text:
using Lora

println("    Testing basic EmpiricalMCTuner constructors...")

mctuners = [EmpiricalMCTuner(0.85),
  EmpiricalMCTuner(0.85, period=50),
  EmpiricalMCTuner(0.85, maxnsteps=100),
  EmpiricalMCTuner(0.85, targetlen=1.25),
  EmpiricalMCTuner(0.85, verbose=true)]

println("    Testing that EmpiricalMCTuner tuners works with all samplers...")

npars = 3
mcmodel = model(v->-dot(v,v), grad=v->-2v, tensor=v->-2*ones(npars), dtensor=v->zeros(npars, npars, npars),
  init=ones(npars))

for mctuner in mctuners
  mcsamplers = [HMC(0.75), MALA(0.75)]

  for mcsampler in mcsamplers
    mcchain = run(mcmodel, mcsampler, SerialMC(nsteps=5000, burnin=1000), mctuner)
  end
end


text:
"""

    OneMol{FT<:AbstractFloat}

Represent properties of a single molecule, with or without phonons, as an abstract type.

See also: [`OneMol_WithPhonons{FT}`](@ref), [`OneMol_NoPhonons{FT}`](@ref)

"""
abstract type OneMol{FT<:AbstractFloat} end;

"""

    OneMol_WithPhonons{FT} <: OneMol{FT}

Represent the properties of a single molecule of ``N`` atoms, with phonons.

`numatoms` is an integer equal to the number of atoms ``N``.

`Qe`, `Me`, `Be`, `Ke`, `MI`, and `BI` are all `numatoms`-element real vectors, representing
their respective RMB oscillator parameters per atom.

`KI` is a `3*numatoms*numblocks1`-by-`3*numatoms*numblocks2` real matrix, representing the
real-space representation of (usually anisotropic) internuclear couplings among atoms
including couplings across unit cells in periodic systems (if appropriate).

`KIBloch` is a `3*numatoms`-by-`3*numatoms` complex matrix, representing the Bloch-space
representation of (usually anisotropic) internuclear couplings among atoms including
couplings across unit cells in periodic systems (if appropriate); if the system is not
periodic or this molecule has no couplings across unit cells, `KIBloch` is numerically the
same as `KI`.

`molpos0` is a 3-by-`numatoms` matrix where `molpos0[:, p]` represents the coordinates of
the atom labeled ``p``.

`molCOM` is a 3-element real vector containing the unweighted molecular center of mass.

`GvacGGinf` is a `3*numatoms`-by-`3*numatoms` complex matrix representing the Green's
function interaction matrix of the molecule in isolation in vacuum (allocated but
uninitialized until a frequency and wavevector are specified).

`alpha0` is a `numatoms`-element real vector of the atomic polarizabilities (allocated but
uninitialized until a frequency and wavevector are specified).

`Ye0ee` is an `numatoms`-element complex vector of the purely electronic contributions to
the response, useful mainly for vdW computations as a preconditioner (allocated but
uninitialized until a frequency and wavevector are specified).

`Ye0II` is an `numatoms`-element complex vector of most of the nuclear contributions to
the response, useful mainly for vdW computations as a preconditioner (allocated but
uninitialized until a frequency and wavevector are specified).

`numblocks1` & `numblocks2` are the integer numbers of neighboring unit cells that have
nontrivial internuclear couplings to the central unit cell (each is 1 for molecules that are
not periodic in that direction or have no internuclear couplings/covalent bonds with
neighboring unit cells in that direction).

`blocklist1` & `blocklist2` are range vectors indexing each unit cell block (each has the
single element 0 for molecules that are not periodic in that direction or have no
internuclear couplings/covalent bonds with neighboring unit cells in that direction).

See also: [`OneMol{FT}`](@ref), [`OneMol_NoPhonons{FT}`](@ref)

"""
struct OneMol_WithPhonons{FT} <: OneMol{FT}

    Qe::Array{FT, 1}
    Me::Array{FT, 1}
    Be::Array{FT, 1}
    Ke::Array{FT, 1}
    MI::Array{FT, 1}
    BI::Array{FT, 1}
    KI::Array{FT, 2}
    KIBloch::Array{Complex{FT}, 2}

    atompos0::Array{FT, 2}
    numatoms::Integer

    ## alphainv::Array{Complex{FT}, 2}
    GvacGGinf::SharedArray{Complex{FT}, 2}
    GvacGGinfdiag::Array{Complex{FT}, 1}
    alpha0::Array{FT, 1}
    Ye0ee::Array{Complex{FT}, 1}
    Ye0II::Array{Complex{FT}, 1}

    numblocks1::Integer
    blocklist1::UnitRange{<:Integer}
    numblocks2::Integer
    blocklist2::UnitRange{<:Integer}

end

## TO DELETE: OneMol_NoPhonons{FT} <: OneMol{FT} because it no longer
## makes sense for phonon conduction (if phonons are to be ignored,
## use the usual RMB code)
"""

    OneMol_NoPhonons{FT} <: OneMol{FT}

Represent the properties of a single molecule of ``N`` atoms, without phonons.

`numatoms` is an integer equal to the number of atoms ``N``.

`Qe`, `Me`, `Be`, and `Ke` are all `numatoms`-element real vectors.

`molpos0` is a 3-by-`numatoms` matrix where `molpos0[:, p]` represents the coordinates of
atom ``p``.

`molCOM` is a 3-element real vector containing the unweighted molecular center of mass.

`alphainv` is a `numatoms`-element complex vector representing the local inverse
susceptibility of each atom (allocated but uninitialized until a frequency and wavevector
are specified).

`GvacGGinf` is a `3*numatoms`-by-`3*numatoms` complex matrix representing the Green's
function interaction matrix of the molecule in isolation in vacuum (allocated but
uninitialized until a frequency and wavevector are specified).

`alphae` is `1 ./ alphainv`, useful mainly for vdW computations (allocated but uninitialized
until a frequency and wavevector are specified).

`alpha0` is `abs.(alphae)`, representing real atomic polarizabilities (allocated but
uninitialized until a frequency and wavevector are specified).

See also: [`OneMol{FT}`](@ref), [`OneMol_WithPhonons{FT}`](@ref)

"""

struct OneMol_NoPhonons{FT} <: OneMol{FT}

    Qe::Array{FT, 1}
    Me::Array{FT, 1}
    Be::Array{FT, 1}
    Ke::Array{FT, 1}

    atompos0::Array{FT, 2}
    numatoms::Integer

    alphainv::Array{Complex{FT}, 1}
    GvacGGinf::SharedArray{Complex{FT}, 2}
    GvacGGinfdiag::Array{Complex{FT}, 1}
    alpha0::Array{FT, 1}
    alphae::Array{Complex{FT}, 1}

end

"""

    MolSystem{T<:OneMol{<:AbstractFloat}}

Represent the properties of a collection of molecules.

`uniqueMolDataArray` contains data for each unique molecule.

`nummols` is the total number of molecules (unique or not).

`numatomslist` is the list of atoms in every molecule (unique or not).

`cumnumatomslist` is the cumulative sum of the number of atoms in every molecule, to keep
track of overall indexing.

`dupemolinds` is the list of indices tracking duplicate molecules: -1 indicates it is the
first instance of that molecule, while positive index `nn` (in the range `1:nummols`) points
to `nn` as the molecule to be duplicated.

`uniquemolinds` is complementary to `dupemolinds`, listing the index of each unique molecule
in the set of unique molecules, or -1 if it is a duplicate molecule.

See also: [`OneMol{FT}`](@ref), [`OneMol_WithPhonons{FT}`](@ref),
[`OneMol_NoPhonons{FT}`](@ref)

"""
struct MolSystem{T<:OneMol{<:AbstractFloat}}

    uniqueMolDataArray::Array{T, 1}
    molLabels::Array{<:AbstractString, 1}

    nummols::Integer
    numatomslist::Array{<:Integer, 1}
    cumnumatomslist::Array{<:Integer, 1}
    dupemolinds::Array{<:Integer, 1}
    uniquemolinds::Array{<:Integer, 1}

end

function extractKIconst!(KIconst::Array{<:AbstractFloat, 2},
                         molstrarray::Array{<:AbstractString, 1},
                         mollistfilename::AbstractString)

    if (size(KIconst, 1) != 3 || size(KIconst, 2) != 3)
        error("extractKIconst!(): KIconst must be a 3-by-3 matrix");
    end

    fill!(KIconst, zero(eltype(KIconst)));
    curridx = findfirst(map(u -> startswith(u, "KIconst "), molstrarray));
    if (curridx != nothing)
        currKIconst = map(u -> parse(eltype(KIconst), u),
                          split(replace(molstrarray[curridx], "KIconst " => "")));
        currKIconstlen = length(currKIconst);
        if (currKIconstlen == 1)
            for ii=1:3
                KIconst[ii, ii] = abs(currKIconst[1]);
            end
        elseif (currKIconstlen == 3)
            for ii=1:3
                KIconst[ii, ii] = abs(currKIconst[ii]);
            end
        elseif (currKIconstlen == 6 || currKIconstlen == 9)
            tempKIconst = zero(KIconst);
            if (currKIconstlen == 6)

                tempKIconst[1, 1] = currKIconst[1];
                tempKIconst[2, 1] = currKIconst[2];
                tempKIconst[3, 1] = currKIconst[3];
                tempKIconst[1, 2] = currKIconst[2];
                tempKIconst[2, 2] = currKIconst[4];
                tempKIconst[3, 2] = currKIconst[5];
                tempKIconst[1, 3] = currKIconst[3];
                tempKIconst[2, 3] = currKIconst[5];
                tempKIconst[3, 3] = currKIconst[6];

            else
                for ii=1:9
                    tempKIconst[ii] = currKIconst[ii];
                end
            end
            curreigvals, curreigvecs = eigen(Symmetric(tempKIconst));
            KIconst .= curreigvecs * broadcast(*, transpose(curreigvecs),
                                               abs(curreigvals));
        else
            error(string(mollistfilename, ": expected scalar (1)",
                         " or tensor (3, 6, or 9) elements for",
                         " KIconst, got something else"));
        end
    else
        println(string(mollistfilename, ": should optionally",
                       " contain a line of the form: KIconst",
                       " \${KIconstvals...}"));
    end

end

function extractRequiredMolDataVec(datalabel::AbstractString,
                                   molstrarray::Array{<:AbstractString, 1},
                                   mollistfilename::AbstractString, T::DataType=Float64)

    curridx = findfirst(map(u -> startswith(u, string(datalabel, " ")), molstrarray));
    if (curridx == nothing)
        error(string(mollistfilename, ": must contain a line",
                     " of the form: ", datalabel, " \${", datalabel, "filename}"));
    end
    return vec(collect(readdlm(replace(molstrarray[curridx],
               string(datalabel, " ") => ""), ',', T)));

end

function extractOptionalMolDataVec(datalabel::AbstractString, numatoms::Integer,
                                   molstrarray::Array{<:AbstractString, 1},
                                   mollistfilename::AbstractString, T::DataType=Float64)

    datavec = zeros(T, numatoms);
    curridx = findfirst(map(u -> startswith(u, string(datalabel, " ")), molstrarray));
    if (curridx != nothing)
        datavec .= vec(collect(readdlm(replace(molstrarray[curridx],
                      string(datalabel, " ") => ""), ',', T)))
    else
        println(string(mollistfilename, ": should optionally, as",
                       " recommended, contain a line of the form: ", datalabel,
                       " \${", datalabel, "filename}"));
    end
    return datavec;

end

function extractRequiredMolDataMat(datalabel::AbstractString,
                                   molstrarray::Array{<:AbstractString, 1},
                                   mollistfilename::AbstractString, T::DataType=Float64)

    curridx = findfirst(map(u -> startswith(u, string(datalabel, " ")), molstrarray));
    if (curridx == nothing)
        error(string(mollistfilename, ": must contain a line",
                     " of the form: ", datalabel, " \${", datalabel, "filename}"));
    end
    return readdlm(replace(molstrarray[curridx],
                   string(datalabel, " ") => ""), ',', T);

end

function extractBescale(molstrarray::Array{<:AbstractString, 1},
                        mollistfilename::AbstractString, T::DataType=Float64)

    curridx = findfirst(map(u -> startswith(u, "Bescale "), molstrarray));
    Bescale = one(T);
    if (curridx != nothing)
        Bescale = parse(T, replace(molstrarray[curridx], "Bescale " => ""));
    else
        println(string(mollistfilename, ": should optionally",
                       " contain a line of the form: Bescale \${Bescaleval}"));
        println("(default: Bescale is set to 1.0)");
    end
    return Bescale;

end

function extractKIconstatomchoices(molstrarray::Array{<:AbstractString, 1},
                                   mollistfilename::AbstractString, numatoms::Integer)

    curridx1 = findfirst(map(u -> startswith(u, "KIconstatom1 "), molstrarray));
    curridx2 = findfirst(map(u -> startswith(u, "KIconstatom2 "), molstrarray));
    atom1 = 1;
    atom2 = numatoms;

    if (curridx1 != nothing)
        atom1 = parse(Int, replace(molstrarray[curridx1], "KIconstatom1 " => ""));
    else
        println(string(mollistfilename, ": should optionally",
                       " contain a line of the form: KIconstatom1 \${KIconstatom1}"));
        println("(default: KIconstatom1 is set to 1)");
    end
    if (curridx2 != nothing)
        atom2 = parse(Int, replace(molstrarray[curridx2], "KIconstatom2 " => ""));
    else
        println(string(mollistfilename, ": should optionally",
                       " contain a line of the form: KIconstatom2 \${KIconstatom2}"));
        println("(default: KIconstatom2 is set to numatoms)");
    end

    if (atom1 < 1 || atom1 > numatoms)
        atom1 = 1;
        println(string(mollistfilename, ": required range is ",
                       "1 <= KIconstatom1 < KIconstatom2 <= numatoms"));
        println("(KIconstatom1 has been reset to 1)");
    end
    if (atom2 < 1 || atom2 > numatoms)
        atom2 = numatoms;
        println(string(mollistfilename, ": required range is ",
                       "1 <= KIconstatom1 < KIconstatom2 <= numatoms"));
        println("(KIconstatom2 has been reset to numatoms)");
    end
    if (atom2 < atom1)
        tempatom = atom2;
        atom2 = atom1;
        atom1 = tempatom;
        println(string(mollistfilename, ": recommended order is ",
                       "1 <= KIconstatom1 < KIconstatom2 <= numatoms"));
        println("(KIconstatom1 and KIconstatom2 have been switched)");
    end
    if (atom1 == atom2)
      println(string(mollistfilename, ": required order being ",
                     "1 <= KIconstatom1 < KIconstatom2 <= numatoms ",
                     "means KIconstatom1 and KIconstatom2 should be distinct"));
        if (atom2 == 1)
            atom2 = numatoms;
            println("(KIconstatom2 has been reset to numatoms)");
        else
            atom1 = 1;
            println("(KIconstatom1 has been reset to 1)");
        end
    end

    return atom1, atom2;

end

function extractGeneralMolData(molstrarray::Array{<:AbstractString, 1},
                               mollistfilename::AbstractString, T::DataType=Float64)

    Qe = extractRequiredMolDataVec("Qe", molstrarray, mollistfilename, T);
    Me = extractRequiredMolDataVec("Me", molstrarray, mollistfilename, T);
    Ke = extractRequiredMolDataVec("Ke", molstrarray, mollistfilename, T);
    numatoms = length(Qe);
    Be = extractOptionalMolDataVec("Be", numatoms, molstrarray, mollistfilename, T);
    Bescale = extractBescale(molstrarray, mollistfilename, T);
    lmul!(Bescale, Be);
    xyz = extractRequiredMolDataMat("xyz", molstrarray, mollistfilename, T);
    return numatoms, Qe, Me, Be, Ke, xyz;

end

function addKIconsttoKIselectedatoms!(KI::Array{FT, 2}, KIconst::Array{FT, 2},
                                      numatoms::Integer, atom1::Integer,
                                      atom2::Integer) where FT<:AbstractFloat

    if (size(KIconst, 1) != 3 || size(KIconst, 2) != 3)
        error("KIconst must be a 3-by-3 array");
    end
    if (rem(size(KI, 1), 3*numatoms) != 0 || rem(size(KI, 2), 3*numatoms) != 0)
        error("KI must have each dimension be an integer multiple of 3N,",
              " where N is numatoms");
    end

    numblocks1 = div(size(KI, 1), 3*numatoms);
    numblocks2 = div(size(KI, 2), 3*numatoms);

    if (iseven(numblocks1) || iseven(numblocks2))
        error("KI must have each dimension be an integer multiple of 3N,",
              " where N is numatoms, and where each multiple must be odd");
    end

    blocklist1 = div(1 - numblocks1, 2):div(numblocks1 - 1, 2);
    blocklist2 = div(1 - numblocks2, 2):div(numblocks2 - 1, 2);
    if (numblocks1 == 1 && numblocks2 == 1)
        KI[3*atom1-2:3*atom1, 3*atom1-2:3*atom1] .+= KIconst[:, :];
        KI[3*atom2-2:3*atom2, 3*atom2-2:3*atom2] .+= KIconst[:, :];
    else
        KI[3*numatoms*div(numblocks1 - 1, 2)+3*atom1-2:3*numatoms*div(numblocks1 - 1,
                                                                      2)+3*atom1,
           3*numatoms*div(numblocks2 - 1, 2)+3*atom1-2:3*numatoms*div(numblocks2 - 1,
                                                                      2)+3*atom1] .+=
           KIconst[:, :];
    end

    return numblocks1, numblocks2, blocklist1, blocklist2;

end

"""

    read1MolDataWithPhonons(molstrarray::Array{<:AbstractString, 1},
                            mollistfilename::AbstractString, T::DataType=Float64)

Read a section of a configuration file `molstrarray` (an array of strings) corresponding to
molecular data and return the data corresponding to that molecule, with phonons, of type
`OneMol_WithPhonons{T}`. Create data types to use `T` (satisfying `T<:AbstractFloat`) for
storing floating-point numbers, and output any errors with reference to the configuration file
`mollistfilename` (a string).

See also: [`OneMol{FT}`](@ref), [`OneMol_WithPhonons{FT}`](@ref)

"""
function read1MolDataWithPhonons(molstrarray::Array{<:AbstractString, 1},
                                 mollistfilename::AbstractString, T::DataType=Float64)

    if (!(T <: AbstractFloat) || isabstracttype(T))
        error("read1MolDataWithPhonons(): T must be a concrete subtype of AbstractFloat");
    end

    KIconst = zeros(T, 3, 3);
    extractKIconst!(KIconst, molstrarray, mollistfilename);

    numatoms, Qe, Me, Be, Ke, molpos0 = extractGeneralMolData(molstrarray,
                                                              mollistfilename, T);

    MI = extractRequiredMolDataVec("MI", molstrarray, mollistfilename, T);
    BI = extractOptionalMolDataVec("BI", numatoms, molstrarray, mollistfilename, T);

    atom1, atom2 = extractKIconstatomchoices(molstrarray, mollistfilename, numatoms);

    KI = extractRequiredMolDataMat("KI", molstrarray, mollistfilename, T);

    numblocks1, numblocks2, blocklist1, blocklist2 =
      addKIconsttoKIselectedatoms!(KI, KIconst, numatoms, atom1, atom2);

    GvacGG = SharedArray{Complex{T}, 2}(3*numatoms, 3*numatoms);
    fill!(GvacGG, zero(eltype(GvacGG)));

    return OneMol_WithPhonons{T}(Qe, Me, Be, Ke, MI, BI, KI,
                                 zeros(Complex{T}, 3*numatoms, 3*numatoms),
                                 molpos0, numatoms, GvacGG,
                                 zeros(Complex{T}, numatoms), zeros(T, numatoms),
                                 zeros(Complex{T}, numatoms),
                                 zeros(Complex{T}, numatoms),
                                 numblocks1, blocklist1,
                                 numblocks2, blocklist2);
end

## TO DELETE: read1MolDataNoPhonons because it no longer makes sense
## for phonon conduction (if phonons are to be ignored, use the usual
## RMB code)
"""

    read1MolDataNoPhonons(molstrarray::Array{<:AbstractString, 1},
                          mollistfilename::AbstractString, T::DataType=Float64)

Read a section of a configuration file `molstrarray` (an array of strings) corresponding to
molecular data and return the data corresponding to that molecule, without phonons, of type
`OneMol_NoPhonons{T}`. Create data types to use `T` (satisfying `T<:AbstractFloat`) for
storing floating-point numbers, and output any errors with reference to the configuration
file `mollistfilename` (a string).

See also: [`OneMol{FT}`](@ref), [`OneMol_NoPhonons{FT}`](@ref)

"""
function read1MolDataNoPhonons(molstrarray::Array{<:AbstractString, 1},
                               mollistfilename::AbstractString, T::DataType=Float64)

    if (!(T<:AbstractFloat) || isabstracttype(T))
        error("read1MolDataWithPhonons(): T must be a concrete subtype of AbstractFloat");
    end

    numatoms, Qe, Me, Be, Ke, molpos0 = extractGeneralMolData(molstrarray,
                                                              mollistfilename, T);

    ## curridx = findfirst(map(u -> startswith(u, "Qe "), molstrarray));
    ## if (curridx == nothing)
    ##     error(string(mollistfilename, ": must contain a line",
    ##                  " of the form: Qe \${Qefilename}"));
    ## end
    ## Qe = vec(collect(readdlm(replace(molstrarray[curridx], "Qe " => ""), ',', T)));
    ## numatoms = length(Qe);

    ## curridx = findfirst(map(u -> startswith(u, "Me "), molstrarray));
    ## if (curridx == nothing)
    ##     error(string(mollistfilename, ": must contain a line",
    ##                  " of the form: Me \${Mefilename}"));
    ## end
    ## Me = vec(collect(readdlm(replace(molstrarray[curridx], "Me " => ""), ',', T)));

    ## curridx = findfirst(map(u -> startswith(u, "Be "), molstrarray));
    ## Be = zeros(T, numatoms);
    ## if (curridx != nothing)
    ##     Be = vec(collect(readdlm(replace(molstrarray[curridx], "Be " => ""), ',', T)));
    ## else
    ##     println(string(mollistfilename, ": should optionally, as",
    ##                    " recommended, contain a line of the form: Be \${Befilename}"));
    ## end

    ## curridx = findfirst(map(u -> startswith(u, "Ke "), molstrarray));
    ## if (curridx == nothing)
    ##     error(string(mollistfilename, ": must contain a line",
    ##                  " of the form: Ke \${Kefilename}"));
    ## end
    ## Ke = vec(collect(readdlm(replace(molstrarray[curridx], "Ke " => ""), ',', T)));

    ## curridx = findfirst(map(u -> startswith(u, "xyz "), molstrarray));
    ## if (curridx == nothing)
    ##     error(string(mollistfilename, ": must contain a line",
    ##                  " of the form: xyz \${xyzfilename}"));
    ## end
    ## molpos0 = readdlm(replace(molstrarray[curridx], "xyz " => ""), ',', T);

    ## curridx = findfirst(map(u -> startswith(u, "Bescale "), molstrarray));
    ## Bescale = one(T);
    ## if (curridx != nothing)
    ##     Bescale = parse(T, replace(molstrarray[curridx], "Bescale " => ""));
    ## else
    ##     println(string(mollistfilename, ": should optionally",
    ##                    " contain a line of the form: Bescale \${Bescaleval}"));
    ## end
    ## lmul!(Bescale, Be);

    GvacGG = SharedArray{Complex{T}, 2}(3*numatoms, 3*numatoms);
    fill!(GvacGG, zero(eltype(GvacGG)));

    return OneMol_NoPhonons{T}(Qe, Me, Be, Ke, molpos0, numatoms,
                               zeros(Complex{T}, numatoms), GvacGG,
                               zeros(Complex{T}, numatoms), zeros(T, numatoms),
                               zeros(Complex{T}, numatoms));
end

"""

    readMolSystem(mollistfilename::AbstractString,
                  nophonons:Bool, T::DataType=Float64)

Read a configuration file `mollistfilename` (a string) to extract molecular data of type
either `OneMol_NoPhonons{T}`, if `nophonons` is `true`, or `OneMol_WithPhonons{T}`
otherwise, returning information about the composite system of type `MolSystem{molType{T}}`.
Ensure data types use `T` (satisfying `T<:AbstractFloat`) for storing floating-point
numbers.

See also: [`OneMol{T}`](@ref), [`OneMol_WithPhonons{T}`](@ref),
[`OneMol_NoPhonons{T}`](@ref), [`MolSystem{molType{T}}`](@ref)

"""
function readMolSystem(mollistfilename::AbstractString,
                       nophonons::Bool, T::DataType=Float64)

    if (!(T<:AbstractFloat) || isabstracttype(T))
        error("read1MolDataWithPhonons(): T must be a concrete subtype of AbstractFloat");
    end

    myreadfunc = read1MolDataNoPhonons;
    molType = OneMol_NoPhonons;
    if (!nophonons)
        molType = OneMol_WithPhonons;
        myreadfunc = read1MolDataWithPhonons;
    end

    molstrarray = strip.(chomp.(readlines(mollistfilename)));
    molstartinds = findall(map(u -> startswith(u, "MOL "), molstrarray));
    molendinds = findall(map(u -> startswith(u, "ENDMOL"), molstrarray));
    nummols = length(molendinds);
    if (length(molstartinds) != nummols)
        error(string(mollistfilename, ": number of MOL and ENDMOL statements",
                     " should be the same"));
    end
    numatomslist = zeros(Int, nummols);
    molbothinds = vec(transpose(hcat(molstartinds, molendinds)));
    if (molbothinds != sort(molbothinds))
        error(string(mollistfilename, ": MOL and ENDMOL statements",
                     " should alternate in ordered occurrence"));
    end

    molstrarraydupeinds = findall(map(u -> startswith(u, "DUPLICATE "), molstrarray));
    numdupes = length(molstrarraydupeinds);
    numunique = nummols - numdupes;
    dupemolinds = fill(-1, nummols);

    uniqueMolDataArray = Array{molType{T}, 1}(undef, numunique);
    allMolLabelArray = Array{String, 1}(undef, nummols);
    uniquemolcounter = zero(Int);
    uniquemolinds = fill(-1, nummols);
    for nn=1:nummols
        allMolLabelArray[nn] = replace(strip(chomp(molstrarray[molstartinds[nn]])),
                                       "MOL " => "");
        currmolstrarray = strip.(chomp.(molstrarray[molstartinds[nn]+1:molendinds[nn]-1]));
        if (findfirst(map(u -> startswith(u, "DUPLICATE "), currmolstrarray)) == nothing)
            uniquemolcounter += 1;
            uniqueMolDataArray[uniquemolcounter] = myreadfunc(currmolstrarray,
                                                              mollistfilename, T);
            numatomslist[nn] = uniqueMolDataArray[uniquemolcounter].numatoms;
            uniquemolinds[nn] = uniquemolcounter;
        elseif (nn > 1)
            dupemolinds[nn] = findfirst(allMolLabelArray .==
                                        replace(strip(chomp(currmolstrarray[1])),
                                                "DUPLICATE " => ""));
            numatomslist[nn] = uniqueMolDataArray[dupemolinds[nn]].numatoms;
        else
            error(string(mollistfilename, ": the first molecule listed cannot be a",
                         "duplicate"));
        end
    end
    cumnumatomslist = cumsum(numatomslist);

    return MolSystem{molType{T}}(uniqueMolDataArray, allMolLabelArray,
                                 nummols, numatomslist, cumnumatomslist,
                                 dupemolinds, uniquemolinds);

end

ext: jl


text:
import LinearAlgebra: I, inv
import ForneyLab: unsafeCov, unsafeMean, unsafePrecision, VariateType,
				  collectNaiveVariationalNodeInbounds, assembleClamp!, ultimatePartner

include("util.jl")

export ruleVariationalNARMAXOutNPPPPPP,
       ruleVariationalNARMAXIn1PNPPPPP,
       ruleVariationalNARMAXIn2PPNPPPP,
       ruleVariationalNARMAXIn3PPPNPPP,
	   ruleVariationalNARMAXIn4PPPPNPP,
	   ruleVariationalNARMAXIn5PPPPPNP,
	   ruleVariationalNARMAXIn6PPPPPPN

function ruleVariationalNARMAXOutNPPPPPP(ϕ :: Function,
									     marg_y :: Nothing,
									     marg_θ :: ProbabilityDistribution{Multivariate},
                                         marg_u :: ProbabilityDistribution{Univariate},
										 marg_x :: ProbabilityDistribution{Multivariate},
										 marg_z :: ProbabilityDistribution{Multivariate},
										 marg_r :: ProbabilityDistribution{Multivariate},
										 marg_τ :: ProbabilityDistribution{Univariate})

	# Extract moments of beliefs
	mθ = unsafeMean(marg_θ)
	mu = unsafeMean(marg_u)
	mx = unsafeMean(marg_x)
	mz = unsafeMean(marg_z)
	mr = unsafeMean(marg_r)
	mτ = unsafeMean(marg_τ)

	# Basis function applied to moments
	ϕm = ϕ([mu; mx; mz; mr])

	# Set outgoing message
	return Message(Univariate, GaussianMeanPrecision, m=mθ'*ϕm, w=mτ)
end

function ruleVariationalNARMAXIn1PNPPPPP(ϕ :: Function,
									     marg_y :: ProbabilityDistribution{Univariate},
                                         marg_θ :: Nothing,
										 marg_u :: ProbabilityDistribution{Univariate},
									     marg_x :: ProbabilityDistribution{Multivariate},
										 marg_z :: ProbabilityDistribution{Multivariate},
										 marg_r :: ProbabilityDistribution{Multivariate},
										 marg_τ :: ProbabilityDistribution{Univariate})

    # Extract moments of beliefs
	my = unsafeMean(marg_y)
	mu = unsafeMean(marg_u)
	mx = unsafeMean(marg_x)
	mz = unsafeMean(marg_z)
	mr = unsafeMean(marg_r)
	mτ = unsafeMean(marg_τ)

	# Basis function applied to moments
	ϕm = ϕ([mu; mx; mz; mr])

	# Number of coefficients
	num_coeffs = size(ϕm,1)

	# Update parameters
	# Ψ = mτ*ϕm*ϕm'
	Ψ = mτ*(ϕm*ϕm' .+ 1e-8*Matrix{Float64}(I, num_coeffs, num_coeffs))
	ψ = mτ*my*ϕm

    return Message(Multivariate, GaussianWeightedMeanPrecision, xi=ψ, w=Ψ)
end

function ruleVariationalNARMAXIn2PPNPPPP(ϕ :: Function,
									  	 marg_y :: ProbabilityDistribution{Univariate},
									  	 marg_θ :: ProbabilityDistribution{Multivariate},
										 marg_u :: Nothing,
										 marg_x :: ProbabilityDistribution{Multivariate},
										 marg_z :: ProbabilityDistribution{Multivariate},
										 marg_r :: ProbabilityDistribution{Multivariate},
										 marg_τ :: ProbabilityDistribution{Univariate})

	error("Variable is assumed to be observed")
    return Nothing
end

function ruleVariationalNARMAXIn3PPPNPPP(ϕ :: Function,
									  	 marg_y :: ProbabilityDistribution{Univariate},
									  	 marg_θ :: ProbabilityDistribution{Multivariate},
									  	 marg_u :: ProbabilityDistribution{Univariate},
										 marg_x :: Nothing,
										 marg_z :: ProbabilityDistribution{Multivariate},
										 marg_r :: ProbabilityDistribution{Multivariate},
										 marg_τ :: ProbabilityDistribution{Univariate})

	error("Variable is assumed to be observed")
    return Nothing
end

function ruleVariationalNARMAXIn4PPPPNPP(ϕ :: Function,
									  	 marg_y :: ProbabilityDistribution{Univariate},
									  	 marg_θ :: ProbabilityDistribution{Multivariate},
										 marg_u :: ProbabilityDistribution{Univariate},
										 marg_x :: ProbabilityDistribution{Multivariate},
							  	      	 marg_z :: Nothing,
										 marg_r :: ProbabilityDistribution{Multivariate},
										 marg_τ :: ProbabilityDistribution{Univariate})

	error("Variable is assumed to be observed")
    return Nothing
end

function ruleVariationalNARMAXIn5PPPPPNP(ϕ :: Function,
									  	 marg_y :: ProbabilityDistribution{Univariate},
									  	 marg_θ :: ProbabilityDistribution{Multivariate},
										 marg_u :: ProbabilityDistribution{Univariate},
										 marg_x :: ProbabilityDistribution{Multivariate},
							  	      	 marg_z :: ProbabilityDistribution{Multivariate},
										 marg_r :: Nothing,
										 marg_τ :: ProbabilityDistribution{Univariate})

	error("Variable is assumed to be observed")
    return Nothing
end

function ruleVariationalNARMAXIn6PPPPPPN(ϕ :: Function,
									 	 marg_y :: ProbabilityDistribution{Univariate},
										 marg_θ :: ProbabilityDistribution{Multivariate},
										 marg_u :: ProbabilityDistribution{Univariate},
										 marg_x :: ProbabilityDistribution{Multivariate},
										 marg_z :: ProbabilityDistribution{Multivariate},
										 marg_r :: ProbabilityDistribution{Multivariate},
  										 marg_τ :: Nothing)

	# Extract moments of beliefs
	mθ,Vθ = unsafeMeanCov(marg_θ)
	my = unsafeMean(marg_y)
	mu = unsafeMean(marg_u)
	mx = unsafeMean(marg_x)
	mz = unsafeMean(marg_z)
	mr = unsafeMean(marg_r)

	# Basis expansion of moments
	ϕm = ϕ([mu; mx; mz; mr])

	# Update parameters
	α = 3/2.
	β = (my^2 - 2*my*(mθ'*ϕm) + (mθ'*ϕm)^2 + ϕm'*Vθ*ϕm)/2.

	# Set message
	return Message(Univariate, Gamma, a=α, b=β)
end

function collectNaiveVariationalNodeInbounds(node::NAutoRegressiveMovingAverageX, entry::ScheduleEntry)
	inbounds = Any[]

	# Push function to calling signature (ϕ needs to be defined in user scope)
	push!(inbounds, Dict{Symbol, Any}(:ϕ => node.ϕ, :keyword => false))

    target_to_marginal_entry = currentInferenceAlgorithm().target_to_marginal_entry

    for node_interface in entry.interface.node.interfaces
        inbound_interface = ultimatePartner(node_interface)
        if node_interface === entry.interface
            # Ignore marginal of outbound edge
            push!(inbounds, nothing)
        elseif (inbound_interface != nothing) && isa(inbound_interface.node, Clamp)
            # Hard-code marginal of constant node in schedule
            push!(inbounds, assembleClamp!(inbound_interface.node, ProbabilityDistribution))
        else
            # Collect entry from marginal schedule
            push!(inbounds, target_to_marginal_entry[node_interface.edge.variable])
        end
    end

    return inbounds
end

max_stars_repo_name: wmkouw/NARMAX

max_issues_repo_name: wmkouw/NARMAX

max_forks_repo_name: wmkouw/NARMAX


text:
@testset "Utils" begin
  @testset "compose" begin
    df1 = compose(jura, (:Cd, :Cu, :Pb, :Co, :Cr, :Ni, :Zn))
    @test names(df1) == ["X","Y","Rock","Land","coda"]
    @test size(df1) == (359, 5)
    @test df1[1,"coda"] == Composition(Cd=1.74, Cu=25.72, Pb=77.36, Co=9.32, Cr=38.32, Ni=21.32, Zn=92.56)

    df2 = compose(jura, ("Cd", "Cu", "Pb", "Co", "Cr", "Ni", "Zn"))
    @test df1 == df2

    df3 = compose(jura, (:Cd, :Cu, :Pb, :Co, :Cr, :Ni, :Zn) => :Composition)
    df4 = compose(jura, ("Cd", "Cu", "Pb", "Co", "Cr", "Ni", "Zn") => "Composition")
    @test names(df3) == names(df4) == ["X","Y","Rock","Land","Composition"]
  end

  @testset "readcoda" begin
    csv = joinpath(datadir,"jura.csv")
    df = readcoda(csv, codanames=(:Cd, :Cu, :Pb, :Co, :Cr, :Ni, :Zn))
    @test names(df) == ["X","Y","Rock","Land","coda"]
    @test size(df) == (359, 5)
    @test df[1,"coda"] == Composition(Cd=1.74, Cu=25.72, Pb=77.36, Co=9.32, Cr=38.32, Ni=21.32, Zn=92.56)

    df = readcoda(csv, codanames=(:Cd, :Cu, :Pb, :Co, :Cr, :Ni, :Zn) => :Composition)
    @test names(df) == ["X","Y","Rock","Land","Composition"]
  end
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

@doc """
Group InsertBeforeOrAfterStatement actions in decreasing order of
(basic block index, statement index) of the actions; for the same
statement, "after" is before "before".
"""
function group_action(
    insert_before_or_after_statement_actions :: Vector{InsertBeforeOrAfterStatement},
    action                                   :: InsertBeforeOrAfterStatement
)
    for i in 1 : length(insert_before_or_after_statement_actions)
        action1 = insert_before_or_after_statement_actions[i]
        if action.bb == action1.bb
            if action.stmt_idx == action1.stmt_idx
                if !action.before && action1.before
                    # Let the actions for a later statement happen earlier since
                    # they do not impact the positions of earlier statements.
                    insert!(insert_before_or_after_statement_actions, i, action)

                    # No need to try match with other elements: the vector is
                    # ensured to have only one elment for one "before" and/or
                    # one "after" for one specific statement.
                    return
                elseif action.before && !action1.before
                    continue
                else
                    # The two actions are both before or both after the statement.
                    append!(action1.new_stmts, action.new_stmts)
                    return
                end
            elseif action.stmt_idx > action1.stmt_idx
                # Let the actions for a later statement happen earlier since
                # they do not impact the positions of earlier statements.
                insert!(insert_before_or_after_statement_actions, i, action)
                return
            end
        elseif action.bb.label > action.bb.label
            insert!(insert_before_or_after_statement_actions, i, action)
            return
        end
    end
    push!(insert_before_or_after_statement_actions, action)
end

@doc """
Group InsertBeforeLoopHead actions.
"""
function group_action(
    insert_before_loop_head_actions :: Vector{InsertBeforeLoopHead},
    action                          :: InsertBeforeLoopHead
)
    for action1 in insert_before_loop_head_actions
        if action1.loop == action.loop && action1.outside_loop == action.outside_loop
            append!(action1.new_stmts, action.new_stmts)

            # No need to try match with other elements: the vector is ensured
            # to have only one elment for one specific block.
            return
        end
    end
    push!(insert_before_loop_head_actions, action)
end

@doc """
Group InsertOnEdge actions.
"""
function group_action(
    insert_on_edge_actions :: Vector{InsertOnEdge},
    action                 :: InsertOnEdge
)
    for action1 in insert_on_edge_actions
        if action1.from_bb == action.from_bb && action1.to_bb == action.to_bb
            append!(action1.new_stmts, action.new_stmts)

            # No need to try match with other elements: the vector is ensured
            # to have only one elment for one specific edge.
            return
        end
    end
    push!(insert_on_edge_actions, action)
end

@doc """
The index of the statement in the block. Return -1 if the statement is not in
the block.
"""
function index_of_statement(
    bb   :: BasicBlock,
    stmt :: Statement
)
    for idx = 1 : length(bb.statements)
        if bb.statements[idx] == stmt
            return idx
        end
    end
    return -1
end

@doc """
Perform all the actions. Return the transformed AST.
"""
function CFG_transformation(
    actions  :: Vector{Action},
    cfg      :: CFG
)
    # Implement each action. Some important details to care:
    # (1) InsertBeforeLoopHead actions must be done last, because it changes the
    # predecessor-successor relationship, which affects InsertOnEdge (Some edge
    # no longer exists after InsertBeforeLoopHead).
    # (2) Actions that work on the same loop for InsertBeforeLoopHead, or the same
    # Statement for InsertBeforeOrAfterStatement, or the same edge for InsertOnEdge,
    # should be grouped and implemented together

    # Group the actions.
    insert_before_or_after_statement_actions = Vector{InsertBeforeOrAfterStatement}()
    insert_before_loop_head_actions          = Vector{InsertBeforeLoopHead}()
    insert_on_edge_actions                   = Vector{InsertOnEdge}()
    for action in actions
        if typeof(action) == InsertBeforeOrAfterStatement
            group_action(insert_before_or_after_statement_actions, action)
        elseif typeof(action) == InsertBeforeLoopHead
            group_action(insert_before_loop_head_actions, action)
        else
            group_action(insert_on_edge_actions, action)
        end
    end

    # Implement InsertOnEdge actions
    for action in insert_on_edge_actions
        if isempty(action.new_stmts)
            continue
        end

        new_bb, new_goto_stmt = CFGs.insertBetween(cfg, action.from_bb.label, action.to_bb.label)
        new_bb.statements     = action.new_stmts
        if new_goto_stmt != nothing
            push!(new_bb.statements, new_goto_stmt)
        end
    end

    # Implement InsertBeforeOrAfterStatement actions
    for action in insert_before_or_after_statement_actions
        if isempty(action.new_stmts)
            continue
        end

        for i = 1 : length(action.new_stmts)
            if action.before
                insert!(action.bb.statements, action.stmt_idx + i - 1, action.new_stmts[i])
            else
                insert!(action.bb.statements, action.stmt_idx + i, action.new_stmts[i])
            end
        end
    end

    # Implement InsertBeforeLoopHead actions
    for action in insert_before_loop_head_actions
        if isempty(action.new_stmts)
            continue
        end

        new_bb, new_goto_stmt = CFGs.insertBefore(cfg, action.loop.head, action.outside_loop, action.loop.back_edge)
        new_bb.statements     = action.new_stmts
        if new_goto_stmt != nothing
            push!(new_bb.statements, new_goto_stmt)
        end
    end
end

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
function ackley(x, a=20, b=0.2, c=2π)
    d = length(x)
    return -a*exp(-b*sqrt(sum(x.^2)/d)) -
        exp(sum(cos.(c*xi) for xi in x)/d) + a +
        exp(1)
end


text:
## Printing
#==========================================================================================#

"""
    print_info(wing :: AbstractWing)

Print the relevant geometric characteristics of a `HalfWing` or `Wing`.
"""
function print_info(wing :: AbstractWing, head = "")
    labels = [ "Aspect Ratio", "Span (m)", "Area (m²)", "Mean Aerodynamic Chord (m)", "Mean Aerodynamic Center (m)" ]
    wing_info = properties(wing)
    data = Any[ labels wing_info ]
    header = [ head, "Value" ]
    h1 = Highlighter( (data,i,j) -> (j == 1), foreground = :cyan, bold = true)

    pretty_table(data, header, alignment = [:c, :c], highlighters = h1, vlines = :none, formatters = ft_round(8))
end

"""
    print_coefficients(nf_coeffs, ff_coeffs, name = "")

Print a pretty table of the nearfield and farfield coefficients with an optional name.
"""
function print_coefficients(nf_coeffs, ff_coeffs, name = "")
    coeffs = [ ifelse(length(nf_coeffs) == 8, ["CD", "CDv"], []); [ "CDi", "CY", "CL", "Cl", "Cm", "Cn" ] ]
    data = [ coeffs nf_coeffs [ ff_coeffs; fill("—", 3) ] ]
    head = [ name, "Nearfield", "Farfield" ]
    h1 = Highlighter( (data,i,j) -> (j == 1), foreground = :cyan, bold = true)

    pretty_table(data, head, alignment = [:c, :c, :c], highlighters = h1, vlines = :none, formatters = ft_round(8))
end

"""
    print_derivatives(nf_coeffs, ff_coeffs, name = "";
                      axes = "")

Print a pretty table of the aerodynamic coefficients and derivatives with an optional name and axis name.
"""
function print_derivatives(comp, name = ""; axes = "")
    coeffs  = ["CX", "CY", "CZ", "Cℓ", "Cm", "Cn"]
    nf_vars = [ "$name" "Values" "" "" "Derivatives" "" "" ; "" "$axes" "∂/∂α, 1/rad" "∂/∂β, 1/rad" "∂/∂p̄" "∂/∂q̄" "∂/∂r̄" ]
    nf_rows = @views [ coeffs comp[1:6,[1,3,4,5,6,7]] ]

    pretty_table(nf_rows, nf_vars, alignment = :c, header_crayon = Crayon(bold = true), subheader_crayon = Crayon(foreground = :yellow, bold = true), highlighters = Highlighter( (data,i,j) -> (j == 1), foreground = :cyan, bold = true), vlines = :none, formatters = ft_round(8))
end

"""
    print_coefficients(system :: VortexLatticeSystem, name = :aircraft;
                       components = false)

Print a pretty table of the total nearfield and farfield coefficients of a `VortexLatticeSystem` with an optional name.

A named Boolean argument `components` is provided to also enable the printing of any possible components.
"""
function print_coefficients(system :: VortexLatticeSystem, name = :aircraft; components = false)
    if components
        nf_c = nearfield_coefficients(system)
        ff_c = farfield_coefficients(system)
        [ print_coefficients(nf_c[key], ff_c[key], key) for key in keys(system.vortices) ]
        print_coefficients(nearfield(system), farfield(system), name)
    else
        print_coefficients(nearfield(system), farfield(system), name)
    end

    nothing
end


text:
module BinningAnalysis

import Statistics: mean, var
import Base: push!, append!, show, summary

include("binning.jl")
export LogBinner, push!, append!

include("statistics.jl")
export mean, var, varN, tau, std_error
export all_vars, all_varNs, all_taus, all_std_errors
export convergence, has_converged

include("cosmetics.jl")

end # module


text:
@testset "SubspaceQs" begin

    dummy_data = rand(4,10)
    subspaces = [[1,2], [3,4]]
    model = SVDD.SubRandomOCClassifier(dummy_data, subspaces)

    params = Dict{Symbol, Any}()
    params[:subspaces] = subspaces
    params[:scale_fct] = min_max_normalize
    params[:combination_fct] = max

    @testset "initialize_qs" begin
        @testset "TestPQs" begin
            p = copy(params)
            p[:x] = 5
            sub_qs = @inferred initialize_qs(SubspaceQs{TestPQs}, model, dummy_data, p)
            @test sub_qs.single_space_strategy.x == 5
            @test sub_qs.scale_fct == min_max_normalize
            @test sub_qs.combination_fct == max
            for labels in [labelmap(fill(:U, 10)), labelmap(fill(:Lin, 10)), labelmap(fill(:Lout, 10))]
                scores = qs_score(sub_qs, dummy_data, labels)
                @test length(scores) == size(dummy_data, 2)
            end
        end

        @testset "scale-and-combine" begin
            labels = labelmap(fill(:U, 10))
            p = copy(params)
            p[:x] = 5

            sub_qs = @inferred initialize_qs(SubspaceQs{TestPQs}, model, dummy_data, p)
            scores = qs_score(sub_qs, dummy_data, labels)
            @test scores ≈ collect(range(0,1,length=10))

            p[:combination_fct] = +
            sub_qs = @inferred initialize_qs(SubspaceQs{TestPQs}, model, dummy_data, p)
            scores = qs_score(sub_qs, dummy_data, labels)
            @test scores ≈ collect(range(0,1,length=10)) .* 2

            p[:scale_fct] = identity
            sub_qs = @inferred initialize_qs(SubspaceQs{TestPQs}, model, dummy_data, p)
            scores = qs_score(sub_qs, dummy_data, labels)
            @test scores ≈ collect(range(1,10,length=10)) .* 2

            scores = qs_score(sub_qs, rand(4,5), labelmap(fill(:U, 5)))
            @test length(scores) == 5

            scores = qs_score(sub_qs, rand(4,25), labelmap(fill(:U, 25)))
            @test length(scores) == 25

            @test_throws DimensionMismatch qs_score(sub_qs, rand(3, 10), labelmap(fill(:U, 10)))
        end

        labels = labelmap(fill(:U, 10))
        for q in [RandomPQs, RandomOutlierPQs]
            @testset "$q" begin
                sub_qs = @inferred initialize_qs(SubspaceQs{q}, model, dummy_data, params)
            end
        end
        @testset "BoundaryNeighborCombinationPQs" begin
            p = copy(params)
            p[:η] = 0.8
            p[:p] = 0.2
            sub_qs = @inferred initialize_qs(SubspaceQs{BoundaryNeighborCombinationPQs{typeof(model)}}, model, dummy_data, p)
            @test_throws DimensionMismatch qs_score(sub_qs, rand(4,5), labels)
            @test sub_qs.single_space_strategy.p == 0.2
            @test sub_qs.single_space_strategy.η == 0.8
            query_data = rand(4,5)
            sub_qs = @inferred initialize_qs(SubspaceQs{BoundaryNeighborCombinationPQs{typeof(model)}}, model, query_data, p)
            scores = qs_score(sub_qs, query_data, labelmap(fill(:U, 5)))
            @test length(scores) == size(query_data, 2)
        end

        for q in [HighConfidencePQs, DecisionBoundaryPQs]
            p = copy(params)
            p[:subspaces] = [[1,2], [3,4], [1,4]]
            @testset "$q" begin
                sub_qs = @inferred initialize_qs(SubspaceQs{q}, model, dummy_data, p)
                @test_throws DimensionMismatch qs_score(sub_qs, rand(3,10), labels)
                @test_throws DimensionMismatch qs_score(sub_qs, rand(5,10), labels)
                scores = qs_score(sub_qs, rand(4, 20), labelmap(fill(:U, 20)))
                @test length(scores) == 20
                scores = qs_score(sub_qs, rand(4, 5), labelmap(fill(:U, 5)))
                @test length(scores) == 5
            end
        end

        @testset "Not implemented Error" begin
            for q in [MinimumLossPQs, MinimumMarginPQs]
                p = copy(params)
                p[:p_inlier] = 0.1
                sub_qs = initialize_qs(SubspaceQs{q}, model, dummy_data, p)
                @test_throws ErrorException qs_score(sub_qs, rand(4, 10), labelmap(fill(:U, 10)))
            end

            for q in [ExpectedMinimumMarginPQs, ExpectedMaximumEntropyPQs, NeighborhoodBasedPQs]
                sub_qs = initialize_qs(SubspaceQs{q}, model, dummy_data, params)
                @test_throws ErrorException qs_score(sub_qs, rand(4, 10), labelmap(fill(:U, 10)))
            end
        end

        @testset "process_queries!" begin
            data = rand(4,10)
            subspaces = [[1,2], [3,4]]
            pools = fill(:U, size(data,2))

            train = BitArray(vcat(false, trues(8), false))
            gamma_strategy = FixedGammaStrategy(MLKernels.GaussianKernel(2))
            C_strategy = FixedCStrategy(0.1)
            init_strategy = SVDD.SimpleSubspaceStrategy(gamma_strategy, C_strategy, gamma_scope=Val(:Global))
            update_strategy = SVDD.FixedWeightStrategy(42.0, 5.0)

            global_query_ids = [2,4]
            query_labels = [:outlier, :inlier]

            @testset "FullSplit" begin
                split_strategy = DataSplits(train, .~train, FullSplitStrat())
                train_data, train_pools, _ = get_train(split_strategy, data, pools)
                model = SVDD.SubSVDD(train_data, subspaces, train_pools)
                set_param!(model, Dict(:weight_update_strategy => update_strategy))
                initialize!(model, init_strategy)
                data_updated, pools_updated, labels_updated = OneClassActiveLearning.process_queries!(global_query_ids, query_labels, model, split_strategy, data, pools, labels)
                # FullSplit -> first train index and third train index are updated
                @test all(model.v .== [5.0, 1.0, 42.0, 1.0, 1.0, 1.0, 1.0, 1.0])
            end

            @testset "UnlabeledAndLabeledInlierSplit" begin
                split_strategy = DataSplits(train, .~train, UnlabeledAndLabeledInlierSplitStrat())
                train_data, train_pools, _ = get_train(split_strategy, data, pools)
                model = SVDD.SubSVDD(train_data, subspaces, train_pools)
                set_param!(model, Dict(:weight_update_strategy => update_strategy))
                initialize!(model, init_strategy)

                data_updated, pools_updated, labels_updated = OneClassActiveLearning.process_queries!(global_query_ids, query_labels, model, split_strategy, data, pools, labels)
                # FullSplit -> first train observation is removed because it is labled :Lout, second train observation in updated data set is updated
                @test all(model.v .== [1.0, 42.0, 1.0, 1.0, 1.0, 1.0, 1.0])
            end
        end
    end
end


text:
"""
Main module of ModiaBase.

* Developers: Hilding Elmqvist, Mogram AB, Martin Otter, DLR
* First version: December 2020
* License: MIT (expat)

"""
module ModiaBase

const path    = dirname(dirname(@__FILE__))   # Absolute path of package directory
const Version = "0.7.2"
const Date    = "2021-11-23"

#println("\nImporting ModiaBase Version $Version ($Date)")

using Unitful

include("LinearIntegerEquations.jl")

include("BLTandPantelidesUtilities.jl")
using .BLTandPantelidesUtilities

include("BLTandPantelides.jl")
using .BLTandPantelides

include("Differentiate.jl")
using .Differentiate

include("Tearing.jl")

include("Simplify.jl")
using .Simplify

include("Symbolic.jl")
using .Symbolic

include("EquationAndStateInfo.jl")
include("StateSelection.jl")

end

max_stars_repo_name: ModiaSim/ModiaBase

max_issues_repo_name: ModiaSim/ModiaBase

max_forks_repo_name: ModiaSim/ModiaBase


text:
# This file contains functions for creating spatial weights from polygon data

# Auxiliar function to return the points from a vector of polygons, multipolygon, or a mix of the two
function _getpointsPoligon(P::T where T <:AbstractPolygon)::Tuple{Vector{Float64}, Vector{Float64}}
    c = coordinates(P)

    ppointsx =  Array{Float64}(undef, length(c[1]))
    ppointsy =  Array{Float64}(undef, length(c[1]))

    ppointsx = map(a -> a[1], c[1])
    ppointsy = map(a -> a[2], c[1])

    return ppointsx, ppointsy
end

function _getpointsPoligon(P::T where T <:AbstractMultiPolygon)::Tuple{Vector{Float64}, Vector{Float64}}
    c = coordinates(P)

    ppointsx =  Array{Float64}(undef, length(c[1]))
    ppointsy =  Array{Float64}(undef, length(c[1]))

    xy = reduce(vcat,
                map(1:length(c)) do i
                    x = map(a -> a[1], c[i][1])
                    y = map(a -> a[2], c[i][1])

                    hcat(x, y)
                end
            )

    ppointsx = xy[:,1]
    ppointsy = xy[:,2]

    return ppointsx, ppointsy
end

_getpointsPoligon(::Missing) = throw(ArgumentError("Missing geometry"));

# Auxiliar function to return the number of times the points of two polygons hits
function _hits(icoords::Vector{Vector{Float64}}, jcoords::Vector{Vector{Float64}}, critthr::Int64, tol::Float64)::Bool
    nhits = 0

    nicoords = length(icoords[1])
    njcoords = length(jcoords[1])

    for pti in 1:nicoords
        for ptj in 1:njcoords
            @inbounds xd = icoords[1][pti] - jcoords[1][ptj]
            if abs(xd) > tol continue end

            @inbounds yd = icoords[2][pti] - jcoords[2][ptj]
            if abs(yd) > tol continue end

            dist = hypot(xd, yd)
            if dist <= tol
                nhits = nhits + 1
            end

            # Leave if number of hits is larger than the required by the criterion
            if nhits >= critthr
                return true
            end
        end
    end

    return false
end

"""
    polyneigh(P, criterion = :Queen)
Build a spatial weights object from a vector of polygons `P`.

# Optional Arguments
- `criterion=:Queen`: neighbour criterion. `:Queen` or `:Rook`.
- `tol=0.0`: tolerance for polygon contiguity.
"""
function polyneigh(P::Vector{T} where T <:Union{Missing,AbstractPolygon,AbstractMultiPolygon}; criterion::Symbol = :Queen, tol::Float64 = 0.0)::SpatialWeights
    n = length(P)

    xy = _getpointsPoligon.(P)

    x = map(a -> a[1], xy)
    y = map(a -> a[2], xy)

    # Bounding box
    xmin, ymin = minimum.(x) .- tol, minimum.(y) .- tol
    xmax, ymax = maximum.(x) .+ tol, maximum.(y) .+ tol

    BBpols = hcat(xmin, ymin, xmax, ymax, [1:n;])

    # Select candidates if bouning box overlaps or touches
    BBpols = sortslices(BBpols, dims=1)

    candidates = copy.(fill(Int[], n))
    Threads.@threads for i in 1:n-1
        for j in (i+1):n
            if (BBpols[j,1] <= BBpols[i,3]) && # xmin[j] <= xmax[i]
                (BBpols[j,2] <= BBpols[i,4]) && # ymin[j] <= ymax[i]
                (BBpols[j,3] >= BBpols[i,1]) && # xmax[j] >= xmin[i]
                (BBpols[j,4] >= BBpols[i,2]) # ymax[j] >= ymin[i]
                    push!(candidates[Int(BBpols[i,5])], BBpols[j,5])
            end
        end
    end

    # Check neighbours
    neighs = copy.(fill(Int[], n))
    weights = copy.(fill(Float64[], n))
    nneighs = zeros(Int,n)

    if criterion == :Queen critthr = 1 end
    if criterion == :Rook  critthr = 2 end

    # Use multithreading to calculate polygons hits
    neighssingle = copy.(fill(Int[], n))
    Threads.@threads for i in 1:n
        @inbounds for j in candidates[i]
            p1 = [x[i], y[i]]
            p2 = [x[j], y[j]]
            polhits = _hits(p1, p2, critthr, tol)

            if polhits
                push!(neighssingle[i], j)
            end
        end
    end

    # Makes neighbors bilateral
    for i in 1:n
        for j in neighssingle[i]
            push!(neighs[i], j)
            push!(neighs[j], i)

            nneighs[i] += 1
            nneighs[j] += 1
        end
    end

    # Sort neighbours lists and compute weights
    for i in 1:n
        sort!(neighs[i])
        weights[i] = ones(nneighs[i]) ./ nneighs[i]
    end

    SpatialWeights(n, neighs, weights, nneighs, :row)
end

"""
    polyneigh(A, criterion = :Queen)
Build a spatial weights object from table A that contains a geometry column.

# Optional Arguments
- `criterion=:Queen`: neighbour criterion. `:Queen` or `:Rook`.
- `tol=0.0`: tolerance for polygon contiguity.
"""
function polyneigh(A::Any; criterion::Symbol = :Queen, tol::Float64 = 0.0)::SpatialWeights

    istable(A) || throw(ArgumentError("Unknown geometry or not polygon geometry"))

    (:geometry in propertynames(A)) || throw(ArgumentError("table does not have :geometry information"))

    return polyneigh(A.geometry, criterion = criterion, tol = tol)
end


text:
include("exp_set_list_all.jl");

println("Printing simulation");

rho=mono.params[:rho];
println("Domain: $(rho), m=",m);
for (i,state)=enumerate(state_list)
    if (!isnothing(state))
        err=showerr(state,output=false,n=state.params[:target_n]*2)
        graphname=state.params[:graphname];
        println("$graphname $(Float64(err))");
    else
        println("x");
    end
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using GAIO

# This example demonstrates how to get a ~2x speedup
# in your code using your CPU's SIMD capabilities.

N = 3
const σ, ρ, β = 10.0, 28.0, 0.4

# By default, GAIO is set up to accept functions of the form
function f(x)
    # Some map, here we use the Lorenz equation
    dx = (
           σ * x[2] -    σ * x[1],
           ρ * x[1] - x[1] * x[3] - x[2],
        x[1] * x[2] -    β * x[3]
    )
    return dx
end

F(x) = rk4_flow_map(f, x)

# Internally, GAIO calls this function on a set of test points
# within the domain. This means many function calls have to be made.
# If your function only uses "basic" instructions, then it is
# possible to simultaneously apply
# Single Instructions to Multiple Data (SIMD).
# This way less total function calls have to be made, saving time.

# To see which instructions are supported, refer to
# https://github.com/eschnett/SIMD.jl.git

# All we need to do is pass :cpu to the BoxMap command.
center, radius = (0,0,25), (30,30,30)
P = BoxPartition(Box(center, radius), (128,128,128))
G = BoxMap(F, P, :cpu)

x = (sqrt(β*(ρ-1)), sqrt(β*(ρ-1)), ρ-1)
W = unstable_set!(G, P[x])

plot(W)


text:
module Hubbard
include("common.jl")
# include("hilbert.jl")
using ..Hilbert
# include("green.jl")
using ..Green
using LinearAlgebra
# import QuantumStatistics.Basis:tau2matfreq, tau2dlr

function fermiHubbard(t, U, μ, sites, bonds)
    Nsite = length(sites)
    Fock = Hilbert.BinaryFock(Nsite)
    cu⁺ = [Hilbert.creation(Fock, s, UP) for s in sites]
    cd⁺ = [Hilbert.creation(Fock, s, DOWN) for s in sites]
    cu⁻ = [c' for c in cu⁺]
    cd⁻ = [c' for c in cd⁺]
    nu = [cu⁺[s] * cu⁻[s] for s in sites]
    nd = [cd⁺[s] * cd⁻[s] for s in sites]

    K = sum([cu⁺[i] * cu⁻[j] + cd⁺[i] * cd⁻[j] for (i, j) in bonds])
    V = sum(nu[s] * nd[s] for s in sites)
    C = sum(nu[s] + nd[s] for s in sites)

    # show(stdout, "text/plain", Matrix(K))
    # println()

    H = -t * K + U * V - μ * C

    # show(stdout, "text/plain", Matrix(H))
    # println()

    return H, cu⁺, cd⁺
end

function hubbardAtom(type, U, μ, β)
    E = [0.0, -μ, -μ, U - 2μ]
    H = zeros(Float, (4, 4))
    H[diagind(H)] = E

    # |0>=|00>=1, |↑>=|10>=2, |↓>=|01>=3, |↑↓>=|11>=4
    # the first is the forck state for ↑ spin, the second is forck state for ↓

    cpup = zeros(Float, (4, 4))
    cpdown = zeros(Float, (4, 4))

    cpup[2, 1], cpup[4, 3] = 1, 1
    cpdown[3, 1], cpdown[4, 2] = 1, 1
    cmup, cmdown = cpup', cpdown'

    @assert abs(tr(cpup * cmdown)) < 1e-16
    @assert abs(tr(cpdown * cmup)) < 1e-16

    m = Model(β, H, [cpup, cpdown])
    println("U=$U, μ=$μ, β=$β")
    println("Model Hilbert space: $(m.dim)")
    println("The lowest eigen energy per site:\n$(E)")
    return m
end

function hubbardAtom2(type, t, U, μ, β)
    Nsite = 2
    bonds = [(1, 2), (2, 1)]
    sites = [s for s in 1:Nsite]

    if type == :fermi
        H, cu⁺, cd⁺ = fermiHubbard(t, U, μ, sites, bonds)
    else
        @error("Not implemented!")
    end

    m = Model(β, H, hcat(cu⁺, cd⁺))

    println("t=$t, U=$U, μ=$μ, β=$β")
    println("Model Hilbert space: $(m.dim)")
    println("The lowest eigen energy per site:\n$(m.E[1:8] / Nsite)")
    # E, U = eigen(Matrix(H))
    # println(E)
    return m
end

function hubbardAtom4(type, t, U, μ, β)
    Nsite = 4
    bonds = [(1, 2), (2, 1), (2, 3), (3, 2), (3, 4), (4, 3), (1, 4), (4, 1)]
    sites = [s for s in 1:Nsite]

    if type == :fermi
        H, cu⁺, cd⁺ = fermiHubbard(t, U, μ, sites, bonds)
    else
        @error("Not implemented!")
    end

    m = Model(β, H, hcat(cu⁺, cd⁺))

    println("t=$t, U=$U, μ=$μ, β=$β")
    println("Model Hilbert space: $(m.dim)")
    println("The lowest eigen energy per site:\n$(m.E[1:16] / Nsite)")
    # E, U = eigen(Matrix(H))
    # println(m.E[1])
    return m
end

# const m = Model(10.0, 5.0, β=10.0)
# propagator(m)
end


text:
# This is part of Calendars.jl. See the copyright note there.
# ========================= Julian dates ====================

# Day number of the start of the Julian calendar.
# JL-0001-01-01 = CE-0000-12-30  ~ DN#-1 = EpochJulian.
# JL-0001-01-03 = CE-0001-01-01  ~ DN#1.
# JD-1111-11-11 = CE-1111-11-18  ~ DN#0405740
# JD-9999-12-31 = CE-10000-03-13 ~ DN#3652132
const EpochJulian = -1
const ValidYearsJulian = (1, 9999) # by convention

# True if year is a Julian leap year.
function isLeapYearJulian(year::DPart)
    rem(year, 4) == 0
end

# Return the last day of the month for the Julian calendar.
function LastDayOfMonthJulian(year::DPart, month::DPart)
    leap = isLeapYearJulian(year)
    (month == 2 && leap) && return 29
    return MonthLength[month]
end

# Is the date a valid Julian date?
function isValidDateJulian(cd::CDate, warn=true)
    cal, year, month, day = cd
    val = ( cal == JD
        && (ValidYearsJulian[1] <= year && year <= ValidYearsJulian[2])
        && (1 <= month && month <= 12)
        && (1 <= day && day <= LastDayOfMonthJulian(year, month)))
    ! val && warn && @warn(Warning(cd))
    return val
end

# Represents the day number a valid Julian date?
function isValidDateJulian(dn::DPart, warn=true)
    val = EpochJulian <= dn
    ! val && warn && @warn(Warning(dn))
    return val
end

# Return the days this year so far.
function DayOfYearJulian(year::DPart, month::DPart, day::DPart)
    for m in (month - 1):-1:1  # days in prior months this year
        day += LastDayOfMonthJulian(year, m)
    end
    return day - 1
end

# Returns the day number from a valid Julian date.
function DayNumberJulian(year::DPart, month::DPart, day::DPart)
    day = DayOfYearJulian(year, month, day)

    return (EpochJulian      # days elapsed before absolute date 1
        + 365 * (year - 1)   # days in previous years ignoring leap days
        + div(year - 1, 4)   # leap days before this year...
        + day                # days this year
    )
end
DayNumberJulian(d::CDate) = DayNumberJulian(d[2], d[3], d[4])

# Computes the Julian date from the day number.
function DateJulian(dn::DPart)
    ! isValidDateJulian(dn) && return InvalidDate

    # Search forward year by year from approximate year
    year = div(dn + EpochJulian, 366)

    while dn >= DayNumberJulian(year + 1, 1, 1)
        year += 1
    end

    # Search forward month by month from January
    month = 1
    while dn > DayNumberJulian(year, month,
                    LastDayOfMonthJulian(year, month))
        month += 1
    end

    day = dn - DayNumberJulian(year, month, 1) + 1
    return (JD, year, month, day)::CDate
end

# Return the day number of the first day in the given Julian year.
function YearStartJulian(year::DPart)
    return DayNumberJulian(year, 1, 1)
end

# Return the day number of the last day in the given Julian year.
function YearEndJulian(year::DPart)
    return DayNumberJulian(year + 1, 1, 1) - 1
end

# Return the number of months in the given Julian year.
function MonthsInYearJulian(year::DPart)
    return 12
end

# Return the number of days in the given Julian year.
function DaysInYearJulian(year::DPart)
    return YearStartJulian(year + 1) - YearStartJulian(year)
end


text:
# Example for basic 2D modeling:
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

using JUDI, SegyIO, LinearAlgebra, PyCall

# Set up model structure
n = (120, 100)   # (x,y,z) or (x,z)
d = (10., 10.)
o = (0., 0.)

# Velocity [km/s]
v = ones(Float32,n) .+ 0.5f0
v0 = ones(Float32,n) .+ 0.5f0
v[:,Int(round(end/2)):end] .= 3.5f0
rho = (v0 .+ .5f0) ./ 2

# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2
m0 = (1f0 ./ v0).^2
dm = vec(m - m0)

# Setup info and model structure
nsrc = 2	# number of sources
model = Model(n, d, o, m)
model0 = Model(n, d, o, m0)

# Set up receiver geometry
nxrec = 120
xrec = range(50f0, stop=1150f0, length=nxrec)
# yrec = range(0f0, stop=0f0, length=nxrec)
yrec = 0f0
zrec = range(0f0, stop=0f0, length=nxrec)

# receiver sampling and recording time
timeR = 1000f0   # receiver recording time [ms]
dtR = 2f0    # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

# Set up source geometry (cell array with source locations for each shot)
xsrc = convertToCell(range(400f0, stop=800f0, length=nsrc))
ysrc = convertToCell(range(0f0, stop=0f0, length=nsrc))
zsrc = convertToCell(range(200f0, stop=200f0, length=nsrc))

# source sampling and number of time steps
timeS = 1000f0  # ms
dtS = 2f0   # ms

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtS, t=timeS)

# setup wavelet
f0 = 0.01f0     # kHz
wavelet = ricker_wavelet(timeS, dtS, f0)
q = judiVector(srcGeometry, wavelet)

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry, recGeometry, model)
info = Info(prod(n), nsrc, ntComp)

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

h5geo = pyimport("h5geopy._h5geo")

geom_filename = "/home/kerim/Documents/Colada_prj/original JUDI examples/modeling_basic_2d/M_400.h5"
# geom_filename = "/home/kerim/Documents/Colada_prj/original JUDI examples/modeling_basic_2d/M_400_NOT_MAPPED.h5"
geom_name = "M_400"

geomCnt = h5geo.openSeisContainerByName(geom_filename)
if isnothing(geomCnt)
  @error "Unable to open geom container: $geom_filename"
  return
end

geom = geomCnt.openSeis(geom_name)
if isnothing(geom)
  @error "Unable to open geom: $geom_name"
  return
end

geom.removePKeySort("SRCX")
geom.addPKeySort("SRCX")

con = H5SeisCon(seis=geom, pkey="SRCX")
dobs = judiVector(con)

# Scan directory for segy files and create out-of-core data container
container = segy_scan("/home/kerim/Documents/Colada_prj/original JUDI examples/modeling_basic_2d/useful/", "M_",
                     ["GroupX", "GroupY", "RecGroupElevation", "SourceSurfaceElevation", "dt"])
dobs_segyio = judiVector(container)

# temporary set geom to avoid round off mismatch error
# recGeometry = dobs.geometry
dobs.geometry = recGeometry

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

# Write shots as segy files to disk
opt = Options(optimal_checkpointing=false, isic=false, subsampling_factor=2, dt_comp=1.0)

# Setup operators
Pr = judiProjection(info, recGeometry)
F = judiModeling(info, model; options=opt)
F0 = judiModeling(info, model0; options=opt)
Ps = judiProjection(info, srcGeometry)
J = judiJacobian(Pr*F0*adjoint(Ps), q)

# Nonlinear modeling
# dobs = Pr*F*adjoint(Ps)*q

# opt.file_path = "/home/kerim/Documents/Colada_prj/original JUDI examples/modeling_basic_2d/"
# opt.file_name = "MM"

# # Adjoint

pdeFull = Ps*adjoint(F)*adjoint(Pr)
qad = pdeFull*dobs

# Linearized modeling
dD = J*dm
# Adjoint jacobian
rtm = adjoint(J)*dD

# evaluate FWI objective function
f, g = fwi_objective(model0, q, dobs; options=opt)

# TWRI
f, gm, gy = twri_objective(model0, q, dobs, nothing; options=opt, optionswri=TWRIOptions(params=:all))
f, gm = twri_objective(model0, q, dobs, nothing; options=Options(frequencies=[[.009, .011], [.008, .012]]),
                                                 optionswri=TWRIOptions(params=:m))

# evaluate LSRTM objective function
fj, gj = lsrtm_objective(model0, q, dD, dm; options=opt)
fjn, gjn = lsrtm_objective(model0, q, dobs, dm; nlind=true, options=opt)

# By extension, lsrtm_objective is the same as fwi_objecive when `dm` is zero
# And with computing of the residual. Small noise can be seen in the difference
# due to floating point roundoff errors with openMP, but running with
# OMP_NUM_THREAS=1 (no parllelism) produces the exact (difference == 0) same result
# gjn2 == g
fjn2, gjn2 = lsrtm_objective(model0, q, dobs, 0f0.*dm; nlind=true, options=opt)


text:
using Test
using LinearAlgebra
using QuantumPropagators
using QuantumControlBase
using QuantumControlBase.TestUtils
using Zygote

@testset "GradGenerator" begin

    N = 10  # size of Hilbert space
    ρ = 1.0  # spectral radius
    Ĥ₀ = random_hermitian_matrix(N, ρ)
    Ĥ₁ = random_hermitian_matrix(N, ρ)
    Ĥ₂ = random_hermitian_matrix(N, ρ)
    Zero = zeros(ComplexF64, N, N)
    ϵ₁ = t->1.0
    ϵ₂ = t->1.0
    Ĥ_of_t = (Ĥ₀, (Ĥ₁, ϵ₁), (Ĥ₂, ϵ₂))
    Ψ = random_state_vector(N)
    ψ_max = maximum(abs.(Ψ))
    Ψtgt = random_state_vector(N)
    𝕚 = 1im
    dt  = 1.25

    vals_dict = IdDict(ϵ₁ => 1.0, ϵ₂ => 1.0)

    G̃_of_t = TimeDependentGradGenerator(Ĥ_of_t)
    G̃ = evalcontrols(G̃_of_t, vals_dict)
    Ĥ = evalcontrols(Ĥ_of_t, vals_dict)

    Û_Ψ = exp(-𝕚*Ĥ*dt) * Ψ

    Ψ̃ = GradVector(Ψ, 2)
    # did the initialization work?
    @test norm(Ψ̃.state - Ψ) < 1e-14
    @test norm(Ψ̃.grad_states[1]) == 0.0
    @test norm(Ψ̃.grad_states[2]) == 0.0
    @test length(Ψ̃.grad_states) == 2

    @test maximum(abs.(Ψ)) == ψ_max  # is Ψ still exactly the same state?
    wrk = NewtonWrk(Ψ̃)
    newton!(Ψ̃, G̃, dt, wrk)
    Ψ̃_out = copy(Ψ̃)
    @test maximum(abs.(Ψ)) == ψ_max  # is Ψ still exactly the same state?
    @test norm(Ψ̃_out.state - Û_Ψ) < 1e-12  # Ψ̃_out.state correct?
    @test norm(Ψ̃_out.grad_states[1]) > 0
    @test norm(Ψ̃_out.grad_states[2]) > 0
    # did the copy work?
    @test norm(Ψ̃_out.state - Ψ̃.state) == 0.0
    @test norm(Ψ̃_out.grad_states[1] - Ψ̃.grad_states[1]) == 0.0
    @test norm(Ψ̃_out.grad_states[2] - Ψ̃.grad_states[2]) == 0.0

    # reset Ψ̃ to the original input (also to test the resetgradvec! method)
    resetgradvec!(Ψ̃, Ψ)
    # state should be the same as after previous initialization
    @test norm(Ψ̃.state - Ψ) < 1e-14
    @test norm(Ψ̃.grad_states[1]) == 0.0
    @test norm(Ψ̃.grad_states[2]) == 0.0

    @test norm(Ψ̃_out.state - Û_Ψ) < 1e-12  # Ψ̃_out.state still correct?

    ###########################################################################
    # Compare against explicit Grad-Gen
    # This checks whether the application of a GradGenerator to a GradVector an
    # all the linear-algebra methods are implemented correctly

    G̃_full = vcat(hcat(Ĥ,    Zero, Ĥ₁),
                  hcat(Zero, Ĥ,    Ĥ₂),
                  hcat(Zero, Zero, Ĥ))

    Ψ̃_full = vcat(Ψ̃.grad_states[1],
                  Ψ̃.grad_states[2],
                  Ψ̃.state)
    # proper initialization? grad_states should be zero
    @test norm(Ψ̃_full) == norm(Ψ̃.state) == norm(Ψ)

    Ψ̃_out_full = exp(-𝕚 * G̃_full * dt) * Ψ̃_full
    # propagation correct?
    @test norm(Ψ̃_out_full[2N+1:3N] - Û_Ψ) < 1e-12

    # do we get the same results as from newton?
    @test norm(Ψ̃_out_full[2N+1:3N] - Ψ̃_out.state) < 1e-12
    @test norm(Ψ̃_out_full[1:N] - Ψ̃_out.grad_states[1]) < 1e-12
    @test norm(Ψ̃_out_full[N+1:2N] - Ψ̃_out.grad_states[2]) < 1e-12

    ###########################################################################
    # Test custom expprop

    wrk_exp = initpropwrk(Ψ̃, [0, dt], Val(:expprop), G̃)
    propstep!(Ψ̃, G̃, dt, wrk_exp)
    Ψ̃_out_exp = copy(Ψ̃)
    @test norm(wrk_exp.Ψ_full - Ψ̃_full) ≈ 0.0
    @test norm(wrk_exp.G_full - G̃_full) ≈ 0.0
    @test norm(Ψ̃_out_exp - Ψ̃_out) < 1e-12
    resetgradvec!(Ψ̃, Ψ)

    ###########################################################################
    # Compare against explicit split-up Grad-Gens
    # This checks whether the extension of the Gradient-Vector to multiple
    # controls is correct (the literature generally only gives the expression
    # for a single control)

    G̃_full1 = vcat(hcat(Ĥ,    Ĥ₁),
                   hcat(Zero, Ĥ))
    G̃_full2 = vcat(hcat(Ĥ,    Ĥ₂),
                   hcat(Zero, Ĥ))
    Ψ̃_full1 = vcat(Ψ̃.grad_states[1],
                   Ψ̃.state)
    @test maximum(abs.(Ψ)) == ψ_max  # is Ψ still exactly the same state?
    @test norm(Ψ̃_full1) == norm(Ψ)  # initialization correct?
    Ψ̃_full2 = vcat(Ψ̃.grad_states[2],
                   Ψ̃.state)
    @test norm(Ψ̃_full2) == norm(Ψ)  # initialization correct?
    Ψ̃_out_full1 = exp(-𝕚 * G̃_full1 * dt) * Ψ̃_full1
    Ψ̃_out_full2 = exp(-𝕚 * G̃_full2 * dt) * Ψ̃_full2

    # propagation correct?
    @test norm(Ψ̃_out_full1[N+1:2N] - Û_Ψ) < 1e-12
    @test norm(Ψ̃_out_full2[N+1:2N] - Û_Ψ) < 1e-12

    # do we get the same results as with the combined grad-gen?
    @test norm(Ψ̃_out_full1[1:N] - Ψ̃_out_full[1:N]) < 1e-12
    @test norm(Ψ̃_out_full2[1:N] - Ψ̃_out_full[N+1:2N]) < 1e-12
    @test norm(Ψ̃_out_full1[N+1:2N] - Ψ̃_out_full[2N+1:3N]) < 1e-12
    @test norm(Ψ̃_out_full2[N+1:2N] - Ψ̃_out_full[2N+1:3N]) < 1e-12

    ###########################################################################
    # Compare against Zygote
    # This checks whether the gradients are correct, albeit a bit indirectly
    # (since Zygote can only calculate the gradient for a scalar function, in
    # this case the square-modulus of the overlap with a target)

    F_sm(ϵ₁, ϵ₂) = abs(dot(Ψtgt, exp(-1im * (Ĥ₀ + ϵ₁*Ĥ₁ + ϵ₂*Ĥ₂) * dt) * Ψ))^2
    grad_zygote = collect(gradient(F_sm, 1.0, 1.0))

    @test norm(Ψ̃_out.state - Û_Ψ) < 1e-12  # still correct?
    τ = dot(Ψtgt, Û_Ψ)
    # `grad` is gradient of F_sm based on Newton-prop of GradGenerator
    # For ∂F/∂τ see Eq. (3.47) of Phd Thesis of Michael Goerz
    grad = [2*real(conj(τ) * dot(Ψtgt, Ψ̃_out.grad_states[1])),
            2*real(conj(τ) * dot(Ψtgt, Ψ̃_out.grad_states[2]))]
    @test abs(grad_zygote[1] - grad[1]) < 1e-10
    @test abs(grad_zygote[2] - grad[2]) < 1e-10

    ###########################################################################
    # Compare against Taylor series
    # This checks the gradients more directly against an alternative method

    commutator(A, B) = A * B - B * A

    """Evaluate ∂/∂ϵ exp(-𝕚 Ĥ dt) via a Taylor expansion."""
    function U_grad(Ĥ, μ̂, dt)
        # See Eq. (14) in de Fouquieres et. al, JMR 212, 412 (2011)
        Û = exp(-𝕚 * Ĥ * dt)
        converged = false
        Ĉ = μ̂
        terms = [(-𝕚 * dt) * Ĉ]
        n = 2
        while !converged
            Ĉ = commutator(Ĥ, Ĉ)
            term = -((𝕚 * dt)^n / factorial(big(n))) * Ĉ
            push!(terms, term)
            converged = (norm(term) < 1e-12)
            n += 1
        end
        return Û * sum(terms)
    end

    grad_taylor = [U_grad(Ĥ, Ĥ₁, dt) * Ψ,
                   U_grad(Ĥ, Ĥ₂, dt) * Ψ]

    @test norm(Ψ̃_out.grad_states[1] - grad_taylor[1]) < 1e-10
    @test norm(Ψ̃_out.grad_states[2] - grad_taylor[2]) < 1e-10

end


text:
using DelimitedFiles
segdata = readdlm("segments.txt",' ')
#@show segdata

function initdisplay()

    display = Vector(undef,10)

    display[1] =  "abcefg"
    display[2] =  "cf"
    display[3] =  "acdeg"
    display[4] =  "acdfg"
    display[5] =  "bcdf"
    display[6] =  "abdfg"
    display[7] =  "abdefg"
    display[8] =  "acf"
    display[9] =  "abcdefg"
    display[10] = "abcdfg"
    #@show display

    displaylen = zeros(Int,1,10)
    for d in eachindex(display)
        displaylen[d] = length(display[d])
    end

    display,displaylen
end

display,displaylen = initdisplay()

@show displaylen

function findeasydigits(segdata)
    uniqseg = 0
    for segment in 1:size(segdata,1)
        #@show segment,segdata[segment,:]
        foundoutput = false
        for sigdata in segdata[segment,:]
            if foundoutput
                #okay, see if it is a unique display
                siglen = length(sigdata)
                #@show sigdata,siglen
                segmatch = count(==(siglen),displaylen)
                if segmatch == 1
                    uniqseg += 1
                end
            elseif sigdata == "|"
                foundoutput = true
            end
        end
    end
    uniqseg
end

@show findeasydigits(segdata)


text:
export MenuItem, buildmenu

mutable struct MenuItem
    txt::AbstractString
    cb::Function

    function MenuItem(txt,cb)
        new(txt,cb)
    end
end

"""
Build a GtkMenu from an Array of MenuItems,
a parent Menu and user_data that will be passed
to callbacks.

Example:

    buildmenu([
        MenuItem("New File",newMenuItem_activate_cb),
        GtkSeparatorMenuItem,
        MenuItem("Quit",quitMenuItem_activate_cb),
        ],
        fileMenu,
        (main_window)
    )
"""
function buildmenu(items::Array,menu::GtkMenu,user_data)
    for i in items
        if typeof(i) == MenuItem
            mi = GtkMenuItem(i.txt)
            push!(menu,mi)
            signal_connect(i.cb, mi, "activate", Nothing,(),false,user_data)
        elseif i == GtkSeparatorMenuItem
            push!(menu,GtkSeparatorMenuItem())
        end
    end
    showall(menu)
end
buildmenu(items::Array,menuItem::GtkMenuItem,user_data) = buildmenu(items,GtkMenu(menuItem),user_data)
buildmenu(items::Array,user_data) = buildmenu(items,GtkMenu(),user_data)
buildmenu(items::MenuItem,user_data) = buildmenu([items],user_data)

max_stars_repo_name: jonathanBieler/GtkExtensions

max_issues_repo_name: jonathanBieler/GtkExtensions

max_forks_repo_name: jonathanBieler/GtkExtensions


text:
export Camera, PerspectiveCamera, OrthoCamera
export xy_to_uv, uv_to_xy, cast_ray, project_point

abstract type Camera end

mutable struct PerspectiveCamera <: Camera
    position::Vec3d
    up::Vec3d
    right::Vec3d
    dir::Vec3d

    PerspectiveCamera(pos, look_at, up, right) = new(pos, up, right,
                                                     normalize(look_at - pos))
end

mutable struct OrthoCamera <: Camera
    position::Vec3d
    up::Vec3d
    right::Vec3d
    dir::Vec3d

    OrthoCamera(pos, up, right) = new(pos, up, right, right × up)
end

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

"""
    normalize_coord(x, span)

Assuming that `x` runs in the range [1, `span`], apply a linear
transformation to `x` so that the result runs from -1 to +1.
"""
function normalize_coord(x, span)
    span == 1 && return 0.0

    2(x - 1 + 0.5) / (span - 1) - 1
end

xy_to_uv(x, y, width, height) = (width / height * normalize_coord(x, width),
                                 normalize_coord(y, height))

uv_to_xy(u, v, width, height) = (((height / width * u + 1) * (width - 1) + 1) / 2,
                                 ((v + 1) * (height - 1) + 1) / 2)

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

"""
    cast_ray(camera::PerspectiveCamera, x, y, width, height)

Cast a ray starting from the camera along the direction `(x, y)`,
where both `x` and `y` are pure numbers in the range [1, width] × [1,
height]. Return a ray (a `Vec3d` object) of length one.
"""
function cast_ray(camera::PerspectiveCamera, x, y, width, height)
    u, v = xy_to_uv(x, y, width, height)
    dir = camera.dir + u * camera.right + v * camera.up

    Ray(camera.position, dir)
end

"""
    cast_ray(camera::OrthoCamera, x, y, width, height)

Cast a ray starting from the camera along the direction `(x, y)`,
where both `x` and `y` are pure numbers in the range [1, width] × [1,
height]. Return a ray (a `Vec3d` object) of length one.
"""
function cast_ray(camera::OrthoCamera, x, y, width, height)
    pos = (camera.position +
           width / height * normalize_coord(x, width) * camera.right +
           normalize_coord(y, height) * camera.up)

    Ray(pos, camera.dir)
end

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

function project_point(camera::OrthoCamera, pt, width, height)
    # For orthogonal cameras, we trace the ray in the inverse
    # direction: from the point `pt` backwards to the
    # plane. Therefore, the ray would be `Ray(pt, -camera.dir)`.

    t = (pt - camera.position) ⋅ camera.dir / norm2(camera.dir)
    pt_on_screen = pt - t * camera.dir

    # Project the point on a (u, v) coordinate system, where both `u`
    # and `v` are in the range [-1, 1]
    (u, v) = (pt_on_screen ⋅ camera.right, pt_on_screen ⋅ camera.up)
    (x, y) = uv_to_xy(u, v, width, height)

    return (x, y, 1 ≤ x ≤ width && 1 ≤ y ≤ height && t > 0)
end

"""
    project_point(camera::PerspectiveCamera, pt, width, height)

Project a 3D point into the camera frame. Return a 3-element tuple
representing the coordinates ``(x, y)`` of a point ``P`` on the screen
and a Boolean value indicating whether the point ``P`` is in the range
`[1, width] × [1, height]` or not (visibility).

Both ``x`` and ``y`` are floating-point numbers.
"""
function project_point(camera::PerspectiveCamera, pt, width, height)
    # Direction from the camera to the point to project on the screen
    dir = pt - camera.position
    plane_center = camera.position + camera.dir

    # Check that the direction is not perpendicular to the viewing direction
    dotprod = camera.dir ⋅ dir
    abs(dotprod) < 1e-7 && return (0.0, 0.0, false)

    # Distance between P (on the screen plane) and the observer, along `dir`
    t = norm2(camera.dir) / dotprod

    # Point on the screen (3D coordinates)
    pt_on_screen = t * dir - camera.dir

    # Project the point on a (u, v) coordinate system, where both `u`
    # and `v` are in the range [-1, 1]
    (u, v) = (pt_on_screen ⋅ camera.right, pt_on_screen ⋅ camera.up)
    (x, y) = uv_to_xy(u, v, width, height)

    return (x, y, 1 ≤ x ≤ width && 1 ≤ y ≤ height && t > 0)
end


text:
# ---
# title: 196. Delete Duplicate Emails
# id: problem196
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories:
# link: <https://leetcode.com/problems/delete-duplicate-emails/description/>
# hidden: true
# ---
#
# Write a SQL query to **delete** all duplicate email entries in a table named
# `Person`, keeping only unique emails based on its _smallest_ **Id**.
#
#
#
#     +----+------------------+
#     | Id | Email            |
#     +----+------------------+
#     | 1  | john@example.com |
#     | 2  | bob@example.com  |
#     | 3  | john@example.com |
#     +----+------------------+
#     Id is the primary key column for this table.
#
#
# For example, after running your query, the above `Person` table should have
# the following rows:
#
#
#
#     +----+------------------+
#     | Id | Email            |
#     +----+------------------+
#     | 1  | john@example.com |
#     | 2  | bob@example.com  |
#     +----+------------------+
#
#
# **Note:**
#
# Your output is the whole `Person` table after executing your sql. Use `delete`
# statement.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#=
Simple unit Tests
to test formulas
=#

using Test
using Dates

include("financeHelper.jl");

println("#### Test Starting")

initialCapital = 100.00                                                         #The initial capital as float
profit_loss_series_one = [10.0, 10.0, 10.0, 10.0, 10.0];                        #Scenario "No Drawdown"
profit_loss_series_two = [-5.0, -5, 0, 10.0, 10, 10, 10, -5, -5, -5, 10.0];     #Scenario "Drawdown somewhere"
profit_loss_series_three = [-5,-5.0];                                           #Scenario "only negativ trades"
profit_loss_series_four =
    [5.5, -5.5, -5, 15, -5, -5, -5, +15, +5, +5, +5, +5, +5, -10, 5, 5, 5,-5];           #Scenario "two times same DD"

@testset "getNominalReturn" begin
    @test 50 == getNominalReturn(profit_loss_series_one)
    @test 25 == getNominalReturn(profit_loss_series_two)
    @test -10 == getNominalReturn(profit_loss_series_three)
    @test 35 == getNominalReturn(profit_loss_series_four)
end

@testset "getAccumulatedProfitLossSeries" begin
    @test [10.0, 20.0, 30.0, 40.0, 50.0] ==
          getAccumulatedProfitLossSeries(profit_loss_series_one)
    @test [-5, -10, -10, 0, 10, 20, 30, 25, 20, 15, 25] ==
          getAccumulatedProfitLossSeries(profit_loss_series_two)
    @test [-5,-10.0] ==
          getAccumulatedProfitLossSeries(profit_loss_series_three)
    @test [5.5,0,-5,10,5,0,-5,10,15,20,25,30,35,25,30,35,40,35] ==
          getAccumulatedProfitLossSeries(profit_loss_series_four)
end

@testset "getMaxDD" begin
    @test 0 == getMaxDD(profit_loss_series_one)
    @test 15 == getMaxDD(profit_loss_series_two)
    @test 10 == getMaxDD(profit_loss_series_three)
    @test 15 == getMaxDD(profit_loss_series_four)
end

@testset "getMaxDDRate_relativeToLastPeak" begin
    @test Inf == getMaxDDRate_relativeToLastPeak(profit_loss_series_one)
    @test 0.5 == getMaxDDRate_relativeToLastPeak(profit_loss_series_two)
    @test -Inf == getMaxDDRate_relativeToLastPeak(profit_loss_series_three)
    @test 1.5 == getMaxDDRate_relativeToLastPeak(profit_loss_series_four)
end

@testset "getReturnRate" begin
    @test 0.5 == getReturnRate(profit_loss_series_one, initialCapital)
    @test 0.25 == getReturnRate(profit_loss_series_two, initialCapital)
end

@testset "getReturnRatePercentage" begin
    @test 50 == getReturnRatePercentage(profit_loss_series_one, initialCapital)
    @test 25 == getReturnRatePercentage(profit_loss_series_two, initialCapital)
end

@testset "getAnnulizedReturnRateBetweenDatesYearly" begin
    #For full year
    @test 0.5004166315911229 == getAnnulizedReturnRateBetweenDatesYearly(
        profit_loss_series_one,
        initialCapital,
        Date(2018, 1, 1),
        Date(2018, 12, 31),
    )
     #For 2 years
    @test 0.224914948717307 == getAnnulizedReturnRateBetweenDatesYearly(
        profit_loss_series_one,
        initialCapital,
        Date(2018, 1, 1),
        Date(2019, 12, 31),
    )
    #For 10 days
    @test 2.7022969825829207e6 == getAnnulizedReturnRateBetweenDatesYearly(
        profit_loss_series_one,
        initialCapital,
        Date(2018, 1, 1),
        Date(2018, 1, 10),
    )
end

@testset "getDaysBetweenDates_includeLast" begin
    @test 2 == getDaysBetweenDates_includeLast(
        Dates.Date(2018, 01, 01),
        Dates.Date(2018, 01, 02),
    )
    @test 365 == getDaysBetweenDates_includeLast(
        Dates.Date(2018, 01, 01),
        Dates.Date(2018, 12, 31),
    )
    @test 365 == getDaysBetweenDates_includeLast(
        Dates.Date(2019, 01, 01),
        Dates.Date(2019, 12, 31),
    )
    @test 366 == getDaysBetweenDates_includeLast(
        Dates.Date(2020, 01, 01),
        Dates.Date(2020, 12, 31),
    )
    @test 1096 == getDaysBetweenDates_includeLast(
        Dates.Date(2018, 01, 01),
        Dates.Date(2020, 12, 31),
    )
end

println("#### Test Ended")

max_stars_repo_name: RainerNeu/tradeStatistics

max_issues_repo_name: RainerNeu/tradeStatistics

max_forks_repo_name: RainerNeu/tradeStatistics


text:
"""
This struct mean all atomic information other than coordinates.
This contains 5 fields like below.

- resname  : Union{String,  Nothing}
- resid    : Union{Int64,   Nothing}
- atomname : Union{String,  Nothing}
- atomid   : Union{Int64,   Nothing}
- mass     : Union{Float32, Nothing}

The constructor is below.

    function Attribute(;resname = nothing, resid = nothing,
                       atomname = nothing, atomid = nothing, mass = nothing)
"""
mutable struct Attribute

    resname::Union{String, Nothing}
    resid::Union{Int64, Nothing}
    atomname::Union{String, Nothing}
    atomid::Union{Int64, Nothing}
    mass::Union{Float32, Nothing}

    function Attribute(;resname = nothing, resid = nothing,
        atomname = nothing, atomid = nothing, mass = nothing)

        new(resname, resid, atomname, atomid, mass)
    end
end

"""
This struct have all information of corresponding atom.
This contains 2 fields like below.

- coordinate : Reference to [`Coordinate`](@ref) object.
- attribute  : Reference to [`Attribute`](@ref) object.
"""
mutable struct Atom
    coordinate::Coordinate{Float32}
    attribute::Attribute
end

function Atom(coordinate)
    Atom(coordinate, Attribute())
end

"""
This struct corresponding to MD trajectory.
This contains 4 fields like below.

- coordinates : AbstractArray{[`Coordinate`](@ref), 2}
- attributes  : Vector{[`Attribute`](@ref)}
- natom       : Number of atoms.
- nframe      : Number of frames in the trajectory.

The column of coordinates matrix means one snapshot.
The row of coordinates matrix means time series of one atom.

```julia
[ a d g j
  b e h k
  c f i l ]
```

In this case, one snapshot correspond to `[a b c]`.

The constructor is below.

    function Trajectory(coordinates::Array{Coordinate, 2},
                        attributes::Vector{Attribute} = [Attribute() for i=1:length(coordinates)])
"""
mutable struct Trajectory{RealT <: Real}
    coordinates::Array{Coordinate{RealT}}
    attributes::Vector{Attribute}
    natom::Int64
    nframe::Int64

    function Trajectory(coordinates::Array{<:Coordinate{RealT}, 2},
        attributes::Array{Attribute, 1} = [Attribute() for i=1:size(coordinates, 1)]
        ) where RealT <: Real

        new{RealT}(coordinates, attributes, size(coordinates, 1), size(coordinates, 2))
    end
end

"""
This struct correspond to one snapshot of trajectory.
This contains 3 fields like below.

- coordinates : Vector{[`Coordinate`](@ref)}
- attribute   : Vector{[`Attribute`](@ref)}
- natom       : Number of atoms.
"""
mutable struct Frame{RealT <: Real}
    coordinates::Vector{Coordinate{RealT}}
    attributes::Vector{Attribute}
    natom::Int64

    function Frame(coordinates::Vector{Coordinate{RealT}},
        attributes::Array{Attribute, 1} = [Attribute() for i=1:length(coordinates)]
        ) where RealT <: Real

        new{RealT}(coordinates, attributes, length(coordinates))
    end
end

mutable struct DCDMetaInfo
    file_size::Int64
    first_step::Int32
    nstep_save::Int32
    total_step::Int32
    total_unit::Int32
    time_step::Float32
    unitcell_flag::Bool
    version::Int32
    number_of_atom::Int32
    header_size::Int32
end


text:
using Base.Filesystem
using Base.Iterators
using JSON
using CSV: CSV
using Random

blacklist = ["lasagna/3787908",
                   "steak/1340977",
                   "bread_pudding/1375816"]

function make_test_pairs(datasetdir)
    p = joinpath(datasetdir,"meta","classes.txt")
    classes = CSV.read(p,header=false)
    classes = Dict([k_name for k_name in enumerate(convert(Array{String,1},classes[:Column1]))])
    json_file = "test.json"
    klass2path = JSON.parsefile(joinpath(datasetdir, "meta", json_file))
    pairs = []
    for (label, klass) in classes
        paths = klass2path[klass]
        for i in 1:length(paths)
            p = paths[i]
            if p in blacklist
                continue
            end
            push!(pairs, (joinpath(datasetdir, "images", p * ".jpg"),label))
        end
    end
    return pairs
end

function make_trainval_pairs(datasetdir;divide_ratio=0.95)
    p = joinpath(datasetdir,"meta","classes.txt")
    classes = CSV.read(p, header=false)
    classes = Dict([k_name for k_name in enumerate(convert(Array{String,1},classes[:Column1]))])
    json_file = "train.json"
    klass2path = JSON.parsefile(joinpath(datasetdir,"meta",json_file))
    train_pairs = []
    val_pairs = []
    for (label, klass) in classes
        paths = klass2path[klass]
        shuffle!(paths)
        for i in shuffle(1:length(paths))
            p = paths[i]
            if p in blacklist
                continue
            end
            if i<= length(paths) * divide_ratio
                push!(train_pairs, (joinpath(datasetdir,"images",p * ".jpg"),label))
            else
                push!(val_pairs, (joinpath(datasetdir,"images",p * ".jpg"),label))
            end
        end
    end
    return train_pairs, val_pairs
end


text:
struct Bleep
  bar
  baz::Int
  quuz::Float64
end

foo = Bleep("Hello, world", 23, 1.5)

println(typeof(foo))
println(foo)
println(foo.bar)
println(foo.baz)
println(foo.quuz)

function baz()
  x = 5
  x = "foo"
  return x
end

print(baz())

max_stars_repo_name: catseye/Dipple

max_issues_repo_name: catseye/Dipple

max_forks_repo_name: catseye/Dipple


text:
using BinDeps

BinDeps.@setup

path = dirname(@__FILE__)

libbouknap = library_dependency("libbouknap")
libminknap = library_dependency("libminknap")
libmulknap = library_dependency("libmulknap")

pisinger_webpage_uri = "http://hjemmesider.diku.dk/~pisinger"

pisinger_files_root = joinpath(path, "..", "libpisinger")
lib_pisinger_root = joinpath(path, "libpisinger")
lib_pisinger_minknap = joinpath(lib_pisinger_root, "minknap")
lib_minknap_build = joinpath(lib_pisinger_minknap, "build")
lib_pisinger_bouknap = joinpath(lib_pisinger_root, "bouknap")
lib_bouknap_build = joinpath(lib_pisinger_bouknap, "build")
lib_pisinger_mulknap = joinpath(lib_pisinger_root, "mulknap")
lib_mulknap_build = joinpath(lib_pisinger_mulknap, "build")

if Sys.iswindows()
    println("Build for windows.")
    using WinRPM
    println("Installing gcc.")
    WinRPM.install("gcc")

    println("Downloading cmake.")
    makedeplnk = "https://sourceforge.net/projects/gnuwin32/files/make/3.81/make-3.81-dep.zip/download"
    makebinlnk = "https://sourceforge.net/projects/gnuwin32/files/make/3.81/make-3.81-bin.zip/download"

    makedepzip = joinpath(pisinger_files_root, "make-3.81-dep.zip")
    makebinzip = joinpath(pisinger_files_root, "make-3.81-bin.zip")
    makebuilddir = joinpath(pisinger_files_root, "make")

    make = joinpath(makebuilddir, "bin", "make.exe")
    icnv = joinpath(makebuilddir, "bin", "libiconv2.dll")
    intl = joinpath(makebuilddir, "bin", "libintl3.dll")

    provides(
        SimpleBuild,
        (@build_steps begin
            FileDownloader(makedeplnk, makedepzip)
            FileDownloader(makebinlnk, makebinzip)
            FileUnpacker(makedepzip, makebuilddir, joinpath("bin", "libiconv2.dll"))
            FileUnpacker(makebinzip, makebuilddir, joinpath("bin", "make.exe"))
            FileDownloader(
                "$pisinger_webpage_uri/bouknap.c",
                joinpath(lib_pisinger_bouknap, "bouknap.c"),
            )
            CreateDirectory(lib_bouknap_build)
            @build_steps begin
                ChangeDirectory(lib_bouknap_build)
                `cmake ..`
                `$make`
            end
        end),
        libbouknap,
        os = :Windows,
    )

    provides(
        SimpleBuild,
        (@build_steps begin
            FileDownloader(
                "$pisinger_webpage_uri/minknap.c",
                joinpath(lib_pisinger_minknap, "minknap.c"),
            )
            CreateDirectory(lib_minknap_build)
            @build_steps begin
                ChangeDirectory(lib_minknap_build)
                `cmake ..`
                `$make`
                MakeTargets()
            end
        end),
        libminknap,
        os = :Windows,
    )

    provides(
        SimpleBuild,
        (@build_steps begin
            FileDownloader(
                "$pisinger_webpage_uri/mulknap.c",
                joinpath(lib_pisinger_mulknap, "mulknap.c"),
            )
            CreateDirectory(lib_mulknap_build)
            @build_steps begin
                ChangeDirectory(lib_mulknap_build)
                `cmake ..`
                `$make`
                MakeTargets()
            end
        end),
        libmulknap,
        os = :Windows,
    )
else
    println("Build not for windows.")
    provides(
        SimpleBuild,
        (@build_steps begin
            FileDownloader(
                "$pisinger_webpage_uri/bouknap.c",
                joinpath(lib_pisinger_bouknap, "bouknap.c"),
            )
            if Sys.isapple() # Removing deprecated headers for osx
                @build_steps begin
                    `sed -i '' '/values\.h/d' $lib_pisinger_bouknap/bouknap.c`
                    `sed -i '' '/malloc\.h/d' $lib_pisinger_bouknap/bouknap.c`
                end
            end
            CreateDirectory(lib_bouknap_build)
            @build_steps begin
                ChangeDirectory(lib_bouknap_build)
                `cmake ..`
                MakeTargets()
            end
        end),
        libbouknap,
        installed_libpath = lib_bouknap_build,
    )

    provides(
        SimpleBuild,
        (@build_steps begin
            FileDownloader(
                "$pisinger_webpage_uri/minknap.c",
                joinpath(lib_pisinger_minknap, "minknap.c"),
            )
            if Sys.isapple() # Removing deprecated headers for osx
                @build_steps begin
                    `sed -i '' '/values\.h/d' $lib_pisinger_minknap/minknap.c`
                    `sed -i '' '/malloc\.h/d' $lib_pisinger_minknap/minknap.c`
                end
            end
            CreateDirectory(lib_minknap_build)
            @build_steps begin
                ChangeDirectory(lib_minknap_build)
                `cmake ..`
                MakeTargets()
            end
        end),
        libminknap,
        installed_libpath = lib_minknap_build,
    )

    provides(
        SimpleBuild,
        (@build_steps begin
            FileDownloader(
                "$pisinger_webpage_uri/mulknap.c",
                joinpath(lib_pisinger_mulknap, "mulknap.c"),
            )
            if Sys.isapple() # Removing deprecated headers for osx
                @build_steps begin
                    `sed -i '' '/values\.h/d' $lib_pisinger_mulknap/mulknap.c`
                    `sed -i '' '/malloc\.h/d' $lib_pisinger_mulknap/mulknap.c`
                end
            end
            CreateDirectory(lib_mulknap_build)
            @build_steps begin
                ChangeDirectory(lib_mulknap_build)
                `cmake ..`
                MakeTargets()
            end
        end),
        libmulknap,
        installed_libpath = lib_mulknap_build,
    )
end

Sys.iswindows() && push!(BinDeps.defaults, BuildProcess)

BinDeps.@install Dict(
    :libbouknap => :_jl_libbouknap,
    :libminknap => :_jl_libminknap,
    :libmulknap => :_jl_libmulknap,
)

Sys.iswindows() && pop!(BinDeps.defaults)