Commit c55e5b61 authored by Jakob Knollmueller's avatar Jakob Knollmueller

current status for Daniel

parent 222f8620
Pipeline #13396 passed with stage
in 5 minutes and 20 seconds
......@@ -64,30 +64,30 @@ if __name__ == "__main__":
# Setting up power space
p_space = PowerSpace(h_space, logarithmic=False,
distribution_strategy=distribution_strategy, nbin=500)
distribution_strategy=distribution_strategy)#, nbin=5)
# Choosing the prior correlation structure and defining correlation operator
p_spec = (lambda k: (.05 / (k + 1) ** 2))
p_spec = (lambda k: (.05 / (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=lambda z: pow_spec(z)**(1./2),
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 = 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
#Instrument._diagonal.val[64:512-64, 64:512-64] = 0
#Adding a harmonic transformation to the instrument
R = AdjointFFTResponse(fft, Instrument)
noise = 1.
noise = .1
N = DiagonalOperator(s_space, diagonal=noise, bare=True)
n = Field.from_random(domain=s_space,
random_type='normal',
......@@ -98,9 +98,9 @@ if __name__ == "__main__":
d = R(sh) + n
realized_power = log(sh.power_analyze(logarithmic=p_space.config["logarithmic"],
nbin=p_space.config["nbin"])**2)
nbin=p_space.config["nbin"]))
data_power = log(fft(d).power_analyze(logarithmic=p_space.config["logarithmic"],
nbin=p_space.config["nbin"])**2)
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')
......@@ -117,9 +117,9 @@ if __name__ == "__main__":
iteration_limit=3,
callback=convergence_measure)
minimizer2 = VL_BFGS(convergence_tolerance=0,
iteration_limit=50,
iteration_limit=7,
callback=convergence_measure,
max_history_length=10)
max_history_length=3)
# Setting starting position
flat_power = Field(p_space,val=10e-8)
......@@ -127,12 +127,13 @@ if __name__ == "__main__":
# 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(0.5*t0),
S0 = create_power_operator(h_space, power_spectrum=exp(t0),
distribution_strategy=distribution_strategy)
for i in range(500):
S0 = create_power_operator(h_space, power_spectrum=exp(0.5*t0),
S0 = create_power_operator(h_space, power_spectrum=exp(t0),
distribution_strategy=distribution_strategy)
# Initializing the nonlinear Wiener Filter energy
......@@ -144,7 +145,7 @@ 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=1000., samples=1)
power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=100., samples=5)
(power_energy, convergence) = minimizer1(power_energy)
......@@ -152,7 +153,7 @@ if __name__ == "__main__":
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)
plot_parameters(m0,t0,log(sp),realized_power, data_power)
# Transforming fields to position space for plotting
......@@ -183,11 +184,3 @@ if __name__ == "__main__":
if rank == 0:
pl.plot([go.Heatmap(z=m_data)], filename='map.html')
collector = fft(n).power_analyze() **2
for i in range(10000):
n = Field.from_random(domain=s_space,
random_type='normal',
std=sqrt(noise),
mean=0)
collector = fft(n).power_analyze()**2
collector /= 10001.
......@@ -48,12 +48,13 @@ class NonlinearResponse(LinearOperator):
def plot_parameters(m,t, t_real):
m = fft.adjoint_times(m)
x = log(t.domain[0].kindex[1:])
m_data = m.val.get_full_data().real
t_data = t.val.get_full_data().real
t_real_data = t_real.val.get_full_data().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")
pl.plot([go.Scatter(x = x,y=t_data),
go.Scatter(x= x, y=t_real_data[1:])], filename="t.html")
class AdjointFFTResponse(LinearOperator):
def __init__(self, FFT, R, default_spaces=None):
super(AdjointFFTResponse, self).__init__(default_spaces)
......@@ -83,7 +84,7 @@ if __name__ == "__main__":
distribution_strategy = 'not'
full_data = np.genfromtxt("train_data.csv", delimiter = ',')
d = full_data.T[3]
d = full_data.T[2]
d[0] = 0.
d -= d.mean()
d[0] = 0.
......@@ -97,7 +98,7 @@ if __name__ == "__main__":
fft = FFTOperator(s_space)
h_space = fft.target[0]
p_space = PowerSpace(h_space, logarithmic=False,
distribution_strategy=distribution_strategy)
distribution_strategy=distribution_strategy)#, nbin=50)
# Choosing the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.01)
......@@ -133,10 +134,10 @@ if __name__ == "__main__":
iteration_limit=30,
callback=convergence_measure)
minimizer1 = VL_BFGS(convergence_tolerance=0,
iteration_limit=30,
minimizer3 = VL_BFGS(convergence_tolerance=0,
iteration_limit=50,
callback=convergence_measure,
max_history_length=3)
max_history_length=10)
......@@ -144,17 +145,17 @@ if __name__ == "__main__":
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=-10)
t0 = Field(p_space, val=log(1./(1+p_space.kindex)**2))
t0 = Field(p_space,val=-13.)
# t0 = log(sp.copy()**2)
S0 = create_power_operator(h_space, power_spectrum=exp(t0),
distribution_strategy=distribution_strategy)
data_power = log(fft(d).power_analyze(logarithmic=p_space.config["logarithmic"],
nbin=p_space.config["nbin"])**2)
for i in range(100):
nbin=p_space.config["nbin"]))
for i in range(500):
S0 = create_power_operator(h_space, power_spectrum=exp(t0),
distribution_strategy=distribution_strategy)
distribution_strategy=distribution_strategy)
# Initializing the nonlinear Wiener Filter energy
map_energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0)
......@@ -165,16 +166,15 @@ 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=.5, samples=3)
if i > 0:
(power_energy, convergence) = minimizer1(power_energy)
else:
(power_energy, convergence) = minimizer1(power_energy)
power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=100000., samples=20)
(power_energy, convergence) = minimizer3(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,data_power)
plot_parameters(m0, t0, data_power)
# Transforming fields to position space for plotting
......
......@@ -6,24 +6,47 @@ from nifty.sugar import generate_posterior_sample
from nifty import Field, exp
class CriticalPowerEnergy(Energy):
"""The Energy for the Gaussian lognormal case.
"""The Energy of the power spectrum according to the critical filter.
It describes the situation of linear measurement of a
lognormal signal with Gaussian noise and Gaussain signal prior.
It describes the energy of the logarithmic amplitudes of the power spectrum of
a Gaussian random field under reconstruction uncertainty with smoothness and
inverse gamma prior. It is used to infer a correlated signal with unknown correlation
structure.
Parameters
----------
d : Field,
the data.
R : Operator,
The nonlinear response operator, describtion of the measurement process.
N : EndomorphicOperator,
The noise covariance in data space.
S : EndomorphicOperator,
The prior signal covariance in harmonic space.
position : Field,
The current position of this energy.
m : Field,
The map whichs power spectrum has to be inferred
D : EndomorphicOperator,
The curvature of the Gaussian encoding the posterior covariance.
If not specified, the map is assumed to be no reconstruction.
default : None
alpha : float
The spectral prior of the inverse gamma distribution. 1.0 corresponds to
non-informative.
default : 1.0
q : float
The cutoff parameter of the inverse gamma distribution. 0.0 corresponds to
non-informative.
default : 0.0
sigma : float
The parameter of the smoothness prior.
default : ??? None? ???????
logarithmic : boolean
Whether smoothness acts on linear or logarithmic scale.
samples : integer
Number of samples used for the estimation of the uncertainty corrections.
default : 3
w : Field
The contribution from the map with or without uncertainty. It is used
to pass on the result of the costly sampling during the minimization.
default : None
"""
def __init__(self, position, m, D=None, alpha =1.0, q=0., sigma=0., w=None, samples=3):
def __init__(self, position, m, D=None, alpha =1.0, q=0., sigma=0.,
logarithmic = True, samples=3, w=None):
super(CriticalPowerEnergy, self).__init__(position = position)
self.m = m
self.D = D
......@@ -31,7 +54,8 @@ 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.T = SmoothnessOperator(domain=self.position.domain[0], sigma=self.sigma,
logarithmic=logarithmic)
self.rho = self.position.domain[0].rho
self.w = w
if self.w is None:
......@@ -72,16 +96,13 @@ class CriticalPowerEnergy(Energy):
posterior_sample = generate_posterior_sample(m, D)
projected_sample = posterior_sample.power_analyze(
logarithmic=self.position.domain[0].config["logarithmic"],
nbin= self.position.domain[0].config["nbin"],
decompose_power=False)
w += (projected_sample **2) * self.rho
nbin= self.position.domain[0].config["nbin"])
w += (projected_sample) * self.rho
w /= float(samples)
else:
w = m.power_analyze(
logarithmic=self.position.domain[0].config["logarithmic"],
nbin=self.position.domain[0].config["nbin"],
decompose_power=False)
w **= 2
nbin=self.position.domain[0].config["nbin"])
w *= self.rho
return w
......
......@@ -5,10 +5,12 @@ class WienerFilterEnergy(Energy):
"""The Energy for the Wiener filter.
It describes the situation of linear measurement with
Gaussian noise and Gaussain signal prior.
Gaussian noise and Gaussain signal prior with known covariance.
Parameters
----------
position: Field,
The current position.
d : Field,
the data.
R : Operator,
......
from nifty.operators import EndomorphicOperator,\
InvertibleOperatorMixin
InvertibleOperatorMixin,\
DiagonalOperator
class CriticalPowerCurvature(InvertibleOperatorMixin, EndomorphicOperator):
def __init__(self, theta, T, inverter=None, preconditioner=None):
self.theta = theta
self.theta = DiagonalOperator(theta.domain, diagonal=theta)
self.T = T
# if preconditioner is None:
# preconditioner = self.T.times
if preconditioner is None:
preconditioner = self.theta.inverse_times
self._domain = self.theta.domain
super(CriticalPowerCurvature, self).__init__(inverter=inverter,
preconditioner=preconditioner)
......@@ -27,4 +28,4 @@ class CriticalPowerCurvature(InvertibleOperatorMixin, EndomorphicOperator):
# ---Added properties and methods---
def _times(self, x, spaces):
return self.T(x) + self.theta * x
return self.T(x) + self.theta(x)
......@@ -70,7 +70,7 @@ def generate_posterior_sample(mean, covariance):
S = covariance.S
R = covariance.R
N = covariance.N
power = sqrt(S.diagonal().power_analyze())
power = S.diagonal().power_analyze()**.5
mock_signal = power.power_synthesize(real_signal=True)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment