Commit 10966a24 authored by Jakob Knollmueller's avatar Jakob Knollmueller

working on laplace

parent 3c0b82ce
Pipeline #13201 failed with stage
in 6 minutes and 17 seconds
File added
......@@ -8,33 +8,33 @@ from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.rank
np.random.seed(62)
np.random.seed(232)
class NonlinearResponse(LinearOperator):
def __init__(self, FFT, Instrument, function, derivative, default_spaces=None):
super(NonlinearResponse, self).__init__(default_spaces)
def plot_parameters(m,t,t_true, t_real, t_d):
m = fft.adjoint_times(m)
m_data = m.val.get_full_data().real
t_data = t.val.get_full_data().real
t_d_data = t_d.val.get_full_data().real
t_true_data = t_true.val.get_full_data().real
t_real_data = t_real.val.get_full_data().real
pl.plot([go.Heatmap(z=m_data)], filename='map.html')
pl.plot([go.Scatter(y=t_data), go.Scatter(y=t_true_data),
go.Scatter(y=t_real_data), go.Scatter(y=t_d_data)], 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 = Instrument.target
self.function = function
self.derivative = derivative
self.I = Instrument
self._target = R.target
self.R = R
self.FFT = FFT
def _times(self, x, spaces=None):
return self.I(self.function(self.FFT.adjoint_times(x)))
return self.R(self.FFT.adjoint_times(x))
def _adjoint_times(self, x, spaces=None):
return self.FFT(self.function(self.I.adjoint_times(x)))
def derived_times(self, x, position):
position_derivative = self.derivative(self.FFT.adjoint_times(position))
return self.I(position_derivative * self.FFT.adjoint_times(x))
def derived_adjoint_times(self, x, position):
position_derivative = self.derivative(self.FFT.adjoint_times(position))
return self.FFT(position_derivative * self.I.adjoint_times(x))
return self.FFT(self.R.adjoint_times(x))
@property
def domain(self):
return self._domain
......@@ -47,16 +47,6 @@ class NonlinearResponse(LinearOperator):
def unitary(self):
return False
def plot_parameters(m,t,t_true, t_real):
m = fft.adjoint_times(m)
m_data = m.val.get_full_data().real
t_data = t.val.get_full_data().real
t_true_data = t_true.val.get_full_data().real
t_real_data = t_real.val.get_full_data().real
pl.plot([go.Heatmap(z=m_data)], filename='map.html')
pl.plot([go.Scatter(y=t_data), go.Scatter(y=t_true_data),
go.Scatter(y=t_real_data)], filename="t.html")
if __name__ == "__main__":
distribution_strategy = 'not'
......@@ -70,11 +60,11 @@ if __name__ == "__main__":
h_space = fft.target[0]
# Setting up power space
p_space = PowerSpace(h_space, logarithmic = True,
distribution_strategy=distribution_strategy)
p_space = PowerSpace(h_space, logarithmic=False,
distribution_strategy=distribution_strategy, nbin=128)
# Choosing the prior correlation structure and defining correlation operator
pow_spec = (lambda k: (.05 / (k + 1) ** 3))
pow_spec = (lambda k: (.05 / (k + 1) ** 2))
S = create_power_operator(h_space, power_spectrum=pow_spec,
distribution_strategy=distribution_strategy)
......@@ -89,43 +79,11 @@ if __name__ == "__main__":
Instrument = DiagonalOperator(s_space, diagonal=1.)
# Instrument._diagonal.val[200:400, 200:400] = 0
# Choosing nonlinearity
# log-normal model:
function = exp
derivative = exp
# tan-normal model
# def function(x):
# return 0.5 * tanh(x) + 0.5
# def derivative(x):
# return 0.5*(1 - tanh(x)**2)
# no nonlinearity, Wiener Filter
# def function(x):
# return x
# def derivative(x):
# return 1
# small quadratic pertubarion
# def function(x):
# return 0.5*x**2 + x
# def derivative(x):
# return x + 1
# def function(x):
# return 0.9*x**4 +0.2*x**2 + x
# def derivative(x):
# return 0.9*4*x**3 + 0.4*x +1
#
#Adding a harmonic transformation to the instrument
R = NonlinearResponse(fft, Instrument, function, derivative)
noise = .01
R = AdjointFFTResponse(fft, Instrument)
noise = 1.
N = DiagonalOperator(s_space, diagonal=noise, bare=True)
n = Field.from_random(domain=s_space,
random_type='normal',
......@@ -134,66 +92,64 @@ if __name__ == "__main__":
# Creating the mock data
d = R(sh) + n
realized_power = log(sh.power_analyze(logarithmic=p_space.config["logarithmic"])**2)
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"],
nbin=p_space.config["nbin"])**2)
d_data = d.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=d_data)], filename='data.html')
# Choosing the minimization strategy
# minimization strategy
def convergence_measure(a_energy, iteration): # returns current energy
x = a_energy.value
print (x, iteration)
# minimizer1 = SteepestDescent(convergence_tolerance=0,
# iteration_limit=50,
# callback=convergence_measure)
minimizer1 = RelaxedNewton(convergence_tolerance=0,
convergence_level=1,
iteration_limit=5,
iteration_limit=2,
callback=convergence_measure)
# minimizer2 = RelaxedNewton(convergence_tolerance=0,
# convergence_level=1,
# iteration_limit=2,
# callback=convergence_measure)
#
# minimizer1 = VL_BFGS(convergence_tolerance=0,
# iteration_limit=5,
# callback=convergence_measure,
# max_history_length=3)
minimizer2 = VL_BFGS(convergence_tolerance=0,
iteration_limit=50,
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))
# t0 = Field(p_space,val=-8)
# t0 = log(sp.copy()**2)
#t0 = data_power
S0 = create_power_operator(h_space, power_spectrum=exp(t0),
distribution_strategy=distribution_strategy)
for i in range(100):
S0 = create_power_operator(h_space, power_spectrum=exp(t0),
distribution_strategy=distribution_strategy)
# Initializing the nonlinear Wiener Filter energy
map_energy = NonlinearWienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0)
map_energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0)
# Minimization with chosen minimizer
(map_energy, convergence) = minimizer1(map_energy)
map_energy = map_energy.analytic_solution()
# Updating parameters for correlation structure reconstruction
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)
(power_energy, convergence) = minimizer1(power_energy)
if i > 0:
(power_energy, convergence) = minimizer1(power_energy)
else:
(power_energy, convergence) = minimizer2(power_energy)
# Setting new power spectrum
t0 = power_energy.position
plot_parameters(m0,t0,log(sp**2),realized_power)
t0.val[-1] = t0.val[-2]
# Plotting current estimate
plot_parameters(m0,t0,log(sp**2),realized_power, data_power)
# Transforming fields to position space for plotting
......@@ -223,10 +179,3 @@ if __name__ == "__main__":
m_data = m.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=m_data)], filename='map.html')
f_m_data = function(m).val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=f_m_data)], filename='f_map.html')
f_ss_data = function(ss).val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=f_ss_data)], filename='f_ss.html')
from nifty import *
import numpy as np
from nifty import Field,\
EndomorphicOperator,\
PowerSpace
import plotly.offline as pl
import plotly.graph_objs as go
import numpy as np
from nifty import Field, \
EndomorphicOperator, \
PowerSpace
class TestEnergy(Energy):
def __init__(self, position, Op):
super(TestEnergy, self).__init__(position)
self.Op = Op
def at(self, position):
return self.__class__(position=position, Op=self.Op)
@property
def value(self):
return 0.5 * self.position.dot(self.Op(self.position))
@property
def gradient(self):
return self.Op(self.position)
@property
def curvature(self):
curv = CurvOp(self.Op)
return curv
class CurvOp(InvertibleOperatorMixin, EndomorphicOperator):
def __init__(self, Op,inverter=None, preconditioner=None):
self.Op = Op
self._domain = Op.domain
super(CurvOp, self).__init__(inverter=inverter,
preconditioner=preconditioner)
def _times(self, x, spaces):
return self.Op(x)
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=False,
distribution_strategy=distribution_strategy, nbin=16)
# 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 = LogLaplaceOperator(p_space)
T = SmoothnessOperator(p_space,sigma=1.)
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,
callback=convergence_measure,
max_history_length=3)
def explicify(op, domain):
space = domain
d = space.dim
res = np.zeros((d, d))
for i in range(d):
x = np.zeros(d)
x[i] = 1.
f = Field(space, val=x)
res[:, i] = op(f).val
return res
A = explicify(lap.times, p_space)
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")
tt = Field.from_random("normal", domain=t.domain)
print "adjointness"
print t.dot(lap(tt))
print tt.dot(lap.adjoint_times(t))
print "log kindex"
aa = Field(p_space, val=p_space.kindex.copy())
aa.val[0] = 1
print lap(log(aa)).val
print "######################"
print test_energy.position.val
\ No newline at end of file
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 NonlinearResponse(LinearOperator):
def __init__(self, FFT, Instrument, function, derivative, default_spaces=None):
super(NonlinearResponse, self).__init__(default_spaces)
self._domain = FFT.target
self._target = Instrument.target
self.function = function
self.derivative = derivative
self.I = Instrument
self.FFT = FFT
def _times(self, x, spaces=None):
return self.I(self.function(self.FFT.adjoint_times(x)))
def _adjoint_times(self, x, spaces=None):
return self.FFT(self.function(self.I.adjoint_times(x)))
def derived_times(self, x, position):
position_derivative = self.derivative(self.FFT.adjoint_times(position))
return self.I(position_derivative * self.FFT.adjoint_times(x))
def derived_adjoint_times(self, x, position):
position_derivative = self.derivative(self.FFT.adjoint_times(position))
return self.FFT(position_derivative * self.I.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([100,100])
# 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
pow_spec = (lambda k: (1. / (k + 1) ** 3))
S = create_power_operator(h_space, power_spectrum=pow_spec,
distribution_strategy=distribution_strategy)
# Drawing a sample sh from the prior distribution in harmonic space
sp = Field(p_space, val=lambda z: pow_spec(z)**(1./2),
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
# Choosing nonlinearity
#
# function = exp
# derivative = exp
def function(x):
return 0.5 * tanh(x) + 0.5
def derivative(x):
return 0.5*(1. - tanh(x)**2)
# def function(x):
# return 20*x
# def derivative(x):
# return 20.
# def function(x):
# return 0.1*0.5*x**2 + x
# def derivative(x):
# return 0.1*x + 1
# def function(x):
# return 0.9*x**4 +0.2*x**2 + x
# def derivative(x):
# return 0.9*4*x**3 + 0.4*x +1
#
#Adding a harmonic transformation to the instrument
noise = 10.
N = DiagonalOperator(s_space, diagonal=noise, bare=True)
n = Field.from_random(domain=s_space,
random_type='normal',
std=sqrt(noise),
mean=0)
R = NonlinearResponse(fft, Instrument, function, derivative)
# Creating the mock data
d = R(sh) + n
# 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,
convergence_level=1,
iteration_limit=3,
callback=convergence_measure)
# minimizer = VL_BFGS(convergence_tolerance=0,
# iteration_limit=50,
# 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)
# Initializing the Wiener Filter energy
energy = NonlinearWienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S)
# Solving the problem with chosen minimization strategy
(energy, convergence) = minimizer(energy)
# Transforming fields to position space for plotting
ss = fft.adjoint_times(sh)
m = fft.adjoint_times(energy.position)
# Plotting
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')
m_data = m.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=m_data)], filename='map.html')
f_m_data = function(m).val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=f_m_data)], filename='f_map.html')
f_ss_data = function(ss).val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=f_ss_data)], filename='f_ss.html')
......@@ -54,12 +54,36 @@ def plot_parameters(m,t, t_real):
pl.plot([go.Scatter(y=m_data)], filename='map.html')
pl.plot([go.Scatter(y=t_data),
go.Scatter(y=t_real_data)], 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'
full_data = np.genfromtxt("train_data.csv", delimiter = ',')
d = full_data.T[2]
d = full_data.T[3]
d[0] = 0.
d -= d.mean()
d[0] = 0.
......@@ -80,42 +104,9 @@ if __name__ == "__main__":
Instrument = DiagonalOperator(s_space, diagonal=1.)
# Instrument._diagonal.val[200:400, 200:400] = 0
# Choosing nonlinearity
# log-normal model:
# function = exp
# derivative = exp
# tan-normal model
# def function(x):
# return 0.5 * tanh(x) + 0.5
# def derivative(x):
# return 0.5*(1 - tanh(x)**2)
# no nonlinearity, Wiener Filter
def function(x):
return x
def derivative(x):
return 1
# small quadratic pertubarion
# def function(x):
# return 0.5*x**2 + x
# def derivative(x):
# return x + 1
# def function(x):
# return 0.9*x**4 +0.2*x**2 + x
# def derivative(x):
# return 0.9*4*x**3 + 0.4*x +1
#
#Adding a harmonic transformation to the instrument
R = NonlinearResponse(fft, Instrument, function, derivative)
R = AdjointFFTResponse(fft, Instrument)
noise = .1
N = DiagonalOperator(s_space, diagonal=noise, bare=True)
......@@ -139,13 +130,13 @@ if __name__ == "__main__":
callback=convergence_measure)
minimizer2 = RelaxedNewton(convergence_tolerance=0,
convergence_level=1,
iteration_limit=10,
iteration_limit=30,
callback=convergence_measure)
# minimizer1 = VL_BFGS(convergence_tolerance=0,
# iteration_limit=5,
# callback=convergence_measure,
# max_history_length=3)
minimizer1 = VL_BFGS(convergence_tolerance=0,
iteration_limit=30,
callback=convergence_measure,
max_history_length=3)
......@@ -159,25 +150,31 @@ if __name__ == "__main__":
S0 = create_power_operator(h_space, power_spectrum=exp(t0),
distribution_strategy=distribution_strategy)
data_power = fft(d).power_analyze()
data_power = log(fft(d).power_analyze(logarithmic=p_space.config["logarithmic"],
nbin=p_space.config["nbin"])**2)
for i in range(100):
S0 = create_power_operator(h_space, power_spectrum=exp(t0),
distribution_strategy=distribution_strategy)
# Initializing the nonlinear Wiener Filter energy
map_energy = NonlinearWienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0)
map_energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0)
# Minimization with chosen minimizer
(map_energy, convergence) = minimizer1(map_energy)
map_energy = map_energy.analytic_solution()
# Updating parameters for correlation structure reconstruction
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=.5, samples=5)
(power_energy, convergence) = minimizer1(power_energy)
power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=.5, samples=3)
if i > 0:
(power_energy, convergence) = minimizer1(power_energy)
else:
(power_energy, convergence) = minimizer1(power_energy)
# Setting new power spectrum
t0 = power_energy.position
plot_parameters(m0,t0,log(data_power**2))
t0.val[-1] = t0.val[-2]
# Plotting current estimate
plot_parameters(m0,t0,data_power)
# Transforming fields to position space for plotting
......@@ -185,6 +182,7 @@ if __name__ == "__main__":
m = fft.adjoint_times(map_energy.position)
# Plotting
d_data = d.val.get_full_data().real
......
......@@ -61,7 +61,7 @@ if __name__ == "__main__":
sp = Field(p_space, val=lambda z: pow_spec(z)**(1./2),
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)
......
......@@ -167,110 +167,20 @@ class Field(Loggable, Versionable, object):
# ---Powerspectral methods---
def power_analyze(self, spaces=None, logarithmic=False, nbin=None, binbounds=None,
decompose_power=False):
# check if all spaces in `self.domain` are either harmonic or
# power_space instances
for sp in self.domain:
if not sp.harmonic and not isinstance(sp, PowerSpace):
self.logger.info(
"Field has a space in `domain` which is neither "
"harmonic nor a PowerSpace.")
# check if the `spaces` input is valid
spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
if spaces is None:
if len(self.domain) == 1:
spaces = (0,)
else:
raise ValueError(
"Field has multiple spaces as domain "
"but `spaces` is None.")
if len(spaces) == 0:
raise ValueError(
"No space for analysis specified.")
elif len(spaces) > 1:
raise ValueError(
"Conversion of only one space at a time is allowed.")
space_index = spaces[0]
if not self.domain[space_index].harmonic:
raise ValueError(
"The analyzed space must be harmonic.")
# Create the target PowerSpace instance:
# If the associated signal-space field was real, we extract the