allow error for each component of data vector to be specified (#33)

lompe/model/data.py

@@ -124,8 +124,12 @@ def __init__(self, values, coordinates = None, LOS = None, components = 'all', d
             magnetometer data and iweight=1.0 for ionospheric convection measurements. Keep in 
             mind that this weight is directly applied to the a priori inverse data covariance matrix, 
             so the data error is effectively increased by a factor of 1/sqrt(iweight).
-        error: array of same length as values, or float, optional
+        error: either float, or 1D array of same length as values, or 2D array
+                (shape 2,N or 3,N), optional
             Measurement error. Used to calculate the data covariance matrix. Use SI units.
+            If 2D array, it means that the error of each component (east, north, up) is
+            specified separately. If 1D array, the same value is used for each component.
+            If float, the same value is used for all observations and components.
 
         """
 
@@ -214,8 +218,11 @@ def subset(self, indices):
         """
 
         self.values  = np.array(self.values , ndmin = 2)[:, indices].squeeze()
-        self.error = self.error[indices]
-
+        if np.array(self.error, ndmin=2).shape[0] >= 2: # errors for each component
+            self.error = self.error[:,indices]
+        else:
+            self.error = self.error[indices]
+
         for key in self.coords.keys():
             self.coords[key] = self.coords[key][indices]
 

lompe/model/model.py

@@ -312,8 +312,11 @@ def run_inversion(self, l1 = 0, l2 = 0,
                     spatial_weight = np.ones(ds.values.size)
 
                 dimensions = np.array(ds.values, ndmin = 2).shape[0]
-                error = np.tile(ds.error, dimensions)
-
+                if np.array(ds.error, ndmin=2).shape[0]==1:
+                    error = np.tile(ds.error, dimensions)
+                else: #error is different for different components
+                    error = ds.error.flatten()
+
                 w_i = spatial_weight * 1/(error**2) * iweights[ii]
                 if iweights[ii] != 1:
                    print('{}: Measurement uncertainty effectively changed from {} to {}'.format(dtype, np.median(error), np.median(error)/np.sqrt(iweights[ii])))