Commit fd1af103 authored by Martin Reinecke's avatar Martin Reinecke

intermediate

parent 7c50e35f
...@@ -102,8 +102,8 @@ if __name__ == "__main__": ...@@ -102,8 +102,8 @@ if __name__ == "__main__":
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"])) nbin=p_space.config["nbin"]))
d_data = d.val.get_full_data().real d_data = d.val.get_full_data().real
if rank == 0: #if rank == 0:
pl.plot([go.Heatmap(z=d_data)], filename='data.html') # pl.plot([go.Heatmap(z=d_data)], filename='data.html')
# minimization strategy # minimization strategy
...@@ -150,7 +150,7 @@ if __name__ == "__main__": ...@@ -150,7 +150,7 @@ if __name__ == "__main__":
# Plotting current estimate # Plotting current estimate
print i print i
if i%50 == 0: #if i%50 == 0:
plot_parameters(m0,t0,log(sp), data_power) # plot_parameters(m0,t0,log(sp), data_power)
...@@ -16,10 +16,8 @@ ...@@ -16,10 +16,8 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes. # and financially supported by the Studienstiftung des deutschen Volkes.
from .energy import Energy
class LineEnergy:
class LineEnergy(Energy):
""" Evaluates an underlying Energy along a certain line direction. """ Evaluates an underlying Energy along a certain line direction.
Given an Energy class and a line direction, its position is parametrized by Given an Energy class and a line direction, its position is parametrized by
...@@ -27,34 +25,31 @@ class LineEnergy(Energy): ...@@ -27,34 +25,31 @@ class LineEnergy(Energy):
Parameters Parameters
---------- ----------
position : float linepos : float
The step length parameter along the given line direction. Defines the full spatial position of this energy via
self.energy.position = zero_point + linepos*line_direction
energy : Energy energy : Energy
The Energy object which will be evaluated along the given direction. The Energy object which will be evaluated along the given direction.
line_direction : Field linedir : Field
Direction used for line evaluation. Direction used for line evaluation. Does not have to be normalized.
zero_point : Field *optional* offset : float *optional*
Fixing the zero point of the line restriction. Used to memorize this Indirectly defines the zero point of the line via the equation
position in new initializations. By the default the current position energy.position = zero_point + offset*line_direction
of the supplied `energy` instance is used (default : None). (default : 0.).
Attributes Attributes
---------- ----------
position : float linepos : float
The position along the given line direction relative to the zero point. The position along the given line direction relative to the zero point.
value : float value : float
The value of the energy functional at given `position`. The value of the energy functional at the given position
gradient : float dd : float
The gradient of the underlying energy instance along the line direction The directional derivative at the given position
projected on the line direction. linedir : Field
curvature : callable
A positive semi-definite operator or function describing the curvature
of the potential at given `position`.
line_direction : Field
Direction along which the movement is restricted. Does not have to be Direction along which the movement is restricted. Does not have to be
normalized. normalized.
energy : Energy energy : Energy
The underlying Energy at the `position` along the line direction. The underlying Energy at the given position
Raises Raises
------ ------
...@@ -72,45 +67,44 @@ class LineEnergy(Energy): ...@@ -72,45 +67,44 @@ class LineEnergy(Energy):
""" """
def __init__(self, position, energy, line_direction, zero_point=None): def __init__(self, linepos, energy, linedir, offset=0.):
super(LineEnergy, self).__init__(position=position) self._linepos = float(linepos)
self.line_direction = line_direction self._linedir = linedir
if zero_point is None:
zero_point = energy.position
self._zero_point = zero_point
position_on_line = self._zero_point + self.position*line_direction pos = energy.position + (self._linepos-float(offset))*self._linedir
self.energy = energy.at(position=position_on_line) self.energy = energy.at(position=pos)
def at(self, position): def at(self, linepos):
""" Returns LineEnergy at new position, memorizing the zero point. """ Returns LineEnergy at new position, memorizing the zero point.
Parameters Parameters
---------- ----------
position : float linepos : float
Parameter for the new position on the line direction. Parameter for the new position on the line direction.
Returns Returns
------- -------
out : LineEnergy
LineEnergy object at new position with same zero point as `self`. LineEnergy object at new position with same zero point as `self`.
""" """
return self.__class__(position, return self.__class__(linepos,
self.energy, self.energy,
self.line_direction, self.linedir,
zero_point=self._zero_point) offset=self.linepos)
@property @property
def value(self): def value(self):
return self.energy.value return self.energy.value
@property @property
def gradient(self): def linepos(self):
return self.energy.gradient.vdot(self.line_direction) return self._linepos
@property
def linedir(self):
return self._linedir
@property @property
def curvature(self): def dd(self):
return self.energy.curvature return self.energy.gradient.vdot(self.linedir)
...@@ -181,7 +181,7 @@ class Field(Loggable, Versionable, object): ...@@ -181,7 +181,7 @@ class Field(Loggable, Versionable, object):
elif isinstance(val, Field): elif isinstance(val, Field):
distribution_strategy = val.distribution_strategy distribution_strategy = val.distribution_strategy
else: else:
self.logger.debug("distribution_strategy set to default!") #self.logger.debug("distribution_strategy set to default!")
distribution_strategy = gc['default_distribution_strategy'] distribution_strategy = gc['default_distribution_strategy']
elif distribution_strategy not in DISTRIBUTION_STRATEGIES['global']: elif distribution_strategy not in DISTRIBUTION_STRATEGIES['global']:
raise ValueError( raise ValueError(
......
...@@ -121,12 +121,18 @@ class DescentMinimizer(Loggable, object): ...@@ -121,12 +121,18 @@ class DescentMinimizer(Loggable, object):
""" """
print "into line search:"
print " pos: ",energy.position.val[0]
print " ene: ",energy.value
convergence = 0 convergence = 0
f_k_minus_1 = None f_k_minus_1 = None
step_length = 0 step_length = 0
iteration_number = 1 iteration_number = 1
while True: while True:
print "line search next iteration:"
print " pos: ",energy.position.val[0]
print " ene: ",energy.value
if self.callback is not None: if self.callback is not None:
try: try:
self.callback(energy, iteration_number) self.callback(energy, iteration_number)
...@@ -147,7 +153,7 @@ class DescentMinimizer(Loggable, object): ...@@ -147,7 +153,7 @@ class DescentMinimizer(Loggable, object):
# current position is encoded in energy object # current position is encoded in energy object
descent_direction = self.get_descent_direction(energy) descent_direction = self.get_descent_direction(energy)
print "descent direction:",descent_direction.val[0]
# compute the step length, which minimizes energy.value along the # compute the step length, which minimizes energy.value along the
# search direction # search direction
step_length, f_k, new_energy = \ step_length, f_k, new_energy = \
...@@ -155,10 +161,20 @@ class DescentMinimizer(Loggable, object): ...@@ -155,10 +161,20 @@ class DescentMinimizer(Loggable, object):
energy=energy, energy=energy,
pk=descent_direction, pk=descent_direction,
f_k_minus_1=f_k_minus_1) f_k_minus_1=f_k_minus_1)
print "out of wolfe:"
print " old pos: ",energy.position.val[0]
print " old ene: ",energy.value
print " new pos: ",new_energy.position.val[0]
print " new ene: ",new_energy.value
print " f_k: ",f_k
f_k_minus_1 = energy.value f_k_minus_1 = energy.value
print " step length: ", step_length
tx1=energy.position-new_energy.position
print " step length 2: ", (energy.position-new_energy.position).norm()
print " step length 3: ", new_energy.position.val[0]-energy.position.val[0]
# check if new energy value is bigger than old energy value # check if new energy value is bigger than old energy value
if (new_energy.value - energy.value) > 0: if (new_energy.value - energy.value) > 0:
print "Line search algorithm returned a new energy that was larger than the old one. Stopping."
self.logger.info("Line search algorithm returned a new energy " self.logger.info("Line search algorithm returned a new energy "
"that was larger than the old one. Stopping.") "that was larger than the old one. Stopping.")
break break
......
...@@ -66,9 +66,9 @@ class LineSearch(Loggable, object): ...@@ -66,9 +66,9 @@ class LineSearch(Loggable, object):
iteration of the line search procedure. (Default: None) iteration of the line search procedure. (Default: None)
""" """
self.line_energy = LineEnergy(position=0., self.line_energy = LineEnergy(linepos=0.,
energy=energy, energy=energy,
line_direction=pk) linedir=pk)
if f_k_minus_1 is not None: if f_k_minus_1 is not None:
f_k_minus_1 = f_k_minus_1.copy() f_k_minus_1 = f_k_minus_1.copy()
self.f_k_minus_1 = f_k_minus_1 self.f_k_minus_1 = f_k_minus_1
......
...@@ -103,16 +103,14 @@ class LineSearchStrongWolfe(LineSearch): ...@@ -103,16 +103,14 @@ class LineSearchStrongWolfe(LineSearch):
""" """
self._set_line_energy(energy, pk, f_k_minus_1=f_k_minus_1) self._set_line_energy(energy, pk, f_k_minus_1=f_k_minus_1)
c1 = self.c1
c2 = self.c2
max_step_size = self.max_step_size max_step_size = self.max_step_size
max_iterations = self.max_iterations max_iterations = self.max_iterations
# initialize the zero phis # initialize the zero phis
old_phi_0 = self.f_k_minus_1 old_phi_0 = self.f_k_minus_1
energy_0 = self.line_energy.at(0) le_0 = self.line_energy.at(0)
phi_0 = energy_0.value phi_0 = le_0.value
phiprime_0 = energy_0.gradient phiprime_0 = le_0.dd
if phiprime_0 == 0: if phiprime_0 == 0:
self.logger.warn("Flat gradient in search direction.") self.logger.warn("Flat gradient in search direction.")
...@@ -135,45 +133,52 @@ class LineSearchStrongWolfe(LineSearch): ...@@ -135,45 +133,52 @@ class LineSearchStrongWolfe(LineSearch):
# start the minimization loop # start the minimization loop
for i in xrange(max_iterations): for i in xrange(max_iterations):
energy_alpha1 = self.line_energy.at(alpha1) print "a0a1:",alpha0, alpha1
phi_alpha1 = energy_alpha1.value print "line search outer iteration", i
le_alpha1 = self.line_energy.at(alpha1)
print "position:", le_alpha1.energy.position.val[0]
phi_alpha1 = le_alpha1.value
print "energy:", le_alpha1.value
if alpha1 == 0: if alpha1 == 0:
self.logger.warn("Increment size became 0.") self.logger.warn("Increment size became 0.")
alpha_star = 0. alpha_star = 0.
phi_star = phi_0 phi_star = phi_0
energy_star = energy_0 le_star = le_0
break break
if (phi_alpha1 > phi_0 + c1*alpha1*phiprime_0) or \ if (phi_alpha1 > phi_0 + self.c1*alpha1*phiprime_0) or \
((phi_alpha1 >= phi_alpha0) and (i > 1)): ((phi_alpha1 >= phi_alpha0) and (i > 1)):
(alpha_star, phi_star, energy_star) = self._zoom( (alpha_star, phi_star, le_star) = self._zoom(
alpha0, alpha1, alpha0, alpha1,
phi_0, phiprime_0, phi_0, phiprime_0,
phi_alpha0, phi_alpha0,
phiprime_alpha0, phiprime_alpha0,
phi_alpha1, phi_alpha1)
c1, c2)
break break
phiprime_alpha1 = energy_alpha1.gradient phiprime_alpha1 = le_alpha1.dd
if abs(phiprime_alpha1) <= -c2*phiprime_0: if abs(phiprime_alpha1) <= -self.c2*phiprime_0:
alpha_star = alpha1 alpha_star = alpha1
phi_star = phi_alpha1 phi_star = phi_alpha1
energy_star = energy_alpha1 le_star = le_alpha1
break break
if phiprime_alpha1 >= 0: if phiprime_alpha1 >= 0:
(alpha_star, phi_star, energy_star) = self._zoom( (alpha_star, phi_star, le_star) = self._zoom(
alpha1, alpha0, alpha1, alpha0,
phi_0, phiprime_0, phi_0, phiprime_0,
phi_alpha1, phi_alpha1,
phiprime_alpha1, phiprime_alpha1,
phi_alpha0, phi_alpha0)
c1, c2)
break break
# update alphas # update alphas
alpha0, alpha1 = alpha1, min(2*alpha1, max_step_size) alpha0, alpha1 = alpha1, min(2*alpha1, max_step_size)
if alpha1 == max_step_size:
alpha_star = alpha1
phi_star = phi_alpha1
le_star = le_alpha1
break
phi_alpha0 = phi_alpha1 phi_alpha0 = phi_alpha1
phiprime_alpha0 = phiprime_alpha1 phiprime_alpha0 = phiprime_alpha1
# phi_alpha1 = self._phi(alpha1) # phi_alpha1 = self._phi(alpha1)
...@@ -182,17 +187,17 @@ class LineSearchStrongWolfe(LineSearch): ...@@ -182,17 +187,17 @@ class LineSearchStrongWolfe(LineSearch):
# max_iterations was reached # max_iterations was reached
alpha_star = alpha1 alpha_star = alpha1
phi_star = phi_alpha1 phi_star = phi_alpha1
energy_star = energy_alpha1 le_star = le_alpha1
self.logger.error("The line search algorithm did not converge.") self.logger.error("The line search algorithm did not converge.")
# extract the full energy from the line_energy # extract the full energy from the line_energy
energy_star = energy_star.energy energy_star = le_star.energy
direction_length = pk.norm() direction_length = pk.norm()
step_length = alpha_star * direction_length step_length = alpha_star * direction_length
return step_length, phi_star, energy_star return step_length, phi_star, energy_star
def _zoom(self, alpha_lo, alpha_hi, phi_0, phiprime_0, def _zoom(self, alpha_lo, alpha_hi, phi_0, phiprime_0,
phi_lo, phiprime_lo, phi_hi, c1, c2): phi_lo, phiprime_lo, phi_hi):
"""Performs the second stage of the line search algorithm. """Performs the second stage of the line search algorithm.
If the first stage was not successful then the Zoom algorithm tries to If the first stage was not successful then the Zoom algorithm tries to
...@@ -203,7 +208,7 @@ class LineSearchStrongWolfe(LineSearch): ...@@ -203,7 +208,7 @@ class LineSearchStrongWolfe(LineSearch):
---------- ----------
alpha_lo : float alpha_lo : float
The lower boundary for the step length interval. The lower boundary for the step length interval.
alph_hi : float alpha_hi : float
The upper boundary for the step length interval. The upper boundary for the step length interval.
phi_0 : float phi_0 : float
Value of the energy at the starting point of the line search Value of the energy at the starting point of the line search
...@@ -219,10 +224,6 @@ class LineSearchStrongWolfe(LineSearch): ...@@ -219,10 +224,6 @@ class LineSearchStrongWolfe(LineSearch):
phi_hi : float phi_hi : float
Value of the energy if we perform a step of length alpha_hi in Value of the energy if we perform a step of length alpha_hi in
descent direction. descent direction.
c1 : float
Parameter for Armijo condition rule.
c2 : float
Parameter for curvature condition rule.
Returns Returns
------- -------
...@@ -234,6 +235,10 @@ class LineSearchStrongWolfe(LineSearch): ...@@ -234,6 +235,10 @@ class LineSearchStrongWolfe(LineSearch):
The new Energy object on the new position. The new Energy object on the new position.
""" """
print "entering zoom"
print alpha_lo, alpha_hi
print "pos1:",self.line_energy.at(alpha_lo).energy.position.val[0]
print "pos2:",self.line_energy.at(alpha_hi).energy.position.val[0]
max_iterations = self.max_zoom_iterations max_iterations = self.max_zoom_iterations
# define the cubic and quadratic interpolant checks # define the cubic and quadratic interpolant checks
cubic_delta = 0.2 # cubic cubic_delta = 0.2 # cubic
...@@ -262,28 +267,28 @@ class LineSearchStrongWolfe(LineSearch): ...@@ -262,28 +267,28 @@ class LineSearchStrongWolfe(LineSearch):
quad_check = quad_delta * delta_alpha quad_check = quad_delta * delta_alpha
alpha_j = self._quadmin(alpha_lo, phi_lo, phiprime_lo, alpha_j = self._quadmin(alpha_lo, phi_lo, phiprime_lo,
alpha_hi, phi_hi) alpha_hi, phi_hi)
# If quadratic was not successfull, try bisection # If quadratic was not successful, try bisection
if (alpha_j is None) or (alpha_j > b - quad_check) or \ if (alpha_j is None) or (alpha_j > b - quad_check) or \
(alpha_j < a + quad_check): (alpha_j < a + quad_check):
alpha_j = alpha_lo + 0.5*delta_alpha alpha_j = alpha_lo + 0.5*delta_alpha
# Check if the current value of alpha_j is already sufficient # Check if the current value of alpha_j is already sufficient
energy_alphaj = self.line_energy.at(alpha_j) le_alphaj = self.line_energy.at(alpha_j)
phi_alphaj = energy_alphaj.value phi_alphaj = le_alphaj.value
# If the first Wolfe condition is not met replace alpha_hi # If the first Wolfe condition is not met replace alpha_hi
# by alpha_j # by alpha_j
if (phi_alphaj > phi_0 + c1*alpha_j*phiprime_0) or\ if (phi_alphaj > phi_0 + self.c1*alpha_j*phiprime_0) or\
(phi_alphaj >= phi_lo): (phi_alphaj >= phi_lo):
alpha_recent, phi_recent = alpha_hi, phi_hi alpha_recent, phi_recent = alpha_hi, phi_hi
alpha_hi, phi_hi = alpha_j, phi_alphaj alpha_hi, phi_hi = alpha_j, phi_alphaj
else: else:
phiprime_alphaj = energy_alphaj.gradient phiprime_alphaj = le_alphaj.dd
# If the second Wolfe condition is met, return the result # If the second Wolfe condition is met, return the result
if abs(phiprime_alphaj) <= -c2*phiprime_0: if abs(phiprime_alphaj) <= -self.c2*phiprime_0:
alpha_star = alpha_j alpha_star = alpha_j
phi_star = phi_alphaj phi_star = phi_alphaj
energy_star = energy_alphaj le_star = le_alphaj
break break
# If not, check the sign of the slope # If not, check the sign of the slope
if phiprime_alphaj*delta_alpha >= 0: if phiprime_alphaj*delta_alpha >= 0:
...@@ -296,12 +301,12 @@ class LineSearchStrongWolfe(LineSearch): ...@@ -296,12 +301,12 @@ class LineSearchStrongWolfe(LineSearch):
phiprime_alphaj) phiprime_alphaj)
else: else:
alpha_star, phi_star, energy_star = \ alpha_star, phi_star, le_star = \
alpha_j, phi_alphaj, energy_alphaj alpha_j, phi_alphaj, le_alphaj
self.logger.error("The line search algorithm (zoom) did not " self.logger.error("The line search algorithm (zoom) did not "
"converge.") "converge.")
return alpha_star, phi_star, energy_star return alpha_star, phi_star, le_star
def _cubicmin(self, a, fa, fpa, b, fb, c, fc): def _cubicmin(self, a, fa, fpa, b, fb, c, fc):
"""Estimating the minimum with cubic interpolation. """Estimating the minimum with cubic interpolation.
......
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