text:
"""
ODEForwardSensitivityFunction{iip,F,A,Tt,OJ,J,JP,S,PJ,TW,TWt,UF,PF,JC,PJC,Alg,fc,JM,pJM,MM,CV} <: DiffEqBase.AbstractODEFunction{iip}

ODEForwardSensitivityFunction is an internal to the ODEForwardSensitivityProblem which extends the AbstractODEFunction
to be used in an ODEProblem, but defines the tools requires for calculating the extra differential equations associated
with the derivative terms.

ODEForwardSensitivityFunction is not intended to be part of the public API.
"""
struct ODEForwardSensitivityFunction{iip,F,A,Tt,OJ,J,JP,S,PJ,TW,TWt,UF,PF,JC,PJC,Alg,fc,JM,pJM,MM,CV} <: DiffEqBase.AbstractODEFunction{iip}
  f::F
  analytic::A
  tgrad::Tt
  original_jac::OJ
  jac::J
  jac_prototype::JP
  sparsity::S
  paramjac::PJ
  Wfact::TW
  Wfact_t::TWt
  uf::UF
  pf::PF
  J::JM
  pJ::pJM
  jac_config::JC
  paramjac_config::PJC
  alg::Alg
  numparams::Int
  numindvar::Int
  f_cache::fc
  mass_matrix::MM
  isautojacvec::Bool
  isautojacmat::Bool
  colorvec::CV
end
has_original_jac(S) = isdefined(S, :original_jac) && S.jac !== nothing

struct NILSSForwardSensitivityFunction{iip,sensefunType,senseType,MM} <: DiffEqBase.AbstractODEFunction{iip}
  S::sensefunType
  sensealg::senseType
  nus::Int
  mass_matrix::MM
end

function ODEForwardSensitivityFunction(f,analytic,tgrad,original_jac,jac,jac_prototype,sparsity,paramjac,Wfact,Wfact_t,uf,pf,u0,
                                    jac_config,paramjac_config,alg,p,f_cache,mm,
                                    isautojacvec,isautojacmat,colorvec,nus)
  numparams = length(p)
  numindvar = length(u0)
  J = isautojacvec ? nothing : Matrix{eltype(u0)}(undef,numindvar,numindvar)
  pJ = Matrix{eltype(u0)}(undef,numindvar,numparams) # number of funcs size

  sensefun = ODEForwardSensitivityFunction{isinplace(f),typeof(f),typeof(analytic),
                                           typeof(tgrad),typeof(original_jac),
                                           typeof(jac),typeof(jac_prototype),typeof(sparsity),
                                           typeof(paramjac),
                                           typeof(Wfact),typeof(Wfact_t),typeof(uf),
                                           typeof(pf),typeof(jac_config),
                                           typeof(paramjac_config),typeof(alg),
                                           typeof(f_cache),
                                           typeof(J),typeof(pJ),typeof(mm),typeof(f.colorvec)}(
                           f,analytic,tgrad,original_jac,jac,jac_prototype,
                           sparsity,paramjac,Wfact,Wfact_t,uf,pf,J,pJ,
                           jac_config,paramjac_config,alg,
                           numparams,numindvar,f_cache,mm,isautojacvec,isautojacmat,colorvec,
                          )
  if nus!==nothing
    sensefun = NILSSForwardSensitivityFunction{isinplace(f), typeof(sensefun),
                             typeof(alg),typeof(mm)}(sensefun,alg,nus,mm)
  end

  return sensefun
end

function (S::ODEForwardSensitivityFunction)(du,u,p,t)
  y = @view u[1:S.numindvar] # These are the independent variables
  dy = @view du[1:S.numindvar]
  S.f(dy,y,p,t) # Make the first part be the ODE

  # Now do sensitivities
  # Compute the Jacobian

  if !S.isautojacvec && !S.isautojacmat
    if has_original_jac(S)
      S.original_jac(S.J,y,p,t) # Calculate the Jacobian into J
    else
      S.uf.t = t
      jacobian!(S.J, S.uf, y, S.f_cache, S.alg, S.jac_config)
    end
  end

  if DiffEqBase.has_paramjac(S.f)
    S.paramjac(S.pJ,y,p,t) # Calculate the parameter Jacobian into pJ
  else
    S.pf.t = t
    copyto!(S.pf.u,y)
    jacobian!(S.pJ, S.pf, p, S.f_cache, S.alg, S.paramjac_config)
  end

  # Compute the parameter derivatives
  if !S.isautojacvec && !S.isautojacmat
    dp = @view du[reshape(S.numindvar+1:(length(p)+1)*S.numindvar,S.numindvar,length(p))]
    Sj = @view u[reshape(S.numindvar+1:(length(p)+1)*S.numindvar,S.numindvar,length(p))]
    mul!(dp,S.J,Sj)
    DiffEqBase.@.. dp += S.pJ
  elseif S.isautojacmat
    S.uf.t = t
    Sj = @view u[reshape(S.numindvar+1:end,S.numindvar,S.numparams)]
    dp = @view du[reshape(S.numindvar+1:end,S.numindvar,S.numparams)]
    jacobianmat!(dp, S.uf, y, Sj, S.alg, S.jac_config)
    DiffEqBase.@.. dp += S.pJ
  else
    S.uf.t = t
    for i in eachindex(p)
      Sj = @view u[i*S.numindvar+1:(i+1)*S.numindvar]
      dp = @view du[i*S.numindvar+1:(i+1)*S.numindvar]
      jacobianvec!(dp, S.uf, y, Sj, S.alg, S.jac_config)
      dp .+= @view S.pJ[:,i]
    end
  end
  return nothing
end
@deprecate ODELocalSensitivityProblem(args...;kwargs...) ODEForwardSensitivityProblem(args...;kwargs...)

struct ODEForwardSensitivityProblem{iip,A}
  sensealg::A
end

function ODEForwardSensitivityProblem(f::F,args...;kwargs...) where F
  ODEForwardSensitivityProblem(ODEFunction(f),args...;kwargs...)
end

function ODEForwardSensitivityProblem(prob::ODEProblem,alg;kwargs...)
  ODEForwardSensitivityProblem(prob.f,prob.u0,prob.tspan,prob.p,alg;kwargs...)
end

const FORWARD_SENSITIVITY_PARAMETER_COMPATABILITY_MESSAGE =
"""
ODEForwardSensitivityProblem requires being able to solve
a differential equation defined by the parameter struct `p`. Thus while
DifferentialEquations.jl can support any parameter struct type, usage
with ODEForwardSensitivityProblem requires that `p` could be a valid
type for being the initial condition `u0` of an array. This means that
many simple types, such as `Tuple`s and `NamedTuple`s, will work as
parameters in normal contexts but will fail during ODEForwardSensitivityProblem
construction. To work around this issue for complicated cases like nested structs,
look into defining `p` using `AbstractArray` libraries such as RecursiveArrayTools.jl
or ComponentArrays.jl.
"""

struct ForwardSensitivityParameterCompatibilityError <: Exception end

function Base.showerror(io::IO, e::ForwardSensitivityParameterCompatibilityError)
  print(io, FORWARD_SENSITIVITY_PARAMETER_COMPATABILITY_MESSAGE)
end

@doc doc"""
function ODEForwardSensitivityProblem(f::Union{Function,DiffEqBase.AbstractODEFunction},
                                      u0,tspan,p=nothing,
                                      alg::AbstractForwardSensitivityAlgorithm = ForwardSensitivity();
                                      kwargs...)

Local forward sensitivity analysis gives a solution along with a timeseries of
the sensitivities. Thus if one wishes to have a derivative at every possible
time point, directly using the `ODEForwardSensitivityProblem` can be the most
efficient method.

!!! warning

      ODEForwardSensitivityProblem requires being able to solve
      a differential equation defined by the parameter struct `p`. Thus while
      DifferentialEquations.jl can support any parameter struct type, usage
      with ODEForwardSensitivityProblem requires that `p` could be a valid
      type for being the initial condition `u0` of an array. This means that
      many simple types, such as `Tuple`s and `NamedTuple`s, will work as
      parameters in normal contexts but will fail during ODEForwardSensitivityProblem
      construction. To work around this issue for complicated cases like nested structs,
      look into defining `p` using `AbstractArray` libraries such as RecursiveArrayTools.jl
      or ComponentArrays.jl.

### ODEForwardSensitivityProblem Syntax

`ODEForwardSensitivityProblem` is similar to an `ODEProblem`, but takes an
`AbstractForwardSensitivityAlgorithm` that describes how to append the forward sensitivity
equation calculation to the time evolution to simultaneously compute the derivative
of the solution with respect to parameters.

```julia
ODEForwardSensitivityProblem(f::SciMLBase.AbstractODEFunction,u0,
                             tspan,p=nothing,
                             sensealg::AbstractForwardSensitivityAlgorithm = ForwardSensitivity();
                             kwargs...)
```

Once constructed, this problem can be used in `solve` just like any other ODEProblem.
The solution can be deconstructed into the ODE solution and sensitivities parts using the
`extract_local_sensitivities` function, with the following dispatches:

```julia
extract_local_sensitivities(sol, asmatrix::Val=Val(false)) # Decompose the entire time series
extract_local_sensitivities(sol, i::Integer, asmatrix::Val=Val(false)) # Decompose sol[i]
extract_local_sensitivities(sol, t::Union{Number,AbstractVector}, asmatrix::Val=Val(false)) # Decompose sol(t)
```

For information on the mathematics behind these calculations, consult
[the sensitivity math page](@ref sensitivity_math)

### Example using an ODEForwardSensitivityProblem

To define a sensitivity problem, simply use the `ODEForwardSensitivityProblem` type
instead of an ODE type. For example, we generate an ODE with the sensitivity
equations attached for the Lotka-Volterra equations by:

```julia
function f(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + u[1]*u[2]
end

p = [1.5,1.0,3.0]
prob = ODEForwardSensitivityProblem(f,[1.0;1.0],(0.0,10.0),p)
```

This generates a problem which the ODE solvers can solve:

```julia
sol = solve(prob,DP8())
```

Note that the solution is the standard ODE system and the sensitivity system combined.
We can use the following helper functions to extract the sensitivity information:

```julia
x,dp = extract_local_sensitivities(sol)
x,dp = extract_local_sensitivities(sol,i)
x,dp = extract_local_sensitivities(sol,t)
```

In each case, `x` is the ODE values and `dp` is the matrix of sensitivities
The first gives the full timeseries of values and `dp[i]` contains the time series of the
sensitivities of all components of the ODE with respect to `i`th parameter.
The second returns the `i`th time step, while the third
interpolates to calculate the sensitivities at time `t`. For example, if we do:

```julia
x,dp = extract_local_sensitivities(sol)
da = dp[1]
```

then `da` is the timeseries for ``\frac{\partial u(t)}{\partial p}``. We can
plot this

```julia
plot(sol.t,da',lw=3)
```

transposing so that the rows (the timeseries) is plotted.

![Local Sensitivity Solution](https://user-images.githubusercontent.com/1814174/170916167-11d1b5c6-3c3c-439a-92af-d3899e24d2ad.png)

Here we see that there is a periodicity to the sensitivity which matches
the periodicity of the Lotka-Volterra solutions. However, as time goes on the
sensitivity increases. This matches the analysis of Wilkins in Sensitivity
Analysis for Oscillating Dynamical Systems.

We can also quickly see that these values are equivalent to those given by
automatic differentiation and numerical differentiation through the ODE solver:

```julia
using ForwardDiff, Calculus
function test_f(p)
  prob = ODEProblem(f,eltype(p).([1.0,1.0]),eltype(p).((0.0,10.0)),p)
  solve(prob,Vern9(),abstol=1e-14,reltol=1e-14,save_everystep=false)[end]
end

p = [1.5,1.0,3.0]
fd_res = ForwardDiff.jacobian(test_f,p)
calc_res = Calculus.finite_difference_jacobian(test_f,p)
```

Here we just checked the derivative at the end point.

### Internal representation of the Solution

For completeness, we detail the internal representation. When using
ForwardDiffSensitivity, the representation is with `Dual` numbers under the
standard interpretation. The values for the ODE's solution at time `i` are the
`ForwardDiff.value.(sol[i])` portions, and the derivative with respect to
parameter `j` is given by `ForwardDiff.partials.(sol[i])[j]`.

When using ForwardSensitivity, the solution to the ODE are the first `n`
components of the solution. This means we can grab the matrix of solution
values like:

```julia
x = sol[1:sol.prob.indvars,:]
```

Since each sensitivity is a vector of derivatives for each function, the sensitivities
are each of size `sol.prob.indvars`. We can pull out the parameter sensitivities from
the solution as follows:

```julia
da = sol[sol.prob.indvars+1:sol.prob.indvars*2,:]
db = sol[sol.prob.indvars*2+1:sol.prob.indvars*3,:]
dc = sol[sol.prob.indvars*3+1:sol.prob.indvars*4,:]
```

This means that `da[1,i]` is the derivative of the `x(t)` by the parameter `a`
at time `sol.t[i]`. Note that all of the functionality available to ODE solutions
is available in this case, including interpolations and plot recipes (the recipes
will plot the expanded system).
"""
function ODEForwardSensitivityProblem(f::F,u0,
                                    tspan,p=nothing,
                                    alg::ForwardSensitivity = ForwardSensitivity();
                                    nus=nothing, # determine if Nilss is used
                                    w0=nothing,
                                    v0=nothing,
                                    kwargs...) where F<:DiffEqBase.AbstractODEFunction
  isinplace = DiffEqBase.isinplace(f)
  # if there is an analytical Jacobian provided, we are not going to do automatic `jac*vec`
  isautojacmat = get_jacmat(alg)
  isautojacvec = get_jacvec(alg)
  p === nothing && error("You must have parameters to use parameter sensitivity calculations!")

  if !(typeof(p) <: Union{Nothing,SciMLBase.NullParameters,AbstractArray}) || (p isa AbstractArray && !Base.isconcretetype(eltype(p)))
    throw(ForwardSensitivityParameterCompatibilityError())
  end

  uf = DiffEqBase.UJacobianWrapper(f,tspan[1],p)
  pf = DiffEqBase.ParamJacobianWrapper(f,tspan[1],copy(u0))
  if isautojacmat
    if alg_autodiff(alg)
      jac_config_seed = ForwardDiff.Dual{typeof(uf)}.(u0,[ntuple(x -> zero(eltype(u0)), length(p)) for i in eachindex(u0)])
      jac_config_buffer = similar(jac_config_seed)
      jac_config = jac_config_seed, jac_config_buffer
    else
      error("Jacobian matrix products only work with automatic differentiation.")
    end
  elseif isautojacvec
    if alg_autodiff(alg)
      # if we are using automatic `jac*vec`, then we need to use a `jac_config`
      # that is a tuple in the form of `(seed, buffer)`
      jac_config_seed = ForwardDiff.Dual{typeof(jacobianvec!)}.(u0,u0)
      jac_config_buffer = similar(jac_config_seed)
      jac_config = jac_config_seed, jac_config_buffer
    else
      jac_config = (similar(u0),similar(u0))
    end
  elseif DiffEqBase.has_jac(f)
    jac_config = nothing
  else
    jac_config = build_jac_config(alg,uf,u0)
  end

  if DiffEqBase.has_paramjac(f)
    paramjac_config = nothing
  else
    paramjac_config = build_param_jac_config(alg,pf,u0,p)
  end

  # TODO: make it better
  if f.mass_matrix isa UniformScaling
    mm = f.mass_matrix
  else
    nn = size(f.mass_matrix, 1)
    mm = zeros(eltype(f.mass_matrix), (length(p)+1)*nn, (length(p)+1)*nn)
    mm[1:nn, 1:nn] = f.mass_matrix
    for i = 1:length(p)
      mm[i*nn+1:(i+1)nn, i*nn+1:(i+1)nn] = f.mass_matrix
    end
  end

  # TODO: Use user tgrad. iW can be safely ignored here.
  sense = ODEForwardSensitivityFunction(f,f.analytic,nothing,f.jac,nothing,
                                        nothing,nothing,f.paramjac,
                                        nothing,nothing,
                                        uf,pf,u0,jac_config,
                                        paramjac_config,alg,
                                        p,similar(u0),mm,
                                        isautojacvec,isautojacmat,f.colorvec,nus)
  if nus===nothing
    sense_u0 = [u0;zeros(eltype(u0),sense.numindvar*sense.numparams)]
  else
    if w0===nothing && v0===nothing
      sense_u0 = [u0;zeros(eltype(u0),(nus+1)*sense.S.numindvar*sense.S.numparams)]
    else
      sense_u0 = [u0;w0;v0]
    end
  end
  ODEProblem(sense,sense_u0,tspan,p,
             ODEForwardSensitivityProblem{DiffEqBase.isinplace(f),
                                          typeof(alg)}(alg);
             kwargs...)
end

function seed_duals(x::AbstractArray{V},f,
                    ::ForwardDiff.Chunk{N} = ForwardDiff.Chunk(x,typemax(Int64)),
                    ) where {V,T,N}
  seeds = ForwardDiff.construct_seeds(ForwardDiff.Partials{N,V})
  duals = ForwardDiff.Dual{typeof(ForwardDiff.Tag(f,eltype(vec(x))))}.(vec(x),seeds)
end

has_continuous_callback(cb::DiscreteCallback) = false
has_continuous_callback(cb::ContinuousCallback) = true
has_continuous_callback(cb::CallbackSet) = !isempty(cb.continuous_callbacks)

function ODEForwardSensitivityProblem(f::DiffEqBase.AbstractODEFunction,u0,
                                      tspan,p,alg::ForwardDiffSensitivity;
                                      du0=zeros(eltype(u0),length(u0),length(p)), # perturbations of initial condition
                                      dp=I(length(p)), # perturbations of parameters
                                      kwargs...)
  num_sen_par = size(du0,2)
  if num_sen_par != size(dp,2)
    error("Same number of perturbations of initial conditions and parameters required")
  end
  if size(du0,1) != length(u0)
    error("Perturbations for all initial conditions required")
  end
  if size(dp,1) != length(p)
    error("Perturbations for all parameters required")
  end

  pdual = ForwardDiff.Dual{typeof(ForwardDiff.Tag(f,eltype(vec(p))))}.(p, [ntuple(j -> dp[i,j], num_sen_par) for i in eachindex(p)])
  u0dual = ForwardDiff.Dual{typeof(ForwardDiff.Tag(f,eltype(vec(u0))))}.(u0, [ntuple(j -> du0[i,j], num_sen_par) for i in eachindex(u0)])

  if (convert_tspan(alg) === nothing &&
    haskey(kwargs,:callback) && has_continuous_callback(kwargs.callback)
    ) || (convert_tspan(alg) !== nothing && convert_tspan(alg))
    tspandual = convert.(eltype(pdual),tspan)
  else
    tspandual = tspan
  end

  prob_dual = ODEProblem(f,u0dual,tspan,pdual,
                         ODEForwardSensitivityProblem{DiffEqBase.isinplace(f),
                                                      typeof(alg)}(alg);
                         kwargs...)
end

"""
extract_local_sensitivities

Extracts the time series for the local sensitivities from the ODE solution. This requires
that the ODE was defined via `ODEForwardSensitivityProblem`.

```julia
extract_local_sensitivities(sol, asmatrix::Val=Val(false)) # Decompose the entire time series
extract_local_sensitivities(sol, i::Integer, asmatrix::Val=Val(false)) # Decompose sol[i]
extract_local_sensitivities(sol, t::Union{Number,AbstractVector}, asmatrix::Val=Val(false)) # Decompose sol(t)
```
"""
extract_local_sensitivities(sol, asmatrix::Val=Val(false)) = extract_local_sensitivities(sol,sol.prob.problem_type.sensealg, asmatrix)
extract_local_sensitivities(sol, asmatrix::Bool) = extract_local_sensitivities(sol, Val{asmatrix}())
extract_local_sensitivities(sol, i::Integer, asmatrix::Val=Val(false)) = _extract(sol, sol.prob.problem_type.sensealg, sol[i], asmatrix)
extract_local_sensitivities(sol, i::Integer, asmatrix::Bool) = extract_local_sensitivities(sol, i, Val{asmatrix}())
extract_local_sensitivities(sol, t::Union{Number,AbstractVector}, asmatrix::Val=Val(false)) = _extract(sol, sol.prob.problem_type.sensealg, sol(t), asmatrix)
extract_local_sensitivities(sol, t, asmatrix::Bool) = extract_local_sensitivities(sol, t, Val{asmatrix}())
extract_local_sensitivities(tmp, sol, t::Union{Number,AbstractVector}, asmatrix::Val=Val(false)) = _extract(sol, sol.prob.problem_type.sensealg, sol(tmp, t), asmatrix)
extract_local_sensitivities(tmp, sol, t, asmatrix::Bool) = extract_local_sensitivities(tmp, sol, t, Val{asmatrix}())

# Get ODE u vector and sensitivity values from all time points
function extract_local_sensitivities(sol,::ForwardSensitivity, ::Val{false})
  ni = sol.prob.f.numindvar
  u = sol[1:ni, :]
  du = [sol[ni*j+1:ni*(j+1),:] for j in 1:sol.prob.f.numparams]
  return u, du
end

function extract_local_sensitivities(sol,::ForwardDiffSensitivity, ::Val{false})
  u = ForwardDiff.value.(sol)
  du_full = ForwardDiff.partials.(sol)
  firststate = first(du_full)
  firstparam = first(firststate)
  Js = map(1:length(firstparam)) do j
    map(CartesianIndices(du_full)) do II
      du_full[II][j]
    end
  end
  return u, Js
end

function extract_local_sensitivities(sol,::ForwardSensitivity, ::Val{true})
  prob = sol.prob
  ni = prob.f.numindvar
  pn = prob.f.numparams
  jsize = (ni, pn)
  sol[1:ni, :], map(sol.u) do u
    collect(reshape((@view u[ni+1:end]), jsize))
  end
end

function extract_local_sensitivities(sol,::ForwardDiffSensitivity, ::Val{true})
  retu = ForwardDiff.value.(sol)
  jsize = length(sol.u[1]), ForwardDiff.npartials(sol.u[1][1])
  du = map(sol.u) do u
    du_i = similar(retu, jsize)
    for i in eachindex(u)
      du_i[i, :] = ForwardDiff.partials(u[i])
    end
    du_i
  end
  retu, du
end

# Get ODE u vector and sensitivity values from sensitivity problem u vector
function _extract(sol, sensealg::ForwardSensitivity, su::AbstractVector, asmatrix::Val = Val(false))
  u = view(su, 1:sol.prob.f.numindvar)
  du = _extract_du(sol, sensealg, su, asmatrix)
  return u, du
end

function _extract(sol, sensealg::ForwardDiffSensitivity, su::AbstractVector, asmatrix::Val = Val(false))
  u = ForwardDiff.value.(su)
  du = _extract_du(sol, sensealg, su, asmatrix)
  return u, du
end

# Get sensitivity values from sensitivity problem u vector (nested form)
function _extract_du(sol, ::ForwardSensitivity, su::Vector, ::Val{false})
  ni = sol.prob.f.numindvar
  return [view(su, ni*j+1:ni*(j+1)) for j in 1:sol.prob.f.numparams]
end

function _extract_du(sol, ::ForwardDiffSensitivity, su::Vector, ::Val{false})
  du_full = ForwardDiff.partials.(su)
  return [[du_full[i][j] for i in 1:size(du_full,1)] for j in 1:length(du_full[1])]
end

# Get sensitivity values from sensitivity problem u vector (matrix form)
function _extract_du(sol, ::ForwardSensitivity, su::Vector, ::Val{true})
  ni = sol.prob.f.numindvar
  np = sol.prob.f.numparams
  return view(reshape(su, ni, np+1), :, 2:np+1)
end

function _extract_du(sol, ::ForwardDiffSensitivity, su::Vector, ::Val{true})
  du_full = ForwardDiff.partials.(su)
  return [du_full[i][j] for i in 1:size(du_full,1), j in 1:length(du_full[1])]
end

### Bonus Pieces

function SciMLBase.remake(prob::ODEProblem{uType,tType,isinplace,P,F,K,<:ODEForwardSensitivityProblem};
                          f=nothing,tspan=nothing,u0=nothing,p=nothing,kwargs...) where
                          {uType,tType,isinplace,P,F,K}
    _p     = p     === nothing ? prob.p : p
    _f     = f     === nothing ? prob.f.f : f
    _u0    = u0    === nothing ? prob.u0[1:prob.f.numindvar] : u0[1:prob.f.numindvar]
    _tspan = tspan === nothing ? prob.tspan : tspan
    ODEForwardSensitivityProblem(_f,_u0,
                                 _tspan,_p,prob.problem_type.sensealg;
                                 prob.kwargs...,kwargs...)
end
SciMLBase.ODEFunction(f::ODEForwardSensitivityFunction; kwargs...) = f


text:
export arb_series

mutable struct arb_series
    poly::arb_poly
    length::Int
end

###############################################################################
#
#   Constructions
#
###############################################################################

function arb_series(poly::arb_poly)
    arb_series(poly, length(poly))
end

function arb_series(poly)
    arb_series(poly, length(poly))
end

function show(io::IO, x::arb_series)
    len = length(x)
    S = var(parent(x.poly))

    for i in 0:len-1
        if !iszero(x[i])
            if !isone(x[i]) || i == 0
                print(io, x[i])
                print(io, "*")
            end
            if i != 0
                print(io, "$S^$i")
            end
            print(io, " + ")
        end
    end
    print(io, "𝒪($S^$len)")
end

###############################################################################
#
#   Array and iterator interface
#
###############################################################################

function Base.getindex(x::arb_series, i::Int)
    0 <= i <= x.length || throw(BoundsError(x, i))
    return coeff(x.poly, i)
end

Base.getindex(x::arb_series, i::Number) = x[convert(Int, i)]
Base.getindex(x::arb_series, I) = [x[i] for i in I]

Base.eltype(x::arb_series) = arb
Base.iterate(x::arb_series, i = 0) = i > x.length ? nothing : (x[i], i + 1)
Base.indexed_iterate(p::arb_series, i::Int, state=1) = (x[i], i + 1)
Base.broadcastable(x::arb_series) = Ref(x)

function Base.setindex!(x::arb_series, v::arb, i::Int)
    0 <= i <= x.length || throw(BoundsError(x, i))
    setcoeff!(x.poly, i, v)
end

Base.firstindex(x::arb_series) = 0
Base.lastindex(x::arb_series) = x.length - 1

# Basic information
length(x::arb_series) = x.length

# Unary series methods
unary_methods = [# Arithmetic
                 (:inv, "arb_poly_inv_series"),
                 # Transforms
                 (:binomial_transform, "arb_poly_binomial_transform"),
                 # Powers and elementary functions
                 (:sqrt, "arb_poly_sqrt_series"),
                 (:rsqrt, "arb_poly_rsqrt_series"),
                 (:log, "arb_poly_log_series"),
                 (:log1p, "arb_poly_log1p_series"),
                 (:atan, "arb_poly_atan_series"),
                 (:asin, "arb_poly_asin_series"),
                 (:acos, "arb_poly_acos_series"),
                 (:exp, "arb_poly_exp_series"),
                 (:sin, "arb_poly_sin_series"),
                 (:cos, "arb_poly_cos_series"),
                 (:tan, "arb_poly_tan_series"),
                 (:sinpi, "arb_poly_sin_pi_series"),
                 (:cospi, "arb_poly_cos_pi_series"),
                 (:tanpi, "arb_poly_tan_pi_series"),
                 (:sinh, "arb_poly_sinh_series"),
                 (:cosh, "arb_poly_cosh_series"),
                 (:sinc, "arb_poly_sinc_series"),
                 (:sincpi, "arb_poly_sinc_pi_series"),
                 # Gamma function and factorials
                 (:gamma, "arb_poly_gamma_series"),
                 (:rgamma, "arb_poly_rgamma_series"),
                 (:lgamma, "arb_poly_lgamma_series"),
                 (:digamma, "arb_poly_digamma_series"),
                 # Zeta function
                 (:riemann_siegel_theta, "arb_poly_riemann_siegel_theta_series"),
                 (:riemann_siegel_z, "arb_poly_riemann_siegel_z_series"),
                 ]

for (s, f) in unary_methods
    @eval begin
        function ($s)(x::arb_series,
                      n::Integer = length(x))
            y = arb_series(parent(x.poly)(), n)
            ccall(($f, Nemo.libarb), Nothing,
                  (Ref{arb_poly}, Ref{arb_poly}, Int, Int),
                  y.poly, x.poly, n, precision(parent(x.poly)))
            return y
        end
    end
end

# Operators
operator_methods = [# Arithmetic
                    (:+, "arb_poly_add_series"),
                    (:-, "arb_poly_sub_series"),
                    (:*, "arb_poly_mullow"),
                    (:/, "arb_poly_div_series"),
                    # Powers and elementary functions
                    (:^, "arb_poly_pow_series"),
                    ]
for (s, f) in operator_methods
    @eval begin
        function ($s)(x::arb_series,
                      y::arb_series)
            n = max(length(x), length(y))
            z = arb_series(parent(x.poly)(), n)
            ccall(($f, Nemo.libarb), Nothing,
                  (Ref{arb_poly}, Ref{arb_poly}, Ref{arb_poly}, Int, Int),
                  z.poly, x.poly, y.poly, n, precision(parent(x.poly)))
            return z
        end
    end
end

# Binary series methods
binary_methods = [# Arithmetic
                  (:add, "arb_poly_add_series"),
                  (:sub, "arb_poly_sub_series"),
                  (:mul, "arb_poly_mullow"),
                  (:div, "arb_poly_div_series"),
                  # Composition
                  (:compose, "arb_poly_compose_series"),
                  (:revert, "arb_poly_rever_series"),
                  # Powers and elementary functions
                  (:pow, "arb_poly_pow_series"),
                  ]
for (s, f) in binary_methods
    @eval begin
        function ($s)(x::arb_series,
                      y::arb_series,
                      n::Integer = max(length(x), length(y)))
            z = arb_series(parent(x.poly)(), n)
            ccall(($f, Nemo.libarb), Nothing,
                  (Ref{arb_poly}, Ref{arb_poly}, Ref{arb_poly}, Int, Int),
                  z.poly, x.poly, y.poly, n, precision(parent(x.poly)))
            return z
        end
    end
end

# Unary length-preserving non-series methods
unary_length_preserving_methods = [(:-, "arb_poly_neg"),
                                   (:borel_transform, "arb_poly_borel_transform"),
                                   (:inv_borel_transform, "arb_poly_inv_borel_transform"),
                                   ]

for (s, f) in unary_length_preserving_methods
    @eval begin
        function ($s)(x::arb_series)
            y = arb_series(parent(x.poly)(), length(x))
            ccall(($f, Nemo.libarb), Nothing,
                  (Ref{arb_poly}, Ref{arb_poly}, Int),
                  y.poly, x.poly, precision(parent(x.poly)))
            return y
        end
    end
end

# Unary non-length-preserving methods
unary_non_length_preserving_methods = [(:derivative, "arb_poly_derivative"),
                                       (:integral, "arb_poly_integral"),
                                       ]

function derivative(x::arb_series,
                    n::Integer = max(length(x) - 1, 0))
    y = arb_series(parent(x.poly)(), n)
    ccall(("arb_poly_derivative", Nemo.libarb), Nothing,
          (Ref{arb_poly}, Ref{arb_poly}, Int),
          y.poly, x.poly, precision(parent(x.poly)))
    return y
end

function integral(x::arb_series,
                  n::Integer = length(x) + 1)
    y = arb_series(parent(x.poly)(), n)
    ccall(("arb_poly_integral", Nemo.libarb), Nothing,
          (Ref{arb_poly}, Ref{arb_poly}, Int),
          y.poly, x.poly, precision(parent(x.poly)))
    return y
end

# Special
(:^, "arb_poly_pow_ui_trunc_binexp")
(:^, "arb_poly_pow_arb_series")

(:sincos, "arb_poly_sin_cos_series")
function sincos(x::arb_series,
                n::Integer = length(x))
    y = arb_series(parent(x.poly)(), n)
    z = arb_series(parent(x.poly)(), n)
    ccall((:arb_poly_sin_cos_series, Nemo.libarb), Nothing,
          (Ref{arb_poly}, Ref{arb_poly}, Ref{arb_poly}, Int, Int),
          y.poly, z.poly, x.poly, n, precision(parent(x.poly)))
    return (y, z)
end

function atan(y::arb_series,
              x::arb_series,
              n::Integer = max(length(x), length(y)))
    xdy = mul(x, derivative(y), n - 1)
    ydx = mul(y, derivative(x), n - 1)
    x2 = mul(x, x, n - 1)
    y2 = mul(y, y, n - 1)

    z = integral((xdy - ydx)/(x2 + y2), n)
    if n > 0
        z[0] = atan(y[0], x[0])
    end
    return z
end

function zeta(s::arb_series,
              a::arb,
              n::Integer = length(s))
    y = arb_series(parent(s.poly)(), n)
    ccall((:arb_poly_zeta_series, Nemo.libarb), Nothing,
          (Ref{arb_poly}, Ref{arb_poly}, Ref{arb}, Int32, Int, Int),
          y.poly, s.poly, a, 0, n, precision(parent(s.poly)))
    return y
end
                 (:zeta, "arb_poly_zeta_series")

###############################################################################
#
#   Promotions
#
###############################################################################

promote_rule(::Type{arb_series}, ::Type{Float64}) = arb_series

promote_rule(::Type{arb_series}, ::Type{BigFloat}) = arb_series

promote_rule(::Type{arb_series}, ::Type{fmpz}) = arb_series

promote_rule(::Type{arb_series}, ::Type{fmpq}) = arb_series

promote_rule(::Type{arb_series}, ::Type{arb}) = arb_series

promote_rule(::Type{arb_series}, ::Type{fmpz_poly}) = arb_series

promote_rule(::Type{arb_series}, ::Type{fmpq_poly}) = arb_series

promote_rule(::Type{arb_series}, ::Type{T}) where {T <: Integer} = arb_series

promote_rule(::Type{arb_series}, ::Type{Rational{T}}) where T <: Union{Int, BigInt} = arb_series

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

for T in [Integer, fmpz, fmpq, Float64, BigFloat, arb, fmpz_poly, fmpq_poly, arb_poly]
    @eval begin
        +(x::arb_series, y::$T) = x + arb_series(parent(x.poly)(y))

        +(x::$T, y::arb_series) = y + x

        -(x::arb_series, y::$T) = x - arb_series(parent(x.poly)(y))

        -(x::$T, y::arb_series) = arb_series(parent(y.poly)(x)) - y

        *(x::arb_series, y::$T) = x * arb_series(parent(x.poly)(y))

        *(x::$T, y::arb_series) = y * x

        /(x::arb_series, y::$T) = x / arb_series(parent(x.poly)(y))

        /(x::$T, y::arb_series) = arb_series(parent(y.poly)(x)) / y

        ^(x::arb_series, y::$T) = x^arb_series(parent(x.poly)(y))

        ^(x::$T, y::arb_series) = arb_series(parent(y.poly)(x))^y
    end
end

+(x::arb_series, y::Rational{T}) where T <: Union{Int, BigInt} = x + arb_series(parent(x.poly)(y))

+(x::Rational{T}, y::arb_series) where T <: Union{Int, BigInt} = y + x

-(x::arb_series, y::Rational{T}) where T <: Union{Int, BigInt} = x - arb_series(parent(x.poly)(y))

-(x::Rational{T}, y::arb_series) where T <: Union{Int, BigInt} = arb_series(parent(y.poly)(x)) - y

*(x::arb_series, y::Rational{T}) where T <: Union{Int, BigInt} = x * arb_series(parent(x.poly)(y))

*(x::Rational{T}, y::arb_series) where T <: Union{Int, BigInt} = y * x

/(x::arb_series, y::Rational{T}) where T <: Union{Int, BigInt} = x / arb_series(parent(x.poly)(y))

/(x::Rational{T}, y::arb_series) where T <: Union{Int, BigInt} = arb_series(parent(y.poly)(x)) / y

^(x::arb_series, y::Rational{T}) where T <: Union{Int, BigInt} = x^arb_series(parent(x.poly)(y))

^(x::Rational{T}, y::arb_series) where T <: Union{Int, BigInt} = arb_series(parent(y.poly)(x))^y

ext: jl


text:
using CGS
# Period estimates for Dan's sample assuming e=0:
ncand = 16
KIC = ["3218908","3239945","4754460","6551440","8410697","8426957","8505215","8738735","8800954","9306307","10187159","10287723","10321319","10602068","10842718","11709124"]
b = [0.,0.207,0.893,0.75,0.,0.889,0.,0.,0.,0.64,0.,0.8,0.7,0.6027,0.,0.55]
rho = [2.716,2.963,0.8725,1.053,1.201,1.145,3.139,1.388,2.502,1.383,1.553,2.64,1.414,1.762,1.177,1.377]
srho1 = [1.432,0.4687,0.476,0.6332,0.6885,0.7186,0.6771,0.5794,0.4743,0.9142,1.3,0.2721,0.8579,1.09,0.6849,0.5512]
srho2 = [0.3974,0.2218,0.7728,0.5275,0.7425,0.4541,0.2093,0.1768,0.3012,1.01,0.6905,0.4805,0.7137,0.3954,0.5948,0.465]
tdur = [21.45,16.202,15.92,10.85,19.77,39.4,20.06,27.44,15.76,8.499,11.81,9.49,16.84,12.804,35.92,17.75]
rad = [0.514,0.876,0.67,0.282,0.698,1.04,0.277,0.355,0.386,1.22,0.43,0.266,0.163,2.0,0.74,0.83]*RJUPITER/REARTH
srad = [0.092,0.039,0.16,0.093,0.107,0.30,0.017,0.045,0.025,0.49,0.21,0.027,0.046,0.66,0.16,0.12]
#rho  = 3/pi^2*(period*24.*3600.)/GRAV/(tdur*3600.)^3*(1.-b^2)^1.5
fac = pi^2/3/24./3600.*3600.^3*GRAV
p = fac.*rho.*tdur.^3./(1-b.^2).^1.5
p1 = p.*(1.-srho1./rho)
p2 = p.*(1.+srho2./rho)
p[2]=2.9328721*365.25
p1[2]=(2.9328721-2.6e-6)*365.25
p2[2]=(2.9328721+2.6e-6)*365.25
p[5]=2.8688097*365.25
p1[5]=(2.8688097-5.3e-6)*365.25
p2[5]=(2.8688097+5.4e-6)*365.25
p[9]=1.9279957*365.25
p1[9]=(1.9279957-9.1e-6)*365.25
p2[9]=(1.9279957+9.2e-6)*365.25
using PyPlot
clf()
for i=1:ncand
  loglog([p1[i],p2[i]],[rad[i],rad[i]])
  loglog([p[i],p[i]],[rad[i]-srad[i],rad[i]+srad[i]])
#  println(KIC[i]," ",p1[i]/365.25," ",p2[i]/365.25," ",rad[i])
  println(KIC[i]," ",p[i]/365.25," ",rad[i])
end
#plot([780,780],[1,100])
#plot([4*365,4*365],[1,100])
plot([600,600],[1,100])
plot([1000,1000],[1,100])

# Notes on individual systems:
# KIC 3218908 - high probability of being a planet, and inner planets agree in duration/period
#   diagram.  The parameters look reasonable, so this is a likely candidate. [High quality]
#
# KIC 3239945 - high FAP (27%); but it has multiple transits, and is consistent with inner planets.
#   The parameters look well determined.  This is a high quality candidate. [High quality]
#
# KIC 4754460 - Likely eclipse, high b, so probably FP. [False alarm]
#
# KIC 6551440 - Has a most likely period that is quite short.  High b, broad eccentricity.
#  FAP is modest 3%.  Perhaps stellar parameters are off? [Stellar params?]
#
# KIC 8410697 - Multiple transits, so well determined period.  Parameters look well constrained.
#  This is a high quality candidate. [High quality]
#
# KIC 8426957 - FAP is high 20%; high b.  Period is *very* long, which seems unlikely.  Perhaps stellar
#  parameters are off. [Stellar params?]
#
# KIC 8505215 - There is a short-period candidate 99.02 which is labeled as a false-positive (may be a blend).
#  DV summary shows an offset during transit. Dan's candidate looks good, though. [Blend?]
#
# KIC 8738735 - Stellar density inferred from inner two planets is inconsistent with Huber's density.
#  Otherwise the parameters look reasonable. [Stellar params?]
#
# KIC 8800954 - Multiple transits, parameters look good. [High quality]
#
# KIC 9306307 - High impact parameter (multi-modal), high eccentricity. Inferred period for e=0
#  is quite short.  Perhaps stellar parameters are off? [Stellar params?]
#
# KIC 10187159 - Density inferred from inner planet is much higher than Huber's.  Parameters of fit look
#  fine, but peaked at very short period. FAP is 9%. [Stellar params?]
#
# KIC 10287723 - High impact parameter, high eccentricity.  Period peaks at ~few years.  Low FAP (5%).
#  [High quality]
#
# KIC 10321319 - Low FAP; parameters look good. [High quality]
#
# KIC 10602068 - High FAP (~100%). [False alarm]
#
# KIC 10842718 - 10% FAP, but parameters look good. [High quality]
#
# KIC 11709124 - Parameters look good. [High quality]

max_stars_repo_name: dfm/single-transits

max_issues_repo_name: dfm/single-transits

max_forks_repo_name: dfm/single-transits


text:
run_task() = _run_task
isdone() = _run_task.state == :done
interrupt_task() = @async Base.throwto(_run_task, InterruptException())

"""
    serialized_value(x) = nothing

Allow types to be serialized back to GtkREPL, by default sends `nothing`.
GtkREPL will call `display` on these objects.
"""
serialized_value(x) = nothing

# FIXME dirty hack
# see https://github.com/JunoLab/Atom.jl/blob/52ef77a93605cb4ddbf455ba36d1d7e9ae637edc/src/eval.jl#L129
function clean_error_msg(s::String)
    r  = r"^(.*)\s\[\d\] _eval_command_remotely.*$"s
    m = match(r, s)
    !isnothing(m) && return m.captures[1]
    s
end

function trim(s::AbstractString, L::Int)#need to be AbstracString to accept SubString
    if length(s) > L
        return string(s[1:L], "...")
    end
    s
end

function format_output(x)
    io = IOBuffer()
    io = IOContext(io, :display_size=>(20, 20))
    io = IOContext(io, :limit=>true)
    show(io, MIME"text/plain"(), x)
    String(take!(io.io))
end

function eval_command_remotely(cmd::String, eval_in::String)
    @async begin
        mod =  @eval Main $(Meta.parse(eval_in))
        t = @elapsed finalOutput, v = _eval_command_remotely(cmd, mod)
        #TODO check stdout at this point ? Let's wait a bit at least
        sleep(50/1000)
        lock(serialize_lock) do
            serialize(gtkrepl.socket, (EvalDone(
                gtkrepl.console_idx,
                (finalOutput,v),
                t
            ),))
        end
    end
    true
end

function eval_shell_remotely(cmd::String, eval_in::String)
    @async begin
        t = @elapsed finalOutput, v = _eval_shell_remotely(cmd)
        #TODO merge this with eval_command_remotely ?
        lock(serialize_lock) do
            serialize(gtkrepl.socket, (EvalDone(
                gtkrepl.console_idx,
                (finalOutput,v),
                t
            ),))
        end
    end
    true
end

function _eval_command_remotely(cmd::String, eval_in::Module)

    evalout = ""
    v = :()
    try
        v = include_string(eval_in, cmd)
        Core.eval(eval_in, :(ans = $(Expr(:quote, v))))
        evalout = isnothing(v) ? "" : format_output(v)
    catch err
        bt = catch_backtrace()
        evalout = clean_error_msg( sprint(showerror, err, bt) )
    end

    evalout = trim(evalout, 4000)
    finalOutput = evalout == "" ? "" : "$evalout\n"

    v = serialized_value(v)

    return finalOutput, v
end

function _eval_shell_remotely(cmd::String)
    evalout = try
        cmd = Core.eval(Main, :(Base.cmd_gen($(Base.shell_parse(cmd)[1]))) )
        evalout = read(cmd, String)
    catch err
        bt = catch_backtrace()
        evalout = clean_error_msg( sprint(showerror, err, bt) )
    end
    return evalout, nothing
end

function remote_completions(cmd)
    comp, dotpos = completions(cmd, lastindex(cmd))
    comp = REPL.completion_text.(comp)# avoid returning the object itself https://github.com/JuliaLang/julia/issues/33747
    comp, dotpos
end


text:
###################################################
## geometry.jl
## Geometry tools
###################################################

"""
    Return bounding box of polygon/ list of points
"""
function boundingBox{T<:AbstractFloat}(poly::Vector{Tuple{T,T}})
    minX::T = Inf; maxX::T = -Inf; minY::T = Inf; maxY::T = -Inf
    for (x,y) in poly
        minX = min(minX,x)
        maxX = max(maxX,x)
        minY = min(minY,y)
        maxY = max(maxY,y)
    end
    return (minX,maxX,minY,maxY)
end

"""
    Scales a given box (i.e. 5% bigger if scale=1.05)
"""
function scaleBox(box, scale)
    minX, maxX, minY, maxY = box
    lenX = maxX - minX; midX = (minX + maxX)/2
    lenY = maxY - minY; midY = (minY + maxY)/2
    return (midX - lenX*scale/2, midX + lenX*scale/2, midY-lenY*scale/2, midY+lenY*scale/2)
end

"""
    Distance between the two points in meters (true geodesic distance from LLA coordinates)
"""
distanceGeo(pt1::Node, pt2::Node) = distanceGeo(pt1.lon, pt1.lat, pt2.lon, pt2.lat)

function distanceGeo(lon1, lat1, lon2, lat2)
    dLat = toradians(lat2 - lat1)
    dLon = toradians(lon2 - lon1)
    lat1 = toradians(lat1)
    lat2 = toradians(lat2)
    a = sin(dLat/2)^2 + sin(dLon/2)^2 * cos(lat1) * cos(lat2)
    2.0 * atan2(sqrt(a), sqrt(1-a)) * 6373.0 * 1000
end

"""
    Distance between the two points in meters (from coordinates)
"""
distanceCoord(pt1::Node, pt2::Node) = distanceCoord(pt1.x, pt1.y, pt2.x, pt2.y)

function distanceCoord{T<:AbstractFloat}(x1::T, y1::T, x2::T, y2::T)
    return sqrt((x2-x1)^2+(y2-y1)^2)
end

"""
    Check is a point is inside a polygon
"""
function pointInsidePolygon{T<:AbstractFloat}(x::AbstractFloat,y::AbstractFloat,poly::Vector{Tuple{T,T}})
    n = length(poly)
    inside =false

    p1x,p1y = poly[1]
    for i in 0:n
        p2x,p2y = poly[i % n + 1]
        if y > min(p1y,p2y) && y <= max(p1y,p2y) && x <= max(p1x,p2x)
            if p1y != p2y
                xinters = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
            end
            if p1x == p2x || x <= xinters
                inside = !inside
            end
        end
        p1x,p1y = p2x,p2y
    end
    return inside
end

"""
    `toENU`, projects latitude and longitude to ENU coordinate system
"""
function toENU(lon, lat, center)
    enu = ENU(LLA(lat,lon), LLA(center[2],center[1]))
    return enu.east, enu.north
end
toENU(lon, lat, n::Network) = toENU(lon, lat, n.projcenter)
"""
    `updateProjection`
    updates the (x,y) projection of the nodes of the network given their latitude/longitude
"""
function updateProjection!(n::Network)
    bounds = boundingBox([(n.lon,n.lat) for n in n.nodes])
    n.projcenter = ((bounds[2]+bounds[1])/2, (bounds[4]+bounds[3])/2)
    nodes = Array{Node}(length(n.nodes))
    for (i,no) in enumerate(n.nodes)
        x,y = toENU(no.lon,no.lat,n.projcenter)
        nodes[i] = Node(x,y,no.lon,no.lat)
    end
    n.nodes = nodes
    n
end

max_stars_repo_name: andrew-zhm/RoutingNetworks

max_issues_repo_name: andrew-zhm/RoutingNetworks

max_forks_repo_name: andrew-zhm/RoutingNetworks


text:
################################################################################
#
#   PYROCHLORE LATTICE
#
################################################################################

# Implementation
# - implementation 1
# - labels <: Any
function getUnitcellPyrochlore(
            unitcell_type  :: Type{U},
            implementation :: Val{1}
        ) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

    # return a new Unitcell
    return newUnitcell(
        # Type of the unitcell
        U,
        # lattice vectors
        Vector{Float64}[
            Float64[0, 0.5, 0.5],
            Float64[0.5, 0, 0.5],
            Float64[0.5, 0.5, 0]
        ],
        # sites
        S[
            newSite(S, Float64[0., 0., 0.], getDefaultLabelN(LS,1)),
            newSite(S, Float64[0., 0.25, 0.25], getDefaultLabelN(LS,2)),
            newSite(S, Float64[0.25, 0., 0.25], getDefaultLabelN(LS,3)),
            newSite(S, Float64[0.25, 0.25, 0.], getDefaultLabelN(LS,4))
        ],
        # bonds
        B[
            newBond(B, 1,2, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 1,3, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 1,4, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 2,1, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 2,3, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 2,4, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 3,1, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 3,2, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 3,4, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 4,1, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 4,2, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 4,3, getDefaultLabel(LB), (0,0,0)),

            newBond(B, 1,4, getDefaultLabel(LB), (0,0,-1)),
            newBond(B, 4,1, getDefaultLabel(LB), (0,0,+1)),
            newBond(B, 1,2, getDefaultLabel(LB), (-1,0,0)),
            newBond(B, 2,1, getDefaultLabel(LB), (+1,0,0)),
            newBond(B, 1,3, getDefaultLabel(LB), (0,-1,0)),
            newBond(B, 3,1, getDefaultLabel(LB), (0,+1,0)),

            newBond(B, 2,3, getDefaultLabel(LB), (+1,-1,0)),
            newBond(B, 3,2, getDefaultLabel(LB), (-1,+1,0)),
            newBond(B, 2,4, getDefaultLabel(LB), (+1,0,-1)),
            newBond(B, 4,2, getDefaultLabel(LB), (-1,0,+1)),
            newBond(B, 3,4, getDefaultLabel(LB), (0,+1,-1)),
            newBond(B, 4,3, getDefaultLabel(LB), (0,-1,+1))
        ]
    )
end

# Implementation
# - implementation 2
# - labels <: Any
function getUnitcellPyrochlore(
            unitcell_type  :: Type{U},
            implementation :: Val{2}
        ) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

    # return a new Unitcell
    return newUnitcell(
        # Type of the unitcell
        U,
        # lattice vectors
        Vector{Float64}[
            Float64[1., 0, 0],
            Float64[0, 1., 0],
            Float64[0, 0, 1.]
        ],
        # sites
        S[
            newSite(S, Float64[0.25, 0, 0], getDefaultLabelN(LS,1)),
            newSite(S, Float64[0, 0.25, 0], getDefaultLabelN(LS,2)),
            newSite(S, Float64[0, 0, 0.25], getDefaultLabelN(LS,3)),
            newSite(S, Float64[0.25, 0.25, 0.25], getDefaultLabelN(LS,4)),

            newSite(S, Float64[0.75, 0.5, 0], getDefaultLabelN(LS,5)),
            newSite(S, Float64[0.5, 0.75, 0], getDefaultLabelN(LS,6)),
            newSite(S, Float64[0.5, 0.5, 0.25], getDefaultLabelN(LS,7)),
            newSite(S, Float64[0.75, 0.75, 0.25], getDefaultLabelN(LS,8)),

            newSite(S, Float64[0.75, 0, 0.5], getDefaultLabelN(LS,9)),
            newSite(S, Float64[0.5, 0.25, 0.5], getDefaultLabelN(LS,10)),
            newSite(S, Float64[0.5, 0, 0.75], getDefaultLabelN(LS,11)),
            newSite(S, Float64[0.75, 0.25, 0.75], getDefaultLabelN(LS,12)),

            newSite(S, Float64[0.25, 0.5, 0.5], getDefaultLabelN(LS,13)),
            newSite(S, Float64[0, 0.75, 0.5], getDefaultLabelN(LS,14)),
            newSite(S, Float64[0, 0.5, 0.75], getDefaultLabelN(LS,15)),
            newSite(S, Float64[0.25, 0.75, 0.75], getDefaultLabelN(LS,16))
        ],
        # bonds
        B[
            newBond(B, 1,2, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 1,3, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 1,4, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 2,1, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 2,3, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 2,4, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 3,1, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 3,2, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 3,4, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 4,1, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 4,2, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 4,3, getDefaultLabel(LB), (0,0,0)),

            newBond(B, 5,6, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 5,7, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 5,8, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 6,5, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 6,7, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 6,8, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 7,5, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 7,6, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 7,8, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 8,5, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 8,6, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 8,7, getDefaultLabel(LB), (0,0,0)),

            newBond(B, 9,10, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 9,11, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 9,12, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 10,9, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 10,11, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 10,12, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 11,9, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 11,10, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 11,12, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 12,9, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 12,10, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 12,11, getDefaultLabel(LB), (0,0,0)),

            newBond(B, 13,14, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 13,15, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 13,16, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 14,13, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 14,15, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 14,16, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 15,13, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 15,14, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 15,16, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 16,13, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 16,14, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 16,15, getDefaultLabel(LB), (0,0,0)),

            newBond(B, 4,7, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 4,10, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 4,13, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 7,4, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 7,10, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 7,13, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 10,4, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 10,7, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 10,13, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 13,4, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 13,7, getDefaultLabel(LB), (0,0,0)),
            newBond(B, 13,10, getDefaultLabel(LB), (0,0,0)),

            newBond(B, 8,3, getDefaultLabel(LB), (1,1,0)),
            newBond(B, 3,8, getDefaultLabel(LB), (-1,-1,0)),
            newBond(B, 8,14, getDefaultLabel(LB), (1,0,0)),
            newBond(B, 14,8, getDefaultLabel(LB), (-1,0,0)),
            newBond(B, 14,3, getDefaultLabel(LB), (0,1,0)),
            newBond(B, 3,14, getDefaultLabel(LB), (0,-1,0)),
            newBond(B, 8,9, getDefaultLabel(LB), (0,1,0)),
            newBond(B, 9,8, getDefaultLabel(LB), (0,-1,0)),
            newBond(B, 14,9, getDefaultLabel(LB), (-1,1,0)),
            newBond(B, 9,14, getDefaultLabel(LB), (1,-1,0)),
            newBond(B, 3,9, getDefaultLabel(LB), (-1,0,0)),
            newBond(B, 9,3, getDefaultLabel(LB), (1,0,0)),

            newBond(B, 12,15, getDefaultLabel(LB), (1,0,0)),
            newBond(B, 15,12, getDefaultLabel(LB), (-1,0,0)),
            newBond(B, 12,2, getDefaultLabel(LB), (1,0,1)),
            newBond(B, 2,12, getDefaultLabel(LB), (-1,0,-1)),
            newBond(B, 15,2, getDefaultLabel(LB), (0,0,1)),
            newBond(B, 2,15, getDefaultLabel(LB), (0,0,-1)),
            newBond(B, 12,5, getDefaultLabel(LB), (0,0,1)),
            newBond(B, 5,12, getDefaultLabel(LB), (0,0,-1)),
            newBond(B, 2,5, getDefaultLabel(LB), (-1,0,0)),
            newBond(B, 5,2, getDefaultLabel(LB), (1,0,0)),
            newBond(B, 15,5, getDefaultLabel(LB), (-1,0,1)),
            newBond(B, 5,15, getDefaultLabel(LB), (1,0,-1)),

            newBond(B, 16,11, getDefaultLabel(LB), (0,1,0)),
            newBond(B, 11,16, getDefaultLabel(LB), (0,-1,0)),
            newBond(B, 6,11, getDefaultLabel(LB), (0,1,-1)),
            newBond(B, 11,6, getDefaultLabel(LB), (0,-1,1)),
            newBond(B, 16,6, getDefaultLabel(LB), (0,0,1)),
            newBond(B, 6,16, getDefaultLabel(LB), (0,0,-1)),
            newBond(B, 16,1, getDefaultLabel(LB), (0,1,1)),
            newBond(B, 1,16, getDefaultLabel(LB), (0,-1,-1)),
            newBond(B, 1,6, getDefaultLabel(LB), (0,-1,0)),
            newBond(B, 6,1, getDefaultLabel(LB), (0,1,0)),
            newBond(B, 1,11, getDefaultLabel(LB), (0,0,-1)),
            newBond(B, 11,1, getDefaultLabel(LB), (0,0,1))
        ]
    )
end


text:
# ## [Пример №4, краткий, динамическая сетка](@id de_4)
#
# Делаем тоже самое, что в [развернутом примере №3](@ref de_3),
# только с использованием функций сокращений и стандартных параметров.
using NonLinearReactionAdvectionDiffusionWithFrontData
using NonLinearReactionAdvectionDiffusionWithFrontData: shishkin_mesh;

## Задача параметров и решение в коротком виде
a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀, mshfrm_std = NonLinearReactionAdvectionDiffusionWithFrontData.dparams_nonuniform();
u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ, create_mesh = mshfrm_std);
nothing #hide

## Генерация априорной информации
ϕl      = phidetermination(qₙ, ulₘ, Xₙ, N, Tₘ, M);                  # Левый вырожденный корень
ϕr      = phidetermination(qₙ, urₘ, Xₙ, N, Tₘ, M, reverseX = true); # Правый вырожденный корень
ϕ       = Φ(ϕl, ϕr, N, M);                                          # Полуразность вырожденных корней
ϕ       = apply_on_dynamic_mesh(ϕ, XX, N, M);                       # Аппроксимация на переменную сетку
ϕl      = apply_on_dynamic_mesh(ϕl, XX, N, M);                      # Аппроксимация на переменную сетку
ϕr      = apply_on_dynamic_mesh(ϕr, XX, N, M);                      # Аппроксимация на переменную сетку
f1_data = f1(ϕ, u, XX, N, M);                                       # Положение переходного слоя
f2_data = f2(f1_data, u, XX, N, M);                                 # Значение функции на переходном слое
nothing # hide


text:
using Documenter, Literate, NCTiles, Pkg

pth=@__DIR__
lst=("Example1.jl","Example2.jl","Example3.jl","Example4.jl","Example5.jl","Example6.jl")
lstExecute=("Example1.jl","Example2.jl","Example3.jl","Example4.jl","Example5.jl","Example6.jl")
for i in lst
    EXAMPLE = joinpath(pth, "..", "examples", i)
    OUTPUT = joinpath(pth, "src","generated")
    Pkg.activate(joinpath(pth,"..","docs"))
    Literate.markdown(EXAMPLE, OUTPUT, documenter = true)
    cd(pth)
    Pkg.activate(joinpath(pth,"..","docs"))
    tmp=xor(occursin.(i,lstExecute)...)
    Literate.notebook(EXAMPLE, OUTPUT, execute = tmp)
    cd(pth)
    #Literate.notebook(EXAMPLE, OUTPUT, flavor = :pluto)
end

makedocs(
    sitename = "NCTiles",
    format = Documenter.HTML(),
   pages=[
        "Home" => "index.md",
        "Examples" => Any[
            "Guide " => "examples.md",
            "Listing" => [map(s -> "generated/$(s[1:end-2])md",lst)...],
            ],
        "maindocs.md",
        "API.md"],
    modules = [NCTiles]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/gaelforget/NCTiles.jl.git",
)


text:
module AutoGrad
export Param, params, grad, value, @diff
export gradloss, getval
export @primitive, @zerograd, @primitive1, @zerograd1, @primitive2, @zerograd2
export @gcheck

# Set ENV["AUTOGRAD_TIMER"]="true" and Pkg.build("AutoGrad") if you want profiling information in AutoGrad.to
using TimerOutputs, Libdl
const TIMER=haskey(ENV,"AUTOGRAD_TIMER")
const to = TimerOutput()
macro gs(); if !isempty(Libdl.find_library("libcudart")); esc(:(ccall(("cudaDeviceSynchronize","libcudart"),UInt32,()))); end; end
macro timer(name,expr); TIMER ? :(@timeit to $(esc(name)) (a=$(esc(expr));@gs;a)) : esc(expr); end

"""
Usage:

    x = Param([1,2,3])          # The user declares parameters with `Param`
    y = @diff sum(x .* x)       # computes gradients using `@diff`
    grad(y,x) => [2,4,6]        # looks up the gradient of a parameter with `grad`

`Param(x)` returns a struct that acts like `x` but marks it as a parameter you want to
compute gradients with respect to.

`@diff expr` evaluates an expression and returns a struct that contains its value (which
should be a scalar) and gradients with respect to the `Param`s used in the computation.

`grad(y, x)` returns the gradient of a `@diff` result `y` with respect to any parameter
`x::Param`. (`nothing` may be returned if the gradient is 0).

`value(x)` returns the value associated with `x` if `x` is a `Param` or the output of
`@diff`, otherwise returns `x`.

`params(x)` returns an iterator of `Param`s found by a recursive search of object `x`, which
is typically a model or a `@diff` result.

Alternative usage:

    x = [1 2 3]
    f(x) = sum(x .* x)
    f(x) => 14
    grad(f)(x) => [2 4 6]
    gradloss(f)(x) => ([2 4 6], 14)

Given a scalar valued function `f`, `grad(f,argnum=1)` returns another function `g` which
takes the same inputs as `f` and returns the gradient of the output with respect to the
argnum'th argument. `gradloss` is similar except the resulting function also returns f's
output.
"""
AutoGrad, Param, params, :(@diff), grad, value, gradloss

using LinearAlgebra, Statistics, SpecialFunctions
include("core.jl")
include("macros.jl")
include("sparse.jl")
include("getindex.jl")
include("iterate.jl")
include("addto.jl")
include("unbroadcast.jl")
include("base.jl")
include("math.jl")
include("statistics.jl")
include("linearalgebra.jl")
include("specialfunctions.jl")
include("cat.jl")
include("show.jl")
include("params.jl")
include("../test/gradcheck.jl")

"Use `AutoGrad.dir(path...)` to construct paths relative to AutoGrad root."
dir(path...) = joinpath(dirname(dirname(@__FILE__)),path...)

# Deprecations:
@deprecate Rec Value
@deprecate getval value

end # module


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using UnitCommitment, Test, LinearAlgebra
import UnitCommitment: _Violation, _offer, _query

@testset "_ViolationFilter" begin
    instance = UnitCommitment.read_benchmark("test/case14")
    filter = UnitCommitment._ViolationFilter(max_per_line = 1, max_total = 2)

    _offer(
        filter,
        _Violation(
            time = 1,
            monitored_line = instance.lines[1],
            outage_line = nothing,
            amount = 100.0,
        ),
    )
    _offer(
        filter,
        _Violation(
            time = 1,
            monitored_line = instance.lines[1],
            outage_line = instance.lines[1],
            amount = 300.0,
        ),
    )
    _offer(
        filter,
        _Violation(
            time = 1,
            monitored_line = instance.lines[1],
            outage_line = instance.lines[5],
            amount = 500.0,
        ),
    )
    _offer(
        filter,
        _Violation(
            time = 1,
            monitored_line = instance.lines[1],
            outage_line = instance.lines[4],
            amount = 400.0,
        ),
    )
    _offer(
        filter,
        _Violation(
            time = 1,
            monitored_line = instance.lines[2],
            outage_line = instance.lines[1],
            amount = 200.0,
        ),
    )
    _offer(
        filter,
        _Violation(
            time = 1,
            monitored_line = instance.lines[2],
            outage_line = instance.lines[8],
            amount = 100.0,
        ),
    )

    actual = _query(filter)
    expected = [
        _Violation(
            time = 1,
            monitored_line = instance.lines[2],
            outage_line = instance.lines[1],
            amount = 200.0,
        ),
        _Violation(
            time = 1,
            monitored_line = instance.lines[1],
            outage_line = instance.lines[5],
            amount = 500.0,
        ),
    ]
    @test actual == expected
end


text:
function getlane(file)
    lane = match(r"L[0-9][0-9][0-9]", file)
    if isnothing(lane)
        return -1
    else
        return parse(Int, lane.match[2:end])
    end
end

function rundemultiplex(readfile, indexfile, metafile, outfolder, ignorelane=true, bl=6)

    println("[SFQ]\tignorelane         :\t", ignorelane)

    meta = CSV.File(metafile, delim='\t')
    indexes = map(strip, meta.Index)
    samples = map(strip, meta.SampleName)
    indexlength = maximum(length, indexes)

    if !all(indexlength .== length.(indexes))
        error("Not all indexes identical length: \nindexes:$indexes \nindexlengths: $(length.(indexes))")
    end

    if hasproperty(meta, :BarcodeLength)
        barcodelengths = meta.BarcodeLength
    else
        barcodelengths = fill(bl, length(indexes))
        println("[SFQ]\tBarcode lengths not specified, using $bl for all samples")
    end

    if !ignorelane && !hasproperty(meta, :Lane)
        println("[SFQ]\tWarning ignorelane=false and lane not set in metafile, ignoring lane")
        ignorelane = true
    end

    if ignorelane
        demultiplex(readfile, indexfile, indexes, samples, outfolder, barcodelengths, indexlength)
    else
        lane = getlane(readfile)
        ind = meta.Lane .== lane
        println("[SFQ]\tRestricting Lane   :\t", lane, " containing ", sum(ind), " samples, excluding ", sum(.!ind), " samples.")
        demultiplex(readfile, indexfile, indexes[ind], samples[ind], outfolder, barcodelengths[ind], indexlength)
    end

end

function demultiplex(readfile, indexfile, indexes, labels, folder, barcodelengths, il=6)

    println("[SFQ]\tSplitting Reads    :\t$readfile")
    println("[SFQ]\tSplitting Indexes  :\t$indexfile")
    println("[SFQ]\tIndexes            :\t$indexes")
    println("[SFQ]\tLabels             :\t$labels")
    println("[SFQ]\tBarcodeLengths     :\t$barcodelengths")
    println("[SFQ]\tFolder             :\t$folder")

    ### open files
    outdirs = joinpath.(folder, string.("Sample_", labels))
    for d in outdirs
        try
            mkpath(d)
        catch
        end
    end

    files   = joinpath.(outdirs, string.(labels, "_", indexes, "_", basename.(readfile)))

    um_file = joinpath(folder, "Sample_unmatched", "unmatched_index_"*basename(readfile))

    mkpath(dirname(um_file))
    files   = [files ; um_file]
    println("[SFQ]\tFiles            :\t$files")

    streams = GzipCompressorStream.(open.(files, "w"))
    index_dict = Dict(indexes[i] => i for i = 1:length(indexes))
    counts = zeros(Int, length(streams))

    rio = open(readfile)  |> GzipDecompressorStream
    iio = open(indexfile) |> GzipDecompressorStream

    total_reads = 0
    st = time()
    while !eof(rio) && !eof(iio)
        read_id   = readline(rio)
        read_read = readline(rio)
        read_qid  = readline(rio)
        read_qs   = readline(rio)

        index_id   = readline(iio)
        index_read = readline(iio)
        index_qid  = readline(iio)
        index_qs   = readline(iio)

        index, barcode = index_read[1:il], index_read[(il+1):end]
        stream_index = get(index_dict, index, length(files))

        if stream_index < length(files)
            bl = barcodelengths[stream_index]
            barcode = barcode[1:bl] ## chop random barcode at specified length
        end

        ids = split(read_id)

        barcode_text = ifelse(isempty(barcode), "", string(":", barcode))

        println(streams[stream_index], ids[1], barcode_text, " ", ids[2], ":", index)
        println(streams[stream_index], read_read)
        println(streams[stream_index], read_qid)
        println(streams[stream_index], read_qs)

        counts[stream_index] += 1
        total_reads += 1

    end

    close(rio)
    close(iio)
    close.(streams)

    println("########################### Summary ###########################")
    println("File\tCount\tProportion")
    s = sum(counts)
    for (f, c) in zip(files, counts)
        println(f, "\t", c, "\t", c/s)
    end
    println("###############################################################")

    println("[SFQ]\tTotal Reads    :\t", total_reads)
    println("[SFQ]\tComplete in    :\t", time() - st, " seconds.")
end


text:
using Base.Test
using EGR

println("TestM1P1")

a = MinusPlusOneVector([1,0,-1],Set([-1]))
b = MinusPlusOneVector([-1,-1,1])
@test a.field == b.field
@test a.numPlus == b.numPlus
@test_throws ArgumentError MinusPlusOneVector([1,1,-2])
@test_throws ArgumentError MinusPlusOneVector([1.0,1,-2])

println("TestM1P1 successful!")
println()


text:
using OnlinePackage
using Documenter: deploydocs, makedocs

makedocs(
    sitename = "OnlinePackage.jl",
    modules = [OnlinePackage],
    doctest = false
)
deploydocs(repo = "github.com/bramtayl/OnlinePackage.jl.git")


text:
#iterate through the momentums trategies and put them into the dataframe
function momentumstrategies!(panel::DataFrame;
  strategies::Dict=PARAM[:momentumstrategies], Fconditional::NSymbol=nothing)
  for (Fmomentum, strategy) ∈ strategies
    #we will attatch this message to an error if the property names don't line up as expected
    if strategy.type ≡ :longshort
      longshortmomentumstrategy!(panel, strategy; Fmomentum, Fconditional)
    elseif strategy.type ≡ :longonly ||strategy.type ≡ :long
      longonlymomentumstrategy!(panel, strategy; Fmomentum, Fconditional)
    else
      throw("unrecognized momentum strategy type")
    end
  end

end

#this function computes a long-short momentum value
function longshortmomentumstrategy!(panel::DataFrame, strategy;
  Fmomentum=throw("Fmomentum is required"),
  Fconditional::NSymbol,
  )
  local longperiod::DatePeriod
  local shortperiod::DatePeriod

  datemin=minimum(skipmissing(panel.date))
  datemax=maximum(skipmissing(panel.date))
  numdates::Int=panel.date |> unique |> length

  #the procedure for computing momentum can either use a tolerant return-based approach
  #or a strict index based approach
  if :Flret ∈ keys(strategy)
    @assert :Findex ∉ keys(strategy)
    @unpack shortperiod, longperiod, shortfraction, longfraction, Flret = strategy

    Fmeasurelong = Symbol(Flret,replace("$longperiod", " "=>""))
    Fmeasureshort = Symbol(Flret,replace("$shortperiod", " "=>""))

    trailingret!(panel, [longperiod, shortperiod],
      [longfraction, shortfraction],
      Fmeasures=[Fmeasurelong, Fmeasureshort]; Flret, Fconditional)
  elseif :Findex ∈ keys(strategy)
    @assert :Flret ∉ keys(strategy)
    @unpack shortperiod, longperiod, Findex = strategy
      @assert longperiod > shortperiod

    Fmeasurelong = Symbol(Findex,replace("$longperiod", " "=>""))
    Fmeasureshort = Symbol(Findex,replace("$shortperiod", " "=>""))

    stricttrailingret!(panel, [longperiod, shortperiod],
      Fmeasures=[Fmeasurelong, Fmeasureshort]; Findex, Fconditional)
  end

  panel[!, Fmomentum] = panel[!, Fmeasurelong] .- panel[!, Fmeasureshort]

  if PARAM[:testmomentum]
    testpermno::Int = PARAM[:testpermnomult]
    testview::SubDataFrame = view(panel, (panel.permno .% testpermno) .== 0, :)
    if size(testview,1)>20_000 #prevents scenarios where the file is larger than needed
      testview=view(testview, 1:20_000,:)
    end
    testview |> CSV.write("$(PARAM[:testpath])\\$(PARAM[:crspdatalabel])_$(Fconditional)_$Fmomentum.csv")
  end

  select!(panel, Not([Fmeasurelong,Fmeasureshort]))

  return nothing
end

#this function computes a long only momentum value
function longonlymomentumstrategy!(panel::DataFrame, strategy;
  Fmomentum=throw("Fmomentum is required"),
  Fconditional::NSymbol,
  )

  #@info "$Fmomentum: minpoints=$minpoints"

  if :Flret ∈ keys(strategy)
    @assert :Findex ∉ keys(strategy)
    @unpack Flret, period, fraction = strategy
    Fmeasure = Symbol(Flret,replace("$period", " "=>""))

    trailingret!(panel, [calendarinterval], [fraction], Fmeasures=[Fmomentum]; Flret, Fconditional)
  elseif :Findex ∈ keys(strategy)
    @assert :Flret ∉ keys(strategy)
    @unpack Findex, period = strategy
    Fmeasure = Symbol(Findex,replace("$period", " "=>""))

   stricttrailingret!(panel, [period], Fmeasures=[Fmomentum]; Findex, Fconditional)
  else
    throw("return or index not found in momentum strategy $strategy")
  end

  if PARAM[:testmomentum]
    testpermno::Int = PARAM[:testpermnomult]
    testview::SubDataFrame = view(panel, (panel.permno .% testpermno) .== 0, :)
    if size(testview,1)>20_000 #prevents scenarios where the file is larger than needed
      testview=view(testview, 1:20_000,:)
    end
    testview |> CSV.write("$(PARAM[:testpath])\\$(PARAM[:crspdatalabel])_$(Fconditional)_$Fmomentum.csv")
  end

  return nothing
end

#############################trailingret########################
#used to compute the actual trailing returns
hpr(v::AbstractVector{T}) where T = hpr(exp(sum(v))-one(T))
hpr(x::Real) = x
hpr(::Missing) = throw("missings found when computing trailing ret!!!")
function trailingret!(crsp::DataFrame, calendarintervals::Vector{<:DatePeriod},
  fractions::Vector{Float64};
  Flret::Symbol,
  Fconditional::NSymbol,
  Fmeasures::Vector{Symbol})

  for (calendarinterval, fraction, Fmeasure) ∈ zip(calendarintervals, fractions, Fmeasures)
    crsp[!, Fmeasure] = Vector{MFloat64}(undef, size(crsp,1))
    trailingmeasure!(hpr, crsp, calendarinterval, fraction; Ftarget=Flret, Fconditional, Fmeasure)
  end

  return nothing
end

idxret(v::AbstractVector{T}) where T = idxret(v[end]/v[1] - one(T))
idxret(x::Real) = isfinite(x) ? x : idxret(missing)
idxret(::Missing) = throw("missings found when computing trailing index ret!!!")
function stricttrailingret!(crsp::DataFrame, calendarintervals::Vector{<:DatePeriod};
  Findex::Symbol,
  Fconditional::NSymbol,
  Fmeasures::Vector{Symbol})

  for (calendarinterval, Fmeasure) ∈ zip(calendarintervals, Fmeasures)
    crsp[!, Fmeasure] = Vector{MFloat64}(undef, size(crsp,1))
    stricttrailingmeasure!(idxret, crsp, calendarinterval,; Ftarget=Findex, Fconditional, Fmeasure)
  end

  return nothing
end

############post-processing of momentum strategies
function conditionmomentumstrategies!(panel::DataFrame; strategies::Dict=PARAM[:momentumstrategies])
  for (Fmomentum, strategy) ∈ strategies
    #winsorize the momentum
    FWmomentum::Symbol = Symbol(:W, Fmomentum)
    FWXmomentum::Symbol = Symbol(:WX, Fmomentum)
    panel[!,FWmomentum] = winsorizequantile(panel[!, Fmomentum], PARAM[:qwinsorize], twosided=true)
    winsorizewithin!(panel, Fsource=Fmomentum, prop=PARAM[:qwinsorize], Fgroup=:date,
      twosided=true, Fnew=FWXmomentum)
    lagwithin2!(panel, [Fmomentum, FWmomentum, FWXmomentum], :permno, date=:date) #this could be parallelized

    #these wont work in the current setup
    #FWiivolWmomentum = Symbol(:Wiivol, FWmomentum)
    #panel[!,Symbol(:Wiivol, FWmomentum)] =  panel[!, :Wiivol] .* panel[!, FWmomentum]
    #panel[!,Symbol(:LWiivolL, FWmomentum)] =  panel[!, :LWiivol] .* panel[!, Symbol(:L, FWmomentum)]
    #lagwithin2!(panel, [FWiivolWmomentum], :permno, date=:date) #this could be parallelized
  end

  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
mutable struct viewer_link_data_t <: LCMType
    name::String
    robot_num::Int32
    num_geom::Int32
    geom::Vector{viewer_geometry_data_t}
end

@lcmtypesetup(viewer_link_data_t,
    geom => (num_geom,)
)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

abstract type AbstractMCMCStats end
AbstractMCMCStats

Base.convert(::Type{AbstractMCMCCallback}, x::AbstractMCMCStats) = MCMCAppendCallback(x)

MCMCAppendCallback(x::AbstractMCMCStats, nonzero_weights::Bool = true) =
    MCMCAppendCallback(x, 1, get_samples!, nonzero_weights)

struct MCMCNullStats <: AbstractMCMCStats end

Base.push!(stats::MCMCNullStats, sv::DensitySampleVector) = stats

Base.append!(stats::MCMCNullStats, sv::DensitySampleVector) = stats

struct MCMCBasicStats{L<:Real,P<:Real} <: AbstractMCMCStats
    param_stats::BasicMvStatistics{P,FrequencyWeights}
    logtf_stats::BasicUvStatistics{L,FrequencyWeights}
    mode::Vector{P}

    function MCMCBasicStats{L,P}(m::Integer) where {L<:Real,P<:Real}
        param_stats = BasicMvStatistics{P,FrequencyWeights}(m)
        logtf_stats = BasicUvStatistics{L,FrequencyWeights}()
        mode = fill(oob(P), m)

        new{L,P}(
            param_stats,
            logtf_stats,
            mode
        )
    end
end

function MCMCBasicStats(::Type{S}, nparams::Integer) where {P,T,W,S<:DensitySample{P,T,W}}
    SL = promote_type(T, Float64)
    SP = promote_type(P, W, Float64)
    MCMCBasicStats{SL,SP}(nparams)
end

MCMCBasicStats(chain::MCMCIterator) = MCMCBasicStats(sample_type(chain), nparams(chain))

function MCMCBasicStats(sv::DensitySampleVector)
    stats = MCMCBasicStats(eltype(sv), innersize(sv.params, 1))
    append!(stats, sv)
end

function Base.empty!(stats::MCMCBasicStats)
    empty!(stats.param_stats)
    empty!(stats.logtf_stats)
    fill!(stats.mode, oob(eltype(stats.mode)))

    stats
end

function Base.push!(stats::MCMCBasicStats, s::DensitySample)
    push!(stats.param_stats, s.params, s.weight)
    if s.logdensity > stats.logtf_stats.maximum
        stats.mode .= s.params
    end
    push!(stats.logtf_stats, s.logdensity, s.weight)
    stats
end

function Base.append!(stats::MCMCBasicStats, sv::DensitySampleVector)
    for i in eachindex(sv)
        p = sv.params[i]
        w = sv.weight[i]
        l = sv.logdensity[i]
        push!(stats.param_stats, p, w)  # Memory allocation (view)!
        if sv.logdensity[i] > stats.logtf_stats.maximum
            stats.mode .= p  # Memory allocation (view)!
        end
        push!(stats.logtf_stats, l, w)
        stats
    end
    stats
end

nparams(stats::MCMCBasicStats) = stats.param_stats.m

nsamples(stats::MCMCBasicStats) = stats.param_stats.cov.sum_w

function Base.merge!(target::MCMCBasicStats, others::MCMCBasicStats...)
    for x in others
        if (x.logtf_stats.maximum > target.logtf_stats.maximum)
            target.mode .= x.mode
        end
        merge!(target.param_stats, x.param_stats)
        merge!(target.logtf_stats, x.logtf_stats)
    end
    target
end

Base.merge(a::MCMCBasicStats, bs::MCMCBasicStats...) = merge!(deepcopy(a), bs...)

function bat_stats(mcmc_stats::MCMCBasicStats)
    (
        mode = mcmc_stats.mode,
        mean = mcmc_stats.param_stats.mean,
        cov = mcmc_stats.param_stats.cov
    )
end


text:
# Helper functions for parallel external sort

"""
Get new chunck of sorted data.
"""
function parallel_get_new_chunk(input_file_name,
                                input_file_position,
                                size_of_sorted_chunk)
    input_file = open(input_file_name, "r")
    input_buffer_array = Array(Int32, 0)

    # input_file_position < 0 means that the file was read already.
    if input_file_position < 0
         return input_buffer_array, -1
    end

    seek(input_file, input_file_position)
    for j in 1:size_of_sorted_chunk
        if eof(input_file)
            break
        else
            push!(input_buffer_array, read(input_file, Int32))
        end
    end

    if eof(input_file)
        new_input_position = -1
    else
        new_input_position = position(input_file)
    end

    return input_buffer_array, new_input_position
end

max_stars_repo_name: rgaiacs/external-sort

max_issues_repo_name: rgaiacs/external-sort

max_forks_repo_name: rgaiacs/external-sort


text:
#TODO : animate bands

################################################################################
### WRAPPERS
################################################################################
"""
    bands, gaps = band_structure(sim, num_bloch::Int; num_bands=5, zone=:half, interpolation=:cubic)

band structure along path specified by `zone`, which must be one of :reduced, :half, :full,
or an integer or array of integers.

Zone specification is

| zone | number |
|:----:|:----:|
| 3    |  4   |
| 1    |  2   |

Except for reduced, each zone is defined by (_ ↑ ←, ↗ ↓ ↖)
"""
function band_structure(sim::Simulation, num_bloch::Int; num_bands=5, zone=:half,
        parallel=nprocs()>1, interpolation=:cubic, calc_type="band structure ",
        disp_opt=true)

    k_paths, zones, order = band_paths(sim, num_bloch, zone)
    bands = Array{Array{Array{Array{Float64,2},1},1},1}(undef, length(zones))
    if disp_opt
        pg = Progress(length(zones)*length(zones[1])*3*num_bloch, PROGRESS_UPDATE_TIME::Float64, calc_type)
    else
        pg = nothing
    end
    for i ∈ eachindex(bands)
        bands[i] = Array{Array{Array{Float64,2},1},1}(undef, length(zones[i]))
        for j ∈ eachindex(zones[i])
            bands[i][j] = Array{Array{Float64,2},1}(undef, 3)
            for m ∈ 1:3
                kb = (k_paths[zones[i][j][m]][1][order[i][j][m]], k_paths[zones[i][j][m]][2][order[i][j][m]])
                if parallel
                    bands[i][j][m] = real(bands_only_p(sim, kb, num_bands, pg))
                else
                    bands[i][j][m] = real(bands_only(sim, kb, num_bands, pg))
                end
            end
        end
    end
    gaps, bands_itp = find_gaps(bands, interpolation)

    return bands, gaps
end

"""
    bands, gaps, kxy = band_structure(sim, num_bloch::Array; num_bands=5, zone=:half)

band structure as full surface over zone specified by `zone`, which must be one of
:reduced, :half, :full.
"""
function band_structure(sim::Simulation, num_bloch::Array{Int,1}; num_bands=5, zone=:half,
    parallel=nprocs()>1, calc_type="band structure ", disp_opt=true, interpolation=nothing,
    line=nothing)

    if length(num_bloch)==1
        num_bloch = fill(num_bloch[1],2)
    end

    if zone == :reduced
        kx = (0:1/num_bloch[1]:1)*2π/2sim.lat.a
        ky = (0:1/num_bloch[2]:1)*2π*sqrt(sim.lat.sin²θ)/2sim.lat.b
    elseif zone == :half
        kx = (0:1/num_bloch[1]:1)*2π/sim.lat.a
        ky = (0:1/num_bloch[2]:1)*2π*sqrt(sim.lat.sin²θ)/2sim.lat.b
    elseif zone == :full
        kx = (0:1/num_bloch[1]:1)*2π/sim.lat.a
        ky = (0:1/num_bloch[2]:1)*2π*sqrt(sim.lat.sin²θ)/sim.lat.b
    else
        throw(ArgumentError("invalid zone $zone, must be one of :reduced, :half, :full"))
    end
    KAB = bravais_coordinates.(kx,ky',Ref(Bravais(sim.lat; :a=>2π/sim.lat.a, :b=>2π/sim.lat.b)))
    ka = Array{Float64}(undef, size(KAB))
    kb = Array{Float64}(undef, size(KAB))
    for i ∈ eachindex(KAB)
        ka[i] = KAB[i][1]
        kb[i] = KAB[i][2]
    end

    if disp_opt
        pg = Progress(prod(size(ka)), PROGRESS_UPDATE_TIME, calc_type)
    else
        pg = nothing
    end
    if parallel
        bands = bands_only_p(sim, (ka, kb), num_bands, pg)
    else
        bands = bands_only(sim, (ka, kb), num_bands, pg)
    end
    gaps, bands_itp = find_gaps([[[real(bands)]]], interpolation)

    return real(bands), gaps, [kx, ky]
end

function band_structure(sim::Simulation, k_bloch::Tuple{AbstractArray{S,M},AbstractArray{T,M}};
    num_bands=5, parallel=nprocs()>1, interpolation=:cubic,
    calc_type = "band structure ",
    pg::Progress=Progress(prod(size(k_bloch[1])), PROGRESS_UPDATE_TIME, calc_type),
    disp_opt=true) where S where T where M

    if !disp_opt
        pg = nothing
    end
    if parallel
        bands = real(bands_only_p(sim, k_bloch, num_bands, pg))
    else
        bands = real(bands_only(sim, k_bloch, num_bands, pg))
    end
    gaps, bands_itp = find_gaps([[[bands]]], interpolation)

    return bands, gaps, bands_itp
end

function band_structure(sim::Simulation, k_bloch::Tuple{AbstractArray{S,M},AbstractArray{T,M}},
    num_bloch_interpolations::Int, reciprocal_basis=true;
    num_bands=5, parallel=nprocs()>1, interpolation=:cubic,
    calc_type = "band structure ",
    pg::Progress=Progress(num_bloch_interpolations*(length(k_bloch[1])-1), PROGRESS_UPDATE_TIME, calc_type),
    disp_opt=true) where S where T where M

    kas = Array{Float64}(undef, num_bloch_interpolations*(length(k_bloch[1])-1))
    kbs = deepcopy(kas)
    for i ∈ 1:length(k_bloch[1])-1
        if reciprocal_basis
            ka_start = k_bloch[1][i]*2π/sim.lat.a
            kb_start = k_bloch[2][i]*2π/sim.lat.b

            ka_stop = k_bloch[1][i+1]*2π/sim.lat.a
            kb_stop = k_bloch[2][i+1]*2π/sim.lat.b
        else
            ka_start = k_bloch[1][i]
            kb_start = k_bloch[2][i]

            ka_stop = k_bloch[1][i+1]
            kb_stop = k_bloch[2][i+1]
        end
        ka_temp = LinRange(ka_start, ka_stop, num_bloch_interpolations)
        kb_temp = LinRange(kb_start, kb_stop, num_bloch_interpolations)
        for j ∈ 1:num_bloch_interpolations
            kas[(i-1)*num_bloch_interpolations + j] = ka_temp[j]
            kbs[(i-1)*num_bloch_interpolations + j] = kb_temp[j]
        end
    end

    bands, gaps, bands_itp = band_structure(sim, (kas, kbs); num_bands=num_bands, parallel=parallel,
        interpolation=interpolation, calc_type=calc_type, pg=pg, disp_opt=disp_opt)

    return bands, gaps, bands_itp, (kas, kbs)
end

"""
    bands, gaps, ks = band_structure(sim; waveguide, num_bloch=17, num_bands=5, interpolation=:cubic)

bands and gaps of photonic crystal that defines background of `waveguide`
"""
function band_structure(sim::Simulation; waveguide, num_bloch=17, num_bands=5,
    parallel=nprocs()>1, interpolation=:cubic, zone=:half, disp_opt=true)

    wg_sim = extract_waveguide_simulation(sim, waveguide)
    pc_sim = Simulation(2, wg_sim)
    bands, gaps = band_structure(pc_sim, num_bloch; num_bands=num_bands,
        parallel=parallel, interpolation=interpolation, zone=zone,
        calc_type="waveguide band structure ", disp_opt=disp_opt)
    return bands, gaps
end

################################################################################
### CORE
################################################################################
"""
    bands = bands_only(sim, k_bloch::Tuple, num_bands, progress)
"""
function bands_only(
    sim::Simulation,
    k_bloch::Tuple{AbstractArray{S,M},AbstractArray{T,M}},
    num_bands,
    pg=nothing) where S where T where M

    ka_bloch = k_bloch[1]
    kb_bloch = k_bloch[2]

    bands = Array{ComplexF64, M+1}(undef, num_bands, size(ka_bloch)...)
    for i ∈ CartesianIndices(ka_bloch)
        bands[:,i] = eig_k(sim, 0.001, num_bands; ka=ka_bloch[i], kb=kb_bloch[i])[1]

        if typeof(pg)==Progress
            next!(pg)
        end
    end

    return bands
end

"""
    gaps, bands_itp = find_gaps(bands, inerpolation)
"""
function find_gaps(bands::Array{Array{Array{Array{Float64,M},1},1},1}, interpolation) where M

    num_bands = size(bands[1][1][1],1)
    bands_itp = Array{AbstractExtrapolation}(undef, num_bands)

    if M==2 # not surface
        band_tops = Array{Float64}(undef, num_bands)
        band_bottoms = Array{Float64}(undef, num_bands)
        band_tops = zeros(Float64, num_bands)
        band_bottoms = Inf*ones(Float64, num_bands)
        for i ∈ eachindex(bands)
            for j ∈ eachindex(bands[i])
                for m ∈ eachindex(bands[i][j])
                    tops = findmax(real(bands[i][j][m]); dims=2)[2]
                    bottoms = findmin(real(bands[i][j][m]); dims=2)[2]
                    for n ∈ 1:num_bands
                        if interpolation == :linear
                            interpolator = BSpline(Linear())
                        elseif interpolation == :quadratic
                            interpolator = BSpline(Quadratic(Line(OnGrid())))
                        elseif interpolation == :cubic
                            interpolator = BSpline(Cubic(Line(OnGrid())))
                        else
                            throw(ArgumentError("invalid order $(interpolation)"))
                        end
                        bands_itp[n] = extrapolate(interpolate(bands[i][j][m][n,:], interpolator), Reflect())
                        band_tops[n] = max(band_tops[n], bands_itp[n](optimize(x->-bands_itp[n](x[1]), [float(tops[n][2])], BFGS()).minimizer)[1])::Float64
                        band_bottoms[n] = min(band_bottoms[n], bands_itp[n](optimize(x->abs(bands_itp[n](x[1])), [float(bottoms[n][2])], BFGS()).minimizer)[1])::Float64
                    end
                end
            end
        end
        gaps = band_bottoms[2:end]-band_tops[1:end-1]
        gaps_start = band_tops[1:end-1][gaps.>0]
        gaps_stop = band_tops[1:end-1][gaps.>0] + gaps[gaps.>0]
    else #surface
        bands = bands[1][1][1]
        tops = (findmax(bands; dims=(2,3))[1])
        bottoms = findmin(bands; dims=(2,3))[1]
        gaps = bottoms[2:end]-tops[1:end-1]
        gaps_start = tops[1:end-1][gaps.>0]
        gaps_stop = tops[1:end-1][gaps.>0] + gaps[gaps.>0]
    end

    return (gaps_start, gaps_stop), bands_itp
end

################################################################################
### AUXILLIARIES
################################################################################
"""
    k_paths, zones, order = band_paths(sim, num_bloch, zone)

`zone` can be an integer, array, or one of `:reduced`, `:half`, `:full`
"""
function band_paths(sim::Simulation, num_bloch, zone=:half)

    k1 = [0,0]
    k2 = 2π*(sim.lat.v1/2sim.lat.a)
    k3 = 2π*(sim.lat.v1/1sim.lat.a)

    k4 = 2π*(sim.lat.v2/2sim.lat.b)
    k5 = 2π*(sim.lat.v1/2sim.lat.a+sim.lat.v2/2sim.lat.b)
    k6 = 2π*(sim.lat.v1/1sim.lat.a+sim.lat.v2/2sim.lat.b)

    k7 = 2π*(sim.lat.v2/1sim.lat.b)
    k8 = 2π*(sim.lat.v1/2sim.lat.a+sim.lat.v2/1sim.lat.b)
    k9 = 2π*(sim.lat.v1/1sim.lat.a+sim.lat.v2/1sim.lat.b)

    ka_bloch01 = LinRange(k1[1], k2[1], num_bloch)
    ka_bloch02 = LinRange(k2[1], k3[1], num_bloch)

    ka_bloch03 = LinRange(k1[1], k4[1], num_bloch)
    ka_bloch04 = LinRange(k2[1], k5[1], num_bloch)

    ka_bloch05 = LinRange(k1[1], k5[1], num_bloch)
    ka_bloch06 = LinRange(k2[1], k4[1], num_bloch)
    ka_bloch07 = LinRange(k2[1], k6[1], num_bloch)
    ka_bloch08 = LinRange(k3[1], k5[1], num_bloch)

    ka_bloch09 = LinRange(k4[1], k5[1], num_bloch)
    ka_bloch10 = LinRange(k5[1], k6[1], num_bloch)

    ka_bloch11 = LinRange(k4[1], k7[1], num_bloch)
    ka_bloch12 = LinRange(k5[1], k8[1], num_bloch)

    ka_bloch13 = LinRange(k4[1], k8[1], num_bloch)
    ka_bloch14 = LinRange(k5[1], k7[1], num_bloch)
    ka_bloch15 = LinRange(k5[1], k9[1], num_bloch)
    ka_bloch16 = LinRange(k6[1], k8[1], num_bloch)

    kb_bloch01 = LinRange(k1[2], k2[2], num_bloch)
    kb_bloch02 = LinRange(k2[2], k3[2], num_bloch)

    kb_bloch03 = LinRange(k1[2], k4[2], num_bloch)
    kb_bloch04 = LinRange(k2[2], k5[2], num_bloch)

    kb_bloch05 = LinRange(k1[2], k5[2], num_bloch)
    kb_bloch06 = LinRange(k2[2], k4[2], num_bloch)
    kb_bloch07 = LinRange(k2[2], k6[2], num_bloch)
    kb_bloch08 = LinRange(k3[2], k5[2], num_bloch)

    kb_bloch09 = LinRange(k4[2], k5[2], num_bloch)
    kb_bloch10 = LinRange(k5[2], k6[2], num_bloch)

    kb_bloch11 = LinRange(k4[2], k7[2], num_bloch)
    kb_bloch12 = LinRange(k5[2], k8[2], num_bloch)

    kb_bloch13 = LinRange(k4[2], k8[2], num_bloch)
    kb_bloch14 = LinRange(k5[2], k7[2], num_bloch)
    kb_bloch15 = LinRange(k5[2], k9[2], num_bloch)
    kb_bloch16 = LinRange(k6[2], k8[2], num_bloch)

    k_bloch01 = (ka_bloch01, kb_bloch01)
    k_bloch02 = (ka_bloch02, kb_bloch02)
    k_bloch03 = (ka_bloch03, kb_bloch03)
    k_bloch04 = (ka_bloch04, kb_bloch04)
    k_bloch05 = (ka_bloch05, kb_bloch05)
    k_bloch06 = (ka_bloch06, kb_bloch06)
    k_bloch07 = (ka_bloch07, kb_bloch07)
    k_bloch08 = (ka_bloch08, kb_bloch08)
    k_bloch09 = (ka_bloch09, kb_bloch09)
    k_bloch10 = (ka_bloch10, kb_bloch10)
    k_bloch11 = (ka_bloch11, kb_bloch11)
    k_bloch12 = (ka_bloch12, kb_bloch12)
    k_bloch13 = (ka_bloch13, kb_bloch13)
    k_bloch14 = (ka_bloch14, kb_bloch14)
    k_bloch15 = (ka_bloch15, kb_bloch15)
    k_bloch16 = (ka_bloch16, kb_bloch16)

    k_paths = (k_bloch01, k_bloch02, k_bloch03, k_bloch04,
        k_bloch05, k_bloch06, k_bloch07, k_bloch08,
        k_bloch09, k_bloch10, k_bloch11, k_bloch12,
        k_bloch13, k_bloch14, k_bloch15, k_bloch16)

    zones = (((1,4,9),(5,3,6)),
             ((2,4,10),(8,3,7)),
             ((1,12,9),(14,11,13)),
             ((2,12,10),(15,11,16)))

    polarity = (((1,1,-1),(-1,1,-1)),
                ((-1,1,1),(-1,1,-1)),
                ((1,-1,-1),(1,-1,1)),
                ((-1,-1,1),(1,-1,1)))
    if zone==:reduced
        zones = (((1,4,5),),)
        polarity = (((1,1,-1),),)
    elseif zone==:half
        zones = (((1,4,9),(5,3,6)),
                 ((2,4,10),(8,3,7)))
        polarity = (((1,1,-1),(-1,1,-1)),
                    ((-1,1,1),(-1,1,-1)))
    elseif zone==:full
        nothing
    elseif typeof(zone) <: Int
        zones = zones[[zone]]
        polarity = polarity[[zone]]
    elseif typeof(zone) <: Array{Int}
        zones = zones[zone]
        polarity = polarity[zone]
    else
        throw(ArgumentError("invalid zone specification $zone, must be one of :reduced, :half, :full, or an integer or array of integers."))
    end

    order = Array{Array{Array{AbstractArray}}}(undef,length(zones))
    for i ∈ 1:length(order)
        order[i] = Array{Array{AbstractArray}}(undef, length(zones[i]))
        for j ∈ 1:length(zones[i])
            order[i][j] = Array{AbstractArray}(undef, 3)
            for m ∈ 1:3
                if polarity[i][j][m]>0
                    step = 1
                    start = 1
                    stop = num_bloch
                else
                    step = -1
                    start = num_bloch
                    stop = 1
                end
                order[i][j][m] = start:step:stop
            end
        end
    end

    return k_paths, zones, order
end


text:
# query.jl

"""
    @query(expr)::QueryResult

Execute a query with expression
"""
macro query(sym::Symbol)
    query_request(Expr(:block, sym))
end

macro query(expr::Expr)
    query_request(expr)
end

function query_request(expr::Expr)
    if expr.head in [:(=), :kw]
        lhs,rhs = expr.args
        lquote = destchains(lhs)
        rquote = destchains(rhs)
        lhs_issym = isa(lhs, Symbol)
        rhs_issym = isa(rhs, Symbol)
        lsymbol = (lhs,)
        rsymbol = (rhs,)
        quote
            (ldest,lsym) = $lquote
            (rdest,rsym) = $rquote
            endpoint = endpointof(rdest, ldest)
            rvalue = $lhs_issym ? first(rsym) : rsym
            if $rhs_issym
                rvalue = valuate(first($rsymbol))
            end
            if isa(ldest, App) && isa(rdest, App) && (ldest.url != rdest.url)
                if $rhs_issym
                else
                    param = Dict("type"=>"Getter", "lhs"=>wrap_symbol(rsym))
                    verb = "/query"
                    dict = request(endpoint, verb, param)
                    result = QueryResult(dict, endpoint, verb, param)
                    rvalue = result.info.value
                end
            end
            if $lhs_issym
                (firstsym,) = $lsymbol
                Swifter.env[firstsym] = rvalue
            else
                endpoint = endpointof(ldest, rdest)
                param = Dict("type"=>"Setter", "lhs"=>wrap_symbol(lsym), "rhs"=>wrap_symbol(rvalue))
                verb = "/query"
                dict = request(endpoint, verb, param)
                QueryResult(dict, endpoint, verb, param)
            end
        end
    else
        lquote = destchains(expr)
        (sym,) = expr.args
        quot = quote
            (ldest,lsym) = $lquote
            endpoint = endpointof(ldest,nothing)
            param = Dict("type"=>"Getter", "lhs"=>wrap_symbol(lsym))
            verb = "/query"
            dict = request(endpoint, verb, param)
            QueryResult(dict, endpoint, verb, param)
        end
        if 1 == length(expr.args) && isa(sym, Symbol)
            quote
                (ldest,lsym) = $lquote
                (val,) = lsym
                val == sym ? $quot : val
            end
        else
            quot
        end
    end
end

function endpointof(a::Union{Void,App}, b::Union{Void,App})
    if isa(a, App)
        return a
    elseif isa(b, App)
        return b
    else
        return current_endpoint
    end
end

chains(ex::Union{Symbol,Int,Float64,Bool,AbstractString}) = chains(Expr(:block, ex))
chains(expr::Expr) = deserial(chaining(expr, 0, []))

function valuate(sym::Symbol)
    if haskey(Swifter.env, sym)
        Swifter.env[sym]
    elseif isdefined(Main, sym)
        getfield(Main, sym)
    else
        sym
    end
end

destchains(ex::Union{Symbol,Int,Float64,Bool,AbstractString}) = destchains(Expr(:block, ex))

function destchains(expr::Expr)
    symbols = chains(expr)
    firstsym = first(symbols)
    quote
        sym = $symbols[1]
        try
            sym = $(esc(firstsym))
        end
        if isa(sym, Symbol)
            $symbols[1] = isa(sym, Symbol) ? valuate(sym) : sym
        end
        if isa(sym, QueryResult)
            (sym.app, vcat(sym, $symbols[2:end]))
        else
            (nothing, $symbols)
        end
    end
end

function chaining(expr::Expr, depth::Int, symbols::Vector)
    if :call == expr.head
        firstarg = first(expr.args)
        callarg = CallArg("", Any[])
        if isa(firstarg, Symbol)
            callarg.name = firstarg
        else
            for ex in firstarg.args
                if isa(ex, Expr)
                    for el in ex.args[1:end]
                        if isa(el, Expr)
                            chaining(el, depth+1, symbols)
                        elseif isa(el, QuoteNode)
                            push!(symbols, el.value)
                        else
                            if length(ex.args) > 1
                                push!(symbols, el)
                            else
                                callarg.name = el
                            end
                        end
                    end
                elseif isa(ex, QuoteNode)
                    callarg.name = ex.value
                else
                    push!(symbols, ex)
                end
            end
        end
        args = Any[]
        if length(expr.args) > 1
            for ex in expr.args[2:end]
                if isa(ex, Expr)
                    push!(args, chaining(ex, depth+1, Any[]))
                else
                    push!(args, ex)
                end
            end
        end
        callarg.args = args
        push!(symbols, callarg)
    else
        for ex in expr.args
            if isa(ex, Expr)
                chaining(ex, depth+1, symbols)
            elseif isa(ex, QuoteNode)
                push!(symbols, ex.value)
            else
                push!(symbols, ex)
            end
        end
    end
    symbols
end

function deserial(symbols::Vector)
    syms = Any[]
    for sym in symbols
        if isa(sym, CallArg)
            if isempty(sym.args)
                push!(syms, (:call, sym.name))
            else
                push!(syms, (:call, (sym.name, wrap_symbol(sym.args))))
            end
        else
            push!(syms, sym)
        end
    end
    syms
end

wrap_symbol(ex::Union{QueryResult,Symbol,Int,Float64,Bool,AbstractString}) = wrap_symbol(Any[ex])

function wrap_symbol(lhs::Vector)
    vals = Any[]
    for ex in lhs
        if isa(ex, Symbol)
            push!(vals, (:symbol, ex))
        elseif isa(ex, QueryResult)
            if isa(ex.info.address, Void)
                push!(vals, ex.info.value)
            else
                push!(vals, (:address, string(ex.info.address)))
            end
        else
            push!(vals, ex)
        end
    end
    vals
end


text:
# For running tests in a "clean" environment.
module Sandbox end

function doctest(io::IO, document::Documentation, verbose::Bool)
    # Import the module being tested into the sandbox.
    eval(Sandbox, Expr(:toplevel, Expr(:using, symbol("$(document.modname)"))))

    print_with_color(:blue, io, "[DOCTEST] running on module $(document.modname).\n")
    for (obj, entry) in document.entries
        println(io, " • [ENTRY] $(obj)")
        for block in entry.docs.content
            if isa(block, Markdown.Code) # block.language == :julia
                print_with_color(:purple, io, "  ∘ [EVAL]\n")
                success = true
                try
                    if endswith(block.code, "\n")
                        print_with_color(:blue, io, "  ~ [SKIPPED]\n")
                        continue
                    end
                    eval(Sandbox, parse("let\n$(block.code)\nend"))
                catch err
                    success = false
                    print_with_color(:red, io, "  - [FAILURE]\n")
                    println(io, "\n", block.code, "\n")
                    print_with_color(:red, io, string(err), "\n\n")
                end
                if success
                    print_with_color(:green, io, "  + [SUCCESS]\n")
                    verbose && print_with_color(:white, io, "\n", block.code, "\n\n")
                end
            end
        end
    end
end

doctest(io::IO, m::Module, verbose) =
    isdefined(m, METADATA) ? doctest(io, getfield(m, METADATA), verbose) :
    println(io, "No docs found in module $(m).")

@doc """
Test all code blocks in the docstrings in a module `m`.

Run `doctest(Docile)` to see some sample output produced by this
function.
""" {
    :tags => ["testing", "code blocks"],
    :parameters => [
        (:m, "the module to test"),
        (:verbose, "optional switch to view detailed output of each test")
        ]
    } ->
function doctest(m::Module; verbose = false)
    doctest(STDOUT, m, verbose)
end