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

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

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

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

# acetotrophic methanogens methanol
df_amm = df_am[(df_am.methanol_methylCoM_or.==1) , :]
Genome_size_amm = convert(Array{Float64}, df_amm.bp)
V_cell_amm = DEBmicroTrait.genome_size_to_cell_volume(Genome_size_amm)
Min_gen_time_amm = df_amm.minimum_generation_time_hours
rrn_copies_amm = DEBmicroTrait.genome_size_to_rRNA_copy_number(Genome_size_amm)
Gram_stain_amm = repeat(["+"], size(V_cell_amm,1))
SAV_amm = DEBmicroTrait.surface_area_volume_ratio(V_cell_amm./1e-18)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

################################################################################
                        # Acetotrophic methanogens - methanol
N_C = [1]
y_DE = [0.44]./N_C
y_EM = ones(size(V_cell_am,1))
ρ_ps_amm= zeros(size(V_cell_am,1))
r_msm = zeros(size(V_cell_am,1))
for i in 1:size(V_cell_am,1)
    find_ρ(x) = DEBmicroTrait.constrain_transporter_density(x, [V_cell_am[i]], [Min_gen_time_am[i]], [Gram_stain_am[i]], [rrn_copies_am[i]], y_DE[1], N_C[1], [y_EM[i]])
    ρ_p = Roots.find_zero(find_ρ, 1.0)
    ρ_ps_amm[i] = ρ_p[1]
    r_msm[i] = DEBmicroTrait.check_growth_rate(ρ_p[1], [V_cell_am[i]], [Min_gen_time_am[i]], [Gram_stain_am[i]], [rrn_copies_am[i]], y_DE[1], N_C[1], [y_EM[i]])
end
# respiration
y_EV = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies_am)
r_growth = @. (1-1/y_EV)*r_msm
r_main   = DEBmicroTrait.cell_volume_to_specific_maintenance_rate(V_cell_am, Min_gen_time_am, Gram_stain_am)
N_SB_amm = [DEBmicroTrait.transporter_density_to_monomer_uptake_sites([V_cell_am[i]], ρ_ps_amm[i]*ones(1,1), [Min_gen_time_am[i]], [Gram_stain_am[i]])[1] for i in 1:size(V_cell_am,1)]
k2p = 180.0*60^2
r_assim = @. y_DE*N_C*k2p*N_SB_amm
r_resp_amm = r_growth .+ r_main .+r_assim
median(r_resp_amm)
std(r_resp_amm)
# half saturation constants
Molecular_weight = [32.04]
D_S = DEBmicroTrait.aqueous_diffusivity(Molecular_weight)
K_D_amm = DEBmicroTrait.specific_reference_affinity(V_cell_am, reshape(ρ_ps_amm*100, 65,1), D_S)
median(K_D_amm.*N_C*12.011)
std(K_D_amm.*N_C*12.011)
K_D_amm = K_D_amm.*N_C*12.011
################################################################################
r_resp_aceto = vcat(r_resp_am, r_resp_amm)
K_D_aceto = vcat(K_D_am, K_D_amm)

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

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

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

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

maximum(K_D_mo)
################################################################################
                        # aerobic diazotroph
N_C = [6]
y_DE = [2.61]./N_C
y_EM = ones(size(V_cell_aerob_diazo,1))

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

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

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

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

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

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

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

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

################################################################################
                        # Acetate nitrite reductase
N_C = [2]
y_DE = [0.38]./N_C
y_EM = ones(size(V_cell_nitrite_red,1))

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

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

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

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

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

df_out = DataFrame()

df_out.functionalgroup = repeat(["Acetotrophic methanogen"],size(df_am.genome,1))
df_out.rrespg =

df_out.genome = vcat(df_aerob_hetero.genome, df_aerob_hetero.genome, df_ferm.genome, df_am.genome, df_am.genome, df_hmo.genome, df_hmo.genome, df_mo.genome)
df_out.habitat = vcat(df_aerob_hetero.habitat_binnedfrom, df_aerob_hetero.habitat_binnedfrom, df_ferm.habitat_binnedfrom, df_am.habitat_binnedfrom,  df_am.habitat_binnedfrom, df_hmo.habitat_binnedfrom, df_hmo.habitat_binnedfrom, df_mo.habitat_binnedfrom)
df_out.rresp = vcat(r_resp_aerob_hetero, r_resp_ferm, r_resp_aceto, r_resp_hydro, r_resp_mo)
df_out.K = vcat(K_D_aerob_hetero[:,1], K_D_ferm[:,1], K_D_aceto[:,1], K_D_hydro[:,1], K_D_mo[:,1])

functional_groups = vcat(repeat(["Aerobic heterotroph glucose"],size(df_aerob_hetero.genome,1)), repeat(["Aerobic heterotroph acetate"],size(df_aerob_hetero.genome,1)), repeat(["Fermenter"],size(df_ferm.genome,1)), repeat(["Acetotrophic methanogen"],size(df_am.genome,1)), repeat(["Acetotrophic methanogen"],size(df_am.genome,1)),
                        repeat(["Hydrogenotrophic methanogen"],size(df_hmo.genome,1)), repeat(["Hydrogenotrophic methanogen"],size(df_hmo.genome,1)),repeat(["Methanotroph"],size(df_mo.genome,1)))

df_out.functionalgroup = functional_groups

CSV.write("/Users/glmarschmann/Data/Zhen/TraitPredictionsIsogenieGenomes_GA.csv", df_out)
################################################################################

################################################################################
                              # Literature comparison
df_lit     = dropmissing(CSV.read("/Users/glmarschmann/Data/Zhen/LiteratureTraits.csv", DataFrame, missingstring="NaN"))

df_lit_aerob_hetero = filter(x->(x.functionalgroup.=="Aerobic heterotroph"), df_lit)
df_lit_ferm = filter(x->(x.functionalgroup.=="Fermenter"), df_lit)
df_lit_aceto = filter(x->(x.functionalgroup.=="Acetotrophic methanogen"), df_lit)
df_lit_hydro = filter(x->(x.functionalgroup.=="Hydrogenotrophic methanogen"), df_lit)
df_lit_meth = filter(x->(x.functionalgroup.=="Methanotroph"), df_lit)

df_pred_aerob_hetero = filter(x->(x.functionalgroup.=="Aerobic heterotroph"), df_out)
df_pred_ferm = filter(x->(x.functionalgroup.=="Fermenter"), df_out)
df_pred_aceto = filter(x->(x.functionalgroup.=="Acetotrophic methanogen"), df_out)
df_pred_hydro = filter(x->(x.functionalgroup.=="Hydrogenotrophic methanogen"), df_out)
df_pred_meth = filter(x->(x.functionalgroup.=="Methanotroph"), df_out)

HypothesisTests.KruskalWallisTest(df_pred_aerob_hetero.rresp, df_lit_aerob_hetero.rresp)
HypothesisTests.KruskalWallisTest(df_pred_aerob_hetero.K, df_lit_aerob_hetero.K)
HypothesisTests.KruskalWallisTest(df_pred_ferm.rresp, df_lit_ferm.rresp)
HypothesisTests.KruskalWallisTest(df_pred_ferm.K, df_lit_ferm.K)
HypothesisTests.KruskalWallisTest(df_pred_aceto.rresp, df_lit_aceto.rresp)
HypothesisTests.KruskalWallisTest(df_pred_aceto.K, df_lit_aceto.K)
HypothesisTests.KruskalWallisTest(df_pred_hydro.rresp, df_lit_hydro.rresp)
HypothesisTests.KruskalWallisTest(df_pred_hydro.K, df_lit_hydro.K)
HypothesisTests.KruskalWallisTest(df_pred_meth.rresp, df_lit_meth.rresp)
HypothesisTests.KruskalWallisTest(df_pred_meth.K, df_lit_meth.K)

(median(df_pred_aerob_hetero.rresp) - 0.125)/0.125
(median(df_pred_ferm.rresp) - 0.125)/0.125
(median(df_pred_aceto.rresp) - 0.125)/0.125
(median(df_pred_hydro.rresp) - 0.125)/0.125
(median(df_pred_meth.rresp) - 0.375)/0.375

(median(df_pred_aerob_hetero.K) - 12)/12
(median(df_pred_ferm.K) - 12)/12
(median(df_pred_aceto.K) - 12)/12
(median(df_pred_hydro.K) - 0.01)/0.01
(median(df_pred_meth.K) - 0.0012)/0.0012

df_fen_aerob_hetero = filter(x->(x.functionalgroup.=="Aerobic heterotroph" && x.habitat.=="fen"), df_out)
df_fen_ferm = filter(x->(x.functionalgroup.=="Fermenter"  && x.habitat.=="fen"), df_out)
df_fen_aceto = filter(x->(x.functionalgroup.=="Acetotrophic methanogen"  && x.habitat.=="fen"), df_out)
df_fen_hydro = filter(x->(x.functionalgroup.=="Hydrogenotrophic methanogen"  && x.habitat.=="fen"), df_out)
df_fen_meth = filter(x->(x.functionalgroup.=="Methanotroph"  && x.habitat.=="fen"), df_out)

df_bog_aerob_hetero = filter(x->(x.functionalgroup.=="Aerobic heterotroph" && x.habitat.=="bog"), df_out)
df_bog_ferm = filter(x->(x.functionalgroup.=="Fermenter"  && x.habitat.=="bog"), df_out)
df_bog_aceto = filter(x->(x.functionalgroup.=="Acetotrophic methanogen"  && x.habitat.=="bog"), df_out)
df_bog_hydro = filter(x->(x.functionalgroup.=="Hydrogenotrophic methanogen"  && x.habitat.=="bog"), df_out)
df_bog_meth = filter(x->(x.functionalgroup.=="Methanotroph"  && x.habitat.=="bog"), df_out)

df_palsa_aerob_hetero = filter(x->(x.functionalgroup.=="Aerobic heterotroph" && x.habitat.=="palsa"), df_out)
df_palsa_ferm = filter(x->(x.functionalgroup.=="Fermenter"  && x.habitat.=="palsa"), df_out)
df_palsa_aceto = filter(x->(x.functionalgroup.=="Acetotrophic methanogen"  && x.habitat.=="palsa"), df_out)
df_palsa_hydro = filter(x->(x.functionalgroup.=="Hydrogenotrophic methanogen"  && x.habitat.=="palsa"), df_out)
df_palsa_meth = filter(x->(x.functionalgroup.=="Methanotroph"  && x.habitat.=="palsa"), df_out)

HypothesisTests.KruskalWallisTest(df_fen_aerob_hetero.rresp, df_bog_aerob_hetero.rresp)
HypothesisTests.KruskalWallisTest(df_fen_aerob_hetero.rresp, df_palsa_aerob_hetero.rresp)
HypothesisTests.KruskalWallisTest(df_palsa_aerob_hetero.rresp, df_bog_aerob_hetero.rresp)

HypothesisTests.KruskalWallisTest(df_fen_ferm.rresp, df_bog_ferm.rresp)
HypothesisTests.KruskalWallisTest(df_fen_ferm.rresp, df_palsa_ferm.rresp)
HypothesisTests.KruskalWallisTest(df_palsa_ferm.rresp, df_bog_ferm.rresp)

HypothesisTests.KruskalWallisTest(df_fen_aceto.rresp, df_bog_aceto.rresp)
HypothesisTests.KruskalWallisTest(df_fen_aceto.rresp, df_palsa_aceto.rresp)
HypothesisTests.KruskalWallisTest(df_palsa_aceto.rresp, df_bog_aceto.rresp)

HypothesisTests.KruskalWallisTest(df_fen_hydro.rresp, df_bog_hydro.rresp)
HypothesisTests.KruskalWallisTest(df_fen_hydro.rresp, df_palsa_hydro.rresp)
HypothesisTests.KruskalWallisTest(df_palsa_hydro.rresp, df_bog_hydro.rresp)

HypothesisTests.KruskalWallisTest(df_fen_meth.rresp, df_bog_meth.rresp)
HypothesisTests.KruskalWallisTest(df_fen_meth.rresp, df_palsa_meth.rresp)
HypothesisTests.KruskalWallisTest(df_palsa_meth.rresp, df_bog_meth.rresp)

HypothesisTests.KruskalWallisTest(df_fen_aerob_hetero.K, df_bog_aerob_hetero.K)
HypothesisTests.KruskalWallisTest(df_fen_aerob_hetero.K, df_palsa_aerob_hetero.K)
HypothesisTests.KruskalWallisTest(df_palsa_aerob_hetero.K, df_bog_aerob_hetero.K)

HypothesisTests.KruskalWallisTest(df_fen_ferm.K, df_bog_ferm.K)
HypothesisTests.KruskalWallisTest(df_fen_ferm.K, df_palsa_ferm.K)
HypothesisTests.KruskalWallisTest(df_palsa_ferm.K, df_bog_ferm.K)

HypothesisTests.KruskalWallisTest(df_fen_aceto.K, df_bog_aceto.K)
HypothesisTests.KruskalWallisTest(df_fen_aceto.K, df_palsa_aceto.K)
HypothesisTests.KruskalWallisTest(df_palsa_aceto.K, df_bog_aceto.K)

HypothesisTests.KruskalWallisTest(df_fen_hydro.K, df_bog_hydro.K)
HypothesisTests.KruskalWallisTest(df_fen_hydro.K, df_palsa_hydro.K)
HypothesisTests.KruskalWallisTest(df_palsa_hydro.K, df_bog_hydro.K)

HypothesisTests.KruskalWallisTest(df_fen_meth.K, df_bog_meth.K)
HypothesisTests.KruskalWallisTest(df_fen_meth.K, df_palsa_meth.K)
HypothesisTests.KruskalWallisTest(df_palsa_meth.K, df_bog_meth.K)

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
module TestSetfield

using LyceumBase: @set!
include("preamble.jl")

struct Imm
    x
    y
end

mutable struct Mut
    a
    b
end

@testset "nested mutable" begin
    let x = orig = Mut(Imm(Mut(1, 2), 3), 4)
        @set! x.a.x.a = 10
        @test x.a.x.a == 10
        @test x === orig
    end
    let x = Imm(Imm(Imm(1, 2), 3), 4)
        @test_throws ErrorException @set! x.x.x.x = 10
    end
end

@testset "arrays" begin
    let v = [1, 2, Mut(1, 2)], x = orig = Mut((x = v, y = 3), 4)
        @set! x.a.x[3].a = 10
        @test x.a.x[3].a == 10
        @test x === orig
    end
    let v = (1, 2, Imm(1, 2)), x = orig = Imm((x = v, y = 3), 4)
        @test_throws ErrorException @set! x.x.x[3].x = 10
    end
end

@testset "index lenses" begin
    @testset "ConstIndexLens" begin
        let v = [1, 2, 3]
            i = j = 1
            @set! v[$(i + j)] = 20
            @test v[2] == 20
        end
        let v = (1, 2, 3)
            i = j = 1
            @test_throws ErrorException @set! v[$(i + j)] = 20
        end
    end
    @testset "DynamicIndexLens" begin
        let v = [1, 2, 3]
            @set! v[end] = 30
            @test v[3] == 30
        end
        let v = (1, 2, 3)
            @test_throws ErrorException @set! v[end] = 30
        end
    end
end

end  # module

ext: jl


text:
# avoid creating a progress bar with @withprogress if progress logging is disabled
# and add a custom progress logger if the current logger does not seem to be able to handle
# progress logs
macro ifwithprogresslogger(progress, exprs...)
    return quote
        if $progress
            if $hasprogresslevel($Logging.current_logger())
                $ProgressLogging.@withprogress $(exprs...)
            else
                $with_progresslogger($Base.@__MODULE__, $Logging.current_logger()) do
                    $ProgressLogging.@withprogress $(exprs...)
                end
            end
        else
            $(exprs[end])
        end
    end |> esc
end

# improved checks?
function hasprogresslevel(logger)
    return Logging.min_enabled_level(logger) ≤ ProgressLogging.ProgressLevel
end

# filter better, e.g., according to group?
function with_progresslogger(f, _module, logger)
    logger1 = LoggingExtras.EarlyFilteredLogger(progresslogger()) do log
        log._module === _module && log.level == ProgressLogging.ProgressLevel
    end
    logger2 = LoggingExtras.EarlyFilteredLogger(logger) do log
        log._module !== _module || log.level != ProgressLogging.ProgressLevel
    end

    Logging.with_logger(f, LoggingExtras.TeeLogger(logger1, logger2))
end

function progresslogger()
    # detect if code is running under IJulia since TerminalLogger does not work with IJulia
    # https://github.com/JuliaLang/IJulia.jl#detecting-that-code-is-running-under-ijulia
    if (Sys.iswindows() && VERSION < v"1.5.3") || (isdefined(Main, :IJulia) && Main.IJulia.inited)
        return ConsoleProgressMonitor.ProgressLogger()
    else
        return TerminalLoggers.TerminalLogger()
    end
end


text:
using Documenter, Literate
using CUDATutorials

import CUDATutorials.make_dir

function generate_docs()
    outpath = normpath(joinpath(pkgdir(CUDATutorials), "docs", "src"))
    make_dir(outpath)
    tutorial_path = normpath(joinpath(pkgdir(CUDATutorials), "src", "tutorials"))
    for (root, dirs, files) in walkdir(tutorial_path)
        for dir in dirs
            make_dir(joinpath(outpath, dir))
        end

        for file in files
            # TODO: non-regex solution ?
            m = match(Regex("^($tutorial_path)/(.+)"), root)

            if file[(end - 2):end] == ".jl"
                filename = file[begin:(end - 3)]
                Literate.markdown(joinpath(root, file), joinpath(outpath, String(m[2]));
                                    name = filename, execute=true, credit = false, documenter = true)
            else
                cp(joinpath(root, file), joinpath(outpath, String(m[2]), file);
                 force=true)
            end
        end
    end
end

@info "Building Tutorials"
generate_docs()

makedocs(
    sitename = "CUDATutorials.jl",
    authors = "Aditya Puranik",
    doctest = false,
    format = Documenter.HTML(prettyurls = false),
    pages = Any[
        "Home" => "index.md",
        "Introduction" => Any[
            "Introduction" => "introduction/01-Introduction.md",
            "Mandelbrot Set" => "./introduction/02-Mandelbrot_Set.md",
            "Shared Memory" => "./introduction/03-Shared_Memory.md",
            "Reduction" => "introduction/04-Reduction.md"
        ]
    ]
)


text:
SafeRational(num::T, den::T) where {T<:SafeSigned} = num//den
SafeRational(num::T, den::T) where {T<:Signed} = safeint(num)//safeint(den)

SafeRational(num::T1, den::T2) where {T1<:SafeSigned, T2<:SafeSigned} = SafeRational(promote(num, den)...)
SafeRational(num::T1, den::T2) where {T1<:Signed, T2<:Signed} = SafeRational(promote(safeint(num), safeint(den))...)

SafeRational(num::T1, den::T2) where {T1<:Signed, T2<:SafeSigned} = SafeRational(promote(safeint(num), den)...)
SafeRational(num::T1, den::T2) where {T1<:SafeSigned, T2<:Signed} = SafeRational(promote(num, safeint(den))...)

SafeRational(num::T, den::T) where {T<:SafeUnsigned} = num//den
SafeRational(num::T, den::T) where {T<:Unsigned} = safeint(num)//safeint(den)

SafeRational(num::T1, den::T2) where {T1<:SafeUnsigned, T2<:SafeUnsigned} = SafeRational(promote(num, den)...)
SafeRational(num::T1, den::T2) where {T1<:Unsigned, T2<:Unsigned} = SafeRational(promote(num, den)...)

SafeRational(num::T1, den::T2) where {T1<:Unsigned, T2<:SafeUnsigned} = SafeRational(promote(safeint(num), den)...)
SafeRational(num::T1, den::T2) where {T1<:SafeUnsigned, T2<:Unsigned} = SafeRational(promote(num, safeint(den))...)

SafeRational(num::T1, den::T2) where {T1<:SafeSigned, T2<:SafeUnsigned} = SafeRational(promote(num, den)...)
SafeRational(num::T1, den::T2) where {T1<:SafeUnsigned, T2<:SafeSigned} = SafeRational(promote(num, den)...)
SafeRational(num::T1, den::T2) where {T1<:Signed, T2<:Unsigned} = SafeRational(promote(num, den)...)
SafeRational(num::T1, den::T2) where {T1<:Unsigned, T2<:Signed} = SafeRational(promote(num, den)...)
SafeRational(num::T1, den::T2) where {T1<:Unsigned, T2<:SafeSigned} = SafeRational(promote(safeint(num), den)...)
SafeRational(num::T1, den::T2) where {T1<:SafeSigned, T2<:Unsigned} = SafeRational(promote(num, safeint(den))...)
SafeRational(num::T1, den::T2) where {T1<:Signed, T2<:SafeUnsigned} = SafeRational(promote(safeint(num), den)...)
SafeRational(num::T1, den::T2) where {T1<:SafeUnsigned, T2<:Signed} = SafeRational(promote(num, safeint(den))...)


text:
const COLORS = ["red", "white", "blue"]

function dutchsort!(a::Array{ASCIIString,1}, lo=COLORS[1], hi=COLORS[end])
    i = 1
    j = 1
    n = length(a)
    while j <= n
        if a[j] == lo
            a[i], a[j] = a[j], a[i]
            i += 1
            j += 1
        elseif a[j] == hi
            a[j], a[n] = a[n], a[j]
            n -= 1
        else
            j += 1
        end
    end
    return a
end

function dutchsort(a::Array{ASCIIString,1}, lo=COLORS[1], hi=COLORS[end])
    dutchsort!(copy(a), lo, hi)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    one_side_group(G::SimpleGraph)

Let `A` be the adjacency matrix of `G`.
Return the `Set` of `Permutation`s `p` such that `A[p,:]` is an adjacency matrix.
"""
function one_side_group(G::SimpleGraph)::Set{Permutation}
    result = Set{Permutation}()
    n = NV(G)
    A = adjacency(G)
    allow = make_allow(A)
    perms = RPG(n, allow)
    @info "There are $(length(perms)) permutations to consider"
    PM = Progress(length(perms))

    for p in perms
        next!(PM)
        B = A[p, :]
        if is_adjacency(B)
            push!(result, Permutation(p))
        end
    end

    return result
end

export one_side_group


text:
tests = [
    "cauchy",
    "geometric",
    "logistic",
    "lognormal",
    "normal",
    "uniform",
    "weibull"]

println("Running rmath tests:")

for t in tests
    test_fn = joinpath("rmath", "$t.jl")
    println(" * $test_fn")
    include(test_fn)
end


text:
export verifyorder

struct VerifyOrder{I} <: AbstractGorIter{I}
    rows::I
end

"""
    verifyorder(rows)
    rows |> verifyorder

Verify order of genome ordered stream `rows`.

An iterator that checks if its input is ordered by `(elt[1], elt[2])`.
The iterator passes the input rows through, but throws an
`ErrorException` if rows are out of order.
"""
function verifyorder(rows)
    VerifyOrder{typeof(rows)}(rows)
end

function Base.iterate(v::VerifyOrder)
    elt_s = iterate(v.rows)
    elt_s == nothing && return nothing
    elt, s = elt_s
    return elt, elt_s
end

function Base.iterate(v::VerifyOrder, last_s)
    last, s = last_s
    elt_s = iterate(v.rows, s)
    elt_s == nothing && return nothing
    elt, s = elt_s

    if (last.Chrom, last.Pos) > (elt.Chrom, elt.Pos)
        msg = "Out of order: ($(elt.Chrom), $(elt.Pos)) after ($(last.Chrom), $(last.Pos))"
        error(msg)
    end

    return elt, elt_s
end


text:
function search_recommender(request; environment=nothing)
    environment == nothing && @error "No search environment provided for search recommender."
    # Generate new  query and overwrite the original one
    request.query, id = generate_query(request.query, environment.dbdata, recommend_id_key=request.request_id_key)

    # Get the linear id of the entry for which recommendations are sought
    target_entry = db_select_entry(environment.dbdata, id, id_key=request.request_id_key)
    linear_id = isempty(target_entry) ? nothing : getproperty(target_entry, environment.id_key)

    # Search (end exclude original entry)
    search(environment, request; exclude=linear_id)
end

function generate_query(req::AbstractString, dbdata; recommend_id_key=DEFAULT_DB_ID_KEY)
    # Methods for transforming various julia types into search query string values
    __transform_value_for_search(value) = value  # default for non-floats, strings, leave untouched
    __transform_value_for_search(value::AbstractFloat) = "[$(0.9*value), $(1.1*value)]"
    __transform_value_for_search(value::AbstractString) = "\"" * value * "\""  # simple string, wrap in quotes

    # Parse a query generation request
    function __parse_request(req::AbstractString, dbschema; recommend_id_key=DEFAULT_DB_ID_KEY)
        # Get target id and fields from request string
        toks = split(req)
        _target_id, fields = try
            (string(toks[1]), Symbol.(toks[2:end]))
        catch e
            @debug "Could not parse target id, fields: $e"
            ("", Symbol[])
        end
        _target_id_type = typeof(_target_id)

        # Get target id type
        for dse in dbschema
            dse.column === recommend_id_key && (_target_id_type = dse.coltype)
        end

        # Parse target id
        target_id = __parse(_target_id_type, _target_id)
        return target_id, fields
    end

    # Initializations
    dbschema = db_create_schema(dbdata)
    target_id, fields = __parse_request(req, dbschema; recommend_id_key=recommend_id_key)
    columns = getproperty.(dbschema, :column)

    # Filter on id, find record
    target_record = db_select_entry(dbdata, target_id; id_key=recommend_id_key)

    # Extract all values using fields (if empty, use all default query generation fields)
    query_toks = Vector{String}()
    if !isempty(target_record)
        for field in fields
            if field in columns
                value = getproperty(target_record, field)
                value_for_query = __transform_value_for_search(value)
                push!(query_toks, "$field:$(value_for_query)")
            end
        end
    end
    query = join(query_toks, " ")
    return (query=query, id=target_id)
end


text:
@testset "SkewT" begin
  @testset "scalefromaltskewt, skewfromaltskewt" begin
    skew = -3
    scale = 2
    altskew = toaltskew(scale, skew)
    altscale = toaltscale(scale, skew)

    @test isapprox(scalefromaltskewt(altscale, altskew), scale)
    @test isapprox(skewfromaltskewt(altscale, altskew), skew)
  end

  @testset "Inference" begin
    Random.seed!(0)
    st = SkewT(2, 1, 7, -10)
    y = rand(st, 1000)
    model = MixSkewT(y, 1)
    print_model_info(model)
    nsamps, nburn, thin = (2000, 2000, 2)
    callback = make_callback(model, nsamps=nsamps, nburn=nburn, thin=thin)
    spl = make_sampler(model)
    @time chain, metrics = mcmc(spl, nsamps, nburn=nburn, thin=thin, callback=callback)

    # Tests
    @test length(unique(chain)) == nsamps
    @test isapprox(st.loc, mean(getindex.(chain, :mu))[1], atol=0.1)
    @test isapprox(st.scale, mean(getindex.(chain, :sigma))[1], atol=0.1)
    @test isapprox(st.df, mean(getindex.(chain, :nu))[1], atol=2)
    @test isapprox(st.skew, mean(getindex.(chain, :phi))[1], atol=3)
  end
end


text:
function goToWest(r)
    counter = 0

    while !isborder(r, West)
        counter+=1
        move!(r, West)
    end
    putmarker!(r)

    return counter
end

function chooseDirect(r)
    isborder(r, Sud) ? move!(r, Ost) : move!(r, Sud)
end

function putMarkers(r)
    directions = (Sud, Ost, Nord, West, Sud)

    for direct = directions
        while !isborder(r, direct) && !ismarker(r)
            putmarker!(r)
            move!(r, direct)
        end
    end
end

function goToStart(r, counter)
    for _ = 1:counter
        move!(r, Ost)
    end
end

function drawPerim(r)
    counter = goToWest(r)

    chooseDirect(r)

    putMarkers(r)

    goToStart(r, counter)
end

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
"""
    sparsify(M::AbstractArray, ratio::Real, p::Real, max_num_bins::Integer,
        impose_null_spaces=false::Bool)

Compute a `SparseMatrixCSC{Float64, Int64}` `X` for a matrix `M`, that is sparse and spectrally close to `M`,
especially in the lower end of the singular value spectrum.

Wrapper function over all functionality.

Compute a sparse approximation for a matrix `M`. First the sparsity pattern is computed with a
`ratio` in ``[0,1]`` that determines the sparsity (`0` means more sparse, `1` means less)
and with `p` in ``(0, \\infty ]`` that determines which norm is used to find the sparsity pattern.

Next the sparsity pattern is used to find a binning pattern, the non-negative integer `max_num_bins`
determines whats the maximum number of bins is (200-1000 is a reasonable choice). `max_num_bins=0`
means no binning is performed and can lead to significant slowdown.

Then an optimization problem, with constraints given by the binning pattern, is solved to find a
sparse approximation `X` for `M`.

If `impose_null_spaces=true`, another optimization problem is solved that ensures that `X` and
`M` have the same null space.

See also: [`sparsity_pattern`](@ref), [`bin_sparse_matrix!`](@ref).

# Examples
```julia-repl
julia> sparsify([16.99 65; 0.1 17.01], 0.6, 2, 200)
2×2 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
 16.8041  64.2499
   ⋅      16.8041
```
"""
function sparsify(M::AbstractArray{T}, ratio::Real, p::Real, max_num_bins::Integer,
    impose_null_spaces=false::Bool) where{T}

    pinv_M, rnull, lnull = pinv_qr(M)

    # sparsity pattern
    num_near_zero_rows, num_near_zero_cols = near_zero_row_col(M)
    min_per_row = max(0, min(size(rnull)[2] - num_near_zero_cols, size(M, 2)))
    min_per_col = max(0, min(size(lnull)[2] - num_near_zero_rows, size(M, 1)))

    M_id = sparsity_pattern(M, ratio, p, min_per_row, min_per_col)

    # binning pattern
    bin_sparse_matrix!(M, M_id, max_num_bins)

    # minimization
    pinv_MTM = pinv_M * pinv_M'
    pinv_MMT = pinv_M' * pinv_M
    D = M * pinv_MTM + pinv_MMT * M
    X = binned_minimization(M_id, pinv_MTM, pinv_MMT, D)
    impose_null_spaces && sps_impose_action!(X, M, rnull, lnull)
    return X
end

function sparsify(M::AbstractArray{<:Int}, ratio::Real, p::Real, max_num_bins::Integer, impose_null_spaces=false::Bool)
    return sparsify(float(M), ratio, p, max_num_bins, impose_null_spaces)
end

function near_zero_row_col(M::AbstractArray{T}) where {T}
    max_in_col = maximum(M, dims=1)
    max_in_row = maximum(M, dims=2)
    A_max = maximum(max_in_col)

    tol = A_max * eps(T) * 100  # magic constant from matlab

    num_near_zero_cols = sum(max_in_col .<= tol)
    num_near_zero_rows = sum(max_in_row .<= tol)
    return num_near_zero_rows, num_near_zero_cols
end

function near_zero_row_col(M::AbstractArray{<:Int})
    near_zero_row_col(float(M))
end

"""
    binned_minimization(M_id::AbstractSparseMatrixCSC,
    B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix)

Solves a optimization problem of the from

```math
Y A + B Y = D.
```
under the binning constraints given in the `SparseMatrixCSC{Int64, Int64}` `M_id`.

See also: [`sparsify`](@ref).

# Examples
```julia-repl
julia> binned_minimization([16.99 65 64.96; 0.1 17.01 0.1], sparse([1 2 2; 0 1 0]), Matrix(I, 3, 3) , Matrix(I, 2, 2), Matrix(I, 2, 3))
2×3 SparseMatrixCSC{Float64, Int64} with 4 stored entries:
 0.5  0.0  0.0
  ⋅   0.5   ⋅
```
"""
function binned_minimization(M_id::AbstractSparseMatrixCSC,
    B::AbstractMatrix{T}, C::AbstractMatrix{T}, D::AbstractMatrix{T}) where{T}
    m = size(M_id, 1)
    n = size(M_id, 2)

    0 <= minimum(M_id) || throw(DomainError(minimum(M_id), "M_id contains negative integers"))
    size(B)[1] == size(B)[2] == n || throw(ArgumentError("B has wrong size"))
    size(C)[1] == size(C)[2] == m || throw(ArgumentError("C has wrong size"))
    size(D, 1) == size(M_id, 1) && size(D, 2) == size(M_id, 2) || throw(ArgumentError("D must be the same size as M"))

    LS_M, LS_b = system_no_null(M_id, B, C, D)
    LS_M, = pinv_qr(LS_M)

    LS_x = LS_M * LS_b

    Y = unknown_to_matrix(LS_x, M_id)
    return Y
end

function system_no_null(M_id::AbstractSparseMatrixCSC, B::AbstractMatrix{T}, C::AbstractMatrix{T}, D::AbstractMatrix{T}) where{T}
    m = size(M_id, 1)
    n = size(M_id, 2)

    0 <= minimum(M_id) || throw(DomainError(minimum(M_id), "M_id contains negative integers"))
    size(B)[1] == size(B)[2] == n || throw(ArgumentError("B has wrong size"))
    size(C)[1] == size(C)[2] == m || throw(ArgumentError("C has wrong size"))
    size(D, 1) == size(M_id, 1) && size(D, 2) == size(M_id, 2) || throw(ArgumentError("D must be the same size as M"))

    N = length(unique(M_id.nzval))

    LS_A = zeros(T, (N, N))
    LS_b = zeros(T, N)

    for idx in findall(!iszero, M_id)
        J, L = findnz(M_id[Tuple(idx)[1], :])
        k = 1
        for l in L
            @views LS_A[M_id[idx], l] = LS_A[M_id[idx], l] + B[J[k], Tuple(idx)[2]]
            k += 1
        end

        J, L = findnz(M_id[:, Tuple(idx)[2]])
        k = 1
        for l in L
            @views LS_A[M_id[idx], l] = LS_A[M_id[idx], l] + C[Tuple(idx)[1], J[k]]
            k += 1
        end
    end

    @views for i in 1:N
        LS_b[i] = sum(D[M_id .== i])
    end

    return LS_A, LS_b
end

function unknown_to_matrix(x::AbstractVector, M_id::AbstractSparseMatrixCSC)
    m = size(M_id, 1)
    n = size(M_id, 2)

    M_id_i, M_id_j, M_id_v = findnz(M_id)
    return sparse(M_id_i, M_id_j, x[M_id_v], m, n)
end

function sps_impose_action!(Y::AbstractSparseMatrixCSC, M::AbstractArray{T}, R_mat::AbstractMatrix, L_mat::AbstractMatrix,
     tol=eps(T)::Real, max_iters=300::Integer) where{T}
    tol > 0 || throw(DomainError("tol has to be greater than zero"))
    size(M, 1) == size(Y, 1) && size(M, 2) == size(Y, 2) || throw(ArgumentError("Y must be the same size as"))

    m = size(Y, 1)
    n = size(Y, 2)
    nR1, nR2 = size(R_mat)
    nL1, nL2 = size(L_mat)

    n == nR1 || throw(ArgumentError("Y * R_mat not possible"))
    m == nL1 || throw(ArgumentError("Y' * L_mat not possible"))

    # Quantities don't change when iterating
    norm_R_mat = norm(R_mat)
    norm_L_mat = norm(L_mat)

    M_R_mat = M * R_mat
    M_L_mat = M' * L_mat

    nnz_Y = count(!iszero, Y)
    Y_i, Y_j, _ = findnz(Y)

    # Quantities that do change when iterating
    Lag_R = zeros(T, m, nR2)
    Lag_L = zeros(T, n, nL2)

    d_R = Y * R_mat - M_R_mat
    d_L = Y' * L_mat - M_L_mat

    q_R = - d_R
    q_L = - d_L

    d_R_proj_vec = zeros(T, nnz_Y, 1)
    d_L_proj_vec = zeros(T, nnz_Y, 1)
    d_proj_vec = zeros(T, nnz_Y, 1)

    iters = 1
    for i in 1:max_iters
        if norm(d_R) <= tol * n * norm(Y - M) * norm_R_mat  norm(d_L) <= tol * m * norm(Y - M) * norm_L_mat
            break
        end

        @views for k in 1:nnz_Y
            d_R_proj_vec[k] = R_mat[Y_j[k], :]' * d_R[Y_i[k], :]
            d_L_proj_vec[k] = L_mat[Y_i[k], :]' * d_L[Y_j[k], :]
        end

        d_proj_vec = d_R_proj_vec .+ d_L_proj_vec

        alpha = (norm(q_R)^2 + norm(q_L)^2) / norm(d_proj_vec)^2

        Lag_R = Lag_R .+ alpha .* d_R
        Lag_L = Lag_L .+ alpha .* d_L

        @views for k in 1:nnz_Y
            Y[Y_i[k], Y_j[k]] = Y[Y_i[k], Y_j[k]] - alpha * d_proj_vec[k]
        end

        norm2_q_RL_old = norm(q_R)^2 + norm(q_L)^2

        q_R = M_R_mat .- Y * R_mat
        q_L = M_L_mat .- Y' * L_mat

        norm2_q_RL_new = norm(q_R)^2 + norm(q_L)^2

        beta = norm2_q_RL_new / norm2_q_RL_old
        d_R = beta .* d_R .- q_R
        d_L = beta .* d_L .- q_L

        iters = iters + 1
    end
    return Y
end

function sps_impose_action!(Y::AbstractSparseMatrixCSC, M::AbstractArray{<:Int}, R_mat::AbstractMatrix,
    L_mat::AbstractMatrix, tol::Real, max_iters=1000::Integer)
    return sps_impose_action!(Y, float(M), R_mat, L_mat, eps(eltype(M)), max_iters)
end


text:
using Test, Trapla

include("gevrey_test.jl")


text:
module TestParameterLabels

using Mimi
using NamedArrays
using Test

import Mimi:
    reset_compdefs

reset_compdefs()

############################################
#    BASIC TEST - use NamedArrays (1/3)    #
############################################

@defcomp compA begin
    regions = Index()

    y = Variable(index=[time, regions])
    x = Parameter(index=[time, regions])

    function run_timestep(p, v, d, t)
        for r in d.regions
            v.y[t, r] = p.x[t, r]
        end
    end
end

x = Array{Float64}(undef, 20, 3)
for t in 1:20
    x[t,1] = 1
    x[t,2] = 2
    x[t,3] = 3
end

region_labels = ["Region1", "Region2", "Region3"]
time_labels = collect(2015:5:2110)

x2 = NamedArray(Array{Float64}(undef, 20, 3), (time_labels, region_labels), (:time, :regions))
for t in time_labels
    x2[:time => t, :regions => "Region1"] = 1
    x2[:time => t, :regions => "Region2"] = 2
    x2[:time => t, :regions => "Region3"] = 3
end

model1 = Model()
set_dimension!(model1, :time, time_labels)
set_dimension!(model1, :regions, region_labels)
add_comp!(model1, compA)
set_param!(model1, :compA, :x, x)

model2 = Model()
set_dimension!(model2, :time, time_labels)
set_dimension!(model2, :regions, region_labels)
add_comp!(model2, compA)
set_param!(model2, :compA, :x, x2) # should perform parameter dimension check

run(model1)
run(model2)

for t in 1:length(time_labels)
    for r in 1:length(region_labels)
        @test(model1[:compA, :y][t, r] == model2[:compA, :y][t, r])
    end
end

@test(size(getdataframe(model1, :compA, :y)) == (60, 3))

#####################################
#  LARGER MULTIREGIONAL TEST (2/3)  #
#####################################

#GROSS ECONOMY COMPONENT
@defcomp grosseconomy begin
    regions = Index()                           #Note that a regional index is defined here

    YGROSS  = Variable(index=[time, regions])   #Gross output
    K       = Variable(index=[time, regions])   #Capital
    l       = Parameter(index=[time, regions])  #Labor
    tfp     = Parameter(index=[time, regions])  #Total factor productivity
    s       = Parameter(index=[time, regions])  #Savings rate
    depk    = Parameter(index=[regions])        #Depreciation rate on capital - Note that it only has a region index
    k0      = Parameter(index=[regions])        #Initial level of capital
    share   = Parameter()                       #Capital share

    function run_timestep(p, v, d, t)
        #Define an equation for K
        for r in d.regions
            if is_first(t)
                v.K[t,r] = p.k0[r]
            else
                v.K[t,r] = (1 - p.depk[r])^5 * v.K[t-1,r] + v.YGROSS[t-1,r] * p.s[t-1,r] * 5
            end
        end

        #Define an equation for YGROSS
        for r in d.regions
            v.YGROSS[t,r] = p.tfp[t,r] * v.K[t,r]^p.share * p.l[t,r]^(1-p.share)
        end
    end
end

#EMISSIONS COMPONENT
@defcomp emissions begin
    regions     = Index()                           #The regions index must be specified for each component

    E           = Variable(index=[time, regions])   #Total greenhouse gas emissions
    E_Global    = Variable(index=[time])            #Global emissions (sum of regional emissions)
    sigma       = Parameter(index=[time, regions])  #Emissions output ratio
    YGROSS      = Parameter(index=[time, regions])  #Gross output - Note that YGROSS is now a parameter

    function run_timestep(p, v, d, t)
        #Define an eqation for E
        for r in d.regions
            v.E[t,r] = p.YGROSS[t,r] * p.sigma[t,r]
        end

        #Define an equation for E_Global
        for r in d.regions
            v.E_Global[t] = sum(v.E[t,:])
        end
    end
end

#DEFINE ALL THE PARAMETERS
l = Array{Float64}(undef, 20, 3)
for t in 1:20
    l[t,1] = (1. + 0.015)^t *2000
    l[t,2] = (1. + 0.02)^t * 1250
    l[t,3] = (1. + 0.03)^t * 1700
end

tfp = Array{Float64}(undef, 20, 3)
for t in 1:20
    tfp[t,1] = (1 + 0.06)^t * 3.2
    tfp[t,2] = (1 + 0.03)^t * 1.8
    tfp[t,3] = (1 + 0.05)^t * 2.5
end

s = Array{Float64}(undef, 20, 3)
for t in 1:20
    s[t,1] = 0.21
    s[t,2] = 0.15
    s[t,3] = 0.28
end

depk = [0.11, 0.135 ,0.15]
k0   = [50.5, 22., 33.5]

sigma = Array{Float64}(undef, 20, 3)
for t in 1:20
    sigma[t,1] = (1. - 0.05)^t * 0.58
    sigma[t,2] = (1. - 0.04)^t * 0.5
    sigma[t,3] = (1. - 0.045)^t * 0.6
end

#FUNCTION TO RUN MY MODEL
function run_my_model()

    my_model = Model()

    set_dimension!(my_model, :time, collect(2015:5:2110))
    set_dimension!(my_model, :regions, ["Region1", "Region2", "Region3"])  #Note that the regions of your model must be specified here

    add_comp!(my_model, grosseconomy)
    add_comp!(my_model, emissions)

    set_param!(my_model, :grosseconomy, :l, l)
    set_param!(my_model, :grosseconomy, :tfp, tfp)
    set_param!(my_model, :grosseconomy, :s, s)
    set_param!(my_model, :grosseconomy, :depk, depk)
    set_param!(my_model, :grosseconomy, :k0, k0)
    set_param!(my_model, :grosseconomy, :share, 0.3)

    #set parameters for emissions component
    set_param!(my_model, :emissions, :sigma, sigma2)
    connect_param!(my_model, :emissions, :YGROSS, :grosseconomy, :YGROSS)

    run(my_model)
    return(my_model)

end

#DEFINE ALL THE PARAMETERS using NAMEDARRAYS
region_labels = ["Region1", "Region2", "Region3"]
time_labels = collect(2015:5:2110)

l2 = NamedArray(Array{Float64}(undef, 20, 3), (time_labels, region_labels), (:time, :regions))
for t in 1:length(time_labels)
    l2[:time => time_labels[t], :regions => region_labels[1]] = (1. + 0.015)^t *2000
    l2[:time => time_labels[t], :regions => region_labels[2]] = (1. + 0.02)^t * 1250
    l2[:time => time_labels[t], :regions => region_labels[3]] = (1. + 0.03)^t * 1700
end

tfp2 = NamedArray(Array{Float64}(undef, 20, 3), (time_labels, region_labels), (:time, :regions))
for t in 1:length(time_labels)
    tfp2[:time => time_labels[t], :regions => region_labels[1]] = (1 + 0.06)^t * 3.2
    tfp2[:time => time_labels[t], :regions => region_labels[2]] = (1 + 0.03)^t * 1.8
    tfp2[:time => time_labels[t], :regions => region_labels[3]] = (1 + 0.05)^t * 2.5
end

s2 = NamedArray(Array{Float64}(undef, 20, 3), (time_labels, region_labels), (:time, :regions))
for t in 1:length(time_labels)
    s2[t, 1] = 0.21
    s2[t, 2] = 0.15
    s2[t, 3] = 0.28
end

depk2 = NamedArray([0.11, 0.135 ,0.15], (region_labels,), (:regions,))
k02   = NamedArray([50.5, 22., 33.5], (region_labels,), (:regions,))

sigma2 = NamedArray(Array{Float64}(undef, 20, 3), (time_labels, region_labels), (:time, :regions))
for t in 1:length(time_labels)
    sigma2[t, 1] = (1. - 0.05)^t * 0.58
    sigma2[t, 2] = (1. - 0.04)^t * 0.5
    sigma2[t, 3] = (1. - 0.045)^t * 0.6
end

function run_my_model2()

    my_model2 = Model()

    set_dimension!(my_model2, :time, collect(2015:5:2110))
    set_dimension!(my_model2, :regions, ["Region1", "Region2", "Region3"])  #Note that the regions of your model must be specified here

    add_comp!(my_model2, grosseconomy)
    add_comp!(my_model2, emissions)

    set_param!(my_model2, :grosseconomy, :l, l2)
    set_param!(my_model2, :grosseconomy, :tfp, tfp2)
    set_param!(my_model2, :grosseconomy, :s, s2)
    set_param!(my_model2, :grosseconomy, :depk,depk2)
    set_param!(my_model2, :grosseconomy, :k0, k02)
    set_param!(my_model2, :grosseconomy, :share, 0.3)

    #set parameters for emissions component
    set_param!(my_model2, :emissions, :sigma, sigma2)
    connect_param!(my_model2, :emissions, :YGROSS, :grosseconomy, :YGROSS)

    run(my_model2)
    return(my_model2)

end

run1 = run_my_model()
run2 = run_my_model2()

#Check results

for t in 1:length(time_labels)
    for r in 1:length(region_labels)
        @test(run1[:grosseconomy, :YGROSS][t, r] == run2[:grosseconomy, :YGROSS][t, r])
        #println(run1[:grosseconomy, :YGROSS][t, r],", ", run2[:grosseconomy, :YGROSS][t, r])
        @test(run1[:grosseconomy, :K][t, r] == run2[:grosseconomy, :K][t, r])
        @test(run1[:emissions, :E][t, r] == run2[:emissions, :E][t, r])
    end
    @test(run1[:emissions, :E_Global][t] == run2[:emissions, :E_Global][t])
end

######################################################
#  set_param! option with list of dimension names  #
######################################################

model3 = Model()
set_dimension!(model3, :time, collect(2015:5:2110))
set_dimension!(model3, :regions, ["Region1", "Region2", "Region3"])
add_comp!(model3, compA)
set_param!(model3, :compA, :x, x, [:time, :regions])
run(model3)

for t in 1:length(time_labels)
    for r in 1:length(region_labels)
        @test(model1[:compA, :y][t, r] == model3[:compA, :y][t, r])
    end
end

end #module


text:
using StatsFuns
using Test
using InverseFunctions

@testset "inverse" begin
    x = 0.7
    InverseFunctions.test_inverse(norminvcdf, x)
    InverseFunctions.test_inverse(norminvccdf, x)
    InverseFunctions.test_inverse(norminvlogcdf, log(x))
    InverseFunctions.test_inverse(norminvlogccdf, log(x))
end


text:
if isdefined(Base, :LIBEXECDIR)
  const exe7z = joinpath(Sys.BINDIR, Base.LIBEXECDIR, "7z.exe")
else
  const exe7z = joinpath(Sys.BINDIR, "7z.exe")
end
using FFTW
println("num threads set to 1")
FFTW.set_num_threads(1)


text:
const TEST_NAME = "PCG_RXS_M_XS_64_64"

include("common.jl")

using RNG.PCG

r = PCGStateSetseq(UInt64, PCG_RXS_M_XS, (0x018cd83e277674ac, 0x436cd6f2434be066))

test_all(r, 100_000_000)


text:
struct NotImplementedError <: Exception
    msg::String
    NotImplementedError(s) = new(
        """
        [$s] hasn't been implemented.
        """
    )
end

struct UnknownOptionError <: Exception
    msg::String
    UnknownOptionError(s, o) = new(
        """
        [$s] is not recognised as a valid option name for $(typeof(o)).
        """
    )
end

struct OptionValueError <: Exception
    msg::String
    OptionValueError(s, v, r="") = new(
        """
        [$s] value given ($v) didn't meet the expected format$(
            ifelse(isempty(r), ".", " $r."))
        """
    )
end


text:
export mh

immutable Choice
    store::Store
    k::Function
    address
    erp::ERP
    score::Float64
    value
    choicescore::Float64
    fresh::Bool # true when re-sampled.
end

type MH <: Ctx
    choices::Vector{Choice}
    oldchoices::Vector{Choice}
    score::Float64
end

MH() = MH(Choice[], Choice[], 0)

function sample(s::Store, k::Function, address, e::ERP, ctx::MH, force=false)
    # This makes a step run in O(n^2) time. (Though the constants are
    # really huge anyway.) For larger state spaces it might be worth
    # using an ordered dictionary.
    if force || ((ix = findfirst(c -> c.address == address, ctx.oldchoices)) == 0)
        value = rand(e)
        fresh = true
    else
        value = ctx.oldchoices[ix].value
        fresh = false
    end
    choicescore = score(e, value)
    push!(ctx.choices, Choice(s,k,address,e,ctx.score,value,choicescore,fresh))
    # Update score /after/ storing current score in choices in-case of restart.
    ctx.score += choicescore
    k(s,value)
end

function factor(s::Store, k::Function, a, score, ctx::MH)
    ctx.score += score
    k(s,nothing)
end

mhexit(s::Store, value) = value

function mh(store::Store, k::Function, address, comp::Function, numsteps=10)
    global ctx
    counts = Dict{Any,Int64}()
    ctxold, ctx = ctx, MH()
    try
        # TODO: Consider introducing a compound type to hold choices,
        # score, value. It could support indexing so that state[1:j]
        # does the right thing.

        # Initialize.
        oldvalue = trampoline(() -> comp(store, mhexit, address))
        oldscore = ctx.score
        ctx.oldchoices = ctx.choices

        for i in 1:numsteps
            j = rand(1:length(ctx.oldchoices))
            restartchoice = ctx.oldchoices[j]

            # Restore state prior to re-start point.
            ctx.choices = ctx.oldchoices[1:j-1]
            ctx.score = restartchoice.score

            # Re-start.
            value = trampoline() do
                sample(restartchoice.store,
                       restartchoice.k,
                       restartchoice.address,
                       restartchoice.erp,
                       ctx, true)
            end

            if rand() < min(1., exp(log_mh_ratio(ctx.choices, ctx.oldchoices,
                                                 ctx.score, oldscore, j)))
                # Accept proposal.
                ctx.oldchoices = ctx.choices
                oldscore = ctx.score
                oldvalue = value
            end

            # Update counts.
            counts[oldvalue] = get(counts, oldvalue, 0) + 1
        end
    finally
        ctx = ctxold
    end
    k(store, Discrete(counts))
end

function log_mh_ratio(choices, oldchoices, score, oldscore, j)
    fw = -log(length(oldchoices))
    for c in choices[j:end]
        if c.fresh
            fw += c.choicescore
        end
    end
    bw = -log(length(choices))
    # This is also O(n^2) per-step.
    for c in oldchoices[j:end]
        ix = findfirst(d -> d.address == c.address, choices)
        if ix==0 || choices[ix].fresh
            bw += c.choicescore
        end
    end
    score - oldscore + bw - fw
end


text:
reward_states = [[4, 1], [5, 5], [7, 6]]

# MDPs initialization
mdp = PyramidMDP(reward_states=reward_states)

# check implementation of required methods
# @POMDPLinter.show_requirements FiniteHorizonPOMDPs.solve(fhsolver, mdp)

# initialize the solver
# max_iterations: maximum number of iterations value iteration runs for (default is 100)
# belres: the value of Bellman residual used in the solver
solver = ValueIterationSolver(max_iterations=10, belres=1e-3, include_Q=true);
fhsolver = FiniteHorizonSolver()

# Solve Value Iteration
VIPolicy = DiscreteValueIteration.solve(solver, mdp);

# Solve Finite Horizon by Value Iteration
FHPolicy = FiniteHorizonPOMDPs.solve(fhsolver, mdp);

@test typeof(HorizonLength(mdp)) != InfiniteHorizon

# Compare resulting policies
@test all(FiniteHorizonValueIteration.action(FHPolicy, s) == DiscreteValueIteration.action(VIPolicy, s) for s in states(mdp))

# Compare FHMDP and IHMDP states
fh_states = Iterators.flatten([FiniteHorizonPOMDPs.stage_states(mdp, i) for i=1:FiniteHorizonPOMDPs.horizon(mdp) + 1])
ih_states = states(mdp)

z = zip(fh_states, ih_states)

@test all((fh == ih for (fh, ih) in z))


text:
module digit_mods

export interpolate, str_el, ΔeltaAB

function ΔeltaAB(A,B;ABS=false)
    size(A)[1] != size(B)[1] ? error("A has not the same rows as B") :

    return ABS ? [abs((A[i,2]-B[i,2])/B[i,2]) for i in 1:size(A)[1]] : [(A[i,2]-B[i,2])/B[i,2] for i in 1:size(A)[1]]
end

function interpolate(a::Array{<:Real}, x::Array{<:Real};side="y")
    # a is input array with a x and an y axis. VALUES ARE IN ORDER AND NOT RANDOMLY PLACED
    # x is the array which carries the values where we want to interpolate data out of a
    (minimum(size(a)) >= 3) | (length(size(a)) > 2) ?  error("Check input array's dimensions") :
    (length(size(x)) > 2) ?  error("Check interpolation array's dimensions") :
    size(a)[1] < size(a)[2] ? a_len = 2 : a_len = 1 # check whether is a row or column matrix
    if side == "y"
       a_len == 1 ? xa = a[:,1] : xa = a[1,:]
       a_len == 1 ? ya = a[:,2] : ya = a[2,:]
    end

    res = Array{Real}(undef,size(x))

    for k in 1:size(x)[1]
        i = x[k]
        for j in 1:size(a)[a_len] - 1
           ((xa[j] <= i) & (xa[j + 1] >= i)) | ((xa[j] >= i) & (xa[j + 1] <= i))  ?  (res[k] = ya[j] + (ya[j + 1] - ya[j]) / (xa[j + 1] - xa[j]) * (i - xa[j]);break) : continue
        end
    end
    return res
end

function str_el(data;concat=false)
   # empty strings elimination
    for i ∈ 1:size(data)[1]
        pos = []
        for j ∈ 1:size(data[i])[1]
            if !((typeof(data[i][j]) == Float64) || (typeof(data[i][j]) == Int64))
                push!(pos, j)
            end
        end
        for p ∈ 1:size(pos)[1]
            deleteat!(data[i], pos[size(pos)[1] - p + 1])
        end
    end

    ((size(data)[1] == 1) & concat) ? (concat = false; @warn "Data were a single column,thus non concatenation is needed.") :
    if concat
        tmp = data[1]
        for i in 2:size(data)[1]
            length(data[1]) != length(data[i]) ? ( @warn "Column $i is not the same length as the previous. Skipping it.";continue) :
            tmp = hcat(tmp, data[i])
        end
        data = tmp
    else
        data = [i for i in data[1]]
    end
    return data
end
# new VS plot for different meshes against the experiment
export popul
function popul(path::AbstractString,a::AbstractArray)
    file = open(path*"Cp_RAE_up_var.dat")
    lines = readlines(file)
    for i in 1:size(lines)[1]
       a[i,1] = parse(Float32,lines[i])
    end
    close(file);
    file = open(path*"Cp_RAE_down_var.dat");
    lines = readlines(file);
    for i in 1:size(lines)[1]
       a[i,2] = parse(Float32,lines[i])
    end
    close(file)
end

end


text:
include(joinpath("linear", "integrator.jl"))
include(joinpath("linear", "body.jl"))
include(joinpath("linear", "test.jl"))


text:
@concrete struct SimpleOrbit <: AbstractOrbit
    period
    t0
    b
    duration
    speed
    half_period
    ref_time
end

"""
    SimpleOrbit(; period, duration, t0=0, b=0.0)

Circular orbit parameterized by the basic observables of a transiting system.

# Parameters
* `period` - The orbital period of the planets, nominally in days
* `duration` - The duration of the transit, similar units as `period`.
* `t0` - The midpoint time of the reference transit, similar units as `period`
* `b` - The impact parameter of the orbit, unitless
"""
function SimpleOrbit(;period, duration, t0=zero(period), b=0)
    half_period = 0.5 * period
    duration > half_period && error("duration cannot be longer than half the period")
    speed = 2.0 * sqrt(1.0 - b^2) / duration
    ref_time =  t0 - half_period
    SimpleOrbit(period, t0, b, duration, speed, half_period, ref_time)
end

period(orbit::SimpleOrbit) = orbit.period
duration(orbit::SimpleOrbit) = orbit.duration

relative_time(orbit::SimpleOrbit, t) = mod(t - orbit.ref_time, period(orbit)) - orbit.half_period

function relative_position(orbit::SimpleOrbit, t)
    Δt = relative_time(orbit, t)
    x = orbit.speed * Δt
    y = orbit.b
    z = abs(Δt) < 0.5 * duration(orbit) ? one(x) : -one(x)
    return SA[x, y, z]
end

# TODO: if texp, ϵ += 0.5 * texp
#function in_transit(orbit::SimpleOrbit, t)
#    Δt = relative_time.(orbit, t)
#    ϵ = 0.5 * duration(orbit)
#    return findall(x -> abs(x) < ϵ, Δt)
#end
#function in_transit(orbit::SimpleOrbit, t, r)
#    Δt = relative_time.(orbit, t)
#    ϵ = √((1.0 + r)^2.0 - orbit.b^2.0) / orbit.speed
#    return findall(x -> abs(x) < ϵ, Δt)
#end

function flip(orbit::SimpleOrbit, ror)
    t0 = orbit.t0 + orbit.half_period
    b = orbit.b / ror
    speed = orbit.speed / ror
    ref_time = orbit.t0
    return SimpleOrbit(period(orbit), t0, b, duration(orbit), speed, orbit.half_period, ref_time)
end

function Base.show(io::IO, orbit::SimpleOrbit)
    T = duration(orbit)
    P = period(orbit)
    b = orbit.b
    t0 = orbit.t0
    print(io, "SimpleOrbit(P=$P, T=$T, t0=$t0, b=$b)")
end

function Base.show(io::IO, ::MIME"text/plain", orbit::SimpleOrbit)
    T = duration(orbit)
    P = period(orbit)
    b = orbit.b
    t0 = orbit.t0
    print(io,
        "SimpleOrbit\n period: ", P,
        "\n duration: ", T,
        "\n t0: ", t0,
        "\n b: ", b
    )
end


text:
using SafeTestsets
include("common.jl")

@testset "Jacobian Sparsity" begin include("jacobian.jl") end
@testset "Paraboloid example" begin include("paraboloid.jl") end
@testset "PDE with globals" begin include("global_PDE.jl") end

@testset "Hessian sparsity" begin include("hessian.jl") end

@testset "Exploration" begin include("ifsandbuts.jl") end
@testset "Examples" begin include("examples.jl") end


text:
#!/usr/bin/env julia
# Tim Sterne-Weiler 2015

const dir = abspath( splitdir(@__FILE__)[1] )
const ver = readline(open(dir * "/VERSION"))

tic()
println( STDERR, "Whippet $ver loading and compiling... " )

using ArgParse

unshift!( LOAD_PATH, dir * "/../src" )
using Whippet

function parse_cmd()
  s = ArgParseSettings()
  # TODO finish options...
  @add_arg_table s begin
    "filename.fastq[.gz]"
      arg_type = String
      required = true
    "paired_mate.fastq[.gz]"
      arg_type = String
      required = false
    "--index", "-x"
      help     = "Output prefix for saving index 'dir/prefix' (default Whippet/index/graph[.jls])"
      arg_type = String
      default  = fixpath( "$dir/../index/graph" )
    "--out", "-o"
      help     = "Where should the gzipped output go 'dir/prefix'?"
      arg_type = String
      default  = fixpath( "./output" )
    "--seed-len", "-L"
      help     = "Seed length"
      arg_type = Int
      default  = 18
    "--seed-try", "-M"
      help     = "Number of failed seeds to try before giving up"
      arg_type = Int
      default  = 3
    "--seed-tol", "-T"
      help     = "Number of seed hits to tolerate"
      arg_type = Int
      default  = 4
    "--seed-buf", "-B"
      help     = "Ignore this many bases from beginning and end of read for seed"
      arg_type = Int
      default  = 5
    "--seed-inc", "-I"
      help     = "Number of bases to increment seed each iteration"
      arg_type = Int
      default  = 18
    "--pair-range", "-P"
      help     = "Seeds for paired end reads must match within _ bases of one another"
      arg_type = Int
      default  = 5000
    "--mismatches", "-X"
      help     = "Allowable number of mismatches in alignment (counted as 1-10^(-phred/10))"
      arg_type = Int
      default  = 3
    "--score-min", "-S"
      help     = "Minimum percent matching (matches - mismatches) / read_length"
      arg_type = Float64
      default  = 0.7
    "--sam"
      help     = "Should SAM format be sent to stdout?"
      action   = :store_true
    "--biascorrect"
      help     = "Apply fragment GC-content and 5' Sequence bias correction methods for more stable PSI values"
      action   = :store_true
    "--stranded"
      help     = "Is the data strand specific in fwd orientation? If so, increase speed with this flag"
      action   = :store_true
    "--pair-same-strand"
      help     = "Whippet by default tries to align fwd/rev pairs, if your data is fwd/fwd or rev/rev set this flag"
      action   = :store_true
    "--phred-33"
      help     = "Qual string is encoded in Phred+33 integers (default)"
      action   = :store_true
    "--phred-64"
      help     = "Qual string is encoded in Phred+64 integers"
      action   = :store_true
    "--url"
      help     = "When this flag is enabled, filename.fastq[.gz] and paired_mate.fastq[.gz] are considered URLs to download/process on the fly"
      action   = :store_true
    "--ebi"
      help     = "Retrieve FASTQ file(s) from ebi.ac.uk using a single seq run id (ie. ERR1994736). This will automatically work for both single and paired-end SRR ids. (sets --url=true)"
      action   = :store_true
    "--circ"
      help     = "Allow back/circular splicing, this will allow output of `BS`-type lines"
      action   = :store_true
    "--force-gz"
      help     = "Regardless of suffix, consider read input as gzipped"
      action   = :store_true
  end
  return parse_args(s)
end

function main()

   args = parse_cmd()

   println(STDERR, " $( round( toq(), 6 ) ) seconds" )

   indexpath = fixpath( args["index"] )
   indexname = hasextension(indexpath, ".jls") ? indexpath : indexpath * ".jls"
   println(STDERR, "Loading splice graph index... $indexname")
   @timer const lib = open(deserialize, indexname)

   ispaired = args["paired_mate.fastq[.gz]"] != nothing ? true : false

   if args["ebi"]
      ebi_res = ident_to_fastq_url( args["filename.fastq[.gz]"] )
      ispaired = ebi_res.paired
      args["url"] = true
      args["filename.fastq[.gz]"] = "http://" * ebi_res.fastq_1_url
      if ispaired
         args["paired_mate.fastq[.gz]"] = "http://" * ebi_res.fastq_2_url
      end
      ebi_res.success || error("Could not fetch data from ebi.ac.uk!!")
   elseif !args["url"]
      args["filename.fastq[.gz]"] = fixpath(args["filename.fastq[.gz]"])
      if ispaired
         args["paired_mate.fastq[.gz]"] = fixpath(args["paired_mate.fastq[.gz]"])
      end
   end

   const ContainerType = ispaired ? SGAlignPaired : SGAlignSingle

   if args["biascorrect"]
      const CounterType   = JointBiasCounter
      const BiasModelType = JointBiasMod
   else
      const CounterType   = DefaultCounter
      const BiasModelType = DefaultBiasMod
   end

   const param = AlignParam( args, ispaired, kmer=lib.kmer )
   const quant = GraphLibQuant{ContainerType,CounterType}( lib )
   const multi = MultiMapping{ContainerType,CounterType}()
   const mod   = BiasModelType()

   const enc_offset = args["phred-64"] ? 64 : 33

   const parser,response = args["url"] ? make_http_fqparser( args["filename.fastq[.gz]"],
                                         forcegzip=args["force-gz"] ) :
                                         make_fqparser( args["filename.fastq[.gz]"],
                                         forcegzip=args["force-gz"] )

   if ispaired
      const mate_parser,mate_response = args["url"] ? make_http_fqparser( args["paired_mate.fastq[.gz]"],
                                                      forcegzip=args["force-gz"] ) :
                                                      make_fqparser( args["paired_mate.fastq[.gz]"],
                                                      forcegzip=args["force-gz"] )
   end

   if nprocs() > 1
      println(STDERR, "Whippet does not currrently support nprocs() > 1")
      return #TODO: first implementation was too slow, ie too much communication overhead
   else
      if args["url"]
         println(STDERR, "Streaming reads from...")
      else
         println(STDERR, "Processing reads from file...")
      end
      if ispaired
         println(STDERR, "FASTQ_1: " * args["filename.fastq[.gz]"])
         println(STDERR, "FASTQ_2: " * args["paired_mate.fastq[.gz]"])
         @timer mapped,totreads,readlen = process_paired_reads!( parser, mate_parser, param, lib, quant, multi, mod,
                                                             sam=args["sam"], qualoffset=enc_offset,
                                                             response=response, mate_response=mate_response,
                                                             http=args["url"] )
         readlen = round(Int, readlen)
         println(STDERR, "Finished mapping $mapped paired-end reads of length $readlen each out of a total of $totreads mate-pairs...")
      else
         println(STDERR, "FASTQ: " * args["filename.fastq[.gz]"])
         @timer mapped,totreads,readlen = process_reads!( parser, param, lib, quant, multi, mod,
                                                       sam=args["sam"], qualoffset=enc_offset,
                                                       response=response, http=args["url"] )
         readlen = round(Int, readlen)
         println(STDERR, "Finished mapping $mapped single-end reads of length $readlen out of a total of $totreads reads...")
      end
   end

   # TPM_EM
   println(STDERR, "Calculating expression values and MLE of equivalence classes with EM:")
   primer_normalize!( mod )
   primer_adjust!( quant, mod )
   primer_adjust!( multi, mod )
   println(STDERR, "- $( length(multi.map) ) multi-gene mapping read equivalence classes...")
   build_equivalence_classes!( quant, lib, assign_long=true )
   println(STDERR, "- $( length(quant.classes) ) multi-isoform equivalence classes...")
   calculate_tpm!( quant, readlen=readlen )
   @timer iter = gene_em!( quant, multi, sig=1, readlen=readlen, maxit=10000 )
   println(STDERR, "Finished calculating transcripts per million (TpM) after $iter iterations of EM...")
   set_gene_tpm!( quant, lib )
   gc_normalize!( mod, lib, quant )
   gc_adjust!( quant, mod )
   gc_adjust!( multi, mod )

   output_tpm( args["out"], lib, quant )
   output_stats( args["out"] * ".map.gz", lib, quant, param, indexpath, totreads, mapped, Int(round(total_multi(multi))), readlen, ver )
   output_junctions( args["out"] * ".jnc.gz", lib, quant )

   println(STDERR, "Assigning multi-mapping reads based on maximum likelihood estimate..")
   # Now assign multi to edges.
   @timer assign_ambig!( quant, lib, multi )

   effective_lengths!( lib, quant, 1, 0)
   println(STDERR, "Calculating maximum likelihood estimate of events..." )
   @timer process_events( args["out"] * ".psi.gz" , lib, quant, isnodeok=false, iscircok=args["circ"], readlen=readlen )
   println(STDERR, "Whippet $ver done." )
end

@timer main()


text:
export setseed, seed!
import Random: rand!, randn!, Random
using Knet.KnetArrays: KnetArray

rand!(a::KnetArray)=(rand!(CuArray(a)); a)
randn!(a::KnetArray)=(randn!(CuArray(a)); a)

function setseed(x)
    @warn "setseed() is deprecated, please use Random.seed!() and/or CUDA.seed!() instead" maxlog=1
    CUDA.functional() && CUDA.seed!(x)
    Random.seed!(x)
end

function seed!(x)
    @warn "Knet.seed!() is deprecated, please use Random.seed!() and/or CUDA.seed!() instead" maxlog=1
     CUDA.functional() && CUDA.seed!(x)
     Random.seed!(x)
end


text:
# runtests.jl

for s in readdir("."; join=true, sort=true)
    f, e = splitext(s)
    if e == ".pdf" || e == ".vtu"
  try
   rm(s)
        catch
            @warn "Failed to remove $s"
        end
    end
end

for s in [
    "LE5_Z_cantilever_examples",
    "convergence_metaanalysis",
    "cos_2t_press_hyperboloid_free_examples",
    "vis_scordelis_lo_examples",
    "barrel_vault_examples",
    "cos_2t_press_cylinder_fixed_examples",
    "hemisphere_open_examples",
    "scordelis_lo_examples",
    "barrel_w_stiffeners_examples",
    "cos_2t_press_cylinder_free_examples",
    "pinched_cylinder_examples",
    "clamped_hypar_examples",
    "cos_2t_press_hyperboloid_fixed_examples",
    "raasch_examples",
    "ss_circular_plate_udl_examples",
    ]
    include(s * ".jl")
end


text:
abstract type Variable end

mutable struct Parameter <: Variable
  name::String
  info::String
  min::Float64
  max::Float64
  init::Float64
  units::String
  constant::Bool
  bounds::Tuple{Float64, Float64}
  # Post fit information
  samplewidth::Float64
  efficiencies::Dict{Symbol, Float64}
  fit::Float64
  low::Float64
  high::Float64
  likelihood_x::Array{Float64}
  likelihood_y::Array{Float64}
  function Parameter( name::String; kwargs... )
    kwargs = Dict(kwargs)
    info         = get(kwargs, :info, "")
    min          = get(kwargs, :min, -Inf)
    max          = get(kwargs, :max, Inf)
    init         = get(kwargs, :init, 0.0)
    units        = get(kwargs, :units, "")
    constant     = get(kwargs, :constant, false)
    bounds       = get(kwargs, :bounds, (-Inf, Inf) )
    samplewidth  = get(kwargs, :samplewidth, 0.0)
    efficiencies = get(kwargs, :efficiencies, Dict{Symbol, Float64}())
    new( name, info, min, max, init, units, constant, bounds,
         samplewidth, efficiencies,
         0.0, 0.0, 0.0,
         Array{Float64}(undef, 0),
         Array{Float64}(undef, 0),
       )
  end

end

mutable struct Constant <: Variable
  name::String
  info::String
  init::Float64
  units::String
  constant::Bool
  function Constant( name::String, init::Float64;
                     info::String="", units::String="" )
    new( name, info, init, units, true )
  end
end

struct Dataset <: Variable
  name::String
  init::String
  constant::Bool
  #init::Symbol
  function Dataset( name::String ; info::String="" )
    new( name, name, true )
  end
end

mutable struct Observable{T<:Real}
  name::String
  sname::Symbol
  min::T
  max::T
end