text:
# using OrdinaryDiffEq
using ModelingToolkit
using DataDrivenDiffEq
using LinearAlgebra, DiffEqSensitivity, Optim
using DiffEqFlux, Flux
using Plots
using DifferentialEquations
gr()

function lotka(du, u, p, t)
    α, β, γ, δ = p
    du[1] = α*u[1] - β*u[2]*u[1]
    du[2] = γ*u[1]*u[2]  - δ*u[2]
end

# Define the experimental parameter
tspan = (0.0f0,3.0f0)
u0 = Float32[0.44249296,4.6280594]
p_ = Float32[1.3, 0.9, 0.8, 1.8]
prob = ODEProblem(lotka, u0,tspan, p_)
solution = solve(prob, Vern7(), abstol=1e-12, reltol=1e-12, saveat = 0.1)

scatter(solution, alpha = 0.25)
plot!(solution, alpha = 0.5)

# Ideal data
X = Array(solution)
# Add noise to the data
println("Generate noisy data")
Xₙ = X + Float32(1e-3)*randn(eltype(X), size(X))

# Define the neueral network which learns L(x, y, y(t-τ))
# Actually, we do not care about overfitting right now, since we want to
# extract the derivative information without numerical differentiation.
L = FastChain(FastDense(2, 32, tanh),FastDense(32, 32, tanh), FastDense(32, 2))
p = initial_params(L)

function dudt_(u, p,t)
    x, y = u
    z = L(u,p)
    [p_[1]*x + z[1],
    -p_[4]*y + z[2]]
end

prob_nn = ODEProblem(dudt_,u0, tspan, p)
sol_nn = concrete_solve(prob_nn, Tsit5(), u0, p, saveat = solution.t)

plot(solution)
plot!(sol_nn)

function predict(θ)
    Array(concrete_solve(prob_nn, Tsit5(), u0, θ, saveat = solution.t,
                         abstol=1e-6, reltol=1e-6,
                         sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP())))
    # @show Array(concrete_solve(prob_nn, Vern7(), u0, θ, saveat = solution.t,
    #                                           abstol=1e-6, reltol=1e-6,
    #                                           sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP())))
  #   tmp_prob = remake(prob_nn,u0=u0,p=θ)
  # tmp_sol =  solve(tmp_prob, Rodas4(), saveat = solution.t,
  #                                           abstol=1e-6, reltol=1e-6,
  #                                           sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP()))
  # # @show size(Array(tmp_sol))
  # pred =  Array(tmp_sol)
end

# No regularisation right now
function loss(θ)
    pred = predict(θ)
    sum(abs2, Xₙ .- pred), pred
end

# Test
loss(p)

const losses = []

callback(θ,l,pred) = begin
    @show l
    push!(losses, l)
    if length(losses)%50==0
        println("Current loss after $(length(losses)) iterations: $(losses[end])")
    end
    false
end

# First train with ADAM for better convergence
res1 = DiffEqFlux.sciml_train(loss, p, ADAM(0.5), cb=callback, maxiters = 100)

# Train with BFGS
res2 = DiffEqFlux.sciml_train(loss, res1.minimizer, BFGS(initial_stepnorm=0.01), cb=callback, maxiters = 10)

println("Final training loss after $(length(losses)) iterations: $(losses[end])")

# Plot the losses
plot(losses, yaxis = :log, xaxis = :log, xlabel = "Iterations", ylabel = "Loss")

# Plot the data and the approximation
NNsolution = predict(res2.minimizer)
# Trained on noisy data vs real solution
plot(solution.t, NNsolution')
plot!(solution.t, X')

# Ideal derivatives
DX = Array(solution(solution.t, Val{1}))

prob_nn2 = ODEProblem(dudt_,u0, tspan, res2.minimizer)
_sol = solve(prob_nn2, Tsit5())
DX_ = Array(_sol(solution.t, Val{1}))

# The learned derivatives
plot(DX')
plot!(DX_')

# Ideal data
L̄ = [-p_[2]*(X[1,:].*X[2,:])';p_[3]*(X[1,:].*X[2,:])']
# Neural network guess
L̂ = L(Xₙ,res2.minimizer)

scatter(L̄')
plot!(L̂')

# Plot the error
scatter(abs.(L̄-L̂)', yaxis = :log)

## Sparse Identification

# Create a Basis
@variables u[1:2]
# Lots of polynomials
polys = Operation[1]

for i ∈ 1:5
    push!(polys, u[1]^i)
    push!(polys, u[2]^i)
    for j ∈ i:5
        if i != j
            push!(polys, (u[1]^i)*(u[2]^j))
            push!(polys, u[2]^i*u[1]^i)
        end
    end
end

# And some other stuff
h = [cos.(u)...; sin.(u)...; polys...]
basis = Basis(h, u)

# Create an optimizer for the SINDy problem
opt = SR3()
# Create the thresholds which should be used in the search process
λ = exp10.(-7:0.1:3)
# Target function to choose the results from; x = L0 of coefficients and L2-Error of the model
g(x) = x[1] < 1 ? Inf : norm(x, 2)
# Test on original data and without further knowledge
println("SINDy on full ideal, unavailable data")
Ψ = SINDy(Xₙ[:, :], DX[:, :], basis, λ, opt, g = g, maxiter = 10000) # Fail
println(Ψ)
print_equations(Ψ)

# Test on ideal derivative data for unknown function ( not available )
println("SINDy on partial ideal, unavailable data")
Ψ = SINDy(Xₙ[:, 1:end], L̄[:, 1:end], basis, λ,opt, g = g, maxiter = 10000) # Succeed
println(Ψ)
print_equations(Ψ)

# Test on uode derivative data
println("SINDy on learned, partial, available data")
Ψ = SINDy(Xₙ[:, 2:end], L̂[:, 2:end], basis, λ,  opt, g = g, maxiter = 10000, normalize = true, denoise = true) # Succeed
println(Ψ)
print_equations(Ψ)

# Extract the parameter
p̂ = parameters(Ψ)
println("First parameter guess : $(p̂)")

# The parameters are a bit off, but the equations are recovered
# Start another SINDy run to get closer to the ground truth
# Create function
unknown_sys = ODESystem(Ψ)
unknown_eq = ODEFunction(unknown_sys)

# Just the equations
b = Basis((u, p, t)->unknown_eq(u, [1.; 1.], t), u)

# Retune for better parameters -> we could also use DiffEqFlux or other parameter estimation tools here.
Ψf = SINDy(Xₙ[:, 2:end], L̂[:, 2:end], b, STRRidge(0.01), maxiter = 100, convergence_error = 1e-18) # Succeed
println(Ψf)
p̂ = parameters(Ψf)
println("Second parameter guess : $(p̂)")

# Create function
recovered_sys = ODESystem(Ψf)
recovered_eq = ODEFunction(recovered_sys)

# Build a ODE for the estimated system
function dudt(du, u, p, t)
    # Add SINDy Term
    α, δ, β, γ = p
    z = recovered_eq(u, [β; γ], t)
    du[1] = α*u[1] + z[1]
    du[2] = -δ*u[2] + z[2]
end

# Create the approximated problem and solution
ps = [p_[[1,4]]; p̂]
approximate_prob = ODEProblem(dudt, u0, tspan, ps)
approximate_solution = solve(approximate_prob, Tsit5(), saveat = 0.01)

# Plot
plot(solution)
plot!(approximate_solution)

## Simulation

# Look at long term prediction
t_long = (0.0, 50.0)
approximate_prob = ODEProblem(dudt, u0, t_long, ps)
approximate_solution_long = solve(approximate_prob, Tsit5(), saveat = 0.1) # Using higher tolerances here results in exit of julia
plot(approximate_solution_long)

true_prob = ODEProblem(lotka, u0, t_long, p_)
true_solution_long = solve(true_prob, Tsit5(), saveat = approximate_solution_long.t)
plot!(true_solution_long)

using JLD2
@save "knowledge_enhanced_NN.jld2" solution recovered_sys approximate_solution_long NNsolution L true_solution_long Xₙ L̄ L̂
@load "knowledge_enhanced_NN.jld2" solution recovered_sys approximate_solution_long NNsolution L true_solution_long Xₙ L̄ L̂

c1 = 3 # RGBA(174/255,192/255,201/255,1) # Maroon
c2 = :orange # RGBA(132/255,159/255,173/255,1) # Red
c3 = :blue # RGBA(255/255,90/255,0,1) # Orange
c4 = :purple # RGBA(153/255,50/255,204/255,1) # Purple

p1 = plot(0.1:0.1:tspan[end],abs.(Array(solution)[:,2:end] .- NNsolution[:,2:end])' .+ eps(Float32),
          lw = 3, yaxis = :log, title = "Timeseries of UODE Error",
          color = [3 :orange], xlabel = "t",
          label = ["x(t)" "y(t)"],
          titlefont = "Helvetica", legendfont = "Helvetica",
          legend = :topright)

# Plot L₂
p2 = plot(X[1,:], X[2,:], L̂[2,:], lw = 3,
     title = "Neural Network Fit of U2(t)", color = c1,
     label = "Neural Network", xaxis = "x", yaxis="y",
     titlefont = "Helvetica", legendfont = "Helvetica",
     legend = :bottomright)
plot!(X[1,:], X[2,:], L̄[2,:], lw = 3, label = "True Missing Term", color=c2)

p3 = scatter(solution, color = [c1 c2], label = ["x data" "y data"],
             title = "Extrapolated Fit From Short Training Data",
             titlefont = "Helvetica", legendfont = "Helvetica",
             markersize = 5)

plot!(p3,true_solution_long, color = [c1 c2], linestyle = :dot, lw=5, label = ["True x(t)" "True y(t)"])
plot!(p3,approximate_solution_long, color = [c3 c4], lw=1, label = ["Estimated x(t)" "Estimated y(t)"])
plot!(p3,[2.99,3.01],[0.0,maximum(hcat(Array(solution),Array(approximate_solution_long)))],lw=2,color=:black)
annotate!([(1.5,9,text("Training \nData", 10, :center, :top, :black, "Helvetica"))])
l = @layout [grid(1,2)
             grid(1,1)]
plot(p1,p2,p3,layout = l)

savefig("sindy_extrapolation.pdf")


text:
using Test, UUIDs
using PackageAnalyzer
using PackageAnalyzer: parse_project, RegistryEntry
using JLLWrappers
using GitHub

get_libpath() = get(ENV, JLLWrappers.LIBPATH_env, nothing)
const orig_libpath = get_libpath()

const auth = GitHub.AnonymousAuth()

@testset "PackageAnalyzer" begin
    general = general_registry()
    @test isdir(general)
    @test all(p -> isdir(p.path), find_packages())
    @test find_package("julia") ∉ find_packages()
    @test all(p -> isdir(p.path), find_packages("Flux"))
    @test isdir(find_package("Flux").path)
    # Test some properties of the `Measurements` package.  NOTE: they may change
    # in the future!
    measurements = analyze(RegistryEntry(joinpath(general, "M", "Measurements")); auth)
    @test measurements.uuid == UUID("eff96d63-e80a-5855-80a2-b1b0885c5ab7")
    @test measurements.reachable
    @test measurements.docs
    @test measurements.runtests
    @test !measurements.buildkite
    @test !isempty(measurements.lines_of_code)
    # Test results of a couple of packages.  Same caveat as above
    packages = [RegistryEntry(joinpath(general, p...)) for p in (("C", "Cuba"), ("P", "PolynomialRoots"))]
    @test Set(packages) == Set(find_packages("Cuba", "PolynomialRoots")) == Set(find_packages(["Cuba", "PolynomialRoots"]))
    @test packages ⊆ find_packages()
    results = analyze(packages; auth)
    cuba, polyroots = results
    @test length(filter(p -> p.reachable, results)) == 2
    @test length(filter(p -> p.runtests, results)) == 2
    @test cuba.drone
    @test !polyroots.docs # Documentation is in the README!
    # We can also use broadcasting!
    @test Set(results) == Set(analyze.(packages; auth))

    # Test `analyze!` directly
    mktempdir() do root
        measurements2 = analyze!(root, RegistryEntry(joinpath(general, "M", "Measurements")); auth)
        @test isequal(measurements, measurements2)
        @test isdir(joinpath(root, "eff96d63-e80a-5855-80a2-b1b0885c5ab7")) # not cleaned up yet
    end
end

@testset "`analyze`" begin
    for pkg in (analyze(PackageAnalyzer; auth), analyze("https://github.com/giordano/PackageAnalyzer.jl"; auth), analyze(joinpath(@__DIR__, ".."); auth))
        @test pkg.uuid == UUID("e713c705-17e4-4cec-abe0-95bf5bf3e10c")
        @test pkg.reachable == true # default
        @test pkg.docs == true
        @test pkg.runtests == true # here we are!
        @test pkg.github_actions == true
        @test length(pkg.license_files) == 1
        @test pkg.license_files[1].licenses_found == ["MIT"]
        @test pkg.license_files[1].license_filename == "LICENSE"
        @test pkg.license_files[1].license_file_percent_covered > 90
        @test pkg.license_files isa Vector{<:NamedTuple}
        @test keys(pkg.license_files[1]) == (:license_filename, :licenses_found, :license_file_percent_covered)
        @test isempty(pkg.licenses_in_project)
        @test !isempty(pkg.lines_of_code)
        @test pkg.lines_of_code isa Vector{<:NamedTuple}
        @test keys(pkg.lines_of_code[1]) == (:directory, :language, :sublanguage, :files, :code, :comments, :blanks)
        idx = findfirst(row -> row.directory=="src" && row.language==:Julia && row.sublanguage===nothing, pkg.lines_of_code)
        @test idx !== nothing
        @test pkg.lines_of_code[idx].code > 200
    end

    # the tests folder isn't a package!
    # But this helps catch issues in error paths for when things go wrong
    bad_pkg = analyze("."; auth)
    @test bad_pkg.repo == ""
    @test bad_pkg.uuid == UUID(UInt128(0))
    @test !bad_pkg.cirrus
    @test isempty(bad_pkg.license_files)
    @test isempty(bad_pkg.licenses_in_project)

    # The argument is a package name
    pkg = analyze("Pluto"; auth)
    # Just make sure we got the UUID correctly and some statistics are collected.
    @test pkg.uuid == UUID("c3e4b0f8-55cb-11ea-2926-15256bba5781")
    @test !isempty(pkg.license_files)
    @test !isempty(pkg.lines_of_code)
    # The argument looks like a package name but it isn't a registered package
    @test_logs (:error, r"Could not find package in registry") match_mode=:any @test_throws ArgumentError analyze("license_in_project"; auth)
end

@testset "`find_packages` with `analyze`" begin
    results = analyze(find_packages("DataFrames", "Flux"); auth) # this method is threaded
    @test results isa Vector{PackageAnalyzer.Package}
    @test length(results) == 2
    # DataFrames currently has 16k LoC; Flux has 5k. Let's check that they aren't mixed up
    # due to some kind of race condition.
    @test results[1].name == "DataFrames"
    @test PackageAnalyzer.count_julia_loc(results[1], "src") > 14000
    @test PackageAnalyzer.count_docs(results[1]) > 5000
    @test PackageAnalyzer.count_readme(results[1]) > 5

    @test results[2].name == "Flux"
    @test PackageAnalyzer.count_julia_loc(results[2].lines_of_code, "src") < 14000

    results = analyze(find_packages("DataFrames"); auth)
    @test results isa Vector{PackageAnalyzer.Package}
    @test length(results) == 1
    @test results[1].name == "DataFrames"

    result = analyze(find_package("DataFrames"); auth)
    @test result isa PackageAnalyzer.Package
    @test result.name == "DataFrames"

    # Don't error when not finding stdlibs in `find_packages`
    @test_logs find_packages("Dates")
    # But we do emit a warning log for non-existent package `Abc`
    @test_logs (:error,) find_packages("Abc")

    # Check we get a nice error for `find_package`
    @test_throws ArgumentError("Standard library Dates not present in registry") find_package("Dates")
    @test_logs (:error,)  @test_throws ArgumentError("Abc not found in registry") find_package("Abc")
end

@testset "`analyze_path!`" begin
    # we check the error path here; the success path is covered by other tests.
    # This also makes sure trying to clone the repo doesn't prompt for
    # username/password
    result = PackageAnalyzer.analyze_path!(mktempdir(), "https://github.com/giordano/DOES_NOT_EXIST.jl"; auth)
    @test result isa PackageAnalyzer.Package
    @test !result.reachable
    @test isempty(result.name)
end

@testset "`subdir` support" begin
    snoop_core_path = only(find_packages("SnoopCompileCore"))
    snoop_core = analyze(snoop_core_path; auth)
    @test !isempty(snoop_core.subdir)
    @test snoop_core.name == "SnoopCompileCore" # this test would fail if we were parsing the wrong Project.toml (that of SnoopCompile)
    # This tests that we find licenses in the subdir and put them first,
    # and find licenses in the repo dir and put them last.
    @test startswith(snoop_core.license_files[1].license_filename, "SnoopCompileCore")
    @test !startswith(snoop_core.license_files[end].license_filename, "SnoopCompileCore")

    # Let's check we're counting LoC right for subdirectories. We have two ways of counting
    # the Julia code in SnoopCompileCore. One, we can ask SnoopCompile for the lines of Julia code in its
    # top-level directory `SnoopCompileCore`. Two, we can ask SnoopCompileCore for all of it's
    # Julia code.
    snoop_compile = analyze(find_package("SnoopCompile"); auth)
    snoop_compile_count_for_core = sum(row.code for row in snoop_compile.lines_of_code if row.language == :Julia && row.sublanguage===nothing && row.directory == "SnoopCompileCore")
    snoop_core_count = sum(row.code for row in snoop_core.lines_of_code if row.language == :Julia && row.sublanguage===nothing)
    @test snoop_core_count == snoop_compile_count_for_core

    # this package doesn't exist in the repo anymore; let's ensure it doesn't throw
    snoop_compile_analysis = analyze(find_package("SnoopCompileAnalysis"); auth)
    @test isempty(snoop_compile_analysis.lines_of_code)
end

@testset "`parse_project`" begin
    bad_project = (; name="Invalid Project.toml", uuid=UUID(UInt128(0)), licenses_in_project=String[])
    # malformatted TOML file
    @test parse_project("missingquote") == bad_project

    # bad UUID
    @test parse_project("baduuid") == bad_project

    # non-existent folder
    @test parse_project("rstratarstra") == bad_project

    # proper Project.toml
    this_project = (; name = "PackageAnalyzer", uuid = UUID("e713c705-17e4-4cec-abe0-95bf5bf3e10c"), licenses_in_project = String[])
    @test parse_project(joinpath(@__DIR__, "..")) == this_project

    # has `license = "MIT"`
    project_1 = (; name = "PackageAnalyzer", uuid = UUID("e713c705-17e4-4cec-abe0-95bf5bf3e10c"), licenses_in_project=["MIT"])
    @test parse_project("license_in_project") == project_1

    # has `license = ["MIT", "GPL"]`
    project_2 = (; name = "PackageAnalyzer", uuid = UUID("e713c705-17e4-4cec-abe0-95bf5bf3e10c"), licenses_in_project=["MIT", "GPL"])
    @test parse_project("licenses_in_project") == project_2
end

@testset "`show`" begin
    # this is mostly to test that `show` doesn't error
    str = sprint(show, analyze(pkgdir(PackageAnalyzer); auth))
    @test occursin("* uuid: e713c705-17e4-4cec-abe0-95bf5bf3e10c", str)
    @test occursin("* OSI approved: true", str)
end

@testset "Contributors" begin
    if PackageAnalyzer.github_auth() isa GitHub.AnonymousAuth
        @warn "Skipping contributors tests since `PackageAnalyzer.github_auth()` is anonymous"
    else
        pkg = analyze("DataFrames")
        @test pkg.contributors isa Vector{<:NamedTuple}
        @test length(pkg.contributors) > 160 # ==183 right now, and it shouldn't go down...
        @test PackageAnalyzer.count_contributors(pkg) > 150
        @test PackageAnalyzer.count_commits(pkg) > 2000
        @test PackageAnalyzer.count_contributors(pkg; type="Anonymous") > 10
    end
end

@testset "Thread-safety" begin
    # Make sure none of the above commands leaks LD_LIBRARY_PATH.  This test
    # should be executed at the very end of the test suite.
    @test orig_libpath == get_libpath()
end

ext: jl


text:
__yield__(x = nothing) = error("@yield used outside @generator")

function on_yield(ctx)
    @assert length(ctx.args) == 1
    x, = ctx.args
    :($__yield__($(esc(x))) && return)
end

on_yieldfrom(_) = error("@yieldfrom not implemented")

macro generator(ex)
    if isexpr(ex, :function, 2) && isexpr(ex.args[1], :tuple, 1)
        ex.args[1] = Expr(:call, :_, ex.args[1].args...)
        # Update MacroTools.jl after anonymous function is supported:
        # https://github.com/MikeInnes/MacroTools.jl/pull/140
    end
    def = splitdef(ex)

    body0 = def[:body]
    body = ContextualMacros.expandwith(
        __module__,
        body0,
        yield = on_yield,
        yieldfrom = on_yieldfrom,
    )
    def[:body] = body

    funcname = :($Transducers.__foldl__)

    if def[:name] === :_
        @assert isempty(def[:kwargs])
        @assert length(def[:args]) == 1
        @assert @capture(def[:args][1], collection_::typename_)
        traceename = gensym(string(typename, "#tracee"))
        def[:name] = traceename
        return quote
            $Base.@inline $(combinedef(def))
            Base.iterate(x::$typename) = $start_generator($traceename, x)
            Base.iterate(::$typename, state) = state()
            $(define_foldl(__module__, funcname, typename, collection, body0))
            nothing
        end |> esc
    end

    allargs = map([def[:args]; def[:kwargs]]) do x
        if @capture(x, args_...)
            args
        else
            x::Symbol
        end
    end

    structname = gensym(string(def[:name], "#itr"))
    structparams = [gensym("T_$a") for a in allargs]
    structfields = [:($a::$T) for (a, T) in zip(allargs, structparams)]

    traceename = gensym(string(def[:name], "#tracee"))

    def[:body] = :($structname($(allargs...)))
    quote
        struct $structname{$(structparams...)}
            $(structfields...)
        end
        $Base.@inline $traceename($(allargs...)) = $body
        $(combinedef(def))
        $Base.IteratorSize(::Type{<:$structname}) = $Base.SizeUnknown()
        $Base.IteratorEltype(::Type{<:$structname}) = $Base.EltypeUnknown()
        Base.iterate(it::$structname) =
            $start_generator($traceename, $([:(it.$a) for a in allargs]...))
        Base.iterate(::$structname, state) = state()
        $(define_foldl(__module__, funcname, structname, allargs, body0))
        nothing
    end |> esc
end

@dynamo function start_generator(tracee, args...)
    ir = IR(tracee, args...)
    ir === nothing && return :(tracee(args...))
    return make_generator_ir(functional(ir))
end

function _last(xs)
    thing = nothing
    for x in xs
        thing = Some(x)
    end
    return something(thing)
end

_resolveref(x::GlobalRef) = getfield(x.mod, x.name)
_resolveref(@nospecialize(x)) = x
_issameref(@nospecialize(a), @nospecialize(b)) = _resolveref(a) === _resolveref(b)

function make_generator_ir(fir)
    ir = empty(fir)
    for a in arguments(fir)
        argument!(ir, a)
    end
    for (v, st) in fir
        if isexpr(st.expr, :lambda)
            args = st.expr.args
            st = Statement(
                st,
                expr = Expr(:lambda, make_generator_ir(args[1]), args[2:end]...),
            )
        end
        push!(ir, st)
    end

    yields = filter!(collect(ir)) do (_, st)
        isexpr(st.expr, :call) && _issameref(st.expr.args[1], __yield__)
    end
    @assert length(yields) <= 1
    if length(yields) == 1
        # Convert `cond(__yield__(x), t, f)` to `tuple(x, f)`
        let vy, sty, v, st
            (vy, sty), = yields  # `vy = __yield__(x)`
            (v, st) = _last(ir)  # `cond(vy, t, f)`
            @assert _issameref(st.expr.args[1], IRTools.cond)
            @assert st.expr.args[2] == vy
            f = st.expr.args[4]
            ir[v] = xcall(:tuple, sty.expr.args[2], f)
            ir[vy] = Statement(nothing)
            # Note: Throwing away "`t` branch" as it was the dummy
            # branch introduced by `@yield`.
        end
    end

    if isempty(ir)
        return return!(ir, nothing)
    else
        return return!(ir, _last(ir)[1])
    end
end


text:
"""
    FEMMAcoustSurfModule

Module for operations on boundaries of domains to construct system matrices and
system vectors for linear acoustics.
"""
module FEMMAcoustSurfModule

import Base.Complex

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: AbstractFESet, gradN!, nodesperelem, manifdim
import ..MatAcoustFluidModule: MatAcoustFluid, bulkmodulus
import FinEtools.MatModule: massdensity
import FinEtools.IntegDomainModule: IntegDomain, integrationdata, Jacobiansurface
import FinEtools.FieldModule: ndofs, gatherdofnums!, gathervalues_asmat!
import FinEtools.NodalFieldModule: NodalField
import FinEtools.GeneralFieldModule: GeneralField
import FinEtools.AssemblyModule: AbstractSysvecAssembler, AbstractSysmatAssembler, SysmatAssemblerSparseSymm, startassembly!, assemble!, makematrix!, SysmatAssemblerSparse
import FinEtools.FEMMBaseModule: AbstractFEMM
import FinEtools.MatrixUtilityModule: add_mggt_ut_only!, add_nnt_ut_only!, complete_lt!, locjac!
import LinearAlgebra: norm, cross

"""
    FEMMAcoustSurf{S<:AbstractFESet, F<:Function, M, NF<:Function} <: AbstractFEMM

Class for linear acoustics finite element modeling machine.
"""
mutable struct FEMMAcoustSurf{S<:AbstractFESet, F<:Function, M, NF<:Function} <: AbstractFEMM
    integdomain::IntegDomain{S, F} # geometry data
    material::M # material object
    getnormal!::NF # get the  normal to the surface
end

"""
    FEMMAcoustSurf(integdomain::IntegDomain{S, F},  material::M) where {S<:AbstractFESet, F<:Function, M}

Create the FEMM for integrals over the surface. The normal is computed
from the geometry of the surface elements.
"""
function FEMMAcoustSurf(integdomain::IntegDomain{S, F},  material::M) where {S<:AbstractFESet, F<:Function, M}
    function getnormal!(n::FFltVec, loc::FFltMat, J::FFltMat)
        sdim, mdim = size(J);
        if     mdim == 1 # 1-D fe
        	N = [J[2,1],-J[1,1]];
        elseif     mdim == 2 # 2-D fe
        	N = cross(J[:,1],J[:,2])
        else
        	error("Got an incorrect size of tangents");
        end
        N=N/norm(N)
        copyto!(n, N)
        return n;
    end
    return FEMMAcoustSurf(integdomain, material, getnormal!)
end

"""
    acousticABC(self::FEMMAcoustSurf, assembler::A,
      geom::NodalField,
      Pdot::NodalField{T}) where {T<:Number, A<:AbstractSysmatAssembler}

Compute the acoustic ABC (Absorbing Boundary Condition) matrix.

# Arguments
- `self`   =  acoustics model
- `assembler`  =  matrix assembler; must be able to assemble unsymmetric matrix
- `geom` = geometry field
- `Pdot` = rate of the acoustic (perturbation) pressure field
"""
function acousticABC(self::FEMMAcoustSurf, assembler::A, geom::NodalField, Pdot::NodalField{T}) where {T<:Number, A<:AbstractSysmatAssembler}
    fes = self.integdomain.fes
    # Constants
    nfes = count(fes); # number of finite elements in the set
    ndn = ndofs(Pdot); # number of degrees of freedom per node
    nne =  nodesperelem(fes); # number of nodes per element
    sdim =  ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes);     # manifold dimension of the element
    Dedim = ndn*nne;          # dimension of the element matrix
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts, Ns, gradNparams, w, pc  =  integrationdata(self.integdomain);
    # Material
    bulk_modulus  =   bulkmodulus(self.material);
    mass_density  =   massdensity(self.material);
    c  =  sqrt(bulk_modulus/mass_density); # sound speed
    # Prepare assembler and temporaries
    ecoords = fill(zero(FFlt), nne, ndofs(geom)); # array of Element coordinates
    De = fill(zero(FFlt), Dedim, Dedim);                # element matrix -- used as a buffer
    dofnums = fill(zero(FInt), Dedim); # degree of freedom array -- used as a buffer
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- used as a buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- used as a buffer
    startassembly!(assembler, Dedim, Dedim, nfes, Pdot.nfreedofs, Pdot.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom, ecoords, fes.conn[i]);
        fill!(De, 0.0); # Initialize element matrix
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
            Jac = Jacobiansurface(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            ffactor = (Jac/c*w[j])
            add_nnt_ut_only!(De, Ns[j], ffactor)
        end # Loop over quadrature points
        complete_lt!(De)
        gatherdofnums!(Pdot, dofnums, fes.conn[i]);# retrieve degrees of freedom
        assemble!(assembler, De, dofnums, dofnums);# assemble symmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

function acousticABC(self::FEMMAcoustSurf, geom::NodalField, Pdot::NodalField{T}) where {T<:Number}
    # Were we supplied assembler object?  If not make a default.
    assembler  =  SysmatAssemblerSparseSymm();
    return acousticABC(self, assembler, geom, Pdot);
end

"""
    pressure2resultantforce(self::FEMMAcoustSurf, assembler::A,
      geom::NodalField,
      P::NodalField{T},
       Force::Field) where {T<:Number, A<:AbstractSysmatAssembler}

Compute the rectangular coupling matrix that transcribes given pressure
on the surface into the resultant force acting on the surface.

# Arguments
- `self`   =  acoustics model
- `assembler`  =  matrix assembler; must be able to assemble unsymmetric matrix
- `geom` = geometry field
- `P` = acoustic (perturbation) pressure field
- `Force` = field for the force resultant
"""
function pressure2resultantforce(self::FEMMAcoustSurf, assembler::A, geom::NodalField, P::NodalField{T}, Force::GeneralField) where {T<:Number, A<:AbstractSysmatAssembler}
    fes = self.integdomain.fes
    # Constants
    nfes = count(fes); # number of finite elements in the set
    ndn = ndofs(P); # number of degrees of freedom per node
    nne =  nodesperelem(fes); # number of nodes per element
    sdim =  ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes);     # manifold dimension of the element
    edim = ndn*nne;          # dimension of the element matrix
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts, Ns, gradNparams, w, pc  =  integrationdata(self.integdomain);
    transposedNs = [reshape(N, 1, nne) for N in Ns]
    ecoords = fill(zero(FFlt), nne, ndofs(geom)); # array of Element coordinates
    Ge = fill(zero(FFlt), 3, nne); # element coupling matrix -- used as a buffer
    coldofnums = zeros(FInt, 1, edim); # degree of freedom array -- used as a buffer
    rowdofnums = zeros(FInt, 1, 3); # degree of freedom array -- used as a buffer
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- used as a buffer
    n = fill(zero(FFlt), 3) # normal vector -- used as a buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- used as a buffer
    gatherdofnums!(Force, rowdofnums, [1 2 3]);# retrieve degrees of freedom
    startassembly!(assembler, 3, edim, count(fes), 3, P.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom, ecoords, fes.conn[i]);
        fill!(Ge, 0.0); # Initialize element matrix
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
            Jac = Jacobiansurface(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            n = self.getnormal!(n, loc, J);
            ffactor = (Jac*w[j])
            Ge = Ge + (ffactor*n)*transposedNs[j]
        end # Loop over quadrature points
        gatherdofnums!(P, coldofnums, fes.conn[i]);# retrieve degrees of freedom
        assemble!(assembler, Ge, rowdofnums, coldofnums);# assemble unsymmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

function pressure2resultantforce(self::FEMMAcoustSurf, geom::NodalField, P::NodalField{T}, Force::GeneralField) where {T<:Number}
    assembler  =  SysmatAssemblerSparse();
    return pressure2resultantforce(self, assembler, geom, P, Force)
end

"""
    pressure2resultanttorque(self::FEMMAcoustSurf, assembler::A,
      geom::NodalField,
      P::NodalField{T},
      Torque::GeneralField, CG::FFltVec) where {T<:Number, A<:AbstractSysmatAssembler}

Compute the rectangular coupling matrix that transcribes given pressure
on the surface into the resultant torque acting on the surface with respect
to the CG.

# Arguments
- `self`   =  acoustics model
- `assembler`  =  matrix assembler; must be able to assemble unsymmetric matrix
- `geom` = geometry field
- `P` = acoustic (perturbation) pressure field
- `Torque` = field for the torque resultant
"""
function pressure2resultanttorque(self::FEMMAcoustSurf, assembler::A, geom::NodalField, P::NodalField{T}, Torque::GeneralField, CG::FFltVec) where {T<:Number,  A<:AbstractSysmatAssembler}
    fes = self.integdomain.fes
    # Constants
    nfes = count(fes); # number of finite elements in the set
    ndn = ndofs(P); # number of degrees of freedom per node
    nne =  nodesperelem(fes); # number of nodes per element
    sdim =  ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes);     # manifold dimension of the element
    edim = ndn*nne;          # dimension of the element matrix
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts, Ns, gradNparams, w, pc  =  integrationdata(self.integdomain);
    transposedNs = [reshape(N, 1, nne) for N in Ns]
    ecoords = fill(zero(FFlt), nne, ndofs(geom)); # array of Element coordinates
    Ge = fill(zero(FFlt), 3, nne); # element coupling matrix -- used as a buffer
    coldofnums = zeros(FInt, 1, edim); # degree of freedom array -- used as a buffer
    rowdofnums = zeros(FInt, 1, 3); # degree of freedom array -- used as a buffer
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- used as a buffer
    n = fill(zero(FFlt), 3) # normal vector -- used as a buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- used as a buffer
    gatherdofnums!(Torque, rowdofnums, [1 2 3]);# retrieve degrees of freedom
    startassembly!(assembler, 3, edim, count(fes), 3, P.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom, ecoords, fes.conn[i]);
        fill!(Ge, 0.0); # Initialize element matrix
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
            Jac = Jacobiansurface(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            n = self.getnormal!(n, loc, J);
            ffactor = (Jac*w[j])
            Ge = Ge + (ffactor*cross(vec(vec(loc)-CG), n))*transposedNs[j]
        end # Loop over quadrature points
        gatherdofnums!(P, coldofnums, fes.conn[i]);# retrieve degrees of freedom
        assemble!(assembler, Ge, rowdofnums, coldofnums);# assemble unsymmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

function pressure2resultanttorque(self::FEMMAcoustSurf, geom::NodalField, P::NodalField{T}, Torque::GeneralField, CG::FFltVec) where {T<:Number}
    assembler  =  SysmatAssemblerSparse();
    return pressure2resultanttorque(self, assembler, geom, P, Torque, CG)
end

"""
    acousticcouplingpanels(self::FEMMAcoustSurf, geom::NodalField, u::NodalField{T}) where {T}

Compute the acoustic pressure-structure coupling matrix.

The acoustic pressure-nodal force matrix transforms the pressure
distributed along the surface to forces acting on the nodes of the
finite element model. Its transpose transforms displacements (or
velocities, or accelerations) into the normal component of the
displacement (or velocity, or acceleration) along the surface.

# Arguments
- `geom`=geometry field
- `u` = displacement field

!!! note

- `n` = outer normal (pointing into the acoustic medium).
- The pressures along the surface are assumed constant (uniform) along
    each finite element –- panel. The panel pressures are assumed to be
    given the same numbers as the serial numbers of the finite elements in
    the set.
"""
function acousticcouplingpanels(self::FEMMAcoustSurf, assembler::A, geom::NodalField, u::NodalField{T}) where {A<:AbstractSysmatAssembler, T}
    fes = self.integdomain.fes
    # Constants
    nne =  nodesperelem(fes); # number of nodes per element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes);     # manifold dimension of the element
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts, Ns, gradNparams, w, pc  =  integrationdata(self.integdomain);
    transposedNs = [reshape(N, 1, nne) for N in Ns]
    ecoords = fill(zero(FFlt), nne, ndofs(geom)); # array of Element coordinates
    coldofnums = zeros(FInt, 1, 1); # degree of freedom array -- used as a buffer
    rowdofnums = zeros(FInt, 1, sdim*nne); # degree of freedom array -- used as a buffer
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- used as a buffer
    n = fill(zero(FFlt), sdim) # normal vector -- used as a buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- used as a buffer
    Ge = fill(zero(FFlt), sdim*nne, 1); # Element matrix -- used as a buffer
    startassembly!(assembler, size(Ge, 1), size(Ge, 2), count(fes), u.nfreedofs, count(fes));
    for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom, ecoords, fes.conn[i]);
        fill!(Ge, 0.0); # Initialize element matrix
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
            Jac = Jacobiansurface(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            n = self.getnormal!(n, loc, J);
            Ge = Ge + (Jac*w[j])*reshape(reshape(n, sdim, 1)*transposedNs[j], size(Ge, 1), size(Ge, 2))
        end # Loop over quadrature points
        coldofnums[1] = i
        gatherdofnums!(u, rowdofnums, fes.conn[i]);# retrieve degrees of freedom
        assemble!(assembler, Ge, rowdofnums, coldofnums);# assemble unsymmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

function acousticcouplingpanels(self::FEMMAcoustSurf, geom::NodalField, u::NodalField{T}) where {T}
    assembler = SysmatAssemblerSparse(); # The matrix is not symmetric
    return acousticcouplingpanels(self, assembler, geom, u)
end

end


text:
module Issue778

using Gridap

d = mktempdir()

domain = (0,1,0,1,0,1)
cells = (3,3,3)
model = CartesianDiscreteModel(domain,cells,isperiodic=(true,false,true))
model = simplexify(model)
writevtk(model,joinpath(d,"tmp_model_3"))

domain = (0,1,0,1)
cells = (3,3)
model = CartesianDiscreteModel(domain,cells,isperiodic=(true,false))
model = simplexify(model)
writevtk(model,joinpath(d,"tmp_model_2"))

rm(d,recursive=true)

end # module


text:
using DataFrames, Plots, Calculus
pyplot()

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  function euler_ex71(f::Function, y0::Float64, x::FloatRange{Float64})
    y = Vector{Real}(length(x))
    y[1] = y0
    h = x.step/x.divisor
    for i in 2:length(x)
      y[i] = y[i-1] + h * f(x[i-1], y[i-1])
    end
    columns=[]
    push!(columns, collect(x))
    push!(columns, y)
    df = DataFrame(columns, [:x, :y])
  end

  f(x,y) = 3x^2 + 2x
  fy(x) = x^3+x^2+1

  x1 = 1.0:0.25:2.0
  x2 = 1.0:0.1:2.0
  x3 = 0.8:0.1:2.2
  df1 = euler_ex71(f, 3.0, x1)
  df2 = euler_ex71(f, 3.0, x2)

  df1 |> display
  println()
  df2 |> display
  println()

  global p2

  p1 = plot(x3, fy.(x3), ylims=(0.8, 13.0))
  plot!(p1, df1[:x], df1[:y], line=(:dash, 1), marker=(:cross))
  plot!(p1, df2[:x], df2[:y], line=(:dash, 1), marker=(:circle))
  show(p1)
  savefig("euler_ex7.1A.png")
  show(p1)
  gui()
end


text:
function get_solution(mechanism::Mechanism{T}) where {T}
    sol = T[]
    for (i, joints) in enumerate(mechanism.joints)
        λ = joints.impulses[2]
        push!(sol, λ...)
    end
    for (i, body) in enumerate(mechanism.bodies)
        v25 = body.state.vsol[2]
        ω25 = body.state.ωsol[2]
        push!(sol, [v25; ω25]...)
    end
    for (i, contacts) in enumerate(mechanism.contacts)
        s = contacts.impulses_dual[2]
        γ = contacts.impulses[2]
        push!(sol, [s; γ]...)
    end
    return sol
end

function set_solution!(mechanism::Mechanism{T}, sol::AbstractVector) where T
    off = 0
    for (i,joints) in enumerate(mechanism.joints)
        nλ = length(joints)
        λ = sol[off .+ (1:nλ)]; off += nλ
        joints.impulses[2] = λ
    end
    for (i,body) in enumerate(mechanism.bodies)
        nv = 3
        nω = 3
        v25 = sol[off .+ (1:nv)]; off += nv
        ω25 = sol[off .+ (1:nω)]; off += nω
        body.state.vsol[2] = v25
        body.state.ωsol[2] = ω25
    end
    for (i,contacts) in enumerate(mechanism.contacts)
        N = length(contacts)
        N½ = Int(N/2)
        s = sol[off .+ (1:N½)]; off += N½
        γ = sol[off .+ (1:N½)]; off += N½
        contacts.impulses_dual[2] = s
        contacts.impulses[2] = γ
    end
    return nothing
end

function evaluate_residual!(mechanism::Mechanism, data::AbstractVector, sol::AbstractVector)
    system = mechanism.system
    set_data!(mechanism, data)
    set_solution!(mechanism, sol)
    set_entries!(mechanism)
    return full_vector(system)
end

function residual_dimension(mechanism::Mechanism)
    return sum(Vector{Int}(length.(mechanism.joints))) +
    	sum(Vector{Int}(length.(mechanism.bodies))) +
    	sum(Vector{Int}(length.(mechanism.contacts)))
end

function finite_difference_data_jacobian(mechanism::Mechanism, data::AbstractVector,
        sol::AbstractVector; δ = 1e-5, verbose=false)
    mechanism = deepcopy(mechanism)
    Nd = data_dim(mechanism, attjac=false)
    Nr = residual_dimension(mechanism)
    jac = zeros(Nr, Nd)
    set_data!(mechanism, deepcopy(data))
    set_solution!(mechanism, deepcopy(sol))

    for i = 1:Nd
        verbose && println("$i / $ndata")
        datap = deepcopy(data)
        datam = deepcopy(data)
        datap[i] += δ
        datam[i] -= δ
        rp = evaluate_residual!(deepcopy(mechanism), datap, deepcopy(sol))
        rm = evaluate_residual!(deepcopy(mechanism), datam, deepcopy(sol))
        jac[:,i] = (rp - rm) / (2δ)
    end
    return jac
end


text:
"""
    Substitute(; statistic=Impute.defaultstats)

Substitute missing values with a summary statistic over the non-missing values.

# Keyword Arguments
* `statistic`: A summary statistic function to be applied to the non-missing values.
  This function should return a value of the same type as the input data `eltype`.
  If this function isn't passed in then the [`Impute.defaultstats`](@ref) function is used to make
  a best guess.

# Example
```jldoctest
julia> using Statistics; using Impute: Substitute, impute

julia> M = [1.0 2.0 missing missing 5.0; 1.1 2.2 3.3 missing 5.5]
2×5 Matrix{Union{Missing, Float64}}:
 1.0  2.0   missing  missing  5.0
 1.1  2.2  3.3       missing  5.5

julia> impute(M, Substitute(); dims=:rows)
2×5 Matrix{Union{Missing, Float64}}:
 1.0  2.0  2.0  2.0   5.0
 1.1  2.2  3.3  2.75  5.5

julia> impute(M, Substitute(; statistic=mean); dims=:rows)
2×5 Matrix{Union{Missing, Float64}}:
 1.0  2.0  2.66667  2.66667  5.0
 1.1  2.2  3.3      3.025    5.5
```
"""
struct Substitute{F<:Function} <: Imputor
    statistic::F
end

Substitute(; statistic=defaultstats) = Substitute(statistic)

function _impute!(data::AbstractArray{Union{T, Missing}}, imp::Substitute) where T
    mask = .!ismissing.(data)
    # Since most summary statistics will require some data, we throw a debug message when
    # all supplied values are missing
    if any(mask)
        x = imp.statistic(disallowmissing(data[mask]))
        return Base.replace!(data, missing => x)
    else
        @warn(
            "Cannot apply substitution function ($(imp.statistic)) " *
            "when all values are missing"
        )
        return data
    end
end

"""
    WeightedSubstitute(; statistic=Impute.defaultstats, weights)

Substitute missing values with a weighted summary statistic over the non-missing values.

# Keyword Arguments
* `statistic`: A summary statistic function to be applied to the non-missing values.
  This function should return a value of the same type as the input data `eltype`.
  If this function isn't passed in then the [`Impute.defaultstats`](@ref) function is used to make
  a best guess.
* `weights`: A set of statistical weights to pass to the `statistic` function.

# Example
```jldoctest
julia> using Statistics, StatsBase; using Impute: WeightedSubstitute, impute

julia> M = [1.0 2.0 missing missing 5.0; 1.1 2.2 3.3 missing 5.5]
2×5 Matrix{Union{Missing, Float64}}:
 1.0  2.0   missing  missing  5.0
 1.1  2.2  3.3       missing  5.5

julia> wv = weights([0.5, 0.2, 0.3, 0.1, 0.4]);

julia> impute(M, WeightedSubstitute(; weights=wv); dims=:rows)
2×5 Matrix{Union{Missing, Float64}}:
 1.0  2.0  2.75  2.75     5.0
 1.1  2.2  3.3   3.11667  5.5
```
"""
struct WeightedSubstitute{F<:Function, W<:AbstractArray{<:Real}} <: Imputor
    statistic::F
    weights::W
end

function WeightedSubstitute(; statistic=defaultstats, weights)
    return WeightedSubstitute(statistic, weights)
end

function _impute!(data::AbstractArray{Union{T, Missing}}, imp::WeightedSubstitute) where T
    mask = .!ismissing.(data)

    # Since most summary statistics will require some data, we throw a debug message when
    # all supplied values are missing
    if any(mask)
        x = imp.statistic(disallowmissing(data[mask]), imp.weights[mask])
        return Base.replace!(data, missing => x)
    else
        @warn(
            "Cannot apply weighted substitution function ($(imp.statistic)) " *
            "when all values are missing"
        )
        return data
    end
end

# Auxiliary functions defining our default substitution rules

@doc """
    defaultstats(data[, wv])

A set of default substitution rules using either `median` or `mode` based on the `eltype` of
the input `data`. Specific rules are summarized as follows.

* `Bool` elements use `mode`
* `Real` elements use `median`
* `Integer` elements where `nunique(data) / length(data) < 0.25` use `mode`
  (ratings, categorical codings, etc)
* `Integer` elements with mostly unique values use `median`
* `!Number` (non-numeric) elements use `mode` as the safest fallback
""" defaultstats

# If we're operating over Bools then we're probably better off using mode
defaultstats(data::AbstractArray{<:Bool}, args...) = _mode(data, args...)

# If we're operating over Reals then we should probably use the median
defaultstats(data::AbstractArray{<:Real}, args...) = median(data, args...)

# If we're operating over integers with relatively few unique values then we're
# likely dealing with either counts or a categorical coding, so mode is probably
# safer
function defaultstats(data::AbstractArray{T}, args...) where T <: Integer
    threshold = 0.25 * length(data)
    nunique = length(unique(data))
    if nunique < threshold
        return _mode(data, args...)
    else
        return round(T, median(data, args...))
    end
end

# Fallback to mode as many types won't support mean or median anyways
defaultstats(data::AbstractArray, args...) = _mode(data, args...)


text:
#
# Constrain Parameters Sampler
#
# constrains the parameter ranges of the subsampler
#
# note: does not enforce actual parameter, nor estimated parameters, to the constrained range
#

type ConstrainParametersSampler <: ModifyingSampler
	subsampler::Sampler
	paramranges::Vector{Tuple{Float64,Float64}}
	function ConstrainParametersSampler(subsampler::Sampler, constrainprs::Vector{Tuple{Float64,Float64}})
		subsamplerprs = paramranges(subsampler)
		length(subsamplerprs) == length(constrainprs) || error("number of constrained parameter ranges does not match number of subsampler parameters")
		prs = Tuple{Float64,Float64}[]
		for i in 1:length(constrainprs)
			cpr = constrainprs[i]
			cpr[1] <= cpr[2] || error("$(i)th constrained parameter range has lower bound greater than upper bound")
			sspr = subsamplerprs[i]
			pr = (max(sspr[1],cpr[1]), min(sspr[2],cpr[2]))
			pr[1] <= pr[2] || error("$(i)th constrained parameter range lies outside subsampler range")
			push!(prs, pr)
		end
		new(subsampler, prs)
	end
end

paramranges(s::ConstrainParametersSampler) = s.paramranges

function setparams!(s::ConstrainParametersSampler, params)
	checkparamranges(s, params)
	setparams!(s.subsampler, params)
end

minimumsupport(s::ConstrainParametersSampler) = minimumsupport(s.subsampler)


text:
module IdentityCellNumbersTests

using Test
using Gridap
using Gridap.CellValuesGallery

a = [2,3,4,6,1,8,3,5]
cn = CellValueFromArray(a)

l = length(a)
id = IdentityCellNumber(Int,l)
r = collect(1:l)
test_index_cell_number(id,r)

cn2 = reindex(cn,id)

@test cn2 === cn

end # module


text:
using Pkg
pkg"activate ."
pkg"add JuMP"
pkg"add Cbc"
pkg"add Juniper"
pkg"add Ipopt"

using JuMP, Ipopt, Juniper

# Escolhendo o modelo
optimizer = Juniper.Optimizer
nl_solver = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
model = Model(optimizer_with_attributes(optimizer, "nl_solver" => nl_solver))

# Definindo os parâmetros que, por sinal, coloquei aleatórios, pois não sei as entradas
m = 4
n = 4
c = rand(1:4, n,m)
a = rand(0:0, m)
b = rand(1000:1000,n)

# Implementando realmente o modelo
@variable(model, x[i = 1:n, g = 1:m], Bin) #Variável
@NLobjective(model, Max, sum(c[i,j] * x[i,g] * x[j,g] for i = 1:n-1, j = i+1:n, g = 1:m)) #Func. objetivo

#Restrições
@constraint(model, [i = 1:n], sum(x[i,g] for g = 1:m) == 1)
@constraint(model, [i = 1:n, g = 1:m], sum(x[i,g] for g = 1:m) ≤ b[g])
@constraint(model, [i = 1:n, g = 1:m], a[g] <= sum(x[i,g] for g = 1:m))

optimize!(model) # Otimizando o modelo

println("Valor da função objetivo = ", objective_value(model))
round.(value.(x))


text:
using SodShockTube
using Documenter

DocMeta.setdocmeta!(SodShockTube, :DocTestSetup, :(using SodShockTube); recursive=true)

makedocs(;
    modules=[SodShockTube],
    authors="Thomas Marks <marksta@umich.edu> and contributors",
    repo="https://github.com/archermarx/SodShockTube.jl/blob/{commit}{path}#{line}",
    sitename="SodShockTube.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://archermarx.github.io/SodShockTube.jl",
        assets=String[],
    ),
    pages=["Home" => "index.md"],
)

deploydocs(; repo="github.com/archermarx/SodShockTube.jl")

max_stars_repo_name: archermarx/SodShockTube

max_issues_repo_name: archermarx/SodShockTube

max_forks_repo_name: archermarx/SodShockTube


text:
module sNMF

include("sNMF_Exact.jl")
include("sNMF_VB.jl")
include("sNMF_MC.jl")
include("BackwardKernel.jl")

using .sNMF_Exact, .sNMF_VB, .sNMF_MC, .BackwardKernel

export generate, log_marginal
export standard_VB, online_VB
export sParticle, smc_weight, particle_mcmc
export EventQueue, sum, length, eltype, iterate

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
using IGC: year_to_short_string, number_to_char, char_to_number
using IGC: IGCLongFilename, IGCShortFilename, filename

@testset "filename" begin
    @testset "short file name" begin
        @testset "year" begin
            @test year_to_short_string(2019) == "9"
        end

        @testset "number_to_char" begin
            @test number_to_char(0) == '0'
            @test number_to_char(9) == '9'
            @test number_to_char(10) == 'A'
            @test number_to_char(31) == 'V'
            @test number_to_char(35) == 'Z'
        end

        @testset "char_to_number" begin
            @test char_to_number('0') == 0
            @test char_to_number('9') == 9
            @test char_to_number('A') == 10
            @test char_to_number('V') == 31
            @test char_to_number('Z') == 35
        end

        @testset "get short file name" begin
            f = IGCShortFilename(
                Date(2015, 06, 17),
                'C',
                "XXX",
                2
            )
            fname_expected = "56HCXXX2.IGC"
            returned_filename = filename(f)
            @test returned_filename == fname_expected

            parsed_filename = parse(IGCShortFilename, fname_expected, year=2010)
            @test parsed_filename.date == f.date
            @test parsed_filename.manufacturer == f.manufacturer
            @test parsed_filename.serial_id == f.serial_id
            @test parsed_filename.flight_of_day == f.flight_of_day
        end
    end

    @testset "get long file name" begin
        f = IGCLongFilename(
            Date(2015, 06, 17),
            "MMM",
            "XXX",
            2
        )
        fname_expected = "2015-06-17-MMM-XXX-02.IGC"
        returned_filename = filename(f)
        @test returned_filename == fname_expected

        parsed_filename = parse(IGCLongFilename, fname_expected)
        @test parsed_filename.date == f.date
        @test parsed_filename.manufacturer == f.manufacturer
        @test parsed_filename.serial_id == f.serial_id
        @test parsed_filename.flight_of_day == f.flight_of_day
    end
end


text:
function CorrelatedRandomVectorGenerator(arg0::RealMatrix, arg1::jdouble, arg2::NormalizedRandomGenerator)
    return CorrelatedRandomVectorGenerator((RealMatrix, jdouble, NormalizedRandomGenerator), arg0, arg1, arg2)
end

function CorrelatedRandomVectorGenerator(arg0::Vector{jdouble}, arg1::RealMatrix, arg2::jdouble, arg3::NormalizedRandomGenerator)
    return CorrelatedRandomVectorGenerator((Vector{jdouble}, RealMatrix, jdouble, NormalizedRandomGenerator), arg0, arg1, arg2, arg3)
end

function get_generator(obj::CorrelatedRandomVectorGenerator)
    return jcall(obj, "getGenerator", NormalizedRandomGenerator, ())
end

function get_rank(obj::CorrelatedRandomVectorGenerator)
    return jcall(obj, "getRank", jint, ())
end

function get_root_matrix(obj::CorrelatedRandomVectorGenerator)
    return jcall(obj, "getRootMatrix", RealMatrix, ())
end

function next_vector(obj::CorrelatedRandomVectorGenerator)
    return jcall(obj, "nextVector", Vector{jdouble}, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end


text:
using DEBmicroTrait
using CSV, DataFrames
using Roots
using Statistics

df     = CSV.read("/Users/glmarschmann/Data/Zhen/IsogenieGenomes.ecosysguilds.csv", DataFrame)

###################### define functional guilds
################################################################################
# aerobic heterotroph
df_aerob_hetero                 = df[(df.heterotroph_withETC.==1) , :]
Genome_size_aerob_hetero        = convert(Array{Float64}, df_aerob_hetero.bp)
V_cell_aerob_hetero             = DEBmicroTrait.genome_size_to_cell_volume(Genome_size_aerob_hetero)
Min_gen_time_aerob_hetero       = df_aerob_hetero.minimum_generation_time_hours
rrn_copies_aerob_hetero         = DEBmicroTrait.genome_size_to_rRNA_copy_number(Genome_size_aerob_hetero)
Gram_stain_aerob_hetero         = repeat(["+"], size(V_cell_aerob_hetero,1))

# acetotrophic methanogens
df_am = df[(df.acetoclastic_methanogenesis.==1) , :]
Genome_size_am = convert(Array{Float64}, df_am.bp)
V_cell_am = DEBmicroTrait.genome_size_to_cell_volume(Genome_size_am)
Min_gen_time_am = df_am.minimum_generation_time_hours
rrn_copies_am = DEBmicroTrait.genome_size_to_rRNA_copy_number(Genome_size_am)
Gram_stain_am = repeat(["+"], size(V_cell_am,1))
SAV_am = DEBmicroTrait.surface_area_volume_ratio(V_cell_am./1e-18)

# hydrogenotrophic methanogens
df_hmo = df[(df.hydrogenotrophic_methanogenesis.==1) , :]
Genome_size_hmo = convert(Array{Float64}, df_hmo.bp)
V_cell_hmo = DEBmicroTrait.genome_size_to_cell_volume(Genome_size_hmo)
Min_gen_time_hmo = df_hmo.minimum_generation_time_hours
rrn_copies_hmo = DEBmicroTrait.genome_size_to_rRNA_copy_number(Genome_size_hmo)
Gram_stain_hmo = repeat(["+"], size(V_cell_hmo,1))
SAV_hmo = DEBmicroTrait.surface_area_volume_ratio(V_cell_hmo./1e-18)

# methane oxidizers
df_mo = df[(df.methane_oxidation.==1) , :]
Genome_size_mo = convert(Array{Float64}, df_mo.bp)
V_cell_mo = DEBmicroTrait.genome_size_to_cell_volume(Genome_size_mo)
Min_gen_time_mo = df_mo.minimum_generation_time_hours
rrn_copies_mo = DEBmicroTrait.genome_size_to_rRNA_copy_number(Genome_size_mo)
Gram_stain_mo = repeat(["+"], size(V_cell_mo,1))
SAV_mo = DEBmicroTrait.surface_area_volume_ratio(V_cell_mo./1e-18)

# Fermenters
df_ferm = df[(df.heterotroph_withETC.==0) , :]
Genome_size_ferm = convert(Array{Float64}, df_ferm.bp)
V_cell_ferm = DEBmicroTrait.genome_size_to_cell_volume(Genome_size_ferm)
Min_gen_time_ferm = df_ferm.minimum_generation_time_hours
rrn_copies_ferm = DEBmicroTrait.genome_size_to_rRNA_copy_number(Genome_size_ferm)
Gram_stain_ferm = repeat(["+"], size(V_cell_ferm,1))
SAV_ferm = DEBmicroTrait.surface_area_volume_ratio(V_cell_ferm./1e-18)

# aerobic diazotroph
df_aerob_diazo = df[(df.aerobic_diazotroph.==1), :]
Genome_size_aerob_diazo  = convert(Array{Float64}, df_aerob_diazo.bp)
V_cell_aerob_diazo  = DEBmicroTrait.genome_size_to_cell_volume(Genome_size_aerob_diazo)
Min_gen_time_aerob_diazo  = df_aerob_diazo.minimum_generation_time_hours
rrn_copies_aerob_diazo  = DEBmicroTrait.genome_size_to_rRNA_copy_number(Genome_size_aerob_diazo)
Gram_stain_aerob_diazo = repeat(["+"], size(V_cell_aerob_diazo ,1))
SAV_aerob_diazo  = DEBmicroTrait.surface_area_volume_ratio(V_cell_aerob_diazo./1e-18)

# anaerobic diazotroph
df_anaerob_diazo = df[(df.anaerobic_diazotroph.==1), :]
Genome_size_anaerob_diazo  = convert(Array{Float64}, df_anaerob_diazo.bp)
V_cell_anaerob_diazo  = DEBmicroTrait.genome_size_to_cell_volume(Genome_size_anaerob_diazo)
Min_gen_time_anaerob_diazo  = df_anaerob_diazo.minimum_generation_time_hours
rrn_copies_anaerob_diazo  = DEBmicroTrait.genome_size_to_rRNA_copy_number(Genome_size_anaerob_diazo)
Gram_stain_anaerob_diazo = repeat(["+"], size(V_cell_anaerob_diazo ,1))
SAV_aerob_andiazo  = DEBmicroTrait.surface_area_volume_ratio(V_cell_anaerob_diazo./1e-18)

# Nitrate reductase
df_nitrate_red = df[(df.nitrite_nitrate_partial.==1) , :]
Genome_size_nitrate_red = convert(Array{Float64}, df_nitrate_red.bp)
V_cell_nitrate_red = DEBmicroTrait.genome_size_to_cell_volume(Genome_size_nitrate_red)
Min_gen_time_nitrate_red = df_nitrate_red.minimum_generation_time_hours
rrn_copies_nitrate_red = DEBmicroTrait.genome_size_to_rRNA_copy_number(Genome_size_nitrate_red)
Gram_stain_nitrate_red= repeat(["+"], size(V_cell_nitrate_red,1))
SAV_nitrate_red = DEBmicroTrait.surface_area_volume_ratio(V_cell_nitrate_red./1e-18)

# Nitrite reductase
df_nitrite_red = df[(df.nitrite_nitricoxide_partial.==1) , :]
Genome_size_nitrite_red = convert(Array{Float64}, df_nitrite_red.bp)
V_cell_nitrite_red = DEBmicroTrait.genome_size_to_cell_volume(Genome_size_nitrite_red)
Min_gen_time_nitrite_red = df_nitrite_red.minimum_generation_time_hours
rrn_copies_nitrite_red = DEBmicroTrait.genome_size_to_rRNA_copy_number(Genome_size_nitrite_red)
Gram_stain_nitrite_red= repeat(["+"], size(V_cell_nitrite_red,1))
SAV_nitrite_red = DEBmicroTrait.surface_area_volume_ratio(V_cell_nitrite_red./1e-18)

# Nitric oxide reductase
df_nitoxide_red = df[(df.nitricoxide_nitrousoxide_partial.==1) , :]
Genome_size_nitoxide_red  = convert(Array{Float64}, df_nitoxide_red.bp)
V_cell_nitoxide_red  = DEBmicroTrait.genome_size_to_cell_volume(Genome_size_nitoxide_red )
Min_gen_time_nitoxide_red  = df_nitoxide_red.minimum_generation_time_hours
rrn_copies_nitoxide_red  = DEBmicroTrait.genome_size_to_rRNA_copy_number(Genome_size_nitoxide_red )
Gram_stain_nitoxide_red = repeat(["+"], size(V_cell_nitoxide_red ,1))
SAV_nitoxide_red  = DEBmicroTrait.surface_area_volume_ratio(V_cell_nitoxide_red./1e-18)

# Nitrous oxide reductase
df_nitrousoxide_red = df[(df.nitrousoxide_dinitrogen.==1) , :]
Genome_size_nitrousoxide_red  = convert(Array{Float64}, df_nitrousoxide_red.bp)
V_cell_nitrousoxide_red  = DEBmicroTrait.genome_size_to_cell_volume(Genome_size_nitrousoxide_red)
Min_gen_time_nitrousoxide_red  = df_nitrousoxide_red.minimum_generation_time_hours
rrn_copies_nitrousoxide_red  = DEBmicroTrait.genome_size_to_rRNA_copy_number(Genome_size_nitrousoxide_red)
Gram_stain_nitrousoxide_red = repeat(["+"], size(V_cell_nitrousoxide_red ,1))
SAV_nitrousoxide_red  = DEBmicroTrait.surface_area_volume_ratio(V_cell_nitrousoxide_red./1e-18)
################################################################################

################################################################################
                        # Glucose oxidation- aerobic heterotroph
# transporter estimation
N_C  = [6]
y_DE = [2.61]./N_C
y_EM = ones(size(V_cell_aerob_hetero,1))
ρ_ps_aerob_hetero_glucose = zeros(size(V_cell_aerob_hetero,1))
r_ms = zeros(size(V_cell_aerob_hetero,1))
for i in 1:size(V_cell_aerob_hetero,1)
    find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell_aerob_hetero[i]], [Min_gen_time_aerob_hetero[i]], [Gram_stain_aerob_hetero[i]], [rrn_copies_aerob_hetero[i]], y_DE[1], N_C[1], [y_EM[i]])
    ρ_p = Roots.find_zero(find_ρ, 1.0)
    ρ_ps_aerob_hetero_glucose[i] = ρ_p[1]
    r_ms[i] = DEBmicroTrait.check_growth_rate(ρ_p[1], [V_cell_aerob_hetero[i]], [Min_gen_time_aerob_hetero[i]], [Gram_stain_aerob_hetero[i]], [rrn_copies_aerob_hetero[i]], y_DE[1], N_C[1], [y_EM[i]])
end
# respiration
y_EV                           = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies_aerob_hetero)
r_growth                       = @. (1-1/y_EV)*r_ms
r_main                         = DEBmicroTrait.cell_volume_to_specific_maintenance_rate(V_cell_aerob_hetero, Min_gen_time_aerob_hetero, Gram_stain_aerob_hetero)
N_SB_aerob_hetero_glucose      = [DEBmicroTrait.transporter_density_to_monomer_uptake_sites([V_cell_aerob_hetero[i]], ρ_ps_aerob_hetero_glucose[i].*ones(1,1), [Min_gen_time_aerob_hetero[i]], [Gram_stain_aerob_hetero[i]])[1] for i in 1:size(V_cell_aerob_hetero,1)]
k2p                            = 180.0*60^2
r_assim                        = @. y_DE*N_C*k2p*N_SB_aerob_hetero_glucose
r_resp_aerob_hetero_glucose    = r_growth .+ r_main .+r_assim
median(r_resp_aerob_hetero_glucose)
std(r_resp_aerob_hetero_glucose)
# half-saturation constant
Molecular_weight       = [180.16]
D_S                    = DEBmicroTrait.aqueous_diffusivity(Molecular_weight)
K_D_aerob_hetero       = DEBmicroTrait.specific_reference_affinity(V_cell_aerob_hetero, reshape(ρ_ps_aerob_hetero_glucose*100, 714,1), D_S)
median(K_D_aerob_hetero.*N_C.*12.011)
std(K_D_aerob_hetero.*N_C*12.011)
K_D_aerob_hetero_glucose = K_D_aerob_hetero.*N_C.*12.011
                        # Acetate oxidation- aerobic heterotroph
# transporter estimation
N_C  = [2]
y_DE = [0.66]./N_C
y_EM = ones(size(V_cell_aerob_hetero,1))
ρ_ps_aerob_hetero_acetate = zeros(size(V_cell_aerob_hetero,1))
r_ms = zeros(size(V_cell_aerob_hetero,1))
for i in 1:size(V_cell_aerob_hetero,1)
    find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell_aerob_hetero[i]], [Min_gen_time_aerob_hetero[i]], [Gram_stain_aerob_hetero[i]], [rrn_copies_aerob_hetero[i]], y_DE[1], N_C[1], [y_EM[i]])
    ρ_p = Roots.find_zero(find_ρ, 1.0)
    ρ_ps_aerob_hetero_acetate[i] = ρ_p[1]
    r_ms[i] = DEBmicroTrait.check_growth_rate(ρ_p[1], [V_cell_aerob_hetero[i]], [Min_gen_time_aerob_hetero[i]], [Gram_stain_aerob_hetero[i]], [rrn_copies_aerob_hetero[i]], y_DE[1], N_C[1], [y_EM[i]])
end
# respiration
y_EV                           = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies_aerob_hetero)
r_growth                       = @. (1-1/y_EV)*r_ms
r_main                         = DEBmicroTrait.cell_volume_to_specific_maintenance_rate(V_cell_aerob_hetero, Min_gen_time_aerob_hetero, Gram_stain_aerob_hetero)
N_SB_aerob_hetero_acetate      = [DEBmicroTrait.transporter_density_to_monomer_uptake_sites([V_cell_aerob_hetero[i]], ρ_ps_aerob_hetero_acetate[i].*ones(1,1), [Min_gen_time_aerob_hetero[i]], [Gram_stain_aerob_hetero[i]])[1] for i in 1:size(V_cell_aerob_hetero,1)]
k2p                            = 180.0*60^2
r_assim                        = @. y_DE*N_C*k2p*N_SB_aerob_hetero_acetate
r_resp_aerob_hetero_acetate    = r_growth .+ r_main .+r_assim
median(r_resp_aerob_hetero_acetate)
std(r_resp_aerob_hetero_acetate)
# half-saturation constant
Molecular_weight       = [59.04]
D_S                    = DEBmicroTrait.aqueous_diffusivity(Molecular_weight)
K_D_aerob_hetero       = DEBmicroTrait.specific_reference_affinity(V_cell_aerob_hetero, reshape(ρ_ps_aerob_hetero_acetate*100, 714,1), D_S)
median(K_D_aerob_hetero.*N_C.*12.011)
std(K_D_aerob_hetero.*N_C*12.011)
K_D_aerob_hetero_acetate = K_D_aerob_hetero.*N_C.*12.011
################################################################################

################################################################################
                        # Glucose fermentation
N_C = [6]
y_DE = [1.06]./N_C
y_EM = ones(size(V_cell_ferm,1))
ρ_ps_ferm= zeros(size(V_cell_ferm,1))
r_ms = zeros(size(V_cell_ferm,1))
for i in 1:size(V_cell_ferm,1)
    find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell_ferm[i]], [Min_gen_time_ferm[i]], [Gram_stain_ferm[i]], [rrn_copies_ferm[i]], y_DE[1], N_C[1], [y_EM[i]])
    ρ_p = Roots.find_zero(find_ρ, 1.0)
    ρ_ps_ferm[i] = ρ_p[1]
    r_ms[i] = DEBmicroTrait.check_growth_rate(ρ_p[1], [V_cell_ferm[i]], [Min_gen_time_ferm[i]], [Gram_stain_ferm[i]], [rrn_copies_ferm[i]], y_DE[1], N_C[1], [y_EM[i]])
end
# respiration
y_EV = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies_ferm)
r_growth = @. (1-1/y_EV)*r_ms
r_main   = DEBmicroTrait.cell_volume_to_specific_maintenance_rate(V_cell_ferm, Min_gen_time_ferm, Gram_stain_ferm)
N_SB_ferm = [DEBmicroTrait.transporter_density_to_monomer_uptake_sites([V_cell_ferm[i]], ρ_ps_ferm[i]*ones(1,1), [Min_gen_time_ferm[i]], [Gram_stain_ferm[i]])[1] for i in 1:size(V_cell_ferm,1)]
k2p = 1800.0*60^2
r_assim = @. y_DE*N_C*k2p*N_SB_ferm
r_resp_ferm = r_growth .+ r_main .+r_assim
median(r_resp_ferm)
std(r_resp_ferm)
# half saturation constant
Molecular_weight = [180.16]
D_S = DEBmicroTrait.aqueous_diffusivity(Molecular_weight)
K_D_ferm = DEBmicroTrait.specific_reference_affinity(V_cell_ferm, reshape(ρ_ps_ferm*100, 815,1), D_S)
median(K_D_ferm*6*12.011)
std(K_D_ferm*6*12.011)

minimum(K_D_ferm*6*12.011)
maximum(K_D_ferm*6*12.011)
K_D_ferm = K_D_ferm*6*12.011
################################################################################

################################################################################
                        # Acetotrophic methanogens
N_C = [2]
y_DE = [1.49]./N_C
y_EM = ones(size(V_cell_am,1))
ρ_ps_am= zeros(size(V_cell_am,1))
r_ms = zeros(size(V_cell_am,1))
for i in 1:size(V_cell_am,1)
    find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell_am[i]], [Min_gen_time_am[i]], [Gram_stain_am[i]], [rrn_copies_am[i]], y_DE[1], N_C[1], [y_EM[i]])
    ρ_p = Roots.find_zero(find_ρ, 1.0)
    ρ_ps_am[i] = ρ_p[1]
    r_ms[i] = DEBmicroTrait.check_growth_rate(ρ_p[1], [V_cell_am[i]], [Min_gen_time_am[i]], [Gram_stain_am[i]], [rrn_copies_am[i]], y_DE[1], N_C[1], [y_EM[i]])
end
# respiration
y_EV = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies_am)
r_growth = @. (1-1/y_EV)*r_ms
r_main   = DEBmicroTrait.cell_volume_to_specific_maintenance_rate(V_cell_am, Min_gen_time_am, Gram_stain_am)
N_SB_am = [DEBmicroTrait.transporter_density_to_monomer_uptake_sites([V_cell_am[i]], ρ_ps_am[i]*ones(1,1), [Min_gen_time_am[i]], [Gram_stain_am[i]])[1] for i in 1:size(V_cell_am,1)]
k2p = 180.0*60^2
r_assim = @. y_DE*N_C*k2p*N_SB_am
r_resp_am = r_growth .+ r_main .+r_assim
median(r_resp_am)
std(r_resp_am)
mean(r_resp_am[r_resp_am .< 0.3])
median(r_growth)
median(r_assim)
# half saturation constants
Molecular_weight = [59.04]
D_S = DEBmicroTrait.aqueous_diffusivity(Molecular_weight)
K_D_am = DEBmicroTrait.specific_reference_affinity(V_cell_am, reshape(ρ_ps_am*100, 65,1), D_S)
median(K_D_am.*N_C*12.011)
std(K_D_am.*N_C*12.011)
K_D_am = K_D_am.*N_C*12.011
################################################################################

################################################################################
                        # hydrogenotrophic methanogens
N_C = [1]
y_DE = [0.166]./N_C
y_EM = ones(size(V_cell_hmo,1))
ρ_ps_hmo= zeros(size(V_cell_hmo,1))
r_ms = zeros(size(V_cell_hmo,1))
for i in 1:size(V_cell_hmo,1)
    find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell_hmo[i]], [Min_gen_time_hmo[i]], [Gram_stain_hmo[i]], [rrn_copies_hmo[i]], y_DE[1], N_C[1], [y_EM[i]])
    ρ_p = Roots.find_zero(find_ρ, 1.0)
    ρ_ps_hmo[i] = ρ_p[1]
    r_ms[i] = DEBmicroTrait.check_growth_rate(ρ_p[1], [V_cell_hmo[i]], [Min_gen_time_hmo[i]], [Gram_stain_hmo[i]], [rrn_copies_hmo[i]], y_DE[1], N_C[1], [y_EM[i]])
end
# respiration
y_EV = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies_hmo)
r_growth = @. (1-1/y_EV)*r_ms
r_main   = DEBmicroTrait.cell_volume_to_specific_maintenance_rate(V_cell_hmo, Min_gen_time_hmo, Gram_stain_hmo)
N_SB_hmo = [DEBmicroTrait.transporter_density_to_monomer_uptake_sites([V_cell_hmo[i]], ρ_ps_hmo[i]*ones(1,1), [Min_gen_time_hmo[i]], [Gram_stain_hmo[i]])[1] for i in 1:size(V_cell_hmo,1)]
k2p = 180.0*60^2
r_assim = @. y_DE*N_C*k2p*N_SB_hmo
r_resp_hmo = r_growth .+ r_main .+r_assim
median(r_resp_hmo)
std(r_resp_hmo)
# half saturation constant
Molecular_weight = [2.016]
D_S = DEBmicroTrait.aqueous_diffusivity(Molecular_weight)
K_D_hmo = DEBmicroTrait.specific_reference_affinity(V_cell_mo, reshape(ρ_ps_hmo, 69,1), D_S)
median(K_D_hmo.*N_C*12.011)
std(K_D_hmo.*N_C*12.011)

################################################################################
                        # methanotrophs
N_C = [6]
y_DE = [5.83]./N_C
y_EM = ones(size(V_cell_mo,1))

ρ_ps_mo= zeros(size(V_cell_mo,1))
r_ms = zeros(size(V_cell_mo,1))
for i in 1:size(V_cell_mo,1)
    find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell_mo[i]], [Min_gen_time_mo[i]], [Gram_stain_mo[i]], [rrn_copies_mo[i]], y_DE[1], N_C[1], [y_EM[i]])
    ρ_p = Roots.find_zero(find_ρ, 1.0)
    ρ_ps_mo[i] = ρ_p[1]
    r_ms[i] = DEBmicroTrait.check_growth_rate(ρ_p[1], [V_cell_mo[i]], [Min_gen_time_mo[i]], [Gram_stain_mo[i]], [rrn_copies_mo[i]], y_DE[1], N_C[1], [y_EM[i]])
end
# respiration
y_EV = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies_mo)
r_growth = @. (1-1/y_EV)*r_ms
r_main   = DEBmicroTrait.cell_volume_to_specific_maintenance_rate(V_cell_mo, Min_gen_time_mo, Gram_stain_mo)
N_SB_mo = [DEBmicroTrait.transporter_density_to_monomer_uptake_sites([V_cell_mo[i]], ρ_ps_mo[i]*ones(1,1), [Min_gen_time_mo[i]], [Gram_stain_mo[i]])[1] for i in 1:size(V_cell_mo,1)]
k2p = 18.0*60^2
r_assim = @. y_DE*N_C*k2p*N_SB_mo
r_resp_mo = r_growth .+ r_main .+r_assim
median(r_resp_mo)
std(r_resp_mo)
################################################################################################
# affinity
Molecular_weight = [16.043]
D_S = DEBmicroTrait.aqueous_diffusivity(Molecular_weight)
K_D_mo = DEBmicroTrait.specific_reference_affinity(V_cell_mo, reshape(ρ_ps_mo./100, 15,1), D_S)
median(K_D_mo.*1*12.011)
std(K_D_mo*1*12.011)
K_D_mo = K_D_mo*1*12.011
################################################################################
                        # aerobic diazotroph
N_C = [6]
y_DE = [2.61]./N_C
y_EM = ones(size(V_cell_aerob_diazo,1))

ρ_ps_aerob_diazo= zeros(size(V_cell_aerob_diazo,1))
r_ms = zeros(size(V_cell_aerob_diazo,1))
for i in 1:size(V_cell_aerob_diazo,1)
    find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell_aerob_diazo[i]], [Min_gen_time_aerob_diazo[i]], [Gram_stain_aerob_diazo[i]], [rrn_copies_aerob_diazo[i]], y_DE[1], N_C[1], [y_EM[i]])
    ρ_p = Roots.find_zero(find_ρ, 1.0)
    ρ_ps_aerob_diazo[i] = ρ_p[1]
    r_ms[i] = DEBmicroTrait.check_growth_rate(ρ_p[1], [V_cell_aerob_diazo[i]], [Min_gen_time_aerob_diazo[i]], [Gram_stain_aerob_diazo[i]], [rrn_copies_aerob_diazo[i]], y_DE[1], N_C[1], [y_EM[i]])
end
# respiration
y_EV = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies_aerob_diazo)
r_growth = @. (1-1/y_EV)*r_ms
r_main   = DEBmicroTrait.cell_volume_to_specific_maintenance_rate(V_cell_aerob_diazo, Min_gen_time_aerob_diazo, Gram_stain_aerob_diazo)
N_SB_aerob_diazo = [DEBmicroTrait.transporter_density_to_monomer_uptake_sites([V_cell_aerob_diazo[i]], ρ_ps_aerob_diazo[i]*ones(1,1), [Min_gen_time_aerob_diazo[i]], [Gram_stain_aerob_diazo[i]])[1] for i in 1:size(V_cell_aerob_diazo,1)]
k2p = 180.0*60^2
r_assim = @. y_DE*N_C*k2p*N_SB_aerob_diazo
r_resp_aerob_diazo = r_growth .+ r_main .+r_assim
median(r_resp_aerob_diazo)
std(r_resp_aerob_diazo)
################################################################################################
# affinity
Molecular_weight = [180.16]
D_S = DEBmicroTrait.aqueous_diffusivity(Molecular_weight)
K_D_aerob_diazo = DEBmicroTrait.specific_reference_affinity(V_cell_aerob_diazo, reshape(ρ_ps_aerob_diazo, 45,1), D_S)
median(K_D_aerob_diazo.*N_C*12.011)
std(K_D_aerob_diazo*N_C*12.011)
K_D_aerob_diazo=K_D_aerob_diazo*N_C*12.011
################################################################################
                        # anaerobic diazotroph
N_C = [6]
y_DE = [4.36]./N_C
y_EM = ones(size(V_cell_anaerob_diazo,1))

ρ_ps_anaerob_diazo= zeros(size(V_cell_anaerob_diazo,1))
r_ms = zeros(size(V_cell_anaerob_diazo,1))
for i in 1:size(V_cell_anaerob_diazo,1)
    find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell_anaerob_diazo[i]], [Min_gen_time_anaerob_diazo[i]], [Gram_stain_anaerob_diazo[i]], [rrn_copies_anaerob_diazo[i]], y_DE[1], N_C[1], [y_EM[i]])
    ρ_p = Roots.find_zero(find_ρ, 1.0)
    ρ_ps_anaerob_diazo[i] = ρ_p[1]
    r_ms[i] = DEBmicroTrait.check_growth_rate(ρ_p[1], [V_cell_anaerob_diazo[i]], [Min_gen_time_anaerob_diazo[i]], [Gram_stain_anaerob_diazo[i]], [rrn_copies_anaerob_diazo[i]], y_DE[1], N_C[1], [y_EM[i]])
end
# respiration
y_EV = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies_anaerob_diazo)
r_growth = @. (1-1/y_EV)*r_ms
r_main   = DEBmicroTrait.cell_volume_to_specific_maintenance_rate(V_cell_anaerob_diazo, Min_gen_time_anaerob_diazo, Gram_stain_anaerob_diazo)
N_SB_anaerob_diazo = [DEBmicroTrait.transporter_density_to_monomer_uptake_sites([V_cell_anaerob_diazo[i]], ρ_ps_anaerob_diazo[i]*ones(1,1), [Min_gen_time_anaerob_diazo[i]], [Gram_stain_anaerob_diazo[i]])[1] for i in 1:size(V_cell_anaerob_diazo,1)]
k2p = 180.0*60^2
r_assim = @. y_DE*N_C*k2p*N_SB_anaerob_diazo
r_resp_anaerob_diazo = r_growth .+ r_main .+r_assim
median(r_resp_anaerob_diazo)
std(r_resp_anaerob_diazo)
# half saturation constant
Molecular_weight = [180.16]
D_S = DEBmicroTrait.aqueous_diffusivity(Molecular_weight)
K_D_anaerob_diazo = DEBmicroTrait.specific_reference_affinity(V_cell_anaerob_diazo, reshape(ρ_ps_anaerob_diazo, 159,1), D_S)
median(K_D_anaerob_diazo.*N_C*12.011)
std(K_D_anaerob_diazo*N_C*12.011)
K_D_anaerob_diazo=K_D_anaerob_diazo*N_C*12.011
################################################################################

################################################################################
                        # Glucose nitrate reductase
N_C = [6]
y_DE = [2.13]./N_C
y_EM = ones(size(V_cell_nitrate_red,1))

ρ_ps_nitrate_red= zeros(size(V_cell_nitrate_red,1))
r_ms = zeros(size(V_cell_nitrate_red,1))
for i in 1:size(V_cell_nitrate_red,1)
    find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell_nitrate_red[i]], [Min_gen_time_nitrate_red[i]], [Gram_stain_nitrate_red[i]], [rrn_copies_nitrate_red[i]], y_DE[1], N_C[1], [y_EM[i]])
    ρ_p = Roots.find_zero(find_ρ, 1.0)
    ρ_ps_nitrate_red[i] = ρ_p[1]
    r_ms[i] = DEBmicroTrait.check_growth_rate(ρ_p[1], [V_cell_nitrate_red[i]], [Min_gen_time_nitrate_red[i]], [Gram_stain_nitrate_red[i]], [rrn_copies_nitrate_red[i]], y_DE[1], N_C[1], [y_EM[i]])
end
# respiration
y_EV = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies_nitrate_red)
r_growth = @. (1-1/y_EV)*r_ms
r_main   = DEBmicroTrait.cell_volume_to_specific_maintenance_rate(V_cell_nitrate_red, Min_gen_time_nitrate_red, Gram_stain_nitrate_red)
N_SB_nitrate_red = [DEBmicroTrait.transporter_density_to_monomer_uptake_sites([V_cell_nitrate_red[i]], ρ_ps_nitrate_red[i]*ones(1,1), [Min_gen_time_nitrate_red[i]], [Gram_stain_nitrate_red[i]])[1] for i in 1:size(V_cell_nitrate_red,1)]
k2p = 180.0*60^2
r_assim = @. y_DE*N_C*k2p*N_SB_nitrate_red
r_resp_nitrate_red_glucose = r_growth .+ r_main .+r_assim
median(r_resp_nitrate_red_glucose)
std(r_resp_nitrate_red_glucose)
# half saturation constant
Molecular_weight = [180.16]
D_S = DEBmicroTrait.aqueous_diffusivity(Molecular_weight)
K_D_nitrate_red = DEBmicroTrait.specific_reference_affinity(V_cell_nitrate_red, reshape(ρ_ps_nitrate_red, 69,1), D_S)
median(K_D_nitrate_red.*N_C*12.011)
std(K_D_nitrate_red*N_C*12.011)
K_D_nitrate_red_glucose = K_D_nitrate_red*N_C*12.011
################################################################################

################################################################################
                        # Glucose nitrite reductase
N_C = [6]
y_DE = [5.38]./N_C
y_EM = ones(size(V_cell_nitrite_red,1))

ρ_ps_nitrite_red= zeros(size(V_cell_nitrite_red,1))
r_ms = zeros(size(V_cell_nitrite_red,1))
for i in 1:size(V_cell_nitrite_red,1)
    find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell_nitrite_red[i]], [Min_gen_time_nitrite_red[i]], [Gram_stain_nitrite_red[i]], [rrn_copies_nitrite_red[i]], y_DE[1], N_C[1], [y_EM[i]])
    ρ_p = Roots.find_zero(find_ρ, 1.0)
    ρ_ps_nitrite_red[i] = ρ_p[1]
    r_ms[i] = DEBmicroTrait.check_growth_rate(ρ_p[1], [V_cell_nitrite_red[i]], [Min_gen_time_nitrite_red[i]], [Gram_stain_nitrite_red[i]], [rrn_copies_nitrite_red[i]], y_DE[1], N_C[1], [y_EM[i]])
end
# respiration
y_EV = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies_nitrite_red)
r_growth = @. (1-1/y_EV)*r_ms
r_main   = DEBmicroTrait.cell_volume_to_specific_maintenance_rate(V_cell_nitrite_red, Min_gen_time_nitrite_red, Gram_stain_nitrite_red)
N_SB_nitrite_red = [DEBmicroTrait.transporter_density_to_monomer_uptake_sites([V_cell_nitrite_red[i]], ρ_ps_nitrite_red[i]*ones(1,1), [Min_gen_time_nitrite_red[i]], [Gram_stain_nitrite_red[i]])[1] for i in 1:size(V_cell_nitrite_red,1)]
k2p = 180.0*60^2
r_assim = @. y_DE*N_C*k2p*N_SB_nitrite_red
r_resp_nitrite_red_glucose = r_growth .+ r_main .+r_assim
median(r_resp_nitrite_red_glucose)
std(r_resp_nitrite_red_glucose)
# half saturation constant
Molecular_weight = [180.16]
D_S = DEBmicroTrait.aqueous_diffusivity(Molecular_weight)
K_D_nitrite_red = DEBmicroTrait.specific_reference_affinity(V_cell_nitrite_red, reshape(ρ_ps_nitrite_red, 63,1), D_S)
median(K_D_nitrite_red.*N_C*12.011)
std(K_D_nitrite_red*N_C*12.011)
K_D_nitrite_red_glucose = K_D_nitrite_red*N_C*12.011
################################################################################

################################################################################
                        # Glucose nitric oxide reductase
N_C = [6]
y_DE = [5.43]./N_C
y_EM = ones(size(V_cell_nitoxide_red,1))

ρ_ps_nitoxide_red= zeros(size(V_cell_nitoxide_red,1))
r_ms = zeros(size(V_cell_nitoxide_red,1))
for i in 1:size(V_cell_nitoxide_red,1)
    find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell_nitoxide_red[i]], [Min_gen_time_nitoxide_red[i]], [Gram_stain_nitoxide_red[i]], [rrn_copies_nitoxide_red[i]], y_DE[1], N_C[1], [y_EM[i]])
    ρ_p = Roots.find_zero(find_ρ, 1.0)
    ρ_ps_nitoxide_red[i] = ρ_p[1]
    r_ms[i] = DEBmicroTrait.check_growth_rate(ρ_p[1], [V_cell_nitoxide_red[i]], [Min_gen_time_nitoxide_red[i]], [Gram_stain_nitoxide_red[i]], [rrn_copies_nitoxide_red[i]], y_DE[1], N_C[1], [y_EM[i]])
end
# respiration
y_EV = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies_nitoxide_red)
r_growth = @. (1-1/y_EV)*r_ms
r_main   = DEBmicroTrait.cell_volume_to_specific_maintenance_rate(V_cell_nitoxide_red, Min_gen_time_nitoxide_red, Gram_stain_nitoxide_red)
N_SB_nitoxide_red = [DEBmicroTrait.transporter_density_to_monomer_uptake_sites([V_cell_nitoxide_red[i]], ρ_ps_nitoxide_red[i]*ones(1,1), [Min_gen_time_nitoxide_red[i]], [Gram_stain_nitoxide_red[i]])[1] for i in 1:size(V_cell_nitoxide_red,1)]
k2p = 180.0*60^2
r_assim = @. y_DE*N_C*k2p*N_SB_nitoxide_red
r_resp_nitoxide_red_glucose = r_growth .+ r_main .+r_assim
median(r_resp_nitoxide_red_glucose)
std(r_resp_nitoxide_red_glucose)
# half saturation constant
Molecular_weight = [180.16]
D_S = DEBmicroTrait.aqueous_diffusivity(Molecular_weight)
K_D_nitoxide_red = DEBmicroTrait.specific_reference_affinity(V_cell_nitoxide_red, reshape(ρ_ps_nitoxide_red, 200,1), D_S)
median(K_D_nitoxide_red.*N_C*12.011)
std(K_D_nitoxide_red*N_C*12.011)
K_D_nitoxide_red_glucose = K_D_nitoxide_red*N_C*12.011
################################################################################

################################################################################
                        # Glucose nitrous oxide reductase
N_C = [6]
y_DE = [4.19]./N_C
y_EM = ones(size(V_cell_nitrousoxide_red,1))

ρ_ps_nitrousoxide_red= zeros(size(V_cell_nitrousoxide_red,1))
r_ms = zeros(size(V_cell_nitrousoxide_red,1))
for i in 1:size(V_cell_nitrousoxide_red,1)
    find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell_nitrousoxide_red[i]], [Min_gen_time_nitrousoxide_red[i]], [Gram_stain_nitrousoxide_red[i]], [rrn_copies_nitrousoxide_red[i]], y_DE[1], N_C[1], [y_EM[i]])
    ρ_p = Roots.find_zero(find_ρ, 1.0)
    ρ_ps_nitrousoxide_red[i] = ρ_p[1]
    r_ms[i] = DEBmicroTrait.check_growth_rate(ρ_p[1], [V_cell_nitrousoxide_red[i]], [Min_gen_time_nitrousoxide_red[i]], [Gram_stain_nitrousoxide_red[i]], [rrn_copies_nitrousoxide_red[i]], y_DE[1], N_C[1], [y_EM[i]])
end
# respiration
y_EV = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies_nitrousoxide_red)
r_growth = @. (1-1/y_EV)*r_ms
r_main   = DEBmicroTrait.cell_volume_to_specific_maintenance_rate(V_cell_nitrousoxide_red, Min_gen_time_nitrousoxide_red, Gram_stain_nitrousoxide_red)
N_SB_nitrousoxide_red = [DEBmicroTrait.transporter_density_to_monomer_uptake_sites([V_cell_nitrousoxide_red[i]], ρ_ps_nitrousoxide_red[i]*ones(1,1), [Min_gen_time_nitrousoxide_red[i]], [Gram_stain_nitrousoxide_red[i]])[1] for i in 1:size(V_cell_nitrousoxide_red,1)]
k2p = 180.0*60^2
r_assim = @. y_DE*N_C*k2p*N_SB_nitrousoxide_red
r_resp_nitrousoxide_red_glucose = r_growth .+ r_main .+r_assim
median(r_resp_nitrousoxide_red_glucose)
std(r_resp_nitrousoxide_red_glucose)
# half saturation constant
Molecular_weight = [180.16]
D_S = DEBmicroTrait.aqueous_diffusivity(Molecular_weight)
K_D_nitrousoxide_red = DEBmicroTrait.specific_reference_affinity(V_cell_nitrousoxide_red, reshape(ρ_ps_nitrousoxide_red, 37,1), D_S)
median(K_D_nitrousoxide_red.*N_C*12.011)
std(K_D_nitrousoxide_red*N_C*12.011)
K_D_nitrousoxide_red_glucose = K_D_nitrousoxide_red*N_C*12.011
################################################################################

################################################################################
                        # Acetate nitrate reductase
N_C = [2]
y_DE = [0.04]./N_C
y_EM = ones(size(V_cell_nitrate_red,1))

ρ_ps_nitrate_red= zeros(size(V_cell_nitrate_red,1))
r_ms = zeros(size(V_cell_nitrate_red,1))
for i in 1:size(V_cell_nitrate_red,1)
    find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell_nitrate_red[i]], [Min_gen_time_nitrate_red[i]], [Gram_stain_nitrate_red[i]], [rrn_copies_nitrate_red[i]], y_DE[1], N_C[1], [y_EM[i]])
    ρ_p = Roots.find_zero(find_ρ, 1.0)
    ρ_ps_nitrate_red[i] = ρ_p[1]
    r_ms[i] = DEBmicroTrait.check_growth_rate(ρ_p[1], [V_cell_nitrate_red[i]], [Min_gen_time_nitrate_red[i]], [Gram_stain_nitrate_red[i]], [rrn_copies_nitrate_red[i]], y_DE[1], N_C[1], [y_EM[i]])
end
# respiration
y_EV = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies_nitrate_red)
r_growth = @. (1-1/y_EV)*r_ms
r_main   = DEBmicroTrait.cell_volume_to_specific_maintenance_rate(V_cell_nitrate_red, Min_gen_time_nitrate_red, Gram_stain_nitrate_red)
N_SB_nitrate_red = [DEBmicroTrait.transporter_density_to_monomer_uptake_sites([V_cell_nitrate_red[i]], ρ_ps_nitrate_red[i]*ones(1,1), [Min_gen_time_nitrate_red[i]], [Gram_stain_nitrate_red[i]])[1] for i in 1:size(V_cell_nitrate_red,1)]
k2p = 180.0*60^2
r_assim = @. y_DE*N_C*k2p*N_SB_nitrate_red
r_resp_nitrate_red_acetate = r_growth .+ r_main .+r_assim
median(r_resp_nitrate_red_acetate)
std(r_resp_nitrate_red_acetate)
# half saturation constant
Molecular_weight = [59.04]
D_S = DEBmicroTrait.aqueous_diffusivity(Molecular_weight)
K_D_nitrate_red = DEBmicroTrait.specific_reference_affinity(V_cell_nitrate_red, reshape(ρ_ps_nitrate_red, 69,1), D_S)
median(K_D_nitrate_red.*N_C*12.011)
std(K_D_nitrate_red*N_C*12.011)
K_D_nitrate_red_acetate = K_D_nitrate_red*N_C*12.011
################################################################################
                        # Acetate nitrite reductase
N_C = [2]
y_DE = [0.35]./N_C
y_EM = ones(size(V_cell_nitrite_red,1))

ρ_ps_nitrite_red= zeros(size(V_cell_nitrite_red,1))
r_ms = zeros(size(V_cell_nitrite_red,1))
for i in 1:size(V_cell_nitrite_red,1)
    find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell_nitrite_red[i]], [Min_gen_time_nitrite_red[i]], [Gram_stain_nitrite_red[i]], [rrn_copies_nitrite_red[i]], y_DE[1], N_C[1], [y_EM[i]])
    ρ_p = Roots.find_zero(find_ρ, 1.0)
    ρ_ps_nitrite_red[i] = ρ_p[1]
    r_ms[i] = DEBmicroTrait.check_growth_rate(ρ_p[1], [V_cell_nitrite_red[i]], [Min_gen_time_nitrite_red[i]], [Gram_stain_nitrite_red[i]], [rrn_copies_nitrite_red[i]], y_DE[1], N_C[1], [y_EM[i]])
end
# respiration
y_EV = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies_nitrite_red)
r_growth = @. (1-1/y_EV)*r_ms
r_main   = DEBmicroTrait.cell_volume_to_specific_maintenance_rate(V_cell_nitrite_red, Min_gen_time_nitrite_red, Gram_stain_nitrite_red)
N_SB_nitrite_red = [DEBmicroTrait.transporter_density_to_monomer_uptake_sites([V_cell_nitrite_red[i]], ρ_ps_nitrite_red[i]*ones(1,1), [Min_gen_time_nitrite_red[i]], [Gram_stain_nitrite_red[i]])[1] for i in 1:size(V_cell_nitrite_red,1)]
k2p = 180.0*60^2
r_assim = @. y_DE*N_C*k2p*N_SB_nitrite_red
r_resp_nitrite_red_acetate = r_growth .+ r_main .+r_assim
median(r_resp_nitrite_red_acetate)
std(r_resp_nitrite_red_acetate)
# half saturation constant
Molecular_weight = [59.04]
D_S = DEBmicroTrait.aqueous_diffusivity(Molecular_weight)
K_D_nitrite_red = DEBmicroTrait.specific_reference_affinity(V_cell_nitrite_red, reshape(ρ_ps_nitrite_red, 63,1), D_S)
median(K_D_nitrite_red.*N_C*12.011)
std(K_D_nitrite_red*N_C*12.011)
K_D_nitrite_red_acetate = K_D_nitrite_red*N_C*12.011
################################################################################

################################################################################
                        # Acetate nitric oxide reductase
N_C = [2]
y_DE = [0.88]./N_C
y_EM = ones(size(V_cell_nitoxide_red,1))

ρ_ps_nitoxide_red= zeros(size(V_cell_nitoxide_red,1))
r_ms = zeros(size(V_cell_nitoxide_red,1))
for i in 1:size(V_cell_nitoxide_red,1)
    find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell_nitoxide_red[i]], [Min_gen_time_nitoxide_red[i]], [Gram_stain_nitoxide_red[i]], [rrn_copies_nitoxide_red[i]], y_DE[1], N_C[1], [y_EM[i]])
    ρ_p = Roots.find_zero(find_ρ, 1.0)
    ρ_ps_nitoxide_red[i] = ρ_p[1]
    r_ms[i] = DEBmicroTrait.check_growth_rate(ρ_p[1], [V_cell_nitoxide_red[i]], [Min_gen_time_nitoxide_red[i]], [Gram_stain_nitoxide_red[i]], [rrn_copies_nitoxide_red[i]], y_DE[1], N_C[1], [y_EM[i]])
end
# respiration
y_EV = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies_nitoxide_red)
r_growth = @. (1-1/y_EV)*r_ms
r_main   = DEBmicroTrait.cell_volume_to_specific_maintenance_rate(V_cell_nitoxide_red, Min_gen_time_nitoxide_red, Gram_stain_nitoxide_red)
N_SB_nitoxide_red = [DEBmicroTrait.transporter_density_to_monomer_uptake_sites([V_cell_nitoxide_red[i]], ρ_ps_nitoxide_red[i]*ones(1,1), [Min_gen_time_nitoxide_red[i]], [Gram_stain_nitoxide_red[i]])[1] for i in 1:size(V_cell_nitoxide_red,1)]
k2p = 180.0*60^2
r_assim = @. y_DE*N_C*k2p*N_SB_nitoxide_red
r_resp_nitoxide_red_acetate = r_growth .+ r_main .+r_assim
median(r_resp_nitoxide_red_acetate)
std(r_resp_nitoxide_red_acetate)
# half saturation constant
Molecular_weight = [59.04]
D_S = DEBmicroTrait.aqueous_diffusivity(Molecular_weight)
K_D_nitoxide_red = DEBmicroTrait.specific_reference_affinity(V_cell_nitoxide_red, reshape(ρ_ps_nitoxide_red, 200,1), D_S)
median(K_D_nitoxide_red.*N_C*12.011)
std(K_D_nitoxide_red*N_C*12.011)
K_D_nitoxide_red_acetate= K_D_nitoxide_red*N_C*12.011
################################################################################

################################################################################
                        # Acetate nitrous oxide reductase
N_C = [2]
y_DE = [1.14]./N_C
y_EM = ones(size(V_cell_nitrousoxide_red,1))

ρ_ps_nitrousoxide_red= zeros(size(V_cell_nitrousoxide_red,1))
r_ms = zeros(size(V_cell_nitrousoxide_red,1))
for i in 1:size(V_cell_nitrousoxide_red,1)
    find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell_nitrousoxide_red[i]], [Min_gen_time_nitrousoxide_red[i]], [Gram_stain_nitrousoxide_red[i]], [rrn_copies_nitrousoxide_red[i]], y_DE[1], N_C[1], [y_EM[i]])
    ρ_p = Roots.find_zero(find_ρ, 1.0)
    ρ_ps_nitrousoxide_red[i] = ρ_p[1]
    r_ms[i] = DEBmicroTrait.check_growth_rate(ρ_p[1], [V_cell_nitrousoxide_red[i]], [Min_gen_time_nitrousoxide_red[i]], [Gram_stain_nitrousoxide_red[i]], [rrn_copies_nitrousoxide_red[i]], y_DE[1], N_C[1], [y_EM[i]])
end
# respiration
y_EV = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies_nitrousoxide_red)
r_growth = @. (1-1/y_EV)*r_ms
r_main   = DEBmicroTrait.cell_volume_to_specific_maintenance_rate(V_cell_nitrousoxide_red, Min_gen_time_nitrousoxide_red, Gram_stain_nitrousoxide_red)
N_SB_nitrousoxide_red = [DEBmicroTrait.transporter_density_to_monomer_uptake_sites([V_cell_nitrousoxide_red[i]], ρ_ps_nitrousoxide_red[i]*ones(1,1), [Min_gen_time_nitrousoxide_red[i]], [Gram_stain_nitrousoxide_red[i]])[1] for i in 1:size(V_cell_nitrousoxide_red,1)]
k2p = 180.0*60^2
r_assim = @. y_DE*N_C*k2p*N_SB_nitrousoxide_red
r_resp_nitrousoxide_red_acetate = r_growth .+ r_main .+r_assim
median(r_resp_nitrousoxide_red_acetate)
std(r_resp_nitrousoxide_red_acetate)
# half saturation constant
Molecular_weight = [59.04]
D_S = DEBmicroTrait.aqueous_diffusivity(Molecular_weight)
K_D_nitrousoxide_red = DEBmicroTrait.specific_reference_affinity(V_cell_nitrousoxide_red, reshape(ρ_ps_nitrousoxide_red, 37,1), D_S)
median(K_D_nitrousoxide_red.*N_C*12.011)
std(K_D_nitrousoxide_red*N_C*12.011)
K_D_nitrousoxide_red_acetate = K_D_nitrousoxide_red*N_C*12.011
################################################################################

df_out = DataFrame()

df_out.genome = vcat(df_aerob_hetero.genome, df_aerob_hetero.genome, df_ferm.genome, df_am.genome, df_hmo.genome, df_mo.genome, df_aerob_diazo.genome, df_anaerob_diazo.genome)
df_out.habitat = vcat(df_aerob_hetero.habitat_binnedfrom, df_aerob_hetero.habitat_binnedfrom, df_ferm.habitat_binnedfrom, df_am.habitat_binnedfrom, df_hmo.habitat_binnedfrom, df_mo.habitat_binnedfrom, df_aerob_diazo.habitat_binnedfrom, df_anaerob_diazo.habitat_binnedfrom)
df_out.rresp = vcat(r_resp_aerob_hetero_glucose, r_resp_aerob_hetero_acetate, r_resp_ferm, r_resp_am, r_resp_hmo, r_resp_mo, r_resp_aerob_diazo, r_resp_anaerob_diazo)
df_out.K = vcat(K_D_aerob_hetero_glucose[:,1], K_D_aerob_hetero_acetate[:,1], K_D_ferm[:,1], K_D_am[:,1], K_D_hmo[:,1], K_D_mo[:,1], K_D_aerob_diazo[:,1], K_D_anaerob_diazo[:,1])

functional_groups = vcat(repeat(["Aerob_hetero_glucose"],size(df_aerob_hetero.genome,1)), repeat(["Aerob_hetero_acetate"],size(df_aerob_hetero.genome,1)), repeat(["Fermenter"],size(df_ferm.genome,1)), repeat(["Acetotrophic methanogen"],size(df_am.genome,1)),
                        repeat(["Hydrogenotrophic methanogen"],size(df_hmo.genome,1)), repeat(["Methanotroph"],size(df_mo.genome,1)), repeat(["Aerobic diazotroph"],size(df_aerob_diazo.genome,1)), repeat(["Anaerobic diazotroph"],size(df_anaerob_diazo.genome,1)))

df_out.functionalgroup = functional_groups

CSV.write("/Users/glmarschmann/.julia/dev/DEBmicroTrait/zhen/TraitPredictionsIsogenieGenomes.csv", df_out)

df_out = DataFrame()
df_out.genome = vcat(df_nitrate_red.genome, df_nitrate_red.genome, df_nitrite_red.genome, df_nitrite_red.genome, df_nitoxide_red.genome, df_nitoxide_red.genome, df_nitrousoxide_red.genome, df_nitrousoxide_red.genome)
df_out.habitat = vcat(df_nitrate_red.habitat_binnedfrom, df_nitrate_red.habitat_binnedfrom, df_nitrite_red.habitat_binnedfrom, df_nitrite_red.habitat_binnedfrom, df_nitoxide_red.habitat_binnedfrom, df_nitoxide_red.habitat_binnedfrom, df_nitrousoxide_red.habitat_binnedfrom, df_nitrousoxide_red.habitat_binnedfrom)
df_out.rresp = vcat(r_resp_nitrate_red_glucose, r_resp_nitrate_red_acetate, r_resp_nitrite_red_glucose, r_resp_nitrite_red_acetate, r_resp_nitoxide_red_glucose, r_resp_nitoxide_red_acetate, r_resp_nitrousoxide_red_glucose, r_resp_nitrousoxide_red_acetate)
df_out.K = vcat(K_D_nitrate_red_glucose, K_D_nitrate_red_acetate, K_D_nitrite_red_glucose, K_D_nitrite_red_acetate, K_D_nitoxide_red_glucose, K_D_nitoxide_red_acetate, K_D_nitrousoxide_red_glucose, K_D_nitrousoxide_red_acetate)

functional_groups = vcat(repeat(["Glucose_Nitrate"],size(df_nitrate_red,1)), repeat(["Acetate_Nitrate"],size(df_nitrate_red,1)), repeat(["Glucose_Nitrite"],size(df_nitrite_red,1)), repeat(["Acetate_Nitrite"],size(df_nitrite_red,1)), repeat(["Glucose_Nitoxide"],size(df_nitoxide_red,1)),
                        repeat(["Acetate_Nitoxide"],size(df_nitoxide_red,1)), repeat(["Glucose_Nitrous"],size(df_nitrousoxide_red,1)), repeat(["Glucose_Nitrous"],size(df_nitrousoxide_red,1)))

df_out.functionalgroup = functional_groups
CSV.write("/Users/glmarschmann/.julia/dev/DEBmicroTrait/emerge/TraitPredictionsIsogenieGenomes_Denitrifiers.csv", df_out)

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
module Day8

export day8

using CombinedParsers
import CombinedParsers.Regexp: word, whitespace_maybe

mutable struct Instruction
    instr::String
    value::Int
end

Instruction(t::NamedTuple{(:instr,:value), Tuple{String,Int}}) =
    Instruction(t[:instr], t[:value])

function parseline(str)
    instruction = Sequence(2, CombinedParsers.Regexp.whitespace_maybe, CombinedParsers.Regexp.word, CombinedParsers.Regexp.whitespace_maybe)

    line = Sequence(:instr => instruction, :value => Numeric(Int), CombinedParsers.Regexp.whitespace_maybe)

    str |> line
end

function day8(::Val{:a}, filename)
    code = Vector{NamedTuple{(:instr,:value),Tuple{String,Int}}}()
    for line in readlines(filename)
        push!(code, parseline(line))
    end
    visited_lines = Set{Int}()
    current_line = 1
    acc = 0
    while !in(current_line, visited_lines)
        t = code[current_line]
        push!(visited_lines, current_line)
        if t[:instr] == "nop"
            current_line += 1
        elseif t[:instr] == "acc"
            current_line += 1
            acc += t[:value]
        elseif t[:instr] == "jmp"
            current_line += t[:value]
        else
            println("Unknown instruction!")
        end
    end
    acc
end

function day8(::Val{:b}, filename)
    code = Vector{Instruction}()
    for line in readlines(filename)
        t::NamedTuple{(:instr, :value), Tuple{String, Int}} = parseline(line)
        push!(code, Instruction(t))
    end

    for i in initrun(code)
        if code[i].instr == "acc"
            continue
        elseif code[i].instr == "jmp"
            code[i].instr = "nop"
            val = runcode(code)
            if !isnothing(val)
                return val
            else
                code[i].instr = "jmp"
            end
        else
            code[i].instr = "jmp"
            val = runcode(code)
            if !isnothing(val)
                return val
            else
                code[i].instr = "nop"
            end
        end
    end
end

function runcode(code)
    visited_lines = Set{Int}()
    current_line = 1
    last_line = length(code)
    acc = 0
    while !in(current_line, visited_lines)
        t = code[current_line]
        push!(visited_lines, current_line)
        if t.instr == "nop"
            current_line += 1
        elseif t.instr == "acc"
            current_line += 1
            acc += t.value
        elseif t.instr == "jmp"
            current_line += t.value
        else
            println("Unknown instruction!")
        end
        if current_line > last_line
            return acc
        end
    end
    return nothing
end

function initrun(code)
    visited_lines = Set{Int}()
    current_line = 1
    last_line = length(code)
    acc = 0
    while !in(current_line, visited_lines)
        t = code[current_line]
        push!(visited_lines, current_line)
        if t.instr == "nop"
            current_line += 1
        elseif t.instr == "acc"
            current_line += 1
            acc += t.value
        elseif t.instr == "jmp"
            current_line += t.value
        else
            println("Unknown instruction!")
        end
    end
    return visited_lines
end

end


text:
using Test, FinancialMonteCarlo, Statistics;
@show "Black Scholes Model"
S0 = 100.0;
K = 100.0;
r = [0.00, 0.02];
T = 1.0;
d = FinancialMonteCarlo.Curve([0.0, 0.02], T);
d2 = FinancialMonteCarlo.Curve([0.0, 0.01, 0.02], T);
D = 90.0;

Nsim = 10000;
Nstep = 30;
sigma = 0.2;
theta1 = 0.01;
k1 = 0.03;
mc = MonteCarloConfiguration(Nsim, Nstep);
mc1 = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.AntitheticMC());
toll = 0.8;

rfCurve = FinancialMonteCarlo.ImpliedZeroRate(r, T);
rfCurve2 = ZeroRate(r, T);

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
EUBin = BinaryEuropeanOption(T, K)
AMData = AmericanOption(T, K)
AmBin = BinaryEuropeanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianFloatingStrikeData = AsianFloatingStrikeOption(T)
AsianFixedStrikeData = AsianFixedStrikeOption(T, K)
Model = VarianceGammaProcess(sigma, theta1, k1, Underlying(S0, d));
Model2 = VarianceGammaProcess(sigma, theta1, k1, Underlying(S0, d2));

display(Model)
@show FwdPrice = pricer(Model, ZeroRate(r[end]), mc, FwdData);
@show FwdPrice = pricer(Model, rfCurve, mc, FwdData);
@show FwdPrice = pricer(Model, rfCurve2, mc, FwdData);
@show FwdPrice = pricer(Model2, rfCurve, mc, FwdData);
@show FwdPrice = pricer(Model, rfCurve, mc1, FwdData);
@show EuPrice = pricer(Model, rfCurve, mc, EUData);
@show EuBinPrice = pricer(Model, rfCurve, mc, EUBin);
@show AmPrice = pricer(Model, rfCurve, mc, AMData);
@show AmBinPrice = pricer(Model, rfCurve, mc, AmBin);
@show BarrierPrice = pricer(Model, rfCurve, mc, BarrierData);
@show AsianPrice1 = pricer(Model, rfCurve, mc, AsianFloatingStrikeData);
@show AsianPrice2 = pricer(Model, rfCurve, mc, AsianFixedStrikeData);