text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    @stage(def)

Add a stage model generation recipe to `stochasticprogram` using the syntax
```julia
@stage stage stochasticprogram::StochasticProgram = begin
    @parameters param1 param2 ...
    @decision(stochasticprogram, var) ...
    @uncertain ξ
    ... JuMPdef ...
    ...
end
```
where `stage` is the stage number and JuMP syntax is used inside the block to define the stage model. During definition, the stage model is referenced using the same variable name as `stochasticprogram`.

## Examples

The following defines the first stage model given by:
```math
  minimize 100x₁ + 150x₂
    s.t  x₁ + x₂ ≤ 120
         x₁ ≥ 40
         x₂ ≥ 20
```
and the second-stage model given by:
```math
  maximize q₁(ξ)y₁ + q₂(ξ)y₂
    s.t  6y₁ + 10y₂ ≤ 60x₁
         8y₁ + 5y₂ ≤ 60x₂
         0 ≤ y₁ ≤ d₁(ξ)
         0 ≤ y₂ ≤ d₂(ξ)
```
where ``q₁(ξ), q₂(ξ), d₁(ξ), d₂(ξ)`` depend on the scenario ``ξ`` and ``x₁, x₂`` are first stage variables. Two scenarios are added so that two second stage models are generated.

```jldoctest
ξ₁ = @scenario q₁ = 24.0 q₂ = 28.0 d₁ = 500.0 d₂ = 100.0 probability = 0.4
ξ₂ = @scenario q₁ = 28.0 q₂ = 32.0 d₁ = 300.0 d₂ = 300.0 probability = 0.6

sp = StochasticProgram([ξ₁, ξ₂])

@stage 1 sp = begin
    @decision(sp, x₁ >= 40)
    @decision(sp, x₂ >= 20)
    @objective(sp, Min, 100*x₁ + 150*x₂)
    @constraint(sp, x₁ + x₂ <= 120)
end

@stage 2 sp = begin
    @uncertain q₁ q₂ d₁ d₂
    @variable(sp, 0 <= y₁ <= d₁)
    @variable(sp, 0 <= y₂ <= d₂)
    @objective(sp, Max, q₁*y₁ + q₂*y₂)
    @constraint(sp, 6*y₁ + 10*y₂ <= 60*x₁)
    @constraint(sp, 8*y₁ + 5*y₂ <= 80*x₂)
end

# output

Stochastic program with:
 * 2 decision variables
 * 2 scenarios of type Scenario
Solver is default solver

```

See also: [`@parameters`](@ref), [`@decision`](@ref), [`@uncertain`](@ref)
"""
macro stage(stage, args)
    _error(x, str...) = begin
        macroname = Symbol(String(x.args[1])[2:end])
        JuMP._macro_error(macroname, prettify.(x.args[3:end]), x.args[2], str...)
    end
    @capture(args, sp_Symbol = def_) || _error(args, "Invalid syntax. Expected @stage stage sp = begin ... end")
    model_name = sp
    # Flags for error checking
    seen_decision = :(false)
    seen_recourse = :(false)
    # Decision definitions might require parameter calculations,
    # so we first need to extract and save any such lines
    vardefs = Expr(:block)
    lastline = __source__
    decisiondefs = postwalk(def) do x
        if x isa LineNumberNode
            lastline = LineNumberNode(x.line, x.file)
        end
        if  @capture(x, @variable(m_Symbol, variabledef__)) ||
            @capture(x, @known(m_Symbol, variabledef__)) ||
            @capture(x, @constraint(m_Symbol, constdef__)) ||
            @capture(x, @expression(m_Symbol, expressiondef__)) ||
            @capture(x, @objective(m_Symbol, objdef__))
            # Check model name
            if m != model_name
                _error(x, "Inconsistent model name \"$m\", should be \"$model_name\".")
            end
            # Skip any line related to the JuMP model
            return Expr(:block)
        elseif @capture(x, @parameters args__) ||
               @capture(x, @uncertain args__)
            # Skip any line related to stochastics or unhandled @known/@parameter lines
            return Expr(:block)
        elseif @capture(x, @decision m_Symbol args__) || @capture(x, @recourse m_Symbol args__)
            # Check model name
            if m != model_name
                _error(x, "Inconsistent model name \"$m\", should be \"$model_name\".")
            end
            # Handle @decision
            return @q begin
                StochasticPrograms.@_decision known $stage $m $((args)...)
            end
        else
            # Anything else could be required for decision variable construction, and is therefore saved
            return x
        end
    end
    # Next, handle @decision annotations in main definition
    def = postwalk(def) do x
        if @capture(x, @decision args__)
            # Bookkeep for error checking
            seen_decision = :(true)
            return @q begin
                StochasticPrograms.@_decision decision $stage $((args)...)
            end
        end
        if @capture(x, @recourse args__)
            # Bookkeep for error checking
            seen_recourse = :(true)
            return @q begin
                StochasticPrograms.@_decision recourse $stage $((args)...)
            end
        end
        return x
    end
    # Handle parameters
    def = postwalk(def) do x
        if @capture(x, @parameters arg_)
            code = Expr(:block)
            for paramdef in block(prettify(arg)).args
                if @capture(paramdef, key_Symbol = val_) || @capture(paramdef, key_Symbol)
                    push!(code.args, :($key = stage.$key))
                else
                    _error(x, "Incorrect @parameters syntax. Specify parameter names, possibly with a default value.")
                end
            end
            # Extracted parameters might be required for decision variable construction
            pushfirst!(decisiondefs.args, code)
            return code
        elseif @capture(x, @parameters args__)
            code = Expr(:block)
            for paramdef in args
                if @capture(paramdef, key_Symbol = val_) || @capture(paramdef, key_Symbol)
                    push!(code.args, :($key = stage.$key))
                else
                    _error(x, "Incorrect @parameters syntax. Specify parameter names, possibly with a default value.")
                end
            end
            # Extracted parameters might be required for decision variable construction
            pushfirst!(decisiondefs.args, code)
            return code
         elseif @capture(x, @uncertain var_Symbol::t_Symbol)
            stage == 1 && _error(x, "@uncertain declarations cannot be used in the first stage.")
            return :($var::$t = scenario)
        elseif @capture(x, @uncertain vars__ from scenvar_Symbol::t_Symbol)
            stage == 1 && _error(x, "@uncertain declarations cannot be used in the first stage.")
            code = @q begin
                $scenvar::$t = scenario
            end
            for var in vars
                varkey = Meta.quot(var)
                push!(code.args, :($var = getproperty($scenvar, $varkey)))
            end
            return code
        elseif @capture(x, @uncertain vars__ from t_Symbol)
            stage == 1 && _error(x, "@uncertain declarations cannot be used in the first stage.")
            code = @q begin
                scenvar::$t = scenario
            end
            for var in vars
                varkey = Meta.quot(var)
                push!(code.args, :($var = getproperty(scenvar, $varkey)))
            end
            return code
        elseif @capture(x, @uncertain var_Symbol::t_Symbol = scenariodef_)
            stage == 1 && _error(x, "@uncertain declarations cannot be used in the first stage.")
            return @q begin
                $var::$t = scenario
            end
        elseif @capture(x, @uncertain vars__ from scenvar_Symbol::t_Symbol = scenariodef_)
            stage == 1 && _error(x, "@uncertain declarations cannot be used in the first stage.")
            code = @q begin
                $scenvar::$t = scenario
            end
            for var in vars
                varkey = Meta.quot(var)
                push!(code.args, :($var = getproperty($scenvar, $varkey)))
            end
            return code
        elseif @capture(x, @uncertain vars__ from t_Symbol = scenariodef_)
            stage == 1 && _error(x, "@uncertain declarations cannot be used in the first stage.")
            code = @q begin
                scenvar::$t = scenario
            end
            for var in vars
                varkey = Meta.quot(var)
                push!(code.args, :($var = getproperty(scenvar, $varkey)))
            end
            return code
        elseif @capture(x, @uncertain args_) && typeof(args) == Expr && args.head == :ref
            stage == 1 && _error(x, "@uncertain declarations cannot be used in the first stage.")
            scenvar = first(args.args)
            reference = :(Containers.@container $args "uncertain value")
            code = @q begin
                typeof(scenario) <: Scenario{D} where D <: Union{Array, StochasticPrograms.DenseAxisArray} || $_error($(Meta.quot(x)), "@uncertain declarations of type `@uncertain ξ[i = ..., j = ..., ...]` only support scenarios of type `Scenario` with Array or DenseAxisArray as underlying data, given scenario is: \n $scenario.\n Use a matching `@scenario` call to generate scenarios or consider declaring a custom scenario type.")
                size(scenario.data) == size($reference) || $_error($(Meta.quot(x)), "Given scenario \n\n$scenario \n\ndoes not match @uncertain declaration \n\n$($reference).")
                $scenvar = scenario.data
            end
            return code
        elseif @capture(x, @uncertain args_) && typeof(args) == Expr && args.head == :typed_vcat
            stage == 1 && _error(x, "@uncertain declarations cannot be used in the first stage.")
            scenvar = first(args.args)
            reference = :(Containers.@container $args "uncertain value")
            code = @q begin
                typeof(scenario) <: Scenario{D} where D <: StochasticPrograms.SparseAxisArray || $_error($(Meta.quot(x)), "@uncertain declarations of type `@uncertain ξ[i = ..., j = ..., ...; ...]` only support scenarios of type `Scenario` with SparseAxisArray as underlying data, given scenario is: \n $scenario.\n Use a matching `@scenario` call to generate scenarios or consider declaring a custom scenario type.")
                keys(scenario.data.data) == keys($reference.data) || $_error($(Meta.quot(x)), "Given scenario \n\n$scenario \n\ndoes not match @uncertain declaration \n\n$($reference).")
                $scenvar = scenario.data
            end
            return code
        elseif @capture(x, @uncertain args__)
            stage == 1 && _error(x, "@uncertain declarations cannot be used in the first stage.")
            code = @q begin
                typeof(scenario) <: Scenario || $_error($(Meta.quot(x)), "@uncertain declarations of type `@uncertain var1, var2, ...` only support scenarios of type `Scenario`, given scenario is:\n$scenario.\n Use a matching `@scenario` call to generate scenarios or consider declaring a custom scenario type.")
            end
            for var in args
                varkey = Meta.quot(var)
                code = @q begin
                    $code
                    $var = try
                        $var = scenario.data[$varkey]
                    catch err
                        $_error($(Meta.quot(x)), "Given scenario $scenario does not match @uncertain declaration.")
                    end
                end
            end
            return code
        else
            return x
        end
    end
    # The first stage is treated different from the subsequent stages
    generatordefs = if stage == 1
        @q begin
            $(esc(sp)).generator[:stage_1] = ($(esc(model_name))::JuMP.Model, $(esc(:stage))) -> begin
                $(esc(def))
                # Cache sense and objective function
                sense = objective_sense($(esc(model_name)))
                obj = moi_function(objective_function($(esc(model_name))))
                add_stage_objective!($(esc(model_name)), 1, sense, obj)
	            return $(esc(model_name))
            end
        end
    else
        @q begin
            $(esc(sp)).generator[Symbol(:stage_,$stage)] = ($(esc(model_name))::JuMP.Model, $(esc(:stage)), $(esc(:scenario))::AbstractScenario) -> begin
                $(esc(def))
                # Cache sense and objective function
                sense = objective_sense($(esc(model_name)))
                obj = moi_function(objective_function($(esc(model_name))))
                add_stage_objective!($(esc(model_name)), $stage, sense, obj)
	            return $(esc(model_name))
            end
        end
    end
    # Create definition code
    code = @q begin
        isa($(esc(sp)), StochasticProgram) || error("Given object is not a stochastic program.")
        n = num_stages($(esc(sp)))
        $stage > n && error("Cannot specify stage $($stage) for stochastic program with $n stages.")
        $stage == n && $seen_decision && error("@decision declarations cannot be used in the final stage. Consider @recourse instead.")
        if $stage < n
            $seen_recourse && error("@recourse declarations can only be used in the final stage.")
            $(esc(sp)).generator[Symbol(:stage_,$stage,:_decisions)] = ($(esc(model_name))::JuMP.Model, $(esc(:stage))) -> begin
                $(esc(decisiondefs))
	            return $(esc(model_name))
            end
        end
        # Stage model generation code
        $generatordefs
        # Return stochastic program
        $(esc(sp))
    end
    # Return code
    return code
end
"""
    @first_stage(def)

Add a first stage model generation recipe to `stochasticprogram` using the syntax
```julia
@first_stage stochasticprogram::StochasticProgram = begin
    ...
end [defer]
```
where JuMP syntax is used inside the block to define the first stage model. During definition, the first stage model is referenced through the reserved keyword matching the name of `stochasticprogram`.

## Examples

The following defines the first stage model given by:

```math
  minimize 100x₁ + 150x₂
    s.t  x₁ + x₂ ≤ 120
         x₁ ≥ 40
         x₂ ≥ 20
```

```julia
@first_stage sp = begin
    @decision(sp, x₁ >= 40)
    @decision(sp, x₂ >= 20)
    @objective(sp, Min, 100*x₁ + 150*x₂)
    @constraint(sp, x₁ + x₂ <= 120)
end
```

See also: [`@second_stage`](@ref)
"""
macro first_stage(arg)
    @capture(arg, sp_Symbol = def_) || error("Invalid syntax. Expected @first_stage sp = begin ... end")
    return esc(:(@stage 1 $arg))
end
"""
    @second_stage(def)

Add a second stage model generation recipe to `stochasticprogram` using the syntax
```julia
@second_stage stochasticprogram::StochasticProgram = begin
    @known var1 var2 ...
    ...
end
```
where JuMP syntax is used inside the block to define the second stage model. During definition, the second stage model is referenced through the reserved keyword matching the name of `stochasticprogram`.

## Examples

The following defines the second stage model given by:
```math
  minimize q₁(ξ)y₁ + q₂(ξ)y₂
    s.t  6y₁ + 10y₂ ≤ 60x₁
         8y₁ + 5y₂ ≤ 60x₂
         0 ≤ y₁ ≤ d₁(ξ)
         0 ≤ y₂ ≤ d₂(ξ)
```
where ``q₁(ξ), q₂(ξ), d₁(ξ), d₂(ξ)`` depend on the scenario ``ξ`` and ``x₁, x₂`` are first stage variables. Two scenarios are added so that two second stage models are generated.

```julia
@second_stage sp = begin
    @known(sp, x₁, x₂)
    @uncertain q₁ q₂ d₁ d₂
    @variable(sp, 0 <= y₁ <= d₁)
    @variable(sp, 0 <= y₂ <= d₂)
    @objective(sp, Min, q₁*y₁ + q₂*y₂)
    @constraint(sp, 6*y₁ + 10*y₂ <= 60*x₁)
    @constraint(sp, 8*y₁ + 5*y₂ <= 80*x₂)
end
```

See also: [`@first_stage`](@ref)
"""
macro second_stage(arg)
    @capture(arg, sp_Symbol = def_) || error("Invalid syntax. Expected @second_stage sp = begin ... end")
    return esc(@q begin
        @stage 2 $arg
        generate!($sp)
    end)
end
"""
    @parameters(def)

Define the problem parameters in a @stage block
```julia
@parameters param1, param2, ...
```
possibly with default values. Any defined parameter without a default value must be supplied as a keyword argument to [`instantiate`](@ref) when creating models.

## Examples

```julia
@parameters d

@parameters begin
    Crops = [:wheat, :corn, :beets]
    Cost = Dict(:wheat=>150, :corn=>230, :beets=>260)
    Budget = 500
end
```

See also [`@decision`](@ref), [`@uncertain`](@ref), [`@stage`](@ref)
"""
macro parameters(def) @warn "@parameters should be used inside a @stage block." end


text:
using .CUDA
using .CUDA: CuArray, @linearidx, GPUArrays, CUBLAS
using LinearAlgebra

export togpu

function togpu(tn::TensorNetwork)
    TensorNetwork(togpu.(tn.tensors))
end

function togpu(t::LabeledTensor)
    LabeledTensor(CuArray(t.array), t.labels, t.meta)
end

function genperm(I::NTuple{N}, perm::NTuple{N}) where N
    ntuple(d-> (@inbounds return I[perm[d]]), Val(N))
end

function LinearAlgebra.permutedims!(dest::GPUArrays.AbstractGPUArray, src::GPUArrays.AbstractGPUArray, perm)
    LinearAlgebra.permutedims!(dest::GPUArrays.AbstractGPUArray, src::GPUArrays.AbstractGPUArray, Tuple(perm))
end
function LinearAlgebra.permutedims!(dest::GPUArrays.AbstractGPUArray, src::GPUArrays.AbstractGPUArray, perm::NTuple)
    perm isa Tuple || (perm = Tuple(perm))
    size_dest = size(dest)
    size_src = size(src)
    CUDA.gpu_call(vec(dest), vec(src), perm; name="permutedims!") do ctx, dest, src, perm
        i = @linearidx src
        I = l2c(size_src, i)
        @inbounds dest[c2l(size_dest, genperm(I, perm))] = src[i]
        return
    end
    return reshape(dest, size(dest))
end

function LinearAlgebra.rmul!(a::StridedCuArray{<:CUBLAS.CublasFloat}, k::Number)
    vec(a) .*= k
    return a
end

ext: jl


text:
function calc_tderivative!(integrator, cache, dtd1, repeat_step)
  @inbounds begin
    @unpack t,dt,uprev,u,f,p = integrator
    @unpack du2,fsalfirst,dT,tf,linsolve_tmp = cache

    # Time derivative
    if !repeat_step # skip calculation if step is repeated
      if DiffEqBase.has_tgrad(f)
        f.tgrad(dT, uprev, p, t)
      else
        tf.uprev = uprev
        tf.p = p
        derivative!(dT, tf, t, du2, integrator, cache.grad_config)
      end
    end

    f(fsalfirst, uprev, p, t)
    integrator.destats.nf += 1
    @.. linsolve_tmp = fsalfirst + dtd1*dT
  end
end

function calc_tderivative(integrator, cache)
  @unpack t,dt,uprev,u,f,p = integrator

  # Time derivative
  if DiffEqBase.has_tgrad(f)
    dT = f.tgrad(uprev, p, t)
  else
    tf = cache.tf
    tf.u = uprev
    tf.p = p
    dT = derivative(tf, t, integrator)
  end
  dT
end

"""
    calc_J(integrator, cache)

Return a new Jacobian object.

If `integrator.f` has a custom Jacobian update function, then it will be called. Otherwise,
either automatic or finite differencing will be used depending on the `uf` object of the
cache.
"""
function calc_J(integrator, cache)
  @unpack t,uprev,f,p,alg = integrator

  if alg isa DAEAlgorithm
    if DiffEqBase.has_jac(f)
      J = f.jac(duprev, uprev, p, t)
    else
      @unpack uf = cache
      x = zero(uprev)
      J = jacobian(uf, x, integrator)
    end
  else
    if DiffEqBase.has_jac(f)
      J = f.jac(uprev, p, t)
    else
      @unpack uf = cache

      uf.f = nlsolve_f(f, alg)
      uf.p = p
      uf.t = t

      J = jacobian(uf, uprev, integrator)
    end

    integrator.destats.njacs += 1

    if alg isa CompositeAlgorithm
      integrator.eigen_est = opnorm(J, Inf)
    end
  end

  J
end

"""
    calc_J!(J, integrator, cache) -> J

Update the Jacobian object `J`.

If `integrator.f` has a custom Jacobian update function, then it will be called. Otherwise,
either automatic or finite differencing will be used depending on the `cache`.
"""
function calc_J!(J, integrator, cache)
  @unpack t,uprev,f,p,alg = integrator

  if alg isa DAEAlgorithm
    if DiffEqBase.has_jac(f)
      f.jac(J, duprev, uprev, p, t)
    else
      @unpack du1, uf, jac_config = cache
      # using `dz` as temporary array
      x = cache.dz
      fill!(x, zero(eltype(x)))
      jacobian!(J, uf, x, du1, integrator, jac_config)
    end
  else
    if DiffEqBase.has_jac(f)
      f.jac(J, uprev, p, t)
    else
      @unpack du1, uf, jac_config = cache

      uf.f = nlsolve_f(f, alg)
      uf.t = t
      uf.p = p

      jacobian!(J, uf, uprev, du1, integrator, jac_config)
    end

    integrator.destats.njacs += 1

    if alg isa CompositeAlgorithm
      integrator.eigen_est = opnorm(J, Inf)
    end
  end

  return nothing
end

"""
    WOperator(mass_matrix,gamma,J[;transform=false])

A linear operator that represents the W matrix of an ODEProblem, defined as

```math
W = MM - \\gamma J
```

or, if `transform=true`:

```math
W = \\frac{1}{\\gamma}MM - J
```

where `MM` is the mass matrix (a regular `AbstractMatrix` or a `UniformScaling`),
`γ` is a real number proportional to the time step, and `J` is the Jacobian
operator (must be a `AbstractDiffEqLinearOperator`). A `WOperator` can also be
constructed using a `*DEFunction` directly as

    WOperator(f,gamma[;transform=false])

`f` needs to have a jacobian and `jac_prototype`, but the prototype does not need
to be a diffeq operator --- it will automatically be converted to one.

`WOperator` supports lazy `*` and `mul!` operations, the latter utilizing an
internal cache (can be specified in the constructor; default to regular `Vector`).
It supports all of `AbstractDiffEqLinearOperator`'s interface.
"""
mutable struct WOperator{T,
  MType <: Union{UniformScaling,AbstractMatrix},
  GType <: Real,
  JType <: DiffEqBase.AbstractDiffEqLinearOperator{T}
  } <: DiffEqBase.AbstractDiffEqLinearOperator{T}
  mass_matrix::MType
  gamma::GType
  J::JType
  transform::Bool       # true => W = mm/gamma - J; false => W = mm - gamma*J
  inplace::Bool
  _func_cache           # cache used in `mul!`
  _concrete_form         # non-lazy form (matrix/number) of the operator
  WOperator(mass_matrix, gamma, J, inplace; transform=false) = new{eltype(J),typeof(mass_matrix),
    typeof(gamma),typeof(J)}(mass_matrix,gamma,J,transform,inplace,nothing,nothing)
end
function WOperator(f::DiffEqBase.AbstractODEFunction, gamma, inplace; transform=false)
  @assert DiffEqBase.has_jac(f) "f needs to have an associated jacobian"
  if isa(f, Union{SplitFunction, DynamicalODEFunction})
    error("WOperator does not support $(typeof(f)) yet")
  end
  # Convert mass matrix, if needed
  mass_matrix = f.mass_matrix
  if !isa(mass_matrix, Union{AbstractMatrix,UniformScaling})
    mass_matrix = convert(AbstractMatrix, mass_matrix)
  end
  # Convert jacobian, if needed
  J = deepcopy(f.jac_prototype)
  if !isa(J, DiffEqBase.AbstractDiffEqLinearOperator)
    J = DiffEqArrayOperator(J; update_func=f.jac)
  end
  return WOperator(mass_matrix, gamma, J, inplace; transform=transform)
end

set_gamma!(W::WOperator, gamma) = (W.gamma = gamma; W)
DiffEqBase.update_coefficients!(W::WOperator,u,p,t) = (update_coefficients!(W.J,u,p,t); W)
function Base.convert(::Type{AbstractMatrix}, W::WOperator)
  if W._concrete_form === nothing || !W.inplace
    # Allocating
    if W.transform
      W._concrete_form = -W.mass_matrix / W.gamma + convert(AbstractMatrix,W.J)
    else
      W._concrete_form = -W.mass_matrix + W.gamma * convert(AbstractMatrix,W.J)
    end
  else
    # Non-allocating
    _W = W._concrete_form
    J = convert(AbstractMatrix,W.J)
    if W.transform
      if _W isa Diagonal # axpby doesn't specialize on Diagonal matrix
        @inbounds for i in axes(W._concrete_form, 1)
          _W[i, i] = J[i, i] - inv(W.gamma) * W.mass_matrix[i, i]
        end
      else
        copyto!(_W, W.mass_matrix)
        axpby!(one(W.gamma), J, -inv(W.gamma), _W)
      end
    else
      if _W isa Diagonal # axpby doesn't specialize on Diagonal matrix
        @inbounds for i in axes(W._concrete_form, 1)
          _W[i, i] = W.gamma*J[i, i] - W.mass_matrix[i, i]
        end
      else
        copyto!(_W, W.mass_matrix)
        axpby!(W.gamma, J, -one(W.gamma), W._concrete_form)
      end
    end
  end
  W._concrete_form
end
function Base.convert(::Type{Number}, W::WOperator)
  if W.transform
    W._concrete_form = -W.mass_matrix / W.gamma + convert(Number,W.J)
  else
    W._concrete_form = -W.mass_matrix + W.gamma * convert(Number,W.J)
  end
  W._concrete_form
end
Base.size(W::WOperator, args...) = size(W.J, args...)
function Base.getindex(W::WOperator, i::Int)
  if W.transform
    -W.mass_matrix[i] / W.gamma + W.J[i]
  else
    -W.mass_matrix[i] + W.gamma * W.J[i]
  end
end
function Base.getindex(W::WOperator, I::Vararg{Int,N}) where {N}
  if W.transform
    -W.mass_matrix[I...] / W.gamma + W.J[I...]
  else
    -W.mass_matrix[I...] + W.gamma * W.J[I...]
  end
end
function Base.:*(W::WOperator, x::Union{AbstractVecOrMat,Number})
  if W.transform
    (W.mass_matrix*x) / -W.gamma + W.J*x
  else
    -W.mass_matrix*x + W.gamma * (W.J*x)
  end
end
function Base.:\(W::WOperator, x::Union{AbstractVecOrMat,Number})
  if size(W) == () # scalar operator
    convert(Number,W) \ x
  else
    convert(AbstractMatrix,W) \ x
  end
end

function LinearAlgebra.mul!(Y::AbstractVecOrMat, W::WOperator, B::AbstractVecOrMat)
  if W._func_cache === nothing
    # Allocate cache only if needed
    W._func_cache = Vector{eltype(B)}(undef, size(Y, 1))
  end
  if W.transform
    # Compute mass_matrix * B
    if isa(W.mass_matrix, UniformScaling)
      a = -W.mass_matrix.λ / W.gamma
      @.. Y = a * B
    else
      mul!(Y, W.mass_matrix, B)
      lmul!(-1/W.gamma, Y)
    end
    # Compute J * B and add
    mul!(W._func_cache, W.J, B)
    Y .+= W._func_cache
  else
    # Compute mass_matrix * B
    if isa(W.mass_matrix, UniformScaling)
      @.. Y = W.mass_matrix.λ * B
    else
      mul!(Y, W.mass_matrix, B)
    end
    # Compute J * B
    mul!(W._func_cache, W.J, B)
    # Add result
    axpby!(W.gamma, W._func_cache, -one(W.gamma), Y)
  end
end

"""
    islinearfunction(integrator) -> Tuple{Bool,Bool}

return the tuple `(is_linear_wrt_odealg, islinearodefunction)`.
"""
islinearfunction(integrator) = islinearfunction(integrator.f, integrator.alg)

"""
    islinearfunction(f, alg) -> Tuple{Bool,Bool}

return the tuple `(is_linear_wrt_odealg, islinearodefunction)`.
"""
function islinearfunction(f, alg)::Tuple{Bool,Bool}
  isode = f isa ODEFunction && islinear(f.f)
  islin = isode || (alg isa SplitAlgorithms && f isa SplitFunction && islinear(f.f1.f))
  return islin, isode
end

function do_newJW(integrator, alg, nlsolver, repeat_step)::NTuple{2,Bool}
  alg isa DAEAlgorithm && return true, true
  repeat_step && return false, false
  # TODO: RosW
  isnewton(nlsolver) || return true, true
  isfirstcall(nlsolver) && return true, true
  isfs = isfirststage(nlsolver)
  iszero(nlsolver.fast_convergence_cutoff) && return isfs, isfs
  W_iγdt = inv(nlsolver.cache.W_γdt)
  iγdt = inv(nlsolver.γ * integrator.dt)
  smallstepchange = abs(iγdt/W_iγdt - 1) <= get_new_W_γdt_cutoff(nlsolver)
  jbad = nlsolver.status === TryAgain && smallstepchange
  errorfail = integrator.EEst > one(integrator.EEst)
  return jbad, (jbad || (!smallstepchange) || (isfs && errorfail))
end

@noinline _throwWJerror(W, J) = throw(DimensionMismatch("W: $(axes(W)), J: $(axes(J))"))
@noinline _throwWMerror(W, mass_matrix) = throw(DimensionMismatch("W: $(axes(W)), mass matrix: $(axes(mass_matrix))"))
@noinline _throwJMerror(J, mass_matrix) = throw(DimensionMismatch("J: $(axes(J)), mass matrix: $(axes(mass_matrix))"))

@inline function jacobian2W!(W::AbstractMatrix, mass_matrix::MT, dtgamma::Number, J::AbstractMatrix, W_transform::Bool)::Nothing where MT
  # check size and dimension
  iijj = axes(W)
  @boundscheck (iijj === axes(J) && length(iijj) === 2) || _throwWJerror(W, J)
  mass_matrix isa UniformScaling || @boundscheck axes(mass_matrix) === axes(W) || _throwWMerror(W, mass_matrix)
  @inbounds if W_transform
    invdtgamma = inv(dtgamma)
    if MT <: UniformScaling
      copyto!(W, J)
      idxs = diagind(W)
      λ = -mass_matrix.λ
      @.. @view(W[idxs]) = muladd(λ, invdtgamma, @view(J[idxs]))
    else
      @.. W = muladd(-mass_matrix, invdtgamma, J)
    end
  else
    if MT <: UniformScaling
      idxs = diagind(W)
      @.. W = dtgamma*J
      λ = -mass_matrix.λ
      @.. @view(W[idxs]) = @view(W[idxs]) + λ
    else
      @.. W = muladd(dtgamma, J, -mass_matrix)
    end
  end
  return nothing
end

@inline function jacobian2W(mass_matrix::MT, dtgamma::Number, J::AbstractMatrix, W_transform::Bool)::Nothing where MT
  # check size and dimension
  mass_matrix isa UniformScaling || @boundscheck axes(mass_matrix) === axes(J) || _throwJMerror(J, mass_matrix)
  @inbounds if W_transform
    invdtgamma = inv(dtgamma)
    if MT <: UniformScaling
      λ = -mass_matrix.λ
      W = J + (λ * invdtgamma)*I
    else
      W = muladd(-mass_matrix, invdtgamma, J)
    end
  else
    if MT <: UniformScaling
      λ = -mass_matrix.λ
      W = dtgamma*J + λ*I
    else
      W = muladd(dtgamma, J, -mass_matrix)
    end
  end
  return W
end

function calc_W!(W, integrator, nlsolver::Union{Nothing,AbstractNLSolver}, cache, dtgamma, repeat_step, W_transform=false)
  @unpack t,dt,uprev,u,f,p = integrator
  lcache = nlsolver === nothing ? cache : nlsolver.cache
  @unpack J = lcache
  isdae = integrator.alg isa DAEAlgorithm
  alg = unwrap_alg(integrator, true)
  if !isdae
    mass_matrix = integrator.f.mass_matrix
  end
  is_compos = integrator.alg isa CompositeAlgorithm

  # handle Wfact
  if W_transform && DiffEqBase.has_Wfact_t(f)
    f.Wfact_t(W, u, p, dtgamma, t)
    isnewton(nlsolver) && set_W_γdt!(nlsolver, dtgamma)
    is_compos && (integrator.eigen_est = opnorm(LowerTriangular(W), Inf) + inv(dtgamma)) # TODO: better estimate
    return nothing
  elseif !W_transform && DiffEqBase.has_Wfact(f)
    f.Wfact(W, u, p, dtgamma, t)
    isnewton(nlsolver) && set_W_γdt!(nlsolver, dtgamma)
    if is_compos
      opn = opnorm(LowerTriangular(W), Inf)
      integrator.eigen_est = (opn + one(opn)) / dtgamma # TODO: better estimate
    end
    return nothing
  end

  # check if we need to update J or W
  new_jac, new_W = do_newJW(integrator, alg, nlsolver, repeat_step)

  if new_jac && isnewton(lcache)
    lcache.J_t = t
    if isdae
      lcache.uf.α = nlsolver.α
      lcache.uf.invγdt = inv(dtgamma)
      lcache.uf.tmp = nlsolver.tmp
    end
  end

  # calculate W
  if W isa WOperator
    isnewton(nlsolver) || DiffEqBase.update_coefficients!(W,uprev,p,t) # we will call `update_coefficients!` in NLNewton
    W.transform = W_transform; set_gamma!(W, dtgamma)
  else # concrete W using jacobian from `calc_J!`
    islin, isode = islinearfunction(integrator)
    islin ? (J = isode ? f.f : f.f1.f) : ( new_jac && (calc_J!(J, integrator, lcache)) )
    new_W && !(isdae) && jacobian2W!(W, mass_matrix, dtgamma, J, W_transform)
  end
  if isnewton(nlsolver)
    set_new_W!(nlsolver, new_W)
    if new_jac && isdae
      set_W_γdt!(nlsolver, nlsolver.α * inv(dtgamma))
    elseif new_W && !isdae
      set_W_γdt!(nlsolver, dtgamma)
    end
  end
  new_W && (integrator.destats.nw += 1)
  return nothing
end

function calc_W(integrator, cache, dtgamma, repeat_step, W_transform=false)
  @unpack t,uprev,p,f = integrator
  isdae = integrator.alg isa DAEAlgorithm
  if !isdae
    mass_matrix = integrator.f.mass_matrix
  end
  isarray = uprev isa AbstractArray
  # calculate W
  is_compos = integrator.alg isa CompositeAlgorithm
  islin, isode = islinearfunction(integrator)
  if islin
    J = isode ? f.f : f.f1.f # unwrap the Jacobian accordingly
    W = WOperator(mass_matrix, dtgamma, J, false; transform=W_transform)
  elseif DiffEqBase.has_jac(f)
    J = f.jac(uprev, p, t)
    if !isa(J, DiffEqBase.AbstractDiffEqLinearOperator)
      J = DiffEqArrayOperator(J)
    end
    W = WOperator(mass_matrix, dtgamma, J, false; transform=W_transform)
    integrator.destats.nw += 1
  else
    integrator.destats.nw += 1
    J = calc_J(integrator, cache)
    if isdae
      W = J
    else
      W_full = W_transform ? -mass_matrix*inv(dtgamma) + J :
                             -mass_matrix + dtgamma*J
      W = W_full isa Number ? W_full : lu(W_full)
    end
  end
  (W isa WOperator && unwrap_alg(integrator, true) isa NewtonAlgorithm) && (W = DiffEqBase.update_coefficients!(W,uprev,p,t)) # we will call `update_coefficients!` in NLNewton
  is_compos && (integrator.eigen_est = isarray ? opnorm(J, Inf) : abs(J))
  return W
end

function calc_rosenbrock_differentiation!(integrator, cache, dtd1, dtgamma, repeat_step, W_transform)
  calc_tderivative!(integrator, cache, dtd1, repeat_step)
  nlsolver = nothing
  # we need to skip calculating `W` when a step is repeated
  repeat_step || calc_W!(cache.W, integrator, nlsolver, cache, dtgamma, repeat_step, W_transform)
  return nothing
end

# update W matrix (only used in Newton method)
update_W!(integrator, cache, dtgamma, repeat_step) =
  update_W!(cache.nlsolver, integrator, cache, dtgamma, repeat_step)

function update_W!(nlsolver::AbstractNLSolver, integrator, cache::OrdinaryDiffEqMutableCache, dtgamma, repeat_step)
  if isnewton(nlsolver)
    calc_W!(get_W(nlsolver), integrator, nlsolver, cache, dtgamma, repeat_step, true)
  end
  nothing
end

function update_W!(nlsolver::AbstractNLSolver, integrator, cache, dtgamma, repeat_step)
  if isnewton(nlsolver)
    isdae = integrator.alg isa DAEAlgorithm
    new_jac, new_W = true, true
    if isdae && new_jac
      lcache = nlsolver.cache
      lcache.uf.α = nlsolver.α
      lcache.uf.invγdt = inv(dtgamma)
      lcache.uf.tmp = @. nlsolver.tmp
      lcache.uf.uprev = @. integrator.uprev
    end
    nlsolver.cache.W = calc_W(integrator, nlsolver.cache, dtgamma, repeat_step, true)
    #TODO: jacobian reuse for oop
    new_jac && (nlsolver.cache.J_t = integrator.t)
    set_new_W!(nlsolver, new_W)
    if new_jac && isdae
      set_W_γdt!(nlsolver, nlsolver.α * inv(dtgamma))
    elseif new_W && !isdae
      set_W_γdt!(nlsolver, dtgamma)
    end
  end
  nothing
end

function build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,::Val{IIP}) where IIP
  islin, isode = islinearfunction(f, alg)
  if f.jac_prototype isa DiffEqBase.AbstractDiffEqLinearOperator
    W = WOperator(f, dt, IIP)
    J = W.J
  elseif IIP && f.jac_prototype !== nothing
    J = similar(f.jac_prototype)
    W = similar(J)
  elseif islin || (!IIP && DiffEqBase.has_jac(f))
    J = islin ? (isode ? f.f : f.f1.f) : f.jac(uprev, p, t) # unwrap the Jacobian accordingly
    if !isa(J, DiffEqBase.AbstractDiffEqLinearOperator)
      J = DiffEqArrayOperator(J)
    end
    W = WOperator(f.mass_matrix, dt, J, IIP)
  else
    J = false .* _vec(u) .* _vec(u)'
    isdae = alg isa DAEAlgorithm
    W = if isdae
      J
    elseif IIP
      similar(J)
    else
      W = if u isa StaticArray
      lu(J)
      elseif u isa Number
        u
      else
        LU{LinearAlgebra.lutype(uEltypeNoUnits)}(Matrix{uEltypeNoUnits}(undef, 0, 0),
                                                     Vector{LinearAlgebra.BlasInt}(undef, 0),
                                                     zero(LinearAlgebra.BlasInt))
      end
    end # end W
  end
  return J, W
end

build_uf(alg::OrdinaryDiffEqAlgorithm,nf,t,p,::Val{true}) =
  UJacobianWrapper(nf,t,p)
build_uf(alg::OrdinaryDiffEqAlgorithm,nf,t,p,::Val{false}) =
  UDerivativeWrapper(nf,t,p)


text:
using PMI
using Test

@test PMI.initialized() == false

PMI.init()

@test PMI.initialized()

nranks = PMI.get_size()
rank = PMI.get_rank()

@test nranks >= 1
@test 0 <= rank < nranks

PMI.barrier()

@test PMI.get_clique_size() >= 1
@test rank ∈ PMI.get_clique_ranks()

PMI.finalize()


text:
using JSON

export Asset, Async, Sync

"""
    Asset(url)
    Asset(name, url)
    Asset(name => url)
    Asset(filetype, name, url)

A browser asset that can be loaded by WebIO.

The `url` parameter may be either a remote URL or a local filepath (which will
be served via AssetRegistry.jl).
All of the following are valid `url` values.
  * `https://unpkg.com/react@16/umd/react.development.js`
  * `//unpkg.com/react@16/umd/react.development.js`
  * `./path/to/foo.js`

By default, the filetype is guessed as the extension of the specified url
(only `"js"`, `"css"`, and `"html"` are currently supported) but may be
specified if nonstandard extensions are in use.
"""
struct Asset
    filetype::String
    name::Union{Nothing, String}
    url::String
end

Asset(name, url) = Asset(getextension(url), name, url)
Asset(url) = Asset(nothing, url)
Asset(x::Pair) = Asset(string(x[1]), x[2])
Asset(x::Asset) = x
function Asset(spec::Dict)
    Base.depwarn(
        "`Asset(::Dict)` is deprecated, use `Asset(url)` instead.`",
        :webio_asset_dict_constructor,
    )

    spec = Dict(string(k) => v for (k, v) in spec)
    url = get(spec, "url") do
        error("Invalid Asset dict specification: missing key \"url\".")
    end
    filetype = get(spec, "type", getextension(url))
    name = get(spec, "name", nothing)

    return Asset(filetype, name, url)
end

function JSON.lower(asset::Asset)
    Dict(
         "type" => asset.filetype,
         "url"  => dep2url(asset.url),
         "name" => asset.name,
    )
end

function tojs(asset::Asset)
    return js"WebIO.importResource($(JSON.lower(asset)))"
end

"""
    Sync(assets...)

A group ("block") of assets that will be imported sequentially (synchronously)
in the browser.
This is useful when dependencies have side effects that must be executed in
order.
The elements of `assets` may either be [`Asset`](@ref)s themselves, any valid
constructor for an [`Asset`](@ref), or a [`Sync`](@ref) or [`Async`](@ref).

If the imports do not need to be imported sequentially, use [`Async`](@ref)
instead.

# Examples
```julia-repl
julia> WebIO.Sync(Asset("foo.js"), "bar" => "bar.js")
Sync(Asset[Asset("js", nothing, "foo.js"), Asset("js", "bar", "bar.js")])
```
"""
struct Sync
    # This is untyped because we can't define a union type that includes all of
    # Asset, Sync, and Async that is then used within Sync and Async.
    imports::Array{Any}

    Sync(imports::Array) = new([ensure_asset(asset) for asset in imports])
end
Sync(assets...) = Sync([assets...])

"""
    Async(assets...)

A group ("block") of assets that will be imported with no specified order
(asynchronously) in the browser.
This is useful when dependencies can be imported in any order.
The elements of `assets` may either be [`Asset`](@ref)s themselves, any valid
constructor for an [`Asset`](@ref), or a [`Sync`](@ref) or [`Async`](@ref).

If the imports need to be imported sequentially, use [`Sync`](@ref) instead.
"""
struct Async
    # See comment about (lack of) typing in Sync above.
    imports::Array{Any}

    Async(imports::Array) = new([ensure_asset(asset) for asset in imports])
end
Async(assets...) = Async([assets...])

# This allows js"await $(Async([....]))" on a tree of assets!
function tojs(asset::Union{Async,Sync})
    return js"WebIO.importBlock($(lowerassets(asset)))"
end

# The output of lowerassets is initially sent with the Scope
# this should trigger loading of the assets before onmount callbacks
lowerassets(x) = JSON.lower(Asset(x))
lowerassets(x::Asset) = JSON.lower(x)
lowerassets(x::Async) = Dict(
    "type" => "async_block",
    "data" => map(lowerassets, x.imports),
)
lowerassets(x::AbstractArray) = lowerassets(Async(x))
lowerassets(x::Sync) = Dict(
    "type"=>"sync_block",
    "data" => map(lowerassets, x.imports),
)

"""
    ensure_asset(asset)

Ensure that `asset` is a valid [`Asset`](@ref) or [`Async`](@ref) or
[`Sync`](@ref).
If it's not, calls the [`Asset`](@ref) constructor on the argument.
"""
ensure_asset(asset::Union{Asset, Sync, Async}) = asset
ensure_asset(arg) = Asset(arg)

"""
    getextension(path)

Get the file extension of the path.
The extension is defined to be the bit after the last dot, excluding any query
string.

# Examples
```julia-repl
julia> WebIO.getextension("foo.bar.js")
"js"
julia> WebIO.getextension("https://my-cdn.net/foo.bar.css?version=1")
"css"
```
"""
getextension(path) = last(split(first(split(path, "?")), "."))

"""
    islocal(path)

Determine whether or not the specified path is a local filesystem path (and not
a remote resource that is hosted on, for example, a CDN).
"""
islocal(path) = !any(startswith.(path, ("//", "https://", "http://", "ftp://")))

"""
    path2url(path)

Register the specified path with AssetRegistry and return the url that
corresponds to that path.

# Examples
```julia-repl
julia> WebIO.path2url(expanduser("~/Documents/foo.js"))
"/assetserver/bd67c48cbe6388c22f85faef9840ff9a2dfc1df6-foo.js"
```
"""
function path2url(path::AbstractString)
    if startswith(path, "/pkg/")
        Base.depwarn(
            "/pkg/ URLs are deprecated, specify files using their absolute paths.",
            :webio_pkg_urls,
        )
        return path
    elseif startswith(path, "/assetserver/") && haskey(AssetRegistry.registry, path)
        return path
    elseif isfile(abspath(path))
        # TODO: this should be implemented in AssetRegistry, not here.
        path = abspath(path)
        # first lookup to see if any of the file itself or any of the parent
        # directories are registered.
        AssetRegistry.isregistered(path) && AssetRegistry.getkey(path)
        cur_path = path
        while true
            if AssetRegistry.isregistered(cur_path) && isdir(cur_path)
                key = AssetRegistry.getkey(cur_path)
                url = key * "/" * replace(replace(path, cur_path => ""), r"^\/"=>"")
                return url
            end
            cur_path1 = dirname(cur_path)
            if cur_path1 == cur_path
                # this means we have reached root directory,
                # and none of the parents are in registry
                # register the original path uniquely
                return AssetRegistry.register(path)
            end
            cur_path = cur_path1
        end
    end
    error("Invalid dependency (must be a url, file, or AssetRegistry path): $path")
end

"""
    dep2url(dep)

Return the URL where the dependency can be fetched from.
If the dependency is a URL (e.g. hosted at a CDN), the same URL is returned.
Otherwise, the depency is registered with AssetRegistry.jl and served from
there.
"""
function dep2url(dep::AbstractString)
    # if is an url, we are done :)
    islocal(dep) || return dep
    query_parts = split(dep, "?") # remove anything after ?
    file_path = first(query_parts)
    query_part = length(query_parts) >= 2 ? "?" * query_parts[2] : ""
    url = path2url(file_path)
    return string(baseurl[], url, query_part)
end


text:
function generate_logdata()
    logdata_path = mktempdir()

    # generate fake data
    s = now() - Day(30)
    e = now()
    dt = Minute(1)

    dates = s:dt:e
    # sin/cos require numeric values. Pass the millisecond
    # value to the function
    temp = [sin(x.instant.periods.value) / 100 for x in dates]
    hum = [cos(x.instant.periods.value) / 110 for x in dates]

    writedlm(joinpath(logdata_path, "test.log"), [dates temp hum])

    return logdata_path
end


text:
cuproj(state, dims::Union{Vector, NTuple}) = cu.(proj(state, dims))

function tensor(ψ::CuMPS)
    devψ = MPS(ψ)
    t = tensor(devψ)
    cu(t)
end


text:
# # Vlasov-Ampere
#
#md # [`notebook`](@__NBVIEWER_ROOT_URL__notebooks/vlasov-ampere.ipynb)
#
# Compute Landau damping by solving Vlasov-Ampere system.
#
# ```math
#  \frac{∂f}{∂t} + υ \frac{∂f}{∂x} - E(t,x) \frac{∂f}{∂υ} = 0
# ```
#
# ```math
# \frac{∂E}{∂t} = - J = ∫ fυ dυ
# ```
#
# ## Algorithm
#
#  - For each ``j`` compute discrete Fourier transform in ``x`` of
#    ``(x_i,υ_j)`` yielding ``f_k^n(υ_j)``,
#  - For `` k ≂̸ 0 ``
#      * Compute ``f^{n+1}_k(υ_j) = e^{−2iπ k υ Δt/L} f_n^k(υ_j), ``
#      * Compute ``ρ_k^{n+1} = Δ υ ∑_j􏰄 f^{n+1}_k(υ_j), ``
#      * Compute ``E^{n+1}_k = ρ^{n+1}_k L/(2iπkϵ_0), ``
#  - For ``k = 0`` do nothing:
# ```math
# f_{n+1}(υ_j) = f^n_k(υ_j), E^{n+1}_k = E^n_k.
# ```
#  - Perform in2erse discrete Fourier transform of ``E^{n+1}_k`` and for each
#    ``j`` of ``f^{n+1}_k (υ_j)``.
#

import Splittings: advection!, Ampere, UniformMesh
import Splittings: @Strang, compute_rho, compute_e
using Plots, LinearAlgebra

#-

function push_t!( f, fᵀ, mesh1, mesh2, e,  dt)

    advection!( f, fᵀ, mesh1, mesh2, e,  dt, Ampere(), 1 )

end

#-

function push_v!(f, fᵀ, mesh1, mesh2, e, dt)

    advection!( f, fᵀ, mesh1, mesh2, e, dt, Ampere(), 2)

end

#-

function vm1d( n1, n2, x1min, x1max, x2min, x2max , tf, nt)

    mesh1 = UniformMesh(x1min, x1max, n1, endpoint=false)
    mesh2 = UniformMesh(x2min, x2max, n2, endpoint=false)

    x = mesh1.points
    v = mesh2.points
    ϵ, kx = 0.001, 0.5

    f = zeros(Complex{Float64},(n1,n2))
    fᵀ= zeros(Complex{Float64},(n2,n1))

    f .= (1.0.+ϵ*cos.(kx*x))/sqrt(2π) .* transpose(exp.(-0.5*v.*v))
    transpose!(fᵀ,f)

    e = zeros(Complex{Float64},n1)

    ρ  = compute_rho(mesh2, f)
    e .= compute_e(mesh1, ρ)

    nrj = Float64[]

    dt = tf / nt

    for i in 1:nt

	push!(nrj, 0.5*log(sum(real(e).^2)*mesh1.step))

        @Strang(  push_v!(f, fᵀ, mesh1, mesh2, e, dt),
                  push_t!(f, fᵀ, mesh1, mesh2, e, dt)
               )

    end
    nrj
end

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

n1, n2 = 32, 64
x1min, x1max =  0., 4π
x2min, x2max = -6., 6.
tf = 50
nt = 500

t = range(0,stop=tf,length=nt)
nrj = vm1d(n1, n2, x1min, x1max, x2min, x2max, tf, nt)
plot(t, nrj)
plot!(t, -0.1533*t.-5.50)
savefig("va-plot.png"); nothing

using Test                      #src
@testset "Vlasov-Ampere" begin  #src
@test length(nrj) > 0           #src
end                             #src

#md # ![](va-plot.png)


text:
using ImageCore, Colors, ColorVectorSpace
using Test, Statistics

# Different access patterns (getindex)
function mysum_elt_boundscheck(A)
    s = zero(eltype(A))
    for a in A
        s += a
    end
    s
end
function mysum_index_boundscheck(A)
    s = zero(eltype(A))
    for I in eachindex(A)
        s += A[I]
    end
    s
end
function mysum_elt_inbounds(A)
    s = zero(eltype(A))
    @inbounds for a in A
        s += a
    end
    s
end
function mysum_index_inbounds_simd(A)
    s = zero(eltype(A))
    @inbounds @simd for I in eachindex(A)
        s += A[I]
    end
    s
end
# setindex!
function myfill1!(A, val)
    f = convert(eltype(A), val)
    for I in eachindex(A)
        A[I] = f
    end
    A
end
function myfill2!(A, val)
    f = convert(eltype(A), val)
    @inbounds @simd for I in eachindex(A)
        A[I] = f
    end
    A
end

# Rather than using BenchmarkTools (and thus run one test repeatedly,
# accumulating timings), we run the same test interleaving the two
# array types. This is designed to reduce the risk of spurious
# failure, particularly on shared machines like Travis where they may
# get "distracted" by other tasks
function test_getindex(f, ar, cv, n)
    t_ar = Array{Float64}(undef, n)
    t_cv = Array{Float64}(undef, n)
    f_ar = Ref(f(ar))
    f_cv = Ref(f(cv))
    for i = 1:n
        t_ar[i] = (tstart = time(); f_ar[] = f(ar); time()-tstart)
        t_cv[i] = (tstart = time(); f_cv[] = f(cv); time()-tstart)
    end
    median(t_ar), median(t_cv), f_ar
end
function test_setindex(f, ar, cv, n)
    t_ar = Array{Float64}(undef, n)
    t_cv = Array{Float64}(undef, n)
    for i = 1:n
        t_ar[i] = @elapsed f(ar, zero(eltype(ar)))
        t_cv[i] = @elapsed f(cv, zero(eltype(cv)))
    end
    median(t_ar), median(t_cv)
end

ssz = (1000,1000)
c = rand(RGB{Float64}, ssz...)
a = copy(reinterpretc(Float64, c))
vchan = channelview(c)
vcol = colorview(RGB, a)
cc_getindex_funcs = (mysum_elt_boundscheck,
                     mysum_index_boundscheck,
                     mysum_elt_inbounds,
                     mysum_index_inbounds_simd)
cc_setindex_funcs = (myfill1!,
                     myfill2!)
chanvtol = Dict(mysum_index_inbounds_simd => 20,   # @simd doesn't work for ChannelView :(
                mysum_elt_boundscheck => 20,
                myfill1! => 20,                    # crappy setindex! performance
                myfill2! => 20)
chanvdefault = 10
colvtol = Dict(mysum_elt_boundscheck=>5,
               mysum_index_boundscheck=>5)
colvdefault = 3

@info "Benchmark tests are warnings for now"
# @testset "benchmarks" begin
for (suite, testf) in ((cc_getindex_funcs, test_getindex),
                       (cc_setindex_funcs, test_setindex))
    for f in suite
        # Channelview
        t_ar, t_cv = testf(f, a, vchan, 10^2)
        tol = haskey(chanvtol, f) ? chanvtol[f] : chanvdefault
        if t_cv >= tol*t_ar
            @warn "ChannelView: failed on $f, time ratio $(t_cv/t_ar), tol $tol"
        end
        # @test t_cv < tol*t_ar
        # ColorView
        t_ar, t_cv = testf(f, c, vcol, 10^2)
        tol = haskey(colvtol, f) ? colvtol[f] : colvdefault
        if t_cv >= tol*t_ar
            @warn "ColorView: failed on $f, time ratio $(t_cv/t_ar), tol $tol"
        end
        # @test t_cv < tol*t_ar
    end
end
# end


text:
"""
    ExtendedHubbardReal1D(address; u=1.0, v=1.0, t=1.0)

Implements the extended Hubbard model on a one-dimensional chain in real space.

```math
\\hat{H} = -t \\sum_{\\langle i,j\\rangle} a_i^† a_j + \\frac{u}{2}\\sum_i n_i (n_i-1) + v \\sum_{\\langle i,j\\rangle} n_i n_j
```

# Arguments

* `address`: the starting address.
* `u`: on-site interaction parameter
* `v`: the next-neighbor interaction
* `t`: the hopping strength

"""
struct ExtendedHubbardReal1D{TT,A<:AbstractFockAddress,U,V,T} <: AbstractHamiltonian{TT}
    add::A
end

# addr for compatibility.
function ExtendedHubbardReal1D(addr; u=1.0, v=1.0, t=1.0)
    U, V, T = promote(float(u), float(v), float(t))
    return ExtendedHubbardReal1D{typeof(U),typeof(addr),U,V,T}(addr)
end

function Base.show(io::IO, h::ExtendedHubbardReal1D)
    print(io, "ExtendedHubbardReal1D($(h.add); u=$(h.u), v=$(h.v), t=$(h.t))")
end

function starting_address(h::ExtendedHubbardReal1D)
    return getfield(h, :add)
end

LOStructure(::Type{<:ExtendedHubbardReal1D{<:Real}}) = IsHermitian()

Base.getproperty(h::ExtendedHubbardReal1D, s::Symbol) = getproperty(h, Val(s))
Base.getproperty(h::ExtendedHubbardReal1D, ::Val{:add}) = getfield(h, :add)
Base.getproperty(h::ExtendedHubbardReal1D{<:Any,<:Any,U}, ::Val{:u}) where U = U
Base.getproperty(h::ExtendedHubbardReal1D{<:Any,<:Any,<:Any,V}, ::Val{:v}) where V = V
Base.getproperty(h::ExtendedHubbardReal1D{<:Any,<:Any,<:Any,<:Any,T}, ::Val{:t}) where T = T

function num_offdiagonals(::ExtendedHubbardReal1D, address::BoseFS)
    return 2 * num_occupied_modes(address)
end

"""
    extended_bose_hubbard_interaction(address, m)

Compute the on-site product sum_j n_j(n_j-1) and the next neighbour term
sum_j n_j n_{j+1} with periodic boundary conditions.
"""
function extended_bose_hubbard_interaction(b::BoseFS{<:Any,M}) where M
    address = b.bs
    # compute the diagonal matrix element of the Extended Bose Hubbard Hamiltonian
    # currently this ammounts to counting occupation numbers of modes
    bhmmatrixelementint = 0
    ebhmmatrixelementint = 0
    bosonnumber2=0
    bosonnumber1 = trailing_ones(address) # count how many bosons inside
    # surpsingly it is faster to not check whether this is nonzero and do the
    # following operations anyway
    bhmmatrixelementint += bosonnumber1 * (bosonnumber1-1)
    firstbosonnumber = bosonnumber1 #keap on memory the boson number of the first
    #to do the calculation with the last boson
    address >>>= bosonnumber1 # remove the countedmode
    address >>>= 1
    for i in 1:M-1
        # proceed to next occupied mode
        bosonnumber2 = trailing_ones(address) # count how many bosons inside
        # surpsingly it is faster to not check whether this is nonzero and do the
        # following operations anyway
        address >>>= bosonnumber2 # remove the countedmode
        ebhmmatrixelementint += bosonnumber2 * (bosonnumber1)
        bhmmatrixelementint+= bosonnumber2 * (bosonnumber2-1)
        bosonnumber1=bosonnumber2
        address >>>= 1
    end
    ebhmmatrixelementint+= bosonnumber2 * firstbosonnumber  #periodic bondary condition
    return ebhmmatrixelementint , bhmmatrixelementint
end

function diagonal_element(h::ExtendedHubbardReal1D, b::BoseFS)
    ebhinteraction, bhinteraction = extended_bose_hubbard_interaction(b)
    return h.u * bhinteraction / 2 + h.v * ebhinteraction
end

function get_offdiagonal(h::ExtendedHubbardReal1D, add::BoseFS, chosen)
    naddress, onproduct = hopnextneighbour(add, chosen)
    return naddress, - h.t * sqrt(onproduct)
end


text:
"""
    Sherlock(optimizer, ϵ::Float64)

Sherlock combines local and global search to estimate the range of the output node.

# Problem requirement
1. Network: any depth, ReLU activation, single output
2. Input: hpolytope and hyperrectangle
3. Output: hyperrectangle (1d interval)

# Return
`CounterExampleResult` or `ReachabilityResult`

# Method
Local search: solve a linear program to find local optima on a line segment of the piece-wise linear network.
Global search: solve a feasibilty problem using MILP encoding (default calling NSVerify).
- `optimizer` default `GLPKSolverMIP()`
- `ϵ` is the margin for global search, default `0.1`.

# Property
Sound but not complete.

# Reference
[S. Dutta, S. Jha, S. Sanakaranarayanan, and A. Tiwari,
"Output Range Analysis for Deep Neural Networks,"
*ArXiv Preprint ArXiv:1709.09130*, 2017.](https://arxiv.org/abs/1709.09130)

[https://github.com/souradeep-111/sherlock](https://github.com/souradeep-111/sherlock)
"""
@with_kw struct Sherlock <: Solver
    optimizer = GLPK.Optimizer
    ϵ::Float64 = 0.1
end

function solve(solver::Sherlock, problem::Problem)
    isbounded(problem.input) || UnboundedInputError("Temprarily, Sherlock cannot accept unbounded inputs")

    (x_u, u) = output_bound(solver, problem, :max)
    (x_l, l) = output_bound(solver, problem, :min)
    println("bounds: [", l, ", ", u, "]")
    bound = Hyperrectangle(low = [l], high = [u])
    reach = Hyperrectangle(low = [l - solver.ϵ], high = [u + solver.ϵ])

    output = problem.output
    if reach ⊆ output
        return ReachabilityResult(:holds, [reach])
    end
    high(bound) > high(output) && return CounterExampleResult(:violated, x_u)
    low(bound)  < low(output)  && return CounterExampleResult(:violated, x_l)
    return ReachabilityResult(:unknown, [reach])
end

function output_bound(solver::Sherlock, problem::Problem, type::Symbol)
    opt = solver.optimizer
    x = an_element(problem.input)
    while true
        (x, bound) = local_search(problem, x, opt, type)
        bound_ϵ = bound + ifelse(type == :max, solver.ϵ, -solver.ϵ)
        (x_new, bound_new, feasible) = global_search(problem, bound_ϵ, opt, type)
        feasible || return (x, bound)
        (x, bound) = (x_new, bound_new)
    end
end

function local_search(problem::Problem, x::Vector{Float64}, optimizer, type::Symbol)
    nnet = problem.network

    model = Model(optimizer)
    model[:δ] = get_activation(nnet, x)
    z = init_vars(model, nnet, :z, with_input=true)
    add_set_constraint!(model, problem.input, first(z))
    encode_network!(model, nnet, StandardLP())

    # the gradient wrt to the input is a `1xn` matrix,
    # because the output is 1-d. Therefore o[1] is the
    # only entry in g*z₀
    gradient = get_gradient(nnet, x)
    o = gradient * first(z)
    index = ifelse(type == :max, 1, -1)
    @objective(model, Max, index * o[1])
    optimize!(model)

    x_new = value(first(z))
    bound_new = compute_output(nnet, x_new)
    return (x_new, bound_new[1])
end

function global_search(problem::Problem, bound::Float64, optimizer, type::Symbol)
    index = ifelse(type == :max, 1.0, -1.0)
    h = HalfSpace([index], index * bound)
    output_set = HPolytope([h])
    problem_new = Problem(problem.network, problem.input, output_set)
    solver  = NSVerify(optimizer = optimizer)
    result  = solve(solver, problem_new)
    if result.status == :violated
        x = result.counter_example
        bound = compute_output(problem.network, x)
        return (x, bound[1], true)
    else
        return ([], 0.0, false)
    end
end


text:
## solve Lx=b where x and b are dense.  x=b on input, solution on output.
function Base.A_ldiv_B!{T}(L::LowerTriangular{T,SparseMatrixCSC{T}}, x::StridedVector{T})
    Ld = L.data
    n = Base.LinAlg.chksquare(Ld)
    length(x) == n || throw(DimensionMismatch())
    Lp = Ld.colptr; Li = rowvals(Ld); Lx = nonzeros(Ld)
    for j in 1:n
        x[j] /= Lx[Lp[j]]
        for p in (Lp[j] + 1):(Lp[j+1] - 1)
            x[Li[p]] -= Lx[p] * x[j]
        end
    end
    x
end

## solve L'x=b where x and b are dense.  x=b on input, solution on output.
function Ac_ldiv_B!{T}(L::LowerTriangular{T,SparseMatrixCSC{T}}, x::StridedVector{T})
    Ld = L.data
    n = Base.LinAlg.chksquare(Ld)
    length(x) == n || throw(DimensionMismatch())
    Lp = Ld.colptr; Li = rowvals(Ld); Lx = nonzeros(Ld)
    for j in reverse(1:n)
        for p in (Lp[j] + 1):(Lp[j+1] - 1)
            x[j] -= Lx[p] * x[Li[p]]
        end
        x[j] /= Lx[Lp[j]]
    end
    x
end

## solve Ux=b where x and b are dense.  x=b on input, solution on output.
function A_ldiv_B!{T}(U::UpperTriangular{T,SparseMatrixCSC{T}}, x::StridedVector{T})
    n = Base.LinAlg.chksquare(U)
    length(x) == n || throw(DimensionMismatch())
    Ud = U.data
    Up = Ud.colptr; Ui = rowvals(Ud); Ux = nonzeros(Ud)
    for j in reverse(1:n)
        x[j] /= Ux[Up[j+1]-1]
        for p in Up[j]:(Up[j+1]-2)
            x[Ui[p]] -= Ux[p] * x[j]
        end
    end
    x
end

## solve U'x=b where x and b are dense.  x=b on input, solution on output.
function Ac_ldiv_B!{T}(U::UpperTriangular{T,SparseMatrixCSC{T}}, x::StridedVector{T})
    n = Base.LinAlg.chksquare(U)
    length(x) == n || throw(DimensionMismatch)
    Ud = U.data
    Up = Ud.colptr; Ui = rowvals(Ud); Ux = nonzeros(Ud)
    for j in 1:n
        for p in Up[j]:(Up[j+1]-2)
            x[j] -= Ux[p] * x[Ui[p]]
            end
        x[j] /= Ux[Up[j+1]-1]
    end
    x
end

## 1-norm of a sparse matrix = max (sum (abs (A))), largest column sum
function js_norm(A::SparseMatrixCSC)
    Ax  = nonzeros(A)
    nrm = zero(eltype(A))
    for j in 1:size(A,2)
        nrm = max(nrm, sumabs(sub(Ax,nzrange(A,j))))
    end
    nrm
end

## solve Gx=b(:,k), where G is either upper (lo=0) or lower (lo=1) triangular
function js_spsolve{Tv,Ti<:Integer}(G::LowerTriangular{Tv,SparseMatrixCSC{Tv}},
                                   B::SparseMatrixCSC{Tv}, k,
                                   xi::Vector{Ti}, x::Vector{Tv}, pinv::Vector)
    n = Base.LinAlg.chksquare(G)
    Gd = G.data
    Gp = Gd.colptr; Gi = rowvals(Gd); Gx = nonzeros(Gd)
    Bp = B.colptr; Bi = rowvals(B); Bx = nonzeros(B)
    ## must have reach return a 1-based index
    top = js_reach(G, B, k, xi, pinv)   # xi[top..n-1]=Reach(B(:,k))
    ## get xi from reach?
    for p in top:n x[xi[p]] = 0. end    # clear x
    for p in Bp[k]:(Bp[k+1]-1) x[Bi[p]] = Bx[p] end # scatter B
    for px in top:n
        j = xi[px]                       # x[j] is nonzero
        J = pinv[j]                      # j maps to col J of G
        if J < 0 continue end            # column J is empty
        x[j] /= Gx[lo? Gp[J]: Gp[J+1]-1] # x(j) /= G[j,j]
        p = lo ? (Gp[J]+1) : Gp[J]       # lo: L[j,j] 1st entry
        q = lo ? Gp[J+1] : (Gp[J+1]-1)   # up: U[j,j] last entry
        for pp in p:q
            x[Gi[pp]] -= Gx[pp] * x[j]   # x[i] -= G[i,j] * x[j]
        end
    end
    top
end


text:
module Spiderman
using Gumbo
using Requests
using AbstractTrees
import HTTPClient

include("css.jl")
export @css_str, compile_css, collect

include("helpers.jl")
export text, parent, parse

include("http.jl")

end # module


text:
using NESSie.Format
using NESSie.Format: readoff_nodes, readoff_elements

testfiles = ((mktemp()..., (0,0)), # empty file
             (mktemp()..., (3,2))) # dummy file
write(testfiles[2][2], """
OFF
3 2 5
 1.0000000000e+00     0.0000000000e+00     0.0000000000e+00
 0.0000000000e+00     1.0000000000e+00     0.0000000000e+00
 0.0000000000e+00     0.0000000000e+00     1.0000000000e+00
3 0 1 2
3 2 0 1
""")

@testset "readoff" begin
    try
        for (fname, fh, len) in testfiles, fn in ("single", "bulk"), T in testtypes
            seekstart(fh)
            nodes = (); elements = ()
            if fn == "single"
                # check single methods
                readline(fh)
                readline(fh)
                nodes = readoff_nodes(fh, len[1], T)
                elements = readoff_elements(fh, len[2], nodes)
            elseif fn == "bulk"
                # check bulk method
                model = readoff(fh, T)
                nodes, elements = (model.nodes, model.elements)
            end
            # check return types
            @test typeof(nodes) == Vector{Vector{T}}
            @test typeof(elements) == Vector{Triangle{T}}

            # check lengths
            @test length(nodes) == len[1]
            @test length(elements) == len[2]

            # check content
            if fname == testfiles[2][1]
                @test nodes[1] == T[1, 0, 0]
                @test nodes[2] == T[0, 1, 0]
                @test nodes[3] == T[0, 0, 1]
                @test elements[1].v1 === nodes[1]
                @test elements[1].v2 === nodes[2]
                @test elements[1].v3 === nodes[3]
                @test elements[2].v1 === nodes[3]
                @test elements[2].v2 === nodes[1]
                @test elements[2].v3 === nodes[2]
            end
        end
    finally
        [begin close(fh); rm(fn); end for (fn, fh, _) in testfiles]
    end
end


text:
function striplinecomment{T<:String,U<:String}(a::T, cchars::U="#;")
    b = strip(a)
    0 < length(cchars) || return b
    for c in cchars
        r = Regex(@sprintf "\\%c.*" c)
        b = replace(b, r, "")
    end
    strip(b)
end

tests = {"apples, pears # and bananas",
         "apples, pears ; and bananas",
         "  apples, pears & bananas   ",
         " # "}

for t in tests
    s = striplinecomment(t)
    println("Testing \"", t, "\":")
    println("    \"", s, "\"")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
is_func_call(x) = false
is_func_call(x::EXPR{Call}) = true
is_func_call(x::EXPR{InvisBrackets}) = is_func_call(x.args[2])
is_func_call(x::UnaryOpCall) = true
is_func_call(x::BinaryOpCall) = true
function is_func_call(x::BinarySyntaxOpCall)
    if is_decl(x.op)
        return is_func_call(x.arg1)
    else
        return false
    end
end
function is_func_call(x::WhereOpCall)
    return is_func_call(x.arg1)
end

is_assignment(x) = x isa BinarySyntaxOpCall && is_eq(x.op)

# OPERATOR
is_exor(x) = x isa OPERATOR && x.kind == Tokens.EX_OR && x.dot == false
is_decl(x) = x isa OPERATOR && x.kind == Tokens.DECLARATION
is_issubt(x) = x isa OPERATOR && x.kind == Tokens.ISSUBTYPE
is_issupt(x) = x isa OPERATOR && x.kind == Tokens.ISSUPERTYPE
is_and(x) = x isa OPERATOR && x.kind == Tokens.AND && x.dot == false
is_not(x) = x isa OPERATOR && x.kind == Tokens.NOT && x.dot == false
is_plus(x) = x isa OPERATOR && x.kind == Tokens.PLUS && x.dot == false
is_minus(x) = x isa OPERATOR && x.kind == Tokens.MINUS && x.dot == false
is_star(x) = x isa OPERATOR && x.kind == Tokens.STAR && x.dot == false
is_eq(x) = x isa OPERATOR && x.kind == Tokens.EQ && x.dot == false
is_dot(x) = x isa OPERATOR && x.kind == Tokens.DOT
is_ddot(x) = x isa OPERATOR && x.kind == Tokens.DDOT
is_dddot(x) = x isa OPERATOR && x.kind == Tokens.DDDOT
is_pairarrow(x) = x isa OPERATOR && x.kind == Tokens.PAIR_ARROW && x.dot == false
is_in(x) = x isa OPERATOR && x.kind == Tokens.IN && x.dot == false
is_elof(x) = x isa OPERATOR && x.kind == Tokens.ELEMENT_OF && x.dot == false
is_colon(x) = x isa OPERATOR && x.kind == Tokens.COLON
is_prime(x) = x isa OPERATOR && x.kind == Tokens.PRIME
is_cond(x) = x isa OPERATOR && x.kind == Tokens.CONDITIONAL
is_where(x) = x isa OPERATOR && x.kind == Tokens.WHERE
is_anon_func(x) = x isa OPERATOR && x.kind == Tokens.ANON_FUNC

# PUNCTUATION
is_punc(x) = x isa PUNCTUATION &&
    x.kind == Tokens.COMMA &&
    x.kind == Tokens.LPAREN &&
    x.kind == Tokens.RPAREN &&
    x.kind == Tokens.LBRACE &&
    x.kind == Tokens.RBRACE &&
    x.kind == Tokens.LSQUARE &&
    x.kind == Tokens.RSQUARE
is_comma(x) = x isa PUNCTUATION && x.kind == Tokens.COMMA
is_lparen(x) = x isa PUNCTUATION && x.kind == Tokens.LPAREN
is_rparen(x) = x isa PUNCTUATION && x.kind == Tokens.RPAREN
is_lbrace(x) = x isa PUNCTUATION && x.kind == Tokens.LBRACE
is_rbrace(x) = x isa PUNCTUATION && x.kind == Tokens.RBRACE
is_lsquare(x) = x isa PUNCTUATION && x.kind == Tokens.LSQUARE
is_rsquare(x) = x isa PUNCTUATION && x.kind == Tokens.RSQUARE

# KEYWORD
is_if(x) = x isa KEYWORD && x.kind == Tokens.IF
is_module(x) = x isa KEYWORD && x.kind == Tokens.MODULE
is_import(x) = x isa KEYWORD && x.kind == Tokens.IMPORT
is_importall(x) = x isa KEYWORD && x.kind == Tokens.IMPORTALL

# Literals

is_lit_string(x) = x isa LITERAL && (x.kind == Tokens.STRING || x.kind == Tokens.TRIPLE_STRING)

is_valid_iterator(x) = false
is_valid_iterator(x::BinarySyntaxOpCall) = is_eq(x.op)
is_valid_iterator(x::BinaryOpCall) = is_in(x.op) || is_elof(x.op)

_arg_id(x) = x
_arg_id(x::IDENTIFIER) = x
_arg_id(x::EXPR{Quotenode}) = x.args[1]
_arg_id(x::EXPR{Curly}) = _arg_id(x.args[1])
_arg_id(x::EXPR{Kw}) = _arg_id(x.args[1])

function _arg_id(x::UnarySyntaxOpCall)
    if is_dddot(x.arg2)
        return _arg_id(x.arg1)
    else
        return x
    end
end

function _arg_id(x::BinarySyntaxOpCall)
    if is_decl(x.op)
        return _arg_id(x.arg1)
    else
        return x
    end
end
function _arg_id(x::WhereOpCall)
    return _arg_id(x.arg1)
end

function get_where_params(x, params = [])
    return params
end

function get_where_params(x::WhereOpCall, params = [])
    for i = 1:length(x.args)
        a = x.args[i]
        if !(a isa PUNCTUATION)
            param = rem_subtype(a)
            param = rem_curly(param)
            push!(params, str_value(param))
        end
    end
    return params
end

function get_curly_params(x, params = [])
    return params
end
function get_curly_params(x::EXPR{Curly}, params = [])
    for i = 2:length(x.args)
        a = x.args[i]
        if !(a isa PUNCTUATION)
            param = rem_subtype(a)
            push!(params, str_value(param))
        end
    end
    return params
end

function get_sig_params(x, params = [])
    get_where_params(x, params)
    if x isa WhereOpCall && x.arg1 isa WhereOpCall
        get_where_params(x.arg1, params)
    end
    x = rem_where(x)
    x = rem_call(x)
    get_curly_params(x, params)
    return params
end

function rem_subtype(x)
    if x isa BinarySyntaxOpCall && x.op isa OPERATOR && x.op.kind == Tokens.ISSUBTYPE
        return x.arg1
    else
        return x
    end
end

function rem_decl(x)
    if x isa BinarySyntaxOpCall && is_decl(x.op)
        return x.arg1
    else
        return x
    end
end

function rem_curly(x)
    if x isa EXPR{Curly}
        return x.args[1]
    else
        return x
    end
end

function rem_call(x)
    if x isa EXPR{Call}
        return x.args[1]
    else
        return x
    end
end

function rem_where(x)
    if x isa WhereOpCall
        return rem_where(x.arg1)
    else
        return x
    end
end

function rem_invis(x)
    if x isa EXPR{InvisBrackets}
        return x.args[2]
    else
        return x
    end
end

function rem_dddot(x)
    if x isa UnarySyntaxOpCall && is_dddot(x.arg2)
        return x.arg1
    else
        return x
    end
end

function rem_kw(x)
    if x isa EXPR{Kw}
        return x.args[1]
    else
        return x
    end
end

# Definitions
defines_function(x) = false
defines_function(x::EXPR{FunctionDef}) = true
function defines_function(x::BinarySyntaxOpCall)
    if is_eq(x.op)
        sig = x.arg1
        while true
            if sig isa EXPR{Call}
                return true
            elseif sig isa BinarySyntaxOpCall && is_decl(sig.op) || sig isa WhereOpCall
                sig = sig.arg1
            elseif sig isa UnaryOpCall || sig isa UnarySyntaxOpCall# && sig.arg isa BinarySyntaxOpCall && is_decl(sig.arg.op) && sig.arg.arg1 isa EXPR{InvisBrackets}
                return true
            elseif sig isa BinaryOpCall || (sig isa BinarySyntaxOpCall && !(sig.op.kind == Tokens.DOT))
                return true
            else
                return false
            end
        end
    end
    return false
end

defines_macro(x) = false
defines_macro(x::EXPR{Macro}) = true

function defines_datatype(x)
    defines_struct(x) || defines_abstract(x) || defines_primitive(x)
end

defines_struct(x) = x isa EXPR{Struct} || defines_mutable(x)
defines_mutable(x) = x isa EXPR{Mutable}
defines_abstract(x) = x isa EXPR{Abstract}
function defines_primitive(x)
    x isa EXPR{Primitive}
end

defines_module(x) = x isa EXPR{ModuleH} || x isa EXPR{BareModule}

defines_anon_function(x) = x isa BinarySyntaxOpCall && is_anon_func(x.op)

function has_sig(x)
    defines_datatype(x) || defines_function(x) || defines_macro(x) || defines_anon_function(x)
end

"""
    get_sig(x)

Returns the full signature of function, macro and datatype definitions.
Should only be called when has_sig(x) == true.
"""
get_sig(x::EXPR{Struct}) = x.args[2]
get_sig(x::EXPR{Mutable}) = x.args[3]
get_sig(x::EXPR{Abstract}) = length(x.args) == 4 ? x.args[3] : x.args[2]
get_sig(x::EXPR{Primitive}) = x.args[3]
get_sig(x::EXPR{FunctionDef}) = x.args[2]
get_sig(x::EXPR{Macro}) = x.args[2]
function get_sig(x::BinarySyntaxOpCall)
    return x.arg1
end

function get_name(x::EXPR{T}) where T <: Union{Struct,Mutable,Abstract,Primitive}
    sig = get_sig(x)
    sig = rem_where(sig)
    sig = rem_subtype(sig)
    sig = rem_curly(sig)
    return sig
end

function get_name(x::EXPR{T}) where T <: Union{ModuleH,BareModule}
    x.args[2]
end

function get_name(sig)
    sig = rem_where(sig)
    sig = rem_decl(sig)
    sig = rem_call(sig)
    sig = rem_curly(sig)
    sig = rem_invis(sig)
    return sig
end

function get_name(x::EXPR{FunctionDef})
    sig = get_sig(x)
    sig = rem_where(sig)
    sig = rem_decl(sig)
    sig = rem_call(sig)
    sig = rem_curly(sig)
    sig = rem_invis(sig)
    return sig
end

function get_name(x::EXPR{Macro})
    sig = get_sig(x)
    sig = rem_call(sig)
    sig = rem_invis(sig)
    return sig
end

function get_name(x::BinarySyntaxOpCall)
    sig = x.arg1
    if sig isa UnaryOpCall
        return sig.op
    end
    sig = rem_where(sig)
    sig = rem_decl(sig)
    sig = rem_call(sig)
    sig = rem_curly(sig)
    sig = rem_invis(sig)
end

function get_args(x::IDENTIFIER)
    return []
end

function get_args(x)
    if defines_anon_function(x) && !(x.arg1 isa EXPR{TupleH})
        arg = x.arg1
        arg = rem_invis(arg)
        arg = get_arg_name(arg)
        return [arg]
    end
    sig = get_sig(x)
    sig = rem_where(sig)
    sig = rem_decl(sig)
    return get_args(sig)
end

function get_args(sig::EXPR{TupleH})
    args = []
    for i = 2:length(sig.args)
        arg = sig.args[i]
        arg isa PUNCTUATION && continue
        arg isa EXPR{Parameters} && continue
        arg_name = get_arg_name(arg)
        push!(args, arg_name)
    end
    return args
end

function get_args(x::EXPR{Do})
    args = []
    for i = 1:length(x.args[3].args)
        arg = x.args[3].args[i]
        arg isa PUNCTUATION && continue
        arg isa EXPR{Parameters} && continue
        arg_name = get_arg_name(arg)
        push!(args, arg_name)
    end
    return args
end

function get_args(sig::EXPR{Call})
    args = []
    sig = rem_where(sig)
    sig = rem_decl(sig)
    if sig isa EXPR{Call}
        for i = 2:length(sig.args)
            arg = sig.args[i]
            arg isa PUNCTUATION && continue
            if arg isa EXPR{Parameters}
                for j = 1:length(arg.args)
                    parg = arg.args[j]
                    parg isa PUNCTUATION && continue
                    parg_name = get_arg_name(parg)
                    push!(args, parg_name)
                end
            else
                arg_name = get_arg_name(arg)
                push!(args, arg_name)
            end
        end
    else
        error("not sig: $sig")
    end
    return args
end

function get_args(x::EXPR{Struct})
    args = []
    for arg in x.args[3]
        if !defines_function(arg)
            arg = rem_decl(arg)
            push!(args, arg)
        end
    end
    return args
end

function get_args(x::EXPR{Mutable})
    args = []
    for arg in x.args[4]
        if !defines_function(arg)
            arg = rem_decl(arg)
            push!(args, arg)
        end
    end
    return args
end

function get_args(x::EXPR{Flatten})
    get_args(x.args[1])
end

function get_args(x::EXPR{T}) where T <: Union{Generator,Filter}
    args = []
    if x.args[1] isa EXPR{Flatten} || x.args[1] isa EXPR{Generator}
        append!(args, get_args(x.args[1]))
    end

    if x.args[3] isa EXPR{Filter}
        return get_args(x.args[3])
    else
        for i = 3:length(x.args)
            arg = x.args[i]
            if is_valid_iterator(arg)
                arg = rem_decl(arg.arg1)
                arg = flatten_tuple(arg)
                arg = rem_decl.(arg)
                append!(args, arg)
            end
        end
        return args
    end
end

function get_arg_name(arg)
    arg = rem_kw(arg)
    arg = rem_dddot(arg)
    arg = rem_where(arg)
    arg = rem_decl(arg)
    arg = rem_subtype(arg)
    arg = rem_curly(arg)
    arg = rem_invis(arg)
end

function get_arg_type(arg)
    if arg isa BinarySyntaxOpCall && is_decl(arg.op)
        return Expr(arg.arg2)
    else
        return :Any
    end
end

get_body(x::EXPR{ModuleH}) = x.args[3]
get_body(x::EXPR{BareModule}) = x.args[3]
get_body(x::EXPR{If}) = x.args[3]
get_body(x::EXPR{For}) = x.args[3]
get_body(x::EXPR{While}) = x.args[3]
get_body(x::EXPR{FunctionDef}) = x.args[3]
get_body(x::EXPR{Macro}) = x.args[3]
get_body(x::EXPR{Struct}) = x.args[3]
get_body(x::EXPR{Mutable}) = x.args[4]

flatten_tuple(x::EXPR{InvisBrackets}, out = []) = flatten_tuple(x.args[2], out)
function flatten_tuple(x, out = [])
    if x isa EXPR{TupleH}
        for arg in x
            arg isa PUNCTUATION && continue
            flatten_tuple(arg, out)
        end
    else
        push!(out, x)
    end
    return out
end

"""
    get_id(x)

Get the IDENTIFIER name of a variable, possibly in the presence of
type declaration operators.
"""
function get_id(x::BinarySyntaxOpCall)
    if is_issubt(x.op) || is_decl(x.op)
        return get_id(x.arg1)
    else
        return x
    end
end

function get_id(x::WhereOpCall)
    return get_id(x.arg1)
end

function get_id(x::UnarySyntaxOpCall)
    if is_dddot(x.arg2)
        return get_id(x.arg1)
    else
        return x
    end
end

get_id(x::EXPR{Curly}) = get_id(x.args[1])
get_id(x::EXPR{InvisBrackets}) = get_id(x.args[2])
get_id(x) = x

"""
    get_t(x)

Basic inference in the presence of type declarations.
"""
get_t(x) = :Any
function get_t(x::BinarySyntaxOpCall)
    if is_decl(x.op)
        return Expr(x.arg2)
    else
        return :Any
    end
end

infer_t(x) = :Any
function infer_t(x::LITERAL)
    if x.kind == Tokens.INTEGER
        return :Int
    elseif x.kind == Tokens.FLOAT
        return :Float64
    elseif x.kind == Tokens.STRING
        return :String
    elseif x.kind == Tokens.TRIPLE_STRING
        return :String
    elseif x.kind == Tokens.CHAR
        return :Char
    elseif x.kind == Tokens.TRUE || x.kind == Tokens.FALSE
        return :Bool
    elseif x.kind == Tokens.CMD
        return :Cmd
    end
end

infer_t(x::EXPR{Vect}) = :(Array{Any,1})
infer_t(x::EXPR{Vcat}) = :(Array{Any,N})
infer_t(x::EXPR{TypedVcat}) = :(Array{$(Expr(x.args[1])),N})
infer_t(x::EXPR{Hcat}) = :(Array{Any,2})
infer_t(x::EXPR{TypedHcat}) = :(Array{$(Expr(x.args[1])),2})
infer_t(x::EXPR{Quote}) = :Expr
infer_t(x::EXPR{StringH}) = :String
infer_t(x::EXPR{Quotenode}) = :QuoteNode

"""
    contributes_scope(x)
Checks whether the body of `x` is included in the toplevel namespace.
"""
contributes_scope(x) = false
contributes_scope(x::EXPR{FileH}) = true
contributes_scope(x::EXPR{Begin}) = true
contributes_scope(x::EXPR{Block}) = true
contributes_scope(x::EXPR{Const}) = true
contributes_scope(x::EXPR{Global}) = true
contributes_scope(x::EXPR{Local}) = true
contributes_scope(x::EXPR{If}) = true
contributes_scope(x::EXPR{MacroCall}) = true
contributes_scope(x::EXPR{TopLevel}) = true


text:
# __precompile__(true)

module Plotly
using Compat
using Compat: String
using Requests
using JSON
using Reexport: @reexport
import Requests: URI, post

@reexport using PlotlyJS
export post
export set_credentials_file, RemotePlot, download_plot

include("utils.jl")

const api_version = "v2"

const default_kwargs = Dict{Symbol,Any}(:filename=>"Plot from Julia API",
                                        :world_readable=> true)

const default_opts = Dict{Symbol,Any}(:origin => "plot",
                                      :platform => "Julia",
                                      :version => "0.2")

## Taken from https://github.com/johnmyleswhite/Vega.jl/blob/master/src/Vega.jl#L51
# Open a URL in a browser
function openurl(url::String)
    if is_apple()
        run(`open $url`)
    elseif is_windows()
        run(`start $url`)
    elseif is_unix()
        run(`xdg-open $url`)
    end
end

openurl(url::URI) = openurl(string(url))

get_plot_endpoint() = "$(get_config().plotly_domain)/clientresp"

"""
Proxy for a plot stored on the Plotly cloud.
"""
immutable RemotePlot
    url::URI
end
RemotePlot(url::String) = RemotePlot(URI(url))

"""
Display a plot stored in the Plotly cloud in a browser window.
"""
Base.open(p::RemotePlot) = openurl(p.url)

"""
Post a local Plotly plot to the Plotly cloud.

Must be signed in first.
"""
function post(p::Plot; kwargs...)
    creds = get_credentials()
    endpoint = get_plot_endpoint()
    opt = merge(default_kwargs, Dict(:layout => p.layout.fields),
    Dict(kwargs))

    data = merge(default_opts,
    Dict("un" => creds.username,
    "key" => creds.api_key,
    "args" => json(p.data),
    "kwargs" => json(opt)))

    r = post(endpoint, data=data)
    body = parse_response(r)
    return RemotePlot(URI(body["url"]))
end

function Requests.post(l::AbstractLayout, meta_opts=Dict(); meta_kwargs...)
    creds = get_credentials()
    endpoint = get_plot_endpoint()

    meta = merge(meta_opts,
    get_required_params(["filename", "fileopt"], meta_opts),
    Dict(meta_kwargs))
    data = merge(default_opts,
    Dict("un" => creds.username,
         "key" => creds.api_key,
         "args" => json(l),
         "origin" => "layout",
         "kwargs" => json(meta)))

    parse_response(post(endpoint, data=data))
end

post(p::PlotlyJS.SyncPlot; kwargs...) = post(p.plot; kwargs...)

function style(style_opts, meta_opts=Dict(); meta_kwargs...)
    creds = get_credentials()
    endpoint = get_plot_endpoint()

    meta = merge(meta_opts,
    get_required_params(["filename", "fileopt"], meta_opts),
    Dict(meta_kwargs))
    data = merge(default_opts,
    Dict("un" => creds.username,
    "key" => creds.api_key,
    "args" => json([style_opts]),
    "origin" => "style",
    "kwargs" => json(meta_opts)))

    parse_response(post(endpoint, data=data))
end

"""
Transport a plot from the Plotly cloud to a local `Plot` object.

Must be signed in first if the plot is not public.
"""
function Base.download(plot::RemotePlot)
    creds = get_credentials()
    username = creds.username
    api_key = creds.api_key
    lib_version = string(default_opts[:platform], " ", default_opts[:version])
    auth = string("Basic ", base64encode("$username:$api_key"))
    options = Dict("Authorization"=>auth, "Plotly-Client-Platform"=>lib_version)
    original_path = plot.url.path
    if original_path[end] == '/'
        path = original_path[1:end-1]
    else
        path = original_path
    end
    endpoint = URI(plot.url, path="$path.json")
    response = get(endpoint, headers=options)
    local_plot = JSON.parse(Plot, bytestring(response))
    return PlotlyJS.SyncPlot(local_plot)
end

download_plot(url) = download(RemotePlot(url))
download_plot(plot::RemotePlot) = download(plot)

immutable PlotlyError <: Exception
    msg::String
end

function Base.show(io::IO, err::PlotlyError)
    print(io, "Plotly error: $(err.msg)")
end

function parse_response(r)
    body = Requests.json(r)
    if statuscode(r) ≠ 200
        throw(PlotlyError("Non-sucessful status code: $(statuscode(r))"))
    elseif "error" ∈ keys(body) && body["error"] ≠ ""
        throw(PlotlyError(body["error"]))
    elseif "detail" ∈ keys(body) && body["detail"] ≠ ""
        throw(PlotlyError(body["detail"]))
    else
        body
    end
end

end


text:
using Test
import ArchGDAL as AG

@testset "test_geos_operations.jl" begin
    function equivalent_to_wkt(geom::AG.Geometry, wkt::String)
        fromWKT(wkt) do test_geom
            @test toWKT(geom) == toWKT(test_geom)
        end
    end

    @testset "Interpolation along a LineString" begin
        AG.createlinestring([
            (8.0, 1.0),
            (9.0, 1.0),
            (9.0, 2.0),
            (8.0, 2.0),
        ]) do ls
            for (dist, dest) in [
                (1.0, (9, 1)),
                (2.0, (9, 2)),
                (1.5, (9.0, 1.5)),
                (2.5, (8.5, 2.0)),
            ]
                AG.pointalongline(ls, dist) do pt1
                    AG.createpoint(dest) do pt2
                        @test AG.toWKT(pt1) == AG.toWKT(pt2)
                    end
                end
                @test AG.toWKT(AG.pointalongline(ls, dist)) ==
                      AG.toWKT(AG.createpoint(dest))
            end
        end
    end

    @testset "Contains operation" begin
        AG.fromWKT("POLYGON EMPTY") do g1
            AG.fromWKT("POLYGON EMPTY") do g2
                @test AG.contains(g1, g2) == false
                @test AG.contains(g2, g1) == false
            end
        end

        AG.fromWKT("POLYGON((1 1,1 5,5 5,5 1,1 1))") do g1
            AG.fromWKT("POINT(2 2)") do g2
                @test AG.contains(g1, g2) == true
                @test AG.contains(g2, g1) == false
            end
        end

        AG.fromWKT("MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))") do g1
            AG.fromWKT("POLYGON((1 1,1 2,2 2,2 1,1 1))") do g2
                @test AG.contains(g1, g2) == true
                @test AG.contains(g2, g1) == false
            end
        end
    end

    @testset "Convex Hull" begin
        AG.fromWKT(
            "MULTIPOINT (130 240, 130 240, 130 240, 570 240, 570 240, 570 240, 650 240)",
        ) do input
            AG.fromWKT("LINESTRING (130 240, 650 240)") do expected
                AG.convexhull(input) do output
                    @test AG.isempty(output) == false
                    @test AG.toWKT(output) == AG.toWKT(expected)
                end
                @test AG.toWKT(AG.convexhull(input)) == AG.toWKT(expected)
            end
        end
    end

    @testset "Delaunay Triangulation" begin
        AG.fromWKT("POLYGON EMPTY") do g1
            AG.delaunaytriangulation(g1, 0, true) do g2
                @test AG.isempty(g1) == true
                @test AG.isempty(g2) == true
                @test AG.toWKT(g2) == "MULTILINESTRING EMPTY"
            end
            @test AG.toWKT(AG.delaunaytriangulation(g1, 0, true)) ==
                  "MULTILINESTRING EMPTY"
        end
        AG.fromWKT("POINT(0 0)") do g1
            AG.delaunaytriangulation(g1, 0, false) do g2
                @test AG.isempty(g2) == true
                @test AG.toWKT(g2) == "GEOMETRYCOLLECTION EMPTY"
            end
            @test AG.toWKT(AG.delaunaytriangulation(g1, 0, false)) ==
                  "GEOMETRYCOLLECTION EMPTY"
        end
        AG.fromWKT("MULTIPOINT(0 0, 5 0, 10 0)") do g1
            AG.delaunaytriangulation(g1, 0, false) do g2
                @test AG.toWKT(g2) == "GEOMETRYCOLLECTION EMPTY"
            end
            AG.delaunaytriangulation(g1, 0, true) do g2
                @test AG.toWKT(g2) == "MULTILINESTRING ((5 0,10 0),(0 0,5 0))"
            end
        end
        AG.fromWKT("MULTIPOINT(0 0, 10 0, 10 10, 11 10)") do g1
            AG.delaunaytriangulation(g1, 2.0, true) do g2
                @test AG.toWKT(g2) ==
                      "MULTILINESTRING ((0 0,10 10),(0 0,10 0),(10 0,10 10))"
            end
        end
    end

    @testset "Distance" begin
        AG.fromWKT("POINT(10 10)") do g1
            AG.fromWKT("POINT(3 6)") do g2
                @test AG.distance(g1, g2) ≈ 8.06225774829855 atol = 1e-12
            end
        end
    end

    function test_method(
        f::Function,
        wkt1::AbstractString,
        wkt2::AbstractString,
    )
        AG.fromWKT(wkt1) do geom
            f(geom) do result
                @test AG.toWKT(result) == wkt2
            end
            @test AG.toWKT(f(geom)) == wkt2
        end
    end

    @testset "Centroid" begin
        test_method(AG.centroid, "POINT(10 0)", "POINT (10 0)")
        test_method(AG.centroid, "LINESTRING(0 0, 10 0)", "POINT (5 0)")
        test_method(
            AG.centroid,
            "POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))",
            "POINT (5 5)",
        )
        test_method(AG.centroid, "LINESTRING EMPTY", "POINT EMPTY")
    end

    @testset "Point on Surface" begin
        test_method(AG.pointonsurface, "POINT(10 0)", "POINT (10 0)")
        test_method(
            AG.pointonsurface,
            "LINESTRING(0 0, 5 0, 10 0)",
            "POINT (5 0)",
        )
        test_method(
            AG.pointonsurface,
            "POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))",
            "POINT (5 5)",
        )
        test_method(AG.pointonsurface, "LINESTRING EMPTY", "POINT EMPTY")
        test_method(AG.pointonsurface, "LINESTRING(0 0, 0 0)", "POINT (0 0)")
    end

    function test_method(
        f::Function,
        wkt1::AbstractString,
        wkt2::AbstractString,
        wkt3::AbstractString,
    )
        AG.fromWKT(wkt1) do geom1
            AG.fromWKT(wkt2) do geom2
                f(geom1, geom2) do result
                    @test AG.toWKT(result) == wkt3
                end
                @test AG.toWKT(f(geom1, geom2)) == wkt3
            end
        end
    end

    @testset "Intersection" begin
        test_method(
            AG.intersection,
            "POLYGON EMPTY",
            "POLYGON EMPTY",
            "POLYGON EMPTY",
        )
        test_method(
            AG.intersection,
            "POLYGON((1 1,1 5,5 5,5 1,1 1))",
            "POINT(2 2)",
            "POINT (2 2)",
        )
        test_method(
            AG.intersection,
            "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))",
            "POLYGON((-1 1,-1 2,2 2,2 1,-1 1))",
            "POLYGON ((0 2,2 2,2 1,0 1,0 2))",
        )
        test_method(
            AG.intersection,
            "MULTIPOLYGON(((0 0,5 10,10 0,0 0),(1 1,1 2,2 2,2 1,1 1),(100 100,100 102,102 102,102 100,100 100)))",
            "POLYGON((0 1,0 2,10 2,10 1,0 1))",
            "GEOMETRYCOLLECTION (POLYGON ((1 2,1 1,0.5 1.0,1 2)),POLYGON ((9.5 1.0,2 1,2 2,9 2,9.5 1.0)),LINESTRING (1 2,2 2),LINESTRING (2 1,1 1))",
        )
    end

    function test_predicate(f::Function, wkt1, wkt2, result::Bool)
        AG.fromWKT(wkt1) do geom1
            AG.fromWKT(wkt2) do geom2
                @test f(geom1, geom2) == result
            end
        end
    end

    @testset "Intersects" begin
        test_predicate(AG.intersects, "POLYGON EMPTY", "POLYGON EMPTY", false)
        test_predicate(
            AG.intersects,
            "POLYGON((1 1,1 5,5 5,5 1,1 1))",
            "POINT(2 2)",
            true,
        )
        test_predicate(
            AG.intersects,
            "POINT(2 2)",
            "POLYGON((1 1,1 5,5 5,5 1,1 1))",
            true,
        )
        test_predicate(
            AG.intersects,
            "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))",
            "POLYGON((1 1,1 2,2 2,2 1,1 1))",
            true,
        )
        test_predicate(
            AG.intersects,
            "POLYGON((1 1,1 2,2 2,2 1,1 1))",
            "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))",
            true,
        )
    end

    @testset "Within" begin
        test_predicate(AG.within, "POLYGON EMPTY", "POLYGON EMPTY", false)
        test_predicate(
            AG.within,
            "POLYGON((1 1,1 5,5 5,5 1,1 1))",
            "POINT(2 2)",
            false,
        )
        test_predicate(
            AG.within,
            "POINT(2 2)",
            "POLYGON((1 1,1 5,5 5,5 1,1 1))",
            true,
        )
        test_predicate(
            AG.within,
            "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))",
            "POLYGON((1 1,1 2,2 2,2 1,1 1))",
            false,
        )
        test_predicate(
            AG.within,
            "POLYGON((1 1,1 2,2 2,2 1,1 1))",
            "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))",
            true,
        )
    end

    @testset "Simplify" begin
        AG.fromWKT(
            "POLYGON((56.528666666700 25.2101666667, 56.529000000000 25.2105000000, 56.528833333300 25.2103333333, 56.528666666700 25.2101666667))",
        ) do g1
            AG.simplify(g1, 0.0) do g2
                @test AG.toWKT(g2) == "POLYGON EMPTY"
            end
            @test AG.toWKT(AG.simplify(g1, 0.0)) == "POLYGON EMPTY"
        end

        AG.fromWKT(
            "POLYGON((56.528666666700 25.2101666667, 56.529000000000 25.2105000000, 56.528833333300 25.2103333333, 56.528666666700 25.2101666667))",
        ) do g1
            AG.simplifypreservetopology(g1, 43.2) do g2
                @test AG.toWKT(g2) ==
                      "POLYGON ((56.5286666667 25.2101666667,56.529 25.2105,56.5288333333 25.2103333333,56.5286666667 25.2101666667))"
            end
            @test AG.toWKT(AG.simplifypreservetopology(g1, 43.2)) ==
                  "POLYGON ((56.5286666667 25.2101666667,56.529 25.2105,56.5288333333 25.2103333333,56.5286666667 25.2101666667))"
        end
    end
end


text:
import Pkg
Pkg.activate(normpath(@__DIR__))

import ProjectMoana: HdrImage, read_pfm_image, normalize_image, clamp_image, greet, tone_mapping
import Images: save
using ArgParse

function parse_commandline()
    s = ArgParseSettings(description = "This program converts a PFM image into a PNG image. Try me!",
                               usage = "usage: [--help] [--file_in FILE_IN] [--clamp CLAMP]  [--a A] [--gamma Γ] [--file_out FILE_OUT]",
                              epilog = "Let's try again!")

    @add_arg_table s begin
        "--file_in"
            help = "input PFM file name;"
            required = true
        "--clamp"
            help = "specify which clamping to use: Custom (C) or ImageMagick (IM);"
            default = "C"
        "--a"
            help = "a_factor;"
            required = false
            default = 1.
            arg_type = Float64
        "--gamma"
            help = "γ factor, for monitor response;"
            required = false
            default = 1.
            arg_type = Float64
        "--file_out"
            help = "output LDR file name (with extension)."
            required = false
            default = "LDR_out.png"
    end

    return parse_args(s)
end

function main()
    # greet()
    params = parse_commandline()

# firtsly, open the input file
    img = HdrImage(1, 1)
    open(params["file_in"], "r") do inpf
        img = read_pfm_image(inpf)
    end

    println("File $(params["file_in"]) has been read from disk.") # check

# then normalizing and clamping

    if params["clamp"] == "C"
        normalize_image(img, params["a"])
        clamp_image(img)
        tone_mapping(img, params["gamma"])
    end

# saving the image in the output format using Images method

    image = reshape(img.pixels, (img.width,img.height))

    save("$(params["file_out"])",image')

    println("File $(params["file_out"]) has been written to disk.") # check

end

main()


text:
"""
    UpwindDifferences

Implements a finite differences method based on an upwind scheme to solve Hamilton-Jacobi-Bellman (HJB) equations.
"""
module UpwindDifferences

#=
Problem:
    ρv(x) = max_c {r(x, c) + ∂vₓ(x)ẋ(x, c)}
    where:  ẋ(x, c) = g(x, c)
            ẋ(x̲, c) ≥ 0
            ẋ(x̅, c) ≤ 0

Let b(x, ∂v) = argmax_c {r(x, c) + ∂vₓ(x)g(x, c)}

This problem can be efficient solved without the use of a cache.
However the cached implementation allows for vectorization.

Inputs:
    1. x -> a one-dimensional state grid, assumed to be evenly spaced, i.e. a LinRange
    2. n -> x.len

Edge cases:
    - Non standard indexing? Enforce LinRange
    - Small n?
=#

"""Contains the optimal reward vector, `R`, and forward and backward drift, `GF` and `GB`."""
struct UpwindResult{T <: Number}
    R::Vector{T}
    GF::Vector{T}
    GB::Vector{T}
end

"""Compute the forward difference of vector `v` at index i"""
function diff(v, dx, i)
    return (v[i+1] - v[i]) / dx
end

"""
    upwind!(UR::UpwindResult, v, xs, reward, drift, policy, zerodrift)

Update `UR` using upwind finite differences based on discretized value function `v`.

See [`upwind`](@ref) for further details.
`UR` contains 3 vectors: (R) the optimal reward; (GF) the optimal drift when non-negative; (GB) the optimal drift when non-positive.
"""
function upwind!(UR::UpwindResult, v, xs, reward, drift, policy, zerodrift)

    R, GF, GB = UR.R, UR.GF, UR.GB
    dx = step(xs)
    n  = length(xs)
    z  = zero(eltype(GB))

    # Lower state constraint
    dvf = diff(v, dx, 1)
    bf = policy(xs[1], dvf)
    gf = drift(xs[1], bf)
    if gf > z
        R[1] = reward(xs[1], bf)
        GF[1] = gf
    else
        R[1] = reward(xs[1], zerodrift(xs[1]))
        GF[1] = z
    end
    GB[1] = z

    # State interior
    for i in 2:n-1
        dvb = dvf
        dvf = diff(v, dx, i)
        R[i], GF[i], GB[i] = get_interior_upwind(xs[i], dvf, dvb, reward, drift, policy, zerodrift)
    end

    # Upper state constraint
    dvb = diff(v, dx, n-1)
    bb = policy(xs[end], dvb)
    gb = drift(xs[end], bb)
    if gb < z
        R[end] = reward(xs[end], bb)
        GB[end] = gb
    else
        R[end] = reward(xs[end], zerodrift(xs[end]))
        GB[end] = z
    end
    GF[end] = z

    return UR
end

"""
    upwind(v::Vector, xs::LinRange, reward, drift, policy, zerodrift)

Compute the optimal reward, and forward and backward drift, based on value function `v` using an upwind differences scheme.

This function can be used when an upwind finite differences approximation of the value function derivative is required to compute an optimal control solution.
As the resulting drift is required to compute the upwind approximation, it is also returned.

See also [`upwind!`](@ref) for an efficient inplace version.

# Arguments
- `v::Vector`: the discretized value function from which the result will be computed.
- `xs::LinRange`: the (one-dimensional) discretized state-space at which `v` is evaluated.
- `reward(x, c)::Function`: returns the flow reward given state value `x` and control value `c`.
- `drift(x, c)::Function`: returns the drift in the state given current state value `x` and control value `v`.
- `policy(x, dv)::Function`: returns the optimal control value given current state value `x` and value function derivative `dv`.
- zerodrift(x)::Function`: returns the control value which results in zero drift, i.e, `drift(x, zerodrift(x)) = 0`.
"""
function upwind(v, xs, reward, drift, policy, zerodrift)
    UR = UpwindResult(similar(v), similar(v), similar(v))
    upwind!(UR, v, xs, reward, drift, policy, zerodrift)
    return UR
end

function get_interior_upwind(x, dvf, dvb, reward, drift, policy, zerodrift)
    # Kernel of upwinding algorithm
    bf = policy(x, dvf)
    bb = policy(x, dvb)
    gf = drift(x, bf)
    gb = drift(x, bb)
    z = zero(typeof(gf))

    if gf > 0 && gb ≥ 0
        R = reward(x, bf)
        GF = gf
        GB = z
    elseif gb < 0 && gf ≤ 0
        R = reward(x, bb)
        GF = z
        GB = gb
    elseif gf ≤ 0 && gb ≥ 0
        R = reward(x, zerodrift(x))
        GF = z
        GB = z
    else
        rf = reward(x, bf)
        Hf = rf + dvf * gf
        rb = reward(x, bb)
        Hb = rb + dvb * gb
        H0 = reward(x, zerodrift(x))

        if Hf > max(Hb, H0)
            R = rf
            GF = gf
            GB = z
        elseif Hb > H0
            R = rb
            GF = z
            GB = gb
        else
            R = H0
            GF = z
            GB = z
        end
    end
    return R, GF, GB
end

export UpwindResult, upwind!, upwind

end # module