Commit 120462aa authored by Martin Reinecke's avatar Martin Reinecke

working, but too slowly

parent 3934ca0c
...@@ -21,7 +21,7 @@ if __name__ == "__main__": ...@@ -21,7 +21,7 @@ if __name__ == "__main__":
ends = np.random.random((2, n)) ends = np.random.random((2, n))
return starts, ends return starts, ends
nlos = 100 nlos = 1000
starts, ends = make_random_los(nlos) starts, ends = make_random_los(nlos)
R = ift.library.LOSResponse(s_space, starts=starts, ends=ends) R = ift.library.LOSResponse(s_space, starts=starts, ends=ends)
......
...@@ -10,37 +10,35 @@ from .. import dobj ...@@ -10,37 +10,35 @@ from .. import dobj
class Pixellister(object): class Pixellister(object):
Rmax = np.sqrt(3.)*0.5
def __init__(self, start, end): def __init__(self, start, end):
self._start = start self._start = start
self._end = end self._end = end
self._dir = end - start self._dir = end - start
self._dirlen = np.sqrt(np.sum(self._dir*self._dir)) self._dirlen = np.linalg.norm(self._dir)
self._dir *= 1./self._dirlen self._dir *= 1./self._dirlen
def _losdist(self, pos): def _losdist(self, pos):
t1 = pos - self._start t1 = pos - self._start
dotpr = np.sum(self._dir*t1) dotpr = np.dot(self._dir, t1)
if dotpr < 0: if dotpr < 0:
return np.sqrt(np.sum(t1*t1)) return np.linalg.norm(t1)
elif dotpr <= self._dirlen: elif dotpr <= self._dirlen:
t2 = self._start + self._dir*dotpr return np.linalg.norm(pos - self._start - self._dir*dotpr)
d1 = pos-t2
return np.sqrt(np.sum(d1*d1))
else: else:
d3 = pos-self._end return np.linalg.norm(pos-self._end)
return np.sqrt(np.sum(d3*d3))
def build_pixellist(self, lo, hi, pixels): def build_pixellist(self, lo, hi, pixels):
Rmax=5.
diff = hi-lo diff = hi-lo
if np.all(diff==1): if np.all(diff == 1):
dist = self._losdist(lo) dist = self._losdist(lo)
if dist < Rmax: if dist < self.Rmax:
pixels.append((lo, dist)) pixels.append((lo, dist))
return return
mid = (lo+hi-1.)/2. mid = (lo+hi-1.)*0.5
t1 = mid-lo t1 = mid-lo
if self._losdist(mid) < np.sqrt(np.sum(t1*t1)) + Rmax: if self._losdist(mid) < np.linalg.norm(t1) + self.Rmax:
dim = np.argmax(diff) dim = np.argmax(diff)
imid = lo[dim]+diff[dim]//2 imid = lo[dim]+diff[dim]//2
lmid = lo.copy() lmid = lo.copy()
...@@ -50,6 +48,7 @@ class Pixellister(object): ...@@ -50,6 +48,7 @@ class Pixellister(object):
self.build_pixellist(lo, hmid, pixels) self.build_pixellist(lo, hmid, pixels)
self.build_pixellist(lmid, hi, pixels) self.build_pixellist(lmid, hi, pixels)
def _comp_traverse(start, end, shp): def _comp_traverse(start, end, shp):
nlos = start.shape[1] nlos = start.shape[1]
inc = np.full(len(shp), 1) inc = np.full(len(shp), 1)
...@@ -59,10 +58,10 @@ def _comp_traverse(start, end, shp): ...@@ -59,10 +58,10 @@ def _comp_traverse(start, end, shp):
for i in range(nlos): for i in range(nlos):
print i print i
pixels = [] pixels = []
pl = Pixellister(start[:,i], end[:,i]) pl = Pixellister(start[:, i], end[:, i])
pl.build_pixellist(np.array(shp)*0, np.array(shp), pixels) pl.build_pixellist(np.array(shp)*0, np.array(shp), pixels)
pix1 = [np.sum(inc*p[0]) for p in pixels] pix1 = np.array([np.sum(inc*p[0]) for p in pixels])
wgt = [np.exp(-p[1]*p[1]) for p in pixels] wgt = np.array([np.exp(-p[1]*p[1]) for p in pixels])
out[i] = (pix1, wgt) out[i] = (pix1, wgt)
return out return out
...@@ -141,19 +140,20 @@ class LOSResponse(LinearOperator): ...@@ -141,19 +140,20 @@ class LOSResponse(LinearOperator):
cnt = 0 cnt = 0
for i in w_i: for i in w_i:
nval = len(i[1]) nval = len(i[1])
ilos[ofs:ofs+nval] = cnt if nval > 0:
iarr[ofs:ofs+nval] = i[0] ilos[ofs:ofs+nval] = cnt
xwgt[ofs:ofs+nval] = i[1] iarr[ofs:ofs+nval] = i[0]
fullidx = np.unravel_index(i[0], self._local_shape) xwgt[ofs:ofs+nval] = i[1]
tmp = np.zeros(nval, dtype=np.float64) fullidx = np.unravel_index(i[0], self._local_shape)
fct = 1. tmp = np.zeros(nval, dtype=np.float64)
for j in range(ndim): fct = 1.
tmp += (fullidx[j]//boxsz)*fct for j in range(ndim):
fct *= self._local_shape[j] tmp += (fullidx[j]//boxsz)*fct
tmp += cnt/float(nlos) fct *= self._local_shape[j]
tmp += iarr[ofs:ofs+nval]/float(nlos*npix) tmp += cnt/float(nlos)
pri[ofs:ofs+nval] = tmp tmp += iarr[ofs:ofs+nval]/float(nlos*npix)
ofs += nval pri[ofs:ofs+nval] = tmp
ofs += nval
cnt += 1 cnt += 1
xtmp = np.argsort(pri) xtmp = np.argsort(pri)
ilos = ilos[xtmp] ilos = ilos[xtmp]
......
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