Add units back in

nifty4/library/noise_energy.py

@@ -25,23 +25,26 @@ from ..minimization.energy import Energy
 
 class NoiseEnergy(Energy):
     def __init__(self, position, d, xi, D, t, ht, Instrument,
-                 nonlinearity, alpha, q, Projection, samples=3,
-                 xi_sample_list=None, inverter=None):
+                 nonlinearity, alpha, q, Projection, munit=1., sunit=1.,
+                 dunit=1., samples=3, xi_sample_list=None, inverter=None):
         super(NoiseEnergy, self).__init__(position=position)
         self.xi = xi
         self.D = D
         self.d = d
-        self.N = DiagonalOperator(diagonal=exp(self.position))
+        self.N = DiagonalOperator(diagonal=dunit**2 * exp(self.position))
         self.t = t
         self.samples = samples
         self.ht = ht
         self.Instrument = Instrument
         self.nonlinearity = nonlinearity
+        self.munit = munit
+        self.sunit = sunit
+        self.dunit = dunit
 
         self.alpha = alpha
         self.q = q
         self.Projection = Projection
-        self.power = self.Projection.adjoint_times(exp(0.5 * self.t))
+        self.power = self.Projection.adjoint_times(munit * exp(0.5 * self.t))
         if xi_sample_list is None:
             if samples is None or samples == 0:
                 xi_sample_list = [xi]
@@ -51,14 +54,14 @@ class NoiseEnergy(Energy):
         self.xi_sample_list = xi_sample_list
         self.inverter = inverter
 
-        A = Projection.adjoint_times(exp(.5*self.t))
+        A = Projection.adjoint_times(munit * exp(.5 * self.t))  # unit: munit
 
         self._gradient = None
         for sample in self.xi_sample_list:
             map_s = self.ht(A * sample)
 
             residual = self.d - \
-                self.Instrument(self.nonlinearity(map_s))
+                self.Instrument(sunit * self.nonlinearity(map_s))
             lh = .5 * residual.vdot(self.N.inverse_times(residual))
             grad = -.5 * self.N.inverse_times(residual.conjugate()*residual)
 
@@ -81,7 +84,8 @@ class NoiseEnergy(Energy):
         return self.__class__(
             position, self.d, self.xi, self.D, self.t, self.ht,
             self.Instrument, self.nonlinearity, self.alpha, self.q,
-            self.Projection, xi_sample_list=self.xi_sample_list,
+            self.Projection, munit=self.munit, sunit=self.sunit,
+            dunit=self.dunit, xi_sample_list=self.xi_sample_list,
             samples=self.samples, inverter=self.inverter)
 
     @property

nifty4/library/nonlinear_power_curvature.py

@@ -21,11 +21,12 @@ from .response_operators import LinearizedPowerResponse
 
 
 def NonlinearPowerCurvature(tau, ht, Instrument, nonlinearity, Projection, N,
-                            T, xi_sample_list, inverter):
+                            T, xi_sample_list, inverter, munit=1., sunit=1.):
     result = None
     for xi_sample in xi_sample_list:
-        LinearizedResponse = LinearizedPowerResponse(
-            Instrument, nonlinearity, ht, Projection, tau, xi_sample)
+        LinearizedResponse = LinearizedPowerResponse(Instrument, nonlinearity,
+                                                     ht, Projection, tau,
+                                                     xi_sample, munit, sunit)
         op = LinearizedResponse.adjoint*N.inverse*LinearizedResponse
         result = op if result is None else result + op
     result = result*(1./len(xi_sample_list)) + T

nifty4/library/nonlinear_power_energy.py

@@ -53,7 +53,7 @@ class NonlinearPowerEnergy(Energy):
 
     def __init__(self, position, d, N, xi, D, ht, Instrument, nonlinearity,
                  Projection, sigma=0., samples=3, xi_sample_list=None,
-                 inverter=None):
+                 inverter=None, munit=1., sunit=1.):
         super(NonlinearPowerEnergy, self).__init__(position)
         self.xi = xi
         self.D = D
@@ -66,6 +66,8 @@ class NonlinearPowerEnergy(Energy):
         self.nonlinearity = nonlinearity
         self.Projection = Projection
         self.sigma = sigma
+        self.munit = munit
+        self.sunit = sunit
         if xi_sample_list is None:
             if samples is None or samples == 0:
                 xi_sample_list = [xi]
@@ -75,7 +77,7 @@ class NonlinearPowerEnergy(Energy):
         self.xi_sample_list = xi_sample_list
         self.inverter = inverter
 
-        A = Projection.adjoint_times(exp(.5 * position))
+        A = Projection.adjoint_times(munit * exp(.5 * position))  # unit: munit
         map_s = self.ht(A * xi)
         Tpos = self.T(position)
 
@@ -84,10 +86,10 @@ class NonlinearPowerEnergy(Energy):
             map_s = self.ht(A * xi_sample)
             LinR = LinearizedPowerResponse(
                 self.Instrument, self.nonlinearity, self.ht, self.Projection,
-                self.position, xi_sample)
+                self.position, xi_sample, munit=self.munit, sunit=self.sunit)
 
             residual = self.d - \
-                self.Instrument(self.nonlinearity(map_s))
+                self.Instrument(sunit * self.nonlinearity(map_s))
             lh = 0.5 * residual.vdot(self.N.inverse_times(residual))
             grad = LinR.adjoint_times(self.N.inverse_times(residual))
 
@@ -109,6 +111,8 @@ class NonlinearPowerEnergy(Energy):
                               self.Projection, sigma=self.sigma,
                               samples=len(self.xi_sample_list),
                               xi_sample_list=self.xi_sample_list,
+                              munit=self.munit,
+                              sunit=self.sunit,
                               inverter=self.inverter)
 
     @property
@@ -125,4 +129,4 @@ class NonlinearPowerEnergy(Energy):
         return NonlinearPowerCurvature(
             self.position, self.ht, self.Instrument, self.nonlinearity,
             self.Projection, self.N, self.T, self.xi_sample_list,
-            self.inverter)
+            self.inverter, self.munit, self.sunit)

nifty4/library/nonlinear_wiener_filter_energy.py

@@ -24,18 +24,19 @@ from .response_operators import LinearizedSignalResponse
 
 class NonlinearWienerFilterEnergy(Energy):
     def __init__(self, position, d, Instrument, nonlinearity, ht, power, N, S,
-                 inverter=None):
+                 inverter=None, sunit=1.):
         super(NonlinearWienerFilterEnergy, self).__init__(position=position)
         self.d = d
+        self.sunit = sunit
         self.Instrument = Instrument
         self.nonlinearity = nonlinearity
         self.ht = ht
         self.power = power
         m = self.ht(self.power * self.position)
         self.LinearizedResponse = LinearizedSignalResponse(
-            Instrument, nonlinearity, ht, power, m)
+            Instrument, nonlinearity, ht, power, m, sunit)
 
-        residual = d - Instrument(nonlinearity(m))
+        residual = d - Instrument(sunit * nonlinearity(m))
         self.N = N
         self.S = S
         self.inverter = inverter
@@ -48,7 +49,7 @@ class NonlinearWienerFilterEnergy(Energy):
     def at(self, position):
         return self.__class__(position, self.d, self.Instrument,
                               self.nonlinearity, self.ht, self.power, self.N,
-                              self.S, self.inverter)
+                              self.S, self.inverter, self.sunit)
 
     @property
     def value(self):

nifty4/library/response_operators.py

@@ -19,12 +19,12 @@
 from ..field import exp
 
 
-def LinearizedSignalResponse(Instrument, nonlinearity, ht, power, m):
-    return Instrument * nonlinearity.derivative(m) * ht * power
+def LinearizedSignalResponse(Instrument, nonlinearity, ht, power, m, sunit):
+    return sunit*(Instrument*nonlinearity.derivative(m)*ht*power)
 
 
-def LinearizedPowerResponse(Instrument, nonlinearity, ht, Projection, tau, xi):
-    power = exp(0.5*tau)
-    position = ht(Projection.adjoint_times(power)*xi)
+def LinearizedPowerResponse(Instrument, nonlinearity, ht, Projection, tau, xi, munit, sunit):
+    power = exp(0.5 * tau) * munit
+    position = ht(Projection.adjoint_times(power) * xi)
     linearization = nonlinearity.derivative(position)
-    return 0.5*Instrument*linearization*ht*xi*Projection.adjoint*power
+    return sunit*(0.5*Instrument*linearization*ht*xi*Projection.adjoint*power)