text:
function corners(box::IntervalBox{2,T}) where {T}
    xL = [box[1].lo, box[2].lo]
    xR = [box[1].hi, box[2].hi]
    return xL, xR
end

function split_box(box, numsplits)
    return mince(box, numsplits)
end

function min_diam(box)
    return minimum(diam.(box))
end

function interval_arithmetic_sign_search(
    func,
    initialbox,
    tol,
    perturbation,
    numsplits,
)

    rtol = tol * diam(initialbox)

    foundpos = foundneg = breachedtol = false
    queue = [initialbox]

    while !isempty(queue)
        if (foundpos && foundneg)
            break
        else
            box = popfirst!(queue)
            if min_diam(box) < rtol
                breachedtol = true
                break
            end
            funcrange = func(box)
            if inf(funcrange) > -perturbation
                foundpos = true
            elseif sup(funcrange) < perturbation
                foundneg = true
            else
                newboxes = split_box(box, numsplits)
                push!(queue, newboxes...)
            end
        end
    end

    if breachedtol
        return 2
    elseif foundpos && foundneg
        return 0
    elseif foundpos && !foundneg
        return +1
    elseif !foundpos && foundneg
        return -1
    else
        error("Unexpected scenario")
    end
end

function Base.sign(func, box; tol = 1e-3, perturbation = 0.0, numsplits = 2)
    return interval_arithmetic_sign_search(
        func,
        box,
        tol,
        perturbation,
        numsplits,
    )
end
"""
    sign(f, box)
return
- `+1` if `f` is uniformly positive on `int`
- `-1` if `f` is uniformly negative on `int`
- `0` if `f` has at least one zero crossing in `int` (f assumed continuous)
"""
function Base.sign(func, xL, xR; tol = 1e-3, perturbation = 0.0, numsplits = 2)
    box = IntervalBox(xL, xR)
    return sign(
        func,
        box;
        tol = tol,
        perturbation = perturbation,
        numsplits = numsplits,
    )
end

function extremal_coeffs_in_box(poly, xL, xR)

    max_coeff = -Inf
    min_coeff = Inf
    points = interpolation_points(basis(poly))
    dim, npoints = size(points)
    for i = 1:npoints
        p = view(points, :, i)
        if all(xL .<= p .<= xR)
            coeff = poly.coeffs[i]
            min_coeff = min(min_coeff, coeff)
            max_coeff = max(max_coeff, coeff)
        end
    end
    return max_coeff, min_coeff
end

function Base.sign(
    P::InterpolatingPolynomial{1,B},
    xL,
    xR;
    tol = 1e-3,
    perturbation = 0.0,
    numsplits = 2,
) where {B<:LagrangeTensorProductBasis}

    max_coeff, min_coeff = extremal_coeffs_in_box(P, xL, xR)
    if max_coeff > 0 && min_coeff < 0
        return 0
    else
        box = IntervalBox(xL, xR)
        return interval_arithmetic_sign_search(
            P,
            box,
            tol,
            perturbation,
            numsplits,
        )
    end
end

function (IP::InterpolatingPolynomial{1,B,T})(
    box::IntervalBox{2,S},
) where {B,T,S}
    @assert PolynomialBasis.dimension(IP) == 2
    return IP(box[1], box[2])
end

function PolynomialBasis.gradient(
    IP::InterpolatingPolynomial{1,B,T},
    box::IntervalBox{2,S},
) where {B,T,S}
    @assert PolynomialBasis.dimension(IP) == 2
    return gradient(IP, box[1], box[2])
end

function (IP::InterpolatingPolynomial{1,B,T})(
    box::IntervalBox{1,S},
) where {B<:PolynomialBasis.AbstractBasis{1},T,S}
    @assert PolynomialBasis.dimension(IP) == 1
    return IP(box[1])
end

function PolynomialBasis.gradient(
    IP::InterpolatingPolynomial{1,B,T},
    box::IntervalBox{1,S},
) where {B,T,S}
    @assert PolynomialBasis.dimension(IP) == 1
    return gradient(IP, box[1])
end

max_stars_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_issues_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_forks_repo_name: ArjunNarayanan/ImplicitDomainQuadrature


text:
using Turing
using Test

@model gdemo() = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0, sqrt(s))
  1.5 ~ Normal(m, sqrt(s))
  2.0 ~ Normal(m, sqrt(s))
  return s, m
end

N = 500
s1 = Gibbs(N, HMC(10, 0.1, 5, :s, :m))
s2 = Gibbs(N, PG(10, 10, :s, :m))
s3 = Gibbs(N, PG(10, 2, :s), HMC(1, 0.4, 8, :m))
s4 = Gibbs(N, PG(10, 3, :s), HMC(2, 0.4, 8, :m); thin=false)
s5 = Gibbs(N, CSMC(10, 2, :s), HMC(1, 0.4, 8, :m))

c1 = sample(gdemo(), s1)
c2 = sample(gdemo(), s2)
c3 = sample(gdemo(), s3)
c4 = sample(gdemo(), s4)
c5 = sample(gdemo(), s5)

# Very loose bound, only for testing constructor.
for c in [c1, c2, c3 ,c4, c5]
  check_numerical(c, [:s, :m], [49/24, 7/6], eps=1.0)
end

@test length(c4[:s]) == N * (3 + 2)

# Test gid of each samplers
g = Turing.Sampler(s3)

@test g.info[:samplers][1].alg.gid == 1
@test g.info[:samplers][2].alg.gid == 2

ext: jl


text:
export Companion

struct Companion{T} <: AbstractArray{T, 2}
    c :: Vector{T}
end
#From polynomial

using Polynomials
#Generate companion matrix from a polynomial

function Companion(P::Poly{T}) where T
   n = length(P)
   c = Array{T}(undef, n-1)
   d=P.a[n]
   for i=1:n-1
       c[i]=P.a[i]/d
   end
   Companion(c)
end

#Basic property computations
size(C::Companion, r::Int) = (r==1 || r==2) ? length(C.c) :
    throw(ArgumentError("Companion is of rank 2"))

function size(C::Companion)
    n = length(C.c)
    n, n
end

#XXX Inefficient but works
# getindex(C::Companion, i, j) = getindex(Matrix(C), i, j)
# isassigned(C::Companion, i, j) = isassigned(Matrix(C), i, j)
getindex(C::Companion{T}, i::Int, j::Int) where T = (j==length(C.c)) ? -C.c[i] : (i==j+1 ? one(T) : zero(T) )
isassigned(C::Companion, i::Int, j::Int) = (j==length(C.c)) ? isassigned(C.c,i) : true

function Matrix(C::Companion{T}) where T
    M = zeros(T, size(C)...)
    M[:,end]=-C.c
    for i=1:size(C,1)-1
    	M[i+1, i] = one(T)
	end
    M
end

#Linear algebra stuff
function mul!(C::Companion{T}, b::Vector{T}) where T
	x = b[end]
	y = -C.c[1]*x
	b[2:end] = b[1:end-1]-C.c[2:end]*x
	b[1] = y
    b
end
*(C::Companion{T}, b::Vector{T}) where T = mul!(C, copy(b))

function mul!(A::Matrix{T}, C::Companion{T}) where T
	v = Array{T}(undef, size(A,1))
	for i=1:size(A,1)
		v[i] =dot(vec(A[i,:]),-C.c)
	end
	for i=1:size(A,1), j=1:size(A,2)-1
		A[i,j] = A[i,j+1]
	end
	A[:,end] = v
	A
end
*(A::Matrix{T}, C::Companion{T}) where T = mul!(copy(A), C)

function inv(C::Companion{T}) where T
	M = zeros(T, size(C)...)
    for i=1:size(C,1)-1
    	M[i, i+1] = one(T)
	end
	d = M[end, 1] = -one(T)/C.c[1]
	M[1:end-1, 1] = d*C.c[2:end]
    M
end


text:
module TextObjects
using Match
using ..TextUtils
using REPL.LineEdit
export word, line

const LE = LineEdit

abstract type Selection end

"""
An inner selection
"""
struct Inner <: Selection end
"""
An "all" selection
"""
struct A <: Selection end

# abstract type TextObject end
struct TextObject{T<:Selection}

end

"""
Create a text object
"""
function textobject(buf::IOBuffer, name::String)
    m = match(r"^([ai])(.)$", name)
    selection = @match m[1] begin
        "i" => Inner()
        "a" => A()
    end

    object_fn::Function = @match m[2] begin
        "w" => word
        _ => error("text object command not found: $(m[2])")
    end
    selection(object)
end

function inner()

end

"""
  For non-block objects: For the "a" commands: The operator applies to the object
   and the white space after the object. If there is no white space after the object
    or when the cursor was in the white space before the object, the white space before
    the object is included.
"""
function (::A)(object_function) end

"""
For the "inner" commands: If the cursor was on the object, the operator applies to
the object. If the cursor was on white space, the  operator applies to the white
space.
"""
function (::Inner)(object_function) end

"""
    A word consists of a sequence of letters, digits and underscores, or a
sequence of other non-blank characters, separated with white space (spaces,
tabs, <EOL>).  This can be changed with the 'iskeyword' option.  An empty line
is also considered to be a word.
"""
function word(buf::IOBuffer)::Tuple{Int,Int}
    origin = position(buf)
    local start
    while !is_object_start(buf)
        skip(buf, -1)
    end
    start = position(buf)
    seek(buf, origin)

    local endd
    while !is_object_end(buf)
        skip(buf, 1)
    end
    endd = position(buf)
    seek(buf, origin)
    return (start, endd)
end

function WORD(buf::IOBuffer)::Tuple{Int,Int}

end

function line(buf::IOBuffer)::Tuple{Int,Int}
    # find the line start
    mark(buf)
    if eof(buf)
        if position(buf) > 0
            LE.char_move_left(buf)
        end
    end

    while !eof(buf) && position(buf) > 0
        c = peek(buf, Char)
        if is_linebreak(c)
            skip(buf, 1)
            break
        end
        LE.char_move_left(buf)
    end
    start = position(buf)
    if ismarked(buf)
        reset(buf)
    end

    # find the line end
    mark(buf)
    while !eof(buf)
        c = read(buf, Char)
        if is_linebreak(c)
            LE.char_move_left(buf)
            break
        end
    end
    stop = position(buf)
    reset(buf)

    return (start, stop)
end

"""
    Identify the text object surrounding a space
"""
function space(buf::IOBuffer)::UnitRange{Int}

end
end


text:
using TensorKit,SUNRepresentations
matlab2gt(m,n) = SUNRepresentations.SUNIrrep((m+n,n,0))

a = TensorMap(rand,ComplexF64,RepresentationSpace((matlab2gt(1,1) => 3))*RepresentationSpace((matlab2gt(1,1) => 3))*RepresentationSpace((matlab2gt(0,1) => 3)),RepresentationSpace((matlab2gt(1,1) => 3))*RepresentationSpace((matlab2gt(2,0) => 3)))
@show norm(a);
#a = permute(a,(3,2,5),(1,4));
#a = permute(a,(4,2,1),(5,3));
a = permute(a,(1,),(2,3,4,5));
a = permute(a,(1,2,3),(4,5))
@show norm(a);


text:
## ExoplanetsSysSim/src/stellar_table.jl
## (c) 2015 Eric B. Ford

## TODO Aditya: adapt to TESS (maybe add a new simparam variable indicating mission)
## specifically, the rows to be preserved may need changing

# this is a wrapper around a stellar dataframe

module StellarTable
using ExoplanetsSysSim
#using DataArrays
using DataFrames
using CSV
using JLD
using JLD2
using FileIO

#if VERSION >= v"0.5-"
#  import Compat: UTF8String, ASCIIString
#end

export setup_star_table, star_table, num_usable_in_star_table, set_star_table, star_table_has_key, columns_in_star_table

df = DataFrame()

function setup(sim_param::SimParam; force_reread::Bool = false)
  global df
  wf = WindowFunction.setup_window_function(sim_param)
  WindowFunction.setup_OSD_interp(sim_param) #read in osd files so they can be interpolated
  if haskey(sim_param,"read_stellar_catalog") && !force_reread
     return df
     #return data
  end
  stellar_catalog_filename = convert(String,joinpath(dirname(pathof(ExoplanetsSysSim)), "..","data", convert(String,get(sim_param,"stellar_catalog")) ) )
  df = setup(stellar_catalog_filename)
  add_param_fixed(sim_param,"read_stellar_catalog",true)
  add_param_fixed(sim_param,"num_kepler_targets",num_usable_in_star_table())
  return df
end

function setup(filename::String; force_reread::Bool = false)
  global df
  if occursin(r".jld2$",filename) || occursin(r".jld$",filename)
  #if occursin(r".jld$",filename)
  try
    data = load(filename)
    df = data["stellar_catalog"]
    Core.typeassert(df,DataFrame)
  catch
    error(string("# Failed to read stellar catalog >",filename,"< in jld2 format."))
  end
  else
  try
    #df = readtable(filename)
    #df = CSV.read(filename,nullable=true)
    df = CSV.read(filename, allowmissing=:all)
  catch
    error(string("# Failed to read stellar catalog >",filename,"< in ascii format."))
  end

  # See options at: http://exoplanetarchive.ipac.caltech.edu/docs/API_keplerstellar_columns.html
  # Now we read in all CDPP's, so can interpolate to transit duration
  symbols_to_keep = [ :kepid, :mass, :mass_err1, :mass_err2, :radius, :radius_err1, :radius_err2, :dens, :dens_err1, :dens_err2, :rrmscdpp01p5, :rrmscdpp02p0, :rrmscdpp02p5, :rrmscdpp03p0, :rrmscdpp03p5, :rrmscdpp04p5, :rrmscdpp05p0, :rrmscdpp06p0, :rrmscdpp07p5, :rrmscdpp09p0, :rrmscdpp10p5, :rrmscdpp12p0, :rrmscdpp12p5, :rrmscdpp15p0, :cdppslplong, :cdppslpshrt, :dataspan, :dutycycle, :limbdark_coeff1, :limbdark_coeff2, :limbdark_coeff3, :limbdark_coeff4 ]
  # tess setup does not go through here
  delete!(df, [~(x in symbols_to_keep) for x in names(df)])    # delete columns that we won't be using anyway
  # until I can put in actual TIC limb-darkening coefficients, am setting all of them to zero.
  is_usable = [ !any(ismissing.([ df[i,j] for j in 1:size(df,2) ])) for i in 1:size(df,1) ]
  usable = find(is_usable)
  df = df[usable, symbols_to_keep]
  end
    df[!,:wf_id] = map(x->ExoplanetsSysSim.WindowFunction.get_window_function_id(x,use_default_for_unknown=false),df[!,:kepid])
    obs_5q = df[!,:wf_id].!=-1
    #df = df[obs_5q,keys(df.colindex)]
    df = df[obs_5q,names(df)]
    StellarTable.set_star_table(df)
  return df
end

setup_star_table(sim_param::SimParam; force_reread::Bool = false) = setup(sim_param, force_reread=force_reread)
setup_star_table(filename::String) = setup(filename)

function num_usable_in_star_table()
  global df
  return size(df,1)
end

function columns_in_star_table()
  global df
  return names(df)
end

function star_table(i::Integer, sym::Symbol)
  global df
  return df[i,sym]
end

function star_table(i::Integer)
  global df
  return df[i,:]
end

function star_table(i::Integer, sym::Vector{Symbol})
  global df
  return df[i,sym]
end

function star_table(i::Vector{Integer}, sym::Symbol)
  global df
  return df[i,sym]
end

function star_table(i::Vector{Integer}, sym::Vector{Symbol})
  global df
  return df[i,sym]
end

function set_star_table(df2::DataFrame)
  global df
  df = df2
end

function star_table_has_key(s::Symbol)
  global df
  hasproperty(df,s)
end

end # module StellarTable

# using ExoplanetsSysSim.StellarTable

#=
function generate_star_from_table(sim_param::SimParam, id::Integer)  # WARNING:  To be renamed once there's a working/tested version that uses a stellar catalog with GAIA data
  mu_r = StellarTable.star_table(id,:radius)
  sig_r1 = StellarTable.star_table(id,:radius_err1)
  sig_r2 = StellarTable.star_table(id,:radius_err2)
  z = randn()
  r = mu_r + (z>0) ?  z*sig_r1 : z*sig_r2
  m = rand(Normal(r,0.1))::Float64
  while m<0.0
    m = rand(Normal(r,0.1))::Float64
  end
  f = rand(Normal(1.0,0.1))::Float64
  while f<0.0
    f = 1.0+0.1*randn()
  end
  ld = LimbDarkeningParam4thOrder(StellarTable.star_table(id,:limbdark_coeff1), StellarTable.star_table(id,:limbdark_coeff2), StellarTable.star_table(id,:limbdark_coeff3), StellarTable.star_table(id,:limbdark_coeff4) )
  return SingleStar(r,m,f,id,ld)
  #return SingleStar(r,m,f,id)
end
=#

function generate_star_from_table(sim_param::SimParam)
  id = rand(1:StellarTable.num_usable_in_star_table())
  generate_star_from_table(sim_param, id)
end


text:
module Elefridge

    export CRPS, CRPS!

    using StatsBase

    include("helpers.jl")
    include("nzextrema.jl")
    include("linquantarrays.jl")
    include("logquantarrays.jl")
    include("bitentropy.jl")
    include("bitgrooming.jl")
    include("bitrounding.jl")
    include("bitstring.jl")
    include("bitinformation.jl")
    include("bittranspose.jl")
    include("xor_delta.jl")
    include("crps.jl")
end


text:
# Load packages
using HomotopyContinuation, LinearAlgebra

# Set HC to don't use compilation since the computations are all small
set_default_compile(:none)

# Built index matrix for catalecticant space Cat(r,s)

sym_to_vec(L,m) = [L[i, j] for i = 1:m for j = i:m]

function vec_to_sym(l::AbstractVector{T},m::Int64) where T
    k = 0
    L = Matrix{T}(undef, m, m)
    for i in 1:m, j in i:m
        k += 1
        L[i, j] = L[j, i] = l[k]
    end
    L
end

function CatalecticantIndices(r,s,m)
    tf = [sum(I) == r for I in Iterators.product(ntuple(i->0:r, s)...)]
    indx = collect.(collect(Iterators.product(ntuple(i->0:r, s)...))[tf])
    M = hcat([indx for _ in 1:m]...)
    return M+permutedims(M)
end

function CatalecticantSpace(M::AbstractMatrix,A::AbstractMatrix,m::Integer)
    eq = Vector{typeof(A[1,1]+A[1,1])}(undef,0)
    count = Vector{typeof(M[1,1])}(undef,0)
    ## Parsing the upper triangular matrix
    for i in 1:m
        for j in i:m
            ## Comparing with the remaining entries
            for i2 in (i+1):m
                for j2 in i2:m
                    if M[i,j] == M[i2,j2]
                       #forcing it's not compared to some entry already used to compare
                       if ([i,j] in count)==false
                         #push!(eq, A[i,j]-A[i2,j2]);
                         push!(eq, A[j,i]-A[j2,i2]);
                         push!(count,[i,j]);
                       end
                    end
                end
            end
        end
    end
return eq
end

function adjugate(M::AbstractMatrix{T}) where T
    nr, nc = size(M)
    out = similar(M)
    rows = BitArray(ones(Int16,nr))
    cols = BitArray(ones(Int16,nc))
    for r in 1:nr
        for c in r:nc
            rows[r] = 0
            cols[c] = 0
            out[c, r] = (-1)^(c+r)*det(M[rows,cols])
            rows[r] = 1
            cols[c] = 1
        end
    end
    return out
end

max_stars_repo_name: LauraBMo/InvertingCats

max_issues_repo_name: LauraBMo/InvertingCats

max_forks_repo_name: LauraBMo/InvertingCats


text:
superscript(i) = map(repr(i)) do c
    c   ==  '-' ? '\u207b' :
    c   ==  '1' ? '\u00b9' :
    c   ==  '2' ? '\u00b2' :
    c   ==  '3' ? '\u00b3' :
    c   ==  '4' ? '\u2074' :
    c   ==  '5' ? '\u2075' :
    c   ==  '6' ? '\u2076' :
    c   ==  '7' ? '\u2077' :
    c   ==  '8' ? '\u2078' :
    c   ==  '9' ? '\u2079' :
    c   ==  '0' ? '\u2070' :
    error("Unexpected Character")
end

subscript(i) = map(repr(i)) do c
    c == '0' ? '₀' :
    c == '1' ? '₁' :
    c == '2' ? '₂' :
    c == '3' ? '₃' :
    c == '4' ? '₄' :
    c == '5' ? '₅' :
    c == '6' ? '₆' :
    c == '7' ? '₇' :
    c == '8' ? '₈' :
    c == '9' ? '₉' :
    error("Unexpected Character")
end


text:
# ---
# title: 840. Magic Squares In Grid
# id: problem840
# author: Indigo
# date: 2021-07-09
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/magic-squares-in-grid/description/>
# hidden: true
# ---
#
# A `3 x 3` magic square is a `3 x 3` grid filled with distinct numbers
# **from**`1` **to**`9` such that each row, column, and both diagonals all have
# the same sum.
#
# Given a `row x col` `grid` of integers, how many `3 x 3` "magic square"
# subgrids are there?  (Each subgrid is contiguous).
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/11/magic_main.jpg)
#
#
#
#     Input: grid = [[4,3,8,4],[9,5,1,9],[2,7,6,2]]
#     Output: 1
#     Explanation:
#     The following subgrid is a 3 x 3 magic square:
#     ![](https://assets.leetcode.com/uploads/2020/09/11/magic_valid.jpg)
#     while this one is not:
#     ![](https://assets.leetcode.com/uploads/2020/09/11/magic_invalid.jpg)
#     In total, there is only one magic square inside the given grid.
#
#
# **Example 2:**
#
#
#
#     Input: grid = [[8]]
#     Output: 0
#
#
# **Example 3:**
#
#
#
#     Input: grid = [[4,4],[3,3]]
#     Output: 0
#
#
# **Example 4:**
#
#
#
#     Input: grid = [[4,7,8],[9,5,1],[2,3,6]]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `row == grid.length`
#   * `col == grid[i].length`
#   * `1 <= row, col <= 10`
#   * `0 <= grid[i][j] <= 15`
#
#
## @lc code=start
using LeetCode

function num_magic_squares_inside(grid::Matrix{Int})::Int
    r = [15, 15, 15]
    function is_magic_square(grid::AbstractMatrix{Int})
        sum(grid; dims=1) != reshape(r, 1, 3) && return false
        sum(grid; dims=2) != reshape(r, 3, 1) && return false
        sum(grid[i, i] for i in 1:3) == sum(grid[i, 4 - i] for i in 1:3) == 15 && return true
        return false
    end
    res = 0
    for j in 1:size(grid,2)-2, i in 1:size(grid,1)-2
        sub_mat = @view(grid[CartesianIndex(i, j):CartesianIndex(i+2, j+2)])
        all(x -> x in sub_mat, 1:9) && is_magic_square(sub_mat) && (res += 1)
    end
    return res
end
## @lc code=end


text:
!(dirname(@__FILE__) in LOAD_PATH) && push!(LOAD_PATH, dirname(@__FILE__))
using DSP, Test, FilterTestHelpers

#
# Butterworth filter prototype
#
@testset "Butterworth" begin
    # Poles of 20 pole Butterworth filter prototype from MATLAB 2013b (buttap(20))
    matlab_p = [-0.07845909572784487+0.996917333733128im,
        -0.07845909572784487-0.996917333733128im,
        -0.2334453638559053+0.9723699203976767im,
        -0.2334453638559053-0.9723699203976767im,
        -0.3826834323650897+0.9238795325112867im,
        -0.3826834323650897-0.9238795325112867im,
        -0.5224985647159488+0.8526401643540923im,
        -0.5224985647159488-0.8526401643540923im,
        -0.6494480483301835+0.760405965600031im,
        -0.6494480483301835-0.760405965600031im,
        -0.7604059656000308+0.6494480483301838im,
        -0.7604059656000308-0.6494480483301838im,
        -0.8526401643540922+0.5224985647159489im,
        -0.8526401643540922-0.5224985647159489im,
        -0.9238795325112867+0.3826834323650899im,
        -0.9238795325112867-0.3826834323650899im,
        -0.9723699203976766+0.2334453638559055im,
        -0.9723699203976766-0.2334453638559055im,
        -0.996917333733128+0.07845909572784507im,
        -0.996917333733128-0.07845909572784507im]
    matlab_butter = ZeroPoleGain(Float64[], matlab_p, 1)

    # Test that our answers are close to MATLAB's and at least as accurate
    butter = Butterworth(20)
    zpkfilter_eq(butter, matlab_butter)
    zpkfilter_accuracy(butter, matlab_butter, Butterworth(BigFloat, 20))

    # Poles of 19 pole Butterworth filter prototype from MATLAB 2013b (buttap(19))
    matlab_p = [-0.08257934547233227+0.9965844930066698im,
        -0.08257934547233227-0.9965844930066698im,
        -0.2454854871407991+0.9694002659393304im,
        -0.2454854871407991-0.9694002659393304im,
        -0.4016954246529694+0.9157733266550574im,
        -0.4016954246529694-0.9157733266550574im,
        -0.5469481581224267+0.8371664782625287im,
        -0.5469481581224267-0.8371664782625287im,
        -0.6772815716257409+0.7357239106731317im,
        -0.6772815716257409-0.7357239106731317im,
        -0.7891405093963935+0.6142127126896679im,
        -0.7891405093963935-0.6142127126896679im,
        -0.8794737512064892+0.4759473930370733im,
        -0.8794737512064892-0.4759473930370733im,
        -0.9458172417006346+0.3246994692046836im,
        -0.9458172417006346-0.3246994692046836im,
        -0.9863613034027224+0.1645945902807336im,
        -0.9863613034027224-0.1645945902807336im,
        -1+0im]

    # Test that our answers are close to MATLAB's and at least as accurate
    zpkfilter_eq(Butterworth(19), ZeroPoleGain(Float64[], matlab_p, 1))
    zpkfilter_accuracy(Butterworth(19), ZeroPoleGain(Float64[], matlab_p, 1), Butterworth(BigFloat, 19))
end

#
# Chebyshev type I filter
#
@testset "Chebyshev I" begin
    # Poles of 20 pole Butterworth filter prototype with 1 dB passband ripple from MATLAB 2013b:
    #=
        [z, p, k] = cheb1ap(20, 1)
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
    =#
    matlab_p = [-0.005606643513655412+0.9994594480354867im,
        -0.01668187637027696+0.9748494394091176im,
        -0.02734634606668882+0.9262354022446964im,
        -0.03733745796561948+0.854814375477951im,
        -0.04640919773351807+0.7623449818206758im,
        -0.05433818902992066+0.6511041246330369im,
        -0.06092919377429644+0.5238309230062319im,
        -0.06601991955554824+0.3836592655598524im,
        -0.06948501580979141+0.2340406437032877im,
        -0.07123916036725134+0.07865916446238384im,
        -0.07123916036725134-0.07865916446238384im,
        -0.06948501580979141-0.2340406437032877im,
        -0.06601991955554824-0.3836592655598524im,
        -0.06092919377429644-0.5238309230062319im,
        -0.05433818902992066-0.6511041246330369im,
        -0.04640919773351807-0.7623449818206758im,
        -0.03733745796561948-0.854814375477951im,
        -0.02734634606668882-0.9262354022446964im,
        -0.01668187637027696-0.9748494394091176im,
        -0.005606643513655412-0.9994594480354867im]
    matlab_k = 3.748372513504540e-06

    # Test that our answers are close to MATLAB's and at least as accurate
    zpkfilter_eq(Chebyshev1(20, 1), ZeroPoleGain(Float64[], matlab_p, matlab_k))
    zpkfilter_accuracy(Chebyshev1(20, 1), ZeroPoleGain(Float64[], matlab_p, matlab_k), Chebyshev1(BigFloat, 20, 1))

    # Poles of 19 pole Butterworth filter prototype with 1 dB passband ripple from MATLAB 2013b:
    #=
        [z, p, k] = cheb1ap(19, 1)
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    matlab_p = [-0.006212227114275604+0.9994004289494353im,
        -0.0184672279812169+0.9721393904901743im,
        -0.03021849100930104+0.9183609236365976im,
        -0.04114547237637834+0.8395319647744703im,
        -0.05095011251526942+0.7378027624097928im,
        -0.0593649664017747+0.6159482239948341im,
        -0.06616049875775722+0.4772922236864799im,
        -0.07115134517517649+0.3256169357239265im,
        -0.07420136837393095+0.1650596665748249im,
        -0.07522737167197566+0im,
        -0.07420136837393095-0.1650596665748249im,
        -0.07115134517517649-0.3256169357239265im,
        -0.06616049875775722-0.4772922236864799im,
        -0.0593649664017747-0.6159482239948341im,
        -0.05095011251526942-0.7378027624097928im,
        -0.04114547237637834-0.8395319647744703im,
        -0.03021849100930104-0.9183609236365976im,
        -0.0184672279812169-0.9721393904901743im,
        -0.006212227114275604-0.9994004289494353im]
    matlab_k = 7.496745027009062e-06

    # Test that our answers are close to MATLAB's and at least as accurate
    zpkfilter_eq(Chebyshev1(19, 1), ZeroPoleGain(Float64[], matlab_p, matlab_k))
    zpkfilter_accuracy(Chebyshev1(19, 1), ZeroPoleGain(Float64[], matlab_p, matlab_k), Chebyshev1(BigFloat, 19, 1))
end

#
# Chebyshev type II filter
#
@testset "Chebyshev II" begin
    # 20 pole Butterworth filter prototype with 1 dB passband ripple from MATLAB 2013b:
    #=
        [z, p, k] = cheb2ap(20, 1)
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    matlab_z = [0+1.003092198482826im, 0-1.003092198482826im,
        0+1.028415193665209im, 0-1.028415193665209im,
        0+1.082392200292394im, 0-1.082392200292394im,
        0+1.172827696614009im, 0-1.172827696614009im,
        0+1.315086999890785im, 0-1.315086999890785im,
        0+1.539769043222366im, 0-1.539769043222366im,
        0+1.913880855430943im, 0-1.913880855430943im,
        0+2.613125929752753im, 0-2.613125929752753im,
        0+4.283657569731185im, 0-4.283657569731185im,
        0+12.74549484318238im, 0-12.74549484318238im]
    matlab_p = [-0.001929675700544753-1.002788598382701im,
        -0.006034873757627809-1.028072311003105im,
        -0.010957850142167-1.081957626858823im,
        -0.01756373888305481-1.172213900920771im,
        -0.02744258171352171-1.314120759169777im,
        -0.04403133250736632-1.538048178699878im,
        -0.07621942506983267-1.910267531994771im,
        -0.1536690174908018-2.6032737531287im,
        -0.4316587213861887-4.238414904224276im,
        -3.610083779585366-11.62012078608171im,
        -3.610083779585366+11.62012078608171im,
        -0.4316587213861887+4.238414904224276im,
        -0.1536690174908018+2.6032737531287im,
        -0.07621942506983267+1.910267531994771im,
        -0.04403133250736632+1.538048178699878im,
        -0.02744258171352171+1.314120759169777im,
        -0.01756373888305481+1.172213900920771im,
        -0.010957850142167+1.081957626858823im,
        -0.006034873757627809+1.028072311003105im,
        -0.001929675700544753+1.002788598382701im]
    matlab_k = 0.8912509381337452
    matlab_f = ZeroPoleGain(matlab_z, matlab_p, matlab_k)

    # The gain shows a greater discrepancy with the gain of the enhanced
    # precision filter than MATLAB's, but AFAICT our filter is still more
    # accurate: the poles and zeros are each more accurate, and computing
    # the gain in enhanced precision yields the same value.

    f = Chebyshev2(20, 1)
    zpkfilter_eq(f, matlab_f)
    zpkfilter_accuracy(f, matlab_f, Chebyshev2(BigFloat, 20, 1), relerr=4)

    # 19 pole Butterworth filter prototype with 1 dB passband ripple from MATLAB 2013b:
    #=
        [z, p, k] = cheb2ap(19, 1)
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    matlab_z = [0+1.003427212662145im, 0-1.003427212662145im,
    0+1.031565634068626im, 0-1.031565634068626im,
    0+1.091973276457601im, 0-1.091973276457601im,
    0+1.194505544554793im, 0-1.194505544554793im,
    0+1.359205519207709im, 0-1.359205519207709im,
    0+1.628100459889458im, 0-1.628100459889458im,
    0+2.101072544213108im, 0-2.101072544213108im,
    0+3.079770972368366im, 0-3.079770972368366im,
    0+6.075533820974263im, 0-6.075533820974263im]
    matlab_p = [-0.002139218568677465-1.003090263923361im,
        -0.006720708432694556-1.031180124433753im,
        -0.01232195186345456-1.091472452757043im,
        -0.02007306444044342-1.193772338103026im,
        -0.03217434756636293-1.357992987527419im,
        -0.05375960613144559-1.625783486413104im,
        -0.09966387580015423-2.09563676642751im,
        -0.2295218019875045-3.061543704346429im,
        -0.9149742173123645-5.932401750359961im,
        -38.84170469908301+0im,
        -0.9149742173123645+5.932401750359961im,
        -0.2295218019875045+3.061543704346429im,
        -0.09966387580015423+2.09563676642751im,
        -0.05375960613144559+1.625783486413104im,
        -0.03217434756636293+1.357992987527419im,
        -0.02007306444044342+1.193772338103026im,
        -0.01232195186345456+1.091472452757043im,
        -0.006720708432694556+1.031180124433753im,
        -0.002139218568677465+1.003090263923361im]
    matlab_k = 37.33930783884512
    matlab_f = ZeroPoleGain(matlab_z, matlab_p, matlab_k)

    f = Chebyshev2(19, 1)
    zpkfilter_eq(f, matlab_f)
    zpkfilter_accuracy(f, matlab_f, Chebyshev2(BigFloat, 19, 1), relerr=2)
end

#
# Elliptic filter
#
@testset "elliptic" begin
    # 20 pole elliptic filter prototype with 0.1 dB passband ripple and 10
    # dB stopband ripple from MATLAB 2013b:
    #=
        [z, p, k] = ellipap(20, 0.1, 10)
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    matlab_z = [0-1.953252853757711im, 0+1.953252853757711im,
        0-1.069599937693626im, 0+1.069599937693626im,
        0-1.007032209276402im, 0+1.007032209276402im,
        0-1.000730367252677im, 0+1.000730367252677im,
        0-1.000076070678288im, 0+1.000076070678288im,
        0-1.000007925871874im, 0+1.000007925871874im,
        0-1.000000826312751im, 0+1.000000826312751im,
        0-1.000000086632384im, 0+1.000000086632384im,
        0-1.000000009576065im, 0+1.000000009576065im,
        0-1.000000001633984im, 0+1.000000001633984im]
    matlab_p = [-0.6665552331000151-0.9586923119289487im,
        -0.6665552331000151+0.9586923119289487im,
        -0.06598422536270995-1.015686365686964im,
        -0.06598422536270995+1.015686365686964im,
        -0.006773291804375061-1.001821704927991im,
        -0.006773291804375061+1.001821704927991im,
        -0.0007045345890457818-1.000191780070487im,
        -0.0007045345890457818+1.000191780070487im,
        -7.339097059946433e-05-1.000020003037497im,
        -7.339097059946433e-05+1.000020003037497im,
        -7.646273200288023e-06-1.000002084836051im,
        -7.646273200288023e-06+1.000002084836051im,
        -7.966432964216136e-07-1.000000217756124im,
        -7.966432964216136e-07+1.000000217756124im,
        -8.299912436196738e-08-1.000000023227775im,
        -8.299912436196738e-08+1.000000023227775im,
        -8.637784231770464e-09-1.000000002962719im,
        -8.637784231770464e-09+1.000000002962719im,
        -8.070974080610083e-10-1.000000000874029im,
        -8.070974080610083e-10+1.000000000874029im]
    matlab_k = 0.3162277662398871
    matlab_f = ZeroPoleGain(matlab_z, matlab_p, matlab_k)

    f = Elliptic(20, 0.1, 10)
    zpkfilter_eq(f, matlab_f)
    zpkfilter_accuracy(f, matlab_f, Elliptic(BigFloat, 20, 0.1, 10), eps=1e-9)

    # 19 pole elliptic filter prototype with 0.1 dB passband ripple and 10
    # dB stopband ripple from MATLAB 2013b:
    #=
        [z, p, k] = ellipap(19, 0.1, 10)
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    matlab_z = [0-1.232609672486912im, 0+1.232609672486912im,
        0-1.021948255265862im, 0+1.021948255265862im,
        0-1.002264470627064im, 0+1.002264470627064im,
        0-1.000235691798532im, 0+1.000235691798532im,
        0-1.000024555183434im, 0+1.000024555183434im,
        0-1.000002559984621im, 0+1.000002559984621im,
        0-1.000000268394018im, 0+1.000000268394018im,
        0-1.00000002966741im, 0+1.00000002966741im,
        0-1.000000005062212im, 0+1.000000005062212im]
    matlab_p = [-0.2102804655411482-1.034611529620449im,
        -0.2102804655411482+1.034611529620449im,
        -0.0210635609742992-1.005504639644632im,
        -0.0210635609742992+1.005504639644632im,
        -0.002183581313750937-1.00059265523267im,
        -0.002183581313750937+1.00059265523267im,
        -0.0002273805846473938-1.000061954788962im,
        -0.0002273805846473938+1.000061954788962im,
        -2.368886404726722e-05-1.000006458811389im,
        -2.368886404726722e-05+1.000006458811389im,
        -2.468065209676509e-06-1.00000067462382im,
        -2.468065209676509e-06+1.00000067462382im,
        -2.571378798158274e-07-1.000000071961468im,
        -2.571378798158274e-07+1.000000071961468im,
        -2.676054057433227e-08-1.000000009178736im,
        -2.676054057433227e-08+1.000000009178736im,
        -2.500451772444974e-09-1.00000000270781im,
        -2.500451772444974e-09+1.00000000270781im,
        -1.309071549907812+0im]
    matlab_k = 0.9266824319626478
    matlab_f = ZeroPoleGain(matlab_z, matlab_p, matlab_k)

    f = Elliptic(19, 0.1, 10)
    zpkfilter_eq(f, matlab_f)
    zpkfilter_accuracy(f, matlab_f, Elliptic(BigFloat, 19, 0.1, 10), eps=4e-9)
end

#
# Conversion between zpk and tf
#
@testset "ZPK to TF" begin
    # Output of [z, p, k] = buttap(20); [b, a] = zp2tf(z, p, k)
    m_a = [1, 12.74549484318237, 81.22381939879423, 343.6513712403923,
        1081.352361133001, 2687.409807920676, 5468.931438945091,
        9326.061201886809, 13528.36656744904, 16852.27707949905,
        18122.54155403868, 16852.27707949905, 13528.36656744904,
        9326.061201886809, 5468.931438945092, 2687.409807920676,
        1081.352361133001, 343.6513712403923, 81.22381939879423,
        12.74549484318237, 1]

    f = convert(PolynomialRatio, Butterworth(20))
    @test coefb(f) ≈ [1]
    @test coefa(f) ≈ m_a

    # Test that our answers are more accurate than MATLAB's
    accurate_a = coefa(convert(PolynomialRatio, Butterworth(BigFloat, 20)))
    @test coefa(f) ≈ map(Float64, accurate_a)
    @test sum(abs.(coefa(f) - accurate_a)) <= sum(abs.(m_a - accurate_a))
end

#
# Conversion between tf and zpk
#
@testset "TF to ZPK" begin
    f = convert(ZeroPoleGain, convert(PolynomialRatio, Butterworth(20)))
    zpkfilter_eq(f, Butterworth(20), 1e-6)

    # TODO: Make this more accurate!

    # Output of [z, p, k] = buttap(20); [b, a] = zp2tf(z, p, k); tf2zpk(b, a)
    m_p = [-0.07845909573254482+0.9969173337335029im,
        -0.07845909573254482-0.9969173337335029im,
        -0.2334453637958131+0.9723699203918822im,
        -0.2334453637958131-0.9723699203918822im,
        -0.3826834327796701+0.9238795325396184im,
        -0.3826834327796701-0.9238795325396184im,
        -0.5224985628488221+0.8526401643454914im,
        -0.5224985628488221-0.8526401643454914im,
        -0.6494480541398985+0.7604059651905597im,
        -0.6494480541398985-0.7604059651905597im,
        -0.760405952587916+0.6494480502272874im,
        -0.760405952587916-0.6494480502272874im,
        -0.8526401859847815+0.5224985598169277im,
        -0.8526401859847815-0.5224985598169277im,
        -0.9238795057196649+0.3826834415328767im,
        -0.9238795057196649-0.3826834415328767im,
        -0.9969173244841298+0.07845911266921719im,
        -0.9969173244841298-0.07845911266921719im,
        -0.97236994351794+0.2334453500964366im,
        -0.97236994351794-0.2334453500964366im]
    # println(ComplexF64([sort(f.p, lt=lt) - sort(Butterworth(BigFloat, 20).p, lt=lt) sort(m_p, lt=lt) - sort(Butterworth(BigFloat, 20).p, lt=lt)]))
    zpkfilter_accuracy(f, ZeroPoleGain(Float64[], m_p, 1), Butterworth(BigFloat, 20); eps=1e-6, relerr=9)
end

#
# Frequency scaling
#
@testset "freq. scaling" begin
    # Output of:
    #=
        [z, p, k] = butter(20, 0.5, 's')
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    m_z = Float64[]
    m_p = [-0.498458666866501+0.03922954786337574im,
        -0.498458666866501-0.03922954786337574im,
        -0.4861849601989526+0.1167226819283556im,
        -0.4861849601989526-0.1167226819283556im,
        -0.4619397662555457+0.191341716182291im,
        -0.4619397662555457-0.191341716182291im,
        -0.4263200821771486+0.2612492823580798im,
        -0.4263200821771486-0.2612492823580798im,
        -0.3802029827999479+0.3247240241651102im,
        -0.3802029827999479-0.3247240241651102im,
        -0.3247240241650893+0.3802029827999696im,
        -0.3247240241650893-0.3802029827999696im,
        -0.2612492823580033+0.4263200821770542im,
        -0.2612492823580033-0.4263200821770542im,
        -0.1913417161825348+0.461939766255657im,
        -0.1913417161825348-0.461939766255657im,
        -0.1167226819279489+0.486184960198832im,
        -0.1167226819279489-0.486184960198832im,
        -0.03922954786392431+0.4984586668665641im,
        -0.03922954786392431-0.4984586668665641im]
    m_k = 9.5367431640625e-07
    m_f = ZeroPoleGain(m_z, m_p, m_k)

    f = analogfilter(Lowpass(0.5), Butterworth(20))
    zpkfilter_eq(f, m_f)

    # Test that our answers are more accurate than MATLAB's
    accurate_ft = analogfilter(Lowpass(0.5), Butterworth(BigFloat, 20))
    zpkfilter_eq(f, accurate_ft)
    zpkfilter_accuracy(f, m_f, accurate_ft)

    # Output of
    #=
        [z, p, k] = ellip(20, 0.5, 20, 0.3, 's')
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    m_z = [0+0.6269347368248613im, 0-0.6269347368248613im,
        3.244237104018702e-18+0.3275342448365736im,
        3.244237104018702e-18-0.3275342448365736im,
        -6.938893903907228e-18+0.3032929712822612im,
        -6.938893903907228e-18-0.3032929712822612im,
        -1.387778780781446e-17+0.3004076906327037im,
        -1.387778780781446e-17-0.3004076906327037im,
        -4.163336342344337e-17+0.3000506878920809im,
        -4.163336342344337e-17-0.3000506878920809im,
        6.938893903907228e-18+0.3000063060899686im,
        6.938893903907228e-18-0.3000063060899686im,
        -6.938893903907228e-18+0.3000007854095975im,
        -6.938893903907228e-18-0.3000007854095975im,
        -1.387778780781446e-17+0.30000009864114im,
        -1.387778780781446e-17-0.30000009864114im,
        1.387778780781446e-17+0.3000000132267227im,
        1.387778780781446e-17-0.3000000132267227im,
        2.775557561562891e-17+0.3000000027630185im,
        2.775557561562891e-17-0.3000000027630185im]
    m_p = [-0.1398677695014979+0.188138982052175im,
        -0.1398677695014979-0.188138982052175im,
        -0.02345522399724158+0.2890911154818891im,
        -0.02345522399724158-0.2890911154818891im,
        -0.003009381182870904+0.2987256455960404im,
        -0.003009381182870904-0.2987256455960404im,
        -0.0003757590363774529+0.2998428279105927im,
        -0.0003757590363774529-0.2998428279105927im,
        -4.676623056761031e-05+0.299980469734862im,
        -4.676623056761031e-05-0.299980469734862im,
        -5.818096677288886e-06+0.2999975716762768im,
        -5.818096677288886e-06-0.2999975716762768im,
        -7.237817530730828e-07+0.2999996988517393im,
        -7.237817530730828e-07-0.2999996988517393im,
        -9.003660077222531e-08+0.2999999634684986im,
        -9.003660077222531e-08-0.2999999634684986im,
        -1.117941275663448e-08+0.2999999963792452im,
        -1.117941275663448e-08-0.2999999963792452im,
        -1.220070558471953e-09+0.300000000410979im,
        -1.220070558471953e-09-0.300000000410979im]
    m_k = 0.1000000000276181
    m_f = ZeroPoleGain(m_z, m_p, m_k)

    f = convert(ZeroPoleGain, analogfilter(Lowpass(0.3), Elliptic(20, 0.5, 20)))
    zpkfilter_eq(f, m_f)
end

#
# High pass filter construction
#
@testset "highpass" begin
    # Output of
    #=
        [z, p, k] = butter(20, 0.5, 'high', 's')
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    m_z = zeros(ComplexF64, 20)
    m_p = [-0.4984586668666513+0.03922954786370267im,
        -0.4984586668666513-0.03922954786370267im,
        -0.4861849601986916+0.1167226819280185im,
        -0.4861849601986916-0.1167226819280185im,
        -0.4619397662557133+0.191341716182565im,
        -0.4619397662557133-0.191341716182565im,
        -0.4263200821770351+0.2612492823579622im,
        -0.4263200821770351-0.2612492823579622im,
        -0.3802029828000168+0.3247240241650909im,
        -0.3802029828000168-0.3247240241650909im,
        -0.3247240241650923+0.3802029828000145im,
        -0.3247240241650923-0.3802029828000145im,
        -0.03922954786392298+0.4984586668665641im,
        -0.03922954786392298-0.4984586668665641im,
        -0.2612492823579744+0.4263200821770468im,
        -0.2612492823579744-0.4263200821770468im,
        -0.1913417161825449+0.4619397662556436im,
        -0.1913417161825449-0.4619397662556436im,
        -0.1167226819279526+0.4861849601988384im,
        -0.1167226819279526-0.4861849601988384im]
    m_k = 1
    m_f = ZeroPoleGain(m_z, m_p, m_k)

    f = analogfilter(Highpass(0.5), Butterworth(20))
    zpkfilter_eq(f, m_f)

    # Test that our answers are more accurate than MATLAB's
    accurate_ft = analogfilter(Highpass(0.5), Butterworth(BigFloat, 20))
    zpkfilter_eq(f, accurate_ft)
    zpkfilter_accuracy(f, m_f, accurate_ft)

    # Output of
    #=
        [z, p, k] = ellip(20, 0.5, 20, 0.4, 'high', 's')
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    m_z = [4.289315073341702e-17+0.1914074830304432im,
        4.289315073341702e-17-0.1914074830304432im,
        -1.734723475976807e-17+0.3663739040779542im,
        -1.734723475976807e-17-0.3663739040779542im,
        -1.387778780781446e-17+0.3956570424057772im,
        -1.387778780781446e-17-0.3956570424057772im,
        -2.775557561562891e-17+0.3994571502056489im,
        -2.775557561562891e-17-0.3994571502056489im,
        2.775557561562891e-17+0.3999324275609072im,
        2.775557561562891e-17-0.3999324275609072im,
        8.673617379884035e-18+0.3999915920567793im,
        8.673617379884035e-18-0.3999915920567793im,
        4.163336342344337e-17+0.399998952789945im,
        4.163336342344337e-17-0.399998952789945im,
        1.387778780781446e-17+0.3999998684785239im,
        1.387778780781446e-17-0.3999998684785239im,
        -2.775557561562891e-17+0.3999999963159753im,
        -2.775557561562891e-17-0.3999999963159753im,
        6.938893903907228e-18+0.3999999823643708im,
        6.938893903907228e-18-0.3999999823643708im]
    m_p = [-0.3053922020601836+0.4107892635097733im,
        -0.3053922020601836-0.4107892635097733im,
        -0.03345814372376468+0.4123794380387859im,
        -0.03345814372376468-0.4123794380387859im,
        -0.004046405082073851+0.4016656239382425im,
        -0.004046405082073851-0.4016656239382425im,
        -0.0005015366410947575+0.4002090441157958im,
        -0.0005015366410947575-0.4002090441157958im,
        -6.236309209528557e-05+0.400026032326622im,
        -6.236309209528557e-05-0.400026032326622im,
        -7.757587819075262e-06+0.4000032376407234im,
        -7.757587819075262e-06-0.4000032376407234im,
        -9.650442748948551e-07+0.4000004015290894im,
        -9.650442748948551e-07-0.4000004015290894im,
        -1.200488304875513e-07+0.4000000487086385im,
        -1.200488304875513e-07-0.4000000487086385im,
        -1.490588337482723e-08+0.4000000048276729im,
        -1.490588337482723e-08-0.4000000048276729im,
        -1.626761100825824e-09+0.3999999994520274im,
        -1.626761100825824e-09-0.3999999994520274im]
    m_k = 0.9440608762859236
    m_f = ZeroPoleGain(m_z, m_p, m_k)

    f = convert(ZeroPoleGain, analogfilter(Highpass(0.4), Elliptic(20, 0.5, 20)))
    zpkfilter_eq(f, m_f)
end

#
# Band pass filter construction
#
@testset "bandpass" begin
    # Output of
    #=
        [z, p, k] = butter(10, [0.25 0.75], 's')
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    m_z = zeros(ComplexF64, 10)
    m_p = [-0.05852661385905927+0.7442329445709392im,
        -0.05852661385905927-0.7442329445709392im,
        -0.1665549471789192+0.6992510095482473im,
        -0.1665549471789192-0.6992510095482473im,
        -0.2479888914449238+0.615606020752747im,
        -0.2479888914449238-0.615606020752747im,
        -0.2869902139927978+0.5070583156429217im,
        -0.2869902139927978-0.5070583156429217im,
        -0.2738314211701858+0.3979722116506606im,
        -0.2738314211701858-0.3979722116506606im,
        -0.2200127491273796+0.3197549791305416im,
        -0.2200127491273796-0.3197549791305416im,
        -0.1585130481013887+0.280063065773154im,
        -0.1585130481013887-0.280063065773154im,
        -0.1055644991483508+0.2620526301594698im,
        -0.1055644991483508-0.2620526301594698im,
        -0.06044030269085397+0.2537477474540631im,
        -0.06044030269085397-0.2537477474540631im,
        -0.01969061866105566+0.2503887742733703im,
        -0.01969061866105566-0.2503887742733703im]
    m_k = 0.0009765624999999997
    m_f = ZeroPoleGain(m_z, m_p, m_k)

    f = analogfilter(Bandpass(0.25, 0.75), Butterworth(10))
    zpkfilter_eq(f, m_f)

    # Test that our answers are more accurate than MATLAB's
    accurate_ft = analogfilter(Bandpass(0.25, 0.75), Butterworth(BigFloat, 10))
    zpkfilter_eq(f, accurate_ft)
    zpkfilter_accuracy(f, m_f, accurate_ft)

    # Output of
    #=
        [z, p, k] = ellip(10, 0.5, 20, [0.2 0.8], 's')
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    m_z = [1.387778780781446e-16+1.370743857088904im,
        1.387778780781446e-16-1.370743857088904im,
        -5.204170427930421e-18+0.1167249440313361im,
        -5.204170427930421e-18-0.1167249440313361im,
        0+0.8445880526167882im, 0-0.8445880526167882im,
        -1.734723475976807e-18+0.8053338099984824im,
        -1.734723475976807e-18-0.8053338099984824im,
        -1.387778780781446e-17+0.1894414673571001im,
        -1.387778780781446e-17-0.1894414673571001im,
        0+0.1986753790956593im, 0-0.1986753790956593im,
        5.551115123125783e-17+0.800148411646084im,
        5.551115123125783e-17-0.800148411646084im,
        1.108054620280186e-16+0.8007109707175198im,
        1.108054620280186e-16-0.8007109707175198im,
        8.673617379884035e-18+0.199822415142662im,
        8.673617379884035e-18-0.199822415142662im,
        4.85722573273506e-17+0.1999629039703328im,
        4.85722573273506e-17-0.1999629039703328im]
    m_p = [-0.2019763467382846+0.6120628823384283im,
        -0.2019763467382846-0.6120628823384283im,
        -0.03721004055421873+0.7823115667177283im,
        -0.03721004055421873-0.7823115667177283im,
        -0.004810677763557675+0.7980130970318471im,
        -0.004810677763557675-0.7980130970318471im,
        -0.0006000695021327641+0.7998051371015736im,
        -0.0006000695021327641-0.7998051371015736im,
        -6.552554548894252e-05+0.8000220660203423im,
        -6.552554548894252e-05-0.8000220660203423im,
        -0.07779252812598642+0.235740074311183im,
        -0.07779252812598642-0.235740074311183im,
        -0.009705975806501373+0.2040604370920958im,
        -0.009705975806501373-0.2040604370920958im,
        -0.001208621817121024+0.2004906765377495im,
        -0.001208621817121024-0.2004906765377495im,
        -0.0001500903998086778+0.2000486149853388im,
        -0.0001500903998086778-0.2000486149853388im,
        -1.638048261950673e-05+0.1999944823054316im,
        -1.638048261950673e-05-0.1999944823054316im]
    m_k = 0.1000000000000025
    m_f = ZeroPoleGain(m_z, m_p, m_k)

    f = convert(ZeroPoleGain, analogfilter(Bandpass(0.2, 0.8), Elliptic(10, 0.5, 20)))
    zpkfilter_eq(f, m_f)
end

#
# Band stop filter construction
#
@testset "bandstop" begin
    # Output of
    #=
        [z, p, k] = butter(10, [0.25 0.75], 'stop', 's')
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    m_z = [0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im,
        0+0.4330127018922193im, -0-0.4330127018922193im]
    m_p = [-0.05852661385905998+0.7442329445709392im,
        -0.05852661385905998-0.7442329445709392im,
        -0.1665549471789191+0.6992510095482469im,
        -0.1665549471789191-0.6992510095482469im,
        -0.2479888914449218+0.6156060207527428im,
        -0.2479888914449218-0.6156060207527428im,
        -0.2869902139927988+0.5070583156429307im,
        -0.2869902139927988-0.5070583156429307im,
        -0.2738314211701937+0.3979722116506589im,
        -0.2738314211701937-0.3979722116506589im,
        -0.2200127491273721+0.3197549791305355im,
        -0.2200127491273721-0.3197549791305355im,
        -0.1585130481013898+0.2800630657731571im,
        -0.1585130481013898-0.2800630657731571im,
        -0.1055644991483515+0.2620526301594707im,
        -0.1055644991483515-0.2620526301594707im,
        -0.06044030269085383+0.2537477474540633im,
        -0.06044030269085383-0.2537477474540633im,
        -0.01969061866105538+0.2503887742733704im,
        -0.01969061866105538-0.2503887742733704im]
    m_k = 1
    m_f = ZeroPoleGain(m_z, m_p, m_k)

    f = analogfilter(Bandstop(0.25, 0.75), Butterworth(10))
    zpkfilter_eq(f, m_f)

    # Test that our answers are more accurate than MATLAB's
    accurate_ft = analogfilter(Bandstop(0.25, 0.75), Butterworth(BigFloat, 10))
    zpkfilter_eq(f, accurate_ft)
    zpkfilter_accuracy(f, m_f, accurate_ft)

    #=
        [z, p, k] = ellip(10, 0.5, 20, [0.2 0.8], 'stop', 's')
        sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
        sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
        sprintf('%.16g', k)
    =#
    m_z = [-1.402517154063671e-17+0.5685129776423937im,
        -1.402517154063671e-17-0.5685129776423937im,
        -6.156235460644255e-18+0.2814359676774932im,
        -6.156235460644255e-18-0.2814359676774932im,
        -2.775557561562891e-17+0.760016963542146im,
        -2.775557561562891e-17-0.760016963542146im,
        -8.326672684688674e-17+0.7947386855861924im,
        -8.326672684688674e-17-0.7947386855861924im,
        5.551115123125783e-17+0.799851645232653im,
        5.551115123125783e-17-0.799851645232653im,
        8.326672684688674e-17+0.7992903327101535im,
        8.326672684688674e-17-0.7992903327101535im,
        0+0.2105216168522108im, 0-0.2105216168522108im,
        -4.423544863740858e-17+0.2013240363176547im,
        -4.423544863740858e-17-0.2013240363176547im,
        6.938893903907228e-18+0.200037095570968im,
        6.938893903907228e-18-0.200037095570968im,
        5.551115123125783e-17+0.2001775743458425im,
        5.551115123125783e-17-0.2001775743458425im]
    m_p = [-0.3778783653696471+0.7819846151027801im,
        -0.3778783653696471-0.7819846151027801im,
        -0.04044889726976308+0.8144610469677399im,
        -0.04044889726976308-0.8144610469677399im,
        -0.004859687206700769+0.8019364482578354im,
        -0.004859687206700769-0.8019364482578354im,
        -0.0006006726587018951+0.8001940294215477im,
        -0.0006006726587018951-0.8001940294215477im,
        -6.551807382872143e-05+0.799977924147539im,
        -6.551807382872143e-05-0.799977924147539im,
        -0.08015534922525629+0.1658741427311564im,
        -0.08015534922525629-0.1658741427311564im,
        -0.009732315617052817+0.1959655887284332im,
        -0.009732315617052817-0.1959655887284332im,
        -0.00120901710644529+0.1995097303565812im,
        -0.00120901710644529-0.1995097303565812im,
        -0.0001500952639848091+0.1999513917362522im,
        -0.0001500952639848091-0.1999513917362522im,
        -1.638042236460496e-05+0.2000055177738602im,
        -1.638042236460496e-05-0.2000055177738602im]
    m_k = 0.9440608762859234
    m_f = ZeroPoleGain(m_z, m_p, m_k)

    f = analogfilter(Bandstop(0.2, 0.8), Elliptic(10, 0.5, 20))
    zpkfilter_eq(f, m_f)
end

#
# IIR digital filter creation
#
@testset "digital IIR" begin
    # Output of [b, a] = butter(20, 0.5)
    m_b = [8.555155699386467e-06, 0.0001711031139877293,
        0.001625479582883429, 0.009752877497300572,
        0.04144972936352743, 0.1326391339632878,
        0.3315978349082195, 0.663195669816439,
        1.077692963451713, 1.436923951268951,
        1.580616346395846, 1.436923951268951,
        1.077692963451713, 0.663195669816439,
        0.3315978349082195, 0.1326391339632878,
        0.04144972936352743, 0.009752877497300572,
        0.001625479582883429, 0.0001711031139877293,
        8.555155699386467e-06]
    m_a = [1, -9.216626485446244e-15, 2.719352210996272,
        -2.276898740746355e-14, 2.954502383565172,
        -2.181204855972978e-14, 1.663404391798304,
        -9.886771142150971e-15, 0.5272563634534938,
        -1.578399203426284e-15, 0.09588816575455775,
        3.201156931050434e-16, 0.009788250615318749,
        1.651392612127069e-16, 0.0005259945756568433,
        3.034673947897289e-17, 1.306795536212839e-05,
        3.205637858266449e-18, 1.137794123801438e-07,
        9.233765069631315e-20, 1.463813780864114e-10]
    m_f = PolynomialRatio(m_b, m_a)

    f = convert(PolynomialRatio, digitalfilter(Lowpass(0.5), Butterworth(20)))
    tffilter_eq(f, m_f)

    # Output of [b, a] = butter(20, 0.3, 'high')
    m_b = [0.001797932208654781, -0.03595864417309561, 0.3416071196444083,
        -2.04964271786645, 8.710981550932413, -27.87514096298372,
        69.6878524074593, -139.3757048149186, 226.4855203242427,
        -301.980693765657, 332.1787631422226, -301.980693765657,
        226.4855203242427, -139.3757048149186, 69.6878524074593,
        -27.87514096298372, 8.710981550932413, -2.04964271786645,
        0.3416071196444083, -0.03595864417309561, 0.001797932208654781]
    m_a = [1, -7.99378842342338, 32.24776489137069, -86.15266953725099,
        169.613151737271, -260.1393113842562, 321.1838937487486,
        -325.8437287835835, 275.1584488419961, -194.8984651152305,
        116.2180984348268, -58.35420576527059, 24.59604119108612,
        -8.645276014889031, 2.506713994242244, -0.5898102957633149,
        0.1098868139124379, -0.01561185396257133, 0.001590100744530764,
        -0.0001034620065564712, 3.232560226924501e-06]
    m_f = PolynomialRatio(m_b, m_a)

    f = convert(PolynomialRatio, digitalfilter(Highpass(0.3), Butterworth(20)))
    tffilter_eq(f, m_f)

    # Output of [b, a] = butter(7, [0.1 0.4])
    m_b = [0.0009628940476886417, 0, -0.006740258333820491, 0,
        0.02022077500146148, 0, -0.03370129166910246,
        0, 0.03370129166910246, 0, -0.02022077500146148,
        0, 0.006740258333820491, 0, -0.0009628940476886417]
    m_a = [1, -7.763446908073567, 29.02150968091795, -69.51575164764138,
        119.5092528789707, -156.1341322158939, 159.8101214681526,
        -130.0847734566226, 84.56503905088692, -43.67462798637182,
        17.64013826588496, -5.405682857208016, 1.18954580728982,
        -0.1686342472129378, 0.01166272804918956]
    m_f = PolynomialRatio(m_b, m_a)

    f = convert(PolynomialRatio, digitalfilter(Bandpass(0.1, 0.4), Butterworth(7)))
    tffilter_eq(f, m_f)

    # Output of [b, a] = butter(4, [0.2 0.25], 'stop')
    m_b = [0.8142545568862771, -4.968628804996451, 14.62659178092118,
        -26.46885325790832, 32.03454308505385, -26.46885325790833,
        14.62659178092118, -4.968628804996451, 0.8142545568862771]
    m_a = [1, -5.789124865012897, 16.16849931613174, -27.77244867328942,
        31.92014561624699, -25.06018377653858, 13.16458034390183,
        -4.253206810966346, 0.6630104843858931]
    m_f = PolynomialRatio(m_b, m_a)

    f = convert(PolynomialRatio, digitalfilter(Bandstop(0.2, 0.25), Butterworth(4)))
    tffilter_eq(f, m_f)

    # Output of [b, a] = ellip(10, 0.7, 13, 0.4)
    m_b = [0.2797942479261318, -0.546743606764335, 1.704005589668604,
        -2.159911306584612, 3.642222926761125, -3.210370872716824,
        3.642222926761128, -2.159911306584614, 1.704005589668607,
        -0.5467436067643362, 0.2797942479261327]
    m_a = [1, -3.029857317356546, 7.982710263075236, -12.64912738840943,
        17.41532979440429, -17.06917827530357, 14.5191214973975,
        -8.798609366466101, 4.541539167299845, -1.419553904214128, 0.3565809079146795]
    m_f = PolynomialRatio(m_b, m_a)

    f = convert(PolynomialRatio, digitalfilter(Lowpass(0.4), Elliptic(10, 0.7, 13)))
    tffilter_eq(f, m_f)

    # Output of [b, a] = ellip(11, 0.7, 17, 0.2, 'high')
    m_b = [0.5737867105773529, -5.288884465999311, 23.07976318423696,
        -62.6831981311594, 117.4527496444787, -159.209133733912,
        159.2091337339121, -117.4527496444788, 62.68319813115945,
        -23.079763184237, 5.288884465999319, -0.5737867105773541]
    m_a = [1, -8.215565996888206, 32.26875224082856, -79.26375768096621,
        134.6328215735552, -165.4634116073934, 149.7532535390064,
        -99.57881504158784, 47.56023313164873, -15.48547977535943,
        3.075141867250689, -0.2777992862430019]
    m_f = PolynomialRatio(m_b, m_a)

    f = convert(PolynomialRatio, digitalfilter(Highpass(0.2), Elliptic(11, 0.7, 17)))
    tffilter_eq(f, m_f)

    # Output of [b, a] = ellip(10, 0.5, 2, [0.2 0.3])
    m_b = [0.7571867287399742, -10.55194327196835, 73.54073467886568,
        -338.6101060102045, 1150.502868412841, -3057.625219554799,
        6581.673203566225, -11732.90330745444, 17574.06211811366,
        -22320.04109078996, 24158.39342879073, -22320.04109078995,
        17574.06211811365, -11732.90330745443, 6581.673203566216,
        -3057.625219554793, 1150.502868412839, -338.6101060102037,
        73.54073467886548, -10.55194327196831, 0.7571867287399715]
    m_a = [1, -13.88122062422027, 96.34606972585827, -441.7050045576331,
        1494.048255114342, -3952.085137071483, 8465.769048060451,
        -15015.81246941687, 22374.64201809294, -28265.00095948501,
        30424.52053126076, -27950.27646663918, 21879.36223420588,
        -14520.36940748522, 8095.785094725445, -3737.670491911818,
        1397.479471832041, -408.6480215434079, 88.17084042962101,
        -12.56723162073412, 0.8957646298422474]
    m_f = PolynomialRatio(m_b, m_a)

    f = convert(PolynomialRatio, digitalfilter(Bandpass(0.2, 0.3), Elliptic(10, 0.5, 2)))
    tffilter_eq(f, m_f)

    # Output of [b, a] = ellip(15, 1.8, 52, [0.6 0.7], 'stop')
    m_b = [0.4754470525338341, 6.428928219976483, 47.47972043083799,
        245.8985324356427, 989.107195159456, 3259.995646867947,
        9099.649115772192, 21987.91200123768, 46709.48715187636,
        88217.96445932264, 149368.0297431405, 228141.9056208539,
        315792.8309895032, 397458.2884412499, 455878.980369108,
        477133.9777299722, 455878.9803691082, 397458.2884412502,
        315792.8309895035, 228141.9056208542, 149368.0297431408,
        88217.96445932281, 46709.48715187647, 21987.91200123774,
        9099.649115772219, 3259.995646867959, 989.1071951594597,
        245.8985324356436, 47.4797204308382, 6.428928219976513,
        0.475447052533836]
    m_a = [1, 12.89634891736776, 90.77856214451822, 447.9765044361528,
        1716.574771803699, 5388.380608104304, 14321.19101752449,
        32940.21064000877, 66586.59389679749, 119619.3933180238,
        192555.4138276511, 279453.5193544072, 367301.6888247615,
        438620.3975638363, 476900.358709065, 472649.3359368506,
        427104.9678341324, 351676.250633706, 263451.4609235193,
        179107.598840823, 110104.024869051, 60897.74794125407,
        30103.22741782092, 13182.01328029648, 5052.946545447418,
        1668.118080642225, 463.5003850459777, 104.731419523004,
        18.20949251383043, 2.194519514615899, 0.1423868039113367]
    m_f = PolynomialRatio(m_b, m_a)

    f = convert(PolynomialRatio, digitalfilter(Bandstop(0.6, 0.7), Elliptic(15, 1.8, 52)))
    tffilter_eq(f, m_f)
end

#
# Special filters
#
@testset "special" begin
    # Output of [num, den] = iirnotch(0.1, 0.1)
    notch = iirnotch(0.1, 0.1)
    @test notch.b0 ≈ 0.863271264002681
    @test notch.b1 ≈ -1.642039521920206
    @test notch.b2 ≈ 0.863271264002681
    @test notch.a1 ≈ -1.642039521920206
    @test notch.a2 ≈ 0.726542528005361

    notch = iirnotch(60, 1, fs=1000)
    @test notch.b0 ≈ 0.996868235770807
    @test notch.b1 ≈ -1.853729290297204
    @test notch.b2 ≈ 0.996868235770807
    @test notch.a1 ≈ -1.853729290297204
    @test notch.a2 ≈ 0.993736471541615
end

#
# IIR digital filter error conditions
#
@testset "errors" begin
    # zero pole
    @test_throws ErrorException Butterworth(0)
    @test_throws ErrorException Chebyshev1(0, 1)
    @test_throws ErrorException Chebyshev2(0, 1)
    @test_throws ErrorException Elliptic(0, 0.1, 0.5)

    # too much passband ripple
    @test_throws ErrorException Elliptic(2, 0.5, 0.1)

    # negative ripple
    @test_throws ErrorException Chebyshev1(2, -1)
    @test_throws ErrorException Chebyshev2(2, -1)
    @test_throws ErrorException Elliptic(2, 0.1, -0.5)
end

#
# Specification of sample rate
#
@testset "sample rate" begin
    @test Lowpass(20, fs=100) == Lowpass(0.4)
    @test Highpass(20, fs=100) == Highpass(0.4)
    @test Bandpass(20, 21, fs=100) == Bandpass(0.4, 0.42)
    @test Bandstop(20, 21, fs=100) == Bandstop(0.4, 0.42)
end

#
# kaiserord
#
@testset "kaiserord" begin
    kaiserord_jl    = kaiserord(0.1, 60)
    kaiserord_scipy = (74, 1.7994885471673767)    # (n,beta)=kaiserord(60,.1); (n,beta/pi)
    @test kaiserord_jl[1] == kaiserord_scipy[1]
    @test kaiserord_jl[2] ≈ kaiserord_scipy[2]

    kaiserord_jl    = kaiserord(0.1, 21.0 + eps(21.0))
    kaiserord_scipy = (20, 3.087704731700782e-07) # (n,beta)=kaiserord(21.000000000000004,.1); (n,beta/pi)
    @test kaiserord_jl[1] == kaiserord_scipy[1]
    @test kaiserord_jl[2] ≈ kaiserord_scipy[2]

    kaiserord_jl    = kaiserord(0.1, 15)
    kaiserord_scipy = (11, 0.0)                   # (n,beta)=kaiserord(15,.1); (n,beta/pi)
    @test kaiserord_jl[1] == kaiserord_scipy[1]
    @test kaiserord_jl[2] ≈ kaiserord_scipy[2]
end

#
# Window FIR filter taps
#
@testset "window FIR" begin
    winfirtaps_jl    = digitalfilter(Lowpass(0.25; fs=1),FIRWindow(hamming(128), scale=false))
    # firwin(128, 0.25, nyq=.5, scale=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_128_lowpass_fc0.25_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    winfirtaps_jl    = digitalfilter(Lowpass(0.25; fs=1),FIRWindow(hamming(129), scale=false))
    # firwin(129, 0.25, nyq=.5, scale=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_lowpass_fc0.25_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    @test_throws ArgumentError digitalfilter(Highpass(0.25; fs=1),FIRWindow(hamming(128), scale=false))

    winfirtaps_jl    = digitalfilter(Highpass(0.25; fs=1),FIRWindow(hamming(129), scale=false))
    # firwin(129, 0.25, nyq=.5, scale=False, pass_zero=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_highpass_fc0.25_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    winfirtaps_jl    = digitalfilter(Bandpass(0.1, 0.2; fs=1),FIRWindow(hamming(128), scale=false))
    # firwin(128, [0.1, 0.2], nyq=.5, scale=False, pass_zero=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_128_bandpass_fc0.1_0.2_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    winfirtaps_jl    = digitalfilter(Bandpass(0.1, 0.2; fs=1),FIRWindow(hamming(129), scale=false))
    # firwin(129, [0.1, 0.2], nyq=.5, scale=False, pass_zero=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_bandpass_fc0.1_0.2_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    @test_throws ArgumentError digitalfilter(Bandstop(0.1, 0.2; fs=1),FIRWindow(hamming(128), scale=false))

    winfirtaps_jl    = digitalfilter(Bandstop(0.1, 0.2; fs=1),FIRWindow(hamming(129), scale=false))
    # firwin(129, [0.1, 0.2], nyq=.5, scale=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_bandstop_fc0.1_0.2_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ vec(winfirtaps_scipy)

    winfirtaps_jl    = digitalfilter(Lowpass(0.25; fs=1),FIRWindow(hamming(128), scale=true))
    # firwin(128, 0.25, nyq=.5)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_128_lowpass_scaled_fc0.25_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    winfirtaps_jl    = digitalfilter(Lowpass(0.25; fs=1),FIRWindow(hamming(129), scale=true))
    # firwin(129, 0.25, nyq=.5)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_lowpass_scaled_fc0.25_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    winfirtaps_jl    = digitalfilter(Highpass(0.25; fs=1),FIRWindow(hamming(129), scale=true))
    # firwin(129, 0.25, nyq=.5, pass_zero=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_highpass_scaled_fc0.25_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    winfirtaps_jl    = digitalfilter(Bandpass(0.1, 0.2; fs=1),FIRWindow(hamming(128), scale=true))
    # firwin(128, [0.1, 0.2], nyq=.5, pass_zero=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_128_bandpass_scaled_fc0.1_0.2_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    winfirtaps_jl    = digitalfilter(Bandpass(0.1, 0.2; fs=1),FIRWindow(hamming(129), scale=true))
    # firwin(129, [0.1, 0.2], nyq=.5, scale=False, pass_zero=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_bandpass_scaled_fc0.1_0.2_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ winfirtaps_scipy

    winfirtaps_jl    = digitalfilter(Bandstop(0.1, 0.2; fs=1),FIRWindow(hamming(129), scale=true))
    # firwin(129, [0.1, 0.2], nyq=.5, scale=False)
    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_bandstop_scaled_fc0.1_0.2_fs1.0.txt"),'\t')
    @test winfirtaps_jl ≈ vec(winfirtaps_scipy)
end

# make sure typing is preserved
# ie, user passes in F32, should get ZPK{CF32, CF32, F32}
@testset "typing" begin
    @testset "prototypes" begin
        @test_throws MethodError Butterworth(ComplexF64, 4)
        @test_throws MethodError Chebyshev1(ComplexF64, 4, 1.01)
        @test_throws MethodError Chebyshev2(ComplexF64, 4, 1.01)
        @test_throws MethodError Elliptic(ComplexF64, 4, 0.5, 1.01)

        @test_throws MethodError Butterworth(Char, 4)
        @test_throws MethodError Chebyshev1(Char, 4, 1.01)
        @test_throws MethodError Chebyshev2(Char, 4, 1.01)
        @test_throws MethodError Elliptic(Char, 4, 0.5, 1.01)

        @testset "$T" for T in [BigFloat, Float16, Float32, Float64]
            @test Butterworth(T, 4) isa ZeroPoleGain{:s, Complex{T}, Complex{T}, T}
            @test Chebyshev1(T, 5, 1.16) isa ZeroPoleGain{:s, Complex{T}, Complex{T}, T}
            @test Chebyshev2(T, 5, 1.16) isa ZeroPoleGain{:s, Complex{T}, Complex{T}, T}
            @test Elliptic(T, 5, 0.31, 1.16) isa ZeroPoleGain{:s, Complex{T}, Complex{T}, T}
        end
    end

    @testset "transformed types" begin
        @testset "$T & $R" for T in [BigFloat, Float16, Float32, Float64],
                R in [Lowpass(0.4), Highpass(0.6), Bandpass(0.3, 0.6), Bandstop(0.7, 0.8)]
            @test digitalfilter(R, Butterworth(T, 5)) isa ZeroPoleGain{:z, Complex{T}, Complex{T}, T}
            @test digitalfilter(R, Chebyshev1(T, 5, 1.16)) isa ZeroPoleGain{:z, Complex{T}, Complex{T}, T}
            @test digitalfilter(R, Chebyshev2(T, 5, 1.16)) isa ZeroPoleGain{:z, Complex{T}, Complex{T}, T}
            @test digitalfilter(R, Elliptic(T, 5, 0.31, 1.16)) isa ZeroPoleGain{:z, Complex{T}, Complex{T}, T}
        end
    end
end


text:
using Hecke

using Base.Test

for scope in Hecke.ASSERT_SCOPE
  set_assert_level(scope, 3)
end

include("AssAlg.jl")
include("EllCrv.jl")
include("GrpAb.jl")
include("LinearAlgebra.jl")
include("Map.jl")
include("Misc.jl")
include("NfAbs.jl")
include("NfOrd.jl")
include("NfRel.jl")
include("RCF.jl")

# x^5 + 514944*x^2 + 123904 test prime decomposition with this (2 is index divisor and only one prime ideal over 2)


text:
using Statistics

function board(X_0::Array, b::Float64)

    x0 = copy(X_0[1])
    y0 = copy(X_0[2])
    u0 = copy(X_0[3])
    v0 = copy(X_0[4])

    g = 9.81
    c = 0.5
    a = 0.25
    p = (2*π)/10.0
    q = (2*π)/4.0

    H = -a*x0 - b*cos(p*x0)*cos(q*y0)
    H_x = -a + b*p*sin(p*x0)*cos(q*y0)
    H_xx = b*p^2 * cos(p*x0)*cos(q*y0)
    H_y = b*q*cos(p*x0)*sin(q*y0)
    H_yy = b*q^2 * cos(p*x0)*cos(q*y0)
    H_xy = -b*q*p*sin(p*x0)*sin(q*y0)

    F = (g + H_xx*u0^2 + 2*H_xy*u0*v0 + H_yy*v0^2)/(1 + H_x^2 + H_y^2)

    dU = -F*H_x - c*u0
    dV = -F*H_y - c*v0

    return [u0, v0, dU, dV]

end

function runge_kutta_step(f::Function, x0, dt::Float64, t=nothing)

    k1 = f(x0, t) * dt
    k2 = f(x0 + k1/2, t) * dt
    k3 = f(x0 + k2/2, t) * dt
    k4 = f(x0 + k3, t) * dt
    x_new = x0 + (k1 + 2*k2 + 2*k3 + k4)/6

    return x_new
end

function solver(f::Function, x0::Array, y0::Array, u0::Array, v0::Array, dt::Float64, N_t::Int, b = 0.5)
    solutions = []
    for k in 1:length(x0)
        values_one_step = [x0[k], y0[k], u0[k], v0[k]]
        for i in 2:N_t + 1
            values_one_step = runge_kutta_step(f, values_one_step, dt, b)
        end
        solutions = vcat(solutions, values_one_step)
    end

    return solutions
end

function bench(n_sleds, n_time)
    x_init = zeros(n_sleds)
    y_init = rand(n_sleds)
    u_init = zeros(n_sleds)
    v_init = zeros(n_sleds) .+ 3.5

    return solver(board, x_init, y_init, u_init, v_init, 0.01, n_time)
end

n_sleds = 10
n_time = 200

nb_runs = 200

times = zeros(nb_runs)

for irun in 1:nb_runs
    times[irun] = @elapsed bench(n_sleds, n_time)
end

println(median(times))

max_stars_repo_name: fluiddyn/piconumpy

max_issues_repo_name: fluiddyn/piconumpy

max_forks_repo_name: fluiddyn/piconumpy


text:
include("opfdata.jl")
using JuMP
using Ipopt

function acopf_solve(opfmodel, opf_data)

  #
  # Initial point - needed especially for pegase cases
  #
  Pg0,Qg0,Vm0,Va0 = acopf_initialPt_IPOPT(opf_data)
  setvalue(getvariable(opfmodel, :Pg), Pg0)
  setvalue(getvariable(opfmodel, :Qg), Qg0)
  setvalue(getvariable(opfmodel, :Vm), Vm0)
  setvalue(getvariable(opfmodel, :Va), Va0)

  status = solve(opfmodel)

  if status != :Optimal
    println("Could not solve the model to optimality.")
  end
  return opfmodel,status
end

function acopf_model(opf_data)
  #shortcuts for compactness
  lines = opf_data.lines; buses = opf_data.buses; generators = opf_data.generators; baseMVA = opf_data.baseMVA
  busIdx = opf_data.BusIdx; FromLines = opf_data.FromLines; ToLines = opf_data.ToLines; BusGeners = opf_data.BusGenerators;

  nbus  = length(buses); nline = length(lines); ngen  = length(generators)

  #branch admitances
  YffR,YffI,YttR,YttI,YftR,YftI,YtfR,YtfI,YshR,YshI = computeAdmitances(lines, buses, baseMVA)

  #
  # JuMP model now
  #
  opfmodel = Model(solver=IpoptSolver(print_level=0))

  @variable(opfmodel, generators[i].Pmin <= Pg[i=1:ngen] <= generators[i].Pmax)
  @variable(opfmodel, generators[i].Qmin <= Qg[i=1:ngen] <= generators[i].Qmax)

  @variable(opfmodel, buses[i].Vmin <= Vm[i=1:nbus] <= buses[i].Vmax)
  @variable(opfmodel, Va[1:nbus])
  #fix the voltage angle at the reference bus
  setlowerbound(Va[opf_data.bus_ref], buses[opf_data.bus_ref].Va)
  setupperbound(Va[opf_data.bus_ref], buses[opf_data.bus_ref].Va)

  # minimize active power
#  @NLobjective(opfmodel,
#		  Min,
#		  sum{ generators[i].coeff[generators[i].n] +
#		       sum{generators[i].coeff[generators[i].n-k]*(baseMVA*Pg[i])^k, k=1:generators[i].n-1},
#		       i=1:ngen}
#		 )

  @NLobjective(opfmodel, Min, sum{ generators[i].coeff[generators[i].n-2]*(baseMVA*Pg[i])^2
			             +generators[i].coeff[generators[i].n-1]*(baseMVA*Pg[i])
				     +generators[i].coeff[generators[i].n  ], i=1:ngen})

  #
  # power flow balance
  #
  for b in 1:nbus
    #real part
    @NLconstraint(
      opfmodel,
      ( sum{ YffR[l], l in FromLines[b]} + sum{ YttR[l], l in ToLines[b]} + YshR[b] ) * Vm[b]^2
      + sum{ Vm[b]*Vm[busIdx[lines[l].to]]  *( YftR[l]*cos(Va[b]-Va[busIdx[lines[l].to]]  ) + YftI[l]*sin(Va[b]-Va[busIdx[lines[l].to]]  )), l in FromLines[b] }
      + sum{ Vm[b]*Vm[busIdx[lines[l].from]]*( YtfR[l]*cos(Va[b]-Va[busIdx[lines[l].from]]) + YtfI[l]*sin(Va[b]-Va[busIdx[lines[l].from]])), l in ToLines[b]   }
      - ( sum{baseMVA*Pg[g], g in BusGeners[b]} - buses[b].Pd ) / baseMVA      # Sbus part
      ==0)
#  end
#  for b in 1:nbus
    #imaginary part
    @NLconstraint(
      opfmodel,
      ( sum{-YffI[l], l in FromLines[b]} + sum{-YttI[l], l in ToLines[b]} - YshI[b] ) * Vm[b]^2
      + sum{ Vm[b]*Vm[busIdx[lines[l].to]]  *(-YftI[l]*cos(Va[b]-Va[busIdx[lines[l].to]]  ) + YftR[l]*sin(Va[b]-Va[busIdx[lines[l].to]]  )), l in FromLines[b] }
      + sum{ Vm[b]*Vm[busIdx[lines[l].from]]*(-YtfI[l]*cos(Va[b]-Va[busIdx[lines[l].from]]) + YtfR[l]*sin(Va[b]-Va[busIdx[lines[l].from]])), l in ToLines[b]   }
      - ( sum{baseMVA*Qg[g], g in BusGeners[b]} - buses[b].Qd ) / baseMVA      #Sbus part
      ==0)
  end
  #
  # branch/lines flow limits
  #
  nlinelim=0
  for l in 1:nline
    if lines[l].rateA!=0 && lines[l].rateA<1.0e10
      nlinelim += 1
      flowmax=(lines[l].rateA/baseMVA)^2

      #branch apparent power limits (from bus)
      Yff_abs2=YffR[l]^2+YffI[l]^2; Yft_abs2=YftR[l]^2+YftI[l]^2
      Yre=YffR[l]*YftR[l]+YffI[l]*YftI[l]; Yim=-YffR[l]*YftI[l]+YffI[l]*YftR[l]
      @NLconstraint(
        opfmodel,
	Vm[busIdx[lines[l].from]]^2 *
	( Yff_abs2*Vm[busIdx[lines[l].from]]^2 + Yft_abs2*Vm[busIdx[lines[l].to]]^2
	  + 2*Vm[busIdx[lines[l].from]]*Vm[busIdx[lines[l].to]]*(Yre*cos(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])-Yim*sin(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]]))
	)
        - flowmax <=0)
#    end
#  end
#  for l in 1:nline
#    if lines[l].rateA!=0 && lines[l].rateA<1.0e10
#
#      flowmax=(lines[l].rateA/baseMVA)^2

      #branch apparent power limits (to bus)
      Ytf_abs2=YtfR[l]^2+YtfI[l]^2; Ytt_abs2=YttR[l]^2+YttI[l]^2
      Yre=YtfR[l]*YttR[l]+YtfI[l]*YttI[l]; Yim=-YtfR[l]*YttI[l]+YtfI[l]*YttR[l]
      @NLconstraint(
        opfmodel,
	Vm[busIdx[lines[l].to]]^2 *
        ( Ytf_abs2*Vm[busIdx[lines[l].from]]^2 + Ytt_abs2*Vm[busIdx[lines[l].to]]^2
          + 2*Vm[busIdx[lines[l].from]]*Vm[busIdx[lines[l].to]]*(Yre*cos(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])-Yim*sin(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]]))
        )
        - flowmax <=0)
    end
  end

  @printf("Buses: %d  Lines: %d  Generators: %d\n", nbus, nline, ngen)
  println("Lines with limits  ", nlinelim)

  return opfmodel
end

  #######################################################

  #values = zeros(2*nbus+2*ngen)
  ## values[1:2*nbus+2*ngen] = readdlm("/sandbox/petra/work/installs/matpower5.1/vars2.txt")
  #values[1:2*nbus+2*ngen] = readdlm("/sandbox/petra/work/installs/matpower5.1/vars3_9241.txt")
  #d = JuMP.NLPEvaluator(opfmodel)
  #MathProgBase.initialize(d, [:Jac])

  #g = zeros(2*nbus+2*nlinelim)
  #MathProgBase.eval_g(d, g, values)
  #println("f=", MathProgBase.eval_f(d,values))

  #gmat=zeros(2*nbus+2*nlinelim)
  #gmat[1:end] = readdlm("/sandbox/petra/work/installs/matpower5.1/cons3_9241.txt")
  #println("diff: ", norm(gmat-g))

  #println(opfmodel)

  #############################################################

function acopf_outputAll(opfmodel, opf_data)
  #shortcuts for compactness
  lines = opf_data.lines; buses = opf_data.buses; generators = opf_data.generators; baseMVA = opf_data.baseMVA
  busIdx = opf_data.BusIdx; FromLines = opf_data.FromLines; ToLines = opf_data.ToLines; BusGeners = opf_data.BusGenerators;

  nbus  = length(buses); nline = length(lines); ngen  = length(generators)

  # OUTPUTING
  println("Objective value: ", getobjectivevalue(opfmodel), "USD/hr")
  VM=getvalue(getvariable(opfmodel,:Vm)); VA=getvalue(getvariable(opfmodel,:Va))
  PG=getvalue(getvariable(opfmodel,:Pg)); QG=getvalue(getvariable(opfmodel,:Qg))

  println("============================= BUSES ==================================")
  println("  BUS    Vm     Va   |   Pg (MW)    Qg(MVAr) ")   # |    P (MW)     Q (MVAr)")  #|         (load)   ")

  println("                     |     (generation)      ")
  println("----------------------------------------------------------------------")
  for i in 1:nbus
    @printf("%4d | %6.2f  %6.2f | %s  | \n",
	    buses[i].bus_i, VM[i], VA[i]*180/pi,
	    length(BusGeners[i])==0?"   --          --  ":@sprintf("%7.2f     %7.2f", baseMVA*PG[BusGeners[i][1]], baseMVA*QG[BusGeners[i][1]]))
  end
  println("\n")

  within=20 # percentage close to the limits

  nflowlim=0
  for l in 1:nline
    if lines[l].rateA!=0 && lines[l].rateA<1.0e10
      nflowlim += 1
    end
  end

  if nflowlim>0
    println("Number of lines with flow limits: ", nflowlim)

    optvec=zeros(2*nbus+2*ngen)
    optvec[1:ngen]=PG
    optvec[ngen+1:2*ngen]=QG
    optvec[2*ngen+1:2*ngen+nbus]=VM
    optvec[2*ngen+nbus+1:2*ngen+2*nbus]=VA

    d = JuMP.NLPEvaluator(opfmodel)
    MathProgBase.initialize(d, [:Jac])

    consRhs = zeros(2*nbus+2*nflowlim)
    MathProgBase.eval_g(d, consRhs, optvec)

    #println(consRhs)

    @printf("================ Lines within %d %s of flow capacity ===================\n", within, "\%")
    println("Line   From Bus    To Bus    At capacity")

    nlim=1
    for l in 1:nline
      if lines[l].rateA!=0 && lines[l].rateA<1.0e10
        flowmax=(lines[l].rateA/baseMVA)^2
        idx = 2*nbus+nlim

        if( (consRhs[idx]+flowmax)  >= (1-within/100)^2*flowmax )
          @printf("%3d      %3d      %3d        %5.3f%s\n", l, lines[l].from, lines[l].to, 100*sqrt((consRhs[idx]+flowmax)/flowmax), "\%" )
          #@printf("%7.4f   %7.4f    %7.4f \n", consRhs[idx], consRhs[idx]+flowmax,  flowmax)
        end
        nlim += 1
      end
    end
  end

  #println(getvalue(Vm))
  #println(getvalue(Va)*180/pi)

  #println(getvalue(Pg))
  #println(getvalue(Qg))

  return
end

# Compute initial point for IPOPT based on the values provided in the case data
function acopf_initialPt_IPOPT(opfdata)
  Pg=zeros(length(opfdata.generators)); Qg=zeros(length(opfdata.generators)); i=1
  for g in opfdata.generators
    # set the power levels in in between the bounds as suggested by matpower
    # (case data also contains initial values in .Pg and .Qg - not used with IPOPT)
    Pg[i]=0.5*(g.Pmax+g.Pmin)
    Qg[i]=0.5*(g.Qmax+g.Qmin)
    i=i+1
  end
  @assert i-1==length(opfdata.generators)

  Vm=zeros(length(opfdata.buses)); i=1;
  for b in opfdata.buses
    # set the ini val for voltage magnitude in between the bounds
    # (case data contains initials values in Vm - not used with IPOPT)
    Vm[i]=0.5*(b.Vmax+b.Vmin);
    i=i+1
  end
  @assert i-1==length(opfdata.buses)

  # set all angles to the angle of the reference bus
  Va = opfdata.buses[opfdata.bus_ref].Va * ones(length(opfdata.buses))

  return Pg,Qg,Vm,Va
end


text:
export RBTorsion

"""
    RBTorsion(; f1, f2, f3, f4)

A Ryckaert-Bellemans torsion angle between four atoms.
"""
struct RBTorsion{T} <: SpecificInteraction
    f1::T
    f2::T
    f3::T
    f4::T
end

RBTorsion(; f1, f2, f3, f4) = RBTorsion{typeof(f1)}(f1, f2, f3, f4)

@inline @inbounds function force(d::RBTorsion, coords_i, coords_j, coords_k,
                                    coords_l, box_size)
    ab = vector(coords_i, coords_j, box_size)
    bc = vector(coords_j, coords_k, box_size)
    cd = vector(coords_k, coords_l, box_size)
    cross_ab_bc = ab × bc
    cross_bc_cd = bc × cd
    bc_norm = norm(bc)
    θ = atan(
        ustrip(dot(cross_ab_bc × cross_bc_cd, bc / bc_norm)),
        ustrip(dot(cross_ab_bc, cross_bc_cd)),
    )
    dEdθ = (d.f1*sin(θ) - 2*d.f2*sin(2*θ) + 3*d.f3*sin(3*θ)) / 2
    fi =  dEdθ * bc_norm * cross_ab_bc / dot(cross_ab_bc, cross_ab_bc)
    fl = -dEdθ * bc_norm * cross_bc_cd / dot(cross_bc_cd, cross_bc_cd)
    v = (dot(-ab, bc) / bc_norm^2) * fi - (dot(-cd, bc) / bc_norm^2) * fl
    fj =  v - fi
    fk = -v - fl
    return SpecificForce4Atoms(fi, fj, fk, fl)
end

@inline @inbounds function potential_energy(d::RBTorsion, coords_i, coords_j, coords_k,
                                            coords_l, box_size)
    θ = torsion_angle(coords_i, coords_j, coords_k, coords_l, box_size)
    return (d.f1 * (1 + cos(θ)) + d.f2 * (1 - cos(2θ)) + d.f3 * (1 + cos(3θ)) + d.f4) / 2
end


text:
# ---
# title: 839. Similar String Groups
# id: problem839
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Depth-first Search, Union Find, Graph
# link: <https://leetcode.com/problems/similar-string-groups/description/>
# hidden: true
# ---
#
# Two strings `X` and `Y` are similar if we can swap two letters (in different
# positions) of `X`, so that it equals `Y`. Also two strings `X` and `Y` are
# similar if they are equal.
#
# For example, `"tars"` and `"rats"` are similar (swapping at positions `0` and
# `2`), and `"rats"` and `"arts"` are similar, but `"star"` is not similar to
# `"tars"`, `"rats"`, or `"arts"`.
#
# Together, these form two connected groups by similarity: `{"tars", "rats",
# "arts"}` and `{"star"}`.  Notice that `"tars"` and `"arts"` are in the same
# group even though they are not similar.  Formally, each group is such that a
# word is in the group if and only if it is similar to at least one other word
# in the group.
#
# We are given a list `strs` of strings where every string in `strs` is an
# anagram of every other string in `strs`. How many groups are there?
#
#
#
# **Example 1:**
#
#
#
#     Input: strs = ["tars","rats","arts","star"]
#     Output: 2
#
#
# **Example 2:**
#
#
#
#     Input: strs = ["omv","ovm"]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= strs.length <= 100`
#   * `1 <= strs[i].length <= 1000`
#   * `sum(strs[i].length) <= 2 * 104`
#   * `strs[i]` consists of lowercase letters only.
#   * All words in `strs` have the same length and are anagrams of each other.
#
#
## @lc code=start
using LeetCode

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


text:
import AngularMomentumAlgebra: @linearly_combinable

struct MySym
    s::Symbol
end
@linearly_combinable MySym

@testset "Linear combinations" begin
    x = MySym(:x)
    y = MySym(:y)

    w = x + y
    @test length(w) == 2
    @test eltype(w) == (MySym,Int)
    @test w.Ts == [x,y]
    @test w.coeffs == [1,1]

    @test (-x).coeffs == [-1]
    @test length((-x)) == 1

    z = 4x - 5y
    z′ = -z
    z′′ = 2z

    @test z.Ts == [x,y]
    @test z.coeffs == [4, -5]

    @test z′.Ts == [x,y]
    @test z′.coeffs == [-4, 5]

    @test z′′.Ts == [x,y]
    @test z′′.coeffs == [8, -10]

    @test all((z′′/2).coeffs .== z.coeffs)

    @test [(t,c) for (t,c) in z] == [(x,4), (y,-5)]
end


text:
using ALGCM1D
using Test
using Dates

@testset "ALGCM1D.jl" begin
    # Write your tests here.
    stime=DateTime(0,7,1)
    ALGCM1D.run(stime)
end


text:
"""
    ImageClassification(classes, sz[; augmentations, ...]) <: Method{ImageClassificationTask}
    ImageClassification(n, ...)

A [`LearningMethod`](#) for multi-class image classification using softmax probabilities.

`classes` is a vector of the category labels. Alternatively, you can pass an integer.
Images are resized to `sz`.

During training, a random crop is used and `augmentations`, a `DataAugmentation.Transform`
are applied.

### Types

- `input::AbstractMatrix{2, <:Colorant}`: an image
- `target` the class label that the image belongs to
- `x::AbstractArray{Float32, 3}`: a normalized 3D-array with dimensions *height, width, channels*
- `y::AbstractVector{Float32}`: one-hot encoding of category

### Model

- input size: `(sz..., ch, batch)` where `ch` depends on color type `C`.
- output size: `(nclasses, batch)`
"""
mutable struct ImageClassification <: DLPipelines.LearningMethod{ImageClassificationTask}
    sz::Tuple{Int, Int}
    classes::AbstractVector
    projectivetransforms::ProjectiveTransforms
    imagepreprocessing::ImagePreprocessing
end

Base.show(io::IO, method::ImageClassification) = print(
    io, "ImageClassification() with $(length(method.classes)) classes")

function ImageClassification(
        classes::AbstractVector,
        sz = (224, 224);
        augmentations = Identity(),
        means = IMAGENET_MEANS,
        stds = IMAGENET_STDS,
        C = RGB{N0f8},
        T = Float32
    )
    projectivetransforms = ProjectiveTransforms(sz, augmentations = augmentations)
    imagepreprocessing = ImagePreprocessing(means, stds; C = C, T = T)
    ImageClassification(sz, classes, projectivetransforms, imagepreprocessing)
end

ImageClassification(n::Int, args...; kwargs...) = ImageClassification(1:n, args...; kwargs...)

# Core interface implementation

function DLPipelines.encodeinput(
        method::ImageClassification,
        context,
        image)
    imagecropped = run(method.projectivetransforms, context, image)
    x = run(method.imagepreprocessing, context, imagecropped)
    return x
end

function DLPipelines.encodetarget(
        method::ImageClassification,
        context,
        category)
    idx = findfirst(isequal(category), method.classes)
    isnothing(idx) && error("`category` could not be found in `method.classes`.")
    return DataAugmentation.onehot(idx, length(method.classes))
end

function DLPipelines.encodetarget!(
        y::AbstractVector{T},
        method::ImageClassification,
        context,
        category) where T
    fill!(y, zero(T))
    idx = findfirst(isequal(category), method.classes)
    y[idx] = one(T)
    return y
end

DLPipelines.decodeŷ(method::ImageClassification, context, ŷ) = method.classes[argmax(ŷ)]

# Interpretation interface

DLPipelines.interpretinput(::ImageClassification, image) = image

function DLPipelines.interpretx(method::ImageClassification, x)
    return invert(method.imagepreprocessing, x)
end

function DLPipelines.interprettarget(task::ImageClassification, class)
    return "Class $class"
end

# Training interface

function DLPipelines.methodmodel(method::ImageClassification, backbone)
    h, w, ch, b = Flux.outdims(backbone, (method.sz..., 3, 1))
    return Chain(
        backbone,
        Chain(
            AdaptiveMeanPool((1,1)),
            flatten,
            Dense(ch, length(method.classes)),
        )
    )
end

DLPipelines.methodlossfn(::ImageClassification) = Flux.Losses.logitcrossentropy

# Testing interface

function DLPipelines.mockinput(method::ImageClassification)
    inputsz = rand.(UnitRange.(method.sz, method.sz .* 2))
    return rand(RGB{N0f8}, inputsz)
end

function DLPipelines.mocktarget(method::ImageClassification)
    rand(1:length(method.classes))
end

function DLPipelines.mockmodel(method::ImageClassification)
    return xs -> rand(Float32, length(method.classes), size(xs)[end])
end


text:
function alignCoef(coef, wavelet, coe, modwt, inverse, filter)
	(N, J, ser) = size(coef)
    for j=1:J
      shift = wtFilterShift(filter, j, wavelet, coe, modwt)

      shift >= N ? shift = shift - floor(shift/N)*N : shift = shift

      if shift == 0
        coef[:,j,:] = coef[:,j,:]
      end

      if inverse
      	coef[:,j,:] = coef[[(end-shift+1):end,1:(end-shift)],j,:]
      else
      	coef[:,j,:] = coef[[(shift+1):end,1:shift],j,:]
      end
    end
    return coef
end

function align(wt::WaveletTransform, coe::Boolean, inverse::Boolean)
  # Do the alignment only if neccessary
  @assume !inverse == wt.aligned

  # Check what kind of transform do I have.
  wt.filter.transform == "modwt" ? (modwt = true) : (modwt = false)

  # Make a copy of the original object
  wtShifted = deepcopy(wt)

  # Adjust the coefficients
  wtShifted.W = alignCoef(wt.W, true, coe, modwt, inverse, wt.filter)
  wtShifted.V = alignCoef(wt.V, false, coe, modwt, inverse, wt.filter)

  inverse ? wtShifted.aligned = false : wtShifted.aligned = true

  return wtShifted
end


text:
# sensors: want to map from sensor name to vector of data

# geom, body, joint, sensor
# site: xpos, etc
# textures, materials, camera

pm = mj_loadXML(modelfile, C_NULL)
pd = mj_makeData(pm)
m, d = mj.mapmujoco(pm, pd)

#sensors = Dict{Symbol, Range}(Symbol(snsr_names[i]) =>
#                              (m.sensor_adr[i]+1):(m.sensor_adr[i]+m.sensor_dim[i])
#                              for i=1:nsnsr)
sensors = mj.name2range(m, mj.get(m, :nsensor),
                        m.name_sensoradr, m.sensor_adr, m.sensor_dim)

sname = String(copy(m.names))   # See https://github.com/JuliaLang/julia/pull/26093
nsnsr = length(m.name_sensoradr)
idx = m.name_sensoradr[1] + 1
snsr_names = split(sname[idx:end], '\0', limit=(nsnsr+1))[1:nsnsr]
snsrs   = Dict{Symbol, SubArray{Float64}}(Symbol(snsr_names[i]) =>
                                          view(d.sensordata, (m.sensor_adr[i]+1):(m.sensor_adr[i]+m.sensor_dim[i]))
                                          for i=1:nsnsr)

bodies = mj.name2idx(m, mj.get(m, :nbody), m.name_bodyadr) # maps body name to it's index

steps = 1000

t = @elapsed begin
    for i=1:steps
        mj_step(m,d)
    end
end

@info "Time for $steps mj.steps: $t"

@test d.xquat[:, bodies[:world] ] == [1.0, 0.0, 0.0, 0.0]
@test isapprox( d.sensordata[sensors[:accel]][1], 8.524, rtol=1e-3)
@test isapprox( snsrs[:accel][1], 8.524, rtol=1e-3)


text:
using Distributions

include("datastats.jl")

"""
Base type for Naive Bayes models.
Inherited classes should have at least following fields:
 c_counts::Dict{C, Int64} - count of ocurrences of each class
 n_obs::Int64             - total number of observations
"""
abstract type NBModel{C} end

#####################################
#####  Multinomial Naive Bayes  #####
#####################################

mutable struct MultinomialNB{C} <: NBModel{C}
    c_counts::Dict{C, Int64}           # count of ocurrences of each class
    x_counts::Dict{C, Vector{Number}}  # count/sum of occurrences of each var
    x_totals::Vector{Number}           # total occurrences of each var
    n_obs::Int64                       # total number of seen observations
end

"""
Multinomial Naive Bayes classifier

classes : array of objects
    Class names
n_vars : Int64
    Number of variables in observations
alpha : Number (optional, default 1)
    Smoothing parameter. E.g. if alpha equals 1, each variable in each class
    is believed to have 1 observation by default
"""
function MultinomialNB(classes::Vector{C}, n_vars::Int64; alpha=1) where C
    c_counts = Dict(zip(classes, ones(Int64, length(classes)) * alpha))
    x_counts = Dict{C, Vector{Int64}}()
    for c in classes
        x_counts[c] = ones(Int64, n_vars) * alpha
    end
    x_totals = ones(Float64, n_vars) * alpha * length(c_counts)
    MultinomialNB{C}(c_counts, x_counts, x_totals, sum(x_totals))
end

function Base.show(io::IO, m::MultinomialNB)
    print(io, "MultinomialNB($(m.c_counts))")
end

#####################################
######  Gaussian Naive Bayes  #######
#####################################

mutable struct GaussianNB{C} <: NBModel{C}
    c_counts::Dict{C, Int64}           # count of ocurrences of each class
    c_stats::Dict{C, DataStats}        # aggregative data statistics
    gaussians::Dict{C, MvNormal}        # precomputed distribution
    # x_counts::Dict{C, Vector{Number}}  # ?? count/sum of occurrences of each var
    # x_totals::Vector{Number}           # ?? total occurrences of each var
    n_obs::Int64                       # total number of seen observations
end

function GaussianNB(classes::Vector{C}, n_vars::Int64) where C
    c_counts = Dict(zip(classes, zeros(Int64, length(classes))))
    c_stats = Dict(zip(classes, [DataStats(n_vars, 2) for i=1:length(classes)]))
    gaussians = Dict{C, MvNormal}()
    GaussianNB{C}(c_counts, c_stats, gaussians, 0)
end

function Base.show(io::IO, m::GaussianNB)
    print(io, "GaussianNB($(m.c_counts))")
end

#####################################
#####  Hybrid Naive Bayes       #####
#####################################
""" a wrapper around key value pairs for a discrete probability distribution """
struct ePDF{C <: AbstractDict}
    pairs::C
end

""" Constructor of ePDF """
function ePDF(x::AbstractVector{T}) where T <: Integer
    cnts = counts(x)
    ρ = map(Float64, cnts)/sum(cnts)
    ρ[ρ .< eps(Float64)] .= eps(Float64)
    d = Dict{Int, Float64}()
    for (k,v) in zip(StatsBase.span(x), ρ)
        d[k]=v
    end
    return ePDF(d)
end

""" query the ePDF to get the probability of n"""
function probability(P::ePDF, n::Integer)
    if n in keys(P.pairs)
        return P.pairs[n]
    else
        return eps(eltype(values(P.pairs)))
    end
end

"""
Initialize a `HybridNB` model with continuous and/or discrete features

### Constructors
```julia
HybridNB(labels::AbstractVector, kde_names::AbstractVector, discrete_names::AbstractVector)
HybridNB(labels::AbstractVector, kde_names::AbstractVector)
HybridNB(labels::AbstractVector, num_kde::Int, num_discrete::Int)
```
### Arguments
* `labels` : An AbstractVector{Any} of feature labels
* `kde_names` : An AbstractVector{Any} of the names of continuous features
* `discrete_names` : An AbstractVector{Any} of the names of discrete features
* `num_kde` : Number of continuous features
* `num_discrete` : Number of discrete features

"""
struct HybridNB{C <: Integer, N}
    c_kdes::Dict{C, Dict{N, InterpKDE}}
    c_discrete::Dict{C, Dict{N, ePDF}}
    classes::Vector{C}
    priors::Dict{C, Float64}
end

function num_features(m::HybridNB)
    length(m.classes) > 0 || throw("Number of kdes is not defined. There are no valid classes in the model.")
    c = m.classes[1]
    return length(m.c_kdes[c]), length(m.c_discrete[c])
end

num_kdes(m::HybridNB) = num_features(m)[1]
num_discrete(m::HybridNB) = num_features(m)[2]

"""
    HybridNB(labels::Vector{Int64}) -> model_h

    HybridNB(labels::Vector{Int64}, AstractString) -> model_h

A constructor for both types of features
"""
function HybridNB(labels::Vector{C}, ::Type{T}=Symbol) where {C <: Integer, T}
    c_kdes = Dict{C, Dict{T, InterpKDE}}()
    c_discrete = Dict{C, Dict{T, ePDF}}()
    priors = Dict{C, Float64}()
    classes = unique(labels)
    for class in classes
        c_kdes[class] = Dict{T, InterpKDE}()
        c_discrete[class] = Dict{T, ePDF}()
    end
    HybridNB{C, T}(c_kdes, c_discrete, classes, priors)
end

# Initialize with the number of continuous and discrete features
function HybridNB(labels::AbstractVector, num_kde::Int = 0, num_discrete::Int = 0)
    return HybridNB(labels, 1:num_kde, 1:num_discrete)
end

function Base.show(io::IO, m::HybridNB)
    println(io, "HybridNB")
    println(io, "  Classes = $(keys(m.c_kdes))")
end