Commit f5205c8f authored by Martin Reinecke's avatar Martin Reinecke Committed by Philipp Arras
Browse files

more fixes

parent 9523951a
...@@ -70,14 +70,14 @@ class ConjugateGradient(Minimizer): ...@@ -70,14 +70,14 @@ class ConjugateGradient(Minimizer):
r = energy.gradient r = energy.gradient
d = r if preconditioner is None else preconditioner(r) d = r if preconditioner is None else preconditioner(r)
previous_gamma = r.vdot(d).real previous_gamma = r.s_vdot(d).real
if previous_gamma == 0: if previous_gamma == 0:
return energy, controller.CONVERGED return energy, controller.CONVERGED
ii = 0 ii = 0
while True: while True:
q = energy.apply_metric(d) q = energy.apply_metric(d)
curv = d.vdot(q).real curv = d.s_vdot(q).real
if curv == 0.: if curv == 0.:
logger.error("Error: ConjugateGradient: curv==0.") logger.error("Error: ConjugateGradient: curv==0.")
return energy, controller.ERROR return energy, controller.ERROR
...@@ -98,7 +98,7 @@ class ConjugateGradient(Minimizer): ...@@ -98,7 +98,7 @@ class ConjugateGradient(Minimizer):
s = r if preconditioner is None else preconditioner(r) s = r if preconditioner is None else preconditioner(r)
gamma = r.vdot(s).real gamma = r.s_vdot(s).real
if gamma < 0: if gamma < 0:
logger.error( logger.error(
"Positive definiteness of preconditioner violated!") "Positive definiteness of preconditioner violated!")
......
...@@ -227,16 +227,16 @@ class L_BFGS(DescentMinimizer): ...@@ -227,16 +227,16 @@ class L_BFGS(DescentMinimizer):
if nhist > 0: if nhist > 0:
for i in range(k-1, k-nhist-1, -1): for i in range(k-1, k-nhist-1, -1):
idx = i % maxhist idx = i % maxhist
alpha[idx] = s[idx].vdot(p)/s[idx].vdot(y[idx]) alpha[idx] = s[idx].s_vdot(p)/s[idx].s_vdot(y[idx])
p = p - alpha[idx]*y[idx] p = p - alpha[idx]*y[idx]
idx = (k-1) % maxhist idx = (k-1) % maxhist
fact = s[idx].vdot(y[idx]) / y[idx].vdot(y[idx]) fact = s[idx].s_vdot(y[idx]) / y[idx].s_vdot(y[idx])
if fact <= 0.: if fact <= 0.:
logger.error("L-BFGS curvature not positive definite!") logger.error("L-BFGS curvature not positive definite!")
p = p*fact p = p*fact
for i in range(k-nhist, k): for i in range(k-nhist, k):
idx = i % maxhist idx = i % maxhist
beta = y[idx].vdot(p) / s[idx].vdot(y[idx]) beta = y[idx].s_vdot(p) / s[idx].s_vdot(y[idx])
p = p + (alpha[idx]-beta)*s[idx] p = p + (alpha[idx]-beta)*s[idx]
self._lastx = x self._lastx = x
self._lastgrad = gradient self._lastgrad = gradient
...@@ -388,12 +388,12 @@ class _InformationStore(object): ...@@ -388,12 +388,12 @@ class _InformationStore(object):
k1 = (k-1) % mmax k1 = (k-1) % mmax
for i in range(m): for i in range(m):
kmi = (k-m+i) % mmax kmi = (k-m+i) % mmax
self.ss[kmi, k1] = self.ss[k1, kmi] = self.s[kmi].vdot(self.s[k1]) self.ss[kmi, k1] = self.ss[k1, kmi] = self.s[kmi].s_vdot(self.s[k1])
self.yy[kmi, k1] = self.yy[k1, kmi] = self.y[kmi].vdot(self.y[k1]) self.yy[kmi, k1] = self.yy[k1, kmi] = self.y[kmi].s_vdot(self.y[k1])
self.sy[kmi, k1] = self.s[kmi].vdot(self.y[k1]) self.sy[kmi, k1] = self.s[kmi].s_vdot(self.y[k1])
for j in range(m-1): for j in range(m-1):
kmj = (k-m+j) % mmax kmj = (k-m+j) % mmax
self.sy[k1, kmj] = self.s[k1].vdot(self.y[kmj]) self.sy[k1, kmj] = self.s[k1].s_vdot(self.y[kmj])
for i in range(m): for i in range(m):
kmi = (k-m+i) % mmax kmi = (k-m+i) % mmax
...@@ -403,10 +403,10 @@ class _InformationStore(object): ...@@ -403,10 +403,10 @@ class _InformationStore(object):
result[i, m+j] = result[m+j, i] = self.sy[kmi, kmj] result[i, m+j] = result[m+j, i] = self.sy[kmi, kmj]
result[m+i, m+j] = self.yy[kmi, kmj] result[m+i, m+j] = self.yy[kmi, kmj]
sgrad_i = self.s[kmi].vdot(self.last_gradient) sgrad_i = self.s[kmi].s_vdot(self.last_gradient)
result[2*m, i] = result[i, 2*m] = sgrad_i result[2*m, i] = result[i, 2*m] = sgrad_i
ygrad_i = self.y[kmi].vdot(self.last_gradient) ygrad_i = self.y[kmi].s_vdot(self.last_gradient)
result[2*m, m+i] = result[m+i, 2*m] = ygrad_i result[2*m, m+i] = result[m+i, 2*m] = ygrad_i
result[2*m, 2*m] = self.last_gradient.norm() result[2*m, 2*m] = self.last_gradient.norm()
......
...@@ -95,7 +95,7 @@ class LineEnergy(object): ...@@ -95,7 +95,7 @@ class LineEnergy(object):
""" """
float : The directional derivative at the given `position`. float : The directional derivative at the given `position`.
""" """
res = self._energy.gradient.vdot(self._line_direction) res = self._energy.gradient.s_vdot(self._line_direction)
if abs(res.imag) / max(abs(res.real), 1.) > 1e-12: if abs(res.imag) / max(abs(res.real), 1.) > 1e-12:
from ..logger import logger from ..logger import logger
logger.warning("directional derivative has non-negligible " logger.warning("directional derivative has non-negligible "
......
...@@ -75,17 +75,17 @@ class NonlinearCG(Minimizer): ...@@ -75,17 +75,17 @@ class NonlinearCG(Minimizer):
if self._beta_heuristic == 'Hestenes-Stiefel': if self._beta_heuristic == 'Hestenes-Stiefel':
# Eq. (5.46) in Nocedal & Wright. # Eq. (5.46) in Nocedal & Wright.
beta = max(0.0, (grad_new.vdot(grad_new-grad_old) / beta = max(0.0, (grad_new.s_vdot(grad_new-grad_old) /
(grad_new-grad_old).vdot(p)).real) (grad_new-grad_old).s_vdot(p)).real)
elif self._beta_heuristic == 'Polak-Ribiere': elif self._beta_heuristic == 'Polak-Ribiere':
# Eq. (5.44) in Nocedal & Wright. (with (5.45) additionally) # Eq. (5.44) in Nocedal & Wright. (with (5.45) additionally)
beta = max(0.0, (grad_new.vdot(grad_new-grad_old) / beta = max(0.0, (grad_new.s_vdot(grad_new-grad_old) /
(grad_old.vdot(grad_old))).real) (grad_old.s_vdot(grad_old))).real)
elif self._beta_heuristic == 'Fletcher-Reeves': elif self._beta_heuristic == 'Fletcher-Reeves':
# Eq. (5.41a) in Nocedal & Wright. # Eq. (5.41a) in Nocedal & Wright.
beta = (grad_new.vdot(grad_new)/(grad_old.vdot(grad_old))).real beta = (grad_new.s_vdot(grad_new)/(grad_old.s_vdot(grad_old))).real
else: else:
# Eq. (5.49) in Nocedal & Wright. # Eq. (5.49) in Nocedal & Wright.
beta = (grad_new.vdot(grad_new) / beta = (grad_new.s_vdot(grad_new) /
((grad_new-grad_old).vdot(p))).real ((grad_new-grad_old).s_vdot(p))).real
p = beta*p - grad_new p = beta*p - grad_new
...@@ -34,9 +34,9 @@ class QuadraticEnergy(Energy): ...@@ -34,9 +34,9 @@ class QuadraticEnergy(Energy):
else: else:
Ax = self._A(self._position) Ax = self._A(self._position)
self._grad = Ax if b is None else Ax - b self._grad = Ax if b is None else Ax - b
self._value = 0.5*self._position.vdot(Ax) self._value = 0.5*self._position.s_vdot(Ax)
if b is not None: if b is not None:
self._value -= b.vdot(self._position) self._value -= b.s_vdot(self._position)
def at(self, position): def at(self, position):
return QuadraticEnergy(position, self._A, self._b) return QuadraticEnergy(position, self._A, self._b)
......
...@@ -74,7 +74,7 @@ class _MinHelper(object): ...@@ -74,7 +74,7 @@ class _MinHelper(object):
def _update(self, x): def _update(self, x):
pos = _toField(x, self._energy.position) pos = _toField(x, self._energy.position)
if (pos != self._energy.position).any(): if (pos != self._energy.position).s_any():
self._energy = self._energy.at(pos) self._energy = self._energy.at(pos)
def fun(self, x): def fun(self, x):
......
...@@ -169,7 +169,7 @@ class MultiField(object): ...@@ -169,7 +169,7 @@ class MultiField(object):
return (nrm ** ord).sum() ** (1./ord) return (nrm ** ord).sum() ** (1./ord)
# return np.sqrt(np.abs(self.vdot(x=self))) # return np.sqrt(np.abs(self.vdot(x=self)))
def sum(self): def s_sum(self):
"""Computes the sum all field values. """Computes the sum all field values.
Returns Returns
...@@ -177,7 +177,7 @@ class MultiField(object): ...@@ -177,7 +177,7 @@ class MultiField(object):
norm : float norm : float
The sum of the field values. The sum of the field values.
""" """
return utilities.my_sum(map(lambda v: v.sum(), self._val)) return utilities.my_sum(map(lambda v: v.s_sum(), self._val))
@property @property
def size(self): def size(self):
...@@ -207,15 +207,15 @@ class MultiField(object): ...@@ -207,15 +207,15 @@ class MultiField(object):
self._domain, self._domain,
tuple(self._val[i].clip(lmin[i], lmax[i]) for i in range(ncomp))) tuple(self._val[i].clip(lmin[i], lmax[i]) for i in range(ncomp)))
def all(self): def s_all(self):
for v in self._val: for v in self._val:
if not v.all(): if not v.s_all():
return False return False
return True return True
def any(self): def s_any(self):
for v in self._val: for v in self._val:
if v.any(): if v.s_any():
return True return True
return False return False
......
...@@ -26,7 +26,7 @@ dom = ift.makeDomain({"d1": ift.RGSpace(10)}) ...@@ -26,7 +26,7 @@ dom = ift.makeDomain({"d1": ift.RGSpace(10)})
def test_vdot(): def test_vdot():
f1 = ift.from_random("normal", domain=dom, dtype=np.complex128) f1 = ift.from_random("normal", domain=dom, dtype=np.complex128)
f2 = ift.from_random("normal", domain=dom, dtype=np.complex128) f2 = ift.from_random("normal", domain=dom, dtype=np.complex128)
assert_allclose(f1.vdot(f2), np.conj(f2.vdot(f1))) assert_allclose(f1.s_vdot(f2), np.conj(f2.s_vdot(f1)))
def test_func(): def test_func():
...@@ -38,7 +38,7 @@ def test_func(): ...@@ -38,7 +38,7 @@ def test_func():
def test_multifield_field_consistency(): def test_multifield_field_consistency():
f1 = ift.full(dom, 27) f1 = ift.full(dom, 27)
f2 = ift.makeField(dom['d1'], f1['d1'].val) f2 = ift.makeField(dom['d1'], f1['d1'].val)
assert_equal(f1.sum(), f2.sum()) assert_equal(f1.s_sum(), f2.s_sum())
assert_equal(f1.size, f2.size) assert_equal(f1.size, f2.size)
...@@ -64,9 +64,9 @@ def test_blockdiagonal(): ...@@ -64,9 +64,9 @@ def test_blockdiagonal():
assert_equal(type(op2), ift.BlockDiagonalOperator) assert_equal(type(op2), ift.BlockDiagonalOperator)
f1 = op2(ift.full(dom, 1)) f1 = op2(ift.full(dom, 1))
for val in f1.values(): for val in f1.values():
assert_equal((val == 400).all(), True) assert_equal((val == 400).s_all(), True)
op2 = op + op op2 = op + op
assert_equal(type(op2), ift.BlockDiagonalOperator) assert_equal(type(op2), ift.BlockDiagonalOperator)
f1 = op2(ift.full(dom, 1)) f1 = op2(ift.full(dom, 1))
for val in f1.values(): for val in f1.values():
assert_equal((val == 40).all(), True) assert_equal((val == 40).s_all(), True)
...@@ -43,8 +43,8 @@ def test_times_adjoint_times(space1, space2): ...@@ -43,8 +43,8 @@ def test_times_adjoint_times(space1, space2):
rand1 = ift.Field.from_random('normal', domain=(space1, space2)) rand1 = ift.Field.from_random('normal', domain=(space1, space2))
rand2 = ift.Field.from_random('normal', domain=(space1, space2)) rand2 = ift.Field.from_random('normal', domain=(space1, space2))
tt1 = rand2.vdot(op.times(rand1)) tt1 = rand2.s_vdot(op.times(rand1))
tt2 = rand1.vdot(op.adjoint_times(rand2)) tt2 = rand1.s_vdot(op.adjoint_times(rand2))
assert_allclose(tt1, tt2) assert_allclose(tt1, tt2)
......
...@@ -42,8 +42,8 @@ def test_times_adjoint(space): ...@@ -42,8 +42,8 @@ def test_times_adjoint(space):
rand2 = ift.Field.from_random('normal', domain=space) rand2 = ift.Field.from_random('normal', domain=space)
diag = ift.Field.from_random('normal', domain=space) diag = ift.Field.from_random('normal', domain=space)
D = ift.DiagonalOperator(diag) D = ift.DiagonalOperator(diag)
tt1 = rand1.vdot(D.times(rand2)) tt1 = rand1.s_vdot(D.times(rand2))
tt2 = rand2.vdot(D.times(rand1)) tt2 = rand2.s_vdot(D.times(rand1))
assert_allclose(tt1, tt2) assert_allclose(tt1, tt2)
......
...@@ -47,8 +47,8 @@ def test_adjoint_times(space, sigma): ...@@ -47,8 +47,8 @@ def test_adjoint_times(space, sigma):
op = ift.HarmonicSmoothingOperator(space, sigma=sigma) op = ift.HarmonicSmoothingOperator(space, sigma=sigma)
rand1 = ift.Field.from_random('normal', domain=space) rand1 = ift.Field.from_random('normal', domain=space)
rand2 = ift.Field.from_random('normal', domain=space) rand2 = ift.Field.from_random('normal', domain=space)
tt1 = rand1.vdot(op.times(rand2)) tt1 = rand1.s_vdot(op.times(rand2))
tt2 = rand2.vdot(op.adjoint_times(rand1)) tt2 = rand2.s_vdot(op.adjoint_times(rand1))
assert_allclose(tt1, tt2) assert_allclose(tt1, tt2)
...@@ -58,7 +58,7 @@ def test_times(space, sigma): ...@@ -58,7 +58,7 @@ def test_times(space, sigma):
fld[0] = 1. fld[0] = 1.
rand1 = ift.Field.from_raw(space, fld) rand1 = ift.Field.from_raw(space, fld)
tt1 = op.times(rand1) tt1 = op.times(rand1)
assert_allclose(1, tt1.sum()) assert_allclose(1, tt1.s_sum())
@pmp('sz', [128, 256]) @pmp('sz', [128, 256])
...@@ -70,7 +70,7 @@ def test_smooth_regular1(sz, d, sigma, tp): ...@@ -70,7 +70,7 @@ def test_smooth_regular1(sz, d, sigma, tp):
inp = ift.Field.from_random( inp = ift.Field.from_random(
domain=sp, random_type='normal', std=1, mean=4, dtype=tp) domain=sp, random_type='normal', std=1, mean=4, dtype=tp)
out = smo(inp) out = smo(inp)
assert_allclose(inp.sum(), out.sum(), rtol=tol, atol=tol) assert_allclose(inp.s_sum(), out.s_sum(), rtol=tol, atol=tol)
@pmp('sz1', [10, 15]) @pmp('sz1', [10, 15])
...@@ -84,4 +84,4 @@ def test_smooth_regular2(sz1, sz2, d1, d2, sigma, tp): ...@@ -84,4 +84,4 @@ def test_smooth_regular2(sz1, sz2, d1, d2, sigma, tp):
inp = ift.Field.from_random( inp = ift.Field.from_random(
domain=sp, random_type='normal', std=1, mean=4, dtype=tp) domain=sp, random_type='normal', std=1, mean=4, dtype=tp)
out = smo(inp) out = smo(inp)
assert_allclose(inp.sum(), out.sum(), rtol=tol, atol=tol) assert_allclose(inp.s_sum(), out.s_sum(), rtol=tol, atol=tol)
Supports Markdown
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