text:
using OffsetArrays

include("rp1_acoustics.jl")
include("inlinelimiter.jl")

# ===================================================================
function step1(parms::CParam,num_eqn::Int, num_waves::Int, num_ghost::Int, num_aux::Int, mx::Int,
               q::OffsetArray{Float64}, aux::OffsetArray{Float64},
               dx::Float64, dt::Float64, method::Array{Int,1}, mthlim::Array{Int},
               f::OffsetArray{Float64}, wave::OffsetArray{Float64}, s::OffsetArray{Float64},
               amdq::OffsetArray{Float64}, apdq::OffsetArray{Float64}, dtdx::OffsetArray{Float64},
               use_fwave::Bool, rp1)
# ===================================================================

    # Take one time step, updating q.

    # method(1) = 1   ==>  Godunov method
    # method(1) = 2   ==>  Slope limiter method
    # mthlim(p)  controls what limiter is used in the pth family
#
#
#     amdq, apdq, wave, s, and f are used locally:
#
#     amdq(num_eqn, 1-num_ghost:mx+num_ghost) = left-going flux-differences
#     apdq(num_eqn, 1-num_ghost:mx+num_ghost) = right-going flux-differences
#        e.g. amdq(m,i) = m'th component of A^- \Delta q from i'th Riemann
#                         problem (between cells i-1 and i).
#
#     wave(num_eqn, num_waves, 1-num_ghost:mx+num_ghost) = waves from solution of
#                                           Riemann problems,
#            wave(m,mw,i) = mth component of jump in q across
#                           wave in family mw in Riemann problem between
#                           states i-1 and i.
#
#     s(num_waves, 1-num_ghost:mx+num_ghost) = wave speeds,
#            s(m,iw) = speed of wave in family mw in Riemann problem between
#                      states i-1 and i.
#
#     f(num_eqn, 1-num_ghost:mx+num_ghost) = correction fluxes for second order method
#            f(m,i) = mth component of flux at left edge of ith cell
#     --------------------------------------------------------------------
#
#    implicit double precision (a-h,o-z)
#    double precision :: q(num_eqn,1-num_ghost:mx+num_ghost)
#    double precision :: aux(num_aux,1-num_ghost:mx+num_ghost)
#    double precision :: f(num_eqn,1-num_ghost:mx+num_ghost)
#    double precision :: s(num_waves,1-num_ghost:mx+num_ghost)
#    double precision :: wave(num_eqn, num_waves,1-num_ghost:mx+num_ghost)
#    double precision :: amdq(num_eqn,1-num_ghost:mx+num_ghost)
#    double precision :: apdq(num_eqn,1-num_ghost:mx+num_ghost)
#    double precision :: dtdx(1-num_ghost:mx+num_ghost)

    cfl::Float64 = 0.0

#    dtdx = OffsetArray(Float64, 1-num_ghost:mx+num_ghost)
#    integer ::          method(7),mthlim(num_waves)
#    logical ::          use_fwave
#    logical :: limit
#    external :: rp1

#f2py intent(in,out) q
#f2py intent(out) cfl
#f2py intent(in) num_eqn
#f2py intent(in) num_ghost
#f2py intent(in) mx
#f2py optional f, amdq, apdq, dtdx, s, wave

    # check if any limiters are used:
    limit = false
    for mw=1:num_waves
        if mthlim[mw] > 0
            limit = true
            break
        end
    end

    index_capa = method[6]
    for i=1-num_ghost:mx+num_ghost
        if index_capa > 0
            if aux[index_capa,i] <= 0.0
                error("capa must be positive")
            end
            dtdx[i] = dt / (dx*aux[index_capa,i])
        else
            dtdx[i] = dt / dx
        end
    end

    # solve Riemann problem at each interface
    # -----------------------------------------

    rp1(parms,mx,num_eqn,num_waves,num_aux,num_ghost,mx,q,q,aux,aux,wave,s,amdq,apdq)

    # Modify q for Godunov update:
    # Note this may not correspond to a conservative flux-differencing
    # for equations not in conservation form.  It is conservative if
    # amdq + apdq = f(q(i)) - f(q(i-1)).

    # print *,"q before, from fortran",q(24,1)

    for i = 1 : mx+1
        # q(:,i-1) is still in cache from last cycle of i loop, so
        # update it first
        for m = 1 : num_eqn
            q[m,i-1] = q[m,i-1] - dtdx[i-1]*amdq[m,i]
        end
        for m = 1 : num_eqn
            q[m,i] = q[m,i] - dtdx[i]*apdq[m,i]
        end
    end

    # compute maximum wave speed:
    cfl = 0.0
    for i=1:mx+1
        for mw=1:num_waves
            # if s>0 use dtdx(i) to compute CFL,
            # if s<0 use dtdx(i-1) to compute CFL:
            #cfl = dmax1(cfl, dtdx(i)*s(mw,i), -dtdx(i-1)*s(mw,i))
            #cfl = max(cfl, dtdx[i]*s[mw,i], -dtdx[i-1]*s[mw,i])
            s1 = s[mw,i]
            cfl = s1 > 0.0 ? s1*dtdx[i] : -s1*dtdx[i-1]
        end
    end

    if method[2] == 1
        return
    end

    # compute correction fluxes for second order q_{xx} terms:
    #     ----------------------------------------------------------

    # apply limiter to waves:
    if limit
        limiter(mx,num_eqn,num_waves,num_ghost,mx,wave,s,mthlim)
    end

    if !use_fwave
        for i=1:mx+1
            for m = 1:num_eqn
                f[m,i] = 0.0
            end
            dtdxave = 0.5 * (dtdx[i-1] + dtdx[i])
            for mw=1:num_waves
                for m=1:num_eqn
                    f[m,i] += 0.5 * abs(s[mw,i]) * (1.0 - abs(s[mw,i])*dtdxave) * wave[m,mw,i]
                end
            end
        end
    else
        for i=1:mx+1
            for m = 1:num_eqn
                f[m,i] = 0.0
            end
            dtdxave = 0.5 * (dtdx[i-1] + dtdx[i])
            for mw=1:num_waves
                for m=1:num_eqn
                    f[m,i] += 0.5 * sign(s[mw,i]) * (1.0 - abs(s[mw,i])*dtdxave) * wave[m,mw,i]
                end
            end
        end
    end

    # update q by differencing correction fluxes
    # ============================================

    # (Note:  Godunov update has already been performed above)

    for i = 1: mx+1
        for m = 1: num_eqn
            q[m,i] -= dtdx[i] * (f[m,i+1] - f[m,i])
        end
    end

    return cfl

end # subroutine step1


text:
# adapted from http://pcjericks.github.io/py-gdalogr-cookbook/geometry.html
using Base.Test
import GDAL, ArchGDAL; const AG = ArchGDAL

@testset "Create a Point" begin
    # Method 1
    AG.createpoint(1198054.34, 648493.09) do point
        println(AG.toWKT(point))
    end
    AG.createpoint((1198054.34, 648493.09)) do point
        println(AG.toWKT(point))
    end

    # Method 2
    AG.createpoint() do point
        AG.addpoint!(point, 1198054.34, 648493.09)
        println(AG.toWKT(point))
    end

    # Method 3
    AG.creategeom(AG.wkbPoint) do point
        AG.addpoint!(point, 1198054.34, 648493.09)
        println(AG.toWKT(point))
    end

    # Method 4
    point = AG.unsafe_creategeom(AG.wkbPoint)
        AG.addpoint!(point, 1198054.34, 648493.09)
        println(AG.toWKT(point))
    AG.destroy(point)
end

@testset "Create a LineString" begin
    # Method 1
    AG.createlinestring([(1116651.439379124,  637392.6969887456),
                         (1188804.0108498496, 652655.7409537067),
                         (1226730.3625203592, 634155.0816022386),
                         (1281307.30760719,   636467.6640211721)]) do line
        println(AG.toWKT(line))
    end

    # Method 2
    AG.createlinestring() do line
        AG.addpoint!(line, 1116651.439379124,  637392.6969887456)
        AG.addpoint!(line, 1188804.0108498496, 652655.7409537067)
        AG.addpoint!(line, 1226730.3625203592, 634155.0816022386)
        AG.addpoint!(line, 1281307.30760719,   636467.6640211721)
        println(AG.toWKT(line))
    end

    # Method 3
    AG.creategeom(AG.wkbLineString) do line
        AG.addpoint!(line, 1116651.439379124,  637392.6969887456)
        AG.addpoint!(line, 1188804.0108498496, 652655.7409537067)
        AG.addpoint!(line, 1226730.3625203592, 634155.0816022386)
        AG.addpoint!(line, 1281307.30760719,   636467.6640211721)
        println(AG.toWKT(line))
    end

    # Method 4
    line = AG.unsafe_creategeom(AG.wkbLineString)
        AG.addpoint!(line, 1116651.439379124,  637392.6969887456)
        AG.addpoint!(line, 1188804.0108498496, 652655.7409537067)
        AG.addpoint!(line, 1226730.3625203592, 634155.0816022386)
        AG.addpoint!(line, 1281307.30760719,   636467.6640211721)
        println(AG.toWKT(line))
    AG.destroy(line)
end

@testset "Create a Polygon" begin
    # Method 1
    AG.createpolygon([(1179091.1646903288, 712782.8838459781),
                      (1161053.0218226474, 667456.2684348812),
                      (1214704.933941905,  641092.8288590391),
                      (1228580.428455506,  682719.3123998424),
                      (1218405.0658121984, 721108.1805541387),
                      (1179091.1646903288, 712782.8838459781)]) do poly
        println(AG.toWKT(poly))
    end

    # Method 2
    AG.createpolygon() do poly
        ring = AG.unsafe_createlinearring(
                    [(1179091.1646903288, 712782.8838459781),
                     (1161053.0218226474, 667456.2684348812),
                     (1214704.933941905,  641092.8288590391),
                     (1228580.428455506,  682719.3123998424),
                     (1218405.0658121984, 721108.1805541387),
                     (1179091.1646903288, 712782.8838459781)])
        AG.addgeomdirectly!(poly, ring)
    end

    AG.createlinearring(
                    [(1179091.1646903288, 712782.8838459781),
                     (1161053.0218226474, 667456.2684348812),
                     (1214704.933941905,  641092.8288590391),
                     (1228580.428455506,  682719.3123998424),
                     (1218405.0658121984, 721108.1805541387),
                     (1179091.1646903288, 712782.8838459781)]) do ring
        AG.createpolygon() do poly
            AG.addgeom!(poly, ring)
            println(AG.toWKT(poly))
        end
    end

    # Method 3
    AG.creategeom(AG.wkbLinearRing) do ring
        AG.addpoint!(ring, 1179091.1646903288, 712782.8838459781)
        AG.addpoint!(ring, 1161053.0218226474, 667456.2684348812)
        AG.addpoint!(ring, 1214704.933941905, 641092.8288590391)
        AG.addpoint!(ring, 1228580.428455506, 682719.3123998424)
        AG.addpoint!(ring, 1218405.0658121984, 721108.1805541387)
        AG.addpoint!(ring, 1179091.1646903288, 712782.8838459781)

        AG.creategeom(AG.wkbPolygon) do poly
            AG.addgeom!(poly, ring)
            println(AG.toWKT(poly))
        end
    end

    # Method 4
    ring = AG.unsafe_creategeom(AG.wkbLinearRing)
        AG.addpoint!(ring, 1179091.1646903288, 712782.8838459781)
        AG.addpoint!(ring, 1161053.0218226474, 667456.2684348812)
        AG.addpoint!(ring, 1214704.933941905, 641092.8288590391)
        AG.addpoint!(ring, 1228580.428455506, 682719.3123998424)
        AG.addpoint!(ring, 1218405.0658121984, 721108.1805541387)
        AG.addpoint!(ring, 1179091.1646903288, 712782.8838459781)

    poly = AG.unsafe_creategeom(AG.wkbPolygon)
        AG.addgeomdirectly!(poly, ring) # ownership of ring passed to poly
        println(AG.toWKT(poly))
    AG.destroy(poly)
end

@testset "Create a Polygon with holes" begin
    # Method 1
    AG.createpolygon([# outring
                      [(1154115.274565847,  686419.4442701361),
                       (1154115.274565847,  653118.2574374934),
                       (1165678.1866605144, 653118.2574374934),
                       (1165678.1866605144, 686419.4442701361),
                       (1154115.274565847,  686419.4442701361)],
                      # innerring(s)
                      [(1149490.1097279799, 691044.6091080031),
                       (1149490.1097279799, 648030.5761158396),
                       (1191579.1097525698, 648030.5761158396),
                       (1191579.1097525698, 691044.6091080031),
                       (1149490.1097279799, 691044.6091080031)]
                     ]) do poly
        println(AG.toWKT(poly))
    end

    # Method 2
    AG.createpolygon() do poly
        outring = AG.unsafe_creategeom(AG.wkbLinearRing)
            AG.addpoint!(outring, 1154115.274565847, 686419.4442701361)
            AG.addpoint!(outring, 1154115.274565847, 653118.2574374934)
            AG.addpoint!(outring, 1165678.1866605144, 653118.2574374934)
            AG.addpoint!(outring, 1165678.1866605144, 686419.4442701361)
            AG.addpoint!(outring, 1154115.274565847, 686419.4442701361)

        innerring = AG.unsafe_creategeom(AG.wkbLinearRing)
            AG.addpoint!(innerring, 1149490.1097279799, 691044.6091080031)
            AG.addpoint!(innerring, 1149490.1097279799, 648030.5761158396)
            AG.addpoint!(innerring, 1191579.1097525698, 648030.5761158396)
            AG.addpoint!(innerring, 1191579.1097525698, 691044.6091080031)
            AG.addpoint!(innerring, 1149490.1097279799, 691044.6091080031)

        AG.addgeomdirectly!(poly, outring)
        AG.addgeomdirectly!(poly, innerring)
        println(AG.toWKT(poly))
    end

    # Method 3
    AG.creategeom(AG.wkbPolygon) do poly
        outring = AG.unsafe_creategeom(AG.wkbLinearRing)
            AG.addpoint!(outring, 1154115.274565847, 686419.4442701361)
            AG.addpoint!(outring, 1154115.274565847, 653118.2574374934)
            AG.addpoint!(outring, 1165678.1866605144, 653118.2574374934)
            AG.addpoint!(outring, 1165678.1866605144, 686419.4442701361)
            AG.addpoint!(outring, 1154115.274565847, 686419.4442701361)

        innerring = AG.unsafe_creategeom(AG.wkbLinearRing)
            AG.addpoint!(innerring, 1149490.1097279799, 691044.6091080031)
            AG.addpoint!(innerring, 1149490.1097279799, 648030.5761158396)
            AG.addpoint!(innerring, 1191579.1097525698, 648030.5761158396)
            AG.addpoint!(innerring, 1191579.1097525698, 691044.6091080031)
            AG.addpoint!(innerring, 1149490.1097279799, 691044.6091080031)

        AG.addgeomdirectly!(poly, outring)
        AG.addgeomdirectly!(poly, innerring)
        println(AG.toWKT(poly))
    end

    # Method 4
    poly = AG.unsafe_creategeom(AG.wkbPolygon)
        outring = AG.unsafe_creategeom(AG.wkbLinearRing)
            AG.addpoint!(outring, 1154115.274565847, 686419.4442701361)
            AG.addpoint!(outring, 1154115.274565847, 653118.2574374934)
            AG.addpoint!(outring, 1165678.1866605144, 653118.2574374934)
            AG.addpoint!(outring, 1165678.1866605144, 686419.4442701361)
            AG.addpoint!(outring, 1154115.274565847, 686419.4442701361)

        innerring = AG.unsafe_creategeom(AG.wkbLinearRing)
            AG.addpoint!(innerring, 1149490.1097279799, 691044.6091080031)
            AG.addpoint!(innerring, 1149490.1097279799, 648030.5761158396)
            AG.addpoint!(innerring, 1191579.1097525698, 648030.5761158396)
            AG.addpoint!(innerring, 1191579.1097525698, 691044.6091080031)
            AG.addpoint!(innerring, 1149490.1097279799, 691044.6091080031)

        AG.addgeomdirectly!(poly, outring)
        AG.addgeomdirectly!(poly, innerring)
        println(AG.toWKT(poly))
    AG.destroy(poly)
end

@testset "Create a MultiPoint" begin
    # Method 1
    AG.createmultipoint([(1251243.7361610543, 598078.7958668759),
                         (1240605.8570339603, 601778.9277371694),
                         (1250318.7031934808, 606404.0925750365)]) do multipoint
        println(AG.toWKT(multipoint))
    end

    # Method 2
    println("method 2")
    AG.createmultipoint() do multipoint
        point1 = AG.unsafe_createpoint(1251243.7361610543, 598078.7958668759)
        point2 = AG.unsafe_createpoint(1240605.8570339603, 601778.9277371694)
        point3 = AG.unsafe_createpoint(1250318.7031934808, 606404.0925750365)
        AG.addgeomdirectly!(multipoint, point1)
        AG.addgeomdirectly!(multipoint, point2)
        AG.addgeomdirectly!(multipoint, point3)
        println(AG.toWKT(multipoint))
    end

    # Method 3
    AG.creategeom(AG.wkbMultiPoint) do multipoint
        point1 = AG.unsafe_createpoint(1251243.7361610543, 598078.7958668759)
        point2 = AG.unsafe_createpoint(1240605.8570339603, 601778.9277371694)
        point3 = AG.unsafe_createpoint(1250318.7031934808, 606404.0925750365)
        AG.addgeomdirectly!(multipoint, point1)
        AG.addgeomdirectly!(multipoint, point2)
        AG.addgeomdirectly!(multipoint, point3)
        println(AG.toWKT(multipoint))
    end

    # Method 4
    multipoint = AG.unsafe_creategeom(AG.wkbMultiPoint)
        point1 = AG.unsafe_createpoint(1251243.7361610543, 598078.7958668759)
        point2 = AG.unsafe_createpoint(1240605.8570339603, 601778.9277371694)
        point3 = AG.unsafe_createpoint(1250318.7031934808, 606404.0925750365)
        AG.addgeomdirectly!(multipoint, point1)
        AG.addgeomdirectly!(multipoint, point2)
        AG.addgeomdirectly!(multipoint, point3)
        println(AG.toWKT(multipoint))
    AG.destroy(multipoint)
end

@testset "Create a MultiLineString" begin
    # Method 1
    AG.createmultilinestring(Vector{Tuple{Float64,Float64}}[
                              [(1214242.4174581182, 617041.9717021306),
                               (1234593.142744733,  629529.9167643716)],
                              [(1184641.3624957693, 626754.8178616514),
                               (1219792.6152635587, 606866.6090588232)]
                             ]) do multiline
        println(AG.toWKT(multiline))
    end

    # Method 2
    AG.createmultilinestring() do multiline
        AG.addgeomdirectly!(multiline, AG.unsafe_createlinestring(
            [(1214242.4174581182, 617041.9717021306),
             (1234593.142744733,  629529.9167643716)]))
        AG.addgeomdirectly!(multiline, AG.unsafe_createlinestring(
            [(1184641.3624957693, 626754.8178616514),
             (1219792.6152635587, 606866.6090588232)]))
        println(AG.toWKT(multiline))
    end

    # Method 3
    AG.creategeom(AG.wkbMultiLineString) do multiline
        line = AG.unsafe_creategeom(AG.wkbLineString)
        AG.addpoint!(line, 1214242.4174581182, 617041.9717021306)
        AG.addpoint!(line, 1234593.142744733, 629529.9167643716)
        AG.addgeomdirectly!(multiline, line)

        line = AG.unsafe_creategeom(AG.wkbLineString)
        AG.addpoint!(line, 1184641.3624957693, 626754.8178616514)
        AG.addpoint!(line, 1219792.6152635587, 606866.6090588232)
        AG.addgeomdirectly!(multiline, line)

        println(AG.toWKT(multiline))
    end

    # Method 4
    multiline = AG.unsafe_creategeom(AG.wkbMultiLineString)
        line = AG.unsafe_creategeom(AG.wkbLineString)
        AG.addpoint!(line, 1214242.4174581182, 617041.9717021306)
        AG.addpoint!(line, 1234593.142744733, 629529.9167643716)
        AG.addgeomdirectly!(multiline, line)

        line = AG.unsafe_creategeom(AG.wkbLineString)
        AG.addpoint!(line, 1184641.3624957693, 626754.8178616514)
        AG.addpoint!(line, 1219792.6152635587, 606866.6090588232)
        AG.addgeomdirectly!(multiline, line)

        println(AG.toWKT(multiline))
    AG.destroy(multiline)
end

@testset "Create a MultiPolygon" begin
    # Method 1
    AG.createmultipolygon_noholes(Vector{Tuple{Float64,Float64}}[
             [(1204067.0548148106, 634617.5980860253),
              (1204067.0548148106, 620742.1035724243),
              (1215167.4504256917, 620742.1035724243),
              (1215167.4504256917, 634617.5980860253),
              (1204067.0548148106, 634617.5980860253)],
             [(1179553.6811741155, 647105.5431482664),
              (1179553.6811741155, 626292.3013778647),
              (1194354.20865529,   626292.3013778647),
              (1194354.20865529,   647105.5431482664),
              (1179553.6811741155, 647105.5431482664)]
            ]) do multipolygon
        println(AG.toWKT(multipolygon))
    end

    # Method 2
    AG.createmultipolygon() do multipolygon
        poly = AG.unsafe_createpolygon(
            [(1204067.0548148106, 634617.5980860253),
             (1204067.0548148106, 620742.1035724243),
             (1215167.4504256917, 620742.1035724243),
             (1215167.4504256917, 634617.5980860253),
             (1204067.0548148106, 634617.5980860253)])
        AG.addgeomdirectly!(multipolygon,poly)

        poly = AG.unsafe_createpolygon(
            [(1179553.6811741155, 647105.5431482664),
             (1179553.6811741155, 626292.3013778647),
             (1194354.20865529,   626292.3013778647),
             (1194354.20865529,   647105.5431482664),
             (1179553.6811741155, 647105.5431482664)])
        AG.addgeomdirectly!(multipolygon,poly)

        println(AG.toWKT(multipolygon))
    end

    # Method 3
    AG.creategeom(AG.wkbMultiPolygon) do multipolygon
        poly = AG.unsafe_creategeom(AG.wkbPolygon)
            ring = AG.unsafe_creategeom(AG.wkbLinearRing)
                AG.addpoint!(ring, 1204067.0548148106, 634617.5980860253)
                AG.addpoint!(ring, 1204067.0548148106, 620742.1035724243)
                AG.addpoint!(ring, 1215167.4504256917, 620742.1035724243)
                AG.addpoint!(ring, 1215167.4504256917, 634617.5980860253)
                AG.addpoint!(ring, 1204067.0548148106, 634617.5980860253)
            AG.addgeomdirectly!(poly, ring)
        AG.addgeomdirectly!(multipolygon,poly)

        poly = AG.unsafe_creategeom(AG.wkbPolygon)
            ring = AG.unsafe_creategeom(AG.wkbLinearRing)
                AG.addpoint!(ring, 1179553.6811741155, 647105.5431482664)
                AG.addpoint!(ring, 1179553.6811741155, 626292.3013778647)
                AG.addpoint!(ring, 1194354.20865529,   626292.3013778647)
                AG.addpoint!(ring, 1194354.20865529,   647105.5431482664)
                AG.addpoint!(ring, 1179553.6811741155, 647105.5431482664)
            AG.addgeomdirectly!(poly, ring)
        AG.addgeomdirectly!(multipolygon,poly)

        println(AG.toWKT(multipolygon))
    end

    # Method 4
    multipolygon = AG.unsafe_creategeom(AG.wkbMultiPolygon)
        poly = AG.unsafe_creategeom(AG.wkbPolygon)
            ring = AG.unsafe_creategeom(AG.wkbLinearRing)
                AG.addpoint!(ring, 1204067.0548148106, 634617.5980860253)
                AG.addpoint!(ring, 1204067.0548148106, 620742.1035724243)
                AG.addpoint!(ring, 1215167.4504256917, 620742.1035724243)
                AG.addpoint!(ring, 1215167.4504256917, 634617.5980860253)
                AG.addpoint!(ring, 1204067.0548148106, 634617.5980860253)
            AG.addgeomdirectly!(poly, ring)
        AG.addgeomdirectly!(multipolygon,poly)

        poly = AG.unsafe_creategeom(AG.wkbPolygon)
            ring = AG.unsafe_creategeom(AG.wkbLinearRing)
                AG.addpoint!(ring, 1179553.6811741155, 647105.5431482664)
                AG.addpoint!(ring, 1179553.6811741155, 626292.3013778647)
                AG.addpoint!(ring, 1194354.20865529,   626292.3013778647)
                AG.addpoint!(ring, 1194354.20865529,   647105.5431482664)
                AG.addpoint!(ring, 1179553.6811741155, 647105.5431482664)
            AG.addgeomdirectly!(poly, ring)
        AG.addgeomdirectly!(multipolygon,poly)

        println(AG.toWKT(multipolygon))
    AG.destroy(multipolygon)
end

@testset "Create a GeometryCollection" begin
    # Method 2
    AG.creategeomcollection() do geomcol
        for g in [AG.unsafe_createpoint(-122.23, 47.09),
                  AG.unsafe_createlinestring([(-122.60, 47.14),
                                              (-122.48, 47.23)])]
            AG.addgeomdirectly!(geomcol, g)
        end
        println(AG.toWKT(geomcol))
    end

    # Method 3
    AG.creategeom(AG.wkbGeometryCollection) do geomcol
        point = AG.unsafe_creategeom(AG.wkbPoint)
            AG.addpoint!(point, -122.23, 47.09)
        AG.addgeomdirectly!(geomcol, point)

        line = AG.unsafe_creategeom(AG.wkbLineString)
            AG.addpoint!(line, -122.60, 47.14)
            AG.addpoint!(line, -122.48, 47.23)
        AG.addgeomdirectly!(geomcol, line)

        println(AG.toWKT(geomcol))
    end

    # Method 4
    geomcol =  AG.unsafe_creategeom(AG.wkbGeometryCollection)
        point = AG.unsafe_creategeom(AG.wkbPoint)
            AG.addpoint!(point, -122.23, 47.09)
        AG.addgeomdirectly!(geomcol, point)

        line = AG.unsafe_creategeom(AG.wkbLineString)
            AG.addpoint!(line, -122.60, 47.14)
            AG.addpoint!(line, -122.48, 47.23)
        AG.addgeomdirectly!(geomcol, line)

        println(AG.toWKT(geomcol))
    AG.destroy(geomcol)
end

@testset "Create Geometry from WKT" begin
    wkt = "POINT (1120351.5712494177 741921.4223245403)"

    # Method 1
    point = AG.unsafe_fromWKT(wkt)
        println((AG.getx(point, 0), AG.gety(point, 0)))
    AG.destroy(point)

    # Method 2
    AG.fromWKT(wkt) do point
        println((AG.getx(point, 0), AG.gety(point, 0)))
    end
end

@testset "Create Geometry from GeoJSON" begin
    geojson = """{"type":"Point","coordinates":[108420.33,753808.59]}"""

    # Method 1
    point = AG.unsafe_fromJSON(geojson)
        println((AG.getx(point, 0), AG.gety(point, 0)))
    AG.destroy(point)

    # Method 2
    AG.fromJSON(geojson) do point
        println((AG.getx(point, 0), AG.gety(point, 0)))
    end
end

@testset "Create Geometry from GML" begin
    gml = """<gml:Point xmlns:gml="http://www.opengis.net/gml">
             <gml:coordinates>108420.33,753808.59</gml:coordinates>
             </gml:Point>"""

    # Method 1
    point = AG.unsafe_fromGML(gml)
        println((AG.getx(point, 0), AG.gety(point, 0)))
    AG.destroy(point)

    # Method 2
    AG.fromGML(gml) do point
        println((AG.getx(point, 0), AG.gety(point, 0)))
    end
end

@testset "Create Geometry from WKB" begin
    wkb = [0x01,0x01,0x00,0x00,0x00,0x7b,0x14,0xae,0x47,0x45,0x78,0xfa,0x40,
           0xe1,0x7a,0x14,0x2e,0x21,0x01,0x27,0x41]

    # Method 1
    point = AG.unsafe_fromWKB(wkb)
        println((AG.getx(point, 0), AG.gety(point, 0)))
    AG.destroy(point)

    # Method 2
    AG.fromWKB(wkb) do point
        println((AG.getx(point, 0), AG.gety(point, 0)))
    end
end

@testset "Count Points in a Geometry" begin
    wkt = "LINESTRING (1181866.263593049 615654.4222507705, 1205917.1207499576 623979.7189589312, 1227192.8790041457 643405.4112779726, 1224880.2965852122 665143.6860159477)"
    AG.fromWKT(wkt) do geom
        println("Geometry has $(AG.npoint(geom)) points")
    end
end

@testset "Count Geometries in a Geometry" begin
    wkt = "MULTIPOINT (1181866.263593049 615654.4222507705, 1205917.1207499576 623979.7189589312, 1227192.8790041457 643405.4112779726, 1224880.2965852122 665143.6860159477)"
    AG.fromWKT(wkt) do geom
        println("Geometry has $(AG.npoint(geom)) points")
    end
end

@testset "Iterate over Geometries in a Geometry" begin
    wkt = "MULTIPOINT (1181866.263593049 615654.4222507705, 1205917.1207499576 623979.7189589312, 1227192.8790041457 643405.4112779726, 1224880.2965852122 665143.6860159477)"
    AG.fromWKT(wkt) do geom
        for i in 0:(AG.ngeom(geom)-1)
            println("$i). $(AG.toWKT(AG.getgeom(geom,i)))")
        end
    end
end

@testset "Iterate over Points in a Geometry" begin
    wkt = "LINESTRING (1181866.263593049 615654.4222507705, 1205917.1207499576 623979.7189589312, 1227192.8790041457 643405.4112779726, 1224880.2965852122 665143.6860159477)"
    AG.fromWKT(wkt) do geom
        for i in 0:(AG.npoint(geom)-1)
            println("$i). POINT $(AG.getpoint(geom, i))")
        end
    end
end

@testset "Buffer a Geometry" begin
    wkt = "POINT (1198054.34 648493.09)"
    AG.fromWKT(wkt) do pt
        bufferdist = 500
        AG.buffer(pt, bufferdist) do poly
            println("$(AG.toWKT(pt)) buffered by $bufferdist is ")
            println("$(AG.toWKT(poly))")
        end
    end
end

# @testset "Calculate Envelope of a Geometry" begin
#     wkt = "LINESTRING (1181866.263593049 615654.4222507705, 1205917.1207499576 623979.7189589312, 1227192.8790041457 643405.4112779726, 1224880.2965852122 665143.6860159477)"
#     geom = ogr.CreateGeometryFromWkt(wkt)
#     # Get Envelope returns a tuple (minX, maxX, minY, maxY)
#     env = geom.GetEnvelope()
#     print "minX: %d, minY: %d, maxX: %d, maxY: %d" %(env[0],env[2],env[1],env[3])
# end

@testset "Calculate the Area of a Geometry" begin
    wkt = "POLYGON ((1162440.5712740074 672081.4332727483, 1162440.5712740074 647105.5431482664, 1195279.2416228633 647105.5431482664, 1195279.2416228633 672081.4332727483, 1162440.5712740074 672081.4332727483))"
    AG.fromWKT(wkt) do poly
        println("Area = $(AG.geomarea(poly))")
    end
end

@testset "Calculate the Length of a Geometry" begin
    wkt = "LINESTRING (1181866.263593049 615654.4222507705, 1205917.1207499576 623979.7189589312, 1227192.8790041457 643405.4112779726, 1224880.2965852122 665143.6860159477)"
    AG.fromWKT(wkt) do line
        println("Length = $(AG.geomlength(line))")
    end
end

@testset "Get the geometry type (as a string) from a Geometry" begin
    wkts = [
        "POINT (1198054.34 648493.09)",
        "LINESTRING (1181866.263593049 615654.4222507705, 1205917.1207499576 623979.7189589312, 1227192.8790041457 643405.4112779726, 1224880.2965852122 665143.6860159477)",
        "POLYGON ((1162440.5712740074 672081.4332727483, 1162440.5712740074 647105.5431482664, 1195279.2416228633 647105.5431482664, 1195279.2416228633 672081.4332727483, 1162440.5712740074 672081.4332727483))"
    ]
    for wkt in wkts
        AG.fromWKT(wkt) do geom
            println(AG.getgeomname(geom))
        end
    end
end

@testset "Calculate intersection between two Geometries" begin
    wkt1 = "POLYGON ((1208064.271243039 624154.6783778917, 1208064.271243039 601260.9785661874, 1231345.9998651114 601260.9785661874, 1231345.9998651114 624154.6783778917, 1208064.271243039 624154.6783778917))"
    wkt2 = "POLYGON ((1199915.6662253144 633079.3410163528, 1199915.6662253144 614453.958118695, 1219317.1067437078 614453.958118695, 1219317.1067437078 633079.3410163528, 1199915.6662253144 633079.3410163528)))"

    AG.fromWKT(wkt1) do poly1
    AG.fromWKT(wkt2) do poly2
        AG.intersection(poly1, poly2) do poly3
            println(AG.toWKT(poly3))
        end
    end
    end
end

@testset "Calculate union between two Geometries" begin
    wkt1 = "POLYGON ((1208064.271243039 624154.6783778917, 1208064.271243039 601260.9785661874, 1231345.9998651114 601260.9785661874, 1231345.9998651114 624154.6783778917, 1208064.271243039 624154.6783778917))"
    wkt2 = "POLYGON ((1199915.6662253144 633079.3410163528, 1199915.6662253144 614453.958118695, 1219317.1067437078 614453.958118695, 1219317.1067437078 633079.3410163528, 1199915.6662253144 633079.3410163528)))"
    AG.fromWKT(wkt1) do poly1
    AG.fromWKT(wkt2) do poly2
        AG.union(poly1, poly2) do poly3
            println(AG.toWKT(poly1))
            println(AG.toWKT(poly2))
            println(AG.toWKT(poly3))
        end
    end
    end
end

@testset "Write Geometry to GeoJSON|GML|WKT|WKB" begin
    AG.createpolygon([(1179091.1646903288, 712782.8838459781),
                      (1161053.0218226474, 667456.2684348812),
                      (1214704.933941905,  641092.8288590391),
                      (1228580.428455506,  682719.3123998424),
                      (1218405.0658121984, 721108.1805541387),
                      (1179091.1646903288, 712782.8838459781)]) do poly
        println(AG.toJSON(poly))
        println(AG.toGML(poly))
        println(AG.toWKT(poly))
        println(AG.toWKB(poly))
    end
end

@testset "Force polygon to multipolygon" begin
    wkt = "POLYGON ((1179091.164690328761935 712782.883845978067257,1161053.021822647424415 667456.268434881232679,1214704.933941904921085 641092.828859039116651,1228580.428455505985767 682719.312399842427112,1218405.065812198445201 721108.180554138729349,1179091.164690328761935 712782.883845978067257))"
    AG.fromWKT([wkt]) do poly
        println("Before: $(AG.toWKT(poly))")
        if AG.getgeomtype(poly) == AG.wkbPolygon
            AG.forceto(poly, AG.wkbMultiPolygon) do mpoly
                println("After: $mpoly")
            end
        end
    end
end

ext: jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AzureReachabilityReportLatencyInfo <: SwaggerModel
    timeStamp::Any # spec type: Union{ Nothing, DateTime } # spec name: timeStamp
    score::Any # spec type: Union{ Nothing, Int32 } # spec name: score

    function AzureReachabilityReportLatencyInfo(;timeStamp=nothing, score=nothing)
        o = new()
        validate_property(AzureReachabilityReportLatencyInfo, Symbol("timeStamp"), timeStamp)
        setfield!(o, Symbol("timeStamp"), timeStamp)
        validate_property(AzureReachabilityReportLatencyInfo, Symbol("score"), score)
        setfield!(o, Symbol("score"), score)
        o
    end
end # type AzureReachabilityReportLatencyInfo

const _property_map_AzureReachabilityReportLatencyInfo = Dict{Symbol,Symbol}(Symbol("timeStamp")=>Symbol("timeStamp"), Symbol("score")=>Symbol("score"))
const _property_types_AzureReachabilityReportLatencyInfo = Dict{Symbol,String}(Symbol("timeStamp")=>"DateTime", Symbol("score")=>"Int32")
Base.propertynames(::Type{ AzureReachabilityReportLatencyInfo }) = collect(keys(_property_map_AzureReachabilityReportLatencyInfo))
Swagger.property_type(::Type{ AzureReachabilityReportLatencyInfo }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_AzureReachabilityReportLatencyInfo[name]))}
Swagger.field_name(::Type{ AzureReachabilityReportLatencyInfo }, property_name::Symbol) =  _property_map_AzureReachabilityReportLatencyInfo[property_name]

function check_required(o::AzureReachabilityReportLatencyInfo)
    true
end

function validate_property(::Type{ AzureReachabilityReportLatencyInfo }, name::Symbol, val)
    if name === Symbol("score")
        Swagger.validate_param(name, "AzureReachabilityReportLatencyInfo", :maximum, val, 100, false)
        Swagger.validate_param(name, "AzureReachabilityReportLatencyInfo", :minimum, val, 1, false)
    end
end


text:
"""
Estimation of the covariance matrix in the standard simple way given returns.
https://en.wikipedia.org/wiki/Sample_mean_and_covariance
"""
function simple_covariance_given_returns(returns::Array{R,2}) where R<:Real
    N = size(returns)[2]
    mat = zeros(N,N)
    for i in 1:N
        for j in i:N
            if i == j
                mat[i,i] = var(returns[:,i])
            else
                mat[i,j] = cov(returns[:,i], returns[:,j])
            end
        end
    end
    return Hermitian(mat)
end

"""
Estimation of the covariance matrix in the standard simple way given a time grid.
"""
function simple_covariance_given_time_grid(ts::SortedDataFrame, assets::Vector{Symbol}, time_grid::Vector; regularisation::Union{Missing,Symbol} = :covariance_default,
                                           regularisation_params::Dict = Dict(), only_regulise_if_not_PSD::Bool = false)
    dd_compiled = latest_value(ts, time_grid; assets = assets)
    dd = get_returns(dd_compiled; rescale_for_duration = false)

    if nrow(dd) < 1 return make_nan_covariance_matrix(assets, ts.time_period_per_unit) end

    returns = Matrix(dd[:, assets])
    covariance = simple_covariance_given_returns(returns)

    # Regularisation
    dont_regulise = ismissing(regularisation) || (only_regulise_if_not_PSD && is_psd_matrix(covariance))
    covariance = dont_regulise ? covariance : regularise(covariance, ts, assets, regularisation; regularisation_params... )

    # av_spacing
    N = length(time_grid)
    spacing = safe_multiply_period(Statistics.mean(time_grid[2:N] .- time_grid[1:(N-1)]), ts.time_period_per_unit)

    # Packing into a CovarianceMatrix and returning.
    cor, vols = cov2cor_and_vol(covariance, spacing, ts.time_period_per_unit)
    return CovarianceMatrix(cor, vols, assets, ts.time_period_per_unit)
end

"""
    simple_covariance(ts::SortedDataFrame, assets::Vector{Symbol} = get_assets(ts);
                      regularisation::Union{Missing,Symbol} = :covariance_default,
                      regularisation_params::Dict = Dict(), only_regulise_if_not_PSD::Bool = false,
                      time_grid::Union{Missing,Vector} = missing, fixed_spacing::Union{Missing,<:Real} = missing,
                      refresh_times::Bool = false, rough_guess_number_of_intervals::Integer = 5)

Estimation of the covariance matrix in the standard textbook way.
### Inputs
* `ts` - The tick data.
* `assets` - The assets you want to estimate volatilities for.
* `regularisation` - A symbol representing what regularisation technique should be used. If missing no regularisation is performed.
* `regularisation_params` - keyword arguments to be consumed in the regularisation algorithm.
* `only_regulise_if_not_PSD` - Should regularisation only be attempted if the matrix is not psd already.
* `time_grid` - The grid with which to calculate returns.
* `fixed_spacing` - A spacing used to calculate a time grid. Not used if `refresh_times=true`.
* `refresh_times` - Should refresh times be used to estimate covariance.
* `rough_guess_number_of_intervals` - A rough number of intervals to calculate a default spacing. Not used if a `time_grid` or `fixed_spacing` is provided or if `refresh_times=true`.
### Returns
* A `CovarianceMatrix`.

### References
"""
function simple_covariance(ts::SortedDataFrame, assets::Vector{Symbol} = get_assets(ts);
                           regularisation::Union{Missing,Symbol} = :covariance_default,
                           regularisation_params::Dict = Dict(), only_regulise_if_not_PSD::Bool = false,
                           time_grid::Union{Missing,Vector} = missing, fixed_spacing::Union{Missing,<:Real} = missing,
                           refresh_times::Bool = false, rough_guess_number_of_intervals::Integer = 5)
   if ismissing(time_grid)
       time_grid = Vector{eltype(ts.df[:,ts.time])}()
       if refresh_times
           time_grid = get_all_refresh_times(ts, assets)
       elseif !ismissing(fixed_spacing)
            time_grid = collect(minimum(ts.df[:,ts.time]):fixed_spacing:maximum(ts.df[:,ts.time]))
       else
           n_grid = default_spacing(ts; rough_guess_number_of_intervals = rough_guess_number_of_intervals)
           vals = collect(values(n_grid))
           spacing = min( mean(vals[(isnan.(vals) .== false) .& (isinf.(vals) .== false)]) , duration(ts; in_dates_period = false)/20   )
           spacing = isnan(spacing) ? duration(ts; in_dates_period = false)/20 : spacing
           time_grid = collect(minimum(ts.df[:,ts.time]):spacing:maximum(ts.df[:,ts.time]))
       end
   end
   return simple_covariance_given_time_grid(ts, assets, time_grid; regularisation = regularisation, regularisation_params = regularisation_params, only_regulise_if_not_PSD = only_regulise_if_not_PSD)
end


text:
function SEMParser(arg0::JString)
    return SEMParser((JString,), arg0)
end

function SEMParser(arg0::JString, arg1::DataProvidersManager, arg2::TimeScales)
    return SEMParser((JString, DataProvidersManager, TimeScales), arg0, arg1, arg2)
end

function get_almanacs(obj::SEMParser)
    return jcall(obj, "getAlmanacs", List, ())
end

function get_prn_numbers(obj::SEMParser)
    return jcall(obj, "getPRNNumbers", List, ())
end

function get_supported_names(obj::SEMParser)
    return jcall(obj, "getSupportedNames", JString, ())
end

function load_data(obj::SEMParser)
    return jcall(obj, "loadData", void, ())
end

function load_data(obj::SEMParser, arg0::InputStream, arg1::JString)
    return jcall(obj, "loadData", void, (InputStream, JString), arg0, arg1)
end

function still_accepts_data(obj::SEMParser)
    return jcall(obj, "stillAcceptsData", jboolean, ())
end


text:
module ArgumentParse

const KEYWORD_ARGUMENT_PREFIX = "--"
const HELP_ARGUMENTS = ["--help", "-h"]

struct NoDefault end

struct ArgumentParsingError <: Exception
    message::String
end

struct ShowHelp <: Exception
    message::String
end

struct ArgumentSpecification
    name::String
    argument_string::String
    help::String
    arg_type::Type
    required::Bool
    default
    action::Symbol
end

struct ArgumentParser
    program_name::String
    propagate_errors::Bool
    positional_arguments::Vector{ArgumentSpecification}
    keyword_arguments::Dict{String, ArgumentSpecification}

    function ArgumentParser(; program_name=nothing, propagate_errors=false)
        if program_name == nothing
            program_name = (
                Base.source_path() != nothing ? basename(Base.source_path()) : "<PROGRAM>"
            )
        end
        new(
            program_name,
            propagate_errors,
            ArgumentSpecification[],
            Dict{String, ArgumentSpecification}()
        )
    end
end

strip_keyword_argument_prefix(name::String) = name[(length(KEYWORD_ARGUMENT_PREFIX) + 1):end]

function add_argument(
    parser::ArgumentParser, argument_string::String; help="", arg_type=String, required=nothing,
    default=NoDefault, action=:store
)
    is_keyword = startswith(argument_string, KEYWORD_ARGUMENT_PREFIX)
    if required == nothing
        required = !is_keyword
    end
    if required
        @assert default == NoDefault
    end
    if is_keyword
        name = strip_keyword_argument_prefix(argument_string)
    else
        @assert action == :store
        name = argument_string
    end
    if action == :store_true
        @assert default == NoDefault || default == false
        default = false
    end
    specification = ArgumentSpecification(
        name,
        argument_string,
        help,
        arg_type,
        required,
        default,
        action,
    )
    if is_keyword
        parser.keyword_arguments[name] = specification
    else
        push!(parser.positional_arguments, specification)
    end
end

function parse_args(parser::ArgumentParser, arguments::Vector{String})
    try
        return parse_args_helper(parser, Dict{String, Any}(), arguments)
    catch exc
        if parser.propagate_errors
            rethrow()
        elseif isa(exc, ArgumentParsingError)
            println(STDERR, exc.message)
            exit(1)
        elseif isa(exc, ShowHelp)
            println(STDERR, exc.message)
            exit(0)
        else
            rethrow()
        end
    end
end

function parse_args_helper(
    parser::ArgumentParser, parsed_args::Dict{String, Any}, arguments::Vector{String}
)
    while !isempty(arguments)
        next_argument = shift!(arguments)
        remaining_positional_args = get_unfilled_args(parser.positional_arguments, parsed_args)
        if in(next_argument, HELP_ARGUMENTS)
            show_help(parser)
        elseif startswith(next_argument, KEYWORD_ARGUMENT_PREFIX)
            name = strip_keyword_argument_prefix(next_argument)
            if !haskey(parser.keyword_arguments, name)
                parsing_error(parser, @sprintf("Unknown keyword '%s'", next_argument))
            else
                specification = parser.keyword_arguments[name]
                parse_keyword_argument(parser, parsed_args, specification, arguments)
            end
        elseif isempty(remaining_positional_args)
            parsing_error(parser, "Too many arguments")
        else
            specification = remaining_positional_args[1]
            parse_positional_argument(parser, parsed_args, specification, next_argument)
        end
    end

    remaining_required_args = vcat(
        get_unfilled_args(parser.positional_arguments, parsed_args),
        get_unfilled_args(values(parser.keyword_arguments), parsed_args),
    )
    remaining_required_args = filter(spec -> spec.required, remaining_required_args)
    if !isempty(remaining_required_args)
        remaining_names = join(
            [specification.argument_string for specification in remaining_required_args],
            ", ",
        )
        parsing_error(parser, @sprintf("Missing required arguments: %s", remaining_names))
    else
        set_defaults(parser, parsed_args)
        return parsed_args
    end
end

function get_unfilled_args(specifications, parsed_args::Dict{String, Any})
    if isempty(specifications)
        return []
    end
    mapreduce(vcat, specifications) do specification::ArgumentSpecification
        if !haskey(parsed_args, specification.name)
            [specification]
        else
            []
        end
    end
end

function parse_keyword_argument(
    parser::ArgumentParser, parsed_args::Dict{String, Any}, specification::ArgumentSpecification,
    following_arguments::Vector{String}
)
    if specification.action == :store_true
        parsed_args[specification.name] = true
    else
        @assert specification.action == :store
        if isempty(following_arguments)
            parsing_error(
                parser,
                @sprintf("Missing value for argument %s", specification.argument_string),
            )
        else
            raw_value = shift!(following_arguments)
            parsed_args[specification.name] = parse_value(specification, raw_value)
        end
    end
end

function parse_positional_argument(
    parser::ArgumentParser, parsed_args::Dict{String, Any}, specification::ArgumentSpecification,
    raw_value::String
)
    parsed_args[specification.name] = parse_value(specification, raw_value)
end

function parse_value(specification::ArgumentSpecification, raw_value::String)
    if specification.arg_type == String
        return raw_value
    else
        return parse(specification.arg_type, raw_value)
    end
end

function set_defaults(parser::ArgumentParser, parsed_args::Dict{String, Any})
    all_specifications = vcat(
        parser.positional_arguments,
        collect(values(parser.keyword_arguments)),
    )
    for specification in all_specifications
        if !haskey(parsed_args, specification.name) && specification.default != NoDefault
            parsed_args[specification.name] = specification.default
        end
    end
end

function parsing_error(parser::ArgumentParser, message::String)
    throw(ArgumentParsingError(@sprintf("ERROR: %s\n%s", message, format_help(parser))))
end

function show_help(parser::ArgumentParser)
    throw(ShowHelp(format_help(parser)))
end

function format_help(parser::ArgumentParser)
    keyword_args = map(values(parser.keyword_arguments)) do specification
        display = @sprintf(
            "%s%s",
            specification.argument_string,
            specification.action == :store ? @sprintf(" %s", specification.name) : "",
        )
        if !specification.required
            display = "[$display]"
        end
        display
    end
    positional_args = map(parser.positional_arguments) do specification
        if !specification.required
            "[$(specification.argument_string)]"
        else
            specification.argument_string
        end
    end
    display_parts = []
    if length(keyword_args) > 0
        append!(display_parts, keyword_args)
    end
    if length(positional_args) > 0
        append!(display_parts, positional_args)
    end
    lines = [
        @sprintf("Usage: %s %s", parser.program_name, join(display_parts, " "))
    ]
    join(lines, "\n")
end

export ArgumentParser, add_argument, parse_args

end # module ArgumentParse


text:
module TestData
    using Base.Test
    using DataArrays

    #test_context("Data types and NA's")

    #test_group("NA's")
    @assert length(NA) == 1
    @assert size(NA) == ()
    @assert isna(3 == NA)
    @assert isna(NA == 3)
    @assert isna(NA == NA)

    #test_group("DataVector creation")
    dvint = @data [1, 2, NA, 4]
    dvint2 = DataArray([5:8])
    dvint3 = convert(DataArray, 5:8)
    dvflt = @data [1.0, 2, NA, 4]
    dvstr = @data ["one", "two", NA, "four"]
    dvdict = DataArray(Dict, 4) # for issue #199

    @assert isa(dvint, DataVector{Int})
    @assert isa(dvint2, DataVector{Int})
    @assert isa(dvint3, DataVector{Int})
    @assert isa(dvflt, DataVector{Float64})
    @assert isa(dvstr, DataVector{ASCIIString})
    # @test throws_exception(DataArray([5:8], falses(2)), Exception)

    #test_group("PooledDataVector creation")
    pdvstr = @pdata ["one", "one", "two", "two", NA, "one", "one"]
    @assert isa(pdvstr, PooledDataVector{ASCIIString})
    # @test throws_exception(PooledDataVector["one", "one", 9], Exception)
    @assert isequal(PooledDataArray(pdvstr), pdvstr)

    #test_group("PooledDataVector creation with predetermined pool")
    pdvpp = PooledDataArray([1, 2, 2, 3], [1, 2, 3, 4])
    @assert isequal(pdvpp.pool, [1, 2, 3, 4])
    @assert string(pdvpp) == "[1, 2, 2, 3]"
    pdvpp = PooledDataArray([1, 2, 2, 3, 2, 1], [1, 2, 3, 4])
    @assert isequal(pdvpp.pool, [1, 2, 3, 4])
    @assert string(pdvpp) == "[1, 2, 2, 3, 2, 1]"
    pdvpp = PooledDataArray(["one", "two", "two"], ["one", "two", "three"])
    @assert isequal(convert(DataArray, pdvpp), @data(["one", "two", "two"]))
    @assert isequal(levels(pdvpp), @data(["one", "two", "three"]))
    @assert isequal(pdvpp.pool, ["one", "two", "three"])
    @assert string(pdvpp) == "[one, two, two]"
    @assert string(PooledDataArray(["one", "two", "four"],
                                   ["one", "two", "three"])) ==
            "[one, two, NA]"

    #test_group("PooledDataVector utf8 support")
    pdvpp = PooledDataArray([utf8("hello")], [false])
    @assert isa(pdvpp[1], UTF8String)
    pdvpp = PooledDataArray([utf8("hello")])
    @assert isa(pdvpp[1], UTF8String)

    #test_group("DataVector access")
    @assert dvint[1] == 1
    @assert isna(dvint[3])
    @assert isequal(dvflt[3:4], @data([NA, 4.0]))
    @assert isequal(dvint[[true, false, true, false]], @data([1, NA]))
    @assert isequal(dvstr[[1, 2, 1, 4]], @data(["one", "two", "one", "four"]))
    # Indexing produces #undef?
    # @assert isequal(dvstr[[1, 2, 1, 3]], DataVector["one", "two", "one", NA])

    #test_group("PooledDataVector access")
    @assert pdvstr[1] == "one"
    @assert isna(pdvstr[5])
    @assert isequal(pdvstr[1:3], @data(["one", "one", "two"]))
    @assert isequal(pdvstr[[true, false, true, false, true, false, true]],
                    @pdata(["one", "two", NA, "one"]))
    @assert isequal(pdvstr[[1, 3, 1, 2]], @data(["one", "two", "one", "one"]))

    #test_group("DataVector methods")
    @assert size(dvint) == (4,)
    @assert length(dvint) == 4
    @assert sum(isna(dvint)) == 1
    @assert eltype(dvint) == Int

    #test_group("PooledDataVector methods")
    @assert size(pdvstr) == (7,)
    @assert length(pdvstr) == 7
    @assert sum(isna(pdvstr)) == 1
    @assert eltype(pdvstr) == ASCIIString

    #test_group("DataVector operations")
    @assert isequal(dvint .+ 1, DataArray([2, 3, 4, 5], [false, false, true, false]))
    @assert isequal(dvint .* 2, @data([2, 4, NA, 8]))
    @assert isequal(dvint .== 2, @data([false, true, NA, false]))
    @assert isequal(dvint .> 1, @data([false, true, NA, true]))

    #test_group("PooledDataVector operations")
    # @assert isequal(pdvstr .== "two", PooledDataVector[false, false, true, true, NA, false, false])

    #test_group("DataVector to something else")
    @assert all(dropna(dvint) .== [1, 2, 4])
    @assert all(array(dvint, 0) .== [1, 2, 0, 4])
    utf8three = convert(UTF8String, "three")
    asciithree = convert(ASCIIString, "three")
    @assert all(array(dvstr, utf8three) .== ["one", "two", "three", "four"])
    @assert all(array(dvstr, asciithree) .== ["one", "two", "three", "four"])
    @assert all(convert(Vector{Int}, dvint2) .== [5:8])
    @assert all([i + 1 for i in dvint2] .== [6:9])
    @assert all([length(x)::Int for x in dvstr] == [3, 3, 1, 4])
    @assert repr(dvint) == "[1,2,NA,4]"

    #test_group("PooledDataVector to something else")
    @assert all(dropna(pdvstr) .== ["one", "one", "two", "two", "one", "one"])
    @assert all(array(pdvstr, "nine") .== ["one", "one", "two", "two", "nine", "one", "one"])
    @assert all([length(i)::Int for i in pdvstr] .== [3, 3, 3, 3, 1, 3, 3])
    @assert string(pdvstr[1:3]) == "[one, one, two]"

    #test_group("DataVector Filter and Replace")
    @assert isequal(dropna(dvint), [1, 2, 4])
    @assert isequal(array(dvint, 7), [1, 2, 7, 4])
    @assert sum(dropna(dvint)) == 7
    @assert sum(array(dvint, 7)) == 14

    #test_group("PooledDataVector Filter and Replace")
    @assert reduce(string, "", dropna(pdvstr)) == "oneonetwotwooneone"
    @assert reduce(string, "", array(pdvstr, "!")) == "oneonetwotwo!oneone"

    #test_group("DataVector assignment")
    assigntest = @data [1, 2, NA, 4]
    assigntest[1] = 8
    @assert isequal(assigntest, (@data [8, 2, NA, 4]))
    assigntest[1:2] = 9
    @assert isequal(assigntest, (@data [9, 9, NA, 4]))
    assigntest[[1,3]] = 10
    @assert isequal(assigntest, (@data [10, 9, 10, 4]))
    assigntest[[true, false, true, true]] = 11
    @assert isequal(assigntest, (@data [11, 9, 11, 11]))
    assigntest[1:2] = [12, 13]
    @assert isequal(assigntest, (@data [12, 13, 11, 11]))
    assigntest[[1, 4]] = [14, 15]
    @assert isequal(assigntest, (@data [14, 13, 11, 15]))
    assigntest[[true, false, true, false]] = [16, 17]
    @assert isequal(assigntest, (@data [16, 13, 17, 15]))
    assigntest[1] = NA
    @assert isequal(assigntest, (@data [NA, 13, 17, 15]))
    assigntest[[1, 2]] = NA
    @assert isequal(assigntest, (@data [NA, NA, 17, 15]))
    assigntest[[true, false, true, false]] = NA
    @assert isequal(assigntest, (@data [NA, NA, NA, 15]))
    assigntest[1] = 1
    assigntest[2:4] = NA
    @assert isequal(assigntest, (@data [1, NA, NA, NA]))

    #test_group("PooledDataVector assignment")
    ret = (pdvstr[2] = "three")
    @assert ret == "three"
    @assert pdvstr[2] == "three"
    ret = pdvstr[[1,2]] = "two"
    @assert ret == "two"
    @assert pdvstr[2] == "two"
    pdvstr2 = @pdata ["one", "one", "two", "two"]
    ret = (pdvstr2[[true, false, true, false]] = "three")
    @assert ret == "three"
    @assert pdvstr2[1] == "three"
    ret = (pdvstr2[[false, true, false, true]] = ["four", "five"])
    @assert isequal(ret, ["four", "five"])
    @assert isequal(pdvstr2[3:4], (@data ["three", "five"]))
    pdvstr2 = @pdata ["one", "one", "two", "two"]
    ret = (pdvstr2[2:3] = "three")
    @assert ret == "three"
    @assert isequal(pdvstr2[3:4], (@data ["three", "two"]))
    ret = (pdvstr2[2:3] = ["four", "five"])
    @assert ret == ["four", "five"]
    @assert isequal(pdvstr2[1:2], (@data ["one", "four"]))
    pdvstr2 = @pdata ["one", "one", "two", "two", "three"]
    @assert isna(begin pdvstr2[1] = NA end)
    @assert all(isna(begin pdvstr2[[1, 2]] = NA end))
    @assert all(isna(begin pdvstr2[[false, false, true, false, false]] = NA end))
    @assert all(isna(begin pdvstr2[4:5] = NA end))
    @assert all(isna(pdvstr2))

    #test_group("PooledDataVector replace!")
    pdvstr2 = @pdata ["one", "one", "two", "two", "three"]
    @assert replace!(pdvstr2, "two", "four") == "four"
    @assert replace!(pdvstr2, "three", "four") == "four"
    @assert isna(replace!(pdvstr2, "one", NA))
    @assert replace!(pdvstr2, NA, "five") == "five"
    @assert isequal(pdvstr2, (@data ["five", "five", "four", "four", "four"]))
end


text:
@testset "MPS-based search produces correctly results for small instances" for L ∈[2, 3]
    N = L^2
    instance = "$(@__DIR__)/instances/basic/$(N)_001.txt"
    ig = ising_graph(instance)
    rank = Tuple(values(get_prop(ig, :rank)))

    β = 2.
    dβ = β / 5.
    schedule = fill(dβ, Int(ceil(β/dβ)))

    ϵ = 1E-8
    Dcut = 32
    var_ϵ = 1E-8
    max_sweeps = 4
    states = all_states(rank)
    ϱ = gibbs_tensor(ig, β)

    @testset "Sqrt of the Gibbs state (aka state tensor)" begin
        ψ = ones(rank...)
        for σ ∈ states
            for i ∈ 1:N
                h = get_prop(ig, i, :h)
                nbrs = unique_neighbors(ig, i)
                ψ[idx.(σ)...] *= exp(-0.5 * β * h * σ[i])
                for j ∈ nbrs
                    J = get_prop(ig, i, j, :J)
                    ψ[idx.(σ)...] *= exp(-0.5 * β * σ[i] * J * σ[j])
                end
            end
        end
        ρ = abs.(ψ) .^ 2
        rψ = MPS(ψ)
        lψ = MPS(ψ, :left)

        @testset "produces correct Gibbs state" begin
            @test ρ / sum(ρ) ≈ ϱ
        end
        @testset "MPS from the tensor" begin
            @testset "can be right normalized" begin
                @test dot(rψ, rψ) ≈ 1
                @test_nowarn is_right_normalized(rψ)
            end
            @testset "can be left normalized" begin
                @test dot(lψ, lψ) ≈ 1
                @test_nowarn is_left_normalized(lψ)
            end
            @testset "both forms are the same (up to a phase factor)" begin
                vlψ = vec(tensor(lψ))
                vrψ = vec(tensor(rψ))
                vψ = vec(ψ)
                vψ /= norm(vψ)
                @test abs(1 - abs(dot(vlψ, vrψ))) < ϵ
                @test abs(1 - abs(dot(vlψ, vψ))) < ϵ
            end
        end
        @testset "MPS from gates" begin
            Gψ = MPS(ig, Dcut, var_ϵ, max_sweeps, schedule)
            @testset "is built correctly" begin
                @test abs(1 - abs(dot(Gψ, rψ))) < ϵ
            end
            @testset "is normalized" begin
                @test dot(Gψ, Gψ) ≈ 1
                @test_nowarn is_right_normalized(Gψ)
            end
            @testset "has correct links and non-trivial bond dimension" begin
                @test bond_dimension(Gψ) > 1
                @test_nowarn verify_bonds(Gψ)
            end
        end
        @testset "Exact probabilities are calculated correctely" begin
            for σ ∈ states
                p, r = dot(rψ, σ), dot(rψ, proj(σ, rank), rψ)
                @test p ≈ r
                @test ϱ[idx.(σ)...] ≈ p
            end
        end
        @testset "Results from MPS-based search agree with brute-force" begin
            for max_states ∈ [1, N, 2*N, 3*N, N^2-3, N^2-2, N^2-1, N^2]
                states, _ = solve(rψ, max_states)
                en = energy.(states, Ref(ig))
                sp = brute_force(ig, num_states = max_states)
                for (e, ẽ) ∈ zip(sort(en), sp.energies)
                    @test e ≈ ẽ
                end
            end
        end
    end
end


text:
#                                                  ....
#                                 ......''''''....;cll:;,..
#                          .';cllc:::cccc:c::::cccldxdlll:,.
#                         .ldxxxxkdlcccc:ccccccccclodoccllc.
#                       .':oddkxxlxkdollccccccc::::cl:';c:'.
#                       ,ccdxkkxd;dOxddoclllcllllc:::c,....
#                      .;llcclc:;:xxxddooxkkddddol:::::;;;'.
#                    ..,;cloollooddddddddxxkO0Odollcc:;;,'.
#                 .,;lddl:cloooddddddddxkOOOOOkxdolc:;,'..
#          .;:,. .:loxOOdlccoodxxxxxxxdddxxxxddolc;,'...
#        .l0K0ko:coodOK0xol;;,;::ccc:;;;;;;;;,''........    ..:dxl,
#       .l0XX0xolodddkOkdolc:::,'.......................  ..,:oO00Oc
#        :kOkxoloxxkkkdc;;;;;;,,''.........................':lodxkx:.
#        .odooolxxxkkkd,;clcccc;,'.......................',:loool:'
#      .;:ododxkxdkkkkx;;olllccc:;,,,,,,,,'''''''...;c:;;;:clc;.
#      ,ldxxxkkkddxkxxxoclclclcc:::::::cc:::::;,'..cxkoc:;;;.
#      'ldxxxxkdlddxkxkOkxxo::c:::::;;;;;::cl:,'.'ckOOdc;'.....
#      .lddddddxxdkkOxxO0kdxOkl,'.......,:lodooloodkkxo:;:ccllcl:;.
#      .cdoc:oxkOOkkOkxk0kxxdxkdo:.....,dkkxxxdoxxxxxdc;;lddxdddol;.
#       'ld:..:lxdxOOkxkkxxl'okxdc.....'lxkdlcdkOOxk0Ol,.:oddddxdcc.
#        ,:.  .;kOOkkOxxkkx:..,;,.......',,,'ckO0kxO0k:...',oddxocl;..  ........
#            .cxkxolc;:xkkd:................ckkkkxxxdl,....'lddolco:'...........
#            .lololc..ckkxxl'..............cxxddxdddl:'....'ldol:lo:..........
#             .,,;'. .oxxxxd:............,lxdoc;,coo:'.....'ldc;loc,........
#                     .,:::'..   ..  . .'oddol;':ooc;.....,cc:.'ll:'..    ..
#                                       .,;:;,..;c:;'. . 'cl;. .....
#                                                 .      ...
# Glowe.jl - Julia interface to GloVe
#            Written in Julia at 0x0α Research
#			 by Corneliu Cofaru, 2018
#
# Original authors: Jeffrey Pennington, Richard Socher
#					and Christopher D. Manning. 2014

module Glowe

    using LinearAlgebra
    using Statistics
    import Base: show, size

    export vocab_count, cooccur, shuffle, glove,
           WordVectors, wordvectors,
           index, size,
           get_vector,
           vocabulary, in_vocabulary,
           cosine, cosine_similar_words,
           analogy, analogy_words,
           similarity

    include("interface.jl")
    include("wordvectors.jl")

end # module


text:
using Pkg

Pkg.test("ITensors")


text:
using ModelConstructors, Test, Random

@testset "Prior computation" begin
    p = ParameterVector{Float64}(undef, 3)
    p[1] = ModelConstructors.parameter(:a, 0.1, (0., 1.), (0., 1.), ModelConstructors.Untransformed(), Normal(0.5, 1.); fixed = false)
    p[2] = ModelConstructors.parameter(:b, 0., (0., 1.), (0., 1.), ModelConstructors.Untransformed(), Normal(0.5, 1.); fixed = true)
    p[3] = ModelConstructors.parameter(:c, 0.4, (0., 1.), (0., 1.), ModelConstructors.Untransformed(), Normal(0.25, 1.); fixed = false)

    @test ModelConstructors.prior(p) == logpdf(p[1]) + logpdf(p[3])
end

@testset "Posterior computation" begin
    p = ParameterVector{Float64}(undef, 3)
    p[1] = ModelConstructors.parameter(:a, 0.1, (0., 1.), (0., 1.), ModelConstructors.Untransformed(), Normal(0.5, 1.); fixed = false)
    p[2] = ModelConstructors.parameter(:b, 0., (0., 1.), (0., 1.), ModelConstructors.Untransformed(), Normal(0.5, 1.); fixed = true)
    p[3] = ModelConstructors.parameter(:c, 0.4, (0., 1.), (0., 1.), ModelConstructors.Untransformed(), Normal(0.25, 1.); fixed = false)

    loglh = (x, data) -> 2.
    data = rand(2,2)

    @test posterior(loglh, p, data, ϕ_smc = 3.104) == 2. * 3.104 + ModelConstructors.prior(p)

    oldvals = map(x -> x.value, p)
    oldpost = 2. * 3.104 + ModelConstructors.prior(p)
    @test posterior!(loglh, p, oldvals, data, ϕ_smc = 3.104) == oldpost

    newvals = copy(oldvals)
    newvals[1] += 0.4 # moves value to the mean
    newvals[3] -= 0.15  # moves value to the mean
    @test posterior!(loglh, p, newvals, data, ϕ_smc = 3.104) > oldpost
end

@testset "Prior computation, updating, and sampling with regime-switching" begin
    p = ParameterVector{Float64}(undef, 3)
    p[1] = ModelConstructors.parameter(:a, 0.1, (0., 1.), (0., 1.), ModelConstructors.Untransformed(), Normal(0.5, 1.); fixed = false)
    p[2] = ModelConstructors.parameter(:b, 0., (0., 1.), (0., 1.), ModelConstructors.Untransformed(), Normal(0.5, 1.); fixed = true)
    p[3] = ModelConstructors.parameter(:c, 0.4, (0., 1.), (0., 1.), ModelConstructors.Untransformed(), Normal(0.25, 1.); fixed = false)
    p3copy = ModelConstructors.parameter(:c, 0.5, (0., 1.), (0., 1.), ModelConstructors.Untransformed(), Normal(0.25, 1.); fixed = false)

    # Set up parameter switching
    set_regime_val!(p[1], 1, 0.1) # test prior switching
    set_regime_val!(p[1], 2, 0.2)
    set_regime_val!(p[1], 3, 0.3)
    set_regime_prior!(p[1], 1, Normal(0.5, 1.))
    set_regime_prior!(p[1], 2, Normal(0.2, 1.))
    set_regime_prior!(p[1], 3, Normal(0.2, 1.))
    set_regime_val!(p[2], 1, 0.1; override_bounds = true)
    set_regime_val!(p[2], 2, 0.2; override_bounds = true)
    set_regime_fixed!(p[2], 1, true)
    set_regime_fixed!(p[2], 2, true)
    @test set_regime_val!(p[2], 2, 0.2) == 0.2 # trying to reset a fixed value won't do anything
    @test regime_val(p[2], 2) == 0.2
    set_regime_val!(p[3], 1, 0.4)
    set_regime_val!(p[3], 2, 0.5; override_bounds = true)
    set_regime_fixed!(p[3], 1, true)
    set_regime_fixed!(p[3], 2, false)

    @test !ModelConstructors._filter_all_fixed_para(p[1])
    @test ModelConstructors._filter_all_fixed_para(p[2])
    @test !ModelConstructors._filter_all_fixed_para(p[3])
    set_regime_fixed!(p[3], 2, true)
    @test ModelConstructors._filter_all_fixed_para(p[3])
    set_regime_fixed!(p[3], 2, false)
    set_regime_valuebounds!(p[3], 2, (0., 1.)) # set_regime_fixed!(p[3], 2, true) will change the valuebounds

    @test prior(p) == sum([logpdf(p[1]), logpdf(p3copy)])

    p_in = [0.2, 0.1, 0.5, 0.3, 0.9, 0.2, 0.6]
    toggle_regime!(p, 1)
    update!(p, p_in)
    p_out = [0.2, 0.1, 0.4, 0.3, 0.9, 0.2, 0.6]
    @test ModelConstructors.get_values(p) == p_out
    @test regime_val(p[1], 1) == 0.2
    @test regime_val(p[1], 2) == 0.3
    @test regime_val(p[1], 3) == 0.9

    Random.seed!(1793)
    out1 = ModelConstructors.rand_regime_switching(p; toggle = true)
    Random.seed!(1793)
    for para in p
        toggle_regime!(para, 1)
    end
    out2 = rand(p; regime_switching = true, toggle = false)
    update!(p, out1)

    @test out1 == out2
    @test regime_val(p[1], 1) == out1[1]
    @test regime_val(p[1], 2) == out1[4]
    @test regime_val(p[1], 3) == out1[5]
    @test regime_val(p[2], 1) == 0.1
    @test regime_val(p[2], 2) == 0.2
    @test regime_val(p[3], 1) == 0.4
    @test regime_val(p[3], 2) == out1[end]
end

nothing

max_stars_repo_name: FRBNY-DSGE/ModelConstructors

max_issues_repo_name: FRBNY-DSGE/ModelConstructors

max_forks_repo_name: FRBNY-DSGE/ModelConstructors


text:
# Tests whether points project correctly

print("\tTest 1 - #generators == dim == 2...")
cone_1 = FiniteCone([[0.0, 1.0] [1.0, 1.0]])
test_data_1 = [([-1.0, 3.0], [0.0, 3.0]),
               ([-1.0, -1.0], [0.0, 0.0]),
               ([2.0, 0.0], [1.0, 1.0]),
               ([0.25, 1.0], [0.25, 1.0])]
for (x, proj_x) in test_data_1
    @test project(cone_1, x) ≈ proj_x
end
print("success!\n")

print("\tTest 2 - 3 = #generators > dim = 2...")
cone_2 = FiniteCone([[0.0, 1.0] [1.0, 1.0] [0.5, 1.0]])
test_data_2 = [([-1.0, 3.0], [0.0, 3.0]),
               ([-1.0, -1.0], [0.0, 0.0]),
               ([2.0, 0.0], [1.0, 1.0]),
               ([0.25, 1.0], [0.25, 1.0])]
for (x, proj_x) in test_data_2
    @test project(cone_2, x) ≈ proj_x
end
print("success!\n")

print("\tTest 3 - 2 = #generators < dim = 3...")
cone_3 = FiniteCone([[0.0, 0.0, 1.0] [1.0, 1.0, 1.0]])
test_data_3 = [([-1.0, -1.0, 1.0], [0.0, 0.0, 1.0]),
               ([3.0, 0.0, 0.0], [1.0, 1.0, 1.0]),
               ([-1.0, -1.0, -1.0], [0.0, 0.0, 0.0]),
               ([1.0, 2.0, 3.0], [1.5, 1.5, 3.0]),
               ([2.5, 2.5, 6.0], [2.5, 2.5, 6.0])]
for (x, proj_x) in test_data_3
    @test project(cone_3, x) ≈ proj_x
end
print("success!\n")


text:
#
# Combinator API.
#

immutable Combinator
    use::Function
end

Query(C::Combinator) =
    cse(format(C.use(ItQuery(Void) |> setnamespace(DB))))

>>(q::Query, C::Combinator) =
    C.use(q)

(C::Combinator)() = C

convert(::Type{Combinator}, q::Query) =
    Combinator() do q′::Query
        (q′ >> q) |> setnamespace(q′.ns.db)
    end

execute{T}(C::Combinator, val::T=nothing; paramvals...) =
    execute(C.use(ItQuery(T)), val; paramvals...)

# Attribute access.

Field(tag::Symbol) =
    Combinator() do q::Query
        scope = oscope(q)
        maybe_binding = lookup(scope, tag, 0)
        @assert !isnull(maybe_binding) "cannot find attribute $tag in the scope of $(scope.dom)"
        binding = get(maybe_binding)
        (q >> compile(binding, scope, Syntax[])) |> setnamespace(scope.db)
    end

Fields(tag::Symbol, tags′::Symbol...) =
    (Field(tag), (Field(tag′) for tag′ in tags′)...)

convert(::Type{Combinator}, tag::Symbol) =
    Field(tag)

# Composition.

>>(L::Combinator, R::Combinator) =
    Combinator() do q::Query
        (q >> L) >> R
    end

# Null.

const Null =
    Combinator() do q::Query
        q >> NullQuery(domain(output(q)))
    end

# Identity.

const It =
    Combinator() do q::Query
        q
    end

# Decorations.

Decorate(F::Combinator, decor::Pair{Symbol}) =
    Combinator() do q::Query
        (q >> F) |> decorate(decor)
    end

ThenDecorate(decor::Pair{Symbol}) =
    Combinator() do q::Query
        q |> decorate(decor)
    end

ThenTag(tag::Symbol) =
    ThenDecorate(:tag => tag)

convert{T}(::Type{Combinator}, pair::Pair{Symbol,T}) =
    let (tag, F) = pair
        convert(Combinator, F) >> ThenTag(tag)
    end

# Constants.

Const(val) =
    Combinator() do q::Query
        q >> ConstQuery(val)
    end

convert(::Type{Combinator}, val::Union{Void,Bool,Signed,Unsigned,String}) =
    Const(val)

# Scalar functions.

Lift(fn, argtypes::Tuple{Vararg{Type}}, restype::Type, F::Combinator, Gs...) =
    Combinator() do q::Query
        it = ostub(q)
        q >> LiftQuery(fn, argtypes, restype, it >> F, (it >> convert(Combinator, G) for G in Gs)...)
    end

Lift(fn, argtypes::Tuple{Vararg{Type}}, F::Combinator, Gs...) =
    Combinator() do q::Query
        it = ostub(q)
        q >> LiftQuery(fn, argtypes, it >> F, (it >> convert(Combinator, G) for G in Gs)...)
    end

Lift(fn, F::Combinator, Gs...) =
    Combinator() do q::Query
        it = ostub(q)
        q >> LiftQuery(fn, it >> F, (it >> convert(Combinator, G) for G in Gs)...)
    end

.==(F::Combinator, G) = Lift(==, F, G)
.!=(F::Combinator, G) = Lift(!=, F, G)
.<(F::Combinator, G) = Lift(<, F, G)
.<=(F::Combinator, G) = Lift(<=, F, G)
.>(F::Combinator, G) = Lift(>, F, G)
.>=(F::Combinator, G) = Lift(>=, F, G)

(~)(F::Combinator) = Lift(~, F)
(&)(F::Combinator, Gs::Combinator...) = Lift(&, F, Gs...)
(|)(F::Combinator, Gs::Combinator...) = Lift(|, F, Gs...)

.+(F::Combinator, Gs...) = Lift(+, F, Gs...)
.-(F::Combinator, Gs...) = Lift(-, F, Gs...)
.*(F::Combinator, Gs...) = Lift(*, F, Gs...)
./(F::Combinator, Gs...) = Lift(/, F, Gs...)
.÷(F::Combinator, Gs...) = Lift(÷, F, Gs...)

# Aggregates.

Count(F::Combinator) =
    Combinator() do q::Query
        q >> CountQuery(ostub(q) >> F)
    end

const ThenCount =
    Combinator() do q::Query
        CountQuery(q) |> setnamespace(q.ns.db)
    end

Exists(F::Combinator) =
    Combinator() do q::Query
        q >> ExistsQuery(ostub(q) >> F)
    end

const ThenExists =
    Combinator() do q::Query
        ExistsQuery(q) |> setnamespace(q.ns.db)
    end

Aggregate(fn, argtype::Type, restype::Type, haszero::Bool, F::Combinator) =
    Combinator() do q::Query
        q >> AggregateQuery(fn, argtype, restype, haszero, ostub(q) >> F)
    end

Aggregate(fn, argtype::Type, restype::Type, F::Combinator) =
    Combinator() do q::Query
        q >> AggregateQuery(fn, argtype, restype, ostub(q) >> F)
    end

Aggregate(fn, argtype::Type, haszero::Bool, F::Combinator) =
    Combinator() do q::Query
        q >> AggregateQuery(fn, argtype, haszero, ostub(q) >> F)
    end

Aggregate(fn, argtype::Type, F::Combinator) =
    Combinator() do q::Query
        q >> AggregateQuery(fn, argtype, ostub(q) >> F)
    end

Aggregate(fn, haszero::Bool, F::Combinator) =
    Combinator() do q::Query
        q >> AggregateQuery(fn, haszero, ostub(q) >> F)
    end

Aggregate(fn, F::Combinator) =
    Combinator() do q::Query
        q >> AggregateQuery(fn, ostub(q) >> F)
    end

Aggregate(fn, argtype::Type, restype::Type, haszero::Bool) =
    F -> Aggregate(fn, argtype, restype, haszero, convert(Combinator, F))

Aggregate(fn, argtype::Type, restype::Type) =
    F -> Aggregate(fn, argtype, restype, convert(Combinator, F))

Aggregate(fn, argtype::Type, haszero::Bool) =
    F -> Aggregate(fn, argtype, haszero, convert(Combinator, F))

Aggregate(fn, argtype::Type) =
    F -> Aggregate(fn, argtype, convert(Combinator, F))

Aggregate(fn, haszero::Bool) =
    F -> Aggregate(fn, haszero, convert(Combinator, F))

Aggregate(fn) =
    F -> Aggregate(fn, convert(Combinator, F))

ThenAggregate(fn, argtype::Type, restype::Type, haszero::Bool) =
    Combinator() do q::Query
        AggregateQuery(fn, argtype, restype, haszero, q) |> setnamespace(q.ns.db)
    end

ThenAggregate(fn, argtype::Type, restype::Type) =
    Combinator() do q::Query
        AggregateQuery(fn, argtype, restype, q) |> setnamespace(q.ns.db)
    end

ThenAggregate(fn, argtype::Type, haszero::Bool) =
    Combinator() do q::Query
        AggregateQuery(fn, argtype, haszero, q) |> setnamespace(q.ns.db)
    end

ThenAggregate(fn, argtype::Type) =
    Combinator() do q::Query
        AggregateQuery(fn, argtype, q) |> setnamespace(q.ns.db)
    end

ThenAggregate(fn, haszero::Bool) =
    Combinator() do q::Query
        AggregateQuery(fn, haszero, q) |> setnamespace(q.ns.db)
    end

ThenAggregate(fn) =
    Combinator() do q::Query
        AggregateQuery(fn, q) |> setnamespace(q.ns.db)
    end

const AnyOf = Aggregate(any)
const AllOf = Aggregate(all)
const MaxOf = Aggregate(maximum, false)
const MinOf = Aggregate(minimum, false)
const SumOf = Aggregate(sum)
const MeanOf = Aggregate(mean)

const ThenAny = ThenAggregate(any)
const ThenAll = ThenAggregate(all)
const ThenMax = ThenAggregate(maximum, false)
const ThenMin = ThenAggregate(minimum, false)
const ThenSum = ThenAggregate(sum)
const ThenMean = ThenAggregate(mean)

# Filtering.

ThenFilter(P::Combinator) =
    Combinator() do q::Query
        FilterQuery(q, ostub(q) >> P) |> setnamespace(q.ns)
    end

# Sorting.

Sort(F::Combinator) =
    Combinator() do q::Query
        (q >> SortQuery(ostub(q) >> F)) |> setnamespace(q.ns)
    end

ThenSort(Fs::Combinator...) =
    Combinator() do q::Query
        it = ostub(q)
        SortQuery(q, (it >> F for F in Fs)...) |> setnamespace(q.ns)
    end

const Asc = ThenDecorate(:rev => false)
const Desc = ThenDecorate(:rev => true)

# Paginating.

ThenTake(N::Combinator) =
    Combinator() do q::Query
        TakeQuery(q, istub(q) >> N) |> setnamespace(q.ns)
    end

ThenSkip(N::Combinator) =
    Combinator() do q::Query
        SkipQuery(q, istub(q) >> N) |> setnamespace(q.ns)
    end

ThenTake(N) = ThenTake(convert(Combinator, N))
ThenSkip(N) = ThenSkip(convert(Combinator, N))

# Output.

ThenSelect(Fs::Combinator...) =
    Combinator() do q::Query
        it = ostub(q)
        compile(StaticBinding{:select}(), oscope(it), q, Query[it >> F for F in Fs])
    end

ThenSelect(Fs...) = ThenSelect((convert(Combinator, F) for F in Fs)...)

ThenSelect(;Fs...) = ThenSelect(Fs...)

ThenSummarize(Fs::Combinator...) =
    Combinator() do q::Query
        it = istub(q)
        compile(StaticBinding{:select}(), oscope(it), it, Query[q >> F for F in Fs])
    end

ThenSummarize(Fs...) = ThenSummarize((convert(Combinator, F) for F in Fs)...)

ThenSummarize(;Fs...) = ThenSummarize(Fs...)

# Hierarchical queries.

Connect(F::Combinator) =
    Combinator() do q::Query
        it = ostub(q)
        (q >> ConnectQuery(false, it >> F)) |> setnamespace(q.ns.db)
    end

# Grouping and roll-up.

ThenGroup(Fs::Combinator...) =
    Combinator() do q::Query
        it = ostub(q)
        compile(StaticBinding{:group}(), oscope(it), q, Query[it >> F for F in Fs])
    end

ThenGroup(Fs...) = ThenGroup((convert(Combinator, F) for F in Fs)...)

ThenGroup(;Fs...) = ThenGroup(Fs...)

ThenRollUp(Fs::Combinator...) =
    Combinator() do q::Query
        it = ostub(q)
        compile(StaticBinding{:rollup}(), oscope(it), q, Query[it >> F for F in Fs])
    end

ThenRollUp(Fs...) = ThenRollUp((convert(Combinator, F) for F in Fs)...)

ThenRollUp(;Fs...) = ThenRollUp(Fs...)

Unique(F::Combinator) =
    Combinator() do q::Query
        it = ostub(q)
        (q >> UniqueQuery(it >> F)) |> setnamespace(q.ns.db)
    end

const ThenUnique =
    Combinator() do q::Query
        UniqueQuery(q)
    end

# Parameters.

Parameter(tag::Symbol, oty) =
    Combinator() do q::Query
        (q >> SlotQuery(tag, oty)) |> setnamespace(q.ns.db)
    end

Given(Fs::Combinator...) =
    Combinator() do q::Query
        it = istub(q)
        GivenQuery(q, (it >> F for F in Fs)...) |> setnamespace(q.ns.db)
    end

Given(Fs...) = Given((convert(Combinator, F) for F in Fs)...)

Given(;Fs...) = Given(Fs...)

function Using(args::Combinator...)
    Q = args[end]
    Fs = args[1:end-1]
    Combinator() do q::Query
        it = ostub(q)
        r = it >> Q
        q >> GivenQuery(r, (it >> F for F in Fs)...) |> setnamespace(r.ns.db)
    end
end

Using(Fs...) = Using((convert(Combinator, F) for F in Fs)...)

# Before and around.

Before(Fs::Combinator...) =
    Combinator() do q::Query
        it = ostub(q)
        (q >> AroundQuery(domain(output(q)), true, true, false, (it >> F for F in Fs)...)) |> setnamespace(q.ns.db)
    end

Around(Fs::Combinator...) =
    Combinator() do q::Query
        it = ostub(q)
        (q >> AroundQuery(domain(output(q)), true, true, true, (it >> F for F in Fs)...)) |> setnamespace(q.ns.db)
    end

const ThenFrame =
    Combinator() do q::Query
        FrameQuery(q) |> setnamespace(q.ns)
    end

Frame(F::Combinator) =
    Combinator() do q::Query
        it = ostub(q)
        r = it >> F
        q >> FrameQuery(r) |> setnamespace(r.ns)
    end

# Cardinality assertions.

ExpectOne(F::Combinator) =
    Combinator() do q::Query
        it = ostub(q)
        r = it >> F
        q >> ExpectOneQuery(r) |> setnamespace(r.ns.db)
    end

ExpectAtMostOne(F::Combinator) =
    Combinator() do q::Query
        it = ostub(q)
        r = it >> F
        q >> ExpectAtMostOneQuery(r) |> setnamespace(r.ns.db)
    end

ExpectAtLeastOne(F::Combinator) =
    Combinator() do q::Query
        it = ostub(q)
        r = it >> F
        q >> ExpectAtLeastOneQuery(r) |> setnamespace(r.ns.db)
    end

const ThenExpectOne =
    Combinator() do q::Query
        ExpectOneQuery(q) |> setnamespace(q.ns.db)
    end

const ThenExpectAtMostOne =
    Combinator() do q::Query
        ExpectAtMostOneQuery(q) |> setnamespace(q.ns.db)
    end

const ThenExpectAtLeastOne =
    Combinator() do q::Query
        ExpectAtLeastOneQuery(q) |> setnamespace(q.ns.db)
    end

# Merging.

Merge(Fs::Combinator...) =
    Combinator() do q::Query
        it = ostub(q)
        q >> MergeQuery(Query[it >> F for F in Fs]) |> setnamespace(q.ns.db)
    end

# Multiplexing.

ThenElse(F::Combinator, G::Combinator) =
    Combinator() do q::Query
        it = ostub(q)
        IfThenElseQuery(q, it >> F, it >> G)
    end

Switch(T::Combinator, F::Combinator, G::Combinator) =
    Combinator() do q::Query
        it = ostub(q)
        t = it >> T
        f = it >> F
        g = it >> G
        q >> IfThenElseQuery(t, f, g)
    end

Switch(T::Combinator, F::Combinator, G::Combinator, H::Combinator, Ks::Combinator...) =
    Combinator() do q::Query
        it = ostub(q)
        t = it >> T
        f = it >> F
        g = it >> Switch(G, H, Ks...)
        q >> IfThenElseQuery(t, f, g)
    end

ExpectAnythingQuery(base::Query) = ExpectQuery(base, true, true)

ExpectAnything(F::Combinator) =
  Combinator() do q::Query
     it = ostub(q)
     q >> ExpectAnythingQuery(it >> F)
  end

FieldIfExists(tag::Symbol) =
    Combinator() do q::Query
        scope = oscope(q)
        maybe_binding = lookup(scope, tag, 0)
        if isnull(maybe_binding)
            return q
        end
        binding = get(maybe_binding)
        (q >> compile(binding, scope, Syntax[])) |> setnamespace(scope.db)
    end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
using BenchlingIntegration
using Test

@testset "BenchlingIntegration.jl" begin
    # Was Aristotle right?
	@test getA() == :A
end


text:
export GammaShapeRate

import Distributions: Gamma, shape, rate
import SpecialFunctions: loggamma, digamma, gamma
import StatsFuns: log2π

struct GammaShapeRate{T <: Real} <: ContinuousUnivariateDistribution
    a :: T
    b :: T
end

GammaShapeRate(a::Real, b::Real)       = GammaShapeRate(promote(a, b)...)
GammaShapeRate(a::Integer, b::Integer) = GammaShapeRate(float(a), float(b))
GammaShapeRate(a::Real)                = GammaShapeRate(a, one(a))
GammaShapeRate()                       = GammaShapeRate(1.0, 1.0)

Distributions.@distr_support GammaShapeRate 0 Inf

Distributions.support(dist::GammaShapeRate) = Distributions.RealInterval(minimum(dist), maximum(dist))

Distributions.shape(dist::GammaShapeRate)  = dist.a
Distributions.rate(dist::GammaShapeRate)   = dist.b
Distributions.scale(dist::GammaShapeRate)  = inv(dist.b)
Distributions.mean(dist::GammaShapeRate)   = shape(dist) / rate(dist)
Distributions.var(dist::GammaShapeRate)    = shape(dist) / abs2(rate(dist))
Distributions.params(dist::GammaShapeRate) = (shape(dist), rate(dist))

function Distributions.entropy(dist::GammaShapeRate)
    a, b = params(dist)
    return a - log(b) + loggamma(a) + (1-a)*digamma(a)
end

function logmean(dist::GammaShapeRate)
    a, b = params(dist)
    return digamma(a) - log(b)
end

function loggammamean(dist::GammaShapeRate)
    a, b = params(dist)
    return 0.5 * (log2π - (digamma(a) - log(b))) + mean(dist) * (-1 + digamma(a + 1) - log(b))
end

function meanlogmean(dist::GammaShapeRate)
    a, b = params(dist)
    return mean(dist) * (digamma(a + 1) - log(b))
end

Base.eltype(::GammaShapeRate{T}) where T = T

Base.convert(::Type{ GammaShapeRate{T} }, a::Real, b::Real) where { T <: Real } = GammaShapeRate(convert(T, a), convert(T, b))

vague(::Type{ <: GammaShapeRate }) = GammaShapeRate(1.0, tiny)

prod_analytical_rule(::Type{ <: GammaShapeRate }, ::Type{ <: GammaShapeRate }) = ProdAnalyticalRuleAvailable()

function prod(::ProdAnalytical, left::GammaShapeRate, right::GammaShapeRate)
    T = promote_type(eltype(left), eltype(right))
    return GammaShapeRate(shape(left) + shape(right) - one(T), rate(left) + rate(right))
end

Distributions.pdf(dist::GammaShapeRate, x::Real)    = (rate(dist)^shape(dist)) / gamma(shape(dist)) * x^(shape(dist)-1)*exp(-rate(dist) * x)
Distributions.logpdf(dist::GammaShapeRate, x::Real) = shape(dist) * log(rate(dist)) - loggamma(shape(dist)) + (shape(dist)-1)*log(x) - rate(dist)*x


text:
using BenchmarkTools
using Random

f(x) = 3x.^2 + 4x + 7x.^3
fdot(x) = 3 .* x.^2 .+ 4 .* x .+ 7 .* x.^3

function main()

    x = randn(10000)

    # @benchmark f($x)
    # @benchmark fdot($x)
    @benchmark f.($x)
end

main()

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
K(c::Float64, cost::CostParameters) = cost.k1/(1+cost.k2*c)
Φ(vN, vI) = maximum([vN vI])
r(ci::Float64, cj::Float64) = ci >= cj ? 0.0 : cj - ci
r1(c1, c2) = r(c1::Float64, c2::Float64)
r2(c1, c2) = r(c2::Float64, c1::Float64)
number_of_strings(nC::Int64) = div(nC*(nC+1)*(2*nC+1),6)

# Functions for the lowest layer

H1(ic1, ic2, pi, τ, ess)     = pi*Φ(τ[ess.nC].EQs[ic1, ic2, ess[ic1,ic2,ess.nC]].vN1,τ[ess.nC].EQs[ic1, ic2, ess[ic1,ic2,ess.nC]].vI1)
H2(ic1, ic2, pi, τ, ess)     = pi*Φ(τ[ess.nC].EQs[ic1, ic2, ess[ic1,ic2,ess.nC]].vN2,τ[ess.nC].EQs[ic1, ic2, ess[ic1,ic2,ess.nC]].vI2)

# Functions for the upper layers

H1(ic1, ic2, ic, pi, τ, ess) = (1-pi)*Φ(τ[ic+1].EQs[ic1, ic2, ess[ic1,ic2,ic+1]].vN1,τ[ic+1].EQs[ic1, ic2, ess[ic1,ic2,ic+1]].vI1) +
                               pi*Φ(τ[ic].EQs[ic1, ic2, ess[ic1,ic2,ic]].vN1,τ[ic].EQs[ic1, ic2, ess[ic1,ic2,ic]].vI1)
H2(ic1, ic2, ic, pi, τ, ess) = (1-pi)*Φ(τ[ic+1].EQs[ic1, ic2, ess[ic1,ic2,ic+1]].vN2,τ[ic+1].EQs[ic1, ic2, ess[ic1,ic2,ic+1]].vI2) +
                               pi*Φ(τ[ic].EQs[ic1, ic2, ess[ic1,ic2,ic]].vN2,τ[ic].EQs[ic1, ic2, ess[ic1,ic2,ic]].vI2)
# Quadratic solver

function quad(a::Float64, b::Float64, c::Float64)
    d = b^2 - 4*a*c
    if abs(a) < 1e-8
        pstar = [0. ; 1.;-c/b]
    else
        if d < 0
            pstar = [0. ;1.]
        elseif d == 0.
            pstar = [0. ;1. ; -b/2a]
        else
            pstar = [0. ;1. ; (-b - sqrt(d))/(2a); (-b + sqrt(d))/(2a)]
        end
    end
    return pstar
end

function jump!(x::ESS)
    activated = false
    if x[1,1,1] == -1
        throw(error("This ESS has already overflown!"))
    end
    X = Array(Int64,length(x))
    R = 1 # Remainder; first "remainder" is one, since we add one
    @inbounds for i in reverse(1:length(x)) # Slightly inefficient, as we go through R = 0 cases
        X[i] = mod(x.number[i] + R,x.bases[i])
        R = div(x.number[i]+R,x.bases[i])
        if R == 0 && activated == false
            x.active = i
            activated = true
        end
    end
    if R > 0
        println("No more equilibria to check.")
        x.number = -1*ones(Int64,length(x))
    else
        x.number = X
    end
end