Skip to content

Commit

Permalink
Merge branch 'dev'
Browse files Browse the repository at this point in the history
  • Loading branch information
@ikroeker
ikroeker committed Jul 6, 2023
2 parents e984c71 + 7663cec commit 85ee56b
Show file tree
Hide file tree
Showing 5 changed files with 172 additions and 53 deletions.
24 changes: 16 additions & 8 deletions aMRPC/datatools.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -814,21 +815,27 @@ 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


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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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'
Expand Down
97 changes: 57 additions & 40 deletions aMRPC/sobol.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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):
Expand All @@ -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:
Expand All @@ -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):
"""
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)

Expand All @@ -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:
Expand All @@ -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
Expand Down Expand Up @@ -419,18 +430,18 @@ 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:
a_r_cf = rsc_dict[a_mkey]
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:
Expand All @@ -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


Expand Down Expand Up @@ -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
Expand All @@ -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):
"""
Expand Down Expand Up @@ -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
Loading

0 comments on commit 85ee56b

Please sign in to comment.