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

11.89

/UncertainSCI/model_examples.py

# Contains some rudimentary (physical-space) models for testing PCE
# approximations. All these functions support the syntax output = f(p), where p
# is a d-dimensional vector, and output is a vector whose size is the dimension
# of the model output.

import numpy as np
from scipy import sparse
from scipy.sparse import linalg as splinalg
import scipy.optimize

def ishigami_function(a, b):
    """
    Returns the Ishigami function,

    .. math::

      f(p) = (1 + b*p_3^4) \\sin p_1 + a \\sin^2 p_2,

    where :math:`p = (p_1, p_2, p_3)` are random parameters that are all
    uniformly distributed over :math:`[-\\pi, \\pi]`: :math:`p_i \\sim
    \\mathcal{U}([-\\pi, \\pi])`.
    """

    def f(p):
        return np.sin(p[0]) * (1 + b*p[2]**4) + a*np.sin(p[1])**2

    return f

def borehole_function():
    """
    Returns the Borehole function,

    .. math::

      f(p) = g_1(p) / (g_2(p) * g_3(p)),
      g_1(p) = 2\\pi p_3 * (p_4 - p_6)
      g_2(p) = log(p_2/p_1)
      g_3(p) = 1 + 2*p_7*p_3/(g_2(p) * p_1^2 * p_8) + p_3/p_5

    where the 8-dimensional parameters :math:`p` are

      p = (p_1, p_2, p_3, p_4, p_5, p_6, p_7, p_8)
        = (r_w, r,   T_u, H_u, T_l, H_l, L,   K_w)
    """

    def g1(p):
        return 2*np.pi*p[2]*(p[3] - p[5])

    def g2(p):
        return np.log(p[1]/p[0])

    def g3(p, g2val):
        return 1 + 2*p[6]*p[2]/(g2val * p[0]**2 * p[7]) + p[2]/p[4]

    def f(p):
        g2val = g2(p)
        return g1(p)/(g2val * g3(p,g2val))

    return f

def taylor_frequency(p):
    """
    Returns ( \\sum_{j=1}^d p_j^j )
    """

    return np.sum(p**(1 + np.arange(p.size)))


def sine_modulation(left=-1, right=1, N=100):
    """
    For a d-dimensional parameter p, defines the model,

      f(x,p) = sin [ pi * ( \\sum_{j=1}^d p_j^j ) * x ],

    where x is N equispaced points on the interval [left, right].

    Returns a function pointer with the syntax p ----> f(p).
    """

    x = np.linspace(left, right, N)

    return lambda p: np.sin(np.pi * x * taylor_frequency(p))


def mercer_eigenvalues_exponential_kernel(N, a, b):
    """
    For a 1D exponential covariance kernel,

      K(s,t) = exp(-|t-s| / a),     s, t \\in [-b,b],

    computes the first N eigenvalues of the associated Mercer integral
    operator.

    Precisely, computes the first N/2 positive solutions to both of the following
    transcendental equations for w and v:

       1 - a v tan(v b) = 0
       a w + tan(w b) = 0

    The eigenvalues are subsequently defined through these solutions.

    Returns (1) the N eigenvalues lamb, (2) the first ceil(N/2) solutions for
    v, (3) the first floor(N/2) solutions for w.
    """

    assert N > 0 and a > 0 and b > 0

    M = int(np.ceil(N/2))
    w = np.zeros(M)
    v = np.zeros(M)

    # First equation transformed:
    # vt = v b
    #
    #  -(b/a) / vt + tan(vt) = 0

    def f(x):
        return -(b/a)/x + np.tan(x)

    for n in range(M):
        # Compute bracketing interval
        # root somewhere in right-hand part of [2*n-1, 2*n+1]*pi/2 interval
        RH_value = -1
        k = 4
        while RH_value < 0:
            k += 1
            right = (2*n+1)*np.pi/2 - 1/k
            RH_value = f(right)

        # Root can't be on LHS of interval
        if n == 0:
            left = 1/k
            while f(left) > 0:
                k += 1
                left = 1/k
        else:
            left = n*np.pi

        v[n] = scipy.optimize.brentq(f, left, right)

    v /= b

    # Second equation transformed:
    # wt = w b
    #
    #  (a/b) wt + tan(wt) = 0

    def f(x):
        return (a/b)*x + np.tan(x)

    for n in range(M):

        # Compute bracketing interval
        # root somewhere in [2*n+1, 2*n+3]*pi/2
        LH_value = 1
        k = 4
        while LH_value > 0:
            k += 1
            left = (2*n+1)*np.pi/2 + 1/k
            LH_value = f(left)

        # Root can't be on RHS of interval
        right = (n+1)*np.pi

        w[n] = scipy.optimize.brentq(f, left, right)

    w /= b

    if (N % 2) == 1:  # Don't need last root for w
        w = w[:-1]

    lamb = np.zeros(N)
    oddinds = [i for i in range(N) if (i % 2) == 0]  # Well, odd for 1-based indexing
    lamb[oddinds] = 2*a/(1+(a*v)**2)
    eveninds = [i for i in range(N) if (i % 2) == 1]  # even for 1-based indexing
    lamb[eveninds] = 2*a/(1+(a*w)**2)

    return lamb, v, w


def KLE_exponential_covariance_1d(N, a, b, mn):
    """
    Returns a pointer to a function the evaluates an N-term Karhunen-Loeve
    Expansion of a stochastic process with exponential covariance function on a
    bounded interval [-b,b]. Let the GP have the covariance function,

        C(s,t) = exp(-|t-s| / a),

    and mean function given by mn. Then the N-term KLE of the process is given
    by

        K_N(x,P) = mn(x) + \\sum_{n=1}^N P_n sqrt(\\lambda_n) \\phi_n(x),

    where (lambda_n, phi_n) are the leading eigenpairs of the associated Mercer
    kernel. The eigenvalues are computed in
    mercer_eigenvalues_exponential_kernel. The (P_n) are iid standard normal
    Gaussian random variables.

    Returns a function lamb(x,P) that takes in a 1D np.ndarray x and a 1D
    np.ndarray vector P and returns the KLE realization on x for that value of
    P.
    """

    lamb, v, w = mercer_eigenvalues_exponential_kernel(N, a, b)

    efuns = N*[None]
    for i in range(N):
        if (i % 2) == 0:
            i2 = int(i/2)
            efuns[i] = (lambda i2: lambda x: np.cos(v[i2]*x) / np.sqrt(b + np.sin(2*v[i2]*b)/(2*v[i2])))(i2)
        else:
            i2 = int((i-1)/2)
            efuns[i] = (lambda i2: lambda x: np.sin(w[i2]*x) / np.sqrt(b - np.sin(2*w[i2]*b)/(2*w[i2])))(i2)

    def KLE(x, p):
        return mn(x) + np.array([np.sqrt(lamb[i])*efuns[i](x) for i in range(N)]).T @ p

    return KLE


def laplace_ode_diffusion(x, p):
    """ Parameterized diffusion coefficient for 1D ODE

    For a d-dimensional parameter p, the diffusion coefficient a(x,p) has the form

      a(x,p) = pi^2/5 + sum_{j=1}^d p_j * sin(j*pi*(x+1)/2) / j^2,

    which is positive for all x if all values of p lie between [-1,1].
    """

    a_val = np.ones(x.shape)*np.pi**2/5
    for q in range(p.size):
        a_val += p[q] * np.sin((q+1)*np.pi*(x+1)/2)/(q+1)**2
    return a_val


def laplace_grid_x(left, right, N):
    """
    Computes one-dimensional equispaced grid with N points on the interval
    (left, right).
    """
    return np.linspace(left, right, N)


def laplace_ode(left=-1., right=1., N=100, f=None, diffusion=laplace_ode_diffusion):
    """

    Computes the solution to the ODE:

      -d/dx [ a(x,p) d/dx u(x,p) ] = f(x),

    with homogeneous Dirichlet boundary conditions at x = left, x = right.

    For a d-dimensional parameter p, a(x,p) is the function defined in laplace_ode_diffusion.

    Uses an equispaced finite-difference discretization of the ODE.

    """

    assert N > 2

    if f is None:
        def f(x):
            return np.pi**2 * np.cos(np.pi*x)

    x = laplace_grid_x(left, right, N)
    h = x[1] - x[0]
    fx = f(x)

    # Set homogeneous Dirichlet conditions
    fx[0], fx[-1] = 0., 0.

    # i+1/2 points
    xh = x[:-1] + h/2.

    def create_system(p):
        nonlocal x, xh, N
        a = diffusion(xh, p)
        number_nonzeros = 1 + 1 + (N-2)*3
        rows = np.zeros(number_nonzeros, dtype=int)
        cols = np.zeros(number_nonzeros, dtype=int)
        vals = np.zeros(number_nonzeros, dtype=float)

        # Set the homogeneous Dirichlet conditions
        rows[0], cols[0], vals[0] = 0, 0, 1.
        rows[1], cols[1], vals[1] = N-1, N-1, 1.
        ind = 2

        for q in range(1, N-1):
            # Column q-1
            rows[ind], cols[ind], vals[ind] = q, q-1, -a[q-1]
            ind += 1

            # Column q
            rows[ind], cols[ind], vals[ind] = q, q, a[q-1] + a[q]
            ind += 1

            # Column q+1
            rows[ind], cols[ind], vals[ind] = q, q+1, -a[q]
            ind += 1

        A = sparse.csc_matrix((vals, (rows, cols)), shape=(N, N))

        return A

    def solve_system(p):
        nonlocal fx, h
        return splinalg.spsolve(create_system(p), fx*(h**2))

    return lambda p: solve_system(p)

def laplace_ode_1d(Nparams, a=1., b=1., abar=3., N=100):
    """ Generates a 1D Laplace ODE model with parameterized diffusion.

    Define model:

    -d/dx a(x,p) d/dx u(x,p) = f(x)

    over x in [-1,1], where a(x,p) is a parameterized diffusion model:

    a(x,p) = abar(x) + sum_{j=1}^d lambda_j Y_j phi_j(x),

    where d = Nparams, (lambda_j, phi_j) are eigenpairs of the exponential
    covariance kernel,

      K(s,t) = exp(-|s-t|/a).

    The Y_j are modeled as iid random variables.
    """

    abarfun = lambda x: abar*np.ones(np.shape(x))
    KLE = KLE_exponential_covariance_1d(Nparams, a, b, abarfun)

    diffusion = lambda x, p: KLE(x, p)
    x = laplace_grid_x(-b, b, N)

    return x, laplace_ode(left=-b, right=b, N=N, diffusion=diffusion)

def laplace_grid_xy(left, right, N1, down, up, N2):
    """
    Computes two-dimensional tensorial equispaced grid corresponding to the
    tensorization of N1 equispaced points on the interval (left, right) and N2
    equispaced points on the interval (down, up).
    """
    x = np.linspace(left, right, N1)
    y = np.linspace(down, up, N2)

    X, Y = np.meshgrid(x, y)

    return X.flatten(order='C'), Y.flatten(order='C')


def laplace_pde_diffusion(x, p):
    """ Parameterized diffusion coefficient for 2D PDE

    For a d-dimensional parameter p, the diffusion coefficient a(x,p) has the form

      a(x,p) = pi^2/5 + sum_{j=1}^d p_j * sin(j*pi*(x+1)/2) / j^2,

    which is positive for all x if all values of p lie between [-1,1].
    """

    a_val = np.ones(x.shape)*np.pi**2/5
    for q in range(p.size):
        a_val += p[q] * np.sin((q+1)*np.pi*(x+1)/2)/(q+1)**2
    return a_val


def genz_oscillatory(w=0., c=None):
    """
    Returns a pointer to the "oscillatory" Genz test function defined as

       f(p) = \\cos{ 2\\pi w + \\sum_{i=1}^dim c_i p_i }

    where p \\in R^d. The default value for w is 0, and that for c is a
    d-dimensional vector of ones.
    """

    def cos_eval(p):
        nonlocal c
        if c is None:
            c = np.ones(p.size)
        return np.cos(2*np.pi*w + np.dot(c, p))

    return lambda p: cos_eval(p)


if __name__ == "__main__":

    from matplotlib import pyplot as plt
    import scipy as sp

    dim = 5

    a = 3
    b = 1

    def mn(x):
        return np.zeros(x.shape)

    KLE = KLE_exponential_covariance_1d(dim, a, b, mn)

    def diffusion(x, p):
        return np.exp(KLE(x, p))

    left = -1.
    right = 1.
    N = 1000
    model = laplace_ode(left=left, right=right, N=N, diffusion=diffusion)
    x = laplace_grid_x(left, right, N)

    K = 4
    p = K*[None]
    u = K*[None]
    a = K*[None]

    for k in range(K):
        p[k] = np.random.rand(dim)*2 - 1
        # a[k] = laplace_ode_diffusion(x, p[k])
        a[k] = diffusion(x, p[k])
        u[k] = model(p[k])

    for k in range(K):

        row = np.floor(k/2) + 1
        col = k % (K/2) + 1

        index = col + (row-1)*K/2

        plt.subplot(2, K, k+1)
        plt.plot(x, a[k], 'r')
        plt.title('Diffusion coefficient')
        plt.ylim([0, 3.0])

        plt.subplot(2, K, k+1+K)
        plt.plot(x, u[k])
        plt.title('Solution u')
        plt.ylim([-5, 5])

    M = 1000
    U = np.zeros([u[0].size, M])
    for m in range(M):
        U[m, :] = model(np.random.rand(dim)*2 - 1)

    _, svs, _ = np.linalg.svd(U)
    _, r, _ = sp.linalg.qr(U, pivoting=True)

    plt.figure()
    plt.semilogy(svs[:100], 'r')
    plt.semilogy(np.abs(np.diag(r)[:100]), 'b')
    plt.legend(["Singular values", "Orthogonalization residuals"])
    plt.show()

1	# Contains some rudimentary (physical-space) models for testing PCE
2	# approximations. All these functions support the syntax output = f(p), where p
3	# is a d-dimensional vector, and output is a vector whose size is the dimension
4	# of the model output.
5
6	import numpy as np	8✔
7	from scipy import sparse	8✔
8	from scipy.sparse import linalg as splinalg	8✔
9	import scipy.optimize	8✔
10
11	def ishigami_function(a, b):	8✔
12	"""
13	Returns the Ishigami function,
14
15	.. math::
16
17	f(p) = (1 + b*p_3^4) \\sin p_1 + a \\sin^2 p_2,
18
19	where :math:`p = (p_1, p_2, p_3)` are random parameters that are all
20	uniformly distributed over :math:`[-\\pi, \\pi]`: :math:`p_i \\sim
21	\\mathcal{U}([-\\pi, \\pi])`.
22	"""
23
24	def f(p):	×
25	return np.sin(p[0]) * (1 + bp[2]4) + anp.sin(p[1])**2	×
26
27	return f	×
28
29	def borehole_function():	8✔
30	"""
31	Returns the Borehole function,
32
33	.. math::
34
35	f(p) = g_1(p) / (g_2(p) * g_3(p)),
36	g_1(p) = 2\\pi p_3 * (p_4 - p_6)
37	g_2(p) = log(p_2/p_1)
38	g_3(p) = 1 + 2p_7p_3/(g_2(p) * p_1^2 * p_8) + p_3/p_5
39
40	where the 8-dimensional parameters :math:`p` are
41
42	p = (p_1, p_2, p_3, p_4, p_5, p_6, p_7, p_8)
43	= (r_w, r, T_u, H_u, T_l, H_l, L, K_w)
44	"""
45
46	def g1(p):	×
47	return 2np.pip[2]*(p[3] - p[5])	×
48
49	def g2(p):	×
50	return np.log(p[1]/p[0])	×
51
52	def g3(p, g2val):	×
53	return 1 + 2p[6]p[2]/(g2val * p[0]*2 p[7]) + p[2]/p[4]	×
54
55	def f(p):	×
56	g2val = g2(p)	×
57	return g1(p)/(g2val * g3(p,g2val))	×
58
59	return f	×
60
61	def taylor_frequency(p):	8✔
62	"""
63	Returns ( \\sum_{j=1}^d p_j^j )
64	"""
65
66	return np.sum(p**(1 + np.arange(p.size)))	×
67
68
69	def sine_modulation(left=-1, right=1, N=100):	8✔
70	"""
71	For a d-dimensional parameter p, defines the model,
72
73	f(x,p) = sin [ pi * ( \\sum_{j=1}^d p_j^j ) * x ],
74
75	where x is N equispaced points on the interval [left, right].
76
77	Returns a function pointer with the syntax p ----> f(p).
78	"""
79
80	x = np.linspace(left, right, N)	×
81
82	return lambda p: np.sin(np.pi * x * taylor_frequency(p))	×
83
84
85	def mercer_eigenvalues_exponential_kernel(N, a, b):	8✔
86	"""
87	For a 1D exponential covariance kernel,
88
89	K(s,t) = exp(-\|t-s\| / a), s, t \\in [-b,b],
90
91	computes the first N eigenvalues of the associated Mercer integral
92	operator.
93
94	Precisely, computes the first N/2 positive solutions to both of the following
95	transcendental equations for w and v:
96
97	1 - a v tan(v b) = 0
98	a w + tan(w b) = 0
99
100	The eigenvalues are subsequently defined through these solutions.
101
102	Returns (1) the N eigenvalues lamb, (2) the first ceil(N/2) solutions for
103	v, (3) the first floor(N/2) solutions for w.
104	"""
105
106	assert N > 0 and a > 0 and b > 0	×
107
108	M = int(np.ceil(N/2))	×
109	w = np.zeros(M)	×
110	v = np.zeros(M)	×
111
112	# First equation transformed:
113	# vt = v b
114	#
115	# -(b/a) / vt + tan(vt) = 0
116
117	def f(x):	×
118	return -(b/a)/x + np.tan(x)	×
119
120	for n in range(M):	×
121	# Compute bracketing interval
122	# root somewhere in right-hand part of [2n-1, 2n+1]*pi/2 interval
123	RH_value = -1	×
124	k = 4	×
125	while RH_value < 0:	×
126	k += 1	×
127	right = (2n+1)np.pi/2 - 1/k	×
128	RH_value = f(right)	×
129
130	# Root can't be on LHS of interval
131	if n == 0:	×
132	left = 1/k	×
133	while f(left) > 0:	×
134	k += 1	×
135	left = 1/k	×
136	else:
137	left = n*np.pi	×
138
139	v[n] = scipy.optimize.brentq(f, left, right)	×
140
141	v /= b	×
142
143	# Second equation transformed:
144	# wt = w b
145	#
146	# (a/b) wt + tan(wt) = 0
147
148	def f(x):	×
149	return (a/b)*x + np.tan(x)	×
150
151	for n in range(M):	×
152
153	# Compute bracketing interval
154	# root somewhere in [2n+1, 2n+3]*pi/2
155	LH_value = 1	×
156	k = 4	×
157	while LH_value > 0:	×
158	k += 1	×
159	left = (2n+1)np.pi/2 + 1/k	×
160	LH_value = f(left)	×
161
162	# Root can't be on RHS of interval
163	right = (n+1)*np.pi	×
164
165	w[n] = scipy.optimize.brentq(f, left, right)	×
166
167	w /= b	×
168
169	if (N % 2) == 1: # Don't need last root for w	×
170	w = w[:-1]	×
171
172	lamb = np.zeros(N)	×
173	oddinds = [i for i in range(N) if (i % 2) == 0] # Well, odd for 1-based indexing	×
174	lamb[oddinds] = 2a/(1+(av)**2)	×
175	eveninds = [i for i in range(N) if (i % 2) == 1] # even for 1-based indexing	×
176	lamb[eveninds] = 2a/(1+(aw)**2)	×
177
178	return lamb, v, w	×
179
180
181	def KLE_exponential_covariance_1d(N, a, b, mn):	8✔
182	"""
183	Returns a pointer to a function the evaluates an N-term Karhunen-Loeve
184	Expansion of a stochastic process with exponential covariance function on a
185	bounded interval [-b,b]. Let the GP have the covariance function,
186
187	C(s,t) = exp(-\|t-s\| / a),
188
189	and mean function given by mn. Then the N-term KLE of the process is given
190	by
191
192	K_N(x,P) = mn(x) + \\sum_{n=1}^N P_n sqrt(\\lambda_n) \\phi_n(x),
193
194	where (lambda_n, phi_n) are the leading eigenpairs of the associated Mercer
195	kernel. The eigenvalues are computed in
196	mercer_eigenvalues_exponential_kernel. The (P_n) are iid standard normal
197	Gaussian random variables.
198
199	Returns a function lamb(x,P) that takes in a 1D np.ndarray x and a 1D
200	np.ndarray vector P and returns the KLE realization on x for that value of
201	P.
202	"""
203
204	lamb, v, w = mercer_eigenvalues_exponential_kernel(N, a, b)	×
205
206	efuns = N*[None]	×
207	for i in range(N):	×
208	if (i % 2) == 0:	×
209	i2 = int(i/2)	×
210	efuns[i] = (lambda i2: lambda x: np.cos(v[i2]x) / np.sqrt(b + np.sin(2v[i2]b)/(2v[i2])))(i2)	×
211	else:
212	i2 = int((i-1)/2)	×
213	efuns[i] = (lambda i2: lambda x: np.sin(w[i2]x) / np.sqrt(b - np.sin(2w[i2]b)/(2w[i2])))(i2)	×
214
215	def KLE(x, p):	×
216	return mn(x) + np.array([np.sqrt(lamb[i])*efuns[i](x) for i in range(N)]).T @ p	×
217
218	return KLE	×
219
220
221	def laplace_ode_diffusion(x, p):	8✔
222	""" Parameterized diffusion coefficient for 1D ODE
223
224	For a d-dimensional parameter p, the diffusion coefficient a(x,p) has the form
225
226	a(x,p) = pi^2/5 + sum_{j=1}^d p_j * sin(jpi(x+1)/2) / j^2,
227
228	which is positive for all x if all values of p lie between [-1,1].
229	"""
230
231	a_val = np.ones(x.shape)np.pi*2/5	×
232	for q in range(p.size):	×
233	a_val += p[q] * np.sin((q+1)np.pi(x+1)/2)/(q+1)**2	×
234	return a_val	×
235
236
237	def laplace_grid_x(left, right, N):	8✔
238	"""
239	Computes one-dimensional equispaced grid with N points on the interval
240	(left, right).
241	"""
242	return np.linspace(left, right, N)	×
243
244
245	def laplace_ode(left=-1., right=1., N=100, f=None, diffusion=laplace_ode_diffusion):	8✔
246	"""
247
248	Computes the solution to the ODE:
249
250	-d/dx [ a(x,p) d/dx u(x,p) ] = f(x),
251
252	with homogeneous Dirichlet boundary conditions at x = left, x = right.
253
254	For a d-dimensional parameter p, a(x,p) is the function defined in laplace_ode_diffusion.
255
256	Uses an equispaced finite-difference discretization of the ODE.
257
258	"""
259
260	assert N > 2	×
261
262	if f is None:	×
263	def f(x):	×
264	return np.pi*2 np.cos(np.pi*x)	×
265
266	x = laplace_grid_x(left, right, N)	×
267	h = x[1] - x[0]	×
268	fx = f(x)	×
269
270	# Set homogeneous Dirichlet conditions
271	fx[0], fx[-1] = 0., 0.	×
272
273	# i+1/2 points
274	xh = x[:-1] + h/2.	×
275
276	def create_system(p):	×
277	nonlocal x, xh, N
278	a = diffusion(xh, p)	×
279	number_nonzeros = 1 + 1 + (N-2)*3	×
280	rows = np.zeros(number_nonzeros, dtype=int)	×
281	cols = np.zeros(number_nonzeros, dtype=int)	×
282	vals = np.zeros(number_nonzeros, dtype=float)	×
283
284	# Set the homogeneous Dirichlet conditions
285	rows[0], cols[0], vals[0] = 0, 0, 1.	×
286	rows[1], cols[1], vals[1] = N-1, N-1, 1.	×
287	ind = 2	×
288
289	for q in range(1, N-1):	×
290	# Column q-1
291	rows[ind], cols[ind], vals[ind] = q, q-1, -a[q-1]	×
292	ind += 1	×
293
294	# Column q
295	rows[ind], cols[ind], vals[ind] = q, q, a[q-1] + a[q]	×
296	ind += 1	×
297
298	# Column q+1
299	rows[ind], cols[ind], vals[ind] = q, q+1, -a[q]	×
300	ind += 1	×
301
302	A = sparse.csc_matrix((vals, (rows, cols)), shape=(N, N))	×
303
304	return A	×
305
306	def solve_system(p):	×
307	nonlocal fx, h
308	return splinalg.spsolve(create_system(p), fx(h*2))	×
309
310	return lambda p: solve_system(p)	×
311
312	def laplace_ode_1d(Nparams, a=1., b=1., abar=3., N=100):	8✔
313	""" Generates a 1D Laplace ODE model with parameterized diffusion.
314
315	Define model:
316
317	-d/dx a(x,p) d/dx u(x,p) = f(x)
318
319	over x in [-1,1], where a(x,p) is a parameterized diffusion model:
320
321	a(x,p) = abar(x) + sum_{j=1}^d lambda_j Y_j phi_j(x),
322
323	where d = Nparams, (lambda_j, phi_j) are eigenpairs of the exponential
324	covariance kernel,
325
326	K(s,t) = exp(-\|s-t\|/a).
327
328	The Y_j are modeled as iid random variables.
329	"""
330
331	abarfun = lambda x: abar*np.ones(np.shape(x))	×
332	KLE = KLE_exponential_covariance_1d(Nparams, a, b, abarfun)	×
333
334	diffusion = lambda x, p: KLE(x, p)	×
335	x = laplace_grid_x(-b, b, N)	×
336
337	return x, laplace_ode(left=-b, right=b, N=N, diffusion=diffusion)	×
338
339	def laplace_grid_xy(left, right, N1, down, up, N2):	8✔
340	"""
341	Computes two-dimensional tensorial equispaced grid corresponding to the
342	tensorization of N1 equispaced points on the interval (left, right) and N2
343	equispaced points on the interval (down, up).
344	"""
345	x = np.linspace(left, right, N1)	×
346	y = np.linspace(down, up, N2)	×
347
348	X, Y = np.meshgrid(x, y)	×
349
350	return X.flatten(order='C'), Y.flatten(order='C')	×
351
352
353	def laplace_pde_diffusion(x, p):	8✔
354	""" Parameterized diffusion coefficient for 2D PDE
355
356	For a d-dimensional parameter p, the diffusion coefficient a(x,p) has the form
357
358	a(x,p) = pi^2/5 + sum_{j=1}^d p_j * sin(jpi(x+1)/2) / j^2,
359
360	which is positive for all x if all values of p lie between [-1,1].
361	"""
362
363	a_val = np.ones(x.shape)np.pi*2/5	×
364	for q in range(p.size):	×
365	a_val += p[q] * np.sin((q+1)np.pi(x+1)/2)/(q+1)**2	×
366	return a_val	×
367
368
369	def genz_oscillatory(w=0., c=None):	8✔
370	"""
371	Returns a pointer to the "oscillatory" Genz test function defined as
372
373	f(p) = \\cos{ 2\\pi w + \\sum_{i=1}^dim c_i p_i }
374
375	where p \\in R^d. The default value for w is 0, and that for c is a
376	d-dimensional vector of ones.
377	"""
378
379	def cos_eval(p):	8✔
380	nonlocal c
381	if c is None:	8✔
382	c = np.ones(p.size)	×
383	return np.cos(2np.piw + np.dot(c, p))	8✔
384
385	return lambda p: cos_eval(p)	8✔
386
387
388	if __name__ == "__main__":	8✔
389
390	from matplotlib import pyplot as plt	×
391	import scipy as sp	×
392
393	dim = 5	×
394
395	a = 3	×
396	b = 1	×
397
398	def mn(x):	×
399	return np.zeros(x.shape)	×
400
401	KLE = KLE_exponential_covariance_1d(dim, a, b, mn)	×
402
403	def diffusion(x, p):	×
404	return np.exp(KLE(x, p))	×
405
406	left = -1.	×
407	right = 1.	×
408	N = 1000	×
409	model = laplace_ode(left=left, right=right, N=N, diffusion=diffusion)	×
410	x = laplace_grid_x(left, right, N)	×
411
412	K = 4	×
413	p = K*[None]	×
414	u = K*[None]	×
415	a = K*[None]	×
416
417	for k in range(K):	×
418	p[k] = np.random.rand(dim)*2 - 1	×
419	# a[k] = laplace_ode_diffusion(x, p[k])
420	a[k] = diffusion(x, p[k])	×
421	u[k] = model(p[k])	×
422
423	for k in range(K):	×
424
425	row = np.floor(k/2) + 1	×
426	col = k % (K/2) + 1	×
427
428	index = col + (row-1)*K/2	×
429
430	plt.subplot(2, K, k+1)	×
431	plt.plot(x, a[k], 'r')	×
432	plt.title('Diffusion coefficient')	×
433	plt.ylim([0, 3.0])	×
434
435	plt.subplot(2, K, k+1+K)	×
436	plt.plot(x, u[k])	×
437	plt.title('Solution u')	×
438	plt.ylim([-5, 5])	×
439
440	M = 1000	×
441	U = np.zeros([u[0].size, M])	×
442	for m in range(M):	×
443	U[m, :] = model(np.random.rand(dim)*2 - 1)	×
444
445	_, svs, _ = np.linalg.svd(U)	×
446	_, r, _ = sp.linalg.qr(U, pivoting=True)	×
447
448	plt.figure()	×
449	plt.semilogy(svs[:100], 'r')	×
450	plt.semilogy(np.abs(np.diag(r)[:100]), 'b')	×
451	plt.legend(["Singular values", "Orthogonalization residuals"])	×
452	plt.show()	×

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