From 3c31261d145cdab8533e0729e1f621eaa3590f23 Mon Sep 17 00:00:00 2001 From: orbitfold Date: Tue, 13 Aug 2019 13:00:50 +0200 Subject: [PATCH 01/48] Initial campaigndb test --- tests/test_db.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) create mode 100644 tests/test_db.py diff --git a/tests/test_db.py b/tests/test_db.py new file mode 100644 index 000000000..2679203cb --- /dev/null +++ b/tests/test_db.py @@ -0,0 +1,15 @@ +import pytest +import os.path +import easyvvuq +from easyvvuq.constants import default_campaign_prefix +from easyvvuq.db.sql import CampaignDB +from easyvvuq.data_structs import CampaignInfo + +def test_file_created(tmp_path): + info = CampaignInfo( + name='test', + campaign_dir_prefix=default_campaign_prefix, + easyvvuq_version=easyvvuq.__version__, + campaign_dir='.') + campaign = CampaignDB(location='sqlite:///{}/test.sqlite'.format(tmp_path), new_campaign=True, name='test', info=info) + assert(os.path.isfile('{}/test.sqlite'.format(tmp_path))) From 53fd2e886d95c7680887702d174981ea333d5ddf Mon Sep 17 00:00:00 2001 From: orbitfold Date: Wed, 14 Aug 2019 15:35:15 +0200 Subject: [PATCH 02/48] Added a campaign fixture: --- tests/test_db.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/tests/test_db.py b/tests/test_db.py index 2679203cb..415a7a3b3 100644 --- a/tests/test_db.py +++ b/tests/test_db.py @@ -5,11 +5,19 @@ from easyvvuq.db.sql import CampaignDB from easyvvuq.data_structs import CampaignInfo -def test_file_created(tmp_path): +@pytest.fixture +def campaign(tmp_path): info = CampaignInfo( name='test', campaign_dir_prefix=default_campaign_prefix, easyvvuq_version=easyvvuq.__version__, campaign_dir='.') campaign = CampaignDB(location='sqlite:///{}/test.sqlite'.format(tmp_path), new_campaign=True, name='test', info=info) - assert(os.path.isfile('{}/test.sqlite'.format(tmp_path))) + campaign.tmp_path = tmp_path + return campaign + +def test_db_file_created(campaign): + assert(os.path.isfile('{}/test.sqlite'.format(campaign.tmp_path))) + +def test_add_runs(tmp_path): + pass From 51bd623e51c34a1318dfbb00a70747224e144f70 Mon Sep 17 00:00:00 2001 From: wedeling Date: Wed, 14 Aug 2019 21:24:13 +0200 Subject: [PATCH 03/48] adding sparse grid capability --- easyvvuq/analysis/sc_analysis.py | 227 ++++++++++++++++++-- easyvvuq/sampling/stochastic_collocation.py | 122 +++++++++-- 2 files changed, 316 insertions(+), 33 deletions(-) diff --git a/easyvvuq/analysis/sc_analysis.py b/easyvvuq/analysis/sc_analysis.py index 9749fb9a6..9c668a005 100644 --- a/easyvvuq/analysis/sc_analysis.py +++ b/easyvvuq/analysis/sc_analysis.py @@ -94,11 +94,23 @@ def analyse(self, data_frame=None): raise RuntimeError( "No data in data frame passed to analyse element") - qoi_cols = self.qoi_cols - results = {'statistical_moments': {}, - 'sobol_indices': {k: {} for k in qoi_cols}} + 'sobol_indices': {k: {} for k in self.qoi_cols}} + + if self.sampler.sparse == False: + self.analyse_tensor_grid(results, self.qoi_cols, data_frame) + else: + self.analyse_sparse_grid(results, self.qoi_cols, data_frame) + return results + + + def analyse_tensor_grid(self, results, qoi_cols, data_frame): + """ + Perform the post-processing analysis for a full tensor product grid + (sparse = False) + """ + # Chaospy computation of 1D weights xi = [] wi = [] @@ -145,6 +157,104 @@ def analyse(self, data_frame=None): results['sobol_indices'][qoi_k] = self.get_sobol_indices(qoi_k, 'all') return results + + def analyse_sparse_grid(self, results, qoi_cols, data_frame): + """ + Perform the post-processing analysis for a sparse grid + (sparse = True) + """ + #use this above too? + self.xi_d = self.sampler.xi_d + self.xi_1d = self.sampler.xi_1d + self.wi_1d = self.sampler.wi_1d + + self.L = self.sampler.q + self.N = self.sampler.N + + self.Map = {} + for level in range(self.N, self.L+1): + self.Map[level] = self.create_map(self.N, level) + + # Extract output values for each quantity of interest from Dataframe + idx = 0 + samples = {k: {} for k in qoi_cols} + + for run_id in data_frame.run_id.unique(): + for k in qoi_cols: + values = data_frame.loc[data_frame['run_id'] == run_id][k] + samples[k][self.xi_d[idx].tobytes()] = values + idx += 1 + self.samples = samples + + self.N_qoi = samples[qoi_cols[0]][self.xi_d[0].tobytes()].size + + # Extract output values for each quantity of interest from Dataframe + samples2 = {k: [] for k in qoi_cols} + for run_id in data_frame.run_id.unique(): + for k in qoi_cols: + values = data_frame.loc[data_frame['run_id'] == run_id][k].values + samples2[k].append(values) + self.samples2 = samples2 + self.data_frame = data_frame + + self.test = [] + + # Compute descriptive statistics for each quantity of interest + for qoi_k in qoi_cols: + mean_k, var_k = self.get_moments(qoi_k) + std_k = var_k**0.5 + + # compute statistical moments + results['statistical_moments'][qoi_k] = {'mean': mean_k, + 'var': var_k, + 'std': std_k} + # compute all Sobol indices + #results['sobol_indices'][qoi_k] = self.get_sobol_indices(qoi_k, 'all') + + def compute_Q_diff(self, i, qoi, mean_f = []): + + """ + ======================================================================= + For every multi index i = (i1, i2, ..., id), Smolyak sums over + tensor products difference quadrature rules: + (Q^1_{i1} - Q^1_{i1-1}) X ... X (Q^1_{id) - Q^1_{id-1}) + Below this product is expanded into individual tensor products, each + of which is then computed as: + Q^1_{k1} X ... X Q^1_{kd} = sum...sum w_{k1}*...*w{kd}*f(x_{k1},...,x_{kd}) + ======================================================================= + """ + #expand the multi-index indices of the tensor product + #(Q^1_{i1} - Q^1_{i1-1}) X ... X (Q^1_{id) - Q^1_{id-1}) + diff_idx = np.array(list(product(*[[k, -(k-1)] for k in i]))) + + #Delta will be the sum of all expanded tensor products + #Q^1_{k1} X ... X Q^1_{kd} = sum...sum w_{k1}*...*w{kd}*f(x_{k1},...,x_{kd}) + Delta = 0.0 + + for diff in diff_idx: + #Q_0 = 0, so if present do not compute Q^1_{k1} X ... X Q^1_{kd} + if not np.in1d(0, diff): + X_k = [self.xi_1d[n][np.abs(diff)[n]] for n in range(self.N)] + X_k = np.array(list(product(*X_k))) + f_k = np.array([self.samples[qoi][x.tobytes()] for x in X_k]) + W_k = [self.wi_1d[n][np.abs(diff)[n]] for n in range(self.N)] + W_k = np.array(list(product(*W_k))) + W_k = np.prod(W_k, axis=1) + + #do some reshaping to align for multiplication + f_k = f_k.reshape([f_k.shape[0], self.N_qoi]) + W_k = W_k.reshape([W_k.shape[0], 1]) + + if mean_f == []: + Q_prod = np.sum(f_k*W_k, axis=0).T + else: + Q_prod = np.sum((f_k - mean_f)**2, axis=0).T + else: + Q_prod = 0.0 + + Delta += np.sign(np.prod(diff))*Q_prod + + return Delta def get_moments(self, qoi): """Compute the first two statistical moments. @@ -161,20 +271,103 @@ def get_moments(self, qoi): float: Variance of samples, using quad weights """ - - # compute the mean and variance of the code samples, using quad weights - mean_f = np.zeros([self.N_qoi, 1]) - var_f = np.zeros([self.N_qoi, 1]) - - for k in range(self._number_of_samples): - sample_k = (self.samples[qoi][k]).values.reshape([self.N_qoi, 1]) - mean_f += sample_k * self.wi_d[k].prod() - + + if self.sampler.sparse == False: + mean_f = np.zeros([self.N_qoi, 1]) + var_f = np.zeros([self.N_qoi, 1]) + + for k in range(self._number_of_samples): + sample_k = (self.samples[qoi][k]).values.reshape([self.N_qoi, 1]) + mean_f += sample_k * self.wi_d[k]#.prod() + + for k in range(self._number_of_samples): + sample_k = (self.samples[qoi][k]).values.reshape([self.N_qoi, 1]) + var_f += (sample_k - mean_f)**2 * self.wi_d[k]#.prod() + + return mean_f, var_f + else: + + i_norm_le_q = self.sampler.compute_sparse_multi_idx(self.L, self.N) + + #compute quadrature sum_i(Delta_i1 X ... X Delta__iN) + Delta = [] + for i in i_norm_le_q: + Delta.append(self.compute_Q_diff(i, qoi)) + + mean_f = np.sum(Delta, axis = 0) + + #compute quadrature sum_i(Delta_i1 X ... X Delta__iN) + Delta = [] + for i in i_norm_le_q: + Delta.append(self.compute_Q_diff(i, qoi, mean_f)) + + var_f = np.sum(Delta, axis = 0) + + return mean_f, var_f + + def get_sample_array(self, qoi): + """ + Returns an array of all samples of qoi + """ + + tmp = np.zeros([self._number_of_samples, self.N_qoi]) + for k in range(self._number_of_samples): - sample_k = (self.samples[qoi][k]).values.reshape([self.N_qoi, 1]) - var_f += (sample_k - mean_f)**2 * self.wi_d[k].prod() - - return mean_f, var_f + tmp[k, :] = (self.samples[qoi][k]).values + + return tmp + + def create_map(self, N, L): + l_norm_le_L = self.sampler.compute_sparse_multi_idx(L, N) + + k = 0 + Map = {} + + for l in l_norm_le_L: + X_l = [self.xi_1d[n][l[n]] for n in range(N)] + X_l = np.array(list(product(*X_l))) + + for x in X_l: + j = np.where((x == self.xi_d).all(axis=1))[0][0] + Map[k] = {'l':l, 'X':x, 'f':j} + k += 1; + + return Map + + def sparse_surrogate(self, qoi, x, N, L): + + if L < N: + return 0.0 + + surr = np.zeros(self.N_qoi) + + for level in range(N, L+1): + + Delta = np.zeros(self.N_qoi) + Map = self.Map[level] + + for k in Map.keys(): + + q_k = self.samples2[qoi][Map[k]['f']] + s_k = q_k - self.sparse_surrogate(qoi, Map[k]['X'], N, level-1) + + l = Map[k]['l'] + + idx = {} + for n in range(N): + # indices of current collocation point xi_d[k] in 1d + # collocation points + idx[n] = (self.xi_1d[n][l[n]] == Map[k]['X'][n]).nonzero()[0][0] + + lagrange = [] + for n in range(N): + # values of Lagrange polynomials at x + lagrange.append(lagrange_poly(x[n], self.xi_1d[n][l[n]], idx[n])) + + Delta += s_k*np.prod(lagrange) + surr += Delta + + return surr def surrogate(self, qoi, x): """Use the SC expansion as a surrogate. @@ -427,7 +620,7 @@ def lagrange_poly(x, x_i, j): x_i : list or array of float - j : float + j : int Returns ------- diff --git a/easyvvuq/sampling/stochastic_collocation.py b/easyvvuq/sampling/stochastic_collocation.py index 0d9a6dd5c..d02684f58 100644 --- a/easyvvuq/sampling/stochastic_collocation.py +++ b/easyvvuq/sampling/stochastic_collocation.py @@ -1,5 +1,7 @@ from .base import BaseSamplingElement, Vary import chaospy as cp +import numpy as np +from itertools import product, chain import logging __author__ = "Wouter Edeling" @@ -30,11 +32,13 @@ class SCSampler(BaseSamplingElement, sampler_name="sc_sampler"): def __init__(self, vary=None, - polynomial_order=4, + quad_order=4, quadrature_rule="G", - count=0): + count=0, + growth=False, + sparse=False): """ - Create the sampler for the Polynomial Chaos Expansion method. + Create the sampler for the Stochastic Collocation method. Parameters ---------- @@ -45,6 +49,11 @@ def __init__(self, quadrature_rule : char, optional The quadrature method, default is Gaussian "G". + + growth: bool, optional + Sets the growth rule to exponential for Clenshaw Curtis quadrature, + which makes it nested, and therefore more efficient for sparse grids. + Default is False. sparse : bool, optional If True use sparse grid instead of normal tensor product grid, @@ -52,7 +61,7 @@ def __init__(self, """ self.vary = Vary(vary) - self.polynomial_order = polynomial_order + #self.polynomial_order = polynomial_order self.quadrature_rule = quadrature_rule # List of the probability distributions of uncertain parameters @@ -64,16 +73,59 @@ def __init__(self, self.joint_dist = cp.J(*params_distribution) # The quadrature information: order, rule and sparsity - self.quad_order = polynomial_order + #self.quad_order = polynomial_order + self.quad_order = quad_order self.quad_rule = quadrature_rule - # self.quad_sparse = sparse - - # the nodes of the collocation grid - xi_d, _ = cp.generate_quadrature(self.quad_order, - self.joint_dist, - rule=quadrature_rule) - - self.xi_d = xi_d.T + self.sparse = sparse + self.quad_sparse = sparse + self.growth = growth + + if sparse == False: + # the nodes of the collocation grid + xi_d, _ = cp.generate_quadrature(self.quad_order, + self.joint_dist, + rule=quadrature_rule) + self.xi_d = xi_d.T + #sparse grid = a linear combination of tensor products of 1D rules + #of different order. Use chaospy to compute these 1D quadrature rules + else: + #q = level of sparse grid + q = self.quad_order + #N = number of uncertain parameters + N = len(params_distribution) + + #q >= N must hold + if q < N: + print("*************************************************************") + print("Level of sparse grid is lower than the dimension N (# params)") + print("Increase level (quad_order) q such that q >= N") + print("*************************************************************") + import sys; sys.exit() + + #for every dimension (parameter), create a hierachy of 1D + #quadrature rules or increasing order + + self.xi_1d = {} + self.wi_1d = {} + + for n in range(N): + self.xi_1d[n] = {} + self.wi_1d[n] = {} + + for n in range(N): + for i in range(1, quad_order): + xi_i, wi_i = cp.generate_quadrature(i, + params_distribution[n], + rule=self.quad_rule, + growth=self.growth) + self.xi_1d[n][i] = xi_i[0] + self.wi_1d[n][i] = wi_i + + #create sparse grid of dimension N and level q using the 1d + #rules in self.xi_1d + self.xi_d = self.sparse_grid_q_N(q, N) + self.q = q + self.N = N self._number_of_samples = self.xi_d.shape[0] @@ -89,7 +141,7 @@ def __init__(self, self.__next__() def element_version(self): - return "0.3" + return "0.4" def is_finite(self): return True @@ -114,6 +166,44 @@ def is_restartable(self): def get_restart_dict(self): return { "vary": self.vary.serialize(), - "polynomial_order": self.polynomial_order, + "quad_order": self.quad_order, "quadrature_rule": self.quadrature_rule, - "count": self.count} + "count": self.count, + "growth": self.growth, + "sparse": self.sparse} + + """ + ======================= + SPARSE GRID SUBROUTINES + ======================= + """ + + def sparse_grid_q_N(self, q, N): + """ + Compute an isotropic sparse grid H_q_N of N dimensions and level q + """ + #multi-index i, such that |i| <= q + i_norm_le_q = self.compute_sparse_multi_idx(q, N) + + H_q_N = [] + #loop over all multi indices i + for i in i_norm_le_q: + + #compute the tensor product of nodes indexed by i + X_i = [self.xi_1d[n][i[n]] for n in range(N)] + H_q_N.append(list(product(*X_i))) + + #flatten the list of lists + H_q_N = np.array(list(chain(*H_q_N))) + + #return unique nodes + return np.unique(H_q_N, axis=0) + + def compute_sparse_multi_idx(self, q, N): + """ + computes all N dimensional multi-indices i = (i1,...,iN) such that + |i| <= Q. Here |i| is the internal sum of i (i1+...+iN) + """ + P = np.array(list(product(range(1,q+1), repeat=N))) + i_norm_le_q = P[np.where(np.sum(P, axis=1) <= q)[0]] + return i_norm_le_q \ No newline at end of file From 8328434b2006440c832e3e2878f9b7afd84d790e Mon Sep 17 00:00:00 2001 From: wedeling Date: Thu, 15 Aug 2019 15:56:29 +0200 Subject: [PATCH 04/48] added SC UQP --- easyvvuq/analysis/sc_analysis.py | 500 ++++++++++---------- easyvvuq/sampling/stochastic_collocation.py | 75 +-- 2 files changed, 292 insertions(+), 283 deletions(-) diff --git a/easyvvuq/analysis/sc_analysis.py b/easyvvuq/analysis/sc_analysis.py index 9c668a005..e310084dd 100644 --- a/easyvvuq/analysis/sc_analysis.py +++ b/easyvvuq/analysis/sc_analysis.py @@ -62,6 +62,19 @@ def __init__(self, sampler=None, qoi_cols=None): self.output_type = OutputType.SUMMARY self.sampler = sampler self._number_of_samples = sampler._number_of_samples + self.sparse = sampler.sparse + + #determine if quadrature rule is nested, default assumption is False + self.nested = False + rule = self.sampler.quadrature_rule + growth = self.sampler.growth + + #TODO: ARE THERE OTHERS WHICH ARE NESTED? genz_keister SHOULD BE, + #BUT SEEMS NOT IMPLEMENTED IN MY cp VERSION (3.0.5) + if rule == 'leja' or rule == 'guass_patterson': + self.nested = True + elif rule == 'c' or rule == 'clenshaw_curtis' and growth == True: + self.nested = True def element_name(self): """Name for this element for logging purposes""" @@ -94,24 +107,8 @@ def analyse(self, data_frame=None): raise RuntimeError( "No data in data frame passed to analyse element") - results = {'statistical_moments': {}, - 'sobol_indices': {k: {} for k in self.qoi_cols}} - - if self.sampler.sparse == False: - self.analyse_tensor_grid(results, self.qoi_cols, data_frame) - else: - self.analyse_sparse_grid(results, self.qoi_cols, data_frame) - - return results - - - def analyse_tensor_grid(self, results, qoi_cols, data_frame): - """ - Perform the post-processing analysis for a full tensor product grid - (sparse = False) - """ - - # Chaospy computation of 1D weights + # Chaospy computation of 1D weights - USED IN SOBOL SUBROUTINES + #REMOVE LATER IN FAVOUR OF xi_1D, wi_1d xi = [] wi = [] for dist in self.sampler.vary.get_values(): @@ -123,27 +120,52 @@ def analyse_tensor_grid(self, results, qoi_cols, data_frame): self.wi = wi # Compute tensor product nodes and weights - xi_d, self.wi_d = cp.generate_quadrature(order=self.sampler.quad_order, + # REMOVE WHEN MOMENT SUBROUTINE USES SC WEIGHTS + _, self.wi_d = cp.generate_quadrature(order=self.sampler.quad_order, domain=self.sampler.joint_dist, rule=self.sampler.quad_rule) - self.xi_d = xi_d.T + + self.xi_d = self.sampler.xi_d + self.xi_1d = self.sampler.xi_1d + self.wi_1d = self.sampler.wi_1d + + #the maximum level (quad order) of the (sparse) grid + self.L = self.sampler.L + + #the number of uncertain parameters + self.N = self.sampler.N + + #if L < L_min: quadratures and interpolations are zero + #For full tensor grid: there is only one level: L_min = L + if self.sparse == False: + self.L_min = self.L + #For sparse grid: multiple levels, L >= N must hold + else: + self.L_min = self.N + + #per level, map a unique index k to all (level multi indices, colloc points) + #combinations. Will differ for sparse or full tensor grids. + #All interpolation/quadrature subroutines loop over the entries in Map + self.Map = {} + for level in range(self.L_min, self.L+1): + self.Map[level] = self.create_map(self.N, level) # Extract output values for each quantity of interest from Dataframe + qoi_cols = self.qoi_cols samples = {k: [] for k in qoi_cols} for run_id in data_frame.run_id.unique(): for k in qoi_cols: - values = data_frame.loc[data_frame['run_id'] == run_id][k] + values = data_frame.loc[data_frame['run_id'] == run_id][k].values samples[k].append(values) - self.samples = samples - # number of uncertain parameters - self.d = self.xi_d.shape[1] - # size of one code sample self.N_qoi = self.samples[qoi_cols[0]][0].size - # Compute descriptive statistics for each quantity of interest + results = {'statistical_moments': {}, + 'sobol_indices': {k: {} for k in self.qoi_cols}} + + # Compute descriptive statistics for each quantity of interest for qoi_k in qoi_cols: mean_k, var_k = self.get_moments(qoi_k) std_k = var_k**0.5 @@ -152,266 +174,254 @@ def analyse_tensor_grid(self, results, qoi_cols, data_frame): results['statistical_moments'][qoi_k] = {'mean': mean_k, 'var': var_k, 'std': std_k} - - # compute all Sobol indices - results['sobol_indices'][qoi_k] = self.get_sobol_indices(qoi_k, 'all') +# # compute all Sobol indices +# results['sobol_indices'][qoi_k] = self.get_sobol_indices(qoi_k, 'all') return results - def analyse_sparse_grid(self, results, qoi_cols, data_frame): - """ - Perform the post-processing analysis for a sparse grid - (sparse = True) + def create_map(self, N, L): """ - #use this above too? - self.xi_d = self.sampler.xi_d - self.xi_1d = self.sampler.xi_1d - self.wi_1d = self.sampler.wi_1d - - self.L = self.sampler.q - self.N = self.sampler.N + Create a map from a unique integer k to each + (level multi index l, collocation point X) combination. Also + compute the index of X (f) in the global (sparse) grid xi_d - self.Map = {} - for level in range(self.N, self.L+1): - self.Map[level] = self.create_map(self.N, level) - - # Extract output values for each quantity of interest from Dataframe - idx = 0 - samples = {k: {} for k in qoi_cols} + Parameters + ---------- + - N (int) = number of parameters + - L (int) = max level of grid - for run_id in data_frame.run_id.unique(): - for k in qoi_cols: - values = data_frame.loc[data_frame['run_id'] == run_id][k] - samples[k][self.xi_d[idx].tobytes()] = values - idx += 1 - self.samples = samples - - self.N_qoi = samples[qoi_cols[0]][self.xi_d[0].tobytes()].size + Returns + -------- + - Map: a dict for level L containing k, l, X, and f + """ - # Extract output values for each quantity of interest from Dataframe - samples2 = {k: [] for k in qoi_cols} - for run_id in data_frame.run_id.unique(): - for k in qoi_cols: - values = data_frame.loc[data_frame['run_id'] == run_id][k].values - samples2[k].append(values) - self.samples2 = samples2 - self.data_frame = data_frame - - self.test = [] + #unique index + k = 0 + Map = {} - # Compute descriptive statistics for each quantity of interest - for qoi_k in qoi_cols: - mean_k, var_k = self.get_moments(qoi_k) - std_k = var_k**0.5 + #full tensor product + if self.sparse == False: + + l = (np.ones(N)*L).astype('int') - # compute statistical moments - results['statistical_moments'][qoi_k] = {'mean': mean_k, - 'var': var_k, - 'std': std_k} - # compute all Sobol indices - #results['sobol_indices'][qoi_k] = self.get_sobol_indices(qoi_k, 'all') + for x in self.xi_d: + Map[k] = {'l':l, 'X':x, 'f':k} + k += 1 + #sparse grid + else: + + #all sparse grid multi indices l with |l| <= L + l_norm_le_L = self.sampler.compute_sparse_multi_idx(L, N) + + #if nested grids are used, compute the sparse grid of |l| <= L-1 + #which constains all points that will also be present in current grid + if self.nested == True and L > N: + H_Lm1_N = self.sampler.sparse_grid(L-1, N) - def compute_Q_diff(self, i, qoi, mean_f = []): + #loop over all multi indices + for l in l_norm_le_L: + + #colloc point of current level index l + X_l = [self.xi_1d[n][l[n]] for n in range(N)] + X_l = np.array(list(product(*X_l))) + + #if nested, remove points from X_l which are also in the previous sparse + #grid (|l| <= L-1) + if self.nested == True and L > N: + for x in X_l: + if list(x) in H_Lm1_N.tolist(): + X_l = np.delete(X_l, x, axis=0) + print("Grid is nested. Not using point", x, "at level", L) + #for each k, store level index l, collocation point x and index j of + #the code sample / collocation point in global grid xi_d + for x in X_l: + j = np.where((x == self.xi_d).all(axis=1))[0][0] + Map[k] = {'l':l, 'X':x, 'f':j} + k += 1; + + return Map + + def surrogate(self, qoi, x): """ - ======================================================================= - For every multi index i = (i1, i2, ..., id), Smolyak sums over - tensor products difference quadrature rules: - (Q^1_{i1} - Q^1_{i1-1}) X ... X (Q^1_{id) - Q^1_{id-1}) - Below this product is expanded into individual tensor products, each - of which is then computed as: - Q^1_{k1} X ... X Q^1_{kd} = sum...sum w_{k1}*...*w{kd}*f(x_{k1},...,x_{kd}) - ======================================================================= + Use sc_expansion UQP as a surrogate + + Parameters + ---------- + - qoi (str): name of the qoi + + Returns + ------- + the interpolated value of qoi at x (float, (N_qoi,)) + """ - #expand the multi-index indices of the tensor product - #(Q^1_{i1} - Q^1_{i1-1}) X ... X (Q^1_{id) - Q^1_{id-1}) - diff_idx = np.array(list(product(*[[k, -(k-1)] for k in i]))) - #Delta will be the sum of all expanded tensor products - #Q^1_{k1} X ... X Q^1_{kd} = sum...sum w_{k1}*...*w{kd}*f(x_{k1},...,x_{kd}) - Delta = 0.0 + return self.sc_expansion(self.L, 'interpolate', self.samples[qoi], x=x) - for diff in diff_idx: - #Q_0 = 0, so if present do not compute Q^1_{k1} X ... X Q^1_{kd} - if not np.in1d(0, diff): - X_k = [self.xi_1d[n][np.abs(diff)[n]] for n in range(self.N)] - X_k = np.array(list(product(*X_k))) - f_k = np.array([self.samples[qoi][x.tobytes()] for x in X_k]) - W_k = [self.wi_1d[n][np.abs(diff)[n]] for n in range(self.N)] - W_k = np.array(list(product(*W_k))) - W_k = np.prod(W_k, axis=1) - - #do some reshaping to align for multiplication - f_k = f_k.reshape([f_k.shape[0], self.N_qoi]) - W_k = W_k.reshape([W_k.shape[0], 1]) - - if mean_f == []: - Q_prod = np.sum(f_k*W_k, axis=0).T - else: - Q_prod = np.sum((f_k - mean_f)**2, axis=0).T - else: - Q_prod = 0.0 - - Delta += np.sign(np.prod(diff))*Q_prod - - return Delta - def get_moments(self, qoi): - """Compute the first two statistical moments. - + """ + Use sc_expansion UQP as a quadrature method + Parameters ---------- - qoi : str - column name of quantity of interest - + - qoi (str): name of the qoi + Returns ------- - float: - Mean of samples, using quad weights - float: - Variance of samples, using quad weights + - mean and variance of qoi (float (N_qoi,)) + """ - if self.sampler.sparse == False: - mean_f = np.zeros([self.N_qoi, 1]) - var_f = np.zeros([self.N_qoi, 1]) - - for k in range(self._number_of_samples): - sample_k = (self.samples[qoi][k]).values.reshape([self.N_qoi, 1]) - mean_f += sample_k * self.wi_d[k]#.prod() - - for k in range(self._number_of_samples): - sample_k = (self.samples[qoi][k]).values.reshape([self.N_qoi, 1]) - var_f += (sample_k - mean_f)**2 * self.wi_d[k]#.prod() - - return mean_f, var_f - else: - - i_norm_le_q = self.sampler.compute_sparse_multi_idx(self.L, self.N) - - #compute quadrature sum_i(Delta_i1 X ... X Delta__iN) - Delta = [] - for i in i_norm_le_q: - Delta.append(self.compute_Q_diff(i, qoi)) - - mean_f = np.sum(Delta, axis = 0) + #compute mean + mean_f = self.sc_expansion(self.L, 'quadrature', self.samples[qoi]) + + #compute variance + variance_samples = [] + for sample in self.samples[qoi]: + variance_samples.append((sample - mean_f)**2) - #compute quadrature sum_i(Delta_i1 X ... X Delta__iN) - Delta = [] - for i in i_norm_le_q: - Delta.append(self.compute_Q_diff(i, qoi, mean_f)) - - var_f = np.sum(Delta, axis = 0) - - return mean_f, var_f - - def get_sample_array(self, qoi): - """ - Returns an array of all samples of qoi + var_f = self.sc_expansion(self.L, 'quadrature', variance_samples) + + return mean_f, var_f + + def sc_expansion(self, L, goal, samples, **args): """ + ----------------------------------------- + This is the UQ Pattern for the SC method. + ----------------------------------------- - tmp = np.zeros([self._number_of_samples, self.N_qoi]) + Performs interpolation and quadrature for both full and sparse grids. - for k in range(self._number_of_samples): - tmp[k, :] = (self.samples[qoi][k]).values + For a qoi q, it computes the following tensor product: - return tmp - - def create_map(self, N, L): - l_norm_le_L = self.sampler.compute_sparse_multi_idx(L, N) + q \approx \sum_{l\in\Lambda} \Delta_{l}[q](x) - k = 0 - Map = {} - - for l in l_norm_le_L: - X_l = [self.xi_1d[n][l[n]] for n in range(N)] - X_l = np.array(list(product(*X_l))) - - for x in X_l: - j = np.where((x == self.xi_d).all(axis=1))[0][0] - Map[k] = {'l':l, 'X':x, 'f':j} - k += 1; - - return Map + where Delta_{l} is the difference at x between surrogates / quadratues + of level L and L-1. See e.g.: + + Dimitrios Loukrezis et. al., "Assessing the Performance of Leja and + Clenshaw-Curtis Collocation for Computational Electromagnetics with + Random Input Data." + + Parameters + ---------- - def sparse_surrogate(self, qoi, x, N, L): + - x (float (N,)): location in stochastic space at which to eval the surrogate + - L (int): max level of the surrogate - if L < N: + Returns + ------- + + surr (float, (N_qoi,)): the interpolated value of qoi at x if goal = 'interpolate' + the expected value of qoi if goal = 'quadrature' + """ + + #for L < L_min the surrogate is defined as zero + if L < self.L_min: return 0.0 surr = np.zeros(self.N_qoi) - for level in range(N, L+1): + #loop over all levels + for level in range(self.L_min, L+1): - Delta = np.zeros(self.N_qoi) + #contains the level multi-indices (l), colloc points x and samples + #indices f of the (sparse) grid Map = self.Map[level] - + + Delta = np.zeros(self.N_qoi) for k in Map.keys(): - q_k = self.samples2[qoi][Map[k]['f']] - s_k = q_k - self.sparse_surrogate(qoi, Map[k]['X'], N, level-1) + #the current code samples + q_k = samples[Map[k]['f']] + + #the hierarchical surplus between the code output q_k and the + #previous surrogate of level L-1 evaluated at the same location. + #Recursively computed. + s_k = q_k - self.sc_expansion(level-1, goal, samples, x = Map[k]['X']) + #the current level multi index (l_1,...,l_N) l = Map[k]['l'] idx = {} - for n in range(N): - # indices of current collocation point xi_d[k] in 1d - # collocation points + # indices of current collocation point (Map[k]['X'][n]), + # in corresponding 1d colloc points (self.xi_1d[n][l[n]]) + # These are the j of the 1D lagrange polynomials l_j(x), see + # lagrange_poly subroutine + for n in range(self.N): idx[n] = (self.xi_1d[n][l[n]] == Map[k]['X'][n]).nonzero()[0][0] - lagrange = [] - for n in range(N): - # values of Lagrange polynomials at x - lagrange.append(lagrange_poly(x[n], self.xi_1d[n][l[n]], idx[n])) - - Delta += s_k*np.prod(lagrange) + weight = [] + for n in range(self.N): + #interpolate + if goal == 'interpolate': + x = args['x'] + # add values of Lagrange polynomials at x + weight.append(lagrange_poly(x[n], self.xi_1d[n][l[n]], idx[n])) + #integrate + elif goal == 'quadrature': + #add quadrature weights + weight.append(self.wi_1d[n][l[n]][idx[n]]) + + #Delta is the interpolation of the hierarchical surplus + Delta += s_k*np.prod(weight) + surr += Delta return surr - - def surrogate(self, qoi, x): - """Use the SC expansion as a surrogate. - + + def compute_weights(self, rule = "G"): + + test = {} + + params = self.sampler.params_distribution + + for n in range(self.N): + for l in range(self.L_min, self.L+1): + + #TODO: compute integral of lagrange polynomials + xi_l, wi_l = cp.generate_quadrature(l, params[n], rule=rule) + + + def get_sample_array(self, qoi): + """ Parameters ---------- - qoi : str - Column name of quantity of interest - x : - - + - qoi (str): name of quantity of interest + Returns ------- - + - array of all samples of qoi """ - - # interpolated QoI - f_int = np.zeros([self.N_qoi, 1]) - - # list with the 1d collocation points of all uncertain parameters - # [self.all_vars[param]['xi_1d'] for param in self.all_vars.keys()] - C = self.xi - - # loop over all samples + + tmp = np.zeros([self._number_of_samples, self.N_qoi]) + for k in range(self._number_of_samples): - - idx = {} - for i in range(self.d): - # indices of current collocation point xi_d[k] in 1d - # collocation points - idx[i] = (C[i] == self.xi_d[k][i]).nonzero()[0] - - L = [] - for i in range(self.d): - # values of Lagrange polynomials at x - L.append(lagrange_poly(x[i], C[i], idx[i])) - - # current sample - qoi_k = self.samples[qoi][k].values.reshape([self.N_qoi, 1]) # .reshape(self.N_qoi) - - # surrogate: samples interpolated via Lagrange polynomials - f_int += qoi_k * np.prod(L) - - return f_int + tmp[k, :] = (self.samples[qoi][k]) + + return tmp + + def plot_grid(self): + """ + If N = 2 or N = 3 plot the (sparse) grid + """ + import matplotlib.pyplot as plt + + if self.N == 2: + fig = plt.figure() + ax = fig.add_subplot(111, xlabel=r'$x_1$', ylabel=r'$x_2$') + ax.plot(self.xi_d[:,0], self.xi_d[:,1], 'ro') + elif self.N == 3: + from mpl_toolkits.mplot3d import Axes3D + fig = plt.figure() + ax = fig.add_subplot(111, projection='3d', xlabel=r'$x_1$', \ + ylabel=r'$x_2$', zlabel=r'$x_3$') + ax.scatter(self.xi_d[:,0], self.xi_d[:,1], self.xi_d[:,2]) + else: + print('Will only plot for N = 2 or N = 3.') + + plt.tight_layout() + plt.show() # Start SC specific methods @@ -481,7 +491,7 @@ def compute_h(self, qoi, u, u_prime, xi_d_u, xi_d_u_prime, wi_d_u_prime): for i_up in range(S_u_prime): # collocation point to be evaluated - xi_s = np.zeros(self.d) + xi_s = np.zeros(self.N) # add the xi of u (at the correct location k) idx = 0 @@ -517,10 +527,10 @@ def get_sobol_indices(self, qoi, typ): """ # multi indices - U = list(range(self.d)) + U = list(range(self.N)) if typ == 'first_order': - P = list(powerset(U))[0:self.d + 1] + P = list(powerset(U))[0:self.N + 1] elif typ == 'all': # all indices u P = list(powerset(U)) @@ -586,9 +596,9 @@ def get_sobol_indices(self, qoi, typ): # End SC specific methods - def powerset(iterable): - """powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3) + """ + powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3) Taken from: https://docs.python.org/3/library/itertools.html#recipes @@ -605,22 +615,20 @@ def powerset(iterable): s = list(iterable) return chain.from_iterable(combinations(s, r) for r in range(len(s) + 1)) - def lagrange_poly(x, x_i, j): - """Lagrange polynomials used for interpolation + """ + Lagrange polynomials used for interpolation l_j(x) = product(x - x_m / x_j - x_m) with 0 <= m <= k and m !=j - TODO: Complete this docstring - Parameters ---------- - x : float + x : (float), location at which to compute the polynomial - x_i : list or array of float + x_i : list or array of float, nodes of the Lagrange polynomials - j : int + j : int, index of node at which l_j(x_j) = 1 Returns ------- @@ -638,4 +646,4 @@ def lagrange_poly(x, x_i, j): l_j *= nom / denom - return l_j + return l_j \ No newline at end of file diff --git a/easyvvuq/sampling/stochastic_collocation.py b/easyvvuq/sampling/stochastic_collocation.py index d02684f58..acc57d1f5 100644 --- a/easyvvuq/sampling/stochastic_collocation.py +++ b/easyvvuq/sampling/stochastic_collocation.py @@ -79,6 +79,30 @@ def __init__(self, self.sparse = sparse self.quad_sparse = sparse self.growth = growth + self.params_distribution = params_distribution + + #L = level of (sparse) grid + L = self.quad_order + #N = number of uncertain parameters + N = len(params_distribution) + + #for every dimension (parameter), create a hierachy of 1D + #quadrature rules of increasing order + self.xi_1d = {} + self.wi_1d = {} + + for n in range(N): + self.xi_1d[n] = {} + self.wi_1d[n] = {} + + for n in range(N): + for i in range(1, quad_order+1): + xi_i, wi_i = cp.generate_quadrature(i, + params_distribution[n], + rule=self.quad_rule, + growth=self.growth) + self.xi_1d[n][i] = xi_i[0] + self.wi_1d[n][i] = wi_i if sparse == False: # the nodes of the collocation grid @@ -89,44 +113,21 @@ def __init__(self, #sparse grid = a linear combination of tensor products of 1D rules #of different order. Use chaospy to compute these 1D quadrature rules else: - #q = level of sparse grid - q = self.quad_order - #N = number of uncertain parameters - N = len(params_distribution) - #q >= N must hold - if q < N: + #L >= N must hold + if L < N: print("*************************************************************") print("Level of sparse grid is lower than the dimension N (# params)") print("Increase level (quad_order) q such that q >= N") print("*************************************************************") import sys; sys.exit() - - #for every dimension (parameter), create a hierachy of 1D - #quadrature rules or increasing order - - self.xi_1d = {} - self.wi_1d = {} - - for n in range(N): - self.xi_1d[n] = {} - self.wi_1d[n] = {} - - for n in range(N): - for i in range(1, quad_order): - xi_i, wi_i = cp.generate_quadrature(i, - params_distribution[n], - rule=self.quad_rule, - growth=self.growth) - self.xi_1d[n][i] = xi_i[0] - self.wi_1d[n][i] = wi_i #create sparse grid of dimension N and level q using the 1d #rules in self.xi_1d - self.xi_d = self.sparse_grid_q_N(q, N) - self.q = q - self.N = N + self.xi_d = self.sparse_grid(L, N) + self.L = L + self.N = N self._number_of_samples = self.xi_d.shape[0] # Fast forward to specified count, if possible @@ -178,26 +179,26 @@ def get_restart_dict(self): ======================= """ - def sparse_grid_q_N(self, q, N): + def sparse_grid(self, L, N): """ - Compute an isotropic sparse grid H_q_N of N dimensions and level q + Compute an isotropic sparse grid H_L_N of N dimensions and level q """ - #multi-index i, such that |i| <= q - i_norm_le_q = self.compute_sparse_multi_idx(q, N) + #multi-index l, such that |l| <= L + l_norm_le_L = self.compute_sparse_multi_idx(L, N) - H_q_N = [] + H_L_N = [] #loop over all multi indices i - for i in i_norm_le_q: + for i in l_norm_le_L: #compute the tensor product of nodes indexed by i X_i = [self.xi_1d[n][i[n]] for n in range(N)] - H_q_N.append(list(product(*X_i))) + H_L_N.append(list(product(*X_i))) #flatten the list of lists - H_q_N = np.array(list(chain(*H_q_N))) + H_L_N = np.array(list(chain(*H_L_N))) #return unique nodes - return np.unique(H_q_N, axis=0) + return np.unique(H_L_N, axis=0) def compute_sparse_multi_idx(self, q, N): """ From 2ed51bc0af5255176bc2860ac67f69cc18430b13 Mon Sep 17 00:00:00 2001 From: wedeling Date: Fri, 16 Aug 2019 15:44:56 +0200 Subject: [PATCH 05/48] SC weights implemented + SC surrogate speedup --- easyvvuq/analysis/sc_analysis.py | 167 +++++++++++++++----- easyvvuq/sampling/stochastic_collocation.py | 9 +- 2 files changed, 130 insertions(+), 46 deletions(-) diff --git a/easyvvuq/analysis/sc_analysis.py b/easyvvuq/analysis/sc_analysis.py index e310084dd..de5a61fd9 100644 --- a/easyvvuq/analysis/sc_analysis.py +++ b/easyvvuq/analysis/sc_analysis.py @@ -71,10 +71,11 @@ def __init__(self, sampler=None, qoi_cols=None): #TODO: ARE THERE OTHERS WHICH ARE NESTED? genz_keister SHOULD BE, #BUT SEEMS NOT IMPLEMENTED IN MY cp VERSION (3.0.5) - if rule == 'leja' or rule == 'guass_patterson': - self.nested = True - elif rule == 'c' or rule == 'clenshaw_curtis' and growth == True: - self.nested = True + if self.sparse == True: + if rule == 'guass_patterson': + self.nested = True + elif rule == 'C' or rule == 'c' or rule == 'clenshaw_curtis' and growth == True: + self.nested = True def element_name(self): """Name for this element for logging purposes""" @@ -107,6 +108,20 @@ def analyse(self, data_frame=None): raise RuntimeError( "No data in data frame passed to analyse element") + #the maximum level (quad order) of the (sparse) grid + self.L = self.sampler.L + + #the number of uncertain parameters + self.N = self.sampler.N + + #if L < L_min: quadratures and interpolations are zero + #For full tensor grid: there is only one level: L_min = L + if self.sparse == False: + self.L_min = self.L + #For sparse grid: multiple levels, L >= N must hold + else: + self.L_min = self.N + # Chaospy computation of 1D weights - USED IN SOBOL SUBROUTINES #REMOVE LATER IN FAVOUR OF xi_1D, wi_1d xi = [] @@ -119,37 +134,22 @@ def analyse(self, data_frame=None): self.xi = xi self.wi = wi - # Compute tensor product nodes and weights - # REMOVE WHEN MOMENT SUBROUTINE USES SC WEIGHTS - _, self.wi_d = cp.generate_quadrature(order=self.sampler.quad_order, - domain=self.sampler.joint_dist, - rule=self.sampler.quad_rule) - self.xi_d = self.sampler.xi_d self.xi_1d = self.sampler.xi_1d - self.wi_1d = self.sampler.wi_1d - - #the maximum level (quad order) of the (sparse) grid - self.L = self.sampler.L - - #the number of uncertain parameters - self.N = self.sampler.N - - #if L < L_min: quadratures and interpolations are zero - #For full tensor grid: there is only one level: L_min = L - if self.sparse == False: - self.L_min = self.L - #For sparse grid: multiple levels, L >= N must hold - else: - self.L_min = self.N + self.wi_1d = self.compute_SC_weights(rule=self.sampler.quad_rule) #per level, map a unique index k to all (level multi indices, colloc points) #combinations. Will differ for sparse or full tensor grids. #All interpolation/quadrature subroutines loop over the entries in Map self.Map = {} + self.surr_lm1 = {} + for level in range(self.L_min, self.L+1): self.Map[level] = self.create_map(self.N, level) + self.clear_surr_lm1('interpolate') + self.clear_surr_lm1('quadrature') + # Extract output values for each quantity of interest from Dataframe qoi_cols = self.qoi_cols samples = {k: [] for k in qoi_cols} @@ -199,6 +199,8 @@ def create_map(self, N, L): k = 0 Map = {} + print('Creating multi-index map for level', L, '...') + #full tensor product if self.sparse == False: @@ -231,7 +233,8 @@ def create_map(self, N, L): for x in X_l: if list(x) in H_Lm1_N.tolist(): X_l = np.delete(X_l, x, axis=0) - print("Grid is nested. Not using point", x, "at level", L) + print("Grid is nested. Not using point", \ + ["%.4f" %xx for xx in x], "at level", L) #for each k, store level index l, collocation point x and index j of #the code sample / collocation point in global grid xi_d @@ -239,6 +242,8 @@ def create_map(self, N, L): j = np.where((x == self.xi_d).all(axis=1))[0][0] Map[k] = {'l':l, 'X':x, 'f':j} k += 1; + + print('done.') return Map @@ -271,16 +276,19 @@ def get_moments(self, qoi): - mean and variance of qoi (float (N_qoi,)) """ - + goal = 'quadrature' #compute mean - mean_f = self.sc_expansion(self.L, 'quadrature', self.samples[qoi]) + mean_f = self.sc_expansion(self.L, goal, self.samples[qoi]) + + #clear the quadrature results in surr_lm1 + self.clear_surr_lm1(goal) #compute variance variance_samples = [] for sample in self.samples[qoi]: variance_samples.append((sample - mean_f)**2) - var_f = self.sc_expansion(self.L, 'quadrature', variance_samples) + var_f = self.sc_expansion(self.L, goal, variance_samples) return mean_f, var_f @@ -290,7 +298,7 @@ def sc_expansion(self, L, goal, samples, **args): This is the UQ Pattern for the SC method. ----------------------------------------- - Performs interpolation and quadrature for both full and sparse grids. + Can perform interpolation and quadrature for both full and sparse grids. For a qoi q, it computes the following tensor product: @@ -335,11 +343,19 @@ def sc_expansion(self, L, goal, samples, **args): #the current code samples q_k = samples[Map[k]['f']] - #the hierarchical surplus between the code output q_k and the + #the hierarchical surplus (s_k) between the code output q_k and the #previous surrogate of level L-1 evaluated at the same location. #Recursively computed. - s_k = q_k - self.sc_expansion(level-1, goal, samples, x = Map[k]['X']) - + + if k in self.surr_lm1[goal][level]: + #print('surrogate already computed') + surr_lm1 = self.surr_lm1[goal][level][k] + else: + surr_lm1 = self.sc_expansion(level-1, goal, samples, x = Map[k]['X']) + self.surr_lm1[goal][level][k] = surr_lm1 + + s_k = q_k - surr_lm1 + #the current level multi index (l_1,...,l_N) l = Map[k]['l'] @@ -370,18 +386,85 @@ def sc_expansion(self, L, goal, samples, **args): return surr - def compute_weights(self, rule = "G"): + def clear_surr_lm1(self, ID): + """ + Clears the interpolation or quadrature results in surr_lm1[ID]. - test = {} + surr_lm1 is a dictionary used to store surrogate results at + previous level (l-1). Used to avoid recomputing the surrogate + in the recursive sc_expansion subroutine. - params = self.sampler.params_distribution + surr_lm1['interpolation'][l][k] stores the interpolation results + of level l-1 at collocation point X_k - for n in range(self.N): - for l in range(self.L_min, self.L+1): - - #TODO: compute integral of lagrange polynomials - xi_l, wi_l = cp.generate_quadrature(l, params[n], rule=rule) - + Parameters + ---------- + - ID (str): either 'interpolate' or 'quadrature' + + """ + self.surr_lm1[ID] = {} + for level in range(self.L_min, self.L+1): + self.surr_lm1[ID][level] = {} + + def compute_SC_weights(self, rule): + """ + Computes the 1D quadrature weights w_j of the SC expansion: + + w_j = int L_j(x)p(x) dx (1) + + Here L_j is a Lagrange polynomial of the SC expansion. + + Parameters + ---------- + - rule ("str"): chaospy quadrature rule used to compute (1), + + + Returns + ------- + - wi_1d (dict): wi_1d[n][l] gives an array + of quadrature weigths for the n-th parameter at level l. + + IMPORTANT: + If rule is the same as the rule used to compute the SC + collocation points, these weights will equal the weights + computed by chaospy, since L_j(x_k) = 1 when j=k and 0 + for the rest. This is the default setting. + """ + + #no need to recompute weights + if rule == self.sampler.quadrature_rule: + return self.sampler.wi_1d + #recompute weights - generally not used + else: + wi_1d = {} + + params = self.sampler.params_distribution + + for n in range(self.N): + #1d weights for n-th parameter + wi_1d[n] = {} + #loop over all level of collocation method + for level in range(1, self.L+1): + #current SC nodes over dimension n and level + xi_1d = self.xi_1d[n][level] + wi_1d[n][level] = np.zeros(xi_1d.size) + + #generate a quadrature rule to compute the SC weights + xi_quad, wi_quad = cp.generate_quadrature(level, params[n], rule=rule) + xi_quad = xi_quad[0] + + #compute integral of the lagrange polynomial through xi_1d, weighted + #by the input distributions: + #w_j = int L_j(xi) p(xi) dxi j = 1,..,xi_1d.size + for j in range(xi_1d.size): + #values of L_i(xi_quad) + lagrange_quad = np.zeros(xi_quad.size) + for i in range(xi_quad.size): + lagrange_quad[i] = lagrange_poly(xi_quad[i], xi_1d, j) + #quadrature + wi_1d[n][level][j] = np.sum(lagrange_quad*wi_quad) + + return wi_1d def get_sample_array(self, qoi): """ diff --git a/easyvvuq/sampling/stochastic_collocation.py b/easyvvuq/sampling/stochastic_collocation.py index acc57d1f5..6023c892d 100644 --- a/easyvvuq/sampling/stochastic_collocation.py +++ b/easyvvuq/sampling/stochastic_collocation.py @@ -32,7 +32,7 @@ class SCSampler(BaseSamplingElement, sampler_name="sc_sampler"): def __init__(self, vary=None, - quad_order=4, + polynomial_order=4, quadrature_rule="G", count=0, growth=False, @@ -74,7 +74,7 @@ def __init__(self, # The quadrature information: order, rule and sparsity #self.quad_order = polynomial_order - self.quad_order = quad_order + self.quad_order = polynomial_order self.quad_rule = quadrature_rule self.sparse = sparse self.quad_sparse = sparse @@ -96,11 +96,12 @@ def __init__(self, self.wi_1d[n] = {} for n in range(N): - for i in range(1, quad_order+1): + for i in range(1, self.quad_order+1): xi_i, wi_i = cp.generate_quadrature(i, params_distribution[n], rule=self.quad_rule, - growth=self.growth) + growth=self.growth, + normalize=True) self.xi_1d[n][i] = xi_i[0] self.wi_1d[n][i] = wi_i From b029d568a6fbacfc0e8ddb2439cd625012a132c5 Mon Sep 17 00:00:00 2001 From: orbitfold Date: Mon, 19 Aug 2019 16:01:07 +0200 Subject: [PATCH 06/48] Added test_add_runs --- tests/test_db.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/tests/test_db.py b/tests/test_db.py index 415a7a3b3..ee45f79fc 100644 --- a/tests/test_db.py +++ b/tests/test_db.py @@ -3,7 +3,8 @@ import easyvvuq from easyvvuq.constants import default_campaign_prefix from easyvvuq.db.sql import CampaignDB -from easyvvuq.data_structs import CampaignInfo +from easyvvuq.data_structs import CampaignInfo, RunInfo + @pytest.fixture def campaign(tmp_path): @@ -16,8 +17,14 @@ def campaign(tmp_path): campaign.tmp_path = tmp_path return campaign + def test_db_file_created(campaign): assert(os.path.isfile('{}/test.sqlite'.format(campaign.tmp_path))) -def test_add_runs(tmp_path): - pass + +def test_add_runs(campaign): + run1 = RunInfo('run1', 'test', '.', 1, {'a' : 1}, 1, 1) + run2 = RunInfo('run2', 'test', '.', 1, {'a' : 1}, 1, 1) + campaign.add_runs([run1, run2]) + campaign.get_run_status('run1', 'test', 1) + campaign.get_run_status('run2', 'test', 1) From b140bc2ee729e31bf6ed9b46fca6444c7e8c3c18 Mon Sep 17 00:00:00 2001 From: orbitfold Date: Tue, 13 Aug 2019 13:00:50 +0200 Subject: [PATCH 07/48] Initial campaigndb test --- tests/test_db.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) create mode 100644 tests/test_db.py diff --git a/tests/test_db.py b/tests/test_db.py new file mode 100644 index 000000000..2679203cb --- /dev/null +++ b/tests/test_db.py @@ -0,0 +1,15 @@ +import pytest +import os.path +import easyvvuq +from easyvvuq.constants import default_campaign_prefix +from easyvvuq.db.sql import CampaignDB +from easyvvuq.data_structs import CampaignInfo + +def test_file_created(tmp_path): + info = CampaignInfo( + name='test', + campaign_dir_prefix=default_campaign_prefix, + easyvvuq_version=easyvvuq.__version__, + campaign_dir='.') + campaign = CampaignDB(location='sqlite:///{}/test.sqlite'.format(tmp_path), new_campaign=True, name='test', info=info) + assert(os.path.isfile('{}/test.sqlite'.format(tmp_path))) From 05eb2a77555e49648bedf40462f1a536bb85f516 Mon Sep 17 00:00:00 2001 From: orbitfold Date: Wed, 14 Aug 2019 15:35:15 +0200 Subject: [PATCH 08/48] Added a campaign fixture: --- tests/test_db.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/tests/test_db.py b/tests/test_db.py index 2679203cb..415a7a3b3 100644 --- a/tests/test_db.py +++ b/tests/test_db.py @@ -5,11 +5,19 @@ from easyvvuq.db.sql import CampaignDB from easyvvuq.data_structs import CampaignInfo -def test_file_created(tmp_path): +@pytest.fixture +def campaign(tmp_path): info = CampaignInfo( name='test', campaign_dir_prefix=default_campaign_prefix, easyvvuq_version=easyvvuq.__version__, campaign_dir='.') campaign = CampaignDB(location='sqlite:///{}/test.sqlite'.format(tmp_path), new_campaign=True, name='test', info=info) - assert(os.path.isfile('{}/test.sqlite'.format(tmp_path))) + campaign.tmp_path = tmp_path + return campaign + +def test_db_file_created(campaign): + assert(os.path.isfile('{}/test.sqlite'.format(campaign.tmp_path))) + +def test_add_runs(tmp_path): + pass From 43f774d23a64970116808084bb1e15f8baeb72db Mon Sep 17 00:00:00 2001 From: orbitfold Date: Mon, 19 Aug 2019 16:01:07 +0200 Subject: [PATCH 09/48] Added test_add_runs --- tests/test_db.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/tests/test_db.py b/tests/test_db.py index 415a7a3b3..ee45f79fc 100644 --- a/tests/test_db.py +++ b/tests/test_db.py @@ -3,7 +3,8 @@ import easyvvuq from easyvvuq.constants import default_campaign_prefix from easyvvuq.db.sql import CampaignDB -from easyvvuq.data_structs import CampaignInfo +from easyvvuq.data_structs import CampaignInfo, RunInfo + @pytest.fixture def campaign(tmp_path): @@ -16,8 +17,14 @@ def campaign(tmp_path): campaign.tmp_path = tmp_path return campaign + def test_db_file_created(campaign): assert(os.path.isfile('{}/test.sqlite'.format(campaign.tmp_path))) -def test_add_runs(tmp_path): - pass + +def test_add_runs(campaign): + run1 = RunInfo('run1', 'test', '.', 1, {'a' : 1}, 1, 1) + run2 = RunInfo('run2', 'test', '.', 1, {'a' : 1}, 1, 1) + campaign.add_runs([run1, run2]) + campaign.get_run_status('run1', 'test', 1) + campaign.get_run_status('run2', 'test', 1) From ef444191496fab41c25b8d6c0309f1d3e4b96c89 Mon Sep 17 00:00:00 2001 From: orbitfold Date: Wed, 21 Aug 2019 16:00:40 +0200 Subject: [PATCH 10/48] test set statuses --- easyvvuq/db/sql.py | 2 +- tests/test_db.py | 10 +++++++--- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/easyvvuq/db/sql.py b/easyvvuq/db/sql.py index f8479386e..53497cdee 100644 --- a/easyvvuq/db/sql.py +++ b/easyvvuq/db/sql.py @@ -479,7 +479,7 @@ def get_run_status(self, run_name, campaign=None, sampler=None): if campaign: filter_options['campaign'] = campaign if sampler: - filter_options['sample'] = sampler + filter_options['sampler'] = sampler selected = self.session.query(RunTable).filter_by(**filter_options) diff --git a/tests/test_db.py b/tests/test_db.py index ee45f79fc..cf6b68874 100644 --- a/tests/test_db.py +++ b/tests/test_db.py @@ -4,6 +4,7 @@ from easyvvuq.constants import default_campaign_prefix from easyvvuq.db.sql import CampaignDB from easyvvuq.data_structs import CampaignInfo, RunInfo +from easyvvuq.constants import Status @pytest.fixture @@ -22,9 +23,12 @@ def test_db_file_created(campaign): assert(os.path.isfile('{}/test.sqlite'.format(campaign.tmp_path))) -def test_add_runs(campaign): +def test_get_and_set_status(campaign): run1 = RunInfo('run1', 'test', '.', 1, {'a' : 1}, 1, 1) run2 = RunInfo('run2', 'test', '.', 1, {'a' : 1}, 1, 1) campaign.add_runs([run1, run2]) - campaign.get_run_status('run1', 'test', 1) - campaign.get_run_status('run2', 'test', 1) + assert(campaign.get_run_status('Run_1') == Status.NEW) + assert(campaign.get_run_status('Run_2') == Status.NEW) + campaign.set_run_statuses(['Run_1'], Status.ENCODED) + assert(campaign.get_run_status('Run_1') == Status.ENCODED) + From 36059772c3cf54006d0c4cca70c9795949f328ee Mon Sep 17 00:00:00 2001 From: orbitfold Date: Wed, 21 Aug 2019 16:03:56 +0200 Subject: [PATCH 11/48] test set statuses --- tests/test_db.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/tests/test_db.py b/tests/test_db.py index cf6b68874..12b18b0c1 100644 --- a/tests/test_db.py +++ b/tests/test_db.py @@ -32,3 +32,10 @@ def test_get_and_set_status(campaign): campaign.set_run_statuses(['Run_1'], Status.ENCODED) assert(campaign.get_run_status('Run_1') == Status.ENCODED) + +def test_too_many_runs_bug(campaign): + runs = [RunInfo('run', 'test', '.', 1, {'a' : 1}, 1, 1) for _ in range(2000)] + run_names = ['Run_{}'.format(i) for i in range(1, 2001)] + campaign.add_runs(runs) + + From 7fa5967ceb30f93901f239ffa6b278615acd63eb Mon Sep 17 00:00:00 2001 From: orbitfold Date: Wed, 21 Aug 2019 16:08:54 +0200 Subject: [PATCH 12/48] test set statuses --- tests/test_db.py | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/tests/test_db.py b/tests/test_db.py index 12b18b0c1..c40e34fa4 100644 --- a/tests/test_db.py +++ b/tests/test_db.py @@ -24,18 +24,9 @@ def test_db_file_created(campaign): def test_get_and_set_status(campaign): - run1 = RunInfo('run1', 'test', '.', 1, {'a' : 1}, 1, 1) - run2 = RunInfo('run2', 'test', '.', 1, {'a' : 1}, 1, 1) - campaign.add_runs([run1, run2]) - assert(campaign.get_run_status('Run_1') == Status.NEW) - assert(campaign.get_run_status('Run_2') == Status.NEW) - campaign.set_run_statuses(['Run_1'], Status.ENCODED) - assert(campaign.get_run_status('Run_1') == Status.ENCODED) - - -def test_too_many_runs_bug(campaign): runs = [RunInfo('run', 'test', '.', 1, {'a' : 1}, 1, 1) for _ in range(2000)] run_names = ['Run_{}'.format(i) for i in range(1, 2001)] campaign.add_runs(runs) - - + assert(all([campaign.get_run_status(name) == Status.NEW for name in run_names])) + campaign.set_run_statuses(run_names, Status.ENCODED) + assert(all([campaign.get_run_status(name) == Status.ENCODED for name in run_names])) From d80205eb6dcc4ae74f980b3706c60b029a43937c Mon Sep 17 00:00:00 2001 From: orbitfold Date: Tue, 27 Aug 2019 16:31:09 +0200 Subject: [PATCH 13/48] Added test_gun_num_runs --- tests/test_db.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/tests/test_db.py b/tests/test_db.py index c40e34fa4..d135a26fd 100644 --- a/tests/test_db.py +++ b/tests/test_db.py @@ -16,6 +16,9 @@ def campaign(tmp_path): campaign_dir='.') campaign = CampaignDB(location='sqlite:///{}/test.sqlite'.format(tmp_path), new_campaign=True, name='test', info=info) campaign.tmp_path = tmp_path + runs = [RunInfo('run', 'test', '.', 1, {'a' : 1}, 1, 1) for _ in range(1010)] + run_names = ['Run_{}'.format(i) for i in range(1, 1011)] + campaign.add_runs(runs) return campaign @@ -24,9 +27,12 @@ def test_db_file_created(campaign): def test_get_and_set_status(campaign): - runs = [RunInfo('run', 'test', '.', 1, {'a' : 1}, 1, 1) for _ in range(2000)] - run_names = ['Run_{}'.format(i) for i in range(1, 2001)] - campaign.add_runs(runs) + run_names = ['Run_{}'.format(i) for i in range(1, 1011)] assert(all([campaign.get_run_status(name) == Status.NEW for name in run_names])) campaign.set_run_statuses(run_names, Status.ENCODED) assert(all([campaign.get_run_status(name) == Status.ENCODED for name in run_names])) + + +def test_get_num_runs(campaign): + assert(campaign.get_num_runs() == 1010) + From fd1f1bdbddde9cd9252c5acd31e338c23ad62259 Mon Sep 17 00:00:00 2001 From: orbitfold Date: Tue, 27 Aug 2019 16:44:33 +0200 Subject: [PATCH 14/48] Added test_app --- tests/test_db.py | 39 ++++++++++++++++++++++++++++++++++++--- 1 file changed, 36 insertions(+), 3 deletions(-) diff --git a/tests/test_db.py b/tests/test_db.py index d135a26fd..4e355e9ab 100644 --- a/tests/test_db.py +++ b/tests/test_db.py @@ -1,9 +1,9 @@ import pytest import os.path -import easyvvuq +import easyvvuq as uq from easyvvuq.constants import default_campaign_prefix from easyvvuq.db.sql import CampaignDB -from easyvvuq.data_structs import CampaignInfo, RunInfo +from easyvvuq.data_structs import CampaignInfo, RunInfo, AppInfo from easyvvuq.constants import Status @@ -12,13 +12,42 @@ def campaign(tmp_path): info = CampaignInfo( name='test', campaign_dir_prefix=default_campaign_prefix, - easyvvuq_version=easyvvuq.__version__, + easyvvuq_version=uq.__version__, campaign_dir='.') campaign = CampaignDB(location='sqlite:///{}/test.sqlite'.format(tmp_path), new_campaign=True, name='test', info=info) campaign.tmp_path = tmp_path runs = [RunInfo('run', 'test', '.', 1, {'a' : 1}, 1, 1) for _ in range(1010)] run_names = ['Run_{}'.format(i) for i in range(1, 1011)] campaign.add_runs(runs) + app_info = AppInfo('test', uq.ParamsSpecification({ + "temp_init": { + "type": "float", + "min": 0.0, + "max": 100.0, + "default": 95.0}, + "kappa": { + "type": "float", + "min": 0.0, + "max": 0.1, + "default": 0.025}, + "t_env": { + "type": "float", + "min": 0.0, + "max": 40.0, + "default": 15.0}, + "out_file": { + "type": "string", + "default": "output.csv"} + }), + None, + uq.encoders.GenericEncoder( + template_fname='tests/cooling/cooling.template', + delimiter='$', + target_filename='cooling_in.json'), + uq.decoders.SimpleCSV(target_filename='output.csv', + output_columns=["te", "ti"], + header=0)) + campaign.add_app(app_info) return campaign @@ -36,3 +65,7 @@ def test_get_and_set_status(campaign): def test_get_num_runs(campaign): assert(campaign.get_num_runs() == 1010) + +def test_app(campaign): + app_dict = campaign.app('test') + assert(app_dict['name'] == 'test') From 7e659b741a31f0ec2f1c5a5bf07f4697855dcff7 Mon Sep 17 00:00:00 2001 From: David W Wright Date: Tue, 27 Aug 2019 16:33:08 +0100 Subject: [PATCH 15/48] Logo text updates --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 3c37e0017..92fb640b4 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -EassyVVUQ icon +EasyVVUQ icon # EasyVVUQ From 38ced42aad673e09923d4f0e8379bab8acededff Mon Sep 17 00:00:00 2001 From: David W Wright Date: Tue, 27 Aug 2019 17:02:18 +0100 Subject: [PATCH 16/48] Logo redirect --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 92fb640b4..d9b9186a8 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -EasyVVUQ icon +EasyVVUQ icon # EasyVVUQ From 3c679bea0335ae47bddc3fa2895c065e15bc3a95 Mon Sep 17 00:00:00 2001 From: wedeling Date: Tue, 27 Aug 2019 23:19:36 +0200 Subject: [PATCH 17/48] sparse grid interpolation and sobol indices implemented --- easyvvuq/analysis/sc_analysis.py | 386 +++++++++++--------- easyvvuq/sampling/stochastic_collocation.py | 50 ++- 2 files changed, 253 insertions(+), 183 deletions(-) diff --git a/easyvvuq/analysis/sc_analysis.py b/easyvvuq/analysis/sc_analysis.py index de5a61fd9..2d88456b9 100644 --- a/easyvvuq/analysis/sc_analysis.py +++ b/easyvvuq/analysis/sc_analysis.py @@ -64,19 +64,6 @@ def __init__(self, sampler=None, qoi_cols=None): self._number_of_samples = sampler._number_of_samples self.sparse = sampler.sparse - #determine if quadrature rule is nested, default assumption is False - self.nested = False - rule = self.sampler.quadrature_rule - growth = self.sampler.growth - - #TODO: ARE THERE OTHERS WHICH ARE NESTED? genz_keister SHOULD BE, - #BUT SEEMS NOT IMPLEMENTED IN MY cp VERSION (3.0.5) - if self.sparse == True: - if rule == 'guass_patterson': - self.nested = True - elif rule == 'C' or rule == 'c' or rule == 'clenshaw_curtis' and growth == True: - self.nested = True - def element_name(self): """Name for this element for logging purposes""" return "SC_Analysis" @@ -118,21 +105,11 @@ def analyse(self, data_frame=None): #For full tensor grid: there is only one level: L_min = L if self.sparse == False: self.L_min = self.L + self.l_norm = np.array([[self.L, self.L]]) #For sparse grid: multiple levels, L >= N must hold else: - self.L_min = self.N - - # Chaospy computation of 1D weights - USED IN SOBOL SUBROUTINES - #REMOVE LATER IN FAVOUR OF xi_1D, wi_1d - xi = [] - wi = [] - for dist in self.sampler.vary.get_values(): - xi_i, wi_i = cp.generate_quadrature( - self.sampler.quad_order, dist, rule=self.sampler.quad_rule) - xi.append(xi_i.flatten()) - wi.append(wi_i.flatten()) - self.xi = xi - self.wi = wi + self.L_min = 1 + self.l_norm = self.sampler.compute_sparse_multi_idx(self.L, self.N) self.xi_d = self.sampler.xi_d self.xi_1d = self.sampler.xi_1d @@ -144,11 +121,12 @@ def analyse(self, data_frame=None): self.Map = {} self.surr_lm1 = {} + self.foo = [] + for level in range(self.L_min, self.L+1): self.Map[level] = self.create_map(self.N, level) - self.clear_surr_lm1('interpolate') - self.clear_surr_lm1('quadrature') + self.clear_surr_lm1() # Extract output values for each quantity of interest from Dataframe qoi_cols = self.qoi_cols @@ -165,7 +143,7 @@ def analyse(self, data_frame=None): results = {'statistical_moments': {}, 'sobol_indices': {k: {} for k in self.qoi_cols}} - # Compute descriptive statistics for each quantity of interest + #Compute descriptive statistics for each quantity of interest for qoi_k in qoi_cols: mean_k, var_k = self.get_moments(qoi_k) std_k = var_k**0.5 @@ -174,8 +152,8 @@ def analyse(self, data_frame=None): results['statistical_moments'][qoi_k] = {'mean': mean_k, 'var': var_k, 'std': std_k} -# # compute all Sobol indices -# results['sobol_indices'][qoi_k] = self.get_sobol_indices(qoi_k, 'all') + # compute all Sobol indices + results['sobol_indices'][qoi_k] = self.get_sobol_indices(qoi_k, 'all') return results @@ -215,11 +193,6 @@ def create_map(self, N, L): #all sparse grid multi indices l with |l| <= L l_norm_le_L = self.sampler.compute_sparse_multi_idx(L, N) - #if nested grids are used, compute the sparse grid of |l| <= L-1 - #which constains all points that will also be present in current grid - if self.nested == True and L > N: - H_Lm1_N = self.sampler.sparse_grid(L-1, N) - #loop over all multi indices for l in l_norm_le_L: @@ -227,17 +200,7 @@ def create_map(self, N, L): X_l = [self.xi_1d[n][l[n]] for n in range(N)] X_l = np.array(list(product(*X_l))) - #if nested, remove points from X_l which are also in the previous sparse - #grid (|l| <= L-1) - if self.nested == True and L > N: - for x in X_l: - if list(x) in H_Lm1_N.tolist(): - X_l = np.delete(X_l, x, axis=0) - print("Grid is nested. Not using point", \ - ["%.4f" %xx for xx in x], "at level", L) - - #for each k, store level index l, collocation point x and index j of - #the code sample / collocation point in global grid xi_d + for x in X_l: j = np.where((x == self.xi_d).all(axis=1))[0][0] Map[k] = {'l':l, 'X':x, 'f':j} @@ -247,7 +210,7 @@ def create_map(self, N, L): return Map - def surrogate(self, qoi, x): + def surrogate(self, qoi, x, **kwargs): """ Use sc_expansion UQP as a surrogate @@ -261,12 +224,98 @@ def surrogate(self, qoi, x): """ - return self.sc_expansion(self.L, 'interpolate', self.samples[qoi], x=x) + if 'L' in kwargs: + L = kwargs['L'] + else: + L = self.L + + return self.sc_expansion(L, self.samples[qoi], x=x) + + def quadrature(self, qoi, **kwargs): + """ + Computes a (Smolyak) quadrature + + Parameters + ---------- + - qoi (str): name of the qoi + + - samples (optional in kwargs): Default: compute the mean + by setting samples = self.samples. To compute the variance, + set samples = (self.samples - mean)**2 + """ + + + if 'samples' in kwargs: + samples = kwargs['samples'] + else: + samples = self.samples[qoi] + + Delta = np.zeros([self.l_norm.shape[0], self.N_qoi]) + idx = 0 + for l in self.l_norm: + #compute the Delta Q := + #(Q^1_l_1 - Q^1_{l_1 - 1}) X ... X (Q^1_{l_N} - Q^1_{L_N - 1}) + #tensor product + Delta[idx, :] = self.compute_Q_diff(l, samples) + idx += 1 + + quadrature_approx = np.sum(Delta, axis=0) + + return quadrature_approx + + def compute_Q_diff(self, l, samples, **kwargs): - def get_moments(self, qoi): """ - Use sc_expansion UQP as a quadrature method + ======================================================================= + For every multi index l = (l1, l2, ..., ld), Smolyak sums over + tensor products difference quadrature rules: + (Q^1_{l1} - Q^1_{l1-1}) X ... X (Q^1_{lN) - Q^1_{lN-1}) + Below this product is expanded into individual tensor products, each + of which is then computed as: + Q^1_{k1} X ... X Q^1_{kN} = sum...sum w_{k1}*...*w{kN}*f(x_{k1},...,x_{kN}) + ======================================================================= + """ + + #expand the multi-index indices of the tensor product + #(Q^1_{i1} - Q^1_{i1-1}) X ... X (Q^1_{id) - Q^1_{id-1}) + diff_idx = np.array(list(product(*[[k, -(k-1)] for k in l]))) + #Delta will be the sum of all expanded tensor products + #Q^1_{k1} X ... X Q^1_{kd} = sum...sum w_{k1}*...*w{kN}*f(x_{k1},...,x_{kd}) + Delta = 0.0 + + #each diff contains the level indices of to a single + #Q^1_{k1} X ... X Q^1_{kN} product + for diff in diff_idx: + #Q_0 = 0, so if present do not compute Q^1_{k1} X ... X Q^1_{kd} + if not np.min(np.abs(diff)) < self.L_min: + + #compute the tensor product of parameter and weight values + X_k = []; W_k = []; + for n in range(self.N): + X_k.append(self.xi_1d[n][np.abs(diff)[n]]) + W_k.append(self.wi_1d[n][np.abs(diff)[n]]) + + X_k = np.array(list(product(*X_k))) + W_k = np.array(list(product(*W_k))) + W_k = np.prod(W_k, axis=1) + W_k = W_k.reshape([W_k.shape[0], 1]) + + #find corresponding code values + f_k = [] + for x in X_k: + j = np.where((x == self.xi_d).all(axis=1))[0][0] + f_k.append(samples[j]) + f_k = np.array(f_k).reshape([len(X_k), self.N_qoi]) + + #quadrature of Q^1_{k1} X ... X Q^1_{kN} product + Q_prod = np.sum(f_k*W_k, axis=0).T + Delta += np.sign(np.prod(diff))*Q_prod + + return Delta + + def get_moments(self, qoi): + """ Parameters ---------- - qoi (str): name of the qoi @@ -276,29 +325,26 @@ def get_moments(self, qoi): - mean and variance of qoi (float (N_qoi,)) """ - goal = 'quadrature' - #compute mean - mean_f = self.sc_expansion(self.L, goal, self.samples[qoi]) - #clear the quadrature results in surr_lm1 - self.clear_surr_lm1(goal) + #compute mean + mean_f = self.quadrature(qoi) #compute variance variance_samples = [] for sample in self.samples[qoi]: variance_samples.append((sample - mean_f)**2) - - var_f = self.sc_expansion(self.L, goal, variance_samples) + + var_f = self.quadrature(qoi, samples=variance_samples) return mean_f, var_f - def sc_expansion(self, L, goal, samples, **args): + def sc_expansion(self, L, samples, **kwargs): """ ----------------------------------------- This is the UQ Pattern for the SC method. ----------------------------------------- - Can perform interpolation and quadrature for both full and sparse grids. + Can perform interpolation for both full and sparse grids. For a qoi q, it computes the following tensor product: @@ -320,8 +366,7 @@ def sc_expansion(self, L, goal, samples, **args): Returns ------- - surr (float, (N_qoi,)): the interpolated value of qoi at x if goal = 'interpolate' - the expected value of qoi if goal = 'quadrature' + surr (float, (N_qoi,)): the interpolated value of qoi at x """ #for L < L_min the surrogate is defined as zero @@ -330,71 +375,65 @@ def sc_expansion(self, L, goal, samples, **args): surr = np.zeros(self.N_qoi) - #loop over all levels - for level in range(self.L_min, L+1): + #contains the level multi-indices (l), colloc points x and samples + #indices f of the (sparse) grid + Map = self.Map[L] + + for k in Map.keys(): + + #the current code samples + q_k = samples[Map[k]['f']] + + #the current level multi index (l_1,...,l_N) + l = Map[k]['l'] + + #the hierarchical surplus (s_k) between the code output q_k and the + #previous surrogate of level L-1 evaluated at the same location. + #Recursively computed. + + if self.sparse == True: + Lm1 = np.sum(l)-1 + else: + Lm1 = self.L-1 + + if k in self.surr_lm1[L]: +# print('surrogate already computed') + surr_lm1 = self.surr_lm1[L][k] + else: + surr_lm1 = self.sc_expansion(Lm1, samples, x = Map[k]['X']) + self.surr_lm1[L][k] = surr_lm1 + + s_k = q_k - surr_lm1 - #contains the level multi-indices (l), colloc points x and samples - #indices f of the (sparse) grid - Map = self.Map[level] + idx = {} + # indices of current collocation point (Map[k]['X'][n]), + # in corresponding 1d colloc points (self.xi_1d[n][l[n]]) + # These are the j of the 1D lagrange polynomials l_j(x), see + # lagrange_poly subroutine + for n in range(self.N): + idx[n] = (self.xi_1d[n][l[n]] == Map[k]['X'][n]).nonzero()[0][0] - Delta = np.zeros(self.N_qoi) - for k in Map.keys(): - - #the current code samples - q_k = samples[Map[k]['f']] - - #the hierarchical surplus (s_k) between the code output q_k and the - #previous surrogate of level L-1 evaluated at the same location. - #Recursively computed. - - if k in self.surr_lm1[goal][level]: - #print('surrogate already computed') - surr_lm1 = self.surr_lm1[goal][level][k] - else: - surr_lm1 = self.sc_expansion(level-1, goal, samples, x = Map[k]['X']) - self.surr_lm1[goal][level][k] = surr_lm1 - - s_k = q_k - surr_lm1 - - #the current level multi index (l_1,...,l_N) - l = Map[k]['l'] + weight = [] + for n in range(self.N): + #interpolate + x = kwargs['x'] + # add values of Lagrange polynomials at x + weight.append(lagrange_poly(x[n], self.xi_1d[n][l[n]], idx[n])) - idx = {} - # indices of current collocation point (Map[k]['X'][n]), - # in corresponding 1d colloc points (self.xi_1d[n][l[n]]) - # These are the j of the 1D lagrange polynomials l_j(x), see - # lagrange_poly subroutine - for n in range(self.N): - idx[n] = (self.xi_1d[n][l[n]] == Map[k]['X'][n]).nonzero()[0][0] - - weight = [] - for n in range(self.N): - #interpolate - if goal == 'interpolate': - x = args['x'] - # add values of Lagrange polynomials at x - weight.append(lagrange_poly(x[n], self.xi_1d[n][l[n]], idx[n])) - #integrate - elif goal == 'quadrature': - #add quadrature weights - weight.append(self.wi_1d[n][l[n]][idx[n]]) - - #Delta is the interpolation of the hierarchical surplus - Delta += s_k*np.prod(weight) + #Delta is the interpolation of the hierarchical surplus + surr += s_k*np.prod(weight) - surr += Delta - return surr - - def clear_surr_lm1(self, ID): + + def clear_surr_lm1(self): """ - Clears the interpolation or quadrature results in surr_lm1[ID]. + Clears the interpolation results in surr_lm1[ID]. surr_lm1 is a dictionary used to store surrogate results at previous level (l-1). Used to avoid recomputing the surrogate in the recursive sc_expansion subroutine. - surr_lm1['interpolation'][l][k] stores the interpolation results + surr_lm1[l][k] stores the interpolation results of level l-1 at collocation point X_k Parameters @@ -402,9 +441,9 @@ def clear_surr_lm1(self, ID): - ID (str): either 'interpolate' or 'quadrature' """ - self.surr_lm1[ID] = {} + self.surr_lm1 = {} for level in range(self.L_min, self.L+1): - self.surr_lm1[ID][level] = {} + self.surr_lm1[level] = {} def compute_SC_weights(self, rule): """ @@ -510,16 +549,19 @@ def plot_grid(self): @staticmethod def compute_tensor_prod_u(xi, wi, u, u_prime): - """Calculate tensor products with dimension of u + """ + Calculate tensor products of weights and collocation points + with dimension of u and u' Parameters ---------- - xi - wi - u - u_prime + xi (array of floats): 1D colloction points + wi (array of floats): 1D quadrature weights + u (array of int): dimensions + u_prime (array of int): remaining dimensions (u union u' = range(N)) Returns + dict of tensor products of weight and points for dimensions u and u' ------- """ @@ -547,27 +589,39 @@ def compute_tensor_prod_u(xi, wi, u, u_prime): return {'xi_d_u': xi_d_u, 'wi_d_u': wi_d_u, 'xi_d_u_prime': xi_d_u_prime, 'wi_d_u_prime': wi_d_u_prime} - def compute_h(self, qoi, u, u_prime, xi_d_u, xi_d_u_prime, wi_d_u_prime): + def compute_marginal(self, qoi, u, u_prime, diff): """ + Computes a marginal integral of the qoi(x) over the dimension defined + by u_prime, for every x value in dimensions u Parameters ---------- - qoi - u - u_prime - xi_d_u - xi_d_u_prime - wi_d_u_prime - + - qoi (str): name of the quantity of interest + - u (array of int): dimensions which are not integrated + - u_prime (array of int): dimensions which are integrated + - diff (array of int): levels + Returns + - Values of the marginal integral ------- """ + + #1d weights and points of the levels in diff + xi = [self.xi_1d[n][np.abs(diff)[n]] for n in range(self.N)] + wi = [self.wi_1d[n][np.abs(diff)[n]] for n in range(self.N)] + + # compute tensor products and weights in dimension u and u' + tmp = self.compute_tensor_prod_u(xi, wi, u, u_prime) + xi_d_u = tmp['xi_d_u'] + wi_d_u = tmp['wi_d_u'] + xi_d_u_prime = tmp['xi_d_u_prime'] + wi_d_u_prime = tmp['wi_d_u_prime'] S_u = xi_d_u.shape[0] S_u_prime = xi_d_u_prime.shape[0] - # coefficients h = f*w' integrated over u', so cardinality is that of u + # marginals h = f*w' integrated over u', so cardinality is that of u h = {} for i_u in range(S_u): h[i_u] = 0. @@ -588,26 +642,35 @@ def compute_h(self, qoi, u, u_prime, xi_d_u, xi_d_u_prime, wi_d_u_prime): xi_s[k] = xi_d_u_prime[i_up][idx] idx += 1 - # + #find the index of the corresponding code sample tmp = np.prod(self.xi_d == xi_s, axis=1) idx = np.where(tmp == 1)[0][0] - h[i_u] += self.samples[qoi][idx].values.flatten() * \ - wi_d_u_prime[i_up].prod() + + #perform quadrature + q_k = self.samples[qoi][idx].flatten() + h[i_u] += q_k * wi_d_u_prime[i_up].prod() - return h + #return marginal and the weights of dimensions u + return h, wi_d_u - def get_sobol_indices(self, qoi, typ): - """Computes Sobol indices using Stochastic Collocation + def get_sobol_indices(self, qoi, typ = 'first_order'): + """ + Computes Sobol indices using Stochastic Collocation. Method: + Tang (2009), GLOBAL SENSITIVITY ANALYSIS FOR STOCHASTIC COLLOCATION + EXPANSION. Parameters ---------- - qoi - typ + qoi (str): name of the Quantity of Interest for which to compute the indices + typ (str): Default = 'first_order'. 'all' is also possible Returns ------- + Either the first order or all Sobol indices of qoi """ + + print('Computing', typ, 'Sobol indices') # multi indices U = list(range(self.N)) @@ -623,13 +686,6 @@ def get_sobol_indices(self, qoi, typ): mu = mu.flatten() D = D.flatten() - # list of 1D nodes and quad weights - xi = self.xi - wi = self.wi - - # total variance might be zero at some locations, Sobol index not defined there - # idx_gt0 = np.where(D > 0)[0] - # partial variances D_u = {P[0]: mu**2} @@ -639,26 +695,25 @@ def get_sobol_indices(self, qoi, typ): # complement of u u_prime = np.delete(U, u) - - # compute corresponding tensor products and GQ weights - tmp = self.compute_tensor_prod_u(xi, wi, u, u_prime) - xi_d_u = tmp['xi_d_u'] - wi_d_u = tmp['wi_d_u'] - xi_d_u_prime = tmp['xi_d_u_prime'] - wi_d_u_prime = tmp['wi_d_u_prime'] - - # cardinality of u - S_u = xi_d_u.shape[0] - - # h coefficients - h = self.compute_h(qoi, u, u_prime, xi_d_u, xi_d_u_prime, wi_d_u_prime) - - # partial variance D_u[u] = 0.0 - for i_u in range(S_u): - D_u[u] += h[i_u]**2 * wi_d_u[i_u].prod() - D_u[u] = D_u[u].flatten() + + for l in self.l_norm: + + diff_idx = np.array(list(product(*[[k, -(k-1)] for k in l]))) + + for diff in diff_idx: + + if not np.min(np.abs(diff)) < self.L_min: + + #mariginal integral h, integrate over dimensions u' + h, wi_d_u = self.compute_marginal(qoi, u, u_prime, diff) + + #square result and integrate over remaining dimensions u + for i_u in range(wi_d_u.shape[0]): + D_u[u] += np.sign(np.prod(diff))*h[i_u]**2 * wi_d_u[i_u].prod() + D_u[u] = D_u[u].flatten() + # all subsets of u W = list(powerset(u))[0:-1] @@ -675,6 +730,7 @@ def get_sobol_indices(self, qoi, typ): sort.append(sobol[u]) + print('done') return sobol # End SC specific methods diff --git a/easyvvuq/sampling/stochastic_collocation.py b/easyvvuq/sampling/stochastic_collocation.py index 6023c892d..5c22f64d9 100644 --- a/easyvvuq/sampling/stochastic_collocation.py +++ b/easyvvuq/sampling/stochastic_collocation.py @@ -95,15 +95,26 @@ def __init__(self, self.xi_1d[n] = {} self.wi_1d[n] = {} - for n in range(N): - for i in range(1, self.quad_order+1): - xi_i, wi_i = cp.generate_quadrature(i, + if sparse == True: + for n in range(N): + for i in range(1, self.quad_order+1): + xi_i, wi_i = cp.generate_quadrature(i+1, + params_distribution[n], + rule=self.quad_rule, + growth=self.growth, + normalize=True) + + self.xi_1d[n][i] = xi_i[0] + self.wi_1d[n][i] = wi_i + else: + for n in range(N): + xi_i, wi_i = cp.generate_quadrature(self.quad_order, params_distribution[n], rule=self.quad_rule, growth=self.growth, - normalize=True) - self.xi_1d[n][i] = xi_i[0] - self.wi_1d[n][i] = wi_i + normalize=True) + self.xi_1d[n][self.quad_order] = xi_i[0] + self.wi_1d[n][self.quad_order] = wi_i if sparse == False: # the nodes of the collocation grid @@ -119,13 +130,16 @@ def __init__(self, if L < N: print("*************************************************************") print("Level of sparse grid is lower than the dimension N (# params)") - print("Increase level (quad_order) q such that q >= N") + print("Increase level (via polynomial_order) p such that p-1 >= N") print("*************************************************************") import sys; sys.exit() + + #multi-index l, such that |l| <= L + l_norm_le_L = self.compute_sparse_multi_idx(L, N) #create sparse grid of dimension N and level q using the 1d #rules in self.xi_1d - self.xi_d = self.sparse_grid(L, N) + self.xi_d = self.generate_grid(L, N, l_norm_le_L) self.L = L self.N = N @@ -180,20 +194,20 @@ def get_restart_dict(self): ======================= """ - def sparse_grid(self, L, N): - """ - Compute an isotropic sparse grid H_L_N of N dimensions and level q - """ - #multi-index l, such that |l| <= L - l_norm_le_L = self.compute_sparse_multi_idx(L, N) + def generate_grid(self, L, N, l_norm, **kwargs): + if 'dimensions' in kwargs: + dimensions = kwargs['dimensions'] + else: + dimensions = range(N) + H_L_N = [] #loop over all multi indices i - for i in l_norm_le_L: + for l in l_norm: #compute the tensor product of nodes indexed by i - X_i = [self.xi_1d[n][i[n]] for n in range(N)] - H_L_N.append(list(product(*X_i))) + X_l = [self.xi_1d[n][l[n]] for n in dimensions] + H_L_N.append(list(product(*X_l))) #flatten the list of lists H_L_N = np.array(list(chain(*H_L_N))) @@ -206,6 +220,6 @@ def compute_sparse_multi_idx(self, q, N): computes all N dimensional multi-indices i = (i1,...,iN) such that |i| <= Q. Here |i| is the internal sum of i (i1+...+iN) """ - P = np.array(list(product(range(1,q+1), repeat=N))) + P = np.array(list(product(range(1, q+1), repeat=N))) i_norm_le_q = P[np.where(np.sum(P, axis=1) <= q)[0]] return i_norm_le_q \ No newline at end of file From 73bde27ee62520ff904c179376ed66f7a8afd9e4 Mon Sep 17 00:00:00 2001 From: wedeling Date: Tue, 27 Aug 2019 23:24:45 +0200 Subject: [PATCH 18/48] fixed some comments --- easyvvuq/analysis/sc_analysis.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/easyvvuq/analysis/sc_analysis.py b/easyvvuq/analysis/sc_analysis.py index 2d88456b9..aaf82111b 100644 --- a/easyvvuq/analysis/sc_analysis.py +++ b/easyvvuq/analysis/sc_analysis.py @@ -699,8 +699,11 @@ def get_sobol_indices(self, qoi, typ = 'first_order'): for l in self.l_norm: + #expand the multi-index indices of the tensor product + #(Q^1_{i1} - Q^1_{i1-1}) X ... X (Q^1_{id) - Q^1_{id-1}) diff_idx = np.array(list(product(*[[k, -(k-1)] for k in l]))) + #perform analysis on each Q^1_l1 X ... X Q^1_l_N tensor prod for diff in diff_idx: if not np.min(np.abs(diff)) < self.L_min: From e655ed7ef7b29e23d69f13873773be95d598d2b5 Mon Sep 17 00:00:00 2001 From: orbitfold Date: Wed, 28 Aug 2019 12:23:15 +0200 Subject: [PATCH 19/48] Extra tests --- tests/test_db.py | 86 ++++++++++++++++++++++++++++++++---------------- 1 file changed, 57 insertions(+), 29 deletions(-) diff --git a/tests/test_db.py b/tests/test_db.py index 4e355e9ab..9a9c04429 100644 --- a/tests/test_db.py +++ b/tests/test_db.py @@ -8,7 +8,40 @@ @pytest.fixture -def campaign(tmp_path): +def app_info(): + app_info = AppInfo('test', uq.ParamsSpecification({ + "temp_init": { + "type": "float", + "min": 0.0, + "max": 100.0, + "default": 95.0}, + "kappa": { + "type": "float", + "min": 0.0, + "max": 0.1, + "default": 0.025}, + "t_env": { + "type": "float", + "min": 0.0, + "max": 40.0, + "default": 15.0}, + "out_file": { + "type": "string", + "default": "output.csv"} + }), + None, + uq.encoders.GenericEncoder( + template_fname='tests/cooling/cooling.template', + delimiter='$', + target_filename='cooling_in.json'), + uq.decoders.SimpleCSV(target_filename='output.csv', + output_columns=["te", "ti"], + header=0)) + return app_info + + +@pytest.fixture +def campaign(tmp_path, app_info): info = CampaignInfo( name='test', campaign_dir_prefix=default_campaign_prefix, @@ -19,34 +52,6 @@ def campaign(tmp_path): runs = [RunInfo('run', 'test', '.', 1, {'a' : 1}, 1, 1) for _ in range(1010)] run_names = ['Run_{}'.format(i) for i in range(1, 1011)] campaign.add_runs(runs) - app_info = AppInfo('test', uq.ParamsSpecification({ - "temp_init": { - "type": "float", - "min": 0.0, - "max": 100.0, - "default": 95.0}, - "kappa": { - "type": "float", - "min": 0.0, - "max": 0.1, - "default": 0.025}, - "t_env": { - "type": "float", - "min": 0.0, - "max": 40.0, - "default": 15.0}, - "out_file": { - "type": "string", - "default": "output.csv"} - }), - None, - uq.encoders.GenericEncoder( - template_fname='tests/cooling/cooling.template', - delimiter='$', - target_filename='cooling_in.json'), - uq.decoders.SimpleCSV(target_filename='output.csv', - output_columns=["te", "ti"], - header=0)) campaign.add_app(app_info) return campaign @@ -67,5 +72,28 @@ def test_get_num_runs(campaign): def test_app(campaign): + with pytest.raises(RuntimeError): + campaign.app('test_') app_dict = campaign.app('test') assert(app_dict['name'] == 'test') + assert(isinstance(app_dict, dict)) + + +def test_add_app(campaign, app_info): + with pytest.raises(RuntimeError): + campaign.add_app(app_info) + + +def test_campaign(campaign): + assert('test' in campaign.campaigns()) + + +def test_get_campaign_id(campaign): + with pytest.raises(RuntimeError): + campaign.get_campaign_id('test_') + assert(campaign.get_campaign_id('test') == 1) + + +def test_campaign_dir(campaign): + assert(campaign.campaign_dir('test') == default_campaign_prefix) + From e17e106cd663c73fe9b20c52f0534e8e7f3f70e9 Mon Sep 17 00:00:00 2001 From: wedeling Date: Wed, 28 Aug 2019 13:06:40 +0200 Subject: [PATCH 20/48] added sparse grid test --- tests/sc/sc_model.py | 1 + tests/sc/sparse_sc_test.py | 169 +++++++++++++++++++++++++++++++++++++ 2 files changed, 170 insertions(+) create mode 100644 tests/sc/sparse_sc_test.py diff --git a/tests/sc/sc_model.py b/tests/sc/sc_model.py index 1785e089a..1dc0fd78e 100755 --- a/tests/sc/sc_model.py +++ b/tests/sc/sc_model.py @@ -114,6 +114,7 @@ def shape(x, XA1, XA2): # run FEM solver at current value of Pe and f u = solve(Pe, f, nel=300) +u[0] = 0.0 # output csv file header = 'u' diff --git a/tests/sc/sparse_sc_test.py b/tests/sc/sparse_sc_test.py new file mode 100644 index 000000000..35099b4a9 --- /dev/null +++ b/tests/sc/sparse_sc_test.py @@ -0,0 +1,169 @@ +""" +SCRIPT TO TEST THE SPARSE GRID INTERPOLATION, QUADRATURE AND SOBOL INDICES +FOR THE STOCHASTIC COLLOCATION METHOD +""" + +import chaospy as cp +import numpy as np +import easyvvuq as uq +import matplotlib.pyplot as plt +import os + +plt.close('all') + +# author: Wouter Edeling +__license__ = "LGPL" + +#home directory of user +home = os.path.expanduser('~') +HOME = os.path.abspath(os.path.dirname(__file__)) + +# Set up a fresh campaign called "sc" +my_campaign = uq.Campaign(name='sc', work_dir='/tmp') + +# Define parameter space +params = { + "Pe": { + "type": "float", + "min": 1.0, + "max": 2000.0, + "default": 100.0}, + "f": { + "type": "float", + "min": 0.0, + "max": 10.0, + "default": 1.0}, + "out_file": { + "type": "string", + "default": "output.csv"}} + +output_filename = params["out_file"]["default"] +output_columns = ["u"] + +# Create an encoder, decoder and collation element +encoder = uq.encoders.GenericEncoder( + template_fname = HOME + '/sc.template', + delimiter='$', + target_filename='ade_in.json') +decoder = uq.decoders.SimpleCSV(target_filename=output_filename, + output_columns=output_columns, + header=0) + +# Create a collation element for this campaign +collater = uq.collate.AggregateSamples(average=False) +my_campaign.set_collater(collater) + +# Add the SC app (automatically set as current app) +my_campaign.add_app(name="sc", + params=params, + encoder=encoder, + decoder=decoder) + +# Create the sampler +vary = { + "Pe": cp.Uniform(100.0, 200.0), + "f": cp.Uniform(0.9, 1.1) +} + +""" +SPARSE GRID PARAMETERS +---------------------- +- sparse = True: use a Smolyak sparse grid +- growth = True: use an exponential rule for the growth of the number + of 1D collocation points per level. Used to make e.g. clenshaw-curtis + quadrature nested. +""" +my_sampler = uq.sampling.SCSampler(vary=vary, polynomial_order = 4, + quadrature_rule="C", sparse=True, + growth=True) + +# Associate the sampler with the campaign +my_campaign.set_sampler(my_sampler) + +# Will draw all (of the finite set of samples) +my_campaign.draw_samples() +my_campaign.populate_runs_dir() + +# Use this instead to run the samples using EasyVVUQ on the localhost +my_campaign.apply_for_each_run_dir(uq.actions.ExecuteLocal( + "./sc_model.py ade_in.json")) + +my_campaign.collate() + +# Post-processing analysis +analysis = uq.analysis.SCAnalysis(sampler=my_sampler, qoi_cols=output_columns) + +my_campaign.apply_analysis(analysis) + +results = my_campaign.get_last_analysis() + +################################### +# Plot the moments and SC samples # +################################### + +mu = results['statistical_moments']['u']['mean'] +std = results['statistical_moments']['u']['std'] + +x = np.linspace(0, 1, 301) + +fig = plt.figure(figsize=[10, 5]) +ax = fig.add_subplot(121, xlabel='x', ylabel='u', + title=r'code mean +/- standard deviation') +ax.plot(x, mu, 'b', label='mean') +ax.plot(x, mu + std, '--r', label='std-dev') +ax.plot(x, mu - std, '--r') + +##################################### +# Plot the random surrogate samples # +##################################### + +if analysis.element_name() == 'SC_Analysis': + + ax = fig.add_subplot(122, xlabel='x', ylabel='u', + title='Surrogate and code samples') + + xi_mc = analysis.xi_d + n_mc = xi_mc.shape[0] + code_samples = analysis.get_sample_array('u') + + # evaluate the surrogate at these values + print('Evaluating surrogate model', n_mc, 'times') + for i in range(n_mc): + ax.plot(x, analysis.surrogate('u', xi_mc[i]), 'g') + ax.plot(x, code_samples[i], 'r+') + print('done') + + plt.tight_layout() + + analysis.plot_grid() + # + ####################### + ## Plot Sobol indices # + ####################### + + fig = plt.figure() + ax = fig.add_subplot( + 111, + xlabel='x', + ylabel='Sobol indices', + title='spatial dist. Sobol indices, Pe only important in viscous regions') + + lbl = ['Pe', 'f', 'Pe-f interaction'] + idx = 0 + + if analysis.element_name() == 'SC_Analysis': + + for S_i in results['sobol_indices']['u']: + ax.plot(x, results['sobol_indices']['u'][S_i], label=lbl[idx]) + idx += 1 + else: + for S_i in results['sobol_indices']['u'][1]: + ax.plot(x, results['sobol_indices']['u'][1][S_i], label=lbl[idx]) + idx += 1 + + leg = plt.legend(loc=0) + leg.set_draggable(True) + + plt.tight_layout() + +plt.show() \ No newline at end of file From 0cbea8c166a95fc0c3a3ef7bb586952e321c75f4 Mon Sep 17 00:00:00 2001 From: wedeling Date: Wed, 28 Aug 2019 13:34:42 +0200 Subject: [PATCH 21/48] cleanup --- easyvvuq/analysis/sc_analysis.py | 3 +++ easyvvuq/sampling/stochastic_collocation.py | 12 ++++++------ 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/easyvvuq/analysis/sc_analysis.py b/easyvvuq/analysis/sc_analysis.py index aaf82111b..34646f717 100644 --- a/easyvvuq/analysis/sc_analysis.py +++ b/easyvvuq/analysis/sc_analysis.py @@ -111,7 +111,10 @@ def analyse(self, data_frame=None): self.L_min = 1 self.l_norm = self.sampler.compute_sparse_multi_idx(self.L, self.N) + #full tensor grid self.xi_d = self.sampler.xi_d + + #1d weights and points per level self.xi_1d = self.sampler.xi_1d self.wi_1d = self.compute_SC_weights(rule=self.sampler.quad_rule) diff --git a/easyvvuq/sampling/stochastic_collocation.py b/easyvvuq/sampling/stochastic_collocation.py index 5c22f64d9..784ade25c 100644 --- a/easyvvuq/sampling/stochastic_collocation.py +++ b/easyvvuq/sampling/stochastic_collocation.py @@ -215,11 +215,11 @@ def generate_grid(self, L, N, l_norm, **kwargs): #return unique nodes return np.unique(H_L_N, axis=0) - def compute_sparse_multi_idx(self, q, N): + def compute_sparse_multi_idx(self, L, N): """ - computes all N dimensional multi-indices i = (i1,...,iN) such that - |i| <= Q. Here |i| is the internal sum of i (i1+...+iN) + computes all N dimensional multi-indices l = (l1,...,lN) such that + |l| <= Q. Here |l| is the internal sum of i (l1+...+lN) """ - P = np.array(list(product(range(1, q+1), repeat=N))) - i_norm_le_q = P[np.where(np.sum(P, axis=1) <= q)[0]] - return i_norm_le_q \ No newline at end of file + P = np.array(list(product(range(1, L+1), repeat=N))) + l_norm_le_q = P[np.where(np.sum(P, axis=1) <= L)[0]] + return l_norm_le_q \ No newline at end of file From f474c77fdb923b18557f2df4d7dbb0e5f62b9bc6 Mon Sep 17 00:00:00 2001 From: orbitfold Date: Thu, 29 Aug 2019 16:14:51 +0200 Subject: [PATCH 22/48] fix to _get_campaign_info in sql.py --- easyvvuq/db/sql.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/easyvvuq/db/sql.py b/easyvvuq/db/sql.py index 8e6e4c681..06038f9ae 100644 --- a/easyvvuq/db/sql.py +++ b/easyvvuq/db/sql.py @@ -547,7 +547,7 @@ def _get_campaign_info(self, campaign_name=None): else: campaign_info = query.filter_by(name=campaign_name).all() - if campaign_info.count() > 1: + if len(campaign_info) > 1: logger.warning( 'More than one campaign selected - using first one.') elif campaign_info.count() == 0: From 67e8869b8a8886c739fd568a9ce0299799bce09f Mon Sep 17 00:00:00 2001 From: orbitfold Date: Thu, 29 Aug 2019 16:46:41 +0200 Subject: [PATCH 23/48] fix to _get_campaign_info in sql.py (for real this time) --- easyvvuq/db/sql.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/easyvvuq/db/sql.py b/easyvvuq/db/sql.py index 06038f9ae..10d29ad68 100644 --- a/easyvvuq/db/sql.py +++ b/easyvvuq/db/sql.py @@ -547,13 +547,14 @@ def _get_campaign_info(self, campaign_name=None): else: campaign_info = query.filter_by(name=campaign_name).all() - if len(campaign_info) > 1: - logger.warning( - 'More than one campaign selected - using first one.') - elif campaign_info.count() == 0: - message = 'No campaign available.' - logger.critical(message) - raise RuntimeError(message) + if campaign_name is not None: + if len(campaign_info) > 1: + logger.warning( + 'More than one campaign selected - using first one.') + elif len(campaign_info) == 0: + message = 'No campaign available.' + logger.critical(message) + raise RuntimeError(message) return campaign_info.first() From b0f4e6d03651c09f6d9c988d2be524372bcacd0a Mon Sep 17 00:00:00 2001 From: David W Wright Date: Mon, 2 Sep 2019 15:42:17 +0100 Subject: [PATCH 24/48] Remove matplotlib requirement which is no longer used. --- requirements.txt | 1 - setup.py | 4 ++-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/requirements.txt b/requirements.txt index a82834046..c230984bb 100644 --- a/requirements.txt +++ b/requirements.txt @@ -9,4 +9,3 @@ SQLAlchemy sqlalchemy-utils jsonpickle cerberus -matplotlib diff --git a/setup.py b/setup.py index a6c313ffb..efaad70a5 100644 --- a/setup.py +++ b/setup.py @@ -27,8 +27,8 @@ install_requires=['numpy', 'pandas>=0.24', 'scipy', 'pytest', 'SQLAlchemy', 'chaospy', - 'sqlalchemy-utils', 'matplotlib', - 'jsonpickle', 'cerberus', 'SALib'], + 'sqlalchemy-utils', 'jsonpickle', + 'cerberus', 'SALib'], packages=find_packages(), From 50d581ba5701774beebea3b2109fd65f89335a90 Mon Sep 17 00:00:00 2001 From: David W Wright Date: Mon, 2 Sep 2019 15:50:21 +0100 Subject: [PATCH 25/48] Fix whitespace --- setup.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/setup.py b/setup.py index efaad70a5..ea6359c47 100644 --- a/setup.py +++ b/setup.py @@ -2,8 +2,6 @@ from setuptools import setup, find_packages import versioneer - - # read the contents of README file this_directory = path.abspath(path.dirname(__file__)) with open(path.join(this_directory, 'README.md'), encoding='utf-8') as f: From fbaaa148b180345eabf0dfe1c996f8bc213edb59 Mon Sep 17 00:00:00 2001 From: David W Wright Date: Mon, 2 Sep 2019 16:00:45 +0100 Subject: [PATCH 26/48] Contributions file. --- CONTRIBUTIONS.md | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 CONTRIBUTIONS.md diff --git a/CONTRIBUTIONS.md b/CONTRIBUTIONS.md new file mode 100644 index 000000000..6b5e71063 --- /dev/null +++ b/CONTRIBUTIONS.md @@ -0,0 +1,18 @@ +# Credits + +In case you wanted to know who to thank/blame for different features ;) +ORCiD's provided for specificity. + +## Concepts and core code development + +* David W. Wright ([0000-0002-5124-8044](https://orcid.org/0000-0002-5124-8044)) +* Robin A. Richardson + +## Statistical methods implementation + +* Quasi Monte Carlo and Polynomial Chaos Expansion - Jalal Lakhlili +* Stochastic Collocation - Wouter Edeling + +## Testing + +* Vytautas Jancauskas From 3a2b61a62f0a0d2b17b1731069de953b9e8b43f3 Mon Sep 17 00:00:00 2001 From: di73kuj2 Date: Wed, 4 Sep 2019 14:57:43 +0200 Subject: [PATCH 27/48] Fix to the test_campaign_dir test --- easyvvuq/db/sql.py | 3 ++- tests/test_db.py | 6 +++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/easyvvuq/db/sql.py b/easyvvuq/db/sql.py index b6ba7246e..53dc23c44 100644 --- a/easyvvuq/db/sql.py +++ b/easyvvuq/db/sql.py @@ -555,7 +555,8 @@ def _get_campaign_info(self, campaign_name=None): message = 'No campaign available.' logger.critical(message) raise RuntimeError(message) - + return campaign_info[0] + return campaign_info.first() def get_campaign_id(self, name): diff --git a/tests/test_db.py b/tests/test_db.py index 9a9c04429..76bb03253 100644 --- a/tests/test_db.py +++ b/tests/test_db.py @@ -46,9 +46,9 @@ def campaign(tmp_path, app_info): name='test', campaign_dir_prefix=default_campaign_prefix, easyvvuq_version=uq.__version__, - campaign_dir='.') + campaign_dir=str(tmp_path)) campaign = CampaignDB(location='sqlite:///{}/test.sqlite'.format(tmp_path), new_campaign=True, name='test', info=info) - campaign.tmp_path = tmp_path + campaign.tmp_path = str(tmp_path) runs = [RunInfo('run', 'test', '.', 1, {'a' : 1}, 1, 1) for _ in range(1010)] run_names = ['Run_{}'.format(i) for i in range(1, 1011)] campaign.add_runs(runs) @@ -95,5 +95,5 @@ def test_get_campaign_id(campaign): def test_campaign_dir(campaign): - assert(campaign.campaign_dir('test') == default_campaign_prefix) + assert(campaign.campaign_dir('test') == campaign.tmp_path) From 45d3d31de7c9f14dbf145fe7d0500eb3a6f8ed99 Mon Sep 17 00:00:00 2001 From: di73kuj2 Date: Wed, 4 Sep 2019 15:19:40 +0200 Subject: [PATCH 28/48] pep8ing it up --- tests/test_db.py | 64 ++++++++++++++++++++++++------------------------ 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/tests/test_db.py b/tests/test_db.py index 76bb03253..77f07f576 100644 --- a/tests/test_db.py +++ b/tests/test_db.py @@ -10,33 +10,33 @@ @pytest.fixture def app_info(): app_info = AppInfo('test', uq.ParamsSpecification({ - "temp_init": { - "type": "float", - "min": 0.0, - "max": 100.0, - "default": 95.0}, - "kappa": { - "type": "float", - "min": 0.0, - "max": 0.1, - "default": 0.025}, - "t_env": { - "type": "float", - "min": 0.0, - "max": 40.0, - "default": 15.0}, - "out_file": { - "type": "string", - "default": "output.csv"} - }), - None, - uq.encoders.GenericEncoder( - template_fname='tests/cooling/cooling.template', - delimiter='$', - target_filename='cooling_in.json'), - uq.decoders.SimpleCSV(target_filename='output.csv', - output_columns=["te", "ti"], - header=0)) + "temp_init": { + "type": "float", + "min": 0.0, + "max": 100.0, + "default": 95.0}, + "kappa": { + "type": "float", + "min": 0.0, + "max": 0.1, + "default": 0.025}, + "t_env": { + "type": "float", + "min": 0.0, + "max": 40.0, + "default": 15.0}, + "out_file": { + "type": "string", + "default": "output.csv"}}), + None, + uq.encoders.GenericEncoder( + template_fname='tests/cooling/cooling.template', + delimiter='$', + target_filename='cooling_in.json'), + uq.decoders.SimpleCSV( + target_filename='output.csv', + output_columns=["te", "ti"], + header=0)) return app_info @@ -47,9 +47,10 @@ def campaign(tmp_path, app_info): campaign_dir_prefix=default_campaign_prefix, easyvvuq_version=uq.__version__, campaign_dir=str(tmp_path)) - campaign = CampaignDB(location='sqlite:///{}/test.sqlite'.format(tmp_path), new_campaign=True, name='test', info=info) + campaign = CampaignDB(location='sqlite:///{}/test.sqlite'.format(tmp_path), + new_campaign=True, name='test', info=info) campaign.tmp_path = str(tmp_path) - runs = [RunInfo('run', 'test', '.', 1, {'a' : 1}, 1, 1) for _ in range(1010)] + runs = [RunInfo('run', 'test', '.', 1, {'a': 1}, 1, 1) for _ in range(1010)] run_names = ['Run_{}'.format(i) for i in range(1, 1011)] campaign.add_runs(runs) campaign.add_app(app_info) @@ -78,7 +79,7 @@ def test_app(campaign): assert(app_dict['name'] == 'test') assert(isinstance(app_dict, dict)) - + def test_add_app(campaign, app_info): with pytest.raises(RuntimeError): campaign.add_app(app_info) @@ -93,7 +94,6 @@ def test_get_campaign_id(campaign): campaign.get_campaign_id('test_') assert(campaign.get_campaign_id('test') == 1) - + def test_campaign_dir(campaign): assert(campaign.campaign_dir('test') == campaign.tmp_path) - From bb2e463cdd8d5ee21c1c2fc8f59dc5de876eba3c Mon Sep 17 00:00:00 2001 From: Robin Richardson Date: Fri, 6 Sep 2019 15:37:17 +0100 Subject: [PATCH 29/48] Make uniform_distribution return python int instead of numpy.int64 (which was breaking cerberus and serialization) --- easyvvuq/distributions/distributions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/easyvvuq/distributions/distributions.py b/easyvvuq/distributions/distributions.py index a4856ea19..83d10eaba 100644 --- a/easyvvuq/distributions/distributions.py +++ b/easyvvuq/distributions/distributions.py @@ -97,7 +97,7 @@ def sample(self, size=()): Random samples with shape ``size``. """ - return np.random.randint(self.lo, self.up, size) + return [int(i) for i in np.random.randint(self.lo, self.up, size)] # TODO: Convert this to a chaospy style distribution From a22edb99d69509aaa9effe5e80d44264ef15ff0f Mon Sep 17 00:00:00 2001 From: di73kuj2 Date: Mon, 9 Sep 2019 14:38:13 +0200 Subject: [PATCH 30/48] Fix get_campaign_id --- easyvvuq/db/sql.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/easyvvuq/db/sql.py b/easyvvuq/db/sql.py index 53dc23c44..ff1c2889a 100644 --- a/easyvvuq/db/sql.py +++ b/easyvvuq/db/sql.py @@ -576,18 +576,18 @@ def get_campaign_id(self, name): selected = self.session.query( CampaignTable.name.label(name), - CampaignTable.id).all() + CampaignTable.id).filter(CampaignTable.name == name).all() if len(selected) == 0: msg = f"No campaign with name {name} found in campaign database" logger.error(msg) - raise Exception(msg) + raise RuntimeError(msg) if len(selected) > 1: msg = ( f"More than one campaign with name {name} found in" f"campaign database. Database state is compromised." ) logger.error(msg) - raise Exception(msg) + raise RuntimeError(msg) # Return the database ID for the specified campaign return selected[0][1] From 002a5525b3d4bfbf26557deef21f16032a2ae5e9 Mon Sep 17 00:00:00 2001 From: di73kuj2 Date: Mon, 9 Sep 2019 14:42:11 +0200 Subject: [PATCH 31/48] pep8fix --- easyvvuq/db/sql.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/easyvvuq/db/sql.py b/easyvvuq/db/sql.py index ff1c2889a..a9c3af277 100644 --- a/easyvvuq/db/sql.py +++ b/easyvvuq/db/sql.py @@ -556,7 +556,7 @@ def _get_campaign_info(self, campaign_name=None): logger.critical(message) raise RuntimeError(message) return campaign_info[0] - + return campaign_info.first() def get_campaign_id(self, name): From 826e3565e7fae22bd5a8d1ef8e70fc4d4102ce3f Mon Sep 17 00:00:00 2001 From: di73kuj2 Date: Tue, 10 Sep 2019 15:32:48 +0200 Subject: [PATCH 32/48] Initial commit --- tests/test_restart.py | 61 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 tests/test_restart.py diff --git a/tests/test_restart.py b/tests/test_restart.py new file mode 100644 index 000000000..24da395a3 --- /dev/null +++ b/tests/test_restart.py @@ -0,0 +1,61 @@ +import pytest +import easyvvuq as uq +from gauss.decoder_gauss import GaussDecoder + + +@pytest.fixture +def test_restart(tmpdir): + my_campaign = uq.Campaign(name=campaign_name, work_dir=tmpdir, db_type='sql') + params = { + "sigma": { + "type": "float", + "min": 0.0, + "max": 100000.0, + "default": 0.25 + }, + "mu": { + "type": "float", + "min": 0.0, + "max": 100000.0, + "default": 1 + }, + "num_steps": { + "type": "integer", + "min": 0, + "max": 100000, + "default": 10 + }, + "out_file": { + "type": "string", + "default": "output.csv" + }, + "bias": { + "type": "fixture", + "allowed": ["bias1", "bias2"], + "default": "bias1" + } + } + encoder = uq.encoders.GenericEncoder(template_fname='tests/gauss/gauss.template', + target_filename='gauss_in.json') + decoder = GaussDecoder(target_filename=params['out_file']['default']) + my_campaign.add_app(name='gauss', + params=params, + encoder=encoder, + decoder=decoder, + fixtures=None) + my_campaign.set_app('gauss') + collater = uq.collate.AggregateSamples(average=False) + my_campaign.set_collater(collater) + sampler = uq.sampling.RandomSampler(vary=vary) + my_campaign.set_sampler(sampler) + my_campaign.draw_samples(num_samples=2, replicas=2) + my_campaign.populate_runs_dir() + my_campaign.collate() + state_file = tmpdir + "{}_state.json".format('gauss') + my_campaign.save_state(state_file) + my_campaign = None + reloaded_campaign = uq.Campaign(state_file=state_file, work_dir=tmpdir) + reloaded_campaign.set_app('gauss') + reloaded_campaign.draw_samples(num_samples=2, replicas=2) + reloaded_campaign.populate_runs_dir() + From f42bbfbd9983d6afdde79a918e91ab97b0ae7ad0 Mon Sep 17 00:00:00 2001 From: di73kuj2 Date: Tue, 10 Sep 2019 16:40:53 +0200 Subject: [PATCH 33/48] Initial test --- tests/test_restart.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/tests/test_restart.py b/tests/test_restart.py index 24da395a3..827de9d4b 100644 --- a/tests/test_restart.py +++ b/tests/test_restart.py @@ -1,11 +1,12 @@ import pytest import easyvvuq as uq +import chaospy as cp from gauss.decoder_gauss import GaussDecoder @pytest.fixture -def test_restart(tmpdir): - my_campaign = uq.Campaign(name=campaign_name, work_dir=tmpdir, db_type='sql') +def restart(tmpdir): + my_campaign = uq.Campaign(name='gauss', work_dir=tmpdir, db_type='sql') params = { "sigma": { "type": "float", @@ -46,6 +47,9 @@ def test_restart(tmpdir): my_campaign.set_app('gauss') collater = uq.collate.AggregateSamples(average=False) my_campaign.set_collater(collater) + vary = { + "mu": cp.Uniform(1.0, 100.0), + } sampler = uq.sampling.RandomSampler(vary=vary) my_campaign.set_sampler(sampler) my_campaign.draw_samples(num_samples=2, replicas=2) @@ -58,4 +62,7 @@ def test_restart(tmpdir): reloaded_campaign.set_app('gauss') reloaded_campaign.draw_samples(num_samples=2, replicas=2) reloaded_campaign.populate_runs_dir() + return reloaded_campaign +def test_restart(restart): + assert(restart.campaign_db is not None) From c8b0d4021157f88e5d0bd719303ca127a44f805c Mon Sep 17 00:00:00 2001 From: di73kuj2 Date: Tue, 10 Sep 2019 16:43:45 +0200 Subject: [PATCH 34/48] 2 lines pep8 fix --- tests/test_restart.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_restart.py b/tests/test_restart.py index 827de9d4b..57ae6e57e 100644 --- a/tests/test_restart.py +++ b/tests/test_restart.py @@ -64,5 +64,6 @@ def restart(tmpdir): reloaded_campaign.populate_runs_dir() return reloaded_campaign + def test_restart(restart): assert(restart.campaign_db is not None) From 5c5e67da50f6e5b74d3379e2b51cb6d745850985 Mon Sep 17 00:00:00 2001 From: di73kuj2 Date: Wed, 11 Sep 2019 13:55:19 +0200 Subject: [PATCH 35/48] Added test_app --- tests/test_restart.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/tests/test_restart.py b/tests/test_restart.py index 57ae6e57e..34d08fff7 100644 --- a/tests/test_restart.py +++ b/tests/test_restart.py @@ -60,10 +60,17 @@ def restart(tmpdir): my_campaign = None reloaded_campaign = uq.Campaign(state_file=state_file, work_dir=tmpdir) reloaded_campaign.set_app('gauss') - reloaded_campaign.draw_samples(num_samples=2, replicas=2) - reloaded_campaign.populate_runs_dir() + #reloaded_campaign.draw_samples(num_samples=2, replicas=2) + #reloaded_campaign.populate_runs_dir() + reloaded_campaign.params_ = params return reloaded_campaign def test_restart(restart): assert(restart.campaign_db is not None) + + +def test_app(restart): + db = restart.campaign_db + assert(db.app('gauss')['name'] == 'gauss') + assert(db.app('gauss')['params'].params_dict['mu']['max'] == 100000.0) From ad2c626e3ee4343b5a02caa45fa5360eb62f637e Mon Sep 17 00:00:00 2001 From: di73kuj2 Date: Wed, 11 Sep 2019 14:12:46 +0200 Subject: [PATCH 36/48] Added test_runs --- tests/test_restart.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/tests/test_restart.py b/tests/test_restart.py index 34d08fff7..9097abc91 100644 --- a/tests/test_restart.py +++ b/tests/test_restart.py @@ -60,8 +60,6 @@ def restart(tmpdir): my_campaign = None reloaded_campaign = uq.Campaign(state_file=state_file, work_dir=tmpdir) reloaded_campaign.set_app('gauss') - #reloaded_campaign.draw_samples(num_samples=2, replicas=2) - #reloaded_campaign.populate_runs_dir() reloaded_campaign.params_ = params return reloaded_campaign @@ -74,3 +72,10 @@ def test_app(restart): db = restart.campaign_db assert(db.app('gauss')['name'] == 'gauss') assert(db.app('gauss')['params'].params_dict['mu']['max'] == 100000.0) + + +def test_runs(restart): + db = restart.campaign_db + assert(db.get_num_runs() == 4) + restart.draw_samples(num_samples=2, replicas=2) + assert(db.get_num_runs() == 8) From 2a6b1d7385cf6bccf3e8510cdb7788aa5a911dd1 Mon Sep 17 00:00:00 2001 From: Derek Groen Date: Wed, 11 Sep 2019 15:04:43 +0100 Subject: [PATCH 37/48] Fix Pandas CSV reader to read columns properly Fix Pandas CSV reader to read columns instead of replacing them --- easyvvuq/decoders/simple_csv.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/easyvvuq/decoders/simple_csv.py b/easyvvuq/decoders/simple_csv.py index fe3a6672c..fcb8e3a0a 100644 --- a/easyvvuq/decoders/simple_csv.py +++ b/easyvvuq/decoders/simple_csv.py @@ -87,7 +87,7 @@ def parse_sim_output(self, run_info={}): data = pd.read_csv( out_path, - names=self.output_columns, + use_cols=self.output_columns, header=self.header) return data From ac63d3eabd818c4ace0226cc04bb651947750c6a Mon Sep 17 00:00:00 2001 From: Derek Groen Date: Wed, 11 Sep 2019 15:09:24 +0100 Subject: [PATCH 38/48] Update simple_csv.py --- easyvvuq/decoders/simple_csv.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/easyvvuq/decoders/simple_csv.py b/easyvvuq/decoders/simple_csv.py index fcb8e3a0a..fef199bdb 100644 --- a/easyvvuq/decoders/simple_csv.py +++ b/easyvvuq/decoders/simple_csv.py @@ -87,7 +87,7 @@ def parse_sim_output(self, run_info={}): data = pd.read_csv( out_path, - use_cols=self.output_columns, + usecols=self.output_columns, header=self.header) return data From 0288cad82b5f040433fb0f610e54ec3a7523ba8c Mon Sep 17 00:00:00 2001 From: di73kuj2 Date: Wed, 11 Sep 2019 17:39:17 +0200 Subject: [PATCH 39/48] Added test_encoder and test_decoder --- tests/test_restart.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/tests/test_restart.py b/tests/test_restart.py index 9097abc91..8b05a6d29 100644 --- a/tests/test_restart.py +++ b/tests/test_restart.py @@ -79,3 +79,16 @@ def test_runs(restart): assert(db.get_num_runs() == 4) restart.draw_samples(num_samples=2, replicas=2) assert(db.get_num_runs() == 8) + + +def test_encoder(restart): + app = restart.campaign_db.app('gauss') + encoder = uq.encoders.GenericEncoder(template_fname='tests/gauss/gauss.template', + target_filename='gauss_in.json') + assert(app['input_encoder'] == encoder.serialize()) + + +def test_decoder(restart): + app = restart.campaign_db.app('gauss') + decoder = GaussDecoder(target_filename=restart.params_['out_file']['default']) + assert(app['output_decoder'] == decoder.serialize()) From ccc2c4b5a6113fad020e674bec1f09a13afe6d3a Mon Sep 17 00:00:00 2001 From: Robin Richardson Date: Fri, 13 Sep 2019 13:06:39 +0100 Subject: [PATCH 40/48] Fix failing pep8 check --- easyvvuq/analysis/sc_analysis.py | 422 ++++++++++---------- easyvvuq/sampling/stochastic_collocation.py | 77 ++-- tests/sc/sparse_sc_test.py | 34 +- 3 files changed, 267 insertions(+), 266 deletions(-) diff --git a/easyvvuq/analysis/sc_analysis.py b/easyvvuq/analysis/sc_analysis.py index 34646f717..c324c0ddc 100644 --- a/easyvvuq/analysis/sc_analysis.py +++ b/easyvvuq/analysis/sc_analysis.py @@ -95,38 +95,38 @@ def analyse(self, data_frame=None): raise RuntimeError( "No data in data frame passed to analyse element") - #the maximum level (quad order) of the (sparse) grid + # the maximum level (quad order) of the (sparse) grid self.L = self.sampler.L - - #the number of uncertain parameters + + # the number of uncertain parameters self.N = self.sampler.N - - #if L < L_min: quadratures and interpolations are zero - #For full tensor grid: there is only one level: L_min = L - if self.sparse == False: + + # if L < L_min: quadratures and interpolations are zero + # For full tensor grid: there is only one level: L_min = L + if not self.sparse: self.L_min = self.L self.l_norm = np.array([[self.L, self.L]]) - #For sparse grid: multiple levels, L >= N must hold - else: + # For sparse grid: multiple levels, L >= N must hold + else: self.L_min = 1 - self.l_norm = self.sampler.compute_sparse_multi_idx(self.L, self.N) + self.l_norm = self.sampler.compute_sparse_multi_idx(self.L, self.N) - #full tensor grid + # full tensor grid self.xi_d = self.sampler.xi_d - - #1d weights and points per level + + # 1d weights and points per level self.xi_1d = self.sampler.xi_1d self.wi_1d = self.compute_SC_weights(rule=self.sampler.quad_rule) - - #per level, map a unique index k to all (level multi indices, colloc points) - #combinations. Will differ for sparse or full tensor grids. - #All interpolation/quadrature subroutines loop over the entries in Map + + # per level, map a unique index k to all (level multi indices, colloc points) + # combinations. Will differ for sparse or full tensor grids. + # All interpolation/quadrature subroutines loop over the entries in Map self.Map = {} self.surr_lm1 = {} - + self.foo = [] - - for level in range(self.L_min, self.L+1): + + for level in range(self.L_min, self.L + 1): self.Map[level] = self.create_map(self.N, level) self.clear_surr_lm1() @@ -145,8 +145,8 @@ def analyse(self, data_frame=None): results = {'statistical_moments': {}, 'sobol_indices': {k: {} for k in self.qoi_cols}} - - #Compute descriptive statistics for each quantity of interest + + # Compute descriptive statistics for each quantity of interest for qoi_k in qoi_cols: mean_k, var_k = self.get_moments(qoi_k) std_k = var_k**0.5 @@ -159,60 +159,59 @@ def analyse(self, data_frame=None): results['sobol_indices'][qoi_k] = self.get_sobol_indices(qoi_k, 'all') return results - + def create_map(self, N, L): """ - Create a map from a unique integer k to each + Create a map from a unique integer k to each (level multi index l, collocation point X) combination. Also compute the index of X (f) in the global (sparse) grid xi_d - + Parameters ---------- - N (int) = number of parameters - L (int) = max level of grid - + Returns -------- - Map: a dict for level L containing k, l, X, and f """ - #unique index + # unique index k = 0 - Map = {} - + Map = {} + print('Creating multi-index map for level', L, '...') - - #full tensor product - if self.sparse == False: - - l = (np.ones(N)*L).astype('int') - + + # full tensor product + if not self.sparse: + + l = (np.ones(N) * L).astype('int') + for x in self.xi_d: - Map[k] = {'l':l, 'X':x, 'f':k} + Map[k] = {'l': l, 'X': x, 'f': k} k += 1 - #sparse grid + # sparse grid else: - #all sparse grid multi indices l with |l| <= L + # all sparse grid multi indices l with |l| <= L l_norm_le_L = self.sampler.compute_sparse_multi_idx(L, N) - - #loop over all multi indices + + # loop over all multi indices for l in l_norm_le_L: - - #colloc point of current level index l + + # colloc point of current level index l X_l = [self.xi_1d[n][l[n]] for n in range(N)] X_l = np.array(list(product(*X_l))) - - + for x in X_l: j = np.where((x == self.xi_d).all(axis=1))[0][0] - Map[k] = {'l':l, 'X':x, 'f':j} - k += 1; - + Map[k] = {'l': l, 'X': x, 'f': j} + k += 1 + print('done.') return Map - + def surrogate(self, qoi, x, **kwargs): """ Use sc_expansion UQP as a surrogate @@ -220,20 +219,20 @@ def surrogate(self, qoi, x, **kwargs): Parameters ---------- - qoi (str): name of the qoi - + Returns ------- the interpolated value of qoi at x (float, (N_qoi,)) - + """ - + if 'L' in kwargs: L = kwargs['L'] else: L = self.L - + return self.sc_expansion(L, self.samples[qoi], x=x) - + def quadrature(self, qoi, **kwargs): """ Computes a (Smolyak) quadrature @@ -241,60 +240,59 @@ def quadrature(self, qoi, **kwargs): Parameters ---------- - qoi (str): name of the qoi - + - samples (optional in kwargs): Default: compute the mean by setting samples = self.samples. To compute the variance, set samples = (self.samples - mean)**2 """ - - + if 'samples' in kwargs: samples = kwargs['samples'] else: samples = self.samples[qoi] - + Delta = np.zeros([self.l_norm.shape[0], self.N_qoi]) idx = 0 for l in self.l_norm: - #compute the Delta Q := - #(Q^1_l_1 - Q^1_{l_1 - 1}) X ... X (Q^1_{l_N} - Q^1_{L_N - 1}) - #tensor product + # compute the Delta Q := + # (Q^1_l_1 - Q^1_{l_1 - 1}) X ... X (Q^1_{l_N} - Q^1_{L_N - 1}) + # tensor product Delta[idx, :] = self.compute_Q_diff(l, samples) idx += 1 - + quadrature_approx = np.sum(Delta, axis=0) return quadrature_approx def compute_Q_diff(self, l, samples, **kwargs): - """ ======================================================================= - For every multi index l = (l1, l2, ..., ld), Smolyak sums over + For every multi index l = (l1, l2, ..., ld), Smolyak sums over tensor products difference quadrature rules: (Q^1_{l1} - Q^1_{l1-1}) X ... X (Q^1_{lN) - Q^1_{lN-1}) - Below this product is expanded into individual tensor products, each + Below this product is expanded into individual tensor products, each of which is then computed as: - Q^1_{k1} X ... X Q^1_{kN} = sum...sum w_{k1}*...*w{kN}*f(x_{k1},...,x_{kN}) + Q^1_{k1} X ... X Q^1_{kN} = sum...sum w_{k1}*...*w{kN}*f(x_{k1},...,x_{kN}) ======================================================================= """ - - #expand the multi-index indices of the tensor product - #(Q^1_{i1} - Q^1_{i1-1}) X ... X (Q^1_{id) - Q^1_{id-1}) - diff_idx = np.array(list(product(*[[k, -(k-1)] for k in l]))) - - #Delta will be the sum of all expanded tensor products - #Q^1_{k1} X ... X Q^1_{kd} = sum...sum w_{k1}*...*w{kN}*f(x_{k1},...,x_{kd}) + + # expand the multi-index indices of the tensor product + # (Q^1_{i1} - Q^1_{i1-1}) X ... X (Q^1_{id) - Q^1_{id-1}) + diff_idx = np.array(list(product(*[[k, -(k - 1)] for k in l]))) + + # Delta will be the sum of all expanded tensor products + # Q^1_{k1} X ... X Q^1_{kd} = sum...sum w_{k1}*...*w{kN}*f(x_{k1},...,x_{kd}) Delta = 0.0 - - #each diff contains the level indices of to a single - #Q^1_{k1} X ... X Q^1_{kN} product + + # each diff contains the level indices of to a single + # Q^1_{k1} X ... X Q^1_{kN} product for diff in diff_idx: - #Q_0 = 0, so if present do not compute Q^1_{k1} X ... X Q^1_{kd} + # Q_0 = 0, so if present do not compute Q^1_{k1} X ... X Q^1_{kd} if not np.min(np.abs(diff)) < self.L_min: - - #compute the tensor product of parameter and weight values - X_k = []; W_k = []; + + # compute the tensor product of parameter and weight values + X_k = [] + W_k = [] for n in range(self.N): X_k.append(self.xi_1d[n][np.abs(diff)[n]]) W_k.append(self.wi_1d[n][np.abs(diff)[n]]) @@ -304,247 +302,247 @@ def compute_Q_diff(self, l, samples, **kwargs): W_k = np.prod(W_k, axis=1) W_k = W_k.reshape([W_k.shape[0], 1]) - #find corresponding code values + # find corresponding code values f_k = [] for x in X_k: j = np.where((x == self.xi_d).all(axis=1))[0][0] - f_k.append(samples[j]) + f_k.append(samples[j]) f_k = np.array(f_k).reshape([len(X_k), self.N_qoi]) - #quadrature of Q^1_{k1} X ... X Q^1_{kN} product - Q_prod = np.sum(f_k*W_k, axis=0).T - Delta += np.sign(np.prod(diff))*Q_prod - + # quadrature of Q^1_{k1} X ... X Q^1_{kN} product + Q_prod = np.sum(f_k * W_k, axis=0).T + Delta += np.sign(np.prod(diff)) * Q_prod + return Delta - + def get_moments(self, qoi): """ Parameters ---------- - qoi (str): name of the qoi - + Returns ------- - mean and variance of qoi (float (N_qoi,)) - + """ - - #compute mean + + # compute mean mean_f = self.quadrature(qoi) - - #compute variance + + # compute variance variance_samples = [] for sample in self.samples[qoi]: variance_samples.append((sample - mean_f)**2) - + var_f = self.quadrature(qoi, samples=variance_samples) - + return mean_f, var_f - + def sc_expansion(self, L, samples, **kwargs): """ ----------------------------------------- This is the UQ Pattern for the SC method. ----------------------------------------- - + Can perform interpolation for both full and sparse grids. - + For a qoi q, it computes the following tensor product: - + q \approx \sum_{l\in\Lambda} \Delta_{l}[q](x) - + where Delta_{l} is the difference at x between surrogates / quadratues of level L and L-1. See e.g.: - - Dimitrios Loukrezis et. al., "Assessing the Performance of Leja and - Clenshaw-Curtis Collocation for Computational Electromagnetics with + + Dimitrios Loukrezis et. al., "Assessing the Performance of Leja and + Clenshaw-Curtis Collocation for Computational Electromagnetics with Random Input Data." - + Parameters ---------- - + - x (float (N,)): location in stochastic space at which to eval the surrogate - L (int): max level of the surrogate - + Returns ------- - - surr (float, (N_qoi,)): the interpolated value of qoi at x - """ - - #for L < L_min the surrogate is defined as zero + + surr (float, (N_qoi,)): the interpolated value of qoi at x + """ + + # for L < L_min the surrogate is defined as zero if L < self.L_min: return 0.0 surr = np.zeros(self.N_qoi) - #contains the level multi-indices (l), colloc points x and samples - #indices f of the (sparse) grid + # contains the level multi-indices (l), colloc points x and samples + # indices f of the (sparse) grid Map = self.Map[L] - + for k in Map.keys(): - #the current code samples + # the current code samples q_k = samples[Map[k]['f']] - #the current level multi index (l_1,...,l_N) + # the current level multi index (l_1,...,l_N) l = Map[k]['l'] - - #the hierarchical surplus (s_k) between the code output q_k and the - #previous surrogate of level L-1 evaluated at the same location. - #Recursively computed. - - if self.sparse == True: - Lm1 = np.sum(l)-1 + + # the hierarchical surplus (s_k) between the code output q_k and the + # previous surrogate of level L-1 evaluated at the same location. + # Recursively computed. + + if self.sparse: + Lm1 = np.sum(l) - 1 else: - Lm1 = self.L-1 - + Lm1 = self.L - 1 + if k in self.surr_lm1[L]: -# print('surrogate already computed') + # print('surrogate already computed') surr_lm1 = self.surr_lm1[L][k] else: - surr_lm1 = self.sc_expansion(Lm1, samples, x = Map[k]['X']) + surr_lm1 = self.sc_expansion(Lm1, samples, x=Map[k]['X']) self.surr_lm1[L][k] = surr_lm1 s_k = q_k - surr_lm1 - + idx = {} # indices of current collocation point (Map[k]['X'][n]), # in corresponding 1d colloc points (self.xi_1d[n][l[n]]) - # These are the j of the 1D lagrange polynomials l_j(x), see + # These are the j of the 1D lagrange polynomials l_j(x), see # lagrange_poly subroutine for n in range(self.N): idx[n] = (self.xi_1d[n][l[n]] == Map[k]['X'][n]).nonzero()[0][0] weight = [] for n in range(self.N): - #interpolate + # interpolate x = kwargs['x'] # add values of Lagrange polynomials at x weight.append(lagrange_poly(x[n], self.xi_1d[n][l[n]], idx[n])) - - #Delta is the interpolation of the hierarchical surplus - surr += s_k*np.prod(weight) + + # Delta is the interpolation of the hierarchical surplus + surr += s_k * np.prod(weight) return surr - + def clear_surr_lm1(self): """ Clears the interpolation results in surr_lm1[ID]. - - surr_lm1 is a dictionary used to store surrogate results at + + surr_lm1 is a dictionary used to store surrogate results at previous level (l-1). Used to avoid recomputing the surrogate in the recursive sc_expansion subroutine. - + surr_lm1[l][k] stores the interpolation results - of level l-1 at collocation point X_k - + of level l-1 at collocation point X_k + Parameters ---------- - ID (str): either 'interpolate' or 'quadrature' - + """ self.surr_lm1 = {} - for level in range(self.L_min, self.L+1): + for level in range(self.L_min, self.L + 1): self.surr_lm1[level] = {} - + def compute_SC_weights(self, rule): """ Computes the 1D quadrature weights w_j of the SC expansion: - + w_j = int L_j(x)p(x) dx (1) - + Here L_j is a Lagrange polynomial of the SC expansion. - + Parameters ---------- - - rule ("str"): chaospy quadrature rule used to compute (1), - - + - rule ("str"): chaospy quadrature rule used to compute (1), + + Returns ------- - wi_1d (dict): wi_1d[n][l] gives an array of quadrature weigths for the n-th parameter at level l. - + IMPORTANT: - If rule is the same as the rule used to compute the SC + If rule is the same as the rule used to compute the SC collocation points, these weights will equal the weights - computed by chaospy, since L_j(x_k) = 1 when j=k and 0 - for the rest. This is the default setting. + computed by chaospy, since L_j(x_k) = 1 when j=k and 0 + for the rest. This is the default setting. """ - - #no need to recompute weights + + # no need to recompute weights if rule == self.sampler.quadrature_rule: return self.sampler.wi_1d - #recompute weights - generally not used + # recompute weights - generally not used else: wi_1d = {} - + params = self.sampler.params_distribution - + for n in range(self.N): - #1d weights for n-th parameter + # 1d weights for n-th parameter wi_1d[n] = {} - #loop over all level of collocation method - for level in range(1, self.L+1): - #current SC nodes over dimension n and level + # loop over all level of collocation method + for level in range(1, self.L + 1): + # current SC nodes over dimension n and level xi_1d = self.xi_1d[n][level] wi_1d[n][level] = np.zeros(xi_1d.size) - - #generate a quadrature rule to compute the SC weights + + # generate a quadrature rule to compute the SC weights xi_quad, wi_quad = cp.generate_quadrature(level, params[n], rule=rule) xi_quad = xi_quad[0] - - #compute integral of the lagrange polynomial through xi_1d, weighted - #by the input distributions: - #w_j = int L_j(xi) p(xi) dxi j = 1,..,xi_1d.size + + # compute integral of the lagrange polynomial through xi_1d, weighted + # by the input distributions: + # w_j = int L_j(xi) p(xi) dxi j = 1,..,xi_1d.size for j in range(xi_1d.size): - #values of L_i(xi_quad) + # values of L_i(xi_quad) lagrange_quad = np.zeros(xi_quad.size) for i in range(xi_quad.size): lagrange_quad[i] = lagrange_poly(xi_quad[i], xi_1d, j) - #quadrature - wi_1d[n][level][j] = np.sum(lagrange_quad*wi_quad) - + # quadrature + wi_1d[n][level][j] = np.sum(lagrange_quad * wi_quad) + return wi_1d - + def get_sample_array(self, qoi): """ Parameters ---------- - qoi (str): name of quantity of interest - + Returns ------- - array of all samples of qoi """ - + tmp = np.zeros([self._number_of_samples, self.N_qoi]) - + for k in range(self._number_of_samples): tmp[k, :] = (self.samples[qoi][k]) - + return tmp - + def plot_grid(self): """ If N = 2 or N = 3 plot the (sparse) grid """ import matplotlib.pyplot as plt - + if self.N == 2: fig = plt.figure() - ax = fig.add_subplot(111, xlabel=r'$x_1$', ylabel=r'$x_2$') - ax.plot(self.xi_d[:,0], self.xi_d[:,1], 'ro') + ax = fig.add_subplot(111, xlabel=r'$x_1$', ylabel=r'$x_2$') + ax.plot(self.xi_d[:, 0], self.xi_d[:, 1], 'ro') elif self.N == 3: from mpl_toolkits.mplot3d import Axes3D fig = plt.figure() - ax = fig.add_subplot(111, projection='3d', xlabel=r'$x_1$', \ - ylabel=r'$x_2$', zlabel=r'$x_3$') - ax.scatter(self.xi_d[:,0], self.xi_d[:,1], self.xi_d[:,2]) + ax = fig.add_subplot(111, projection='3d', xlabel=r'$x_1$', + ylabel=r'$x_2$', zlabel=r'$x_3$') + ax.scatter(self.xi_d[:, 0], self.xi_d[:, 1], self.xi_d[:, 2]) else: print('Will only plot for N = 2 or N = 3.') - + plt.tight_layout() plt.show() @@ -595,7 +593,7 @@ def compute_tensor_prod_u(xi, wi, u, u_prime): def compute_marginal(self, qoi, u, u_prime, diff): """ Computes a marginal integral of the qoi(x) over the dimension defined - by u_prime, for every x value in dimensions u + by u_prime, for every x value in dimensions u Parameters ---------- @@ -603,17 +601,17 @@ def compute_marginal(self, qoi, u, u_prime, diff): - u (array of int): dimensions which are not integrated - u_prime (array of int): dimensions which are integrated - diff (array of int): levels - + Returns - Values of the marginal integral ------- """ - - #1d weights and points of the levels in diff + + # 1d weights and points of the levels in diff xi = [self.xi_1d[n][np.abs(diff)[n]] for n in range(self.N)] wi = [self.wi_1d[n][np.abs(diff)[n]] for n in range(self.N)] - + # compute tensor products and weights in dimension u and u' tmp = self.compute_tensor_prod_u(xi, wi, u, u_prime) xi_d_u = tmp['xi_d_u'] @@ -645,21 +643,21 @@ def compute_marginal(self, qoi, u, u_prime, diff): xi_s[k] = xi_d_u_prime[i_up][idx] idx += 1 - #find the index of the corresponding code sample + # find the index of the corresponding code sample tmp = np.prod(self.xi_d == xi_s, axis=1) idx = np.where(tmp == 1)[0][0] - - #perform quadrature + + # perform quadrature q_k = self.samples[qoi][idx].flatten() h[i_u] += q_k * wi_d_u_prime[i_up].prod() - #return marginal and the weights of dimensions u + # return marginal and the weights of dimensions u return h, wi_d_u - def get_sobol_indices(self, qoi, typ = 'first_order'): + def get_sobol_indices(self, qoi, typ='first_order'): """ Computes Sobol indices using Stochastic Collocation. Method: - Tang (2009), GLOBAL SENSITIVITY ANALYSIS FOR STOCHASTIC COLLOCATION + Tang (2009), GLOBAL SENSITIVITY ANALYSIS FOR STOCHASTIC COLLOCATION EXPANSION. Parameters @@ -672,7 +670,7 @@ def get_sobol_indices(self, qoi, typ = 'first_order'): Either the first order or all Sobol indices of qoi """ - + print('Computing', typ, 'Sobol indices') # multi indices @@ -699,27 +697,27 @@ def get_sobol_indices(self, qoi, typ = 'first_order'): # complement of u u_prime = np.delete(U, u) D_u[u] = 0.0 - + for l in self.l_norm: - - #expand the multi-index indices of the tensor product - #(Q^1_{i1} - Q^1_{i1-1}) X ... X (Q^1_{id) - Q^1_{id-1}) - diff_idx = np.array(list(product(*[[k, -(k-1)] for k in l]))) - - #perform analysis on each Q^1_l1 X ... X Q^1_l_N tensor prod + + # expand the multi-index indices of the tensor product + # (Q^1_{i1} - Q^1_{i1-1}) X ... X (Q^1_{id) - Q^1_{id-1}) + diff_idx = np.array(list(product(*[[k, -(k - 1)] for k in l]))) + + # perform analysis on each Q^1_l1 X ... X Q^1_l_N tensor prod for diff in diff_idx: - + if not np.min(np.abs(diff)) < self.L_min: - - #mariginal integral h, integrate over dimensions u' + + # mariginal integral h, integrate over dimensions u' h, wi_d_u = self.compute_marginal(qoi, u, u_prime, diff) - - #square result and integrate over remaining dimensions u - for i_u in range(wi_d_u.shape[0]): - D_u[u] += np.sign(np.prod(diff))*h[i_u]**2 * wi_d_u[i_u].prod() + + # square result and integrate over remaining dimensions u + for i_u in range(wi_d_u.shape[0]): + D_u[u] += np.sign(np.prod(diff)) * h[i_u]**2 * wi_d_u[i_u].prod() D_u[u] = D_u[u].flatten() - + # all subsets of u W = list(powerset(u))[0:-1] @@ -741,6 +739,7 @@ def get_sobol_indices(self, qoi, typ = 'first_order'): # End SC specific methods + def powerset(iterable): """ powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3) @@ -760,6 +759,7 @@ def powerset(iterable): s = list(iterable) return chain.from_iterable(combinations(s, r) for r in range(len(s) + 1)) + def lagrange_poly(x, x_i, j): """ Lagrange polynomials used for interpolation @@ -791,4 +791,4 @@ def lagrange_poly(x, x_i, j): l_j *= nom / denom - return l_j \ No newline at end of file + return l_j diff --git a/easyvvuq/sampling/stochastic_collocation.py b/easyvvuq/sampling/stochastic_collocation.py index 784ade25c..886495fc4 100644 --- a/easyvvuq/sampling/stochastic_collocation.py +++ b/easyvvuq/sampling/stochastic_collocation.py @@ -49,7 +49,7 @@ def __init__(self, quadrature_rule : char, optional The quadrature method, default is Gaussian "G". - + growth: bool, optional Sets the growth rule to exponential for Clenshaw Curtis quadrature, which makes it nested, and therefore more efficient for sparse grids. @@ -80,14 +80,14 @@ def __init__(self, self.quad_sparse = sparse self.growth = growth self.params_distribution = params_distribution - - #L = level of (sparse) grid + + # L = level of (sparse) grid L = self.quad_order - #N = number of uncertain parameters + # N = number of uncertain parameters N = len(params_distribution) - #for every dimension (parameter), create a hierachy of 1D - #quadrature rules of increasing order + # for every dimension (parameter), create a hierachy of 1D + # quadrature rules of increasing order self.xi_1d = {} self.wi_1d = {} @@ -95,49 +95,50 @@ def __init__(self, self.xi_1d[n] = {} self.wi_1d[n] = {} - if sparse == True: + if sparse: for n in range(N): - for i in range(1, self.quad_order+1): - xi_i, wi_i = cp.generate_quadrature(i+1, - params_distribution[n], - rule=self.quad_rule, + for i in range(1, self.quad_order + 1): + xi_i, wi_i = cp.generate_quadrature(i + 1, + params_distribution[n], + rule=self.quad_rule, growth=self.growth, normalize=True) - + self.xi_1d[n][i] = xi_i[0] self.wi_1d[n][i] = wi_i else: for n in range(N): - xi_i, wi_i = cp.generate_quadrature(self.quad_order, - params_distribution[n], - rule=self.quad_rule, + xi_i, wi_i = cp.generate_quadrature(self.quad_order, + params_distribution[n], + rule=self.quad_rule, growth=self.growth, - normalize=True) + normalize=True) self.xi_1d[n][self.quad_order] = xi_i[0] self.wi_1d[n][self.quad_order] = wi_i - if sparse == False: + if not sparse: # the nodes of the collocation grid xi_d, _ = cp.generate_quadrature(self.quad_order, self.joint_dist, rule=quadrature_rule) self.xi_d = xi_d.T - #sparse grid = a linear combination of tensor products of 1D rules - #of different order. Use chaospy to compute these 1D quadrature rules + # sparse grid = a linear combination of tensor products of 1D rules + # of different order. Use chaospy to compute these 1D quadrature rules else: - #L >= N must hold + # L >= N must hold if L < N: print("*************************************************************") print("Level of sparse grid is lower than the dimension N (# params)") print("Increase level (via polynomial_order) p such that p-1 >= N") print("*************************************************************") - import sys; sys.exit() + import sys + sys.exit() - #multi-index l, such that |l| <= L + # multi-index l, such that |l| <= L l_norm_le_L = self.compute_sparse_multi_idx(L, N) - - #create sparse grid of dimension N and level q using the 1d + + # create sparse grid of dimension N and level q using the 1d #rules in self.xi_1d self.xi_d = self.generate_grid(L, N, l_norm_le_L) @@ -187,7 +188,7 @@ def get_restart_dict(self): "count": self.count, "growth": self.growth, "sparse": self.sparse} - + """ ======================= SPARSE GRID SUBROUTINES @@ -195,31 +196,31 @@ def get_restart_dict(self): """ def generate_grid(self, L, N, l_norm, **kwargs): - + if 'dimensions' in kwargs: dimensions = kwargs['dimensions'] else: dimensions = range(N) - - H_L_N = [] - #loop over all multi indices i + + H_L_N = [] + # loop over all multi indices i for l in l_norm: - - #compute the tensor product of nodes indexed by i + + # compute the tensor product of nodes indexed by i X_l = [self.xi_1d[n][l[n]] for n in dimensions] H_L_N.append(list(product(*X_l))) - - #flatten the list of lists + + # flatten the list of lists H_L_N = np.array(list(chain(*H_L_N))) - - #return unique nodes + + # return unique nodes return np.unique(H_L_N, axis=0) - + def compute_sparse_multi_idx(self, L, N): """ computes all N dimensional multi-indices l = (l1,...,lN) such that |l| <= Q. Here |l| is the internal sum of i (l1+...+lN) """ - P = np.array(list(product(range(1, L+1), repeat=N))) + P = np.array(list(product(range(1, L + 1), repeat=N))) l_norm_le_q = P[np.where(np.sum(P, axis=1) <= L)[0]] - return l_norm_le_q \ No newline at end of file + return l_norm_le_q diff --git a/tests/sc/sparse_sc_test.py b/tests/sc/sparse_sc_test.py index 35099b4a9..70ef5fd36 100644 --- a/tests/sc/sparse_sc_test.py +++ b/tests/sc/sparse_sc_test.py @@ -14,7 +14,7 @@ # author: Wouter Edeling __license__ = "LGPL" -#home directory of user +# home directory of user home = os.path.expanduser('~') HOME = os.path.abspath(os.path.dirname(__file__)) @@ -40,9 +40,9 @@ output_filename = params["out_file"]["default"] output_columns = ["u"] -# Create an encoder, decoder and collation element +# Create an encoder, decoder and collation element encoder = uq.encoders.GenericEncoder( - template_fname = HOME + '/sc.template', + template_fname=HOME + '/sc.template', delimiter='$', target_filename='ade_in.json') decoder = uq.decoders.SimpleCSV(target_filename=output_filename, @@ -73,8 +73,8 @@ of 1D collocation points per level. Used to make e.g. clenshaw-curtis quadrature nested. """ -my_sampler = uq.sampling.SCSampler(vary=vary, polynomial_order = 4, - quadrature_rule="C", sparse=True, +my_sampler = uq.sampling.SCSampler(vary=vary, polynomial_order=4, + quadrature_rule="C", sparse=True, growth=True) # Associate the sampler with the campaign @@ -121,38 +121,38 @@ ax = fig.add_subplot(122, xlabel='x', ylabel='u', title='Surrogate and code samples') - + xi_mc = analysis.xi_d - n_mc = xi_mc.shape[0] + n_mc = xi_mc.shape[0] code_samples = analysis.get_sample_array('u') - + # evaluate the surrogate at these values print('Evaluating surrogate model', n_mc, 'times') for i in range(n_mc): ax.plot(x, analysis.surrogate('u', xi_mc[i]), 'g') ax.plot(x, code_samples[i], 'r+') print('done') - + plt.tight_layout() - + analysis.plot_grid() # ####################### ## Plot Sobol indices # ####################### - + fig = plt.figure() ax = fig.add_subplot( 111, xlabel='x', ylabel='Sobol indices', title='spatial dist. Sobol indices, Pe only important in viscous regions') - + lbl = ['Pe', 'f', 'Pe-f interaction'] idx = 0 - + if analysis.element_name() == 'SC_Analysis': - + for S_i in results['sobol_indices']['u']: ax.plot(x, results['sobol_indices']['u'][S_i], label=lbl[idx]) idx += 1 @@ -160,10 +160,10 @@ for S_i in results['sobol_indices']['u'][1]: ax.plot(x, results['sobol_indices']['u'][1][S_i], label=lbl[idx]) idx += 1 - + leg = plt.legend(loc=0) leg.set_draggable(True) - + plt.tight_layout() -plt.show() \ No newline at end of file +plt.show() From 1e195febb40f1872755b25381c7774b3fea6b70a Mon Sep 17 00:00:00 2001 From: Robin Richardson Date: Fri, 13 Sep 2019 13:11:13 +0100 Subject: [PATCH 41/48] Fix sparse SC test (wrong sc_model.py path for tests) --- tests/sc/sparse_sc_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/sc/sparse_sc_test.py b/tests/sc/sparse_sc_test.py index 70ef5fd36..52fad60d9 100644 --- a/tests/sc/sparse_sc_test.py +++ b/tests/sc/sparse_sc_test.py @@ -86,7 +86,7 @@ # Use this instead to run the samples using EasyVVUQ on the localhost my_campaign.apply_for_each_run_dir(uq.actions.ExecuteLocal( - "./sc_model.py ade_in.json")) + "./tests/sc/sc_model.py ade_in.json")) my_campaign.collate() From e51cc939e8483690525eb34d39df17dd8af135a5 Mon Sep 17 00:00:00 2001 From: di73kuj2 Date: Fri, 13 Sep 2019 14:12:01 +0200 Subject: [PATCH 42/48] Fix csv file header format in the test simulations --- tests/cannonsim/src/cannonsim.cc | 2 +- tests/cooling/cooling_model.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/cannonsim/src/cannonsim.cc b/tests/cannonsim/src/cannonsim.cc index 1db850669..584732e04 100644 --- a/tests/cannonsim/src/cannonsim.cc +++ b/tests/cannonsim/src/cannonsim.cc @@ -100,7 +100,7 @@ int main(int argc, char **argv) exit(1); } - outfile << "# Dist,lastvx,lastvy\n"; + outfile << "Dist,lastvx,lastvy\n"; outfile << out_dist << "," << out_vx << "," << out_vy << "\n"; outfile.close(); diff --git a/tests/cooling/cooling_model.py b/tests/cooling/cooling_model.py index 94ea846cf..7975569c6 100755 --- a/tests/cooling/cooling_model.py +++ b/tests/cooling/cooling_model.py @@ -41,7 +41,7 @@ def f(T, time, kappa, T_env): te = model(t, temp0, kappa, t_env) ti = -te # output csv file -header = 'te, ti' +header = 'te,ti' np.savetxt(output_filename, np.c_[te, ti], delimiter=",", comments='', header=header) From 0df54c6878ea8897090209a8f145f62c8774b6c8 Mon Sep 17 00:00:00 2001 From: Robin Richardson Date: Fri, 13 Sep 2019 13:15:45 +0100 Subject: [PATCH 43/48] Fix pep8 for comment in SC test --- tests/sc/sparse_sc_test.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/sc/sparse_sc_test.py b/tests/sc/sparse_sc_test.py index 52fad60d9..229ebaba8 100644 --- a/tests/sc/sparse_sc_test.py +++ b/tests/sc/sparse_sc_test.py @@ -136,9 +136,9 @@ plt.tight_layout() analysis.plot_grid() - # + ####################### - ## Plot Sobol indices # + # Plot Sobol indices # ####################### fig = plt.figure() From 90cb506b70b5b97a9fe5a06eda2c0714d7ccfe85 Mon Sep 17 00:00:00 2001 From: Robin Richardson Date: Fri, 13 Sep 2019 13:24:34 +0100 Subject: [PATCH 44/48] Try forcing chaospy to be 3.0.7 to fix SC issue (no keyword normalize) --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index a82834046..74d1d7ade 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,7 @@ numpy pandas>=0.24 scipy -chaospy>=3.0.7 +chaospy==3.0.7 SALib pytest pytest-pep8 From 2e7e9cbdd0e545810a46bbc3731f4a25600e3304 Mon Sep 17 00:00:00 2001 From: Robin Richardson Date: Fri, 13 Sep 2019 13:58:19 +0100 Subject: [PATCH 45/48] Remove deprecated normalize keyword from generate_quadrature() --- easyvvuq/sampling/stochastic_collocation.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/easyvvuq/sampling/stochastic_collocation.py b/easyvvuq/sampling/stochastic_collocation.py index 886495fc4..fbcb00c58 100644 --- a/easyvvuq/sampling/stochastic_collocation.py +++ b/easyvvuq/sampling/stochastic_collocation.py @@ -101,8 +101,8 @@ def __init__(self, xi_i, wi_i = cp.generate_quadrature(i + 1, params_distribution[n], rule=self.quad_rule, - growth=self.growth, - normalize=True) + growth=self.growth) +# normalize=True) self.xi_1d[n][i] = xi_i[0] self.wi_1d[n][i] = wi_i @@ -111,8 +111,8 @@ def __init__(self, xi_i, wi_i = cp.generate_quadrature(self.quad_order, params_distribution[n], rule=self.quad_rule, - growth=self.growth, - normalize=True) + growth=self.growth) +# normalize=True) self.xi_1d[n][self.quad_order] = xi_i[0] self.wi_1d[n][self.quad_order] = wi_i From da05351b040f06c3a76b2b1ecf18fb1a186f0cb0 Mon Sep 17 00:00:00 2001 From: Robin Richardson Date: Fri, 13 Sep 2019 14:00:44 +0100 Subject: [PATCH 46/48] Fix mismatch between pdf() signature in uniform_integer relative to the CP Dist base class --- easyvvuq/distributions/distributions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/easyvvuq/distributions/distributions.py b/easyvvuq/distributions/distributions.py index 83d10eaba..6139fbd88 100644 --- a/easyvvuq/distributions/distributions.py +++ b/easyvvuq/distributions/distributions.py @@ -68,7 +68,7 @@ def bnd(self): return self.lo, self.up - def pdf(self): + def pdf(self, step=1e-07): """ Probability density function. """ From c2277f847da5876f9dcf36219380964ad6ace016 Mon Sep 17 00:00:00 2001 From: Robin Richardson Date: Fri, 13 Sep 2019 14:14:01 +0100 Subject: [PATCH 47/48] Add specific versions to every dependency in requirements.txt to try to avoid stable versions being broken by version increments --- requirements.txt | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/requirements.txt b/requirements.txt index c230984bb..932b0c193 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,11 +1,11 @@ -numpy -pandas>=0.24 -scipy -chaospy>=3.0.7 -SALib -pytest -pytest-pep8 -SQLAlchemy -sqlalchemy-utils -jsonpickle -cerberus +numpy==1.16.2 +pandas==0.25 +scipy==1.3.1 +chaospy==3.0.7 +SALib==1.3.8 +pytest==4.3.1 +pytest-pep8==1.0.6 +SQLAlchemy==1.3.8 +sqlalchemy-utils==0.34.2 +jsonpickle==1.2 +cerberus==1.3.1 From 75a57016c39ea53d0ec5f6534ccd69756fbad045 Mon Sep 17 00:00:00 2001 From: Robin Richardson Date: Fri, 13 Sep 2019 14:30:52 +0100 Subject: [PATCH 48/48] Fixed broken deserialization of SC sampler (quad_order -> polynomial_order inconsistency) --- easyvvuq/sampling/stochastic_collocation.py | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/easyvvuq/sampling/stochastic_collocation.py b/easyvvuq/sampling/stochastic_collocation.py index fbcb00c58..e6b18c789 100644 --- a/easyvvuq/sampling/stochastic_collocation.py +++ b/easyvvuq/sampling/stochastic_collocation.py @@ -61,7 +61,6 @@ def __init__(self, """ self.vary = Vary(vary) - #self.polynomial_order = polynomial_order self.quadrature_rule = quadrature_rule # List of the probability distributions of uncertain parameters @@ -73,8 +72,7 @@ def __init__(self, self.joint_dist = cp.J(*params_distribution) # The quadrature information: order, rule and sparsity - #self.quad_order = polynomial_order - self.quad_order = polynomial_order + self.polynomial_order = polynomial_order self.quad_rule = quadrature_rule self.sparse = sparse self.quad_sparse = sparse @@ -82,7 +80,7 @@ def __init__(self, self.params_distribution = params_distribution # L = level of (sparse) grid - L = self.quad_order + L = self.polynomial_order # N = number of uncertain parameters N = len(params_distribution) @@ -97,28 +95,26 @@ def __init__(self, if sparse: for n in range(N): - for i in range(1, self.quad_order + 1): + for i in range(1, self.polynomial_order + 1): xi_i, wi_i = cp.generate_quadrature(i + 1, params_distribution[n], rule=self.quad_rule, growth=self.growth) -# normalize=True) self.xi_1d[n][i] = xi_i[0] self.wi_1d[n][i] = wi_i else: for n in range(N): - xi_i, wi_i = cp.generate_quadrature(self.quad_order, + xi_i, wi_i = cp.generate_quadrature(self.polynomial_order, params_distribution[n], rule=self.quad_rule, growth=self.growth) -# normalize=True) - self.xi_1d[n][self.quad_order] = xi_i[0] - self.wi_1d[n][self.quad_order] = wi_i + self.xi_1d[n][self.polynomial_order] = xi_i[0] + self.wi_1d[n][self.polynomial_order] = wi_i if not sparse: # the nodes of the collocation grid - xi_d, _ = cp.generate_quadrature(self.quad_order, + xi_d, _ = cp.generate_quadrature(self.polynomial_order, self.joint_dist, rule=quadrature_rule) self.xi_d = xi_d.T @@ -183,7 +179,7 @@ def is_restartable(self): def get_restart_dict(self): return { "vary": self.vary.serialize(), - "quad_order": self.quad_order, + "polynomial_order": self.polynomial_order, "quadrature_rule": self.quadrature_rule, "count": self.count, "growth": self.growth,