text:
#---------------------------------------------------------
# # [River discharge](@id river-discharge)
#---------------------------------------------------------

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/5_river_discharge.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/5_river_discharge.ipynb)

#md # !!! note
#md #     All the AIBECS tutorials and how-to guides are available as Jupyter notebooks.
#md #     You can view them with [nbviewer](https://nbviewer.jupyter.org/)or execute them online with [binder](https://mybinder.org/) by clicking on the badges above!
#md #     (Note that binder can be slow to launch and its memory caps can be a problem when running.)

# In this tutorial we will simulate a fictitious radioactive tracer that is injected into the ocean by the 200 largest rivers (by estimated discharge).
# The 200 major rivers dataset from [*Dai and Trenberth* (2002)](https://rda.ucar.edu/datasets/ds551.0/) is available from within the AIBECS.
# Once "born", our ficitious tracer decays with a parameter timescale $\tau$ as it flows through ocean basins.
#
# The 3D tracer equation is:
#
# $$\left[\frac{\partial}{\partial t} + \nabla \cdot (\boldsymbol{u} + \mathbf{K}\nabla)\right] x = s_\mathsf{rivers} - x / \tau$$
#
# where $\nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \nabla \right]$ represents the ocean circulation transport.
# ([Tracer transport operators are described in the documentation](@ref tracer-transport-operators).)
# The riverine source of the tracer is $s_\mathsf{rivers}$, and $x / \tau$ is the decay rate.

# In AIBECS, we must recast this equation in the generic form
#
# $$\left[\frac{\partial}{\partial t} + \mathbf{T}(\boldsymbol{p})\right] \boldsymbol{x} = \boldsymbol{G}(\boldsymbol{x},\boldsymbol{p}).$$

# We start by telling Julia that we want to use the AIBECS and the OCIM2 transport matrix for the ocean circulation.

using AIBECS
grd, T_OCIM2 = OCIM2.load()

# For the radioactive decay, we simply use

function decay(x, p)
    @unpack τ = p
    return x / τ
end

# To build the river sources, we will load the geographic locations and discharge (in m³ s⁻¹) from the [*Dai and Trenberth* (2017) dataset](https://rda.ucar.edu/datasets/ds551.0/index.html#!description).

RIVERS = Rivers.load()

# This is an array of rivers, for which the type `River{T}` contains the river's name, lat–lon coordinates, and discharge in m³ s⁻¹.
# For example, the first river is the Amazon

r = RIVERS[1]

# We can check the locations with

using Plots
scatter([r.lon for r in RIVERS], [r.lat for r in RIVERS],
        zcolor=log10.(uconvert.(NoUnits, [r.VFR for r in RIVERS] / u"m^3/s")),
        colorbartitle="log₁₀(discharge / (1 m³ s⁻¹))")

# We can regrid these into the OCIM2 grid and return the corresponding vector with

rivers = regrid(RIVERS, grd)

# (Note this regridding uses [NearestNeighbors.jl](https://github.com/KristofferC/NearestNeighbors.jl) to assign a wet box as the mouth of each river, which sometimes is not exactly the real location of the river mouth.)

# We control the global magnitude of the river discharge, $\sigma$ (in mol s⁻¹), by making it a parameter of our model.
# For that, we separate the river source
#
# $$s_\mathsf{rivers} = \sigma s_0$$
#
# into global magnitude ($\sigma$) and spatial pattern ($s_0$).
#
# Since $\int s_0 \mathrm{d}V = 1$, `s_0` can be computed by normalizing `rivers`.
# In Julia/AIBECS, this can be done by dividing `rivers` by the dot product `v ⋅ rivers` (or `v'rivers` in matrix form).
# (`v ⋅ x` is the discrete equivalent of the volume integral $\int x \mathrm{d}V$.)

v = vector_of_volumes(grd)
s_0 = rivers / (v'rivers)
function s_rivers(p)
    @unpack σ = p
    return σ * ustrip.(s_0) # we must remove the units in AIBECS here :(
end

# We then write the generic $\boldsymbol{G}$ function, which is

G_radiorivers(x,p) = s_rivers(p) - decay(x,p)

# ##### Parameters

# We specify some initial values for the parameters and also include units.

import AIBECS: @units, units
import AIBECS: @initial_value, initial_value
@initial_value @units struct RadioRiversParameters{U} <: AbstractParameters{U}
    τ::U | 5.0 | u"yr"
    σ::U | 1.0 | u"Gmol/yr"
end

# Finally, thanks to the initial values we provided, we can instantiate the parameter vector succinctly as

p = RadioRiversParameters()

# We generate the state function `F` that gives the tendencies,

F = AIBECSFunction(T_OCIM2, G_radiorivers)

# generate the steady-state problem `prob`,

nb = sum(iswet(grd))
x = ones(nb) # initial guess
prob = SteadyStateProblem(F, x, p)

# and solve it

sol = solve(prob, CTKAlg()).u * u"mol/m^3"

# Let's now run some visualizations using the plot recipes.
# Taking a horizontal slice of the 3D field at 200m gives

cmap = :viridis
plothorizontalslice(sol, grd, zunit=u"μmol/m^3", depth=200, color=cmap)

# and at 500m:

plothorizontalslice(sol, grd, zunit=u"μmol/m^3", depth=500, color=cmap, clim=(0,0.05))

# Or we can change the timescale and watch the tracer fill the oceans:

p = RadioRiversParameters(τ = 50.0u"yr")
prob = SteadyStateProblem(F, x, p)
sol_τ50 = solve(prob, CTKAlg()).u * u"mol/m^3"
plothorizontalslice(sol_τ50, grd, zunit=u"μmol/m^3", depth=500, color=cmap, clim=(0,1))

# Point-wise sources like the rivers used here can be problematic numerically because these
# can generate large gradients and numerical noise with the given spatial discretization:

cmap = :RdYlBu_4
plothorizontalslice(sol_τ50, grd, zunit=u"μmol/m^3", depth=0, color=cmap, clim=(-10,10))

# Note, for an example of numerical noise, the negative values appearing in red.
# Hence, it might be wise to smooth the 3D field of the river sources.
# We can do this using the `smooth_operator` function
# (`smooth_operator` creates an operator that can be applied to smooth the AIBECS vector
# along the stencil of the given transport matrix and in a volume-conservative way).
# And we can check that the negative values disappear after smoothing of the sources:

S = smooth_operator(grd, T_OCIM2)
s_0 = S * (S * s_0) # smooth river sources twice
sol_smooth = solve(prob, CTKAlg()).u * u"mol/m^3"
plothorizontalslice(sol_smooth, grd, zunit=u"μmol/m^3", depth=0, color=cmap, clim=(-10,10))

# Note that such numerical noise generally dampens out with distance and depth

cmap = :viridis
plot(
    plothorizontalslice(sol_τ50, grd, zunit=u"μmol/m^3", depth=200, color=cmap, clim=(0,2)),
    plothorizontalslice(sol_smooth, grd, zunit=u"μmol/m^3", depth=200, color=cmap, clim=(0,2)),
    layout=(2,1),
    size=(600,600)
)

# Nevertheless, it can be good to reduce noise as much as possible!


text:
module GOESRLevel0

using Reexport

include("Rice/Rice.jl")
@reexport using .Rice

"""
Descibes the type of the scene

| name           | value  | Description                              |
|:-------------- |:------ |:---------------------------------------- |
| `disk`         | `0x00` | Whole visible disk                       |
| `conus`        | `0x01` | CONUS region                             |
| `meso1`        | `0x02` | Meso 1 region                            |
| `meso2`        | `0x03` | Meso 2 region                            |
| `starlook_vis` | `0x05` | VIS starlook calibration                 |
| `starlook_ir`  | `0x0a` | IR starlook calibration                  |
| `cal_ir`       | `0x0b` | IR calibration (black body calibration?) |
| `spacelook`    | `0x0d` | space look calibration                   |

"""
@enum SceneType::UInt8 begin
    disk = 0x00
    conus = 0x01
    meso1 = 0x02
    meso2 = 0x03
    starlook_vis = 0x05
    starlook_ir = 0x0a
    cal_ir = 0x0b
    spacelook = 0x0d
end

export SceneType

end # module

ext: jl


text:
using Test

using Dates
using DelimitedFiles
using LinearAlgebra
using Printf
using ReferenceFrameRotations
using SatelliteToolbox
using SatelliteToolbox.SGP4
using StaticArrays

@testset "Atmospheric Models" begin
    cd("./earth/atmospheric_models/expatmosphere/")
    include("./earth/atmospheric_models/expatmosphere/expatmosphere.jl")
    cd("../../../")
    cd("./earth/atmospheric_models/nrlmsise00/")
    include("./earth/atmospheric_models/nrlmsise00/nrlmsise00.jl")
    cd("../../../")
    cd("./earth/atmospheric_models/jb2008/")
    include("./earth/atmospheric_models/jb2008/jb2008.jl")
    cd("../../../")
    cd("./earth/atmospheric_models/jr1971/")
    include("./earth/atmospheric_models/jr1971/jr1971.jl")
    cd("../../../")
end
println("")

@testset "Geomagnetic Field Models" begin
    cd("./earth/geomagnetic_field_models/")
    include("./earth/geomagnetic_field_models/dipole.jl")
    include("./earth/geomagnetic_field_models/igrf.jl")
    cd("../../")
end
println("")

@testset "Gravity Field Models" begin
    cd("./earth/gravity_models/")
    include("./earth/gravity_models/gravity_models.jl")
    cd("../../")
end
println("")

@testset "Space Indices" begin
    cd("./earth/space_indices/")
    include("./earth/space_indices/space_indices.jl")
    cd("../../")
end
println("")

@testset "General orbit functions" begin
    include("./orbit/general.jl")
end
println("")

@testset "Orbit propagators" begin
    cd("./orbit/propagators")
    include("./orbit/propagators/sgp4.jl")
    include("./orbit/propagators/twobody.jl")
    cd("../../")
end
println("")

@testset "Coordinate transformations" begin
    cd("./transformations/fk5/")
    include("./transformations/fk5/fk5.jl")
    cd("../../")
    cd("./transformations/iau2006/")
    include("./transformations/iau2006/iau2006.jl")
    cd("../../")
    cd("./transformations/teme/")
    include("./transformations/teme/teme.jl")
    cd("../../")
    cd("./transformations/")
    include("./transformations/ecef_to_ecef.jl")
    include("./transformations/ecef_to_eci.jl")
    include("./transformations/eci_to_ecef.jl")
    include("./transformations/eci_to_eci.jl")
    include("./transformations/geodetic_geocentric.jl")
    include("./transformations/orbit_elements.jl")
    include("./transformations/sv_ecef_to_ecef.jl")
    include("./transformations/sv_ecef_to_eci.jl")
    include("./transformations/sv_eci_to_ecef.jl")
    include("./transformations/sv_eci_to_eci.jl")
    cd("../")
end
println("")

@testset "Functions related with time" begin
     include("./time/time.jl")
end
println("")

@testset "Functions related with the Sun" begin
    include("./sun/equation_of_time.jl")
    include("./sun/sun_position.jl")
    include("./sun/sun_velocity.jl")
end
println("")

@testset "Functions related with the Moon" begin
    include("./moon/moon_position.jl")
end
println("")

@testset "Submodules" begin
    cd("./submodules/")
    include("./submodules/submodules.jl")
    cd("../")
end


text:
# ---
# title: 1105. Filling Bookcase Shelves
# id: problem1105
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/filling-bookcase-shelves/description/>
# hidden: true
# ---
#
# We have a sequence of `books`: the `i`-th book has thickness `books[i][0]` and
# height `books[i][1]`.
#
# We want to place these books **in order**  onto bookcase shelves that have
# total width `shelf_width`.
#
# We choose some of the books to place on this shelf (such that the sum of their
# thickness is `<= shelf_width`), then build another level of shelf of the
# bookcase so that the total height of the bookcase has increased by the maximum
# height of the books we just put down.  We repeat this process until there are
# no more books to place.
#
# Note again that at each step of the above process, _the order of the books we
# place is the same order as the given sequence of books_.  For example, if we
# have an ordered list of 5 books, we might place the first and second book onto
# the first shelf, the third book on the second shelf, and the fourth and fifth
# book on the last shelf.
#
# Return the minimum possible height that the total bookshelf can be after
# placing shelves in this manner.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/06/24/shelves.png)
#
#
#
#     Input: books = [[1,1],[2,3],[2,3],[1,1],[1,1],[1,1],[1,2]], shelf_width = 4
#     Output: 6
#     Explanation:
#     The sum of the heights of the 3 shelves are 1 + 3 + 2 = 6.
#     Notice that book number 2 does not have to be on the first shelf.
#
#
#
#
# **Constraints:**
#
#   * `1 <= books.length <= 1000`
#   * `1 <= books[i][0] <= shelf_width <= 1000`
#   * `1 <= books[i][1] <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export reshape_tensor, tensor_identity

"""
    function reshape_tensor(tensor, indices)

Permutes and reshapes tensor to move and merge indices
"""
function reshape_tensor(tensor, indices)
    index_order = Array{Int, 1}()
    dims = Array{Int, 1}()
    for x in indices
        dim = 1
        for y in x
            push!(index_order, y)
            dim *= size(tensor)[y]
        end
        push!(dims, dim)
    end
    tensor = permutedims(tensor, index_order)
    reshape(tensor, Tuple(dims))
end

function tensor_identity(m, n, d)
    id = zeros(m, n, d, d)
    id_small = Diagonal(ones(d))
    for i in 1:m
        for j in 1:n
            if i == j
                id[i, j, :, :] = id_small
            end
        end
    end
    id
end


text:
module SurfaceQG

export
  Problem,
  set_b!,
  updatevars!,

  kinetic_energy,
  buoyancy_variance,
  buoyancy_dissipation,
  buoyancy_work,
  buoyancy_advection,
  kinetic_energy_advection

  using
    CUDA,
    Reexport,
    DocStringExtensions

@reexport using FourierFlows

using LinearAlgebra: mul!, ldiv!
using FourierFlows: parsevalsum

nothingfunction(args...) = nothing

"""
    Problem(dev::Device=CPU();
                 nx = 256,
                 Lx = 2π,
                 ny = nx,
                 Ly = Lx,
                  ν = 0,
                 nν = 1,
                 dt = 0.01,
            stepper = "RK4",
              calcF = nothingfunction,
         stochastic = false,
  aliased_fraction = 1/3,
                  T = Float64)

Construct a surface quasi-geostrophic `problem` on device `dev`.

Keyword arguments
=================
  - `dev`: (required) `CPU()` or `GPU()`; computer architecture used to time-step `problem`.
  - `nx`: Number of grid points in ``x``-domain.
  - `ny`: Number of grid points in ``y``-domain.
  - `Lx`: Extent of the ``x``-domain.
  - `Ly`: Extent of the ``y``-domain.
  - `ν`: Small-scale (hyper)-viscosity coefficient.
  - `nν`: (Hyper)-viscosity order, `nν```≥ 1``.
  - `dt`: Time-step.
  - `stepper`: Time-stepping method.
  - `calcF`: Function that calculates the Fourier transform of the forcing, ``F̂``.
  - `stochastic`: `true` or `false`; boolean denoting whether `calcF` is temporally stochastic.
  - `aliased_fraction`: the fraction of high-wavenumbers that are zero-ed out by `dealias!()`.
  - `T`: `Float32` or `Float64`; floating point type used for `problem` data.
"""
function Problem(dev::Device=CPU();
  # Numerical parameters
                nx = 256,
                ny = nx,
                Lx = 2π,
                Ly = Lx,
  # Hyper-viscosity parameters
                 ν = 0,
                nν = 1,
  # Timestepper and equation options
                dt = 0.01,
           stepper = "RK4",
             calcF = nothingfunction,
        stochastic = false,
  # Float type and dealiasing
  aliased_fraction = 1/3,
                 T = Float64)

  grid = TwoDGrid(dev, nx, Lx, ny, Ly; aliased_fraction=aliased_fraction, T=T)

  params = Params{T}(ν, nν, calcF)

  vars = calcF == nothingfunction ? DecayingVars(dev, grid) : (stochastic ? StochasticForcedVars(dev, grid) : ForcedVars(dev, grid))

  equation = Equation(params, grid)

  return FourierFlows.Problem(equation, stepper, dt, grid, vars, params, dev)
end

# ----------
# Parameters
# ----------

"""
    Params{T}(ν, nν, calcF!)

A struct containing the parameters for Surface QG dynamics. Included are:

$(TYPEDFIELDS)
"""
struct Params{T} <: AbstractParams
    "buoyancy (hyper)-viscosity coefficient"
       ν :: T
    "buoyancy (hyper)-viscosity order"
      nν :: Int
    "function that calculates the Fourier transform of the forcing, ``F̂``"
  calcF! :: Function
end

Params(ν, nν) = Params(ν, nν, nothingfunction)

# ---------
# Equations
# ---------

"""
    Equation(params, grid)

Return the `equation` for surface QG dynamics with `params` and `grid`. The linear
opeartor ``L`` includes (hyper)-viscosity of order ``n_ν`` with coefficient ``ν``,

```math
L = - ν |𝐤|^{2 n_ν} .
```

Plain old viscocity corresponds to ``n_ν=1``.

The nonlinear term is computed via function `calcN!()`.
"""
function Equation(params::Params, grid::AbstractGrid)
  L = @. - params.ν * grid.Krsq^params.nν
  CUDA.@allowscalar L[1, 1] = 0

  return FourierFlows.Equation(L, calcN!, grid)
end

# ----
# Vars
# ----

abstract type SurfaceQGVars <: AbstractVars end

"""
    Vars{Aphys, Atrans, F, P}(b, u, v, bh, uh, vh, Fh, prevsol)

The variables for surface QG problem:

$(FIELDS)
"""
struct Vars{Aphys, Atrans, F, P} <: SurfaceQGVars
    "buoyancy"
        b :: Aphys
    "x-component of velocity"
        u :: Aphys
    "y-component of velocity"
        v :: Aphys
    "Fourier transform of buoyancy"
       bh :: Atrans
    "Fourier transform of x-component of velocity"
       uh :: Atrans
    "Fourier transform of y-component of velocity"
       vh :: Atrans
    "Fourier transform of forcing"
       Fh :: F
    "`sol` at previous time-step"
  prevsol :: P
end

const DecayingVars = Vars{<:AbstractArray, <:AbstractArray, Nothing, Nothing}
const ForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, Nothing}
const StochasticForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, <:AbstractArray}

"""
    DecayingVars(dev, grid)

Return the `vars` for unforced surface QG dynamics on device `dev` and with `grid`.
"""
function DecayingVars(::Dev, grid::AbstractGrid) where Dev
  T = eltype(grid)
  @devzeros Dev T (grid.nx, grid.ny) b u v
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) bh uh vh

  return Vars(b, u, v, bh, uh, vh, nothing, nothing)
end

"""
    ForcedVars(dev, grid)

Return the vars for forced surface QG dynamics on device `dev` and with `grid`.
"""
function ForcedVars(dev::Dev, grid) where Dev
  T = eltype(grid)
  @devzeros Dev T (grid.nx, grid.ny) b u v
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) bh uh vh Fh

  return Vars(b, u, v, bh, uh, vh, Fh, nothing)
end

"""
    StochasticForcedVars(dev, grid)

Return the `vars` for stochastically forced surface QG dynamics on device `dev` and with `grid`.
"""
function StochasticForcedVars(dev::Dev, grid) where Dev
  T = eltype(grid)

  @devzeros Dev T (grid.nx, grid.ny) b u v
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) bh uh vh Fh prevsol

  return Vars(b, u, v, bh, uh, vh, Fh, prevsol)
end

# -------
# Solvers
# -------

"""
    calcN_advection(N, sol, t, clock, vars, params, grid)

Calculate the Fourier transform of the advection term, ``- 𝖩(ψ, b)`` in conservative
form, i.e., ``- ∂_x[(∂_y ψ)b] - ∂_y[(∂_x ψ)b]`` and store it in `N`:

```math
N = - \\widehat{𝖩(ψ, b)} = - i k_x \\widehat{u b} - i k_y \\widehat{v b} .
```
"""
function calcN_advection!(N, sol, t, clock, vars, params, grid)
  @. vars.bh = sol
  @. vars.uh =   im * grid.l  * sqrt(grid.invKrsq) * sol
  @. vars.vh = - im * grid.kr * sqrt(grid.invKrsq) * sol

  ldiv!(vars.u, grid.rfftplan, vars.uh)
  ldiv!(vars.v, grid.rfftplan, vars.vh)
  ldiv!(vars.b, grid.rfftplan, vars.bh)

  ub, ubh = vars.u, vars.uh         # use vars.u, vars.uh as scratch variables
  vb, vbh = vars.v, vars.vh         # use vars.v, vars.vh as scratch variables

  @. ub *= vars.b                   # u*b
  @. vb *= vars.b                   # v*b

  mul!(ubh, grid.rfftplan, ub)      # \hat{u*b}
  mul!(vbh, grid.rfftplan, vb)      # \hat{v*b}

  @. N = - im * grid.kr * ubh - im * grid.l * vbh

  return nothing
end

"""
    calcN!(N, sol, t, clock, vars, params, grid)

Calculate the nonlinear term, that is the advection term and the forcing,

```math
N = - \\widehat{𝖩(ψ, b)} + F̂ .
```
"""
function calcN!(N, sol, t, clock, vars, params, grid)
  dealias!(sol, grid)

  calcN_advection!(N, sol, t, clock, vars, params, grid)

  addforcing!(N, sol, t, clock, vars, params, grid)

  return nothing
end

"""
    addforcing!(N, sol, t, clock, vars, params, grid)

When the problem includes forcing, calculate the forcing term ``F̂`` and add it to the
nonlinear term ``N``.
"""
addforcing!(N, sol, t, clock, vars::Vars, params, grid) = nothing

function addforcing!(N, sol, t, clock, vars::ForcedVars, params, grid)
  params.calcF!(vars.Fh, sol, t, clock, vars, params, grid)
  @. N += vars.Fh

  return nothing
end

function addforcing!(N, sol, t, clock, vars::StochasticForcedVars, params, grid)
  if t == clock.t # not a substep
    @. vars.prevsol = sol # sol at previous time-step is needed to compute budgets for stochastic forcing
    params.calcF!(vars.Fh, sol, t, clock, vars, params, grid)
  end

  @. N += vars.Fh

  return nothing
end

# ----------------
# Helper functions
# ----------------

"""
    updatevars!(prob)

Update variables in `vars` with solution in `sol`.
"""
function updatevars!(prob)
  vars, grid, sol = prob.vars, prob.grid, prob.sol

  dealias!(sol, grid)

  @. vars.bh = sol
  @. vars.uh =   im * grid.l  * sqrt(grid.invKrsq) * sol
  @. vars.vh = - im * grid.kr * sqrt(grid.invKrsq) * sol

  ldiv!(vars.b, grid.rfftplan, deepcopy(vars.bh))
  ldiv!(vars.u, grid.rfftplan, deepcopy(vars.uh))
  ldiv!(vars.v, grid.rfftplan, deepcopy(vars.vh))

  return nothing
end

"""
    set_b!(prob, b)

Set the solution `sol` as the transform of `b` and update all variables.
"""
function set_b!(prob, b)
  mul!(prob.sol, prob.grid.rfftplan, b)
  CUDA.@allowscalar prob.sol[1, 1] = 0 # zero domain average

  updatevars!(prob)

  return nothing
end

"""
    kinetic_energy(prob)

Return the domain-averaged surface kinetic energy. Since ``u² + v² = |{\\bf ∇} ψ|²``, we get
```math
\\int \\frac1{2} |{\\bf ∇} ψ|² \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} \\frac1{2} |𝐤|² |ψ̂|² .
```
In SQG, this is identical to half the domain-averaged surface buoyancy variance.
"""
@inline function kinetic_energy(prob)
  sol, vars, grid = prob.sol, prob.vars, prob.grid

  ψh = vars.uh                     # use vars.uh as scratch variable
  kinetic_energyh = vars.bh        # use vars.bh as scratch variable

  @. ψh = sqrt(grid.invKrsq) * sol
  @. kinetic_energyh = 1 / 2 * grid.Krsq * abs2(ψh)

  return 1 / (grid.Lx * grid.Ly) * parsevalsum(kinetic_energyh, grid)
end

"""
    buoyancy_variance(prob)

Return the buoyancy variance,
```math
\\int b² \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} |b̂|² .
```
In SQG, this is identical to the velocity variance (i.e., twice the domain-averaged kinetic
energy).
"""
@inline function buoyancy_variance(prob)
  sol, grid = prob.sol, prob.grid

  return 1 / (grid.Lx * grid.Ly) * parsevalsum(abs2.(sol), grid)
end

"""
    buoyancy_dissipation(prob)

Return the domain-averaged dissipation rate of surface buoyancy variance due
to small scale (hyper)-viscosity,
```math
2 ν (-1)^{n_ν} \\int b ∇^{2n_ν} b \\frac{𝖽x 𝖽y}{L_x L_y} = - 2 ν \\sum_{𝐤} |𝐤|^{2n_ν} |b̂|² ,
```
where ``ν`` the (hyper)-viscosity coefficient ``ν`` and ``nν`` the (hyper)-viscosity order.
In SQG, this is identical to twice the rate of kinetic energy dissipation
"""
@inline function buoyancy_dissipation(prob)
  sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid
  buoyancy_dissipationh = vars.uh         # use vars.uh as scratch variable

  @. buoyancy_dissipationh = 2 * params.ν * grid.Krsq^params.nν * abs2(sol)
  CUDA.@allowscalar buoyancy_dissipationh[1, 1] = 0

  return 1 / (grid.Lx * grid.Ly) * parsevalsum(buoyancy_dissipationh, grid)
end

"""
    buoyancy_work(prob)
    buoyancy_work(sol, vars, grid)

Return the domain-averaged rate of work of buoyancy variance by the forcing,
```math
\\int 2 b F \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} 2 b̂ F̂^* .
```
"""
@inline function buoyancy_work(sol, vars::ForcedVars, grid)
  buoyancy_workh = vars.uh         # use vars.uh as scratch variable

  @. buoyancy_workh =  2 * sol * conj(vars.Fh)    # 2*b̂*conj(f̂)
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(buoyancy_workh, grid)
end

@inline function buoyancy_work(sol, vars::StochasticForcedVars, grid)
  buoyancy_workh = vars.uh         # use vars.uh as scratch variable

  @. buoyancy_workh =  (vars.prevsol + sol) * conj(vars.Fh) # Stratonovich
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(buoyancy_workh, grid)
end

@inline buoyancy_work(prob) = buoyancy_work(prob.sol, prob.vars, prob.grid)

end # module


text:
using Test, VaderSentiment

struct TestTuple
    sentence::AbstractString
    original_sentiment::Dict
end

test_set =  [TestTuple("VADER is smart, handsome, and funny.", Dict("neg"=> 0.0, "neu"=> 0.254, "pos"=> 0.746, "compound"=> 0.8316)),
             TestTuple("VADER is smart, handsome, and funny!", Dict("neg"=> 0.0, "neu"=> 0.248, "pos"=> 0.752, "compound"=> 0.8439)),
             TestTuple("VADER is very smart, handsome, and funny.", Dict("neg"=> 0.0, "neu"=> 0.299, "pos"=> 0.701, "compound"=> 0.8545)),
             TestTuple("VADER is VERY SMART, handsome, and FUNNY.", Dict("neg"=> 0.0, "neu"=> 0.246, "pos"=> 0.754, "compound"=> 0.9227)),
             TestTuple("VADER is VERY SMART, handsome, and FUNNY!!!", Dict("neg"=> 0.0, "neu"=> 0.233, "pos"=> 0.767, "compound"=> 0.9342)),
             TestTuple("VADER is VERY SMART, uber handsome, and FRIGGIN FUNNY!!!", Dict("neg"=> 0.0, "neu"=> 0.294, "pos"=> 0.706, "compound"=> 0.9469)),
             TestTuple("VADER is not smart, handsome, nor funny.", Dict("neg"=> 0.646, "neu"=> 0.354, "pos"=> 0.0, "compound"=> -0.7424)),
             TestTuple("The book was good.", Dict("neg"=> 0.0, "neu"=> 0.508, "pos"=> 0.492, "compound"=> 0.4404)),
             TestTuple("At least it isn't a horrible book.", Dict("neg"=> 0.0, "neu"=> 0.678, "pos"=> 0.322, "compound"=> 0.431)),
             TestTuple("The book was only kind of good.", Dict("neg"=> 0.0, "neu"=> 0.697, "pos"=> 0.303, "compound"=> 0.3832)),
             TestTuple("The plot was good, but the characters are uncompelling and the dialog is not great.", Dict("neg"=> 0.327, "neu"=> 0.579, "pos"=> 0.094, "compound"=> -0.7042)),
             TestTuple("Today SUX!", Dict("neg"=> 0.779, "neu"=> 0.221, "pos"=> 0.0, "compound"=> -0.5461)),
             TestTuple("Today only kinda sux! But I'll get by, lol", Dict("neg"=> 0.127, "neu"=> 0.556, "pos"=> 0.317, "compound"=> 0.5249)),
             TestTuple("Make sure you :) or :D today!", Dict("neg"=> 0.0, "neu"=> 0.294, "pos"=> 0.706, "compound"=> 0.8633)),
             TestTuple("Catch utf-8 emoji such as such as 💘 and 💋 and 😁", Dict("neg"=> 0.0, "neu"=> 0.615, "pos"=> 0.385, "compound"=> 0.875)),
             TestTuple("Not bad at all", Dict("neg"=> 0.0, "neu"=> 0.513, "pos"=> 0.487, "compound"=> 0.431)),
             TestTuple("Not GREATLY bad at all!!!!!", Dict("neg"=> 0.0, "neu"=> 0.456, "pos"=> 0.544, "compound"=> 0.6982)),
             TestTuple("Not GREATLY bad at all??", Dict("neg"=> 0.0, "neu"=> 0.502, "pos"=> 0.498, "compound"=> 0.6084)),
             TestTuple("Not GREATLY bad at all????", Dict("neg"=> 0.0, "neu"=> 0.467, "pos"=> 0.533, "compound"=> 0.6777)),
             TestTuple("Not least GREATLY bad at all????", Dict("neg"=> 0.0, "neu"=> 0.523, "pos"=> 0.477, "compound"=> 0.6777)),
             TestTuple("Not GREATLY least bad at all????", Dict("neg"=> 0.436, "neu"=> 0.564, "pos"=> 0.0, "compound"=> -0.5944)),
             TestTuple("Catch utf-8 emoji such as such as💘 and 💋 and 😁", Dict("neg"=> 0.0, "neu"=> 0.615, "pos"=> 0.385, "compound"=> 0.875)),
             TestTuple("least bad at all????", Dict("neg"=> 0.0, "neu"=> 0.441, "pos"=> 0.559, "compound"=> 0.5873)),
             TestTuple("No not GREATLY least bad at all????", Dict("neg"=> 0.548, "neu"=> 0.452, "pos"=> 0.0, "compound"=> -0.7238)),
             TestTuple("The book was only kind of no good.", Dict("neg"=> 0.381, "neu"=> 0.619, "pos"=> 0.0, "compound"=> -0.4017)),
             TestTuple("The book was only kind of bad ass good.", Dict("neg"=> 0.246, "neu"=> 0.422, "pos"=> 0.331, "compound"=> 0.0534)),
             TestTuple("The book was only kind of never so bad ass good.", Dict("neg"=> 0.0, "neu"=> 0.516, "pos"=> 0.484, "compound"=> 0.7579)),
             TestTuple("The book was only kind of without doubt bad ass good.", Dict("neg"=> 0.0, "neu"=> 0.412, "pos"=> 0.588, "compound"=> 0.8406)),
             TestTuple("The book was only kind of badn't ass good.", Dict("neg"=> 0.196, "neu"=> 0.571, "pos"=> 0.233, "compound"=> 0.1139)),
             TestTuple("I am a big happy cat.", Dict("neg"=> 0.0, "neu"=> 0.575, "pos"=> 0.425, "compound"=> 0.5719))]

analyzer = VaderSentiment.SentimentIntensityAnalyzer

@testset "Test Vader Sentiment" begin
    for test_case in test_set
        @test (analyzer(test_case.sentence).polarity_scores) == test_case.original_sentiment end end

max_stars_repo_name: nusretipek/JuliaVader

max_issues_repo_name: nusretipek/JuliaVader


text:
#=

  Euler #40 in Julia.

  Problem 40:
  """
  An irrational decimal fraction is created by concatenating the positive integers:

  0.123456789101112131415161718192021...

  It can be seen that the 12th digit of the fractional part is 1.

  If dn represents the nth digit of the fractional part, find the
  value of the following expression.

  d1 × d10 × d100 × d1000 × d10000 × d100000 × d1000000
  """

  This Julia model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/Julia_progs/

=#

include("Euler.jl")

# 0.00941826s
function euler40a()
    i = 1
    dlen = 1
    p = 1
    index = 10; # Index = 10, 100, 1000, ..., 1000000
    while dlen <= 1000000
        i += 1
        istr = string(i)
        istrlen = length(istr)
        if dlen+istrlen>=index
            p *= parse(Int,istr[index-dlen])
            index *= 10
        end
        dlen += istrlen
    end

    return p
end

run_euler(euler40a)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module LSH

####
#
## LSH.jl
#
# This package implements Locality Sensitive Hashing as desribed in [AI06].
# The implementation is guided by the description in [AM04]
#
## Basic Description of the algorithm
#
## Notations and conventions
#
# - Let 𝒫 be the data set we wish to query where each point is of
#   type T𝒫.
#
# - Let H <: HashFunction be the datastructure holding the necessary
#   data to perform a hash of data point yielding a value of type RH.
#   In particular, we need the invariant that
#   ∀ p, typeof(p) = T𝒫 -> H(t)::RH.
#   Implementations for RH === Uint32 are provided in this package.

## References
#
# [AI06] Alexandr Andoni and Piotr Indyk
#   "Near-Optimal Hashing Algorithms for Near Neighbor Problem
#    in High Dimensions"
#
# [AM04] Alexandr Andoni, Mayur Datar, Nicole Immorlica, Piotr Indyk,
#           and Vahab Mirrokni
#   "Locality-Sensitive Hashing Scheme Based on p-Stable Distributions"
#
##

using DataStructures
using ProgressMeter
import DataStructures:  MultiHashDict, AbstractHashDict, HashDict, Unordered, Ordered

abstract HashFunction

export dimension, HashFunction, LSHTable

# ============================= GroupingSet ====================================

import Base: push!, getindex, get
import DataStructures: rehash_item, hashindex

# Groups all elements added to it via `push!(x, k)` into a vector utilizing
# a hash function as the grouping key. All elements inserted into the vector
# so far can be retrieved by looking up any item.

immutable PointGroupingSet{V,GroupingF,RH,HashF,DeriveF,DK,T𝒫} <: AbstractHashDict{V,Vector{V},Unordered}
    group::GroupingF
    𝒫::Vector{T𝒫}
    dict::MultiHashDict{RH,V,Unordered,HashF,DeriveF,DK}
end

function PointGroupingSet{V,GroupingF,RH,HashF,DeriveF,DK,T𝒫}(
        𝒫::Vector{T𝒫},::Type{V},groupF::GroupingF,::Type{RH},hash::HashF,derive::DeriveF,::Type{DK})
    g = PointGroupingSet{V,GroupingF,RH,HashF,DeriveF,DK,T𝒫}(
        groupF,𝒫,MultiHashDict{RH,V,Unordered,HashF,DeriveF,DK}(hash,derive))
    sizehint(g.dict,2*length(𝒫))
    g
end

function push!{V,GroupingF,RH,HashF,DeriveF,DK,T𝒫}(g::PointGroupingSet{V,GroupingF,RH,HashF,DeriveF,DK,T𝒫},
        item_no, p=nothing)
    if is(p,nothing)
        DataStructures._push!(g,g.dict,g.group(g.𝒫[item_no])=>V(item_no))
    else
        DataStructures._push!(g,g.dict,g.group(p)=>V(item_no))
    end
end

rehash(g::PointGroupingSet, sz) = DataStructures._rehash(g,g.dict,sz)

# We can find out what the hash value with this finger print was by passing
# any element of the group back throught the hash
function rehash_item(g::PointGroupingSet, k, item, sz)
    new_key = g.group(g.𝒫[item])
    @assert isequal(k,g.dict.derive(new_key))
    ind = hashindex(g.dict, new_key, sz)
    ind
end

function getindex(g::PointGroupingSet, q)
    key = g.group(q)
    g.dict[key]
end

# ========================= LSH hashindex and derive ===========================

# In the notation of [AM04], we have hashindex == t_1 and derive == t_2

const P = UInt32(UInt64(2)^32 - 5)

immutable ModPHash{RT} <: HashFunction
    # TODO: Make NTuple
    r::Vector{Uint32}
end
dimension(T::ModPHash) = length(T.r)

function check_length(T::ModPHash, N)
    if length(T.r) != N
        error("""Hash function has $(length(T.r)) coefficients, but
                 was called on an input of dimension $N""")
    end
end

function call{RT}(T::ModPHash{RT},z::Vector{Int32})
    check_length(T,length(z))
    result::Uint32 = Uint32(0)
    for i in 1:length(z)
        r = z[i]*Int64(T.r[i])
        result = ((result + (r % P)) % P) % Uint32
    end
    result % RT
end

function call{RT}(T::ModPHash{RT},z::Int32)
    check_length(T,1)
    (widemul(z,T.r[1]) % P) % RT
end

# ============================== LSHTable ======================================

import Base: push!

# The main data structure on which we may perform query
# contains a collection of k BucketTables
immutable LSHTable{H,RH,DK,V,T𝒫,HashF,DeriveF}
    R::Float64
    𝒫::Vector{T𝒫}
    tables::Vector{PointGroupingSet{V,H,RH,HashF,DeriveF,DK,T𝒫}}
end

# Given a vector of hashes and a dataset, construct an LSHTable
function LSHTable{H,T𝒫,RH,DK}(R::Float64, hashes::Vector{H}, 𝒫::Vector{T𝒫},
        ::Type{RH}, ::Type{DK}; progress = false)

    if progress
        p = Progress(length(𝒫), 1)
    end

    @show typeof(𝒫)

    # Set up the tables.
    tables = [
        (pgs = PointGroupingSet(
            𝒫,
            Int32,
            h,
            RH,
            ModPHash{Uint32}(rand(Uint32,dimension(h))), # t_1 in [AM04]
            ModPHash{DK}(rand(Uint32,dimension(h))), # t_2 in [AM04]
            DK
        ); pgs)
        for h in hashes ]

    T = LSHTable(R,𝒫,tables)

    for (i,point) in enumerate(𝒫)
        push!(T,i,point)
        progress && next!(p)
    end

    T
end

function push!(T::LSHTable,i,v)
    v = precompute(T.tables[1].group,v)
    for t in T.tables
        push!(t,i,v)
    end
end

function LSHTable{H,T𝒫}(R::Float64,hashes::Vector{H}, datapoints::Vector{T𝒫}; progress = false)
    LSHTable(R,hashes,datapoints,RH(H),Uint32; progress=progress)
end

using Distances

function getindex{H,RH,DK,V,T𝒫,HashP,DeriveP}(
        T::LSHTable{H,RH,DK,V,T𝒫,HashP,DeriveP},q::T𝒫)
    results = ObjectIdDict()
    tried = falses(length(T.𝒫))
    qp = precompute(T.tables[1].group,q)
    for table in T.tables
        for p in table[qp]
            if !tried[p]
                point = T.𝒫[p]
                if (euclidean(point,q) <= T.R)
                    results[point] = p
                end
                tried[p] = true
            end
        end
    end
    #@show length(tried)
    collect(keys(results))
end

# =============================== LSH Family ===================================
using Distributions

import Distributions: Gaussian, Uniform
import Base: rand

# Implements the LSH family given in [AM04]

immutable AM04HashFamily
    d::Int
    w::Float64
    R::Float64
end

# Represents an h_{a,b} from [AM04]
immutable AM04Hash
    a::Vector{Float64}
    b::Float64
    R::Float64
end

function rand(family::AM04HashFamily)
    #println("Creating new hash")
    AM04Hash(
        rand(Gaussian(0,1.0/family.w),family.d),
        rand(Uniform(0,1)),
        family.R
    )
end

call(h::AM04Hash,v) = floor(Int32,(dot(h.a,v)/h.R + h.b))

RH(::Type{AM04Hash}) = Int32

# The g_i = {h_{1,i}...h_{1,k}}
immutable HashCollection{H}
    hs::Vector{AM04Hash}
end

RH{H}(::Type{HashCollection{H}}) = Vector{H}
dimension(c::HashCollection) = length(c.hs)

call(c::HashCollection,v) = [ h(v) for h in c.hs ]

function createHashesAM04(d,w,k,L,R)
    family = AM04HashFamily(d,w,R)
    [ HashCollection{AM04Hash}([rand(family) for _ = 1:k]) for _ = 1:L ]
end

# The hashing scheme define in section 3.4.1 of the E2LSH manual

immutable CompositeHashCollection{H}
    first::Uint8
    second::Uint8
    hs::Vector{HashCollection{H}}
end

RH{H}(::Type{CompositeHashCollection{H}}) = Vector{H}
dimension(c::CompositeHashCollection) = 2*dimension(first(c.hs))

immutable PrecomputedHashes{H}
    hs::Vector{HashCollection{H}}
    computed::Vector{Vector{Int32}}
end

call(c::CompositeHashCollection,v) = vcat(c.hs[c.first](v),c.hs[c.second](v))
function call(c::CompositeHashCollection,v::PrecomputedHashes)
    @assert is(c.hs,v.hs)
    vcat(v.computed[c.first],v.computed[c.second])
end

function precompute{H}(hs::Vector{HashCollection{H}},q)
    PrecomputedHashes(hs,[h(q) for h in hs])
end
precompute(h::CompositeHashCollection,q) = precompute(h.hs,q)
precompute(_,q) = q

function createUHashesAM04(d,w,k,L,R,m)
    family = AM04HashFamily(d,w,R)
    hs = [ HashCollection{AM04Hash}([rand(family) for _ = 1:(k/2)]) for _ = 1:m ]
    r = Array(CompositeHashCollection{AM04Hash},div(m*(m-1),2))
    z = 1
    for i = 1:m
        for j = (i+1):m
            r[z] = CompositeHashCollection{AM04Hash}(i,j,hs)
            z += 1
        end
    end
    @assert z-1 == length(r)
    r
end

# Implements the LSH family given in [AI06]

end # module


text:
module PyCall

export pyinitialize, pyfinalize, pycall, pyimport, pybuiltin, PyObject,
       pysym, PyPtr, pyincref, pydecref, pyversion, PyArray, PyArray_Info,
       pyerr_check, pyerr_clear, pytype_query, PyAny, @pyimport, PyWrapper,
       PyDict, pyisinstance, pywrap, @pykw, pytypeof, pyeval, pyhassym,
       PyVector, pystring, pyraise

import Base.size, Base.ndims, Base.similar, Base.copy, Base.ref, Base.assign,
       Base.stride, Base.convert, Base.pointer, Base.summary, Base.convert,
       Base.show, Base.has, Base.keys, Base.values, Base.eltype, Base.get,
       Base.delete!, Base.empty!, Base.length, Base.isempty, Base.start,
       Base.done, Base.next, Base.filter!, Base.hash, Base.delete!, Base.pop!

#########################################################################

# Mirror of C PyObject struct (for non-debugging Python builds).
# We won't actually access these fields directly; we'll use the Python
# C API for everything.  However, we need to define a unique Ptr type
# for PyObject*, and we might as well define the actual struct layout
# while we're at it.
immutable PyObject_struct
    ob_refcnt::Int
    ob_type::Ptr{Void}
end

typealias PyPtr Ptr{PyObject_struct} # type for PythonObject* in ccall

#########################################################################

# Global configuration variables.  Note that, since Julia does not allow us
# to attach type annotations to globals, we need to annotate these explicitly
# as initialized::Bool and libpython::Ptr{Void} when we use them.
initialized = false # whether Python is initialized
finalized = false # whether Python has been finalized
libpython = C_NULL # Python shared library (from dlopen)

pysym(func::Symbol) = dlsym(libpython::Ptr{Void}, func)
pysym_e(func::Symbol) = dlsym_e(libpython::Ptr{Void}, func)
pyhassym(func::Symbol) = pysym_e(func) != C_NULL

# call pysym_e on the arguments and return the first non-NULL result
function pysym_e(funcs...)
    for func in funcs
        p = pysym_e(func)
        if p != C_NULL
            return p
        end
    end
    return C_NULL
end

# Macro version of pysym to cache dlsym lookup (thanks to vtjnash)
macro pysym(func)
    z = gensym(string(func))
    @eval global $z = C_NULL
    quote begin
        global $z
        if $z::Ptr{Void} == C_NULL
            $z::Ptr{Void} = dlsym(libpython::Ptr{Void}, $(esc(func)))
        end
        $z::Ptr{Void}
    end end
end

# Macro version of pyinitialize() to inline initialized? check
macro pyinitialize()
    :(initialized::Bool ? nothing : pyinitialize())
end

#########################################################################
# Wrapper around Python's C PyObject* type, with hooks to Python reference
# counting and conversion routines to/from C and Julia types.

type PyObject
    o::PyPtr # the actual PyObject*
    function PyObject(o::PyPtr)
        po = new(o)
        finalizer(po, pydecref)
        return po
    end
    PyObject() = PyObject(convert(PyPtr, C_NULL))
end

function pydecref(o::PyObject)
    if initialized::Bool # don't decref after pyfinalize!
        ccall((@pysym :Py_DecRef), Void, (PyPtr,), o.o)
    end
    o.o = convert(PyPtr, C_NULL)
    o
end

function pyincref(o::PyObject)
    ccall((@pysym :Py_IncRef), Void, (PyPtr,), o)
    o
end

# doing an incref *before* creating a PyObject may safer in the
# case of borrowed references, to ensure that no exception or interrupt
# induces a double decref.
function pyincref(o::PyPtr)
    ccall((@pysym :Py_IncRef), Void, (PyPtr,), o)
    PyObject(o)
end

pyisinstance(o::PyObject, t::PyObject) =
  t.o != C_NULL && ccall((@pysym :PyObject_IsInstance), Cint, (PyPtr,PyPtr), o, t.o) == 1

pyisinstance(o::PyObject, t::Ptr{Void}) =
  t != C_NULL && ccall((@pysym :PyObject_IsInstance), Cint, (PyPtr,PyPtr), o, t) == 1

pyquery(q::Ptr{Void}, o::PyObject) =
  ccall(q, Cint, (PyPtr,), o) == 1

pytypeof(o::PyObject) = o.o == C_NULL ? throw(ArgumentError("NULL PyObjects have no Python type")) : pycall(TypeType, PyObject, o)

# conversion to pass PyObject as ccall arguments:
convert(::Type{PyPtr}, po::PyObject) = po.o

# use constructor for generic conversions to PyObject
convert(::Type{PyObject}, o) = PyObject(o)
PyObject(o::PyObject) = o

#########################################################################

inspect = PyObject() # inspect module, needed for module introspection
builtin = PyObject() # __builtin__ module, needed for pybuiltin

# Python has zillions of types that a function be, in addition to the FunctionType
# in the C API.  We have to obtain these at runtime and cache them in globals
BuiltinFunctionType = PyObject()
TypeType = PyObject() # "type" function
MethodType = PyObject()
MethodWrapperType = PyObject()
# also WrapperDescriptorType = type(list.__add__) and
#      MethodDescriptorType = type(list.append) ... is it worth detecting these?

# special function type used in NumPy and SciPy (if available)
ufuncType = PyObject()

# cache Python None
pynothing = PyObject()

# Python 2/3 compatibility: cache dlsym for renamed functions
pystring_fromstring = C_NULL
pystring_asstring = C_NULL
pystring_type = C_NULL
pyint_type = C_NULL
pyint_from_size_t = C_NULL
pyint_from_ssize_t = C_NULL
pyint_as_ssize_t = C_NULL

# cache ctypes.c_void_p type and function if available
c_void_p_Type = PyObject()
py_void_p = p::Ptr -> PyObject(uint(p))

# PyCObject_Check and PyCapsule_CheckExact are actually macros
# that check against PyCObject_Type and PyCapsule_Type global variables,
# which we cache if they are available:
PyCObject_Type = C_NULL
PyCapsule_Type = C_NULL

# Py_SetProgramName needs its argument to persist as long as Python does
pyprogramname = bytestring("")

# PyUnicode_* may actually be a #define for another symbol,
# so we cache the correct dlsym
PyUnicode_AsUTF8String = C_NULL
PyUnicode_DecodeUTF8 = C_NULL

# whether to use unicode for strings by default, ala Python 3
pyunicode_literals = false

# traceback.format_tb function, for show(PyError)
format_traceback = PyObject()

# cache whether Python hash values are Clong (Python < 3.2) or Int (>= 3.2)
pyhashlong = false

# Return the Python version as a Julia VersionNumber
pyversion = v"0"

# low-level initialization, given a pointer to dlopen result on libpython,
# or C_NULL if python symbols are in the global namespace:
# initialize the Python interpreter (no-op on subsequent calls)
function pyinitialize(libpy::Ptr{Void})
    global initialized
    global finalized
    global libpython
    global pyprogramname
    global inspect
    global BuiltinFunctionType
    global TypeType
    global MethodType
    global MethodWrapperType
    global ufuncType
    global PyUnicode_AsUTF8String
    global PyUnicode_DecodeUTF8
    global pyunicode_literals
    global pynothing
    global pyhashlong
    global pystring_fromstring
    global pystring_asstring
    global pystring_type
    global pyint_type
    global pyint_from_size_t
    global pyint_from_ssize_t
    global pyint_as_ssize_t
    global pyversion
    global PyCObject_Type
    global PyCapsule_Type
    global c_void_p_Type
    global py_void_p
    global format_traceback
    if !initialized::Bool
        if finalized::Bool
            # From the Py_Finalize documentation:
            #    "Some extensions may not work properly if their
            #     initialization routine is called more than once; this
            #     can happen if an application calls Py_Initialize()
            #     and Py_Finalize() more than once."
            # For example, numpy and scipy seem to crash if this is done.
            error("Calling pyinitialize after pyfinalize is not supported")
        end
        libpython::Ptr{Void} = libpy == C_NULL ? ccall(:jl_load_dynamic_library, Ptr{Void}, (Ptr{Uint8},Cuint), C_NULL, 0) : libpy
        already_inited = 0 != ccall((@pysym :Py_IsInitialized), Cint, ())
        if !already_inited
            if !isempty(pyprogramname::ASCIIString)
                ccall((@pysym :Py_SetProgramName), Void, (Ptr{Uint8},),
                      pyprogramname::ASCIIString)
            end
            ccall((@pysym :Py_InitializeEx), Void, (Cint,), 0)
        end
        initialized::Bool = true
        inspect::PyObject = pyimport("inspect")
        types = pyimport("types")
        BuiltinFunctionType::PyObject = types["BuiltinFunctionType"]
        TypeType::PyObject = pybuiltin("type")
        MethodType::PyObject = types["MethodType"]
        MethodWrapperType::PyObject = pytypeof(PyObject(PyObject[])["__add__"])
        try
            ufuncType = pyimport("numpy")["ufunc"]
        catch
            ufuncType = PyObject() # NumPy not available
        end
        PyUnicode_AsUTF8String::Ptr{Void} =
          pysym_e(:PyUnicode_AsUTF8String,
                  :PyUnicodeUCS4_AsUTF8String,
                  :PyUnicodeUCS2_AsUTF8String)
        PyUnicode_DecodeUTF8::Ptr{Void} =
          pysym_e(:PyUnicode_DecodeUTF8,
                  :PyUnicodeUCS4_DecodeUTF8,
                  :PyUnicodeUCS2_DecodeUTF8)
        pynothing = pyincref(convert(PyPtr, pysym(:_Py_NoneStruct)))
        if pyhassym(:PyString_FromString)
            pystring_fromstring::Ptr{Void} = pysym(:PyString_FromString)
            pystring_asstring::Ptr{Void} = pysym(:PyString_AsString)
            pystring_type::Ptr{Void} = pysym(:PyString_Type)
        else
            pystring_fromstring::Ptr{Void} = pysym(:PyBytes_FromString)
            pystring_asstring::Ptr{Void} = pysym(:PyBytes_AsString)
            pystring_type::Ptr{Void} = pysym(:PyBytes_Type)
        end
        if pyhassym(:PyInt_Type)
            pyint_type::Ptr{Void} = pysym(:PyInt_Type)
            pyint_from_size_t::Ptr{Void} = pysym(:PyInt_FromSize_t)
            pyint_from_ssize_t::Ptr{Void} = pysym(:PyInt_FromSsize_t)
            pyint_as_ssize_t::Ptr{Void} = pysym(:PyInt_AsSsize_t)
        else
            pyint_type::Ptr{Void} = pysym(:PyLong_Type)
            pyint_from_size_t::Ptr{Void} = pysym(:PyLong_FromSize_t)
            pyint_from_ssize_t::Ptr{Void} = pysym(:PyLong_FromSsize_t)
            pyint_as_ssize_t::Ptr{Void} = pysym(:PyLong_AsSsize_t)
        end
        PyCObject_Type::Ptr{Void} = pysym_e(:PyCObject_Type)
        PyCapsule_Type::Ptr{Void} = pysym_e(:PyCapsule_Type)
        try
            c_void_p_Type::PyObject = pyimport("ctypes")["c_void_p"]
            py_void_p::Function = p::Ptr -> pycall(c_void_p_Type::PyObject, PyObject, uint(p))
        catch # fallback to CObject
            pycobject_new = pysym(:PyCObject_FromVoidPtr)
            py_void_p::Function = p::Ptr -> PyObject(ccall(pycobject_new, PyPtr, (Ptr{Void}, Ptr{Void}), p, C_NULL))
        end
        pyversion::VersionNumber =
          VersionNumber(convert((Int,Int,Int,String,Int),
                                pyimport("sys")["version_info"])[1:3]...)
        pyhashlong::Bool = pyversion::VersionNumber < v"3.2"
        pyunicode_literals::Bool = pyversion::VersionNumber >= v"3.0"
        format_traceback::PyObject = pyimport("traceback")["format_tb"]
        pyexc_initialize()
        if !already_inited
            # some modules (e.g. IPython) expect sys.argv to be set
            argv_s = bytestring("")
            argv = convert(Ptr{Uint8}, argv_s)
            ccall(pysym(:PySys_SetArgvEx), Void, (Cint,Ptr{Ptr{Uint8}},Cint),
                  1, &argv, 0)
        end
    end
    return
end

pyconfigvar(python::String, var::String) = chomp(readall(`$python -c "import distutils.sysconfig; print(distutils.sysconfig.get_config_var('$var'))"`))

function libpython_name(python::String)
    lib = pyconfigvar(python, "LDLIBRARY")
    @osx_only if lib[1] != '/'; lib = pyconfigvar(python, "PYTHONFRAMEWORKPREFIX")*"/"*lib end
    lib
end

# initialize the Python interpreter (no-op on subsequent calls)
function pyinitialize(python::String)
    global initialized
    global pyprogramname
    if !initialized::Bool
        libpy = try
            dlopen(python, RTLD_LAZY|RTLD_DEEPBIND|RTLD_GLOBAL)
          catch
            dlopen(libpython_name(python), RTLD_LAZY|RTLD_DEEPBIND|RTLD_GLOBAL)
          end
        pyprogramname::ASCIIString = bytestring(python)
        pyinitialize(libpy)
    end
    return
end

pyinitialize() = pyinitialize("python") # default Python executable name
libpython_name() = libpython_name("python") # ditto

# end the Python interpreter and free associated memory
function pyfinalize()
    global initialized
    global finalized
    global libpython
    global inspect
    global builtin
    global BuiltinFunctionType
    global TypeType
    global MethodType
    global MethodWrapperType
    global ufuncType
    global pynothing
    global c_void_p_Type
    global format_traceback
    if initialized::Bool
        pydecref(ufuncType::PyObject)
        npyfinalize()
        pydecref(format_traceback::PyObject)
        pydecref(c_void_p_Type::PyObject)
        pydecref(pynothing::PyObject)
        pydecref(BuiltinFunctionType::PyObject)
        pydecref(TypeType::PyObject)
        pydecref(MethodType::PyObject)
        pydecref(MethodWrapperType::PyObject)
        pydecref(builtin::PyObject)
        pydecref(inspect::PyObject)
        pyexc_finalize()
        pygc_finalize()
        gc() # collect/decref any remaining PyObjects
        ccall((@pysym :Py_Finalize), Void, ())
        dlclose(libpython::Ptr{Void})
        libpython::Ptr{Void} = C_NULL
        initialized::Bool = false
        finalized::Bool = true
    end
    return
end

#########################################################################

include("exception.jl")

#########################################################################

include("gc.jl")

# make a PyObject that embeds a reference to keep, to prevent Julia
# from garbage-collecting keep until o is finalized.
PyObject(o::PyPtr, keep::Any) = pyembed(PyObject(o), keep)

#########################################################################

include("conversions.jl")

include("pytype.jl")

include("callback.jl")

#########################################################################
# Pretty-printing PyObject

function pystring(o::PyObject)
    if o.o == C_NULL
        return "NULL"
    else
        s = ccall((@pysym :PyObject_Repr), PyPtr, (PyPtr,), o)
        if (s == C_NULL)
            pyerr_clear()
            s = ccall((@pysym :PyObject_Str), PyPtr, (PyPtr,), o)
            if (s == C_NULL)
                pyerr_clear()
                return string(o.o)
            end
        end
        return convert(String, PyObject(s))
    end
end

function show(io::IO, o::PyObject)
    print(io, "PyObject $(pystring(o))")
end

#########################################################################
# computing hashes of PyObjects

pysalt = hash("PyCall.PyObject") # "salt" to mix in to PyObject hashes

function hash(o::PyObject)
    if o.o == C_NULL
        bitmix(pysalt::Uint, hash(C_NULL))
    elseif is_pyjlwrap(o)
        # call native Julia hash directly on wrapped Julia objects,
        # since on 64-bit Windows the Python 2.x hash is only 32 bits
        bitmix(pysalt::Uint, hash(unsafe_pyjlwrap_to_objref(o.o)))
    else
        h = pyhashlong::Bool ? # changed to Py_hash_t in Python 3.2
               ccall((@pysym :PyObject_Hash), Clong, (PyPtr,), o) :
               ccall((@pysym :PyObject_Hash), Int, (PyPtr,), o)
        if h == -1 # error
            pyerr_clear()
            return bitmix(pysalt::Uint, hash(o.o))
        end
        bitmix(pysalt::Uint, uint(h))
    end
end

#########################################################################
# For o::PyObject, make o["foo"] and o[:foo] equivalent to o.foo in Python,
# with the former returning an raw PyObject and the latter giving the PyAny
# conversion.

function ref(o::PyObject, s::String)
    if (o.o == C_NULL)
        throw(ArgumentError("ref of NULL PyObject"))
    end
    p = ccall((@pysym :PyObject_GetAttrString), PyPtr,
              (PyPtr, Ptr{Uint8}), o, bytestring(s))
    if p == C_NULL
        pyerr_clear()
        throw(KeyError(s))
    end
    return PyObject(p)
end

ref(o::PyObject, s::Symbol) = convert(PyAny, ref(o, string(s)))

function assign(o::PyObject, v, s::String)
    if (o.o == C_NULL)
        throw(ArgumentError("assign of NULL PyObject"))
    end
    if -1 == ccall((@pysym :PyObject_SetAttrString), Cint,
                   (PyPtr, Ptr{Uint8}, PyPtr), o, bytestring(s), PyObject(v))
        pyerr_clear()
        throw(KeyError(s))
    end
    o
end

assign(o::PyObject, v, s::Symbol) = assign(o, v, string(s))

#########################################################################
# Create anonymous composite w = pywrap(o) wrapping the object o
# and providing access to o's members (converted to PyAny) as w.member.

abstract PyWrapper

# still provide w["foo"] low-level access to unconverted members:
ref(w::PyWrapper, s) = ref(w.___jl_PyCall_PyObject___, s)

typesymbol(T::DataType) = T.name.name
typesymbol(T) = :Any # punt

function pywrap(o::PyObject)
    @pyinitialize
    members = convert(Vector{(String,PyObject)},
                      pycall(inspect["getmembers"], PyObject, o))
    tname = gensym("PyCall_PyWrapper")
    @eval begin
        $(Expr(:type, true, Expr(:<:, tname, :PyWrapper),
               Expr(:block, :(___jl_PyCall_PyObject___::PyObject),
                    map(m -> Expr(:(::), symbol(m[1]),
                                  typesymbol(pytype_query(m[2]))),
                        members)...)))
        $(Expr(:call, tname, o,
               [ :(convert(PyAny, $(members[i][2])))
                for i = 1:length(members) ]...))
    end
end

#########################################################################

pyimport(name::String) =
    PyObject(@pycheckn ccall((@pysym :PyImport_ImportModule), PyPtr,
                             (Ptr{Uint8},), bytestring(name)))

pyimport(name::Symbol) = pyimport(string(name))

# convert expressions like :math or :(scipy.special) into module name strings
modulename(s::Symbol) = string(s)
function modulename(e::Expr)
    if e.head == :.
        string(modulename(e.args[1]), :., modulename(e.args[2]))
    elseif e.head == :quote
        modulename(e.args...)
    else
        throw(ArgumentError("invalid module"))
    end
end

macro pyimport(name, optional_varname...)
    mname = modulename(name)
    len = length(optional_varname)
    Name = len > 0 && (len != 2 || optional_varname[1] != :as) ?
      throw(ArgumentError("usage @pyimport module [as name]")) :
      (len == 2 ? optional_varname[2] :
       typeof(name) == Symbol ? name :
       throw(ArgumentError("$mname is not a valid module variable name, use @pyimport $mname as <name>")))
    quote
        $(esc(Name)) = pywrap(pyimport($mname))
        nothing
    end
end

#########################################################################

# look up a global builtin
function pybuiltin(name)
    global builtin
    if (builtin::PyObject).o == C_NULL
        builtin::PyObject = try
            pyimport("__builtin__")
        catch
            pyimport("builtins") # renamed in Python 3
        end
    end
    (builtin::PyObject)[name]
end

#########################################################################
# Keyword arguments are just passed to pycall as @pykw kw1=val1 kw2=val2...
# in the last pycall argument.  The @pykw macro converts this into a special
# dictionary object for passing to pycall.
#
# Note that I don't use the Options package because a lot of the things
# like typechecking and defaults don't make sense the pycall context,
# so this can be a lot simpler.  All Python needs is a dictionary.

# In order for PyCall to differentiate between the keyword dictionary
# and an ordinary dictionary parameter, we wrap it in a special type.
type PyKW
    d::Dict{String,Any} # dictionary of (keyword => value) pairs.
end
PyObject(kw::PyKW) = isempty(kw.d) ? PyObject() : PyObject(kw.d)

# from a single :(x = y) expression make a single :(x => y) expression
function pykw1(ex::Expr)
    if ex.head == :(=) && length(ex.args) == 2 && typeof(ex.args[1]) == Symbol
        Expr(:(=>), string(ex.args[1]), ex.args[2])
    else
        throw(ArgumentError("invalid keyword expression $ex"))
    end
end

macro pykw(args...)
    :(PyKW($(Expr(:typed_dict, :(String=>Any), map(pykw1, args)...))))
end

#########################################################################

typealias TypeTuple Union(Type,NTuple{Type})

function pycall(o::PyObject, returntype::TypeTuple, args...)
    oargs = map(PyObject, args)
    nargs = length(args)
    if nargs > 0 && isa(args[end], PyKW)
        kw = PyObject(args[end])
        nargs -= 1
    else
        kw = PyObject()
    end
    arg = PyObject(@pycheckn ccall((@pysym :PyTuple_New), PyPtr, (Int,),
                                   nargs))
    for i = 1:nargs
        @pycheckzi ccall((@pysym :PyTuple_SetItem), Cint, (PyPtr,Int,PyPtr),
                         arg, i-1, oargs[i])
        pyincref(oargs[i]) # PyTuple_SetItem steals the reference
    end
    ret = PyObject(@pycheckni ccall((@pysym :PyObject_Call), PyPtr,
                                    (PyPtr,PyPtr,PyPtr), o, arg, kw))
    jret = convert(returntype, ret)
    return jret
end

#########################################################################

const Py_eval_input = 258 # from Python.h
const pyeval_fname = bytestring("PyCall.jl") # filename for pyeval

# evaluate a python string, returning PyObject, given a dictionary
# (string/symbol => value) of local variables to use in the expression
function pyeval_(s::String, locals::PyDict)
    sb = bytestring(s) # use temp var to prevent gc before we are done with o
    o = PyObject(@pycheckn ccall((@pysym :Py_CompileString), PyPtr,
                                  (Ptr{Uint8}, Ptr{Uint8}, Cint),
                                  sb, pyeval_fname, Py_eval_input))
    main = @pycheckni ccall((@pysym :PyImport_AddModule),
                           PyPtr, (Ptr{Uint8},),
                           bytestring("__main__"))
    maindict = @pycheckni ccall((@pysym :PyModule_GetDict), PyPtr, (PyPtr,),
                                main)
    PyObject(@pycheckni ccall((@pysym :PyEval_EvalCode),
                              PyPtr, (PyPtr, PyPtr, PyPtr),
                              o, maindict, locals))
end

pyeval(s::String, locals::PyDict, returntype::TypeTuple) =
   convert(returntype, pyeval_(s, locals))
pyeval(s::String, locals::PyDict) = pyeval(s, locals, PyAny)
pyeval(s::String, locals::Associative, returntype::TypeTuple) =
   pyeval(s, PyDict(locals), returntype)
pyeval(s::String, locals::Associative) =
   pyeval(s, PyDict(locals), PyAny)
pyeval(s::String) = pyeval(s, PyDict())

#########################################################################

end # module PyCall


text:
VERSION < v"0.7.0-beta2.199" && __precompile__()

"""
The CondaBinDeps module provides access to the [conda](http://conda.pydata.org/) packages
manager as a BinDeps provider, to install binary
dependencies of other Julia packages.

To use Anaconda as a binary provider for BinDeps, the `CondaBinDeps.Manager` type is proposed. A
small example looks like this:

```julia
# Declare dependency
using BinDeps
@BinDeps.setup
netcdf = library_dependency("netcdf", aliases = ["libnetcdf","libnetcdf4"])

using CondaBinDeps
#  Use alternative conda channel.
CondaBinDeps.Conda.add_channel("my_channel")
provides(CondaBinDeps.Manager, "libnetcdf", netcdf)
```
"""
module CondaBinDeps
import Conda
using BinDeps

include("bindeps_conda.jl")

end


text:
@json struct User
    account_id::String
    account_status::String
    created_on::DateTime
    display_name::String
    has_2fa_enabled::Bool
    is_staff::Bool
    links::NamedTuple
    location::String
    nickname::String
    type::String
    username::String
    uuid::UUID
    website::String
end
@json_struct User

@json struct WorkspaceMembership
    links::NamedTuple
    user::User
end
@json_struct WorkspaceMembership

endpoint(::BitbucketAPI, ::typeof(get_user)) = Endpoint(:GET, "/user")
endpoint(::BitbucketAPI, ::typeof(get_user), name::AStr) =
    Endpoint(:GET, "/workspaces/$name/members",
             query = Dict(
                 :q=> "user.username=\"$name\""
             )
             )
@not_implemented(::BitbucketAPI, ::typeof(get_user), ::UUID)
postprocessor(::BitbucketAPI, ::typeof(get_user)) = DoSomething() do resp
    if match(r".*/workspaces/", resp.request.target) !== nothing
        StructTypes.constructfrom(Page{get_user, WorkspaceMembership}, JSON2.read(String(resp.body))).values[1].user
    else
        StructTypes.constructfrom(User, JSON2.read(String(resp.body)))
    end
end
into(::BitbucketAPI, ::typeof(get_user)) = User

endpoint(api::BitbucketAPI, ::typeof(get_users), workspace::AStr) = Endpoint(:GET, "/workspaces/$workspace/members")
into(::BitbucketAPI, ::typeof(get_users)) = Page{get_users, User}

@not_implemented(api::BitbucketAPI, ::typeof(update_user))
@not_implemented(api::BitbucketAPI, ::typeof(update_user), id::UUID)

@not_implemented(api::BitbucketAPI, ::typeof(create_user))

@not_implemented(api::BitbucketAPI, ::typeof(delete_user), id::UUID)


text:
using DimensionalData, Extents, Test

a = zeros(X(10.0:10.0:100.0), Y(0.1:0.1:1.0))
st = DimStack((a = a, b=a[Y=1]))
@test extent(st) == extent(a) == extent(dims(a)) == Extent(X=(10.0, 100.0), Y=(0.1, 1.0))
@test extent(st, X) == extent(a, X) == extent(dims(a), X) == Extent(; X=(10.0, 100.0))
@test dims(extent(st)) == dims(extent(a)) == (X((10.0, 100.0)), Y((0.1, 1.0)))
@test dims(extent(a), Y) == Y((0.1, 1.0))


text:
import ThorinDistributions
using Test

@testset "UnivariateGammaConvolution.jl tests" begin

    dist = ThorinDistributions.UnivariateGammaConvolution([1,0.5, 3.7],[4,2, 10])

    @testset "functions are working and consistant with past selves." begin
        spl = [
            44.21448353310021,
            40.67778979017133,
            38.67640776714249,
            23.525094134635108,
            23.38644058822338,
            29.498656708345194,
            35.340080215403574,
            56.380727733998995,
            25.52972573138749,
            51.20505953884155
            ]
            x = pdf.(dist,spl)
            x_fixed = [
            0.018757112919480665,
            0.020559416575509366,
            0.02140738452673978,
            0.019147239099563968,
            0.019031193841015707,
            0.02236341049555192,
            0.02240803708392265,
            0.011765753186399592,
            0.020617044066218452,
            0.014700161862923838
            ]

            y = cdf.(dist,spl)
            y_fixed = [
            0.6082512562725385,
            0.5386338347308577,
            0.49661228006226665,
            0.16557008552908928,
            0.16292327691403613,
            0.29129135121181804,
            0.4233565994115752,
            0.7937252726673726,
            0.20549205556263866,
            0.725339366143653
            ]

        @test x ≈ x_fixed
        @test y ≈ y_fixed
    end

    @testset "results of pdf and cdf are consistant with R coga" begin
        # We check against COGA from R on this simple example.
        # > library(coga)
        # > options(digits=20)
        # > coga::dcoga(1:10,c(1,0.5,3.7),1/c(4,2,10))
        #  [1] 9.1697877166357581e-07 1.4306872316747290e-05 6.6823198132184901e-05
        #  [4] 1.9070137722253317e-04 4.1580841085495259e-04 7.6517030290220020e-04
        #  [7] 1.2530854177965063e-03 1.8849283018035432e-03 2.6579560778065099e-03
        # [10] 3.5626340599333940e-03
        # > coga::pcoga(1:10,c(1,0.5,3.7),1/c(4,2,10))
        #  [1] 1.8115256558630823e-07 5.8062605413247208e-06 4.1780076813316473e-05
        #  [4] 1.6327100296424765e-04 4.5700519074037604e-04 1.0364052739079920e-03
        #  [7] 2.0336306415139858e-03 3.5906428800766086e-03 5.8506245787791631e-03
        # [10] 8.9504969486005927e-03
        # >

        values = Float64.(1:10)
        p = pdf.(dist,values)
        c = cdf.(dist,values)
        pᵣ = [
        9.1697877166357581e-07,
        1.4306872316747290e-05,
        6.6823198132184901e-05,
        1.9070137722253317e-04,
        4.1580841085495259e-04,
        7.6517030290220020e-04,
        1.2530854177965063e-03,
        1.8849283018035432e-03,
        2.6579560778065099e-03,
        3.5626340599333940e-03
        ]
        cᵣ = [
        1.8115256558630823e-07,
        5.8062605413247208e-06,
        4.1780076813316473e-05,
        1.6327100296424765e-04,
        4.5700519074037604e-04,
        1.0364052739079920e-03,
        2.0336306415139858e-03,
        3.5906428800766086e-03,
        5.8506245787791631e-03,
        8.9504969486005927e-03
        ]
        @test round.(p,digits=10) ≈ round.(pᵣ,digits=10)
        @test round.(c,digits=10) ≈ round.(cᵣ,digits=10)
    end

end


text:
module Sp3Wrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const BufferedReader = @jimport java.io.BufferedReader
const Class = @jimport java.lang.Class
const Function = @jimport java.util.function.Function
const InputStream = @jimport java.io.InputStream
const JString = @jimport java.lang.String
const Map = @jimport java.util.Map
const Object = @jimport java.lang.Object
const SP3File = @jimport org.orekit.files.sp3.SP3File
const SP3File_SP3FileType = @jimport org.orekit.files.sp3.SP3File$SP3FileType
const SP3File_SP3OrbitType = @jimport org.orekit.files.sp3.SP3File$SP3OrbitType
const SP3File_TimeSystem = @jimport org.orekit.files.sp3.SP3File$TimeSystem
const SP3Parser = @jimport org.orekit.files.sp3.SP3Parser
const TimeScales = @jimport org.orekit.time.TimeScales

export SP3File
export SP3Parser
export add_satellite
export contains_satellite
export equals
export get_agency
export get_class
export get_coordinate_system
export get_data_used
export get_day_fraction
export get_epoch
export get_epoch_interval
export get_gps_week
export get_julian_day
export get_number_of_epochs
export get_orbit_type
export get_orbit_type_key
export get_satellite_count
export get_satellites
export get_seconds_of_week
export get_time_system
export get_type
export hash_code
export notify
export notify_all
export parse
export to_string
export wait

include("s_p3_file.jl")
include("s_p3_parser.jl")

end


text:
function get_inflation(outfile = "temp.xlsx")
  url = "https://www.philadelphiafed.org/-/media/research-and-data/real-time-center/survey-of-professional-forecasters/historical-data/inflation.xlsx?la=en"

  urldownload(url, parser = identity, save_raw = outfile)
  df = DataFrame(XLSX.readtable(outfile, "INFLATION")...)
  df = ifelse.(df .== "#N/A", missing, df)

  outfile == "temp.xlsx" ? rm(outfile) : nothing

  return df
end


text:
module Landau2P
export advance_distribution!

function advance_distribution!()
end

end


text:
export sample_free_goal, sample_free!

function sample_free_goal(P::MPProblem)
    v = sample_goal(P.goal, P.SS)
    while !is_free_state(v, P.CC, P.SS)
        v = sample_goal(P.goal, P.SS)
    end
    return v
end

function sample_free!(P::MPProblem, N::Integer, ensure_goal::Bool = true; goal_bias = 0.0, ensure_goal_ct = 5)
    N <= 0 && return volume(P.SS)  # TODO: do something principled about replanning with fewer samples (see note in fmtstar!)
    V = P.V.V
    W = Array{eltype(V)}(N)
    if length(V) > 0 && V[1] == P.init
        sample_count = 0
    else
        sample_count = 1
        W[1] = P.init
    end
    attempts = 0
    successes = 0
    for v in repeatedly(() -> sample_space(P.SS)) # TODO: batch sampling; even pre-allocated perhaps
        attempts += 1
        if is_free_state(v, P.CC, P.SS)
            sample_count += 1
            successes += 1
            if 0.0 < goal_bias && rand() < goal_bias
                W[sample_count] = sample_free_goal(P)
            else
                W[sample_count] = v
            end
            if sample_count == N
                break
            end
        end
    end
    if ensure_goal
        for i in 1:min(ensure_goal_ct, N-1)
            W[N+1-i] = sample_free_goal(P)
        end
    end
    P.V = addpoints(P.V, W)
    return volume(P.SS) #ci(BinomialTest(successes, attempts), .05, tail = :left)[2] * volume(P.SS)
end


text:
using Base.Test

using PEGParser
import PEGParser: Node, StandardCache
using Compat

function test_string1_ncache()
  @grammar grammar begin
    start = "a"
  end

  (ast, pos, error) = parse(grammar, "a")
  @test error === nothing
  @test length(ast.children) == 0
  @test ast.value == "a"
  @test ast.sym === :start
  @test pos == 2
end

function test_string1_cache()
  @grammar grammar begin
    start = "a"
  end

  (ast, pos, error) = parse(grammar, "a", cache=StandardCache())
  @test error === nothing
  @test length(ast.children) == 0
  @test ast.value == "a"
  @test ast.sym === :start
  @test pos == 2
end

function test_or1_ncache()
  @grammar grammar begin
    start = "a" | "b" | "c"
  end

  (ast, pos, error) = parse(grammar, "a")
  @test ast.value == "a"
  @test pos == 2

  (ast, pos, error) = parse(grammar, "b")
  @test ast.value == "b"
  @test pos == 2

  (ast, pos, error) = parse(grammar, "c")
  @test ast.value == "c"
  @test pos == 2

  (ast, pos, error) = parse(grammar, "d")
  @test error !== nothing
  @test pos == 1
end

function test_or1_cache()
  @grammar grammar begin
    start = "a" | "b" | "c"
  end

  (ast, pos, error) = parse(grammar, "a", cache=StandardCache())
  @test ast.value == "a"
  @test pos == 2

  (ast, pos, error) = parse(grammar, "b", cache=StandardCache())
  @test ast.value == "b"
  @test pos == 2

  (ast, pos, error) = parse(grammar, "c", cache=StandardCache())
  @test ast.value == "c"
  @test pos == 2

  (ast, pos, error) = parse(grammar, "d", cache=StandardCache())
  @test error !== nothing
  @test pos == 1
end

function test_reference1_ncache()
  @grammar grammar begin
    start = a
    a = "a"
  end

  (ast, pos, error) = parse(grammar, "a")
  @test ast.value == "a"

  (ast, pos, error) = parse(grammar, "b")
  @test error !== nothing
end

function test_reference1_cache()
  @grammar grammar begin
    start = a
    a = "a"
  end

  (ast, pos, error) = parse(grammar, "a", cache=StandardCache())
  @test ast.value == "a"

  (ast, pos, error) = parse(grammar, "b", cache=StandardCache())
  @test error !== nothing
end

function test_or2_ncache()
  @grammar grammar begin
    start = a | b
    a = "a"
    b = "b"
  end

  (ast, pos, error) = parse(grammar, "a")
  @test ast.value == "a"
  @test pos == 2

  (ast, pos, error) = parse(grammar, "b")
  @test ast.value == "b"
  @test pos == 2

  (ast, pos, error) = parse(grammar, "c")
  @test error !== nothing
end

function test_or2_cache()
  @grammar grammar begin
    start = a | b
    a = "a"
    b = "b"
  end

  (ast, pos, error) = parse(grammar, "a", cache=StandardCache())
  @test ast.value == "a"
  @test pos == 2

  (ast, pos, error) = parse(grammar, "b", cache=StandardCache())
  @test ast.value == "b"
  @test pos == 2

  (ast, pos, error) = parse(grammar, "c", cache=StandardCache())
  @test error !== nothing
end

function test_and1_ncache()
  @grammar grammar begin
    start = a + b
    a = "a"
    b = "b"
  end

  (ast, pos, error) = parse(grammar, "ab")
  @test ast.value == "ab"
  @test length(ast.children) == 2
  @test ast.children[1].value == "a"
  @test ast.children[2].value == "b"
end

function test_and1_cache()
  @grammar grammar begin
    start = a + b
    a = "a"
    b = "b"
  end

  (ast, pos, error) = parse(grammar, "ab", cache=StandardCache())
  @test ast.value == "ab"
  @test length(ast.children) == 2
  @test ast.children[1].value == "a"
  @test ast.children[2].value == "b"
end

function test_andor_ncache()
  @grammar grammar begin
    start = a | b
    a = "a" + "b"
    b = "a" + "c"
  end

  (ast, pos, error) = parse(grammar, "ab")
  @test ast.value == "ab"

  (ast, pos, error) = parse(grammar, "ac")
  @test ast.value == "ac"
end

function test_andor_cache()
  @grammar grammar begin
    start = a | b
    a = "a" + "b"
    b = "a" + "c"
  end

  (ast, pos, error) = parse(grammar, "ab", cache=StandardCache())
  @test ast.value == "ab"

  (ast, pos, error) = parse(grammar, "ac", cache=StandardCache())
  @test ast.value == "ac"
end

function test_oneormore_ncache()
  @grammar grammar begin
    start = +("a")
  end

  (ast, pos, error) = parse(grammar, "")
  @test error !== nothing

  (ast, pos, error) = parse(grammar, "a")
  @test ast.value == "a"

  (ast, pos, error) = parse(grammar, "aa")
  @test ast.value == "aa"

  (ast, pos, error) = parse(grammar, "aaa")
  @test ast.value == "aaa"
end

function test_oneormore_cache()
  @grammar grammar begin
    start = +("a")
  end

  (ast, pos, error) = parse(grammar, "", cache=StandardCache())
  @test error !== nothing

  (ast, pos, error) = parse(grammar, "a", cache=StandardCache())
  @test ast.value == "a"

  (ast, pos, error) = parse(grammar, "aa", cache=StandardCache())
  @test ast.value == "aa"

  (ast, pos, error) = parse(grammar, "aaa", cache=StandardCache())
  @test ast.value == "aaa"
end

function test_zeroormore_ncache()
  @grammar grammar begin
    start = *("a")
  end

  (ast, pos, error) = parse(grammar, "b")
  @test ast === nothing

  (ast, pos, error) = parse(grammar, "a")
  @test ast.value == "a"

  (ast, pos, error) = parse(grammar, "aa")
  @test ast.value == "aa"

  (ast, pos, error) = parse(grammar, "aaa")
  @test ast.value == "aaa"
end

function test_zeroormore_cache()
  @grammar grammar begin
    start = *("a")
  end

  (ast, pos, error) = parse(grammar, "b", cache=StandardCache())
  @test ast === nothing

  (ast, pos, error) = parse(grammar, "a", cache=StandardCache())
  @test ast.value == "a"

  (ast, pos, error) = parse(grammar, "aa", cache=StandardCache())
  @test ast.value == "aa"

  (ast, pos, error) = parse(grammar, "aaa", cache=StandardCache())
  @test ast.value == "aaa"
end

function test_optional_ncache()
  @grammar grammar begin
    start = "a" + ?("b")
  end

  (ast, pos, error) = parse(grammar, "a")
  @test ast.value == "a"

  (ast, pos, error) = parse(grammar, "ab")
  @test ast.value == "ab"
end

function test_optional_cache()
  @grammar grammar begin
    start = "a" + ?("b")
  end

  (ast, pos, error) = parse(grammar, "a", cache=StandardCache())
  @test ast.value == "a"

  (ast, pos, error) = parse(grammar, "ab", cache=StandardCache())
  @test ast.value == "ab"
end

function test_regex_ncache()
  @grammar grammar begin
    start = r"[0-9]+"
  end

  (ast, pos, error) = parse(grammar, "0")
  @test ast.value == "0"

  (ast, pos, error) = parse(grammar, "42")
  @test ast.value == "42"
end

function test_regex_cache()
  @grammar grammar begin
    start = r"[0-9]+"
  end

  (ast, pos, error) = parse(grammar, "0", cache=StandardCache())
  @test ast.value == "0"

  (ast, pos, error) = parse(grammar, "42", cache=StandardCache())
  @test ast.value == "42"
end

function test_ambigious_ncache()
  @grammar grammar begin
    # ambigious on "c"
    start = "a" + ((Y + "e") | (?("b") + Y + "f"))
    Y = "c" + "d"
  end

  # first branch
  (ast, pos, error) = parse(grammar, "acde")
  @test ast.value == "acde"

  # second branch; skip optional "b"
  (ast, pos, error) = parse(grammar, "acdf")
  @test ast.value == "acdf"

  # second branch, use optional "b"
  (ast, pos, error) = parse(grammar, "abcdf")
  @test ast.value == "abcdf"
end

function test_ambigious_cache()
  @grammar grammar begin
    # ambigious on "c"
    start = "a" + ((Y + "e") | (?("b") + Y + "f"))
    Y = "c" + "d"
  end

  # first branch
  (ast, pos, error) = parse(grammar, "acde", cache=StandardCache())
  @test ast.value == "acde"

  # second branch; skip optional "b"
  (ast, pos, error) = parse(grammar, "acdf", cache=StandardCache())
  @test ast.value == "acdf"

  # We should traverse the first branch, bail because 'f' doesn't match
  # and then go up the second branch. At that point cache should be used.
  # Once a method to trace calls is created, can test this.
  # second branch, use optional "b"

  (ast, pos, error) = parse(grammar, "abcdf", cache=StandardCache())
  @test ast.value == "abcdf"
end

function test_list_ncache()
  @grammar grammar begin
    start = list(number, ",")
    number = r"[0-9]"
  end

  (ast, pos, error) = parse(grammar, "1,2,3")
  @test error === nothing
  @test ast.value == "1,2,3"

  (ast, pos, error) = parse(grammar, "")
  @test error !== nothing
end

function test_list0_ncache()
  @grammar grammar begin
    start = list(number, ",", 0)
    number = r"[0-9]"
  end

  (ast, pos, error) = parse(grammar, "1,2,3")
  @test error === nothing
  @test ast.value == "1,2,3"

  (ast, pos, error) = parse(grammar, "")
  @test error === nothing
  @test length(ast.children) == 0
end

function test_integer_parser()
  @grammar grammar begin
    start = integer
  end

  (ast, pos, error) = parse(grammar, "5")
  @test error === nothing
  @test ast.ruleType === IntegerRule
  @test ast.value == "5"
end

function test_float_parser()
  @grammar grammar begin
    start = float
  end

  (ast, pos, error) = parse(grammar, "5.123")
  @test error === nothing
  @test ast.ruleType === FloatRule
  @test ast.value == "5.123"
end

function test_semantic_action()
  @grammar grammar begin
    start = list(number, COMMA)
    number = (-SPACE + float){ parse(Float64, _1.value) } |
             (-SPACE + integer){ parse(Int64, _1.value) }
    COMMA = SPACE + ","
    SPACE = r"[ \t\n]*"
  end

  (ast, pos, error) = parse(grammar, "1,2.0,3,4.5")
  @test error === nothing
  @test length(ast.children) == 4
  @test ast.children[1] == 1
  @test ast.children[2] == 2.0
  @test ast.children[3] == 3
  @test ast.children[4] == 4.5
end

test_string1_ncache()
test_string1_cache()
test_or1_ncache()
test_or1_cache()
test_reference1_ncache()
test_reference1_cache()
test_or2_ncache()
test_or2_cache()
test_and1_ncache()
test_and1_cache()
test_andor_ncache()
test_andor_cache()
test_oneormore_ncache()
test_oneormore_cache()
test_zeroormore_ncache()
test_zeroormore_cache()
test_optional_ncache()
test_optional_cache()
test_regex_ncache()
test_regex_cache()
test_ambigious_ncache()
test_ambigious_cache()
test_list_ncache()
test_list0_ncache()

test_integer_parser()
test_float_parser()
test_semantic_action()


text:
# Type
"""
	(type) Ray{T<:AnyComplex} in the complex plane
Each `Ray` type is parameterized according to the common type of its complex input arguments.
"""
struct Ray{T<:AnyComplex} <: AbstractCurve
	base::T
	angle::AbstractFloat
	reverse::Bool
	function Ray{T}(a,d,rev=false) where T<:AnyComplex
		new(a,mod2pi(d),rev)
	end
end

# Untyped constructor
"""
	Ray(a,θ,reverse=false)
Construct the ray starting at `a` and extending to infinity at the angle `θ`. If `reverse` is true, the ray is considered to extend from infinity to `a` at angle `-θ`.
"""
function Ray(a::Number,θ::Real,rev=false)
	a = complex(float(a))
	Ray{typeof(a)}(a,float(θ),rev)
end

# Complex type converters
for ctype in [:Spherical,:Polar,:Complex]
	@eval begin
		function $ctype(R::Ray{T}) where T<:AnyComplex
			Ray($ctype(R.base),R.angle,R.reverse)
		end
	end
end

# Required methods
arclength(R::Ray) = Inf
function point(R::Ray{T},t::Real) where T
	if R.reverse
		t = 1-t
	end
	# avoid NaNs
	if t==0
		R.base
	elseif t==1
		T(Inf)
	else
		R.base + t/(1-t)*exp(complex(0,R.angle))
	end
end
(C::Ray)(t::Real) = point(C,t)

"""
	arg(R::Ray,z)
Find the parameter argument `t` such that `R(t)==z` is true.

This gives undefined results if `z` is not actually on the ray.
"""
function arg(R::Ray,z::Number)
	if isinf(z)
		t = 1
	else
		δ = abs(z - R.base)
		t = δ / (1+δ)
	end
	return R.reverse ? 1-t : t
end

function unittangent(R::Ray{T},t::Real=0) where T <: AnyComplex
	τ = T( exp(complex(0,R.angle)) )
	R.reverse ? -τ : τ
end

function tangent(R::Ray{T},t::Real) where T <: AnyComplex
	τ = exp(complex(0,R.angle))
	R.reverse ? T(-τ/t^2) : T(τ/(1-t)^2)
end

+(R::Ray,z::Number) = Ray(R.base+z,R.angle,R.reverse)
-(R::Ray) = Ray(-R.base,mod2pi(R.angle+pi),R.reverse)
*(R::Ray,z::Number) = Ray(z*R.base,mod2pi(R.angle+sign(z)),R.reverse)

"""
	inv(R)
Invert the ray `R` through the origin. In general the inverse is an `Arc`.
"""
function inv(R::Ray)
	w = 1 ./ point(R,[0,0.5,1])
	Arc(w...)
end

# Other methods
isfinite(::Ray) = false
conj(R::Ray) = Ray(conj(R.base),-R.angle,R.reverse)
reverse(R::Ray) = Ray(R.base,R.angle,!R.reverse)

"""
	isapprox(R1::Ray,R2::Ray; tol=<default>)
	R1 ≈ R2
Determine if `R1` and `R2` represent the same segment, irrespective of the type or values of its parameters. Identity is determined by agreement within `tol`, which is interpreted as the weaker of absolute and relative differences.
"""
function isapprox(R1::Ray,R2::Ray;tol=DEFAULT[:tol])
	return isapprox(R1.base,R2.base,rtol=tol,atol=tol) && (abs(mod2pi(R1.angle-R2.angle)) < tol)
end

"""
	isleft(z,R::Ray)
Determine whether the number `z` lies "to the left" of ray `R`. This means that the angle it makes with `tangent(R)` is in the interval (0,π).

Note that `isleft` and `isright` are *not* logical opposites; a point on the (extended) ray should give `false` in both cases.
"""
function isleft(z::Number,R::Ray)
	a,b = point(R,[0.2,0.8])  # accounts for reversal
	(real(b)-real(a)) * (imag(z)-imag(a)) > (real(z)-real(a)) * (imag(b)-imag(a))
end

"""
	isright(z,R::Ray)
Determine whether the number `z` lies "to the right" of ray `R`. This means that the angle it makes with `tangent(R)` is in the interval (-π,0).

Note that `isleft` and `isright` are *not* logical opposites; a point on the (extended) ray should give `false` in both cases.
"""
function isright(z::Number,R::Ray)
	a,b = point(R,[0.2,0.8])  # accounts for reversal
	(real(b)-real(a)) * (imag(z)-imag(a)) < (real(z)-real(a)) * (imag(b)-imag(a))
end

"""
	dist(z,R::Ray)
Compute the distance from number `z` to the ray `R`.
"""
dist(z::Number,R::Ray) = abs(z - closest(z,R))

"""
	closest(z,R::Ray)

Find the point on ray `R` that lies closest to `z`.
"""
function closest(z::Number,R::Ray)
	# translate and rotate to positive Re axis
	s = exp(complex(0,R.angle))
	ζ = (z-R.base)/s
	R.base + max(real(ζ),0)*s
end

sign(R::Ray) = unittangent(R)

# Display methods
function show(io::IO,R::Ray{T}) where {T}
	print(IOContext(io,:compact=>true),"Ray(",R.base,",",R.angle,",",R.reverse,")")
end

function show(io::IO,::MIME"text/plain",R::Ray{T}) where {T}
	if R.reverse
		print(io,"Ray from ∞ to ",R.base," at angle ",R.angle)
	else
		print(io,"Ray from ",R.base," to ∞ at angle ",R.angle)
	end
end

# in the plane, use two (finite) points for plotting
plotdata(R::Ray{T}) where T<:Union{Complex,Polar} = R.reverse ? [R(0.3),R.base] : [R.base,R(0.7)]


text:
export SampledInteraction, IntegratedInteraction
export sampled_interaction, integrated_interaction, compute_interaction
export SI, II

@doc raw"""
    sampled_interaction([t,] x, W′, ys)

Computes ``-(W' * \dot\rho)(t, x)``, which is the sampled approximation of ``-(W' * \rho)(t, x)``,
where ``\rho`` is the piecewise-constant density associated to the particles `ys`.

It `t` is omitted, then `W′(x)` is assumed independent of time.

See also [`integrated_interaction`](@ref), [`compute_interaction`](@ref).
"""
function sampled_interaction end

function sampled_interaction(x::Real, Wprime, ys::AbstractVector{<:Real})
    -sum(Wprime(x - y) for y in ys) / (length(ys) - 1)
end

function sampled_interaction(t::Real, x::Real, Wprime, ys::AbstractVector{<:Real})
    -sum(Wprime(t, x - y) for y in ys) / (length(ys) - 1)
end

function sampled_interaction(x::Real; Wprime, particles::AbstractVector{<:Real})
    -sum(Wprime(x - p) for p in particles) / (length(particles) - 1)
end

function sampled_interaction(t::Real, x::Real; Wprime, particles::AbstractVector{<:Real})
    -sum(Wprime(t, x - p) for p in particles) / (length(particles) - 1)
end

@doc raw"""
    integrated_interaction([t,] x, W, ys[, dens_diff])

Computes ``-(W' * \rho)(t, x) = -(W * \rho')(t,x)``, where ``\rho`` is the piecewise-constant
density associated to the particles `ys`.

It `t` is omitted, then `W(x)` is assumed independent of time.

!!! note

    To ensure the correctness of the computation, `dens_diff` must coincide with `diff(pwc_density(ys))`.
    It can be pre-computed and passed explicitly to allow reuse (as an optimization).

See also [`sampled_interaction`](@ref), [`compute_interaction`](@ref).
"""
function integrated_interaction end

function integrated_interaction(x::Real, W, ys::AbstractVector{<:Real}, dens_diff::AbstractVector{<:Real}=diff(pwc_density(ys)))
    -sum(i -> dens_diff[i] * W(x - ys[i]), eachindex(ys))
end

function integrated_interaction(t::Real, x::Real, W, ys::AbstractVector{<:Real}, dens_diff::AbstractVector{<:Real}=diff(pwc_density(ys)))
    -sum(i -> dens_diff[i] * W(t, x - ys[i]), eachindex(ys))
end

"""
    SampledInteraction(W′)

Represents the sampled interaction induced by the function `W′`.

The interaction can be evaluated with [`compute_interaction`](@ref).

See also [`IntegratedInteraction`](@ref).
"""
struct SampledInteraction{TWprime}
    Wprime::TWprime
end

"""
    IntegratedInteraction(W)

Represents the integrated interaction induced by the function `W`.

The interaction can be evaluated with [`compute_interaction`](@ref).

See also [`SampledInteraction`](@ref).
"""
struct IntegratedInteraction{TW}
    W::TW
end

"""
    compute_interaction(t, x, interaction, ys[, dens_diff])

Evaluates the interaction generated by the particles `ys` at `(t, x)`.

`interaction` can be a `SampledInteraction`, an `IntegratedInteraction`, or any other
form of interaction that provides this interface.

See also [`SampledInteraction`](@ref), [`IntegratedInteraction`](@ref),
[`sampled_interaction`](@ref), [`integrated_interaction`](@ref).
"""
function compute_interaction end

function compute_interaction(t::Real, x::Real, int::SampledInteraction, ys::AbstractVector{<:Real}, dens_diff::AbstractVector{<:Real}=Float16[])
    sampled_interaction(t, x, int.Wprime, ys)
end

function compute_interaction(t::Real, x::Real, int::IntegratedInteraction, ys::AbstractVector{<:Real}, dens_diff::AbstractVector{<:Real}=diff(pwc_density(ys)))
    integrated_interaction(t, x, int.W, ys, dens_diff)
end

const SI = SampledInteraction
const II = IntegratedInteraction


text:
using ModelingToolkit
using Test

# `Expr`, `Number` -> `Operation`
@Param a
@Param b
@Unknown x(t) y()
@test convert(Expression, 2) == 2
expr = :(-inv(2sqrt(+($a, $b))))
op   = Operation(-, [Operation(inv,
                [Operation(*, [2, Operation(sqrt,
                              [Operation(+, [a, b])])])])])
@test convert(Expression, expr) == op
expr1 = :($x^($y-1))
op1   = Operation(^, [x, Operation(-, [y, 1])])
@test convert(Expression, expr1) == op1

@test_throws ArgumentError convert(Expression, :([a, b]))
@test_throws ArgumentError convert(Expression, :(a ? b : c))


text:
include("C:/Users/stravaglino3/Downloads/julia_code/julia_code/model_fit_global.jl")
include("C:/Users/stravaglino3/Downloads/julia_code/models.jl")
include("C:/Users/stravaglino3/Downloads/julia_code/helper_functions.jl")

f = (du,u,p,t) -> bimolecular_diss!(du,u,p,t,dens,cons) # enclose constants

# u0 = [0.0;0.0]
u₀ = [1.0]          # Initial condition
tspan = (0.0,10.0)  # Simulation interval
p = 0.5      # equation parameter. p = [k₁, k₋₁]
dens = []
cons = []

# Setup the ODE problem, then solve
prob = ODEProblem(f,u₀,tspan,p)
sol = solve(prob,Tsit5())
# sol = solve(prob,Tsit5())
plot(sol)
dataset1 = generate_data(sol)
# plot(sol.t,log.(Array(sol))')
#
# fs = collect(0:20)
# ks = bell_diss.(0.5,10,fs)
# plot(fs,ks)
model_opt2 = (1.0e-6,10.0)
pinit = 0.0
t = collect(range(0,stop=10,length=101))
rates,hes,l = optimization(bimolecular_diss!,[vec(dataset1)],
                            [[]],[[]],model_opt2,[t])

Juno.@enter optimization(bimolecular_diss!,[vec(dataset1)],
            [[]],[[]],model_opt2,[t])


text:
################################################################################
# Comparison of Zig-Zag for diffusion bridges with tailored Poisson rates
# and with adapted Poisson rate. Reference: https://arxiv.org/abs/2001.05889.
################################################################################

using ZigZagBoomerang, SparseArrays, LinearAlgebra
#using CairoMakie
include("faberschauder.jl")
const ZZB = ZigZagBoomerang

# Drift
const α = 1.5
const L = 7
const T = 50.0
b(x) = α * sin(x)
# First derivative
b′(x) = α * cos(x)
# Second derivative
b″(x) = -α * sin(x)

n = (2 << L) + 1
ξ0 = 0randn(n)
u, v = -π, 3π  # initial and fianl point
ξ0[1] = u / sqrt(T)
ξ0[end] = v / sqrt(T)
c = ones(n)
c[end] = c[1] = 0.0
θ0 = rand((-1.0, 1.0), n)
θ0[end] = θ0[1] = 0.0 # fix final point
T′ = 30000.0 # final clock of the pdmp

Γ = sparse(1.0I, n, n)
#trace, (t, ξ, θ), (acc, num) = @time pdmp(∇ϕ!, 0.0, ξ0, θ0, T, 10.0, Boomerang(Γ, ξ0*0, 0.1; ρ=0.9), 1, L, adapt=false);
#trace, (t, ξ, θ), (acc, num) = @time pdmp(∇ϕ, 0.0, ξ0, rand((-1.0, 1.0), n), T, 40.0*ones(n), ZigZag(Γ, ξ0*0), 5, L, adapt=false);
trace, (t, ξ, θ), (acc, num), c = @time spdmp(∇ϕmoving, 0.0, ξ0, θ0, T′, c, ZigZag(Γ, ξ0 * 0),
                        SelfMoving(), L, T, adapt = true);

@show acc/num
################################################################################
# Overloafing Poisson times in order to have tighter upperbounds
################################################################################
struct MyBound
    c::Float64
end
function ZZB.adapt!(b::Vector{MyBound}, i, x)
    b[i] = MyBound(b[i].c * x)
end

"""
    poisson_time(a, b, c, u)
Obtaining waiting time for inhomogeneous Poisson Process
with rate of the form λ(t) = a + (b + c*t)^+, where `c`,`a`> 0 ,`b` ∈ R, `u` uniform random variable
"""
function ZZB.poisson_time((a, b, c)::NTuple{3}, u = rand()) # formula (22)
    if b > 0
        return (-(b + a) + sqrt((b + a)^2 - 2.0 * c * log(u))) / c # positive solution of quadratic equation c*0.5 x^2 + (b + a) x + log(u) = 0
    elseif a * b / c <= log(u)
        return -log(u) / a
    else
        return (-(a + b) + sqrt((a + b)^2 - 2.0 * c * (b * b * 0.5 / c + log(u)))) / c    # positive solution of quadratic equation c*0.5 x^2 + (b + a) x + log(u) + b*b*0.5/c = 0
    end
end

ZZB.sλ̄((a, b, c)::NTuple{3}, Δt) = a + ZZB.pos(b + c * Δt)
"""
    abc(G, i, x, θ, c, Flow)

Returns the constant term `a` and linear term `b` when computing the Poisson times
from the upper upper bounding rates λᵢ(t) = max(a + b*t)^2. The factors `a` and `b`
can be function of the current position `x`, velocity `θ`, tuning parameter `c` and
the Graph `G`
"""
function ZZB.ab(G, i, x, θ, c::Vector{MyBound}, F::ZigZag)
        if i == 1
            a = c[i].c + T^(1.5)*0.5*(α^2 + α) * abs(θ[i])  # initial point
            b1 = θ[i]*(x[i] - x[end])
            b2  θ[i]*(θ[i] - θ[end])
        elseif i == (2 << L) + 1
            a = c[i].c + T^(1.5)*0.5*(α^2 + α) * abs(θ[i])  # final point
            b1 = θ[i]*(x[i] - x[1])
            b2  θ[i]*(θ[i] - θ[1])
        else
            l = lvl(i, L)
            a = c[i].c + T^(1.5) / 2^((L - l) * 1.5 + 2) * (α^2 + α) * abs(θ[i]) # formula (22)
            b1 = x[i] * θ[i]
            b2 = θ[i] * θ[i]
        end
        return a, b1, b2
    end
end

c = [MyBound(0.0) for i in 1:n]
trace, (t, ξ, θ), (acc, num), c = @time spdmp(∇ϕmoving, 0.0, ξ0, θ0, T′, c, ZigZag(Γ, ξ0 * 0),
                        SelfMoving(), L, T, adapt = false);

@show acc/num
ts, ξs = splitpairs(discretize(trace, T′/n))
S = T*(0:n)/(n+1)

#using CairoMakie
p1 = lines(S, [dotψ(ξ, s, L, T) for s in S], linewidth=0.3)
for ξ in ξs[1:5:end]
    lines!(p1, S, [dotψ(ξ, s, L, T) for s in S], linewidth=0.3)
end
display(p1)


text:
include("runBatchCompareSim.jl")
include("linspecer.jl")
using Printf

function compareSims(A::Array{Float64,2}, As::Array{Array{Float64,2},1}, C::Array{Float64,2}, muPrior::Array{Float64,1}, endT::Number, dt::Number,
                    sigmas::Array{Float64,1}, simNames::Array{String,1})
    set_zero_subnormals(true)
    variationNum = length(As)
    colors = linspecer(variationNum)
    static = StaticWorld(A, C, muPrior, endT, dt)
    # figure()
    for v = 1:variationNum
        println(string("Variation #", v))
        plotOpts = PlotOpts(simNames[v], colors[v,:])
        simOpts = SimOpts(sigmas, 1000)
        @time runBatchSim(plotOpts, static, simOpts, As)
    end
    legend()

    return
end

function compareSims(As::Array{Array{Float64,2},1}, C::Array{Float64,2}, muPrior::Array{Float64,1}, sigmas::Array{Float64,1}, simNames::Array{String,1})
    set_zero_subnormals(true)
    variationNum = length(As)
    colors = linspecer(variationNum)
    static = StaticWorld(A, C, muPrior)
    # figure()
    for v = 1:variationNum
        println(string("Variation #", v))
        plotOpts = PlotOpts(simNames[v], colors[v,:])
        simOpts = SimOpts(sigmas, 1000)
        @time runBatchSim(plotOpts, static, simOpts, As)
    end
    legend()

    return
end

max_stars_repo_name: rhoadesScholar/Drugowitsch-Rotation

max_issues_repo_name: rhoadesScholar/Drugowitsch-Rotation

max_forks_repo_name: rhoadesScholar/Drugowitsch-Rotation


text:
global vis = nothing

function init_vis()
    global vis
    if vis == nothing
        vis = Visualizer()
        open(vis)
    end
    vis
end

function visualize(m::AbstractMesh...)

    # first visual
    vis = init_vis()
    delete!(vis) # clear visualizer
    for mesh in m
        setobject!(vis, mesh)
    end
end

visualize(m::AbstractVector) = visualize(m...)


text:
function close_eqs(sys::MomentEquations, closure_exp::OrderedDict,
                   closure::OrderedDict, polynorm::Bool)

    # TODO: improve performance

    closed_eqs = Equation[]
    for eq in sys.odes.eqs

        closed_rhs = substitute(eq.rhs, closure_exp)
        # apply binomial expansion on the expressions
        if polynorm # depending on the functional form
            closed_rhs = expand(closed_rhs)
        else
            closed_rhs = expand_expr(closed_rhs)
        end
        push!(closed_eqs, Equation(eq.lhs, closed_rhs))
    end

    iv = get_iv(sys.odes)
    ps = get_ps(sys.odes)

    vars = sys.odes.states[1:(length(sys.iter_1)+length(sys.iter_m))]

    odename = Symbol(nameof(sys), "_closed")
    odes = ODESystem(closed_eqs, iv, vars, ps; name=odename)

    ClosedMomentEquations(odes, closure, sys)

end

"""
    moment_closure(sys::MomentEquations, closure::String, binary_vars::Array{Int,1}=Int[])

Given `MomentEquations`, apply the specified moment closure approximation and return
the [`ClosedMomentEquations`](@ref).

The supported `closure` options are:
- [`"zero"`](@ref zero_closure)
- [`"normal"`](@ref normal_closure)
- [`"log-normal"`](@ref log-normal_closure)
- [`"poisson"`](@ref poisson_closure)
- [`"gamma"`](@ref gamma_closure)
- [`"derivative matching"`](@ref derivative_matching)
- [`"conditional gaussian"`](@ref conditional_gaussian_closure)
- [`"conditional derivative matching"`](@ref conditional_derivative_matching)

# Notes
- `binary_vars` *must* be specified for conditional closures as an array of indices of all species
  (as in [`speciesmap`](@ref)) which molecule number is a Bernoulli variable. Although not necessary
  for other closures, specifying `binary_vars` is recommended as the properties of Bernoulli variables
  will be used to remove the redundant moment equations and simplify the expressions, which can
  significantly improve numerical stability.
"""
function moment_closure(sys::MomentEquations, closure::String, binary_vars::Array{Int,1}=Int[])

    if closure == "zero"
        zero_closure(sys, binary_vars)
    elseif closure == "normal"
        normal_closure(sys, binary_vars)
    elseif closure == "log-normal"
        log_normal_closure(sys, binary_vars)
    elseif closure == "poisson"
        poisson_closure(sys, binary_vars)
    elseif closure == "gamma"
        gamma_closure(sys, binary_vars)
    elseif closure == "derivative matching"
        derivative_matching(sys, binary_vars)
    elseif closure == "conditional gaussian"
        conditional_gaussian_closure(sys, binary_vars)
    elseif closure == "conditional derivative matching"
        conditional_derivative_matching(sys, binary_vars)
    else
        error(closure*" closure does not exist")
    end

end