text:
"""

    TransData{FT<:AbstractFloat}

Store information relevant to the particular geometrical configuration, not intrinsic to
molecules themselves, involving rigid geometrical transformations (translations and
rotations); this depends on the number of molecules `nummols` from
`MolSystem{T<:OneMol{<:AbstractFloat}}`.

`numTrans` is the number of transformations ``t`` (an integer).

`transArray` is a `3*nummols`-by-`numTrans` real array where each 3-element vector
`transArray[3*nn-2:3*nn, tt]` is the center-of-mass translation vector (relative to the
baseline) for molecule `nn` in transformation labeled `tt`.

`rotArray` is a `3*nummols`-by-`3*numTrans` real array where each 3-by-3 block
`rotArray[3*nn-2:3*nn, 3*tt-2:3*tt]` is the center-of-mass rotation matrix (relative to the
baseline) for molecule `nn` in transformation labeled `tt`.

`transLabels` is a `numTrans`-element vector of strings labeling each transformation `tt`.

`changedFromBefore` is an `nummols`-by-`numTrans` `BitArray` such that
`changedFromBefore[nn, tt]` is `true` if molecule `nn` is transformed in transformation `tt`
relative to transformation `tt - 1`, `false` otherwise.

`allAtomPos` is the horizontal concatenation of atomic coordinates (initially unmodified,
later modified, by transformations) across all molecules.

`baselineTrans` is a `3*nummols`-element real vector where each 3-element vector
`baselineTrans[3*nn-2:3*nn]` is the baseline center of mass coordinate for molecule `nn`.

baselineRot is a `3*nummols`-by-3 real matrix where each 3-by-3 block
`baselineRot[3*nn-2:3*nn,:]` is the baseline center of mass rotation matrix for the
orientation of molecule `nn` (using the alibi convention).
"""
struct TransData{FT<:AbstractFloat}

    transArray::Array{FT, 2}
    rotArray::Array{FT, 2}
    transLabels::Array{<:AbstractString, 1}
    changedFromBefore::BitArray{2}
    numTrans::Integer
    allAtomPos::Array{FT, 2}
    baselineTrans::Array{FT, 1}
    baselineRot::Array{FT, 2}

end

"""

    readTransData(myMolData::MolSystem{<:OneMol{FT}},
                  translistfilename::AbstractString="") where FT<:AbstractFloat

Construct TransData from existing molecular data `myMolData` and a file of name
`translistfilename` containing information about geometric transformations; if no file is
specified, produce a default set of transformations (no displacements or rotations at all).

See also: [`MolSystem{<:OneMol{FT}}`](@ref)

"""
function readTransData(myMolData::MolSystem{<:OneMol{FT}},
                       translistfilename::AbstractString="") where FT<:AbstractFloat

    ## default: no transformation file specified
    if (translistfilename == "")
        numtrans = 1;
        allmoltrans = zeros(FT, numtrans, 6*myMolData.nummols);
        allmoltranslabels = ["DEFAULT"];
        changedfrombefore = falses(myMolData.nummols, numtrans);
        transArray = zeros(FT, 3*myMolData.nummols, numtrans);
        rotArray = repeat(Matrix{FT}(I, 3, 3), myMolData.nummols, numtrans);

        allAtomPos = Array{FT,2}(undef, 3, myMolData.cumnumatomslist[end]);
        concatenateAtomPos!(myMolData, allAtomPos);

        baselineTrans = zeros(FT, 3*myMolData.nummols);
        baselineRot = repeat(Matrix{FT}(I, 3, 3), myMolData.nummols, 1);

        myTransData = TransData{FT}(transArray, rotArray, allmoltranslabels,
                                    changedfrombefore, numtrans, allAtomPos,
                                    baselineTrans, baselineRot);
        initializeAtomPos!(myTransData, myMolData);
        return myTransData;
    end

    # read file and format each line to get rid of extraneous spaces and things like that
    transstrarray = strip.(chomp.(readlines(translistfilename)));
    transstartinds = findall(map(u -> startswith(u, "TRANS "), transstrarray));
    transendinds = findall(transstrarray .== "ENDTRANS");
    numtrans = length(transendinds);
    if (length(transstartinds) != numtrans)
        error(string(translistfilename, ": number of TRANS and ENDTRANS statements",
                     " should be the same"));
    end
    transbothinds = vec(transpose(hcat(transstartinds, transendinds)));
    if (transbothinds != sort(transbothinds))
        error(string(translistfilename, ": TRANS and ENDTRANS statements should",
                     " alternate in ordered occurrence"));
    end

    # populate the lists of transformation labels, translation vectors, and Euler angles
    allmoltranslabels = Array{String, 1}(undef, numtrans);
    allmoltrans = zeros(FT, 6*myMolData.nummols + 1, numtrans);
    for tt=1:numtrans
        allmoltranslabels[tt] = replace(strip(chomp(transstrarray[transstartinds[tt]])),
                                        "TRANS " => "");
        currtransstrarray = strip.(chomp.(transstrarray[transstartinds[tt]+1:transendinds[tt]-1]));
        for currtransstr in currtransstrarray
            currtransstrsplit = split(currtransstr);
            currmollabel = currtransstrsplit[1];
            currmolind = findfirst(myMolData.molLabels .== currmollabel);

            if (currmolind != nothing)
                ## assemble transformation list in the following order for each row in a transformation:
                ## [phi_1, theta_1, psi_1, ..., phi_N, theta_N, psi_N, x_1, y_1, z_1, ..., x_N, y_N, z_N]
                ## Note: all Euler angles should be given in the transformation file as ALIBI transformations in DEGREES
                allmoltrans[3*myMolData.nummols+3*currmolind-2,
                            tt] = parse(FT,currtransstrsplit[2]);
                allmoltrans[3*myMolData.nummols+3*currmolind-1,
                            tt] = parse(FT, currtransstrsplit[3]);
                allmoltrans[3*myMolData.nummols+3*currmolind,
                            tt] = parse(FT, currtransstrsplit[4]);
                if (length(currtransstrsplit) == 7)
                    allmoltrans[3*currmolind-2, tt] = parse(FT, currtransstrsplit[5]);
                    allmoltrans[3*currmolind-1, tt] = parse(FT, currtransstrsplit[6]);
                    allmoltrans[3*currmolind, tt] = parse(FT, currtransstrsplit[7]);
                end
            end
        end
        allmoltrans[end, tt] = tt;
    end

    # populate the lists of baseline translation vectors and Euler angles
    baselinestartind = findfirst(map(u -> startswith(u, "BASELINE"), transstrarray));
    baselineendind = findfirst(map(u -> startswith(u, "ENDBASELINE"), transstrarray));
    baselinetrans = zeros(FT, 6*myMolData.nummols);

    if (baselinestartind != nothing && baselineendind != nothing)
        baselinestrarray = strip.(chomp.(transstrarray[(baselinestartind+
                                                        1):(baselineendind-1)]));
        if (length(baselinestrarray) > 0)
            for currbaselinestr in baselinestrarray
                currbaselinestrsplit = split(currbaselinestr);
                currmollabel = currbaselinestrsplit[1];
                currmolind = findfirst(myMolData.molLabels .== currmollabel);

                if (currmolind != nothing)
                    ## Note: all Euler angles should be given in the transformation file in DEGREES
                    baselinetrans[3*myMolData.nummols+3*currmolind-2] = parse(FT, currbaselinestrsplit[2]);
                    baselinetrans[3*myMolData.nummols+3*currmolind-1] = parse(FT, currbaselinestrsplit[3]);
                    baselinetrans[3*myMolData.nummols+3*currmolind] = parse(FT, currbaselinestrsplit[4]);
                    if (length(currbaselinestrsplit) == 7)
                        baselinetrans[3*currmolind-2] = parse(FT, currbaselinestrsplit[5]);
                        baselinetrans[3*currmolind-1] = parse(FT, currbaselinestrsplit[6]);
                        baselinetrans[3*currmolind] = parse(FT, currbaselinestrsplit[7]);
                    end
                end
            end
        end
    end

    # sort such that the displacement of the last molecule changes first
    # between transformation labels, then the displacement of the
    # next-to-last molecule, and so on, then the rotation of the last
    # molecule, and so on
    allmoltranssorted = sortslices(allmoltrans, dims=2);
    sortindvec = round.(Int, vec(allmoltranssorted[end, :]));
    permute!(allmoltranslabels, sortindvec);
    allmoltrans = allmoltranssorted[1:end-1, :];

    # convert Euler angles into (alibi) rotation matrices
    allmoltranslations = deepcopy(allmoltrans[3*myMolData.nummols+1:end, :]);
    allmolrotations = Array{FT, 2}(undef, 3*myMolData.nummols, 3*numtrans);
    for tt=1:numtrans
        for nn=1:myMolData.nummols
            if (allmoltrans[3*nn-2:3*nn, tt] == zeros(FT, 3))
                allmolrotations[3*nn-2:3*nn, 3*tt-2:3*tt] = Matrix{FT}(I, 3, 3);
            else
                cosphi = cosd(allmoltrans[3*nn-2, tt]);
                sinphi = sind(allmoltrans[3*nn-2, tt]);
                costheta = cosd(allmoltrans[3*nn-1, tt]);
                sintheta = sind(allmoltrans[3*nn-1, tt]);
                cospsi = cosd(allmoltrans[3*nn, tt]);
                sinpsi = sind(allmoltrans[3*nn, tt]);

                allmolrotations[3*nn-2, 3*tt-2] = costheta*cospsi;
                allmolrotations[3*nn-2, 3*tt-1] = costheta*sinpsi;
                allmolrotations[3*nn-2, 3*tt] = -1.0*sintheta;
                allmolrotations[3*nn-1, 3*tt-2] = sinphi*sintheta*cospsi - cosphi*sinpsi;
                allmolrotations[3*nn-1, 3*tt-1] = sinphi*sintheta*sinpsi + cosphi*cospsi;
                allmolrotations[3*nn-1, 3*tt] = sinphi*costheta;
                allmolrotations[3*nn, 3*tt-2] = cosphi*sintheta*cospsi + sinphi*sinpsi;
                allmolrotations[3*nn, 3*tt-1] = cosphi*sintheta*sinpsi - sinphi*cospsi;
                allmolrotations[3*nn, 3*tt] = cosphi*costheta;
            end
        end
    end

    baselinetranslations = deepcopy(baselinetrans[3*myMolData.nummols+1:end]);
    baselinerotations = Array{FT, 2}(undef, 3*myMolData.nummols, 3);
    for nn=1:myMolData.nummols
        if (baselinetrans[3*nn-2:3*nn] == zeros(FT, 3))
            baselinerotations[3*nn-2:3*nn, :] = Matrix{FT}(I, 3, 3);
        else
            cosphi = cosd(baselinetrans[3*nn-2]);
            sinphi = sind(baselinetrans[3*nn-2]);
            costheta = cosd(baselinetrans[3*nn-1]);
            sintheta = sind(baselinetrans[3*nn-1]);
            cospsi = cosd(baselinetrans[3*nn]);
            sinpsi = sind(baselinetrans[3*nn]);

            baselinerotations[3*nn-2, 1] = costheta*cospsi;
            baselinerotations[3*nn-2, 2] = costheta*sinpsi;
            baselinerotations[3*nn-2, 3] = -1.0*sintheta;
            baselinerotations[3*nn-1, 1] = sinphi*sintheta*cospsi - cosphi*sinpsi;
            baselinerotations[3*nn-1, 2] = sinphi*sintheta*sinpsi + cosphi*cospsi;
            baselinerotations[3*nn-1, 3] = sinphi*costheta;
            baselinerotations[3*nn, 1] = cosphi*sintheta*cospsi + sinphi*sinpsi;
            baselinerotations[3*nn, 2] = cosphi*sintheta*sinpsi - sinphi*cospsi;
            baselinerotations[3*nn, 3] = cosphi*costheta;
        end
    end

    # after sorting, determine if each successive transformation affects a given molecule
    changedfrombefore = trues(myMolData.nummols, numtrans);
    for tt=2:numtrans
        for currmolind=1:myMolData.nummols
            if (allmoltrans[3*myMolData.nummols+3*currmolind-2, tt] ==
                allmoltrans[3*myMolData.nummols+3*currmolind-2, tt-1] &&
                allmoltrans[3*myMolData.nummols+3*currmolind-1, tt] ==
                allmoltrans[3*myMolData.nummols+3*currmolind-1, tt-1] &&
                allmoltrans[3*myMolData.nummols+3*currmolind, tt] ==
                allmoltrans[3*myMolData.nummols+3*currmolind, tt-1] &&
                allmoltrans[3*currmolind-2, tt] == allmoltrans[3*currmolind-2, tt-1] &&
                allmoltrans[3*currmolind-1, tt] == allmoltrans[3*currmolind-1, tt-1] &&
                allmoltrans[3*currmolind, tt] == allmoltrans[3*currmolind, tt-1])

                changedfrombefore[currmolind, tt] = false;

            end
        end
    end

    # horizontally concatenate all atomic coordinates over all (not just unique) molecules
    allAtomPos = Array{FT, 2}(undef, 3, myMolData.cumnumatomslist[end]);

    allAtomPos[:, 1:myMolData.numatomslist[1]] .=
    myMolData.uniqueMolDataArray[1].atompos0[:, :];

    for nn=2:myMolData.nummols

        currmolDupeInd = myMolData.dupemolinds[nn];
        if (currmolDupeInd == -1)
            allAtomPos[:,
                       myMolData.cumnumatomslist[nn-1]+1:myMolData.cumnumatomslist[nn]] .=
            myMolData.uniqueMolDataArray[myMolData.uniquemolinds[nn]].atompos0[:, :];
        else
            allAtomPos[:,
                       myMolData.cumnumatomslist[nn-1]+1:myMolData.cumnumatomslist[nn]] .=
            myMolData.uniqueMolDataArray[myMolData.dupemolinds[nn]].atompos0[:, :];
        end

    end

    # create and return the TransData
    myTransData = TransData{FT}(allmoltranslations, allmolrotations,
                                allmoltranslabels, changedfrombefore,
                                numtrans, allAtomPos, baselinetranslations,
                                baselinerotations);
    initializeAtomPos!(myTransData, myMolData);
    return myTransData;

end

"""

    initializeAtomPos!(myTransData::TransData{FT},
                       myMolData::MolSystem{<:OneMol{FT}}) where FT<:AbstractFloat

Initialize atomic coordinates across all molecules by first shifting coordinates of all
molecules to have centers of mass of 0, then rotating all molecular coordinates according to
baseline orientations of each, then translating to desired baseline centers of mass,
modifying `myTransData` in-place.

See also: [`transformAtomPos!`](@ref), [`untransformAtomPos!`](@ref)

"""
function initializeAtomPos!(myTransData::TransData{FT},
                            myMolData::MolSystem{<:OneMol{FT}}) where FT<:AbstractFloat
    startatom = 1;
    endatom = myMolData.numatomslist[1];
    for nn=1:myMolData.nummols
        if (nn > 1)
            startatom = myMolData.cumnumatomslist[nn-1] + 1;
            endatom = myMolData.cumnumatomslist[nn];
        end
        makeCOM0!(myTransData.allAtomPos, startatom, endatom);
        if (myTransData.baselineRot[3*nn-2:3*nn, :] != Matrix{FT}(I, 3, 3))
            rotatecoords!(myTransData.allAtomPos, myTransData.baselineRot[3*nn-2:3*nn, :],
                          startatom, endatom);
        end
        if (myTransData.baselineTrans[3*nn-2:3*nn] != zeros(FT, 3))
            translateCOM!(myTransData.allAtomPos, myTransData.baselineTrans[3*nn-2:3*nn],
                          startatom, endatom);
        end
    end
end

"""

    transformAtomPos!(myTransData::TransData{FT},
                      myMolData::MolSystem{<:OneMol{FT}},
                      tt::Integer) where FT<:AbstractFloat

Transform atomic coordinates across all molecules by first shifting coordinates of all
molecules to have centers of mass of 0, then rotating all molecular coordinates according to
orientations corresponding to transformation `tt` relative to the baseline, then translating
to desired centers of mass corresponding to transformation `tt` relative to the baseline,
modifying `myTransData` in-place.

See also: [`initializeAtomPos!`](@ref), [`untransformAtomPos!`](@ref)

"""
function transformAtomPos!(myTransData::TransData{FT},
                           myMolData::MolSystem{<:OneMol{FT}},
                           tt::Integer) where FT<:AbstractFloat
    startatom = 1;
    endatom = myMolData.numatomslist[1];
    for nn=1:myMolData.nummols
        if (nn > 1)
            startatom = myMolData.cumnumatomslist[nn-1] + 1;
            endatom = myMolData.cumnumatomslist[nn];
        end
        if (myTransData.baselineTrans[3*nn-2:3*nn] != zeros(FT, 3))
            makeCOM0!(myTransData.allAtomPos, startatom, endatom);
        end
        if (myTransData.rotArray[3*nn-2:3*nn, 3*tt-2:3*tt] != Matrix{FT}(I, 3, 3))
            rotatecoords!(myTransData.allAtomPos,
                          myTransData.rotArray[3*nn-2:3*nn, 3*tt-2:3*tt],
                          startatom, endatom);
        end
        if (myTransData.baselineTrans[3*nn-2:3*nn] .+
            myTransData.transArray[3*nn-2:3*nn, tt] != zeros(FT, 3))
            translateCOM!(myTransData.allAtomPos,
                          myTransData.baselineTrans[3*nn-2:3*nn] .+
                          myTransData.transArray[3*nn-2:3*nn, tt],
                          startatom, endatom);
        end
    end
end

"""

    untransformAtomPos!(myTransData::TransData{FT},
                        myMolData::MolSystem{<:OneMol{FT}},
                        tt::Integer) where FT<:AbstractFloat

Undo transformation of atomic coordinates across all molecules by first shifting coordinates
of all molecules to have centers of mass of 0, then rotating all molecular coordinates
according to orientations corresponding to the opposites of transformation `tt` to return to
the baseline, then translating back to the desired baseline centers of mass, modifying
`myTransData` in-place.

See also: [`initializeAtomPos!`](@ref), [`transformAtomPos!`](@ref)

"""
function untransformAtomPos!(myTransData::TransData{FT},
                             myMolData::MolSystem{<:OneMol{FT}},
                             tt::Integer) where FT<:AbstractFloat
    startatom = 1;
    endatom = myMolData.numatomslist[1];
    for nn=1:myMolData.nummols
        if (nn > 1)
            startatom = myMolData.cumnumatomslist[nn-1] + 1;
            endatom = myMolData.cumnumatomslist[nn];
        end
        if (myTransData.baselineTrans[3*nn-2:3*nn] .+
            myTransData.transArray[3*nn-2:3*nn, tt] != zeros(FT, 3))
            makeCOM0!(myTransData.allAtomPos, startatom, endatom);
        end
        if (myTransData.rotArray[3*nn-2:3*nn, 3*tt-2:3*tt] != Matrix{FT}(I, 3, 3))
            rotatecoords!(myTransData.allAtomPos,
                          transpose(myTransData.rotArray[3*nn-2:3*nn,3*tt-2:3*tt]),
                          startatom, endatom);
        end
        if (myTransData.baselineTrans[3*nn-2:3*nn] != zeros(FT, 3))
            translateCOM!(myTransData.allAtomPos,
                          myTransData.baselineTrans[3*nn-2:3*nn],
                          startatom, endatom);
        end
    end
end


text:
#! /usr/bin/env julia

using SerialPorts

function terminal(args)

    # arg parsing
    if length(args) < 4
        println("usage: listener.jl --port <port> --baud <baud>")
        exit(0)
    else
        p_i = findfirst(args, "--port")
        port = args[p_i+1]
        b_i = findfirst(args, "--baud")
        baud = parse(Int,args[b_i+1])
    end

    ser = SerialPort(port, baud)

    # setup handler at exit to close the serial port
    atexit_handler() = (println("exiting..."); close(ser))
    atexit(atexit_handler)

    println("Listening to $port @ $(baud)baud...")
    while true
        print(readavailable(ser))
    end
end

terminal(ARGS)

ext: jl


text:
###################################################################
#  Axis related definitions & functions
###################################################################

config_axis_spec = Dict(
  :axisWidth       => Real,
  :layer           => AbstractString,
  :offset          => Real,
  :grid            => Bool,
  :labels          => Bool,
  :labelAngle      => Real,
  :labelAlign      => AbstractString,
  :labelBaseline   => AbstractString,
  :labelMaxLength  => Real,
  :shortTimeLabels => Bool,
  :subdivide       => Real,
  :ticks           => Real,
  :tickPadding     => Real,
  :tickSize        => Real,
  :tickSizeMajor   => Real,
  :tickSizeMinor   => Real,
  :tickSizeEnd     => Real,
  :titleOffset     => Real,
  :titleMaxLength  => Real,
  :characterWidth  => Real
)

axis_help = """
General

- `axisWidth` (`Real`) : Width of the axis line
- `layer` (`AbstractString`) : A string indicating if the axis (and any gridlines) should be placed above or below the data marks.
- `offset` (`Real`) : The offset, in pixels, by which to displace the axis from the edge of the enclosing group or data rectangle.

Grid

- `grid` (`Bool`) : A flag indicate if gridlines should be created in addition to ticks. If `grid` is unspecified, the default value is `true` for ROW and COL. For X and Y, the default value is `true` for quantitative and time fields and `false` otherwise.

Labels

- `labels` (`Bool`) : Enable or disable labels.
- `labelAngle` (`Real`) : The rotation angle of the axis labels.
- `labelAlign` (`AbstractString`) : Text alignment for the Label.
- `labelBaseline` (`AbstractString`) : Text baseline for the label.
- `labelMaxLength` (`Real`) : Truncate labels that are too long.
- `shortTimeLabels` (`Bool`) : Whether month and day names should be abbreviated.

Ticks

- `subdivide` (`Real`) : If provided, sets the number of minor ticks between major ticks (the value 9 results in decimal subdivision). Only applicable for axes visualizing quantitative scales.
- `ticks` (`Real`) : A desired number of ticks, for axes visualizing quantitative scales. The resulting number may be different so that values are "nice" (multiples of 2, 5, 10) and lie within the underlying scale's range.
- `tickPadding` (`Real`) : The padding, in pixels, between ticks and text labels.
- `tickSize` (`Real`) : The size, in pixels, of major, minor and end ticks.
- `tickSizeMajor` (`Real`) : The size, in pixels, of major ticks.
- `tickSizeMinor` (`Real`) : The size, in pixels, of minor ticks.
- `tickSizeEnd` (`Real`) : The size, in pixels, of end ticks.

Title

- `titleOffset` (`Real`) : A title offset value for the axis.
- `titleMaxLength` (`Real`) : Max length for axis title if the title is automatically generated from the field's description. By default, this is automatically based on cell size and characterWidth property.
- `characterWidth` (`Real`) : Character width for automatically determining title max length.
"""

########### helper function for encoding()  ########
"""
Sets axis properties for the `axis` argument within `encoding()` :

```julia
encoding_x(..., axis=axis(subdivide=3), ...)
```

$axis_help
"""
axis(;properties...) = _mkdict(config_axis_spec, properties)

########### global level function ########
"""
Sets axis properties

`config_axis(;properties...)`

$axis_help
"""
config_axis(;properties...) = _mkvis((:config, :axis), config_axis_spec, properties)


text:
type Rnd <: Op; rgen; Rnd(;rgen=Uniform(0,1), o...)=new(rgen); end
Kenv.kdef(:rnd,Rnd)
ninputs(::Rnd)=0
canoverwrite(::Rnd)=false
back_reads_x(::Rnd)=false
back_reads_y(::Rnd)=false
infersize(::Rnd,y)=tuple(y)
forw(r::Rnd, y; o...)=rgen!(r.rgen, y)
back(r::Rnd, dy; o...)=error("rnd has no backward pass")


text:
#  Copyright 2016, Liye Zhang and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuSO
#
# See http://github.com/zhangliye/JuSO.jl
#############################################################################

function test_ut1()
  println("I am test_ut1()")
  return 1
end


text:
#= Valuation — Hard-coded Rates =#
@testset "Valuation" begin
    rates_a = ExchangeRateTable(:USD => 1.0, :CAD => 0.75)
    rates_b = ExchangeRateTable(Dict(
        :USD => 1.0,
        :EUR => 1.3,
        :GBP => 1.5))
    rates_c = ExchangeRateTable(Date(2015, 12, 02), Dict(
        :USD => 1.0,
        :EUR => 1.3,
        :GBP => 1.5,
        :CAD => 0.7,
        :JPY => 0.01))
    rates_d = ExchangeRateTable(
        Date(2015, 12, 02),
        :USD => 1.0,
        :CAD => 0.75)

    @test valuate(rates_a, :CAD, 21USD) == 28CAD
    @test valuate(rates_a, :CAD, Basket([21USD, 10CAD])) == 38CAD
    @test valuate(rates_b, :USD, 10USD) == 10USD
    @test valuate(rates_b, :EUR, 0.13USD) == 0.1EUR
    @test valuate(rates_c, :USD, Basket([USD, EUR, GBP])) == 3.8USD
    @test valuate(rates_c, :EUR, 100CAD) == 53.85EUR
    @test valuate(rates_c, :JPY, 1USD) == 100JPY
    @test valuate(rates_c, :JPY, -1USD) == -100JPY
    @test valuate(rates_c, :USD, Basket([200JPY, EUR])) == 3.3USD
    @test valuate(rates_c, :JPY, 0USD) == 0JPY
    @test valuate(rates_d, :CAD, 3.14CAD) == 3.14CAD
    @test occursin(":USD=>1.0", string(rates_d))
    @test occursin(":CAD=>0.75", string(rates_d))

    @test_throws KeyError valuate(rates_a, :JPY, 100USD)

    # Valuation at custom precision
    @testset "At Custom Precision" begin
        precisecad = Monetary(:CAD, 10000; precision=4)
        @test valuate(rates_d, Monetary{:USD, Int, 4}, CAD) ==
            Monetary(:USD, 7500; precision=4)
        @test valuate(rates_d, Monetary{:USD}, precisecad) == 0.75USD
    end

    # Valuation with custom representation
    @testset "With Custom Representation" begin
        bigcad = Monetary(:CAD, BigInt(100))
        @test valuate(rates_d, Monetary{:USD, BigInt, 4}, CAD) ==
            Monetary(:USD, BigInt(7500); precision=4)
        @test valuate(rates_d, :USD, bigcad) == 0.75USD
    end

    # Valuation to convert similar currencies
    @testset "Convert Similar Currencies" begin
        bigusd = Monetary(:USD; storage=BigInt)
        @test USD ≡ valuate(Dict(:USD => 1.0), :USD, bigusd)
        @test bigusd == valuate(Dict(:USD => 1.0), typeof(bigusd), USD)
    end
end

#= Valuation — ECB data =#
#= BROKEN TODO
@testset "ECB data" begin
    rates_e = ecbrates()  # most recent
    rates_f = ecbrates(Date(2015, 08, 05))  # fixed date
    rates_g = ecbrates(Date(2015, 08, 05))  # same as rates_f

    # test cache
    @test rates_e ≡ ecbrates()
    @test rates_f ≡ rates_g

    # test recent rates object
    @test rates_e.date + Dates.Day(7) > Date(now())
    @test !isempty(rates_e)
    @test isa(rates_e, ExchangeRateTable)
    @test valuate(rates_e, :USD, 1USD) == 1USD
    @test valuate(rates_e, :EUR, 20EUR) == 20EUR
    @test haskey(rates_e, :JPY)
    @test length(rates_e) > 10

    # test fixed rates object
    @test isa(rates_f, ExchangeRateTable)
    @test rates_f.date == Date(2015, 08, 05)
    @test valuate(rates_f, :AUD, 100USD) == 135.58AUD
end
=#


text:
"""
    lag(v::AbstractArray, n = 1; default = missing)

Return a `ShiftedArray` object, with underlying data `v`. The second argument gives the amount
to shift in each dimension. If it is an integer, it is assumed to refer to the first dimension.
`default` specifies a default value when you are out of bounds.

# Examples

```jldoctest lag
julia> v = [1, 3, 5, 4];

julia> lag(v)
4-element ShiftedArrays.ShiftedArray{Int64,Missing,1,Array{Int64,1}}:
  missing
 1
 3
 5

julia> w = 1:2:9
1:2:9

julia> s = lag(w, 2)
5-element ShiftedArrays.ShiftedArray{Int64,Missing,1,StepRange{Int64,Int64}}:
  missing
  missing
 1
 3
 5

julia> copy(s)
5-element Array{Union{Int64, Missing},1}:
  missing
  missing
 1
 3
 5

julia> v = reshape(1:16, 4, 4);

julia> s = lag(v, (0, 2))
4×4 ShiftedArrays.ShiftedArray{Int64,Missing,2,Base.ReshapedArray{Int64,2,UnitRange{Int64},Tuple{}}}:
 missing  missing  1  5
 missing  missing  2  6
 missing  missing  3  7
 missing  missing  4  8
```
"""
lag(v::AbstractArray, n = 1; default = missing) = ShiftedArray(v, n; default = default)

"""
    lead(v::AbstractArray, n = 1; default = missing)

Return a `ShiftedArray` object, with underlying data `v`. The second argument gives the amount
to shift negatively in each dimension. If it is an integer, it is assumed to refer
to the first dimension. `default` specifies a default value when you are out of bounds.

# Examples

```jldoctest lead
julia> v = [1, 3, 5, 4];

julia> lead(v)
4-element ShiftedArrays.ShiftedArray{Int64,Missing,1,Array{Int64,1}}:
 3
 5
 4
  missing

julia> w = 1:2:9
1:2:9

julia> s = lead(w, 2)
5-element ShiftedArrays.ShiftedArray{Int64,Missing,1,StepRange{Int64,Int64}}:
 5
 7
 9
  missing
  missing

julia> copy(s)
5-element Array{Union{Int64, Missing},1}:
 5
 7
 9
  missing
  missing

julia> v = reshape(1:16, 4, 4);

julia> s = lead(v, (0, 2))
4×4 ShiftedArrays.ShiftedArray{Int64,Missing,2,Base.ReshapedArray{Int64,2,UnitRange{Int64},Tuple{}}}:
  9  13  missing  missing
 10  14  missing  missing
 11  15  missing  missing
 12  16  missing  missing
```
"""
lead(v::AbstractArray, n = 1; default = missing) = ShiftedArray(v, map(-, n); default = default)


text:
function temp_pkg_dir(fn::Function)
    local env_dir
    local old_load_path
    local old_depot_path
    try
        old_load_path = copy(LOAD_PATH)
        old_depot_path = copy(DEPOT_PATH)
        empty!(LOAD_PATH)
        empty!(DEPOT_PATH)
        mktempdir() do env_dir
            mktempdir() do depot_dir
                pushfirst!(LOAD_PATH, env_dir)
                pushfirst!(DEPOT_PATH, depot_dir)
                # Add the standard library paths back
                vers = "v$(VERSION.major).$(VERSION.minor)"
                push!(LOAD_PATH, abspath(Sys.BINDIR, "..", "local", "share", "julia", "site", vers))
                push!(LOAD_PATH, abspath(Sys.BINDIR, "..", "share", "julia", "site", vers))
                fn(env_dir)
            end
        end
    finally
        empty!(LOAD_PATH)
        empty!(DEPOT_PATH)
        append!(LOAD_PATH, old_load_path)
        append!(DEPOT_PATH, old_depot_path)
    end
end

function cd_tempdir(f)
    mktempdir() do tmp
        cd(tmp) do
            f(tmp)
        end
    end
end

isinstalled(pkg) = Base.locate_package(Base.PkgId(pkg.uuid, pkg.name)) !== nothing

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
""" Computes divided differences of tabulated data n points x and
    for each of the n points, up to m derivatives  dy
"""
function DividedDifference2(x,dy)  # generate array of divided differences
    n, m = size(dy)                # n data points, m derivatives (0 to m-1)
    dd = zeros(typeof(x[1]),n*m,n*m)               # matrix of divided differences
    z = zeros(typeof(x[1]),n*m)
    k=1;
    for i = 1:n                    # n data points
        for j = 1:m                # m derivatives (0 to m-1) at each point
            k = (i-1)*m + j        # row index
            z[k] = x[i]
            dd[k,1] = dy[i,1]      # 0th divided difference in first column
            for l = 2:k            # column index for the remaining columns
                if dd[k,l-1] == dd[k-1,l-1]  # left and top-left neighbors are repeated
                    if l<=m
                        dd[k,l] = dy[i,l] / factorial(l-1)
                    else dd[k,l] = big(0.0)
                    end
                else
                    dd[k,l] = (dd[k,l-1] - dd[k-1,l-1]) / (z[k]-z[k-l+1])
                end
            end
        end
    end
    return dd
end


text:
# include("../../src/SFFM.jl")
using LinearAlgebra, Plots, JLD2, GLM, SFFM

## define the model(s)
include("exampleModelDef.jl")

## load sims
@load pwd()*"/examples/paperNumerics/dump/sims.jld2" sims

## analytic X distribution for comparison
# construction
Ψₓ = SFFM.PsiFunX( approxModel)
ξₓ = SFFM.MakeXiX( approxModel, Ψₓ)
pₓ, πₓ, Πₓ, Kₓ = SFFM.StationaryDistributionX( approxModel, Ψₓ, ξₓ)

## section 4.5: error for approximation of stationary distribution of X
Basis = "lagrange"
Δs = [1.6; 0.8; 0.4; 0.2; 0.1; 0.05]
NBasesRange = [1; 2; 3; 4]

πnorms = zeros(length(Δs), length(NBasesRange))
Ψnorms = zeros(length(Δs), length(NBasesRange))
times = Array{Any}(undef, length(Δs), length(NBasesRange))
gctimes = Array{Any}(undef, length(Δs), length(NBasesRange))
mems = Array{Any}(undef, length(Δs), length(NBasesRange))
alloc = Array{Any}(undef, length(Δs), length(NBasesRange))
approxSpec = Array{Any}(undef, length(Δs), length(NBasesRange))

for d = 1:length(Δs), n = 1:length(NBasesRange)
    if d + 2 * n <= 10
        # define the mesh for each iteration
        nBases = NBasesRange[n]
        Δtemp = Δs[d]
        println("mesh details; h = ", Δtemp, " NBases = ", nBases)

        # collect time and memory stats
        ~, times[d, n], mems[d, n], gctimes[d, n], alloc[d, n] = @timed begin
            Nodes = collect(approxBounds[1, 1]:Δtemp:approxBounds[1, 2])
            mesh = SFFM.DGMesh(
                approxModel,
                Nodes,
                nBases,
                Basis = Basis,
            )
            approxSpec[d, n] = (Δtemp, nBases, SFFM.TotalNBases(mesh) * SFFM.NPhases(approxModel))

            # compute the marginal via DG
            All = SFFM.MakeAll( approxModel, mesh, approxType = "projection")
            Ψ = SFFM.PsiFun( All.D)
            # @btime SFFM.PsiFun( $All.D)

            # the distribution of X when Y first returns to 0
            ξ = SFFM.MakeXi(All.B, Ψ)
            # @btime SFFM.MakeXi($All.B.BDict, $Ψ)

            marginalX, p, K = SFFM.MakeLimitDistMatrices(
                All.B,
                All.D,
                All.R.RDict,
                Ψ,
                ξ,
                mesh,
                approxModel,
            )
            # @btime SFFM.MakeLimitDistMatrices(
            #     $All.B.BDict,
            #     $All.D,
            #     $All.R.RDict,
            #     $Ψ,
            #     $ξ,
            #     $mesh,
            #     $approxModel,
            # )
        end
        # convert marginalX to a distribution for analysis
        DGStationaryDist = SFFM.Coeffs2Dist(
            approxModel,
            mesh,
            marginalX,
            SFFM.SFFMProbability,
        )

        ## stationary distirbution stuff
        # evaluate the analytic result given the mesh
        analyticX = SFFM.SFFMProbability(
            [pₓ[:]; 0; 0],
            Πₓ(Matrix(mesh.Nodes[2:end]')) - Πₓ(Matrix(mesh.Nodes[1:end-1]')),
            mesh.Nodes[1:end-1] + SFFM.Δ(mesh) / 2,
        )

        # save them
        πnorms[d, n] = SFFM.starSeminorm( DGStationaryDist, analyticX)

        ## now do Ψ
        ## turn sims into a cdf
        simprobs =
            SFFM.Sims2Dist( simModel, mesh, sims, SFFM.SFFMProbability)

        # do DG for Ψ
        theNodes = SFFM.CellNodes(mesh)[:, convert(Int, ceil(5 / Δtemp))]
        basisValues = zeros(length(theNodes))
        for n = 1:length(theNodes)
            basisValues[n] =
                prod(5.0 .- theNodes[[1:n-1; n+1:end]]) ./
                prod(theNodes[n] .- theNodes[[1:n-1; n+1:end]])
        end
        initpm = [
            zeros(sum(approxModel.C .<= 0)) # LHS point mass
            zeros(sum(approxModel.C .>= 0)) # RHS point mass
        ]
        V=SFFM.vandermonde(NBases(mesh))
        initprobs = zeros(Float64, SFFM.NBases(mesh), SFFM.NIntervals(mesh), SFFM.NPhases(approxModel))
        initprobs[:, convert(Int, ceil(5 / Δtemp)), 3] =
            basisValues' * V.V * V.V' .* 2 / Δtemp
        initdist =
            SFFM.SFFMDensity(initpm, initprobs, SFFM.CellNodes(mesh)) # convert to a distribution object so we can apply Dist2Coeffs
        # convert to Coeffs α in the DG context
        x0 = SFFM.Dist2Coeffs( approxModel, mesh, initdist)
        # the initial condition on Ψ is restricted to + states so find the + states
        plusIdx = [
            mesh.Fil["p+"]
            repeat(mesh.Fil["+"]', SFFM.NBases(mesh), 1)[:]
            mesh.Fil["q+"]
        ]
        # get the elements of x0 in + states only
        x0 = x0[plusIdx]'

        # compute x-distribution at the time when Y returns to 0
        w = x0 * Ψ
        # this can occur in - states only, so find the - states
        minusIdx = [
            mesh.Fil["p-"]
            repeat(mesh.Fil["-"]', SFFM.NBases(mesh), 1)[:]
            mesh.Fil["q-"]
        ]
        # then map to the whole state space for analysis
        z = zeros(
            Float64,
            SFFM.NBases(mesh) * SFFM.NIntervals(mesh) * SFFM.NPhases(approxModel) +
            sum(approxModel.C .<= 0) +
            sum(approxModel.C .>= 0),
        )
        z[minusIdx] = w

        # convert to a distribution object
        DGΨProbs = SFFM.Coeffs2Dist(
            approxModel,
            mesh,
            z,
            SFFM.SFFMProbability,
        )

        Ψnorms[d, n] = SFFM.starSeminorm( DGΨProbs, simprobs)
    else
        πnorms[d, n] = NaN
        Ψnorms[d, n] = NaN
        times[d, n], mems[d, n], gctimes[d, n], alloc[d, n] = NaN, NaN, NaN, NaN
    end
end

let p = plot()
    types = [:solid,:dash,:dashdot,:dashdotdot]
    for n = 1:length(NBasesRange)
        p = plot!(
            Δs,
            πnorms[:, n],
            xaxis = :log,
            yaxis = :log,
            xlabel = "h",
            ylabel = "error",
            label = "N_k = " * string(NBasesRange[n]),
            legend = :outertopright,
            linestyle = types[n],
            seriestype = :line,
            markershape = :auto,
            title = "Error",
            grid = false,
        )
    end
    display(p)
    # savefig(pwd()*"/examples/paperNumerics/dump/piErrorVsDelta.png")
end

let p = plot()
    types = [:solid,:dash,:dashdot,:dashdotdot,:solid,:dash]
    for d = 1:length(Δs)
        p = plot!(
            NBasesRange,
            πnorms[d, :],
            yaxis = :log,
            xlabel = "N_k",
            ylabel = "error",
            label = "h = " * string(Δs[d]),
            legend = :outertopright,
            linestyle = types[d],
            seriestype = :line,
            markershape = :auto,
            title = "Error",
            grid = false,
        )
    end
    display(p)
    # savefig(pwd()*"/examples/paperNumerics/dump/piErrorVsNBases.png")
end

begin
    n = 1
    data = log.(πnorms[:,n])
    ind = .!isnan.(data)
    data = data[ind]
    println("rate of convergence for ", n, " basis functions is ")
    display(
        lm([ones(length(Δs[ind])) log.(Δs[ind])], data)
    )
    println("")

    n = 2
    data = log.(πnorms[:,n])
    ind = .!isnan.(data)
    data = data[ind]
    println("rate of convergence for ", n, " basis functions is ")
    display(
        lm([ones(length(Δs[ind])) log.(Δs[ind])], data)
    )
    println("")

    n = 3
    data = log.(πnorms[:,n])
    ind = .!isnan.(data)
    ind[4] = false
    data = data[ind]
    data = data
    println("rate of convergence for ", n, " basis functions is ")
    display(
        lm([ones(length(Δs[ind])) log.(Δs[ind])], data)
    )
    println("")
end

for d = 1:4
    data = log.(πnorms[d,:])
    ind = .!isnan.(data)
    data = data[ind]
    println("rate of convergence for h = ", Δs[d])
    display(
        lm([ones(length(NBasesRange[ind])) NBasesRange[ind]], data)
    )
    println("")
end

mp =
    (log.(πnorms[1:end-1, :]) .- log.(πnorms[2:end, :])) ./
    (log.(Δs[1:end-1]) - log.(Δs[2:end]))
np =
    (log.(πnorms[:, 1:end-1]) .- log.(πnorms[:, 2:end])) ./
    ((NBasesRange[1:end-1]) - (NBasesRange[2:end]))'

let q = plot()
    for n = 1:length(NBasesRange)
        q = plot!(
            Δs,
            Ψnorms[:, n],
            xaxis = :log,
            yaxis = :log,
            xlabel = "h",
            ylabel = "error",
            label = "NBases = " * string(NBasesRange[n]),
            legend = :outertopright,
            seriestype = :line,
            markershape = :diamond,
            title = "Error for Ψ",
        )
    end
    display(q)
    # savefig(pwd()*"/examples/paperNumerics/dump/psiErrorVsDelta.png")
end

mq =
    (log.(Ψnorms[1:end-1, :]) .- log.(Ψnorms[2:end, :])) ./
    (log.(Δs[1:end-1]) - log.(Δs[2:end]))

stats = Array{Any}(undef,length(Δs)+1,length(NBasesRange)+1)
for n in 1:length(NBasesRange), d in 1:length(Δs)
    if d + 2 * n <= 10
        stats[d+1,n+1] = (πnorms[d,n], times[d,n], mems[d,n]/(1024*1024), approxSpec[d,n][3])
    else
        stats[d+1,n+1] = (0,0,0,0)
    end
end
stats[1,2:end] = NBasesRange
stats[2:end,1] = Δs
display(stats)
using NumericIO
innards = stats[2:end,2:end]
asciiexponentfmt = NumericIO.IOFormattingExpNum(
	"\\times10^{", false, '+', '-', NumericIO.ASCII_SUPERSCRIPT_NUMERALS
)
fmt = NumericIO.IOFormattingReal(asciiexponentfmt,
	ndigits=3, decpos=0, decfloating=true, eng=false, minus='-', inf="Inf"
)

begin
    println("\\begin{tabular}{ | c | l | l | l | l | l | }")
    println("\\hline")
    println("\\multirow{2}{*}{\\(\\Delta\\)} & \\multirow{2}{*}{}& \\multicolumn{4}{c|}{Number of basis functions, \\(N_k\\)} \\\\ \\cline{3-6}
          &      & 1            & 2           & 3           & 4           \\\\\\hline"
        )
    for d in 1:length(Δs)
        println(Δs[d], " &")
        println("
        \\(\\begin{array}{l}
            \\mbox{error}
            \\\\ \\mbox{time (sec)}
            \\\\ \\mbox{memory (MB)}
            \\\\ n_\\varphi
        \\end{array}\\)"
        )
        for n in 1:length(NBasesRange)
            println("&")
            if d + 2 * n <= 10
                println("
                \\(\\begin{array}{l}
                    ", formatted(innards[d,n][1],fmt),"}
                    \\\\",round(innards[d,n][2],sigdigits=3),"
                    \\\\",round(innards[d,n][3],sigdigits=3),"
                    \\\\",convert(Int,round(innards[d,n][4],digits=0)),"
                \\end{array}\\)"
                )
            else
                println(" - ")
            end
        end
        println("\\\\ \\hline")
    end
    println("\\end{tabular}")
end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
module POMDPIDMModel

using Random
using Distributions
using LinearAlgebra

using POMDPs
using POMDPModelTools
using POMDPPolicies
using POMDPSimulators
using BeliefUpdaters
using DiscreteValueIteration

using Geo
using Plots
using Measures

import Base: round, ceil
import POMDPs: update, stateindex
import POMDPModelTools: update_info
import QMDP: solve, QMDPSolver

export
    CarSt,
    Sts,
    RoutesGenerator,
    DrivePOMDP,
    solve,
    Simulator,
    Visualisation,
    stateindex,
    SSRound,
    sv_boundry!,
    IDMtransit,
    show_DB,
    action,
    IDMtransit

include("Structs.jl")
include("Route.jl")
include("IDM.jl")
include("NormalDistribution.jl")
include("DPConstructor.jl")

include("Action.jl")
include("State.jl")
include("Transition.jl")
include("Observation.jl")
include("Reward.jl")
include("QMDPMod.jl")
include("Interface.jl")

include("Visualisation.jl")

end # module


text:
run(`mkfifo /tmp/pwd-fifo`)
prompt = "sudo need your password, please run [read -s p && echo \$p > /tmp/pwd-fifo] in another shell.\n"
sudo = @async run(`bash -c "sudo -S -p '' touch /tmp/tx < /tmp/pwd-fifo"`)
print(prompt)
wait(sudo)


text:
module tictactoe

using LinearAlgebra
using DelimitedFiles

include("create_game.jl")
include("minimax.jl")

export playgame, playgame_ai

end

max_stars_repo_name: HariharanJayashankar/tictactoe

max_issues_repo_name: HariharanJayashankar/tictactoe

max_forks_repo_name: HariharanJayashankar/tictactoe


text:
function dimComb(dims::Dim)
  Set(filter(x -> x[1] < x[2], collect(Iterators.product(Dim(1):(dims-1), Dim(2):dims))))
end

max_stars_repo_name: gabysbrain/multid-slicing

max_issues_repo_name: gabysbrain/multid-slicing


text:
abstract type Clock <: System end
timeunit(::Type{<:Clock}) = u"hr"
@system Clock{timeunit = timeunit(Clock)} begin
    context ~ ::Nothing
    config ~ ::Config(override)
    init => 0 ~ preserve(unit=timeunit, parameter)
    step => 1 ~ preserve(unit=timeunit, parameter)
    time => nothing ~ advance(init=init, step=step, unit=timeunit)
    tick => nothing ~ advance::int
end

abstract type DailyClock <: Clock end
timeunit(::Type{<:DailyClock}) = u"d"
@system DailyClock{timeunit = timeunit(DailyClock)}(Clock) <: Clock

timeunit(c::C) where {C<:Clock} = timeunit(C)

export Clock


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

using Markdown
using InteractiveUtils

# ╔═╡ 89edd424-485c-11ec-09eb-1712ed5a188b
a = 1+2

# ╔═╡ 086e8c68-657c-42e4-b7e0-f2c4c51106d3
b=a+1

# ╔═╡ ca5bc530-09d8-44e8-a19e-6e9b4f9ca307
This is not code, it is **markdown**

# ╔═╡ Cell order:
# ╠═89edd424-485c-11ec-09eb-1712ed5a188b
# ╠═086e8c68-657c-42e4-b7e0-f2c4c51106d3
# ╠═ca5bc530-09d8-44e8-a19e-6e9b4f9ca307

max_stars_repo_name: sylvaticus/JuliaMLCourse

max_issues_repo_name: sylvaticus/JuliaMLCourse

max_forks_repo_name: sylvaticus/JuliaMLCourse


text:
module Estimator

import Eive.CGA: cga

function eive(
    dirtyx::Array{T,1},
    y::Array{T,1},
    otherx::Union{Nothing, Array{T,2}, Array{T, 1}},
    popsize::Int = 50,
    numdummies::Int = 10,
) where {T<:Real}

    if isnothing(otherx)
        eivewithoutotherx(dirtyx, y, popsize, numdummies)
    else
        return eivewithotherx(dirtyx, y, otherx, popsize, numdummies)
    end
end

function eivewithotherx(
    dirtyx::Array{T,1},
    y::Array{T,1},
    otherx::Union{Array{T,2}, Array{T, 1}},
    popsize::Int = 50,
    numdummies::Int = 10,
) where {T<:Real}

    n = length(dirtyx)
    myones = ones(Float64, n)
    chsize = n * numdummies

    function costfn(bits::Array{Int,1})
        auxX = reshape(bits, n, numdummies)
        betas = auxX \ dirtyx
        cleanX = auxX * betas

        X = hcat(myones, cleanX, otherx)

        outerbetas = X \ y
        res = y .- X * outerbetas
        return sum(res .^ 2.0)
    end

    finalbits = cga(chsize = chsize, costfunction = costfn, popsize = popsize)

    auxX = reshape(finalbits, n, numdummies)
    betas = auxX \ dirtyx
    cleanX = auxX * betas

    X = hcat(myones, cleanX, otherx)

    outerbetas = X \ y

    return outerbetas
end

function eivewithoutotherx(
    dirtyx::Array{T,1},
    y::Array{T,1},
    popsize::Int = 50,
    numdummies::Int = 10,
) where {T<:Real}

    n = length(dirtyx)
    myones = ones(Float64, n)
    chsize = n * numdummies

    function costfn(bits::Array{Int,1})
        auxX = reshape(bits, n, numdummies)
        betas = auxX \ dirtyx
        cleanX = auxX * betas

        X = hcat(myones, cleanX)

        outerbetas = X \ y
        res = y .- X * outerbetas
        return sum(res .^ 2.0)
    end

    finalbits = cga(chsize = chsize, costfunction = costfn, popsize = popsize)

    auxX = reshape(finalbits, n, numdummies)
    betas = auxX \ dirtyx
    cleanX = auxX * betas

    X = hcat(myones, cleanX)

    outerbetas = X \ y

    return outerbetas
end

export eive

end # end of module


text:
function read_navheader(f::IOStream,
            nav_hstrings_set::Set{String}=nav_hstrings_set,
        )::Dict{String,String}
    header_hs = Dict{String,String}()
    if get_header!(header_hs, f, nav_hstrings_set)
        haskey(header_hs, ohsVersion) ||
        throw(error("header error, not $ohsVersion"))

        return header_hs
    end
    throw(error("read header error!"))
end

function read_navheader(filename::String)
    open(filename, "r") do f
        return read_navheader(f)
    end
end

function parse_navheader(header_hs::Dict{String,String})

    version     = parse_nhsVersion(get(header_hs, nhsVersion, ""))
    runBy       = parse_nhsRunBy(get(header_hs, nhsRunBy, ""))
    comment     = parse_nhsComment(get(header_hs, nhsComment, ""))
    ionoCorr    = parse_nhsIonoCorr(get(header_hs, nhsIonoCorr, ""))
    corrSysTime = get(header_hs, nhsCorrSysTime, "")
    deltaUTC    = get(header_hs, nhsDeltaUTC, "")
    dUTC        = get(header_hs, nhsDUTC, "")
    if version[1] < 3
        if version[4] == sstGlonass
            timeSysCorr = parse_nhsCorrSysTime(corrSysTime) # R2.11GPS
        elseif version[4] == sstGPS
            timeSysCorr = parse_nhsDeltaUTC(deltaUTC)
        elseif version[4] == sstGeosync
            timeSysCorr = parse_nhsDUTC(dUTC) # R2.11GEO

        end
    else
        timeSysCorr = parse_nhsTimeSysCorr(get(header_hs, nhsTimeSysCorr, ""))
    end
    leapSeconds = parse_nhsLeapSeconds(get(header_hs, nhsLeapSeconds, "")) # R2.10GLO
    ionAlpha    = parse_nhsIonAlpha(get(header_hs, nhsIonAlpha, "")) # R2.11
    ionBeta     = parse_nhsIonBeta(get(header_hs, nhsIonBeta, ""))
    eoH         = parse_nhsEoH(get(header_hs, nhsEoH, ""))

    RinexNavHeader(version,runBy,comment,
                ionoCorr,timeSysCorr,
                leapSeconds,corrSysTime,deltaUTC,dUTC,
                ionAlpha,ionBeta,eoH )
end
function parse_navd(str::SubString)
    if str[16] == ' ' return 0.0 end

    # b = parse(Float64, str[1:15])
    # e = parse(Float64, str[17:19])
    # b * 10^e
    parse(Float64, replace(str, 'D' => 'e'))
end

function read_navline(line::SubString, d::Array{Float64,1})
    len = length(line)
    start = 1
    # @show line
    while start < len
        push!(d, parse_navd(line[start:start + 18]))
        start += 19
    end
    return
end

const pvaSSTs = Set([sstGlonass,sstGeosync])

function parse_2navdt(sbdt::SubString)
    # @show sbdt
    y  = parse(Int64, sbdt[1:2])
    if y < 70
        y += 2000
    else
        y += 1900
    end
    m  = parse(Int64, sbdt[4:5])
    d  = parse(Int64, sbdt[7:8])
    H  = parse(Int64, sbdt[10:11])
    M  = parse(Int64, sbdt[13:14])
    S  = parse(Float64, sbdt[15:19])
    ds = Millisecond(0)
    if S >= 60.0
        ds = Millisecond(S * 1000)
        S = 0.0
    end
    DateTime(y, m, d, H, M, S) + ds
end

function parse_3navdt(sbdt::SubString)
    # @show sbdt
    y  = parse(Int64, sbdt[1:4])

    m  = parse(Int64, sbdt[6:7])
    d  = parse(Int64, sbdt[9:10])
    H  = parse(Int64, sbdt[12:13])
    M  = parse(Int64, sbdt[15:16])
    S  = parse(Int64, sbdt[18:19])
    ds = Second(0)
    if S >= 60
        ds = Second(S)
        S = 0
    end
    DateTime(y, m, d, H, M, S) + ds
end

function read_navdata(nh::RinexNavHeader, f::IOStream, fc::IOStream)

    if nh.version[1] > 2.99
        s = sstUnknown
        index_plus = 1
        check_index = [4,9,12,15,18,21]
    else
        s = nh.version[4]
        index_plus = 0
        check_index = [3, 6, 9, 12, 15, 18]
    end
    index_health = 25
    rnd = Dict{SatID,EphEpochStore}()

    while !eof(f)
        time_existed = false
        line = readline(f)
        for i in check_index
            line[i] == ' ' ||
                throw(error("no valid epoch line:\n|$line|"))
        end
        if index_plus == 0
            id = parse(Int64, line[1:2])
            si = SatID(s, id)
            dto = parse_2navdt(SubString(line, 4, 22), )
        else
            si = SatID(line[1:3])
            s = si.ss
            dto = parse_3navdt(SubString(line, 5, 23), )
        end

        # @show stt2tst[s]
        # @show diff2UTC(dto, stt2tst[s])
        dt = dto + diff2UTC(dto, stt2tst[s])
        if in(s, pvaSSTs)
            n = 3
            index_health = 7
        else
            n = 7
            index_health = 25
        end
        if haskey(rnd, si)
            # @show si "haskey rnd si"
            ees = rnd[si]
            if haskey(ees.eph, dt)
                # @show ees.keys, dt, line
                # @show dt "haskey ees.eph dt"
                # @show line
                existed_arr = ees.eph[dt].d
                time_existed = true
                # throw(error("duplicate epoch? $si-$dt"))
            end
            # @show dt "this dt"
        else
            ees = EphEpochStore(si, Array{DateTime,1}(),
                            Dict{DateTime,EphEpoch}(),)
            rnd[si] = ees
        end
        arr = Array{Float64,1}()
        sizehint!(arr, 3 + n * 4)
        read_navline(SubString(line, 23 + index_plus,
                                        79 + index_plus),
                    arr, )
        for i in 1:n
            line = readline(f)
            len = length(line)
            read_navline(SubString(line, 4 + index_plus, len),
                            arr, )
        end

        arr[index_health] == 0 || continue
        if second(dto) != 0 continue end

        if s == sstGlonass
            if minute(dto) != 15 && minute(dto) != 45
                continue
            end
        elseif s == sstGalileo
            if minute(dto) % 10 != 0 continue end
        end

        if time_existed
            # @show si "time_existed"
            # @show dt "time_existed"
            # @show ees.eph[dt].d
            time_existed = false
            for (iien, (iex, inow)) in enumerate(zip(existed_arr, arr, ))
                if iex == inow continue end
                iienstr = @sprintf("%02d",iien)
                warnline = "$si $dto NO.$iienstr data different: $iex, $inow"
                write(fc, warnline, "\n")
                if s == sstGlonass
                    # @warn warnline
                    continue
                end
                if s == sstGalileo
                    # @warn warnline
                    continue
                #  Test this by comparing Toe values. OrbitEphStore.cpp Line 276
                end
                if s == sstBeiDou
                    # @warn warnline
                    continue
                end
                # throw(error(warnline))
                @warn warnline
            end
            continue
        end
        push!(ees.keys, dt)
        sort!(ees.keys)
        ees.eph[dt] = EphEpoch(si, dt, arr)
    end
    # @show rnd
    RinexNavData(rnd)
end

function load_nav(filename::String)
    open(filename, "r") do f
        check_file = string(filename, ".check")
        open(check_file, "w") do fc
            nhd = read_navheader(f)
            nh  = parse_navheader(nhd)
            nd  = read_navdata(nh, f, fc)
            return RinexNav(basename(filename), nh, nd)
        end
    end
end

function merge_navs(navs::Array{RinexNav,1})

    data = Dict{SatID,EphEpochStore}()
    mheaders = Array{RinexNavHeader,1}()

    for inav in navs
        push!(mheaders, inav.header)
        for (isat, ees) in inav.data.data
            if haskey(data, isat)
                this_data = data[isat]
                this_keys = this_data.keys
                this_eph = this_data.eph
            else
                data[isat] = EphEpochStore(ees.si, ees.keys, ees.eph)
                continue
            end
            for (idt, iee) in ees.eph
                if haskey(this_eph, idt)
                    continue
                end
                this_eph[idt] = iee
                push!(this_keys, idt)
            end
            sort!(this_keys)
        end
    end
    mheaders, data
end

function load_navs(nav_files::Array{String,1})
    navs = Array{RinexNav,1}(undef, length(nav_files))
    for (i, this_nav_file) in enumerate(nav_files)
        navs[i] = load_nav(this_nav_file)
    end
    merge_navs(navs)
end

function get_freq(navdata::Dict{SatID,EphEpochStore},
    basefrequencies::Dict{SatelliteSystemType,Dict{CarrierBand,Tuple{Float64,Float64}}}=basefrequencies)
    freq = Dict{SatID,Dict{CarrierBand,Float64}}()
    for (isat, ees) in navdata

        freqn = 0.0
        if isat.ss == sstGlonass
            freqng = NaN
            for (idt, iee) in ees.eph
                if isnan(freqng) freqng = iee.d[11]
                elseif freqng == iee.d[11] continue
                else @warn string("$isat frequency number is different? $freqng ", iee[11])
                end
            end
            if !isnan(freqng) freqn = freqng end
        end
        this_freq = Dict{CarrierBand,Float64}()
        for (icb, ibf) in basefrequencies[isat.ss]
            this_freq[icb] = ibf[1]  + ibf[2] * freqn
        end
        freq[isat]  = this_freq
    end
    freq
end


text:
const TEST_NAME = "AESNI1x128"

include("common.jl")

using RandomNumbers.Random123

r = AESNI1x(123)

test_all(r, 100_000_000)


text:
"""
        game(decisionMatrix; verbose = false)

    Solve the zero-sum game.

# Arguments:
 - `decisionMatrix::Array{Float,2}`: n × m matrix of gain values for row player.
 - `verbose::Bool`: Logical value indicating whether the optimization logs are shown or not. Default is false.

# Description
game() solves zero-sum games using a gain matrix designed for the row player. The game matrix has n rows
and m columns for n and m strategies of row and column players, respectively.

# Output
- `::GameResult`: GameResult object that holds mixed strategy probabilities for row player and the value of game.
If a pure strategy exists for the row player, than the propabilities vector is a one-hat vector.

# Examples
```julia-repl
julia> # Rock & Paper & Scissors game
julia> mat = [0 -1 1; 1 0 -1; -1 1 0]
3×3 Array{Int64,2}:
  0  -1   1
  1   0  -1
 -1   1   0

julia> dm = makeDecisionMatrix(mat)
3×3 DataFrame
 Row │ Crt1     Crt2     Crt3
     │ Float64  Float64  Float64
─────┼───────────────────────────
   1 │     0.0     -1.0      1.0
   2 │     1.0      0.0     -1.0
   3 │    -1.0      1.0      0.0

julia> result = game(dm);

julia> result.row_player_probabilities
3-element Array{Float64,1}:
 0.3333333333333333
 0.33333333333333337
 0.3333333333333333

julia> result.value
0.0

```

# References
Zhou, Hai-Jun. "The rock–paper–scissors game." Contemporary Physics 57.2 (2016): 151-163.
"""
    function game(decisionMatrix::DataFrame; verbose::Bool=false)::GameResult

    newDecisionMatrix = similar(decisionMatrix)

    nrow, ncol = size(decisionMatrix)

    minmat = minimum(decisionMatrix)
    modified = false

    if minmat < 0
        modified = true
        newDecisionMatrix = decisionMatrix .- minmat
    end

    # dm = convert(Matrix, newDecisionMatrix)
    dm = Matrix(newDecisionMatrix)

    model = Model(GLPK.Optimizer);
    MOI.set(model, MOI.Silent(), !verbose)

    @variable(model, x[1:nrow])
    @objective(model, Min, sum(x))

    for i in 1:nrow
        @constraint(model, sum(x[1:nrow] .* dm[:, i]) >= 1)
    end

    for i in 1:nrow
        @constraint(model, x[i] >= 0)
    end

    optimize!(model);

    values = JuMP.value.(x)

    gamevalue = 1 / objective_value(model)
    p = gamevalue * values

    if modified
        gamevalue -= abs(minmat)
    end

result = GameResult(
        p,
        gamevalue
    )

    return result

end

max_stars_repo_name: jbytecode/jMCDM

max_issues_repo_name: jbytecode/jMCDM

max_forks_repo_name: jbytecode/jMCDM


text:
function jt08()
#
#  test polynomial suggested by Jenkins and Traub
#
    p = poly([-1.0;-1;-1;-1;-1])
    z = [-1 5]
    p, PolyZeros(z)
end


text:
using ExaTron
using Test
using Random
using LinearAlgebra
using SparseArrays
using StatsBase
using CUDA

@testset "Test ExaTron" begin
    if has_cuda_gpu()
        @testset "gputest" begin
            include("gputest.jl")
        end
    end
    @testset "qptests" begin
        include("qptest.jl")
    end
    @testset "densetest" begin
        include("densetest.jl")
    end
    @testset "proxal_wrapper" begin
        include("proxal_wrapper.jl")
    end
    @testset "admmtest" begin
        include("admmtest.jl")
    end
end


text:
import Sobol.next!

"""
    A `NumberGenerator` must always return a vector of uniforms when called with the function `next!(::NumberGenerator)`.
    It can contain whatever it wants as long as it does that method.
"""
abstract type NumberGenerator end

"""
    A NumberGenerator wrapper for `SobolSeq`. This makes quasirandom numbers
"""
struct SobolGen <: NumberGenerator
    seq_::SobolSeq
end

"""
    A NumberGenerator wrapper for `MersenneTwister`. This makes pseudorandom numbers.
"""
struct Mersenne <: NumberGenerator
    twister_::MersenneTwister
    number_of_itos_::UInt
end

"""
    next!(number_generator::NumberGenerator)
This extracts a random draw given a number generator struct (like a `SobolGen` or a `Mersenne`).

"""
function next!(number_generator::SobolGen)
    return next!(number_generator.seq_)
end
function next!(number_generator::Mersenne)
    return rand(number_generator.twister_, number_generator.number_of_itos_)
end


text:
using Pkg, Test, DistStat

type=[Float64,Float32]

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    using CUDA
    ArrayType = CuArray
else
    ArrayType = Array
end

for T in type

 data =ArrayType{T}(reshape(collect(1:42),6,7))
 data_dist1 = distribute(data)
 data_dist2 = distribute(ArrayType{T}(transpose(data)))
 cols1=data_dist1.partitioning[DistStat.Rank()+1][2]
 cols2=data_dist2.partitioning[DistStat.Rank()+1][2]

 @test data_dist1.localarray==data[:,cols1]
 @test data_dist2.localarray==(ArrayType{T}(transpose(data)))[:,cols2]

end


text:
using Test

using Oceananigans
using Oceananigans.Units
using Oceananigans.Simulations:
    stop_iteration_exceeded, stop_time_exceeded, wall_time_limit_exceeded,
    TimeStepWizard, new_time_step, reset!

using Dates: DateTime

include("utils_for_runtests.jl")

archs = test_architectures()

function wall_time_step_wizard_tests(arch)
    grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1))
    model = NonhydrostaticModel(grid=grid)

    Δx = grid.Δxᶜᵃᵃ
    CFL = 0.45
    u₀ = 7
    Δt = 2.5
    model.velocities.u[1, 1, 1] = u₀

    wizard = TimeStepWizard(cfl=CFL, max_change=Inf, min_change=0)
    Δt = new_time_step(Δt, wizard, model)
    @test Δt ≈ CFL * Δx / u₀

    wizard = TimeStepWizard(cfl=CFL, max_change=Inf, min_change=0.75)
    Δt = new_time_step(1.0, wizard, model)
    @test Δt ≈ 0.75

    wizard = TimeStepWizard(cfl=CFL, max_change=Inf, min_change=0, min_Δt=1.99)
    Δt = new_time_step(Δt, wizard, model)
    @test Δt ≈ 1.99

    model.velocities.u[1, 1, 1] = u₀/100

    wizard = TimeStepWizard(cfl=CFL, max_change=1.1, min_change=0)
    Δt = new_time_step(1.0, wizard, model)
    @test Δt ≈ 1.1

    wizard = TimeStepWizard(cfl=CFL, max_change=Inf, min_change=0, max_Δt=3.99)
    Δt = new_time_step(Δt, wizard, model)
    @test Δt ≈ 3.99

    grid_stretched = RectilinearGrid(arch,
                                    size = (1, 1, 1),
                                     x = (0, 1),
                                     y = (0, 1),
                                     z = z -> z,
                                     halo = (1, 1, 1))

    model = NonhydrostaticModel(grid=grid_stretched)

    Δx = grid_stretched.Δxᶜᵃᵃ
    CFL = 0.45
    u₀ = 7
    Δt = 2.5
    model.velocities.u .= u₀

    wizard = TimeStepWizard(cfl=CFL, max_change=Inf, min_change=0)
    Δt = new_time_step(Δt, wizard, model)
    @test Δt ≈ CFL * Δx / u₀

    return nothing
end

function run_basic_simulation_tests(arch)
    grid  = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1))
    model = NonhydrostaticModel(grid=grid)
    simulation = Simulation(model, Δt=3, stop_iteration=1)

    # Just make sure we can construct a simulation without any errors.
    @test simulation isa Simulation

    simulation.running = true
    stop_iteration_exceeded(simulation)
    @test simulation.running

    run!(simulation)

    # Just make sure run! executes without any errors.
    @test simulation isa Simulation

    # Some basic tests
    simulation.running = true
    stop_iteration_exceeded(simulation)
    @test !(simulation.running)

    @test model.clock.time ≈ simulation.Δt
    @test model.clock.iteration == 1
    @test simulation.run_wall_time > 0

    simulation.running = true
    stop_time_exceeded(simulation)
    @test simulation.running

    simulation.running = true
    simulation.stop_time = 1e-12 # less than the current time.
    stop_time_exceeded(simulation)
    @test !(simulation.running)

    simulation.running = true
    wall_time_limit_exceeded(simulation)
    @test simulation.running

    simulation.running = true
    simulation.wall_time_limit = 1e-12
    wall_time_limit_exceeded(simulation)
    @test !(simulation.running)

    # Test that simulation stops at `stop_iteration`.
    reset!(simulation)
    simulation.stop_iteration = 3
    run!(simulation)

    @test simulation.model.clock.iteration == 3

    # Test that simulation stops at `stop_time`.
    reset!(simulation)
    simulation.stop_time = 20.20
    run!(simulation)

    @test simulation.model.clock.time ≈ 20.20

    # Test that we can run a simulation with TimeStepWizard
    reset!(simulation)
    simulation.stop_iteration = 2

    wizard = TimeStepWizard(cfl=0.1)
    simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(1))

    run!(simulation)

    @test simulation.callbacks[:wizard].func isa TimeStepWizard

    # Test time-step alignment with callbacks
    reset!(simulation)
    simulation.stop_time = 2.0
    simulation.Δt = 1.0

    called_at = Float64[]
    schedule = TimeInterval(0.31)
    capture_call_time(sim) = push!(called_at, sim.model.clock.time)
    simulation.callbacks[:tester] = Callback(capture_call_time, schedule)
    run!(simulation)

    @show called_at
    @test all(called_at .≈ 0.0:schedule.interval:simulation.stop_time)

    return nothing
end

function run_simulation_date_tests(arch, start_time, stop_time, Δt)
    grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1))

    clock = Clock(time=start_time)
    model = NonhydrostaticModel(grid=grid, clock=clock)
    simulation = Simulation(model, Δt=Δt, stop_time=stop_time)

    @test model.clock.time == start_time
    @test simulation.stop_time == stop_time

    run!(simulation)

    @test model.clock.time == stop_time
    @test simulation.stop_time == stop_time

    return nothing
end

function run_nan_checker_test(arch; erroring)
    grid = RectilinearGrid(arch, size=(4, 2, 1), extent=(1, 1, 1))
    model = NonhydrostaticModel(grid=grid)
    simulation = Simulation(model, Δt=1, stop_iteration=1)
    model.velocities.u[1, 1, 1] = NaN
    erroring && erroring_NaNChecker!(simulation)

    if erroring
        @test_throws ErrorException run!(simulation)
    else
        run!(simulation)
        @test model.clock.iteration == 0 # simulation did not run
    end

    return nothing
end

@testset "Time step wizard" begin
    for arch in archs
        @info "Testing time step wizard [$(typeof(arch))]..."
        wall_time_step_wizard_tests(arch)
    end
end

@testset "Simulations" begin
    for arch in archs
        @info "Testing simulations [$(typeof(arch))]..."
        run_basic_simulation_tests(arch)

        @testset "NaN Checker [$(typeof(arch))]" begin
            @info "  Testing NaN Checker [$(typeof(arch))]..."
            run_nan_checker_test(arch, erroring=true)
            run_nan_checker_test(arch, erroring=false)
        end

        @info "Testing simulations with DateTime [$(typeof(arch))]..."
        run_simulation_date_tests(arch, 0.0, 1.0, 0.3)
        run_simulation_date_tests(arch, DateTime(2020), DateTime(2021), 100days)
        run_simulation_date_tests(arch, TimeDate(2020), TimeDate(2021), 100days)
    end
end


text:
function update_feed!(
    state::Tuple{AbstractGraph, AbstractArray}, agent_idx::Integer,
    config::Config
)
    graph, agent_list = state
    this_agent = agent_list[agent_idx]
    unique!(this_agent.feed)
    deleted_posts = Integer[]
    for (index, post) in enumerate(this_agent.feed)
        if post.weight == -1 || !(post.source_agent in inneighbors(graph, agent_idx))
            push!(deleted_posts, index)
        else
            post.weight = 0.5 * post.weight
        end
    end
    deleteat!(this_agent.feed, deleted_posts)
    sort!(this_agent.feed, lt=<, rev=true)
    if length(this_agent.feed) > 15
        this_agent.feed = this_agent.feed[1:15]
    end

    for post in this_agent.feed
        if !(agent_idx in post.seen_by)
            push!(post.seen_by, agent_idx)
        end
    end

    return state
end

# suppress output of include()
;


text:
using CartesianGeneticProgramming
using Cambrian
using ArcadeLearningEnvironment
using ArgParse
import Cambrian.mutate
import Random

```
Playing Atari games using classic CGP on RAM values

If an individual is provided using --ind, an evaluation loop with rendering will
be performed. Otherwise, an evolution is launched. By default this uses seed=0
for each evaluation for a deterministic environment, but this can be removed for
a stochastic and more realistic result.
```

s = ArgParseSettings()
@add_arg_table s begin
    "--cfg"
    help = "configuration script"
    default = "cfg/atari_ram.yaml"
    "--game"
    help = "game rom name"
    default = "centipede"
    "--seed"
    help = "random seed for evolution"
    arg_type = Int
    default = 0
    "--ind"
    help = "individual for evaluation"
    arg_type = String
    default = ""
end
args = parse_args(ARGS, s)

function play_atari(ind::CGPInd, rom_name::String; seed=0, max_frames=18000, render=false)
    ale = ALE_new()
    setInt(ale, "random_seed", Cint(seed))
    if render
        setBool(ale, "display_screen", true)
        setBool(ale, "sound", true)
    end
    loadROM(ale, rom_name)
    actions = getMinimalActionSet(ale)
    reward = 0.0
    frames = 0
    while ~game_over(ale)
        ram = getRAM(ale) ./ typemax(UInt8)
        output = process(ind, ram)
        action = actions[argmax(output)]
        reward += act(ale, action)
        frames += 1
        if frames > max_frames
            break
        end
    end
    ALE_del(ale)
    [reward]
end

ale = ALE_new()
loadROM(ale, args["game"])
n_in = length(getRAM(ale))
n_out = length(getMinimalActionSet(ale))
ALE_del(ale)

cfg = get_config(args["cfg"]; game=args["game"], n_in=n_in, n_out=n_out)
Random.seed!(args["seed"])

if length(args["ind"]) > 0
    ind = CGPInd(cfg, read(args["ind"], String))
    reward = play_atari(ind, args["game"]; seed=args["seed"], render=true)
    println(reward)
else
    mutate(i::CGPInd) = goldman_mutate(cfg, i)
    fit(i::CGPInd) = play_atari(i, cfg.game)
    e = CGPEvolution(cfg, fit)
    run!(e)
end


text:
using Interpolations;

function showGif(field,title_in="Solar Forcing of the Earth: t = ",minval=-46.073107944313556,maxval=32.52216883281448)
    World = read_geography(joinpath(@__DIR__,"..", "input", "The_World_Outline.dat"))*(-80)
    (width, height) = size(World) # NO hinzugefügt
    Plots.@gif for t in 1:size(field,3)
        title = title_in * lpad(string(floor(t*365/48)), 3, ' ') * " days"

        data = field[:,:,t]
        data = transpose(data)
        Plots.heatmap(LinRange(-180,180,width),LinRange(-90,90,height),data,
            clims=(minval,maxval),
            yflip=true,
            title=title,
            xlabel="longitude [°]",ylabel="latitude [°]")
        Plots.contour!(LinRange(-180,180,width),LinRange(-90,90,height),transpose(World),xlim=(-180,180), ylim=(-90, 90),levels=-40:40,seriesalpha=0.05)
        fpath = "png/field/" * lpad(string(t), 3, '0') * ".png"
        println(fpath)
        Plots.savefig(fpath)
    end
end

# This is a dirty function. What was it again?
# First make it work, then make it pretty ... something like that.(?)
function apply_richardson_projection(data)

    lat = [0 , 5 , 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90] * pi/180.0
    X   = [1.0000, 0.9986, 0.9954, 0.99  , 0.9822, 0.973 , 0.96  , 0.9427, 0.9216, 0.8962, 0.8679, 0.835 , 0.7986, 0.7597, 0.7186, 0.6732, 0.6213, 0.5722, 0.5322]
    Y   = [0.0000, 0.0620, 0.1240, 0.1860, 0.2480, 0.3100, 0.3720, 0.4340, 0.4958, 0.5571, 0.6176, 0.6769, 0.7346, 0.7903, 0.8435, 0.8936, 0.9394, 0.9761, 1.0000]

    itpl_X = LinearInterpolation(lat,X);
    itpl_Y = LinearInterpolation(lat,Y);

    calc_x = (radius,longitude,latitude) -> 0.8487 * radius .* itpl_X.(abs.(latitude)) .* longitude;
    calc_y = (radius,longitude,latitude) -> 1.3523 * radius .* itpl_Y.(abs.(latitude)) .* 0.2536 .* sign.(latitude)

    nx = size(data,1)
    ny = size(data,2)

    longitude = LinRange(-pi,pi,nx) .* ones(ny)'
    latitude  = ones(nx) .* LinRange(-pi/2,pi/2,ny)'

    x = calc_x(1.0,longitude,latitude);
    y = calc_y(1.0,longitude,latitude);
    z = reverse(reverse(data),dims=1)

    x = 180.0*(x ./ maximum(x));
    y =  90.0*(y ./ maximum(y));

    return x,y,z

end


text:
import Base.(==)
#=
== : get a new Assoc where all of the elements of input Assoc mataches the given Element.
=#
function ==(A::Assoc, E::Union{AbstractString,Number})
    tarIndex = searchsortedfirst(A.val,E)
    if (isa(E,Number) & (size(Val(A),1)==1) & (Val(A)[1] == 1.0)  )
        tarIndex = E
    end
    rowkey, colkey, valkey = findnz(A.A)
    mapping = find( x-> x == tarIndex, valkey)
#    rowkey , colkey = unique(rowkey[mapping]), unique(colkey[mapping])
#    sort!(rowkey)
#    sort!(colkey)
    return A[rowkey[mapping],colkey[mapping]]
end

==(E::Union{AbstractString,Number},A::Assoc) = (A == E)

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
module Hydrogen_FEM
    include("hydrogen_fem_module.jl")
    using LinearAlgebra
    using Match
    using MKL
    using .Hydrogen_FEM_module

    const NODE_TOTAL = 10000

    function construct()
        param = Hydrogen_FEM_module.Hydrogen_FEM_param("eigenfunc.csv", "eigenval.csv", NODE_TOTAL, NODE_TOTAL - 1, 1000.0, 0.0)
        val = Hydrogen_FEM_module.Hydrogen_FEM_variables(
            Symmetric(zeros(param.ELE_TOTAL, param.ELE_TOTAL)),
            Array{Float64}(undef, param.ELE_TOTAL),
            Array{Float64, 3}(undef, param.ELE_TOTAL, 2, 2),
            Array{Float64, 3}(undef, param.ELE_TOTAL, 2, 2),
            Array{Int64, 2}(undef, param.ELE_TOTAL, 2),
            Array{Float64, 2}(undef, param.ELE_TOTAL, 2),
            Array{Float64}(undef, param.NODE_TOTAL),
            Array{Float64}(undef, param.NODE_TOTAL),
            Symmetric(zeros(param.ELE_TOTAL, param.ELE_TOTAL)))

        return param, val
    end

    function do_run(param, val)
        # データの生成
        make_data!(param, val)

        # 要素行列の生成
        make_element_matrix!(param, val)

        # 全体行列を生成
        hg_tmp, ug_tmp = make_global_matrix(param, val)

        # 境界条件処理を行う
        boundary_conditions!(param, val, hg_tmp, ug_tmp)

        # 一般化固有値問題を解く
        eigenval, phi = eigen!(val.hg, val.ug)

        # 基底状態の固有ベクトルを取り出す
        val.phi = @view(phi[:,1])

        # 固有ベクトルの要素数を増やす
        resize!(val.phi, NODE_TOTAL)

        # 端点rcでの値を0にする
        val.phi[NODE_TOTAL] = 0.0

        # 固有ベクトル（波動関数）を規格化
        normalize!(val)

        return eigenval
    end

    save_result(param, eigenval, val) = let
        open(param.EIGENFUNC_FILENAME, "w" ) do fp
            for i = 1:length(val.phi)
                println(fp, @sprintf "%.14f, %.14f, %.14f" val.node_r_glo[i] val.phi[i] 2.0 * exp(- val.node_r_glo[i]))
            end
        end

        open(param.EIGENVAL_FILENAME, "w" ) do fp
            for i = 1:length(eigenval)
                println(fp, @sprintf "%d, %.14f, %.14f" i eigenval[i] -0.5 * 1.0 / float(i * i))
            end
        end
    end

    function boundary_conditions!(param, val, hg_tmp, ug_tmp)
        @inbounds for i = 1:param.ELE_TOTAL
            for j = i - 1:i + 1
                if j != 0 && j != param.NODE_TOTAL
                    # 左辺の全体行列のN行とN列を削る
                    val.hg.data[j, i] = hg_tmp.data[j, i]

                    # 右辺の全体行列のN行とN列を削る
                    val.ug.data[j, i] = ug_tmp.data[j, i]
                end
            end
        end
    end

    get_A_matrix_element(e, le, p, q) = let
        ed = float(e - 1)
        @match p begin
            1 =>
                @match q begin
                    1 => return  0.5 * le * (ed * ed + ed + 1.0 / 3.0) - le * le * (ed / 3.0 + 1.0 / 12.0)
                    2 => return -0.5 * le * (ed * ed + ed + 1.0 / 3.0) - le * le * (ed / 6.0 + 1.0 / 12.0)
                    _ => return 0.0
                end

            2 =>
                @match q begin
                    2 => return 0.5 * le * (ed * ed + ed + 1.0 / 3.0) - le * le * (ed / 3.0 + 0.25)
                    _ => return 0.0
                end

            _ => return 0.0
        end
    end

    get_B_matrix_element(e, le, p, q) = let
        ed = float(e - 1)
        @match p begin
            1 =>
                @match q begin
                    1 => return le * le * le * (ed * ed / 3.0 + ed / 6.0 + 1.0 / 30.0)
                    2 => return le * le * le * (ed * ed / 6.0 + ed / 6.0 + 0.05)
                    _ => return 0.0
                end

            2 =>
                @match q begin
                    2 => return le * le * le * (ed * ed / 3.0 + ed / 2.0 + 0.2)
                    _ => return 0.0
                end

            _ => return 0.0
        end
    end

    function make_data!(param, val)
        # Global節点のx座標を定義(R_MIN～R_MAX）
        dr = (param.R_MAX - param.R_MIN) / float(param.ELE_TOTAL)
        @inbounds for i = 0:param.NODE_TOTAL - 1
            # 計算領域を等分割
            val.node_r_glo[i + 1] = param.R_MIN + float(i) * dr
        end

        @inbounds for e = 1:param.ELE_TOTAL
            val.node_num_seg[e, 1] = e
            val.node_num_seg[e, 2] = e + 1
        end

        @inbounds for e = 1:param.ELE_TOTAL
            for i = 1:2
                val.node_r_ele[e, i] = val.node_r_glo[val.node_num_seg[e, i]]
            end
        end
    end

    function make_element_matrix!(param, val)
        # 各線分要素の長さを計算
        @inbounds for e = 1:param.ELE_TOTAL
            val.length[e] = abs(val.node_r_ele[e, 2] - val.node_r_ele[e, 1])
        end

        # 要素行列の各成分を計算
        @inbounds for e = 1:param.ELE_TOTAL
            le = val.length[e]
            for j = 1:2
                for i = 1:j
                    val.mat_A_ele[e, i, j] = get_A_matrix_element(e, le, i, j)
                    val.mat_B_ele[e, i, j] = get_B_matrix_element(e, le, i, j)
                    end
                end
            end
        end

    function make_global_matrix(param, val)
        hg_tmp = Symmetric(zeros(param.NODE_TOTAL, param.NODE_TOTAL))
        ug_tmp = Symmetric(zeros(param.NODE_TOTAL, param.NODE_TOTAL))

        @inbounds for e = 1:param.ELE_TOTAL
            for j = 1:2
                for i = 1:j
                    hg_tmp.data[val.node_num_seg[e, i], val.node_num_seg[e, j]] += val.mat_A_ele[e, i, j]
                    ug_tmp.data[val.node_num_seg[e, i], val.node_num_seg[e, j]] += val.mat_B_ele[e, i, j]
                end
            end
        end

        return hg_tmp, ug_tmp
    end

    function normalize!(val)
        sum = 0.0
        max = length(val.phi) - 2

        # Simpsonの公式によって数値積分する
        @inbounds @simd for i = 1:2:max
            f0 = val.phi[i] * val.phi[i] * val.node_r_glo[i] * val.node_r_glo[i]
            f1 = val.phi[i + 1] * val.phi[i + 1] * val.node_r_glo[i + 1] * val.node_r_glo[i + 1]
            f2 = val.phi[i + 2] * val.phi[i + 2] * val.node_r_glo[i + 2] * val.node_r_glo[i + 2]

            sum += (f0 + 4.0 * f1 + f2)
        end

        val.phi = map(x -> abs(x / sqrt(sum * val.length[1] / 3.0)), val.phi)
    end
end


text:
hsv = [rand(1:360); rand(2)...]
rgb = hsvtorgb(hsv)

@test length(rgb) == 3

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
"""
Perdew, Burke, Ernzerhof (DOI: 10.1103/PhysRevLett.77.3865)
"""
function energy_per_particle(::Val{:gga_c_pbe}, ρ, σ)
    β = 0.06672455060314922
    γ = (1 - log(2)) / π^2
    energy_per_particle_c_pbe(ρ, σ; β, γ)
end

"""
Xu, Goddard 2004 (DOI 10.1063/1.1771632)
"""
function energy_per_particle(::Val{:gga_c_xpbe}, ρ, σ)
    β = 0.089809  # Fitted constants, Table I
    α = 0.197363  # Fitted constants, Table I
    γ = β^2 / 2α
    energy_per_particle_c_pbe(ρ, σ; β, γ)
end

"""
Perdew, Ruzsinszky, Csonka and others 2008 (DOI 10.1103/physrevlett.100.136406)
"""
function energy_per_particle(::Val{:gga_c_pbe_sol}, ρ, σ)
    β = 0.046  # Page 3, left column below figure 1
    γ = (1 - log(2)) / π^2
    energy_per_particle_c_pbe(ρ, σ; β, γ)
end

"""
Constantin, Fabiano, Laricchia 2011 (DOI 10.1103/physrevlett.106.186406)
"""
function energy_per_particle(::Val{:gga_c_apbe}, ρ, σ)
    μ = 0.260   # p. 1, right column, bottom
    β = 3μ / π^2
    γ = (1 - log(2)) / π^2  # like in PBE
    energy_per_particle_c_pbe(ρ, σ; β, γ)
end

"""
del Campo, Gazqez, Trickey and others 2012 (DOI 10.1063/1.3691197)
"""
function energy_per_particle(::Val{:gga_c_pbe_mol}, ρ, σ)
    # β made to cancel self-interaction error in hydrogen
    β = 0.08384             # p. 4, right column, first paragraph
    γ = (1 - log(2)) / π^2  # like in PBE
    energy_per_particle_c_pbe(ρ, σ; β, γ)
end

"""
Sarmiento-Perez, Silvana, Marques 2015 (DOI 10.1021/acs.jctc.5b00529)
"""
function energy_per_particle(::Val{:gga_c_pbefe}, ρ, σ)
    β = 0.043                    # Fitted constants, Table I
    γ = 0.031090690869654895034  # Fitted constants, Table I
    energy_per_particle_c_pbe(ρ, σ; β, γ)
end

function energy_per_particle_c_pbe(ρ, σ; β, γ)
    # Spin-scaling factor with ζ spin polarization.
    # Yue Wang and John P. Perdew. Phys. Rev. B 43, 8911 (1991).
    # DOI 10.1103/PhysRevB.43.8911
    ϕ(ζ) = ((1+ζ)^(2/3) + (1-ζ)^(2/3))/2  # == 1 for non-spin-polarised

    # ε = UEG correlation energy per particle
    A(ε, ϕ³) = β/γ / expm1(-ε / (γ * ϕ³))  # (8)
    function H(ε, t², ϕ³)  # (7)
        At² = A(ε, ϕ³) * t²
        γ * ϕ³ * log(1 + β/γ * t² * (1 + At²) / (1 + At² + (At²)^2))
    end

    phi = #= ϕ(ζ) =# 1.0
    ε = energy_per_particle(Val(:lda_c_pw), ρ; improved_constants=true)
    t² = (1/12 * 3^(5/6) * π^(1/6))^2 * σ / (phi^2 * ρ^(7/3))  # page 2, left column, top
    ε + H(ε, t², phi)
end


text:
using RobotDynamics
using StaticArrays
using Test
using BenchmarkTools
using ForwardDiff
using FiniteDiff
using LinearAlgebra
using Random
const RD = RobotDynamics

using RobotDynamics: @autodiff, state, control, KnotPoint, getdata, getstate, getcontrol, getparams,
                     FiniteDifference, ForwardAD, StaticReturn, InPlace

using RobotDynamics: evaluate, evaluate!, jacobian!
function test_allocs(fun)
    n,m,p = RD.dims(fun)
    z_ = @SVector randn(n+m)
    z = KnotPoint{n,m}(z_,1.0,0.1)
    test_allocs(fun, z)
end
function test_allocs(fun, z)
    n = RD.state_dim(z)
    m = RD.control_dim(z)
    z_ = RD.getdata(z)
    p = RD.output_dim(fun)
    w = RD.input_dim(fun)
    y = zeros(p)
    J = zeros(p,w)
    allocs = 0
    allocs += @allocated evaluate(fun, z)
    allocs += @allocated evaluate!(fun, y, z)
    allocs += @allocated jacobian!(RobotDynamics.InPlace(), RobotDynamics.ForwardAD(), fun, J, y, z)
    allocs += @allocated jacobian!(RobotDynamics.InPlace(), RobotDynamics.FiniteDifference(), fun, J, y, z)
    allocs += @allocated jacobian!(RobotDynamics.StaticReturn(), RobotDynamics.ForwardAD(), fun, J, y, z)
    allocs += @allocated jacobian!(RobotDynamics.StaticReturn(), RobotDynamics.FiniteDifference(), fun, J, y, z)

    # check for zero allocations with normal array inputs for in place methods
    z = KnotPoint{n,m}(Vector(z_),1.0,0.1)
    allocs += @allocated evaluate(fun, z)
    allocs += @allocated evaluate!(fun, y, z)
    allocs += @allocated jacobian!(RobotDynamics.InPlace(), RobotDynamics.ForwardAD(), fun, J, y, z)
    allocs += @allocated jacobian!(RobotDynamics.InPlace(), RobotDynamics.FiniteDifference(), fun, J, y, z)
    allocs += @allocated jacobian!(RobotDynamics.StaticReturn(), RobotDynamics.FiniteDifference(), fun, J, y, z)
end

function test_fun(fun)
    Random.seed!(1)
    @test RD.dims(fun) == (2,2,3)

    z_ = @SVector randn(4)
    z = KnotPoint{2,2}(z_,1.0,0.1)
    y = zeros(3)
    J = zeros(3,4)
    J0 = zeros(3,4)

    evaluate!(fun, y, z)
    @test y ≈ evaluate(fun, z)

    jacobian!(fun, J0, y, z)

    jacobian!(RobotDynamics.InPlace(), RobotDynamics.UserDefined(), fun, J, y, z)
    @test J ≈ J0
    jacobian!(RobotDynamics.StaticReturn(), RobotDynamics.UserDefined(), fun, J, y, z)
    @test J ≈ J0

    jacobian!(RobotDynamics.InPlace(), RobotDynamics.ForwardAD(), fun, J, y, z)
    @test J ≈ J0
    jacobian!(RobotDynamics.StaticReturn(), RobotDynamics.ForwardAD(), fun, J, y, z)
    @test J ≈ J0

    jacobian!(RobotDynamics.InPlace(), RobotDynamics.FiniteDifference(), fun, J, y, z)
    @test J ≈ J0 atol = 1e-6
    jacobian!(RobotDynamics.StaticReturn(), RobotDynamics.FiniteDifference(), fun, J, y, z)
    @test J ≈ J0 atol = 1e-6

    run_alloc_tests && (@test test_allocs(fun) == 0)
end

##
struct TestFun0{CH} <: RobotDynamics.AbstractFunction
    cfg::ForwardDiff.JacobianConfig{Nothing, Float64, CH, Tuple{Vector{ForwardDiff.Dual{Nothing, Float64, CH}}, Vector{ForwardDiff.Dual{Nothing, Float64, CH}}}}
    cache::FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}, Nothing, Val{:forward}(), Float64}
    function TestFun0()
        n,m,p = 2,2,3
        cfg = ForwardDiff.JacobianConfig(nothing, zeros(p), zeros(n+m))
        cache = FiniteDiff.JacobianCache(zeros(n+m), zeros(p))
        new{length(cfg.seeds)}(cfg, cache)
    end
end
function RobotDynamics.evaluate(::TestFun0, x, u, p)
    return SA[cos(x[1]) * u[1], sin(x[2]^2 * x[1]) * u[2], exp(x[2] + x[1]/10)] * p[1]
end
function RobotDynamics.evaluate!(::TestFun0, y, x, u, p)
    y[1] = cos(x[1])  * u[1]
    y[2] = sin(x[2]^2 * x[1]) * u[2]
    y[3] = exp(x[2] + x[1]/10)
    y .*= p[1]
    return nothing
end
RobotDynamics.state_dim(::TestFun0) = 2
RobotDynamics.control_dim(::TestFun0) = 2
RobotDynamics.output_dim(::TestFun0) = 3

function RobotDynamics.jacobian!(::TestFun0, J, y, x, u, p)
    J .= 0
    J[1,1] = -sin(x[1]) * u[1]
    J[1,3] = cos(x[1])
    J[2,1] = x[2]^2 * cos(x[2]^2 * x[1]) * u[2]
    J[2,2] = 2 * x[1] * x[2] * cos(x[2]^2 * x[1]) * u[2]
    J[2,4] = sin(x[2]^2 * x[1])
    J[3,1] = exp(x[2] + x[1]/10) / 10
    J[3,2] = exp(x[2] + x[1]/10)
    J .*= p[1]
    return nothing
end

function RobotDynamics.jacobian!(::StaticReturn, ::ForwardAD, fun::TestFun0, J, y, z)
    f(_z) = RobotDynamics.evaluate(fun, getstate(z, _z), getcontrol(z, _z), getparams(z))
    J .= ForwardDiff.jacobian(f, getdata(z))
    return nothing
end

function RobotDynamics.jacobian!(::InPlace, ::ForwardAD, fun::TestFun0, J, y, z)
    f!(_y,_z) = RobotDynamics.evaluate!(fun, _y, getstate(z, _z), getcontrol(z, _z), getparams(z))
    ForwardDiff.jacobian!(J, f!, y, getdata(z), fun.cfg)
    return nothing
end

function RobotDynamics.jacobian!(::StaticReturn, ::FiniteDifference, fun::TestFun0, J, y, z)
    f!(_y,_z) = _y .= RobotDynamics.evaluate(fun, getstate(z, _z), getcontrol(z, _z), getparams(z))
    FiniteDiff.finite_difference_jacobian!(J, f!, getdata(z), fun.cache)
    return nothing
end

function RobotDynamics.jacobian!(::InPlace, ::FiniteDifference, fun::TestFun0, J, y, z)
    f!(_y,_z) = RobotDynamics.evaluate!(fun, _y, getstate(z, _z), getcontrol(z, _z), getparams(z))
    FiniteDiff.finite_difference_jacobian!(J, f!, getdata(z), fun.cache)
    return nothing
end

fun = TestFun0()
n,m,p = RD.dims(fun)
x = @SVector randn(n)
u = @SVector randn(m)
t = 1.2
dt = 0.1

zs = KnotPoint{n,m}([x;u],t,dt)
z = KnotPoint{n,m}(Vector([x;u]),t,dt)
z_ = copy(z.z)
@test getstate(zs, z_) isa SVector{n}
@test getcontrol(zs, z_) isa SVector{m}
@test getstate(z, z_) isa SubArray
@test getcontrol(z, z_) isa SubArray

RD.evaluate(fun, zs)
y = zeros(3)
RD.evaluate!(fun, y, z)

test_fun(fun)

##############################
# Autogen With Inner Constructor
##############################
@autodiff struct TestFun <: RobotDynamics.AbstractFunction
    a::Int
    function TestFun()
        new(1.0)
    end
end

RobotDynamics.state_dim(::TestFun) = 2
RobotDynamics.control_dim(::TestFun) = 2
RobotDynamics.output_dim(::TestFun) = 3

function RobotDynamics.evaluate(::TestFun, x, u)
    return SA[sin(x[1]), cos(x[2]), u[1] * exp(u[2])]
end

function RobotDynamics.evaluate!(::TestFun, y, x, u)
    y[1] = sin(x[1])
    y[2] = cos(x[2])
    y[3] = u[1] * exp(u[2])
end

function RobotDynamics.jacobian!(::TestFun, J, y, x, u)
    J .= 0
    J[1,1] = cos(x[1])
    J[2,2] = -sin(x[2])
    J[3,3] = exp(u[2])
    J[3,4] = u[1] * exp(u[2])
    return nothing
end

fun = TestFun()
@test :cfg ∈ fieldnames(typeof(fun))
@test :cache ∈ fieldnames(typeof(fun))
@test fun.a == 1.0
test_fun(fun)

##############################
# With Input Parameters
##############################
@autodiff struct TestFunTime <: RobotDynamics.AbstractFunction end

function RobotDynamics.evaluate(::TestFunTime, x, u, t)
    return SA[sin(x[1]) * t.t, cos(x[2]) * t.t, u[1] * exp(u[2])]
end

function RobotDynamics.evaluate!(::TestFunTime, y, x, u, t)
    y[1] = sin(x[1]) * t.t
    y[2] = cos(x[2]) * t.t
    y[3] = u[1] * exp(u[2])
end

function RobotDynamics.jacobian!(::TestFunTime, J, y, x, u, t)
    J .= 0
    J[1,1] = cos(x[1]) * t.t
    J[2,2] = -sin(x[2]) * t.t
    J[3,3] = exp(u[2])
    J[3,4] = u[1] * exp(u[2])
    return nothing
end
RobotDynamics.state_dim(::TestFunTime) = 2
RobotDynamics.control_dim(::TestFunTime) = 2
RobotDynamics.output_dim(::TestFunTime) = 3

fun = TestFunTime()
z_ = @SVector zeros(4)
z = KnotPoint{2,2}(z_, 1.2, 0.1)
y = zeros(3)
J = zeros(3,4)
out = RobotDynamics.evaluate(fun, z)
RobotDynamics.evaluate!(fun, y, z)
@test y ≈ out
@test out[2] == z.t

jacobian!(RobotDynamics.StaticReturn(), RobotDynamics.ForwardAD(), fun, J, y, z)
@test J[1,1] == z.t
z.t = 2.4
jacobian!(RobotDynamics.StaticReturn(), RobotDynamics.ForwardAD(), fun, J, y, z)
@test J[1,1] == z.t

test_fun(fun)

##############################
# With Type Parameters
##############################
@autodiff struct TestFunParam{T} <: RobotDynamics.AbstractFunction
    a::T
    function TestFunParam(a::T) where T
        new{T}(a)
    end
end

RobotDynamics.state_dim(::TestFunParam) = 2
RobotDynamics.control_dim(::TestFunParam) = 2
RobotDynamics.output_dim(::TestFunParam) = 3

function RobotDynamics.evaluate(::TestFunParam, x, u)
    return SA[sin(x[1]), cos(x[2]), u[1] * exp(u[2])]
end

function RobotDynamics.evaluate!(::TestFunParam, y, x, u)
    y[1] = sin(x[1])
    y[2] = cos(x[2])
    y[3] = u[1] * exp(u[2])
end

function RobotDynamics.jacobian!(::TestFunParam, J, y, x, u)
    J .= 0
    J[1,1] = cos(x[1])
    J[2,2] = -sin(x[2])
    J[3,3] = exp(u[2])
    J[3,4] = u[1] * exp(u[2])
    return nothing
end

fun = TestFunParam(1.0)
@test fun.a == 1.0
test_fun(fun)

##############################
# Without inner constructor and inherited
##############################
abstract type TestFunBase <: RobotDynamics.AbstractFunction end

@autodiff struct TestFunInherited <: TestFunBase end
RobotDynamics.state_dim(::TestFunBase) = 2
RobotDynamics.control_dim(::TestFunBase) = 2
RobotDynamics.output_dim(::TestFunBase) = 3

function RobotDynamics.evaluate(::TestFunBase, x, u)
    return SA[sin(x[1]), cos(x[2]), u[1] * exp(u[2])]
end

function RobotDynamics.evaluate!(::TestFunBase, y, x, u)
    y[1] = sin(x[1])
    y[2] = cos(x[2])
    y[3] = u[1] * exp(u[2])
end

function RobotDynamics.jacobian!(::TestFunBase, J, y, x, u)
    J .= 0
    J[1,1] = cos(x[1])
    J[2,2] = -sin(x[2])
    J[3,3] = exp(u[2])
    J[3,4] = u[1] * exp(u[2])
    return nothing
end

fun = TestFunInherited()
test_fun(fun)

##############################
# One-liner inner constructor
##############################
@autodiff struct TestInner <: TestFunBase
    a::Int
    TestInner(a) = new(a)
end

fun = TestInner(1)
test_fun(fun)

##############################
## State-Only
##############################
@autodiff struct TestFunStateOnly <: RobotDynamics.AbstractFunction end

RD.state_dim(::TestFunStateOnly) = 3
RD.output_dim(::TestFunStateOnly) = 2
RD.functioninputs(::TestFunStateOnly) = RD.StateOnly()

function RD.evaluate(::TestFunStateOnly, x::RD.DataVector)
    return SA[cos(x[1] - x[3]), sin(x[2]^2 * x[1])]
end

function RD.evaluate!(::TestFunStateOnly, y, x::RD.DataVector)
    y[1] = cos(x[1] - x[3])
    y[2] = sin(x[2]^2 * x[1])
    return nothing
end

function RD.jacobian!(::TestFunStateOnly, J, y, x::RD.DataVector)
    J .= 0
    J[1,1] = -sin(x[1] - x[3])
    J[1,3] = sin(x[1] - x[3])
    J[2,1] = cos(x[2]^2 * x[1]) * x[2]^2
    J[2,2] = cos(x[2]^2 * x[1]) * 2 * x[2] * x[1]
    return nothing
end

fun = TestFunStateOnly()
@test RD.functioninputs(fun) == RD.StateOnly()

n = RD.state_dim(fun)
p = RD.output_dim(fun)
x = @SVector randn(n)
u = @SVector randn(p)
y = zeros(p)
J = zeros(p,n)
J0 = copy(J)
fun = TestFunStateOnly()
zs = RD.KnotPoint(x,u,0.0,0.1)
z = RD.KnotPoint(Vector(x),Vector(u),0.0,0.1)

evaluate!(fun, y, z)
@test y ≈ evaluate(fun, z) ≈ evaluate(fun, x)
evaluate!(fun, y, x)
@test y ≈ evaluate(fun, z) ≈ evaluate(fun, x)

jacobian!(fun, J0, y, x)
jacobian!(RD.InPlace(), RD.UserDefined(), fun, J, y, z)
@test J ≈ J0
jacobian!(RD.StaticReturn(), RD.UserDefined(), fun, J, y, z)
@test J ≈ J0

jacobian!(RD.StaticReturn(), RD.ForwardAD(), fun, J, y, z)
@test J ≈ J0
jacobian!(RD.InPlace(), RD.ForwardAD(), fun, J, y, z)
@test J ≈ J0
jacobian!(RD.StaticReturn(), RD.FiniteDifference(), fun, J, y, z)
@test J ≈ J0 atol=1e-5
jacobian!(RD.InPlace(), RD.FiniteDifference(), fun, J, y, z)
@test J ≈ J0 atol=1e-5

@test RD.getinput(RD.functioninputs(fun), z) === state(z)
run_alloc_tests && @test test_allocs(fun, zs) == 0

##############################
## Control Only
##############################
@autodiff struct TestFunControlOnly <: RobotDynamics.AbstractFunction end

RD.control_dim(::TestFunControlOnly) = 2
RD.output_dim(::TestFunControlOnly) = 3
RD.functioninputs(::TestFunControlOnly) = RD.ControlOnly()

function RD.evaluate(::TestFunControlOnly, u::RD.DataVector)
    SA[exp(10*u[1]/u[2]), 1/u[1], (u[1] - u[2])^2]
end

function RD.evaluate!(::TestFunControlOnly, y, u::RD.DataVector)
    y[1] = exp(10*u[1]/u[2])
    y[2] = 1/u[1]
    y[3] = (u[1] - u[2])^2
    return nothing
end

function RD.jacobian!(::TestFunControlOnly, J, y, u::RD.DataVector)
    J .= 0
    e = exp(10*u[1] / u[2])
    J[1,1] = 10/u[2] * e
    J[1,2] = -10*u[1]/(u[2]^2) * e
    J[2,1] = -1/u[1]^2
    J[3,1] = 2*(u[1] - u[2])
    J[3,2] = -2*(u[1] - u[2])
    return nothing
end

fun = TestFunControlOnly()

m = RD.control_dim(fun)
p = RD.output_dim(fun)
x = @SVector randn(4)
u = @SVector randn(m)
y = zeros(p)
J = zeros(p,m)
J0 = copy(J)
zs = RD.KnotPoint(x,u,0.0,0.1)
z = RD.KnotPoint(Vector(x),Vector(u),0.0,0.1)

evaluate!(fun, y, z)
@test y ≈ evaluate(fun, z) ≈ evaluate(fun, u)
evaluate!(fun, y, u)
@test y ≈ evaluate(fun, z) ≈ evaluate(fun, u)

jacobian!(fun, J0, y, u)
jacobian!(RD.InPlace(), RD.UserDefined(), fun, J, y, z)
@test J ≈ J0
J
jacobian!(RD.StaticReturn(), RD.UserDefined(), fun, J, y, z)
@test J ≈ J0

jacobian!(RD.StaticReturn(), RD.ForwardAD(), fun, J, y, z)
@test J ≈ J0
jacobian!(RD.InPlace(), RD.ForwardAD(), fun, J, y, z)
@test J ≈ J0
jacobian!(RD.StaticReturn(), RD.FiniteDifference(), fun, J, y, z)
@test J ≈ J0 rtol=1e-2
jacobian!(RD.InPlace(), RD.FiniteDifference(), fun, J, y, z)
@test J ≈ J0 rtol=1e-2

@test RD.getinput(RD.functioninputs(fun), z) === control(z)
run_alloc_tests && @test test_allocs(fun, zs) == 0


text:
#!/usr/bin/env julia

using DocOpt

import Celeste.ParallelRun: BoundingBox, infer_box, infer_field
import Celeste.Stripe82Score: score_field_disk, score_object_disk
import Celeste.SDSSIO: RunCamcolField

const DOC =
"""Run Celeste.

Usage:
  celeste.jl [options] infer-box <ramin> <ramax> <decmin> <decmax> <outdir>
  celeste.jl [options] infer-field <run> <camcol> <field> <outdir>
  celeste.jl [options] infer-object <run> <camcol> <field> <objid> <outdir>
  celeste.jl [options] score-field <run> <camcol> <field> <resultdir> <truthfile>
  celeste.jl [options] score-object <run> <camcol> <field> <objid> <resultdir> <truthfile>
  celeste.jl -h | --help
  celeste.jl --version

Options:
  -h, --help         Show this screen.
  --version          Show the version.
  --logging=<level>  Level for the Lumberjack package (debug, info, warn, error). [default: info]

The `infer-box` subcommand runs Celeste on all sources in the given
RA/Dec range, using all available images.

The `infer-field` subcommand runs Celeste on all sources in the given
field, using only that field.

The `score-field` subcommand is not yet implemented for the new API.
"""

function main()
    args = docopt(DOC, version=v"0.1.0", options_first=true)
#   TODO: re-enable selective logging by level
#    set_logging_level(args["--logging"])
    stagedir = ENV["CELESTE_STAGE_DIR"]
    if args["infer-box"]
        ramin = parse(Float64, args["<ramin>"])
        ramax = parse(Float64, args["<ramax>"])
        decmin = parse(Float64, args["<decmin>"])
        decmax = parse(Float64, args["<decmax>"])
        box = BoundingBox(ramin, ramax, decmin, decmax)
        outdir = args["<outdir>"]
        infer_box(box, stagedir, outdir)
    else
        run = parse(Int, args["<run>"])
        camcol = parse(Int, args["<camcol>"])
        field = parse(Int, args["<field>"])
        rcf = RunCamcolField(run, camcol, field)
        if args["infer-field"]
            outdir = args["<outdir>"]
            infer_field(rcf, stagedir, outdir)
        elseif args["infer-object"]
            outdir = args["<outdir>"]
            objid = args["<objid>"]
            infer_field(rcf, stagedir, outdir; objid=objid)
        elseif args["score-field"]
            resultdir = args["<resultdir>"]
            truthfile = args["<truthfile>"]
            score_field_disk(rcf, resultdir, truthfile, stagedir)
        elseif args["score-object"]
            objid = args["<objid>"]
            resultdir = args["<resultdir>"]
            truthfile = args["<truthfile>"]
            score_object_disk(rcf, objid, resultdir, truthfile, stagedir)
        end
    end
end

main()