text:
# This script runs all the tests and create a table for comparison.

using RAFF
using Printf
using DelimitedFiles

fpath = "../test_problems/"

getType = begin

    d = Dict(
        "C" => "cubic",
        "L" => "linear",
        "LOG" => "logistic",
        "CI" => "circle",
        "E" => "expon"
    )

    (s) -> d[s]

end

function run_tests(;kwargs...)

    tfp = open("/tmp/table.csv", "w")

    @printf(tfp, """
| Name | Dim. | N Points | N Outl. | Found | Correct | Time (s) | Status | Solution |
| ---- | ---- | -------- | ------- | ----- | ------- | -------- | ------ | -------- |
""")

    for fname in readdir(fpath)

        occursin("sol", fname) && continue

        m = match(r"[^1-9]*", fname)

        (m == nothing) && continue

        mtype = getType(m.match)

        open(fpath * fname) do fp

            global N = parse(Int, readline(fp))

            global data = readdlm(fp)

        end

        n, model, = RAFF.model_list[mtype]

        np, = size(data)

        # Retrieve the true outliers
        outliers = findall(data[:, end] .== 1.0)

        # Run RAFF
        rsol, t, = @timed raff(model, data[:, 1:end - 1], n; kwargs...)

        # Save problem data

        name = match(r"[^.]*", fname).match

        @printf(tfp, "| %5s | %3d | %6d | %3d | ", name, n, np, length(outliers))

        nexact = 0

        for i in rsol.outliers

            (i in outliers) && (nexact += 1)

        end

        @printf(tfp, "%3d | %3d | %10.4f | %d | ",
                length(rsol.outliers), nexact, t, rsol.status)

        # Save best solution

        @printf(tfp, "[")

        for i = 1:n - 1

            @printf(tfp, "%10.3e, ", rsol.solution[i])

        end

        @printf(tfp, "%10.3e] |\n", rsol.solution[n])

    end

    close(tfp)

end


text:
"""
    Z = get_design_matrix(x, nested=[]; maxlev=0, cov=false)
    Z = get_design_matrix(x, y)
    Z = get_design_matrix(X; maxlev=[], cov=[])

Returns a design matrix accoding to class variables or covariates in a vector `x` or
matrix `X`. For a vector `x`, optional argument `dested` defines `x` as the nested
covariate within a class. Option `maxval` specifies the number of levels in this effect.
Omitting `maxval`, the function uses the maximum integer as the maximum level
in `x` or `X`. Option `cov` specifies the input is covariate and passes it through
the output. The matrix `X` doesn't support nested covariates.

If two vectors `x` and `y` are provided, the function returns the design matrix for the intercation between the two effects; the main effects are not included in the result.
Also, in this case, the function does not accept any additional options.

```juliadoctest
julia> x = [1,2,1,1,2]
julia> X = get_design_matrix(x)
5×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0
 1.0  0.0
 1.0  0.0
 0.0  1.0
```
"""
function get_design_matrix(x::Vector{<:Real}; nested::Vector{<:Integer}=Int[], maxlev=0, cov=false)
   n = length(x)

   # no change
   if length(nested)==0 && cov
      return x[:,:]
   end

   # maximum level
   if maxlev==0
      if length(nested)>0
         lv = maximum(nested)
      else
         lv = maximum(x)
      end
   else
      lv = maxlev
   end
   # fill the value
   Z = zeros(Float64,n,lv)
   if length(nested)>0
      # nested variable
      for i=1:n
         loc = Int(nested[i])
         if loc < 1 || loc > lv
            error("invalid category code: $(nested[i])")
         end
         Z[i,loc] = x[i]
      end
   else
      # non-nested case
      for i=1:n
         loc = Int(x[i])
         if loc < 1 || loc > lv
            error("invalid category code: $(x[i])")
         end
         Z[i,loc] = 1.0
      end
   end
   return Z
end

function get_design_matrix(X::Matrix{<:Real}; maxlev::Vector{<:Integer}=Int[], cov::Vector{Bool}=Bool[])
   # check
   (n,neff)=size(X)
   if length(maxlev)>0 && length(maxlev)!=neff
      error("argument maxlev size mismatch")
   end
   if length(cov)>0 && length(cov)!=neff
      error("argument cov size mismatch")
   end
   if length(cov)>0
      covar=cov
   else
      covar=zeros(Bool,neff)
   end
   # levels
   if length(maxlev)>0
      if minimum(maxlev)<0
         error("argument maxlev with 0 or negative value")
      end
      lv = maxlev
   else
      lv = zeros(Int,neff)
      for i=1:neff
         if covar[i]
            lv[i] = 1
         else
            lv[i] = maximum(X[:,i])
         end
      end
   end
   # build Y
   ncol = sum(lv)
   Z = zeros(Float64,n,ncol)
   collst = 0
   for i=1:neff
      colfst = collst + 1
      collst = colfst + lv[i] - 1
      Z[:,colfst:collst] .= get_design_matrix(X[:,i],maxlev=lv[i],cov=covar[i])
   end
   return Z
end

function get_design_matrix(x1::Vector{<:Real},x2::Vector{<:Real})
   X1 = get_design_matrix(x1)
   X2 = get_design_matrix(x2)
   return get_interaction_matrix(X1,X2)
end

"""
    X = get_interaction_matrix(X1, X2)

Returns a design matrix for the interaction between effects.
The design matrices for the two main effects (`X1` and `X2`) should be provided.
   Note that this function will not check if the supplied matrices are really the design matrices - it just creats all the combination of column products among the input matrices.
"""
function get_interaction_matrix(X1::Matrix{<:Real}, X2::Matrix{<:Real})
   if size(X1,1) != size(X2,1)
      throw(DimensionMismatch("input matrices X1 and X2"))
   end
   m = size(X1,1)
   n1 = size(X1,2)
   n2 = size(X2,2)
   Z = zeros(m,n1*n2)
   k = 0
   for i=1:n1
      for j=1:n2
         k = k + 1
         Z[:,k] = X1[:,i] .* X2[:,j]
      end
   end
   return Z
end

"""
     M = directsum(X...)

Returns the direct sum of supplied matrices.
"""
function directsum(X::Matrix...)
    num = length(X)
    mX = Int[ size(x, 1) for x in X ]
    nX = Int[ size(x, 2) for x in X ]
    m = sum(mX)
    n = sum(nX)
    Tv = promote_type(map(x->eltype(x), X)...)

    M =	  zeros(Tv,m,n)
    prev_i = 0
    prev_j = 0
    for i=1:num
       ifst = prev_i + 1
       ilst = ifst + mX[i] - 1
       jfst = prev_j + 1
       jlst = jfst + nX[i] - 1
       M[ifst:ilst,jfst:jlst] = X[i]
       prev_i =		      ilst
       prev_j =		      jlst
    end
    return M
end

ext: jl


text:
using AWSS3
using Test
using Dates
using AWSCore
using Retry
using HTTP
using FilePathsBase
using FilePathsBase.TestPaths
using LinearAlgebra  # for norm S3Path tests
using UUIDs: uuid4

AWSCore.set_debug_level(0)

@testset "AWSS3.jl" begin
    include("s3path.jl")
    include("awss3.jl")
end


text:
using .DynamicHMC: DynamicHMC

function DynamicHMC.GaussianKineticEnergy(M⁻¹::WoodburyPDMat)
    return DynamicHMC.GaussianKineticEnergy(M⁻¹, WoodburyLeftInvFactor(M⁻¹))
end

function DynamicHMC.kinetic_energy(
    κ::DynamicHMC.GaussianKineticEnergy{<:WoodburyPDMat}, p, q=nothing
)
    return PDMats.quad(κ.M⁻¹, p) / 2
end

# utility object so we can use DynamicHMC.GaussianKineticEnergy
struct WoodburyLeftInvFactor{T<:Real,TW<:PDMats.AbstractPDMat{T}} <: AbstractMatrix{T}
    A::TW
end

Base.size(L::WoodburyLeftInvFactor, dims::Int...) = size(L.A, dims...)

function LinearAlgebra.mul!(
    r::StridedVecOrMat{T}, L::WoodburyLeftInvFactor{T}, x::StridedVecOrMat{T}
) where {T}
    return invunwhiten!(r, L.A, x)
end

function Base.Matrix(L::WoodburyLeftInvFactor)
    W = L.A
    n, m = size(W.B)
    k = min(m, n)
    Lmat = zeros(n, n)
    Lmat[diagind(Lmat)] .= true
    @views ldiv!(W.UC, Lmat[1:k, 1:k])
    lmul!(W.Q, Lmat)
    ldiv!(W.UA, Lmat)
    return Lmat
end

function Base.AbstractMatrix{T}(L::WoodburyLeftInvFactor) where {T}
    return WoodburyLeftInvFactor(AbstractMatrix{T}(L.A))
end

function Base.getindex(L::WoodburyLeftInvFactor, i::Int, j::Int)
    n = size(L, 2)
    v = zeros(n)
    v[j] = 1
    mul!(v, L, v)
    return v[i]
end


text:
# For testing bug found on issue # 52
function solve_vaf_sdp(factory::JuMP.OptimizerFactory)
    model = Model(factory)
    @variable(model, x)
    @constraint(model, [4x x; x 4x] - ones(2, 2) in PSDCone())
    @objective(model, Min, x)
    optimize!(model)
    return JuMP.objective_value(model) # 0.6
end

@testset "JuMP dualize" begin
    @testset "direct_mode" begin
        JuMP_model = JuMP.direct_model(MOIU.MockOptimizer(MOIU.Model{Float64}()))
        err = ErrorException("Dualization does not support solvers in DIRECT mode")
        @test_throws ErrorException dualize(JuMP_model)
    end
    @testset "dualize JuMP models" begin
        # Only testing with a small LP
        for i in eachindex(primal_linear_factory)
            JuMP_model = JuMP.Model()
            MOI.copy_to(JuMP.backend(JuMP_model), lp1_test())
            dual_JuMP_model = dualize(JuMP_model)
            @test backend(dual_JuMP_model).state == MOIU.NO_OPTIMIZER
            dual_JuMP_model = dualize(JuMP_model, primal_linear_factory[i])
            @test backend(dual_JuMP_model).state == MOIU.EMPTY_OPTIMIZER
            @test MOI.get(backend(dual_JuMP_model), MOI.SolverName()) ==
                MOI.get(primal_linear_optimizer[i], MOI.SolverName())#"GLPK"
        end
    end
    @testset "set_dot on different sets" begin
        for i in eachindex(primal_power_cone_factory)
            primal_scs = solve_vaf_sdp(primal_power_cone_factory[i])
            dual_scs = solve_vaf_sdp(dual_power_cone_factory[i])
            @test isapprox(primal_scs, dual_scs; atol = 1e-3)
        end
        for i in eachindex(primal_conic_factory)
            primal_csdp = solve_vaf_sdp(primal_conic_factory[i])
            dual_csdp = solve_vaf_sdp(dual_conic_factory[i])
            @test isapprox(primal_csdp, dual_csdp; atol = 1e-6)
        end
    end
end


text:
module samplet4
using StaticArrays
function samplet4mesh()
    xyz = [
    0.0 0.0 0.0
    1.0 0.0 0.0
    2.0 0.0 0.0
    3.0 0.0 0.0
    0.0 2.0 0.0
    1.0 2.0 0.0
    2.0 2.0 0.0
    3.0 2.0 0.0
    0.0 4.0 0.0
    1.0 4.0 0.0
    2.0 4.0 0.0
    3.0 4.0 0.0
    0.0 6.0 0.0
    1.0 6.0 0.0
    2.0 6.0 0.0
    3.0 6.0 0.0
    0.0 8.0 0.0
    1.0 8.0 0.0
    2.0 8.0 0.0
    3.0 8.0 0.0
    0.0 0.0 2.5
    1.0 0.0 2.5
    2.0 0.0 2.5
    3.0 0.0 2.5
    0.0 2.0 2.5
    1.0 2.0 2.5
    2.0 2.0 2.5
    3.0 2.0 2.5
    0.0 4.0 2.5
    1.0 4.0 2.5
    2.0 4.0 2.5
    3.0 4.0 2.5
    0.0 6.0 2.5
    1.0 6.0 2.5
    2.0 6.0 2.5
    3.0 6.0 2.5
    0.0 8.0 2.5
    1.0 8.0 2.5
    2.0 8.0 2.5
    3.0 8.0 2.5
    0.0 0.0 5.0
    1.0 0.0 5.0
    2.0 0.0 5.0
    3.0 0.0 5.0
    0.0 2.0 5.0
    1.0 2.0 5.0
    2.0 2.0 5.0
    3.0 2.0 5.0
    0.0 4.0 5.0
    1.0 4.0 5.0
    2.0 4.0 5.0
    3.0 4.0 5.0
    0.0 6.0 5.0
    1.0 6.0 5.0
    2.0 6.0 5.0
    3.0 6.0 5.0
    0.0 8.0 5.0
    1.0 8.0 5.0
    2.0 8.0 5.0
    3.0 8.0 5.0
    ]
    c = [1 25 21 22
    6 5 2 26
    26 2 22 25
    5 26 25 2
    2 1 22 25
    5 25 1 2
    21 45 41 42
    26 25 22 46
    46 22 42 45
    25 46 45 22
    22 21 42 45
    25 45 21 22
    5 29 25 26
    10 9 6 30
    30 6 26 29
    9 30 29 6
    6 5 26 29
    9 29 5 6
    25 49 45 46
    30 29 26 50
    50 26 46 49
    29 50 49 26
    26 25 46 49
    29 49 25 26
    9 33 29 30
    14 13 10 34
    34 10 30 33
    13 34 33 10
    10 9 30 33
    13 33 9 10
    29 53 49 50
    34 33 30 54
    54 30 50 53
    33 54 53 30
    30 29 50 53
    33 53 29 30
    13 37 33 34
    18 17 14 38
    38 14 34 37
    17 38 37 14
    14 13 34 37
    17 37 13 14
    33 57 53 54
    38 37 34 58
    58 34 54 57
    37 58 57 34
    34 33 54 57
    37 57 33 34
    2 26 22 23
    7 6 3 27
    27 3 23 26
    6 27 26 3
    3 2 23 26
    6 26 2 3
    22 46 42 43
    27 26 23 47
    47 23 43 46
    26 47 46 23
    23 22 43 46
    26 46 22 23
    6 30 26 27
    11 10 7 31
    31 7 27 30
    10 31 30 7
    7 6 27 30
    10 30 6 7
    26 50 46 47
    31 30 27 51
    51 27 47 50
    30 51 50 27
    27 26 47 50
    30 50 26 27
    10 34 30 31
    15 14 11 35
    35 11 31 34
    14 35 34 11
    11 10 31 34
    14 34 10 11
    30 54 50 51
    35 34 31 55
    55 31 51 54
    34 55 54 31
    31 30 51 54
    34 54 30 31
    14 38 34 35
    19 18 15 39
    39 15 35 38
    18 39 38 15
    15 14 35 38
    18 38 14 15
    34 58 54 55
    39 38 35 59
    59 35 55 58
    38 59 58 35
    35 34 55 58
    38 58 34 35
    3 27 23 24
    8 7 4 28
    28 4 24 27
    7 28 27 4
    4 3 24 27
    7 27 3 4
    23 47 43 44
    28 27 24 48
    48 24 44 47
    27 48 47 24
    24 23 44 47
    27 47 23 24
    7 31 27 28
    12 11 8 32
    32 8 28 31
    11 32 31 8
    8 7 28 31
    11 31 7 8
    27 51 47 48
    32 31 28 52
    52 28 48 51
    31 52 51 28
    28 27 48 51
    31 51 27 28
    11 35 31 32
    16 15 12 36
    36 12 32 35
    15 36 35 12
    12 11 32 35
    15 35 11 12
    31 55 51 52
    36 35 32 56
    56 32 52 55
    35 56 55 32
    32 31 52 55
    35 55 31 32
    15 39 35 36
    20 19 16 40
    40 16 36 39
    19 40 39 16
    16 15 36 39
    19 39 15 16
    35 59 55 56
    40 39 36 60
    60 36 56 59
    39 60 59 36
    36 35 56 59
    39 59 35 36
    ]
    return xyz, c
end
end


text:
depwarn2(old::AbstractString, new::AbstractString, funcsym::Symbol) =
  Base.depwarn("$old is deprecated, use $new instead.", funcsym)

## Deprecated at 0.11.2

@deprecate logit(x::AbstractArray) logit.(x)
@deprecate invlogit(x::AbstractArray) invlogit.(x)


text:
import Base: .<, .>, .<=, .>=, .==, .!=
export .<, .>, .<=, .>=, .==, .!=, and, or, maxOf, minOf

import Base: .+, .-, .*, ./, .\, .%
export .+, .-, .*, ./, .\, .%

export addAssign!,subAssign!,divAssign!,mulAssign!,modAssign!

immutable Binary <: AFImpl
    le::Ptr{Void}
    lt::Ptr{Void}
    ge::Ptr{Void}
    gt::Ptr{Void}
    eq::Ptr{Void}
    neq::Ptr{Void}
    and::Ptr{Void}
    or::Ptr{Void}
    add::Ptr{Void}
    sub::Ptr{Void}
    mul::Ptr{Void}
    div::Ptr{Void}
    maxOf::Ptr{Void}
    minOf::Ptr{Void}
    mod::Ptr{Void}

    function Binary(ptr)
        new(
            Libdl.dlsym(ptr, :af_le),
            Libdl.dlsym(ptr, :af_lt),
            Libdl.dlsym(ptr, :af_ge),
            Libdl.dlsym(ptr, :af_gt),
            Libdl.dlsym(ptr, :af_eq),
            Libdl.dlsym(ptr, :af_neq),
            Libdl.dlsym(ptr, :af_and),
            Libdl.dlsym(ptr, :af_or),
            Libdl.dlsym(ptr, :af_add),
            Libdl.dlsym(ptr, :af_sub),
            Libdl.dlsym(ptr, :af_mul),
            Libdl.dlsym(ptr, :af_div),
            Libdl.dlsym(ptr, :af_maxof),
            Libdl.dlsym(ptr, :af_minof),
            Libdl.dlsym(ptr, :af_mod)
        )
    end
end

macro binOp(op, cFunc)
    quote
        function $(esc(op)){B}(lhs::AFArray{B}, rhs::AFArray{B})
            verifyAccess(lhs)
            verifyAccess(rhs)
            af = lhs.af
            result = af.results.ptr
            err = ccall(af.binary.$cFunc,
                Cint, (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Void}, Bool),
                result, lhs.ptr, rhs.ptr, af.batch)
            assertErr(err)
            AFArray{B}(af, result[])
        end

        function $(esc(op)){B}(lhs::AFArray{B}, rhsConst::Number)
            verifyAccess(lhs)
            af = lhs.af
            result = af.results.ptr
            rhs = constant(af, rhsConst, size(lhs)...)
            try
                err = ccall(af.binary.$cFunc,
                    Cint, (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Void}, Bool),
                    result, lhs.ptr, rhs.ptr, af.batch)
                assertErr(err)
                AFArray{B}(af, result[])
            finally
                release!(rhs)
            end
        end

        function $(esc(op)){B}(lhsConst::Number, rhs::AFArray{B})
            verifyAccess(rhs)
            af = rhs.af
            result = af.results.ptr
            lhs = constant(af, lhsConst, size(rhs)...)
            try
                err = ccall(af.binary.$cFunc,
                    Cint, (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Void}, Bool),
                    result, lhs.ptr, rhs.ptr, af.batch)
                assertErr(err)
                AFArray{B}(af, result[])
            finally
                release!(lhs)
            end
        end
    end
end

@binOp(.<, lt)
@binOp(.<=, le)
@binOp(.>, gt)
@binOp(.>=, ge)
@binOp(.==, eq)
@binOp(.!=, neq)
@binOp(and, and)
@binOp(or, or)

@binOp(.+, add)
@binOp(.-, sub)
@binOp(.*, mul)
@binOp(./, div)
@binOp(.%, mod)
@binOp(maxOf, maxOf)
@binOp(minOf, minOf)

function .\{B}(lhs::AFArray{B}, rhs::AFArray{B})
    rhs ./ lhs
end

function .\(lhs::AFArray, rhs::Number)
    rhs ./ lhs
end

function .\(lhs::Number, rhs::AFArray)
    rhs ./ lhs
end

macro binOpAssign(op, cFunc)
    quote
        function $(esc(op)){B}(lhs::AFArray{B}, rhs::AFArray{B})
            verifyAccess(lhs)
            verifyAccess(rhs)
            af = lhs.af
            result = af.results.ptr
            err = ccall(af.binary.$cFunc,
                Cint, (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Void}, Bool),
                result, lhs.ptr, rhs.ptr, af.batch)
            assertErr(err)
            updatePtr(lhs, result[])
        end

        function $(esc(op)){B}(lhs::AFArray{B}, rhsConst::Number)
            verifyAccess(lhs)
            af = lhs.af
            result = af.results.ptr
            rhs = constant(af, rhsConst, size(lhs)...)
            try
                err = ccall(af.binary.$cFunc,
                    Cint, (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Void}, Bool),
                    result, lhs.ptr, rhs.ptr, af.batch)
                assertErr(err)
                updatePtr(lhs, result[])
            finally
                release!(rhs)
            end
        end

        function $(esc(op)){B}(lhsConst::Number, rhs::AFArray{B})
            verifyAccess(rhs)
            af = rhs.af
            result = af.results.ptr
            lhs = constant(af, lhsConst, size(rhs)...)
            try
                err = ccall(af.binary.$cFunc,
                    Cint, (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Void}, Bool),
                    result, lhs.ptr, rhs.ptr, af.batch)
                assertErr(err)
                updatePtr(lhs, result[])
            finally
                release!(lhs)
            end
        end
    end
end

function updatePtr(arr::AFArray, ptr::Ptr{Void})
    if ptr != arr.ptr
        release!(arr.af, arr.ptr)
        arr.ptr = ptr
    end
end

@binOpAssign(addAssign!, add)
@binOpAssign(subAssign!, sub)
@binOpAssign(mulAssign!, mul)
@binOpAssign(divAssign!, div)
@binOpAssign(modAssign!, mod)

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
module KhepriUnity
using KhepriBase
using Sockets

# functions that need specialization
include(khepribase_interface_file())

include("Unity.jl")

function __init__()
   add_current_backend(unity)
end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 6fc66216-fa46-11eb-27bc-8fcb9143502c
begin
    import Pkg
    Pkg.activate(Base.current_project())

	using IntrinsicDecoherence
end

# ╔═╡ 46d9d011-ad41-49fb-aea0-8a55e93d97b5
begin
	using PyPlot
	using LaTeXStrings
	ion()
	rc("text", usetex = true)
end

# ╔═╡ 198df15f-890d-44f6-bb13-d9cde21c0e8a
begin
	const α = 4
	const r = 0.3
	const ω = 4
end;

# ╔═╡ be4e4005-6372-4d95-8c40-d9eb7d6a91f4
tlist = range(0, stop = 7, length = 500);

# ╔═╡ 1af1ee37-55e5-4f9f-b7bc-896a52ffdc8f
begin
	ctrs = [quads_a(α, 0, 0, 0.7, ω, 10, tlist[τ]) for τ in 1:length(tlist)];
	tsq1 = [quads_a(α, r, 0.7, 0.7, ω, 10, tlist[τ]) for τ in 1:length(tlist)];
	tsq2 = [quads_a(α, r, 1.4, 0.7, ω, 10, tlist[τ]) for τ in 1:length(tlist)];
	tsq3 = [quads_a(α, r, 2.1, 0.7, ω, 10, tlist[τ]) for τ in 1:length(tlist)];
end;

# ╔═╡ 0ce11282-0612-4236-91ea-8fd0a17b1fc6
begin
	figreim1 = figure(figsize = (10, 6))
	axreim1 = gca()

	#axreim1.plot(tlist, real(ctrs), "y--", lw = 2, alpha = 0.5)
	axreim1.plot(tlist, real(tsq1), "r", lw = 3, label = L"\theta = 0.7")
	axreim1.plot(tlist, real(tsq2), "g", lw = 3, label = L"\theta = 1.4")
	axreim1.plot(tlist, real(tsq3), "b", lw = 3, label = L"\theta = 2.1")

	#axreim1.set_ylim([-9, 9])

	#axreim1.text(1,4, L"r = 0.3", fontsize = 20)

	axreim1.legend(fontsize = 18)
	axreim1.set_ylabel(L"\textrm{Amplitude}", fontsize = 20)
	axreim1.set_xlabel(L"t\; [s]", fontsize = 18)
	setp(axreim1.get_xticklabels(), fontsize = 18)
	setp(axreim1.get_yticklabels(), fontsize = 18)
	axreim1.set_title(L"Expectation value of $\langle\hat{a}^{\dagger}+\hat{a}\rangle$", fontsize = 18)

	tight_layout()

	savefig("quads.pdf", dpi = 600)
end

# ╔═╡ e0492707-23ea-44bf-8530-4fe0df014f76
begin
	ctrS = [nums_a(α, 0, 0, 0.7, ω, 10, tlist[τ]) for τ in 1:length(tlist)];
	tSq1 = [nums_a(α, r, 0.7, 0.7, ω, 10, tlist[τ]) for τ in 1:length(tlist)];
	tSq2 = [nums_a(α, r, 1.4, 0.7, ω, 10, tlist[τ]) for τ in 1:length(tlist)];
	tSq3 = [nums_a(α, r, 2.1, 0.7, ω, 10, tlist[τ]) for τ in 1:length(tlist)];
end;

# ╔═╡ 0c1b4982-ad06-44d2-89b7-e4fab429521b
begin
	figreim2 = figure(figsize = (10, 6))
	axreim2 = gca()

	axreim2.plot(tlist, real(ctrS), "y--", lw = 2, alpha = 0.5)
	axreim2.plot(tlist, real(tSq1), "r", lw = 3, label = L"\theta = 0.7")
	axreim2.plot(tlist, real(tSq2), "g", lw = 3, label = L"\theta = 1.4")
	axreim2.plot(tlist, real(tSq3), "b", lw = 3, label = L"\theta = 2.1")

	#axreim1.set_ylim([-9, 9])

	#axreim2.text(5, 17.2, L"r = \theta = 0", fontsize = 16)
	axreim2.annotate(L"r = \theta = 0", xy=(5, 17.5), xytext=(5.5, 16), arrowprops = Dict("arrowstyle"=>"->"), fontsize = 14)

	axreim2.legend(fontsize = 18, loc = 1)
	axreim2.set_ylabel(L"\textrm{Number of Photons}", fontsize = 20)
	axreim2.set_xlabel(L"t\; [s]", fontsize = 18)
	setp(axreim2.get_xticklabels(), fontsize = 18)
	setp(axreim2.get_yticklabels(), fontsize = 18)
	axreim2.set_title(L"Expectation value of $\langle\hat{a}^{\dagger}\hat{a}\rangle$", fontsize = 18)

	tight_layout()

	savefig("nums.pdf", dpi = 600)
end

# ╔═╡ Cell order:
# ╠═6fc66216-fa46-11eb-27bc-8fcb9143502c
# ╠═46d9d011-ad41-49fb-aea0-8a55e93d97b5
# ╠═198df15f-890d-44f6-bb13-d9cde21c0e8a
# ╠═be4e4005-6372-4d95-8c40-d9eb7d6a91f4
# ╠═1af1ee37-55e5-4f9f-b7bc-896a52ffdc8f
# ╠═0ce11282-0612-4236-91ea-8fd0a17b1fc6
# ╠═e0492707-23ea-44bf-8530-4fe0df014f76
# ╠═0c1b4982-ad06-44d2-89b7-e4fab429521b

max_stars_repo_name: rurz/IntrinsicDecoherence

max_issues_repo_name: rurz/IntrinsicDecoherence

max_forks_repo_name: rurz/IntrinsicDecoherence


text:
using OrdinaryDiffEq
using HCubature

function ode_observations!(du::AbstractArray{Float64},
        u::AbstractArray{Float64}, fit::PhaseType, t::Float64)
    # dc = T * C + t * a
    a = fit.π' * exp(fit.T * t)
    du[:] = vec(fit.T * reshape(u, fit.p, fit.p) + fit.t * a)
end

function ode_censored!(du::AbstractArray{Float64}, u::AbstractArray{Float64},
        fit::PhaseType, t::Float64)
    # dc = T * C + 1 * a
    a = fit.π' * exp(fit.T * t)
    du[:] = vec(fit.T * reshape(u, fit.p, fit.p) + ones(fit.p) * a)
end

function c_integrand(u, fit, y)
    # Compute the two vector terms in the integrand
    first = fit.π' * exp(fit.T * u)
    second = exp(fit.T * (y-u)) * fit.t

    # Construct the matrix of integrands using outer product
    # then reshape it to a vector.
    C = second * first
    vec(C)
end

function d_integrand(u, fit, y)
    # Compute the two vector terms in the integrand
    first = fit.π' * exp(fit.T * u)
    second = exp(fit.T * (y-u)) * ones(fit.p)

    # Construct the matrix of integrands using outer product
    # then reshape it to a vector.
    D = second * first
    vec(D)
end

function e_step_observed_ode!(s::Sample, fit::PhaseType,
        Bs::AbstractArray{Float64}, Zs::AbstractArray{Float64},
        Ns::AbstractArray{Float64})
    # Setup initial conditions.
    p = fit.p
    u0 = zeros(p*p)

    # Run the ODE solver.
    prob = ODEProblem(ode_observations!, u0, (0.0, maximum(s.obs)), fit)
    sol = solve(prob, Tsit5())
    Id = Matrix{Float64}(I, p, p)

    for k = 1:length(s.obs)
        weight = s.obsweight[k]

        expTy = exp(fit.T * s.obs[k])
        a = transpose(fit.π' * expTy)
        b = expTy * fit.t

        u = sol(s.obs[k])
        C = reshape(u, p, p)

        if minimum(C) < 0
            (C,err) = hquadrature(x -> c_integrand(x, fit, s.obs[k]), 0,
            s.obs[k], rtol=1e-1, maxevals=500)
            C = reshape(C, p, p)
        end

        denom = fit.π' * b
        Bs[:] = Bs[:] + weight * (fit.π .* b) / denom
        Zs[:] = Zs[:] + weight * diag(C) / denom
        Ns[:,1:p] = Ns[:,1:p] + weight * (fit.T .* C' .* (1 .- Id)) / denom
        Ns[:,p+1] = Ns[:,end] + weight * (fit.t .* a) / denom
    end
end

function e_step_censored_ode!(s::Sample, fit::PhaseType,
        Bs::AbstractArray{Float64}, Zs::AbstractArray{Float64},
        Ns::AbstractArray{Float64})
    # Setup initial conditions.
    p = fit.p
    u0 = zeros(p*p)

    # Should tell the ODE solver to evaluate all right and interval censored
    # points.
    allcens = vcat(s.cens, vec(s.int))

    # Run the ODE solver.
    prob = ODEProblem(ode_censored!, u0, (0.0, maximum(allcens)), fit)
    sol = solve(prob, Tsit5())

    for k = 1:length(s.cens)
        weight = s.censweight[k]

        h = exp(fit.T * s.cens[k]) * ones(p)

        u = sol(s.cens[k])
        D = reshape(u, p, p)

        if minimum(D) < 0
            (D, err) = hquadrature(x -> d_integrand(x, fit, s.cens[k]), 0,
                s.cens[k], rtol=1e-1, maxevals=500)
            D = reshape(D, p, p)
        end

        denom = fit.π' * h
        Bs[:] = Bs[:] + weight * (fit.π .* h) / denom
        Zs[:] = Zs[:] + weight * diag(D) / denom
        Ns[:,1:p] = Ns[:,1:p] + weight * (fit.T .* transpose(D) .*
            (1 .- Matrix{Float64}(I, p, p))) / denom
    end

    a = 0; h = 0; g = 0; D = 0; expTRight = 0

    for k = 1:size(s.int, 1)
        weight = s.intweight[k]
        left, right = s.int[k,:]

        if k > 1 && left == s.int[k-1,2]
            a_left = a
            h_left = h
            g_left = g
            D_left = D
            expTLeft = expTRight
        else
            expTLeft = exp(fit.T * left)
            a_left = transpose(fit.π' * expTLeft)
            h_left = expTLeft * ones(p)
            g_left = try
                        transpose(fit.T) \ (a_left - fit.π)
                    catch
                        transpose(transpose(a_left - fit.π) * pinv(fit.T))
                    end

            u_left = sol(left)
            D_left = reshape(u_left, p, p)
            if minimum(D_left) < 0 && iscoxian(fit)
                minD = minimum(D_left)
                (D_left, err) = hquadrature(x -> d_integrand(x, fit, left), 0,
                    left, rtol=1e-1, maxevals=500)
                D_left = reshape(D_left, p, p)
            end
        end

        expTRight = exp(fit.T * right)
        a = transpose(fit.π' * expTRight)
        h = expTRight * ones(p)
        g = try
                transpose(fit.T) \ (a - fit.π)
            catch
                transpose(transpose(a - fit.π) * pinv(fit.T))
            end

        u_right = sol(right)
        D = reshape(u_right, p, p)
        if minimum(D) < 0 && iscoxian(fit)
            (D, err) = hquadrature(x -> d_integrand(x, fit, right), 0, right,
                rtol=1e-1, maxevals=500)
            D = reshape(D, p, p)
        end

        T_legal = fit.T .> 0

        denom = fit.π' * (expTLeft - expTRight) * ones(fit.p)

        Bs[:] = Bs[:] + weight * (fit.π .* (h_left - h)) / denom
        Zs[:] = Zs[:] + weight * (g - g_left + diag(D_left) - diag(D)) / denom
        Ns[:,1:p] = Ns[:,1:p] + weight * T_legal .* (fit.T .* ((g - g_left) .+
            transpose(D_left - D))) / denom
        Ns[:,p+1] = Ns[:,p+1] + weight * (fit.t .* (g - g_left)) / denom
    end
end


text:
using Distributions
using LinearAlgebra
using Random
using StatsFuns: softplus
using Roots # for inverse

################################################################################
#                            Planar and Radial Flows                           #
#             Ref: Variational Inference with Normalizing Flows,               #
#               D. Rezende, S. Mohamed(2015) arXiv:1505.05770                  #
################################################################################

###############
# PlanarLayer #
###############

mutable struct PlanarLayer{T1,T2} <: Bijector{1}
    w::T1
    u::T1
    b::T2
end

function get_u_hat(u, w)
    # To preserve invertibility
    return (
        u + (planar_flow_m(w' * u) - w' * u)[1]
        * w / (norm(w[:,1],2) ^ 2)
    )   # from A.1
end

function PlanarLayer(dims::Int, container=Array)
    w = container(randn(dims, 1))
    u = container(randn(dims, 1))
    b = container(randn(1))
    return PlanarLayer(w, u, b)
end

planar_flow_m(x) = -1 .+ softplus.(x)   # for planar flow from A.1
dtanh(x) = 1 .- (tanh.(x)) .^ 2         # for planar flow
ψ(z, w, b) = dtanh(w' * z .+ b) .* w    # for planar flow from eq(11)

# An internal version of transform that returns intermediate variables
function _transform(flow::PlanarLayer, z)
    u_hat = get_u_hat(flow.u, flow.w)
    transformed = z + u_hat * tanh.(flow.w' * z .+ flow.b) # from eq(10)
    return (transformed=transformed, u_hat=u_hat)
end

(b::PlanarLayer)(z) = _transform(b, z).transformed

function _forward(flow::PlanarLayer, z)
    transformed, u_hat = _transform(flow, z)
    # Compute log_det_jacobian
    psi = ψ(z, flow.w, flow.b)
    log_det_jacobian = log.(abs.(1.0 .+ psi' * u_hat))          # from eq(12)
    return (rv=transformed, logabsdetjac=vec(log_det_jacobian)) # from eq(10)
end

forward(flow::PlanarLayer, z::AbstractMatrix) = _forward(flow, z)

function forward(flow::PlanarLayer, z::AbstractVector{<: Real})
    res = _forward(flow, z)
    return (rv=vec(res.rv), logabsdetjac=res.logabsdetjac[1])
end

function (ib::Inversed{<: PlanarLayer})(y::AbstractMatrix{<: Real})
    flow = ib.orig
    u_hat = get_u_hat(flow.u, flow.w)
    # Define the objective functional; implemented with reference from A.1
    f(y) = alpha -> (flow.w' * y)[1] - alpha - (flow.w' * u_hat)[1] * tanh(alpha+flow.b[1])
    # Run solver
    alphas_ = [find_zero(f(y[:,i:i]), 0.0, Order16()) for i in 1:size(y, 2)]
    alphas = alphas_'
    z_para = (flow.w ./ norm(flow.w,2)) * alphas
    z_per = y - z_para - u_hat * tanh.(flow.w' * z_para .+ flow.b)

    return z_para + z_per
end

function (ib::Inversed{<: PlanarLayer})(y::AbstractVector{<: Real})
    return vec(ib(reshape(y, (length(y), 1))))
end

logabsdetjac(flow::PlanarLayer, x) = forward(flow, x).logabsdetjac
isclosedform(b::Inversed{<:PlanarLayer}) = false

###############
# RadialLayer #
###############

# FIXME: using `TrackedArray` for the parameters, we end up with
# nested tracked structures; don't want this.
mutable struct RadialLayer{T1,T2} <: Bijector{1}
    α_::T1
    β::T1
    z_0::T2
end

function RadialLayer(dims::Int, container=Array)
    α_ = container(randn(1))
    β = container(randn(1))
    z_0 = container(randn(dims, 1))
    return RadialLayer(α_, β, z_0)
end

h(α, r) = 1 ./ (α .+ r)     # for radial flow from eq(14)
dh(α, r) = - h(α, r) .^ 2   # for radial flow; derivative of h()

# An internal version of transform that returns intermediate variables
function _transform(flow::RadialLayer, z)
    α = softplus(flow.α_[1])            # from A.2
    β_hat = -α + softplus(flow.β[1])    # from A.2
    r = sqrt.(sum((z .- flow.z_0).^2; dims = 1))
    transformed = z + β_hat .* h(α, r) .* (z .- flow.z_0)   # from eq(14)
    return (transformed=transformed, α=α, β_hat=β_hat, r=r)
end

(b::RadialLayer)(z::AbstractMatrix{<:Real}) = _transform(b, z).transformed
(b::RadialLayer)(z::AbstractVector{<:Real}) = vec(_transform(b, z).transformed)

function _forward(flow::RadialLayer, z)
    transformed, α, β_hat, r = _transform(flow, z)
    # Compute log_det_jacobian
    d = size(flow.z_0, 1)
    h_ = h(α, r)
    log_det_jacobian = @. (
        (d - 1) * log(1.0 + β_hat * h_)
        + log(1.0 +  β_hat * h_ + β_hat * (- h_ ^ 2) * r)
    )   # from eq(14)
    return (rv=transformed, logabsdetjac=vec(log_det_jacobian))
end

forward(flow::RadialLayer, z) = _forward(flow, z)

function forward(flow::RadialLayer, z::AbstractVector{<:Real})
    res = _forward(flow, z)
    return (rv=vec(res.rv), logabsdetjac=res.logabsdetjac[1])
end

function (ib::Inversed{<:RadialLayer})(y)
    flow = ib.orig
    α = softplus(flow.α_[1])            # from A.2
    β_hat = - α + softplus(flow.β[1])   # from A.2
    # Define the objective functional
    f(y) = r -> norm(y - flow.z_0, 2) - r * (1 + β_hat / (α + r))   # from eq(26)
    # Run solver
    rs = [find_zero(f(y[:,i:i]), 0.0, Order16()) for i in 1:size(y, 2)]'    # from A.2
    z_hat = (y .- flow.z_0) ./ (rs .* (1 .+ β_hat ./ (α .+ rs)))            # from eq(25)
    z = flow.z_0 .+ rs .* z_hat # from A.2
    return z
end

function (ib::Inversed{<: RadialLayer})(y::AbstractVector{<:Real})
    return vec(ib(reshape(y, (length(y), 1))))
end

logabsdetjac(flow::RadialLayer, x) = forward(flow, x).logabsdetjac
isclosedform(b::Inversed{<:RadialLayer}) = false


text:
a{T<:Real}(x::T,t::T)=one(x)*one(t)
b{T<:Real}(x::T,t::T)=zero(x)*zero(t)
c{T<:Real}(x::T,t::T)=zero(x)*zero(t)

s{T<:Real}(t::T)=exp(t)
sp{T<:Real}(t::T)=exp(t)
u{T<:Real}(x::T,t::T)=exp(-pi^2*t)*cos(pi*x)
ux{T<:Real}(x::T,t::T)=pi*exp(-pi^2*t)*sin(pi*x)
ut{T<:Real}(x::T,t::T)=-pi^2*exp(-pi^2*t)*sin(pi*x)

γf{T<:Real}(x::T,t::T)=one(x)*one(t)

max_stars_repo_name: jgoldfar/isp-julia-mwe

max_issues_repo_name: jgoldfar/isp-julia-mwe

max_forks_repo_name: jgoldfar/isp-julia-mwe


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

########################################
# I/O
df_isolates     = CSV.read("/Users/glmarschmann/.julia/dev/DEBmicroTrait/final_manuscript/files/isolates2traits.csv", DataFrame, missingstring="N/A")
df_metabolites  = CSV.read("/Users/glmarschmann/.julia/dev/DEBmicroTrait/final_manuscript/files/exudation_properties.csv", DataFrame, missingstring="N/A")
########################################

########################################
# isolate traits
V_cell          = DEBmicroTrait.genome_size_to_cell_volume(convert(Array{Float64,1}, df_isolates.Genome_size))
Min_gen_time    = df_isolates.Min_gen_time
Gram_stain      = convert(Array{String,1}, df_isolates.gram_stain)
rrn_copies      = convert(Array{Float64,1}, df_isolates.rRNA_genes)

z_sugars        = reshape(convert(Array{Float64,1}, df_isolates.z_sugars./df_isolates.Genome_size*1e6),1,39)
z_organics      = reshape(convert(Array{Float64,1}, df_isolates.z_organic_acids./df_isolates.Genome_size*1e6),1,39)
z_aminos        = reshape(convert(Array{Float64,1}, df_isolates.z_amino_acids./df_isolates.Genome_size*1e6), 1,39)
z_fattys        = reshape(convert(Array{Float64,1}, df_isolates.z_fatty_acids./df_isolates.Genome_size*1e6),1,39)
z_nucleos       = reshape(convert(Array{Float64,1}, df_isolates.z_nucleotides./df_isolates.Genome_size*1e6),1,39)
z_auxins        = reshape(convert(Array{Float64,1}, df_isolates.z_auxins./df_isolates.Genome_size*1e6),1,39)
genome_distr    = vcat(z_sugars, z_organics, z_aminos, z_fattys, z_nucleos, z_auxins)

y_EM            = ones(size(V_cell,1))
########################################

########################################
# calc transporter density
ρ_ps            = zeros(size(df_metabolites.Name,1), size(V_cell,1))
y_DEs           = zeros(size(df_metabolites.Name,1), size(V_cell,1))

for j in 1:size(df_metabolites.Name,1)
    if df_metabolites.Ontology[j] == "Sugars"
        for i in 1:size(V_cell,1)
            #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
            find_ρ(x) = DEBmicroTrait.constrain_transporter_density_cost(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_metabolites.Formula[j])
            ρ_p = Roots.find_zero(find_ρ, 1.0)
            closure = genome_distr[:,i]./sum(genome_distr[:,i])
            ρ_ps[j,i] = ρ_p[1].*closure[1]
            y_DE = DEBmicroTrait.yield_transporter_density_cost(ρ_ps[j,i], [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_metabolites.Formula[j])
            y_DEs[j,i] = y_DE[1]
        end
    elseif df_metabolites.Ontology[j] == "Organic acids"
        for i in 1:size(V_cell,1)
            #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
            find_ρ(x) = DEBmicroTrait.constrain_transporter_density_cost(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_metabolites.Formula[j])
            ρ_p = Roots.find_zero(find_ρ, 1.0)
            closure = genome_distr[:,i]./sum(genome_distr[:,i])
            ρ_ps[j,i] = ρ_p[1].*closure[2]
            y_DE = DEBmicroTrait.yield_transporter_density_cost(ρ_ps[j,i], [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_metabolites.Formula[j])
            y_DEs[j,i] = y_DE[1]
        end
    elseif df_metabolites.Ontology[j] == "Amino acids"
        for i in 1:size(V_cell,1)
            #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
            find_ρ(x) = DEBmicroTrait.constrain_transporter_density_cost(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_metabolites.Formula[j])
            ρ_p = Roots.find_zero(find_ρ, 1.0)
            closure = genome_distr[:,i]./sum(genome_distr[:,i])
            ρ_ps[j,i] = ρ_p[1].*closure[3]
            y_DE = DEBmicroTrait.yield_transporter_density_cost(ρ_ps[j,i], [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_metabolites.Formula[j])
            y_DEs[j,i] = y_DE[1]
        end
    elseif df_metabolites.Ontology[j] == "Fatty acids"
        for i in 1:size(V_cell,1)
            #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
            find_ρ(x) = DEBmicroTrait.constrain_transporter_density_cost(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_metabolites.Formula[j])
            ρ_p = Roots.find_zero(find_ρ, 1.0)
            closure = genome_distr[:,i]./sum(genome_distr[:,i])
            ρ_ps[j,i] = ρ_p[1].*closure[4]
            y_DE = DEBmicroTrait.yield_transporter_density_cost(ρ_ps[j,i], [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_metabolites.Formula[j])
            y_DEs[j,i] = y_DE[1]
        end
    elseif df_metabolites.Ontology[j] == "Nucleotides"
        for i in 1:size(V_cell,1)
            #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
            find_ρ(x) = DEBmicroTrait.constrain_transporter_density_cost(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_metabolites.Formula[j])
            ρ_p = Roots.find_zero(find_ρ, 1.0)
            closure = genome_distr[:,i]./sum(genome_distr[:,i])
            ρ_ps[j,i] = ρ_p[1].*closure[5]
            y_DE = DEBmicroTrait.yield_transporter_density_cost(ρ_ps[j,i], [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_metabolites.Formula[j])
            y_DEs[j,i] = y_DE[1]
        end
    else
        for i in 1:size(V_cell,1)
            #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]],[rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
            find_ρ(x) = DEBmicroTrait.constrain_transporter_density_cost(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_metabolites.Formula[j])
            ρ_p = Roots.find_zero(find_ρ, 1.0)
            closure = genome_distr[:,i]./sum(genome_distr[:,i])
            ρ_ps[j,i] = ρ_p[1].*closure[6]
            y_DE = DEBmicroTrait.yield_transporter_density_cost(ρ_ps[j,i], [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_metabolites.Formula[j])
            y_DEs[j,i] = y_DE[1]
        end
    end
end

ρ_ps[ρ_ps.==0.0] .= 1e-12

N_C = zeros(size(df_metabolites.Name,1))
for i in 1:size(df_metabolites.Name,1)
    elementstring = df_metabolites.Formula[i]
    N_C[i] = DEBmicroTrait.extract_composition(elementstring)[1]
end

N_SB              = DEBmicroTrait.transporter_density_to_monomer_uptake_sites(V_cell, ρ_ps, Min_gen_time, Gram_stain)
Vmax              = @. 180.0*60^2*N_SB.*N_C

D_S               = DEBmicroTrait.aqueous_diffusivity(df_metabolites.Molecular_weight)
K_D               = DEBmicroTrait.specific_reference_affinity(V_cell, ρ_ps, D_S)

a_s               = Vmax./K_D
# I/O
df_assimilation = DataFrame()
df_assimilation.transporter_density = vec(ρ_ps)
df_assimilation.Vmax = vec(Vmax)
df_assimilation.KD = vec(K_D)
df_assimilation.affinity = vec(Vmax./K_D)

df_assimilation.ontology = repeat(df_metabolites.Ontology, size(V_cell,1))
response = Array{String}(undef, (size(df_metabolites.Name,1),size(V_cell,1)))
for i in 1:size(V_cell,1)
     response[:,i] .= df_isolates.Rhizosphere_response[i]
 end

df_assimilation.response = vec(response)
genome_size = Array{Int}(undef, (size(df_metabolites.Name,1),size(V_cell,1)))
for i in 1:size(V_cell,1)
     genome_size[:,i] .= df_isolates.Genome_size[i]
end
df_assimilation.genome_size = vec(genome_size)

cell_size = Array{Float64}(undef, (size(df_metabolites.Name,1),size(V_cell,1)))
for i in 1:size(V_cell,1)
     cell_size[:,i] .= V_cell[i]
end
df_assimilation.cell_size = vec(cell_size)
min_gt = Array{Float64}(undef, (size(df_metabolites.Name,1),size(V_cell,1)))
for i in 1:size(V_cell,1)
     min_gt[:,i] .= df_isolates.Min_gen_time[i]
end
df_assimilation.mingt = vec(min_gt)

SAV = Array{Float64}(undef, (size(df_metabolites.Name,1),size(V_cell,1)))
for i in 1:size(V_cell,1)
     SAV[:,i] .= DEBmicroTrait.surface_area_volume_ratio([V_cell[i]])
end
df_assimilation.sav = vec(SAV)

df_assimilation.yield = vec(abs.(y_DEs))

CSV.write("/Users/glmarschmann/.julia/dev/DEBmicroTrait/final_manuscript/files/isolates_assimilation.csv", df_assimilation)
save("/Users/glmarschmann/.julia/dev/DEBmicroTrait/final_manuscript/files/isolates_assimilation.jld", "rho", ρ_ps, "NSB", N_SB, "KD", K_D, "yEM", y_EM, "yDE", y_DEs, "NC", N_C)

# statistics
using HypothesisTests
#sugars
df_p_sugars = filter(x->(x.ontology.=="Sugars"&&x.response.=="positive") , df_assimilation)
df_n_sugars = filter(x->(x.ontology.=="Sugars"&&x.response.=="negative") , df_assimilation)
df_u_sugars = filter(x->(x.ontology.=="Sugars"&&x.response.=="undefined") , df_assimilation)
ρ_sugars    = KruskalWallisTest(df_p_sugars.transporter_density, df_n_sugars.transporter_density)
ρ_sugars_p  = median(df_p_sugars.transporter_density)
ρ_sugars_n  = median(df_n_sugars.transporter_density)
ρ_sugars_u  = median(df_u_sugars.transporter_density)
Vmax_sugars    = KruskalWallisTest(df_p_sugars.Vmax, df_n_sugars.Vmax)
Vmax_sugars_p  = median(df_p_sugars.Vmax)
Vmax_sugars_n  = median(df_n_sugars.Vmax)
Vmax_sugars_u  = median(df_u_sugars.Vmax)
K_sugars    = KruskalWallisTest(df_p_sugars.KD, df_n_sugars.KD)
K_sugars_p  = median(df_p_sugars.KD)
K_sugars_n  = median(df_n_sugars.KD)
K_sugars_u  = median(df_u_sugars.KD)
a_sugars    = KruskalWallisTest(df_p_sugars.Vmax./df_p_sugars.KD, df_n_sugars.Vmax./df_n_sugars.KD)
a_sugars_p  = median(df_p_sugars.Vmax./df_p_sugars.KD)
a_sugars_n  = median(df_n_sugars.Vmax./df_n_sugars.KD)

#organic acids
df_p_organics = filter(x->(x.ontology.=="Organic acids"&&x.response.=="positive") , df_assimilation)
df_n_organics = filter(x->(x.ontology.=="Organic acids"&&x.response.=="negative") , df_assimilation)
df_u_organics = filter(x->(x.ontology.=="Organic acids"&&x.response.=="undefined") , df_assimilation)
ρ_organics    = KruskalWallisTest(df_p_organics.transporter_density, df_n_organics.transporter_density)
ρ_organics_p  = median(df_p_organics.transporter_density)
ρ_organics_n  = median(df_n_organics.transporter_density)
ρ_organics_u  = median(df_u_organics.transporter_density)
Vmax_organics    = KruskalWallisTest(df_p_organics.Vmax, df_n_organics.Vmax)
Vmax_organics_p  = median(df_p_organics.Vmax)
Vmax_organics_n  = median(df_n_organics.Vmax)
Vmax_organics_u  = median(df_u_organics.Vmax)
K_organics    = KruskalWallisTest(df_p_organics.KD, df_n_organics.KD)
K_organics_p  = median(df_p_organics.KD)
K_organics_n  = median(df_n_organics.KD)
K_organics_u  = median(df_u_organics.KD)
a_organics    = KruskalWallisTest(df_p_organics.Vmax./df_p_organics.KD, df_n_organics.Vmax./df_n_organics.KD)
a_organics_p  = median(df_p_organics.Vmax./df_p_organics.KD)
a_organics_n  = median(df_n_organics.Vmax./df_n_organics.KD)

#amino acids
df_p_aminos = filter(x->(x.ontology.=="Amino acids"&&x.response.=="positive") , df_assimilation)
df_n_aminos = filter(x->(x.ontology.=="Amino acids"&&x.response.=="negative") , df_assimilation)
df_u_aminos = filter(x->(x.ontology.=="Amino acids"&&x.response.=="undefined") , df_assimilation)
ρ_aminos    = KruskalWallisTest(df_p_aminos.transporter_density, df_n_aminos.transporter_density,  df_u_aminos.transporter_density)
ρ_aminos_p  = median(df_p_aminos.transporter_density)
ρ_aminos_n  = median(df_n_aminos.transporter_density)
ρ_aminos_u  = median(df_u_aminos.transporter_density)
Vmax_aminos    = KruskalWallisTest(df_p_aminos.Vmax, df_n_aminos.Vmax)
Vmax_aminos_p  = median(df_p_aminos.Vmax)
Vmax_aminos_n  = median(df_n_aminos.Vmax)
Vmax_aminos_u  = median(df_u_aminos.Vmax)
K_aminos    = KruskalWallisTest(df_p_aminos.KD, df_n_aminos.KD)
K_aminos_p  = median(df_p_aminos.KD)
K_aminos_n  = median(df_n_aminos.KD)
K_aminos_u  = median(df_u_aminos.KD)
a_aminos    = KruskalWallisTest(df_p_aminos.Vmax./df_p_aminos.KD, df_n_aminos.Vmax./df_n_aminos.KD)
a_aminos_p  = median(df_p_aminos.Vmax./df_p_aminos.KD)
a_aminos_n  = median(df_n_aminos.Vmax./df_n_aminos.KD)

#fatty acids
df_p_fattys = filter(x->(x.ontology.=="Fatty acids"&&x.response.=="positive") , df_assimilation)
df_n_fattys = filter(x->(x.ontology.=="Fatty acids"&&x.response.=="negative") , df_assimilation)
df_u_fattys = filter(x->(x.ontology.=="Fatty acids"&&x.response.=="undefined") , df_assimilation)
ρ_fattys    = KruskalWallisTest(df_p_fattys.transporter_density, df_n_fattys.transporter_density, df_u_fattys.transporter_density)
ρ_fattys_p  = median(df_p_fattys.transporter_density)
ρ_fattys_n  = median(df_n_fattys.transporter_density)
ρ_fattys_u  = median(df_u_fattys.transporter_density)
Vmax_fattys    = KruskalWallisTest(df_p_fattys.Vmax, df_n_fattys.Vmax)
Vmax_fattys_p  = median(df_p_fattys.Vmax)
Vmax_fattys_n  = median(df_n_fattys.Vmax)
Vmax_fattys_u  = median(df_u_fattys.Vmax)
K_fattys    = KruskalWallisTest(df_p_fattys.KD, df_n_fattys.KD)
K_fattys_p  = median(df_p_fattys.KD)
K_fattys_n  = median(df_n_fattys.KD)
K_fattys_u  = median(df_u_fattys.KD)
a_fattys    = KruskalWallisTest(df_p_fattys.Vmax./df_p_fattys.KD, df_n_fattys.Vmax./df_n_fattys.KD)
a_fattys_p  = median(df_p_fattys.Vmax./df_p_fattys.KD)
a_fattys_n  = median(df_n_fattys.Vmax./df_n_fattys.KD)

# Nucleotides
df_p_nucleos = filter(x->(x.ontology.=="Nucleotides"&&x.response.=="positive") , df_assimilation)
df_n_nucleos = filter(x->(x.ontology.=="Nucleotides"&&x.response.=="negative") , df_assimilation)
df_u_nucleos = filter(x->(x.ontology.=="Nucleotides"&&x.response.=="undefined") , df_assimilation)
ρ_nucleos    = KruskalWallisTest(df_p_nucleos.transporter_density, df_n_nucleos.transporter_density, df_u_nucleos.transporter_density)
ρ_nucleos_p  = median(df_p_nucleos.transporter_density)
ρ_nucleos_n  = median(df_n_nucleos.transporter_density)
ρ_nucleos_u  = median(df_u_nucleos.transporter_density)
Vmax_nucleos    = KruskalWallisTest(df_p_nucleos.Vmax, df_n_nucleos.Vmax)
Vmax_nucleos_p  = median(df_p_nucleos.Vmax)
Vmax_nucleos_n  = median(df_n_nucleos.Vmax)
Vmax_nucleos_u  = median(df_u_nucleos.Vmax)
K_nucleos    = KruskalWallisTest(df_p_nucleos.KD, df_n_nucleos.KD)
K_nucleos_p  = median(df_p_nucleos.KD)
K_nucleos_n  = median(df_n_nucleos.KD)
K_nucleos_u  = median(df_u_nucleos.KD)
a_nucleos    = KruskalWallisTest(df_p_nucleos.Vmax./df_p_nucleos.KD, df_n_nucleos.Vmax./df_n_nucleos.KD)
a_nucleos_p  = median(df_p_nucleos.Vmax./df_p_nucleos.KD)
a_nucleos_n  = median(df_n_nucleos.Vmax./df_n_nucleos.KD)

#Auxins
df_p_auxins = filter(x->(x.ontology.=="Auxins"&&x.response.=="positive") , df_assimilation)
df_n_auxins = filter(x->(x.ontology.=="Auxins"&&x.response.=="negative") , df_assimilation)
df_u_auxins = filter(x->(x.ontology.=="Auxins"&&x.response.=="undefined") , df_assimilation)
ρ_auxins    = KruskalWallisTest(df_p_auxins.transporter_density, df_n_auxins.transporter_density)
ρ_auxins_p  = median(df_p_auxins.transporter_density)
ρ_auxins_n  = median(df_n_auxins.transporter_density)
ρ_auxins_u  = median(df_u_auxins.transporter_density)
Vmax_auxins    = KruskalWallisTest(df_p_auxins.Vmax, df_n_auxins.Vmax)
Vmax_auxins_p  = median(df_p_auxins.Vmax)
Vmax_auxins_n  = median(df_n_auxins.Vmax)
Vmax_auxins_u  = median(df_u_auxins.Vmax)
K_auxins    = KruskalWallisTest(df_p_auxins.KD, df_n_auxins.KD)
K_auxins_p  = median(df_p_auxins.KD)
K_auxins_n  = median(df_n_auxins.KD)
K_auxins_u  = median(df_u_auxins.KD)
a_auxins    = KruskalWallisTest(df_p_auxins.Vmax./df_p_auxins.KD, df_n_auxins.Vmax./df_n_auxins.KD)
a_auxins_p  = median(df_p_auxins.Vmax./df_p_auxins.KD)
a_auxins_n  = median(df_n_auxins.Vmax./df_n_auxins.KD)

# yields
df_y_sugars = filter(x->(x.ontology.=="Sugars"), df_assimilation)
y_sugars    = @. 1/df_y_sugars.yield
y_sugar_med = median(y_sugars)
df_y_organics = filter(x->(x.ontology.=="Organic acids"), df_assimilation)
y_organics     =  @. 1/df_y_organics.yield
y_organics_med = median(y_organics)
df_y_aminos = filter(x->(x.ontology.=="Amino acids"), df_assimilation)
y_aminos     =  @. 1/df_y_aminos.yield
y_aminos_med = median(y_aminos)
df_y_fattys = filter(x->(x.ontology.=="Fatty acids"), df_assimilation)
y_fattys     =  @. 1/df_y_fattys.yield
y_fattys_med = median(y_fattys)
df_y_nucleos = filter(x->(x.ontology.=="Nucleotides"), df_assimilation)
y_nucleos     =  @. 1/df_y_nucleos.yield
y_nucloes_med = median(y_nucleos)
df_y_auxins = filter(x->(x.ontology.=="Auxins"), df_assimilation)
y_auxins     =  @. 1/df_y_auxins.yield
y_auxins_med = median(y_auxins)

kw_sugars_organics = KruskalWallisTest(y_sugars, y_organics)
kw_sugars_aminos = KruskalWallisTest(y_sugars, y_aminos)
kw_sugars_fattys = KruskalWallisTest(y_sugars, y_fattys)
kw_sugars_auxins = KruskalWallisTest(y_sugars, y_auxins)
kw_sugars_nucleos = KruskalWallisTest(y_sugars, y_nucleos)

kw_organics_aminos =  KruskalWallisTest(y_organics, y_aminos)
kw_organics_fattys =  KruskalWallisTest(y_organics, y_fattys)
kw_organics_auxins =  KruskalWallisTest(y_organics, y_auxins)
kw_organics_nucleos =  KruskalWallisTest(y_organics, y_nucleos)

kw_aminos_fattys =  KruskalWallisTest(y_aminos, y_fattys)
kw_aminos_auxins =  KruskalWallisTest(y_aminos, y_auxins)
kw_aminos_nucleos =  KruskalWallisTest(y_aminos, y_nucleos)

kw_fattys_auxins =  KruskalWallisTest(y_auxins, y_fattys)
kw_fattys_nucleos =  KruskalWallisTest(y_nucleos, y_fattys)

kw_nucleos_auxins = KruskalWallisTest(y_nucleos, y_auxins)

################################################################################
# Aromatic organic acids

df_aromatics = filter(x->(x.Name.=="nicotinic acid" || x.Name.=="shikimic acid" ||
                          x.Name.=="shikimic acid" || x.Name.=="salicylic acid" ||
                          x.Name.=="thymidine" || x.Name.=="threonic acid" ||
                          x.Name.=="indole-3-acetic acid" || x.Name.=="deoxyguanosine" ||
                          x.Name.=="cytidine" || x.Name.=="choline" ||
                          x.Name.=="3-Dehydroshikimic acid" || x.Name.=="trans-cinnamic acid" ||
                          x.Name.=="aminobutyric acid" || x.Name.=="L-pyroglutamic acid") , df_metabolites)

########################################
# isolate traits
V_cell          = DEBmicroTrait.genome_size_to_cell_volume(convert(Array{Float64,1}, df_isolates.Genome_size))
Min_gen_time    = df_isolates.Min_gen_time
Gram_stain      = convert(Array{String,1}, df_isolates.gram_stain)
rrn_copies      = convert(Array{Float64,1}, df_isolates.rRNA_genes)

z_sugars        = reshape(convert(Array{Float64,1}, df_isolates.z_sugars./df_isolates.Genome_size*1e6),1,39)
z_organics      = reshape(convert(Array{Float64,1}, df_isolates.z_organic_acids./df_isolates.Genome_size*1e6),1,39)
z_aminos        = reshape(convert(Array{Float64,1}, df_isolates.z_amino_acids./df_isolates.Genome_size*1e6), 1,39)
z_fattys        = reshape(convert(Array{Float64,1}, df_isolates.z_fatty_acids./df_isolates.Genome_size*1e6),1,39)
z_nucleos       = reshape(convert(Array{Float64,1}, df_isolates.z_nucleotides./df_isolates.Genome_size*1e6),1,39)
z_auxins        = reshape(convert(Array{Float64,1}, df_isolates.z_auxins./df_isolates.Genome_size*1e6),1,39)
genome_distr    = vcat(z_sugars, z_organics, z_aminos, z_fattys, z_nucleos, z_auxins)

y_EM            = ones(size(V_cell,1))
########################################

########################################
# calc transporter density
ρ_ps            = zeros(size(df_aromatics.Name,1), size(V_cell,1))
y_DEs           = zeros(size(df_aromatics.Name,1), size(V_cell,1))

for j in 1:size(df_aromatics.Name,1)
  if df_aromatics.Ontology[j] == "Sugars"
      for i in 1:size(V_cell,1)
          #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
          find_ρ(x) = DEBmicroTrait.constrain_transporter_density_cost(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_aromatics.Formula[j])
          ρ_p = Roots.find_zero(find_ρ, 1.0)
          closure = genome_distr[:,i]./sum(genome_distr[:,i])
          ρ_ps[j,i] = ρ_p[1].*closure[1]
          y_DE = DEBmicroTrait.yield_transporter_density_cost(ρ_ps[j,i], [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_aromatics.Formula[j])
          y_DEs[j,i] = y_DE[1]
      end
  elseif df_aromatics.Ontology[j] == "Organic acids"
      for i in 1:size(V_cell,1)
          #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
          find_ρ(x) = DEBmicroTrait.constrain_transporter_density_cost(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_aromatics.Formula[j])
          ρ_p = Roots.find_zero(find_ρ, 1.0)
          closure = genome_distr[:,i]./sum(genome_distr[:,i])
          ρ_ps[j,i] = ρ_p[1].*closure[2]
          y_DE = DEBmicroTrait.yield_transporter_density_cost(ρ_ps[j,i], [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_aromatics.Formula[j])
          y_DEs[j,i] = y_DE[1]
      end
  elseif df_aromatics.Ontology[j] == "Amino acids"
      for i in 1:size(V_cell,1)
          #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
          find_ρ(x) = DEBmicroTrait.constrain_transporter_density_cost(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_aromatics.Formula[j])
          ρ_p = Roots.find_zero(find_ρ, 1.0)
          closure = genome_distr[:,i]./sum(genome_distr[:,i])
          ρ_ps[j,i] = ρ_p[1].*closure[3]
          y_DE = DEBmicroTrait.yield_transporter_density_cost(ρ_ps[j,i], [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_aromatics.Formula[j])
          y_DEs[j,i] = y_DE[1]
      end
  elseif df_aromatics.Ontology[j] == "Fatty acids"
      for i in 1:size(V_cell,1)
          #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
          find_ρ(x) = DEBmicroTrait.constrain_transporter_density_cost(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_aromatics.Formula[j])
          ρ_p = Roots.find_zero(find_ρ, 1.0)
          closure = genome_distr[:,i]./sum(genome_distr[:,i])
          ρ_ps[j,i] = ρ_p[1].*closure[4]
          y_DE = DEBmicroTrait.yield_transporter_density_cost(ρ_ps[j,i], [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_aromatics.Formula[j])
          y_DEs[j,i] = y_DE[1]
      end
  elseif df_aromatics.Ontology[j] == "Nucleotides"
      for i in 1:size(V_cell,1)
          #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
          find_ρ(x) = DEBmicroTrait.constrain_transporter_density_cost(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_aromatics.Formula[j])
          ρ_p = Roots.find_zero(find_ρ, 1.0)
          closure = genome_distr[:,i]./sum(genome_distr[:,i])
          ρ_ps[j,i] = ρ_p[1].*closure[5]
          y_DE = DEBmicroTrait.yield_transporter_density_cost(ρ_ps[j,i], [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_aromatics.Formula[j])
          y_DEs[j,i] = y_DE[1]
      end
  else
      for i in 1:size(V_cell,1)
          #find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]],[rrn_copies[i]], y_DE[j], N_C[j], [y_EM[i]])
          find_ρ(x) = DEBmicroTrait.constrain_transporter_density_cost(x, [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_aromatics.Formula[j])
          ρ_p = Roots.find_zero(find_ρ, 1.0)
          closure = genome_distr[:,i]./sum(genome_distr[:,i])
          ρ_ps[j,i] = ρ_p[1].*closure[6]
          y_DE = DEBmicroTrait.yield_transporter_density_cost(ρ_ps[j,i], [V_cell[i]], [Min_gen_time[i]], [Gram_stain[i]], [rrn_copies[i]], [y_EM[i]], df_aromatics.Formula[j])
          y_DEs[j,i] = y_DE[1]
      end
  end
end

ρ_ps[ρ_ps.==0.0] .= 1e-12

N_C = zeros(size(df_aromatics.Name,1))
for i in 1:size(df_aromatics.Name,1)
  elementstring = df_aromatics.Formula[i]
  N_C[i] = DEBmicroTrait.extract_composition(elementstring)[1]
end

N_SB              = DEBmicroTrait.transporter_density_to_monomer_uptake_sites(V_cell, ρ_ps, Min_gen_time, Gram_stain)
Vmax              = @. 180.0*60^2*N_SB.*N_C

D_S               = DEBmicroTrait.aqueous_diffusivity(df_aromatics.Molecular_weight)
K_D               = DEBmicroTrait.specific_reference_affinity(V_cell, ρ_ps, D_S)

df_aromatics_out = DataFrame()
df_aromatics_out.Vmax = vec(Vmax)
df_aromatics_out.KD = vec(K_D)
df_aromatics_out.affinity = vec(Vmax./K_D)

df_aromatics_out.ontology = repeat(df_aromatics.Ontology, size(V_cell,1))
response = Array{String}(undef, (size(df_aromatics.Name,1),size(V_cell,1)))
for i in 1:size(V_cell,1)
     response[:,i] .= df_isolates.Rhizosphere_response[i]
 end

df_aromatics_out.response = vec(response)

df_aromatics_out.yield = vec(abs.(y_DEs))
df_aromatics_out.name = repeat(df_aromatics.Name, size(V_cell,1))

df_aromatics_out.transporter_density = vec(ρ_ps)

CSV.write("/Users/glmarschmann/.julia/dev/DEBmicroTrait/final_manuscript/files/isolates_assimilation_aromatics.csv", df_aromatics_out)

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
function evaluation_performed(obj::KalmanObserver, arg0::KalmanEstimation)
    return jcall(obj, "evaluationPerformed", void, (KalmanEstimation,), arg0)
end


text:
@testset "arraymath" begin
    @testset "inv(::Matrix{$T})" for T in (Float64, ComplexF64)
        N = 3
        B = generate_well_conditioned_matrix(T, N)
        frule_test(inv, (B, randn(T, N, N)))
        rrule_test(inv, randn(T, N, N), (B, randn(T, N, N)))
    end

    @testset "*: $T" for T in (Float64, ComplexF64)
        ⋆(a) = round.(5*randn(T, a))  # Helper to generate nice random values
        ⋆(a, b) = ⋆((a, b))  # matrix
        ⋆() = only(⋆(()))  # scalar

        ⋆₂(a) = (⋆(a), ⋆(a)) # Helper to generate random matrix and its cotangent
        ⋆₂(a, b) = ⋆₂((a, b))  #matrix
        ⋆₂() = ⋆₂(())  # scalar

        @testset "Scalar-Array $dims" for dims in ((3,), (5,4), (2, 3, 4, 5))
            rrule_test(*, ⋆(dims), ⋆₂(), ⋆₂(dims))
            rrule_test(*, ⋆(dims), ⋆₂(dims), ⋆₂())
        end

        @testset "AbstractMatrix-AbstractVector n=$n, m=$m" for n in (2, 3), m in (4, 5)
            @testset "Array" begin
                rrule_test(*, ⋆(n), n ⋆₂ m, ⋆₂(m))
            end
        end

        @testset "AbstractVector-AbstractMatrix n=$n, m=$m" for n in (2, 3), m in (4, 5)
            @testset "Array" begin
                rrule_test(*, n ⋆ m, ⋆₂(n), 1 ⋆₂ m)
            end
        end

        @testset "AbstractMatrix-AbstractMatrix" begin
            @testset "Matrix * Matrix n=$n, m=$m, p=$p" for n in (2, 5), m in (2, 4), p in (2, 3)
                @testset "Array" begin
                    rrule_test(*, n⋆p, (n⋆₂m), (m⋆₂p))
                end

                @testset "SubArray - $indexname" for (indexname, m_index) in (
                    ("fast", :), ("slow", Ref(m:-1:1))
                )
                    rrule_test(*, n⋆p, view.(n⋆₂m, :, m_index), view.(m⋆₂p, m_index, :))
                    rrule_test(*, n⋆p, n⋆₂m, view.(m⋆₂p, m_index, :))
                    rrule_test(*, n⋆p, view.(n⋆₂m, :, m_index), m⋆₂p)
                end

                @testset "Adjoints and Transposes" begin
                    rrule_test(*, n⋆p, Transpose.(m⋆₂n), Transpose.(p⋆₂m))
                    rrule_test(*, n⋆p, Adjoint.(m⋆₂n), Adjoint.(p⋆₂m))

                    rrule_test(*, n⋆p, Transpose.(m⋆₂n), (m⋆₂p))
                    rrule_test(*, n⋆p, Adjoint.(m⋆₂n), (m⋆₂p))

                    rrule_test(*, n⋆p, (n⋆₂m), Transpose.(p⋆₂m))
                    rrule_test(*, n⋆p, (n⋆₂m), Adjoint.(p⋆₂m))
                end
            end
        end

        @testset "Covector * Vector n=$n" for n in (3, 5)
            @testset "$f" for f in (adjoint, transpose)
                # This should be same as dot product and give a scalar
                rrule_test(*, ⋆(), f.(⋆₂(n)), ⋆₂(n))
            end
        end
    end

    @testset "$f" for f in (/, \)
        @testset "Matrix" begin
            for n in 3:5, m in 3:5
                A = randn(m, n)
                B = randn(m, n)
                Ȳ = randn(size(f(A, B)))
                rrule_test(f, Ȳ, (A, randn(m, n)), (B, randn(m, n)))
            end
        end
        @testset "Vector" begin
            x = randn(10)
            y = randn(10)
            ȳ = randn(size(f(x, y))...)
            rrule_test(f, ȳ, (x, randn(10)), (y, randn(10)))
        end
        if f == (\)
            @testset "Matrix $f Vector" begin
                X = randn(10, 4)
                y = randn(10)
                ȳ = randn(size(f(X, y))...)
                rrule_test(f, ȳ, (X, randn(size(X))), (y, randn(10)))
            end
            @testset "Vector $f Matrix" begin
                x = randn(10)
                Y = randn(10, 4)
                ȳ = randn(size(f(x, Y))...)
                rrule_test(f, ȳ, (x, randn(size(x))), (Y, randn(size(Y))))
            end
        end
    end
    @testset "/ and \\ Scalar-AbstractArray" begin
        A = randn(3, 4, 5)
        Ā = randn(3, 4, 5)
        Ȳ = randn(3, 4, 5)
        rrule_test(/, Ȳ, (A, Ā), (7.2, 2.3))
        rrule_test(\, Ȳ, (7.2, 2.3), (A, Ā))
    end

    @testset "negation" begin
        A = randn(4, 4)
        Ā = randn(4, 4)
        Ȳ = randn(4, 4)
        rrule_test(-, Ȳ, (A, Ā))
        rrule_test(-, Diagonal(Ȳ), (Diagonal(A), Diagonal(Ā)))
    end
end


text:
# Uses new formulation from limbdark paper.
# Plot of the flux of a linearly limb-darkened star with u_1=1.
# Normalized so that uneclipsed star has a flux of unity.
using Limbdark
import Limbdark: s2
#include("../../../src/s2.jl")
#include("../../../src/define_constants.jl")
include("../../../test/loglinspace.jl")

# Run some tests:

using PyPlot
nb = 401; nr = 401
b=linearspace(0.0,2,nr)
r=linearspace(0.0,2,nr)
fgrid = zeros(Float64,nr,nb)
# Normalization constant (which is uneclipsed s_1):
norm = 3/(2pi)
for ib=1:nb
#  if mod(ib,10) == 0
#    println("Finished: ",ib/nr*100,"%")
#  end
  for ir=1:nr
    fgrid[ir,ib]=norm*s2(r[ir],b[ib])
  end
end

clf()
img=imshow(fgrid,interpolation="nearest",origin="lower",extent=[0,2,0,2],
           cmap="plasma")
plot([0,1],[1,2],"k-")
plot([1,2],[0,1],"k-")
plot([0,0.38],[0,0.38],"k-")
plot([0.5,2],[0.5,2],"k-")
plot([0,1],[1,0],"k-")

annotate(L"$r = 1 + b$", xy=(0.5, 1.5), xytext=(0, 0), xycoords="data",
         textcoords="offset points", ha="center", va="bottom",
         fontsize=12, color="w", rotation=45)

annotate(L"$r = b$", xy=(1, 1), xytext=(0, 0), xycoords="data",
         textcoords="offset points", ha="center", va="bottom",
         fontsize=12, color="k", rotation=45)

annotate(L"$r = b - 1$", xy=(1.5, 0.5), xytext=(0, 0), xycoords="data",
         textcoords="offset points", ha="center", va="bottom",
         fontsize=12, color="k", rotation=45)

annotate(L"$r = 1 - b$", xy=(0.5, 0.5), xytext=(-9, -9), xycoords="data",
         textcoords="offset points", ha="center", va="center",
         fontsize=12, color="k", rotation=-45)

xlabel("Impact parameter")
ylabel("Radius ratio")
colorbar(img)
savefig("transit_linear.pdf", bbox_inches="tight")

max_issues_repo_name: rodluger/limbdark


text:
import ..GPU

"""
    postprocess(yolomat::Array{Float32}, conf_thresh::Float32, iou_thresh::Float32)

Post processing function.
Confidence score threshold to select correct predictions. Recommended : 0.3
IoU threshold to remove unnecessary predictions: Recommended:0.3
"""
function postprocess(yolomat::Array{Float32},settings::Settings; conf_thresh::T = 0.3, iou_thresh::T = 0.3) where {T<:AbstractFloat}
    im_w = settings.image_shape[1]
    im_h = settings.image_shape[2]
    num_images = size(yolomat,4)
    all_detections = Vector{Vector{PredictLabel}}(undef, num_images)
    RATE = 32
    for i in 1:num_images
        all_detections[i] = PredictLabel[]
        for cy in 1:13
            for cx in 1:13
                @views for b in 1:5
                    channel = (b-1)*(settings.num_classes + 5)
                    tx = yolomat[cy,cx,channel+1,i]
                    ty = yolomat[cy,cx,channel+2,i]
                    tw = yolomat[cy,cx,channel+3,i]
                    th = yolomat[cy,cx,channel+4,i]
                    tc = yolomat[cy,cx,channel+5,i]
                    x = (sigmoid(tx) + cx-1) * RATE
                    y = (sigmoid(ty) + cy-1) * RATE
                    w = exp(tw) * settings.anchors[b][1] * RATE
                    h = exp(th) * settings.anchors[b][2] * RATE
                    conf = sigmoid(tc)
                    classScores = yolomat[cy,cx,channel+6:channel+25,i]
                    classScores = softmax(classScores)
                    classNo = argmax(classScores)
                    bestScore = classScores[classNo]
                    classConfidenceScore = conf * bestScore
                    if classConfidenceScore > conf_thresh
                        bbox = BBOX(
                            x = max(0.0,x-w/2)/im_w,
                            y = max(0.0,y-h/2)/im_h,
                            w = min(w,im_w)/im_w,
                            h = min(h,im_h)/im_h)
                        p = PredictLabel(bbox = bbox, class = classNo, conf = classConfidenceScore)
                        push!(all_detections[i], p)
                    end
                end
            end
        end
    end
    nonMaxSupression!.(all_detections, iou_thresh)
    if length(all_detections) == 1
        return all_detections[1]
    else
        return all_detections
    end
end
function postprocess(yolomat::KnetArray{Float32}, settings::Settings; conf_thresh::T = 0.3, iou_thresh::T = 0.3) where {T<:AbstractFloat}
    postprocess(Array{Float32}(yolomat), settings; conf_thresh = conf_thresh, iou_thresh = iou_thresh)
end
"""
    nonMaxSupression!(detections::Vector{PredictLabel}, iou_thresh::Float32)

Removes the predictions overlapping.
"""
function nonMaxSupression!(detections::Vector{PredictLabel}, iou_thresh::T) where {T<:AbstractFloat}
    sort!(detections, by = x -> x.conf, rev = true)
    for i = 1:length(detections)
        k = i + 1
        @views while k <= length(detections)
            iou = ioumatch(detections[i].bbox, detections[k].bbox)
            if iou > iou_thresh && (detections[i].class == detections[k].class)
                deleteat!(detections, k)
                k -= 1
            end
            k += 1
        end
    end
end

"""
    ioumatch(bbox1::BBOX, bbox2::BBOX)

Calculates IoU score (overlapping rate)
"""
function ioumatch(bbox1::BBOX, bbox2::BBOX)
    r1 = bbox1.x + bbox1.w
    l1 = bbox1.x
    t1 = bbox1.y
    b1 = bbox1.y + bbox1.h
    r2 = bbox2.x + bbox2.w
    l2 = bbox2.x
    t2 = bbox2.y
    b2 = bbox2.y + bbox2.h
    a = min(r1, r2)
    b = max(t1, t2)
    c = max(l1, l2)
    d = min(b1, b2)
    intersec = (d - b) * (a - c)
    return intersec / (bbox1.w * bbox1.h + bbox2.w * bbox2.h - intersec)
end

imgslice_to_rgbimg(imgslice) = collect(colorview(RGB, permutedims(imgslice,[3,1,2])))

## TODO Convert the remainer

# #Calculates accuracy for Voc Dataset
# acc(model,data,conf_thresh,iou_thresh,iou,predictions) =(acc!(model,args,conf_thresh,iou_thresh,iou,predictions) for args in data)
#
# function acc!(model,args,conf_thresh,iou_thresh,iou,predictions)
#     yolomat = model(args[1])
#     yolomat = postprocess(yolomat,conf_thresh,iou_thresh)
#     check = zeros(length(args[2][1])-2)
#     sort!(yolomat,by = x-> x[6],rev=true)
#     for k in 1:length(yolomat)
#         tp,loc = istrue(yolomat[k],args[2][1][3:length(args[2][1])],check,iou)
#         push!(predictions[numsdic[yolomat[k][5]]],(tp,yolomat[k][6]))
#         if tp
#             check[loc] = 1
#         end
#     end
# end
# #conf_thresh => confidence score threshold. 0.0 for calculating accuracy
# #iou_thresh => intersection over union threshold. if 2 images overlap more than this threshold, one of them is removed
# #iou => intersection over union. True positive threshold
# function accuracy(model,data,conf_thresh,iou_thresh,iou)
#     predictions = Dict("aeroplane"=>[],"bicycle"=>[],"bird"=>[], "boat"=>[],
#                     "bottle"=>[],"bus"=>[],"car"=>[],"cat"=>[],"chair"=>[],
#                     "cow"=>[],"diningtable"=>[],"dog"=>[],"horse"=>[],"motorbike"=>[],
#                     "person"=>[],"pottedplant"=>[],"sheep"=>[],"sofa"=>[],"train"=>[],"tvmonitor"=>[])
#     apdic= Dict("aeroplane"=>0.0,"bicycle"=>0.0,"bird"=>0.0, "boat"=>0.0,
#                     "bottle"=>0.0,"bus"=>0.0,"car"=>0.0,"cat"=>0.0,"chair"=>0.0,
#                     "cow"=>0.0,"diningtable"=>0.0,"dog"=>0.0,"horse"=>0.0,"motorbike"=>0.0,
#                     "person"=>0.0,"pottedplant"=>0.0,"sheep"=>0.0,"sofa"=>0.0,"train"=>0.0,"tvmonitor"=>0.0)
#
#     println("Calculating accuracy...")
#     progress!(acc(model,data,conf_thresh,iou_thresh,iou,predictions))
#     for key in keys(predictions)
#         sort!(predictions[key], by = x ->x[2],rev = true)
#         tp = 0
#         fp = 0
#         total = totaldic[key]
#         preRecall = []
#         p = predictions[key]
#         for i in 1:length(p)
#             if p[i][1]
#                 tp = tp+1
#             else
#                 fp = fp+1
#             end
#             if total==0
#                 push!(preRecall,[0,0])
#             end
#             push!(preRecall,[tp/(tp+fp),tp/total])
#         end
#         #smooth process
#         rightMax = preRecall[length(preRecall)][1]
#         location = length(preRecall)-1
#         while(location >= 1)
#             if preRecall[location][1] > rightMax
#                 rightMax = preRecall[location][1]
#             else
#                 preRecall[location][1] = rightMax
#             end
#             location = location -1
#         end
#         #make calculation
#         sum = 0
#         for i in 2:length(preRecall)
#             sum = sum + (preRecall[i][2]-preRecall[i-1][2]) * preRecall[i][1]
#         end
#             apdic[key] = sum
#     end
#     println("Calculated")
#     return apdic
# end
#
# #Checks if given prediction is true positive or false negative
# function istrue(prediction,labels,check,iou)
#     min = iou
#     result = false
#     location = length(labels) + 1
#     for i in 1:length(labels)
#         if prediction[5] == namesdic[labels[i][5]] && check[i] == 0 && ioumatch(prediction[1],prediction[2],prediction[3],prediction[4],labels[i][1],labels[i][2],labels[i][3],labels[i][4]) > min
#             min = ioumatch(prediction[1],prediction[2],prediction[3],prediction[4],labels[i][1],labels[i][2],labels[i][3],labels[i][4])
#             result = true
#             location = i
#         end
#     end
#     return result ,location
# end
#
# function calculatemean(dict)
#     sum = 0
#     number = 0
#     for key in keys(dict)
#         sum = sum + dict[key]
#         number = number + 1
#     end
#     return sum/number
# end


text:
using LightGraphs
using LightGraphsIO

@testset "Types" begin
    dg = complete_digraph(5)
    g = complete_graph(5)

    # DOT
    @test typeof(DOTFormat()) <: AbstractGraphFormat
    @test load_graph(1, DOTFormat) == 1
    @test_broken write_graph(g, DOTFormat) == 1

    # EdgeList
    @test typeof(EdgeListFormat()) <: AbstractGraphFormat
    @test load_graph(1, EdgeListFormat) == 2
    @test write_graph(g, EdgeListFormat) == 2

    # GML
    @test typeof(GMLFormat()) <: AbstractGraphFormat
    @test load_graph(1, GMLFormat) == 3
    @test write_graph(g, GMLFormat) == 3

    # NET
    @test typeof(NETFormat()) <: AbstractGraphFormat
    @test load_graph(1, NETFormat) == 4
    @test write_graph(g, NETFormat) == 4

    # GEXF
    @test typeof(GEXFFormat()) <: AbstractGraphFormat
    @test load_graph(1, GEXFFormat) == 5
    @test write_graph(g, GEXFFormat) == 5

    # Graph6
    @test typeof(Graph6Format()) <: AbstractGraphFormat
    @test load_graph(1, Graph6Format) == 6
    @test write_graph(g, Graph6Format) == 6

    # GraphML
    @test typeof(GraphMLFormat()) <: AbstractGraphFormat
    @test load_graph(1, GraphMLFormat) == 7
    @test write_graph(g, GraphMLFormat) == 7

    # CDF
    @test typeof(CDFFormat()) <: AbstractGraphFormat
    @test load_graph(1, CDFFormat) == 8
    @test_broken write_graph(g, CDFFormat) == 8

    # LGCompressed
    @test typeof(LGCompressedFormat()) <: AbstractGraphFormat
    @test load_graph(1, LGCompressedFormat) == 9
    @test write_graph(g, LGCompressedFormat) == 9
end


text:
# include("AlgAss/AlgAssTypes.jl")
include("AlgAss/AbsAlgAss.jl")
include("AlgAss/AlgAss.jl")
include("AlgAss/AlgGrp.jl")
include("AlgAss/AlgMat.jl")
include("AlgAss/Elem.jl")
include("AlgAss/AlgMatElem.jl")
include("AlgAss/Map.jl")
include("AlgAss/Ideal.jl")


text:
using FixedPointNumbers, Test

function test_op(fun::F, ::Type{T}, fx, fy, fxf, fyf, tol) where {F,T}
    # Make sure that the result is representable
    (typemin(T) <= fun(fxf, fyf) <= typemax(T)) || return nothing
    @assert abs(fun(fx, fy) - convert(T, fun(fxf, fyf))) <= tol
    @assert abs(convert(Float64, fun(fx, fy)) - fun(fxf, fyf)) <= tol
end

function test_fixed(::Type{T}, f) where {T}
    values = [-10:0.01:10; -180:.01:-160; 160:.01:180]
    tol = 2.0^-f
    for x in values
        # Ignore values outside the representable range
        # typemin <, otherwise for -(-0.5)  > typemax
        if !(typemin(T) < x <= typemax(T))
            continue
        end
        # isinteger(x) && @show x
        fx = convert(T,x)
        @test convert(T,convert(Float64, fx)) == fx
        @test convert(T,convert(Float64, -fx)) == -fx
        @test convert(Float64, -fx) == -convert(Float64, fx)

        fxf = convert(Float64, fx)

        rx = convert(Rational{BigInt},fx)
        @assert isequal(fx,rx) == isequal(hash(fx),hash(rx))

        for y in values
            if !(typemin(T) < y <= typemax(T))
                continue
            end

            fy = convert(T,y)
            fyf = convert(Float64, fy)

            @assert fx==fy || x!=y
            @assert fx<fy  || x>=y
            @assert fx<=fy || x>y

            test_op(+, T, fx, fy, fxf, fyf, tol)
            test_op(-, T, fx, fy, fxf, fyf, tol)
            test_op(*, T, fx, fy, fxf, fyf, tol)
            fy != 0 && test_op(/, T, fx, fy, fxf, fyf, tol)

            @assert isequal(fx,fy) == isequal(hash(fx),hash(fy))
        end
    end
end

@testset "conversion" begin
    @test isapprox(convert(Fixed{Int8,7}, 0.8), 0.797, atol=0.001)
    @test isapprox(convert(Fixed{Int8,7}, 0.9), 0.898, atol=0.001)
    @test_throws InexactError convert(Fixed{Int8, 7}, 0.999)
    @test_throws InexactError convert(Fixed{Int8, 7}, 1.0)
    @test_throws InexactError convert(Fixed{Int8, 7}, 1)
    @test_throws InexactError convert(Fixed{Int8, 7}, 2)
    @test_throws InexactError convert(Fixed{Int8, 7}, 128)
end

@testset "test_fixed" begin
    for (TI, f) in [(Int8, 8), (Int16, 8), (Int16, 10), (Int32, 16)]
        T = Fixed{TI,f}
        # println("  Testing $T")
        test_fixed(T, f)
    end
end

@testset "modulus" begin
    T = Fixed{Int8,7}
    for i = -1.0:0.1:typemax(T)
        @test i % T === T(i)
    end
    @test ( 1.5 % T).i == round(Int,  1.5*128) % Int8
    @test (-0.3 % T).i == round(Int, -0.3*128) % Int8

    T = Fixed{Int16,9}
    for i = -64.0:0.1:typemax(T)
        @test i % T === T(i)
    end
    @test ( 65.2 % T).i == round(Int,  65.2*512) % Int16
    @test (-67.2 % T).i == round(Int, -67.2*512) % Int16
end

@testset "testapprox" begin
    for T in [Fixed{Int8,7}, Fixed{Int16,8}, Fixed{Int16,10}]
        testapprox(T)  # defined in ufixed.jl
    end
end

@testset "reductions" begin
    F8 = Fixed{Int8,8}
    a = F8[0.498, 0.1]
    acmp = Float64(a[1]) + Float64(a[2])
    @test sum(a) == acmp
    @test sum(a, dims=1) == [acmp]

    F6 = Fixed{Int8,6}
    a = F6[1.2, 1.4]
    acmp = Float64(a[1])*Float64(a[2])
    @test prod(a) == acmp
    @test prod(a, dims=1) == [acmp]
end

@testset "convert result type" begin
    x = Fixed{Int8,8}(0.3)
    for T in (Float16, Float32, Float64, BigFloat)
        y = convert(T, x)
        @test isa(y, T)
    end
end

@testset "Integer conversions" begin
    @test convert(Int, Q1f6(1)) === 1
    @test convert(Integer, Q1f6(1)) === Int8(1)
end

@testset "Floating-point conversions" begin
    @test isa(float(one(Fixed{Int8,6})),   Float32)
    @test isa(float(one(Fixed{Int32,18})), Float64)
    @test isa(float(one(Fixed{Int32,25})), Float64)
end

@testset "Show" begin
    x = Fixed{Int32,5}(0.25)
    iob = IOBuffer()
    show(iob, x)
    str = String(take!(iob))
    @test str == "0.25Q26f5"
    @test eval(Meta.parse(str)) == x
end

@testset "rand" begin
    for T in (Fixed{Int8,8}, Fixed{Int16,8}, Fixed{Int16,10}, Fixed{Int32,16})
        a = rand(T)
        @test isa(a, T)
        a = rand(T, (3, 5))
        @test ndims(a) == 2 && eltype(a) == T
        @test size(a) == (3,5)
    end
end

@testset "floatmin" begin
    # issue #79
    @test floatmin(Q11f4) == Q11f4(0.06)
end

@testset "Disambiguation constructors" begin
    @test_throws ArgumentError Fixed{Int32,16}('a')
    @test_throws InexactError  Fixed{Int32,16}(complex(1.0, 1.0))
    @test Fixed{Int32,16}(complex(1.0, 0.0))        == 1
    @test Fixed{Int32,16}(Base.TwicePrecision(1.0, 0.0)) == 1
end

@testset "fractional fixed-point numbers" begin
    # test all-fractional fixed-point numbers (issue #104)
    for (T, f) in ((Int8, 7),
                 (Int16, 15),
                 (Int32, 31),
                 (Int64, 63))
        tmax = typemax(Fixed{T, f})
        @test tmax == BigInt(typemax(T)) / BigInt(2)^f
        tol = (tmax + BigFloat(1.0)) / (sizeof(T) * 8)
        for x in range(-1, stop=BigFloat(tmax)-tol, length=50)
            @test abs(Fixed{T, f}(x) - x) <= tol
        end
    end
end

@testset "Promotion within Fixed" begin
    @test @inferred(promote(Q0f7(0.25), Q0f7(0.75))) ===
        (Q0f7(0.25), Q0f7(0.75))
    @test @inferred(promote(Fixed{Int16,3}(0.25), Fixed{Int8,3}(0.875))) ===
        (Fixed{Int16,3}(0.25), Fixed{Int16,3}(0.875))
    @test @inferred(promote(Fixed{Int8,6}(0.125), Fixed{Int8,4}(0.75))) ===
        (Fixed{Int16,6}(0.125), Fixed{Int16,6}(0.75))

    @test Fixed{Int16,15}(-1)   == Fixed{Int8,7}(-1)
    @test Fixed{Int16,15}(0.25) == Fixed{Int8,7}(0.25)
    @test Fixed{Int16,7}(-1)   == Fixed{Int8,7}(-1)
    @test Fixed{Int16,7}(0.25) == Fixed{Int8,7}(0.25)
    @test Fixed{Int16,15}(-1)   == Fixed{Int8,5}(-1)
    @test Fixed{Int16,15}(5/32) == Fixed{Int8,5}(5/32)
    @test Fixed{Int16,3}(-1)   == Fixed{Int8,5}(-1)
    @test Fixed{Int16,3}(0.25) == Fixed{Int8,5}(0.25)

    @test promote_type(Q0f7,Float32,Int) == Float32
    @test promote_type(Q0f7,Int,Float32) == Float32
    @test promote_type(Int,Q0f7,Float32) == Float32
    @test promote_type(Int,Float32,Q0f7) == Float32
    @test promote_type(Float32,Int,Q0f7) == Float32
    @test promote_type(Float32,Q0f7,Int) == Float32
    @test promote_type(Q0f7,Q1f6,Q2f5,Q3f4,Q4f3,Q5f2) == Fixed{Int128,7}
end


text:
using Dates
using OrderedCollections
using Setfield
using DataFrames

import Agtor: Climate

@with_kw mutable struct Crop <: AgComponent
    name::String
    crop_type::String

    # Growth pattern
    growth_stages::NamedTuple
    coef_stages::NamedTuple

    # Crop properties
    plant_date::Union{Date, AgParameter}
    yield_per_ha::Union{Float64, AgParameter}
    price_per_yield::Union{Float64, AgParameter}
    variable_cost_per_ha::Union{Float64, AgParameter}
    water_use_ML_per_ha::Union{Float64, AgParameter}
    root_depth_m::Union{Float64, AgParameter}
    et_coef::Union{Float64, AgParameter}
    wue_coef::Union{Float64, AgParameter}
    rainfall_threshold::Union{Float64, AgParameter}
    ssm_coef::Union{Float64, AgParameter}
    effective_root_zone::Union{Float64, AgParameter}
    harvest_date::Union{Date, AgParameter}
    harvest_offset::Union{Day, Int64, AgParameter}
    naive_crop_income::Float64

    function Crop(name::String, crop_type::String, growth_stages::Dict, start_dt::Date; props...)::Crop
        plant_date = props[:plant_date]
        yield_per_ha = props[:yield_per_ha]
        price_per_yield = props[:price_per_yield]
        variable_cost_per_ha = props[:variable_cost_per_ha]
        water_use_ML_per_ha = props[:water_use_ML_per_ha]
        root_depth_m = props[:root_depth_m]
        et_coef = props[:et_coef]
        wue_coef = props[:wue_coef]
        rainfall_threshold = props[:rainfall_threshold]
        ssm_coef = props[:ssm_coef]
        effective_root_zone = props[:effective_root_zone]

        start_year = year(start_dt)
        sow_date = Date("$(start_year)-$(plant_date)")

        h_day = 0
        offset = 0
        start_date = sow_date
        g_stages = NamedTuple()
        coef_stages = NamedTuple()

        @inbounds for (k, v) in growth_stages
            offset = v[:stage_length]
            end_of_stage = start_date + Dates.Day(offset)
            @set! g_stages[Symbol(k)] = (
                start = start_date,
                var"end" = end_of_stage,
                stage_length = offset
            )

            @set! coef_stages[k] = NamedTuple{Tuple(keys(v))}(values(v))

            h_day += offset
            start_date = end_of_stage + Dates.Day(offset+1)

        end

        harvest_offset = Dates.Day(h_day)
        harvest_date = sow_date + harvest_offset

        c = new(name, crop_type, g_stages, coef_stages, sow_date,
                yield_per_ha, price_per_yield, variable_cost_per_ha, water_use_ML_per_ha,
                root_depth_m, et_coef, wue_coef, rainfall_threshold,
                ssm_coef, effective_root_zone, harvest_date, harvest_offset)

        # variable cost does not include water usage costs (added later in optimized allocation step)
        c.naive_crop_income = (c.price_per_yield * c.yield_per_ha) - c.variable_cost_per_ha

        return c
    end
end

"""Update growth stages with corresponding dates from given sowing date."""
function update_stages!(c::Crop, dt::Date)::Nothing
    stages::NamedTuple = c.growth_stages::NamedTuple
    start_date::Date = Date(yearmonthday(dt)...)
    @inbounds for (k, v) in pairs(c.growth_stages)
        offset::Int64 = v[:stage_length]
        end_of_stage = start_date + Dates.Day(offset)
        @set! stages[k] = (
            start = start_date,
            var"end" = end_of_stage,
            stage_length = offset
        )

        start_date = end_of_stage + Dates.Day(offset+1)
    end

    c.growth_stages = stages

    return nothing
end

function get_stage_coefs(c::Crop, dt::Date)::NamedTuple
    @inbounds for (k::Symbol, v::NamedTuple) in pairs(c.growth_stages::NamedTuple)
        s::Date = v[:start]
        e::Date = v[:end]

        # if in season...
        if (s <= dt <= e)
            return c.coef_stages[k]
        end
    end

    # Not in season so just return initial growth stage
    return c.coef_stages[:initial]
end

"""Naive estimation of net income."""
function estimate_income_per_ha(c::Crop)::Float64
    return (c.price_per_yield * c.yield_per_ha) - c.variable_cost_per_ha
end

function subtotal_costs(c::Crop, year::Int64)::Float64
    # cost of production is handled by factoring in
    # water application costs and other maintenance costs.
    # The variable_cost_per_ha is only used to inform estimates.
    return 0.0
end

function create(spec::Dict, start_dt::Date)::Crop
    data = deepcopy(spec)
    _ = pop!(data, :component)
    return Crop(data[:name], data[:crop_type], data[:growth_stages], start_dt; data...)
end

"""Get all in-season dates for a given crop."""
function in_season_dates(c::Climate, crop::Crop)::DataFrame
    d = c.data
    sow_d = monthday(crop.plant_date)
    harvest_d = monthday(crop.harvest_date)
    mds = monthday.(d[:Date])

    return d[(mds .>= [sow_d]) .| (mds .< [harvest_d]), :]
end

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
#%%
using GdMn2O5
using Parameters
using Base.Threads
using DataFrames
using FileIO
using DrWatson
using Plots
using Query
using Statistics
#%%
df = load(datadir("sims.jld2"))["df"]

function window_pbs(pbs, hmags)
    window_pbs = Float64[]
    window_hmags = Float64[]
    count = 0
    for (h, p) in zip(hmags, pbs)
        if 1 <= h <= 2.5
            push!(window_hmags, h)
            push!(window_pbs, p)
        end
    end
    return window_pbs, window_hmags
end

function Pbs_within_half(m, windowPbs, upper_m, lower_m)
    out = Float64[]
    ids = Int[]
    for (i, p) in enumerate(windowPbs)
        if (p < m && (m - p)/3 < p - lower_m)  || ( p >= m && p - m < upper_m - p)
            push!(out, p)
            push!(ids, i)
        end
    end
    return out, ids
end

function sanitize_Pb!(Pbs)
    m = mean(Pbs)
    upper_branch = filter(x-> x >= m, Pbs)
    lower_branch = filter(x-> x <  m, Pbs)
    upper_m = mean(upper_branch)
    lower_m = mean(lower_branch)
    for i in reverse(1:length(Pbs)-1)
        if (Pbs[i] > upper_m && Pbs[i+1] > upper_m && abs((Pbs[i] - upper_m) - (Pbs[i+1] - upper_m)) > 1) || (Pbs[i] < lower_m && Pbs[i+1] < lower_m && abs((lower_m - Pbs[i])-(lower_m - Pbs[i+1])) > 1)
            Pbs[i] = Pbs[i+1]
        end
    end
    return Pbs
end

function m_lower_upper(points)
    m = mean(points)
    return (m, mean(filter(x->x<m, points)), mean(filter(x->x>=m, points)))
end

#%%

t=df |>
    @groupby({_.J1, _.J2, _.K_Gd, _.K_L , _.γ ,_.dipole_strength})|>
    @map(x-> begin
        ih1 = findfirst(y -> y == deg2rad(0.0), x.Hϕ)
        ih2 = findfirst(y -> y == deg2rad(6), x.Hϕ)

        pbs1 = sanitize_Pb!(x.Pbs[ih1])
        pbs2 = sanitize_Pb!(x.Pbs[ih2])
        m1, lower_m1, upper_m1 = m_lower_upper(pbs1)
        m2, lower_m2, upper_m2 = m_lower_upper(pbs2)

        wpbs1, window_mags1 = window_pbs(pbs1, x.Hmags[1])
        wpbs2, window_mags2 = window_pbs(pbs2, x.Hmags[1])
        # sanitize_Pb!(wpbs1)
        # sanitize_Pb!(wpbs2)
        if !isempty(wpbs1) &&  !isempty(wpbs2)
            pbs_in_half1, ids1 = Pbs_within_half(mean(wpbs1), wpbs1, upper_m1, lower_m1)
            pbs_in_half2, ids2 = Pbs_within_half(mean(wpbs2), wpbs2, upper_m2, lower_m2)
            N = length(pbs_in_half2) - length(pbs_in_half1)
            Hmags1 = window_mags1[ids1]
            Hmags2 = window_mags2[ids2]
        else
            pbs_in_half1, ids1 = Float64[], Int[]
            pbs_in_half2, ids2 = Float64[], Int[]
            Hmags1 = Float64[]
            Hmags2 = Float64[]
            N = 0
        end

        window_m1, window_lower_m1, window_upper_m1 = m_lower_upper(wpbs1)
        window_value = (window_upper_m1 - window_m1) + (window_m1 - window_lower_m1)
        window_rank = window_value > 10 ? 1 : 0

        return (J1             =x.J1[1],
                J2             =x.J2[1],
                K_Gd           =x.K_Gd[1],
                K_L            =x.K_L[1],
                γ              =x.γ[1],
                χ_Gd           = x.χ_Gd[1],
                χ_L           = x.χ_L[1],
                g           = x.g[1],
                α           = x.α[1],
                dipole_strength=x.dipole_strength[1],
                Pbs1           =x.Pbs[ih1],
                Pbs2           =x.Pbs[ih2],
                pbs_in_half1   = pbs_in_half1,
                pbs_in_half2   = pbs_in_half2,
                Hmags1   = Hmags1,
                Hmags2   = Hmags2,
                m1 = m1,
                m2 = m2,
                Hmags          = x.Hmags[1],
                N              =window_rank*N,
                ϕ1s_1 = x.ϕ1s[ih1],
                ϕ1s_2 = x.ϕ1s[ih2],
                ϕ2s_1 = x.ϕ2s[ih1],
                ϕ2s_2 = x.ϕ2s[ih2],
                ϕ3s_1 = x.ϕ3s[ih1],
                ϕ3s_2 = x.ϕ3s[ih2],
                ϕ4s_1 = x.ϕ4s[ih1],
                ϕ4s_2 = x.ϕ4s[ih2],
                ϕ5s_1 = x.ϕ5s[ih1],
                ϕ5s_2 = x.ϕ5s[ih2],
                ϕ6s_1 = x.ϕ6s[ih1],
                ϕ6s_2 = x.ϕ6s[ih2],
                ϕ7s_1 = x.ϕ7s[ih1],
                ϕ7s_2 = x.ϕ7s[ih2],
                ϕ8s_1 = x.ϕ8s[ih1],
                ϕ8s_2 = x.ϕ8s[ih2],
                L1s_1 = x.L1ϕs[ih1],
                L1s_2 = x.L1ϕs[ih2],
                L2s_1 = x.L2ϕs[ih1],
                L2s_2 = x.L2ϕs[ih2],
                )
    end) |>
    @orderby_descending(_.N) |> DataFrame
#%%
# using CSV
# CSV.write(datadir("sims_processed.csv"), t[1:10, :])

# Dict(pairs(t))

# function plot_row(t, i)
#     p1 = scatter(t[i, :Hmags], t[i, :Pbs1], color=:blue)
#     scatter!(p1, t[i, :Hmags1], t[i, :pbs_in_half1], color=:red)
#     plot!(p1, x -> t[i, :m1], color=:green)
#     p2 = scatter(t[i, :Hmags], t[i, :Pbs2], color=:blue)
#     scatter!(p2, t[i, :Hmags2], t[i, :pbs_in_half2], color=:red)
#     plot!(p2, x -> t[i, :m2], color=:green)
#     # p3 = scatter(t[i, :Hmags], t[i, :ϕ1s_1], color=:green, title="ϕs1")
#     # p4 = scatter(t[i, :Hmags], t[i, :ϕ5s_1], color=:purple, title="ϕs5")
#     # p5 = scatter(t[i, :Hmags], t[i, :L1s_1], color=:purple, title="Lϕs1")
#     # plot(p1, p2, p3, p4, p5)
#     plot(p1, p2)
# end

# plot_row(t,50)

# #%%
# using VegaLite

# t |> @take(20) |>
#     @vlplot(:point, x=:K_Gd, y=:γ, color=:dipole_strength, column=:N)

# plot_row(t, 1)
# #%%

# t[:, :N]
# plot(scatter(t[1:20, :K_Gd]./t[1:20,:dipole_strength]))
# plot(scatter((t[1:20, :J1] ./ t[1:20,:γ]).*(t[1:20,:K_Gd]./t[1:20,:K_L])))

# (t[1:20, :J1] ./ t[1:20,:γ]).*(t[1:20,:K_Gd]./t[1:20,:K_L])

# t[1:4, :]
# t[2, :]

# # {1,2,3}
# # df = hcat(df, df |>  |> DataFrame)
# # df = hcat(df, df |> @map({N_within_half_Pb=length(_.within_half_Pb)}) |> @orderby(_.N_within_half_Pb) |> DataFrame)

# # df[:, :Hϕ]
# # df1 = df |> @filter(_.Hϕ == deg2rad(0)) |> DataFrame
# # df2 = df |> @filter(_.Hϕ == deg2rad(6)) |> DataFrame
# # df3 = df |> @filter(_.Hϕ == deg2rad(8)) |> DataFrame

# t[end,:]

# #%%

# df[end, :N_within_half_Pb]
# df
# scatter(df[end-10, :Hmags], df[end-10, :Pbs])
# sortperm(df[1, :Hmags])
# scatter(df[1, :Hmags], df[1, :Pbs])
# for i = 8000:8200
#     display(scatter(df[i, :Hmags], df[i, :Pbs], title="$i"))
# end

# findall(x -> x==6, df[1, :Hmags])
# df[1, :Pbs][50]
# scatter(df[:, :K_Gd])

# #%%


text:
using Trixi, OrdinaryDiffEq

surface_flux = FluxLaxFriedrichs()
volume_flux = flux_ranocha
dg = DGMulti(polydeg = 3, element_type = Line(), approximation_type = Polynomial(),
             surface_integral = SurfaceIntegralWeakForm(surface_flux),
             volume_integral = VolumeIntegralFluxDifferencing(volume_flux))

equations = CompressibleEulerEquations1D(1.4)

initial_condition = initial_condition_convergence_test
source_terms = source_terms_convergence_test

cells_per_dimension = (8, )
vertex_coordinates, EToV = StartUpDG.uniform_mesh(dg.basis.elementType, cells_per_dimension...)
mesh = VertexMappedMesh(vertex_coordinates, EToV, dg, is_periodic=(true,))
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg;
                                    source_terms=source_terms)

tspan = (0.0, 1.1)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()
alive_callback = AliveCallback(alive_interval=10)
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt = 0.5 * estimate_dt(mesh, dg), save_everystep=false, callback=callbacks);

summary_callback() # print the timer summary


text:
n = rand(1:7)
@test matrixdepot("sampling", n) == matrixdepot("sampling", Float64, n)

A = matrixdepot("sampling", n)
v = sort(eigvals(A))
@test_approx_eq_eps v [0: n - 1;] 1e-10
println("'sampling' passed test...")