There is a maintenance of MPCDF Gitlab on Thursday, April 22st 2020, 9:00 am CEST - Expect some service interruptions during this time

Commit 57cbbc96 authored by Martin Reinecke's avatar Martin Reinecke

Merge branch 'do_not_shadow_dir' into 'NIFTy_7'

Do not shadow native `dir`; rename to `direction`

See merge request !596
parents dffac406 de7697cc
Pipeline #93546 passed with stages
in 24 minutes and 27 seconds
......@@ -291,22 +291,22 @@ def _performance_check(op, pos, raise_on_fail):
def _get_acceptable_location(op, loc, lin):
if not np.isfinite(lin.val.s_sum()):
raise ValueError('Initial value must be finite')
dir = from_random(loc.domain, dtype=loc.dtype)
dirder = lin.jac(dir)
direction = from_random(loc.domain, dtype=loc.dtype)
dirder = lin.jac(direction)
if dirder.norm() == 0:
dir = dir * (lin.val.norm()*1e-5)
direction = direction * (lin.val.norm() * 1e-5)
else:
dir = dir * (lin.val.norm()*1e-5/dirder.norm())
direction = direction * (lin.val.norm() * 1e-5 / dirder.norm())
# Find a step length that leads to a "reasonable" location
for i in range(50):
try:
loc2 = loc+dir
loc2 = loc + direction
lin2 = op(Linearization.make_var(loc2, lin.want_metric))
if np.isfinite(lin2.val.s_sum()) and abs(lin2.val.s_sum()) < 1e20:
break
except FloatingPointError:
pass
dir = dir*0.5
direction = direction * 0.5
else:
raise ValueError("could not find a reasonable initial step")
return loc2, lin2
......@@ -368,21 +368,21 @@ def _jac_vs_finite_differences(op, loc, tol, ntries, only_r_differentiable):
for _ in range(ntries):
lin = op(Linearization.make_var(loc))
loc2, lin2 = _get_acceptable_location(op, loc, lin)
dir = loc2-loc
direction = loc2 - loc
locnext = loc2
dirnorm = dir.norm()
dirnorm = direction.norm()
hist = []
for i in range(50):
locmid = loc + 0.5*dir
locmid = loc + 0.5 * direction
linmid = op(Linearization.make_var(locmid))
dirder = linmid.jac(dir)
numgrad = (lin2.val-lin.val)
dirder = linmid.jac(direction)
numgrad = (lin2.val - lin.val)
xtol = tol * dirder.norm() / np.sqrt(dirder.size)
hist.append((numgrad-dirder).norm())
# print(len(hist),hist[-1])
if (abs(numgrad-dirder) <= xtol).s_all():
hist.append((numgrad - dirder).norm())
# print(len(hist),hist[-1])
if (abs(numgrad - dirder) <= xtol).s_all():
break
dir = dir*0.5
direction = direction * 0.5
dirnorm *= 0.5
loc2, lin2 = locmid, linmid
else:
......
......@@ -42,11 +42,11 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, sig, erf):
out = [None]*nlos
for i in range(nlos):
dir = end[:, i]-start[:, i]
dirx = np.where(dir == 0., 1e-12, dir)
d0 = np.where(dir == 0., ((start[:, i] > 0)-0.5)*1e12,
direction = end[:, i]-start[:, i]
dirx = np.where(direction == 0., 1e-12, direction)
d0 = np.where(direction == 0., ((start[:, i] > 0)-0.5)*1e12,
-start[:, i]/dirx)
d1 = np.where(dir == 0., ((start[:, i] < pmax)-0.5)*-1e12,
d1 = np.where(direction == 0., ((start[:, i] < pmax)-0.5)*-1e12,
(pmax-start[:, i])/dirx)
(dmin, dmax) = (np.minimum(d0, d1), np.maximum(d0, d1))
dmin = dmin.max()
......@@ -61,18 +61,18 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, sig, erf):
out[i] = (np.full(0, 0, dtype=np.int64), np.full(0, 0.))
continue
# determine coordinates of first cell crossing
c_first = np.ceil(start[:, i]+dir*dmin)
c_first = np.where(dir > 0., c_first, c_first-1.)
c_first = np.ceil(start[:, i]+direction*dmin)
c_first = np.where(direction > 0., c_first, c_first-1.)
c_first = (c_first-start[:, i])/dirx
pos1 = np.asarray((start[:, i]+dmin*dir), dtype=np.int)
pos1 = np.asarray((start[:, i]+dmin*direction), dtype=np.int)
pos1 = np.sum(pos1*inc)
cdist = np.empty(0, dtype=np.float64)
add = np.empty(0, dtype=np.int)
for j in range(ndim):
if dir[j] != 0:
step = inc[j] if dir[j] > 0 else -inc[j]
if direction[j] != 0:
step = inc[j] if direction[j] > 0 else -inc[j]
tmp = np.arange(start=c_first[j], stop=dmax,
step=abs(1./dir[j]))
step=abs(1./direction[j]))
cdist = np.append(cdist, tmp)
add = np.append(add, np.full(len(tmp), step, dtype=np.int64))
idx = np.argsort(cdist)
......@@ -80,7 +80,7 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, sig, erf):
add = add[idx]
cdist = np.append(np.full(1, dmin), cdist)
cdist = np.append(cdist, np.full(1, dmax))
corfac = np.linalg.norm(dir*dist)
corfac = np.linalg.norm(direction*dist)
cdist *= corfac
wgt = np.diff(cdist)
mdist = 0.5*(cdist[:-1]+cdist[1:])
......
......@@ -119,18 +119,18 @@ class Energy(metaclass=NiftyMeta):
"""
raise NotImplementedError
def longest_step(self, dir):
"""Returns the longest allowed step size along `dir`
def longest_step(self, direction):
"""Returns the longest allowed step size along `direction`
Parameters
----------
dir : Field
direction : Field
the search direction
Returns
-------
float or None
the longest allowed step when starting from `self.position` along
`dir`. If None, the step size is not limited.
`direction`. If None, the step size is not limited.
"""
return None
......@@ -242,7 +242,7 @@ class GaussianEnergy(EnergyOperator):
not. Note that for a complex Gaussian the inverse_covariance is
.. math ::
(<ff^dagger>)^{-1}_P(f)/2,
where the additional factor of 2 is necessary because the
where the additional factor of 2 is necessary because the
domain of s has double as many dimensions as in the real case.
Note
......
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