Commit 2d176740 authored by Martin Reinecke's avatar Martin Reinecke

merge master

parents 90af04d0 f64c8f54
Pipeline #16806 failed with stage
in 13 minutes and 12 seconds
......@@ -23,9 +23,9 @@ def plot_parameters(m, t, p, p_d):
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.Heatmap(z=m)], filename='map.html', auto_open=False)
pl.plot([go.Scatter(x=x, y=t), go.Scatter(x=x, y=p),
go.Scatter(x=x, y=p_d)], filename="t.html")
go.Scatter(x=x, y=p_d)], filename="t.html", auto_open=False)
class AdjointFFTResponse(LinearOperator):
......@@ -106,7 +106,7 @@ if __name__ == "__main__":
data_power = log(fft(d).power_analyze(binbounds=p_space.binbounds))
d_data = d.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=d_data)], filename='data.html')
pl.plot([go.Heatmap(z=d_data)], filename='data.html', auto_open=False)
# Minimization strategy
def convergence_measure(a_energy, iteration): # returns current energy
......
# -*- coding: utf-8 -*-
import numpy as np
import nifty as ift
from nifty import plotting
from keepers import Repository
if __name__ == "__main__":
signal_to_noise = 1.5 # The signal to noise ratioa
# Setting up parameters |\label{code:wf_parameters}|
correlation_length_1 = 1. # Typical distance over which the field is correlated
field_variance_1 = 2. # Variance of field in position space
response_sigma_1 = 0.05 # Smoothing length of response (in same unit as L)
def power_spectrum_1(k): # note: field_variance**2 = a*k_0/4.
a = 4 * correlation_length_1 * field_variance_1**2
return a / (1 + k * correlation_length_1) ** 4.
# Setting up the geometry |\label{code:wf_geometry}|
L_1 = 2. # Total side-length of the domain
N_pixels_1 = 512 # Grid resolution (pixels per axis)
signal_space_1 = ift.RGSpace([N_pixels_1], distances=L_1/N_pixels_1)
harmonic_space_1 = ift.FFTOperator.get_default_codomain(signal_space_1)
fft_1 = ift.FFTOperator(harmonic_space_1, target=signal_space_1,
domain_dtype=np.complex, target_dtype=np.complex)
power_space_1 = ift.PowerSpace(harmonic_space_1, distribution_strategy='fftw')
mock_power_1 = ift.Field(power_space_1, val=power_spectrum_1,
distribution_strategy='not')
# Setting up parameters |\label{code:wf_parameters}|
correlation_length_2 = 1. # Typical distance over which the field is correlated
field_variance_2 = 2. # Variance of field in position space
response_sigma_2 = 0.01 # Smoothing length of response (in same unit as L)
def power_spectrum_2(k): # note: field_variance**2 = a*k_0/4.
a = 4 * correlation_length_2 * field_variance_2**2
return a / (1 + k * correlation_length_2) ** 2.5
# Setting up the geometry |\label{code:wf_geometry}|
L_2 = 2. # Total side-length of the domain
N_pixels_2 = 512 # Grid resolution (pixels per axis)
signal_space_2 = ift.RGSpace([N_pixels_2], distances=L_2/N_pixels_2)
harmonic_space_2 = ift.FFTOperator.get_default_codomain(signal_space_2)
fft_2 = ift.FFTOperator(harmonic_space_2, target=signal_space_2,
domain_dtype=np.complex, target_dtype=np.complex)
power_space_2 = ift.PowerSpace(harmonic_space_2, distribution_strategy='not')
mock_power_2 = ift.Field(power_space_2, val=power_spectrum_2,
distribution_strategy='not')
fft = ift.ComposedOperator((fft_1, fft_2))
mock_power = ift.Field(domain=(power_space_1, power_space_2),
val=np.outer(mock_power_1.val.get_full_data(),
mock_power_2.val.get_full_data()),
distribution_strategy='not')
diagonal = mock_power.power_synthesize(spaces=(0, 1), mean=1, std=0,
real_signal=False,
distribution_strategy='fftw')**2
S = ift.DiagonalOperator(domain=(harmonic_space_1, harmonic_space_2),
diagonal=diagonal)
np.random.seed(10)
mock_signal = fft(mock_power.power_synthesize(real_signal=True,
distribution_strategy='fftw'))
# Setting up a exemplary response
N1_10 = int(N_pixels_1/10)
mask_1 = ift.Field(signal_space_1, val=1., distribution_strategy='fftw')
mask_1.val[N1_10*7:N1_10*9] = 0.
N2_10 = int(N_pixels_2/10)
mask_2 = ift.Field(signal_space_2, val=1., distribution_strategy='not')
mask_2.val[N2_10*7:N2_10*9] = 0.
R = ift.ResponseOperator((signal_space_1, signal_space_2),
sigma=(response_sigma_1, response_sigma_2),
exposure=(mask_1, mask_2)) #|\label{code:wf_response}|
data_domain = R.target
R_harmonic = ift.ComposedOperator([fft, R], default_spaces=(0, 1, 0, 1))
# Setting up the noise covariance and drawing a random noise realization
N = ift.DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise,
bare=True,
distribution_strategy='fftw')
noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise),
mean=0,
distribution_strategy='fftw')
data = R(mock_signal) + noise #|\label{code:wf_mock_data}|
# Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data))
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R_harmonic)
wiener_curvature._InvertibleOperatorMixin__inverter.convergence_tolerance = 1e-3
m_k = wiener_curvature.inverse_times(j) #|\label{code:wf_wiener_filter}|
m = fft(m_k)
# Probing the variance
class Proby(ift.DiagonalProberMixin, ift.Prober): pass
proby = Proby((signal_space_1, signal_space_2), probe_count=100)
proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z))))
# sm = SmoothingOperator(signal_space, sigma=0.02)
# variance = sm(proby.diagonal.weight(-1))
variance = proby.diagonal.weight(-1)
repo = Repository('repo_100.h5')
repo.add(mock_signal, 'mock_signal')
repo.add(data, 'data')
repo.add(m, 'm')
repo.add(variance, 'variance')
repo.commit()
plot_space = ift.RGSpace((N_pixels_1, N_pixels_2))
plotter = plotting.RG2DPlotter(color_map=plotting.colormaps.PlankCmap())
plotter.figure.xaxis = ift.plotting.Axis(label='Pixel Index')
plotter.figure.yaxis = ift.plotting.Axis(label='Pixel Index')
plotter.plot.zmin = 0.
plotter.plot.zmax = 3.
sm = ift.SmoothingOperator(plot_space, sigma=0.03)
plotter(ift.log(ift.sqrt(sm(ift.Field(plot_space, val=variance.val.real)))), path='uncertainty.html')
plotter.plot.zmin = np.real(mock_signal.min());
plotter.plot.zmax = np.real(mock_signal.max());
plotter(ift.Field(plot_space, val=mock_signal.val.real), path='mock_signal.html')
plotter(ift.Field(plot_space, val=data.val.get_full_data().real), path = 'data.html')
plotter(ift.Field(plot_space, val=m.val.real), path = 'map.html')
# -*- coding: utf-8 -*-
import nifty as ift
from nifty import plotting
import numpy as np
from keepers import Repository
if __name__ == "__main__":
ift.nifty_configuration['default_distribution_strategy'] = 'fftw'
# Setting up parameters |\label{code:wf_parameters}|
correlation_length_scale = 1. # Typical distance over which the field is correlated
fluctuation_scale = 2. # Variance of field in position space
response_sigma = 0.05 # Smoothing length of response (in same unit as L)
signal_to_noise = 1.5 # The signal to noise ratio
np.random.seed(43) # Fixing the random seed
def power_spectrum(k): # Defining the power spectrum
a = 4 * correlation_length_scale * fluctuation_scale**2
return a / (1 + (k * correlation_length_scale)**2) ** 2
# Setting up the geometry |\label{code:wf_geometry}|
L = 2. # Total side-length of the domain
N_pixels = 512 # Grid resolution (pixels per axis)
signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
harmonic_space = ift.FFTOperator.get_default_codomain(signal_space)
fft = ift.FFTOperator(harmonic_space, target=signal_space, target_dtype=np.float)
power_space = ift.PowerSpace(harmonic_space)
# Creating the mock signal |\label{code:wf_mock_signal}|
S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum)
mock_power = ift.Field(power_space, val=power_spectrum)
mock_signal = fft(mock_power.power_synthesize(real_signal=True))
# Setting up an exemplary response
mask = ift.Field(signal_space, val=1.)
N10 = int(N_pixels/10)
mask.val[N10*5:N10*9, N10*5:N10*9] = 0.
R = ift.ResponseOperator(signal_space, sigma=(response_sigma,), exposure=(mask,)) #|\label{code:wf_response}|
data_domain = R.target[0]
R_harmonic = ift.ComposedOperator([fft, R], default_spaces=[0, 0])
# Setting up the noise covariance and drawing a random noise realization
N = ift.DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise, bare=True)
noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(mock_signal) + noise #|\label{code:wf_mock_data}|
# Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data))
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R_harmonic)
m_k = wiener_curvature.inverse_times(j) #|\label{code:wf_wiener_filter}|
m = fft(m_k)
# Probing the uncertainty |\label{code:wf_uncertainty_probing}|
class Proby(ift.DiagonalProberMixin, ift.Prober): pass
proby = Proby(signal_space, probe_count=800)
proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z)))) #|\label{code:wf_variance_fft_wrap}|
sm = ift.SmoothingOperator(signal_space, sigma=0.03)
variance = ift.sqrt(sm(proby.diagonal.weight(-1))) #|\label{code:wf_variance_weighting}|
repo = Repository('repo_800.h5')
repo.add(mock_signal, 'mock_signal')
repo.add(data, 'data')
repo.add(m, 'm')
repo.add(variance, 'variance')
repo.commit()
# Plotting #|\label{code:wf_plotting}|
plotter = plotting.RG2DPlotter(color_map=plotting.colormaps.PlankCmap())
plotter.figure.xaxis = ift.plotting.Axis(label='Pixel Index')
plotter.figure.yaxis = ift.plotting.Axis(label='Pixel Index')
plotter.plot.zmax = variance.max(); plotter.plot.zmin = 0
plotter(variance, path = 'uncertainty.html')
plotter.plot.zmax = mock_signal.max(); plotter.plot.zmin = mock_signal.min()
plotter(mock_signal, path='mock_signal.html')
plotter(ift.Field(signal_space, val=data.val), path='data.html')
plotter(m, path='map.html')
......@@ -2,22 +2,23 @@ import numpy as np
from nifty import RGSpace, PowerSpace, Field, FFTOperator, ComposedOperator,\
DiagonalOperator, ResponseOperator, plotting,\
create_power_operator
create_power_operator, nifty_configuration
from nifty.library import WienerFilterCurvature
if __name__ == "__main__":
distribution_strategy = 'not'
nifty_configuration['default_distribution_strategy'] = 'fftw'
nifty_configuration['harmonic_rg_base'] = 'real'
# Setting up variable parameters
# Typical distance over which the field is correlated
correlation_length = 0.01
correlation_length = 0.05
# Variance of field in position space sqrt(<|s_x|^2>)
field_variance = 2.
# smoothing length of response (in same unit as L)
response_sigma = 0.1
response_sigma = 0.01
# The signal to noise ratio
signal_to_noise = 0.7
......@@ -36,19 +37,17 @@ if __name__ == "__main__":
signal_space = RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
harmonic_space = FFTOperator.get_default_codomain(signal_space)
fft = FFTOperator(harmonic_space, target=signal_space,
domain_dtype=np.complex, target_dtype=np.float)
power_space = PowerSpace(harmonic_space,
distribution_strategy=distribution_strategy)
fft = FFTOperator(harmonic_space, target=signal_space)
power_space = PowerSpace(harmonic_space)
# Creating the mock data
S = create_power_operator(harmonic_space, power_spectrum=power_spectrum,
distribution_strategy=distribution_strategy)
S = create_power_operator(harmonic_space, power_spectrum=power_spectrum)
mock_power = Field(power_space, val=power_spectrum,
distribution_strategy=distribution_strategy)
mock_power = Field(power_space, val=power_spectrum)
np.random.seed(43)
mock_harmonic = mock_power.power_synthesize(real_signal=True)
if nifty_configuration['harmonic_rg_base'] == 'real':
mock_harmonic = mock_harmonic.real
mock_signal = fft(mock_harmonic)
R = ResponseOperator(signal_space, sigma=(response_sigma,))
......@@ -73,9 +72,11 @@ if __name__ == "__main__":
m_s = fft(m)
plotter = plotting.RG2DPlotter()
plotter.title = 'mock_signal.html';
plotter(mock_signal)
plotter.title = 'data.html'
plotter(Field(signal_space,
val=data.val.get_full_data().reshape(signal_space.shape)))
plotter.title = 'map.html'; plotter(m_s)
\ No newline at end of file
plotter.path = 'mock_signal.html'
plotter(mock_signal.real)
plotter.path = 'data.html'
plotter(Field(
signal_space,
val=data.val.get_full_data().real.reshape(signal_space.shape)))
plotter.path = 'map.html'
plotter(m_s.real)
......@@ -10,6 +10,7 @@ rank = comm.rank
np.random.seed(42)
class AdjointFFTResponse(LinearOperator):
def __init__(self, FFT, R, default_spaces=None):
super(AdjointFFTResponse, self).__init__(default_spaces)
......@@ -23,6 +24,7 @@ class AdjointFFTResponse(LinearOperator):
def _adjoint_times(self, x, spaces=None):
return self.FFT(self.R.adjoint_times(x))
@property
def domain(self):
return self._domain
......@@ -36,13 +38,12 @@ class AdjointFFTResponse(LinearOperator):
return False
if __name__ == "__main__":
distribution_strategy = 'not'
# Set up position space
s_space = RGSpace([128,128])
s_space = RGSpace([128, 128])
# s_space = HPSpace(32)
# Define harmonic transformation and associated harmonic space
......@@ -52,7 +53,8 @@ if __name__ == "__main__":
# Setting up power space
p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
# Choosing the prior correlation structure and defining correlation operator
# 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,
......@@ -69,7 +71,7 @@ if __name__ == "__main__":
Instrument = DiagonalOperator(s_space, diagonal=1.)
# Instrument._diagonal.val[200:400, 200:400] = 0
#Adding a harmonic transformation to the instrument
# 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)
......@@ -84,9 +86,9 @@ if __name__ == "__main__":
# Choosing the minimization strategy
def convergence_measure(energy, iteration): # returns current energy
def convergence_measure(energy, iteration): # returns current energy
x = energy.value
print (x, iteration)
print(x, iteration)
# minimizer = SteepestDescent(convergence_tolerance=0,
# iteration_limit=50,
......@@ -109,20 +111,19 @@ if __name__ == "__main__":
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)
energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S)
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)
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 = sugar.generate_posterior_sample(m0, D0)
sample_variance += sample**2
sample_mean += sample
variance = sample_variance/n_samples - (sample_mean/n_samples)
......@@ -23,7 +23,7 @@ from nifty.field import Field
__all__ = ['cos', 'sin', 'cosh', 'sinh', 'tan', 'tanh', 'arccos', 'arcsin',
'arccosh', 'arcsinh', 'arctan', 'arctanh', 'sqrt', 'exp', 'log',
'conjugate', 'clipped_exp']
'conjugate', 'clipped_exp', 'limited_exp']
def _math_helper(x, function):
......@@ -99,6 +99,19 @@ def clipped_exp(x):
return _math_helper(x, lambda z: np.exp(np.minimum(200, z)))
def limited_exp(x):
thr = 200
expthr = np.exp(thr)
return _math_helper(x, lambda z: _limited_exp_helper(z, thr, expthr))
def _limited_exp_helper(x, thr, expthr):
mask = (x > thr)
result = np.exp(x)
result[mask] = ((1-thr) + x[mask])*expthr
return result
def log(x, base=None):
result = _math_helper(x, np.log)
if base is not None:
......
......@@ -70,11 +70,18 @@ variable_default_distribution_strategy = keepers.Variable(
if z == 'fftw' else True),
genus='str')
variable_harmonic_rg_base = keepers.Variable(
'harmonic_rg_base',
['real', 'complex'],
lambda z: z in ['real', 'complex'],
genus='str')
nifty_configuration = keepers.get_Configuration(
name='NIFTy',
variables=[variable_fft_module,
variable_default_field_dtype,
variable_default_distribution_strategy],
variable_default_distribution_strategy,
variable_harmonic_rg_base],
file_name='NIFTy.conf',
search_paths=[os.path.expanduser('~') + "/.config/nifty/",
os.path.expanduser('~') + "/.config/",
......
......@@ -66,12 +66,9 @@ class Energy(Loggable, object):
__metaclass__ = NiftyMeta
def __init__(self, position):
super(Energy, self).__init__()
self._cache = {}
try:
position = position.copy()
except AttributeError:
pass
self._position = position
self._position = position.copy()
def at(self, position):
""" Initializes and returns a new Energy object at the new position.
......
......@@ -68,12 +68,16 @@ class LineEnergy(object):
"""
def __init__(self, line_position, energy, line_direction, offset=0.):
super(LineEnergy, self).__init__()
self._line_position = float(line_position)
self._line_direction = line_direction
pos = energy.position \
+ (self._line_position-float(offset))*self._line_direction
self.energy = energy.at(position=pos)
if self._line_position==float(offset):
self.energy = energy
else:
pos = energy.position \
+ (self._line_position-float(offset))*self._line_direction
self.energy = energy.at(position=pos)
def at(self, line_position):
""" Returns LineEnergy at new position, memorizing the zero point.
......
......@@ -19,7 +19,6 @@
from __future__ import division
import ast
import itertools
import numpy as np
from keepers import Versionable,\
......@@ -406,7 +405,6 @@ class Field(Loggable, Versionable, object):
distribution_strategy=distribution_strategy,
logarithmic=logarithmic, nbin=nbin,
binbounds=binbounds)
power_spectrum = cls._calculate_power_spectrum(
field_val=work_field.val,
pdomain=power_domain,
......@@ -437,6 +435,7 @@ class Field(Loggable, Versionable, object):
target_shape=field_val.shape,
target_strategy=field_val.distribution_strategy,
axes=axes)
power_spectrum = pindex.bincount(weights=field_val,
axis=axes)
rho = pdomain.rho
......@@ -462,14 +461,14 @@ class Field(Loggable, Versionable, object):
"A slicing distributor shall not be reshaped to "
"something non-sliced.")
semiscaled_shape = [1, ] * len(target_shape)
for i in axes:
semiscaled_shape[i] = target_shape[i]
semiscaled_local_shape = [1, ] * len(target_shape)
for i in range(len(axes)):
semiscaled_local_shape[axes[i]] = pindex.local_shape[i]
local_data = pindex.get_local_data(copy=False)
semiscaled_local_data = local_data.reshape(semiscaled_shape)
semiscaled_local_data = local_data.reshape(semiscaled_local_shape)
result_obj = pindex.copy_empty(global_shape=target_shape,
distribution_strategy=target_strategy)
result_obj.set_full_data(semiscaled_local_data, copy=False)
result_obj.data[:] = semiscaled_local_data
return result_obj
......@@ -478,7 +477,7 @@ class Field(Loggable, Versionable, object):
""" Yields a sampled field with `self`**2 as its power spectrum.
This method draws a Gaussian random field in the harmonic partner
domain of this fields domains, using this field as power spectrum.
domain of this field's domains, using this field as power spectrum.
Parameters
----------
......@@ -697,7 +696,7 @@ class Field(Loggable, Versionable, object):
# ---Properties---
def set_val(self, new_val=None, copy=False):
""" Sets the fields distributed_data_object.
""" Sets the field's distributed_data_object.
Parameters
----------
......@@ -870,7 +869,7 @@ class Field(Loggable, Versionable, object):
dtype : type
The datatype the output shall have. This can be used to override
the fields dtype.
the field's dtype.
Returns
-------
......
......@@ -7,7 +7,7 @@ class WienerFilterEnergy(Energy):
"""The Energy for the Wiener filter.
It covers the case of linear measurement with
Gaussian noise and Gaussain signal prior with known covariance.
Gaussian noise and Gaussian signal prior with known covariance.
Parameters
----------
......
......@@ -123,7 +123,6 @@ class DescentMinimizer(Loggable, object):
convergence = 0
f_k_minus_1 = None
step_length = 0
iteration_number = 1
while True:
......@@ -150,7 +149,7 @@ class DescentMinimizer(Loggable, object):
# compute the step length, which minimizes energy.value along the
# search direction
try:
step_length, f_k, new_energy = \
new_energy = \
self.line_searcher.perform_line_search(
energy=energy,
pk=descent_direction,
......@@ -160,12 +159,10 @@ class DescentMinimizer(Loggable, object):
"Stopping because of RuntimeError in line-search")
break
if f_k_minus_1 is None:
delta = 1e30
else:
delta = (abs(f_k-f_k_minus_1) /
max(abs(f_k), abs(f_k_minus_1), 1.))
f_k_minus_1 = energy.value
f_k = new_energy.value
delta = (abs(f_k-f_k_minus_1) /
max(abs(f_k), abs(f_k_minus_1), 1.))
# check if new energy value is bigger than old energy value
if (new_energy.value - energy.value) > 0:
self.logger.info("Line search algorithm returned a new energy "
......@@ -174,9 +171,9 @@ class DescentMinimizer(Loggable, object):
energy = new_energy
# check convergence
self.logger.debug("Iteration:%08u step_length=%3.1E "
self.logger.debug("Iteration:%08u "
"delta=%3.1E energy=%3.1E" %
(iteration_number, step_length, delta,
(iteration_number, delta,
energy.value))
if delta == 0:
convergence = self.convergence_level + 2
......
......@@ -34,8 +34,6 @@ class LineSearch(Loggable, object):
----------
line_energy : LineEnergy Object
LineEnergy object from which we can extract energy at a specific point.
f_k_minus_1 : Field
Value of the field at the k-1 iteration of the line search procedure.
preferred_initial_step_size : float
Initial guess for the step length.
......@@ -44,32 +42,8 @@ class LineSearch(Loggable, object):
__metaclass__ = abc.ABCMeta
def __init__(self):
self.line_energy = None
self.f_k_minus_1 = None
self.preferred_initial_step_size = None
def _set_line_energy(self, energy, pk, f_k_minus_1=None):
"""Set the coordinates for a new line search.
Parameters
----------
energy : Energy object
Energy object from which we can calculate the energy, gradient and
curvature at a specific point.
pk : Field
Unit vector pointing into the search direction.
f_k_minus_1 : float
Value of the fuction (energy) which will be minimized at the k-1
iteration of the line search procedure. (Default: None)
"""
self.line_energy = LineEnergy(line_position=0.,
energy=energy,
line_direction=pk)
if f_k_minus_1 is not None:
f_k_minus_1 = f_k_minus_1.copy()
self.f_k_minus_1 = f_k_minus_1
@abc.abstractmethod
def perform_line_search(self, energy, pk, f_k_minus_1=None):
raise NotImplementedError
......@@ -19,12 +19,13 @@
import numpy as np
from .line_search import LineSearch
from ...energies import LineEnergy
class LineSearchStrongWolfe(LineSearch):
"""Class for finding a step size that satisfies the strong Wolfe conditions.
Algorithm contains two stages. It begins whit a trial step length and
Algorithm contains two stages. It begins with a trial step length and
keeps increasing it until it finds an acceptable step length or an
interval. If it does not satisfy the Wolfe conditions, it performs the Zoom
algorithm (second stage). By interpolating it decreases the size of the
......@@ -77,8 +78,8 @@ class LineSearchStrongWolfe(LineSearch):
"""Performs the first stage of the algorithm.
It starts with a trial step size and it keeps increasing it until it
satisfy the strong Wolf conditions. It also performs the descent and
returns the optimal step length and the new enrgy.
satisfies the strong Wolf conditions. It also performs the descent and
returns the optimal step length and the new energy.
Parameters
----------
......@@ -86,29 +87,22 @@ class LineSearchStrongWolfe(LineSearch):
Energy object from which we will calculate the energy and the
gradient at a specific point.
pk : Field
Unit vector pointing into the search direction.
Vector pointing into the search direction.
f_k_minus_1 : float
Value of the fuction (which is being minimized) at the k-1
iteration of the line search procedure. (Default: None)
Returns
-------
alpha_star : float
The optimal step length in the descent direction.
phi_star : float
Value of the energy after the performed descent.
energy_star : Energy object
The new Energy object on the new position.
"""
self._set_line_energy(energy, pk, f_k_minus_1=f_k_minus_1)
max_step_size = self.max_step_size
max_iterations = self.max_iterations
le_0 = LineEnergy(0., energy, pk, 0.)