# Soften symmetry operator requirements

XMLWordPrintable

#### Details

• Type: Story
• Status: Done
• Resolution: Done
• Fix Version/s: None
• Component/s:
• Labels:
None
• Story Points:
1
• Sprint:
DRP S17-4
• Team:
Data Release Production

#### Description

The current symmetry operator enforces strict symmetry, meaning that all pixels that do not have a symmetric partner in the footprint are penalized. To ease this requirement, which may be necessary if some of the flux for one (or more) of the objects lies outside of the footprint, we allow the user to use a less stringent penalty value.

#### Activity

Hide
Fred Moolekamp added a comment -

For simplicity this work is not being done on a separate branch, but instead will be included in DM-9172. The code for the new symmetry operator appears below:

 def getPeakSymmetry(shape, px, py, fillValue=0):  """Build the operator to symmetrize a the intensities for a single row  """  center = (np.array(shape)-1)/2.0  # If the peak is centered at the middle of the footprint,  # make the entire footprint symmetric  if px==center[1] and py==center[0]:  return scipy.sparse.coo_matrix(np.fliplr(np.eye(shape[0]*shape[1])))    # Otherwise, find the bounding box that contains the minimum number of pixels needed to symmetrize  if py<(shape[0]-1)/2.:  ymin = 0  ymax = 2*py+1  elif py>(shape[0]-1)/2.:  ymin = 2*py-shape[0]+1  ymax = shape[0]  else:  ymin = 0  ymax = shape[0]  if px<(shape[1]-1)/2.:  xmin = 0  xmax = 2*px+1  elif px>(shape[1]-1)/2.:  xmin = 2*px-shape[1]+1  xmax = shape[1]  else:  xmin = 0  xmax = shape[1]    fpHeight, fpWidth = shape  fpSize = fpWidth*fpHeight  tWidth = xmax-xmin  tHeight = ymax-ymin  extraWidth = fpWidth-tWidth  pixels = (tHeight-1)*fpWidth+tWidth    # This is the block of the matrix that symmetrizes intensities at the peak position  subOp = np.eye(pixels, pixels)  for i in range(0,tHeight-1):  for j in range(extraWidth):  idx = (i+1)*tWidth+(i*extraWidth)+j  subOp[idx, idx] = fillValue  subOp = np.fliplr(subOp)    smin = ymin*fpWidth+xmin  smax = (ymax-1)*fpWidth+xmax  if fillValue!=0:  symmetryOp = np.identity(fpSize)*fillValue  else:  symmetryOp = np.zeros((fpSize, fpSize))  symmetryOp[smin:smax,smin:smax] = subOp    # Return a sparse matrix, which greatly speeds up the processing  return scipy.sparse.coo_matrix(symmetryOp)   def getPeakSymmetryOp(shape, px, py, fillValue=0):  """Operator to calculate the difference from the symmetric intensities  """  symOp = getPeakSymmetry(shape, px, py, fillValue)  diffOp = scipy.sparse.identity(symOp.shape[0])-symOp  # In cases where the symmetry operator is very small (eg. a small isolated source)  # scipy doesn't return a sparse matrix, so we test whether or not the matrix is sparse  # and if it is, use a sparse matrix that works best with the proximal operators.  if hasattr(diffOp, "tocoo"):  diffOp = diffOp.tocoo()  return diffOp 

Show
Fred Moolekamp added a comment - For simplicity this work is not being done on a separate branch, but instead will be included in DM-9172 . The code for the new symmetry operator appears below: def getPeakSymmetry(shape, px, py, fillValue = 0 ): """Build the operator to symmetrize a the intensities for a single row """ center = (np.array(shape) - 1 ) / 2.0 # If the peak is centered at the middle of the footprint, # make the entire footprint symmetric if px = = center[ 1 ] and py = = center[ 0 ]: return scipy.sparse.coo_matrix(np.fliplr(np.eye(shape[ 0 ] * shape[ 1 ])))   # Otherwise, find the bounding box that contains the minimum number of pixels needed to symmetrize if py<(shape[ 0 ] - 1 ) / 2. : ymin = 0 ymax = 2 * py + 1 elif py>(shape[ 0 ] - 1 ) / 2. : ymin = 2 * py - shape[ 0 ] + 1 ymax = shape[ 0 ] else : ymin = 0 ymax = shape[ 0 ] if px<(shape[ 1 ] - 1 ) / 2. : xmin = 0 xmax = 2 * px + 1 elif px>(shape[ 1 ] - 1 ) / 2. : xmin = 2 * px - shape[ 1 ] + 1 xmax = shape[ 1 ] else : xmin = 0 xmax = shape[ 1 ]   fpHeight, fpWidth = shape fpSize = fpWidth * fpHeight tWidth = xmax - xmin tHeight = ymax - ymin extraWidth = fpWidth - tWidth pixels = (tHeight - 1 ) * fpWidth + tWidth   # This is the block of the matrix that symmetrizes intensities at the peak position subOp = np.eye(pixels, pixels) for i in range ( 0 ,tHeight - 1 ): for j in range (extraWidth): idx = (i + 1 ) * tWidth + (i * extraWidth) + j subOp[idx, idx] = fillValue subOp = np.fliplr(subOp)   smin = ymin * fpWidth + xmin smax = (ymax - 1 ) * fpWidth + xmax if fillValue! = 0 : symmetryOp = np.identity(fpSize) * fillValue else : symmetryOp = np.zeros((fpSize, fpSize)) symmetryOp[smin:smax,smin:smax] = subOp   # Return a sparse matrix, which greatly speeds up the processing return scipy.sparse.coo_matrix(symmetryOp)   def getPeakSymmetryOp(shape, px, py, fillValue = 0 ): """Operator to calculate the difference from the symmetric intensities """ symOp = getPeakSymmetry(shape, px, py, fillValue) diffOp = scipy.sparse.identity(symOp.shape[ 0 ]) - symOp # In cases where the symmetry operator is very small (eg. a small isolated source) # scipy doesn't return a sparse matrix, so we test whether or not the matrix is sparse # and if it is, use a sparse matrix that works best with the proximal operators. if hasattr (diffOp, "tocoo" ): diffOp = diffOp.tocoo() return diffOp
Hide
Fred Moolekamp added a comment -

Hi Peter,
Would you mind taking a look at the new symmetry operator and testing it out?

It requires you to add a new optional parameter fillValue=0 to nmf_deblender, and call getPeakSymmetryOp using C = getPeakSymmetryOp((N,M), px, py, fillValue). fillValue is the inverse of the penalty to use, so fillValue=0 is the current implementation that strictly enforces symmetry, while fillValue=1 doesn't penalize pixels outside of the symmetric patch at all.

Let me know if you have any problems or questions.

-Fred

Show
Fred Moolekamp added a comment - Hi Peter, Would you mind taking a look at the new symmetry operator and testing it out? It requires you to add a new optional parameter fillValue=0 to nmf_deblender , and call getPeakSymmetryOp using C = getPeakSymmetryOp((N,M), px, py, fillValue) . fillValue is the inverse of the penalty to use, so fillValue=0 is the current implementation that strictly enforces symmetry, while fillValue=1 doesn't penalize pixels outside of the symmetric patch at all. Let me know if you have any problems or questions. -Fred
Hide
Peter Melchior added a comment -

Works as expected. Interesting to note: fillValue can be some float, but it's only effective in a narrow range around 1.

Show
Peter Melchior added a comment - Works as expected. Interesting to note: fillValue can be some float, but it's only effective in a narrow range around 1.

#### People

Assignee:
Fred Moolekamp
Reporter:
Fred Moolekamp
Reviewers:
Peter Melchior
Watchers:
Fred Moolekamp, Peter Melchior