in syne_tune/optimizer/schedulers/searchers/bayesopt/gpautograd/learncurve/issm.py [0:0]
def issm_likelihood_slow_computations(
targets: List[np.ndarray], issm_params: Dict, r_min: int, r_max: int,
skip_c_d: bool = False,
profiler: Optional[SimpleProfiler] = None) -> Dict:
"""
Naive implementation of `issm_likelihood_computations`, which does not
require precomputations, but is much slower. Here, results are computed
one datapoint at a time, instead of en bulk.
This code is used in unit testing, and called from `sample_posterior_joint`.
"""
num_configs = len(targets)
num_res = r_max + 1 - r_min
assert num_configs > 0, "targets must not be empty"
assert num_res > 0, f"r_min = {r_min} must be <= r_max = {r_max}"
compute_wtw = targets[0].shape[1] == 1
alphas = _flatvec(issm_params['alpha'])
betas = _flatvec(issm_params['beta'])
gamma = issm_params['gamma']
n = getval(alphas.shape[0])
assert n == num_configs, f"alpha.size = {n} != {num_configs}"
n = getval(betas.shape[0])
assert n == num_configs, f"beta.size = {n} != {num_configs}"
# Outer loop over configurations
c_lst = []
d_lst = []
vtv_lst = []
wtv_lst = []
wtw_lst = []
num_data = 0
for i, ymat in enumerate(targets):
alpha = alphas[i]
alpha_m1 = alpha - 1.0
beta = betas[i]
ydim = ymat.shape[0]
if profiler is not None:
profiler.start('issm_part1')
num_data += ydim
r_obs = r_min + ydim # Observed in range(r_min, r_obs)
assert 0 < ydim <= num_res,\
f"len(y[{i}]) = {ydim}, num_res = {num_res}"
if not skip_c_d:
# c_i, d_i
if ydim < num_res:
lrvec = anp.array(
[np.log(r) for r in range(r_obs, r_max + 1)]) *\
alpha_m1 + beta
c_scal = alpha * anp.exp(logsumexp(lrvec))
d_scal = anp.square(gamma * alpha) * anp.exp(
logsumexp(lrvec * 2.0))
c_lst.append(c_scal)
d_lst.append(d_scal)
else:
c_lst.append(0.0)
d_lst.append(0.0)
# Inner loop for v_i, w_i
if profiler is not None:
profiler.stop('issm_part1')
profiler.start('issm_part2')
yprev = ymat[-1].reshape((1, -1)) # y_{j-1} (vector)
vprev = 1.0 # v_{j-1} (scalar)
wprev = yprev # w_{j-1} (row vector)
vtv = vprev * vprev # scalar
wtv = wprev * vprev # row vector
if compute_wtw:
wtw = wprev * wprev # shape (1, 1)
for j in range(1, ydim):
ycurr = ymat[ydim - j - 1].reshape((1, -1)) # y_j (row vector)
# a_{j-1}
ascal = alpha * anp.exp(np.log(r_obs - j) * alpha_m1 + beta)
escal = gamma * ascal + 1.0
vcurr = escal * vprev # v_j
wcurr = escal * wprev + ycurr - yprev + ascal # w_j
vtv = vtv + vcurr * vcurr
wtv = wtv + wcurr * vcurr
if compute_wtw:
wtw = wtw + wcurr * wcurr
yprev = ycurr
vprev = vcurr
wprev = wcurr
vtv_lst.append(vtv)
wtv_lst.append(wtv)
if compute_wtw:
assert wtw.shape == (1, 1)
wtw_lst.append(wtw.item())
if profiler is not None:
profiler.stop('issm_part2')
# Compile results
result = {
'num_data': num_data,
'vtv': anp.array(vtv_lst),
'wtv': anp.vstack(wtv_lst)}
if compute_wtw:
result['wtw'] = anp.array(wtw_lst)
if not skip_c_d:
result['c'] = anp.array(c_lst)
result['d'] = anp.array(d_lst)
return result