text:
using Bigleaf
using Test, StableRNGs
using Pipe: @pipe
using DataFrames, Dates, TimeZones
using Tables: columntable
using StaticArrays
using Statistics, StatsBase
using Missings

isapproxm(args...; kwargs...) = isapprox(args...; kwargs...);
isapproxm(::Missing, ::Missing; kwargs...) = true

tha48 = DataFrame(
    datetime = ZonedDateTime("2014-06-01T00:00:00+00:00","yyyy-mm-ddTHH:MM:SSzzzz") .+
        (0:47) .* Minute(30),
    ustar = Union{Missing, Float64}[0.5400000214576721, 0.4900000095367432, 0.4799999892711639, 0.449999988079071, 0.5099999904632568, 0.4600000083446503, 0.5400000214576721, 0.5199999809265137, 0.4600000083446503, 0.4099999964237213, 0.5299999713897705, 0.5099999904632568, 0.5199999809265137, 0.4099999964237213, 0.5400000214576721, 0.4600000083446503, 0.5099999904632568, 0.5, 0.6299999952316284, 0.5799999833106995, 0.6800000071525574, 0.6499999761581421, 0.7400000095367432, 0.7300000190734863, 0.7699999809265137, 0.7400000095367432, 0.7799999713897705, 0.8700000047683716, 0.7799999713897705, 0.6700000166893005, 0.5600000023841858, 0.6600000262260437, 0.7099999785423279, 0.6100000143051147, 0.5899999737739563, 0.6100000143051147, 0.5299999713897705, 0.449999988079071, 0.3400000035762787, 0.3300000131130219, 0.3300000131130219, 0.3700000047683716, 0.3799999952316284, 0.3799999952316284, 0.3600000143051147, 0.3700000047683716, 0.3600000143051147, 0.3400000035762787],
    Tair = [11.88000011444092, 11.67000007629395, 11.1899995803833, 10.80000019073486, 10.67000007629395, 10.13000011444092, 9.779999732971191, 9.5, 9.09000015258789, 8.800000190734863, 8.6899995803833, 8.899999618530273, 9.430000305175781, 9.710000038146973, 10.55000019073486, 11.19999980926514, 11.72999954223633, 12.68000030517578, 13.1899995803833, 13.46000003814697, 14.1899995803833, 14.73999977111816, 14.65999984741211, 14.8100004196167, 15.02999973297119, 14.98999977111816, 14.77999973297119, 15.35000038146973, 15.5, 15.61999988555908, 15.47000026702881, 16.20000076293945, 16.13999938964844, 15.56999969482422, 15.14999961853027, 15.18000030517578, 14.81999969482422, 14.60000038146973, 14.10000038146973, 13.23999977111816, 12.65999984741211, 12.46000003814697, 12.09000015258789, 11.85999965667725, 11.72000026702881, 11.64000034332275, 11.52999973297119, 11.43000030517578],
    pressure = [97.63999938964844, 97.62999725341797, 97.61000061035156, 97.61000061035156, 97.61000061035156, 97.61000061035156, 97.61000061035156, 97.61000061035156, 97.62999725341797, 97.66000366210938, 97.66999816894531, 97.69000244140625, 97.69000244140625, 97.69999694824219, 97.70999908447266, 97.69999694824219, 97.7300033569336, 97.72000122070312, 97.70999908447266, 97.69999694824219, 97.69999694824219, 97.69999694824219, 97.69999694824219, 97.70999908447266, 97.70999908447266, 97.72000122070312, 97.70999908447266, 97.70999908447266, 97.69999694824219, 97.68000030517578, 97.68000030517578, 97.66999816894531, 97.68000030517578, 97.66999816894531, 97.66999816894531, 97.66000366210938, 97.6500015258789, 97.63999938964844, 97.63999938964844, 97.6500015258789, 97.66000366210938, 97.66999816894531, 97.69000244140625, 97.68000030517578, 97.68000030517578, 97.69000244140625, 97.69000244140625, 97.69000244140625],
    H = Union{Missing, Float64}[-68.18000030517578, -48.54000091552734, -59.09999847412109, -60.11000061035156, -59.40000152587891, -48.91999816894531, -53.52999877929688, -40.75, -29.28000068664551, -33.36999893188477, -9.050000190734863, 2.400000095367432, 38.29999923706055, 38.45000076293945, 85.08000183105469, 108.7799987792969, 176.6799926757812, 210.6900024414062, 230.0800018310547, 227.4199981689453, 321.6499938964844, 358.8999938964844, 329.9299926757812, 317.3099975585938, 375.1900024414062, 336.489990234375, 260.7099914550781, 360.6400146484375, 264.2200012207031, 185.3200073242188, 123.7699966430664, 236.0200042724609, 207.5200042724609, 128.3399963378906, 66.05999755859375, 79.6500015258789, 44.93999862670898, 16.22999954223633, -22.35000038146973, -29.10000038146973, -47.79999923706055, -63.61999893188477, -61.38000106811523, -65.11000061035156, -59.77000045776367, -60.52000045776367, -49.02999877929688, -53.45000076293945],
    LE = Union{Missing, Float64}[9.9399995803833, 5.269999980926514, 3.980000019073486, missing, missing, missing, missing, missing, 5.190000057220459, 10.11999988555908, 26.27000045776367, 27.3700008392334, 23.52000045776367, 6.639999866485596, 52.2599983215332, 54.02000045776367, 92.27999877929688, 100.9199981689453, 111.7200012207031, 223.1000061035156, 200.7400054931641, 185.0500030517578, 170.9799957275391, 183.4900054931641, 187.6900024414062, 223.0299987792969, 161.6600036621094, 160.9700012207031, 153.7599945068359, 77.20999908447266, 118.4300003051758, 31.86000061035156, 116.4100036621094, 84.98999786376953, 62.2400016784668, 68.12999725341797, 29.88999938964844, 41.65000152587891, 28.77000045776367, 16.53000068664551, 16.04000091552734, 10.67000007629395, 2.880000114440918, 5.369999885559082, 1.820000052452087, 1.580000042915344, 4.139999866485596, 4.880000114440918],
    Rn = [-86.48999786376953, -84.19999694824219, -81.80999755859375, -77.9000015258789, -79.08999633789062, -79.08999633789062, -80.13999938964844, -73.36000061035156, -66.63999938964844, -37.93999862670898, -25.25, 78.0, 113.2399978637695, 121.9499969482422, 230.6199951171875, 302.1700134277344, 394.5400085449219, 453.8800048828125, 518.2100219726562, 609.9000244140625, 693.4099731445312, 729.1400146484375, 767.3599853515625, 778.1699829101562, 778.5599975585938, 778.239990234375, 606.7899780273438, 724.239990234375, 524.8499755859375, 425.1499938964844, 399.75, 515.6400146484375, 371.9100036621094, 304.1300048828125, 161.6799926757812, 176.8000030517578, 108.5299987792969, 35.54999923706055, -25.29999923706055, -72.19999694824219, -91.4800033569336, -92.30000305175781, -92.12000274658203, -90.70999908447266, -90.12000274658203, -89.77999877929688, -87.5, -86.76000213623047],
    G = [-4.934999942779541, -5.085000038146973, -5.135000228881836, -5.210000038146973, -5.28000020980835, -5.34499979019165, -5.5, -5.804999828338623, -6.239999771118164, -6.550000190734863, -6.385000228881836, -6.039999961853027, -5.034999847412109, -4.070000171661377, -3.019999980926514, -1.475000023841858, 0.5550000071525574, 1.534999966621399, 2.630000114440918, 7.550000190734863, 22.06500053405762, 17.09499931335449, 18.93000030517578, 15.5649995803833, 16.90500068664551, 12.55000019073486, 30.14999961853027, 23.71500015258789, 14.81999969482422, 15.06999969482422, 9.925000190734863, 6.525000095367432, 4.864999771118164, 5.394999980926514, 4.920000076293945, 3.289999961853027, 2.345000028610229, 1.559999942779541, 0.6499999761581421, -0.4550000131130219, -1.669999957084656, -2.795000076293945, -3.674999952316284, -4.295000076293945, -4.775000095367432, -5.190000057220459, -5.364999771118164, -5.440000057220459],
    VPD = [0.5745999813079834, 0.5633999824523925, 0.513700008392334, 0.4560999870300293, 0.4184000015258789, 0.3608999967575073, 0.3033999919891358, 0.2822999954223633, 0.2641999959945679, 0.2621000051498413, 0.2773000001907349, 0.3075999975204468, 0.3339999914169312, 0.3270999908447266, 0.373799991607666, 0.4267000198364258, 0.4580999851226807, 0.5138000011444092, 0.5421999931335449, 0.6954999923706054, 0.8619000434875488, 1.010499954223633, 1.014999961853027, 1.075800037384033, 1.090100002288818, 1.094499969482422, 0.9982000350952148, 1.085700035095215, 1.094299983978271, 1.09379997253418, 1.096599960327148, 1.164799976348877, 1.048600006103516, 0.9133999824523926, 0.8229999542236328, 0.8409000396728515, 0.7889999866485595, 0.7486999988555908, 0.6928999900817872, 0.601800012588501, 0.5798999786376953, 0.5870999813079834, 0.5499000072479248, 0.5432000160217285, 0.5427000045776367, 0.5219999790191651, 0.5199999809265137, 0.511299991607666],
    wind = Union{Missing, Float64}[4.210000038146973, 4.460000038146973, 4.539999961853027, 4.079999923706055, 3.950000047683716, 4.019999980926514, 3.269999980926514, 3.480000019073486, 3.240000009536743, 3.420000076293945, 3.440000057220459, 2.970000028610229, 3.349999904632568, 2.660000085830688, 2.779999971389771, 2.329999923706055, 2.220000028610229, 2.160000085830688, 2.519999980926514, 2.980000019073486, 2.359999895095825, 2.420000076293945, 2.920000076293945, 3.359999895095825, 2.759999990463257, 3.279999971389771, 3.410000085830688, 3.480000019073486, 3.039999961853027, 3.950000047683716, 2.700000047683716, 2.829999923706055, 2.730000019073486, 3.640000104904175, 2.75, 3.299999952316284, 2.630000114440918, 2.049999952316284, 1.879999995231628, 2.130000114440918, 2.400000095367432, 2.619999885559082, 2.640000104904175, 2.589999914169312, 2.589999914169312, 2.819999933242798, 2.779999971389771, 2.660000085830688],
    #
    # from Su2001 of entire tha-dataset (see tmpf below)
    Gb_h = Union{Missing, Float64}[0.1822661462615249, 0.1750029796115819, 0.16439640538330097, 0.15115204867322993, 0.17536657313260987, 0.1625005163806663, 0.18840537219804904, 0.18879242029415302, 0.17645201287966256, 0.1538975543389697, 0.21088990351067222, 0.2159981919177606, 0.26124601608024695, 0.2773512635093451, 0.31590536172758477, 0.40377826077805495, 0.4767947347202397, 0.56854692623082, 0.42750700682547943, 0.4686532838500691, 0.4713224944023968, 0.5367616432054471, 0.4390590733192645, 0.4365982669032558, 0.45028995760447604, 0.44342264360628164, 0.3870369356275578, 0.4091200789143718, 0.38965757992457295, 0.3741648297018856, 0.3598397227673043, 0.41828517646506513, 0.37765531002288344, 0.346427027859903, 0.2926818717409249, 0.3033011874663078, 0.2735927127688745, 0.23189051804555683, 0.13765984163452955, 0.1230228546581407, 0.10194295317288511, 0.11062842016540558, 0.11745785563008848, 0.11461165845075606, 0.107937932508708, 0.11269011080949713, 0.11690251453131727, 0.1020485594839973],
    Ga_h = Union{Missing, Float64}[0.050190555728311555, 0.04116957600962841, 0.03877814728450016, 0.03736361775914192, 0.047872524670215494, 0.03975837173989957, 0.06052647712603172, 0.05504593321238934, 0.04766632363024636, 0.03725385088675965, 0.05886451561135729, 0.06231168613359922, 0.06166420574946217, 0.05146825750242876, 0.07874566459614012, 0.07414026084897143, 0.09405110045117922, 0.09616429374872118, 0.1150966610946552, 0.09097297096976756, 0.13839881444493318, 0.1317378026984531, 0.13140678312765053, 0.11633916710154416, 0.1454360006411824, 0.12128624346147493, 0.12212101975187602, 0.14200569128979465, 0.1322214598437316, 0.0871694750533624, 0.08780626502144728, 0.11251759597304214, 0.12401547627348973, 0.07893328170976276, 0.08836491703204545, 0.08219874492644172, 0.07681770620097074, 0.06927204867424476, 0.042503890501350394, 0.036116991360095725, 0.03139917116540786, 0.03548953025999413, 0.03731866634395528, 0.03750740234735876, 0.034189028583981027, 0.03392950934600214, 0.03332805495857941, 0.030478849592838437],
);

i_tmp = () -> begin
    Dl=0.01
    df = aerodynamic_conductance!(copy(df, copycols=false); Gb_model=Su2001(),
        Dl, LAI=thal.LAI, zh=thal.zh, zr=thal.zr);
    @show df.Gb_h[1:48];
    @show df.Ga_h[1:48];
end

thal = (
    LAI = 7.6,   # leaf area index
    zh  = 26.5,  # average vegetation height (m)
    zr  = 42,    # sensor height (m)
    Dl  = 0.01,  # leaf characteristic dimension (m)
)

@testset "Bigleaf" begin
    @testset "bigleaf_constants" begin
        include("bigleaf_constants.jl")
    end
    @testset "util" begin
        include("util.jl")
    end
    @testset "unit_conversions" begin
        include("unit_conversions.jl")
    end
    @testset "filter_data" begin
        include("filter_data.jl")
    end
    @testset "meteorological_variables" begin
        include("meteorological_variables.jl")
    end
    @testset "sun_position" begin
        include("sun_position.jl")
    end
    @testset "potential_radiation" begin
        include("potential_radiation.jl")
    end
    @testset "stability_correction" begin
        include("stability_correction.jl")
    end
    @testset "surface_roughness" begin
        include("surface_roughness.jl")
    end
    @testset "boundary_layer_conductance" begin
        include("boundary_layer_conductance.jl")
    end
    @testset "aerodynamic_conductance" begin
        include("aerodynamic_conductance.jl")
    end
    @testset "surface_conductance" begin
        include("surface_conductance.jl")
    end
    @testset "evapotranspiration" begin
        include("evapotranspiration.jl")
    end
end


text:
module Walk1DMCTSVisualize

using Plots, Walk1DMDP, POMDPToolbox, MCTS

export Walk1DMCTSVis

mutable struct Walk1DMCTSVis
    mdp::Walk1D
    s::Walk1DState
    actions::Vector{Vector{Float64}}
    next_states::Vector{Vector{Float64}}
    ns::Vector{Int}
    radii::Vector{Vector{Float64}}
end
Walk1DMCTSVis(mdp::Walk1D, s::Walk1DState) = Walk1DMCTSVis(mdp, s, Vector{Float64}[], Int[], Vector{Float64}[], Float64[])

function circle(x,y,r)
    c = Shape(Plots.partialcircle(0,2π,20,r))
    translate!(c, x, y)
    c
end

@recipe function f(vis::Walk1DMCTSVis, index::Int, h::SimHistory)
    m = vis.mdp
    xlim := (-1, (m.p.threshx+1))
    ylim := (-m.p.threshx-1, m.p.threshx+1)
	aspect_ratio := :equal
    @series begin
        label := "actions"
        alpha := 0.25
        x = state_hist(h)[end-1].t+1
        y = vis.next_states[index]
        [circle(x,y,r) for (y,r) in zip(y, vis.radii[index])]
    end
    @series begin
        seriestype := :line
        color := :black
        label := "thresh+"
        x = [-2.0, m.p.t_max]
        y = [m.p.threshx, m.p.threshx]
        x, y
    end
    @series begin
        seriestype := :line
        color := :black
        label := "thresh-"
        x = [-2.0, m.p.t_max]
        y = [-m.p.threshx, -m.p.threshx]
        x, y
    end
    @series begin
        label := "path"
        color := :blue
        x = [s.t for s in state_hist(h)[1:end-1]]
        y = [s.x for s in state_hist(h)[1:end-1]]
        x, y
    end
    @series begin
        seriestype := :scatter
        label := "current position"
        t = state_hist(h)[end-1].t
        x = state_hist(h)[end-1].x
        if abs(x) < abs(m.p.threshx)
            color := :blue
        else
            color := :red
        end
        [t], [x]
    end
end

function MCTS.notify_listener(vis::Walk1DMCTSVis,dsb::DSBPlanner,s,a,sp,r,snode,sanode,spnode)
    if s == vis.s
        tree = get(dsb.tree)
        sol = dsb.solver
        actions = [tree.a_labels[c] for c in tree.children[snode]]
        if isempty(vis.actions) || (actions != vis.actions[end])  #only push the deltas
            push!(vis.actions, actions)
            push!(vis.next_states, [tree.s_labels[tree.transitions[c][1][1]].x for c in tree.children[snode]])  #deterministic transition
            push!(vis.ns, tree.total_n[snode])
            push!(vis.radii, [sol.r0_action/tree.total_n[snode]^sol.lambda_action for action in actions])
        end
    end
end
function MCTS.notify_listener(vis::Walk1DMCTSVis,asb::ASBPlanner,s,a,sp,r,snode,sanode,spnode)
    if s == vis.s
        tree = get(asb.tree)
        sol = asb.solver
        actions = [tree.a_labels[c] for c in tree.children[snode]]
        if isempty(vis.actions) || (actions != vis.actions[end])  #only push the deltas
            push!(vis.actions, actions)
            push!(vis.next_states, [tree.s_labels[tree.transitions[c][1][1]].x for c in tree.children[snode]])  #deterministic transition
            push!(vis.ns, tree.total_n[snode])
            push!(vis.radii, [tree.a_radius[x] for x in tree.children[snode]])
        end
    end
end

end # module

ext: jl


text:
import SchumakerSpline.evaluate
import SchumakerSpline.evaluate_integral

"""
    right_integral(f::UnivariateFunction, left::Real)
    right_integral(f::UnivariateFunction, left::Union{Date,DateTime})

This calculates the integral of a function from a left limit and returns it as a UnivariateFunction.
So if you were to then evaluate this integral function at a point x then you would
get the integral between left and x.

If a `Date`, `DateTime` is input then it is first converted to a scalar with the `years_from_global_base` function.

### Inputs
* `f` - A `UnivariateFunction`.
* `left` - A left limit (scalar)
### Returns
* A `UnivariateFunction`.
"""
function right_integral(f::UnivariateFunction, left::Real)
    indef_int = indefinite_integral(f)
    left_constant = evaluate(indef_int, left)
    return indef_int - left_constant
end

function right_integral(f::UnivariateFunction, left::Union{Date,DateTime})
    left_float = years_from_global_base(left)
    return right_integral(f, left_float)
end

function right_integral(f::Piecewise_Function, left::R) where R<:Real
    whole_number_of_intervals = length(f.starts_)
    which_interval_contains_left = searchsortedlast(f.starts_, left)
    number_of_intervals = whole_number_of_intervals - which_interval_contains_left + 1
    starts_    =  Array{R,1}(undef, number_of_intervals)
    starts_[1] = left
    first_indefinite_integral = indefinite_integral(f.functions_[which_interval_contains_left])
    functions_    =  Array{UnivariateFunction}(undef, number_of_intervals)
    if typeof(first_indefinite_integral) == UnivariateFunctions.Undefined_Function
        return Undefined_Function()
    end
    functions_[1] = first_indefinite_integral - evaluate(first_indefinite_integral, left)
    if number_of_intervals == 1
        return Piecewise_Function(starts_, functions_)
    end
    displacement  = which_interval_contains_left - 1
    for i in 2:number_of_intervals
        starts_[i]    = f.starts_[i + displacement]
        indef_int = indefinite_integral(f.functions_[i + displacement])
        if typeof(indef_int) == UnivariateFunctions.Undefined_Function
            return Piecewise_Function(starts_[1:i], vcat(functions_[1:(i-1)], Undefined_Function()  ) )
        end
        value_at_previous_end = evaluate(functions_[i-1], starts_[i])
        functions_[i] = indef_int - evaluate(indef_int, starts_[i]) + value_at_previous_end
    end
    return Piecewise_Function(starts_, functions_)
end

"""
    left_integral(f::UnivariateFunction, right::Real)
    left_integral(f::UnivariateFunction, right::Union{Date,DateTime})

This calculates the integral of a function from a right limit and returns it as a UnivariateFunction.
So if you were to then evaluate this integral function at a point x then you would
get the integral between right and x.

If a `Date`, `DateTime` is input then it is first converted to a scalar with the `years_from_global_base` function.

### Inputs
* `f` - A `UnivariateFunction`.
* `right` - A right limit (scalar)
### Returns
* A `UnivariateFunction`.
"""
function left_integral(f::UnivariateFunction, right::Real)
    indef_int = indefinite_integral(f)
    right_constant = evaluate(indef_int, right)
    return right_constant - indef_int
end

function left_integral(f::UnivariateFunction, right::Union{Date,DateTime})
    right_float = years_from_global_base(right)
    return left_integral(f, right_float)
end

function left_integral(f::Piecewise_Function, right::R) where R<:Real
    whole_number_of_intervals = length(f.starts_)
    which_interval_contains_right = searchsortedlast(f.starts_, right)
    number_of_intervals = which_interval_contains_right
    starts_    =  Array{R,1}(undef, number_of_intervals)
    starts_[number_of_intervals] = f.starts_[which_interval_contains_right]
    last_indefinite_integral = indefinite_integral(f.functions_[which_interval_contains_right])
    functions_    =  Array{UnivariateFunction}(undef, number_of_intervals)
    if typeof(last_indefinite_integral) == UnivariateFunctions.Undefined_Function
        return Undefined_Function()
    end
    functions_[number_of_intervals] = evaluate(last_indefinite_integral, right) - last_indefinite_integral
    if number_of_intervals == 1
        return Piecewise_Function(starts_, functions_)
    end
    for i in reverse(1:(which_interval_contains_right-1))
        starts_[i]    = f.starts_[i]
        value_at_next_start = evaluate(functions_[i+1], starts_[i+1])
        indef_int = indefinite_integral(f.functions_[i])
        if typeof(indef_int) == UnivariateFunctions.Undefined_Function
            return Piecewise_Function(starts_[(i+1):which_interval_contains_right], functions_[(i+1):which_interval_contains_right]   )
        end
        functions_[i] = evaluate(indef_int, starts_[i+1] ) - indef_int + value_at_next_start
    end
    return Piecewise_Function(starts_, functions_)
end

"""
    evaluate_integral(f::UnivariateFunction,left::Real, right::Real)
    evaluate_integral(f::UnivariateFunction,left::Union{Date,DateTime}, right::Union{Date,DateTime})

This calculates the integral of a function from a left limit to a right limit and returns a scalar.

If a `Date`, `DateTime` is input then it is first converted to a scalar with the `years_from_global_base` function.

### Inputs
* `f` - A `UnivariateFunction`.
* `left` - A left limit (scalar)
* `right` - A right limit (scalar)
### Returns
* A `Real`.
"""
function evaluate_integral(f::UnivariateFunction,left::Real, right::Real)
    indef_int  = indefinite_integral(f)
    left_eval  = evaluate(indef_int, left)
    right_eval = evaluate(indef_int, right)
    return (right_eval - left_eval)
end
function evaluate_integral(f::UnivariateFunction,left::Union{Date,DateTime}, right::Union{Date,DateTime})
    left_as_float  = years_from_global_base(left)
    right_as_float = years_from_global_base(right)
    return evaluate_integral(f, left_as_float, right_as_float)
end

function evaluate_integral(f::Piecewise_Function,left::Real, right::Real)
    return evaluate(right_integral(f, left), right)
end


text:
module straight_duct_examples
using FinEtools
using Gaston

function straight_duct_Q8_example()
    t0  =  time()

    rho = 1.21*phun("kg/m^3");# mass density
    c  = 343.0*phun("m/s");# sound speed
    bulk =  c^2*rho;
    omega =  54.5901*phun("rev/s")
    vn0 =  -1.0*phun("m/s")
    Lx = 10.0*phun("m");# length of the box, millimeters
    Ly = 1.0*phun("m"); # length of the box, millimeters
    n = 20;#number of elements along the length

    println("""

    Straight duct with anechoic termination.
    Example from Boundary element acoustics: Fundamentals and computer codes, TW Wu, page 44.
    Both real and imaginary components of the pressure should have amplitude of
    rho*c = $(rho*c).

    Two-dimensional model. Quadratic qquadrilateral mesh.
    """)

    fens,fes  =  Q8block(Lx,Ly,n,2); # Mesh
    bfes  =  meshboundary(fes)
    L0 = selectelem(fens,bfes,facing = true, direction = [-1.0 0.0])
    L10 = selectelem(fens,bfes,facing = true, direction = [+1.0 0.0])
    nLx = selectnode(fens,box = [0.0 Lx  0.0 0.0], inflate = Lx/1.0e5)

    geom  =  NodalField(fens.xyz)
    P  =  NodalField(zeros(ComplexF64,size(fens.xyz,1),1))

    numberdofs!(P)

    material = MatAcoustFluid(bulk,rho)
    femm  =  FEMMAcoust(IntegData(fes, GaussRule(2, 3)), material)

    S  =  acousticstiffness(femm, geom, P);
    C  =  acousticmass(femm, geom, P);

    E10femm  =  FEMMAcoustSurf(IntegData(subset(bfes,L10),GaussRule(1, 3)), material)
    D  =  acousticABC(E10femm, geom, P);

    E0femm  =  FEMMBase(IntegData(subset(bfes,L0), GaussRule(1,  3)))
    fi  =  ForceIntensity(-1.0im*omega*rho*vn0);
    F  =  distribloads(E0femm, geom, P, fi, 2);

    p = (-omega^2*S +omega*1.0im*D + C)\F
    scattersysvec!(P, p[:])

    println("Pressure amplitude bounds: ")
    println("  real $(minimum(real(P.values)))/$(maximum(real(P.values)))")
    println("  imag $(minimum(imag(P.values)))/$(maximum(imag(P.values)))")

    println("Total time elapsed  =  ",time() - t0,"s")

    File  =   "straight_duct_Q8.vtk"
    scalars = real(P.values);
    vtkexportmesh(File, connasarray(fes), geom.values, FinEtools.MeshExportModule.Q8;  scalars = [("Pressure", scalars)])
    @async run(`"paraview.exe" $File`)

    # @pyimport matplotlib.pyplot as plt
    # plt.style[:use]("seaborn-whitegrid")
    # fig = plt.figure()
    # ax = plt.axes()
    # ix = sortperm(geom.values[nLx,1])
    # ax[:plot](geom.values[nLx,1][ix], real(P.values)[nLx][ix], marker=:o, color = :blue, label = "real")
    # ax[:plot](geom.values[nLx,1][ix], imag(P.values)[nLx][ix], marker=:d,  color = :red, label  =  "imag")
    # plt.legend()
    # plt.show()
    f = figure()
    set(axis="normal", plotstyle="linespoints", linewidth=2, pointsize = 2, color = "black", xlabel = "x", ylabel = "Pressure", grid="on", title = "")
    ix = sortperm(geom.values[nLx,1])
    plot(geom.values[nLx,1][ix], real(P.values)[nLx][ix], legend = "real", marker = "ecircle")
    plot!(geom.values[nLx,1][ix], imag(P.values)[nLx][ix], legend = "imag", color = "red", marker = "edmd")
    figure(f)

    true

end # straight_duct_Q8_example

function straight_duct_T3_example()
    t0  =  time()

    rho = 1.21*phun("kg/m^3");# mass density
    c  = 343.0*phun("m/s");# sound speed
    bulk =  c^2*rho;
    omega =  54.5901*phun("rev/s")
    vn0 =  -1.0*phun("m/s")
    Lx = 10.0*phun("m");# length of the box, millimeters
    Ly = 1.0*phun("m"); # length of the box, millimeters
    n = 20;#number of elements along the length

    println("""

    Straight duct with anechoic termination.
    Example from Boundary element acoustics: Fundamentals and computer codes, TW Wu, page 44.
    Both real and imaginary components of the pressure should have amplitude of
    rho*c = $(rho*c).

    Triangle mesh.
    """)

    fens,fes  =  T3block(Lx,Ly,n,2); # Mesh
    bfes  =  meshboundary(fes)
    L0 = selectelem(fens,bfes,facing = true, direction = [-1.0 0.0])
    L10 = selectelem(fens,bfes,facing = true, direction = [+1.0 0.0])
    nLx = selectnode(fens,box = [0.0 Lx  0.0 0.0], inflate = Lx/1.0e5)

    geom  =  NodalField(fens.xyz)
    P  =  NodalField(zeros(ComplexF64,size(fens.xyz,1),1))

    numberdofs!(P)

    material = MatAcoustFluid(bulk,rho)
    femm  =  FEMMAcoust(IntegData(fes, TriRule(1)), material)

    S  =  acousticstiffness(femm, geom, P);
    C  =  acousticmass(femm, geom, P);

    E10femm  =  FEMMAcoustSurf(IntegData(subset(bfes,L10),GaussRule(1, 2)), material)
    D  =  acousticABC(E10femm, geom, P);

    E0femm  =  FEMMBase(IntegData(subset(bfes,L0), GaussRule(1,  2)))
    fi  =  ForceIntensity(-1.0im*omega*rho*vn0);
    F  =  distribloads(E0femm, geom, P, fi, 2);

    p = (-omega^2*S +omega*1.0im*D + C)\F
    scattersysvec!(P, p[:])

    println("Pressure amplitude bounds: ")
    println("  real $(minimum(real(P.values)))/$(maximum(real(P.values)))")
    println("  imag $(minimum(imag(P.values)))/$(maximum(imag(P.values)))")

    println("Total time elapsed  =  ",time() - t0,"s")

    File  =   "straight_duct.vtk"
    scalars = real(P.values);
    vtkexportmesh(File, connasarray(fes), geom.values, FinEtools.MeshExportModule.T3; scalars = [("Pressure", scalars)])
    @async run(`"paraview.exe" $File`)

    # @pyimport matplotlib.pyplot as plt
    # plt.style[:use]("seaborn-whitegrid")
    # fig = plt.figure()
    # ax = plt.axes()
    # ix = sortperm(geom.values[nLx,1])
    # ax[:plot](geom.values[nLx,1][ix], real(P.values)[nLx][ix], marker=:o, color = :blue, label = "real")
    # ax[:plot](geom.values[nLx,1][ix], imag(P.values)[nLx][ix], marker=:d,  color = :red, label  =  "imag")
    # plt.legend()
    # plt.show()
    # fig = figure()
    # figure(fig)
    set(axis="normal", plotstyle="linespoints", linewidth=2, pointsize = 2, color = "black", xlabel = "x", ylabel = "Pressure", grid="on", title = "")
    f = figure()
    ix = sortperm(geom.values[nLx,1])
    plot(geom.values[nLx,1][ix], real(P.values)[nLx][ix], legend = "real", marker = "ecircle")
    plot!(geom.values[nLx,1][ix], imag(P.values)[nLx][ix], legend = "imag", color = "red", marker = "edmd")
    figure(f)

    true

end # straight_duct_T3_example

function allrun()
    println("#####################################################")
    println("# straight_duct_Q8_example ")
    straight_duct_Q8_example()
    println("#####################################################")
    println("# straight_duct_T3_example ")
    straight_duct_T3_example()
end # function allrun

end # module straight_duct_examples


text:
export LocalComm

#This class is to create a stand in for Tpetra's local maps

"""
    LocalComm(::Comm{GID, PID, LID})

Creates a comm object that creates an error when inter-process communication is attempted, but still allows access to the correct process ID information
"""
struct LocalComm{GID <: Integer, PID <: Integer, LID <: Integer} <: Comm{GID, PID, LID}
    original::Comm{GID, PID, LID}
end

function barrier(comm::LocalComm)
    throw(InvalidStateError("Cannot call barrier on a local comm"))
end

function broadcastAll(comm::LocalComm, v::AbstractArray, root::Integer)
    throw(InvalidStateError("Cannot call broadcastAll on a local comm"))
end

function gatherAll(comm::LocalComm, v::AbstractArray)
    throw(InvalidStateError("Cannot call gatherAll on a local comm"))
end

function sumAll(comm::LocalComm, v::AbstractArray)
    throw(InvalidStateError("Cannot call sumAll on a local comm"))
end

function maxAll(comm::LocalComm, v::AbstractArray)
    throw(InvalidStateError("Cannot call maxAll on a local comm"))
end

function minAll(comm::LocalComm, v::AbstractArray)
    throw(InvalidStateError("Cannot call minAll on a local comm"))
end

function scanSum(comm::LocalComm, v::AbstractArray)
    throw(InvalidStateError("Cannot call scanSum on a local comm"))
end

function myPid(comm::LocalComm)
    myPid(comm.original)
end

function numProc(comm::LocalComm)
    numProc(comm.original)
end

function createDistributor(comm::LocalComm)
    throw(InvalidStateError("Cannot call createDistributor on a local comm"))
end


text:
#----------------------------------------
#-- Real Manhattan network, from OpenStreetMap data
#----------------------------------------

immutable Coordinates
  x::Float64
  y::Float64
end

type Manhattan <: TaxiProblem
  network::Network
  roadTime::SparseMatrixCSC{Float64, Int}
  roadCost::SparseMatrixCSC{Float64, Int}
  custs::Array{Customer,1}
  taxis::Array{Taxi,1}
  nTime::Int
  waitingCost::Float64
  sp::ShortPaths

  #--------------
  #Specific attributes
  distances::SparseMatrixCSC{Float64, Int}
  positions::Vector{Coordinates}
  tStart::DateTime
  tEnd::DateTime
  function Manhattan(;sp=false, emptyType=false)
    if emptyType
      c = new()
      return c
    else
      c = new()
      data = load("../src/Cities/Manhattan/manhattan.jld")
      c.network   = data["network"]
      c.distances = data["distances"]
      c.roadTime  = data["timings"]
      c.roadCost  = c.roadTime/100 #temporary
      c.positions = [Coordinates(i,j) for (i,j) in data["positions"]]
      if sp
          c.sp = shortestPaths(c.network, c.roadTime, c.roadCost)
      end
      c.custs = Customer[]
      c.taxis = Taxi[]
      c.nTime = 0
      return c
    end
  end
end

#Output the graph vizualization to pdf file (see GraphViz library)
function drawNetwork(pb::Manhattan, name::AbstractString = "graph")
  stdin, proc = open(`neato -n2 -Tpdf -o ../outputs/$name.pdf`, "w")
  write(stdin, "digraph  citygraph {\n")
  for i in vertices(pb.network)
    write(stdin, "$i [\"pos\"=\"$(pb.positions[i].x),$(pb.positions[i].y)!\"]\n")
  end
  for i in vertices(pb.network), j in out_neighbors(pb.network,i)
    write(stdin, "$i -> $j\n")
  end
  write(stdin, "}\n")
  close(stdin)
end

max_stars_repo_name: sebmart/NYC-timing

max_issues_repo_name: sebmart/NYC-timing

max_forks_repo_name: sebmart/NYC-timing


text:
"""
    DAMMplot(Tₛ::Array{Float64, 1}, θ::Array{Float64, 1}, Rₛ::Array{Float64, 1}, r::Int64)

Plot scatter of data and fitted DAMM surface

# Example
```julia-repl
julia> Tₛ = collect(15.0:2.5:40.0)
julia> θ = collect(0.2:0.05:0.7)
julia> Rₛ = [1.0, 1.2, 1.5, 2.0, 2.7, 3.8, 4.9, 6.7, 4.1, 2.0, 0.4]
julia> r = 10
julia> fig = DAMMplot(Tₛ, θ, Rₛ, r)
```
"""
function DAMMplot(Tₛ::Array{Float64, 1}, θ::Array{Float64, 1}, Rₛ::Array{Float64, 1}, r::Int64)
  fig = Figure()
  ax3D = Axis3(fig[1, 1])
  out = DAMMmat(Tₛ, θ, Rₛ, r)
  ax3D.xlabel = to_latex("T_{soil} (°C)");
  ax3D.ylabel = to_latex("\\theta (m^3 m^{-3})");
  ax3D.zlabel = to_latex("R_{soil} (\\mumol m^{-2} s^{-1})");
  data3D = Vec3f.(Tₛ, θ, Rₛ)
  p3D = scatter!(ax3D, data3D, markersize = 2500, strokewidth = 3,
	color = Rₛ, colormap = Reverse(:Spectral))
  s3D = surface!(ax3D, out.x, out.y, out.DAMM_Matrix, colormap = Reverse(:Spectral),
	transparency = true, alpha = 0.1, shading = false)
  w3D = wireframe!(ax3D, out.x, out.y, out.DAMM_Matrix, overdraw = true,
	transparency = true, color = (:black, 0.1));
  Colorbar(fig[1, 2], colormap = Reverse(:Spectral),
	   label = to_latex("R_{soil} (\\mumol m^{-2} s^{-1})"));
  fig
  return fig
end


text:
using DynamicGrids, Test
using DynamicGrids: inbounds, isinbounds, _cyclic_index, SimData, _unwrap, ismasked

@testset "boundary boundary checks are working" begin
    @testset "inbounds with Remove() returns index and false for an boundaryed index" begin
        @test inbounds((1, 1), (4, 5), Remove()) == ((1,1),true)
        @test inbounds((2, 3), (4, 5), Remove()) == ((2,3),true)
        @test inbounds((4, 5), (4, 5), Remove()) == ((4,5),true)
        @test inbounds((-3, -100), (4, 5), Remove()) == ((-3,-100),false)
        @test inbounds((0, 0), (4, 5), Remove()) == ((0,0),false)
        @test inbounds((2, 3), (3, 2), Remove()) == ((2,3),false)
        @test inbounds((2, 3), (1, 4), Remove()) == ((2,3),false)
        @test inbounds((200, 300), (2, 3), Remove()) == ((200,300),false)
    end
    @testset "inbounds with Wrap() returns new index and true for an boundaryed index" begin
        @test inbounds((-2,3), (10, 10), Wrap()) == ((8,3),true)
        @test inbounds((2,0), (10, 10), Wrap()) == ((2,10),true)
        @test inbounds((22,0), (10, 10), Wrap()) == ((2,10),true)
        @test inbounds((-22,0), (10, 10), Wrap()) == ((8,10),true)
    end
    @testset "isinbounds" begin
        @test isinbounds((4, 5), (4, 5)) == true
        @test isinbounds((200, 300), (2, 3)) == false
        @test isinbounds((-22,0), (10, 10)) == false
    end
    @testset "boundscheck objects" begin
        output = ArrayOutput(zeros(Int, 10, 10); tspan=1:10)
        sd = SimData(output.extent, Ruleset())
        @test inbounds((5, 5), sd) == ((5, 5), true)
        @test inbounds((5, 5), first(sd)) == ((5, 5), true)
        @test inbounds((12, 5), sd) == ((12, 5), false)
        sd_wrap = SimData(output.extent, Ruleset(; boundary=Wrap()))
        @test inbounds((5, 5), sd_wrap) == ((5, 5), true)
        @test inbounds((12, 5), sd_wrap) == ((2, 5), true)
        @test inbounds((12, 5), first(sd_wrap)) == ((2, 5), true)
    end
end

@testset "isinferred" begin
    @testset "unstable conditional" begin
        rule = let threshold = 20
            Cell() do x
                x > 1 ? 2 : 0.0
            end
        end
        output = ArrayOutput(rand(Int, 10, 10); tspan=1:10)
        @test_throws ErrorException isinferred(output, rule)
    end

    @testset "return type" begin
        rule = Neighbors{:a,:a}(Moore{1}(zeros(Bool, 3, 3))) do hood, x
            round(Int, x + sum(hood))
        end
        output = ArrayOutput((a=rand(Int, 10, 10),); tspan=1:10)
        @test isinferred(output, rule)
        output = ArrayOutput((a=rand(Bool, 10, 10),); tspan=1:10)
        @test_throws ErrorException isinferred(output, rule)
    end

    @testset "let blocks" begin
        a = 0.7
        rule = SetCell() do data, index, x
            add!(first(data), round(Int, a + x), index...)
        end
        output = ArrayOutput(zeros(Int, 10, 10); tspan=1:10)
        @test_throws ErrorException isinferred(output, Ruleset(rule))
        a = 0.7
        rule = let a = a
            SetCell() do data, index, x
                add!(first(data), round(Int, a), index...)
            end
        end
        output = ArrayOutput(zeros(Int, 10, 10); tspan=1:10)
        @test isinferred(output, Ruleset(rule))
    end

end

@testset "ismasked" begin
    output = ArrayOutput(zeros(2, 2); mask=Bool[1 0; 0 1], tspan=1:10)
    sd = SimData(output.extent, Ruleset())
    @test ismasked(sd, 1, 2) == true
    @test ismasked(sd, 2, 2) == false
    output_nomask = ArrayOutput(zeros(2, 2); tspan=1:10)
    sd = SimData(output_nomask.extent, Ruleset())
    @test ismasked(sd, 1, 2) == false
end

@testset "unwrap" begin
    @test _unwrap(1) == 1
    @test _unwrap(Val(:a)) == :a
    @test _unwrap(Aux(:a)) == :a
    @test _unwrap(Grid(:a)) == :a
    @test _unwrap(Aux{:x}) == :x
    @test _unwrap(Grid{:x}) == :x
    @test _unwrap(Val{:x}) == :x
end


text:
using Test
using ForestBiometrics

@test typeof(sdi_chart(500,8)) == Plots.Plot{Plots.GRBackend}
@test typeof(sdi_chart(500,8;maxsdi=600)) == Plots.Plot{Plots.GRBackend}


text:
@testset "Group wrapped in ArrayManifold" begin
    G = SpecialOrthogonal(3)
    M = Rotations(3)
    AG = ArrayManifold(G)
    @test base_group(AG) === G
    @test Manifolds.is_decorator_group(AG) === Val(true)

    eg = Matrix{Float64}(I, 3, 3)
    ω = [[1.0, 2.0, 3.0], [3.0, 2.0, 1.0]]
    x, y = [exp(M, eg, hat(M, eg, ωi)) for ωi in ω]
    v = hat(M, eg, [-1.0, 2.0, 0.5])

    e = Identity(AG)
    @test e === Identity(G)
    x2, y2 = ArrayMPoint(x), ArrayMPoint(y)
    v2 = ArrayTVector(v)

    @test identity(AG, x2) isa ArrayMPoint
    @test isapprox(G, identity(AG, x2).value, identity(G, x))
    @test identity(AG, e) === e
    @test_throws DomainError identity(AG, Identity(TranslationGroup(3)))

    eg = similar(x2)
    identity!(AG, eg, x2)
    @test isapprox(G, eg.value, identity(G, x))
    eg = similar(x2)
    identity!(AG, eg, e)
    @test isapprox(G, eg.value, identity(G, x))

    @test inv(AG, x2) isa ArrayMPoint
    @test isapprox(G, inv(AG, x2).value, inv(G, x))
    @test inv(AG, e) === e
    @test_throws DomainError inv(AG, Identity(TranslationGroup(3)))

    xinv = similar(x2)
    inv!(AG, xinv, x2)
    @test isapprox(G, xinv.value, inv(G, x))
    eg = similar(x2)
    inv!(AG, eg, e)
    @test isapprox(G, eg.value, e)

    @test compose(AG, x2, y2) isa ArrayMPoint
    @test isapprox(G, compose(AG, x2, y2).value, compose(G, x, y))
    @test compose(AG, x2, e) === x2
    @test compose(AG, e, x2) === x2

    xy = similar(x2)
    compose!(AG, xy, x2, y2)
    @test isapprox(G, xy.value, compose(G, x, y))
    xy = similar(x2)
    compose!(AG, xy, e, y2)
    @test isapprox(G, xy.value, compose(G, e, y))
    xy = similar(x2)
    compose!(AG, xy, x2, e)
    @test isapprox(G, xy.value, compose(G, x, e))

    for conv in (LeftAction(), RightAction())
        @test translate(AG, x2, y2, conv) isa ArrayMPoint
        @test isapprox(G, translate(AG, x2, y2, conv).value, translate(G, x, y, conv))

        xy = similar(x2)
        translate!(AG, xy, x2, y2, conv)
        @test isapprox(G, xy.value, translate(G, x, y, conv))

        @test inverse_translate(AG, x2, y2, conv) isa ArrayMPoint
        @test isapprox(G, inverse_translate(AG, x2, y2, conv).value, inverse_translate(G, x, y, conv))

        xinvy = similar(x2)
        inverse_translate!(AG, xinvy, x2, y2, conv)
        @test isapprox(G, xinvy.value, inverse_translate(G, x, y, conv))
    end

    for conv in (LeftAction(), RightAction())
        @test translate_diff(AG, y2, x2, v2, conv; atol = 1e-10) isa ArrayTVector
        @test isapprox(G, translate_diff(AG, y2, x2, v2, conv; atol = 1e-10).value, translate_diff(G, y, x, v, conv))

        vout = similar(v2)
        translate_diff!(AG, vout, y2, x2, v2, conv; atol = 1e-10)
        @test isapprox(vout.value, translate_diff(G, y, x, v, conv))

        @test inverse_translate_diff(AG, y2, x2, v2, conv; atol = 1e-10) isa ArrayTVector
        @test isapprox(G, inverse_translate_diff(AG, y2, x2, v2, conv; atol = 1e-10).value, inverse_translate_diff(G, y, x, v, conv))

        vout = similar(v2)
        inverse_translate_diff!(AG, vout, y2, x2, v2, conv; atol = 1e-10)
        @test isapprox(vout.value, inverse_translate_diff(G, y, x, v, conv))
    end
end


text:
const atom_names = Dict(
    "H"=>1,
    "He"=>2,
    "Li"=>3,
    "Be"=>4,
    "B"=>5,
    "C"=>6,
    "N"=>7,
    "O"=>8,
    "F"=>9,
    "Ar"=>18
)

const m_electron = 9.10938291E-31
const m_u = 1.660539040E-27
const m_au = m_u/m_electron
const proton_mass = 1836.0

const atom_masses = Dict(
    "H"=>1.008*m_au,
    "He"=>4.0026*m_au,
    "Li"=>6.94*m_au,
    "Be"=>9.0122*m_au,
    "B"=>10.81*m_au,
    "C"=>12.011*m_au,
    "N"=>14.007*m_au,
    "O"=>15.999*m_au,
    "F"=>18.998*m_au,
    "Ne"=>20.180*m_au,
    "Ar"=>39.948*m_au
)


text:
using Test
using UncertainData
using Distributions
using StaticArrays
using StatsBase
using KernelDensity

#############################################
# UncertainValues module
#############################################
include("uncertain_values/test_assign_distributions.jl")
include("uncertain_values/test_uncertain_values.jl")
include("uncertain_values/test_minmax.jl")
include("uncertain_values/test_certain_values.jl")
include("uncertain_values/test_merging.jl")

include("uncertain_values/populations/test_UncertainScalarPopulation.jl")
include("uncertain_values/populations/test_ConstrainedUncertainScalarPopulation.jl")

#############################################
# UncertainDatasets module
#############################################
include("uncertain_datasets/test_uncertain_datasets.jl")

#############################################
# Resampling uncertain values
#############################################
include("resampling/test_resampling_uncertain_values.jl")

#############################################
# Resampling uncertain tuples
#############################################
include("resampling/test_resampling_uncertain_tuples.jl")

#############################################
# Resampling uncertain datasets
#############################################

# Resampling schemes
include("resampling/resampling_schemes/test_BinnedResampling.jl")
include("resampling/resampling_schemes/test_BinnedWeightedResampling.jl")
include("resampling/resampling_schemes/test_BinnedMeanResampling.jl")
include("resampling/resampling_schemes/test_BinnedMeanWeightedResampling.jl")

include("resampling/resampling_schemes/test_ConstrainedIndexValueResampling.jl")
include("resampling/resampling_schemes/test_ConstrainedValueResampling.jl")
include("resampling/resampling_schemes/test_SequentialResampling.jl")
include("resampling/resampling_schemes/test_SequentialInterpolatedResampling.jl")
include("resampling/resampling_schemes/test_RandomSequences.jl")

# Define an array of uncertain values `uvals` that we can construct datasets from.
include("resampling/define_uncertain_values.jl")

# Resampling vectors of uncertain values
#---------------------------------------
include("resampling/uncertain_vectors/test_resampling_vectors.jl")
include("resampling/uncertain_vectors/test_resampling_vectors_apply_funcs.jl")
include("resampling/uncertain_vectors/test_resampling_vectors_constraints.jl")

# Resampling uncertain datasets
#-------------------------------
include("resampling/uncertain_datasets/test_resampling_datasets.jl")

# Resampling uncertain value datasets
#-------------------------------------
include("resampling/uncertain_datasets/test_resampling_datasets_uncertainvaluedataset.jl")
include("resampling/uncertain_datasets/test_resampling_datasets_uncertainvaluedataset_apply_funcs.jl")
include("resampling/uncertain_datasets/test_resampling_datasets_uncertainvaluedataset_constraints.jl")

# Resampling uncertain index datasets
#-------------------------------------
include("resampling/uncertain_datasets/test_resampling_datasets_uncertainindexdataset.jl")
include("resampling/uncertain_datasets/test_resampling_datasets_uncertainindexdataset_apply_funcs.jl")
include("resampling/uncertain_datasets/test_resampling_datasets_uncertainindexdataset_constraints.jl")

# Resampling uncertain index-value datasets
#-------------------------------------
include("resampling/uncertain_datasets/test_resampling_uncertainindexvaluedataset.jl")
include("resampling/uncertain_datasets/test_resampling_with_schemes.jl")

# Resampling uncertain vectors
#-----------------------------
include("resampling/uncertain_vectors/test_resampling_vectors.jl")

# Special resampling constraints
#-----------------------------
include("resampling/uncertain_datasets/sequential/test_resampling_sequential_increasing.jl")
include("resampling/uncertain_datasets/sequential/test_resampling_sequential_decreasing.jl")

# Resampling inplace.
#-----------------------------
include("resampling/test_resampling_inplace.jl")

#############################################
# Resampling uncertain datasets element-wise
#############################################
include("resampling/uncertain_datasets/test_resampling_abstractuncertainvaluedataset_elwise.jl")

############################################
# Interpolation and binning
#############################################
include("generic_interpolation/test_findall_nan_chunks.jl")
include("generic_interpolation/test_interpolate_nans.jl")
include("resampling/uncertain_datasets/test_interpolation.jl")
include("binning/test_binning.jl")

# Interpolation with resampling schemes
include("resampling/binning/test_bin_BinnedResampling.jl")
include("resampling/binning/test_bin_BinnedWeightedResampling.jl")

#############################
# Mathematics
#############################
include("mathematics/test_mathematics.jl")

############################################
# Uncertain statistics
#############################################
include("uncertain_statistics/uncertain_values/test_core_stats_values_point_estimates.jl")
include("uncertain_statistics/uncertain_values/test_core_stats_values_pair_estimates.jl")

include("uncertain_statistics/uncertain_datasets/test_core_stats_datasets_single_dataset_estimates.jl")
include("uncertain_statistics/uncertain_datasets/test_core_stats_datasets_pairwise_estimates.jl")

include("uncertain_statistics/test_hypothesistests.jl")
#include("uncertain_statistics/test_hypothesistests_timeseries.jl")

#########################################
# Resampling with constraints and models
#########################################

# #####################
# # Resampling schemes
# #####################
# include("resampling/resampling_schemes/test_ConstrainedResampling.jl")

# #######################
# # Sampling constraints
# #######################
include("sampling_constraints/test_sampling_constraints.jl")
include("sampling_constraints/test_constrain_certainvalue.jl")
include("sampling_constraints/test_constrain_population.jl")
include("sampling_constraints/test_constrain_uncertainvalues.jl")
include("sampling_constraints/test_constrain_uncertainvalues_kde.jl")
include("sampling_constraints/test_constrain_uncertaindatasets.jl")
include("sampling_constraints/test_constrain_uncertainvaluedatasets.jl")
include("sampling_constraints/test_constrain_uncertainindexdatasets.jl")

include("sampling_constraints/test_constrain_with_schemes.jl")


text:
# Simple general-purpose queue.

type EmptyQueue <: Exception
end

type UnderfullQueue <: Exception
end

type QueueNode{T}
    value::T
    next::(@compat Union{QueueNode, Void})
end

type Queue{T}
    front::(@compat Union{QueueNode{T}, Void})
    back::(@compat Union{QueueNode{T}, Void})
    length::UInt64

    function Queue()
        new(nothing, nothing, 0)
    end
end

isempty(q::Queue) = q.front === nothing
length(q::Queue) = q.length

function enqueue!{T}(q::Queue{T}, value)
    newback = QueueNode{T}(convert(T, value), nothing)
    if q.back === nothing
        q.front = q.back = newback
    else
        q.back.next = newback
        q.back = newback
    end
    q.length += 1
    q.back.value
end

# Enqueue into kth place.
function enqueue!{T}(q::Queue{T}, value, k)
    newnode = QueueNode{T}(convert(T, value), nothing)
    if k == 0
        newnode.next = q.front
        q.front = newnode
    elseif k > q.length
        throw(UnderfullQueue())
    else
        front = q.front
        for _ in 1:(k-1)
            front = front.next
        end
        newnode = front.next
        front.next = newnode
    end
    if newnode.next === nothing
        q.back = newnode
    end
    q.length += 1
    newnode.value
end

function dequeue!(q::Queue)
    if q.front === nothing
        throw(EmptyQueue())
    else
        value = q.front.value
        q.front = q.front.next
        if q.front === nothing
            q.back = nothing
        end
        q.length -= 1
        value
    end
end

function peek(q::Queue)
    if q.front === nothing
        throw(EmptyQueue())
    else
        return q.front.value
    end
end


text:
module BenchMissingArgmax

using BenchmarkTools
using Transducers

function manual_missing_argmax(xs)
    v = typemin(Base.nonmissingtype(eltype(xs)))
    n = 0
    @inbounds @simd for i in eachindex(xs)
        x = xs[i]
        if x !== missing
            if x > v
                v = x
                n = i
            end
        end
    end
    return n, v
end

random_missings(n = 10^3, th=2) = [abs(x) > th ? missing : x for x in randn(n)]

argext_step(should_update) =
    (old, (index, value)) ->
        if old === nothing || should_update(old[2], value)
            (index, value)
        else
            old
        end

xf_scanext(should_update) = Scan(argext_step(should_update), nothing)

xf_argmax = Enumerate() |>
    OfType(Tuple{Integer, Number}) |>
    xf_scanext(<)

suite = BenchmarkGroup()

let xs = randn(10)
    @assert manual_missing_argmax(xs)[1] == argmax(xs)
end

let xs = random_missings()
    @assert manual_missing_argmax(xs) == foldl(right, eduction(xf_argmax, xs))
end

suite["xf"] = @benchmarkable(
    foldl(right, ed; init=nothing),
    setup=(ed = eduction(xf_argmax, random_missings())))

rf = reducingfunction(
    xf_argmax,
    right;
    # simd = true,
)
suite["rf"] = @benchmarkable(
    transduce($rf, nothing, xs),
    setup=(xs = random_missings()))

suite["man"] = @benchmarkable(
    manual_missing_argmax(xs),
    setup=(xs = random_missings()))

end  # module
BenchMissingArgmax.suite


text:
######### StanVariational Bernoulli example  ###########

using StanVariational

bernoulli_model = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

bernoulli_data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

# Keep tmpdir across multiple runs to prevent re-compilation
tmpdir = joinpath(@__DIR__, "tmp")

sm = VariationalModel(
  "bernoulli", bernoulli_model; tmpdir = tmpdir)

rc = stan_variational(sm; data=bernoulli_data)

if success(rc)

  (samples, names) = read_variational(sm)

  # Show the same output in DataFrame format
  sdf = read_summary(sm)
  display(sdf)
  println()

  # Retrieve mean value of theta from the summary
  sdf[sdf.parameters .== :theta, :mean]

end


text:
module SpeedMapping
	using LinearAlgebra, ForwardDiff
	include("main.jl")
	export speedmapping
end


text:
using Pkg
using Test

@testset "Pkg.jl" begin
    # Write your tests here.
end


text:
export continuity

Space(d::IntervalOrSegment) = Chebyshev(d)
Space(d::FullSpace{<:Real}) = Chebyshev(Line())

Fun(::typeof(identity), d::IntervalOrSegment{T}) where {T<:Number} =
    Fun(Chebyshev(d), [mean(d), complexlength(d)/2])

## Calculus

# the default domain space is higher to avoid negative ultraspherical spaces
Integral(d::IntervalOrSegment,n::Integer) = Integral(Ultraspherical(1,d),n)

for Func in (:DefiniteIntegral,:DefiniteLineIntegral)
    @eval begin
        #TODO: this may be misleading
        $Func(d::IntervalOrSegment) = $Func(JacobiWeight(-0.5,-0.5,Chebyshev(d)))
        function $Func(α::Number,β::Number,d::IntervalOrSegment)
            @assert α == β
            @assert round(Int,α+.5) == α+.5
            @assert round(Int,α+.5) >= 0
            $Func(JacobiWeight(α,β,Ultraspherical(round(Int,α+.5),d)))
        end
        $Func(α::Number,β::Number) = $Func(α,β,ChebyshevInterval())
    end
end


text:
using POMDPs
using ContinuousPOMDPTreeSearchExperiments
using POMDPToolbox
using POMCP
using POMCPOW
using LightDarkPOMDPs
using ProgressMeter
using PmapProgressMeter
using ParticleFilters
using JLD
using CPUTime
using Powseeker

N = 500

@everywhere begin
    using Powseeker
    using POMDPs
    using ContinuousPOMDPTreeSearchExperiments
    using POMDPToolbox
    using POMCP
    using POMCPOW
    using ParticleFilters
    using CPUTime

    pomdp = PowseekerPOMDP()
    p_rng = MersenneTwister(4)
    solvers = Dict{String, Union{Solver,Policy}}(

        "bt_100_500" => begin
            rng3 = copy(p_rng)
            action_gen = GPSFirst(rng3)
            rollout_policy = Downhill(pomdp)
            tree_queries = 500
            node_updater = ObsAdaptiveSRFilter(pomdp, LowVarianceResampler(100), 0.05, rng3)
            solver = POMCPDPWSolver(next_action=action_gen,
                                    tree_queries=tree_queries,
                                    c = exp(35.0),
                                    max_depth=mdp(pomdp).duration+2,
                                    k_action=8.0,
                                    alpha_action=1/20,
                                    k_observation=8.0,
                                    alpha_observation=1/10,
                                    estimate_value=FORollout(rollout_policy),
                                    node_sr_belief_updater=node_updater,
                                    rng=rng3
                                   )
        end,

        "pomdpdpw_10k" => begin
            rng3 = copy(p_rng)
            action_gen = GPSFirst(rng3)
            rollout_policy = Downhill(pomdp)
            tree_queries = 10_000
            solver = POMCPDPWSolver(next_action=action_gen,
                                    tree_queries=tree_queries,
                                    c = exp(35.0),
                                    max_depth=mdp(pomdp).duration+2,
                                    k_action=8.0,
                                    alpha_action=1/20,
                                    k_observation=8.0,
                                    alpha_observation=1/10,
                                    estimate_value=FORollout(rollout_policy),
                                    # node_sr_belief_updater=node_updater,
                                    rng=rng3
                                   )
        end,

        "pomcpow_10k" => begin
            rng3 = copy(p_rng)
            action_gen = GPSFirst(rng3)
            rollout_policy = Downhill(pomdp)
            tree_queries = 10_000
            solver = POMCPOWSolver(next_action=action_gen,
                                    tree_queries=tree_queries,
                                    criterion=MaxUCB(exp(35)),
                                    final_criterion=MaxTries(),
                                    max_depth=mdp(pomdp).duration+2,
                                    k_action=8.0,
                                    alpha_action=1/20,
                                    k_observation=8.0,
                                    alpha_observation=1/10,
                                    estimate_value=FORollout(rollout_policy),
                                    rng=rng3
                                   )
        end,

        "pomcpow_10k_so" => begin
            rng3 = copy(p_rng)
            action_gen = GPSFirst(rng3)
            rollout_policy = SkiOver(pomdp)
            tree_queries = 10_000
            solver = POMCPOWSolver(next_action=action_gen,
                                    tree_queries=tree_queries,
                                    criterion=MaxUCB(exp(35)),
                                    final_criterion=MaxTries(),
                                    max_depth=mdp(pomdp).duration+2,
                                    k_action=8.0,
                                    alpha_action=1/20,
                                    k_observation=8.0,
                                    alpha_observation=1/10,
                                    estimate_value=FORollout(rollout_policy),
                                    rng=rng3
                                   )
        end,

        "pomcpow_50k" => begin
            rng3 = copy(p_rng)
            action_gen = GPSFirst(rng3)
            rollout_policy = Downhill(pomdp)
            tree_queries = 50_000
            solver = POMCPOWSolver(next_action=action_gen,
                                    tree_queries=tree_queries,
                                    criterion=MaxUCB(exp(35)),
                                    final_criterion=MaxTries(),
                                    max_depth=mdp(pomdp).duration+2,
                                    k_action=8.0,
                                    alpha_action=1/20,
                                    k_observation=8.0,
                                    alpha_observation=1/10,
                                    estimate_value=FORollout(rollout_policy),
                                    rng=rng3
                                   )
        end,

    )
end

solver_keys = collect(keys(solvers))
rewards = Dict{String, AbstractVector{Float64}}()
# state_hists = Dict{String, AbstractVector{AbstractVector{state_type(pomdp)}}}()
times = Dict{String, AbstractVector{Float64}}()
steps = Dict{String, AbstractVector{Int}}()

for (j, sk) in enumerate(solver_keys)
    s_rewards = SharedArray(Float64, N)
    # s_hists = SharedArray(AbstractVector{state_type(pomdp)}, N)
    s_times = SharedArray(Float64, N)
    s_steps = SharedArray(Float64, N)
    prog = Progress(N, desc="$sk ($j of $(length(solver_keys)))...")
    pmap(prog, 1:N) do i
        sim_rng = MersenneTwister(i)
        up_rng = MersenneTwister(i+100_000)
        policy = solve(solvers[sk], pomdp)
        up = ObsAdaptiveParticleFilter(deepcopy(pomdp), LowVarianceResampler(10_000), 0.05, up_rng)
        sim = HistoryRecorder(max_steps=40, rng=sim_rng)
        s_times[i] = @CPUelapsed hist = simulate(sim, deepcopy(pomdp), policy, up)
        s_rewards[i] = discounted_reward(hist)
        s_steps[i] = n_steps(hist)
    end
    rewards[sk] = sdata(s_rewards)
    @show mean(rewards[sk])
    # state_hists[sk] = sdata(s_hists)
    times[sk] = sdata(s_times)
    steps[sk] = sdata(s_steps)
end

for k in solver_keys
    println("$k mean: $(mean(rewards[k])) sem: $(std(rewards[k])/sqrt(N))")
    println("$k time: $(mean(times[k]))")
end

filename = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "data", "compare_$(Dates.format(now(), "E_d_u_HH_MM")).jld")
println("saving to $filename...")
@save(filename, solver_keys, rewards, times, steps)
println("done.")


text:
using ModernGL, GLFW
using MicroLogging

# GLFW callbacks
error_callback(error::Cint, description::Ptr{GLchar}) = error("GLFW ERROR: code $error msg: $description")
function window_size_callback(window::GLFW.Window, width::Cint, height::Cint)
    global glfwWidth = width
    global glfwHeight = height
    return nothing
end

function create_window_context(width, height)
    @static if is_apple()
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 4)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 1)
        GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
        GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, true)
    else
        GLFW.DefaultWindowHints()
    end
    # create window
    window = GLFW.CreateWindow(width, height, "Initializing...")
    window == C_NULL && error("could not open window with GLFW3.")
    # set callbacks
    GLFW.SetErrorCallback(error_callback)
    GLFW.SetWindowSizeCallback(window, window_size_callback)
    GLFW.MakeContextCurrent(window)
    GLFW.WindowHint(GLFW.SAMPLES, 4)
    return window
end

let previousTime = time()
    frameCount = 0
    global function updatefps!(window::GLFW.Window)
        currentTime = time()
        elapsedTime = currentTime - previousTime
        if elapsedTime > 0.25
            previousTime = currentTime
            fps = frameCount / elapsedTime
            s = @sprintf "RealSense @ fps: %.2f" fps
            GLFW.SetWindowTitle(window, s)
            frameCount = 0
        end
        frameCount = frameCount + 1
    end
end

function mipmap()
    glGenerateMipmap(GL_TEXTURE_2D)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_LINEAR)
end

# shader program
function shaderinfo(shaderID::GLuint)
    maxLength = Ref{GLsizei}(0)
    glGetShaderiv(shaderID, GL_INFO_LOG_LENGTH, maxLength)
    actualLength = Ref{GLsizei}(0)
    log = Vector{GLchar}(maxLength[])
    glGetShaderInfoLog(shaderID, maxLength[], actualLength, log)
    @info "Shader $shaderID info log: \n $(String(log))"
end

function compileshader(source::AbstractString, shaderType::GLenum)
    @info "Compiling GLSL shader from $source..."
    shader = readstring(source)
    shaderID = glCreateShader(shaderType)
    glShaderSource(shaderID, 1, Ptr{GLchar}[pointer(shader)], C_NULL)
    glCompileShader(shaderID)
    shaderTypeName = GLENUM(shaderType).name
    # get shader compile status
    compileResult = Ref{GLint}(-1)
    glGetShaderiv(shaderID, GL_COMPILE_STATUS, compileResult)
    shaderinfo(shaderID)
    compileResult[] == GL_FALSE && error("$shaderTypeName $shaderID did not compile!")
    @info "$shaderTypeName $shaderID compiled successfully."
    return shaderID
end

function programinfo(programID::GLuint)
    maxLength = Ref{GLsizei}(0)
    glGetProgramiv(programID, GL_INFO_LOG_LENGTH, maxLength)
    actualLength = Ref{GLsizei}(0)
    log = Vector{GLchar}(maxLength[])
    glGetProgramInfoLog(programID, maxLength[], actualLength, log)
    @info "Program $programID info log: \n $(String(log))"
end

function program_verbose_info(programID::GLuint)
    result = Ref{GLint}(-1)
    @debug "Program $programID verbose info:"
    glGetProgramiv(programID, GL_LINK_STATUS, result)
    @debug "GL_LINK_STATUS = $(result[])"
    glGetProgramiv(programID, GL_ATTACHED_SHADERS, result)
    @debug "GL_ATTACHED_SHADERS = $(result[])"
    glGetProgramiv(programID, GL_ACTIVE_ATTRIBUTES, result)
    @debug "GL_ACTIVE_ATTRIBUTES = $(result[])"
    for i in 0:result[]-1
        maxLength = Ref{GLsizei}(0)
        glGetProgramiv(programID, GL_ACTIVE_ATTRIBUTE_MAX_LENGTH, maxLength)
        actualLength = Ref{GLsizei}(0)
        attributeSize = Ref{GLint}(0)
        attributeType = Ref{GLenum}(0)
        name = Vector{GLchar}(maxLength[])
        glGetActiveAttrib(programID, i, maxLength[], actualLength, attributeSize, attributeType, name)
        attributeTypeName = GLENUM(attributeType[]).name
        if attributeSize[] > 1
            for j = 1:attributeSize[]
                longName = @sprintf "%s[%i]" name j
                location = glGetAttribLocation(programID, longName)
                @debug "  $attributeTypeName $(String(longName)) at location $location"
            end
        else
            location = glGetAttribLocation(programID, name)
            @debug "  $attributeTypeName $(String(name)) at location $location"
        end
    end
    glGetProgramiv(programID, GL_ACTIVE_UNIFORMS, result)
    @debug "GL_ACTIVE_UNIFORMS = $(result[])"
    for i in 0:result[]-1
        maxLength = Ref{GLsizei}(0)
        glGetProgramiv(programID, GL_ACTIVE_UNIFORM_MAX_LENGTH, maxLength)
        actualLength = Ref{GLsizei}(0)
        attributeSize = Ref{GLint}(0)
        attributeType = Ref{GLenum}(0)
        name = Vector{GLchar}(maxLength[])
        glGetActiveUniform(programID, i, maxLength[], actualLength, attributeSize, attributeType, name)
        attributeTypeName = GLENUM(attributeType[]).name
        if attributeSize[] > 1
            for j = 1:attributeSize[]
                longName = @sprintf "%s[%i]" name j
                location = glGetUniformLocation(programID, longName)
                @debug "  $attributeTypeName $(String(longName)) at location $location"
            end
        else
            location = glGetUniformLocation(programID, name)
            @debug "  $attributeTypeName $(String(name)) at location $location"
        end
    end
    programinfo(programID)
end

function validateprogram(programID::GLuint)
    validResult = Ref{GLint}(-1)
    glValidateProgram(programID)
    glGetProgramiv(programID, GL_VALIDATE_STATUS, validResult)
    @info "Program $programID: GL_VALIDATE_STATUS = $(validResult[])"
    program_verbose_info(programID)
    validResult[] == GL_FALSE && return false
    @info "Program $programID created successfully."
    return true
end

function createprogram(vertID::GLuint, fragID::GLuint)
    programID = glCreateProgram()
    @info "Creating program $programID..."
    @info "Attaching shaders: $vertID and $fragID..."
    glAttachShader(programID, vertID)
    glAttachShader(programID, fragID)
    glLinkProgram(programID)
    # checkout programe linking status
    linkingResult = Ref{GLint}(-1)
    glGetProgramiv(programID, GL_LINK_STATUS, linkingResult)
    linkingResult[] == GL_FALSE && @error "Program $programID did not created!"
    @assert validateprogram(programID)
    glDeleteShader(vertID)
    glDeleteShader(fragID)
    return programID
end

function createprogram(vertSource::AbstractString, fragSource::AbstractString)
    vertexShader = compileshader(vertSource, GL_VERTEX_SHADER)
    fragmentShader = compileshader(fragSource, GL_FRAGMENT_SHADER)
    return createprogram(vertexShader, fragmentShader)
end


text:
using CircuitModelZoo
using Test

@testset "API Tests" begin

    @test isfile(zoo_file("clts", "2.clt"))
    @test isfile(zoo_cnf_file("8.cnf"))
    @test isfile(zoo_dnf_file("8.dnf"))
    @test isfile(zoo_sdd_file("random.sdd"))
    @test isfile(zoo_vtree_file("random.vtree"))
    #@test isfile(zoo_jpc_file("test.jpc")) no jpc at the moment in the zoo
    @test isfile(zoo_psdd_file("little_4var.psdd"))

end


text:
# Adapted from:
# SOSDEMO2 --- Lyapunov Function Search
# Section 3.2 of SOSTOOLS User's Manual

@testset "SOSDEMO2 with $(factory.constructor)" for factory in sdp_factories
    @polyvar x[1:3]

    # Constructing the vector field dx/dt = f
    f = [-x[1]^3-x[1]*x[3]^2,
         -x[2]-x[1]^2*x[2],
    -x[3]+3*x[1]^2*x[3]-3*x[3]/(x[3]^2+1)]

    m = SOSModel(factory)

    # The Lyapunov function V(x):
    Z = vec(x).^2
    @variable m V Poly(Z)

    @constraint m V >= sum(x.^2)

    # dV/dx*(x[3]^2+1)*f <= 0
    P = dot(differentiate(V, x), f).num # the denominator is x[3]^2+1
    @constraint m P <= 0

    JuMP.optimize!(m)

    @test JuMP.primal_status(m) == MOI.FeasiblePoint

    @test iszero(removemonomials(JuMP.value(V), Z))
end


text:
# single position and single velocity
function __tracegeodesics(
    m::AbstractMetricParams{T},
    init_pos::AbstractVector{T},
    init_vel::AbstractVector{T},
    time_domain::Tuple{T,T},
    solver;
    callback,
    kwargs...
) where {T}
    prob = integrator_problem(m, init_pos, init_vel, time_domain)
    cbs = create_callback_set(m, callback)
    solve_geodesic(solver, prob; callback = cbs, kwargs...)
end

# columnar of positions and velocity
function __tracegeodesics(
    m::AbstractMetricParams{T},
    init_positions::AbstractArray{T,2},
    init_velocities::AbstractArray{T,2},
    time_domain::Tuple{T,T},
    solver;
    ensemble = EnsembleThreads(),
    callback,
    kwargs...
) where {T}
    prob = integrator_problem(
        m,
        @view(init_positions[:, 1]),
        @view(init_velocities[:, 1]),
        time_domain
    )
    ens_prob = EnsembleProblem(
        prob,
        prob_func = (prob, i, repeat) -> integrator_problem(
            m,
            @view(init_positions[:, i]),
            @view(init_velocities[:, i]),
            time_domain
        ),
        safetycopy = false
    )

    cbs = create_callback_set(m, callback)
    solve_geodesic(
        solver,
        ens_prob,
        ensemble;
        trajectories = size(init_velocities, 2),
        callback = cbs,
        kwargs...
    )
end

# indexables
function __tracegeodesics(
    m::AbstractMetricParams{T},
    init_positions,
    init_velocities,
    time_domain::Tuple{T,T},
    solver;
    ensemble = EnsembleThreads(),
    callback,
    kwargs...
) where {T}
    prob = integrator_problem(m, init_positions[1], init_velocities[1], time_domain)
    ens_prob = EnsembleProblem(
        prob,
        prob_func = (prob, i, repeat) ->
            integrator_problem(m, init_positions[i], init_velocities[i], time_domain),
        safetycopy = false
    )

    cbs = create_callback_set(m, callback)
    solve_geodesic(
        solver,
        ens_prob,
        ensemble;
        trajectories = length(init_velocities),
        callback = cbs,
        kwargs...
    )
end

function solve_geodesic(solver, prob, ensemble; solver_opts...)
    solve(prob, solver, ensemble; solver_opts...)
end

function solve_geodesic(solver, prob; solver_opts...)
    solve(prob, solver, ; solver_opts...)
end


text:
module Polyhedra

using SparseArrays
using LinearAlgebra

export Polyhedron

abstract type Library end
abstract type Polyhedron{T} end

using JuMP
export optimizer_with_attributes

coefficient_type(::Union{AbstractVector{T}, Type{<:AbstractVector{T}}}) where T = T
similar_type(::Type{<:Vector}, ::Int, ::Type{T}) where T = Vector{T}
similar_type(::Type{SparseVector{S, IT}}, ::Int, ::Type{T}) where {S, IT, T} = SparseVector{T, IT}

emptymatrix(::Type{MT}, m, n) where {MT<:AbstractMatrix} = MT(undef, m, n)
emptymatrix(::Type{SparseMatrixCSC{T, Int}}, m, n) where T = spzeros(T, m, n)
similar_type(::Type{<:Matrix}, ::Type{T}) where T = Matrix{T}
similar_type(::Type{SparseMatrixCSC{S, I}}, ::Type{T}) where {S, I, T} = SparseMatrixCSC{T, I}

# Interface/Definitions
include("dimension.jl")
include("elements.jl")
include("comp.jl")
include("representation.jl")
include("indices.jl")
include("incidence.jl")
include("triangulation.jl")
include("iterators.jl")
include("polyhedron.jl")

vvectortype(rep::Rep) = vvectortype(typeof(rep))
hvectortype(rep::Rep) = hvectortype(typeof(rep))
# TODO Only define vectortype for the type for Polyhedron v0.4
function vectortype(RepT::Type{<:Polyhedron})
    @assert hvectortype(RepT) == vvectortype(RepT)
    hvectortype(RepT)
end
vectortype(HRepT::Type{<:HRepresentation}) = hvectortype(HRepT)
vectortype(VRepT::Type{<:VRepresentation}) = vvectortype(VRepT)
vectortype(p::Rep) = vectortype(typeof(p))

vectortype(::Type{<:AbstractSparseArray{T}}) where T = SparseVector{T, Int}
vectortype(::Type{<:AbstractMatrix{T}}) where T = Vector{T}
vectortype(::Type{StaticArrays.SArray{Tuple{N, M}, T, 2, NM}}) where {N, M, T, NM} = StaticArrays.SVector{M, T}

hmatrixtype(RepT::Type{<:HRep}, T::Type) = matrixtype(similar_type(hvectortype(RepT), T))
vmatrixtype(RepT::Type{<:VRep}, T::Type) = matrixtype(similar_type(vvectortype(RepT), T))
matrixtype(::Type{<:AbstractVector{T}}) where T = Matrix{T}
matrixtype(::Type{<:AbstractSparseVector{T}}) where T = SparseMatrixCSC{T, Int}

function similar_type(::Type{ET}, ::Type{Tout}) where {Tout, ET<:Union{HRepElement, VRepElement, Rep}}
    similar_type(ET, FullDim(ET), Tout)
end
function similar_type(::Type{ET}, d::FullDim) where {ET<:Union{HRepElement, VRepElement, Rep}}
    similar_type(ET, d, coefficient_type(ET))
end

# Operations
include("repop.jl")
include("center.jl")
include("repelemop.jl")
include("aff.jl")
include("linearity.jl")
include("redundancy.jl")
include("projection.jl")
include("extended.jl")

# Implementations of representations
include("vecrep.jl")
include("mixedrep.jl")
include("lphrep.jl")
include("jump.jl")
include("matrep.jl")
include("liftedrep.jl")
include("doubledescription.jl") # FIXME move it after projection.jl once it stops depending on LiftedRep
include("interval.jl") # 1D polyhedron
include("planar.jl")
include("defaultlibrary.jl")

include("fulldim.jl")

# Optimization
include("opt.jl")
include("polyhedra_to_lp_bridge.jl")
include("vrep_optimizer.jl")
include("default.jl")
include("projection_opt.jl")

# Visualization
include("show.jl")
include("recipe.jl")
include("decompose.jl")

end # module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

function GitReference(repo::GitRepo, refname::AbstractString)
    ref_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_reference_lookup, :libgit2), Cint,
                  (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Cstring),
                   ref_ptr_ptr, repo.ptr, refname)
    return GitReference(repo, ref_ptr_ptr[])
end

function GitReference(repo::GitRepo, obj_oid::GitHash, refname::AbstractString = Consts.HEAD_FILE;
                      force::Bool=false, msg::AbstractString="")
    ref_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_reference_create, :libgit2), Cint,
                  (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Ptr{UInt8}, Ptr{GitHash}, Cint, Cstring),
                   ref_ptr_ptr, repo.ptr, refname, Ref(obj_oid), Cint(force),
                   isempty(msg) ? C_NULL : msg)
    return GitReference(repo, ref_ptr_ptr[])
end

"""
    LibGit2.isorphan(repo::GitRepo)

Check if the current branch is an "orphan" branch, i.e. has no commits. The first commit
to this branch will have no parents.
"""
function isorphan(repo::GitRepo)
    r = @check ccall((:git_repository_head_unborn, :libgit2), Cint,
                     (Ptr{Cvoid},), repo.ptr)
    r != 0
end

"""
    LibGit2.head(repo::GitRepo) -> GitReference

Return a `GitReference` to the current HEAD of `repo`.
"""
function head(repo::GitRepo)
    head_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_repository_head, :libgit2), Cint,
                  (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}), head_ptr_ptr, repo.ptr)
    return GitReference(repo, head_ptr_ptr[])
end

"""
    LibGit2.shortname(ref::GitReference)

Return a shortened version of the name of `ref` that's
"human-readable".

```julia-repl
julia> repo = LibGit2.GitRepo(path_to_repo);

julia> branch_ref = LibGit2.head(repo);

julia> LibGit2.name(branch_ref)
"refs/heads/master"

julia> LibGit2.shortname(branch_ref)
"master"
```
"""
function shortname(ref::GitReference)
    isempty(ref) && return ""
    Base.@gc_preserve ref begin
        name_ptr = ccall((:git_reference_shorthand, :libgit2), Cstring, (Ptr{Cvoid},), ref.ptr)
        name_ptr == C_NULL && return ""
        name = unsafe_string(name_ptr)
    end
    return name
end

"""
    LibGit2.reftype(ref::GitReference) -> Cint

Return a `Cint` corresponding to the type of `ref`:
  * `0` if the reference is invalid
  * `1` if the reference is an object id
  * `2` if the reference is symbolic
"""
function reftype(ref::GitReference)
    return ccall((:git_reference_type, :libgit2), Cint, (Ptr{Cvoid},), ref.ptr)
end

"""
    LibGit2.fullname(ref::GitReference)

Return the name of the reference pointed to by the
symbolic reference `ref`. If `ref` is not a symbolic
reference, return an empty string.
"""
function fullname(ref::GitReference)
    isempty(ref) && return ""
    reftype(ref) == Consts.REF_OID && return ""
    Base.@gc_preserve ref begin
        rname = ccall((:git_reference_symbolic_target, :libgit2), Cstring, (Ptr{Cvoid},), ref.ptr)
        rname == C_NULL && return ""
        name = unsafe_string(rname)
    end
    return name
end

"""
    LibGit2.name(ref::GitReference)

Return the full name of `ref`.
"""
function name(ref::GitReference)
    isempty(ref) && return ""
    Base.@gc_preserve ref begin
        name_ptr = ccall((:git_reference_name, :libgit2), Cstring, (Ptr{Cvoid},), ref.ptr)
        name_ptr == C_NULL && return ""
        name = unsafe_string(name_ptr)
    end
    return name
end

function branch(ref::GitReference)
    isempty(ref) && return ""
    str_ptr_ptr = Ref{Cstring}()
    Base.@gc_preserve ref begin
        @check ccall((:git_branch_name, :libgit2), Cint,
                      (Ptr{Cstring}, Ptr{Cvoid},), str_ptr_ptr, ref.ptr)
        str = unsafe_string(str_ptr_ptr[])
    end
    return str
end

function ishead(ref::GitReference)
    isempty(ref) && return false
    err = ccall((:git_branch_is_head, :libgit2), Cint,
                  (Ptr{Cvoid},), ref.ptr)
    return err == 1
end

function isbranch(ref::GitReference)
    isempty(ref) && return false
    err = ccall((:git_reference_is_branch, :libgit2), Cint,
                  (Ptr{Cvoid},), ref.ptr)
    return err == 1
end

function istag(ref::GitReference)
    isempty(ref) && return false
    err = ccall((:git_reference_is_tag, :libgit2), Cint,
                  (Ptr{Cvoid},), ref.ptr)
    return err == 1
end

function isremote(ref::GitReference)
    isempty(ref) && return false
    err = ccall((:git_reference_is_remote, :libgit2), Cint,
                  (Ptr{Cvoid},), ref.ptr)
    return err == 1
end

function Base.show(io::IO, ref::GitReference)
    println(io, "GitReference:")
    if isremote(ref)
        println(io, "Remote with name ", name(ref))
    elseif isbranch(ref)
        println(io, "Branch with name ", name(ref))
        if ishead(ref)
            println(io, "Branch is HEAD.")
        else
            println(io, "Branch is not HEAD.")
        end
    elseif istag(ref)
        println(io, "Tag with name ", name(ref))
    end
end

"""
    peel([T,] ref::GitReference)

Recursively peel `ref` until an object of type `T` is obtained. If no `T` is provided,
then `ref` will be peeled until an object other than a [`GitTag`](@ref) is obtained.

- A `GitTag` will be peeled to the object it references.
- A [`GitCommit`](@ref) will be peeled to a [`GitTree`](@ref).

!!! note
    Only annotated tags can be peeled to `GitTag` objects. Lightweight tags (the default)
    are references under `refs/tags/` which point directly to `GitCommit` objects.
"""
function peel(::Type{T}, ref::GitReference) where T<:GitObject
    obj_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_reference_peel, :libgit2), Cint,
                 (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Cint), obj_ptr_ptr, ref.ptr, Consts.OBJECT(T))
    return T(ref.owner, obj_ptr_ptr[])
end
peel(ref::GitReference) = peel(GitObject, ref)

"""
    LibGit2.ref_list(repo::GitRepo) -> Vector{String}

Get a list of all reference names in the `repo` repository.
"""
function ref_list(repo::GitRepo)
    sa_ref = Ref(StrArrayStruct())
    @check ccall((:git_reference_list, :libgit2), Cint,
                      (Ptr{StrArrayStruct}, Ptr{Cvoid}), sa_ref, repo.ptr)
    res = convert(Vector{String}, sa_ref[])
    free(sa_ref)
    res
end

"""
    LibGit2.create_branch(repo::GitRepo, bname::AbstractString, commit_obj::GitCommit; force::Bool=false)

Create a new branch in the repository `repo` with name `bname`, which
points to commit `commit_obj` (which has to be part of `repo`). If
`force` is `true`, overwrite an existing branch named `bname` if it
exists. If `force` is `false` and a branch already exists named `bname`,
this function will throw an error.
"""
function create_branch(repo::GitRepo,
                       bname::AbstractString,
                       commit_obj::GitCommit;
                       force::Bool=false)
    ref_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_branch_create, :libgit2), Cint,
                  (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Cstring, Ptr{Cvoid}, Cint),
                   ref_ptr_ptr, repo.ptr, bname, commit_obj.ptr, Cint(force))
    return GitReference(repo, ref_ptr_ptr[])
end

"""
    LibGit2.delete_branch(branch::GitReference)

Delete the branch pointed to by `branch`.
"""
function delete_branch(branch::GitReference)
    @check ccall((:git_branch_delete, :libgit2), Cint, (Ptr{Cvoid},), branch.ptr)
end

"""
    LibGit2.head!(repo::GitRepo, ref::GitReference) -> GitReference

Set the HEAD of `repo` to the object pointed to by `ref`.
"""
function head!(repo::GitRepo, ref::GitReference)
    ref_name = name(ref)
    @check ccall((:git_repository_set_head, :libgit2), Cint,
                  (Ptr{Cvoid}, Cstring), repo.ptr, ref_name)
    return ref
end

"""
    lookup_branch(repo::GitRepo, branch_name::AbstractString, remote::Bool=false) -> Union{GitReference, Nothing}

Determine if the branch specified by `branch_name` exists in the repository `repo`.
If `remote` is `true`, `repo` is assumed to be a remote git repository. Otherwise, it
is part of the local filesystem.

Return either a `GitReference` to the requested branch
if it exists, or [`nothing`](@ref) if not.
"""
function lookup_branch(repo::GitRepo,
                       branch_name::AbstractString,
                       remote::Bool=false)
    ref_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    branch_type = remote ? Consts.BRANCH_REMOTE : Consts.BRANCH_LOCAL
    err = ccall((:git_branch_lookup, :libgit2), Cint,
                 (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Ptr{UInt8}, Cint),
                  ref_ptr_ptr, repo.ptr, branch_name, branch_type)
    if err != Int(Error.GIT_OK)
        if err == Int(Error.ENOTFOUND)
            return nothing
        end
        if ref_ptr_ptr[] != C_NULL
            close(GitReference(repo, ref_ptr_ptr[]))
        end
        throw(Error.GitError(err))
    end
    return GitReference(repo, ref_ptr_ptr[])
end

"""
    upstream(ref::GitReference) -> Union{GitReference, Nothing}

Determine if the branch containing `ref` has a specified upstream branch.

Return either a `GitReference` to the upstream branch if it exists,
or [`nothing`](@ref) if the requested branch does not have an upstream counterpart.
"""
function upstream(ref::GitReference)
    isempty(ref) && return nothing
    ref_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    err = ccall((:git_branch_upstream, :libgit2), Cint,
                  (Ref{Ptr{Cvoid}}, Ptr{Cvoid},), ref_ptr_ptr, ref.ptr)
    if err != Int(Error.GIT_OK)
        if err == Int(Error.ENOTFOUND)
            return nothing
        end
        if ref_ptr_ptr[] != C_NULL
            close(GitReference(ref.owner, ref_ptr_ptr[]))
        end
        throw(Error.GitError(err))
    end
    return GitReference(ref.owner, ref_ptr_ptr[])
end

repository(ref::GitReference) = ref.owner

function target!(ref::GitReference, new_oid::GitHash; msg::AbstractString="")
    ref_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_reference_set_target, :libgit2), Cint,
             (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Ptr{GitHash}, Cstring),
             ref_ptr_ptr, ref.ptr, Ref(new_oid), isempty(msg) ? C_NULL : msg)
    return GitReference(ref.owner, ref_ptr_ptr[])
end

function GitBranchIter(repo::GitRepo, flags::Cint=Cint(Consts.BRANCH_LOCAL))
    bi_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_branch_iterator_new, :libgit2), Cint,
                  (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Cint), bi_ptr, repo.ptr, flags)
    return GitBranchIter(repo, bi_ptr[])
end

function Base.start(bi::GitBranchIter)
    ref_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    btype = Ref{Cint}()
    err = ccall((:git_branch_next, :libgit2), Cint,
                 (Ptr{Ptr{Cvoid}}, Ptr{Cint}, Ptr{Cvoid}),
                  ref_ptr_ptr, btype, bi.ptr)
    err != Int(Error.GIT_OK) && return (nothing, -1, true)
    return (GitReference(bi.owner, ref_ptr_ptr[]), btype[], false)
end

Base.done(bi::GitBranchIter, state) = Bool(state[3])

function Base.next(bi::GitBranchIter, state)
    ref_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    btype = Ref{Cint}()
    err = ccall((:git_branch_next, :libgit2), Cint,
                 (Ptr{Ptr{Cvoid}}, Ptr{Cint}, Ptr{Cvoid}),
                  ref_ptr_ptr, btype, bi.ptr)
    err != Int(Error.GIT_OK) && return (state[1:2], (nothing, -1, true))
    return (state[1:2], (GitReference(bi.owner, ref_ptr_ptr[]), btype[], false))
end

Base.iteratorsize(::Type{GitBranchIter}) = Base.SizeUnknown()

function Base.map(f::Function, bi::GitBranchIter)
    res = nothing
    s = start(bi)
    while !done(bi, s)
        val = f(s[1:2])
        if res === nothing
            res = Vector{typeof(val)}()
        end
        push!(res, val)
        val, s = next(bi, s)
    end
    return res
end

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
using DataFrames, UrlDownload, Dash, DashHtmlComponents, DashCoreComponents

df3 = DataFrame(urldownload("https://raw.githubusercontent.com/plotly/datasets/master/2011_us_ag_exports.csv"))

function generate_table(dataframe, max_rows = 10)
    html_table([
        html_thead(html_tr([html_th(col) for col in names(df3)])),
        html_tbody([
            html_tr([html_td(dataframe[r, c]) for c in names(dataframe)]) for r = 1:min(nrow(dataframe), max_rows)
        ]),
    ])
end

app = dash()

app.layout = html_div() do
    html_h4("US Agriculture Exports (2011)"),
    generate_table(df3, 10)
end

run_server(app, "0.0.0.0", debug=true)

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
if haskey(ENV, "BARON_EXEC")
    path = ENV["BARON_EXEC"]
else
    error("Unable to locate BARON executable. Make sure the solver has been separately downloaded, and that you properly set the BARON_EXEC environment variable.")
end

open(joinpath(@__DIR__, "path.jl"), "w") do io
    write(io, """const baron_exec = $(repr(path))\n""")
end


text:
#==============================================================================#
# SQS/test/runtests.jl
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

using AWSSQS
using Base.Test

AWSCore.set_debug_level(1)

#-------------------------------------------------------------------------------
# Load credentials...
#-------------------------------------------------------------------------------

aws = AWSCore.aws_config(region="ap-southeast-2")

#-------------------------------------------------------------------------------
# SQS tests
#-------------------------------------------------------------------------------

for q in sqs_list_queues(aws, "ocaws-jl-test-queue-")
    sqs_delete_queue(q)
end

test_queue = "ocaws-jl-test-queue-" * lowercase(Dates.format(now(Dates.UTC),
                                                            "yyyymmddTHHMMSSZ"))
@testset "test queue" begin
qa = sqs_create_queue(aws, test_queue)

qb = sqs_get_queue(aws, test_queue)

@test qa[:resource] == qb[:resource]

sqs_send_message(qa, "Hello!")

m = sqs_receive_message(qa)
@test m[:message] == "Hello!"

sqs_delete_message(qa, m)
sqs_flush(qa)

msgs = repeat(["test message"], outer = 10)
sqs_send_message_batch(qa, msgs) # this do not work
for i in 1:10
    m = sqs_receive_message(qa)
    @test m[:message] == "test message"
    sqs_delete_message(qa, m)
end
sqs_flush(qa)

info = sqs_get_queue_attributes(qa)
@test info["ApproximateNumberOfMessages"] == "0"
@test sqs_count(qa) == 0

sqs_delete_queue(qa)

end
#==============================================================================#
# End of file.
#==============================================================================#


text:
using Test
import GDAL
import ArchGDAL; const AG = ArchGDAL

@testset "Test methods for rasterband" begin
    AG.read("data/utmsmall.tif") do dataset
        ds_result = """
        GDAL Dataset (Driver: GTiff/GeoTIFF)
        File(s):
          data/utmsmall.tif

        Dataset (width x height): 100 x 100 (pixels)
        Number of raster bands: 1
          [GA_ReadOnly] Band 1 (Gray): 100 x 100 (UInt8)
        """
        @test sprint(print, dataset) == ds_result
        rb = AG.getband(dataset, 1)
        sprint(print, rb) == """
        [GA_ReadOnly] Band 1 (Gray): 100 x 100 (UInt8)
        blocksize: 100×81, nodata: -1.0e10, units: 1.0px + 0.0
        overviews: """
        @test sprint(print, AG.getdataset(rb)) == ds_result

        @test AG.getunittype(rb) == ""
        AG.setunittype!(rb,"ft")
        @test AG.getunittype(rb) == "ft"
        AG.setunittype!(rb,"")
        @test AG.getunittype(rb) == ""

        @test AG.getoffset(rb) == 0
        AG.setoffset!(rb, 10)
        @test AG.getoffset(rb) ≈ 10
        AG.setoffset!(rb, 0)
        @test AG.getoffset(rb) ≈ 0

        @test AG.getscale(rb) == 1
        AG.setscale!(rb, 0.5)
        @test AG.getscale(rb) ≈ 0.5
        AG.setscale!(rb, 2)
        @test AG.getscale(rb) ≈ 2
        AG.setscale!(rb, 1)
        @test AG.getscale(rb) ≈ 1

        @test AG.getnodatavalue(rb) ≈ -1e10
        AG.setnodatavalue!(rb, -100)
        @test AG.getnodatavalue(rb) ≈ -100
        AG.deletenodatavalue!(rb)
        @test AG.getnodatavalue(rb) ≈ -1e10

        AG.copy(dataset) do dest
            destband = AG.getband(dest, 1)
            AG.copywholeraster!(rb, destband)
            @test sprint(print, destband) == """
            [GA_Update] Band 1 (Gray): 100 x 100 (UInt8)
                blocksize: 100×81, nodata: -1.0e10, units: 1.0px + 0.0
                overviews: """
            @test AG.noverview(destband) == 0
            AG.buildoverviews!(dest, Cint[2, 4, 8])
            @test AG.noverview(destband) == 3
            @test sprint(print, destband) == """
            [GA_Update] Band 1 (Gray): 100 x 100 (UInt8)
                blocksize: 100×81, nodata: -1.0e10, units: 1.0px + 0.0
                overviews: (0) 50x50 (1) 25x25 (2) 13x13
                           """
            @test AG.getcolorinterp(destband) == GDAL.GCI_GrayIndex
            AG.setcolorinterp!(destband, GDAL.GCI_RedBand)
            @test AG.getcolorinterp(destband) == GDAL.GCI_RedBand

            @test sprint(print, AG.sampleoverview(destband, 100)) == """
            [GA_Update] Band 1 (Gray): 13 x 13 (UInt8)
                blocksize: 128×128, nodata: -1.0e10, units: 1.0px + 0.0
                overviews: """
            @test sprint(print, AG.sampleoverview(destband, 200)) == """
            [GA_Update] Band 1 (Gray): 25 x 25 (UInt8)
                blocksize: 128×128, nodata: -1.0e10, units: 1.0px + 0.0
                overviews: """
            @test sprint(print, AG.sampleoverview(destband, 500)) == """
            [GA_Update] Band 1 (Gray): 25 x 25 (UInt8)
                blocksize: 128×128, nodata: -1.0e10, units: 1.0px + 0.0
                overviews: """
            AG.sampleoverview(destband, 1000) do result
                @test sprint(print, result) == """
                [GA_Update] Band 1 (Gray): 50 x 50 (UInt8)
                    blocksize: 128×128, nodata: -1.0e10, units: 1.0px + 0.0
                    overviews: """
            end
            @test sprint(print, AG.getmaskband(destband)) == """
            [GA_ReadOnly] Band 0 (Undefined): 100 x 100 (UInt8)
                blocksize: 100×81, nodata: -1.0e10, units: 1.0px + 0.0
                overviews: """
            @test AG.maskflags(destband) == 1
            AG.createmaskband!(destband, 3)
            AG.getmaskband(destband) do maskband
                @test sprint(print, maskband) == """
                [GA_Update] Band 1 (Gray): 100 x 100 (UInt8)
                    blocksize: 100×81, nodata: -1.0e10, units: 1.0px + 0.0
                    overviews: """
            end
            @test AG.maskflags(destband) == 3
            AG.fillraster!(destband, 3)
            AG.setcategorynames!(destband, ["foo","bar"])
            @test AG.getcategorynames(destband) == ["foo", "bar"]

            AG.getoverview(destband, 0) do overview
                AG.regenerateoverviews!(destband, [
                    overview,
                    AG.getoverview(destband, 2)
                ])
            end

            AG.createcolortable(GDAL.GPI_RGB) do ct
                AG.createcolorramp!(ct,
                    128, GDAL.GDALColorEntry(0,0,0,0),
                    255, GDAL.GDALColorEntry(0,0,255,0)
                )
                AG.setcolortable!(destband, ct)
                @test AG.ncolorentry(ct) == 256
                AG.getcolortable(destband) do ct2
                    @test AG.ncolorentry(ct) == AG.ncolorentry(ct2)
                end
                AG.clearcolortable!(destband)

                AG.createRAT(ct) do rat
                    AG.setdefaultRAT!(destband, rat)
                    @test AG.getdefaultRAT(destband) != C_NULL
                end
            end
        end
    end
end

# untested
# setcategorynames!(rasterband, names)
# getcolortable(band) do C_NULL


text:
using Documenter, PolyhedralRelaxations
using DocumenterTools: Themes

Themes.compile(
    joinpath(@__DIR__, "src/assets/pr-light.scss"),
    joinpath(@__DIR__, "src/assets/themes/documenter-light.css"),
)

Themes.compile(
    joinpath(@__DIR__, "src/assets/pr-dark.scss"),
    joinpath(@__DIR__, "src/assets/themes/documenter-dark.css"),
)

makedocs(
    sitename = "PolyhedralRelaxations",
    format = Documenter.HTML(
        mathengine = Documenter.MathJax(),
        assets=[asset("https://fonts.googleapis.com/css?family=Montserrat|Source+Code+Pro&display=swap", class=:css)],
        prettyurls = get(ENV, "CI", nothing) == "true",
    ),
    strict = true,
    authors = "Sujeevraja Sanjeevi and Kaarthik Sundar",
    pages = [
        "Introduction" => "index.md",
        "API" => "api.md",
        "Interfacing with JuMP" => "interface.md",
        "Reference" => "reference.md",
    ],
)

deploydocs(repo = "github.com/sujeevraja/PolyhedralRelaxations.jl.git")