text:
using ImageFiltering

@testset "MS-SSIM" begin
    @info "test: MS-SSIM"

    # constructors
    iqi_a = MSSSIM()
    iqi_b = MSSSIM(KernelFactors.gaussian(1.5, 11))
    iqi_c = MSSSIM(KernelFactors.gaussian(1.5, 11), (0.0448, 0.2856, 0.3001, 0.2363, 0.1333))
    iqi_d = MSSSIM(KernelFactors.gaussian(1.5, 11), (
        (0.0448, 0.0448, 0.0448),
        (0.2856, 0.2856, 0.2856),
        (0.3001, 0.3001, 0.3001),
        (0.2363, 0.2363, 0.2363),
        (0.1333, 0.1333, 0.1333),
    ))
    iqi_e = MSSSIM(; num_scales=5)
    iqi_f = @suppress_err MSSSIM(KernelFactors.gaussian(1.5, 11), 2 .* (0.0448, 0.2856, 0.3001, 0.2363, 0.1333))
    @test iqi_a == iqi_b == iqi_c == iqi_d == iqi_e == iqi_f

    output = @capture_err begin
        MSSSIM(KernelFactors.gaussian(1.5, 11), (0.0448, 0.2856, 0.3001))
    end
    @test output == ""
    output = @capture_err begin
        MSSSIM(KernelFactors.gaussian(1.5, 11), (0.1, 0.3, 0.6); num_scales=1)
    end
    @test occursin("truncate MS-SSIM weights to scale", output)
    @test occursin("normalize MS-SSIM weights so that (∑α, ∑β, ∑γ) == (1.0, 1.0, 1.0)", output)

    # Gray test images
    img1 = testimage("fabio_gray_512")
    img2 = testimage("lena_gray_512")
    # Tensorflow, pytorch_msssim and MATLAB original implementation: 0.1072
    # We have strong reason to suspect that MSSSIM in MATLAB R2020a is incorrect
    # https://github.com/JuliaImages/ImageQualityIndexes.jl/pull/19#issuecomment-655208537
    @test assess_msssim(img1, img2) ≈ 0.11109 atol=1e-4
    @test assess(MSSSIM(), img1, img2) ≈ MSSSIM()(img1, img2)

    # RGB test images
    # MS-SSIM in MATLAB R2020a and its original implementation only support Gray images
    img3 = testimage("mandril_color")
    img4 = testimage("lena_color_512")
    @test assess_msssim(img3, img4) ≈ 0.08516 atol=1e-4 # pytorch_msssim: 0.0786
    @test assess(MSSSIM(), img3, img4) ≈ MSSSIM()(img3, img4)

    # Varying αᵢ, βᵢ, γᵢ
    iqi_e = @suppress_err MSSSIM(KernelFactors.gaussian(1.5, 11), (
        (0.0448, 0.2856, 0.3001),
        (0.3001, 0.0448, 0.2856),
        (0.2856, 0.3001, 0.0448),
    ))
    @test iqi_e(img3, img4) ≈ 0.13109 atol=1e-4
    # non standard parameters, result may differ from other implementations

    # Comparing with SSIM
    iqi_f = MSSSIM(KernelFactors.gaussian(1.5, 11), (1, ))
    @test iqi_f(img3, img4) ≈ SSIM()(img3, img4) atol=1e-3

    # images should have the same size, but their axes can differs
    @test assess_msssim(img1, OffsetArray(img2, -2, -2)) == assess_msssim(img1, img2)
    @test_throws ArgumentError assess_msssim(img1, restrict(img2))

    # RGB - Gray
    @test assess_msssim(img3, Gray.(img4)) == assess_msssim(img3, RGB.(Gray.(img4)))
    @test assess_msssim(img3, img4) ≠ assess_msssim(channelview(img3), channelview(img4))

    # Gray type tests
    type_list = generate_test_types([Float32, N0f8], [Gray])
    A = rand(128,128)
    B = rand(128,128)

    for T in type_list

        a = A .|> T
        b = B .|> T

        @test assess_msssim(a, b) == assess(MSSSIM(), a, b) == MSSSIM()(a, b)
        @test assess_msssim(a, a) ≈ 1.0 atol=1e-4
        @test assess_msssim(a, b; num_scales=1) == assess(MSSSIM(num_scales=1), a, b) == MSSSIM(num_scales=1)(a, b)
    end

    # check if numbers are treated the same as gray colorants
    a = Gray.(A)
    b = Gray.(B)
    @test assess_msssim(a, b) ≈ assess_msssim(A, B)

    # RGB type tests
    type_list = generate_test_types([Float32, N0f8], [RGB])
    A = [RGB(0.0, 0.0, 0.0) RGB(0.0, 1.0, 0.0) RGB(0.0, 1.0, 1.0)
        RGB(0.0, 0.0, 1.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 1.0, 0.0)
        RGB(1.0, 1.0, 1.0) RGB(1.0, 0.0, 1.0) RGB(0.0, 0.0, 0.0)]
    B = [RGB(0.0, 0.0, 0.0) RGB(0.0, 0.0, 1.0) RGB(1.0, 1.0, 1.0)
        RGB(0.0, 1.0, 0.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 0.0, 1.0)
        RGB(0.0, 1.0, 1.0) RGB(1.0, 1.0, 0.0) RGB(0.0, 0.0, 0.0)]

    for T in type_list

        a = A .|> T
        b = B .|> T

        @test assess_msssim(a, b) == assess(MSSSIM(), a, b) == MSSSIM()(a, b)
        @test assess_msssim(a, a) ≈ 1.0 atol=1e-4
        @test assess_msssim(a, b; num_scales=1) == assess(MSSSIM(num_scales=1), a, b) == MSSSIM(num_scales=1)(a, b)
    end

    # Other Color3 type tests
    type_list = generate_test_types([Float32], [Lab, HSV])
    A = [RGB(0.0, 0.0, 0.0) RGB(0.0, 1.0, 0.0) RGB(0.0, 1.0, 1.0)
        RGB(0.0, 0.0, 1.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 1.0, 0.0)
        RGB(1.0, 1.0, 1.0) RGB(1.0, 0.0, 1.0) RGB(0.0, 0.0, 0.0)]
    B = [RGB(0.0, 0.0, 0.0) RGB(0.0, 0.0, 1.0) RGB(1.0, 1.0, 1.0)
        RGB(0.0, 1.0, 0.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 0.0, 1.0)
        RGB(0.0, 1.0, 1.0) RGB(1.0, 1.0, 0.0) RGB(0.0, 0.0, 0.0)]
    for T in type_list
        a = A .|> T
        b = B .|> T

        @test_nowarn assess_msssim(A, b), assess_msssim(a, B)

        @test assess_msssim(a, b) == assess(MSSSIM(), a, b) == MSSSIM()(a, b)
        @test assess_msssim(a, a) ≈ 1.0 atol=1e-4
        @test assess_msssim(a, b; num_scales=1) == assess(MSSSIM(num_scales=1), a, b) == MSSSIM(num_scales=1)(a, b)

    end
    @test assess_msssim(A, B) ≈ assess_msssim(Lab.(A), B) atol=1e-4

    type_list = generate_test_types([Float32, N0f8], [RGB, BGR])
    test_cross_type(MSSSIM(), A, B, type_list)
end


text:
module SCF

#Julia Libraries/Modules
using Printf
using LinearAlgebra

#LearnHatreeFock.jl Modules
using TypesBasis,TypesParticles
using CoulombExchange

export runscf

@doc raw"""
function scf()

description: Self-consistent field soultion approach.

The Hatree-Fock operator for non-interacting electrons in an orthonormal basis set is give as:

$H_o C = \epsilon S C$ where $C,epsilon$ are the eigenvalues,eigenvectors and $S$ is the overlap matrix.

The first step is to find the eigenvalues and vectors of the overlap matrix. This will be used to transform from
atomic centered basis-set to a orthonormal molecular basis-set, for example, $H_{MO}=
"""
function runscf(numelec::Int,overlap::Array{Float64,2},
             ho::Array{Float64,2},elecelecrepul::Array{Float64,4};
             maxcycle=30::Int,tolerance=1.0e-6::Float64)

    convergflag = false;

    #Orthogonal overlap matrix, i.e., molecular orbitals
    oevals,oevecs = eigen(overlap);
    oevalsmat = Diagonal(oevals);

    #Try to invert eigenvalues matrix assuming non-singular, if singular
    # promote datatype ComplexF64
    invoevalsmat = returnevalsinvhalf(oevalsmat);
    orthoverlap = oevecs*invoevalsmat*adjoint(oevecs);
    orthoverlapstar = adjoint(orthoverlap);

    #Fock matrix construction
    fock = copy(ho);
    fockprime = orthoverlapstar*fock*orthoverlap;

    #HF Roothaan eq.
    #Find molecular orbital single electron solutions
    evals,evecsprime = eigen(fockprime);

    #Go back to atomic orbitals
    evecs = orthoverlap*evecsprime;

    #Sort lowest eigen values of for atomic orbitals
    sorteigen!(evals,evecs);

    #Build electron density for Hatree potential?
    density = builddensity(evecs,numelec)

    totalenergy = zeros(Float64,maxcycle+1);
    #SCF cycle
    for cycle=2:maxcycle+1

        #Get exchange and fock potential
        exchangepot = getcoulexchange(density,elecelecrepul)
        fock = ho + exchangepot;

        #find solutions
        fockprime = orthoverlapstar*fock*orthoverlap;
        evals,evecsprime = eigen(fockprime);
        evecs = orthoverlap*evecsprime;

        sorteigen!(evals,evecs);

        density = builddensity(evecs,numelec);

        totalenergy[cycle] = getfockenergy(ho,fock,density);

        #Check conditions
        delta = abs(totalenergy[cycle]-totalenergy[cycle-1]);
        printscf(cycle-1,delta,totalenergy[cycle])
        if delta < tolerance
            return totalenergy[cycle]
        end

    end #cycle=2:maxcycle+1

    return totalenergy[end];

end #runscf

@doc raw"""
function returnevalsinvhalf(evals)

Try to invert eigenvalues matrix assuming non-singular, if singular promote
datatype ComplexF64

"""
function returnevalsinvhalf(evals)
    try
        invevalsmat = evals^-0.5e0;
    catch
        invevalsmat = (convert(Array{ComplexF64},evals))^(-0.5e0);
    end
end #returnevalsinvhalf

@doc raw"""
function builddensity(evecs::Array{Number,2},n::)

Restricted HF density

"""
function builddensity(evecs::Matrix,numelec::Int)
    #Restricted HF density
    nhalf = Int(numelec/2);
    esize = size(evecs)[1];

    if typeof(evecs) == Array{Complex,2}
        density = zeros(Complex,esize,esize)
    else
        density = zeros(esize,esize)
    end
    for n=1:esize
        for m=1:esize
            for i=1:nhalf
                density[n,m] += evecs[n,i]*evecs[m,i];
            end #i=1:nhalf
        end
    end #n=1:esize
    return density

end #builddensity

@doc raw"""
function sortevals!(evals::Array,numelec::Int)

description: sort the eignevalues and vectors.

NOTES: I would like to rework this so that its more
    succinct.
"""
function sorteigen!(evals::Vector{T},evecs::Matrix{T}) where {T<:Real}
    ncol = size(evecs)[2];

    #Make a shallow copy and force local scope
    #CHECK: I'm not sure what the most Julianic way is to do this.
    local sortedevals = copy(evals);
    local sortedevecs = copy(evecs);

    #Get sorted indexes for eigenvals
    sortedindex = sortperm(evals);
    evals[:] = sortedevals[sortedindex];

    #Get sorted eigenvectors
    for i=1:ncol
        sortedevecs[:,i] = evecs[:,sortedindex[i]];
    end

    evecs[:,:] = sortedevecs[:,:];

end #sortevals

function sorteigen!(evals::Vector{T},evecs::Matrix{T}) where {T<:Complex}
    ncol = size(evecs)[2];

    #Make a shallow copy and force local scope
    #CHECK: I'm not sure what the most Julianic way is to do this.
    local sortedevals = copy(evals);
    local sortedevecs = copy(evecs);

    magevals = adjoint(evals)*evals;
    println(magevals);
    #Get sorted indexes for eigenvals
    sortedindex = sortperm(real(magevals));
    evals[:] = sortedevals[sortedindex];

    #Get sorted eigenvectors
    for i=1:ncol
        sortedevecs[:,i] = evecs[:,sortedindex[i]];
    end

    evecs[:,:] = sortedevecs[:,:];

end #sortevals

@doc raw"""
sorteigen()

TODO
"""
function sorteigen(evals::Vector{T},evecs::Matrix{T}) where {T<:Real}
    p = sortperm(evals);
    return evals[p], evecs[p,:]
end

@doc raw"""
function getfockenergy()

"""
function getfockenergy(ho::Array,fock::Array,density::Array)
    numbasis = size(density)[1];
    energy = 0.00e0;
    for n=1:numbasis
        for m=1:numbasis
            energy += density[n,m] * (ho[n,m]+fock[n,m]);
        end
    end#n=1:numbasis
    return energy
end #getfockenergy

@doc raw"""
function printscf()

"""
function printscf(iteration::Int,delta::Float64,energy::Float64)
    if iteration < 2
        @printf("                   HF SCF                         \n");
        @printf("--------------------------------------------------\n");
        @printf("Iteration         ΔE [Ha]          HF-Energy [Ha] \n");
        @printf("--------------------------------------------------\n");
        @printf("%i                %0.6f           %0.6e   \n",
            iteration,delta,energy);
    else
        @printf("%i                %0.6f           %0.6e   \n",
            iteration,delta,energy);
    end
end #printscf

end #SCF

ext: jl


text:
using Test
using ThermodynamicsBase

struct mything <: AbstractEquationOfState
    name::AbstractString
end
struct mything2 <: AbstractPhaseState
    name::AbstractString
end
struct mything3 <: AbstractProperty
    name::AbstractString
end

Base.show(a::mything) = a.name
Base.show(b::mything2) = b.name
Base.show(c::mything3) = c.name

#function thermocalc(st::AbstractPhaseState, eos::AbstractEquationOfState, p::AbstractProperty)
#println("We're sorry, the property $show(p) isn't implemented for $show(st) , using the $show(eos) equation of state.") end
@test thermocalc(mything2("Generic Phase State"), mything("Generic EOS"), mything3("Generic Property")) == "We're sorry, the property Generic Property isn't implemented for Generic State , using the Generic EOS equation of state."


text:
#Función de inflación de exclusión dinámica

#InflationDynamicExclusion.jl - Definición de exclusión dinámica

"""
    InflationDynamicExclusion <: InflationFunction

Función para computar la inflación de exclusión dinámica dado el factor inferior (`lower_factor`) y factor superior (`upper_factor`).

## Utilización

    function (inflfn::InflationDynamicExclusion)(base::VarCPIBase)

Define cómo opera `InflationDynamicExclusion` sobre un objeto de tipo VarCPIBase.

    function (inflfn::InflationDynamicExclusion)(cs::CountryStructure)

Define cómo opera `InflationDynamicExclusion` sobre un objeto de tipo CountryStructure.

### Ejemplo
Cálculo el recorte dinámico simétrico (2, 2) la distribución de variaciones intermensuales ponderadas.

```julia-repl
julia> dynExfn = InflationDynamicExclusion(2, 2)
(::InflationDynamicExclusion) (generic function with 5 methods)
julia>dynExfn(gtdata) #gtdata es de tipo UniformCountryStructure
231×1 Matrix{Float32}:
 7.0480227
 7.3734045
 7.6767564
 7.6933146
 7.7317834
 ⋮
 1.0896564
 1.1332393
 1.0305643
 1.1091232
```
"""
Base.@kwdef struct InflationDynamicExclusion <: InflationFunction
    lower_factor::Float32
    upper_factor::Float32
end

# Método para recibir argumentos como par en una tupla.
InflationDynamicExclusion(factors::Tuple{Real, Real}) = InflationDynamicExclusion(
    convert(Float32, factors[1]),
    convert(Float32, factors[2])
)

# Método para recibir argumentos como par en una lista.
function (InflationDynamicExclusion)(factor_vec::Vector{<:Real})
    length(factor_vec) != 2 && return @error "Dimensión incorrectal del vector"
    InflationDynamicExclusion(
        convert(Float32, factor_vec[1]),
        convert(Float32, factor_vec[2])
    )
end

"""
    measure_name(inflfn::InflationDynamicExclusion)

Indica qué medida se utiliza para una instancia de una función de inflación.

# Ejemplo
```julia-repl
julia> dynExfn = InflationDynamicExclusion(2, 2)
julia> measure_name(dynExfn)
"Inflación de exclusión dinámica (2.0, 2.0)"
```
"""
function CPIDataBase.measure_name(inflfn::InflationDynamicExclusion)
    round_lower_factor, round_upper_factor = string.(
        round.(
            [inflfn.lower_factor, inflfn.upper_factor], digits = 2
        )
    )
    "Inflación de exclusión dinámica ($(round_lower_factor), $(round_upper_factor))"
end

"""
    measure_tag(inflfn::InflationDynamicExclusion)

Indica qué medida se utiliza para una instancia de una función de inflación.

# Ejemplo
```julia-repl
julia> dynExfn = InflationDynamicExclusion(2, 2)
julia> measure_tag(dynExfn)
"DynEx(2.0,2.0)"
```
"""
# Extendemos `params`, que devuelve los parámetros de la medida de inflación
CPIDataBase.params(inflfn::InflationDynamicExclusion) = (inflfn.lower_factor, inflfn.upper_factor)

# Las funciones sobre VarCPIBase resumen en variaciones intermensuales
function (inflfn::InflationDynamicExclusion)(base::VarCPIBase)
    lower_factor = inflfn.lower_factor
    upper_factor = inflfn.upper_factor

    std_v = std(base.v, dims = 2)
    mean_v = mean(base.v, dims = 2)

    dynEx_filter = (mean_v - (lower_factor .* std_v)) .<=
        base.v .<=
        (mean_v + (upper_factor .* std_v))

    dynEx_w = base.w' .* dynEx_filter
    dynEx_w = dynEx_w ./ sum(dynEx_w, dims = 2)

    dynEx_v = sum(
        base.v .* dynEx_w,
        dims = 2
    ) |> vec
end


text:
using CurveProximityQueries
using IntervalArithmetic
using Plots

pgfplots()

include("EllipseBounds.jl")
using .EllipseBounds

curve = Bernstein([[-0.6, 0.0], [-0.6, 0.0], [0.1, -0.4], [0.4, 1.6], [-0.2, -1.6], [-0.1, 0.4], [0.8, 0.0], [0.8, 0.0]])
ell = Ellipse(curve, 0.2..0.8)

plot(ell,
    linecolor=RGBA(102/255,130/225,223/255,0.5),
    fillcolor=RGBA(102/255,130/225,223/255,0.5),
    )
plot!(curve, 0.0..0.2,
      linewidth=2.0,
      linestyle=:dash,
      linecolor=Gray(0.5),
     )
plot!(curve, 0.2..0.8,
      linewidth=2.0,
      linecolor=:black,
     )
plot!(curve, 0.8..1.0,
      linewidth=2.0,
      linestyle=:dash,
      linecolor=Gray(0.5),
     )
scatter!(curve(0.2), markersize=4, markercolor=:black)
scatter!(curve(0.8), markersize=4, markercolor=:black)
annotate!(curve(0.5)[1]+0.03, curve(0.5)[2]+0.17, "\$\\mathcal{U}_{\\mathcal{Q}}\$")
annotate!(curve(0.2)[1], curve(0.2)[2]-0.07, "\$\\psi(\\alpha)\$")
annotate!(curve(0.8)[1], curve(0.8)[2]+0.07, "\$\\psi(\\beta)\$")
savefig("out/fig3.pdf")


text:
## Initialize
import Pkg
Pkg.activate("/Users/arreyes/.julia/environments/v1.5/")
using Plots, BenchmarkTools
gr()
Pkg.activate("/Users/arreyes/.julia/dev/RamirezReyes_Yang_SpontaneousCyclogenesis/")
using AvailablePotentialEnergyFramework
using Statistics:mean
using NCDatasets
import ImageSegmentation: seeded_region_growing, labels_map, region_adjacency_graph
## Read test file
test_file_with_cyclones = "/Users/arreyes/Documents/Research/Developement/testNetCDF/f5e-4_2km_1000km_homoRad_2d_timesteps.nc"
test_file_with_cyclones_3d = "/Users/arreyes/Documents/Research/Developement/testNetCDF/f5e-4_2km_1000km_homoRad_3d_timesteps.nc"
ds = Dataset(test_file_with_cyclones);
ds_3d = Dataset(test_file_with_cyclones_3d);
pres = ds_3d["PP"][:,:,:,:] :: Array{Float32,4};
U = ds_3d["U"][:,:,:,:] :: Array{Float32,4};
V = ds_3d["V"][:,:,:,:] :: Array{Float32,4};
close(ds)
close(ds_3d)
windspeed = sqrt.(U.*U .+ V.*V);

pres_anomaly = pres[:,:,1,:] .- mean(pres[:,:,1,:],dims=(1,2));

## Test center detection
anim = @animate for t in 1:size(windspeed,4)
    centers = @views AvailablePotentialEnergyFramework.findcyclonecenters_aspressureminima(pres_anomaly[:,:,t],-5,2000);
    centers_x= [centers[i][1] for i in 1:length(centers)];
    centers_y = [centers[i][2] for i in 1:length(centers)];
    p1 = @views heatmap(permutedims(windspeed[:,:,1,t]),title="Wind speed (m/s)",c=:viridis)
    p2 = @views heatmap(permutedims(pres_anomaly[:,:,t]),c=:viridis,title = "Pressure anomaly (hPa) and centers , timestep = $t")
    scatter!(p2,centers_x,centers_y,markersize=10,markershape=:cross,c=:red)
    plot(p1,p2,size=(800,380))
end
gif(anim, "testresults/center_detection.gif", fps = 1)

## Test cyclone detection and segmentation

anim = @animate for t in 1:size(windspeed,4)
    centers_and_labels,cyclones = AvailablePotentialEnergyFramework.detect_cyclones(pres[:,:,1,t],-5,2000)
    labels_map2 = copy(labels_map(cyclones));
    labels_map2[labels_map2 .== 1000] .= 0;
    p1 = @views heatmap(permutedims(pres_anomaly[:,:,t]),c=:viridis,title = "Pressure anomaly (hPa)")
    p2 = heatmap(permutedims(labels_map2),reuse=false,title="Segmented cyclones timestep = $t")
    plot(p1,p2,size=(800,380))
end
gif(anim, "testresults/cyclone_segmentation.gif", fps = 1)

#Test azimuthal averaging

radiusbins = 1000:2000:512000;
anim = @animate for t in 1:size(windspeed,4)
    centers_and_labels,cyclones = AvailablePotentialEnergyFramework.detect_cyclones(pres[:,:,1,t],-5,2000)
    if t == 1
        global presaverage = AvailablePotentialEnergyFramework.azimuthalaverage_allcyclones(radiusbins,pres[:,:,1,t],cyclones,centers_and_labels,2000)
        global windspeedaverage = AvailablePotentialEnergyFramework.azimuthalaverage_allcyclones(radiusbins,windspeed[:,:,:,t],cyclones,centers_and_labels,2000)
    else
        presaverage = 0.5(presaverage .+ AvailablePotentialEnergyFramework.azimuthalaverage_allcyclones(radiusbins,pres[:,:,1,t],cyclones,centers_and_labels,2000))
        windspeedaverage = 0.5(windspeedaverage .+ AvailablePotentialEnergyFramework.azimuthalaverage_allcyclones(radiusbins,windspeed[:,:,:,t],cyclones,centers_and_labels,2000))
    end
    p1 = @views heatmap(permutedims(pres_anomaly[:,:,t]),c=:viridis,title = "Pressure anomaly (hPa)")
    p2 = plot(radiusbins[1:end-1]./1e3,presaverage./100,title="Azimuthally averaged surface pressure anomaly (hPa) ")
    p3 = heatmap(radiusbins[1:end-1]./1e3,1:80,windspeedaverage,title="Azimuthally averaged surface windspeed (m/s) t = $t")
    plot(p1,p2,p3,size=(1200,390))
end
gif(anim, "testresults/azimuthal_averages_alltimesteps.gif", fps = 1)


text:
# using IntervalArithmetic
using DifferentialEquations
using LinearAlgebra

function Flo(du,u,p,t)
    h=BigFloat(".125")
    A11=(-1) *(1+4. *h^2+4. *h*cos(2. *t)).^(-1/2).*sin(2. *t)
    A12=(1+4. *h^2+(1+4. *h^2+4. *h*cos(2. *t)).^(1/2)+2. *h*cos(2. *t).*(2+(1+4. *h^2+ 4. *h*cos(2. *t)).^(1/2))).^(-1) *((-1) *cos(2. *t).*(1+(-8).*h^2+(1+4. *h^2+4. *h*cos(2. *t)).^(1/2))+h*(1+(-4).*h+16. *h^2+2. *(1+ 4. *h^2+4. *h*cos(2. *t)).^(1/2)+(-8).*h*(1+4. *h^2+4. *h*cos(2. * t)).^(1/2)+(-1) *cos(4. *t)))
    A21=(1+4. *h^2+(1+4. *h^2+4. *h*cos(2. *t) ).^(1/2)+2. *h*cos(2. *t).*(2+(1+4. *h^2+4. *h*cos(2. *t)).^(1/2))).^(-1) *(h+4. *h^2. *(1+4. *h)+2. *h*(1+4. *h).*(1+4. *h^2+4. *h*cos(2. *t)).^(1/2)+(-1) *cos(2. *t).*(1+(-8).*h^2+(1+4. *h^2+4. *h*cos( 2. *t)).^(1/2))+(-1) *h*cos(4. *t))
    A22=(1+4. *h^2+4. *h*cos(2. *t)).^( .-1/2).*sin(2. *t)

    du[1]=A11*u[1]+A12*u[2]
    du[2]=A21*u[1]+A22*u[2]
end
#

T_end=BigFloat(pi)
tspan=(BigFloat("0"),T_end)
relT=BigFloat("1e-72")
absT=relT*1e-2

u0=[BigFloat("1"),BigFloat("0")]
prob = ODEProblem(Flo,u0,tspan)
t,u1=solve(prob, Feagin14(),reltol=relT,abstol=absT,maxiters=1e20)
diff=sqrt((u1[1,end]-u1[1,1])^2+(u1[2,end]-u1[2,1])^2)
println(diff)
#
# u0=[BigFloat("0"),BigFloat("1")]
# prob = ODEProblem(Flo,u0,tspan)
# t,u2=solve(prob, RK4(),reltol=relT,abstol=absT,maxiters=1e15)
#
# monodromy_matrix=[u1[1,end] u2[1,end]; u1[2,end] u2[2,end]]
# discriminant=abs(tr(monodromy_matrix))-2
# println(discriminant)
# Flo1=eigvals(monodromy_matrix)[1]
# Flo2=eigvals(monodromy_matrix)[2]


text:
push!(LOAD_PATH, "./src")

using Pkg
Pkg.activate("MLJ_env", shared=true)

using Test
using YaVocabulary

@testset "Basics" begin
  vlst = ["B", "A", "G", "E", "A", "M", "P", "M"]
  vocab = Vocabulary(vlst)

  @test size(vocab) == (length ∘ unique)(vlst)
  @test get_id(vocab, "A") == 2
  @test one_hot_encode(vocab, "B") == [1., 0., 0., 0., 0., 0.]
  @test isnothing(get_id(vocab, "z"))
  @test get_word(vocab, 6) == "P"
  @test isnothing(get_word(vocab, 10))

  add!(vocab, "Y")
  @test size(vocab) == (length ∘ unique)(vlst) + 1
  @test one_hot_encode(vocab, "B") == [1., 0., 0., 0., 0., 0., 0.]
  @test one_hot_encode(vocab, "Y") == [0., 0., 0., 0., 0., 0., 1.]
  @test isnothing(one_hot_encode(vocab, "z"))
end

max_stars_repo_path: Data Science From Scratch/test/ya_vocabulary_test.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
using Test
using Statistics
using JLD2

using Oceananigans
using Oceananigans.Architectures: array_type

using Oceananigans.Fields: location

function generate_some_interesting_simulation_data(Nx, Ny, Nz; architecture=CPU())
    grid = RegularRectilinearGrid(size=(Nx, Ny, Nz), extent=(64, 64, 32))

    T_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(5e-5), bottom = GradientBoundaryCondition(0.01))
    u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(-3e-4))

    @inline Qˢ(x, y, t, S, evaporation_rate) = - evaporation_rate * S
    evaporation_bc = FluxBoundaryCondition(Qˢ, field_dependencies=:S, parameters=3e-7)
    S_bcs = FieldBoundaryConditions(top=evaporation_bc)

    model = NonhydrostaticModel(
               architecture = architecture,
                       grid = grid,
        boundary_conditions = (u=u_bcs, T=T_bcs, S=S_bcs)
    )

    dTdz = 0.01
    Tᵢ(x, y, z) = 20 + dTdz * z + 1e-6 * randn()
    uᵢ(x, y, z) = 1e-3 * randn()
    set!(model, u=uᵢ, w=uᵢ, T=Tᵢ, S=35)

    wizard = TimeStepWizard(cfl=1.0, Δt=10.0, max_change=1.1, max_Δt=1minute)
    simulation = Simulation(model, Δt=wizard, stop_time=2minutes)

    u, v, w = model.velocities

    computed_fields = (
        b = BuoyancyField(model),
        ζ = ComputedField(∂x(v) - ∂y(u)),
        ke = ComputedField(√(u^2 + v^2))
    )

    fields_to_output = merge(model.velocities, model.tracers, computed_fields)

    simulation.output_writers[:jld2_3d_with_halos] =
        JLD2OutputWriter(model, fields_to_output,
                  prefix = "test_3d_output_with_halos",
            field_slicer = FieldSlicer(with_halos=true),
                schedule = TimeInterval(30seconds),
                   force = true)

    profiles = NamedTuple{keys(fields_to_output)}(AveragedField(f, dims=(1, 2)) for f in fields_to_output)

    simulation.output_writers[:jld2_1d_with_halos] =
        JLD2OutputWriter(model, profiles,
                  prefix = "test_1d_output_with_halos",
            field_slicer = FieldSlicer(with_halos=true),
                schedule = TimeInterval(30seconds),
                   force = true)

    run!(simulation)

    return nothing
end

@testset "OutputReaders" begin
    @info "Testing output readers..."

    Nx, Ny, Nz = 16, 10, 5
    generate_some_interesting_simulation_data(Nx, Ny, Nz)
    Nt = 5

    filepath3d = "test_3d_output_with_halos.jld2"
    filepath1d = "test_1d_output_with_halos.jld2"

    for arch in archs
        @testset "FieldTimeSeries{InMemory} [$(typeof(arch))]" begin
            @info "  Testing FieldTimeSeries{InMemory} [$(typeof(arch))]..."

            ## 3D Fields

            u3 = FieldTimeSeries(filepath3d, "u", architecture=arch)
            v3 = FieldTimeSeries(filepath3d, "v", architecture=arch)
            w3 = FieldTimeSeries(filepath3d, "w", architecture=arch)
            T3 = FieldTimeSeries(filepath3d, "T", architecture=arch)
            b3 = FieldTimeSeries(filepath3d, "b", architecture=arch)
            ζ3 = FieldTimeSeries(filepath3d, "ζ", architecture=arch)

            @test location(u3) == (Face, Center, Center)
            @test location(v3) == (Center, Face, Center)
            @test location(w3) == (Center, Center, Face)
            @test location(T3) == (Center, Center, Center)
            @test location(b3) == (Center, Center, Center)
            @test location(ζ3) == (Face, Face, Center)

            @test size(u3) == (Nx, Ny, Nz,   Nt)
            @test size(v3) == (Nx, Ny, Nz,   Nt)
            @test size(w3) == (Nx, Ny, Nz+1, Nt)
            @test size(T3) == (Nx, Ny, Nz,   Nt)
            @test size(b3) == (Nx, Ny, Nz,   Nt)
            @test size(ζ3) == (Nx, Ny, Nz,   Nt)

            ArrayType = array_type(arch)
            for fts in (u3, v3, w3, T3, b3, ζ3)
                @test fts.data.parent isa ArrayType
            end

            if arch isa CPU
                @test u3[1, 2, 3, 4] isa Number
                @test u3[1] isa Field
                @test v3[2] isa Field
            end

            ## 1D AveragedFields

            u1 = FieldTimeSeries(filepath1d, "u", architecture=arch)
            v1 = FieldTimeSeries(filepath1d, "v", architecture=arch)
            w1 = FieldTimeSeries(filepath1d, "w", architecture=arch)
            T1 = FieldTimeSeries(filepath1d, "T", architecture=arch)
            b1 = FieldTimeSeries(filepath1d, "b", architecture=arch)
            ζ1 = FieldTimeSeries(filepath1d, "ζ", architecture=arch)

            @test location(u1) == (Nothing, Nothing, Center)
            @test location(v1) == (Nothing, Nothing, Center)
            @test location(w1) == (Nothing, Nothing, Face)
            @test location(T1) == (Nothing, Nothing, Center)
            @test location(b1) == (Nothing, Nothing, Center)
            @test location(ζ1) == (Nothing, Nothing, Center)

            @test size(u1) == (1, 1, Nz,   Nt)
            @test size(v1) == (1, 1, Nz,   Nt)
            @test size(w1) == (1, 1, Nz+1, Nt)
            @test size(T1) == (1, 1, Nz,   Nt)
            @test size(b1) == (1, 1, Nz,   Nt)
            @test size(ζ1) == (1, 1, Nz,   Nt)

            for fts in (u1, v1, w1, T1, b1, ζ1)
                @test fts.data.parent isa ArrayType
            end

            if arch isa CPU
                @test u1[1, 1, 3, 4] isa Number
                @test u1[1] isa Field
                @test v1[2] isa Field
            end
        end
    end

    for arch in archs
        @testset "FieldTimeSeries{OnDisk} [$(typeof(arch))]" begin
            @info "  Testing FieldTimeSeries{OnDisk} [$(typeof(arch))]..."

            ArrayType = array_type(arch)

            ζ = FieldTimeSeries(filepath3d, "ζ", backend=OnDisk(), architecture=arch)
            @test location(ζ) == (Face, Face, Center)
            @test size(ζ) == (Nx, Ny, Nz, Nt)
            @test ζ[1] isa Field
            @test ζ[2] isa Field
            @test ζ[1].data.parent isa ArrayType

            b = FieldTimeSeries(filepath1d, "b", backend=OnDisk(), architecture=arch)
            @test location(b) == (Nothing, Nothing, Center)
            @test size(b) == (1, 1, Nz, Nt)
            @test b[1] isa Field
            @test b[2] isa Field
        end
    end

    for arch in archs
        @testset "FieldTimeSeries{InMemory} reductions" begin
            @info "  Testing FieldTimeSeries{InMemory} reductions..."

            for name in ("u", "v", "w", "T", "b", "ζ"), fun in (sum, mean, maximum, minimum)
                f = FieldTimeSeries(filepath3d, name, architecture=CPU())

                ε = eps(maximum(abs, f.data.parent))

                val1 = fun(f)
                val2 = fun([fun(f[n]) for n in 1:Nt])

                @test val1 ≈ val2 atol=4ε
            end
        end
    end

    for Backend in [InMemory, OnDisk]
        @testset "FieldDataset{$Backend}" begin
            @info "  Testing FieldDataset{$Backend}..."

            ds = FieldDataset(filepath3d, backend=Backend())

            @test ds isa FieldDataset
            @test length(keys(ds.fields)) == 8
            @test ds["u"] isa FieldTimeSeries
            @test ds["v"][1] isa Field
            @test ds["T"][2] isa Field
        end
    end

    rm("test_3d_output_with_halos.jld2")
    rm("test_1d_output_with_halos.jld2")
end


text:
module FooBuild

using Test
import BasicSandbox

@test 2 == 2
BasicSandbox.greet()

end


text:
using Heuristics, StaticArrays

# Particle Type
N = 10
T = typeof(10.0)
pp = Heuristics.Particle{T,N}(undef)

# Testing types
@test pp.x isa SizedArray{Tuple{N},T,1,1}
@test pp.v isa SizedArray{Tuple{N},T,1,1}
@test pp.p isa SizedArray{Tuple{N},T,1,1}
@test pp.fx == 0.0
@test pp.fp == 0.0
@test pp.fx isa T
@test pp.fp isa T

tempX = rand(N)
tempV = rand(N)
tempP = rand(N)

pp.x .= tempX
pp.v .= tempV
pp.p .= tempP
pp.fx = 5;
pp.fp = 6;

@test pp.fx == 5.0
@test pp.fp == 6.0
@test pp.fx isa T
@test pp.fp isa T
@test pp.x == tempX
@test pp.v == tempV
@test pp.p == tempP
@test pp.x isa SizedArray{Tuple{N},T,1,1}
@test pp.v isa SizedArray{Tuple{N},T,1,1}
@test pp.p isa SizedArray{Tuple{N},T,1,1}

# Testing length(p::Particle)
@test length(pp) == N


text:
export calculate_press_grad, calculate_press_grad_all

"""
    port_dist(port1, port2, port2channels, portdists)
    port_dist(port1, port2, res)

Distance between two pressure ports
"""
port_dist(port1, port2, res::RcRes) = port_dist(port1, port2, res.ex)
function port_dist(port1, port2, ex::RcExp)
    @unpack portdists, port2channels = ex.setup
    port_dist(port1, port2, port2channels, portdists)
end
function port_dist(port1, port2, port2channels, portdists)
    ports = keys(port2channels)
    slot1, slot2 = find(x->x==port1, ports)[1], find(x->x==port2, ports)[1]
    deltaL = sum(portdists[slot1:slot2-1])
    return deltaL
end

"""
    calculate_press_grad(res, port1, port2)

Calculates the pressure gradient over the ports:

"""
calculate_press_grad(res::RcRes, port1, port2) =
    calculate_press_grad(res.ex, res.resl)
function calculate_press_grad(ex::RcExp, resl::ExpLabViewResults, port1, port2)
    @unpack press, press_std = resl
    #press, press_std = resl.press_raw, resl.press_std_raw
    deltaL = port_dist(port1, port2, ex)
    grad, grad_std = calculate_press_grad(press[port1], press[port2], press_std[port1], press_std[port2], deltaL)
    return grad, grad_std, (port1=>port2, deltaL)
end
function calculate_press_grad(press1, press2, press1_std, press2_std, deltaL)
    grad = (press1-press2)./deltaL
    grad_std = sqrt.(press1_std.^2 + press2_std.^2)./deltaL
    return grad, grad_std
end
calculate_press_grad(res::RcExp) = calculate_press_grad(res.ex, res.resl)
function calculate_press_grad(ex::RcExp, resl::ExpLabViewResults)
    @unpack port2channels = ex.setup
    # select ice1, ice2 or in4
    nums = [8,9,4:-1:1;]
    port1 = 0
    for n in nums
        if haskey(port2channels, Port.PressPorts(n))
            port1 = Port.PressPorts(n)
            break
        end
    end
    # select ice4, ice3 or out1
    nums = [11,10,4:7;]
    port2 = 0
    for n in nums
        if haskey(port2channels, Port.PressPorts(n))
            port2 = Port.PressPorts(n)
            break
        end
    end
    calculate_press_grad(res, port1, port2)
end

calculate_press_grad_all(res::RcRes) = calculate_press_grad_all(res.ex, res.resl)
function calculate_press_grad_all(ex::RcExp, resl::ExpLabViewResults,
                                  ports1 = collect(keys(ex.setup.port2channels))[1:end-1],
                                  ports2 = collect(keys(ex.setup.port2channels))[2:end])
    grads = Float64[]
    grads_std = Float64[]
    ports = []
    for (p1,p2) in zip(ports1, ports2)
        pd, spd, ps = calculate_press_grad(ex, resl, p1, p2)
        append!(grads, pd)
        append!(grads_std, spd)
        push!(ports, ps)
    end
    nl = length(ports1)
    return reshape(grads, length(grads)÷nl, nl),
           reshape(grads_std, length(grads_std)÷nl, nl),
           ports
end

calc_f(grad_p, D, Q) =  (grad_p * pi^2 * D^5)/(8*Q^2*rho)
function calc_f(res::RcRes, port1, port2, inds=nothing)
    @unpack resl, resi = res
    @unpack trun, irun = resl
    if inds!=nothing
        irun = irun[inds]
        trun = trun[inds]
    end
    grad_p, grad_p_std, ports = calculate_press_grad(res, port1, port2)
    grad_p, grad_p_std = grad_p[irun], grad_p_std[irun]
    Q = resl.Q[irun]
    D, D_std = get_time_series(resi, trun)
    calc_f.(grad_p, D, Q)
end

function Re_nr(D,Q)
    A = pi/4*D^2
    D*Q/(nu*A)
end
function Re_nr(res::RcRes, inds=nothing)
    @unpack resl, resi = res
    @unpack trun, irun = resl
    if inds!=nothing
        irun = irun[inds]
        trun = trun[inds]
    end
    D, D_std = get_time_series(resi, trun)
    Q = resl.Q[irun]
    Re_nr.(D,Q)
end

"""
    f_haaland(Re, epsilon, Dh)
    f_haaland(Re, rel_roughness)

D-W f factor by Haaland 1983.  Epsilon is the size of the roughness,
Dh hyd-diameter.
"""
f_haaland(Re, epsilon, Dh) = f_haaland(Re, epsilon/Dh)
function f_haaland(Re, rel_roughness)
    tmp = -1.8 * log( ((rel_roughness/3.7)^1.11 + 6.9/Re) )
    1/tmp^2
end

### stacking
function stack_it(vec, step)
    halfstep = step÷2
    out = eltype(vec)[]
    for i = halfstep+1:step:length(vec)-halfstep
        push!(out, mean(vec[i-halfstep:i+halfstep]))
    end
    out
end
function stack_it_time(time, step)
    halfstep = step÷2
    out = eltype(time)[]
    for i = halfstep+1:step:length(time)-halfstep
        push!(out, time[i])
    end
    out
end

### purge high std
"""
    find_low_std(press_std, ports, irun=1:length(first(press_std)), thresh=100)
    -> inds

Find indices of low std.
"""
function find_low_std(press_std, ports, irun=1:length(first(press_std)), thresh=100)
    out = press_std[ports[1]][irun].<thresh
    for i=2:length(ports)
        out = (press_std[ports[i]][irun].<thresh) .& out
    end
    find(out)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ fdb00d02-cb5b-49e0-aa5b-af108ab078f4
begin
	using PlutoUI
end

# ╔═╡ deb978d1-54a6-4d1d-9722-0a44778c4a94
md"""
### Discrete State Models for Transcriptional Regulation

A few lectures ago, we introduced balance equations for messenger RNA (mRNA) and protein abundance.
The material balance equations governing the abundance of the ith mRNA species ($m_{i}$) and the ith protein species ($p_{i}$) are given by (in cellmass specific units in a growing population):

$$\begin{eqnarray}
\frac{dm_{i}}{dt} & = & r_{X,i}u_{i} - (\theta_{m,i}+\mu)m_{i}+\lambda_{i}\\
\frac{dp_{i}}{dt} & = & r_{L,i}w_{i} - (\theta_{p,i}+\mu)p_{i}
\end{eqnarray}$$

The $r_{X, i}$ denote the _kinetic limit_ of transcription (if a gene were `ON`, this is the rate of production that you would expect), while the $r_{L, i}$ denote the _kinetic limit_ of translation. The second set of terms represents the degradation and dilution due to growth. The term $\lambda_{i}$ denotes the rate of unregulated or background expression (units: conc/time). Finally, the $u_{i}\in\left[0,1\right]$ and $w_{i}\in\left[0,1\right]$ terms describe the input controlling the transcription of gene $i$, and the translation of mRNA $i$, respectively. There is a __huge__ variety of approaches to formulate these control functions (many are super interesting)!

In this lecture, we'll take a deeper dive into the modeling and analysis of transcription (TX) and translation (TL). In particular, we will:

* Introduce the discrete state promoter model formulation for the $u\left(\dots\right)$ function
* Explore gene expression dynamics in a population of single _Escherichia coli_  cells
* Revisit the mRNA and protein balance equations to reveal their shameful dark secret

"""

# ╔═╡ 7bef424f-8e94-4bdd-8c13-09a08d54ab7b
md"""
### Some fun: Cybernetic models
Cybernetic models, developed over the last 30 years by [Ramkrishna and coworkers](https://engineering.purdue.edu/ChE/people/ptProfile?resource_id=11806), use ideas from Microeconomics (allocation of a scarce resource) to predict global gene expression:

* [Kompala DS, Ramkrishna D, Jansen NB, Tsao GT. Investigation of bacterial growth on mixed substrates: experimental evaluation of cybernetic models. Biotechnol Bioeng. 1986 Jul;28(7):1044-55. doi: 10.1002/bit.260280715. PMID: 18555426.](https://pubmed.ncbi.nlm.nih.gov/18555426/)
* [Effective Dynamic Models of Metabolic Networks Michael Vilkhovoy, Mason Minot, Jeffrey D. Varner bioRxiv 047316; doi: https://doi.org/10.1101/047316](https://www.biorxiv.org/content/10.1101/047316v2)
"""

# ╔═╡ 9aff53c6-0d8e-431f-b690-72c9b97be740
md"""
### Discrete state model for promoter functions

Suppose a promoter $P$ can exits in one of $s=1,2\dots,\mathcal{S}$ possible _disjoint_ microstates, where each microstate $s$ has some pseudo energy $\epsilon_{s}$. Some microstates will lead to expression (the ability to produce an mRNA molecule by transcription), while others do not. For each microstate $s$, let's assign a pseudo energy $\epsilon_{s}$, where by definition $\epsilon_{1}=0$; we assume the base state has zero energy. Next, suppose the probability that promoter $P$ is in microstate $s$ follows a [Boltzmann distribution](https://en.wikipedia.org/wiki/Boltzmann_distribution) which says:

$$p_{i} = \frac{1}{Z} \times f_{i}\exp\left(-\beta\epsilon_{i}\right)\qquad{i=1,2,\dots,\mathcal{S}}$$

where $p_{i}$ denotes the probability that promoter $P$ is in microstate $i=1,2,\dots,\mathcal{S}$, $f_{i}$ denotes a state-specific factor $f_{i}\in\left[0,1\right]$, $\beta$ denotes the [thermodynamic beta](https://en.wikipedia.org/wiki/Thermodynamic_beta) and $Z$ denotes a normalization factor (called the [Partiton function](https://en.wikipedia.org/wiki/Partition_function_(statistical_mechanics)) in the statistical physics community). We can find $Z$ using the summation law of discrete probolity e.g.,  $\sum_{s}p_{s} = 1$ which gives:

$$Z = \sum_{s=1}^{\mathcal{S}}f_{i}\exp\left(-\beta\epsilon_{i}\right)$$

which gives:

$$p_{i} = \frac{f_{i}\exp\left(-\beta\epsilon_{i}\right)}{\displaystyle \sum_{s=1}^{\mathcal{S}}f_{i}\exp\left(-\beta\epsilon_{i}\right)}\qquad{i=1,2,\dots,\mathcal{S}}$$.

Finally, we relate the probability that promoter $P$ is in microstate $s$ back to the $\bar{u}\left(\dots\right)$ control function by computing the overall probability that the desired event happens, e.g., promoter $P$ undergoes transcription. We know if $\Omega = \left\{1,2,\dots,\mathcal{S}\right\}$, then we can define the subset $\mathcal{A}\subseteq\Omega$ in which the desired event happens (in this case transcription). Given $\mathcal{A}$, the $\bar{u}\left(\dots\right)$ function becomes:

$$\bar{u}=\sum_{s\in{\mathcal{A}}}p_{s}$$

where the control function $\bar{u}$ is the sum of a regulated $u\left(\dots\right)$ and unregulated $u^{\dagger}\left(\dots\right)$ component:

$$\bar{u} = u + u^{\dagger}$$

##### Example

* [Moon TS, Lou C, Tamsir A, Stanton BC, Voigt CA. Genetic programs constructed from layered logic gates in single cells. Nature. 2012;491(7423):249-253. doi:10.1038/nature11516](https://www.ncbi.nlm.nih.gov/labs/pmc/articles/PMC3904217/)
"""

# ╔═╡ 8494b59a-bcb4-4836-a3dd-5a72d4060d7b
md"""
### A microscopic view of gene expression

Typical transcription and translation models (such as those we'll use in this class) are population-averaged descriptions of gene expression dynamics. However, is this what is happening inside a single cell? For example, is there significant cell to cell variation in the _Escherichia coli_ cells we are using to make the product `XYZ`? Yes, there is. However, can we still use population-averaged models? It depends.

##### Example:

* [Golding I, Paulsson J, Zawilski SM, Cox EC. Real-time kinetics of gene activity in individual bacteria. Cell. 2005 Dec 16;123(6):1025-36. doi: 10.1016/j.cell.2005.09.031. PMID: 16360033.](https://pubmed.ncbi.nlm.nih.gov/16360033/)
"""

# ╔═╡ ffddbade-a162-4ebc-a81c-3ff672ae3a2c
md"""

### Revisit: mRNA and protein balance equations

The mRNA and protein balance equations shown above hold a dark secret; they assume a well-mixed population of cells. The work of Golding et al. suggests that the well-mixed assumption is not valid. Does it matter?

##### Derivation
Let's see where the well-mixed assumption comes into the derivation of the mRNA and protein balance equations. Consider the integral form of the mol balance around the ith mRNA $m_{i}$ written with respect to some volume basis $\beta$ in a _batch_ system (no convective or diffusive transport terms):

$$\frac{d}{dt}\left(\int_{\Omega}m_{i}d\beta\right) = \int_{\Omega}\left(
\sum_{j=1}^{\mathcal{R}}\sigma_{ij}\hat{r}_{j}\right)d\beta\qquad{i=1,2,\dots,\mathcal{M}}$$

The term on the left-hand side is the accumulation term (in the domain $\Omega$). In contrast, the right-hand side is the reaction term (in the domain $\Omega$) where $\sigma_{ij}$ denotes the stoichiometric coefficient for species $i$ in reaction $j$, $\mathcal{R}$ represents the number of reactions, and $\hat{r}_{j}$ denotes the kinetics of reaction $j$ per unit $\beta$.

###### Well mixed assumption
Let's assume the system has no $\beta$ dependence, i.e., the system is well mixed such that we can pull $m_{i}$ out of the integral:

$$\frac{d}{dt}\left(\int_{\Omega}m_{i}d\beta\right) \simeq \frac{d}{dt}\left(\langle m_{i} \rangle\beta\right)$$

Similarly, we can pull the out the reaction terms to give:

$$\int_{\Omega}\left(
\sum_{j=1}^{\mathcal{R}}\sigma_{ij}\hat{r}_{j}\right)d\beta \simeq
\langle \sum_{j=1}^{\mathcal{R}}\sigma_{ij}\hat{r}_{j} \rangle\beta$$

where $\langle\cdot\rangle$ denotes an averaged quantity over the domain $\Omega$. Thus, the mRNA balance equation is given by:

$$\frac{d}{dt}\left(\langle m_{i} \rangle\beta\right) = \langle \sum_{j=1}^{\mathcal{R}}\sigma_{ij}\hat{r}_{j} \rangle\beta\qquad{i=1,2,\dots,\mathcal{M}}$$

Expanding the accumulation term and simplifying gives:

$$\frac{d\langle m_{i}\rangle}{dt} = \langle \sum_{j=1}^{\mathcal{R}}\sigma_{ij}\hat{r}_{j} \rangle -
\frac{\langle m_{i} \rangle}{\beta}\frac{d\beta}{dt}\qquad{i=1,2,\dots,\mathcal{M}}$$

When we use cell mass specific units, we assume $\beta = XV_{R}$, where $X$ denotes the cell mass abundance (units: gDW/L) in the culture, and $V_{R}$ (units: L) represents the culture volume. Thus, for a _batch_ culture of cells we have:

$$\frac{d\langle m_{i}\rangle}{dt} = \langle \sum_{j=1}^{\mathcal{R}}\sigma_{ij}\hat{r}_{j} \rangle -
\mu\langle m_{i} \rangle\qquad{i=1,2,\dots,\mathcal{M}}$$

where $\mu$ denotes the specific growth rate (units: hr$^{-1}$). Finally, for mRNA we know that we have both regulated and unregulated production of transcript,  and a transcript degradation term, thus:

$$\frac{dm_{i}}{dt} = r_{X,i}u_{i} - (\theta_{m,i}+\mu)m_{i}+\lambda_{i}\qquad{i=1,2,\dots,\mathcal{M}}$$

where we have dropped the $\langle \cdot \rangle$ notation and substituted our TX and degradation reactions.
"""

# ╔═╡ ba43f178-01c8-4801-b750-32b84b923707
md"""
### Summary and conclusions

In this lecture we:

* Introduced discrete state promoter models used them to capture averaged expression dynamics for logic circuits
* Explored gene expression dynamics in a population of single _Escherichia coli_ cells
* Derived the population-averaged mRNA and protein balance equations in cell mass specific units
"""

# ╔═╡ f25cda13-83fd-470a-b755-ec7b4877b19d
md"""
### Next time

During the next lecture, we'll discuss:
* numerical and experimental examples of deterministic and stochastic gene expression
* models of kinetics of transcription and translation

"""

# ╔═╡ 9a85e011-6188-41f0-b309-b0a869ec5e9d
TableOfContents(title="📚 Table of Contents", indent=true, depth=5, aside=true)

# ╔═╡ f5cafa0c-4173-4354-a873-0953e9bcfcb1
html"""
<script>
	// initialize -
	var section = 0;
	var subsection = 0;
	var subsubsection = 0;
	var headers = document.querySelectorAll('h3, h5, h6');

	// main loop -
	for (var i=0; i < headers.length; i++) {

		var header = headers[i];
	    var text = header.innerText;
	    var original = header.getAttribute("text-original");
	    if (original === null) {

			// Save original header text
	        header.setAttribute("text-original", text);
	    } else {

			// Replace with original text before adding section number
	        text = header.getAttribute("text-original");
	    }

	    var numbering = "";
	    switch (header.tagName) {
	        case 'H3':
	            section += 1;
	            numbering = section + ".";
	            subsection = 0;
	            break;
	        case 'H5':
	            subsection += 1;
	            numbering = section + "." + subsection;
	            break;
			case 'H6':
	            subsubsection += 1;
	            numbering = section + "." + subsection + "." + subsubsection;
	            break;
	    }
		// update the header text
		header.innerText = numbering + " " + text;
	};
</script>"""

# ╔═╡ 248cbb68-9ec6-11ec-05b5-abe80f6c2b53
html"""
<style>
main {
    max-width: 860px;
    width: 70%;
    margin: auto;
    font-family: "Roboto, monospace";
}

a {
    color: blue;
    text-decoration: none;
}
</style>"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.37"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.3"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "bf0a1121af131d9974241ba53f601211e9303a9e"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.37"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─deb978d1-54a6-4d1d-9722-0a44778c4a94
# ╟─7bef424f-8e94-4bdd-8c13-09a08d54ab7b
# ╟─9aff53c6-0d8e-431f-b690-72c9b97be740
# ╟─8494b59a-bcb4-4836-a3dd-5a72d4060d7b
# ╟─ffddbade-a162-4ebc-a81c-3ff672ae3a2c
# ╟─ba43f178-01c8-4801-b750-32b84b923707
# ╟─f25cda13-83fd-470a-b755-ec7b4877b19d
# ╠═9a85e011-6188-41f0-b309-b0a869ec5e9d
# ╠═fdb00d02-cb5b-49e0-aa5b-af108ab078f4
# ╟─f5cafa0c-4173-4354-a873-0953e9bcfcb1
# ╟─248cbb68-9ec6-11ec-05b5-abe80f6c2b53
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
@testset "c api memory leak test" begin # Issue https://github.com/JuliaIO/VideoIO.jl/issues/246

    if(Sys.islinux())  # TODO: find a method to get cross platform memory usage, see: https://discourse.julialang.org/t/how-to-get-current-julia-process-memory-usage/41734/4

        function get_memory_usage()
            open("/proc/$(getpid())/statm") do io
                split(read(io, String))[1]
            end
        end

        file = joinpath(videodir, "annie_oakley.ogg")

        @testset "open file test" begin
            check_size = 10
            usage_vec = Vector{String}(undef, check_size)

            for i in 1:check_size

                f = VideoIO.openvideo(file)
                close(f)
                GC.gc()

                usage_vec[i] = get_memory_usage()
            end

            println(usage_vec)

            @test usage_vec[end-1] == usage_vec[end]
        end

        @testset "open and read file test" begin
            check_size = 10
            usage_vec = Vector{String}(undef, check_size)

            for i in 1:check_size

                f = VideoIO.openvideo(file)
                img = read(f)
                close(f)
                GC.gc()

                usage_vec[i] = get_memory_usage()
            end

            println(usage_vec)

            @test usage_vec[end-1] == usage_vec[end]
        end
    end
end


text:
#!/usr/bin/env julia
include("Common.jl")
using Common
using RealNeuralNetworks.SWCs
using RealNeuralNetworks.Neurons
using HDF5
#using Plots
# use pyplot backend for vectorized output
#Plots.pyplot()

function main()
    args = parse_commandline()
    numSegmentingPointsList = Vector{Int}()
    totalPathLengthList = Vector{Float64}()
    @sync begin
        for fileName in readdir( args["swcdir"] )
            if contains(fileName, ".swc")
                @async begin
                    @time swc = SWCs.load( joinpath( args["swcdir"], fileName ) )
                    totalPathLength = SWCs.get_total_length( swc )
                    println("total path length of cell $(fileName): " *
                            "$(totalPathLength/1000) micron")
                    push!(totalPathLengthList, totalPathLength)

                    neuron = Neuron(swc)
                    numSegmentingPoints = Neurons.get_num_branching_points(neuron)
                    println("number of branching points of cell $(fileName): $numSegmentingPoints")
                    push!(numSegmentingPointsList, numSegmentingPoints)
                end
            end
        end
    end
    h5write("statistics.h5",    "totalPathLengthList",      totalPathLengthList)
    h5write("statistics.h5",    "numSegmentingPointsList",   numSegmentingPointsList)
    #p = histogram( lengthList )
    # display(p)
end

@time main()