text:
################################################################
"""
````
integrate(system,F,U; boundary=false)
````

Integrate node function (same signature as reaction or storage)
 `F` of  solution vector region-wise over domain or boundary.
The result is  `nspec x nregion` vector.
"""
function integrate(system::AbstractSystem{Tv,Ti,Tm},F::Function,U::AbstractMatrix{Tu}; boundary=false) where {Tu,Tv,Ti,Tm}
    grid=system.grid
    data=system.physics.data
    nspecies=num_species(system)
    res=zeros(Tu,nspecies)
    node=Node{Tv,Ti}(system)
    nodeparams=(node,)
    bnode=BNode{Tv,Ti}(system)
    nodeparams=(bnode,)
    if isdata(data)
        nodeparams=(node,data,)
    end

    if boundary

        geom=grid[BFaceGeometries][1]
        bfaceregions=grid[BFaceRegions]
        nbfaceregions=maximum(bfaceregions)
        integral=zeros(Tu,nspecies,nbfaceregions)

        for ibface=1:num_bfaces(grid)
            for inode=1:num_nodes(geom)
                _fill!(bnode,inode,ibface)
                res.=zero(Tv)
                @views F(res,U[:,bnode.index],nodeparams...)
                for ispec=1:nspecies
                    if system.node_dof[ispec,bnode.index]==ispec
                        integral[ispec,bnode.region]+=system.bfacenodefactors[inode,ibface]*res[ispec]
                    end
                end
            end
        end
    else
        geom=grid[CellGeometries][1]
        cellnodes=grid[CellNodes]
        cellregions=grid[CellRegions]
        ncellregions=maximum(cellregions)
        integral=zeros(Tu,nspecies,ncellregions)

        for icell=1:num_cells(grid)
            for inode=1:num_nodes(geom)
                _fill!(node,inode,icell)
                res.=zero(Tv)
                @views F(res,U[:,node.index],nodeparams...)
                for ispec=1:nspecies
                    if system.node_dof[ispec,node.index]==ispec
                        integral[ispec,node.region]+=system.cellnodefactors[inode,icell]*res[ispec]
                    end
                end
            end
        end
    end

    return integral
end

############################################################################
"""
$(SIGNATURES)

Reconstruction of  edge flux as  vector function  on the nodes  of the
triangulation.  The result  can be seen as a  piecewiesw linear vector
function in the FEM space spanned by the discretization nodes exhibiting
the flux density.

The reconstruction is based on the  "magic formula"
R. Eymard, T. Gallouet, R. Herbin, IMA Journal of Numerical Analysis (2006)
26, 326−353, Lemma 2.4 .

The return value is a `dim x nspec x nnodes` vector. The flux of species i
can  e.g. plotted via GridVisualize.vectorplot.

Example:
```julia
    ispec=3
    vis=GridVisualizer(Plotter=Plotter)
    scalarplot!(vis,grid,solution[ispec,:],clear=true,colormap=:summer)
    vectorplot!(vis,grid,nf[:,ispec,:],clear=false)
    reveal(vis)
```

CAVEAT: there is a possible unsolved problem with the values at domain
corners in the code. Please see any potential boundary artifacts as a manifestation
of this issue and report them.
"""
function nodeflux(system::AbstractSystem{Tv,Ti},U::AbstractArray{Tu,2}) where {Tu,Tv,Ti}
    grid=system.grid
    dim=dim_space(grid)
    nnodes=num_nodes(grid)
    nspecies=num_species(system)
    nodeflux=zeros(Tu,dim,nspecies,nnodes)
    edgeflux=zeros(Tu,nspecies)
    xsigma=grid[VoronoiFaceCenters]
    coord=grid[Coordinates]
    nodevol=zeros(Tv,nnodes)
    cellnodes=grid[CellNodes]

    physics=system.physics
    node=Node{Tv,Ti}(system)
    bnode=BNode{Tv,Ti}(system)
    edge=Edge{Tv,Ti}(system)
    bedge=BEdge{Tv,Ti}(system)
    @create_physics_wrappers(physics,node,bnode,edge,bedge)

    UKL=Array{Tu,1}(undef,2*nspecies)
    geom=grid[CellGeometries][1]

    for icell=1:num_cells(grid)
        for iedge=1:num_edges(geom)
            _fill!(edge,iedge,icell)
            K=edge.node[1]
            L=edge.node[2]
            fac=system.celledgefactors[iedge,icell]
            @views UKL[1:nspecies].=U[:,edge.node[1]]
            @views UKL[nspecies+1:2*nspecies].=U[:,edge.node[2]]
            edgeflux.=zero(Tv)
            fluxwrap(edgeflux,UKL)
            for ispec=1:nspecies
                if system.node_dof[ispec,K]==ispec && system.node_dof[ispec,L]==ispec
                    @views nodeflux[:,ispec,K].+=fac*edgeflux[ispec]*(xsigma[:,edge.index]-coord[:,K])
                    @views nodeflux[:,ispec,L].-=fac*edgeflux[ispec]*(xsigma[:,edge.index]-coord[:,L])
                end
            end
        end

        for inode=1:num_nodes(geom)
            nodevol[cellnodes[inode,icell]]+=system.cellnodefactors[inode,icell]
        end
    end
    for inode=1:nnodes
        @views nodeflux[:,:,inode]/=nodevol[inode]
    end
    nodeflux
end


text:
# A symmetry operation (SymOp) is a couple (W, w) of a
# unitary (in cartesian coordinates, but not in reduced coordinates)
# matrix W and a translation w such that, for each atom of
# type A at position a, W a + w is also an atom of type A.
# This induces an operator
# (Uu)(x) = u(W x + w)
# or in Fourier space
# (Uu)(G) = e^{-i G τ} u(S^-1 G)
# with
# S = W'
# τ = -W^-1 w
# (all these formulas are valid both in reduced and cartesian coordinates)

# Time-reversal symmetries are the anti-unitaries
# (Uu)(x) = conj(u(Wx+w))
# or in Fourier space
# (Uu)(G) = e^{i G τ} conj(u(-S^-1 G))

# Tolerance to consider two atomic positions as equal (in relative coordinates)
const SYMMETRY_TOLERANCE = 1e-5

is_approx_integer(r; tol=SYMMETRY_TOLERANCE) = all(ri -> abs(ri - round(ri)) ≤ tol, r)

struct SymOp{T <: Real}
    # (Uu)(x) = u(W x + w) in real space
    W::Mat3{Int}
    w::Vec3{T}

    # (Uu)(G) = e^{-i G τ} u(S^-1 G) in reciprocal space
    S::Mat3{Int}
    τ::Vec3{T}

    function SymOp(W, w)
        w = mod.(w, 1)
        S = W'
        τ = -W \w
        new{eltype(τ)}(W, w, S, τ)
    end
end

Base.:(==)(op1::SymOp, op2::SymOp) = op1.W == op2.W && op1.w == op2.w
function Base.isapprox(op1::SymOp, op2::SymOp; atol=SYMMETRY_TOLERANCE)
    op1.W == op2.W && is_approx_integer(op1.w - op2.w; tol=atol)
end
Base.one(::Type{SymOp}) = SymOp(Mat3{Int}(I), Vec3(zeros(Bool, 3)))
Base.one(::SymOp) = one(SymOp)

# group composition and inverse.
function Base.:*(op1::SymOp, op2::SymOp)
    W = op1.W * op2.W
    w = op1.w + op1.W * op2.w
    SymOp(W, w)
end
Base.inv(op::SymOp) = SymOp(inv(op.W), -op.W\op.w)

function check_group(symops::Vector; kwargs...)
    is_approx_in_symops(s1) = any(s -> isapprox(s, s1; kwargs...), symops)
    is_approx_in_symops(one(SymOp)) || error("check_group: no identity element")
    for s in symops
        if !is_approx_in_symops(inv(s))
            error("check_group: symop $s with inverse $(inv(s)) is not in the group")
        end
        for s2 in symops
            if !is_approx_in_symops(s*s2) || !is_approx_in_symops(s2*s)
                error("check_group: product is not stable")
            end
        end
    end
    symops
end

ext: jl


text:
# CXXRecordDecl
function getCanonicalDecl(x::AbstractCXXRecordDecl)
    @check_ptrs x
    return CXCXXRecordDecl(clang_CXXRecordDecl_getCanonicalDecl(x))
end

function getPreviousDecl(x::AbstractCXXRecordDecl)
    @check_ptrs x
    return CXCXXRecordDecl(clang_CXXRecordDecl_getPreviousDecl(x))
end

function getMostRecentDecl(x::AbstractCXXRecordDecl)
    @check_ptrs x
    return CXCXXRecordDecl(clang_CXXRecordDecl_getMostRecentDecl(x))
end

function getMostRecentNonInjectedDecl(x::AbstractCXXRecordDecl)
    @check_ptrs x
    return CXCXXRecordDecl(clang_CXXRecordDecl_getMostRecentNonInjectedDecl(x))
end

function getDefinition(x::AbstractCXXRecordDecl)
    @check_ptrs x
    return CXCXXRecordDecl(clang_CXXRecordDecl_getDefinition(x))
end

function hasDefinition(x::AbstractCXXRecordDecl)
    @check_ptrs x
    return clang_CXXRecordDecl_hasDefinition(x)
end

function isLambda(x::AbstractCXXRecordDecl)
    @check_ptrs x
    return clang_CXXRecordDecl_isLambda(x)
end

function isGenericLambda(x::AbstractCXXRecordDecl)
    @check_ptrs x
    return clang_CXXRecordDecl_isGenericLambda(x)
end

function getGenericLambdaTemplateParameterList(x::AbstractCXXRecordDecl)
    @check_ptrs x
    return TemplateParameterList(clang_CXXRecordDecl_getGenericLambdaTemplateParameterList(x))
end

function isAggregate(x::AbstractCXXRecordDecl)
    @check_ptrs x
    return clang_CXXRecordDecl_isAggregate(x)
end

function isPOD(x::AbstractCXXRecordDecl)
    @check_ptrs x
    return clang_CXXRecordDecl_isPOD(x)
end

function isCLike(x::AbstractCXXRecordDecl)
    @check_ptrs x
    return clang_CXXRecordDecl_isCLike(x)
end

function isEmpty(x::AbstractCXXRecordDecl)
    @check_ptrs x
    return clang_CXXRecordDecl_isEmpty(x)
end


text:
@testset "Check pressure calculation" begin

    function compute_ẋ!(ẋ, x, t)
        # Zero the velocity
        reset_velocity!(ẋ)
        # Compute the self-induced velocity of the system
        self_induce_velocity!(ẋ, x, t)
    end

    atol_AD  = 100*eps()
    atol  = 1e-6

    Δt = sqrt(eps())

    δ = 1e-4

    Δs = 0.1
    ss = -3.0:Δs:3.0
    sensors = complex.(ss)
    Ny = length(ss)

    # Test with nothing
    testsys = vcat(Vortex.Blob(im, 0.0, δ),
                   Vortex.Blob(-im, 0.0, δ))

    press = pressure(sensors, deepcopy(testsys), 0.0)
    # Expected value 0.0 everywhere
    @test isapprox(press, zeros(Ny), atol = atol)

    # Test with one vortex (parallel translation)
    testsys = vcat(Vortex.Blob(im, 1.0, δ),
         Vortex.Blob(-im, -1.0, δ))
    press_FD = pressure_FD(sensors, deepcopy(testsys), 0.0, Δt)

    testϕ₋ = real.(complexpotential(sensors, testsys))
    testẋ = allocate_velocity(testsys)
    compute_ẋ!(testẋ, testsys, 0.0)
    testsys₊ = deepcopy(testsys)
    advect!(testsys₊, testsys₊, testẋ, Δt)
    testϕ₊ = real.(complexpotential(sensors, testsys₊))
    @test isapprox(press_FD, -0.5*abs2.(induce_velocity(sensors, testsys , 0.0)) - (testϕ₊-testϕ₋)/Δt,
                   atol = atol)

    # Test with vortices and sources
    testsys = (vcat(Vortex.Blob(im, 1.0, δ),
                    Vortex.Blob(-im, -1.0, δ)),
               vcat(Source.Blob(2*im,  1.0, δ),
                    Source.Blob(-2*im,  1.0, δ)))
    press_FD = pressure_FD(sensors, deepcopy(testsys), 0.0, Δt)

    testϕ₋ = real.(complexpotential(sensors, (testsys[1:2])))
    testẋ = allocate_velocity(testsys)
    compute_ẋ!(testẋ, testsys, 0.0)
    testsys₊ = deepcopy(testsys)
    advect!(testsys₊, testsys₊, testẋ, Δt)
    testϕ₊ = real.(complexpotential(sensors, (testsys₊[1:2])))
    @test isapprox(press_FD, (-0.5*abs2.(induce_velocity(sensors, testsys, 0.0)) - (testϕ₊-testϕ₋)/Δt), atol = atol)

    # Test with vortices and sources and freestream
    testsys = (vcat(Vortex.Blob(im, 1.0, δ),
               Vortex.Blob(-im, -1.0, δ)),
               vcat(Source.Blob(2*im,  1.0, δ),
               Source.Blob(-2*im,  1.0, δ)))
    press_FD = pressure_FD(sensors, deepcopy(testsys), 0.0, Δt)

    testϕ₋ = real.(complexpotential(sensors, (testsys[1:2])))
    testẋ = allocate_velocity(testsys)
    compute_ẋ!(testẋ, testsys, 0.0)
    testsys₊ = deepcopy(testsys)
    advect!(testsys₊, testsys₊, testẋ, Δt)
    testϕ₊ = real.(complexpotential(sensors, (testsys₊[1:2])))
    @test isapprox(press_FD, (-0.5*abs2.(induce_velocity(sensors, testsys, 0.0)) - (testϕ₊-testϕ₋)/Δt), atol = atol)

    # Frozen advected pressure field by a pair of two vortices
    testsys = (vcat(Vortex.Blob(im, 1.0, 0.0),
           Vortex.Blob(-im, -1.0, 0.0)),
           vcat(Source.Blob(2*im,  0.0, 0.0),
           Source.Blob(-2*im,  0.0, 0.0)))

    press_analytical(x) = -0.5*(1/(π*(x^2+1.0)))^2 + 1/(4*π^2*(x^2 + 1.0))
    press = pressure(complex(-10.0:0.01:10.0), deepcopy(testsys), 0.0)
    press_FD = pressure_FD(complex(-10.0:0.01:10.0), deepcopy(testsys), 0.0, Δt);

    @test isapprox(press, press_analytical.(-10.0:0.01:10.0), atol = atol_AD)
    @test isapprox(press_FD, press_analytical.(-10.0:0.01:10.0), atol = atol)
end


text:
function faceconnectivity(p,t,f,f2t,elemtype,prdexpr,porder)

dim = size(p,1);
nprd = size(prdexpr,1);
nf = size(f2t,2);

philocvl,philocfc,~,~,perm = localbasis(porder,dim,elemtype);
philocvl = philocvl';
philocfc = philocfc';
face = getelemface(dim,elemtype);

npf = size(perm,1);
npe = size(philocvl,2);

facecon = zeros(Int,npf,2,nf);
for i = 1:nf
    e1 = f2t[1,i];
    l1 = f2t[2,i];
    e2 = f2t[3,i];
    l2 = f2t[4,i];

    f1 = t[face[:,l1],e1];
    pf = p[:,f1]*philocfc; # nodal points on face i
    pf = round.(pf,digits=8);

    t1 = t[:,e1];
    p1 = p[:,t1]*philocvl;
    p1 = p1[:,perm[:,l1]];
    p1 = round.(p1,digits=8);

    if f[l1,e1]<0 # periodic face
        pbnd = -f[l1,e1];
        for k=1:nprd
            if (prdexpr[k,1]==pbnd)
                pf = prdexpr[k,2](pf);
                p1 = prdexpr[k,2](p1);
                break;
            elseif (prdexpr[k,3]==pbnd)
                pf = prdexpr[k,4](pf);
                p1 = prdexpr[k,4](p1);
                break;
            end
        end
    end

    j1 = xiny(p1,pf,0);
    facecon[:,1,i] = (e1-1)*npe .+ perm[j1,l1];

    if e2>0 # face i is an interior face
        t2 = t[:,e2];
        p2 = p[:,t2]*philocvl;
        p2 = p2[:,perm[:,l2]];
        p2 = round.(p2,digits=8);
        if f[l2,e2]<0 # periodic face
            pbnd = -f[l2,e2];
            if (prdexpr[k,1]==pbnd)
                p2 = prdexpr[k,2](p2);
            elseif (prdexpr[k,3]==pbnd)
                p2 = prdexpr[k,4](p2);
            end
        end
        j2 = xiny(p2,pf,0);
        facecon[:,2,i] = (e2-1)*npe .+ perm[j2,l2];
    else
        facecon[:,2,i] = facecon[:,1,i];
    end
end

return facecon

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
# Note:
# This module is meant to be removed from ComoniconGUI when Comonicon itself is ready.

module Parse

using ExprTools
using Comonicon.Types
using Comonicon.CodeGen: prettify, pushmaybe!
using Comonicon.Parse
using Comonicon.Parse: xcall, checksum
using Comonicon.Parse: parse_args, parse_kwargs, parse_function
using Comonicon.Parse: create_casted_commands, casted_commands, set_cmd!

using ComoniconGUI
using ComoniconGUI.CodeGen

include("utils.jl")
include("cast.jl")

export @cast, @main,
        iscached,
        cachefile,
        create_cache,
        enable_cache,
        disable_cache

end


text:
# Copyright 2021 The CLVR Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# This file contains the Results struct's definition and its dispatches.

"""
Defines the progress of execution at each logging step.
* iterations: Number of iterations elapsed.
* times: Elapsed times since start of execution (in seconds).
* fvaluegaps: Primal and dual objective value gaps.
* metricLPs: The computed values of the LP metric.
"""
struct Results
    iterations::Vector{Float64}
    times::Vector{Float64}
    fvaluegaps::Vector{Float64}
    metricLPs::Vector{Float64}

    function Results()
        new(Array{Int64}([]), Array{Float64}([]), Array{Float64}([]), Array{Float64}([]))
    end
end

"""Append execution measures to Results."""
function logresult!(r::Results, currentiter, elapsedtime, fvaluegap, metricLP)

    push!(r.iterations, currentiter)
    push!(r.times, elapsedtime)
    push!(r.fvaluegaps, fvaluegap)
    push!(r.metricLPs, metricLP)
    return
end

max_stars_repo_name: ericlincc/Efficient-GLP

max_issues_repo_name: ericlincc/Efficient-GLP

max_forks_repo_name: ericlincc/Efficient-GLP


text:
using Flux
export LayerDef, create_model_arch, output_dims, calculate_final_outdims, create_model

mutable struct Model{T}
    paths::T
end

Model(paths...) = Model(paths)

Flux.@functor Model

(m::Model)(x::AbstractArray) = map(f -> f(x), m.paths)

function build_Model(path, input_dims, nclasses)
    model_defs = JSON.parse(String(read(path)))
    out = []
    for key in keys(model_defs)
        m = eval(Meta.parse(model_defs[key][1]))
        latent_dims = Flux.outputsize(m, tuple(vcat(input_dims, 1)...))
        push!(out, Chain(m, flatten, Dense(prod(latent_dims), nclasses)))
    end
    return Model(tuple(out...))
end


text:
# indicator of the Cartesian product of real binary sets

export IndBinary

"""
**Indicator of the product of binary sets**

    IndBinary(low, up)

Returns the indicator function of the set
```math
S = \\{ x : x_i = low_i\\ \\text{or}\\ x_i = up_i \\},
```
Parameters `low` and `up` can be either scalars or arrays of the same dimension as the space.
"""
struct IndBinary{T, S} <: ProximableFunction
    low::T
    high::S
end

is_set(f::IndBinary) = true

IndBinary() = IndBinary(0.0, 1.0)

IndBinary_low(f::IndBinary{<: Number, S}, i) where S = f.low
IndBinary_low(f::IndBinary{T, S}, i) where {T, S} = f.low[i]
IndBinary_high(f::IndBinary{T, <: Number}, i) where T = f.high
IndBinary_high(f::IndBinary{T, S}, i) where {T, S} = f.high[i]

function (f::IndBinary)(x::AbstractArray{C}) where {R <: Real, C <: RealOrComplex{R}}
    for k in eachindex(x)
        if x[k] != IndBinary_low(f, k) && x[k] != IndBinary_high(f, k)
            return R(Inf)
        end
    end
    return R(0)
end

function prox!(y::AbstractArray{T}, f::IndBinary, x::AbstractArray{T}, gamma::R=R(1)) where {R <: Real, T <: Union{R, Complex{R}}}
    for k in eachindex(x)
        low = IndBinary_low(f, k)
        high = IndBinary_high(f, k)
        if abs(x[k] - low) < abs(x[k] - high)
            y[k] = low
        else
            y[k] = high
        end
    end
    return R(0)
end

fun_name(f::IndBinary) = "indicator of binary array"
fun_dom(f::IndBinary) = "AbstractArray{Real}"

function prox_naive(f::IndBinary, x::AbstractArray{T}, gamma::R=R(1)) where {R <: Real, T <: Union{R, Complex{R}}}
    distlow = abs.(x .- f.low)
    disthigh = abs.(x .- f.high)
    indlow = distlow .< disthigh
    indhigh = distlow .>= disthigh
    y = f.low.*indlow + f.high.*indhigh
    return y, R(0)
end


text:
const ALWB_URI = URI(scheme="http", host="www.bom.gov.au", path="/jsp/awra/thredds/fileServer/AWRACMS")

abstract type DataMode end

"""
    Values <: DataMode

Get the dataset as regular measured values.
"""
struct Values <: DataMode end

"""
    Deciles <: DataMode

Get the dataset in relative deciles.
"""
struct Deciles <: DataMode end

# Docs below
struct ALWB{M<:DataMode,D<:Union{Day,Month,Year}} <: RasterDataSource end

layers(::Type{<:ALWB}) = (
    :rain_day, :s0_pct, :ss_pct, :sd_pct, :sm_pct, :qtot, :etot,
    :e0, :ma_wet, :pen_pet, :fao_pet, :asce_pet, :msl_wet, :dd
)

# Days are in 1 year nc files
date_step(::Type{<:ALWB{<:Any,Day}}) = Year(1)
# Months and years are in single files
date_step(::Type{<:ALWB{<:Any,Month}}) = Year(100)
date_step(::Type{<:ALWB{<:Any,Year}}) = Year(100)

has_constant_dims(::Type{<:ALWB}) = false

@doc """
    ALWB{Union{Deciles,Values},Union{Day,Month,Year}} <: RasterDataSource

Data from the Australian Landscape Water Balance (ALWB) data source.

See: [www.bom.gov.au/water/landscape](http://www.bom.gov.au/water/landscape)

The dataset contains NetCDF files. They have a time dimension so that multiple
dates are stored in each file.

The available layers are: `$(layers(ALWB))`, available in daily, monthly and
annual resolutions, and as `Values` or relative `Deciles`.

`getraster` for `ALWB` must use a `date` keyword to specify the date to download.

See the [`getraster`](@ref) docs for implementation details.
""" ALWB

# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/rain_day_2017.nc
# Precipiation = "rain_day"

# SoilMoisture_Upper = "s0_pct"
# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/ss_pct_2017.nc
# SoilMoisture_Lower = "ss_pct"
# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/sd_pct_2017.nc
# SoilMoisture_Deep = "sd_pct"
# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/sm_pct_2017.nc
# SoilMoisture_RootZone = "sm_pct" # http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/qtot_2017.nc # Runoff = "qtot"
# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/etot_2017.nc
# Evapotrans_Actual = "etot"
# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/e0_2017.nc
# Evapotrans_Potential_Landscape = "e0"
# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/ma_wet_2017.nc
# Evapotrans_Potential_Areal = "ma_wet"
# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/pen_pet_2017.nc
# Evapotrans_Potential_SyntheticPan = "pen_pet"

# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/fao_pet_2017.nc
# Evapotrans_RefCrop_Short = "fao_pet"
# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/asce_pet_2017.nc
# Evapotrans_RefCrop_Tall = "asce_pet"
# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/etot_2017.nc

# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/msl_wet_2017.nc
# Evaporation_OpenWater = "msl_wet"

# http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/dd_2017.nc
# DeepDrainage = "dd"

"""
    getraster(source::Type{<:ALWB{Union{Deciles,Values},Union{Day,Month,Year}}}, [layer]; date)

Download [`ALWB`](@ref) weather data from
[www.bom.gov.au/water/landscape](http://www.bom.gov.au/water/landscape) as values or
deciles with timesteps of `Day`, `Month` or `Year`.

# Arguments
- `layer`: `Symbol` or `Tuple` of `Symbol` from `$(layers(ALWB))`. Without a
    `layer` argument, all layers will be downloaded, and a `NamedTuple` of paths returned.

# Keywords
- `date`: a `DateTime`, `AbstractVector` of `DateTime` or a `Tuple` of start and end dates.
    For multiple dates, a `Vector` of multiple filenames will be returned.
    ALWB is available with a daily, monthly, and yearly, timestep.

# Example
This will return the file containing annual averages, including your date:

```julia
julia> getraster(ALWB{Values,Year}, :ss_pct; date=Date(2001, 2))
"/your/RASTERDATASOURCES_PATH/ALWB/values/month/ss_pct.nc"
```

Returns the filepath/s of the downloaded or pre-existing files.
"""
function getraster(T::Type{<:ALWB}, layers::Union{Tuple,Symbol}; date)
     _getraster(T, layers, date)
end

function _getraster(T::Type{<:ALWB{M,P}}, layers, dates::Tuple) where {M,P}
    _getraster(T, layers, date_sequence(T, dates))
end
function _getraster(T::Type{<:ALWB}, layers, dates::AbstractArray)
    _getraster.(T, Ref(layers), dates)
end
function _getraster(T::Type{<:ALWB}, layers::Tuple, date::Dates.TimeType)
    _map_layers(T, layers, date)
end
function _getraster(T::Type{<:ALWB}, layer::Symbol, date::Dates.TimeType)
    _check_layer(T, layer)
    mkpath(rasterpath(T))
    url = rasterurl(T, layer; date=date)
    path = rasterpath(T, layer; date=date)
    _maybe_download(url, path)
    path
end

rastername(T::Type{<:ALWB{M,P}}, layer; date) where {M,P} =
    string(layer, _pathsegment(P, date), ".nc")

rasterpath(::Type{ALWB}) = joinpath(rasterpath(), "ALWB")
rasterpath(::Type{ALWB{M,P}}) where {M,P} =
    joinpath(joinpath(rasterpath(), "ALWB"), map(_pathsegment, (M, P))...)
rasterpath(T::Type{<:ALWB}, layer; date=nothing) =
    joinpath(rasterpath(T), rastername(T, layer; date))

rasterurl(T::Type{<:ALWB{M,P}}, layer; date) where {M,P} =
    joinpath(ALWB_URI, _pathsegments(T)..., rastername(T, layer; date))

# Utility methods

_pathsegments(::Type{ALWB{M,P}}) where {M,P} = _pathsegment(M), _pathsegment(P)
_pathsegment(::Type{Values}) = "values"
_pathsegment(::Type{Deciles}) = "deciles"
_pathsegment(::Type{Day}) = "day"
_pathsegment(::Type{Month}) = "month"
_pathsegment(::Type{Year}) = "year"
# Days are in whole-year files
_pathsegment(::Type{Day}, date) = "_" * string(year(date))
# Months and years are all in one file
_pathsegment(::Type{<:Union{Year,Month}}, date) = ""

max_forks_repo_name: cesaraustralia/RasterDataSources


text:
using Test

location(s::Symbol) = (s === :u ? (Face, Center, Center) :
                       s === :v ? (Center, Face, Center) :
                       s === :w ? (Center, Center, Face) :
                                  (Center, Center, Center))

print_min_max_mean(ψ, name="") =
    @info @sprintf("%s min: %.9e, max: %.9e, mean: %.9e", name, minimum(ψ), maximum(ψ), mean(ψ))

function extract_two_solutions(analytical_solution, filename; name=:u)
    grid = RectilinearGrid(filename)
    iters = iterations(filename)
    loc = location(name)

    ψ_raw = field_data(filename, name, iters[end])

    ψ_raw = ψ_raw[
                  1 : grid.Nx + 2 * grid.Hx,
                  1 : grid.Ny + 2 * grid.Hy,
                  1 : grid.Nz + 2 * grid.Hz
                 ]

    tx, ty, tz = size(ψ_raw)
    ψ_data = OffsetArray(ψ_raw, 0:tx-1, 0:ty-1, 0:tz-1)
    ψ_simulation = Field{loc[1], loc[2], loc[3]}(ψ_data, grid, FieldBoundaryConditions(grid, loc))

    x, y, z = nodes(ψ_simulation, reshape=true)

    ψ_simulation = interior(ψ_simulation)

    t = iteration_time(filename, iters[end])
    ψ_analytical = analytical_solution.(x, y, z, t)

    print_min_max_mean(ψ_simulation, "simulation")
    print_min_max_mean(ψ_analytical, "analytical")

    return ψ_simulation, ψ_analytical
end

function compute_error(u_simulation, u_analytical)
    # Convert any CuArrays to Arrays.
    u_simulation = Array(u_simulation)
    u_analytical = Array(u_analytical)

    absolute_error = @. abs(u_simulation - u_analytical)
    absolute_truth = abs.(u_analytical)
    L₁ = mean(absolute_error)
    L₂ = mean(absolute_error.^2)
    L∞ = maximum(absolute_error)

    return (L₁=L₁, L∞=L∞)
end

function compute_error(analytical_solution, filename::String; kwargs...)
    u_simulation, u_analytical = extract_two_solutions(analytical_solution, filename; kwargs...)
    return compute_error(u_simulation, u_analytical)
end

compute_errors(analytical_solution, filenames::String...; kwargs...) =
    [compute_error(analytical_solution, filename; kwargs...) for filename in filenames]

extract_sizes(filenames...) = [size(RectilinearGrid(filename)) for filename in filenames]

function test_rate_of_convergence(error, N; name="", Ntest=N[end], expected, atol)
    i = searchsortedfirst(N, Ntest)
    ROC = log10(error[i-1] / error[i]) / log10(N[i-1] / N[i])
    @info (name == "" ? "" : name * " ") * "rate of convergence = $ROC (expected ≈ $expected, atol = $atol)"
    @test isapprox(ROC, expected, atol=atol)
    return ROC
end


text:
using Random

function tick(rng::AbstractRNG)

end

function run_model()
    println("hi")
end


text:
using TimeWarp
using TimeWarp.Datasets
using Plots
gr()

data,labels = ucr_traindata("50words")

plot([ fastdtw(data[:,i], data[:,i+10], radius)[1] for radius=1:100, i=1:10 ], show=true)
xlabel!("radius")


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm, SparseArrays

sz=[3, 3]

# V,(VV,EV,FV,CV) = Lar.cuboidGrid([3,3,3], true)
# V,(VV,EV,FV,CV) = Lar.cuboidGrid([sz[1], sz[2], sz[3]], false)
# V,(VV,EV,FV,CV) = Lar.cuboidGrid(sz, true)
V,(VV,EV,FV) = Lar.cuboidGrid(sz, true)

# Plasm.view( Plasm.numbering(.6)((V,[VV, EV, FV])) )
# Plasm.view( Plasm.numbering(.6)((V,[VV, EV, FV])) )
Plasm.view( Plasm.numbering(.6)((V,[VV, EV])) )


text:
function Base.:+{N,ES}(lhs::Valid{N,ES}, rhs::Valid{N,ES})

  (isempty(lhs)    || isempty(rhs))    && (return Valid{N,ES}(∅))
  (isallreals(lhs) || isallreals(rhs)) && (return Valid{N,ES}(ℝp))

  roundsinf_lhs::Bool = roundsinf(lhs)
  roundsinf_rhs::Bool = roundsinf(rhs)

  (roundsinf_lhs && roundsinf_rhs) && (return Valid{N,ES}(ℝp))

  result = ((@lower rhs) + (@lower lhs)) → ((@upper rhs) + (@upper lhs))

  if (roundsinf_lhs || roundsinf_rhs)
    # let's make sure our result still rounds inf, this is a property which is
    # invariant under addition.  Losing this property suggests that the answer
    # should be recast as "allreals."  While we're at it, check to see if the
    # answer ends now "touch", which makes them "allreals".
    roundsinf(result) || (return Valid{N,ES}(ℝp))
  end

  return result
end

Base.:-{N,ES}(lhs::Valid{N,ES}, rhs::Valid{N,ES}) = lhs + (-rhs)
Base.:-{N,ES}(x::Valid{N,ES}) = Valid{N,ES}(-x.upper, -x.lower)

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
using LinearAlgebra
using StaticArrays
using Distributions
using PyPlot
using TimerOutputs
const Vec2f = SVector{2, Float64}
const Vec4f = SVector{4, Float64}
const VecFrame = SVector{18, Float64}

include("src/geometry/geometry.jl")
include("src/world.jl")
include("src/double_integrator.jl")
include("src/fmt_cvae.jl")
include("read_scenario.jl")
include("predict.jl")
include("PyCallModel.jl")
include("PyCallPredict.jl")

lane1_c = 7.231027419583393
lane2_c = 18.84918099161813
lane3_c = 29.990037979406054
lane4_c = 40.99582664476257
lane5_c = 53.099915981198585
lane_c = [7.231027419583393, 18.84918099161813, 29.990037979406054,
          40.99582664476257, 53.099915981198585]
lane12 = (lane1_c + lane2_c)/2
lane23 = (lane2_c + lane3_c)/2
lane34 = (lane3_c + lane4_c)/2
lane45 = (lane4_c + lane5_c)/2
lane5edge = lane5_c + (lane5_c - lane45)

#us101_1_cluster.7 no case 35
#us101_1_cluster.9 no case 28
#us101_3_cluster.9 no case 3

#us101_2_cluster.9 case 45 invalid sample = 500

#us101_2_cluster.9 case 7 not easy to reach sample = 100
#us101_2_cluster.9 case 28, 38 not easy to reach sample = 100
#us101_2_cluster.9 case 45 not reach sample = 100
#us101_2_cluster.10 case 30 not easy to reach sample = 100

#us101_3_cluster.12 case 11 not easy to reach sample = 100

#sample = 50:
#us101_1_cluster.7 no case 20, 25
for case_id in 1:100
    print("case id: ",case_id, "\n")
    sleep(3)
    case_name = "./us101_2/cluster_10/test_case_file"*string(case_id)*".csv"
    host_file_name = "./us101_2/cluster_10/test_host_file"*string(case_id)*".csv"
    global scenario, start_frame = read_scenario(case_name)
    global lane1_c, lane2_c, lane3_c, lane4_c, lane5_c, lane12, lane23, lane34, lane45, lane5edge

    global init_goal = get_init_goal(host_file_name)

    global s_init = init_goal[1]
    global s_goal = init_goal[4]
    global n_goal = 2

    frame_id = 0
    curr_t = 0
    global lane_id = 0
    global left_prob = 0
    global keep_prob = 0
    global right_prob = 0

    result_f = open("./us101_2_cluster_100/cluster10/test_"*string(case_id)*".txt", "a+")

    # global to = TimerOutput()
    reset_timer!()

    while frame_id < 40
        println("frame id:", frame_id)
        global to
        global start_frame
        global s_goal
        global isOK = true
        global new_start = false
        if frame_id == 0
            s_init = init_goal[1]
            s_goal = init_goal[4]
            if 0 < s_init[1] < lane12
                global lane_id = 1
            elseif lane12 < s_init[1] < lane23
                global lane_id = 2
            elseif lane23 < s_init[1] < lane34
                global lane_id = 3
            elseif lane34 < s_init[1] < lane45
                global lane_id = 4
            elseif lane45 < s_init[1] < lane5edge
                global lane_id = 5
            end
            # println(length(frame))
            global prev_frame = get_frame(scenario, frame_id-1, start_frame)
            global frame1 = get_frame(scenario, frame_id, start_frame)
            #global frame2, frame3, frame4, frame5, frame6, frame7, frame8, frame9, frame10 = predict_action(prev_frame, frame1, frame_id+start_frame)
            global frame2 = get_frame(scenario, frame_id+1, start_frame)
            global frame3 = get_frame(scenario, frame_id+2, start_frame)
            global frame4 = get_frame(scenario, frame_id+3, start_frame)
            global frame5 = get_frame(scenario, frame_id+4, start_frame)
            global frame6 = get_frame(scenario, frame_id+5, start_frame)
            global frame7 = get_frame(scenario, frame_id+6, start_frame)
            global frame8 = get_frame(scenario, frame_id+7, start_frame)
            global frame9 = get_frame(scenario, frame_id+8, start_frame)
            global frame10 = get_frame(scenario, frame_id+9, start_frame)

            center = Vec2f(s_init[1], s_init[2])
            # println(length(frame))
            global neighbors1, front_neighbors1 = find_neighbors(frame1, center, lane_id, 150.0)
            global neighbors2, front_neighbors2 = find_neighbors(frame2, center, lane_id, 150.0)
            global neighbors3, front_neighbors3 = find_neighbors(frame3, center, lane_id, 150.0)
            global neighbors4, front_neighbors4 = find_neighbors(frame4, center, lane_id, 150.0)
            global neighbors5, front_neighbors5 = find_neighbors(frame5, center, lane_id, 150.0)
            global neighbors6, front_neighbors6 = find_neighbors(frame6, center, lane_id, 150.0)
            global neighbors7, front_neighbors7 = find_neighbors(frame7, center, lane_id, 150.0)
            global neighbors8, front_neighbors8 = find_neighbors(frame8, center, lane_id, 150.0)
            global neighbors9, front_neighbors9 = find_neighbors(frame9, center, lane_id, 150.0)
            global neighbors10, front_neighbors10 = find_neighbors(frame10, center, lane_id, 150.0)
            # println(length(neighbors))
            obstacle_set = set_obstacles(neighbors1)
        else
            global traj, future_wy, reach = get_next_position(cum_time, waypoints, 0.3, s_goal)
            if !reach
                # println("trajvalid1")
                trajvalid1 = check_traj(cum_time, waypoints, 0.1, frame_id+1, s_goal, scenario, start_frame)
                # println("trajvalid2")
                trajvalid2 = check_traj(cum_time, waypoints, 0.2, frame_id+2, s_goal, scenario, start_frame)
                # println("trajvalid3")
                trajvalid3 = check_traj(cum_time, waypoints, 0.3, frame_id+3, s_goal, scenario, start_frame)
                # println("trajvalid4")
                # trajvalid4 = check_traj(cum_time, waypoints, 0.4, frame_id+4, s_goal, scenario, start_frame)
                # trajvalid5 = check_traj(cum_time, waypoints, 0.5, frame_id+5, s_goal, scenario, start_frame)
                # trajvalid6 = check_traj(cum_time, waypoints, 0.6, frame_id+6, s_goal, scenario, start_frame)
                # trajvalid7 = check_traj(cum_time, waypoints, 0.7, frame_id+7, s_goal, scenario, start_frame)
                # trajvalid8 = check_traj(cum_time, waypoints, 0.8, frame_id+8, s_goal, scenario, start_frame)
                #trajvalid9 = check_traj(cum_time, waypoints, 0.9, frame_id+9, s_goal, scenario, start_frame)
                #trajvalid10 = check_traj(cum_time, waypoints, 1.0, frame_id+10, s_goal, scenario, start_frame)
                if !trajvalid1
                    println("traj1 invalid")
                end
                if !trajvalid2
                    println("traj2 invalid")
                end
                if !trajvalid3
                    println("traj3 invalid")
                end
                # if !trajvalid4
                #     println("traj4 invalid")
                # end
                # if !trajvalid5
                #     println("traj5 invalid")
                # end
                # if !trajvalid6
                #     println("traj6 invalid")
                # end
                # if !trajvalid7
                #     println("traj7 invalid")
                # end
                # if !trajvalid8
                #     println("traj8 invalid")
                # end
                #if !trajvalid9
                #    println("traj9 invalid")
                #end
                #if !trajvalid10
                #    println("traj10 invalid")
                #end
                # if !(trajvalid1 && trajvalid2&& trajvalid3)
                #    println("collision!!!!!")
                #    exit()
                # end
            end

            if reach
                println("congratulations!!!successful！！！")
                for p in traj
                    write(result_f, string(p[1]))
                    write(result_f, ',')
                    write(result_f, string(p[2]))
                    write(result_f, ',')
                    write(result_f, string(p[3]))
                    write(result_f, ',')
                    write(result_f, string(p[4]))
                    write(result_f, '\n')
                end
                print_timer()
                break
            end
            global s_init = traj[end]
            print("init:", s_init)
            print("goal:", s_goal)
            println("goal distance:", sqrt((s_init[1]-s_goal[1])^2+(s_init[2]-s_goal[2])^2))
            if sqrt((s_init[1]-s_goal[1])^2+(s_init[2]-s_goal[2])^2+(s_init[3]-s_goal[3])^2+(s_init[4]-s_goal[4])^2) < 10
                println("congratulations!!!successful！！！")
                break
            end
            global prev_frame = get_frame(scenario, frame_id-1, start_frame)
            global frame1 = get_frame(scenario, frame_id, start_frame)
            # global frame2, frame3, frame4, frame5, frame6, frame7, frame8, frame9, frame10 = predict_action(prev_frame, frame1, frame_id+start_frame)
            global frame2 = get_frame(scenario, frame_id+1, start_frame)
            global frame3 = get_frame(scenario, frame_id+2, start_frame)
            global frame4 = get_frame(scenario, frame_id+3, start_frame)
            global frame5 = get_frame(scenario, frame_id+4, start_frame)
            global frame6 = get_frame(scenario, frame_id+5, start_frame)
            global frame7 = get_frame(scenario, frame_id+6, start_frame)
            global frame8 = get_frame(scenario, frame_id+7, start_frame)
            global frame9 = get_frame(scenario, frame_id+8, start_frame)
            global frame10 = get_frame(scenario, frame_id+9, start_frame)

            global center = Vec2f(s_init[1], s_init[2])
            if 0 < s_init[1] < lane12
                global lane_id = 1
            elseif lane12 < s_init[1] < lane23
                global lane_id = 2
            elseif lane23 < s_init[1] < lane34
                global lane_id = 3
            elseif lane34 < s_init[1] < lane45
                global lane_id = 4
            elseif lane45 < s_init[1] < lane5edge
                global lane_id = 5
            end
            global neighbors1, front_neighbors1 = find_neighbors(frame1, center, lane_id, 150.0)
            global neighbors2, front_neighbors2 = find_neighbors(frame2, center, lane_id, 150.0)
            global neighbors3, front_neighbors3 = find_neighbors(frame3, center, lane_id, 150.0)
            global neighbors4, front_neighbors4 = find_neighbors(frame4, center, lane_id, 150.0)
            global neighbors5, front_neighbors5 = find_neighbors(frame5, center, lane_id, 150.0)
            global neighbors6, front_neighbors6 = find_neighbors(frame6, center, lane_id, 150.0)
            global neighbors7, front_neighbors7 = find_neighbors(frame7, center, lane_id, 150.0)
            global neighbors8, front_neighbors8 = find_neighbors(frame8, center, lane_id, 150.0)
            global neighbors9, front_neighbors9 = find_neighbors(frame9, center, lane_id, 150.0)
            global neighbors10, front_neighbors10 = find_neighbors(frame10, center, lane_id, 150.0)
            global obstacle_set = set_obstacles(neighbors1)
            global x_min, x_max, v_min, v_max = get_boundary(neighbors1, center, lane_id, s_goal)
            global W = World(x_min, x_max, v_min, v_max, obstacle_set)
            if isValid(W, future_wy) && !new_start && length(future_wy) > 10 && false
                global cum_time = cal_time(future_wy)
                global waypoints = future_wy
                global M = zeros(4, length(waypoints))
                for i in 1:length(waypoints)
                    global M[:, i] = waypoints[i]
                end
                clf()
                show(W)
                scatter(s_init[1], s_init[2], c=:blue, s=50)
                scatter(s_goal[1], s_goal[2], c=:blue, s=50)
                plot(M[1, :], M[2, :], c=:blue, linewidth=1.5)
                savefig("./fig/test-cvae-real-data3-2-3-with-pre/real_data-"*string(frame_id)*".png")
                for p in traj
                    write(result_f, string(p[1]))
                    write(result_f, ',')
                    write(result_f, string(p[2]))
                    write(result_f, ',')
                    write(result_f, string(p[3]))
                    write(result_f, ',')
                    write(result_f, string(p[4]))
                    write(result_f, '\n')
                end
                frame_id += 3
                continue
            end
            for p in traj
                write(result_f, string(p[1]))
                write(result_f, ',')
                write(result_f, string(p[2]))
                write(result_f, ',')
                write(result_f, string(p[3]))
                write(result_f, ',')
                write(result_f, string(p[4]))
                write(result_f, '\n')
            end
        end

        lane1_v = 0
        lane2_v = 0
        lane3_v = 0
        lane4_v = 0
        lane5_v = 0
        for v in neighbors1
            next_lane_id = predict_vehicle(scenario, center, convert(Int64, v[1]), frame_id+5, start_frame)
            if next_lane_id == 1
                lane1_v += 1
            elseif next_lane_id == 2
                lane2_v += 1
            elseif next_lane_id == 3
                lane3_v += 1
            elseif next_lane_id == 4
            lane4_v += 1
            elseif next_lane_id == 5
                lane5_v += 1
            end
        end
        if lane_id == 1
            global keep_prob, right_prob = softmax(1/lane1_v, 1/lane2_v)
            global left_prob = 0
        elseif lane_id == 2
            global left_prob, keep_prob, right_prob = softmax(1/lane1_v, 1/lane2_v, 1/lane3_v)
        elseif lane_id == 3
            global left_prob, keep_prob, right_prob = softmax(1/lane2_v, 1/lane3_v, 1/lane4_v)
        elseif lane_id == 4
            global left_prob, keep_prob, right_prob = softmax(1/lane3_v, 1/lane4_v, 1/lane5_v)
        elseif lane_id == 5
            global left_prob, keep_prob = softmax(1/lane4_v, 1/lane5_v)
            global right_prob = 0
        end

        # obstacle_set = set_obstacles(neighbors)
        if left_prob == 0
            left_prob = 0.0
        end
        if keep_prob == 0
            keep_prob = 0.0
        end
        if right_prob == 0
            right_prob = 0.0
        end
        # println(left_prob)
        # println(keep_prob)
        # println(right_prob)
        x_min, x_max, v_min, v_max = get_boundary(neighbors1, center, lane_id, s_goal)
        W = World(x_min, x_max, v_min, v_max, obstacle_set)
        # show(W)
        # scatter(s_init[1], s_init[2])
        # scatter(s_goal[1], s_goal[2])
        # savefig("./us101_2_cluster_100/cluster10/test_0_world.png")
        # exit()
        prob = [left_prob, keep_prob, right_prob]

        # println("goal:", s_goal)

        Nsample = 1000
        # right_prob = 1/3
        # left_prob = 1/3
        # keep_prob = 1/3
        println("init:", s_init)
        println("goal:", s_goal)
        frames = [frame1, frame2, frame3, frame4, frame5, frame6, frame7, frame8, frame9, frame10]
        fmt = @time FMTree(s_init, s_goal, Nsample, W, right_prob, left_prob, keep_prob, lane_id, frames)
        println(length(fmt.Pset))
        with_savefig = false
        # println("fmt.s_init:", fmt.s_init)
        @timeit "solve "*string(frame_id) idx_solution = solve(fmt, frame_id, with_savefig)
        clf()
        show(fmt)

        t = Float64[]
        waypoints = Vec4f[]
        for idx in idx_solution
            s0 = fmt.Pset[fmt.parent[idx]]
            s1 = fmt.Pset[idx]
            tau = fmt.time[idx]
            println("s0:", s0)
            println("s1:", s1)
            # println(tau)
            points = show_trajectory(s0, s1, tau, 20, :blue, 1.5)
            if idx != 1
                waypoints = vcat(points, waypoints)
            end
            # print(waypoints)
            scatter(s0[1], s0[2], c=:black, s=10)
            # scatter(s1[1], s1[2], c=:red, s=10)
        end
        global waypoints
        global cum_time = cal_time(waypoints)
        # println(fmt.itr)
        # println(t)
        savefig("./us101_2_cluster_100/cluster10/test_case_"*string(case_id)*"_"*string(frame_id)*".png")
        frame_id += 3
        # break
        # print_timer()
    end
    open("./us101_2_cluster_100/cluster10/mytime"*string(case_id)*".txt","a+") do io
        print_timer(io, TimerOutputs.get_defaulttimer(), allocations=false)
    end
end
exit()


text:
module StaticValues

import Base: splitprec, truncbits, truncmask, twiceprecision, TwicePrecision, canonicalize2
import Base: eltype, promote_eltype, values, log10, isfinite, zero, iszero,
             abs, abs2, isless, max, min, div, rem, promote_rule, @pure
import Base: getindex, length, first, step, last, firstindex, lastindex, size, eltype

export SInt128, SInt16, SInt32, SInt64, SInt, SInt8,
       SUInt128, SUInt64, SUInt, SUInt32, SUInt16, SUInt8,
       SBool,
       SFloat16, SFloat32, SFloat64,
       AbstractInt64, AbstractInt32, AbstractInt16, AbstractInt8,
       AbstractUInt64, AbstractUInt32, AbstractUInt16, AbstractUInt8,
       AbstractFloat64, AbstractFloat32, AbstractFloat16,
       SSigned, SUnsigned, SInteger, SFloat, SReal, SNumber,
       SVal, TPVal,
       SOne, SZero,
       SChar, SSymbol,
       UnitSRange, UnitMRange,
#       OneToSRange,
       StaticStepRangeLen,
       StepSRangeLen,
       StepSRange, StepMRange,
       LinSRange, LinMRange,
       OffsetSRange, OffsetMRange,
       CartesianSIndices, CartesianMIndices, LinearSIndices, LinearMIndices,
       NamedRange,
       AbstractIndices, SIndices, MIndices

struct Dynamic end

const BaseUnsigned = Union{UInt128,UInt16,UInt32,UInt64,UInt8}
const BaseSigned = Union{Int128,Int16,Int32,Int64,Int8}
const BaseInteger = Union{Int128,Int16,Int32,Int64,Int8,
                          UInt128,UInt16,UInt32,UInt64,UInt8,Bool}
const BaseFloat = Union{Float64,Float32,Float16}
const BaseReal = Union{Int128,Int16,Int32,Int64,Int8,
                       UInt128,UInt16,UInt32,UInt64,UInt8,
                       Bool,Float64,Float32,Float16,Rational,Irrational}
const BaseNumber = Union{Int128,Int16,Int32,Int64,Int8,
                         UInt128,UInt16,UInt32,UInt64,UInt8,
                         Bool,Float64,Float32,Float16,
                         Rational,Irrational,
                         Complex}
const BaseAny = Union{<:BaseNumber,Char,TwicePrecision}

include("deffunctions.jl")
include("signed.jl")
include("unsigned.jl")
include("int.jl")
include("float.jl")
include("rational.jl")
include("real.jl")
include("complex.jl")
include("char.jl")

include("types.jl")
include("intfuncs.jl")
include("floatfuncs.jl")

include("twiceprecision.jl")

include("functions.jl")

#const IEEESFloat = Union{Float16,Float32,Float64,SFloat16,SFloat32,SFloat64}

const F_or_FF = Union{<:AbstractFloat, Tuple{<:AbstractFloat,<:AbstractFloat}}

f64(x::BaseFloat) = Float64(x)
f64(x::Tuple{<:BaseFloat,<:BaseFloat}) = Float64(x[1]) + Float64(x[2])
f64(x::SFloat) = SFloat64(x)
f64(x::Tuple{<:SFloat,<:SFloat}) = SFloat64(x[1]) + SFloat64(x[2])

tp64(x::BaseNumber) = TwicePrecision{Float64}(x)
tp64(x::Tuple{<:BaseNumber,<:BaseNumber}) = TwicePrecision{Float64}(x[1], x[2])
tp64(x::SReal) = TPVal(ofeltype(Float64, x))
tp64(x::Tuple{H,L}) where {H<:SReal,L<:SReal} = TPVal(Float64, x[1]::H, x[2]::L)
ofeltype(::Type{T}, x::BaseNumber) where T = convert(eltype(T), x)

include("staticrange.jl")
include("staticstartrange.jl")
include("unitrange.jl")
include("staticordinalrange.jl")
include("steprange.jl")
include("linrange.jl")
include("steprangelen.jl")
include("srange.jl")
include("colon.jl")
include("rangefuncs.jl")
include("offsetrange.jl")
include("staticrangewrapper.jl")
include("namedrange.jl")
include("indices.jl")

"""
    isstatic(x) -> Bool

Determine whether a type is static.

# Examples
```jldoctest
julia> isstatic(SOne)
true

julia> isstatic(1)
false
```
"""
function isstatic end

# TODO examples, and ensure all ranges have this
"""
    tostatic(x)

Returns the statically defined version of a type.
"""
function tostatic end

#include("math.jl")

end


text:
using MPSKit,TensorKit,Test

#we pass this function to dmrg, it gets executed once per iteration
function finalize(iter,state,ham,pars)
    final_bonddim = 12;

    upperbound = max_Ds(state);
    shouldincrease = reduce((a,i) -> a && dim(virtualspace(state,i))>=final_bonddim || upperbound[i+1] == dim(virtualspace(state,i)),1:length(state),init=true);

    if (shouldincrease)
        (state,pars) = changebonds(state, ham, OptimalExpand(),pars);
    end

    return (state,pars,true)
end

#defining the hamiltonian
th = nonsym_ising_ham(lambda = 4.0);
szt = TensorMap([1 0;0 -1],ℂ^2,ℂ^2)

#finite mps (dmrg)
ts = FiniteMPS(fill(TensorMap(rand,ComplexF64,ℂ^1*ℂ^2,ℂ^1),10));
(ts,pars,_)=find_groundstate(ts,th,Dmrg(finalize=finalize));

szval_finite= sum(expectation_value(ts,szt))/length(ts)
@test szval_finite ≈ 0 atol=1e-12

#twosite dmrg
ts = FiniteMPS(fill(TensorMap(rand,ComplexF64,ℂ^1*ℂ^2,ℂ^1),10));
(ts,pars,_)=find_groundstate(ts,th,Dmrg2(trscheme = truncdim(15)));

szval_finite= sum(expectation_value(ts,szt))/length(ts)
@test szval_finite ≈ 0 atol=1e-12

#uniform mps
ts=InfiniteMPS([ℂ^2],[ℂ^50]);
(ts,pars,_)=find_groundstate(ts,th,Vumps(maxiter=400));

szval_infinite=@tensor ts.AC[1][1,2,3]*szt[4,2]*conj(ts.AC[1][1,4,3])
@test szval_infinite ≈ 0 atol=1e-12

# gradient optimisation using the Grassmann manifold structure
ts = InfiniteMPS([ℂ^2], [ℂ^5]);
(ts, pars, _) = find_groundstate(ts, th, GradientGrassmann(maxiter=400));

szval_infinite = @tensor ts.AC[1][1,2,3]*szt[4,2]*conj(ts.AC[1][1,4,3])
@test szval_infinite ≈ 0 atol=1e-12


text:
abstract type AbstractID end
Base.UInt(id::AbstractID) = id.value
idstr(id::AbstractID) = string(UInt(id); base = 16, pad = Sys.WORD_SIZE >> 2)
idprefix(::AbstractID) = ""

Base.print(io::IO, id::AbstractID) = print(io, idprefix(id), idstr(id))

function Base.show(io::IO, id::AbstractID)
    if get(io, :typeinfo, Any) <: Union{typeof(id),Missing,Nothing}
        show(io, UInt(id))
    else
        invoke(show, Tuple{IO,Any}, io, id)
    end
end

struct TaskID <: AbstractID
    value::UInt
end

TaskID(task::Task) = TaskID(UInt(pointer_from_objref(task)))
idprefix(::TaskID) = "task:"

struct ObjectID <: AbstractID
    value::UInt

    function ObjectID(x)
        if ismutable(x)
            new(UInt(pointer_from_objref(x)))
        else
            new(objectid(x))
        end
    end
end

Historic.taskid(task::Task = current_task()) = TaskID(task)
Historic.objid(x) = ObjectID(x)

const LAST_UNIQUE_ID = UInt64[]
unique_id_index(tid = Threads.threadid()) = tid * 8

function init_unique_id()
    resize!(LAST_UNIQUE_ID, Threads.nthreads() * 8)
    fill!(LAST_UNIQUE_ID, 0)
    for tid in 1:Threads.nthreads()
        LAST_UNIQUE_ID[unique_id_index(tid)] = tid
    end
end

Historic.uniqueid() = LAST_UNIQUE_ID[unique_id_index()] += Threads.nthreads()


text:
"""
Three-stage Runge-Kutta method of Zhao and Wei with parameter C=4
[Math. Meth. Appl. Sci. 2014, 37 1042–1071, figure 1b]
"""
@generated function timestep!(op::CollisionOperator{M,N,ℳ,𝒩,RT,CT}, u₀::Matrix{RT}, u₁::Matrix{RT}, Δt::RT) where {M,N,ℳ,𝒩,RT,CT}
    local u::Vector{Matrix{RT}} = [zeros(RT,M,N), zeros(RT,M,N), zeros(RT,M,N)]
    local f::Vector{Matrix{RT}} = [zeros(RT,M,N), zeros(RT,M,N), zeros(RT,M,N)]

    local a::Vector{RT} = [0.0, 0.5, 1.0]
    local b::Vector{RT} = [1/6, 2/3, 1/6]
    local c::Vector{RT} = [0.0, 0.5, 1.0]

    quote
        @assert size(u₀) == size(u₁) == (M,N)

        $u[1] .= u₀
        collision_operator!(op, $u[1], $f[1])

        update_field!($f[1], u₀, $u[2], Δt * $a[2])
        collision_operator!(op, $u[2], $f[2])

        update_field!($f[2], u₀, $u[3], Δt * $a[3])
        collision_operator!(op, $u[3], $f[3])

        update_field!($f[1], u₀, u₁, Δt * $b[1])
        update_field!($f[2], u₁, u₁, Δt * $b[2])
        update_field!($f[3], u₁, u₁, Δt * $b[3])

        u₁
    end
end

function update_field!(f::Matrix{RT}, u₀::Matrix{RT}, u₁::Matrix{RT}, fac::RT) where {RT}
    @assert size(f) == size(u₀) == size(u₁)
    @inbounds for j in 1:size(u₁,2)
        for i in 1:size(u₁,1)
            u₁[i,j] = u₀[i,j] + fac * f[i,j]
        end
    end
end


text:
module CleanUpPullRequestPreviews

import GitHub
import HTTP

export remove_old_previews

struct AlwaysAssertionError <: Exception
    msg::String
end

function always_assert(cond::Bool, message::String)
    cond || throw(AlwaysAssertionError(message))
    return nothing
end

struct PRandPath
    pr::Int
    path::String
end

struct Config
    email_address::String
    git_command::String
    my_regex::Regex
    num_samples::Int
    push_to_origin::Bool
    repo_main::String
    repo_previews::String
    repo_previews_branch::String
    sample_delay_seconds::Int
    username::String
end

function Config(;
                email_address::AbstractString = "41898282+github-actions[bot]@users.noreply.github.com",
                git_command::AbstractString = "git",
                my_regex::Regex,
                num_samples::Integer = 3,
                push_to_origin::Bool = true,
                repo_main::AbstractString,
                repo_previews::AbstractString,
                repo_previews_branch::AbstractString,
                sample_delay_seconds::Integer = 60,
                username::AbstractString = "github-actions[bot]",
                )
    result = Config(
        email_address,
        git_command,
        my_regex,
        num_samples,
        push_to_origin,
        repo_main,
        repo_previews,
        repo_previews_branch,
        sample_delay_seconds,
        username,
    )
    return result
end

"""
    remove_old_previews

Remove old pull request previews.
"""
function remove_old_previews(;
                             api::GitHub.GitHubAPI = GitHub.GitHubWebAPI(HTTP.URI("https://api.github.com")),
                             auth::GitHub.Authorization = GitHub.authenticate(ENV["GITHUB_TOKEN"]),
                             kwargs...)
    config = Config(; kwargs...)
    remove_old_previews(api, auth, config)::Nothing
    return nothing
end

"""
    remove_old_previews

Remove old pull request previews.
"""
function remove_old_previews(api::GitHub.GitHubAPI,
                             auth::GitHub.Authorization,
                             config::Config)
    original_directory::String = pwd()::String
    clone_directory::String = _git_clone(api, auth, config)::String
    prs_and_paths::Vector{PRandPath} = _get_prs_and_paths(api, auth, config, clone_directory)::Vector{PRandPath}
    pr_is_open::Dict{Int, Bool} = _pr_is_open(api, auth, config, prs_and_paths)::Dict{Int, Bool}
    _delete_preview_directories(api, auth, config, clone_directory, prs_and_paths, pr_is_open)::Nothing
    _git_add_commit_push(api, auth, config, clone_directory)::Nothing
    cd(original_directory)
    rm(clone_directory; force = true, recursive = true)
    return nothing
end

function _git_clone(api::GitHub.GitHubAPI,
                    auth::GitHub.Authorization,
                    config::Config)::String
    parent_dir = mktempdir()
    atexit(() -> rm(parent_dir; force = true, recursive = true))
    original_directory = pwd()
    cd(parent_dir)
    run(`$(config.git_command) clone $(config.repo_previews) CLONEDIRECTORY`)
    clone_directory::String = joinpath(parent_dir, "CLONEDIRECTORY")::String
    cd(clone_directory)
    run(Cmd(String[config.git_command, "config", "user.name", strip(config.username)]))
    run(Cmd(String[config.git_command, "config", "user.email", strip(config.email_address)]))
    run(`$(config.git_command) checkout $(config.repo_previews_branch)`)
    cd(original_directory)
    return clone_directory
end

function _get_prs_and_paths(api::GitHub.GitHubAPI,
                            auth::GitHub.Authorization,
                            config::Config,
                            root_directory::AbstractString)::Vector{PRandPath}
    result::Vector{PRandPath} = Vector{PRandPath}(undef, 0)::Vector{PRandPath}
    original_directory = pwd()
    cd(root_directory)
    for (root, dirs, files) in walkdir(".")
        for dir in dirs
            path = joinpath(root, dir)
            m = match(config.my_regex, path)
            if m !== nothing
                pr = parse(Int, m[1])
                push!(result, PRandPath(pr, path))
            end
        end
    end
    cd(original_directory)
    return result
end

function _pr_is_open(api::GitHub.GitHubAPI,
                     auth::GitHub.Authorization,
                     config::Config,
                     prs_and_paths)::Dict{Int, Bool}
    list_of_all_pr_numbers::Vector{Int} = Int[pr_and_path.pr for pr_and_path in prs_and_paths]::Vector{Int}
    unique!(list_of_all_pr_numbers)
    sort!(list_of_all_pr_numbers)

    intermediate = Vector{Dict{Int, Bool}}(undef, config.num_samples)
    for i = 1:config.num_samples
        intermediate[i] = Dict{Int, Bool}()
    end

    for i = 1:config.num_samples
        if i != 1
            @info("Waiting for $(config.sample_delay_seconds) second(s)...")
            sleep(config.sample_delay_seconds)
        end
        for pr_number in list_of_all_pr_numbers
            @info("Getting state of PR # $(pr_number)")
            github_pr = GitHub.pull_request(api,
                                            config.repo_main,
                                            pr_number;
                                            auth = auth)
            github_pr_is_open = github_pr.state != "closed"
            intermediate[i][pr_number] = github_pr_is_open
        end
    end

    result::Dict{Int, Bool} = Dict{Int, Bool}()::Dict{Int, Bool}

    for pr_number in list_of_all_pr_numbers
        pr_intermediates = Vector{Bool}(undef, config.num_samples)
        for i = 1:config.num_samples
            pr_intermediates[i] = intermediate[i][pr_number]
        end
        result[pr_number] = any(pr_intermediates)
        @info("PR # $(pr_number) is open: $(result[pr_number])")
    end

    return result
end

function _git_add_commit_push(api::GitHub.GitHubAPI,
                              auth::GitHub.Authorization,
                              config::Config,
                              clone_directory::AbstractString)::Nothing
    original_directory = pwd()
    cd(clone_directory)
    run(`$(config.git_command) add -A`)

    p = run(pipeline(`$(config.git_command) commit -m "Automated commit created by CleanUpPullRequestPreviews.jl"`; stdout=stdout, stderr=stderr); wait = false)
    wait(p)
    @info("Commit succeeded = $(success(p))")

    if config.push_to_origin
        run(`$(config.git_command) push origin --all`)
    end

    cd(original_directory)
    return nothing
end

function _delete_preview_directories(api::GitHub.GitHubAPI,
                                     auth::GitHub.Authorization,
                                     config::Config,
                                     clone_directory::AbstractString,
                                     prs_and_paths::Vector{PRandPath},
                                     pr_is_open::Dict{Int, Bool})::Nothing
    original_directory = pwd()
    for pr_and_path in prs_and_paths
        pr = pr_and_path.pr
        path = pr_and_path.path
        this_pr_is_open = pr_is_open[pr]
        if !this_pr_is_open
            cd(clone_directory)
            @info("Removing \"$(path)\"")
            rm(path; force = true, recursive = true)
        end
    end
    cd(original_directory)
    return nothing
end

end  # end module CleanUpPullRequestPreviews


text:
function fitmm!(
    gcm::GLMCopulaVCModel,
    solver=Ipopt.IpoptSolver(print_level=0)
    )
    optm = MathProgBase.NonlinearModel(solver)
    lb = fill(-Inf, gcm.p)
    ub = fill( Inf, gcm.p)
    MathProgBase.loadproblem!(optm, gcm.p, 0, lb, ub, Float64[], Float64[], :Max, gcm)
    MathProgBase.setwarmstart!(optm, gcm.β)
    MathProgBase.optimize!(optm)
    optstat = MathProgBase.status(optm)
    optstat == :Optimal || @warn("Optimization unsuccesful; got $optstat")
    copy_par!(gcm, MathProgBase.getsolution(optm))
    loglikelihood!(gcm)
    # gcm
    nothing
end

function MathProgBase.initialize(
    gcm::GLMCopulaVCModel,
    requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in [:Grad, :Hess])
            error("Unsupported feature $feat")
        end
    end
end

MathProgBase.features_available(gcm::GLMCopulaVCModel) = [:Grad, :Hess]

function MathProgBase.eval_f(
    gcm::GLMCopulaVCModel,
    par::Vector)
    copy_par!(gcm, par)
    # update nuisance parameter
    if typeof(gcm.data[1].d) <: NegativeBinomial
        new_d = update_r!(gcm)
        @show new_d
    end
    # maximize σ2 and τ at current β using MM
    update_θ!(gcm)
    @show gcm.θ
    # evaluate loglikelihood
    loglikelihood!(gcm, false, false)
end

function MathProgBase.eval_grad_f(
    gcm::GLMCopulaVCModel,
    grad::Vector,
    par::Vector)
    # println("reached eval_grad_f")
    copy_par!(gcm, par)
    # maximize σ2 and τ at current β using MM
    @show gcm.β
    # update nuisance parameter
    if typeof(gcm.data[1].d) <: NegativeBinomial
        new_d = update_r!(gcm)
        @show new_d
    end
    update_θ!(gcm)
    @show gcm.θ
    # evaluate gradient
    logl = loglikelihood!(gcm, true, false)
    copyto!(grad, gcm.∇β)
    nothing
end

function copy_par!(
    gcm::GLMCopulaVCModel,
    par::Vector)
    copyto!(gcm.β, par)
    par
end

function MathProgBase.hesslag_structure(gcm::GLMCopulaVCModel)
    Iidx = Vector{Int}(undef, (gcm.p * (gcm.p + 1)) >> 1)
    Jidx = similar(Iidx)
    ct = 1
    for j in 1:gcm.p
        for i in j:gcm.p
            Iidx[ct] = i
            Jidx[ct] = j
            ct += 1
        end
    end
    Iidx, Jidx
end

function MathProgBase.eval_hesslag(
    gcm::GLMCopulaVCModel,
    H::Vector{T},
    par::Vector{T},
    σ::T,
    μ::Vector{T}) where {T}
    # println("reached eval_hesslag")
    copy_par!(gcm, par)
    if typeof(gcm.data[1].d) <: NegativeBinomial
        new_d = update_r!(gcm)
        @show new_d
    end
    # maximize σ2 and τ at current β using MM
    update_θ!(gcm)
    # evaluate Hessian
    loglikelihood!(gcm, true, true)
    # copy Hessian elements into H
    ct = 1
    for j in 1:gcm.p
        for i in j:gcm.p
            H[ct] = gcm.Hβ[i, j]
            ct += 1
        end
    end
    H .*= σ
end


text:
import JuliaDB: select, Not
import PrettyTables: pretty_table, tf_html_simple, TextFormat, tf_mysql
const _TF_COMPACT_ = TextFormat(up_right_corner = '─',
                   up_left_corner      = '─',
                   bottom_left_corner  = '─',
                   bottom_right_corner = '─',
                   up_intersection     = '─',
                   left_intersection   = '─',
                   right_intersection  = '─',
                   middle_intersection = '─',
                   bottom_intersection = '─',
                   column              = ' ',
                   row                 = '─',
                   hlines              = [:begin,:header,:end],
                   vlines              = :all);

function pretty_table(quran::AbstractQuran; kwargs...)
    if quran isa Chapter
        if quran.numbers isa Int64
            tbl = select(quran.data, Not(:chapter))
        else
            tbl = quran.data
        end
    elseif quran isa Verse
        if quran.chapters isa Int64 && quran.numbers isa Int64
            tbl = select(quran.data, Not(:chapter, :verse))
        elseif quran.chapters isa Int64 && quran.numbers isa UnitRange{Int64}
            tbl = select(quran.data, Not(:chapter))
        elseif quran.chapters isa UnitRange{Int64} && quran.numbers isa Int64
            tbl = quran.data
        else
            tbl = quran.data
        end
    elseif quran isa Word
        if quran.chapters isa Int64 && quran.verses isa Int64 && quran.numbers isa Int64
            tbl = select(quran.data, Not(:chapter, :verse, :word))
        elseif quran.chapters isa Int64 && quran.verses isa Int64 && quran.numbers isa UnitRange{Int64}
            tbl = select(quran.data, Not(:chapter, :verse))
        elseif quran.chapters isa Int64 && quran.verses isa UnitRange{Int64} && quran.numbers isa Int64
            tbl = select(quran.data, Not(:chapter))
        elseif quran.chapters isa Int64 && quran.verses isa UnitRange{Int64} && quran.numbers isa UnitRange{Int64}
            tbl = select(quran.data, Not(:chapter))
        elseif quran.chapters isa UnitRange{Int64} && quran.verses isa Int64 && quran.numbers isa Int64
            tbl = quran.data
        elseif quran.chapters isa UnitRange{Int64} && quran.verses isa Int64 && quran.numbers isa UnitRange{Int64}
            tbl = quran.data
        elseif quran.chapters isa UnitRange{Int64} && quran.verses isa UnitRange{Int64} && quran.numbers isa Int64
            tbl = quran.data
        else
            tbl = quran.data
        end
    elseif quran isa Part
        if quran.chapters isa Int64 && quran.verses isa Int64 && quran.words isa Int64 && quran.numbers isa Int64
            tbl = select(quran.data, Not(:chapter, :verse, :word, :part))
        elseif quran.chapters isa Int64 && quran.verses isa Int64 && quran.words isa Int64 && quran.numbers isa UnitRange{Int64}
            tbl = select(quran.data, Not(:chapter, :verse, :word))
        elseif quran.chapters isa Int64 && quran.verses isa Int64 && quran.words isa UnitRange{Int64} && quran.numbers isa UnitRange{Int64}
            tbl = select(quran.data, Not(:chapter, :verse))
        elseif quran.chapters isa Int64 && quran.verses isa Int64 && quran.words isa UnitRange{Int64} && quran.numbers isa Int64
            tbl = select(quran.data, Not(:chapter, :verse))
        elseif quran.chapters isa Int64 && quran.verses isa UnitRange{Int64} && quran.words isa UnitRange{Int64} && quran.numbers isa UnitRange{Int64}
            tbl = select(quran.data, Not(:chapter, :verse))
        elseif quran.chapters isa Int64 && quran.verses isa UnitRange{Int64} && quran.words isa UnitRange{Int64} && quran.numbers isa Int64
            tbl = select(quran.data, Not(:chapter, :verse))
        elseif quran.chapters isa Int64 && quran.verses isa UnitRange{Int64} && quran.words isa Int64 && quran.numbers isa UnitRange{Int64}
            tbl = select(quran.data, Not(:chapter))
        elseif quran.chapters isa UnitRange{Int64} && quran.verses isa UnitRange{Int64} && quran.words isa UnitRange{Int64} && quran.numbers isa UnitRange{Int64}
            tbl = quran.data
        elseif quran.chapters isa UnitRange{Int64} && quran.verses isa UnitRange{Int64} && quran.words isa UnitRange{Int64} && quran.numbers isa Int64
            tbl = quran.data
        elseif quran.chapters isa UnitRange{Int64} && quran.verses isa UnitRange{Int64} && quran.words isa Int64 && quran.numbers isa UnitRange{Int64}
            tbl = quran.data
        elseif quran.chapters isa UnitRange{Int64} && quran.verses isa Int64 && quran.words isa UnitRange{Int64} && quran.numbers isa UnitRange{Int64}
            tbl = quran.data
        elseif quran.chapters isa UnitRange{Int64} && quran.verses isa UnitRange{Int64} && quran.words isa Int64 && quran.numbers isa Int64
            tbl = quran.data
        elseif quran.chapters isa Int64 && quran.verses isa Int64 && quran.words isa UnitRange{Int64} && quran.numbers isa UnitRange{Int64}
            tbl = select(quran.data, Not(:chapter, :verse))
        elseif quran.chapters isa Int64 && quran.verses isa UnitRange{Int64} && quran.words isa UnitRange{Int64} && quran.numbers isa Int64
            tbl = select(quran.data, Not(:chapter))
        elseif quran.chapters isa UnitRange{Int64} && quran.verses isa Int64 && quran.words isa Int64 && quran.numbers isa UnitRange{Int64}
            tbl = quran.data
        else
            tbl = quran.data
        end
    else
        tbl = quran.data
    end
    pretty_table(tbl; kwargs...)
end

function description(feat::Union{Lemma,Root,Special})
    println(typeof(feat), ":")
    println(" └ data: ", feat.data)
end

function description(feat::AbstractPartOfSpeech)
    println(typeof(feat), ":"); j = 1
    fields = fieldnames(typeof(feat))
    for f in fields
        if j < length(fields)
            if length(fields) > 1
                println(" ├ ", f, ": ", getproperty(feat, f))
            else
                println(" └ ", f, ": ", getproperty(feat, f))
            end
        else
            println(" └ ", f, ": ", getproperty(feat, f))
        end
        j += 1
    end
end

function description(feat::Union{Prefix,Stem,Suffix})
    println(typeof(feat))
    println(join(repeat(["─"], length(string(typeof(feat))))))
    if length(fieldnames(typeof(feat))) < 3
        description(feat.pos)
    else
        description(feat.pos)
        for f in feat.feats
            description(f)
        end
    end
end

"""
    @desc(expr)

Extract the detailed description of a `AbstractFeature`.

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps);
julia> tnzldata = table(tnzl);
julia> feat = parse(Features, select(crpsdata.data, :features)[53])
julia> @desc feat
Stem
────
Negative:
 ├ data: NEG
 ├ desc: Negative particle
 └ ar_label: حرف نفي
Lemma:
 └ data: laA
Special:
 └ data: <in~
```
"""
macro desc(expr)
    esc(quote
        try
            return description(eval($expr))
        catch
            return missing
        end
    end)
end

"""
    @data(expr)

Extract the data property object.

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps);
julia> tnzldata = table(tnzl);
julia> feat = parse(Features, select(crpsdata.data, :features)[53])
julia> @data feat
:NEG
```
"""
macro data(expr)
    return esc(Meta.parse(string(expr) * ".data"))
end

function Base.show(io::IO, t::AbstractPartOfSpeech)
    print(io, t.data)
end

function Base.show(io::IO, t::SimpleEncoder)
    println(io, typeof(t), ":")
    println(io, " └ encoder: ", t.encode)
end

function Base.show(io::IO, m::MetaData)
    println(io, m.title, " v", m.version)
    println(io, "Copyright (C) ", m.year, " ", m.author)
    println(io, m.license)
    println(io, m.website, "\n")
    println(io, m.description)
end

function Base.show(io::IO, ::MIME"text/plain", c::CorpusRaw)
    print(io, "(CorpusRaw) ")
    show(io, "text/plain", c.data)
end

function Base.show(io::IO, ::MIME"text/plain", c::TanzilRaw)
    print(io, "(TanzilRaw) ")
    show(io, "text/plain", c.data)
end

function Base.show(io::IO, m::CorpusData)
    println(io, m.meta.title)
    println(io, "(C) ", m.meta.year, " ", m.meta.author, "\n")
    println(io, m.data)
end

function Base.show(io::IO, m::TanzilData)
    println(io, m.meta.title)
    println(io, "(C) ", m.meta.year, " ", m.meta.author, "\n")
    println(io, m.data)
end

function Base.show(io::IO, quran::Chapter)
    chapterlab = ChapterLabel()
    if quran.numbers isa Int64
        println(io, "Chapter ", quran.numbers, ": ", chapterlab.arabic[quran.numbers], " (", chapterlab.english[quran.numbers], ")\n")
        println(io, select(quran.data, Not(:chapter)))
    elseif quran.numbers isa Array{Int64}
        println(io, "Chapters: "); j = 1
        for i in quran.numbers
            if j < length(quran.numbers)
                if length(quran.numbers) > 1
                    print(io, " ├")
                    println(io, lpad(i, 3), " (", chapterlab.arabic[i], "-", chapterlab.english[i], ") ")
                else
                    print(io, " └")
                    print(io, i, " (", chapterlab.arabic[i], "-", chapterlab.english[i], ") ")
                end
            else
                print(io, " └")
                print(io, lpad(i, 3), " (", chapterlab.arabic[i], "-", chapterlab.english[i], ")")
            end
            j += 1
        end
        println(io, "\n")
        println(io, quran.data)
    else
        println(io, "Chapter ", quran.numbers.start, "-", quran.numbers.stop, ": ",
            chapterlab.arabic[quran.numbers.start], "-", chapterlab.arabic[quran.numbers.stop],
            " (", chapterlab.english[quran.numbers.start], "-", chapterlab.english[quran.numbers.stop], ")\n")
        println(io, quran.data)
    end
end

function Base.show(io::IO, quran::Union{Verse,Word,Part})
    if quran isa Union{Word,Part}
        verse = quran.verses
    else
        verse = quran.numbers
    end
    chapterlab = ChapterLabel()
    if quran.chapters isa Int64 && verse isa Int64
        println(io, "Chapter ", quran.chapters, " ", chapterlab.arabic[quran.chapters], " (", chapterlab.english[quran.chapters], ")")
        println(io, "Verse ", verse, "\n")
        println(io, select(quran.data, Not(:chapter, :verse)))
    elseif quran.chapters isa Int64 && verse isa UnitRange{Int64}
        println(io, "Chapter ", quran.chapters, " ", chapterlab.arabic[quran.chapters], " (", chapterlab.english[quran.chapters], ")")
        println(io, "Verses ", verse.start, "-", verse.stop, "\n")
        println(io, select(quran.data, Not(:chapter)))
    elseif quran.chapters isa Int64 && verse isa Array{Int64,1}
        println(io, "Chapter ", quran.chapters, " ", chapterlab.arabic[quran.chapters], " (", chapterlab.english[quran.chapters], ")")
        if length(verse) > 1
            print(io, "Verses ")
        else
            print(io, "Verse ")
        end
        j = 1
        for i in verse
            if j < length(verse)
                print(io, i, ", ")
            else
                println(io, i)
            end
            j += 1
        end
        println(io, "\n", quran.data)
    elseif quran.chapters isa UnitRange{Int64} && verse isa Int64
        println(io, "Chapters ", quran.chapters.start, "-", quran.chapters.stop, ": ",
            chapterlab.arabic[quran.chapters.start], "-", chapterlab.arabic[quran.chapters.stop],
            " (", chapterlab.english[quran.chapters.start], "-", chapterlab.english[quran.chapters.stop], ")")
        println(io, "Verse ", verse, "\n")
        println(io, quran.data)
    elseif quran.chapters isa UnitRange{Int64} && verse isa Array{Int64,1}
        println(io, "Chapters ", quran.chapters.start, "-", quran.chapters.stop, ": ",
            chapterlab.arabic[quran.chapters.start], "-", chapterlab.arabic[quran.chapters.stop],
            " (", chapterlab.english[quran.chapters.start], "-", chapterlab.english[quran.chapters.stop], ")")
        if length(verse) > 1
            print(io, "Verses ")
        else
            print(io, "Verse ")
        end
        j = 1
        for i in verse
            if j < length(verse)
                print(io, i, ", ")
            else
                println(io, i)
            end
            j += 1
        end
        println(io, "\n", quran.data)
    elseif quran.chapters isa UnitRange{Int64} && verse isa UnitRange{Int64}
        println(io, "Chapters ", quran.chapters.start, "-", quran.chapters.stop, ": ",
            chapterlab.arabic[quran.chapters.start], "-", chapterlab.arabic[quran.chapters.stop],
            " (", chapterlab.english[quran.chapters.start], "-", chapterlab.english[quran.chapters.stop], ")")
        println(io, "Verses ", verse.start, "-", verse.stop, "\n")
        println(io, quran.data)
    elseif quran.chapters isa Array{Int64,1} && verse isa UnitRange{Int64}
        if length(quran.chapters) > 1
            println(io, "Chapters: ")
        else
            println(io, "Chapter: ")
        end

        j = 1
        for i in quran.chapters
            if j < length(quran.chapters)
                if length(quran.chapters) > 1
                    print(io, " ├")
                    println(io, lpad(i, 3), " (", chapterlab.arabic[i], "-", chapterlab.english[i], ") ")
                else
                    print(io, " └")
                    print(io, i, " (", chapterlab.arabic[i], "-", chapterlab.english[i], ") ")
                end
            else
                print(io, " └")
                print(io, lpad(i, 3), " (", chapterlab.arabic[i], "-", chapterlab.english[i], ")")
            end
            j += 1
        end
        println(io, "\nVerses ", verse.start, "-", verse.stop, "\n")
        println(io, quran.data)
    elseif quran.chapters isa Array{Int64,1} && verse isa Int64
        if length(quran.chapters) > 1
            println(io, "Chapters: ")
        else
            println(io, "Chapter: ")
        end

        j = 1
        for i in quran.chapters
            if j < length(quran.chapters)
                if length(quran.chapters) > 1
                    print(io, " ├")
                    println(io, lpad(i, 3), " (", chapterlab.arabic[i], "-", chapterlab.english[i], ") ")
                else
                    print(io, " └")
                    print(io, i, " (", chapterlab.arabic[i], "-", chapterlab.english[i], ") ")
                end
            else
                print(io, " └")
                print(io, lpad(i, 3), " (", chapterlab.arabic[i], "-", chapterlab.english[i], ")")
            end
            j += 1
        end
        println(io, "\nVerse ", verse, "\n")
        println(io, quran.data)
    else
        if length(quran.chapters) > 1
            println(io, "Chapters: ")
        else
            println(io, "Chapter: ")
        end

        j = 1
        for i in quran.chapters
            if j < length(quran.chapters)
                if length(quran.chapters) > 1
                    print(io, " ├")
                    println(io, lpad(i, 3), " (", chapterlab.arabic[i], "-", chapterlab.english[i], ") ")
                else
                    print(io, " └")
                    print(io, i, " (", chapterlab.arabic[i], "-", chapterlab.english[i], ") ")
                end
            else
                print(io, " └")
                print(io, lpad(i, 3), " (", chapterlab.arabic[i], "-", chapterlab.english[i], ")")
            end
            j += 1
        end
        if length(verse) > 1
            print(io, "\nVerses ")
        else
            print(io, "\nVerse ")
        end
        j = 1
        for i in verse
            if j < length(verse)
                print(io, i, ", ")
            else
                println(io, i)
            end
            j += 1
        end
        println(io, "\n", quran.data)
    end
end


text:
module particle

import Base.==
import ..operators

using ..bases, ..states, ..operators, ..operators_sparse, ..operators_lazy

export PositionBasis, MomentumBasis,
        gaussianstate,
        spacing, samplepoints,
        positionoperator, momentumoperator, laplace_x, laplace_p, FFTOperator

"""
Basis for a particle in real space.

For simplicity periodic boundaries are assumed which means that
the rightmost point defined by xmax is not included in the basis
but is defined to be the same as xmin.

Arguments
---------

xmin
    Left boundary.
xmax
    Right boundary.
N
    Number of discrete points in the basis.
"""
type PositionBasis <: Basis
    shape::Vector{Int}
    xmin::Float64
    xmax::Float64
    N::Int
    PositionBasis(xmin::Real, xmax::Real, N::Int) = new([N], xmin, xmax, N)
end

"""
Basis for a particle in momentum space.

For simplicity periodic boundaries are assumed which means that
pmax is not included in the basis but is defined to be the same as pmin.

Arguments
---------

pmin
    Left boundary.
pmax
    Right boundary.
N
    Number of discrete points in the basis.
"""
type MomentumBasis <: Basis
    shape::Vector{Int}
    pmin::Float64
    pmax::Float64
    N::Int
    MomentumBasis(pmin::Real, pmax::Real, N::Int) = new([N], pmin, pmax, N)
end

"""
Create a PositionBasis from a MomentumBasis.

Since periodic boundary conditions are assumed the values for :math:`x_{min}`
and :math:`x_{max}` are more or less arbitrary and are chosen to be
:math:`-\\pi/dp` and :math:`\\pi/dp` with :math:`dp=(p_{max}-p_{min})/N`.
"""
PositionBasis(b::MomentumBasis) = (dp = (b.pmax - b.pmin)/b.N; PositionBasis(-pi/dp, pi/dp, b.N))

"""
Create a MomentumBasis from a PositionBasis.

Since periodic boundary conditions are assumed the values for :math:`p_{min}`
and :math:`p_{max}` are more or less arbitrary and are chosen to be
:math:`-\\pi/dx` and :math:`\\pi/dx` with :math:`dx=(x_{max}-x_{min})/N`.
"""
MomentumBasis(b::PositionBasis) = (dx = (b.xmax - b.xmin)/b.N; MomentumBasis(-pi/dx, pi/dx, b.N))

==(b1::PositionBasis, b2::PositionBasis) = b1.xmin==b2.xmin && b1.xmax==b2.xmax && b1.N==b2.N
==(b1::MomentumBasis, b2::MomentumBasis) = b1.pmin==b2.pmin && b1.pmax==b2.pmax && b1.N==b2.N

"""
Create a Gaussian state around x0 and p0 with width sigma in real space.

This implementation is based on the definition

.. math::

    \\Psi(x) = \\frac{\\sqrt{\\Delta x}}{\\pi^{1/4}\\sqrt{\\sigma}}
                e^{i p_0 (x-x_0) - \\frac{(x-x_0)^2}{2 \\sigma^2}}

which has the properties

* :math:`|\\Psi(p)|^2 = 1`
* :math:`\\langle p \\rangle = p_0`
* :math:`\\langle x \\rangle = x_0`
* :math:`\\mathrm{Var}(x) = \\frac{\\sigma^2}{2}`
* :math:`\\mathrm{Var}(p) = \\frac{1}{2 \\sigma^2}`

and is connected to the MomentumBasis based implementation via a
Fourier-transformation

.. math::

    \\Psi(p) = \\frac{1}{\\sqrt{2\\pi}}
                \\int_{-\\infty}^{\\infty} e^{-ipx}\\Psi(x) \\mathrm{d}x

Due to the numerically necessary discretization an additional scaling
factor :math:`\\sqrt{\\Delta x}` is used so that
:math:`\\Psi_i = \\sqrt{\\Delta x} \\Psi(x_i)` and the resulting state
is normalized to

.. math::

    \\sum_i |\\Psi_i|^2 = 1

Arguments
---------
b
    PositionBasis.
x0
    Center of the Gaussian in position space.
p0
    Center of the Gaussian in momentum space.
sigma
    Width of the Gaussian.
"""
function gaussianstate(b::PositionBasis, x0::Real, p0::Real, sigma::Real)
    psi = Ket(b)
    dx = spacing(b)
    alpha = 1./(pi^(1/4)*sqrt(sigma))*sqrt(dx)
    x = b.xmin
    for i=1:b.N
        psi.data[i] = alpha*exp(1im*p0*(x-x0) - (x-x0)^2/(2*sigma^2))
        x += dx
    end
    return psi
end

"""
Create a Gaussian state around x0 and p0 with width sigma in momentum space.

This implementation is based on the definition

.. math::

    \\Psi(p) = \\frac{\\sqrt{\\sigma} \\sqrt{\\Delta x}}{\\pi^{1/4}}
                e^{-i x_0 p - \\frac{1}{2}(p-p_0)^2 \\sigma^2}

which has the properties

* :math:`|\\Psi(p)|^2 = 1`
* :math:`\\langle p \\rangle = p_0`
* :math:`\\langle x \\rangle = x_0`
* :math:`\\mathrm{Var}(x) = \\frac{\\sigma^2}{2}`
* :math:`\\mathrm{Var}(p) = \\frac{1}{2 \\sigma^2}`

and is connected to the MomentumBasis based implementation via a
Fourier-transformation

.. math::

    \\Psi(x) = \\frac{1}{\\sqrt{2\\pi}}
                    \\int_{-\\infty}^{\\infty} e^{ipx}\\Psi(p) \\mathrm{d}p

Due to the numerically necessary discretization an additional scaling
factor :math:`\\sqrt{\\Delta x}` is used so that
:math:`\\Psi_i = \\sqrt{\\Delta p} \\Psi(p_i)` and the resulting state
is normalized to

.. math::

    \\sum_i |\\Psi_i|^2 = 1

Arguments
---------
b
    MomentumBasis.
x0
    Center of the Gaussian in position space.
p0
    Center of the Gaussian in momentum space.
sigma
    Width of the Gaussian.
"""
function gaussianstate(b::MomentumBasis, x0::Real, p0::Real, sigma::Real)
    psi = Ket(b)
    dp = spacing(b)
    alpha = sqrt(sigma)/pi^(1/4)*sqrt(dp)
    p = b.pmin
    for i=1:b.N
        psi.data[i] = alpha*exp(-1im*x0*p - (p-p0)^2/2*sigma^2)
        p += dp
    end
    return psi
end

"""
Distance between two adjacent points of the real space basis.
"""
spacing(b::MomentumBasis) = (b.pmax - b.pmin)/b.N

"""
Momentum difference between two adjacent points of the momentum basis.
"""
spacing(b::PositionBasis) = (b.xmax - b.xmin)/b.N

"""
x values of the real space basis.
"""
samplepoints(b::PositionBasis) = (dx = spacing(b); Float64[b.xmin + i*dx for i=0:b.N-1])

"""
p values of the momentum basis.
"""
samplepoints(b::MomentumBasis) = (dp = spacing(b); Float64[b.pmin + i*dp for i=0:b.N-1])

"""
Position operator in real space.
"""
positionoperator(b::PositionBasis) = SparseOperator(b, spdiagm(complex(samplepoints(b)), 0, length(b), length(b)))

"""
Position operator in momentum space.
"""
function positionoperator(b::MomentumBasis)
    p_op = DenseOperator(b)
    u = 1im/(12*spacing(b))
    for i=1:b.N-2
        p_op.data[i,i+2] = -u
        p_op.data[i,i+1] = 8*u
        p_op.data[i+1,i] = -8*u
        p_op.data[i+2,i] = u
    end
    return SparseOperator(p_op)
end

"""
Momentum operator in momentum space.
"""
momentumoperator(b::MomentumBasis) = SparseOperator(b, spdiagm(complex(samplepoints(b)), 0, length(b), length(b)))

"""
Momentum operator in real space.
"""
function momentumoperator(b::PositionBasis)
    p_op = DenseOperator(b)
    u = -1im/(12*spacing(b))
    for i=1:b.N-2
        p_op.data[i,i+2] = -u
        p_op.data[i,i+1] = 8*u
        p_op.data[i+1,i] = -8*u
        p_op.data[i+2,i] = u
    end
    return SparseOperator(p_op)
end

"""
Second x-derivative in real space.
"""
function laplace_x(b::PositionBasis)
    x_op = DenseOperator(b)
    u = 1/spacing(b)^2
    for i=1:b.N-1
        x_op.data[i+1,i] = u
        x_op.data[i,i] = -2*u
        x_op.data[i,i+1] = u
    end
    x_op.data[b.N,b.N] = -2*u
    return SparseOperator(x_op)
end

"""
Second x-derivative in momentum space.
"""
laplace_x(b::MomentumBasis) = SparseOperator(b, spdiagm(complex(samplepoints(b)).^2, 0, length(b), length(b)))

"""
Second p-derivative in momentum space.
"""
function laplace_p(b::MomentumBasis)
    p_op = DenseOperator(b)
    u = 1/spacing(b)^2
    for i=1:b.N-1
        p_op.data[i+1,i] = u
        p_op.data[i,i] = -2*u
        p_op.data[i,i+1] = u
    end
    p_op.data[b.N,b.N] = -2*u
    return SparseOperator(p_op)
end

"""
Second p-derivative in real space.
"""
laplace_p(b::PositionBasis) = SparseOperator(b, spdiagm(complex(samplepoints(b)).^2, 0, length(b), length(b)))

type FFTOperator <: LazyOperator
    basis_l::Basis
    basis_r::Basis
    fft_l!
    fft_r!
    mul_before::Vector{Complex128}
    mul_after::Vector{Complex128}
end

"""
DenseOperator performing a fast Fourier transformation of a state.
"""
function FFTOperator(basis_l::MomentumBasis, basis_r::PositionBasis)
    Lx = (basis_r.xmax - basis_r.xmin)
    dp = spacing(basis_l)
    dx = spacing(basis_r)
    if basis_l.N != basis_r.N || abs(2*pi/dp - Lx)/Lx > 1e-12
        throw(IncompatibleBases())
    end
    x = zeros(Complex128, length(basis_r))
    mul_before = exp(-1im*basis_l.pmin*(samplepoints(basis_r)-basis_r.xmin))
    mul_after = exp(-1im*basis_r.xmin*samplepoints(basis_l))/sqrt(basis_r.N)
    FFTOperator(basis_l, basis_r, plan_bfft!(x), plan_fft!(x), mul_before, mul_after)
end

function FFTOperator(basis_l::PositionBasis, basis_r::MomentumBasis)
    Lx = (basis_l.xmax - basis_l.xmin)
    dp = spacing(basis_r)
    dx = spacing(basis_l)
    if basis_l.N != basis_r.N || abs(2*pi/dp - Lx)/Lx > 1e-12
        throw(IncompatibleBases())
    end
    x = zeros(Complex128, length(basis_r))
    mul_before = exp(1im*basis_l.xmin*(samplepoints(basis_r)-basis_r.pmin))
    mul_after = exp(1im*basis_r.pmin*samplepoints(basis_l))/sqrt(basis_r.N)
    FFTOperator(basis_l, basis_r, plan_fft!(x), plan_bfft!(x), mul_before, mul_after)
end

operators.dagger(op::FFTOperator) = FFTOperator(op.basis_r, op.basis_l)

function operators.gemv!{T<:Complex}(alpha::T, M::FFTOperator, b::Ket, beta::T, result::Ket)
    N::Int = M.basis_r.N
    if beta==Complex(0.)
        @inbounds for i=1:N
            result.data[i] = M.mul_before[i] * b.data[i]
        end
        M.fft_r! * result.data
        @inbounds for i=1:N
            result.data[i] *= M.mul_after[i] * alpha
        end
    else
        psi_ = Ket(M.basis_l, deepcopy(b.data))
        @inbounds for i=1:N
            psi_.data[i] *= M.mul_before[i]
        end
        M.fft_r! * psi_.data
        @inbounds for i=1:N
            result.data[i] = beta*result.data[i] + alpha * psi_.data[i] * M.mul_after[i]
        end
    end
    nothing
end

function operators.gemv!{T<:Complex}(alpha::T, b::Bra, M::FFTOperator, beta::T, result::Bra)
    N::Int = M.basis_l.N
    if beta==Complex(0.)
        @inbounds for i=1:N
            result.data[i] = conj(M.mul_after[i]) * conj(b.data[i])
        end
        M.fft_l! * result.data
        @inbounds for i=1:N
            result.data[i] = conj(result.data[i]) * M.mul_before[i] * alpha
        end
    else
        psi_ = Bra(M.basis_r, conj(b.data))
        @inbounds for i=1:N
            psi_.data[i] *= conj(M.mul_after[i])
        end
        M.fft_l! * psi_.data
        @inbounds for i=1:N
            result.data[i] = beta*result.data[i] + alpha * conj(psi_.data[i]) * M.mul_before[i]
        end
    end
    nothing
end

end # module