//# SynthesisImagerVi2.cc: Implementation of SynthesisImager.h
//# Copyright (C) 1997-2021
//# Associated Universities, Inc. Washington DC, USA.
//# This library is free software; you can redistribute it and/or modify it
//# under the terms of the GNU General Public License as published by
//# the Free Software Foundation; either version 3 of the License, or (at your
//# option) any later version.
//#
//# This library is distributed in the hope that it will be useful, but WITHOUT
//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
//# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
//# License for more details.
//#
//# https://www.gnu.org/licenses/
//#
//# You should have received a copy of the GNU  General Public License
//# along with this library; if not, write to the Free Software Foundation,
//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
//#
//# Queries concerning CASA should be submitted at
//#        https://help.nrao.edu
//#
//#        Postal address: CASA Project Manager 
//#                        National Radio Astronomy Observatory
//#                        520 Edgemont Road
//#                        Charlottesville, VA 22903-2475 USA
//#
//#
//# $Id$

#define CFC_VERBOSE false /* Control the verbosity when building CFCache. */

#include <casacore/casa/Exceptions/Error.h>
#include <iostream>
#include <sstream>

#include <casacore/casa/Arrays/Matrix.h>
#include <casacore/casa/Arrays/ArrayMath.h>
#include <casacore/casa/Arrays/ArrayLogical.h>


#include <casacore/casa/Logging.h>
#include <casacore/casa/Logging/LogIO.h>
#include <casacore/casa/Logging/LogMessage.h>
#include <casacore/casa/Logging/LogSink.h>
#include <casacore/casa/Logging/LogMessage.h>
#include <casacore/casa/System/ProgressMeter.h>

#include <casacore/casa/OS/DirectoryIterator.h>
#include <casacore/casa/OS/File.h>
#include <casacore/casa/OS/HostInfo.h>
#include <casacore/casa/OS/Path.h>
//#include <casa/OS/Memory.h>

#include <casacore/lattices/LRegions/LCBox.h>

#include <casacore/measures/Measures/MeasTable.h>

#include <casacore/ms/MeasurementSets/MSHistoryHandler.h>
#include <casacore/ms/MeasurementSets/MeasurementSet.h>
#include <casacore/ms/MSSel/MSSelection.h>


#include <synthesis/ImagerObjects/SynthesisImagerVi2.h>

#include <synthesis/ImagerObjects/SynthesisUtilMethods.h>
#include <synthesis/ImagerObjects/SIImageStore.h>
#include <synthesis/ImagerObjects/SIImageStoreMultiTerm.h>
#include <synthesis/ImagerObjects/CubeMajorCycleAlgorithm.h>
#include <synthesis/ImagerObjects/CubeMakeImageAlgorithm.h>
#include <synthesis/MeasurementEquations/VPManager.h>
#include <imageanalysis/Utilities/SpectralImageUtil.h>
#include <msvis/MSVis/MSUtil.h>
#include <msvis/MSVis/VisSetUtil.h>
#include <msvis/MSVis/VisImagingWeight.h>

#include <synthesis/TransformMachines2/GridFT.h>
#include <synthesis/TransformMachines2/WPConvFunc.h>
#include <synthesis/TransformMachines2/WProjectFT.h>
#include <synthesis/TransformMachines2/VisModelData.h>
#include <synthesis/TransformMachines2/AWProjectFT.h>
#include <synthesis/TransformMachines2/HetArrayConvFunc.h>
#include <synthesis/TransformMachines2/MosaicFTNew.h>
#include <synthesis/TransformMachines2/AWPLPG.h>
#include <synthesis/TransformMachines2/MultiTermFTNew.h>
#include <synthesis/TransformMachines2/AWProjectFT.h>
#include <synthesis/TransformMachines2/AWProjectWBFTHPG.h>
#include <synthesis/TransformMachines2/AWVisResamplerHPG.h>
#include <synthesis/TransformMachines2/AWConvFunc.h>
#include <synthesis/TransformMachines2/SDGrid.h>
#include <synthesis/TransformMachines/WProjectFT.h>
#include <synthesis/TransformMachines2/BriggsCubeWeightor.h>
#include <casacore/casa/Utilities/Regex.h>
#include <casacore/casa/OS/Directory.h>
#include <msvis/MSVis/VisibilityIteratorImpl2.h>
//#include <casadbus/utilities/BusAccess.h>
//#include <casadbus/session/DBusSession.h>

#include <sys/types.h>
#include <unistd.h>
#include <iomanip>

#include <chrono>
#include <thread>
#include <synthesis/Parallel/Applicator.h>

#ifdef USE_HPG
#include <hpg/hpg.hpp>
#include <dlfcn.h>
#endif


    using namespace std;

using namespace casacore;

namespace casa { //# NAMESPACE CASA - BEGIN
  extern Applicator applicator;
  SynthesisImagerVi2::SynthesisImagerVi2() : SynthesisImager(), vi_p(0), fselections_p(nullptr), imparsVec_p(0), gridparsVec_p(0) {
	/*cerr << "is applicator initialized " << applicator.initialized() << endl;
	if(!applicator.initialized()){
	  int argc=1;
        char **argv;
        casa::applicator.init ( argc, argv );
        cerr << "controller ?" <<  applicator.isController() <<  " worker? " <<  applicator.isWorker() <<  " numprocs " << applicator.numProcs() <<  endl;
	}
	*/
    mss_p.resize(0);
  }

  SynthesisImagerVi2::~SynthesisImagerVi2(){
    for (uInt k=0; k < mss_p.nelements(); ++k){
      if(mss_p[k])
	delete mss_p[k];
    }
     
    //if(gridparsVec_p[0].ftmachine=="awphpg")
    //  hpg::finalize();
      SynthesisUtilMethods::getResource("End Run");
}

  Bool SynthesisImagerVi2::selectData(const SynthesisParamsSelect& selpars){
 LogIO os( LogOrigin("SynthesisImagerVi2","selectData",WHERE) );
 Bool retval=True;

    SynthesisUtilMethods::getResource("Start Run");
    
    try
      {


	MeasurementSet thisms;
	{ ///Table system seems to have a bug when running in multi-process as the source table disappears
	  /// temporarily when other processes are updating it 
	  uInt exceptCounter=0;
	  
	  while(true){
	    try{
	      //Respect the readonly flag...necessary for multi-process access
	      thisms=MeasurementSet(selpars.msname, TableLock(TableLock::UserNoReadLocking),
				    selpars.readonly ? Table::Old : Table::Update);
	      break;
	    }
	    catch(AipsError &x){
	      
	      String mes=x.getMesg();
	      if(mes.contains("FilebufIO::readBlock") || mes.contains("SOURCE")){
		std::this_thread::sleep_for(std::chrono::milliseconds(50));
		os << LogIO::WARN << "#####CATCHING a sleep because "<< mes<< LogIO::POST;
	      }
	      else
		throw(AipsError("Error in selectdata: "+mes));
	      
	      if(exceptCounter > 100){
		throw(AipsError("Error in selectdata got 100 of this exeception: "+mes));
		
	      }
	      
	    }
	    ++exceptCounter;
	  }
	}//End of work around for table disappearing bug
	

    
    useScratch_p=selpars.usescratch;
    readOnly_p = selpars.readonly;
    lockMS(thisms);	
    thisms.setMemoryResidentSubtables (MrsEligibility::defaultEligible());

    
    //    cout << "**************** usescr : " << useScratch_p << "     readonly : " << readOnly_p << endl;
    //if you want to use scratch col...make sure they are there
    ///Need to clear this in first pass only...child processes or loops for cube ///should skip it
    if(!(impars_p.mode.contains("cube")) || ((impars_p.mode.contains("cube")) && doingCubeGridding_p)){
      if(selpars.usescratch && !selpars.readonly){
	VisSetUtil::addScrCols(thisms, true, false, true, false);
	refim::VisModelData::clearModel(thisms);
      }
    ////TESTOO
    //Int CPUID;
	//MPI_Comm_rank(MPI_COMM_WORLD, &CPUID);
    //cerr  << " SELPARS " << selpars.toRecord()  << endl;
      if(!selpars.incrmodel && !selpars.usescratch && !selpars.readonly)
	refim::VisModelData::clearModel(thisms, selpars.field, selpars.spw);
    }//end of first pass if for cubes
    // unlockMSs();


    os << "MS : " << selpars.msname << " | ";

    //Some MSSelection 
    //If everything is empty (which is valid) it will throw an exception..below
    //So make sure the main defaults are not empy i.e field and spw
    MSSelection thisSelection;
    if(selpars.field != ""){
      thisSelection.setFieldExpr(selpars.field);
      os << "Selecting on fields : " << selpars.field << " | " ;//LogIO::POST;
    }else
      thisSelection.setFieldExpr("*");
    if(selpars.spw != ""){
	thisSelection.setSpwExpr(selpars.spw);
	os << "Selecting on spw :"<< selpars.spw  << " | " ;//LogIO::POST;
    }else
      thisSelection.setSpwExpr("*");
    
    if(selpars.antenna != ""){
      Vector<String> antNames(1, selpars.antenna);
      // thisSelection.setAntennaExpr(MSSelection::nameExprStr( antNames));
      thisSelection.setAntennaExpr(selpars.antenna);
      os << "Selecting on antenna names : " << selpars.antenna << " | " ;//LogIO::POST;
	
    }            
    if(selpars.timestr != ""){
	thisSelection.setTimeExpr(selpars.timestr);
	os << "Selecting on time range : " << selpars.timestr << " | " ;//LogIO::POST;	
      }
    if(selpars.uvdist != ""){
      thisSelection.setUvDistExpr(selpars.uvdist);
      os << "Selecting on uvdist : " << selpars.uvdist << " | " ;//LogIO::POST;	
    }
    if(selpars.scan != ""){
      thisSelection.setScanExpr(selpars.scan);
      os << "Selecting on scan : " << selpars.scan << " | " ;//LogIO::POST;	
    }
    if(selpars.obs != ""){
      thisSelection.setObservationExpr(selpars.obs);
      os << "Selecting on Observation Expr : " << selpars.obs << " | " ;//LogIO::POST;	
    }
    if(selpars.state != ""){
      thisSelection.setStateExpr(selpars.state);
      os << "Selecting on Scan Intent/State : " << selpars.state << " | " ;//LogIO::POST;	
    }
    // if(selpars.taql != ""){
    // 	thisSelection.setTaQLExpr(selpars.taql);
    // 	os << "Selecting via TaQL : " << selpars.taql << " | " ;//LogIO::POST;	
    // }
    os << "[Opened " << (readOnly_p?"in readonly mode":(useScratch_p?"with scratch model column":"with virtual model column"))  << "]" << LogIO::POST;
    TableExprNode exprNode=thisSelection.toTableExprNode(&thisms);
    if(!(exprNode.isNull()))
      {
	
    
	MeasurementSet thisMSSelected0 = MeasurementSet(thisms(exprNode));
	mss_p.resize(mss_p.nelements()+1, false, true);
	if(selpars.taql != "")
	  {
	    MSSelection mss0;
	    mss0.setTaQLExpr(selpars.taql);
	    os << "Selecting via TaQL : " << selpars.taql << " | " ;//LogIO::POST;	

	    TableExprNode tenWithTaQL=mss0.toTableExprNode(&thisMSSelected0);
	    MeasurementSet thisMSSelected1 = MeasurementSet(thisMSSelected0(tenWithTaQL));
	    //mss4vi_p[mss4vi_p.nelements()-1]=MeasurementSet(thisms(exprNode));
	    mss_p[mss_p.nelements()-1]=new MeasurementSet(thisMSSelected1);
	  }
	else
	  mss_p[mss_p.nelements()-1]=new MeasurementSet(thisMSSelected0);
	  
	os << "  NRows selected : " << (mss_p[mss_p.nelements()-1])->nrow() << LogIO::POST;
	unlockMSs();
      }
    else{
      throw(AipsError("Selection for given MS "+selpars.msname+" is invalid"));
    }
    if((mss_p[mss_p.nelements()-1])->nrow() ==0){
      delete mss_p[mss_p.nelements()-1];
      mss_p.resize(mss_p.nelements()-1, True, True);
      if(mss_p.nelements()==0)
	throw(AipsError("Data selection ended with 0 rows"));
      //Sill have some valid ms's so return false and do not proceed to do 
      //channel selection
      unlockMSs();
      return False;
    }


    
    ///Channel selection
    {
      if(!fselections_p) fselections_p=new FrequencySelections();
      Matrix<Int> chanlist = thisSelection.getChanList(mss_p[mss_p.nelements()-1]);
      Matrix<Double> freqList=thisSelection.getChanFreqList(mss_p[mss_p.nelements()-1]);
      //cerr << std::setprecision(12) << "FreqList " << freqList << endl;
      IPosition shape = freqList.shape();
      uInt nSelections = shape[0];
      ///temporary variable as we carry that for tunechunk...till we get rid of it
      selFreqFrame_p=selpars.freqframe;
      Bool ignoreframe=False;
      MFrequency::Types freqFrame=MFrequency::castType(MSColumns(*mss_p[mss_p.nelements()-1]).spectralWindow().measFreqRef()(Int(chanlist(0,0))));
  
      if(selpars.freqframe == MFrequency::REST ||selpars.freqframe == MFrequency::Undefined){	
	selFreqFrame_p=freqFrame;
	ignoreframe=True;
      }
      if(selpars.freqbeg==""){
	   // Going round the problem of CAS-8829
        /*vi::FrequencySelectionUsingChannels channelSelector;

        channelSelector.add(thisSelection, mss_p[mss_p.nelements()-1]);

        fselections_p.add(channelSelector);
        */
        ////////////////////////////
        Double lowfreq;
        Double topfreq;
      Vector<Int> fieldList=thisSelection.getFieldList(mss_p[mss_p.nelements()-1]);
	 // cerr << "chanlist " << chanlist.column(0) << "\n fieldList " << fieldList << endl;
        
	//cerr << "selpars.freqframe " << selpars.freqframe << endl;
        vi::FrequencySelectionUsingFrame channelSelector(selFreqFrame_p);
	///temporary variable as we carry that for tunechunk
		
	Bool selectionValid=False;
    	  for(uInt k=0; k < nSelections; ++k){
	    Bool thisSpwSelValid=False;
	    //The getChanfreqList is wrong for beg and end..going round that too.
	    Vector<Double> freqies=MSColumns(*mss_p[mss_p.nelements()-1]).spectralWindow().chanFreq()(Int(chanlist(k,0)));
	    Vector<Double> chanwidth=MSColumns(*mss_p[mss_p.nelements()-1]).spectralWindow().chanWidth()(Int(chanlist(k,0)));
            
	    if(freqList(k,3) < 0.0){
	      topfreq=freqies(chanlist(k,1));//-chanwidth(chanlist(k,1))/2.0;
	      lowfreq=freqies(chanlist(k,2));//+chanwidth(chanlist(k,2))/2.0;
	      //lowfreq=freqList(k,2); //+freqList(k,3)/2.0;
	      //topfreq=freqList(k, 1); //-freqList(k,3)/2.0;
	    }
	    else{
	      lowfreq=freqies(chanlist(k,1));//-chanwidth(chanlist(k,1))/2.0;
	      topfreq=freqies(chanlist(k,2));//+chanwidth(chanlist(k,2))/2.0;
	      //lowfreq=freqList(k,1);//-freqList(k,3)/2.0;
	      //topfreq=freqList(k, 2);//+freqList(k,3)/2.0;
	    }
	    
	    if(!ignoreframe){
		
			
          //cerr << "begin " << lowfreq << "  " << topfreq << endl; 
	      //vi::VisibilityIterator2 tmpvi(mss_p, vi::SortColumns(), false); 
	      //VisBufferUtil::getFreqRangeFromRange(lowfreq, topfreq,  freqFrame, lowfreq,  topfreq, tmpvi, selFreqFrame_p);
	      //    lockMS(*(const_cast<MeasurementSet*> (mss_p[mss_p.nelements()-1])));
              if(MSUtil::getFreqRangeInSpw( lowfreq,
				  topfreq, Vector<Int>(1,chanlist(k,0)), Vector<Int>(1,chanlist(k,1)),
				  Vector<Int>(1, chanlist(k,2)-chanlist(k,1)+1),
				 *mss_p[mss_p.nelements()-1] , 
				  selFreqFrame_p,
						 fieldList, False))
                {
                  selectionValid=True;
                  thisSpwSelValid=True;
                }
	      //  unlockMSs();
		    
	    }
	    
	    if(thisSpwSelValid || ignoreframe){
	      andFreqSelection(mss_p.nelements()-1, Int(freqList(k,0)), lowfreq, topfreq, selFreqFrame_p);
	      andChanSelection(mss_p.nelements()-1, Int(chanlist(k,0)), Int(chanlist(k,1)),Int(chanlist(k,2)));
	    }
          }
	  if(! (selectionValid && !ignoreframe)){
	    //os << "Did not match spw selection in the selected ms " << LogIO::WARN << LogIO::POST;
	    retval=False;
	  }
	    //fselections_p->add(channelSelector);
          //////////////////////////////////
      }
      else{

	//////More workaroung CAS-8829
	//MFrequency::Types freqFrame=MFrequency::castType(MSColumns(*mss_p[mss_p.nelements()-1]).spectralWindow().measFreqRef()(Int(freqList(0,0))));
    	  Quantity freq;
    	  Quantity::read(freq, selpars.freqbeg);
    	  Double lowfreq=freq.getValue("Hz");
    	  Quantity::read(freq, selpars.freqend);
    	  Double topfreq=freq.getValue("Hz");
    	 
	  ////Work aroun CAS-8829
	  // if(vi_p) 
	    //VisBufferUtil::getFreqRangeFromRange(lowfreq, topfreq,  selpars.freqframe, lowfreq,  topfreq, *vi_p, freqFrame);
	  //cerr << "lowFreq "<< lowfreq << " topfreq " << topfreq << endl;
	  //vi::FrequencySelectionUsingFrame channelSelector((vi_p ? freqFrame :selpars.freqframe));
	  //vi::FrequencySelectionUsingFrame channelSelector(selpars.freqframe);
    	  for(uInt k=0; k < nSelections; ++k){
            //channelSelector.add(Int(freqList(k,0)), lowfreq, topfreq);
	    //andFreqSelection((mss_p.nelements()-1), Int(freqList(k,0)), lowfreq, topfreq, vi_p ?freqFrame : selpars.freqframe);
	    andFreqSelection((mss_p.nelements()-1), Int(freqList(k,0)), lowfreq, topfreq, selFreqFrame_p);
          }
    	  //fselections_p->add(channelSelector);

      }
    }//End of channel selection
          
    writeAccess_p=writeAccess_p && !selpars.readonly;
    createVisSet(writeAccess_p);

    /////// Remove this when the new vi/vb is able to get the full freq range.
    mssFreqSel_p.resize();
    mssFreqSel_p  = thisSelection.getChanFreqList(NULL,true);
   
    //// Set the data column on which to operate
    // TT: added checks for the requested data column existace 
    //    cout << "Using col : " << selpars.datacolumn << endl;
    if( selpars.datacolumn.contains("data") || selpars.datacolumn.contains("obs") ) 
      {    if( thisms.tableDesc().isColumn("DATA") ) { datacol_p = FTMachine::OBSERVED; }
           else { os << LogIO::SEVERE <<"DATA column does not exist" << LogIO::EXCEPTION;}
      }
    else if( selpars.datacolumn.contains("corr") ) { datacol_p = FTMachine::CORRECTED; }
    else { os << LogIO::WARN << "Invalid data column : " << selpars.datacolumn << ". Using corrected (or observed if corrected doesn't exist)" << LogIO::POST;  datacol_p =  FTMachine::CORRECTED;}


    dataSel_p.resize(dataSel_p.nelements()+1, true);
    dataSel_p[dataSel_p.nelements()-1]=selpars;
	//cerr << "AT THE end of DATASEL " << selpars.toRecord() << endl;

    unlockMSs();
      }
    catch(AipsError &x)
      {
	unlockMSs();
	throw( AipsError("Error in selectData() : "+x.getMesg()) );
      }

    return retval;



  }
void SynthesisImagerVi2::andChanSelection(const Int msId, const Int spwId, const Int startchan, const Int endchan){

	map<Int, Vector<Int> > spwsel;
	auto it=channelSelections_p.find(msId);
	if(it !=channelSelections_p.end())
		spwsel=it->second;
	auto hasspw=spwsel.find(spwId);
	Vector<Int>chansel(2,-1);
	if(hasspw != spwsel.end()){
		chansel.resize();
		chansel=hasspw->second;
	}
	Int nchan=endchan-startchan+1;
	if(chansel(1)== -1)
		chansel(1)=startchan;
	if(chansel(1) >= startchan){
	  if(nchan > (chansel(1)-startchan+chansel(0))){
			chansel(0)=nchan;
	  }
	  else{
			chansel(0)=chansel(1)-startchan+chansel(0);
	  }
	  chansel(1)=startchan;
	}
	else{
		if((chansel(0) -(startchan - chansel(1)+1)) < nchan){	
		  chansel(0)=nchan+(startchan-chansel(1));
		}
	}
	spwsel[spwId]=chansel;
	channelSelections_p[msId]=spwsel;
	//cerr << "chansel "<< channelSelections_p << endl;
	
}
  void SynthesisImagerVi2::andFreqSelection(const Int msId, const Int spwId,  const Double freqBeg, const Double freqEnd, const MFrequency::Types frame){
    
    
    Int key=msId;
   
    Bool isDefined=False;
    FrequencySelectionUsingFrame frameSel(frame);
    for (uInt k =0; k<freqBegs_p.size(); ++k){ 
      //cerr <<freqBegs_p[k].first  << " == " << key << " && " << freqSpws_p[k].second<< " == " << spwId << " && " << freqBeg << " < " << freqEnds_p[k].second<< " && " << freqEnd << " > " << freqBegs_p[k].second << endl;
	if((freqBegs_p[k].first == key || key <0 ) && (freqSpws_p[k].second==spwId || spwId <0)  && (freqBeg < freqEnds_p[k].second) && (freqEnd > freqBegs_p[k].second)){
	isDefined=True;
	//cerr << k << " inside freqBegs " << freqBegs_p[k].second << "  " << freqBeg << endl;  
	if(freqBegs_p[k].second < freqBeg)
	  freqBegs_p[k].second=freqBeg;
	if(freqEnds_p[k].second > freqEnd)
	  freqEnds_p[k].second=freqEnd;
	if(msId < 0) key=freqBegs_p[k].first;
	//cerr << "modified " <<  freqBegs_p[k].second << "   "  <<  freqEnds_p[k].second << endl;
      }
	///add only those that have the same msid
	if(freqBegs_p[k].first == key){
	  //cerr << "added " << k << " freqBegs " << freqBegs_p[k].second << "  " << freqEnds_p[k].second << endl;  
	  frameSel.add(freqSpws_p[k].second ,  freqBegs_p[k].second, freqEnds_p[k].second);
	}
    }
    if(!isDefined && msId >=0){
      //cerr << "undefined " << key << " freqBegs "  << freqBeg << "  " << freqEnd << endl;  
      freqBegs_p.push_back(make_pair(key, freqBeg));
      freqEnds_p.push_back(make_pair(key, freqEnd));
      freqSpws_p.push_back(make_pair(key, spwId));
      frameSel.add(spwId,  freqBeg, freqEnd);
    }
    CountedPtr<vi::FrequencySelections> copyFsels=fselections_p->clone();
    uInt nMSs=copyFsels->size() <=msId ? msId+1 : copyFsels->size();
    //cerr << "nms " << nMSs << endl;
    fselections_p=new FrequencySelections();
    for (uInt k=0;  k < nMSs ; ++k){
      if(k==uInt(key)){
	fselections_p->add(frameSel);
	//cerr <<"adding framesel " << frameSel.toString() << endl;
      }
      else{
	const FrequencySelectionUsingFrame& thissel= static_cast<const FrequencySelectionUsingFrame &> (copyFsels->get(k));
	//cerr <<"framesel orig " << thissel.toString() << endl;
	fselections_p->add(thissel);

      }
    }
    
 

  }

  void SynthesisImagerVi2::tuneChunk(const Int gmap) {

    CoordinateSystem cs=itsMappers.imageStore(gmap)->getCSys();
    IPosition imshape=itsMappers.imageStore(gmap)->getShape();
    /////For some reason imagestore returns 0 channel image sometimes
    ////
    if(imshape(3) < 1) {
      return;
    }
   
    Double minFreq=SpectralImageUtil::worldFreq(cs, 0.0);
    Double maxFreq=SpectralImageUtil::worldFreq(cs,imshape(3)-1);
   
    if(maxFreq < minFreq){
      Double tmp=minFreq;
      minFreq=maxFreq;
      maxFreq=tmp;
    }
    
    Int spectralIndex=cs.findCoordinate(Coordinate::SPECTRAL);
    SpectralCoordinate spectralCoord=cs.spectralCoordinate(spectralIndex);
    maxFreq+=fabs(spectralCoord.increment()(0))/2.0;
    minFreq-=fabs(spectralCoord.increment()(0))/2.0;
    if(minFreq < 0.0) minFreq=0.0;
    MFrequency::Types intype=spectralCoord.frequencySystem(True);
    
    if(!VisBufferUtil::getFreqRangeFromRange(minFreq, maxFreq,  intype, minFreq,  maxFreq, *vi_p, selFreqFrame_p)){
      //Do not retune if conversion did not happen
      return;
    }
      
    maxFreq+=3.0*fabs(spectralCoord.increment()(0))/2.0;
    minFreq-=3.0*fabs(spectralCoord.increment()(0))/2.0;
    if(minFreq < 0.0) minFreq=0.0;
    
    auto copyFreqBegs=freqBegs_p;
    auto copyFreqEnds=freqEnds_p;
    auto copyFreqSpws=  freqSpws_p;
    ///////////////TESTOO
    //cerr << std::setprecision(12) << "AFTER maxFreq " << maxFreq << "  minFreq " << minFreq << endl;
    //for (Int k =0 ; k < (fselections_p->size()) ; ++k){
    //  cerr << k << (fselections_p->get(k)).toString() << endl;
    // }
    ///////////////////////////////////////	   
    ///TESTOO
    // andFreqSelection(-1, -1, minFreq, maxFreq, MFrequency::TOPO); 
    andFreqSelection(-1, -1, minFreq, maxFreq, selFreqFrame_p);
    
    vi_p->setFrequencySelection (*fselections_p);

    freqBegs_p=copyFreqBegs;
    freqEnds_p=copyFreqEnds;
    freqSpws_p=copyFreqSpws;
  }


Bool SynthesisImagerVi2::defineImage(
        SynthesisParamsImage& impars,
        const SynthesisParamsGrid& gridpars)
  {
    LogIO os( LogOrigin("SynthesisImagerVi2", "defineImage", WHERE) );

    if (mss_p.nelements() == 0)
      os << "SelectData has to be run before defineImage" << LogIO::EXCEPTION;

    CoordinateSystem csys;
    CountedPtr<refim::FTMachine> ftm, iftm;
    impars_p = impars;
    gridpars_p = gridpars; 

    try {
      os << "Define image coordinates for [" << impars.imageName << "] : "
         << LogIO::POST;
	os << "Define image coordinates for [" << impars.imageName << "] : " << LogIO::POST;
	//    cerr <<  "DEFIM " <<  gridpars_p.ftmachine <<  endl;
	//    cerr <<  "###### gridpars compute " <<  gridpars.computePAStep <<  "   " <<  gridpars_p.computePAStep <<  endl;
	csys = impars_p.buildCoordinateSystem( *vi_p, channelSelections_p, mss_p );
	//use the location defined for coordinates frame;
	mLocation_p=impars_p.obslocation;
	IPosition imshape = impars_p.shp();

	
         os << "Impars: start " << impars_p.start << LogIO::POST;
         os << "Shape: " << imshape
         << " Spectral: " << csys.spectralCoordinate().referenceValue()
         << " at " << csys.spectralCoordinate().referencePixel()
         << " with increment " << csys.spectralCoordinate().increment()
         << LogIO::POST;

	if( (itsMappers.nMappers()==0) || 
	    (impars_p.imsize[0]*impars_p.imsize[1] > itsMaxShape[0]*itsMaxShape[1]))
	  {
	    itsMaxShape=imshape;
	    itsMaxCoordSys=csys;
	  }
        itsNchan = imshape[3];
        itsCsysRec = impars_p.getcsys();
	/*
	os << "Define image  [" << impars.imageName << "] : nchan : " << impars.nchan 
	   //<< ", freqstart:" << impars.freqStart.getValue() << impars.freqStart.getUnit() 
	   << ", start:" << impars.start
	   <<  ", imsize:" << impars.imsize 
	   << ", cellsize: [" << impars.cellsize[0].getValue() << impars.cellsize[0].getUnit() 
	   << " , " << impars.cellsize[1].getValue() << impars.cellsize[1].getUnit() 
	   << LogIO::POST;
	*/
        // phasecenter
        if (impars_p.phaseCenterFieldId == -1) {
          // user-specified
          phaseCenter_p = impars_p.phaseCenter;
        } else if (impars_p.phaseCenterFieldId >= 0) {
          // FIELD_ID
          auto const msobj = mss_p[0];
          MSFieldColumns msfield(msobj->field());
          phaseCenter_p=msfield.phaseDirMeas(impars_p.phaseCenterFieldId);
        } else {
          // use default FIELD_ID (0)
          auto const msobj = mss_p[0];
          MSFieldColumns msfield(msobj->field());
          phaseCenter_p=msfield.phaseDirMeas(0);
        }

      }
    catch(AipsError &x)
      {
	os << "Error in building Coordinate System and Image Shape: " << x.getMesg() << LogIO::EXCEPTION;
      }

	
    try
      {
	os << "Set Gridding options for [" << impars_p.imageName << "] with ftmachine : " << gridpars.ftmachine << LogIO::POST;

	itsVpTable=gridpars.vpTable;

	itsMakeVP= ( gridpars.ftmachine.contains("mosaicft") ||
                     (gridpars.ftmachine.at(0,3)=="awp") )?False:True;
       
      createFTMachine(
        ftm, iftm, gridpars.ftmachine, impars_p.nTaylorTerms, gridpars.mType,
        gridpars.facets, gridpars.wprojplanes,
        gridpars.padding, gridpars.useAutoCorr, gridpars.useDoublePrec,
        gridpars.convFunc,
        gridpars.aTermOn, gridpars.psTermOn, gridpars.mTermOn,
        gridpars.wbAWP, gridpars.cfCache, gridpars.usePointing, gridpars.pointingOffsetSigDev.tovector(),
        gridpars.doPBCorr, gridpars.conjBeams,
        gridpars.computePAStep, gridpars.rotatePAStep,
        gridpars.interpolation, impars_p.freqFrameValid, 1000000000, 16, impars_p.stokes,
        impars_p.imageName,
        gridpars.pointingDirCol, gridpars.convertFirst, gridpars.skyPosThreshold,
        gridpars.convSupport, gridpars.truncateSize, gridpars.gwidth, gridpars.jwidth,
        gridpars.minWeight, gridpars.clipMinMax, impars_p.pseudoi
      );

      }
    catch (AipsError &x) {
      os << "Error in setting up FTMachine(): " << x.getMesg() << LogIO::EXCEPTION;
    }


    try {
      appendToMapperList(
        impars_p.imageName, csys, impars_p.shp(),
        ftm, iftm,
        gridpars.distance, gridpars.facets, gridpars.chanchunks, impars_p.overwrite,
        gridpars.mType, gridpars.padding, impars_p.nTaylorTerms, impars_p.startModel
      );

      imageDefined_p = true;
    }
    catch(AipsError &x) {
      os << "Error in adding Mapper: " + x.getMesg() << LogIO::EXCEPTION;
    }

    imparsVec_p.resize(imparsVec_p.nelements()+1, true);
    imparsVec_p[imparsVec_p.nelements()-1] = impars_p;
    ///For now cannot deal with cube and mtmfs in C++ parallel mode
    if (imparsVec_p[0].deconvolver == "mtmfs") setCubeGridding(False);
    // cerr << "DECONV " << imparsVec_p[0].deconvolver
    //      << " cube gridding " << doingCubeGridding_p << endl;
    gridparsVec_p.resize(gridparsVec_p.nelements()+1, true);
    gridparsVec_p[imparsVec_p.nelements()-1] = gridpars_p;
    // For now as awproject does not work with the c++ mpi cube gridding
    // make sure it works the old way as mfs
    // if ( gridparsVec_p[0].ftmachine.contains("awproject") )
    //   setCubeGridding(False);
    itsMakeVP =
      (  gridparsVec_p[0].ftmachine.contains("mosaicft") or
        (gridparsVec_p[0].ftmachine.at(0,3) == "awp")
      ) ? False : True;
    return true;
  }


Bool SynthesisImagerVi2::defineImage(
        CountedPtr<SIImageStore> imstor,
        SynthesisParamsImage& impars,
        const SynthesisParamsGrid& gridpars)
  {
    gridpars_p=gridpars;
    Int id = itsMappers.nMappers();
    CoordinateSystem csys = imstor->getCSys();
    IPosition imshape = imstor->getShape();
    Int nx=imshape[0], ny=imshape[1];
    if ( (id==0) || (nx*ny > itsMaxShape[0]*itsMaxShape[1]) ) {
      itsMaxShape=imshape;
      itsMaxCoordSys=csys;
    }
    mLocation_p=impars.obslocation;
    // phasecenter
    if (impars.phaseCenterFieldId == -1) { // user-specified
      phaseCenter_p = impars.phaseCenter;
    }
    else if (impars.phaseCenterFieldId >= 0) { // FIELD_ID
      auto const msobj = mss_p[0];
      MSFieldColumns msfield(msobj->field());
      phaseCenter_p = msfield.phaseDirMeas(impars.phaseCenterFieldId);
    }
    else { // use default FIELD_ID (0)
      auto const msobj = mss_p[0];
      MSFieldColumns msfield(msobj->field());
      phaseCenter_p = msfield.phaseDirMeas(0);
    }

    itsVpTable = gridpars.vpTable;
    itsMakeVP = (  gridpars.ftmachine.contains("mosaicft") or
                  (gridpars.ftmachine.at(0,3) == "awp")
                ) ? False : True;
    CountedPtr<refim::FTMachine> ftm, iftm;

    createFTMachine(
      ftm, iftm, gridpars.ftmachine, impars.nTaylorTerms, gridpars.mType,
      gridpars.facets, gridpars.wprojplanes,
      gridpars.padding,gridpars.useAutoCorr,gridpars.useDoublePrec,
      gridpars.convFunc,
      gridpars.aTermOn,gridpars.psTermOn, gridpars.mTermOn,
      gridpars.wbAWP,gridpars.cfCache,gridpars.usePointing,
      gridpars.pointingOffsetSigDev.tovector(),
      gridpars.doPBCorr,gridpars.conjBeams,
      gridpars.computePAStep,gridpars.rotatePAStep,
      gridpars.interpolation, impars.freqFrameValid, 1000000000, 16, impars.stokes,
      impars.imageName,
      gridpars.pointingDirCol, gridpars.convertFirst, gridpars.skyPosThreshold,
      gridpars.convSupport, gridpars.truncateSize, gridpars.gwidth, gridpars.jwidth,
      gridpars.minWeight, gridpars.clipMinMax, impars.pseudoi);

  if (gridpars.facets >1) {
    // Make and connect the list.
    Block<CountedPtr<SIImageStore> > imstorList =
      createFacetImageStoreList( imstor, gridpars.facets );
    for( uInt facet=0; facet<imstorList.nelements(); facet++) {
      CountedPtr<refim::FTMachine> new_ftm, new_iftm;
      if (facet == 0) {
        new_ftm = ftm;
        new_iftm = iftm;
      }
      else {
        new_ftm = ftm->cloneFTM();
        new_iftm = iftm->cloneFTM();
      }
      itsMappers.addMapper(
        createSIMapper( gridpars.mType, imstorList[facet], new_ftm, new_iftm)
      );
    }
  }
  else {
    itsMappers.addMapper(
      createSIMapper( gridpars.mType, imstor, ftm, iftm)
    );
  }
  impars_p = impars;
  gridpars_p = gridpars;
  imageDefined_p = true;

  imparsVec_p.resize(imparsVec_p.nelements()+1, true);
  imparsVec_p[imparsVec_p.nelements()-1] = impars_p;

  gridparsVec_p.resize(gridparsVec_p.nelements()+1, true);
  gridparsVec_p[gridparsVec_p.nelements()-1] = gridpars_p;

  return true;
}

Bool SynthesisImagerVi2::defineImage(CountedPtr<SIImageStore> imstor, 
				    const String& ftmachine)
  {
    CountedPtr<refim::FTMachine> ftm, iftm;

    // The following call to createFTMachine() uses the
    // following defaults
    //
    // facets=1, wprojplane=1, padding=1.0, useAutocorr=false, 
    // useDoublePrec=true, gridFunction=String("SF")
    //
    createFTMachine(ftm, iftm, ftmachine);
    
    Int id=itsMappers.nMappers();
    CoordinateSystem csys =imstor->getCSys();
    IPosition imshape=imstor->getShape();
    Int nx=imshape[0], ny=imshape[1];
    if( (id==0) || (nx*ny > itsMaxShape[0]*itsMaxShape[1]))
      {
	itsMaxShape=imshape;
	itsMaxCoordSys=csys;
      }

    itsMappers.addMapper(  createSIMapper( "default", imstor, ftm, iftm ) );
    imageDefined_p=true;
    return true;
  }
Bool SynthesisImagerVi2::defineImage(CountedPtr<SIImageStore> imstor, 
				    const Record& ftmRec, const Record& iftmRec)
  {
    CountedPtr<refim::FTMachine> ftm, iftm;
	ftm=refim::VisModelData::NEW_FT(ftmRec);
	iftm=refim::VisModelData::NEW_FT(iftmRec);
    
    Int id=itsMappers.nMappers();
    CoordinateSystem csys =imstor->getCSys();
    IPosition imshape=imstor->getShape();
    Int nx=imshape[0], ny=imshape[1];
    if( (id==0) || (nx*ny > itsMaxShape[0]*itsMaxShape[1]))
      {
	itsMaxShape=imshape;
	itsMaxCoordSys=csys;
      }

    itsMappers.addMapper(  createSIMapper( "default", imstor, ftm, iftm, id ) );
    imageDefined_p=true;
    return true;
  }
 Bool SynthesisImagerVi2::weight(const Record& inrec){
	String type, rmode, filtertype;
	Quantity noise, fieldofview,filterbmaj,filterbmin, filterbpa;  
	Double robust, fracBW;
	Int npixels;
	Bool multiField, useCubeBriggs;
	SynthesisUtilMethods::getFromWeightRecord(type, rmode,noise, robust,fieldofview,npixels, multiField, useCubeBriggs,
				  filtertype, filterbmaj,filterbmin, filterbpa, fracBW, inrec);
	return weight(type, rmode,noise, robust,fieldofview,npixels, multiField, useCubeBriggs,
				  filtertype, filterbmaj,filterbmin, filterbpa, fracBW );
				
	 
 }
 Bool SynthesisImagerVi2::weight(const String& type, const String& rmode,
			       const Quantity& noise, const Double robust,
			       const Quantity& fieldofview,
				 const Int npixels, const Bool multiField, const Bool useCubeBriggs,
			       const String& filtertype, const Quantity& filterbmaj,
			       const Quantity& filterbmin, const Quantity& filterbpa, Double fracBW)
  {
      LogIO os(LogOrigin("SynthesisImagerVi2", "weight()", WHERE));
      
      if(rmode=="bwtaper") //See CAS-13021 for bwtaper algorithm details
      {
          if(fracBW == 0.0)
          {
              Double minFreq = 0.0;
              Double maxFreq = 0.0;
              
              if(itsMaxShape(3) < 1) {
                cout << "SynthesisImagerVi2::weight Only one channel in image " << endl;
              }
              else{
                  minFreq=abs(SpectralImageUtil::worldFreq(itsMaxCoordSys, 0.0));
                  maxFreq=abs(SpectralImageUtil::worldFreq(itsMaxCoordSys,itsMaxShape(3)-1));

                  if(maxFreq < minFreq){
                    Double tmp=minFreq;
                    minFreq=maxFreq;
                    maxFreq=tmp;
                  }
                  
                  if((maxFreq != 0.0) || (minFreq != 0.0)) fracBW = 2*(maxFreq - minFreq)/(maxFreq + minFreq);
                  
                  os << LogIO::NORMAL << " Fractional bandwidth used by briggsbwtaper " << fracBW << endl;  //<< LogIO::POST;
                  
              }
          }
      }
      
	weightParams_p=SynthesisUtilMethods::fillWeightRecord(type, rmode,noise, robust,fieldofview,
				 npixels, multiField, useCubeBriggs,filtertype, filterbmaj,filterbmin, filterbpa, fracBW);

       try {
    	//Int nx=itsMaxShape[0];
    	//Int ny=itsMaxShape[1];
        

         ///////////////////////
	 Quantity cellx=Quantity(itsMaxCoordSys.increment()[0], itsMaxCoordSys.worldAxisUnits()[0]);
	 Quantity celly=Quantity(itsMaxCoordSys.increment()[1], itsMaxCoordSys.worldAxisUnits()[1]);
	 os << LogIO::NORMAL // Loglevel INFO
	    << "Set imaging weights : " ; //<< LogIO::POST;
	 
	 if (type=="natural") {
	   os << LogIO::NORMAL // Loglevel INFO
	      << "Natural weighting" << LogIO::POST;
	   imwgt_p=VisImagingWeight("natural");
	 }
      else if (type=="radial") {
	os << "Radial weighting" << LogIO::POST;
    	  imwgt_p=VisImagingWeight("radial");
      }
      else{
	vi_p->originChunks();
	vi_p->origin();
    	  if(!imageDefined_p)
    		  throw(AipsError("Need to define image"));
    	  Int nx=itsMaxShape[0];
    	  Int ny=itsMaxShape[1];
    	  Quantity cellx=Quantity(itsMaxCoordSys.increment()[0], itsMaxCoordSys.worldAxisUnits()[0]);
    	  Quantity celly=Quantity(itsMaxCoordSys.increment()[1], itsMaxCoordSys.worldAxisUnits()[1]);
    	  if(type=="superuniform"){
    		  if(!imageDefined_p) throw(AipsError("Please define image first"));
    		  Int actualNpix=npixels;
    		  if(actualNpix <=0)
    			  actualNpix=3;
    		  os << LogIO::NORMAL // Loglevel INFO
    				  << "SuperUniform weighting over a square cell spanning ["
    				  << -actualNpix
    				  << ", " << actualNpix << "] in the uv plane" << LogIO::POST;
    		  imwgt_p=VisImagingWeight(*vi_p, rmode, noise, robust, nx,
    				  ny, cellx, celly, actualNpix,
    				  actualNpix, multiField);
    	  }
    	  else if ((type=="robust")||(type=="uniform")||(type=="briggs")) {
   		  if(!imageDefined_p) throw(AipsError("Please define image first"));
    		  Quantity actualFieldOfView_x(fieldofview), actualFieldOfView_y(fieldofview) ;
    		  Int actualNPixels_x(npixels),actualNPixels_y(npixels) ;
    		  String wtype;
    		  if(type=="briggs") {
    			  wtype = "Briggs";
    		  }
    		  else {
    			  wtype = "Uniform";
    		  }
    		  if(actualFieldOfView_x.get().getValue()==0.0&&actualNPixels_x==0) {
    			  actualNPixels_x=nx;
    			  actualFieldOfView_x=Quantity(actualNPixels_x*cellx.get("rad").getValue(),"rad");
    			  actualNPixels_y=ny;
    			  actualFieldOfView_y=Quantity(actualNPixels_y*celly.get("rad").getValue(),"rad");
    			  os << LogIO::NORMAL // Loglevel INFO
    					  << wtype
    					  << " weighting: sidelobes will be suppressed over full image"
    					  << LogIO::POST;
    		  }
    		  else if(actualFieldOfView_x.get().getValue()>0.0&&actualNPixels_x==0) {
    			  actualNPixels_x=nx;
    			  actualNPixels_y=ny;
    			  os << LogIO::NORMAL // Loglevel INFO
    					  << wtype
    					  << " weighting: sidelobes will be suppressed over specified field of view: "
			                  << actualFieldOfView_x.get("arcsec").getValue() << " arcsec by " 
			                  << actualFieldOfView_y.get("arcsec").getValue()  << " arcsec" << LogIO::POST;
    		  }
    		  else if(actualFieldOfView_x.get().getValue()==0.0&&actualNPixels_x>0) {
    			  actualFieldOfView_x=Quantity(actualNPixels_x*cellx.get("rad").getValue(),"rad");
    			  actualFieldOfView_y=Quantity(actualNPixels_y*celly.get("rad").getValue(),"rad");
    			  os << LogIO::NORMAL // Loglevel INFO
    					  << wtype
    					  << " weighting: sidelobes will be suppressed over full image field of view: "
			                  << actualFieldOfView_x.get("arcsec").getValue() << " arcsec by " 
    					  << actualFieldOfView_y.get("arcsec").getValue() << " arcsec" << LogIO::POST;
    		  }
    		  else {
    			  os << LogIO::NORMAL // Loglevel INFO
    					  << wtype
    					  << " weighting: sidelobes will be suppressed over specified field of view: "
			                  << actualFieldOfView_x.get("arcsec").getValue() << " arcsec by " 
    					  << actualFieldOfView_y.get("arcsec").getValue() << " arcsec" << LogIO::POST;
    		  }
    		  os << LogIO::DEBUG1
		     << "Weighting used " << actualNPixels_x << " by " << actualNPixels_y << " uv pixels."
		     << LogIO::POST;
    		  Quantity actualCellSize_x(actualFieldOfView_x.get("rad").getValue()/actualNPixels_x, "rad");
    		  Quantity actualCellSize_y(actualFieldOfView_y.get("rad").getValue()/actualNPixels_y, "rad");


		  //		  cerr << "rmode " << rmode << " noise " << noise << " robust " << robust << " npixels " << actualNPixels << " cellsize " << actualCellSize << " multifield " << multiField << endl;
		  //		  Timer timer;
		  //timer.mark();
		  //Construct imwgt_p with old vi for now if old vi is in use as constructing with vi2 is slower
		  //Determine if any image is cube
		  if(useCubeBriggs){
		    String outstr=String("Doing spectral cube Briggs weighting formula --  " + rmode + (rmode=="abs" ? " with estimated noise "+ String::toString(noise.getValue())+noise.getUnit()  : "")); 
		    os << outstr << LogIO::POST;
		    //VisImagingWeight nat("natural");
		    //vi_p->useImagingWeight(nat);
		    if(rmode=="abs" && robust==0.0 && noise.getValue()==0.0)
		      throw(AipsError("Absolute Briggs formula does not allow for robust 0 and estimated noise per visibility 0"));
                
            CountedPtr<refim::BriggsCubeWeightor> bwgt=new refim::BriggsCubeWeightor(wtype=="Uniform" ? "none" : rmode, noise, robust,fracBW, npixels, multiField);
            for (Int k=0; k < itsMappers.nMappers(); ++k){
                      itsMappers.getFTM2(k)->setBriggsCubeWeight(bwgt);
                    }
              
              
		  }
		  else
		  {
		    imwgt_p=VisImagingWeight(*vi_p, wtype=="Uniform" ? "none" : rmode, noise, robust,
					     actualNPixels_x, actualNPixels_y, actualCellSize_x,
					     actualCellSize_y, 0, 0, multiField);
		  }

		  /*
		  if(rvi_p !=NULL){
		    imwgt_p=VisImagingWeight(*rvi_p, rmode, noise, robust,
                                 actualNPixels, actualNPixels, actualCellSize,
                                 actualCellSize, 0, 0, multiField);
		  }
		  else{
		    ////This is slower by orders of magnitude as of 2014/06/25
		    imwgt_p=VisImagingWeight(*vi_p, rmode, noise, robust,
                                 actualNPixels, actualNPixels, actualCellSize,
                                 actualCellSize, 0, 0, multiField);
		  }
		  */
		    //timer.show("After making visweight ");

    	  }
    	  else {
    		  //this->unlock();
    		  os << LogIO::SEVERE << "Unknown weighting " << type
    				  << LogIO::EXCEPTION;
    		  return false;
    	  }
      }
	 
	 //// UV-Tapering
	 //cout << "Taper type : " << filtertype << " : " << (filtertype=="gaussian") <<  endl;
	 if( filtertype == "gaussian" ) {
	   //	   os << "Setting uv-taper" << LogIO::POST;
	   imwgt_p.setFilter( filtertype,  filterbmaj, filterbmin, filterbpa );
	 }
	 vi_p->useImagingWeight(imwgt_p);
      ///////////////////////////////
	 
	     SynthesisUtilMethods::getResource("Set Weighting");

	 ///	 return true;
	 
       }
       catch(AipsError &x)
	 {
	   throw( AipsError("Error in Weighting : "+x.getMesg()) );
	 }
       
       return true;
  }

void SynthesisImagerVi2::appendToMapperList(String imagename,  
					   CoordinateSystem& csys, 
					   IPosition imshape,
					    CountedPtr<refim::FTMachine>& ftm,
					    CountedPtr<refim::FTMachine>& iftm,
					   Quantity distance,
					   Int facets,
					   Int chanchunks,
					   const Bool overwrite,
					   String mappertype,
					   Float padding,
					   uInt ntaylorterms,
					   const Vector<String> &startmodel)
    {
      LogIO log_l(LogOrigin("SynthesisImagerVi2", "appendToMapperList(ftm)"));
      //---------------------------------------------
      // Some checks..
      
      if(facets > 1 && itsMappers.nMappers() > 0)
	log_l << "Facetted image has to be the first of multifields" << LogIO::EXCEPTION;

     TcleanProcessingInfo procInfo;
     CompositeNumber cn(uInt(imshape[0] * 2));
     // heuristic factors multiplied to imshape based on gridder
     size_t fudge_factor = 15;
     if (ftm->name()=="MosaicFTNew") {
         fudge_factor = 20;
     }
     else if (ftm->name()=="GridFT") {
         fudge_factor = 9;
     }
     std::tie(procInfo, std::ignore, std::ignore) =
         nSubCubeFitInMemory(fudge_factor, imshape, padding);

     // chanchunks auto-calculation block, for now still here for awproject (CAS-12204)
     if(chanchunks<1)
	{
	  log_l << "Automatically calculated chanchunks";
	  log_l << " using imshape : " << imshape << LogIO::POST;

	  // Do calculation here.
	  // This runs once per image field (for multi-field imaging)
	  // This runs once per cube partition, and will see only its own partition's shape
		//chanchunks=1;

                chanchunks = procInfo.chnchnks;

		/*log_l << "Required memory " << required_mem / nlocal_procs / 1024. / 1024. / 1024.
                 << "\nAvailable memory " << memory_avail / 1024. / 1024 / 1024.
                 << " (rc: memory fraction " << usr_memfrac << "% rc memory " << usr_mem / 1024.
                 << ")\n" << nlocal_procs << " other processes on node\n"
                 << "Setting chanchunks to " << chanchunks << LogIO::POST;
		*/
	}
	//record this in gridpars_p
	gridpars_p.chanchunks=chanchunks;
      if( imshape.nelements()==4 && imshape[3]<chanchunks )
	{
	  log_l << LogIO::WARN << "An image with " << imshape[3] << " channel(s) cannot be divided into " << chanchunks << " chunks. Please set chanchunks=1 or choose chanchunks<nchan." << LogIO::EXCEPTION;
	}

      if(chanchunks > 1 && itsMappers.nMappers() > 0)
	log_l << "Channel chunking is currently not supported with multi(outlier)-fields. Please submit a feature request if needed." << LogIO::EXCEPTION;

      if(chanchunks > 1) itsDataLoopPerMapper=true;
      
      AlwaysAssert( ( ( ! (ftm->name()=="MosaicFTNew" && mappertype=="imagemosaic") )  && 
		      ( ! (ftm->name()=="AWProjectWBFT" && mappertype=="imagemosaic") )) ,
		    AipsError );
      //---------------------------------------------

      // Create the ImageStore object
      CountedPtr<SIImageStore> imstor;
      MSColumns msc(*(mss_p[0]));
      imstor = createIMStore(imagename, csys, imshape, overwrite,msc, mappertype, ntaylorterms, distance, procInfo, facets, iftm->useWeightImage(), startmodel );

      // Create the Mappers
      if( facets<2 && chanchunks<2) // One facet. Just add the above imagestore to the mapper list.
	{
	  itsMappers.addMapper(  createSIMapper( mappertype, imstor, ftm, iftm, ntaylorterms) );
	}
      else // This field is facetted. Make a list of reference imstores, and add all to the mapper list.
	{

	  if ( facets>1 && chanchunks==1 )
	    {
	      // Make and connect the list.
	      Block<CountedPtr<SIImageStore> > imstorList = createFacetImageStoreList( imstor, facets );
	      for( uInt facet=0; facet<imstorList.nelements(); facet++)
		{
		  CountedPtr<refim::FTMachine> new_ftm, new_iftm;
		  if(facet==0){ new_ftm = ftm;  new_iftm = iftm; }
		  else{ new_ftm=ftm->cloneFTM();  new_iftm=iftm->cloneFTM(); }
// 		  imstorList[facet]->setDataPolFrame(imstor->getDataPolFrame());
		  itsMappers.addMapper(createSIMapper( mappertype, imstorList[facet], new_ftm, new_iftm, ntaylorterms));
		}
	    }// facets
	  else if ( facets==1 && chanchunks>1 )
	    {
	      // Make and connect the list.
	      Block<CountedPtr<SIImageStore> > imstorList = createChanChunkImageStoreList( imstor, chanchunks );
	      for( uInt chunk=0; chunk<imstorList.nelements(); chunk++)
		{
		  
		  CountedPtr<refim::FTMachine> new_ftm, new_iftm;
		  if(chunk==0){ 
		    new_ftm = ftm;  
		    new_iftm = iftm; }
		  else{ 
		    new_ftm=ftm->cloneFTM();  
		    new_iftm=iftm->cloneFTM(); }
		  imstorList[chunk]->setDataPolFrame(imstor->getDataPolFrame());
		  itsMappers.addMapper(createSIMapper( mappertype, imstorList[chunk], new_ftm, new_iftm, ntaylorterms));
		}
	    }// chanchunks
	  else
	    {
	      throw( AipsError("Error in requesting "+String::toString(facets)+" facets on a side with " + String::toString(chanchunks) + " channel chunks.  Support for faceting along with channel chunking is not yet available. Please submit a feature-request if you need multiple facets as well as chanchunks. ") );
	    }

	}// facets or chunks

    }

  /////////////////////////
  /**
   * Calculations of memory required / available -> nchunks .
   *
   * Returns a tuple with a TcleanProcessingInfo, vector of start channels per subchunk,
   * vector of end channels.
   */
  std::tuple<TcleanProcessingInfo, Vector<Int>, Vector<Int> > SynthesisImagerVi2::nSubCubeFitInMemory(const Int fudge_factor, const IPosition& imshape, const Float padding){
	LogIO log_l(LogOrigin("SynthesisImagerVi2", "nSubCubeFitInMemory"));

	Double required_mem = fudge_factor * sizeof(Float);
	int nsubcube=1;
	CompositeNumber cn(uInt(imshape[0] * 2));
	for (size_t i = 0; i < imshape.nelements(); i++) {
			// gridding pads image and increases to composite number
			if (i < 2) {
				required_mem *= cn.nextLargerEven(Int(padding*Float(imshape[i])-0.5));
			}
			else {
				required_mem *= imshape[i];
			}
	}

	// get number of tclean processes running on the same machine
	size_t nlocal_procs = 1;
	if (getenv("OMPI_COMM_WORLD_LOCAL_SIZE")) {
		std::stringstream ss(getenv("OMPI_COMM_WORLD_LOCAL_SIZE"));
		ss >> nlocal_procs;
	}
        //cerr << "NUM_PROC " << nlocal_procs << endl;
	// assumes all processes need the same amount of memory
	required_mem *= nlocal_procs;
	Double usr_memfrac, usr_mem;
	AipsrcValue<Double>::find(usr_memfrac, "system.resources.memfrac", 80.);
	AipsrcValue<Double>::find(usr_mem, "system.resources.memory", -1.);
	Double memory_avail;
	if (usr_mem > 0.) {
		memory_avail = usr_mem * 1024. * 1024.;
	}
	else {
	    memory_avail = Double(HostInfo::memoryFree()) * (usr_memfrac / 100.) * 1024.;
	}
	// compute required chanchunks to fit into the available memory
	nsubcube = (int)std::ceil((Double)required_mem / memory_avail);
        Int nworkers= applicator.numProcs() < 2 ? 1 : applicator.numProcs()-1;
        if((nsubcube/nworkers) >1 && nworkers !=1){
          nsubcube=(Int(std::floor(Float(nsubcube)/Float(nworkers)))+1)*nworkers;

        }
	if (imshape.nelements() == 4 && imshape[3] < nsubcube) {
		nsubcube = imshape[3];
              // TODO make chanchunks a divisor of nchannels?
	}
	nsubcube = nsubcube < 1 ? 1 : nsubcube;

        
        if( (imshape[3] >= nworkers) && (nsubcube < nworkers)){
          nsubcube=nworkers;
          ///TESTOO
          //if(imshape[3] > 2*nworkers)
          //  nsubcube=2*nworkers;

        }
        else if(imshape[3] < (applicator.numProcs()-1)){
		nsubcube=imshape[3]; 
        }
        cerr << "FTM name " << gridpars_p.ftmachine << endl;
        if(gridpars_p.ftmachine.at(0,4)=="awph")
          nsubcube=imshape[3];
        cerr << "SUBCUBE " << nsubcube << " shp  " << imshape[3] << endl;
        
	Int chunksize=imshape[3]/nsubcube;
	Int rem=imshape[3] % nsubcube;
	//case of nchan < numprocs
	if(chunksize==0 && rem > 0){
		nsubcube=rem;
		chunksize=1;
		rem=0;
	}
	///Avoid an extra chunk with 1 channel as it cause bumps in linear interpolation
	///See CAS-12625
        /*
	while((rem==1) && (chunksize >1)){
		nsubcube +=1;
		chunksize=imshape[3]/nsubcube;
		rem=imshape[3] % nsubcube;
	}
	if(rem !=0) ++nsubcube;
        . */
	Vector<Int> start(nsubcube,0);
	Vector<Int> end(nsubcube,chunksize-1);
        if(rem >0){
          end(0)+=1;
          --rem;
        }
	for (Int k=1; k < nsubcube; ++k){
		start(k)=end(k-1)+1;
                //	end(k)=((k !=nsubcube-1) || rem==0)? (start(k)+chunksize-1) : (start(k)+rem-1);
                end(k)=(start(k)+chunksize-1);
                if(rem > 0){
                  end(k)+=1;
                  --rem;
                }
	}

        // print mem related info to log
        const float toGB = 1024.0 * 1024.0 * 1024.0;
        std::ostringstream usr_mem_msg;
        if (usr_mem > 0.) {
            usr_mem_msg << usr_mem / 1024.;
        } else {
            usr_mem_msg << "-";
        }
        std::ostringstream oss;
        oss << setprecision(4);
        oss << "Required memory: " << required_mem / toGB
            << " GB. Available mem.: " << memory_avail / toGB
            << " GB (rc, mem. fraction: " << usr_memfrac
            << "%, memory: " << usr_mem_msg.str()
            << ") => Subcubes: " << nsubcube
            << ". Processes on node: " << nlocal_procs << ".\n";
        log_l << oss.str() << LogIO::POST;

        TcleanProcessingInfo procInfo;
        procInfo.mpiprocs = nlocal_procs;
        procInfo.chnchnks = nsubcube;
        procInfo.memavail = memory_avail / toGB;
        procInfo.memreq = required_mem / toGB;
        return make_tuple(procInfo, start, end);
  }
  
 void SynthesisImagerVi2::runMajorCycleCube( const Bool dopsf, 
				      const Record inpcontrol) {
	LogIO os( LogOrigin("SynthesisImagerVi2","runMajorCycleCube",WHERE) );		  
	if(dopsf){
	  runCubeGridding(True);
	  ///Store the beamsets in ImageInfo
	  for(Int k=0; k < itsMappers.nMappers(); ++k){
	   
	    (itsMappers.imageStore(k))->getBeamSet();
	  }
	}
	else
		runCubeGridding(False, inpcontrol);
	
			  
			  
	}
 void SynthesisImagerVi2::runMajorCycle(const Bool dopsf, 
				      const Bool savemodel)
  {
    LogIO os( LogOrigin("SynthesisImagerVi2","runMajorCycle",WHERE) );

    //    cout << "Savemodel : " << savemodel << "   readonly : " << readOnly_p << "   usescratch : " << useScratch_p << endl;
    Bool savemodelcolumn = savemodel && !readOnly_p && useScratch_p;
    Bool savevirtualmodel = savemodel && !readOnly_p && !useScratch_p;

    if( savemodelcolumn ) os << "Saving model column" << LogIO::POST;
    if( savevirtualmodel ) os << "Saving virtual model" << LogIO::POST;

    SynthesisUtilMethods::getResource("Start Major Cycle");

    itsMappers.checkOverlappingModels("blank");

    {
      vi::VisBuffer2* vb=vi_p->getVisBuffer();
      vi_p->originChunks();
      vi_p->origin();
      /////////////////////////////////////
      /////////////////
      if(gridparsVec_p[0].ftmachine=="awphpg"){
        //reset the ftmachine as it keeps dying on second run
        resetAWPHPG();
      }
      //////////////////////////////////
      Double numcoh=0;
      for (uInt k=0; k< mss_p.nelements(); ++k)
	numcoh+=Double(mss_p[k]->nrow());
      ProgressMeter pm(1.0, numcoh, 
			 dopsf?"Gridding Weights and PSF":"Major Cycle", "","","",true);
      rownr_t cohDone=0;


    	if(!dopsf)itsMappers.initializeDegrid(*vb);
    	itsMappers.initializeGrid(*vi_p,dopsf);
	SynthesisUtilMethods::getResource("After initGrid for all mappers");
        ////Under some peculiar selection criterion and low channel ms  vb2 seems to return more channels than in spw
        
        if (gridparsVec_p[0].ftmachine.at(0,3) != "awp")
        {
          vi_p->originChunks();
          vi_p->origin();
          Int nchannow=vb->nChannels();
          Int spwnow=vb->spectralWindows()[0];
          Int nchaninms=MSColumns(vb->ms()).spectralWindow().numChan()(spwnow);
          //cerr << "chans " << nchaninms << "   " << nchannow << endl;
         
          if (nchaninms < nchannow){
            cerr << "NCHANS ms" << nchaninms << " now " << nchannow << " spw " << spwnow << "   " << vb->spectralWindows()[0] << endl;
            throw(AipsError("A nasty Visbuffer2 error occured...wait for CNGI"));
          }
        }
          //////
    	for (vi_p->originChunks(); vi_p->moreChunks();vi_p->nextChunk())
    	{

	  for (vi_p->origin(); vi_p->more(); vi_p->next())
    		{
		  //if (SynthesisUtilMethods::validate(*vb)==SynthesisUtilMethods::NOVALIDROWS) break; // No valid rows in this VB
		  //		  cerr << "nRows "<< vb->nRow() << "   " << max(vb->visCube()) <<  endl;
		  if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
		    {
    			if(!dopsf) {
			  { Cube<Complex> mod(vb->nCorrelations(), vb->nChannels(), vb->nRows(), Complex(0.0));
			    vb->setVisCubeModel(mod); 
			  }
                          if(gridparsVec_p[0].ftmachine !="awphpg"){
                            itsMappers.degrid(*vb, savevirtualmodel );
                            if(savemodelcolumn && writeAccess_p ){	
                              const_cast<MeasurementSet& >((vi_p->ms())).lock(true);
                              vi_p->writeVisModel(vb->visCubeModel());
                              const_cast<MeasurementSet& >((vi_p->ms())).unlock();
                            }
			    //static_cast<VisibilityIteratorImpl2 *> (vi_p->getImpl())->writeVisModel(vb->visCubeModel());

			    // Cube<Complex> tt=vb->visCubeModel();
			    // tt = 20.0;
			    // cout << "Vis:" << tt << endl;
			    // static_cast<VisibilityIteratorImpl2 *> (vi_p->getImpl())->writeVisModel(tt);
			  }
    			}
    			itsMappers.grid(*vb, dopsf, (refim::FTMachine::Type)datacol_p);
			
			cohDone += vb->nRows();
			pm.update(Double(cohDone));
		    }
    		}
    	}

	// cerr << "VI2 data: " << cohDone << endl;
	// exit(0);
    	//cerr << "IN SYNTHE_IMA" << endl;
    	//VisModelData::listModel(rvi_p->getMeasurementSet());
	SynthesisUtilMethods::getResource("Before finalize for all mappers");
    	if(!dopsf) itsMappers.finalizeDegrid(*vb);
    	itsMappers.finalizeGrid(*vb, dopsf);

    }

    itsMappers.checkOverlappingModels("restore");

    unlockMSs();

    SynthesisUtilMethods::getResource("End Major Cycle");

  }// end runMajorCycle



  void SynthesisImagerVi2::resetAWPHPG(){
    CountedPtr<refim::FTMachine> theFT=nullptr;
    CountedPtr<refim::FTMachine> theIFT=nullptr;
    SynthesisParamsGrid *gp= &(gridparsVec_p[0]);
    createAWPFTMachine(theFT, theIFT, gp->ftmachine, gp->facets,
                       gp->wprojplanes, 
                       gp->padding, gp->useAutoCorr, gp->useDoublePrec,
                       gp->convFunc, //gridFunction,
                       gp->aTermOn, gp->psTermOn, gp->mTermOn, gp->wbAWP,
                       gp->cfCache, gp->usePointing, (gp->pointingOffsetSigDev).tovector(),
                       gp->doPBCorr, gp->conjBeams, gp->computePAStep,
                       gp->rotatePAStep, 1000000000,16,imparsVec_p[0].imageName);
    itsMappers.setFTM2(0, theIFT, true);
    itsMappers.setFTM2(0, theFT, false);
    loadMosaicSensitivity();

  }
  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

  /// The mapper loop is outside the data iterator loop.
  /// This is for cases where the image size is large compared to the RAM and
  /// where data I/O is the relatively minor cost.
  void SynthesisImagerVi2::runMajorCycle2(const Bool dopsf, 
				      const Bool savemodel)
  {
    LogIO os( LogOrigin("SynthesisImagerVi2","runMajorCycle2",WHERE) );

    //    cout << "Savemodel : " << savemodel << "   readonly : " << readOnly_p << "   usescratch : " << useScratch_p << endl;

    Bool savemodelcolumn = savemodel && !readOnly_p && useScratch_p;
    Bool savevirtualmodel = savemodel && !readOnly_p && !useScratch_p;

    if( savemodelcolumn ) os << "Saving model column" << LogIO::POST;
    if( savevirtualmodel ) os << "Saving virtual model" << LogIO::POST;

    itsMappers.checkOverlappingModels("blank");

    Bool resetModel=False;
    if( savemodelcolumn && writeAccess_p)
      {
	resetModel=True;
	os << "Iterating through the model column to reset it to zero" << LogIO::POST;
	vi::VisBuffer2* vb=vi_p->getVisBuffer();
    	vi_p->originChunks();
    	vi_p->origin();
	Double numcoh=0;
	for (uInt k=0; k< mss_p.nelements(); ++k)
	  numcoh+=Double(mss_p[k]->nrow());
	ProgressMeter pm(1.0, numcoh, 
			 dopsf?"Seting model column to zero":"pre-Major Cycle", "","","",True);
        rownr_t cohDone=0;
    	for (vi_p->originChunks(); vi_p->moreChunks();vi_p->nextChunk())
	  {
	    
	    for (vi_p->origin(); vi_p->more(); vi_p->next())
	      {
		if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
		  {
		    { Cube<Complex> mod(vb->nCorrelations(), vb->nChannels(), vb->nRows(), Complex(0.0));
			    vb->setVisCubeModel(mod); 
		    }
		    vi_p->writeVisModel(vb->visCubeModel());
		    
		  }
		cohDone += vb->nRows();;
		pm.update(Double(cohDone));
	      }
	  }
      }// setting model to zero

    //Need to inialize the the forward ft machine to save the virtual model on first pass of each ms.
    if(!dopsf && savevirtualmodel){
      vi::VisBuffer2* vb=vi_p->getVisBuffer();
      vi_p->originChunks();
      vi_p->origin();
      itsMappers.initializeDegrid(*vb, -1);
    }
    
    for(Int gmap=0;gmap<itsMappers.nMappers();gmap++)
       {
	 os << "Running major cycle for chunk : " << gmap << LogIO::POST;

	 SynthesisUtilMethods::getResource("Start Major Cycle for mapper"+String::toString(gmap));
	 CountedPtr<vi::FrequencySelections> copyFsels=fselections_p->clone();
	 ///CAS-12132  create a new visiter for each chunk
	 createVisSet(writeAccess_p);
	 ////////////////////////
	 vi::VisBuffer2* vb=vi_p->getVisBuffer();
	 /// Careful where tunechunk 
	 tuneChunk(gmap);

	 vi_p->originChunks();
	 vi_p->origin();

	 Double numcoh=0;
	 for (uInt k=0; k< mss_p.nelements(); ++k)
	   numcoh+=Double(mss_p[k]->nrow());


	 ProgressMeter pm(1.0, numcoh, 
			  dopsf?"Gridding Weights and PSF":"Major Cycle", "","","",true);
	rownr_t cohDone=0;


	itsMappers.getFTM2(gmap, False)->reset();
	itsMappers.getFTM2(gmap, True)->reset();

    	if(!dopsf){
	  itsMappers.initializeDegrid(*vb, gmap);
		  //itsMappers.getMapper(gmap)->initializeDegrid(*vb);
	}
	itsMappers.initializeGrid(*vi_p,dopsf, gmap);
		//itsMappers.getMapper(gmap)->initializeGrid(*vb,dopsf);

	SynthesisUtilMethods::getResource("After initialize for mapper"+String::toString(gmap));
	Int iterNum=0;

	
    	for (vi_p->originChunks(); vi_p->moreChunks();vi_p->nextChunk())
    	{

	  for (vi_p->origin(); vi_p->more(); vi_p->next())
	    {
	      //if (SynthesisUtilMethods::validate(*vb)==SynthesisUtilMethods::NOVALIDROWS) break; // No valid rows in this VB
	      //		  cerr << "nRows "<< vb->nRow() << "   " << max(vb->visCube()) <<  endl;
	      if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
		{
		  
		  if(!dopsf) {
		    if(resetModel==False) 
		      { 
			Cube<Complex> mod(vb->nCorrelations(), vb->nChannels(), vb->nRows(), Complex(0.0));
			vb->setVisCubeModel(mod); 
		      }
		    itsMappers.degrid(*vb, savevirtualmodel, gmap );
		    //itsMappers.getMapper(gmap)->degrid(*vb); //, savevirtualmodel );
		    if(savemodelcolumn && writeAccess_p ){
		      vi_p->writeVisModel(vb->visCubeModel());
		      //vi_p->writeBackChanges(vb);
		      // static_cast<VisibilityIteratorImpl2 *> (vi_p->getImpl())->writeVisModel(vb->visCubeModel());
		    }

		  }
		  itsMappers.grid(*vb, dopsf, (refim::FTMachine::Type)(datacol_p), gmap);
		  //itsMappers.getMapper(gmap)->grid(*vb, dopsf, datacol_p);
		  cohDone += vb->nRows();
		  ++iterNum;
		  pm.update(Double(cohDone));
		}
	    }
    	}
    	//cerr << "IN SYNTHE_IMA" << endl;
    	//VisModelData::listModel(rvi_p->getMeasurementSet());

	SynthesisUtilMethods::getResource("Before finalize for mapper"+String::toString(gmap));
	
    	if(!dopsf) 
	  {
	    itsMappers.finalizeDegrid(*vb,gmap);
	    //itsMappers.getMapper(gmap)->finalizeDegrid();
	  }
	itsMappers.finalizeGrid(*vb, dopsf,gmap);
    	//itsMappers.getMapper(gmap)->finalizeGrid(*vb, dopsf);
	
	//	itsMappers.getMapper(gmap)->releaseImageLocks();
	itsMappers.getMapper(gmap)->imageStore()->releaseComplexGrids();        
	
	SynthesisUtilMethods::getResource("End Major Cycle for mapper"+String::toString(gmap));
	fselections_p=copyFsels;
       }// end of mapper loop
    ///CAS-12132  create a new visiter for each chunk
    createVisSet(writeAccess_p);
    ////////////////////////
    //////vi_p->setFrequencySelection(*fselections_p);

    itsMappers.checkOverlappingModels("restore");

    unlockMSs();

    SynthesisUtilMethods::getResource("End Major Cycle");

  }// end runMajorCycle2

//////////////////////////////
 
  ////////////////////////////////////////////////////////////////////////////////////////////////

  bool SynthesisImagerVi2::runCubeGridding(Bool dopsf, const Record inpcontrol){

	LogIO logger(LogOrigin("SynthesisImagerVi2", "runCubeGridding", WHERE));

	  //dummy for now as init is overloaded on this signature
        int argc=1;
        char **argv=nullptr;
        casa::applicator.init ( argc, argv );
	  //For serial or master only
	if ( applicator.isController() ) {
		CubeMajorCycleAlgorithm cmc;
		//casa::applicator.defineAlgorithm(cmc);
		//Initialize everything to get the setup as serial
		{
                
			vi_p->originChunks();
			vi_p->origin();
		
		}
		///Break things into chunks for parallelization or memory availabbility
		Vector<Int> startchan;
		Vector<Int> endchan; 
		Int numchunks;
		Int fudge_factor = 15;
		if ((itsMappers.getFTM2(0))->name()=="MosaicFTNew") {
			fudge_factor = 20;
		}
		else if ((itsMappers.getFTM2(0))->name()=="GridFT") {
			fudge_factor = 9;
		}
                TcleanProcessingInfo procInfo;
		std::tie(procInfo, startchan, endchan)=nSubCubeFitInMemory(fudge_factor, itsMaxShape, gridpars_p.padding);
                numchunks = procInfo.chnchnks;
		////TESTOO
		//numchunks=2;
		//startchan.resize(2);startchan[0]=0; startchan[1]=2;
		//endchan.resize(2); endchan[0]=1; endchan[1]=2;
		
		/////END TESTOO
		//cerr << "NUMCHUNKS " << numchunks << " start " <<  startchan << " end " << endchan << endl;
		Record controlRecord=inpcontrol;
		//For now just field 0 but should loop over all
		///This is to pass in explicit model, residual names etc
		controlRecord.define("nfields", Int(imparsVec_p.nelements()));
        //CountedPtr<SIImageStore> imstor = imageStore ( 0 );
        // checking that psf,  residual and sumwt is allDone
        //cerr << "shapes "  <<  imstor->residual()->shape() <<  " " <<  imstor->sumwt()->shape() <<  endl;
		if(!dopsf){
		        
		  //controlRecord.define("lastcycle",  savemodel);
		  controlRecord.define("nmajorcycles", nMajorCycles);
			Vector<String> modelnames(Int(imparsVec_p.nelements()),"");
			for(uInt k=0; k < imparsVec_p.nelements(); ++k){
				Int imageStoreId=k;
				if(k>0){
					if(gridparsVec_p[0].chanchunks >1 && uInt(itsMappers.nMappers()) >=    (gridparsVec_p[0].chanchunks + gridparsVec_p.nelements()))
						imageStoreId+=gridparsVec_p[0].chanchunks-1;
					if(gridparsVec_p[0].facets >1)
						imageStoreId+=gridparsVec_p[0].facets*gridparsVec_p[0].facets-1;
				}
				if((itsMappers.imageStore(imageStoreId))->hasModel()){
					modelnames(k)=imparsVec_p[k].imageName+".model";
					(itsMappers.imageStore(k))->model()->unlock();
					controlRecord.define("modelnames", modelnames);
				}
			}
			
		}
		Vector<String> weightnames(Int(imparsVec_p.nelements()),"");
                Vector<String> pbnames(Int(imparsVec_p.nelements()),"");
		for(uInt k=0; k < imparsVec_p.nelements(); ++k){
			Int imageStoreId=k;
			if(k>0){
					if(gridparsVec_p[0].chanchunks >1 && uInt(itsMappers.nMappers()) >=    (gridparsVec_p[0].chanchunks + gridparsVec_p.nelements()))
						imageStoreId+=gridparsVec_p[0].chanchunks-1;
					if(gridparsVec_p[0].facets >1)
						imageStoreId+=gridparsVec_p[0].facets*gridparsVec_p[0].facets-1;
				}
			//cerr << "FTMachine name " << (itsMappers.getFTM2(imageStoreId))->name() << endl;
			if((itsMappers.getFTM2(imageStoreId))->useWeightImage()){
                          //cerr << "Mosaic weight image " << max(itsMappers.imageStore(imageStoreId)->weight(k)->get()) << endl;
				weightnames(k)=itsMappers.imageStore(imageStoreId)->weight(k)->name();
                                
                                
			}
			if(itsMakeVP){
			  pbnames(k)=itsMappers.imageStore(imageStoreId)->pb(k)->name();
                           (itsMappers.imageStore(imageStoreId)->pb(k))->unlock();
                        }
		}
		controlRecord.define("weightnames", weightnames);
                controlRecord.define("pbnames", pbnames);
        // Tell the child processes not to do the dividebyweight process as this is done
		// tell each child to do the normars stuff
		controlRecord.define("dividebyweight",  True);
                controlRecord.defineRecord("normpars", normpars_p); 
		///Let's see if no per chan weight density was chosen
		String weightdensityimage=getWeightDensity();
		if(weightdensityimage != "")
			controlRecord.define("weightdensity", weightdensityimage);
        
		Record vecSelParsRec, vecImParsRec, vecGridParsRec;
		Vector<Int>vecRange(2);
		for (uInt k = 0; k < dataSel_p.nelements(); ++k) {
			Record selparsRec = dataSel_p[k].toRecord();
			vecSelParsRec.defineRecord(String::toString(k), selparsRec);
		}
		for (uInt k=0; k < imparsVec_p.nelements(); ++k){
			Record imparsRec = imparsVec_p[k].toRecord();
			//need to send polrep
			imparsRec.define("polrep", Int((itsMappers.imageStore(k))->getDataPolFrame()));
                        //need to send movingSourceName if any
                        imparsRec.define("movingsource", movingSource_p);
			Record gridparsRec = gridparsVec_p[k].toRecord();
			/* Might need this to pass the state of the global ftmachines...test for parallel when needed
			String err;
			Record ftmrec, iftmrec;
			if(!( (itsMappers.getFTM2(k)->toRecord(err, iftmrec,False)) && (itsMappers.getFTM2(k, false)->toRecord(err, ftmrec,False))))
				throw(AipsError("FTMachines serialization failed"));
			gridparsRec.defineRecord("ftmachine", ftmrec);
			gridparsRec.defineRecord("iftmachine", iftmrec);
			*/
			vecImParsRec.defineRecord(String::toString(k), imparsRec);
			vecGridParsRec.defineRecord(String::toString(k), gridparsRec);
		}
		String workingdir="";
		//Int numchan=(dopsf) ? imstor->psf()->shape()[3] : imstor->residual()->shape() [3];
		//copy the imageinfo of main image here
		if(dopsf)
		  cubePsfImageInfo_p=(itsMappers.imageStore(0)->psf())->imageInfo();
        for(Int k=0; k < itsMappers.nMappers(); ++k){
	 
			if(dopsf){
			   
				for(uInt j =0; j <(itsMappers.imageStore(k)->getNTaylorTerms(true)); ++j){
                                  ///TESTOO
                                  //(itsMappers.imageStore(k))->psf(j)->set(0.0);
                                  /////////
					(itsMappers.imageStore(k))->psf(j)->unlock();
                                        (itsMappers.imageStore(k))->pb()->unlock();
				}
                                
			}
			else{
				for(uInt j =0; j <(itsMappers.imageStore(k)->getNTaylorTerms(false)); ++j){
                                  /////////TESTOO
                                  //(itsMappers.imageStore(k))->residual(j)->set(0.0);
                                    ///////
					(itsMappers.imageStore(k))->residual(j)->unlock();
				}
			}
			for(uInt j =0; j <(itsMappers.imageStore(k)->getNTaylorTerms(true)); ++j){
			//cerr << k << " type " << (itsMappers.imageStore(k))->sumwt(j)->imageType() << " name " << (itsMappers.imageStore(k))->sumwt(j)->name() << endl;
				
			       
				Path namewgt( (itsMappers.imageStore(k))->sumwt(j)->name());
				workingdir=namewgt.dirName();
                                ///TESTOO
                                //(itsMappers.imageStore(k))->sumwt(j)->set(0.0);
                                ////
				(itsMappers.imageStore(k))->sumwt(j)->unlock();
				//(itsMappers.imageStore(k))->releaseLocks();
			}
                        (itsMappers.imageStore(k))->releaseLocks();   

	}		
		//Send the working directory as the child and master may be at different places
		
		controlRecord.define("workingdirectory", workingdir);
		// For now this contains lastcycle if necessary in the future this
		// should come from the master control record
        
                //Int numprocs = applicator.numProcs();
        //cerr << "Number of procs: " << numprocs << endl;
        Int rank ( 0 );
        Bool assigned; //(casa::casa::applicator.nextAvailProcess(pwrite, rank));
        Bool allDone ( false );
        Vector<Bool> retvals(numchunks, False);
        Int indexofretval=0;
        for ( Int k=0; k < numchunks; ++k ) {
            assigned=casa::applicator.nextAvailProcess ( cmc, rank );
            //cerr << "assigned "<< assigned << endl;
            while ( !assigned ) {
                rank = casa::applicator.nextProcessDone ( cmc, allDone );
                //cerr << "while rank " << rank << endl;
                Bool status;
		Record returnRec;
		casa::applicator.get(returnRec);
                casa::applicator.get ( status );
                retvals(indexofretval)=status;
		if(dopsf)
		  updateImageBeamSet(returnRec);
		if(returnRec.isDefined("tempfilenames")){
		  std::vector<String> b=returnRec.asArrayString("tempfilenames").tovector();
		  tempFileNames_p.insert(std::end(tempFileNames_p), std::begin(b), std::end(b));
		}
		  
                ++indexofretval;
                if ( status )
                  //cerr << k << " rank " << rank << " successful " << endl;
                  cerr << "" ;
                else
                    logger << LogIO::SEVERE << k << " rank " << rank << " failed " << LogIO::POST;
                assigned = casa::applicator.nextAvailProcess ( cmc, rank );

            }
            std::this_thread::sleep_for(std::chrono::milliseconds(50));
            ///send process info
            // put data sel params #1
            applicator.put ( vecSelParsRec );
            // put image sel params #2
            applicator.put ( vecImParsRec );
            // put gridders params #3
            applicator.put ( vecGridParsRec );
            // put which channel to process #4
			vecRange(0)=startchan(k);
			vecRange(1)=endchan(k);
            applicator.put ( vecRange );
            // psf or residual CubeMajorCycleAlgorithm #5
            applicator.put ( dopsf );
            // store modelvis and other controls #6
            applicator.put ( controlRecord );
			// weighting scheme #7
			applicator.put( weightParams_p);
            /// Tell worker to process it
            applicator.apply ( cmc );

        }
        // Wait for all outstanding processes to return
        rank = casa::applicator.nextProcessDone ( cmc, allDone );
        while ( !allDone ) {
            Bool status;
	    Record returnRec;
	    casa::applicator.get(returnRec);
            casa::applicator.get ( status );
	    if(dopsf)
	      updateImageBeamSet(returnRec);
	    if(returnRec.isDefined("tempfilenames")){
	      std::vector<String> b=returnRec.asArrayString("tempfilenames").tovector();
	      tempFileNames_p.insert(std::end(tempFileNames_p), std::begin(b), std::end(b));
	    }
            retvals(indexofretval)=status;
            ++indexofretval;
            if ( status )
              //cerr << "remainder rank " << rank << " successful " << endl;
              cerr << "";
            else
                logger << LogIO::SEVERE << "remainder rank " << rank << " failed " << LogIO::POST;

            rank = casa::applicator.nextProcessDone ( cmc, allDone );
			if(casa::applicator.isSerial())
				allDone=true;
        }
        if(anyEQ(retvals, False)){
          //cerr << retvals << endl;
          ostringstream oss;
          oss << "One or more  of the cube section failed in de/gridding. Return values for "
              "the sections: " << retvals;
          throw(AipsError(oss));
        }
        if(!dopsf && normpars_p.isDefined("pblimit") && (normpars_p.asFloat("pblimit") > 0.0) ){
          try{
	    SIImageStore::copyMask(itsMappers.imageStore(0)->pb(), itsMappers.imageStore(0)->residual());
            (itsMappers.imageStore(0))->residual()->unlock();
            //(itsMappers.imageStore(0)->pb())->pixelMask().unlock();
            (itsMappers.imageStore(0))->pb()->unlock();
          }
          catch(AipsError &x) {
            if(!String(x.getMesg()).contains("T/F"))
              throw(AipsError(x.getMesg()));
	    else{
	      logger << LogIO::WARN << "Error : " << x.getMesg() << LogIO::POST;
	      //cout << "x.getMesg() " << endl;
	    }
            ///ignore copy mask error and proceed as this happens with interactive
          }
        }
	else{
	  LatticeLocker lock1 (*(itsMappers.imageStore(0)->psf()), FileLocker::Write);
	  itsMappers.imageStore(0)->psf()->setImageInfo(cubePsfImageInfo_p);
	  itsMappers.imageStore(0)->psf()->unlock();
          (itsMappers.imageStore(0))->pb()->unlock();
	}

        }  
	  
  
	return true;
  
  }
  /////////////////////////
  void SynthesisImagerVi2::predictModel(){
    LogIO os( LogOrigin("SynthesisImagerVi2","predictModel ",WHERE) );

    os << "---------------------------------------------------- Predict Model ---------------------------------------------" << LogIO::POST;
    
    Bool savemodelcolumn = !readOnly_p && useScratch_p;
    Bool savevirtualmodel = !readOnly_p && !useScratch_p;
    //os<<"PREDICTMODEL: readOnly_p=="<<readOnly_p<<" useScratch_p=="<<useScratch_p<<LogIO::POST;
    if( savemodelcolumn ) os << "Saving model column" << LogIO::POST;
    if( savevirtualmodel ) os << "Saving virtual model" << LogIO::POST;

    itsMappers.checkOverlappingModels("blank");


    {
      vi::VisBuffer2* vb = vi_p->getVisBuffer();;
      vi_p->originChunks();
      vi_p->origin();
      Double numberCoh=0;
      for (uInt k=0; k< mss_p.nelements(); ++k)
	numberCoh+=Double(mss_p[k]->nrow());

      ProgressMeter pm(1.0, numberCoh, "Predict Model", "","","",true);
      rownr_t cohDone=0;

      itsMappers.initializeDegrid(*vb);
      for (vi_p->originChunks(); vi_p->moreChunks();vi_p->nextChunk())
	{
	  
	  for (vi_p->origin(); vi_p->more(); vi_p->next())
	    {
	      //if (SynthesisUtilMethods::validate(*vb)==SynthesisUtilMethods::NOVALIDROWS) break; //No valid rows in this MS
	      //if !usescratch ...just save
	      vb->setVisCubeModel(Complex(0.0, 0.0));
	      itsMappers.degrid(*vb, savevirtualmodel);

	      if(savemodelcolumn && writeAccess_p ){

		const_cast<MeasurementSet& >((vi_p->ms())).lock(true);

		vi_p->writeVisModel(vb->visCubeModel());

	      //cerr << "nRows "<< vb->nRows() << "   " << max(vb->visCubeModel()) <<  endl;
		const_cast<MeasurementSet& >((vi_p->ms())).unlock();
	      }

	      cohDone += vb->nRows();
	      pm.update(Double(cohDone));

	    }
	}
      itsMappers.finalizeDegrid(*vb);
    }

    itsMappers.checkOverlappingModels("restore");
    itsMappers.releaseImageLocks();
    unlockMSs();
   
  }// end of predictModel

  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  void SynthesisImagerVi2::makeSdImage(Bool dopsf)
  {
    LogIO os( LogOrigin("SynthesisImagerVi2","makeSdImage",WHERE) );

    // Bool dopsf=false;
    if(datacol_p==FTMachine::PSF) dopsf=true;

    {
      vi::VisBuffer2* vb = vi_p->getVisBuffer();;
      vi_p->originChunks();
      vi_p->origin();

      Double numberCoh=0;
      for (uInt k=0; k< mss_p.nelements(); ++k)
        numberCoh+=Double(mss_p[k]->nrow());

      ProgressMeter pm(1.0, numberCoh, "Predict Model", "","","",true);
      rownr_t cohDone=0;

      itsMappers.initializeGrid(*vi_p,dopsf);
      for (vi_p->originChunks(); vi_p->moreChunks(); vi_p->nextChunk())
      {
        if (vi_p->getImpl()->isNewMs()) {
          itsMappers.handleNewMs(vi_p->ms());
        }
        for (vi_p->origin(); vi_p->more(); vi_p->next())
        {
          itsMappers.grid(*vb, dopsf, (refim::FTMachine::Type)datacol_p);
          cohDone += vb->nRows();
          pm.update(Double(cohDone));
        }
      }
      itsMappers.finalizeGrid(*vb, dopsf);

    }

    unlockMSs();

  }// end makeSDImage
  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  void SynthesisImagerVi2::makeImage(String type, const String& imagename, const String& complexImage, const Int whichModel)
  {
    LogIO os( LogOrigin("SynthesisImager","makeImage",WHERE) );

    
    refim::FTMachine::Type seType(refim::FTMachine::OBSERVED);
    if(type=="observed") {
      seType=refim::FTMachine::OBSERVED;
      os << LogIO::NORMAL // Loglevel INFO
         << "Making dirty image from " << type << " data "
	 << LogIO::POST;
    }
    else if (type=="model") {
      seType=refim::FTMachine::MODEL;
      os << LogIO::NORMAL // Loglevel INFO
         << "Making dirty image from " << type << " data "
	 << LogIO::POST;
    }
    else if (type=="corrected") {
      seType=refim::FTMachine::CORRECTED;
      os << LogIO::NORMAL // Loglevel INFO
         << "Making dirty image from " << type << " data "
	 << LogIO::POST;
    }
    else if (type=="psf") {
      seType=refim::FTMachine::PSF;
      os << "Making point spread function "
	 << LogIO::POST;
    }
    else if (type=="residual") {
      seType=refim::FTMachine::RESIDUAL;
      os << LogIO::NORMAL // Loglevel INFO
         << "Making dirty image from " << type << " data "
	 << LogIO::POST;
    }
    else if (type=="singledish-observed") {
      seType=refim::FTMachine::OBSERVED;
      os << LogIO::NORMAL 
         << "Making single dish image from observed data" << LogIO::POST;
    }
    else if (type=="singledish") {
      seType=refim::FTMachine::CORRECTED;
      os << LogIO::NORMAL 
         << "Making single dish image from corrected data" << LogIO::POST;
    }
    else if (type=="coverage") {
      seType=refim::FTMachine::COVERAGE;
      os << LogIO::NORMAL 
         << "Making single dish coverage function "
	 << LogIO::POST;
    }
    else if (type=="holography") {
      seType=refim::FTMachine::CORRECTED;
      os << LogIO::NORMAL
         << "Making complex holographic image from corrected data "
	 << LogIO::POST;
    }
    else if (type=="holography-observed") {
      seType=refim::FTMachine::OBSERVED;
      os << LogIO::NORMAL 
         << "Making complex holographic image from observed data "
	 << LogIO::POST;
    }

    String imageName=(itsMappers.imageStore(whichModel))->getName();
    String cImageName(complexImage);
    if(complexImage=="") {
      cImageName=imageName+".compleximage";
    }
    Bool keepComplexImage=(complexImage!="")||(type.contains("holography"));
    //cerr << "name " << (itsMappers.getFTM2(whichModel))->name() << endl;
 if(((itsMappers.getFTM2(whichModel))->name())!="MultiTermFTNew"){
   ////Non multiterm case    
	//cerr << "whichModel " << whichModel << " itsMappers num " << itsMappers.nMappers() << " shape " << (itsMappers.imageStore(whichModel))->getShape() << endl;
	///the SIImageStore sometimes return 0 shape...
	CoordinateSystem csys=itsMaxCoordSys;
	IPosition shp=itsMaxShape;
	if((itsMappers.imageStore(whichModel))->getShape().product()!=0 ){
		csys=	(itsMappers.imageStore(whichModel))->getCSys();
		shp=(itsMappers.imageStore(whichModel))->getShape();
	}
	itsMappers.releaseImageLocks();
    PagedImage<Float> theImage( shp, csys, imagename);
    PagedImage<Complex> cImageImage(theImage.shape(),
				    theImage.coordinates(),
				    cImageName);
    if(!keepComplexImage)
      cImageImage.table().markForDelete();


    Matrix<Float> weight;
    if(cImageImage.shape()[3] > 1){
		cImageImage.set(Complex(0.0));
		cImageImage.tempClose();
		makeComplexCubeImage(cImageName, seType, whichModel);
	}
	else{
    (itsMappers.getFTM2(whichModel))->makeImage(seType, *vi_p, cImageImage, weight);
	}
    if(seType==refim::FTMachine::PSF){
       StokesImageUtil::ToStokesPSF(theImage, cImageImage);
       StokesImageUtil::normalizePSF(theImage);
    }
    else{
      StokesImageUtil::To(theImage, cImageImage);
    }
 }
 else{
   ///Multiterm
   //refim::MultiTermFTNew *theft=static_cast<refim::MultiTermFTNew *>( (itsMappers.getFTM2(whichModel))->cloneFTM());
   refim::MultiTermFTNew *theft=static_cast<refim::MultiTermFTNew *>( (itsMappers.getFTM2(whichModel)).get());
   Int ntaylor= seType==refim::FTMachine::PSF ? theft->psfNTerms() : theft->nTerms();
   if(ntaylor<2)
     throw(AipsError("some issue with muti term setting "));
   Vector<CountedPtr<ImageInterface<Float> > >theImage(ntaylor);
   Vector<CountedPtr<ImageInterface<Complex> > >cImageImage(ntaylor);
   Vector<CountedPtr<Matrix<Float> > >weight(ntaylor);
   for (Int taylor=0; taylor < ntaylor; ++taylor){
     theImage[taylor]=new PagedImage<Float>((itsMappers.imageStore(whichModel))->getShape(), (itsMappers.imageStore(whichModel))->getCSys(), imagename+".tt"+String::toString(taylor));
     cImageImage[taylor]=new PagedImage<Complex> (theImage[taylor]->shape(),
						  theImage[taylor]->coordinates(),
						  cImageName+".tt"+String::toString(taylor));
      if(!keepComplexImage)
	static_cast<PagedImage<Complex> *> (cImageImage[taylor].get())->table().markForDelete();
      weight[taylor]=new Matrix<Float>();

   }
   theft->makeMTImages(seType, *vi_p, cImageImage, weight);
   Float maxpsf=1.0;
   for (Int taylor=0; taylor < ntaylor; ++taylor){
     if(seType==refim::FTMachine::PSF){
       StokesImageUtil::ToStokesPSF(*(theImage[taylor]), *(cImageImage[taylor]));
       if(taylor==0){
	 maxpsf=StokesImageUtil::normalizePSF(*theImage[taylor]);
	 //cerr << "maxpsf " << maxpsf << endl;
       }
       else{
	 ///divide by max;
	 (*theImage[taylor]).copyData((LatticeExpr<Float>)((*theImage[taylor])/maxpsf));
       }
    }
    else{
      StokesImageUtil::To(*(theImage[taylor]), *(cImageImage[taylor]));
    }
   }
   //delete theft;
     
 }
    unlockMSs();

  }// end makeImage
  /////////////////////////////////////////////////////

void SynthesisImagerVi2::makeComplexCubeImage(const String& cimage, const refim::FTMachine::Type imtype, const Int whichmodel){
	CubeMakeImageAlgorithm *cmi=new CubeMakeImageAlgorithm();
	// Dummies for using the right signature for init
	Path cimpath(cimage);
	String cimname=cimpath.absoluteName();
	//cerr << "ABSOLUTE path " << cimname << endl;
	Int argc = 1;
	char **argv = nullptr;
	casa::applicator.init(argc, argv);
	if(applicator.isController())
    {
		Record vecSelParsRec;
		for (uInt k = 0; k < dataSel_p.nelements(); ++k) {
			Record selparsRec = dataSel_p[k].toRecord();
			vecSelParsRec.defineRecord(String::toString(k), selparsRec);
		}
		Record imparsRec = imparsVec_p[whichmodel].toRecord();
		//cerr << "which model " << whichmodel << " record " << imparsRec << endl;
		Record gridparsRec = gridparsVec_p[whichmodel].toRecord();
        ///Break things into chunks for parallelization or memory availabbility
        Vector<Int> startchan;
        Vector<Int> endchan;
        Int numchunks;
        Int fudge_factor = 15;
        if ((itsMappers.getFTM2(0))->name()=="MosaicFTNew") {
            fudge_factor = 20;
        }
        else if ((itsMappers.getFTM2(0))->name()=="GridFT") {
            fudge_factor = 9;
        }
        TcleanProcessingInfo  procInfo;
        std::tie(procInfo, startchan, endchan)=nSubCubeFitInMemory(fudge_factor, itsMaxShape, gridpars_p.padding);
        numchunks = procInfo.chnchnks;
        
		Int imageType=Int(imtype);
		Int rank(0);
		Bool assigned; 
		Bool allDone(false);
		Vector<Int> chanRange(2);
		for (Int k=0; k < numchunks; ++k) {
			assigned=casa::applicator.nextAvailProcess ( *cmi, rank );
            //cerr << "assigned "<< assigned << endl;
            while ( !assigned ) {
                rank = casa::applicator.nextProcessDone ( *cmi, allDone );
                //cerr << "while rank " << rank << endl;
                Bool status;
                casa::applicator.get(status);
                /*if ( status )
                  cerr << k << " rank " << rank << " successful " << endl;
                else
                    cerr << k << " rank " << rank << " failed " << endl;
                */
                assigned = casa::applicator.nextAvailProcess ( *cmi, rank );

            }
            applicator.put(vecSelParsRec);
			// put image sel params #2
			applicator.put(imparsRec);
			// put gridder params #3
			applicator.put(gridparsRec);
			// put which channel to process #4
			chanRange(0)=startchan(k);
			chanRange(1)=endchan(k);
			applicator.put(chanRange);
			//Type of image #5
			applicator.put(imageType);
			// weighting parameters #6
			applicator.put( weightParams_p);
			// complex imagename #7
			applicator.put(cimname);
			applicator.apply(*cmi);
		}
		// Wait for all outstanding processes to return
        rank = casa::applicator.nextProcessDone(*cmi, allDone);
        while (!allDone) {
            Bool status;
            casa::applicator.get(status);
            /*
            if(status)
                cerr << "remainder rank " << rank << " successful " << endl;
            else
                cerr << "remainder rank " << rank << " failed " << endl;
            */
            rank = casa::applicator.nextProcessDone(*cmi, allDone);
			if(casa::applicator.isSerial())
				allDone=true;
        }
    }//applicator controller section
}



CountedPtr<SIMapper> SynthesisImagerVi2::createSIMapper(String mappertype,  
							   CountedPtr<SIImageStore> imagestore,
							CountedPtr<refim::FTMachine> ftmachine,
							CountedPtr<refim::FTMachine> iftmachine,
						       uInt /*ntaylorterms*/)
  {
    LogIO os( LogOrigin("SynthesisImagerVi2","createSIMapper",WHERE) );
    
    CountedPtr<SIMapper> localMapper;

    try
      {
	
	if( mappertype == "default" || mappertype == "multiterm" )
	  {
	    localMapper = new SIMapper( imagestore, ftmachine, iftmachine );
	  }
	else if( mappertype == "imagemosaic") // || mappertype == "mtimagemosaic" )
	  {
	    localMapper = new SIMapperImageMosaic( imagestore, ftmachine, iftmachine );
	  }
	else
	  {
	    throw(AipsError("Unknown mapper type : " + mappertype));
	  }

      }
    catch(AipsError &x) {
	throw(AipsError("Error in createSIMapper : " + x.getMesg() ) );
      }
    return localMapper;
  }
  
void SynthesisImagerVi2::lockMS(MeasurementSet& thisms){
  thisms.lock(!readOnly_p);
    thisms.antenna().lock(false);
	thisms.dataDescription().lock(false);
    thisms.feed().lock(false);
    thisms.field().lock(false);
    thisms.observation().lock(false);
    thisms.polarization().lock(false);
    thisms.processor().lock(false);
	thisms.spectralWindow().lock(false);
    thisms.state().lock(false);
    if(!thisms.doppler().isNull())
      thisms.doppler().lock(false);
    if(!thisms.flagCmd().isNull())
      thisms.flagCmd().lock(false);
    if(!thisms.freqOffset().isNull())
      thisms.freqOffset().lock(false);
	//True here as we can write in that
    if(!thisms.history().isNull())
    // thisms.history().lock(!readOnly_p);
      thisms.history().lock(false);
    if(!thisms.pointing().isNull())
      thisms.pointing().lock(false);
	//we write virtual model here
    if(!thisms.source().isNull())
      thisms.source().lock(!readOnly_p);
    if(!thisms.sysCal().isNull())
      thisms.sysCal().lock(false);
    if(!thisms.weather().isNull())
      thisms.weather().lock(false);
	
}
///////////////////////
  
  Vector<SynthesisParamsSelect> SynthesisImagerVi2::tuneSelectData(){
    vi_p->originChunks();
    vi_p->origin();
    vi::VisBuffer2* vb=vi_p->getVisBuffer();
    Int spwnow=vb->spectralWindows()[0];
    Int nchaninms=MSColumns(vb->ms()).spectralWindow().numChan()(spwnow);
    //For some small number of channels in the ms vi/vb2 retuning the selection
    //will sometimes return more channels than what is in the ms...this is a
    //kludge here to bypass it...mostly seen in test_tclean
    /// write to the test !!  till someboody fixes this is vi2 or wait for cngi
    //if savescratch column we have tune...otherwise some channel may be 0
    // when chunking or in parallel
    //cerr << "nchanims " << nchaninms << endl;
    if(nchaninms <30 && !(!readOnly_p && useScratch_p))
      return dataSel_p;
    
    return SynthesisImager::tuneSelectData();
  }
//////////////////////
  void SynthesisImagerVi2::lockMSs(){
    for(uInt i=0;i<mss_p.nelements();i++)
      { 
       
	MeasurementSet *ms_l = 	const_cast<MeasurementSet* >(mss_p[i]);
	lockMS(*ms_l);
      }

  }
void SynthesisImagerVi2::unlockMSs()
  {
    LogIO os( LogOrigin("SynthesisImagerVi2","unlockMSs",WHERE) );
    for(uInt i=0;i<mss_p.nelements();i++)
      { 
	os << LogIO::NORMAL2 << "Unlocking : " << (mss_p[i])->tableName() << LogIO::POST;
	MeasurementSet *ms_l = 	const_cast<MeasurementSet* >(mss_p[i]);
	ms_l->unlock(); 
	ms_l->antenna().unlock();
	ms_l->dataDescription().unlock();
	ms_l->feed().unlock();
	ms_l->field().unlock();
	ms_l->observation().unlock();
	ms_l->polarization().unlock();
	ms_l->processor().unlock();
	ms_l->spectralWindow().unlock();
	ms_l->state().unlock();
	//
	// Unlock the optional sub-tables as well, if they are present
	//
	if(!(ms_l->source().isNull()))     ms_l->source().unlock();
	if(!(ms_l->doppler().isNull()))    ms_l->doppler().unlock();
	if(!(ms_l->flagCmd().isNull()))    ms_l->flagCmd().unlock();
	if(!(ms_l->freqOffset().isNull())) ms_l->freqOffset().unlock();
	if(!(ms_l->history().isNull()))    ms_l->history().unlock();
	if(!(ms_l->pointing().isNull()))   ms_l->pointing().unlock();
	if(!(ms_l->sysCal().isNull()))     ms_l->sysCal().unlock();
	if(!(ms_l->weather().isNull()))    ms_l->weather().unlock();
      }
  }
  void SynthesisImagerVi2::createFTMachine(
            CountedPtr<refim::FTMachine>& theFT,
            CountedPtr<refim::FTMachine>& theIFT,
            const String& ftname,
            const uInt nTaylorTerms,
            const String mType,
            const Int facets,            //=1
            //------------------------------
            const Int wprojplane,        //=1,
            const Float padding,         //=1.0,
            const Bool useAutocorr,      //=false,
            const Bool useDoublePrec,    //=true,
            const String gridFunction,   //=String("SF"),
            //------------------------------
            const Bool aTermOn,          //= true,
            const Bool psTermOn,         //= true,
            const Bool mTermOn,          //= false,
            const Bool wbAWP,            //= true,
            const String cfCache,        //= "",
            const Bool usePointing,      //= false,
            // const Vector<Float> pointingOffsetSigDev, //= 10.0,
            const vector<float> pointingOffsetSigDev, // = {10,10}
            const Bool doPBCorr,         //= true,
            const Bool conjBeams,        //= true,
            const Float computePAStep,   //=360.0
            const Float rotatePAStep,    //=5.0
            const String interpolation,  //="linear"
            const Bool freqFrameValid,   //=true
            const Int cache,             //=1000000000,
            const Int tile,              //=16
            const String stokes,         //=I
            const String imageNamePrefix,
            //---------------------------
            const String &pointingDirCol,
            const String &convertFirst,
            const Float skyPosThreshold,
            const Int convSupport,
            const Quantity &truncateSize,
            const Quantity &gwidth,
            const Quantity &jwidth,
            const Float minWeight,
            const Bool clipMinMax,
            const Bool pseudoI
           )

  {
    LogIO os( LogOrigin("SynthesisImagerVi2","createFTMachine",WHERE));


    if (ftname == "gridft") {
      if (facets >1) {
        theFT = new refim::GridFT(
            cache, tile,
            gridFunction, mLocation_p, phaseCenter_p, padding,
            useAutocorr, useDoublePrec
        );
        theIFT = new refim::GridFT(
            cache, tile,
            gridFunction, mLocation_p, phaseCenter_p, padding,
            useAutocorr, useDoublePrec
        );
      }
      else {
        theFT = new refim::GridFT(
            cache, tile,
            gridFunction, mLocation_p, padding,
            useAutocorr, useDoublePrec
        );
        theIFT = new refim::GridFT(
            cache, tile,
            gridFunction, mLocation_p, padding,
            useAutocorr, useDoublePrec
        );
      }
    }
    else if (ftname == "wprojectft") {
      Double maxW = -1.0;
      Double minW = -1.0;
      Double rmsW = -1.0;
      if (wprojplane < 1)
        casa::refim::WProjectFT::wStat(*vi_p, minW, maxW, rmsW);
      if (facets > 1) {
        theFT = new refim::WProjectFT(wprojplane, phaseCenter_p, mLocation_p,
         cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
        theIFT = new refim::WProjectFT(wprojplane, phaseCenter_p, mLocation_p,
          cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
      }
      else {
        theFT = new refim::WProjectFT(wprojplane, mLocation_p,
          cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
        theIFT = new refim::WProjectFT(wprojplane, mLocation_p,
          cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
      }
      CountedPtr<refim::WPConvFunc> sharedconvFunc =
        static_cast<refim::WProjectFT &>(*theFT).getConvFunc();
      //static_cast<WProjectFT &>(*theFT).setConvFunc(sharedconvFunc);
      static_cast<refim::WProjectFT &>(*theIFT).setConvFunc(sharedconvFunc);
    }

    else if ( ftname == "mosaic" || ftname== "mosft" || ftname == "mosaicft" || ftname== "MosaicFT" || ftname == "awp2"){

      createMosFTMachine(theFT, theIFT, padding, useAutocorr, useDoublePrec, rotatePAStep, stokes, conjBeams);
    } 
    else if ((ftname.at(0,3)=="awp") || (ftname== "mawprojectft") || (ftname == "protoft")) {
      createAWPFTMachine(theFT, theIFT, ftname, facets, wprojplane, 
			 padding, useAutocorr, useDoublePrec, gridFunction,
			 aTermOn, psTermOn, mTermOn, wbAWP, cfCache, 
			 usePointing, pointingOffsetSigDev, doPBCorr, conjBeams, computePAStep,
			 rotatePAStep, cache,tile,imageNamePrefix);
    }

    else if ( ftname == "sd" ) {
      createSDFTMachine(
        theFT, theIFT,
        pointingDirCol, convertFirst,
        skyPosThreshold, doPBCorr, rotatePAStep,
        gridFunction, convSupport, truncateSize, gwidth, jwidth,
        minWeight, clipMinMax, cache, tile, stokes
      );
    }
    else {
      throw( AipsError( "Invalid FTMachine name : " + ftname ) );
    }

    // Now, clone and pack the chosen FT into a MultiTermFT if needed.
    if( mType == "multiterm" ) {
      AlwaysAssert( nTaylorTerms>=1 , AipsError );

      CountedPtr<refim::FTMachine> theMTFT =
        new refim::MultiTermFTNew( theFT, nTaylorTerms, true /*forward*/ );
      CountedPtr<refim::FTMachine> theMTIFT =
        new refim::MultiTermFTNew( theIFT, nTaylorTerms, false /*forward*/ );

      theFT = theMTFT;
      theIFT = theMTIFT;
    }

    // Now, set the SkyJones if needed, and if not internally generated.
    if ( mType == "imagemosaic" and
        ( ftname != "awprojectft" and
          ftname != "mawprojectft" and
          ftname != "proroft" ) ) {
      CountedPtr<refim::SkyJones> vp;
      MSColumns msc(*(mss_p[0]));
      Quantity parang(0.0,"deg");
      Quantity skyposthreshold(0.0,"deg");
      vp = new refim::VPSkyJones(
        msc, true,  parang, BeamSquint::NONE, skyposthreshold
      );

      Vector<CountedPtr<refim::SkyJones> > skyJonesList(1);
      skyJonesList(0) = vp;
      theFT->setSkyJones( skyJonesList );
      theIFT->setSkyJones( skyJonesList );
    }

    // For mode=cubedata, set the freq frame to invalid.
    // get this info from buildCoordSystem
    // theFT->setSpw( tspws, false );
    // theIFT->setSpw( tspws, false );
    theFT->setFrameValidity( freqFrameValid );
    theIFT->setFrameValidity( freqFrameValid );

    // Set interpolation mode
    theFT->setFreqInterpolation( interpolation );
    theIFT->setFreqInterpolation( interpolation );

    // Set tracking of moving source if any
    if (movingSource_p != "") {
      theFT->setMovingSource(movingSource_p);
      theIFT->setMovingSource(movingSource_p);
    }
    /* vi_p has chanselection now
    //channel selections from spw param
    theFT->setSpwChanSelection(chanSel_p);
    theIFT->setSpwChanSelection(chanSel_p);
    */

    // Set pseudo-I if requested.
    if (pseudoI == true) {
      os << "Turning on Pseudo-I gridding" << LogIO::POST;
      theFT->setPseudoIStokes(true);
      theIFT->setPseudoIStokes(true);
    }

  }

  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  void SynthesisImagerVi2::createAWPFTMachine(CountedPtr<refim::FTMachine>& theFT, CountedPtr<refim::FTMachine>& theIFT, 
					      const String& ftmName,
					      const Int,// facets,            //=1
					      //------------------------------
					      const Int wprojPlane,        //=1,
					      const Float,// padding,         //=1.0,
					      const Bool,// useAutocorr,      //=false,
					      const Bool useDoublePrec,    //=true,
					      const String,// gridFunction,   //=String("SF"),
					      //------------------------------
					      const Bool aTermOn,          //= true,
					      const Bool psTermOn,         //= true,
					      const Bool mTermOn,          //= false,
					      const Bool wbAWP,            //= true,
					      const String cfCache,        //= "",
					      const Bool usePointing,       //= false,
					      const vector<Float> pointingOffsetSigDev,//={10,10},
					      const Bool doPBCorr,         //= true,
					      const Bool conjBeams,        //= true,
					      const Float computePAStep,   //=360.0
					      const Float rotatePAStep,    //=5.0
					      const Int cache,             //=1000000000,
					      const Int tile,               //=16
					      const String imageNamePrefix
					      )

  {
    LogIO os( LogOrigin("SynthesisImagerVi2","createAWPFTMachine",WHERE));

    if (wprojPlane<=1)
      {
	os << LogIO::NORMAL
	   << "You are using wprojplanes=1. Doing co-planar imaging (no w-projection needed)" 
	   << LogIO::POST;
	os << LogIO::NORMAL << "Performing WBA-Projection" << LogIO::POST; // Loglevel PROGRESS
      }
    // if((wprojPlane>1)&&(wprojPlane<64)) 
    //   {
    // 	os << LogIO::WARN
    // 	   << "No. of w-planes set too low for W projection - recommend at least 128"
    // 	   << LogIO::POST;
    // 	os << LogIO::NORMAL << "Performing WBAW-Projection" << LogIO::POST; // Loglevel PROGRESS
    //   }

    // CountedPtr<ATerm> apertureFunction = createTelescopeATerm(mss4vi_p[0], aTermOn);
    // CountedPtr<PSTerm> psTerm = new PSTerm();
    // CountedPtr<WTerm> wTerm = new WTerm();
    
    // //
    // // Selectively switch off CFTerms.
    // //
    // if (aTermOn == false) {apertureFunction->setOpCode(CFTerms::NOOP);}
    // if (psTermOn == false) psTerm->setOpCode(CFTerms::NOOP);

    // //
    // // Construct the CF object with appropriate CFTerms.
    // //
    // CountedPtr<ConvolutionFunction> tt;
    // tt = AWProjectFT::makeCFObject(aTermOn, psTermOn, true, mTermOn, wbAWP);
    // CountedPtr<ConvolutionFunction> awConvFunc;
    // //    awConvFunc = new AWConvFunc(apertureFunction,psTerm,wTerm, !wbAWP);
    // if ((ftmName=="mawprojectft") || (mTermOn))
    //   awConvFunc = new AWConvFuncEPJones(apertureFunction,psTerm,wTerm,wbAWP);
    // else
    //   awConvFunc = new AWConvFunc(apertureFunction,psTerm,wTerm,wbAWP);

    MSObservationColumns msoc((mss_p[0])->observation());
    String telescopeName=msoc.telescopeName()(0);
    CountedPtr<refim::ConvolutionFunction> awConvFunc = refim::AWProjectFT::makeCFObject(telescopeName, 
									   aTermOn,
									   psTermOn, (wprojPlane > 1),
									   mTermOn, wbAWP, conjBeams);

    //
    // Construct the appropriate re-sampler.
    //
    CountedPtr<refim::VisibilityResamplerBase> visResampler;
    //    if (ftmName=="protoft") visResampler = new ProtoVR();
    //elsef
    //cerr << "############ Starting " << ftmName << endl;
    if (ftmName=="awphpg")
      {  //hpg::initialize();
#ifndef USE_HPG
        throw(AipsError("Code has not been built with gpu libraries"));
#else
      try {
        // test for cuda
        void *cudalib = dlopen("libcuda.so", RTLD_LAZY);
        if (!cudalib)
          throw(AipsError("Cannot run hpg on this machine"));
        else {
          dlclose(cudalib);
        }
        if (!hpg::is_initialized())
          hpg::initialize();
      } catch (...) {
        throw(AipsError("Trying to use GPU code with the wrong GPU or no GPU"));
      }
      cerr << "DOING HPG" << endl;
      visResampler = new refim::AWVisResamplerHPG(false);
      visResampler->setModelImage("");
#endif
      }
    else
      visResampler = new refim::AWVisResampler();
    //    CountedPtr<VisibilityResamplerBase> visResampler = new VisibilityResampler();

    //
    // Construct and initialize the CF cache object.
    //


    // CountedPtr<CFCache> cfCacheObj = new CFCache();
    // cfCacheObj->setCacheDir(cfCache.data());
    // //    cerr << "Setting wtImagePrefix to " << imageNamePrefix.c_str() << endl;
    // cfCacheObj->setWtImagePrefix(imageNamePrefix.c_str());
    // cfCacheObj->initCache2();

    CountedPtr<refim::CFCache> cfCacheObj;
      

    //
    // Finally construct the FTMachine with the CFCache, ConvFunc and
    // Re-sampler objects.  
    //
    Float pbLimit_l=1e-3;
    if(ftmName=="awphpg"){
      theFT=new refim::AWProjectWBFTHPG(wprojPlane, cache/2, 
					   cfCacheObj, awConvFunc,
					   visResampler,
					   /*true */usePointing, pointingOffsetSigDev, doPBCorr,
					   tile, computePAStep, pbLimit_l, true,conjBeams,
					   useDoublePrec);
      theFT->setPBReady(true);
    }
    else{
      theFT = new refim::AWProjectWBFT(wprojPlane, cache/2,
				       cfCacheObj, awConvFunc,
				       visResampler,
				       /*true */usePointing, pointingOffsetSigDev ,doPBCorr,
				       tile, computePAStep, pbLimit_l, true,conjBeams,
				       useDoublePrec);
    }
  if(ftmName != "awphpg"){
    cfCacheObj = new refim::CFCache();
    cfCacheObj->setCacheDir(cfCache.data());
    // Get the LAZYFILL setting from the user configuration.  If not
    // found, default to False.
    //
    // With lazy fill ON, CFCache loads the required CFs on-demand
    // from the disk.  And periodically triggers garbage collection to
    // release CFs that aren't required immediately.
    //cfCacheObj->setLazyFill(refim::SynthesisUtils::getenv("CFCache.LAZYFILL",1)==1);
    cfCacheObj->setLazyFill(False);

    //    cerr << "Setting wtImagePrefix to " << imageNamePrefix.c_str() << endl;
    cfCacheObj->setWtImagePrefix(imageNamePrefix.c_str());
    cfCacheObj->initCache2(CFC_VERBOSE);

    theFT->setCFCache(cfCacheObj);
    Quantity rotateOTF(rotatePAStep,"deg");
    static_cast<refim::AWProjectWBFT &>(*theFT).setPAIncrement(Quantity(computePAStep,"deg"),rotateOTF);
    static_cast<refim::AWProjectWBFT &>(*theFT).setObservatoryLocation(mLocation_p);
  }

   
   

    // theIFT = new AWProjectWBFT(wprojPlane, cache/2, 
    // 			       cfCacheObj, awConvFunc, 
    // 			       visResampler,
    // 			       /*true */usePointing, doPBCorr, 
    // 			       tile, computePAStep, pbLimit_l, true,conjBeams,
    // 			       useDoublePrec);

    // static_cast<AWProjectWBFT &>(*theIFT).setObservatoryLocation(mLocation_p);
    // static_cast<AWProjectWBFT &>(*theIFT).setPAIncrement(Quantity(computePAStep,"deg"),rotateOTF);
    if(ftmName=="awphpg"){
      /// the gridder and degridder (for hpg) are the same except it needs to load the modelimage
      static_cast<refim::AWProjectWBFTHPG &>(*theFT).setObservatoryLocation(mLocation_p);
      //static_cast<refim::AWProjectWBFTHPG &>(*theFT).setPAIncrement(Quantity(computePAStep,"deg"),rotateOTF);
      theIFT = theFT;
    }
    else
      theIFT = new refim::AWProjectWBFT(static_cast<refim::AWProjectWBFT &>(*theFT));

    os << "Sending frequency selection information " <<  mssFreqSel_p  <<  " to AWP FTM." << LogIO::POST;
    theFT->setSpwFreqSelection( mssFreqSel_p );
    theIFT->setSpwFreqSelection( mssFreqSel_p );

  }

  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

  void SynthesisImagerVi2:: createMosFTMachine(CountedPtr<refim::FTMachine>& theFT,CountedPtr<refim::FTMachine>&  theIFT, const Float /*padding*/, const Bool useAutoCorr, const Bool useDoublePrec, const Float rotatePAStep, const String stokes, const Bool doConjBeams){
    
    LogIO os(LogOrigin("SynthesisImagerVi2", "createMosFTMachine",WHERE));
   
    MSColumns msc(vi_p->ms());
    String telescop=msc.observation().telescopeName()(0);
    Bool multiTel=False;
    Int msid=0;
     for(vi_p->originChunks(); vi_p->moreChunks(); vi_p->nextChunk()){
       if(((vi_p->getVisBuffer())->msId() != msid) && telescop !=  MSColumns(vi_p->ms()).observation().telescopeName()(0)){
	 msid=(vi_p->getVisBuffer())->msId();
	 multiTel=True;
       }
     }
    vi_p->originChunks();
  
  

    PBMath::CommonPB kpb;
    Record rec;
    getVPRecord( rec, kpb, telescop );
   

    if(rec.empty()){os << LogIO::SEVERE << "Cannot proceed with mosaicft gridder without a valid PB model" << LogIO::POST; }
    
    /*
   VPManager *vpman=VPManager::Instance();
    PBMath::CommonPB kpb;
    PBMath::enumerateCommonPB(telescop, kpb);
    Record rec;
    vpman->getvp(rec, telescop);
    */

   refim::VPSkyJones* vps= nullptr;
   //cerr << "rec " << rec << " kpb " << kpb << endl;
   //cerr <<  "createMOs ftname " <<  gridpars_p.ftmachine <<  endl;
   if (!gridpars_p.ftmachine.contains("mos")) {
     cerr <<  "PASTERP " <<  rotatePAStep <<  "   " <<  gridpars_p.computePAStep <<  endl;
     bool dosquint = (gridpars_p.computePAStep < 180);       //anything beneath 180 deg ...you are not serious about squint correction  
    //  TESTOO
    //dosquint = False;
    ///////
    
    cerr <<  "Doing AWPLPG" <<   " wprojplanes " << gridpars_p.wprojplanes << endl;
     theFT = new refim::AWPLPG(vps , gridpars_p.wprojplanes, dosquint, gridpars_p.computePAStep*(C::pi)/180.0, mLocation_p, stokes, useAutoCorr, useDoublePrec, gridpars_p.usePointing);
     theIFT = new refim::AWPLPG(vps , gridpars_p.wprojplanes, dosquint, gridpars_p.computePAStep*(C::pi)/180.0, mLocation_p, stokes, useAutoCorr, useDoublePrec, gridpars_p.usePointing);
     CountedPtr<refim::SimplePBConvFunc> mospb=new refim::HetArrayConvFunc();
      static_cast<refim::AWPLPG &>(*theFT).setConvFunc(mospb);
      static_cast<refim::AWPLPG &>(*theIFT).setConvFunc(mospb);
      
   }
   else{
    if(rec.asString("name")=="COMMONPB" && kpb !=PBMath::UNKNOWN ){
      vps= new refim::VPSkyJones(msc, true, Quantity(rotatePAStep, "deg"), BeamSquint::GOFIGURE, Quantity(360.0, "deg"));
      /////Don't know which parameter has pb threshold cutoff that the user want 
      ////leaving at default
      ////vps.setThreshold(minPB);
      
    }
    else{
      PBMath myPB(rec);
      String whichPBMath;
      PBMathInterface::namePBClass(myPB.whichPBClass(), whichPBMath);
      os  << "Using the PB defined by " << whichPBMath << " for beam calculation for telescope " << telescop << LogIO::POST;
      vps= new refim::VPSkyJones(telescop, myPB, Quantity(rotatePAStep, "deg"), BeamSquint::GOFIGURE, Quantity(360.0, "deg"));
      //kpb=PBMath::DEFAULT;
    }
   
    theFT = new refim::MosaicFTNew(vps, mLocation_p, stokes, 1000000000, 16, useAutoCorr, 
				   useDoublePrec, doConjBeams, gridpars_p.usePointing);
    PBMathInterface::PBClass pbtype=((kpb==PBMath::EVLA) || multiTel)? PBMathInterface::COMMONPB: PBMathInterface::AIRY;
    if(rec.asString("name")=="IMAGE")
       pbtype=PBMathInterface::IMAGE;
    ///Use Heterogenous array mode for the following
    if((kpb == PBMath::UNKNOWN) || (kpb==PBMath::OVRO) || (kpb==PBMath::ACA)
       || (kpb==PBMath::ALMA) || (kpb==PBMath::EVLA) || multiTel){
      CountedPtr<refim::SimplePBConvFunc> mospb=new refim::HetArrayConvFunc(pbtype, itsVpTable);
      static_cast<refim::MosaicFTNew &>(*theFT).setConvFunc(mospb);
    }
    ///////////////////make sure both FTMachine share the same conv functions.
    theIFT= new refim::MosaicFTNew(static_cast<refim::MosaicFTNew &>(*theFT));
   }
  }
  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

  void SynthesisImagerVi2::createSDFTMachine(CountedPtr<refim::FTMachine>& theFT,
      CountedPtr<refim::FTMachine>& theIFT,
      const String &pointingDirCol,
      const String &convertFirst,
      const Float skyPosThreshold,
      const Bool /*doPBCorr*/,
      const Float rotatePAStep,
      const String& gridFunction,
      const Int convSupport,
      const Quantity& truncateSize,
      const Quantity& gwidth,
      const Quantity& jwidth,
      const Float minWeight,
      const Bool clipMinMax,
      const Int cache,
      const Int tile,
      const String &stokes,
      const Bool pseudoI) {
//    // member variable itsVPTable is VP table name
    LogIO os(LogOrigin("SynthesisImagerVi2", "createSDFTMachine", WHERE));
    os << LogIO::NORMAL // Loglevel INFO
       << "Performing single dish gridding..." << LogIO::POST;
    os << LogIO::NORMAL1 // gridFunction is too cryptic for most users.
       << "with convolution function " << gridFunction << LogIO::POST;

    // Now make the Single Dish Gridding
    os << LogIO::NORMAL // Loglevel INFO
       << "Gridding will use specified common tangent point:" << LogIO::POST;
    os << LogIO::NORMAL << SynthesisUtilMethods::asComprehensibleDirectionString(phaseCenter_p)
        << LogIO::POST; // Loglevel INFO
    String const gridfunclower = downcase(gridFunction);
    if(gridfunclower=="pb") {
      refim::SkyJones *vp = nullptr;
      MSColumns msc(*mss_p[0]);
      // todo: NONE is OK?
      BeamSquint::SquintType squintType = BeamSquint::NONE;
      Quantity skyPosThresholdQuant((Double)skyPosThreshold+180.0, "deg");
      Quantity parAngleStepQuant((Double)rotatePAStep, "deg");
      if (itsVpTable.empty()) {
        os << LogIO::NORMAL // Loglevel INFO
            << "Using defaults for primary beams used in gridding" << LogIO::POST;
        vp=new refim::VPSkyJones(msc, true, parAngleStepQuant, squintType,
            skyPosThresholdQuant);
      } else {
        os << LogIO::NORMAL // Loglevel INFO
            << "Using VP as defined in " << itsVpTable <<  LogIO::POST;
        Table vpTable( itsVpTable );
        vp=new refim::VPSkyJones(msc, vpTable, parAngleStepQuant, squintType,
            skyPosThresholdQuant);
      }
      theFT = new refim::SDGrid(mLocation_p, *vp, cache/2, tile, gridfunclower,
          convSupport, minWeight, clipMinMax);
      theIFT = new refim::SDGrid(mLocation_p, *vp, cache/2, tile, gridfunclower,
          convSupport, minWeight, clipMinMax);
    } else if (gridfunclower=="gauss" || gridfunclower=="gjinc") {
      DirectionCoordinate dirCoord = itsMaxCoordSys.directionCoordinate();
      Vector<String> units = dirCoord.worldAxisUnits();
      Vector<Double> increments = dirCoord.increment();
      Quantity cellx(increments[0], units[0]);
      Quantity celly(increments[1], units[1]);
      if (cellx != celly &&
          ((!truncateSize.getUnit().empty()||truncateSize.getUnit()=="pixel")
              || (!gwidth.getUnit().empty()||gwidth.getUnit()=="pixel")
              || (!jwidth.getUnit().empty()||jwidth.getUnit()=="pixel"))) {
        os << LogIO::WARN
            << "The " << gridFunction << " gridding doesn't support non-square grid." << endl
            << "Result may be wrong." << LogIO::POST;
      }
      Float truncateValue, gwidthValue, jwidthValue;
      if (truncateSize.getUnit().empty() || truncateSize.getUnit()=="pixel")
        truncateValue = truncateSize.getValue();
      else
        truncateValue = truncateSize.getValue("rad")/celly.getValue("rad");
      if (gwidth.getUnit().empty() || gwidth.getUnit()=="pixel")
        gwidthValue = gwidth.getValue();
      else
        gwidthValue = gwidth.getValue("rad")/celly.getValue("rad");
      if (jwidth.getUnit().empty() || jwidth.getUnit()=="pixel")
        jwidthValue = jwidth.getValue();
      else
        jwidthValue = jwidth.getValue("rad")/celly.getValue("rad");
      theFT = new refim::SDGrid(mLocation_p, cache/2, tile, gridfunclower,
                        truncateValue, gwidthValue, jwidthValue, minWeight, clipMinMax);
      theIFT = new refim::SDGrid(mLocation_p, cache/2, tile, gridfunclower,
                        truncateValue, gwidthValue, jwidthValue, minWeight, clipMinMax);
    }
    else {
      theFT = new refim::SDGrid(mLocation_p, cache/2, tile, gridfunclower,
                        convSupport, minWeight, clipMinMax);
      theIFT = new refim::SDGrid(mLocation_p, cache/2, tile, gridfunclower,
                        convSupport, minWeight, clipMinMax);
    }
    theFT->setPointingDirColumn(pointingDirCol);
    static_cast<refim::SDGrid*>(theFT.get())->setConvertFirst(convertFirst);
    static_cast<refim::SDGrid*>(theIFT.get())->setConvertFirst(convertFirst);

    // turn on Pseudo Stokes mode if necessary
    if (pseudoI || stokes == "XX" || stokes == "YY" || stokes == "XXYY"
        || stokes == "RR" || stokes == "LL" || stokes == "RRLL") {
      theFT->setPseudoIStokes(True);
      theIFT->setPseudoIStokes(True);
    }
  }
  
  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  //// Get/Set Weight Grid.... write to disk and read

  /// todo : do for full mapper list, and taylor terms.
  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  Bool SynthesisImagerVi2::setWeightDensity( const String& weightimagename)
  {
    LogIO os(LogOrigin("SynthesisImagerVi2", "setWeightDensity()", WHERE));
    try
      {
	if(weightimagename.size() !=0){
	  Table::isReadable(weightimagename, True);
	  PagedImage<Float> im(weightimagename);
	  imwgt_p=VisImagingWeight(im);
          im.unlock();
	}
	else{
          ////In memory weight densities is being deprecated...we should get rid of this bit
	  Int ndensities=1;
	  if((itsMappers.imageStore(0)->gridwt()->nelements())==5)
	    ndensities=(itsMappers.imageStore(0)->gridwt())->shape()[4];
	  Int nx=(itsMappers.imageStore(0)->gridwt())->shape()[0];
	  Int ny=(itsMappers.imageStore(0)->gridwt())->shape()[1];
	  Block<Matrix<Float> > densitymatrices(ndensities);
	  if(((itsMappers.imageStore(0)->gridwt())->shape().nelements())==5){
	    IPosition blc(Vector<Int>(5,0));
	    for (uInt fid=0;fid<densitymatrices.nelements();fid++)
	      {
		densitymatrices[fid].resize();
		Array<Float> lala;
		blc[4]=fid;
		itsMappers.imageStore(0)->gridwt()->getSlice(lala, blc, IPosition(5, nx, ny,1,1,1), True);
		densitymatrices[fid].reference( lala.reform(IPosition(2, nx, ny)));
	      }
	  }
	  else{
	    Array<Float> lala;
	    itsMappers.imageStore(0)->gridwt()->get(lala, True);
	    densitymatrices[0].reference( lala.reform(IPosition(2, nx, ny)));
	  }
	  imwgt_p.setWeightDensity( densitymatrices );
	}

	vi_p->useImagingWeight(imwgt_p);
	itsMappers.releaseImageLocks();

      }
    catch (AipsError &x)
      {
	throw(AipsError("In setWeightDensity : "+x.getMesg()));
      }
    return true;
  }
  
  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  void SynthesisImagerVi2::createVisSet(const Bool /*writeAccess*/)
  {
    LogIO os( LogOrigin("SynthesisImagerVi2","createVisSet",WHERE) );
    cerr << "mss_p num" << mss_p.nelements() <<  " sel  " << fselections_p->size() << endl;
    lockMSs();
    if(mss_p.nelements() > uInt(fselections_p->size()) && (fselections_p->size() !=0)){
      throw(AipsError("Discrepancy between Number of MSs and Frequency selections"));
    }

    Block<Int> col;
    col.resize(4);
    col[0] = MS::ARRAY_ID;
    col[1] = MS::DATA_DESC_ID;
    col[2] = MS::FIELD_ID;
    col[3] = MS::TIME;
    vi::SortColumns sc(col, false);
  vi_p = new VisibilityIterator2(mss_p, sc, true); //writeAccess);

    if(fselections_p->size() !=0){
      CountedPtr<vi::FrequencySelections> tmpfselections=new FrequencySelections();
      //Temporary fix till we get rid of old vi and we can get rid of tuneSelect
      if(uInt(fselections_p->size()) > mss_p.nelements()){
	for(uInt k=0 ; k <  mss_p.nelements(); ++k){
	  tmpfselections->add(fselections_p->get(k));
	}
      }
      else{
	tmpfselections=fselections_p;
      }
      ////end of fix for tuneSelectdata 
      vi_p->setFrequencySelection (*tmpfselections);

    }
    //
    vi_p->originChunks();
    vi_p->origin();
    ////make sure to use the latest imaging weight scheme
    vi_p->useImagingWeight(imwgt_p);
    unlockMSs();
  }// end of createVisSet

  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  // Method to run the AWProjectFT machine to seperate the CFCache
  // construction from imaging.  This is done by splitting the
  // operation in two steps: (1) run the FTM in "dry" mode to create a
  // "blank" CFCache, and (2) re-load the "blank" CFCache and "fill"
  // it.
  //
  // If someone can get me (SB) out of the horrible statc_casts in the
  // code below, I will be most grateful (we are out of it! :-)).
  //
  void SynthesisImagerVi2::dryGridding(const Vector<String>& cfList)
  {
    LogIO os( LogOrigin("SynthesisImagerVi2","dryGridding",WHERE) );

    rownr_t cohDone=0;
    Int whichFTM=0;
    (void)cfList;
    // If not an AWProject-class FTM, make this call a NoOp.  Might be
    // useful to extend it to other projection FTMs -- but later.
    String ftmName = ((*(itsMappers.getFTM2(whichFTM)))).name();

    if (!((itsMappers.getFTM2(whichFTM,true))->isUsingCFCache())) return;

    os << "---------------------------------------------------- Dry Gridding ---------------------------------------------" << LogIO::POST;

    //
    // Go through the entire MS in "dry" mode to set up a "blank"
    // CFCache.  This is done by setting the AWPWBFT in dryrun mode
    // and gridding.  The process of gridding emits CFCache, which
    // will be "blank" in a dry run.
    {
      vi::VisBuffer2* vb=vi_p->getVisBuffer();
      vi_p->originChunks();
      vi_p->origin();
      Double numberCoh=0;
      for (uInt k=0; k< mss_p.nelements(); ++k)
	numberCoh+=Double(mss_p[k]->nrow());

      ProgressMeter pm(1.0, numberCoh, "dryGridding", "","","",true);

      itsMappers.initializeGrid(*vi_p);
    
      // Set the gridder (iFTM) to run in dry-gridding mode
      (itsMappers.getFTM2(whichFTM,true))->setDryRun(true);

      Bool aTermIsOff=False;
      {
	CountedPtr<refim::FTMachine> ftm=itsMappers.getFTM2(whichFTM,True);
	CountedPtr<refim::ConvolutionFunction> cf=ftm->getAWConvFunc();
	aTermIsOff = cf->getTerm("ATerm")->isNoOp();
      }

      os << "Making a \"blank\" CFCache"
	 << (aTermIsOff?" (without the A-Term)":"")
	 << LogIO::WARN << LogIO::POST;

      // Step through the MS.  This triggers the logic in the Gridder
      // to determine all the CFs that will be required.  These empty
      // CFs are written to the CFCache which can then be filled via
      // a call to fillCFCache().
      for (vi_p->originChunks(); vi_p->moreChunks();vi_p->nextChunk())
	{
	  for (vi_p->origin(); vi_p->more(); vi_p->next())
	    {
	      if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS) 
		{
		  itsMappers.grid(*vb, true, refim::FTMachine::OBSERVED, whichFTM);
		  cohDone += vb->nRows();
		  pm.update(Double(cohDone));
		  // If there is no term that depends on time, don't iterate over the entire data base
		  if (aTermIsOff) break;
		}
	    }
	}
    }
    if (cohDone == 0) os << "No valid rows found in dryGridding." << LogIO::EXCEPTION << LogIO::POST;
    // Unset the dry-gridding mode.
    (itsMappers.getFTM2(whichFTM,true))->setDryRun(false);

    //itsMappers.checkOverlappingModels("restore");
    unlockMSs();
    //fillCFCache(cfList);
  }
  //
  // Re-load the CFCache from the disk using the supplied list of CFs
  // (as cfList).  Then extract the ConvFunc object (which was setup
  // in the FTM) and call it's makeConvFunction2() to fill the CF.
  // Finally, unset the dry-run mode of the FTM.
  //
  void SynthesisImagerVi2::fillCFCache(const Vector<String>& cfList,
				       const String& ftmName,
				       const String& cfcPath,
				       const Bool& psTermOn,
				       const Bool& aTermOn,
				       const Bool& conjBeams)
    {
      LogIO os( LogOrigin("SynthesisImagerVi2","fillCFCache",WHERE) );
      // If not an AWProject-class FTM, make this call a NoOp.  Might be
      // useful to extend it to other projection FTMs -- but later.
      // String ftmName = ((*(itsMappers.getFTM(whichFTM)))).name();

      if ( !(ftmName.at(0,3)=="awp") &&
	  !ftmName.contains("multitermftnew")) return;
      //if (!ftmName.contains("awproject")) return;
      
      os << "---------------------------------------------------- fillCFCache ---------------------------------------------" << LogIO::POST;

      //String cfcPath = itsMappers.getFTM(whichFTM)->getCacheDir();
      //String imageNamePrefix=itsMappers.getFTM(whichFTM)->getCFCache()->getWtImagePrefix();

      //cerr << "Path = " << path << endl;

      // CountedPtr<AWProjectWBFTNew> tmpFT = new AWProjectWBFTNew(static_cast<AWProjectWBFTNew &> (*(itsMappers.getFTM(whichFTM))));
      cerr << "@@@@IN making CFCache" << endl; 

      Float dPA=360.0,selectedPA=2*360.0;
      if (cfList.nelements() > 0)
	{
	  CountedPtr<refim::CFCache> cfCacheObj = new refim::CFCache();
	  //Vector<String> wtCFList; wtCFList.resize(cfList.nelements());
	  //for (Int i=0; i<wtCFList.nelements(); i++) wtCFList[i] = "WT"+cfList[i];
	  //Directory dir(path);
	  Vector<String> cfList_p=cfList;//dir.find(Regex(Regex::fromPattern("CFS*")));
	  Vector<String> wtCFList_p;
	  wtCFList_p.resize(cfList_p.nelements());
	  for (Int i=0; i<(Int)wtCFList_p.nelements(); i++) wtCFList_p[i]="WT"+cfList_p[i];

	  //cerr << cfList_p << endl;
      	  cfCacheObj->setCacheDir(cfcPath.data());

	  os << "Re-loading the \"blank\" CFCache for filling" << LogIO::WARN << LogIO::POST;

      	  cfCacheObj->initCacheFromList2(cfcPath, cfList_p, wtCFList_p,
      					 selectedPA, dPA,CFC_VERBOSE);

	  // tmpFT->setCFCache(cfCacheObj);
	  Vector<Double> uvScale, uvOffset;
	  Matrix<Double> vbFreqSelection;
	  CountedPtr<refim::CFStore2> cfs2 = CountedPtr<refim::CFStore2>(&cfCacheObj->memCache2_p[0],false);//new CFStore2;
	  CountedPtr<refim::CFStore2> cfwts2 =  CountedPtr<refim::CFStore2>(&cfCacheObj->memCacheWt2_p[0],false);//new CFStore2;

	  //
	  // Get whichFTM from itsMappers (SIMapperCollection) and
	  // cast it as AWProjectWBFTNew.  Then get the ConvFunc from
	  // the FTM and cast it as AWConvFunc.  Finally call
	  // AWConvFunc::makeConvFunction2().
	  //
	  // (static_cast<AWConvFunc &> 
	  //  (*(static_cast<AWProjectWBFTNew &> (*(itsMappers.getFTM(whichFTM)))).getAWConvFunc())
	  //  ).makeConvFunction2(String(path), uvScale, uvOffset, vbFreqSelection,
	  // 		       *cfs2, *cfwts2);

	  // This is a global methond in AWConvFunc.  Does not require
	  // FTM to be constructed (which can be expensive in terms of
	  // memory footprint).
	  refim::AWConvFunc::makeConvFunction2(String(cfcPath), uvScale, uvOffset, vbFreqSelection,
					       *cfs2, *cfwts2, psTermOn, aTermOn, conjBeams);
      	}
      //cerr << "Mem used = " << itsMappers.getFTM(whichFTM)->getCFCache()->memCache2_p[0].memUsage() << endl;
      //(static_cast<AWProjectWBFTNew &> (*(itsMappers.getFTM(whichFTM)))).getCFCache()->initCache2();
    }

  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  void SynthesisImagerVi2::reloadCFCache()
  {
      LogIO os( LogOrigin("SynthesisImagerVi2","reloadCFCache",WHERE) );
      Int whichFTM=0; 
      CountedPtr<refim::FTMachine> ftm=itsMappers.getFTM2(whichFTM,true);

      // Proceed only if FMTs uses the CFCache mechanism. The first FTM
      // in the Mapper is used to make this decision.  Not sure if the
      // framework pipes allow other FTMs in SIMapper to be
      // fundamentally different. If it does, and if that is
      // triggered, the current decision may be insufficient.
      if (!(ftm->isUsingCFCache())) return; // Better check than checking against FTM name
      
      os << "-------------------------------------------- Re-load CFCache ---------------------------------------------" << LogIO::POST;

      // Following code that distinguishes between MultiTermFTM and
      // all others should ideally be replaced with a polymorphic
      // solution.  I.e. all FTMs should have a working getFTM2(bool)
      // method.  This is required since MultiTermFTM is a container
      // FTM and it's getFTM2() returns the internal (per-MTMFS term)
      // FTMs.  Non-container FTMs must return a pointer to
      // themselves.  The if-else below is because attempts to make
      // AWProjectFT::getFTM2() work have failed.
      //
      // Control reaches this stage only if the isUsingCFCache() test
      // above return True.  The only FTMs what will pass that test
      // for now are AWProjectFT (and its derivatives) and
      // MultiTermFTM if it is constructed with AWP.
      //
      CountedPtr<refim::CFCache> cfc;
      if (ftm->name().contains("MultiTerm")) cfc = ftm->getFTM2(true)->getCFCache();
      else                                   cfc = ftm->getCFCache();
      cfc->setLazyFill((refim::SynthesisUtils::getenv("CFCache.LAZYFILL",1)==1));
      cfc->initCache2();


      // String path,imageNamePrefix;
      // if (ftm->name().contains("MultiTerm"))
      // 	{
      // 	  path = ftm->getFTM2(true)->getCacheDir();
      // 	  imageNamePrefix = ftm->getFTM2(true)->getCFCache()->getWtImagePrefix();
      // 	}
      // else
      // 	{
      // 	  path = ftm->getCacheDir();
      // 	  imageNamePrefix = ftm->getCFCache()->getWtImagePrefix();
      // 	}
	

      // CountedPtr<refim::CFCache> cfCacheObj = new refim::CFCache();
      // cfCacheObj->setCacheDir(path.c_str());
      // cfCacheObj->setWtImagePrefix(imageNamePrefix.c_str());
      // cfCacheObj->setLazyFill((refim::SynthesisUtils::getenv("CFCache.LAZYFILL",1)==1));
      // cfCacheObj->initCache2();

      // // This assumes the itsMappers is always SIMapperCollection.
      // for (whichFTM = 0; whichFTM < itsMappers.nMappers(); whichFTM++)
      // 	{
      // 	  CountedPtr<refim::FTMachine> ifftm=itsMappers.getFTM2(whichFTM,true),
      // 	    fftm=itsMappers.getFTM2(whichFTM,false);
	
      // 	  ifftm->setCFCache(cfCacheObj,true);
      // 	  fftm->setCFCache(cfCacheObj,true);
      // 	}
  }
  //////////////////

   bool  SynthesisImagerVi2::makeMosaicSensitivity(){
     ///We will bother with the first image. As A projection style gridding
     ///usually is done on that first image.
     /// if necessary in the future we will need to migrate this to SIMapper to
     /// do it for all fields if multiple fields are A-projected. 
     if(!itsMappers.getFTM2(0))
       return False;
     /////////////////
       /////////////////
     if(gridparsVec_p[0].ftmachine=="awphpg"){
       //reset the ftmachine as it keeps dying on second run
       resetAWPHPG();
     }

          

    vi::VisBuffer2* vb=vi_p->getVisBuffer();
     vi_p->originChunks();
     vi_p->origin();
     Double numcoh=0;
      for (uInt k=0; k< mss_p.nelements(); ++k)
	numcoh+=Double(mss_p[k]->nrow());
      ProgressMeter pm(1.0, numcoh, 
                          "Gridding Weights for PB", "","","",true);
      rownr_t cohDone=0;
      

      ///This will initialize weight grid too.
      itsMappers.initializeGrid(*vi_p,True);
      //These 2 lines are for AWProj ftmachines
      itsMappers.getFTM2(0)->setPBReady(false);
      itsMappers.getFTM2(0)->setFTMType(casa::refim::FTMachine::WEIGHT);


      for (vi_p->originChunks(); vi_p->moreChunks();vi_p->nextChunk())
    	{
          
	  for (vi_p->origin(); vi_p->more(); vi_p->next())
            {
              if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
		    {
                      itsMappers.getFTM2(0)->gridImgWeights(*vb);//This just calls AWP::put();

                      cohDone += vb->nRows();
                      pm.update(Double(cohDone));
		    }
    		}
    	}
      //now load the images in weight and sumwt
      itsMappers.getFTM2(0)-> finalizeToWeightImage(*vb, imageStore(0));  
      itsMappers.getFTM2(0)->setPBReady(false);

      //cerr << "@@@@@@@MAKING PB " << endl;
      return True;
     

   }
  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  bool  SynthesisImagerVi2::loadMosaicSensitivity(){
    if(!itsMappers.getFTM2(0))
      return False;
    String ftmname=itsMappers.getFTM2(0)->name();
    //cerr << "########Trying to load PB" << endl;
    if(ftmname.contains("Mosaic") || ftmname.contains("AWProjectWB")){
      //sumwt has been calcuated
      Bool donesumwt=(max(itsMappers.imageStore(0)->sumwt()->get()) > 0.0);
      if(donesumwt){
        IPosition shp=itsMappers.imageStore(0)->weight()->shape();
        CoordinateSystem cs=itsMappers.imageStore(0)->weight()->coordinates();
        CountedPtr<ImageInterface<Float> > wgtim=new TempImage<Float>(shp, cs);
        wgtim->copyData(*(itsMappers.imageStore(0)->weight()));
        (const_cast<CountedPtr<refim::FTMachine>& >(itsMappers.getFTM2(0,False)))->setWeightImage(*wgtim);
        (const_cast<CountedPtr<refim::FTMachine>& >(itsMappers.getFTM2(0,True)))->setWeightImage(*wgtim);
        //cerr <<"@@@@@@@@LOADING PB" << endl;
        return true;

      }


    }
    return false;
  }
  /////////////////////////////////////////////////
  Record SynthesisImagerVi2::apparentSensitivity() 
  {
    LogIO os(LogOrigin("SynthesisImagerVi2", "apparentSensitivity()", WHERE));
    
    Record outrec;
    try {

      os << LogIO::NORMAL // Loglevel INFO
	 << "Calculating apparent sensitivity from MS weights, as modified by gridding weight function"
	 << LogIO::POST;
      os << LogIO::NORMAL // Loglevel INFO
	 << "(assuming that MS weights have correct scale and units)"
	 << LogIO::POST;
      
      Double sumNatWt=0.0;
      Double sumGridWt=0.0;
      Double sumGridWt2OverNatWt=0.0;
    
      Float iNatWt(0.0),iGridWt(0.0);
      
      vi::VisBuffer2* vb = vi_p->getVisBuffer();
      vi_p->originChunks();
      vi_p->origin();
      
      // Discover if weightSpectrum non-trivially available
      Bool doWtSp=vi_p->weightSpectrumExists();

      //////
      for(vi_p->originChunks(); vi_p->moreChunks(); vi_p->nextChunk())
	{
	  for (vi_p->origin(); vi_p->more(); vi_p->next())
	    {
	      auto nRow=vb->nRows();
	      const Vector<Bool>& rowFlags(vb->flagRow());

	      const Vector<Int>& a1(vb->antenna1()), a2(vb->antenna2());

              // Extract weights correctly (use WEIGHT_SPECTRUM, if available)
	      Int nCorr=vb->nCorrelations();
              Matrix<Float> wtm;
              Cube<Float> wtc;
              if (doWtSp)
                // WS available [nCorr,nChan,nRow]
                wtc.reference(vb->weightSpectrum());       
              else {
                // WS UNavailable weight()[nCorr,nRow] --> [nCorr,nChan,nRow]
                wtc.reference(vb->weight().reform(IPosition(3,nCorr,1,nRow)));  // unchan'd weight as single-chan
              }
	      Int nChanWt=wtc.shape()(1);  // Might be 1 (no WtSp)

	      Cube<Bool> flagCube(vb->flagCube());
	      for (rownr_t row=0; row<nRow; row++) {
		if (!rowFlags(row) && a1(row)!=a2(row)) {  // exclude ACs

		  for (Int ich=0;ich<vb->nChannels();++ich) {
		    if( !flagCube(0,ich,row) && !flagCube(nCorr-1,ich,row)) {  // p-hands unflagged

		      // Accumulate relevant info

		      // Simple sum of p-hand for now
		      iNatWt=wtc(0,ich%nChanWt,row)+wtc(nCorr-1,ich%nChanWt,row);

		      iGridWt=2.0f*vb->imagingWeight()(ich,row);

		      if (iGridWt>0.0 && iNatWt>0.0) {
			sumNatWt+=(iNatWt);
			sumGridWt+=(iGridWt);
			sumGridWt2OverNatWt+=(iGridWt*iGridWt/iNatWt);
		      }
		    }
		  }
		}
	      } // row
	    } // vb
	} // chunks
      
      if (sumNatWt==0.0) {
	os << "Cannot calculate sensitivity: sum of selected natural weights is zero" << LogIO::EXCEPTION;
      }
      if (sumGridWt==0.0) {
	os << "Cannot calculate sensitivity: sum of gridded weights is zero" << LogIO::EXCEPTION;
      }

      Double effSensitivity = sqrt(sumGridWt2OverNatWt)/sumGridWt;

      Double natSensitivity = 1.0/sqrt(sumNatWt);
      Double relToNat=effSensitivity/natSensitivity;

      os << LogIO::NORMAL << "RMS Point source sensitivity  : " // Loglevel INFO
	 << effSensitivity      //  << " Jy/beam"       // actually, units are arbitrary
	 << LogIO::POST;
      os << LogIO::NORMAL // Loglevel INFO
	 << "Relative to natural weighting : " << relToNat << LogIO::POST;

      // Fill output Record
      outrec.define("relToNat",relToNat);
      outrec.define("effSens",effSensitivity);

    } catch (AipsError x) {
      throw(x);
      return outrec;
    } 
    return outrec;

  }
  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////

  Bool SynthesisImagerVi2::makePB()
  {
      LogIO os( LogOrigin("SynthesisImagerVi2","makePB",WHERE) );

      if( itsMakeVP==False )
	{
          if( ((itsMappers.getFTM2(0))->name())!="MultiTermFTNew")
            if(!loadMosaicSensitivity()){
              if(!makeMosaicSensitivity())
                throw(AipsError("Problem with making/loading sensitivity image for A -projection gridder"));
            }
	}
      else
	{
	  Bool doDefaultVP = itsVpTable.length()>0 ? False : True;

	  CoordinateSystem coordsys=itsMappers.imageStore(0)->getCSys();
	  String telescope=coordsys.obsInfo().telescope();
	  
	  if (doDefaultVP) {
	    
	    MSAntennaColumns ac(mss_p[0]->antenna());
	    Double dishDiam=ac.dishDiameter()(0);
	    if(!allEQ(ac.dishDiameter().getColumn(), dishDiam))
	      os << LogIO::WARN
		 << "The MS has multiple antenna diameters ..PB could be wrong "
		 << LogIO::POST;
	    return makePBImage( telescope, False, dishDiam);
	  }
	  else{
	    return makePBImage(telescope );	
	  }
	  
	}
 
      return False;
  }

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

  Bool SynthesisImagerVi2::makePrimaryBeam(PBMath& pbMath)
  {
    LogIO os( LogOrigin("SynthesisImagerVi2","makePrimaryBeam",WHERE) );

    os << "vi2 : Evaluating Primary Beam model onto image grid(s)" << LogIO::POST;

    itsMappers.initPB();

    vi::VisBuffer2* vb = vi_p->getVisBuffer();
    vi_p->originChunks();
    vi_p->origin();
    std::map<Int, std::set<Int>> fieldsDone;
    VisBufferUtil vbU(*vb);
    ///////if tracking a moving source
    MDirection origMovingDir;
    MDirection newPhaseCenter;
    Bool trackBeam=getMovingDirection(*vb, origMovingDir, True);
    //////
    for(vi_p->originChunks(); vi_p->moreChunks(); vi_p->nextChunk())
      {
	for (vi_p->origin(); vi_p->more(); vi_p->next())
	  {
	    Bool fieldDone=False;
	    if(fieldsDone.count(vb->msId() >0)){
	      fieldDone=fieldDone || (fieldsDone[vb->msId()].count(vb->fieldId()(0)) > 0);
	    }
	    else{
	      fieldsDone[vb->msId()]=std::set<int>();
	    }
	    if(!fieldDone){
	      fieldsDone[vb->msId()].insert(vb->fieldId()(0));
	      if(trackBeam){
		MDirection newMovingDir;
		getMovingDirection(*vb, newMovingDir);
		//newPhaseCenter=vb->phaseCenter();
                newPhaseCenter=vbU.getPhaseCenter(*vb);
		newPhaseCenter.shift(MVDirection(-newMovingDir.getAngle()+origMovingDir.getAngle()), False);
	      }
	      itsMappers.addPB(*vb,pbMath, newPhaseCenter, trackBeam);
	      
	    }
	  }
      }
    itsMappers.releaseImageLocks();
    unlockMSs();

    return True;
  }// end makePB

  Bool SynthesisImagerVi2::getMovingDirection(const vi::VisBuffer2& vb,  MDirection& outDir, const Bool useImageEpoch){
    MDirection movingDir;
    Bool trackBeam=False;
      
    MeasFrame mFrame(MEpoch(Quantity(vb.time()(0), "s"), MSColumns(vb.ms()).timeMeas()(0).getRef()), mLocation_p);
    if(useImageEpoch){
      mFrame.resetEpoch((itsMappers.imageStore(0))->getCSys().obsInfo().obsDate());

    }
    if(movingSource_p != ""){
      MDirection::Types refType;
      trackBeam=True;
      if(Table::isReadable(movingSource_p, False)){
	//seems to be a table so assuming ephemerides table
	Table laTable(movingSource_p);
	Path leSentier(movingSource_p);
	MeasComet laComet(laTable, leSentier.absoluteName());
	movingDir.setRefString("COMET");
	mFrame.set(laComet);
      }
      ///if not a table 
      else  if(casacore::MDirection::getType(refType, movingSource_p)){
	if(refType > casacore::MDirection::N_Types && refType < casacore::MDirection:: N_Planets ){
	  ///A known planet
	  movingDir.setRefString(movingSource_p);
	}
      }
      else if(upcase(movingSource_p)=="TRACKFIELD"){
        VisBufferUtil vbU(vb);
	movingDir=vbU.getEphemDir(vb, MEpoch(mFrame.epoch()).get("s").getValue());
      }
      else{
	throw(AipsError("Erroneous tracking direction set to make pb"));
      }
      MDirection::Ref outref1(MDirection::AZEL, mFrame);
      MDirection tmphazel=MDirection::Convert(movingDir, outref1)();
      MDirection::Ref outref(vb.phaseCenter().getRef().getType(), mFrame);
      outDir=MDirection::Convert(tmphazel, outref)();  
    }
    else{
      outDir=vb.phaseCenter();
      trackBeam=False;
    }
      return trackBeam;


  }
  CountedPtr<vi::VisibilityIterator2> SynthesisImagerVi2::getVi(){
	return vi_p;  
  }
  CountedPtr<refim::FTMachine> SynthesisImagerVi2::getFTM(const Int fid, Bool ift){
	if(itsMappers.nMappers() <=fid)
		throw(AipsError("Mappers have not been set for field "+String::toString(fid)));
	return (itsMappers.getFTM2(fid, ift));
	  
  }
  void SynthesisImagerVi2::updateImageBeamSet(Record& returnRec){
    if(returnRec.isDefined("imageid") && returnRec.asInt("imageid")==0){
      //ImageInfo ii=(itsMappers.imageStore(0)->psf())->imageInfo();
      Vector<Int> chanRange(0);
      if(returnRec.isDefined("chanrange"))
	chanRange=returnRec.asArrayInt("chanrange");
      Int npol=(itsMappers.imageStore(0)->getShape())(2);
      Int nchan=(itsMappers.imageStore(0)->getShape())(3);
      if(chanRange.nelements() ==2 && chanRange(1) < nchan){

	ImageBeamSet iibeamset=cubePsfImageInfo_p.getBeamSet();
	Matrix<GaussianBeam> mbeams=iibeamset.getBeams();
	if(mbeams.nelements()==0){
	  mbeams.resize(itsMappers.imageStore(0)->getShape()(3), itsMappers.imageStore(0)->getShape()(2));
	  mbeams.set(GaussianBeam(Vector<Quantity>(3, Quantity(1e-12, "arcsec"))));
	}
	Cube<Float> recbeams(returnRec.asArrayFloat("beams"));
	for(Int c=chanRange[0]; c <= chanRange[1]; ++c){
	  for(Int p=0; p < npol; ++p){
	    mbeams(c,p)=GaussianBeam(Quantity(recbeams(c-chanRange[0], p, 0),"arcsec"), Quantity(recbeams(c-chanRange[0], p, 1),"arcsec"), Quantity(recbeams(c-chanRange[0], p, 2),"deg"));

	  }
	}
	cubePsfImageInfo_p.setBeams(ImageBeamSet(mbeams));

      }
      //itsMappers.imageStore(0)->psf()->setImageInfo(ii);
      //itsMappers.imageStore(0)->psf()->unlock();

      
    }
  }
  bool SynthesisImagerVi2::inithpg(){
    bool rstat = false;
#ifdef USE_HPG
    try {
      if (!hpg::is_initialized())
        hpg::initialize();
      auto devices = hpg::devices();
      // cerr << "DEvices " << devices << endl;
      rstat = true;
    } catch (...) {
      rstat = false;
      throw(AipsError("Trying to use GPU code with the wrong GPU or no GPU"));
    }
#endif
    return rstat;
  }
  bool SynthesisImagerVi2::hpg_enabled(){
#ifdef USE_HPG
    return true;
#endif
    return false;
  }

} //# NAMESPACE CASA - END