text:
"""

    assemble1MolGvacinf!(mySingleMolData::OneMol{FT}, myPeriodicData::PeriodicData{FT},
                         freq::Union{FT, Complex{FT}},
                         blochk::Array{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

Assemble the Green's function interaction matrix of a single molecule in vacuum,
`mySingleMolData.GvacGGinf`, in-place, at frequency `freq` and wavevector `blochk` using
additional information from `myPeriodicData` to speed up the computation. In particular,
if the periodic dimensionality is 0, then this matrix is symmetric (whether `freq` is
imaginary or not). Even if the periodic dimensionality is not 0, if `freq` is strictly
imaginary, then this matrix is Hermitian. Both of those cases allow for computing only the
lower triangle of 3-by-3 blocks. Otherwise, the full matrix must be assembled.

!!! warning

    This function depends on `mySingleMolData.alpha0` having been initialized to the desired
    atomic susceptibilities at the same `freq` and `blochk`, but does not check this.

!!! note

    This function multiplies the Green's function tensor for every pair of atoms by the
    charges of each associated electronic oscillator.

"""
function assemble1MolGvacinf!(mySingleMolData::OneMol{FT}, myPeriodicData::PeriodicData{FT},
                              freq::Union{FT, Complex{FT}},
                              blochk::Array{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

    fill!(mySingleMolData.GvacGGinf, zero(eltype(mySingleMolData.GvacGGinf)));
    Gtemp = Array{eltype(mySingleMolData.GvacGGinf), 2}(undef, 3, 3);
    if (myPeriodicData.numdims == 0)

        ## @inbounds @sync @distributed for qq=1:mySingleMolData.numatoms
        @inbounds for qq=1:mySingleMolData.numatoms
            @inbounds for pp=qq:mySingleMolData.numatoms

                GFVACGG!(Gtemp, 1, 1, freq, view(mySingleMolData.atompos0, 1:3, pp),
                        view(mySingleMolData.atompos0, 1:3, qq),
                        mySingleMolData.alpha0[pp], mySingleMolData.alpha0[qq]);
                Gtemp .= mySingleMolData.Qe[pp] * mySingleMolData.Qe[qq] * Gtemp;
                mySingleMolData.GvacGGinf[3*pp-2:3*pp, 3*qq-2:3*qq] .= Gtemp;
                mySingleMolData.GvacGGinf[3*qq-2:3*qq, 3*pp-2:3*pp] .= transpose(Gtemp);

            end
        end

    elseif (real(freq) == zero(real(freq)))

        ## @inbounds @sync @distributed for qq=1:mySingleMolData.numatoms
        @inbounds for qq=1:mySingleMolData.numatoms
            @inbounds for pp=qq:mySingleMolData.numatoms

                GFVACGG!(Gtemp, 1, 1, freq, blochk, view(mySingleMolData.atompos0, 1:3, pp),
                         view(mySingleMolData.atompos0, 1:3, qq),
                         myPeriodicData.latticevecs, myPeriodicData.reciprocalvecs,
                         mySingleMolData.alpha0[pp], mySingleMolData.alpha0[qq],
                         myPeriodicData.numdims);
                Gtemp .= mySingleMolData.Qe[pp] * mySingleMolData.Qe[qq] * Gtemp;
                mySingleMolData.GvacGGinf[3*pp-2:3*pp, 3*qq-2:3*qq] .= Gtemp;
                mySingleMolData.GvacGGinf[3*qq-2:3*qq, 3*pp-2:3*pp] .= adjoint(Gtemp);
            end
        end

    else

        ## @inbounds @sync @distributed for qq=1:mySingleMolData.numatoms
        @inbounds for qq=1:mySingleMolData.numatoms
            @inbounds for pp=1:mySingleMolData.numatoms

                GFVACGG!(Gtemp, 1, 1, freq, blochk,
                         view(mySingleMolData.atompos0, 1:3, pp),
                         view(mySingleMolData.atompos0, 1:3, qq),
                         myPeriodicData.latticevecs, myPeriodicData.reciprocalvecs,
                         mySingleMolData.alpha0[pp], mySingleMolData.alpha0[qq],
                         myPeriodicData.numdims);
                Gtemp .= mySingleMolData.Qe[pp] * mySingleMolData.Qe[qq] * Gtemp;
                mySingleMolData.GvacGGinf[3*pp-2:3*pp, 3*qq-2:3*qq] .= Gtemp;

            end
        end

    end

end

"""

    assemble1MolGvacinfdiag!(mySingleMolData::OneMol{FT},
                             freq::Union{FT, Complex{FT}}) where FT<:AbstractFloat

Assemble the scalar coincident Green's function interaction matrix elements of a single
molecule in vacuum, `mySingleMolData.GvacGGinfdiag`, in-place, at frequency `freq`.

!!! warning

    This function depends on `mySingleMolData.alpha0` having been initialized to the desired
    atomic susceptibilities at the same `freq` (and wavevector), but does not check this.

!!! note

    This function multiplies the diagonal Green's function scalar for every atom by the
    charge squared of each associated electronic oscillator.

"""
function assemble1MolGvacinfdiag!(mySingleMolData::OneMol{FT},
                                  freq::Union{FT, Complex{FT}}) where FT<:AbstractFloat

    mySingleMolData.GvacGGinfdiag .= map(u -> GFVACGGcoincident(freq, u, u),
                                         mySingleMolData.alpha0) .* mySingleMolData.Qe.^2;

end

## TO DELETE: subtract1MolalphainvGvac!() appears to be deprecated, not used anywhere
"""

    subtract1MolalphainvGvac!(myMolData::OneMol{<:AbstractFloat})

Subtract the Green's function interaction matrix of a single molecule in vacuum,
`myMolData.GvacGGinf`, in-place, from its inverse susceptibility `myMolData.alphainv`.
Separately dispatch cases with versus without phonons.

!!! warning

This function depends on `myMolData.alphainv` having been initialized to the desired
molecular susceptibility at the same `freq` and `blochk`, but does not check this.

"""
function subtract1MolalphainvGvac!(myMolData::OneMol_NoPhonons{<:AbstractFloat})

    subtractfromdiag!(myMolData.GvacGGinf, vecNto3N(myMolData.alphainv));

end
function subtract1MolalphainvGvac!(myMolData::OneMol_WithPhonons{<:AbstractFloat})

    myMolData.GvacGGinf[:, :] .= myMolData.alphainv[:, :] .- myMolData.GvacGGinf[:, :];

end

"""

    assembleGscaMolBlock!(A::SharedArray{Complex{FT}, 2},
                          myAllMolData::MolSystem{<:OneMol{FT}},
                          myTransData::TransData{FT}, nn1::Integer, nn2::Integer,
                          startidx1::Integer, startidx2::Integer,
                          myPeriodicData::PeriodicData{FT},
                          GFSCAGG!::Function, freq::Union{FT, Complex{FT}},
                          blochk::AbstractArray{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

Assemble the scattering Green's function interaction matrix (evaluated via the function
`GFSCAGG!`) between molecules `nn1` and `nn2` (sources evaluated in the latter, fields
evaluated in the former), at frequency `freq` and wavevector `blochk`, using atomic
polarizability and coordinate data respectively from `myAllMolData` and `myTransData`, and
add the result in-place to the matrix block
`A[startidx1:endidx1, startidx2:endidx2]` where `endidx1 - startidx1 + 1` is 3 times the
number of atoms in molecule `nn1`, and likewise `endidx2 - startidx2 + 1` is 3 times the
number of atoms in molecule `nn2`. If `nn1` and `nn2` are the same, the matrix
block is symmetric if the periodic dimensionality is 0, or Hermitian if the periodic
dimensionality is nonzero but `freq` is purely imaginary, so only the lower triangle of
blocks is assembled. Otherwise, the full set of matrix blocks is assembled.

!!! note

    If `GFSCAGG!` is the scattering Green's function of vacuum (returning a 3-by-3 block of
    zeros), namely `GFVACGG_sca!` or `GFVACGG0_sca!`, this function does nothing, saving on
    computations.

!!! warning

    Currently, only `GFPECGG!` and `GFPECGG0!` are supported for `GFSCAGG!`, but this is not
    explicitly checked, and may change in future versions of this code.

!!! warning

    This function depends on `mySingleMolData.alpha0` having been initialized to the desired
    atomic susceptibilities at the same `freq` and `blochk`, as well as
    `myTransData.allAtomPos` having been modified to reflect the desired transformation, but
    does not check these.

!!! note

    This function multiplies the Green's function tensor for every pair of atoms by the
    charges of each associated electronic oscillator.

"""
function assembleGscaMolBlock!(A::AbstractArray{Complex{FT}, 2},
                               myAllMolData::MolSystem{<:OneMol{FT}},
                               myTransData::TransData{FT}, nn1::Integer, nn2::Integer,
                               startidx1::Integer, startidx2::Integer,
                               myPeriodicData::PeriodicData{FT},
                               GFSCAGG!::Function, freq::Union{FT, Complex{FT}},
                               blochk::AbstractArray{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

    if ((GFSCAGG!) != (GFVACGG_sca!) && (GFSCAGG!) != (GFVACGG0_sca!)) # only do anything if we aren't in vacuum

        currmolind1 = myAllMolData.dupemolinds[nn1];
        currmolind2 = myAllMolData.dupemolinds[nn2];
        if (currmolind1 < 0)
            currmolind1 = myAllMolData.uniquemolinds[nn1];
        end
        if (currmolind2 < 0)
            currmolind2 = myAllMolData.uniquemolinds[nn2];
        end

        if (size(A, 1) < 3*myAllMolData.numatomslist[nn1])
            error("A (currently with ", size(A, 1), " rows) should have at least ",
                  3*myAllMolData.numatomslist[nn1], " rows in assembleGenvMolBlock!()");
        end
        if (size(A, 2) < 3*myAllMolData.numatomslist[nn2])
            error("A (currently with ", size(A, 2), " columns) should have at least ",
                  3*myAllMolData.numatomslist[nn2], " columns in assembleGenvMolBlock!()");
        end
        if (startidx1 < 1 || startidx1 > size(A, 1) - 3*myAllMolData.numatomslist[nn1] + 1)
            error("assembleGenvMolBlock!(): startidx1 (currently ", startidx1,
                  ") must be between 1 and ",
                  size(A, 1) - 3*myAllMolData.numatomslist[nn1] + 1,
                  " as A has ", size(A, 1), " rows and the block will have ",
                  3*myAllMolData.numatomslist[nn1], " rows");
        end
        if (startidx2 < 1 || startidx2 > size(A, 2) - 3*myAllMolData.numatomslist[nn2] + 1)
            error("assembleGenvMolBlock!(): startidx2 (currently ", startidx2,
                  ") must be between 1 and ",
                  size(A, 2) - 3*myAllMolData.numatomslist[nn2] + 1,
                  " as A has ", size(A, 2), " columns and the block will have ",
                  3*myAllMolData.numatomslist[nn2], " columns");
        end

        Gtemp = Array{eltype(A), 2}(undef, 3, 3);
        fill!(Gtemp, zero(eltype(Gtemp)));
        allatomoffset1 = 0;
        if (nn1 > 1)
            allatomoffset1 = myAllMolData.cumnumatomslist[nn1-1];
        end
        allatomoffset2 = 0;
        if (nn2 > 1)
            allatomoffset2 = myAllMolData.cumnumatomslist[nn2-1];
        end
        posfield = zeros(eltype(myTransData.allAtomPos), 3);
        possource = zeros(eltype(myTransData.allAtomPos), 3);

        if (nn1 == nn2 && myPeriodicData.numdims == 0)

            @inbounds for qq=1:myAllMolData.numatomslist[currmolind2]
                possource .= myTransData.allAtomPos[1:3, allatomoffset2 + qq];
                posfield .= myTransData.allAtomPos[1:3, allatomoffset1 + qq];
                GFSCAGG!(Gtemp, 1, 1, freq, posfield, possource,
                         myAllMolData.uniqueMolDataArray[currmolind1].alpha0[qq],
                         myAllMolData.uniqueMolDataArray[currmolind2].alpha0[qq]);
                sympart!(Gtemp);
                Gtemp .= myAllMolData.uniqueMolDataArray[currmolind1].Qe[qq] *
                         myAllMolData.uniqueMolDataArray[currmolind2].Qe[qq] *
                         Gtemp;
                A[(startidx1-1).+(3*qq-2:3*qq), (startidx2-1).+(3*qq-2:3*qq)] .+=
                    Gtemp;
            end

            @inbounds for qq=1:myAllMolData.numatomslist[currmolind2]
                possource .= myTransData.allAtomPos[1:3, allatomoffset2 + qq];
                @inbounds for pp=(qq+1):myAllMolData.numatomslist[currmolind1]
                    posfield .= myTransData.allAtomPos[1:3, allatomoffset1 + pp];
                    GFSCAGG!(Gtemp, 1, 1, freq, posfield, possource,
                             myAllMolData.uniqueMolDataArray[currmolind1].alpha0[pp],
                             myAllMolData.uniqueMolDataArray[currmolind2].alpha0[qq]);
                    Gtemp .= myAllMolData.uniqueMolDataArray[currmolind1].Qe[pp] *
                             myAllMolData.uniqueMolDataArray[currmolind2].Qe[qq] *
                             Gtemp;
                    A[(startidx1-1).+(3*pp-2:3*pp), (startidx2-1).+(3*qq-2:3*qq)] .+=
                    Gtemp;
                    A[(startidx2-1).+(3*qq-2:3*qq), (startidx1-1).+(3*pp-2:3*pp)] .+=
                    transpose(Gtemp);
                end
            end

        elseif (nn1 == nn2 && real(freq) == zero(real(freq)))

            @inbounds for qq=1:myAllMolData.numatomslist[currmolind2]
                possource .= myTransData.allAtomPos[1:3, allatomoffset2 + qq];
                posfield .= myTransData.allAtomPos[1:3, allatomoffset1 + qq];
                GFSCAGG!(Gtemp, 1, 1, freq, blochk, posfield, possource,
                         myPeriodicData.latticevecs, myPeriodicData.reciprocalvecs,
                         myAllMolData.uniqueMolDataArray[currmolind1].alpha0[qq],
                         myAllMolData.uniqueMolDataArray[currmolind2].alpha0[qq],
                         myPeriodicData.numdims);
                hermpart!(Gtemp);
                Gtemp .= myAllMolData.uniqueMolDataArray[currmolind1].Qe[qq] *
                         myAllMolData.uniqueMolDataArray[currmolind2].Qe[qq] *
                         Gtemp;
                A[(startidx1-1).+(3*qq-2:3*qq), (startidx2-1).+(3*qq-2:3*qq)] .+=
                    Gtemp;
            end

            @inbounds for qq=1:myAllMolData.numatomslist[currmolind2]
                possource .= myTransData.allAtomPos[1:3, allatomoffset2 + qq];
                @inbounds for pp=(qq+1):myAllMolData.numatomslist[currmolind1]
                    posfield .= myTransData.allAtomPos[1:3, allatomoffset1 + pp];
                    GFSCAGG!(Gtemp, 1, 1, freq, blochk, posfield, possource,
                             myPeriodicData.latticevecs, myPeriodicData.reciprocalvecs,
                             myAllMolData.uniqueMolDataArray[currmolind1].alpha0[pp],
                             myAllMolData.uniqueMolDataArray[currmolind2].alpha0[qq],
                             myPeriodicData.numdims);
                    Gtemp .= myAllMolData.uniqueMolDataArray[currmolind1].Qe[pp] *
                             myAllMolData.uniqueMolDataArray[currmolind2].Qe[qq] *
                             Gtemp;
                    A[(startidx1-1).+(3*pp-2:3*pp), (startidx2-1).+(3*qq-2:3*qq)] .+=
                    Gtemp;
                    A[(startidx2-1).+(3*qq-2:3*qq), (startidx1-1).+(3*pp-2:3*pp)] .+=
                    adjoint(Gtemp);
                end
            end

        else

            @inbounds for qq=1:myAllMolData.numatomslist[currmolind2]
                possource .= myTransData.allAtomPos[1:3, allatomoffset2 + qq];
                @inbounds for pp=1:myAllMolData.numatomslist[currmolind1]
                    posfield .= myTransData.allAtomPos[1:3, allatomoffset1 + pp];
                    GFSCAGG!(Gtemp, 1, 1, freq, blochk, posfield, possource,
                             myPeriodicData.latticevecs, myPeriodicData.reciprocalvecs,
                             myAllMolData.uniqueMolDataArray[currmolind1].alpha0[pp],
                             myAllMolData.uniqueMolDataArray[currmolind2].alpha0[qq],
                             myPeriodicData.numdims);
                    Gtemp .= myAllMolData.uniqueMolDataArray[currmolind1].Qe[pp] *
                             myAllMolData.uniqueMolDataArray[currmolind2].Qe[qq] *
                             Gtemp;
                    A[(startidx1-1).+(3*pp-2:3*pp), (startidx2-1).+(3*qq-2:3*qq)] .+=
                    Gtemp;
                end
            end

        end

    end

end

"""

    assembleGvacMolBlock!(A::AbstractArray{Complex{FT}, 2},
                          myAllMolData::MolSystem{<:OneMol{FT}},
                          myTransData::TransData{FT}, nn1::Integer, nn2::Integer,
                          startidx1::Integer, startidx2::Integer,
                          myPeriodicData::PeriodicData{FT},
                          freq::Union{FT, Complex{FT}},
                          blochk::AbstractArray{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

Assemble the vacuum Green's function interaction matrix (evaluated via the function
`GFVACGG!`) between molecules `nn1` and `nn2` (sources evaluated in the latter, fields
evaluated in the former), at frequency `freq` and wavevector `blochk`, using atomic
polarizability and coordinate data respectively from `myAllMolData` and `myTransData`, and
add the result in-place to the matrix block
`A[startidx1:endidx1, startidx2:endidx2]` where `endidx1 - startidx1 + 1` is 3 times the
number of atoms in molecule `nn1`, and likewise `endidx2 - startidx2 + 1` is 3 times the
number of atoms in molecule `nn2`. If `nn1` and `nn2` are the same, the matrix
block is symmetric if the periodic dimensionality is 0, or Hermitian if the periodic
dimensionality is nonzero but `freq` is purely imaginary, so only the lower triangle of
blocks is assembled. Otherwise, the full set of matrix blocks is assembled.

!!! tip

    Use `assemble1MolGvacinf!` when `nn1` and `nn2` are the same, then stamp that into `A`
    as needed. Only use this when `nn1` and `nn2` are different.

!!! warning

    This function depends on `mySingleMolData.alpha0` having been initialized to the desired
    atomic susceptibilities at the same `freq` and `blochk`, as well as
    `myTransData.allAtomPos` having been modified to reflect the desired transformation, but
    does not check these.

!!! note

    This function multiplies the Green's function tensor for every pair of atoms by the
    charges of each associated electronic oscillator.

"""
function assembleGvacMolBlock!(A::AbstractArray{Complex{FT}, 2},
                               myAllMolData::MolSystem{<:OneMol{FT}},
                               myTransData::TransData{FT},
                               nn1::Integer, nn2::Integer,
                               startidx1::Integer,
                               startidx2::Integer,
                               myPeriodicData::PeriodicData{FT},
                               freq::Union{FT, Complex{FT}},
                               blochk::AbstractArray{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

    currmolind1 = myAllMolData.dupemolinds[nn1];
    if (currmolind1 < 0)
        currmolind1 = myAllMolData.uniquemolinds[nn1];
    end
    currmolind2 = myAllMolData.dupemolinds[nn2];
    if (currmolind2 < 0)
        currmolind2 = myAllMolData.uniquemolinds[nn2];
    end

    if (size(A, 1) < 3*myAllMolData.numatomslist[nn1])
        error("A (currently with ", size(A, 1), " rows) should have at least ",
              3*myAllMolData.numatomslist[nn1], " rows in assembleGenvMolBlock!()");
    end
    if (size(A, 2) < 3*myAllMolData.numatomslist[nn2])
        error("A (currently with ", size(A, 2), " rows) should have at least ",
              3*myAllMolData.numatomslist[nn2], " rows in assembleGenvMolBlock!()");
    end
    if (startidx1 < 1 || startidx1 > size(A, 1) - 3*myAllMolData.numatomslist[nn1] + 1)
        error("assembleGenvMolBlock!(): startidx1 (currently ", startidx1,
              ") must be between 1 and ",
              size(A, 1) - 3*myAllMolData.numatomslist[nn1] + 1,
              " as A has ", size(A, 1), " rows and the block will have ",
              3*myAllMolData.numatomslist[nn1], " rows");
    end
    if (startidx2 < 1 || startidx2 > size(A, 2) - 3*myAllMolData.numatomslist[nn2] + 1)
        error("assembleGenvMolBlock!(): startidx2 (currently ", startidx2,
              ") must be between 1 and ",
              size(A, 2) - 3*myAllMolData.numatomslist[nn2] + 1,
              " as A has ", size(A, 2), " columns and the block will have ",
              3*myAllMolData.numatomslist[nn2], " columns");
    end

    Gtemp = Array{eltype(A), 2}(undef, 3, 3);
    fill!(Gtemp, zero(eltype(Gtemp)));
    allatomoffset1 = 0;
    if (nn1 > 1)
        allatomoffset1 = myAllMolData.cumnumatomslist[nn1-1];
    end
    allatomoffset2 = 0;
    if (nn2 > 1)
        allatomoffset2 = myAllMolData.cumnumatomslist[nn2-1];
    end
    posfield = zeros(eltype(myTransData.allAtomPos), 3);
    possource = zeros(eltype(myTransData.allAtomPos), 3);

    if (nn1 == nn2 && myPeriodicData.numdims == 0)

        @inbounds for qq=1:myAllMolData.numatomslist[currmolind2]
            possource .= myTransData.allAtomPos[1:3, allatomoffset2 + qq];
            @inbounds for pp=(qq+1):myAllMolData.numatomslist[currmolind1]
                posfield .= myTransData.allAtomPos[1:3, allatomoffset1 + pp];
                GFVACGG!(Gtemp, 1, 1, freq, posfield, possource,
                         myAllMolData.uniqueMolDataArray[currmolind1].alpha0[pp],
                         myAllMolData.uniqueMolDataArray[currmolind2].alpha0[qq]);
                Gtemp .= myAllMolData.uniqueMolDataArray[currmolind1].Qe[pp] *
                         myAllMolData.uniqueMolDataArray[currmolind2].Qe[qq] *
                         Gtemp;
                A[(startidx1-1).+(3*pp-2:3*pp), (startidx2-1).+(3*qq-2:3*qq)] .+=
                Gtemp;
                A[(startidx2-1).+(3*qq-2:3*qq), (startidx1-1).+(3*pp-2:3*pp)] .+=
                transpose(Gtemp);
            end
        end

    elseif (nn1 == nn2 && real(freq) == zero(real(freq)))

        @inbounds for qq=1:myAllMolData.numatomslist[currmolind2]
            possource .= myTransData.allAtomPos[1:3, allatomoffset2 + qq];
            @inbounds for pp=(qq+1):myAllMolData.numatomslist[currmolind1]
                posfield .= myTransData.allAtomPos[1:3, allatomoffset1 + pp];
                GFVACGG!(Gtemp, 1, 1, freq, blochk, posfield, possource,
                         myPeriodicData.latticevecs, myPeriodicData.reciprocalvecs,
                         myAllMolData.uniqueMolDataArray[currmolind1].alpha0[pp],
                         myAllMolData.uniqueMolDataArray[currmolind2].alpha0[qq],
                         myPeriodicData.numdims);
                Gtemp .= myAllMolData.uniqueMolDataArray[currmolind1].Qe[pp] *
                         myAllMolData.uniqueMolDataArray[currmolind2].Qe[qq] *
                         Gtemp;
                A[(startidx1-1).+(3*pp-2:3*pp), (startidx2-1).+(3*qq-2:3*qq)] .+=
                Gtemp;
                A[(startidx2-1).+(3*qq-2:3*qq), (startidx1-1).+(3*pp-2:3*pp)] .+=
                adjoint(Gtemp);
            end
        end

    else

        @inbounds for qq=1:myAllMolData.numatomslist[currmolind2]
            possource .= myTransData.allAtomPos[1:3, allatomoffset2 + qq];
            @inbounds for pp=1:myAllMolData.numatomslist[currmolind1]
                posfield .= myTransData.allAtomPos[1:3, allatomoffset1 + pp];
                GFVACGG!(Gtemp, 1, 1, freq, blochk, posfield, possource,
                         myPeriodicData.latticevecs, myPeriodicData.reciprocalvecs,
                         myAllMolData.uniqueMolDataArray[currmolind1].alpha0[pp],
                         myAllMolData.uniqueMolDataArray[currmolind2].alpha0[qq],
                         myPeriodicData.numdims);
                Gtemp .= myAllMolData.uniqueMolDataArray[currmolind1].Qe[pp] *
                         myAllMolData.uniqueMolDataArray[currmolind2].Qe[qq] *
                         Gtemp;
                A[(startidx1-1).+(3*pp-2:3*pp), (startidx2-1).+(3*qq-2:3*qq)] .+=
                Gtemp;
            end
        end

    end

end


text:
module FEM

using DocStringExtensions
using StaticArrays
using LinearAlgebra
using SparseArrays

using WaveProp
using WaveProp.Mesh
using WaveProp.Geometry
using WaveProp.Integration
using WaveProp.Interpolation

include("elementary_matrix.jl")
include("dof_numbering.jl")

export elementary_matrix
export DOFNumbering, LocalDOFNumbering, dofs, elt2dof, local_dof_numbering
export assembly

end

ext: jl

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
# convert a Reaction to a ChemicalHyperEdge
function Base.convert(::Type{ChemicalHyperEdge}, r::Catalyst.Reaction)
    t_subs = eltype(r.substrates)
    t_prods = eltype(r.products)
    t = (isequal(t_subs, t_prods) ? t_subs : Simulacrum.get_lowest_subtype(t_subs, t_prods))
    ChemicalHyperEdge(SpeciesSet(Vector{t}(r.substrates), r.substoich), SpeciesSet(Vector{t}(r.products), r.prodstoich), r.rate)
end
function Base.convert(::Type{ChemicalHyperEdge}, rs::AbstractVector{T}) where {T<:Catalyst.Reaction}
    convert.(ChemicalHyperEdge, rs)
end

# now these functions may be defined via the conversion defined above
isunimolecular(r::Catalyst.Reaction) = isunimolecular(convert(ChemicalHyperEdge, r))
ismultimolecular(r::Catalyst.Reaction) = !(isunimolecular(r))

# convert a ChemicalHyperEdge to a Reaction
# note that rate(che)[1] is needed -- not ideal
function Base.convert(::Type{Catalyst.Reaction}, che::ChemicalHyperEdge)
    Reaction(rate(che)[1], src(che), tgt(che), src_stoich(che), tgt_stoich(che))
end
function Base.convert(::Type{Catalyst.Reaction}, ches::AbstractVector{T}) where {T<:ChemicalHyperEdge}
    convert.(Catalyst.Reaction, ches)
end

# convert a ReactionSystem to a ChemicalHyperGraph
function Base.convert(::Type{ChemicalHyperGraph}, rsys::Catalyst.ReactionSystem)
    ChemicalHyperGraph(convert(ChemicalHyperEdge, equations(rsys)))
end

# return reactions from a chemical hypergraph
catalyst_reactions(chx::ChemicalHyperGraph) = convert(Catalyst.Reaction, hyperedges(chx))

# convert a ChemicalHyperGraph to a ReactionSystem
function Base.convert(::Type{Catalyst.ReactionSystem}, chx::ChemicalHyperGraph; p = [], iv = Symbolics.variable(:t))
    Catalyst.ReactionSystem(catalyst_reactions(chx), iv, vertices(chx), p, name = :rsys)
end

# overload the ChemicalHyperGraph constructor to Reaction and ReactionSystem types
ChemicalHyperGraph(r::Catalyst.Reaction) = ChemicalHyperGraph(convert(ChemicalHyperEdge, r))
function ChemicalHyperGraph(rs::AbstractVector{T}) where {T<:Catalyst.Reaction}
    ChemicalHyperGraph(convert(ChemicalHyperEdge, rs))
end
ChemicalHyperGraph(rsys::Catalyst.ReactionSystem) = convert(ChemicalHyperGraph, rsys)


text:
# This file is part of Kpax3. License is MIT.

function biased_random_walk!(data::Matrix{UInt8},
                             priorR::PriorRowPartition,
                             priorC::PriorColPartition,
                             settings::KSettings,
                             support::MCMCSupport,
                             state::AminoAcidState)
  k = state.k

  i = StatsBase.sample(1:support.n)
  hi = state.R[i]
  hj = 0

  v = 0
  if state.v[hi] > 1
    if k > 1
      v = StatsBase.sample(1:k)

      if v == k
        # move i to a new cluster
        k += 1
      end
    else
      # move i to a new cluster
      k = 2
    end
  else
    # move i to another cluster
    k -= 1
    v = StatsBase.sample(1:k)
  end

  if k == state.k
    hj = state.cl[v] < hi ? state.cl[v] : state.cl[v + 1]

    initsupportbrwmove!(i, hi, hj, data, support, state)

    support.lograR = logratiopriorrowmove(state.v[hi], state.v[hj], priorR)

    updatelogmargliki!(priorC, support)
    updatelogmarglikj!(priorC, support)

    logmarglikbrwmove!(state.cl, state.k, hi, hj, priorC, support)

    ratio = exp(support.lograR + support.logmlikcandidate - support.logmlik)

    if ratio >= 1 || ((ratio > 0) && (rand() <= ratio))
      performbrwmove!(i, hi, hj, support, state)
    end
  elseif k > state.k
    initsupportbrwsplit!(k, i, hi, data, priorC, settings, support, state)

    support.lograR = logratiopriorrowsplit(k, state.v[hi], priorR)

    updatelogmargliki!(priorC, support)
    updatelogmarglikj!(priorC, support)

    logmarglikbrwsplit!(state.cl, state.k, hi, priorC, support)

    ratio = exp(support.lograR + support.logmlikcandidate - support.logmlik)

    if ratio >= 1 || ((ratio > 0) && (rand() <= ratio))
      performbrwsplit!(i, hi, settings, support, state)
    end
  else
    hj = state.cl[v] < hi ? state.cl[v] : state.cl[v + 1]

    initsupportbrwmerge!(i, hj, data, support, state)

    support.lograR = logratiopriorrowmerge(k, state.v[hj], priorR)

    updatelogmarglikj!(priorC, support)

    logmarglikbrwmerge!(state.cl, state.k, hi, hj, priorC, support)

    ratio = exp(support.lograR + support.logmlikcandidate - support.logmlik)

    if ratio >= 1 || ((ratio > 0) && (rand() <= ratio))
      performbrwmerge!(i, hi, hj, settings, support, state)
    end
  end

  nothing
end

function initsupportbrwmerge!(i::Int,
                              hj::Int,
                              data::Matrix{UInt8},
                              support::MCMCSupport,
                              state::AminoAcidState)
  support.vj = state.v[hj] + 1

  for b in 1:support.m
    support.nj[b] = state.n1s[hj, b] + float(data[b, i])
  end

  nothing
end

function initsupportbrwmove!(i::Int,
                             hi::Int,
                             hj::Int,
                             data::Matrix{UInt8},
                             support::MCMCSupport,
                             state::AminoAcidState)
  support.vi = state.v[hi] - 1
  support.vj = state.v[hj] + 1

  for b in 1:support.m
    support.ni[b] = state.n1s[hi, b] - float(data[b, i])
    support.nj[b] = state.n1s[hj, b] + float(data[b, i])
  end

  idx = 0
  for a in 1:state.v[hi]
    if state.unit[hi][a] != i
      idx += 1
      support.ui[idx] = state.unit[hi][a]
    end
  end

  nothing
end

function initsupportbrwsplit!(k::Int,
                              i::Int,
                              hi::Int,
                              data::Matrix{UInt8},
                              priorC::AminoAcidPriorCol,
                              settings::KSettings,
                              support::MCMCSupport,
                              state::AminoAcidState)
  resizesupport!(support, k)

  support.vi = state.v[hi] - 1
  support.vj = 1

  for b in 1:support.m
    support.nj[b] = float(data[b, i])
    support.ni[b] = state.n1s[hi, b] - support.nj[b]
  end

  idx = 0
  for a in 1:state.v[hi]
    if state.unit[hi][a] != i
      idx += 1
      support.ui[idx] = state.unit[hi][a]
    end
  end

  nothing
end

function performbrwmove!(i::Int,
                         hi::Int,
                         hj::Int,
                         support::MCMCSupport,
                         state::AminoAcidState)
  state.R[i] = hj

  state.v[hi] = support.vi
  state.v[hj] = support.vj

  for b in 1:support.m
    state.n1s[hi, b] = support.ni[b]
    state.n1s[hj, b] = support.nj[b]

    support.lp[1, hi, b] = support.lpi[1, b]
    support.lp[2, hi, b] = support.lpi[2, b]
    support.lp[3, hi, b] = support.lpi[3, b]
    support.lp[4, hi, b] = support.lpi[4, b]

    support.lp[1, hj, b] = support.lpj[1, b]
    support.lp[2, hj, b] = support.lpj[2, b]
    support.lp[3, hj, b] = support.lpj[3, b]
    support.lp[4, hj, b] = support.lpj[4, b]
  end

  copyto!(state.unit[hi], 1, support.ui, 1, support.vi)

  resize!(state.unit[hj], state.v[hj])
  state.unit[hj][state.v[hj]] = i

  state.logpR += support.lograR

  support.logmlik = support.logmlikcandidate

  nothing
end

function performbrwsplit!(i::Int,
                          hi::Int,
                          settings::KSettings,
                          support::MCMCSupport,
                          state::AminoAcidState)
  resizestate!(state, state.k + 1, settings)

  hj = findfirst(state.emptycluster)

  state.R[i] = hj

  state.emptycluster[hj] = false

  support.k = 0
  for a in 1:length(state.emptycluster)
    if !state.emptycluster[a]
      support.k += 1
      state.cl[support.k] = a
    end
  end

  state.k = support.k

  state.v[hi] = support.vi
  state.v[hj] = support.vj

  for b in 1:support.m
    state.n1s[hi, b] = support.ni[b]
    state.n1s[hj, b] = support.nj[b]

    support.lp[1, hi, b] = support.lpi[1, b]
    support.lp[2, hi, b] = support.lpi[2, b]
    support.lp[3, hi, b] = support.lpi[3, b]
    support.lp[4, hi, b] = support.lpi[4, b]

    support.lp[1, hj, b] = support.lpj[1, b]
    support.lp[2, hj, b] = support.lpj[2, b]
    support.lp[3, hj, b] = support.lpj[3, b]
    support.lp[4, hj, b] = support.lpj[4, b]
  end

  copyto!(state.unit[hi], 1, support.ui, 1, support.vi)

  resize!(state.unit[hj], 1)
  state.unit[hj][1] = i

  state.logpR += support.lograR

  support.logmlik = support.logmlikcandidate

  nothing
end

function performbrwmerge!(i::Int,
                          hi::Int,
                          hj::Int,
                          settings::KSettings,
                          support::MCMCSupport,
                          state::AminoAcidState)
  state.R[i] = hj

  state.emptycluster[hi] = true

  support.k = 0
  for a in 1:length(state.emptycluster)
    if !state.emptycluster[a]
      support.k += 1
      state.cl[support.k] = a
    end
  end

  state.k = support.k

  state.v[hj] = support.vj

  for b in 1:support.m
    state.n1s[hj, b] = support.nj[b]

    support.lp[1, hj, b] = support.lpj[1, b]
    support.lp[2, hj, b] = support.lpj[2, b]
    support.lp[3, hj, b] = support.lpj[3, b]
    support.lp[4, hj, b] = support.lpj[4, b]
  end

  resize!(state.unit[hj], state.v[hj])
  state.unit[hj][state.v[hj]] = i

  state.logpR += support.lograR

  support.logmlik = support.logmlikcandidate

  nothing
end


text:
#!/usr/bin/env julia

using Colors
using Compose

function golden_rect(n::Int)
    if n == 0 return context() end
    c = compose(context(), rectangle(), fill(LCHab(90, 80, 70 - 15n)))
    compose(c, context(units=UnitBox(0, -1/φ, 1h, 1/φ), rotation=Rotation(π/2, 0, 1)),  golden_rect(n - 1))
end

draw(SVG("golden_rect.svg", φ * 3inch, 3inch),
    compose(golden_rect(10), fill(nothing), stroke(colorant"white"), linewidth(0.2mm)))


text:
using MPI, Test
include(joinpath("..", "..", "testhelpers.jl"))

@testset "DGmethods" begin
    runmpi(joinpath(@__DIR__, "vars_test.jl"))
    runmpi(joinpath(@__DIR__, "integral_test.jl"))
    runmpi(joinpath(@__DIR__, "integral_test_sphere.jl"))
    runmpi(joinpath(@__DIR__, "Euler/isentropicvortex.jl"))
    runmpi(joinpath(@__DIR__, "Euler/isentropicvortex_imex.jl"))
    runmpi(joinpath(@__DIR__, "Euler/isentropicvortex_multirate.jl"))
    runmpi(joinpath(
        @__DIR__,
        "advection_diffusion/pseudo1D_advection_diffusion.jl",
    ))
    runmpi(joinpath(
        @__DIR__,
        "advection_diffusion/advection_diffusion_model_1dimex_bgmres.jl",
    ))
    runmpi(joinpath(@__DIR__, "compressible_Navier_Stokes/ref_state.jl"))
    runmpi(joinpath(@__DIR__, "horizontal_integral_test.jl"))
    runmpi(joinpath(@__DIR__, "courant.jl"), ntasks = 2)
end


text:
#=
Definitions of the dfe (distribution of fitness effects) functions which define
   selection coefficients in terms of Gamma distributions.
Split off from nn_poplist.jl  on 5/14/18.
=#
using Distributions

function dfe_deleterious( x::Int64; alpha::Float64=0.2, beta::Float64=0.5 )
  global dist_deleterious
  if !isdefined(Main,:dist_deleterious)
    dist_deleterious = Distributions.Gamma(alpha,beta)
  end
  return max(0.1,1.0-rand(dist_deleterious))
end

function dfe_advantageous( x::Int64; alpha::Float64=1.0, beta::Float64=0.5 )
  global dist_advantageous
  if !isdefined(Main,:dist_advantageous)
    dist_advantageous = Distributions.Gamma(alpha,beta)
  end
  return 1.0+rand(dist_advantageous)
end

@doc """ function dfe_mixed( x::Int64; adv_probability::Float64=0.2, alpha_disadv::Float64=0.2, alpha_adv::Float64=1.0, beta_disadv::Float64=1.0, beta_adv::Float64=0.01 )
Mixed advantageous and deleterious.
"""

function dfe_mixed( x::Int64; adv_probability::Float64=0.2, alpha_disadv::Float64=0.2, alpha_adv::Float64=1.0, beta_disadv::Float64=1.0, beta_adv::Float64=0.05 )
  global dist_deleterious
  if !isdefined(Main,:dist_deleterious)
    dist_deleterious = Distributions.Gamma(alpha_disadv,beta_disadv)
  end
  global dist_advantageous
  if !isdefined(Main,:dist_advantageous)
    dist_advantageous = Distributions.Gamma(alpha_adv,beta_adv)
  end
  rand() < adv_probability ? 1.0+rand(dist_advantageous) : max(0.01,1.0-rand(dist_deleterious))
end

@doc """ function dfe_mod( )
Increased fitness for every x which is divisible by modulaus.
The increase in fitness is by a multiple of fit_inc.
"""
function dfe_mod( x::Int64; modulus::Int64=5, fit_inc::Float64=1.2 )
  if x % modulus == 0
	  return fit_inc
  else
    return 1.0  # flat fitness for no selection
  end
end

@doc """function dfe_neutral( x::Int64 )
"""
function dfe_neutral( x::Int64 )
  return 1.0   # flat fitness for no selection
end

@doc """function dfe_fixed( )
"""
function dfe_fixed( x::Int64; s::Float64=0.0 )
  return 1.0 + s   # fitness corresponding to selection coefficient s
end

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
# ---------------------------------------------------------------------------
# empirical CDF

@recipe function f(ecdf::StatsBase.ECDF)
    seriestype := :steppost
    legend --> :topleft
    x = [ecdf.sorted_values[1]; ecdf.sorted_values]
    if :weights in propertynames(ecdf) && !isempty(ecdf.weights)
         # support StatsBase versions >v0.32.0
        y = [0; cumsum(ecdf.weights) ./ sum(ecdf.weights)]
    else
        y = range(0, 1; length = length(x))
    end
    x, y
end

@userplot ECDFPlot
recipetype(::Val{:ecdfplot}, args...) = ECDFPlot(args)
@recipe function f(p::ECDFPlot)
    x = p.args[1]
    if !isa(x, StatsBase.ECDF)
        x = StatsBase.ecdf(x)
    end
    x
end


text:
"""
DFS
---
    compute depth first search distances and returns the distance (d), the discover (dt),
    the finish time(ft), and the predecessor array (pred) in the tuple (d,dt,ft, pred).\n
    pred[i] = 0 if vertex i is in a component not reachable from u and i != u.

Functions
---------
- (d,dt,ft,pred) = dfs(A::MatrixNetwork,u::Int64,full::Int64,target::Int64)
- (d,dt,ft,pred) = dfs(A::MatrixNetwork,u::Int64)
- (d,dt,ft,pred) = dfs{T}(A::SparseMatrixCSC{T,Int64},u::Int64,full::Int64,target::Int64)
- (d,dt,ft,pred) = dfs{T}(A::SparseMatrixCSC{T,Int64},u::Int64)
- (d,dt,ft,pred) = dfs(ei::Vector{Int64},ej::Vector{Int64},u::Int64,full::Int64,target::Int64)
- (d,dt,ft,pred) = dfs(ei::Vector{Int64},ej::Vector{Int64},u::Int64)

Example
-------
~~~
A = load_matrix_network("dfs_example")
(d,dt,ft,pred)  = dfs(A,1)
~~~
"""

function dfs(A::MatrixNetwork,u::Int64,full::Int64,target::Int64)

    (rp,ci) = (A.rp,A.ci)
    n=length(rp)-1
    d=-1*ones(Int64,n)
    dt=-1*ones(Int64,n)
    ft=-1*ones(Int64,n)
    pred=zeros(Int64,n)
    rs=zeros(Int64,2*n)
    rss=0 # recursion stack holds two nums (v,ri)

    #start dfs at u
    t=0
    targethit=0
    for i=1:n
        if i==1
            v=u
        else
            v=mod(u+i-1,n)+1
            if d[v]>0
                continue
            end
        end
        d[v]=0
        dt[v]=t
        t=t+1
        ri=rp[v]
        rss=rss+1
        rs[2*rss-1]=v
        rs[2*rss]=ri # add v to the stack

        while rss>0
            v=rs[2*rss-1]
            ri=rs[2*rss]
            rss=rss-1 # pop v from the stack
            if v==target || targethit == 1
                ri=rp[v+1]
                targethit=1 # end the algorithm if v is the target
            end
            while ri<rp[v+1] #Traversing thru all its adjancencies
                w=ci[ri]
                ri=ri+1
                if d[w]<0
                    d[w]=d[v]+1
                    pred[w]=v
                    rss=rss+1
                    rs[2*rss-1]=v
                    rs[2*rss]=ri # add v to the stack
                    v=w
                    ri=rp[w]
                    dt[v]=t
                    t=t+1
                    continue #discover a new vertex!
                end
            end
            ft[v]=t
            t=t+1 # finish with v
        end
        if full == 0
            break
        end
    end
    return (d,dt,ft,pred)
end

############################
### Additional functions ###
############################

function dfs(A::MatrixNetwork,u::Int64)
    return dfs(A,u,0,0);
end

## CSC sparse matrices:
function dfs{T}(A::SparseMatrixCSC{T,Int64},u::Int64)
    return dfs(MatrixNetwork(A),u)
end

function dfs{T}(A::SparseMatrixCSC{T,Int64},u::Int64,full::Int64,target::Int64)
    return dfs(MatrixNetwork(A),u,full,target)
end

## Triplet Format:
function dfs(ei::Vector{Int64},ej::Vector{Int64},u::Int64)
    return dfs(MatrixNetwork(ei,ej),u)
end

function dfs(ei::Vector{Int64},ej::Vector{Int64},u::Int64,full::Int64,target::Int64)
    return dfs(MatrixNetwork(ei,ej),u,full,target)
end

## CSR sparse matrices:
function dfs{T}(rp::Vector{Int64},ci::Vector{Int64},vals::Vector{T},n::Int64,u::Int64)
    return dfs(MatrixNetwork(n,rp,ci,vals),u)
end

function dfs{T}(rp::Vector{Int64},ci::Vector{Int64},vals::Vector{T},n::Int64,u::Int64,full::Int64,target::Int64)
    return dfs(MatrixNetwork(n,rp,ci,vals),u,full,target)
end


text:
using FEASTSolver
using LinearAlgebra
using Profile
using NonlinearEigenproblems: PEP, compute_Mder, contour_block_SS, polyeig, compute_Mlincomb, compute_MM

function info(Λ, X, residuals, c, r)
    inside(x) = in_contour(x, c, r)
    in_eig = Λ[inside.(Λ)]
    in_res = residuals[inside.(Λ)]
    print("number inside : ")
    println(size(Λ[inside.(Λ)])[1])
    if sum(inside.(Λ)) > 0
        in_res_conv = in_res[in_res .<= 1e-5]
        in_eig_conv = in_eig[in_res .<= 1e-5]
        print("number inside converged : ")
        println(size(in_eig_conv)[1])
        print("max res inside: ")
        println(maximum(residuals[inside.(Λ)]))
        if size(in_res_conv, 1) > 0
            print("max res inside non spurious: ")
            println(maximum(in_res_conv))
        end
    end
    # display(in_eig_conv)
end

function bf_arrays()
    N = diagm(-1 => [1,1,1,1,1,1,1])
    Mh0 = 1/6*(4*I + N + N')
    Mh1 = N - N'
    Mh2 = -1*(2*I - N - N')
    Mh3 = Mh1
    Mh4 = -Mh2
    c = [0.6 1.3; 1.3 0.1; 0.1 1.2; 1.0 1.0; 1.2 1.0]
    I8 = Matrix(I, 8, 8)
    M0 = c[1,1] * kron(I8, Mh0) + c[1,2] * kron(Mh0, I8)
    M1 = c[2,1] * kron(I8, Mh1) + c[2,2] * kron(Mh1, I8)
    M2 = c[3,1] * kron(I8, Mh2) + c[3,2] * kron(Mh2, I8)
    M3 = c[4,1] * kron(I8, Mh3) + c[4,2] * kron(Mh3, I8)
    M4 = c[5,1] * kron(I8, Mh4) + c[5,2] * kron(Mh4, I8)
    return [M0, M1, M2, M3, M4]
end

function bf()
    N = diagm(-1 => [1,1,1,1,1,1,1])
    Mh0 = 1/6*(4*I + N + N')
    Mh1 = N - N'
    Mh2 = -1*(2*I - N - N')
    Mh3 = Mh1
    Mh4 = -Mh2
    c = [0.6 1.3; 1.3 0.1; 0.1 1.2; 1.0 1.0; 1.2 1.0]
    I8 = Matrix(I, 8, 8)
    M0 = c[1,1] * kron(I8, Mh0) + c[1,2] * kron(Mh0, I8)
    M1 = c[2,1] * kron(I8, Mh1) + c[2,2] * kron(Mh1, I8)
    M2 = c[3,1] * kron(I8, Mh2) + c[3,2] * kron(Mh2, I8)
    M3 = c[4,1] * kron(I8, Mh3) + c[4,2] * kron(Mh3, I8)
    M4 = c[5,1] * kron(I8, Mh4) + c[5,2] * kron(Mh4, I8)
    function butterfly(A0, A1, A2, A3, A4)
        f(z) =  z^4*A4 + z^3*A3 + z^2*A2 + z*A1 + A0
    end

    return butterfly(M0, M1, M2, M3, M4)
end

C = complex(1.0,1.0)
R = 0.5

T = bf()
e, v, res = nlfeast!(T, rand(ComplexF64,64,30), 2^2, 0, c=complex(1.0,1.0), r=0.5)
e, v, res = @timev nlfeast!(T, rand(ComplexF64,64,20), 2^9, 0, c=complex(1.0,1.0), r=0.5, debug=false, ϵ=10e-16, spurious=5e-3)
# e, v, res = nlfeast_moments!(T, rand(ComplexF64,64,15), 2^4, 40, moments=2, c=complex(1.0,1.0), r=0.5, debug=true, ϵ=10e-16, spurious=5e-3)
# # e, v, res = companion(bf_arrays())
# # e, v, res = @profile nlfeast!(T, rand(ComplexF64,64,30), 2^4, 20, c=complex(1.0,1.0), r=0.5, debug=true, ϵ=10e-16)
# # Profile.print()
# # e, v, res = nlfeast_it!(T, rand(ComplexF64,64,30), 2^3, 0, c=complex(1.0,1.0), r=0.5)
# # e, v, res = @timev nlfeast_it!(T, rand(ComplexF64,64,30), 2^3, 20, c=complex(1.0,1.0), r=0.5, debug=true, ϵ=10e-14)
# # display(e)
# inside(x) = in_contour(x, complex(1.0, 1.0), 0.5)
# in_eig = e[inside.(e)]
# in_res = res[inside.(e)]
# in_res_conv = in_res[in_res .<= 5e-3]
# in_eig_conv = in_eig[in_res .<= 5e-3]
# print("\nmax res inside: ")
# println(maximum(res[inside.(e)]))
# # print("max res inside non spurious: ")
# # println(maximum(in_res_conv))
# print("number inside : ")
# println(size(e[inside.(e)])[1])
# # print("number inside converged : ")
# # println(size(in_eig_conv)[1])
# display(in_eig)
# println()
# display(in_res)

# inside(x) = in_contour(x, complex(1.0, 1.0), 0.5)
#
# nep = PEP(bf_arrays())
# # e,v = polyeig(nep)
# e, v = contour_block_SS(nep, radius=R, σ=C, N=2^8, k=30,K=8)
# res = zeros(size(e,1))
# for i=1:size(e,1)
#     res[i] = norm(compute_Mlincomb(nep, e[i], v[:,i]))
# end
#
# info(e, v, res, C, R)
# in_eig = e[inside.(e)]
# display(in_eig)
#
# println()
#
# e, v, res = block_SS!(T, rand(ComplexF64,64,30), 2^6, 4, c=C, r=R)
# # e, v, res = nlfeast_moments_SS!(T, rand(ComplexF64,64,16), 2^6, 1, c=C, r=R, moments=2, debug=true, ϵ=10e-16, spurious=1e-3)
# info(e, v, res, C, R)
# inside(x) = in_contour(x, complex(1.0, 1.0), 0.5)
# in_eig = e[inside.(e)]
# in_res = res[inside.(e)]
# display(in_eig)
# display(in_res)


text:
# Basic goto tests

function goto_test1()
    @goto a
    return false
    @label a
    return true
end
@test goto_test1()

@test_throws ErrorException eval(
    quote
        function goto_test2()
            @goto a
            @label a
            @label a
            return
        end
    end)

@test_throws ErrorException eval(
    quote
        function goto_test3()
            @goto a
            return
        end
    end)

@test_throws ErrorException eval(
    quote
        function goto_test4()
            @goto a
            try
                @label a
            catch
            end
        end
    end)

# test that labels in macros are reassigned
macro goto_test5_macro()
    @label a
end

@test_throws ErrorException eval(
    quote
        function goto_test5()
            @goto a
            @goto_test5_macro
            return
        end
    end)

@test_throws ErrorException eval(
    quote
        function goto_test6()
            try
                @goto a
            finally
            end
            @label a
            return
        end
    end)

function goto_test6()
    @goto a
    @label a
end

@test goto_test6() == nothing

function goto_test7(x)
    @label a
    if x
        @goto a
    end
end

@test goto_test7(false) == nothing

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
# TODO: could also expose the channel options if someone wanted
# to define those; can probably leave for later until feedback has
# been received.

# Specifications: RSS 2.0 -- https://cyber.harvard.edu/rss/rss.html#sampleFiles
# steps:
# 0. check if the relevant variables are defined otherwise don't generate the RSS
# 1. is there an RSS file?
#  --> remove it and create a new one (bc items may have been updated)
# 2. go over all pages
#  --> is there a rss var?
#  NO  --> skip
#  YES --> recuperate the `fd_ctime` and `fd_mtime` and add a RSS channel object
# 3. save the file

struct RSSItem
    # -- required fields
    title::String
    link::String
    description::String  # note: should not contain

    # -- optional fields
    author::String       # note: should be a valid email
    category::String
    comments::String     # note: should be a valid URL
    enclosure::String
    # guid == link
    pubDate::Date        # note: should respect RFC822 (https://www.w3.org/Protocols/rfc822/)
end

const RSS_DICT = LittleDict{String,RSSItem}()

"""
    jor(a, b)

Convenience function for fallback fields.
"""
jor(a::String, b::String) = ifelse(isempty(locvar(a)), locvar(b), locvar(a))

"""
    remove_html_ps

Convenience function to remove
 and
 in RSS description (not supposed to
happen).
"""
remove_html_ps(s::String)::String = replace(s, r"</?p>" => "")

"""
$SIGNATURES

RSS should not contain relative links so this finds relative links and prepends
them with the canonical link.
"""
fix_relative_links(s::String, link::String) =
    replace(s, r"(href|src)\s*?=\s*?\"\/" => SubstitutionString("\\1=\"$link"))

"""
$SIGNATURES

Create an `RSSItem` out of the provided fields defined in the page vars.
"""
function add_rss_item()
    link  = url_curpage()
    title = jor("rss_title", "title")
    descr = jor("rss", "rss_description")

    descr = fd2html(descr; internal=true) |> remove_html_ps

    author    = locvar(:rss_author)
    category  = locvar(:rss_category)
    comments  = locvar(:rss_comments)
    enclosure = locvar(:rss_enclosure)

    pubDate = locvar(:rss_pubdate)
    if pubDate == Date(1)
        pubDate = locvar(:date)
        if !isa(pubDate, Date) || pubDate == Date(1)
            pubDate = Date(locvar(:fd_mtime_raw))
        end
    end

    # warning for title which should really be defined
    isnothing(title) && (title = "")
    isempty(title)   && print_warning("""
        An RSS description was found but without title for page '$link'.
        """)

    res = RSS_DICT[link] = RSSItem(title, link, descr, author,
                                   category, comments, enclosure, pubDate)
    return res
end

"""
$SIGNATURES

Extract the entries from RSS_DICT and assemble the RSS. If the dictionary is empty, nothing
is generated.
"""
function rss_generator()::Nothing
    # is there anything to go in the RSS feed?
    isempty(RSS_DICT) && return nothing

    # are the basic defs there? otherwise warn and break
    rss_title = globvar("website_title")
    rss_descr = globvar("website_descr")
    rss_link  = globvar("website_url")

    if any(isempty, (rss_title, rss_descr, rss_link))
        print_warning("""
            RSS items were found but the RSS feed is improperly described:
            at least one of the following variables have not been defined in
            your 'config.md': 'website_title', 'website_descr', 'website_url'.
            The feed will not be (re)generated.
            \nRelevant pointer:
            $POINTER_PV
            """)
        return nothing
    end

    endswith(rss_link, "/") || (rss_link *= "/")
    rss_descr = fd2html(rss_descr; internal=true) |> remove_html_ps

    # is there an RSS file already? if so remove it
    rss_path = joinpath(PATHS[:site], "feed.xml")
    isfile(rss_path) && rm(rss_path)

    # create a buffer which will correspond to the output
    rss_buff = IOBuffer()
    write(rss_buff,
        """
        <rss version="2.0" xmlns:atom="http://www.w3.org/2005/Atom">
        <channel>
          <title>$rss_title</title>
          <description><![CDATA[$(fix_relative_links(rss_descr, rss_link))]]></description>
          <link>$rss_link</link>
          <atom:link href="$(rss_link)feed.xml" rel="self" type="application/rss+xml" />
        """)
    # loop over items
    for (k, v) in RSS_DICT
        full_link = rss_link
        if startswith(v.link, "/")
            full_link *= v.link[2:end]
        else
            full_link *= v.link
        end
        write(rss_buff,
          """
            <item>
              <title>$(v.title)</title>
              <link>$(full_link)</link>
              <description><![CDATA[$(fix_relative_links(v.description, rss_link))
Read more]]></description>
          """)
        for elem in (:author, :category, :comments, :enclosure)
            e = getproperty(v, elem)
            isempty(e) || write(rss_buff,
              """
                  <$elem>$e</$elem>
              """)
        end
        write(rss_buff,
          """
              <guid>$(full_link)</guid>
              <pubDate>$(Dates.format(v.pubDate, "e, d u Y")) 00:00:00 UT</pubDate>
            </item>
          """)
    end
    # finalize
    write(rss_buff,
        """
        </channel>
        </rss>
        """)
    write(rss_path, take!(rss_buff))

    return nothing
end


text:
"""
    Torque(;name)

Input signal acting as external torque on a flange
"""
function Torque(;name, use_support=false)
    @named partial_element = PartialElementaryOneFlangeAndSupport2(use_support=use_support)
    @unpack flange = partial_element
    @named tau = RealInput() # Accelerating torque acting at flange (= -flange.tau)
    eqs = [flange.tau ~ -tau.u]
    return extend(ODESystem(eqs, t, [], []; name=name, systems=[tau]), partial_element)
end


text:
module Newton

# Write your package code here.
using Statistics, LinearAlgebra, ForwardDiff

function newtonroot(f,f′;x₀,tol=1E-7,maxiter=1000)
    println("the first dispach")
    x=zeros(maxiter+1)
    x_old=x₀
    i=1
    norm_dif=Inf
    while i <= maxiter &&  norm_dif>=tol
        x_new=x_old - f(x_old)/f′(x_old)
        i+=1
        norm_dif=norm(x_new-x_old)
        x_old=x_new
    end
    if i == maxiter+1
        return nothing
    else
        return (iteration=i, norm_dif=norm_dif, x=x_old)
    end
end

function newtonroot(f;x₀,tol = 1E-7,maxiter = 1000)
    println("the second dispach")
    x=zeros(maxiter+1)
    x_old=x₀
    i=1
    norm_dif=Inf
    f_prime(x)=ForwardDiff.derivative(f,x)
    while i <= maxiter &&  norm_dif>=tol
        x_new=x_old - f(x_old)/f_prime(x_old)
        i+=1
        norm_dif=norm(x_new-x_old)
        x_old=x_new
    end
    if i == maxiter+1
        return nothing
    else
        return (iteration=i, norm_dif=norm_dif, x=x_old)
    end
end

export newtonroot
end


text:
# Helper functions
sameSense(pra::Int, ra::Int) = mod(pra,2)==mod(ra,2)
downSense(ra::Int) = (ra>0) .& (mod(ra,2)==1)
upSense(ra::Int) = (ra>0) .& (mod(ra,2)==0)

# Reward function for VerticalCAS MDP
function POMDPs.reward(mdp::VerticalCAS_MDP, s::stateType, ra::actType)
    h = s[1]; vown = s[2]; vint = s[3]; pra = s[4]; resp = s[5];
    tau = mdp.currentTau
    r = 0.0
    sep = abs(h)
    closure = abs(vint-vown)
    crossing = ((h<0) .& downSense(ra)) .| ((h>0) .& upSense(ra))
    deltaVown = 0.0
    corrective = false
    preventitive = false
    weakening = false
    strengthening = false
    reversal = false

    if ra>COC
        if pra>COC
            reversal = mod(pra,2)!=mod(ra,2)
        end
        if !reversal
            weakening = pra>ra
            strengthening = pra<ra
        end
        vLow, vHigh = mdp.velRanges[ra]
        corrective = (vown>vLow) .& (vown < vHigh)
        preventitive = !corrective
        if corrective
            if downSense(ra)
                deltaVown = abs(vLow-vown)
            else
                deltaVown = abs(vHigh-vown)
            end
        end
    end
    lolo = (ra==DNC) .| (ra==DND)

    if (sep<=175) .& (tau==0)
        r-=1.0
    end

    if mdp.allowedTrans[pra][ra+1]==0
        r-=1.0
    end

    if crossing
        if preventitive
            r-=1.0
        end
        if sep>500
            r-=0.01
        end
    end
    if corrective
        r-=1e-5
        if (sep>650) .& (closure<2000.0/60.0)
            r-= 0.1
        end
        if (sep>1000) .& (closure<4000.0/60.0)
            r-=0.03
        end
    elseif preventitive
        if (sep>650) .& (closure<2000.0/60.0)
            r-=0.01
        end
    end
    if reversal
        r-= 8e-3
    end
    if strengthening
        r-=5e-3
    end
    if weakening
        r-=1e-3
    end
    if lolo
        r-=1e-4
        if closure > 3000.0/60.0
            r-=5e-4
        end
    elseif (ra!=COC) .& (closure > 3000.0/60.0)
        r-=1.5e-3
    end
    if closure < 3000.0/60.0
        r-=2.3e-3
    end
    if ra==COC
        r+=1e-9
    else
        r-=3e-5*deltaVown
        if closure > 3000.0/60.0
            r-=1.5e-3
        end
    end
    return r
end

max_stars_repo_name: alisamiii/VerticalCAS

max_issues_repo_name: alisamiii/VerticalCAS

max_forks_repo_name: bethelhall/CoCoCAS


text:
__precompile__()

module Timeout

using Compat.Distributed
using Compat.Dates
using Compat: @warn

export ptimeout, @ptimeout

if isdefined(Base, :_UVError)
    using Base: _UVError
else
    const _UVError = Base.UVError
end

if isdefined(Base, :SIGTERM)
    using Base: SIGTERM
else
    const SIGTERM = 15 # just taking the value in 1.0
end

"""
    ptimeout(f, limit; worker=1, poll=0.5, verbose=true)

Run the given function `f` in a separate process on worker `worker` for a maximum time
of `limit`, which can be a number of seconds or a `Dates.Period`. If the time limit is
reached, the remote process will be interrupted, first by sending SIGINTs and then by
sending a SIGTERM to forcibly kill the process if the interrupts are ineffective. If `f`
completed without timing out, `true` is returned, otherwise `false`.

The keyword argument `worker` selects the worker by process ID from the pool of worker
processes known to Julia, i.e. that which is returned by `workers()`. The chosen worker
cannot be the current process.

`poll` is the number of seconds to wait before rechecking whether `f` has finished
executing on the remote worker. For long running jobs, you may want to set this higher
than the default 0.5, which will poll every half second.

If `verbose` is `true`, a warning will be logged by the calling process before attempting
to interrupt.

!!! warn
    The return value refers ONLY to whether the computation timed out. Computation of
    `f` may have terminated due to an error that's unrelated to an interrupt sent by
    `ptimeout`, and this will NOT be reflected in the return value.
"""
function ptimeout(f::Function, secs::Real; worker=1, poll=0.5, verbose=true)
    nprocs() > 1 || throw(ArgumentError("No worker processes available"))
    worker in workers() || throw(ArgumentError("Unknown worker process ID: $worker"))
    worker == myid() && throw(ArgumentError("Can't run ptimeout on the current process"))
    poll > 0 || throw(ArgumentError("Can't poll every $poll seconds"))

    # We need the worker process to be on the same host as the calling process, otherwise
    # sending a SIGTERM to the result of getpid might kill off something local
    if gethostname() != remotecall_fetch(gethostname, worker)
        throw(ArgumentError("Can't run ptimeout with a worker on a different host"))
    end

    # Now start by getting the OS process ID for the worker so that we have something to
    # forcibly kill if need be
    ospid = remotecall_fetch(getpid, worker)

    # Run the function on the given worker, with a channel for communicating with the
    # process so that checking isready won't block
    channel = Channel(1)
    @async put!(channel, remotecall_fetch(f, worker))
    timedwait(()->isready(channel), float(secs), pollint=float(poll))
    isready(channel) && return true
    verbose && @warn "Time limit for computation exceeded. Interrupting..."
    patience = 10
    while !isready(channel) && (patience -= 1) > 0
        interrupt(worker)
    end
    # If our interrupts didn't work, forcibly kill the process
    if !isready(channel)
        rc = ccall(:uv_kill, Cint, (Cint, Cint), ospid, SIGTERM)
        rc == 0 || throw(_UVError("kill", rc))
    end
    close(channel)
    false
end

ptimeout(f::Function, time::Period; kwargs...) =
    ptimeout(f, Dates.value(convert(Second, time)); kwargs...)

"""
    @ptimeout worker limit expr

Wrap the given expression `expr` in a function and run it on a remote worker up to the
given amount of time, `limit`. This macro is a thin convenience wrapper around the
[`ptimeout`](@ref) function.
"""
macro ptimeout(worker, limit, expr)
    :(ptimeout(()->$expr, $limit, worker=$worker))
end

end # module


text:
# Demonstrate how to optimize L2 penalty parameter with LogisticOptTools
# We follow the method introduced in:
# Optimizing for Generalization in Machine Learning with Cross-Validation Gradients
# S Barratt, R Sharma
# arXiv preprint arXiv:1805.07072

using Optim
using LogisticOptTools
using Plots

const LOT = LogisticOptTools

# Select dataset where to optimize l2 penalty parameter
SVM_DATASET = "diabetes"
dataset = LOT.LogitData(SVM_DATASET)

# Instantiate optimizer. We have an outer and an inner optimizers there.
options = Optim.Options(iterations=250, g_tol=1e-5, show_trace=false)
optimizer = LOT.L2PenaltyOptimizer(3, false, false, BFGS(), LBFGS(), options)

# Set initial parameter
γ0 = 1.0
# Optimize (be careful, problem is not necessary convex)!
l♯, γ♯ = LOT.fit!(optimizer, dataset.X, dataset.y, γ0)

# Show effectiveness of optimal parameter found
costs = []
# Range of γ to explore
t = range(-10, step=.1, stop=10.)
for i in  t
    γ = 10^i
    c, dc = LOT.crossvalid_loss(dataset.X, dataset.y, γ, penopt)
    push!(costs, c)
end

plot(t, costs)
vline!([log10(γ♯)])
xlabel!("ℓ2 parameter (log scale)")
ylabel!("Costs")


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule DynarePreprocessor_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("DynarePreprocessor")
JLLWrappers.@generate_main_file("DynarePreprocessor", UUID("23afba7c-24e5-5ee2-bc2c-b42e07f0492a"))
end  # module DynarePreprocessor_jll


text:
using ClimateMachine
ClimateMachine.init()
using ClimateMachine.Atmos
using ClimateMachine.ConfigTypes
using ClimateMachine.DGmethods
using ClimateMachine.DGmethods.NumericalFluxes
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Topologies
using ClimateMachine.MoistThermodynamics
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates
import ClimateMachine.Atmos:
    MoistureModel,
    temperature,
    pressure,
    soundspeed,
    total_specific_enthalpy,
    thermo_state

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using LinearAlgebra
using MPI
using StaticArrays
using Test

include("mms_solution_generated.jl")

"""
    MMSDryModel

Assumes the moisture components is in the dry limit.
"""
struct MMSDryModel <: MoistureModel end

function total_specific_enthalpy(
    bl::AtmosModel,
    moist::MMSDryModel,
    state::Vars,
    aux::Vars,
)
    zero(eltype(state))
end
function thermo_state(
    bl::AtmosModel,
    moist::MMSDryModel,
    state::Vars,
    aux::Vars,
)
    PS = typeof(bl.param_set)
    return PhaseDry{eltype(state), PS}(
        bl.param_set,
        internal_energy(bl, state, aux),
        state.ρ,
    )
end
function pressure(bl::AtmosModel, moist::MMSDryModel, state::Vars, aux::Vars)
    T = eltype(state)
    γ = T(7) / T(5)
    ρinv = 1 / state.ρ
    return (γ - 1) * (state.ρe - ρinv / 2 * sum(abs2, state.ρu))
end
function soundspeed(bl::AtmosModel, moist::MMSDryModel, state::Vars, aux::Vars)
    T = eltype(state)
    γ = T(7) / T(5)
    ρinv = 1 / state.ρ
    p = pressure(bl, bl.moisture, state, aux)
    sqrt(ρinv * γ * p)
end

function mms3_init_state!(bl, state::Vars, aux::Vars, (x1, x2, x3), t)
    state.ρ = ρ_g(t, x1, x2, x3, Val(3))
    state.ρu = SVector(
        U_g(t, x1, x2, x3, Val(3)),
        V_g(t, x1, x2, x3, Val(3)),
        W_g(t, x1, x2, x3, Val(3)),
    )
    state.ρe = E_g(t, x1, x2, x3, Val(3))
end

function mms3_source!(
    bl,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)
    x1, x2, x3 = aux.coord
    source.ρ = Sρ_g(t, x1, x2, x3, Val(3))
    source.ρu = SVector(
        SU_g(t, x1, x2, x3, Val(3)),
        SV_g(t, x1, x2, x3, Val(3)),
        SW_g(t, x1, x2, x3, Val(3)),
    )
    source.ρe = SE_g(t, x1, x2, x3, Val(3))
end

function main()
    FT = Float64

    # DG polynomial order
    N = 4

    t0 = FT(0)
    timeend = FT(1)

    ode_dt = 0.00125
    nsteps = ceil(Int64, timeend / ode_dt)
    ode_dt = timeend / nsteps

    ode_solver = ClimateMachine.ExplicitSolverType(
        solver_method = LSRK54CarpenterKennedy,
    )

    expected_result = FT(3.403104838700577e-02)

    model = AtmosModel{FT}(
        AtmosLESConfigType,
        param_set;
        orientation = NoOrientation(),
        ref_state = NoReferenceState(),
        turbulence = ConstantViscosityWithDivergence(FT(μ_exact)),
        moisture = MMSDryModel(),
        source = mms3_source!,
        boundarycondition = InitStateBC(),
        init_state_conservative = mms3_init_state!,
    )

    brickrange = (
        range(FT(0); length = 5, stop = 1),
        range(FT(0); length = 5, stop = 1),
        range(FT(0); length = 5, stop = 1),
    )
    topl = BrickTopology(
        MPI.COMM_WORLD,
        brickrange,
        periodicity = (false, false, false),
    )
    warpfun =
        (x1, x2, x3) -> begin
            (
                x1 + (x1 - 1 / 2) * cos(2 * π * x2 * x3) / 4,
                x2 + exp(sin(2π * (x1 * x2 + x3))) / 20,
                x3 + x1 / 4 + x2^2 / 2 + sin(x1 * x2 * x3),
            )
        end
    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ClimateMachine.array_type(),
        polynomialorder = N,
        meshwarp = warpfun,
    )
    driver_config = ClimateMachine.DriverConfiguration(
        AtmosLESConfigType(),
        "MMS3",
        N,
        FT,
        ClimateMachine.array_type(),
        ode_solver,
        model,
        MPI.COMM_WORLD,
        grid,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
        ClimateMachine.AtmosLESSpecificInfo(),
    )

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_solver_type = ode_solver,
        ode_dt = ode_dt,
    )
    dg = DGModel(
        model,
        grid,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    ClimateMachine.invoke!(solver_config)

    Qe = init_ode_state(dg, timeend)
    result = euclidean_distance(solver_config.Q, Qe)
    @test result ≈ expected_result
end

main()


text:
@testset "Moses Problems" begin
    (A, B, H, K, a, e, r, t, w, x, y, z, ) = @variables A B H K a e r t w x y z

    #= ::Package:: =#

    #= ::Title:: =#
    #=Joel*Moses - Symbolic*Integration*Ph.D. Thesis*(1967)=#

    #= ::Section::Closed:: =#
    #=Chapter*2 - How*SIN*differs*from*SAINT=#

    @test_int [cot(x)^4, x, 3, x + cot(x) - cot(x)^3/3]
    @test_int [1/(x^4*(1 + x^2)), x, 3, -(1/(3*x^3)) + 1/x + atan(x)]
    @test_int [(x^2 + x)/sqrt(x), x, 2, (2*x^(3/2))/3 + (2*x^(5/2))/5]
    @test_int [cos(x), x, 1, sin(x)]
    @test_int [x*ℯ^x^2, x, 1, ℯ^x^2/2]
    @test_int [tan(x)*sec(x)^2, x, 2, sec(x)^2/2]
    @test_int [x*sqrt(1 + x^2), x, 1, (1/3)*(1 + x^2)^(3/2)]
    @test_int [sin(x)*ℯ^x, x, 1, (-(1/2))*ℯ^x*cos(x) + (1/2)*ℯ^x*sin(x)]

    #= ::Section::Closed:: =#
    #=Chapter*3 - SCHATCHEN - A*Matching*Program*for*Algebraic*Expressions=#

    @test_int [csc(x)^2*(cos(x)/sin(x)^2), x, 2, (-(1/3))*csc(x)^3]

    #= ::Section::Closed:: =#
    #=Chapter*4 - The*First*Stage*of*sin=#

    @test_int [sin(ℯ^x), x, 2, SinIntegral(ℯ^x)]
    @test_int [sin(y)/y, y, 1, SinIntegral(y)]

    @test_int [sin(x) + ℯ^x, x, 3, ℯ^x - cos(x)]
    @test_int [ℯ^x^2 + 2*x^2*ℯ^x^2, x, 4, ℯ^x^2*x]
    @test_int [(x + ℯ^x)^2, x, 5, -2*ℯ^x + ℯ^(2*x)/2 + 2*ℯ^x*x + x^3/3]
    @test_int [x^2 + 2*ℯ^x + ℯ^(2*x), x, 3, 2*ℯ^x + ℯ^(2*x)/2 + x^3/3]

    @test_int [sin(x)*cos(x), x, 2, sin(x)^2/2]
    @test_int [x*ℯ^x^2, x, 1, ℯ^x^2/2]
    @test_int [x*sqrt(1 + x^2), x, 1, (1/3)*(1 + x^2)^(3/2)]
    @test_int [ℯ^x/(1 + ℯ^x), x, 2, log(1 + ℯ^x)]
    @test_int [x^(3/2), x, 1, (2*x^(5/2))/5]
    @test_int [cos(2*x + 3), x, 1, (1/2)*sin(3 + 2*x)]
    @test_int [2*y*z*ℯ^(2*x), x, 2, ℯ^(2*x)*y*z]
    @test_int [cos(ℯ^x)^2*sin(ℯ^x)*ℯ^x, x, 3, (-(1/3))*cos(ℯ^x)^3]

    #= ::Section::Closed:: =#
    #=Chapter*4 - The*Second*Stage*of*sin=#

    @test_int [x*sqrt(x + 1), x, 2, (-(2/3))*(1 + x)^(3/2) + (2/5)*(1 + x)^(5/2)]
    @test_int [1/(x^4 - 1), x, 3, -(atan(x)/2) - atanh(x)/2]

    #= ::Subsection::Closed:: =#
    #=Method*1)  Elementary*function*of*exponentials=#

    @test_int [ℯ^x/(2 + 3*ℯ^(2*x)), x, 2, atan(sqrt(3/2)*ℯ^x)/sqrt(6)]
    @test_int [ℯ^(2*x)/(A + B*ℯ^(4*x)), x, 2, atan((sqrt(B)*ℯ^(2*x))/sqrt(A))/(2*sqrt(A)*sqrt(B))]
    @test_int [ℯ^(x + 1)/(1 + ℯ^x), x, 3, ℯ*log(1 + ℯ^x)]
    @test_int [10^x*ℯ^x, x, 1, (10*ℯ)^x/(1 + log(10))]

    #= ::Subsection::Closed:: =#
    #=Method*2)  Substitution*for*an*integral*power=#

    @test_int [x^3*sin(x^2), x, 3, (-(1/2))*x^2*cos(x^2) + sin(x^2)/2]
    @test_int [x^7/(x^12 + 1), x, 7, -(atan((1 - 2*x^4)/sqrt(3))/(4*sqrt(3))) - (1/12)*log(1 + x^4) + (1/24)*log(1 - x^4 + x^8)]
    @test_int [x^(3*a)*sin(x^(2*a)), x, 3, (im*x^(1 + 3*a)*Gamma((1/2)*(3 + 1/a), (-im)*x^(2*a)))/(((-im)*x^(2*a))^((1 + 3*a)/(2*a))*(4*a)) - (im*x^(1 + 3*a)*Gamma((1/2)*(3 + 1/a), im*x^(2*a)))/((im*x^(2*a))^((1 + 3*a)/(2*a))*(4*a))]

    #= ::Subsection::Closed:: =#
    #=Method*3)  Substitution*for*a*rational*root*of*a*linear*function*of*x=#

    @test_int [cos(sqrt(x)), x, 3, 2*cos(sqrt(x)) + 2*sqrt(x)*sin(sqrt(x))]
    @test_int [x*sqrt(x + 1), x, 2, (-(2/3))*(1 + x)^(3/2) + (2/5)*(1 + x)^(5/2)]
    @test_int [1/(x^(1/2) + x^(1/3)), x, 4, 6*x^(1/6) - 3*x^(1/3) + 2*sqrt(x) - 6*log(1 + x^(1/6))]
    @test_int [sqrt((x + 1)/(2*x + 3)), x, 4, (1/2)*sqrt(1 + x)*sqrt(3 + 2*x) - asinh(sqrt(2)*sqrt(1 + x))/(2*sqrt(2))]

    #= ::Subsection::Closed:: =#
    #=Method*4)  Binomial - Chebyschev=#

    @test_int [x^4/(1 - x^2)^(5/2), x, 3, x^3/(3*(1 - x^2)^(3/2)) - x/sqrt(1 - x^2) + asin(x)]
    @test_int [x^(1/2)*(1 + x)^(5/2), x, 6, (5/64)*sqrt(x)*sqrt(1 + x) + (5/32)*x^(3/2)*sqrt(1 + x) + (5/24)*x^(3/2)*(1 + x)^(3/2) + (1/4)*x^(3/2)*(1 + x)^(5/2) - (5*asinh(sqrt(x)))/64]

    #= ::Subsection::Closed:: =#
    #=Method*5)  Arctrigonometric*substitutions=#

    @test_int [x^4/(1 - x^2)^(5/2), x, 3, x^3/(3*(1 - x^2)^(3/2)) - x/sqrt(1 - x^2) + asin(x)]
    @test_int [sqrt(A^2 + B^2 - B^2*y^2)/(1 - y^2), y, 5, B*atan((B*y)/sqrt(A^2 + B^2 - B^2*y^2)) + A*atanh((A*y)/sqrt(A^2 + B^2 - B^2*y^2))]

    #= ::Subsection::Closed:: =#
    #=Method*6)  Elementary*function*of*trigonometric*functions=#

    @test_int [sin(x)^2, x, 2, x/2 - (1/2)*cos(x)*sin(x)]
    @test_int [sqrt(A^2 + B^2*sin(x)^2)/sin(x), x, 6, (-B)*atan((B*cos(x))/sqrt(A^2 + B^2*sin(x)^2)) - A*atanh((A*cos(x))/sqrt(A^2 + B^2*sin(x)^2)), (-B)*atan((B*cos(x))/sqrt(A^2 + B^2 - B^2*cos(x)^2)) - A*atanh((A*cos(x))/sqrt(A^2 + B^2 - B^2*cos(x)^2))]
    @test_int [1/(1 + cos(x)), x, 1, sin(x)/(1 + cos(x))]

    #= ::Subsection::Closed:: =#
    #=Method*7)  Rational*function*times*an*exponential=#

    @test_int [x*ℯ^x, x, 2, -ℯ^x + ℯ^x*x]
    @test_int [(x/(x + 1)^2)*ℯ^x, x, 1, ℯ^x/(1 + x)]
    @test_int [(1 + 2*x^2)*ℯ^x^2, x, 5, ℯ^x^2*x]
    @test_int [ℯ^x^2, x, 1, (1/2)*sqrt(pi)*erfi(x)]
    @test_int [ℯ^x/x, x, 1, ExpIntegralEi(x)]

    #= ::Subsection::Closed:: =#
    #=Method*8)  Rational*functions=#

    @test_int [x/(x^3 + 1), x, 6, -(atan((1 - 2*x)/sqrt(3))/sqrt(3)) - (1/3)*log(1 + x) + (1/6)*log(1 - x + x^2)]
    @test_int [1/(x^6 - 1), x, 10, -(atan((sqrt(3)*x)/(1 - x^2))/(2*sqrt(3))) - atanh(x)/3 - (1/6)*atanh(x/(1 + x^2)), atan((1 - 2*x)/sqrt(3))/(2*sqrt(3)) - atan((1 + 2*x)/sqrt(3))/(2*sqrt(3)) - atanh(x)/3 + (1/12)*log(1 - x + x^2) - (1/12)*log(1 + x + x^2)]
    @test_int [1/((B^2 - A^2)*x^2 - A^2*B^2 + A^4), x, 1, atanh(x/A)/(A*(A^2 - B^2))]

    #= ::Subsection::Closed:: =#
    #=Method*9)  Rational*function*times*a*log*or*arctrigonometric*function=#

    @test_int [x*log(x), x, 1, -(x^2/4) + (1/2)*x^2*log(x)]
    @test_int [x^2*asin(x), x, 4, sqrt(1 - x^2)/3 - (1/9)*(1 - x^2)^(3/2) + (1/3)*x^3*asin(x)]
    @test_int [1/(x^2 + 2*x + 1), x, 2, -(1/(1 + x))]

    #= ::Subsection::Closed:: =#
    #=Method*10)  Rational*function*times*an*elementary*function*of*log(a+b*x)=#

    @test_int [log(x)/(log(x) + 1)^2, x, 7, x/(1 + log(x))]
    @test_int [1/(x*(1 + log(x)^2)), x, 2, atan(log(x))]
    @test_int [1/log(x), x, 1, LogIntegral(x)]

    #= ::Subsection::Closed:: =#
    #=Method*11)  Expansion*of*the*integrand=#

    @test_int [x*(cos(x) + sin(x)), x, 6, cos(x) - x*cos(x) + sin(x) + x*sin(x)]
    @test_int [(x + ℯ^x)/ℯ^x, x, 4, -ℯ^(-x) + x - x/ℯ^x]
    @test_int [x*(1 + ℯ^x)^2, x, 6, -2*ℯ^x - ℯ^(2*x)/4 + 2*ℯ^x*x + (1/2)*ℯ^(2*x)*x + x^2/2]

    #= ::Section::Closed:: =#
    #=Chapter*4 - The*Third*Stage*of*sin=#

    @test_int [x*cos(x), x, 2, cos(x) + x*sin(x)]
    @test_int [cos(sqrt(x)), x, 3, 2*cos(sqrt(x)) + 2*sqrt(x)*sin(sqrt(x))]

    #= ::Subsection::Closed:: =#
    #=The*Integration-by-Parts*Methods=#

    @test_int [x*cos(x), x, 2, cos(x) + x*sin(x)]
    @test_int [x*log(x)^2, x, 2, x^2/4 - (1/2)*x^2*log(x) + (1/2)*x^2*log(x)^2]

    #= ::Subsection::Closed:: =#
    #=The*Derivative-divides*Method=#

    @test_int [cos(x)*(1 + sin(x)^3), x, 2, sin(x) + sin(x)^4/4]
    @test_int [1/(x*(1 + log(x)^2)), x, 2, atan(log(x))]
    @test_int [1/(sqrt(1 - x^2)*(1 + asin(x)^2)), x, 2, atan(asin(x))]
    @test_int [sin(x)/(sin(x) + cos(x)), x, 2, x/2 - (1/2)*log(cos(x) + sin(x))]

    #= ::Subsection::Closed:: =#
    #=An*Example*of*SIN's*Performance=#

    @test_int [-sqrt(A^2 + B^2*(1 - y^2))/(1 - y^2), y, 6, (-B)*atan((B*y)/sqrt(A^2 + B^2 - B^2*y^2)) - A*atanh((A*y)/sqrt(A^2 + B^2 - B^2*y^2))]
    @test_int [(-(A^2 + B^2))*(cos(z)^2/(B*(1 - ((A^2 + B^2)/B^2)*sin(z)^2))), z, 5, (-B)*z - A*atanh((A*tan(z))/B)]
    @test_int [-(A^2 + B^2)/(B*(1 + w^2)^2*(1 - ((A^2 + B^2)/B^2)*(w^2/(1 + w^2)))), w, 6, (-B)*atan(w) - A*atanh((A*w)/B)]
    @test_int [(-B)*((A^2 + B^2)/((1 + w^2)*(B^2 - A^2*w^2))), w, 4, (-B)*atan(w) - A*atanh((A*w)/B)]

    #= ::Subsection::Closed:: =#
    #=SAINT*and*SIN*solutions*of*the*same*problem=#

    @test_int [x^4/(1 - x^2)^(5/2), x, 3, x^3/(3*(1 - x^2)^(3/2)) - x/sqrt(1 - x^2) + asin(x)]
    @test_int [sin(y)^4/cos(y)^4, y, 3, y - tan(y) + tan(y)^3/3]
    @test_int [z^4/(1 + z^2), z, 3, -z + z^3/3 + atan(z)]

    #= ::Section::Closed:: =#
    #=Chapter*5 - The*Edge*Heuristic=#

    @test_int [(2*x^2 + 1)*ℯ^x^2, x, 5, ℯ^x^2*x]
    @test_int [((2*x^6 + 5*x^4 + x^3 + 4*x^2 + 1)/(x^2 + 1)^2)*ℯ^x^2, x, 10, ℯ^x^2*x + ℯ^x^2/(2*(1 + x^2))]
    @test_int [1/ℯ^1/ℯ^x, x, 1, -ℯ^(-1 - x)]
    @test_int [(x + 1/x)*log(x), x, 5, -(x^2/4) + (1/2)*x^2*log(x) + log(x)^2/2]
    @test_int [x/(1 + x^4), x, 2, atan(x^2)/2]
    @test_int [x^5/(1 + x^4), x, 3, x^2/2 - atan(x^2)/2]
    @test_int [1/(1 + tan(x)^2), x, 3, x/2 + (1/2)*cos(x)*sin(x)]
    @test_int [x^4/(1 - x^2)^(5/2), x, 3, x^3/(3*(1 - x^2)^(3/2)) - x/sqrt(1 - x^2) + asin(x)]
    @test_int [-x^2/(1 - x^2)^(3/2), x, 2, -(x/sqrt(1 - x^2)) + asin(x)]
    @test_int [sin(x)*ℯ^x, x, 1, (-(1/2))*ℯ^x*cos(x) + (1/2)*ℯ^x*sin(x)]

    #= ::Section::Closed:: =#
    #=Appendix*C - Slagle's*Thesis*Integration*Problems=#

    @test_int [1/x, x, 1, log(x)]
    @test_int [sec(2*t)/(1 + sec(t)^2 + 3*tan(t)), t, 4, (-(1/12))*log(cos(t) - sin(t)) - (1/4)*log(cos(t) + sin(t)) + (1/3)*log(2*cos(t) + sin(t)) - 1/(2*(1 + tan(t)))]
    @test_int [1/sec(x)^2, x, 2, x/2 + (1/2)*cos(x)*sin(x)]
    @test_int [(x^2 + 1)/sqrt(x), x, 2, 2*sqrt(x) + (2*x^(5/2))/5]
    @test_int [x/sqrt(x^2 + 2*x + 5), x, 3, sqrt(5 + 2*x + x^2) - asinh((1 + x)/2)]
    @test_int [sin(x)^2*cos(x), x, 2, sin(x)^3/3]
    @test_int [ℯ^x/(1 + ℯ^x), x, 2, log(1 + ℯ^x)]
    @test_int [ℯ^(2*x)/(1 + ℯ^x), x, 3, ℯ^x - log(1 + ℯ^x)]
    @test_int [1/(1 - cos(x)), x, 1, -(sin(x)/(1 - cos(x)))]
    @test_int [tan(x)*sec(x)^2, x, 2, sec(x)^2/2]
    @test_int [x*log(x), x, 1, -(x^2/4) + (1/2)*x^2*log(x)]
    @test_int [sin(x)*cos(x), x, 2, sin(x)^2/2]
    @test_int [(x + 1)/sqrt(2*x - x^2), x, 3, -sqrt(2*x - x^2) - 2*asin(1 - x)]
    @test_int [2*(ℯ^x/(2 + 3*ℯ^(2*x))), x, 3, sqrt(2/3)*atan(sqrt(3/2)*ℯ^x)]
    @test_int [x^4/(1 - x^2)^(5/2), x, 3, x^3/(3*(1 - x^2)^(3/2)) - x/sqrt(1 - x^2) + asin(x)]
    @test_int [ℯ^(6*x)/(ℯ^(4*x) + 1), x, 3, ℯ^(2*x)/2 - (1/2)*atan(ℯ^(2*x))]
    @test_int [log(2 + 3*x^2), x, 3, -2*x + 2*sqrt(2/3)*atan(sqrt(3/2)*x) + x*log(2 + 3*x^2)]

    #= ::Section::Closed:: =#
    #=Appendix*D - MacIntosh*Integration*Problems=#

    @test_int [1/(r*sqrt(2*H*r^2 - a^2)), x, 1, x/(r*sqrt(-a^2 + 2*H*r^2))]
    @test_int [1/(r*sqrt(2*H*r^2 - a^2 - e^2)), x, 1, x/(r*sqrt(-a^2 - e^2 + 2*H*r^2))]
    @test_int [1/(r*sqrt(2*H*r^2 - a^2 - 2*K*r^4)), x, 1, x/(r*sqrt(-a^2 + 2*H*r^2 - 2*K*r^4))]
    @test_int [1/(r*sqrt(2*H*r^2 - a^2 - e^2 - 2*K*r^4)), x, 1, x/(r*sqrt(-a^2 - e^2 + 2*H*r^2 - 2*K*r^4))]
    @test_int [1/(r*sqrt(2*H*r^2 - a^2 - 2*K*r)), x, 1, x/(r*sqrt(-a^2 - 2*r*(K - H*r)))]
    @test_int [1/(r*sqrt(2*H*r^2 - a^2 - e^2 - 2*K*r)), x, 1, If(13>=8, x/(r*sqrt(-a^2 - e^2 - 2*r*(K - H*r))), x/(r*sqrt(-a^2 - e^2 - 2*K*r + 2*H*r^2)))]
    @test_int [r/sqrt(2*ℯ*r^2 - a^2), x, 1, (r*x)/sqrt(-a^2 + 2*ℯ*r^2)]
    @test_int [r/sqrt(2*ℯ*r^2 - a^2 - e^2), x, 1, (r*x)/sqrt(-a^2 - e^2 + 2*ℯ*r^2)]
    @test_int [r/sqrt(2*ℯ*r^2 - a^2 - 2*K*r^4), x, 1, (r*x)/sqrt(-a^2 + 2*ℯ*r^2 - 2*K*r^4)]
    @test_int [r/sqrt(2*ℯ*r^2 - a^2 - e^2 - 2*K*r^4), x, 1, (r*x)/sqrt(-a^2 - e^2 + 2*ℯ*r^2 - 2*K*r^4)]
    @test_int [r/sqrt(2*H*r^2 - a^2 - e^2 - 2*K*r), x, 1, If(13>=8, (r*x)/sqrt(-a^2 - e^2 - 2*r*(K - H*r)), (r*x)/sqrt(-a^2 - e^2 - 2*K*r + 2*H*r^2))]
    end

max_stars_repo_path: test/0 Independent test suites/Moses Problems.jl


text:
@testset "learn" begin
    bper = BPELearner([joinpath(dirname(@__FILE__), "data/corpus.en")], 1000)
    learn!(bper)
    bper_result = emit(bper)
    open(joinpath(dirname(@__FILE__), "data/bpe.ref")) do ref
        _v, lines = Iterators.peel(readlines(ref))
        result = map((x)->Tuple(split(x)), lines)
        @test length(bper.bpe) == length(result)

        for (lr, res) ∈ zip(bper_result, result)
            @test lr == res
        end
    end

    bpefile = joinpath(dirname(@__FILE__), "data/bpe.out")
    emit(bper, bpefile; comment = "this is for testing...")

    open(bpefile) do bf
        _h = readline(bf)
        for (i, line) ∈ enumerate(eachline(bf))
            pair = Tuple(split(line, " "))
            @test pair == bper_result[i]
        end
    end
end


text:
# https://github.com/goedman/Stan.jl/blob/master/Examples/Mamba/Binormal/binormal.jl

using CmdStan, Turing, ContinuousBenchmarks

include(joinpath(ContinuousBenchmarks.STAN_MODELS_DIR, "binormal-stan.model.jl"))
include(joinpath(ContinuousBenchmarks.STAN_MODELS_DIR, "binormal.model.jl"))

@tbenchmark(HMC(20, 0.5, 5), binormal, ())
bench_res = @tbenchmark(HMC(2000, 0.5, 5), binormal, ())
chn = sample(binormal(), HMC(2000,0.5,5))

describe(chn)

# turing_time = sum(chn[:elapsed])

logd = Dict(
  "name" => "Binormal (sampling from the prior)",
  "engine" => bench_res[1],
  "time" => bench_res[2],
  "mem" => bench_res[3],
  "turing" => Dict("y[1]" => mean(bench_res[4][:y])[1], "y[2]" => mean(bench_res[4][:y])[2])
)

logd = build_logd("Binormal: sampling from the prior", bench_res...)

#include(joinpath(ContinuousBenchmarks.BENCH_DIR, "binormal-stan.run.jl"))

# logd["stan"] = Dict("s" => mean(s_stan), "m" => mean(m_stan))
#logd["time_stan"] = get_stan_time("binormal")

# println("CmdStan time  : $stan_time")
# println("Turing time: $turing_time")

print_log(logd)
send_log(logd)


text:
#### TurbConvs

export TurbConv
export get_ϕ_ψ

struct TurbConv{G, SVQ, SVT, TD, OD}
    grid::G
    q::SVQ
    q_new::SVQ
    tendencies::SVQ
    aux::SVT
    tri_diag::TD
    aux_O2::Dict
    output_dir::OD
end

"""
    UniformGridParams{FT,IT}
Grid parameters used to initialize a grid.
"""
Base.@kwdef struct UniformGridParams{FT, IT}
    z_min::FT
    z_max::FT
    n_elems::IT
end

import ..FiniteDifferenceGrids: UniformGrid

UniformGrid(gp::UniformGridParams) = UniformGrid(gp.z_min, gp.z_max, gp.n_elems)

"""
    get_ϕ_ψ(ϕ)

Convenience function to get individual
variables names from co-variance.
"""
function get_ϕ_ψ(ϕ)
    if ϕ == :cv_q_tot
        _ϕ = :q_tot
        _ψ = :q_tot
    elseif ϕ == :cv_θ_liq
        _ϕ = :θ_liq
        _ψ = :θ_liq
    elseif ϕ == :cv_θ_liq_q_tot
        _ϕ = :q_tot
        _ψ = :θ_liq
    end
    return _ϕ, _ψ
end

function TurbConv(params, case::Case, output_dir::String)
    @unpack N_subdomains = params
    n_ud = N_subdomains - 2

    grid = UniformGrid(params[:UniformGridParams])

    #! format: off

    domain_set = DomainSet(gm = 1, en = 1, ud = n_ud)

    unkowns = (
        (:a             , DomainSubSet(gm = true, en = true, ud = true)),
        (:w             , DomainSubSet(gm = true, en = true, ud = true)),
        (:q_tot         , DomainSubSet(gm = true, en = true, ud = true)),
        (:θ_liq         , DomainSubSet(gm = true, en = true, ud = true)),
        (:tke           , DomainSubSet(en = true)),
        (:u             , DomainSubSet(gm = true, en = true, ud = true)),
        (:v             , DomainSubSet(gm = true, en = true, ud = true)),
    )

    aux_vars = (
        (:ρ_0           , DomainSubSet(gm = true)),
        (:p_0           , DomainSubSet(gm = true)),
        (:T             , DomainSubSet(gm = true, en = true, ud = true)),
        (:q_liq         , DomainSubSet(gm = true, en = true, ud = true)),
        (:HVSD_a        , DomainSubSet(ud = true)),
        (:HVSD_w        , DomainSubSet(ud = true)),
        (:buoy          , DomainSubSet(gm = true, en = true, ud = true)),
        (:nh_press      , DomainSubSet(ud = true)),
        (:∇buoyancy     , DomainSubSet(gm = true)),
        (:δ_model       , DomainSubSet(ud = true)),
        (:ε_model       , DomainSubSet(ud = true)),
        (:l_mix         , DomainSubSet(gm = true)),
        (:K_m           , DomainSubSet(gm = true)),
        (:K_h           , DomainSubSet(gm = true)),
        (:α_0           , DomainSubSet(gm = true)),
        (:q_tot_dry     , DomainSubSet(gm = true)),
        (:θ_dry         , DomainSubSet(gm = true)),
        (:subsidence    , DomainSubSet(gm = true)),
        (:ug            , DomainSubSet(gm = true)),
        (:vg            , DomainSubSet(gm = true)),
        (:dTdt          , DomainSubSet(gm = true)),
        (:dqtdt         , DomainSubSet(gm = true)),
        (:t_cloudy      , DomainSubSet(gm = true)),
        (:q_vap_cloudy  , DomainSubSet(gm = true)),
        (:q_tot_cloudy  , DomainSubSet(gm = true)),
        (:θ_cloudy      , DomainSubSet(gm = true)),
        (:CF            , DomainSubSet(gm = true)),
        (:mf_θ_liq      , DomainSubSet(gm = true)),
        (:mf_q_tot      , DomainSubSet(gm = true)),
        (:mf_tend_θ_liq , DomainSubSet(gm = true)),
        (:mf_tend_q_tot , DomainSubSet(gm = true)),
        (:mf_aux        , DomainSubSet(ud = true)),
        (:θ_ρ           , DomainSubSet(gm = true)),
        (:u_forcing     , DomainSubSet(gm = true)),
        (:v_forcing     , DomainSubSet(gm = true)),
        (:θ_liq_forcing , DomainSubSet(gm = true)),
        (:q_tot_forcing , DomainSubSet(gm = true)),
    )

    tri_diag_vars = (
        (:a             , DomainSubSet(gm = true)),
        (:b             , DomainSubSet(gm = true)),
        (:c             , DomainSubSet(gm = true)),
        (:f             , DomainSubSet(gm = true)),
        (:ρaK           , DomainSubSet(gm = true)),
    )

    q_2MO_vars = (
        (:entr_gain     , DomainSubSet(gm = true)),
        (:buoy          , DomainSubSet(gm = true)),
        (:press         , DomainSubSet(gm = true)),
        (:shear         , DomainSubSet(gm = true)),
        (:interdomain   , DomainSubSet(gm = true)),
    )
    #! format: on

    q = StateVec(unkowns, grid, domain_set)
    aux = StateVec(aux_vars, grid, domain_set)
    tri_diag = Dict{Symbol, StateVec}()
    tri_diag[:tke] = StateVec(tri_diag_vars, grid, domain_set)
    tri_diag[:q_tot] = StateVec(tri_diag_vars, grid, domain_set)
    tri_diag[:θ_liq] = StateVec(tri_diag_vars, grid, domain_set)
    q_new = deepcopy(q)
    tendencies = deepcopy(q)
    aux_O2 = Dict()
    aux_O2[:tke] = StateVec(q_2MO_vars, grid, domain_set)

    turb_conv =
        TurbConv(grid, q, q_new, tendencies, aux, tri_diag, aux_O2, output_dir)
end


text:
abstract type BaseObj end

abstract type LeafObj <: BaseObj end

abstract type NodeObj <: BaseObj end

abstract type BaseWorld <: NodeObj end

mutable struct NoNode <: NodeObj
end

get_id(obj::BaseObj) = obj.id

function set_parent(child::BaseObj, parent::NodeObj = NoNode())
	if !isa(child.parent, NoNode)
		removing_child(parent, child)
		@assert child.parent.children[get_id(child)] == child
		delete!(child.parent.children, get_id(child))
	end
	child.parent = parent
	if !isa(parent, NoNode)
		@assert !haskey(parent.children, get_id(child))
		parent.children[get_id(child)] = child
		added_child(parent, child)
	end
end

added_child(parent::NodeObj, child::BaseObj) = nothing
removing_child(parent::NodeObj, child::BaseObj) = nothing

function top(obj::BaseObj)
	t = obj
	while !isa(t.parent, NoNode)
		t = t.parent
	end
	t
end

function init_children(engine::Engine, parent::NodeObj, def::Dict{Symbol, Any})
	if haskey(def, :children)
		children = def[:children]
		for c in children
			child = load_def(engine, c)
			set_parent(child, parent)
		end
	end
	for (id, child) in parent.children
		inited(child)
	end
	inited(parent)
	parent
end

inited(obj::BaseObj) = nothing

for_children(f::Function, leaf::LeafObj, preorder::Bool = false) = f(leaf)
function for_children(f::Function, node::NodeObj, preorder::Bool = false)
	preorder && f(node)
	for (id, child) in node.children
		for_children(f, child, preorder)
	end
	!preorder && f(node)
end

render(leaf::LeafObj, engine::Engine) = nothing
get_local_bound(obj::LeafObj) = Empty{Float32}()

function get_local_transform(node::NodeObj)
	!haskey(node.children, :spatial) && return Matrix{Float32}(I, 4, 4)
	get_local_transform(node.children[:spatial])
end

function get_local_bound(node::NodeObj)
	!haskey(obj.children, :spatial) && return Empty{Float32}()
	get_local_bound(obj.children[:spatial])
end

function calc_local_bound(node::NodeObj, leafOnly::Bool = true)
	bound = Empty{Float32}()
	for (id, c) in node.children
		leafOnly && !isa(c, LeafObj) && continue
		localBound = get_local_bound(c)
		isa(localBound, Empty) && continue
		if isa(bound, Empty)
			bound = localBound
		else
			union!(bound, localBound)
		end
	end
	bound
end

function get_world_transform(node::NodeObj)
	!haskey(node.children, :spatial) && return Matrix{Float32}(I, 4, 4)
	get_world_transform(node.children[:spatial])
end

function get_world_bound(obj::NodeObj)
	!haskey(obj.children, :spatial) && return Empty{Float32}()
	get_world_bound(obj.children[:spatial])
end

import .Octree.getbound
getbound(obj::NodeObj) = convert(AABB{Float32}, get_world_bound(obj))

mutable struct Object <: NodeObj
	id::Symbol
	parent::NodeObj
	children::Dict{Symbol, BaseObj}

	Object(id::Symbol) = new(id, NoNode(), Dict{Symbol, BaseObj}())
end

function init(engine::Engine, ::Type{Object}, def::Dict{Symbol, Any})::Object
	obj = Object(get_id!(def))
	init_children(engine, obj, def)
	inited(obj)
	obj
end

function inited(obj::Object)
	if haskey(obj.children, :spatial)
		bound = calc_local_bound(obj)
		set_local_bound(obj.children[:spatial], bound)
	end
end

function added_child(parent::NodeObj, obj::Object)
	world = top(parent)
	!isa(world, BaseWorld) && return
	for_children(obj, true) do child
		add_to_world(world, child)
	end
end

function removing_child(parent::NodeObj, obj::Object)
	world = top(parent)
	!isa(world, BaseWorld) && return
	for_children(obj, false) do child
		remove_from_world(world, child)
	end
end

function render(obj::Object, engine::Engine)
	for (id, child) in obj.children
		if isa(child, LeafObj)
			render(child, engine)
		end
	end
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

try
    # Download a copy of https://www.ngdc.noaa.gov/geomag/data/poles/NP.xy in order to fix
    # issue #6 (https://github.com/JuliaAstro/AstroLib.jl/issues/6).
    download("https://bintray.com/giordano/AstroLib.jl-Data/download_file?file_path=NP.xy-2017-03-08",
             "NP.xy")
catch
    @warn("""Could not download file `NP.xy', you will not be able to use
"geo2mag" and "mag2geo" functions, but you can use all other routines.""")
end


text:
module TestMisc

using Test
using ThreadPools

include("util.jl")

@testset "miscellaneous" begin

    @testset "issue25" begin
        a = []
        @bthreads for x in 1:3
            if x == 2
                continue
            end
            push!(a, x)
        end
        @test sort(a) == [1,3]

        a = []
        @qthreads for x in 1:3
            if x == 2
                continue
            end
            push!(a, x)
        end
        @test sort(a) == [1,3]

        a = []
        @qbthreads for x in 1:3
            if x == 2
                continue
            end
            push!(a, x)
        end
        @test sort(a) == [1,3]
    end

end

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -14.06559781669544 <= q <= 13.37751741699121)

                     add_NL_constraint(m, :(swish(-0.7556338259975872 + -0.2850076752045205*swish(0.7432962759699406 + 0.616365769922441*swish(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + -0.9563519707270385*swish(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.21036883565748976*swish(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + -0.5819748667402531*swish(0.6367931081099356 + 0.08862555157308405*swish(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.5030828808953012*swish(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + -0.4316150506981855*swish(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + 0.4871037750371179*swish(0.12597161846431115 + 0.10770543162761159*swish(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.46799585735795457*swish(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.45029698767672954*swish(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5])))) + swish(0.9186311390475046 + -0.8559712370251158*swish(0.7432962759699406 + 0.616365769922441*swish(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + -0.9563519707270385*swish(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.21036883565748976*swish(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + -0.2512416329037399*swish(0.6367931081099356 + 0.08862555157308405*swish(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.5030828808953012*swish(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + -0.4316150506981855*swish(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + 0.7592758779705528*swish(0.12597161846431115 + 0.10770543162761159*swish(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.46799585735795457*swish(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.45029698767672954*swish(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5])))) + swish(-0.7151630404444922 + -0.15417874070690862*swish(0.7432962759699406 + 0.616365769922441*swish(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + -0.9563519707270385*swish(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.21036883565748976*swish(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + -0.016895941959414795*swish(0.6367931081099356 + 0.08862555157308405*swish(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.5030828808953012*swish(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + -0.4316150506981855*swish(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + -0.17047547651198514*swish(0.12597161846431115 + 0.10770543162761159*swish(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.46799585735795457*swish(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.45029698767672954*swish(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions