text:
## Evaluation

getindex{J<:Jacobi}(op::ConcreteEvaluation{J,Bool},k::Integer) =
    op[k:k][1]

getindex{J<:Jacobi}(op::ConcreteEvaluation{J},k::Integer) =
    op[k:k][1]

function getindex{J<:Jacobi}(op::ConcreteEvaluation{J,Bool},kr::Range)
    @assert op.order <= 2
    sp=op.space
    a=sp.a;b=sp.b
    x=op.x

    if op.order == 0
        jacobip(kr-1,a,b,x?1.0:-1.0)
    elseif op.order == 1&& !x && b==0
        d=domain(op)
        @assert isa(d,Interval)
        Float64[tocanonicalD(d,d.a)*.5*(a+k)*(k-1)*(-1)^k for k=kr]
    elseif op.order == 1
        d=domain(op)
        @assert isa(d,Interval)
        if kr[1]==1 && kr[end] ≥ 2
            0.5*tocanonicalD(d,d.a)*(a+b+kr).*[0.;jacobip(0:kr[end]-2,1+a,1+b,x?1.:-1.)]
        elseif kr[1]==1  # kr[end] ≤ 1
            zeros(eltype(op),length(kr))
        else
            0.5*tocanonicalD(d,d.a)*(a+b+kr).*jacobip(kr-1,1+a,1+b,x?1.:-1.)
        end
    elseif op.order == 2
        @assert !x && b==0
        @assert domain(op)==Interval()
        Float64[-.125*(a+k)*(a+k+1)*(k-2)*(k-1)*(-1)^k for k=kr]
    end
end
function getindex{J<:Jacobi}(op::ConcreteEvaluation{J,Float64},kr::Range)
    @assert op.order == 0
    jacobip(kr-1,op.space.a,op.space.b,tocanonical(domain(op),op.x))
end

## Derivative

Derivative(J::Jacobi,k::Integer)=k==1?ConcreteDerivative(J,1):DerivativeWrapper(TimesOperator(Derivative(Jacobi(J.a+1,J.b+1,J.domain),k-1),ConcreteDerivative(J,1)),k)

rangespace{J<:Jacobi}(D::ConcreteDerivative{J})=Jacobi(D.space.a+D.order,D.space.b+D.order,domain(D))
bandinds{J<:Jacobi}(D::ConcreteDerivative{J})=0,D.order

getindex{J<:Jacobi}(T::ConcreteDerivative{J},k::Integer,j::Integer) =
    j==k+1? (k+1+T.space.a+T.space.b)./complexlength(domain(T)) : zero(eltype(T))

function Derivative{T,DDD<:Interval}(S::WeightedJacobi{T,DDD})
    if S.α>0 && S.β>0 && S.α==S.space.b && S.β==S.space.a
        ConcreteDerivative(S,1)
    else
        jacobiweightDerivative(S)
    end
end

bandinds{T,DDD<:Interval}(D::ConcreteDerivative{WeightedJacobi{T,DDD}})=-1,0
rangespace{T,DDD<:Interval}(D::ConcreteDerivative{WeightedJacobi{T,DDD}})=WeightedJacobi(domainspace(D).α-1,domainspace(D).β-1,domain(D))

getindex{T,DDD<:Interval}(D::ConcreteDerivative{WeightedJacobi{T,DDD}},k::Integer,j::Integer) =
    j==k-1? -4(k-1)./complexlength(domain(D)) : zero(eltype(D))

## Integral

function Integral(J::Jacobi,k::Integer)
    if k > 1
        Q=Integral(J,1)
        IntegralWrapper(TimesOperator(Integral(rangespace(Q),k-1),Q),k)
    elseif J.a > 0 && J.b > 0   # we have a simple definition
        ConcreteIntegral(J,1)
    else   # convert and then integrate
        sp=Jacobi(J.a+1,J.b+1,domain(J))
        C=Conversion(J,sp)
        Q=Integral(sp,1)
        IntegralWrapper(TimesOperator(Q,C),1)
    end
end

rangespace{J<:Jacobi}(D::ConcreteIntegral{J})=Jacobi(D.space.a-D.order,D.space.b-D.order,domain(D))
bandinds{J<:Jacobi}(D::ConcreteIntegral{J})=-D.order,0

function getindex{J<:Jacobi}(T::ConcreteIntegral{J},k::Integer,j::Integer)
    @assert T.order==1
    if k≥2 && j==k-1
        complexlength(domain(T))./(k+T.space.a+T.space.b-2)
    else
        zero(eltype(J))
    end
end

## Volterra Integral operator

Volterra(d::Interval) = Volterra(Legendre(d))
function Volterra(S::Jacobi,order::Integer)
    @assert S.a == S.b == 0.0
    @assert order==1
    ConcreteVolterra(S,order)
end

rangespace{J<:Jacobi}(V::ConcreteVolterra{J})=Jacobi(0.0,-1.0,domain(V))
bandinds{J<:Jacobi}(V::ConcreteVolterra{J})=-1,0

function getindex{J<:Jacobi}(V::ConcreteVolterra{J},k::Integer,j::Integer)
    d=domain(V)
    C = 0.5(d.b-d.a)
    if k≥2
        if j==k-1
            C/(k-1.5)
        elseif j==k
            -C/(k-0.5)
        else
            zero(eltype(V))
        end
    else
        zero(eltype(V))
    end
end

## Conversion
# We can only increment by a or b by one, so the following
# multiplies conversion operators to handle otherwise

function Conversion(L::Jacobi,M::Jacobi)
    if isapproxinteger(L.a-M.a) && isapproxinteger(L.b-M.b)
        dm=domain(M)
        D=typeof(dm)
        if isapprox(M.a,L.a) && isapprox(M.b,L.b)
            ConversionWrapper(eye(L))
        elseif (isapprox(M.b,L.b+1) && isapprox(M.a,L.a)) || (isapprox(M.b,L.b) && isapprox(M.a,L.a+1))
            ConcreteConversion(L,M)
        elseif M.b > L.b+1
            ConversionWrapper(TimesOperator(Conversion(Jacobi(M.a,M.b-1,dm),M),Conversion(L,Jacobi(M.a,M.b-1,dm))))
        else  #if M.a >= L.a+1
            ConversionWrapper(TimesOperator(Conversion(Jacobi(M.a-1,M.b,dm),M),Conversion(L,Jacobi(M.a-1,M.b,dm))))
        end
    elseif L.a ≈ L.b ≈ 0. && M.a ≈ M.b ≈ 0.5
        Conversion(L,Ultraspherical(L),Ultraspherical(M),M)
    elseif L.a ≈ L.b ≈ 0. && M.a ≈ M.b ≈ -0.5
        Conversion(L,Ultraspherical(L),Chebyshev(M),M)
    elseif L.a ≈ L.b ≈ -0.5 && M.a ≈ M.b ≈ 0.5
        Conversion(L,Chebyshev(L),Ultraspherical(M),M)
    else # L.a - M.a ≈ L.b - M.b
        error("Implement for $L ↦ $M")
    end
end

bandinds{J<:Jacobi}(C::ConcreteConversion{J,J})=(0,1)

function Base.getindex{J<:Jacobi,T}(C::ConcreteConversion{J,J,T},k::Integer,j::Integer)
    L=C.domainspace
    if L.b+1==C.rangespace.b
        if j==k
            k==1?T(1):T((L.a+L.b+k)/(L.a+L.b+2k-1))
        elseif j==k+1
            T((L.a+k)./(L.a+L.b+2k+1))
        else
            zero(T)
        end
    elseif L.a+1==C.rangespace.a
        if j==k
            k==1?T(1):T((L.a+L.b+k)/(L.a+L.b+2k-1))
        elseif j==k+1
            T(-(L.b+k)./(L.a+L.b+2k+1))
        else
            zero(T)
        end
    else
        error("Not implemented")
    end
end

# return the space that has banded Conversion to the other
function conversion_rule(A::Jacobi,B::Jacobi)
    if !isapproxinteger(A.a-B.a) || !isapproxinteger(A.b-B.b)
        NoSpace()
    else
        Jacobi(min(A.a,B.a),min(A.b,B.b),domain(A))
    end
end

## Ultraspherical Conversion

# Assume m is compatible

function Conversion(A::PolynomialSpace,B::Jacobi)
    J = Jacobi(A)
    J == B ? ConcreteConversion(A,B) :
             ConversionWrapper(TimesOperator(Conversion(J,B),Conversion(A,J)))
end

function Conversion(A::Jacobi,B::PolynomialSpace)
    J = Jacobi(B)
    J == A ? ConcreteConversion(A,B) :
             ConversionWrapper(TimesOperator(Conversion(J,B),Conversion(A,J)))
end

function Conversion(A::Jacobi,B::Chebyshev)
    if A.a == A.b == -0.5
        ConcreteConversion(A,B)
    elseif A.a == A.b == 0
        ConversionWrapper(
            SpaceOperator(
                Conversion(Ultraspherical(1//2),B),
                A,B))
    elseif A.a == A.b
        US = Ultraspherical(A)
        ConversionWrapper(Conversion(US,B)*Conversion(A,US))
    else
        J = Jacobi(B)
        Conversion(J,B)*Conversion(A,J)
    end
end

function Conversion(A::Chebyshev,B::Jacobi)
    if B.a == B.b == -0.5
        ConcreteConversion(A,B)
    elseif B.a == B.b == 0
        ConversionWrapper(
            SpaceOperator(
                Conversion(A,Ultraspherical(1//2,domain(B))),
                A,B))
    elseif B.a == B.b
        US = Ultraspherical(B)
        ConversionWrapper(Conversion(US,B)*Conversion(A,US))
    else
        J = Jacobi(A)
        Conversion(J,B)*Conversion(A,J)
    end
end

function Conversion(A::Jacobi,B::Ultraspherical)
    if A.a == A.b == -0.5
        ConversionWrapper(Conversion(Chebyshev(domain(A)),B)*Conversion(A,Chebyshev(domain(A))))
    elseif A.a == A.b == order(B)-0.5
        ConcreteConversion(A,B)
    elseif A.a == A.b == 0
        ConversionWrapper(
            SpaceOperator(
                Conversion(Ultraspherical(1//2),B),
                A,B))
    elseif A.a == A.b
        US = Ultraspherical(A)
        ConversionWrapper(Conversion(US,B)*Conversion(A,US))
    else
        J = Jacobi(B)
        Conversion(J,B)*Conversion(A,J)
    end
end

function Conversion(A::Ultraspherical,B::Jacobi)
    if B.a == B.b == -0.5
        ConversionWrapper(Conversion(Chebyshev(domain(A)),B)*Conversion(A,Chebyshev(domain(A))))
    elseif B.a == B.b == order(A)-0.5
        ConcreteConversion(A,B)
    elseif B.a == B.b == 0
        ConversionWrapper(
            SpaceOperator(
                Conversion(A,Ultraspherical(1//2,domain(B))),
                A,B))
    elseif B.a == B.b
        US = Ultraspherical(B)
        ConversionWrapper(Conversion(US,B)*Conversion(A,US))
    else
        J = Jacobi(A)
        Conversion(J,B)*Conversion(A,J)
    end
end

bandinds{US<:Chebyshev,J<:Jacobi}(C::ConcreteConversion{US,J}) = 0,0
bandinds{US<:Chebyshev,J<:Jacobi}(C::ConcreteConversion{J,US}) = 0,0

bandinds{US<:Ultraspherical,J<:Jacobi}(C::ConcreteConversion{US,J}) = 0,0
bandinds{US<:Ultraspherical,J<:Jacobi}(C::ConcreteConversion{J,US}) = 0,0

function getindex{J<:Jacobi,CC<:Chebyshev,T}(C::ConcreteConversion{CC,J,T},k::Integer,j::Integer)
    if j==k
        one(T)/jacobip(k-1,-one(T)/2,-one(T)/2,one(T))
    else
        zero(T)
    end
end

function getindex{J<:Jacobi,CC<:Chebyshev,T}(C::ConcreteConversion{J,CC,T},k::Integer,j::Integer)
    if j==k
        jacobip(k-1,-one(T)/2,-one(T)/2,one(T))
    else
        zero(T)
    end
end

function getindex{US<:Ultraspherical,J<:Jacobi,T}(C::ConcreteConversion{US,J,T},k::Integer,j::Integer)
    if j==k
        S=rangespace(C)
        jp=jacobip(k-1,S.a,S.b,one(T))
        um=Evaluation(setcanonicaldomain(domainspace(C)),one(T))[k]
        um/jp
    else
        zero(T)
    end
end

function getindex{US<:Ultraspherical,J<:Jacobi,T}(C::ConcreteConversion{J,US,T},k::Integer,j::Integer)
    if j==k
        S=domainspace(C)
        jp=jacobip(k-1,S.a,S.b,one(T))
        um=Evaluation(setcanonicaldomain(rangespace(C)),one(T))[k]
        jp/um
    else
        zero(T)
    end
end

union_rule(A::Jacobi,B::Jacobi)=conversion_type(A,B)
maxspace_rule(A::Jacobi,B::Jacobi) = Jacobi(max(A.a,B.a),max(A.b,B.b),domain(A))

for (OPrule,OP) in ((:conversion_rule,:conversion_type),(:maxspace_rule,:maxspace),(:union_rule,:(Base.union)))
    @eval begin
        function $OPrule(A::Chebyshev,B::Jacobi)
            if isapprox(B.a,-0.5)&&isapprox(B.b,-0.5)
                # the spaces are the same
                A
            else
                $OP(Jacobi(A),B)
            end
        end
        function $OPrule(A::Ultraspherical,B::Jacobi)
            m=order(A)
            if isapprox(B.a,m-0.5)&&isapprox(B.b,m-0.5)
                # the spaces are the same
                A
            else
                $OP(Jacobi(A),B)
            end
        end
    end
end

hasconversion(a::Jacobi,b::Chebyshev) = hasconversion(a,Jacobi(b))
hasconversion(a::Chebyshev,b::Jacobi) = hasconversion(Jacobi(a),b)

hasconversion(a::Jacobi,b::Ultraspherical) = hasconversion(a,Jacobi(b))
hasconversion(a::Ultraspherical,b::Jacobi) = hasconversion(Jacobi(a),b)

## Special Multiplication
# special multiplication operators exist when multiplying by
# (1+x) or (1-x) by _decreasing_ the parameter.  Thus the

## <: IntervalDomain avoids a julia bug
function Multiplication{C<:ConstantSpace,DD<:IntervalDomain}(f::Fun{JacobiWeight{C,DD}},S::Jacobi)
    # this implements (1+x)*P and (1-x)*P special case
    # see DLMF (18.9.6)
    d=domain(f)
    if ((space(f).α==1 && space(f).β==0 && S.b >0) ||
                        (space(f).α==0 && space(f).β==1 && S.a >0))
        ConcreteMultiplication(f,S)
    elseif isapproxinteger(space(f).α) && space(f).α ≥ 1 && S.b >0
        # decrement α and multiply again
        M=Multiplication(f.coefficients[1]*jacobiweight(1.,0.,d),S)
        MultiplicationWrapper(f,Multiplication(jacobiweight(space(f).α-1,space(f).β,d),rangespace(M))*M)
    elseif isapproxinteger(space(f).β) && space(f).β ≥ 1 && S.a >0
        # decrement β and multiply again
        M=Multiplication(f.coefficients[1]*jacobiweight(0.,1.,d),S)
        MultiplicationWrapper(f,Multiplication(jacobiweight(space(f).α,space(f).β-1,d),rangespace(M))*M)
    else
# default JacobiWeight
        M=Multiplication(Fun(f.coefficients,space(f).space),S)
        rsp=JacobiWeight(space(f).α,space(f).β,rangespace(M))
        MultiplicationWrapper(f,SpaceOperator(M,S,rsp))
    end
end

Multiplication{C<:ConstantSpace,DD<:IntervalDomain}(f::Fun{JacobiWeight{C,DD}},S::Union{Ultraspherical,Chebyshev}) =
    MultiplicationWrapper(f,Multiplication(f,Jacobi(S))*Conversion(S,Jacobi(S)))

function rangespace{J<:Jacobi,C<:ConstantSpace,DD<:IntervalDomain}(M::ConcreteMultiplication{JacobiWeight{C,DD},J})
    S=domainspace(M)
    if space(M.f).α==1
        # multiply by (1+x)
        Jacobi(S.a,S.b-1,domain(S))
    elseif space(M.f).β == 1
        # multiply by (1-x)
        Jacobi(S.a-1,S.b,domain(S))
    else
        error("Not implemented")
    end
end

bandinds{J<:Jacobi,C<:ConstantSpace,DD<:IntervalDomain}(::ConcreteMultiplication{JacobiWeight{C,DD},J})=-1,0

function getindex{J<:Jacobi,C<:ConstantSpace,DD<:IntervalDomain}(M::ConcreteMultiplication{JacobiWeight{C,DD},J},k::Integer,j::Integer)
    @assert ncoefficients(M.f)==1
    a,b=domainspace(M).a,domainspace(M).b
    c=M.f.coefficients[1]
    if space(M.f).α==1
        @assert space(M.f).β==0
        # multiply by (1+x)
        if j==k
            c*2(k+b-1)/(2k+a+b-1)
        elseif k > 1 && j==k-1
            c*(2k-2)/(2k+a+b-3)
        else
            zero(eltype(M))
        end
    elseif space(M.f).β == 1
        @assert space(M.f).α==0
        # multiply by (1-x)
        if j==k
            c*2(k+a-1)/(2k+a+b-1)
        elseif k > 1 && j==k-1
            -c*(2k-2)/(2k+a+b-3)
        else
            zero(eltype(M))
        end
    else
        error("Not implemented")
    end
end

# We can exploit the special multiplication to construct a Conversion

for FUNC in (:maxspace_rule,:union_rule,:hasconversion)
    @eval function $FUNC{T,DD<:Interval}(A::WeightedJacobi{T,DD},B::Jacobi)
        if A.α==A.β+1 && A.space.b>0
            $FUNC(Jacobi(A.space.a,A.space.b-1,domain(A)),B)
        elseif A.β==A.α+1 && A.space.a>0
            $FUNC(Jacobi(A.space.a-1,A.space.b,domain(A)),B)
        else
            $FUNC(A,JacobiWeight(0.,0.,B))
        end
    end
end

# represents [b+(1+z)*d/dz] (false) or [a-(1-z)*d/dz] (true)
immutable JacobiSD{T} <:Operator{T}
    lr::Bool
    S::Jacobi
end

JacobiSD(lr,S)=JacobiSD{Float64}(lr,S)

Base.convert{T}(::Type{Operator{T}},SD::JacobiSD)=JacobiSD{T}(SD.lr,SD.S)

domain(op::JacobiSD)=domain(op.S)
domainspace(op::JacobiSD)=op.S
rangespace(op::JacobiSD)=op.lr?Jacobi(op.S.a-1,op.S.b+1,domain(op.S)):Jacobi(op.S.a+1,op.S.b-1,domain(op.S))
bandinds(::JacobiSD)=0,0

function getindex(op::JacobiSD,A,k::Integer,j::Integer)
    m=op.lr?op.S.a:op.S.b
    if k==j
        k+m-1
    else
        zero(eltype(op))
    end
end


text:
@testset "867.transpose-matrix.jl" begin
    @test transpose_matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) == [[1, 4, 7], [2, 5, 8], [3, 6, 9]]
    @test transpose_matrix([[1, 2, 3], [4, 5, 6]]) == [[1, 4], [2, 5], [3, 6]]
end

ext: jl


text:
# linear phase FIR filter design that optimizes maximum error
# in the frequency domain

#=============================================

Translated from C code in scipy into Julia.
Tom Krauss, 2018.

SCIPY license:
Copyright (c) 2001, 2002 Enthought, Inc.
All rights reserved.

Copyright (c) 2003-2017 SciPy Developers.
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

  a. Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
  b. Redistributions in binary form must reproduce the above copyright
     notice, this list of conditions and the following disclaimer in the
     documentation and/or other materials provided with the distribution.
  c. Neither the name of Enthought nor the names of the SciPy Developers
     may be used to endorse or promote products derived from this software
     without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS
BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.

--- end of scipy license

The remez and associated functions were extracted from sigtoolsmodule.c,
with the following header comment:

SIGTOOLS module by Travis Oliphant

Copyright 2005 Travis Oliphant
Permission to use, copy, modify, and distribute this software without fee
is granted under the SciPy License.

C CODE BANNER

/********************************************************
 *
 *  Code taken from remez.c by Erik Kvaleberg which was
 *    converted from an original FORTRAN by:
 *
 * AUTHORS: JAMES H. MCCLELLAN
 *
 *         DEPARTMENT OF ELECTRICAL ENGINEERING AND COMPUTER SCIENCE
 *         MASSACHUSETTS INSTITUTE OF TECHNOLOGY
 *         CAMBRIDGE, MASS. 02139
 *
 *         THOMAS W. PARKS
 *         DEPARTMENT OF ELECTRICAL ENGINEERING
 *         RICE UNIVERSITY
 *         HOUSTON, TEXAS 77001
 *
 *         LAWRENCE R. RABINER
 *         BELL LABORATORIES
 *         MURRAY HILL, NEW JERSEY 07974
 *
 *
 *  Adaptation to C by
 *      egil kvaleberg
 *      husebybakken 14a
 *      0379 oslo, norway
 *  Email:
 *      egil@kvaleberg.no
 *  Web:
 *      http://www.kvaleberg.com/
 *
 *********************************************************/

=============================================#

# RemezFilterType:
#    Type I and II symmetric linear phase: neg==0   (filter_type==bandpass)
#    Type III and IV negative symmetric linear phase: neg==1   (filter_type==hilbert or differentiator)
@enum RemezFilterType filter_type_bandpass=1 filter_type_differentiator=2 filter_type_hilbert=3

"""/*
lagrange_interp(k::Integer, n::Integer, m::Integer, x::AbstractVector)

CALCULATE THE LAGRANGE INTERPOLATION COEFFICIENTS
"""
function lagrange_interp(k::Integer, n::Integer, m::Integer, x::AbstractVector)
    retval = 1.0
    q = x[k]
    for l = 1 : m
        for j = l : m : n
            if j != k
                retval *= 2.0 * (q - x[j])
            end
        end
    end
    1.0 / retval
end

"""
build_grid(numtaps, band_defs, Hz, grid_density, neg)
return "grid" and "des" and "wt" arrays
"""
function build_grid(nfilt, band_defs, Hz, grid_density, neg)
    nodd = isodd(nfilt)
    nfcns = nfilt ÷ 2
    if nodd && !neg
        nfcns = nfcns + 1
    end

    #
    # SET UP THE DENSE GRID. THE NUMBER OF POINTS IN THE GRID
    # IS (FILTER LENGTH + 1)*GRID DENSITY/2
    #
    delf = 0.5 / (grid_density * nfcns)

    flimlow = neg ? delf : 0.0
    flimhigh = (neg == nodd) ? 0.5 - delf : 0.5
    function normalize_banddef_entry(b)
        # normalize and clamp band-edges
        fl = convert(Float64, clamp(b.first[1] / Hz, flimlow, flimhigh))
        fu = convert(Float64, clamp(b.first[2] / Hz, flimlow, flimhigh))
        # make sure desired and weight are functions
        if !isa(b.second, Tuple{Any, Any})
            desired = b.second
            weight = 1.0
        else
            desired = b.second[1]
            weight = b.second[2]
        end
        if isa(desired, Real)
            let d = desired
                desired = _ -> d
            end
        end
        if isa(weight, Real)
            let w = weight
                weight = _ -> w
            end
        end
        return Pair{Tuple{Float64,Float64},Tuple{Any,Any}}((fl, fu), (desired, weight))
    end
    normalized_band_defs = normalize_banddef_entry.(band_defs)

    local ngrid
    # work around JuliaLang/julia#15276
    let delf=delf
        ngrid = sum(max(length(band_def.first[1]:delf:band_def.first[2]), 1) for band_def in normalized_band_defs)#::Int
    end

    grid = zeros(Float64, ngrid)  # the array of frequencies, between 0 and 0.5
    des = zeros(Float64, ngrid)   # the desired function on the grid
    wt = zeros(Float64, ngrid)    # array of weights

    j = 1

    #
    # CALCULATE THE DESIRED MAGNITUDE RESPONSE AND THE WEIGHT
    # FUNCTION ON THE GRID
    #
    # and
    #
    # SET UP A NEW APPROXIMATION PROBLEM WHICH IS EQUIVALENT
    # TO THE ORIGINAL PROBLEM
    #
    for band_def in normalized_band_defs
        flow = band_def.first[1]
        fup = band_def.first[2]
        # outline inner loop to have it type-stable (band_def.second is Tuple{Any,Any})
        j = _buildgrid!(grid, des, wt, j, [(flow:delf:fup)[1:end-1]; fup], neg, nodd, Hz,
                        band_def.second)
    end
    @assert ngrid == j - 1

    return grid, des, wt
end

function _buildgrid!(grid, des, wt, j, fs, neg, nodd, Hz, des_wt)
    for f in fs
        change = neg ? (nodd ? sinpi(2f) : sinpi(f)) : (nodd ? 1.0 : cospi(f))
        grid[j] = cospi(2f)
        des[j] = des_wt[1](f * Hz) / change
        wt[j] = des_wt[2](f * Hz) * change
        j += 1
    end
    return j
end

"""
function freq_eval(xf, x::AbstractVector, y::AbstractVector, ad::AbstractVector)
-----------------------------------------------------------------------
 FUNCTION: freq_eval (gee)
  FUNCTION TO EVALUATE THE FREQUENCY RESPONSE USING THE
  LAGRANGE INTERPOLATION FORMULA IN THE BARYCENTRIC FORM
-----------------------------------------------------------------------
"""
function freq_eval(xf, x::AbstractVector, y::AbstractVector, ad::AbstractVector)
    d = 0.0
    p = 0.0

    for j = 1:length(ad)
        c = ad[j] / (xf - x[j])
        d += c
        p += c * y[j]
    end

    p/d
end

function initialize_y(dev::Float64, nz::Integer, iext::AbstractArray, des::AbstractArray, wt::AbstractArray, y::AbstractArray)
    nu = 1
    if dev > 0.0
        nu = -1
    end
    dev = -nu * dev
    k = nu
    for j = 1 : nz
        l = iext[j]
        y[j] = des[l] + k * dev / wt[l]
        k = -k
    end
    nu, dev
end

#=========
Banner from C code

-----------------------------------------------------------------------
 SUBROUTINE: remez
  THIS SUBROUTINE IMPLEMENTS THE REMEZ EXCHANGE ALGORITHM
  FOR THE WEIGHTED CHEBYSHEV APPROXIMATION OF A CONTINUOUS
  FUNCTION WITH A SUM OF COSINES.  INPUTS TO THE SUBROUTINE
  ARE A DENSE GRID WHICH REPLACES THE FREQUENCY AXIS, THE
  DESIRED FUNCTION ON THIS GRID, THE WEIGHT FUNCTION ON THE
  GRID, THE NUMBER OF COSINES, AND AN INITIAL GUESS OF THE
  EXTREMAL FREQUENCIES.  THE PROGRAM MINIMIZES THE CHEBYSHEV
  ERROR BY DETERMINING THE BSMINEST LOCATION OF THE EXTREMAL
  FREQUENCIES (POINTS OF MAXIMUM ERROR) AND THEN CALCULATES
  THE COEFFICIENTS OF THE BEST APPROXIMATION.
-----------------------------------------------------------------------
=========#
"""
    remez(numtaps::Integer, band_defs;
          Hz::Real=1.0,
          neg::Bool=false,
          maxiter::Integer=25,
          grid_density::Integer=16)

Calculate the minimax optimal filter using the Remez exchange algorithm [^McClellan1973a] [^McClellan1973b].

This is the simplified API that accepts just 2 required arguments (numtaps, band_defs).
For a scipy compatible version see the 3 arguments version (numtaps, bands, desired).

Calculate the filter-coefficients for the finite impulse response
(FIR) filter whose transfer function minimizes the maximum error
between the desired gain and the realized gain in the specified
frequency bands using the Remez exchange algorithm.

# Arguments
- `numtaps::Integer`: The desired number of taps in the filter.
    The number of taps is the number of terms in the filter, or the filter
    order plus one.
- `bands_defs`: A sequence of band definitions.
    This sequence defines the bands. Each entry is a pair. The pair's
    first item is a tuple of band edges (low, high). The pair's second item
    defines the desired response and weight in that band. The weight is optional
    and defaults to 1.0. Both the desired response and weight may be either scalars
    or functions. If a function, the function should accept a real frequency and
    return the real desired response or real weight. Examples:
    + LPF with unity weights. `[(0, 0.475) => 1, (0.5, 1.0) => 0]`
    + LPF with weight of 2 in the stop band. `[(0, 0.475) => (1, 1), (0.5, 1.0) => (0, 2)]`
    + BPF with unity weights. `[(0, 0.375) => 0, (0.4, 0.5) => 1, (0.525, 1.0) => 0]`
    + Hilbert transformer. `[(0.1, 0.95) => 1]; neg=true`
    + Differentiator. `[(0.01, 0.99) => (f -> f/2, f -> 1/f)]; neg=true`
- `Hz::Real`: The sampling frequency in Hz. Default is 1.
- `neg::Bool`: Whether the filter has negative symmetry or not. Default is false.
    If false, the filter is even-symmetric. If true, the filter is odd-symmetric.
    neg=true means that h[n]=-h[end+1-n]; neg=false means that h[n]=h[end+1-n].
- `maxiter::Integer`: (optional)
    Maximum number of iterations of the algorithm. Default is 25.
- `grid_density:Integer`: (optional)
    Grid density. The dense grid used in `remez` is of size
    `(numtaps + 1) * grid_density`. Default is 16.

# Returns
- `h::Array{Float64,1}`: A rank-1 array containing the coefficients of the optimal
    (in a minimax sense) filter.

[^McClellan1973a]:
J. H. McClellan and T. W. Parks, A unified approach to the
design of optimum FIR linear phase digital filters,
IEEE Trans. Circuit Theory, vol. CT-20, pp. 697-701, 1973.

[^McClellan1973b]:
J. H. McClellan, T. W. Parks and L. R. Rabiner, A Computer
Program for Designing Optimum FIR Linear Phase Digital
Filters, IEEE Trans. Audio Electroacoust., vol. AU-21,
pp. 506-525, 1973.

# Examples
Construct a length 35 filter with a passband at 0.15-0.4 Hz
(desired response of 1), and stop bands at 0-0.1 Hz and 0.45-0.5 Hz
(desired response of 0). Note: the behavior in the frequency ranges between
those bands - the transition bands - is unspecified.

```jldoctest; setup = :(using DSP)
julia> bpass = remez(35, [(0, 0.1)=>0, (0.15, 0.4)=>1, (0.45, 0.5)=>0]);
```

You can trade-off maximum error achieved for transition bandwidth.
The wider the transition bands, the lower the maximum error in the
bands specified. Here is a bandpass filter with the same passband, but
wider transition bands.

```jldoctest; setup = :(using DSP)
julia> bpass2 = remez(35, [(0, 0.08)=>0, (0.15, 0.4)=>1, (0.47, 0.5)=>0]);
```

Here we compute the frequency responses and plot them in dB.

```julia-repl
julia> using PyPlot
julia> b = DSP.Filters.PolynomialRatio(bpass, [1.0])
julia> b2 = DSP.Filters.PolynomialRatio(bpass2, [1.0])
julia> f = range(0, stop=0.5, length=1000)
julia> plot(f, 20*log10.(abs.(freqresp(b,f,1.0))))
julia> plot(f, 20*log10.(abs.(freqresp(b2,f,1.0))))
julia> grid()
```

# Examples from the unittests - standard (even) symmetry.

Length 151 LPF (Low Pass Filter).
```jldoctest; setup = :(using DSP)
julia> h = remez(151, [(0, 0.475) => 1, (0.5, 1.0) => 0]; Hz=2.0);
```

Length 152 LPF. Non-default "weight" input.
```jldoctest; setup = :(using DSP)
julia> h = remez(152, [(0, 0.475) => (1, 1), (0.5, 1.0) => (0, 2)]; Hz=2.0);
```

Length 51 HPF (High Pass Filter).
```jldoctest; setup = :(using DSP)
julia> h = remez(51, [(0, 0.75) => 0, (0.8, 1.0) => 1]; Hz=2.0);
```

Length 180 BPF (Band Pass Filter).
```jldoctest; setup = :(using DSP)
julia> h = remez(180, [(0, 0.375) => 0, (0.4, 0.5) => 1, (0.525, 1.0) => 0]; Hz=2.0, maxiter=30);
```

# Examples from the unittests - Odd-symmetric filters - hilbert and differentiators type.
Even length - has a much better approximation since the response is not constrained to 0 at
the nyquist frequency.  Length 20 Hilbert transformer.
```jldoctest; setup = :(using DSP)
julia> h = remez(20, [(0.1, 0.95) => 1]; neg=true, Hz=2.0);
```

Length 21 Hilbert transformer.
```jldoctest; setup = :(using DSP)
julia> h = remez(21, [(0.1, 0.95) => 1]; neg=true, Hz=2.0);
```

Length 200 differentiator.
```jldoctest; setup = :(using DSP)
julia> h = remez(200, [(0.01, 0.99) => (f -> f/2, f -> 1/f)]; neg=true, Hz=2.0);
```

Length 201 differentiator.
```jldoctest; setup = :(using DSP)
julia> h = remez(201, [(0.05, 0.95) => (f -> f/2, f -> 1/f)]; neg=true, Hz=2.0);
```

Inverse sinc filter - custom response function
```julia-repl
julia> L = 64; Fs = 4800*L;
julia> passband_response_function = f -> (f==0) ? 1.0 : abs.((π*f/4800) ./ sin.(π*f/4800));
julia> h = remez(201, [(    0.0, 2880.0) => (passband_response_function, 1.0),
                (10000.0,   Fs/2) => (0.0, 100.0)]; Hz=Fs);
```
"""
function remez(numtaps::Integer, band_defs;
                Hz::Real=1.0,
                neg::Bool=false,
                maxiter::Integer=25,
                grid_density::Integer=16)
    all(b.first[1] <= b.first[2] for b in band_defs) ||
        throw(ArgumentError("lower band edge higher then upper band edge"))
    all(band_defs[i].first[2] < band_defs[i+1].first[1] for i in 1:length(band_defs)-1) ||
        throw(ArgumentError("band edges is not monotonically increasing"))
    (0 <= band_defs[1].first[1]) && (band_defs[end].first[2] <= 0.5*Hz) ||
        throw(ArgumentError("band edges must be between 0 and `Hz`/2"))

    grid, des, wt = build_grid(numtaps, band_defs, Hz, grid_density, neg)

    nfilt = numtaps
    ngrid = length(grid)

    nodd = isodd(nfilt)   # boolean: "nodd" means filter length is odd
    nfcns = numtaps ÷ 2   # integer divide
    if nodd && !neg
        nfcns = nfcns + 1
    end

    nz  = nfcns+1
    nzz = nfcns+2
    iext = zeros(Int64, nzz)   # indices of extremals
    x = zeros(Float64, nzz)
    y = zeros(Float64, nz)

    for j = 1:nz
        iext[j] = (j-1)*(ngrid-1) ÷ nfcns + 1
    end

    dev = 0.0     # deviation from the desired function,
                  # that is, the amount of "ripple" on the extremal set
    devl = -1.0   # deviation on last iteration
    niter = 0
    ad = zeros(Float64, nz)

    jet = ((nfcns-1) ÷ 15) + 1

    while true

        #
        # Start next iteration
        #
      @label L100
        iext[nzz] = ngrid + 1
        niter += 1
        if niter > maxiter
            @warn("remez() iteration count exceeds maxiter = $maxiter, filter is not converged; try increasing maxiter")
            # the filter is returned in its current, unconverged state.
            break
        end

        x[1:nz] = grid[iext[1:nz]]

        for j = 1 : nz
            ad[j] = lagrange_interp(j, nz, jet, x)
        end

        dnum = 0.0
        dden = 0.0
        k = 1
        for j = 1 : nz
            l = iext[j]
            dnum += ad[j] * des[l]
            dden += k * ad[j] / wt[l]
            k = -k
        end
        dev = dnum / dden

        fill!(y, 0.0)
        nu, dev = initialize_y(dev, nz, iext, des, wt, y)

        if dev <= devl
            # finished
            throw(ErrorException("remez() - failure to converge at iteration $niter, try reducing transition band width"))
        end
        devl = dev

        #
        # SEARCH FOR THE EXTREMAL FREQUENCIES OF THE BEST APPROXIMATION
        #

        # Between here and L370, the extremal index set is updated in a loop
        # roughly over the index "j" - although the logic is complicated as
        # the extremal set may grow or shrink in an iteration.
        # j - the index of the current extremal being updated
        # nz - the number of cosines in the approximation (including the constant term).
        #      nz = nfcns + 1 where nfcns = nfilt / 2, and
        #      nfilt is the filter length or number of taps.
        #      For example, for a length 15 filter, nfcns = 7 and nz = 8.
        # jchgne - the number of extremal indices that changed this iteration
        jchnge = 0
        k1 = iext[1]
        knz = iext[nz]
        klow = 0
        nut = -nu
        j = 1

      @label L200
        j == nzz && (ynz = comp)   # equivalent to "if (j == nzz) ynz = comp; end"
        j >= nzz && @goto L300
        kup = iext[j+1]
        l = iext[j]+1
        nut = -nut
        j == 2 && (y1 = comp)
        comp = dev
        l >= kup && @goto L220
        err = (freq_eval(grid[l], x, y, ad) - des[l]) * wt[l]
        nut*err <= comp && @goto L220
        comp = nut * err
      @label L210
        l += 1; l >= kup && @goto L215
        err = (freq_eval(grid[l], x, y, ad) - des[l]) * wt[l]
        nut*err <= comp && @goto L215
        comp = nut * err
        @goto L210

      @label L215
        iext[j] = l - 1; j += 1
        klow = l - 1
        jchnge += 1
        @goto L200

      @label L220
        l -= 1
      @label L225
        l -= 1; l <= klow && @goto L250
        err = (freq_eval(grid[l], x, y, ad) - des[l]) * wt[l]
        nut*err > comp && @goto L230
        jchnge <= 0 && @goto L225
        @goto L260

      @label L230
        comp = nut * err
      @label L235
        l -= 1; l <= klow && @goto L240
        err = (freq_eval(grid[l], x, y, ad) - des[l]) * wt[l]
        nut*err <= comp && @goto L240
        comp = nut * err
        @goto L235
      @label L240
        klow = iext[j]
        iext[j] = l+1
        j += 1
        jchnge += 1
        @goto L200

      @label L250
        l = iext[j]+1
        jchnge > 0 && @goto L215

      @label L255
        l += 1; l >= kup && @goto L260
        err = (freq_eval(grid[l], x, y, ad) - des[l]) * wt[l]
        nut*err <= comp && @goto L255
        comp = nut * err

        @goto L210
      @label L260
        klow = iext[j]; j += 1
        @goto L200

      @label L300
        j > nzz && @goto L320
        k1 > iext[1] && (k1 = iext[1])
        knz < iext[nz] && (knz = iext[nz])
        nut1 = nut
        nut = -nu
        l = 0
        kup = k1
        comp = ynz*(1.00001)
        luck = 1
      @label L310
        l += 1; l >= kup && @goto L315
        err = (freq_eval(grid[l], x, y, ad) - des[l]) * wt[l]
        nut*err <= comp && @goto L310
        comp =  nut * err
        j = nzz
        @goto L210

      @label L315
        luck = 6
        @goto L325

      @label L320
        luck > 9 && @goto L350
        comp > y1 && (y1 = comp)
        k1 = iext[nzz]
      @label L325
        l = ngrid+1
        klow = knz
        nut = -nut1
        comp = y1*(1.00001)
      @label L330
        l -= 1; l <= klow && @goto L340
        err = (freq_eval(grid[l], x, y, ad) - des[l]) * wt[l]
        nut*err <= comp && @goto L330
        j = nzz
        comp =  nut * err
        luck = luck + 10
        @goto L235
      @label L340
        luck == 6 && @goto L370
        for j = 1 : nfcns
            iext[nzz-j] = iext[nz-j]
        end
        iext[1] = k1
        @goto L100
      @label L350
        for j = 1:nz
            iext[j] = iext[j+1]
        end

        @goto L100
      @label L370

        if jchnge <= 0  # we are done if none of the extremal indices changed
            break
        end
    end  # while

    #
    #    CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION
    #    USING THE INVERSE DISCRETE FOURIER TRANSFORM
    #

    a = zeros(Float64, nfcns)   # frequency response on evenly spaced grid
    p = zeros(Float64, nfcns)
    q = zeros(Float64, nfcns-2)
    alpha = zeros(Float64, nzz)   # return vector

    fsh = 1.0e-06
    x[nzz] = -2.0
    delf = 1 / (2*nfcns - 1)
    l = 1

    # Boolean for "kkk" in C code.
    full_grid = (band_defs[1].first[1] == 0.0 && band_defs[end].first[2] == 0.5*Hz) || (nfcns <= 3)
    if !full_grid
        aa    = 2.0/(grid[1]-grid[ngrid])
        bb    = -(grid[1]+grid[ngrid])/(grid[1]-grid[ngrid])
    end

    # Fill in "a" array with the frequency response on an evenly
    # spaced set of frequencies. Care is taken to use "y[l]" -
    # an already computed response on one of the extremals "l" -
    # if the extremal is equal to the frequency ft. If no y[l]
    # matches, a[j] is computed using freq_eval.
    for j = 1 : nfcns
        xt = cospi(2 * (j - 1) * delf)
        if !full_grid
            xt = (xt-bb)/aa
        end
        if (l > 1)
            l = l-1
        end
        while x[l]-xt >= fsh
            l += 1
        end
        if xt-x[l] < fsh
            a[j] = y[l]
        else
            a[j] = freq_eval(xt, x, y, ad)
        end
    end

    nm1 = nfcns - 1   # nm1 => "nfcns minus 1"

    for j = 1 : nfcns
        dtemp = 0.0
        for k = 1 : nm1
            dtemp += a[k+1] * cospi(2 * (j-1) * delf * k)
        end
        alpha[j] = 2.0 * dtemp + a[1]
    end

    for j = 2 : nfcns
        alpha[j] *= 2.0 * delf
    end
    alpha[1] *= delf

    if !full_grid
        p[1] = 2.0*alpha[nfcns]*bb+alpha[nm1]
        p[2] = 2.0*aa*alpha[nfcns]
        q[1] = alpha[nfcns-2]-alpha[nfcns]
        for j = 2 : nm1
            if j >= nm1
                aa *= 0.5
                bb *= 0.5
            end
            p[j+1] = 0.0
            for k = 1 : j
                a[k] = p[k]
                p[k] = 2.0 * bb * a[k]
            end
            p[2] += a[1] * 2.0 *aa
            for k = 1 : j-1
                p[k] += q[k] + aa * a[k+1]
            end
            for k = 3 : j+1
                p[k] += aa * a[k-1]
            end

            if j != nm1
                for k = 1 : j
                    q[k] = -a[k]
                end
                q[1] += alpha[nfcns - 1 - j]
            end
        end
        for j = 1 : nfcns
            alpha[j] = p[j]
        end
    end

    if nfcns <= 3
        alpha[nfcns+1] = alpha[nfcns+2] = 0.0
    end

    #
    # CALCULATE THE IMPULSE RESPONSE.
    #
    h = zeros(Float64, nfilt)
    if !neg
        if nodd
            for j = 1 : nm1
                h[j] = 0.5 * alpha[nz-j]
            end
            h[nfcns] = alpha[1]
        else
            h[1] = 0.25 * alpha[nfcns]
            for j = 2 : nm1
                h[j] = 0.25 * (alpha[nz-j] + alpha[nfcns+2-j])
            end
            h[nfcns] = 0.5*alpha[1] + 0.25*alpha[2]
        end
    else
        if nodd
            h[1] = 0.25 * alpha[nfcns]
            h[2] = 0.25 * alpha[nm1]
            for j = 1 : nm1
                h[j] = 0.25 * (alpha[nz-j] - alpha[nfcns+3-j])
            end
            h[nfcns] = 0.5 * alpha[1] - 0.25 * alpha[3]
            h[nz] = 0.0
        else
            h[1] = 0.25 * alpha[nfcns]
            for j = 2 : nm1
                h[j] = 0.25 * (alpha[nz-j] - alpha[nfcns+2-j])
            end
            h[nfcns] = 0.5 * alpha[1] - 0.25 * alpha[2]
        end
    end

    for j = 1 : nfcns
        k = nfilt + 1 - j
        if !neg
           h[k] = h[j]
        else
           h[k] = -h[j]
        end
    end
    if neg && nodd
        h[nz] = 0.0
    end

    return h
end

"""
    remez(numtaps::Integer,
          bands::Vector,
          desired::Vector;
          weight::Vector=[],
          Hz::Real=1.0,
          filter_type::RemezFilterType=filter_type_bandpass,
          maxiter::Integer=25,
          grid_density::Integer=16)

This is the scipy compatible version that requires 3 arguments (numtaps, bands, desired).
For a simplified API, see the 2 argument version (numtaps, band_defs). The filters
designed are equivalent, the inputs are just specified in a different way.
Below the arguments and examples are described that differ from the simplified
API version.

# Arguments
- `bands::Vector`: A monotonic sequence containing the band edges in Hz.
    All elements must be non-negative and less than half the sampling
    frequency as given by `Hz`.
- `desired::Vector`:A sequence half the size of bands containing the desired
    gain in each of the specified bands.
- `weight::Vector`: (optional)
    A relative weighting to give to each band region. The length of
    `weight` has to be half the length of `bands`.
- `filter_type::RemezFilterType`: Default is `filter_type_bandpass`.
    The type of filter:
    +  `filter_type_bandpass` : flat response in bands. This is the default.
    +  `filter_type_differentiator` : frequency proportional response in bands.
        Odd symetric as in `filter_type_hilbert` case, but with a linear sloping
        desired response.
    +  `filter_type_hilbert` : filter with odd symmetry, that is, type III
                  (for even order) or type IV (for odd order)
                  linear phase filters.

# Examples
Compare the examples with the simplified API and the Scipy API.
Each of the following blocks first designs a filter using the
simplified (recommended) API, and then designs the same filter
using the Scipy-compatible API.

```jldoctest; setup = :(using DSP)
julia> bpass = remez(35, [(0, 0.1)=>0, (0.15, 0.4)=>1, (0.45, 0.5)=>0]);

julia> bpass = remez(35, [0, 0.1, 0.15, 0.4, 0.45, 0.5], [0, 1, 0]);

```

```jldoctest; setup = :(using DSP)
julia> bpass2 = remez(35, [(0, 0.08)=>0, (0.15, 0.4)=>1, (0.47, 0.5)=>0]);

julia> bpass2 = remez(35, [0, 0.08, 0.15, 0.4, 0.47, 0.5], [0, 1, 0]);

```

```jldoctest; setup = :(using DSP)
julia> h = remez(20, [(0.1, 0.95) => 1]; neg=true, Hz=2.0);

julia> h = remez(20, [0.1, 0.95], [1]; filter_type=filter_type_hilbert, Hz=2.0);

```

```jldoctest; setup = :(using DSP)
julia> h = remez(200, [(0.01, 0.99) => (f -> f/2, f -> 1/f)]; neg=true, Hz=2.0);

julia> h = remez(200, [0.01, 0.99], [1]; filter_type=filter_type_differentiator, Hz=2.0);

```
"""
function remez(numtaps::Integer, bands::Vector, desired::Vector;
               weight::Vector=fill(1.0, length(desired)),
               Hz::Real=1.0,
               filter_type::RemezFilterType=filter_type_bandpass,
               kwargs...)
    issorted(bands) || throw(ArgumentError("`bands` is not monotonically increasing"))
    length(bands) == 2length(desired) ||
        throw(ArgumentError("`desired` must be half the length of `bands`."))
    length(bands) == 2length(weight) ||
        throw(ArgumentError("`weight` must be half the length of `bands`."))
    band_ranges = [(bands[i], bands[i+1]) for i in 1:2:length(bands)]
    if filter_type == filter_type_differentiator
        eff = [f -> d*f/Hz for d in desired]
        wate = [d > 0.0001 ? f -> w/(f/Hz) : f -> w for (w, d) in zip(weight, desired)]
    else
        eff = desired
        wate = weight
    end
    band_defs = [r => (d, w) for (r, d, w) in zip(band_ranges, eff, wate)]
    neg = filter_type in (filter_type_hilbert, filter_type_differentiator)
    return remez(numtaps, band_defs; Hz=Hz, neg=neg, kwargs...)
end


text:
"""
plot3d(arg1::Array; kwargs...)

reads (x,y,z) triplets from files [or standard input] and generates PostScript code that will plot lines,
polygons, or symbols at those locations in 3-D.

Full option list at [`psxyz`](http://gmt.soest.hawaii.edu/doc/latest/psxyz.html)

Parameters
----------

- **A** : **straight_lines** : -- Str --

    By default, geographic line segments are drawn as great circle arcs. To draw them as straight lines, use the -A flag.
- $(GMT.opt_J)
- $(GMT.opt_Jz)
- $(GMT.opt_R)
- $(GMT.opt_B)
- **C** : **color** : -- Str --

    Give a CPT or specify -Ccolor1,color2[,color3,...] to build a linear continuous CPT from those colors automatically.
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/psxyz.html#c)
- **D** : **offset** : -- Str --

    Offset the plot symbol or line locations by the given amounts dx/dy.
- **E** : **error_bars** : -- Str --

    Draw symmetrical error bars.
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/psxyz.html#e)
- **F** : **conn** : **connection** : -- Str --

    Alter the way points are connected
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/psxyz.html#f)
- **G** : **fill** : **markerfacecolor** : **MarkerFaceColor** : -- Str --

    Select color or pattern for filling of symbols or polygons. BUT WARN: the alias 'fill' will set the
    color of polygons OR symbols but not the two together. If your plot has polygons and symbols, use
    'fill' for the polygons and 'markerfacecolor' for filling the symbols. Same applyies for W bellow
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/psxyz.html#g)
- **I** : **intens** : -- Str or number --

    Use the supplied intens value (in the [-1 1] range) to modulate the fill color by simulating illumination.
- **L** : **closed_polygon** : -- Str --

    Force closed polygons.
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/psxyz.html#l)
- **N** : **no_clip** : --- Str or [] --

    Do NOT clip symbols that fall outside map border
- $(GMT.opt_P)
- **S** : **symbol** : **marker** : **Marker** : -- Str --

    Plot symbols (including vectors, pie slices, fronts, decorated or quoted lines).
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/psxyz.html#s)
    Alternatively select a sub-set of symbols using the aliases: **marker** or **Marker** and values:

    + **-**, **x_dash**
    + **+**, **plus**
    + **a**, *, **star**
    + **c**, **circle**
    + **d**, **diamond**
    + **g**, **octagon**
    + **h**, **hexagon**
    + **i**, **v**, **inverted_tri**
    + **n**, **pentagon**
    + **p**, **.**, **point**
    + **r**, **rectangle**
    + **s**, **square**
    + **t**, **^**, **triangle**
    + **x**, **cross**
    + **y**, **y_dash**
- **W** : **line_attribs** : **markeredgecolor** : **MarkerEdgeColor** : -- Str --
    Set pen attributes for lines or the outline of symbols
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/psxyz.html#w)
    WARNING: the pen attributes will set the pen of polygons OR symbols but not the two together.
    If your plot has polygons and symbols, use **W** or **line_attribs** for the polygons and
    **markeredgecolor** or **MarkerEdgeColor** for filling the symbols. Similar to S above.
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_a)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)
"""
plot3d(arg1::Array; extra="", K=false, O=false, first=true, kw...) =
    GMT.xyz(extra, arg1; caller="plot3d", K=K, O=O, first=first, kw...)
plot3d!(arg1::Array; extra="", K=false, O=false, first=true, kw...) =
    GMT.xyz(extra, arg1; caller="plot3d", K=true, O=true, first=false, kw...)

# -----------------------------------------------------------------------------------------------------
plot3d(arg1::GMTdataset; extra="", K=false, O=false, first=true, kw...) =
	GMT.xyz(extra, arg1; caller="plot3d", K=K, O=O, first=first, kw...)
plot3d!(arg1::GMTdataset; extra="", K=false, O=false, first=true, kw...) =
	GMT.xyz(extra, arg1; caller="plot3d", K=true, O=true, first=false, kw...)
plot3d(arg1::Array{GMTdataset,1}; extra="", K=false, O=false, first=true, kw...) =
	GMT.xyz(extra, arg1; caller="plot3d", K=K, O=O, first=first, kw...)
plot3d!(arg1::Array{GMTdataset,1}; extra="", K=false, O=false, first=true, kw...) =
	GMT.xyz(extra, arg1; caller="plot3d", K=true, O=true, first=false, kw...)

# ------------------------------------------------------------------------------------------------------
plot3d(arg1::String; extra="", K=false, O=false, first=true, kw...) =
	GMT.xyz(extra, []; caller="plot3d", K=K, O=O, first=first, kw...)
plot3d!(arg1::String; extra="", K=false, O=false, first=true, kw...) =
	GMT.xyz(extra, []; caller="plot3d", K=true, O=true, first=false, kw...)

# ------------------------------------------------------------------------------------------------------
function plot3d(arg1::Array, arg2::Array; extra="", K=false, O=false, first=true, kw...)
	arg = hcat(arg1, arg2)
	GMT.xyz(extra, arg; caller="plot3d", K=K, O=O, first=first, kw...)
end
function plot3d!(arg1::Array, arg2::Array; extra="", K=false, O=false, first=true, kw...)
	arg = hcat(arg1, arg2)
	GMT.xyz(extra, arg; caller="plot3d", K=true, O=true, first=false, kw...)
end
# ------------------------------------------------------------------------------------------------------


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

using Base.Test
using TimerOutputs
const to = TimerOutput()

test_files = String[]
push!(test_files, "test_create_basis.jl")
push!(test_files, "test_lagrange_elements.jl")
push!(test_files, "test_nurbs_elements.jl")
push!(test_files, "test_math.jl")

@testset "FEMBasis.jl" begin
    for fn in test_files
        timeit(to, fn) do
            include(fn)
        end
    end
end
println()
println("Test statistics:")
println(to)


text:
using QRupdate
using LinearOperators

macro fprintf(args...)
    :( f()=@printf $(args...); f() )
end

type OptimResults
    # Statistics
    itns::Int
    time::Float64
    stat::Symbol
    exitmsg::AbstractString
    nprodA::Int
    nprodAt::Int
    numtrim::Int

    # State
    active::Vector{Int}
    state::Vector{Int}
    y::Vector{Float64}
    S::Matrix{Float64}
    R::Matrix{Float64}
end

# Constructor
function OptimResults(m::Int, n::Int)
    active = zeros(Int,0)
    state = zeros(Int,n)
    y = zeros(m)
    S = zeros(m,0)
    R = zeros(0,0)
    OptimResults(0,0,:exit_unknown,"Exit Unknown",0,0,0,
                 active, state, y, S, R)
end

"""
Solve the optimization problem

   DP:       minimize_y   - b'y  +  1/2 lambda y'y
             subject to   bl <= A'y <= bu
using given A,b,bl,bu,lambda.  When bl = -e, bu = e = ones(n,1),
DP is the dual Basis Pursuit problem

   BPdual:   maximize_y   b'y  -  1/2 lambda y'y
             subject to   ||A'y||_inf  <=  1.

General call to BPdual:
(active,state,x,y,S,R,info) = BPdual(A,b,bl,bu,lambda,active,state,y,S,R)

EXAMPLE
=======
Generate data:
[A,b,lambda] = gendata(1);

Solve BPdual:
[active,state,x,y,S,info] = BPdual(A,b,-1,1,lambda);

Resolve with a different lambda
lambda = lambda / 10;
[active,state,x,y,S,R,info] = BPdual(A,b,-1,1,lambda,active,state,y,S,R);

INPUT
=====
A          is an m-by-n linear operator.
b          is an m-vector.
lambda     is a nonnegative scalar.
bl, bu     are n-vectors.
active,state,y,S,R  may be empty (i.e., equal to [])
           or output from BPdual with some previous value of lambda.

OUTPUT
======
active     is an nact-vector of indices j with nact = length(active),
           listing which of the constraints bl(j) <= A(:,j)'y <= bu(j)
           is active.
state      is an n-vector giving the state of each dual constraint.
           state(j) = -2  if         Aj'*y < bl(j)  infeas below
                    = -1  if         Aj'*y = bl(j)  active below
                    = +1  if         Aj'*y = bu(j)  active above
                    = +2  if         Aj'*y > bu(j)  infeas above
                    =  0  if bl(j) < Aj'*y < bu(j)  inactive
x          is an nact-vector of primal solution values.
           The n-vector xx = zeros(n,1), xx(active) = x
           solves the primal BP problem (see below).
y          is an m-vector of BPdual solution values.
S          is the submatrix A(:,active).
R          is the Cholesky factorization of S.
info       is a structure with the following components:
           .time    total solution time (seconds)
           .stat    = 0  solution is optimal
                    = 1  too many iterations
                    = 2  current S is rank-deficient
                    = 3  dual infeasible point
           .itn     number of iterations
           .exitmsg exit message
           .nprodA  number of products with A
           .nprodAt number of products with A'
           .rNorm   2-norm of the final residual r = b - A(:,active)*x
           .xNorm   1-norm of the primal solution
           .numtrim number of constraints trimmed from final working set

BPdual Toolbox
Copyright 2008, Michael P. Friedlander and Michael A. Saunders
http://www.cs.ubc.ca/labs/scl/bpdual

The primal BP problem may be written

BP:        min_{x,y} ||x||_1 +  1/2 lambda||y||_2^2
              st         Ax     +      lambda  y  =  b.

If lambda > 0, we have the residual vector r = lambda y = b - Ax,
and BP is equivalent to the original BP-denoising problem

   BPDN:      min_{x,y} lambda ||x||_1 +  1/2 ||r||_2^2
              st         Ax     +      r  =  b.

09 Jun 2007: First version of BPdual.
             Michael Friedlander and Michael Saunders.
             A can be dense or sparse (but not an operator).
             S = A(:,active) holds active columns.
             Least-squares subproblems solved by x = S\g.
29 Aug 2007: "A" is now an operator.
09 Nov 2007: Generalize the l1 dual constraints from -e <= A'y <= e
             to bl <= A'y <= bu.
"""
function bpdual(# Positional arguments
                A::Union{AbstractMatrix,AbstractLinearOperator},
                b::Vector,
                λin::Real,
                bl::Vector,
                bu::Vector;

                # Keyword arguments
                inform::OptimResults = OptimResults(size(A)...),
                coldstart::Bool = true,
                trim::Int = 0,
                itnMax::Int = 10*maximum(size(A)),
                feaTol::Real = 5e-05,
                λmin::Real = √eps(1.),
                optTol::Real = 1e-05,
                gapTol::Real = 1e-06,
                pivTol::Real = 1e-12,
                actMax::Real = Inf,
                loglevel::Int = 1)

    time0 = time()

    # ------------------------------------------------------------------
    # Check input arguments.
    # ------------------------------------------------------------------
    m, n = size(A)

    active = inform.active
    state = inform.state
    y = inform.y
    S = inform.S
    R = inform.R

    coldstart = isempty(S)

    # ------------------------------------------------------------------
    # Grab input options and set defaults where needed.
    # ------------------------------------------------------------------
    tieTol = feaTol + 1e-8        # Perturbation to break ties
    λ = max(λin, λmin)

    # ------------------------------------------------------------------
    # Print log header.
    # ------------------------------------------------------------------
    if loglevel > 0
        @fprintf("\n")
        @fprintf(" %s\n",repeat("=",80))
        @fprintf(" %s\n","BPdual")
        @fprintf(" %s\n",repeat("=",80))
        @fprintf(" %-20s: %8i %5s"    ,"No. rows"          ,m       ,"")
        @fprintf(" %-20s: %8.2e\n"    ,"λ"                 ,λ          )
        @fprintf(" %-20s: %8i %5s"    ,"No. columns"       ,n       ,"")
        @fprintf(" %-20s: %8.2e\n"    ,"Optimality tol"    ,optTol     )
        @fprintf(" %-20s: %8i %5s"    ,"Maximum iterations",itnMax  ,"")
        @fprintf(" %-20s: %8.2e\n"    ,"Duality tol"       ,gapTol     )
        @fprintf(" %-20s: %8i %5s"    ,"Support trimming"  ,trim,    "")
        @fprintf(" %-20s: %8.2e\n"    ,"Pivot tol"         ,pivTol     )
    end

    # ------------------------------------------------------------------
    # Initialize local variables.
    # ------------------------------------------------------------------
    EXIT_STATUS =
        Dict(:exit_optimal=>"optimal solution found -- full Newton step",
             :exit_singular_ls=>"singular least-squares subproblem",
             :exit_iterations=>"too many iterations",
             :exit_infeasible=>"dual infeasible point",
             :exit_requested=>"user requested exit",
             :exit_actmax=>"max no. of active constraints reached",
             :exit_small_gap=>"optimal solution found -- small duality gap",
             :exit_unknown=>"unknown exit")
    eFlag = :exit_unknown
    itn       = 0
    step      = 0.0
    p         = 0                         # index of added constraint
    q         = 0                         # index of deleted constraint
    svar      = ""                        # ...and their string value.
    r         = zeros(m)
    zerovec   = zeros(n)
    nBrks     = 0
    numtrim   = 0
    nprodA    = 0
    nprodAt   = 0

    # ------------------------------------------------------------------
    # Cold/warm-start initialization.
    # ------------------------------------------------------------------
    if coldstart
        x = zeros(0)
        z = zeros(n)

    else
        # Trim active-set components that correspond to huge bounds,
        # and then find the min-norm change to y that puts the active
        # constraints exactly on their bounds.
        g = b - λ*y              # Compute steepest-descent dir
        triminf!(active,state,S,R,bl,bu,g)
        nact = length(active)
        x = zeros(nact)
        y = restorefeas(y, active, state, S, R, bl, bu)
        z = A'*y; nprodAt = nprodAt+1

    end

    # Make sure initial point is feasible.
    # bl ≤ z ≤ bu  -->  bl - z ≤ 0
    #              -->  z - bu ≤ 0
    for j in 1:n
        infeasible = (bl[j] - z[j] > feaTol) || (z[j] - bu[j] > feaTol)
        if infeasible
            eFlag = :exit_infeasible
        end
    end

    # ------------------------------------------------------------------
    # Main loop.
    # ------------------------------------------------------------------
    while true

        g = b - λ*y            # Steepest-descent direction

        condS = condest(R)

        if condS > 1e+10
            eFlag = :exit_singular_ls
            # Pad x with enough zeros to make it compatible with S.
            resize!(x, size(S,2))

        else
            x, dy = newtonstep(S, R, g, x, λ)

        end
        r = b - S*x

        # ------------------------------------------------------------------
        # Print to log.
        # ------------------------------------------------------------------
        yNorm = norm(y,2)
        rNorm = norm(r,2)
        xNorm = norm(x,1)

        pObj, dObj, rGap = objectives(x,y,active,b,bl,bu,λ,rNorm,yNorm)
        nact = length(active)   # size of the working set

        if loglevel > 0
            if q != 0
                svar = @sprintf("%8i%1s", q, "-")
            elseif p != 0
                svar = @sprintf("%8i%1s", p ," ")
            else
                svar = @sprintf("%8s%1s"," "," ")
            end
            if mod(itn,50)==0
                @fprintf("\n %4s  %8s %9s %7s %11s %11s %11s %7s %7s\n",
                            "Itn","step","Add/Drop","Active","rNorm2",
                            "xNorm1","Objective","RelGap","condS")
            end
            @fprintf(
            " %4i  %8.1e %8s %7i %11.4e %11.4e %11.4e %7.1e %7.1e\n",
                    itn, step, svar, nact, rNorm, xNorm,
                    dObj, rGap, condS)
        end

        # ------------------------------------------------------------------
        # Check exit conditions.
        # ------------------------------------------------------------------
        if eFlag == :exit_unknown && rGap < gapTol && nact > 0
            eFlag = :exit_small_gap
        end
        if eFlag == :exit_unknown && itn >= itnMax
            eFlag = :exit_iterations
        end
        if eFlag == :exit_unknown && nact >= actMax
            eFlag = :exit_actmax
        end

        # If this is an optimal solution, trim multipliers before exiting.
        if eFlag == :exit_optimal || eFlag == :exit_small_gap
            if trim == 1
                # Optimal trimming. λ0 may be different from λ.
                # Recompute gradient just in case.
                g = b - λin*y
                trimx!(x,S,R,active,state,g,b,λ,feaTol,optTol,loglevel)
                numtrim = nact - length(active)
                nact = length(active)

            elseif trim == 2
                # Threshold trimming.
                # Not yet implemented.
            end
        end

        # Act on any live exit conditions.
        if eFlag != :exit_unknown
            break
        end

        # --------------------------------------------------------------
        # New iteration starts here.
        # --------------------------------------------------------------
        itn += 1
        p = q = 0

        # --------------------------------------------------------------
        # Compute search direction dz.
        # --------------------------------------------------------------
        if norm(dy,Inf) < eps(1.0)       # save a mat-vec prod if dy is 0
            dz = zeros(n)
        else
            dz = A'*dy
            nprodAt += 1
        end

        # ---------------------------------------------
        # Find step to the nearest inactive constraint
        # ---------------------------------------------
        pL, stepL, pU, stepU = step_to_bnd(z, dz, bl, bu, state, tieTol, pivTol)

        step = min(1.0, stepL, stepU)
        hitBnd = step < 1.0
        if hitBnd                    # We bumped into a new constraint
            if step == stepL
                p = pL
                state[p] = -1
            else
                p = pU
                state[p] = +1
            end
        end

        y += step*dy                 # Update dual variables.
        if mod(itn,50)==0
            z = A'*y                 # Occasionally set z = A'y directly.
            nprodAt += 1
        else
            z += step*dz
        end

        if hitBnd
            zerovec[p] = 1            # Extract a = A(:,p)
            a = A*zerovec
            nprodA += 1
            zerovec[p] = 0
            R = qraddcol(S,R,a)       # Update R
            S = [S a]                 # Expand S, active
            push!(active, p)
            push!(x, 0)

        else
            #------------------------------------------------------
            # step = 1.  We moved to a minimizer for the current S.
            # See if we can delete an active constraint.
            #------------------------------------------------------
            drop = false
            if length(active) > 0
                dropl = (state[active].==-1) & (x .> +optTol)
                dropu = (state[active].==+1) & (x .< -optTol)
                dropa = dropl | dropu
                drop = any(dropa)
            end

            # Delete an active constraint.
            if drop
                nact = length(active)
                qa = indmax(abs(x.*dropa))
                q = active[qa]
                state[q] = 0
                S = S[:, 1:nact .!= qa]
                deleteat!(active, qa)
                deleteat!(x, qa)
                R = qrdelcol(R,qa)

            else
                eFlag = :exit_optimal

            end

        end # if step < 1

  end # while true

  tottime = time() - time0

  if loglevel > 0
      @fprintf("\n EXIT BPdual -- %s\n\n",EXIT_STATUS[eFlag])
      @fprintf(" %-20s: %8i %5s","No. significant nnz",sparsity(x),"")
      @fprintf(" %-20s: %8i\n","Products with A",nprodA)
      @fprintf(" %-20s: %8i %5s","No. trimmed nnz",numtrim,"")
      @fprintf(" %-20s: %8i\n","Products with At",nprodAt)
      @fprintf(" %-20s: %8.1e %5s","Solution time (sec)",tottime,"")
      @fprintf("\n")
  end

   # Gather exit data.
   return (x, r, OptimResults(itn,
                              tottime,
                              eFlag,
                              EXIT_STATUS[eFlag],
                              nprodA,
                              nprodAt,
                              numtrim,
                              active, state, y, S, R)
           )
end # function bpdual

# ----------------------------------------------------------------------
# Find step to first bound.
# ----------------------------------------------------------------------
function step_to_bnd(z, dz, bl, bu, state, tieTol, pivTol)

    n = length(dz)
    stepL = stepLtie = stepU = stepUtie = Inf
    pL = pU = 0

    for i in 1:n

        sL = bl[i] -  z[i]
        sU =  z[i] - bu[i]

        # skip variables that aren't free.
        state[i] != 0 && continue

        # variable moves to LOWER bound.
        if dz[i] < -pivTol
            tmp = (sL-tieTol)/dz[i]
            if tmp < stepLtie
                stepLtie = tmp
                stepL = max(0, sL/dz[i])
                pL = i
            end
        end

        # variable moves to UPPER bound.
        if dz[i] > +pivTol
            tmp = -(sU-tieTol)/dz[i]
            if tmp < stepUtie
                stepUtie = tmp
                stepU = max(0, -sU/dz[i])
                pU = i
            end
        end

    end # if

    return (pL, stepL, pU, stepU)

end

# ----------------------------------------------------------------------
# Compute the primal and dual objective values, and the duality gap:
#
#    DP:  minimize_y   - b'y  +  1/2 λ y'y
#         subject to   bl <= A'y <= bu
#
#    PP:  minimize_x   bl'neg(x) + bu'pos(x) + 1/2 λ y'y
#         subject to   Ax + λ y = b.
# ----------------------------------------------------------------------
function objectives(x, y, active,
                    b, bl, bu,
                    λ, rNorm, yNorm)

    bigNum = 1e20

    if isempty(x)
        blx = 0.
        bux = 0.
    else
        blx = bl[active]
        blx[blx .< -bigNum] = 0
        blx = dot(blx, min(x,0))

        bux = bu[active]
        bux[bux .>  bigNum] = 0
        bux = dot(bux, max(x,0))
    end

    if λ > eps(1.0)
        pObj = blx + bux + rNorm^2/2    # primal objective
    else
        pObj = blx + bux
    end
    dObj  = λ*yNorm^2/2 - dot(b,y)         # dual   objective
    maxpd = max(1, pObj, dObj)
    rGap  = abs(pObj+dObj)/maxpd               # relative duality gap

    return (pObj, dObj, rGap)

end # function objectives

# ----------------------------------------------------------------------
# Compute a Newton step.  This is a step to a minimizer of the EQP
#
#   min   g'dy + 1/2 λ dy'dy   subj to   S'dy = 0.
#
# The optimality conditions are given by
#
#   [ -λ I   S ] [dy] = [ h ], with  h = b - lam y - S x,
#   [   S'     ] [dx]   [ 0 ]
#
# where x is an estimate of the Lagrange multiplier. Thus, dx solves
# min ||S dx - h||. On input, g = b - lam y. Alternatively, solve the LS problem
# min ||S  x - g||.
# ----------------------------------------------------------------------
function newtonstep(S, R, g, x, λ)
    m, n = size(S)
    if m==0 || n==0
        dy = g/λ            # Steepest descent
        return (x, dy)
    end
    x, dr = csne(R,S,g)     # LS problem
    if m > n                # Overdetermined system
        dy = dr/λ           # dy is the scaled residual
    else                    # System is square or underdetermined
        dy = zeros(m)       # Anticipate that the residual is 0
    end
    return (x, dy)
end # function newtonstep

# ----------------------------------------------------------------------
# Compute the infeasibility of z relative to bounds bl/bu.
#
#  (bl - z) <= 0  implies  z is   feasible wrt to lower bound
#           >  0  ...           infeasible ...
#  (z - bu) <= 0  implies  z is   feasible wrt to lower bound.
#           >  0  ...           infeasible ...
function infeasibilities(bl::Vector, bu::Vector, z::Vector)
    sL = bl - z
    sU = z - bu
    return (sL, sU)
end

# ----------------------------------------------------------------------
# Trim working constraints with "infinite" bounds.
function triminf!(active::Vector, state::Vector, S::Matrix, R::Matrix,
                  bl::Vector, bu::Vector, g::Vector, b::Vector)

    bigbnd = 1e10
    nact = length(active)

    # Generate a list of constraints to trim.
    tlistbl = find( state[active] .== -1 & bl[active] .< -bigbnd )
    tlistbu = find( state[active] .== +1 & bu[active] .> +bigbnd )
    tlist   = [tlistbl; tlistbu]

    if isempty(tlist)
        return
    end

    for q in tlist

        qa = active[q]          # Index of active constraint to delete
        nact = nact - 1
        S = S[:,1:size(S,2) .!= qa] # Delete column from S
        deleteat!(active, qa)   # Delete index from active set
        R = qrdelcol(R, qa)     # Recompute new QR factorization
        state[q] = 0            # Mark constraint as free
        x = csne(R, S, g)       # Recompute multipliers

        rNorm = norm(b - S*x)
        xNorm = norm(x, 1)
        @fprintf("%4i %4i %4i %10.3e %10.e %10.e %10.e\n",
                k,q,nact,bl[q],bu[q],rNorm,xNorm)
    end

end

# ----------------------------------------------------------------------
# n = sparsity(x,threshold) returns the number of elements in x that
# contain most of the weight of the vector.
function sparsity(x::Vector, threshold::Real = 0.9995)
    threshold = clamp(threshold, 0, 1)
    x = sort(abs(x), rev=true)
    tsum = sum(x)
    csum = 0.0
    j = 0
    while j < length(x)
        csum += x[j+1] / tsum
        if csum >= threshold
            j += 1
        end
    end
    return j
end

# ----------------------------------------------------------------------
# Estimate condition number of R'R.
function condest(R::Matrix)
    if isempty(R)
        condS = 1.0
    else
        n = size(R,1)
        rmin = +Inf
        rmax = -Inf
        for i in 1:n
            rii = R[i,i]
            rmin = rii < rmin ? rii : rmin
            rmax = rii > rmax ? rii : rmax
        end
        condS = rmax / rmin
    end
    return condS
end # function condest

# ----------------------------------------------------------------------
# Trim unneeded constraints from the active set.
#
# `trimx` assumes that the current active set is optimal, ie, 1. x has
# the correct sign pattern, and 2. z := A'y is feasible.
#
# Keep trimming the active set until one of these conditions is
# violated. Condition 2 isn't checked directly. Instead, we simply check
# if the resulting change in y, dy, is small.  It would be more
# appropriate to check that dz := A'dy is small, but we want to avoid
# incurring additional products with A. (Implicitly we're assuming that
# A has a small norm.)
function trimx!(x,S,R,active,state,g,b,λ,featol,opttol,loglevel)

    k = 0
    nact = length(active)
    xabs = abs(x)
    xmin, qa = findmin(xabs)
    gNorm = norm(g,Inf)

    while xmin < opttol

        e = sign(x.*(xabs .> opttol)) # Signs of significant multipliers
        q = active[qa] # Index of the corresponding constraint
        a = S[:,qa]    # Save the col from S in case we need to add it back.
        xsmall = x[qa] # Value of candidate multiplier

        # Trim quantities related to the small multiplier.
        deleteat!(e, qa)
        deleteat!(active, qa)
        S = S[:,1:nact .!= qa]
        R = qrdelcol(R, qa)

        # Recompute the remaining multipliers and their signs.
        (x,dy) = csne(R,S,g)           # min ||g - Sx||_2
        xabs = abs(x)
        et = sign(x.*(xabs .> opttol))
        dyNorm = norm(dy,Inf)/λ        # dy = (g - Sx) / lambda

        # Check if the trimmed active set is still optimal
        if any( et .!= e ) || (dyNorm / max(1,gNorm) > featol)
            R = qraddcol(S,R,a)
            S = [S a]
            push!(active, q)
            x = csne(R, S, g)
            break
        end

        if loglevel > 0 && mod(k,50)==0
            @fprintf("\n %4s  %8s %8s  %7s  %10s  %10s  %10s\n",
                    "Itn","xSmall","Add/Drop","Active","rNorm2",
                    "xNorm1","dyNorm")
        end

        # The trimmed x is still optimal.
        k += 1
        nact -= 1
        state[q] = 0               # Mark this constraint as free.
        rNorm = norm(b - S*x)
        xNorm = norm(x,1)

        # Grab the next canddate multiplier.
        xmin, qa = findmin(xabs)

        # Logging.
        if loglevel > 0
            @fprintf(" %4iT %8.1e %8i- %7i  %10.4e  %10.4e  %10.4e\n",
                    k,xsmall,q,nact,rNorm,xNorm,dyNorm)
        end

    end # while

    return (x,S,R,active,state)

end # function trimx


text:
using Augmentor, ImageTransformations, MNIST, ISICArchive, FixedPointNumbers, Colors, FileIO

mkpath(fname) = joinpath(dirname(@__FILE__()), fname)
srand(123)

toimg(A) = imresize(reinterpret(Gray{N0f8}, UInt8.(A)), 28*4, 28*4)

input_img  = MNIST.trainimage(19)
output_img = augment(input_img, Rotate180())
save(mkpath("readthedocs/background_mnist_in.png"),  toimg(input_img))
save(mkpath("readthedocs/background_mnist_out.png"), toimg(output_img))

input_img  = get(ImageThumbnailRequest(id = "5592ac599fc3c13155a57a85"))
output_img = augment(input_img, Rotate180())
save(mkpath("readthedocs/background_isic_in.png"),  input_img)
save(mkpath("readthedocs/background_isic_out.png"), output_img)

max_stars_repo_name: JuliaTagBot/FileStorage

max_issues_repo_name: JuliaTagBot/FileStorage

max_forks_repo_name: JuliaTagBot/FileStorage


text:
using ModelingToolkit, AbstractTrees, Test

include("../examples/rc_model.jl")

io = IOBuffer()
print_tree(io, rc_model)
ser = String(take!(io))
str = """rc_model
      ├─ resistor
      │  ├─ p
      │  └─ n
      ├─ capacitor
      │  ├─ p
      │  └─ n
      ├─ source
      │  ├─ p
      │  └─ n
      └─ ground
         └─ g
      """
@test strip(ser) == strip(str)


text:
module PlotlyV2Tests

using Plotly

using Base.Test

statuscode(res) = res.status

M = Plotly  # simplify name

@testset "utils" begin
    @test M.basic_auth("foo", "pass") == "Basic Zm9vOnBhc3M="
end

@testset "api -- search" begin
    @test 200 == statuscode(M.search_list_raw("public:false flow"))
end

@testset "api -- files" begin
    # TODO: untested: file_content, file_drop_reference
    # NOTE: Not allowed to test these ones: file_copy, file_permanent_delete

    # TODO: come back to this and change to a plot once I have a plot here!!
    fid = "sglyon:257"  # plot `_jl_testing/myplot`

    get_fid_funcs = [
        M.file_retrieve_raw, M.file_path_raw, M.file_image_raw,
        M.file_sources_raw
    ]
    for func in get_fid_funcs
        @test 200 == statuscode(func(fid))
    end

    res_star_raw = M.file_star_raw(fid)
    @test 200 == statuscode(res_star_raw)
    res = M.get_json_data(res_star_raw)
    @test length(res["stars"]["results"]) > 0
    @test 204 == statuscode(M.file_remove_star(fid))

    @test 200 == statuscode(M.file_lookup_raw("_jl_testing/mygrid"))
    res = M.file_lookup("_jl_testing/mygrid")
    @test res["filename"] == "mygrid"

    @test 200 == statuscode(M.file_trash_raw(fid))
    @test 200 == statuscode(M.file_restore_raw(fid))
    @test 200 == statuscode(M.file_partial_update_raw(fid, world_readable=true))
    @test 200 == statuscode(M.file_update_raw(fid, world_readable=false))
end

@testset "api -- grids" begin
    # TODO: need to test grid_drop_reference, grid_permanent_delete, grid_upload

    fid = "sglyon:251"
    uid1 = "d6434f"
    uid2 = "8c7b60"

    res_retrieve = M.grid_retrieve_raw(fid)
    @test Dict(res_retrieve.headers)["Content-Type"] == "application/json"
    @test 200 == statuscode(res_retrieve)
    json_retrieve = M.get_json_data(res_retrieve)

    @test 200 == statuscode(M.grid_content_raw(fid))

    data = Dict(
        :cols => Dict(
            "first column" => Dict("data" => ["a", "b", "c"], "order" => 0),
            "second column" => Dict("data" => 1:3, "order" => 1),
        )
    )
    del_fn = join(rand('a':'z', 10), "")
    res = M.grid_create_raw(data, parent_path="_jl_testing", filename=del_fn)
    json_delete = M.get_json_data(res)
    fid_delete = json_delete["file"]["fid"]
    uid1_delete = json_delete["file"]["cols"][1]["uid"]
    uid2_delete = json_delete["file"]["cols"][2]["uid"]
    @test statuscode(res) == 201

    @test 201 == statuscode(M.grid_row(fid_delete, rand(2, 2)))

    new_col1 = [Dict("data" => ["d", "e", "f"])]
    res_put_col = M.grid_put_col_raw(fid_delete, new_col1, uid1_delete)

    cols_post = [Dict("name" => "third column", "data" => 4:6)]
    M.grid_post_col(fid_delete, cols_post)

    res2 = M.grid_destroy(fid_delete)
    @test statuscode(res2) == 204

    res_col = M.grid_get_col_raw(fid, uid1)
    @test 200 == statuscode(res_col)
    col_json = M.get_json_data(res_col)
    @test col_json["cols"][1]["data"] == ["a", "b", "c"]

    res_col = M.grid_get_col_raw(fid, uid1, uid2)
    @test 200 == statuscode(res_col)
    col_json = M.get_json_data(res_col)
    @test col_json["cols"][1]["data"] == ["a", "b", "c"]
    @test col_json["cols"][2]["data"] == [1, 2, 3]

    @test 200 == statuscode(M.grid_trash_raw(fid))
    @test 200 == statuscode(M.grid_restore_raw(fid))
    @test 200 == statuscode(M.grid_partial_update_raw(fid, world_readable=true))
    @test 200 == statuscode(M.grid_update_raw(fid, world_readable=false))

    res_lookup = M.grid_lookup_raw("_jl_testing/mygrid")
    @test 200 == statuscode(res_lookup)
    lookup_json = M.get_json_data(res_lookup)
    @test lookup_json["filename"] == "mygrid"
end

@testset "api -- plots" begin
    @test 200 == statuscode(M.plot_list_raw())
    @test 200 == statuscode(M.plot_feed_raw())

    fid = "sglyon:257"

    # test plot_create
    figure = Dict(
        "data" => [Dict(
            "type" => "bar", "xsrc" => "sglyon:251:d6434f", "ysrc" => "sglyon:251:8c7b60"
        )]
    )
    del_fn = join(rand('a':'z', 10), "")
    res_del = M.plot_create_raw(figure, parent_path="_jl_testing", filename=del_fn)
    @test 201 == statuscode(res_del)
    del_json = M.get_json_data(res_del)
    del_fid = del_json["file"]["fid"]

    res_update = M.plot_update_raw(del_fid, world_readable=true)
    @test 200 == statuscode(res_update)
    @test M.get_json_data(res_update)["world_readable"] == true

    res_p_update = M.plot_partial_update_raw(del_fid, world_readable=false)
    @test 200 == statuscode(res_p_update)
    @test M.get_json_data(res_p_update)["world_readable"] == false

    M.file_trash(del_fid)

    res_detail = M.plot_detail_raw(fid)
    @test 200 == statuscode(res_detail)
    @test M.get_json_data(res_detail)["filename"] == "myplot"

    res_content = M.plot_content_raw(fid)
    @test 200 == statuscode(res_content)
    @test M.get_json_data(res_content) == figure
end

@testset "api -- extras" begin
    # TODO
end
@testset "api -- folders" begin
    # TODO
end
@testset "api -- images" begin
    # TODO
end
@testset "api -- comments" begin
    # TODO
end
@testset "api -- plot-schema" begin
    # TODO
end

end


text:
module Domains2D

using ..Boundaries
using ..DielectricFunctions
using ..PumpFunctions
using ..Dispersions
using ..Lattices
using ..Points
using ..Shapes
using RecipesBase

import ..Symmetric, ..Unsymmetric
import LinearAlgebra.norm

import ..LatticeDomain

function LatticeDomain(
            boundary::Boundary{2},
            lattice::Lattice{2,Cartesian},
            n::Number = 1;
            type::Symbol = :generic,
            name::Symbol = :anonymous,
            fit::Bool = boundary.shape.ϕ ≈ lattice.ϕ)

    bnd = boundary
    lat = lattice

    if (typeof(bnd.shape) <: AbstractRectangle) && fit
        nx = ceil(Int,bnd.shape.a/lat.dx)
        ny = ceil(Int,bnd.shape.b/lat.dy)
        dx = bnd.shape.a/nx
        dy = bnd.shape.a/ny
        lat = Lattices.Cartesian(dx, dy; origin=bnd.shape.origin)
        CLASS = Symmetric()
    else
        CLASS = Unsymmetric()
    end

    return lattice_domain_classy(boundary,lattice,n,CLASS;type=type,name=name)
end

function LatticeDomain(
            boundary::Boundary{2},
            lattice::Lattice{2,Polar},
            n::Number = 1;
            type::Symbol = :generic,
            name::Symbol = :anonymous,
            fit::Bool = boundary.shape.origin ≈ lattice.origin)

    bnd = boundary
    lat = lattice

    if (typeof(bnd.shape) <: AbstractUniformDisk) && fit
        nr = ceil(Int,bnd.shape.r/lat.dr)
        nϕ = ceil(Int,2π/lat.dϕ)
        dr = bnd.shape.r/nr
        dϕ = 2π/nϕ
        lat = Lattices.Polar(dr, dϕ; origin=bnd.shape.origin)
        CLASS = Symmetric()
    else
        CLASS = Unsymmetric()
    end

    return lattice_domain_classy(boundary,lattice,n,CLASS;type=type,name=name)
end

function lattice_domain_classy(
            boundary::Boundary{2},
            lattice::Lattice{2},
            n::Number,
            ::CLASS;
            type::Symbol = :generic,
            name::Symbol = :anonymous
            ) where CLASS

    bnd = boundary
    lat = lattice

    ijminmax = _generate_spanning_indices(bnd,lat)
    p = _generate_points(bnd,lat,ijminmax...)
    interior = _generate_interior(bnd,lat,p)
    translations = _generate_translations(bnd,lat,ijminmax...)
    surface = _generate_surface!(interior,bnd,lat,p,translations...)
    nnxm, nnxp, nnym, nnyp = _generate_neighbor_indices(bnd,lat,interior)
    corner = _generate_corner(bnd,lat,p,surface,interior)

    inds_keep = findall(interior)
    x = p[inds_keep]
    interior = interior[inds_keep]
    surface = surface[inds_keep]
    bulk = interior .& .!surface
    corner = corner[inds_keep]
    nnxm = nnxm[inds_keep]
    nnxp = nnxp[inds_keep]
    nnym = nnym[inds_keep]
    nnyp = nnyp[inds_keep]
    indices = map(z->z+CartesianIndex(ijminmax[1]-1,ijminmax[3]-1),inds_keep)

    return LatticeDomain{CLASS}(bnd, lat, n, x, indices, interior,
                            bulk, surface, corner, (nnxm,nnym), (nnxp,nnyp); type=type, name=name)
end

function Base.propertynames(::LatticeDomain{2}, private=false)
    if private
        return fieldnames(LatticeDomain)
    else
        return (:boundary, :lattice, :shape, :n, :ε, :type, :name, :x, :interior, :surface, :corner, :bulk)
    end
end

################################################################################
# DOMAIN BUILDING UTITLITIES

# cartesian spanning indices
function _generate_spanning_indices(bnd::Boundary{2}, lat::Lattice{2,Cartesian})
    i,j = latticeindex(lat,bnd.shape.origin)
    imin, imax = floor(Int,i), ceil(Int,i)
    jmin, jmax = floor(Int,j), ceil(Int,j)

    flag_xmin = flag_xmax = flag_ymin = flag_ymax = true
    init_flag = false
    while flag_ymin || flag_ymax || flag_xmin|| flag_xmax
        I = imin:imax
        J = jmin:jmax
        PBool = BitArray(undef,2length(I)+2length(J))
        count1 = 0
        count2 = 0
        temp = Vector{Float64}(undef,2)
        for i ∈ I
            count2 += 1
            P = lat[i,J[1]]
            PBool[1+count1+i-I[1]] = bnd.shape(P)
        end
        !init_flag || any(PBool[(count1+1):count2]) ? (jmin -= 1; flag_ymin=true) : flag_ymin=false
        !init_flag ? init_flag = any(PBool[(count1+1):count2]) : nothing
        count1 = count2
        for i ∈ I
            count2 += 1
            P = lat[i,J[end]]
            PBool[1+count1+i-I[1]] = bnd.shape(P)
        end
        !init_flag || any(PBool[(count1+1):count2]) ? (jmax += 1; flag_ymax=true) : flag_ymax=false
        !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
        count1 = count2
        for j ∈ J
            count2 += 1
            P = lat[I[1],j]
            PBool[1+count1+j-J[1]] = bnd.shape(P)
        end
        !init_flag || any(PBool[(count1+1):count2]) ? (imin -= 1; flag_xmin=true) : flag_xmin=false
        !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
        count1 = count2
        for j ∈ J
            count2 += 1
            P = lat[I[end],j]
            PBool[2length(I) + length(J) + j+1-J[1]] = bnd.shape(P)
        end
        !init_flag || any(PBool[(count1+1):count2]) ? (imax += 1; flag_xmax=true) : flag_xmax=false
        !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
    end
    return imin-2, imax+2, jmin-2, jmax+2
end

# polar spanning indices
function _generate_spanning_indices(bnd::Boundary, lat::Lattice{2,Polar})
    @assert typeof(bnd.shape)<:Union{Circle,Annulus} "polar coordinates only for circular outermost region"
    @assert lat.x0==bnd.shape.x0 && lat.y0==bnd.shape.y0 "polar lattice origin $((lat.x0,lat.y0)) must be same as circle origin $((bnd.shape.x0,bnd.shape.y0))"
    # i,j = get_lattice_index(lat,bnd.shape.x0,bnd.shape.y0)
    # imin, imax = floor(Int,i), ceil(Int,i)
    # jmin, jmax = floor(Int,j), ceil(Int,j)
    # imin ≤ 0 ? imin = 0 : nothing
    # jmin ≤ 0 ? jmin = 0 : nothing
    imin = 0
    imax = 1
    jmin = 0
    jmax = floor(Int,2π/lat.dθ)

    flag_rmin = flag_rmax = flag_θmin = flag_θmax = true
    init_flag = false
    while flag_rmin || flag_rmax || flag_θmin|| flag_θmax
        I = imin:imax
        J = jmin:jmax
        XYBool = BitArray(undef,2length(I)+2length(J))
        count1 = 0
        count2 = 0
        temp = Array{Float64}(undef,2)
        for i ∈ I
            count2 += 1
            X, Y = lat(i,J[1])
            XYBool[1+count1+i-I[1]] = bnd.shape(X,Y)
        end
        (!init_flag || any(XYBool[(count1+1):count2])) && jmin≥1 ? (jmin -= 1; flag_θmin=true) : flag_θmin=false
        !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
        count1 = count2
        for i ∈ I
            count2 += 1
            X,Y = lat(i,J[end])
            XYBool[1+count1+i-I[1]] = bnd.shape(X,Y)
        end
        (!init_flag || any(XYBool[(count1+1):count2])) && (jmax+1)*lat.dθ≤2π ? (jmax += 1; flag_θmax=true) : flag_θmax=false
        !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
        count1 = count2
        for j ∈ J
            count2 += 1
            X,Y = lat(I[1],j)
            XYBool[1+count1+j-J[1]] = bnd.shape(X,Y)
        end
        (!init_flag || any(XYBool[(count1+1):count2])) && imin≥1 ? (imin -= 1; flag_rmin=true) : flag_rmin=false
        !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
        count1 = count2
        for j ∈ J
            count2 += 1
            X,Y = lat(I[end],j)
            XYBool[2length(I) + length(J) + j+1-J[1]] = bnd.shape(X,Y)
        end
        !init_flag || any(XYBool[(count1+1):count2]) ? (imax += 1; flag_rmax=true) : flag_rmax=false
        !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
    end
    return imin, imax, jmin, jmax
end

function _generate_points(bnd::Boundary{2},lat::Lattice{2},imin,imax,jmin,jmax)
    I, J = imin:imax, jmin:jmax
    return lat[[i for i ∈ I, j ∈ J],[j for i ∈ I, j ∈ J]]
end

function _generate_interior(bnd::Boundary{2},lat::Lattice{2},p)
    interior = BitArray(undef,size(p)...)
    for k ∈ CartesianIndices(p) interior[k] = bnd.shape(p[k]) end
    return interior
end

function _generate_translations(bnd::Boundary{2},lat::Lattice{2},imin,imax,jmin,jmax)
    I, J = imin:imax, jmin:jmax
    p1 = Matrix{Point{2,Cartesian}}(undef,length(I),length(J))
    p2 = Matrix{Point{2,Cartesian}}(undef,length(I),length(J))
    p3 = Matrix{Point{2,Cartesian}}(undef,length(I),length(J))
    p4 = Matrix{Point{2,Cartesian}}(undef,length(I),length(J))
    for k ∈ CartesianIndices(p1)
        i,j = k[1],k[2]
        p1[k] = Cartesian(lat[I[i]-1,J[j]])
        p2[k] = Cartesian(lat[I[i]+1,J[j]])
        p3[k] = Cartesian(lat[I[i],J[j]-1])
        p4[k] = Cartesian(lat[I[i],J[j]+1])
    end
    return p1, p2, p3, p4
end

function _generate_surface!(interior::BitArray,bnd::Boundary{2},lat::Lattice{2},p,t1,t2,t3,t4)
    surface = falses(size(interior)...)
    surface_temp = BitArray(undef,size(interior)...)
    for k ∈ CartesianIndices(interior)
        surface_temp[k] =   interior[k] &&
                            (!bnd.shape(t1[k]) ||
                             !bnd.shape(t2[k]) ||
                             !bnd.shape(t3[k]) ||
                             !bnd.shape(t4[k])
                            )
        if surface_temp[k]
            # n,t,d = normal_distance(bnd.shape,x[k],y[k])
            # if d[1]<min(lat.dx*SURFACE_BUFFER_FACTOR,lat.dy*SURFACE_BUFFER_FACTOR)
                # interior[k] = false
                # surface[k] = false
                # surface[k+CartesianIndex( 1, 0)]=true
                # surface[k+CartesianIndex(-1, 0)]=true
                # surface[k+CartesianIndex( 0, 1)]=true
                # surface[k+CartesianIndex( 0,-1)]=true
            # else
                surface[k] = true
            # end
        end
    end
    return surface
end

function _generate_neighbor_indices(bnd::Boundary{2},lat::Lattice{2,Cartesian},interior)
    NN = LinearIndices(interior) - reshape(cumsum(.!interior[:]),size(interior)...)
    nnxm = vcat(zeros(Int,1,size(interior,2)),NN[1:end-1,:])
    nnxp = vcat(NN[2:end,:],zeros(Int,1,size(interior,2)))
    nnym = hcat(zeros(Int,size(interior,1),1),NN[:,1:end-1])
    nnyp = hcat(NN[:,2:end],zeros(Int,size(interior,1),1))
    return nnxm,nnxp,nnym,nnyp
end

function _generate_neighbor_indices(bnd::Boundary{2},lat::Lattice{2,Polar},interior)
    NN = LinearIndices(interior) - reshape(cumsum(.!interior[:]),size(interior)...)
    nnrm = vcat(zeros(Int,1,size(interior,2)),NN[1:end-1,:])
    nnrp = vcat(NN[2:end,:],zeros(Int,1,size(interior,2)))
    nnθm = hcat(NN[:,end],NN[:,1:end-1])
    nnθp = hcat(NN[:,2:end],NN[:,1])
    return nnrm,nnrp,nnθm,nnθp
end

function _generate_corner(bnd::Boundary{2},lat::Lattice{2,Cartesian},p,surface,interior)
    corner = falses(size(surface)...)
    c = bnd.shape.corners
    dx = hypot(lat.dx,lat.dy)#*(1+SURFACE_BUFFER_FACTOR)
    if length(c)>0
        for k ∈ CartesianIndices(surface)
            if surface[k]
                d = map(z->norm(p[k]-z),c)
                corner[k] = minimum(d) ≤ dx+eps(dx)
                P1 = p[k+CartesianIndex( 2, 0)]
                P2 = p[k+CartesianIndex(-2, 0)]
                P3 = p[k+CartesianIndex( 0, 2)]
                P4 = p[k+CartesianIndex( 0,-2)]
                count = sum(
                            (!bnd.shape(P1) && !interior[k+CartesianIndex( 1, 0)],
                             !bnd.shape(P2) && !interior[k+CartesianIndex(-1, 0)],
                             !bnd.shape(P3) && !interior[k+CartesianIndex( 0, 1)],
                             !bnd.shape(P4) && !interior[k+CartesianIndex( 0,-1)],)
                            )
                corner[k] = corner[k] && (count>1)
            end
        end
    end
    return corner
end

_generate_corner(bnd::Boundary{2},lat::Lattice{2,Polar},p,surface,interior) = falses(size(surface)...)

################################################################################
# PLOTTING

import ...Common.MARKERSIZE_SCALE
import ...Common.MARKERSHAPE

@recipe function f(d::LatticeDomain{2})
    aspect_ratio --> 1
    legend --> false
    @series begin
        seriestype --> :scatter
        markersize --> MARKERSIZE_SCALE/sqrt(length(d.x))
        markerstrokealpha --> 0
        shape --> MARKERSHAPE
        d.x[d.bulk]
    end
    @series begin
        seriestype --> :scatter
        markersize --> MARKERSIZE_SCALE/sqrt(length(d.x))
        markerstrokealpha --> 0
        shape --> MARKERSHAPE
        d.x[d.surface]
    end
    @series begin
        seriestype --> :scatter
        markersize --> MARKERSIZE_SCALE/sqrt(length(d.x))
        markerstrokealpha --> 0
        shape --> MARKERSHAPE
        d.x[d.corner]
    end
    # @series d.boundary
end

end

using .Domains2D


text:
__precompile__()

module Simpler

# package code goes here
include("simple.jl")

end # module


text:
println("Dg) Test of the PhotoIonizationFluores module with ASF from an internally generated initial-, intermediate and final-state multiplets.")
#
wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=JAC.Radial.Grid(true), nuclearModel=Nuclear.Model(26.),
                        initialConfigs=[Configuration("1s^2 2s"), Configuration("1s^2 2p")],
                        intermediateConfigs=[Configuration("1s 2s^2 2p"), Configuration("1s 2s 2p^2") ],
                        finalConfigs  =[Configuration("1s^2 2s^2"), Configuration("1s^2 2s 2p") ],
                        process = PhotoIonFluor(),
                        processSettings=PhotoIonizationFluores.Settings()  )

wb = perform(wa)

## [E1, M1], [JAC.UseCoulomb, JAC.UseBabushkin], false, false, false, true, true, Tuple{Int64,Int64,Int64}[(1,1,0), (1,2,0), (1,2,0)], 0., 0., 0., CoulombInteraction()


text:
""" Misc math functions used in FreeFermionGutzwiller.jl"""

function get_Determinant(r,states)
    """ Computes slater determinant |ϕ_j(r_i)|"""
    mat = states[r,:]
    return det(mat)
end

function get_Inverse_Matrix(r,states)
    """ Compute the inverse of the matrix ϕ_j(r_i) """
    mat = states[r,:]
    return inv(mat)
end

using LinearAlgebra
function det_ratio_factor(A_inv,u,v)
    """ Compute ratio of determinants using Matrix Determinant Lemma"""
    # extension for higher rank update
    correction = hcat(v'*A_inv*u) # turns scalar output into 1 by 1 matrix ¯\_(ツ)_/¯

    one = Matrix{Float64}(I,size(correction))
    return det(one+ correction)
end

function sherman_morrison_inverse_update(A_inv,u,v)
    """ Compute (A')^-1 where A' = A + u v' """
    return A_inv - (A_inv*u*v'*A_inv)/(1+v'*A_inv*u)
end


text:
include("parse.jl")

#Value Stuff

struct closV <: Value
    parms::Array
    body::exprC
    env::Env
end

struct numV <: Value
    val::Number
end

struct strV <: Value
    str::String
end

struct primV <: Value
    op::String
end

struct errorV <: Value
end

struct boolV <: Value
    val::Bool
end

#extends an environment
function extend_env(env, parms, vals)
    nenv = env.bindings
    for i in range(1, length(parms))
        prepend!(nenv, [Binding(parms[i], vals[i])])
    end
    return Env(nenv)
end

#looks up an id in an environment
function env_lookup(find, env)
    for binding in env.bindings
        if cmp(binding.name, find) == 0
            return binding.val
        end
    end
    error("TULI5: name not found in Environment")
end

#fucntion that handles all the two variable operations (primV's)
function binop(op, a1, a2)
    if cmp("+", op) == 0
        if typeof(a1) != numV || typeof(a2) != numV
            error("TULI5 invalid + call")
        end
        return numV(a1.val + a2.val)
    elseif cmp("-", op) == 0
        if typeof(a1) != numV || typeof(a2) != numV
            error("TULI5 invalid - call")
        end
        return numV(a1.val - a2.val)
    elseif cmp("*", op) == 0
        if typeof(a1) != numV || typeof(a2) != numV
            error("TULI5 invalid * call")
        end
        return numV(a1.val * a2.val)
    elseif cmp("/", op) == 0
        if typeof(a1) != numV || typeof(a2) != numV
            error("TULI5 invalid / call")
        end
        if a2.val == 0
            error("TULI5 division by 0")
        end
        return numV(a1.val / a2.val)
    elseif cmp("<=", op) == 0
        if typeof(a1) != numV || typeof(a2) != numV
            error("TULI5 invalid <= call")
        end
        return boolV(a1.val <= a2.val)
    elseif cmp("equal?", op) == 0
        return boolV(a1 == a2)
    end
end

#serializes values into strings
function serialize(Value)
    if typeof(Value) == numV
        return string(Value.val)
    elseif typeof(Value) == boolV
        if Value.val
            return "true"
        else
            return "false"
        end
    elseif typeof(Value) == strV
        return Value.str
    elseif typeof(Value) == closV
        return "<procedure>"
    elseif typeof(Value) == primV
        return "<primative>"
    else
        error("TULI5 unserializable value ", Value)
    end
end

#top interp function
function top_interp(expr)
    return serialize(interp(parse(expr), top_level_env))
end

#Interp function returns a value
function interp(exprC, env)
    if typeof(exprC) == numC
        return numV(exprC.val)
    elseif typeof(exprC) == strC
        return strV(exprC.str)
    elseif typeof(exprC) == ifC
        interped = interp(exprC.val, env)
        if typeof(interped) == boolV
            if interped.val
                return interp(exprC.then, env)
            else
                return interp(exprC.elsee, env)
            end
        else
            error("Invalid TULI5 if statement")
        end
    elseif typeof(exprC) == idC
        return env_lookup(exprC.id, env)
    elseif typeof(exprC) == lamC
        return closV(exprC.parms, exprC.body, env)
    elseif typeof(exprC) == appC
        interped = interp(exprC.fun, env)
        args = exprC.args
        if typeof(interped) == closV
            arg_vals = []
            for arg in args
                append!(arg_vals, [interp(arg, env)])
            end
            if length(interped.parms) == length(arg_vals)
                interp(interped.body, extend_env(interped.env, interped.parms, arg_vals))
            else
                error("number of TULI5 parameters and arguments do not match")
            end
        elseif typeof(interped) == primV
            if length(args) == 2
                binop(interped.op, interp(args[1], env), interp(args[2], env))
            else
                error("Invalid num of args for builtin TULI5 function")
            end
        elseif typeof(interped) == errorV
            if length(args) == 1
                error("TULI 5 User error ", serialize(args[1]))
            else
                error("Invalid num of args for error TULI5 function")
            end
        else
            error("not a TULI5 Function ", exprC)
        end
    end
end

top_level_env = Env(
    [
        Binding("+", primV("+")),
        Binding("-", primV("-")),
        Binding("*", primV("*")),
        Binding("/", primV("/")),
        Binding("<=", primV("<=")),
        Binding("equal?", primV("equal?")),
        Binding("true", boolV(true)),
        Binding("false", boolV(false)),
        Binding("error", errorV())
    ]
)

println(top_interp(4))
println(top_interp(["+", 3, 4]))
println(top_interp(["-", 3, 4]))
println(top_interp(["/", 3, 4]))
println(top_interp(["*", 3, 4]))
println(top_interp([["fn" ["x", "y"] ["*", "x", "y"]], [3, 4]]))

max_stars_repo_name: nhausmna/julia-hello-world

max_issues_repo_name: nhausmna/julia-hello-world

max_forks_repo_name: nhausmna/julia-hello-world