Commit cd7dfa8c authored by Martin Reinecke's avatar Martin Reinecke

explicitly disallow in-place operations on Fields and MultiFields

parent 0be9da82
...@@ -141,9 +141,7 @@ class RGSpace(StructuredDomain): ...@@ -141,9 +141,7 @@ class RGSpace(StructuredDomain):
@staticmethod @staticmethod
def _kernel(x, sigma): def _kernel(x, sigma):
from ..sugar import exp from ..sugar import exp
tmp = x*x return exp(x*x * (-2.*np.pi*np.pi*sigma*sigma))
tmp *= -2.*np.pi*np.pi*sigma*sigma
return exp(tmp)
def get_fft_smoothing_kernel_function(self, sigma): def get_fft_smoothing_kernel_function(self, sigma):
if (not self.harmonic): if (not self.harmonic):
......
...@@ -31,7 +31,7 @@ def _get_acceptable_model(M): ...@@ -31,7 +31,7 @@ def _get_acceptable_model(M):
raise ValueError('Initial Model value must be finite') raise ValueError('Initial Model value must be finite')
dir = from_random("normal", M.position.domain) dir = from_random("normal", M.position.domain)
dirder = M.jacobian(dir) dirder = M.jacobian(dir)
dir *= val/(dirder).norm()*1e-5 dir = dir * val * (1e-5/dirder.norm())
# Find a step length that leads to a "reasonable" Model # Find a step length that leads to a "reasonable" Model
for i in range(50): for i in range(50):
try: try:
...@@ -40,7 +40,7 @@ def _get_acceptable_model(M): ...@@ -40,7 +40,7 @@ def _get_acceptable_model(M):
break break
except FloatingPointError: except FloatingPointError:
pass pass
dir *= 0.5 dir = dir*0.5
else: else:
raise ValueError("could not find a reasonable initial step") raise ValueError("could not find a reasonable initial step")
return M2 return M2
...@@ -52,7 +52,7 @@ def _get_acceptable_energy(E): ...@@ -52,7 +52,7 @@ def _get_acceptable_energy(E):
raise ValueError('Initial Energy must be finite') raise ValueError('Initial Energy must be finite')
dir = from_random("normal", E.position.domain) dir = from_random("normal", E.position.domain)
dirder = E.gradient.vdot(dir) dirder = E.gradient.vdot(dir)
dir *= np.abs(val)/np.abs(dirder)*1e-5 dir = dir * (np.abs(val)/np.abs(dirder)*1e-5)
# Find a step length that leads to a "reasonable" energy # Find a step length that leads to a "reasonable" energy
for i in range(50): for i in range(50):
try: try:
...@@ -61,7 +61,7 @@ def _get_acceptable_energy(E): ...@@ -61,7 +61,7 @@ def _get_acceptable_energy(E):
break break
except FloatingPointError: except FloatingPointError:
pass pass
dir *= 0.5 dir = dir*0.5
else: else:
raise ValueError("could not find a reasonable initial step") raise ValueError("could not find a reasonable initial step")
return E2 return E2
...@@ -92,7 +92,7 @@ def check_value_gradient_consistency(E, tol=1e-8, ntries=100): ...@@ -92,7 +92,7 @@ def check_value_gradient_consistency(E, tol=1e-8, ntries=100):
xtol = tol*Emid.gradient_norm xtol = tol*Emid.gradient_norm
if abs(numgrad-dirder) < xtol: if abs(numgrad-dirder) < xtol:
break break
dir *= 0.5 dir = dir*0.5
dirnorm *= 0.5 dirnorm *= 0.5
E2 = Emid E2 = Emid
else: else:
...@@ -117,7 +117,7 @@ def check_value_gradient_metric_consistency(E, tol=1e-8, ntries=100): ...@@ -117,7 +117,7 @@ def check_value_gradient_metric_consistency(E, tol=1e-8, ntries=100):
if abs((E2.value-val)/dirnorm - dirder) < xtol and \ if abs((E2.value-val)/dirnorm - dirder) < xtol and \
(abs((E2.gradient-E.gradient)/dirnorm-dgrad) < xtol).all(): (abs((E2.gradient-E.gradient)/dirnorm-dgrad) < xtol).all():
break break
dir *= 0.5 dir = dir*0.5
dirnorm *= 0.5 dirnorm *= 0.5
E2 = Emid E2 = Emid
else: else:
......
...@@ -670,3 +670,12 @@ for op in ["__add__", "__radd__", ...@@ -670,3 +670,12 @@ for op in ["__add__", "__radd__",
return NotImplemented return NotImplemented
return func2 return func2
setattr(Field, op, func(op)) setattr(Field, op, func(op))
for op in ["__iadd__", "__isub__", "__imul__", "__idiv__",
"__itruediv__", "__ifloordiv__", "__ipow__"]:
def func(op):
def func2(self, other):
raise TypeError(
"In-place operations are deliberately not supported")
return func2
setattr(Field, op, func(op))
...@@ -84,8 +84,7 @@ class ConjugateGradient(Minimizer): ...@@ -84,8 +84,7 @@ class ConjugateGradient(Minimizer):
logger.error("Error: ConjugateGradient: alpha<0.") logger.error("Error: ConjugateGradient: alpha<0.")
return energy, controller.ERROR return energy, controller.ERROR
q *= -alpha r = r - q*alpha
r = r + q
energy = energy.at_with_grad(energy.position - alpha*d, r) energy = energy.at_with_grad(energy.position - alpha*d, r)
...@@ -103,7 +102,6 @@ class ConjugateGradient(Minimizer): ...@@ -103,7 +102,6 @@ class ConjugateGradient(Minimizer):
if status != controller.CONTINUE: if status != controller.CONTINUE:
return energy, status return energy, status
d *= max(0, gamma/previous_gamma) d = d * max(0, gamma/previous_gamma) + s
d += s
previous_gamma = gamma previous_gamma = gamma
...@@ -61,16 +61,16 @@ class L_BFGS(DescentMinimizer): ...@@ -61,16 +61,16 @@ class L_BFGS(DescentMinimizer):
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].vdot(p)/s[idx].vdot(y[idx])
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].vdot(y[idx]) / y[idx].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 *= 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].vdot(p) / s[idx].vdot(y[idx])
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
self._k += 1 self._k += 1
......
...@@ -67,7 +67,7 @@ class VL_BFGS(DescentMinimizer): ...@@ -67,7 +67,7 @@ class VL_BFGS(DescentMinimizer):
descent_direction = delta[0] * b[0] descent_direction = delta[0] * b[0]
for i in range(1, len(delta)): for i in range(1, len(delta)):
descent_direction += delta[i] * b[i] descent_direction = descent_direction + delta[i]*b[i]
return descent_direction return descent_direction
......
...@@ -220,3 +220,12 @@ for op in ["__add__", "__radd__", ...@@ -220,3 +220,12 @@ for op in ["__add__", "__radd__",
return MultiField(result_val) return MultiField(result_val)
return func2 return func2
setattr(MultiField, op, func(op)) setattr(MultiField, op, func(op))
for op in ["__iadd__", "__isub__", "__imul__", "__idiv__",
"__itruediv__", "__ifloordiv__", "__ipow__"]:
def func(op):
def func2(self, other):
raise TypeError(
"In-place operations are deliberately not supported")
return func2
setattr(MultiField, op, func(op))
...@@ -185,7 +185,6 @@ class DiagonalOperator(EndomorphicOperator): ...@@ -185,7 +185,6 @@ class DiagonalOperator(EndomorphicOperator):
res = Field.from_random(random_type="normal", domain=self._domain, res = Field.from_random(random_type="normal", domain=self._domain,
dtype=dtype) dtype=dtype)
if from_inverse: if from_inverse:
res /= np.sqrt(self._ldiag) return res/np.sqrt(self._ldiag)
else: else:
res *= np.sqrt(self._ldiag) return res*np.sqrt(self._ldiag)
return res
...@@ -118,10 +118,7 @@ class FFTOperator(LinearOperator): ...@@ -118,10 +118,7 @@ class FFTOperator(LinearOperator):
fct = self._domain[self._space].scalar_dvol fct = self._domain[self._space].scalar_dvol
else: else:
fct = self._target[self._space].scalar_dvol fct = self._target[self._space].scalar_dvol
if fct != 1: return Tval if fct == 1 else Tval*fct
Tval *= fct
return Tval
@property @property
def domain(self): def domain(self):
......
...@@ -172,14 +172,15 @@ class SumOperator(LinearOperator): ...@@ -172,14 +172,15 @@ class SumOperator(LinearOperator):
def apply(self, x, mode): def apply(self, x, mode):
self._check_mode(mode) self._check_mode(mode)
for i, op in enumerate(self._ops): res = None
if i == 0: for op, neg in zip(self._ops, self._neg):
res = -op.apply(x, mode) if self._neg[i] else op.apply(x, mode) if res is None:
res = -op.apply(x, mode) if neg else op.apply(x, mode)
else: else:
if self._neg[i]: if neg:
res -= op.apply(x, mode) res = res - op.apply(x, mode)
else: else:
res += op.apply(x, mode) res = res + op.apply(x, mode)
return res return res
def draw_sample(self, from_inverse=False, dtype=np.float64): def draw_sample(self, from_inverse=False, dtype=np.float64):
......
...@@ -31,9 +31,9 @@ class StatCalculator(object): ...@@ -31,9 +31,9 @@ class StatCalculator(object):
self._M2 = 0.*value self._M2 = 0.*value
else: else:
delta = value - self._mean delta = value - self._mean
self._mean += delta*(1./self._count) self._mean = self.mean + delta*(1./self._count)
delta2 = value - self._mean delta2 = value - self._mean
self._M2 += delta*delta2 self._M2 = self._M2 + delta*delta2
@property @property
def mean(self): def mean(self):
......
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