 News |
FAQ
|
| Search
|
Home
|
| Getting Started | Documentation | Glish | Learn More | Programming | Contact Us |
|
| Version 1.9 Build 1556 |
|
Here we give some examples. These should be considered illustrative rather than production quality code. However the general framework will stand.
The following reads a MeasurementSet and an image model from disk and calibrates the data, first for a simple G-term every 5 minutes and then for a D-term every 12 hours.
// Make a SkyEquation
SkyEquation se;
// Read the VisSet from disk
VisSet vs("3c84.MS");
se.setVisSet(vs);
// Create an ImageSkyModel from an image on disk
ImageSkyModel ism(Image<StokesVector>("3c84"));
se.setSkyModel(ism);
// Use DFT for forward transform, FFT for reverse
IFFTCoh ift;
DFTCoh ft;
se.setFTCoh(ft);
se.setIFTCoh(ift);
// Predict the visibility set
se.predict();
// Make a Vis equation
VisEquation se();
// Solve for calibration of G matrix every 5 minutes
GJones gj(vs, 5*60);
ve.solve(gj);
ve.setVisJones(gj);
// Solve for calibration of D matrix every 12 hours
DJones dj(vs, 12*60*60);
ve.solve(dj);
ve.setVisJones(dj);
Next we give a simple two-dimensional Clean algorithm that finds peaks (actually the peak maximum eigenvalue of the coherence matrix) and subtracts each longwindedly.
// Clean solver
Bool HogbomCleanImageSkyModel::solve(SkyEquation& se) {
// Make a local copy of the Sky Equation so we can change
// some of the entries
SkyEquation lse(se);
// Tell the SkyEquation to use this SkyModel
lse.setSkyModel(*this);
// Predict visibility for this Sky Model
lse.predict ();
// Find the gradients
lse.gradientsChiSquared(*this);
// Make the residual image
makeNewtonRaphsonStep(1.0);
Int nx=image_.shape()(0);
Int ny=image_.shape()(1);
Matrix<StokesVector> limageStep(imageStep_.shape());
Matrix<StokesVector> limage(image_.shape());
Matrix<Float> lpsf(psf_.shape());
IPosition ipStart(image_.ndim(), 0);
IPosition ipStride(image_.ndim(), 1);
imageStep_.getSlice(limageStep, ipStart, imageStep_.shape(), ipStride);
image_.getSlice(limage, ipStart, image_.shape(), ipStride);
psf_.getSlice(lpsf, IPosition(psf_.ndim(),0), psf_.shape(),
IPosition(psf_.ndim(),1));
cout<<"Center of PSF = "<<lpsf(nx/2,ny/2)<<endl;
// Iterate
Float max=0.0;
for (uInt iter=0;iter<numberIterations();iter++) {
// Now find peak in residual image
StokesVector maxVal;
Int ix, iy;
Int px=nx/4;
Int py=ny/4;
max=0.0;
for (iy=ny/4;iy<3*ny/4-1;iy++) {
for (ix=nx/4;ix<3*nx/4-1;ix++) {
if(abs(limageStep(ix,iy).maxEigenValue())>max) {
px=ix;
py=iy;
maxVal=limageStep(ix,iy);
max=abs(limageStep(ix,iy).maxEigenValue());
}
}
}
AlwaysAssert(max>0.0, AipsError);
AlwaysAssert(px>=nx/4&&px<3*nx/4-1,AipsError);
AlwaysAssert(py>=ny/4&&py<3*ny/4-1,AipsError);
// Output ten lines of information if run to the end
Int cycle;
cycle=numberIterations()/10;
if(iter==0||(iter%cycle)==0) {
cout<<"Iteration "<<iter<<" peak is "<<maxVal<<" at "<<px<<","<<py<<endl;
}
if(max<threshold()) {
cout<<"Converged"<<endl;
cout<<"Final iteration "<<iter<<" peak is "<<maxVal<<" at "<<px<<","<<py<<endl;
break;
};
// Add the scaled peak to the current image
StokesVector pv=gain()*maxVal;
limage(px,py)+=pv;
// Subtract the scaled PSF from the residual image
for (iy=ny/4;iy<3*ny/4-1;iy++) {
for (ix=nx/4;ix<3*nx/4-1;ix++) {
limageStep(ix,iy)-=pv*lpsf(nx/2+ix-px,ny/2+iy-py);
}
}
};
imageStep_.putSlice(limageStep, ipStart, ipStride);
image_.putSlice(limage, ipStart, ipStride);
if(max>threshold()) {
cout<<"Did not converge in "<<numberIterations()<<" iterations"<<endl;
return(False);
}
else {
return(True);
};
};
Notes:
.
Next we give an example of the solution for a single fixed ``D'' term per antenna. This actually solves for a general matrix with no assumed symmetries.
// Solve for the Jones matrix. Updates the VisJones thus found.
// Also inserts it into the MeasurementEquation thus it is not const.
Bool TimeVarVisJones::solve (VisEquation& ve)
{
// Make a local copy of the Vis Equation
VisEquation lve(ve);
AlwaysAssert(gain()>0.0,AipsError);
AlwaysAssert(numberIterations()>0,AipsError);
AlwaysAssert(tolerance()>0.0,AipsError);
// Mask for gradient computations
Matrix<Bool> required(2,2);
setMask(required);
// Set the Jones matrix in the Measurement Equation
lve.setVisJones(*this);
// Count number of failed solutions
Int failed=0;
// Tell the MeasurementEquation to use the internal (chunked) VisSet
lve.setVisSet(intvs_);
// Reset (chunked) iterators to start of data
intvs_.observedCoherence().originChunks();
intvs_.correctedCoherence().originChunks();
intvs_.modelCoherence().originChunks();
// Now iterate over the solution intervals
for (currentSlot_=0; currentSlot_<nInterval_; currentSlot_++) {
cout<<"Finding solution for slot "<<currentSlot_<<endl;
// Find gradient and Hessian
lve.gradientsChiSquared(required,*this);
// Assess fit
Float sumwt=0.0;
Float currentFit=0.0;
currentFit=chisq_;
cout<<"Initial fit = "<<currentFit<<endl;
Float originalFit=currentFit;
Float previousFit=currentFit;
// Iterate
for (uInt iter=0;iter<numberIterations();iter++) {
// update antenna gains from gradients
updateAntGain();
updateIntGain();
// Find gradient and Hessian
lve.gradientsChiSquared(required,*this);
// Assess fit
previousFit=currentFit;
currentFit=chisq_;
cout<<"Iteration "<<iter<<" fit = "<<currentFit<<endl;
// Stop?
if(abs(currentFit-previousFit)<tolerance()*originalFit) break;
}
if(abs(currentFit-previousFit)>tolerance()*originalFit) failed++;
// Advance the iterators to the next interval
intvs_.observedCoherence().nextChunk();
intvs_.modelCoherence().nextChunk();
intvs_.correctedCoherence().nextChunk();
}
if(failed>0) {
cout<<"Did not converge in "<<numberIterations()<<" iterations"<<
"for "<<failed<<" time intervals"<<endl;
return(False);
}
else {
return(True);
}
}
Notes: