minor changes to major change

Author: jverzani

README.md

@@ -22,6 +22,7 @@ julia> using Polynomials
 
 * `Polynomial` - Standard polynomials
 * `ChebyshevT` - Chebyshev polynomials of the first kind
+* `Bersntein` - Bernstein polynomials of degree n
 
 #### Construction and Evaluation
 
@@ -130,8 +131,8 @@ julia> roots(Polynomial([1, 0, -1]))
 
 julia> roots(Polynomial([1, 0, 1]))
 2-element Array{Complex{Float64},1}:
- 0.0+1.0im
- 0.0-1.0im
+ 0.0 - 1.0im
+ 0.0 + 1.0im
 
 julia> roots(Polynomial([0, 0, 1]))
 2-element Array{Float64,1}:
@@ -146,11 +147,11 @@ Fit a polynomial (of degree `deg`) to `x` and `y` using a least-squares approxim
 ```julia
 julia> xs = 0:4; ys = @. exp(-xs) + sin(xs);
 
-julia> fit(xs, ys)
-Polynomial(1.0000000000000016 + 0.059334723072240664*x + 0.39589720602859824*x^2 - 0.2845598112184312*x^3 + 0.03867830809692903*x^4)
+julia> fit(xs, ys) |> x -> round(x, digits=4)
+Polynomial(1.0 + 0.0593*x + 0.3959*x^2 - 0.2846*x^3 + 0.0387*x^4)
 
-julia> fit(ChebyshevT, xs, ys, deg=2)
-ChebyshevT([0.541280671210034, -0.8990834124779993, -0.4237852336242923])
+julia> fit(ChebyshevT, xs, ys, deg=2) |> x -> round(x, digits=4)
+ChebyshevT(0.5413⋅T_0(x) - 0.8991⋅T_1(x) - 0.4238⋅T_2(x))
 ```
 
 Visual example:
@@ -166,15 +167,16 @@ Polynomial objects also have other methods:
 
 * `coeffs`: returns the entire coefficient vector
 
-* `degree`: returns the polynomial degree, `length` is 1 plus the degree
+* `degree`: returns the polynomial degree, `length` is number of stored coefficients
 
-* `variable`: returns the polynomial symbol as a degree 1 polynomial
+* `variable`: returns the polynomial symbol as polynomial in the underlying type
 
 * `norm`: find the `p`-norm of a polynomial
 
-* `conj`: finds the conjugate of a polynomial over a complex fiel
+* `conj`: finds the conjugate of a polynomial over a complex field
 
 * `truncate`: set to 0 all small terms in a polynomial;
+
 * `chop` chops off any small leading values that may arise due to floating point operations.
 
 * `gcd`: greatest common divisor of two polynomials.

docs/src/extending.md

@@ -6,7 +6,7 @@ The [`AbstractPolynomial`](@ref) type was made to be extended via a rich interfa
 AbstractPolynomial
 ```
 
-A polynomial's  coefficients  are  relative to some *basis*. The `Polynomial` type relates coefficients  `[a0, a1,  ..., an]`, say,  to the  polynomial  `a0 +  a1*x + a2*x^  + ... +  an*x^n`,  through the natural  basis  `1,  x,  x^2, ..., x^n`.  New p olynomial  types typically represent the polynomial through a different  basis. For example,  `CheyshevT` uses a basis  `T_0=1, T_1=x,  T_2=2x^2-1,  ...,  T_n  =  2xT_{n-1} - T_{n-2}`.  For this type  the  coefficients  `[a0,a1,...,an]` are associated with  the polynomial  `a0*T0  + a1*T_1 +  ...  +  an*T_n`.
+A polynomial's  coefficients  are  relative to some *basis*. The `Polynomial` type relates coefficients  `[a0, a1,  ..., an]`, say,  to the  polynomial  `a0 +  a1*x + a2*x^  + ... +  an*x^n`,  through the standard  basis  `1,  x,  x^2, ..., x^n`.  New polynomial  types typically represent the polynomial through a different  basis. For example,  `CheyshevT` uses a basis  `T_0=1, T_1=x,  T_2=2x^2-1,  ...,  T_n  =  2xT_{n-1} - T_{n-2}`.  For this type  the  coefficients  `[a0,a1,...,an]` are associated with  the polynomial  `a0*T0  + a1*T_1 +  ...  +  an*T_n`.
 
 To implement a new polynomial type, `P`, the following methods should
 be implemented.

docs/src/index.md

@@ -104,8 +104,8 @@ ERROR: Polynomials must have same variable
 ### Integrals and Derivatives
 
 Integrate the polynomial `p` term by term, optionally adding constant
-term `C`. The order of the resulting polynomial is one higher than the
-order of `p`.
+term `C`. The degree of the resulting polynomial is one higher than the
+degree of `p`.
 
 ```jldoctest
 julia> integrate(Polynomial([1, 0, -1]))
@@ -115,8 +115,8 @@ julia> integrate(Polynomial([1, 0, -1]), 2)
 Polynomial(2.0 + 1.0*x - 0.3333333333333333*x^3)
 ```
 
-Differentiate the polynomial `p` term by term. The order of the
-resulting polynomial is one lower than the order of `p`.
+Differentiate the polynomial `p` term by term. The degree of the
+resulting polynomial is one lower than the degree of `p`.
 
 ```jldoctest
 julia> derivative(Polynomial([1, 3, -1]))
@@ -137,8 +137,8 @@ julia> roots(Polynomial([1, 0, -1]))
 
 julia> roots(Polynomial([1, 0, 1]))
 2-element Array{Complex{Float64},1}:
- -0.0 + 1.0im
   0.0 - 1.0im
+  0.0 + 1.0im
 
 julia> roots(Polynomial([0, 0, 1]))
 2-element Array{Float64,1}:
@@ -148,13 +148,13 @@ julia> roots(Polynomial([0, 0, 1]))
 
 ### Fitting arbitrary data
 
-Fit a polynomial (of order `deg`) to `x` and `y` using a least-squares approximation.
+Fit a polynomial (of degree `deg`) to `x` and `y` using polynomial interpolation or a (weighted) least-squares approximation.
 
 ```@example
 using Plots, Polynomials
 xs = range(0, 10, length=10)
 ys = exp.(-xs)
-f = fit(xs, ys)
+f = fit(xs, ys)  # fit(xs, ys, k)  for fitting a kth degreee polynomial
 
 scatter(xs, ys, label="Data");
 plot!(f, extrema(xs)..., label="Fit");
@@ -163,9 +163,47 @@ savefig("polyfit.svg"); nothing # hide
 
 ![](polyfit.svg)
 
+### Other bases
+
+A polynomial, e.g. `a_0 + a_1 x + a_2 x^2 + ... + a_n x^n`, can be seen as a collection of coefficients, `[a_0, a_1, ..., a_n]`, relative to some polynomial basis. The most  familiar basis being  the standard one: `1`, `x`, `x^2`, ...  Alternative bases are possible.  The `ChebyshevT` polynomials are  implemented, as an example. Instead of `Polynomial`  or  `Polynomial{T}`, `ChebyshevT` or  `ChebyshevT{T}` constructors are used:
+
+```jldoctest
+julia> p1 = ChebyshevT([1.0, 2.0, 3.0])
+ChebyshevT(1.0⋅T_0(x) + 2.0⋅T_1(x) + 3.0⋅T_2(x))
+
+julia> p2 = ChebyshevT{Float64}([0, 1, 2])
+ChebyshevT(1.0⋅T_1(x) + 2.0⋅T_2(x))
+
+julia> p1 + p2
+ChebyshevT(1.0⋅T_0(x) + 3.0⋅T_1(x) + 5.0⋅T_2(x))
+
+julia> p1 * p2
+ChebyshevT(4.0⋅T_0(x) + 4.5⋅T_1(x) + 3.0⋅T_2(x) + 3.5⋅T_3(x) + 3.0⋅T_4(x))
+
+julia> derivative(p1)
+ChebyshevT(2.0⋅T_0(x) + 12.0⋅T_1(x))
+
+julia> integrate(p2)
+ChebyshevT(0.25⋅T_0(x) - 1.0⋅T_1(x) + 0.25⋅T_2(x) + 0.3333333333333333⋅T_3(x))
+
+julia> convert(Polynomial, p1)
+Polynomial(-2.0 + 2.0*x + 6.0*x^2)
+
+julia> convert(ChebyshevT, Polynomial([1.0, 2,  3]))
+ChebyshevT(2.5⋅T_0(x) + 2.0⋅T_1(x) + 1.5⋅T_2(x))
+```
+
+
 ### Iteration
 
-A polynomial, e.g. `a_0 + a_1 x + a_2 x^2 + ... + a_n x^n` can be seen as a collection of coefficients, `[a_0, a_1, ..., a_n]`, reelative to some polynomial basis. The most  familiar basis being  `1`, `x`, `x^2`, ... If the basis is implicit, then a polynomial is just a vector. Vectors or 1-based, but polynomial types are 0-based, for purposes of indexing (e.g. `getindex`, `setindex!`, `eachindex`). Iteration over a polynomial steps through the basis vectors, e.g. `a_0`, `a_1 x`, ...
+If its basis is implicit, then a polynomial may be  seen as just a vector of  coefficients. Vectors or 1-based, but, for convenience, polynomial types are 0-based, for purposes of indexing (e.g. `getindex`, `setindex!`, `eachindex`). Iteration over a polynomial steps through the basis vectors, e.g. `a_0`, `a_1*x`, ...
+
+```jldoctest
+julia> as = [1,2,3,4,5]; p = Polynomial(as);
+
+julia> as[3], p[2], collect(p)[3]
+(3, 3, Polynomial(3*x^2))
+```
 
 ## Related Packages
 

docs/src/polynomials/bernstein.md

@@ -7,7 +7,7 @@ end
 ```
 
 
-The [Bernstein polynomials](https://en.wikipedia.org/wiki/Bernstein_polynomial) are a family of polynomials defined on the interval `[0,1]`. For each `n` there are `n+1` polynomials, given by: `b(n,nu) = choose(n,nu) x^nu * (1-x)^(n-nu)` for `nu` in `0:n`. Together, these form a basis that can represent any polynomial of degree `n` or less through a linear combination.
+The [Bernstein polynomials](https://en.wikipedia.org/wiki/Bernstein_polynomial) are a family of polynomials defined on the interval `[0,1]`. For each `n` there are `n+1` polynomials, given by: `b(n,nu) = choose(n,nu) x^nu * (1-x)^(n-nu)` for `nu` in `0:n`. Together, these form a basis that can represent any polynomial of degree `n` or less through a unique linear combination.
 
 
 ## Bernstein(N,T)

docs/src/polynomials/chebyshev.md

@@ -7,7 +7,7 @@ end
 ```
 
 
-The [Chebyshev polynomials](https://en.wikipedia.org/wiki/Chebyshev_polynomials) are two sequences of polynomials, `T_n` and `U_n`. The Chebyshev polynomials of the first kind, `T_n`, can be defined by the recurrence relation `T_0(x)=1`, `T_1(x)=x`, and `T_{n+1}(x) = 2xT_n{x}-T_{n-1}(x)`. The Chebyshev polynomioals of the second kind, `U_n(x)`, can be defined by `U_0(x)=1`, `U_1(x)=2x`, and `U_{n+1}(x) = 2xU_n(x) - U_{n-1}(x)`. Both `T_n` and `U_n` have degree `n`, and any polynomial  of  degree `n` may be written as a linear combination of the polynomials `T_0`, `T_1`, ..., `T_n` (similarly with `U`).
+The [Chebyshev polynomials](https://en.wikipedia.org/wiki/Chebyshev_polynomials) are two sequences of polynomials, `T_n` and `U_n`. The Chebyshev polynomials of the first kind, `T_n`, can be defined by the recurrence relation `T_0(x)=1`, `T_1(x)=x`, and `T_{n+1}(x) = 2xT_n{x}-T_{n-1}(x)`. The Chebyshev polynomioals of the second kind, `U_n(x)`, can be defined by `U_0(x)=1`, `U_1(x)=2x`, and `U_{n+1}(x) = 2xU_n(x) - U_{n-1}(x)`. Both `T_n` and `U_n` have degree `n`, and any polynomial  of  degree `n` may be uniquely written as a linear combination of the polynomials `T_0`, `T_1`, ..., `T_n` (similarly with `U`).
 
 
 ## First Kind

docs/src/reference.md

@@ -91,6 +91,7 @@ plot(::AbstractPolynomial, a, b; kwds...)
 will plot the polynomial within the range `[a, b]`.
 
 ### Example: The Polynomials.jl logo
+
 ```@example
 using Plots, Polynomials
 # T1, T2, T3, and T4:

src/Polynomials.jl

@@ -11,10 +11,10 @@ include("contrib.jl")
 # Polynomials
 include("polynomials/Polynomial.jl")
 include("polynomials/ChebyshevT.jl")
-include("polynomials/ChebyshevU.jl")
+#include("polynomials/ChebyshevU.jl")
 include("polynomials/Bernstein.jl")
 
-include("polynomials/Poly.jl") # Deprecated -> Will be removed
+include("polynomials/Poly.jl") # to be deprecated, then removed
 include("pade.jl")
 include("compat.jl") # Where we keep deprecations
 

src/common.jl

@@ -60,39 +60,40 @@ fromroots(A::AbstractMatrix{T}; var::SymbolLike = :x) where {T <: Number} =
     fromroots(Polynomial, eigvals(A), var = var)
 
 """
-    fit(x, y; [weights], deg=length(x) - 1, var=:x)
-    fit(::Type{<:AbstractPolynomial}, x, y; [weights], deg=length(x)-1, var=:x)
+    fit(x, y, deg=length(x) - 1; [weights], var=:x)
+    fit(::Type{<:AbstractPolynomial}, x, y, deg=length(x)-1; [weights], var=:x)
 
-Fit the given data as a polynomial type with the given degree. Uses linear least squares. When weights are given, as either a `Number`, `Vector` or `Matrix`, will use weighted linear least squares. The default polynomial type is [`Polynomial`](@ref). This will automatically scale your data to the [`domain`](@ref) of the polynomial type using [`mapdomain`](@ref)
+Fit the given data as a polynomial type with the given degree. Uses linear least squares. When weights are given, as either a `Number`, `Vector` or `Matrix`, will use weighted linear least squares. The default polynomial type is [`Polynomial`](@ref). This will automatically scale your data to the [`domain`](@ref) of the polynomial type using [`mapdomain`](@ref). To specify a different range to scale to, specify `domain=(a,b)` where `a <= minimum(xs) && maximum(xs) <= b`.
 """
 function fit(P::Type{<:AbstractPolynomial},
-    x::AbstractVector{T},
-    y::AbstractVector{T};
-    weights = nothing,
-    deg::Integer = length(x) - 1,
-    var = :x,) where {T}
-    x = mapdomain(P, x)
+             x::AbstractVector{T},
+             y::AbstractVector{T},
+             deg::Integer=length(x) - 1;
+             domain=(minimum(x), maximum(x)),
+             weights = nothing,
+             var = :x,) where {T}
+    x = mapdomain(P, first(domain),last(domain)).(x)
     vand = vander(P, x, deg)
     if weights !== nothing
         coeffs = _wlstsq(vand, y, weights)
     else
-        coeffs = pinv(vand) * y
+        coeffs = vand \ y #pinv(vand) * y
     end
     return P(T.(coeffs), var)
 end
 
 fit(P::Type{<:AbstractPolynomial},
     x,
-    y;
+    y,
+    deg::Integer = length(x) - 1;
     weights = nothing,
-    deg::Integer = length(x) - 1,
-    var = :x,) = fit(P, promote(collect(x), collect(y))...; weights = weights, deg = deg, var = var)
+    var = :x,) = fit(P, promote(collect(x), collect(y))..., deg; weights = weights, var = var)
 
 fit(x::AbstractVector,
-    y::AbstractVector;
+    y::AbstractVector,
+    deg::Integer = length(x) - 1;
     weights = nothing,
-    deg::Integer = length(x) - 1,
-    var = :x,) = fit(Polynomial, x, y; weights = weights, deg = deg, var = var)
+    var = :x,) = fit(Polynomial, x, y, deg; weights = weights, var = var)
 
 # Weighted linear least squares
 _wlstsq(vand, y, W::Number) = _wlstsq(vand, y, fill!(similar(y), W))
@@ -353,6 +354,11 @@ function mapdomain(P::Type{<:AbstractPolynomial}, x::AbstractArray)
 end
 mapdomain(::P, x::AbstractArray) where {P <: AbstractPolynomial} = mapdomain(P, x)
 
+function mapdomain(P::Type{<:AbstractPolynomial}, a::Number, b::Number)
+    a, b = a < b ? (a,b) : (b,a)
+    x -> mapdomain(P, [a,x,b])[2]
+end
+
 #=
 indexing =#
 Base.firstindex(p::AbstractPolynomial) = 0
@@ -362,10 +368,14 @@ Base.broadcastable(p::AbstractPolynomial) = Ref(p)
 
 # basis
 # return the kth basis polynomial for the given polynomial type, e.g. x^k for Polynomial{T}
-function basis(p::P, k::Int) where {P <: AbstractPolynomial}
-    zs = zeros(eltype(p), k+1)
-    zs[k+1] = 1
-    P(zs, p.var)
+function basis(p::P, k::Int) where {P<:AbstractPolynomial}
+    basis(P, k)
+end
+
+function basis(::Type{P}, k::Int; var=:x) where {P <: AbstractPolynomial}
+    zs = zeros(Int, k+1)
+    zs[end] = 1
+    P(zs, var)
 end
 
 # iteration

src/compat.jl

@@ -19,8 +19,8 @@ end
 polyint(p::AbstractPolynomial, C = 0) = integrate(p, C)
 polyint(p::AbstractPolynomial, a, b) = integrate(p, a, b)
 polyder(p::AbstractPolynomial, ord = 1) = derivative(p, ord)
-polyfit(x, y, n = length(x) - 1) = fit(Polynomial, x, y; deg = n)
-polyfit(x, y, sym::Symbol) = fit(Polynomial, x, y; var = sym)
+polyfit(x, y, n = length(x) - 1, sym=:x) = fit(Polynomial, x, y, n; var = sym)
+polyfit(x, y, sym::Symbol) = fit(Polynomial, x, y, var = sym)
 
 padeval(PQ::Pade, x::Number) = PQ(x)
 padeval(PQ::Pade, x) = PQ.(x)