text:
struct TrimConditions
    # TODO: improve typing
    tas::Number  # m/s
    ψ::Number  # rad
    position::AbstractArray  # [x, y, z] earth (m)
    ψ_dot::Number  # rad/s
    γ::Number  # rad
    aircraft::Aircraft
    atmosphere::Atmosphere
    gravity::Gravity
end

struct TrimSolution
    # TODO: improve typing
    α::Number
    β::Number
    controls::AbstractArray
    θ::Number
    ϕ::Number
    p::Number
    q::Number
    r::Number
end

function trim(
    dss_guess, controls_guess, aircraft, atmosphere, gravity, γ=0.0, ψ_dot=0.0;
    show_trace=false,
    ftol=1e-16,
    iterations=5000
)

    # TRIMMING SOLUTION
    sol_gues = [
        get_α(dss_guess),  # alpha (rad)
        get_β(dss_guess),  # beta (rad)
        controls_guess..., # thtl  (0-1), el (deg), ail (deg), rdr (deg)
    ]

    # CONSTS
    trim_conditions = TrimConditions(
        get_tas(dss_guess),  # TAS (m/s)
        get_euler_angles(dss_guess)[1],  # psi (rad)
        get_earth_position(dss_guess),  # north, east, down (m)
        ψ_dot,  # ψ_dot (rad/s)
        γ,  # γ (rad)
        aircraft,
        atmosphere,
        gravity,
    )

    ds_type = typeof(dss_guess)

    f_opt(sol) = trim_cost_function(ds_type, sol, trim_conditions; full_output=false)

    result = nlsolve(
        f_opt, sol_gues;
        ftol=ftol, show_trace=show_trace, iterations=iterations
    )

     if show_trace
            println(result)
     end

    # TODO: return a DSStateDot instead of x
    sol = result.zero
    dssd, controls, outputs, cost = trim_cost_function(
        ds_type, sol, trim_conditions; full_output=true
    )
    return dssd, controls, outputs, cost
end

function trim_cost_function(ds::Type{T}, sol, trim_conditions; full_output=false) where {T<:DSState}

    tas = trim_conditions.tas
    ψ = trim_conditions.ψ
    x, y, z = trim_conditions.position
    ψ_dot = trim_conditions.ψ_dot
    γ = trim_conditions.γ
    aircraft = trim_conditions.aircraft
    atmosphere = trim_conditions.atmosphere
    gravity = trim_conditions.gravity

    gd = get_gravity_accel(gravity)

    ϕ, θ, p, q, r = apply_trimmer_constrains(
        tas,  # m/s
        sol[1],  # α (rad)
        sol[2],  # β (rad)
        γ,  # rad
        ψ_dot,  # rad/s
        gd,  # m/s²
    )

    trim_solution = TrimSolution(
        sol[1],  # α (rad)
        sol[2],  # β (rad)
        sol[3:6],  # controls  -> thtl  (0-1), el (deg), ail (deg), rdr (deg)
        θ,  # rad
        ϕ,  # rad
        p,  # rad/s
        q,  # rad/s
        r,  # rad/s
    )
    # Construct state vector
    dss = ds(trim_conditions, trim_solution, aircraft)
    dssd, outputs = f(time, dss, trim_solution.controls, aircraft, atmosphere, gravity)

    x_dot = get_xdot(dssd)
    cost = [x_dot[1:3]..., x_dot[7:9]...]

    if full_output
        return dssd, trim_solution.controls, outputs, cost
    else
        return cost
    end
end

function apply_trimmer_constrains(tas, α, β, γ, ψ_dot, gd)
    # Coordinated turn bank --> phi
    ϕ = coordinated_turn_bank(ψ_dot, α, β, tas, γ, gd)
    # Climb -> theta
    θ = rate_of_climb_constrain_no_wind(γ, α, β, ϕ)
    # Angular kinemtic -> p, q, r
    p, q, r = ψθϕ_dot_2_pqr(ψ_dot, 0, 0, θ, ϕ)

    return ϕ, θ, p, q, r
end

function SixDOFAeroEuler{T}(tc::TrimConditions, ts::TrimSolution, ac::Aircraft) where {T}
    x = [
        tc.tas,
        ts.α,
        ts.β,
        ts.ϕ,
        ts.θ,
        tc.ψ,
        ts.p,
        ts.q,
        ts.r,
        tc.position...,
        tgear(ac, ts.controls[1])
    ]
    return SixDOFAeroEuler(x)
end

function SixDOFBodyEuler{T}(tc::TrimConditions, ts::TrimSolution, ac::Aircraft) where {T}
    x = [
        wind2body(tc.tas, 0, 0, ts.α, ts.β)...,
        ts.ϕ,
        ts.θ,
        tc.ψ,
        ts.p,
        ts.q,
        ts.r,
        tc.position...,
        tgear(ac, ts.controls[1])
    ]
    return SixDOFBodyEuler(x)
end


text:
#### Document type

struct _XMLDocStruct  # use the same layout as C
    # common part
    _private::Ptr{Cvoid}
    nodetype::Cint
    name::Xstr
    children::Xptr
    last::Xptr
    parent::Xptr
    next::Xptr
    prev::Xptr
    doc::Xptr

    # specific part
    compression::Cint
    standalone::Cint
    intsubset::Xptr
    extsubset::Xptr
    oldns::Xptr
    version::Xstr
    encoding::Xstr

    ids::Ptr{Cvoid}
    refs::Ptr{Cvoid}
    url::Xstr
    charset::Cint
    dict::Xstr
    psvi::Ptr{Cvoid}
    parseflags::Cint
    properties::Cint
end

mutable struct XMLDocument
    ptr::Xptr
    _struct::_XMLDocStruct

    function XMLDocument(ptr::Xptr)
        s::_XMLDocStruct = unsafe_load(convert(Ptr{_XMLDocStruct}, ptr))

        # validate integrity
        @assert s.nodetype == XML_DOCUMENT_NODE
        @assert s.doc == ptr

        new(ptr, s)
    end

    function XMLDocument()
        # create an empty document
        ptr = ccall((:xmlNewDoc,libxml2), Xptr, (Cstring,), "1.0")
        XMLDocument(ptr)
    end
end

version(xdoc::XMLDocument) = unsafe_string(xdoc._struct.version)
encoding(xdoc::XMLDocument) = unsafe_string(xdoc._struct.encoding)
compression(xdoc::XMLDocument) = Int(xdoc._struct.compression)
standalone(xdoc::XMLDocument) = Int(xdoc._struct.standalone)

function root(xdoc::XMLDocument)
    pr = ccall((:xmlDocGetRootElement,libxml2), Xptr, (Xptr,), xdoc.ptr)
    pr != C_NULL || throw(XMLNoRootError())
    XMLElement(pr)
end

#### construction & free

function free(xdoc::XMLDocument)
    ccall((:xmlFreeDoc,libxml2), Cvoid, (Xptr,), xdoc.ptr)
    xdoc.ptr = C_NULL
end

function set_root(xdoc::XMLDocument, xroot::XMLElement)
    ccall((:xmlDocSetRootElement,libxml2), Xptr, (Xptr, Xptr), xdoc.ptr, xroot.node.ptr)
end

function create_root(xdoc::XMLDocument, name::AbstractString)
    xroot = new_element(name)
    set_root(xdoc, xroot)
    return xroot
end

#### parse and free

function _check_result(p)
    p != C_NULL || throw(XMLParseError("Failure in parsing an XML file."))
    XMLDocument(p)
end

parse_file(filename::AbstractString) =
    _check_result(ccall((:xmlParseFile,libxml2), Xptr, (Cstring,), filename))

parse_file(filename::AbstractString, encoding, options::Integer) =
    _check_result(ccall((:xmlReadFile,libxml2), Xptr, (Cstring, Ptr{Cchar}, Cint),
        filename, encoding, options))

parse_string(s::AbstractString) =
    _check_result(ccall((:xmlParseMemory,libxml2), Xptr, (Xstr, Cint), s, sizeof(s) + 1))

#### output

function save_file(xdoc::XMLDocument, filename::AbstractString; encoding::AbstractString="utf-8")
    ret = ccall((:xmlSaveFormatFileEnc,libxml2), Cint,
                (Cstring, Xptr, Cstring, Cint),
                filename, xdoc.ptr, encoding, 1)
    ret < 0 && throw(XMLWriteError("Failed to save XML to file $filename"))
    Int(ret)  # number of bytes written
end

function Base.string(xdoc::XMLDocument; encoding::AbstractString="utf-8")
    buf_out = Vector{Xstr}(undef, 1)
    len_out = Vector{Cint}(undef, 1)
    ccall((:xmlDocDumpFormatMemoryEnc,libxml2), Cvoid,
          (Xptr, Ptr{Xstr}, Ptr{Cint}, Cstring, Cint),
          xdoc.ptr, buf_out, len_out, encoding, 1)
    _xcopystr(buf_out[1])
end

Base.show(io::IO, xdoc::XMLDocument) = print(io, string(xdoc))

ext: jl


text:
# Calculates and plots the learning/test curves as we increase the number of grids points in the volfrac X Radius space.
function learning_saturation_curve_grid_volfrac_rad_crossvalidation(train::Array{StatisticalMoments{Float64}}, test::Array{StatisticalMoments{Float64}}, fm::FeatureMap, options, kernelname::AbstractString, output_type::AbstractString,kernelparams, crossnum,lambdas, numparts, dataname::AbstractString, datatestname::AbstractString)

  pairs = unique([ [td.label[2], td.label[1] ] for td in train]); # pairs of (radius, volume)
  lpr = length(pairs);
  batchsize = convert(Int64,floor(lpr/convert(Float64,numparts)));
  Rsqds = zeros(numparts);
  #Reshuffling the data
  s = sample(1:lpr, lpr,replace=false);
  pairs = pairs[s];   # Reshuffle the data
  # Get best cross valid lambdas and kernelparameter
  lambda, kernelparam = get_crossvalid_parameters(fm,kernelname,output_type, options, dataname);
  for k =1:numparts #Adding on more pairs in each iteration
    trainselection = 1:batchsize*k;  println("Training with ",batchsize*k, " points" )
    if(k==numparts) trainselection = 1:lpr end # get all data on last batch
    trainpairs = pairs[trainselection];
    #Get the testpair coresponding to the crossnum th section
    trainbatch = filter(td -> ([td.label[2],td.label[1]] in trainpairs), train);
    #
    lambda, kernelparam=     cross_valid_inner(trainbatch, fm, options, kernelname, kernelparams,crossnum,lambdas, output_type, false, dataname)
    ml = setup_ml(trainbatch, fm,  kernelname,options, output_type);
    fit_L2(ml, options,  kernelname, kernelparam, lambda);
    Rsqd = test_linear(test,datatestname, fm,  ml,options,output_type);
    Rsqds[k] =Rsqd;
  end
  ml = setup_ml(train, fm,  kernelname,options, output_type);
  savename ="$(ml.name)-learncurve-crossvalid"
  default_path = "../data/"
  println("saving the learning curves in $(default_path)$(savename).jld")
  save("$(default_path)$(savename).jld", "Rsqds",Rsqds)
  # Plotting learning saturation curve
  plot((1:1:numparts)./numparts ,Rsqds,line = :line, xlabel = "fraction of train data",
  ylabel = string(L"$R^2$ of the ", output_type),label = "", ylims = (-1,1), grid = false)
  #grid("off")
  savefig("../figures/$(savename).pdf");
end


text:
# Generates a tablaue of weights and points for Gauss Legendre quadrature
# on a line
using FastGaussQuadrature

const MAX_ORDER = 5
open(joinpath(@__DIR__, "..", "src", "quaddata.jl"), "w") do f
    println(f, "const QUAD_DATA = [")
    for g in 1:MAX_ORDER
        points, weights = gausslegendre(g)
        println(f, "[(", join(points, ", "), ",),")
        println(f, "(", join(weights, ", "), ",)],")
    end
    println(f, "]")
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 4bfcd930-7aa0-11eb-2c61-71d83e5a7efc
begin
	using Images
	using ImageMagick
	using Colors
	#using HypertextLiteral
	vine = load("neko_vignetto.jpg")
end

# ╔═╡ b9d33ef0-7a9f-11eb-3aa6-3d8dc5749a94
begin
	progname = "Julia"
	notebookname = "Pluto Notebook"
end

# ╔═╡ 938783e0-7aa0-11eb-193a-efc9d1424da5
begin
	(h, w) = size(vine)
	wmodify = (w ÷ 2) - 300
	vine_face = vine[1:h, wmodify:w]
end

# ╔═╡ f85d9250-7aa0-11eb-3958-c184816057ae
size(vine)

# ╔═╡ 3db90b80-7aa2-11eb-0448-7d886fdb4367
[
	vine_face reverse(vine_face, dims = 2)
	reverse(vine_face, dims = 1) reverse(reverse(vine_face, dims = 1), dims = 2)
]

# ╔═╡ 71fe20f0-7aa3-11eb-179a-99bfc4165276
function redify(color)
	return RGB(0, color.g, 0)
end

# ╔═╡ 12ff3820-7b78-11eb-0f2f-27f9776a33ef
[RGB(redval, 0.0, blueval) for redval in 0:0.1:1, blueval in 0:0.1:1]

# ╔═╡ 8d443ed0-7aa3-11eb-2888-0920888050e6
begin
	takecolor = RGB(0.98, 0.15, 0.67)
	[takecolor, redify.(takecolor)]
end

# ╔═╡ ab3afdc0-7aa3-11eb-16ab-177c4031bfdd
redify.(vine_face)

# ╔═╡ b78d7b60-7aa4-11eb-204f-cbc675fcd835
function monochrome(color)
	RGB(color.r, color.g, color.b)
	if color.b < 0.50
		return RGB(0, 0, 0)
	else
		return RGB(1, 1, 1)
	end
end

# ╔═╡ 70693e30-7aa5-11eb-0b7e-e553fd22275b
monochrome.(vine_face)

# ╔═╡ 9d4db480-7aa5-11eb-1a9e-21c7c6960b5f
function grayscale(color)
	RGB(color.r, color.g, color.b)
	μ = (color.r + color.g + color.b) / 3
	return RGB(μ, μ, μ)
end

# ╔═╡ e6106a00-7aa5-11eb-275b-e97981411116
grayscale.(vine_face)

# ╔═╡ 81da4a60-7b7c-11eb-2ece-a73c0c7b4498
rand(5)

# ╔═╡ 6d60c900-7b7d-11eb-20c9-2f0d414e96e6
function inverse(color)
	RGB(color.r, color.g, color.b)
	return RGB(1 - color.r, 1 - color.g, 1 - color.b)
end

# ╔═╡ c525a840-7b7d-11eb-1e4d-ab0e27dce733
inverse.(vine_face)

# ╔═╡ 49ef6c30-7b80-11eb-1266-13c00983662e
function distgauss(x, μ, σ)
	return (1 / (sqrt(2π) * σ)) * exp(-(1 / 2) * ((x - μ) / σ)^2)
end

# ╔═╡ 4f9ce470-7b7e-11eb-3fb7-4707db8ecee2
function noise(color)
	RGB(color.r, color.g, color.b)
	η = rand(3)
	return RGB(color.r + η[1], color.g + η[1], color.b + η[1])
end

# ╔═╡ 837371e0-7b85-11eb-0f9c-dbb24f6f7ab3

# ╔═╡ 75126e50-7b7e-11eb-0a8f-332434b49591
noise.(vine_face)

# ╔═╡ Cell order:
# ╠═b9d33ef0-7a9f-11eb-3aa6-3d8dc5749a94
# ╠═4bfcd930-7aa0-11eb-2c61-71d83e5a7efc
# ╠═938783e0-7aa0-11eb-193a-efc9d1424da5
# ╠═f85d9250-7aa0-11eb-3958-c184816057ae
# ╠═3db90b80-7aa2-11eb-0448-7d886fdb4367
# ╠═71fe20f0-7aa3-11eb-179a-99bfc4165276
# ╠═12ff3820-7b78-11eb-0f2f-27f9776a33ef
# ╠═8d443ed0-7aa3-11eb-2888-0920888050e6
# ╠═ab3afdc0-7aa3-11eb-16ab-177c4031bfdd
# ╠═b78d7b60-7aa4-11eb-204f-cbc675fcd835
# ╠═70693e30-7aa5-11eb-0b7e-e553fd22275b
# ╠═9d4db480-7aa5-11eb-1a9e-21c7c6960b5f
# ╠═e6106a00-7aa5-11eb-275b-e97981411116
# ╠═81da4a60-7b7c-11eb-2ece-a73c0c7b4498
# ╠═6d60c900-7b7d-11eb-20c9-2f0d414e96e6
# ╠═c525a840-7b7d-11eb-1e4d-ab0e27dce733
# ╠═49ef6c30-7b80-11eb-1266-13c00983662e
# ╠═4f9ce470-7b7e-11eb-3fb7-4707db8ecee2
# ╠═837371e0-7b85-11eb-0f9c-dbb24f6f7ab3
# ╠═75126e50-7b7e-11eb-0a8f-332434b49591

max_stars_repo_name: m-RezaFahlevi/Mathematica-Works

max_issues_repo_name: m-RezaFahlevi/Mathematica-Works

max_forks_repo_name: m-RezaFahlevi/Mathematica-Works


text:
isdefined(Base, :__precompile__) && __precompile__()

"""
Module for constructing self-contained power system objects.
"""
module PowerSystems

#################################################################################
# Exports

export System
export Bus
export Arc
export AggregationTopology
export Area
export LoadZone
export get_aggregation_topology_accessor

export PowerSystemType
export Component
export Device
export Branch
export StaticInjection
export ACBranch
export Line
export MonitoredLine
export DCBranch
export HVDCLine
export VSCDCLine
export Transformer2W
export TapTransformer
export PhaseShiftingTransformer

export ThreePartCost
export TwoPartCost
export VariableCost
export MultiStartCost
export get_slopes
export get_breakpoint_upperbounds

export Generator
export HydroGen
export HydroDispatch
export HydroEnergyReservoir

export RenewableGen
export RenewableFix
export RenewableDispatch

export ThermalGen
export ThermalStandard
export ThermalMultiStart

export ElectricLoad
export StaticLoad
export PowerLoad
export PowerLoadPF
export FixedAdmittance
export ControllableLoad
export InterruptibleLoad

export Storage
export GenericBattery

export DynamicInjection
export DynamicGenerator
export RegulationDevice

#AVR Exports
export AVR
export AVRFixed
export AVRSimple
export AVRTypeI
export AVRTypeII

#Machine Exports
export Machine
export BaseMachine
export RoundRotorMachine
export SalientPoleMachine
export RoundRotorQuadratic
export SalientPoleQuadratic
export RoundRotorExponential
export SalientPoleExponential
export OneDOneQMachine
export MarconatoMachine
export SimpleMarconatoMachine
export AndersonFouadMachine
export SimpleAFMachine
export FullMachine
export SimpleFullMachine

#PSS Exports
export PSS
export PSSFixed
export PSSFixed

#Shaft Exports
export SingleMass
export FiveMassShaft

#TG Exports
export TurbineGov
export TGFixed
export TGTypeI
export TGTypeII

export DynamicInverter
# Converter Exports
export Converter
export AverageConverter

# DC Source Exports
export DCSource
export FixedDCSource

# Filter Exports
export Filter
export LCLFilter

# FrequencyEstimator Exports
export FrequencyEstimator
export KauraPLL

# Outer Control Exports
export OuterControl
export VirtualInertia
export ReactivePowerDroop

# InnerControl Export
export InnerControl
export CurrentControl

export Source

export Service
export Reserve
export ReserveDirection
export ReserveUp
export ReserveDown
export StaticReserve
export VariableReserve
export AGC
export ReserveDemandCurve
export Transfer

export PTDF
export Ybus
export LODF
export GeneratorCostModels
export AngleUnits
export BusTypes
export LoadModels
export PrimeMovers
export ThermalFuels

export Forecast
export Deterministic
export Probabilistic
export ScenarioBased
export PiecewiseFunction

export get_dynamic_components

export solve_powerflow!

export parse_file
export add_forecasts!
export add_forecast!
export remove_forecast!
export clear_forecasts!
export add_component!
export remove_component!
export remove_components!
export clear_components!
export add_service!
export remove_service!
export clear_services!
export has_service
export has_forecasts
export get_buses
export get_components_in_aggregation_topology
export get_aggregation_topology_mapping
export get_contributing_devices
export get_contributing_device_mapping
export ServiceContributingDevices
export ServiceContributingDevicesKey
export ServiceContributingDevicesMapping
export are_forecasts_contiguous
export generate_initial_times
export get_component
export get_components
export get_components_by_name
export get_available_components
export get_forecast_labels
export get_forecast_initial_times
export get_forecast
export get_forecast_values
export get_forecasts_horizon
export get_forecasts_initial_time
export get_forecasts_interval
export get_forecasts_resolution
export get_horizon
export get_initial_time
export get_resolution
export get_data
export iterate_components
export iterate_forecasts
export make_forecasts
export get_bus_numbers
export get_name
export get_basepower
export get_frequency
export to_json
export check_forecast_consistency
export validate_forecast_consistency
export clear_ext
export convert_component!
export set_area!
export set_load_zone!
export TamuSystem
export PowerModelsData
export set_dynamic_injector!
export set_static_injector!
export get_V_ref
export get_P_ref
export get_Q_ref
export get_saturation_coeffs
export set_droop!
export set_participation_factor!
export set_inertia!
export set_reserve_limit_up!
export set_reserve_limit_dn!
export set_cost!
export get_droop
export get_inertia
export get_reserve_limit_up
export get_reserve_limit_dn
export get_participation_factor
export get_cost

#################################################################################
# Imports

import SparseArrays
import LinearAlgebra: LAPACK.getri!
import LinearAlgebra: LAPACK.getrf!
import LinearAlgebra: BLAS.gemm
import LinearAlgebra
import Dates
import TimeSeries
import DataFrames
import JSON
import JSON2
import CSV
import YAML
import UUIDs
import Base.to_index

import InfrastructureSystems
import InfrastructureSystems:
    Components,
    Deterministic,
    Probabilistic,
    Forecast,
    ScenarioBased,
    PiecewiseFunction,
    InfrastructureSystemsType,
    InfrastructureSystemsInternal,
    FlattenIteratorWrapper,
    LazyDictFromIterator,
    DataFormatError,
    InvalidRange,
    InvalidValue,
    get_data,
    get_horizon,
    get_initial_time,
    get_resolution,
    get_name

const IS = InfrastructureSystems

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

using DocStringExtensions

@template (FUNCTIONS, METHODS) = """
                                 $(SIGNATURES)
                                 $(DOCSTRING)
                                 """

#################################################################################
# Includes

"""
Supertype for all PowerSystems types.
All subtypes must include a InfrastructureSystemsInternal member.
Subtypes should call InfrastructureSystemsInternal() by default, but also must
provide a constructor that allows existing values to be deserialized.
"""
abstract type PowerSystemType <: IS.InfrastructureSystemsType end

abstract type Component <: PowerSystemType end
# supertype for "devices" (bus, line, etc.)
abstract type Device <: Component end
abstract type DeviceParameter <: PowerSystemType end

include("common.jl")
include("models/static_models.jl")
include("models/dynamic_models.jl")
include("models/injection.jl")

# Include utilities
include("utils/IO/base_checks.jl")

# PowerSystems models
include("models/topological_elements.jl")
include("models/branches.jl")
include("models/operational_cost.jl")
#include("models/network.jl")

# Static types
include("models/services.jl")
include("models/reserves.jl")
include("models/generation.jl")
include("models/storage.jl")
include("models/loads.jl")
include("models/regulation_device.jl")
include("models/dynamic_generator_components.jl")
include("models/dynamic_inverter_components.jl")

# Include all auto-generated structs.
include("models/generated/includes.jl")

#Methods for devices
include("models/devices.jl")

# Dynamic Composed types
include("models/dynamic_generator.jl")
include("models/dynamic_inverter.jl")
include("models/dynamic_machines.jl")

include("models/supplemental_constructors.jl")
include("models/supplemental_accessors.jl")

# Definitions of PowerSystem
include("base.jl")

#Data Checks
include("utils/IO/system_checks.jl")
include("utils/IO/branchdata_checks.jl")

# network calculations
include("utils/network_calculations/common.jl")
include("utils/network_calculations/ybus_calculations.jl")
include("utils/network_calculations/ptdf_calculations.jl")
include("utils/network_calculations/lodf_calculations.jl")

#PowerFlow
include("utils/power_flow/power_flow.jl")

# Include Parsing files
include("parsers/common.jl")
include("parsers/enums.jl")
include("parsers/pm_io.jl")
include("parsers/im_io.jl")
include("parsers/power_system_table_data.jl")
include("parsers/power_models_data.jl")
include("parsers/TAMU_data.jl")

# Better printing
include("utils/print.jl")

include("models/serialization.jl")

# Download test data
include("utils/data.jl")
import .UtilsData: TestData

end # module


text:
module _XPermEn
export XPermEn
using StatsBase: fit, Histogram
using Combinatorics: permutations

    """
        XPerm = XPermEn(Sig)

    Returns the cross-permuation entropy estimates (`XPerm`) estimated betweeen
    the data sequences contained in `Sig` using the default parameters:
    embedding dimension = 3, time delay = 1, logarithm = base 2,

        XPerm = XPermEn(Sig::AbstractArray{T,2} where T<:Real; m::Int=3, tau::Int=1, Logx::Real=exp(1))

    Returns the permutation entropy estimates (`XPerm`) estimated between the
    data sequences contained in `Sig` using the specified 'keyword' arguments:

    # Arguments:
    `m`     - Embedding Dimension, an integer > 2   [default: 3]  \n
            **Note: XPerm is undefined for embedding dimensions < 3.**\n
    `tau`   - Time Delay, a positive integer        [default: 1]    \n
    `Logx`  - Logarithm base, a positive scalar     [default: 2]
            ** enter 0 for natural log.**    \n

    # See also `PermEn`, `XApEn`, `XSampEn`, `XFuzzEn`, `XMSEn`

    # References:
        [1] Wenbin Shi, Pengjian Shang, and Aijing Lin,
            "The coupling analysis of stock market indices based on
            cross-permutation entropy."
            Nonlinear Dynamics
            79.4 (2015): 2439-2447.

    """
    function XPermEn(Sig::AbstractArray{T,2} where T<:Real; m::Int=3, tau::Int=1, Logx::Real=exp(1))

    (size(Sig,2) > size(Sig,1)) ? Sig = transpose(Sig) : nothing

    N = size(Sig,1)
    (N>10 && size(Sig,2)==2) ? nothing :  error("Sig:   must be a 2-columns matrix")
    (m > 2) ? nothing : error("m:     must be an integer > 1")
    (tau>0) ? nothing : error("tau:   must be an integer > 0")
    (Logx>0) ? nothing : error("Logx:     must be a positive number > 0")

    S1 = Sig[:,1]; S2 = Sig[:,2];
    N = length(S1)-(m-1)*tau
    Sx1 = zeros(N,m)
    Sx2 = zeros(N,m)
    for k = 1:m
        Sx1[:,k] = S1[1+(k-1)*tau:N+(k-1)*tau]
        Sx2[:,k] = S2[1+(k-1)*tau:N+(k-1)*tau]
    end

    Temp = sortind(Sx1[1:N,1:m])
    Gx = zeros(N,m)
    for k = 1:N
        Gx[k,:] = Sx2[k,Temp[k,:]]
    end

    Kt = zeros(m-2,m-2,N)
    for k = 1:m-2
        for j = k+1:m-1
            G1 = Gx[:,j+1] .- Gx[:,k]
            G2 = Gx[:,k] .- Gx[:,j]
            Kt[k,j-1,:] = (G1.*G2 .> 0)
        end
    end

    Di = sum(Kt,dims=(1,2))[:]
    Ppi = fit(Histogram, Di, -.5:((m-2)*(m-1) + 1)/2).weights/N
    Ppi = Ppi[Ppi.!=0]
    XPerm = -sum(Ppi.*log.(Logx,Ppi))
    if round(sum(Ppi),digits=6)!=1
        @warn("Potential error with probability calculation")
    end

    return XPerm
    end

    function sortind(X)
        Y = zeros(Int, size(X))
        for k = 1:length(X[:,1])
            Y[k,:] = sortperm(X[k,:])
        end
        return Y
    end
end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
"""
    predict(y, x, xb, xt, hyp, ∇)

    predictive distribution for SPGP given hyperparameters and pseudo-inputs

    y -- training targets (N x 1)
    x -- training inputs (N x dim)
    xb -- pseudo inputs (n x dim)
    xt -- test inputs (Nt x dim)
    hyp -- hyperparameters (including noise)
           for Gaussian covariance: (dim+2 x 1)

          hyp(1:dim) = log( b )
          hyp(dim+1) = log( c )
          hyp(dim+2) = log( sig )

         where cov = c * exp[-0.5 * sum_d b_d * (x_d - x'_d)^2]
                         + sig * delta(x,x')

    ∇ -- OPTIONAL jitter (default 1e-6)
    μ -- predictive mean
    σ² -- predictive variance of latent function
          (add noise if full variance is required)
"""
function predict(y, x, xb, xt, hyp, ∇=1e-6)
    n_observ, n_dim = size(x, 1), size(x, 2)
    n_pseudo = size(xb, 1)
    n_train = size(xt, 1)

    σ = exp(hyp[end])

    K = kernel(xb, xb, hyp) + ∇ * I
    L = cholesky(Hermitian(K)).L
    K = kernel(xb, x, hyp)
    V = L \ K

    @inline function kdiag(x, hyp)
        c = exp(hyp[end - 1])
        Kd = fill(c, size(x, 1), 1)

        return Kd
    end

    ɛ = 1 .+ (kdiag(x, hyp) - transpose(sum(V.^2, dims=1))) / σ

    V = V ./ repeat(transpose(sqrt.(ɛ)), n_pseudo, 1);
    y = y ./ sqrt.(ɛ);
    Lm = transpose(cholesky(σ * I + V * transpose(V)).U);
    bet = Lm \ (V * y);

    K = kernel(xb, xt, hyp);
    lst = L \ K;
    lmst = Lm \ lst;
    μ = transpose((transpose(bet) * lmst));

    σ² = kdiag(xt, hyp) .- transpose(sum(lst.^2, dims=1)) .+ σ .* transpose(sum(lmst.^2, dims=1));

    return μ, σ²
end


text:
# This file includes test set for modulator

@testset "ModulatorTestset" begin

# Fields
@test :scheme in fieldnames(Modulator)
@test :Ep in fieldnames(Modulator)

# Construction
modulator = Modulator(PSK(4), 10.)
modulator = Modulator(ASK(4))
@test modulator.Ep == 1.

# Modulation
M = 4
bits = [0, 0, 1, 1, 0, 1, 1, 0, 0, 1]
coding = GrayCoding(M)
modulator = Modulator(PSK(M))
tx = bits |> coding |> modulator
@test tx[1] ≈ [1, 0] / sqrt(2)
@test tx[2] ≈ [-1, 0] / sqrt(2)
@test tx[3] ≈ [0, 1] / sqrt(2)
@test tx[4] ≈ [0, -1] / sqrt(2)
@test tx[5] ≈ [0, 1] / sqrt(2)

end # testset


text:
# ------------------------------------------------------------------- #
# CES Producer
# ------------------------------------------------------------------- #
immutable CESProducer <: AbstractProducer
    η::Float64
    ν::Float64
end

# call(p::CESProducer, k, l) = (p.η*k.^p.ν + (1-p.η)*l.^(p.ν)).^(1./p.ν)
Base.produce(p::CESProducer, k, l) = (p.η*k.^p.ν + (1-p.η)*l.^(p.ν)).^(1./p.ν)
f_k(p::CESProducer, k, l) = p.η * (produce(p, k, l) ./ k).^(1.-p.ν)

# ------------------------------------------------------------------- #
# Cobb-Douglas Producer
# ------------------------------------------------------------------- #
immutable CobbDouglassProducer <: AbstractProducer
    α_k::Float64
    α_l::Float64
end

is_crs(i::CobbDouglassProducer) = i.α_k == 1. - i.α_l

max_stars_repo_name: NYUEcon/GrowthModels

max_issues_repo_name: NYUEcon/GrowthModels

max_forks_repo_name: NYUEcon/GrowthModels


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

using Markdown
using InteractiveUtils

# ╔═╡ 211f1a10-061d-452a-b061-aa5eb5f07428
begin
    using Statistics
    using DataFrames
    using CSV
    using Random
	using PlantBiophysics
end

# ╔═╡ 994a4062-39cf-49b7-b19e-63ff90e9a574
md"""

# _PlantBiophysics.jl_ benchmark

The main objective of this notebook is to compare the computational times of `PlantBiophysics.jl` against the `plantecophys` R package and the `LeafGasExchange.jl` Julia package from the `Cropbox.jl` framework. The comparison follows three steps:
- create an N-large basis of random conditions.
- benchmark the computational time of the three packages via similar functions (_i.e._ photosynthesis-stomatal conductance-energy balance coupled model for C3 leaves): `energy_balance`, `photosynEB` and `simulate` with `ModelC3MD`.
- compare the results with plots and statistics.

For running the models you'll need to use the following packages:
```julia
begin
    using CairoMakie
    using BenchmarkTools
	using PlantBiophysics
    using Cropbox
    using LeafGasExchange
    using RCall
end
```

!!! note
	You need to have R installed on your system to use RCall.

!!! warning
	This notebook does not perform the benchmark by itself for the obvious reason that it takes forever to run, and because there is an overhead cause by Pluto (benchmark and reactive are not a good mix). Instead, it shows the outputs of [a script from this repository](https://github.com/VEZY/PlantBiophysics-paper/blob/main/tutorials/PlantBiophysics_performance_noPluto.jl) that implements the code shown here. If you want to perform the benchmark by yourself, you can run this script from the command line. The versions used for the above dependencies are available in the [Project.toml](https://github.com/VEZY/PlantBiophysics-paper/blob/main/tutorials/Project.toml) of the repository.
"""

# ╔═╡ 9c3f432b-4d1f-4765-9989-17e983cee5f7
md"""

## Parameters

### Benchmark parameters

You'll find below the main parameters of the benchmark. In few words, each package runs a simulation for `N` different time-steps `microbenchmark_steps` times repeated `microbenchmark_evals` times. We make `N` different simulations because the simulation duration can vary depending on the inputs due to iterative computations in the code, *i.e.* different initial conditions can make the algorithms converge more or less rapidly.
"""

# ╔═╡ e7e23262-1813-4fe5-80fb-2be20cb54f89
begin
    Random.seed!(1) # Set random seed
    microbenchmark_steps = 100 # Number of times the microbenchmark is run
    microbenchmark_evals = 1 # N. times each sample is run to be sure of the output
    N = 100 # Number of timesteps simulated for each microbenchmark step
end

# ╔═╡ 65d00bf6-8f9a-40b4-947a-fd46ddc89753
md"""

### Random input simulation dataset

We create possible ranges for input parameters. These ranges where chosen so all of the three packages don't return errors during computation (plantecophys has issues with low temperatures).

- Ta: air temperature ($°C$)
- Wind: wind speed ($m.s^{-1}$)
- P: ambient pressure ($kPa$)
- Rh: relative humidity (between 0 and 1)
- Ca: air CO₂ concentration ($ppm$)
- Jmax: potential rate of electron transport ($\mu mol_{CO2}.m^{-2}.s^{-1}$)
- Vmax: maximum rate of Rubisco activity ($\mu mol_{CO2}.m^{-2}.s^{-1}$)
- Rd: mitochondrial respiration in the light at reference temperature ($\mu mol_{CO2}.m^{-2}.s^{-1}$)
- TPU: triose phosphate utilization-limited photosynthesis rate ($\mu mol_{CO2}.m^{-2}.s^{-1}$)
- Rs: short-wave net radiation ($W.m^{-1}$)
- skyF: Sun-visible fraction of the leaf (between 0 and 1)
- d: characteristic length ($m$)
- g0: residual stomatal conductance ($mol_{CO2}.m^{-2}.s^{-1}$)
- g1: slope of the stomatal conductance relationship.
"""

# ╔═╡ acc685f9-7ae8-483d-982e-ff45ccd9e860
begin
    # Create the ranges of input parameters
    length_range = 10000
    Rs = range(10, 500, length=length_range)
    Ta = range(18, 40, length=length_range)
    Wind = range(0.5, 20, length=length_range)
    P = range(90, 101, length=length_range)
    Rh = range(0.1, 0.98, length=length_range)
    Ca = range(360, 900, length=length_range)
    skyF = range(0.0, 1.0, length=length_range)
    d = range(0.001, 0.5, length=length_range)
    Jmax = range(200.0, 300.0, length=length_range)
    Vmax = range(150.0, 250.0, length=length_range)
    Rd = range(0.3, 2.0, length=length_range)
    TPU = range(5.0, 20.0, length=length_range)
    g0 = range(0.001, 2.0, length=length_range)
    g1 = range(0.5, 15.0, length=length_range)
    vars = hcat([Ta, Wind, P, Rh, Ca, Jmax, Vmax, Rd, Rs, skyF, d, TPU, g0, g1])
    nothing
end

# ╔═╡ b8126c39-8c6c-490d-9992-8922e88f8857
md"We then sample `N` conditions from the given ranges:"

# ╔═╡ a0476d0d-64d9-4457-b78d-41519e38e859
begin
    set = [rand.(vars) for i in 1:N]
    set = reshape(vcat(set...), (length(set[1]), length(set)))'
    name = ["T", "Wind", "P", "Rh", "Ca", "JMaxRef", "VcMaxRef", "RdRef", "Rs", "sky_fraction", "d", "TPURef", "g0", "g1"]
    set = DataFrame(set, name)
    @. set[!, :vpd] = e_sat(set.T) - vapor_pressure(set.T, set.Rh)
    @. set[!, :PPFD] = set.Rs * 0.48 * 4.57
    set
end

# ╔═╡ ae00e44e-d864-4d04-b2e8-fac510ce44bb
md"
## 2. Benchmarking
"

# ╔═╡ 8419c830-8e35-493b-91b2-adfc97e30a61
md"""
##### plantecophys

Preparing R to make the benchmark:

```julia

R\"\"\"
if(!require("plantecophys")){
	install.packages("plantecophys", repos = "https://cloud.r-project.org")
}
if(!require("microbenchmark")){
	install.packages("microbenchmark", repos = "https://cloud.r-project.org")
}
\"\"\"

# Make variables available to the R session
@rput set N microbenchmark_steps

```

Making the benchmark:

```julia
R\"\"\"
	# Define the function call in a function that takes a list as input to limit DataFrame overhead
	function_EB <- function(input) {
		PhotosynEB(
			Tair = input$Tair, VPD = input$VPD, Wind = input$Wind,
			Wleaf = input$Wleaf,Ca = input$Ca,  StomatalRatio = 1,
			LeafAbs = input$LeafAbs, gsmodel = "BBOpti", g0 = input$g0, g1 = input$g1,
			alpha = 0.24, theta = 0.7, Jmax = input$Jmax,
			Vcmax = input$Vcmax, TPU = input$TPU, Rd = input$Rd,
			RH = input$RH, PPFD=input$PPFD, Patm = input$Patm
		)
	}

	time_PE = c()
	for(i in seq_len(N)){
		# Put the inputs into a vector to limit dataframe overhead:
		input = list(
			Tair = set$T[i], VPD = set$vpd[i], Wind = set$Wind[i], Wleaf = set$d[i],
			Ca = set$Ca[i], LeafAbs = set$sky_fraction[i], g0 = set$g0[i], g1 = set$g1[i],
			Jmax = set$JMaxRef[i], Vcmax = set$VcMaxRef[i], TPU = set$TPURef[i],
			Rd = set$RdRef[i], RH = set$Rh[i]*100, PPFD=set$PPFD[i],Patm = set$P[i]
		)

		m = microbenchmark(function_EB(input), times = microbenchmark_steps)

		time_PE = append(time_PE,m$time * 10e-9) # transform in seconds
	}
\"\"\"

@rget time_PE
```

"""

# ╔═╡ 4dad9b75-3a03-4da1-9aa8-fea02fb48a97
md"""

##### LeafGasExchange.jl

Note that we benchmark `LeafGasExchange.jl` with the `nounit` flag to make a faire comparison with `PlantBiophysics.jl` in case computing units takes time (it shouldn't much).

```julia
time_LG = []
n_lg = fill(0, N)
for i in 1:N
	config = :Weather => (
		PFD=set.PPFD[i],
		CO2=set.Ca[i],
		RH=set.Rh[i] * 100,
		T_air=set.T[i],
		wind=set.Wind[i],
		P_air=set.P[i],
		g0=set.g0[i],
		g1=set.g1[i],
		Vcmax=set.VcMaxRef[i],
		Jmax=set.JMaxRef[i],
		Rd=set.RdRef[i],
		TPU=set.TPURef[i]
	)
	b_LG = @benchmark simulate($ModelC3MD; config=$config) evals = microbenchmark_evals samples = microbenchmark_steps
	append!(time_LG, b_LG.times .* 1e-9) # transform in seconds
	n_lg[i] = 1
end
```

"""

# ╔═╡ 60f769b9-bf8c-461c-9739-f31ca43b27b8
md"""
##### PlantBiophysics.jl

Benchmarking `PlantBiophysics.jl`:

```julia
time_PB = []
for i in 1:N
	leaf = LeafModels(
		energy=Monteith(),
		photosynthesis=
		Fvcb(
			VcMaxRef=set.VcMaxRef[i],
			JMaxRef=set.JMaxRef[i],
			RdRef=set.RdRef[i],
			TPURef=set.TPURef[i]
		),
		stomatal_conductance=Medlyn(set.g0[i], set.g1[i]),
		Rₛ=set.Rs[i], sky_fraction=set.sky_fraction[i], PPFD=set.PPFD[i], d=set.d[i]
	)
	meteo = Atmosphere(T=set.T[i], Wind=set.Wind[i], P=set.P[i], Rh=set.Rh[i], Cₐ=set.Ca[i])

	b_PB = @benchmark energy_balance!($leaf, $meteo) evals = microbenchmark_evals samples = microbenchmark_steps
	append!(time_PB, b_PB.times .* 1e-9) # transform in seconds
end
```
"""

# ╔═╡ 87cc52dd-9878-422d-a10e-d491db2a04c6
md"
## 3. Comparison
"

# ╔═╡ 8294a065-5285-4221-8553-e38ca23257f2
md"""
##### Statistics

We compute here basic statistics, _i.e._ mean, median, min, max, standard deviation. You can access it using `statsPE.min` for example.

```julia
statsPB = basic_stat(time_PB)
statsPE = basic_stat(time_PE)
statsLG = basic_stat(time_LG)

factorPE = mean(time_PE) / mean(time_PB)
factorLG = mean(time_LG) / mean(time_PB)

# Write overall timings:
df = DataFrame(
	[getfield(j, i) for i in fieldnames(StatResults), j in [statsPB, statsPE, statsLG]],
	["PlantBiophysics", "plantecophys", "LeafGasExchange"]
)
insertcols!(df, 1, :Stat => [fieldnames(StatResults)...])
CSV.write("benchmark.csv", df)

# Write timing for each sample:
CSV.write("benchmark_full.csv",
	DataFrame(
		"package" => vcat(
			[
				repeat([i.first], length(i.second)) for i in [
					"PlantBiophysics" => time_PB,
					"plantecophys" => time_PE,
					"LeafGasExchange" => time_LG
				]
			]...
		),
		"sample_time" => vcat(time_PB, time_PE, time_LG)
	)
)
```
"""

# ╔═╡ 59ce5f30-3b74-402d-919d-b1b3d1c76788
CSV.read("benchmark.csv", DataFrame)

# ╔═╡ 12497dfb-87ff-4978-801a-c78cd118b032
md"
##### Histogram plotting
"

# ╔═╡ 84a4b3b4-23d8-4b12-9204-110d3ba4724b
md"""
We here display the computational time histogram of each package on the same scale in order to compare them: `PlantBiophysics.jl` (soon), `LeafGasExchange.jl` (a) and `plantecophys`(b). y-axis represents the density (_i.e._ reaching 0.3 means that 30% of the computed times are in this bar). Orange zone represents the interval [mean - standard deviation; mean + standard deviation]. Red dashed line represents the mean. Note that x-axis is logarithmic.

```julia
fig = plot_benchmark_Makie(statsPB, statsPE, statsLG, time_PB, time_PE, time_LG)
    save("benchmark_each_time_steps.png", fig, px_per_unit=3)
```

![](https://github.com/VEZY/PlantBiophysics-paper/blob/main/tutorials/benchmark_each_time_steps.png?raw=true)

!!! note
	This is the plot from the latest commit on <https://github.com/VEZY/PlantBiophysics-paper/>. If you want to make your own benchmarking, run the script that was used to perform it, but careful, it takes a long time to perform!
"""

# ╔═╡ b2958ece-3a4c-497e-88e3-afc45d8f3688
md"""
# References
"""

# ╔═╡ 22ec1978-aaee-4fa1-b739-1246a6aeda20
"""
	StatResults(
		mean::AbstractFloat
    	median::AbstractFloat
    	stddev::AbstractFloat
    	min::AbstractFloat
    	max::AbstractFloat
	)

Structure to hold basic statistics of model performance.
"""
struct StatResults
    mean::AbstractFloat
    median::AbstractFloat
    stddev::AbstractFloat
    min::AbstractFloat
    max::AbstractFloat
end

# ╔═╡ e56f829e-4b83-457a-8e90-a08e071cd5f3
begin
    Base.show(io::IO, m::StatResults) = print(io, m.mean, "(±", m.stddev, ')')

    function Base.show(io::IO, ::MIME"text/plain", m::StatResults)
        print(
            io,
            "Benchmark:",
            "\nMean time -> ", m.mean, " ± ", m.stddev,
            "\nMedian time -> ", m.median,
            "\nMinimum time -> ", m.min,
            "\nMaximum time -> ", m.max,
        )
    end

    md"""
    **Base.show**

    	Base.show(io::IO, m::StatResults)
    	Base.show(io::IO, ::MIME"text/plain", m::StatResults)

    Add a show method for our `StatResults` type.
    """
end

# ╔═╡ 30a52bf1-ff77-4f9b-a87a-e9dd8ad0c1c4
"""
	basic_stat(df)

Compute basic statistics from the benchmarking
"""
function basic_stat(df)
    m = mean(df)
    med = median(df)
    std = Statistics.std(df)
    min = findmin(df)[1]
    max = findmax(df)[1]
    return stat_res(m, med, std, min, max)
end

# ╔═╡ ac330756-f1ac-431c-9e60-7d346e06aa1f
function plot_benchmark_Makie(statsPB, statsPE, statsLG, time_PB, time_PE, time_LG)
    size_inches = (6.7, 5)
    size_pt = 72 .* size_inches
    bins = 220
    noto_sans = assetpath("fonts", "NotoSans-Regular.ttf")
    fig = Figure(backgroundcolor=RGBf(1, 1, 1), resolution=size_pt, font=noto_sans, fontsize=10)
    ep = 1e-9
    extr = extrema(vcat(time_PB, time_PE, time_LG))
    interval = (extr[1] * 1e-1, extr[2])

    ga = fig[1, 1] = GridLayout()

    axa = Axis(ga[1, 1], title="(a) PlantBiophysics.jl", xscale=log10, titlealign=:left, titlesize=10)
    stddevi = poly!(
        axa,
        Rect(max(ep, statsPB.mean - statsPB.stddev), 0.0, 2 * statsPB.stddev, 1),
        color=(:orange, 0.3), yautolimits=false
    )
    moy = vlines!(axa, statsPB.mean; color=:red, linewidth=3, linestyle=:dot)
    hist!(axa, time_PB, normalization=:probability, bins=bins)
    # h = axa.finallimits[].widths[2]
    axislegend(axa, [stddevi, moy], ["95% confidence interval", "Mean"], "", position=:rb,
        orientation=:vertical, labelsize=8, framevisible=false)
    xlims!(axa, interval)

    axb = Axis(ga[2, 1], title="(b) plantecophys", xscale=log10, ylabel="Density", titlealign=:left, titlesize=10)
    stddevi = poly!(
        axb,
        Rect(statsPE.mean - statsPE.stddev, 0.0, 2 * statsPE.stddev, 1),
        color=(:orange, 0.3), yautolimits=false
    )
    vlines!(axb, statsPE.mean; color=:red, linewidth=3, linestyle=:dot)
    hist!(axb, time_PE, normalization=:probability, bins=bins)
    xlims!(axb, interval)

    # axc = Axis(ga[3, 1], title="(c) LeafGasExchange.jl", yminorticks=IntervalsBetween(10),
    #     xscale=log10, xminorticks=IntervalsBetween(10), yminorgridvisible=true, yminorticksvisible=true,
    #     xminorgridvisible=true, xminorticksvisible=true, xlabel="Time (s)")
    axc = Axis(ga[3, 1], title="(c) LeafGasExchange.jl", xscale=log10, xlabel="Time (s)", titlealign=:left, titlesize=10)
    stddevi = poly!(
        axc,
        Rect(statsLG.mean - statsLG.stddev, 0.0, 2 * statsLG.stddev, 1),
        color=(:orange, 0.3), yautolimits=false
    )
    vlines!(axc, statsLG.mean; color=:red, linewidth=3, linestyle=:dot)
    hist!(axc, time_LG, normalization=:probability, bins=bins)
    xlims!(axc, interval)

    rowgap!(ga, 7)
    hidexdecorations!(axa, grid=false)
    hidexdecorations!(axb, grid=false)
    fig
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
PlantBiophysics = "7ae8fcfa-76ad-4ec6-9ea7-5f8f5e2d6ec9"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[compat]
CSV = "~0.10.4"
DataFrames = "~1.3.4"
PlantBiophysics = "~0.3.0"
"""

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

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractFFTs]]
deps = ["ChainRulesCore", "LinearAlgebra"]
git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.1.0"

[[deps.AbstractTrees]]
git-tree-sha1 = "03e0550477d86222521d254b741d470ba17ea0b5"
uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
version = "0.3.4"

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

[[deps.ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "62e51b39331de8911e4a7ff6f5aaf38a5f4cc0ae"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.2.0"

[[deps.ArrayInterface]]
deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "81f0cb60dc994ca17f68d9fb7c942a5ae70d9ee4"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "5.0.8"

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

[[deps.BSON]]
git-tree-sha1 = "306bb5574b0c1c56d7e1207581516c557d105cad"
uuid = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
version = "0.3.5"

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

[[deps.BinaryProvider]]
deps = ["Libdl", "Logging", "SHA"]
git-tree-sha1 = "ecdec412a9abc8db54c0efc5548c64dfce072058"
uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
version = "0.5.10"

[[deps.CSV]]
deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings"]
git-tree-sha1 = "873fb188a4b9d76549b81465b1f75c82aaf59238"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.10.4"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.14.0"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[deps.CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[deps.CommonSubexpressions]]
deps = ["MacroTools", "Test"]
git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
version = "0.3.0"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "b153278a25dd42c65abbf4e62344f9d22e59191b"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.43.0"

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

[[deps.CovarianceEstimation]]
deps = ["LinearAlgebra", "Statistics", "StatsBase"]
git-tree-sha1 = "a3e070133acab996660d31dcf479ea42849e368f"
uuid = "587fd27a-f159-11e8-2dae-1979310e6154"
version = "0.2.7"

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[deps.DataAPI]]
git-tree-sha1 = "fb5f5316dd3fd4c5e7c30a24d50643b73e37cd40"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.10.0"

[[deps.DataDeps]]
deps = ["BinaryProvider", "HTTP", "Libdl", "Reexport", "SHA", "p7zip_jll"]
git-tree-sha1 = "4f0e41ff461d42cfc62ff0de4f1cd44c6e6b3771"
uuid = "124859b0-ceae-595e-8997-d05f6a7a8dfe"
version = "0.7.7"

[[deps.DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "daa21eb85147f72e41f6352a57fccea377e310a9"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.3.4"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

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

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

[[deps.DiffResults]]
deps = ["StaticArrays"]
git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805"
uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
version = "1.0.3"

[[deps.DiffRules]]
deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "28d605d9a0ac17118fe2c5e9ce0fbb76c3ceb120"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.11.0"

[[deps.Distances]]
deps = ["LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "3258d0659f812acde79e8a74b11f17ac06d0ca04"
uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
version = "0.10.7"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Distributions]]
deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "8a6b49396a4058771c5c072239b2e0a76e2e898c"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.58"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

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

[[deps.EzXML]]
deps = ["Printf", "XML2_jll"]
git-tree-sha1 = "0fa3b52a04a4e210aeb1626def9c90df3ae65268"
uuid = "8f5d6c58-4d21-5cfd-889c-e3ad7ee6a615"
version = "1.1.0"

[[deps.FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "9267e5f50b0e12fdfd5a2455534345c4cf2c7f7a"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.14.0"

[[deps.FilePathsBase]]
deps = ["Compat", "Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "129b104185df66e408edd6625d480b7f9e9823a0"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.18"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "246621d23d1f43e3b9c368bf3b72b2331a27c286"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.13.2"

[[deps.FiniteDiff]]
deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"]
git-tree-sha1 = "51c8f36c81badaa0e9ec405dcbabaf345ed18c84"
uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
version = "2.11.1"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.ForwardDiff]]
deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
git-tree-sha1 = "89cc49bf5819f0a10a7a3c38885e7c7ee048de57"
uuid = "f6369f11-7733-5829-9624-2563aa707210"
version = "0.10.29"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[deps.Graphs]]
deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "57c021de207e234108a6f1454003120a1bf350c4"
uuid = "86223c79-3864-5bf0-83f7-82e725a168b6"
version = "1.6.0"

[[deps.HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "0fa77022fe4b511826b39c894c90daf5fce3334a"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.17"

[[deps.IfElse]]
git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.1"

[[deps.Impute]]
deps = ["BSON", "CSV", "DataDeps", "Distances", "IterTools", "LinearAlgebra", "Missings", "NamedDims", "NearestNeighbors", "Random", "Statistics", "StatsBase", "TableOperations", "Tables"]
git-tree-sha1 = "d3fb6342d94030706ad31f05c23514962c29296c"
uuid = "f7bf1975-0170-51b9-8c5f-a992d46b9575"
version = "0.6.8"

[[deps.Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[deps.IniFile]]
git-tree-sha1 = "f550e6e32074c939295eb5ea6de31849ac2c9625"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.1"

[[deps.InlineStrings]]
deps = ["Parsers"]
git-tree-sha1 = "61feba885fac3a407465726d0c330b3055df897f"
uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
version = "1.1.2"

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

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.3"

[[deps.InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.IterTools]]
git-tree-sha1 = "fa6287a4469f5e048d763df38279ee729fbd44e5"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.4.0"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JLD2]]
deps = ["FileIO", "MacroTools", "Mmap", "OrderedCollections", "Pkg", "Printf", "Reexport", "TranscodingStreams", "UUIDs"]
git-tree-sha1 = "81b9477b49402b47fbe7f7ae0b252077f53e4a08"
uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
version = "0.4.22"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

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

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

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

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

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

[[deps.Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

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

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "76c987446e8d555677f064aaac1145c4c17662f8"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.14"

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

[[deps.LsqFit]]
deps = ["Distributions", "ForwardDiff", "LinearAlgebra", "NLSolversBase", "OptimBase", "Random", "StatsBase"]
git-tree-sha1 = "91aa1442e63a77f101aff01dec5a821a17f43922"
uuid = "2fda8390-95c7-5789-9bda-21331edee243"
version = "0.12.1"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

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

[[deps.MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

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

[[deps.MetaGraphsNext]]
deps = ["Graphs", "JLD2"]
git-tree-sha1 = "f8e0351036278130f6bf966cb903e5fddb93778c"
uuid = "fa8bd995-216d-47f1-8a91-f3b68fbeb377"
version = "0.2.0"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

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

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

[[deps.MultiScaleTreeGraph]]
deps = ["AbstractTrees", "DataFrames", "DelimitedFiles", "Graphs", "MetaGraphsNext", "MutableNamedTuples", "OrderedCollections", "Printf", "SHA", "XLSX"]
git-tree-sha1 = "d0dc30c70217c03211d690a6e16f77552aa58fb0"
uuid = "dd4a991b-8a45-4075-bede-262ee62d5583"
version = "0.5.0"

[[deps.MutableNamedTuples]]
git-tree-sha1 = "bd551bdbe872d8f132bb0c5fb6e6dec236fa01f9"
uuid = "af6c499f-54b4-48cc-bbd2-094bba7533c7"
version = "0.1.1"

[[deps.NLSolversBase]]
deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"]
git-tree-sha1 = "50310f934e55e5ca3912fb941dec199b49ca9b68"
uuid = "d41bc354-129a-5804-8e4c-c37616107c6c"
version = "7.8.2"

[[deps.NaNMath]]
git-tree-sha1 = "b086b7ea07f8e38cf122f5016af580881ac914fe"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.7"

[[deps.NamedDims]]
deps = ["AbstractFFTs", "ChainRulesCore", "CovarianceEstimation", "LinearAlgebra", "Pkg", "Requires", "Statistics"]
git-tree-sha1 = "0856b62716585eb90cc1dada226ac9eab5f69aa5"
uuid = "356022a1-0364-5f58-8944-0da4b18d706f"
version = "0.2.47"

[[deps.NearestNeighbors]]
deps = ["Distances", "StaticArrays"]
git-tree-sha1 = "ded92de95031d4a8c61dfb6ba9adb6f1d8016ddd"
uuid = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
version = "0.4.10"

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

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

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.OptimBase]]
deps = ["NLSolversBase", "Printf", "Reexport"]
git-tree-sha1 = "9cb1fee807b599b5f803809e85c81b582d2009d6"
uuid = "87e2bd06-a317-5318-96d9-3ecbac512eee"
version = "2.0.2"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "3114946c67ef9925204cc024a73c9e679cebe0d7"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.8"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "1285416549ccfcdf0c50d4997a94331e88d68413"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.3.1"

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

[[deps.PlantBiophysics]]
deps = ["CSV", "DataFrames", "Dates", "Impute", "LsqFit", "MultiScaleTreeGraph", "MutableNamedTuples", "OrderedCollections", "RecipesBase", "Statistics", "Test", "YAML"]
git-tree-sha1 = "c8891016722498da99581b983f753c4d44f86ab6"
uuid = "7ae8fcfa-76ad-4ec6-9ea7-5f8f5e2d6ec9"
version = "0.3.0"

[[deps.PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "a6062fe4063cdafe78f4a0a81cfffb89721b30e7"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.2"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.3.0"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "dfb54c4e414caa595a1f2ed759b160f5a3ddcba5"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.3.1"

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

[[deps.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

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

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

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

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

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

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

[[deps.SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "6a2f7d70512d205ca8c7ee31bfa9f142fe74310c"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.12"

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

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

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

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

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

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.4"

[[deps.Static]]
deps = ["IfElse"]
git-tree-sha1 = "5309da1cdef03e95b73cd3251ac3a39f887da53e"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.6.4"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "cd56bf18ed715e8b09f06ef8c6b781e6cdc49911"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.4.4"

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

[[deps.StatsAPI]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "c82aaa13b44ea00134f8c9c89819477bd3986ecd"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.3.0"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8977b17906b0a1cc74ab2e3a05faa16cf08a8291"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.16"

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "5950925ff997ed6fb3e985dcce8eb1ba42a0bbe7"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.18"

[[deps.StringEncodings]]
deps = ["Libiconv_jll"]
git-tree-sha1 = "50ccd5ddb00d19392577902f0079267a72c5ab04"
uuid = "69024149-9ee7-55f6-a4c4-859efe599b68"
version = "0.3.5"

[[deps.SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

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

[[deps.TableOperations]]
deps = ["SentinelArrays", "Tables", "Test"]
git-tree-sha1 = "e383c87cf2a1dc41fa30c093b2a19877c83e1bc1"
uuid = "ab02a1b2-a7df-11e8-156e-fb1833f50b87"
version = "1.2.0"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
git-tree-sha1 = "5ce79ce186cc678bbb5c5681ca3379d1ddae11a1"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.7.0"

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

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

[[deps.TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[deps.URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

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

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

[[deps.WeakRefStrings]]
deps = ["DataAPI", "InlineStrings", "Parsers"]
git-tree-sha1 = "b1be2855ed9ed8eac54e5caff2afcdb442d52c23"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "1.4.2"

[[deps.XLSX]]
deps = ["Dates", "EzXML", "Printf", "Tables", "ZipFile"]
git-tree-sha1 = "7fa8618da5c27fdab2ceebdff1da8918c8cd8b5d"
uuid = "fdbf4ff8-1666-58a4-91e7-1b58723a45e0"
version = "0.7.10"

[[deps.XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[deps.YAML]]
deps = ["Base64", "Dates", "Printf", "StringEncodings"]
git-tree-sha1 = "3c6e8b9f5cdaaa21340f841653942e1a6b6561e5"
uuid = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
version = "0.4.7"

[[deps.ZipFile]]
deps = ["Libdl", "Printf", "Zlib_jll"]
git-tree-sha1 = "3593e69e469d2111389a9bd06bac1f3d730ac6de"
uuid = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
version = "0.9.4"

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

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

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

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

# ╔═╡ Cell order:
# ╠═211f1a10-061d-452a-b061-aa5eb5f07428
# ╟─994a4062-39cf-49b7-b19e-63ff90e9a574
# ╟─9c3f432b-4d1f-4765-9989-17e983cee5f7
# ╠═e7e23262-1813-4fe5-80fb-2be20cb54f89
# ╟─65d00bf6-8f9a-40b4-947a-fd46ddc89753
# ╠═acc685f9-7ae8-483d-982e-ff45ccd9e860
# ╟─b8126c39-8c6c-490d-9992-8922e88f8857
# ╟─a0476d0d-64d9-4457-b78d-41519e38e859
# ╟─ae00e44e-d864-4d04-b2e8-fac510ce44bb
# ╟─8419c830-8e35-493b-91b2-adfc97e30a61
# ╟─4dad9b75-3a03-4da1-9aa8-fea02fb48a97
# ╟─60f769b9-bf8c-461c-9739-f31ca43b27b8
# ╟─87cc52dd-9878-422d-a10e-d491db2a04c6
# ╟─8294a065-5285-4221-8553-e38ca23257f2
# ╟─59ce5f30-3b74-402d-919d-b1b3d1c76788
# ╟─12497dfb-87ff-4978-801a-c78cd118b032
# ╟─84a4b3b4-23d8-4b12-9204-110d3ba4724b
# ╟─b2958ece-3a4c-497e-88e3-afc45d8f3688
# ╟─22ec1978-aaee-4fa1-b739-1246a6aeda20
# ╟─e56f829e-4b83-457a-8e90-a08e071cd5f3
# ╟─30a52bf1-ff77-4f9b-a87a-e9dd8ad0c1c4
# ╟─ac330756-f1ac-431c-9e60-7d346e06aa1f
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: VEZY/PlantBiophysics-paper

max_issues_repo_name: VEZY/PlantBiophysics-paper

max_forks_repo_name: VEZY/PlantBiophysics-paper


text: Genie.secret_token!("82310a8b3018ec3407426246a2684948f037e7971c7f6a3d35aefc4277a4e144")


text:
using Plots

func1(R) = R/2
func2(R) = R/2 + [1 ; 0]
func3(R) = R/2 + [0 ; 1]
P = 200
num = 500000
x = []
y = []
for i in 1:num
    point = [rand() ; rand()]
    for j in 1:P
        point = rand([func1(point), func2(point), func3(point)])
    end
    push!(x, point[1])
    push!(y, point[2])
end

scatter(x,y,markersize = 0.000001, legend = false, border=:none, dpi=500)
savefig("C:\\Users\\Yaghoub\\Documents\\GitHub\\A-few-fractals-in-Julia\\Fractals-Fig\\SRT.png")

max_stars_repo_name: shahmari/A-few-fractals-in-Julia

max_issues_repo_name: shahmari/A-few-fractals-in-Julia

max_forks_repo_name: shahmari/A-few-fractals-in-Julia