293 def solve(self, science, fringes):
294 """Solve for the scale factors with iterative clipping.
295
296 Parameters
297 ----------
298 science : `numpy.array`
299 Array of measured science image values at each of the
300 positions supplied.
301 fringes : `numpy.array`
302 Array of measured fringe values at each of the positions
303 supplied.
304
305 Returns
306 -------
307 solution : `np.array`
308 Fringe solution amplitudes for each input fringe frame.
309 rms : `float`
310 RMS error for the fit solution for this exposure.
311 """
312 import lsstDebug
314
315 origNum = len(science)
316
317 def emptyResult(msg=""):
318 """Generate an empty result for return to the user
319
320 There are no good pixels; doesn't matter what we return.
321 """
322 self.log.warning("Unable to solve for fringes: no good pixels%s", msg)
323 out = [0]
324 if len(fringes) > 1:
325 out = out*len(fringes)
326 return numpy.array(out), numpy.nan
327
328 good = numpy.where(numpy.logical_and(numpy.isfinite(science), numpy.any(numpy.isfinite(fringes), 1)))
329 science = science[good]
330 fringes = fringes[good]
331 oldNum = len(science)
332 if oldNum == 0:
333 return emptyResult()
334
335
336
337 good = select(science, self.config.clip)
338 for ff in range(fringes.shape[1]):
339 good &= select(fringes[:, ff], self.config.clip)
340 science = science[good]
341 fringes = fringes[good]
342 oldNum = len(science)
343 if oldNum == 0:
344 return emptyResult(" after initial rejection")
345
346 for i in range(self.config.iterations):
347 solution = self._solve(science, fringes)
348 resid = science - numpy.sum(solution*fringes, 1)
349 rms = stdev(resid)
350 good = numpy.logical_not(abs(resid) > self.config.clip*rms)
351 self.log.debug("Iteration %d: RMS=%f numGood=%d", i, rms, good.sum())
352 self.log.debug("Solution %d: %s", i, solution)
353 newNum = good.sum()
354 if newNum == 0:
355 return emptyResult(" after %d rejection iterations" % i)
356
357 if doPlot:
358 import matplotlib.pyplot as plot
359 for j in range(fringes.shape[1]):
360 fig = plot.figure(j)
361 fig.clf()
362 try:
363 fig.canvas._tkcanvas._root().lift()
364 except Exception:
365 pass
366 ax = fig.add_subplot(1, 1, 1)
367 adjust = science.copy()
368 others =
set(range(fringes.shape[1]))
369 others.discard(j)
370 for k in others:
371 adjust -= solution[k]*fringes[:, k]
372 ax.plot(fringes[:, j], adjust, 'r.')
373 xmin = fringes[:, j].
min()
374 xmax = fringes[:, j].
max()
375 ymin = solution[j]*xmin
376 ymax = solution[j]*xmax
377 ax.plot([xmin, xmax], [ymin, ymax], 'b-')
378 ax.set_title("Fringe %d: %f" % (j, solution[j]))
379 ax.set_xlabel("Fringe amplitude")
380 ax.set_ylabel("Science amplitude")
381 ax.set_autoscale_on(False)
382 ax.set_xbound(lower=xmin, upper=xmax)
383 ax.set_ybound(lower=ymin, upper=ymax)
384 fig.show()
385 while True:
386 ans = input("Enter or c to continue [chp]").lower()
387 if ans in ("", "c",):
388 break
389 if ans in ("p",):
390 import pdb
391 pdb.set_trace()
392 elif ans in ("h", ):
393 print("h[elp] c[ontinue] p[db]")
394
395 if newNum == oldNum:
396
397 break
398 oldNum = newNum
399 good = numpy.where(good)
400 science = science[good]
401 fringes = fringes[good]
402
403
404 solution = self._solve(science, fringes)
405 self.log.info("Fringe solution: %s RMS: %f Good: %d/%d", solution, rms, len(science), origNum)
406 return solution, rms
407
daf::base::PropertySet * set