text:
# all the message handlers for AUV example

## middleware handler (callback) functions

function lbl_hdlr(channel::String, msgdata::pose_t, dfg::AbstractDFG, dashboard::Dict, lblDict::Dict, lbllog)
  msgtime = unix2datetime(msgdata.utime*1e-6)
  dashboard[:lastMsgTime] = msgtime

  # TODO -- use buffer instead (legacy)
  if length(lblDict) == 0
    lblDict[dashboard[:lastMsgTime]] = (msgdata.pos[1:2] |> deepcopy)
  end

  # add to the ciruclar buffer
  push!( dashboard[:lblBuffer], (msgtime, (msgdata.pos[1:2] |> deepcopy), Bool[false;]) )

  # also record in log for later printing
  println(lbllog, "$msgtime, $(msgdata.pos[1]),$(msgdata.pos[2])")

  # check if any lbl can be added as prior to factor graph
  for lbl in dashboard[:lblBuffer]
    # might have been added in the past
    lbl[3][1] ? continue : nothing
    # find nearest in time
    corr = findVariableNearTimestamp(dfg, lbl[1], r"x\d")  #tags=[:POSE;])
    # is there a match and its real close in time (assume 50Hz max rate)
    if 1 == length(corr) && abs(corr[1][2]) < Millisecond(10)
      # set the reference data in that variable
      setVariableRefence!(dfg, corr[1][1][1], reshape(lbl[2],:,1), refKey=:lbl)
      lbl[3][1] = true
    end
  end
  nothing
end

function mag_hdlr(channel::String, msgdata::pose_t, dfg::AbstractDFG, dashboard::Dict, magDict::Dict)
  msgtime = unix2datetime(msgdata.utime*1e-6)
  dashboard[:lastMsgTime] = msgtime
  ea = TU.convert(Euler, Quaternion(msgdata.orientation[1],msgdata.orientation[2:4]) )

  # add first element only
  # TODO use buffer instead (legacy)
  if length(magDict) == 0
    magDict[dashboard[:lastMsgTime]] = ea.Y
  end

  push!( dashboard[:magBuffer], (msgtime, ea.Y, Bool[false;]) )

  # check if any mag can be added as prior to factor graph
  for mag in dashboard[:magBuffer]
    # might have been added in the past
    mag[3][1] ? continue : nothing
    corr = findVariableNearTimestamp(dfg, mag[1], r"x\d")  #tags=[:POSE;])
    # is there a match and its real close in time (assume 50Hz max rate)
    if 1 == length(corr) && abs(corr[1][2]) < Millisecond(10)
      # double check that the variable does not already have a mag prior on it
      if !hasTagsNeighbors(dfg, corr[1][1][1], [:MAGNETOMETER;])
        # add the mag measurement as a prior
        @show "add mag"
        magPrior = PartialPriorYawPose2(Normal(mag[2], dashboard[:magNoise]))
        addFactor!(dfg, corr[1][1][1:1], magPrior, autoinit=false, solvable=1, tags=[:MAGNETOMETER;])
        # set flag that this has already been added to the factor graph
        mag[3][1] = true
      end
    end
  end

  nothing
end

function range_hdlr(channel::String, msgdata::raw_t, dfg::AbstractDFG, dashboard::Dict)
  msgtime = unix2datetime(msgdata.utime*1e-6)
  dashboard[:lastMsgTime] = msgtime
  # at some point there was a bug where range data was not in lcmlog as expected -- it was not understood, hence the assert
  @assert msgdata.length >= 130000

  # add msg to buffer based on RANGESTRIDE
  data = Array{Float64,2}(undef, 0,0)
  dashboard[:rangeCount] += 1
  if dashboard[:rangeCount] % dashboard[:RANGESTRIDE] == 0
    dashboard[:lastRangeMsg] = deepcopy(msgdata)
    data = reshape(reinterpret(Float64, msgdata.data), :, 2)
    # Tuple of msgtime, data, ifadded
    push!( dashboard[:rangesBuffer], (msgtime, data, Bool[false;]) )
    dashboard[:rangeCount] = 0
  end

  ## should a range message be added to the factor graph?
  # should be at least one
  if 0 < length(dashboard[:rangesBuffer])
    lastPose = getLastPoses(dfg, number=1)[1]
    lastT = getTimestamp(getVariable(dfg, lastPose))
    # check latest pose against all ranges in buffer
    rantim = map(x->(dashboard[:rangesBuffer][x][1],dashboard[:rangesBuffer][x][3][1]), 1:length(dashboard[:rangesBuffer]))
    @show ind, mdt, corrs = findClosestTimestamp([(lastT, lastPose);], rantim)
    # there should be only one correspondence
    @assert corrs <= 1
    # check if already added, add and set flag
    if corrs == 1 && mdt < Millisecond(100)
      rangeIdx = ind[2]
      if !dashboard[:rangesBuffer][rangeIdx][3][1]
        bera = Pose2Point2Range(AliasingScalarSampler(dashboard[:rangesBuffer][rangeIdx][2][:,1],dashboard[:rangesBuffer][rangeIdx][2][:,2], SNRfloor=dashboard[:SNRfloor]))
        addFactor!(dfg, [lastPose; :l1], bera, solvable=0, autoinit=false)
        fSym = lastPose*"l1f1"
        setTimestamp!(getFactor(dfg, fSym), msgtime)
        dashboard[:rangesBuffer][rangeIdx][3][1] = true
        # check that ranges are not lost owing to clock drift
        deltaClk = lastT - dashboard[:rangesBuffer][rangeIdx][1]
        @assert abs(deltaClk) < Millisecond(100)

        # set the factors as solvable (landmark if not there yet too)
        addLm = getSolvable(dfg, :l1) == 1 ? Symbol[] : [:l1; :l1f1]
        put!(dashboard[:solveSettings].solvables, [fSym;addLm])
      end
    end
  end
  nothing
end

function pose_hdlr(channel::String,
                   msgdata::pose_t,
                   dfg::AbstractDFG,
                   dashboard::Dict,
                   timinglog,
                   drtlog,
                   drolog,
                   dirodolog,
                   rawodolog,
                   alltht )
  #
  global ODOSCALE
  t0 = time_ns()
  odoT = unix2datetime(msgdata.utime*1e-6)
  dashboard[:lastMsgTime] = odoT # unix2datetime(msgdata.utime*1e-6)
  # "accumulateDiscreteLocalFrame! on all RTTs" |> println
  scodo = ODOSCALE.*msgdata.pos[1:3]
  for (vsym, drt) in dashboard[:drtMpp]
    for symFct in intersect(ls(dfg, MutablePose2Pose2Gaussian), ls(dfg,vsym))
      drtFct = getFactorType(dfg, symFct)
      accumulateDiscreteLocalFrame!(drtFct, scodo, dashboard[:Qc_odo], 1.0, Phik=SE2(scodo)) # drt
      println(alltht, "$odoT, $(length(dashboard[:drtMpp])), $vsym, $(symFct), $(scodo[1]), $(scodo[2]), $(scodo[3]), $(drtFct.Zij.μ[1]), $(drtFct.Zij.μ[2]), $(drtFct.Zij.μ[3])")
    end
  end
  dt_acc = (time_ns()-t0)/1e9

  # get message time

  # write the latest DRT solution to file
  drtFnc = string(dashboard[:solveSettings].drtCurrent[1], dashboard[:solveSettings].drtCurrent[2], "f1") |> Symbol
  val = accumulateFactorMeans(dfg, [drtFnc;])
  drtFncMu = getFactorType(dfg, drtFnc).Zij.μ
  println(drtlog, "$odoT, $(drtFnc), $(dashboard[:lastPose]), $(val[1]), $(val[2]), $(val[3]), $(scodo[1]), $(scodo[2]), $(scodo[3]), $(collect(keys(dashboard[:drtMpp]))), $(drtFncMu)")
  drval = accumulateFactorMeans(dfg, [:x0drt_0f1;])
  println(drolog, "$odoT, $(dashboard[:lastPose]), $(drval[1]), $(drval[2]), $(drval[3])")
  dt_drt = (time_ns()-t0)/1e9

  # check if a new pose and factor must be added?
  if (dashboard[:odoTime] + dashboard[:poseRate]) <= odoT
    # update odo time
    dashboard[:odoTime] = odoT

    # get factor to next pose
    nPose = nextPose(dashboard[:lastPose])
    @info "pose_hdlr, adding new pose $nPose, at $odoT.  DRTs $(collect(keys(dashboard[:drtMpp]))), drtCurrent=$(dashboard[:solveSettings].drtCurrent)"

    # get drt from last pose
    lastPoseDrt = intersect(ls(dfg, dashboard[:lastPose]), ls(dfg, MutablePose2Pose2Gaussian))[1]
    drtLast = getFactorType(dfg, lastPoseDrt)
    # drtLast = dashboard[:drtMpp][dashboard[:lastPose]]

    # add new variable and factor to the graph
    fctsym = duplicateToStandardFactorVariable(Pose2Pose2,
                                               drtLast,
                                               dfg,
                                               dashboard[:lastPose],
                                               nPose,
                                               solvable=0,
                                               autoinit=false,
                                               cov=dashboard[:odoCov]  )
    #
    # set timestamp to msg times
    # setTimestamp!(getVariable(dfg, nPose), odoT)
    setTimestamp!(dfg, nPose, odoT)
    # delete drt unless used by real time prediction
    # TODO assuming stride ends on a 0
    poseStrideTail = sortDFG(ls(dfg, r"x\d+9\b|x9\b", solvable=0))
    poseStrideTail = 3 < length(poseStrideTail) ? poseStrideTail[end-3:end] : poseStrideTail
    for vkey in keys(dashboard[:drtMpp])
      if vkey != dashboard[:lastPose] && vkey != :x0 && !(vkey in poseStrideTail) # dashboard[:drtCurrent][1]
      # if dashboard[:lastPose] != :x0 && !(dashboard[:lastPose] in poseStrideTail) # dashboard[:drtCurrent][1]
        delete!(dashboard[:drtMpp], vkey)
        # delete!(dashboard[:drtMpp], dashboard[:lastPose])
      end
    end

    # update last pose to new pose
    dashboard[:lastPose] = nPose

    # create new drt on new variable
    nDrt = Symbol(string("drt_",string(nPose)[2:end]))
    addVariable!(dfg, nDrt, Pose2, solvable=0)
    drec1 = MutablePose2Pose2Gaussian(MvNormal(zeros(3), Matrix{Float64}(LinearAlgebra.I, 3,3)))
    addFactor!(dfg, [nPose; nDrt], drec1, solvable=0, autoinit=false)
    dashboard[:drtMpp][nPose] = drec1
    put!(dashboard[:solveSettings].solvables, [nPose; fctsym])
    dashboard[:poseStride] += 1

    if dashboard[:doDelay] && 10 <= dashboard[:poseStride]
      # how far to escalate -- :canTakePoses is de-escalated by solver in manageSolveTree()
      # if dashboard[:solveSettings].canTakePoses == HSMHandling
      #   dashboard[:solveSettings].canTakePoses = HSMOverlapHandling
      # elseif dashboard[:solveSettings].canTakePoses == HSMOverlapHandling
      #   dashboard[:solveSettings].canTakePoses = HSMBlocking
      # end
      @info "new solve token $nPose"
      put!(dashboard[:solveSettings].poseSolveToken, nPose)
      dashboard[:poseStride] = 0
    end

    # also store current factor graph if part of plotting request
    if dashboard[:savePlotting]
      dfg_ = deepcopy(dfg)
      @spawn saveDFG(dfg_, joinpath(getLogPath(dfg_), "fg_$nPose"))
    end
  end
  dt_pose = (time_ns()-t0)/1e9

  ## Separate odometry check
  # mutable pose2pose2 factor holding odometry only reference
  # mpp = dashboard[:drtOdoRef][3]
  # accumulate odo
  accumulateDiscreteLocalFrame!(getFactorType(dfg, :x0drt_reff1), scodo[1:3], dashboard[:Qc_odo], 1.0, Phik=SE2(scodo[1:3]))
  # write result to file
  drval = accumulateFactorMeans(dfg, [:x0drt_reff1;])
  println(dirodolog, "$odoT, $(dashboard[:lastPose]), $(drval[1]), $(drval[2]), $(drval[3])")
  println(rawodolog, "$odoT, $(dashboard[:lastPose]), $(scodo[1]), $(scodo[2]), $(scodo[3])")

  dt_total = (time_ns()-t0)/1e9
  println(timinglog, "$odoT, $(dashboard[:lastPose]), $t0, $dt_acc, $dt_drt, $dt_pose, $dt_total")

  nothing
end

#


text:
module TestDoctest

using StaticStorages
using Documenter: doctest
using Test

test_doctest() = doctest(StaticStorages, manual = false)

end  # module

ext: jl


text:
using Compose
import Colors
typealias ComposeColor Union{Colors.Color,
                             Colors.AlphaColor,
                             Colors.AbstractString}

"""
Given an adjacency matrix and two vectors of X and Y coordinates, returns
a Compose tree of the graph layout

Arguments:
    adj_matrix       Adjacency matrix of some type. Non-zero of the eltype
                     of the matrix is used to determine if a link exists,
                     but currently no sense of magnitude
    locs_x, locs_y   Locations of the nodes. Can be any units you want,
                     but will be normalized and centered anyway
    labels           Optional. Labels for the vertices. Default: Any[]
    nodefillc        Color to fill the nodes with. Default: #AAAAFF
    nodestrokec      Color for the nodes stroke. Default: #BBBBBB
    edgestrokec      Color for the edge strokes. Default: #BBBBBB
    arrowlengthfrac  Fraction of line length to use for arrows.
                     Set to 0 for no arrows. Default: 0.1
    angleoffset      angular width in radians for the arrows. Default: π/9 (20 degrees).

"""
function compose_layout_adj{S, T<:Real}(
    adj_matrix::Array{S,2},
    locs_x::Vector{T}, locs_y::Vector{T};
    labels::Vector=Any[],
    filename::AbstractString="",
    labelc::ComposeColor="#000000",
    nodefillc::ComposeColor="#AAAAFF",
    nodestrokec::ComposeColor="#BBBBBB",
    edgestrokec::ComposeColor="#BBBBBB",
    labelsize::Real=4.0,
    arrowlengthfrac::Real=0.1,
    angleoffset=20.0/180.0*π)

    length(locs_x) != length(locs_y) && error("Vectors must be same length")
    const N = length(locs_x)
    if length(labels) != N && length(labels) != 0
        error("Must have one label per node (or none)")
    end

    # Scale to unit square
    min_x, max_x = minimum(locs_x), maximum(locs_x)
    min_y, max_y = minimum(locs_y), maximum(locs_y)
    function scaler(z, a, b)
        2.0*((z - a)/(b - a)) - 1.0
    end
    map!(z -> scaler(z, min_x, max_x), locs_x)
    map!(z -> scaler(z, min_y, max_y), locs_y)

    # Determine sizes
    const NODESIZE    = 0.25/sqrt(N)
    const LINEWIDTH   = 3.0/sqrt(N)
    const ARROWLENGTH = LINEWIDTH * arrowlengthfrac

    # Create lines and arrow heads
    lines = Any[]
    for i = 1:N
        for j = 1:N
            i == j && continue
            if adj_matrix[i,j] != zero(eltype(adj_matrix))
                push!(lines, lineij(locs_x, locs_y, i,j, NODESIZE, ARROWLENGTH, angleoffset))
            end
        end
    end

    # Create nodes
    nodes = [circle(locs_x[i],locs_y[i],NODESIZE) for i=1:N]

    # Create labels (if wanted)
    texts = length(labels) == N ?
        [text(locs_x[i],locs_y[i],labels[i],hcenter,vcenter) for i=1:N] : Any[]

    compose(context(units=UnitBox(-1.2,-1.2,+2.4,+2.4)),
        compose(context(), texts..., fill(labelc), stroke(nothing), fontsize(labelsize)),
        compose(context(), nodes..., fill(nodefillc), stroke(nodestrokec)),
        compose(context(), lines..., stroke(edgestrokec), linewidth(LINEWIDTH))
    )
end

function arrowcoords(θ, endx, endy, arrowlength, angleoffset=20.0/180.0*π)
    arr1x = endx - arrowlength*cos(θ+angleoffset)
    arr1y = endy - arrowlength*sin(θ+angleoffset)
    arr2x = endx - arrowlength*cos(θ-angleoffset)
    arr2y = endy - arrowlength*sin(θ-angleoffset)
    return (arr1x, arr1y), (arr2x, arr2y)
end

function lineij(locs_x, locs_y, i, j, NODESIZE, ARROWLENGTH, angleoffset)
    Δx = locs_x[j] - locs_x[i]
    Δy = locs_y[j] - locs_y[i]
    d  = sqrt(Δx^2 + Δy^2)
    θ  = atan2(Δy,Δx)
    endx  = locs_x[i] + (d-NODESIZE)*1.00*cos(θ)
    endy  = locs_y[i] + (d-NODESIZE)*1.00*sin(θ)
    if ARROWLENGTH > 0.0
        arr1, arr2 = arrowcoords(θ, endx, endy, ARROWLENGTH, angleoffset)
        composenode = compose(
                context(),
                line([(locs_x[i], locs_y[i]), (endx, endy)]),
                line([arr1, (endx, endy)]),
                line([arr2, (endx, endy)])
            )
    else
        composenode = compose(
                context(),
                line([(locs_x[i], locs_y[i]), (endx, endy)])
            )
    end
    return composenode
end

"""
Given an adjacency matrix and two vectors of X and Y coordinates, returns
an SVG of the graph layout.

Arguments:
    adj_matrix       Adjacency matrix of some type. Non-zero of the eltype
                     of the matrix is used to determine if a link exists,
                     but currently no sense of magnitude
    locs_x, locs_y   Locations of the nodes. Can be any units you want,
                     but will be normalized and centered anyway
    labels           Optional. Labels for the vertices. Default: Any[]
    filename         Optional. Output filename for SVG. If blank, just
                     tries to draw it anyway, which will display in IJulia
    nodefillc        Color to fill the nodes with. Default: #AAAAFF
    nodestrokec      Color for the nodes stroke. Default: #BBBBBB
    edgestrokec      Color for the edge strokes. Default: #BBBBBB
    arrowlengthfrac  Fraction of line length to use for arrows.
                     Set to 0 for no arrows. Default: 0.1
    angleoffset      angular width in radians for the arrows. Default: π/9 (20 degrees).
"""
function draw_layout_adj{S, T<:Real}(
    adj_matrix::Array{S,2},
    locs_x::Vector{T}, locs_y::Vector{T};
    labels::Vector=Any[],
    filename::AbstractString="",
    labelc::ComposeColor="#000000",
    nodefillc::ComposeColor="#AAAAFF",
    nodestrokec::ComposeColor="#BBBBBB",
    edgestrokec::ComposeColor="#BBBBBB",
    labelsize::Real=4.0,
    arrowlengthfrac::Real=0.1,
    angleoffset=20.0/180.0*π)

    draw(filename == "" ? SVG(8inch, 8inch) : SVG(filename, 8inch, 8inch),
        compose_layout_adj(adj_matrix, locs_x, locs_y, labels=labels,
            labelc=labelc, nodefillc=nodefillc, nodestrokec=nodestrokec,
            edgestrokec=edgestrokec, labelsize=labelsize,
            arrowlengthfrac=arrowlengthfrac, angleoffset=angleoffset)
    )
end


text:
N = 10
D1 = speye(N-1,N) - spdiagm(ones(N-1),1,N-1,N)
D = [ kron(D1, speye(N)); kron(speye(N), D1) ]
i, j = N*1 + 2, N*7+7
b = zeros(N^2); b[i], b[j] = 1, -1
v = (D' * D) \ b
v[i] - v[j]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
immutable ForwardPass
    net::ForwardNet
    input::Array{Symbol}
    output::Array{Symbol}
    activation_order::Vector{Int} # in order
end
function Base.show(io::IO, forwardpass::ForwardPass)
    net = forwardpass.net

    if isempty(forwardpass.activation_order)
        print(io, "ForwardPass: (empty) ", forwardpass.input, " → ", forwardpass.output)
    else
        print(io, "ForwardPass: ", name(net[forwardpass.activation_order[1]]))
        for i in 2 : length(forwardpass.activation_order)
            print(io, " → ", name(net[forwardpass.activation_order[i]]))
        end
    end
end
Base.show(forwardpass::ForwardPass) = print_forward_pass(STDOUT, forwardpass)

function calc_forwardpass(net::ForwardNet, input::Array{Symbol}, output::Array{Symbol})
    #=
    1 - get a topological sort of net
    2 - run through nodes in topologogical order and activate those that are:
            # have at least one ancestor in input and one descendent in output
    3 - will need to run forward! on all active nodes, in topological order
    =#

    input_indeces = Set{Int}()
    input_dijkstra = Dict{Int, LightGraphs.DijkstraState{Int}}()
    for name in input
        index = indexof(net, name)
        push!(input_indeces, index)
        input_dijkstra[index] = dijkstra_shortest_paths(net.dag, index)
    end

    output_indeces = Set{Int}()
    for name in output
        push!(output_indeces, indexof(net, name))
    end

    activation_order = Int[]
    for i in topological_sort_by_dfs(net.dag)

        if !isa(net[i], Variable)

            add_it = false

            for parent in input_indeces

                if parent == i
                    add_it = true
                elseif !in(i, output_indeces) && input_dijkstra[parent].dists[i] != typemax(Int) # is descendent of parent
                    dijkstra = dijkstra_shortest_paths(net.dag, i)

                    for child in output_indeces
                        if dijkstra.dists[child] != typemax(Int)
                            # child is a descendent of our node
                            add_it = true
                        end
                    end

                end
            end

            if add_it
                push!(activation_order, i)
            end
        end
    end

    ForwardPass(net, input, output, activation_order)
end
calc_forwardpass(net::ForwardPass, input::Symbol, output::Array{Symbol}) = calc_forwardpass(net, [input], output)
calc_forwardpass(net::ForwardPass, input::Array{Symbol}, output::Symbol) = calc_forwardpass(net, input, [output])
calc_forwardpass(net::ForwardPass, input::Symbol,        output::Symbol) = calc_forwardpass(net, [input], [output])

function forward!(forwardpass::ForwardPass)
    for index in forwardpass.activation_order
        forward!(forwardpass.net.nodes[index])
    end
    forwardpass
end


text:
@testset "Estimators" begin
  # floating point tolerance
  tol = 100eps()

  # create some data
  dim = 3; nobs = 10
  pset = PointSet(10*rand(dim, nobs))
  data = georef((z=rand(nobs),), pset)

  γ = GaussianVariogram(sill=1., range=1., nugget=0.)
  simkrig = SimpleKriging(data, γ, mean(data[:z]))
  ordkrig = OrdinaryKriging(data, γ)
  unikrig = UniversalKriging(data, γ, 1)
  drikrig = ExternalDriftKriging(data, γ, [x->1.])

  # Kriging is an interpolator
  for j in 1:nobs
    SKestimate, SKvar = predict(simkrig, :z, pset[j])
    OKestimate, OKvar = predict(ordkrig, :z, pset[j])
    UKestimate, UKvar = predict(unikrig, :z, pset[j])
    DKestimate, DKvar = predict(drikrig, :z, pset[j])

    # estimate checks
    @test SKestimate ≈ data[:z][j]
    @test OKestimate ≈ data[:z][j]
    @test UKestimate ≈ data[:z][j]
    @test DKestimate ≈ data[:z][j]

    # variance checks
    @test SKvar + tol ≥ 0
    @test OKvar + tol ≥ 0
    @test UKvar + tol ≥ 0
    @test DKvar + tol ≥ 0
    @test SKvar ≤ OKvar + tol
  end

  # save results on a particular location pₒ
  pₒ = rand(Point3)
  SKestimate, SKvar = predict(simkrig, :z, pₒ)
  OKestimate, OKvar = predict(ordkrig, :z, pₒ)
  UKestimate, UKvar = predict(unikrig, :z, pₒ)
  DKestimate, DKvar = predict(drikrig, :z, pₒ)

  # Kriging is translation-invariant
  h = rand(Vec3)
  pset_h = PointSet([pset[i] + h for i in 1:nelements(pset)])
  data_h = georef((z=data[:z],), pset_h)
  simkrig_h = SimpleKriging(data_h, γ, mean(data_h[:z]))
  ordkrig_h = OrdinaryKriging(data_h, γ)
  unikrig_h = UniversalKriging(data_h, γ, 1)
  drikrig_h = ExternalDriftKriging(data_h, γ, [x->1.])
  SKestimate_h, SKvar_h = predict(simkrig_h, :z, pₒ + h)
  OKestimate_h, OKvar_h = predict(ordkrig_h, :z, pₒ + h)
  UKestimate_h, UKvar_h = predict(unikrig_h, :z, pₒ + h)
  DKestimate_h, DKvar_h = predict(drikrig_h, :z, pₒ + h)
  @test SKestimate_h ≈ SKestimate
  @test SKvar_h ≈ SKvar
  @test OKestimate_h ≈ OKestimate
  @test OKvar_h ≈ OKvar
  @test UKestimate_h ≈ UKestimate
  @test UKvar_h ≈ UKvar
  @test DKestimate_h ≈ DKestimate
  @test DKvar_h ≈ DKvar

  # Kriging estimate is invariant under covariance scaling
  # Kriging variance is multiplied by the same factor
  α = rand()
  γ_α = GaussianVariogram(sill=α, range=1., nugget=0.)
  simkrig_α = SimpleKriging(data, γ_α, mean(data[:z]))
  ordkrig_α = OrdinaryKriging(data, γ_α)
  unikrig_α = UniversalKriging(data, γ_α, 1)
  drikrig_α = ExternalDriftKriging(data, γ_α, [x->1.])
  SKestimate_α, SKvar_α = predict(simkrig_α, :z, pₒ)
  OKestimate_α, OKvar_α = predict(ordkrig_α, :z, pₒ)
  UKestimate_α, UKvar_α = predict(unikrig_α, :z, pₒ)
  DKestimate_α, DKvar_α = predict(drikrig_α, :z, pₒ)
  @test SKestimate_α ≈ SKestimate
  @test SKvar_α ≈ α*SKvar
  @test OKestimate_α ≈ OKestimate
  @test OKvar_α ≈ α*OKvar
  @test UKestimate_α ≈ UKestimate
  @test UKvar_α ≈ α*UKvar
  @test DKestimate_α ≈ DKestimate
  @test DKvar_α ≈ α*DKvar

  # Kriging variance is a function of data configuration, not data values
  δ = rand(nobs)
  data_δ = georef((z=data[:z].+δ,), pset)
  simkrig_δ = SimpleKriging(data_δ, γ, mean(data_δ[:z]))
  ordkrig_δ = OrdinaryKriging(data_δ, γ)
  unikrig_δ = UniversalKriging(data_δ, γ, 1)
  drikrig_δ = ExternalDriftKriging(data_δ, γ, [x->1.])
  SKestimate_δ, SKvar_δ = predict(simkrig_δ, :z, pₒ)
  OKestimate_δ, OKvar_δ = predict(ordkrig_δ, :z, pₒ)
  UKestimate_δ, UKvar_δ = predict(unikrig_δ, :z, pₒ)
  DKestimate_δ, DKvar_δ = predict(drikrig_δ, :z, pₒ)
  @test SKvar_δ ≈ SKvar
  @test OKvar_δ ≈ OKvar
  @test UKvar_δ ≈ UKvar
  @test DKvar_δ ≈ DKvar

  # Ordinary Kriging ≡ Universal Kriging with 0th degree drift
  unikrig_0th = UniversalKriging(data, γ, 0)
  OKestimate, OKvar = predict(ordkrig, :z, pₒ)
  UKestimate, UKvar = predict(unikrig_0th, :z, pₒ)
  @test OKestimate ≈ UKestimate
  @test OKvar ≈ UKvar

  # Ordinary Kriging ≡ Kriging with constant external drift
  driftkrig_const = ExternalDriftKriging(data, γ, [x->1.])
  OKestimate, OKvar = predict(ordkrig, :z, pₒ)
  DKestimate, DKvar = predict(driftkrig_const, :z, pₒ)
  @test OKestimate ≈ DKestimate
  @test OKvar ≈ DKvar

  # Non-stationary variograms are allowed
  γ_ns = PowerVariogram()
  ordkrig_ns = OrdinaryKriging(data, γ_ns)
  unikrig_ns = UniversalKriging(data, γ_ns, 1)
  drikrig_ns = ExternalDriftKriging(data, γ_ns, [x->1.])
  for j in 1:nobs
    OKestimate, OKvar = predict(ordkrig_ns, :z, pset[j])
    UKestimate, UKvar = predict(unikrig_ns, :z, pset[j])
    DKestimate, DKvar = predict(drikrig_ns, :z, pset[j])

    # estimate checks
    @test OKestimate ≈ data[:z][j]
    @test UKestimate ≈ data[:z][j]
    @test DKestimate ≈ data[:z][j]

    # variance checks
    @test OKvar + tol ≥ 0
    @test UKvar + tol ≥ 0
    @test DKvar + tol ≥ 0
  end

  # Floating point precision checks
  X_f         = rand(Float32, dim, nobs)
  z_f         = rand(Float32, nobs)
  X_d         = Float64.(X_f)
  z_d         = Float64.(z_f)
  pset_f      = PointSet(X_f)
  data_f      = georef((z=z_f,), pset_f)
  pset_d      = PointSet(X_d)
  data_d      = georef((z=z_d,), pset_d)
  pₒ_f        = rand(Point{dim,Float32})
  pₒ_d        = convert(Point{dim,Float64}, pₒ_f)
  γ_f         = GaussianVariogram(sill=1f0, range=1f0, nugget=0f0)
  simkrig_f   = SimpleKriging(data_f, γ_f, mean(data_f[:z]))
  ordkrig_f   = OrdinaryKriging(data_f, γ_f)
  unikrig_f   = UniversalKriging(data_f, γ_f, 1)
  drikrig_f   = ExternalDriftKriging(data_f, γ_f, [x->1f0])
  γ_d         = GaussianVariogram(sill=1., range=1., nugget=0.)
  simkrig_d   = SimpleKriging(data_d, γ_d, mean(data_d[:z]))
  ordkrig_d   = OrdinaryKriging(data_d, γ_d)
  unikrig_d   = UniversalKriging(data_d, γ_d, 1)
  drikrig_d = ExternalDriftKriging(data_d, γ_d, [x->1.])
  SKestimate_f, SKvar_f = predict(simkrig_f, :z, pₒ_f)
  OKestimate_f, OKvar_f = predict(ordkrig_f, :z, pₒ_f)
  UKestimate_f, UKvar_f = predict(unikrig_f, :z, pₒ_f)
  DKestimate_f, DKvar_f = predict(drikrig_f, :z, pₒ_f)
  SKestimate_d, SKvar_d = predict(simkrig_d, :z, pₒ_d)
  OKestimate_d, OKvar_d = predict(ordkrig_d, :z, pₒ_d)
  UKestimate_d, UKvar_d = predict(unikrig_d, :z, pₒ_d)
  DKestimate_d, DKvar_d = predict(drikrig_d, :z, pₒ_d)
  @test isapprox(SKestimate_f, SKestimate_d, atol=1e-4)
  @test isapprox(SKvar_f, SKvar_d, atol=1e-4)
  @test isapprox(OKestimate_f, OKestimate_d, atol=1e-4)
  @test isapprox(OKvar_f, OKvar_d, atol=1e-4)
  @test isapprox(UKestimate_f, UKestimate_d, atol=1e-4)
  @test isapprox(UKvar_f, UKvar_d, atol=1e-4)
  @test isapprox(DKestimate_f, DKestimate_d, atol=1e-4)
  @test isapprox(DKvar_f, DKvar_d, atol=1e-4)

  # ------------------
  # CHANGE OF SUPPORT
  # ------------------

  # create some data
  dim = 2; nobs = 10
  pset = PointSet(10*rand(dim, nobs))
  data = georef((z=rand(nobs),), pset)

  # basic estimators
  γ = GaussianVariogram(sill=1., range=1., nugget=0.)
  simkrig = SimpleKriging(data, γ, mean(data[:z]))
  ordkrig = OrdinaryKriging(data, γ)
  unikrig = UniversalKriging(data, γ, 1)
  drikrig = ExternalDriftKriging(data, γ, [x->1.])

  # prediction on a quadrangle
  uₒ = Quadrangle((0.,0.), (1.,0.), (1.,1.), (0.,1.))
  _, SKvar = predict(simkrig, :z, uₒ)
  _, OKvar = predict(ordkrig, :z, uₒ)
  _, UKvar = predict(unikrig, :z, uₒ)
  _, DKvar = predict(drikrig, :z, uₒ)

  # variance checks
  @test SKvar + tol ≥ 0
  @test OKvar + tol ≥ 0
  @test UKvar + tol ≥ 0
  @test DKvar + tol ≥ 0
  @test SKvar ≤ OKvar + tol
end


text:
module PatternFolds

import Base.length, Base.rand

# exports
export PatternFold, VectorFold
export pattern, gap, folds
export length, unfold

# includes
include("common.jl")
include("vector.jl")

end


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

# Instant arithmetic
(+)(x::Instant) = x
(-)(x::T, y::T) where {T<:Instant} = x.periods - y.periods

# TimeType arithmetic
(+)(x::TimeType) = x
(-)(x::T, y::T) where {T<:TimeType} = x.instant - y.instant

# Date-Time arithmetic
"""
    dt::Date + t::Time -> DateTime

The addition of a `Date` with a `Time` produces a `DateTime`. The hour, minute, second, and millisecond parts of
the `Time` are used along with the year, month, and day of the `Date` to create the new `DateTime`.
Non-zero microseconds or nanoseconds in the `Time` type will result in an `InexactError` being thrown.
"""
(+)(dt::Date, t::Time) = DateTime(dt ,t)
(+)(t::Time, dt::Date) = DateTime(dt, t)

# TimeType-Year arithmetic
function (+)(dt::DateTime, y::Year)
    oy, m, d = yearmonthday(dt); ny = oy + value(y); ld = daysinmonth(ny, m)
    return DateTime(ny, m, d <= ld ? d : ld, hour(dt), minute(dt), second(dt), millisecond(dt))
end
function (+)(dt::Date,y::Year)
    oy, m, d = yearmonthday(dt); ny = oy + value(y); ld = daysinmonth(ny, m)
    return Date(ny, m, d <= ld ? d : ld)
end
function (-)(dt::DateTime,y::Year)
    oy, m, d = yearmonthday(dt); ny = oy - value(y); ld = daysinmonth(ny, m)
    return DateTime(ny, m, d <= ld ? d : ld, hour(dt), minute(dt), second(dt), millisecond(dt))
end
function (-)(dt::Date,y::Year)
    oy, m, d = yearmonthday(dt); ny = oy - value(y); ld = daysinmonth(ny, m)
    return Date(ny, m, d <= ld ? d : ld)
end

# TimeType-Month arithmetic
# monthwrap adds two months with wraparound behavior (i.e. 12 + 1 == 1)
monthwrap(m1, m2) = (v = mod1(m1 + m2, 12); return v < 0 ? 12 + v : v)
# yearwrap takes a starting year/month and a month to add and returns
# the resulting year with wraparound behavior (i.e. 2000-12 + 1 == 2001)
yearwrap(y, m1, m2) = y + fld(m1 + m2 - 1, 12)

function (+)(dt::DateTime, z::Month)
    y,m,d = yearmonthday(dt)
    ny = yearwrap(y, m, value(z))
    mm = monthwrap(m, value(z)); ld = daysinmonth(ny, mm)
    return DateTime(ny, mm, d <= ld ? d : ld, hour(dt), minute(dt), second(dt), millisecond(dt))
end

function (+)(dt::Date, z::Month)
    y,m,d = yearmonthday(dt)
    ny = yearwrap(y, m, value(z))
    mm = monthwrap(m, value(z)); ld = daysinmonth(ny, mm)
    return Date(ny, mm, d <= ld ? d : ld)
end
function (-)(dt::DateTime, z::Month)
    y,m,d = yearmonthday(dt)
    ny = yearwrap(y, m, -value(z))
    mm = monthwrap(m, -value(z)); ld = daysinmonth(ny, mm)
    return DateTime(ny, mm, d <= ld ? d : ld, hour(dt), minute(dt), second(dt), millisecond(dt))
end
function (-)(dt::Date, z::Month)
    y,m,d = yearmonthday(dt)
    ny = yearwrap(y, m, -value(z))
    mm = monthwrap(m, -value(z)); ld = daysinmonth(ny, mm)
    return Date(ny, mm, d <= ld ? d : ld)
end

(+)(x::Date, y::Quarter) = x + Month(y)
(-)(x::Date, y::Quarter) = x - Month(y)
(+)(x::DateTime, y::Quarter) = x + Month(y)
(-)(x::DateTime, y::Quarter) = x - Month(y)
(+)(x::Date, y::Week) = return Date(UTD(value(x) + 7 * value(y)))
(-)(x::Date, y::Week) = return Date(UTD(value(x) - 7 * value(y)))
(+)(x::Date, y::Day)  = return Date(UTD(value(x) + value(y)))
(-)(x::Date, y::Day)  = return Date(UTD(value(x) - value(y)))
(+)(x::DateTime, y::Period) = return DateTime(UTM(value(x) + toms(y)))
(-)(x::DateTime, y::Period) = return DateTime(UTM(value(x) - toms(y)))
(+)(x::Time, y::TimePeriod) = return Time(Nanosecond(value(x) + tons(y)))
(-)(x::Time, y::TimePeriod) = return Time(Nanosecond(value(x) - tons(y)))
(+)(y::Period, x::TimeType) = x + y

# AbstractArray{TimeType}, AbstractArray{TimeType}
(-)(x::OrdinalRange{T}, y::OrdinalRange{T}) where {T<:TimeType} = Vector(x) - Vector(y)
(-)(x::AbstractRange{T}, y::AbstractRange{T}) where {T<:TimeType} = Vector(x) - Vector(y)

# Allow dates, times, and time zones to broadcast as unwrapped scalars
Base.Broadcast.broadcastable(x::AbstractTime) = Ref(x)
Base.Broadcast.broadcastable(x::TimeZone) = Ref(x)


text:
"""
    StudentsTLikelihood(σ²)

Student's T likelihood with `ν` degrees of freedom and variance `σ²` variance.
This is to be used if we assume that the uncertainity associated with the data
follows a Student's T distribution.

```math
    p(y|f) = StudentsT(y | f, σ², ν)
```
On calling, this would return a Student's T distribution with mean `f`, variance σ²
and `ν` degrees of freedom.
"""
struct StudentsTLikelihood{T<:Real}
    ν::Int
    σ²::T
end

@functor StudentsTLikelihood (σ²,)

(l::StudentsTLikelihood)(f::Real) = LocationScale(f, sqrt(l.σ²), TDist(l.ν))

(l::StudentsTLikelihood)(fs::AbstractVector{<:Real}) = Product(LocationScale.(fs, sqrt(l.σ²), TDist(l.ν)))


text:
# Statistical model for analyzing Radial Velocity observations
# Assumes existance of data in arrays times, obs, sigma_obs
# Assumes existance of functions: loglikelihood(theta), logprior_noise(theta), and is_valid_noise(theta)

export logprior_planets, logprior, logtarget, negative_logtarget
export set_times, set_obs, set_sigma_obs, num_noise_param

# Observational data to compare model to
global times # ::Array{Float64,1}
global obs # ::Array{Float64,1}
global sigma_obs #::Array{Float64,1}

# Functions to set global data within module
function set_times(t::Array{Float64,1}) global times = t   end
function set_obs(o::Array{Float64,1})  global obs = o   end
function set_sigma_obs(so::Array{Float64,1}) global sigma_obs = so  end

function is_valid(p::Vector)
   is_valid_noise(p) && is_valid_planets(p)
end

function logprior_period(P::Real)
  const min_period = 1.0
  const max_period = 10000.0
  const norm = -log(max_period/min_period)
  if !(min_period<=P<=max_period)
     return -Inf
  end
  return norm-log(P)
end

function logprior_amplitude(K::Real)
  const min_K = 0.0
  const max_K = 999.0
  const K0 = 1.0
  const norm = -log(K0*log1p(maxK/min_K))
  if !(min_K<=K<=max_K)
     return -Inf
  end
  return norm-log1p(K/K0)
end

function logprior_planets(p::Vector)
  if !is_valid_planets(p) return -Inf end  # prempt model evaluation
  num_pl = num_planets(p)
  logp = zero(eltype(p))
  if num_pl <= 0 return logp end
  if num_param_per_planet == 3
     logp -= num_pl*log(2pi)
  elseif num_param_per_planet == 5
     logp -= 2*num_pl*log(2pi)
  end
  for plid in 1:num_pl
    P::eltype(p) = extract_period(p,plid=plid)
    K::eltype(p) = extract_amplitude(p,plid=plid)
    logp += logprior_period(P) + logprior_amplitude(K)
  end
  return logp::eltype(p)
end

function logprior_offset(C::Real)
  const max_C = 1000.0
  if !(-max_C<=C<=max_C)
     return -Inf
  end
  return -log(2*max_C)
end

logprior_offset(p::Vector) = logprior_offset(extract_rvoffset(p))

function logprior(p::Vector)
  logprior_noise(p) + logprior_offset(p) + logprior_planets(p)
end

function logtarget(p::Vector)
  val = logprior(p)
  if val==-Inf
     return val
  end
  val += loglikelihood(p)
  return val
end

negative_logtarget(p::Vector) = -logtarget(p)


text:
using LinearAlgebra, CUDA, SparseArrays

export MatFacModel

mutable struct MatFacModel

    X::AbstractMatrix            # KxM "instance factor" matrix
    Y::AbstractMatrix            # KxN "feature factor" matrix

    instance_reg_mats::AbstractVector{AbstractMatrix}  # K x (M x M) (usually sparse)
    feature_reg_mats::AbstractVector{AbstractMatrix}  # K x (N x N) (usually sparse)

    instance_offset::AbstractVector # instance-wise "intercept" terms
    instance_offset_reg::AbstractMatrix # regularizer matrix (usually sparse)

    feature_offset::AbstractVector  # feature-wise "intercept" terms
    feature_offset_reg::AbstractMatrix # regularizer matrix (usually sparse)

    feature_precision::AbstractVector # feature-specific precisions

    instance_covariate_coeff::Union{Nothing,AbstractMatrix} # Regression coefficients
                                                          # for sample covariates
    instance_covariate_coeff_reg::AbstractMatrix

    losses::AbstractVector       # N-dim vector of feature-specific losses

    factor_importances::AbstractVector # K-dim vector of factor importances
end

function MatFacModel(instance_reg_mats::AbstractVector,
                     feature_reg_mats::AbstractVector,
                     losses::AbstractVector;
                     instance_offset_reg::Union{Nothing,AbstractMatrix}=nothing,
                     feature_offset_reg::Union{Nothing,AbstractMatrix}=nothing,
                     instance_covariate_coeff_reg::Union{Nothing,AbstractMatrix}=nothing,
                     K::Union{Nothing,Integer}=nothing)

    M = size(instance_reg_mats[1],1)
    N = size(feature_reg_mats[1],1)

    if K == nothing
        K = max(length(instance_reg_mats),
                length(feature_reg_mats))
    else
        @assert K >= max(length(instance_reg_mats),
                         length(feature_reg_mats))
    end

    X = 0.01 .* randn(K, M) ./ sqrt(K)
    Y = 0.01 .* randn(K, N)

    instance_offset = 0.01 .* randn(M)
    feature_offset = 0.01 .* randn(N)

    if instance_offset_reg == nothing
        instance_offset_reg = spzeros(M,M)
    end
    if feature_offset_reg == nothing
        feature_offset_reg = spzeros(N,N)
    end
    if instance_covariate_coeff_reg == nothing
        instance_covariate_coeff_reg = spzeros(N,N)
    end

    feature_precision = ones(N)

    instance_covariate_coeff = nothing

    factor_importances = zeros(K)

    return MatFacModel(X, Y, instance_reg_mats,
                             feature_reg_mats,
                             instance_offset,
                             instance_offset_reg,
                             feature_offset,
                             feature_offset_reg,
                             feature_precision,
                             instance_covariate_coeff,
                             instance_covariate_coeff_reg,
                             losses,
                             factor_importances
                      )
end


text:
using Reexport

@reexport module TriangulationEstimators

    import ..TransferOperator
    import ..TransferOperatorGenerator
    import ..TransferOperatorApproximation
    import ..transopergenerator
    import ..transferoperator
    import ..isboundarycondition
    import ..invariantmeasure
    import ..InvariantDistribution

    using Requires
    function __init__()
        @require Simplices="d5428e67-3037-59ba-9ab1-57a04f0a3b6a" begin
            export TriangulationBasedTransferOperator

            """ Supertype of all triangulation-based transfer operator estimators """
            abstract type TriangulationBasedTransferOperator <: TransferOperator end

            # Requirements and useful types
            include("simplex_types/simplex_types.jl")
            include("delaunay_triangulations/DelaunayTriangulations.jl")

            # Triangulation estimators
            include("common.jl")
            include("point/SimplexPoint.jl")
            include("exact/SimplexExact.jl")
            include("invariant_measure_simplex.jl")
        end
    end

end


text:
using Revise
using AdFem
using PyCall
using LinearAlgebra
using PyPlot
using SparseArrays
using MAT
using ADCMEKit
np = pyimport("numpy")

stepsize = 1
if length(ARGS)==1
  global stepsize = parse(Int64, ARGS[1])
end
@info stepsize

mode = "training"

## alpha-scheme
β = 1/4; γ = 1/2
a = b = 0.1    # damping

n = 10
m = 2n
h = 0.01
NT = 500
Δt = 2.0/NT
ηmax = 1
ηmin = 0.5

obs_idx = collect(1:stepsize:m+1)

bdedge = bcedge("right", m, n, h)
bdnode = bcnode("lower", m, n, h)

# λ = Variable(1.0)
# μ = Variable(1.0)
# invη = Variable(1.0)

function eta_model(idx)
  if idx == 1
    out = ηmin * ones(n)
    out[1:div(n,3)] .= ηmax
    out
  elseif idx==2
    out = ηmin * ones(n)
    out[1:div(n,3)] .= ηmax
    out[2div(n,3):end] .= ηmax
    out
  end

end

function visualize_inv_eta(X, k)
    x = LinRange(0.5h,m*h, m)
    y = LinRange(0.5h,n*h, n)
    V = zeros(m, n)
    for i = 1:m
        for j = 1:n
            elem = (j-1)*m + i
            V[i, j] = mean(X[4(elem-1)+1:4elem])
        end
    end
    close("all")
    pcolormesh(x, y, V'/50.0, vmin=ηmin-(ηmax-ηmin)/4, vmax=ηmax+(ηmax-ηmin)/4)
    colorbar(shrink=0.5)
    xlabel("x")
    ylabel("y")
    # title("Iteration = $k")
    axis("scaled")
    gca().invert_yaxis()
    if k == "true"
      title("True Model")
      savefig("true.png")
      return
    end
    k_ = string(k)
    k_ = reduce(*, "0" for i = 1:3-length(k_))*k_
    title("Iteration = $k_")
    savefig("iter$k_.png")
end

λ = constant(2.0)
μ = constant(0.2)
if mode=="data"
  global invη_var = constant(eta_model(2))
  invη = reshape(repeat(invη_var, 1, 4m), (-1,))
  global invη *= 50.0
else
    global invη_var = Variable((ηmin + ηmax)/2*ones(n))
    invη_ = reshape(repeat(invη_var, 1, 4m), (-1,))
    # invη_ = constant(eta_model(1))
    global invη = 50.0*invη_
end

fn_G = invη->begin
  G = tensor([1/Δt+2/3*μ*invη -μ/3*invη 0.0
    -μ/3*invη 1/Δt+2/3*μ*invη 0.0
    0.0 0.0 1/Δt+μ*invη])
  invG = inv(G)
end
invG = map(fn_G, invη)
S = tensor([2μ/Δt+λ/Δt λ/Δt 0.0
    λ/Δt 2μ/Δt+λ/Δt 0.0
    0.0 0.0 μ/Δt])
H = invG*S

M = compute_fem_mass_matrix1(m, n, h)
Zero = spzeros((m+1)*(n+1), (m+1)*(n+1))
M = SparseTensor([M Zero;Zero M])

K = compute_fem_stiffness_matrix(H, m, n, h)
C = a*M + b*K # damping matrix
L = M + γ*Δt*C + β*Δt^2*K
L, Lbd = fem_impose_Dirichlet_boundary_condition_experimental(L, bdnode, m, n, h)
L = factorize(L)

a = TensorArray(NT+1); a = write(a, 1, zeros(2(m+1)*(n+1))|>constant)
v = TensorArray(NT+1); v = write(v, 1, zeros(2(m+1)*(n+1))|>constant)
d = TensorArray(NT+1); d = write(d, 1, zeros(2(m+1)*(n+1))|>constant)
U = TensorArray(NT+1); U = write(U, 1, zeros(2(m+1)*(n+1))|>constant)
Sigma = TensorArray(NT+1); Sigma = write(Sigma, 1, zeros(4*m*n, 3)|>constant)
Varepsilon = TensorArray(NT+1); Varepsilon = write(Varepsilon, 1,zeros(4*m*n, 3)|>constant)

Forces = zeros(NT, 2(m+1)*(n+1))
for i = 1:NT
  T = eval_f_on_boundary_edge((x,y)->0.1, bdedge, m, n, h)

  # if i>=NT÷2
  #   T *= 0.0
  # end
  T = [-T T]
#   T = [T T]
  rhs = compute_fem_traction_term(T, bdedge, m, n, h)

#   if i*Δt>0.5
#     rhs = zero(rhs)
#   end
  Forces[i, :] = rhs
end
Forces = constant(Forces)

function condition(i, tas...)
  i <= NT
end

function body(i, tas...)
  a_, v_, d_, U_, Sigma_, Varepsilon_ = tas
  a = read(a_, i)
  v = read(v_, i)
  d = read(d_, i)
  U = read(U_, i)
  Sigma = read(Sigma_, i)
  Varepsilon = read(Varepsilon_, i)

  res = batch_matmul(invG/Δt, Sigma)
  F = compute_strain_energy_term(res, m, n, h) - K * U
  rhs = Forces[i] - F

  td = d + Δt*v + Δt^2/2*(1-2β)*a
  tv = v + (1-γ)*Δt*a
  rhs = rhs - C*tv - K*td
  rhs = scatter_update(rhs, constant([bdnode; bdnode.+(m+1)*(n+1)]), constant(zeros(2*length(bdnode))))

  ## alpha-scheme
  a = L\rhs # bottleneck
  d = td + β*Δt^2*a
  v = tv + γ*Δt*a
  U_new = d

  Varepsilon_new = eval_strain_on_gauss_pts(U_new, m, n, h)

  res2 = batch_matmul(invG * S, Varepsilon_new-Varepsilon)
  Sigma_new = res +  res2

  i+1, write(a_, i+1, a), write(v_, i+1, v), write(d_, i+1, d), write(U_, i+1, U_new),
        write(Sigma_, i+1, Sigma_new), write(Varepsilon_, i+1, Varepsilon_new)
end

i = constant(1, dtype=Int32)
_, _, _, _, u, sigma, varepsilon = while_loop(condition, body,
                  [i, a, v, d, U, Sigma, Varepsilon])

U = stack(u)
Sigma = stack(sigma)
Varepsilon = stack(varepsilon)

if mode!="data"
  data = matread("viscoelasticity.mat")
  global Uval,Sigmaval, Varepsilonval = data["U"], data["Sigma"], data["Varepsilon"]
  U.set_shape((NT+1, size(U, 2)))
  idx0 = 1:4m*n
  Sigma = map(x->x[idx0,:], Sigma)
  global loss = sum((U[:, obs_idx] - Uval[:, obs_idx])^2)
end

if !isdir(string(stepsize));mkdir(string(stepsize)); end
sess = Session(); init(sess)

cb = (v, i, l)->begin
  println("[$i] loss = $l")
  if i=="true" || mod(i,20)==0
    inv_eta = v[1]
    matwrite("$stepsize/eta$i.mat", Dict("eta"=>inv_eta))
  end
end

if mode=="data"
  Uval,Sigmaval, Varepsilonval = run(sess, [U, Sigma, Varepsilon])
  matwrite("viscoelasticity.mat", Dict("U"=>Uval, "Sigma"=>Sigmaval, "Varepsilon"=>Varepsilonval))

  # p = visualize_von_mises_stress(Sigmaval[1:5:end,:,:], m, n, h); saveanim(p, "space_s.gif")
  # p = visualize_displacement(Uval[1:5:end,:], m, n, h); saveanim(p, "space_u.gif")

  # visualize_von_mises_stress(Sigmaval[end,:,:], m, n, h); savefig("space_s.pdf")
  # visualize_displacement(Uval[end,:], m, n, h); savefig("space_u.pdf")

  visualize_inv_eta(run(sess, invη), "true")
  # cb([run(sess, invη)], "true", 0.0)
  error("Stop!")
end

v_ = []
i_ = []
l_ = []

@info run(sess, loss)
loss_ = BFGS!(sess, loss*1e10, vars=[invη],
    callback=cb, var_to_bounds=Dict(invη_var=>(0.1,2.0)))

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
using CheckedArithmeticCore
using BenchmarkTools
using Base.Checked # TODO: re-export

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1

xs = Dict{Type, Matrix}()
ys = Dict{Type, Matrix}()
zs = Dict{Type, Matrix}()

eltypes = (Int8, UInt8, Int16, UInt16, Int32, UInt32)
for T in eltypes
    push!(xs, T => rand(T, 1000, 1000))
    push!(ys, T => rand(T, 1000, 1000))
    push!(zs, T => zeros(T, 1000, 1000))
end

SUITE = BenchmarkGroup()
SUITE["add"] = BenchmarkGroup([],
    "wrapping" => BenchmarkGroup(),
    "saturating" => BenchmarkGroup(),
    "checked" => BenchmarkGroup(),
)

for T in eltypes
    x = xs[T]::Matrix{T}
    y = ys[T]::Matrix{T}
    z = zs[T]::Matrix{T}
    t = string(T)
    SUITE["add"]["checked"   ][t] = @benchmarkable    checked_add.($x, $z)
end


text:
"""
``RotationAngle3D`` is an abstract type to be used
as the first argument for functions.
It represents the rotation angle
in a 2D plane.
"""
abstract RotationAngle2D


text:
using BenchmarkTools
using Random

include("../examples/run.jl")

Random.seed!(7)
a = @benchmark run()
display(a)


text:
include("solenoid.jl")
include("clvs.jl")
include("lss.jl")
using Test
function test_dsolenoid()
	s = [1.,4]
	m = 10
	u_trj = rand(m,3)
	eps = 1.e-5
	du_x = (solenoid(u_trj .+ eps.*[1 0. 0.],
						s, 1)[2,:,:] -
				solenoid(u_trj .- eps.*[1 0. 0.],
						s, 1)[2,:,:])./(2*eps)
	du_y = (solenoid(u_trj .+ eps.*[0 1. 0.],
						s, 1)[2,:,:] -
				solenoid(u_trj .- eps.*[0 1. 0.],
						s, 1)[2,:,:])./(2*eps)
	du_z = (solenoid(u_trj .+ eps.*[0 0. 1.],
						s, 1)[2,:,:] -
				solenoid(u_trj .- eps.*[0 0. 1.],
						s, 1)[2,:,:])./(2*eps)
	du_fd = reshape(collect([du_x; du_y; du_z]),
						3, 3, m)
	du = dsolenoid(u_trj, s)
	@assert all(isapprox.(du_fd, du, rtol=1.e-2)) == 1
end
function test_les()
	s = [1., 4.]
	m = 1
	n = 25000
	u0 = rand(3,m)
	u_trj = lorenz63(u0, s, n)[:,:,1]
	du_trj = dlorenz63(u_trj, s)
	du_trj = permutedims(du_trj,[2,3,1])
	les, clv_trj = clvs(du_trj,3)
	println(les/0.005)
	@assert all(isapprox.(les/0.005, [0.91, 0.,
		-14.572], rtol=1.e-1)) == 1
end
function test_perturbation()
	s = [1.,4.]
	n = 10
	eps = 1.e-4
	u = rand(n,3)
	dFds_fd = (solenoid(u, s .+ eps.*[1,0], 1)[2,:,:] -
			   solenoid(u, s .- eps.*[1,0], 1)[2,:,:])/
			   (2*eps)
	dFds = perturbation(u, s)
	@test all(isapprox.(dFds_fd - dFds, 0., atol=1.e-8)) ==
	1
end
function test_lss()
	s = [3.,4.]
    m = 1
    n = 1000
    n_runup = 1000
    d = 3
    d_u = 1

	u0 = rand(m,d)
	u_init = transpose(solenoid(u0, s, n_runup
								)[end,:,:])
   	n_samples = 10
	dJds = zeros(n_samples)
	vsh = zeros(d,n+1)
	f = zeros(d,n+1)
	u = zeros(n+1,d)
	for i=1:n_samples
		println("Starting LSS, sample ", i)
		u_trj = solenoid(u_init, s, n)[:,:,1]
		du_trj = dsolenoid(u_trj, s)
		X = perturbation(u_trj,s) #ith col in T_{u_{i+1}} M
		J = objective(u_trj,s)
		dJ = dobjective(u_trj, s)
		v, dJds[i] = lss(u_trj, du_trj, X, f, J, dJ, s, d_u)
		vsh .= v
		u .= u_trj
		u_init .= reshape(u_trj[end,:],1,3)
		println(dJds[i])
	end
	# check with analytical shadowing direction
	# which is [r,0,0]
	n_spinup = 100
	u = u[n_spinup+1:end,:]
	vsh = vsh[:, n_spinup+1:end]
	n = n - n_spinup
	x, y, z = view(u,:,1), view(u,:,2),
				  view(u,:,3)
	r = sqrt.(x.*x .+ y.*y .+ z.*z)
	t_dir = [-y./r  x./r  zeros(n+1)]'
	r_dir = [x./r   y./r  zeros(n+1)]'
	prj_theta = [dot(t_dir[:,i], vsh[:,i]) for i=1:n]
	prj_z = [dot([0,0,1.], vsh[:,i]) for i=1:n]
	#@test all(isapprox.(prj_theta, 0., atol=0.01)) == 1
	#@test all(isapprox.(prj_z, 0., atol=0.01)) == 1
	@test isapprox((sum(dJds)/n_samples),1.0,rtol=0.1)
	return vsh, dJds
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
using TrafficAssignment
using Base.Test

ta_data = load_ta_network("Anaheim")
ta_data = load_ta_network("Barcelona")
ta_data = load_ta_network("Chicago Sketch")
ta_data = load_ta_network("Winnipeg")

data_time = time()
ta_data = load_ta_network("Sioux Falls")

println("Data Loading Completed, time:", time() - data_time, " seconds")

link_flow, link_travel_time, objective = ta_frank_wolfe(ta_data, method=:bfw, step=:exact, log=:on, tol=1e-3, max_iter_no=5)

link_flow, link_travel_time, objective = ta_frank_wolfe(ta_data, method=:cfw, step=:newton, log=:on, tol=1e-3, max_iter_no=5)

link_flow, link_travel_time, objective = ta_frank_wolfe(ta_data, method=:fw, step=:exact, log=:on, tol=1e-3, max_iter_no=5)
@assert abs( objective - 4.963799502172654e6 ) < 1e6


text:
# Default styles for Grammar of Graphics like interface
import ColorBrewer

const palette_dict = OrderedDict(
    :color => ColorBrewer.palette("Set1", 8),
    :markershape => [:diamond, :circle, :triangle, :star5],
    :linestyle => [:solid, :dash, :dot, :dashdot],
    :linewidth => [1,4,2,3],
    :markersize => [3,9,5,7]
)

const styles = collect(keys(palette_dict))

function get_style(t::SelectedData, by = splitby(t))
    grpd = groupby(first, t.table, by)
    args_dict = Dict{Symbol, Any}()
    for (key, val) in t.categoricalstyle
        s = unique(sort(columns(grpd, key)))
        d = Dict(zip(s, 1:length(s)))
        args_dict[val] = reshape(palette_dict[val][map(t->getindex(d, t), columns(grpd, key))], 1, :)
    end
    args_dict
end


text:
module GaussLegendre
    using FastGaussQuadrature

    function constract(integtablenum)
        return gausslegendre(integtablenum)
    end

    function gl_integ(f, x1, x2, vh_val)
        xm = 0.5 * (x1 + x2)
        xr = 0.5 * (x2 - x1)

        s = sum(i -> vh_val.w[i] * f(xm + xr * vh_val.x[i]), eachindex(vh_val.x))

        return s * xr
    end
end


text:
using Curiosity
using MinimalRLCore
using Statistics
using GVFHordes
using StatsBase
using JLD2
using Plots
const MCU = Curiosity.MountainCarUtils
const MCC = Curiosity.MountainCarConst

import ..get_normalized_state

include("../experiments/mountain_car.jl")

#Hacky this function should be called from environment
function get_normalized_state(env::MountainCar,state)
    pos_limit = MCC.pos_limit
    vel_limit = MCC.vel_limit
    return Float32[(state[1] - pos_limit[1])/(pos_limit[2] - pos_limit[1]),
                   (state[2] - vel_limit[1])/(vel_limit[2] - vel_limit[1])]
end

StatsBase.sample(p::GVFHordes.GVFParamFuncs.FunctionalPolicy, s, actions) =
    sample(Weights([p.func(;state_t = s, action_t = a) for a in actions]))

StatsBase.sample(rng, p::GVFHordes.GVFParamFuncs.AbstractPolicy, s, actions) =
    sample(Weights([get(p;state_t = s, action_t = a) for a in actions]))

function gen_start_state_eval_set()
    parsed = MountainCarExperiment.default_args()
    # horde = MountainCarExperiment.get_horde(parsed)
    # gvfs = [MCU.steps_to_wall_gvf(), MCU.steps_to_goal_gvf()]
    policies = MCU.get_policies(Dict("learned_policy_names" => ["Wall","Goal"] , "learned_policy" => true))
    gvfs = [MCU.steps_to_wall_gvf(policies[1]),MCU.steps_to_goal_gvf(policies[2])]

    normalized = true
    env = MountainCar(0.0,0.0, normalized)
    start_states = []

    pos_step_interval = 0.025
    vel_step_interval = 0.005
    # possible_pos = collect(0:0.05:1)
    # possible_vels = collect(0:0.1:1)
    possible_pos = collect(MCC.pos_initial_range[1]:pos_step_interval:MCC.pos_initial_range[2])
    possible_vels = collect(MCC.vel_initial_range[1]:vel_step_interval:MCC.vel_initial_range[2])

    states = []
    for p in possible_pos
        for v in possible_vels
            push!(states,get_normalized_state(env,[p,v]))
        end
    end
    # states = [[0.45,0.3]]
    actions = rand(1:3,length(states))
    gvf_rets = Array{Float64, 2}(undef, length(gvfs),length(states))

    num_returns = 1
    γ_thresh=1e-6

    for (gvf_i,gvf) in enumerate(gvfs)
        rets = monte_carlo_returns(env, gvf, states, actions,num_returns, γ_thresh)
        episode_lengths = [log.(0.99,r) for r in rets]
        println("Mean Episode length: ", mean(rets))
        rets = mean(rets, dims = 2)
        rets = collect(Iterators.flatten(rets))
        # scatter(x,rets, legend=false, ylabel="Cumulant Val", xlabel="Starting X Pos", title = "GVF: $( gvf_i)")
        # savefig("./MC_gvf_$(gvf_i).png")
        gvf_rets[gvf_i,:] = rets
    end
    return gvf_rets,states,actions
end

function gen_learned_eval_set()
    parsed = MountainCarExperiment.default_args()
    # horde = MountainCarExperiment.get_horde(parsed)
    # gvfs = [MCU.steps_to_wall_gvf(), MCU.steps_to_goal_gvf()]
    policies = MCU.get_policies(Dict("learned_policy_names" => ["Wall","Goal"] , "learned_policy" => true))
    gvfs = [MCU.steps_to_wall_gvf(policies[1]),MCU.steps_to_goal_gvf(policies[2])]

    normalized = true
    env = MountainCar(0.0,0.0, normalized)
    start_states = []

    pos_step_interval = 0.05
    vel_step_interval = 0.01
    possible_pos = collect(0:0.05:1)
    possible_vels = collect(0:0.1:1)
    # possible_pos = collect(MCC.pos_limit[1]:pos_step_interval:MCC.pos_limit[2])
    # possible_vels = collect(MCC.vel_limit[1]:vel_step_interval,MCC.vel_limit[2])

    states = []
    for p in possible_pos
        for v in possible_vels
            push!(states,[p,v])
        end
    end
    # states = [[0.45,0.3]]
    actions = rand(1:3,length(states))
    gvf_rets = Array{Float64, 2}(undef, length(gvfs),length(states))

    num_returns = 1
    γ_thresh=1e-6

    for (gvf_i,gvf) in enumerate(gvfs)
        println("here in gvf loop")
        rets = monte_carlo_returns(env, gvf, states, actions,num_returns, γ_thresh)
        @show size(rets)
        @show rets[1],rets[2],rets[3]
        episode_lengths = [log.(0.99,r) for r in rets]
        println("Mean Episode length: ", mean(rets))
        rets = mean(rets, dims = 2)
        rets = collect(Iterators.flatten(rets))
        # scatter(x,rets, legend=false, ylabel="Cumulant Val", xlabel="Starting X Pos", title = "GVF: $( gvf_i)")
        # savefig("./MC_gvf_$(gvf_i).png")
        gvf_rets[gvf_i,:] = rets
    end
    return gvf_rets,states,actions
end

function plot_eval_set(states,rets)
    x = [xy[1] for xy in states]
    y = [xy[2] for xy in states]

    s1 = scatter(x,y,rets[1,:], xlabel="Position", ylabel="Velocity", title = "To Back Wall")
    s2 = scatter(x,y,rets[2,:], xlabel="Position", ylabel="Velocity", title = "To Top Hill")
    plot([s1,s2]...)
end

function save_learned_policy(rets,states,actions)
    MCLearnedEvalSet = Dict()
    MCLearnedEvalSet["ests"] = rets
    MCLearnedEvalSet["actions"] = actions
    MCLearnedEvalSet["states"] = states
    @save "./src/data/MCLearnedEvalSet.jld2" MCLearnedEvalSet
end

function gen_eval_set()
    parsed = MountainCarExperiment.default_args()
    # horde = MountainCarExperiment.get_horde(parsed)
    gvfs = [MCU.steps_to_wall_gvf(), MCU.steps_to_goal_gvf()]

    normalized = true
    env = MountainCar(0.0,0.0, normalized)

    num_returns = 1
    γ_thresh=1e-6

    task_gvf = Curiosity.MountainCarUtils.task_gvf()
    # gvfs = [horde.gvfs..., task_gvf]
    # gvfs = [horde.gvfs..., task_gvf]

    num_start_states = 400
    gvf_rets = Array{Float64, 2}(undef, length(gvfs),num_start_states)
    start_states = []

    for i in 1:num_start_states
        MinimalRLCore.reset!(env)
        s = MinimalRLCore.get_state(env)
        push!(start_states,s)
    end
    possible_actions = MinimalRLCore.get_actions(env)
    action_probs = ones(length(possible_actions)) / length(possible_actions)
    start_actions = [sample(possible_actions, Weights(action_probs)) for i in 1:num_start_states]

    for (gvf_i,gvf) in enumerate(gvfs)
        rets = monte_carlo_returns(env, gvf, start_states, start_actions,num_returns, γ_thresh)

        x = [x for (x,y) in start_states]
        y = [y for (x,y) in start_states]

        rets = mean(rets, dims = 2)
        rets = collect(Iterators.flatten(rets))
        # scatter(x,rets, legend=false, ylabel="Cumulant Val", xlabel="Starting X Pos", title = "GVF: $( gvf_i)")
        # savefig("./MC_gvf_$(gvf_i).png")
        gvf_rets[gvf_i,:] = rets

    end
    return start_states, start_actions, gvf_rets
end

function save_data(rets,obs,actions)
    MCEvalSet = Dict()
    MCEvalSet["ests"] = rets
    MCEvalSet["actions"] = actions
    MCEvalSet["states"] = obs
    @save "./src/data/MCEvalSet.jld2" MCEvalSet
end

# using Plots
# pyplot()
# for (gvf_i,gvf) in enumerate(gvfs)
# # Create heatmap version
# increment = 0.01
# vel_limit = (0, 1)
# pos_initial_range = (0.4, 0.9)
#
# ys = collect(vel_limit[1]:increment:vel_limit[2])
# xs = collect(pos_initial_range[1]:increment:pos_initial_range[2])
#
# start_states = [[x,y] for x in xs for y in ys]
# rets = monte_carlo_returns(env, gvf, start_states, num_returns, γ_thresh)
# rets = collect(Iterators.flatten(rets))
# ret_heatmap = reshape(rets,length(xs),length(ys))
#
# surface(xs,ys,ret_heatmap', xlabel=" Starting X position", ylabel = "Starting Velocity", title = "GVF: $( gvf_i)")
# savefig("./MC_gvf_$(gvf_i).png")
#
# end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
""" Julia interface to the PyNeid API for accessing the NEID archive at NExScI """
__precompile__()
module NeidArchive
  using PyCall
  using Dates
  #using CSV, DataFrames

  const PyNeid = PyNULL()
  const archive = PyNULL()
  const valid_datalevel = [ "l0", "l1", "l2", "eng", "solarl0", "solarl1", "solarl2", "solareng"]
  const valid_format = [ "ipac", "votable", "csv", "tsv"]
  const default_format = "csv"
  const datefmt = DateFormat("yyyy-mm-dd HH:MM:SS")

  function __init__()
      copy!(PyNeid, pyimport("pyneid.neid"))
      copy!(archive, PyNeid.Neid)
  end

  #= Broken:  Code doesn't get an updated token.
  function login(;userid::String, password::String, token::String)
    neid.login(userid=userid, password=password, token=token)
    return token
  end
  =#

  """
      login(;userid, password, cookiepath)
  Login to the NEID archive at NExScI.  All three named arguments are required.  Login credentials stored in cookiepath.
  """
  function login(;userid::String="", password::String="", cookiepath::String)
    if length(userid)>=1 && length(password)>=1
       archive.login(userid=userid, password=password, cookiepath=cookiepath) #, debugfile="./archive.debug")
    else
       archive.login(cookiepath=cookiepath)
    end
 end

  """
      query_criteria(param::Dict{String,String}; cookiepath, [outpath] )
  Query the NEID archive using constraints in param.  Named argument cookiepath is required.
  """
  function query_criteria(param::Dict{String,String}; cookiepath::String, outpath::String=".", format::String = default_format)
    archive.query_criteria(param, cookiepath=cookiepath, format=format, outpath=outpath)
  end

  """
      query(param::Dict{String,String}; cookiepath, [outpath] )
  Shorthaned for query_criteria """
  query = query_criteria

# TODO: Implement query_adql, query_datetime, query_object, query_piname, query_program, query_qobject

  """
      download(filename, datalevel; cookiepath, [format, outdir, start_row, end_row] )
  Download files returned from a query of the NEID archive.
  Named argument cookiepath is required.
  """
  function download(filename::String, datalevel::String; format::String = default_format, outdir::String=".", cookiepath::String = "", start_row::Integer=0, end_row::Integer=1000)
    @assert datalevel ∈ valid_datalevel
    @assert format ∈ valid_format
    @assert 0 <= start_row <= end_row
    @assert end_row-start_row <= 10000

    if length(cookiepath) >= 1
      if start_row==0 && end_row==1000
         println("# Calling pyNEID archive.download without start_row or end_row")
         archive.download(filename, datalevel, format, outdir, cookiepath=cookiepath)
      else
         println("# Calling pyNEID archive.download with start_row=",start_row," end_row=",end_row)
         archive.download(filename, datalevel, format, outdir, cookiepath=cookiepath, start_row=start_row, end_row=end_row)
      end
    else
      if start_row==0 && end_row==1000
         println("# Calling pyNEID archive.download without start_row or end_row")
         archive.download(filename, datalevel, format, outdir)
      else
         println("# Calling pyNEID archive.download with start_row=",start_row," end_row=",end_row)
         archive.download(filename, datalevel, format, outdir, start_row=start_row, end_row=end_row)
      end
    end
  end

  #=
  # TODO: See if there's a clever way to pass all extra optional args to pyneid.  Maybe something like
  function download(filename::String, datalevel::String; format::String = default_format, outdir::String="."), kwargs...)
    archive.download(filename, datalevel, format, outdir, kwargs...)
  end
  =#

  """ Create a date range string in PyNeid format from start to stop."""
  function datetime_range(start::DateTime, stop::DateTime)
    Dates.format(start,datefmt) * "/" * Dates.format(stop,datefmt)
  end

  function datetime_range(start::Date, stop::Date)
    datetime_range(DateTime(start), DateTime(year(stop),month(stop),day(stop),23,59,59))
  end

  """ Create a date range string in PyNeid format for one day."""
  function datetime_one_day(day::Date)
    datetime_range(day, day)
  end

  """ Create a date range string in PyNeid format for 24 hours starting at given DateTime."""
  function datetime_one_day(start::DateTime)
    datetime_range(start, start+Dates.Day(1))
  end

  """ Create a date range string in PyNeid format for 24 hours starting at 15:00 on given Date."""
  function datetime_one_day_solar(start::Date)
    datetime_one_day(DateTime(start)+Dates.Hour(15))
  end

  """ Create a date range string in PyNeid format starting at 00:00:00 on start date."""
  function datetime_range_after(start::DateTime )
    Dates.format(start,datefmt) * "/"
  end

  function datetime_range_after(start::Date )
    datetime_range_after(DateTime(start))
  end

  """ Create a date range string in PyNeid format ending at 23:59:59 on stop date."""
  function datetime_range_before(stop::DateTime)
    datetime_range_before(DateTime(year(stop),month(stop),day(stop),23,59,59))
  end

end


text:
using Distributions
import Distributions: cdf, pdf

@doc raw"""
Weibull count model.

**Weibull count model** is a discrete distribution obtained as the *count model with inter-arrival times assumed to be i.i.d. Weibull distributed* [1, Section 3.1].

In words, a Weibull count model models the *number of events* which occur within some time `t`, where the time *between* each of the events are (i.i.d.) Weibull distributed.

The PDF is given by the following formula
```math
P \big( I(t) = i \big) = \sum_{j = 1}^{\infty} \frac{(-1)^{j + i} (\lambda t^c)^j a_j^i}{\Gamma(cj + 1)}, \quad i = 0, 1, 2, \dots
```
where
```math
\alpha_j^0 = \frac{\Gamma(cj + 1)}{\Gamma(j + 1)}, \quad j = 0, 1, 2, \dots
```
and
```math
\alpha_j^{i + 1} = \sum_{m = i}^{j - 1} \alpha_m^i \frac{\Gamma(cj - cm + 1)}{\Gamma(j - m + 1)}, \quad i = 0, 1, 2, \dots, \quad j = i + 1, i + 2, \dots
```
The derivation can be found in [1, Section 3.1].

## Notes
The computation of the PDF, CDF, etc. for a Weibull counting model requires approximation of an infinite series, which is highly unstable. We therefore make use of `BigInt` and `BigFloat` for all computations under the hood.
Most of the methods related to `WeibullCountModel` are dependent on the series approximation, and will therefore throw a `WeibullCountModel.DivergedError` upon failing to converge.

## References
[1] Adrian, M., Bradlow, E., Fader, P., & McShane, B., Count models based on weibull interarrival times, CoRR, (),  (2013).

"""
struct WeibullCountModel <: Distribution{Univariate, Discrete}
    λ::BigFloat
    c::BigFloat
    t::BigFloat
end

WeibullCountModel(λ::Number, c::Number, t::Number) = WeibullCountModel(BigFloat(λ), BigFloat(c), BigFloat(1.0))
WeibullCountModel(λ::Number, c::Number) = WeibullCountModel(λ, c, 1.0)

# PDF
pdf(d::WeibullCountModel, x::Array{BigInt}) = begin
    # extract the actual result of the computation
    weibull_count_pdf.(x, d.λ, d.c, d.t)[1]
end
pdf(d::WeibullCountModel, x::AbstractArray) = pdf.(d, BigInt.(x))

pdf(d::WeibullCountModel, x::BigInt) = weibull_count_pdf(x, d.λ, d.c)

# make explicit to prevent ambiguity from `Distributions.jl`
for T ∈ [Int64, Integer, Real, Number]
    eval(quote pdf(d::WeibullCountModel, x::$T) = pdf(d, BigInt(x)) end)
end

# CDF
cdf(d::WeibullCountModel, x::Array{BigInt}) = begin
    # extract the actual result of the computation
    weibull_count_cdf.(x, d.λ, d.c)[1]
end
cdf(d::WeibullCountModel, x::AbstractArray) = cdf(d, BigInt.(x))

cdf(d::WeibullCountModel, x::BigInt) = weibull_count_cdf(x, d.λ, d.c)
# make explicit to prevent ambiguity from `Distributions.jl`
for T ∈ [Int64, Integer, Real, Number]
    eval(quote cdf(d::WeibullCountModel, x::$T) = cdf(d, BigInt(x)) end)
end


text:
module Chipmunk

import Base: step, sleep
import SFML

dlsym = Base.Libdl.dlsym

function __init__()
    try
        deps = Pkg.dir("Chipmunk")*"/deps"
        global const libchipmunk = Libdl.dlopen("$deps/libchipmunk", Libdl.RTLD_GLOBAL)
        global const libchipmunkjl = Libdl.dlopen("$deps/libchipmunkjl", Libdl.RTLD_GLOBAL)
    catch exception
        println("Error: Could not find libchipmunk")
        println("Please rebuild")
    end
end

use_debug_draw = true
try
    global const debug_circle = SFML.CircleShape()
    global debug_line = SFML.Line(SFML.Vector2f(0, 0), SFML.Vector2f(0, 0), 0)
    global const debug_polygon = SFML.ConvexShape()
    global const debug_dot = SFML.CircleShape()
catch exception
    warn("SFML could not be loaded.")
    warn("Debug draw will be unavailable")
    use_debug_draw = false
end

include("julia/chipmunk.jl")
include("julia/chipmunk_types.jl")
include("julia/cpBody.jl")
include("julia/cpSpace.jl")
include("julia/cpShape.jl")
include("julia/cpPolyShape.jl")
include("julia/cpConstraint.jl")
include("julia/cpSimpleMotor.jl")
include("julia/cpPivotJoint.jl")
include("julia/cpPinJoint.jl")
include("julia/cpDebugDraw.jl")
include("julia/cpCollisionHandler.jl")

end


text:
####################################################
# Extract key/value combinations from a saved JLD dictionary
function extract(d)
    """
    Extracts variables in a dictionary into the workspace
    """
    expr = quote end
    for (k, v) in d
        push!(expr.args, :($(Symbol(k)) = $v))
    end
    eval(expr)
    return
end

####################################################
# Reduces the size of the variable dictionary to time t
function reduce_dict(exog_states, node)
    """
    Pulls out the current time exogenous variables
    """
    exog_states_t = Dict{Symbol,Float64}()

    for (key,value) in exog_states
        exog_states_t[key] = exog_states[key][node]
    end

    return exog_states_t

end

####################################################
# Gets bounds for simulated trajectories
function get_bounds(state, percentile)
"""
Recovers percentile bounds for simulated state paths
"""
    out = vcat([[quantile(hcat(state...)[row, :], percentile) quantile(hcat(state...)[row, :], 1 - percentile)] for row = 1:size(hcat(state...), 1)]...)
    return (out[:, 1], out[:, 2])
end

####################################################
# Extracts data from saved runs for use in welfare calculations
function extract_welfare_data(saved_data)
    """
    Pull out approximation space, coefficients, initial state, and parameters from
    a saved model
    """
    extract(saved_data)
    model = space
    coeff = c_store[1]
    if space[:d] < 4
        initial_state = [mparams[:kstart], mparams[:Mstart], mparams[:lossstart]]
    else
        initial_state = [mparams[:kstart], mparams[:Mstart], mparams[:expstart], mparams[:varstart], mparams[:lossstart]]
    end

    exog_states = simulate_exog_vars(mparams)

    return model, coeff, initial_state, mparams, exog_states
end

####################################################
# Compares ex ante welfare between two saved models
function ex_ante_welfare(model_1, model_2)
    """
    Compare welfare in terms of bge and pv lump sum
    """
    space_1, coeff_1, initial_state_1, mparams_1, exog_states = extract_welfare_data(model_1)
    space_2, coeff_2, initial_state_2, mparams_2, exog_states = extract_welfare_data(model_2)

    # Balanced growth equivalent
    bge = (-100*(1 .- (funeval(coeff_2, space_2, initial_state_2)[1].*mparams_2[:scaling_factor]/
        (funeval(coeff_1, space_1, initial_state_1)[1].*mparams_1[:scaling_factor]))^(1/(1-2))))[1]

    # Marginal utility function
    MU(x,t) = x.^(mparams[:rho] - 1)/mparams_1[:scaling_factor]

    # Lump sum benefit in 2005 (billion $)
    # -(V_1-V_2)/(MU)
    # A_0 L_0 adjusts MU from effective terms into real terms
    lump = (-(funeval(coeff_1, space_1, initial_state_1)[1].*mparams_1[:scaling_factor] .-
        funeval(coeff_2, space_2, initial_state_2)[1].*mparams_2[:scaling_factor])./
        MU(consumption[1], 1).*1e3.*(mparams_1[:A0].*mparams_1[:L0]))[1]
    return bge, lump
end

####################################################
# Computes ex-post welfare along a set of simulated trajectories
function ex_post_welfare(model, horizons, vars_in, mparams_in)
    """
    Computes ex-post sum of utility
    """
    # Extract consumption paths
    extract(model)
    cons = hcat(consumption...)
    cons = cons[:, 1:mparams_in[:num_runs]]

    # Present value utility fct in non-effective terms
    util_fct(x, t) = prod(vars_in[:betaeff][1:t]).*(x.^mparams_in[:rho])./mparams_in[:rho]/mparams[:scaling_factor]

    # Build utility array
    util = hcat([util_fct(cons[t, :], t) for t = 1:size(cons, 1)]...)'

    # Sum up over horizon length and take average utility across simulations
    avg_util = Vector{Float64}(undef, length(horizons))

    for (h_index, horizon) in enumerate(horizons)
        avg_util[h_index] = mean(sum(util[1:horizon, :], dims = 1))
    end

    return avg_util

end

max_stars_repo_name: irudik/optimal-climate-policy-aej

max_issues_repo_name: irudik/optimal-climate-policy-aej

max_forks_repo_name: irudik/optimal-climate-policy-aej


text:
@testset "BinomialLikelihood" begin
    lik = BinomialLikelihood(5)
    test_interface(lik, SqExponentialKernel(), rand(10))
end


text:
using IntervalArithmetic, IntervalRootFinding, SpecialFunctions
setformat(:full)

include("IntervalSpecialFunctions.jl")

# Complementary error function
function Phi_c(x::T) where {T}
    invsqrt2 = inv(sqrt(T(2)))
    erfc(x * invsqrt2)/2
end

# Standard normal p.d.f.
function phi(x::T) where {T}
    invsqrt2π = inv(sqrt(2*T(π)))
    exp(-abs2(x)/2) * invsqrt2π
end

# The Mills ratio
M(x) = Phi_c(x)/phi(x)

# The function to maximize
f(x) = x - x^2 * M(x)

# The first derivative of f
f_p(x) = 1 + x^2  - x*(2+x^2)*M(x)

# The second derivative of f
f_pp(x) = x^3 + 4*x - M(x)*(2+5*x^2+x^4)

rs = roots(f_p, f_pp, 0.0..3.0, Newton, 1e-15)
x0_interval = interval(rs[1])

println("We see there is a unique root:")
println(rs)
println("The root x0 occurs in the interval:")
println(x0_interval)
println("The value of f at x0 is in the interval:")
println(f(x0_interval))
println("The value of f' at x0 is in the interval:")
println(f_p(x0_interval))
println("The value of f'' at x0 is in the interval:")
println(f_pp(x0_interval))


text:
#include("../get_home_dir.jl")

using DataFrames
using CSV
using PyPlot
using Random
using StatsBase

#include(string(PROJECT_HOME,"/assim_code/time_averaging.jl"))

rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams");
rcParams["lines.linewidth"] = 1;
rcParams["font.size"] = 24;

mpa2mm = 10^6/9.8;
gravity_fac = (18.5+9)/2 * 1000/mpa2mm;

pmean = CSV.read("chain3_range/lwp_means.csv",DataFrame, header=0);
p25 = CSV.read("chain3_range/lwp_lower.csv",DataFrame,header=0) ;
p75 = CSV.read("chain3_range/lwp_upper.csv",DataFrame,header=0) ;

ptrue = CSV.read("chain3_range/lwp_true.csv",DataFrame,header=0) ;

colors_list = ["tab:blue","tab:green","tab:orange","tab:red"];
label_list = ["All obs.", "1 AM/PM obs.", "6 AM/PM obs."];

year_x = collect(1:(size(pmean)[1]))/365 .+ 2006;

figure()

plot(year_x, ptrue[:,1],color="k",linewidth=3, label="Truth")
for i in 1:3
	plot(year_x, pmean[:,i],color=colors_list[i],linewidth=3, label=label_list[i])
end

for i in 1:3
	fill_between(year_x, p25[:,i], p75[:,i] , color=colors_list[i],alpha=0.33)
end

legend()
xlabel("Time (years)")
ylabel("Leaf water potential (MPa)")
#xlim(2011, 2013.5)

#=

spmean = CSV.read("swp_means.csv",DataFrame, header=0);
sp25 = CSV.read("swp_lower.csv",DataFrame,header=0) ;
sp75 = CSV.read("swp_upper.csv",DataFrame,header=0) ;
sptrue = CSV.read("swp_true.csv",DataFrame,header=0) ;
#=
plot(year_x,ptrue[:,1],"k", label="True leaf")
#plot(year_x[1:(end-1)],(pmean[1:(end-1),1] + pmean[2:end,1])/2)
plot(year_x , sptrue[:,1] .- gravity_fac,"r:", label="True soil")

plot(year_x,pmean[:,1],"b", label="All leaf")
#plot(year_x[1:(end-1)],(pmean[1:(end-1),1] + pmean[2:end,1])/2)
plot(year_x , spmean[:,1] .- gravity_fac,":" , color="orange",label="All soil")
=#

figure()

plot(year_x, sptrue[:,1],color="k",linewidth=3, label="Truth")
for i in 1:3
	plot(year_x, spmean[:,i],color=colors_list[i],linewidth=3, label=label_list[i])
end

for i in 1:3
	fill_between(year_x, sp25[:,i], sp75[:,i] , color=colors_list[i],alpha=0.33)
end

legend()
xlabel("Time (years)")
ylabel("Soil water potential (MPa)")
#xlim(2011, 2013.5)

swmean = CSV.read("smc_means.csv",DataFrame, header=0);
sw25 = CSV.read("smc_lower.csv",DataFrame,header=0) ;
sw75 = CSV.read("smc_upper.csv",DataFrame,header=0) ;
swtrue = CSV.read("smc_true.csv",DataFrame,header=0) ;

figure()

plot(year_x, swtrue[:,1],color="k",linewidth=3, label="Truth")
for i in 1:3
	plot(year_x, swmean[:,i],color=colors_list[i],linewidth=3, label=label_list[i])
end

for i in 1:3
	fill_between(year_x, sw25[:,i], sw75[:,i] , color=colors_list[i],alpha=0.33)
end

legend()
xlabel("Time (years)")
ylabel("Soil moisture")
#xlim(2011, 2013.5)

=#

#need to start at 0.3 soil moisture in year 3


text:
"""Umm al-Qura calendar constants."""
module ummalqura

GREGORIAN_RANGE = ((1924, 8, 1), (2077, 11, 16))
"""Inclusive range of supported Gregorian dates (year, month and day)."""

HIJRI_RANGE = ((1343, 1, 1), (1500, 12, 30))
"""Inclusive range of supported Hijri dates (year, month and day)."""

HIJRI_OFFSET = 1342 * 12
"""Total Hijri months elapsed before the beginning of Hijri range."""

# insert types
MONTH_STARTS = (
    23999,
    24029,
    24058,
    24088,
    24118,
    24147,
    24177,
    24207,
    24237,
    24265,
    24295,
    24325,
    24355,
    24384,
    24413,
    24443,
    24472,
    24502,
    24531,
    24561,
    24590,
    24620,
    24649,
    24679,
    24708,
    24738,
    24767,
    24797,
    24826,
    24857,
    24886,
    24916,
    24944,
    24974,
    25004,
    25033,
    25063,
    25092,
    25121,
    25151,
    25181,
    25210,
    25240,
    25270,
    25299,
    25328,
    25358,
    25388,
    25417,
    25446,
    25475,
    25505,
    25535,
    25564,
    25594,
    25624,
    25653,
    25683,
    25713,
    25742,
    25771,
    25801,
    25830,
    25860,
    25889,
    25919,
    25948,
    25978,
    26008,
    26037,
    26066,
    26097,
    26125,
    26155,
    26184,
    26214,
    26243,
    26273,
    26302,
    26332,
    26362,
    26392,
    26420,
    26451,
    26481,
    26510,
    26540,
    26569,
    26599,
    26628,
    26657,
    26687,
    26716,
    26746,
    26776,
    26805,
    26835,
    26865,
    26894,
    26924,
    26953,
    26983,
    27012,
    27041,
    27070,
    27100,
    27130,
    27159,
    27189,
    27219,
    27248,
    27278,
    27307,
    27337,
    27366,
    27396,
    27425,
    27455,
    27484,
    27514,
    27543,
    27573,
    27602,
    27632,
    27662,
    27691,
    27721,
    27750,
    27780,
    27810,
    27839,
    27868,
    27898,
    27927,
    27957,
    27986,
    28016,
    28045,
    28075,
    28105,
    28134,
    28164,
    28193,
    28223,
    28252,
    28282,
    28311,
    28341,
    28370,
    28400,
    28429,
    28459,
    28488,
    28518,
    28548,
    28577,
    28607,
    28636,
    28665,
    28695,
    28724,
    28754,
    28783,
    28813,
    28843,
    28872,
    28901,
    28931,
    28960,
    28990,
    29019,
    29049,
    29078,
    29108,
    29137,
    29167,
    29196,
    29226,
    29255,
    29285,
    29315,
    29345,
    29375,
    29404,
    29434,
    29463,
    29492,
    29522,
    29551,
    29580,
    29610,
    29640,
    29669,
    29699,
    29729,
    29759,
    29788,
    29818,
    29847,
    29876,
    29906,
    29935,
    29964,
    29994,
    30023,
    30053,
    30082,
    30112,
    30141,
    30171,
    30200,
    30230,
    30259,
    30289,
    30318,
    30348,
    30378,
    30408,
    30437,
    30467,
    30496,
    30526,
    30555,
    30585,
    30614,
    30644,
    30673,
    30703,
    30732,
    30762,
    30791,
    30821,
    30850,
    30880,
    30909,
    30939,
    30968,
    30998,
    31027,
    31057,
    31086,
    31116,
    31145,
    31175,
    31204,
    31234,
    31263,
    31293,
    31322,
    31352,
    31381,
    31411,
    31441,
    31471,
    31500,
    31530,
    31559,
    31589,
    31618,
    31648,
    31676,
    31706,
    31736,
    31766,
    31795,
    31825,
    31854,
    31884,
    31913,
    31943,
    31972,
    32002,
    32031,
    32061,
    32090,
    32120,
    32150,
    32180,
    32209,
    32239,
    32268,
    32298,
    32327,
    32357,
    32386,
    32416,
    32445,
    32475,
    32504,
    32534,
    32563,
    32593,
    32622,
    32652,
    32681,
    32711,
    32740,
    32770,
    32799,
    32829,
    32858,
    32888,
    32917,
    32947,
    32976,
    33006,
    33035,
    33065,
    33094,
    33124,
    33153,
    33183,
    33213,
    33243,
    33272,
    33302,
    33331,
    33361,
    33390,
    33420,
    33450,
    33479,
    33509,
    33539,
    33568,
    33598,
    33627,
    33657,
    33686,
    33716,
    33745,
    33775,
    33804,
    33834,
    33863,
    33893,
    33922,
    33952,
    33981,
    34011,
    34040,
    34069,
    34099,
    34128,
    34158,
    34187,
    34217,
    34247,
    34277,
    34306,
    34336,
    34365,
    34395,
    34424,
    34454,
    34483,
    34512,
    34542,
    34571,
    34601,
    34631,
    34660,
    34690,
    34719,
    34749,
    34778,
    34808,
    34837,
    34867,
    34896,
    34926,
    34955,
    34985,
    35015,
    35044,
    35074,
    35103,
    35133,
    35162,
    35192,
    35222,
    35251,
    35280,
    35310,
    35340,
    35370,
    35399,
    35429,
    35458,
    35488,
    35517,
    35547,
    35576,
    35605,
    35635,
    35665,
    35694,
    35723,
    35753,
    35782,
    35811,
    35841,
    35871,
    35901,
    35930,
    35960,
    35989,
    36019,
    36048,
    36078,
    36107,
    36136,
    36166,
    36195,
    36225,
    36254,
    36284,
    36314,
    36343,
    36373,
    36403,
    36433,
    36462,
    36492,
    36521,
    36551,
    36580,
    36610,
    36639,
    36669,
    36698,
    36728,
    36757,
    36786,
    36816,
    36845,
    36875,
    36904,
    36934,
    36963,
    36993,
    37022,
    37052,
    37081,
    37111,
    37141,
    37170,
    37200,
    37229,
    37259,
    37288,
    37318,
    37347,
    37377,
    37406,
    37436,
    37465,
    37495,
    37524,
    37554,
    37584,
    37613,
    37643,
    37672,
    37701,
    37731,
    37760,
    37790,
    37819,
    37849,
    37878,
    37908,
    37938,
    37967,
    37997,
    38027,
    38056,
    38085,
    38115,
    38144,
    38174,
    38203,
    38233,
    38262,
    38292,
    38322,
    38351,
    38381,
    38410,
    38440,
    38469,
    38499,
    38528,
    38558,
    38587,
    38617,
    38646,
    38676,
    38705,
    38735,
    38764,
    38794,
    38823,
    38853,
    38882,
    38912,
    38941,
    38971,
    39001,
    39030,
    39059,
    39089,
    39118,
    39148,
    39178,
    39208,
    39237,
    39267,
    39297,
    39326,
    39355,
    39385,
    39414,
    39444,
    39473,
    39503,
    39532,
    39562,
    39592,
    39621,
    39650,
    39680,
    39709,
    39739,
    39768,
    39798,
    39827,
    39857,
    39886,
    39916,
    39946,
    39975,
    40005,
    40035,
    40064,
    40094,
    40123,
    40153,
    40182,
    40212,
    40241,
    40271,
    40300,
    40330,
    40359,
    40389,
    40418,
    40448,
    40477,
    40507,
    40536,
    40566,
    40595,
    40625,
    40655,
    40685,
    40714,
    40744,
    40773,
    40803,
    40832,
    40862,
    40892,
    40921,
    40951,
    40980,
    41009,
    41039,
    41068,
    41098,
    41127,
    41157,
    41186,
    41216,
    41245,
    41275,
    41304,
    41334,
    41364,
    41393,
    41422,
    41452,
    41481,
    41511,
    41540,
    41570,
    41599,
    41629,
    41658,
    41688,
    41718,
    41748,
    41777,
    41807,
    41836,
    41865,
    41894,
    41924,
    41953,
    41983,
    42012,
    42042,
    42072,
    42102,
    42131,
    42161,
    42190,
    42220,
    42249,
    42279,
    42308,
    42337,
    42367,
    42397,
    42426,
    42456,
    42485,
    42515,
    42545,
    42574,
    42604,
    42633,
    42662,
    42692,
    42721,
    42751,
    42780,
    42810,
    42839,
    42869,
    42899,
    42929,
    42958,
    42988,
    43017,
    43046,
    43076,
    43105,
    43135,
    43164,
    43194,
    43223,
    43253,
    43283,
    43312,
    43342,
    43371,
    43401,
    43430,
    43460,
    43489,
    43519,
    43548,
    43578,
    43607,
    43637,
    43666,
    43696,
    43726,
    43755,
    43785,
    43814,
    43844,
    43873,
    43903,
    43932,
    43962,
    43991,
    44021,
    44050,
    44080,
    44109,
    44139,
    44169,
    44198,
    44228,
    44258,
    44287,
    44317,
    44346,
    44375,
    44405,
    44434,
    44464,
    44493,
    44523,
    44553,
    44582,
    44612,
    44641,
    44671,
    44700,
    44730,
    44759,
    44788,
    44818,
    44847,
    44877,
    44906,
    44936,
    44966,
    44996,
    45025,
    45055,
    45084,
    45114,
    45143,
    45172,
    45202,
    45231,
    45261,
    45290,
    45320,
    45350,
    45380,
    45409,
    45439,
    45468,
    45498,
    45527,
    45556,
    45586,
    45615,
    45644,
    45674,
    45704,
    45733,
    45763,
    45793,
    45823,
    45852,
    45882,
    45911,
    45940,
    45970,
    45999,
    46028,
    46058,
    46088,
    46117,
    46147,
    46177,
    46206,
    46236,
    46265,
    46295,
    46324,
    46354,
    46383,
    46413,
    46442,
    46472,
    46501,
    46531,
    46560,
    46590,
    46620,
    46649,
    46679,
    46708,
    46738,
    46767,
    46797,
    46826,
    46856,
    46885,
    46915,
    46944,
    46974,
    47003,
    47033,
    47063,
    47092,
    47122,
    47151,
    47181,
    47210,
    47240,
    47269,
    47298,
    47328,
    47357,
    47387,
    47417,
    47446,
    47476,
    47506,
    47535,
    47565,
    47594,
    47624,
    47653,
    47682,
    47712,
    47741,
    47771,
    47800,
    47830,
    47860,
    47890,
    47919,
    47949,
    47978,
    48008,
    48037,
    48066,
    48096,
    48125,
    48155,
    48184,
    48214,
    48244,
    48273,
    48303,
    48333,
    48362,
    48392,
    48421,
    48450,
    48480,
    48509,
    48538,
    48568,
    48598,
    48627,
    48657,
    48687,
    48717,
    48746,
    48776,
    48805,
    48834,
    48864,
    48893,
    48922,
    48952,
    48982,
    49011,
    49041,
    49071,
    49100,
    49130,
    49160,
    49189,
    49218,
    49248,
    49277,
    49306,
    49336,
    49365,
    49395,
    49425,
    49455,
    49484,
    49514,
    49543,
    49573,
    49602,
    49632,
    49661,
    49690,
    49720,
    49749,
    49779,
    49809,
    49838,
    49868,
    49898,
    49927,
    49957,
    49986,
    50016,
    50045,
    50075,
    50104,
    50133,
    50163,
    50192,
    50222,
    50252,
    50281,
    50311,
    50340,
    50370,
    50400,
    50429,
    50459,
    50488,
    50518,
    50547,
    50576,
    50606,
    50635,
    50665,
    50694,
    50724,
    50754,
    50784,
    50813,
    50843,
    50872,
    50902,
    50931,
    50960,
    50990,
    51019,
    51049,
    51078,
    51108,
    51138,
    51167,
    51197,
    51227,
    51256,
    51286,
    51315,
    51345,
    51374,
    51403,
    51433,
    51462,
    51492,
    51522,
    51552,
    51582,
    51611,
    51641,
    51670,
    51699,
    51729,
    51758,
    51787,
    51816,
    51846,
    51876,
    51906,
    51936,
    51965,
    51995,
    52025,
    52054,
    52083,
    52113,
    52142,
    52171,
    52200,
    52230,
    52260,
    52290,
    52319,
    52349,
    52379,
    52408,
    52438,
    52467,
    52497,
    52526,
    52555,
    52585,
    52614,
    52644,
    52673,
    52703,
    52733,
    52762,
    52792,
    52822,
    52851,
    52881,
    52910,
    52939,
    52969,
    52998,
    53028,
    53057,
    53087,
    53116,
    53146,
    53176,
    53205,
    53235,
    53264,
    53294,
    53324,
    53353,
    53383,
    53412,
    53441,
    53471,
    53500,
    53530,
    53559,
    53589,
    53619,
    53648,
    53678,
    53708,
    53737,
    53767,
    53796,
    53825,
    53855,
    53884,
    53914,
    53943,
    53973,
    54003,
    54032,
    54062,
    54092,
    54121,
    54151,
    54180,
    54209,
    54239,
    54268,
    54297,
    54327,
    54357,
    54387,
    54416,
    54446,
    54476,
    54505,
    54535,
    54564,
    54593,
    54623,
    54652,
    54681,
    54711,
    54741,
    54770,
    54800,
    54830,
    54859,
    54889,
    54919,
    54948,
    54977,
    55007,
    55036,
    55066,
    55095,
    55125,
    55154,
    55184,
    55213,
    55243,
    55273,
    55302,
    55332,
    55361,
    55391,
    55420,
    55450,
    55479,
    55508,
    55538,
    55567,
    55597,
    55627,
    55657,
    55686,
    55716,
    55745,
    55775,
    55804,
    55834,
    55863,
    55892,
    55922,
    55951,
    55981,
    56011,
    56040,
    56070,
    56100,
    56129,
    56159,
    56188,
    56218,
    56247,
    56276,
    56306,
    56335,
    56365,
    56394,
    56424,
    56454,
    56483,
    56513,
    56543,
    56572,
    56601,
    56631,
    56660,
    56690,
    56719,
    56749,
    56778,
    56808,
    56837,
    56867,
    56897,
    56926,
    56956,
    56985,
    57015,
    57044,
    57074,
    57103,
    57133,
    57162,
    57192,
    57221,
    57251,
    57280,
    57310,
    57340,
    57369,
    57399,
    57429,
    57458,
    57487,
    57517,
    57546,
    57576,
    57605,
    57634,
    57664,
    57694,
    57723,
    57753,
    57783,
    57813,
    57842,
    57871,
    57901,
    57930,
    57959,
    57989,
    58018,
    58048,
    58077,
    58107,
    58137,
    58167,
    58196,
    58226,
    58255,
    58285,
    58314,
    58343,
    58373,
    58402,
    58432,
    58461,
    58491,
    58521,
    58551,
    58580,
    58610,
    58639,
    58669,
    58698,
    58727,
    58757,
    58786,
    58816,
    58845,
    58875,
    58905,
    58934,
    58964,
    58994,
    59023,
    59053,
    59082,
    59111,
    59141,
    59170,
    59200,
    59229,
    59259,
    59288,
    59318,
    59348,
    59377,
    59407,
    59436,
    59466,
    59495,
    59525,
    59554,
    59584,
    59613,
    59643,
    59672,
    59702,
    59731,
    59761,
    59791,
    59820,
    59850,
    59879,
    59909,
    59939,
    59968,
    59997,
    60027,
    60056,
    60086,
    60115,
    60145,
    60174,
    60204,
    60234,
    60264,
    60293,
    60323,
    60352,
    60381,
    60411,
    60440,
    60469,
    60499,
    60528,
    60558,
    60588,
    60618,
    60647,
    60677,
    60707,
    60736,
    60765,
    60795,
    60824,
    60853,
    60883,
    60912,
    60942,
    60972,
    61002,
    61031,
    61061,
    61090,
    61120,
    61149,
    61179,
    61208,
    61237,
    61267,
    61296,
    61326,
    61356,
    61385,
    61415,
    61445,
    61474,
    61504,
    61533,
    61563,
    61592,
    61621,
    61651,
    61680,
    61710,
    61739,
    61769,
    61799,
    61828,
    61858,
    61888,
    61917,
    61947,
    61976,
    62006,
    62035,
    62064,
    62094,
    62123,
    62153,
    62182,
    62212,
    62242,
    62271,
    62301,
    62331,
    62360,
    62390,
    62419,
    62448,
    62478,
    62507,
    62537,
    62566,
    62596,
    62625,
    62655,
    62685,
    62715,
    62744,
    62774,
    62803,
    62832,
    62862,
    62891,
    62921,
    62950,
    62980,
    63009,
    63039,
    63069,
    63099,
    63128,
    63157,
    63187,
    63216,
    63246,
    63275,
    63305,
    63334,
    63363,
    63393,
    63423,
    63453,
    63482,
    63512,
    63541,
    63571,
    63600,
    63630,
    63659,
    63689,
    63718,
    63747,
    63777,
    63807,
    63836,
    63866,
    63895,
    63925,
    63955,
    63984,
    64014,
    64043,
    64073,
    64102,
    64131,
    64161,
    64190,
    64220,
    64249,
    64279,
    64309,
    64339,
    64368,
    64398,
    64427,
    64457,
    64486,
    64515,
    64545,
    64574,
    64603,
    64633,
    64663,
    64692,
    64722,
    64752,
    64782,
    64811,
    64841,
    64870,
    64899,
    64929,
    64958,
    64987,
    65017,
    65047,
    65076,
    65106,
    65136,
    65166,
    65195,
    65225,
    65254,
    65283,
    65313,
    65342,
    65371,
    65401,
    65431,
    65460,
    65490,
    65520,
    65549,
    65579,
    65608,
    65638,
    65667,
    65697,
    65726,
    65755,
    65785,
    65815,
    65844,
    65874,
    65903,
    65933,
    65963,
    65992,
    66022,
    66051,
    66081,
    66110,
    66140,
    66169,
    66199,
    66228,
    66258,
    66287,
    66317,
    66346,
    66376,
    66405,
    66435,
    66465,
    66494,
    66524,
    66553,
    66583,
    66612,
    66641,
    66671,
    66700,
    66730,
    66760,
    66789,
    66819,
    66849,
    66878,
    66908,
    66937,
    66967,
    66996,
    67025,
    67055,
    67084,
    67114,
    67143,
    67173,
    67203,
    67233,
    67262,
    67292,
    67321,
    67351,
    67380,
    67409,
    67439,
    67468,
    67497,
    67527,
    67557,
    67587,
    67617,
    67646,
    67676,
    67705,
    67735,
    67764,
    67793,
    67823,
    67852,
    67882,
    67911,
    67941,
    67971,
    68000,
    68030,
    68060,
    68089,
    68119,
    68148,
    68177,
    68207,
    68236,
    68266,
    68295,
    68325,
    68354,
    68384,
    68414,
    68443,
    68473,
    68502,
    68532,
    68561,
    68591,
    68620,
    68650,
    68679,
    68708,
    68738,
    68768,
    68797,
    68827,
    68857,
    68886,
    68916,
    68946,
    68975,
    69004,
    69034,
    69063,
    69092,
    69122,
    69152,
    69181,
    69211,
    69240,
    69270,
    69300,
    69330,
    69359,
    69388,
    69418,
    69447,
    69476,
    69506,
    69535,
    69565,
    69595,
    69624,
    69654,
    69684,
    69713,
    69743,
    69772,
    69802,
    69831,
    69861,
    69890,
    69919,
    69949,
    69978,
    70008,
    70038,
    70067,
    70097,
    70126,
    70156,
    70186,
    70215,
    70245,
    70274,
    70303,
    70333,
    70362,
    70392,
    70421,
    70451,
    70481,
    70510,
    70540,
    70570,
    70599,
    70629,
    70658,
    70687,
    70717,
    70746,
    70776,
    70805,
    70835,
    70864,
    70894,
    70924,
    70954,
    70983,
    71013,
    71042,
    71071,
    71101,
    71130,
    71159,
    71189,
    71218,
    71248,
    71278,
    71308,
    71337,
    71367,
    71397,
    71426,
    71455,
    71485,
    71514,
    71543,
    71573,
    71602,
    71632,
    71662,
    71691,
    71721,
    71751,
    71781,
    71810,
    71839,
    71869,
    71898,
    71927,
    71957,
    71986,
    72016,
    72046,
    72075,
    72105,
    72135,
    72164,
    72194,
    72223,
    72253,
    72282,
    72311,
    72341,
    72370,
    72400,
    72429,
    72459,
    72489,
    72518,
    72548,
    72577,
    72607,
    72637,
    72666,
    72695,
    72725,
    72754,
    72784,
    72813,
    72843,
    72872,
    72902,
    72931,
    72961,
    72991,
    73020,
    73050,
    73080,
    73109,
    73139,
    73168,
    73197,
    73227,
    73256,
    73286,
    73315,
    73345,
    73375,
    73404,
    73434,
    73464,
    73493,
    73523,
    73552,
    73581,
    73611,
    73640,
    73669,
    73699,
    73729,
    73758,
    73788,
    73818,
    73848,
    73877,
    73907,
    73936,
    73965,
    73995,
    74024,
    74053,
    74083,
    74113,
    74142,
    74172,
    74202,
    74231,
    74261,
    74291,
    74320,
    74349,
    74379,
    74408,
    74437,
    74467,
    74497,
    74526,
    74556,
    74585,
    74615,
    74645,
    74675,
    74704,
    74733,
    74763,
    74792,
    74822,
    74851,
    74881,
    74910,
    74940,
    74969,
    74999,
    75029,
    75058,
    75088,
    75117,
    75147,
    75176,
    75206,
    75235,
    75264,
    75294,
    75323,
    75353,
    75383,
    75412,
    75442,
    75472,
    75501,
    75531,
    75560,
    75590,
    75619,
    75648,
    75678,
    75707,
    75737,
    75766,
    75796,
    75826,
    75856,
    75885,
    75915,
    75944,
    75974,
    76003,
    76032,
    76062,
    76091,
    76121,
    76150,
    76180,
    76210,
    76239,
    76269,
    76299,
    76328,
    76358,
    76387,
    76416,
    76446,
    76475,
    76505,
    76534,
    76564,
    76593,
    76623,
    76653,
    76682,
    76712,
    76741,
    76771,
    76801,
    76830,
    76859,
    76889,
    76918,
    76948,
    76977,
    77007,
    77036,
    77066,
    77096,
    77125,
    77155,
    77185,
    77214,
    77243,
    77273,
    77302,
    77332,
    77361,
    77390,
    77420,
    77450,
    77479,
    77509,
    77539,
    77569,
    77598,
    77627,
    77657,
    77686,
    77715,
    77745,
    77774,
    77804,
    77833,
    77863,
    77893,
    77923,
    77952,
    77982,
    78011,
    78041,
    78070,
    78099,
    78129,
    78158,
    78188,
    78217,
    78247,
    78277,
    78307,
    78336,
    78366,
    78395,
    78425,
    78454,
    78483,
    78513,
    78542,
    78572,
    78601,
    78631,
    78661,
    78690,
    78720,
    78750,
    78779,
    78808,
    78838,
    78867,
    78897,
    78926,
    78956,
    78985,
    79015,
    79044,
    79074,
    79104,
    79133,
    79163,
    79192,
    79222,
    79251,
    79281,
    79310,
    79340,
    79369,
    79399,
    79428,
    79458,
    79487,
    79517,
    79546,
    79576,
    79606,
    79635,
    79665,
    79695,
    79724,
    79753,
    79783,
    79812,
    79841,
    79871,
    79900,
    79930,
    79960,
    79990,
)
"""Ordered list of Reduced Julian Day (RJD) numbers for the beginning of supported
Hijri months.
"""
end

max_stars_repo_name: ceferisbarov/hijri-converter

max_issues_repo_name: ceferisbarov/hijri-converter


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule HOHQMesh_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("HOHQMesh")
JLLWrappers.@generate_main_file("HOHQMesh", UUID("1d5cbd98-5122-5a8a-bea1-c186d986ee7f"))
end  # module HOHQMesh_jll


text:
using Test

# https://discourse.julialang.org/t/getting-data-directly-from-a-website/12241/6

fNames = ["train-images-idx3-ubyte","train-labels-idx1-ubyte","t10k-images-idx3-ubyte","t10k-labels-idx1"]
origPath = ["http://yann.lecun.com/exdb/mnist/"]
destPath = joinpath(dirname(Base.find_package("BetaML")),"..","test","data","mnist")

iris     = readdlm(joinpath(dirname(Base.find_package("BetaML")),"..","test","data","iris.csv"),',',skipstart=1)

http://yann.lecun.com/exdb/mnist/train-images-idx3-ubyte.gz
http://yann.lecun.com/exdb/mnist/train-labels-idx1-ubyte.gz
http://yann.lecun.com/exdb/mnist/t10k-images-idx3-ubyte.gz
http://yann.lecun.com/exdb/mnist/t10k-labels-idx1-ubyte.gz

r = HTTP.get(origPath*fNames[1]*".gz", cookies=true);

using HTTP, GZip, IDX # https://github.com/jlegare/IDX.git
r = HTTP.get("http://yann.lecun.com/exdb/mnist/train-images-idx3-ubyte.gz", cookies=true);
destPath = joinpath(dirname(Base.find_package("BetaML")),"..","test","data","minst")
zippedFile = joinpath(destPath,"test.gz")
unZippedFile = joinpath(destPath,"test.idx3")
open(zippedFile,"w") do f
    write(f,String(r.body))
end
fh = GZip.open(zippedFile)
open(unZippedFile,"w") do f
    write(f,read(fh))
end
train_set = load(unZippedFile)
img1 = train_set[3][:,:,1]

train_set = load(read(fh))

a = "aaa"
b = "bb"
c = a*b


text:
module AccurateArithmetic

export
       four, two, negone, onehalf, onequarter,
       fourx, twox, halfx, quarterx,
#=
       These functions appear in the literature, named similarly.
=#
       two_sum, two_diff, two_prod,
       fast_two_sum, fast_two_diff,
       ufp, ulp, splitting, extractscalar,
#=
     These functions evaluate as multipart floats.
     The trailing '_' signifies this extra accuracy.
=#
       add_, sub_, add_hilo_, sub_hilo_,
       mul_, inv_, dvi_,
       powr2_, root2_, # powr3_, root3_,
       sum_,
       add_2, add_3, add_4, add_5,
       add_hilo_2, add_hilo_3, add_hilo_4, add_hilo_5,
       sub_2, sub_3,
       sub_hilo_2, sub_hilo_3,
       powr2_2, #powr3_2, powr3_3,
       root2_2, #root3_2,
       mul_2, mul_3,
       inv_2, dvi_2,
       fma_, fma_3, fma_2, fms_, fms_3, fms_2

if VERSION >= v"0.7.0-"
    import Base.IEEEFloat
else
    const IEEEFloat = Union{Float64, Float32, Float16}
end

import Base: signbit, flipsign, copysign, abs, inv, (+), (-), (*), (/)

include("constants.jl")
include("twoarith.jl")
include("renormalize.jl")
include("errorfree.jl")
include("faithful.jl")
include("nearfaithful.jl")
include("compensated.jl")
include("pairarith.jl")

end # module AccurateAritmetic