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

97.69

/src/smol_util.jl

# Unrelated to the actual grid, just stuff we use

struct Permuter{T<:AbstractVector}
    a::T
    len::Int
    # sort so we can deal with repeated elements
    Permuter{T}(a::T) where {T} = new{T}(sort(a), length(a))
end
Permuter(a::T) where {T<:AbstractVector} = Permuter{T}(a)

Base.start(p::Permuter) = p.len
Base.done(p::Permuter, i::Int) = i == 0
function Base.next(p::Permuter, i::Int)
    while true
        i -= 1

        if p.a[i] < p.a[i+1]
            j = p.len

            while p.a[i] >= p.a[j]
                j -= 1
            end

            p.a[i], p.a[j] = p.a[j], p.a[i]  # swap(p.a[j], p.a[i])
            t = p.a[(i + 1):end]
            reverse!(t)
            p.a[(i + 1):end] = t

            return copy(p.a), p.len
        end

        if i == 1
            reverse!(p.a)
            return copy(p.a), 0
        end
    end
end

function cartprod(arrs, out=Array{eltype(arrs[1])}(
                                prod([length(a) for a in arrs]),
                                length(arrs))
                                )
    sz = Int[length(a) for a in arrs]
    k = 1
    for v in product(arrs...)
        i = 1
        for el in v
            out[k, i] = el
            i += 1
        end
        k += 1
    end

    return out
end

# Building the actual grid

function m_i(i::Int)
    i < 0 && error("DomainError: i must be positive")
    i == 0 ? 0 : i == 1 ? 1 : 2^(i - 1) + 1
end

function cheby2n(x::AbstractArray{T}, n::Int, kind::Int=1) where T<:Number
    out = Array{T}(size(x)..., n+1)
    cheby2n!(out, x, n, kind)
end

"""
Computes first `n` Chebychev polynomials of the `kind` kind
evaluated at each point in `x` and places them in `out`. The trailing dimension
of `out` indexes the chebyshev polynomials. All inner dimensions correspond to
points in `x`.
"""
function cheby2n!(out::AbstractArray{T}, x::AbstractArray{T,N},
                  n::Int, kind::Int=1) where {T<:Number,N}
    if size(out) != tuple(size(x)..., n+1)
        error("out must have dimensions $(tuple(size(x)..., n+1))")
    end
    
    if n == 0
        out .= one(T)
        return out
    end

    R = CartesianRange(size(x))
    # fill first element with ones
    @inbounds @simd for I in R
        out[I, 1] = one(T)
        out[I, 2] = kind*x[I]
    end

    @inbounds for i in 3:n+1
        @simd for I in R
            out[I, i] = 2x[I] * out[I, i - 1] - out[I, i - 2]
        end
    end
    out
end

"""
Finds the set `S_n` , which is the `n`th Smolyak set of Chebychev extrema
"""
function s_n(n::Int)
    n < 1 && error("DomainError: n must be positive")

    if n == 1
        return [0.0]
    end

    m = m_i(n)
    j = 1:m
    pts = cos.(pi .* (j .- 1.0) ./ (m .- 1.0))
    @inbounds @simd for i in eachindex(pts)
        pts[i] = abs(pts[i]) < 1e-12 ? 0.0 : pts[i]
    end

    pts
end

doc"""
Finds all of the unidimensional disjoint sets of Chebychev extrema that are
used to construct the grid.  It improves on past algorithms by noting  that
$A_{n} = S_{n}$ [evens] except for $A_1= \{0\}$  and $A_2 = \{-1, 1\}$.
Additionally, $A_{n} = A_{n+1}$ [odds] This prevents the calculation of these
nodes repeatedly. Thus we only need to calculate biggest of the S_n's to build
the sequence of $A_n$ 's

See section 3.2 of the paper...
"""
function a_chain(n::Int)
    sn = s_n(n)
    a = Dict{Int,Vector{Float64}}()
    sizehint!(a, n)

    # These are constant and don't follow the pattern.
    a[1] = [0.0]
    a[2] = [-1.0, 1.0]

    for i in n:-1:3
        a[i] = sn[2:2:end]
        sn = sn[1:2:end]
    end

    a
end



doc"""
For each number in 1 to `n`, compute the Smolyak indices for the corresponding
basis functions. This is the $n$ in $\phi_n$. The output is A dictionary whose
keys are the Smolyak index `n` and values are ranges containing all basis
polynomial subscripts for that Smolyak index
"""
function phi_chain(n::Int)
    max_ind = m_i(n)
    phi = Dict{Int, UnitRange{Int64}}()
    phi[1] = 1:1
    phi[2] = 2:3
    low_ind = 4  # current lower index

    for i = 3:n
        high_ind = m_i(i)
        phi[i] = low_ind:high_ind
        low_ind = high_ind + 1
    end

    phi
end

## ---------------------- ##
#- Construction Utilities -#
## ---------------------- ##

doc"""
    smol_inds(d::Int, mu::Int)

Finds all of the indices that satisfy the requirement that $d \leq \sum_{i=1}^d
\leq d + \mu$.
"""
function smol_inds(d::Int, mu::Int)

    p_vals = 1:(mu+1)

    # PERF: size_hint here if it is slow
    poss_inds = Vector{Int}[]

    for el in with_replacement_combinations(p_vals, d)
        if d < sum(el) <= d + mu
            push!(poss_inds, el)
        end
    end

    # PERF: size_hint here if it is slow
    true_inds = Vector{Int}[ones(Int, d)]  # we will always have (1, 1, ...,  1)
    for val in poss_inds
        for el in Permuter(val)
            push!(true_inds, el)
        end
    end

    return true_inds
end

doc"""
    smol_inds(d::Int, mu::AbstractVector{Int})

Finds all of the indices that satisfy the requirement that $d \leq \sum_{i=1}^d
\leq d + \mu_i$.

This is the anisotropic version of the method that allows mu to vary for each
dimension
"""
function smol_inds(d::Int, mu::AbstractVector{Int})
    # Compute indices needed for anisotropic smolyak grid given number of
    # dimensions d and a vector of mu parameters mu

    length(mu) != d &&  error("ValueError: mu must have d elements.")

    mu_max = maximum(mu)
    mup1 = mu + 1

    p_vals = 1:(mu_max+1)

    poss_inds = Vector{Int}[]

    for el in with_replacement_combinations(p_vals, d)
        if d < sum(el) <= d + mu_max
            push!(poss_inds, el)
        end
    end

    true_inds = Vector{Int}[ones(Int64, d)]
    for val in poss_inds
        for el in Permuter(val)
            if all(el .<= mup1)
                push!(true_inds, el)
            end
        end
    end

    return true_inds
end

"""
Build indices specifying all the Cartesian products of Chebychev polynomials
needed to build Smolyak polynomial
"""
function poly_inds(d::Int, mu::IntSorV, inds::Vector{Vector{Int}}=smol_inds(d, mu))::Matrix{Int}
    phi_n = phi_chain(maximum(mu) + 1)
    vcat([cartprod([phi_n[i] for i in el]) for el in inds]...)
end

"""
Use disjoint Smolyak sets to construct Smolyak grid of degree `d` and density
parameter `mu`
"""
function build_grid(d::Int, mu::IntSorV, inds::Vector{Vector{Int}}=smol_inds(d, mu))::Matrix{Float64}
    An = a_chain(maximum(mu) + 1)  # use maximum in case mu is Vector
    vcat([cartprod([An[i] for i in el]) for el in inds]...)
end


"""
Compute the matrix `B(pts)` from equation 22 in JMMV 2013. This is the basis
matrix
"""
function build_B!(out::AbstractMatrix{T}, d::Int, mu::IntSorV,
                  pts::Matrix{Float64}, b_inds::Matrix{Int64}) where T
    # check dimensions
    npolys = size(b_inds, 1)
    npts = size(pts, 1)
    size(out) == (npts, npolys) || error("Out should be size $((npts, npolys))")

    # fill out with ones so tensor product below works
    fill!(out, one(T))

    # compute all the chebyshev polynomials we'll need
    Ts = cheby2n(pts, m_i(maximum(mu) + 1))

    @inbounds for ind in 1:npolys, k in 1:d
        b = b_inds[ind, k]
        for i in 1:npts
            out[i, ind] *= Ts[i, k, b]
        end
    end

    return out
end

function build_B(d::Int, mu::IntSorV, pts::Matrix{Float64}, b_inds::Matrix{Int64})
    build_B!(Array{Float64}(size(pts, 1), size(b_inds, 1)), d, mu, pts, b_inds)
end

function dom2cube!(out::AbstractMatrix, pts::AbstractMatrix,
                   lb::AbstractVector, ub::AbstractVector)
    d = length(lb)
    n = size(pts, 1)

    size(out) == (n, d) || error("out should be $((n, d))")

    @inbounds for i_d in 1:d
        center = lb[i_d] + (ub[i_d] - lb[i_d])/2
        radius = (ub[i_d] - lb[i_d])/2
        @simd for i_n in 1:n
            out[i_n, i_d] = (pts[i_n, i_d] - center)/radius
        end
    end

    out
end

function cube2dom!(out::AbstractMatrix, pts::AbstractMatrix,
                   lb::AbstractVector, ub::AbstractVector)
    d = length(lb)
    n = size(pts, 1)

    size(out) == (n, d) || error("out should be $((n, d))")

    @inbounds for i_d in 1:d
        center = lb[i_d] + (ub[i_d] - lb[i_d])/2
        radius = (ub[i_d] - lb[i_d])/2
        @simd for i_n in 1:n
            out[i_n, i_d] = center + pts[i_n, i_d]*radius
        end
    end

    out
end

for f in [:dom2cube!, :cube2dom!]
    no_bang = Symbol(string(f)[1:end-1])
    @eval $(no_bang)(pts::AbstractMatrix{T}, lb::AbstractVector, ub::AbstractVector) where {T} =
        $(f)(Array{T}(size(pts, 1), length(lb)), pts, lb, ub)
end

1	# Unrelated to the actual grid, just stuff we use
2
3	struct Permuter{T<:AbstractVector}
4	a::T
5	len::Int
6	# sort so we can deal with repeated elements
7	Permuter{T}(a::T) where {T} = new{T}(sort(a), length(a))	482×
8	end
9	Permuter(a::T) where {T<:AbstractVector} = Permuter{T}(a)	482×
10
11	Base.start(p::Permuter) = p.len	483×
12	Base.done(p::Permuter, i::Int) = i == 0	4,042,579×
13	function Base.next(p::Permuter, i::Int)
14	while true	4,042,097×
15	i -= 1	6,951,229×
16
17	if p.a[i] < p.a[i+1]	6,951,229×
18	j = p.len	4,041,615×
19
20	while p.a[i] >= p.a[j]	4,041,615×
21	j -= 1	1,457,038×
22	end
23
24	p.a[i], p.a[j] = p.a[j], p.a[i] # swap(p.a[j], p.a[i])	4,041,615×
25	t = p.a[(i + 1):end]	4,041,615×
26	reverse!(t)	4,041,615×
27	p.a[(i + 1):end] = t	4,041,615×
28
29	return copy(p.a), p.len	4,041,615×
30	end
31
32	if i == 1	2,909,614×
33	reverse!(p.a)	482×
34	return copy(p.a), 0	482×
35	end
36	end
37	end
38
39	function cartprod(arrs, out=Array{eltype(arrs[1])}(
40	prod([length(a) for a in arrs]),
41	length(arrs))
42	)
43	sz = Int[length(a) for a in arrs]	3,030×
44	k = 1	1,515×
45	for v in product(arrs...)	1,515×
46	i = 1	36,041×
47	for el in v	36,041×
48	out[k, i] = el	194,829×
49	i += 1	194,829×
50	end
51	k += 1	36,041×
52	end
53
54	return out	1,515×
55	end
56
57	# Building the actual grid
58
59	function m_i(i::Int)
60	i < 0 && error("DomainError: i must be positive")	286×
61	i == 0 ? 0 : i == 1 ? 1 : 2^(i - 1) + 1	276×
62	end
63
64	function cheby2n(x::AbstractArray{T}, n::Int, kind::Int=1) where T<:Number
65	out = Array{T}(size(x)..., n+1)	194×
66	cheby2n!(out, x, n, kind)	97×
67	end
68
69	"""
70	Computes first `n` Chebychev polynomials of the `kind` kind
71	evaluated at each point in `x` and places them in `out`. The trailing dimension
72	of `out` indexes the chebyshev polynomials. All inner dimensions correspond to
73	points in `x`.
74	"""
75	function cheby2n!(out::AbstractArray{T}, x::AbstractArray{T,N},
76	n::Int, kind::Int=1) where {T<:Number,N}
77	if size(out) != tuple(size(x)..., n+1)	116×
78	error("out must have dimensions $(tuple(size(x)..., n+1))")	!
79	end
80
81	if n == 0	116×
NEW 82	out .= one(T)	!
NEW 83	return out	!
84	end
85
86	R = CartesianRange(size(x))	116×
87	# fill first element with ones
88	@inbounds @simd for I in R	116×
89	out[I, 1] = one(T)	115,595×
90	out[I, 2] = kind*x[I]	115,595×
91	end
92
93	@inbounds for i in 3:n+1	116×
94	@simd for I in R	1,018×
95	out[I, i] = 2x[I] * out[I, i - 1] - out[I, i - 2]	561,440×
96	end
97	end
98	out	116×
99	end
100
101	"""
102	Finds the set `S_n` , which is the `n`th Smolyak set of Chebychev extrema
103	"""
104	function s_n(n::Int)
105	n < 1 && error("DomainError: n must be positive")	61×
106
107	if n == 1	61×
108	return [0.0]	2×
109	end
110
111	m = m_i(n)	59×
112	j = 1:m	59×
113	pts = cos.(pi .* (j .- 1.0) ./ (m .- 1.0))	59×
114	@inbounds @simd for i in eachindex(pts)	59×
115	pts[i] = abs(pts[i]) < 1e-12 ? 0.0 : pts[i]	1,483×
116	end
117
118	pts	59×
119	end
120
121	doc"""
122	Finds all of the unidimensional disjoint sets of Chebychev extrema that are
123	used to construct the grid. It improves on past algorithms by noting that
124	$A_{n} = S_{n}$ [evens] except for $A_1= \{0\}$ and $A_2 = \{-1, 1\}$.
125	Additionally, $A_{n} = A_{n+1}$ [odds] This prevents the calculation of these
126	nodes repeatedly. Thus we only need to calculate biggest of the S_n's to build
127	the sequence of $A_n$ 's
128
129	See section 3.2 of the paper...
130	"""
131	function a_chain(n::Int)
132	sn = s_n(n)	45×
133	a = Dict{Int,Vector{Float64}}()	45×
134	sizehint!(a, n)	45×
135
136	# These are constant and don't follow the pattern.
137	a[1] = [0.0]	45×
138	a[2] = [-1.0, 1.0]	45×
139
140	for i in n:-1:3	45×
141	a[i] = sn[2:2:end]	73×
142	sn = sn[1:2:end]	73×
143	end
144
145	a	45×
146	end
147
148
149
150	doc"""
151	For each number in 1 to `n`, compute the Smolyak indices for the corresponding
152	basis functions. This is the $n$ in $\phi_n$. The output is A dictionary whose
153	keys are the Smolyak index `n` and values are ranges containing all basis
154	polynomial subscripts for that Smolyak index
155	"""
156	function phi_chain(n::Int)
157	max_ind = m_i(n)	48×
158	phi = Dict{Int, UnitRange{Int64}}()	48×
159	phi[1] = 1:1	48×
160	phi[2] = 2:3	48×
161	low_ind = 4 # current lower index	48×
162
163	for i = 3:n	48×
164	high_ind = m_i(i)	99×
165	phi[i] = low_ind:high_ind	99×
166	low_ind = high_ind + 1	99×
167	end
168
169	phi	48×
170	end
171
172	## ---------------------- ##
173	#- Construction Utilities -#
174	## ---------------------- ##
175
176	doc"""
177	smol_inds(d::Int, mu::Int)
178
179	Finds all of the indices that satisfy the requirement that $d \leq \sum_{i=1}^d
180	\leq d + \mu$.
181	"""
182	function smol_inds(d::Int, mu::Int)
183
184	p_vals = 1:(mu+1)	83×
185
186	# PERF: size_hint here if it is slow
187	poss_inds = Vector{Int}[]	83×
188
189	for el in with_replacement_combinations(p_vals, d)	83×
190	if d < sum(el) <= d + mu	2,584×
191	push!(poss_inds, el)	290×
192	end
193	end
194
195	# PERF: size_hint here if it is slow
196	true_inds = Vector{Int}[ones(Int, d)] # we will always have (1, 1, ..., 1)	83×
197	for val in poss_inds	83×
198	for el in Permuter(val)	290×
199	push!(true_inds, el)	2,501×
200	end
201	end
202
203	return true_inds	83×
204	end
205
206	doc"""
207	smol_inds(d::Int, mu::AbstractVector{Int})
208
209	Finds all of the indices that satisfy the requirement that $d \leq \sum_{i=1}^d
210	\leq d + \mu_i$.
211
212	This is the anisotropic version of the method that allows mu to vary for each
213	dimension
214	"""
215	function smol_inds(d::Int, mu::AbstractVector{Int})
216	# Compute indices needed for anisotropic smolyak grid given number of
217	# dimensions d and a vector of mu parameters mu
218
219	length(mu) != d && error("ValueError: mu must have d elements.")	43×
220
221	mu_max = maximum(mu)	40×
222	mup1 = mu + 1	40×
223
224	p_vals = 1:(mu_max+1)	40×
225
226	poss_inds = Vector{Int}[]	40×
227
228	for el in with_replacement_combinations(p_vals, d)	40×
229	if d < sum(el) <= d + mu_max	1,716×
230	push!(poss_inds, el)	181×
231	end
232	end
233
234	true_inds = Vector{Int}[ones(Int64, d)]	40×
235	for val in poss_inds	40×
236	for el in Permuter(val)	181×
237	if all(el .<= mup1)	1,676×
238	push!(true_inds, el)	1,611×
239	end
240	end
241	end
242
243	return true_inds	40×
244	end
245
246	"""
247	Build indices specifying all the Cartesian products of Chebychev polynomials
248	needed to build Smolyak polynomial
249	"""
250	function poly_inds(d::Int, mu::IntSorV, inds::Vector{Vector{Int}}=smol_inds(d, mu))::Matrix{Int}
251	phi_n = phi_chain(maximum(mu) + 1)	64×
252	vcat([cartprod([phi_n[i] for i in el]) for el in inds]...)	37×
253	end
254
255	"""
256	Use disjoint Smolyak sets to construct Smolyak grid of degree `d` and density
257	parameter `mu`
258	"""
259	function build_grid(d::Int, mu::IntSorV, inds::Vector{Vector{Int}}=smol_inds(d, mu))::Matrix{Float64}
260	An = a_chain(maximum(mu) + 1) # use maximum in case mu is Vector	63×
261	vcat([cartprod([An[i] for i in el]) for el in inds]...)	39×
262	end
263
264
265	"""
266	Compute the matrix `B(pts)` from equation 22 in JMMV 2013. This is the basis
267	matrix
268	"""
269	function build_B!(out::AbstractMatrix{T}, d::Int, mu::IntSorV,
270	pts::Matrix{Float64}, b_inds::Matrix{Int64}) where T
271	# check dimensions
272	npolys = size(b_inds, 1)	59×
273	npts = size(pts, 1)	59×
274	size(out) == (npts, npolys) \|\| error("Out should be size $((npts, npolys))")	59×
275
276	# fill out with ones so tensor product below works
277	fill!(out, one(T))	59×
278
279	# compute all the chebyshev polynomials we'll need
280	Ts = cheby2n(pts, m_i(maximum(mu) + 1))	59×
281
282	@inbounds for ind in 1:npolys, k in 1:d	59×
283	b = b_inds[ind, k]	8,829×
284	for i in 1:npts	8,829×
285	out[i, ind] *= Ts[i, k, b]	10,972,375×
286	end
287	end
288
289	return out	59×
290	end
291
292	function build_B(d::Int, mu::IntSorV, pts::Matrix{Float64}, b_inds::Matrix{Int64})
293	build_B!(Array{Float64}(size(pts, 1), size(b_inds, 1)), d, mu, pts, b_inds)	27×
294	end
295
296	function dom2cube!(out::AbstractMatrix, pts::AbstractMatrix,
297	lb::AbstractVector, ub::AbstractVector)
298	d = length(lb)	27×
299	n = size(pts, 1)	27×
300
301	size(out) == (n, d) \|\| error("out should be $((n, d))")	27×
302
303	@inbounds for i_d in 1:d	27×
304	center = lb[i_d] + (ub[i_d] - lb[i_d])/2	57×
305	radius = (ub[i_d] - lb[i_d])/2	57×
306	@simd for i_n in 1:n	57×
307	out[i_n, i_d] = (pts[i_n, i_d] - center)/radius	2,518×
308	end
309	end
310
311	out	27×
312	end
313
314	function cube2dom!(out::AbstractMatrix, pts::AbstractMatrix,
315	lb::AbstractVector, ub::AbstractVector)
316	d = length(lb)	32×
317	n = size(pts, 1)	32×
318
319	size(out) == (n, d) \|\| error("out should be $((n, d))")	32×
320
321	@inbounds for i_d in 1:d	32×
322	center = lb[i_d] + (ub[i_d] - lb[i_d])/2	68×
323	radius = (ub[i_d] - lb[i_d])/2	68×
324	@simd for i_n in 1:n	68×
325	out[i_n, i_d] = center + pts[i_n, i_d]*radius	3,005×
326	end
327	end
328
329	out	32×
330	end
331
332	for f in [:dom2cube!, :cube2dom!]
333	no_bang = Symbol(string(f)[1:end-1])
334	@eval $(no_bang)(pts::AbstractMatrix{T}, lb::AbstractVector, ub::AbstractVector) where {T} =	32×
335	$(f)(Array{T}(size(pts, 1), length(lb)), pts, lb, ub)
336	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