From 61b1f73bcc19121ed2ecb9f19ffe67564169530b Mon Sep 17 00:00:00 2001 From: "Knollmueller, Jakob (kjako)" <jakob@knollmueller.de> Date: Fri, 2 Jun 2017 14:30:38 +0200 Subject: [PATCH] critical filtering seems to work --- .gitignore | 1 + demos/critical_filtering.py | 33 +++++++++---------- demos/laplace_testing.py | 14 ++++---- .../energy_library/critical_power_energy.py | 25 +++++++------- .../critical_power_curvature.py | 8 ++--- .../operator_library/smoothness_operator.py | 4 +-- 6 files changed, 41 insertions(+), 44 deletions(-) diff --git a/.gitignore b/.gitignore index 311909456..1716a17d0 100644 --- a/.gitignore +++ b/.gitignore @@ -7,6 +7,7 @@ *.egg *.egg-info *~ +*.DS_Store *.csv .git/ .svn/ diff --git a/demos/critical_filtering.py b/demos/critical_filtering.py index 6d157fc1a..42360bda0 100644 --- a/demos/critical_filtering.py +++ b/demos/critical_filtering.py @@ -8,7 +8,7 @@ from mpi4py import MPI comm = MPI.COMM_WORLD rank = comm.rank -np.random.seed(232) +np.random.seed(42) def plot_parameters(m,t,t_true, t_real, t_d): @@ -60,8 +60,8 @@ if __name__ == "__main__": h_space = fft.target[0] # Setting up power space - p_space = PowerSpace(h_space, logarithmic=False, - distribution_strategy=distribution_strategy, nbin=128) + p_space = PowerSpace(h_space, logarithmic=True, + distribution_strategy=distribution_strategy, nbin=120) # Choosing the prior correlation structure and defining correlation operator pow_spec = (lambda k: (.05 / (k + 1) ** 2)) @@ -83,7 +83,7 @@ if __name__ == "__main__": #Adding a harmonic transformation to the instrument R = AdjointFFTResponse(fft, Instrument) - noise = 1. + noise = 5. N = DiagonalOperator(s_space, diagonal=noise, bare=True) n = Field.from_random(domain=s_space, random_type='normal', @@ -95,7 +95,7 @@ if __name__ == "__main__": realized_power = log(sh.power_analyze(logarithmic=p_space.config["logarithmic"], nbin=p_space.config["nbin"])**2) - data_power = log(fft(d).power_analyze(logarithmic=p_space.config["logarithmic"], + data_power = log(fft(d).power_analyze(logarithmic=p_space.config["logarithmic"], nbin=p_space.config["nbin"])**2) d_data = d.val.get_full_data().real if rank == 0: @@ -110,24 +110,24 @@ if __name__ == "__main__": minimizer1 = RelaxedNewton(convergence_tolerance=0, convergence_level=1, - iteration_limit=2, + iteration_limit=3, callback=convergence_measure) minimizer2 = VL_BFGS(convergence_tolerance=0, iteration_limit=50, callback=convergence_measure, - max_history_length=3) + max_history_length=10) # 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)) - #t0 = data_power + # t0 = Field(p_space, val=log(1./(1+p_space.kindex)**2)) + t0 = data_power- 1. S0 = create_power_operator(h_space, power_spectrum=exp(t0), distribution_strategy=distribution_strategy) - for i in range(100): + for i in range(500): S0 = create_power_operator(h_space, power_spectrum=exp(t0), distribution_strategy=distribution_strategy) @@ -140,14 +140,13 @@ if __name__ == "__main__": m0 = map_energy.position D0 = map_energy.curvature # Initializing the power energy with updated parameters - power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=10., samples=3) - if i > 0: - (power_energy, convergence) = minimizer1(power_energy) - else: - (power_energy, convergence) = minimizer2(power_energy) + power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=1000., samples=1) + + (power_energy, convergence) = minimizer1(power_energy) + # Setting new power spectrum - t0 = power_energy.position - t0.val[-1] = t0.val[-2] + t0.val = power_energy.position.val.real + # t0.val[-1] = t0.val[-2] # Plotting current estimate plot_parameters(m0,t0,log(sp**2),realized_power, data_power) diff --git a/demos/laplace_testing.py b/demos/laplace_testing.py index 3577720c8..db8d1e7b7 100644 --- a/demos/laplace_testing.py +++ b/demos/laplace_testing.py @@ -58,23 +58,23 @@ if __name__ == "__main__": # Setting up power space p_space = PowerSpace(h_space, logarithmic=False, - distribution_strategy=distribution_strategy, nbin=70) + distribution_strategy=distribution_strategy, nbin=200) # Choosing the prior correlation structure and defining correlation operator pow_spec = (lambda k: (.05 / (k + 1) ** 2)) # t = Field(p_space, val=pow_spec) t= Field.from_random("normal", domain=p_space) - lap = LaplaceOperator(p_space) - T = SmoothnessOperator(p_space,sigma=1.) + lap = LaplaceOperator(p_space, logarithmic=True) + T = SmoothnessOperator(p_space,sigma=1., logarithmic=True) test_energy = TestEnergy(t,T) def convergence_measure(a_energy, iteration): # returns current energy x = a_energy.value print (x, iteration) minimizer1 = VL_BFGS(convergence_tolerance=0, - iteration_limit=1000, + iteration_limit=10, callback=convergence_measure, - max_history_length=3) + max_history_length=10) def explicify(op, domain): @@ -91,7 +91,7 @@ if __name__ == "__main__": B = explicify(lap.adjoint_times, p_space) test_energy,convergence = minimizer1(test_energy) data = test_energy.position.val.get_full_data() - pl.plot([go.Scatter(x=log(p_space.kindex)[1:], y=data[1:])], filename="t.html") + pl.plot([go.Scatter(x=(p_space.kindex)[1:], y=data[1:])], filename="t.html") tt = Field.from_random("normal", domain=t.domain) print "adjointness" print t.dot(lap(tt)) @@ -100,6 +100,6 @@ if __name__ == "__main__": aa = Field(p_space, val=p_space.kindex.copy()) aa.val[0] = 1 - print lap(log(aa)).val + print lap(log(aa)**2).val print "######################" print test_energy.position.val \ No newline at end of file diff --git a/nifty/library/energy_library/critical_power_energy.py b/nifty/library/energy_library/critical_power_energy.py index 81c75c005..46edc70b6 100644 --- a/nifty/library/energy_library/critical_power_energy.py +++ b/nifty/library/energy_library/critical_power_energy.py @@ -1,6 +1,6 @@ from nifty.energies.energy import Energy from nifty.library.operator_library import CriticalPowerCurvature,\ - LaplaceOperator + SmoothnessOperator from nifty.sugar import generate_posterior_sample from nifty import Field, exp @@ -31,14 +31,12 @@ class CriticalPowerEnergy(Energy): self.sigma = sigma self.alpha = Field(self.position.domain, val=alpha) self.q = Field(self.position.domain, val=q) - #self.T = SmoothnessOperator(domain=self.position.domain[0], sigma=self.sigma) - self.Laplace = LaplaceOperator(self.position.domain[0]) - self.roughness = self.Laplace(self.position) + self.T = SmoothnessOperator(domain=self.position.domain[0], sigma=self.sigma) self.rho = self.position.domain[0].rho self.w = w if self.w is None: self.w = self._calculate_w(self.m, self.D, self.samples) - self.theta = exp(-self.position) * (self.q + self.w / 2.) + self.theta = (exp(-self.position) * (self.q + self.w / 2.)) def at(self, position): return self.__class__(position, self.m, D=self.D, @@ -49,22 +47,22 @@ class CriticalPowerEnergy(Energy): @property def value(self): - energy = exp(-self.position).dot(self.q + self.w / 2.) - energy += self.position.dot(self.alpha - 1. + self.rho / 2.) - energy += 0.5 * self.roughness.dot(self.roughness) / self.sigma**2 + energy = exp(-self.position).dot(self.q + self.w / 2., bare= True) + energy += self.position.dot(self.alpha - 1. + self.rho / 2., bare=True) + energy += 0.5 * self.position.dot(self.T(self.position)) return energy.real @property def gradient(self): - gradient = - self.theta - gradient += self.alpha - 1. + self.rho / 2. - gradient += self.Laplace(self.roughness) / self.sigma **2 + gradient = - self.theta.weight(-1) + gradient += (self.alpha - 1. + self.rho / 2.).weight(-1) + gradient += self.T(self.position) + gradient.val = gradient.val.real return gradient @property def curvature(self): - curvature = CriticalPowerCurvature(theta=self.theta, Laplace = self.Laplace, - sigma = self.sigma) + curvature = CriticalPowerCurvature(theta=self.theta.weight(-1), T=self.T) return curvature def _calculate_w(self, m, D, samples): @@ -81,6 +79,7 @@ class CriticalPowerEnergy(Energy): else: w = m.power_analyze( logarithmic=self.position.domain[0].config["logarithmic"], + nbin=self.position.domain[0].config["nbin"], decompose_power=False) w **= 2 w *= self.rho diff --git a/nifty/library/operator_library/critical_power_curvature.py b/nifty/library/operator_library/critical_power_curvature.py index ece3cb58b..3f5478071 100644 --- a/nifty/library/operator_library/critical_power_curvature.py +++ b/nifty/library/operator_library/critical_power_curvature.py @@ -3,11 +3,10 @@ from nifty.operators import EndomorphicOperator,\ class CriticalPowerCurvature(InvertibleOperatorMixin, EndomorphicOperator): - def __init__(self, theta, Laplace, sigma, inverter=None, preconditioner=None): + def __init__(self, theta, T, inverter=None, preconditioner=None): self.theta = theta - self.Laplace = Laplace - self.sigma = sigma + self.T = T # if preconditioner is None: # preconditioner = self.T.times self._domain = self.theta.domain @@ -28,5 +27,4 @@ class CriticalPowerCurvature(InvertibleOperatorMixin, EndomorphicOperator): # ---Added properties and methods--- def _times(self, x, spaces): - return self.Laplace.adjoint_times(self.Laplace(x)) / self.sigma ** 2 \ - + self.theta * x + return self.T(x) + self.theta * x diff --git a/nifty/library/operator_library/smoothness_operator.py b/nifty/library/operator_library/smoothness_operator.py index fca9d91e2..50e62f29a 100644 --- a/nifty/library/operator_library/smoothness_operator.py +++ b/nifty/library/operator_library/smoothness_operator.py @@ -6,7 +6,7 @@ from laplace_operator import LaplaceOperator class SmoothnessOperator(EndomorphicOperator): - def __init__(self, domain, sigma=1., + def __init__(self, domain, sigma=1.,logarithmic = True, default_spaces=None): super(SmoothnessOperator, self).__init__(default_spaces) @@ -21,7 +21,7 @@ class SmoothnessOperator(EndomorphicOperator): self._sigma = sigma - self._Laplace = LaplaceOperator(domain=self._domain[0]) + self._Laplace = LaplaceOperator(domain=self._domain[0], logarithmic=logarithmic) """ SmoothnessOperator -- GitLab