Repository: fjosw/pyerrors

Issue:
Suspicious dvalues for mean of multiple ensembles
### Describe the issue:

When we have an observable that is calculated on multiple ensembles (with similar results and errors) I would expect the the error of the whole mean to decrease with the number of independent ensembles added. This does not seem to be the case. 
Is there a reason, why the error should increase?

### Code example:

```python
def make_test_obs(number_of_ensembles:int):
    samples=[[1.1,0.9,0.9,1.0,1.0,1.1,1.1,0.9]]*number_of_ensembles
    names=[f"ensemble_{i}" for i in range(number_of_ensembles)]
    o=pe.Obs(samples,names)
    o.gm()
    return o
lsp=[1,10,100,500,1000]
plt.plot(lsp,[make_test_obs(n).dvalue for n in lsp],"bo"),
plt.show()
```

### Error message:

```shell

```

### Runtime information:

system    	 Linux
python    	 3.12.7
pyerrors  	 2.13.0
numpy     	 2.2.2
scipy     	 1.15.1
matplotlib	 3.10.0
pandas    	 2.2.3

Hints from issue comments:
@s-kuberski: Hi,
I am not sure that your example does make sense. What kind of observable would be constructed like this? 

If you would have multiple replica (i.e., you would assume that each distribution has the same mean, variance and autocorrelation), then your expectation is fulfilled and the error decreases with $\sqrt{N}$, as you can check by including a ``|`` in the name of the ensembles. This is because we know that the distributions have to match:

```python
def make_test_obs(number_of_ensembles:int):
    samples=[[1.1,0.9,0.9,1.0,1.0,1.1,1.1,0.9]]*number_of_ensembles
    names=[f"ensemble_|{i}" for i in range(number_of_ensembles)]
    o=pe.Obs(samples,names)
    o.gm()
    return o
lsp=[1,10,100,500,1000]
plt.plot(lsp,[make_test_obs(n).dvalue for n in lsp],"bo"),
plt.show()
```

As another way to properly construct the observable when you have ensembles with matching parameters but different autocorrelation functions, you could perform an average. For example, you could weight each ensemble with the number of configurations, if the autocorrelation times are very similar. For the example here, this also fulfills your expectations

``` python
def make_test_obs(number_of_ensembles:int):
    o = sum([pe.Obs([[1.1,0.9,0.9,1.0,1.0,1.1,1.1,0.9]], [f"ensemble_{i}"]) for i in range(number_of_ensembles)]) / number_of_ensembles
    o.gm()
    return o
lsp=[1,10,100,500,1000]
plt.plot(lsp,[make_test_obs(n).dvalue for n in lsp],"bo"),
plt.show()
```
Here, I put in by hand the knowledge, that the observables on the independent ensembles have the same central value.

So I guess the question is if it should be allowed to construct the `Obs` in the way you do it in your example, right?@s-kuberski: To summarize, in your example, you are adding an independent contribution to the error with each additional ensemble. The resulting error of your observable should thus increase as $\sqrt{N}$ because you have $\sigma = \sqrt{\sigma_0^2 + \sigma_1^2 + \dots + \sigma_{N-1}^2}$.@JanNeuendorf: Yes, I did not use replicas on purpose and I agree that it appears to work correctly for them. 
I do not get, why I would expect it to increase with $\sqrt{N}$.
If I set up $N$ ensembles and perform $N$ measurements with similar errors, I would expect $\frac{1}{\sqrt{N}}$ because I have independent samples.  You can argue for this in the following way:
Lets say that my individual ensembles all have a $\sigma$ of around 1 (and they actually come from uncorrelated samples with a nice Gaussian distribution).
So when I ask the program for their individual dvalue, I get a meaningful answer. 
"You will find 68% of values in this range around the mean".
When I now combine 100 of those ensemble, the dvalue would gain a factor of 10. 
So what does this error mean? @JanNeuendorf: > What kind of observable would be constructed like this?

I do not have a real example (someone else told me about weird errors when using multiple ensembles and I just tried to build a minimal example). 
But I can conceive of a scenario where I would use this. 
Let`s say I have two ensembles. They are completely independent and only differ in some parameter that I expect not to matter for what I am calculating. I have samples from both and I want to combine them. I could construct individual obs and average them. This is probably the intended way. But I would have assumed that building a single Obs with all the data in it would have the same effect. 



> So I guess the question is if it should be allowed to construct the Obs in the way you do it in your example, right?

I guess so. But I can construct it this way and the mean value does correspond to the average, so why would I expect the errors to work different from an average?@s-kuberski: Hi,

> I do not get, why I would expect it to increase with $\sqrt{N}$

Think about what happens when we evaluate the contributions to the uncertainty, even in the case where we completely neglect autocorrelation and can just assume independent measurements. For each of the $N$ ensembles, you will get a contribution to the uncertainty of the observable that scales with $\sqrt{n_i}$, which is the number of configurations in this ensemble. So by adding an ensemble with fluctuations, you increase the uncertainty such that
$$\sigma \propto \sqrt{\sum_{i=0}^{N - 1} \sigma_i^2}\,\quad \text{with} \quad \sigma_i \propto \sqrt{n_i}\,.$$

Going to the case of your example, where all ensembles are identical, then $\sigma_i = \sigma^\ast$ and 
$$\sigma = \sqrt{N} \sigma^\ast$$
so you get a trivial relation for the uncertainty.


> If I set up $N$ ensembles [...] When I now combine 100 of those ensemble, the dvalue would gain a factor of 10.

All of this relies on the assumption that you have the same distribution underlying each of the ensembles. How should the code know that this is the case? How should it distinguish this case from the ususal case where two ensembles contribute meaningfully with their own uncertainty? 
The proper way to handle this are replica. However, when the autocorrelation functions differ between ensembles that have the same parameters of the Lagrangian, there is no one-fits-all way to combine the data. The users have to do it themselves. 

> I do not have a real example (someone else told me about weird errors when using multiple ensembles and I just tried to build a minimal example). 

Maybe we should really prevent this way of creating `Obs`. I have to think about it to see if it would have negative implications. We use the feature internally, e.g. here https://github.com/fjosw/pyerrors/blob/7eabd68c5f19ad299eb677314be17985f410ee0a/pyerrors/obs.py#L1347 so we would have to find a way to do it internally without allowing it when calling `Obs()`. The user should really never take primary data from different ensembles to naively construct an `Obs`.

> This is probably the intended way. But I would have assumed that building a single Obs with all the data in it would have the same effect.

The point is also, that any such combination is not free of ambiguities. It is really a step that should be performed by the user.

> But I can construct it this way and the mean value does correspond to the average, so why would I expect the errors to work different from an average?

That's another argument to enforce some kind of barrier to prevent users from doing the wrong thing. When there is an ambiguity that can lead to unexpected behavior, then we have to remove it.
@JanNeuendorf: I think the fundamental problem is that the `.value` of an `Obs` created this way corresponds to the (weighted) mean of all ensembles, but its `.dvalue` does not. This is extremely confusing.
 So I think you are right that this should be prevented. 
I can not come up with something that is not a breaking change. 
- One Idea could be to still allow creating `Obs` like this. After all, there is nothing wrong with that by itself. 
But then there is some flag that does not allow you to call the gamma method on it. 
- Or, even simpler, we prevent initializing an `Obs` with multiple (non-replica) ensembles and have a second method `Obs.ensemble_average(...)` that does a correct average. 
 @s-kuberski: I thought about raising an exception if multiple ensembles are used but the mean values are not set explicitly. Warnings are tricky because they could be raised in internal routines and this should never happen, if the use case is sane.

Raising an exception would still allow us to create new Obs internally, e.g. in 

https://github.com/fjosw/pyerrors/blob/7eabd68c5f19ad299eb677314be17985f410ee0a/pyerrors/obs.py#L1347

In the reweighting routine, 

https://github.com/fjosw/pyerrors/blob/7eabd68c5f19ad299eb677314be17985f410ee0a/pyerrors/obs.py#L1389

the means are not passed explicitly, but I think that in principle there one should not be allowed to reweight observables that are based on several ensembles, don't you think? At least I don't see how reweighting would be defined there, because the gauge average is changed by this procedure and this average should be fixed before one combines ensembles. So having multiple ensembles in this routine should probably raise an exception anyways.

A more difficult case is the `correlate` routine, where a new `Obs` is created in 
https://github.com/fjosw/pyerrors/blob/7eabd68c5f19ad299eb677314be17985f410ee0a/pyerrors/obs.py#L1463-L1469
and in principle it is allowed to have observables that are based on several ensembles. I currently cannot come up with a case where this should be the case, but we would have to be careful here. I have to think about it, maybe also this case boils down to a creation of an `Obs` from primary data on several ensembles (one could have done the correlation on the primary data before creating the observables) and then it should not be allowed. However, @fjosw has written the sentence 
https://github.com/fjosw/pyerrors/blob/7eabd68c5f19ad299eb677314be17985f410ee0a/pyerrors/obs.py#L1445
in the description of the function. I would interpret this as having to have to identical sets of configurations is not strictly necessary (not that one has to have a single ensemble), but maybe I am missing something here.

In `merge_obs`
https://github.com/fjosw/pyerrors/blob/7eabd68c5f19ad299eb677314be17985f410ee0a/pyerrors/obs.py#L1757
we do basically the same thing that would be done when initializing an `Obs`. I thought of this routine to be used to combine data from several replica. If we think that one should not create `Obs` with different ensembles, then an error should also be thrown in `merge_obs`, so this would be fine (one could do it explicitly, here).

The two other cases where some work would have to be done (but can be done, I think) is the JSON input in 
https://github.com/fjosw/pyerrors/blob/7eabd68c5f19ad299eb677314be17985f410ee0a/pyerrors/input/json.py#L139
and
https://github.com/fjosw/pyerrors/blob/7eabd68c5f19ad299eb677314be17985f410ee0a/pyerrors/input/json.py#L334


What do you think? Concerning a routine to merge data from several ensembles: One could come up with something like this. There is the naive average where one weights each ensemble with the number of configs, as an alternative other one could weight with the number of effective configs, taking autocorrelation into account. If you think that this would be a useful routine, we could add it. However, we would have to very specfic about how to use it in which cases. I also realized that the `merge_obs` routine is not perfectly well documented.
@JanNeuendorf: I think that sounds very reasonable.
Maybe there could just be an extra boolean argument to the `__init__()` method that just means "Don't throw an exception for multiple ensembles; I know what I am doing." And then this is used in the internal calls. 

> The two other cases where some work would have to be done (but can be done, I think) is the JSON input in

Ah, the JSON input might be a bit of a problem. Not reading previously valid files is probably a No-Go. 
I think we just have to accept Obs from previous versions as they are. @s-kuberski: I've implemented what I had in mind in #258. For users who don't do anything wrong, nothing should change when we would decide to use this. We can discuss whether this goes in the right direction.  

Solution patch:
diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml
index ff5d8223..a4c27116 100644
--- a/.github/workflows/pytest.yml
+++ b/.github/workflows/pytest.yml
@@ -20,7 +20,9 @@ jobs:
         python-version: ["3.9", "3.10", "3.11", "3.12", "3.13"]
         include:
           - os: macos-latest
-            python-version: "3.10"
+            python-version: "3.12"
+          - os: ubuntu-24.04-arm
+            python-version: "3.12"
 
     steps:
       - name: Checkout source
diff --git a/pyerrors/fits.py b/pyerrors/fits.py
index 89f73e5b..675bdca6 100644
--- a/pyerrors/fits.py
+++ b/pyerrors/fits.py
@@ -365,6 +365,8 @@ def chisqfunc_uncorr(p):
             if (chol_inv[1] != key_ls):
                 raise ValueError('The keys of inverse covariance matrix are not the same or do not appear in the same order as the x and y values.')
             chol_inv = chol_inv[0]
+            if np.any(np.diag(chol_inv) <= 0) or (not np.all(chol_inv == np.tril(chol_inv))):
+                raise ValueError('The inverse covariance matrix inv_chol_cov_matrix[0] has to be a lower triangular matrix constructed from a Cholesky decomposition.')
         else:
             corr = covariance(y_all, correlation=True, **kwargs)
             inverrdiag = np.diag(1 / np.asarray(dy_f))
diff --git a/pyerrors/input/dobs.py b/pyerrors/input/dobs.py
index aea9b7a9..6907ec3c 100644
--- a/pyerrors/input/dobs.py
+++ b/pyerrors/input/dobs.py
@@ -529,7 +529,8 @@ def import_dobs_string(content, full_output=False, separator_insertion=True):
                 deltas.append(repdeltas)
                 idl.append(repidl)
 
-        res.append(Obs(deltas, obs_names, idl=idl))
+        obsmeans = [np.average(deltas[j]) for j in range(len(deltas))]
+        res.append(Obs([np.array(deltas[j]) - obsmeans[j] for j in range(len(obsmeans))], obs_names, idl=idl, means=obsmeans))
         res[-1]._value = mean[i]
     _check(len(e_names) == ne)
 
diff --git a/pyerrors/input/json.py b/pyerrors/input/json.py
index ca3fb0d2..a0d3304a 100644
--- a/pyerrors/input/json.py
+++ b/pyerrors/input/json.py
@@ -133,10 +133,11 @@ def _nan_Obs_like(obs):
         names = []
         idl = []
         for key, value in obs.idl.items():
-            samples.append([np.nan] * len(value))
+            samples.append(np.array([np.nan] * len(value)))
             names.append(key)
             idl.append(value)
-        my_obs = Obs(samples, names, idl)
+        my_obs = Obs(samples, names, idl, means=[np.nan for n in names])
+        my_obs._value = np.nan
         my_obs._covobs = obs._covobs
         for name in obs._covobs:
             my_obs.names.append(name)
@@ -331,7 +332,8 @@ def get_Obs_from_dict(o):
         cd = _gen_covobsd_from_cdatad(o.get('cdata', {}))
 
         if od:
-            ret = Obs([[ddi[0] + values[0] for ddi in di] for di in od['deltas']], od['names'], idl=od['idl'])
+            r_offsets = [np.average([ddi[0] for ddi in di]) for di in od['deltas']]
+            ret = Obs([np.array([ddi[0] for ddi in od['deltas'][i]]) - r_offsets[i] for i in range(len(od['deltas']))], od['names'], idl=od['idl'], means=[ro + values[0] for ro in r_offsets])
             ret._value = values[0]
         else:
             ret = Obs([], [], means=[])
@@ -356,7 +358,8 @@ def get_List_from_dict(o):
         taglist = o.get('tag', layout * [None])
         for i in range(layout):
             if od:
-                ret.append(Obs([list(di[:, i] + values[i]) for di in od['deltas']], od['names'], idl=od['idl']))
+                r_offsets = np.array([np.average(di[:, i]) for di in od['deltas']])
+                ret.append(Obs([od['deltas'][j][:, i] - r_offsets[j] for j in range(len(od['deltas']))], od['names'], idl=od['idl'], means=[ro + values[i] for ro in r_offsets]))
                 ret[-1]._value = values[i]
             else:
                 ret.append(Obs([], [], means=[]))
@@ -383,7 +386,8 @@ def get_Array_from_dict(o):
         taglist = o.get('tag', N * [None])
         for i in range(N):
             if od:
-                ret.append(Obs([di[:, i] + values[i] for di in od['deltas']], od['names'], idl=od['idl']))
+                r_offsets = np.array([np.average(di[:, i]) for di in od['deltas']])
+                ret.append(Obs([od['deltas'][j][:, i] - r_offsets[j] for j in range(len(od['deltas']))], od['names'], idl=od['idl'], means=[ro + values[i] for ro in r_offsets]))
                 ret[-1]._value = values[i]
             else:
                 ret.append(Obs([], [], means=[]))
diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 623c37fd..87591cd9 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -82,6 +82,8 @@ def __init__(self, samples, names, idl=None, **kwargs):
                     raise ValueError('Names are not unique.')
                 if not all(isinstance(x, str) for x in names):
                     raise TypeError('All names have to be strings.')
+                if len(set([o.split('|')[0] for o in names])) > 1:
+                    raise ValueError('Cannot initialize Obs based on multiple ensembles. Please average separate Obs from each ensemble.')
             else:
                 if not isinstance(names[0], str):
                     raise TypeError('All names have to be strings.')
@@ -1407,6 +1409,8 @@ def reweight(weight, obs, **kwargs):
             raise ValueError('Error: Not possible to reweight an Obs that contains covobs!')
         if not set(obs[i].names).issubset(weight.names):
             raise ValueError('Error: Ensembles do not fit')
+        if len(obs[i].mc_names) > 1 or len(weight.mc_names) > 1:
+            raise ValueError('Error: Cannot reweight an Obs that contains multiple ensembles.')
         for name in obs[i].names:
             if not set(obs[i].idl[name]).issubset(weight.idl[name]):
                 raise ValueError('obs[%d] has to be defined on a subset of the configs in weight.idl[%s]!' % (i, name))
@@ -1442,9 +1446,12 @@ def correlate(obs_a, obs_b):
     -----
     Keep in mind to only correlate primary observables which have not been reweighted
     yet. The reweighting has to be applied after correlating the observables.
-    Currently only works if ensembles are identical (this is not strictly necessary).
+    Only works if a single ensemble is present in the Obs.
+    Currently only works if ensemble content is identical (this is not strictly necessary).
     """
 
+    if len(obs_a.mc_names) > 1 or len(obs_b.mc_names) > 1:
+        raise ValueError('Error: Cannot correlate Obs that contain multiple ensembles.')
     if sorted(obs_a.names) != sorted(obs_b.names):
         raise ValueError(f"Ensembles do not fit {set(sorted(obs_a.names)) ^ set(sorted(obs_b.names))}")
     if len(obs_a.cov_names) or len(obs_b.cov_names):
@@ -1755,7 +1762,11 @@ def import_bootstrap(boots, name, random_numbers):
 
 
 def merge_obs(list_of_obs):
-    """Combine all observables in list_of_obs into one new observable
+    """Combine all observables in list_of_obs into one new observable.
+    This allows to merge Obs that have been computed on multiple replica
+    of the same ensemble.
+    If you like to merge Obs that are based on several ensembles, please
+    average them yourself.
 
     Parameters
     ----------


Test patch:
diff --git a/tests/fits_test.py b/tests/fits_test.py
index 2eeb6a49..283ff6a2 100644
--- a/tests/fits_test.py
+++ b/tests/fits_test.py
@@ -223,6 +223,9 @@ def fitf(p, x):
             diff_inv_cov_combined_fit.gamma_method()
             assert(diff_inv_cov_combined_fit.is_zero(atol=1e-12))
 
+        with pytest.raises(ValueError):
+            pe.least_squares(x_dict, data_dict, fitf_dict,  correlated_fit = True, inv_chol_cov_matrix = [corr,chol_inv_keys_combined_fit])
+
 def test_least_squares_invalid_inv_cov_matrix_input():
     xvals = []
     yvals = []
diff --git a/tests/json_io_test.py b/tests/json_io_test.py
index dafaaa41..a9263691 100644
--- a/tests/json_io_test.py
+++ b/tests/json_io_test.py
@@ -12,7 +12,7 @@ def test_jsonio():
     o = pe.pseudo_Obs(1.0, .2, 'one')
     o2 = pe.pseudo_Obs(0.5, .1, 'two|r1')
     o3 = pe.pseudo_Obs(0.5, .1, 'two|r2')
-    o4 = pe.merge_obs([o2, o3])
+    o4 = pe.merge_obs([o2, o3, pe.pseudo_Obs(0.5, .1, 'two|r3', samples=3221)])
     otag = 'This has been merged!'
     o4.tag = otag
     do = o - .2 * o4
@@ -101,8 +101,8 @@ def test_json_string_reconstruction():
 
 
 def test_json_corr_io():
-    my_list = [pe.Obs([np.random.normal(1.0, 0.1, 100)], ['ens1']) for o in range(8)]
-    rw_list = pe.reweight(pe.Obs([np.random.normal(1.0, 0.1, 100)], ['ens1']), my_list)
+    my_list = [pe.Obs([np.random.normal(1.0, 0.1, 100), np.random.normal(1.0, 0.1, 321)], ['ens1|r1', 'ens1|r2'], idl=[range(1, 201, 2), range(321)]) for o in range(8)]
+    rw_list = pe.reweight(pe.Obs([np.random.normal(1.0, 0.1, 100), np.random.normal(1.0, 0.1, 321)], ['ens1|r1', 'ens1|r2'], idl=[range(1, 201, 2), range(321)]), my_list)
 
     for obs_list in [my_list, rw_list]:
         for tag in [None, "test"]:
@@ -111,40 +111,51 @@ def test_json_corr_io():
                 for corr_tag in [None, 'my_Corr_tag']:
                     for prange in [None, [3, 6]]:
                         for gap in [False, True]:
-                            my_corr = pe.Corr(obs_list, padding=[pad, pad], prange=prange)
-                            my_corr.tag = corr_tag
-                            if gap:
-                                my_corr.content[4] = None
-                            pe.input.json.dump_to_json(my_corr, 'corr')
-                            recover = pe.input.json.load_json('corr')
-                            os.remove('corr.json.gz')
-                            assert np.all([o.is_zero() for o in [x for x in (my_corr - recover) if x is not None]])
-                            for index, entry in enumerate(my_corr):
-                                if entry is None:
-                                    assert recover[index] is None
-                            assert my_corr.tag == recover.tag
-                            assert my_corr.prange == recover.prange
-                            assert my_corr.reweighted == recover.reweighted
+                            for mult in [1., pe.cov_Obs([12.22, 1.21], [.212**2, .11**2], 'renorm')[0]]:
+                                my_corr = mult * pe.Corr(obs_list, padding=[pad, pad], prange=prange)
+                                my_corr.tag = corr_tag
+                                if gap:
+                                    my_corr.content[4] = None
+                                pe.input.json.dump_to_json(my_corr, 'corr')
+                                recover = pe.input.json.load_json('corr')
+                                os.remove('corr.json.gz')
+                                assert np.all([o.is_zero() for o in [x for x in (my_corr - recover) if x is not None]])
+                                for index, entry in enumerate(my_corr):
+                                    if entry is None:
+                                        assert recover[index] is None
+                                assert my_corr.tag == recover.tag
+                                assert my_corr.prange == recover.prange
+                                assert my_corr.reweighted == recover.reweighted
 
 
 def test_json_corr_2d_io():
-    obs_list = [np.array([[pe.pseudo_Obs(1.0 + i, 0.1 * i, 'test'), pe.pseudo_Obs(0.0, 0.1 * i, 'test')], [pe.pseudo_Obs(0.0, 0.1 * i, 'test'), pe.pseudo_Obs(1.0 + i, 0.1 * i, 'test')]]) for i in range(4)]
+    obs_list = [np.array([
+        [
+         pe.merge_obs([pe.pseudo_Obs(1.0 + i, 0.1 * i, 'test|r2'), pe.pseudo_Obs(1.0 + i, 0.1 * i, 'test|r1', samples=321)]), 
+         pe.merge_obs([pe.pseudo_Obs(0.0, 0.1 * i, 'test|r2'), pe.pseudo_Obs(0.0, 0.1 * i, 'test|r1', samples=321)]),
+        ], 
+        [
+         pe.merge_obs([pe.pseudo_Obs(0.0, 0.1 * i, 'test|r2'), pe.pseudo_Obs(0.0, 0.1 * i, 'test|r1', samples=321),]),
+         pe.merge_obs([pe.pseudo_Obs(1.0 + i, 0.1 * i, 'test|r2'), pe.pseudo_Obs(1.0 + i, 0.1 * i, 'test|r1', samples=321)]),
+        ],
+        ]) for i in range(4)]
 
     for tag in [None, "test"]:
         obs_list[3][0, 1].tag = tag
         for padding in [0, 1]:
             for prange in [None, [3, 6]]:
-                my_corr = pe.Corr(obs_list, padding=[padding, padding], prange=prange)
-                my_corr.tag = tag
-                pe.input.json.dump_to_json(my_corr, 'corr')
-                recover = pe.input.json.load_json('corr')
-                os.remove('corr.json.gz')
-                assert np.all([np.all([o.is_zero() for o in q]) for q in [x.ravel() for x in (my_corr - recover) if x is not None]])
-                for index, entry in enumerate(my_corr):
-                    if entry is None:
-                        assert recover[index] is None
-                assert my_corr.tag == recover.tag
-                assert my_corr.prange == recover.prange
+                for mult in [1., pe.cov_Obs([12.22, 1.21], [.212**2, .11**2], 'renorm')[0]]:
+                    my_corr = mult * pe.Corr(obs_list, padding=[padding, padding], prange=prange)
+                    my_corr.tag = tag
+                    pe.input.json.dump_to_json(my_corr, 'corr')
+                    recover = pe.input.json.load_json('corr')
+                    os.remove('corr.json.gz')
+                    assert np.all([np.all([o.is_zero() for o in q]) for q in [x.ravel() for x in (my_corr - recover) if x is not None]])
+                    for index, entry in enumerate(my_corr):
+                        if entry is None:
+                            assert recover[index] is None
+                    assert my_corr.tag == recover.tag
+                    assert my_corr.prange == recover.prange
 
 
 def test_json_dict_io():
@@ -211,6 +222,7 @@ def list_check_obs(l1, l2):
             'd': pe.pseudo_Obs(.01, .001, 'testd', samples=10) * pe.cov_Obs(1, .01, 'cov1'),
             'se': None,
             'sf': 1.2,
+            'k': pe.cov_Obs(.1, .001**2, 'cov') * pe.merge_obs([pe.pseudo_Obs(1.0, 0.1, 'test|r2'), pe.pseudo_Obs(1.0, 0.1, 'test|r1', samples=321)]),
         }
     }
 
@@ -314,7 +326,7 @@ def test_dobsio():
 
     o2 = pe.pseudo_Obs(0.5, .1, 'two|r1')
     o3 = pe.pseudo_Obs(0.5, .1, 'two|r2')
-    o4 = pe.merge_obs([o2, o3])
+    o4 = pe.merge_obs([o2, o3, pe.pseudo_Obs(0.5, .1, 'two|r3', samples=3221)])
     otag = 'This has been merged!'
     o4.tag = otag
     do = o - .2 * o4
@@ -328,7 +340,7 @@ def test_dobsio():
     o5 /= co2[0]
     o5.tag = 2 * otag
 
-    tt1 = pe.Obs([np.random.rand(100), np.random.rand(100)], ['t|r1', 't|r2'], idl=[range(2, 202, 2), range(22, 222, 2)])
+    tt1 = pe.Obs([np.random.rand(100), np.random.rand(102)], ['t|r1', 't|r2'], idl=[range(2, 202, 2), range(22, 226, 2)])
     tt3 = pe.Obs([np.random.rand(102)], ['qe|r1'])
 
     tt = tt1 + tt3
@@ -337,7 +349,7 @@ def test_dobsio():
 
     tt4 = pe.Obs([np.random.rand(100), np.random.rand(100)], ['t|r1', 't|r2'], idl=[range(1, 101, 1), range(2, 202, 2)])
 
-    ol = [o2, o3, o4, do, o5, tt, tt4, np.log(tt4 / o5**2), np.exp(o5 + np.log(co3 / tt3 + o4) / tt)]
+    ol = [o2, o3, o4, do, o5, tt, tt4, np.log(tt4 / o5**2), np.exp(o5 + np.log(co3 / tt3 + o4) / tt), o4.reweight(o4)]
     print(ol)
     fname = 'test_rw'
 
@@ -362,9 +374,12 @@ def test_dobsio():
 
 
 def test_reconstruct_non_linear_r_obs(tmp_path):
-    to = pe.Obs([np.random.rand(500), np.random.rand(500), np.random.rand(111)],
-                ["e|r1", "e|r2", "my_new_ensemble_54^£$|8'[@124435%6^7&()~#"],
-                idl=[range(1, 501), range(0, 500), range(1, 999, 9)])
+    to = (
+        pe.Obs([np.random.rand(500), np.random.rand(1200)],
+                ["e|r1", "e|r2", ],
+                idl=[range(1, 501), range(0, 1200)])
+        + pe.Obs([np.random.rand(111)], ["my_new_ensemble_54^£$|8'[@124435%6^7&()~#"], idl=[range(1, 999, 9)])
+        )
     to = np.log(to ** 2) / to
     to.dump((tmp_path / "test_equality").as_posix())
     ro = pe.input.json.load_json((tmp_path / "test_equality").as_posix())
@@ -372,9 +387,12 @@ def test_reconstruct_non_linear_r_obs(tmp_path):
 
 
 def test_reconstruct_non_linear_r_obs_list(tmp_path):
-    to = pe.Obs([np.random.rand(500), np.random.rand(500), np.random.rand(111)],
-                ["e|r1", "e|r2", "my_new_ensemble_54^£$|8'[@124435%6^7&()~#"],
-                idl=[range(1, 501), range(0, 500), range(1, 999, 9)])
+    to = (
+        pe.Obs([np.random.rand(500), np.random.rand(1200)],
+                ["e|r1", "e|r2", ],
+                idl=[range(1, 501), range(0, 1200)])
+        + pe.Obs([np.random.rand(111)], ["my_new_ensemble_54^£$|8'[@124435%6^7&()~#"], idl=[range(1, 999, 9)])
+        )
     to = np.log(to ** 2) / to
     for to_list in [[to, to, to], np.array([to, to, to])]:
         pe.input.json.dump_to_json(to_list, (tmp_path / "test_equality_list").as_posix())
diff --git a/tests/linalg_test.py b/tests/linalg_test.py
index 4fb952d3..9323cfcf 100644
--- a/tests/linalg_test.py
+++ b/tests/linalg_test.py
@@ -34,7 +34,7 @@ def test_matmul():
             my_list = []
             length = 100 + np.random.randint(200)
             for i in range(dim ** 2):
-                my_list.append(pe.Obs([np.random.rand(length), np.random.rand(length + 1)], ['t1', 't2']))
+                my_list.append(pe.Obs([np.random.rand(length)], ['t1']) + pe.Obs([np.random.rand(length + 1)], ['t2']))
             my_array = const * np.array(my_list).reshape((dim, dim))
             tt = pe.linalg.matmul(my_array, my_array) - my_array @ my_array
             for t, e in np.ndenumerate(tt):
@@ -43,8 +43,8 @@ def test_matmul():
             my_list = []
             length = 100 + np.random.randint(200)
             for i in range(dim ** 2):
-                my_list.append(pe.CObs(pe.Obs([np.random.rand(length), np.random.rand(length + 1)], ['t1', 't2']),
-                                       pe.Obs([np.random.rand(length), np.random.rand(length + 1)], ['t1', 't2'])))
+                my_list.append(pe.CObs(pe.Obs([np.random.rand(length)], ['t1']) + pe.Obs([np.random.rand(length + 1)], ['t2']),
+                                       pe.Obs([np.random.rand(length)], ['t1']) + pe.Obs([np.random.rand(length + 1)], ['t2'])))
             my_array = np.array(my_list).reshape((dim, dim)) * const
             tt = pe.linalg.matmul(my_array, my_array) - my_array @ my_array
             for t, e in np.ndenumerate(tt):
@@ -151,7 +151,7 @@ def test_multi_dot():
         my_list = []
         length = 1000 + np.random.randint(200)
         for i in range(dim ** 2):
-            my_list.append(pe.Obs([np.random.rand(length), np.random.rand(length + 1)], ['t1', 't2']))
+            my_list.append(pe.Obs([np.random.rand(length)], ['t1']) + pe.Obs([np.random.rand(length + 1)], ['t2']))
         my_array = pe.cov_Obs(1.0, 0.002, 'cov') * np.array(my_list).reshape((dim, dim))
         tt = pe.linalg.matmul(my_array, my_array, my_array, my_array) - my_array @ my_array @ my_array @ my_array
         for t, e in np.ndenumerate(tt):
@@ -160,8 +160,8 @@ def test_multi_dot():
         my_list = []
         length = 1000 + np.random.randint(200)
         for i in range(dim ** 2):
-            my_list.append(pe.CObs(pe.Obs([np.random.rand(length), np.random.rand(length + 1)], ['t1', 't2']),
-                                   pe.Obs([np.random.rand(length), np.random.rand(length + 1)], ['t1', 't2'])))
+            my_list.append(pe.CObs(pe.Obs([np.random.rand(length)], ['t1']) + pe.Obs([np.random.rand(length + 1)], ['t2']),
+                                   pe.Obs([np.random.rand(length)], ['t1']) + pe.Obs([np.random.rand(length + 1)], ['t2'])))
         my_array = np.array(my_list).reshape((dim, dim)) * pe.cov_Obs(1.0, 0.002, 'cov')
         tt = pe.linalg.matmul(my_array, my_array, my_array, my_array) - my_array @ my_array @ my_array @ my_array
         for t, e in np.ndenumerate(tt):
@@ -209,7 +209,7 @@ def test_irregular_matrix_inverse():
     for idl in [range(8, 508, 10), range(250, 273), [2, 8, 19, 20, 78, 99, 828, 10548979]]:
         irregular_array = []
         for i in range(dim ** 2):
-            irregular_array.append(pe.Obs([np.random.normal(1.1, 0.2, len(idl)), np.random.normal(0.25, 0.1, 10)], ['ens1', 'ens2'], idl=[idl, range(1, 11)]))
+            irregular_array.append(pe.Obs([np.random.normal(1.1, 0.2, len(idl))], ['ens1'], idl=[idl]) + pe.Obs([np.random.normal(0.25, 0.1, 10)], ['ens2'], idl=[range(1, 11)]))
         irregular_matrix = np.array(irregular_array).reshape((dim, dim)) * pe.cov_Obs(1.0, 0.002, 'cov') * pe.pseudo_Obs(1.0, 0.002, 'ens2|r23')
 
         invertible_irregular_matrix = np.identity(dim) + irregular_matrix @ irregular_matrix.T
diff --git a/tests/obs_test.py b/tests/obs_test.py
index 2c642ad4..546a4bfd 100644
--- a/tests/obs_test.py
+++ b/tests/obs_test.py
@@ -333,7 +333,7 @@ def test_derived_observables():
 
 def test_multi_ens():
     names = ['A0', 'A1|r001', 'A1|r002']
-    test_obs = pe.Obs([np.random.rand(50), np.random.rand(50), np.random.rand(50)], names)
+    test_obs = pe.Obs([np.random.rand(50)], names[:1]) + pe.Obs([np.random.rand(50), np.random.rand(50)], names[1:])
     assert test_obs.e_names == ['A0', 'A1']
     assert test_obs.e_content['A0'] == ['A0']
     assert test_obs.e_content['A1'] == ['A1|r001', 'A1|r002']
@@ -345,6 +345,9 @@ def test_multi_ens():
         ensembles.append(str(i))
     assert my_sum.e_names == sorted(ensembles)
 
+    with pytest.raises(ValueError):
+        test_obs = pe.Obs([np.random.rand(50), np.random.rand(50), np.random.rand(50)], names)
+
 
 def test_multi_ens2():
     names = ['ens', 'e', 'en', 'e|r010', 'E|er', 'ens|', 'Ens|34', 'ens|r548984654ez4e3t34terh']
@@ -498,18 +501,25 @@ def test_reweighting():
     with pytest.raises(ValueError):
         pe.reweight(my_irregular_obs, [my_obs])
 
+    my_merged_obs = my_obs + pe.Obs([np.random.rand(1000)], ['q'])
+    with pytest.raises(ValueError):
+        pe.reweight(my_merged_obs, [my_merged_obs])
+
 
 def test_merge_obs():
-    my_obs1 = pe.Obs([np.random.rand(100)], ['t'])
-    my_obs2 = pe.Obs([np.random.rand(100)], ['q'], idl=[range(1, 200, 2)])
+    my_obs1 = pe.Obs([np.random.normal(1, .1, 100)], ['t|1'])
+    my_obs2 = pe.Obs([np.random.normal(1, .1, 100)], ['t|2'], idl=[range(1, 200, 2)])
     merged = pe.merge_obs([my_obs1, my_obs2])
-    diff = merged - my_obs2 - my_obs1
-    assert diff == -(my_obs1.value + my_obs2.value) / 2
+    diff = merged - (my_obs2 + my_obs1) / 2
+    assert np.isclose(0, diff.value, atol=1e-16)
     with pytest.raises(ValueError):
         pe.merge_obs([my_obs1, my_obs1])
     my_covobs = pe.cov_Obs(1.0, 0.003, 'cov')
     with pytest.raises(ValueError):
         pe.merge_obs([my_obs1, my_covobs])
+    my_obs3 = pe.Obs([np.random.rand(100)], ['q|2'], idl=[range(1, 200, 2)])
+    with pytest.raises(ValueError):
+        pe.merge_obs([my_obs1, my_obs3])
 
 
 
@@ -542,6 +552,9 @@ def test_correlate():
     my_obs6 = pe.Obs([np.random.rand(100)], ['t'], idl=[range(5, 505, 5)])
     corr3 = pe.correlate(my_obs5, my_obs6)
     assert my_obs5.idl == corr3.idl
+    my_obs7 = pe.Obs([np.random.rand(99)], ['q'])
+    with pytest.raises(ValueError):
+        pe.correlate(my_obs1, my_obs7)
 
     my_new_obs = pe.Obs([np.random.rand(100)], ['q3'])
     with pytest.raises(ValueError):
@@ -681,14 +694,14 @@ def test_gamma_method_irregular():
     assert (a.dvalue - 5 * a.ddvalue < expe and expe < a.dvalue + 5 * a.ddvalue)
 
     arr2 = np.random.normal(1, .2, size=N)
-    afull = pe.Obs([arr, arr2], ['a1', 'a2'])
+    afull = pe.Obs([arr], ['a1']) + pe.Obs([arr2], ['a2'])
 
     configs = np.ones_like(arr2)
     for i in np.random.uniform(0, len(arr2), size=int(.8*N)):
         configs[int(i)] = 0
     zero_arr2 = [arr2[i] for i in range(len(arr2)) if not configs[i] == 0]
     idx2 = [i + 1 for i in range(len(configs)) if configs[i] == 1]
-    a = pe.Obs([zero_arr, zero_arr2], ['a1', 'a2'], idl=[idx, idx2])
+    a = pe.Obs([zero_arr], ['a1'], idl=[idx]) + pe.Obs([zero_arr2], ['a2'], idl=[idx2])
 
     afull.gamma_method()
     a.gamma_method()
@@ -1022,7 +1035,7 @@ def test_correlation_intersection_of_idls():
 
 
 def test_covariance_non_identical_objects():
-    obs1 = pe.Obs([np.random.normal(1.0, 0.1, 1000), np.random.normal(1.0, 0.1, 1000), np.random.normal(1.0, 0.1, 732)], ["ens|r1", "ens|r2", "ens2"])
+    obs1 = pe.Obs([np.random.normal(1.0, 0.1, 1000), np.random.normal(1.0, 0.1, 1000)], ["ens|r1", "ens|r2"]) + pe.Obs([np.random.normal(1.0, 0.1, 732)], ['ens2'])
     obs1.gamma_method()
     obs2 = obs1 + 1e-18
     obs2.gamma_method()
@@ -1106,6 +1119,9 @@ def test_reweight_method():
     obs1 = pe.pseudo_Obs(0.2, 0.01, 'test')
     rw = pe.pseudo_Obs(0.999, 0.001, 'test')
     assert obs1.reweight(rw) == pe.reweight(rw, [obs1])[0]
+    rw2 = pe.pseudo_Obs(0.999, 0.001, 'test2')
+    with pytest.raises(ValueError):
+        obs1.reweight(rw2)
 
 
 def test_jackknife():