LCOV - code coverage report
Current view: top level - synthesis/MeasurementEquations - MultiTermMatrixCleaner.cc (source / functions) Hit Total Coverage
Test: casacpp_coverage.info Lines: 0 534 0.0 %
Date: 2024-10-04 18:58:15 Functions: 0 30 0.0 %

          Line data    Source code
       1             : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
       2             : //# Associated Universities, Inc. Washington DC, USA.
       3             : //#
       4             : //# This library is free software; you can redistribute it and/or modify it
       5             : //# under the terms of the GNU Library General Public License as published by
       6             : //# the Free Software Foundation; either version 2 of the License, or (at your
       7             : //# option) any later version.
       8             : //#
       9             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      10             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      11             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      12             : //# License for more details.
      13             : //#
      14             : //# You should have received a copy of the GNU Library General Public License
      15             : //# along with this library; if not, write to the Free Software Foundation,
      16             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      17             : //#
      18             : //# Correspondence concerning AIPS++ should be addressed as follows:
      19             : //#        Internet email: casa-feedback@nrao.edu.
      20             : //#        Postal address: AIPS++ Project Office
      21             : //#                        National Radio Astronomy Observatory
      22             : //#                        520 Edgemont Road
      23             : //#                        Charlottesville, VA 22903-2475 USA
      24             : //#
      25             : //# $Id: MultiTermMatrixCleaner.cc 13656 2010-12-04 02:08:02 UrvashiRV$
      26             : //
      27             : // Implementation of the paper "A multi-scale multi-frequency deconvolution algorithm
      28             : // for synthesis imaging in radio interferometry" by Rau and Cornwell
      29             : // https://www.aanda.org/articles/aa/pdf/2011/08/aa17104-11.pdf
      30             : 
      31             : // Same include list as in MatrixCleaner.cc
      32             : #include <casacore/casa/Arrays/Matrix.h>
      33             : #include <casacore/casa/Arrays/Cube.h>
      34             : #include <casacore/casa/Arrays/ArrayMath.h>
      35             : #include <casacore/casa/Arrays/MatrixMath.h>
      36             : #include <casacore/casa/IO/ArrayIO.h>
      37             : #include <casacore/casa/BasicMath/Math.h>
      38             : #include <casacore/casa/BasicSL/Complex.h>
      39             : #include <casacore/casa/Logging/LogIO.h>
      40             : #include <casacore/casa/OS/File.h>
      41             : #include <casacore/casa/Containers/Record.h>
      42             : 
      43             : #include <casacore/lattices/LRegions/LCBox.h>
      44             : #include <casacore/casa/Arrays/Slicer.h>
      45             : #include <casacore/scimath/Mathematics/FFTServer.h>
      46             : #include <casacore/casa/OS/HostInfo.h>
      47             : #include <casacore/casa/Arrays/ArrayError.h>
      48             : #include <casacore/casa/Arrays/ArrayIter.h>
      49             : #include <casacore/casa/Arrays/VectorIter.h>
      50             : 
      51             : 
      52             : #include <casacore/casa/Utilities/GenSort.h>
      53             : #include <casacore/casa/BasicSL/String.h>
      54             : #include <casacore/casa/Utilities/Assert.h>
      55             : #include <casacore/casa/Utilities/Fallible.h>
      56             : 
      57             : #include <casacore/casa/BasicSL/Constants.h>
      58             : 
      59             : #include <casacore/casa/Logging/LogSink.h>
      60             : #include <casacore/casa/Logging/LogMessage.h>
      61             : 
      62             : #include <synthesis/MeasurementEquations/MatrixCleaner.h>
      63             : #include <casacore/coordinates/Coordinates/TabularCoordinate.h>
      64             : 
      65             : // Additional include files
      66             : #include <casacore/lattices/Lattices/SubLattice.h>
      67             : #include <casacore/scimath/Mathematics/MatrixMathLA.h>
      68             : #include <casacore/images/Images/PagedImage.h>
      69             : 
      70             : #include<synthesis/MeasurementEquations/MultiTermMatrixCleaner.h>
      71             : 
      72             : using namespace casacore;
      73             : namespace casa { //# NAMESPACE CASA - BEGIN
      74             : 
      75             : #define MIN(a,b) ((a)<=(b) ? (a) : (b))
      76             : #define MAX(a,b) ((a)>=(b) ? (a) : (b))
      77             :         
      78           0 :   MultiTermMatrixCleaner::MultiTermMatrixCleaner():
      79             :     MatrixCleaner(),
      80           0 :     ntaylor_p(0),psfntaylor_p(0),nscales_p(0),nx_p(0),ny_p(0),totalIters_p(0),
      81           0 :     maxscaleindex_p(0), globalmaxpos_p(IPosition(0)),
      82           0 :     /*donePSF_p(false),*/donePSP_p(false),doneCONV_p(false),memoryMB_p(0),
      83           0 :     adbg(false)
      84           0 :   { }
      85             : 
      86             :   /*
      87             : MultiTermMatrixCleaner::MultiTermMatrixCleaner(const MultiTermMatrixCleaner & other):
      88             :   MatrixCleaner(),ntaylor_p(other.ntaylor_p)  //And others... minus some...
      89             :     {
      90             :     }
      91             : 
      92             : MultiTermMatrixCleaner &MultiTermMatrixCleaner:: operator=(const MultiTermMatrixCleaner & other) 
      93             : {
      94             :   operator=(other);
      95             :   return *this;
      96             : }
      97             :   */
      98             : 
      99           0 :  MultiTermMatrixCleaner::~MultiTermMatrixCleaner()
     100             : {
     101             :   //cout << "In MultiTermMatrixCleaner destructor" << endl;
     102           0 : }
     103             : 
     104             : 
     105           0 : Bool MultiTermMatrixCleaner::setscales(const Vector<Float> & scales)
     106             : {
     107           0 :    nscales_p = scales.nelements();
     108           0 :    scaleSizes_p.resize();
     109           0 :    scaleSizes_p = scales;
     110           0 :    totalScaleFlux_p.resize(nscales_p);
     111           0 :    totalScaleFlux_p.set(0.0);
     112           0 :    maxScaleVal_p.resize(nscales_p);
     113           0 :    maxScaleVal_p.set(0.0);
     114           0 :    maxScalePos_p.resize(nscales_p);
     115           0 :    globalmaxval_p=-1e+10;
     116           0 :    globalmaxpos_p=IPosition(2,0);
     117           0 :    maxscaleindex_p=0;
     118           0 :    return true;
     119             : }
     120             : 
     121             : 
     122           0 : Bool MultiTermMatrixCleaner::setntaylorterms(const int & nterms)
     123             : {
     124           0 :    ntaylor_p = nterms;
     125           0 :    psfntaylor_p = 2*nterms-1;
     126           0 :    totalTaylorFlux_p.resize(ntaylor_p);
     127           0 :    totalTaylorFlux_p.set(0.0);
     128           0 :    return true;
     129             : }
     130             : 
     131             : // Allocate memory, based on nscales and ntaylor
     132           0 : Bool MultiTermMatrixCleaner::initialise(Int nx, Int ny)
     133             : {
     134           0 :   LogIO os(LogOrigin("MultiTermMatrixCleaner", "initialise()", WHERE));
     135             :   
     136             :   /* Verify Image Shapes */
     137             :   //  nx_p = model.shape(0);
     138           0 :   nx_p = nx;
     139           0 :   ny_p = ny;
     140             :   
     141           0 :   if(adbg) os << "Checking shapes" << LogIO::POST;
     142             :   
     143             :   /* Verify nscales_p and ntaylor_p */
     144           0 :   AlwaysAssert(nscales_p>0, AipsError);
     145           0 :   AlwaysAssert(ntaylor_p>0, AipsError);
     146           0 :   AlwaysAssert(scaleSizes_p.nelements()==static_cast<uInt>(nscales_p), AipsError);
     147             : 
     148           0 :   if(adbg) os << "Verify Scale sizes" << LogIO::POST;
     149           0 :   verifyScaleSizes();
     150             : 
     151             :   /* Calculate PSF/scale support - after verifying scale sizes and before allocating memory  */
     152             : 
     153             :   // PSF support : main lobe width x how many times this width the patch should cover
     154           0 :   Float maxscalesize = scaleSizes_p[nscales_p-1];
     155             : 
     156             :   //// N times the size of the main lobe of the PSF at the max scale size
     157           0 :   Float psfbeam = 4.0;
     158           0 :   Float nbeams = 20.0;
     159           0 :   Int psupport = findBeamPatch(maxscalesize, nx_p, ny_p, psfbeam, nbeams);
     160             :   // //----------------------Encapsulated in the findBeamPatch() method from here...-----------
     161             :   // Int psupport = (Int) ( sqrt( psfbeam*psfbeam + maxscalesize*maxscalesize ) * nbeams  );
     162             : 
     163             :   // // At least this big...
     164             :   // if(psupport < psfbeam*nbeams ) psupport = static_cast<Int>(psfbeam*nbeams);
     165             : 
     166             :   // // Not too big...
     167             :   // if(psupport > nx_p || psupport > ny_p)   psupport = MIN(nx_p,ny_p);
     168             : 
     169             :   // // Make it even.
     170             :   // if (psupport%2 != 0) psupport -= 1;
     171             :   // //----------------------...till here------------------------------------
     172             : 
     173             : 
     174             :   // Full image
     175             :   //psfsupport_p = IPosition(2,MIN(nx_p,ny_p),MIN(nx_p,ny_p));
     176             :   //psfpeak_p = IPosition(2,nx_p/2, ny_p/2);
     177             : 
     178             :   // Inner quater image
     179             :   //psfsupport_p = IPosition(2,MIN(nx_p/2,ny_p/2),MIN(nx_p/2,ny_p/2));
     180             :   //psfpeak_p = IPosition(2,nx_p/4, ny_p/4);
     181             : 
     182             :   // Generic support-size.
     183           0 :   psfsupport_p = IPosition(2,psupport,psupport);
     184           0 :   psfpeak_p = IPosition(2,psupport/2, psupport/2);
     185             : 
     186           0 :   os << "Using a PSF patch of " << psupport << " pixels on each side for minor-cycle updates." << endl;
     187             : 
     188             :   /* Force the scale images to be */
     189             :   
     190             :   
     191           0 :   if(adbg) os << "Start allocating mem" << LogIO::POST;
     192             : 
     193             :   /* Allocate memory for many many TempMatrices. */
     194             :   /* Check on a return value of -1 - for insuffucient memory */
     195           0 :   if( allocateMemory() == -1 ) return false;
     196             : 
     197             :   /* FFTServer */
     198           0 :   fftcomplex = FFTServer<Float,Complex>(IPosition(2,nx_p,ny_p));
     199             :   
     200             :   /* Create the scale functions and their FTs */
     201           0 :   setupScaleFunctions();
     202             : 
     203           0 :   totalIters_p=0;
     204           0 :   prev_max_p = 1e+10;
     205           0 :   min_max_p = 1e+10;
     206           0 :   rmaxval_p = -1.0;
     207             :   
     208           0 :   if(adbg) os << "Finished initializing MultiTermMatrixCleaner" << LogIO::POST;
     209           0 :   return true;
     210           0 : }
     211             : 
     212             : 
     213           0 : Bool MultiTermMatrixCleaner::buildImagePatches()
     214             : {
     215             : 
     216           0 :    gip = IPosition(2,nx_p,ny_p);  
     217             : 
     218             :    /* The update region. */
     219           0 :    IPosition inc(2,1,1);
     220           0 :    blc_p = IPosition(globalmaxpos_p - psfsupport_p/2);
     221           0 :    trc_p = IPosition(globalmaxpos_p + psfsupport_p/2 - IPosition(2,1,1));
     222             :    //cout << "residual box 1 : " << blc << trc << endl;
     223           0 :    LCBox::verify(blc_p, trc_p, inc, gip);
     224             :    //cout << "residual box 2 : " << blc << trc << endl;
     225             : 
     226             :    
     227             :    /* Shifted region, with the psf at the globalmaxpos_p. */
     228           0 :    blcPsf_p = IPosition(blc_p + psfpeak_p - globalmaxpos_p); // OLD
     229           0 :    trcPsf_p = IPosition(trc_p + psfpeak_p - globalmaxpos_p); // OLD
     230             :    //cout << "psf box 1 : " << blcPsf << trcPsf << endl;
     231           0 :    LCBox::verify(blcPsf_p, trcPsf_p, inc, psfsupport_p); // NEW
     232             :    //cout << "psf box 2 : " << blcPsf << trcPsf << endl;
     233             : 
     234             :    /* Reconcile box sizes/locations with the image size */
     235           0 :    makeBoxesSameSize(blc_p,trc_p,blcPsf_p,trcPsf_p);
     236             : 
     237             :    //cout << "maxpos : " << globalmaxpos_p <<  "  blc : " << blc_p << " trc : " << trc_p << "  blcPsf : " << blcPsf_p << " trcPsf : " << trcPsf_p << endl;
     238             : 
     239           0 :    return true;
     240           0 : }
     241             : 
     242             : 
     243           0 : Bool MultiTermMatrixCleaner::setpsf(int order, Matrix<Float> & psf)
     244             : {
     245           0 :         AlwaysAssert((order>=(int)0 && order<(int)vecPsfFT_p.nelements()), AipsError);
     246           0 :         if(order==0) AlwaysAssert(validatePsf(psf), AipsError);
     247             :         // Directly store the FFT of the PSFs.
     248             :         // The FT'd matrix is reshaped automatically.
     249           0 :         fftcomplex.fft0(vecPsfFT_p[order],psf,false);
     250             :         //fftcomplex.flip(vecPsfFT_p[order],false,false);
     251             :         //        Matrix<Float> temp(real(vecPsfFT_p[order]));
     252             :         //        writeMatrixToDisk("psfFT_"+String::toString(order), temp);
     253           0 :         return true;
     254             : }
     255             : 
     256             : /* Input : Dirty Images */
     257           0 : Bool MultiTermMatrixCleaner::setresidual(int order, Matrix<Float> & dirty)
     258             : {
     259           0 :         AlwaysAssert((order>=(int)0 && order<(int)vecDirty_p.nelements()), AipsError);
     260           0 :         vecDirty_p[order].reference(dirty);
     261           0 :   return true;
     262             : }
     263             : 
     264             : /* Input : Model Component Image */
     265             : //TODO : This is an extra copy of the model image. Why not hold it by reference ???
     266           0 : Bool MultiTermMatrixCleaner::setmodel(int order, Matrix<Float> & model)
     267             : {
     268           0 :         AlwaysAssert((order>=(int)0 && order<(int)vecModel_p.nelements()), AipsError);
     269           0 :         vecModel_p[order].assign(model);
     270           0 :         vecInitialModel_p[order].assign(model);
     271           0 :         totalTaylorFlux_p[order] = (sum(vecModel_p[order]));
     272           0 :   return true;
     273             : }
     274             : 
     275             : /* Input : Mask */
     276             : 
     277           0 : Bool MultiTermMatrixCleaner::setmask(Matrix<Float> & mask)
     278             : {
     279           0 :   if(itsMask.null()) 
     280           0 :        itsMask = new Matrix<Float>(mask);  // it's a counted ptr
     281             :   else
     282             :     {
     283           0 :       AlwaysAssert(itsMask->shape()==mask.shape(), AipsError);
     284           0 :       (*itsMask).assign(mask);
     285             :     }
     286           0 :   return true;
     287             : }
     288             : 
     289             : /* Output : Model Component Image */
     290             : 
     291           0 : Bool MultiTermMatrixCleaner::getmodel(int order, Matrix<Float> & model)
     292             : {
     293           0 :         AlwaysAssert((order>=(int)0 && order<(int)vecModel_p.nelements()), AipsError);
     294           0 :         model.assign(vecModel_p[order]);
     295           0 :   return true;
     296             : }
     297             : 
     298             : /* Output Residual Image */
     299             : 
     300             : 
     301           0 : Bool MultiTermMatrixCleaner::getresidual(int order, Matrix<Float> & residual)
     302             : {
     303           0 :         AlwaysAssert((order>=(int)0 && order<(int)vecDirty_p.nelements()), AipsError);
     304             :         //AlwaysAssert(residual, AipsError);
     305           0 :         residual.assign(vecDirty_p[order]);
     306             :         ///residual.assign( matR_p[IND2(order,0)]  ); // This is the one that gets updated during iters.
     307             : 
     308           0 :   return true;
     309             : }
     310             : 
     311             : 
     312             : 
     313             :  /* Output Hessian matrix */
     314           0 : Bool MultiTermMatrixCleaner::getinvhessian(Matrix<Double> & invhessian)
     315             : {
     316           0 :         invhessian.resize((invMatA_p[0]).shape());
     317           0 :         invhessian = (invMatA_p[0]); 
     318           0 :   return true;
     319             : }
     320             : 
     321             : 
     322             : /* Do the deconvolution */
     323           0 : Int MultiTermMatrixCleaner::mtclean(Int maxniter, Float stopfraction, Float inputgain, Float userthreshold)
     324             : { 
     325           0 :   LogIO os(LogOrigin("MultiTermMatrixCleaner", "mtclean()", WHERE));
     326           0 :   if(adbg)os << "SOLVER for Multi-Frequency Synthesis deconvolution" << LogIO::POST;
     327             : 
     328             :   /* Initialize some variables */  
     329             : 
     330           0 :  maxniter_p = maxniter;
     331           0 :  stopfraction_p=  stopfraction;
     332           0 :  inputgain_p=inputgain;
     333           0 :  userthreshold_p=userthreshold;
     334             : 
     335             :  //cout << "MTMC : User threshold : " << userthreshold_p << endl;
     336             : 
     337             : 
     338           0 :   Int convergedflag = 0;
     339           0 :   Float fluxlimit = -1;
     340           0 :   Float loopgain = 0.5; 
     341           0 :   if(inputgain>(Float)0.0) loopgain=inputgain;
     342           0 :   Int criterion=5; // This chooses among a list of 'penalty functions'
     343           0 :   itercount_p=0;
     344           0 :   Int iterdone = 0;
     345             :  
     346             :   /* Compute all convolutions and the matrix A */
     347             :   /* If matrix is not invertible, return ! */
     348             :   /* This is done only once - for the first major cycle. Null operation otherwise */
     349           0 :   if( doneCONV_p == false )
     350             :     {
     351           0 :       if( computeHessianPeak() == -2 )
     352           0 :         return -2;
     353             :     }
     354             : 
     355             :   /* Set up the Mask image */
     356             :   /* This must be after the Hessian is computed.. */
     357           0 :   setupUserMask();
     358             :    
     359             :   /* Compute the convolutions of the current residual with all PSFs and scales */
     360           0 :   os << "Calculating convolutions of residual images with scales " << LogIO::POST;
     361           0 :   computeRHS();
     362             : 
     363             :   /* Check if peak RHS is already less than the threshold */
     364             :   /* This step computes the flux limit also - when called here... */
     365           0 :   convergedflag = checkConvergence(criterion,fluxlimit,loopgain);
     366           0 :   if(convergedflag==1) return 0;
     367             :   
     368             :   /* Compute the flux limits that determine the depth of the minor cycles. */
     369             :   //  computeFluxLimit(fluxlimit,thresh);
     370             :   
     371             :  /********************** START MINOR CYCLE ITERATIONS ***********************/
     372             :   //os << "Doing deconvolution iterations..." << LogIO::POST;
     373             :   
     374           0 :   for(itercount_p=0;itercount_p<maxniter_p;itercount_p++)
     375             :   {
     376           0 :       globalmaxval_p=-1e+10;
     377             : 
     378             :       // if itercount_p==0, set the blc/trc to the full image size.
     379             :       // For later iterations, it will get updated according to globalmaxpos_p and psfsupport_p
     380             :       //                              when buildImagePatches() is called.
     381           0 :       if(itercount_p==0)
     382             :         {
     383           0 :           blc_p = IPosition(2,0,0);
     384           0 :           trc_p = IPosition(2,nx_p-1,ny_p-1);
     385             :         }
     386             : 
     387           0 :       Int scale=0;
     388           0 :       Int ntaylor=ntaylor_p;
     389           0 :       IPosition blc(blc_p), trc(trc_p);
     390             :       //OMP//      #pragma omp parallel default(shared) private(scale) firstprivate(ntaylor,criterion,blc,trc)
     391             :        { 
     392             :          //OMP//         #pragma omp for 
     393           0 :           for(scale=0;scale<nscales_p;scale++)
     394             :           {
     395             :             /* Solve the matrix eqn for all pixels */
     396           0 :             solveMatrixEqn(ntaylor,scale,blc,trc);
     397             :             /* Choose a component from the list of solutions. Record the location for each scale.*/
     398             :             // Calculate penalty function
     399             :             // Find max across all pixels.
     400           0 :             chooseComponent(ntaylor, scale,criterion,blc,trc);
     401             :           }
     402             :        }//end pragma omp
     403             : 
     404             :        /* Find the best component over all scales */
     405           0 :        for(Int scale=0;scale<nscales_p;scale++)
     406             :        {
     407           0 :           if((maxScaleVal_p[scale]*scaleBias_p[scale]) > globalmaxval_p)
     408             :           {
     409           0 :             globalmaxval_p = maxScaleVal_p[scale]*scaleBias_p[scale];
     410           0 :             globalmaxpos_p = maxScalePos_p[scale];
     411           0 :             maxscaleindex_p = scale;
     412             :           }
     413             :        }
     414             : 
     415             :        /* Update the image and psf patch sizes according to the 
     416             :           current globalmaxval and psfsupport.
     417             :           This patch is over which the next-iteration's solveMatrixEqn is computed */
     418           0 :        buildImagePatches();
     419             : 
     420             : 
     421             :         /* Update the current solution by this chosen step */
     422           0 :         updateModelAndRHS(loopgain);
     423             :       
     424             :         /* Increment iteration count */
     425           0 :         totalIters_p++;
     426           0 :         iterdone++;
     427             :   
     428             :         /* Check for convergence */
     429           0 :         convergedflag = checkConvergence(criterion,fluxlimit,loopgain);
     430             : 
     431             :         /* Reached stopping threshold for this major cycle */
     432             :         /* Break out of minor-cycle loop, but signal that major cycles should continue */
     433           0 :        if(convergedflag)
     434             :         {
     435           0 :           break;
     436             :         }
     437             : 
     438           0 :   }//end minor cycle
     439           0 :   if(convergedflag==0)
     440             :     {
     441           0 :       os << "Reached max number of iterations for this minor cycle " << LogIO::POST;
     442             :     }
     443             : 
     444             :   /* returning because of non-convergence */
     445           0 :   if(convergedflag == -1)
     446             :     {
     447           0 :       os << "Stopping minor cycle iterations because of possible divergence." << LogIO::POST;
     448             :       //return(-1);
     449             :     }
     450             :   
     451             :   /********************** END MINOR CYCLE ITERATIONS ***********************/           
     452             :   
     453             :   /* Print out flux counts so far */
     454             :   //if(adbg)
     455             :   {
     456           0 :     os << "Total flux by scale :"; 
     457           0 :     for(Int scale=0;scale<nscales_p;scale++) os << "  [" << scaleSizes_p[scale] << "]: " << totalScaleFlux_p[scale] ;
     458           0 :      os << " (in this run) " << LogIO::POST;
     459           0 :      os << "Total flux by Taylor coefficient :";
     460           0 :      for(Int taylor=0;taylor<ntaylor_p;taylor++) os << "  [" << taylor << "]: " << totalTaylorFlux_p[taylor] ;
     461           0 :      os << LogIO::POST;
     462             :      //     for(Int taylor=0;taylor<ntaylor_p;taylor++) os << "Taylor " << taylor << " has total flux = " << totalTaylorFlux_p[taylor] << LogIO::POST;
     463             :   }
     464             : 
     465             :   /* Write out the multiscale model images - to disk */ 
     466             :   //  for(Int scale=0;scale<nscales_p;scale++)
     467             :   //  writeMatrixToDisk("modelimage_scale_"+String::toString(scale) , vecScaleModel_p[scale]);
     468             : 
     469             : 
     470             :   // At this stage, vecDirty_p contains the original residual images from the start of minorcycle iterations,
     471             :   // Here, we need to update vecDirty_p with image-domain residuals that account for the model 
     472             :   // components that have been found in this call to mtclean(). 
     473             :   //  Math :  Update residual :  Ires_t = Ires_t -  sum_q [ Ipsf_tq * Im_q ]
     474             :   //              where  Ipsf_tq = sum_nu ( w^{t+q} Ipsf_nu   =   Ipsf_{00, 11, 12/21}
     475           0 :   for(Int taylor1=0;taylor1<ntaylor_p;taylor1++)
     476             :     {
     477           0 :       for(Int taylor2=0;taylor2<ntaylor_p;taylor2++)
     478             :         {
     479           0 :           Matrix<Float>newModel;
     480           0 :           newModel = vecModel_p[taylor2] - vecInitialModel_p[taylor2];
     481             :           // Convolve model with psf
     482           0 :           Matrix<Complex> modelFT;  
     483           0 :           fftcomplex.fft0(modelFT, newModel, false);
     484           0 :           modelFT= ( vecPsfFT_p[taylor1+taylor2]  ) * modelFT;
     485           0 :           Matrix<Float> smoothMod(vecDirty_p[taylor1].shape());
     486           0 :           fftcomplex.fft0(smoothMod, modelFT, false);
     487           0 :           fftcomplex.flip(smoothMod, false, false);
     488             :           // Subtract from initial residual (vecDirty_p)
     489           0 :           vecDirty_p[taylor1] =  vecDirty_p[taylor1]  - smoothMod;
     490           0 :         }
     491             :     }
     492             :   // Note : In the restore step, vecDirty_p is over-written and reused in computeprincipalsolution.
     493             : 
     494             :   /* Return the number of minor-cycle iterations completed in this run */
     495           0 :   return iterdone;
     496           0 : }
     497             : 
     498             : /* Indexing Rules... */
     499           0 : Int MultiTermMatrixCleaner::IND2(Int taylor, Int scale)
     500             : {
     501           0 :         return  taylor * nscales_p + scale;
     502             : }
     503             : 
     504           0 : Int MultiTermMatrixCleaner::IND4(Int taylor1, Int taylor2, Int scale1, Int scale2)
     505             : {
     506           0 :         Int tt1=taylor1;
     507           0 :         Int tt2=taylor2;
     508           0 :         Int ts1=scale1;
     509           0 :         Int ts2=scale2;
     510           0 :         scale1 = MAX(ts1,ts2);
     511           0 :         scale2 = MIN(ts1,ts2);
     512           0 :         taylor1 = MAX(tt1,tt2);
     513           0 :         taylor2 = MIN(tt1,tt2);
     514           0 :         Int totscale = nscales_p*(nscales_p+1)/2;
     515           0 :         return ((taylor1*(taylor1+1)/2)+taylor2)*totscale + ((scale1*(scale1+1)/2)+scale2);
     516             : }
     517             : 
     518             : 
     519             :  /* Check if scale sizes are appropriate to the image size 
     520             :      If some scales are too big, ignore them.
     521             :      Reset nscales_p BEFORE all arrays are allocated */
     522           0 : Int MultiTermMatrixCleaner::verifyScaleSizes()
     523             : {
     524           0 :   LogIO os(LogOrigin("MultiTermMatrixCleaner", "verifyScaleSizes()", WHERE));
     525           0 :   Vector<Int> scaleflags(nscales_p); scaleflags=0;
     526           0 :   for(Int scale=0;scale<nscales_p;scale++)
     527             :    {
     528           0 :      if(scaleSizes_p[scale] > nx_p/2  || scaleSizes_p[scale] > ny_p/2)
     529             :        {
     530           0 :          scaleflags[scale]=1;
     531           0 :          os << LogIO::WARN << "Scale size of " << scaleSizes_p[scale] << " pixels is too big for an image size of " << nx_p << " x " << ny_p << " pixels. This scale size will be ignored.  " << LogIO::POST;
     532             :        }
     533             :    }
     534           0 :   if(sum(scaleflags)>0) // prune the scalevector and change nscales_p
     535             :     {
     536           0 :       Vector<Float> newscalesizes(nscales_p - sum(scaleflags));
     537           0 :       uInt i=0;
     538           0 :       for(Int scale=0;scale<nscales_p;scale++)
     539             :         {
     540           0 :           if(!scaleflags[scale])
     541             :           {
     542           0 :              AlwaysAssert(i<newscalesizes.nelements(),AipsError);
     543           0 :              newscalesizes[i]=scaleSizes_p[scale];
     544           0 :              i++;
     545             :           }
     546             :         }
     547           0 :           AlwaysAssert(i==newscalesizes.nelements(), AipsError);
     548           0 :           scaleSizes_p.assign(newscalesizes);
     549           0 :           nscales_p = newscalesizes.nelements();
     550           0 :           totalScaleFlux_p.resize(nscales_p,true);
     551           0 :           maxScaleVal_p.resize(nscales_p,true);
     552           0 :           maxScalePos_p.resize(nscales_p,true);
     553             : 
     554           0 :     }
     555           0 :   os << "Scale sizes to be used for deconvolution : " << scaleSizes_p << LogIO::POST;
     556             : 
     557           0 :   return scaleSizes_p.shape()(0);
     558           0 : }
     559             : 
     560             : /*************************************
     561             :  *          Allocate Memory
     562             :  *************************************/
     563             : 
     564           0 : Int MultiTermMatrixCleaner::allocateMemory()
     565             : {
     566           0 :        LogIO os(LogOrigin("MultiTermMatrixCleaner", "allocateMemory()", WHERE));
     567             :         // Define max memory usage. (half of available);
     568             :         
     569           0 :         Int ntotal4d = (nscales_p*(nscales_p+1)/2) * (ntaylor_p*(ntaylor_p+1)/2);
     570             : 
     571           0 :         Int nHess = ntotal4d;                    // cubeA
     572           0 :         Int ntempfull =  1                          // scratch
     573           0 :                                + nscales_p * 3      // vecScales, vecScaleMasks, vecWork_p, ////vecScaleModel_p
     574           0 :                                + ntaylor_p   // vecModel (vecDirty is a ref)
     575           0 :                                + 2 * nscales_p * ntaylor_p; // matR and matCoeff
     576             : 
     577           0 :         Int ntemphalf = 2  // scratch
     578           0 :                                 + nscales_p  // vecSCalesFT
     579           0 :                                 + psfntaylor_p; // vecPsfFT
     580             : 
     581           0 :         Int numMB = static_cast<Int>(( nx_p*ny_p*4*(ntempfull + ntemphalf/2.0)  + psfsupport_p[0]*psfsupport_p[1]*nHess  )/(1024*1024));
     582           0 :         memoryMB_p = Double(HostInfo::memoryTotal()/1024);
     583             : 
     584           0 :         if(adbg) os << "This algorithm needs to allocate " << numMB << " MBytes." << LogIO::POST;
     585           0 :         if (numMB > 0.75*memoryMB_p) 
     586             :        { 
     587           0 :          os << LogIO::WARN << "This algorithm needs to allocate " << numMB << " MBytes for " << ntempfull + ntemphalf << " Matrices. " << LogIO::POST;
     588           0 :            os << LogIO::WARN << "Available memory for this process is  " << memoryMB_p << " MB. Please reduce imsize/nscales/nterms and try again " << LogIO::POST;
     589           0 :            return -1;
     590             :        }
     591             : 
     592             :         // shape of all matrices
     593           0 :         gip = IPosition(2,nx_p,ny_p);  
     594           0 :         IPosition tgip(2,ntaylor_p,ntaylor_p);
     595             :         
     596             :         // Small HessianPeak matrix to be inverted for each point..
     597           0 :         matA_p.resize(nscales_p); invMatA_p.resize(nscales_p);
     598           0 :         for(Int i=0;i<nscales_p;i++)
     599             :         {
     600           0 :           matA_p[i].resize(tgip); 
     601           0 :           invMatA_p[i].resize(tgip);
     602             :         }
     603             :         
     604             :         // I_D 
     605           0 :         dirtyFT_p.resize(); 
     606             :         
     607             :         // Temporary work-holder
     608           0 :         cWork_p.resize(); 
     609             :         //tWork_p.resize(gip);
     610             :         
     611             :         // Scales 
     612           0 :         vecScales_p.resize(nscales_p);
     613           0 :         vecScalesFT_p.resize(nscales_p);
     614           0 :         vecScaleMasks_p.resize(nscales_p);
     615           0 :         vecWork_p.resize(nscales_p);
     616             :         //        vecScaleModel_p.resize(nscales_p);
     617           0 :         for(Int i=0;i<nscales_p;i++) 
     618             :         {
     619             :           //vecScales_p[i].resize(gip); // OLD
     620           0 :           vecScales_p[i].resize(psfsupport_p)
     621             : ;
     622           0 :           vecScalesFT_p[i].resize();
     623           0 :           vecScaleMasks_p[i].resize(gip);
     624           0 :           vecWork_p[i].resize(gip);
     625             :           //          vecScaleModel_p[i].resize(gip);
     626             :           //          vecScaleModel_p[i] = 0.0;
     627             :         }
     628             :         
     629             :         // Psfs and Models 
     630           0 :         vecPsfFT_p.resize(psfntaylor_p);
     631           0 :         for(Int i=0;i<psfntaylor_p;i++) 
     632             :         {
     633           0 :                   vecPsfFT_p[i].resize();
     634             :         }
     635             :         
     636             :         // Dirty/Residual Images
     637           0 :         vecDirty_p.resize(ntaylor_p);
     638           0 :         vecModel_p.resize(ntaylor_p);
     639           0 :         vecInitialModel_p.resize(ntaylor_p);
     640           0 :         for(Int i=0;i<ntaylor_p;i++) 
     641             :         {
     642           0 :           vecDirty_p[i].resize();
     643           0 :           vecModel_p[i].resize(gip);
     644           0 :           vecInitialModel_p[i].resize(gip);
     645             :         }
     646             :         
     647             :         // (Psf * Scales) * (Psf * Scales)
     648           0 :         cubeA_p.resize(ntotal4d);
     649             : 
     650           0 :         for(Int i=0;i<ntotal4d;i++) 
     651             :         {
     652             :           //cubeA_p[i].resize(gip); // OLD
     653           0 :           cubeA_p[i].resize(psfsupport_p); // NEW
     654             :         }
     655             :         
     656             :         // I_D * (Psf * Scales)
     657           0 :         matR_p.resize(ntaylor_p*nscales_p);
     658             :         // Coefficients to be solved for.
     659           0 :         matCoeffs_p.resize(ntaylor_p*nscales_p);
     660             :         
     661           0 :         for(Int i=0;i<ntaylor_p*nscales_p;i++) 
     662             :         {
     663           0 :                   matR_p[i].resize(gip); 
     664           0 :                   matCoeffs_p[i].resize(gip);
     665           0 :                   matCoeffs_p[i] = 0.0;
     666             :         }
     667             :   
     668           0 :         return 0;
     669           0 : }
     670             : 
     671             : /* Make masks for each scale size */
     672             : /* At the end, force a scale-dependent border */
     673           0 : Int MultiTermMatrixCleaner::setupUserMask()
     674             : {
     675           0 :     if(itsMask.null())
     676             :     {
     677             :            /* Make a mask the full size, for all scales */
     678           0 :            for(Int scale=0;scale<nscales_p;scale++)
     679             :             {
     680           0 :               vecScaleMasks_p[scale] = 1.0;
     681             :             }
     682             :     }
     683             :     else
     684             :     {
     685             :       /* Compute the convolutions of the input mask with each scale size */
     686           0 :       if(adbg) os << "Convolving user mask with scales, and using 0.1 as a threshold." << LogIO::POST;
     687           0 :       Matrix<Complex> maskft;
     688           0 :       fftcomplex.fft0(maskft , *itsMask , false);
     689             : 
     690           0 :       for(Int scale=0;scale<nscales_p;scale++)
     691             :       {
     692           0 :         AlwaysAssert(maskft.shape() == vecScalesFT_p[scale].shape(),AipsError);
     693           0 :         cWork_p.assign(maskft * vecScalesFT_p[scale]);
     694           0 :         fftcomplex.fft0( vecScaleMasks_p[scale] , cWork_p , false  );
     695           0 :         fftcomplex.flip( vecScaleMasks_p[scale] , false , false);
     696             :         
     697           0 :         for (Int j=0 ; j < (itsMask->shape())(1); ++j)
     698           0 :         for (Int k =0 ; k < (itsMask->shape())(0); ++k)
     699             :         {
     700           0 :           if(itsMaskThreshold > (Float)0.0)  // if negative, use the continuous mask
     701             :                 {
     702           0 :                   (vecScaleMasks_p[scale])(k,j) =  (vecScaleMasks_p[scale])(k,j) > (Float)0.1 ? 1.0 : 0.0;
     703             :                 }
     704             :         }
     705             :         
     706             :         //             writeMatrixToDisk("scalemask."+String::toString(scale), vecScaleMasks_p[scale]);
     707             :       }// end of for scale
     708             : 
     709             :           /* TO DO... maybe :
     710             :              Map some local variables to those in the parant MatrixCleaner class 
     711             :              so that its functions can be used. Later, clean this up by directly using
     712             :              the MatrixCleaner variables.  Currently not possible because the 
     713             :              function that creates scales and populates itsScaleXfrs (setscales), 
     714             :              does a lot of extra stuff that we don't need here. */
     715             : 
     716             :           /* Must reuse itsNscales, itsScaleSizes, itsScales, itsScaleXfrs, itsScaleMasks, itsScalesValid 
     717             :               and eliminate : vecScales_p, scaleSizes_p, vecScalesFT_p, vecScaleMasks_p */
     718             :           /*
     719             :           itsScaleXfrs.resize(nscales_p);
     720             :           itsNscales = nscales_p;
     721             :           for(Int scale=0;scale<nscales_p;scale++)
     722             :             {
     723             :               itsScaleXfrs[scale].reference(vecScalesFT_p[scale]);
     724             :             }
     725             :           itsScalesValid=true;
     726             : 
     727             :           makeScaleMasks();
     728             :           AlwaysAssert(nscales_p==itsScaleMasks.nelements(),AipsError);
     729             : 
     730             :           for(Int scale=0;scale<nscales_p;scale++)
     731             :             {
     732             :               cout << "Shapes for scale : " << scale << itsScaleMasks[scale] << endl;
     733             :               vecScaleMasks_p[scale].reference(itsScaleMasks[scale]);
     734             :             }
     735             :           */
     736             : 
     737           0 :     }// end of if else
     738             : 
     739             : 
     740             :  /* Set the edges of the masks according to the scale size */
     741             :  // Set the values OUTSIDE the box to zero....
     742             :   // "The vecScaleMasks_p are used when finding the peakres for each iteration. If we subtract
     743             :   // the (tt/scale convolved) psf from too close to the edge of the image, we could add a hard edge
     744             :   // to the image, which would result in ringing in an fft. Therefore, we use a border to ensure
     745             :   // that we don't pick a peakres position too close to the edge." - approximately what Urvashi said
     746             :   // Start at 1: don't add a 1-pixel border to vecScaleMasks_p[0], because that mask is used to find
     747             :   // the peakres in checkConvergence(), which could result in the minor cycle disagreeing with the
     748             :   // major cycle about what the peakres value is if the peakres is at the edge.
     749           0 :   for(Int scale=1;scale<nscales_p;scale++)
     750             :   {
     751           0 :       Int border = (Int)(scaleSizes_p[scale]*1.5);
     752             :       // bottom
     753           0 :       IPosition blc1(2, 0 , 0 );
     754           0 :       IPosition trc1(2,nx_p-1, border );
     755           0 :       IPosition inc1(2, 1);
     756           0 :       LCBox::verify(blc1,trc1,inc1,vecScaleMasks_p[scale].shape());
     757           0 :       (vecScaleMasks_p[scale])(blc1,trc1) = 0.0;
     758             :       // top
     759           0 :       blc1[0]=0; blc1[1]=ny_p-border-1;
     760           0 :       trc1[0]=nx_p-1; trc1[1]=ny_p-1;
     761           0 :       LCBox::verify(blc1,trc1,inc1,vecScaleMasks_p[scale].shape());
     762           0 :       (vecScaleMasks_p[scale])(blc1,trc1) = 0.0;
     763             :       // left
     764           0 :       blc1[0]=0; blc1[1]=border;
     765           0 :       trc1[0]=border; trc1[1]=ny_p-border-1;
     766           0 :       LCBox::verify(blc1,trc1,inc1,vecScaleMasks_p[scale].shape());
     767           0 :       (vecScaleMasks_p[scale])(blc1,trc1) = 0.0;
     768             :       // right
     769           0 :       blc1[0]=nx_p-border-1; blc1[1]=border;
     770           0 :       trc1[0]=nx_p; trc1[1]=ny_p-border-1;
     771           0 :       LCBox::verify(blc1,trc1,inc1,vecScaleMasks_p[scale].shape());
     772           0 :       (vecScaleMasks_p[scale])(blc1,trc1) = 0.0;
     773             : 
     774             :       //      writeMatrixToDisk("maskforscale_"+String::toString(scale),vecScaleMasks_p[scale]);
     775           0 :   }
     776             : 
     777           0 :    return 0;
     778             : }/* end of setupUserMask() */
     779             : 
     780             : 
     781             : /***************************************
     782             :  *  Set up the Blobs of various scales.
     783             :  ****************************************/
     784             : 
     785           0 : Int MultiTermMatrixCleaner::setupScaleFunctions()
     786             : {
     787           0 :   LogIO os(LogOrigin("MultiTermMatrixCleaner", "setupScaleFunctions", WHERE));
     788             :         // Set the scale sizes
     789           0 :         if(scaleSizes_p.nelements()==0)
     790             :         {
     791           0 :                 scaleSizes_p.resize(nscales_p);
     792           0 :                 Float scaleInc = 2.0;
     793           0 :                 scaleSizes_p[0] = 0.0;
     794           0 :                 for (Int scale=1; scale<nscales_p;scale++) 
     795             :                 {
     796           0 :                         scaleSizes_p[scale] = scaleInc * pow(10.0, (Float(scale)-2.0)/2.0) ;
     797             :                 }  
     798             :         }
     799             : 
     800             :         /* SCALE BIAS : Can get rid of this entirely. Leaving it in, just in case scalebiases are needed later */       
     801           0 :         scaleBias_p.resize(nscales_p);
     802           0 :         totalScaleFlux_p.resize(nscales_p);
     803             :         //Float prefScale=2.0;
     804             :         //Float fac=6.0;
     805           0 :         if(nscales_p>1)
     806             :         {
     807           0 :                 for(Int scale=0;scale<nscales_p;scale++) 
     808             :                 {
     809           0 :                   scaleBias_p[scale] = 1 - itsSmallScaleBias * scaleSizes_p[scale]/scaleSizes_p(nscales_p-1);
     810             :                         //scaleBias_p[scale] = 1 - 0.4 * scaleSizes_p[scale]/scaleSizes_p(nscales_p-1);
     811             :                   //scaleBias_p[scale] = 1.0;
     812             :                         //////scaleBias_p[scale] = pow((Float)scale/fac,prefScale)*exp(-1.0*scale/fac)/(pow(prefScale/fac,prefScale)*exp(-1.0*prefScale/fac));
     813             :                         //scaleBias_p[scale] = pow((Float)(scale+1)/fac,prefScale)*exp(-1.0*(scale+1)/fac);
     814           0 :                   os << "scale " << scale+1 << " = " << scaleSizes_p(scale) << " pixels with bias = " << scaleBias_p[scale] << LogIO::POST;
     815           0 :                   totalScaleFlux_p[scale]=0.0;
     816             :                 }
     817             :         }
     818           0 :         else scaleBias_p[0]=1.0;
     819             :  
     820             :         // Compute the scaled blobs - prolate spheroids with tapering and truncation
     821             :         // vecScales_p, scaleSizes_p, vecScalesFT_p
     822           0 :         if(!donePSP_p) 
     823             :         {
     824             :                 // Compute h(s1), h(s2),... depending on the number of scales chosen.
     825             :                 // NSCALES = 1;
     826           0 :                 if(adbg) os << "Calculating scales and their FTs " << LogIO::POST;
     827             :                         
     828           0 :                 for (Int scale=0; scale<nscales_p;scale++) 
     829             :                 {
     830             :                   /*
     831             :                         // First make the scale
     832             :                         makeScale(vecScales_p[scale], scaleSizes_p(scale));
     833             :                         // Now store the XFR (shape of FT is set automatically)
     834             :                         fftcomplex.fft0(vecScalesFT_p[scale] , vecScales_p[scale] , false);
     835             :                         //fftcomplex.flip(vecScalesFT_p[scale] , false , false);
     836             :                   */
     837             :                         // First make the scale
     838           0 :                         makeScale(vecWork_p[0], scaleSizes_p(scale));
     839             :                         // Now store the XFR (shape of FT is set automatically)
     840           0 :                         fftcomplex.fft0(vecScalesFT_p[scale] , vecWork_p[0] , false);
     841             :                         // Copy the scale onto the smaller vecScales image.
     842           0 :                         IPosition immid(2,nx_p/2, ny_p/2);
     843           0 :                         Matrix<Float> psfpatch = ( vecWork_p[0] ) (immid-psfsupport_p/2,immid+psfsupport_p/2-IPosition(2,1,1));  
     844           0 :                         vecScales_p[scale] = psfpatch; 
     845             : 
     846           0 :                 }
     847           0 :                 donePSP_p=true;
     848             :         }
     849             :         
     850           0 :         return 0;
     851           0 : }/* end of setupBlobs() */
     852             : 
     853             : 
     854             : 
     855             : /***************************************
     856             :  *  Compute convolutions and the A matrix.
     857             :  ****************************************/
     858             : 
     859           0 : Int MultiTermMatrixCleaner::computeHessianPeak()
     860             : {
     861           0 :   LogIO os(LogOrigin("MultiTermMatrixCleaner", "computeHessianPeak", WHERE));
     862           0 :    gip = IPosition(2,nx_p,ny_p);  
     863             :    
     864           0 :    if(!doneCONV_p)
     865             :    {  
     866             :       // Compute the convolutions of the smoothed psfs with each other.
     867             :       // Compute Assxx
     868             :       // Compute A100, A101, A102
     869             :       //         A110, A111, A112
     870             :       //         A120, A121, A122   for h(s1)
     871             :       // Compute A200, A201, A202
     872             :       //         A210, A211, A212
     873             :       //         A220, A221, A222   for h(s2)
     874             :       //... depending on the number of scales chosen
     875             : 
     876             :       // (PSF * scale) * (PSF * scale) -> cubeA_p [nx_p,ny_p,ntaylor,ntaylor,nscales]
     877           0 :       os << "Calculating PSF and Scale convolutions " << LogIO::POST;
     878           0 :       for (Int taylor1=0; taylor1<ntaylor_p;taylor1++) 
     879           0 :       for (Int taylor2=0; taylor2<=taylor1;taylor2++) 
     880           0 :       for (Int scale1=0; scale1<nscales_p;scale1++) 
     881           0 :       for (Int scale2=0; scale2<=scale1;scale2++) 
     882             :       {
     883           0 :         Int ttay1 = taylor1+taylor2;
     884             :         
     885             :         // CALC Hess : Calculate  PSF_(t1+t2)  * scale_1 * scale 2
     886             : 
     887           0 :         cWork_p.assign( (vecPsfFT_p[ttay1]) *(vecScalesFT_p[scale1])*(vecScalesFT_p[scale2]) );
     888             : 
     889           0 :         Bool usepatch=true;
     890           0 :         if(usepatch)
     891             :           {
     892           0 :             fftcomplex.fft0( vecWork_p[0]  , cWork_p , false  );
     893           0 :             Matrix<Float> psfpatch = ( vecWork_p[0] ) (itsPositionPeakPsf-psfsupport_p/2,itsPositionPeakPsf+psfsupport_p/2-IPosition(2,1,1));  
     894           0 :             cubeA_p[IND4(taylor1,taylor2,scale1,scale2)] = psfpatch; 
     895           0 :           }
     896             :         else
     897             :           {
     898           0 :             fftcomplex.fft0( cubeA_p[IND4(taylor1,taylor2,scale1,scale2)]  , cWork_p , false  );
     899             :           }
     900             :         
     901             :         //writeMatrixToDisk("psfconv_t_"+String::toString(taylor1)+"-"+String::toString(taylor2)+"_s_"+String::toString(scale1)+"-"+String::toString(scale2)+".im", cubeA_p[IND4(taylor1,taylor2,scale1,scale2)] );
     902             :       }   
     903             : 
     904             :       // Construct A, invA for each scale.
     905           0 :       if(itsPositionPeakPsf != IPosition(2,(nx_p/2),(ny_p/2)))
     906           0 :         os << LogIO::WARN << "The PSF peak is not at the center of the image..." << LogIO::POST;
     907             : 
     908           0 :       Int stopnow=false;
     909           0 :       for (Int scale=0; scale<nscales_p;scale++) 
     910             :       {
     911             :               // Fill up A
     912           0 :               for (Int taylor1=0; taylor1<ntaylor_p;taylor1++) 
     913           0 :               for (Int taylor2=0; taylor2<ntaylor_p;taylor2++) 
     914             :               {
     915             :                 //(matA_p[scale])(taylor1,taylor2) = (cubeA_p[IND4(taylor1,taylor2,scale,scale)])(itsPositionPeakPsf); // OLD
     916           0 :                 (matA_p[scale])(taylor1,taylor2) = (cubeA_p[IND4(taylor1,taylor2,scale,scale)])(psfpeak_p); // NEW
     917             :                 /* Check for exact zeros ON MAIN DIAGONAL. Usually indicative of error */
     918           0 :                 if( taylor1==taylor2 &&
     919           0 :                     fabs( (matA_p[scale])(taylor1,taylor2) )  == 0.0 ) stopnow = true;
     920             :               }
     921             :               
     922           0 :               if(stopnow)
     923             :               {
     924           0 :                 os << "Multi-Term Hessian has exact zeros on its main diagonal. Not proceeding further." << LogIO::WARN << endl;
     925           0 :                 os << "The Matrix [A] is : " << (matA_p[scale]) << LogIO::POST;
     926           0 :                 return -2;
     927             :               }
     928             : 
     929             :               /* If all elements are non-zero, check by brute-force if the rows/cols 
     930             :                  are nearly linearly dependant, making the matrix nearly singular. */
     931           0 :                Vector<Float> ratios(ntaylor_p);
     932           0 :                Float tsum=0.0;
     933           0 :                for(Int taylor1=0; taylor1<ntaylor_p-1; taylor1++)
     934             :                {
     935           0 :                    for(Int taylor2=0; taylor2<ntaylor_p; taylor2++)
     936           0 :                      ratios[taylor2] = (matA_p[scale])(taylor1,taylor2) / (matA_p[scale])(taylor1+1,taylor2);
     937           0 :                    tsum=0.0;
     938           0 :                    for(Int taylor2=0; taylor2<ntaylor_p-1; taylor2++)
     939           0 :                       tsum += fabs(ratios[taylor2] - ratios[taylor2+1]);
     940           0 :                    tsum /= ntaylor_p-1;
     941           0 :                    if(tsum < 1e-04) stopnow=true;
     942             : 
     943             :                    //cout << "Ratios for row " << taylor1 << " are " << ratios << endl;
     944             :                    //cout << "tsum : " << tsum << endl;
     945             :                }
     946             : 
     947           0 :               if(stopnow)
     948             :               {
     949           0 :                 os << "Multi-Term Hessian has linearly-dependent rows/cols. Not proceeding further." << LogIO::WARN << endl;
     950           0 :                 os << "The Matrix [A] is : " << (matA_p[scale]) << LogIO::POST;
     951           0 :                 return -2;
     952             :               }
     953             :                   
     954             :               /* Try to invert the matrix. If it fails, return. 
     955             :                  The invertSymPosDef will check if it's positive definite or not. 
     956             :                  By construction, it should be pos-def. */
     957             :               // Compute inv(A) 
     958             :               // Use MatrixMathLA::invert
     959             :               // or Use invertSymPosDef...
     960             :               //
     961             :               try
     962             :               {
     963           0 :                      Double deter=0.0;
     964             :                      //invert((*invMatA_p[scale]),deter,(*matA_p[scale]));
     965             :                      //os << "Matrix Inverse : inv(A) : " << (*invMatA_p[scale]) << LogIO::POST;
     966             :               
     967             :                      //if(adbg) 
     968           0 :                      os << "The Matrix [H] for " << scaleSizes_p[scale] << " pixel scale is : " << (matA_p[scale]) << LogIO::POST;
     969           0 :                      invertSymPosDef((invMatA_p[scale]),deter,(matA_p[scale]));
     970           0 :                      if(adbg) os << "Lapack Cholesky Decomposition Inverse of [A] is : " << (invMatA_p[scale]) << LogIO::POST;
     971             :                      //if(adbg)os << "A matrix determinant : " << deter << LogIO::POST;
     972             :                      //if(fabs(deter) < 1.0e-08) os << "SINGULAR MATRIX !! STOP!! " << LogIO::EXCEPTION;
     973             :               }
     974           0 :               catch(AipsError &x)
     975             :               {
     976           0 :                       os << "The Matrix [A] is : " << (matA_p[scale]) << LogIO::POST;
     977           0 :                       os << "Cannot Invert matrix : " << x.getMesg() << LogIO::WARN;
     978           0 :                       return -2;
     979           0 :               }
     980             : 
     981             :               /* Add more checks based on the condition number of these matrices */
     982             :               
     983           0 :       }// end for scales
     984             :       
     985           0 :       doneCONV_p=true;
     986             :    } 
     987             :    
     988           0 :    return 0;
     989           0 : }/* end of computeHessianPeak() */
     990             : 
     991             : 
     992             : /***************************************
     993             :  *  Compute convolutions of the residual images ( all specs ) with scales.
     994             :  *  --> the Right-Hand-Side of the matrix equation.
     995             :  ****************************************/
     996             : 
     997           0 : Int MultiTermMatrixCleaner::computeRHS()
     998             : {
     999             :   //  LogIO os(LogOrigin("MultiTermMatrixCleaner", "computeRHS()", WHERE));
    1000           0 :         IPosition blc1(2,0,0);
    1001           0 :         IPosition trc1(2,nx_p,ny_p);
    1002           0 :         IPosition inc1(2, 1);
    1003             :         
    1004             :         /* Compute R10 = I_D*B10, R11 = I_D*B11, R12 = I_D*B12
    1005             :          * Compute R20 = I_D*B20, R21 = I_D*B21, R22 = I_D*B22
    1006             :          * ... depending on the number of scales chosen.
    1007             :          */
    1008             : 
    1009             :         /* I_D * (PSF * scale) -> matR_p [nx_p,ny_p,ntaylor,nscales] */
    1010           0 :         for (Int taylor=0; taylor<ntaylor_p;taylor++) 
    1011             :         {
    1012             :            /* Compute FT of dirty image */
    1013           0 :           fftcomplex.fft0( dirtyFT_p , vecDirty_p[taylor] , false );
    1014             :            
    1015           0 :            for (Int scale=0; scale<nscales_p;scale++) 
    1016             :            {
    1017             : 
    1018             :              // CALC RHS :  Calculate   Dirty_t  * scale_s
    1019             : 
    1020             :                 // Let cWork get resized if needed. Force matR to have the right shape
    1021             :                 ////                cWork_p.assign( (dirtyFT_p)*(vecPsfFT_p[0])*(vecScalesFT_p[scale]) );
    1022           0 :                 cWork_p.assign( (dirtyFT_p)*(vecScalesFT_p[scale]) );
    1023           0 :                 fftcomplex.fft0( matR_p[IND2(taylor,scale)] , cWork_p , false );
    1024           0 :                 fftcomplex.flip(  matR_p[IND2(taylor,scale)] , false , false );
    1025             :            }
    1026             :            //      writeMatrixToDisk("resid_"+String::toString(taylor)+".im", matR_p[IND2(taylor,0)] );
    1027             :         }
    1028             :         
    1029           0 :         return 0;
    1030           0 : }/* end of computeRHS() */
    1031             : 
    1032             : /***************************************
    1033             :  *  Compute  flux limit for minor cycles
    1034             :  ****************************************/
    1035             : 
    1036           0 : Int MultiTermMatrixCleaner::computeFluxLimit(Float &fluxlimit, Float threshold)
    1037             : {
    1038             :   // LogIO os(LogOrigin("MultiTermMatrixCleaner", "computeFluxLimit", WHERE));
    1039             : 
    1040             :   // For now, since this calculation is being done outside.....
    1041           0 :     fluxlimit = threshold;
    1042             : 
    1043             :     /// fluxlimit = MAX(threshold , itsGain * max....
    1044             : 
    1045             :   // Later, implement the equivalent of an iteration-based cycle-speedup - if required.
    1046             : 
    1047           0 :         return 0;
    1048             : }/* end of computeFluxLimit() */
    1049             : 
    1050             : 
    1051             : 
    1052             : /***************************************
    1053             :  *  Solve the matrix eqn for each point in the lattice.
    1054             : Note : This function is called within the 'scale' omp/pragma loop. Needs to be thread-safe
    1055             :  ****************************************/
    1056           0 : Int MultiTermMatrixCleaner::solveMatrixEqn(Int ntaylor, Int scale, IPosition blc, IPosition trc)
    1057             : {
    1058             :   
    1059           0 :         for(Int taylor1=0;taylor1<ntaylor;taylor1++)
    1060             :         {
    1061           0 :              Matrix<Float> coeffs = (matCoeffs_p[IND2(taylor1,scale)])(blc,trc);  
    1062           0 :              coeffs = 0.0;
    1063           0 :              for(Int taylor2=0;taylor2<ntaylor;taylor2++)
    1064             :              {
    1065           0 :                Matrix<Float> rhs = (matR_p[IND2(taylor2,scale)])(blc,trc);
    1066           0 :                coeffs = coeffs +  ((Float)(invMatA_p[scale])(taylor1,taylor2))* rhs;
    1067           0 :              }
    1068           0 :         }
    1069             :   
    1070             :   /* Solve for the coefficients, one scale at at time*/
    1071             :         /*      
    1072             :         for(Int taylor1=0;taylor1<ntaylor;taylor1++)
    1073             :         {
    1074             :              (matCoeffs_p[IND2(taylor1,scale)]) = 0.0; 
    1075             :              for(Int taylor2=0;taylor2<ntaylor;taylor2++)
    1076             :              {
    1077             :                   matCoeffs_p[IND2(taylor1,scale)] = matCoeffs_p[IND2(taylor1,scale)] + ((Float)(invMatA_p[scale])(taylor1,taylor2))*(matR_p[IND2(taylor2,scale)]);
    1078             :              }
    1079             :         }
    1080             :         */
    1081           0 :         return 0;
    1082             : }/* end of solveMatrixEqn() */
    1083             :         
    1084             : /***************************************
    1085             :  * Find the 'peak' and its location. Fill this into maxScaleVal_p, maxScalePos_p  
    1086             :  *  Compute the update direction
    1087             :  *  Options are
    1088             :  *  (1) Original implementation : something related to the derivative of Chi-Sq. Not sure why it ever worked.
    1089             :  *  (2) Peak residual in the 00 image
    1090             :  *  (3) Peak t_0 component (for each scale)
    1091             :  *  (4) Derivative of chi-square
    1092             :  *  (5) Similar to (1), but with the PSF replaced by the peak of the PSF only.
    1093             : Note : This function is called within the 'scale' omp/pragma loop. Needs to be thread-safe
    1094             :  ****************************************/
    1095             : 
    1096           0 :   Int MultiTermMatrixCleaner::chooseComponent(Int ntaylor, Int scale, Int /*criterion*/, IPosition blc, IPosition trc)
    1097             : {
    1098             : 
    1099             : 
    1100           0 :   Matrix<Float> work = ( vecWork_p[scale] )(blc,trc);  
    1101             : 
    1102           0 :   work = 0.0;
    1103           0 :   for(Int taylor1=0;taylor1<ntaylor;taylor1++)
    1104             :     {
    1105           0 :       Matrix<Float> coeffs1 = (matCoeffs_p[IND2(taylor1,scale)])(blc,trc);
    1106           0 :       Matrix<Float> resid = (matR_p[IND2(taylor1,scale)])(blc,trc);
    1107             : 
    1108             :       /*
    1109             :       work = work + (Float)2.0 * coeffs1 * resid;
    1110             :       for(Int taylor2=0;taylor2<ntaylor;taylor2++)
    1111             :         {
    1112             :           Matrix<Float> coeffs2 = (matCoeffs_p[IND2(taylor2,scale)])(blc,trc);
    1113             :           work = work - (Float)((matA_p[scale])(taylor1,taylor2)) * coeffs1 * coeffs2;
    1114             :         }
    1115             :       */
    1116             : 
    1117           0 :       work = work +  coeffs1 * resid;
    1118             : 
    1119           0 :     }
    1120           0 :   findMaxAbsMask(vecWork_p[scale],vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
    1121             :   
    1122             :   /*
    1123             :     
    1124             :   vecWork_p[scale] = 0.0;
    1125             :   for(Int taylor1=0;taylor1<ntaylor;taylor1++)
    1126             :     {
    1127             :       vecWork_p[scale] = vecWork_p[scale] + (Float)2.0  * (  (matCoeffs_p[IND2(taylor1,scale)])  *  (matR_p[IND2(taylor1,scale)])  );
    1128             :       for(Int taylor2=0;taylor2<ntaylor;taylor2++)
    1129             :         vecWork_p[scale] = vecWork_p[scale] - (Float)((matA_p[scale])(taylor1,taylor2)) * matCoeffs_p[IND2(taylor1,scale)] * matCoeffs_p[IND2(taylor2,scale)] ;
    1130             :     }
    1131             :   findMaxAbsMask(vecWork_p[scale],vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
    1132             :   */
    1133             : 
    1134             : 
    1135             :   //    switch(criterion)
    1136             :   //    {
    1137             :        //     case 1 : /* For each scale, find the maximum chi-square derivative (maybe) */
    1138             :          //       {
    1139             :          /*
    1140             :                /// Code block using a private matrix
    1141             :                Matrix<Float> ttWork_p((matR_p[IND2(0,0)]).shape());
    1142             :                ttWork_p = 0.0;
    1143             :                 for(Int taylor1=0;taylor1<ntaylor_p;taylor1++)
    1144             :                 {
    1145             :                       ttWork_p = ttWork_p + (Float)2.0  * (  (matCoeffs_p[IND2(taylor1,scale)])  *  (matR_p[IND2(taylor1,scale)])  );
    1146             :                       for(Int taylor2=0;taylor2<ntaylor_p;taylor2++)
    1147             :                              ttWork_p = ttWork_p -  matCoeffs_p[IND2(taylor1,scale)] * matCoeffs_p[IND2(taylor2,scale)] * cubeA_p[IND4(taylor1,taylor2,scale,scale)];
    1148             :                 }
    1149             :                 findMaxAbsMask(ttWork_p,vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
    1150             :          */
    1151             :          /*
    1152             :                /// Code block using pre-allocated matrices
    1153             :                vecWork_p[scale] = 0.0;
    1154             :                 for(Int taylor1=0;taylor1<ntaylor_p;taylor1++)
    1155             :                 {
    1156             :                       vecWork_p[scale] = vecWork_p[scale] + (Float)2.0  * (  (matCoeffs_p[IND2(taylor1,scale)])  *  (matR_p[IND2(taylor1,scale)])  );
    1157             :                       for(Int taylor2=0;taylor2<ntaylor_p;taylor2++)
    1158             :                              vecWork_p[scale] = vecWork_p[scale] -  matCoeffs_p[IND2(taylor1,scale)] * matCoeffs_p[IND2(taylor2,scale)] * cubeA_p[IND4(taylor1,taylor2,scale,scale)];
    1159             :                 }
    1160             :                 findMaxAbsMask(vecWork_p[scale],vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
    1161             :          */
    1162             :        //       }
    1163             :        //       break;
    1164             :        /*
    1165             :      case 2 : // For each scale, find the peak residual 
    1166             :        {
    1167             :                Float norm = sqrt((1.0/(matA_p[scale])(0,0)));
    1168             :                findMaxAbsMask( norm*matR_p[IND2(0,scale)] , vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
    1169             :        }
    1170             :        break;
    1171             :      case 3: // For each scale, find the peak t_0 coefficnent 
    1172             :        {
    1173             :                findMaxAbsMask( matCoeffs_p[IND2(0,scale)] , vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
    1174             :        }
    1175             :        break;
    1176             :      case 4: // For each scale find the max chisq derivative (diag approx for all hessians) - BAD
    1177             :        {
    1178             :                Matrix<Float> ttWork_p((matR_p[IND2(0,0)]).shape());
    1179             :                 ttWork_p = 0.0;
    1180             :                 for(Int taylor1=0;taylor1<ntaylor_p;taylor1++)
    1181             :                 {
    1182             :                     ttWork_p = ttWork_p + (Float)(2.0  * (matA_p[scale])(taylor1,taylor1) ) *  (matR_p[IND2(taylor1,scale)])  ;
    1183             :                     for(Int taylor2=0;taylor2<ntaylor_p;taylor2++)
    1184             :                         ttWork_p = ttWork_p - (Float)(2.0*( (matA_p[scale])(taylor1,taylor1) * (matA_p[scale])(taylor1,taylor2))) * matCoeffs_p[IND2(taylor2,scale)] ; 
    1185             :                 }
    1186             :                 findMaxAbsMask(ttWork_p,vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
    1187             :        }
    1188             :        break;
    1189             :        */
    1190             :   //     case 5 : /* For each scale, same as 1, but use only psf peaks */
    1191             :   //      {
    1192             :          /*    
    1193             :                /// Code block using a private matrix
    1194             :                Matrix<Float> ttWork_p((matR_p[IND2(0,0)]).shape());
    1195             :                ttWork_p = 0.0;
    1196             :                 for(Int taylor1=0;taylor1<ntaylor_p;taylor1++)
    1197             :                 {
    1198             :                       ttWork_p = ttWork_p + (Float)2.0  * (  (matCoeffs_p[IND2(taylor1,scale)])  *  (matR_p[IND2(taylor1,scale)])  );
    1199             :                       for(Int taylor2=0;taylor2<ntaylor_p;taylor2++)
    1200             :                         ttWork_p = ttWork_p - (Float)((matA_p[scale])(taylor1,taylor2)) * matCoeffs_p[IND2(taylor1,scale)] * matCoeffs_p[IND2(taylor2,scale)] ;
    1201             :                 }
    1202             :                 findMaxAbsMask(ttWork_p,vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
    1203             :          */
    1204             :   /*            
    1205             :          /// Code block using pre-allocated matrices
    1206             :                 vecWork_p[scale] = 0.0;
    1207             :                 for(Int taylor1=0;taylor1<ntaylor;taylor1++)
    1208             :                 {
    1209             :                   vecWork_p[scale] = vecWork_p[scale] + (Float)2.0  * (  (matCoeffs_p[IND2(taylor1,scale)])  *  (matR_p[IND2(taylor1,scale)])  );
    1210             :                   for(Int taylor2=0;taylor2<ntaylor;taylor2++)
    1211             :                         vecWork_p[scale] = vecWork_p[scale] - (Float)((matA_p[scale])(taylor1,taylor2)) * matCoeffs_p[IND2(taylor1,scale)] * matCoeffs_p[IND2(taylor2,scale)] ;
    1212             :                 }
    1213             :                 findMaxAbsMask(vecWork_p[scale],vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
    1214             :                 
    1215             : 
    1216             :        }
    1217             :        break;
    1218             :   */
    1219             :        /*
    1220             :      case 6 : // chi-square for this scale = sum of abs of residual images for taylor terms 
    1221             :        {
    1222             : 
    1223             :          /// Code block using pre-allocated matrices
    1224             :                 vecWork_p[scale] = 0.0;
    1225             :                 for(Int taylor1=0;taylor1<ntaylor;taylor1++)
    1226             :                 {
    1227             :                   vecWork_p[scale] = vecWork_p[scale] + (  abs(matR_p[IND2(taylor1,scale)])  );
    1228             :                 }
    1229             :                 findMaxAbsMask(vecWork_p[scale],vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
    1230             :                 
    1231             : 
    1232             :        }
    1233             :        break;
    1234             :        */
    1235             :        //     default:
    1236             :        //       os << LogIO::SEVERE << "Internal error : Unknown option for type of update direction" << LogIO::POST;
    1237             :        //     }
    1238             :  
    1239             : 
    1240           0 :         return 0;
    1241           0 : }/* end of chooseComponent() */
    1242             : 
    1243             : /* Update the RHS vector - Called from 'updateModelandRHS'.
    1244             : Note : This function is called within the 'scale' omp/pragma loop. Needs to be thread-safe for scales.
    1245             :  */
    1246           0 : Int MultiTermMatrixCleaner::updateRHS(Int ntaylor, Int scale, Float loopgain, Vector<Float> coeffs, IPosition blc, IPosition trc, IPosition blcPsf, IPosition trcPsf)
    1247             : {
    1248           0 :     for(Int taylor1=0;taylor1<ntaylor;taylor1++)
    1249             :     {
    1250           0 :       Matrix<Float> residSub = (matR_p[IND2(taylor1,scale)])(blc,trc);  
    1251           0 :            for(Int taylor2=0;taylor2<ntaylor;taylor2++)
    1252             :            {
    1253           0 :              Matrix<Float> smoothSub = (cubeA_p[IND4(taylor1,taylor2,scale,maxscaleindex_p)])(blcPsf,trcPsf);
    1254           0 :              residSub -= smoothSub * loopgain * coeffs[taylor2];
    1255             :              //      residSub = residSub - smoothSub * loopgain * (matCoeffs_p[IND2(taylor2,maxscaleindex_p)])(globalmaxpos_p);
    1256           0 :            }
    1257           0 :     }
    1258             :     //
    1259           0 :     return 0;
    1260             : }/* end of updateRHS */
    1261             : 
    1262             : 
    1263             : /***************************************
    1264             :  *  Update the model images and the convolved residuals
    1265             : TODO : The current assumption is that only one scale is chosen at a time,
    1266             :                However the chi-sq derivative can be computed across scales too
    1267             :                and the update could be all scales at once for the 'best' location.
    1268             :  ****************************************/
    1269           0 : Int MultiTermMatrixCleaner::updateModelAndRHS(Float loopgain)
    1270             : {
    1271             : 
    1272             :    /* Max support size for all updates is the full image size */
    1273             :    // blc, trc :model image -> needs to be centred on the component location
    1274             :    // blcPsf : psf or scale-blob -> centred on the psf image center (peak).
    1275             : 
    1276             :   /////   gip = IPosition(2,nx_p,ny_p);  
    1277             : 
    1278             :    //IPosition support(2,nx_p,ny_p); // OLD
    1279             :    //IPosition psfpeak(itsPositionPeakPsf);
    1280             :    /* The update region. */
    1281             :    //IPosition inc(2,1,1);
    1282             :    //IPosition blc(globalmaxpos_p - support/2);
    1283             :    //IPosition trc(globalmaxpos_p + support/2 - IPosition(2,1,1));
    1284             :    //LCBox::verify(blc, trc, inc, gip);
    1285             :    /* Shifted region, with the psf at the globalmaxpos_p. */
    1286             :    //IPosition blcPsf(blc + psfpeak - globalmaxpos_p); // OLD
    1287             :    //IPosition trcPsf(trc + psfpeak - globalmaxpos_p); // OLD
    1288             :    ///LCBox::verify(blcPsf, trcPsf, inc, gip); // OLD
    1289             : 
    1290             : 
    1291             :    /* The update region. */
    1292             :    /////IPosition inc(2,1,1);
    1293             :    /////IPosition blc(globalmaxpos_p - psfsupport_p/2);
    1294             :   /////   IPosition trc(globalmaxpos_p + psfsupport_p/2 - IPosition(2,1,1));
    1295             :   /////   LCBox::verify(blc, trc, inc, gip);
    1296             : 
    1297             :    
    1298             :    /* Shifted region, with the psf at the globalmaxpos_p. */
    1299             :   /////   IPosition blcPsf(blc + psfpeak_p - globalmaxpos_p); // OLD
    1300             :   /////   IPosition trcPsf(trc + psfpeak_p - globalmaxpos_p); // OLD
    1301             :   /////   LCBox::verify(blcPsf, trcPsf, inc, psfsupport_p); // NEW
    1302             : 
    1303             :    /* Reconcile box sizes/locations with the image size */
    1304             :   /////   makeBoxesSameSize(blc,trc,blcPsf,trcPsf);
    1305             : 
    1306             :   //buildImagePatches();
    1307             : 
    1308             :    //UUU   
    1309             :    /*
    1310             :    cout << "Source location : " << globalmaxpos_p << endl;
    1311             :    cout << "region around peak residual : " << blc << trc << endl;
    1312             :    cout << "around the PSF peak : " << blcPsf << trcPsf << endl;
    1313             :    cout << "around the Scale blob : " << blcScale << trcScale << endl;
    1314             :    */
    1315             : 
    1316             :  
    1317             :    /* Update the model images */
    1318             :    ///   Matrix<Float> scaleSub = (vecScales_p[maxscaleindex_p])(blcPsf,trcPsf);  // OLD
    1319             :    /// Matrix<Float> scaleSub = (vecScales_p[maxscaleindex_p])(blcScale, trcScale);  // NEW
    1320           0 :    Matrix<Float> scaleSub = (vecScales_p[maxscaleindex_p])(blcPsf_p, trcPsf_p);  // NEWER (same size as psf)
    1321           0 :    for(Int taylor=0;taylor<ntaylor_p;taylor++)
    1322             :    {
    1323           0 :      Matrix<Float> modelSub = (vecModel_p[taylor])(blc_p,trc_p); 
    1324           0 :      modelSub += scaleSub * loopgain * (matCoeffs_p[IND2(taylor,maxscaleindex_p)])(globalmaxpos_p);
    1325           0 :    }
    1326             : 
    1327             :    /* Update the convolved residuals */
    1328           0 :    Vector<Float> coeffs(ntaylor_p);
    1329           0 :    for(Int taylor=0;taylor<ntaylor_p;taylor++) 
    1330           0 :           coeffs[taylor] = (matCoeffs_p[IND2(taylor,maxscaleindex_p)])(globalmaxpos_p);
    1331             :    
    1332             :    Int scale;
    1333           0 :    Int ntaylor=ntaylor_p;
    1334           0 :    IPosition blc(blc_p), trc(trc_p), blcPsf(blcPsf_p), trcPsf(trcPsf_p);
    1335             :    //OMP// #pragma omp parallel default(shared) private(scale) firstprivate(ntaylor,loopgain,coeffs,blc,trc,blcPsf,trcPsf)
    1336             :   { 
    1337             :     //OMP// #pragma omp for 
    1338           0 :     for(scale=0;scale<nscales_p;scale++)
    1339             :    {
    1340           0 :      updateRHS(ntaylor,scale, loopgain, coeffs, blc, trc, blcPsf, trcPsf);
    1341             :    }
    1342             :   }//End pragma parallel  
    1343             : 
    1344             :    /* Update flux counters */
    1345           0 :    for(Int taylor=0;taylor<ntaylor_p;taylor++)
    1346             :    {
    1347           0 :            totalTaylorFlux_p[taylor] += loopgain*(matCoeffs_p[IND2(taylor,maxscaleindex_p)])(globalmaxpos_p);
    1348             :    }
    1349           0 :    totalScaleFlux_p[maxscaleindex_p] += loopgain*(matCoeffs_p[IND2(0,maxscaleindex_p)])(globalmaxpos_p);
    1350             :    
    1351           0 :    return 0;
    1352           0 : }/* end of updateModelAndRHS() */
    1353             : 
    1354             : /* ................ */
    1355           0 : Int MultiTermMatrixCleaner::checkConvergence(Int /*criterion*/, Float &fluxlimit, Float &loopgain)
    1356             : {
    1357           0 :     Float rmaxval=0.0;
    1358             :     
    1359             :     /* Use the maximum residual (current), to compare against the convergence threshold */
    1360           0 :     Float maxres=0.0;
    1361           0 :     IPosition maxrespos;
    1362             : 
    1363           0 :     findMaxAbsMask((matR_p[IND2(0,0)]),vecScaleMasks_p[0],maxres,maxrespos);
    1364           0 :     Float norma = (1.0/(matA_p[0])(0,0));
    1365           0 :     rmaxval = abs(maxres*norma);
    1366           0 :     rmaxval_p = fabs(rmaxval);
    1367             : 
    1368             :     /* // Calc the max residual across all scales....
    1369             :     Int maxscale=0;
    1370             :     Float maxscaleresidual=0.0;
    1371             :     for (Int scale =0; scale<nscales_p; scale++)
    1372             :       {
    1373             :         findMaxAbsMask((matR_p[IND2(0,scale)]),vecScaleMasks_p[scale],maxres,maxrespos);
    1374             :         if ( maxscaleresidual < maxres )
    1375             :           {
    1376             :             maxscaleresidual = maxres;
    1377             :             maxscale = scale;
    1378             :           }
    1379             :       }
    1380             :     Float norma = (1.0/(matA_p[maxscale])(0,0));
    1381             :     rmaxval = abs(maxscaleresidual*norma);
    1382             :     */
    1383             : 
    1384             :     /* Check for convergence */
    1385           0 :     Int convergedflag = 0;
    1386             :     ///    cout << "MTFT::checkconvergence : maxval : " << fabs(rmaxval) << "  userthreshold : " << userthreshold_p << "    fluxlimit : " << fluxlimit << endl;
    1387           0 :     if( fabs(rmaxval) < MAX(userthreshold_p,fluxlimit) ) 
    1388             :     {
    1389           0 :       LogIO os(LogOrigin("MultiTermMatrixCleaner", "mtclean()", WHERE));
    1390           0 :       os << "Reached stopping threshold at iteration " << totalIters_p << ". Peak residual " << fabs(rmaxval) ;
    1391           0 :       if( ! itsMask.null() ){os << " (within mask) " ;}
    1392           0 :         os << LogIO::POST;
    1393           0 :        convergedflag=1; 
    1394           0 :     }
    1395             : 
    1396             :     /* Levenberg-Macquart-like change in step size */
    1397           0 :     if(itercount_p>1 && inputgain_p<=(Float)0.0)
    1398             :    { 
    1399             :      
    1400           0 :         if (globalmaxval_p < prev_max_p) 
    1401           0 :             loopgain = loopgain * 1.5;
    1402             :         else 
    1403           0 :             loopgain = loopgain / 1.5;
    1404             : 
    1405           0 :         loopgain = MIN((1-stopfraction_p),loopgain);
    1406           0 :         loopgain = MIN((Float)0.6,loopgain);
    1407             :         
    1408             :         /* detect divergence by approximately 10 consecutive increases in maxval */
    1409           0 :         if(loopgain < (Float)0.01)
    1410             :         {
    1411           0 :           LogIO os(LogOrigin("MultiTermMatrixCleaner", "mtclean()", WHERE));
    1412           0 :           os << "Not converging any more. May be diverging. Stopping minor cycles at iteration " << totalIters_p << ". Peak residual " << fabs(rmaxval) << LogIO::POST;
    1413           0 :           convergedflag=-1; 
    1414           0 :          }
    1415             :      
    1416             :     /* Stop if there is divergence : 200% increase from the current minimum value */
    1417           0 :         if( fabs(  (min_max_p-globalmaxval_p)/min_max_p ) > (Float)2.0 )
    1418             :       {
    1419           0 :         LogIO os(LogOrigin("MultiTermMatrixCleaner", "mtclean()", WHERE));
    1420           0 :         os << "Diverging.... Stopping minor cycles at iteration " << totalIters_p << ". Peak residual " << fabs(rmaxval) << " Max: " << globalmaxval_p << LogIO::POST;
    1421           0 :         convergedflag=-1;
    1422           0 :       }
    1423             : 
    1424             :     // Stop if the maxval has changed by less than 5% in 5 iterations 
    1425             :     // --- this is similar to saying less than 1% change per iteration.... same as a loopgain < 0.01 
    1426             :         
    1427             :         //      if( abs(prev5_max - globalmaxval_p) < 0.01*abs(prev5_max) )
    1428             :         //  { 
    1429             :         //    os << "Not converging any more. Flattening out. Stopping minor cycles at iteration " << totalIters_p << ". Peak residual " << fabs(rmaxval) << LogIO::POST;
    1430             :         //    convergedflag=-1;
    1431             :         //  }
    1432             : 
    1433             : 
    1434             :         
    1435             :     }// end of if(itercount_p>1)
    1436             : 
    1437             :     /* Store current max value - to use in the next iteration */    
    1438           0 :     prev_max_p = globalmaxval_p;
    1439             :     //    if(itercount_p%5 == 0) 
    1440             :     //     prev5_max=globalmaxval_p;
    1441           0 :     min_max_p = MIN(min_max_p,abs(globalmaxval_p));
    1442             : 
    1443             :     /* Stop, if there are negatives on the largest scale in the Io image */
    1444             :     //if(nscales_p>1 && maxscaleindex_p == nscales_p-2)
    1445             :     //  if((*matCoeffs_p[IND2(0,maxscaleindex_p)]).getAt(globalmaxpos_p) < 0.0)
    1446             :     //  {converged = false;break;}
    1447             : 
    1448             :     /* Detect divergence, and signal it.... */
    1449             :     // TODO
    1450             : 
    1451             :     /* Print out coefficients for a few iterations */
    1452           0 :     if(convergedflag==0)
    1453             :       {
    1454           0 :     if(fluxlimit==-1) 
    1455             :     {
    1456           0 :          fluxlimit = rmaxval * stopfraction_p;
    1457             :          //        os << "Peak convolved residual : " << rmaxval << "    Minor cycle stopping threshold : " << itsThreshold.getValue("Jy")  << LogIO::POST;
    1458           0 :          LogIO os(LogOrigin("MultiTermMatrixCleaner", "mtclean()", WHERE));
    1459           0 :          os << "Peak convolved residual" ;
    1460           0 :          if( ! itsMask.null() ){os << " (within mask) ";}
    1461           0 :          os << " : " << rmaxval;
    1462           0 :          if( fluxlimit > 0.0 ){ os << "  : Minor cycle stopping threshold : " << fluxlimit;}
    1463           0 :          os << LogIO::POST;
    1464           0 :     }
    1465             :     else
    1466             :     {
    1467             :       //       if(1)
    1468           0 :         if( totalIters_p==maxniter_p || (adbg==(Bool)true) || maxniter_p < (int)5 || (totalIters_p%(Int)20==0) )
    1469             :        {
    1470             :          
    1471           0 :             os << "[" << totalIters_p << "] Res: " << rmaxval << " Max: " << globalmaxval_p;
    1472           0 :             os << " Gain: " << loopgain;
    1473             :              ////        os << "[" << totalIters_p << "] Res: " << rmaxval;
    1474             :             //os << "[" << totalIters_p << "] Max: " << globalmaxval_p;
    1475           0 :             os << " Pos: " <<  globalmaxpos_p << " Scale: " << scaleSizes_p[maxscaleindex_p];
    1476           0 :             os << " Coeffs: ";
    1477           0 :             for(Int taylor=0;taylor<ntaylor_p;taylor++)
    1478           0 :                os << (matCoeffs_p[IND2(taylor,maxscaleindex_p)])(globalmaxpos_p) << "  ";
    1479           0 :             if(adbg)
    1480             :               {
    1481           0 :               os << " OrigRes: ";
    1482           0 :               for(Int taylor=0;taylor<ntaylor_p;taylor++)
    1483           0 :                    os << (matR_p[IND2(taylor,maxscaleindex_p)])(globalmaxpos_p) << "  ";
    1484             :               }
    1485           0 :             os << LogIO::POST;
    1486             : 
    1487             :             //UUU cout << "minor " << totalIters_p << " " << (vecModel_p[0])( IPosition(2,nx_p/2, ny_p/2) ) << " " << (vecModel_p[1])( IPosition(2,nx_p/2, ny_p/2) ) << "  " << (matR_p[0])( IPosition(2,nx_p/2, ny_p/2) ) <<  endl;
    1488             : 
    1489             :         }
    1490             :     }
    1491             : 
    1492             :       }// if converged-flag is still 0
    1493             : 
    1494             : 
    1495           0 :     return convergedflag;
    1496             : 
    1497           0 : }/* end of checkConvergence */
    1498             : 
    1499             : 
    1500             : 
    1501             : /* Save a matrix to disk */
    1502           0 : Int MultiTermMatrixCleaner::writeMatrixToDisk(String imagename, Matrix<Float>& themat)
    1503             : {
    1504           0 :   TabularCoordinate tab1(0, 1.0, 0, String("deg"), String("axis1"));
    1505           0 :   TabularCoordinate tab2(0, 1.0, 0, String("deg"), String("axis2"));
    1506           0 :       CoordinateSystem csys;
    1507           0 :       csys.addCoordinate(tab1);
    1508           0 :       csys.addCoordinate(tab2);
    1509           0 :       PagedImage<Float> limage(themat.shape(), csys, imagename);
    1510           0 :       limage.put(themat);
    1511             :       //return value does not seemed to be used so will make compiler happy
    1512           0 :       return 1;
    1513           0 : }
    1514             : 
    1515             : 
    1516             : /* Compute principal solution in-place on the list of residual images ( vecDirty ) 
    1517             :     -- Call MTMC::setresidual() repeatedly to fill in final residuals before computing
    1518             :        the principal solution. This does the same as solveMatrixEquation(), but 
    1519             :        stores the results in-place in the residual images */
    1520           0 : Bool MultiTermMatrixCleaner::computeprincipalsolution()
    1521             : {
    1522           0 :     LogIO os(LogOrigin("MultiTermMatrixCleaner", "computeprincipalsolution()", WHERE));
    1523             : 
    1524           0 :     os << "MTMC :: Computing principal solution on residuals" << LogIO::POST;
    1525             : 
    1526             :     /// If this is being called with niter=0, the Hessian won't exist. So, make it.
    1527           0 :     if( doneCONV_p == false )
    1528             :     {
    1529           0 :         if( computeHessianPeak() == -2 )
    1530           0 :             return false;
    1531             :     }
    1532             : 
    1533           0 :     AlwaysAssert((vecDirty_p.nelements()>0), AipsError);
    1534             : 
    1535             :     /* Solve for the coefficients at the delta-function scale*/
    1536           0 :     for(Int taylor1=0;taylor1<ntaylor_p;taylor1++)
    1537             :     {
    1538           0 :         (matCoeffs_p[IND2(taylor1,0)]) = 0.0; 
    1539           0 :         for(Int taylor2=0;taylor2<ntaylor_p;taylor2++)
    1540             :         {
    1541           0 :             matCoeffs_p[IND2(taylor1,0)] = matCoeffs_p[IND2(taylor1,0)] + ((Float)(invMatA_p[0])(taylor1,taylor2))*(vecDirty_p[taylor2]);
    1542             :         }
    1543             :     }
    1544             :         
    1545             :     /* Copy this into the residual vector */
    1546           0 :     for(Int taylor=0; taylor<ntaylor_p;taylor++)
    1547             :     {
    1548           0 :         vecDirty_p[taylor].assign(matCoeffs_p[IND2(taylor,0)]);
    1549             :     }
    1550             : 
    1551           0 :     return true;
    1552           0 : }
    1553             : 
    1554             : 
    1555             : 
    1556             : } //# NAMESPACE CASA - END

Generated by: LCOV version 1.16