Merge pull request #48 from JanNeuendorf/develop

Changes in the correlator init, GEVP and docstrings & tag format includes dict and prange

pyerrors/correlators.py

@@ -28,8 +28,8 @@ class Corr:
 
         Parameters
         ----------
-        data_input : list
+        data_input : list or array
-            list of Obs or list of arrays of Obs.
+            list of Obs or list of arrays of Obs or array of Corrs
         padding : list, optional
             List with two entries where the first labels the padding
             at the front of the correlator and the second the padding
@@ -39,25 +39,56 @@ class Corr:
             region indentified for this correlator.
         """
 
-        if not isinstance(data_input, list):
+        if isinstance(data_input, np.ndarray):  # Input is an array of Corrs
-            raise TypeError('Corr__init__ expects a list of timeslices.')
 
-        if all([(isinstance(item, Obs) or isinstance(item, CObs)) or item is None for item in data_input]):
+            # This only works, if the array fulfills the conditions below
-            _assert_equal_properties([o for o in data_input if o is not None])
+            if not len(data_input.shape) == 2 and data_input.shape[0] == data_input.shape[1]:
-            self.content = [np.asarray([item]) if item is not None else None for item in data_input]
+                raise Exception("Incompatible array shape")
-            self.N = 1
+            if not all([isinstance(item, Corr) for item in data_input.flatten()]):
+                raise Exception("If the input is an array, its elements must be of type pe.Corr")
+            if not all([item.N == 1 for item in data_input.flatten()]):
+                raise Exception("Can only construct matrix correlator from single valued correlators")
+            if not len(set([item.T for item in data_input.flatten()])) == 1:
+                raise Exception("All input Correlators must be defined over the same timeslices.")
 
-        elif all([isinstance(item, np.ndarray) or item is None for item in data_input]) and any([isinstance(item, np.ndarray) for item in data_input]):
+            T = data_input[0, 0].T
-            self.content = data_input
+            N = data_input.shape[0]
+            input_as_list = []
+            for t in range(T):
+                if any([(item.content[t][0] is None) for item in data_input.flatten()]):
+                    if not all([(item.content[t][0] is None) for item in data_input.flatten()]):
+                        warnings.warn("Input ill-defined at different timeslices. Conversion leads to data loss!", RuntimeWarning)
+                    input_as_list.append(None)
+                else:
+                    array_at_timeslace = np.empty([N, N], dtype="object")
+                    for i in range(N):
+                        for j in range(N):
+                            array_at_timeslace[i, j] = data_input[i, j][t]
+                    input_as_list.append(array_at_timeslace)
+            data_input = input_as_list
 
-            noNull = [a for a in self.content if not (a is None)]  # To check if the matrices are correct for all undefined elements
+        if isinstance(data_input, list):
-            self.N = noNull[0].shape[0]
+
-            if self.N > 1 and noNull[0].shape[0] != noNull[0].shape[1]:
+            if all([(isinstance(item, Obs) or isinstance(item, CObs)) for item in data_input]):
-                raise Exception("Smearing matrices are not NxN")
+                _assert_equal_properties(data_input)
-            if (not all([item.shape == noNull[0].shape for item in noNull])):
+                self.content = [np.asarray([item]) for item in data_input]
-                raise Exception("Items in data_input are not of identical shape." + str(noNull))
+            if all([(isinstance(item, Obs) or isinstance(item, CObs)) or item is None for item in data_input]):
+                _assert_equal_properties([o for o in data_input if o is not None])
+                self.content = [np.asarray([item]) if item is not None else None for item in data_input]
+                self.N = 1
+
+            elif all([isinstance(item, np.ndarray) or item is None for item in data_input]) and any([isinstance(item, np.ndarray) for item in data_input]):
+                self.content = data_input
+                noNull = [a for a in self.content if not (a is None)]  # To check if the matrices are correct for all undefined elements
+                self.N = noNull[0].shape[0]
+                if self.N > 1 and noNull[0].shape[0] != noNull[0].shape[1]:
+                    raise Exception("Smearing matrices are not NxN")
+                if (not all([item.shape == noNull[0].shape for item in noNull])):
+                    raise Exception("Items in data_input are not of identical shape." + str(noNull))
+            else:
+                raise Exception("data_input contains item of wrong type")
         else:
-            raise Exception("data_input contains item of wrong type")
+            raise Exception("Data input was not given as list or correct array")
 
         self.tag = None
 
@@ -214,8 +245,27 @@ class Corr:
     # There are two ways, the GEVP metod can be called.
     # 1. return_list=False will return a single eigenvector, normalized according to V*C(t_0)*V=1
     # 2. return_list=True  will return a new eigenvector for every timeslice. The time t_s is used to order the vectors according to. arXiv:2004.10472 [hep-lat]
-    def GEVP(self, t0, ts, state=0, sorting="Eigenvalue", return_list=False):
+    def GEVP(self, t0, ts=None, state=0, sorted_list=None):
-        if not return_list:
+        """Solve the general eigenvalue problem on the current correlator
+
+        Parameters
+        ----------
+        t0 : int
+            The time t0 for G(t)v= lambda G(t_0)v
+        ts : int
+            fixed time G(t_s)v= lambda G(t_0)v  if return_list=False
+            If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
+        state : int
+            The state one is interested in ordered by energy. The lowest state is zero.
+        sorted list : string
+            if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned.
+             "Eigenvalue"  -  The eigenvector is chosen according to which einvenvalue it belongs individually on every timeslice.
+             "Eigenvector" -  Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state.
+                                 The referense state is identified by its eigenvalue at t=ts
+        """
+        if sorted_list is None:
+            if (ts is None):
+                raise Exception("ts is required if return_list=False")
             if (self.content[t0] is None) or (self.content[ts] is None):
                 raise Exception("Corr not defined at t0/ts")
             G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double")
@@ -227,7 +277,8 @@ class Corr:
             sp_vecs = _GEVP_solver(Gt, G0)
             sp_vec = sp_vecs[state]
             return sp_vec
-        if return_list:
+        else:
+
             all_vecs = []
             for t in range(self.T):
                 try:
@@ -238,14 +289,16 @@ class Corr:
                             Gt[i, j] = self.content[t][i, j].value
 
                     sp_vecs = _GEVP_solver(Gt, G0)
-                    if sorting == "Eigenvalue":
+                    if sorted_list == "Eigenvalue":
                         sp_vec = sp_vecs[state]
                         all_vecs.append(sp_vec)
                     else:
                         all_vecs.append(sp_vecs)
                 except Exception:
                     all_vecs.append(None)
-            if sorting == "Eigenvector":
+            if sorted_list == "Eigenvector":
+                if (ts is None):
+                    raise Exception("ts is required for the Eigenvector sorting method.")
                 all_vecs = _sort_vectors(all_vecs, ts)
                 all_vecs = [a[state] for a in all_vecs]
 

pyerrors/input/json.py

@@ -186,6 +186,11 @@ def create_json_string(ol, description='', indent=1):
             dat['tag'].append(corr_meta_data)
         else:
             dat['tag'] = [corr_meta_data]
+        taglist = dat['tag']
+        dat['tag'] = {}  # tag is now a dictionary, that contains the previous taglist in the key "tag"
+        dat['tag']['tag'] = taglist
+        if my_corr.prange is not None:
+            dat['tag']['prange'] = my_corr.prange
         return dat
 
     if not isinstance(ol, list):
@@ -395,7 +400,19 @@ def import_json_string(json_string, verbose=True, full_output=False):
         return np.reshape(ret, layout)
 
     def get_Corr_from_dict(o):
-        taglist = o.get('tag')
+        if isinstance(o.get('tag'), list):  # supports the old way
+            taglist = o.get('tag')  # This had to be modified to get the taglist from the dictionary
+            temp_prange = None
+        elif isinstance(o.get('tag'), dict):
+            tagdic = o.get('tag')
+            taglist = tagdic['tag']
+            if 'prange' in tagdic:
+                temp_prange = tagdic['prange']
+            else:
+                temp_prange = None
+        else:
+            raise Exception("The tag is not a list or dict")
+
         corr_tag = taglist[-1]
         tmp_o = o
         tmp_o['tag'] = taglist[:-1]
@@ -405,6 +422,8 @@ def import_json_string(json_string, verbose=True, full_output=False):
         my_corr = Corr([None if np.isnan(o.ravel()[0].value) else o for o in list(dat)])
         if corr_tag != 'None':
             my_corr.tag = corr_tag
+
+        my_corr.prange = temp_prange
         return my_corr
 
     json_dict = json.loads(json_string)