CEMemModel.h

Classes

CEMemModel -- Implements the Cornwell & Evans MEM Algorithm on Lattices (full description)

class CEMemModel: public LinearModel< Lattice<Float> >

Types

enum GradientType

H = 0
Entropy
C = 1
Chi_sq
F = 2
Flux
J = 3
Objective function J

Interface

Public Members
CEMemModel(Entropy &ent, Lattice<Float> & model, uInt nIntegrations = 10, Float sigma = 0.001, Float targetFlux = 1.0, Bool useFluxConstraint = False, Bool initializeModel = True, Bool imagePlane = False)
CEMemModel(Entropy &ent, Lattice<Float> & model, Lattice<Float> & prior, uInt nIntegrations = 10, Float sigma = 0.001, Float targetFlux = 1.0, Bool useFluxConstraint = False, Bool initializeModel = True, Bool imagePlane = False)
CEMemModel(Entropy &ent, Lattice<Float> & model, Lattice<Float> & prior, Lattice<Float> & mask, uInt nIntegrations = 10, Float sigma = 0.001, Float targetFlux = 1.0, Bool useFluxConstraint = False, Bool initializeModel = True, Bool imagePlane = False)
virtual ~CEMemModel()
void state()
Bool solve(ResidualEquation<Lattice<Float> > & eqn)
void setEntropy(Entropy &myEntropy )
void getEntropy(Entropy &myEntropy )
Lattice<Float>& getModel() const
void setModel(const Lattice<Float> & model)
Lattice<Float>& getPrior() const
void setPrior(Lattice<Float> & prior)
Lattice<Float>& getMask() const
void setMask(Lattice<Float> & mask)
Float getAlpha() const
Float getBeta() const
void setAlpha(Float alpha)
void setBeta(Float beta)
Float getTolerance()
void setTolerance(Float x)
Float getQ()
void setQ(Float x)
Float getGain()
void setGain(Float x)
Float getMaxNormGrad()
void setMaxNormGrad(Float x)
Int getInitialNumberIterations()
void setInitialNumberIterations(Int x)
void setThreshold (const Float x )
void setThresholdSpeedup (const Float iter)
Float getThreshold()
virtual void setProgress(CEMemProgress& ccp)
virtual CEMemProgress& getProgress()
Int numberIterations ()
void setCycleFlux(Float x)
Float getCycleFlux()
Protected Members
void oneIteration()
Bool applyMask( Lattice<Float> & lat )
void formEntropy()
void formGDG()
void formGDGStep()
void formGDS()
void entropyType(String &str)
void relaxMin()
Bool testConvergence()
CEMemModel()
void letEntropyDie()
Bool initStuff()
void changeAlphaBeta ()
void initializeAlphaBeta()
void updateAlphaBeta()
Bool checkImage(const Lattice<Float> *)
Bool checkImages(const Lattice<Float> *one, const Lattice<Float> *other)
Bool ok()
void calculateStep()
void takeStep(Float wt1, Float wt2)
Float formFlux()
Bool testConvergenceThreshold()

Description

Review Status

Date Reviewed:
yyyy/mm/dd

Prerequisite

Etymology

This class is called CEMemModel because it uses the Cornwell and Evans MEM algorithm to deconvolve the model.

Synopsis

This class is used to perform the Cornwell and Evans MEM Algorithm on an Array. Only the deconvolved model of the sky are directly stored by this class. The point spread function (psf) and convolved (dirty) image are stored in a companion class which is must be derived from ResidualEquation.

The deconvolution works like this. The user constructs a CEMemModel by specifying an initial model of the sky. This can by be one,two,three... dimensional depending on the dimension of the psf (see below). The user then constructs a class which implements the forward equation between the model and the dirty image. Typically this will be the ConvolutionEquation class, although any class which has a ResidualEquation interface will work (but perhaps very slowly, as the ConvolutionEquation class has member functions optimised for CLEAN and MEM)

The user then calls the solve() function (with the appropriate equation class as an arguement), and this class will perform the MEM deconvolution. The various MEM parameters are set (prior to calling solve) using the functions derived from the Iterate class, in particular setNumberIterations().

The solve() function does not return either the deconvolved model or the residuals. The solved model can be obtained using the getModel() function (derived from ArrayModel()) and the residual can be obtained using the residual() member function of the Convolution/Residual Equation Class.

The size and shape of the model used in this class MUST be the same as the convolved data (Dirty Image), stored in the companion ResidualEquation Class. However the model (and convolved data) can have more dimensions than the psf, as well as a different size (either larger or smaller). When the dimensionality is different the deconvolution is done independendtly in each "plane" of the model. (Note this has not been implemented yet but is relatively simple to do if necessary).

StokesVectors are not yet implemented.

A companion class to this one is MaskedCEMemModel. This provides the same functionality but is used with MaskedArrays which indicate which regions of the NewtonRaphson (residual) image to apply when forming the step image (MaskedCEMemModel is not yet implemented).

Example

    Matrix<Float> psf(12,12), dirty(10,10), initialModel(10,10);
    ...put appropriate values into psf, dirty, & initialModel....
    CEMemModel<Float> deconvolvedModel(initialModel); 
    ConvolutionEquation convEqn(psf, dirty);
    deconvolvedModel.setSigma(0.001); 
    deconvolvedModel.setTargetFlux(-2.500); 
    deconvolvedModel.setNumberIterations(30);
    Bool convWorked = deconvolvedModel.solve(convEqn);
    Array<Float> finalModel, residuals;
    if (convWorked){
      finalModel = deconvolvedModel.getModel();
      ConvEqn.residual(deconvolvedModel, finalResidual);
    }
    

Motivation

This class is needed to deconvolve extended images. In time, the MEM algorithm will be a principle player in the mosaicing stuff.

Template Type Argument Requirements (T)

To Do

Member Description

CEMemModel(Entropy &ent, Lattice<Float> & model, uInt nIntegrations = 10, Float sigma = 0.001, Float targetFlux = 1.0, Bool useFluxConstraint = False, Bool initializeModel = True, Bool imagePlane = False)

Construct the CEMemModel object and initialise the model.

CEMemModel(Entropy &ent, Lattice<Float> & model, Lattice<Float> & prior, uInt nIntegrations = 10, Float sigma = 0.001, Float targetFlux = 1.0, Bool useFluxConstraint = False, Bool initializeModel = True, Bool imagePlane = False)

Construct the CEMemModel object, initialise the model and Prior images.

CEMemModel(Entropy &ent, Lattice<Float> & model, Lattice<Float> & prior, Lattice<Float> & mask, uInt nIntegrations = 10, Float sigma = 0.001, Float targetFlux = 1.0, Bool useFluxConstraint = False, Bool initializeModel = True, Bool imagePlane = False)

Construct the CEMemModel object, initialise the model, Prior, and mask images.

virtual ~CEMemModel()

destructor

void state()

solve the convolution equation returns True if converged

Gives information about the state of the CEMem (ie, using mask image, using prior image; more work here!)

Bool solve(ResidualEquation<Lattice<Float> > & eqn)

This needs to be "ResidualEquation", using LatConvEquation as polymorphism is broken

void setEntropy(Entropy &myEntropy )

Set and get various state images and classes

set or get the Entropy class

Lattice<Float>& getModel() const

Set and get various state images and classes

set or get the Model image

Lattice<Float>& getPrior() const

Set and get various state images and classes

set or get the Prior image

Lattice<Float>& getMask() const

Set and get various state images and classes

set or get the Mask image

void getEntropy(Entropy &myEntropy )
void setModel(const Lattice<Float> & model)
void setPrior(Lattice<Float> & prior)
void setMask(Lattice<Float> & mask)

Set and get various state images and classes

Float getAlpha() const
Float getBeta() const
void setAlpha(Float alpha)
void setBeta(Float beta)

set and get alpha and beta

Float getTolerance()
void setTolerance(Float x)
Float getQ()
void setQ(Float x)
Float getGain()
void setGain(Float x)
Float getMaxNormGrad()
void setMaxNormGrad(Float x)
Int getInitialNumberIterations()
void setInitialNumberIterations(Int x)

Set various controlling parameters (The most popular controlling parameters are set in the constructor)

void setThreshold (const Float x )

The convergence can also be in terms of the maximum residual (ie, for use in stopping in a major cycle).

void setThresholdSpeedup (const Float iter)

Thresh doubles in iter iterations

Float getThreshold()

virtual void setProgress(CEMemProgress& ccp)
virtual CEMemProgress& getProgress()

Set/get the progress display

Int numberIterations ()

return the number of iterations completed

void setCycleFlux(Float x)

if this MEMModel is constructed in a MF loop, we may need to increment the flux by the last iterations flux

Float getCycleFlux()

void oneIteration()

Perform One Iteration

Bool applyMask( Lattice<Float> & lat )

apply mask to a lattice; returns True if mask is available, False if not

void formEntropy()
void formGDG()
void formGDGStep()
void formGDS()
void entropyType(String &str)
void relaxMin()
Bool testConvergence()

Helper functions that interface with Entropy routines Access to the entropy should be through this interface; the points at which the Entropy is mentioned is then limited to these lines right here, and to the constructors which set the Entropy. The entropy should not ever change the private

CEMemModel()

protected generic constrcutor: DON'T USE IT!

void letEntropyDie()

Set entropy pointer to zero: called by Entropy's destructor

Bool initStuff()

initialize itsStep and itsResidual and other stuff

void initializeAlphaBeta()

controls how to change Alpha and Beta

initialize Alpha-Beta (called by changeAlphaBeta)

void updateAlphaBeta()

controls how to change Alpha and Beta

update Alpha-Beta (called by changeAlphaBeta)

void changeAlphaBeta ()

controls how to change Alpha and Beta

Bool checkImages(const Lattice<Float> *one, const Lattice<Float> *other)

Generic utility functions

check that the lattices and the underlying tile sizes are consistent

Bool ok()

Generic utility functions

check that all is well in Denmark: ensure all images are the same size, ensure we have an entropy, ensure state variables have reasonable values

Bool checkImage(const Lattice<Float> *)

Generic utility functions

void takeStep(Float wt1, Float wt2)

Helper functions for oneIteration:

take one step: clipped addition of wt1*itsModel + wt2*itsStep

Float formFlux()

Helper functions for oneIteration:

Calculate the flux, itsModMin, and itsModMax

void calculateStep()

Helper functions for oneIteration:

Bool testConvergenceThreshold()

Determine if the peak residual is less than the getThreshold()

enum GradientType

Enumerate the different Gradient subscript types