quantecon / compecon.jl Goto Github PK

We need to make it clearer that in order to use a higher degree spline you need to produce y data for x = nodes(basis)[1].

Trying to use the initial nodes passed as the knot vector will not work because we pad the endpoints internally to preserve continuity at the borders.

Since Julia handles loops very well, we should see if a 1D golden method solver looped over the vector function is faster than using the Matlabish vector golden search

If you have a Basis{N} N>1 and you try to construct a BasisStructure at a single point, you cannot represent that point as a Vector{Float64} of length N.

Test case to reproduce:

julia> using CompEcon

julia> basis = Basis(ChebParams(10, -1, 1), ChebParams(5, 2, 4))
2 dimensional Basis on the hypercube formed by [-1.0,2.0] × [1.0,4.0].
Basis families are CompEcon.Cheb × CompEcon.Cheb


julia> pt = [0.5, 2.1]
2-element Array{Float64,1}:
 0.5
 2.1

julia> BasisStructure(basis, Direct(), pt)
ERROR: BoundsError: attempt to access 2-element Array{Float64,1}:
 0.5
 2.1
  at index [Colon(),2]
 in throw_boundserror at abstractarray.jl:155
 in call at /Users/sglyon/.julia/v0.4/CompEcon/src/basis_structure.jl:160
 in call at /Users/sglyon/.julia/v0.4/CompEcon/src/basis_structure.jl:142

julia> BasisStructure(basis, Direct(), pt')
BasisStructure{CompEcon.Direct} of order [0 0]

Not an issue per se, but I have two questions about the way you define types:

• Why do you use immutable types? It sounds like immutable types are passed by values, which does not sound efficient. From the documentation:

  " An object with an immutable type is passed around (both in assignment statements and in function calls) by copying, whereas a mutable type is passed around by reference."

• I think the documentation advises to avoid abstract types in field definitions.
  http://julia.readthedocs.org/en/latest/manual/faq/#how-do-abstract-or-ambiguous-fields-in-types-interact-with-the-compiler.

We might want to write a custom sparse data-type for use inside BasisStructure.vals when we have a finite element basis (spline, lin). Here a rough sketch of something that might work

immutable BasisAtPoint{Tv,Ti}
    f_vals::Vector{T}
    inds::UnitRange{Ti}
end

where length(f_vals) == length(inds) --> true and length(f_vals) <= degree_of_spline+1.

We could then collect all of these in a Vector{BasisAtPoint{Tv,Ti}}. When we eventually need to do the matrix operation to evaluate the interpoland we would simply need to do

# let bs.vals <: Vector{BasisAtPoint{Float64,Int}}, for example
function evaluate{Tv}(coefs::Vector{Float64}, bs::BasisStructure{Expanded})
    npoints = length(bs)
    out = Array(Tv, npoints)
    @inbounds for i=1:npoints
        out[i] = dot(bs.vals[i], bs.inds[i], coefs, bs.inds[i])
    end
    out
end

To make this work with minimal change to the current code base (not necessarily something we need to aim for), we would need to define getindex, row_kron, and * for this new type and make Vector{BasisAtPoint} <: AbstractMatrix

memory

This would mean that for each point in the domain where we want to construct a basis structure we would store only:

• the non-zero basis function evaluations at that point (at most k+1 numbers, usually Float64)
• The 2 integers that are start and stop for the UnitRange.

The best we could do using either CSR or CSC is to store the same non-zero basis functions, plus the same number of integers for the indval, plus at least one integer for the indptr.

We can gain memory efficiency here by leveraging the fact that we know the non-zero elements of each row are consecutive. We might be able to achieve the same memory lightness with some block sparse matrix type, but we don't even have the beginnings of that yet. This new data type might be a first step in that direction

Hi @spencerlyon2

how are you doing with this? I can see you are being very active in many places, Interpolations.jl for example, what's your feeling about this project here? I think it would be great to have a unifying framework a la CompEcon. I'm not saying this is necessarily the best way to achieve this, but it is a solid start. I like that one can construct different interpolations by just giving a different keyword.
What about pulling your Smolyak package into this here? This is nothing else than a different way of constructing a basis function, so should fit well with the setup here.

We should definitely be able to speed this up!

I don't think we have any tests

Right now the Sparse matrix version of row_kron simply constructs a sparse matrix and hand off to the dense version of the function. This duplicates a lot of work, we could save some time by only dealing directly with the non-zero elements

There is a lot of flexibility in how we allow users to construct/evaluate a Basis/BasisStructure ect.

This comes at a cost.

We really should try to specialize more on input types and clean things up...

A quick look at @code_warntype for most of our methods reveals the current status.

Provide docs/usage examples for matlab users

We have this already in the various convert methods for BasisStructure, but it would be great to expose another method to combine the output of funbase.

The theme of the new methods would be to make it easy for the user to combine correctly and hard for them to do it incorrectly.

I tested the working example (which I named by ChebshevApp.jl) and have the following error messages. My Julia version is 1.5.0

`julia> include("ChebshevApp.jl")
[ Info: Precompiling CompEcon [1375f997-1772-542d-b6fa-a8ee39124556]
WARNING: Method definition dot(Any, Any, Any) in module LinearAlgebra at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\generic.jl:924 overwritten in module Optim at C:\Users\user.julia\packages\Optim\Agd3B\src\multivariate\precon.jl:23.
** incremental compilation may be fatally broken for this module **

WARNING: Method definition dot(Any, Any, Any) in module LinearAlgebra at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\generic.jl:924 overwritten in module Optim at C:\Users\user.julia\packages\Optim\Agd3B\src\multivariate\precon.jl:23.
** incremental compilation may be fatally broken for this module **

WARNING: Method definition dot(Any, Any, Any) in module LinearAlgebra at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\generic.jl:924 overwritten in module Optim at C:\Users\user.julia\packages\Optim\Agd3B\src\multivariate\precon.jl:23.
** incremental compilation may be fatally broken for this module **

WARNING: Method definition dot(Any, Any, Any) in module LinearAlgebra at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\generic.jl:924 overwritten in module Optim at C:\Users\user.julia\packages\Optim\Agd3B\src\multivariate\precon.jl:23.
** incremental compilation may be fatally broken for this module **

ERROR: LoadError: MethodError: no method matching exp(::Array{Float64,1})
Closest candidates are:
exp(::Float16) at math.jl:1144
exp(::Missing) at math.jl:1197
exp(::BigFloat) at mpfr.jl:603
...
Stacktrace:
[1] f(::Array{Float64,1}) at C:\Users\user\Documents\Julia\ChebshevApp.jl:3
[2] funfitf(::Basis{1,Tuple{ChebParams{Float64}}}, ::Function) at C:\Users\user.julia\packages\BasisMatrices\PZ1uM\src\interp.jl:68
[3] funfitf(::Dict{Symbol,Any}, ::Function) at C:\Users\user.julia\packages\CompEcon\9POAs\src\core.jl:134
[4] top-level scope at C:\Users\user\Documents\Julia\ChebshevApp.jl:16
[5] include(::String) at .\client.jl:457
[6] top-level scope at none:1
in expression starting at C:\Users\user\Documents\Julia\ChebshevApp.jl:16`

Test case:

using CompEcon
basis = Basis(ChebParams(10, -1, 1), ChebParams(5, 2, 4));
grid = nodes(basis)[1];
y = grid[:, 1] - grid[:, 2];  # simple function x - y
c = funfitxy(basis, grid, y)[1];
pt = [0.5, 2.1];
bs = BasisStructure(basis, Direct(), pt);
funeval(c, bs)

REPL example:

julia> using CompEcon

julia> basis = Basis(ChebParams(10, -1, 1), ChebParams(5, 2, 4));


julia> grid = nodes(basis)[1];

julia> y = grid[:, 1] - grid[:, 2];  # simple function x - y

julia> c = funfitxy(basis, grid, y)[1];

julia> pt = [0.5, 2.1];

julia> bs = BasisStructure(basis, Direct(), pt);

julia> funeval(c, bs)
0-dimensional Array{Float64,0}:
-1.6

Hi,

Thanks for this great package!

Since this package use ReExport in it, it holds back the version of JLD2 to v0.4.3. This old version of JLD2 conflicts with the latest version of Julia (actually, this old version of JLD2 does not work well in any Julia 1.8x). Could you please fix this issue? Many thanks!

It turns out that if x1 and x2 are 1-dimensional, but the f function produces 2-dimensional vectors, then we get this error:

ERROR: BoundsError()
in getindex_bool_1d at array.jl:285
in getindex at array.jl:299
in golden_method at /home/ap/.julia/v0.3/CompEcon/src/optimization.jl:19

Could you explain a little bit more the order argument in BasisStructure and in funeval? In the documentation, you write "The order field keeps track of what order of derivative or integral the arrays inside vals correspond to". Why is it a matrix?
What's the best way to evaluate the zero order, first order, second order derivative of a set of Chebyshev coefficients on a given grid? I'm doing the following.

basis = Basis(Cheb, 50, 0.0, 1.0)
bs0 = BasisStructure(basis, Direct(), nodes(basis)[1], 0) 
bs1 = BasisStructure(basis, Direct(), nodes(basis)[1], 1) 
bs2 = BasisStructure(basis, Direct(), nodes(basis)[1], 2) 
funeval(coef, bs0, CompEcon._check_order(n, 0))
funeval(coef, bs1, CompEcon._check_order(n, 1))
funeval(coef, bs2, CompEcon._check_order(n, 2))

Is there a chance that this could be registered in Metadata? I have a bunch of stuff that relies on this, and it would be super handy if I could simply add this to a REQUIRE file and wouldn't have to add manual install instructions to a README.

The current version of golden_method doesn't work with functions that use scalars.

f(x::Real) = x^2

julia> golden_method(f, -1., 1.)
ERROR: MethodError: `*` has no method matching *(::Array{Float64,1}, ::Array{Float64,1})
Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...)
  *{T<:Union{Complex{Float32},Float64,Float32,Complex{Float64}},S}(::Union{DenseArray{T<:Union{Complex{Float32},Float64,Float32,Complex{Float64}},2},SubArray{T<:Union{Complex{Float32},Float64,Float32,Complex{Float64}},2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD}}, ::Union{SubArray{S,1,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD},DenseArray{S,1}})
  *{TA,TB}(::Base.LinAlg.AbstractTriangular{TA,S<:AbstractArray{T,2}}, ::Union{SubArray{TB,1,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD},DenseArray{TB,2},SubArray{TB,2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD},DenseArray{TB,1}})
  ...
 in power_by_squaring at intfuncs.jl:78
 in f at none:1
 in golden_method at /home/chase/.julia/v0.4/CompEcon.jl/src/optimization.jl:10
 in golden_method at /home/chase/.julia/v0.4/CompEcon.jl/src/optimization.jl:37

There are a number of things that could be sped up here. Most of them revolve around allocating less

• We could cache augbreaks based on minorder
• We can avoid calling vec so many times in the call to sparse
• We could try to not even use the sparse method and directly call the SparseMatrixCSC constructor
• Change things like collect(1:m)'' to reshape(1:m, 1, m)
• Don't make B a cell we know what it will really be (Vector{SparseMatrixCSC{Float64,Int64}})
• Roll out loops instead of any [:, ind] = stuff[:, ind]
• Use things like @inbounds or @simd where appropriate

I'm sure there are more, but that is what I came up with from a first look at the method through performance lenses

It looks like tests pass on v0.7, but there are some deprecation warnings, I print them below. Is there someone working on this at the moment? If not, I could have a try, but some guidance would be welcome.

┌ Warning: Deprecated syntax `parametric method syntax get_bformat{T <: BasisMatrix{Direct}}(b::T)` around C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:12.
│ Use `get_bformat(b::T) where T <: BasisMatrix{Direct}` instead.
└ @ C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:12
┌ Warning: Deprecated syntax `parametric method syntax get_bformat{T <: BasisMatrix{Expanded}}(b::T)` around C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:13.
│ Use `get_bformat(b::T) where T <: BasisMatrix{Expanded}` instead.
└ @ C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:13
┌ Warning: Deprecated syntax `parametric method syntax get_bformat{T <: BasisMatrix{Tensor}}(b::T)` around C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:14.
│ Use `get_bformat(b::T) where T <: BasisMatrix{Tensor}` instead.
└ @ C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:14
WARNING: Base.Test is deprecated, run `using Test` instead
  likely near C:\Users\Ilja\.julia\dev\CompEcon\test\runtests.jl:5
┌ Warning: `linspace(start, stop, length::Integer)` is deprecated, use `range(start, stop=stop,
length=length)` instead.
│   caller = top-level scope at none:0
└ @ Core none:0
┌ Warning: `linspace(start, stop, length::Integer)` is deprecated, use `range(start, stop=stop,
length=length)` instead.
│   caller = top-level scope at none:0
└ @ Core none:0
Test Summary:            | Pass  Total
Test original API compat |   10     10
┌ Warning: `linspace(start, stop, length::Integer)` is deprecated, use `range(start, stop=stop,
length=length)` instead.
│   caller = top-level scope at none:0
└ @ Core none:0
┌ Warning: `Array{T}(m::Int) where T` is deprecated, use `Array{T}(undef, m)` instead.
│   caller = fundef(::Tuple{Symbol,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Int64,Int64}) at core.jl:76
└ @ CompEcon C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:76
┌ Warning: `Array{T}(m::Int) where T` is deprecated, use `Array{T}(undef, m)` instead.
│   caller = fundef(::Tuple{Symbol,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Int64,Int64}) at core.jl:77
└ @ CompEcon C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:77
┌ Warning: `Array{T}(m::Int) where T` is deprecated, use `Array{T}(undef, m)` instead.
│   caller = fundef(::Tuple{Symbol,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Int64,Int64}) at core.jl:79
└ @ CompEcon C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:79
┌ Warning: `linspace(start, stop, length::Integer)` is deprecated, use `range(start, stop=stop,
length=length)` instead.
│   caller = top-level scope at none:0
└ @ Core none:0
┌ Warning: `Array{T}(m::Int) where T` is deprecated, use `Array{T}(undef, m)` instead.
│   caller = fundef(::Tuple{Symbol,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Int64}) at core.jl:76
└ @ CompEcon C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:76
┌ Warning: `Array{T}(m::Int) where T` is deprecated, use `Array{T}(undef, m)` instead.
│   caller = fundef(::Tuple{Symbol,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Int64}) at core.jl:77
└ @ CompEcon C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:77
┌ Warning: `Array{T}(m::Int) where T` is deprecated, use `Array{T}(undef, m)` instead.
│   caller = fundef(::Tuple{Symbol,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Int64}) at core.jl:79
└ @ CompEcon C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:79
┌ Warning: `Array{T}(m::Int) where T` is deprecated, use `Array{T}(undef, m)` instead.
│   caller = fundef(::Tuple{Symbol,Int64,Float64,Float64}) at core.jl:76
└ @ CompEcon C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:76
┌ Warning: `Array{T}(m::Int) where T` is deprecated, use `Array{T}(undef, m)` instead.
│   caller = fundef(::Tuple{Symbol,Int64,Float64,Float64}) at core.jl:77
└ @ CompEcon C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:77
┌ Warning: `Array{T}(m::Int) where T` is deprecated, use `Array{T}(undef, m)` instead.
│   caller = fundef(::Tuple{Symbol,Int64,Float64,Float64}) at core.jl:79
└ @ CompEcon C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:79
┌ Warning: `Array{T}(m::Int, n::Int) where T` is deprecated, use `Array{T}(undef, m, n)` instead.
│   caller = fundefn(::Symbol, ::Array{Int64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64) at core.jl:107
└ @ CompEcon C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:107
┌ Warning: `Array{T}(m::Int) where T` is deprecated, use `Array{T}(undef, m)` instead.
│   caller = fundef(::Array{Any,1}, ::Vararg{Array{Any,1},N} where N) at core.jl:76
└ @ CompEcon C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:76
┌ Warning: `Array{T}(m::Int) where T` is deprecated, use `Array{T}(undef, m)` instead.
│   caller = fundef(::Array{Any,1}, ::Vararg{Array{Any,1},N} where N) at core.jl:77
└ @ CompEcon C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:77
┌ Warning: `Array{T}(m::Int) where T` is deprecated, use `Array{T}(undef, m)` instead.
│   caller = fundef(::Array{Any,1}, ::Vararg{Array{Any,1},N} where N) at core.jl:79
└ @ CompEcon C:\Users\Ilja\.julia\dev\CompEcon\src\core.jl:79
┌ Warning: Deprecated syntax `implicit assignment to global variable `i``.
│ Use `global i` instead.
└ @ none:0
´´´

Before we add functionality or go for performance optimizations, we need tests!

Some of this has already been started, but we are far from being done here.

List of files we need tests for:

• src/basis.jl:
• src/basis_structure.jl:
• src/cheb.jl:
• src/interp.jl:
• src/lin.jl:
• src/optimization.jl:
• src/quad.jl: Done
• src/spline.jl:
• src/util.jl: For row_kron tests see here
• src/original/*.jl: These should be very simple and just test that if users call code with the original Matlab API, the code runs. Correctness checks should all happen in the files inside src/

This is to start a discussion about whether or not we should transfer this repository in to the QuantEcon organization.

I've been contacted recently by more than one source asking about the status of this code and the relationship to quantecon. Parties have suggested that with the repo in my name, with me being the primary contributor, there is a fair amount of uncertainty regarding the future maintenance of the code. If the repo were housed under the QuantEcon umbrella, that might send a stronger signal that we have at least medium run goals to keep this package going.

I use this code in my research, so I plan on maintaining it going forward. But, unless someone asks there's no way to send that signal.

I'd be happy to move this over. Before I do that I'd like to write up some quick docs. And make sure the rest of the team is on board.

cc @jstac, @cc7768, @mmcky, @albep, @vgregory757, @albop