diff --git a/aMRPC/datatools.py b/aMRPC/datatools.py index a3acba1..2b491ef 100644 --- a/aMRPC/datatools.py +++ b/aMRPC/datatools.py @@ -13,6 +13,7 @@ import pandas as pd import numpy as np from scipy.linalg import lstsq +# from multiprocessing import freeze_support from . import polytools as pt from . import utils as u # from . import wavetools as wt @@ -751,7 +752,7 @@ def gen_phi_as(pol_vals, alpha_mask): return phi -#@njit +# @njit def gen_cov_mx_4lh(phi, s_sigma_n, s_sigma_p): """ generates covariance etc. matrixes for likelihood @@ -814,13 +815,19 @@ def gen_cov_mx_4lh_noex(phi, s_sigma_n, s_sigma_p): -------- cov_mx_inv inverse cov. matrix """ -# if s_sigma_n < 1e-42 or s_sigma_p < 1e-42: -# print("sigma n/p", s_sigma_n, s_sigma_p) -# return np.nan, np.nan Q = np.eye(phi.shape[0]) * s_sigma_n R = np.eye(phi.shape[1]) * s_sigma_p cov_mx = phi @ R @ phi.T + Q - cov_mx_inv = np.linalg.pinv(cov_mx) + # try: + # cov_mx_inv = np.linalg.pinv(cov_mx) + # except: + # inverse according to eq. (A.9) in Rasmussen and Williams + Q_inv = np.ascontiguousarray(np.eye(phi.shape[0]) / s_sigma_n) + R_inv = np.ascontiguousarray(np.eye(phi.shape[1]) / s_sigma_p) + P = phi.T @ Q_inv @ phi + R_inv + P_inv = np.ascontiguousarray(np.linalg.pinv(P)) + # inverse according to eq. (A.9) in Rasmussen and Williams + cov_mx_inv = Q_inv - Q_inv @ phi @ P_inv @ phi.T @ Q_inv return cov_mx, cov_mx_inv @@ -828,7 +835,7 @@ def sample_amrpc_rec(samples, mk_list, alphas, f_cfs, f_cov_mx, npc_dict, nrb_dict, mk2sid, alpha_masks=None, **kwargs): """ - Generates function reconstruction + Samples function reconstruction (slow version) f(sample, x) = sum_(p in alphas) f_cfs(sample_mk, p,x)*pol(alpha_p, sample) Parameters @@ -1080,7 +1087,7 @@ def pc_eval_mv_par(pcfs, X): dtype=np.float64) -def sample_amprc_cfs(mk_list, alphas, f_cfs, f_cov_mx, +def sample_amrpc_cfs(mk_list, alphas, f_cfs, f_cov_mx, mk2sid, alpha_masks=None, **kwargs): """ samples the polynomial-coeficients f_cfs(sample_mk, p, x) for @@ -1484,7 +1491,7 @@ def gen_amrpc_dec_ls_mask_aux(data, sids, pol_vals, alpha_mask, cov_mask, #v_ls, resid, rank, sigma = np.linalg.lstsq( # Phi, data[sids, idx_x], rcond=None) # LS - output v_ls, _, _, _ = np.linalg.lstsq(phi, rs_data, - rcond=-1) # LS - output + rcond=-1) # LS - output else: v_ls = np.ravel(data[0, dt_idx_x]/phi) @@ -1860,6 +1867,7 @@ def update_pol_vals_on_samples(samples_updated, new_samples_cnt, pol_vals, def main(): + # freeze_support() """ some tests """ # data location # url = '../data/InputParameters.txt' diff --git a/aMRPC/sobol.py b/aMRPC/sobol.py index 2545152..736be11 100644 --- a/aMRPC/sobol.py +++ b/aMRPC/sobol.py @@ -16,6 +16,8 @@ import aMRPC.utils as u # Sobol idx for aPC + + def sobol_idx_pc(pc_coefs, alphas, idx_set, eps=1e-15): """ computes Sobol indexes from polynomial chaos (PC) expansion @@ -53,7 +55,7 @@ def sobol_idx_pc(pc_coefs, alphas, idx_set, eps=1e-15): srcs = list(range(dim)) assert max(idx_set) < dim not_in_idx_set = list(set(srcs)-set(idx_set)) - #sobol = 0 + # sobol = 0 for pidx in range(p_max): alpha = alphas[pidx, :] chk_in = alpha[idx_set].min() > 0 @@ -65,6 +67,7 @@ def sobol_idx_pc(pc_coefs, alphas, idx_set, eps=1e-15): sobol += pc_coefs[pidx]**2 / var_pc return sobol + def gen_idx_subsets(dim): """ Generates all possible source combination for dim-sources @@ -85,6 +88,7 @@ def gen_idx_subsets(dim): for t in it.combinations(items, length)] return set(sub_idx) + def gen_sidx_subsets(sidx): """ Generates all possible source / index combination of sidx @@ -105,6 +109,7 @@ def gen_sidx_subsets(sidx): for t in it.combinations(items, length)] return set(sub_idx) + def sobol_tot_sen_pc(pc_coefs, alphas, src_idxs, idx_list): """ Computes total sensitivity coeficients for polynomial-chaos (PC) expansion @@ -133,6 +138,7 @@ def sobol_tot_sen_pc(pc_coefs, alphas, src_idxs, idx_list): tot += sobol_idx_pc(pc_coefs, alphas, list(idx_it)) return tot + def sobol_tot_sens(sob_dict, src_idxs, idx_list): """ total sensitivity @@ -159,6 +165,7 @@ def sobol_tot_sens(sob_dict, src_idxs, idx_list): tot += sob_dict[idx_it] return tot + def sobol_idx_amrpc_helper(pc_coefs, rsc_dict, mk2sid, alphas, idx_list, eps=1e-15): """ Helper function for computation of Sobol indexes for aMR-PC expansion @@ -190,7 +197,7 @@ def sobol_idx_amrpc_helper(pc_coefs, rsc_dict, mk2sid, alphas, idx_list, eps=1e- srcs = list(range(dim)) assert max(idx_list) < dim not_in_idx_set = list(set(srcs)-set(idx_list)) - #loc_rsc_cf = 2**(- len(not_in_idx_set)*4) + # loc_rsc_cf = 2**(- len(not_in_idx_set)*4) # omit normalization for var <=eps var_thresh = var <= eps if max(var_thresh): @@ -206,24 +213,24 @@ def sobol_idx_amrpc_helper(pc_coefs, rsc_dict, mk2sid, alphas, idx_list, eps=1e- for mkey, sids in mk2sid.items(): loc_pc = pc_coefs[sids[0], :] - sobol_mk = 0 #loc_pc[0]**2 + sobol_mk = 0 # loc_pc[0]**2 r_cf = rsc_dict[mkey] c_cf = u.gen_corr_rcf(mkey, not_in_idx_set) for a_mkey, a_sids in mk2sid.items(): loc_pc_a = pc_coefs[a_sids[0], :] a_r_cf = rsc_dict[a_mkey] a_c_cf = u.gen_corr_rcf(a_mkey, not_in_idx_set) - #idx_diff = u.multi_key_intersect_srcs(mkey, a_mkey) - #if set(idx_diff) <= set(idx_set): - #mk_chk = mkey == a_mkey + # idx_diff = u.multi_key_intersect_srcs(mkey, a_mkey) + # if set(idx_diff) <= set(idx_set): + # mk_chk = mkey == a_mkey # sobol_mk += loc_pc[0]**2 - #else: + # else: mk_chk = u.compare_multi_key_for_idx(mkey, a_mkey, idx_list) if mk_chk: - #sobol_mk += loc_pc[0] * loc_pc_a[0] + # sobol_mk += loc_pc[0] * loc_pc_a[0] for pidx in range(p_max): alpha = alphas[pidx, :] - #chk_in = alpha[idx_list].min() > 0 + # chk_in = alpha[idx_list].min() > 0 if not not_in_idx_set: chk_out = True else: @@ -234,9 +241,10 @@ def sobol_idx_amrpc_helper(pc_coefs, rsc_dict, mk2sid, alphas, idx_list, eps=1e- sobol_mk *= math.sqrt(r_cf * c_cf) sobol_ns += sobol_mk sobol_ns -= mean**2 - #print(mean, var, 1/rsc_dict[mkey]) + # print(mean, var, 1/rsc_dict[mkey]) return sobol_ns / var - #return sobol_ns + # return sobol_ns + def sobol_idx_amrpc(sobol_dict, idx_set): """ @@ -257,17 +265,18 @@ def sobol_idx_amrpc(sobol_dict, idx_set): """ idx_set_len = len(idx_set) - ret_val = 0 #sobol_dict[idx_set] + ret_val = 0 # sobol_dict[idx_set] if idx_set_len > 1: sub_idx = gen_sidx_subsets(idx_set) - #print(idx_set, sub_idx) + # print(idx_set, sub_idx) for idx in sub_idx: - #print(((-1)**(len(idx_set-idx))), idx_set, idx) + # print(((-1)**(len(idx_set-idx))), idx_set, idx) ret_val += ((-1)**(len(idx_set-idx)))*sobol_dict[idx] else: ret_val = sobol_dict[idx_set] return ret_val + def sobol_idx_amrpc_jj(pc_coefs, rsc_dict, mk2sid, alphas, idx_list, eps=1e-15): """ Generates auxiliary values for aux. dict. for computing of Sobol coeficients @@ -298,6 +307,7 @@ def sobol_idx_amrpc_jj(pc_coefs, rsc_dict, mk2sid, alphas, idx_list, eps=1e-15): else: return sobol_idx_amrpc_jj_4s(pc_coefs, rsc_dict, mk2sid, alphas, idx_list, eps) + def sobol_idx_amrpc_jj_4s(pc_coefs, rsc_dict, mk2sid, alphas, idx_list, eps=1e-15): """ Generates auxiliary values for aux. dict. for computing of Sobol coeficients @@ -344,10 +354,10 @@ def sobol_idx_amrpc_jj_4s(pc_coefs, rsc_dict, mk2sid, alphas, idx_list, eps=1e-1 for mkey, sids in mk2sid.items(): loc_pc = pc_coefs[sids[0], :] - sobol_mk = 0 #loc_pc[0]**2 + sobol_mk = 0 # loc_pc[0]**2 r_cf = rsc_dict[mkey] c_cf = u.gen_corr_rcf(mkey, not_in_idx_set) - #c_cf = u.gen_corr_rcf(mkey, not_in_idx_set) + # c_cf = u.gen_corr_rcf(mkey, not_in_idx_set) for a_mkey, a_sids in mk2sid.items(): mk_chk = u.compare_multi_key_for_idx(mkey, a_mkey, idx_list) @@ -357,7 +367,7 @@ def sobol_idx_amrpc_jj_4s(pc_coefs, rsc_dict, mk2sid, alphas, idx_list, eps=1e-1 a_c_cf = u.gen_corr_rcf(a_mkey, not_in_idx_set) for pidx in range(p_max): alpha = alphas[pidx, :] - #chk_in = alpha[idx_list].min() > 0 + # chk_in = alpha[idx_list].min() > 0 if not not_in_idx_set: chk_out = True else: @@ -370,9 +380,10 @@ def sobol_idx_amrpc_jj_4s(pc_coefs, rsc_dict, mk2sid, alphas, idx_list, eps=1e-1 sobol_mk *= math.sqrt(r_cf*c_cf) sobol_ns += sobol_mk sobol_ns -= mean**2 - #var_ths = var >= eps + # var_ths = var >= eps return sobol_ns / var + def sobol_idx_amrpc_jj_4mkey(pc_coefs, rsc_dict, alphas, idx_list, eps=1e-15): """ Generates auxiliary values for aux. dict. for computing of Sobol coeficients @@ -419,10 +430,10 @@ def sobol_idx_amrpc_jj_4mkey(pc_coefs, rsc_dict, alphas, idx_list, eps=1e-15): sobol_ns = np.zeros(mean.shape) for mkey, loc_pc in pc_coefs.items(): - sobol_mk = 0 #loc_pc[0]**2 + sobol_mk = 0 # loc_pc[0]**2 r_cf = rsc_dict[mkey] c_cf = u.gen_corr_rcf(mkey, not_in_idx_set) - #c_cf = u.gen_corr_rcf(mkey, not_in_idx_set) + # c_cf = u.gen_corr_rcf(mkey, not_in_idx_set) for a_mkey, a_loc_pc in pc_coefs.items(): mk_chk = u.compare_multi_key_for_idx(mkey, a_mkey, idx_list) if mk_chk: @@ -430,7 +441,7 @@ def sobol_idx_amrpc_jj_4mkey(pc_coefs, rsc_dict, alphas, idx_list, eps=1e-15): a_c_cf = u.gen_corr_rcf(a_mkey, not_in_idx_set) for pidx in range(p_max): alpha = alphas[pidx, :] - #chk_in = alpha[idx_list].min() > 0 + # chk_in = alpha[idx_list].min() > 0 if not not_in_idx_set: chk_out = True else: @@ -443,7 +454,7 @@ def sobol_idx_amrpc_jj_4mkey(pc_coefs, rsc_dict, alphas, idx_list, eps=1e-15): sobol_mk *= math.sqrt(r_cf*c_cf) sobol_ns += sobol_mk sobol_ns -= mean**2 - #var_ths = var >= eps + # var_ths = var >= eps return sobol_ns / var @@ -478,12 +489,13 @@ def gen_sobol_amrpc_dict(pc_coefs, rsc_dict, mk2sid, alphas, idx_list, eps=1e-15 sub_idx = gen_sidx_subsets(set(idx_list)) sobol_dict = {} for j in sub_idx: - #print(j) + # print(j) sobol_dict[j] = sobol_idx_amrpc_jj(pc_coefs, rsc_dict, mk2sid, alphas, list(j), eps) return sobol_dict + def sobol_idx_amrpc_comb(help_sobol_dict, idx_set, tmp_sobol_dict): """ Computes Sobol indexes for aMR-PC expansion using dictionary provided by @@ -505,26 +517,28 @@ def sobol_idx_amrpc_comb(help_sobol_dict, idx_set, tmp_sobol_dict): idx_set_len = len(idx_set) ret_val = help_sobol_dict[idx_set] if idx_set_len > 1: - #sub_idx = gen_sidx_subsets(list(idx_set)) - idx_set + # sub_idx = gen_sidx_subsets(list(idx_set)) - idx_set items = list(idx_set) sub_idx = [frozenset(t) for length in range(1, idx_set_len) - for t in it.combinations(items, length)] - #print(idx_set, sub_idx) + for t in it.combinations(items, length)] + # print(idx_set, sub_idx) for j_idx in sub_idx: j_len = len(j_idx) - #print('jj:',j_idx, help_sobol_dict[j_idx], j_len) + # print('jj:',j_idx, help_sobol_dict[j_idx], j_len) if j_idx not in tmp_sobol_dict.keys(): if j_len == 1: tmp_sobol_dict[j_idx] = help_sobol_dict[j_idx] else: _, tmp_sobol_dict = sobol_idx_amrpc_comb(help_sobol_dict, - j_idx, tmp_sobol_dict) + j_idx, + tmp_sobol_dict) ret_val -= tmp_sobol_dict[j_idx] - #print('ret=', ret_val) + # print('ret=', ret_val) tmp_sobol_dict[idx_set] = ret_val return ret_val, tmp_sobol_dict + def sobol_idx_amrpc_dynamic(idx_set, pc_coefs, rsc_dict, mk2sid, alphas, sobol_dict, help_sobol_dict, eps=1e-15): """ @@ -572,25 +586,28 @@ def sobol_idx_amrpc_dynamic(idx_set, pc_coefs, rsc_dict, mk2sid, alphas, items = list(idx_set) sub_idx = [frozenset(t) for length in range(1, idx_set_len) for t in it.combinations(items, length)] - #print(idx_set, sub_idx) + # print(idx_set, sub_idx) for j_idx in sub_idx: j_len = len(j_idx) - #print('jj:',j_idx, help_sobol_dict[j_idx], j_len) + # print('jj:',j_idx, help_sobol_dict[j_idx], j_len) if j_idx not in sobol_dict.keys(): if j_len == 1: if j_idx not in help_sobol_dict.keys(): - help_sobol_dict[j_idx] = sobol_idx_amrpc_jj(pc_coefs, rsc_dict, - mk2sid, alphas, - list(j_idx), eps) + help_sobol_dict[j_idx] = sobol_idx_amrpc_jj(pc_coefs, + rsc_dict, + mk2sid, + alphas, + list(j_idx), + eps) sobol_dict[j_idx] = help_sobol_dict[j_idx] else: _, sobol_dict, help_sobol_dict = sobol_idx_amrpc_dynamic(j_idx, pc_coefs, - rsc_dict, mk2sid, - alphas, - sobol_dict, - help_sobol_dict, - eps) + rsc_dict, mk2sid, + alphas, + sobol_dict, + help_sobol_dict, + eps) ret_val -= sobol_dict[j_idx] - #print('ret=', ret_val) + # print('ret=', ret_val) sobol_dict[idx_set] = ret_val return ret_val, sobol_dict, help_sobol_dict diff --git a/aMRPC/stat.py b/aMRPC/stat.py index 07d0c4c..4fe4c9c 100755 --- a/aMRPC/stat.py +++ b/aMRPC/stat.py @@ -362,6 +362,41 @@ def lbme_norm_response(observation, response_surfaces, covariance_matrix, return llh_cf + np.log(mean_lh) if mean_lh > 0 else -np.inf +def lbme_norm_response_inv(observation, response_surfaces, + covariance_matrix_inv, llh_cf, eps=0.0): + """ + Computes log(BME) (Bayesian Model Evidence) + + Parameters + ---------- + observation : np.array + observation trajectory. + response_surfaces : np.array + surrogate / model response. + covariance_matrix_inv : np.array + inverse covariance matrix. + llh_cf : np.array + log-likelihood coefficient + Returns + ------- + float + BME. + + """ + n, m = response_surfaces.shape[0:2] + if m == len(observation): + sample_cnt = n + else: + sample_cnt = m + response_surfaces = response_surfaces.T + + mean_lh = np.mean(np.array([cmp_norm_likelihood_core_inv(observation, + response_surfaces[sample, :], + covariance_matrix_inv) + eps + for sample in range(sample_cnt)], dtype=np.float64)) + return llh_cf + np.log(mean_lh) if mean_lh > 0 else -np.inf + + def d_kl_norm_prior_response(observation, response_surfaces, covariance_matrix, eps=0): """ @@ -411,6 +446,46 @@ def d_kl_norm_prior_response(observation, response_surfaces, covariance_matrix, return np.mean(llhs[mask]) - np.log(bme) if bme > 0 else np.nan +def d_kl_norm_prior_response_inv(observation, response_surfaces, + covariance_matrix_inv, + eps=0): + """ + computes Kullback-Leibler divergence + + Parameters + ---------- + observation : np.array + observation. + response_surfaces : np.array + surrogate / model response. + covariance_matrix_inv : np.array + invserse covariance matrix. + **kwargs : dict + eps - error bad condition compensation. + + Returns + ------- + float. + + """ + # eps = kwargs.get('eps', 0) + n, m = response_surfaces.shape + if m == len(observation): + sample_cnt = n + else: + sample_cnt = m + response_surfaces = response_surfaces.T + + llhs = np.array([cmp_log_likelihood_core_sinv(observation, + response_surfaces[sample, :], + covariance_matrix_inv) + for sample in range(sample_cnt)]) + + mask = llhs - llhs.max() >= np.log(np.random.uniform(0, 1, llhs.shape)) + bme = np.exp(llhs).mean() + eps + + return np.mean(llhs[mask]) - np.log(bme) if bme > 0 else np.nan + # if NJM: # jit_module(nopython=True, error_model="numpy") diff --git a/aMRPC/utils.py b/aMRPC/utils.py index 5620c6f..d4dded5 100644 --- a/aMRPC/utils.py +++ b/aMRPC/utils.py @@ -77,7 +77,7 @@ def gen_midx_mask(alphas, no_max): are below ( <=)no_max """ p_cnt = alphas.shape[0] - a_mask = np.zeros(p_cnt, dtype=np.bool8) + a_mask = np.empty(p_cnt, dtype=np.bool8) for i in range(p_cnt): a_mask[i] = alphas[i, :].sum() <= no_max return a_mask @@ -100,6 +100,24 @@ def gen_midx_mask_part(alphas, no_min, no_max, idx_set): return a_mask +@jit(nopython=True, nogil=True) +def gen_midx_mask_hyp(alphas, no_max, p_norm): + """ + generates a mask for alphas, such that all P-Norms of + multi-index polynomial degrees + are below ( <=) no_max + see Hyperbolic trunction in + Sparse Polynomial Chaos Expansions: Literature Survey and Benchmark + Nora Lüthen, Stefano Marelli, and Bruno Sudret + """ + p_cnt = alphas.shape[0] + a_mask = np.empty(p_cnt, dtype=np.bool8) + for i in range(p_cnt): + i_arr = alphas[i, :].astype(np.float32) + a_mask[i] = np.linalg.norm(i_arr, p_norm) <= no_max + return a_mask + + def gen_nri_range(nrs): """ generates an array with Nri-entries diff --git a/setup.py b/setup.py index 5b4501f..b9c8a85 100644 --- a/setup.py +++ b/setup.py @@ -4,22 +4,23 @@ long_description = fh.read() setuptools.setup( - name="aMRPC-pkg", - version="1.0.1", + name="aMRPC", + version="1.0.2", author="Ilja Kröker", author_email="ilja.kroeker@iws.uni-stuttgart.de", description="aMRPC python implementation", keywords="arbirtray multi-resolution polynomial chaos, multi-wavelet", long_description=long_description, #long_description_content_type="text/markdown", - url="https://git.iws.uni-stuttgart.de/ikroeker/ik_amr-pc", + #url="https://git.iws.uni-stuttgart.de/ikroeker/ik_amr-pc", + url="https://github.com/ikroeker/aMR-PC", packages=setuptools.find_packages(), classifiers=[ "Programming Language :: Python :: 3", "License :: OSI Approved :: MIT License", "Operating System :: OS Independent", ], - python_requires='>=3.6', + python_requires='>=3.10', requires=[ "numpy", "pandas", "scipy", "math", "pickle", "os", "itertools", "numba"