merge from master

.gitignore

@@ -2,6 +2,11 @@
 setup.cfg
 .idea
 .DS_Store
+*.pyc
+*.html
+.document
+.svn/
+*.csv
 
 # from https://github.com/github/gitignore/blob/master/Python.gitignore
 
@@ -95,4 +100,4 @@ ENV/
 .spyderproject
 
 # Rope project settings
-.ropeproject
\ No newline at end of file
+.ropeproject

demos/critical_filtering.py

+
+from nifty import *
+from nifty.library.wiener_filter import WienerFilterEnergy
+from nifty.library.critical_filter import CriticalPowerEnergy
+import plotly.offline as pl
+import plotly.graph_objs as go
+
+from mpi4py import MPI
+comm = MPI.COMM_WORLD
+rank = comm.rank
+
+np.random.seed(42)
+
+
+def plot_parameters(m,t,p, p_d):
+
+    x = log(t.domain[0].kindex)
+    m = fft.adjoint_times(m)
+    m = m.val.get_full_data().real
+    t = t.val.get_full_data().real
+    p = p.val.get_full_data().real
+    p_d = p_d.val.get_full_data().real
+    pl.plot([go.Heatmap(z=m)], filename='map.html')
+    pl.plot([go.Scatter(x=x,y=t), go.Scatter(x=x ,y=p), go.Scatter(x=x, y=p_d)], filename="t.html")
+
+
+class AdjointFFTResponse(LinearOperator):
+    def __init__(self, FFT, R, default_spaces=None):
+        super(AdjointFFTResponse, self).__init__(default_spaces)
+        self._domain = FFT.target
+        self._target = R.target
+        self.R = R
+        self.FFT = FFT
+
+    def _times(self, x, spaces=None):
+        return self.R(self.FFT.adjoint_times(x))
+
+    def _adjoint_times(self, x, spaces=None):
+        return self.FFT(self.R.adjoint_times(x))
+    @property
+    def domain(self):
+        return self._domain
+
+    @property
+    def target(self):
+        return self._target
+
+    @property
+    def unitary(self):
+        return False
+
+if __name__ == "__main__":
+
+    distribution_strategy = 'not'
+
+    # Set up position space
+    s_space = RGSpace([128,128])
+    # s_space = HPSpace(32)
+
+    # Define harmonic transformation and associated harmonic space
+    fft = FFTOperator(s_space)
+    h_space = fft.target[0]
+
+    # Setting up power space
+    p_space = PowerSpace(h_space, logarithmic=True,
+                         distribution_strategy=distribution_strategy)
+
+    # Choosing the prior correlation structure and defining correlation operator
+    p_spec = (lambda k: (.5 / (k + 1) ** 3))
+    S = create_power_operator(h_space, power_spectrum=p_spec,
+                              distribution_strategy=distribution_strategy)
+
+    # Drawing a sample sh from the prior distribution in harmonic space
+    sp = Field(p_space,  val=p_spec,
+               distribution_strategy=distribution_strategy)
+    sh = sp.power_synthesize(real_signal=True)
+
+
+    # Choosing the measurement instrument
+    # Instrument = SmoothingOperator(s_space, sigma=0.01)
+    Instrument = DiagonalOperator(s_space, diagonal=1.)
+    # Instrument._diagonal.val[200:400, 200:400] = 0
+    #Instrument._diagonal.val[64:512-64, 64:512-64] = 0
+
+
+    #Adding a harmonic transformation to the instrument
+    R = AdjointFFTResponse(fft, Instrument)
+
+    noise = 1.
+    N = DiagonalOperator(s_space, diagonal=noise, bare=True)
+    n = Field.from_random(domain=s_space,
+                          random_type='normal',
+                          std=sqrt(noise),
+                          mean=0)
+
+    # Creating the mock data
+    d = R(sh) + n
+
+    # The information source
+    j = R.adjoint_times(N.inverse_times(d))
+    realized_power = log(sh.power_analyze(logarithmic=p_space.config["logarithmic"],
+                                          nbin=p_space.config["nbin"]))
+    data_power = log(fft(d).power_analyze(logarithmic=p_space.config["logarithmic"],
+                                          nbin=p_space.config["nbin"]))
+    d_data = d.val.get_full_data().real
+    if rank == 0:
+        pl.plot([go.Heatmap(z=d_data)], filename='data.html')
+
+    #  minimization strategy
+
+    def convergence_measure(a_energy, iteration): # returns current energy
+        x = a_energy.value
+        print (x, iteration)
+
+
+    minimizer1 = RelaxedNewton(convergence_tolerance=10e-2,
+                              convergence_level=2,
+                              iteration_limit=3,
+                              callback=convergence_measure)
+
+    minimizer2 = VL_BFGS(convergence_tolerance=0,
+                       iteration_limit=7,
+                       callback=convergence_measure,
+                       max_history_length=3)
+
+    # Setting starting position
+    flat_power = Field(p_space,val=10e-8)
+    m0 = flat_power.power_synthesize(real_signal=True)
+
+    t0 = Field(p_space, val=log(1./(1+p_space.kindex)**2))
+
+
+
+    for i in range(500):
+        S0 = create_power_operator(h_space, power_spectrum=exp(t0),
+                              distribution_strategy=distribution_strategy)
+
+        # Initializing the  nonlinear Wiener Filter energy
+        map_energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0, inverter=inverter)
+        # Solving the Wiener Filter analytically
+        D0 = map_energy.curvature
+        m0 = D0.inverse_times(j)
+        # Initializing the power energy with updated parameters
+        power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, smoothness_prior=10., samples=3)
+
+        (power_energy, convergence) = minimizer1(power_energy)
+
+
+        # Setting new power spectrum
+        t0.val  = power_energy.position.val.real
+
+        # Plotting current estimate
+        print i
+        if i%50 == 0:
+            plot_parameters(m0,t0,log(sp), data_power)
+
+

demos/wiener_filter_advanced.py

+
+from nifty import *
+
+import plotly.offline as pl
+import plotly.graph_objs as go
+
+from mpi4py import MPI
+comm = MPI.COMM_WORLD
+rank = comm.rank
+
+np.random.seed(42)
+
+class AdjointFFTResponse(LinearOperator):
+    def __init__(self, FFT, R, default_spaces=None):
+        super(AdjointFFTResponse, self).__init__(default_spaces)
+        self._domain = FFT.target
+        self._target = R.target
+        self.R = R
+        self.FFT = FFT
+
+    def _times(self, x, spaces=None):
+        return self.R(self.FFT.adjoint_times(x))
+
+    def _adjoint_times(self, x, spaces=None):
+        return self.FFT(self.R.adjoint_times(x))
+    @property
+    def domain(self):
+        return self._domain
+
+    @property
+    def target(self):
+        return self._target
+
+    @property
+    def unitary(self):
+        return False
+
+
+
+if __name__ == "__main__":
+
+    distribution_strategy = 'not'
+
+    # Set up position space
+    s_space = RGSpace([128,128])
+    # s_space = HPSpace(32)
+
+    # Define harmonic transformation and associated harmonic space
+    fft = FFTOperator(s_space)
+    h_space = fft.target[0]
+
+    # Setting up power space
+    p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
+
+    # Choosing the prior correlation structure and defining correlation operator
+    p_spec = (lambda k: (42 / (k + 1) ** 3))
+
+    S = create_power_operator(h_space, power_spectrum=p_spec,
+                              distribution_strategy=distribution_strategy)
+
+    # Drawing a sample sh from the prior distribution in harmonic space
+    sp = Field(p_space, val=p_spec,
+               distribution_strategy=distribution_strategy)
+    sh = sp.power_synthesize(real_signal=True)
+    ss = fft.adjoint_times(sh)
+
+    # Choosing the measurement instrument
+    # Instrument = SmoothingOperator(s_space, sigma=0.05)
+    Instrument = DiagonalOperator(s_space, diagonal=1.)
+#    Instrument._diagonal.val[200:400, 200:400] = 0
+
+    #Adding a harmonic transformation to the instrument
+    R = AdjointFFTResponse(fft, Instrument)
+    signal_to_noise = 1.
+    N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
+    n = Field.from_random(domain=s_space,
+                          random_type='normal',
+                          std=ss.std()/np.sqrt(signal_to_noise),
+                          mean=0)
+
+    # Creating the mock data
+    d = R(sh) + n
+    j = R.adjoint_times(N.inverse_times(d))
+
+    # Choosing the minimization strategy
+
+    def convergence_measure(energy, iteration): # returns current energy
+        x = energy.value
+        print (x, iteration)
+
+#    minimizer = SteepestDescent(convergence_tolerance=0,
+#                                iteration_limit=50,
+#                                callback=convergence_measure)
+
+    minimizer = RelaxedNewton(convergence_tolerance=0,
+                              iteration_limit=1,
+                              callback=convergence_measure)
+    #
+    # minimizer = VL_BFGS(convergence_tolerance=0,
+    #                    iteration_limit=50,
+    #                    callback=convergence_measure,
+    #                    max_history_length=3)
+    #
+
+    inverter = ConjugateGradient(convergence_level=3,
+                                 convergence_tolerance=10e-5,
+                                 preconditioner=None)
+    # Setting starting position
+    m0 = Field(h_space, val=.0)
+
+    # Initializing the Wiener Filter energy
+    energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S, inverter=inverter)
+    D0 = energy.curvature
+
+    # Solving the problem analytically
+    m0 = D0.inverse_times(j)
+
+    sample_variance = Field(sh.domain,val=0. + 0j)
+    sample_mean = Field(sh.domain,val=0. + 0j)
+
+    # sampling the uncertainty map
+    n_samples = 1
+    for i in range(n_samples):
+        sample = sugar.generate_posterior_sample(m0,D0)
+        sample_variance += sample**2
+        sample_mean += sample
+    variance = sample_variance/n_samples - (sample_mean/n_samples)
+

demos/wiener_filter.py → demos/wiener_filter_easy.py

@@ -20,7 +20,7 @@ if __name__ == "__main__":
     correlation_length = 0.1
     #variance of field in position space sqrt(<|s_x|^2>) (in unit of s)
     field_variance = 2.
-    #smoothing length that response (in same unit as L)
+    #smoothing length of response (in same unit as L)
     response_sigma = 0.1
 
     #defining resolution (pixels per dimension)
@@ -63,15 +63,8 @@ if __name__ == "__main__":
     d = R(ss) + n
 
     # Wiener filter
+
     j = R.adjoint_times(N.inverse_times(d))
     D = PropagatorOperator(S=S, N=N, R=R)
 
     m = D(j)
-
-    d_data = d.val.get_full_data().real
-    m_data = m.val.get_full_data().real
-    ss_data = ss.val.get_full_data().real
-#    if rank == 0:
-#        pl.plot([go.Heatmap(z=d_data)], filename='data.html')
-#        pl.plot([go.Heatmap(z=m_data)], filename='map.html')
-#        pl.plot([go.Heatmap(z=ss_data)], filename='map_orig.html')

demos/wiener_filter_hamiltonian.py

-from __future__ import division
-from __future__ import print_function
-from builtins import object
-
-from nifty import *
-
-import plotly.offline as pl
-import plotly.graph_objs as go
-
-from mpi4py import MPI
-comm = MPI.COMM_WORLD
-rank = comm.rank
-
-np.random.seed(42)
-
-class WienerFilterEnergy(Energy):
-    def __init__(self, position, D, j):
-        # in principle not necessary, but useful in order to make the signature
-        # explicit
-        super(WienerFilterEnergy, self).__init__(position)
-        self.D = D
-        self.j = j
-
-    def at(self, position):
-        return self.__class__(position, D=self.D, j=self.j)
-
-    @property
-    def value(self):
-        D_inv_x = self.D_inverse_x()
-        H = 0.5 * D_inv_x.vdot(self.position) - self.j.vdot(self.position)
-        return H.real
-
-    @property
-    def gradient(self):
-        D_inv_x = self.D_inverse_x()
-        g = D_inv_x - self.j
-        return_g = g.copy_empty(dtype=np.float)
-        return_g.val = g.val.real
-        return return_g
-
-    @property
-    def curvature(self):
-        class Dummy(object):
-            def __init__(self, x):
-                self.x = x
-            def inverse_times(self, *args, **kwargs):
-                return self.x.times(*args, **kwargs)
-        my_dummy = Dummy(self.D)
-        return my_dummy
-
-
-    @memo
-    def D_inverse_x(self):
-        return D.inverse_times(self.position)
-
-
-if __name__ == "__main__":
-
-    distribution_strategy = 'not'
-
-    # Set up spaces and fft transformation
-    s_space = RGSpace([512, 512])
-    fft = FFTOperator(s_space)
-    h_space = fft.target[0]
-    p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
-
-    # create the field instances and power operator
-    pow_spec = (lambda k: (42. / (k + 1) ** 3))
-    S = create_power_operator(h_space, power_spectrum=pow_spec,
-                              distribution_strategy=distribution_strategy)
-
-    sp = Field(p_space, val=lambda z: pow_spec(z)**0.5,
-               distribution_strategy=distribution_strategy)
-    sh = sp.power_synthesize(real_signal=True)
-    ss = fft.inverse_times(sh)
-
-    # model the measurement process
-    R = SmoothingOperator.make(s_space, sigma=0.01)
-#    R = DiagonalOperator(s_space, diagonal=1.)
-#    R._diagonal.val[200:400, 200:400] = 0
-
-    signal_to_noise = 1
-    N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
-    n = Field.from_random(domain=s_space,
-                          random_type='normal',
-                          std=ss.std()/np.sqrt(signal_to_noise),
-                          mean=0)
-
-    # create mock data
-    d = R(ss) + n
-
-    # set up reconstruction objects
-    j = R.adjoint_times(N.inverse_times(d))
-    D = PropagatorOperator(S=S, N=N, R=R)
-
-    def distance_measure(energy, iteration):
-        x = energy.position
-        print(iteration, (x-ss).norm()/ss.norm().real)
-
-#    minimizer = SteepestDescent(convergence_tolerance=0,
-#                                iteration_limit=50,
-#                                callback=distance_measure)
-
-    minimizer = RelaxedNewton(convergence_tolerance=0,
-                              iteration_limit=2,
-                              callback=distance_measure)
-
-#    minimizer = VL_BFGS(convergence_tolerance=0,
-#                        iteration_limit=50,
-#                        callback=distance_measure,
-#                        max_history_length=3)
-
-    m0 = Field(s_space, val=1.)
-
-    energy = WienerFilterEnergy(position=m0, D=D, j=j)
-
-    (energy, convergence) = minimizer(energy)
-
-    m = energy.position
-
-    d_data = d.val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Heatmap(z=d_data)], filename='data.html')
-
-
-    ss_data = ss.val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Heatmap(z=ss_data)], filename='ss.html')
-
-    sh_data = sh.val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Heatmap(z=sh_data)], filename='sh.html')
-
-    j_data = j.val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Heatmap(z=j_data)], filename='j.html')
-
-    jabs_data = np.abs(j.val.get_full_data())
-    jphase_data = np.angle(j.val.get_full_data())
-    if rank == 0:
-        pl.plot([go.Heatmap(z=jabs_data)], filename='j_abs.html')
-        pl.plot([go.Heatmap(z=jphase_data)], filename='j_phase.html')
-
-    m_data = m.val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Heatmap(z=m_data)], filename='map.html')
-
-#    grad_data = grad.val.get_full_data().real
-#    if rank == 0:
-#        pl.plot([go.Heatmap(z=grad_data)], filename='grad.html')

demos/wiener_filter_harmonic.py

-from __future__ import division
-from builtins import range
-from nifty import *
-from mpi4py import MPI
-import plotly.offline as py
-import plotly.graph_objs as go
-
-
-comm = MPI.COMM_WORLD
-rank = comm.rank
-
-
-def plot_maps(x, name):
-
-    trace = [None]*len(x)
-
-    keys = list(x.keys())
-    field = x[keys[0]]
-    domain = field.domain[0]
-    shape = len(domain.shape)
-    max_n = domain.shape[0]*domain.distances[0]
-    step = domain.distances[0]
-    x_axis = np.arange(0, max_n, step)
-
-    if shape == 1:
-        for ii in range(len(x)):
-            trace[ii] = go.Scatter(x= x_axis, y=x[keys[ii]].val.get_full_data(), name=keys[ii])
-        fig = go.Figure(data=trace)
-
-        py.plot(fig, filename=name)
-
-    elif shape == 2:
-        for ii in range(len(x)):
-            py.plot([go.Heatmap(z=x[keys[ii]].val.get_full_data())], filename=keys[ii])
-    else:
-        raise TypeError("Only 1D and 2D field plots are supported")
-
-def plot_power(x, name):
-
-    layout = go.Layout(
-        xaxis=dict(
-            type='log',
-            autorange=True
-        ),
-        yaxis=dict(
-            type='log',
-            autorange=True
-        )
-    )
-
-    trace = [None]*len(x)
-
-    keys = list(x.keys())
-    field = x[keys[0]]
-    domain = field.domain[0]
-    x_axis = domain.kindex
-
-    for ii in range(len(x)):
-        trace[ii] = go.Scatter(x= x_axis, y=x[keys[ii]].val.get_full_data(), name=keys[ii])
-
-    fig = go.Figure(data=trace, layout=layout)
-    py.plot(fig, filename=name)
-
-np.random.seed(42)
-
-
-
-if __name__ == "__main__":
-
-    distribution_strategy = 'not'
-
-    # setting spaces
-    npix = np.array([500])  # number of pixels
-    total_volume = 1.  # total length
-
-    # setting signal parameters
-    lambda_s = .05  # signal correlation length
-    sigma_s = 10.  # signal variance
-
-
-    #setting response operator parameters
-    length_convolution = .025
-    exposure = 1.
-
-    # calculating parameters
-    k_0 = 4. / (2 * np.pi * lambda_s)
-    a_s = sigma_s ** 2. * lambda_s * total_volume
-
-    # creation of spaces
-    # x1 = RGSpace([npix,npix], distances=total_volume / npix,
-    #              zerocenter=False)
-    # k1 = RGRGTransformation.get_codomain(x1)
-
-    x1 = HPSpace(64)
-    k1 = HPLMTransformation.get_codomain(x1)
-    p1 = PowerSpace(harmonic_partner=k1, logarithmic=False)
-
-    # creating Power Operator with given spectrum
-    spec = (lambda k: a_s / (1 + (k / k_0) ** 2) ** 2)
-    p_field = Field(p1, val=spec)
-    S_op = create_power_operator(k1, spec)
-
-    # creating FFT-Operator and Response-Operator with Gaussian convolution
-    # adjust dtype_target probperly
-    Fft_op = FFTOperator(domain=x1, target=k1,
-                        domain_dtype=np.float64,
-                        target_dtype=np.float64)
-    R_op = ResponseOperator(x1, sigma=[length_convolution],
-                            exposure=[exposure])
-
-    # drawing a random field
-    sk = p_field.power_synthesize(real_power=True, mean=0.)
-    s = Fft_op.adjoint_times(sk)
-
-    signal_to_noise = 1
-    N_op = DiagonalOperator(R_op.target, diagonal=s.var()/signal_to_noise, bare=True)
-    n = Field.from_random(domain=R_op.target,
-                          random_type='normal',
-                          std=s.std()/np.sqrt(signal_to_noise),
-                          mean=0)
-
-    d = R_op(s) + n
-
-    # Wiener filter
-    j = Fft_op.times(R_op.adjoint_times(N_op.inverse_times(d)))
-    D = HarmonicPropagatorOperator(S=S_op, N=N_op, R=R_op)
-
-    mk = D(j)
-
-    m = Fft_op.adjoint_times(mk)
-
-    # z={}
-    # z["signal"] = s
-    # z["reconstructed_map"] = m
-    # z["data"] = d
-    # z["lambda"] = R_op(s)
-    #