4659511182

Build Type

github

Committed by GitHub

Commit Message

Merge 3320693ee into 52de2ab68

Pull Request Pull Request #108: Sphinx build

Run Details

1589 of 3151 relevant lines covered (50.43%)

4.46 hits per line

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

66.67

/UncertainSCI/utils/quad.py

import numpy as np

from UncertainSCI.families import jacobi_recurrence_values,\
                                  jacobi_weight_normalized

from UncertainSCI.opoly1d import linear_modification, quadratic_modification
from UncertainSCI.opoly1d import gauss_quadrature_driver

# from UncertainSCI.utils.compute_subintervals import compute_subintervals


def compute_subintervals(a, b, singularity_list):
    """
    Returns an M x 4 numpy array, where each row contains the left-hand point,
    right-hand point, left-singularity strength, and right-singularity
    strength.
    """

    # Tolerance for resolving internal versus boundary singularities.
    tol = 1e-12
    singularities = np.array([entry[0] for entry in singularity_list])
    strength_left = np.array([entry[1] for entry in singularity_list])
    strength_right = np.array([entry[2] for entry in singularity_list])

    # We can discard any singularities that lie to the left of a or the right
    # of b
    discard = []
    for (ind, s) in enumerate(singularities):
        if s < a-tol or s > b+tol:
            discard.append(ind)

    singularities = np.delete(singularities, discard)
    strength_left = np.delete(strength_left, discard)
    strength_right = np.delete(strength_right, discard)

    # Sort remaining valid singularities
    order = np.argsort(singularities)
    singularities = singularities[order]
    strength_left = strength_left[order]
    strength_right = strength_right[order]

    # Make sure there aren't doubly-specified singularities
    if np.any(np.diff(singularities) < tol):
        raise ValueError("Overlapping singularities were specified. \
                          Singularities must be unique")

    S = singularities.size

    if S > 0:

        # Extend the singularities lists if we need to add interval endpoints
        a_sing = np.abs(singularities[0] - a) <= tol
        b_sing = np.abs(singularities[-1] - b) <= tol

        # Figure out if singularities match endpoints
        if not b_sing:
            singularities = np.hstack([singularities, b])
            strength_left = np.hstack([strength_left, 0])
            strength_right = np.hstack([strength_right, 0])  # Doesn't matter
        if not a_sing:
            singularities = np.hstack([a, singularities])
            strength_left = np.hstack([0, strength_left])  # Doesn't matter
            strength_right = np.hstack([0, strength_right])  # Doesn't matter

        # Use the singularities lists to identify subintervals
        S = singularities.size
        subintervals = np.zeros([S-1, 4])
        for q in range(S-1):
            subintervals[q, :] = [singularities[q], singularities[q+1],
                                  strength_right[q], strength_left[q+1]]

    else:

        subintervals = np.zeros([1, 4])
        subintervals[0, :] = [a, b, 0, 0]

    return subintervals


def gq_modification(integrand, a, b, N, roots=np.zeros(0),
                    quadroots=np.zeros(0), Nmax=100,
                    gamma=(0, 0), leading_coefficient=1.):
    """
    Uses Gaussian quadrature to approximate

      \\int_a^b q(x) integrand(x) dx,

    where integrand is an input function, and q is a polynomial defined as,

      q(x) = leading_coefficient * prod_{j=1}^R (x - roots[j-1]) *
                                   prod_{q=1}^Q (x - quadroots[j-1])**2,

    where R = roots.size and Q = quadroots.size, and it is assumed that no
    entries of roots lie inside the interval (a,b).  The coefficient
    leading_coefficient can be negative.

    The procedure rewrites the integral using the input gamma as

      \\int_a^b (integrand(x) / w(x)) q(x) w(x) dx,

    where w is a mapped Jacobi weight function:

        w(x) = v(r(x)),    gamma = (alpha, beta)

    where r(x) affinely maps (a,b) to (-1,1), and v(r) is the (alpha,beta)
    Jacobi weight function in UncertainSCI.families.jacobi_weight_normalized.

    The procedure then performs measure modifications on w, absorbing q into
    the measure. An N-point Gaussian quadrature rule with respect to this
    modified measure is used to integrate (integrand / w).
    """

    assert (a < b) and (N > 0)

    # If the gamma parameters are > 0, then we write
    #   gamma[0] = gam + G,
    # where gam \in (-1,0) and G is a positive integer. And similarly for
    # gamma[1].
    #
    # Then we take the Jacobi weight to be w associated with gam,
    # and set N, Nmax += G
    for ind in range(2):
        if gamma[ind] > 0:
            G = np.ceil(gamma[ind])
            gam = gamma[ind] - G
            N += int(G)
            Nmax += int(G)
            gamma[ind] = gam

    assert (gamma[0] <= 0.) and (gamma[1] <= 0.)
    Nmax = max(Nmax, N)

    R = roots.size
    Q = quadroots.size

    # Map everything to [-1,1], and scale by Jacobian afterward
    jac = 2/(b-a)

    def map_to_standard(x):  # Maps [a,b] to [-1,1]
        return jac*(x - a) - 1
    # map_to_standard = lambda x: jac*(x - a) - 1

    def map_to_ab(x):  # Maps [-1,1], to [a,b]
        return (x + 1) / jac + a
    # map_to_ab = lambda x: (x+1)/jac + a

    # Recurrence coefficients for appropriate Jacobi probability measure
    ab = jacobi_recurrence_values(Nmax+R+2*Q, gamma[0], gamma[1])
    ab[0, 1] = 1.

    # The sign of q is determined by how many zeros lie to the right of (a,b)
    sgn = (-1)**(np.sum(roots >= b) % 2) if R > 0 else 1.

    # Modify ab by the zeros of q
    # Simultaneously gradually scale by leading_coefficient and jacobians
    sgn *= np.sign(leading_coefficient)
    if (R+2*Q) > 0:
        # Scale each linear modification by C = leading_coefficient**(1/(R+2Q))
        # There is also a jac factor since (x - x0) = 1/jac * (r-r0)
        # Finally, we need the square root of this because ab[0,1] scales like
        # the square root of an integral.
        C = np.exp(np.log(np.abs(leading_coefficient))/(R+2*Q))/jac
        Csqrt = np.sqrt(C)

        for x0 in quadroots:
            ab = quadratic_modification(ab, map_to_standard(x0))
            ab[0, 1] *= C
        for x0 in roots:
            ab = linear_modification(ab, map_to_standard(x0))
            ab[0, 1] *= Csqrt
    else:
        ab[0, 1] *= np.abs(leading_coefficient)

    x, w = gauss_quadrature_driver(ab, N)

    # Gauss quadrature
    integral = np.sum(w * integrand(map_to_ab(x)) /
                      jacobi_weight_normalized(x, gamma[0], gamma[1]))

    # Jacobian factor for the affine map and sign:
    integral *= sgn/jac

    return integral


def gq_modification_adaptive(integrand, a, b, N, N_step=10, tol=1e-12,
                             **kwargs):
    s = gq_modification(integrand, a, b, N, **kwargs)
    s_new = gq_modification(integrand, a, b, N=N + N_step, **kwargs)

    while np.abs(s - s_new) > tol:
        s = s_new
        N += N_step
        s_new = gq_modification(integrand, a, b, N=N, **kwargs)
    return s_new


def gq_modification_composite(integrand, a, b, N,
                              subintervals=np.zeros([0, 4]), adaptive=True,
                              **kwargs):
    """
    Uses a composite quadrature rule where each subinterval uses
    gq_modification to integrate. The integral is split into
    subintervals.shape[0] integrals.

    subintervals is an M x 4 array, where each row has the form

      [left, right, left_sing, right_sing]

    with each subinterval being [left, right], with the appropriate left- and
    right-hand singularity strengths.

    If subintervals is empty, performs integration over an interval specified
    by

      [a, b, 0, 0]

    Typical keyword arguments that should be input to this function are:

        - quadroots
        - roots
        - leading_coefficient

    See gq_modification for a description of these inputs.
    """
    integral = 0.

    if subintervals.shape[0] == 0:
        subintervals = np.zeros([1, 4])
        subintervals[0, :] = [a, b, 0, 0]

    for q in range(subintervals.shape[0]):
        gamma = [subintervals[q, 3], subintervals[q, 2]]
        if adaptive:
            integral += gq_modification_adaptive(integrand, subintervals[q, 0],
                                                 subintervals[q, 1], N,
                                                 gamma=gamma, **kwargs)
        else:
            integral += gq_modification(integrand, subintervals[q, 0],
                                        subintervals[q, 1], N, gamma=gamma,
                                        **kwargs)

    return integral


def gq_modification_unbounded_composite(integrand, a, b, N, singularity_list,
                                        adaptive=True, tol=1e-12, step=1,
                                        **kwargs):
    """
    Functional as the same as util.quad.gq_modification_composite but with a
    unbounded interval [a,b] I expect this method lying in the util.quad, but
    the issue is:

    1. Subintervals are changed when extending the intervals that we're
    integrating on, thus subintervals cannot be a argument, instead, We use
    singularity_list.

    2. Method compute_subintervals in the composite is required to obtain the
    subintervals for different interval [a,b]

    So We temporarily put this method here, maybe change this later for bettter
    consideration.
    """

    if a == -np.inf and b == np.inf:
        # l = -1.; r = 1.
        le, r = -1., 1.
        subintervals = compute_subintervals(le, r, singularity_list)
        integral = gq_modification_composite(integrand, le, r, N, subintervals,
                                             adaptive, **kwargs)

        integral_new = 1.
        while np.abs(integral_new) > tol:
            r = le
            le = r - step
            subintervals = compute_subintervals(le, r, singularity_list)
            integral_new = gq_modification_composite(integrand, le, r, N,
                                                     subintervals, adaptive,
                                                     **kwargs)
            integral += integral_new

        le, r = -1., 1.
        # l = -1.; r = 1.
        integral_new = 1.
        while np.abs(integral_new) > tol:
            le = r
            r = le + step
            subintervals = compute_subintervals(le, r, singularity_list)
            integral_new = gq_modification_composite(integrand, le, r, N,
                                                     subintervals, adaptive,
                                                     **kwargs)
            integral += integral_new

    elif a == -np.inf:
        r = b
        le = b - step
        subintervals = compute_subintervals(le, r, singularity_list)
        integral = gq_modification_composite(integrand, le, r, N, subintervals,
                                             adaptive, **kwargs)
        integral_new = 1.
        while np.abs(integral_new) > tol:
            r = le
            le = r - step
            subintervals = compute_subintervals(le, r, singularity_list)
            integral_new = gq_modification_composite(integrand, le, r, N,
                                                     subintervals, adaptive,
                                                     **kwargs)
            integral += integral_new

    elif b == np.inf:
        le = a
        r = a + step
        subintervals = compute_subintervals(le, r, singularity_list)
        integral = gq_modification_composite(integrand, le, r, N, subintervals,
                                             adaptive, **kwargs)
        integral_new = 1.
        while np.abs(integral_new) > tol:
            le = r
            r = le + step
            subintervals = compute_subintervals(le, r, singularity_list)
            integral_new = gq_modification_composite(integrand, le, r, N,
                                                     subintervals, adaptive,
                                                     **kwargs)
            integral += integral_new

    return integral

1	import numpy as np	8✔
2
3	from UncertainSCI.families import jacobi_recurrence_values,\	8✔
4	jacobi_weight_normalized
5
6	from UncertainSCI.opoly1d import linear_modification, quadratic_modification	8✔
7	from UncertainSCI.opoly1d import gauss_quadrature_driver	8✔
8
9	# from UncertainSCI.utils.compute_subintervals import compute_subintervals
10
11
12	def compute_subintervals(a, b, singularity_list):	8✔
13	"""
14	Returns an M x 4 numpy array, where each row contains the left-hand point,
15	right-hand point, left-singularity strength, and right-singularity
16	strength.
17	"""
18
19	# Tolerance for resolving internal versus boundary singularities.
20	tol = 1e-12	8✔
21	singularities = np.array([entry[0] for entry in singularity_list])	8✔
22	strength_left = np.array([entry[1] for entry in singularity_list])	8✔
23	strength_right = np.array([entry[2] for entry in singularity_list])	8✔
24
25	# We can discard any singularities that lie to the left of a or the right
26	# of b
27	discard = []	8✔
28	for (ind, s) in enumerate(singularities):	8✔
29	if s < a-tol or s > b+tol:	×
30	discard.append(ind)	×
31
32	singularities = np.delete(singularities, discard)	8✔
33	strength_left = np.delete(strength_left, discard)	8✔
34	strength_right = np.delete(strength_right, discard)	8✔
35
36	# Sort remaining valid singularities
37	order = np.argsort(singularities)	8✔
38	singularities = singularities[order]	8✔
39	strength_left = strength_left[order]	8✔
40	strength_right = strength_right[order]	8✔
41
42	# Make sure there aren't doubly-specified singularities
43	if np.any(np.diff(singularities) < tol):	8✔
44	raise ValueError("Overlapping singularities were specified. \	×
45	Singularities must be unique")
46
47	S = singularities.size	8✔
48
49	if S > 0:	8✔
50
51	# Extend the singularities lists if we need to add interval endpoints
52	a_sing = np.abs(singularities[0] - a) <= tol	×
53	b_sing = np.abs(singularities[-1] - b) <= tol	×
54
55	# Figure out if singularities match endpoints
56	if not b_sing:	×
57	singularities = np.hstack([singularities, b])	×
58	strength_left = np.hstack([strength_left, 0])	×
59	strength_right = np.hstack([strength_right, 0]) # Doesn't matter	×
60	if not a_sing:	×
61	singularities = np.hstack([a, singularities])	×
62	strength_left = np.hstack([0, strength_left]) # Doesn't matter	×
63	strength_right = np.hstack([0, strength_right]) # Doesn't matter	×
64
65	# Use the singularities lists to identify subintervals
66	S = singularities.size	×
67	subintervals = np.zeros([S-1, 4])	×
68	for q in range(S-1):	×
69	subintervals[q, :] = [singularities[q], singularities[q+1],	×
70	strength_right[q], strength_left[q+1]]
71
72	else:
73
74	subintervals = np.zeros([1, 4])	8✔
75	subintervals[0, :] = [a, b, 0, 0]	8✔
76
77	return subintervals	8✔
78
79
80	def gq_modification(integrand, a, b, N, roots=np.zeros(0),	8✔
81	quadroots=np.zeros(0), Nmax=100,
82	gamma=(0, 0), leading_coefficient=1.):
83	"""
84	Uses Gaussian quadrature to approximate
85
86	\\int_a^b q(x) integrand(x) dx,
87
88	where integrand is an input function, and q is a polynomial defined as,
89
90	q(x) = leading_coefficient * prod_{j=1}^R (x - roots[j-1]) *
91	prod_{q=1}^Q (x - quadroots[j-1])**2,
92
93	where R = roots.size and Q = quadroots.size, and it is assumed that no
94	entries of roots lie inside the interval (a,b). The coefficient
95	leading_coefficient can be negative.
96
97	The procedure rewrites the integral using the input gamma as
98
99	\\int_a^b (integrand(x) / w(x)) q(x) w(x) dx,
100
101	where w is a mapped Jacobi weight function:
102
103	w(x) = v(r(x)), gamma = (alpha, beta)
104
105	where r(x) affinely maps (a,b) to (-1,1), and v(r) is the (alpha,beta)
106	Jacobi weight function in UncertainSCI.families.jacobi_weight_normalized.
107
108	The procedure then performs measure modifications on w, absorbing q into
109	the measure. An N-point Gaussian quadrature rule with respect to this
110	modified measure is used to integrate (integrand / w).
111	"""
112
113	assert (a < b) and (N > 0)	8✔
114
115	# If the gamma parameters are > 0, then we write
116	# gamma[0] = gam + G,
117	# where gam \in (-1,0) and G is a positive integer. And similarly for
118	# gamma[1].
119	#
120	# Then we take the Jacobi weight to be w associated with gam,
121	# and set N, Nmax += G
122	for ind in range(2):	8✔
123	if gamma[ind] > 0:	8✔
124	G = np.ceil(gamma[ind])	8✔
125	gam = gamma[ind] - G	8✔
126	N += int(G)	8✔
127	Nmax += int(G)	8✔
128	gamma[ind] = gam	8✔
129
130	assert (gamma[0] <= 0.) and (gamma[1] <= 0.)	8✔
131	Nmax = max(Nmax, N)	8✔
132
133	R = roots.size	8✔
134	Q = quadroots.size	8✔
135
136	# Map everything to [-1,1], and scale by Jacobian afterward
137	jac = 2/(b-a)	8✔
138
139	def map_to_standard(x): # Maps [a,b] to [-1,1]	8✔
140	return jac*(x - a) - 1	8✔
141	# map_to_standard = lambda x: jac*(x - a) - 1
142
143	def map_to_ab(x): # Maps [-1,1], to [a,b]	8✔
144	return (x + 1) / jac + a	8✔
145	# map_to_ab = lambda x: (x+1)/jac + a
146
147	# Recurrence coefficients for appropriate Jacobi probability measure
148	ab = jacobi_recurrence_values(Nmax+R+2*Q, gamma[0], gamma[1])	8✔
149	ab[0, 1] = 1.	8✔
150
151	# The sign of q is determined by how many zeros lie to the right of (a,b)
152	sgn = (-1)**(np.sum(roots >= b) % 2) if R > 0 else 1.	8✔
153
154	# Modify ab by the zeros of q
155	# Simultaneously gradually scale by leading_coefficient and jacobians
156	sgn *= np.sign(leading_coefficient)	8✔
157	if (R+2*Q) > 0:	8✔
158	# Scale each linear modification by C = leading_coefficient**(1/(R+2Q))
159	# There is also a jac factor since (x - x0) = 1/jac * (r-r0)
160	# Finally, we need the square root of this because ab[0,1] scales like
161	# the square root of an integral.
162	C = np.exp(np.log(np.abs(leading_coefficient))/(R+2*Q))/jac	8✔
163	Csqrt = np.sqrt(C)	8✔
164
165	for x0 in quadroots:	8✔
166	ab = quadratic_modification(ab, map_to_standard(x0))	8✔
167	ab[0, 1] *= C	8✔
168	for x0 in roots:	8✔
169	ab = linear_modification(ab, map_to_standard(x0))	8✔
170	ab[0, 1] *= Csqrt	8✔
171	else:
172	ab[0, 1] *= np.abs(leading_coefficient)	8✔
173
174	x, w = gauss_quadrature_driver(ab, N)	8✔
175
176	# Gauss quadrature
177	integral = np.sum(w * integrand(map_to_ab(x)) /	8✔
178	jacobi_weight_normalized(x, gamma[0], gamma[1]))
179
180	# Jacobian factor for the affine map and sign:
181	integral *= sgn/jac	8✔
182
183	return integral	8✔
184
185
186	def gq_modification_adaptive(integrand, a, b, N, N_step=10, tol=1e-12,	8✔
187	**kwargs):
188	s = gq_modification(integrand, a, b, N, **kwargs)	8✔
189	s_new = gq_modification(integrand, a, b, N=N + N_step, **kwargs)	8✔
190
191	while np.abs(s - s_new) > tol:	8✔
192	s = s_new	×
193	N += N_step	×
194	s_new = gq_modification(integrand, a, b, N=N, **kwargs)	×
195	return s_new	8✔
196
197
198	def gq_modification_composite(integrand, a, b, N,	8✔
199	subintervals=np.zeros([0, 4]), adaptive=True,
200	**kwargs):
201	"""
202	Uses a composite quadrature rule where each subinterval uses
203	gq_modification to integrate. The integral is split into
204	subintervals.shape[0] integrals.
205
206	subintervals is an M x 4 array, where each row has the form
207
208	[left, right, left_sing, right_sing]
209
210	with each subinterval being [left, right], with the appropriate left- and
211	right-hand singularity strengths.
212
213	If subintervals is empty, performs integration over an interval specified
214	by
215
216	[a, b, 0, 0]
217
218	Typical keyword arguments that should be input to this function are:
219
220	- quadroots
221	- roots
222	- leading_coefficient
223
224	See gq_modification for a description of these inputs.
225	"""
226	integral = 0.	8✔
227
228	if subintervals.shape[0] == 0:	8✔
229	subintervals = np.zeros([1, 4])	×
230	subintervals[0, :] = [a, b, 0, 0]	×
231
232	for q in range(subintervals.shape[0]):	8✔
233	gamma = [subintervals[q, 3], subintervals[q, 2]]	8✔
234	if adaptive:	8✔
235	integral += gq_modification_adaptive(integrand, subintervals[q, 0],	8✔
236	subintervals[q, 1], N,
237	gamma=gamma, **kwargs)
238	else:
239	integral += gq_modification(integrand, subintervals[q, 0],	×
240	subintervals[q, 1], N, gamma=gamma,
241	**kwargs)
242
243	return integral	8✔
244
245
246	def gq_modification_unbounded_composite(integrand, a, b, N, singularity_list,	8✔
247	adaptive=True, tol=1e-12, step=1,
248	**kwargs):
249	"""
250	Functional as the same as util.quad.gq_modification_composite but with a
251	unbounded interval [a,b] I expect this method lying in the util.quad, but
252	the issue is:
253
254	1. Subintervals are changed when extending the intervals that we're
255	integrating on, thus subintervals cannot be a argument, instead, We use
256	singularity_list.
257
258	2. Method compute_subintervals in the composite is required to obtain the
259	subintervals for different interval [a,b]
260
261	So We temporarily put this method here, maybe change this later for bettter
262	consideration.
263	"""
264
265	if a == -np.inf and b == np.inf:	8✔
266	# l = -1.; r = 1.
267	le, r = -1., 1.	8✔
268	subintervals = compute_subintervals(le, r, singularity_list)	8✔
269	integral = gq_modification_composite(integrand, le, r, N, subintervals,	8✔
270	adaptive, **kwargs)
271
272	integral_new = 1.	8✔
273	while np.abs(integral_new) > tol:	8✔
274	r = le	8✔
275	le = r - step	8✔
276	subintervals = compute_subintervals(le, r, singularity_list)	8✔
277	integral_new = gq_modification_composite(integrand, le, r, N,	8✔
278	subintervals, adaptive,
279	**kwargs)
280	integral += integral_new	8✔
281
282	le, r = -1., 1.	8✔
283	# l = -1.; r = 1.
284	integral_new = 1.	8✔
285	while np.abs(integral_new) > tol:	8✔
286	le = r	8✔
287	r = le + step	8✔
288	subintervals = compute_subintervals(le, r, singularity_list)	8✔
289	integral_new = gq_modification_composite(integrand, le, r, N,	8✔
290	subintervals, adaptive,
291	**kwargs)
292	integral += integral_new	8✔
293
294	elif a == -np.inf:	×
295	r = b	×
296	le = b - step	×
297	subintervals = compute_subintervals(le, r, singularity_list)	×
298	integral = gq_modification_composite(integrand, le, r, N, subintervals,	×
299	adaptive, **kwargs)
300	integral_new = 1.	×
301	while np.abs(integral_new) > tol:	×
302	r = le	×
303	le = r - step	×
304	subintervals = compute_subintervals(le, r, singularity_list)	×
305	integral_new = gq_modification_composite(integrand, le, r, N,	×
306	subintervals, adaptive,
307	**kwargs)
308	integral += integral_new	×
309
310	elif b == np.inf:	×
311	le = a	×
312	r = a + step	×
313	subintervals = compute_subintervals(le, r, singularity_list)	×
314	integral = gq_modification_composite(integrand, le, r, N, subintervals,	×
315	adaptive, **kwargs)
316	integral_new = 1.	×
317	while np.abs(integral_new) > tol:	×
318	le = r	×
319	r = le + step	×
320	subintervals = compute_subintervals(le, r, singularity_list)	×
321	integral_new = gq_modification_composite(integrand, le, r, N,	×
322	subintervals, adaptive,
323	**kwargs)
324	integral += integral_new	×
325
326	return integral	8✔

SCIInstitute / UncertainSCI / 4659511182

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