text:
@doc raw"""
    Plane3

An object `h` of the data type [`Plane3`](@ref) is an oriented plane in the
three-dimensional Euclidean space ``Ε^3``.

It is defined by the set of points with Cartesian coordinates ``(x,y,z)`` that
satisfy the plane equation

```math
    h:\; a\, x + b\, y + c\, z + d = 0.
```

The plane splits ``Ε^3`` in a *positive* and a *negative side*. A point `p` with
Cartesian coordinates ``(px, py, pz)`` is on the positive side of `h`, iff
``a\, px + b\, py + c\, pz + d > 0``. It is on the negative side, iff
``a\, px + b\, py + c\, pz + d < 0``.
"""
Plane3

@doc raw"""
    Plane3(a::Real, b::Real, c::Real, d::Real)

Create a plane `h` defined by the equation ``a\, px + b\, py + c\, pz + d = 0``.

Notice that `h` is degenerate if ``a = b = c = 0``.
"""
Plane3(a::Real, b::Real, c::Real, d::Real) =
    Plane3(convert(FT, a), convert(FT, b),
           convert(FT, c), convert(FT, d))

"""
    Plane3(p::Point3, q::Point3, r::Point3)

Create a plane `h` passing through the points `p`, `q`, and `r`.

The plane is oriented such that `p`, `q` and `r` are oriented in a positive
sense (that is counterclockwise) when seen form the positive side of `h`. Notice
that `h` is degenerate if the points are collinear.
"""
Plane3(p::Point3, q::Point3, r::Point3)

"""
    Plane3(p::Point3, v::Vector3)

Introduces a plane `h` that passes through point `p` and that is orthogonal to
`v`.
"""
Plane3(p::Point3, v::Vector3)

"""
    Plane3(p::Point3, d::Direction3)

Introduces a plane `h` that passes through point `p` and that has an orthogonal
direction equal to `d`.
"""
Plane3(p::Point3, d::Direction3)

"""
    Plane3(l::Line3, p::Point3)

Introduces a plane `h` that is defined through the three points `point(l, 0)`,
`point(l, 1)` and `p`.
"""
Plane3(l::Line3, p::Point3)

"""
    Plane3(r::Ray3, p::Point3)

Introduces a plane `h` that is defined through the three points `point(r, 0)`,
`point(r, 1)` and `p`.
"""
Plane3(r::Ray3, p::Point3)

"""
    Plane3(s::Segment3, p::Point3)

Introduces a plane `h` that is defined through the three points `source(s)`,
`target(s)` and `p`.
"""
Plane3(s::Segment3, p::Point3)

"""
    Plane3(c::Circle3)

Introduces a plane `h` that is defined as the plane containing the circle `c`.
"""
Plane3(c::Circle3)

"""
    ==(h₁::Plane3, h₂::Plane3)

Test for equality: two planes are equal, iff they have a non empty intersection
and the same orientation.
"""
==(h₁::Plane3, h₂::Plane3)

"""
    a(h::Plane3)

Returns the first coefficient of `h`.
"""
a(h::Plane3)

"""
    b(h::Plane3)

Returns the second coefficient of `h`.
"""
b(h::Plane3)

"""
    c(h::Plane3)

Returns the third coefficient of `h`.
"""
c(h::Plane3)

"""
    d(h::Plane3)

Returns the fourth coefficient of `h`.
"""
d(h::Plane3)

"""
    perpendicular_line(h::Plane3, p::Point3)

Returns the line that is perpendicular to `h` and that passes through point `p`.

The line is oriented from the negative to the positive side of `h`.
"""
perpendicular_line(h::Plane3, p::Point3)

"""
    projection(h::Plane3, p::Point3)

Returns the orthogonal projection of `p` on `h`.
"""
projection(h::Plane3, p::Point3)

"""
    opposite(h::Plane3)

Returns the plane with opposite orientation.
"""
opposite(h::Plane3)

"""
    point(h::Plane3)

Returns an arbitrary point on `h`.
"""
point(h::Plane3)

"""
    orthogonal_vector(h::Plane3)

Returns a vector that is orthogonal to `h` and that is directed to the positive
side of `h`.
"""
orthogonal_vector(h::Plane3)

"""
    orthogonal_direction(h::Plane3)

Returns a direction that is orthogonal to `h` and that is directed to the positive
side of `h`.
"""
orthogonal_direction(h::Plane3)

"""
    base1(h::Plane3)

Returns a vector orthogonal to [`orthogonal_vector()`](@ref).
"""
base1(h::Plane3)

"""
    base2(h::Plane3)

Returns a vector that is both orthogonal to [`base1()`](@ref), and to
[`orthogonal_vector()`](@ref), and such that the result of `orientation(
point(h), point(h) + base1(h), point(h) +
base2(h), point(h) + orthogonal_vector(h) )` is positive.
"""
base2(h::Plane3)

"""
    to_2d(h::Plane3, p::Point3)

Returns the image point of the projection of `p` under an affine transformation,
which maps `h` onto the ``xy``-plane, with the ``z``-coordinate removed.
"""
to_2d(h::Plane3, p::Point3)

"""
    to_3d(h::Plane3, p::Point2)

Returns a point `q`, such that `to_2d(h, to_3d(h, p))` is equal to `p`.
"""
to_3d(h::Plane3, p::Point2)

"""
    oriented_side(h::Plane3, p::Point3)

Returns either [`ON_ORIENTED_BOUNDARY`](@ref), or the constant
[`ON_POSITIVE_SIDE`](@ref), or the constant [`ON_NEGATIVE_SIDE`](@ref),
determined by the position of `p` relative to the oriented plane `h`.
"""
oriented_side(h::Plane3, p::Point3)

"""
    has_on(h::Plane3, p::Point3)

Returns `true`, iff `p` properly lies on `h`.
"""
has_on(h::Plane3, p::Point3)

"""
    has_on_positive_side(h::Plane3, p::Point3)

Returns `true`, iff `p` lies on the positive side of `h`.
"""
has_on_positive_side(h::Plane3, p::Point3)

"""
    has_on_negative_side(h::Plane3, p::Point3)

Returns `true`, iff `p` lies on the negative side of `h`.
"""
has_on_negative_side(h::Plane3, p::Point3)

"""
    has_on(h::Plane3, l::Line3)

Returns `true`, iff `l` properly lies on `h`.
"""
has_on(h::Plane3, l::Line3)

"""
    has_on(h::Plane3, c::Circle3)

Returns `true`, iff `c` properly lies on `h`.
"""
has_on(h::Plane3, c::Circle3)

"""
    is_degenerate(h::Plane3)

Plane `h` is degenerate, if the coefficients `a`, `b`, and `c` of the plane
equation are zero.
"""
is_degenerate(h::Plane3)

"""
    transform(h::Plane3, t::AffTransformation3)

Returns the plane obtained by applying `t` on a point of `h` and the
orthogonal direction of `h`.
"""
transform(h::Plane3, t::AffTransformation3)


text:
using HDF5, JLD

function savedata(savepath, psfst, skyst, algost, admmst, toolst)

    JLD.save("imagerec.jld","admmst.x",admmst.x)

    JLD.save(savepath,"algost.lastiter",algost.lastiter,
                      "algost.nfreq",algost.nfreq,
                      "toolst.err",toolst.err,
                      "toolst.tol1",toolst.tol1,
                      "toolst.tol2",toolst.tol2,
                      "toolst.tol3",toolst.tol3,
                      "toolst.tol4",toolst.tol4,
                      "toolst.tol5",toolst.tol5,
                      "skyst.mydata",skyst.mydata,
                      "skyst.sky",skyst.sky,
                      "psfst.nu",psfst.nu,
                      "psfst.nu0",psfst.nu0,
                      "admmst",admmst)

                    #   "skyst.noise",skyst.noise)
end

#   "toolst.snr",toolst.snr,
#   "toolst.errorrec",toolst.errorrec,
#   "toolst.errorest",toolst.errorest,
#   "toolst.errorraw",toolst.errorraw,
# "admmst.x",admmst.x
# "admmst.t",admmst.t,
# "admmst.taut",admmst.taut,
# "admmst.v",admmst.v,
# "admmst.tauv",admmst.tauv,
# "admmst.s",admmst.s,
# "admmst.taus",admmst.taus,
# "admmst.p",admmst.p,
# "admmst.taup",admmst.taup

ext: jl


text:
using Test
using Distributions
using KernelDensity

import KernelDensity: kernel_dist, default_bandwidth, kde_boundary, kde_range, tabulate

for D in [Tuple{Normal,Normal}, Tuple{Uniform,Uniform}, Tuple{Logistic,Logistic},
          (Normal, Normal), (Uniform, Uniform), (Logistic, Logistic)]
    d = KernelDensity.kernel_dist(D,(0.5,0.5))
    dx,dy = d
    @test mean(dx) == 0.0
    @test mean(dy) == 0.0
    @test std(dx) ≈ 0.5
    @test std(dy) ≈ 0.5
end

r = kde_range((-2.0,2.0), 128)
@test step(r) > 0

for X in ([0.0], [0.0,0.0], [0.0,0.5], [-0.5:0.1:0.5;])
    w = default_bandwidth(X)
    @test w > 0
    lo, hi = kde_boundary(X,w)
    @test lo < hi
    kr = kde_range((lo,hi), 10)
    @test step(kr) > 0

    for D in (Normal, )
        k1 = tabulate((X,X),(r,r))
        @test isa(k1,BivariateKDE)
        @test size(k1.density) == (length(k1.x), length(k1.y))
        @test all(k1.density .>= 0.0)
        @test sum(k1.density)*step(k1.x)*step(k1.y) ≈ 1.0

        k2 = KernelDensity.conv(k1,kernel_dist(Tuple{D,D}, (0.1,0.1)))
        @test isa(k2,BivariateKDE)
        @test size(k2.density) == (length(k2.x), length(k2.y))
        @test all(k2.density .>= 0.0)
        @test sum(k2.density)*step(k2.x)*step(k2.y) ≈ 1.0

        k3 = kde((X,X);kernel=D)
        @test isa(k3,BivariateKDE)
        @test size(k3.density) == (length(k3.x), length(k3.y))
        @test all(k3.density .>= 0.0)
        @test sum(k3.density)*step(k3.x)*step(k3.y) ≈ 1.0

        k4 = kde((X,X),(r,r);kernel=D)
        @test isa(k4,BivariateKDE)
        @test size(k4.density) == (length(k4.x), length(k4.y))
        @test all(k4.density .>= 0.0)
        @test sum(k4.density)*step(k4.x)*step(k4.y) ≈ 1.0

        k5 = kde([X X];kernel=D)
        @test isa(k5,BivariateKDE)
        @test size(k5.density) == (length(k5.x), length(k5.y))
        @test all(k5.density .>= 0.0)
        @test sum(k5.density)*step(k5.x)*step(k5.y) ≈ 1.0

        k6 = kde([X X],(r,r);kernel=D, weights=fill(1.0/length(X),length(X)))
        @test k4.density ≈ k6.density
    end
end

k11 = kde([0.0 0.0; 1.0 1.0], (r,r), bandwidth=(1,1), weights=[0,1])
k12 = kde([1.0 1.0], (r,r), bandwidth=(1,1))
@test k11.density ≈ k12.density


text:
export  excelPopPlot

function excelPopPlot( mpSim::MPsim, filename::AbstractString,
    sheetname::AbstractString="Node population plots" )

    filename = string( filename, endswith( filename, ".xlsx" ) ? "" : ".xlsx" )

    if !ispath( filename )
        @warn string( "File '", filename,
            "' not found. Can't generate plots/reports." )
        return
    end  # if !ispath( filename )

    plotsToMake = Dict{Tuple, Vector{String}}()

    XLSX.openxlsx( filename ) do xf
        if !XLSX.hassheet( xf, sheetname )
            @warn string( "Sheet '", sheetname,
                "' not found in file. Can't generate plots/reports." )
            return
        end  # if !XLSX.hassheet( xf, sheetname )

        xs = xf[sheetname]
        nToProcess = xs["B4"]

        for ii in (1:nToProcess) .+ 6
            # The name of the node to plot.
            node = xs[ii, 1]
            node = node isa Missing ? "active" : node

            # The time resolution.
            timeRes = xs[ii, 2]

            if !( timeRes isa Real ) || ( timeRes <= 0 )
                continue
            end  # if !( timeRes isa Real ) || ...

            # The type of the flux plots.
            fluxType = lowercase( string( xs[ii, 3] ) )
            fluxType = fluxType == "line" ? :normal :
                ( fluxType == "stacked" ? :stacked : :none )

            # The type of the composition plots.
            compType = :none

            if haskey( mpSim.compoundNodeList, node )
                if fluxType === :none
                    compType = lowercase( string( xs[ii, 4] ) )
                    compType = compType == "line" ? :normal :
                        ( compType == "stacked" ? :stacked : :none )
                    else
                    compType = fluxType
                end  # if fluxType === :none
            end  # if haskey( mpSim.compoundNodeList, node )

            # Show plot in browser.
            showPlot = lowercase( string( xs[ii, 5] ) ) == "yes"

            # Save plot.
            savePlot = lowercase( string( xs[ii, 6] ) ) == "yes"

            # Generate report.
            makeReport = lowercase( string( xs[ii, 7] ) ) == "yes"

            plotkey = (timeRes, fluxType, compType, showPlot, savePlot,
                makeReport)

            if haskey( plotsToMake, plotkey )
                push!( plotsToMake[plotkey], node )
            else
                plotsToMake[plotkey] = [node]
            end  # if haskey( plotsToMake, plotkey )
        end  # for ii in 1:nToProcess
    end  # XLSX.openxlsx( filename ) do xf

    plotdir = filename[1:(end-5)]

    if !ispath( plotdir )
        mkpath( plotdir )
    end  # if !ispath( plotdir )

    dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
    isFirst = true
    reportname = string( plotdir, "/node report ", dateStr )

    for plotkey in keys( plotsToMake )
        if !any( [plotkey[4], plotkey[5], plotkey[6]] )
            continue
        end  # if !any( [plotkey[4], plotkey[5], plotkey[6]] )

        dateStr = string( " (", Dates.format( now(), "yyyymmdd HHMMSS" ), ";",
            plotkey[1], ")" )
        makePlot = plotkey[4] || plotkey[5]
        makeReport = plotkey[6]

        if plotkey[2] === :none
            # Generate the report.
            tStart = now()
            popReport = nodePopReport( mpSim, plotkey[1],
                plotsToMake[plotkey]... )
            reportGenerationTime = ( now() - tStart ).value / 1000.0

            if isempty( popReport )
                continue
            end  # if isempty( popReport )

            plotname = string( plotdir, "/node population", dateStr )

            if makePlot
                generatePopPlot( popReport, 12.0, plotkey[4], plotkey[5],
                    plotname )
            end  # makePlot

            if makeReport
                writeNodePopData( popReport, mpSim, 12.0, reportGenerationTime,
                    reportname, isFirst )
                isFirst = false
            end  # if makeReport

            if plotkey[3] !== :none
                nodes = filter( node -> haskey( mpSim.compoundNodeList, node ),
                    plotsToMake[plotkey] )

                # Generate the report.
                tStart = now()
                compReport = nodeCompositionReport( mpSim, plotkey[1],
                    nodes... )
                reportGenerationTime = ( now() - tStart ).value / 1000.0

                if isempty( compReport )
                    continue
                end  # if isempty( compReport )

                plotFunction = plotTypes[plotkey[3]][2]
                plotname = string( plotdir, "/node composition", dateStr,
                    ".svg" )

                if makePlot
                    generateCompPlot( compReport, nodes, 12.0, plotkey[4],
                        plotkey[5], plotname, plotFunction )
                end  # if makePlot

                if makeReport
                    writeNodeCompositionData( compReport, nodes, mpSim, 12.0,
                        reportGenerationTime, reportname, false )
                end  # if makeReport
            end  # if plotkey[3] !== :none
        else
            # Generate the report.
            tStart = now()
            popReport, fluxReports = nodeEvolutionReport( mpSim, plotkey[1],
                plotsToMake[plotkey]... )
            reportGenerationTime = ( now() - tStart ).value / 1000.0

            if isempty( popReport ) || isempty( fluxReports )
                return
            end  # if isempty( popReport ) || ...

            plotFunction = plotTypes[plotkey[2]][1]
            plotFunction2 = plotTypes[plotkey[2]][2]
            plotname = string( plotdir, "/node evolution", dateStr, ".svg" )

            if makePlot
                generateEvolutionPlot( popReport, mpSim, fluxReports,
                    Tuple( plotsToMake[plotkey] ), 12.0, plotkey[4], plotkey[5],
                    plotname, plotFunction, plotFunction2 )
            end  # if makePlot

            if makeReport
                writeNodeEvolutionData( (popReport, fluxReports), mpSim, 12.0,
                    reportGenerationTime, reportname, isFirst )
                isFirst = false
            end  # if makeReport
        end  # if plotkey[2] === :none
    end  # for plotkey in keys( plotsToMake )

end  # excelPopPlot( mpSim, filename )


text:
"""
    LaTeXSVG(latex::AbstractString, svg::String; standalone::Bool=false)

This type contains the LaTeX code to be rendered, the preamble, and the rendered SVG string. Both [`latexsvg`](@ref) and [`@Lsvg_str`](@ref) return objects of this type.

`LaTeXSVG` objects contain 4 fields:
- `latex` contains the LaTeX input.
- `standalone` indicates whether `latex` is a complete LaTeX document. If `standalone=true` then `pre` is empty.
- `pre` is a vector of strings containing the preamble, which is recorded from the global preamble when a `LaTeXSVG` object is contructed, namely when `latexsvg` or `@Lsvg_str` is called.
- `svg` contains the SVG string.

You can access these fields with the usual dot syntax.

!!! note
    You should never need to contruct a `LaTeXSVG` object yourself.

!!! note
    If you are using VSCode, by default the color of the SVG displayed in the plot pane adapts to the color scheme of VSCode (black for light mode and white for dark mode.) Rest assured that this adaptive color is not written into the SVG itself. Also, if you defined custom colors in your LaTeX, it should be displayed (and saved to file) faithfully.
"""
struct LaTeXSVG
    latex::AbstractString
    standalone::Bool
    pre::Vector{String}
    svg::String
    function LaTeXSVG(latex::AbstractString, svg::String; standalone::Bool=false)
        if standalone
            return new(latex, true, String[], svg)
        else
            return new(latex, false, deepcopy(preamble()), svg)
        end
    end
end

function _shorten(line::AbstractString, maxlength::Int)
    if length(line) <= maxlength
        return line
    else
        return line[begin:(begin+maxlength-4)] * "..."
    end
end

function Base.print(io::IO, svg::LaTeXSVG)
    print(io, "LaTeXSVG\n\n    SVG: ")

    delim = Sys.iswindows() ? "\r\n" : "\n"
    svg_split = split(svg.svg, delim, limit=8)
    line_length = length(svg_split[begin + 1])
    print(io, svg_split[begin], "\n")
    for line in svg_split[begin + 1:end - 1]
        print(io, "         ", _shorten(line, line_length), "\n")
    end
    print(io, "         .........\n")

    print(io, "\n    LaTeX: ")
    latex_split = split(svg.latex, "\n")
    print(io, latex_split[begin], "\n")
    for line in latex_split[begin + 1:end]
        print(io, "           $(line)\n")
    end

    if !svg.standalone
        print(io, "\n    preamble: $(svg.pre[begin])")
        for pre in svg.pre[begin + 1:end]
            print(io, "\n              " * pre)
        end
    end
end

Base.show(io::IO, ::MIME"text/plain", svg::LaTeXSVG) = print(io, svg)
Base.show(io::IO, ::MIME"image/svg+xml", svg::LaTeXSVG) = write(io, svg.svg)


text:
using FLoops
using Folds
using FoldsKernelAbstractions
using Referenceables
using Test

function inc!(xs, ex = nothing)
    @floop ex for x in referenceable(xs)
        x[] += 1
    end
    return xs
end

function test_size(
    device;
    groupsize,
    basesize = 2,
    extraitems = 0,
    n = groupsize * basesize + extraitems,
)
    ex = KAEx(device; groupsize, basesize)

    @test Folds.sum(1:n, ex) == sum(1:n)
    @test Folds.reduce(xor, 1:n, ex; init = 0) == reduce(xor, 1:n; init = 0)

    xs = FoldsKernelAbstractions.arrayfor(device)(1:n)
    @test collect(inc!(xs, ex)) == 2:n+1
end

function test_sweep_sizes(device)
    @testset for groupsize in [1, 2, 8, 64, 256],
        basesize in [2, 5, 16, 19],
        extraitems in [0, 1, 3]

        test_size(device; groupsize, basesize, extraitems)
    end
    @testset "$kwds" for kwds in [
        (groupsize = 4, n = 3),
        (groupsize = 4, n = 2),
        (groupsize = 4, n = 1),
        (groupsize = 4, basesize = 10, n = 30),
        (groupsize = 4, basesize = 10, n = 20),
        (groupsize = 4, basesize = 10, n = 10),
    ]
        test_size(device; kwds...)
    end
end

function test_all(device)
    @testset "sweep sizes" begin
        test_sweep_sizes(device)
    end
end


text:
struct SFixedGainDivHandler end

gethandler(::Bool, ::Type{typeof(÷)}, ::Type{<:SBitstreamLike}, ::Type{<:Real}) = SFixedGainDivHandler()
init_state(::SFixedGainDivHandler) = (id = 0,)

function (handler::SFixedGainDivHandler)(buffer, netlist, state, inputs, outputs)
    # update netlist with inputs
    setsigned!(netlist, inputs[1], true)

    # compute output size
    outsize = netsize(outputs[1])

    # add output net to netlist
    setsigned!(netlist, outputs[1], true)

    write(buffer, """
        $stdcomment
        // BEGIN fdiv$(state.id)
        stoch_fixed_gain_div_mat #(
                .GAIN($(name(inputs[2]))),
                .NUM_ROWS($(outsize[1])),
                .NUM_COLS($(outsize[2]))
            ) fdiv$(state.id)_p (
                .CLK(CLK),
                .nRST(nRST),
                .A($(name(inputs[1]))_p),
                .Y($(name(outputs[1]))_p)
            );
        stoch_fixed_gain_div_mat #(
                .GAIN($(name(inputs[2]))),
                .NUM_ROWS($(outsize[1])),
                .NUM_COLS($(outsize[2]))
            ) fdiv$(state.id)_m (
                .CLK(CLK),
                .nRST(nRST),
                .A($(name(inputs[1]))_m),
                .Y($(name(outputs[1]))_m)
            );
        // END fdiv$(state.id)
        \n""")

    return buffer, (id = state.id + 1,)
end


text:
struct Point{T}
    x::T
    y::T
    z::Union{T, Nothing}
    name::String
end

Point(x::T, y::T, name::String) where {T} = Point(x, y, nothing, name)

function moveto(sp, point::Point; safe_height=55, wait=false)
    println("Moving to $(point.name)")
    write(sp, "G1 Z$(safe_height)\n")
    write(sp, "G1 X$(point.x) Y$(point.y)\n")
    (!isnothing(point.z)) && write(sp, "G1 Z$(point.z)\n")

    wait && wait_for_key("press any key to continue")
end

function create_wells(topleft; well_spacing = 9.0, well_depth = 11.0)
    well_names = map(x->"$(x[1])$(x[2])", Iterators.product('A':'H', 1:12))

    wells = Array{Point, 2}(undef, size(well_names))

    for I in CartesianIndices(well_names)
        wells[I] = Point(
            plate.x + well_spacing * (I[2] - 1),
            plate.y - well_spacing * (I[1] - 1),
            plate.z - well_depth,
            well_names[I]
        )
    end

    return wells
end

wait_for_key(prompt) = (print(stdout, prompt); read(stdin, 1); nothing)

max_stars_repo_name: tlnagy/fractionation

max_issues_repo_name: tlnagy/fractionation

max_forks_repo_name: tlnagy/fractionation


text:
using Distributions
using MultidimensionalTables

module T

using MultidimensionalTables
using DataFrames

println("test 1")

begin
  lar = @larr(a=repmat(1:100,10000), b=rand(1000000))
  #lar = @larr(a=repmat(1:1000,10000), b=rand(10000000))
  df = DataFrame(lar)

  @time r1 = @select(lar, by[:a], m=length(_))
  @time r2 = by(df, :a) do df
    DataFrame(m=length(df))
  end
  @time r1 = @select(lar, by[:a], m=sum(_b))
  @time r2 = by(df, :a) do df
    DataFrame(m=sum(df[:a]))
  end
end

println("test filter")
for n in [10_000_000, 100_000_000][1:1]
  for m in [100, 10_000, 1_000_000][1:1]
    d = darr(x=sample(1:m, n, replace=true), y=rand(n))
    dm = darr(x=sample(1:m, m))
    @show @time @select(d, where[_x .>= 10], where[_x .< 20])
    @show @time @select(d, where[_x .>= Nullable(10)], where[_x .< Nullable(20)])
  end
end

println("test sort")
for n in [10_000_000, 100_000_000][1:1]
  for m in [100, 10_000, 1_000_000][1:1]
    d = darr(x=sample(1:m, n, replace=true), y=rand(n))
    dm = darr(x=sample(1:m, m))
    @show @time sort(d, 1, :x)
    @show @time sort(d, 1, :y)
  end
end

println("test new column")
for n in [10_000_000, 100_000_000][1:1]
  for m in [100, 10_000, 1_000_000][1:1]
    d = darr(x=sample(1:m, n, replace=true), y=rand(n))
    dm = darr(x=sample(1:m, m))
    @show @time @update(d, y2=2*_y)
  end
end

println("test aggregation")
for n in [10_000_000, 100_000_000][1:1]
  for m in [100, 10_000, 1_000_000][1:1]
    d = darr(x=sample(1:m, n, replace=true), y=rand(n))
    dm = darr(x=sample(1:m, m))
    @show @time @select(d, by[:x], ym=mean(_y))
  end
end

println("test join")
for n in [10_000_000, 100_000_000][1:1]
  for m in [100, 10_000, 1_000_000][1:1]
    d = darr(x=sample(1:m, n, replace=true), y=rand(n))
    dm = larr(axis1=darr(x=sample(1:m, m, replace=false)), n=1:m)
    @show @time innerjoin(d, dm, 1)
    @show @time leftjoin(d, dm, 1)
  end
end
#@show(@allocated r1 = @select(lar, by[:a], m=length(_)))
#@show(@allocated r2 = by(df, :a) do df
#  DataFrame(m=length(df))
#end)
#@show(@allocated r1 = @select(lar, by[:a], m=sum(_b)))
#@show(@allocated r2 = by(df, :a) do df
#  DataFrame(m=sum(df[:a]))
#end)

end


text:
@testset "policies" begin
    @testset "DeterministicPolicy" begin
        p = DeterministicPolicy([1,2,1], 2)
        @test p(1) == 1
        @test p(2) == 2
        @test p(3) == 1

        @test p(1, Val(:dist)) == [1.0, 0.0]
        @test p(2, Val(:dist)) == [0.0, 1.0]
        @test p(3, Val(:dist)) == [1.0, 0.0]

        @test p(3, 1) == 1.0
        @test p(3, 2) == 0.0

        update!(p, 3, 2)

        @test p(3) == 2
        @test p(3, Val(:dist)) == [0.0, 1.0]
        @test p(3, 1) == 0.0
        @test p(3, 2) == 1.0
    end

    @testset "EpsilonGreedyPolicy" begin
        ϵ, nactions, best_actions = 0.1, 2, [1, 2, 1]
        p = EpsilonGreedyPolicy(best_actions, nactions, ϵ)
        seed!(123)
        for (s, a) in enumerate(best_actions)
            @test isapprox(sum((i) -> p(s) == a ? 1 : 0, 1:1000) / 1000, 1 - ϵ + ϵ / nactions, atol=0.05)
        end

        @test p(1, Val(:dist)) ≈ [0.95, 0.05]
        @test p(2, Val(:dist)) ≈ [0.05, 0.95]
        @test p(3, Val(:dist)) ≈ [0.95, 0.05]

        @test p(1, 1) ≈ 0.95
        @test p(1, 2) ≈ 0.05

        update!(p, 1, 2)
        @test p(1, Val(:dist)) ≈ [0.05, 0.95]
        @test p(1, 1) ≈ 0.05
        @test p(1, 2) ≈ 0.95
    end

    @testset "FunctionalPolicy" begin
        p = FunctionalPolicy(s -> 1)
        @test p(nothing) == 1
    end

    @testset "LinearPolicy" begin
        s1 = [1. 0. 1.; 0. 1. 0.]
        s2 = [0. 1. 1.; 1. 1. 1.]
        features = zeros(2, 2, 3)
        features[1, :, :] .= s1
        features[2, :, :] .= s2
        weights = [1., 1., 1.]
        p = LinearPolicy(features, weights)

        seed!(123)
        action_counts = countmap(p(1) for _ in 1:1000)
        prob = [ℯ^2/(ℯ^2+ℯ), ℯ/(ℯ^2+ℯ)]

        @test isapprox([action_counts[1]/1000, action_counts[2]/1000], prob, atol=0.1)
        @test isapprox(p(1, 1), prob[1])
        @test isapprox(p(1, Val(:dist)), prob)

        update!(p, -1.)
        @test isapprox(p(1, Val(:dist)), [0.5, 0.5])
    end

    @testset "RandomPolicy" begin
        prob = [0.3, 0.7]
        p = RandomPolicy(prob)

        seed!(123)
        action_counts = countmap(p(1) for _ in 1:1000)
        @test isapprox([action_counts[1]/1000, action_counts[2]/1000], prob, atol=0.1)
        @test p(1, Val(:dist)) == prob
        @test p(1, 1) == prob[1]
        @test p(1, 2) == prob[2]

        prob = [0.3  0.7; 0.7 0.3]
        p = RandomPolicy(prob)

        seed!(123)
        action_counts = countmap(p(1) for _ in 1:1000)
        @test isapprox([action_counts[1]/1000, action_counts[2]/1000], prob[1, :], atol=0.1)
        action_counts = countmap(p(2) for _ in 1:1000)
        @test isapprox([action_counts[1]/1000, action_counts[2]/1000], prob[2, :], atol=0.1)

        @test p(1, Val(:dist)) == prob[1, :]
        @test p(2, Val(:dist)) == prob[2, :]
        @test p(1, 1) == prob[1, 1]
        @test p(1, 2) == prob[1, 2]
        @test p(2, 1) == prob[2, 1]
        @test p(2, 2) == prob[2, 2]
    end
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 867b290c-6bd9-11eb-1afb-47757ee13936
begin
	using PlutoUI
	using NativeSVG
end

# ╔═╡ 66ddc134-6bea-11eb-3bde-0fea30c863af
let
	t = Turtle()
	include("../src/chap05.jl")
	penup(t)
	forward(t, -200)
	turn(t, 90)
	forward(t, 60)
	turn(t,-90)
	pendown(t)
	koch(t, 400)
	Drawing(t, 720, 130)
end

# ╔═╡ 4cc6b8ac-6bd9-11eb-1d23-f5d9caf02cd3
md"""# Conditionals and Recursion

The main topic of this chapter is the `if` statement, which executes different code depending on the state of the program. But first I want to introduce two new operators: floor division and modulus."""

# ╔═╡ 81ea78c0-6bd9-11eb-1d3b-5f67fdf090f4
md"""## Floor Division and Modulus

The *floor division* operator, `÷` (**`\div TAB`**), divides two numbers and rounds down to an integer. For example, suppose the running time of a movie is 105 minutes. You might want to know how long that is in hours. Conventional division returns a floating-point number:

```julia
julia> minutes = 105
105
julia> minutes / 60
1.75
```

But we don’t normally write hours with decimal points. Floor division returns the integer number of hours, rounding down:

```julia
julia> hours = minutes ÷ 60
1
```

To get the remainder, you could subtract one hour in minutes:

```julia
julia> remainder = minutes - hours * 60
45
```

An alternative is to use the *modulus operator*, `%`, which divides two numbers and returns the remainder:

```julia
julia> remainder = minutes % 60
45
```

!!! tip
    The modulus operator is more useful than it seems. For example, you can check whether one number is divisible by another—if `x % y` is `0`, then `x` is divisible by `y`.

    Also, you can extract the rightmost digit or digits from a number. For example, `x % 10` yields the rightmost digit of an integer `x` (in base 10). Similarly, `x % 100` yields the last two digits.
"""

# ╔═╡ fa7e737c-6bd9-11eb-3cd0-d5b1d2ae7af1
md"""## Boolean Expressions

A *Boolean expression* is an expression that is either true or false. The following examples use the operator `==`, which compares two operands and produces `true` if they are equal and `false` otherwise:

```julia
julia> 5 == 5
true
julia> 5 == 6
false
```

`true` and `false` are special values that belong to the type `Bool`; they are not strings:

```julia
julia> typeof(true)
Bool
julia> typeof(false)
Bool
```

The `==` operator is one of the relational operators (operators that compare their operands). The others are:

```julia
x!=y # x is not equal to y
x≠y  # (\ne TAB)
x>y  # x is greater than y
x<y  # x is less than y
x>=y # x is greater than or equal to y x≥y # (\ge TAB)
x<=y # x is less than or equal to y
x≤y  # (\le TAB)
```

!!! danger
    Although these operations are probably familiar to you, the Julia symbols are different from the mathematical symbols. A common error is to use a single equals sign (`=`) instead of a double equals sign (`==`). Remember that `=` is an assignment operator and `==` is a relational operator. There is no such thing as `=<` or `=>`.
"""

# ╔═╡ 01a13616-6bdb-11eb-1b8b-8bdda36e623f
md"""## Logical Operators

There are three *logical operators*: `&&` (and), `||` (or), and `!` (not). The semantics (meaning) of these operators is similar to their meaning in English. For example, `x > 0 && x < 10` is true only if `x` is greater than `0` *and* less than `10`.

`n % 2 == 0 || n % 3 == 0` is true if *either or both* of the conditions is true; that is, if the number is divisible by 2 or 3.

Both `&&` and `||` associate to the right (i.e., grouped from the right), but `&&` has higher precedence than `||` does.

Finally, the `!` operator negates a Boolean expression, so `!(x > y)` is true if `x > y` is false; that is, if `x` is less than or equal to `y`.
"""

# ╔═╡ 7478967a-6bdb-11eb-03f7-81f894751f6d
md"""## Conditional Execution

In order to write useful programs, we almost always need the ability to check conditions and change the behavior of the program accordingly. *Conditional statements* give us this ability. The simplest form is the `if` statement:

```julia
if x > 0
	println("x is positive")
end
```

The Boolean expression after `if` is called the *condition*. If it is true, the indented statement runs. If not, nothing happens.

`if` statements have the same structure as function definitions: a header followed by a body terminated with the keyword `end`. Statements like this are called *compound statements*.

There is no limit on the number of statements that can appear in the body. Occasionally, it is useful to have a body with no statements (usually as a placeholder for code you haven’t written yet):

```julia
if x < 0
	# TODO: need to handle negative values!
end
```
"""

# ╔═╡ e75638d2-6bdb-11eb-0a7e-93446ea41c2d
md"""## Alternative Execution

A second form of the `if` statement is “alternative execution,” in which there are two possibilities and the condition determines which one runs. The syntax looks like this:

```julia
if x % 2 == 0
	println("x is even")
else
	println("x is odd")
end
```

If the remainder when `x `is divided by 2 is 0, then we know that `x` is even, and the program displays an appropriate message. If the condition is false, the second set of statements runs. Since the condition must be true or false, exactly one of the alternatives will run. The alternatives are called *branches*, because they are branches in the flow of execution.
"""

# ╔═╡ 482d8390-6bdc-11eb-3e3c-33b9f4f244b9
md"""## Chained Conditionals

Sometimes there are more than two possibilities and we need more than two branches. One way to express a computation like that is using a *chained conditional*:

```julia
if x < y
	println("x is less than y")
elseif x > y
	println("x is greater than y")
else
	println("x and y are equal")
end
```

Again, exactly one branch will run. There is no limit on the number of `elseif` statements. If there is an `else` clause, it has to be at the end, but there doesn’t have to be one:

```julia
if choice == "a"
	draw_a()
elseif choice == "b"
	draw_b()
elseif choice == "c"
	draw_c()
end
```

Each condition is checked in order. If the first is false, the next is checked, and so on. If one of them is true, the corresponding branch runs and the statement ends. Even if more than one condition is true, only the first true branch runs.
"""

# ╔═╡ 9a3080ec-6be4-11eb-32eb-f162987e0fd6
md"""
!!! languages
    In Python, `if` statements have the same structure as function definitions: a header followed by an indented body.

    ```python
    if x < y:
        print('x is less than y')
    elif x > y:
        print('x is greater than y')
    else:
        print('x and y are equal')

    ```

    In the C programming language, we find a similar `if` statement:

    ```c
    if (x < y) {
        printf("x is less than y\n");
    }
    else if (x > y) {
        printf("x is greater than y\n");
    }
    else {
        printf("x and y are equal\n");
    }
    ```
"""

# ╔═╡ c0f31858-6bdc-11eb-00f8-9b4fb5bebb3e
md"""## Nested Conditionals

One conditional can also be nested within another. We could have written the example in the previous section like this:

```julia
if x == y
	println("x and y are equal")
else
	if x < y
		println("x is less than y")
	else
		println("x is greater than y")
	end
end
```

The outer conditional contains two branches. The first branch contains a simple statement. The second branch contains another `if` statement, which has two branches of its own. Those two branches are both simple statements, although they could have been conditional statements as well.

Although the (noncompulsory) indentation of the statements makes the structure apparent, *nested conditionals* become difficult to read very quickly. It is a good idea to avoid them when you can.

Logical operators often provide a way to simplify nested conditional statements. For example, we can rewrite the following code using a single conditional:

```julia
if 0 < x
	if x < 10
		println("x is a positive single-digit number.")
	end
end
```

The print statement runs only if we make it past both conditionals, so we can get the same effect with the `&&` operator:

```julia
if 0 < x && x < 10
	println("x is a positive single-digit number.")
end
```

For this kind of condition, Julia provides a more concise syntax:

```julia
if 0 < x < 10
	println("x is a positive single-digit number.")
end
```
"""

# ╔═╡ 81410cd0-6bdd-11eb-3ac7-21c132df1a56
md"""## Recursion

It is legal for one function to call another; it is also legal for a function to call itself. It may not be obvious why that is a good thing, but it turns out to be one of the most magical things a program can do. For example, look at the following function:

```julia
function countdown(n)
	if n ≤ 0
		println("Blastoff!")
	else
		print(n, " ")
		countdown(n-1)
	end
end
```

If `n` is `0` or negative, it outputs the word `"Blastoff!"`. Otherwise, it outputs `n` and then calls a function named `countdown`—itself—passing `n-1` as an argument.

What happens if we call this function like this?

```julia
julia> countdown(3)
3 2 1 Blastoff!
```

* The execution of `countdown` begins with `n = 3`, and since `n` is greater than `0`, it outputs the value `3`, and then calls itself...
  * The execution of `countdown` begins with `n = 2`, and since `n` is greater than `0`, it outputs the value `2`, and then calls itself...

    * The execution of `countdown` begins with `n = 1`, and since `n` is greater than `0`, it outputs the value `1`, and then calls itself...

      * The execution of `countdown` begins with `n = 0`, and since `n` is not greater than `0`, it outputs the word `"Blastoff!"` and then returns.

    * The countdown that got `n = 1` returns.

  * The countdown that got `n = 2` returns.

* The countdown that got `n = 3` returns. And then you’re back in `Main`.

A function that calls itself is *recursive*; the process of executing it is called *recursion*.

As another example, we can write a function that prints a string ``n`` times:

```julia
function printn(s, n)
	if n ≤ 0
		return
	end
	println(s)
	printn(s, n-1)
end
```

If `n <= 0` the *`return` statement* exits the function. The flow of execution immediately returns to the caller, and the remaining lines of the function don’t run.

The rest of the function is similar to countdown: it displays s and then calls itself to display `s` `n-1` additional times. So, the number of lines of output is `1 + (n - 1)`, which adds up to `n`.

For simple examples like this, it is probably easier to use a `for` loop. But we will see examples later that are hard to write with a for loop and easy to write with recursion, so it is good to start early.
"""

# ╔═╡ 481f2708-6bde-11eb-2e2f-059722a2e288
md"""## Stack Diagrams for Recursive Functions

In “Stack Diagrams”, we used a stack diagram to represent the state of a program during a function call. The same kind of diagram can help us interpret a recursive function.

Every time a function gets called, Julia creates a frame to contain the function’s local variables and parameters. For a recursive function, there might be more than one frame on the stack at the same time.
"""

# ╔═╡ e7bc0384-6bdf-11eb-23d5-1f8454e713f2
Drawing(width=720, height=200) do
	text(x=320, y=30, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("Main")
	end
    rect(x=360, y=10, width=100, height=30, fill="rgb(242, 242, 242)", stroke="black")
	text(x=320, y=70, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("countdown")
	end
    rect(x=360, y=50, width=100, height=30, fill="rgb(242, 242, 242)", stroke="black")
	text(x=380, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("n")
	end
	text(x=400, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("->")
	end
	text(x=430, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("3")
	end
	text(x=320, y=110, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("countdown")
	end
    rect(x=360, y=90, width=100, height=30, fill="rgb(242, 242, 242)", stroke="black")
	text(x=380, y=110, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("n")
	end
	text(x=400, y=110, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("->")
	end
	text(x=430, y=110, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("2")
	end
	text(x=320, y=150, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("countdown")
	end
    rect(x=360, y=130, width=100, height=30, fill="rgb(242, 242, 242)", stroke="black")
	text(x=380, y=150, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("n")
	end
	text(x=400, y=150, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("->")
	end
	text(x=430, y=150, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	text(x=320, y=190, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("countdown")
	end
    rect(x=360, y=170, width=100, height=30, fill="rgb(242, 242, 242)", stroke="black")
	text(x=380, y=190, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("n")
	end
	text(x=400, y=190, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("->")
	end
	text(x=430, y=190, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("0")
	end
end

# ╔═╡ 85d66ec4-6be0-11eb-32c2-f7068bc06a62
md"*Figure 5-1. Stack diagram.*"

# ╔═╡ 0658bd3c-6be0-11eb-21bd-d753ae791777
md"""As usual, at the top of the stack is the frame for `Main`. It is empty because we did not create any variables in `Main` or pass any arguments to it.

The four `countdown` frames have different values for the parameter `n`. The bottom of the stack, where `n = 0`, is called the *base case*. It does not make a recursive call, so there are no more frames."""

# ╔═╡ c7a395b0-6be1-11eb-2d6c-8d15b718faa9
md"""#### Exercise 5-1

As an exercise, draw a stack diagram for printn called with `s = "Hello"` and `n = 2`. Then write a function called `do_n` that takes a function object and a number, `n`, as arguments, and that calls the given function `n` times.
"""

# ╔═╡ f94724b0-6be1-11eb-2183-1d97bab81764
md"""## Infinite Recursion

If a recursion never reaches a base case, it goes on making recursive calls forever, and the program never terminates. This is known as *infinite recursion*, and it is generally not a good idea. Here is a minimal program with an infinite recursion:

```julia
function recurse()
	recurse()
end
```

In most programming environments, a program with an infinite recursion does not really run forever. Julia reports an error message when the maximum recursion depth is reached:

```julia
julia> recurse()
ERROR: StackOverflowError:
Stacktrace:
     [1] recurse() at ./REPL[1]:2 (repeats 80000 times)
```

This stacktrace is a little bigger than the one we saw in the previous chapter. When the error occurs, there are 80,000 `recurse` frames on the stack!

If you encounter an infinite recursion by accident, review your function to confirm that there is a base case that does not make a recursive call. And if there is a base case, check whether you are guaranteed to reach it.
"""

# ╔═╡ b5bbb0ac-6be2-11eb-33fd-7971aa3a4c24
md"""## Keyboard Input

The programs we have written so far accept no input from the user. They just do the same thing every time.

Julia provides a built-in function called readline that stops the program and waits for the user to type something. When the user presses Return or Enter, the program resumes and readline returns what the user typed as a string:

```julia
julia> text = readline()
What are you waiting for?
"What are you waiting for?"
```

Before getting input from the user, it is a good idea to print a prompt telling her what to type:

```julia
julia> print("What...is your name? "); readline()
What...is your name? Arthur, King of the Britons!
"Arthur, King of the Britons!"
```

A semicolon (`;`) allows you to put multiple statements on the same line. In the REPL, only the last statement returns its value. semicolon (;) allows you to put multiple statements on the same line. In the REPL, only the last statement returns its value.

If you expect the user to type an integer, you can try to convert the return value to `Int64`:

```julia
julia> println("What...is the airspeed velocity of an unladen swallow?"); speed = readline()
What...is the airspeed velocity of an unladen swallow?
42
"42"
julia> parse(Int64, speed)
42
```

But if the user types something other than a string of digits, you get an error:

```julia
julia> println("What...is the airspeed velocity of an unladen swallow? "); speed = readline();
What...is the airspeed velocity of an unladen swallow?
What do you mean, an African or a European swallow?
julia> parse(Int64, speed)
ERROR: ArgumentError: invalid base 10 digit 'W' in "What do you mean, an African or a European swallow?"
[...]
```

We will see how to handle this kind of error later.

!!! info
    `readline` does not function in a Pluto notebook. However, a variable can be bound directly to a text input field:

    ```julia
    @bind x html"<input type=text>"
    ```

    Entering text in the text field will update the `x` variable and all top level statement using `x` will be reexecuted.

    All html input type are available.
"""

# ╔═╡ 96c3d1d8-6be3-11eb-1a7a-51786b2935ee
md"""## Debugging

When a syntax or runtime error occurs, the error message contains a lot of information. This can be overwhelming. The most useful parts are usually:

* What kind of error it was
* Where it occurred

Syntax errors are usually easy to find, but there are a few gotchas. In general, error messages indicate where the problem was discovered, but the actual error might be earlier in the code, sometimes on a previous line.

The same is true of runtime errors. Suppose you are trying to compute a signal-to- noise ratio in decibels. The formula is:

```math
\textit{SNR}_\mathrm{db} = 10 \log_{10} \frac{P_\mathrm{signal}}{P_\mathrm{noise}}
```

In Julia, you might write something like this:

```julia
signal_power = 9
noise_power = 10
ratio = signal_power ÷ noise_power
decibels = 10 * log10(ratio)
print(decibels)
```

And you’d get:

```julia
-Inf
```

This is not the result you expected.

To find the error, it might be useful to print the value of `ratio`, which turns out to be `0`. The problem is in line 3, which uses floor division instead of floatingpoint division.

!!! warning
    You should take the time to read error messages carefully, but don’t assume that everything they say is correct.
"""

# ╔═╡ acd667ba-6be3-11eb-00a0-55632ca86594
md"""## Glossary

*floor division*:
An operator, denoted `÷`, that divides two numbers and rounds down (toward negative infinity) to an integer.

*modulus operator*:
An operator, denoted with a percent sign (`%`), that works on integers and returns the remainder when one number is divided by another.

*Boolean expression*:
An expression whose value is either `true` or `false`.

*relational operator*:
One of the operators that compares its operands: `==`, `≠` (`!=`), `>`, `<`, `≥` (`>=`), and `≤` (`<=`).

*logical operator*:
One of the operators that combines Boolean expressions: `&&` (and), `||` (or), and `!` (not).

*conditional statement*:
A statement that controls the flow of execution depending on some condition.

*condition*:
The Boolean expression in a conditional statement that determines which branch runs.

*compound statement*:
A statement that consists of a header and a body. The body is terminated with the keyword `end`.

*branch*:
One of the alternative sequences of statements in a conditional statement.

*chained conditional*:
A conditional statement with a series of alternative branches.

*nested conditional*:
A conditional statement that appears in one of the branches of another conditional statement.

*recursive function*:
A function that calls itself.

*recursion*:
The process of calling the function that is currently executing.

*return statement*:
A statement that causes a function to end immediately and return to the caller.

*base case*:
A conditional branch in a recursive function that does not make a recursive call.

*infinite recursion*:
A recursion that doesn’t have a base case, or never reaches it. Eventually, an infinite recursion causes a runtime error.
"""

# ╔═╡ 81e8c81e-6be7-11eb-1052-434aee7931cf
md"""## Exercises

#### Exercise 5-2

The function `time` returns the current Greenwich Mean Time in seconds since “the epoch,” which is an arbitrary time used as a reference point. On Unix systems, the epoch is 1 January 1970:

```julia
julia> time()
1.612992834821836e9
```

Write a script that reads the current time and converts it to a time of day in hours, minutes, and seconds, plus the number of days since the epoch.
"""

# ╔═╡ bdd9c6c0-6be7-11eb-1cf0-65fbaf11abd0
md"""#### Exercise 5-3

Fermat’s Last Theorem says that there are no positive integers ``a``, ``b``, and ``c`` such that:

```math
a^n + b^n = c^n
```

for any value of ``n`` greater than ``2``.

1. Write a function named `checkfermat` that takes four parameters—`a`, `b`, `c`, and `n`—and checks to see if Fermat’s theorem holds. If `n` is greater than `2` and `a^n + b^n == c^n` the program should print, “Holy smokes, Fermat was wrong!” Otherwise, the program should print, “No, that doesn’t work.”

2. Write a function that prompts the user to input values for `a`, `b`, `c`, and `n`, converts them to integers, and uses `checkfermat` to check whether they violate Fermat’s theorem.
"""

# ╔═╡ 3a3ae2be-6be8-11eb-1fdc-db6aadd33a76
md"""#### Exercise 5-4

If you are given three sticks, you may or may not be able to arrange them in a triangle. For example, if one of the sticks is 12 inches long and the other two are 1 inch long, you will not be able to get the short sticks to meet in the middle. For any three lengths, there is a simple test to see if it is possible to form a triangle:

*If any of the three lengths is greater than the sum of the other two, then you cannot form a triangle. Otherwise, you can. (If the sum of two lengths equals the third, they form what is called a “degenerate” triangle.)*

1. Write a function named `istriangle` that takes three integers as arguments and prints either “Yes” or “No,” depending on whether you can or cannot form a triangle from sticks with the given lengths.

2. Write a function that prompts the user to input three stick lengths, converts them to integers, and uses `istriangle` to check whether sticks with the given lengths can form a triangle.
"""

# ╔═╡ 894fdcfe-6be8-11eb-021a-df6b152e0a2d
md"""#### Exercise 5-5

What is the output of the following program? Draw a stack diagram that shows the state of the program when it prints the result.

```julia
function recurse(n, s)
	if n == 0
		println(s)
	else
		recurse(n-1, n+s)
	end
end

recurse(3, 0)
```

1. What would happen if you called this function like this: `recurse(-1, 0)`?

2. Write a docstring that explains everything someone would need to know in order to use this function (and nothing else).
"""

# ╔═╡ 052ec49c-6bea-11eb-325f-d924ddacca6f
md"""
!!! info
    The following exercises use the `NativeSVG` module, described in Chapter 4.

#### Exercise 5-6

Read the following function and see if you can figure out what it does. Then run it and see if you got it right.

```julia
function draw(t, length, n)
	if n == 0
		return
	end
	angle = 50
	forward(t, length*n)
	turn(t, -angle)
	draw(t, length, n-1)
	turn(t, 2*angle)
	draw(t, length, n-1)
	turn(t, -angle)
	forward(t, -length*n)
end
```
"""

# ╔═╡ 5609252e-6bea-11eb-073c-0db3c1ab2da7
md"""#### Exercise 5-7

The Koch curve is a fractal that looks something like Figure 5-2.
"""

# ╔═╡ a70a425e-6beb-11eb-08e3-5d02826766e1
md"*Figure 5-2. A Koch curve.*"

# ╔═╡ b7655da0-6beb-11eb-21e4-a3372fe53c50
md"""To draw a Koch curve with length ``x``, all you have to do is:

1. Draw a Koch curve with length ``\frac{x}{3}``.
2. Turn left ``60`` degrees.
3. Draw a Koch curve with length ``\frac{x}{3}``.
4. Turn right ``120`` degrees.
5. Draw a Koch curve with length ``\frac{x}{3}``.
6. Turn left ``60`` degrees.
7. Draw a Koch curve with length ``\frac{x}{3}``.

The exception is if ``x`` is less than ``3``: in that case, you can just draw a straight line with length ``x``.

1. Write a function called `koch` that takes a turtle and a length as parameters, and that uses the turtle to draw a Koch curve with the given length.

2. Write a function called `snowflake` that draws three Koch curves to make the out‐ line of a snowflake.

3. The Koch curve can be generalized in several ways. See [https://en.wikipedia.org/ wiki/Koch_snowflake](https://en.wikipedia.org/ wiki/Koch_snowflake) for examples and implement your favorite.
"""

# ╔═╡ Cell order:
# ╟─867b290c-6bd9-11eb-1afb-47757ee13936
# ╟─4cc6b8ac-6bd9-11eb-1d23-f5d9caf02cd3
# ╟─81ea78c0-6bd9-11eb-1d3b-5f67fdf090f4
# ╟─fa7e737c-6bd9-11eb-3cd0-d5b1d2ae7af1
# ╟─01a13616-6bdb-11eb-1b8b-8bdda36e623f
# ╟─7478967a-6bdb-11eb-03f7-81f894751f6d
# ╟─e75638d2-6bdb-11eb-0a7e-93446ea41c2d
# ╟─482d8390-6bdc-11eb-3e3c-33b9f4f244b9
# ╟─9a3080ec-6be4-11eb-32eb-f162987e0fd6
# ╟─c0f31858-6bdc-11eb-00f8-9b4fb5bebb3e
# ╟─81410cd0-6bdd-11eb-3ac7-21c132df1a56
# ╟─481f2708-6bde-11eb-2e2f-059722a2e288
# ╟─e7bc0384-6bdf-11eb-23d5-1f8454e713f2
# ╟─85d66ec4-6be0-11eb-32c2-f7068bc06a62
# ╟─0658bd3c-6be0-11eb-21bd-d753ae791777
# ╟─c7a395b0-6be1-11eb-2d6c-8d15b718faa9
# ╟─f94724b0-6be1-11eb-2183-1d97bab81764
# ╟─b5bbb0ac-6be2-11eb-33fd-7971aa3a4c24
# ╟─96c3d1d8-6be3-11eb-1a7a-51786b2935ee
# ╟─acd667ba-6be3-11eb-00a0-55632ca86594
# ╟─81e8c81e-6be7-11eb-1052-434aee7931cf
# ╟─bdd9c6c0-6be7-11eb-1cf0-65fbaf11abd0
# ╟─3a3ae2be-6be8-11eb-1fdc-db6aadd33a76
# ╟─894fdcfe-6be8-11eb-021a-df6b152e0a2d
# ╟─052ec49c-6bea-11eb-325f-d924ddacca6f
# ╟─5609252e-6bea-11eb-073c-0db3c1ab2da7
# ╟─66ddc134-6bea-11eb-3bde-0fea30c863af
# ╟─a70a425e-6beb-11eb-08e3-5d02826766e1
# ╟─b7655da0-6beb-11eb-21e4-a3372fe53c50


text:
function sumBetweenElements(a::Array{Int, 1}, min::Integer, max::Integer, left = 1, right = length(a))
	if left == right
		if min ≤ a[left] ≤ max
			return a[left]
		else
			return 0
		end
	end
	med = (left+right) ÷ 2
	return sumBetweenElements(a, min, max, left, med) +
		   sumBetweenElements(a, min, max, med+1, right)
end


text:
function zad34(d, el1, el2, el3)
    A = zeros(d,d)
    for i in 1:d
        for j in 1:d
            if i + j == d+1
                A[i,j] = el1
            elseif i + j < d+1
                A[i,j] = el2
            else
                A[i,j] = el3
            end
        end
    end

    return A
end


text:
import UncertainData: UncertainDataset, UncertainIndexDataset, UncertainValueDataset, resample

surrogate(x::UncertainDataset, method::Surrogate) = surrogate(resample(x), method)
surrogate(x::UncertainValueDataset, method::Surrogate) = surrogate(resample(x), method)
surrogate(x::UncertainIndexDataset, method::Surrogate) = surrogate(resample(x), method)


text:
#-----------------------------------------------------------------------
#		EVALUATION SCRIPT
#           It is recommended that you run this script with 4 threads
#           enabled.
#-----------------------------------------------------------------------

using LinearAlgebra, Serialization, BenchmarkTools

#include("../src/squirrel.jl")
include("../src/metric.jl")

g  	= metric.g          # Gordon metric with standard parameters
gk 	= metric.ge         # Kerr-Schild with Earth parameters

Neval	= 1000          # Number of test cases to evaluate
Nsamp	= 1000          # Number of test cases in generated sample file

nb      = 24            # Number of steps for Broyden solver
tol	    = 1e-10         # Tolerance for ODE solver (OrdinaryDiffEq.jl)
tpfl    = Float64       # Floating point precision to use
ξ	    = 2e1           # Threshold for outlier detection code

sufx 	= ""                # Filename string suffix
dir 	= "../res/"         # Results directory
Nes	    = string(Neval)     # Turn Neval into a string
Nfs	    = string(Nsamp)     # Turn Nfile into a string

#-----------------------------------------------------------------------
#		DESERIALIZE GENERATED SAMPLE FILES (ANALOGUE GEOMETRY)
#-----------------------------------------------------------------------

pfx	= "tct" ;
tctloc	= dir*pfx*"-"*Nfs*sufx ;

tct	= Serialization.deserialize(tctloc) ;

tc	= squirrel.seval.tup2tc( tct ) ;

#-----------------------------------------------------------------------
#		DESERIALIZE GENERATED SAMPLE FILES (KERR GEOMETRY)
#-----------------------------------------------------------------------

pfx	= "tck" ;
tckloc	= dir*pfx*"-"*Nfs*sufx ;

tctk	= Serialization.deserialize(tckloc) ;

tck	 = squirrel.seval.tup2tc( tctk ) ;

#-----------------------------------------------------------------------
ne	= 5     # Number of emission points to consider
#-----------------------------------------------------------------------

# The following are parameters for the perturbation model
h0a  = tpfl[  0  ;  4  ;  8  ;  12  ;  16  ]
σa   = tpfl[  2  ; 1.5 ; 1.8 ;  1.7 ;  1.5 ]
h0i  = tpfl[  150  ;  200  ;  250  ;  300  ;  350  ]
σi   = tpfl[  21   ;  15   ;  18   ;  21   ;  10   ]

pfx  = "td"         # Prefix for filename

#-----------------------------------------------------------------------

# Fractional uncertainty of atmospheric n due to temperature and
# pressure uncertainties
δ1	= 0.001

# Fractional uncertainty of ionospheric n
δ2	= 0.10

# Atmospheric and ionospheric perturbation model
Patm 	= h->metric.P(h,h0a,σa)
Pion 	= h->metric.P(h,h0i,σi)

# Perturbed metric
gp 	    = x->metric.gp(x,δ1,δ2,Patm,Pion)

# Run evaluation function
tdL = squirrel.seval.main(tc,squirrel.locator,gp,Neval,tpfl,tol,ξ,nb,ne)

sfx	= "n"*string(ne)*"p"*string(Int(round(δ2*100)))    # Filename suffix

tdLtup = squirrel.seval.td2tup( tdL )

Serialization.serialize(dir*pfx*"-"*Nes*"-"*sfx*sufx,tdLtup)

#-----------------------------------------------------------------------

δ1	= 0.001
δ2	= 0.01
gp 	= x->metric.gp(x,δ1,δ2)
tdS = squirrel.seval.main(tc,squirrel.locator,gp,Neval,tpfl,tol,ξ,nb,ne)

sfx	= "n"*string(ne)*"p"*string(Int(round(δ2*100)))

tdStup = squirrel.seval.td2tup( tdS )

Serialization.serialize(dir*pfx*"-"*Nes*"-"*sfx*sufx,tdStup)

#-----------------------------------------------------------------------

td0 = squirrel.seval.main(tc,squirrel.locator,g,Neval,tpfl,tol,ξ,nb,ne)

sfx	= "n"*string(ne)*"p0"

td0tup = squirrel.seval.td2tup( td0 )

Serialization.serialize(dir*pfx*"-"*Nes*"-"*sfx*sufx,td0tup)

#-----------------------------------------------------------------------

tdk = squirrel.seval.main(tck,squirrel.locator,gk,Neval,tpfl,tol,ξ,nb,
                          ne)

sfx	= "n"*string(ne)*"k"

tdktup = squirrel.seval.td2tup( tdk )

Serialization.serialize(dir*pfx*"-"*Nes*"-"*sfx*sufx,tdktup)

#-----------------------------------------------------------------------


text:
using Random, LightGraphs, SimpleWeightedGraphs

@testset "Testing local_search.jl" begin

    @testset "Testing pearson_correlation" begin

        # Check equal
        b1 = [i/6 for i in 1:6]
        b2 = [i/6 for i in 1:6]

        @test AbrkgaObop.pearson_correlation(b1,b2,6) == 1.0

        # Check reverse
        b3 = reverse(b1)

        @test AbrkgaObop.pearson_correlation(b1,b3,6) ≈ -1.0
    end;

    @testset "Testing weighted_label_propagation" begin

        Random.seed!(0)

        keys = [0.16666666666666666 0.159 0.84;
                0.3333333333333333  0.32  0.82;
                0.5                 0.49  0.65;
                0.6666666666666666  0.65  0.49;
                0.8333333333333334  0.82  0.32;
                1.0                 0.84  0.159]

        # Create weighted graph
        sigma = 0.7                 # Used to make graph sparse

        n_itens = size(keys,1)
        n_voters = size(keys,2)

        graph = SimpleWeightedGraph(n_voters)
        for i in 1:n_voters-1
            for j in i+1:n_voters
                weight = AbrkgaObop.pearson_correlation(keys[:,i],keys[:,j],n_itens)
                if weight >= sigma
                    add_edge!(graph, i, j, weight)
                end
            end
        end

        # Return the index of the chromosomes that the LS will be applied
        index_local_search = AbrkgaObop.weighted_label_propagation(graph)

        @test index_local_search == [1,3]

    end

    @testset "Testing objective_partial" begin

        instance = read_dataset("datasets/movie.toi")

        bucket = [0 0;
                  0 0;
                  0 0;
                  0 0;
                  0 0;
                  0 0]
        keys = [0.16666666666666666 0.159;
                0.3333333333333333  0.32;
                0.5                 0.49;
                0.6666666666666666  0.65;
                0.8333333333333334  0.82;
                1.0                 0.84]
        fitness = [0,0]
        AbrkgaObop.decoder!(1,bucket,keys)
        AbrkgaObop.objective!(1,bucket,fitness,instance)

        # Check i -> j
        fitness_from_i_j = AbrkgaObop.objective_partial(bucket[:,1], fitness[1], 1, 1, 6, instance)
        @test fitness_from_i_j == 32

        # Check j -> i
        fitness_from_j_i = AbrkgaObop.objective_partial(bucket[:,1], fitness_from_i_j, 1, 6, 1, instance)
        @test fitness_from_j_i == fitness[1]

        # Compare with another bucket with same bucket order
        AbrkgaObop.decoder!(2,bucket,keys)
        AbrkgaObop.objective!(2,bucket,fitness,instance)
        @test fitness_from_i_j == AbrkgaObop.objective_partial(bucket[:,2], fitness[2], 1, 1, 6, instance)

    end

    @testset "Testing local_search" begin

        Random.seed!(0)

        instance = read_dataset("datasets/movie.toi")

        bucket = [0 0;
                  0 0;
                  0 0;
                  0 0;
                  0 0;
                  0 0]
        keys = [0.16666666666666666 0.159;
                0.3333333333333333  0.32;
                0.5                 0.49;
                0.6666666666666666  0.65;
                0.8333333333333334  0.82;
                1.0                 0.84]
        fitness = [0,0]

        statistics = AbrkgaObop.OBOPStatistics([],[],[],[],[],[],0.0,0.0,0,0)

        AbrkgaObop.decoder!(1,bucket,keys)
        AbrkgaObop.objective!(1,bucket,fitness,instance)

        best_solution = AbrkgaObop.OBOPSolution(instance.total_itens,time())

        AbrkgaObop.local_search!(1,keys,bucket,fitness,instance,best_solution,1,statistics)

        @test best_solution.objective == 72
        @test best_solution.bucket == [1,1,1,2,2,3]

    end

    @testset "Testing local_search_parallel" begin

        Random.seed!(0)

        instance = read_dataset("datasets/movie.toi")

        bucket = [0 0;
                  0 0;
                  0 0;
                  0 0;
                  0 0;
                  0 0]
        keys = [0.16666666666666666 0.159;
                0.3333333333333333  0.32;
                0.5                 0.49;
                0.6666666666666666  0.65;
                0.8333333333333334  0.82;
                1.0                 0.84]
        fitness = [0,0]

        statistics = AbrkgaObop.OBOPStatistics([],[],[],[],[],[],0.0,0.0,0,0)

        AbrkgaObop.decoder!(1,bucket,keys)
        AbrkgaObop.objective!(1,bucket,fitness,instance)

        best_solution = AbrkgaObop.OBOPSolution(instance.total_itens,time())

        AbrkgaObop.local_search_parallel!(1,keys,bucket,fitness,instance,best_solution,1,statistics)

        @test best_solution.objective == 72
        @test best_solution.bucket == [1,1,1,2,2,3]

    end

    #=
    @testset "Testing clustering_search" begin

        Random.seed!(0)

        instance = read_dataset("datasets/movie.toi")

        bucket = [0 0;
                  0 0;
                  0 0;
                  0 0;
                  0 0;
                  0 0]
        keys = [0.16666666666666666 0.159;
                0.3333333333333333  0.32;
                0.5                 0.49;
                0.6666666666666666  0.65;
                0.8333333333333334  0.82;
                1.0                 0.84]
        fitness = [0,0]
        AbrkgaObop.decoder!(1,bucket,keys)
        AbrkgaObop.objective!(1,bucket,fitness,instance)
        AbrkgaObop.decoder!(2,bucket,keys)
        AbrkgaObop.objective!(2,bucket,fitness,instance)

        best_solution = AbrkgaObop.OBOPSolution(instance.total_itens,time())

        statistics = AbrkgaObop.OBOPStatistics([],[],[],[],[],[],0.0,0.0,0,0)

        AbrkgaObop.clustering_search(2,
                                     [1,2],
                                     keys,
                                     bucket,
                                     fitness,
                                     instance,
                                     best_solution,
                                     1,
                                     statistics)

        @test best_solution.objective == 72
        @test best_solution.bucket == [1,1,1,2,2,3]

    end
    =#

end

max_stars_repo_name: LuizHNLorena/AbrkgaObop

max_issues_repo_name: LuizHNLorena/AbrkgaObop

max_forks_repo_name: LuizHNLorena/AbrkgaObop


text:
struct SubProblem{T <: AbstractFloat, A <: AbstractVector, S <: LQSolver}
    id::Int
    π::T
    solver::S
    c::A
    x::A
    y::A
    ρ::A
    optimvector::A

    function (::Type{SubProblem})(model::JuMP.Model,
                                  id::Integer,
                                  π::AbstractFloat,
                                  xdim::Integer,
                                  optimsolver::MPB.AbstractMathProgSolver)
        solver = LQSolver(model,optimsolver)
        solver()
        optimvector = getsolution(solver)
        x₀ = optimvector[1:xdim]
        y₀ = optimvector[xdim+1:end]
        T = promote_type(eltype(optimvector),Float32)
        c_ = convert(AbstractVector{T},JuMP.prepAffObjective(model))
        c_ *= model.objSense == :Min ? 1 : -1
        x₀_ = convert(AbstractVector{T},x₀)
        y₀_ = convert(AbstractVector{T},y₀)
        optimvector_ = convert(AbstractVector{T},optimvector)
        A = typeof(x₀_)
        subproblem = new{T,A,typeof(solver)}(id,
                                             π,
                                             solver,
                                             c_,
                                             x₀_,
                                             y₀_,
                                             zero(x₀_),
                                             optimvector_)
        return subproblem
    end

    function (::Type{SubProblem})(model::JuMP.Model,
                                  id::Integer,
                                  π::AbstractFloat,
                                  x₀::AbstractVector,
                                  y₀::AbstractVector,
                                  optimsolver::MPB.AbstractMathProgSolver)
        T = promote_type(eltype(x₀),eltype(y₀),Float32)
        c_ = convert(AbstractVector{T},JuMP.prepAffObjective(model))
        c_ *= model.objSense == :Min ? 1 : -1
        x₀_ = convert(AbstractVector{T},x₀)
        y₀_ = convert(AbstractVector{T},y₀)
        A = typeof(x₀_)
        solver = LQSolver(model,optimsolver)
        subproblem = new{T,A,typeof(solver)}(id,
                                             π,
                                             solver,
                                             c_,
                                             x₀_,
                                             y₀_,
                                             zero(x₀_),
                                             [x₀_...,y₀_...])
        return subproblem
    end
end

function update_subproblem!(subproblem::SubProblem, ξ::AbstractVector, r::AbstractFloat)
    subproblem.ρ[:] = subproblem.ρ + r*(subproblem.x - ξ)
end
update_subproblems!(subproblems::Vector{<:SubProblem}, ξ::AbstractVector, r::AbstractFloat) = map(prob -> update_subproblem!(prob,ξ,r), subproblems)

function reformulate_subproblem!(subproblem::SubProblem, ξ::AbstractVector, r::AbstractFloat)
    model = subproblem.solver.lqmodel
    # Linear part
    c = copy(subproblem.c)
    c[1:length(ξ)] += subproblem.ρ
    c[1:length(ξ)] -= r*ξ
    MPB.setobj!(model, c)
    # Quadratic part
    qidx = collect(1:length(subproblem.optimvector))
    qval = zeros(length(subproblem.optimvector))
    qval[1:length(ξ)] .= r
    if applicable(MPB.setquadobj!, model, qidx, qidx, qval)
        MPB.setquadobj!(model, qidx, qidx, qval)
    else
        error("Setting a quadratic penalty requires a solver that handles quadratic objectives")
    end
end
reformulate_subproblems!(subproblems::Vector{<:SubProblem}, ξ::AbstractVector, r::AbstractFloat) = map(prob -> reformulate_subproblem!(prob,ξ,r), subproblems)

function get_objective_value(subproblem::SubProblem)
    return subproblem.π*subproblem.c⋅subproblem.optimvector
end

function get_solution(subproblem::SubProblem)
    return copy(subproblem.y), getredcosts(subproblem.solver)[length(subproblem.x)+1:end], getduals(subproblem.solver)
end

function (subproblem::SubProblem)()
    subproblem.solver(subproblem.optimvector)
    solvestatus = status(subproblem.solver)
    if solvestatus == :Optimal
        xdim = length(subproblem.x)
        subproblem.optimvector[:] = getsolution(subproblem.solver)
        subproblem.x[:] = subproblem.optimvector[1:xdim]
        subproblem.y[:] = subproblem.optimvector[xdim+1:end]
        return get_objective_value(subproblem)
    elseif solvestatus == :Infeasible
        return Inf
    elseif solvestatus == :Unbounded
        return -Inf
    else
        error(@sprintf("Subproblem %d was not solved properly, returned status code: %s", subproblem.id, string(solvestatus)))
    end
end


text:
export
    ThreePassBP,
    # The export of three_pass_BP is kept for backward compatibility,
    # but new implementations should use ThreePassBP.
    three_pass_BP

"""
    ThreePassBP

    An instant algorithm that runs three passes of belief propagation.
"""
struct ThreePassBP <: BP
    default_range_size::Int
    ThreePassBP(drs = 10) = new(drs)
end

function debugparams(sf, range, prob, varname, fname, name)
    Dict(:type=>(output_type(sf) <: Real ? :cont : :discrete),
         :numBins=>min(10, length(range)),
         :range=>range,
         :prob=>prob,
         :varname=>varname,
         :fname=>fname,
         :name=>name)
end

function three_pass_BP(runtime::InstantRuntime)
    for node in get_nodes(get_network(runtime))
        remove_messages!(runtime, node, :pi_message)
        remove_messages!(runtime, node, :lambda_message)
    end
    run_bp(ThreePassBP(), runtime)
end

function run_bp(::ThreePassBP, runtime::InstantRuntime)
    network = get_network(runtime)
    variables = [v for v in topsort(get_initial_graph(network)) if v isa Variable]
    ranges = Dict{Symbol, Array{Any, 1}}()
    for node in get_nodes(network)
        inst = current_instance(runtime, node)
        ranges[node.name] = get_range(runtime, inst)
    end

    for var in variables
        for par in get_parents(network, var)
            if par isa Variable
                rng = ranges[par.name]
                initial_message = Cat(rng,  ones(length(rng)))
                set_message!(runtime, par, var, :pi_message, initial_message)
            end
        end
    end

    for var in variables
        _forwardstep(runtime, var, ranges, true)
    end

    for var in reverse(variables)
        _backstep(runtime, var, ranges)
    end

    for var in variables
        _forwardstep(runtime, var, ranges, false)
    end
end


text:
# This file defines the HistoryEntry type. This type marks when a change in the
#   status of a personnel member occurred.
# Note that the entry does NOT contain the parameter that is changed, because
#   that information is already contained in the representation of the personnel
#   member.

# The HistoryEntry type requires no additional types.
requiredTypes = [ ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

export HistoryEntry
type HistoryEntry
    # The timestamp of the entry.
    timestamp::Float64

    # The new value of the parameter.
    newVal

    function HistoryEntry( t::T, val ) where T <: Real
        newEntry = new()
        newEntry.timestamp = t
        newEntry.newVal = val
        return newEntry
    end  # HistoryEntry( t, par, newVal )
end  # type HistoryEntry


text:
__precompile__()

module CxxWrap

using Compat

export wrap_modules, wrap_module, wrap_module_types, wrap_module_functions, safe_cfunction, load_module, ptrunion, CppEnum, ConstPtr, ConstArray, gcprotect, gcunprotect

# Convert path if it contains lib prefix on windows
function lib_path(so_path::AbstractString)
  path_copy = so_path
  @static if is_windows()
    basedir, libname = splitdir(so_path)
    libdir_suffix = Sys.WORD_SIZE == 32 ? "32" : ""
    if startswith(libname, "lib") && !isfile(so_path)
      path_copy = joinpath(basedir*libdir_suffix, libname[4:end])
    end
  end
  return path_copy
end

const depsfile = joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl")
if !isfile(depsfile)
  error("$depsfile not found, CxxWrap did not build properly")
end
include(depsfile)
const jlcxx_path = _l_jlcxx

# Base type for wrapped C++ types
@compat abstract type CppAny end
@compat abstract type CppBits <: CppAny end
@compat abstract type CppDisplay <: Display end
@compat abstract type CppArray{T,N} <: AbstractArray{T,N} end
@compat abstract type CppAssociative{K,V} <: Associative{K,V} end

# Enum type interface
@compat abstract type CppEnum end
Base.convert(::Type{Int32}, x::CppEnum) = reinterpret(Int32, x)
import Base: +, |
+{T <: CppEnum}(a::T, b::T) = reinterpret(T, Int32(a) + Int32(b))
|{T <: CppEnum}(a::T, b::T) = reinterpret(T, Int32(a) | Int32(b))

cxxdowncast(x) = error("No downcast for type $(supertype(typeof(x))). Did you specialize SuperType to enable automatic downcasting?")

"""
Base class for smart pointers
"""
@compat abstract type SmartPointer{T} <: CppAny end

"""
Concrete smart pointer implementation. PT is a hash for the pointer types, DerefPtr a pointer to the dereference function,
ConstructPtr a function pointer to construct from a compatible smart pointer type (e.g. weak_ptr from shared_ptr)
CastPtr is a function pointer to cast to the direct base class
"""
type SmartPointerWithDeref{T,PT,DerefPtr,ConstructPtr,CastPtr} <: SmartPointer{T}
  ptr::Ptr{Void}
end

reference_type(t::Type) = Any

@generated function dereference_smart_pointer{T,PT,DerefPtr,ConstructPtr,CastPtr,DereferencedT}(p::SmartPointerWithDeref{T,PT,DerefPtr,ConstructPtr,CastPtr}, ::Type{DereferencedT})
  quote
    ccall(DerefPtr, $DereferencedT, (Ptr{Void},), p.ptr)
  end
end

function Base.getindex{T,PT,DerefPtr,ConstructPtr,CastPtr}(p::SmartPointerWithDeref{T,PT,DerefPtr,ConstructPtr,CastPtr})::reference_type(T)
  return dereference_smart_pointer(p, reference_type(T))
end

# No conversion if source and target type are identical
Base.convert{T,PT,DerefPtr,ConstructPtr,CastPtr}(::Type{SmartPointerWithDeref{T,PT,DerefPtr,ConstructPtr,CastPtr}}, p::SmartPointerWithDeref{T,PT,DerefPtr,ConstructPtr,CastPtr}) = p

# Construct from a related pointer, e.g. a std::weak_ptr from std::shared_ptr
function Base.convert{T,PT1,PT2,B1,B2,B3,D1,D2,D3}(::Type{SmartPointerWithDeref{T,PT1,B1,B2,B3}}, p::SmartPointerWithDeref{T,PT2,D1,D2,D3})
  ccall(B2, Any, (Any,), p)
end

# Construct from a related pointer, downcasting the type to base class
function Base.convert{BaseT,DerivedT,PT1,PT2,B1,B2,B3,D1,D2,D3}(t::Type{SmartPointerWithDeref{BaseT,PT1,B1,B2,B3}}, p::SmartPointerWithDeref{DerivedT,PT2,D1,D2,D3})
  if !(DerivedT <: BaseT)
    error("$DerivedT does not inherit from $BaseT in smart pointer convert")
  end
  # First convert to base type
  base_p = ccall(D3, Any, (Ptr{Void},), p.ptr)
  return convert(t, base_p)
end

# Cast to base type enclosed in same pointer type
function Base.convert{BaseT,DerivedT,PT,B1,B2,B3,D1,D2,D3}(::Type{SmartPointerWithDeref{BaseT,PT,B1,B2,B3}}, p::SmartPointerWithDeref{DerivedT,PT,D1,D2,D3})
  if !(DerivedT <: BaseT)
    error("$DerivedT does not inherit from $BaseT in smart pointer convert")
  end
  return convert(SmartPointerWithDeref{BaseT,PT,B1,B2,B3}, ccall(D3, Any, (Ptr{Void},), p.ptr))
end

immutable StrictlyTypedNumber{NumberT}
  value::NumberT
end
Base.convert{NumberT}(::Type{StrictlyTypedNumber{NumberT}}, n::NumberT) = StrictlyTypedNumber{NumberT}(n)

immutable ConstPtr{T} <: CppBits
  ptr::Ptr{T}
end

immutable ConstArray{T,N} <: CppArray{T,N}
  ptr::ConstPtr{T}
  size::NTuple{N,Int}
end

ConstArray{T,N}(ptr::ConstPtr{T}, args::Vararg{Int,N}) = ConstArray{T,N}(ptr, (args...))

@compat Base.IndexStyle(::ConstArray) = IndexLinear()
Base.size(arr::ConstArray) = arr.size
Base.getindex(arr::ConstArray, i::Integer) = unsafe_load(arr.ptr.ptr, i)

# Encapsulate information about a function
type CppFunctionInfo
  name::Any
  argument_types::Array{Type,1}
  reference_argument_types::Array{Type,1}
  return_type::Type
  function_pointer::Ptr{Void}
  thunk_pointer::Ptr{Void}
end

function __init__()
  @static if is_windows()
    Libdl.dlopen(jlcxx_path, Libdl.RTLD_GLOBAL)
  end
  ccall((:initialize, jlcxx_path), Void, (Any, Any, Any), CxxWrap, CppAny, CppFunctionInfo)
end

# Helper function to create a Module in a parent module
function create_module(name::String, parent::Module)
  return Core.eval(parent, :(module $(Symbol(name)) end))
end

# Load the modules in the shared library located at the given path
function load_modules(path::AbstractString, parent_module)
  module_lib = Libdl.dlopen(path, Libdl.RTLD_GLOBAL)
  registry = ccall((:create_registry, jlcxx_path), Ptr{Void}, (Any,Any), parent_module, nothing)
  ccall(Libdl.dlsym(module_lib, "register_julia_modules"), Void, (Ptr{Void},), registry)
  return registry
end

"""
Load a single module, to be wrapped in wrapped_module
"""
function load_module(path::AbstractString, wrapped_module)
  module_lib = Libdl.dlopen(path, Libdl.RTLD_GLOBAL)
  registry = ccall((:create_registry, jlcxx_path), Ptr{Void}, (Any,Any), module_parent(wrapped_module), wrapped_module)
  ccall(Libdl.dlsym(module_lib, "register_julia_modules"), Void, (Ptr{Void},), registry)
  return registry
end

function get_modules(registry::Ptr{Void})
  ccall((:get_modules, jlcxx_path), Any, (Ptr{Void},), registry)
end

function get_wrapped_module(registry::Ptr{Void})
  ccall((:get_wrapped_module, jlcxx_path), Any, (Ptr{Void},), registry)
end

function get_module_functions(registry::Ptr{Void})
  ccall((:get_module_functions, jlcxx_path), Any, (Ptr{Void},), registry)
end

function bind_constants(registry::Ptr{Void}, m::Module)
  ccall((:bind_module_constants, jlcxx_path), Void, (Ptr{Void},Any), registry, m)
end

function exported_symbols(registry::Ptr{Void}, modname::AbstractString)
  ccall((:get_exported_symbols, jlcxx_path), Any, (Ptr{Void},AbstractString), registry, modname)
end

function reference_types(registry::Ptr{Void}, modname::AbstractString)
  ccall((:get_reference_types, jlcxx_path), Any, (Ptr{Void},AbstractString), registry, modname)
end

function allocated_types(registry::Ptr{Void}, modname::AbstractString)
  ccall((:get_allocated_types, jlcxx_path), Any, (Ptr{Void},AbstractString), registry, modname)
end

"""
Protect a variable from garbage collection by adding it to the global array kept by CxxWrap
"""
function gcprotect(x)
  ccall((:gcprotect, jlcxx_path), Void, (Any,), x)
end

"""
Unprotect a variable from garbage collection by removing it from the global array kept by CxxWrap
"""
function gcunprotect(x)
  ccall((:gcunprotect, jlcxx_path), Void, (Any,), x)
end

# Interpreted as a constructor for Julia  > 0.5
type ConstructorFname
  _type::DataType
end

# Interpreted as an operator call overload
type CallOpOverload
  _type::DataType
end

process_fname(fn::Symbol) = fn
process_fname(fn::ConstructorFname) = :(::$(Type{fn._type}))
function process_fname(fn::CallOpOverload)
  return :(arg1::$(fn._type))
end

make_func_declaration(fn, argmap) = :($(process_fname(fn))($(argmap...)))
function make_func_declaration(fn::CallOpOverload, argmap)
  return :($(process_fname(fn))($((argmap[2:end])...)))
end

# By default, no argument overloading happens
argument_overloads(t::Type) = Type[]
@static if Int != Cint
  argument_overloads(t::Type{Cint}) = [Int]
end
@static if UInt != Cuint
  argument_overloads(t::Type{Cuint}) = [UInt, Int]
else
  argument_overloads(t::Type{Cuint}) = [Int]
end
argument_overloads(t::Type{Float64}) = [Int, Irrational]
function argument_overloads(t::Type{Array{AbstractString,1}})
  return [Array{String,1}]
end
argument_overloads{T <: Number}(t::Type{Ptr{T}}) = [Array{T,1}]

"""
Create a Union containing the type and a smart pointer to any type derived from it
"""
function ptrunion{T}(::Type{T})
  @compat result{T2 <: T} = Union{T2, SmartPointer{T2}}
  return result
end

smart_pointer_type(t::Type) = t
smart_pointer_type{T <: CppAny}(x::Type{T}) = ptrunion(x)
smart_pointer_type{T <: CppArray}(x::Type{T}) = ptrunion(x)
smart_pointer_type{T <: CppAssociative}(x::Type{T}) = ptrunion(x)

function smart_pointer_type{T,PT,DerefPtr,ConstructPtr,CastPtr}(::Type{SmartPointerWithDeref{T,PT,DerefPtr,ConstructPtr,CastPtr}})
  @compat result{T2 <: T} = SmartPointer{T2}
  return result
end

map_julia_arg_type(t::Type) = Union{smart_pointer_type(t),argument_overloads(t)...}
map_julia_arg_type{T}(a::Type{StrictlyTypedNumber{T}}) = T

# Build the expression to wrap the given function
function build_function_expression(func::CppFunctionInfo)
  # Arguments and types
  argtypes = func.argument_types
  argsymbols = map((i) -> Symbol(:arg,i[1]), enumerate(argtypes))

  # Function pointer
  fpointer = func.function_pointer
  assert(fpointer != C_NULL)

  # Thunk
  thunk = func.thunk_pointer

  map_c_arg_type(t::Type) = t
  map_c_arg_type{T <: AbstractString}(::Type{Array{T,1}}) = Any
  map_c_arg_type(::Type{Type}) = Any
  map_c_arg_type{T <: Tuple}(::Type{T}) = Any
  map_c_arg_type{T,N}(::Type{ConstArray{T,N}}) = Any
  map_c_arg_type{T <: SmartPointer}(::Type{T}) = Any

  map_return_type(t) = map_c_arg_type(t)
  map_return_type{T}(t::Type{Ref{T}}) = Ptr{T}

  # Build the types for the ccall argument list
  c_arg_types = [map_c_arg_type(t) for t in func.reference_argument_types]
  return_type = map_return_type(func.return_type)

  converted_args = ([:(Base.cconvert($t,$a)) for (t,a) in zip(func.reference_argument_types,argsymbols)]...)

  # Build the final call expression
  call_exp = nothing
  if thunk == C_NULL
    call_exp = :(ccall($fpointer, $return_type, ($(c_arg_types...),), $(converted_args...))) # Direct pointer call
  else
    call_exp = :(ccall($fpointer, $return_type, (Ptr{Void}, $(c_arg_types...)), $thunk, $(converted_args...))) # use thunk (= std::function)
  end
  assert(call_exp != nothing)

  # Build an array of arg1::Type1... expressions
  function argmap(signature)
    result = Expr[]
    for (t, s) in zip(signature, argsymbols)
      push!(result, :($s::$(map_julia_arg_type(t))))
    end
    return result
  end

  function_expression = :($(make_func_declaration(func.name, argmap(argtypes)))::$(func.return_type) = $call_exp)
  return function_expression
end

function wrap_reference_converters(registry, julia_mod)
  mod_name = string(module_name(julia_mod))
  reftypes = reference_types(registry, mod_name)
  alloctypes = allocated_types(registry, mod_name)
  for (rt, at) in zip(reftypes, alloctypes)
    st = supertype(at)
    Core.eval(Base, :(cconvert(::Type{$rt}, x::$rt) = x))
    Core.eval(Base, :(cconvert(::Type{$rt}, x::$at) = unsafe_load(reinterpret(Ptr{$rt}, pointer_from_objref(x)))))
    Core.eval(Base, :(cconvert{T <: $st}(t::Type{$rt}, x::T) = Base.cconvert(t, $(cxxdowncast)(x))))
    Core.eval(Base, :(cconvert(::Type{$rt}, x::$(SmartPointer{st})) = x[]))
    Core.eval(Base, :(cconvert{T <: $st}(t::Type{$rt}, x::$(SmartPointer){T}) = Base.cconvert(t, $(cxxdowncast)(x[]))))
    Core.eval(CxxWrap, :(reference_type(::Type{$st}) = $rt))
  end
end

# Wrap functions from the cpp module to the passed julia module
function wrap_functions(functions, julia_mod)
  basenames = Set([
    :getindex,
    :setindex!,
    :convert,
    :deepcopy_internal,
    :+,
    :*,
    :(==)
  ])

  cxxnames = Set([
    :cxxdowncast
  ])

  for func in functions
    if func.name ∈ basenames
      Core.eval(Base, build_function_expression(func))
    elseif func.name ∈ cxxnames
      Core.eval(CxxWrap, build_function_expression(func))
    else
      Core.eval(julia_mod, build_function_expression(func))
    end
  end
end

# Create modules defined in the given library, wrapping all their functions and types
function wrap_modules(registry::Ptr{Void})
  jl_modules = get_modules(registry)
  for jl_mod in jl_modules
    bind_constants(registry, jl_mod)
  end

  module_functions = get_module_functions(registry)
  for (jl_mod, mod_functions) in zip(jl_modules, module_functions)
    wrap_reference_converters(registry, jl_mod)
    wrap_functions(mod_functions, jl_mod)
  end

  for jl_mod in jl_modules
    exps = [Symbol(s) for s in exported_symbols(registry, string(jl_mod))]
    Core.eval(jl_mod, :(export $(exps...)))
  end
end

# Wrap modules in the given path
function wrap_modules(so_path::AbstractString, parent_mod=Main)
  registry = CxxWrap.load_modules(lib_path(so_path), parent_mod)
  wrap_modules(registry)
end

# Place the types for the module with the name corresponding to the current module name in the current module
function wrap_module_types(registry)
  jlmod = get_wrapped_module(registry)
  wanted_name = string(module_name(jlmod))
  bind_constants(registry, jlmod)

  exps = [Symbol(s) for s in exported_symbols(registry, wanted_name)]
  Core.eval(jlmod, :(export $(exps...)))
end

function wrap_module_functions(registry)
  modules = get_modules(registry)
  @assert length(modules) == 1
  mod_idx = 1

  if mod_idx == 0
    error("Module $wanted_name not found in C++")
  end

  jlmod = get_wrapped_module(registry)
  module_functions = get_module_functions(registry)
  wrap_reference_converters(registry, jlmod)
  wrap_functions(module_functions[mod_idx], jlmod)
end

# Place the functions and types into the current module
function wrap_module(registry::Ptr{Void})
  wrap_module_types(registry)
  wrap_module_functions(registry)
end

"""
Wrap C++ types into module m. Assumes the library adds only a single module with the same name as m.
"""
function wrap_module(so_path::AbstractString, m::Module)
  registry = CxxWrap.load_module(lib_path(so_path), m)
  wrap_module(registry)
end

immutable SafeCFunction
  fptr::Ptr{Void}
  return_type::Type
  argtypes::Array{Type,1}
end

safe_cfunction(f::Function, rt::Type, args::Tuple) = SafeCFunction(cfunction(f, rt, args), rt, [t for t in args])

wstring_to_julia(p::Ptr{Cwchar_t}, L::Int) = transcode(String, unsafe_wrap(Array, p, L))
wstring_to_cpp(s::String) = transcode(Cwchar_t, s)

end # module


text:
using SafeTestsets

@safetestset "Nystrom mesh tests" begin include("nystrommesh_test.jl") end

@safetestset "Pde tests" begin include("pde_test.jl") end

@safetestset "Kernels tests" begin include("kernels_test.jl") end

@safetestset "Density tests" begin include("density_test.jl") end

@safetestset "Potential tests" begin include("potential_test.jl") end

@safetestset "Operator tests" begin include("integraloperator_test.jl") end

@safetestset "Default assemble" begin include("assemble_test.jl") end

@safetestset "Density interpolation" begin include("dim_test.jl") end

@safetestset "Local density interpolation" begin include("ldim_test.jl") end

@safetestset "Martensen-Kussmaul tests" begin include("martensen_kussmaul_test.jl") end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
##############
# TM1 API #
##############

"""
Represents the API to interact with, either an actual TM1 instance,
or a mock API for testing purposes
"""
abstract type TM1API end

struct TM1WebAPI <: TM1API
  # could maybe change this to take host, port and ssl and construct the URI
  endpoint::HTTP.URI
end

const DEFAULT_API = TM1WebAPI(HTTP.URI("http://localhost:12345/api/v1/"))

using Base.Meta

"""
    @api_default function f(api, args...)
    ...
    end

For a method taking an `api` argument, add a new method without the `api` argument
that just calls the method with DEFAULT_API.
"""
macro api_default(func)
  call = func.args[1]
  has_kwargs = isexpr(call.args[2], :parameters)
  newcall = Expr(
    :call,
    call.args[1],
    (
      has_kwargs ? [Expr(:parameters, Expr(:..., :kwargs)); call.args[4:end]] :
      call.args[3:end]
    )...,
  )
  argnames = map(has_kwargs ? call.args[4:end] : call.args[3:end]) do expr
    isexpr(expr, :kw) && (expr = expr.args[1])
    isexpr(expr, Symbol("::")) && return expr.args[1]
    @assert isa(expr, Symbol)
    return expr
  end
  esc(Expr(
    :toplevel,
    :(Base.@__doc__ $func),
    Expr(
      :function,
      newcall,
      Expr(:block, :($(call.args[1])(DEFAULT_API, $(argnames...); kwargs...))),
    ),
  ))
end

####################
# Default API URIs #
####################

function api_uri(api::TM1WebAPI, path)
  merge(api.endpoint, path = api.endpoint.path * HTTP.URIs.escapepath(path))
end

function tm1_request(
  api::TM1API,
  request_method,
  endpoint;
  auth = AnonymousAuth(),
  handle_error = true,
  headers = Dict(),
  params = Dict(),
  allowredirects = true,
)
  authenticate_headers!(headers, auth)
  params = tm12json(params)
  api_endpoint = api_uri(api, endpoint)
  _headers = convert(Dict{String,String}, headers)
  !haskey(_headers, "User-Agent") && (_headers["User-Agent"] = "TM1-jl")
  if request_method == HTTP.get
    r = request_method(
      merge(api_endpoint, query = params),
      _headers,
      redirect = allowredirects,
      status_exception = false,
      idle_timeout = 20,
    )
  else
    r = request_method(
      string(api_endpoint),
      _headers,
      JSON.json(params),
      redirect = allowredirects,
      status_exception = false,
      idle_timeout = 20,
    )
  end
  handle_error && handle_response_error(r)
  return r
end

tm1_get(api::TM1API, endpoint = ""; options...) =
  tm1_request(api, HTTP.get, endpoint; options...)
tm1_post(api::TM1API, endpoint = ""; options...) =
  tm1_request(api, HTTP.post, endpoint; options...)
tm1_put(api::TM1API, endpoint = ""; options...) =
  tm1_request(api, HTTP.put, endpoint; options...)
tm1_delete(api::TM1API, endpoint = ""; options...) =
  tm1_request(api, HTTP.delete, endpoint; options...)
tm1_patch(api::TM1API, endpoint = ""; options...) =
  tm1_request(api, HTTP.patch, endpoint; options...)

tm1_get_json(api::TM1API, endpoint = ""; options...) =
  JSON.parse(HTTP.payload(tm1_get(api, endpoint; options...), String))
tm1_post_json(api::TM1API, endpoint = ""; options...) =
  JSON.parse(HTTP.payload(tm1_post(api, endpoint; options...), String))
tm1_put_json(api::TM1API, endpoint = ""; options...) =
  JSON.parse(HTTP.payload(tm1_put(api, endpoint; options...), String))
tm1_delete_json(api::TM1API, endpoint = ""; options...) =
  JSON.parse(HTTP.payload(tm1_delete(api, endpoint; options...), String))
tm1_patch_json(api::TM1API, endpoint = ""; options...) =
  JSON.parse(HTTP.payload(tm1_patch(api, endpoint; options...), String))

##############
# Pagination #
##############

has_page_links(r) = HTTP.hasheader(r, "Link")
get_page_links(r) = split(HTTP.header(r, "Link"), ",")

function find_page_link(links, rel)
  relstr = "rel=\"$(rel)\""
  for i = 1:length(links)
    if occursin(relstr, links[i])
      return i
    end
  end
  return 0
end

extract_page_url(link) = match(r"<.*?>", link).match[2:end-1]

function tm1_paged_get(
  api,
  endpoint;
  page_limit = Inf,
  start_page = "",
  handle_error = true,
  auth = AnonymousAuth(),
  headers = Dict(),
  params = Dict(),
  options...,
)
  authenticate_headers!(headers, auth)
  _headers = convert(Dict{String,String}, headers)
  !haskey(_headers, "User-Agent") && (_headers["User-Agent"] = "TM1-jl")
  if isempty(start_page)
    r = tm1_get(
      api,
      endpoint;
      handle_error = handle_error,
      headers = _headers,
      params = params,
      auth = auth,
      options...,
    )
  else
    @assert isempty(params) "`start_page` kwarg is incompatible with `params` kwarg"
    r = HTTP.get(start_page, headers = _headers)
  end
  results = HTTP.Response[r]
  page_data = Dict{String,String}()
  if has_page_links(r)
    page_count = 1
    while page_count < page_limit
      links = get_page_links(r)
      next_index = find_page_link(links, "next")
      next_index == 0 && break
      r = HTTP.get(extract_page_url(links[next_index]), headers = _headers)
      handle_error && handle_response_error(r)
      push!(results, r)
      page_count += 1
    end
    links = get_page_links(r)
    for page in ("next", "last", "first", "prev")
      page_index = find_page_link(links, page)
      if page_index != 0
        page_data[page] = extract_page_url(links[page_index])
      end
    end
  end
  return results, page_data
end

# for APIs which return just a list
function tm1_get_paged_json(api, endpoint = ""; options...)
  results, page_data = tm1_paged_get(api, endpoint; options...)
  return mapreduce(r -> JSON.parse(HTTP.payload(r, String)), vcat, results), page_data
end

# for APIs which return a Dict(key => list, "total_count" => count)
function tm1_get_paged_json(api, endpoint, key; options...)
  results, page_data = tm1_paged_get(api, endpoint; options...)
  local total_count
  list = mapreduce(vcat, results) do r
    dict = JSON.parse(HTTP.payload(r, String))
    total_count = dict["total_count"]
    dict[key]
  end
  list, page_data, total_count
end

##################
# Error Handling #
##################

function handle_response_error(r::HTTP.Response)
  if r.status >= 400
    message, docs_url, errors = "", "", ""
    body = HTTP.payload(r, String)
    try
      data = JSON.parse(body)
      message = get(data, "message", "")
      docs_url = get(data, "documentation_url", "")
      errors = get(data, "errors", "")
    catch
    end
    error(
      "Error found in TM1 reponse:\n",
      "\tStatus Code: $(r.status)\n",
      (
        (isempty(message) && isempty(errors)) ? ("\tBody: $body",) :
        ("\tMessage: $message\n", "\tDocs URL: $docs_url\n", "\tErrors: $errors")
      )...,
    )
  end
end