text:
using Test, Onda, Dates, MsgPack

@testset "round trip" begin
    mktempdir() do root
        # generate a test dataset
        dataset = Dataset(joinpath(root, "test"); create=true)
        @test dataset isa Dataset
        @test isdir(dataset.path)
        @test isdir(joinpath(dataset.path, "samples"))
        duration_in_seconds = Second(10)
        duration_in_nanoseconds = Nanosecond(duration_in_seconds)
        uuid, recording = create_recording!(dataset)
        Ts = (UInt8, UInt16, UInt32, UInt64, Int8, Int16, Int32, Int64)
        sample_rate = 50.5
        signals = Dict(Symbol(:x, i) => Signal(Symbol.([:a, :b, :c], i), Nanosecond(0),
                                               duration_in_nanoseconds, Symbol(:unit, i),
                                               0.25, i, T, sample_rate, Symbol("lpcm.zst"),
                                               nothing) for (i, T) in enumerate(Ts))
        samples = Dict(k => Samples(v, true, rand(v.sample_type, 3, sample_count(v)))
                       for (k, v) in signals)
        for (name, s) in samples
            @test channel_count(s) == length(s.signal.channel_names)
            @test channel_count(s.signal) == length(s.signal.channel_names)
            @test sample_count(s) == size(s.data, 2)
            @test sizeof_samples(s.signal) == sizeof(s.data)
            @test encode(s) === s
            tmp = similar(s.data)
            encode!(tmp, s)
            @test tmp == s.data
            d = decode(s)
            @test decode(d) === d
            tmp = similar(d.data)
            decode!(tmp, d)
            @test tmp == d.data
            tmp = similar(d.data)
            decode!(tmp, s)
            @test tmp == d.data
            tmp = similar(d.data)
            decode!(tmp, s.signal.sample_resolution_in_unit, s.signal.sample_offset_in_unit,
                    s.data)
            @test tmp == d.data
            @test d.data == (s.data .* s.signal.sample_resolution_in_unit .+
                   s.signal.sample_offset_in_unit)
            if sizeof(s.signal.sample_type) >= 8
                # decoding from 64-bit to floating point is fairly lossy
                tmp = similar(s.data)
                @test isapprox(encode(d, nothing).data, s.data; rtol=10)
                encode!(tmp, d, nothing)
                @test isapprox(tmp, s.data; rtol=10)
                @test isapprox(encode(d, missing).data, s.data; rtol=10)
                encode!(tmp, d, missing)
                @test isapprox(tmp, s.data; rtol=10)
            else
                tmp = similar(s.data)
                encode!(tmp, d, nothing)
                @test tmp == s.data
                @test encode(d, nothing).data == s.data
                encode!(tmp, d, missing)
                @test isapprox(tmp, s.data; rtol=1)
                @test isapprox(encode(d, missing).data, s.data; rtol=1)
            end
            chs = s.signal.channel_names
            i = 93
            @test s[:, i].data == s.data[:, i:i]
            t = Onda.time_from_index(s.signal.sample_rate, i)
            t2 = t + Nanosecond(Second(3))
            j = Onda.index_from_time(s.signal.sample_rate, t2) - 1
            for ch_inds in (:, 1:2, 2:3, 1:3, [3, 1], [2, 3, 1], [1, 2, 3])
                @test s[chs[ch_inds], t].data == s[ch_inds, i].data
                @test s[chs[ch_inds], TimeSpan(t, t2)].data == s.data[ch_inds, i:j]
                @test s[chs[ch_inds], i:j].data == s.data[ch_inds, i:j]
                @test s[ch_inds, t].data == s[ch_inds, i].data
                @test s[ch_inds, TimeSpan(t, t2)].data == s.data[ch_inds, i:j]
                @test s[ch_inds, i:j].data == s.data[ch_inds, i:j]
            end
            @test size(s[:, TimeSpan(0, Second(1))].data, 2) == floor(s.signal.sample_rate)
            for i in 1:length(chs)
                @test channel(s, chs[i]) == i
                @test channel(s, i) == chs[i]
                @test channel(s.signal, chs[i]) == i
                @test channel(s.signal, i) == chs[i]
            end
            @test duration(s) == duration_in_seconds
            @test s[:, TimeSpan(0, duration(s))].data == s.data
            store!(dataset, uuid, name, s)
        end
        save_recordings_file(dataset)

        # read back in the test dataset, add some annotations
        old_dataset = dataset
        dataset = Dataset(joinpath(root, "test"))
        @test length(dataset.recordings) == 1
        uuid, recording = first(dataset.recordings)
        x1 = load(dataset, uuid, :x1)
        @test x1.signal == signals[:x1]
        xs = load(dataset, uuid, (:x3, :x2))
        @test xs[:x3].signal == signals[:x3]
        @test xs[:x2].signal == signals[:x2]
        xs = load(dataset, uuid)
        span = TimeSpan(Second(1), Second(2))
        xs_span = load(dataset, uuid, span)
        for (name, s) in samples
            xi = xs[name]
            @test xi.signal == signals[name]
            @test xi.encoded
            @test xi.data == s.data
            @test xi[:, span].data == xs_span[name].data
        end
        for i in 1:3
            annotate!(recording,
                      Annotation("value_$i", Nanosecond(i),
                                 Nanosecond(i + rand(1:1000000))))
        end
        save_recordings_file(dataset)

        # read back in annotations
        old_uuid = uuid
        old_recording = recording
        old_dataset = dataset
        dataset = Dataset(joinpath(root, "test"))
        uuid, recording = first(dataset.recordings)
        @test old_recording == recording
        delete!(dataset.recordings, uuid)
        uuid, recording = create_recording!(dataset)
        foreach(x -> annotate!(recording, x), old_recording.annotations)
        foreach(x -> store!(dataset, uuid, x, load(old_dataset, old_uuid, x)),
                keys(old_recording.signals))
        merge!(dataset, old_dataset; only_recordings=true)
        @test length(dataset.recordings) == 2
        r1 = dataset.recordings[old_uuid]
        @test r1 == old_recording
        r2 = dataset.recordings[uuid]
        @test r2 == recording
        @test old_uuid != uuid
        @test r1.signals == r2.signals
        @test r1.annotations == r2.annotations

        old_duration = duration(r2)
        new_duration = old_duration + Nanosecond(1)
        r2signals = set_span!(r2, TimeSpan(Nanosecond(0), new_duration))
        @test keys(r2signals) == keys(r2.signals)
        @test all(duration.(values(r2signals)) .== new_duration)
        set_span!(r2, TimeSpan(Nanosecond(0), old_duration))

        r = dataset.recordings[uuid]
        original_signals_length = length(r.signals)
        signal_name, signal = first(r.signals)
        signal_samples = load(dataset, uuid, signal_name)
        signal_samples_path = samples_path(dataset, uuid, signal_name)
        delete!(dataset, uuid, signal_name)
        @test r === dataset.recordings[uuid]
        @test length(r.signals) == (original_signals_length - 1)
        @test !haskey(r.signals, signal_name)
        @test !isfile(signal_samples_path)
        store!(dataset, uuid, signal_name, signal_samples)

        # read back everything, but without assuming an order on the metadata
        dataset = Dataset(joinpath(root, "test"))
        Onda.write_recordings_file(dataset.path,
                                   Onda.Header(dataset.header.onda_format_version, false),
                                   dataset.recordings)
        dataset = Dataset(joinpath(root, "test"))
        @test Dict(old_uuid => old_recording) == dataset.recordings
        delete!(dataset, old_uuid)
        save_recordings_file(dataset)

        # read back the dataset that should now be empty
        dataset = Dataset(joinpath(root, "test"))
        @test isempty(dataset.recordings)
        @test !isdir(joinpath(dataset.path, "samples", string(old_uuid)))

        # make sure samples directory is appropriately created if not present
        no_samples_path = joinpath(root, "no_samples_dir.onda")
        mkdir(no_samples_path)
        cp(joinpath(dataset.path, "recordings.msgpack.zst"),
           joinpath(no_samples_path, "recordings.msgpack.zst"))
        Dataset(no_samples_path; create=false)
        return @test isdir(joinpath(no_samples_path, "samples"))
    end
end

@testset "Error conditions" begin
    mktempdir() do root
        mkdir(joinpath(root, "i_exist.onda"))
        touch(joinpath(root, "i_exist.onda", "memes"))
        @test_throws ArgumentError Dataset(joinpath(root, "i_exist.onda"); create=true)

        dataset = Dataset(joinpath(root, "okay.onda"); create=true)
        uuid, recording = create_recording!(dataset)
        signal = Signal([:a], Nanosecond(0), Nanosecond(Second(10)), :mv, 0.25, 0.0, Int8,
                        100, Symbol("lpcm.zst"), nothing)
        @test_throws DimensionMismatch Samples(signal, true, rand(Int8, 2, 10))
        @test_throws ArgumentError Samples(signal, true, rand(Float32, 1, 10))
        samples = Samples(signal, true, rand(Int8, 1, 10 * 100))
        @test_throws ArgumentError store!(dataset, uuid, Symbol("***HI***"), samples)
        store!(dataset, uuid, :name_okay, samples)
        @test_throws ArgumentError store!(dataset, uuid, :name_okay, samples;
                                          overwrite=false)

        @test_throws ArgumentError Annotation("hi", Nanosecond(20), Nanosecond(4))

        mkdir(joinpath(root, "other.onda"))
        other = Dataset(joinpath(root, "other.onda"); create=true)  # Using existing empty directory
        create_recording!(other, uuid)
        @test_throws ArgumentError create_recording!(other, uuid)
        store!(other, uuid, :cool_stuff, samples)
        @test_throws ErrorException merge!(dataset, other; only_recordings=false)
        return @test_throws ArgumentError merge!(dataset, other; only_recordings=true)
    end
end

@testset "upgrade_onda_format_from_v0_2_to_v0_3!" begin
    mktempdir() do new_path
        old_path = joinpath(@__DIR__, "old_test_v0_2.onda")
        cp(old_path, new_path; force=true)
        dataset = Onda.upgrade_onda_format_from_v0_2_to_v0_3!(new_path,
                                                              (k, v) -> string(k, '.', v))
        @test dataset.path == new_path
        @test dataset.header.onda_format_version == v"0.3.0"
        @test dataset.header.ordered_keys
        old_recordings = MsgPack.unpack(Onda.zstd_decompress(read(joinpath(old_path,
                                                                           "recordings.msgpack.zst"))))[2]
        new_customs = MsgPack.unpack(Onda.zstd_decompress(read(joinpath(new_path,
                                                                        "recordings_custom.msgpack.zst"))))
        @test length(dataset.recordings) == 1
        @test length(new_customs) == 1
        uuid = first(keys(dataset.recordings))
        recording = first(values(dataset.recordings))
        old_recording = first(values(old_recordings))
        @test string(uuid) == first(keys(new_customs)) == first(keys(old_recordings))
        @test first(values(new_customs)) == old_recording["custom"]
        sorted_annotations = sort(collect(recording.annotations); by=first)
        sorted_old_annotations = sort(old_recording["annotations"];
                                      by=(x -> x["start_nanosecond"]))
        @test length(sorted_annotations) == length(sorted_old_annotations)
        for (ann, old_ann) in zip(sorted_annotations, sorted_old_annotations)
            @test ann.value == string(old_ann["key"], '.', old_ann["value"])
            @test ann.start_nanosecond.value == old_ann["start_nanosecond"]
            @test ann.stop_nanosecond.value == old_ann["stop_nanosecond"]
        end
        old_signals = old_recording["signals"]
        @test keys(recording.signals) == Set(Symbol.(keys(old_signals)))
        return for (signal_name, signal) in recording.signals
            old_signal = old_signals[string(signal_name)]
            @test signal.channel_names == Symbol.(old_signal["channel_names"])
            @test signal.start_nanosecond == Nanosecond(0)
            @test signal.stop_nanosecond ==
                  Nanosecond(old_recording["duration_in_nanoseconds"])
            @test signal.sample_unit == Symbol(old_signal["sample_unit"])
            @test signal.sample_resolution_in_unit ==
                  old_signal["sample_resolution_in_unit"]
            @test signal.sample_offset_in_unit == 0.0
            @test signal.sample_type ==
                  Onda.julia_type_from_onda_sample_type(old_signal["sample_type"])
            @test signal.sample_rate == old_signal["sample_rate"]
            @test signal.file_extension == Symbol(old_signal["file_extension"])
            @test signal.file_options == old_signal["file_options"]
        end
    end
end


text:
module ReindexTests

using Test
using FillArrays
using Gridap.Arrays
using Gridap.Arrays: AppliedArray

a = [1,2,3]
b = reindex(a,[3,2,1])
for i=1:3
    b[i]=i
end
@test b == [1,2,3]
@test a == [3,2,1]

a = [[1,2,4,5],[2,4,6,7],[4,3,5,1],[2,3]]
b = [3,1,2]

c = reindex(a,b)

d = a[b]
test_array(c,d)

a = Fill(30.0,10)
b = [3,1,2]

c = reindex(a,b)
@test isa(c,Fill)

d = a[b]
test_array(c,d)

a = Fill(1.0,5,5)
b = [13 23; 15 25]

c = reindex(a,b)
@test isa(c,Fill)

d = a[b]
test_array(c,d)

a = CompressedArray([30,40,10,20,30],[1,2,3,5,3,1,4,2])
b = [3,1,2]

c = reindex(a,b)
@test isa(c,CompressedArray)

d = a[b]
test_array(c,d)

x = [1,2,3,5,3,1,4,2]
a = apply(-,x)
b = [3,1,2]

c = reindex(a,b)
@test isa(c,AppliedArray)

d = [a[bi] for bi in b]
test_array(c,d)

d = a[b]
test_array(c,d)

x = [[1,2,4,5],[2,4,6,7],[4,3,5,1],[2,3]]
a = Table(x)
b = [3,1,2]
c = reindex(a,b)
d = a[b]
test_array(c,d)

end # module

ext: jl


text:
module StippleLatex

import Genie
import Stipple

using Stipple.Reexport

export latex, @latex, @latex!!

const DEFAULT_WRAPPER = Genie.Renderer.Html.template
const COMPONENTS = ["'vue-katex'" => :vueKatex]

#===#

function deps() :: String
  Genie.Router.route("/js/stipplekatex/katex.min.css") do
    Genie.Renderer.WebRenderable(
      read(joinpath(@__DIR__, "..", "files", "katex.min.css"), String),
      :css) |> Genie.Renderer.respond
  end
  Genie.Router.route("/js/stipplekatex/katex.min.js") do
    Genie.Renderer.WebRenderable(
      read(joinpath(@__DIR__, "..", "files", "katex.min.js"), String),
      :javascript) |> Genie.Renderer.respond
  end
  Genie.Router.route("/js/stipplekatex/auto-render.min.js") do
    Genie.Renderer.WebRenderable(
      read(joinpath(@__DIR__, "..", "files", "auto-render.min.js"), String),
      :javascript) |> Genie.Renderer.respond
  end
  Genie.Router.route("/js/stipplekatex/vue-katex.min.js") do
    Genie.Renderer.WebRenderable(
      read(joinpath(@__DIR__, "..", "files", "vue-katex.min.js"), String),
      :javascript) |> Genie.Renderer.respond
  end

  Genie.Router.route("/js/stipplekatex/fonts/KaTeX_Main-Regular.woff2") do
    Genie.Renderer.WebRenderable(
      read(joinpath(@__DIR__, "..", "files", "fonts", "KaTeX_Main-Regular.woff2"), String),
      :fontwoff2) |> Genie.Renderer.respond
  end
  Genie.Router.route("/js/stipplekatex/fonts/KaTeX_Math-Italic.woff2") do
    Genie.Renderer.WebRenderable(
      read(joinpath(@__DIR__, "..", "files", "fonts", "KaTeX_Math-Italic.woff2"), String),
      :fontwoff2) |> Genie.Renderer.respond
  end

  string(
    Genie.Renderer.Html.script(src="$(Genie.config.base_path)js/stipplekatex/katex.min.js"),
    Genie.Renderer.Html.script(src="$(Genie.config.base_path)js/stipplekatex/auto-render.min.js", onload="renderMathInElement(document.body);"),
    Genie.Renderer.Html.script(src="$(Genie.config.base_path)js/stipplekatex/vue-katex.min.js", defer=true),
    Genie.Renderer.Html.link(href="$(Genie.config.base_path)js/stipplekatex/katex.min.css", rel="stylesheet"),
  )
end

#===#

function __init__()
  Stipple.DEPS[@__MODULE__] = deps
end

#===#

function latex(content::String = "",
                args...;
                wrap::Function = DEFAULT_WRAPPER,
                expression::String = "",
                auto::Bool = true,
                display::Bool = false,
                throw_on_error::Bool = false,
                error_color::String = "#CC0000",
                max_size::String="Infinity",
                max_expand::Union{String,Number} = 1000,
                allowed_protocols::Vector{String} = String[],
                strict::Union{Bool,String} = "warn",
                kwargs...)
  isempty(expression) || (content = expression)
  options = Dict(:displayMode => display,
                  :throwOnError => throw_on_error,
                  :errorColor => error_color,
                  :maxSize => max_size,
                  :maxExpand => max_expand,
                  :allowedProtocols => allowed_protocols,
                  :strict => strict)
  arguments = "{ expression: '$(content)', options: $(Stipple.JSONParser.json(options)) }"
  arguments = replace(arguments, '"'=>"'")

  wrap() do
    auto ?
      Genie.Renderer.Html.div(v__katex!!auto = arguments, args...; kwargs...) :
      Genie.Renderer.Html.div(v__katex = arguments, args...; kwargs...)
  end
end

macro latex(expr)
  quote
    "v-katex:auto='$($(esc(expr)))'"
  end
end

macro latex!!(expr)
  quote
    "v-katex:display='$($(esc(expr)))'"
  end
end

end


text:
#
println("Hb) Tests of the StrongField module to calculate energy and momentum distributions of photoelectron in SFA.")
#

dataName        = "Neon_StrongFieldIonization" #If files are exported below (plot, data), all filenames will include dataName

#########################################################
#-------------------PHYSICAL SETTINGS--------------------
#########################################################

#_____________1. Atomic target______________________________________________

#An initial and final level for the strong-field ionization process will be computed
#In the strong-field computation, the highest lying one-electron  orbital will be taken as "active electron"

#Specify the initial atomic and final ionic states
initialName     = "Li ground-state"
nuclearModel    = Nuclear.Model(3.39336) #Adjust the nuclear charge so that the ionization potential of the selected active electron fits (Li: 3.39336, Ne: 10.463599, Ar: 18.50435055, Kr: 36.529677, Xe: 54.5504443)
refConfigInitial= [Configuration("[He] 2s")]

finalName       = "He ground-state"
refConfigFinal  = [Configuration("[He]")]

hydrogenic      = false             #True: Use hydrogenic wave functions for the initial state (quantum numbers n,l,m and the ionization potential are taken from initialLevel)
hydrogenic1s    = false             #True (AND hydrogenic = true): Same as hydrogenic, but quantum numbers n,l,m of initial level are ignored and a 1s state is used instead (ionization potential is taken from initialLevel)

#___________________________________________________________________________

#_____________2. Laser field and electron continuum_________________________

#---Define laser beam parameters---------------
wavelength      = 800.                                                          #nm
intensity       = 5e14                                                          #W/cm^2
CEP             = 0.                                                            #carrier-envelope phase 0 <= CEP < 2*pi
envelope        = Pulse.SinSquaredEnvelope(8)                                   #InfiniteEnvelope(), SinSquaredEnvelope(np), GaussianEnvelope(np) with the number of optical cycles np
polarization    = Basics.RightElliptical(0.36)                                  #RightCircular(), LeftCircular(),
                                                                                #RightElliptical(0 <= ellipticity <= 1), LeftElliptical(0 <= ellipticity <= 1)
                                                                                #NOTE: Elliptical polarization is only implemented for SinSquaredEnvelope so far

#---Choose electron continuum in the laser field-------------------
volkov          = StrongField.FreeVolkov()                                #FreeVolkov(), CoulombVolkov(Z) with the charge of the ion Z, or DistortedVolkov()

#-----AUTOMATIC: derived quantities and settings-------------------
omega           = convertUnits("energy: from wavelength [nm] to atomic", wavelength)
intensity       = convertUnits("intensity: from W/cm^2 to atomic", intensity)
A0              = Pulse.computeFieldAmplitude(intensity, omega)
beam            = Pulse.PlaneWaveBeam(A0, omega, CEP)                            #Arguments: A0, omega, CEP #Only PlaneWaveBeam possible

#___________________________________________________________________________

#_____________3. Observable that will be computed__________________________

observable      = StrongField.SfaEnergyDistribution(pi/2, pi/2, 10, 10*omega)          #Arguments: theta, phi, number of energies, maximum energy
#observable      = StrongField.SfaMomentumDistribution(pi/2, 30, 100, 10*omega)         #Arguments: theta, number of azimuthal angles, number of energies, maximum energy
#observable      = StrongField.SfaAzimuthalAngularDistribution(pi/2, 100, 2.21*omega)   #Arguments: theta, number of azimuthal angles, energy
#observable      = StrongField.SfaPolarAngularDistribution(pi/2, 100, 2.5*omega)        #Arguments: phi, number of azimuthal angles, energy

#___________________________________________________________________________

#########################################################
#-----------------COMPUTATIONAL SETTINGS-----------------
#########################################################

asfSettings     = AsfSettings(AsfSettings(), generateScf=true)
rGrid           = Radial.Grid(Radial.Grid(false), rnt = 4.0e-6, h = 5.0e-2, hp = 2.0e-2, rbox = 20.0) #Radial grid for atomic computation and matrix elements

multipoles      = [E1]              #Only "dipole" ( [E1] ) SFA implemented so far
gauges          = [UseCoulomb]      #Only Coulomb (velocity) gauge implemented so far
mAverage        = true              #True: Average over initialState projections mj or ml and sum over final state spin projection msp
coupledBasis    = true              #True: compute in |j l mj> angular momentum basis; False: compute in |l ml> angular momentum basis

printAmplitudes = true              #True: Print all requested amplitudes during computation

exportData      = true              #True: Probabilities will be saved to dataName-observable-j.csv as (E_p,P) or (p_x,p_y,P)
exportPlot      = true              #True: Energy distributions or Angular distributions will be plotted and saved to dataName-observable.pdf
plotEnergyScale = "omega"           #Energy scaling used in plot ("atomic" or "omega")
plotProbabilityScale = "linear"     #Energy scaling used in plot ("linear" or "log" - "log" here means log10)
exportRadialWF  = true              #True: Radial wave functions of initial states used in StrongField computations will be saved to dataName-initial_radial_wave_function-j.csv
plotRadialWF    = true              #True: Radial wave functions of initial states used in StrongField computations will be plotted and saved to dataName-initial_radial_wave_function.pdf

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

##############################################################
#-----------------COMPUTATION AND DATA EXPORT-----------------
##############################################################

#_____________1. Atomic computations: Compute the initial and final levels for the SFA computation_____
AtomicComp              = Atomic.Computation(Atomic.Computation(), name=initialName, grid=rGrid, nuclearModel=nuclearModel, configs=refConfigInitial,  asfSettings=asfSettings )
AtomicData              = perform(AtomicComp, output=true)
initialLevel            = AtomicData["multiplet:"].levels[1]
#
AtomicComp              = Atomic.Computation(Atomic.Computation(), name=finalName, grid=rGrid, nuclearModel=nuclearModel, configs=refConfigFinal,  asfSettings=asfSettings )
AtomicData              = perform(AtomicComp, output=true)
finalLevel              = AtomicData["multiplet:"].levels[1]

#_____________2. StrongField computation(s)___________
sfaSettings         = StrongField.Settings(multipoles, gauges, printAmplitudes, coupledBasis, hydrogenic, hydrogenic1s, mAverage)
SFIComp1            = StrongField.Computation(observable, nuclearModel, rGrid, initialLevel, finalLevel, beam, envelope, polarization, volkov, sfaSettings)
SFIData1            = StrongField.perform(SFIComp1, output=true)

compArray           = [SFIComp1] #Multiple computations can be combined here into an array and will be plotted and exported below
dataArray           = [SFIData1]

#_____________3. Export data and generate plots___________

if exportData
    StrongField.exportData(compArray, dataArray, dataName)
end

if exportPlot
    StrongField.plot(compArray, dataArray, plotEnergyScale, plotProbabilityScale, dataName)
end

if exportRadialWF
    StrongField.exportRadialWavefunctions( compArray, dataName, plotRadialWF )
end


text:
#=
The functions here are intended to generate AdmittanceModels.Circuit objects
from ANSYS Q3D plain text output files containing RLGC parameters.

Here's a small example of the expected Q3D file format, valid for ANSYS Electronics Desktop 2016:

################################################################################################################
DesignVariation:\$DummyParam1='1' \$DummyParam2='1mm' \$dummy_param3='100e19pF' \$dummy_param4='1.2eUnits'
Setup1:LastAdaptive
Problem Type:C
C Units:farad, G Units:sie
Reduce Matrix:Original
Frequency: 4E+09 Hz

Capacitance Matrix
	net_name_1	net_name_2
net_name_1	1.1924E-13	-1.3563E-16
net_name_2	-1.3563E-16	3.7607E-13

Conductance Matrix
	net_name_1	net_name_2
net_name_1	2.4489E-09	-2.3722E-12
net_name_2	-2.3722E-12	7.7123E-09

Capacitance Matrix Coupling Coefficient
	net_name_1	net_name_2
net_name_1	1	0.00064047	0.00018631
net_name_2	0.00064047	1	0.00010129

Conductance Matrix Coupling Coefficient
	net_name_1	net_name_2
net_name_1	1	0.00054586	0.00017696
net_name_2	0.00054586	1	9.2987E-05
################################################################################################################

Each data block like \"Capacitance Matrix\" must start with a line enumerating
the N net names. These are the column labels for the given matrix. The matrix
rows are specified directly below this line. Each of the N row lines start with
a net name to label the row, followed by N space-separated float values for the
matrix elements. We expect these to be strings `parse(Float64, _)` can handle.
=#

using LinearAlgebra: diagind

"""
    parse_value_and_unit(s::AbstractString)

Parse strings like "9.9e-9mm" => (9.9e-9, "mm").
"""
function parse_value_and_unit(s::AbstractString)
    # Three capture groups in the regex intended to match
    # (±X.X)(e-X)(unit)
    m = match(r"(\-?[\d|\.]+)(e-?[\d]+)?(\D*)", s)
    exp_str = m[2] == nothing ? "" : m[2]
    num_str = m[1] * exp_str
    num = parse(Float64, num_str)
    unit_str = String(m[3])
    (num, unit_str)
end

const matrix_types = Dict(
    :capacitance => "Capacitance Matrix",
    :conductance => "Conductance Matrix",
    # :dc_inductance => "DC Inductance Matrix",
    # :dc_resistance => "DC Resistance Matrix",
    # :ac_inductance => "AC Inductance Matrix",
    # :ac_resistance => "AC Resistance Matrix",
    # :dc_plus_ac_resistance => "DCPlusAC Resistance Matrix",
)

const unit_names = Dict(
    :capacitance => "C Units",
    :conductance => "G Units"
)

const units_to_multipliers = Dict(
    "fF" => 1e-15,
    "pF" => 1e-12,
    "nF" => 1e-9,
    "uF" => 1e-6,
    "mF" => 1e-3,
    "farad" => 1.0,
    "mSie" => 1e-3,
    "sie" => 1.0
)

"""
    parse_q3d_txt(filepath, matrix_type::Symbol)

Parses a plain text file generated by ANSYS Q3D containing RLGC simulation output.

# Args
- `q3d_file`: the path to the text file.
- `matrix_type`: A `Symbol` indicating the type of data to read. Currently,
  symbols in $(collect(keys(matrix_types))) are the only available choices,
  and only the units in $(collect(keys(units_to_multipliers))) are handled.

# Returns
- A tuple `(design_variation, net_names, matrix)`, where:
    - `design_variation` is a `Dict` of design variables
    - `net_names` is a `Vector` of net names (strings)
    - `matrix` is the requested matrix
"""
function parse_q3d_txt(q3d_file, matrix_type::Symbol)
    !haskey(matrix_types, matrix_type) && error("matrix type not implemented.")

    matrix_name = matrix_types[matrix_type]
    unit_name = unit_names[matrix_type]

    # Windows uses \r\n for newlines, where Linux just uses \n.
    # We work around that difference here by deleting all carriage return (\r)
    # characters from the file.
    linesep = "\n"
    file_text = replace(read(q3d_file, String), "\r" => "")
    file_chunks = split(file_text, linesep * linesep)
    @assert length(file_chunks) > 1 "expected more than one block in Q3D plain text file."

    is_chunk(header) = chunk -> startswith(chunk, header)
    function get_chunk(header)
        # the chunk may be missing, guard against errors
        idx = findfirst(is_chunk(header), file_chunks)
        !isnothing(idx) && return file_chunks[idx]
        return nothing
    end
    get_chunk_line(chunk, line_num) = split(chunk, linesep)[line_num]

    # The following parses lines like:
    # DesignVariation:var1='12' var2='1e-05' var3='123um'
    # into Dict("var1" => (12, ""),
    #           "var2" => (1e-5, ""),
    #           "var3" => (123, "um"))
    design_variation_chunk = get_chunk("DesignVariation")
    design_variation = if !isnothing(design_variation_chunk)
        _, design_variation_data_str = split(get_chunk_line(design_variation_chunk, 1), ":")
        design_strings = split(design_variation_data_str)
        parse_design_kv((k, v)) = String(k) => parse_value_and_unit(strip(v, ['\'']))
        Dict(parse_design_kv.(split.(design_strings, ['='])))
    else
        Dict{String, Tuple{Float64, String}}()
    end

    local unit_multiplier
    for l in split(file_chunks[1], linesep)
        if occursin("Units", l)
            unit = Dict(split.(split(l, ", "), ':'))[unit_name]
            !haskey(units_to_multipliers, unit) && error("unit $unit not implemented.")
            unit_multiplier = units_to_multipliers[unit]
            break
        end
    end
    (@isdefined unit_multiplier) || error("units not given in Q3D plain text file.")

    # Don't forget to include the line separator here.
    # We typically have matrix_name = "Capacitance Matrix"
    # But, these files can have a heading like "Capacitance Matrix Coupling Coefficient"
    matrix_chunk = get_chunk(matrix_name * linesep)
    net_names = String.(split(get_chunk_line(matrix_chunk, 2)))

    get_matrix_row(row_idx) = parse.(Float64,
        split(get_chunk_line(matrix_chunk, 2 + row_idx))[2:end]) .* unit_multiplier
    matrix = reduce(hcat, get_matrix_row.(1:length(net_names)))

    (design_variation, net_names, matrix)
end

"""
    Circuit(q3d_file; matrix_types = [:capacitance])

Return an `AdmittanceModels.Circuit` from a plain text file generated by
ANSYS Q3D containing RLGC simulation output.

# Args
- `q3d_file`: the path to the text file.
- `matrix_type`: A list of symbols indicating which type of matrix data to read.
  Currently, symbols in $(collect(keys(matrix_types))) are the only available
  choices, and only the units in $(collect(keys(units_to_multipliers))) are handled.
"""
function Circuit(q3d_file; matrix_types = [:capacitance])
    parse_dict = Dict(t => parse_q3d_txt(q3d_file, t) for t ∈ matrix_types)

    # Check that all matrix types are defined over the same set of nets
    # `take_only` asserts the iterator it's passed has only one element
    nets = take_only(collect(Set(map(t -> t[2], values(parse_dict)))))
    matrix_dict = Dict(t => parse_dict[t][3] for t ∈ matrix_types)

    zero_mat() = zeros(length(nets), length(nets))

    # The C and G matrices expected in AdmittanceModels.Circuit are the same
    # as those parsed from Q3D, except with all diagonal entries set to 0 and
    # all values made positive.
    function prep_matrix(matrix_type)
        x = get(matrix_dict, matrix_type, zero_mat())
        x[diagind(x)] .= 0.
        abs.(x)
    end

    k = zero_mat()
    g = prep_matrix(:conductance)
    c = prep_matrix(:capacitance)
    return Circuit(k, g, c, nets)
end

function take_only(xs)
    @assert length(xs) == 1 "Expected $xs to have one element."
    return xs[1]
end


text:
module NB5I

using NBodyIPs.FastPolys

using StaticArrays

# Primary invariants for NBody=5and deg=7
 # : definitions at the beginning of the file
const P1_1 = (1,2,3,4,5,6,7,8,9,10,)

const P2_1 = (1,1,1,2,2,3,1,1,1,2,3,4,2,2,3,4,3,4,5,5,6,5,5,6,7,6,7,8,8,9,)
const P2_2 = (2,3,4,3,4,4,5,6,7,5,6,7,8,9,8,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,)

const P3_1 = (1,2,3,4,5,6,7,8,9,10,)

const P4_1 = (1,1,1,2,1,1,1,2,2,3,4,3,4,2,3,4,5,5,6,7,)
const P4_2 = (2,2,3,3,5,5,6,5,5,6,7,6,7,8,8,9,6,8,8,9,)
const P4_3 = (3,4,4,4,6,7,7,8,9,8,9,10,10,9,10,10,7,9,10,10,)

const P5_1 = (1,2,3,4,5,6,7,8,9,10,)

const P6_1 = (1,1,2,3,4,)
const P6_2 = (2,5,5,6,7,)
const P6_3 = (3,6,8,8,9,)
const P6_4 = (4,7,9,10,10,)

const P7_1 = (1,2,3,4,5,6,7,8,9,10,)

const P8_1 = (1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,3,3,1,1,1,2,2,3,4,3,4,2,3,4,)
const P8_2 = (2,2,2,2,2,2,2,3,4,2,3,4,3,4,3,4,4,4,5,5,5,5,5,5,5,6,6,5,6,7,)
const P8_3 = (3,3,3,3,3,3,5,5,5,5,6,7,5,5,6,7,6,7,6,6,6,6,7,6,7,7,7,8,8,8,)
const P8_4 = (4,4,4,4,4,4,6,6,6,8,8,9,8,8,8,9,8,9,7,7,7,8,8,8,9,8,9,9,9,9,)
const P8_5 = (5,6,7,8,9,10,7,7,7,9,10,10,9,9,10,10,10,10,8,9,10,9,9,10,10,10,10,10,10,10,)

const P9_1 = (1,2,3,4,5,6,7,8,9,10,)

const P10_1 = (1,2,3,4,5,6,7,8,9,10,)

# Irreducible secondaries for NBody=5and deg=7
 # : definitions at the beginning of the file
const IS1_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS1_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS2_1 = (1,1,1,2,2,3,5,5,6,8,)
const IS2_2 = (2,3,4,3,4,4,6,7,7,9,)
const IS2_3 = (5,6,7,8,9,10,8,9,10,10,)

const IS3_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS3_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS4_1 = (1,1,1,2,2,3,1,1,1,2,3,4,2,2,3,4,3,4,5,5,6,5,5,6,7,6,7,8,8,9,)
const IS4_2 = (2,3,4,3,4,4,5,6,7,5,6,7,8,9,8,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,)

const IS5_1 = (1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,8,9,10,5,5,6,7,6,7,8,9,10,8,9,10,)
const IS5_2 = (2,3,4,1,1,1,3,4,2,2,4,3,1,1,1,2,2,3,6,7,5,5,7,6,5,5,6,9,8,8,)
const IS5_3 = (5,6,7,5,6,7,8,9,8,9,10,10,2,3,4,3,4,4,8,9,8,9,10,10,6,7,7,10,10,9,)

const IS6_1 = (1,1,1,2,2,3,4,3,4,2,3,4,1,1,1,2,2,3,4,3,4,2,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,6,7,5,6,7,8,9,8,9,10,10,)
const IS6_2 = (2,2,3,1,1,1,1,1,1,3,2,2,5,5,6,5,5,6,7,6,7,8,8,9,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,5,5,8,8,9,5,5,6,7,6,7,)
const IS6_3 = (3,4,4,3,4,2,2,4,3,4,4,3,6,7,7,8,9,8,9,10,10,9,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,7,7,6,9,10,10,9,8,10,10,8,9,)

const IS7_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,8,9,8,9,10,10,)
const IS7_2 = (3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,4,3,4,3,2,2,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,7,6,7,6,5,5,6,7,5,5,7,6,)
const IS7_3 = (5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,8,9,8,9,10,10,9,8,10,10,8,9,)

const IS8_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS8_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS9_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS9_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS10_1 = (1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,8,9,10,5,5,6,7,6,7,8,9,10,8,9,10,)
const IS10_2 = (2,3,4,1,1,1,3,4,2,2,4,3,1,1,1,2,2,3,6,7,5,5,7,6,5,5,6,9,8,8,)
const IS10_3 = (5,6,7,5,6,7,8,9,8,9,10,10,2,3,4,3,4,4,8,9,8,9,10,10,6,7,7,10,10,9,)

const IS11_1 = (1,1,1,2,2,3,1,1,1,2,3,4,2,2,3,4,3,4,5,5,6,5,5,6,7,6,7,8,8,9,)
const IS11_2 = (2,3,4,3,4,4,5,6,7,5,6,7,8,9,8,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,)
const IS11_3 = (5,6,7,8,9,10,2,3,4,1,1,1,3,4,2,2,4,3,8,9,10,6,7,5,5,7,6,10,9,8,)

const IS12_1 = (1,1,1,2,2,3,4,3,4,2,3,4,1,1,1,2,2,3,4,3,4,2,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,6,7,5,6,7,8,9,8,9,10,10,)
const IS12_2 = (2,2,3,1,1,1,1,1,1,3,2,2,5,5,6,5,5,6,7,6,7,8,8,9,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,5,5,8,8,9,5,5,6,7,6,7,)
const IS12_3 = (3,4,4,3,4,2,2,4,3,4,4,3,6,7,7,8,9,8,9,10,10,9,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,7,7,6,9,10,10,9,8,10,10,8,9,)

const IS13_1 = (1,1,1,1,1,1,2,2,3,2,2,3,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,5,5,6,7,6,7,8,8,9,5,5,6,5,5,6,7,6,7,8,8,9,)
const IS13_2 = (2,2,3,4,3,4,3,4,4,3,4,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,)
const IS13_3 = (3,4,2,2,4,3,1,1,1,4,3,2,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,1,1,1,2,2,3,4,3,4,2,3,4,7,6,5,9,8,10,10,8,9,5,6,7,)

const IS14_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,8,9,8,9,10,10,)
const IS14_2 = (3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,4,3,4,3,2,2,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,7,6,7,6,5,5,6,7,5,5,7,6,)
const IS14_3 = (5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,8,9,8,9,10,10,9,8,10,10,8,9,)

const IS15_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS15_2 = (2,2,2,2,3,3,1,1,1,1,1,1,1,1,1,1,1,1,3,3,2,2,2,2,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,3,3,2,2,2,2,3,3,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,6,5,5,5,5,6,7,5,5,7,6,5,5,5,5,6,6,5,5,6,7,6,7,)
const IS15_3 = (3,4,3,4,4,4,3,4,2,2,4,3,3,4,2,2,4,3,4,4,4,3,4,3,5,5,5,5,6,6,5,5,6,7,6,7,5,5,6,7,6,7,8,8,8,9,8,9,2,2,3,4,3,4,2,2,3,4,3,4,3,4,3,4,4,4,3,4,3,4,4,4,6,7,5,5,7,6,6,7,5,5,7,6,5,5,5,5,6,6,9,8,9,8,8,8,7,7,7,6,7,6,8,8,8,9,8,9,6,7,6,7,7,7,9,8,9,8,8,8,)
const IS15_4 = (5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,6,7,7,7,8,9,8,9,10,10,8,9,8,9,10,10,9,9,10,10,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,8,9,8,9,10,10,8,9,8,9,10,10,6,7,6,7,7,7,10,10,10,10,9,9,8,9,8,9,10,10,9,9,10,10,10,10,9,8,10,10,8,9,10,10,10,10,9,9,)

const IS16_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS16_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS17_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS17_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS18_1 = (1,1,1,2,2,3,1,1,1,2,3,4,2,2,3,4,3,4,5,5,6,5,5,6,7,6,7,8,8,9,)
const IS18_2 = (2,3,4,3,4,4,5,6,7,5,6,7,8,9,8,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,)

const IS19_1 = (1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,8,9,10,5,5,6,7,6,7,8,9,10,8,9,10,)
const IS19_2 = (2,3,4,1,1,1,3,4,2,2,4,3,1,1,1,2,2,3,6,7,5,5,7,6,5,5,6,9,8,8,)
const IS19_3 = (5,6,7,5,6,7,8,9,8,9,10,10,2,3,4,3,4,4,8,9,8,9,10,10,6,7,7,10,10,9,)

const IS20_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS20_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)
const IS20_3 = (5,6,7,5,6,7,8,9,8,9,10,10,2,3,4,1,1,1,3,4,2,2,4,3,2,3,4,1,1,1,3,4,2,2,4,3,8,9,8,9,10,10,6,7,5,5,7,6,6,7,5,5,7,6,10,10,9,8,9,8,)

const IS21_1 = (1,1,1,2,2,3,5,5,6,8,)
const IS21_2 = (2,3,4,3,4,4,6,7,7,9,)
const IS21_3 = (5,6,7,8,9,10,8,9,10,10,)

const IS22_1 = (1,1,1,2,2,3,4,3,4,2,3,4,1,1,1,2,2,3,4,3,4,2,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,6,7,5,6,7,8,9,8,9,10,10,)
const IS22_2 = (2,2,3,1,1,1,1,1,1,3,2,2,5,5,6,5,5,6,7,6,7,8,8,9,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,5,5,8,8,9,5,5,6,7,6,7,)
const IS22_3 = (3,4,4,3,4,2,2,4,3,4,4,3,6,7,7,8,9,8,9,10,10,9,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,7,7,6,9,10,10,9,8,10,10,8,9,)

const IS23_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS23_2 = (2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)
const IS23_3 = (3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,4,3,4,3,2,2,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,7,6,7,6,5,5,9,8,10,10,8,9,9,8,10,10,8,9,5,5,6,7,6,7,)

const IS24_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,8,9,8,9,10,10,)
const IS24_2 = (3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,4,3,4,3,2,2,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,7,6,7,6,5,5,6,7,5,5,7,6,)
const IS24_3 = (5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,8,9,8,9,10,10,9,8,10,10,8,9,)

const IS25_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS25_2 = (2,2,2,2,3,3,1,1,1,1,1,1,1,1,1,1,1,1,3,3,2,2,2,2,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,3,3,2,2,2,2,3,3,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,6,5,5,5,5,6,7,5,5,7,6,5,5,5,5,6,6,5,5,6,7,6,7,)
const IS25_3 = (3,4,3,4,4,4,3,4,2,2,4,3,3,4,2,2,4,3,4,4,4,3,4,3,5,5,5,5,6,6,5,5,6,7,6,7,5,5,6,7,6,7,8,8,8,9,8,9,2,2,3,4,3,4,2,2,3,4,3,4,3,4,3,4,4,4,3,4,3,4,4,4,6,7,5,5,7,6,6,7,5,5,7,6,5,5,5,5,6,6,9,8,9,8,8,8,7,7,7,6,7,6,8,8,8,9,8,9,6,7,6,7,7,7,9,8,9,8,8,8,)
const IS25_4 = (5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,6,7,7,7,8,9,8,9,10,10,8,9,8,9,10,10,9,9,10,10,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,8,9,8,9,10,10,8,9,8,9,10,10,6,7,6,7,7,7,10,10,10,10,9,9,8,9,8,9,10,10,9,9,10,10,10,10,9,8,10,10,8,9,10,10,10,10,9,9,)

const IS26_1 = (1,1,1,1,1,1,2,2,3,2,2,3,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,5,5,6,7,6,7,8,8,9,5,5,6,5,5,6,7,6,7,8,8,9,)
const IS26_2 = (2,2,3,4,3,4,3,4,4,3,4,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,)
const IS26_3 = (3,4,2,2,4,3,1,1,1,4,3,2,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,2,2,3,4,3,4,2,3,4,7,6,5,6,7,5,5,7,6,5,6,7,)
const IS26_4 = (5,5,6,7,6,7,8,9,10,8,9,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,8,9,10,6,7,5,5,7,6,10,9,8,8,9,10,9,8,10,10,8,9,10,9,8,)

const IS27_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS27_2 = (5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,8,9,8,9,10,10,6,7,5,5,7,6,6,7,5,5,7,6,10,10,9,8,9,8,8,9,8,9,10,10,6,7,5,5,7,6,6,7,5,5,7,6,10,10,9,8,9,8,)
const IS27_3 = (3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,4,3,4,3,2,2,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,7,6,7,6,5,5,9,8,10,10,8,9,9,8,10,10,8,9,5,5,6,7,6,7,)

const IS28_1 = (1,2,3,4,5,6,7,8,9,10,)

const IS29_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS29_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS30_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS30_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS31_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS31_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS32_1 = (1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,8,9,10,5,5,6,7,6,7,8,9,10,8,9,10,)
const IS32_2 = (2,3,4,1,1,1,3,4,2,2,4,3,1,1,1,2,2,3,6,7,5,5,7,6,5,5,6,9,8,8,)
const IS32_3 = (5,6,7,5,6,7,8,9,8,9,10,10,2,3,4,3,4,4,8,9,8,9,10,10,6,7,7,10,10,9,)

const IS33_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS33_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)
const IS33_3 = (5,6,7,5,6,7,8,9,8,9,10,10,2,3,4,1,1,1,3,4,2,2,4,3,2,3,4,1,1,1,3,4,2,2,4,3,8,9,8,9,10,10,6,7,5,5,7,6,6,7,5,5,7,6,10,10,9,8,9,8,)

const IS34_1 = (1,1,1,2,2,3,1,1,1,2,3,4,2,2,3,4,3,4,5,5,6,5,5,6,7,6,7,8,8,9,)
const IS34_2 = (2,3,4,3,4,4,5,6,7,5,6,7,8,9,8,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,)
const IS34_3 = (5,6,7,8,9,10,2,3,4,1,1,1,3,4,2,2,4,3,8,9,10,6,7,5,5,7,6,10,9,8,)

const IS35_1 = (1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,8,9,10,5,5,6,7,6,7,8,9,10,8,9,10,)
const IS35_2 = (2,3,4,1,1,1,3,4,2,2,4,3,1,1,1,2,2,3,6,7,5,5,7,6,5,5,6,9,8,8,)
const IS35_3 = (5,6,7,5,6,7,8,9,8,9,10,10,2,3,4,3,4,4,8,9,8,9,10,10,6,7,7,10,10,9,)

const IS36_1 = (1,1,1,2,2,3,4,3,4,2,3,4,1,1,1,2,2,3,4,3,4,2,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,6,7,5,6,7,8,9,8,9,10,10,)
const IS36_2 = (2,2,3,1,1,1,1,1,1,3,2,2,5,5,6,5,5,6,7,6,7,8,8,9,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,5,5,8,8,9,5,5,6,7,6,7,)
const IS36_3 = (3,4,4,3,4,2,2,4,3,4,4,3,6,7,7,8,9,8,9,10,10,9,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,7,7,6,9,10,10,9,8,10,10,8,9,)

const IS37_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS37_2 = (2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)
const IS37_3 = (3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,4,3,4,3,2,2,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,7,6,7,6,5,5,9,8,10,10,8,9,9,8,10,10,8,9,5,5,6,7,6,7,)

const IS38_1 = (1,1,1,1,1,1,2,2,3,2,2,3,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,5,5,6,7,6,7,8,8,9,5,5,6,5,5,6,7,6,7,8,8,9,)
const IS38_2 = (2,2,3,4,3,4,3,4,4,3,4,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,)
const IS38_3 = (3,4,2,2,4,3,1,1,1,4,3,2,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,1,1,1,2,2,3,4,3,4,2,3,4,7,6,5,9,8,10,10,8,9,5,6,7,)

const IS39_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,8,9,8,9,10,10,)
const IS39_2 = (3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,4,3,4,3,2,2,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,7,6,7,6,5,5,6,7,5,5,7,6,)
const IS39_3 = (5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,8,9,8,9,10,10,9,8,10,10,8,9,)

const IS40_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS40_2 = (2,2,2,2,3,3,1,1,1,1,1,1,1,1,1,1,1,1,3,3,2,2,2,2,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,3,3,2,2,2,2,3,3,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,6,5,5,5,5,6,7,5,5,7,6,5,5,5,5,6,6,5,5,6,7,6,7,)
const IS40_3 = (3,4,3,4,4,4,3,4,2,2,4,3,3,4,2,2,4,3,4,4,4,3,4,3,5,5,5,5,6,6,5,5,6,7,6,7,5,5,6,7,6,7,8,8,8,9,8,9,2,2,3,4,3,4,2,2,3,4,3,4,3,4,3,4,4,4,3,4,3,4,4,4,6,7,5,5,7,6,6,7,5,5,7,6,5,5,5,5,6,6,9,8,9,8,8,8,7,7,7,6,7,6,8,8,8,9,8,9,6,7,6,7,7,7,9,8,9,8,8,8,)
const IS40_4 = (5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,6,7,7,7,8,9,8,9,10,10,8,9,8,9,10,10,9,9,10,10,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,8,9,8,9,10,10,8,9,8,9,10,10,6,7,6,7,7,7,10,10,10,10,9,9,8,9,8,9,10,10,9,9,10,10,10,10,9,8,10,10,8,9,10,10,10,10,9,9,)

# Primary invariants for NBody=5and deg=7
 # : definitions of the types at the beginning of the file
const P1 = Val((P1_1,))
const P2 = Val((P2_1,P2_2,))
const P3 = Val((P3_1,))
const P4 = Val((P4_1,P4_2,P4_3,))
const P5 = Val((P5_1,))
const P6 = Val((P6_1,P6_2,P6_3,P6_4,))
const P7 = Val((P7_1,))
const P8 = Val((P8_1,P8_2,P8_3,P8_4,P8_5,))
const P9 = Val((P9_1,))
const P10 = Val((P10_1,))
# Irreducible secondaries for NBody=5and deg=7
 # : definitions of the types at the beginning of the file
const IS1 = Val((IS1_1,IS1_2,))
const IS2 = Val((IS2_1,IS2_2,IS2_3,))
const IS3 = Val((IS3_1,IS3_2,))
const IS4 = Val((IS4_1,IS4_2,))
const IS5 = Val((IS5_1,IS5_2,IS5_3,))
const IS6 = Val((IS6_1,IS6_2,IS6_3,))
const IS7 = Val((IS7_1,IS7_2,IS7_3,))
const IS8 = Val((IS8_1,IS8_2,))
const IS9 = Val((IS9_1,IS9_2,))
const IS10 = Val((IS10_1,IS10_2,IS10_3,))
const IS11 = Val((IS11_1,IS11_2,IS11_3,))
const IS12 = Val((IS12_1,IS12_2,IS12_3,))
const IS13 = Val((IS13_1,IS13_2,IS13_3,))
const IS14 = Val((IS14_1,IS14_2,IS14_3,))
const IS15 = Val((IS15_1,IS15_2,IS15_3,IS15_4,))
const IS16 = Val((IS16_1,IS16_2,))
const IS17 = Val((IS17_1,IS17_2,))
const IS18 = Val((IS18_1,IS18_2,))
const IS19 = Val((IS19_1,IS19_2,IS19_3,))
const IS20 = Val((IS20_1,IS20_2,IS20_3,))
const IS21 = Val((IS21_1,IS21_2,IS21_3,))
const IS22 = Val((IS22_1,IS22_2,IS22_3,))
const IS23 = Val((IS23_1,IS23_2,IS23_3,))
const IS24 = Val((IS24_1,IS24_2,IS24_3,))
const IS25 = Val((IS25_1,IS25_2,IS25_3,IS25_4,))
const IS26 = Val((IS26_1,IS26_2,IS26_3,IS26_4,))
const IS27 = Val((IS27_1,IS27_2,IS27_3,))
const IS28 = Val((IS28_1,))
const IS29 = Val((IS29_1,IS29_2,))
const IS30 = Val((IS30_1,IS30_2,))
const IS31 = Val((IS31_1,IS31_2,))
const IS32 = Val((IS32_1,IS32_2,IS32_3,))
const IS33 = Val((IS33_1,IS33_2,IS33_3,))
const IS34 = Val((IS34_1,IS34_2,IS34_3,))
const IS35 = Val((IS35_1,IS35_2,IS35_3,))
const IS36 = Val((IS36_1,IS36_2,IS36_3,))
const IS37 = Val((IS37_1,IS37_2,IS37_3,))
const IS38 = Val((IS38_1,IS38_2,IS38_3,))
const IS39 = Val((IS39_1,IS39_2,IS39_3,))
const IS40 = Val((IS40_1,IS40_2,IS40_3,IS40_4,))

function invariants_gen(x1::SVector{10, T}) where {T}
   x2 = x1.*x1
   x3 = x2.*x1
   x4 = x3.*x1
   x5 = x4.*x1
   x6 = x5.*x1
   x7 = x6.*x1
   #------------------------------------------------
   # Primaries
   #------------------------------------------------

# Primary invariants for NBody=5and deg=7
 # : what goes in the function for the evaluation
P1 = fpoly((x1,) , NB5I.P1)
P2 = fpoly((x1,x1,) , NB5I.P2)
P3 = fpoly((x2,) , NB5I.P3)
P4 = fpoly((x1,x1,x1,) , NB5I.P4)
P5 = fpoly((x3,) , NB5I.P5)
P6 = fpoly((x1,x1,x1,x1,) , NB5I.P6)
P7 = fpoly((x4,) , NB5I.P7)
P8 = fpoly((x1,x1,x1,x1,x1,) , NB5I.P8)
P9 = fpoly((x5,) , NB5I.P9)
P10 = fpoly((x6,) , NB5I.P10)

   #------------------------------------------------
   # Irreducible secondaries
   #------------------------------------------------

# Irreducible secondaries for NBody=5and deg=7
 # : what goes in the function for the evaluation
IS1 = fpoly((x2,x1,) , NB5I.IS1)
IS2 = fpoly((x1,x1,x1,) , NB5I.IS2)
IS3 = fpoly((x3,x1,) , NB5I.IS3)
IS4 = fpoly((x2,x2,) , NB5I.IS4)
IS5 = fpoly((x2,x1,x1,) , NB5I.IS5)
IS6 = fpoly((x2,x1,x1,) , NB5I.IS6)
IS7 = fpoly((x2,x1,x1,) , NB5I.IS7)
IS8 = fpoly((x4,x1,) , NB5I.IS8)
IS9 = fpoly((x3,x2,) , NB5I.IS9)
IS10 = fpoly((x3,x1,x1,) , NB5I.IS10)
IS11 = fpoly((x2,x2,x1,) , NB5I.IS11)
IS12 = fpoly((x3,x1,x1,) , NB5I.IS12)
IS13 = fpoly((x2,x2,x1,) , NB5I.IS13)
IS14 = fpoly((x3,x1,x1,) , NB5I.IS14)
IS15 = fpoly((x2,x1,x1,x1,) , NB5I.IS15)
IS16 = fpoly((x5,x1,) , NB5I.IS16)
IS17 = fpoly((x4,x2,) , NB5I.IS17)
IS18 = fpoly((x3,x3,) , NB5I.IS18)
IS19 = fpoly((x4,x1,x1,) , NB5I.IS19)
IS20 = fpoly((x3,x2,x1,) , NB5I.IS20)
IS21 = fpoly((x2,x2,x2,) , NB5I.IS21)
IS22 = fpoly((x4,x1,x1,) , NB5I.IS22)
IS23 = fpoly((x3,x2,x1,) , NB5I.IS23)
IS24 = fpoly((x4,x1,x1,) , NB5I.IS24)
IS25 = fpoly((x3,x1,x1,x1,) , NB5I.IS25)
IS26 = fpoly((x2,x2,x1,x1,) , NB5I.IS26)
IS27 = fpoly((x3,x2,x1,) , NB5I.IS27)
IS28 = fpoly((x7,) , NB5I.IS28)
IS29 = fpoly((x6,x1,) , NB5I.IS29)
IS30 = fpoly((x5,x2,) , NB5I.IS30)
IS31 = fpoly((x4,x3,) , NB5I.IS31)
IS32 = fpoly((x5,x1,x1,) , NB5I.IS32)
IS33 = fpoly((x4,x2,x1,) , NB5I.IS33)
IS34 = fpoly((x3,x3,x1,) , NB5I.IS34)
IS35 = fpoly((x3,x2,x2,) , NB5I.IS35)
IS36 = fpoly((x5,x1,x1,) , NB5I.IS36)
IS37 = fpoly((x4,x2,x1,) , NB5I.IS37)
IS38 = fpoly((x3,x3,x1,) , NB5I.IS38)
IS39 = fpoly((x5,x1,x1,) , NB5I.IS39)
IS40 = fpoly((x4,x1,x1,x1,) , NB5I.IS40)

   #------------------------------------------------
   # All secondaries
   #------------------------------------------------

 SEC1  = 1
 SEC2  = IS1
 SEC3  = IS2
 SEC4  = IS3
 SEC5  = IS4
 SEC6  = IS5
 SEC7  = IS6
 SEC8  = IS7
 SEC9  = IS8
 SEC10  = IS9
 SEC11  = IS10
 SEC12  = IS11
 SEC13  = IS12
 SEC14  = IS13
 SEC15  = IS14
 SEC16  = IS15
 SEC17  = IS1 *IS1
 SEC18  = IS1 *IS2
 SEC19  = IS2 *IS2
 SEC20  = IS16
 SEC21  = IS17
 SEC22  = IS18
 SEC23  = IS19
 SEC24  = IS20
 SEC25  = IS21
 SEC26  = IS22
 SEC27  = IS23
 SEC28  = IS24
 SEC29  = IS25
 SEC30  = IS26
 SEC31  = IS27
 SEC32  = IS1 *IS3
 SEC33  = IS2 *IS3
 SEC34  = IS1 *IS4
 SEC35  = IS2 *IS4
 SEC36  = IS1 *IS5
 SEC37  = IS2 *IS5
 SEC38  = IS1 *IS6
 SEC39  = IS2 *IS6
 SEC40  = IS1 *IS7
 SEC41  = IS2 *IS7
 SEC42  = IS28
 SEC43  = IS29
 SEC44  = IS30
 SEC45  = IS31
 SEC46  = IS32
 SEC47  = IS33
 SEC48  = IS34
 SEC49  = IS35
 SEC50  = IS36
 SEC51  = IS37
 SEC52  = IS38
 SEC53  = IS39
 SEC54  = IS40

return (@SVector [P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,]), (@SVector [SEC1,SEC2,SEC3,SEC4,SEC5,SEC6,SEC7,SEC8,SEC9,SEC10,SEC11,SEC12,SEC13,SEC14,SEC15,SEC16,SEC17,SEC18,SEC19,SEC20,SEC21,SEC22,SEC23,SEC24,SEC25,SEC26,SEC27,SEC28,SEC29,SEC30,SEC31,SEC32,SEC33,SEC34,SEC35,SEC36,SEC37,SEC38,SEC39,SEC40,SEC41,SEC42,SEC43,SEC44,SEC45,SEC46,SEC47,SEC48,SEC49,SEC50,SEC51,SEC52,SEC53,SEC54,])
 end

function invariants_d_gen(x1::SVector{10, T}) where {T}
   x2 = x1.*x1
   x3 = x2.*x1
   x4 = x3.*x1
   x5 = x4.*x1
   x6 = x5.*x1
   x7 = x6.*x1

   dx1 = @SVector ones(10)
   dx2 = 2 * x1
   dx3 = 3 * x2
   dx4 = 4 * x3
   dx5 = 5 * x4
   dx6 = 6 * x5
   dx7 = 7 * x6
   #------------------------------------------------
   # Primaries
   #------------------------------------------------

# Primary invariants for NBody=5and deg=7
 # : what goes in the function for the derivatives
dP1 = fpoly_d((x1,),(dx1,) , NB5I.P1)
dP2 = fpoly_d((x1,x1,),(dx1,dx1,) , NB5I.P2)
dP3 = fpoly_d((x2,),(dx2,) , NB5I.P3)
dP4 = fpoly_d((x1,x1,x1,),(dx1,dx1,dx1,) , NB5I.P4)
dP5 = fpoly_d((x3,),(dx3,) , NB5I.P5)
dP6 = fpoly_d((x1,x1,x1,x1,),(dx1,dx1,dx1,dx1,) , NB5I.P6)
dP7 = fpoly_d((x4,),(dx4,) , NB5I.P7)
dP8 = fpoly_d((x1,x1,x1,x1,x1,),(dx1,dx1,dx1,dx1,dx1,) , NB5I.P8)
dP9 = fpoly_d((x5,),(dx5,) , NB5I.P9)
dP10 = fpoly_d((x6,),(dx6,) , NB5I.P10)

   #------------------------------------------------
   # Irreducible secondaries
   #------------------------------------------------

# Irreducible secondaries for NBody=5and deg=7
 # : what goes in the function for the evaluation
IS1 = fpoly((x2,x1,) , NB5I.IS1)
IS2 = fpoly((x1,x1,x1,) , NB5I.IS2)
IS3 = fpoly((x3,x1,) , NB5I.IS3)
IS4 = fpoly((x2,x2,) , NB5I.IS4)
IS5 = fpoly((x2,x1,x1,) , NB5I.IS5)
IS6 = fpoly((x2,x1,x1,) , NB5I.IS6)
IS7 = fpoly((x2,x1,x1,) , NB5I.IS7)
IS8 = fpoly((x4,x1,) , NB5I.IS8)
IS9 = fpoly((x3,x2,) , NB5I.IS9)
IS10 = fpoly((x3,x1,x1,) , NB5I.IS10)
IS11 = fpoly((x2,x2,x1,) , NB5I.IS11)
IS12 = fpoly((x3,x1,x1,) , NB5I.IS12)
IS13 = fpoly((x2,x2,x1,) , NB5I.IS13)
IS14 = fpoly((x3,x1,x1,) , NB5I.IS14)
IS15 = fpoly((x2,x1,x1,x1,) , NB5I.IS15)
IS16 = fpoly((x5,x1,) , NB5I.IS16)
IS17 = fpoly((x4,x2,) , NB5I.IS17)
IS18 = fpoly((x3,x3,) , NB5I.IS18)
IS19 = fpoly((x4,x1,x1,) , NB5I.IS19)
IS20 = fpoly((x3,x2,x1,) , NB5I.IS20)
IS21 = fpoly((x2,x2,x2,) , NB5I.IS21)
IS22 = fpoly((x4,x1,x1,) , NB5I.IS22)
IS23 = fpoly((x3,x2,x1,) , NB5I.IS23)
IS24 = fpoly((x4,x1,x1,) , NB5I.IS24)
IS25 = fpoly((x3,x1,x1,x1,) , NB5I.IS25)
IS26 = fpoly((x2,x2,x1,x1,) , NB5I.IS26)
IS27 = fpoly((x3,x2,x1,) , NB5I.IS27)
IS28 = fpoly((x7,) , NB5I.IS28)
IS29 = fpoly((x6,x1,) , NB5I.IS29)
IS30 = fpoly((x5,x2,) , NB5I.IS30)
IS31 = fpoly((x4,x3,) , NB5I.IS31)
IS32 = fpoly((x5,x1,x1,) , NB5I.IS32)
IS33 = fpoly((x4,x2,x1,) , NB5I.IS33)
IS34 = fpoly((x3,x3,x1,) , NB5I.IS34)
IS35 = fpoly((x3,x2,x2,) , NB5I.IS35)
IS36 = fpoly((x5,x1,x1,) , NB5I.IS36)
IS37 = fpoly((x4,x2,x1,) , NB5I.IS37)
IS38 = fpoly((x3,x3,x1,) , NB5I.IS38)
IS39 = fpoly((x5,x1,x1,) , NB5I.IS39)
IS40 = fpoly((x4,x1,x1,x1,) , NB5I.IS40)

# Irreducible secondaries for NBody=5and deg=7
 # : what goes in the function for the derivatives
dIS1 = fpoly_d((x2,x1,),(dx2,dx1,) , NB5I.IS1)
dIS2 = fpoly_d((x1,x1,x1,),(dx1,dx1,dx1,) , NB5I.IS2)
dIS3 = fpoly_d((x3,x1,),(dx3,dx1,) , NB5I.IS3)
dIS4 = fpoly_d((x2,x2,),(dx2,dx2,) , NB5I.IS4)
dIS5 = fpoly_d((x2,x1,x1,),(dx2,dx1,dx1,) , NB5I.IS5)
dIS6 = fpoly_d((x2,x1,x1,),(dx2,dx1,dx1,) , NB5I.IS6)
dIS7 = fpoly_d((x2,x1,x1,),(dx2,dx1,dx1,) , NB5I.IS7)
dIS8 = fpoly_d((x4,x1,),(dx4,dx1,) , NB5I.IS8)
dIS9 = fpoly_d((x3,x2,),(dx3,dx2,) , NB5I.IS9)
dIS10 = fpoly_d((x3,x1,x1,),(dx3,dx1,dx1,) , NB5I.IS10)
dIS11 = fpoly_d((x2,x2,x1,),(dx2,dx2,dx1,) , NB5I.IS11)
dIS12 = fpoly_d((x3,x1,x1,),(dx3,dx1,dx1,) , NB5I.IS12)
dIS13 = fpoly_d((x2,x2,x1,),(dx2,dx2,dx1,) , NB5I.IS13)
dIS14 = fpoly_d((x3,x1,x1,),(dx3,dx1,dx1,) , NB5I.IS14)
dIS15 = fpoly_d((x2,x1,x1,x1,),(dx2,dx1,dx1,dx1,) , NB5I.IS15)
dIS16 = fpoly_d((x5,x1,),(dx5,dx1,) , NB5I.IS16)
dIS17 = fpoly_d((x4,x2,),(dx4,dx2,) , NB5I.IS17)
dIS18 = fpoly_d((x3,x3,),(dx3,dx3,) , NB5I.IS18)
dIS19 = fpoly_d((x4,x1,x1,),(dx4,dx1,dx1,) , NB5I.IS19)
dIS20 = fpoly_d((x3,x2,x1,),(dx3,dx2,dx1,) , NB5I.IS20)
dIS21 = fpoly_d((x2,x2,x2,),(dx2,dx2,dx2,) , NB5I.IS21)
dIS22 = fpoly_d((x4,x1,x1,),(dx4,dx1,dx1,) , NB5I.IS22)
dIS23 = fpoly_d((x3,x2,x1,),(dx3,dx2,dx1,) , NB5I.IS23)
dIS24 = fpoly_d((x4,x1,x1,),(dx4,dx1,dx1,) , NB5I.IS24)
dIS25 = fpoly_d((x3,x1,x1,x1,),(dx3,dx1,dx1,dx1,) , NB5I.IS25)
dIS26 = fpoly_d((x2,x2,x1,x1,),(dx2,dx2,dx1,dx1,) , NB5I.IS26)
dIS27 = fpoly_d((x3,x2,x1,),(dx3,dx2,dx1,) , NB5I.IS27)
dIS28 = fpoly_d((x7,),(dx7,) , NB5I.IS28)
dIS29 = fpoly_d((x6,x1,),(dx6,dx1,) , NB5I.IS29)
dIS30 = fpoly_d((x5,x2,),(dx5,dx2,) , NB5I.IS30)
dIS31 = fpoly_d((x4,x3,),(dx4,dx3,) , NB5I.IS31)
dIS32 = fpoly_d((x5,x1,x1,),(dx5,dx1,dx1,) , NB5I.IS32)
dIS33 = fpoly_d((x4,x2,x1,),(dx4,dx2,dx1,) , NB5I.IS33)
dIS34 = fpoly_d((x3,x3,x1,),(dx3,dx3,dx1,) , NB5I.IS34)
dIS35 = fpoly_d((x3,x2,x2,),(dx3,dx2,dx2,) , NB5I.IS35)
dIS36 = fpoly_d((x5,x1,x1,),(dx5,dx1,dx1,) , NB5I.IS36)
dIS37 = fpoly_d((x4,x2,x1,),(dx4,dx2,dx1,) , NB5I.IS37)
dIS38 = fpoly_d((x3,x3,x1,),(dx3,dx3,dx1,) , NB5I.IS38)
dIS39 = fpoly_d((x5,x1,x1,),(dx5,dx1,dx1,) , NB5I.IS39)
dIS40 = fpoly_d((x4,x1,x1,x1,),(dx4,dx1,dx1,dx1,) , NB5I.IS40)

   #------------------------------------------------
   # All secondaries
   #------------------------------------------------

 dSEC1   = @SVector zeros(10)
 dSEC2  = dIS1
 dSEC3  = dIS2
 dSEC4  = dIS3
 dSEC5  = dIS4
 dSEC6  = dIS5
 dSEC7  = dIS6
 dSEC8  = dIS7
 dSEC9  = dIS8
 dSEC10  = dIS9
 dSEC11  = dIS10
 dSEC12  = dIS11
 dSEC13  = dIS12
 dSEC14  = dIS13
 dSEC15  = dIS14
 dSEC16  = dIS15
 dSEC17   =  dIS1 *IS1 + IS1 *dIS1
 dSEC18   =  dIS1 *IS2 + IS1 *dIS2
 dSEC19   =  dIS2 *IS2 + IS2 *dIS2
 dSEC20  = dIS16
 dSEC21  = dIS17
 dSEC22  = dIS18
 dSEC23  = dIS19
 dSEC24  = dIS20
 dSEC25  = dIS21
 dSEC26  = dIS22
 dSEC27  = dIS23
 dSEC28  = dIS24
 dSEC29  = dIS25
 dSEC30  = dIS26
 dSEC31  = dIS27
 dSEC32   =  dIS1 *IS3 + IS1 *dIS3
 dSEC33   =  dIS2 *IS3 + IS2 *dIS3
 dSEC34   =  dIS1 *IS4 + IS1 *dIS4
 dSEC35   =  dIS2 *IS4 + IS2 *dIS4
 dSEC36   =  dIS1 *IS5 + IS1 *dIS5
 dSEC37   =  dIS2 *IS5 + IS2 *dIS5
 dSEC38   =  dIS1 *IS6 + IS1 *dIS6
 dSEC39   =  dIS2 *IS6 + IS2 *dIS6
 dSEC40   =  dIS1 *IS7 + IS1 *dIS7
 dSEC41   =  dIS2 *IS7 + IS2 *dIS7
 dSEC42  = dIS28
 dSEC43  = dIS29
 dSEC44  = dIS30
 dSEC45  = dIS31
 dSEC46  = dIS32
 dSEC47  = dIS33
 dSEC48  = dIS34
 dSEC49  = dIS35
 dSEC50  = dIS36
 dSEC51  = dIS37
 dSEC52  = dIS38
 dSEC53  = dIS39
 dSEC54  = dIS40

return (@SVector [dP1,dP2,dP3,dP4,dP5,dP6,dP7,dP8,dP9,dP10,]),
       (@SVector [dSEC1,dSEC2,dSEC3,dSEC4,dSEC5,dSEC6,dSEC7,dSEC8,dSEC9,dSEC10,dSEC11,dSEC12,dSEC13,dSEC14,dSEC15,dSEC16,dSEC17,dSEC18,dSEC19,dSEC20,dSEC21,dSEC22,dSEC23,dSEC24,dSEC25,dSEC26,dSEC27,dSEC28,dSEC29,dSEC30,dSEC31,dSEC32,dSEC33,dSEC34,dSEC35,dSEC36,dSEC37,dSEC38,dSEC39,dSEC40,dSEC41,dSEC42,dSEC43,dSEC44,dSEC45,dSEC46,dSEC47,dSEC48,dSEC49,dSEC50,dSEC51,dSEC52,dSEC53,dSEC54,])

 end

function invariants_ed_gen(x1::SVector{10, T}) where {T}
   x2 = x1.*x1
   x3 = x2.*x1
   x4 = x3.*x1
   x5 = x4.*x1
   x6 = x5.*x1
   x7 = x6.*x1

   dx1 = @SVector ones(10)
   dx2 = 2 * x1
   dx3 = 3 * x2
   dx4 = 4 * x3
   dx5 = 5 * x4
   dx6 = 6 * x5
   dx7 = 7 * x6
   #------------------------------------------------
   # Primaries
   #------------------------------------------------

# Primary invariants for NBody=5and deg=7
 # : what goes in the function for the evaluation and derivatives
P1, dP1 = fpoly_ed((x1,),(dx1,) , NB5I.P1)
P2, dP2 = fpoly_ed((x1,x1,),(dx1,dx1,) , NB5I.P2)
P3, dP3 = fpoly_ed((x2,),(dx2,) , NB5I.P3)
P4, dP4 = fpoly_ed((x1,x1,x1,),(dx1,dx1,dx1,) , NB5I.P4)
P5, dP5 = fpoly_ed((x3,),(dx3,) , NB5I.P5)
P6, dP6 = fpoly_ed((x1,x1,x1,x1,),(dx1,dx1,dx1,dx1,) , NB5I.P6)
P7, dP7 = fpoly_ed((x4,),(dx4,) , NB5I.P7)
P8, dP8 = fpoly_ed((x1,x1,x1,x1,x1,),(dx1,dx1,dx1,dx1,dx1,) , NB5I.P8)
P9, dP9 = fpoly_ed((x5,),(dx5,) , NB5I.P9)
P10, dP10 = fpoly_ed((x6,),(dx6,) , NB5I.P10)

   #------------------------------------------------
   # Irreducible secondaries
   #------------------------------------------------

# Irreducible secondaries for NBody=5and deg=7
 # : what goes in the function for the evaluation and derivatives
IS1, dIS1 = fpoly_ed((x2,x1,),(dx2,dx1,) , NB5I.IS1)
IS2, dIS2 = fpoly_ed((x1,x1,x1,),(dx1,dx1,dx1,) , NB5I.IS2)
IS3, dIS3 = fpoly_ed((x3,x1,),(dx3,dx1,) , NB5I.IS3)
IS4, dIS4 = fpoly_ed((x2,x2,),(dx2,dx2,) , NB5I.IS4)
IS5, dIS5 = fpoly_ed((x2,x1,x1,),(dx2,dx1,dx1,) , NB5I.IS5)
IS6, dIS6 = fpoly_ed((x2,x1,x1,),(dx2,dx1,dx1,) , NB5I.IS6)
IS7, dIS7 = fpoly_ed((x2,x1,x1,),(dx2,dx1,dx1,) , NB5I.IS7)
IS8, dIS8 = fpoly_ed((x4,x1,),(dx4,dx1,) , NB5I.IS8)
IS9, dIS9 = fpoly_ed((x3,x2,),(dx3,dx2,) , NB5I.IS9)
IS10, dIS10 = fpoly_ed((x3,x1,x1,),(dx3,dx1,dx1,) , NB5I.IS10)
IS11, dIS11 = fpoly_ed((x2,x2,x1,),(dx2,dx2,dx1,) , NB5I.IS11)
IS12, dIS12 = fpoly_ed((x3,x1,x1,),(dx3,dx1,dx1,) , NB5I.IS12)
IS13, dIS13 = fpoly_ed((x2,x2,x1,),(dx2,dx2,dx1,) , NB5I.IS13)
IS14, dIS14 = fpoly_ed((x3,x1,x1,),(dx3,dx1,dx1,) , NB5I.IS14)
IS15, dIS15 = fpoly_ed((x2,x1,x1,x1,),(dx2,dx1,dx1,dx1,) , NB5I.IS15)
IS16, dIS16 = fpoly_ed((x5,x1,),(dx5,dx1,) , NB5I.IS16)
IS17, dIS17 = fpoly_ed((x4,x2,),(dx4,dx2,) , NB5I.IS17)
IS18, dIS18 = fpoly_ed((x3,x3,),(dx3,dx3,) , NB5I.IS18)
IS19, dIS19 = fpoly_ed((x4,x1,x1,),(dx4,dx1,dx1,) , NB5I.IS19)
IS20, dIS20 = fpoly_ed((x3,x2,x1,),(dx3,dx2,dx1,) , NB5I.IS20)
IS21, dIS21 = fpoly_ed((x2,x2,x2,),(dx2,dx2,dx2,) , NB5I.IS21)
IS22, dIS22 = fpoly_ed((x4,x1,x1,),(dx4,dx1,dx1,) , NB5I.IS22)
IS23, dIS23 = fpoly_ed((x3,x2,x1,),(dx3,dx2,dx1,) , NB5I.IS23)
IS24, dIS24 = fpoly_ed((x4,x1,x1,),(dx4,dx1,dx1,) , NB5I.IS24)
IS25, dIS25 = fpoly_ed((x3,x1,x1,x1,),(dx3,dx1,dx1,dx1,) , NB5I.IS25)
IS26, dIS26 = fpoly_ed((x2,x2,x1,x1,),(dx2,dx2,dx1,dx1,) , NB5I.IS26)
IS27, dIS27 = fpoly_ed((x3,x2,x1,),(dx3,dx2,dx1,) , NB5I.IS27)
IS28, dIS28 = fpoly_ed((x7,),(dx7,) , NB5I.IS28)
IS29, dIS29 = fpoly_ed((x6,x1,),(dx6,dx1,) , NB5I.IS29)
IS30, dIS30 = fpoly_ed((x5,x2,),(dx5,dx2,) , NB5I.IS30)
IS31, dIS31 = fpoly_ed((x4,x3,),(dx4,dx3,) , NB5I.IS31)
IS32, dIS32 = fpoly_ed((x5,x1,x1,),(dx5,dx1,dx1,) , NB5I.IS32)
IS33, dIS33 = fpoly_ed((x4,x2,x1,),(dx4,dx2,dx1,) , NB5I.IS33)
IS34, dIS34 = fpoly_ed((x3,x3,x1,),(dx3,dx3,dx1,) , NB5I.IS34)
IS35, dIS35 = fpoly_ed((x3,x2,x2,),(dx3,dx2,dx2,) , NB5I.IS35)
IS36, dIS36 = fpoly_ed((x5,x1,x1,),(dx5,dx1,dx1,) , NB5I.IS36)
IS37, dIS37 = fpoly_ed((x4,x2,x1,),(dx4,dx2,dx1,) , NB5I.IS37)
IS38, dIS38 = fpoly_ed((x3,x3,x1,),(dx3,dx3,dx1,) , NB5I.IS38)
IS39, dIS39 = fpoly_ed((x5,x1,x1,),(dx5,dx1,dx1,) , NB5I.IS39)
IS40, dIS40 = fpoly_ed((x4,x1,x1,x1,),(dx4,dx1,dx1,dx1,) , NB5I.IS40)

   #------------------------------------------------
   # All secondaries
   #------------------------------------------------

 SEC1  = 1
 SEC2  = IS1
 SEC3  = IS2
 SEC4  = IS3
 SEC5  = IS4
 SEC6  = IS5
 SEC7  = IS6
 SEC8  = IS7
 SEC9  = IS8
 SEC10  = IS9
 SEC11  = IS10
 SEC12  = IS11
 SEC13  = IS12
 SEC14  = IS13
 SEC15  = IS14
 SEC16  = IS15
 SEC17  = IS1 *IS1
 SEC18  = IS1 *IS2
 SEC19  = IS2 *IS2
 SEC20  = IS16
 SEC21  = IS17
 SEC22  = IS18
 SEC23  = IS19
 SEC24  = IS20
 SEC25  = IS21
 SEC26  = IS22
 SEC27  = IS23
 SEC28  = IS24
 SEC29  = IS25
 SEC30  = IS26
 SEC31  = IS27
 SEC32  = IS1 *IS3
 SEC33  = IS2 *IS3
 SEC34  = IS1 *IS4
 SEC35  = IS2 *IS4
 SEC36  = IS1 *IS5
 SEC37  = IS2 *IS5
 SEC38  = IS1 *IS6
 SEC39  = IS2 *IS6
 SEC40  = IS1 *IS7
 SEC41  = IS2 *IS7
 SEC42  = IS28
 SEC43  = IS29
 SEC44  = IS30
 SEC45  = IS31
 SEC46  = IS32
 SEC47  = IS33
 SEC48  = IS34
 SEC49  = IS35
 SEC50  = IS36
 SEC51  = IS37
 SEC52  = IS38
 SEC53  = IS39
 SEC54  = IS40

 dSEC1   = @SVector zeros(10)
 dSEC2  = dIS1
 dSEC3  = dIS2
 dSEC4  = dIS3
 dSEC5  = dIS4
 dSEC6  = dIS5
 dSEC7  = dIS6
 dSEC8  = dIS7
 dSEC9  = dIS8
 dSEC10  = dIS9
 dSEC11  = dIS10
 dSEC12  = dIS11
 dSEC13  = dIS12
 dSEC14  = dIS13
 dSEC15  = dIS14
 dSEC16  = dIS15
 dSEC17   =  dIS1 *IS1 + IS1 *dIS1
 dSEC18   =  dIS1 *IS2 + IS1 *dIS2
 dSEC19   =  dIS2 *IS2 + IS2 *dIS2
 dSEC20  = dIS16
 dSEC21  = dIS17
 dSEC22  = dIS18
 dSEC23  = dIS19
 dSEC24  = dIS20
 dSEC25  = dIS21
 dSEC26  = dIS22
 dSEC27  = dIS23
 dSEC28  = dIS24
 dSEC29  = dIS25
 dSEC30  = dIS26
 dSEC31  = dIS27
 dSEC32   =  dIS1 *IS3 + IS1 *dIS3
 dSEC33   =  dIS2 *IS3 + IS2 *dIS3
 dSEC34   =  dIS1 *IS4 + IS1 *dIS4
 dSEC35   =  dIS2 *IS4 + IS2 *dIS4
 dSEC36   =  dIS1 *IS5 + IS1 *dIS5
 dSEC37   =  dIS2 *IS5 + IS2 *dIS5
 dSEC38   =  dIS1 *IS6 + IS1 *dIS6
 dSEC39   =  dIS2 *IS6 + IS2 *dIS6
 dSEC40   =  dIS1 *IS7 + IS1 *dIS7
 dSEC41   =  dIS2 *IS7 + IS2 *dIS7
 dSEC42  = dIS28
 dSEC43  = dIS29
 dSEC44  = dIS30
 dSEC45  = dIS31
 dSEC46  = dIS32
 dSEC47  = dIS33
 dSEC48  = dIS34
 dSEC49  = dIS35
 dSEC50  = dIS36
 dSEC51  = dIS37
 dSEC52  = dIS38
 dSEC53  = dIS39
 dSEC54  = dIS40

return (@SVector [P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,]),
       (@SVector [SEC1,SEC2,SEC3,SEC4,SEC5,SEC6,SEC7,SEC8,SEC9,SEC10,SEC11,SEC12,SEC13,SEC14,SEC15,SEC16,SEC17,SEC18,SEC19,SEC20,SEC21,SEC22,SEC23,SEC24,SEC25,SEC26,SEC27,SEC28,SEC29,SEC30,SEC31,SEC32,SEC33,SEC34,SEC35,SEC36,SEC37,SEC38,SEC39,SEC40,SEC41,SEC42,SEC43,SEC44,SEC45,SEC46,SEC47,SEC48,SEC49,SEC50,SEC51,SEC52,SEC53,SEC54,]),
       (@SVector [dP1,dP2,dP3,dP4,dP5,dP6,dP7,dP8,dP9,dP10,]),
       (@SVector [dSEC1,dSEC2,dSEC3,dSEC4,dSEC5,dSEC6,dSEC7,dSEC8,dSEC9,dSEC10,dSEC11,dSEC12,dSEC13,dSEC14,dSEC15,dSEC16,dSEC17,dSEC18,dSEC19,dSEC20,dSEC21,dSEC22,dSEC23,dSEC24,dSEC25,dSEC26,dSEC27,dSEC28,dSEC29,dSEC30,dSEC31,dSEC32,dSEC33,dSEC34,dSEC35,dSEC36,dSEC37,dSEC38,dSEC39,dSEC40,dSEC41,dSEC42,dSEC43,dSEC44,dSEC45,dSEC46,dSEC47,dSEC48,dSEC49,dSEC50,dSEC51,dSEC52,dSEC53,dSEC54,])
 end

end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Filters data of 'n_before' seconds before CPA and all data after CPA.  If CPA occurs prior to 't_min',
sample uniformly between 't_min' and 't_max'.
"""
module FilterNMACInfo

export script_dasc, script_libcas098small, script_libcas0100star, script_libcas0100llcem
export remove_cpa, find_cpa, filter_nmac_info

using FilterRecords
using DataFrameSets
using DataFrames
using Datasets

##############################
#Add an entry here

const DASC_FILTERED = Pkg.dir("Datasets/data/dascfilt")
const LIBCAS098SMALL_FILTERED = Pkg.dir("Datasets/data/libcas098smallfilt")
const LIBCAS0100STAR_FILTERED = Pkg.dir("Datasets/data/libcas0100starfilt")
const LIBCAS0100LLCEM_FILTERED = Pkg.dir("Datasets/data/libcas0100llcemfilt")

script_dasc() = remove_cpa("dasc", DASC_FILTERED; t_min=35, n_before=5)
script_libcas098small() = remove_cpa("libcas098small", LIBCAS098SMALL_FILTERED; t_min=35, n_before=5)
script_libcas0100star() = remove_cpa("libcas0100star", LIBCAS0100STAR_FILTERED; t_min=35, n_before=5)
script_libcas0100llcem() = remove_cpa("libcas0100llcem", LIBCAS0100LLCEM_FILTERED; t_min=35, n_before=5)

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

function cpa_metric(r::DataFrameRow)
  return r[:abs_alt_diff] + r[:horizontal_range] / 5
end

#returns row of cpa
function find_cpa(D::DataFrame)
  minval = realmax(Float64)
  minrow = 0
  for (i, r) in enumerate(eachrow(D))
    val = cpa_metric(r)
    if val < minval
      minval = val
      minrow = i
    end
  end
  return minrow
end

function remove_cpa(dataname::AbstractString, outdir::AbstractString; kwargs...)
  Ds = dataset(dataname)
  remove_cpa!(Ds; kwargs...)
  save_csvs(outdir, Ds)
end

function remove_cpa!(Ds::DFSet; kwargs...)
  for D in getrecords(Ds)
    remove_cpa!(D; kwargs...)
  end
end

function remove_cpa!(D::DataFrame; t_min::Int64=35, n_before::Int64=5,
    rng::AbstractRNG=MersenneTwister(0))
  i_cpa = find_cpa(D)

  #not in acceptable range, select randomly
  if !(t_min <= i_cpa <= nrow(D))
    i_cpa = rand(rng, t_min:nrow(D))
  end
  deleterows!(D, (i_cpa - n_before):nrow(D))
end

############
#deprecated
function isnmac(row::DataFrameRow)
  return row[:horizontal_range] <= 500.0 && row[:abs_alt_diff] <= 100.0
end

function filter_nmac_info(f::Function, dataname::AbstractString, outdir::AbstractString)
  Ds = dataset(dataname)
  filter_by_bool!(f, Ds, n_before=5, n_after=typemax(Int64))

  save_csvs(outdir, Ds)
end

#dasc_script() = filter_nmac_info(isnmac, "dasc", DASC_FILTERED)
#libcas098small_script() = filter_nmac_info(isnmac, "libcas098small", LIBCAS098SMALL_FILTERED)
############

end #module


text:
@inline function initialize!(integrator,cache::StochasticCompositeCache,f=integrator.f)
  cache.current = cache.choice_function(integrator)
  initialize!(integrator,cache.caches[cache.current])
end

@inline function perform_step!(integrator,cache::StochasticCompositeCache,f=integrator.f)
  perform_step!(integrator,cache.caches[cache.current],f)
end

@inline choose_algorithm!(integrator,cache::StochasticDiffEqCache) = nothing
@inline function choose_algorithm!(integrator,cache::StochasticCompositeCache)
  new_current = cache.choice_function(integrator)
  if new_current != cache.current
    initialize!(integrator,cache.caches[new_current])
    reset_alg_dependent_opts!(integrator,integrator.alg.algs[cache.current],integrator.alg.algs[new_current])
    transfer_cache!(integrator,integrator.cache.caches[cache.current],integrator.cache.caches[new_current])
    cache.current = new_current
  end
end

"""
If no user default, then this will change the default to the defaults
for the second algorithm.
"""
@inline function reset_alg_dependent_opts!(integrator,alg1,alg2)
  integrator.dtchangeable = isdtchangeable(alg2)
  if integrator.opts.adaptive == isadaptive(alg1)
    integrator.opts.adaptive = isadaptive(alg2)
  end
  if integrator.opts.qmin == qmin_default(alg1)
    integrator.opts.qmin = qmin_default(alg2)
  end
  if integrator.opts.qmax == qmax_default(alg1)
    integrator.opts.qmax == qmax_default(alg2)
  end
  if integrator.opts.beta2 == beta2_default(alg1)
    integrator.opts.beta2 = beta2_default(alg2)
  end
  if integrator.opts.beta1 == beta1_default(alg1,integrator.opts.beta2)
    integrator.opts.beta1 = beta1_default(alg2,integrator.opts.beta2)
  end
end

# Write how to transfer the cache variables from one cache to the other
# Example: send the history variables from one multistep method to another

@inline transfer_cache!(integrator,alg1,alg2) = nothing


text:
module mmesh1
using LibFEM
using LibFEM.Mesh
using Test
function test()
    @test Mesh.P1.manifdim == 0
    @test Mesh.L2.nfacets == 2
    true
end
end
using .mmesh1
mmesh1.test()

module mmesh2
using StaticArrays
using LibFEM
using LibFEM.Mesh: Vertices, nvertices, coordinates
using Test
function test()
    xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0; 1900.0 0.0; 0.0 400.0; 633.3333333333334 400.0; 1266.6666666666667 400.0; 1900.0 400.0; 0.0 800.0; 633.3333333333334 800.0; 1266.6666666666667 800.0; 1900.0 800.0]
    v =  Vertices(xyz)
    @test nvertices(v) == 12
    x = coordinates(v, SVector{2}([2, 4]))
    @test x[1] == SVector{2}([633.3333333333334 0.0])
    true
end
end
using .mmesh2
mmesh2.test()

module mmesh3
using StaticArrays
using LibFEM
using LibFEM.Mesh: L2, Q4, Shapes, connectivity, manifdim, nnodes, nfacets, facettype
using Test
function test()
    shapetype = Q4
    c = [(1, 2, 6, 5), (5, 6, 10, 9), (2, 3, 7, 6), (6, 7, 11, 10), (3, 4, 8, 7), (7, 8, 12, 11)]
    cc = [SVector{shapetype.nnodes}(c[idx]) for idx in 1:length(c)]
    shapes = Shapes(shapetype, cc)
    @test connectivity(shapes, 3)[3] == 7
    @test connectivity(shapes, SVector{2}([2, 4]))[1][4] == 9
    @test manifdim(shapes) == 2
    @test nnodes(shapes) == 4
    @test facettype(shapes) == L2
    @test nfacets(shapes) == 4
    true
end
end
using .mmesh3
mmesh3.test()


text:
export code_spir, code_hsa, src_hsail

function to_tuple_type(t::ANY)
    if isa(t,Tuple) || isa(t,AbstractArray) || isa(t,SimpleVector)
        t = Tuple{t...}
    end
    if isa(t,Type) && t<:Tuple
        if !all(p->(isa(p,Type)||isa(p,TypeVar)), t.parameters)
            error("argument tuple type must contain only types")
        end
    else
        error("expected tuple type")
    end
    t
end

# display ir of the generated spir function
function src_spir(f::Function, types::ANY)
    t = to_tuple_type(types)
    llvmf = ccall(:jl_get_spirf, Ptr{Void}, (Any, Any), f, t)

    if llvmf == C_NULL
        error("no method found for the specified argument types")
    end

    str = ccall(:jl_dump_function_ir, Any, (Ptr{Void},), llvmf)::ByteString

    return str
end

code_spir(io::IO, f::Function, types::ANY) = print(io, src_spir(f,types))

code_spir(f::ANY, t::ANY) = code_spir(STDOUT, f, t)

# display ir after linking to the hsa intrinsics
function src_hsa(f::Function, types::ANY)
    t = to_tuple_type(types)
    llvmf = ccall(:jl_get_hsailf, Ptr{Void}, (Any, Any, Bool), f, t, false)

    if llvmf == C_NULL
        error("no method found for the specified argument types")
    end

    str = ccall(:jl_dump_function_ir, Any, (Ptr{Void},), llvmf)::ByteString

    return str
end

code_hsa(io::IO, f::Function, types::ANY) = print(io, src_hsa(f,types))

code_hsa(f::ANY, t::ANY) = code_hsa(STDOUT, f, t)

# output final hsail code
function src_hsail(f::Function, types::ANY)
    t = to_tuple_type(types)
    llvmf = ccall(:jl_get_hsailf, Ptr{Void}, (Any, Any), f, t)

    if llvmf == C_NULL
        error("no method found for the specified argument types")
    end

    str = ccall(:jl_dump_function_hsail, Any, (Ptr{Void},), llvmf)::ByteString

    return str
end


text:
### GUI controls for play forward/backward, up/down, and frame stepping ###

module Navigation

using Tk, Compat

## Type for holding GUI state
# This specifies a particular 2d slice from a possibly-4D image
type NavigationState
    # Dimensions:
    zmax::Int          # number of frames in z, set to 1 if only 2 spatial dims
    tmax::Int          # number of frames in t, set to 1 if only a single image
    z::Int             # current position in z-stack
    t::Int             # current moment in time
    # Other state data:
    timer              # nothing if not playing, Timer if we are
    fps::Float64       # playback speed in frames per second
end

NavigationState(zmax::Integer, tmax::Integer, z::Integer, t::Integer) = NavigationState(@compat(Int(zmax)), @compat(Int(tmax)), @compat(Int(z)), @compat(Int(t)), nothing, 30.0)
NavigationState(zmax::Integer, tmax::Integer) = NavigationState(zmax, tmax, 1, 1)

function stop_playing!(state::NavigationState)
    if state.timer != nothing
        close(state.timer)
        state.timer = nothing
    end
end

## Type for holding "handles" to GUI controls
type NavigationControls
    stepup                            # z buttons...
    stepdown
    playup
    playdown
    stepback                          # t buttons...
    stepfwd
    playback
    playfwd
    stop
    editz                             # edit boxes
    editt
    textz                             # static text (information)
    textt
    scalez                            # scale (slider) widgets
    scalet
end
NavigationControls() = NavigationControls(nothing, nothing, nothing, nothing,
                                          nothing, nothing, nothing, nothing,
                                          nothing, nothing, nothing, nothing,
                                          nothing, nothing, nothing)

# f is a TkFrame
function init_navigation!(f, ctrls::NavigationControls, state::NavigationState, showframe::Function)
    btnsz, pad = widget_size()
    stop = trues(btnsz)
    mask = copy(stop)
    stop[[1,btnsz[1]],:] = false
    stop[:,[1,btnsz[2]]] = false
    bkg = "gray70"
    icon = Tk.Image(stop, mask, bkg, "black")
    ctrls.stop = Button(f, icon)
    bind(ctrls.stop, "command", path -> stop_playing!(state))
    local zindex
    local tindex
    local stopindex
    havez = state.zmax > 1
    havet = state.tmax > 1
    zindex = 1:6
    stopindex = 7
    tindex = 8:13
    if !havez
        stopindex = 1
        tindex = 2:7
    end
    grid(ctrls.stop, 1, stopindex, padx=3*pad, pady=pad)
    win = toplevel(f)
    if havez || havet
        bind(win, "<space>", path->stop_playing!(state))
    end
    if havez
        callback = (path->stepz(1,ctrls,state,showframe), path->playz(1,ctrls,state,showframe),
            path->playz(-1,ctrls,state,showframe), path->stepz(-1,ctrls,state,showframe),
            path->setz(ctrls,state,showframe), path->scalez(ctrls,state,showframe))
        ctrls.stepup, ctrls.playup, ctrls.playdown, ctrls.stepdown, ctrls.textz, ctrls.editz, ctrls.scalez =
            addbuttons(f, btnsz, bkg, pad, zindex, "z", callback, 1:state.zmax)
        bind(win, "<Alt-Up>", path->stepz(1,ctrls,state,showframe))
        bind(win, "<Alt-Down>", path->stepz(-1,ctrls,state,showframe))
        bind(win, "<Alt-Shift-Up>", path->playz(1,ctrls,state,showframe))
        bind(win, "<Alt-Shift-Down>", path->playz(-1,ctrls,state,showframe))
        updatez(ctrls, state)
    end
    if havet
        callback = (path->stept(-1,ctrls,state,showframe), path->playt(-1,ctrls,state,showframe),
            path->playt(1,ctrls,state,showframe), path->stept(1,ctrls,state,showframe),
            path->sett(ctrls,state,showframe), path->scalet(ctrls,state,showframe))
        ctrls.stepback, ctrls.playback, ctrls.playfwd, ctrls.stepfwd, ctrls.textt, ctrls.editt, ctrls.scalet =
            addbuttons(f, btnsz, bkg, pad, tindex, "t", callback, 1:state.tmax)
        bind(win, "<Alt-Right>", path->stept(1,ctrls,state,showframe))
        bind(win, "<Alt-Left>", path->stept(-1,ctrls,state,showframe))
        bind(win, "<Alt-Shift-Right>", path->playt(1,ctrls,state,showframe))
        bind(win, "<Alt-Shift-Left>", path->playt(-1,ctrls,state,showframe))
        updatet(ctrls, state)
    end
    # Context menu for settings
    menu = Menu(f)
    menu_fps = menu_add(menu, "Playback speed...", path -> set_fps!(state))
    tk_popup(f, menu)
end

# GUI to set the frame rate
function set_fps!(state::NavigationState)
    win = Toplevel()
    f = Frame(win)
    pack(f, expand=true, fill="both")

    l = Label(f, "Frames per second:")
    e = Entry(f, width=5)
    set_value(e, string(state.fps))
    ok = Button(f, "OK")
    cancel = Button(f, "Cancel")

    grid(l, 1, 1)
    grid(e, 1, 2, pady=5)
    grid(cancel, 2, 1)
    grid(ok, 2, 2)

    function set_close!(state::NavigationState)
        try
            fps = float64(get_value(e))
            state.fps = fps
            destroy(win)
        catch
            set_value(e, string(state.fps))
        end
    end
    bind(e, "<Return>", path->set_close!(state))
    bind(ok, "command", path->set_close!(state))
    bind(cancel, "command", path->destroy(win))
end

function widget_size()
    btnsz = OS_NAME == :Darwin ? (13, 13) : (21, 21)
    pad = 5
    return btnsz, pad
end

# Functions for drawing icons
function arrowheads(sz, vert::Bool)
    datasm = icondata(sz, 0.5)
    datalg = icondata(sz, 0.8)
    if vert
        return datasm[:,end:-1:1], datalg[:,end:-1:1], datalg, datasm
    else
        datasm = datasm'
        datalg = datalg'
        return datasm[end:-1:1,:], datalg[end:-1:1,:], datalg, datasm
    end
end

function icondata(iconsize, frac)
    center = ceil(Integer, iconsize[1]/2)
    data = Bool[ 2abs(i-center)< iconsize[2]-(j-1)/frac for i = 1:iconsize[1], j = 1:iconsize[2] ]
    data .== true
end

# index contains the grid position of each object
# orientation is "t" or "z"
# callback is an array of 5 entries, the 5th being the edit box
function addbuttons(f, sz, bkg, pad, index, orientation, callback, rng)
    rotflag = orientation == "z"
    ctrl = Array(Any, 7)
    ctrl[1], ctrl[2], ctrl[3], ctrl[4] = arrowheads(sz, rotflag)
    mask = trues(sz)
    const color = ("black", "green", "green", "black")
    ibutton = [1,2,5,6]
    for i = 1:4
        icon = Tk.Image(ctrl[i], mask, bkg, color[i])
        b = Button(f, icon)
        grid(b, 1, index[ibutton[i]], padx=pad, pady=pad)
        bind(b, "command", callback[i])
        ctrl[i] = b
    end
    ctrl[5] = Label(f, orientation*":")
    grid(ctrl[5], 1, index[3], padx=pad, pady=pad)
    ctrl[6] = Entry(f, "1")
    configure(ctrl[6], width=5)
    grid(ctrl[6], 1, index[4], padx=pad, pady=pad)
    bind(ctrl[6], "<Return>", callback[5])
    ctrl[7] = Slider(f, rng)
    grid(ctrl[7], 2, index, sticky="ew", padx=pad)
    bind(ctrl[7], "command", callback[6])
    tuple(ctrl...)
end

function updatez(ctrls, state)
    set_value(ctrls.editz, string(state.z))
    set_value(ctrls.scalez, state.z)
    enabledown = state.z > 1
    set_enabled(ctrls.stepdown, enabledown)
    set_enabled(ctrls.playdown, enabledown)
    enableup = state.z < state.zmax
    set_enabled(ctrls.stepup, enableup)
    set_enabled(ctrls.playup, enableup)
end

function updatet(ctrls, state)
    set_value(ctrls.editt, string(state.t))
    set_value(ctrls.scalet, state.t)
    enableback = state.t > 1
    set_enabled(ctrls.stepback, enableback)
    set_enabled(ctrls.playback, enableback)
    enablefwd = state.t < state.tmax
    set_enabled(ctrls.stepfwd, enablefwd)
    set_enabled(ctrls.playfwd, enablefwd)
end

function incrementz(inc, ctrls, state, showframe)
    state.z += inc
    updatez(ctrls, state)
    showframe(state)
end

function stepz(inc, ctrls, state, showframe)
    if 1 <= state.z+inc <= state.zmax
        incrementz(inc, ctrls, state, showframe)
    else
        stop_playing!(state)
    end
end

function playz(inc, ctrls, state, showframe)
    if !(state.fps > 0)
        error("Frame rate is not positive")
    end
    stop_playing!(state)
    dt = 1/state.fps
    state.timer = Timer(timer -> stepz(inc, ctrls, state, showframe), dt, dt)
end

function setz(ctrls,state, showframe)
    zstr = get_value(ctrls.editz)
    try
        val = parse(Int, zstr)
        state.z = val
        updatez(ctrls, state)
        showframe(state)
    catch
        updatez(ctrls, state)
    end
end

function scalez(ctrls, state, showframe)
    state.z = round(Int, get_value(ctrls.scalez))
    updatez(ctrls, state)
    showframe(state)
end

function incrementt(inc, ctrls, state, showframe)
    state.t += inc
    updatet(ctrls, state)
    showframe(state)
end

function stept(inc, ctrls, state, showframe)
    if 1 <= state.t+inc <= state.tmax
        incrementt(inc, ctrls, state, showframe)
    else
        stop_playing!(state)
    end
end

function playt(inc, ctrls, state, showframe)
    if !(state.fps > 0)
        error("Frame rate is not positive")
    end
    stop_playing!(state)
    dt = 1/state.fps
    state.timer = Timer(timer -> stept(inc, ctrls, state, showframe), dt, dt)
end

function sett(ctrls,state, showframe)
    tstr = get_value(ctrls.editt)
    try
        val = parse(Int, tstr)
        state.t = val
        updatet(ctrls, state)
        showframe(state)
    catch
        updatet(ctrls, state)
    end
end

function scalet(ctrls, state, showframe)
    state.t = round(Int, get_value(ctrls.scalet))
    updatet(ctrls, state)
    showframe(state)
end

export NavigationState,
    NavigationControls,
    init_navigation!

end


text:
using SAMTools
using Test

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


text:
@testset "grad: simple" begin
    args3 = (rand(4, 3), rand(4), rand(3))
    args4 = (rand(4, 3), rand(4), rand(3), rand(4))
    @test gradcheck((W, b, x) -> sum(W * x .+ b), args3...)
    @test gradcheck((W, b, x) -> sum(tanh.(W * x .+ b)), args3...)
    @test gradcheck((W, b, x, y) -> sum(abs2.(tanh.(W * x .+ b) .- y)), args4...)
end