text:
import MD5
import GZip
import Tar
import Arrow
import DataFrames

# URL for the taxonomy dump
const ncbi_ftp = "https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/new_taxdump/"
const archive = ncbi_ftp * "new_taxdump.tar.gz"
const checksum = archive * ".md5"

if !haskey(ENV, "NCBITAXONOMY_PATH")
    @warn """
    The environmental varialbe NCBITAXONOMY_PATH is not set, so the tables will
    be stored in the package path. This is not ideal, and you really should set
    the NCBITAXONOMY_PATH.
    """
end
const taxpath = get(ENV, "NCBITAXONOMY_PATH", joinpath(@__DIR__, "..", "deps"))
ispath(taxpath) || mkpath(taxpath)

chk_file = download(checksum)
chk = split(readlines(chk_file)[1], " ")[1]
@info "Checksum of the most recent NCBI taxonomy: $(chk)"

function download_dump(url, chk, dest)
    @info "Downloading the taxonomy data from $(url)"
    if ispath(joinpath(taxpath, dest))
        @info "Removing the previous version of the taxonomy"
        rm(joinpath(taxpath, dest); force=true, recursive=true)
        mkpath(joinpath(taxpath, dest))
    else
        mkpath(joinpath(taxpath, dest))
    end
    arc = download(url)
    vrf = bytes2hex(open(MD5.md5, arc))
    vrf == chk || throw(ErrorException("Wrong checksum for the NCBI taxonomy archive file - unable to download"))
    write(joinpath(taxpath, ".checksum"), vrf)
    Tar.extract(GZip.open(arc), joinpath(taxpath, dest))
end

# The next block is about making sure that we don't download something that has
# not changed when we build the package. The taxonomy dump is not gigantic, but
# there is no need to get it over and over again.
if !isfile(joinpath(taxpath, ".checksum"))
    @info "No local taxonomy checksum found"
    download_dump(archive, chk, "dump")
else
    local_chk = readline(joinpath(taxpath, ".checksum"))
    if local_chk != chk
        @info "Local and remote checksum do not match"
        download_dump(archive, chk, "dump")
    else
        @info "Local taxonomy dump ($(local_chk)) is up to date"
    end
end

@info "Materializing the taxonomy"

# This is the separator used for fields in all the files
tables = joinpath(taxpath, "tables")
ispath(tables) || mkpath(tables)

# Utility functions

include(joinpath(@__DIR__, "..", "src", "types.jl"))

function _class_to_enum(c::T) where {T <: String}
    c = replace(c, " " => "_")
    c = replace(c, "-" => "_")
    return getproperty(@__MODULE__, Symbol("class_$(c)"))
end

"""
This function is responsible for converting a raw line from the taxonomy dump to
a series of correctly typed values. It is quite clearly not the most elegant
version, but it works based on the content of the taxonomy files. As this is
only called at build time, we think that the performance cost is not a big
issue.
"""
function _materialize_data(::Type{T}, v) where {T}
    if v != ""
        T <: Number && return parse(T, v)
        T <: Union{Bool,Missing} && return parse(Bool, v)
        T <: Union{Int,Missing} && return parse(Int, v)
        T <: Symbol && return Symbol(v)
        T <: NCBINameClass && return _class_to_enum(v)
        return v
    else
        return missing
    end
end

function _build_arrow_file(df, dump_file)
    names_columns_types = [eltype(col) for col in eachcol(df)]
    for names_line in readlines(dump_file)
        t = String.(split(names_line, r"\t\|\t?")[1:(end - 1)])
        r = Vector{Any}(undef, length(t))
        for (i, e) in enumerate(t)
            r[i] = _materialize_data(names_columns_types[i], e)
        end
        push!(df, tuple(r...))
    end
    return df
end

# Get the data

@info "Building the names file"
ncbi_names_file_in = joinpath(taxpath, "dump", "names.dmp")
ncbi_names_file_out = joinpath(taxpath, "tables", "names.arrow")
ncbi_names = DataFrames.DataFrame(tax_id=Int[], name=String[], unique_name=Union{String,Missing}[], class=NCBINameClass[])
names_df = _build_arrow_file(ncbi_names, ncbi_names_file_in)
names_df.class = Int.(names_df.class)
Arrow.write(ncbi_names_file_out, names_df)
names_df = nothing
GC.gc()

@info "Building the division file"
ncbi_division_file_in = joinpath(taxpath, "dump", "division.dmp")
ncbi_division_file_out = joinpath(taxpath, "tables", "division.arrow")
ncbi_division = DataFrames.DataFrame(division_id=Int[], division_code=Symbol[], division_name=Symbol[], comments=Union{String,Missing}[])
division_df = _build_arrow_file(ncbi_division, ncbi_division_file_in)
Arrow.write(ncbi_division_file_out, division_df)
division_df = nothing
GC.gc()

@info "Building the nodes file"
ncbi_nodes_file_in = joinpath(taxpath, "dump", "nodes.dmp")
ncbi_nodes_file_out = joinpath(taxpath, "tables", "nodes.arrow")
ncbi_nodes = DataFrames.DataFrame(
    tax_id=Int[], parent_tax_id=Int[],
    rank=Symbol[],
    embl=Union{String,Missing}[],
    division_id=Int[], inherited_div=Union{Bool,Missing}[],
    genetic_code_id=Int[], inherited_gc=Union{Bool,Missing}[],
    mitochondrial_genetic_code_id=Union{Int,Missing}[], inherited_mgc=Union{Bool,Missing}[],
    genbank_hidden=Union{Bool,Missing}[],
    hidden_subtree=Union{Bool,Missing}[],
    comments=Union{String,Missing}[],
    plastid_genetic_code_id=Union{Int,Missing}[], inherited_pgc=Union{Bool,Missing}[],
    specified_species=Union{Bool,Missing}[],
    hydrogenosome_code_id=Union{Int,Missing}[], inherited_hgc=Union{Bool,Missing}[]
    )
nodes_df = _build_arrow_file(ncbi_nodes, ncbi_nodes_file_in)
Arrow.write(ncbi_nodes_file_out, nodes_df)
nodes_df = nothing
GC.gc()


text:
using Test
using PETScX

@testset "lib tests" begin
    # This should error, not library Float64, Float64 is possible!
    @test_throws ErrorException PETScX.getpetsclib(Float64, Float64)

    for petsclib in PETScX.petsclibs
        # Get the library type
        PetscLib = typeof(petsclib)

        # Pull out the parameters
        PetscScalar = PetscLib.parameters[1]
        PetscInt = PetscLib.parameters[2]
        PetscReal = real(PetscScalar)

        # Make sure we can find the library correctly
        @test petsclib === PETScX.getpetsclib(PetscScalar, PetscInt)

        # Check the type accessors
        @test PetscScalar === PETScX.scalartype(petsclib)
        @test PetscScalar === PETScX.scalartype(PetscLib)
        @test PetscInt === PETScX.inttype(petsclib)
        @test PetscInt === PETScX.inttype(PetscLib)
        @test PetscReal === PETScX.realtype(petsclib)
        @test PetscReal === PETScX.realtype(PetscLib)
    end
end

ext: jl


text:
"""
This module supplies some functions for testing of the implementation.
Look at `?` for

* `fdtest` : general finite-difference test
* `fdtest_R2R`: fd test for F : ℝ → ℝ
"""
module Testing

using JuLIP: AbstractCalculator, AbstractAtoms, energy, grad, forces,
         constraint, calculator, set_positions!, dofs, NullConstraint,
         mat, vecs, positions, rattle!
using JuLIP.Potentials: PairPotential, evaluate, evaluate_d, @D
using JuLIP.Constraints: FixedCell

export fdtest

"""
generic finite-difference test for scalar F

* `fdtest(F::Function, dF::Function, x)`

* `fdtest(V::PairPotential, r::Vector{Float64})`

TODO: complete documentation
"""
function fdtest(F::Function, dF::Function, x; verbose=true)
   errors = Float64[]
   E = F(x)
   dE = dF(x)
   # loop through finite-difference step-lengths
   @printf("---------|----------- \n")
   @printf("    h    | error \n")
   @printf("---------|----------- \n")
   for p = 2:11
      h = 0.1^p
      dEh = zeros(dE)
      for n = 1:length(dE)
         x[n] += h
         dEh[n] = (F(x) - E) / h
         x[n] -= h
      end
      push!(errors, vecnorm(dE - dEh, Inf))
      @printf(" %1.1e | %4.2e  \n", h, errors[end])
   end
   @printf("---------|----------- \n")
   if minimum(errors) <= 1e-3 * maximum(errors)
      println("passed")
      return true
   else
      warn("""It seems the finite-difference test has failed, which indicates
      that there is an inconsistency between the function and gradient
      evaluation. Please double-check this manually / visually. (It is
      also possible that the function being tested is poorly scaled.)""")
      return false
   end
end

"finite-difference test for a function V : ℝ → ℝ"
function fdtest_R2R(F::Function, dF::Function, x::Vector{Float64};
                     verbose=true)
   errors = Float64[]
   E = [ F(t) for t in x ]
   dE = [ dF(t) for t in x ]
   # loop through finite-difference step-lengths
   if verbose
      @printf("---------|----------- \n")
      @printf("    h    | error \n")
      @printf("---------|----------- \n")
   end
   for p = 2:11
      h = 0.1^p
      dEh = ([F(t+h) for t in x ] - E) / h
      push!(errors, vecnorm(dE - dEh, Inf))
      if verbose
         @printf(" %1.1e | %4.2e  \n", h, errors[end])
      end
   end
   if verbose
      @printf("---------|----------- \n")
   end
   if minimum(errors) <= 1e-3 * maximum(errors[1:2])
      println("passed")
      return true
   else
      warn("""is seems the finite-difference test has failed, which indicates
            that there is an inconsistency between the function and gradient
            evaluation. Please double-check this manually / visually. (It is
            also possible that the function being tested is poorly scaled.)""")
      return false
   end
end

fdtest(V::PairPotential, r::AbstractVector; kwargs...) =
               fdtest_R2R(s -> V(s), s -> (@D V(s)), collect(r); kwargs...)

function fdtest(calc::AbstractCalculator, at::AbstractAtoms; verbose=true)
   X0 = copy(positions(at))
   # perturb atom positions a bit to get out of equilibrium states
   at = rattle!(at, 0.02)
   # if no constraint is attached, then attach the NullConstraint
   cons = constraint(at)
   if typeof(cons) == NullConstraint
      cons = FixedCell(at)
   end
   # call the actual FD test
   fdtest( x-> energy(calc, set_positions!(cons, at, x)),
           x-> mat(grad(calc, set_positions!(cons, at, x)))[:],
           dofs(at, cons)
         )
   # restore original atom positions
   set_positions!(at, X0)
end

end


text:
# ---
# title: 922. Sort Array By Parity II
# id: problem922
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array, Sort
# link: <https://leetcode.com/problems/sort-array-by-parity-ii/description/>
# hidden: true
# ---
#
# Given an array `A` of non-negative integers, half of the integers in A are
# odd, and half of the integers are even.
#
# Sort the array so that whenever `A[i]` is odd, `i` is odd; and whenever `A[i]`
# is even, `i` is even.
#
# You may return any answer array that satisfies this condition.
#
#
#
# **Example 1:**
#
#
#
#     Input: [4,2,5,7]
#     Output: [4,5,2,7]
#     Explanation: [4,7,2,5], [2,5,4,7], [2,7,4,5] would also have been accepted.
#
#
#
#
# **Note:**
#
#   1. `2 <= A.length <= 20000`
#   2. `A.length % 2 == 0`
#   3. `0 <= A[i] <= 1000`
#
#
#
#
## @lc code=start
using LeetCode

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


text:
using AtomicTables
using Test

@testset "AtomicTables.jl" begin
t = AtomicTables.fakedata(100)
end


text:
function Base.binomial{T<:FloatingPoint}(n::T, k::T)
    exp(lfact(n) - lfact(n - k) - lfact(k))
end

function Base.factorial{T<:FloatingPoint}(n::T, k::T)
    exp(lfact(n) - lfact(k))
end

⊞{T<:Real}(n::T, k::T) = binomial(n, k)
⊠{T<:Real}(n::T, k::T) = factorial(n, n-k)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# consider using the project of this directory for tested versions
using PGFPlotsX, KernelDensity, HighestDensityRegions, Distributions

d = MixtureModel([Normal(-1, 1), Normal(1, 0.2)], [0.6, 0.4])
x = rand(d, 2000)
k = kde(x)
xgrid = range(minimum(x); stop = maximum(x), length = 1000)
ygrid = pdf.(Ref(k), xgrid)
threshold = first(hdr_thresholds([0.5], ygrid))

p = @pgf Axis({ xlabel = "x", ylabel = "kernel density", width = "10cm", height = "8cm" },
              Plot({ very_thick, mesh, point_meta = "explicit" , shader = "interp"},
                   Coordinates(xgrid, ygrid; meta = Int.(threshold .< ygrid))))
pgfsave("1d-density.svg", p)


text:
#! /usr/bin/env julia

using PackageCompiler
import Pkg

out_path = length(ARGS) > 1 ? ARGS[1] : "Baysor";
baysor_path = length(ARGS) > 2 ? ARGS[2] : dirname(@__DIR__);

Base.LOAD_PATH .= baysor_path
create_app(baysor_path, out_path; precompile_execution_file="$(baysor_path)/bin/precompile.jl")

max_stars_repo_name: repushko/Baysor

max_issues_repo_name: repushko/Baysor

max_forks_repo_name: repushko/Baysor


text:
using Test

include("./vector.jl")


text:
module VectorClouds

using LinearAlgebra
using NearestNeighbors
using ..Utils.Mathes

const VectorCloud = Matrix

export VectorCloud

@inline function to_kd_tree(self::VectorCloud; leafsize::Integer=1)
    if isempty(self)
        return nothing
    else
        return KDTree(self[1:3, :], leafsize=leafsize)
    end
end

@inline function to_kd_tree(self::Nothing)
    return nothing
end

"""
    query(self::VectorCloud{T}, targetVectorCloud::VectorCloud;
                targetTree::KDTree = to_kd_tree(targetVectorCloud))
query a VectorCloud to a target KDTree
"""
function query(self::VectorCloud{T}, targetVectorCloud::VectorCloud{T};
               targetTree::KDTree=to_kd_tree(targetVectorCloud)) where T
    queryVectorCloud = self

    idxs, dists = knn(targetTree, queryVectorCloud[1:3, :], 1, false)

    # distance matrix
    ret = Matrix{T}(undef, 2, size(queryVectorCloud, 2))

    for (i, nodeIndexList) in idxs |> enumerate
        targetNodeIndex = nodeIndexList[1]
        # physical distance
        dist = dists[i][1]

        queryVector = queryVectorCloud[4:6, i]
        targetVector = targetVectorCloud[4:6, targetNodeIndex]
        # absolute dot product
        adp = abs( dot( queryVector, targetVector ) )

        ret[:, i] = [dist, adp]
    end
    ret
end

end # end of module VectorClouds


text:
"""
Attribute
"""
struct MangalAttribute
    id::Int64
    name:: AbstractString
    description::AbstractString
    unit::Union{Missing,AbstractString}
end

"""
Reference
"""
struct MangalReference
    id::Int64
    year::Union{Missing,Int64}
    doi::Union{String,Missing}
    jstor::Union{String,Missing}
    pmid::Union{String,Missing}
    bibtex::Union{String,Missing}
    paper::Union{String,Missing}
    data::Union{String,Missing}
end

"""
A `MangalDataset` identifies a collection of networks, possibly containing
a single element. A dataset is identified by its `id` or `name` (both of
which are *unique*).

`name` (`AbstractString`): a unique name describing the dataset.

`public` (`Bool`): indicates whether the dataset details are available to others
than its owner.

`reference` (`Union{Int64,Nothing}`) (*optional*): a reference to the `id` of
the `MangalReference`, or `nothing` if there is no associated reference for this
dataset.

`user` (`Int64`): `id` of the user who added the dataset to the database. This
is *not necessarily* the author of the dataset, see `reference` (and the same
field in the `MangalNetwork`) to get the actual authorship.

`description` (`AbstractString`): a free-form description of the dataset.
"""
struct MangalDataset
    id::Int64
    public::Bool
    name::AbstractString
    created::DateTime
    updated::DateTime
    reference::Union{MangalReference,Missing}
    user::Int64
    description::AbstractString
end

"""
A `MangalNetwork` is a wrapper around *nodes* (and not around interactions, for
reasons which are really not worth mentioning here, but see the documentation
for `MangalNode` for some hints).

`name` (`AbstractString`): a unique name describing the network.

`dataset` (`Int64`): the unique id of the `MangalDataset` to which the network
belongs.

`public` (`Bool`): indicates whether the network details are available to others
than its owner.

`date` (`DateTime`): date and time at which the network was sampled.

`position` (`AbstractGeometry`): the location at which the network was sampled.
This can be any sort of geospatial construct, most notably points *or* polygons.

`complete` (`Bool`): indicates whether the network was sampled completely, or is
a collection of interactions with possible gaps.

`reference` (`Union{Int64,Nothing}`) (*optional*): a reference to the `id` of
the `MangalReference`, or `nothing` if there is no associated reference for this
network.

`user` (`Int64`): `id` of the user who added the network to the database. This
is *not necessarily* the author of the network, see `reference` to get the
actual authorship.

`description` (`AbstractString`): a free-form description of the network.
"""
struct MangalNetwork
    id::Int64
    public::Bool
    name::AbstractString
    date::Union{DateTime,Missing}
    position::Union{AbstractGeometry,Missing}
    created::DateTime
    updated::DateTime
    user::Union{Int64,Missing}
    description::AbstractString
    complete::Bool
    dataset::Union{MangalDataset,Missing}
end

"""
Reference taxon (unique identifier of network nodes)
"""
struct MangalReferenceTaxon
    id::Int64
    name::AbstractString
    bold::Union{Int64,Missing}
    tsn::Union{Int64,Missing}
    ncbi::Union{Int64,Missing}
    eol::Union{Int64,Missing}
    gbif::Union{Int64,Missing}
    created::DateTime
    updated::DateTime
end

"""
Node in a network

The `taxon` field is a `MangalReferenceTaxon` object, so that one can, for
example, query the TSN identifier of a node through `object.taxon.tsn`.

This approach has been chosen because (i) names of nodes in networks can be non
unique and (ii) nodes within the same networks can refer to various taxonomic
levels. As an example, if a network has four distinct nodes identified as
`Ascariasis sp.`, they will represent four nodes in the networks, but map onto
the same `MangalReferenceTaxon` (representing the entire *Ascariasis* genus).
This approach provides a seemless representation of the same taxon across
different networks, but also of the same taxon *within* networks.
"""
struct MangalNode
    id::Int64
    name::AbstractString
    created::DateTime
    updated::DateTime
    taxon::Union{MangalReferenceTaxon,Missing}
end

"""
Interaction
"""
struct MangalInteraction
    id::Int64
    network::MangalNetwork
    from::MangalNode
    to::MangalNode
    date::Union{DateTime,Missing}
    position::Union{AbstractGeometry,Missing}
    directed::Bool
    interaction::Symbol
    method::Union{AbstractString,Missing}
    strength::Union{Number,Missing}
    created::DateTime
    updated::DateTime
    attribute::Union{Missing,MangalAttribute}
end


text:
# borrowed 2-point stencil macros from ParallelStencil.FiniteDifferences; extended to higher-order differences

import ..ParallelStencil: INDICES, WITHIN_DOC
iz, iy, ix = INDICES[1], INDICES[2], INDICES[3]
@static if (_fd.order == 2)
    izi, iyi, ixi = :($iz + 1), :($iy + 1), :($ix + 1)
elseif (_fd.order == 4)
    izi, iyi, ixi = :($iz + 3), :($iy + 3), :($ix + 3)
elseif (_fd.order == 6)
    izi, iyi, ixi = :($iz + 5), :($iy + 5), :($ix + 5)
elseif (_fd.order == 8)
    izi, iyi, ixi = :($iz + 7), :($iy + 7), :($ix + 7)
end

@static if (_fd.order == 2)

    macro within(macroname::String, A::Symbol)
        if macroname == "@all"
            esc(:($iz <= size($A, 1) && $iy <= size($A, 2) && $ix <= size($A, 3)))
        elseif macroname == "@inn"
            esc(
                :(
                    $iz <= size($A, 1) - 2 &&
                    $iy <= size($A, 2) - 2 &&
                    $ix <= size($A, 3) - 2
                ),
            )
        else
            error(
                "unkown macroname: $macroname. If you want to add your own assignement macros, overwrite the macro 'within(macroname::String, A::Symbol)'; to still use the exising macro within as well call ParallelStencil.FiniteDifferences{1|2|3}D.@within(macroname, A) at the end.",
            )
        end
    end

    macro d_za(A::Symbol)
        esc(:($A[$iz+1, $iy, $ix] - $A[$iz, $iy, $ix]))
    end
    macro d_ya(A::Symbol)
        esc(:($A[$iz, $iy+1, $ix] - $A[$iz, $iy, $ix]))
    end
    macro d_xa(A::Symbol)
        esc(:($A[$iz, $iy, $ix+1] - $A[$iz, $iy, $ix]))
    end
    macro d_zi(A::Symbol)
        esc(:($A[$iz+1, $iyi, $ixi] - $A[$iz, $iyi, $ixi]))
    end
    macro d_yi(A::Symbol)
        esc(:($A[$izi, $iy+1, $ixi] - $A[$izi, $iy, $ixi]))
    end
    macro d_xi(A::Symbol)
        esc(:($A[$izi, $iyi, $ix+1] - $A[$izi, $iyi, $ix]))
    end
elseif (_fd.order == 4)

    macro within(macroname::String, A::Symbol)
        if macroname == "@all"
            esc(:($iz <= size($A, 1) && $iy <= size($A, 2) && $ix <= size($A, 3)))
        elseif macroname == "@inn"
            esc(
                :(
                    $iz <= size($A, 1) - 6 &&
                    $iy <= size($A, 2) - 6 &&
                    $ix <= size($A, 3) - 6
                ),
            )
        else
            error(
                "unkown macroname: $macroname. If you want to add your own assignement macros, overwrite the macro 'within(macroname::String, A::Symbol)'; to still use the exising macro within as well call ParallelStencil.FiniteDifferences{1|2|3}D.@within(macroname, A) at the end.",
            )
        end
    end

    macro d_za(A::Symbol)
        esc(
            :(
                $A[$iz+2, $iy, $ix] * 27.0 - $A[$iz+1, $iy, $ix] * 27.0 +
                $A[$iz, $iy, $ix] - $A[$iz+3, $iy, $ix]
            ),
        )
    end
    macro d_ya(A::Symbol)
        esc(
            :(
                $A[$iz, $iy+2, $ix] * 27.0 - $A[$iz, $iy+1, $ix] * 27.0 +
                $A[$iz, $iy, $ix] - $A[$iz, $iy+3, $ix]
            ),
        )
    end
    macro d_xa(A::Symbol)
        esc(
            :(
                $A[$iz, $iy, $ix+2] * 27.0 - $A[$iz, $iy, $ix+1] * 27.0 +
                $A[$iz, $iy, $ix] - $A[$iz, $iy, $ix+3]
            ),
        )
    end
    macro d_zi(A::Symbol)
        esc(
            :(
                $A[$iz+2, $iyi, $ixi] * 27.0 - $A[$iz+1, $iyi, $ixi] * 27.0 +
                $A[$iz, $iyi, $ixi] - $A[$iz+3, $iyi, $ixi]
            ),
        )
    end
    macro d_yi(A::Symbol)
        esc(
            :(
                $A[$izi, $iy+2, $ixi] * 27.0 - $A[$izi, $iy+1, $ixi] * 27.0 +
                $A[$izi, $iy, $ixi] - $A[$izi, $iy+3, $ixi]
            ),
        )
    end
    macro d_xi(A::Symbol)
        esc(
            :(
                $A[$izi, $iyi, $ix+2] * 27.0 - $A[$izi, $iyi, $ix+1] * 27.0 +
                $A[$izi, $iyi, $ix] - $A[$izi, $iyi, $ix+3]
            ),
        )
    end
elseif (_fd.order == 6)
    macro within(macroname::String, A::Symbol)
        if macroname == "@all"
            esc(:($iz <= size($A, 1) && $ix <= size($A, 2)))
        elseif macroname == "@inn"
            esc(
                :(
                    $iz <= size($A, 1) - 10 &&
                    $iy <= size($A, 2) - 10 &&
                    $ix <= size($A, 3) - 10
                ),
            )
        else
            error(
                "unkown macroname: $macroname. If you want to add your own assignement macros, overwrite the macro 'within(macroname::String, A::Symbol)'; to still use the exising macro within as well call ParallelStencil.FiniteDifferences{1|2|3}D.@within(macroname, A) at the end.",
            )
        end
    end

    macro d_za(A::Symbol)
        esc(
            :(
                $A[$iz+3,  $iy,  $ix] * 2250.0 - $A[$iz+2,  $iy,  $ix] * 2250.0 +
                $A[$iz+1,  $iy,  $ix] * 125.0 - $A[$iz+4,  $iy,  $ix] * 125.0 +
                $A[$iz+5,  $iy,  $ix] * 9.0 - $A[$iz,  $iy,  $ix] * 9.0
            ),
        )
    end
    macro d_ya(A::Symbol)
        esc(
            :(
                $A[$iz,  $iy+3,  $ix] * 2250.0 - $A[$iz,  $iy+2,  $ix] * 2250.0 +
                $A[$iz,  $iy+1,  $ix] * 125.0 - $A[$iz,  $iy+4,  $ix] * 125.0 +
                $A[$iz,  $iy+5,  $ix] * 9.0 - $A[$iz,  $iy,  $ix] * 9.0
            ),
        )
    end
    macro d_xa(A::Symbol)
        esc(
            :(
                $A[$iz,  $iy,  $ix+3] * 2250.0 - $A[$iz,  $iy,  $ix+2] * 2250.0 +
                $A[$iz,  $iy,  $ix+1] * 125.0 - $A[$iz,  $iy,  $ix+4] * 125.0 +
                $A[$iz,  $iy,  $ix+5] * 9.0 - $A[$iz,  $iy,  $ix] * 9.0
            ),
        )
    end
    macro d_zi(A::Symbol)
        esc(
            :(
                $A[$iz+3,  $iyi,  $ixi] * 2250.0 - $A[$iz+2,  $iyi,  $ixi] * 2250.0 +
                $A[$iz+1,  $iyi,  $ixi] * 125.0 - $A[$iz+4,  $iyi,  $ixi] * 125.0 +
                $A[$iz+5,  $iyi,  $ixi] * 9.0 - $A[$iz,  $iyi,  $ixi] * 9.0
            ),
        )
    end
    macro d_yi(A::Symbol)
        esc(
            :(
                $A[$izi,  $iy+3,  $ixi] * 2250.0 - $A[$izi,  $iy+2,  $ixi] * 2250.0 +
                $A[$izi,  $iy+1,  $ixi] * 125.0 - $A[$izi,  $iy+4,  $ixi] * 125.0 +
                $A[$izi,  $iy+5,  $ixi] * 9.0 - $A[$izi,  $iy,  $ixi] * 9.0
            ),
        )
    end
    macro d_xi(A::Symbol)
        esc(
            :(
                $A[$izi,  $iyi,  $ix+3] * 2250.0 - $A[$izi,  $iyi,  $ix+2] * 2250.0 +
                $A[$izi,  $iyi,  $ix+1] * 125.0 - $A[$izi,  $iyi,  $ix+4] * 125.0 +
                $A[$izi,  $iyi,  $ix+5] * 9.0 - $A[$izi,  $iyi,  $ix] * 9.0
            ),
        )
    end
elseif (_fd.order == 8)
    macro within(macroname::String, A::Symbol)
        if macroname == "@all"
            esc(:($iz <= size($A, 1) && $ix <= size($A, 2)))
        elseif macroname == "@inn"
            esc(
                :(
                    $iz <= size($A, 1) - 14 &&
                    $iy <= size($A, 2) - 14 &&
                    $ix <= size($A, 3) - 14
                ),
            )
        else
            error(
                "unkown macroname: $macroname. If you want to add your own assignement macros, overwrite the macro 'within(macroname::String, A::Symbol)'; to still use the exising macro within as well call ParallelStencil.FiniteDifferences{1|2|3}D.@within(macroname, A) at the end.",
            )
        end
    end

    macro d_za(A::Symbol)
        esc(
            :(
                $A[$iz+4,  $iy,  $ix] * 128625.0 - $A[$iz+3,  $iy,  $ix] * 128625.0 +
                $A[$iz+2,  $iy,  $ix] * 8575.0 - $A[$iz+5,  $iy,  $ix] * 8575.0 +
                $A[$iz+6,  $iy,  $ix] * 1029.0 - $A[$iz+1,  $iy,  $ix] * 1029.0 +
                $A[$iz,  $iy,  $ix] * 75.0 - $A[$iz+7,  $iy,  $ix] * 75.0
            ),
        )
    end
    macro d_ya(A::Symbol)
        esc(
            :(
                $A[$iz,  $iy+4,  $ix] * 128625.0 - $A[$iz,  $iy+3,  $ix] * 128625.0 +
                $A[$iz,  $iy+2,  $ix] * 8575.0 - $A[$iz,  $iy+5,  $ix] * 8575.0 +
                $A[$iz,  $iy+6,  $ix] * 1029.0 - $A[$iz,  $iy+1,  $ix] * 1029.0 +
                $A[$iz,  $iy,  $ix] * 75.0 - $A[$iz,  $iy+7,  $ix] * 75.0
            ),
        )
    end
    macro d_xa(A::Symbol)
        esc(
            :(
                $A[$iz,  $iy,  $ix+4] * 128625.0 - $A[$iz,  $iy,  $ix+3] * 128625.0 +
                $A[$iz,  $iy,  $ix+2] * 8575.0 - $A[$iz,  $iy,  $ix+5] * 8575.0 +
                $A[$iz,  $iy,  $ix+6] * 1029.0 - $A[$iz,  $iy,  $ix+1] * 1029.0 +
                $A[$iz,  $iy,  $ix] * 75.0 - $A[$iz,  $iy,  $ix+7] * 75.0
            ),
        )
    end
    macro d_zi(A::Symbol)
        esc(
            :(
                $A[$iz+4,  $iyi,  $ixi] * 128625.0 - $A[$iz+3,  $iyi,  $ixi] * 128625.0 +
                $A[$iz+2,  $iyi,  $ixi] * 8575.0 - $A[$iz+5,  $iyi,  $ixi] * 8575.0 +
                $A[$iz+6,  $iyi,  $ixi] * 1029.0 - $A[$iz+1,  $iyi,  $ixi] * 1029.0 +
                $A[$iz,  $iyi,  $ixi] * 75.0 - $A[$iz+7,  $iyi,  $ixi] * 75.0
            ),
        )
    end
    macro d_yi(A::Symbol)
        esc(
            :(
                $A[$izi,  $iy+4,  $ixi] * 128625.0 - $A[$izi,  $iy+3,  $ixi] * 128625.0 +
                $A[$izi,  $iy+2,  $ixi] * 8575.0 - $A[$izi,  $iy+5,  $ixi] * 8575.0 +
                $A[$izi,  $iy+6,  $ixi] * 1029.0 - $A[$izi,  $iy+1,  $ixi] * 1029.0 +
                $A[$izi,  $iy,  $ixi] * 75.0 - $A[$izi,  $iy+7,  $ixi] * 75.0
            ),
        )
    end
    macro d_xi(A::Symbol)
        esc(
            :(
                $A[$izi,  $iyi,  $ix+4] * 128625.0 - $A[$izi,  $iyi,  $ix+3] * 128625.0 +
                $A[$izi,  $iyi,  $ix+2] * 8575.0 - $A[$izi,  $iyi,  $ix+5] * 8575.0 +
                $A[$izi,  $iyi,  $ix+6] * 1029.0 - $A[$izi,  $iyi,  $ix+1] * 1029.0 +
                $A[$izi,  $iyi,  $ix] * 75.0 - $A[$izi,  $iyi,  $ix+7] * 75.0
            ),
        )
    end
end

macro all(A::Symbol)
    esc(:($A[$iz, $iy, $ix]))
end
macro inn(A::Symbol)
    esc(:($A[$izi, $iyi, $ixi]))
end
macro av(A::Symbol)
    esc(
        :(
            (
                $A[$iz, $iy, $ix] +
                $A[$iz+1, $iy, $ix] +
                $A[$iz+1, $iy+1, $ix] +
                $A[$iz+1, $iy+1, $ix+1] +
                $A[$iz, $iy+1, $ix+1] +
                $A[$iz, $iy, $ix+1] +
                $A[$iz+1, $iy, $ix+1] +
                $A[$iz, $iy+1, $ix]
            ) * 0.125
        ),
    )
end
macro av_zi(A::Symbol)
    esc(:(($A[$iz, $iyi, $ixi] + $A[$iz+1, $iyi, $ixi]) * 0.5))
end
macro av_yi(A::Symbol)
    esc(:(($A[$izi, $iy, $ixi] + $A[$izi, $iy+1, $ixi]) * 0.5))
end
macro av_xi(A::Symbol)
    esc(:(($A[$izi, $iyi, $ix] + $A[$izi, $iyi, $ix+1]) * 0.5))
end
macro av_yzi(A::Symbol)
    esc(
        :(
            (
                $A[$iz, $iy, $ixi] +
                $A[$iz+1, $iy, $ixi] +
                $A[$iz, $iy+1, $ixi] +
                $A[$iz+1, $iy+1, $ixi]
            ) * 0.25
        ),
    )
end
macro av_xzi(A::Symbol)
    esc(
        :(
            (
                $A[$iz, $iyi, $ix] +
                $A[$iz+1, $iyi, $ix] +
                $A[$iz, $iyi, $ix+1] +
                $A[$iz+1, $iyi, $ix+1]
            ) * 0.25
        ),
    )
end
macro av_xyi(A::Symbol)
    esc(
        :(
            (
                $A[$izi, $iy, $ix] +
                $A[$izi, $iy+1, $ix] +
                $A[$izi, $iy, $ix+1] +
                $A[$izi, $iy+1, $ix+1]
            ) * 0.25
        ),
    )
end


text:
"""
**Return an interpreted occurrence given its key**

	occurrence(key::Union{String, Integer})

The key can be given as a string or as an integer.
"""
function occurrence(key::Union{String, Integer})
  occ_url = gbifurl * "occurrence/" * string(key)
  occ_key_req = HTTP.get(occ_url)
  return GBIFRecord(JSON.parse(String(occ_key_req.body)))
end

"""
**Retrieve latest occurrences**

    occurrences()

This function will return the latest occurrences -- usually 20, but this is
entirely determined by the server default page size. This is mostly useful to
get a few results rapidly for illustration purposes.
"""
function occurrences()
  occ_s_url = gbifurl * "occurrence/search"
  occ_s_req = HTTP.get(occ_s_url)
  if occ_s_req.status == 200
    body = JSON.parse(String(occ_s_req.body))
    occ = map(GBIFRecord, body["results"])
    maxocc = body["count"] > 200000 ? 200000 : body["count"]
    return GBIFRecords(
      body["offset"],
      maxocc,
      nothing,
      occ,
      ones(Bool, length(occ))
    )
  else
    @error "Non-OK status returned: $(occ_s_req.status)"
  end
end

"""
**Retrieve latest occurrences based on a query**

    occurrences(q::Dict)

Returns occurrences that correspond to a filter, given in `q` as a dictionary.
When first called, this function will return the latest 20 hits (or whichever
default page size GBIF uses). Future occurrences can be queried with `next!` or
`complete!`.
"""
function occurrences(q::Dict)
  check_records_parameters!(q)
  occ_s_url = gbifurl * "occurrence/search"
  occ_s_req = HTTP.get(occ_s_url, query=q)
  if occ_s_req.status == 200
    body = JSON.parse(String(occ_s_req.body))
    occ = map(GBIFRecord, body["results"])
    maxocc = body["count"] > 200000 ? 200000 : body["count"]
    return GBIFRecords(
      body["offset"],
      maxocc,
      q,
      occ,
      ones(Bool, length(occ))
    )
  end
end

"""
**Retrieve latest occurrences for a taxon based on a query**

    occurrences(t::GBIFTaxon, q::Dict)

Returns occurrences that correspond to a filter, given in `q` as a dictionary.
When first called, this function will return the latest 20 hits (or whichever
default page size GBIF uses). Future occurrences can be queried with `next!` or
`complete!`.
"""
function occurrences(t::GBIFTaxon, q::Dict)
	levels = [:kingdom, :phylum, :class, :order, :family, :genus, :species]
	for l in levels
		if getfield(t, l) !== nothing
			q[String(l)*"Key"] = getfield(t, l).second
		end
	end
	return occurrences(q)
end

"""
**Retrieve latest occurrences for a taxon**

    occurrences(t::GBIFTaxon)

Returns occurrences that correspond to a filter, given in `q` as a dictionary.
When first called, this function will return the latest 20 hits (or whichever
default page size GBIF uses). Future occurrences can be queried with `next!` or
`complete!`.
"""
function occurrences(t::GBIFTaxon)
	return occurrences(t, Dict())
end


text:
# This file is part of the TaylorSeries.jl Julia package, MIT license
#
# Luis Benet & David P. Sanders
# UNAM
#
# MIT Expat license
#

## Constructors ##
doc"""
    Taylor1{T<:Number} <: Number

DataType for polynomial expansions in one independent variable.

**Fields:**

- `coeffs :: Array{T,1}` Expansion coefficients; the $i$-th
component is the coefficient of degree $i-1$ of the expansion.
- `order  :: Int64` Maximum order (degree) of the polynomial.
"""
immutable Taylor1{T<:Number} <: Number
    coeffs :: Array{T,1}
    order :: Int

    ## Inner constructor ##
    function Taylor1(coeffs::Array{T,1}, order::Int)
        lencoef = length(coeffs)
        order = max(order, lencoef-1)
        if order == lencoef-1
            return new(coeffs, order)
        else
            resize!(coeffs, order+1)
            for i = lencoef+1:order+1
                coeffs[i] = zero(T)
            end
            return new(coeffs, order)
        end
    end
end

## Outer constructors ##
Taylor1{T<:Number}(x::Taylor1{T}, order::Int) = Taylor1{T}(x.coeffs, order)
Taylor1{T<:Number}(x::Taylor1{T}) = x
Taylor1{T<:Number}(coeffs::Array{T,1}, order::Int) = Taylor1{T}(coeffs, order)
Taylor1{T<:Number}(coeffs::Array{T,1}) = Taylor1{T}(coeffs, length(coeffs)-1)
Taylor1{T<:Number}(x::T, order::Int) = Taylor1{T}([x], order)
# Taylor1{T<:Number}(x::T) = Taylor1{T}([x], 0)

# Shortcut to define Taylor1 independent variables
"""
    Taylor1(T, [order=1])
    Taylor1([order=1])

Shortcut to define the independent variable of a `Taylor1{T}` polynomial of
given `order`. If the type `T` is ommitted, `Float64` is assumend.
"""
Taylor1{T<:Number}(::Type{T}, order::Int=1) = Taylor1{T}( [zero(T), one(T)], order)
Taylor1(order::Int=1) = Taylor1(Float64, order)

## get_coeff ##
"""
    get_coeff(a, n)

Return the coefficient of order `n::Int` of a `a::Taylor1` polynomial.
"""
get_coeff(a::Taylor1, n::Int) = (@assert 0 <= n <= a.order+1;
    return a.coeffs[n+1])

## Type, length ##
eltype{T<:Number}(::Taylor1{T}) = T
length{T<:Number}(a::Taylor1{T}) = a.order

## Conversion and promotion rules ##
convert{T<:Number}(::Type{Taylor1{T}}, a::Taylor1) =
    Taylor1(convert(Array{T,1}, a.coeffs), a.order)
function convert{T<:Integer, S<:AbstractFloat}(::Type{Taylor1{Rational{T}}},
    a::Taylor1{S})
    v = Array(Rational{T}, length(a.coeffs))
    for i in eachindex(v)
        # v[i] = convert(Rational{T}, a.coeffs[i])
        v[i] = rationalize(a.coeffs[i], tol=eps(one(S)))
    end
    Taylor1(v)
end
convert{T<:Number}(::Type{Taylor1{T}}, b::Array{T,1}) = Taylor1(b, length(b)-1)
convert{T<:Number, S<:Number}(::Type{Taylor1{T}}, b::Array{S,1}) =
    Taylor1(convert(Array{T,1},b))
convert{T<:Number, S<:Number}(::Type{Taylor1{T}}, b::S) = Taylor1([convert(T,b)], 0)
convert{T<:Number}(::Type{Taylor1{T}}, b::T) = Taylor1([b], 0)

promote_rule{T<:Number}(::Type{Taylor1{T}}, ::Type{Taylor1{T}}) = Taylor1{T}
promote_rule{T<:Number, S<:Number}(::Type{Taylor1{T}}, ::Type{Taylor1{S}}) =
    Taylor1{promote_type(T, S)}
promote_rule{T<:Number}(::Type{Taylor1{T}}, ::Type{Array{T,1}}) = Taylor1{T}
promote_rule{T<:Number, S<:Number}(::Type{Taylor1{T}}, ::Type{Array{S,1}}) =
    Taylor1{promote_type(T, S)}
promote_rule{T<:Number}(::Type{Taylor1{T}}, ::Type{T}) = Taylor1{T}
promote_rule{T<:Number, S<:Number}(::Type{Taylor1{T}}, ::Type{S}) =
    Taylor1{promote_type(T, S)}

## Auxiliary function ##
function firstnonzero{T<:Number}(ac::Vector{T})
    nonzero::Int = length(ac)
    for i in eachindex(ac)
        if ac[i] != zero(T)
            nonzero = i-1
            break
        end
    end
    nonzero
end
firstnonzero{T<:Number}(a::Taylor1{T}) = firstnonzero(a.coeffs)

function fixshape{T<:Number}(a::Taylor1{T}, b::Taylor1{T})
    if a.order == b.order
        return a, b
    elseif a.order < b.order
        return Taylor1(a, b.order), b
    end
    return a, Taylor1(b, a.order)
end

function fixshape{T<:Number, S<:Number}(a::Taylor1{T}, b::Taylor1{S})
    fixshape(promote(a,b)...)
end

## real, imag, conj and ctranspose ##
for f in (:real, :imag, :conj)
    @eval ($f){T<:Number}(a::Taylor1{T}) = Taylor1(($f)(a.coeffs), a.order)
end
ctranspose{T<:Number}(a::Taylor1{T}) = conj(a)

## zero and one ##
zero{T<:Number}(a::Taylor1{T}) = Taylor1(zero(T), a.order)
one{T<:Number}(a::Taylor1{T}) = Taylor1(one(T), a.order)

## Equality ##
function ==(a::Taylor1, b::Taylor1)
    a, b = fixshape(a, b)
    return a.coeffs == b.coeffs
end

# Tests `isinf` and `isnan` for all polynomial coefficients
for f in (:isinf, :isnan)
    @eval begin
        function ($f)(a::Taylor1)
            test = false
            for i in eachindex(a.coeffs)
                @inbounds test = ($f)(a.coeffs[i])
                test && break
            end
            test
        end
    end
end

## Addition and substraction ##
for f in (:+, :-)
    @eval begin
        function ($f)(a::Taylor1, b::Taylor1)
            a, b = fixshape(a, b)
            v = similar(a.coeffs)
            @simd for i in eachindex(a.coeffs)
                @inbounds v[i] = ($f)(a.coeffs[i], b.coeffs[i])
            end
            return Taylor1(v, a.order)
        end
        function ($f)(a::Taylor1)
            v = similar(a.coeffs)
            @simd for i in eachindex(v)
                @inbounds v[i] = ($f)(a.coeffs[i])
            end
            return Taylor1(v, a.order)
        end
        function ($f)(a::Taylor1, b::Union{Real,Complex})
            @inbounds aux = ($f)(a.coeffs[1], b)
            v = Array(typeof(aux), length(a.coeffs))
            @simd for i in eachindex(v)
                @inbounds v[i] = a.coeffs[i]
            end
            @inbounds v[1] = aux
            Taylor1(v, a.order)
        end
        function ($f)(a::Union{Real,Complex}, b::Taylor1)
            @inbounds aux = ($f)(a, b.coeffs[1])
            v = Array(typeof(aux), length(b.coeffs))
            @simd for i in eachindex(v)
                @inbounds v[i] = ($f)(b.coeffs[i])
            end
            @inbounds v[1] = aux
            Taylor1(v, b.order)
        end
    end
end

## Multiplication ##
*(a::Bool, b::Taylor1) = *(promote(a,b)...)
*(a::Taylor1, b::Bool) = b*a
function *(a::Union{Real,Complex}, b::Taylor1)
    @inbounds aux = a * b.coeffs[1]
    v = Array(typeof(aux), length(b.coeffs))
    @simd for i in eachindex(v)
        @inbounds v[i] = a * b.coeffs[i]
    end
    Taylor1(v, b.order)
end
*(a::Taylor1, b::Union{Real,Complex}) = b * a
doc"""
```
*(a, b)
```

Return the Taylor expansion of $a \cdot b$, of order `max(a.order,b.order)`, for
`a::Taylor1`, `b::Taylor1` polynomials.

For details on making the Taylor expansion, see [`TaylorSeries.mulHomogCoef`](@ref).
"""
function *(a::Taylor1, b::Taylor1)
    a, b = fixshape(a, b)
    coeffs = similar(a.coeffs)
    @inbounds coeffs[1] = a.coeffs[1] * b.coeffs[1]
    @inbounds for k = 1:a.order
        coeffs[k+1] = mulHomogCoef(k, a.coeffs, b.coeffs)
    end
    Taylor1(coeffs, a.order)
end

# Homogeneous coefficient for the multiplication
doc"""
    mulHomogCoef(kcoef, ac, bc)

Compute the `k`-th expansion coefficient of $c = a\cdot b$ given by

\begin{equation*}
c_k = \sum_{j=0}^k a_j b_{k-j},
\end{equation*}

with $a$ and $b$ `Taylor1` polynomials.

Inputs are the `kcoef`-th coefficient, and the vectors of the expansion coefficients
`ac` and `bc`, corresponding respectively to `a` and `b`.
"""
function mulHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, bc::Array{T,1})
    kcoef == 0 && return ac[1] * bc[1]
    coefhomog = zero(T)
    @inbounds for i = 0:kcoef
        coefhomog += ac[i+1] * bc[kcoef-i+1]
    end
    coefhomog
end

## Division ##
/{T<:Real}(a::Taylor1, b::T) = a * inv(b)
/{T<:Complex}(a::Taylor1, b::T) = a * inv(b)
doc"""
```
/(a, b)
```

Return the Taylor expansion of $a/b$, of order `max(a.order,b.order)`, for
`a::Taylor1`, `b::Taylor1` polynomials.

For details on making the Taylor expansion, see [`TaylorSeries.divHomogCoef`](@ref).
"""
function /(a::Taylor1, b::Taylor1)
    a, b = fixshape(a, b)
    # order and coefficient of first factorized term
    orddivfact, cdivfact = divfactorization(a, b)
    T = typeof(cdivfact)
    v1 = convert(Array{T,1}, a.coeffs)
    v2 = convert(Array{T,1}, b.coeffs)
    coeffs = zeros(T, a.order+1)
    @inbounds coeffs[1] = cdivfact
    @inbounds for k = orddivfact+1:a.order
        coeffs[k-orddivfact+1] = divHomogCoef(k, v1, v2, coeffs, orddivfact)
    end
    Taylor1(coeffs, a.order)
end

function divfactorization(a1::Taylor1, b1::Taylor1)
    # order of first factorized term; a1 and b1 assumed to be of the same order
    a1nz = firstnonzero(a1)
    b1nz = firstnonzero(b1)
    orddivfact = min(a1nz, b1nz)
    if orddivfact > a1.order
        orddivfact = a1.order
    end
    cdivfact = a1.coeffs[orddivfact+1] / b1.coeffs[orddivfact+1]

    # Is the polynomial factorizable?
    if isinf(cdivfact) || isnan(cdivfact)
        throw(ArgumentError(
        """Division does not define a Taylor1 polynomial
        or its first non-zero coefficient is Inf/NaN.
        Order k=$(orddivfact) => coeff[$(orddivfact+1)]=$(cdivfact)."""))
    end

    return orddivfact, cdivfact
end

# Homogeneous coefficient for the division
doc"""
    divHomogCoef(kcoef, ac, bc, coeffs, ordfact)

Compute the `k-th` expansion coefficient of $c = a / b$ given by

\begin{equation*}
c_k =  \frac{1}{b_0} (a_k - \sum_{j=0}^{k-1} c_j b_{k-j}),
\end{equation*}

with $a$ and $b$ `Taylor1` polynomials.

Inputs are the `kcoef`-th coefficient, the vectors of the expansion coefficients
`ac` and `bc`, corresponding respectively to `a` and `b`, the
already calculated expansion coefficients `coeffs` of `c`, and `ordfact`
which is the order of the factorized term of the denominator,
whenever `b_0` is zero.
"""
function divHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, bc::Array{T,1},
    coeffs::Array{T,1}, ordfact::Int)
    #
    @inbounds kcoef == ordfact && return ac[ordfact+1] / bc[ordfact+1]
    coefhomog = mulHomogCoef(kcoef, coeffs, bc)
    @inbounds coefhomog = (ac[kcoef+1]-coefhomog) / bc[ordfact+1]
    coefhomog
end

## Division functions: rem and mod
for op in (:mod, :rem)
    @eval begin
        function ($op){T<:Real, S<:Real}(a::Taylor1{T}, x::S)
            coeffs = copy(a.coeffs)
            @inbounds coeffs[1] = ($op)(a.coeffs[1], x)
            return Taylor1(coeffs, a.order)
        end
    end
end

function mod2pi{T<:Real}(a::Taylor1{T})
    coeffs = copy(a.coeffs)
    @inbounds coeffs[1] = mod2pi( a.coeffs[1] )
    return Taylor1( coeffs, a.order)
end

## abs function ##
"""
    abs(a)

Return `a` or `-a` depending on the 0-th order coefficient of the
`Taylor1` polynomial `a`.
If `a.coeffs[1]` is zero, an `ArgumentError` is thrown.
"""
function abs{T<:Real}(a::Taylor1{T})
    if a.coeffs[1] > zero(T)
        return a
    elseif a.coeffs[1] < zero(T)
        return -a
    else
        throw(ArgumentError(
        """The 0th order Taylor1 coefficient must be non-zero
        (abs(x) is not differentiable at x=0)."""))
    end
end

## Int power ##
doc"""
```
^(a, x)
```

Return the Taylor expansion of $a^x$ for `a::Taylor1` polynomial and `x::Number`.
If `x` is non integer and the 0-th order coefficient is zero, an
`ArgumentError` is thrown.
"""
function ^{T<:Number}(a::Taylor1{T}, n::Integer)
    n == 0 && return one(a)
    n == 1 && return a
    n == 2 && return square(a)
    n < 0 && return inv( a^(-n) )
    return power_by_squaring(a, n)
end

function ^{T<:Integer}(a::Taylor1{T}, n::Integer)
    n == 0 && return one(a)
    n == 1 && return a
    n == 2 && return square(a)
    n < 0 && throw(DomainError())
    return power_by_squaring(a, n)
end

## power_by_squaring; slightly modified from base/intfuncs.jl
## Licensed under MIT "Expat"
function power_by_squaring(x::Taylor1, p::Integer)
    p == 1 && return copy(x)
    p == 0 && return one(x)
    p == 2 && return square(x)
    t = trailing_zeros(p) + 1
    p >>= t

    while (t -= 1) > 0
        x = square(x)
    end

    y = x
    while p > 0
        t = trailing_zeros(p) + 1
        p >>= t
        while (t -= 1) >= 0
            x = square(x)
        end
        y *= x
    end

    return y
end

## Rational power ##
^{T<:Integer}(a::Taylor1,x::Rational{T}) = a^(x.num/x.den)

^(a::Taylor1, b::Taylor1) = exp( b*log(a) )

## Real power ##
function ^{S<:Real}(a::Taylor1, x::S)
    x == zero(x) && return one(a)
    x == one(x)/2 && return sqrt(a)
    isinteger(x) && return a^round(Int,x)

    # First non-zero coefficient
    l0nz = firstnonzero(a)
    l0nz > a.order && return zero(a)

    # The first non-zero coefficient of the result; must be integer
    lnull = x*l0nz
    !isinteger(lnull) &&
        throw(ArgumentError(
        """The 0th order Taylor1 coefficient must be non-zero
        to raise the Taylor1 polynomial to a non-integer exponent."""))

    # Reaching this point, it is possible to implement the power of the Taylor1
    # polynomial. The last l0nz coefficients are set to zero.
    lnull = trunc(Int,lnull)
    #l0nz > 0 && warn("The last k=$(l0nz) Taylor1 coefficients ARE SET to 0.")
    @inbounds aux = (a.coeffs[l0nz+1])^x
    T = typeof(aux)
    v = convert(Array{T,1}, a.coeffs)
    coeffs = zeros(T, a.order+1)
    @inbounds coeffs[lnull+1] = aux
    k0 = lnull+l0nz
    @inbounds for k = k0+1:a.order
        coeffs[k-l0nz+1] = powHomogCoef(k, v, x, coeffs, l0nz)
    end

    Taylor1(coeffs,a.order)
end
^{T<:Complex}(a::Taylor1, x::T) = exp( x*log(a) )

# Homogeneous coefficients for real power
doc"""
    powHomogCoef(kcoef, ac, x, coeffs, knull)

Compute the `k-th` expansion coefficient of $c = a^x$, given by

\begin{equation*}
c_k = \frac{1}{k a_0} \sum_{j=0}^{k-1} ((k-j)x -j) - j)a_{k-j} c_j,
\end{equation*}

with $a$ a `Taylor1` polynomial, and `x` a number.

Inputs are the `kcoef`-th coefficient, the vector of the expansion coefficients
`ac` of `a`, the exponent `x`, the already calculated expansion
coefficients `coeffs` of `c`, and `knull`
which is the order of the first non-zero coefficient of `a`.
"""
function powHomogCoef{T<:Number, S<:Real}(kcoef::Int, ac::Array{T,1}, x::S,
    coeffs::Array{T,1}, knull::Int)

    kcoef == knull && return (ac[knull+1])^x
    coefhomog = zero(T)
    for i = 0:kcoef-knull-1
        aux = x*(kcoef-i)-i
        @inbounds coefhomog += aux*ac[kcoef-i+1]*coeffs[i+1]
    end
    aux = kcoef - knull*(x+1)
    @inbounds coefhomog = coefhomog / (aux*ac[knull+1])

    coefhomog
end

## Square ##
function square(a::Taylor1)
    coeffs = Array(eltype(a), a.order+1)
    coeffs[1] = a.coeffs[1]^2
    @inbounds for k = 1:a.order
        coeffs[k+1] = squareHomogCoef(k, a.coeffs)
    end
    Taylor1(coeffs,a.order)
end

# Homogeneous coefficients for square
doc"""
    squareHomogCoef(kcoef, ac)

Compute the `k-th` expansion coefficient of $c = a^2$, given by

\begin{eqnarray*}
c_k &=& 2 \sum_{j=0}^{(k-1)/2} a_{k-j} a_j, \text{ if $k$ is odd, or }  \\\\
c_k &=& 2 \sum_{j=0}^{(k-2)/2} a_{k-j} a_j + (a_{k/2})^2, \text{ if $k$ is even, }
\end{eqnarray*}

with $a$ a `Taylor1` polynomial.

Inputs are the `kcoef`-th coefficient and the vector of the expansion coefficients
`ac` of `a`.
"""
function squareHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1})
    kcoef == 0 && return ac[1]^2
    coefhomog = zero(T)
    kodd = kcoef%2
    kend = div(kcoef - 2 + kodd, 2)
    @inbounds for i = 0:kend
        coefhomog += ac[i+1]*ac[kcoef-i+1]
    end
    coefhomog = convert(T,2) * coefhomog
    if kodd == 0
        @inbounds coefhomog += ac[div(kcoef,2)+1]^2
    end
    coefhomog
end

## Square root ##
doc"""
    sqrt(a)

Return  the Taylor expansion of $\sqrt(a)$, of order `a.order`,
for `a::Taylor1` polynomial.
If the first non-vanishing coefficient of `a` is an odd power, and
`ArgumentError` is thrown.

For details on making the Taylor expansion, see [`TaylorSeries.sqrtHomogCoef`](@ref).
"""
function sqrt(a::Taylor1)
    # First non-zero coefficient
    l0nz = firstnonzero(a)
    if l0nz > a.order
        return zero(a)
    elseif l0nz%2 == 1 # l0nz must be pair
        throw(ArgumentError(
        """First non-vanishing Taylor1 coefficient must correspond
        to an **even power** in order to expand `sqrt` around 0."""))
    end

    # Reaching this point, it is possible to implement the sqrt of the Taylor1 polynomial.
    # The last l0nz coefficients are set to zero.
    ##l0nz > 0 && warn("The last k=$(l0nz) Taylor1 coefficients ARE SET to 0.")
    lnull = div(l0nz, 2)
    @inbounds aux = sqrt(a.coeffs[l0nz+1])
    T = typeof(aux)
    v = convert(Array{T,1}, a.coeffs)
    coeffs = zeros(T, a.order+1)
    @inbounds coeffs[lnull+1] = aux
    @inbounds for k = lnull+1:a.order-l0nz
        coeffs[k+1] = sqrtHomogCoef(k, v, coeffs, lnull)
    end
    Taylor1(coeffs, a.order)
end

# Homogeneous coefficients for the square-root
doc"""
    sqrtHomogCoef(kcoef, ac, coeffs, knull)

Compute the `k-th` expansion coefficient of $c = \sqrt(a)$, given by

\begin{eqnarray*}
c_k &=& \frac{1}{2 c_0} ( a_k - 2 \sum_{j=0}^{(k-1)/2} c_{k-j}c_j),
\text{ if $k$ is odd, or } \\\\
c_k &=& \frac{1}{2 c_0} ( a_k - 2 \sum_{j=0}^{(k-2)/2} c_{k-j}c_j) - (c_{k/2})^2,
\text{ if $k$ is even,}
\end{eqnarray*}

with $a$ a `Taylor1` polynomial.

Inputs are the `kcoef`-th coefficient, the vector of the expansion coefficients
`ac` of `a`, the already calculated expansion coefficients `coeffs` of `c`,
and `knull`, which is half of the order of the first non-zero coefficient of `a`.
"""
function sqrtHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, coeffs::Array{T,1}, knull::Int)
    kcoef == knull && return sqrt(ac[2*knull+1])
    coefhomog = zero(T)
    kodd = (kcoef - knull)%2
    kend = div(kcoef - knull - 2 + kodd, 2)
    two = convert(T,2)
    @inbounds for i = knull+1:knull+kend
        coefhomog += coeffs[i+1]*coeffs[kcoef+knull-i+1]
    end
    @inbounds aux = ac[kcoef+knull+1] - two*coefhomog

    if kodd == 0
        @inbounds aux = aux - (coeffs[kend+knull+2])^2
    end
    @inbounds coefhomog = aux / (two*coeffs[knull+1])

    coefhomog
end

## Exp ##
doc"""
    exp(a)

Return the Taylor expansion of $e^a$, of order `a.order`, for
`a::Taylor1` polynomial.

For details on making the Taylor expansion, see [`TaylorSeries.expHomogCoef`](@ref).
"""
function exp(a::Taylor1)
    @inbounds aux = exp( a.coeffs[1] )
    T = typeof(aux)
    v = convert(Array{T,1}, a.coeffs)
    coeffs = similar(v)
    @inbounds coeffs[1] = aux
    @inbounds for k = 1:a.order
        coeffs[k+1] = expHomogCoef(k, v, coeffs)
    end
    Taylor1( coeffs, a.order )
end

# Homogeneous coefficients for exp
doc"""
    expHomogCoef(kcoef, ac, coeffs)

Compute the `k-th` expansion coefficient of $c = \exp(a)$ given by

\begin{equation*}
c_k = \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} c_j,
\end{equation*}

with $a$ a `Taylor1` polynomial.

Inputs are the `kcoef`-th coefficient, the vector of the expansion coefficients
`ac` of $a$, and the already calculated expansion coefficients `coeffs` of `c`.
"""
function expHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, coeffs::Array{T,1})
    kcoef == 0 && return exp(ac[1])
    coefhomog = zero(T)
    @inbounds for i = 0:kcoef-1
        coefhomog += (kcoef-i) * ac[kcoef-i+1] * coeffs[i+1]
    end
    coefhomog = coefhomog/kcoef
    coefhomog
end

## Log ##
doc"""
    log(a)

Return the Taylor expansion of $\log(a)$, of order `a.order`, for `a::Taylor1` polynomial.

For details on making the Taylor expansion, see [`TaylorSeries.logHomogCoef`](@ref).
"""
function log(a::Taylor1)
    ( firstnonzero(a)>0 ) && throw(
        ArgumentError("Impossible to expand `log` around 0."))
    @inbounds aux = log( a.coeffs[1] )
    T = typeof(aux)
    ac = convert(Array{T,1}, a.coeffs)
    coeffs = similar(ac)
    @inbounds coeffs[1] = aux
    @inbounds for k = 1:a.order
        coeffs[k+1] = logHomogCoef(k, ac, coeffs)
    end
    Taylor1( coeffs, a.order )
end

# Homogeneous coefficients for log
doc"""
    logHomogCoef(kcoef, ac, coeffs)

Compute the `k-th` expansion coefficient of $c = \log(a)$, given by

\begin{equation*}
c_k = \frac{1}{a_0} (a_k - \frac{1}{k} \sum_{j=0}^{k-1} j a_{k-j} c_j ),
\end{equation*}

with $a$ a `Taylor1` polynomial.

Inputs are the `kcoef`-th coefficient, the vector of the expansion coefficients
`ac` of `a`, and the already calculated expansion coefficients `coeffs` of `c`.
"""
function logHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, coeffs::Array{T,1})
    kcoef == 0 && return log( ac[1] )
    coefhomog = zero(T)
    @inbounds for i = 1:kcoef-1
        coefhomog += (kcoef-i) * ac[i+1] * coeffs[kcoef-i+1]
    end
    @inbounds coefhomog = (ac[kcoef+1] -coefhomog/kcoef) / ac[1]
    coefhomog
end

### TRIGONOMETRIC FUNCTIONS ###

## Sin ##
doc"""
    sin(a)

Return the Taylor expansion of $\sin(a)$, of order `a.order`, for
`a::Taylor1` polynomial.

For details on making the Taylor expansion, see [`TaylorSeries.sincosHomogCoef`](@ref).
"""
sin(a::Taylor1) = sincos(a)[1]

## Cos ##
doc"""
    cos(a)

Return the Taylor expansion of $\cos(a)$, of order `a.order`,
for `a::Taylor1` polynomial

For details on making the Taylor expansion, see [`TaylorSeries.sincosHomogCoef`](@ref).
"""
cos(a::Taylor1) = sincos(a)[2]

## Sin and Cos ##
function sincos(a::Taylor1)
    @inbounds aux = sin( a.coeffs[1] )
    T = typeof(aux)
    v = convert(Array{T,1}, a.coeffs)
    sincoeffs = similar(v)
    coscoeffs = similar(v)
    @inbounds sincoeffs[1] = aux
    @inbounds coscoeffs[1] = cos( a.coeffs[1] )
    @inbounds for k = 1:a.order
        sincoeffs[k+1], coscoeffs[k+1] = sincosHomogCoef(k, v, sincoeffs, coscoeffs)
    end
    return Taylor1(sincoeffs, a.order), Taylor1(coscoeffs, a.order)
end

# Homogeneous coefficients for sincos
doc"""
    sincosHomogCoef(kcoef, ac, scoeffs, ccoeffs)

Compute the `k-th` expansion coefficient of $s = \sin(a)$ and $c=\cos(a)$
simultaneously given by

\begin{eqnarray*}
s_k &=& \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} c_j \\\ \\

c_k &=& -\frac{1}{k}\sum_{j=0}^{k-1} (k-j) a_{k-j} s_j
\end{eqnarray*}

with $a$ a `Taylor1` polynomial.

Inputs are the `kcoef`-th coefficient, the vector of the expansion coefficients
`ac` of `a`, and the already calculated expansion coefficients `scoeffs`
and `ccoeffs` of `sin(a)` and `cos(a)`, respectvely.
"""
function sincosHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1},
    scoeffs::Array{T,1}, ccoeffs::Array{T,1})

    kcoef == 0 && return sin( ac[1] ), cos( ac[1] )
    sincoefhom = zero(T)
    coscoefhom = zero(T)

    @inbounds for i = 1:kcoef
        x = i * ac[i+1]
        sincoefhom += x * ccoeffs[kcoef-i+1]
        coscoefhom -= x * scoeffs[kcoef-i+1]
    end

    sincoefhom = sincoefhom/kcoef
    coscoefhom = coscoefhom/kcoef
    return sincoefhom, coscoefhom
end

## tan ##
doc"""
    tan(a)

Return the Taylor expansion of $\tan(a)$, of order `a.order`, for
`a::Taylor1` polynomial.

For details on making the Taylor expansion, see [`TaylorSeries.tanHomogCoef`](@ref).
"""
function tan(a::Taylor1)
    aux = tan( a.coeffs[1] )
    T = typeof(aux)
    v = convert(Array{T,1}, a.coeffs)
    coeffs = similar(v)
    coeffst2 = similar(v)
    @inbounds coeffs[1] = aux
    @inbounds coeffst2[1] = aux^2
    @inbounds for k = 1:a.order
        coeffs[k+1] = tanHomogCoef(k, v, coeffst2)
        coeffst2[k+1] = squareHomogCoef(k, coeffs)
    end
    Taylor1( coeffs, a.order )
end

# Homogeneous coefficients for tan
doc"""
    tanHomogCoef(kcoef, ac, coeffst2)

Compute the `k-th` expansion coefficient of $c = \tan(a)$ given by

\begin{equation*}
c_k = a_k + \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} p_j,
\end{equation*}

with $a$ a `Taylor1` polynomial and $p = c^2$.

Inputs are the `kcoef`-th coefficient, the vector of the expansion coefficients
`ac` of `a`, and the already calculated expansion coefficients `coeffst2`
of `c^2`.
"""
function tanHomogCoef{T<:Number}(kcoef::Int,ac::Array{T,1},coeffst2::Array{T,1})
    kcoef == 0 && return tan( ac[1] )
    coefhomog = zero(T)
    @inbounds for i = 0:kcoef-1
        coefhomog += (kcoef-i)*ac[kcoef-i+1]*coeffst2[i+1]
    end
    @inbounds coefhomog = ac[kcoef+1] + coefhomog/kcoef
    coefhomog
end

### INVERSE TRIGONOMETRIC FUNCTIONS ###

## Arcsin ##
doc"""
    asin(a)

Return the Taylor expansion of $\arcsin(a)$, of order `a.order`, for
`a::Taylor1` polynomial.

For details on making the Taylor expansion, see [`TaylorSeries.asinHomogCoef`](@ref).
"""
function asin(a::Taylor1)
    @inbounds aux = asin( a.coeffs[1] )
    T = typeof(aux)
    ac = convert(Array{T,1}, a.coeffs)
    ac[1]^2 == one(T) && throw(ArgumentError("""
        Recursion formula diverges due to vanishing `sqrt`."""))
    rc = sqrt(one(T) - a^2).coeffs
    coeffs = zeros(ac)
    @inbounds coeffs[1] = aux
    @inbounds for k in 1:a.order
        coeffs[k+1] = asinHomogCoef(k, ac, rc, coeffs)
    end
    Taylor1( coeffs, a.order )
end

# Homogeneous coefficients for arcsin
doc"""
    asinHomogCoef(kcoef, ac, rc, coeffs)

Compute the `k-th` expansion coefficient of $s = \arcsin(a)$ given by

\begin{equation*}
s_k = \frac{1}{ \sqrt{r_0} } \big( a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} s_j \big),
\end{equation*}

with $a$ a `Taylor1` polynomial and $r = \sqrt{1 - a^2}$.

Inputs are the `kcoef`-th coefficient, the vector of the expansion coefficients
`ac` of `a`, the already calculated expansion coefficients `rc`
of $r$ (see above), and the already calculated expansion coefficients
`coeffs` of `asin(a)`.
"""
function asinHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, rc::Array{T,1},
    coeffs::Array{T,1})
    kcoef == 0 && return asin( ac[1] )
    coefhomog = zero(T)
    @inbounds for i in 1:kcoef-1
        coefhomog += (kcoef-i) * rc[i+1] * coeffs[kcoef-i+1]
    end
    @inbounds coefhomog = (ac[kcoef+1] - coefhomog/kcoef) / rc[1]
    coefhomog
end

## Arccos ##
doc"""
    acos(a)

Return the Taylor expansion of $\arccos(a)$, of order `a.order`, for
`a::Taylor1` polynomial.

For details on making the Taylor expansion, see [`TaylorSeries.acosHomogCoef`](@ref).
"""
function acos(a::Taylor1)
    @inbounds aux = asin( a.coeffs[1] )
    T = typeof(aux)
    ac = convert(Array{T,1}, a.coeffs)
    ac[1]^2 == one(T) && throw(ArgumentError("""
        Recursion formula diverges due to vanishing `sqrt`."""))
    rc = sqrt(one(T) - a^2).coeffs
    coeffs = zeros(ac)
    @inbounds coeffs[1] = aux
    @inbounds for k in 1:a.order
        coeffs[k+1] = asinHomogCoef(k , ac, rc, coeffs)
    end
    @inbounds coeffs[1] = -acos( a.coeffs[1] )
    Taylor1( -coeffs, a.order )
end

# Homogeneous coefficients for arccos
doc"""
    acosHomogCoef(kcoef, ac, rc, coeffs)

Compute the `k-th` expansion coefficient of $c = \arccos(a)$ given by

\begin{equation*}
c_k = - \frac{1}{ r_0 } \big( a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} c_j \big),
\end{equation*}

with $a$ a `Taylor1` polynomial and $r = \sqrt{1 - a^2}$.

Inputs are the `kcoef`-th coefficient, the vector of the expansion coefficients
`ac` of `a`, the already calculated expansion coefficients `rc`
of $r$ (see above), and the already calculated expansion coefficients
`coeffs` of `acos(a)`.
"""
function acosHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, rc::Array{T,1}, coeffs::Array{T,1})
    kcoef == 0 && return asin( ac[1] )
    coefhomog = zero(T)
    @inbounds for i in 1:kcoef-1
        coefhomog += (kcoef-i) * rc[i+1] * coeffs[kcoef-i+1]
    end
    @inbounds coefhomog = -(ac[kcoef+1] - coefhomog/kcoef) / rc[1]
    coefhomog
end

## Arctan
doc"""
    atan(a)

Return the Taylor expansion of $\arctan(a)$, of order `a.order`, for
`a::Taylor1` polynomial.

For details on making the Taylor expansion, see [`TaylorSeries.atanHomogCoef`](@ref).
"""
function atan(a::Taylor1)
    @inbounds aux = atan( a.coeffs[1] )
    T = typeof(aux)
    rc = (one(T) + a^2).coeffs
    rc[1] == zero(T) && throw(ArgumentError("""
        Recursion formula has a pole."""))
    ac = convert(Array{T,1}, a.coeffs)
    coeffs = zeros(ac)
    @inbounds coeffs[1] = aux
    @inbounds for k in 1:a.order
        coeffs[k+1] = atanHomogCoef(k , ac, rc, coeffs)
    end
    Taylor1( coeffs, a.order )
end

# Homogeneous coefficients for arctan
doc"""
    atanHomogCoef(kcoef, ac, rc, coeffs)

Compute the `k-th` expansion coefficient of $c = \arctan(a)$ given by

\begin{equation*}
t_k = \frac{1}{r_0}(a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} t_j) ,
\end{equation*}

with $a$ a `Taylor1` polynomial and $r = 1 + a^2$.

Inputs are the `kcoef`-th coefficient, the vector of the expansion coefficients
`ac` of `a`, the already calculated expansion coefficients `rc`
of $r$ (see above), and the already calculated expansion coefficients
`coeffs` of `asin(a)`.
"""
function atanHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, rc::Array{T,1}, coeffs::Array{T,1})
    kcoef == 0 && return atan( ac[1] )
    coefhomog = zero(T)
    @inbounds for i in 1:kcoef-1
        coefhomog += (kcoef-i) * rc[i+1] * coeffs[kcoef-i+1]
    end
    @inbounds coefhomog = (ac[kcoef+1] - coefhomog/kcoef) / rc[1]
    coefhomog
end

### HYPERBOLIC FUNCTIONS ###

## Sinh ##
doc"""
    sinh(a)

Return the Taylor expansion of $\sinh(a)$, of order `a.order`, for
`a::Taylor1` polynomial.

For details on making the Taylor expansion, see [`TaylorSeries.sinhcoshHomogCoef`](@ref).
"""
sinh(a::Taylor1) = sinhcosh(a)[1]

## Cosh ##
doc"""
    cosh(a)

Return the Taylor expansion of $\cosh(a)$, of order `a.order`,
for `a::Taylor1` polynomial

For details on making the Taylor expansion, see [`TaylorSeries.sinhcoshHomogCoef`](@ref).
"""
cosh(a::Taylor1) = sinhcosh(a)[2]

## Sinh and Cosh ##
function sinhcosh(a::Taylor1)
    @inbounds aux = sinh( a.coeffs[1] )
    T = typeof(aux)
    v = convert(Array{T,1}, a.coeffs)
    sinhcoeffs = similar(v)
    coshcoeffs = similar(v)
    @inbounds sinhcoeffs[1] = aux
    @inbounds coshcoeffs[1] = cosh( a.coeffs[1] )
    @inbounds for k = 1:a.order
        sinhcoeffs[k+1], coshcoeffs[k+1] = sinhcoshHomogCoef(k, v, sinhcoeffs, coshcoeffs)
    end
    return Taylor1(sinhcoeffs, a.order), Taylor1(coshcoeffs, a.order)
end

# Homogeneous coefficients for sinhcosh
doc"""
    sinhcoshHomogCoef(kcoef, ac, scoeffs, ccoeffs)

Compute the `k-th` expansion coefficient of $sh = \sinh(a)$ and $ch=\cosh(a)$
simultaneously given by

\begin{eqnarray*}
sh_k &=& \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} ch_j \\\ \\

ch_k &=& \frac{1}{k}\sum_{j=0}^{k-1} (k-j) a_{k-j} sh_j
\end{eqnarray*}

with $a$ a `Taylor1` polynomial.

Inputs are the `kcoef`-th coefficient, the vector of the expansion coefficients
`ac` of `a`, and the already calculated expansion coefficients `shcoeffs`
and `chcoeffs` of `sinh(a)` and `cosh(a)`, respectvely.
"""
function sinhcoshHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1},
    shcoeffs::Array{T,1}, chcoeffs::Array{T,1})

    kcoef == 0 && return sinh( ac[1] ), cosh( ac[1] )
    sinhcoefhom = zero(T)
    coshcoefhom = zero(T)

    @inbounds for i = 1:kcoef
        x = i * ac[i+1]
        sinhcoefhom += x * chcoeffs[kcoef-i+1]
        coshcoefhom += x * shcoeffs[kcoef-i+1]
    end

    sinhcoefhom = sinhcoefhom/kcoef
    coshcoefhom = coshcoefhom/kcoef
    return sinhcoefhom, coshcoefhom
end

## tanh ##
doc"""
    tanh(a)

Return the Taylor expansion of $\tanh(a)$, of order `a.order`, for
`a::Taylor1` polynomial.

For details on making the Taylor expansion, see [`TaylorSeries.tanhHomogCoef`](@ref).
"""
function tanh(a::Taylor1)
    aux = tanh( a.coeffs[1] )
    T = typeof(aux)
    v = convert(Array{T,1}, a.coeffs)
    coeffs = similar(v)
    coeffst2 = similar(v)
    @inbounds coeffs[1] = aux
    @inbounds coeffst2[1] = aux^2
    @inbounds for k = 1:a.order
        coeffs[k+1] = tanhHomogCoef(k, v, coeffst2)
        coeffst2[k+1] = squareHomogCoef(k, coeffs)
    end
    Taylor1( coeffs, a.order )
end

# Homogeneous coefficients for tanh
doc"""
    tanhHomogCoef(kcoef, ac, coeffst2)

Compute the `k-th` expansion coefficient of $c = \tanh(a)$ given by

\begin{equation*}
th_k = a_k - \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} p_j,
\end{equation*}

with $a$ a `Taylor1` polynomial and $p = th^2$.

Inputs are the `kcoef`-th coefficient, the vector of the expansion coefficients
`ac` of `a`, and the already calculated expansion coefficients `coeffst2`
of `th^2`.
"""
function tanhHomogCoef{T<:Number}(kcoef::Int,ac::Array{T,1},coeffst2::Array{T,1})
    kcoef == 0 && return tanh( ac[1] )
    coefhomog = zero(T)
    @inbounds for i = 0:kcoef-1
        coefhomog += (kcoef-i)*ac[kcoef-i+1]*coeffst2[i+1]
    end
    @inbounds coefhomog = ac[kcoef+1] - coefhomog/kcoef
    coefhomog
end

## Differentiating ##
"""
    derivative(a)

Return the `Taylor1` polynomial of the differential of `a::Taylor1`.
The last coefficient is set to zero.
"""
function derivative(a::Taylor1)
    coeffs = zero(a.coeffs)
    @inbounds coeffs[1] = a.coeffs[2]
    @inbounds for i = 1:a.order
        coeffs[i] = i*a.coeffs[i+1]
    end
    return Taylor1(coeffs, a.order)
end

"""
    derivative(n, a)

Return the value of the `n`-th derivative of the polynomial `a`.
"""
function derivative{T<:Number}(n::Int, a::Taylor1{T})
    @assert a.order >= n >= 0
    factorial( widen(n) ) * a.coeffs[n+1] :: T
end

## Integrating ##
"""
    integrate(a, [x])

Return the integral of `a::Taylor1`. The constant of integration
(0-th order coefficient) is set to `x`, which is zero if ommitted.
"""
function integrate{T<:Number, S<:Number}(a::Taylor1{T}, x::S)
    R = promote_type(T, typeof(a.coeffs[1] / 1), S)
    coeffs = zeros(R, a.order+1)
    @inbounds for i = 1:a.order
        coeffs[i+1] = a.coeffs[i] / i
    end
    @inbounds coeffs[1] = convert(R, x)
    return Taylor1(coeffs, a.order)
end
integrate{T<:Number}(a::Taylor1{T}) = integrate(a, zero(T))

## Evaluating ##
"""
    evaluate(a, [dx])

Evaluate a `Taylor1` polynomial using Horner's rule (hand coded). If `dx` is
ommitted, its value is considered as zero.
"""
function evaluate{T<:Number}(a::Taylor1{T}, dx::T)
    @inbounds suma = a.coeffs[end]
    @inbounds for k = a.order:-1:1
        suma = suma*dx + a.coeffs[k]
    end
    suma
end
function evaluate{T<:Number,S<:Number}(a::Taylor1{T}, dx::S)
    R = promote_type(T,S)
    @inbounds suma = convert(R, a.coeffs[end])
    @inbounds for k = a.order:-1:1
        suma = suma*dx + a.coeffs[k]
    end
    suma
end
evaluate{T<:Number}(a::Taylor1{T}) = a.coeffs[1]

doc"""
    evaluate(x, δt)

Evaluates each element of `x::Array{Taylor1{T},1}`, representing
the dependent variables of an ODE, at *time* δt.
"""
function evaluate{T<:Number, S<:Number}(x::Array{Taylor1{T},1}, δt::S)
    R = promote_type(T,S)
    return evaluate(convert(Array{Taylor1{R},1},x), convert(R,δt))
end
function evaluate{T<:Number}(x::Array{Taylor1{T},1}, δt::T)
    xnew = Array{T}( length(x) )
    evaluate!(x, δt, xnew)
    return xnew
end
evaluate{T<:Number}(a::Array{Taylor1{T},1}) = evaluate(a, zero(T))

doc"""
    evaluate!(x, δt, x0)

Evaluates each element of `x::Array{Taylor1{T},1}`,
representing the Taylor expansion for the dependent variables
of an ODE at *time* δt. It updates the vector `x0` with the
computed values.
"""
function evaluate!{T<:Number}(x::Array{Taylor1{T},1}, δt::T, x0::Array{T,1})
    @assert length(x) == length(x0)
    @inbounds for i in eachindex(x)
        x0[i] = evaluate( x[i], δt )
    end
    nothing
end

"""
    evaluate(a, x)

Substitute `x::Taylor1` as independent variable in a `a::Taylor1` polynomial.
"""
function evaluate{T<:Number,S<:Number}(a::Taylor1{T}, x::Taylor1{S})
    a, x = fixshape(a, x)
    @inbounds suma = a.coeffs[end]
    @inbounds for k = a.order:-1:1
        suma = suma*x + a.coeffs[k]
    end
    suma
end

"""
    A_mul_B!(Y, A, B)

Multiply A*B and save the result in Y.
"""
function A_mul_B!{T<:Number}(y::Vector{Taylor1{T}}, a::Union{Matrix{T},SparseMatrixCSC{T}},
    b::Vector{Taylor1{T}})

    n, k = size(a)
    @assert (length(y)== n && length(b)== k)

    # determine the maximal order of b
    order = maximum([b1.order for b1 in b])

    # Use matrices of coefficients (of proper size) and A_mul_B!
    B = zeros(T, k, order+1)
    for i = 1:k
        @inbounds ord = b[i].order
        @inbounds for j = 1:ord+1
            B[i,j] = b[i].coeffs[j]
        end
    end
    Y = Array(T, n, order+1)
    A_mul_B!(Y, a, B)
    @inbounds for i = 1:n
        y[i] = Taylor1( collect(Y[i,:]), order)
    end

    return y
end


text:
module PUDLData

using DataDeps

export pudl_datasets

# `consts` providing info for `datadeps/release*.jl`
include("v1_0_0.jl")
include("v1_1_0.jl")
include("version_info.jl")
include("datadep_helpers.jl")

function __init__()
    datadep_dir = joinpath(@__DIR__, "datadeps")
    for filename in readdir(datadep_dir)
        include(joinpath(datadep_dir, filename))
    end
end

end  # module


text:
# This file was generated, do not modify it. # hide
High = ifelse.(carseats.Sales .<= 8, "No", "Yes") |> categorical;
carseats[!, :High] = High;


text:
using ClusterValidityIndices
using Test
using Logging
using Random

# Set the log level
LogLevel(Logging.Info)

# Common code tests for coverage
@testset "common.jl" begin
    # Load the common methods
    include("../src/common.jl")
    # Load the data
    data_path = "../data/correct_partition.csv"
    data, labels = get_cvi_data(data_path)
     # Sort it into sequential order
    data, labels = sort_cvi_data(data, labels)
    # Verify that the labels are now monotonic
    @test issorted(labels)
end

# Test all of the provided example CVI scripts
@testset "Example Scripts" begin
    # Switch to the top for execution because our scripts point to the datasets
    # relative to themselves, not relative to the test dir
    @info "Test directory" pwd()
    cd("../")
    @info "Switching working directory to top for running scripts" pwd()

    # Run scripts
    include("../src/examples/combined.jl")
    include("../src/examples/db.jl")
    include("../src/examples/ps.jl")
    include("../src/examples/xb.jl")
end


text:
# this is all of the tests - not everything is included in runtests.jl because it depends on unregistered packages

include("build.jl")

include("runtests.jl")
include("lightdark.jl")


text:
function quadform(x::Value, A::AbstractExpr; assume_psd=false)
    return x' * A * x
end

"""
    is_psd(A; tol)

Check whether `A` is positive semi-definite by computing a LDLᵀ factorization of `A + tol*I`
"""
function is_psd(A; tol=sqrt(eps(float(real(eltype(A))))))
    T = eltype(A)
    # If `A` is neither a Matrix nor SparseMatrixCSC, we do the following:
    # * sparse fallback if the arithmetic is supported
    # * dense fallack otherwise
    if T <: AbstractFloat || T <: LinearAlgebra.BlasFloat
        return is_psd(sparse(A); tol=tol)
    else
        return is_psd(Matrix(A); tol=tol)
    end
end
is_psd(A::SparseMatrixCSC{Complex{T}}; tol::T=sqrt(eps(T))) where{T<:LinearAlgebra.BlasReal} = isposdef(A + tol*I)
function is_psd(A::SparseMatrixCSC{T}; tol::T=sqrt(eps(T))) where{T<:AbstractFloat}
    # LDLFactorizations requires the input matrix to only have the upper triangle.
    A_ = Symmetric(sparse(UpperTriangular(A)) + tol*I)
    try
        F = ldl(A_)
        d = F.D.diag
        return minimum(d) >= 0
    catch err
        # If the matrix could not be factorized, then it is not PSD
        isa(err, LDLFactorizations.SQDException) && return false
        rethrow()  # Something else happened
    end
end
is_psd(A::Matrix; tol=sqrt(eps(float(real(eltype(A)))))) = isposdef(A + tol*I)

function quadform(x::AbstractExpr, A::Value; assume_psd=false)
    if length(size(A)) != 2 || size(A, 1) != size(A, 2)
        error("Quadratic form only takes square matrices")
    end
    if !ishermitian(A)
        error("Quadratic form only defined for Hermitian matrices")
    end
    if assume_psd
        factor = 1
    else
        if is_psd(A)
            factor = 1
        elseif is_psd(-A)
            factor = -1
        else
            error("Quadratic forms supported only for semidefinite matrices")
        end
    end

    P = sqrt(Hermitian(factor * A))
    return factor * square(norm2(P * x))
end

"""
    quadform(x::AbstractExpr, A::AbstractExpr; assume_psd=false)

Represents `x' * A * x` where either:

* `x` is a vector-valued variable and `A` is a positive semidefinite or negative semidefinite matrix (and in particular Hermitian or real symmetric). If `assume_psd=true`, then
`A` will be assumed to be positive semidefinite. Otherwise, `Convex.is_psd` will be used to check if `A` is positive semidefinite or negative semidefinite.
* or `A` is a matrix-valued variable and `x` is a vector.

"""
function quadform(x::AbstractExpr, A::AbstractExpr; assume_psd=false)
    if vexity(x) == ConstVexity()
        return quadform(evaluate(x), A; assume_psd=assume_psd)
    elseif vexity(A) == ConstVexity()
        return quadform(x, evaluate(A); assume_psd=assume_psd)
    else
        error("Either `x` or `A` must be constant in `quadform(x,A)`.")
    end
end


text:
include("sndintpr.jl")
Scheme.repl()

max_stars_repo_name: carl-gjerris/pseudo-scheme-interpreter

max_issues_repo_name: carl-gjerris/pseudo-scheme-interpreter

max_forks_repo_name: carl-gjerris/pseudo-scheme-interpreter


text:
using .C2BasisDefinition: C1, C1Norm

function opnormbound(N::Type{<:C1}, M, B)
    n = size(M, 2)

    est = 0

    for i in 1:n
        v = zeros(n)
        v[i] = 1
        z = C1Norm(B, v)
        est+= C1Norm(B, M*v)/z
    end
    return est.hi
end

function norms_of_powers_basis(B::Basis, m::Integer, Q::DiscretizedOperator, f::AbstractArray;
    normv0::Real=-1., #used as "missing" value
    normQ::Real=-1.,
    normE::Real=-1.,
    normEF::Real=-1.,
    normIEF::Real=-1.,
    normN::Real=-1.)

    @assert eltype(f) <: Interval
    T = typeof(zero(eltype(Q.L)).hi) # gets "Float64" from Q.L
    n = size(Q.L, 1)
    M = mid.(Q.L)
    R = radius.(Q.L)
    #δ = opnormbound(N, R)
    #γz = gamma(T, max_nonzeros_per_row(Q.L))
    #γn = gamma(T, n+3) # not n+2 like in the paper, because we wish to allow for f to be the result of rounding
    #ϵ = zero(T)
    midf = mid.(f)

    norm(v) = C2BasisDefinition.C1Norm(B, v)

    norms = zeros(m)

    S = zeros((n, m))

    for v in AverageZero(B)
        norm_0 = norm(v)
        v/= norm_0.lo
        for i in 1:m
            w = M*v
            v = w - Q.e * (midf*w)[1]
            S[:, i]+= abs.(v)
        end
    end

    return [norm(S[:, i]) for i in 1:m]

end

max_stars_repo_name: orkolorko/ResponseManneville

max_issues_repo_name: orkolorko/ResponseManneville

max_forks_repo_name: orkolorko/ResponseManneville


text:
function show(b::UInt64)
    border = "\n   " * '▄'^18 * "\n"
    @printf(Crayon(reset=true, foreground=:dark_gray), "%s", border)
    dark = false
    for rank = 8:-1:1
        @printf(Crayon(reset=true, foreground=:dark_gray), " %s █", rank)
        for file = 1:8
            i = (8-rank)*8+file
            if b & INT2UINT[i] != EMPTY
                @printf(Crayon(reset=true, foreground=:red, bold=true,
                               background=(dark ? :dark_gray : :light_gray)), "%s", 'x')
            else
                @printf(Crayon(reset=true, foreground=(dark ? :dark_gray : :light_gray)), "█")
            end
            @printf(Crayon(reset=true, foreground=(dark ? :dark_gray : :light_gray)), "█")
            file != 8 ? dark = !dark : continue
        end
        @printf(Crayon(reset=true, foreground=:dark_gray), "█")
        rank != 1 ? @printf("\n") : continue
    end
    lower_border = "\n   " * '▀'^18 * "\n"
    @printf("%s", lower_border)
    @printf("    a b c d e f g h\n")
end

function getSymbol(cs::ChessSet, sq::UInt64)
    if cs.P & sq != EMPTY
        return 'o'
    elseif cs.N & sq != EMPTY
        return 'ζ'
    elseif cs.B & sq != EMPTY
        return 'Δ'
    elseif cs.R & sq != EMPTY
        return 'Π'
    elseif cs.Q & sq != EMPTY
        return 'Ψ'
    elseif cs.K & sq != EMPTY
        return '+'
    end
    return ' '
end

function Base.show(io::IO, b::Board)
    frame = :dark_gray
    dark_square = :blue
    light_square = :cyan
    border = "\n   " * '▄'^18 * "\n"
    @printf(Crayon(reset=true, foreground=frame), "%s", border)
    dark = false
    for rank = 8:-1:1
        @printf(Crayon(reset=true, foreground=frame), " %s █", rank)
        for file = 1:8
            i = (8-rank)*8+file
            if b.taken & INT2UINT[i] != EMPTY
                c = getSymbol(b.white, INT2UINT[i])
                if c != ' '
                    color = :red
                else
                    color = :dark_gray
                    c = getSymbol(b.black, INT2UINT[i])
                end
                @printf(Crayon(reset=true, foreground=color, bold=true,
                               background=(dark ? dark_square : light_square)), "%s", c)
            else
                @printf(Crayon(reset=true, foreground=(dark ? dark_square : light_square)), "█")
            end
            @printf(Crayon(reset=true, foreground=(dark ? dark_square : light_square)), "█")
            file != 8 ? dark = !dark : continue
        end
        @printf(Crayon(reset=true, foreground=frame), "█")
        rank != 1 ? @printf("\n") : continue
    end
    lower_border = "\n   " * '▀'^18 * "\n"
    @printf("%s", lower_border)
    @printf("    a b c d e f g h\n")
end

function plainPrint(b::Board)
    border = "\n   " * '▄'^18 * "\n"
    @printf("%s", border)
    dark = false
    for rank = 8:-1:1
        @printf(" %s █", rank)
        for file = 1:8
            i = (8-rank)*8+file
            if b.taken & INT2UINT[i] != EMPTY
                c = getSymbol(b.white, INT2UINT[i])
                if c != ' '
                    color = :red
                else
                    color = :dark_gray
                    c = getSymbol(b.black, INT2UINT[i])
                end
                @printf("%s", c)
            else
                @printf("⋅")
            end
            @printf(" ")
            file != 8 ? dark = !dark : continue
        end
        @printf("█")
        rank != 1 ? @printf("\n") : continue
    end
    lower_border = "\n   " * '▀'^18 * "\n"
    @printf("%s", lower_border)
    @printf("    a b c d e f g h\n")
end

function loadFen(fen::String)
    board, active_str, castling_str, enpassant, halfmove, fullmove = split(fen)
    active = active_str == "w"
    castling = UInt8(0)
    if castling_str != "-"
        if occursin("K", castling_str); castling |= CK end
        if occursin("Q", castling_str); castling |= CQ end
        if occursin("k", castling_str); castling |= Ck end
        if occursin("q", castling_str); castling |= Cq end
    end
    enpassant != "-" ? enpassant = PGN2UINT[enpassant] : enpassant = EMPTY

    pieces = Dict{Char,UInt64}()
    for c in ['P', 'N', 'B', 'R', 'Q', 'K']
        push!(pieces, c=>EMPTY)
        push!(pieces, lowercase(c)=>EMPTY)
    end

    i = 1
    str_i = 1
    while str_i <= length(board)
        c = board[str_i]
        if isnumeric(c)
            i += parse(Int64, c)
        elseif c != '/'
            pieces[c] |= INT2UINT[i]
            i += 1
        end
        str_i += 1
    end

    whites = EMPTY
    blacks = EMPTY
    for c in ['P', 'N', 'B', 'R', 'Q', 'K']
        whites |= pieces[c]
        blacks |= pieces[lowercase(c)]
    end

    white = ChessSet(pieces['P'], pieces['N'], pieces['B'],
                     pieces['R'], pieces['Q'], pieces['K'],
                     whites)
    black = ChessSet(pieces['p'], pieces['n'], pieces['b'],
                     pieces['r'], pieces['q'], pieces['k'],
                     blacks)

    taken = EMPTY
    taken = whites |= blacks

    halfmove = parse(Int64, halfmove)
    fullmove = parse(Int64, fullmove)

    return Board(white, black, taken, active, castling, enpassant, halfmove, fullmove)
end

function setBoard()
    return loadFen("rnbqkbnr/pppppppp/8/8/8/8/PPPPPPPP/RNBQKBNR w KQkq - 0 1")
end

function getTypeAt(pieces::ChessSet, square::UInt64)
    if pieces.P & square != EMPTY
        return :pawn
    elseif pieces.N & square != EMPTY
        return :knight
    elseif pieces.B & square != EMPTY
        return :bishop
    elseif pieces.R & square != EMPTY
        return :rook
    elseif pieces.Q & square != EMPTY
        return :queen
    elseif pieces.K & square != EMPTY
        return :king
    else
        return :none
    end
end

function bitshift(p::UInt64)
    rx = 64-trailing_zeros(p)
    return INT2UINT[rx], abs(leading_zeros(p)+1-rx)
end


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

# some test strings
astr = "Hello, world.\n"
u8str = "∀ ε > 0, ∃ δ > 0: |x-y| < δ ⇒ |f(x)-f(y)| < ε"

# I think these should give error on 4 also, and "" is not treated
# consistently with SubString("",1,1), nor with Char[]
for ind in (0, 5)
    @test_throws BoundsError search("foo", SubString("",1,1), ind)
    @test_throws BoundsError rsearch("foo", SubString("",1,1), ind)
    @test_throws BoundsError searchindex("foo", SubString("",1,1), ind)
    @test_throws BoundsError rsearchindex("foo", SubString("",1,1), ind)
end

# Note: the commented out tests will be enabled after fixes to make
# sure that search/rsearch/searchindex/rsearchindex are consistent
# no matter what type of AbstractString the second argument is
@test_throws BoundsError search("foo", Char[], 0)
@test_throws BoundsError search("foo", Char[], 5)
# @test_throws BoundsError rsearch("foo", Char[], 0)
@test_throws BoundsError rsearch("foo", Char[], 5)

# @test_throws BoundsError searchindex("foo", Char[], 0)
# @test_throws BoundsError searchindex("foo", Char[], 5)
# @test_throws BoundsError rsearchindex("foo", Char[], 0)
# @test_throws BoundsError rsearchindex("foo", Char[], 5)

# @test_throws ErrorException in("foobar","bar")
@test_throws BoundsError search(b"\x1\x2",0x1,0)
@test rsearchindex(b"foo",b"o",0) == 0
@test rsearchindex(SubString("",1,1),SubString("",1,1)) == 1

@test search(b"foo",'o') == 2
@test rsearch(b"foo",'o') == 3
@test search(b"foó",'ó') == 3
@test rsearch(b"foó",'ó') == 3

# ascii search
for str in [astr, GenericString(astr)]
    @test_throws BoundsError search(str, 'z', 0)
    @test_throws BoundsError search(str, '∀', 0)
    @test search(str, 'x') == 0
    @test search(str, '\0') == 0
    @test search(str, '\u80') == 0
    @test search(str, '∀') == 0
    @test search(str, 'H') == 1
    @test search(str, 'l') == 3
    @test search(str, 'l', 4) == 4
    @test search(str, 'l', 5) == 11
    @test search(str, 'l', 12) == 0
    @test search(str, ',') == 6
    @test search(str, ',', 7) == 0
    @test search(str, '\n') == 14
    @test search(str, '\n', 15) == 0
    @test_throws BoundsError search(str, 'ε', nextind(str,endof(str))+1)
    @test_throws BoundsError search(str, 'a', nextind(str,endof(str))+1)
end

# ascii rsearch
for str in [astr]
    @test rsearch(str, 'x') == 0
    @test rsearch(str, '\0') == 0
    @test rsearch(str, '\u80') == 0
    @test rsearch(str, '∀') == 0
    @test rsearch(str, 'H') == 1
    @test rsearch(str, 'H', 0) == 0
    @test rsearch(str, 'l') == 11
    @test rsearch(str, 'l', 5) == 4
    @test rsearch(str, 'l', 4) == 4
    @test rsearch(str, 'l', 3) == 3
    @test rsearch(str, 'l', 2) == 0
    @test rsearch(str, ',') == 6
    @test rsearch(str, ',', 5) == 0
    @test rsearch(str, '\n') == 14
end

# utf-8 search
for str in (u8str, GenericString(u8str))
    @test_throws BoundsError search(str, 'z', 0)
    @test_throws BoundsError search(str, '∀', 0)
    @test search(str, 'z') == 0
    @test search(str, '\0') == 0
    @test search(str, '\u80') == 0
    @test search(str, '∄') == 0
    @test search(str, '∀') == 1
    @test_throws UnicodeError search(str, '∀', 2)
    @test search(str, '∀', 4) == 0
    @test search(str, '∃') == 13
    @test_throws UnicodeError search(str, '∃', 15)
    @test search(str, '∃', 16) == 0
    @test search(str, 'x') == 26
    @test search(str, 'x', 27) == 43
    @test search(str, 'x', 44) == 0
    @test search(str, 'δ') == 17
    @test_throws UnicodeError search(str, 'δ', 18)
    @test search(str, 'δ', nextind(str,17)) == 33
    @test search(str, 'δ', nextind(str,33)) == 0
    @test search(str, 'ε') == 5
    @test search(str, 'ε', nextind(str,5)) == 54
    @test search(str, 'ε', nextind(str,54)) == 0
    @test search(str, 'ε', nextind(str,endof(str))) == 0
    @test search(str, 'a', nextind(str,endof(str))) == 0
    @test_throws BoundsError search(str, 'ε', nextind(str,endof(str))+1)
    @test_throws BoundsError search(str, 'a', nextind(str,endof(str))+1)
end

# utf-8 rsearch
for str in [u8str]
    @test rsearch(str, 'z') == 0
    @test rsearch(str, '\0') == 0
    @test rsearch(str, '\u80') == 0
    @test rsearch(str, '∄') == 0
    @test rsearch(str, '∀') == 1
    @test rsearch(str, '∀', 0) == 0
    @test rsearch(str, '∃') == 13
    @test rsearch(str, '∃', 14) == 13
    @test rsearch(str, '∃', 13) == 13
    @test rsearch(str, '∃', 12) == 0
    @test rsearch(str, 'x') == 43
    @test rsearch(str, 'x', 42) == 26
    @test rsearch(str, 'x', 25) == 0
    @test rsearch(str, 'δ') == 33
    @test rsearch(str, 'δ', 32) == 17
    @test rsearch(str, 'δ', 16) == 0
    @test rsearch(str, 'ε') == 54
    @test rsearch(str, 'ε', 53) == 5
    @test rsearch(str, 'ε', 4) == 0
end

# string search with a single-char string
@test search(astr, "x") == 0:-1
@test search(astr, "H") == 1:1
@test search(astr, "H", 2) == 0:-1
@test search(astr, "l") == 3:3
@test search(astr, "l", 4) == 4:4
@test search(astr, "l", 5) == 11:11
@test search(astr, "l", 12) == 0:-1
@test search(astr, "\n") == 14:14
@test search(astr, "\n", 15) == 0:-1

@test search(u8str, "z") == 0:-1
@test search(u8str, "∄") == 0:-1
@test search(u8str, "∀") == 1:1
@test search(u8str, "∀", 4) == 0:-1
@test search(u8str, "∃") == 13:13
@test search(u8str, "∃", 16) == 0:-1
@test search(u8str, "x") == 26:26
@test search(u8str, "x", 27) == 43:43
@test search(u8str, "x", 44) == 0:-1
@test search(u8str, "ε") == 5:5
@test search(u8str, "ε", 7) == 54:54
@test search(u8str, "ε", 56) == 0:-1

# string rsearch with a single-char string
@test rsearch(astr, "x") == 0:-1
@test rsearch(astr, "H") == 1:1
@test rsearch(astr, "H", 2) == 1:1
@test rsearch(astr, "H", 0) == 0:-1
@test rsearch(astr, "l") == 11:11
@test rsearch(astr, "l", 10) == 4:4
@test rsearch(astr, "l", 4) == 4:4
@test rsearch(astr, "l", 3) == 3:3
@test rsearch(astr, "l", 2) == 0:-1
@test rsearch(astr, "\n") == 14:14
@test rsearch(astr, "\n", 13) == 0:-1

@test rsearch(u8str, "z") == 0:-1
@test rsearch(u8str, "∄") == 0:-1
@test rsearch(u8str, "∀") == 1:1
@test rsearch(u8str, "∀", 0) == 0:-1
#TODO: setting the limit in the middle of a wide char
#      makes search fail but rsearch succeed.
#      Should rsearch fail as well?
#@test rsearch(u8str, "∀", 2) == 0:-1 # gives 1:3
@test rsearch(u8str, "∃") == 13:13
@test rsearch(u8str, "∃", 12) == 0:-1
@test rsearch(u8str, "x") == 43:43
@test rsearch(u8str, "x", 42) == 26:26
@test rsearch(u8str, "x", 25) == 0:-1
@test rsearch(u8str, "ε") == 54:54
@test rsearch(u8str, "ε", 53) == 5:5
@test rsearch(u8str, "ε", 4) == 0:-1

# string search with a single-char regex
@test search(astr, r"x") == 0:-1
@test search(astr, r"H") == 1:1
@test search(astr, r"H", 2) == 0:-1
@test search(astr, r"l") == 3:3
@test search(astr, r"l", 4) == 4:4
@test search(astr, r"l", 5) == 11:11
@test search(astr, r"l", 12) == 0:-1
@test search(astr, r"\n") == 14:14
@test search(astr, r"\n", 15) == 0:-1
@test search(u8str, r"z") == 0:-1
@test search(u8str, r"∄") == 0:-1
@test search(u8str, r"∀") == 1:1
@test search(u8str, r"∀", 4) == 0:-1
@test search(u8str, r"∀") == search(u8str, r"\u2200")
@test search(u8str, r"∀", 4) == search(u8str, r"\u2200", 4)
@test search(u8str, r"∃") == 13:13
@test search(u8str, r"∃", 16) == 0:-1
@test search(u8str, r"x") == 26:26
@test search(u8str, r"x", 27) == 43:43
@test search(u8str, r"x", 44) == 0:-1
@test search(u8str, r"ε") == 5:5
@test search(u8str, r"ε", 7) == 54:54
@test search(u8str, r"ε", 56) == 0:-1
for i = 1:endof(astr)
    @test search(astr, r"."s, i) == i:i
end
for i = 1:endof(u8str)
    if isvalid(u8str,i)
        @test search(u8str, r"."s, i) == i:i
    end
end

# string search with a zero-char string
for i = 1:endof(astr)
    @test search(astr, "", i) == i:i-1
end
for i = 1:endof(u8str)
    @test search(u8str, "", i) == i:i-1
end
@test search("", "") == 1:0

# string rsearch with a zero-char string
for i = 1:endof(astr)
    @test rsearch(astr, "", i) == i:i-1
end
for i = 1:endof(u8str)
    @test rsearch(u8str, "", i) == i:i-1
end
@test rsearch("", "") == 1:0

# string search with a zero-char regex
for i = 1:endof(astr)
    @test search(astr, r"", i) == i:i-1
end
for i = 1:endof(u8str)
    # TODO: should regex search fast-forward invalid indices?
    if isvalid(u8str,i)
        @test search(u8str, r""s, i) == i:i-1
    end
end

# string search with a two-char string literal
@test search("foo,bar,baz", "xx") == 0:-1
@test search("foo,bar,baz", "fo") == 1:2
@test search("foo,bar,baz", "fo", 3) == 0:-1
@test search("foo,bar,baz", "oo") == 2:3
@test search("foo,bar,baz", "oo", 4) == 0:-1
@test search("foo,bar,baz", "o,") == 3:4
@test search("foo,bar,baz", "o,", 5) == 0:-1
@test search("foo,bar,baz", ",b") == 4:5
@test search("foo,bar,baz", ",b", 6) == 8:9
@test search("foo,bar,baz", ",b", 10) == 0:-1
@test search("foo,bar,baz", "az") == 10:11
@test search("foo,bar,baz", "az", 12) == 0:-1

# issue #9365
# string search with a two-char UTF-8 (2 byte) string literal
@test search("ééé", "éé") == 1:3
@test search("ééé", "éé", 1) == 1:3
# string search with a two-char UTF-8 (3 byte) string literal
@test search("€€€", "€€") == 1:4
@test search("€€€", "€€", 1) == 1:4
# string search with a two-char UTF-8 (4 byte) string literal
@test search("\U1f596\U1f596\U1f596", "\U1f596\U1f596") == 1:5
@test search("\U1f596\U1f596\U1f596", "\U1f596\U1f596", 1) == 1:5

# string search with a two-char UTF-8 (2 byte) string literal
@test search("éé", "éé") == 1:3
@test search("éé", "éé", 1) == 1:3
# string search with a two-char UTF-8 (3 byte) string literal
@test search("€€", "€€") == 1:4
@test search("€€", "€€", 1) == 1:4
# string search with a two-char UTF-8 (4 byte) string literal
@test search("\U1f596\U1f596", "\U1f596\U1f596") == 1:5
@test search("\U1f596\U1f596", "\U1f596\U1f596", 1) == 1:5

# string rsearch with a two-char UTF-8 (2 byte) string literal
@test rsearch("ééé", "éé") == 3:5
@test rsearch("ééé", "éé", endof("ééé")) == 3:5
# string rsearch with a two-char UTF-8 (3 byte) string literal
@test rsearch("€€€", "€€") == 4:7
@test rsearch("€€€", "€€", endof("€€€")) == 4:7
# string rsearch with a two-char UTF-8 (4 byte) string literal
@test rsearch("\U1f596\U1f596\U1f596", "\U1f596\U1f596") == 5:9
@test rsearch("\U1f596\U1f596\U1f596", "\U1f596\U1f596", endof("\U1f596\U1f596\U1f596")) == 5:9

# string rsearch with a two-char UTF-8 (2 byte) string literal
@test rsearch("éé", "éé") == 1:3        # should really be 1:4!
@test rsearch("éé", "éé", endof("ééé")) == 1:3
# string search with a two-char UTF-8 (3 byte) string literal
@test rsearch("€€", "€€") == 1:4        # should really be 1:6!
@test rsearch("€€", "€€", endof("€€€")) == 1:4
# string search with a two-char UTF-8 (4 byte) string literal
@test rsearch("\U1f596\U1f596", "\U1f596\U1f596") == 1:5        # should really be 1:8!
@test rsearch("\U1f596\U1f596", "\U1f596\U1f596", endof("\U1f596\U1f596\U1f596")) == 1:5

# string rsearch with a two-char string literal
@test rsearch("foo,bar,baz", "xx") == 0:-1
@test rsearch("foo,bar,baz", "fo") == 1:2
@test rsearch("foo,bar,baz", "fo", 1) == 0:-1
@test rsearch("foo,bar,baz", "oo") == 2:3
@test rsearch("foo,bar,baz", "oo", 2) == 0:-1
@test rsearch("foo,bar,baz", "o,") == 3:4
@test rsearch("foo,bar,baz", "o,", 1) == 0:-1
@test rsearch("foo,bar,baz", ",b") == 8:9
@test rsearch("foo,bar,baz", ",b", 6) == 4:5
@test rsearch("foo,bar,baz", ",b", 3) == 0:-1
@test rsearch("foo,bar,baz", "az") == 10:11
@test rsearch("foo,bar,baz", "az", 10) == 0:-1

# array rsearch
@test rsearch(UInt8[1,2,3],UInt8[2,3],3) == 2:3
@test rsearch(UInt8[1,2,3],UInt8[2,3],1) == 0:-1

# string search with a two-char regex
@test search("foo,bar,baz", r"xx") == 0:-1
@test search("foo,bar,baz", r"fo") == 1:2
@test search("foo,bar,baz", r"fo", 3) == 0:-1
@test search("foo,bar,baz", r"oo") == 2:3
@test search("foo,bar,baz", r"oo", 4) == 0:-1
@test search("foo,bar,baz", r"o,") == 3:4
@test search("foo,bar,baz", r"o,", 5) == 0:-1
@test search("foo,bar,baz", r",b") == 4:5
@test search("foo,bar,baz", r",b", 6) == 8:9
@test search("foo,bar,baz", r",b", 10) == 0:-1
@test search("foo,bar,baz", r"az") == 10:11
@test search("foo,bar,baz", r"az", 12) == 0:-1

@test searchindex("foo", 'o') == 2
@test searchindex("foo", 'o', 3) == 3

# string searchindex with a two-char UTF-8 (2 byte) string literal
@test searchindex("ééé", "éé") == 1
@test searchindex("ééé", "éé", 1) == 1
# string searchindex with a two-char UTF-8 (3 byte) string literal
@test searchindex("€€€", "€€") == 1
@test searchindex("€€€", "€€", 1) == 1
# string searchindex with a two-char UTF-8 (4 byte) string literal
@test searchindex("\U1f596\U1f596\U1f596", "\U1f596\U1f596") == 1
@test searchindex("\U1f596\U1f596\U1f596", "\U1f596\U1f596", 1) == 1

# string searchindex with a two-char UTF-8 (2 byte) string literal
@test searchindex("éé", "éé") == 1
@test searchindex("éé", "éé", 1) == 1
# string searchindex with a two-char UTF-8 (3 byte) string literal
@test searchindex("€€", "€€") == 1
@test searchindex("€€", "€€", 1) == 1
# string searchindex with a two-char UTF-8 (4 byte) string literal
@test searchindex("\U1f596\U1f596", "\U1f596\U1f596") == 1
@test searchindex("\U1f596\U1f596", "\U1f596\U1f596", 1) == 1

# string rsearchindex with a two-char UTF-8 (2 byte) string literal
@test rsearchindex("ééé", "éé") == 3
@test rsearchindex("ééé", "éé", endof("ééé")) == 3
# string rsearchindex with a two-char UTF-8 (3 byte) string literal
@test rsearchindex("€€€", "€€") == 4
@test rsearchindex("€€€", "€€", endof("€€€")) == 4
# string rsearchindex with a two-char UTF-8 (4 byte) string literal
@test rsearchindex("\U1f596\U1f596\U1f596", "\U1f596\U1f596") == 5
@test rsearchindex("\U1f596\U1f596\U1f596", "\U1f596\U1f596", endof("\U1f596\U1f596\U1f596")) == 5

# string rsearchindex with a two-char UTF-8 (2 byte) string literal
@test rsearchindex("éé", "éé") == 1
@test rsearchindex("éé", "éé", endof("ééé")) == 1
# string searchindex with a two-char UTF-8 (3 byte) string literal
@test rsearchindex("€€", "€€") == 1
@test rsearchindex("€€", "€€", endof("€€€")) == 1
# string searchindex with a two-char UTF-8 (4 byte) string literal
@test rsearchindex("\U1f596\U1f596", "\U1f596\U1f596") == 1
@test rsearchindex("\U1f596\U1f596", "\U1f596\U1f596", endof("\U1f596\U1f596\U1f596")) == 1

@test_throws ErrorException "ab" ∈ "abc"