text:
"""
Description
===========

The `Model` type represents a solution in the genetic algorithm.

Fields
======

- `model::Matrix`: The underlying model in matrix form. Each row is a layer with
  depth to the layer's top in the first column and the layer's resistivity in
  the second.
- `N::Integer`: The number of layers in the model.
- `depth_lims::Array{LayerBC}`: The depth constraints for each layer.
- `res_lims::Array{LayerBC}`: The resistivity constraints for each layer.
- `fitness::Real`: The fitness of the model in the genetic algorithm. A -ve
  fitness indicates the fitness needs to be recalculated.
"""
type Model
    model::Matrix
    N::Integer
    depth_lims::Array{LayerBC}
    res_lims::Array{LayerBC}
    fitness::Real

    function Model(m::Matrix, z::Array{LayerBC}, r::Array{LayerBC})
        if (z[1].min != 0 || z[1].max != 0)
            error("Range of depths for first layer must be (0,0)!")
        end
        if size(m)[1] != length(z) || size(m)[1] != length(r)
            error("Dimensions of zParams or rParams does not match number of layers")
        end

        new(m, size(m)[1], z, r, -1)
    end
end

"""
Description
===========

Create an instance of the `Model` type with a randomly generated model within a
set of constraints.

Arguments
=========

- `N::Integer`: The number of layers to use in the model.
- `depth_lims::Array{LayerBC}`: The depth constraints for each layer.
- `res_lims::Array{LayerBC}`: The resistivity constraints for each layer.

Returns
=======

- `::Model`: An instance of the `Model` type.

Notes
=====

- **Throws an error** if `length(depth_lims) != length(res_lims)`.
- **Throws an error** if the first element of `depth_lims != LayerBC(0,0)`.
"""
function rand_model(N::Integer, depth_lims::Array{LayerBC}, res_lims::Array{LayerBC})
    if depth_lims[1].min != 0 || depth_lims[1].max != 0
        error("Range of depths for first layer must be (0,0)!")
    end

    if length(depth_lims) != N || length(res_lims) != N
        error("Dimensions of depth_lims or res_lims does not match number of layers")
    end

    model = zeros(N,2)
    model[1,2] = rand(res_lims[1].min:res_lims[1].max)
    for n = 2:N
        layer_depth_lim = depth_lims[n]
        layer_res_lim = res_lims[n]
        model[n,1] = rand(layer_depth_lim.min:layer_depth_lim.max)
        model[n,2] = rand(layer_res_lim.min:layer_res_lim.max)
    end

    Model(model, depth_lims, res_lims)
end

"""
Description
===========

Calculate the fitness of a `Model` instance by calculating the RMS misfit
between the model's forward response and some data.

Arguments
=========

- `c::Model`: `Model` Instance to find the fitness of. The instance's `fitness`
  field will be mutated.
- `data::Matrix`: An N×5 matrix containing the data to find the misfit of. The
  column content should be:
    1. Period (s)
    2. Apparent resistivity (Ωm)
    3. Phase (°)
    4. Apparent resistivity uncertainty
    5. Phase uncertainty

Side Effects
============

- Updates the `fitness` field of the instance `c::Model` with the newly
  calculated misfit.

Returns
=======

- `nothing::Void`
"""
function calcfitness!(c::Model, data::Matrix)
    if c.fitness >= 0
        # Fitness is already calculated
        return
    end

    fs = data[:,1].^-1
    ρ, Φ = calc_response(c.model, fs)

    RMSρ = calc_rms(data[:,2], ρ, data[:,4])
    RMSΦ = calc_rms(data[:,3], Φ, data[:,5])

    c.fitness = sqrt(RMSρ^2 + RMSΦ^2)

    # Are the model depths sorted from lowest to highest? Apply a
    # penalty if they aren't
    if !issorted(c.model[:,1])
        c.fitness += 50
    end

    # Do any of the layers in the current model exceed the search
    # boundaries? Apply a penalty proportional to how much they exceed
    # them if they do.
    depth_penalty = 2
    res_penalty = 2
    for N in 1:c.N
        depth_lims = c.depth_lims[N]
        res_lims = c.res_lims[N]
        layer = c.model[N,:]

        if layer[1] < depth_lims.min
            c.fitness += abs(depth_lims.min - layer[1]) * depth_penalty
        elseif layer[1] > depth_lims.max
            c.fitness += abs(layer[1] - depth_lims.max) * depth_penalty
        end
        if layer[2] < res_lims.min
            c.fitness += abs(res_lims.min - layer[2]) * res_penalty
        elseif layer[2] > res_lims.max
            c.fitness += abs(layer[2] - res_lims.max) * res_penalty
        end
    end

    return
end

"""
Description
===========

Crossover two parent `Model` instances to produce two children using SBX.

Arguments
=========

- `a::Model`: The first parent.
- `b::Model`: The second parent.

Returns
=======

- `Tuple{Model,Model}`: The two offsprings of the parents.

Notes
=====

- The offsprings' fitness fields will not be initialised.
"""
function crossover(a::Model, b::Model, η::Integer)
    if η < 1
        error("Distribution parameter (η) must be a non-negative integer")
    end

    ca = deepcopy(a)
    cb = deepcopy(b)

    # Void the RMS calculations for the children
    ca.fitness = -1
    cb.fitness = -1

    for n = 1:a.N
        layer_depth_lims = a.depth_lims[n]
        layer_res_lims = a.res_lims[n]
        # Crossover depths
        if a.model[n,1] < b.model[n,1]
            βa, βb = boundedsbx(a.model[n,1],
                                b.model[n,1],
                                η,
                                layer_depth_lims.min,
                                layer_depth_lims.max)
            ca.model[n,1] = 0.5*(1+βa)*a.model[n,1] + 0.5 * (1-βa) * b.model[n,1]
            cb.model[n,1] = 0.5*(1-βb)*a.model[n,1] + 0.5 * (1+βb) * b.model[n,1]
        end

        # Crossover resistivities
        if a.model[n,2] < b.model[n,2]
            βa, βb = boundedsbx(a.model[n,2],
                                b.model[n,2],
                                η,
                                layer_res_lims.min,
                                layer_res_lims.max)
            ca.model[n,2] = 0.5*(1+βa)*a.model[n,2] + 0.5 * (1-βa) * b.model[n,2]
            cb.model[n,2] = 0.5*(1-βb)*a.model[n,2] + 0.5 * (1+βb) * b.model[n,2]
        end
    end

    # Force first layer to have zero depth
    ca.model[1,1] = 0
    cb.model[1,1] = 0

    ca.model = sortrows(ca.model)
    cb.model = sortrows(cb.model)

    return (ca, cb)
end

"""
Description
===========

Calculate the β values for the bounded SBX of two values.

Arguments
=========

- `p1::Real`: First parent value to crossover.
- `p2::Real`: Second parent value to crossover.
- `η::Integer`: Distribution parameter for the SBX.
- `a::Real`: Lower bound of search space.
- `b::Real`: Upper bound of search space.

Returns
=======

- `Tuple{Real,Real}`: The calculated values of β(u,a) and β(u,b).

Notes
=====

- **Throws an error** if `a > b` or if `p1 > p2`.
- Implementation is based on report by Deb, K. and Jain, H. (#2011017)
"""
function boundedsbx(p1::Real, p2::Real, η::Integer, a::Real, b::Real)
    if a > b
        error("boundedsbx: lower bound a was greater than b")
    end

    if !(p1 < p2)
        error("boundedsbx: Parent 1 was greater than parent 2")
    end
    βa = 1 + (p1 - a) / (p2 - p1)
    βb = 1 + (b - p2) / (p2 - p1)
    γa = 1 / (2 * βa^(η+1))
    γb = 1 / (2 * βb^(η+1))

    βua, βub = 0, 0

    ua = rand()
    if ua <= 0.5/(1-γa)
        βua = (2ua*(1-γa))^(1/(η+1))
    else
        βua = (1 / (2*(1-ua*(1-γa))))^(1/(η+1))
    end

    ub = rand()
    if ub <= 0.5/(1-γb)
        βub = (2ub*(1-γb))^(1/(η+1))
    else
        βub = (1 / (2*(1-ub*(1-γb))))^(1/(η+1))
    end

    return (βua, βub)
end

"""
Description
===========

Mutate a `Model` instance. For each depth and resistivity value in the model,
pick a random number. If it is less than the probability of mutation, pick a new
value within the constraints of the model.

Arguments
=========

- `c::Model`: Model to mutate.
- `Pm::Real`: Probability of mutation. Should be in the range [0,1].

Side Effects
============

- If a mutation occurs, `c`'s `model` field will be altered.

Returns
=======

- `nothing::Void`
"""
function mutate!(c::Model, Pm::Real)
    for n in 1:c.N
        # Mutate depth
        if rand() < Pm
            # Do nothing to the depth if this is the top layer
            if n != 1
                layer_depth_lims = c.depth_lims[n]
                c.model[n,1] = rand(layer_depth_lims.min:layer_depth_lims.max)
                sortrows(c.model)
                c.fitness = -1
            end
        end

        # Mutate resistivity
        if rand() < Pm
            layer_res_params = c.res_lims[n]
            c.model[n,2] = rand(layer_res_params.min:layer_res_params.max)
            sortrows(c.model)
            c.fitness = -1
        end
    end

    return
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## dummy stub for https://github.com/JuliaBinaryWrappers/LibSSH2_jll.jl

module LibSSH2_jll

using Libdl
using MbedTLS_jll

const PATH_list = String[]
const LIBPATH_list = String[]

export libssh2

# These get calculated in __init__()
PATH = Ref("")
LIBPATH = Ref("")
artifact_dir = ""
libssh2_handle = C_NULL
libssh2_path = ""

if Sys.iswindows()
    const libssh2 = "libssh2.dll"
elseif Sys.isapple()
    const libssh2 = "@rpath/libssh2.1.dylib"
else
    const libssh2 = "libssh2.so.1"
end

function __init__()
    global artifact_dir = dirname(Sys.BINDIR)
    global LIBPATH[] = joinpath(Sys.BINDIR, Base.LIBDIR, "julia")
    global libssh2_handle = dlopen(libssh2)
    global libssh2_path = dlpath(libssh2_handle)
end

is_available() = true

end  # module LibSSH2_jll

ext: jl

max_stars_repo_name: saberdy/julia

max_issues_repo_name: saberdy/julia

max_forks_repo_name: saberdy/julia


text:
module Atmos

module Turbulence

""" Smagorinsky Coefficient [dimensionless] """
function C_smag end

""" Drag coefficient """
function C_drag end

end # module Turbulence

end # module Atmos

max_stars_repo_name: trontrytel/CLIMA

max_issues_repo_name: trontrytel/CLIMA

max_forks_repo_name: trontrytel/CLIMA


text:
# generate all possible combinations
const VALS = Int8.(1:10)

function gen_all_solved_comb()
    all_comb = Iterators.product([VALS for i in 1:4]...) |> collect ;
    move_combine_compact!.(all_comb)
end

const MOVE_LOOKUP = gen_all_comb()

function move_lookup!(x)
    x .= MOVE_LOOKUP[x[1], x[2], x[3], x[4]]
    x
end


text:
function extensiblefunction!(@nospecialize(f::Function), varargs...)
    return extensiblefunction!(f, varargs)
end

function extensiblefunction!(@nospecialize(f::Function),
                             @nospecialize(varargs::Tuple))
    global _registry_extensibleunion_to_genericfunctions
    global _registry_genericfunctions_to_extensibleunions
    if !haskey(_registry_genericfunctions_to_extensibleunions, f)
        _registry_genericfunctions_to_extensibleunions[f] = Set{Any}()
    end
    for i = 1:length(varargs)
        if isextensibleunion(varargs[i])
            push!(_registry_extensibleunion_to_genericfunctions[varargs[i]],
                  f)
            push!(_registry_genericfunctions_to_extensibleunions[f],
                  varargs[i])
        else
            throw(ArgumentError(
                "Argument is not a registered extensible union."))
        end
    end
    _update_all_methods_for_extensiblefunction!(f)
    return f
end

function isextensiblefunction(@nospecialize(f::Function))
    global _registry_genericfunctions_to_extensibleunions
    return haskey(_registry_genericfunctions_to_extensibleunions, f)
end

macro extensible(fdef)
    @assert fdef.head == :(=) || fdef.head == :(:function)
    ext_unions = Symbol[]
    name  = fdef.args[1].args[1]
    fargs = fdef.args[1].args[2:end]
    for (i, arg) in enumerate(fargs)
        if arg isa Expr && arg.head == :(::)
            T = arg.args[end]
            if T isa Expr && T.head == :braces && length(T.args) == 1
                push!(ext_unions, T.args[1])
                fdef.args[1].args[1+i].args[end] = T.args[1]
            end
        end
    end
    fdef
    quote
        $(fdef)
        extensiblefunction!($name, $(ext_unions...))
    end |> esc
end


text:
function polynomial_features(data_in::Array; dim::Int=2,
                             subset::Vector=collect(1:size(data_in, dim)), degree::Int=2,
                             include_bias::Bool=true, interaction_only::Bool=false)

    # Data type of the elements of the input array
    T = eltype(data_in)

    # If features are in columns and examples are in rows
    if dim == 2

        # Size of original input matrix
        num_obs, num_feats_in = size(data_in)

        # Create matrix of the chosen subset of features that will generate the new polynomial features
        data_sub = data_in[:, subset]

        # How many original features are chosen to generate the polynomial features
        num_feats_sub = length(subset)

        # If we're only generating the interaction terms (i.e. only single power cross terms)
        if interaction_only

            # Number of columns needed for the new features
            num_feats_poly = 0
            for deg = 2:degree
                num_feats_poly += binomial(num_feats_sub, deg)
            end

            # Initialise design matrix with:
            # column for bias if wanted i.e. if include_bias = true
            # columns for original features i.e. data_in
            # columns for the polynomial features that will be generated (initialised to zero)
            if include_bias
                data_out = [ones(T, num_obs) data_in zeros(T, num_obs, num_feats_poly)]
            else
                data_out = [data_in zeros(T, num_obs, num_feats_poly)]
            end

            # Generate polynomial features and add to data_out
            ind = (include_bias ? 1 : 0) + num_feats_in + 1        # column number to start placing the new features
            for deg = 2:degree                                     # generate features for each power - 2:degree
                for c in combinations(collect(1:num_feats_sub), deg)      # create combinations iterator
                    p = zeros(Int, num_feats_sub)                  # zeros vector of length = number of features in the subset
                    p[c] = 1                                       # change this combination's elements to 1
                    data_out[:, ind] = prod(data_sub .^ p', 2)     # use the power vector to create the new feature
                    ind += 1                                       # move one column to the right to place next new feature
                end
            end

        # If we're generating all the terms
        else

            # Number of columns needed for the new features
            num_feats_poly = 0
            for deg = 2:degree
                num_feats_poly += binomial(num_feats_sub + deg - 1, num_feats_sub - 1)
            end

            # Initialise design matrix with:
            # column for bias if wanted i.e. if include_bias = true
            # columns for original features i.e. data_in
            # columns for the polynomial features that will be generated (initialised to zero)
            if include_bias
                data_out = [ones(T, num_obs) data_in zeros(T, num_obs, num_feats_poly)]
            else
                data_out = [data_in zeros(T, num_obs, num_feats_poly)]
            end

            # Generate polynomial features and add to data_out
            ind = (include_bias ? 1 : 0) + num_feats_in + 1                                 # column number to start placing the new features
            for deg = 2:degree                                                              # generate features for each power - 2:degree
                tmp = [(num_feats_sub + deg); zeros(Int, num_feats_sub)]                    # temp vector used to convert format of combination
                for c in combinations(collect(1:(num_feats_sub + deg - 1)), num_feats_sub - 1)     # create combinations iterator
                    tmp[num_feats_sub:-1:2] = c                                             # place combination into tmp vector
                    p = -diff(tmp)-1                                                        # vector of powers
                    data_out[:, ind] = prod(data_sub .^ p', 2)                              # use the power vector to create the new feature
                    ind += 1                                                                # move one column to the right to place next new feature
                end
            end

        end

    # If features are in rows and examples are in cols
    elseif dim == 1

        # Size of original input matrix
        num_feats_in, num_obs = size(data_in)

        # Create matrix of the chosen subset of features that will generate the new polynomial features
        data_sub = data_in[subset, :]

        # How many original features are chosen to generate the polynomial features
        num_feats_sub = length(subset)

        # If we're only generating the interaction terms (i.e. only single power cross terms)
        if interaction_only

            # Number of rows needed for the new features
            num_feats_poly = 0
            for deg = 2:degree
                num_feats_poly += binomial(num_feats_sub, deg)
            end

            # Initialise design matrix with:
            # row for bias if wanted i.e. if include_bias = true
            # rows for original features i.e. data_in
            # rows for the polynomial features that will be generated (initialised to zero)
            if include_bias
                data_out = [ones(T, 1, num_obs) ; data_in ; zeros(T, num_feats_poly, num_obs)]
            else
                data_out = [data_in ; zeros(T, num_feats_poly, num_obs)]
            end

            # Generate polynomial features and add to data_out
            ind = (include_bias ? 1 : 0) + num_feats_in + 1       # row number to start placing the new features
            for deg = 2:degree                                    # generate features for each power - 2:degree
                for c in combinations(collect(1:num_feats_sub), deg)     # create combinations iterator
                    p = zeros(Int, num_feats_sub)                 # zeros vector of length = number of features in the subset
                    p[c] = 1                                      # change this combination's elements to 1
                    data_out[ind, :] = prod(data_sub .^ p, 1)     # use the power vector to create the new feature
                    ind += 1                                      # move one column to the right to place next new feature
                end
            end

        # If we're generating all the terms
        else

            # Number of columns needed for the new features
            num_feats_poly = 0
            for deg = 2:degree
                num_feats_poly += binomial(num_feats_sub + deg - 1, num_feats_sub - 1)
            end

            # Initialise design matrix with:
            # row for bias if wanted i.e. if include_bias = true
            # rows for original features i.e. data_in
            # rows for the polynomial features that will be generated (initialised to zero)
            if include_bias
                data_out = [ones(T, 1, num_obs) ; data_in ; zeros(T, num_feats_poly, num_obs)]
            else
                data_out = [data_in ; zeros(T, num_feats_poly, num_obs)]
            end

            # Generate polynomial features and add to data_out
            ind = (include_bias ? 1 : 0) + num_feats_in + 1                                 # row number to start placing the new features
            for deg = 2:degree                                                              # generate features for each power - 2:degree
                tmp = [(num_feats_sub + deg); zeros(Int, num_feats_sub)]                    # temp vector used to convert format of combination
                for c in combinations(collect(1:(num_feats_sub + deg - 1)), num_feats_sub - 1)     # create combinations iterator
                    tmp[num_feats_sub:-1:2] = c                                             # place combination into tmp vector
                    p = -diff(tmp)-1                                                        # vector of powers
                    data_out[ind, :] = prod(data_sub .^ p, 1)                               # use the power vector to create the new feature
                    ind += 1                                                                # move one column to the right to place next new feature
                end
            end

        end

    # If dim entered is not in range i.e. must be 1 or 2
    else throw(DimensionMismatch("Dimension must be 1 (features are rows) or 2, the default, (features are cols)."))

    end

    return data_out
end


text:
###############################################################################
#
# Diff function to minimize by ConstrainedRootSolvers
# Useful for Ballberry, Leuning, and Medlyn models
#
###############################################################################
"""
    solution_diff!(x::FT,
                canopyi::CanopyLayer{FT},
                hs::LeafHydraulics{FT},
                svc::AbstractSoilVC{FT},
                psoil::FT,
                swc::FT,
                envir::AirLayer{FT},
                sm::OptimizationStomatalModel{FT},
                bt::AbstractBetaFunction{FT},
                mode::GlcDrive,
                ind::Int
    ) where {FT<:AbstractFloat}
    solution_diff!(x::FT,
                canopyi::CanopyLayer{FT},
                hs::LeafHydraulics{FT},
                envir::AirLayer{FT},
                sm::AbstractStomatalModel{FT},
                mode::AbstractDrive,
                ind::Int
    ) where {FT<:AbstractFloat}

Calculate the difference to be minimized for a given
- `x` Assumed leaf diffusive conductance or stomatal conductance, depending on
    `mode`
- `canopyi`[`CanopyLayer`](@ref) type struct
- `hs` Leaf hydraulic system
- `psoil` Soil water potential `[MPa]`
- `swc` Soil water content
- `envir`[`AirLayer`] type struct
- `sm` [`EmpiricalStomatalModel`](@ref) or [`OptimizationStomatalModel`](@ref)
- `bt` [`AbstractBetaFunction`](@ref) type struct
- `mode` [`GlcDrive`](@ref) or [`GswDrive`](@ref) mode
- `ind` Nth leaf in the canopy layer

The former function works for all empirical stomatal models, and the latter
    works for all optimization based models.
"""
function solution_diff!(
            x::FT,
            canopyi::CanopyLayer{FT},
            hs::LeafHydraulics{FT},
            svc::AbstractSoilVC{FT},
            psoil::FT,
            swc::FT,
            envir::AirLayer{FT},
            sm::EmpiricalStomatalModel{FT},
            bt::AbstractBetaG{FT},
            mode::GlcDrive,
            ind::Int
) where {FT<:AbstractFloat}
    # unpack variables
    @unpack ps = canopyi;
    g_bc  = canopyi.g_bc[ind];
    g_m   = canopyi.g_m[ind];

    # update photosynthesis for ps
    leaf_photosynthesis!(ps, envir, GCO₂Mode(), x);

    # calculate g_sw from stomatal model
    β    = β_factor(hs, svc, bt, hs.p_element[end], psoil, swc);
    g_md = stomatal_conductance(sm, ps, envir, β);
    g_md = min(canopyi.g_max, g_md);

    # calculate model predicted g_lc
    g_lm = 1 / (FT(1.6)/g_md + 1/g_bc + 1/g_m);

    return g_lm - x
end

function solution_diff!(
            x::FT,
            canopyi::CanopyLayer{FT},
            hs::LeafHydraulics{FT},
            svc::AbstractSoilVC{FT},
            psoil::FT,
            swc::FT,
            envir::AirLayer{FT},
            sm::EmpiricalStomatalModel{FT},
            bt::AbstractBetaV{FT},
            mode::GlcDrive,
            ind::Int
) where {FT<:AbstractFloat}
    # unpack variables
    ps   = canopyi.ps;
    g_bc = canopyi.g_bc[ind];
    g_m  = canopyi.g_m[ind];

    # update photosynthesis for ps
    leaf_photosynthesis!(ps, envir, GCO₂Mode(), x);

    # make beta correction over the photosynthesis system
    β    = β_factor(hs, svc, bt, hs.p_element[end], psoil, swc);
    _rat = ps.PSM.v_cmax25_ww * β / ps.PSM.v_cmax25;

    # make sure it does not conflict with update_leaf_TP!
    if typeof(ps.PSM) <: C4VJPModel
        ps.PSM.r_d25    *= _rat;
        ps.PSM.r_d      *= _rat;
        ps.PSM.v_cmax25 *= _rat;
        ps.PSM.v_cmax   *= _rat;
        ps.PSM.v_pmax25 *= _rat;
        ps.PSM.v_pmax   *= _rat;
    elseif typeof(ps.PSM) <: C3VJPModel
        ps.PSM.j_max25  *= _rat;
        ps.PSM.j_max    *= _rat;
        ps.PSM.r_d25    *= _rat;
        ps.PSM.r_d      *= _rat;
        ps.PSM.v_cmax25 *= _rat;
        ps.PSM.v_cmax   *= _rat;
    else
        ps.PSM.r_d25    *= _rat;
        ps.PSM.r_d      *= _rat;
        ps.PSM.v_cmax25 *= _rat;
        ps.PSM.v_cmax   *= _rat;
        ps.PSM.v_qmax   *= _rat;
    end;

    # calculate g_sw from stomatal model
    g_md = stomatal_conductance(sm, ps, envir, FT(1));
    g_md = min(canopyi.g_max, g_md);

    # calculate model predicted g_lc
    g_lm = 1 / (FT(1.6)/g_md + 1/g_bc + 1/g_m);

    return g_lm - x
end

function solution_diff!(
            x::FT,
            canopyi::CanopyLayer{FT},
            hs::LeafHydraulics{FT},
            envir::AirLayer{FT},
            sm::OSMEller,
            mode::GlcDrive,
            ind::Int
) where {FT<:AbstractFloat}
    # unpack variables
    @unpack g_max, g_min, p_sat, ps = canopyi;
    @unpack P_AIR, p_H₂O = envir;
    g_bc = canopyi.g_bc[ind];
    g_bw = canopyi.g_bw[ind];
    g_m  = canopyi.g_m[ind];

    # calculate the limit of g_lc to ensure x in the physiological range
    _glh = 1 / (1/g_bc + FT(1.6)/g_max + 1/g_m);
    _gll = 1 / (1/g_bc + FT(1.6)/g_min + 1/g_m);

    if x > _glh
        x = _glh + FT(1e-4);
        canopyi.g_sw[ind] = 2 * g_max;

        return FT(0)
    elseif x< _gll
        x = _gll - FT(1e-4);
        canopyi.g_sw[ind] = FT(0);

        return FT(0)
    else
        # update photosynthesis from x-FT(1e-3)
        g_lc = x - FT(1e-3);
        leaf_photosynthesis!(ps, envir, GCO₂Mode(), g_lc);

        g_sc = 1 / ( 1/g_lc - 1/g_m - 1/g_bc );
        g_sw = g_sc * FT(1.6);
        g_lw = 1 / (1/g_sw + 1/g_bw);
        a_1  = canopyi.ps.PSM.a_net;
        e_1  = g_lw * (p_sat - p_H₂O) / P_AIR;
        k_1  = xylem_risk(hs, e_1);

        # update photosynthesis from x
        g_lc = x;
        leaf_photosynthesis!(ps, envir, GCO₂Mode(), g_lc);

        g_sc = 1 / ( 1/g_lc - 1/g_m - 1/g_bc );
        g_sw = g_sc * FT(1.6);
        g_lw = 1 / (1/g_sw + 1/g_bw);
        a_2  = canopyi.ps.PSM.a_net;
        e_2  = g_lw * (p_sat - p_H₂O) / P_AIR;
        k_2  = xylem_risk(hs, e_2);

        ∂A∂E = (a_2 - a_1) / (e_2 - e_1);
        ∂Θ∂E = (k_1 - k_2) / (e_2 - e_1) * a_2 / k_2;
        diff = ∂A∂E - ∂Θ∂E;

        #= used for debugging
        @show x;
        @show ∂A∂E;
        @show ∂Θ∂E;
        println();
        #sleep(0.1);
        # =#

        return diff
    end
end

function solution_diff!(
            x::FT,
            canopyi::CanopyLayer{FT},
            hs::LeafHydraulics{FT},
            envir::AirLayer{FT},
            sm::OSMSperry{FT},
            mode::GlcDrive,
            ind::Int
) where {FT<:AbstractFloat}
    # unpack variables
    @unpack g_max, g_min, kr_max, p_sat, ps = canopyi;
    @unpack P_AIR, p_H₂O = envir;
    a_max = canopyi.a_max[ind];
    g_bc  = canopyi.g_bc[ind];
    g_bw  = canopyi.g_bw[ind];
    g_m   = canopyi.g_m[ind];

    # calculate the limit of g_lc to ensure x in the physiological range
    _glh = 1 / (1/g_bc + FT(1.6)/g_max + 1/g_m);
    _gll = 1 / (1/g_bc + FT(1.6)/g_min + 1/g_m);

    if x > _glh
        x = _glh + FT(1e-4);
        canopyi.g_sw[ind] = 2 * g_max;

        return FT(0)
    elseif x< _gll
        x = _gll - FT(1e-4);
        canopyi.g_sw[ind] = FT(0);

        return FT(0)
    else
        # update photosynthesis from x-FT(1e-3)
        g_lc = x - FT(1e-3);
        leaf_photosynthesis!(ps, envir, GCO₂Mode(), g_lc);

        g_sc = 1 / ( 1/g_lc - 1/g_m - 1/g_bc );
        g_sw = g_sc * FT(1.6);
        g_lw = 1 / (1/g_sw + 1/g_bw);
        a_1  = canopyi.ps.PSM.a_net;
        e_1  = g_lw * (p_sat - p_H₂O) / P_AIR;
        k_1  = xylem_risk(hs, e_1);

        # update photosynthesis from x
        g_lc = x;
        leaf_photosynthesis!(ps, envir, GCO₂Mode(), g_lc);

        g_sc = 1 / ( 1/g_lc - 1/g_m - 1/g_bc );
        g_sw = g_sc * FT(1.6);
        g_lw = 1 / (1/g_sw + 1/g_bw);
        a_2  = canopyi.ps.PSM.a_net;
        e_2  = g_lw * (p_sat - p_H₂O) / P_AIR;
        k_2  = xylem_risk(hs, e_2);

        ∂A∂E = (a_2 - a_1) / (e_2 - e_1);
        ∂Θ∂E = (k_1 - k_2) / (e_2 - e_1) * a_max / kr_max;
        diff = ∂A∂E - ∂Θ∂E;

        #= used for debugging
        @show x;
        @show ∂A∂E;
        @show ∂Θ∂E;
        println();
        #sleep(0.1);
        # =#

        return diff
    end
end

function solution_diff!(
            x::FT,
            canopyi::CanopyLayer{FT},
            hs::LeafHydraulics{FT},
            envir::AirLayer{FT},
            sm::OSMWang{FT},
            mode::GlcDrive,
            ind::Int
) where {FT<:AbstractFloat}
    # unpack variables
    @unpack ec, g_max, g_min, p_sat, ps = canopyi;
    @unpack P_AIR, p_H₂O = envir;
    g_bc  = canopyi.g_bc[ind];
    g_bw  = canopyi.g_bw[ind];
    g_m   = canopyi.g_m[ind];

    # calculate the limit of g_lc to ensure x in the physiological range
    _glh = 1 / (1/g_bc + FT(1.6)/g_max + 1/g_m);
    _gll = 1 / (1/g_bc + FT(1.6)/g_min + 1/g_m);

    if x > _glh
        x = _glh + FT(1e-4);
        canopyi.g_sw[ind] = 2 * g_max;

        return FT(0)
    elseif x< _gll
        x = _gll - FT(1e-4);
        canopyi.g_sw[ind] = FT(0);

        return FT(0)
    else
        # update photosynthesis from x-FT(1e-3)
        g_lc = x - FT(1e-3);
        leaf_photosynthesis!(ps, envir, GCO₂Mode(), g_lc);

        g_sc = 1 / ( 1/g_lc - 1/g_m - 1/g_bc );
        g_sw = g_sc * FT(1.6);
        g_lw = 1 / (1/g_sw + 1/g_bw);
        a_1  = canopyi.ps.PSM.a_net;
        e_1  = g_lw * (p_sat - p_H₂O) / P_AIR;

        # update photosynthesis from x
        g_lc = x;
        leaf_photosynthesis!(ps, envir, GCO₂Mode(), g_lc);

        g_sc = 1 / ( 1/g_lc - 1/g_m - 1/g_bc );
        g_sw = g_sc * FT(1.6);
        g_lw = 1 / (1/g_sw + 1/g_bw);
        a_2  = canopyi.ps.PSM.a_net;
        e_2  = g_lw * (p_sat - p_H₂O) / P_AIR;

        ∂A∂E = (a_2 - a_1) / (e_2 - e_1);
        ∂Θ∂E = a_2 / max(ec - e_2, FT(1e-7));
        diff = ∂A∂E - ∂Θ∂E;

        #= used for debugging
        @show x;
        @show ∂A∂E;
        @show ∂Θ∂E;
        println();
        #sleep(0.1);
        # =#

        return diff
    end
end

function solution_diff!(
            x::FT,
            canopyi::CanopyLayer{FT},
            hs::LeafHydraulics{FT},
            envir::AirLayer{FT},
            sm::OSMWAP{FT},
            mode::GlcDrive,
            ind::Int
) where {FT<:AbstractFloat}
    # unpack variables
    @unpack g_max, g_min, p_sat, ps = canopyi;
    @unpack P_AIR, p_H₂O = envir;
    g_bc = canopyi.g_bc[ind];
    g_bw = canopyi.g_bw[ind];
    g_m  = canopyi.g_m[ind];

    # calculate the limit of g_lc to ensure x in the physiological range
    _glh = 1 / (1/g_bc + FT(1.6)/g_max + 1/g_m);
    _gll = 1 / (1/g_bc + FT(1.6)/g_min + 1/g_m);

    if x > _glh
        x = _glh + FT(1e-4);
        canopyi.g_sw[ind] = 2 * g_max;

        return FT(0)
    elseif x< _gll
        x = _gll - FT(1e-4);
        canopyi.g_sw[ind] = FT(0);

        return FT(0)
    else
        # update photosynthesis from x-FT(1e-3)
        g_lc = x - FT(1e-3);
        leaf_photosynthesis!(ps, envir, GCO₂Mode(), g_lc);

        g_sc = 1 / ( 1/g_lc - 1/g_m - 1/g_bc );
        g_sw = g_sc * FT(1.6);
        g_lw = 1 / (1/g_sw + 1/g_bw);
        a_1  = canopyi.ps.PSM.a_net;
        e_1  = g_lw * (p_sat - p_H₂O) / P_AIR;
        p_1  = end_pressure(hs, e_1);

        # update photosynthesis from x
        g_lc = x;
        leaf_photosynthesis!(ps, envir, GCO₂Mode(), g_lc);

        g_sc = 1 / ( 1/g_lc - 1/g_m - 1/g_bc );
        g_sw = g_sc * FT(1.6);
        g_lw = 1 / (1/g_sw + 1/g_bw);
        a_2  = canopyi.ps.PSM.a_net;
        e_2  = g_lw * (p_sat - p_H₂O) / P_AIR;
        p_2  = end_pressure(hs, e_2);

        ∂A∂E = (a_2 - a_1) / (e_2 - e_1);
        ∂Θ∂E = (-2 * sm.a * p_2 + sm.b) / (e_2 - e_1) * (p_1 - p_2);
        diff = ∂A∂E - ∂Θ∂E;

        #= used for debugging
        @show x;
        @show ∂A∂E;
        @show ∂Θ∂E;
        println();
        #sleep(0.1);
        # =#

        return diff
    end
end

function solution_diff!(
            x::FT,
            canopyi::CanopyLayer{FT},
            hs::LeafHydraulics{FT},
            envir::AirLayer{FT},
            sm::OSMWAPMod{FT},
            mode::GlcDrive,
            ind::Int
) where {FT<:AbstractFloat}
    # unpack variables
    @unpack g_max, g_min, p_sat, ps = canopyi;
    @unpack P_AIR, p_H₂O = envir;
    g_bc = canopyi.g_bc[ind];
    g_bw = canopyi.g_bw[ind];
    g_m  = canopyi.g_m[ind];

    # calculate the limit of g_lc to ensure x in the physiological range
    _glh = 1 / (1/g_bc + FT(1.6)/g_max + 1/g_m);
    _gll = 1 / (1/g_bc + FT(1.6)/g_min + 1/g_m);

    if x > _glh
        x = _glh + FT(1e-4);
        canopyi.g_sw[ind] = 2 * g_max;

        return FT(0)
    elseif x< _gll
        x = _gll - FT(1e-4);
        canopyi.g_sw[ind] = FT(0);

        return FT(0)
    else
        # update photosynthesis from x-FT(1e-3)
        g_lc = x - FT(1e-3);
        leaf_photosynthesis!(ps, envir, GCO₂Mode(), g_lc);

        g_sc = 1 / ( 1/g_lc - 1/g_m - 1/g_bc );
        g_sw = g_sc * FT(1.6);
        g_lw = 1 / (1/g_sw + 1/g_bw);
        a_1  = canopyi.ps.PSM.a_net;
        e_1  = g_lw * (p_sat - p_H₂O) / P_AIR;
        p_1  = end_pressure(hs, e_1);

        # update photosynthesis from x
        g_lc = x;
        leaf_photosynthesis!(ps, envir, GCO₂Mode(), g_lc);

        g_sc = 1 / ( 1/g_lc - 1/g_m - 1/g_bc );
        g_sw = g_sc * FT(1.6);
        g_lw = 1 / (1/g_sw + 1/g_bw);
        a_2  = canopyi.ps.PSM.a_net;
        e_2  = g_lw * (p_sat - p_H₂O) / P_AIR;
        p_2  = end_pressure(hs, e_2);

        ∂A∂E = (a_2 - a_1) / (e_2 - e_1);
        ∂Θ∂E = (-2 * sm.a * p_2 * a_2) / (e_2 - e_1) * (p_1 - p_2);
        diff = ∂A∂E - ∂Θ∂E;

        #= used for debugging
        @show x;
        @show ∂A∂E;
        @show ∂Θ∂E;
        println();
        #sleep(0.1);
        # =#

        return diff
    end
end

function solution_diff!(
            x::FT,
            canopyi::CanopyLayer{FT},
            hs::LeafHydraulics{FT},
            envir::AirLayer{FT},
            sm::OptimizationStomatalModel{FT},
            mode::GswDrive,
            ind::Int
) where {FT<:AbstractFloat}
    # gsw has already been guaranteed to be in the range
    g_bc = canopyi.g_bc[ind];
    g_bw = canopyi.g_bw[ind];
    g_m  = canopyi.g_m[ind];
    g_lc = 1 / (1/g_bc + FT(1.6)/x + 1/g_m);

    return solution_diff!(g_lc, canopyi, hs, envir, sm, GlcDrive(), ind)
end


text:
function Base.:*(H::Hᶜ{ITensor}, v::ITensor)
  ∑h = H.∑h
  Hᴸ = H.Hᴸ
  Hᴿ = H.Hᴿ
  ψ = H.ψ
  ψ′ = dag(ψ)'
  n = H.n
  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)
  δˡ = δ(l[n], l′[n])
  δˡ⁻¹ = δ(l[n - 1], l′[n - 1])
  δʳ = δ(r[n], r′[n])
  δʳ⁺¹ = δ(r[n + 1], r′[n + 1])
  Hᶜᴸv = v * Hᴸ[n] * dag(δʳ)
  Hᶜᴿv = v * δˡ * Hᴿ[n]
  Hᶜʰv =
    v * ψ.AL[n] * δˡ⁻¹ * ψ′.AL[n] * ∑h[(n, n + 1)] * ψ.AR[n + 1] * dag(δʳ⁺¹) * ψ′.AR[n + 1]
  Hᶜv = Hᶜᴸv + Hᶜʰv + Hᶜᴿv
  return Hᶜv * dag(δˡ) * δʳ
end

function Base.:*(H::Hᴬᶜ{ITensor}, v::ITensor)
  ∑h = H.∑h
  Hᴸ = H.Hᴸ
  Hᴿ = H.Hᴿ
  ψ = H.ψ
  ψ′ = dag(ψ)'
  Nsites = nsites(ψ)
  range_∑h = nrange(ψ, ∑h[1])
  n = H.n
  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)
  s = siteinds(only, ψ)
  s′ = siteinds(only, ψ′)

  δˢ(n) = δ(s[n], s′[n])
  δˡ(n) = δ(l[n], l′[n])
  δʳ(n) = δ(r[n], r′[n])

  Hᴬᶜᴸv = v * Hᴸ[n - 1] * dag(δˢ(n)) * dag(δʳ(n))
  Hᴬᶜᴿv = v * δˡ(n - 1) * dag(δˢ(n)) * Hᴿ[n]
  Hᴬᶜʰ¹v = v * ψ.AL[n - 1] * δˡ(n - 2) * ψ′.AL[n - 1] * ∑h[(n - 1, n)] * dag(δʳ(n))
  Hᴬᶜʰ²v = v * δˡ(n - 1) * ψ.AR[n + 1] * dag(δʳ(n + 1)) * ψ′.AR[n + 1] * ∑h[(n, n + 1)]
  Hᴬᶜv = Hᴬᶜᴸv + Hᴬᶜʰ¹v + Hᴬᶜʰ²v + Hᴬᶜᴿv
  return Hᴬᶜv * dag(δˡ(n - 1)) * δˢ(n) * δʳ(n)
end

function tdvp_iteration_sequential(
  solver::Function,
  ∑h::InfiniteSum{ITensor},
  ψ::InfiniteCanonicalMPS;
  (ϵᴸ!)=fill(1e-15, nsites(ψ)),
  (ϵᴿ!)=fill(1e-15, nsites(ψ)),
  time_step,
  solver_tol=(x -> x / 100),
  eager=true,
)
  Nsites = nsites(ψ)
  ϵᵖʳᵉˢ = max(maximum(ϵᴸ!), maximum(ϵᴿ!))
  _solver_tol = solver_tol(ϵᵖʳᵉˢ)
  ψᴴ = dag(ψ)
  ψ′ = ψᴴ'
  # XXX: make this prime the center sites
  ψ̃ = prime(linkinds, ψᴴ)

  # TODO: replace with linkinds(ψ)
  l = CelledVector([commoninds(ψ.AL[n], ψ.AL[n + 1]) for n in 1:Nsites])
  l′ = CelledVector([commoninds(ψ′.AL[n], ψ′.AL[n + 1]) for n in 1:Nsites])
  r = CelledVector([commoninds(ψ.AR[n], ψ.AR[n + 1]) for n in 1:Nsites])
  r′ = CelledVector([commoninds(ψ′.AR[n], ψ′.AR[n + 1]) for n in 1:Nsites])

  ψ = copy(ψ)
  C̃ = InfiniteMPS(Vector{ITensor}(undef, Nsites), translator(ψ))
  Ãᶜ = InfiniteMPS(Vector{ITensor}(undef, Nsites), translator(ψ))
  Ãᴸ = InfiniteMPS(Vector{ITensor}(undef, Nsites), translator(ψ))
  Ãᴿ = InfiniteMPS(Vector{ITensor}(undef, Nsites), translator(ψ))
  eᴸ = Vector{Float64}(undef, Nsites)
  eᴿ = Vector{Float64}(undef, Nsites)
  for n in 1:Nsites
    hᴸ = InfiniteMPS([
      δ(only(l[k - 2]), only(l′[k - 2])) *
      ψ.AL[k - 1] *
      ψ.AL[k] *
      ∑h[(k - 1, k)] *
      ψ′.AL[k - 1] *
      ψ′.AL[k] for k in 1:Nsites
    ])
    hᴿ = InfiniteMPS([
      δ(only(dag(r[k + 2])), only(dag(r′[k + 2]))) *
      ψ.AR[k + 2] *
      ψ.AR[k + 1] *
      ∑h[(k + 1, k + 2)] *
      ψ′.AR[k + 2] *
      ψ′.AR[k + 1] for k in 1:Nsites
    ])
    eᴸ = [
      (hᴸ[k] * ψ.C[k] * δ(only(dag(r[k])), only(dag(r′[k]))) * ψ′.C[k])[] for k in 1:Nsites
    ]
    eᴿ = [(hᴿ[k] * ψ.C[k] * δ(only(l[k]), only(l′[k])) * ψ′.C[k])[] for k in 1:Nsites]
    for k in 1:Nsites
      # TODO: remove `denseblocks` once BlockSparse + DiagBlockSparse is supported
      hᴸ[k] -= eᴸ[k] * denseblocks(δ(inds(hᴸ[k])))
      hᴿ[k] -= eᴿ[k] * denseblocks(δ(inds(hᴿ[k])))
    end

    function left_environment_cell(ψ, ψ̃, hᴸ, n)
      Nsites = nsites(ψ)
      𝕙ᴸ = copy(hᴸ)
      for k in reverse((n - Nsites + 2):n)
        𝕙ᴸ[k] = 𝕙ᴸ[k - 1] * ψ.AL[k] * ψ̃.AL[k] + 𝕙ᴸ[k]
      end
      return 𝕙ᴸ[n]
    end

    #for k in 2:Nsites
    #  hᴸ[k] = hᴸ[k - 1] * ψ.AL[k] * ψ̃.AL[k] + hᴸ[k]
    #end
    𝕙ᴸ = copy(hᴸ)
    for k in 1:Nsites
      𝕙ᴸ[k] = left_environment_cell(ψ, ψ̃, hᴸ, k)
    end
    Hᴸ = left_environment(hᴸ, 𝕙ᴸ, ψ; tol=_solver_tol)
    for k in 2:Nsites
      hᴿ[k] = hᴿ[k + 1] * ψ.AR[k + 1] * ψ̃.AR[k + 1] + hᴿ[k]
    end
    Hᴿ = right_environment(hᴿ, ψ; tol=_solver_tol)

    Cvalsₙ₋₁, Cvecsₙ₋₁, Cinfoₙ₋₁ = solver(
      Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n - 1), time_step, ψ.C[n - 1], _solver_tol, eager
    )
    Cvalsₙ, Cvecsₙ, Cinfoₙ = solver(
      Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n), time_step, ψ.C[n], _solver_tol, eager
    )
    Avalsₙ, Avecsₙ, Ainfoₙ = solver(
      Hᴬᶜ(∑h, Hᴸ, Hᴿ, ψ, n), time_step, ψ.AL[n] * ψ.C[n], _solver_tol, eager
    )

    C̃[n - 1] = Cvecsₙ₋₁
    C̃[n] = Cvecsₙ
    Ãᶜ[n] = Avecsₙ

    function ortho_overlap(AC, C)
      AL, _ = polar(AC * dag(C), uniqueinds(AC, C))
      return noprime(AL)
    end

    function ortho_polar(AC, C)
      UAC, _ = polar(AC, uniqueinds(AC, C))
      UC, _ = polar(C, commoninds(C, AC))
      return noprime(UAC) * noprime(dag(UC))
    end

    Ãᴸ[n] = ortho_polar(Ãᶜ[n], C̃[n])
    Ãᴿ[n] = ortho_polar(Ãᶜ[n], C̃[n - 1])

    # Update state for next iteration
    #ψ = InfiniteCanonicalMPS(Ãᴸ, C̃, Ãᴿ)
    ψ.AL[n] = Ãᴸ[n]
    ψ.AR[n] = Ãᴿ[n]
    ψ.C[n - 1] = C̃[n - 1]
    ψ.C[n] = C̃[n]
    ψᴴ = dag(ψ)
    ψ′ = ψᴴ'
    # XXX: make this prime the center sites
    ψ̃ = prime(linkinds, ψᴴ)

    # TODO: replace with linkinds(ψ)
    l = CelledVector([commoninds(ψ.AL[n], ψ.AL[n + 1]) for n in 1:Nsites])
    l′ = CelledVector([commoninds(ψ′.AL[n], ψ′.AL[n + 1]) for n in 1:Nsites])
    r = CelledVector([commoninds(ψ.AR[n], ψ.AR[n + 1]) for n in 1:Nsites])
    r′ = CelledVector([commoninds(ψ′.AR[n], ψ′.AR[n + 1]) for n in 1:Nsites])
  end
  for n in 1:Nsites
    ϵᴸ![n] = norm(Ãᶜ[n] - Ãᴸ[n] * C̃[n])
    ϵᴿ![n] = norm(Ãᶜ[n] - C̃[n - 1] * Ãᴿ[n])
  end
  return ψ, (eᴸ, eᴿ)
end

function tdvp_iteration_parallel(
  solver::Function,
  ∑h::InfiniteSum{ITensor},
  ψ::InfiniteCanonicalMPS;
  (ϵᴸ!)=fill(1e-15, nsites(ψ)),
  (ϵᴿ!)=fill(1e-15, nsites(ψ)),
  time_step,
  solver_tol=(x -> x / 100),
  eager,
)
  Nsites = nsites(ψ)
  range_∑h = nrange(ψ, ∑h[1])
  ϵᵖʳᵉˢ = max(maximum(ϵᴸ!), maximum(ϵᴿ!))
  _solver_tol = solver_tol(ϵᵖʳᵉˢ)
  ψᴴ = dag(ψ)
  ψ′ = ψᴴ'
  # XXX: make this prime the center sites
  ψ̃ = prime(linkinds, ψᴴ)

  # TODO: replace with linkinds(ψ)
  l = CelledVector([commoninds(ψ.AL[n], ψ.AL[n + 1]) for n in 1:Nsites])
  l′ = CelledVector([commoninds(ψ′.AL[n], ψ′.AL[n + 1]) for n in 1:Nsites])
  r = CelledVector([commoninds(ψ.AR[n], ψ.AR[n + 1]) for n in 1:Nsites])
  r′ = CelledVector([commoninds(ψ′.AR[n], ψ′.AR[n + 1]) for n in 1:Nsites])

  hᴸ = InfiniteMPS([
    δ(only(l[n - 2]), only(l′[n - 2])) *
    ψ.AL[n - 1] *
    ψ.AL[n] *
    ∑h[(n - 1, n)] *
    ψ′.AL[n - 1] *
    ψ′.AL[n] for n in 1:Nsites
  ])

  hᴿ = InfiniteMPS([
    δ(only(dag(r[n + 2])), only(dag(r′[n + 2]))) *
    ψ.AR[n + 2] *
    ψ.AR[n + 1] *
    ∑h[(n + 1, n + 2)] *
    ψ′.AR[n + 2] *
    ψ′.AR[n + 1] for n in 1:Nsites
  ])

  eᴸ = [
    (hᴸ[n] * ψ.C[n] * δ(only(dag(r[n])), only(dag(r′[n]))) * ψ′.C[n])[] for n in 1:Nsites
  ]
  eᴿ = [(hᴿ[n] * ψ.C[n] * δ(only(l[n]), only(l′[n])) * ψ′.C[n])[] for n in 1:Nsites]

  for n in 1:Nsites
    # TODO: use these instead, for now can't subtract
    # BlockSparse and DiagBlockSparse tensors
    #hᴸ[n] -= eᴸ[n] * δ(inds(hᴸ[n]))
    #hᴿ[n] -= eᴿ[n] * δ(inds(hᴿ[n]))
    hᴸ[n] -= eᴸ[n] * denseblocks(δ(inds(hᴸ[n])))
    hᴿ[n] -= eᴿ[n] * denseblocks(δ(inds(hᴿ[n])))
  end

  # Sum the Hamiltonian terms in the unit cell
  function left_environment_cell(ψ, ψ̃, hᴸ, n)
    Nsites = nsites(ψ)
    𝕙ᴸ = copy(hᴸ)
    for k in reverse((n - Nsites + 2):n)
      𝕙ᴸ[k] = 𝕙ᴸ[k - 1] * ψ.AL[k] * ψ̃.AL[k] + 𝕙ᴸ[k]
    end
    return 𝕙ᴸ[n]
  end

  #for k in 2:Nsites
  #  hᴸ[k] = hᴸ[k - 1] * ψ.AL[k] * ψ̃.AL[k] + hᴸ[k]
  #end
  𝕙ᴸ = copy(hᴸ)
  for k in 1:Nsites
    𝕙ᴸ[k] = left_environment_cell(ψ, ψ̃, hᴸ, k)
  end
  Hᴸ = left_environment(hᴸ, 𝕙ᴸ, ψ; tol=_solver_tol)

  for n in 2:Nsites
    hᴿ[n] = hᴿ[n + 1] * ψ.AR[n + 1] * ψ̃.AR[n + 1] + hᴿ[n]
  end
  Hᴿ = right_environment(hᴿ, ψ; tol=_solver_tol)

  C̃ = InfiniteMPS(Vector{ITensor}(undef, Nsites), translator(ψ))
  Ãᶜ = InfiniteMPS(Vector{ITensor}(undef, Nsites), translator(ψ))
  for n in 1:Nsites
    Cvalsₙ, Cvecsₙ, Cinfoₙ = solver(
      Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n), time_step, ψ.C[n], _solver_tol, eager
    )
    Avalsₙ, Avecsₙ, Ainfoₙ = solver(
      Hᴬᶜ(∑h, Hᴸ, Hᴿ, ψ, n), time_step, ψ.AL[n] * ψ.C[n], _solver_tol, eager
    )

    C̃[n] = Cvecsₙ
    Ãᶜ[n] = Avecsₙ
  end

  Ãᴸ = InfiniteMPS(Vector{ITensor}(undef, Nsites), translator(ψ))
  Ãᴿ = InfiniteMPS(Vector{ITensor}(undef, Nsites), translator(ψ))
  for n in 1:Nsites
    Ãᴸ[n] = ortho_polar(Ãᶜ[n], C̃[n])
    Ãᴿ[n] = ortho_polar(Ãᶜ[n], C̃[n - 1])
  end

  for n in 1:Nsites
    ϵᴸ![n] = norm(Ãᶜ[n] - Ãᴸ[n] * C̃[n])
    ϵᴿ![n] = norm(Ãᶜ[n] - C̃[n - 1] * Ãᴿ[n])
  end
  return InfiniteCanonicalMPS(Ãᴸ, C̃, Ãᴿ), (eᴸ, eᴿ)
end


text:
"""
    count_nucleotides(strand)

The frequency of each nucleotide within `strand` as a dictionary.

Invalid strands raise a `DomainError`.

"""
function count_nucleotides(strand)
    d = Dict('A' => 0, 'C' => 0, 'G' => 0, 'T' => 0)
    for c in strand
        try d[c] += 1 catch; throw(DomainError("Invalid strand")) end
    end
    return d
end

max_stars_repo_name: aadimator/exercism

max_issues_repo_name: aadimator/exercism

max_forks_repo_name: aadimator/exercism


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

function _check_reshape_varshape_compat(trgshape::AbstractValueShape, srcshape::AbstractValueShape)
    if !(srcshape <= trgshape)
        throw(ArgumentError("Shape of variate doesn't match target variate trgshape, with variate of type $(typeof(v)) and expected trgshape $(trgshape)"))
    end
end
function ChainRulesCore.rrule(::typeof(_check_reshape_varshape_compat), trgshape::AbstractValueShape, srcshape::AbstractValueShape)
    _check_reshape_varshape_compat_pullback(ΔΩ) = (NoTangent(), NoTangent(), ZeroTangent())
    return _check_reshape_varshape_compat(trgshape, srcshape), _check_reshape_varshape_compat_pullback
end

function _check_reshape_ndof_compat(trgshape::AbstractValueShape, srcshape::AbstractValueShape)
    ndof_trg = totalndof(trgshape)
    ndof_src = totalndof(srcshape)
    if ndof_trg != ndof_src
        throw(ArgumentError("Source varshape has ($ndof_src) degrees of freedom but target varshape has $ndof_trg"))
    end
end
function ChainRulesCore.rrule(::typeof(_check_reshape_ndof_compat), trgshape::AbstractValueShape, srcshape::AbstractValueShape)
    _check_reshape_ndof_compat_pullback(ΔΩ) = (NoTangent(), NoTangent(), ZeroTangent())
    return _check_reshape_ndof_compat(trgshape, srcshape), _check_reshape_ndof_compat_pullback
end

reshape_variate(trgshape::Missing, srcshape::Any, v::Any) = throw(ArgumentError("Can't reshape variate without trgshape information"))

function reshape_variate(trgshape::VS, srcshape::VS, v::Any) where {VS<:AbstractValueShape}
    _check_reshape_varshape_compat(trgshape, srcshape)
    v
end

function reshape_variate(trgshape::ScalarShape{Real}, srcshape::ScalarShape{Real}, v::Real)
    v
end

function reshape_variate(trgshape::ArrayShape{<:Real,1}, srcshape::ArrayShape{<:Real,1}, v::AbstractVector{<:Real})
    _check_reshape_varshape_compat(trgshape, srcshape)
    v
end

function reshape_variate(trgshape::ArrayShape{<:Real,1}, srcshape::AbstractValueShape, v::Any)
    _check_reshape_ndof_compat(trgshape, srcshape)
    inverse(srcshape)(v)
end

function reshape_variate(trgshape::AbstractValueShape, srcshape::ArrayShape{<:Real,1}, v::AbstractVector{<:Real})
    _check_reshape_ndof_compat(trgshape, srcshape)
    trgshape(v)
end

function reshape_variates(trgshape::VS, srcshape::VS, vs::AbstractVector{<:Any}) where {VS<:AbstractValueShape}
    _check_reshape_varshape_compat(trgshape, srcshape)
    vs
end

function reshape_variates(trgshape::ScalarShape{Real}, srcshape::ScalarShape{Real}, vs::AbstractVector{Real})
    vs
end

function reshape_variates(trgshape::ArrayShape{<:Real,1}, srcshape::ArrayShape{<:Real,1}, vs::AbstractVector{<:AbstractVector{<:Real}})
    _check_reshape_varshape_compat(trgshape, srcshape)
    vs
end

function reshape_variates(trgshape::ArrayShape{<:Real,1}, srcshape::AbstractValueShape, vs::AbstractVector{<:AbstractVector{<:Real}})
    _check_reshape_ndof_compat(trgshape, srcshape)
    inverse(srcshape).(vs)
end

function reshape_variates(trgshape::AbstractValueShape, srcshape::ArrayShape{<:Real,1}, vs::AbstractVector{<:AbstractVector{<:Real}})
    _check_reshape_ndof_compat(trgshape, srcshape)
    trgshape(vs)
end

function check_variate end

function ChainRulesCore.rrule(::typeof(check_variate), trgshape::Any, v::Any)
    result = check_variate(trgshape, v)
    _check_variate_pullback(ΔΩ) = (NoTangent(), NoTangent(), ZeroTangent())
    return result, _check_variate_pullback
end

check_variate(trgshape::ScalarShape{Real}, v::Real) = nothing

function check_variate(trgshape::ArrayShape{<:Real,N}, v::AbstractArray{<:Real,N}) where N
    ndof = length(eachindex(v))
    ndof_expected = totalndof(trgshape)
    if ndof != ndof_expected
        throw(ArgumentError("Invalid length ($ndof) of variate, target trgshape  $(trgshape) has $ndof_expected degrees of freedom"))
    end
    nothing
end

function check_variate(trgshape::NamedTupleShape, v::ShapedAsNT)
    if !(valshape(v) <= trgshape)
        throw(ArgumentError("Shape of variate incompatible with target variate trgshape, with variate of type $(typeof(v)) and expected trgshape $(trgshape)"))
    end
    nothing
end

function check_variate(trgshape::NamedTupleShape{names}, v::NamedTuple{names}) where names
    nothing
end

function check_variate(trgshape::Any, v::Any)
    throw(ArgumentError("Shape of variate incompatible with target variate trgshape, with variate of type $(typeof(v)) and expected trgshape $(trgshape)"))
end

function check_variate(trgshape::Missing, v::Any)
    throw(ArgumentError("Cannot evaluate without value trgshape information"))
end


text:
__precompile__(true)

module EarthEngine

using Logging
using PyCall
using OrderedCollections
# import Core: String

const module_dir = @__DIR__

# include the Base module methods to extend
include(joinpath(module_dir,"extendfuncs.jl"))

# include the functions to support interactive mapping
# include(joinpath(module_dir,"eemap.jl"))

# define a constant pointer to the ee object
const ee = PyNULL()

version() = VersionNumber(ee.__version__)

const pre_type_map = []

# Maps a python object corresponding to an EE class to a Julia type which
# wraps that class.
const type_map = OrderedDict{PyObject,DataType}()

abstract type AbstractEEObject end

PyCall.PyObject(x::AbstractEEObject) = x.pyo

"""
    @pytype name class

Macro for creating a Julia Type that wraps a PyObject class
"""
macro pytype(name, class)
    quote
        struct $(name) <: AbstractEEObject
            pyo::PyObject
            $(esc(name))(pyo::PyObject) = new(pyo)
            function $(esc(name))(args...; kwargs...)
                ee_method = ($class)()
                new(pycall(ee_method, PyObject, args...; kwargs...))
            end
        end

        # This won't work until PyCall is updated to support
        # the Julia 1.0 iteration protocol.
        function Base.iterate(x::$name, state...)
            res = Base.iterate(x.pyo, state...)
            if res === nothing
                return nothing
            else
                value, state = res
                return ee_wrap(value), state
            end
        end

        push!(pre_type_map, ($class, $name))
    end
end

"""
    ee_wrap(pyo::PyObject)

Function for wrapping a Python object defined in the type map
"""
function ee_wrap(pyo::PyObject)
    for (pyt, pyv) in type_map
        pyt === nothing && continue
        if pyisinstance(pyo, pyt)
            return pyv(pyo)
        end
    end
    return convert(PyAny, pyo)
end

quot(x) = Expr(:quote, x)

ee_wrap(x::Union{AbstractArray,Tuple}) = [ee_wrap(_) for _ in x]

ee_wrap(pyo) = pyo

fix_arg(x) = x

pyattr(class, method) = pyattr(class, method, method)

"""
    pyattr(class method orig_method)

Function for creating a method signature for a Julia Type
In Python world method(T::class) is analagous to class.method()
"""
function pyattr(class, jl_method, py_method)
    quote
        function $(esc(jl_method))(pyt::$class, args...; kwargs...)
            new_args = fix_arg.(args)
            method = pyt.pyo.$(string(py_method))
            pyo = pycall(method, PyObject, new_args...; kwargs...)
            wrapped = ee_wrap(pyo)
        end
    end
end

"""
    @pyattr class method

Macro for creating a method signature for a Julia Type
In Python world method(T::class) is analagous to class.method()
"""
macro pyattr(class, method)
    pyattr(class, method)
end

"""
    @pyattr class method orig_method

Macro for creating a method signature for a Julia Type
In Python world method(T::class) is analagous to class.method()
This will create a new method name which calls the orig_method
"""
macro pyattr(class, method, orig_method)
    pyattr(class, method, orig_method)
end

"""
    pyattr_set(types, methods...)

For each Julia type `T<:AbstractEEObject` in `types` and each method `m` in `methods`,
define a new function `m(t::T, args...)` that delegates to the underlying
pyobject wrapped by `t`.
"""
function pyattr_set(classes, methods...)
    for class in classes
        for method in methods
            @eval @pyattr($class, $method)
        end
    end
end

"""
    eefunc(method)

Function to wrap Julia defined method on the Python side. This is needed when applying
functions to specific types such as EE.List.
"""
function eefunc(method)
    # extract out the callable attribute from python object
    return PyObject(method)."__call__"
end

"""
    @eefunc method

Macro to wrap Julia defined method on the Python side. This is needed when applying
functions to specific types such as EE.List.
"""
macro eefunc(method)
    quote
        eefunc($(esc(method)))
    end
end

"""
    @eefunc method type

Macro to wrap Julia defined method with a specific type signature on the Python side.
This is needed when applying functions with type annotations using `map()`
"""
macro eefunc(method, type)
    quote
        wrapped(pyo::PyObject) = $(esc(method))($(esc(type))(pyo))
        eefunc(wrapped)
    end
end

# constant set to hold unique method names from the EE Python API
const ee_exports = OrderedSet()
# specify a list of unique names to not export
const excludeexports = [:String]

"""
    Initialize(args...; kwargs...)

Function to initialize an EarthEngine session (analagous to ee.Initialize() from
the Python API). Accepts arguments and keywords from the Python ee.Initialize()
function. This function also dynamically builds the EE API and creates the methods
with signatures for each EE Type.
"""
function Initialize(args...; kwargs...)
    try
        copy!(ee, pyimport("ee"))
    catch err
        error(
            "The `earthengine-api` package could not be imported. You must install the Python earthengine-api before using this package. The error was $err",
        )
    end

    try
        ee.Initialize(args...; kwargs...)

    catch
        @warn "Could not initialize an `ee` session. Trying authentication workflow..."

        try
            Authenticate()
            ee.Initialize()
        catch
            error(
                "Could not initialize an `ee` session or run authentication workflow. Please authenticate manually using the earthengine-api CLI (i.e. `\$ earthengine authenticate`",
            )
        end
    end

    # pull in the types and dynamically wrap things after initialization
    include(joinpath(module_dir,"eetypes.jl"))
    include(joinpath(module_dir,"eefuncs.jl"))

    # export the EE methods
    for f in Symbol.(ee_exports)
        if ~(f in excludeexports)
            @eval export $f
        end
    end

end

"""
    Authenticate()

Function to execute the EarthEngine authetication workflow (analgous to
ee.Authenticate() in the Python API). This function should only be executed
once if the EE API has not be used before.
"""
function Authenticate(args...; kwargs...)
    try
        ee.Authenticate(args...; kwargs...)
    catch err
        @error "Could not authenticate `ee`...see following error: $err"
    end
end

# create an abreviated variable for the module
const EE = EarthEngine

export EE, ee, Initialize, @eefunc

end # module


text:
using SymbolicUtils: Symbolic

"""
    @register_symbolic(expr, define_promotion = true, Ts = [Num, Symbolic, Real])

Overload appropriate methods so that Symbolics can stop tracing into the
registered function. If `define_promotion` is true, then a promotion method in
the form of
```julia
SymbolicUtils.promote_symtype(::typeof(f_registered), args...) = Real # or the annotated return type
```
is defined for the register function. Note that when defining multiple register
overloads for one function, all the rest of the registers must set
`define_promotion` to `false` except for the first one, to avoid method
overwritting.

# Examples
```julia
@register_symbolic foo(x, y)
@register_symbolic foo(x, y::Bool) false # do not overload a duplicate promotion rule
@register_symbolic goo(x, y::Int) # `y` is not overloaded to take symbolic objects
@register_symbolic hoo(x, y)::Int # `hoo` returns `Int`
```
"""
macro register_symbolic(expr, define_promotion = true, Ts = [])
    if expr.head === :(::)
        ret_type = expr.args[2]
        expr = expr.args[1]
    else
        ret_type = Real
    end

    @assert expr.head === :call

    f = expr.args[1]
    args = expr.args[2:end]

    types = map(args) do x
        if x isa Symbol
            :(($Real, $wrapper_type($Real), $Symbolic{<:$Real}))
        elseif Meta.isexpr(x, :(::))
            T = x.args[2]
            :($has_symwrapper($T) ?
              ($T, $Symbolic{<:$T}, $wrapper_type($T),) :
              ($T, $Symbolic{<:$T}))
        else
            error("Invalid argument format $x")
        end
    end

    eval_method = :(@eval function $f($(Expr(:$, :(s...))),)
                        args = map($unwrap, [$(Expr(:$, :(s_syms...)))])
                        if !any(x->$issym(x) || $istree(x), args)
                            $f(args...)
                        else
                            $Term{$ret_type}($f, args)
                        end |> $wrap
                    end)
    verbose = false
    mod, fname = f isa Expr && f.head == :(.) ? f.args : (:(@__MODULE__), QuoteNode(f))
    Ts = Symbol("##__Ts")
    quote
        $Ts = [Tuple{x...} for x in Iterators.product($(types...),)
                if !all(x->Symbolics.has_symwrapper(x) &&
                        Symbolics.wrapper_type(x) != x, x)]
        if $verbose
            println("Candidates")
            map(println, $Ts)
        end

        for sig in $Ts
            s = map(((i,T,),)->Expr(:(::), Symbol("arg", i), T), enumerate(sig.parameters))
            s_syms = map(x->x.args[1], s)
            $eval_method
        end
        if $define_promotion
            (::$typeof($promote_symtype))(::$typeof($f), args...) = $ret_type
        end
    end |> esc
end

Base.@deprecate_binding var"@register" var"@register_symbolic"

# Ensure that Num that get @registered from outside the ModelingToolkit
# module can work without having to bring in the associated function into the
# ModelingToolkit namespace. We basically store information about functions
# registered at runtime in a ModelingToolkit variable,
# `registered_external_functions`. It's not pretty, but we are limited by the
# way GeneralizedGenerated builds a function (adding "ModelingToolkit" to every
# function call).
# ---
const registered_external_functions = Dict{Symbol,Module}()
function inject_registered_module_functions(expr)
    MacroTools.postwalk(expr) do x
        # Find all function calls in the expression and extract the function
        # name and calling module.
        MacroTools.@capture(x, f_module_.f_name_(xs__))
        if isnothing(f_module)
            MacroTools.@capture(x, f_name_(xs__))
        end

        if !isnothing(f_name)
            # Set the calling module to the module that registered it.
            mod = get(registered_external_functions, f_name, f_module)
            if !isnothing(mod) && mod != Base
                x.args[1] = :($mod.$f_name)
            end
        end

        return x
    end
end


text:
struct DFP{T1} <: QuasiNewton{T1}
    approx::T1
end
DFP(; inverse = true) = DFP(inverse ? Inverse() : Direct())
hasprecon(::DFP) = NoPrecon()

summary(dfp::DFP{Inverse}) = "Inverse DFP"
summary(dfp::DFP{Direct}) = "Direct DFP"

# function update!(scheme::DFP, B::Inverse, s, y)
#    B.A = B.A + s*s'/dot(s, y) - B.A*y*y'*B.A/(y'*B.A*y)
# end
# function update!(scheme::DFP, B::Direct, s, y)
#    B.A = (I - y*s'/dot(y, s))*B.A*(I - s*y'/dot(y, s)) + y*y'/dot(y, s)
# end

function update(scheme::DFP{<:Inverse}, H, s, y)
    σ = dot(s, y)
    ρ = inv(σ)
    H = H + ρ * s * s' - H * (y * y') * H / (y' * H * y)
    H
end
function update(scheme::DFP{<:Direct}, B, s, y)
    σ = dot(s, y)
    ρ = inv(σ)

    C = (I - ρ * y * s')
    B = C * B * C' + ρ * y * y'
    B
end
function update!(scheme::DFP{<:Inverse}, H, s, y)
    σ = dot(s, y)
    ρ = inv(σ)

    H .+= ρ * s * s' - H * (y * y') * H / (y' * H * y)
    H
end
function update!(scheme::DFP{<:Direct}, B, s, y)
    σ = dot(s, y)
    ρ = inv(σ)

    C = (I - ρ * y * s')
    B .= C * B * C' + ρ * y * y'

    B
end
update!(scheme::DFP{<:Inverse}, A::UniformScaling, s, y) = update(scheme, A, s, y)
update!(scheme::DFP{<:Direct}, A::UniformScaling, s, y) = update(scheme, A, s, y)

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
# CUDA forward compatibility package as taken from the datacenter driver
#
# - https://docs.nvidia.com/deploy/cuda-compatibility/index.html#forward-compatibility-title
# - https://docs.nvidia.com/datacenter/tesla/index.html
# - https://www.nvidia.com/Download/index.aspx

using BinaryBuilder, Pkg

include("../../../fancy_toys.jl")

name = "CUDA_compat"
version = v"11.6"

cuda_version = "$(version.major)-$(version.minor)"
driver_version = "510.39.01"
build = 1

sources_linux_x86 = [
    FileSource("https://developer.download.nvidia.com/compute/cuda/repos/rhel8/x86_64/cuda-compat-$(cuda_version)-$(driver_version)-$(build).x86_64.rpm",
               "348bb6cd66042bd90bc58f7dd909962c3bc9c3eb7a4ef75aa3178aad9fc00d5a", "compat.rpm")
]
sources_linux_ppc64le = [
    FileSource("https://developer.download.nvidia.com/compute/cuda/repos/rhel8/ppc64le/cuda-compat-$(cuda_version)-$(driver_version)-$(build).ppc64le.rpm",
               "7407dec3754ac416f6e10b5448756c463e9ab068f588bd30afbfc8fd600c3b45", "compat.rpm")
]
sources_linux_aarch64 = [
    FileSource("https://developer.download.nvidia.com/compute/cuda/repos/rhel8/sbsa/cuda-compat-$(cuda_version)-$(driver_version)-$(build).aarch64.rpm",
               "7d8eacda7955677839694ffb02f6d1444e2c66a0ff68327e5a82d2f0d8d994ac", "compat.rpm")
]

dependencies = []

script = raw"""
apk add rpm2cpio
rpm2cpio compat.rpm | cpio -idmv

mkdir -p ${libdir}

mv usr/local/cuda-*/compat/* ${libdir}
"""

products = [
    LibraryProduct("libcuda", :libcuda),
    LibraryProduct("libnvidia-ptxjitcompiler", :libnvidia_ptxjitcompiler),
]

non_reg_ARGS = filter(arg -> arg != "--register", ARGS)

if should_build_platform("x86_64-linux-gnu")
    build_tarballs(non_reg_ARGS, name, version, sources_linux_x86, script,
                   [Platform("x86_64", "linux")], products, dependencies;
                   lazy_artifacts=true, skip_audit=true)
end

if should_build_platform("powerpc64le-linux-gnu")
    build_tarballs(non_reg_ARGS, name, version, sources_linux_ppc64le, script,
                   [Platform("powerpc64le", "linux")], products, dependencies;
                   lazy_artifacts=true, skip_audit=true)
end

if should_build_platform("aarch64-linux-gnu")
    build_tarballs(ARGS, name, version, sources_linux_aarch64, script,
                   [Platform("aarch64", "linux")], products, dependencies;
                   lazy_artifacts=true, skip_audit=true)
end

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
var1(x) = length(x) > 1 ? var(x) : 0.0
std1(x) = length(x) > 1 ? std(x) : 0.0
sem1(x) = length(x) > 1 ? sem(x) : 0.0
variation1(x) = length(x) > 1 ? variation(x) : 0.0
skewness1(x) = length(x) > 1 ? skewness(x) : 1.0
kurtosis1(x) = length(x) > 1 ? kurtosis(x) : -1.0

for T1 in (:T, :(Union{Missing,T}))
    for (R,F) in ((:runmin, :minimum), (:runmax, :maximum),
                  (:runmean, :mean), (:runmedian, :median),
                  (:runvar, :var1), (:runstd, :std1),
                  (:runskewness, :skewness1), (:runkurtosis, :kurtosis1),
                  (:runsem, :sem1),
                  (:runmad, :mad_not_normalized),
                  (:runmad_normalized, :mad_normalized),
                  (:runvariation, :variation1))
        @eval begin
            $R(data::V, windowspan::Int) where {T, V<:AbstractVector{$T1}} =
                running($F, data, windowspan)
            $R(data::V, windowspan::Int, weighting::AbstractVector{S}) where {T, V<:AbstractVector{$T1}, S} =
                running($F, data, windowspan, weighting)
            $R(data::V, windowspan::Int, weighting::AbstractWeights) where {T, V<:AbstractVector{$T1}} =
                running($F, data, windowspan, weighting.values)
        end
    end
end

runcor(data1::V1, data2::V2, windowspan::Int) where {T, V1<:Union{AbstractVector{T}, AbstractVector{Union{Missing,T}}},
                                                         V2<:Union{AbstractVector{T},AbstractVector{Union{Missing,T}}}} =
    running(cor, data1, data2, windowspan, 1)

runcov(data1::V1, data2::V2, windowspan::Int) where {T, V1<:Union{AbstractVector{T}, AbstractVector{Union{Missing,T}}},
                                                         V2<:Union{AbstractVector{T},AbstractVector{Union{Missing,T}}}} =
    running(cov, data1, data2, windowspan, 0)

runcor(data1::V1, data2::V2, windowspan::Int, weighting::AbstractVector{S}) where {S, T, V1<:Union{AbstractVector{T}, AbstractVector{Union{Missing,T}}},
                                                                                V2<:Union{AbstractVector{T}, AbstractVector{Union{Missing,T}}}} =
    running(cor, data1, data2, windowspan, weighting, 1)

runcov(data1::V1, data2::V2, windowspan::Int, weighting::AbstractVector{S}) where {S, T, V1<:Union{AbstractVector{T}, AbstractVector{Union{Missing,T}}},
                                                                                V2<:Union{AbstractVector{T}, AbstractVector{Union{Missing,T}}}} =
    running(cov, data1, data2, windowspan, weighting, 0)


text:
# """
#     MatrixProductOperator

# `lx` is the length of the MPO or the number of sites.  `d` is the
# physical on-site dimension, `dims` is the dimenions of bonds, starting
# and ending with 1.  `tensors` is the local tensors, in four
# dimenionsional tensor format, where the indeces are ordered starting
# from the left one using the counterclockwise convention for the
# indeces, so left, down, right, up.

# """
# struct MatrixProductOperator{T<:RLorCX}
#     lx :: Int
#     d :: Int
#     dims :: Vector{Int64}
#     tensors :: Vector{Array{T,4}}
# end

# struct InfMatrixProductOperator{T<:Number, OP<:Union{Array{T, 4}, SymTensor{T, 4}}}
#     ly :: Int
#     ds :: Vector{Int}
#     Ws :: Vector{OP}
# end
# ### constructors
# ################

function xxz_mpo(T::DataType, lx::Int64, d::Int64, delta::Float64=1.0) #where {T<:RLorCX}

    Sz, Sp, Sm = spinoperators(1/2)
    I2 = I(2)

    W = zeros(T, 5, d, 5, d)
    W[1,:,1,:] = I2
    W[2,:,1,:] = 0.5 * Sp
    W[3,:,1,:] = 0.5 * Sm
    W[4,:,1,:] = delta * Sz

    W[5,:,2,:] = Sm
    W[5,:,3,:] = Sp
    W[5,:,4,:] = Sz
    W[5,:,5,:] = I2

    mpo = MPO{T}()
    push!(mpo, W[5:5,:,:,:])
    for site=2:lx-1
        push!(mpo, W)
    end
    push!(mpo, W[:,:,1:1,:])

    mpo
end

function xxzlong_mpo(T::DataType, lx::Int64, d::Int64, delta::Float64=1.0, r::Float64=0.5)
    Sz, Sp, Sm = spinoperators(1/2)
    I2 = I(2)

    W = zeros(T, 5, d, 5, d)
    W[1,:,1,:] = I2
    W[2,:,1,:] = 0.5 * Sp
    W[3,:,1,:] = 0.5 * Sm
    W[4,:,1,:] = delta * Sz

    W[5,:,2,:] = Sm
    W[5,:,3,:] = Sp
    W[5,:,4,:] = Sz
    W[5,:,5,:] = I2

    W[2,:,2,:] = r.*I2
    W[3,:,3,:] = r.*I2
    W[4,:,4,:] = r.*I2

    mpo = MPO{T}()
    push!(mpo, W[5:5,:,:,:])
    for site=2:lx-1
        push!(mpo, W)
    end
    push!(mpo, W[:,:,1:1,:])

    mpo
end

function qitf_mpo(T::DataType, lx::Int64, d::Int64,
                  g::Float64=1.0, J::Float64=-1.0) #where {T<:RLorCX}

    sz, sp ,sm = spinoperators(1/2)
    X = sp + sm
    Z = 2 * sz

    W = zeros(T, 3, d, 3, d)
    W[1,:,1,:] = I(2)
    W[2,:,1,:] = J * Z
    W[3,:,1,:] = g * X

    W[3,:,2,:] = Z
    W[3,:,3,:] = I(2)

    mpo = MPO{T}()
    push!(mpo, W[3:3,:,:,:])
    for site=2:lx-1
        push!(mpo, W)
    end
    push!(mpo, W[:,:,1:1,:])

    mpo
end

function j1j2_mpo(lx::Int64, j1::T, j2::T, d::Int64=2) where {T<:Number}
    W = zeros(T, 8, d, 8, d)
    sz, sp, sm = spinoperators(1/2)
    Id = I(d)
    W[1,:,1,:] = Id
    W[2,:,1,:] = sz
    W[3,:,1,:] = sp
    W[4,:,1,:] = sm

    W[5,:,2,:] = Id
    W[6,:,3,:] = Id
    W[7,:,4,:] = Id

    W[8,:,2,:] = j1 * sz
    W[8,:,3,:] = j1 * 0.5 * sm
    W[8,:,4,:] = j1 * 0.5 * sp
    W[8,:,5,:] = j2 * sz
    W[8,:,6,:] = j2 * 0.5 * sm
    W[8,:,7,:] = j2 * 0.5 * sp
    W[8,:,8,:] = Id

    mpo = MPO{T}()
    push!(mpo, W[8:8,:,:,:])
    for site=2:lx-1
        push!(mpo, W)
    end
    push!(mpo, W[:,:,1:1,:])

    mpo
end

# ### conversions
# ###############

# function convert(::Type{MatrixProductOperator{ComplexF64}},
#                  mpo::MatrixProductOperator{Float64})
#     MatrixProductOperator{ComplexF64}(mpo.lx, mpo.d, mpo.dims,
#                                       convert(Vector{Array{ComplexF64,4}}, mpo.tensors))
# end

# function mpo2hamiltonian(mpo::MatrixProductOperator)
#     lx = mpo.lx
#     d = mpo.d
#     mpo.d^lx > 10000 && error("model is too large for explicit Hamiltonian!")
#     L = mpo.tensors[1]
#     for two in mpo.tensors[2:end]
#         @tensor L[l,o1,o2,r,o1',o2'] := L[l,o1,m,o1'] * two[m, o2, r, o2']
#         dd = size(L,2)*d
#         L = reshape(L, 1, dd, size(two,3), dd)
#     end
#     reshape(L, d^lx, d^lx)
# end

# function reducempo!(mpo::MatrixProductOperator)
#     A = mpo.tensors[1]
#     for l=1:mpo.lx-1
#         U, S, Vt = svdtrunc(reshape(permutedims(A, [1,2,4,3]),
#                                      prod(size(A)[[1,2,4]]), size(A, 3)))
#         mpo.dims[l+1] = size(S, 1)
#         mpo.tensors[l] = permutedims(reshape(U*S,
#                                                size(A)[[1,2,4]]..., size(S, 1)),
#                                      [1,2,4,3])
#         @tensor A[l,o,r,o'] := Vt[l, m] * mpo.tensors[l+1][m,o,r,o']
#     end

#     for l=mpo.lx:-1:2
#         U, S, Vt = svdtrunc(reshape(A, size(A, 1), prod(size(A)[2:4])))
#         mpo.dims[l] = size(S, 1)
#         mpo.tensors[l] = reshape(S*Vt, size(S, 1), size(A)[2:4]...)

#         @tensor A[l,o,r,o'] := mpo.tensors[l-1][l,o,m,o'] *  U[m, r]
#     end
#     mpo.tensors[1] = A
#     mpo
# end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
#####################################################################
# IMPLEMENTACIÓN DEL ALGORITMOS DE INICIALIZACIÓN DE PESOS
# Basado en capítulo 7, algoritmo 6
#   CASO ENTRADA REAL (una dimensión) SALIDA REAL (una dimensión)
#####################################################################
export nn_from_data
"""
nn_from_data(X_train::Vector,Y_train::Vector, n::Int, M=10)::AbstractOneLayerNeuralNetwork
    Devuelve una red neuronal de entrada de una dimensión y
    de salida una dimensión con los pesos ya inicializados
    de acorte a los conjuntos de entrenamiento.
    `n` es el número de neuronas en la capa oculta.
    El tamaño de entrada  y salida de la red neuronal vienen determinados por
    por los propios datos de entranamiento.
    El tamaño de entrada se corresponde con el número de atributos de X_train (su número de columnas)
    y a salida con el número de columnas de `Y_train`.

    M es una constante que depende de la función de activación,
    por lo visto en teoría M=10 funciona para todas las

"""
function nn_from_data(X_train::Vector,Y_train::Vector, n::Int, M=10)::AbstractOneLayerNeuralNetwork
    entry_dimension =  1
    output_dimension = 1

    nodes = []
    index = 1
    tam = 0
    y_values = zeros(n)

    while tam < n && index <= n
        if !(X_train[index] in nodes)
            append!(nodes, X_train[index] )
            tam += 1
            y_values[tam] = Y_train[index]
        end
        index += 1
    end
    ordered_index = sortperm(nodes)
    # Matrices de la red neuronal
    # A = n x 1
    # S = n x 1
    # B = 1 x n
    A = zeros(Float64, n, entry_dimension)
    S = zeros(Float64, n)
    B = zeros(Float64, output_dimension, n)

    # valores iniciales
    x_a = nodes[ordered_index[1]]
    y_a = y_values[ordered_index[1]]
    # Función afín constantemente Y_1
    S[1]= M
    B[1] = y_a

    # Cálculo del resto de neuronas
    for (index,key) in collect(Iterators.zip(2:n, ordered_index[2:n]))
        x_s = nodes[key]
        y_s = y_values[key]

        A[index] = 2M / (x_s - x_a)
        S[index] = M - x_s * A[index]
        B[index] = y_s - y_a

        x_a = x_s
        y_a = y_s

    end
    return FromMatrixNN(S,A,B)
end


text:
module BurglaryModel

include("../src/Jaynes.jl")
using .Jaynes
using Gen

jmodel = @jaynes function burglary_model()::Nothing
    burglary ~ Bernoulli(0.05)
    burglary ? disabled ~ Bernoulli(0.1) : disabled = false
    !disabled ? alarm ~ Bernoulli(burglary ? 0.94 : 0.01) : alarm = false
    call ~ Bernoulli(alarm ? 0.70 : 0.05)
    return nothing
end (NoStatic)

display(jmodel)

# Inference.
chm = choicemap((:call, true))
trs, lnws, lmle = importance_sampling(jmodel, (), chm, 5000)
est = sum(map(zip(lnws, trs)) do (lnw, tr)
              (tr[:burglary] ? 1 : 0) * exp(lnw)
          end)
println(est)

end # module


text:
using QuantumOptics

b = SpinBasis(1//2);
sx = 1/2*sigmax(b);
sy = 1/2*sigmay(b);
sz = 1/2*sigmaz(b);
XX = sx ⊗ sx ;
YY = sy ⊗ sy ;
ZZ = sz ⊗ sz;

Jxy=1
Jz=1

function basis(N)
    spinchain=b
    for ispin in 2:N
        spinchain = spinchain ⊗ b

    end
    return spinchain; # same local
end

function CrtModelH(N)
    Htwoint=Jxy*(XX+YY)+Jz*ZZ
    HXXZ = embed(basis(N), [1,2] , Htwoint)
    for s in 2:N-1
        HXXZ += embed(basis(N), [s,s+1] , Htwoint)
    end
    return HXXZ;
end

s5z = embed(basis(10) , [5] , sz);
s6z = embed(basis(10) , [6] , sz);
sup = spinup(b);
sdo = spindown(b);
psi0 = sup
for ispin in 2:10
    global psi0
    if ispin == 5
       psi0 = psi0 ⊗ sdo;
    else
       psi0 = psi0 ⊗ sup;
    end
end

times = 0:.1:12;
tout, psi_t = timeevolution.schroedinger(times, psi0, CrtModelH(10));

m5z = expect(s5z, psi_t);
m6z = expect(s6z, psi_t);

using PyPlot

PyPlot.plot(tout, m5z,linestyle = ":",linewidth = 5)
PyPlot.plot(tout, m6z,linestyle = ":",linewidth = 5)

max_stars_repo_name: BlackBodyRadiation/JupyterTNN

max_issues_repo_name: BlackBodyRadiation/JupyterTNN

max_forks_repo_name: BlackBodyRadiation/JupyterTNN


text:
export Fixed, Periodic
########################
# Abstract Boundaries
########################

"""
    LatticeTag

Abstract type for lattice properties can be
determined in compile time.
"""
abstract type LatticeTag <: AbstractTag end

"""
    Boundary <: LatticeTag

Abstract type for boundary conditions.
"""
abstract type Boundary <: LatticeTag end

"""
    Periodic <: Boundary

Periodic boundary tag.
"""
abstract type Periodic <: Boundary end

"""
    Fixed <: Boundary

Fixed boundary tag.
"""
abstract type Fixed <: Boundary end

# Square Bonds Iterator
export Vertical, Horizontal, UpRight, UpLeft

abstract type LatticeIterTag end
abstract type Vertical{N} <: LatticeIterTag end
abstract type Horizontal{N} <: LatticeIterTag end
abstract type UpRight{N} <: LatticeIterTag end
abstract type UpLeft{N} <: LatticeIterTag end

export AbstractRegion, Region
# Regions

"""
    AbstractRegion

abstract region.
"""
abstract type AbstractRegion end

"""
    Region{ID} <: AbstractRegion

denotes the `ID`th Region.
"""
abstract type Region{ID} <: AbstractRegion end


text:
using Test
using ToyECDSA

const key = genkey()
const msg = "Bill says this is an elliptic curve digital signature algorithm."
const altered = "Bill says this isn't an elliptic curve digital signature algorithm."

publickey, c, d = ECDSA_sign(msg, key)

@test isverifiedECDSA(msg, publickey, c, d)) == true

@test isverifiedECDSA(altered, publickey, c, d) == false


text:
using Plots
pyplot()

#---------------------------------------------------------------------#
#This code plots the Stability of for the 1D Wave Equation.
#Written by F.X. Giraldo on 03/2019
#           Department of Applied Maths
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#Variables:
#---------------------------------------------------------------------#
#method=input(" 1=FE & Upwind \n 2=FE & Centered \n 3=FE & Downwind \n 4=BE & Upwind \n 5=BE & Centered \n 6=BE & Downwind \n 7=TR & Upwind \n 8=TR & Centered \n 9=TR & Downwind \n Enter Method: ")
method=1

n=100
dtheta=2*pi/n
theta=0:dtheta:2*pi
if (method <= 3)
    cmax=1
elseif [method >= 4]
    cmax=10
end
m=10*cmax
dc=cmax/m
c=0:dc:cmax
em1=exp.(-im*theta)
ep1=exp.(+im*theta)

a=ones(size(em1))
#Plot Stability of Specific Equations
xmatrix=Array{Float64}(undef,n,m+2)
ymatrix=Array{Float64}(undef,n,m+2)
if method == 1 #FE + Upwind
    for j=1:m+1
        z=a - c[j]*(a - em1) #1st order upwind [explicit 1st order in time]
        xmatrix[:,j]=real(z)
        ymatrix[:,j]=imag(z)
    end
elseif method == 2 #FE + Centered
    for j=2:m+1
        z=a - 0.5*c[j]*(ep1 - em1); #2nd order centered [explicit 1st order in time]
        xmatrix[:,j]=real(z)
        ymatrix[:,j]=imag(z)
    end
elseif method == 3 #FE + Downwind
    for j=2:m+1
        z=a - c[j]*(ep1 - a); #1st order downwind [explicit 1st order in time]
        xmatrix[:,j]=real(z)
        ymatrix[:,j]=imag(z)
    end
elseif method == 4 #BE + Upwind
    for j=2:m+1
        z=a/(a + c[j]*(a - em1)) #1st order uwpind, (implicit 1st order in time)
        xmatrix[:,j]=real(z)
        ymatrix[:,j]=imag(z)
    end
elseif method == 5 #BE + Centered
    for j=2:m+1
        z=a/(a + 0.5*c[j]*(ep1 - em1)) #2nd order centered, (implicit 1st order in time)
        xmatrix[:,j]=real(z)
        ymatrix[:,j]=imag(z)
    end
elseif method == 6 #BE + Downwind
    for j=2:m+1
        z=a/(a + c[j]*(ep1 - a)) #1st order downwind, (implicit 1st order in time)
        xmatrix[:,j]=real(z)
        ymatrix[:,j]=imag(z)
    end
elseif method == 7 #TR + Upwind
    for j=2:m+1
        z=(a - 0.5*c[j]*(a - em1))./(a + 0.5*c[j]*(a - em1)) #1st order upwind [Trapezoidal 2nd order in time]
        xmatrix[:,j]=real(z)
        ymatrix[:,j]=imag(z)
    end
elseif method == 8 #TR + Centered
    for j=2:m+1
        z=(a - 0.25*c[j]*(ep1 - em1))./(a + 0.25*c[j]*(ep1 - em1)) #2nd order centered [TR 2nd order in time]
        xmatrix[:,j]=real(z)
        ymatrix[:,j]=imag(z)
    end
elseif method == 9 #TR + Downwind
    for j=2:m+1
        z=(a - 0.5*c[j]*(ep1 - a))./(a + 0.5*c[j]*(ep1 - a)) #1st order downwind [TR 2nd order in time]
        xmatrix[:,j]=real(z)
        ymatrix[:,j]=imag(z)
    end
end

#Plot Unit Circle
z=ep1
xmatrix[:,m+2]=real(z)
ymatrix[:,m+2]=imag(z)
#plot_handle=plot(xmatrix,ymatrix,leg=false,w=2,aspect_ratio=1)
plot(xmatrix,ymatrix,leg=false,w=1,line=(:dash,1),aspect_ratio=1)
plot!(real(z),imag(z),leg=false,w=1,color=[:black],aspect_ratio=1)
xaxis!("Re(z)")
yaxis!("Im(z)")
title!("Von Neumann Stability Analysis")
#display(plot(plot_handle))


text:
using SymbolicUtils: Rewriters

const KEEP = typemin(Int)

function alias_elimination(sys)
    # FIXME: update `structure` too
    #sys = flatten(sys)
    #s = get_structure(sys)
    #if !(s isa SystemStructure)
    sys = initialize_system_structure(sys)
    s = structure(sys)
    #end
    is_linear_equations, eadj, cadj = find_linear_equations(sys)

    v_eliminated, v_types, n_null_vars, degenerate_equations, linear_equations = alias_eliminate_graph(
        s, is_linear_equations, eadj, cadj
    )

    s = structure(sys)
    @unpack fullvars, graph = s

    n_reduced_states = length(v_eliminated) - n_null_vars
    reduced_states = similar(v_eliminated, Any, n_reduced_states)
    subs = OrderedDict()
    if n_reduced_states > 0
        for (i, v) in enumerate(@view v_eliminated[n_null_vars+1:end])
            reduced_states[i] = fullvars[v]
            subs[fullvars[v]] = iszeroterm(v_types, v) ? 0.0 :
                                isalias(v_types, v) ? fullvars[alias(v_types, v)] :
                                -fullvars[negalias(v_types, v)]
        end
    end

    dels = Set{Int}()
    eqs = copy(equations(sys))
    for (ei, e) in enumerate(linear_equations)
        vs = 𝑠neighbors(graph, e)
        if isempty(vs)
            push!(dels, e)
        else
            rhs = 0
            for vj in eachindex(vs)
                var = fullvars[vs[vj]]
                rhs += cadj[ei][vj] * var
            end
            eqs[e] = 0 ~ rhs
        end
    end
    dels = sort(collect(dels))
    deleteat!(eqs, dels)

    dict = Dict(subs)
    for (ieq, eq) in enumerate(eqs)
        if !isdiffeq(eq) && !_iszero(eq.lhs)
            eq = 0 ~ eq.rhs - eq.lhs
        end
        eqs[ieq] = eq.lhs ~ fixpoint_sub(eq.rhs, dict)
    end

    newstates = []
    sts = states(sys)
    for j in eachindex(fullvars)
        if isirreducible(v_types, j)
            isdervar(s, j) || push!(newstates, fullvars[j])
        end
    end

    @set! sys.eqs = eqs
    @set! sys.states = newstates
    @set! sys.reduced_states = [get_reduced_states(sys); reduced_states]
    @set! sys.observed = [get_observed(sys); [lhs ~ rhs for (lhs, rhs) in pairs(subs)]]
    @set! sys.structure = nothing
    return sys
end

function alias_eliminate_graph(s::SystemStructure, is_linear_equations, eadj, cadj)
    @unpack graph, varassoc = s
    invvarassoc = inverse_mapping(varassoc)

    old_cadj = map(copy, cadj)

    is_not_potential_state = iszero.(varassoc)
    is_linear_variables = copy(is_not_potential_state)
    for i in 𝑠vertices(graph); is_linear_equations[i] && continue
        for j in 𝑠neighbors(graph, i)
            is_linear_variables[j] = false
        end
    end
    solvable_variables = findall(is_linear_variables)

    linear_equations = findall(is_linear_equations)

    rank1 = bareiss!(
        (eadj, cadj),
        old_cadj, linear_equations, is_linear_variables, 1
       )

    v_solved = [eadj[i][1] for i in 1:rank1]
    v_eliminated = setdiff(solvable_variables, v_solved)
    n_null_vars = length(v_eliminated)

    v_types = fill(KEEP, ndsts(graph))
    for v in v_eliminated
        v_types[v] = 0
    end

    rank2 = bareiss!(
        (eadj, cadj),
        old_cadj, linear_equations, is_not_potential_state, rank1+1
       )

    rank3 = bareiss!(
        (eadj, cadj),
        old_cadj, linear_equations, nothing, rank2+1
       )

    # kind of like the backward substitution
    for ei in reverse(1:rank2)
        locally_structure_simplify!(
                                    (eadj[ei], old_cadj[ei]),
                                    invvarassoc, v_eliminated, v_types
                                   )
    end

    reduced = false
    for ei in 1:rank2
        if length(cadj[ei]) >= length(old_cadj[ei])
            cadj[ei] = old_cadj[ei]
        else
            # MEMORY ALIAS of a vector
            eadj[ei] = 𝑠neighbors(graph, linear_equations[ei])
            reduced |= locally_structure_simplify!(
                                                   (eadj[ei], cadj[ei]),
                                                   invvarassoc, v_eliminated, v_types
                                                  )
        end
    end

    while reduced
        for ei in 1:rank2
            if !isempty(eadj[ei])
                reduced |= locally_structure_simplify!(
                                                       (eadj[ei], cadj[ei]),
                                                       invvarassoc, v_eliminated, v_types
                                                      )
                reduced && break # go back to the begining of equations
            end
        end
    end

    for ei in rank2+1:length(linear_equations)
        eadj[ei] = old_cadj[ei]
    end

    for (ei, e) in enumerate(linear_equations)
        graph.fadjlist[e] = eadj[ei]
    end

    degenerate_equations = rank3 < length(linear_equations) ? linear_equations[rank3+1:end] : Int[]
    return v_eliminated, v_types, n_null_vars, degenerate_equations, linear_equations
end

iszeroterm(v_types, v) = v_types[v] == 0
isirreducible(v_types, v) = v_types[v] == KEEP
isalias(v_types, v) = v_types[v] > 0 && !isirreducible(v_types, v)
alias(v_types, v) = v_types[v]
negalias(v_types, v) = -v_types[v]

function locally_structure_simplify!(
        (vars, coeffs),
        invvarassoc, v_eliminated, v_types
       )
    while length(vars) > 1 && any(!isequal(KEEP), (v_types[v] for v in @view vars[2:end]))
        for vj in 2:length(vars)
            v = vars[vj]
            if isirreducible(v_types, v)
                continue
            elseif iszeroterm(v_types, v)
                deleteat!(vars, vj)
                deleteat!(coeffs, vj)
                break
            else
                coeff = coeffs[vj]
                if isalias(v_types, v)
                    v = alias(v_types, v)
                else
                    v = negalias(v_types, v)
                    coeff = -coeff
                end

                has_v = false
                for vi in 2:length(vars)
                    (vi !== vj && vars[vi] == v) || continue
                    has_v = true
                    c = (coeffs[vi] += coeff)
                    if c == 0
                        if vi < vj
                            deleteat!(vars, [vi, vj])
                            deleteat!(coeffs, [vi, vj])
                        else
                            deleteat!(vars, [vj, vi])
                            deleteat!(coeffs, [vj, vi])
                        end
                    end
                    break
                end # for vi

                if has_v
                    break
                else
                    vars[vj] = v
                    coeffs[vj] = coeff
                end # if
            end # else
        end # for
    end # while

    v = first(vars)
    if invvarassoc[v] == 0
        if length(vars) == 1
            push!(v_eliminated, v)
            v_types[v] = 0
            empty!(vars); empty!(coeffs)
            return true
        elseif length(vars) == 2 && abs(coeffs[1]) == abs(coeffs[2])
            if (coeffs[1] > 0 && coeffs[2] < 0) || (coeffs[1] < 0 && coeffs[2] > 0)
                # positive alias
                push!(v_eliminated, v)
                v_types[v] = vars[2]
            else
                # negative alias
                push!(v_eliminated, v)
                v_types[v] = -vars[2]
            end
            empty!(vars); empty!(coeffs)
            return true
        end
    end
    return false
end

"""
$(SIGNATURES)

Use Bareiss algorithm to compute the nullspace of an integer matrix exactly.
"""
function bareiss!(
        (eadj, cadj),
        old_cadj, linear_equations, is_linear_variables, offset
       )
    m = length(eadj)
    # v = eadj[ei][vj]
    v = ei = vj = 0
    pivot = last_pivot = 1
    tmp_incidence = Int[]
    tmp_coeffs = Int[]

    for k in offset:m
        ###
        ### Pivoting:
        ###
        ei, vj = find_first_linear_variable(eadj, k:m, is_linear_variables, isequal(1))
        if vj == 0
            ei, vj = find_first_linear_variable(eadj, k:m, is_linear_variables, isequal(2))
        end
        if vj == 0
            ei, vj = find_first_linear_variable(eadj, k:m, is_linear_variables, _->true)
        end

        if vj > 0 # has a pivot
            pivot = old_cadj[ei][vj]
            deleteat!(old_cadj[ei] , vj)
            v = eadj[ei][vj]
            deleteat!(eadj[ei], vj)
            if ei != k
                swap!(cadj, ei, k)
                swap!(old_cadj, ei, k)
                swap!(eadj, ei, k)
                swap!(linear_equations, ei, k)
            end
        else # rank deficient
            return k-1
        end

        for ei in k+1:m
            # elimate `v`
            coeff = 0
            ivars = eadj[ei]
            vj = findfirst(isequal(v), ivars)
            if vj === nothing # `v` is not in in `e`
                continue
            else # remove `v`
                coeff = old_cadj[ei][vj]
                deleteat!(old_cadj[ei], vj)
                deleteat!(eadj[ei], vj)
            end

            # the pivot row
            kvars = eadj[k]
            kcoeffs = old_cadj[k]
            # the elimination target
            ivars = eadj[ei]
            icoeffs = old_cadj[ei]

            empty!(tmp_incidence)
            empty!(tmp_coeffs)
            vars = union(ivars, kvars)

            for v in vars
                ck = getcoeff(kvars, kcoeffs, v)
                ci = getcoeff(ivars, icoeffs, v)
                ci = (pivot*ci - coeff*ck) ÷ last_pivot
                if ci !== 0
                    push!(tmp_incidence, v)
                    push!(tmp_coeffs, ci)
                end
            end

            eadj[ei], tmp_incidence = tmp_incidence, eadj[ei]
            old_cadj[ei], tmp_coeffs = tmp_coeffs, old_cadj[ei]
        end
        last_pivot = pivot
        # add `v` in the front of the `k`-th equation
        pushfirst!(eadj[k], v)
        pushfirst!(old_cadj[k], pivot)
    end

    return m # fully ranked
end

swap!(v, i, j) = v[i], v[j] = v[j], v[i]

function getcoeff(vars, coeffs, var)
    for (vj, v) in enumerate(vars)
        v == var && return coeffs[vj]
    end
    return 0
end

"""
$(SIGNATURES)

Find the first linear variable such that `𝑠neighbors(adj, i)[j]` is true given
the `constraint`.
"""
@inline function find_first_linear_variable(
        eadj,
        range,
        mask,
        constraint,
    )
    for i in range
        vertices = eadj[i]
        if constraint(length(vertices))
            for (j, v) in enumerate(vertices)
                (mask === nothing || mask[v]) && return i, j
            end
        end
    end
    return 0, 0
end

function inverse_mapping(assoc)
    invassoc = zeros(Int, length(assoc))
    for (i, v) in enumerate(assoc)
        v <= 0 && continue
        invassoc[v] = i
    end
    return invassoc
end

"""
$(SIGNATURES)

Use Kahn's algorithm to topologically sort observed equations.

Example:
```julia
julia> @variables t x(t) y(t) z(t) k(t)
(t, x(t), y(t), z(t), k(t))

julia> eqs = [
           x ~ y + z
           z ~ 2
           y ~ 2z + k
       ];

julia> ModelingToolkit.topsort_equations(eqs, [x, y, z, k])
3-element Vector{Equation}:
 Equation(z(t), 2)
 Equation(y(t), k(t) + 2z(t))
 Equation(x(t), y(t) + z(t))
```
"""
function topsort_equations(eqs, states; check=true)
    graph, assigns = observed2graph(eqs, states)
    neqs = length(eqs)
    degrees = zeros(Int, neqs)

    for 𝑠eq in 1:length(eqs); var = assigns[𝑠eq]
        for 𝑑eq in 𝑑neighbors(graph, var)
            # 𝑠eq => 𝑑eq
            degrees[𝑑eq] += 1
        end
    end

    q = Queue{Int}(neqs)
    for (i, d) in enumerate(degrees)
        d == 0 && enqueue!(q, i)
    end

    idx = 0
    ordered_eqs = similar(eqs, 0); sizehint!(ordered_eqs, neqs)
    while !isempty(q)
        𝑠eq = dequeue!(q)
        idx+=1
        push!(ordered_eqs, eqs[𝑠eq])
        var = assigns[𝑠eq]
        for 𝑑eq in 𝑑neighbors(graph, var)
            degree = degrees[𝑑eq] = degrees[𝑑eq] - 1
            degree == 0 && enqueue!(q, 𝑑eq)
        end
    end

    (check && idx != neqs) && throw(ArgumentError("The equations have at least one cycle."))

    return ordered_eqs
end

function observed2graph(eqs, states)
    graph = BipartiteGraph(length(eqs), length(states))
    v2j = Dict(states .=> 1:length(states))

    # `assigns: eq -> var`, `eq` defines `var`
    assigns = similar(eqs, Int)

    for (i, eq) in enumerate(eqs)
        lhs_j = get(v2j, eq.lhs, nothing)
        lhs_j === nothing && throw(ArgumentError("The lhs $(eq.lhs) of $eq, doesn't appear in states."))
        assigns[i] = lhs_j
        vs = vars(eq.rhs)
        for v in vs
            j = get(v2j, v, nothing)
            j !== nothing && add_edge!(graph, i, j)
        end
    end

    return graph, assigns
end

function fixpoint_sub(x, dict)
    y = substitute(x, dict)
    while !isequal(x, y)
        y = x
        x = substitute(y, dict)
    end

    return x
end

function substitute_aliases(eqs, dict)
    sub = Base.Fix2(fixpoint_sub, dict)
    map(eq->eq.lhs ~ sub(eq.rhs), eqs)
end


text:
function get_statuses_mentions_timeline(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/statuses/mentions_timeline.json", options))
end

function get_statuses_user_timeline(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/statuses/user_timeline.json", options))
end

function get_statuses_home_timeline(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/statuses/home_timeline.json", options))
end

function get_statuses_retweets_of_me(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://dev.twitter.com/rest/reference/get/statuses/retweets_of_me", options))
end

# Requires id
function get_statuses_retweets_id(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/statuses/retweets/$id.json", options))
end

# Requires id
function get_statuses_show_id(auth::TwAuth, options::Dict, id::Int64)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/statuses/show.json", options))
end

# Requires id
function get_statuses_oembed(auth::TwAuth, options::Dict, format="json")
  format = lowercase(format)
  if format == "json"
    return   Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/statuses/oembed.$format", options))
  elseif format == "xml"
    return get_oauth(auth, "https://api.twitter.com/1.1/statuses/oembed.$format", options)
  else
    throw(ArgumentError("Format must be either 'json' or 'xml'. User entered '$format'"))
  end
end

# Requires id
function get_statuses_retweeters_ids(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/statuses/retweeters/ids.json", options))
end

# Requires id
function get_statuses_lookup(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/statuses/lookup.json", options))
end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter