Merge branch 'develop' into feature/eliminate_derived_array

This commit is contained in:
Fabian Joswig 2021-12-06 15:30:14 +00:00
commit 3ae6959bac
6 changed files with 72 additions and 22 deletions

View file

@ -301,8 +301,7 @@ def total_least_squares(x, y, func, silent=False, **kwargs):
result = []
for i in range(n_parms):
result.append(derived_observable(lambda x, **kwargs: x[0], list(x.ravel()) + list(y), man_grad=list(deriv_x[i]) + list(deriv_y[i])))
result[-1]._value = out.beta[i]
result.append(derived_observable(lambda my_var, **kwargs: my_var[0] / x.ravel()[0].value * out.beta[i], list(x.ravel()) + list(y), man_grad=list(deriv_x[i]) + list(deriv_y[i])))
output.fit_parameters = result + const_par
@ -419,8 +418,7 @@ def _prior_fit(x, y, func, priors, silent=False, **kwargs):
result = []
for i in range(n_parms):
result.append(derived_observable(lambda x, **kwargs: x[0], list(y) + list(loc_priors), man_grad=list(deriv[i])))
result[-1]._value = params[i]
result.append(derived_observable(lambda x, **kwargs: x[0] / y[0].value * params[i], list(y) + list(loc_priors), man_grad=list(deriv[i])))
output.fit_parameters = result
output.chisquare = chisqfunc(np.asarray(params))
@ -614,8 +612,7 @@ def _standard_fit(x, y, func, silent=False, **kwargs):
result = []
for i in range(n_parms):
result.append(derived_observable(lambda x, **kwargs: x[0], list(y), man_grad=list(deriv[i])))
result[-1]._value = fit_result.x[i]
result.append(derived_observable(lambda x, **kwargs: x[0] / y[0].value * fit_result.x[i], list(y), man_grad=list(deriv[i])))
output.fit_parameters = result + const_par

View file

@ -74,10 +74,15 @@ class Obs:
if idl is not None:
if len(idl) != len(names):
raise Exception('Length of idl incompatible with samples and names.')
if len(names) != len(set(names)):
raise Exception('names are not unique.')
if not all(isinstance(x, str) for x in names):
raise TypeError('All names have to be strings.')
name_length = len(names)
if name_length > 1:
if name_length != len(set(names)):
raise Exception('names are not unique.')
if not all(isinstance(x, str) for x in names):
raise TypeError('All names have to be strings.')
else:
if not isinstance(names[0], str):
raise TypeError('All names have to be strings.')
if min(len(x) for x in samples) <= 4:
raise Exception('Samples have to have at least 5 entries.')
@ -623,9 +628,9 @@ class Obs:
name = self.names[0]
full_data = self.deltas[name] + self.r_values[name]
n = full_data.size
mean = np.mean(full_data)
mean = self.value
tmp_jacks = np.zeros(n + 1)
tmp_jacks[0] = self.value
tmp_jacks[0] = mean
tmp_jacks[1:] = (n * mean - full_data) / (n - 1)
return tmp_jacks
@ -1259,22 +1264,22 @@ def reweight(weight, obs, **kwargs):
for i in range(len(obs)):
if len(obs[i].cov_names):
raise Exception('Error: Not possible to reweight an Obs that contains covobs!')
if sorted(weight.names) != sorted(obs[i].names):
if not set(obs[i].names).issubset(weight.names):
raise Exception('Error: Ensembles do not fit')
for name in weight.names:
for name in obs[i].names:
if not set(obs[i].idl[name]).issubset(weight.idl[name]):
raise Exception('obs[%d] has to be defined on a subset of the configs in weight.idl[%s]!' % (i, name))
new_samples = []
w_deltas = {}
for name in sorted(weight.names):
for name in sorted(obs[i].names):
w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name])
new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name]))
tmp_obs = Obs(new_samples, sorted(weight.names), idl=[obs[i].idl[name] for name in sorted(weight.names)])
tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
if kwargs.get('all_configs'):
new_weight = weight
else:
new_weight = Obs([w_deltas[name] + weight.r_values[name] for name in sorted(weight.names)], sorted(weight.names), idl=[obs[i].idl[name] for name in sorted(weight.names)])
new_weight = Obs([w_deltas[name] + weight.r_values[name] for name in sorted(obs[i].names)], sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
result.append(derived_observable(lambda x, **kwargs: x[0] / x[1], [tmp_obs, new_weight], **kwargs))
result[-1].reweighted = True

View file

@ -33,6 +33,5 @@ def find_root(d, func, guess=1.0, **kwargs):
da = jacobian(lambda u, v: func(v, u))(d.value, root[0])
deriv = - da / dx
res = derived_observable(lambda x, **kwargs: x[0], [d], man_grad=[deriv])
res._value = root[0]
res = derived_observable(lambda x, **kwargs: x[0] / d.value * root[0], [d], man_grad=[deriv])
return res

View file

@ -52,7 +52,7 @@ def test_least_squares():
outc = pe.least_squares(x, oyc, func)
betac = outc.fit_parameters
for i in range(2):
betac[i].gamma_method(S=1.0)
assert math.isclose(betac[i].value, popt[i], abs_tol=1e-5)
@ -97,7 +97,7 @@ def test_least_squares():
return p[1] * np.exp(-p[0] * x)
fitp = pe.least_squares(x, data, fitf, expected_chisquare=True)
fitpc = pe.least_squares(x, data, fitf, correlated_fit=True)
for i in range(2):
diff = fitp[i] - fitpc[i]
@ -170,7 +170,7 @@ def test_total_least_squares():
diffc = outc.fit_parameters[0] - betac[0]
assert(diffc / betac[0] < 1e-3 * betac[0].dvalue)
assert((outc.fit_parameters[1] - betac[1]).is_zero())
outc = pe.total_least_squares(oxc, oy, func)
betac = outc.fit_parameters
@ -208,3 +208,41 @@ def test_odr_derivatives():
tfit = pe.fits.fit_general(x, y, func, base_step=0.1, step_ratio=1.1, num_steps=20)
assert np.abs(np.max(np.array(list(fit1[1].deltas.values()))
- np.array(list(tfit[1].deltas.values())))) < 10e-8
def test_r_value_persistence():
def f(a, x):
return a[0] + a[1] * x
a = pe.pseudo_Obs(1.1, .1, 'a')
assert np.isclose(a.value, a.r_values['a'])
a_2 = a ** 2
assert np.isclose(a_2.value, a_2.r_values['a'])
b = pe.pseudo_Obs(2.1, .2, 'b')
y = [a, b]
[o.gamma_method() for o in y]
fitp = pe.fits.least_squares([1, 2], y, f)
assert np.isclose(fitp[0].value, fitp[0].r_values['a'])
assert np.isclose(fitp[0].value, fitp[0].r_values['b'])
assert np.isclose(fitp[1].value, fitp[1].r_values['a'])
assert np.isclose(fitp[1].value, fitp[1].r_values['b'])
fitp = pe.fits.total_least_squares(y, y, f)
assert np.isclose(fitp[0].value, fitp[0].r_values['a'])
assert np.isclose(fitp[0].value, fitp[0].r_values['b'])
assert np.isclose(fitp[1].value, fitp[1].r_values['a'])
assert np.isclose(fitp[1].value, fitp[1].r_values['b'])
fitp = pe.fits.least_squares([1, 2], y, f, priors=y)
assert np.isclose(fitp[0].value, fitp[0].r_values['a'])
assert np.isclose(fitp[0].value, fitp[0].r_values['b'])
assert np.isclose(fitp[1].value, fitp[1].r_values['a'])
assert np.isclose(fitp[1].value, fitp[1].r_values['b'])

View file

@ -381,6 +381,16 @@ def test_merge_obs():
assert diff == -(my_obs1.value + my_obs2.value) / 2
def test_merge_obs_r_values():
a1 = pe.pseudo_Obs(1.1, .1, 'a|1')
a2 = pe.pseudo_Obs(1.2, .1, 'a|2')
a = pe.merge_obs([a1, a2])
assert np.isclose(a.r_values['a|1'], a1.value)
assert np.isclose(a.r_values['a|2'], a2.value)
assert np.isclose(a.value, np.mean([a1.value, a2.value]))
def test_correlate():
my_obs1 = pe.Obs([np.random.rand(100)], ['t'])
my_obs2 = pe.Obs([np.random.rand(100)], ['t'])

View file

@ -15,6 +15,7 @@ def test_root_linear():
my_root = pe.roots.find_root(my_obs, root_function)
assert np.isclose(my_root.value, value)
assert np.isclose(my_root.value, my_root.r_values['t'])
difference = my_obs - my_root
assert difference.is_zero()