text:
using MemoizedMethods, Test, Pkg

@test_throws LoadError eval(:(@memoize))
@test_throws LoadError eval(:(@memoize () = ()))

# you can't use test_throws in macros
arun = 0
@memoize function memadd(x::Int, y::Int)::Int
    global arun += 1
    print("evaluating memadd $x $y\n")
    return x + y
end
@test memadd(1, 2) == 3
@test memadd(1, 4) == 5
@test memadd(1, 2) == 3
@test arun == 2

run = 0
@memoize function simple(a)
    global run += 1
    a
end
@test simple(5) == 5
@test run == 1
@test simple(5) == 5
@test run == 1
@test simple(6) == 6
@test run == 2
@test simple(6) == 6
@test run == 2

map(forget!, methods(simple))
@test simple(6) == 6
@test run == 3
@test simple(6) == 6
@test run == 3

run = 0
lambda = @memoize (a) -> begin
    global run += 1
    a
end
@test lambda(5) == 5
@test run == 1
@test lambda(5) == 5
@test run == 1
@test lambda(6) == 6
@test run == 2
@test lambda(6) == 6
@test run == 2

forget!(lambda, Tuple{Any})
@test lambda(6) == 6
@test run == 3
@test lambda(6) == 6
@test run == 3

run = 0
@memoize function typed(a::Int)
    global run += 1
    a
end
@test typed(5) == 5
@test run == 1
@test typed(5) == 5
@test run == 1
@test typed(6) == 6
@test run == 2
@test typed(6) == 6
@test run == 2

run = 0
@memoize function default(a = 2)
    global run += 1
    a
end
@test default() == 2
@test run == 1
@test default() == 2
@test run == 1
@test default(2) == 2
@test run == 1
@test default(6) == 6
@test run == 2
@test default(6) == 6
@test run == 2

run = 0
@memoize function default_typed(a::Int = 2)
    global run += 1
    a
end
@test default_typed() == 2
@test run == 1
@test default_typed() == 2
@test run == 1
@test default_typed(2) == 2
@test run == 1
@test default_typed(6) == 6
@test run == 2
@test default_typed(6) == 6
@test run == 2

run = 0
@memoize function kw(; a = 2)
    global run += 1
    a
end
@test kw() == 2
@test run == 1
@test kw() == 2
@test run == 1
@test kw(a = 2) == 2
@test run == 1
@test kw(a = 6) == 6
@test run == 2
@test kw(a = 6) == 6
@test run == 2

run = 0
@memoize function kw_typed(; a::Int = 2)
    global run += 1
    a
end
@test kw_typed() == 2
@test run == 1
@test kw_typed() == 2
@test run == 1
@test kw_typed(a = 2) == 2
@test run == 1
@test kw_typed(a = 6) == 6
@test run == 2
@test kw_typed(a = 6) == 6
@test run == 2

run = 0
@memoize function default_kw(a = 1; b = 2)
    global run += 1
    (a, b)
end
@test default_kw() == (1, 2)
@test run == 1
@test default_kw() == (1, 2)
@test run == 1
@test default_kw(1, b = 2) == (1, 2)
@test run == 1
@test default_kw(1, b = 3) == (1, 3)
@test run == 2
@test default_kw(1, b = 3) == (1, 3)
@test run == 2
@test default_kw(2, b = 3) == (2, 3)
@test run == 3
@test default_kw(2, b = 3) == (2, 3)
@test run == 3

run = 0
@memoize function default_kw_typed(a::Int = 1; b::Int = 2)
    global run += 1
    (a, b)
end
@test default_kw_typed() == (1, 2)
@test run == 1
@test default_kw_typed() == (1, 2)
@test run == 1
@test default_kw_typed(1, b = 2) == (1, 2)
@test run == 1
@test default_kw_typed(1, b = 3) == (1, 3)
@test run == 2
@test default_kw_typed(1, b = 3) == (1, 3)
@test run == 2
@test default_kw_typed(2, b = 3) == (2, 3)
@test run == 3
@test default_kw_typed(2, b = 3) == (2, 3)
@test run == 3

run = 0
@memoize function required_kw(; a)
    global run += 1
    a
end
@test required_kw(a = 1) == 1
@test run == 1
@test required_kw(a = 1) == 1
@test run == 1
@test required_kw(a = 2) == 2
@test run == 2
@test required_kw(a = 2) == 2
@test run == 2
@test_throws UndefKeywordError required_kw()

run = 0
@memoize function ellipsis(a, b...)
    global run += 1
    (a, b)
end
@test ellipsis(1) == (1, ())
@test run == 1
@test ellipsis(1) == (1, ())
@test run == 1
@test ellipsis(1, 2, 3) == (1, (2, 3))
@test run == 2
@test ellipsis(1, 2, 3) == (1, (2, 3))
@test run == 2

run = 0
@memoize Dict() function kw_ellipsis(;a...)
    global run += 1
    a
end
@test isempty(kw_ellipsis())
@test run == 1
@test isempty(kw_ellipsis())
@test run == 1
@test kw_ellipsis(a = 1) == pairs((a = 1,))
@test run == 2
@test kw_ellipsis(a = 1) == pairs((a = 1,))
@test run == 2
@test kw_ellipsis(a = 1, b = 2) == pairs((a = 1, b = 2))
@test run == 3
@test kw_ellipsis(a = 1, b = 2) == pairs((a = 1, b = 2))
@test run == 3

run = 0
@memoize function multiple_dispatch(a::Int)
    global run += 1
    1
end
@memoize function multiple_dispatch(a::Float64)
    global run += 1
    2
end
@test multiple_dispatch(1) == 1
@test run == 1
@test multiple_dispatch(1) == 1
@test run == 1
@test multiple_dispatch(1.0) == 2
@test run == 2
@test multiple_dispatch(1.0) == 2
@test run == 2

run = 0
@memoize function where_clause(a::T) where T
    global run += 1
    T
end
@test where_clause(1) == Int
@test run == 1
@test where_clause(1) == Int
@test run == 1

function outer()
    run = 0
    @memoize function inner(x)
        run += 1
        x
    end
    @memoize function inner(x, y)
        run += 1
        x + y
    end
    @test inner(5) == 5
    @test run == 1
    @test inner(5) == 5
    @test run == 1
    @test inner(6) == 6
    @test run == 2
    @test inner(6) == 6
    @test run == 2
    @test inner(5, 1) == 6
    @test run == 3
    @test inner(5, 1) == 6
    @test run == 3
end
outer()
@test !@isdefined inner

function outer_multi(y)
    run = 0
    @memoize function inner(x)
        run += 1
        (x + 1, y, run)
    end
    #note that calling inner here would result in an error,
    #since both definitions of inner are evaluated before the
    #body of outer_overwrite runs, and the cache for the second definition
    #of inner has not been set up yet.
    @memoize function inner(x)
        run += 1
        (x, y, run)
    end
    @test inner(5) == (5, y, 1)
    @test run == 1
    @test inner(5) == (5, y, 1)
    @test run == 1
    @test inner(6) == (6, y, 2)
    @test run == 2
    @memoize function inner(x::String)
        run += 1
        (x, y, run)
    end
    return inner
end

inner_1 = outer_multi(7)
inner_2 = outer_multi(42)
@test inner_1(5) == (5, 7, 1)
@test inner_1(6) == (6, 7, 2)
@test inner_1(7) == (7, 7, 3)
@test inner_1("hello") == ("hello", 7, 4)
@test inner_2(7) == (7, 42, 3)
@test inner_2(5) == (5, 42, 1)
@test inner_2(6) == (6, 42, 2)
@test inner_2("goodbye") == ("goodbye", 42, 4)
@test inner_2("hello") == ("hello", 42, 5)
forget!(inner_1, Tuple{Any})
@test inner_1(5) == (5, 7, 1)
@test inner_2(6) == (6, 42, 2)
@test inner_1("hello") == ("hello", 7, 4)

trait_function(a, ::Bool) = (-a,)
run = 0
@memoize function trait_function(a, ::Int)
    global run += 1
    (a,)
end
@test trait_function(1, true) == (-1,)
@test run == 0
@test trait_function(2, true) == (-2,)
@test run == 0
@test trait_function(1, 1) == (1,)
@test run == 1
@test trait_function(1, 2) == (1,)
@test run == 1
@test trait_function(2, 2) == (2,)
@test run == 2
@test trait_function(2, 2) == (2,)
@test run == 2

run = 0
@memoize function trait_params(a, ::T) where {T}
    global run += 1
    (a, T)
end
@test trait_params(1, true) == (1, Bool)
@test run == 1
@test trait_params(1, false) == (1, Bool)
@test run == 1
@test trait_params(2, true) == (2, Bool)
@test run == 2
@test trait_params(2, false) == (2, Bool)
@test run == 2
@test trait_params(1, 3) == (1, Int)
@test run == 3
@test trait_params(1, 4) == (1, Int)
@test run == 3

run = 0
struct callable_object
    a
end
@memoize function (o::callable_object)(b)
    global run += 1
    (o.a, b)
end
@test callable_object(1)(2) == (1, 2)
@test run == 1
@test callable_object(1)(2) == (1, 2)
@test run == 1
@test callable_object(1)(3) == (1, 3)
@test run == 2
@test callable_object(1)(3) == (1, 3)
@test run == 2
@test callable_object(2)(3) == (2, 3)
@test run == 3
@test callable_object(2)(3) == (2, 3)
@test run == 3

run = 0
struct callable_trait_object{T}
    a::T
end
@memoize function (::callable_trait_object{T})(b) where {T}
    global run += 1
    (T, b)
end
@test callable_trait_object(1)(2) == (Int, 2)
@test run == 1
@test callable_trait_object(2)(2) == (Int, 2)
@test run == 1
@test callable_trait_object(false)(2) == (Bool, 2)
@test run == 2
@test callable_trait_object(true)(3) == (Bool, 3)
@test run == 3
@test callable_trait_object(1)(3) == (Int, 3)
@test run == 4
@test callable_trait_object(2)(3) == (Int, 3)
@test run == 4

run = 0
struct callable_type{T}
    a::T
end
@memoize function callable_type{T}(b) where {T}
    global run += 1
    (T, b)
end
@test callable_type{Int}(2) == (Int, 2)
@test run == 1
@test callable_type{Int}(2) == (Int, 2)
@test run == 1
@test callable_type{Int}(3) == (Int, 3)
@test run == 2
@test callable_type{Int}(3) == (Int, 3)
@test run == 2
@test callable_type{Bool}(3) == (Bool, 3)
@test run == 3
@test callable_type{Bool}(3) == (Bool, 3)
@test run == 3

run = 0
struct Constructable
    x
    @memoize function Constructable(a)
        global run += 1
        new(a)
    end
end
@test Constructable(5).x == 5
@test run == 1
@test Constructable(5).x == 5
@test run == 1
@test Constructable(6).x == 6
@test run == 2
@test Constructable(6).x == 6
@test run == 2
@test memories(Constructable, Tuple{Any}) isa IdDict

genrun = 0
@memoize function genspec(a)
    global genrun += 1
    a + 1
end
specrun = 0
@test genspec(5) == 6
@test genrun == 1
@test specrun == 0
@memoize function genspec(a::Int)
    global specrun += 1
    a + 2
end
@test genspec(5) == 7
@test genrun == 1
@test specrun == 1
@test genspec(5) == 7
@test genrun == 1
@test specrun == 1
@test genspec(true) == 2
@test genrun == 2
@test specrun == 1
@test invoke(genspec, Tuple{Any}, 5) == 6
@test genrun == 2
@test specrun == 1

map(forget!, methods(genspec, Tuple{Int}))
@test genspec(5) == 7
@test genrun == 2
@test specrun == 2

@memoize function typeinf(x)
    x + 1
end
@inferred typeinf(1)
@inferred typeinf(1.0)

finalized = false
@memoize function method_rewrite()
    x = []
    finalizer(x->(global finalized; finalized = true), x)
    x
end
method_rewrite()
@memoize function method_rewrite(x) end
GC.gc()
@test !finalized
@memoize function method_rewrite() end
GC.gc()
@test finalized

run = 0
""" documented function """
@memoize function documented_function(a)
    global run += 1
    a
end
@test strip(string(@doc documented_function)) == "documented function"
@test documented_function(1) == 1
@test run == 1
@test documented_function(1) == 1
@test run == 1
@test documented_function(2) == 2
@test run == 2
@test documented_function(2) == 2
@test run == 2

run = 0
@memoize function vararg_func(list::Vararg{Tuple{Int,Int}})
    global run += 1
    return list[1]
end
@test vararg_func((1,1), (1,1)) == (1,1)
@test run == 1
@test vararg_func((1,1), (1,1)) == (1,1)
@test run == 1
@test vararg_func((1,1), (1,2)) == (1,1)
@test run == 2
@test vararg_func((1,1), (1,2)) == (1,1)
@test run == 2

module MemoizeTest
using Test
using MemoizedMethods

const MyDict = Dict

run = 0
@memoize MyDict() function custom_dict(a)
    global run += 1
    a
end
@test custom_dict(1) == 1
@test run == 1
@test custom_dict(1) == 1
@test run == 1
@test custom_dict(2) == 2
@test run == 2
@test custom_dict(2) == 2
@test run == 2
@test custom_dict(1) == 1
@test run == 2

end # module

using .MemoizeTest
using .MemoizeTest: custom_dict

map(forget!, methods(custom_dict))
@test custom_dict(1) == 1
@test MemoizeTest.run == 3
@test custom_dict(1) == 1
@test MemoizeTest.run == 3

map(forget!, methods(MemoizeTest.custom_dict))
@test custom_dict(1) == 1
@test MemoizeTest.run == 4

if VERSION >= v"1.4.2"
    Pkg.activate(temp=true)
    Pkg.develop(path=joinpath(@__DIR__, "TestPrecompile"))
    using TestPrecompile

    @test TestPrecompile.run == 1
    @test TestPrecompile.forgetful(1)
    @test TestPrecompile.run == 1

    Pkg.develop(path=joinpath(@__DIR__, "TestPrecompile2"))
    using TestPrecompile2

    @test TestPrecompile.run == 1
    @test TestPrecompile.forgetful(2)
    @test TestPrecompile.run == 2
end

run = 0
@memoize Dict{Tuple{String},Int}() function dict_call(a::String)::Int
    global run += 1
    length(a)
end
@test dict_call("a") == 1
@test run == 1
@test dict_call("a") == 1
@test run == 1
@test dict_call("bb") == 2
@test run == 2
@test dict_call("bb") == 2

run = 0
@memoize Dict{__Key__,__Value__}() function auto_dict_call(a::String)::Int
    global run += 1
    length(a)
end
@test auto_dict_call("a") == 1
@test run == 1
@test auto_dict_call("a") == 1
@test run == 1
@test auto_dict_call("bb") == 2
@test run == 2
@test auto_dict_call("bb") == 2
@test run == 2
@test memories(auto_dict_call, Tuple{String}) isa Dict{Tuple{String}, Int}
forget!(auto_dict_call, Tuple{String})
@test length(memories(auto_dict_call, Tuple{String})) == 0

@memoize non_allocating(x) = x+1
non_allocating(10)
@test @allocated(non_allocating(10)) == 0


text:
using Flux
using Flux: gradient
using Flux.Optimise: update!
using DelimitedFiles, Statistics
using Parameters: @with_kw

# This replicates the housing data example from the Knet.jl readme. Although we
# could have reused more of Flux (see the mnist example), the library's
# abstractions are very lightweight and don't force you into any particular
# strategy.

# Struct to define hyperparameters
@with_kw mutable struct Hyperparams
    lr::Float64 = 0.1		# learning rate
    split_ratio::Float64 = 0.1	# Train Test split ratio, define percentage of data to be used as Test data
end

function get_processed_data(args)
    isfile("housing.data") ||
        download("https://raw.githubusercontent.com/MikeInnes/notebooks/master/housing.data",
            "housing.data")

    rawdata = readdlm("housing.data")'

    # The last feature is our target -- the price of the house.
    split_ratio = args.split_ratio # For the train test split

    x = rawdata[1:13,:]
    y = rawdata[14:14,:]

    # Normalise the data
    x = (x .- mean(x, dims = 2)) ./ std(x, dims = 2)

    # Split into train and test sets
    split_index = floor(Int,size(x,2)*split_ratio)
    x_train = x[:,1:split_index]
    y_train = y[:,1:split_index]
    x_test = x[:,split_index+1:size(x,2)]
    y_test = y[:,split_index+1:size(x,2)]

    train_data = (x_train, y_train)
    test_data = (x_test, y_test)

    return train_data,test_data
end

# Struct to define model
mutable struct model
    W::AbstractArray
    b::AbstractVector
end

# Function to predict output from given parameters
predict(x, m) = m.W*x .+ m.b

# Define the mean squared error function to be used in the loss
# function. An implementation is also available in the Flux package
# (https://fluxml.ai/Flux.jl/stable/models/losses/#Flux.Losses.mse).
meansquarederror(ŷ, y) = sum((ŷ .- y).^2)/size(y, 2)

function train(; kws...)
    # Initialize the Hyperparamters
    args = Hyperparams(; kws...)

    # Load the data
    (x_train,y_train),(x_test,y_test) = get_processed_data(args)

    # The model
    m = model((randn(1,13)),[0.])

    loss(x, y) = meansquarederror(predict(x, m), y)

    ## Training
    η = args.lr
    θ = params(m.W, m.b)

    for i = 1:500
        g = gradient(() -> loss(x_train, y_train), θ)
        for x in θ
            update!(x, g[x]*η)
        end
        if i%100==0
            @show loss(x_train, y_train)
        end
    end

    # Predict the RMSE on the test set
    err = meansquarederror(predict(x_test, m),y_test)
    println(err)
end
##===============mmy
using DifferentialEquations, Plots
## Setup ODE
function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end
u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob = ODEProblem(lotka_volterra,u0,tspan,p)
# Verify ODE solution
sol = solve(prob,Tsit5())
plot(sol)

# Generate data from the ODE
sol = solve(prob,Tsit5(),saveat=0.1)
A = sol[1,:] # length 101 vector
t = 0:0.1:10.0
scatter!(t,A)

using Flux, DiffEqFlux
# Build a neural network that sets the cost as the difference
# from the generated data and 1
p = [2.2, 1.0, 2.0, 0.4]
p = Flux.params(p) # Initial Parameter Vector

function predict_rd() # Our 1-layer neural network
  diffeq_rd(p,prob,Tsit5(),saveat=0.1)[1,:]
end
loss_rd() = sum(abs2,x-1 for x in predict_rd()) # loss function

# Optimize the parameters so the ODE's solution stays near 1
data = Iterators.repeated((), 100)
opt = ADAM(0.1)
Flux.train!(loss_rd, [p], data, opt, cb = cb)

##=====Poisson
using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux
@parameters x y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
# 2D PDE
eq  = Dxx(u(x,y)) + Dyy(u(x,y)) ~ -sin(pi*x)*sin(pi*y)
# Boundary conditions
bcs = [u(0,y) ~ 0.f0, u(1,y) ~ -sin(pi*1)*sin(pi*y),
       u(x,0) ~ 0.f0, u(x,1) ~ -sin(pi*x)*sin(pi*1)]
# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0), y ∈ IntervalDomain(0.0,1.0)]
# Discretization
dx = 0.1
# Neural network
dim = 2 # number of dimensions
chain = FastChain(FastDense(dim,16,Flux.σ),FastDense(16,16,Flux.σ),
                  FastDense(16,1))
discret = PhysicsInformedNN(chain, GridTraining(dx)) # PINNs model
pde_system = PDESystem(eq,bcs,domains,[x,y],[u])
prob = discretize(pde_system,discret)
res = GalacticOptim.solve(prob, Optim.BFGS(); maxiters=1000)
phi = discret.phi

cd(@__DIR__)
train()

ext: jl


text:
module ModelParameters
# Use the README as the module docs
@doc let
    path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    read(path, String)
end ModelParameters

import AbstractNumbers,
       ConstructionBase,
       ConstructionBaseExtras,
       Flatten,
       PrettyTables,
       Tables

using Setfield

export AbstractModel, Model, StaticModel

export AbstractParam, Param

export params, printparams, stripparams, update, update!, withunits, stripunits

include("interface.jl")
include("param.jl")
include("model.jl")
include("tables.jl")

end


text:
using Flux.Data
using Base.Test

@test cmudict()["CATASTROPHE"] == :[K,AH0,T,AE1,S,T,R,AH0,F,IY0].args

@test length(CMUDict.phones()) == 39

@test length(CMUDict.symbols()) == 84


text:
function fixargs(args::Vector{Any}, stub::Any)
    n = length(args)
    data = Array(Any, n)
    na = Array(Any, n)
    for i in 1:n
        if args[i] == :NA
            data[i] = stub
            na[i] = true
        else
            data[i] = args[i]
            na[i] = false
        end
    end
    return data, na
end

# We assume that data has at least one "value" that isn't NA
function findstub_vector(ex::Expr)
    if ex.args[1] != :NA
        return ex.args[1]
    end
    for i in 2:length(ex.args)
        if ex.args[i] != :NA
            return ex.args[i]
        end
    end
    error("Type of literal data cannot be inferred")
end

# We assume that data has at least one "value" that isn't NA
function findstub_matrix(ex::Expr)
    if ex.args[1].args[1] != :NA
        return ex.args[1].args[1]
    end
    nrows = length(ex.args)
    for row in 1:nrows
        subex = ex.args[row]
        for i in 1:length(subex.args)
            if subex.args[i] != :NA
                return subex.args[i]
            end
        end
    end
    error("Type of literal data cannot be inferred")
end

function parsevector(ex::Expr)
    stub = findstub_vector(ex)
    data, na = fixargs(ex.args, stub)
    return Expr(:vcat, data...), Expr(:vcat, na...)
end

function parsematrix(ex::Expr)
    stub = findstub_matrix(ex)
    nrows = length(ex.args)
    datarows = Array(Expr, nrows)
    narows = Array(Expr, nrows)
    for row in 1:nrows
        data, na = fixargs(ex.args[row].args, stub)
        datarows[row] = Expr(:row, data...)
        narows[row] = Expr(:row, na...)
    end
    return Expr(:vcat, datarows...), Expr(:vcat, narows...)
end

function parsedata(ex::Expr)
    if length(ex.args) == 0
        return :([]), Expr(:call, :Array, :Bool, 0)
    end
    if isa(ex.args[1], Expr) && ex.args[1].head == :row
        return parsematrix(ex)
    else
        return parsevector(ex)
    end
end

macro data(ex)
    if ex.head != :vcat
        return Expr(:call, :DataArray, esc(ex))
    end
    dataexpr, naexpr = parsedata(ex)
    return Expr(:call, :DataArray, esc(dataexpr), esc(naexpr))
end

macro pdata(ex)
    if ex.head != :vcat
        return Expr(:call, :PooledDataArray, esc(ex))
    end
    dataexpr, naexpr = parsedata(ex)
    return Expr(:call, :PooledDataArray, esc(dataexpr), esc(naexpr))
end


text:
using DEBTrait
using CSV, DataFrames, Statistics
using DifferentialEquations, Plots

n_polymers              = 0
n_monomers              = 6
n_microbes              = 1
n_enzymes               = 1
n_minerals              = 0

df_isolates             = CSV.read("/Users/glmarschmann/.julia/dev/DEBTrait/data/microbes2traits.csv", DataFrame, missingstring = "N/A")
df_monomers             = CSV.read("/Users/glmarschmann/.julia/dev/DEBTrait/files/thermodynamic_props_class.csv",DataFrame, missingstring = "N/A")

isolate_id              = ["positive", "negative"]
monomer_id              = collect(1:6)
monomer_concentration   = [1e-7, 5e-7, 1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 1.0]
#
BGE                     = zeros(2, 6, 15)
r                       = zeros(2, 6, 15)
m_E                     = zeros(2, 6, 15)
x                       = zeros(2, 6, 15)
#
function condition(u,t,integrator) # Event when event_f(u,t) == 0
  u[1] < 1e-8
end
affect!(integrator) = terminate!(integrator)
cb = DiscreteCallback(condition,affect!)

plt = plot()

for i in 1:2
    id_isolate              = isolate_id[i]
    df_p                    = filter(row -> row[:Rhizosphere_response] == id_isolate, df_isolates)
    #
    L_DNA                   = median(df_p.Genome_size)*ones(n_microbes)
    rRNA                    = median(df_p.rRNA_genes)*ones(n_microbes)
    min_gentime             = median(df_p.Min_gen_time)*ones(n_microbes)
    #
    z_sugars_p              = median(convert(Array{Float64,1}, df_p.z_sugars))*ones(1,1)
    z_organics_p            = median(convert(Array{Float64,1}, df_p.z_organic_acids))*ones(1,1)
    z_aminos_p              = median(convert(Array{Float64,1}, df_p.z_amino_acids))*ones(1,1)
    z_fattys_p              = median(convert(Array{Float64,1}, df_p.z_fatty_acids))*ones(1,1)
    z_nucleos_p             = median(convert(Array{Float64,1}, df_p.z_nucleotides))*ones(1,1)
    z_auxins_p              = median(convert(Array{Float64,1}, df_p.z_auxins))*ones(1,1)
    genome_distr_p          = vcat(z_sugars_p, z_organics_p, z_aminos_p, z_fattys_p, z_nucleos_p, z_auxins_p)
    #
    for j in 1:6
        n_monomers              = 1
        id_monomer              = monomer_id[j]
        Isolate                 = BaseGenome(L_DNA, rRNA, min_gentime, genome_distr_p[id_monomer]*ones(n_monomers,n_microbes))
        Monomer                 = BaseMonomer(df_monomers.D_S[id_monomer]*ones(n_monomers), df_monomers.N_C[id_monomer]*ones(n_monomers), df_monomers.N_N[id_monomer]*ones(n_monomers), df_monomers.y_DE[id_monomer]*ones(n_monomers))
        Soil                    = BaseSoil(0.0, 0.0, 1.0, 1*ones(n_microbes))
        #
        MetabolismC             = PMetabolismC{BaseGenome, Array{Float64,1}}(Isolate)
        MetabolismC.y_EM        = ones(n_microbes)
        MetabolismC.y_EX        = ones(n_microbes).*1.4
        MetabolismC.α           = zeros(n_microbes)

        if df_p.Rhizosphere_response[1] .== "positive"
            MetabolismC.α .= 0.01
        else
            MetabolismC.α .= 0.3
        end
        #
        AssimilationC           = PAssimilation(Isolate, Monomer, Soil, MetabolismC)
        #
        Setup                   = SetupBatchC(Isolate, n_polymers, n_monomers, n_microbes, n_enzymes, n_minerals)
        p                       = ParamsBatchC(Setup, Monomer, MetabolismC, AssimilationC)
        #########################
        p.n_polymers              = 0
        p.n_monomers              = 1
        p.n_microbes              = 1
        p.n_enzymes               = 1
        p.n_minerals              = 0
        #
        u0 = zeros(p.dim)
        for k in 1:15
            init_D0 = monomer_concentration[k]
            u0[1+p.n_polymers:p.n_polymers+p.n_monomers] = init_D0*ones(p.n_monomers)
            u0[1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes] = p.Setup.u0_reserves*1e-6
            u0[1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes] = p.Setup.u0_structures*1e-6
            #
            tspan = (0.0,  50.0)
            prob  = ODEProblem(rhs_BatchC!, u0, tspan, p)
            sol   = solve(prob, alg_hints=[:stiff], callback=cb)
            #
            BGE_tmp = [DEBTrait.calc_BGE(zeros(p.dim), sol.u[i], p, 0.0)[1] for i in 1:size(sol.t,1)]
            if isempty(BGE_tmp[BGE_tmp.>1e-8])
                BGE[i,j,k] = 0.0
            else
                BGE[i,j,k] = median(BGE_tmp[BGE_tmp.>1e-8])
            end
            r_tmp   = [DEBTrait.calc_growth_rate(sol.u[i], p)[1] for i in 1:size(sol.t,1)]
            r[i,j,k] = median(r_tmp)
            E = sol.u[end][1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes]
            V = sol.u[end][1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes]
            m_E_tmp = E./V
            m_E_tmp[m_E_tmp.>=1.0] .= 0.0
            m_E[i,j,k] = m_E_tmp[1]
            x_tmp = [DEBTrait.calc_growth_production(sol.u[i], p)[1][1] for i in 1:size(sol.t,1)]
            x[i,j,k] = median(x_tmp)
        end
    end
end

p = plot(legend=false)
for j in 1:6
    scatter!(monomer_concentration, m_E[1,j,:], xscale=:log10, color="green")
    scatter!(monomer_concentration, m_E[2,j,:], xscale=:log10, color="red")
    display(p)
end

xlabel!("Substrate concentration [mol/m3]")
ylabel!("reserve density [mol/mol]")

id_positive = 1
r_med_sugar = m_E[id_positive,1,:]
r_med_organics = m_E[id_positive,2,:]
r_med_aminos = m_E[id_positive,3,:]
r_med_fattys = m_E[id_positive,4,:]
r_med_nucleos = m_E[id_positive,5,:]
r_med_auxins = m_E[id_positive,6,:]
r_med_pos = (r_med_sugar+r_med_organics+r_med_aminos+r_med_fattys+r_med_nucleos+r_med_auxins)./6
scatter!(monomer_concentration, r_med_pos, markersize=8, color="green")

id_positive = 2
r_med_sugar = m_E[id_positive,1,:]
r_med_organics = m_E[id_positive,2,:]
r_med_aminos = m_E[id_positive,3,:]
r_med_fattys = m_E[id_positive,4,:]
r_med_nucleos = m_E[id_positive,5,:]
r_med_auxins = m_E[id_positive,6,:]
r_med_pos = (r_med_sugar+r_med_organics+r_med_aminos+r_med_fattys+r_med_nucleos+r_med_auxins)./6
scatter!(monomer_concentration, r_med_pos, markersize=8, color="red")


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

include("abstract_variate_transform.jl")
include("generic_transforms.jl")
include("distribution_transform.jl")
include("variate_transform_chain.jl")
include("named_tuple_transform.jl")


text:
using Pkg
 Pkg.activate(".")

verbose = true
using Revise
 if verbose   println("# Loading RvSpecML")    end
using RvSpectML
 if verbose   println("# Loading other packages")    end
 using DataFrames, Query, Statistics

# USER:  The default paths that specify where datafiles can be entered here or overridden in examples/data_paths.jl
target_subdir = "20190918"   # USER: Replace with directory of your choice
 fits_target_str = "Solar"
 output_dir = "examples/output"
 default_paths_to_search = [pwd(),"examples",joinpath(pkgdir(RvSpectML),"examples"), "/gpfs/group/ebf11/default/ebf11/neid_solar"]

make_plots = false
 write_ccf_to_csv = true
 write_rvs_to_csv = true
 write_template_to_csv = false
 write_spectral_grid_to_jld2 = true
 write_dcpca_to_csv = false

has_loaded_data = false
 has_computed_ccfs = false
 has_computed_rvs = false
 has_computed_tempalte = false
 has_computed_dcpca = false

if !has_loaded_data
   df_files = make_manifest(target_subdir, NEID )
   eval(code_to_include_param_jl())
   df_files_use

   if verbose println("# Reading in FITS files.")  end
  @time neid_data = map(NEID.read_solar_data,eachrow(df_files_use))
  has_loaded_data = true
end

if !has_computed_ccfs
 if verbose println("# Extracting order list timeseries from spectra.")  end
 order_list_timeseries = RvSpectML.make_order_list_timeseries(neid_data)
 order_list_timeseries = RvSpectML.filter_bad_chunks(order_list_timeseries,verbose=true)
 RvSpectML.normalize_spectra!(order_list_timeseries,neid_data);
end

begin
 if verbose println("# Reading line list for CCF: ", linelist_for_ccf_filename, ".")  end
 lambda_range_with_good_data = get_λ_range(neid_data)
 espresso_filename = joinpath(pkgdir(RvSpectML),"data","masks",linelist_for_ccf_filename)
 espresso_df = RvSpectML.read_linelist_espresso(espresso_filename)
 line_list_df = NEID.filter_line_list(espresso_df,first(neid_data).inst)
end
begin
    # Compute CCF's & measure RVs
 if verbose println("# Computing CCF.")  end
 mask_shape = RvSpectML.CCF.TopHatCCFMask(order_list_timeseries.inst, scale_factor=tophap_ccf_mask_scale_factor)
 line_list = RvSpectML.CCF.BasicLineList(line_list_df.lambda, line_list_df.weight)
 ccf_plan = RvSpectML.CCF.BasicCCFPlan(mask_shape = mask_shape, line_list=line_list, midpoint=ccf_mid_velocity)
 v_grid = RvSpectML.CCF.calc_ccf_v_grid(ccf_plan)
 @time ccfs = RvSpectML.CCF.calc_ccf_chunklist_timeseries(order_list_timeseries, ccf_plan)
 # Write CCFs to file
 if write_ccf_to_csv && false
    using CSV
    CSV.write(joinpath(output_dir,target_subdir * "_ccfs.csv"),Tables.table(ccfs',header=Symbol.(v_grid)))
 end
 has_computed_ccfs = true
end

if make_plots
   using Plots
   t_idx = 1
   using Plots
   plot(v_grid,ccfs[:,t_idx],label=:none)
   xlabel!("v (m/s)")
   ylabel!("CCF")
 end

if !has_computed_rvs
 if verbose println("# Measuring RVs from CCF.")  end
 #rvs_ccf_gauss = [ RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs[:,i],fit_type = "gaussian") for i in 1:length(order_list_timeseries) ]
 rvs_ccf_gauss = RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs,fit_type = "gaussian")
 # Store estimated RVs in metadata
 map(i->order_list_timeseries.metadata[i][:rv_est] = rvs_ccf_gauss[i]-mean(rvs_ccf_gauss), 1:length(order_list_timeseries) )
 if write_rvs_to_csv
    using CSV
    CSV.write(joinpath(output_dir,target_subdir * "_rvs_ccf.csv"),DataFrame("Time [MJD]"=>order_list_timeseries.times,"CCF RV [m/s]"=>rvs_ccf_gauss))
 end
 has_computed_rvs = true
end

#=
using DataFrames, Query, Statistics
vald_filename = joinpath(ancilary_solar_data_path,"VALD_Fe1_DP_rejectTelluricSlope0.0_badLineFilterESPRESSO-strict-NEID-BIS_overlapcutoff6e-05_depthcutoff0.05_allowBlends0_wavesReiners_depthssolar_nbin1depth0.mas")
 vald_df = RvSpectML.read_mask_vald(vald_filename)

λmin, λmax = get_λ_range(neid_data)
 orig_chunk_list_df = vald_df |>
    @filter(λmin <= _.lambda ) |>
    @filter( _.lambda < λmax ) |>
#    @filter( _.lambda < 6000.0 ) |>                       # Avoid tellurics at redder wavelengths
#    @filter( _.lambda >6157 || _.lambda < 6155  ) |>   # Avoid "line" w/ large variability
   DataFrame

find_overlapping_chunks(orig_chunk_list_df)
chunk_list_df = RvSpectML.merge_chunks(orig_chunk_list_df)
@assert find_overlapping_chunks(chunk_list_df) == nothing

# TODO:  FIX:  Problem with ESPRESSO mask having lots of lines getting meged into humongous chunks that don't fit into the order!
chunk_list_timeseries = RvSpectML.make_chunk_list_timeseries(neid_data,chunk_list_df)

# Check that no NaN's included
(chunk_list_timeseries, chunk_list_df) = RvSpectML.filter_bad_chunks(chunk_list_timeseries,chunk_list_df)
println(size(chunk_list_df), " vs ", num_chunks(chunk_list_timeseries) )

chunk_grids = map(c->RvSpectML.make_grid_for_chunk(chunk_list_timeseries,c), 1:num_chunks(chunk_list_timeseries) )
spectra_matrix = RvSpectML.pack_chunk_list_timeseries_to_matrix(chunk_list_timeseries,chunk_grids, alg=:Sinc, oversample_factor=1)
RvSpectML.make_template_spectra(chunk_list_timeseries)

=#

if !has_computed_tempalte
   if verbose println("# Making template spectra.")  end
   @time ( spectral_orders_matrix, f_mean, var_mean, deriv, deriv2, order_grids )  = RvSpectML.make_template_spectra(order_list_timeseries)

   if write_template_to_csv
      using CSV
      CSV.write(joinpath(output_dir,target_subdir * "_template.csv"),DataFrame("λ"=>spectral_orders_matrix.λ,"flux_template"=>f_mean,"var"=>var_mean, "dfluxdlnλ_template"=>deriv,"d²fluxdlnλ²_template"=>deriv2))
   end
   if write_spectral_grid_to_jld2
      using JLD2, FileIO
      save(joinpath(output_dir,target_subdir * "_matrix.jld2"), Dict("λ"=>spectral_orders_matrix.λ,"spectra"=>spectral_orders_matrix.flux,"var_spectra"=>spectral_orders_matrix.var,"flux_template"=>f_mean,"var"=>var_mean, "dfluxdlnλ_template"=>deriv,"d²fluxdlnλ²_template"=>deriv2))
   end
end

if !has_computed_dcpca
   if verbose println("# Performing Doppler constrained PCA analysis.")  end
   using MultivariateStats
   dcpca_out, M = RvSpectML.DCPCA.doppler_constrained_pca(spectral_orders_matrix.flux, deriv, rvs_ccf_gauss)
   frac_var_explained = 1.0.-cumsum(principalvars(M))./tvar(M)
   println("# Fraction of variance explained = ", frac_var_explained[1:min(5,length(frac_var_explained))])
   if write_dcpca_to_csv
      using CSV
      CSV.write(joinpath(output_dir,target_subdir * "_dcpca_basis.csv"),  Tables.table(M.proj) )
      CSV.write(joinpath(output_dir,target_subdir * "_dcpca_scores.csv"), Tables.table(dcpca_out) )
   end
end

#=
using MultivariateStats
function doppler_constrained_pca(flux::AbstractArray{T1,2}, deriv::AbstractVector{T2}, rvs::AbstractVector{T3} )  where { T1<:Real, T2<:Real, T3<:Real  }#, times::AbstractVector, rvs::AbstractVector)
  fm_perp = RvSpectML.compute_spectra_perp_doppler_shift(flux,deriv, rvs )
  idx_good_obs = 1:length(rvs)
  M = fit(PCA, fm_perp[:,idx_good_obs]; maxoutdim=10)
  pca_out = MultivariateStats.transform(M,fm_perp[:,idx_good_obs])
  return pca_out, M
end

=#
if make_plots
  using Plots
  # Set parameters for plotting analysis
  plt_order = 42
  plt_order_pix = 4500:5000
  plt = scatter(frac_var_explained, xlabel="Number of PCs", ylabel="Frac Variance Unexplained")
  display(plt)
end

if make_plots
  plt_order = 40
  RvSpectML.plot_basis_vectors(order_grids, f_mean, deriv, M.proj, idx_plt = spectral_orders_matrix.chunk_map[plt_order], num_basis=3 )
  #xlims!(5761.5,5766)
  xlims!(5315,5330)
end

if make_plots
  RvSpectML.plot_basis_scores(order_list_timeseries.times, rvs_ccf_gauss, dcpca_out, num_basis=3 )
end

if make_plots
  RvSpectML.plot_basis_scores_cor( rvs_ccf_gauss, dcpca_out, num_basis=3)
end

has_loaded_data = false
has_computed_ccfs = false
has_computed_rvs = false
has_computed_tempalte = false
has_computed_dcpca = false


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

function args_morespecific(a, b)
    sp = (ccall(:jl_type_morespecific, Cint, (Any,Any), a, b) != 0)
    if sp  # make sure morespecific(a,b) implies !morespecific(b,a)
        @test ccall(:jl_type_morespecific, Cint, (Any,Any), b, a) == 0
    end
    return sp
end

# issue #8652
let
    a  = Tuple{Type{T1}, T1} where T1<:Integer
    b2 = Tuple{Type{T2}, Integer} where T2<:Integer
    @test args_morespecific(a, b2)
    @test !args_morespecific(b2, a)
    a  = Tuple{Type{T1}, Ptr{T1}} where T1<:Integer
    b2 = Tuple{Type{T2}, Ptr{Integer}} where T2<:Integer
    @test args_morespecific(a, b2)
    @test !args_morespecific(b2, a)
end

# issue #11534
let
    t1 = Tuple{AbstractArray, Tuple{Vararg{Base.RangeIndex}}}
    t2 = Tuple{Array, T} where T<:Tuple{Vararg{Base.RangeIndex}}
    @test !args_morespecific(t1, t2)
    @test  args_morespecific(t2, t1)
end

let
    a = Tuple{Array{T,N}, Vararg{Int,N}} where T where N
    b = Tuple{Array,Int}
    @test  args_morespecific(a, b)
    @test !args_morespecific(b, a)
    a = Tuple{Array, Vararg{Int,N}} where N
    @test !args_morespecific(a, b)
    @test  args_morespecific(b, a)
end

# another specificity issue
_z_z_z_(x, y) = 1
_z_z_z_(::Int, ::Int, ::Vector) = 2
_z_z_z_(::Int, c...) = 3
@test _z_z_z_(1, 1, []) == 2

@test  args_morespecific(Tuple{T,Vararg{T}} where T<:Number,  Tuple{Number,Number,Vararg{Number}})
@test !args_morespecific(Tuple{Number,Number,Vararg{Number}}, Tuple{T,Vararg{T}} where T<:Number)

@test args_morespecific(Tuple{Array{T} where T<:Union{Float32,Float64,ComplexF32,ComplexF64}, Any},
                        Tuple{Array{T} where T<:Real, Any})

@test  args_morespecific(Tuple{1,T} where T, Tuple{Any})
@test  args_morespecific(Tuple{T} where T, Tuple{T,T} where T)
@test !args_morespecific(Type{T} where T<:Integer, Type{Any})

# issue #21016
@test args_morespecific(Tuple{IO, Core.TypeofBottom}, Tuple{IO, Type{T}} where T<:Number)

# issue #21382
@test args_morespecific(Tuple{Type{Pair{A,B} where B}} where A, Tuple{DataType})
@test args_morespecific(Tuple{Union{Int,String},Type{Pair{A,B} where B}} where A, Tuple{Integer,UnionAll})

# PR #21750
let A = Tuple{Any, Tuple{Vararg{Integer,N} where N}},
    B = Tuple{Any, Tuple{Any}},
    C = Tuple{Any, Tuple{}}
    @test args_morespecific(A, B)
    @test args_morespecific(C, A)
    @test args_morespecific(C, B)
end

# with bound varargs
_bound_vararg_specificity_1(::Type{Array{T,N}}, d::Vararg{Int, N}) where {T,N} = 0
_bound_vararg_specificity_1(::Type{Array{T,1}}, d::Int) where {T} = 1
@test _bound_vararg_specificity_1(Array{Int,1}, 1) == 1
@test _bound_vararg_specificity_1(Array{Int,2}, 1, 1) == 0

# issue #21710
@test args_morespecific(Tuple{Array}, Tuple{AbstractVector})
@test args_morespecific(Tuple{Matrix}, Tuple{AbstractVector})

# Method specificity
begin
    local f, A
    f(dims::Tuple{}, A::AbstractArray{T,0}) where {T} = 1
    f(dims::NTuple{N,Int}, A::AbstractArray{T,N}) where {T,N} = 2
    f(dims::NTuple{M,Int}, A::AbstractArray{T,N}) where {T,M,N} = 3
    A = zeros(2,2)
    @test f((1,2,3), A) == 3
    @test f((1,2), A) == 2
    @test f((), reshape([1])) == 1
    f(dims::NTuple{N,Int}, A::AbstractArray{T,N}) where {T,N} = 4
    @test f((1,2), A) == 4
    @test f((1,2,3), A) == 3
end

# a method specificity issue
c99991(::Type{T},x::T) where {T} = 0
c99991(::Type{UnitRange{T}},x::StepRangeLen{T}) where {T} = 1
c99991(::Type{UnitRange{T}},x::AbstractRange{T}) where {T} = 2
@test c99991(UnitRange{Float64}, 1.0:2.0) == 1
@test c99991(UnitRange{Int}, 1:2) == 2

# issue #17016, method specificity involving vararg tuples
T_17016{N} = Tuple{Any,Any,Vararg{Any,N}}
f17016(f, t::T_17016) = 0
f17016(f, t1::Tuple) = 1
@test f17016(0, (1,2,3)) == 0

@test !args_morespecific(Tuple{Type{Any}, Any}, Tuple{Type{T}, Any} where T<:VecElement)
@test !args_morespecific((Tuple{Type{T}, Any} where T<:VecElement), Tuple{Type{Any}, Any})

@test !args_morespecific(Tuple{Type{T}, Tuple{Any, Vararg{Any, N} where N}} where T<:Tuple{Any, Vararg{Any, N} where N},
                         Tuple{Type{Any}, Any})
@test !args_morespecific(Tuple{Type{T}, Tuple{Any, Vararg{Any, N} where N}} where T<:Tuple{Any, Vararg{Any, N} where N},
                         Tuple{Type{Tuple}, Tuple})
@test !args_morespecific(Tuple{Type{T}, T} where T<:Tuple{Any, Vararg{Any, N} where N},
                         Tuple{Type{T}, Any} where T<:VecElement)

@test args_morespecific(Tuple{Any, Tuple{}, Tuple{}}, Tuple{Any, Tuple{Any}})
@test args_morespecific(Tuple{Any, Tuple{Any}, Tuple{Any}}, Tuple{Any, Tuple{Any, Any}})
@test args_morespecific(Tuple{Any, Vararg{Tuple{}, N} where N}, Tuple{Any, Tuple{Any}})

@test  args_morespecific(Tuple{T, T} where T<:AbstractFloat, Tuple{T, T, T} where T<:AbstractFloat)
@test  args_morespecific(Tuple{T, Real, T} where T<:AbstractFloat, Tuple{T, T} where T<:Real)
@test  args_morespecific(Tuple{Real, Real}, Tuple{T, T, T} where T <: Real)
@test !args_morespecific(Tuple{Real, Real, Real}, Tuple{T, T, T} where T <: Real)
@test !args_morespecific(Tuple{Real, Real, Vararg{Real}}, Tuple{T, T, T} where T <: Real)
@test  args_morespecific(Tuple{Real, Real, Vararg{Int}}, Tuple{T, T, T} where T <: Real)

@test  args_morespecific(Tuple{Type{Base.Some{T}}} where T, Tuple{Type{T}, Any} where T)
@test !args_morespecific(Tuple{Type{Base.Some{T}}, T} where T, Tuple{Type{Base.Some{T}}} where T)

@test  args_morespecific(Tuple{Union{Base.StepRange{T, S} where S, Base.StepRangeLen{T, T, S} where S},
                               Union{Base.StepRange{T, S} where S, Base.StepRangeLen{T, T, S} where S}} where T,
                         Tuple{T, T} where T<:Union{Base.StepRangeLen, Base.LinRange})

@test args_morespecific(Tuple{Type{Tuple}, Any, Any},
                        Tuple{Type{Tuple{Vararg{E, N} where N}}, Any, Any} where E)

@test args_morespecific(Tuple{Type{Tuple{}}, Tuple{}},
                        Tuple{Type{T}, T} where T<:Tuple{Any, Vararg{Any, N} where N})

@test args_morespecific(Tuple{Type{CartesianIndex{N}}} where N,
                        Tuple{Type{CartesianIndex{N}},Vararg{Int,N}} where N)

# issue #22164
let A = Tuple{Type{D},D} where D<:Pair,
    B = Tuple{Type{Any}, Any},
    C = Tuple{Type{Pair}, Pair}
    @test  args_morespecific(C, A)
    @test !args_morespecific(A, B)
    @test !args_morespecific(C, B)
end

# issue #22338
let A = Tuple{Ref, Tuple{T}} where T,
    B = Tuple{Ref{T}, Tuple{Vararg{T}}} where T,
    C = Tuple{Ref{T}, Tuple{T}} where T
    @test  args_morespecific(C, A)
    @test  args_morespecific(C, B)
    @test !args_morespecific(A, B)
    @test !args_morespecific(B, A)
end

# issue #22339
let A = Tuple{T, Array{T, 1}} where T,
    B = Tuple{T} where T,
    C = Tuple{T} where T<:AbstractFloat
    @test args_morespecific(B, A)
    @test args_morespecific(C, B)
    @test args_morespecific(C, A)
end

# issue #22908
f22908(::Union) = 2
f22908(::Type{Union{Int, Float32}}) = 1
@test f22908(Union{Int, Float32}) == 1

let x = Type{Union{Tuple{T}, Tuple{Ptr{T}, Ptr{T}, Any}} where T},
    y = Type{Union{Tuple{T}, Tuple{Array{T, N} where N, Any, Array{T, N} where N, Any, Any}} where T}
    @test !args_morespecific(x, y)
    @test !args_morespecific(y, x)
    @test !args_morespecific(x.parameters[1], y.parameters[1])
    @test !args_morespecific(y.parameters[1], x.parameters[1])
end

let A = Tuple{Array{T,N}, Vararg{Int,N}} where {T,N},
    B = Tuple{Array, Int},
    C = Tuple{AbstractArray, Int, Array}
    @test args_morespecific(A, B)
    @test args_morespecific(B, C)
    @test args_morespecific(A, C)
end

# transitivity issue found in #26915
let A = Tuple{Vector, AbstractVector},
    B = Tuple{AbstractVecOrMat{T}, AbstractVector{T}} where T,
    C = Tuple{AbstractVecOrMat{T}, AbstractVecOrMat{T}} where T
    @test args_morespecific(A, B)
    @test args_morespecific(B, C)
    @test args_morespecific(A, C)
end

# issue #27361
f27361(::M) where M <: Tuple{2} = nothing
f27361(::M) where M <: Tuple{3} = nothing
@test length(methods(f27361)) == 2

# specificity of TypeofBottom
@test args_morespecific(Tuple{Core.TypeofBottom}, Tuple{DataType})
@test args_morespecific(Tuple{Core.TypeofBottom}, Tuple{Type{<:Tuple}})

max_stars_repo_name: helgee/julia

max_issues_repo_name: helgee/julia

max_forks_repo_name: helgee/julia


text:
# This file is a part of RadiationDetectorDSP.jl, licensed under the MIT License (MIT).

"""
    add_rect_pulse!(samples::AbstractSamples, start::Integer, pulselen::Integer, amplitude::Real = 1.0)

Add a rectangular pulse to `samples`.
"""
function add_rect_pulse!(samples::AbstractSamples, start::Integer, pulselen::Integer, amplitude::Real = 1.0)
    T = eltype(samples)
    A = view(samples, start:(start + pulselen - 1))
    A .= A .+ T(amplitude)
    samples
end
export add_rect_pulse!

"""
    gen_rect_pulse(tracelen::Integer, start::Integer, pulselen::Integer, amplitude::Real = 1.0)

Generate a rectangular pulse.
"""
function gen_rect_pulse(tracelen::Integer, start::Integer, pulselen::Integer, amplitude::Real = 1.0)
    add_rect_pulse!(fill(zero(amplitude), tracelen), start, pulselen, amplitude)
end
export gen_rect_pulse

# ToDo: add_gaussian_noise!(samples::AbstractSamples, sigma::AbstractFloat = 1.0) = ...


text:
using RecipesBase: apply_recipe
using Statistics

APERTURES = [
    CircularAperture(3, 3, 3),
    CircularAnnulus(3, 3, 2, 4),
    EllipticalAperture(3, 3, 4, 2, 45),
    EllipticalAnnulus(3, 3, 3, 4, 2, -26),
    RectangularAperture(3, 3, 3, 4, 15),
    RectangularAnnulus(3, 3, 3, 4, 2, -5)
]

@testset "Aperture Plots - $(typeof(ap))" for ap in APERTURES
    recipes = apply_recipe(Dict{Symbol,Any}(), ap)
    for rec in recipes
        @test getfield(rec, 1) == Dict{Symbol,Any}(
            :seriescolor => :red,
            :fillcolor => nothing,
            :linecolor => :match,
            :label => "",
            :seriestype => :shape,
            :aspect_ratio => :equal)

        # test to make sure our position is correct (should be +0.5 the given (x,y))
        points = rec.args[1]
        x_tot = y_tot = 0
        for point in points
            x_tot += point[1]
            y_tot += point[2]
        end
        N = length(points)
        x_mean = x_tot / N
        y_mean = y_tot / N
        @test x_mean ≈ 3.5 atol = x_mean / sqrt(N)
        @test y_mean ≈ 3.5 atol = y_mean / sqrt(N)
    end
end


text:
import Base.display

export Result, Responses, Plefka, Parameters, tspan, status_update

"""
    struct Result

    Struct with the result of the chemical species.

    Contains:
    time    Vector of time steps
    data    Matrix of mean copy numbers for each species

"""
struct Result
    time::Vector{Float64}
    data::Array{Float64,2}
end

function Base.display(res::Result)
	display([res.time'; res.data])
end

"""
    struct Responses

    Struct with the result of the response functions.

    Contains:
    time    Vector of time steps
    data    3-dim Array with species, time, time

"""
struct Responses
    time::Vector{Float64}
    data::Array{Float64,3}
end

function Base.display(resp::Responses)
	display(permutedims(resp.data, (2,3,1)))
end

"""
	struct Fields

	Struct with all auxiliary fields
	Differentiate between Plefka theory!

	Contains:
	hatTheta1,
	hatTheta2,
	hatR1,
	hatR2

"""
mutable struct Fields_lin1
	hatTheta1::Array{Float64,2}
end
function Base.display(fields::Fields_lin1)
	display(fields.hatTheta1)
end

mutable struct Fields_lin2
	hatTheta1::Array{Float64,2}
	hatTheta2::Array{Float64,2}
end
function Base.display(fields::Fields_lin2)
	display(fields.hatTheta1)
	display(fields.hatTheta2)
end

mutable struct Fields_quad1
	hatTheta1::Array{Float64,2}
	hatR1::Array{Float64,2}
end
function Base.display(fields::Fields_quad1)
	display(fields.hatTheta1)
	display(fields.hatR1)
end

mutable struct Fields_quad2
	hatTheta1::Array{Float64,2}
	hatTheta2::Array{Float64,2}
	hatR1::Array{Float64,2}
	hatR2::Array{Float64,3}
end
function Base.display(fields::Fields_quad2)
	display(fields.hatTheta1)
	display(fields.hatTheta2)
	display(fields.hatR1)
	display(permutedims(fields.hatR2, (2,3,1)))
end

"""
	struct Plefka

	Struct with all information about the Plefka approximation.

	Contains:
	alpha		number between 0 and 1, typically == 1
	alphaOrder	order of alpha expansion, either 1 or 2
	orderParam	either 'linear' or 'quad'

"""
struct Plefka
	α::Float64
	alphaOrder::Int64
	orderParam::String
	function Plefka(α::Float64, alphaOrder::Int64, orderParam::String)
		if α < 0.0 || α > 1.0
			error("PLEFKA: invalid alpha value! Between 0 and 1!")
		end
		if alphaOrder != 1 && alphaOrder != 2
			error("PLEFKA: invalid alpha order! Either 1 or 2!")
		end
		if orderParam != "linear" && orderParam != "quad"
			error("PLEFKA: invalid order parameter! Either 'linear' or 'quad'!")
		end
		new(α, alphaOrder, orderParam)
	end
end

"""
	mutable struct Parameters

	inherits all important information about the
	parameters of the system

	x0		initial condition
	k		reaction rates
	s_i		stoichiometric products
	r_i		stoichiometric reactants
	t_init	initial time
	t_final	final time
	delta_t	time discretization

"""
mutable struct Parameters
	x0::Vector{Float64}
	k::Vector{Vector{Float64}}
	s_i::Matrix{Int64}
	r_i::Matrix{Int64}
	t_init::Float64
	t_final::Float64
	delta_t::Float64
end

"""
	tspan(p::Parameters)

	return an array of the time grid
	according to the time discretization
	in the Parameters set p.

"""
function tspan(p::Parameters)
	return collect(p.t_init:p.delta_t:p.t_final)
end

"""
	status_update(t, len_ts)

	Prints the status in percent from [0..100%] in 10% steps

	t		loop index
	len_ts	length of the loop

	Put in a loop!

"""
function status_update(t, len_ts)

	if t == 1
		println("[  0%]")
	elseif t%(len_ts÷10) == 0 && t < 0.91*len_ts
		println("[ ", round(Int, 100*t/len_ts), "%]")
	elseif t == len_ts-1
		println("[100%]")
	end

end

max_stars_repo_name: mkurjahn/CRN

max_issues_repo_name: mkurjahn/CRN

max_forks_repo_name: mkurjahn/CRN


text:
@testset "Math Intrinsics" begin
    for intr in AMDGPUnative.MATH_INTRINSICS
        jlintr = intr.jlname
        if intr.isbroken || !(isdefined(Base, jlintr) || isdefined(SpecialFunctions, jlintr))
            @test_skip "$jlintr()"
            continue
        end
        modname = (isdefined(Base, jlintr) ? Base : SpecialFunctions)
        # FIXME: Handle all input and output args
        T = intr.inp_args[1]
        intr_kern = Symbol("intr_$(jlintr)_$T")
        @eval begin
            function $intr_kern(a)
                i = threadIdx().x
                a[i] = AMDGPUnative.$jlintr(a[i])
                return nothing
            end
            dims = (8,)
            a = rand($T, dims)
            if $(QuoteNode(jlintr)) == :acosh
                a .+= one($T)
            end
            d_a = HSAArray(a)
            len = prod(dims)

            @debug begin
                aspace = AMDGPUnative.AS.Global
                arrdims = ndims(a)
                arrT = ROCDeviceArray{Float32,arrdims,aspace}
                @debug "LLVM IR"
                AMDGPUnative.code_llvm($intr_kern, Tuple{arrT}; kernel=true)
                @debug "GCN Device Code"
                AMDGPUnative.code_gcn($intr_kern, Tuple{arrT}; kernel=true)
                ""
            end

            wait(@roc groupsize=len $intr_kern(d_a))
            _a = Array(d_a)
            @test $modname.$jlintr.(a) ≈ _a
        end
    end
end


text:
struct LSSSchur{wBType,wEType,BType,EType}
  wBinv::wBType
  wEinv::wEType
  B::BType
  E::EType
end

struct LSSSensitivityFunction{iip,F,A,J,JP,S,PJ,UF,PF,JC,PJC,Alg,fc,JM,pJM,MM,CV,
     PGPU,PGPP,CONFU,CONGP,DG} <: DiffEqBase.AbstractODEFunction{iip}
  f::F
  analytic::A
  jac::J
  jac_prototype::JP
  sparsity::S
  paramjac::PJ
  uf::UF
  pf::PF
  J::JM
  pJ::pJM
  jac_config::JC
  paramjac_config::PJC
  alg::Alg
  numparams::Int
  numindvar::Int
  f_cache::fc
  mass_matrix::MM
  colorvec::CV
  pgpu::PGPU
  pgpp::PGPP
  pgpu_config::CONFU
  pgpp_config::CONGP
  dg_val::DG
end

function LSSSensitivityFunction(sensealg,f,analytic,jac,jac_prototype,sparsity,paramjac,u0,
                                    alg,p,f_cache,mm,
                                    colorvec,tspan,g,dg)

  uf = DiffEqBase.UJacobianWrapper(f,tspan[1],p)
  pf = DiffEqBase.ParamJacobianWrapper(f,tspan[1],copy(u0))

  if DiffEqBase.has_jac(f)
    jac_config = nothing
  else
    jac_config = build_jac_config(sensealg,uf,u0)
  end

  if DiffEqBase.has_paramjac(f)
    paramjac_config = nothing
  else
    paramjac_config = build_param_jac_config(sensealg,pf,u0,p)
  end
  numparams = length(p)
  numindvar = length(u0)
  J = Matrix{eltype(u0)}(undef,numindvar,numindvar)
  pJ = Matrix{eltype(u0)}(undef,numindvar,numparams) # number of funcs size

  # compute gradients of objective
  if dg != nothing
    pgpu = nothing
    pgpp = nothing
    pgpu_config = nothing
    pgpp_config = nothing
    if dg isa Tuple && length(dg) == 2
      dg_val = (similar(u0, numindvar),similar(u0, numparams))
      dg_val[1] .= false
      dg_val[2] .= false
    else
      dg_val = similar(u0, numindvar) # number of funcs size
      dg_val .= false
    end
  else
    pgpu = UGradientWrapper(g,tspan[1],p) # ∂g∂u
    pgpp = ParamGradientWrapper(g,tspan[1],u0) #∂g∂p
    pgpu_config = build_grad_config(sensealg,pgpu,u0,tspan[1])
    pgpp_config = build_grad_config(sensealg,pgpp,p,tspan[1])
    dg_val = (similar(u0, numindvar),similar(u0, numparams))
    dg_val[1] .= false
    dg_val[2] .= false
  end

  LSSSensitivityFunction{isinplace(f),typeof(f),typeof(analytic),
                             typeof(jac),typeof(jac_prototype),typeof(sparsity),
                             typeof(paramjac),
                             typeof(uf),
                             typeof(pf),typeof(jac_config),
                             typeof(paramjac_config),typeof(alg),
                             typeof(f_cache),
                             typeof(J),typeof(pJ),typeof(mm),typeof(f.colorvec),
                             typeof(pgpu),typeof(pgpp),typeof(pgpu_config),typeof(pgpp_config),typeof(dg_val)}(
                             f,analytic,jac,jac_prototype,sparsity,paramjac,uf,pf,J,pJ,
                             jac_config,paramjac_config,alg,
                             numparams,numindvar,f_cache,mm,colorvec,
                             pgpu,pgpp,pgpu_config,pgpp_config,dg_val)
end

struct ForwardLSSProblem{A,C,solType,dtType,umidType,dudtType,SType,Ftype,bType,ηType,wType,vType,windowType,
    ΔtType,G0,G,DG,resType}
  sensealg::A
  diffcache::C
  sol::solType
  dt::dtType
  umid::umidType
  dudt::dudtType
  S::SType
  F::Ftype
  b::bType
  η::ηType
  w::wType
  v::vType
  window::windowType
  Δt::ΔtType
  Nt::Int
  g0::G0
  g::G
  dg::DG
  res::resType
end

function ForwardLSSProblem(sol, sensealg::ForwardLSS, g, dg = nothing;
                            kwargs...)

  @unpack f, p, u0, tspan = sol.prob
  isinplace = DiffEqBase.isinplace(f)

  p == nothing && error("You must have parameters to use parameter sensitivity calculations!")
  !(sol.u isa AbstractVector) && error("`u` has to be an AbstractVector.")

  sense = LSSSensitivityFunction(sensealg,f,f.analytic,f.jac,
                                     f.jac_prototype,f.sparsity,f.paramjac,
                                     u0,sensealg,
                                     p,similar(u0),f.mass_matrix,
                                     f.colorvec,
                                     tspan,g,dg)

  @unpack numparams, numindvar = sense
  Nt = length(sol.t)
  Ndt = Nt-one(Nt)

  # pre-allocate variables
  dt = similar(sol.t, Ndt)
  umid = Matrix{eltype(u0)}(undef,numindvar,Ndt)
  dudt = Matrix{eltype(u0)}(undef,numindvar,Ndt)
  # compute their values
  discretize_ref_trajectory!(dt, umid, dudt, sol, Ndt)

  S = LSSSchur(dt,u0,numindvar,Nt,Ndt,sensealg.alpha)

  if sensealg.alpha isa Number
    η = similar(dt,Ndt)
    window = nothing
    g0 = g(u0,p,tspan[1])
  else
    η = nothing
    window = similar(dt,Nt)
    g0 = nothing
  end

  b = Matrix{eltype(u0)}(undef,numindvar*Ndt,numparams)
  w = similar(dt,numindvar*Ndt)
  v = similar(dt,numindvar*Nt)

  Δt = tspan[2] - tspan[1]
  wB!(S,Δt,Nt,numindvar,dt)
  wE!(S,Δt,dt,sensealg.alpha)

  B!(S,dt,umid,sense,sensealg)
  E!(S,dudt,sensealg.alpha)
  F = SchurLU(S)

  res = similar(u0, numparams)

  ForwardLSSProblem{typeof(sensealg),typeof(sense),typeof(sol),typeof(dt),
    typeof(umid),typeof(dudt),
    typeof(S),typeof(F),typeof(b),typeof(η),typeof(w),typeof(v),typeof(window),typeof(Δt),
    typeof(g0),typeof(g),typeof(dg),typeof(res)}(sensealg,sense,sol,dt,umid,dudt,S,F,b,η,w,v,
    window,Δt,Nt,g0,g,dg,res)
end

function LSSSchur(dt,u0,numindvar,Nt,Ndt,alpha)
  wBinv = similar(dt,numindvar*Nt)
  if alpha isa Number
    wEinv = similar(dt,Ndt)
    E = Matrix{eltype(u0)}(undef,numindvar*Ndt,Ndt)
  else
    wEinv = nothing
    E = nothing
  end
  B = Matrix{eltype(u0)}(undef,numindvar*Ndt,numindvar*Nt)

  LSSSchur(wBinv,wEinv,B,E)
end

# compute discretized reference trajectory
function discretize_ref_trajectory!(dt, umid, dudt, sol, Ndt)
  for i=1:Ndt
    tr = sol.t[i+1]
    tl = sol.t[i]
    ur = sol.u[i+1]
    ul = sol.u[i]
    dt[i] = tr-tl
    copyto!((@view umid[:,i]), (ur + ul)/2)
    copyto!((@view dudt[:,i]), (ur - ul)/dt[i])
  end
  return nothing
end

function wB!(S::LSSSchur,Δt,Nt,numindvar,dt)
  @unpack wBinv = S
  fill!(wBinv, one(Δt))
  dim = numindvar * Nt
  tmp = @view wBinv[1:numindvar]
  tmp ./= dt[1]
  tmp = @view wBinv[dim-2:end]
  tmp ./= dt[end]
  for indx = 2:Nt-1
    tmp = @view wBinv[(indx-1)*numindvar+1:indx*numindvar]
    tmp ./= (dt[indx]+dt[indx-1])
  end

  wBinv .*= 2*Δt
  return nothing
end

wE!(S::LSSSchur,Δt,dt,alpha::Union{CosWindowing,Cos2Windowing}) = nothing

function wE!(S::LSSSchur,Δt,dt,alpha)
  @unpack wEinv = S
  @. wEinv = Δt/(alpha^2*dt)
  return nothing
end

function B!(S::LSSSchur,dt,umid,sense,sensealg)
  @unpack B = S
  @unpack f,J,uf,numindvar,f_cache,jac_config = sense

  fill!(B, zero(eltype(J)))

  for (i,u) in enumerate(eachcol(umid))
    if DiffEqBase.has_jac(f)
      f.jac(J,u,uf.p,uf.t) # Calculate the Jacobian into J
    else
      jacobian!(J, uf, u, f_cache, sensealg, jac_config)
    end
    B0 = @view B[(i-1)*numindvar+1:i*numindvar,i*numindvar+1:(i+1)*numindvar]
    B1 = @view B[(i-1)*numindvar+1:i*numindvar,(i-1)*numindvar+1:i*numindvar]
    B0 .+= I/dt[i] - J/2
    B1 .+= -I/dt[i] -J/2
  end
  return nothing
end

E!(S::LSSSchur,dudt,alpha::Union{CosWindowing,Cos2Windowing}) = nothing

function E!(S::LSSSchur,dudt,alpha)
  @unpack E = S
  numindvar, Ndt = size(dudt)
  for i=1:Ndt
    tmp = @view E[(i-1)*numindvar+1:i*numindvar,i]
    copyto!(tmp, (@view dudt[:,i]))
  end
  return nothing
end

# compute Schur
function SchurLU(S::LSSSchur)
  @unpack B, E, wBinv, wEinv = S
  Smat = B*Diagonal(wBinv)*B'
  (wEinv !== nothing) && (Smat .+= E*Diagonal(wEinv)*E')
  F = lu!(Smat)
  return F
end

function b!(b, prob::ForwardLSSProblem)
  @unpack diffcache, umid, sensealg = prob
  @unpack f, f_cache, pJ, pf, paramjac_config, uf, numindvar = diffcache

  for (i,u) in enumerate(eachcol(umid))
    if DiffEqBase.has_paramjac(f)
      f.paramjac(pJ, u, uf.p, pf.t)
    else
      pf.u = u
      jacobian!(pJ, pf, uf.p, f_cache, sensealg, paramjac_config)
    end
    tmp = @view b[(i-1)*numindvar+1:i*numindvar,:]
    copyto!(tmp, pJ)
  end
  return nothing
end

function shadow_forward(prob::ForwardLSSProblem; t0skip=zero(prob.Δt), t1skip=zero(prob.Δt))
  shadow_forward(prob,prob.sensealg,prob.sensealg.alpha,t0skip,t1skip)
end

function shadow_forward(prob::ForwardLSSProblem,sensealg::ForwardLSS,alpha::Number,t0skip,t1skip)
  @unpack sol, S, F, window, Δt, diffcache, b, w, v, η, res, g, g0, dg, umid = prob
  @unpack wBinv, wEinv, B, E = S
  @unpack dg_val, pgpu, pgpu_config, pgpp, pgpp_config, numparams, numindvar, uf = diffcache

  n0 = searchsortedfirst(sol.t, sol.t[1]+t0skip)
  n1 = searchsortedfirst(sol.t, sol.t[end]-t1skip)

  b!(b,prob)

  ures = @view sol.u[n0:n1]
  umidres = @view umid[:,n0:n1-1]

  # reset
  res .*=false

  for i=1:numparams
    #running average
    g0 *= false
    bpar = @view b[:,i]
    w .= F\bpar
    v .= Diagonal(wBinv)*(B'*w)
    η .= Diagonal(wEinv)*(E'*w)

    ηres = @view η[n0:n1-1]

    for (j,u) in enumerate(ures)
      vtmp = @view v[(n0+j-2)*numindvar+1:(n0+j-1)*numindvar]
      #  final gradient result for ith parameter
      if dg===nothing
        if dg_val isa Tuple
          DiffEqSensitivity.gradient!(dg_val[1], pgpu, u, sensealg,pgpu_config)
          DiffEqSensitivity.gradient!(dg_val[2], pgpp, uf.p, sensealg,pgpp_config)
          res[i] += dot(dg_val[1],vtmp)
          res[i] += dg_val[2][i]
        else
          DiffEqSensitivity.gradient!(dg_val, pgpu, u, sensealg,pgpu_config)
          res[i] += dot(dg_val,vtmp)
        end
      else
        if dg_val isa Tuple
          dg[1](dg_val[1],u,uf.p,nothing,n0+j-1)
          dg[2](dg_val[2],u,uf.p,nothing,n0+j-1)
          res[i] -= dot(dg_val[1],vtmp)
          res[i] -= dg_val[2][i]
        else
          dg(dg_val,u,uf.p,nothing,n0+j-1)
          res[i] -= dot(dg_val,vtmp)
        end
      end
    end
    # mean value
    res[i] = res[i]/(n1-n0+1)

    for (j,u) in enumerate(eachcol(umidres))
      # compute objective
      gtmp = g(u,uf.p,nothing)
      g0 += gtmp
      res[i] -= ηres[j]*gtmp/(n1-n0)
    end
    res[i] = res[i] + sum(ηres)*g0/(n1-n0)^2

  end
  return res
end

function shadow_forward(prob::ForwardLSSProblem,sensealg::ForwardLSS,alpha::CosWindowing,t0skip,t1skip)
  @unpack sol, S, F, window, Δt, diffcache, b, w, v, dg, res = prob
  @unpack wBinv, B = S
  @unpack dg_val, pgpu, pgpu_config, pgpp, pgpp_config, numparams, numindvar, uf = diffcache

  b!(b,prob)

  # windowing (cos)
  @. window = (sol.t-sol.t[1])*convert(eltype(Δt),2*pi/Δt)
  @. window = one(eltype(window)) - cos(window)
  window ./= sum(window)

  res .*= false

  for i=1:numparams
    bpar = @view b[:,i]
    w .= F\bpar
    v .= Diagonal(wBinv)*(B'*w)
    for (j,u) in enumerate(sol.u)
      vtmp = @view v[(j-1)*numindvar+1:j*numindvar]
      #  final gradient result for ith parameter
      if dg===nothing
        if dg_val isa Tuple
          DiffEqSensitivity.gradient!(dg_val[1], pgpu, u, sensealg,pgpu_config)
          DiffEqSensitivity.gradient!(dg_val[2], pgpp, uf.p, sensealg,pgpp_config)
          res[i] += dot(dg_val[1],vtmp)*window[j]
          res[i] += dg_val[2][i]*window[j]
        else
          DiffEqSensitivity.gradient!(dg_val, pgpu, u, sensealg,pgpu_config)
          res[i] += dot(dg_val,vtmp)*window[j]
        end
      else
        if dg_val isa Tuple
          dg[1](dg_val[1],u,uf.p,nothing,j)
          dg[2](dg_val[2],u,uf.p,nothing,j)
          res[i] -= dot(dg_val[1],vtmp)*window[j]
          res[i] -= dg_val[2][i]*window[j]
        else
          dg(dg_val,u,uf.p,nothing,j)
          res[i] -= dot(dg_val,vtmp)*window[j]
        end
      end
    end
  end
  return res
end

function shadow_forward(prob::ForwardLSSProblem,sensealg::ForwardLSS,alpha::Cos2Windowing,t0skip,t1skip)
    @unpack sol, S, F, window, Δt, diffcache, b, w, v, dg, res = prob
    @unpack wBinv, B = S
    @unpack dg_val, pgpu, pgpu_config, pgpp, pgpp_config, numparams, numindvar, uf = diffcache

    b!(b,prob)

    res .*= false

    # windowing cos2
    @. window = (sol.t-sol.t[1])*convert(eltype(Δt),2*pi/Δt)
    @. window = (one(eltype(window)) - cos(window))^2
    window ./= sum(window)

    for i=1:numparams
      bpar = @view b[:,i]
      w .= F\bpar
      v .= Diagonal(wBinv)*(B'*w)
      for (j,u) in enumerate(sol.u)
        vtmp = @view v[(j-1)*numindvar+1:j*numindvar]
        #  final gradient result for ith parameter
        if dg===nothing
          if dg_val isa Tuple
            DiffEqSensitivity.gradient!(dg_val[1], pgpu, u, sensealg,pgpu_config)
            DiffEqSensitivity.gradient!(dg_val[2], pgpp, uf.p, sensealg,pgpp_config)
            res[i] += dot(dg_val[1],vtmp)*window[j]
            res[i] += dg_val[2][i]*window[j]
          else
            DiffEqSensitivity.gradient!(dg_val, pgpu, u, sensealg,pgpu_config)
            res[i] += dot(dg_val,vtmp)*window[j]
          end
        else
          if dg_val isa Tuple
            dg[1](dg_val[1],u,uf.p,nothing,j)
            dg[2](dg_val[2],u,uf.p,nothing,j)
            res[i] -= dot(dg_val[1],vtmp)*window[j]
            res[i] -= dg_val[2][i]*window[j]
          else
            dg(dg_val,u,uf.p,nothing,j)
            res[i] -= dot(dg_val,vtmp)*window[j]
          end
        end
      end
    end
    return res
end

struct AdjointLSSProblem{A,C,solType,dtType,umidType,dudtType,SType,FType,hType,bType,wType,
    ΔtType,G0,G,DG,resType}
  sensealg::A
  diffcache::C
  sol::solType
  dt::dtType
  umid::umidType
  dudt::dudtType
  S::SType
  F::FType
  h::hType
  b::bType
  wa::wType
  Δt::ΔtType
  Nt::Int
  g0::G0
  g::G
  dg::DG
  res::resType
end

function AdjointLSSProblem(sol, sensealg::AdjointLSS, g, dg = nothing;
                            kwargs...)

  @unpack f, p, u0, tspan = sol.prob
  isinplace = DiffEqBase.isinplace(f)

  p == nothing && error("You must have parameters to use parameter sensitivity calculations!")
  !(sol.u isa AbstractVector) && error("`u` has to be an AbstractVector.")

  sense = LSSSensitivityFunction(sensealg,f,f.analytic,f.jac,
                                     f.jac_prototype,f.sparsity,f.paramjac,
                                     u0,sensealg,
                                     p,similar(u0),f.mass_matrix,
                                     f.colorvec,
                                     tspan,g,dg)

  @unpack numparams, numindvar = sense
  Nt = length(sol.t)
  Ndt = Nt-one(Nt)

  # pre-allocate variables
  dt = similar(sol.t, Ndt)
  umid = Matrix{eltype(u0)}(undef,numindvar,Ndt)
  dudt = Matrix{eltype(u0)}(undef,numindvar,Ndt)
  # compute their values
  discretize_ref_trajectory!(dt, umid, dudt, sol, Ndt)

  S = LSSSchur(dt,u0,numindvar,Nt,Ndt,sensealg.alpha)

  if sensealg.alpha isa Number
    g0 = g(u0,p,tspan[1])
  else
    g0 = nothing
  end

  b = Vector{eltype(u0)}(undef,numindvar*Ndt)
  h = Vector{eltype(u0)}(undef,Ndt)
  wa = similar(dt,numindvar*Ndt)

  Δt = tspan[2] - tspan[1]
  wB!(S,Δt,Nt,numindvar,dt)
  wE!(S,Δt,dt,sensealg.alpha)

  B!(S,dt,umid,sense,sensealg)
  E!(S,dudt,sensealg.alpha)
  F = SchurLU(S)
  wBcorrect!(S,sol,g,Nt,sense,sensealg,dg)

  h!(h,g0,g,umid,p,S.wEinv)

  res = similar(u0, numparams)

  AdjointLSSProblem{typeof(sensealg),typeof(sense),typeof(sol),typeof(dt),
    typeof(umid),typeof(dudt),
    typeof(S),typeof(F),typeof(h),typeof(b),typeof(wa),typeof(Δt),
    typeof(g0),typeof(g),typeof(dg),typeof(res)}(sensealg,sense,sol,dt,umid,dudt,S,F,h,b,wa,
    Δt,Nt,g0,g,dg,res)
end

function h!(h,g0,g,u,p,wEinv)

  for (j,uj) in enumerate(eachcol(u))
    # compute objective
    h[j] = g(uj,p,nothing)
  end
  h .= -(h .- mean(h)) / (size(u)[2])

  @. h = wEinv*h

  return nothing
end

function wBcorrect!(S,sol,g,Nt,sense,sensealg,dg)
  @unpack dg_val, pgpu, pgpu_config, numparams, numindvar, uf = sense
  @unpack wBinv = S

  for (i,u) in enumerate(sol.u)
    _wBinv = @view wBinv[(i-1)*numindvar+1:i*numindvar]
    if dg === nothing
      if dg_val isa Tuple
        DiffEqSensitivity.gradient!(dg_val[1], pgpu, u, sensealg, pgpu_config)
        @. _wBinv = _wBinv*dg_val[1]/Nt
      else
        DiffEqSensitivity.gradient!(dg_val, pgpu, u, sensealg, pgpu_config)
        @. _wBinv = _wBinv*dg_val/Nt
      end
    else
      if dg_val isa Tuple
        dg[1](dg_val[1],u,uf.p,nothing,i)
        @. _wBinv = -_wBinv*dg_val[1]/Nt
      else
        dg(dg_val,u,uf.p,nothing,i)
        @. _wBinv = -_wBinv*dg_val/Nt
      end
    end
  end
  return nothing
end

function shadow_adjoint(prob::AdjointLSSProblem; t0skip=zero(prob.Δt), t1skip=zero(prob.Δt))
  shadow_adjoint(prob,prob.sensealg,prob.sensealg.alpha,t0skip,t1skip)
end

function shadow_adjoint(prob::AdjointLSSProblem,sensealg::AdjointLSS,alpha::Number,t0skip,t1skip)
  @unpack sol, S, F, Δt, diffcache, h, b, wa, res, g, g0, dg, umid = prob
  @unpack wBinv, B, E = S
  @unpack dg_val, pgpp, pgpp_config, numparams, numindvar, uf, f, f_cache, pJ, pf, paramjac_config = diffcache

  b .= E*h + B*wBinv
  wa .= F\b

  n0 = searchsortedfirst(sol.t, sol.t[1]+t0skip)
  n1 = searchsortedfirst(sol.t, sol.t[end]-t1skip)

  umidres = @view umid[:,n0:n1-1]
  wares = @view wa[(n0-1)*numindvar+1:(n1-1)*numindvar]

  # reset
  res .*= false

  if dg_val isa Tuple
    for (j,u) in enumerate(eachcol(umidres))
      if dg === nothing
        DiffEqSensitivity.gradient!(dg_val[2], pgpp, uf.p, sensealg, pgpp_config)
        @. res += dg_val[2]
      else
        dg[2](dg_val[2],u,uf.p,nothing,n0+j-1)
        @. res -= dg_val[2]
      end
    end
    res ./= (size(umidres)[2])
  end

  for (j,u) in enumerate(eachcol(umidres))
    _wares = @view wares[(j-1)*numindvar+1:j*numindvar]
    if DiffEqBase.has_paramjac(f)
      f.paramjac(pJ, u, uf.p, pf.t)
    else
      pf.u = u
      jacobian!(pJ, pf, uf.p, f_cache, sensealg, paramjac_config)
    end

    res .+= pJ'*_wares
  end

  return res
end

check_for_g(sensealg::Union{ForwardLSS,AdjointLSS},g)=((sensealg.alpha isa Number && g===nothing) && error("Time dilation needs explicit knowledge of g. Either pass `g` as a kwarg or use ForwardLSS/AdjointLSS with windowing."))


text:
# Builders for regions.
push!(r::Region, b::Block) = MLIR.API.mlirRegionAppendOwnedBlock(r, b)
push_block!(r::Region, b::Block) = push!(r, b)

# Builders for blocks.
push!(b::Block, op::Operation) = MLIR.API.mlirBlockAppendOwnedOperation(b, unwrap(op))
push_operation!(b::Block, op::Operation) = push!(b, op)

# Builders for state.
push!(state::OperationState, n::Int, results::Vector{Type}) = MLIR.API.mlirOperationStateAddResults(Ref(state), n, results)
push!(state::OperationState, n::Int, operands::Vector{Value}) = MLIR.API.mlirOperationStateAddOperands(Ref(state), n, operands)
push!(state::OperationState, n::Int, regions::Vector{Region}) = MLIR.API.mlirOperationStateAddOwnedRegions(Ref(state), n, regions)
push!(state::OperationState, n::Int, successors::Vector{Block}) = MLIR.API.mlirOperationStateAddSuccessors(Ref(state), n, successors)
push!(state::OperationState, n::Int, attrs::Vector{NamedAttribute}) = MLIR.API.mlirOperationStateAddAttributes(Ref(state), n, attrs)
push!(state::OperationState, args::Vector) = push!(state, length(args), args)
push!(state::OperationState, arg::T) where T = push!(state, T[arg])


text:
module NUSFSpectra

using Dates, PyCall, JLD2, FileIO, DelimitedFiles

import Base: *, -, +, /

function __init__()
    push!(pyimport("sys")["path"], joinpath(@__DIR__, "../pythonmodules/pySfgProcess/"))
    global dfg = pyimport("dfg")
    global pscalib = pyimport("pscalib")
    global spectrum = pyimport("spectrum")
    global winspec = pyimport("winspec")
end

mutable struct NUSFSpectrum{T<:Number,N} <: AbstractArray{T,N}
    s::Array{T,N}
    wavenumbers::Array{T,1}
    pyobjs::Array{PyObject,1}
    name::String
    dates::Array{DateTime,1}

    """
    Inner constructor. Access the PyObjects data to put in wavelength and signal
    arrays directly into the object.
    """
    function NUSFSpectrum(pyobjs, name, dates)
        firstspectrum = pyobjs[1][:dfgTruncatedSum]
        spectrumlength = length(firstspectrum)
        if length(pyobjs) > 1
            s = Array{eltype(firstspectrum), 2}(undef, spectrumlength, length(pyobjs))
        elseif length(pyobjs) == 1
            s = Array{eltype(firstspectrum), 1}(undef, spectrumlength)
        else
            error("2D Spectra not supported.")
        end

        s[:,1] = firstspectrum
        for i = 2:length(pyobjs)
            s[:,i] = pyobjs[i][:dfgTruncatedSum]
        end

        wavenumbers = pyobjs[1][:fullwn]

        new{eltype(s), ndims(s)}(s, wavenumbers, pyobjs, name, dates)
    end
end

Base.size(s::NUSFSpectrum) = size(s.s)
Base.getindex(s::NUSFSpectrum, i::Int) = getindex(s.s, i)
Base.getindex(s::NUSFSpectrum{T,N}, I::Vararg{Int, N}) where {N,T} = getindex(s.s, I...)
Base.setindex!(s::NUSFSpectrum, v::Number, i::Int) = setindex!(s.s, v, i)
Base.setindex!(s::NUSFSpectrum{T,N}, v::Number, I::Vararg{Int, N}) where {N,T} = setindex!(s.s, v, I...)
Base.copy(s::NUSFSpectrum) = NUSFSpectrum(copy(s.id), copy(s.s))
Base.stride(s::NUSFSpectrum, k) = size(s, 1)
Base.strides(s::NUSFSpectrum) = (1, size(s, 1))
Base.convert(T::Array{Number}, s::NUSFSpectrum) = s.s |> T

+(s::NUSFSpectrum, t::NUSFSpectrum) = s.s + t.s
-(s::NUSFSpectrum, t::NUSFSpectrum) = s.s - t.s
*(s::NUSFSpectrum, a::Number)     = s.s * a
*(a::Number, s::NUSFSpectrum)     = a * s.s
/(s::NUSFSpectrum, a::Number)     = s.s / a
/(a::Number, s::NUSFSpectrum)     = a ./ s.s

function calibrate(filepath, peak1=(2825, 2860), peak2=(3045, 3070))
    calibfile = filepath

    ps = pscalib[:PScalib](calibfile)
    ps[:plot]()

    #starting and ending values of segment to fit
    val1 = peak1[1]
    val2 = peak1[2]

    #fit it with gaussian, calculate shift
    ps[:fitPeak](val1,val2,0)

    #starting and ending values of segment to fit
    val1 = peak2[1]
    val2 = peak2[2]

    #fit it with gaussian, calculate shift
    ps[:fitPeak](val1,val2,1)

    calib = ps[:evaluateShift]()
    println("Shift: $calib wavenumbers")

    path = splitdir(calibfile)[1]

    save(joinpath(path, "calibration.jld2"), "calib", calib)
    println("""Saved calibration to $(joinpath(path, "calibration.jld2"))""")

    return calib
end

function batchprocess(path, calibfile=joinpath(path, "calibration.jld2");
                        trunc=1e-6, doplot=false,
                        savefile=joinpath(path, "data.jld2"))
    (root, dirs, files) = first(walkdir(path))
    filter!(x -> !startswith(x, "."), dirs)
    paths = joinpath.(root, dirs)

    subpaths = Array{String}[]
    for i in 1:length(paths)
        root, dirs, files = first(walkdir(paths[i]))
        push!(subpaths, joinpath.(paths[i], dirs))
    end

    specs = NUSFSpectrum[]
    for p in subpaths
        numruns = length(p)
        dates = Array{DateTime}(undef, numruns)
        pyobjs = Array{PyObject}(undef, numruns)
        spectra = Array{Array{Float64,1},1}(undef, numruns)
        for (i,s) in enumerate(p)
            # p are the arrays of subpaths, so s are the individual runs
            specinfo = processgold(s, calibfile, doplot=doplot, trunc=trunc)
            dates[i] = specinfo[2][1]
            pyobjs[i] = specinfo[3]
        end

        name = split(p[1], "/")[end-1]

        ## Order spectra by date
        i = sortperm(dates)
        pyobjs = pyobjs[i]
        dates = dates[i]

        push!(specs, NUSFSpectrum(pyobjs, name, dates))
    end

    return specs

    splitext(savefile)[end] != ".jld2" && (savefile *= ".jld2")

    save(savefile, "spectra", specs)

    @info "Saved the following spectra to $savefile:"
    for n in getfield.(specs, :name) println(n) end

    return specs
end

function processgold(path, calibrationfile; trunc=1e-6, doplot=false)

    calib = load(calibrationfile, "calib")
    path
    # get all the .SPE files
    filenames = filter!(x -> splitext(x)[2] == ".SPE", readdir(path))
    filepaths = joinpath.(path, filenames)

    #get datetime of the files
    datetimes = Array{DateTime}(undef, length(filepaths))
    for i = 1:length(filepaths)
        f = winspec[:SpeFile](filepaths[i])
        datetimestr = f[:header][:date] * f[:header][:ExperimentTimeLocal]
        datetimes[i] = DateTime(datetimestr, dateformat"dduuuyyyyHHMMSS")
    end

    #create object, loads each sample and background DFG
    gold = spectrum[:Spectrum](path, shift=calib)

    #subtract the appropriate background DFG from each sample DFG
    gold[:subtractBGs]()

    #plot the imported sample DFGs
    doplot && gold[:plotDFGs]()

    #plot the imported background DFGs
    doplot && gold[:plotBGs]()

    #pad the dfgs with zeros so they align and can be summed up
    gold[:padDFGs]()

    #find the indices of where the reference signal falls off to guide truncation
    gold[:findTruncateIndices](trunc)

    #truncate
    gold[:truncateFullDFGs](gold)

    #plot the truncated DFGs
    doplot && gold[:plotTruncatedDFGs]()

    #sum the truncated DFGs
    gold[:sumTruncatedDFGs]()

    #plot the summed, truncated DFGS
    doplot && gold[:plotSumTruncatedDFG]()

    wavenumbers = gold[:fullwn]
    signal = gold[:dfgTruncatedSum]
    filename = "processed_spectrum.csv"
    savepath = joinpath(path, filename)

    writedlm(savepath, [wavenumbers signal])

    return ([wavenumbers signal], datetimes, gold)

end

export NUSFSpectrum, batchprocess, processgold, calibrate

end # module

max_stars_repo_name: MLackner/NUSFSpectra

max_issues_repo_name: MLackner/NUSFSpectra

max_forks_repo_name: MLackner/NUSFSpectra