QuantEcon / BasisMatrices.jl / 149

Committed 12 Feb 2018 - 18:36 coverage: 91.01%. First build

Build # 149

Build Type

Pull #48

travis-ci

Committed by

web-flow

Commit Message

TEST: fix testing error

Pull Request Pull Request #48: Avoid bounds violation with 0-degree Cheb. polynomial.

Run Details

2 of 6 new or added lines in 2 files covered. (33.33%)

1225 of 1346 relevant lines covered (91.01%)

486283.95 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

82.08

/src/cheb.jl

# --------------- #
# Chebyshev Basis #
# --------------- #

struct Cheb <: BasisFamily end

mutable struct ChebParams{T<:Number} <: BasisParams
    n::Int
    a::T
    b::T

    function ChebParams{T}(n::Int, a::T, b::T) where T
        n <= 0 && error("n must be positive")
        a >= b && error("left endpoint (a) must be less than right end point (b)")
        new{T}(n, a, b)
    end
end

ChebParams(n::Int, a::T, b::T) where {T<:Number} = ChebParams{T}(n, a, b)
ChebParams(n::Int, a::T, b::T) where {T<:Integer} = ChebParams(n, Float64(a), Float64(b))

## BasisParams interface
# define these methods on the type, the instance version is defined over
# BasisParams
family(::Type{T}) where {T<:ChebParams} = Cheb
family_name(::Type{T}) where {T<:ChebParams} = "Cheb"
@generated Base.eltype(::Type{T}) where {T<:ChebParams} = T.parameters[1]

# methods that only make sense for instances
Base.min(cp::ChebParams) = cp.a
Base.max(cp::ChebParams) = cp.b
Base.length(cp::ChebParams) = cp.n

function Base.show(io::IO, p::ChebParams)
    m = string("Chebyshev interpoland parameters with ",
               "$(p.n) basis functions from $(p.a), $(p.b)")
    print(io, m)
end

function nodes(p::ChebParams{T}, ::Type{Val{0}}) where T
    s = (p.b-p.a) / 2  # 21
    m = (p.b+p.a) / 2  # 22
    half = convert(T, 1/2)
    k = π*(half:(p.n - half))  # 25
    m .- cos.(k ./ p.n) .* s  # 26
end

function nodes(p::ChebParams, ::Type{Val{1}})
    x = nodes(p, Val{0})
    aa = x[1]
    bb = x[end]
    c1 = (bb*p.a - aa*p.b)/(bb-aa)
    c2 = (p.b.-p.a)/(bb-aa)

    @inbounds @simd for ix in eachindex(x)
        x[ix] = c1 + c2 * x[ix]
    end
    x
end

function nodes(p::ChebParams{T}, ::Type{Val{2}}) where T
    s = (p.b-p.a) / 2  # 21
    m = (p.b+p.a) / 2  # 22
    k = pi*(zero(T):(p.n - one(T)))  # 33
    m .- cos.(k ./ (p.n-1)) .* s  # 34
end

# chebnode.m -- DONE
nodes(p::ChebParams, nodetype::Integer=0) = nodes(p, Val{min(2, nodetype)})

function derivative_op(p::ChebParams, x, order=1)
    n, a, b = p.n, p.a, p.b
    if order > 0
        # TODO: figure out some caching mechanism that avoids globals
        D = Array{SparseMatrixCSC{basis_eltype(p, x),Int64}}(max(2, order)) # 49
        i = repmat(1:n', 1, n)
        j = i'  # 50

        # 51
        inds = find((rem.(i + j, 2) .== 1) .& (j .> i))
        r, c = similar(inds), similar(inds)
        for ix in 1:length(inds)
            r[ix], c[ix] = ind2sub((n, n), inds[ix])
        end

        d = sparse(r, c, (4/(b-a)) * (vec(j[1, c])-1), n-1, n)  # 52
        d[1, :] ./= 2  # 53
        D[1] = d  # 54
        for ii in 2:max(2, order)
            D[ii] = d[1:n-ii, 1:n-ii+1] * D[ii-1]  # 56
        end
    else
        D = Array{SparseMatrixCSC{basis_eltype(p, x),Int64}}(abs(order))  # 64
        nn = n - order  # 65
        z = (0.25 * (b - a)) ./(1:nn)  # 66
        d = sparse(vcat(1:nn, 1:nn-2), vcat(1:nn, 3:nn), vcat(z, -z[1:nn-2]),
                   nn, nn)  # 67
        d[1, 1] *= 2  # 68
        d0 = ((-1).^(0:nn-1)') .* sum(d, 1)  # 69
        D[1] = sparse(vcat(d0[1:n]', d[1:n, 1:n]))  # 70
        for ii=-2:-1:order
            ind = 1:n-ii-1
            D[-ii] = sparse(vcat(d0[ind]', d[ind, ind]) * D[-ii-1])
        end
    end
    D, ChebParams(n-order, a, b)
end

function evalbase(p::ChebParams, x::AbstractArray=nodes(p, 1), order::Int=0, nodetype::Int=1)
    n, a, b = p.n, p.a, p.b
    minorder = min(0, order)  # 30

    # compute 0-order basis
    local bas::Matrix{basis_eltype(p,x)}  # stupid type stability...
    if nodetype == 0
        temp = ((n-0.5):-1:0.5)''  # 41
        bas = cos.((pi./n).*temp.*(0:(n-1-minorder))')  # 42
    else
        bas = evalbasex(ChebParams(n-minorder, a, b), x)  # 44
    end

    if order != 0
        D = derivative_op(p, x, order)[1]
        B = bas[:, 1:n-order]*D[abs(order)]
    else
        B = bas
    end

    return B
end

function evalbase(p::ChebParams, x::AbstractArray, order::AbstractVector{Int}, nodetype::Int=1)
    n, a, b = p.n, p.a, p.b
    minorder = min(0, minimum(order))  # 30
    maxorder = maximum(order)

    # compute 0-order basis
    # if nodetype == 0
    #     temp = ((n-0.5):-1:0.5)''  # 41
    #     bas = convert(Matrix{basis_eltype(p, x)},
    #         cos.((pi./n).*temp.*(0:(n-1-minorder))')
    #     )  # 42
    # else
        bas = evalbasex(ChebParams(n-minorder, a, b), x)  # 44
    # end

    B = Array{Matrix{basis_eltype(p, x)}}(length(order))
    if maxorder > 0 D = derivative_op(p, x, maxorder)[1] end
    if minorder < 0 I = derivative_op(p, x, minorder)[1] end

    for ii=1:length(order)
        if order[ii] == 0
            B[ii] = bas[:, 1:n]
        elseif order[ii] > 0
            B[ii] = bas[:, 1:n-order[ii]] * D[order[ii]]
        else
            B[ii] = bas[:, 1:n-order[ii]] * I[-order[ii]]
        end
    end

    return B
end

_unscale(p::ChebParams, x::T) where {T<:Number} = (2/(p.b-p.a)) * (x-(p.a+p.b)/2)

function evalbasex!(out::AbstractMatrix, z::AbstractVector{T},
                    p::ChebParams, x::AbstractVector{T}) where T<:Number
    if size(out) != (size(x, 1), p.n)
        throw(DimensionMismatch("out must be (size(x, 1), p.n)"))
    end

    if size(z) != size(x)
        throw(DimensionMismatch("z must be same size as x"))
    end
    
    if p.n == 1
        out .= 1.0
        return out
    end

    # Note: for julia 0.6+ we can do z .= _unscale.(p, x)
    z .= _unscale.([p], x)
    m = length(z)

    @inbounds out[:, 1] = 1.0
    @inbounds out[:, 2] = z

    scale!(z, 2.0)

    @inbounds for j in 3:p.n
        @simd for i in 1:m
            out[i, j] = z[i] * out[i, j-1] - out[i, j-2]
        end
    end
    out
end

function evalbasex!(out::AbstractMatrix, p::ChebParams, x::AbstractArray)
    z = similar(x)
    evalbasex!(out, z, p, x)
end

function evalbasex(p::ChebParams, x::AbstractArray)
    m = size(x, 1)
    n = p.n
    bas = Array{basis_eltype(p, x)}(m, n)
    evalbasex!(bas, p, x)
end

1	# --------------- #
2	# Chebyshev Basis #
3	# --------------- #
4
5	struct Cheb <: BasisFamily end	5×
6
7	mutable struct ChebParams{T<:Number} <: BasisParams
8	n::Int
9	a::T
10	b::T
11
12	function ChebParams{T}(n::Int, a::T, b::T) where T
13	n <= 0 && error("n must be positive")	165×
14	a >= b && error("left endpoint (a) must be less than right end point (b)")	165×
15	new{T}(n, a, b)	165×
16	end
17	end
18
19	ChebParams(n::Int, a::T, b::T) where {T<:Number} = ChebParams{T}(n, a, b)	165×
20	ChebParams(n::Int, a::T, b::T) where {T<:Integer} = ChebParams(n, Float64(a), Float64(b))	4×
21
22	## BasisParams interface
23	# define these methods on the type, the instance version is defined over
24	# BasisParams
25	family(::Type{T}) where {T<:ChebParams} = Cheb	2×
26	family_name(::Type{T}) where {T<:ChebParams} = "Cheb"	5×
27	@generated Base.eltype(::Type{T}) where {T<:ChebParams} = T.parameters[1]	414×
28
29	# methods that only make sense for instances
30	Base.min(cp::ChebParams) = cp.a	5×
31	Base.max(cp::ChebParams) = cp.b	5×
32	Base.length(cp::ChebParams) = cp.n	9×
33
34	function Base.show(io::IO, p::ChebParams)
35	m = string("Chebyshev interpoland parameters with ",	1×
36	"$(p.n) basis functions from $(p.a), $(p.b)")
37	print(io, m)	1×
38	end
39
40	function nodes(p::ChebParams{T}, ::Type{Val{0}}) where T
41	s = (p.b-p.a) / 2 # 21	23×
42	m = (p.b+p.a) / 2 # 22	23×
43	half = convert(T, 1/2)	23×
44	k = π*(half:(p.n - half)) # 25	23×
45	m .- cos.(k ./ p.n) .* s # 26	23×
46	end
47
48	function nodes(p::ChebParams, ::Type{Val{1}})
49	x = nodes(p, Val{0})	2×
50	aa = x[1]	2×
51	bb = x[end]	2×
52	c1 = (bbp.a - aap.b)/(bb-aa)	2×
53	c2 = (p.b.-p.a)/(bb-aa)	2×
54
55	@inbounds @simd for ix in eachindex(x)	2×
56	x[ix] = c1 + c2 * x[ix]	17×
57	end
58	x	2×
59	end
60
61	function nodes(p::ChebParams{T}, ::Type{Val{2}}) where T
62	s = (p.b-p.a) / 2 # 21	2×
63	m = (p.b+p.a) / 2 # 22	2×
64	k = pi*(zero(T):(p.n - one(T))) # 33	2×
65	m .- cos.(k ./ (p.n-1)) .* s # 34	2×
66	end
67
68	# chebnode.m -- DONE
69	nodes(p::ChebParams, nodetype::Integer=0) = nodes(p, Val{min(2, nodetype)})	46×
70
71	function derivative_op(p::ChebParams, x, order=1)
72	n, a, b = p.n, p.a, p.b	6×
73	if order > 0	6×
74	# TODO: figure out some caching mechanism that avoids globals
75	D = Array{SparseMatrixCSC{basis_eltype(p, x),Int64}}(max(2, order)) # 49	5×
76	i = repmat(1:n', 1, n)	5×
77	j = i' # 50	5×
78
79	# 51
80	inds = find((rem.(i + j, 2) .== 1) .& (j .> i))	5×
81	r, c = similar(inds), similar(inds)	5×
82	for ix in 1:length(inds)	5×
83	r[ix], c[ix] = ind2sub((n, n), inds[ix])	337×
84	end
85
86	d = sparse(r, c, (4/(b-a)) * (vec(j[1, c])-1), n-1, n) # 52	5×
87	d[1, :] ./= 2 # 53	5×
88	D[1] = d # 54	5×
89	for ii in 2:max(2, order)	5×
90	D[ii] = d[1:n-ii, 1:n-ii+1] * D[ii-1] # 56	5×
91	end
92	else
93	D = Array{SparseMatrixCSC{basis_eltype(p, x),Int64}}(abs(order)) # 64	1×
94	nn = n - order # 65	1×
95	z = (0.25 * (b - a)) ./(1:nn) # 66	1×
96	d = sparse(vcat(1:nn, 1:nn-2), vcat(1:nn, 3:nn), vcat(z, -z[1:nn-2]),	1×
97	nn, nn) # 67
98	d[1, 1] *= 2 # 68	1×
99	d0 = ((-1).^(0:nn-1)') .* sum(d, 1) # 69	1×
100	D[1] = sparse(vcat(d0[1:n]', d[1:n, 1:n])) # 70	1×
101	for ii=-2:-1:order	1×
102	ind = 1:n-ii-1	!
103	D[-ii] = sparse(vcat(d0[ind]', d[ind, ind]) * D[-ii-1])	!
104	end
105	end
106	D, ChebParams(n-order, a, b)	6×
107	end
108
109	function evalbase(p::ChebParams, x::AbstractArray=nodes(p, 1), order::Int=0, nodetype::Int=1)
110	n, a, b = p.n, p.a, p.b	253×
111	minorder = min(0, order) # 30	127×
112
113	# compute 0-order basis
114	local bas::Matrix{basis_eltype(p,x)} # stupid type stability...	127×
115	if nodetype == 0	127×
116	temp = ((n-0.5):-1:0.5)'' # 41	1×
117	bas = cos.((pi./n).temp.(0:(n-1-minorder))') # 42	1×
118	else
119	bas = evalbasex(ChebParams(n-minorder, a, b), x) # 44	126×
120	end
121
122	if order != 0	127×
123	D = derivative_op(p, x, order)[1]	4×
124	B = bas[:, 1:n-order]*D[abs(order)]	4×
125	else
126	B = bas	123×
127	end
128
129	return B	127×
130	end
131
132	function evalbase(p::ChebParams, x::AbstractArray, order::AbstractVector{Int}, nodetype::Int=1)
133	n, a, b = p.n, p.a, p.b	!
134	minorder = min(0, minimum(order)) # 30	!
135	maxorder = maximum(order)	!
136
137	# compute 0-order basis
138	# if nodetype == 0
139	# temp = ((n-0.5):-1:0.5)'' # 41
140	# bas = convert(Matrix{basis_eltype(p, x)},
141	# cos.((pi./n).temp.(0:(n-1-minorder))')
142	# ) # 42
143	# else
144	bas = evalbasex(ChebParams(n-minorder, a, b), x) # 44	!
145	# end
146
147	B = Array{Matrix{basis_eltype(p, x)}}(length(order))	!
148	if maxorder > 0 D = derivative_op(p, x, maxorder)[1] end	!
149	if minorder < 0 I = derivative_op(p, x, minorder)[1] end	!
150
151	for ii=1:length(order)	!
152	if order[ii] == 0	!
153	B[ii] = bas[:, 1:n]	!
154	elseif order[ii] > 0	!
155	B[ii] = bas[:, 1:n-order[ii]] * D[order[ii]]	!
156	else
157	B[ii] = bas[:, 1:n-order[ii]] * I[-order[ii]]	!
158	end
159	end
160
161	return B	!
162	end
163
164	_unscale(p::ChebParams, x::T) where {T<:Number} = (2/(p.b-p.a)) * (x-(p.a+p.b)/2)	14,807×
165
166	function evalbasex!(out::AbstractMatrix, z::AbstractVector{T},
167	p::ChebParams, x::AbstractVector{T}) where T<:Number
168	if size(out) != (size(x, 1), p.n)	128×
169	throw(DimensionMismatch("out must be (size(x, 1), p.n)"))	2×
170	end
171
172	if size(z) != size(x)	126×
173	throw(DimensionMismatch("z must be same size as x"))	!
174	end
175
176	if p.n == 1	126×
NEW 177	out .= 1.0	!
NEW 178	return out	!
179	end
180
181	# Note: for julia 0.6+ we can do z .= _unscale.(p, x)
182	z .= _unscale.([p], x)	126×
183	m = length(z)	126×
184
185	@inbounds out[:, 1] = 1.0	126×
186	@inbounds out[:, 2] = z	126×
187
188	scale!(z, 2.0)	126×
189
190	@inbounds for j in 3:p.n	126×
191	@simd for i in 1:m	2,178×
192	out[i, j] = z[i] * out[i, j-1] - out[i, j-2]	187,410×
193	end
194	end
195	out	126×
196	end
197
198	function evalbasex!(out::AbstractMatrix, p::ChebParams, x::AbstractArray)
199	z = similar(x)	126×
200	evalbasex!(out, z, p, x)	126×
201	end
202
203	function evalbasex(p::ChebParams, x::AbstractArray)
204	m = size(x, 1)	126×
205	n = p.n	126×
206	bas = Array{basis_eltype(p, x)}(m, n)	126×
207	evalbasex!(bas, p, x)	126×
208	end

QuantEcon / BasisMatrices.jl / 149

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous