Merge a5ff91d into ea33dbe

.travis.yml

@@ -1,6 +1,5 @@
 language: julia
 julia:
-    - 0.4
     - 0.5
     - nightly
 matrix:

REQUIRE

@@ -1,3 +1,5 @@
-julia 0.4
+julia 0.5
 QuantEcon
 Compat 0.9
+Combinatorics
+Iterators

demo/basis_mat_formats.jl

@@ -24,7 +24,7 @@ c1 = kron(bmt.vals[2], bmt.vals[1]) \ f
 # calling funfitxy with the BasisMatrix instance will use optimized version of
 # that operation that never constructs full kronecker product
 c2 = funfitxy(basis, bmt, f)[1]
-@assert maxabs(c1 - c2) < 1e-14
+@assert maximum(abs, c1 - c2) < 1e-14
 
 ## Direct form
 bmd = BasisMatrix(basis, Direct(), S)
@@ -42,21 +42,21 @@ bmd = BasisMatrix(basis, Direct(), S)
 # Now to get coefficients we will need to compute row_kron(bm.vals[2], bm.vals[2]) \ f
 c3 = row_kron(bmd.vals[2], bmd.vals[1]) \ f
 
-@assert maxabs(c3 - c2) < 1e-14
+@assert maximum(abs, c3 - c2) < 1e-14
 
 # calling funfitxy with the BasisMatrix instance will use optimized version of
 # that operation that never constructs row wise kronecker product
 c4 = funfitxy(basis, bmd, f)[1]
 
-@assert maxabs(c4 - c2) < 1e-14
+@assert maximum(abs, c4 - c2) < 1e-14
 
 ## Expanded form
 bme = BasisMatrix(basis, Expanded(), S)
 
 # to get coefficients we just need to do bme.vals[1] \ f
 c5 = bme.vals[1] \ f
 
-@assert maxabs(c5 - c2) < 1e-14
+@assert maximum(abs, c5 - c2) < 1e-14
 
 # the reason this worked is that Expanded form computes the fully expanded
 # version of the basis matrix. This means that the expansion down columns and
@@ -65,7 +65,7 @@ c5 = bme.vals[1] \ f
 # the funfitxy method for the Expanded BasisMatrix does exactly this operation
 c6 = funfitxy(basis, bme, f)[1]
 
-@assert maxabs(c6 - c2) < 1e-14
+@assert maximum(abs, c6 - c2) < 1e-14
 
 
 #=
@@ -107,33 +107,33 @@ rows of the matrix argument.
     hcat([repeat(bmt.vals[2][:, i], inner=length(x)) for i in 1:length(y)]...)
     )
 
-@assert maxabs(bmd.vals[1] - Φd1[1]) < 1e-14
-@assert maxabs(bmd.vals[2] - Φd1[2]) < 1e-14
+@assert maximum(abs, bmd.vals[1] - Φd1[1]) < 1e-14
+@assert maximum(abs, bmd.vals[2] - Φd1[2]) < 1e-14
 
 ## Tensor -> expanded
 # we've already seen this one as just kron(bmt.vals[2], bmt.vals[1]) (think
 # about how we computed the c1 and c5)
 Φe1 = kron(bmt.vals[2], bmt.vals[1])
-@assert maxabs(Φe1 - bme.vals[1]) < 1e-14
+@assert maximum(abs, Φe1 - bme.vals[1]) < 1e-14
 
 ## Direct -> Expanded
 # We've also seen that this is row_kron(bmd.vals[2], bmd.vals[1]) (check c3 and
 # c5)
 Φe2 = row_kron(bmd.vals[2], bmd.vals[1])
-@assert maxabs(Φe2 - bme.vals[1]) < 1e-14
+@assert maximum(abs, Φe2 - bme.vals[1]) < 1e-14
 
 # The julia methods `convert(Format, bm, [order])` do these operations for you,
 # but return a fully formed BasisMatrix object instead of just a plain Julia
 # AbstractMatrix subtype
 Φd2 = convert(Direct, bmt)
-@assert maxabs(bmd.vals[1] - Φd2.vals[1]) < 1e-14
-@assert maxabs(bmd.vals[2] - Φd2.vals[2]) < 1e-14
+@assert maximum(abs, bmd.vals[1] - Φd2.vals[1]) < 1e-14
+@assert maximum(abs, bmd.vals[2] - Φd2.vals[2]) < 1e-14
 
 Φe3 = convert(Expanded, bmt)
-@assert maxabs(Φe3.vals[1] - bme.vals[1]) < 1e-14
+@assert maximum(abs, Φe3.vals[1] - bme.vals[1]) < 1e-14
 
 Φe4 = convert(Expanded, bmd)
-@assert maxabs(Φe4.vals[1] - bme.vals[1]) < 1e-14
+@assert maximum(abs, Φe4.vals[1] - bme.vals[1]) < 1e-14
 
 ### Evaluation off nodes
 
@@ -158,19 +158,19 @@ ft1 = kron(bmt.vals[2], Φ1_x2)*c1
 
 # the approximation actually doesn't fit func extremely well, so we will have a
 # relatively loose sense of success in our interpolation
-@assert maxabs(ft1 - f2) < 1.0
+@assert maximum(abs, ft1 - f2) < 1.0
 
 # if we didn't want to form the full `kron` we could have construct a
 # BasisMatrix instance by hand (NOTE: to do this we had to do a relatively ugly
 # thing and allocate vals first, then populate it and we had to pass type
 # params to the BasisMatrix constructor. We can get around this, we just need
 # to be more clever in how we accept arguments to BasisMatrix).
-bmt2_vals = Array(typeof(Φ1_x2), 1, 2)
+bmt2_vals = Array{typeof(Φ1_x2)}(1, 2)
 bmt2_vals[1] = Φ1_x2; bmt2_vals[2] = bmt.vals[2]
 bmt2 = BasisMatrix{Tensor,typeof(Φ1_x2)}([0 0], bmt2_vals)
 ft2 = funeval(c1, bmt2)
 
-@assert maxabs(ft1 - ft2) < 1e-14
+@assert maximum(abs, ft1 - ft2) < 1e-14
 
 ## Using Direct form
 # to use direct form, we first need to construct the expanded grid on x2 and y
@@ -179,12 +179,12 @@ bmd2 = BasisMatrix(basis, Direct(), S2)
 
 # now we can evaluate using row_kron(bmd2.vals[2], bmd2.vals[1]) * c
 fd1 = row_kron(bmd2.vals[2], bmd2.vals[1]) * c1
-@assert maxabs(fd1 - ft1) < 1e-14
+@assert maximum(abs, fd1 - ft1) < 1e-14
 
 # we could also use the funeval method which would help us avoid building that
 # full row_kron matrix
 fd2 = funeval(c1, bmd2)
-@assert maxabs(fd2 - ft1) < 1e-14
+@assert maximum(abs, fd2 - ft1) < 1e-14
 
 
 ## Using Expanded form
@@ -194,8 +194,8 @@ bme2 = BasisMatrix(basis, Expanded(), S2)
 # evaluation now is juse bme2.vals[2] * c
 fe1 = bme2.vals[1]*c1
 
-@assert maxabs(fe1 - ft1) < 1e-14
+@assert maximum(abs, fe1 - ft1) < 1e-14
 
 # Here the funeval version is identical to the operation from above
 fe2 = funeval(c1, bme2)
-@assert maxabs(fe2 - ft1) < 1e-14
+@assert maximum(abs, fe2 - ft1) < 1e-14

src/BasisMatrices.jl

@@ -2,43 +2,7 @@ __precompile__()
 
 module BasisMatrices
 
-#=
-TODO: still need to write fund, minterp
-TODO: also need splidop, lindop
-TODO: funeval fails for scalar input and does weird thing for 1-element
-      vector input
-
-TODO: rethink the `order::Matrix` thing. Maybe a convention we can adopt is
-the following:
-
-- If calling evalbase on a `*Params` type, we will say that order must be of
-type `Integer` or `T<:Integer, Vector{T}`. The implementation for the `Vector`
-case will just loop over each element and then call the `Integer` method,
-which will have all the code
-- If calling `evalbase` on a `Basis{1}` type, the semantics will be the same
-as for the `*Params` case
-- If calling `evalbase` on a `Basis{N>1}` type, then `order` can either be a
-`T<:Integer Vector{T}` or `T<:Integer Matrix{T}`. If it is a vector, we call
-the `Int` evalbase for the `*Params` along each dimension, then combine as
-necessary. If it is a `Matrix` then we say that each column is the orders that
-should be computed for each dimension and we farm it out that way.
-
-So method signatures would be
-
-```
-evalbase(p::ChebParam, order::Int) => Matrix{Float64}
-evalbase(p::SplineParam, order::Int) => SparseMatrixCSC{Float64}
-evalbase(p::LinParam, order::Int) => SparseMatrixCSC{Float64}
-
-evalbase(p::ChebParam, order::Vector{Int}) => Vector{Matrix{Float64}}
-evalbase(p::SplineParam, order::Vector{Int}) => Vector{SparseMatrixCSC{Float64}}
-evalbase(p::LinParam, order::Vector{Int}) => Vector{SparseMatrixCSC{Float64}}
-
-evalbase(b::Basis{1}, order::Vector{Int}) => .... can't remember
-evalbase(b::Basis{N}, order::Matrix{Int}) => Vector{... can't remembers}
-```
-=#
-
+# TODO: still need to write fund, minterp
 
 #=
 Note that each subtype of `BT<:BasisFamily` (with associated `PT<:BasisParam`)
@@ -56,14 +20,17 @@ nodes(::PT)
 =#
 
 import Base: ==, *, \
+using Base.Cartesian
 
 using QuantEcon: gridmake, gridmake!, ckron, fix, fix!
 
-import Compat
+using Compat
+using Combinatorics: with_replacement_combinations
+using Iterators: product
 
 # types
-export BasisFamily, Cheb, Lin, Spline, Basis,
-       BasisParams, ChebParams, LinParams, SplineParams,
+export BasisFamily, Cheb, Lin, Spline, Basis, Smolyak,
+       BasisParams, ChebParams, LinParams, SplineParams, SmolyakParams,
        AbstractBasisMatrixRep, Tensor, Expanded, Direct,
        BasisMatrix, Interpoland, SplineSparse, RowKron
 
@@ -72,19 +39,62 @@ export nodes, get_coefs, funfitxy, funfitf, funeval, evalbase,
        derivative_op, row_kron, evaluate, fit!, update_coefs!,
        complete_polynomial, complete_polynomial!, n_complete
 
+#re-exports
+export gridmake, gridmake!, ckron
+
+abstract BasisFamily
+abstract BasisParams
+typealias TensorX Union{Tuple{Vararg{AbstractVector}},AbstractVector{TypeVar(:TV,AbstractVector)}}
+typealias IntSorV Union{Int, AbstractVector{Int}}
+
 include("util.jl")
 include("spline_sparse.jl")
-include("basis.jl")
-include("basis_structure.jl")
-include("interp.jl")
 
-# include the rest of the Julian API
+# include the families
+
+# BasisParams interface
+Base.issparse{T<:BasisParams}(::Type{T}) = false
+Base.ndims(::BasisParams) = 1
+for f in [:family, :family_name, :(Base.issparse), :(Base.eltype)]
+    @eval $(f){T<:BasisParams}(::T) = $(f)(T)
+end
 include("cheb.jl")
-include("spline.jl")
 include("lin.jl")
-
-# include comlpete
+include("spline.jl")
 include("complete.jl")
+include("smolyak.jl")
+
+# now some more interface methods that only make sense once we have defined
+# the subtypes
+
+# give the type of the `vals` field based on the family type parameter of the
+# corresponding basis. `Spline` and `Lin` use sparse, `Cheb` uses dense
+# a hybrid must fall back to a generic AbstractMatrix{Float64}
+# the default is just a dense matrix
+bmat_type{TP<:BasisParams}(::Type{TP}) = Matrix{eltype(TP)}
+bmat_type{TP<:BasisParams,T2}(::Type{T2}, ::Type{TP}) = Matrix{eltype(TP)}
+
+# default to SparseMatrixCSC
+bmat_type{T}(::Union{Type{LinParams{T}},Type{SplineParams{T}}}) = SparseMatrixCSC{eltype(T),Int}
+bmat_type{T,T2}(::Type{T2}, ::Union{Type{LinParams{T}},Type{SplineParams{T}}}) = SparseMatrixCSC{eltype(T),Int}
+
+# specialize for SplineSparse
+bmat_type{T,T2<:SplineSparse}(::Type{T2}, ::Union{Type{LinParams{T}},Type{SplineParams{T}}}) =
+    SplineSparse{eltype(T),Int}
+
+# add methods to instances
+bmat_type{T<:BasisParams}(::T) = bmat_type(T)
+bmat_type{TF<:BasisParams,T2}(ss::Type{T2}, ::TF) = bmat_type(T2, TF)
+
+# default method for evalbase with extra type hint is to just ignore the extra
+# type hint
+evalbase{T}(::Type{T}, bp::BasisParams, x, order) = evalbase(bp, x, order)
+
+
+# include other
+include("basis.jl")
+include("basis_structure.jl")
+include("interp.jl")
 
 
 # deprecations