text:
"""

    vdWenergy(myARGS)

Extract command-line arguments for computing the vdW interaction energy among a collection
of molecules, then compute the vdW energy integrand in imaginary frequency appropriately.

!!! note

    Any frequencies will be converted to their absolute values times `1im`. Frequencies are
    not allowed to be exactly zero.

!!! note

    This computes the integrand ``\\ln(\\det(\\mathbb{T}_{\\infty} \\mathbb{T}^{-1}))`` at
    specified imaginary frequencies, without doing any integration or Matsubara summation.
    This is typically a negative dimensionless quantity, and no factors of ``2\\pi`` or
    ``\\hbar`` are included.

!!! tip

    This function can be embarrassingly parallelized over frequency and wavevector. However,
    the integrand at each frequency and wavevector is calculated for every transformation,
    as this maximizes the efficiency of the program by virtue of reusing quantities
    independent of geometric transformations.

!!! tip

    The recommended way to do frequency integration or Matsubara summation is to supply a
    list of frequencies (typically not too many), and then interpolate the integrand to do
    integration or Matsubara summation after the fact.

"""
function vdWenergy(myARGS)

    mollistfilename = readRequiredFilenameARG(myARGS, "mollistfilename", "AARMBEM_vdW.jl");
    outfilename = readRequiredFilenameARG(myARGS, "outfilename", "AARMBEM_vdW.jl");

    translistfilename = readOptionalFilenameARG(myARGS, "translistfilename");
    periodicfilename = readOptionalFilenameARG(myARGS, "periodicfilename");

    Genvstr = readGenvstr(myARGS);

    nophonons, outfilename = readOptionalBoolcondARG(myARGS, "nophonons", outfilename);
    nonretarded, outfilename = readOptionalBoolcondARG(myARGS, "nonretarded", outfilename);

    FT = readOptionalFloatType(myARGS);

    myAllMolData = readMolSystem(mollistfilename, nophonons, FT);
    myTransData = readTransData(myAllMolData, translistfilename);
    myPeriodicData = readPeriodicData(FT, periodicfilename);

    freqlist = Array{FT, 1}(undef, 0);
    readfreqlist!(freqlist, myARGS, "AARMBEM_vdW.jl");

    klist = readklist(myARGS, "AARMBEM_vdW.jl", myPeriodicData.numdims, FT);

    freqklist = hcat(vec(repeat(transpose(freqlist), size(klist, 1), 1)),
                     repeat(klist, length(freqlist), 1));

    vdWenergy_general(outfilename, myAllMolData, myTransData, myPeriodicData, Genvstr,
                      nonretarded, freqklist)

end

"""

    vdWenergy_general(outfilename, myAllMolData, myTransData, myPeriodicData, Genvstr,
                      nonretarded, freqklist)

Compute the actual vdW energy integrand and print it to `outfilename`.

"""
function vdWenergy_general(outfilename::AbstractString,
                           myAllMolData::MolSystem{<:OneMol{FT}},
                           myTransData::TransData{FT},
                           myPeriodicData::PeriodicData{FT},
                           Genvstr::AbstractString, nonretarded::Bool,
                           freqklist::Array{FT, 2}) where FT<:AbstractFloat

    GFSCAGG! = GFVACGG_sca!;
    if (Genvstr == "PEC")
        GFSCAGG! = GFPECGG!;
    end

    totalalphaeTinvmat = SharedArray{Complex{FT}, 2}(3*myAllMolData.cumnumatomslist[end],
                                                     3*myAllMolData.cumnumatomslist[end]);
    fill!(totalalphaeTinvmat, zero(eltype(totalalphaeTinvmat)));
    integrandvaclist = zeros(FT, myAllMolData.nummols);
    myintegrandvac = zero(FT);
    integrandtot = zero(FT);

    for ll=1:size(freqklist, 1)

        currfreq = 1im * abs(freqklist[ll, 1]);
        currfreqG = currfreq;
        if (nonretarded)
            currfreqG = zero(currfreqG);
        end
        currk = vec(freqklist[ll, 2:4]);
        fill!(totalalphaeTinvmat, zero(eltype(totalalphaeTinvmat)));

        for nn=1:myAllMolData.nummols

            currumolind = myAllMolData.uniquemolinds[nn];
            if (currumolind > -1)

                constructMolalphaealphainv!(myAllMolData.uniqueMolDataArray[currumolind],
                                            myPeriodicData, currfreq, currk);

                constructMolalphaeGvacinf!(myAllMolData.uniqueMolDataArray[currumolind],
                                           myPeriodicData, currfreqG, currk);

                subtract1MolalphainvGvac!(myAllMolData.uniqueMolDataArray[currumolind]);

                integrandvaclist[nn] =
                logabsdet(myAllMolData.uniqueMolDataArray[currumolind].GvacGGinf)[1];

            end

        end

        myintegrandvac = integrandvac(myAllMolData, integrandvaclist);

        for tt=1:myTransData.numTrans

            transformAtomPos!(myTransData, myAllMolData, tt);

            for nn2=1:myAllMolData.nummols

                startidx2 = 1;
                if (nn2 > 1)
                    startidx2 = 3*myAllMolData.cumnumatomslist[nn2-1] + 1;
                end
                currmolind2 = myAllMolData.uniquemolinds[nn2];
                if (currmolind2 < 0)
                    currmolind2 = myAllMolData.dupemolinds[nn2];
                end

                if (myTransData.changedFromBefore[nn2, tt])
                    assemblealphaeTinvDiagBlock!(totalalphaeTinvmat, myAllMolData,
                                                 myTransData, nn2, currmolind2, startidx2,
                                                 tt, myPeriodicData, GFSCAGG!, currfreqG,
                                                 currk);
                end

            end

            for nn2=1:myAllMolData.nummols

                startidx2 = 1;
                if (nn2 > 1)
                    startidx2 = 3*myAllMolData.cumnumatomslist[nn2-1] + 1;
                end

                for nn1=(nn2+1):myAllMolData.nummols

                    startidx1 = 3*myAllMolData.cumnumatomslist[nn1-1] + 1;

                    if (myTransData.changedFromBefore[nn2, tt] ||
                        myTransData.changedFromBefore[nn1, tt])

                        assemblealphaeTinvOffDiagBlock!(totalalphaeTinvmat, myAllMolData,
                                                        myTransData, nn1, nn2, startidx1,
                                                        startidx2, myPeriodicData, GFSCAGG!,
                                                        currfreqG, currk);
                    end
                end

            end

            integrandtot = logabsdet(totalalphaeTinvmat)[1] - myintegrandvac;

            outfile = open(outfilename, "a+");
            if (myPeriodicData.numdims == 0)
                @printf(outfile, "%s %.16e %.16e\n", myTransData.transLabels[tt],
                        imag(currfreq), integrandtot);
            else
                @printf(outfile, "%s %.16e %.16e %.16e %.16e %.16e\n",
                        myTransData.transLabels[tt], imag(currfreq), currk[1], currk[2],
                        currk[3], integrandtot);
                @printf(outfile, "%s %.16e %.16e %.16e %.16e %.16e\n",
                        myTransData.transLabels[tt], imag(currfreq), -1*currk[1],
                        -1*currk[2], -1*currk[3], integrandtot);
            end
            close(outfile);

            untransformAtomPos!(myTransData, myAllMolData, tt);

        end

    end

end

"""

    constructMolalphaealphainv!(myMolData::OneMol{FT},
                                myPeriodicData::PeriodicData{FT},
                                freq::Union{FT, Complex{FT}},
                                k::Array{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

Construct the product of the electronic susceptibility with the full inverse molecular
susceptibility `myMolData.alphainv` in-place at frequency `freq` and wavevector `k`, the
latter of which must be a 3-element real vector. With phonons, `myMolData.alphainv` will
be a matrix, the vector of atomic polarizabilities `myMolData.alpha0` will depend on the
elements of `myMolData.alphainv`, and `myMolData.alphae` will be independent of these two.
Without phonons, all three will be simply related to each other. This is more efficient
than calling `constructMolAlpha!` and then performing a `broadcast` multiplication involving
`myMolData.alphae` multiplying the rows of `myMolData.alphainv`.

"""
function constructMolalphaealphainv!(myMolData::OneMol_NoPhonons{FT},
                                     myPeriodicData::PeriodicData{FT},
                                     freq::Union{FT, Complex{FT}},
                                     k::Array{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

    setindex!(myMolData.alphae,
              (myMolData.Qe.^2) ./ (myMolData.Ke .- (1im.*freq.*myMolData.Be .+
                                freq^2 .* myMolData.Me)),
              1:myMolData.numatoms);
    setindex!(myMolData.alpha0, abs.(myMolData.alphae), 1:myMolData.numatoms);
    fill!(myMolData.alphainv, one(eltype(myMolData.alphainv)));

end
function constructMolalphaealphainv!(myMolData::OneMol_WithPhonons{FT},
                                     myPeriodicData::PeriodicData{FT},
                                     freq::Union{FT, Complex{FT}},
                                     k::Array{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

    Ke3N = vecNto3N(myMolData.Ke);
    Qe3N = vecNto3N(myMolData.Qe);

    setindex!(myMolData.alphainv, KIk(myMolData, myPeriodicData, k),
              1:3*myMolData.numatoms, 1:3*myMolData.numatoms);

    addtodiag!(myMolData.alphainv,
               Ke3N .- vecNto3N(1im.*freq.*myMolData.BI .+ freq^2 .* myMolData.MI));

    setindex!(myMolData.alphainv, inv(myMolData.alphainv),
              1:3*myMolData.numatoms, 1:3*myMolData.numatoms);

    mulMatDiagleft!(myMolData.alphainv, Ke3N);
    mulMatDiagright!(myMolData.alphainv, Ke3N);

    alpha = deepcopy(myMolData.alphainv);
    subtractfromdiag!(alpha,
                      Ke3N .- vecNto3N(1im.*freq.*myMolData.Be .+ freq^2 .* myMolData.Me));

    mulMatDiagright!(alpha, 1 ./ Qe3N);
    mulMatDiagleft!(alpha, 1 ./ Qe3N);

    alpha = inv(alpha);
    alphacarttr = Array{eltype(alpha), 2}(undef, myMolData.numatoms, myMolData.numatoms);
    for qq=1:myMolData.numatoms
        for pp=1:myMolData.numatoms
            alphacarttr[pp,qq] = tr(alpha[3*pp-2:3*pp, 3*qq-2:3*qq]);
        end
    end
    setindex!(myMolData.alpha0, vec(abs.(sum(alphacarttr, dims=2))), 1:myMolData.numatoms);
    lmul!(1/3, myMolData.alpha0);

    mulMatDiagright!(myMolData.alphainv, 1 ./ Qe3N);

    mulMatDiagleft!(myMolData.alphainv,
                    Qe3N  ./
                    (Ke3N .- vecNto3N(1im.*freq.*myMolData.Be .+
                                      freq^2 .* myMolData.Me)));

    subtractfromdiag!(myMolData.alphainv, ones(eltype(myMolData.alphainv),
                                               3*myMolData.numatoms));

    setindex!(myMolData.alphae,
              myMolData.Qe.^2 ./ (myMolData.Ke .-
                                  (1im.*freq.*myMolData.Be .+
                                   freq^2 .* myMolData.Me)),
              1:myMolData.numatoms);

end

"""

    constructMolalphaeGvacinf!(myMolData::OneMol{FT},
                               myPeriodicData::PeriodicData{FT},
                               currfreqG::Union{FT, Complex{FT}},
                               currk::Array{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

Assemble the Green's function interaction matrix of a molecule in vacuum
`myMolData.GvacGGinf`, then multiplies on the left by its electronic polarizabilities
in-place, at frequency `currfreqG` and wavevector `currk` passing `myMolData` and
`myPeriodicData` to the subsidiary functions.

"""
function constructMolalphaeGvacinf!(myMolData::OneMol{FT},
                                    myPeriodicData::PeriodicData{FT},
                                    currfreqG::Union{FT, Complex{FT}},
                                    currk::Array{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

    assemble1MolGvacinf!(myMolData, myPeriodicData, currfreqG, currk);
    mulMatDiagleft!(myMolData.GvacGGinf, vecNto3N(myMolData.alphae));

end

"""

    mulalphaeGBlock!(A::AbstractArray{Complex{FT}, 2},
                     myAllMolData::MolSystem{<:OneMol{FT}},
                     nn1::Integer, startidx1::Integer,
                     startidx2::Integer,
                     endidx2::Integer) where FT<:AbstractFloat

Assuming `A` represents a Green's function interaction matrix, multiply its block
`A[startidx1-1+3*myAllMolData.numatomslist[nn1], startidx2:endidx2]` on the left
by the electronic polarizabilities of molecule `nn1` from `myAllMolData`.

"""
function mulalphaeGBlock!(A::AbstractArray{Complex{FT}, 2},
                          myAllMolData::MolSystem{<:OneMol{FT}},
                          nn1::Integer, startidx1::Integer,
                          startidx2::Integer,
                          endidx2::Integer) where FT<:AbstractFloat

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

    if (size(A, 1) < 3*myAllMolData.numatomslist[nn1])
        error("A (currently with ", size(A, 1), " rows) should have at least ",
              3*myAllMolData.numatomslist[nn1], " rows in mulalphaeGBlock!()");
    end

    if (startidx1 < 1 || startidx1 > size(A, 1) - 3*myAllMolData.numatomslist[nn1] + 1)
        error("mulalphaeGBlock!(): startidx1 (currently ", startidx1,
              ") must be between 1 and ", size(A, 1) - 3*myAllMolData.numatomslist[nn1] + 1,
              " as A has ", size(A, 1), " rows and the block will have ",
              3*myAllMolData.numatomslist[nn1], " rows");
    end
    if (startidx2 < 1 || startidx2 > size(A, 2))
        error("mulalphaeGBlock!(): startidx2 (currently ", startidx2,
              ") must be between 1 and ", size(A, 2),
              " as A has ", size(A, 2), " columns");
    end
    if (endidx2 < 1 || endidx2 > size(A, 2))
        error("mulalphaeGBlock!(): endidx2 (currently ", startidx2,
              ") must be between 1 and ", size(A, 2),
              " as A has ", size(A, 2), " columns");
    end
    if (startidx2 > endidx2)
        error("mulalphaeGBlock!(): startidx2 (currently ", startidx2,
              ") must not be larger than endidx2 (currently ", endidx2, ")");
    end

    tempview = view(A, (startidx1-1).+1:3*myAllMolData.numatomslist[nn1], startidx2:endidx2);
    mulMatDiagleft!(tempview, vecNto3N(myAllMolData.uniqueMolDataArray[currmolind1].alphae));

end

"""

    integrandvac(myAllMolData::MolSystem{<:OneMol{FT}},
                 integrandvaclist::Array{FT, 1}) where FT<:AbstractFloat

Populate the list of vacuum contributions to the vdW integrand `integrandvaclist` with
contributions from duplicate molecules where they exist, then return the sum yielding the
total contribution to the integrand in vacuum.

"""
function integrandvac(myAllMolData::MolSystem{<:OneMol{FT}},
                      integrandvaclist::Array{FT, 1}) where FT<:AbstractFloat

    myintegrandvac = zero(eltype(integrandvaclist));
    for nn=1:myAllMolData.nummols
        if (myAllMolData.dupemolinds[nn] > 0)
            integrandvaclist[nn] = integrandvaclist[myAllMolData.dupemolinds[nn]];
        end
        myintegrandvac += integrandvaclist[nn];
    end

    return myintegrandvac;

end

"""

    assemblealphaeTinvDiagBlock!(totalalphaeTinvmat::SharedArray{Complex{FT}, 2},
                                 myAllMolData::MolSystem{<:OneMol{FT}},
                                 myTransData::TransData{FT}, nn2::Integer,
                                 currmolind2::Integer, startidx2::Integer,
                                 tt::Integer, myPeriodicData::PeriodicData{FT},
                                 GFSCAGG!::Function,
                                 currfreqG::Union{FT, Complex{FT}},
                                 currk::Array{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

Assemble the molecular diagonal block of the scattering Green's function interaction matrix
of molecule `nn2` (whose unique molecular index `currmolind2` is precomputed and passed as
an argument) at frequency `currfreqG` and wavevector `currk` using `GFSCAGG!`, multiply on
the left by the electronic polarizabilities, then add the vacuum contribution which has
already been computed as the electronic polarizabilities multiplied by the inverse
T-operator of molecule `nn2`, performing any required matrix rotations corresponding to
transformation `tt` for molecule `nn2`, and stamp that into
`totalalphaeTinvmat[(startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]), (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2])]`.
Finally, undo any rotation corresponding to transformation `tt` in the inverse T-operator of
molecule `nn2` for reuse with future transformations.

"""
function assemblealphaeTinvDiagBlock!(totalalphaeTinvmat::SharedArray{Complex{FT}, 2},
                                      myAllMolData::MolSystem{<:OneMol{FT}},
                                      myTransData::TransData{FT}, nn2::Integer,
                                      currmolind2::Integer, startidx2::Integer,
                                      tt::Integer, myPeriodicData::PeriodicData{FT},
                                      GFSCAGG!::Function,
                                      currfreqG::Union{FT, Complex{FT}},
                                      currk::Array{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

    tempview = view(totalalphaeTinvmat, (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]),
                    (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]));

    fill!(tempview, zero(eltype(tempview)));

    assembleGscaMolBlock!(tempview, myAllMolData, myTransData, nn2, nn2, 1, 1,
                          myPeriodicData, GFSCAGG!, currfreqG, currk);

    lmul!(-1, tempview);
    mulalphaeGBlock!(tempview, myAllMolData, nn2, 1, 1, 3*myAllMolData.numatomslist[nn2]);

    rotate3Nmat!(myAllMolData.uniqueMolDataArray[currmolind2].GvacGGinf,
                 myTransData.rotArray[3*nn2-2:3*nn2, 3*tt-2:3*tt] *
                 myTransData.baselineRot[3*nn2-2:3*nn2, :]);

    tempview .+= myAllMolData.uniqueMolDataArray[currmolind2].GvacGGinf[:, :];

    rotate3Nmat!(myAllMolData.uniqueMolDataArray[currmolind2].GvacGGinf,
                 transpose(myTransData.rotArray[3*nn2-2:3*nn2, 3*tt-2:3*tt] *
                           myTransData.baselineRot[3*nn2-2:3*nn2, :]));

end

"""

    assemblealphaeTinvOffDiagBlock!(totalalphaeTinvmat::SharedArray{Complex{FT}, 2},
                                    myAllMolData::MolSystem{<:OneMol{FT}},
                                    myTransData::TransData{FT},
                                    nn1::Integer, nn2::Integer,
                                    startidx1::Integer,
                                    startidx2::Integer,
                                    myPeriodicData::PeriodicData{FT},
                                    GFSCAGG!::Function,
                                    currfreqG::Union{FT, Complex{FT}},
                                    currk::Array{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

Assemble the molecular off-diagonal block of the scattering Green's function interaction
matrix of molecules `nn2` as the source and `nn1` as the field (whose respective unique
molecular indices `currmolind1` and `currmolind2` are precomputed and passed as arguments)
at frequency `currfreqG` and wavevector `currk` using `GFSCAGG!`, multiply on
the left by the electronic polarizabilities, add the vacuum contribution, and stamp the
negative of that into
`totalalphaeTinvmat[(startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]), (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2])]`.
Do the same for the adjoint, swapping `nn2` and `nn1` without performing further Green's
function computations.

"""
function assemblealphaeTinvOffDiagBlock!(totalalphaeTinvmat::SharedArray{Complex{FT}, 2},
                                         myAllMolData::MolSystem{<:OneMol{FT}},
                                         myTransData::TransData{FT},
                                         nn1::Integer, nn2::Integer,
                                         startidx1::Integer,
                                         startidx2::Integer,
                                         myPeriodicData::PeriodicData{FT},
                                         GFSCAGG!::Function,
                                         currfreqG::Union{FT, Complex{FT}},
                                         currk::Array{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

    tempview = view(totalalphaeTinvmat, (startidx1-1).+(1:3*myAllMolData.numatomslist[nn1]),
                    (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]));
    fill!(tempview, zero(eltype(tempview)));

    assembleGvacMolBlock!(tempview, myAllMolData, myTransData, nn1, nn2, 1, 1,
                          myPeriodicData, currfreqG, currk);
    assembleGscaMolBlock!(tempview, myAllMolData, myTransData, nn1, nn2, 1, 1,
                          myPeriodicData, GFSCAGG!, currfreqG, currk);
    lmul!(-1, tempview);

    totalalphaeTinvmat[(startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]),
                       (startidx1-1).+(1:3*myAllMolData.numatomslist[nn1])] .=
    adjoint(totalalphaeTinvmat[(startidx1-1).+(1:3*myAllMolData.numatomslist[nn1]),
                               (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2])])

    mulalphaeGBlock!(tempview, myAllMolData, nn1, 1, 1, 3*myAllMolData.numatomslist[nn2]);
    tempview = view(totalalphaeTinvmat, (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]),
                    (startidx1-1).+(1:3*myAllMolData.numatomslist[nn1]));
    mulalphaeGBlock!(tempview, myAllMolData, nn2, 1, 1, 3*myAllMolData.numatomslist[nn1]);

end


text:
using JuMP, GLPK

const no_strat = vec([1.0])

function strat_vec(l::Int64, i::Int64)
    to_return = zeros(l)
    @inbounds to_return[i] = 1.0
    return to_return
end

minmax(R::Matrix{Float64}, m::Int64) =
    @inbounds mapreduce(x -> maximum(R[:, x]), min, 1:m)
maxmin(R::Matrix{Float64}, n::Int64) =
    @inbounds mapreduce(x -> minimum(R[x, :]), max, 1:n)
findminmax(R::Matrix{Float64}, n::Int64) =
    @inbounds strat_vec(n, argmax(map(x -> minimum(R[x, :]), 1:n)))
findmaxmin(R::Matrix{Float64}, m::Int64) =
    @inbounds strat_vec(m, argmin(map(x -> maximum(R[:, x]), 1:m)))

struct NashResult
    payoff::Float64
    row_strategy::Vector{Float64}
    column_strategy::Vector{Float64}
end

function nash(R::Matrix{Float64})
    n, m = size(R)

    # Check if we have to do linear programming
    n == 1 && return NashResult(
        minimum(R), no_strat, strat_vec(m, argmin(R)[2]))
    m == 1 && return NashResult(
        maximum(R), strat_vec(n, argmax(R)[1]), no_strat)
    minmax(R, m) == maxmin(R, n) && return NashResult(
        minmax(R, m), findminmax(R, n), findmaxmin(R, m))

    # Set up model and payoff
    model = direct_model(GLPK.Optimizer())
    @variable(model, z)
    @objective(model, Max, 1.0 * z)

    # Solve for row player
    @variable(model, x[1:n], lower_bound = 0.0)
    @constraint(model, c1, x' * R .>= z)
    @constraint(model, sum(x) == 1.0)

    optimize!(model)

    return NashResult(JuMP.value(z), JuMP.value.(x), vec(shadow_price.(c1)))
end

function SM(state::DynamicState, static_state::StaticState, depth::Int64;
    allow_nothing::Bool = false, allow_overfarming::Bool = false,
    sim_to_end::Bool = false)
    A, B = get_possible_decisions(state, static_state,
        allow_nothing = allow_nothing, allow_overfarming = allow_overfarming)

    (Base.ctpop_int(A) == 0x00 || Base.ctpop_int(B) == 0x00 || depth == 0) &&
        return sim_to_end ?
            NashResult(sum(battle_scores(state, static_state, 100)
                / 100) - 0.5, no_strat, no_strat) :
            NashResult(battle_score(state, static_state) - 0.5,
                no_strat, no_strat)

    payoffs = zeros(Float64, Base.ctpop_int(A), Base.ctpop_int(B))

    for i = 0x01:Base.ctpop_int(A), j = 0x01:Base.ctpop_int(B)
        turn_output = play_turn(state, static_state, get_decisions(A, B, i, j))
        if turn_output.odds == 1.0
            @inbounds payoffs[i, j] = SM(turn_output.next_state_1,
                static_state, depth - 1,
                allow_nothing = allow_nothing,
                allow_overfarming = allow_overfarming,
                sim_to_end = sim_to_end).payoff
        else
            @inbounds payoffs[i, j] = turn_output.odds *
                SM(turn_output.next_state_1,
                    static_state, depth - 1,
                    allow_nothing = allow_nothing,
                    allow_overfarming = allow_overfarming,
                    sim_to_end = sim_to_end).payoff + (1 - turn_output.odds) *
                SM(turn_output.next_state_2, static_state, depth - 1,
                    allow_nothing = allow_nothing,
                    allow_overfarming = allow_overfarming,
                    sim_to_end = sim_to_end).payoff
        end
    end
    return nash(payoffs)
end

function solve_battle(s::DynamicState, static_s::StaticState, depth::Int64;
    allow_nothing::Bool = false, allow_overfarming::Bool = false,
    sim_to_end::Bool = false)
    value = 0.0
    decision = 0, 0
    strat = Strategy([], [], [], [])
    while true
        A, B = get_possible_decisions(s, static_s,
            allow_nothing = allow_nothing,
            allow_overfarming = allow_overfarming)

        (Base.ctpop_int(A) == 0x00 || Base.ctpop_int(B) == 0x00) &&
            return value, strat
        if Base.ctpop_int(A) == 0x01 && Base.ctpop_int(B) == 0x01
            decision = get_decisions(A, B, 0x01, 0x01)
        else
            nash_result = SM(s, static_s, depth, sim_to_end = sim_to_end)
            value = nash_result.payoff
            d1, d2 = rand(rb_rng), rand(rb_rng)
            decision1, decision2 = UInt8(length(nash_result.row_strategy)),
                UInt8(length(nash_result.column_strategy))
            j = 0.0
            for i = 0x01:(UInt8(length(nash_result.row_strategy)) - 0x01)
                @inbounds j += nash_result.row_strategy[i]
                if d1 < j
                    decision1 = i
                    break
                end
            end
            j = 0.0
            for i = 0x01:(UInt8(length(nash_result.column_strategy)) - 0x01)
                @inbounds j += nash_result.column_strategy[i]
                if d2 < j
                    decision2 = i
                    break
                end
            end
            decision = get_decisions(A, B, decision1, decision2)
        end
        s = select_next_state(play_turn(s, static_s, decision))
        push!(strat.decisions, decision)
        push!(strat.scores, value + 0.5)
        push!(strat.activeMons, get_active(s))
        push!(strat.hps, ((get_hp(s[0x01][0x01]), get_hp(s[0x01][0x02]),
                           get_hp(s[0x01][0x03])), (get_hp(s[0x02][0x01]),
                           get_hp(s[0x02][0x02]), get_hp(s[0x02][0x03]))))
    end
end

ext: jl


text:
using LoggedDicts
using Base.Test

if !isfile("ld_test.log")
    run(`touch $(joinpath(dirname(@__FILE__), "ld_test.log"))`)
end
if !isfile("loggeddict.test")
    run(`touch $(joinpath(dirname(@__FILE__), "loggeddict.test"))`)
end

# Create and populate a LoggedDict
ld = LoggedDict("my_ld", "ld_test.log")      # Init empty LoggedDict, with logs entries recorded in "mydict.log"
set!(ld, "key1", "some_value")
set!(ld, "key2", 2)                          # ld["key2"] equals 2
@test get(ld, "key2") == 2
set!(ld, "key2", "key21", [1.2, 3.4])        # ld["key2"] equals Dict("key21" => [1.2, 3.4]), overwrites previous value of 2
set!(ld, "key3", "key31", Set([1, 2, 3]))    # ld["key3"] equals Dict("key31" => Set([1, 2, 3]))
set!(ld, "key3", "key32", 32)                # ld["key3"] equals Dict("key31" => Set([1, 2, 3]), "key32" => 32)
pop!(ld, "key3", "key31", 2)                 # ld["key3"] equals Dict("key31" => Set([1, 3]), "key32" => 32)
push!(ld, "key3", "key31", 4)                # ld["key3"] equals Dict("key31" => Set([1, 3, 4]), "key32" => 32)
set!(ld, "key4", Dict("key41" => 41, "key42" => 42))    # ld["key4"] equals Dict("key41" => 4, "key42" => 42)

# Test some values
@test get(ld, "key1") == "some_value"
@test get(ld, "key2") == Dict("key21" => [1.2, 3.4])
@test get(ld, "key3") == Dict("key31" => Set([1, 3, 4]), "key32" => 32)
@test get(ld, "key4") == Dict("key41" => 41, "key42" => 42)
@test_throws ErrorException get(ld, "key5", "key51")

# Test: delete! and haskey
delete!(ld, "key1")
delete!(ld, "key3", "key32")
@test haskey(ld, "key1") == false
@test haskey(ld, "key3") == true
@test haskey(ld, "key3", "key32") == false

# Write the LoggedDict to disk
write_logged_dict("loggeddict.test", ld)

# Read the LoggedDict from disk
ld = ""
ld = read_logged_dict("loggeddict.test")

# Retest the values
@test get(ld, "key2") == Dict("key21" => [1.2, 3.4])
@test get(ld, "key3") == Dict("key31" => Set([1, 3, 4]))
@test get(ld, "key4") == Dict("key41" => 41, "key42" => 42)

# Test LoggedDict with logging off
ld = LoggedDict("my_ld", "ld_test.log", true)      # Init empty LoggedDict, with logs entries recorded in "mydict.log"
set!(ld, "key1", "some_value")
@test get(ld, "key1") == "some_value"

### EOF


text:
using PupperSim
using Test

@testset "PupperSim.jl" begin
    @test size(names(PupperSim)) == (4,)
    @test names(PupperSim)[1] == :PupperSim
    @test names(PupperSim)[2] == :loadmodel
    @test names(PupperSim)[3] == :pupper
    @test names(PupperSim)[4] == :simulate
end

@testset "image_flip.jl" begin
A = [1 2 3; 4 5 6; 7 8 9]
PupperSim.hflip!(A)
@test A == [3 2 1; 6 5 4; 9 8 7]
PupperSim.vflip!(A)
@test A == [9 8 7; 6 5 4; 3 2 1]
end

@testset "pupper" begin
    robot = pupper()
    @test string(robot.command) == "QuadrupedController.Command([0.4, 0.0], 0.0, -0.06, 0.1, 0.0, 0, false, false, false)"
end

@testset "PupperSim.RobotCmd" begin
    @test string(PupperSim.TURN_LEFT) == "TURN_LEFT"
    @test Int(PupperSim.TURN_LEFT) == 13
end

using GLFW
const GLFW_MOD_CAPS_LOCK = 0x0010   # Caps Lock key is enabled
const GLFW_MOD_NUM_LOCK  = 0x0020   # Num Lock key is enabled

@testset "PupperSim.keypadcmd" begin
    @test PupperSim.keypadcmd(GLFW.KEY_KP_ADD, Int32(0))                    == PupperSim.CYCLE_HOP
    @test PupperSim.keypadcmd(GLFW.KEY_KP_ADD, Int32(GLFW.MOD_SHIFT))       == PupperSim.CYCLE_HOP
    @test PupperSim.keypadcmd(GLFW.KEY_KP_ADD, Int32(GLFW.MOD_CONTROL))     == PupperSim.NO_COMMAND
    @test PupperSim.keypadcmd(GLFW.KEY_KP_ADD, Int32(GLFW.MOD_ALT))         == PupperSim.NO_COMMAND
    @test PupperSim.keypadcmd(GLFW.KEY_KP_ADD, Int32(GLFW.MOD_SUPER))       == PupperSim.CYCLE_HOP
    @test PupperSim.keypadcmd(GLFW.KEY_KP_ADD, Int32(GLFW_MOD_CAPS_LOCK))   == PupperSim.CYCLE_HOP
    @test PupperSim.keypadcmd(GLFW.KEY_KP_ADD, Int32(GLFW_MOD_NUM_LOCK))    == PupperSim.CYCLE_HOP

    @test PupperSim.keypadcmd(GLFW.KEY_KP_7, Int32(0))                      == PupperSim.INCREASE_HEIGHT
    @test PupperSim.keypadcmd(GLFW.KEY_KP_7, Int32(GLFW.MOD_SHIFT))         == PupperSim.NO_COMMAND
    @test PupperSim.keypadcmd(GLFW.KEY_KP_7, Int32(GLFW.MOD_CONTROL))       == PupperSim.NO_COMMAND
    @test PupperSim.keypadcmd(GLFW.KEY_KP_7, Int32(GLFW.MOD_ALT))           == PupperSim.NO_COMMAND
    @test PupperSim.keypadcmd(GLFW.KEY_KP_7, Int32(GLFW.MOD_SUPER))         == PupperSim.NO_COMMAND
    @test PupperSim.keypadcmd(GLFW.KEY_KP_7, Int32(GLFW_MOD_CAPS_LOCK))     == PupperSim.INCREASE_HEIGHT
    @test PupperSim.keypadcmd(GLFW.KEY_KP_7, Int32(GLFW_MOD_NUM_LOCK))      == PupperSim.NO_COMMAND
end

@testset "PupperSim.keyboardcmd" begin
    @test PupperSim.keyboardcmd(GLFW.KEY_PAGE_UP, Int32(0))                 == PupperSim.NO_COMMAND
    @test PupperSim.keyboardcmd(GLFW.KEY_PAGE_UP, Int32(GLFW.MOD_SHIFT))    == PupperSim.INCREASE_VELOCITY
    @test PupperSim.keyboardcmd(GLFW.KEY_PAGE_UP, Int32(GLFW.MOD_CONTROL))  == PupperSim.NO_COMMAND
    @test PupperSim.keyboardcmd(GLFW.KEY_PAGE_UP, Int32(GLFW.MOD_ALT))      == PupperSim.NO_COMMAND
    @test PupperSim.keyboardcmd(GLFW.KEY_PAGE_UP, Int32(GLFW.MOD_SUPER))    == PupperSim.NO_COMMAND
    @test PupperSim.keyboardcmd(GLFW.KEY_PAGE_UP, Int32(GLFW_MOD_CAPS_LOCK))== PupperSim.NO_COMMAND
    @test PupperSim.keyboardcmd(GLFW.KEY_PAGE_UP, Int32(GLFW_MOD_NUM_LOCK)) == PupperSim.NO_COMMAND

    @test PupperSim.keyboardcmd(GLFW.KEY_PERIOD, Int32(0))                  == PupperSim.NO_COMMAND
    @test PupperSim.keyboardcmd(GLFW.KEY_PERIOD, Int32(GLFW.MOD_SHIFT))     == PupperSim.TURN_RIGHT
    @test PupperSim.keyboardcmd(GLFW.KEY_PERIOD, Int32(GLFW.MOD_CONTROL))   == PupperSim.NO_COMMAND
    @test PupperSim.keyboardcmd(GLFW.KEY_PERIOD, Int32(GLFW.MOD_ALT))       == PupperSim.NO_COMMAND
    @test PupperSim.keyboardcmd(GLFW.KEY_PERIOD, Int32(GLFW.MOD_SUPER))     == PupperSim.NO_COMMAND
    @test PupperSim.keyboardcmd(GLFW.KEY_PERIOD, Int32(GLFW_MOD_CAPS_LOCK)) == PupperSim.NO_COMMAND
    @test PupperSim.keyboardcmd(GLFW.KEY_PERIOD, Int32(GLFW_MOD_NUM_LOCK))  == PupperSim.NO_COMMAND
end


text:
struct VectorOfBlockArrayCoo{T,N,B,X,Z} <: AbstractVector{T}
  blocks::B
  blockids::Vector{NTuple{N,Int}}
  axes::X
  ptrs::Array{Int,N}
  zero_blocks::Z
  function VectorOfBlockArrayCoo(
    blocks::Tuple,
    blockids::Vector{NTuple{N,Int}},
    axes::AbstractArray{<:NTuple{N}},
    ptrs::Array{Int,N} = _compute_ptrs(blockids,testitem(axes)),
    zero_blocks::Tuple=_compute_zero_blocks_array(blocks,ptrs,axes)) where N

    msg = "Trying to build a VectorOfBlockArrayCoo with repeated blocks"
    @assert _no_has_repeaded_blocks(blockids,ptrs) msg

    if length(blocks) != 0
      @assert all(map(length,blocks) .==  length(first(blocks)) )
      blocks_i = collect(map(testitem,blocks))
      zero_blocks_i = collect(eltype(blocks_i),map(testitem,zero_blocks))
    else
      @assert length(zero_blocks) !=0
      zero_blocks_i = collect(map(testitem,zero_blocks))
      blocks_i = collect(eltype(zero_blocks_i),map(testitem,blocks))
    end

    B = typeof(blocks)
    X = typeof(axes)
    axes_i = testitem(axes)
    t = BlockArrayCoo(blocks_i,blockids,axes_i,ptrs,zero_blocks_i)
    T = typeof(t)
    Z = typeof(zero_blocks)
    new{T,N,B,X,Z}(blocks,blockids,axes,ptrs,zero_blocks)
  end
end

function VectorOfBlockArrayCoo(
  blocks::Tuple,
  blockids::Vector{NTuple{N,Int}},
  axes::AbstractArray{<:NTuple{N}},
  zero_blocks::Tuple) where N

  ptrs = _compute_ptrs(blockids,testitem(axes))
  VectorOfBlockArrayCoo(blocks,blockids,axes,ptrs,zero_blocks)
end

const VectorOfBlockVectorCoo = VectorOfBlockArrayCoo{T,1} where T
const VectorOfBlockMatrixCoo = VectorOfBlockArrayCoo{T,2} where T

function _no_has_repeaded_blocks(blockids::Vector{NTuple{N,Int}},ptrs) where N
  maxblocks = size(ptrs)
  touched = zeros(Int,maxblocks)
  for c in blockids
    touched[c...] += 1
  end
  all( touched .<= 1 )
end

function _compute_zero_blocks_array(blocks,ptrs,axs)
  A = eltype(first(blocks))
  cis = CartesianIndices(ptrs)
  zero_blocks = []
  for ci in cis
    p = ptrs[ci]
    if p<0
      block = _zeros_like(A,Tuple(ci),axs)
      push!(zero_blocks,block)
    end
  end
  Tuple(zero_blocks)
end

function zeros_like(a::VectorOfBlockArrayCoo{T,N} where T) where N
  ptrs = copy(a.ptrs)
  for k in 1:length(ptrs)
    ptrs[k] = -k
  end
  zero_blocks = _compute_zero_blocks_array(a.blocks,ptrs,a.axes)
  blocks = ()
  blockids = NTuple{N,Int}[]
  VectorOfBlockArrayCoo(blocks,blockids,a.axes,ptrs,zero_blocks)
end

function _zeros_like(::Type{A},I,axs) where A
  block = apply(axs) do a
    laxs = map( local_range, a, I)
    zeros_like(A,laxs)
  end
  block
end

function _zeros_like(::Type{<:BlockArrayCoo{T,N,A} where T},I,axs) where {N,A}
  axsI = apply(axs) do a
    map(local_range, a, I)
  end
  blocks = ()
  blockids = NTuple{N,Int}[]
  ptrs = zeros(Int,map(x->blocksize(x,1),testitem(axsI)))
  for k in 1:length(ptrs)
    ptrs[k] = -k
  end
  zero_blocks = []
  cis = CartesianIndices(ptrs)
  for ci in cis
    block = _zeros_like(A,Tuple(ci),axsI)
    push!(zero_blocks,block)
  end
  VectorOfBlockArrayCoo(blocks,blockids,axsI,ptrs,Tuple(zero_blocks))
end

BlockArrays.blocksize(a::VectorOfBlockArrayCoo) = blocksize(testitem(a))

BlockArrays.blocksize(a::VectorOfBlockArrayCoo,i::Integer) = blocksize(testitem(a),i)

function BlockArrays.eachblock(a::VectorOfBlockArrayCoo)
  cis = CartesianIndices(blocksize(a))
  blocks = map(ci->Block(Tuple(ci)),cis)
  ( a[block] for block in blocks  )
end

Base.IndexStyle(::Type{<:VectorOfBlockArrayCoo}) = IndexLinear()

function Base.size(a::VectorOfBlockArrayCoo)
  if length(a.blocks) != 0
    (length(first(a.blocks)),)
  else
    (length(first(a.zero_blocks)),)
  end
end

function array_cache(a::VectorOfBlockArrayCoo)

  if length(a.blocks) != 0
    blocks_i = collect(map(testitem,a.blocks))
    zero_blocks_i = collect(eltype(blocks_i),map(testitem,a.zero_blocks))
  else
    zero_blocks_i = collect(map(testitem,a.zero_blocks))
    blocks_i = collect(eltype(zero_blocks_i),map(testitem,a.blocks))
  end

  ca = array_cache(a.axes)
  cb = array_caches(a.blocks...)
  cz = array_caches(a.zero_blocks...)
  (blocks_i,zero_blocks_i,ca,cb,cz)
end

@inline function getindex!(cache,a::VectorOfBlockArrayCoo,i::Integer)
  blocks_i, zero_blocks_i, ca, cb, cz = cache
  axes_i = getindex!(ca,a.axes,i)
  blocks_i .= getitems!(cb,a.blocks,i)
  zero_blocks_i .= getitems!(cz,a.zero_blocks,i)
  BlockArrayCoo(blocks_i,a.blockids,axes_i,a.ptrs,zero_blocks_i)
end

function Base.getindex(a::VectorOfBlockArrayCoo,i::Integer)
  cache = array_cache(a)
  getindex!(cache,a,i)
end

function Base.getindex(a::VectorOfBlockArrayCoo,b::Block)
  _get_block_index(a,b)
end

function Base.getindex(a::VectorOfBlockVectorCoo,b::Block{1})
  _get_block_index(a,b)
end

function _get_block_index(a,b)
  i = convert(Tuple,b)
  p = a.ptrs[i...]
  if p>0
    a.blocks[p]
  else
    a.zero_blocks[-p]
  end
end

function Base.getindex(a::VectorOfBlockArrayCoo,b::Block{1}...)
  a[Block(map(Int,b)...)]
end

function is_zero_block(a::VectorOfBlockArrayCoo,i::Integer...)
  p = a.ptrs[i...]
  @assert p != 0
  p < 0
end

function apply(::typeof(transpose),a::VectorOfBlockMatrixCoo)
  blocks = [ apply(transpose,block) for block in a.blocks ]
  zero_blocks = [ apply(transpose,block) for block in a.zero_blocks ]
  blockids = [ (j,i) for (i,j) in a.blockids ]
  axs = apply( ax->(ax[2],ax[1]), a.axes)
  ptrs = collect(Transpose(a.ptrs))
  VectorOfBlockArrayCoo(Tuple(blocks),blockids,axs,ptrs,Tuple(zero_blocks))
end


text:
using ImgHdr
using ImgHdr: imgtype
using Test: @testset, @test, @test_throws

@testset "test_images_types" begin
    png = "images/example.png"
    gif = "images/example.gif"

    @test ImgHdr.isa_png(png)
    @test imgtype(png) == "png"
    @test imgtype(png) != "gif"

    @test ImgHdr.isa_gif(gif)
    @test imgtype(gif) == "gif"
    @test imgtype(gif) != "png"

    jpeg = "images/example.jpeg"
    @test ImgHdr.isa_jpeg(jpeg)
    @test imgtype(jpeg) == "jpeg"
    @test imgtype(jpeg) != "pgn"

    exr = "images/example.exr"
    @test ImgHdr.isa_exr(exr)
    @test imgtype(exr) == "exr"
    @test imgtype(exr) != "gif"

    tiff = "images/example.tiff"
    @test ImgHdr.isa_tiff(tiff)
    @test imgtype(tiff) == "tiff"
    @test imgtype(tiff) != "png"

    webp = "images/example.webp"
    @test ImgHdr.isa_webp(webp)
    @test imgtype(webp) == "webp"
    @test imgtype(webp) != "gif"

    pbm = "images/example.pbm"
    @test ImgHdr.isa_pbm(pbm)
    @test imgtype(pbm) == "pbm"
    @test imgtype(pbm) != "png"

    xbm = "images/example.xbm"
    @test ImgHdr.isa_xbm(xbm)
    @test imgtype(xbm) == "xbm"
    @test imgtype(xbm) != "gif"

    bmp = "images/example.bmp"
    @test ImgHdr.isa_bmp(bmp)
    @test imgtype(bmp) == "bmp"
    @test imgtype(bmp) != "png"

    pgm = "images/example.pgm"
    @test ImgHdr.isa_pgm(pgm)
    @test imgtype(pgm) == "pgm"
    @test imgtype(pgm) != "gif"

    ppm = "images/example.ppm"
    @test ImgHdr.isa_ppm(ppm)
    @test imgtype(ppm) == "ppm"
    @test imgtype(ppm) != "png"

    rast = "images/example.ras"
    @test ImgHdr.isa_rast(rast)
    @test imgtype(rast) == "rast"
    @test imgtype(rast) != "gif"

    rgb = "images/example.rgb"
    @test ImgHdr.isa_rgb(rgb)
    @test imgtype(rgb) == "rgb"
    @test imgtype(rgb) != "png"

    not = "images/not-a-image.txt"
    @test ImgHdr.isa_ppm(not) == false
    @test imgtype(not) == nothing
end


text:
module Bandits
import Base: >=, <=, +, /

export
    Bandit, StaticBandit, staticbandit, pull, regret,
    Policy, GreedyPolicy, EpsilonGreedyPolicy, ExploreThenExploit, UCB1, ThompsonSampling, choose,
    Agent, BasicAgent, BetaBernoulliAgent, NormalAgent, randomSample, choose, observe,
    Action, BanditStats, simulate, aggregate_simulate

include("utils.jl")
include("actions.jl")
include("bandit.jl")
include("policy.jl")
include("agent.jl")
include("simulate.jl")

end


text:
export energy, train, get_gradients
using Flux: Optimise

function get_gradients(qpeps::QPEPSMachine, model)
    return map(qpeps.runtime.rotblocks) do r
        r.theta += π/2
        E₊ = energy(qpeps, model)
        r.theta -= π
        E₋ = energy(qpeps, model)
        r.theta += π/2
        0.5(E₊ - E₋)
    end
end

function train(config, model; maxiter=200, optimizer=Optimise.ADAM(0.1), nbatch=1024, use_cuda=true)
    @assert nspins(config) == nspins(model)
    #@assert config.nv+config.nrepeat == model.size[1]
    reg0 = zero_state(nqubits(config); nbatch=nbatch)
    use_cuda && (reg0 = reg0 |> cu)
    qpeps = QPEPSMachine(config, reg0)
    circuit = qpeps.runtime.circuit
    rotblocks = qpeps.runtime.rotblocks
    dispatch!(circuit, :random)
    println("E0 = $(energy(qpeps, model))")
    flush(stdout)

    history = Float64[]
    params = parameters(circuit)
    println("Number of parameters is $(length(params))")
    for i in 1:maxiter
        grad = get_gradients(qpeps, model)
        Optimise.update!(optimizer, params, grad)
        dispatch!.(rotblocks, params)
        push!(history, energy(qpeps, model))
        println("Iter $i, E = $(history[end])")
        flush(stdout)
    end
    qpeps, history
end


text:
# -----------------------------------------------
# Total Welfare
# -----------------------------------------------

@defcomp fund_welfare begin

    regions  = Index()                                   # 16 regions in FUND.

    ρ                 = Parameter()                      # Pure rate of time preference.
    η                 = Parameter()                      # Elasticity of marginal utility of consumption.
    welfare_year      = Parameter{Int}()                      # First year to begin counting total discounted welfare over time.
    pop               = Parameter(index=[time, regions]) # Regional population levels (millions of people).
    c                 = Parameter(index=[time, regions]) # Regional total consumption levels ($).

    UTILITY           = Variable(index=[time])           # Total economic welfare over model time horizon (using the same name as the default RICE welfare component).
    cpc               = Variable(index=[time, regions])  # Per capita consumption levels ($).

    function run_timestep(p, v, d, t)

        # Calculate regional per capita consumption levels.
        for r in d.regions
            v.cpc[t,r] = p.c[t,r] / p.pop[t,r] / 1000000.0
        end

        # Calculate total welfare for each period. If η = 1, the social welfare function becomes log(x).
        # Note: FUND starts a spin-up period in 1950 to initialize the model, so should not track welfare in early periods.
        if gettime(t) < p.welfare_year
            v.UTILITY[t] = 0.0
        else
            # Calculate total welfare over time.
            if p.η == 1.0
                v.UTILITY[t] = v.UTILITY[t-1] + sum(log.(v.cpc[t,:]) .* p.pop[t,:]) / (1.0 + p.ρ)^(gettime(t) - p.welfare_year)
            else
                v.UTILITY[t] = v.UTILITY[t-1] + sum((v.cpc[t,:] .^ (1.0 - p.η)) ./ (1.0 - p.η) .* p.pop[t,:]) / (1.0 + p.ρ)^(gettime(t) - p.welfare_year)
            end
        end
    end
end

max_stars_repo_name: Environment-Research/Utilitarianism

max_issues_repo_name: Environment-Research/Utilitarianism

max_forks_repo_name: Environment-Research/Utilitarianism


text:
#! /usr/bin/julia

# Rosetta Code, Non-continuous subsequences

iscontseq(n::Integer) = count_zeros(n) == leading_zeros(n) + trailing_zeros(n)
iscontseq(n::BigInt) = !ismatch(r"0", rstrip(bin(n), '0'))

function makeint2seq(n::Integer)
    const idex = collect(1:n)
    function int2seq(m::Integer)
        d = digits(m, 2, n)
        idex[d .== 1]
    end
    return int2seq
end

immutable NCSubSeq{T<:Integer}
    n::T
end

type NCSubState{T<:Integer}
    m::T
    m2s::Function
end

Base.length(a::NCSubSeq) = 2^a.n - div(a.n*(a.n+1), 2) - 1
Base.start(a::NCSubSeq) = NCSubState(5, makeint2seq(a.n))
Base.done(a::NCSubSeq, as::NCSubState) = 2^a.n-3 < as.m

function Base.next(a::NCSubSeq, as::NCSubState)
    s = as.m2s(as.m)
    as.m += 1
    while iscontseq(as.m)
        as.m += 1
    end
    return (s, as)
end

n = 4
print("Testing NCSubSeq for ", n, " items:\n   ")
for a in NCSubSeq(n)
    print(" ", a)
end
println()

s = "Rosetta"
cs = split(s, "")
m = 10
n = length(NCSubSeq(length(s))) - m
println()
println("The first and last ", m, " NC sub-sequences of \"", s, "\":")
for (i,a) in enumerate(NCSubSeq(length(cs)))
    i <= m || n < i || continue
    println(@sprintf "%6d %s" i join(cs[a], ""))
    i == m || continue
    println("    .. ......")
end

t = {}
append!(t, collect(1:10))
append!(t, collect(20:10:40))
append!(t, big(50):50:200)
println()
println("Numbers of NC sub-sequences of a given length:")
for i in t
    println(@sprintf("%7d => ", i), length(NCSubSeq(i)))
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
##############################################################################
##
## Partial
## http://chairnerd.seatgeek.com/fuzzywuzzy-fuzzy-string-matching-in-python/
##
##############################################################################
type Partial{T <: PreMetric} <: PreMetric
    dist::T
end

# general
function compare(dist::Partial, s1::AbstractStringorGraphemeIterator, s2::AbstractStringorGraphemeIterator, len1::Integer, len2::Integer)
    len1 == len2 && return compare(dist.dist, s1, s2, len1, len2)
    len1 == 0 && return compare(dist.dist, "", "", 0, 0)
    iter = QGramIterator(s2, len2, len1)
    state = start(iter)
    s, state = next(iter, state)
    out = compare(dist.dist, s1, s)
    while !done(iter, state)
        s, state = next(iter, state)
        curr = compare(dist.dist, s1, s)
        out = max(out, curr)
    end
    return out
end

# Specialization for RatcliffObershelp distance
# Code follows https://github.com/seatgeek/fuzzywuzzy/blob/master/fuzzywuzzy/fuzz.py
function compare(dist::Partial{RatcliffObershelp}, s1::AbstractStringorGraphemeIterator, s2::AbstractStringorGraphemeIterator, len1::Integer, len2::Integer)
    len1 == len2 && return compare(dist.dist, s1, s2, len1, len2)
    out = 0.0
    result = matching_blocks(s1, s2)
    for r in result
        # here I difffer from fuzz.py by making sure the substring of s2 has length len1
        s2_start = r[2] - r[1] + 1
        s2_end = s2_start + len1 - 1
        if s2_start <= 0
            s2_end += 1 - s2_start
            s2_start += 1 - s2_start
        elseif s2_end > len2
            s2_start += len2 - s2_end
            s2_end += len2 - s2_end
        end
        i2_start =  chr2ind(s2, s2_start)
        i2_end = s2_end == len2 ? endof(s2) : (chr2ind(s2, s2_end + 1) - 1)
        curr = compare(RatcliffObershelp(), s1, SubString(s2, i2_start, i2_end), len1, len1)
        out = max(out, curr)
    end
    return out
end

##############################################################################
##
## TokenSort
## http://chairnerd.seatgeek.com/fuzzywuzzy-fuzzy-string-matching-in-python/
##
##############################################################################
type TokenSort{T <: PreMetric} <: PreMetric
    dist::T
end

function compare(dist::TokenSort, s1::AbstractStringorGraphemeIterator, s2::AbstractStringorGraphemeIterator,
    len1::Integer, len2::Integer)
    if search(s1, Base._default_delims) > 0
        s1 = iterator(typeof(s1), join(sort!(split(s1)), " "))
    end
    if search(s2, Base._default_delims) > 0
        s2 = iterator(typeof(s2), join(sort!(split(s2)), " "))
    end
    @show s1, s2
    compare(dist.dist, s1, s2)
end

##############################################################################
##
## TokenSet
## http://chairnerd.seatgeek.com/fuzzywuzzy-fuzzy-string-matching-in-python/
##
##############################################################################
type TokenSet{T <: PreMetric} <: PreMetric
    dist::T
end

function compare(dist::TokenSet, s1::AbstractStringorGraphemeIterator, s2::AbstractStringorGraphemeIterator,
    len1::Integer, len2::Integer)
    v0, v1, v2 = _separate!(split(s1), split(s2))
    s0 = iterator(typeof(s1), join(v0, " "))
    s1 = iterator(typeof(s1), join(chain(v0, v1), " "))
    s2 = iterator(typeof(s2), join(chain(v0, v2), " "))
    if isempty(s0)
        # otherwise compare(dist, "", "a")== 1.0
        compare(dist.dist, s1, s2)
    else
        max(compare(dist.dist, s0, s1),
            compare(dist.dist, s1, s2),
            compare(dist.dist, s0, s2))
    end
end

# separate 2 vectors in intersection, setdiff1, setdiff2 (all sorted)
function _separate!(v1::Vector, v2::Vector)
    sort!(v1)
    sort!(v2)
    out = eltype(v1)[]
    start = 1
    i1 = 0
    while i1 < length(v1)
        i1 += 1
        x = v1[i1]
        i2 = searchsortedfirst(v2, x, start, length(v2), Base.Forward)
        i2 > length(v2) && break
        if i2 > 0 && v2[i2] == x
            deleteat!(v1, i1)
            deleteat!(v2, i2)
            push!(out, x)
            i1 -= 1
            start = i2
        end
    end
    return out, v1, v2
end

##############################################################################
##
## TokenMax
##
##############################################################################
type TokenMax{T <: PreMetric} <: PreMetric
    dist::T
end

function compare(dist::TokenMax, s1::AbstractStringorGraphemeIterator, s2::AbstractStringorGraphemeIterator,
    len1::Integer, len2::Integer)
    base = compare(dist.dist, s1, s2, len1, len2)
    unbase_scale = 0.95
    # if one string is much much shorter than the other
    if len2 >= 1.5 * len1
        # if strings are of dissimilar length, use partials
        partial = compare(Partial(dist.dist), s1, s2, len1, len2)
        ptsor = compare(TokenSort(Partial(dist.dist)), s1, s2, len1, len2)
        ptser = compare(TokenSet(Partial(dist.dist)), s1, s2, len1, len2)
        partial_scale = len2 > (8 * len1) ? 0.6 : 0.9
        return max(base,
                partial * partial_scale,
                ptsor * unbase_scale * partial_scale,
                ptser * unbase_scale * partial_scale)
    else
        ptsor = compare(TokenSort(dist.dist), s1, s2, len1, len2)
        ptser = compare(TokenSet(dist.dist), s1, s2, len1, len2)
        return max(base,
                ptsor * unbase_scale,
                ptser * unbase_scale)
    end
end


text:
module CoordsModule

using UserExceptionModule: UserError
using SpecialDirectionCodesModule

abstract type Coords end

function is_ocean(coords::Coords)
  throw(UserError())
end

function is_outflow(coords::Coords)
  throw(UserError())
end

function is_truesink(coords::Coords)
  throw(UserError())
end

function is_lake(coords::Coords)
  throw(UserError())
end

function are_valid(coords::Coords)
  throw(UserError())
end

struct LatLonCoords <: Coords
  lat::Int64
  lon::Int64
end

struct Generic1DCoords <: Coords
  index::Int64
end

function is_ocean(coords::LatLonCoords)
  return coords.lat == coord_base_indicator_ocean &&
         coords.lon == coord_base_indicator_ocean
end

function is_outflow(coords::LatLonCoords)
  return coords.lat == coord_base_indicator_outflow &&
         coords.lon == coord_base_indicator_outflow
end

function is_truesink(coords::LatLonCoords)
  return coords.lat == coord_base_indicator_truesink &&
         coords.lon == coord_base_indicator_truesink
end

function is_lake(coords::LatLonCoords)
  return coords.lat == coord_base_indicator_lake &&
         coords.lon == coord_base_indicator_lake
end

function is_ocean(coords::Generic1DCoords)
  return coords.index == coord_base_indicator_ocean
end

function is_outflow(coords::Generic1DCoords)
  return coords.index == coord_base_indicator_outflow
end

function is_truesink(coords::Generic1DCoords)
  return coords.index == coord_base_indicator_truesink
end

function is_lake(coords::Generic1DCoords)
  return coords.index == coord_base_indicator_lake
end

struct DirectionIndicator
  coords::Coords
end

function get_next_cell_coords(flow_direction::DirectionIndicator,
                              coords::Coords)
  return flow_direction.coords
end

function is_ocean(flow_direction::DirectionIndicator)
  return is_ocean(flow_direction.coords)
end

function is_outflow(flow_direction::DirectionIndicator)
  return is_outflow(flow_direction.coords)
end

function is_truesink(flow_direction::DirectionIndicator)
  return is_truesink(flow_direction.coords)
end

function is_lake(flow_direction::DirectionIndicator)
  return is_lake(flow_direction.coords)
end

abstract type SectionCoords end

struct LatLonSectionCoords <: SectionCoords
  min_lat::Int64
  max_lat::Int64
  min_lon::Int64
  max_lon::Int64
end

end

max_stars_repo_name: ThomasRiddick/DynamicHD

max_issues_repo_name: ThomasRiddick/DynamicHD

max_forks_repo_name: ThomasRiddick/DynamicHD


text:
using Distributed
using JLD2

include("../BioClim/src/required.jl")
include("expandlayers_vulpes_vulpes.jl")

#find probability of species occuring together
#create array, probs, using pred data
probs = copy(newlayers[1].grid)
#for all prey data
@time for j in 2:length(newlayers)     # iterate through all preys
    for k in 1:length(newlayers[1].grid)   # iterate through all coordinates
        #if the probability of overlap is not NaN (ie. overlap) AND if prey occurence is not NaN value
        if !isnan(probs[k]) && !isnan(newlayers[j][k])
            #then multiply the probability of sp 1,2,3... occurence by the next sp
            probs[k] = probs[k] .* newlayers[j][k]
        end
    end
end
probs

#check if probability of occurence changed (thus overlap in species occurence) -- if newlayers[1] value changes or not
presence = probs .== newlayers[1].grid
#if no change, change to NaN
probs[presence] .= NaN
#filter only occurence data (ie. not NaN values)
filter(!isnan, probs)

# Convert to SDMLayer
map_layer = SDMLayer(probs, newlayers[1].left, newlayers[1].right, newlayers[1].bottom, newlayers[1].top)

plot_layer = plotSDM(map_layer)
plot_layer = plot!(title="Probability of species overlap: vulpes vulpes")
plot_layer

savefig("../SDM/figures/sdm_map-vulpes_probability_america.pdf")

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
## PARTE 2
cd(@__DIR__())
using DataStructures: Queue, enqueue!, nil, cons, list, Cons
maskrx = r"mask = (.*)"
memrx = r"mem\[([0-9]+)\] = ([0-9]+)"
function parse_mask_v2(address::T, s::AbstractString) where {T<:Unsigned}
    b1 = one(T)
    n = length(s)
    address′ = foldl(1:n; init=address) do a, i
        s[n-i+1] == '1' || return a
        a | b1 << (i - 1)
    end
    foldl(1:n; init=list(address′)) do addrs, i
        s[n-i+1] == 'X' || return addrs
        foldl(addrs; init=nil(T)) do xs, x
            mask = b1 << (i - 1)
            y = x | mask
            z = x & ~mask
            cons(z, cons(y, xs))
        end
    end
end
function parse_mem(::Type{T}, address, value) where {T<:Unsigned}
    parse(T, address) => parse(T, value)
end
function set_mem!(tape::Dict, mem, mask)
    isnothing(mem) && return
    address, value = mem
    addrs = parse_mask_v2(address, mask)
    foreach(addrs) do a
        tape[a] = value
    end
end
function run_program_v2(input, ::Type{T} = UInt64) where {T<:Unsigned}
    tape = Dict{T,Int}()
    mem = nothing
    mask = nothing
    foreach(eachline(input)) do line
        if (m = match(maskrx, line); !isnothing(m))
            mask = m[1]
        elseif (m = match(memrx, line); !isnothing(m))
            mem = parse_mem(T, m[1], m[2])
            set_mem!(tape, mem, mask)
        end
    end
    sum(values(tape))
end

##
@show run_program_v2("input.txt", UInt64);


text:
##################
#  Society Type  #
##################

abstract type Society{N, K, T} end

###########################
#  Society Redefinitions  #
###########################

length(soc::Society)       = length(soc.agents)
size(soc::Society)         = (length(soc.agents), length(soc.agents[1]))
size(soc::Society, i::Int) = size(soc)[i]

getindex(soc::Society, i::Int)                     = soc.agents[i]
getindex(soc::Society, c::Colon)                   = soc.agents[:]
getindex(soc::Society, c1::Colon, c2::Colon)       = soc.links[:, :]
endof(soc::Society)                                = soc.agents[end]

getindex(soc::Society, i::Int, j::Int)             = soc.links[i, j]
getindex(soc::Society, i::Int, r::UnitRange)       = soc.links[i, r]
getindex(soc::Society, r::UnitRange, j::Int)       = soc.links[r, j]
getindex(soc::Society, r::UnitRange, s::UnitRange) = soc.links[r, s]
getindex(soc::Society, i::Int, c::Colon)           = soc.links[i, c]
getindex(soc::Society, c::Colon, j::Int)           = soc.links[c, j]
getindex(soc::Society, c::Colon, r::UnitRange)     = soc.links[c, r]
getindex(soc::Society, r::UnitRange, c::Colon)     = soc.links[r, c]

setindex!(soc::Society, x::MoralVector, i::Int64) =
    soc.agents[i] = x

start(soc::Society)        = 1
done(soc::Society, s::Int) = s > length(soc)
next(soc::Society, s::Int) = (soc[s], s+1)

#########################
#  Society Definitions  #
#########################

"""
    agents(soc::Society)

Returns the agents vector of a Society soc
"""
agents(soc::Society)           = soc[:]

"""
    agents(soc::Society, i::Int)

Returns the ith agent of a Society soc
"""
agents(soc::Society, i::Int) = soc[i]

"""
    interactions(soc::Society)

Returns the links of a Society soc
"""
interactions(soc::Society)     = soc.links

"""
    insertagent!{N,K,T}(soc::Society{N,K,T}, x::MoralVector{K,T}, i::Int)

Insert new agent x at position i in Society soc
"""
insertagent!{N, K, T}(soc::Society{N,K,T}, x::MoralVector{K,T}, i::Int) =
    soc.agents[i] = x

"""
    hamiltonian{N,K,T}(soc::Society{N,K,T}, x::MoralVector{K,T})

Evaluates the hamiltonian, that is, the total social cost of a Society given some MoralVector x
"""
function hamiltonian{N,K,T}(soc::Society{N,K,T}, x::MoralVector{K,T})
    hamil = 0

    for j in 1:N for i in 1:N
        Rij = soc[i, j]

        if Rij != 0
            hamil += Rij*cogcost(soc,i,j,x)
        end
    end end

    return hamil
end

"""
    magnetization{N,K,T}(soc::Society{N,K,T}, x::MoralVector{K,T})

Evaluates the magnetization of a Society, that is, the total sum of the MoralVectors opinions
"""
function magnetization{N,K,T}(soc::Society{N,K,T}, x::MoralVector{K,T})
    mag = 0.

    for agi in agents(soc)
        mag += agi ⋅ x
    end

    return mag / N
end

"""
    believeness{N,K,T}(soc::Society{N,K,T}, x::MoralVector{K,T})

Evaluates the believeness of a Society, that is, the sum of the absolute values of MoralVectors opinions
"""
function believeness{N,K,T}(soc::Society{N,K,T}, x::MoralVector{K,T})
    bel = 0.

    for agi in agents(soc)
        bel += abs(agi ⋅ x)
    end

    return bel / N
end

"""
    quadrupole{N,K,T}(soc::Society{N,K,T}, x::MoralVector{K,T})

Evaluates the quadrupole of a Society, that is, something very nasty and obscure:

* \frac{1}{N^2 K} \sum_{\langle ij \rangle} R_{ij} j_i h_j
"""
function quadrupole{N,K,T}(soc::Society{N,K,T}, x::MoralVector{K,T})
    qd = 0.

    for i in 1:N for j in (i+1):N
        qd += soc[i, j] * (soc[i] ⋅ x) * (soc[j] ⋅ x)
    end end

    return qd / N^2
end

"""
    consensus{N,K,T}(soc::Society{N,K,T}, i::Int, j::Int)

Evaluates the consensus of MoralVectors i and j within a Society, that is, the between the agents inner representations
"""
function consensus{N,K,T}(soc::Society{N,K,T}, i::Int, j::Int)
    if 0 < i <= N && 0 < j <= N
        return soc[i] ⋅ soc[j]
    else
        error("The indices given are out of range")
    end
end

"""
    consensus{N,K,T}(soc::Society{N,K,T})

Evaluates the consensus of a Society, that is, the matrix of the correlations between the MoralVectors inner representations
"""
function consensus{N,K,T}(soc::Society{N,K,T}; format = :none)
    ψ = zeros(N, N)

    for i in 1:N for j in (i+1):N
        ψ[i, j] = consensus(soc, i, j)
    end end

    if format == :ordered
        return ψ[ψ .!= 0.]
    end
    return ψ
end


text:
# Minimally adapted from https://github.com/scipy/scipy/blob/b5d8bab88af61d61de09641243848df63380a67f/scipy/special/_spherical_bessel.pxd#L74
function sphericalbesselj(n::Integer, x::T) where {T<:AbstractFloat}
    isnan(x) && return x
    n < 0 && throw(DomainError(n, "n must be non-negative"))
    isinf(x) && return zero(T)
    iszero(x) && return iszero(n) ? one(T) : zero(T)

    if n > 0 && n ≥ x
        # We don't use SpecialFunctions.sphericalbesselj directly because it errors out on large x
        return T(SpecialFunctions.sphericalbesselj(n, x))
    else
        return _sphericalbesselj(n, x)
    end
end

function _sphericalbesselj(n::Integer, x::T) where {T<:AbstractFloat}
    invx = inv(x)
    s0 = sin(x) * invx
    iszero(n) && return s0

    s1 = (s0 - cos(x)) * invx
    isone(n) && return s1

    sn = zero(T)
    for idx in 2:n
        sn = (2idx - 1) * invx * s1 - s0
        s0 = s1
        s1 = sn
    end
    return sn
end

# Adapted from https://github.com/numpy/numpy/blob/4adc87dff15a247e417d50f10cc4def8e1c17a03/numpy/polynomial/legendre.py#L832-L914
function legval(x, c::AbstractVector)
    nd = length(c)
    nd ≥ 2 || return last(c)

    c0, c1 = @inbounds c[nd - 1], c[nd]
    @inbounds for j in (nd - 2):-1:1
        k = j / (j + 1)
        c0, c1 = c[j] - c1 * k, c0 + c1 * x * (k + 1)
    end
    return c0 + c1 * x
end

# Adapted from https://github.com/numpy/numpy/blob/4adc87dff15a247e417d50f10cc4def8e1c17a03/numpy/polynomial/legendre.py#L1126-L1176
"""
legvander(x, deg)

Pseudo-Vandermonde matrix of degree `deg`.
"""
function legvander(x::AbstractVector{T}, deg::Integer) where {T}
    deg ≥ 0 || throw(DomainError(deg, "degree needs to be non-negative"))

    vsize = (length(x), deg + 1)
    v = Matrix{T}(undef, vsize...)

    # Use forward recursion to generate the entries. This is not as accurate
    # as reverse recursion in this application but it is more efficient.
    @inbounds begin
        v[:, 1] .= one(T)
        if deg > 0
            v[:, 2] .= x
            for i in 2:deg
                invi = inv(i)
                @views @. v[:, i + 1] = v[:, i] * x * (2 - invi) - v[:, i - 1] * (1 - invi)
            end
        end
    end

    return v
end

# Adapted from https://github.com/numpy/numpy/blob/4adc87dff15a247e417d50f10cc4def8e1c17a03/numpy/polynomial/legendre.py#L612-L701
"""
    legder
"""
function legder(c::AbstractMatrix{T}) where {T}
    c = copy(c)
    n, m = size(c)
    n -= 1
    der = Matrix{T}(undef, n, m)
    @views @inbounds for j in n:-1:2
        @. der[j, :] = (2j - 1) * c[j + 1, :]
        @. c[j - 1, :] += c[j + 1, :]
    end
    @views @inbounds @. der[1, :] = c[2, :]
    return der
end


text:
#######
# Tag #
#######

struct Tag{F,H} end

# Here, we could've just as easily used `hash`; however, this
# is unsafe/undefined behavior if `hash(::Type{V})` is overloaded
# in a module loaded after ForwardDiff. Thus, we instead use
# `hash(Symbol(V))`, which is somewhat safer since it's far less
# likely that somebody would overwrite the Base definition for
# `Symbol(::DataType)` or `hash(::Symbol)`.
@generated function Tag(::Type{F}, ::Type{V}) where {F,V}
    H = hash(Symbol(V))
    return quote
        $(Expr(:meta, :inline))
        Tag{F,$H}()
    end
end

##################
# AbstractConfig #
##################

abstract type AbstractConfig{T,N} end

struct ConfigMismatchError{F,G,H} <: Exception
    f::F
    cfg::AbstractConfig{Tag{G,H}}
end

function Base.showerror(io::IO, e::ConfigMismatchError{F,G}) where {F,G}
    print(io, "The provided configuration (of type $(typeof(e.cfg))) was constructed for a",
              " function other than the current target function. ForwardDiff cannot safely",
              " perform differentiation in this context; see the following issue for details:",
              " https://github.com/JuliaDiff/ForwardDiff.jl/issues/83. You can resolve this",
              " problem by constructing and using a configuration with the appropriate target",
              " function, e.g. `ForwardDiff.GradientConfig($(e.f), x)`")
end

Base.copy(cfg::AbstractConfig) = deepcopy(cfg)

Base.eltype(cfg::AbstractConfig) = eltype(typeof(cfg))

@inline chunksize(::AbstractConfig{T,N}) where {T,N} = N

####################
# DerivativeConfig #
####################

struct DerivativeConfig{T,D} <: AbstractConfig{T,1}
    duals::D
end

"""
    ForwardDiff.DerivativeConfig(f!, y::AbstractArray, x::AbstractArray)

Return a `DerivativeConfig` instance based on the type of `f!`, and the types/shapes of the
output vector `y` and the input vector `x`.

The returned `DerivativeConfig` instance contains all the work buffers required by
`ForwardDiff.derivative` and `ForwardDiff.derivative!` when the target function takes the form
`f!(y, x)`.

If `f!` is `nothing` instead of the actual target function, then the returned instance can
be used with any target function. However, this will reduce ForwardDiff's ability to catch
and prevent perturbation confusion (see https://github.com/JuliaDiff/ForwardDiff.jl/issues/83).

This constructor does not store/modify `y` or `x`.
"""
function DerivativeConfig(::F,
                          y::AbstractArray{Y},
                          x::X,
                          ::T = Tag(F, X)) where {F,X<:Real,Y<:Real,T}
    duals = similar(y, Dual{T,Y,1})
    return DerivativeConfig{T,typeof(duals)}(duals)
end

Base.eltype(::Type{DerivativeConfig{T,D}}) where {T,D} = eltype(D)

##################
# GradientConfig #
##################

struct GradientConfig{T,V,N,D} <: AbstractConfig{T,N}
    seeds::NTuple{N,Partials{N,V}}
    duals::D
end

"""
    ForwardDiff.GradientConfig(f, x::AbstractArray, chunk::Chunk = Chunk(x))

Return a `GradientConfig` instance based on the type of `f` and type/shape of the input
vector `x`.

The returned `GradientConfig` instance contains all the work buffers required by
`ForwardDiff.gradient` and `ForwardDiff.gradient!`.

If `f` is `nothing` instead of the actual target function, then the returned instance can
be used with any target function. However, this will reduce ForwardDiff's ability to catch
and prevent perturbation confusion (see https://github.com/JuliaDiff/ForwardDiff.jl/issues/83).

This constructor does not store/modify `x`.
"""
function GradientConfig(::F,
                        x::AbstractArray{V},
                        ::Chunk{N} = Chunk(x),
                        ::T = Tag(F, V)) where {F,V,N,T}
    seeds = construct_seeds(Partials{N,V})
    duals = similar(x, Dual{T,V,N})
    return GradientConfig{T,V,N,typeof(duals)}(seeds, duals)
end

Base.eltype(::Type{GradientConfig{T,V,N,D}}) where {T,V,N,D} = Dual{T,V,N}

##################
# JacobianConfig #
##################

struct JacobianConfig{T,V,N,D} <: AbstractConfig{T,N}
    seeds::NTuple{N,Partials{N,V}}
    duals::D
end

"""
    ForwardDiff.JacobianConfig(f, x::AbstractArray, chunk::Chunk = Chunk(x))

Return a `JacobianConfig` instance based on the type of `f` and type/shape of the input
vector `x`.

The returned `JacobianConfig` instance contains all the work buffers required by
`ForwardDiff.jacobian` and `ForwardDiff.jacobian!` when the target function takes the form
`f(x)`.

If `f` is `nothing` instead of the actual target function, then the returned instance can
be used with any target function. However, this will reduce ForwardDiff's ability to catch
and prevent perturbation confusion (see https://github.com/JuliaDiff/ForwardDiff.jl/issues/83).

This constructor does not store/modify `x`.
"""
function JacobianConfig(::F,
                        x::AbstractArray{V},
                        ::Chunk{N} = Chunk(x),
                        ::T = Tag(F, V)) where {F,V,N,T}
    seeds = construct_seeds(Partials{N,V})
    duals = similar(x, Dual{T,V,N})
    return JacobianConfig{T,V,N,typeof(duals)}(seeds, duals)
end

"""
    ForwardDiff.JacobianConfig(f!, y::AbstractArray, x::AbstractArray, chunk::Chunk = Chunk(x))

Return a `JacobianConfig` instance based on the type of `f!`, and the types/shapes of the
output vector `y` and the input vector `x`.

The returned `JacobianConfig` instance contains all the work buffers required by
`ForwardDiff.jacobian` and `ForwardDiff.jacobian!` when the target function takes the form
`f!(y, x)`.

If `f!` is `nothing` instead of the actual target function, then the returned instance can
be used with any target function. However, this will reduce ForwardDiff's ability to catch
and prevent perturbation confusion (see https://github.com/JuliaDiff/ForwardDiff.jl/issues/83).

This constructor does not store/modify `y` or `x`.
"""
function JacobianConfig(::F,
                        y::AbstractArray{Y},
                        x::AbstractArray{X},
                        ::Chunk{N} = Chunk(x),
                        ::T = Tag(F, X)) where {F,Y,X,N,T}
    seeds = construct_seeds(Partials{N,X})
    yduals = similar(y, Dual{T,Y,N})
    xduals = similar(x, Dual{T,X,N})
    duals = (yduals, xduals)
    return JacobianConfig{T,X,N,typeof(duals)}(seeds, duals)
end

Base.eltype(::Type{JacobianConfig{T,V,N,D}}) where {T,V,N,D} = Dual{T,V,N}

#################
# HessianConfig #
#################

struct HessianConfig{T,V,N,D,H,DJ} <: AbstractConfig{T,N}
    jacobian_config::JacobianConfig{Tag{Void,H},V,N,DJ}
    gradient_config::GradientConfig{T,Dual{Tag{Void,H},V,N},N,D}
end

"""
    ForwardDiff.HessianConfig(f, x::AbstractArray, chunk::Chunk = Chunk(x))

Return a `HessianConfig` instance based on the type of `f` and type/shape of the input
vector `x`.

The returned `HessianConfig` instance contains all the work buffers required by
`ForwardDiff.hessian` and `ForwardDiff.hessian!`. For the latter, the buffers are
configured for the case where the `result` argument is an `AbstractArray`. If
it is a `DiffResult`, the `HessianConfig` should instead be constructed via
`ForwardDiff.HessianConfig(f, result, x, chunk)`.

If `f` is `nothing` instead of the actual target function, then the returned instance can
be used with any target function. However, this will reduce ForwardDiff's ability to catch
and prevent perturbation confusion (see https://github.com/JuliaDiff/ForwardDiff.jl/issues/83).

This constructor does not store/modify `x`.
"""
function HessianConfig(f::F,
                       x::AbstractArray{V},
                       chunk::Chunk = Chunk(x),
                       tag::Tag = Tag(F, Dual{Void,V,0})) where {F,V}
    jacobian_config = JacobianConfig(nothing, x, chunk)
    gradient_config = GradientConfig(f, jacobian_config.duals, chunk, tag)
    return HessianConfig(jacobian_config, gradient_config)
end

"""
    ForwardDiff.HessianConfig(f, result::DiffResult, x::AbstractArray, chunk::Chunk = Chunk(x))

Return a `HessianConfig` instance based on the type of `f`, types/storage in `result`, and
type/shape of the input vector `x`.

The returned `HessianConfig` instance contains all the work buffers required by
`ForwardDiff.hessian!` for the case where the `result` argument is an `DiffResult`.

If `f` is `nothing` instead of the actual target function, then the returned instance can
be used with any target function. However, this will reduce ForwardDiff's ability to catch
and prevent perturbation confusion (see https://github.com/JuliaDiff/ForwardDiff.jl/issues/83).

This constructor does not store/modify `x`.
"""
function HessianConfig(f::F,
                       result::DiffResult,
                       x::AbstractArray{V},
                       chunk::Chunk = Chunk(x),
                       tag::Tag = Tag(F, Dual{Void,V,0})) where {F,V}
    jacobian_config = JacobianConfig(nothing, DiffBase.gradient(result), x, chunk)
    gradient_config = GradientConfig(f, jacobian_config.duals[2], chunk, tag)
    return HessianConfig(jacobian_config, gradient_config)
end

Base.eltype(::Type{HessianConfig{T,V,N,D,H,DJ}}) where {T,V,N,D,H,DJ} = Dual{T,Dual{Tag{Void,H},V,N},N}


text:
module Day12

using AdventOfCode2019

function day12(input::String = readInput(joinpath(@__DIR__, "..", "data", "day12.txt")))
    stripChar = (s, r) -> replace(s, Regex("[$r]") => "")
    coordStringList = split(stripChar(input, "<>=xyz,"))
    length(coordStringList) % 3 == 0 || throw(AssertionError("Invalid input"))
    coords = parse.(Int, coordStringList)
    positions = Array{Array{Int,1},1}()
    for i = 1:3:length(coords)
        push!(positions, Array{Int,1}([coords[i], coords[i+1], coords[i+2]]))
    end

    velocities = [zeros(Int, 3) for i=1:length(positions)]

    initial_positions = deepcopy(positions)
    initial_velocities = deepcopy(velocities)

    # Part 1
    nsteps = 1000
    for i = 1:nsteps
        _timestep!(positions, velocities)
    end
    total_energy = _total_energy(positions, velocities)

    # Part 2
    positions = deepcopy(initial_positions)
    velocities = deepcopy(initial_velocities)

    cycles = zeros(Int, 3)
    n = 0
    while iszero(cycles[1]) || iszero(cycles[2]) || iszero(cycles[3])
        _timestep!(positions, velocities)
        n += 1
        for d in findall(x->x==0, cycles)
            if _equals(initial_positions, positions, d) && _equals(initial_velocities, velocities, d)
                cycles[d] = n
            end
        end
    end
    return [total_energy, lcm(cycles)]
end

function _equals(s1, s2, dim)
    length(s1) == length(s2) || return false
    for (x, y) in zip(s1, s2)
        x[dim] == y[dim] || return false
    end
    return true
end

function _update_velocities!(positions, velocities)
    for i = 1:length(positions)
        for j = i+1:length(positions)
            for k = 1:3
                if positions[i][k] > positions[j][k]
                    velocities[i][k] -= 1
                    velocities[j][k] += 1
                elseif positions[i][k] < positions[j][k]
                    velocities[i][k] += 1
                    velocities[j][k] -= 1
                end
            end
        end
    end
end

function _update_positions!(positions, velocities)
    positions .+= velocities
end

function _timestep!(positions, velocities)
    _update_velocities!(positions, velocities)
    _update_positions!(positions, velocities)
end

function _potential_energies(positions)
    return [abs.(pos) |> sum for pos in positions]
end

function _kinetic_energies(velocities)
    return [abs.(vel) |> sum for vel in velocities]
end

function _total_energies(positions, velocities)
    return _potential_energies(positions) .* _kinetic_energies(velocities)
end

function _total_energy(positions, velocities)
    return _total_energies(positions, velocities) |> sum
end

end # module


text:
using Plots

# Model
include_model("walker")

model = free_time_model(model)

# Visualize
# - Pkg.add any external deps from visualize.jl
include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)

"""
Useful methods

joint-space inertia matrix:
	M_func(model, q)

position (px, pz) of foot 1 heel:
	kinematics_2(model,
		get_q⁺(x), body = :calf_1, mode = :ee)

position (px, pz) of foot 2 heel:
	kinematics_2(model,
		get_q⁺(x), body = :calf_2, mode = :ee)

position (px, pz) of foot 1 toe:
		kinematics_3(model,
	        get_q⁺(x), body = :foot_1, mode = :ee)

position (px, pz) of foot 2 toe:
		kinematics_3(model,
	        get_q⁺(x), body = :foot_2, mode = :ee)

jacobian to foot 1 heel
	jacobian_2(model, get_q⁺(x), body = :calf_1, mode = :ee)

jacobian to foot 2 heel
	jacobian_2(model, get_q⁺(x), body = :calf_2, mode = :ee)

jacobian to foot 1 toe
	jacobian_3(model, get_q⁺(x), body = :foot_1, mode = :ee)

jacobian to foot 2 toe
	jacobian_3(model, get_q⁺(x), body = :foot_2, mode = :ee)

gravity compensating torques
	initial_torque(model, q, h)[model.idx_u] # returns torques for system
	initial_torque(model, q, h) # returns torques for system and contact forces

	-NOTE: this method is iterative and is not differentiable (yet)
		-TODO: implicit function theorem to compute derivatives

contact models
	- maximum-dissipation principle:
		contact_constraints(model, T)
	- no slip:
		contact_no_slip_constraints(T)

loop constraint
	- constrain variables at two indices to be the same:
		loop_constraints(model, x_idx, t1_idx, t2_idx)
"""

function get_q⁺(x)
	view(x, model.nq .+ (1:model.nq))
end

function ellipse_traj(x_start, x_goal, z, T)
	dist = x_goal - x_start
	a = 0.5 * dist
	b = z

	z̄ = 0.0

	x = range(x_start, stop = x_goal, length = T)

	z = sqrt.(max.(0.0, (b^2) * (1.0 .- ((x .- (x_start + a)).^2.0) / (a^2.0))))

	return x, z
end

# Horizon
T = 31
T_fix = 7

# Time step
tf = 0.5
h = tf / (T - 1)

# Permutation matrix
perm = @SMatrix [1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
                 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
				 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0;
				 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0;
				 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0;
				 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0;
				 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0;
				 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0;
				 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0]

# Configurations
# 1: x pos
# 2: z pos
# 3: torso angle (rel. to downward vertical)
# 4: thigh 1 angle (rel. to downward vertical)
# 5: calf 1 (rel. to downward vertical)
# 6: thigh 2 (rel. to downward vertical)
# 7: calf 2 (rel. to downward vertical)
# 8: foot 1 (rel. to downward vertical)
# 9: foot 2 (rel. to downward vertical)
function initial_configuration(model, θ_torso, θ_thigh_1, θ_leg_1, θ_thigh_2)
    q1 = zeros(model.nq)
    q1[3] = θ_torso
    q1[4] = θ_thigh_1
    q1[5] = θ_leg_1
    z1 = model.l_thigh1 * cos(q1[4]) + model.l_calf1 * cos(q1[5])
    q1[6] = θ_thigh_2
    q1[7] = -1.0 * acos((z1 - model.l_thigh2 * cos(q1[6])) / model.l_calf2)
	# q1[6] = θ_thigh_1
	# q1[7] = θ_leg_1
	q1[2] = z1
	q1[8] = pi / 2.0
	q1[9] = pi / 2.0

    # p1 = kinematics_2(model, q1, body = :calf_1, mode = :ee)
    # p2 = kinematics_2(model, q1, body = :calf_2, mode = :ee)
    # @show stride = abs(p1[1] - p2[1])
	#
    # q1[1] = -0.75#-1.0 * p1[1]
	#
    # qM = copy(q1)
    # qM[4] = q1[6]
    # qM[5] = q1[7]
    # qM[6] = q1[4]
    # qM[7] = q1[5]
    # qM[1] = abs(p2[1])
	#
    # pM_1 = kinematics_2(model, qM, body = :calf_1, mode = :ee)
    # pM_2 = kinematics_2(model, qM, body = :calf_2, mode = :ee)
	#
    # qT = copy(q1)
    # qT[1] = 0.75#q1[1] + 1.0#2 * stride
	#
    # pT_1 = kinematics_2(model, qT, body = :calf_1, mode = :ee)
    # pT_2 = kinematics_2(model, qT, body = :calf_2, mode = :ee)

    return q1
end

q1 = initial_configuration(model, -pi / 100.0, pi / 20.0, -pi / 20.0, -pi / 40.0)
pf1 = kinematics_3(model, q1, body = :foot_1, mode = :com)
pf2 = kinematics_3(model, q1, body = :foot_2, mode = :com)

strd = 2 * (pf1 - pf2)[1]

qT = Array(perm) * copy(q1)
qT[1] += 0.5 * strd
q_ref = linear_interpolation(q1, qT, T+1)
x_ref = configuration_to_state(q_ref)
visualize!(vis, model, q_ref, Δt = h)

zh = 0.075
xf1_el, zf1_el = ellipse_traj(pf1[1], pf1[1] + strd, zh, T - T_fix)
xf2_el, zf2_el = ellipse_traj(pf2[1], pf2[1] + strd, zh, T - T_fix)

zf2 = [[zf2_el[1] for t = 1:T_fix]..., zf2_el...]
xf2 = [[xf2_el[1] for t = 1:T_fix]..., xf2_el...]
zf1 = [zf1_el[1] for t = 1:T]
xf1 = [xf1_el[1] for t = 1:T]

p1_ref = [[xf1[t]; zf1[t]] for t = 1:T]
p2_ref = [[xf2[t]; zf2[t]] for t = 1:T]

plot(hcat(p1_ref...)', legend = :topleft)
plot!(hcat(p2_ref...)')

# Control
# u = (τ1..7, λ1..4, β1..8, ψ1..4, η1...8, s1)
# τ1: torso angle
# τ2: thigh 1 angle
# τ3: calf 1
# τ4: thigh 2
# τ5: calf 2
# τ6: foot 1
# τ7: foot 2

# ul <= u <= uu
# u1 = initial_torque(model, q1, h)[model.idx_u] # gravity compensation for current q
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= Inf
_uu[end] = 2.0 * h
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf
_ul[end] = 0.25 * h
ul, uu = control_bounds(model, T, _ul, _uu)

qL = [-Inf; -Inf; q1[3] - pi / 50.0; q1[4:end] .- pi / 6.0; -Inf; -Inf; q1[3] - pi / 50.0; q1[4:end] .- pi / 6.0]
qU = [Inf; q1[2] + 0.001; 0.0; q1[4:end] .+ pi / 6.0; Inf; Inf; q1[3:end] .+ pi / 6.0]
qL[8] = q1[8] - pi / 10.0
qL[9] = q1[9] - pi / 10.0
qL[9 + 8] = q1[8] - pi / 10.0
qL[9 + 9] = q1[9] - pi / 10.0

qU[8] = q1[8] + pi / 10.0
qU[9] = q1[9] + pi / 10.0
qU[9 + 8] = q1[8] + pi / 10.0
qU[9 + 9] = q1[9] + pi / 10.0

xl, xu = state_bounds(model, T,
    qL, qU,
    x1 = [Inf * ones(model.nq); q1],
	xT = [Inf * ones(model.nq); qT[1]; Inf * ones(model.nq - 1)])

# Objective
include_objective(["velocity", "nonlinear_stage", "control_velocity"])

x0 = configuration_to_state(q_ref)

# penalty on slack variable
obj_penalty = PenaltyObjective(1.0e5, model.m-1)

# quadratic tracking objective
# Σ (x - xref)' Q (x - x_ref) + (u - u_ref)' R (u - u_ref)
q_penalty = 1.0e-1 * ones(model.nq)

x_penalty = [q_penalty; q_penalty]
obj_control = quadratic_time_tracking_objective(
    [Diagonal(x_penalty) for t = 1:T],
    [Diagonal([1.0e-4 * ones(model.nu)..., 1.0 * ones(model.nc)..., 1.0 * ones(model.nb)..., 1.0e-8 * ones(model.m - model.nu - model.nc - model.nb - 1)..., 0.0]) for t = 1:T-1],
    [x_ref[end] for t = 1:T],
    [[zeros(model.nu); zeros(model.m - model.nu)] for t = 1:T-1],
	1.0)

obj_ctrl_vel = control_velocity_objective(Diagonal([1.0e-1 * ones(model.nu); 1.0 * ones(model.nc + model.nb); zeros(model.m - model.nu - model.nc - model.nb)]))

# quadratic velocity penalty
q_v = 1.0e-2 * ones(model.nq)
obj_velocity = velocity_objective(
    [Diagonal(q_v) for t = 1:T-1],
    model.nq,
    h = h)

function l_foot_height(x, u, t)
	J = 0.0
	q1 = view(x, 1:9)
	# q2 = view(x, 9 .+ (1:9))

	if true
		J += 10000.0 * sum((p2_ref[t] - kinematics_3(model, q1, body = :foot_2, mode = :toe)).^2.0)
		J += 10000.0 * sum((p2_ref[t] - kinematics_3(model, q1, body = :foot_2, mode = :heel)).^2.0)
		J += 10000.0 * (q1[3] - qT[3])^2.0

	end

	return J
end

l_foot_height(x) = l_foot_height(x, nothing, T)

obj_foot_height = nonlinear_stage_objective(l_foot_height, l_foot_height)

obj = MultiObjective([obj_penalty,
                      obj_control,
                      obj_velocity,
					  obj_ctrl_vel,
					  obj_foot_height])
# Constraints
include_constraints(["contact", "loop", "free_time", "stage"])

function pinned1!(c, x, u, t)
    q1 = view(x, 1:9)
    c[1:2] = p1_ref[t] - kinematics_3(model, q1, body = :foot_1, mode = :com)
	nothing
end

function pinned2!(c, x, u, t)
    q1 = view(x, 1:9)
	c[1:2] = p2_ref[t] - kinematics_3(model, q1, body = :foot_2, mode = :com)
	nothing
end

n_stage = 2
t_idx1 = vcat([1:T]...)
t_idx2 = vcat([t for t = 1:T_fix]..., T)

con_pinned1 = stage_constraints(pinned1!, n_stage, (1:0), t_idx1)
con_pinned2 = stage_constraints(pinned2!, n_stage, (1:0), t_idx2)

con_loop = loop_constraints(model, collect([(2:9)...,(11:18)...]), 1, T, perm = Array(cat(perm, perm, dims = (1:2))))
con_contact = contact_constraints(model, T)
con_free_time = free_time_constraints(T)
con = multiple_constraints([con_contact,
	con_free_time,
	con_loop,
	con_pinned1,
	con_pinned2])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# trajectory initialization
u0 = [[1.0e-2 * randn(model.nu); 0.01 * randn(model.m - model.nu - 1); h] for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

# Solve
# NOTE: run multiple times to get good trajectory
@time z̄, info = solve(prob, copy(z0),
    nlp = :ipopt,
	max_iter = 2000,
    tol = 1.0e-3, c_tol = 1.0e-3, mapl = 5,
    time_limit = 60 * 3)

@show check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
_tf, _t, h̄ = get_time(ū)

q = state_to_configuration(x̄)
u = [u[model.idx_u] for u in ū]
γ = [u[model.idx_λ] for u in ū]
b = [u[model.idx_b] for u in ū]
ψ = [u[model.idx_ψ] for u in ū]
η = [u[model.idx_η] for u in ū]
h̄ = mean(h̄)
#
# norm(hcat([γ[t] .* ϕ_func(model, q[t+2]) for t = 1:T-1]...),Inf)
# norm(hcat([b[t] .* η[t] for t = 1:T-1]...), Inf)
# function friction_cone(model::Flamingo, λ, b)
# 	# @show model.μ
# 	return @SVector [model.μ * λ[1] - sum(b[1:2]),
# 					 model.μ * λ[2] - sum(b[3:4]),
# 					 model.μ * λ[3] - sum(b[5:6]),
# 					 model.μ * λ[4] - sum(b[7:8])]
# end
# norm(hcat([ψ[t] .* friction_cone(model, γ[t], b[t]) for t = 1:T-1]...), Inf)
#
# maximum([norm(fd(model, x̄[t+1], x̄[t], ū[t], zeros(model.d), h̄, t)) for t = 1:T-1])
# (q[end][1] - q[1][1]) / _tf
#
# vis = Visualizer()
# render(vis)
# visualize!(vis, model,
# 	[[x̄[1][1:model.nq] for i = 1:10]...,
# 	 state_to_configuration(x̄)...,
# 	 [x̄[end][model.nq .+ (1:model.nq)] for i = 1:10]...], Δt = h̄[1])
#
#
# plot(hcat(u...)', linetype = :steppost)
# plot(hcat([γt[1:2] for γt in γ]...)', linetype = :steppost)
# plot!(hcat([γt[3:4] for γt in γ]...)', linetype = :steppost)
#
# plot!(hcat([bt[1:2] for bt in b]...)', linetype = :steppost)
# plot!(hcat([bt[3:4] for bt in b]...)', linetype = :steppost)
# plot!(hcat([bt[5:6] for bt in b]...)', linetype = :steppost)
# plot!(hcat([bt[7:8] for bt in b]...)', linetype = :steppost)
#
_pf1 = [kinematics_3(model, qt, body = :foot_1, mode = :com) for qt in q]
_pf2 = [kinematics_3(model, qt, body = :foot_2, mode = :com) for qt in q]

plot(hcat(p1_ref...)', title = "foot 1",
 	legend = :topleft, label = ["x" "z"], color = :black, width = 2.0)
plot!(hcat(_pf1...)', legend = :topleft, color = :red, width = 1.0)

plot(hcat(p2_ref...)', title = "foot 2",
	legend = :topleft, label = ["x" "z"], color = :black, width = 2.0)
plot!(hcat(_pf2...)', title = "foot 2", legend = :bottomright,
	color = :red, width = 1.0)
#
# # function get_q_viz(q̄; N = 4)
# # 	q_viz = [q̄...]
# # 	shift_vec = zeros(model.nq)
# # 	shift_vec[1] = q̄[end][1]
# # 	for i = 1:N
# # 		q_update = [q + shift_vec for q in q̄[2:end]]
# # 		push!(q_viz, q_update...)
# # 		shift_vec[1] = q_update[end][1]
# # 	end
# #
# # 	return q_viz
# # end
#
# # q_viz = get_q_viz(q)
#
#
perm4 = [0.0 0.0 1.0 0.0;
         0.0 0.0 0.0 1.0;
		 1.0 0.0 0.0 0.0;
		 0.0 1.0 0.0 0.0]

perm7 = [0.0 0.0 1.0 0.0 0.0 0.0;
		 0.0 0.0 0.0 1.0 0.0 0.0;
		 1.0 0.0 0.0 0.0 0.0 0.0;
		 0.0 1.0 0.0 0.0 0.0 0.0;
		 0.0 0.0 0.0 0.0 0.0 1.0;
		 0.0 0.0 0.0 0.0 1.0 0.0]

perm8 = [0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0;
         0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0;
		 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0;
		 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0;
		 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
		 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0;
		 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0;
		 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0]

function mirror_gait(q, u, γ, b, ψ, η, T)
	qm = [deepcopy(q)...]
	um = [deepcopy(u)...]
	γm = [deepcopy(γ)...]
	bm = [deepcopy(b)...]
	ψm = [deepcopy(ψ)...]
	ηm = [deepcopy(η)...]

	stride = zero(qm[1])
	@show stride[1] = q[T+1][1] - q[2][1]
	@show 0.5 * strd
	for t = 1:T-1
		push!(qm, Array(perm) * q[t+2] + stride)
		push!(um, perm7 * u[t])
		push!(γm, perm4 * γ[t])
		push!(bm, perm8 * b[t])
		push!(ψm, perm4 * ψ[t])
		push!(ηm, perm8 * η[t])
	end

	return qm, um, γm, bm, ψm, ηm
end

hm = h̄
μm = model.μ
qm, um, γm, bm, ψm, ηm = mirror_gait(q, u, γ, b, ψ, η, T)

@save joinpath(@__DIR__, "biped_mirror_gait.jld2") qm um γm bm ψm ηm μm hm
#
# plot(hcat(q...)', color = :black, width = 2.0, label = "")
# plot!(hcat(qm...)', color = :red, width = 1.0, label = "")
#
# plot(hcat(u...)', color = :black, width = 2.0, label = "", linetype = :steppost)
# plot!(hcat(um...)', color = :red, width = 1.0, label = "", linetype = :steppost)
#
# plot(hcat(γ...)', color = :black, width = 2.0, label = "", linetype = :steppost)
# plot!(hcat(γm...)', color = :red, width = 1.0, label = "", linetype = :steppost)
#
# plot(hcat(b...)', color = :black, width = 2.0, label = "", linetype = :steppost)
# plot!(hcat(bm...)', color = :red, width = 1.0, label = "", linetype = :steppost)
#
vis = Visualizer()
render(vis)
visualize!(vis, model,
	qm,
	Δt = h̄[1])

#
# xm = configuration_to_state(qm)
# Tm = length(um) + 1
# norm(hcat([γm[t] .* ϕ_func(model, qm[t+2]) for t = 1:T-1]...),Inf)
# norm(hcat([bm[t] .* ηm[t] for t = 1:Tm-1]...), Inf)
# function friction_cone(model::Walker, λ, b)
# 	# @show model.μ
# 	return @SVector [model.μ * λ[1] - sum(b[1:2]),
# 					 model.μ * λ[2] - sum(b[3:4]),
# 					 model.μ * λ[3] - sum(b[5:6]),
# 					 model.μ * λ[4] - sum(b[7:8])]
# end
# norm(hcat([ψm[t] .* friction_cone(model, γm[t], bm[t]) for t = 1:Tm-1]...), Inf)
#
# maximum([norm(fd(model, xm[t+1], xm[t], [um[t]; γm[t]; bm[t]; ψm[t]; ηm[t]; 0.0; hm], zeros(model.d), hm, t)) for t = 1:Tm-1])

max_issues_repo_name: thowell/DirectMotionPlanning


text:
"""
  perp(u,n)

  cross product of manifold's outward unit normal and vector-valued field u
"""
function perp(u,n)
   n×u
end
const ⟂ = perp

"""
Outward unit normal of sphere
"""
function normal_unit_sphere(θϕ)
 θ,ϕ = θϕ
 n=spherical_to_cartesian_matrix(VectorValue(θ,ϕ,1.0))⋅VectorValue(0,0,1)
 n
end

"""
Jacobian of parametric unit sphere
"""
function J_unit_sphere(θϕ)
  θ,ϕ = θϕ
  TensorValue{3,2}(-sin(θ)*cos(ϕ), cos(θ)*cos(ϕ),      0,
                   -sin(ϕ)*cos(θ),-sin(ϕ)*sin(θ), cos(ϕ))
end

"""
First fundamental form parametric unit sphere (a.k.a. metric / Grammian)
"""
function G_unit_sphere(θϕ)
  Jθϕ=J_unit_sphere(θϕ)
  transpose(Jθϕ)⋅Jθϕ
end

"""
Spherical divergence on unit sphere
IMPORTANT NOTE:
* The input function "v" and output function
  take values on the parametric space
"""
function divergence_unit_sphere(v)
  function tmp(θϕ)
     Gθϕ=G_unit_sphere(θϕ)
     Jθϕ=J_unit_sphere(θϕ)
     function f(θϕ)
        sqrt(det(Gθϕ))*inv(Gθϕ)⋅transpose(Jθϕ)⋅(v(θϕ))
     end
     1.0/sqrt(det(Gθϕ))*(∇⋅(f))(θϕ)
  end
end

"""
Spherical Laplacian on unit sphere
IMPORTANT NOTE:
  * The input function "v" and output function
    take values on the parametric space
"""
function laplacian_unit_sphere(v)
  # function tmp(θϕ)
  #    Gθϕ=G_unit_sphere(θϕ)
  #    Jθϕ=J_unit_sphere(θϕ)
  #    function f(θϕ)
  #       sqrt(det(Gθϕ))*inv(Gθϕ)⋅(∇(v)(θϕ))
  #    end
  #    1.0/sqrt(det(Gθϕ))*(∇⋅(f))(θϕ)
  # end
  function tmp(θϕ)
    divergence_unit_sphere(gradient_unit_sphere(v))(θϕ)
  end
end

"""
Spherical gradient on unit sphere
IMPORTANT NOTE:
* The input function "v" and output function
  take values on the parametric space
"""
function gradient_unit_sphere(v)
  function tmp(θϕ)
     Gθϕ=G_unit_sphere(θϕ)
     Jθϕ=J_unit_sphere(θϕ)
     Jθϕ⋅inv(Gθϕ)⋅(∇(v))(θϕ)
  end
end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
using Plasm, ViewerGL, LinearAlgebra
GL = ViewerGL
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

#include("")
store = [];
scaling = 1.5;
V,(VV,EV,FV,CV) = Lar.cuboid([0.25,0.25,0.25],true,[-0.25,-0.25,-0.25]);
mybox = (V,CV,FV,EV);
for k=1:5
	size = rand()*scaling
	scale = Lar.s(size,size,size)
	transl = Lar.t(rand(3)...)
	alpha = 2*pi*rand()
	rx = Lar.r(alpha,0,0); ry = Lar.r(0,alpha,0); rz = Lar.r(0,0,alpha)
	rot = rx * ry * rz
	str = Lar.Struct([ transl, scale, rot, mybox ])
	obj = Lar.struct2lar(str)
	vs = obj[1]
	diag = LinearAlgebra.norm(vs[:,8]-vs[:,1])
	if diag > 1/5
		push!(store, obj)
	end
end

str = Lar.Struct(store);
V,CV,FV,EV = Lar.struct2lar(str);
# V = Plasm.normalize3D(V) TODO:  solve MethodError bug

open("/tmp/lar.txt", "w") do f
	write(f, "V = $V\n\n")
	write(f, "CV = $CV\n\n")
	write(f, "FV = $FV\n\n")
	write(f, "EV = $EV\n\n")
	close(f)
end

GL.VIEW([ GL.GLPol(V,CV, GL.COLORS[2], 0.1) ]);

function testarrangement(V,CV,FV,EV)
		cop_EV = Lar.coboundary_0(EV::Lar.Cells);
		cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
		W = convert(Lar.Points, V');

		V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement(
				W::Lar.Points, cop_EV::Lar.ChainOp, cop_FE::Lar.ChainOp);

		V = convert(Lar.Points, V');
		V,CVs,FVs,EVs = Lar.pols2tria(V, copEV, copFE, copCF) # whole assembly
		GL.VIEW(GL.GLExplode(V,FVs,1.1,1.1,1.1,99,1));
		GL.VIEW(GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1));
		GL.VIEW(GL.GLExplode(V,CVs,1,1,1,99,0.2));
end

testarrangement(V,CV,FV,EV);


text:
using Test
using Boltzmann
using BenchmarkTools

"""
Runs a series of basic smoke and benchmark tests on the interface of any RBM.

Default methods tested:

* `fit(model, X)`
* `transform(model, X)`
* `generate(model, X)`

Args:
* model::Union{AbstractRBM, Net} - the model to test

Optional:
* opts::Dict - the context passed when running fit. (default=DEFAULT_CONTEXT)
* n_obs::Int - total number of observations in generated dataset (default=1000)

Returns: BenchmarkGroup
"""
function benchmark!(rbm::AbstractRBM{T,V,H}, suite::BenchmarkGroup; opts::Dict=DefaultContext(), input_size=-1, debug=false, n_obs=(1000, 5000, 10000)) where {T,V,H}
    tname = String(typeof(rbm).name)
    suite[tname, String(T.name)] = BenchmarkGroup()

    if debug && isa(ctx[:reporter], TestReporter)
        ctx[:reporter].log = true
    end

    for m in n_obs
        ctx = deepcopy(opts)
        n_hid, n_vis = size(rbm.W)

        if input_size > 0
            n_vis = input_size
        end

        DX = generate_dataset(T, n_vis; n_obs=n_obs)
        SX = generate_dataset(T, n_vis; n_obs=n_obs, sparsity=0.3)

        for X in (DX, SX)
            suite[tname]["fit", size(X), String(typeof(X).name)] = @benchmarkable fit($rbm, $X, $ctx)
            suite[tname]["transform", size(X), String(typeof(X).name)] = @benchmarkable transform($rbm, $X)
            suite[tname]["generate", size(X), String(typeof(X).name)] = @benchmarkable generate($rbm, $X)
        end
    end

    return suite
end

#=
Available test groups are:
* :unit,
* :integration,
* :acceptance,
* :benchmark,
=#
const DEFAULT_TEST_GROUPS = Set([:unit, :integration])
TEST_GROUPS = DEFAULT_TEST_GROUPS

# If the JULIA_TEST_GROUPS ENV variable is set
# then use the inte
if haskey(ENV, "JULIA_TEST_GROUPS")
    TEST_GROUPS = Set(
        Symbol.(split(ENV["JULIA_TEST_GROUPS"]))
    )
end

tests = [
    "utils",
    "rbm",
    "nets",
    "conditional",
]

for t in tests
    println("\n\nRunning $t tests...\n")
    include(string("test", t, ".jl"))
end


text:
module ElasticFusion

import Images

include("impl.jl")

export Fusion, process_frame!, render, save_ply

using ElasticFusionLib_jll

end # module