//# FringeJones.cc: Implementation of FringeJones
//# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003,2011
//# 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 Library General Public License as published by
//# the Free Software Foundation; either version 2 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 Library General Public
//# License for more details.
//#
//# You should have received a copy of the GNU Library General Public License
//# along with this library; if not, write to the Free Software Foundation,
//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
//#
//# Correspondence concerning AIPS++ should be addressed as follows:
//#        Internet email: aips2-request@nrao.edu.
//#        Postal address: AIPS++ Project Office
//#                        National Radio Astronomy Observatory
//#                        520 Edgemont Road
//#                        Charlottesville, VA 22903-2475 USA
//#

#include <synthesis/MeasurementComponents/FringeJones.h>
#include <msvis/MSVis/VisBuffer.h>
#include <msvis/MSVis/VisBuffAccumulator.h>
#include <casacore/ms/MeasurementSets/MSColumns.h>
#include <synthesis/CalTables/CTIter.h>
#include <synthesis/MeasurementEquations/VisEquation.h>  // *
#include <synthesis/MeasurementComponents/SolveDataBuffer.h>
#include <synthesis/MeasurementComponents/MSMetaInfoForCal.h>
#include <casacore/lattices/Lattices/ArrayLattice.h>
#include <casacore/lattices/LatticeMath/LatticeFFT.h>
#include <casacore/scimath/Mathematics/FFTServer.h>

#include <casacore/casa/Arrays/ArrayMath.h>
#include <casacore/casa/Arrays/MatrixMath.h>
#include <casacore/casa/Arrays/ArrayLogical.h>
#include <casacore/casa/BasicSL/String.h>
#include <casacore/casa/Utilities/Assert.h>
#include <casacore/casa/Exceptions/Error.h>
#include <casacore/casa/System/Aipsrc.h>

#include <sstream>

#include <casacore/measures/Measures/MCBaseline.h>
#include <casacore/measures/Measures/MDirection.h>
#include <casacore/measures/Measures/MEpoch.h>
#include <casacore/measures/Measures/MeasTable.h>

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

#include <casacore/casa/Arrays/MaskedArray.h>
#include <casacore/casa/Arrays/MaskArrMath.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_spblas.h>
#include <gsl/gsl_multilarge_nlinear.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_linalg.h>
#include <iomanip>                // needed for setprecision

// DEVDEBUG gates the development debugging information to standard
// error; it should be set to 0 for production.
#define DEVDEBUG false
#define KDISPSCALE 1e6

using namespace casa::vi;
using namespace casacore;

namespace casa { //# NAMESPACE CASA - BEGIN

// Start of GSL compliant solver

void my_gsl_error_handler (const char * reason,
                           const char * file,
                           int line,
                           int gsl_errno) {
    throw(AipsError(reason));
} 

  
class AuxParamBundle {
public:
    SDBList &sdbs;    
    size_t nCalls;
private:
    // We make sure there are no copy or default constructors to
    // preserve the integrity of our reference member.
    AuxParamBundle();
    AuxParamBundle(AuxParamBundle const&);
    AuxParamBundle const& operator=(AuxParamBundle const&);

    size_t refant;
    size_t nCorrelations;
    size_t corrStep;
    Double t0;
    Double reftime;
    std::map< Int, std::set< Int > > activeAntennas;
    std::map< Int, Int > antennaIndexMap;
    // Can't I just have a vector, which maps indices to values anyway?
    std::vector<bool> parameterFlags;
    Int nParams;
    std::map< Int, Int > parameterMap;
    Int activeCorr;
public:
    AuxParamBundle(SDBList& sdbs_, size_t refant, const std::map< Int, std::set<Int> >& activeAntennas_, Vector<Bool> paramActive) :
        sdbs(sdbs_),
        nCalls(0),
        refant(refant),
        nCorrelations(sdbs.nCorrelations() > 1 ? 2 : 1),
        corrStep(sdbs.nCorrelations() > 2 ? 3 : 1),
        activeAntennas(activeAntennas_),
        parameterFlags(),
        parameterMap(),
        activeCorr(-1)
        // corrStep(3)
        {
            Int last_index = sdbs.nSDB() - 1 ;
            t0 = sdbs(0).time()(0);
            Double tlast = sdbs(last_index).time()(0);
            reftime = 0.5*(t0 + tlast);

            parameterFlags = paramActive.tovector();
            Int j = 0; // The CASA parameter index (0=peculiar phase, 1=delay, 2=rate, 3=dispersive)
            Int i = 0; // the Least Squares parameter vector index, depending on what's being solved for
            for (auto p=parameterFlags.begin(); p!=parameterFlags.end(); p++) {
                if (*p) {
                    parameterMap.insert(std::pair<Int, Int>(j, i));
                    i++;
                }
                j++;
            }
            if (i==0) {
                throw(AipsError("No parameters specified!"));
            }
            nParams = i; // There's always at least one parameter!
            // cerr << "AuxParamBundle reftime " << reftime << " t0 " << t0 <<" dt " << tlast - t0 << endl;
        }

    Int nParameters() {
        return nParams;
    }
    Double get_t0() {
        return t0;
    }
    Double
    get_ref_time() {
        return reftime;
    }
    size_t
    get_num_corrs() {
        //return sdbs.nCorrelations() > 1 ? 2 : 1;
        return nCorrelations;
    }
    size_t
    get_num_antennas() {
        if (activeCorr < 0) {
            throw(AipsError("Correlation out of range."));
        }
        std::set< Int > ants = activeAntennas.find(activeCorr)->second;
        return (size_t) ants.size();
    }
    size_t
    get_max_antenna_index() {
        if (activeCorr < 0) {
            throw(AipsError("Correlation out of range."));
        }
        return *(activeAntennas.find(activeCorr)->second.rbegin());
    }
    // Sometimes there is Int, sometimes size_t; the following ones are casacore::Int.
    Int
    get_num_data_points() {
        Int nTotalRows = 0;
        for (Int i = 0; i != sdbs.nSDB(); i++) {
            nTotalRows += sdbs(i).nRows();
        }
        return nTotalRows * sdbs.nChannels();
    }
    Int
    get_actual_num_data_points() {
        Int nTotalRows = 0;
        for (Int i = 0; i != sdbs.nSDB(); i++) {
            SolveDataBuffer& s (sdbs(i));
            for (Int irow=0; irow!=s.nRows(); irow++) {
                if (s.flagRow()(irow)) continue;
                nTotalRows++;
            }
        }
        return nTotalRows * sdbs.nChannels();
    }
    size_t 
    get_data_corr_index(size_t icorr) {
        if (icorr > nCorrelations) {
            throw(AipsError("Correlation out of range."));
        }
        size_t dcorr = icorr * corrStep;
        return dcorr;
    }
    bool
    isActive(size_t iant) {
        std::set<Int> ants = activeAntennas.find(activeCorr)->second;
        if (iant == refant) return true;
        else return (ants.find(iant) != ants.end());
    }
    Int
    get_param_corr_param_index(size_t iant0, size_t ipar) {
        if (iant0 == refant) return -1;
        int iant1 = antennaIndexMap[iant0];
        if (iant1 > antennaIndexMap[refant]) {
            iant1 -= 1;
        }
        int ipar1;
        auto p = parameterMap.find(ipar);
        if (p==parameterMap.end()) {
            ipar1 = -1;
        }
        else {
            ipar1 = (iant1 * nParameters()) + p->second;
        }
        return ipar1;
    }
    size_t
    get_active_corr() {
        return activeCorr;
    }
    void
    set_active_corr(size_t icorr) {
        activeCorr = icorr;
        antennaIndexMap.clear();
        Int i = 0;
        std::set<Int>::iterator it;
        std::set<Int> ants = activeAntennas.find(activeCorr)->second;
        for (it = ants.begin(); it != ants.end(); it++) {
            antennaIndexMap[*it] = i++;
        }
    }
};



void
print_baselines(std::set<std::pair< Int, Int > > baselines) {
    cerr << "Baselines encountered ";
    std::set<std::pair< Int, Int > >::iterator it;
    for (it=baselines.begin(); it != baselines.end(); ++it) {
        cerr << "(" << it->first << ", " << it->second << ") ";
    }
    cerr << endl;
}


int
expb_f(const gsl_vector *param, void *d, gsl_vector *f)
{
    AuxParamBundle *bundle = (AuxParamBundle *)d;
    SDBList& sdbs = bundle->sdbs;
    Double refTime = bundle->get_t0();
    
    gsl_vector_set_zero(f);
    //    Vector<Double> freqs = sdbs.freqs();

    const Double reffreq0=sdbs(0).freqs()(0);  // First freq in first SDB
    
    size_t count = 0; // This is the master index.

    Double sumwt = 0.0;
    Double xi_squared = 0.0;

    for (Int ibuf=0; ibuf < sdbs.nSDB(); ibuf++) {
        SolveDataBuffer& s (sdbs(ibuf));
        if (!s.Ok()) continue;

	const Vector<Double> freqs(s.freqs()); // This ibuf's freqs
        Float fmin_ = min(freqs);
        Float fmax = max(freqs);

        const Cube<Complex>& v(s.visCubeCorrected());
        const Cube<Bool>& fl(s.flagCube());
        const Cube<Float>& weights(s.weightSpectrum());
           
        for (Int irow=0; irow!=s.nRows(); irow++) {
            if (s.flagRow()(irow)) continue;

            Int ant1(s.antenna1()(irow));
            Int ant2(s.antenna2()(irow));
            if (!bundle->isActive(ant1) || !bundle->isActive(ant2))
                continue;            
            if (ant1==ant2) continue;

            // VisBuffer.h seems to suggest that a vb.visCube may have shape
            // (nCorr(), nChannel(), nRow())
            size_t icorr0 = bundle->get_active_corr();
            size_t dcorr = bundle->get_data_corr_index(icorr0);
            // We also need to get the right parameters for this,
            // polarization (icorr is an encoding of the
            // polarization of the correlation products).
            
            Double phi0, tau, r, disp;
            {
                Int i;
                Double phi0_1, tau1, r1, disp1;
                Double phi0_2, tau2, r2, disp2;
                
                phi0_1 = ((i = bundle->get_param_corr_param_index(ant1, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
                tau1   = ((i = bundle->get_param_corr_param_index(ant1, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
                r1     = ((i = bundle->get_param_corr_param_index(ant1, 2))>=0) ? gsl_vector_get(param, i) : 0.0;
                disp1  = ((i = bundle->get_param_corr_param_index(ant1, 3))>=0) ? gsl_vector_get(param, i) : 0.0;
                
                phi0_2 = ((i = bundle->get_param_corr_param_index(ant2, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
                tau2   = ((i = bundle->get_param_corr_param_index(ant2, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
                r2     = ((i = bundle->get_param_corr_param_index(ant2, 2))>=0) ? gsl_vector_get(param, i) : 0.0;
                disp2  = ((i = bundle->get_param_corr_param_index(ant2, 3))>=0) ? gsl_vector_get(param, i) : 0.0;
                
                phi0 = phi0_2 - phi0_1;
                tau = tau2 - tau1;
                r = r2-r1;
                disp = disp2-disp1;
            }
            
            for (size_t ichan = 0; ichan != v.ncolumn(); ichan++) {
                if (fl(dcorr, ichan, irow)) continue;

                Float freq = freqs(ichan);
                Float k_disp = KDISPSCALE*C::_2pi*(1.0/freq + (freq-fmin_-fmax)/(fmin_*fmax));

                Complex vis = v(dcorr, ichan, irow);
                Double w0 = weights(dcorr, ichan, irow);
                // FIXME: what should we use to scale the weights?
                // Double weightScale = norm(vis);
                // Double weightScale = abs(vis);
                // Double weightScale = 1;
                // Double weightScale = 1/sqrt(w0); // Actually AIPS 0, not AIPS 1!
                // Double weightScale = pow(w0, -0.75); // AIPS 2
                //  Double weightScale = 1/w0; // AIPS 3
                // Double weightScale = norm(vis); // Casa 1, I guess
                Double w = sqrt(w0);
                sumwt += w*w;
                if (fabs(w) < FLT_EPSILON) continue;
                // We have to turn the delay back into seconds from nanoseconds.
                // Freq difference is in Hz, which comes out typically as 1e6 bands
                //Double wDf = C::_2pi*(freqs(ichan) - freqs(0))*1e-9;
                Double wDf = C::_2pi*(freqs(ichan) - reffreq0)*1e-9;
                //
                Double t1 = s.time()(0);
                // FIXME: Remind me why we *do* scale wDf with 1e-9
                // but do *not* do that with ref_freq?
                // I have a theory which is mine:
                // this is because tau is in nanoseconds.
                //Double ref_freq = freqs(0);
                //Double wDt = C::_2pi*(t1 - refTime) * ref_freq; 
                Double wDt = C::_2pi*(t1 - refTime) * reffreq0; 

                Double mtheta = -(phi0 + tau*wDf + r*wDt + disp*k_disp); 
                Double vtheta = arg(vis);

                Double c_r = w*(cos(mtheta) - cos(vtheta));
                Double c_i = w*(sin(mtheta)  - sin(vtheta));
                gsl_vector_set(f, count, c_r);
                gsl_vector_set(f, count+1, c_i);

                count += 2;
                xi_squared += c_r*c_r + c_i*c_i;
            }
        }
    }
    // cerr << "Residual xi-squared = " << xi_squared << endl;
    return GSL_SUCCESS;
}

    
int
expb_df(CBLAS_TRANSPOSE_t TransJ, const gsl_vector* param, const gsl_vector *u, void *bundle_, gsl_vector *v, gsl_matrix *JTJ)
{

    // param is the current vector for which we're finding the jacobian.
    // if TransJ is true, evaluate J^T u and store in v.
    // Also store J^T . J in lower half of JTJ.
    std::set <std::pair < Int, Int> > baselines;
    AuxParamBundle *bundle = (AuxParamBundle *)bundle_;

    SDBList& sdbs = bundle->sdbs;
    const Double reffreq0=sdbs(0).freqs()(0);  // First freq in first SDB

    size_t count = 0; // This is the master index.

    gsl_vector_set_zero(v);
    gsl_matrix_set_zero(JTJ);
    
    Double refTime = bundle->get_t0();

    for (Int ibuf=0; ibuf < sdbs.nSDB(); ibuf++) {
        SolveDataBuffer& s (sdbs(ibuf));
        if (!s.Ok()) continue;

	const Vector<Double>& freqs(s.freqs()); // This ibuf's freqs
        Float fmin_ = min(freqs);
        Float fmax = max(freqs);
        const Cube<Complex>& vis(s.visCubeCorrected());
        const Cube<Bool>& fl(s.flagCube());
        const Cube<Float>& weights(s.weightSpectrum());

        Double t1 = s.time()(0);
        for (Int irow=0; irow!=s.nRows(); irow++) {
            if (s.flagRow()(irow)) continue;

            Int ant1(s.antenna1()(irow));
            Int ant2(s.antenna2()(irow));
            
            if (ant1==ant2) continue;
            if (!bundle->isActive(ant1) || !bundle->isActive(ant2)) {
                continue;
            }

            // VisBuffer.h seems to suggest that a vb.visCube may have shape
            // (nCorr(), nChannel(), nRow()) 

            size_t icorr0 = bundle->get_active_corr();
            size_t dcorr = bundle->get_data_corr_index(icorr0);
            // We also need to get the right parameters for this
            // polarization (icorr is an encoding of the
            // polarization of the correlation products).
            
            Double phi0, tau, r, disp;
            {
                Int i;
                Double phi0_1, tau1, r1, disp1;
                Double phi0_2, tau2, r2, disp2;
                
                phi0_1 = ((i = bundle->get_param_corr_param_index(ant1, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
                tau1   = ((i = bundle->get_param_corr_param_index(ant1, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
                r1     = ((i = bundle->get_param_corr_param_index(ant1, 2))>=0) ? gsl_vector_get(param, i) : 0.0;
                disp1  = ((i = bundle->get_param_corr_param_index(ant1, 3))>=0) ? gsl_vector_get(param, i) : 0.0;
                
                phi0_2 = ((i = bundle->get_param_corr_param_index(ant2, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
                tau2   = ((i = bundle->get_param_corr_param_index(ant2, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
                r2     = ((i = bundle->get_param_corr_param_index(ant2, 2))>=0) ? gsl_vector_get(param, i) : 0.0;
                disp2  = ((i = bundle->get_param_corr_param_index(ant2, 3))>=0) ? gsl_vector_get(param, i) : 0.0;
                
                phi0 = phi0_2 - phi0_1;
                tau = tau2 - tau1;
                r = r2-r1;
                disp = disp2-disp1;
            }

                
            //Double ref_freq = freqs(0); 
            //Double wDt = C::_2pi*(t1 - refTime) * ref_freq; 
            Double wDt = C::_2pi*(t1 - refTime) * reffreq0; 
            bool found_data = false;

            for (size_t ichan = 0; ichan != vis.ncolumn(); ichan++) {
                if (fl(dcorr, ichan, irow)) continue;
                Double w0 = weights(dcorr, ichan, irow);
                Double w = sqrt(w0);
                if (fabs(w) < FLT_EPSILON) continue;
                found_data = true;

                Float freq = freqs(ichan);
                Float k_disp = KDISPSCALE*C::_2pi*(1.0/freq + (freq-fmin_-fmax)/(fmin_*fmax));
                    
                // Add a 1e-9 factor because tau parameter is in nanoseconds.
                //Double wDf = C::_2pi*(freqs(ichan) - freqs(0))*1e-9;
                Double wDf = C::_2pi*(freqs(ichan) - reffreq0)*1e-9;
                //
                Double mtheta = -(phi0 + tau*wDf + r*wDt + disp*k_disp);
                Double ws = sin(mtheta);
                Double wc = cos(mtheta);

                Double p0 = 1.0;
                Double p1 = wDf;
                Double p2 = wDt;
                Double p3 = k_disp;
                
                Vector<Double> dterm2(4);
                dterm2(0) = -p0;
                dterm2(1) = -p1;
                dterm2(2) = -p2;
                dterm2(3) = -p3;

                Vector<Double> dterm1(4);
                dterm1(0) = p0;
                dterm1(1) = p1;
                dterm1(2) = p2;
                dterm1(3) = p3;
                
                /* 
                   What we want to express is just:
                   J[count + 0, iparam2 + 0] = w*-ws*-1.0; 
                   J[count + 1, iparam2 + 0] = w*+wc*-1.0;
                   J[count + 0, iparam2 + 1] = w*-ws*-wDf;
                   J[count + 1, iparam2 + 1] = w*+wc*-wDf;
                   J[count + 0, iparam2 + 2] = w*-ws*-wDt;
                   J[count + 1, iparam2 + 2] = w*+wc*-wDt;
                   
                   But in the GSL multilarge framework we have to
                   be ready to calculate either J*u for a given u
                   or J^T*u, depending on the flag TransJ, and we also have to fill in the 
                   
                   v[iparam + ...] = J[count + ..., iparam + ...] * u[iparam + ...]
                   
                   or
                   
                   v[iparam + ...] = J^T[iparam + ..., count + ...] * u[count + ...]
                   
                   <https://www.gnu.org/software/gsl/doc/html/nls.html#c.gsl_multifit_nlinear_default_parameters>
                   
                   "Additionally, the normal equations matrix J^T J should be stored in the lower half of JTJ."
                   
                   So we should also use
                   JTJ[iparam + ..., iparam + ...] += J^T[iparam + ..., count + ...] J[count + ..., iparam + ...] 
                   
                */
                if (TransJ==CblasNoTrans) {
                    for (Int di=0; di<4; di++) {
                        Int i;
                        if ((i = bundle->get_param_corr_param_index(ant2, di))>=0) {
                            (*gsl_vector_ptr(v, count + 0)) += (w*-ws*dterm2(di)) * gsl_vector_get(u, i);
                            (*gsl_vector_ptr(v, count + 1)) += (w*+wc*dterm2(di)) * gsl_vector_get(u, i);
                        }
                        if ((i = bundle->get_param_corr_param_index(ant1, di))>=0) {
                            (*gsl_vector_ptr(v, count + 0)) += gsl_vector_get(u, i) * (w*-ws*dterm1(di));
                            (*gsl_vector_ptr(v, count + 1)) += gsl_vector_get(u, i) * (w*+wc*dterm1(di));
                        }
                    }
                } else {
                    for (Int di=0; di<4; di++) {
                        Int i;
                        if ((i = bundle->get_param_corr_param_index(ant2, di))>=0) {
                            (*gsl_vector_ptr(v, i)) += (w*-ws*dterm2(di)) * gsl_vector_get(u, count + 0);
                            (*gsl_vector_ptr(v, i)) += (w*+wc*dterm2(di)) * gsl_vector_get(u, count + 1);
                        }
                        if ((i = bundle->get_param_corr_param_index(ant1, di))>=0) {
                            (*gsl_vector_ptr(v, i)) += gsl_vector_get(u, count + 0) * (w*-ws*dterm1(di));
                            (*gsl_vector_ptr(v, i)) += gsl_vector_get(u, count + 1) * (w*+wc*dterm1(di));
                        }
                    }
                }
                if (JTJ) {
                    Int i, j;
                    Double wterm = (-ws) * (-ws) + (+wc) * (+wc);
                    if (fabs(1-wterm) > 1e-15)
                        throw AipsError("Insufficiently at one");
                    for (Int di=0; di<4; di++) {
                        for (Int dj=0; dj<=di; dj++) {
                            if (((i = bundle->get_param_corr_param_index(ant2, di))>=0) &&
                                ((j = bundle->get_param_corr_param_index(ant2, dj))>=0)) {
                                (*gsl_matrix_ptr(JTJ, i, j)) += w0*dterm2(di)*dterm2(dj);
                            }
                            if (((i = bundle->get_param_corr_param_index(ant1, di))>=0) &&
                                ((j = bundle->get_param_corr_param_index(ant1, dj))>=0)) {
                                (*gsl_matrix_ptr(JTJ, i, j)) += w0*dterm1(di)*dterm1(dj);
                            }
                        }
                    }
                    // iant1 != iant2, so we don't have to worry about collisions
                    for (Int di=0; di<4; di++) {
                        for (Int dj=0; dj<4; dj++) {
                            Int i0, j0;
                            if (((i0 = bundle->get_param_corr_param_index(ant1, di))>=0) &&
                                ((j0 = bundle->get_param_corr_param_index(ant2, dj))>=0)) {
                                Int i1 = max(i0, j0);
                                Int j1 = min(i0, j0);
                                (*gsl_matrix_ptr(JTJ, i1, j1)) += w0*dterm2(di)*dterm1(dj);
                            }
                        }
                    }
                    count += 2;
                } // if JTJ
            } // loop over channels
            if (found_data) {
                std::pair<Int, Int> antpair = std::make_pair(ant1, ant2);
                bool newBaseline = (baselines.find(antpair) == baselines.end());
                if (newBaseline) {
                    baselines.insert(antpair);
                    // cerr << "paramFlagging for antenna "<< ant1 << ": ";
                    // for (size_t di=0; di<4; di++) {
                    //     cerr << (bundle->get_param_corr_param_index(ant1, di)>=0) << " ";
                    // }
                    // cerr << endl;
                    // cerr << "indices for antenna "<< ant1 << ": ";
                    // if (bundle->get_param_corr_param_index(ant1, 0) >= 0) { 
                    //     for (size_t di=0; di<4; di++) {
                    //         cerr << bundle->get_param_corr_param_index(ant1, di) << " ";
                    //     }                    
                    //     cerr << endl;
                    // }
                    // cerr << "phi0 " << phi0 << " tau " << tau << " r " << r << endl;
                }
            }
        }
    }
    if (DEVDEBUG && 0) {
        print_baselines(baselines);
        cerr << "count " << count << endl;
        cerr << "v = ";
        for (size_t i=0; i!=v->size; i++) {
            cerr << gsl_vector_get(v, i) << " ";
        }
        cerr << endl;
        // if (JTJ) {
        //     cerr <<"JTJ " << std::scientific << endl;
        //     for (size_t i=0; i!=JTJ->size1; i++) {
        //         for (size_t j=0; j!=JTJ->size2; j++) {
        //             cerr << gsl_matrix_get(JTJ, i, j) << " ";
        //         }
        //         cerr << endl;
        //     }
        //     cerr << endl;
        // }
    }
    return GSL_SUCCESS;
}
    


int
expb_hess(gsl_vector *param, AuxParamBundle *bundle, gsl_matrix *hess, Double xi_squared, gsl_vector *snr_vector, LogIO& logSink)
{
    // We calculate the diagonal for the hessian as used by AIPS for
    // the signal to noise. The AIPS formulation, by Fred Schwab, is a
    // hand-rolled routine that solves a different problem to ours: by
    // using a triangular matrix for the Jacobian (requiring antenna i<
    // antenna j) the J^T * J term vanishes throughout and the Hessian
    // of the Xi^2 functional *only* includes the second-order
    // derivative terms, which are usually neglected.
    //
    // This is very clever, but it also means different covariance and
    // information matrices, and therefore a different SNR.  Here we
    // use a generic least squares solver but we cheat slightly and use
    // the AIPS form for the Hessian and SNR calculations.
    //
    // FIXME: Is there any compelling reason to use gsl_vectors for
    // this, given that we're not really hooked in to the gsl
    // least squares framework by the time we do this?
    
    SDBList& sdbs = bundle->sdbs;
    Double refTime = bundle->get_t0();

    // Dimensions of (num_antennas); is the same dimension as
    // param vector here.
    gsl_matrix_set_zero(hess);

    const Double reffreq0=sdbs(0).freqs()(0);  // First freq in first SDB

    size_t nobs = 0;
    Double sumwt = 0;
    size_t numpar = param->size;

    for (Int ibuf=0; ibuf < sdbs.nSDB(); ibuf++)
    {
        SolveDataBuffer& s (sdbs(ibuf));
        if (!s.Ok()) continue;

	const Vector<Double> freqs(s.freqs()); // This ibuf's freqs
        Float fmin_ = min(freqs);
        Float fmax = max(freqs);
            
        const Cube<Complex>& v(s.visCubeCorrected());
        const Cube<Bool>& fl(s.flagCube());
        const Cube<Float>& weights(s.weightSpectrum());
           
            
        for (Int irow=0; irow!=s.nRows(); irow++) {
            if (s.flagRow()(irow)) continue;

            Int ant1(s.antenna1()(irow));
            Int ant2(s.antenna2()(irow));
            if (!bundle->isActive(ant1) || !bundle->isActive(ant2))
                continue;            
            if (ant1==ant2) continue;

            // VisBuffer.h seems to suggest that a vb.visCube may have shape
            // (nCorr(), nChannel(), nRow())
            size_t icorr0 = bundle->get_active_corr();
            size_t dcorr = bundle->get_data_corr_index(icorr0);
            // We also need to get the right parameters for this,
            // polarization (icorr is an encoding of the
            // polarization of the correlation products).
            
            Double phi0, tau, r, disp;
            {
                Int i;
                Double phi0_1, tau1, r1, disp1;
                Double phi0_2, tau2, r2, disp2;
                
                phi0_1 = ((i = bundle->get_param_corr_param_index(ant1, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
                tau1   = ((i = bundle->get_param_corr_param_index(ant1, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
                r1     = ((i = bundle->get_param_corr_param_index(ant1, 2))>=0) ? gsl_vector_get(param, i) : 0.0;
                disp1  = ((i = bundle->get_param_corr_param_index(ant1, 3))>=0) ? gsl_vector_get(param, i) : 0.0;

                
                phi0_2 = ((i = bundle->get_param_corr_param_index(ant2, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
                tau2   = ((i = bundle->get_param_corr_param_index(ant2, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
                r2     = ((i = bundle->get_param_corr_param_index(ant2, 2))>=0) ? gsl_vector_get(param, i) : 0.0;  
                disp2  = ((i = bundle->get_param_corr_param_index(ant2, 3))>=0) ? gsl_vector_get(param, i) : 0.0;
              
                phi0 = phi0_2 - phi0_1;
                tau = tau2 - tau1;
                r = r2-r1;
                disp = disp2-disp1;
            }

            for (size_t ichan = 0; ichan != v.ncolumn(); ichan++) {
                if (fl(dcorr, ichan, irow)) continue;
                Complex vis = v(dcorr, ichan, irow);
                // Fixme: this isn't a square root.
                Double w0 = weights(dcorr, ichan, irow);
                sumwt += w0;
                Double w = w0;
                nobs++;
                if (fabs(w) < FLT_EPSILON) continue;
                
                // We have to turn the delay back into seconds from nanoseconds.
                // Freq difference is in Hz, which comes out typically as 1e6 bands
                //Double wDf = C::_2pi*(freqs(ichan) - freqs(0))*1e-9;
                Double wDf = C::_2pi*(freqs(ichan) - reffreq0)*1e-9;
                //
                Double t1 = s.time()(0);

                //Double ref_freq = freqs(0);
                //Double wDt = C::_2pi*(t1 - refTime) * ref_freq;
                Double wDt = C::_2pi*(t1 - refTime) * reffreq0; 

                Float freq = freqs(ichan);
                Float k_disp = KDISPSCALE*C::_2pi*(1.0/freq + (freq-fmin_-fmax)/(fmin_*fmax));

                Double mtheta = -(phi0 + tau*wDf + r*wDt + disp*k_disp); 
                Double vtheta = arg(vis);

                // Hold on a minute though! 
                Double cx = w*cos(vtheta - mtheta);

                Matrix<Double> dterm(4,4);
                dterm(0, 0) = cx;
                dterm(0, 1) = wDf*cx;
                dterm(0, 2) = wDt*cx;
                dterm(0, 3) = k_disp*dterm(0, 1);
                dterm(1, 1) = wDf*dterm(0, 1);
                dterm(1, 2) = wDt*dterm(0, 1);
                dterm(1, 3) = wDf*dterm(0, 3);
                dterm(2, 2) = wDt*dterm(1, 2);
                dterm(2, 3) = wDt*dterm(1, 3);
                dterm(3, 3) = k_disp*dterm(2, 3);

                // Symmetry terms:
                dterm(1, 0) = dterm(0, 1);
                dterm(2, 0) = dterm(0, 2);
                dterm(3, 0) = dterm(0, 3);
                dterm(2, 1) = dterm(1, 2);
                dterm(3, 1) = dterm(1, 3);
                dterm(3, 2) = dterm(2, 3);
                


                for (Int di=0; di<4; di++) {
                    for (Int dj=0; dj<4; dj++) {
                        Int i, j;
                        if (((i = bundle->get_param_corr_param_index(ant1, di))>=0) &&
                            ((j = bundle->get_param_corr_param_index(ant1, dj))>=0)) {
                            *gsl_matrix_ptr(hess, i, j) += dterm(di, dj);
                        }
                        // Exactly the same logic, but with antenna2
                        if (((i = bundle->get_param_corr_param_index(ant2, di))>=0) &&
                            ((j = bundle->get_param_corr_param_index(ant2, dj))>=0)) {
                            *gsl_matrix_ptr(hess, i, j) += dterm(di, dj);
                        }
                    }
                }
                // FIXME: Not just diagonal terms any more!
                // Note that some of these are not in the lower
                // triangular part, even though they are copied
                // faithfully from AIPS which thinks it is filling
                // a triangular matrix and handles symmetry
                // later. Unless I've missed something (again).
                for (Int di=0;  di<4; di++) {
                    for (Int dj=0; dj<4; dj++) {
                        Int i, j;
                        if (((i = bundle->get_param_corr_param_index(ant1, di))>=0) &&
                            ((j = bundle->get_param_corr_param_index(ant2, dj))>=0)) {
                            *gsl_matrix_ptr(hess, i, j) -= dterm(di, dj);
                            *gsl_matrix_ptr(hess, j, i) -= dterm(dj, di);
                        } // if
                    } // dj
                } // di
            } // ichan
        } // irow
    } // ibuff
    // s is more tricky: it is the xi^2 term from exp_f
    
    xi_squared = max(xi_squared, 1e-25);

    if (DEVDEBUG && 0) {
        cerr << "The matrix is" << endl;
        cerr << setprecision(3) << scientific;
        for (size_t i = 0; i != hess->size1; i++) {
            for (size_t j = 0; j < hess->size2; j++) {
                cerr << gsl_matrix_get(hess,i,j) << " ";
            }
            cerr << endl;
        }
        cerr << endl;
        // str.unsetf(cerr:floatfield);
        cerr << std::defaultfloat;
    }
    //
    size_t hsize = hess->size1;
    int signum;
    gsl_permutation *perm = gsl_permutation_alloc (hsize);
    gsl_matrix *lu = gsl_matrix_alloc(hsize, hsize);
    gsl_matrix *inv_hess = gsl_matrix_alloc(hsize, hsize);

    gsl_linalg_LU_decomp(hess, perm, &signum);
    Double det = gsl_linalg_LU_det(hess, signum);
    if ((fabs(det) < GSL_DBL_EPSILON) || std::isnan(det)) {
        logSink << "Hessian matrix singular (determinant=" << det << "); setting signal-to-noise ratio to zero." << LogIO::POST;
       // Singular matrix; fill snrs with zero.
        for (size_t i=0; i < hess->size1; i+=bundle->nParameters()) {
            Double snr = 0;
            gsl_vector_set(snr_vector, i, snr);
        }
    }
    else {
        gsl_linalg_LU_invert(hess, perm, inv_hess);
    
        Double sigma2 = xi_squared / (nobs - numpar) * nobs / sumwt;
        // cerr << "xi_squared " << xi_squared << " Nobs " << nobs << " sumwt " << sumwt << " sigma2 " << sigma2 << endl;
        for (size_t i=0; i < hess->size1; i+=bundle->nParameters()) {
            Double h = gsl_matrix_get(inv_hess, i, i);
            Double snr0 = sqrt(sigma2*h*0.5);
            snr0 = min(snr0, 9999.999);
            Double snr = (snr0 > 1e-20) ? 1.0/snr0 : snr0;
            gsl_vector_set(snr_vector, i, snr);
        }
    }
    gsl_matrix_free(lu);
    gsl_matrix_free(inv_hess);
    gsl_permutation_free(perm);
    // SNR[i], according to aips, is 1/sqrt(sigma2*hess(i1,i1)*0.5);
    // Note that in AIPS i1 ranges over 1..NANT
    // We use 1 as a success code.
    return 1;
}


Int
findRefAntWithData(SDBList& sdbs, Vector<Int>& refAntList, Int prtlev) {
    std::set<Int> activeAntennas;
    for (Int ibuf=0; ibuf != sdbs.nSDB(); ibuf++) {
        SolveDataBuffer& s(sdbs(ibuf));
        if (!s.Ok())
            continue;
        Cube<Bool> fl = s.flagCube();
        for (Int irow=0; irow!=s.nRows(); irow++) {
            if (s.flagRow()(irow))
                continue;
            Int a1(s.antenna1()(irow));
            Int a2(s.antenna2()(irow));
            // Not using irow
            Matrix<Bool> flr = fl.xyPlane(irow);
            if (!allTrue(flr)) {
                activeAntennas.insert(a1);
                activeAntennas.insert(a2);
            }
        }
    }
    if (prtlev > 2) {
        cout << "[FringeJones.cc::findRefAntWithData] refantlist " << refAntList << endl;
        cout << "[FringeJones.cc::findRefAntWithData] activeAntennas: ";
        std::copy(
            activeAntennas.begin(),
            activeAntennas.end(),
            std::ostream_iterator<Int>(std::cout, " ")
            );
        cout << endl;
    }
    Int refAnt = -1;
    for (Vector<Int>::ConstIteratorSTL a = refAntList.begin(); a != refAntList.end(); a++) {
        if (activeAntennas.find(*a) != activeAntennas.end()) {
            if (prtlev > 2)
                cout << "[FringeJones.cc::findRefAntWithData] We are choosing refant " << *a << endl;
            refAnt = *a;
            break;
        } else {
            if (prtlev > 2)
                cout << "[FringeJones.cc::findRefAntWithData] No data for refant " << *a << endl;
        }
    }
    return refAnt;
}


// Stolen from SolveDataBuffer
void
aggregateTimeCentroid(SDBList& sdbs, Int refAnt, std::map<Int, Double>& aggregateTime) {
    // Weighted average of SDBs' timeCentroids
    std::map<Int, Double> aggregateWeight;
    for (Int i=0; i < sdbs.nSDB(); ++i) {
        SolveDataBuffer& s = sdbs(i);
        Vector<Double> wtvD;
        // Sum over correlations and channels to get a vector of weights for each row
        Vector<Float> wtv(partialSums(s.weightSpectrum(), IPosition(2, 0, 1)));
        wtvD.resize(wtv.nelements());
        convertArray(wtvD, wtv);
        for (Int j=0; j<s.nRows(); j++) {
            Int a1 = s.antenna1()(j);
            Int a2 = s.antenna2()(j);
            Int ant;
            if (a1 == refAnt) { ant = a2; }
            else if (a2 == refAnt) { ant = a1; }
            else continue;
            Double w = wtv(j);
            aggregateTime[ant] += w*s.timeCentroid()(j);
            aggregateWeight[ant] += w;
        }
    }
    for (auto it=aggregateTime.begin(); it!=aggregateTime.end(); ++it) {
        Int a = it->first;
        it->second /= aggregateWeight[a];
    }

}


void
print_gsl_vector(gsl_vector *v)
{
    const size_t n = v->size;
    for (size_t i=0; i!=n; i++) {
        cerr << gsl_vector_get(v, i) << " ";
        if (i>0 && (i % 4)==3) cerr << endl;
    }
    cerr << endl;
}

void
print_max_gsl3(gsl_vector *v)
{
    double phi_max = 0.0;
    double del_max = 0.0;
    double rat_max = 0.0;
        
    const size_t n = v->size;
    for (size_t i=0; i!=n/3; i++) {
        if (fabs(gsl_vector_get(v, 3*i+0)) > fabs(phi_max)) phi_max = gsl_vector_get(v, 3*i+0);
        if (fabs(gsl_vector_get(v, 3*i+1)) > fabs(del_max)) del_max = gsl_vector_get(v, 3*i+1);
        if (fabs(gsl_vector_get(v, 3*i+2)) > fabs(rat_max)) rat_max = gsl_vector_get(v, 3*i+2);
    }
    cerr << "phi_max " << phi_max << " del_max " << del_max << " rat_max " << rat_max << endl;
}



/*
gsl-2.4/multilarge_nlinear/fdf.c defines gsl_multilarge_nlinear_driver,
which I have butchered for my purposes here into
not_gsl_multilarge_nlinear_driver(). We still iterate the nonlinear
least squares solver until completion, but we adopt a convergence
criterion copied from AIPS.

Inputs: maxiter  - maximum iterations to allow
        w        - workspace

Additionally I've removed the info parameter, and I may yet regret it.

Originally:
        info     - (output) info flag on why iteration terminated
                   1 = stopped due to small step size ||dx|
                   2 = stopped due to small gradient
                   3 = stopped due to small change in f
                   GSL_ETOLX = ||dx|| has converged to within machine
                               precision (and xtol is too small)
                   GSL_ETOLG = ||g||_inf is smaller than machine
                               precision (gtol is too small)
                   GSL_ETOLF = change in ||f|| is smaller than machine
                               precision (ftol is too small)

Return:
GSL_SUCCESS if converged
GSL_MAXITER if maxiter exceeded without converging
GSL_ENOPROG if no accepted step found on first iteration
*/

int
least_squares_inner_driver (const size_t maxiter,
                            gsl_multilarge_nlinear_workspace * w,
                            AuxParamBundle *bundle)
{
  int status;
  size_t iter = 0;
  /* call user callback function prior to any iterations
   * with initial system state */
  Double s;
  Double last_s = 1.0e30;
  Bool converged = false;
  do  {
      status = gsl_multilarge_nlinear_iterate (w);
      /*
       * If the solver reports no progress on the first iteration,
       * then it didn't find a single step to reduce the
       * cost function and more iterations won't help so return.
       *
       * If we get a no progress flag on subsequent iterations,
       * it means we did find a good step in a previous iteration,
       * so continue iterating since the solver has now reset
       * mu to its initial value.
       */
      if (status == GSL_ENOPROG && iter == 0) {
          return GSL_EMAXITER;
      }

      Double fnorm = gsl_blas_dnrm2(w->f);      
      s = 0.5 * fnorm * fnorm;
      if (DEVDEBUG) {
          cerr << "Iter: " << iter << " ";
          cerr << "Parameters: " << endl;
          print_gsl_vector(w->x);
          print_max_gsl3(w->dx);
          cerr << "1 - s/last_s=" << 1 - s/last_s << endl;
      }
      ++iter;
      if ((1 - s/last_s < 5e-6) && (iter > 1)) converged = true;
      last_s = s;
      /* old test for convergence:
         status = not_gsl_multilarge_nlinear_test(xtol, gtol, ftol, info, w); */
  } while (!converged && iter < maxiter);
  /*
   * the following error codes mean that the solution has converged
   * to within machine precision, so record the error code in info
   * and return success
   */
  if (status == GSL_ETOLF || status == GSL_ETOLX || status == GSL_ETOLG)
  {
      status = GSL_SUCCESS;
  }
  /* check if max iterations reached */
  if (iter >= maxiter && status != GSL_SUCCESS)
      status = GSL_EMAXITER;  return status;
} /* gsl_multilarge_nlinear_driver() */




void
least_squares_driver(SDBList& sdbs, Matrix<Float>& casa_param, Matrix<Bool>& casa_flags, Matrix<Float>& casa_snr, Int refant,
                     const std::map< Int, std::set<Int> >& activeAntennas, Int maxits, Vector<Bool> paramActive, LogIO& logSink) {
    // The variable casa_param is the Casa calibration framework's RParam matrix; we transcribe our results into it only at the end.
    // n below is number of variables,
    // p is number of parameters

    gsl_error_handler_t *old_handler = gsl_set_error_handler (&my_gsl_error_handler);

    // We could pass in an AuxParamBundle instead I guess?
    AuxParamBundle bundle(sdbs, refant, activeAntennas, paramActive);
    for (size_t icor=0; icor != bundle.get_num_corrs(); icor++) {
        bundle.set_active_corr(icor);
        if (bundle.get_num_antennas() == 0) {
            logSink << "No antennas for correlation " << icor << "; not running least-squares solver." << LogIO::POST;
            continue;
        }
        if (bundle.get_num_antennas() == 1) {
            logSink << "No baselines for correlation " << icor << "; not running least-squares solver." << LogIO::POST;
            continue;
        }
        // Four parameters for every antenna, with dispersion
        size_t p = bundle.nParameters() * (bundle.get_num_antennas() - 1);
        // We need to store complex visibilities in a real matrix so we
        // just store real and imaginary components separately.
        size_t n = 2 * bundle.get_actual_num_data_points();

        if (DEVDEBUG) {
            cerr << "bundle.nParameters() " << bundle.nParameters()
                 << " bundle.get_num_antennas() " <<bundle.get_num_antennas() << endl;
            cerr << "p " << p << " n " << n << endl;
        }
        // Parameters for the least-squares solver.
        // param_tol sets roughly the number of decimal places accuracy you want in the answer;
        // I feel that 3 is probably plenty for fringe fitting.
        // param_tol is not used
        //const double param_tol = 1.0e-3;
        // gtol is not used
        // const double gtol = pow(GSL_DBL_EPSILON, 1.0/3.0);
        // ftol is not used
        // const double ftol = 1.0e-20;   
        const size_t max_iter = maxits ;

        const gsl_multilarge_nlinear_type *T = gsl_multilarge_nlinear_trust;
        gsl_multilarge_nlinear_parameters params = gsl_multilarge_nlinear_default_parameters();
        // the Moré scaling is the best equipped to handle very different scales;
        //it should be the best choice to accommodate dispersive terms of O(f)
        params.scale = gsl_multilarge_nlinear_scale_more;
        params.trs = gsl_multilarge_nlinear_trs_lm;
        params.solver = gsl_multilarge_nlinear_solver_cholesky;
        gsl_multilarge_nlinear_workspace *w = gsl_multilarge_nlinear_alloc(T, &params, n, p);
        gsl_multilarge_nlinear_fdf f;

        f.f = &expb_f;
        /* Can't set to NULL for finite-difference Jacobian in multilarge case. */
        f.df =  &expb_df;   
        f.n = n;    /* number of data points */
        f.p = p;    /* number of parameters */
        f.params = &bundle;

        
        // Our original param is a matrix of (3*nCorr, nElem).
        // We have to transcribe it to a vector.

        gsl_vector *gp = gsl_vector_alloc(p);
        gsl_vector_set_zero(gp);

        // We transcribe Casa parameters into gsl vector format, as required by the solver.
        for (size_t iant=0; iant != bundle.get_max_antenna_index()+1; iant++) {
            if (!bundle.isActive(iant)) {
                continue;
            }
            for (int di=0; di<4; di++) {
                Int param_ind = bundle.get_param_corr_param_index(iant, di);
                if (param_ind < 0) continue;
                gsl_vector_set(gp, param_ind, casa_param(4*icor + di, iant));
            }
        }
        gsl_vector *gp_orig = gsl_vector_alloc(p);
        // Keep a copy of original parameters
        gsl_vector_memcpy (gp_orig, gp);
        // initialise workspace
        gsl_multilarge_nlinear_init(gp, &f, w);
    
        // compute initial residual norm */
        gsl_vector *res_f = gsl_multilarge_nlinear_residual(w);

        int info;
        int status = least_squares_inner_driver(max_iter, w, &bundle);
        double chi1 = gsl_blas_dnrm2(res_f);
        
        gsl_vector_sub(gp_orig, w->x);
        gsl_vector *res = gsl_multilarge_nlinear_position(w);
        
        // We transcribe values back from gsl_vector to the param matrix
        
        gsl_matrix *hess = gsl_matrix_alloc(p,p);
        gsl_vector *snr_vector = gsl_vector_alloc(p);
        gsl_matrix_set_zero(hess);
        gsl_vector_set_zero(snr_vector);
        expb_hess(gp, &bundle, hess, chi1*chi1, snr_vector, logSink);
        
        // Transcribe parameters back into CASA arrays
        for (size_t iant=0; iant != bundle.get_max_antenna_index()+1; iant++) {
            if (!bundle.isActive(iant)) continue;
            Int iparam = bundle.get_param_corr_param_index(iant, 0);
            if (iparam<0) {
                continue;
            }
            if (DEVDEBUG) {
                logSink << "iparam " << iparam << endl;
                logSink << "Old values for ant " << iant << " correlation " << icor 
                        << " delay " << casa_param(4*icor + 1, iant) << " ns "
                        << " rate " << casa_param(4*icor + 2, iant)
                        << " angle " << casa_param(4*icor + 0, iant)
                        << endl;
                logSink << "New values for ant " << iant << " correlation " << icor << ":";
                int i;
                if ((i=bundle.get_param_corr_param_index(iant, 0))>=0) {
                    logSink << " Angle " << gsl_vector_get(res, i);
                }
                if ((i=bundle.get_param_corr_param_index(iant, 1))>=0) {
                    logSink << " delay " << gsl_vector_get(res, i) << " ns ";
                }
                if ((i=bundle.get_param_corr_param_index(iant, 2))>=0) {
                    logSink << " rate " << gsl_vector_get(res, i);
                }
                logSink << "." << LogIO::POST;
            }
            if (status==GSL_SUCCESS || status==GSL_EMAXITER) {
                // Current policy is to assume that exceeding max
                // number of iterations is not a deal-breaker, leave it
                // to SNR calculation to decide if the results are
                // useful.
                for (size_t di=0; di<4; di++) {
                    int i=bundle.get_param_corr_param_index(iant, di);
                    int i0 = bundle.get_param_corr_param_index(iant, 0);
                    if (i>=0) {
                        casa_param(4*icor + di, iant) = gsl_vector_get(res, i);
                        casa_snr(4*icor + di, iant) = gsl_vector_get(snr_vector, i0);
                    } else {
                        casa_param(4*icor + di, iant) = 0.0;
                        casa_snr(4*icor + di, iant) = 0.0;
                    }
                }
            } else { // gsl solver failed; flag data
                logSink << "Least-squares solver failed to converge; flagging" << endl;
                for (size_t di=0; di<4; di++) {
                    casa_flags(4*icor + di, iant) = false;
                }
            }
        }

        logSink <<  "Least squares complete for correlation " << icor
                << " after " <<  gsl_multilarge_nlinear_niter(w) << " iterations." << LogIO::POST;

            // << "reason for stopping: " << ((info == 1) ? "small step size" : "small gradient") << endl
            // << "initial |f(x)| = " << chi0 << endl
            // << "final   |f(x)| = " << chi1 << endl
            // << "final step taken = " << diffsize 

        if (DEVDEBUG) {
            cerr << "casa_param " << casa_param << endl;

            switch (info) {
            case 1:
                logSink << "Small step size." << endl;
                break;
            case 2:
                logSink << "Flatness." << endl;
            }
            logSink << LogIO::POST;
        }
	gsl_set_error_handler(old_handler);
        gsl_vector_free(gp);
        gsl_vector_free(gp_orig);
        gsl_vector_free(snr_vector);
        gsl_matrix_free(hess);
        gsl_multilarge_nlinear_free(w);
    }
}

    


// **********************************************************
//  CTRateAwareTimeInterp1 Implementations
//

CTRateAwareTimeInterp1::CTRateAwareTimeInterp1(NewCalTable& ct,
					       const casacore::String& timetype,
					       casacore::Array<Float>& result,
					       casacore::Array<Bool>& rflag) :
  CTTimeInterp1(ct,timetype,result,rflag)
{}

// Destructor (nothing to do locally)
CTRateAwareTimeInterp1::~CTRateAwareTimeInterp1() {}

Bool CTRateAwareTimeInterp1::interpolate(Double newtime) {
  // Call generic first
  if (CTTimeInterp1::interpolate(newtime)) {
    // Only if generic yields new calibration
    // NB: lastWasExact_=exact in generic
    applyPhaseRate(timeType().contains("nearest") || lastWasExact_);
    return true;
  }
  else
    // No change
    return false;

}

// Do the phase rate math
void CTRateAwareTimeInterp1::applyPhaseRate(Bool single)
{

  Int ispw=mcols_p->spwId()(0);  // should only be one (sliced ct_)!
  MSSpectralWindow msSpw(ct_.spectralWindow());
  MSSpWindowColumns msCol(msSpw);
  //Vector<Double> refFreqs;
  //msCol.refFrequency().getColumn(refFreqs,True);

  Vector<Double> freqs;
  msCol.chanFreq().get(ispw,freqs,True);  // should only be 1
  Double centroidFreq=freqs(0);

  // cout << "time = " << (currTime_ - timeRef_) << endl;

  if (single) {
    for (Int ipol=0;ipol<2;ipol++) {
      Double dtime=(currTime_-timeRef_)-timelist_(currIdx_);
      Double phase=result_(IPosition(2,ipol*4,0));
      Double rate=result_(IPosition(2,ipol*4+2,0));
      phase+=2.0*C::pi*rate*centroidFreq*dtime;
      result_(IPosition(2,ipol*4,0))=phase;
    }
  } else {
    Vector<uInt> rows(2); indgen(rows); rows+=uInt(currIdx_);
    Cube<Float> r(mcols_p->fparamArray("",rows));

    Vector<Double> dtime(2);
    dtime(0)=(currTime_-timeRef_)-timelist_(currIdx_);
    dtime(1)=(currTime_-timeRef_)-timelist_(currIdx_+1);
    Double wt=dtime(1) / (dtime(1)-dtime(0));


    for (Int ipol=0;ipol<2;ipol++) {
      Vector<Double> phase(2), rate(2);
      phase(0)=r.xyPlane(0)(IPosition(2,ipol*4,0));
      phase(1)=r.xyPlane(1)(IPosition(2,ipol*4,0));
      rate(0)=r.xyPlane(0)(IPosition(2,ipol*4+2,0));
      rate(1)=r.xyPlane(1)(IPosition(2,ipol*4+2,0));

      phase(0)+=2.0*C::pi*rate(0)*centroidFreq*dtime(0);
      phase(1)+=2.0*C::pi*rate(1)*centroidFreq*dtime(1);

      Vector<Complex> ph(2);
      ph(0)=Complex(cos(phase(0)),sin(phase(0)));
      ph(1)=Complex(cos(phase(1)),sin(phase(1)));
      ph(0)=Float(wt)*ph(0) + Float(1.0-wt)*ph(1);
      result_(IPosition(2,ipol*4,0))=arg(ph(0));
    }
  }
}




// **********************************************************
//  FringeJones Implementations
//
FringeJones::FringeJones(VisSet& vs) :
    VisCal(vs),             // virtual base
    VisMueller(vs),         // virtual base
    GJones(vs)       // immediate parent
{
    if (prtlev()>2) cout << "FringeJones::FringeJones(vs)" << endl;
}

FringeJones::FringeJones(String msname,Int MSnAnt,Int MSnSpw) :
    VisCal(msname,MSnAnt,MSnSpw),             // virtual base
    VisMueller(msname,MSnAnt,MSnSpw),         // virtual base
    GJones(msname,MSnAnt,MSnSpw)    // immediate parent
{
    if (prtlev()>2) cout << "FringeJones::FringeJones(msname,MSnAnt,MSnSpw)" << endl;
}

FringeJones::FringeJones(const MSMetaInfoForCal& msmc) :
    VisCal(msmc),             // virtual base
    VisMueller(msmc),         // virtual base
    GJones(msmc)    // immediate parent
{
    if (prtlev()>2) cout << "FringeJones::FringeJones(msmc)" << endl;
}

FringeJones::FringeJones(Int nAnt) :
    VisCal(nAnt), 
    VisMueller(nAnt),
    GJones(nAnt)
{
    if (prtlev()>2) cout << "FringeJones::FringeJones(nAnt)" << endl;
}

FringeJones::~FringeJones() {
    if (prtlev()>2) cout << "FringeJones::~FringeJones()" << endl;
}

void FringeJones::setApply(const Record& apply) {
    // Call parent to do conventional things
    GJones::setApply(apply);

    if (calWt()) 
        logSink() << " (" << this->typeName() << ": Enforcing calWt()=false for phase/delay-like terms)" << LogIO::POST;

   // Enforce calWt() = false for delays
    calWt()=false;

    /*
    // Extract per-spw ref Freq for phase(delay) calculation
    //  from the CalTable
    // TBD:  revise as per refFreq decisions
    MSSpectralWindow msSpw(ct_->spectralWindow());
    MSSpWindowColumns msCol(msSpw);
    msCol.refFrequency().getColumn(KrefFreqs_,true);
    KrefFreqs_/=1.0e9;  // in GHz

    /// Re-assign KrefFreq_ according spwmap (if any)
    if (spwMap().nelements()>0) {
        Vector<Double> tmpfreqs;
        tmpfreqs.assign(KrefFreqs_);
        for (uInt ispw=0;ispw<spwMap().nelements();++ispw)
            if (spwMap()(ispw)>-1)
                KrefFreqs_(ispw)=tmpfreqs(spwMap()(ispw));
    }
    */

    // Use the "physical" (centroid) frequency, per spw 
    MSSpectralWindow msSpw(ct_->spectralWindow());
    MSSpWindowColumns msCol(msSpw);
    Vector<Double> tmpfreqs;
    Vector<Double> chanfreq;
    tmpfreqs.resize(msSpw.nrow());
    for (rownr_t ispw=0;ispw<msSpw.nrow();++ispw) {
      if (ispw < msSpw.nrow()) {
	msCol.chanFreq().get(ispw,chanfreq,true);  // reshape, if nec.
	Int nch=chanfreq.nelements();
	tmpfreqs(ispw)=chanfreq(nch/2);
      }
    }

    KrefFreqs_.resize(nSpw()); KrefFreqs_.set(0.0);
    for (rownr_t ispw=0;ispw<(rownr_t)nSpw();++ispw) {
        if (ispw < tmpfreqs.nelements())
            KrefFreqs_(ispw)=tmpfreqs(ispw);
    }

    /// Re-assign KrefFreq_ according spwmap (if any)
    if (spwMap().nelements()>0) {
      for (uInt ispw=0;ispw<spwMap().nelements();++ispw)
	if (spwMap()(ispw)>-1 && ispw < tmpfreqs.nelements())
	  KrefFreqs_(ispw)=tmpfreqs(spwMap()(ispw));
    }

    KrefFreqs_/=1.0e9;  // in GHz
}

void FringeJones::setApply() {
  // This was omitted in copying KJones. It shouldn't have been.
    
  // Call parent to do conventional things
  GJones::setApply();

  // Enforce calWt() = false for delays
  calWt()=false;

  // Set the ref freqs to something usable
  KrefFreqs_.resize(nSpw());
  KrefFreqs_.set(5.0);

}


void FringeJones::setCallib(const Record& callib,
                            const MeasurementSet& selms) {

    // Call parent to do conventional things
    SolvableVisCal::setCallib(callib,selms);

    /*
    if (calWt()) 
        logSink() << " (" << this->typeName() << ": Enforcing calWt()=false for phase/delay-like terms)" << LogIO::POST;
    */
    // Enforce calWt() = false for delays
    calWt()=false;

    /*
    // Extract per-spw ref Freq for phase(delay) calculation
    //  from the CalTable 
   KrefFreqs_.assign(cpp_->refFreqIn());
    KrefFreqs_/=1.0e9;  // in GHz

    // Re-assign KrefFreq_ according spwmap (if any)
    if (spwMap().nelements()>0) {
        Vector<Double> tmpfreqs;
        tmpfreqs.assign(KrefFreqs_);
        for (uInt ispw=0;ispw<spwMap().nelements();++ispw)
            if (spwMap()(ispw)>-1)
                KrefFreqs_(ispw)=tmpfreqs(spwMap()(ispw));
    }
    */

    // Use the "physical" (centroid) frequency, per spw 
    KrefFreqs_.resize(nSpw());
    for (Int ispw=0;ispw<nSpw();++ispw) {
      const Vector<Double>& f(cpp_->freqIn(ispw));
      Int nf=f.nelements();
      KrefFreqs_[ispw]=f[nf/2];  // center (usually this will be same as [0])
    }
    KrefFreqs_/=1.0e9;  // In GHz

    // Re-assign KrefFreq_ according spwmap (if any)
    if (spwMap().nelements()>0) {
      Vector<Double> tmpfreqs;
      tmpfreqs.assign(KrefFreqs_);
      for (uInt ispw=0;ispw<spwMap().nelements();++ispw)
	if (spwMap()(ispw)>-1)
	  KrefFreqs_(ispw)=tmpfreqs(spwMap()(ispw));
    }
    
}

void FringeJones::setSolve(const Record& solve) {

    // Call parent to do conventional things
    GJones::setSolve(solve);
    preavg() = -DBL_MAX;

    // refant isn't properly set until selfSolveOne.  We set it to a
    // known value here so that it can be checked in debugging code.
    refant() = -1;
    if (prtlev() > 2) {
        cout << "Before GJones::setSolve" << endl
             << "FringeJones::setSolve()" <<endl
             << "FringeJones::refant() = "<< refant() <<endl
             << "FringeJones::refantlist() = "<< refantlist() <<endl;
    }
    if (solve.isDefined("zerorates")) {
        zeroRates() = solve.asBool("zerorates");
    }
    if (solve.isDefined("globalsolve")) {
        globalSolve() = solve.asBool("globalsolve");
    }
    if (solve.isDefined("delaywindow")) {
        Array<Double> dw = solve.asArrayDouble("delaywindow");
        delayWindow() = dw;
    } else {
        cerr << "No delay window!" << endl;
    }
    if (solve.isDefined("ratewindow")) {
        rateWindow() = solve.asArrayDouble("ratewindow");
    } else {
        cerr << "No rate window!" << endl;
    }
    if (solve.isDefined("niter")) {
        maxits() = solve.asInt("niter");
    }
    if (solve.isDefined("paramactive")) {
        paramActive() = solve.asArrayBool("paramactive");
    }
    if (solve.isDefined("concatspws")) {
        concatSPWs() = solve.asBool("concatspws");
    }
    if (solve.isDefined("corrcomb")) {
      //cerr << "FringeJones::setsolve() Corrcomb is set! To:"
      //     << solve.asString("corrcomb")
      //     << endl;
        corrcomb() = solve.asString("corrcomb");
    }
}

// Note: this was previously omitted
void FringeJones::specify(const Record& specify) {

  return SolvableVisCal::specify(specify);

}

void FringeJones::calcAllJones() {

  if (prtlev()>6) cout << "       FringeJones::calcAllJones()" << endl;

  // Should handle OK flags in this method, and only
  //  do Jones calc if OK

  Vector<Complex> oneJones;
  Vector<Bool> oneJOK;
  Vector<Float> onePar;
  Vector<Bool> onePOK;

  ArrayIterator<Complex> Jiter(currJElem(),1);
  ArrayIterator<Bool>    JOKiter(currJElemOK(),1);
  ArrayIterator<Float>   Piter(currRPar(),1);
  ArrayIterator<Bool>    POKiter(currParOK(),1);

  Double phase;

  for (Int iant=0; iant<nAnt(); iant++) {
    onePar.reference(Piter.array());
    onePOK.reference(POKiter.array());

    Float fmin_ = min(currFreq());
    Float fmax = max(currFreq());
    for (Int ich=0; ich<nChanMat(); ich++) {
      
      oneJones.reference(Jiter.array());
      oneJOK.reference(JOKiter.array());

      for (Int ipar=0;ipar<nPar();ipar+=4) {
	if (onePOK(ipar)) {
          Float freq = currFreq()(ich);
          Float k_disp = 1e-9*KDISPSCALE*C::_2pi*(1.0/freq + (freq-fmin_-fmax)/(fmin_*fmax));
                          
	  phase=onePar(ipar);
	  phase+=2.0*C::pi*onePar(ipar+1)*
	    (currFreq()(ich)-KrefFreqs_(currSpw()));
	  phase+=2.0*C::pi*onePar(ipar+2)*KrefFreqs_(currSpw())*1e9*
	    (currTime() - refTime());
          Float dph_d = onePar(ipar+3) * k_disp;
          if (DEVDEBUG && 0) {
              cerr << "fmin_ " << fmin_ << " fmax " << fmax << " k_disp " << k_disp
                   << " param " << onePar(ipar+3) << " dph_d " << dph_d << endl;
          }
          phase+=dph_d;
	  oneJones(ipar/4)=Complex(cos(phase),sin(phase));
	  oneJOK(ipar/4)=True;
	} else {
	  oneJOK(ipar/4)=False;
	}
      }
      // Advance iterators
      Jiter.next();
      JOKiter.next();
    }
    // Step to next antenns's pars
    Piter.next();
    POKiter.next();
  }
}


void
FringeJones::calculateSNR(Int nCorr, DelayRateFFT& drf) {
    Matrix<Float> sRP(solveRPar().nonDegenerate(1));
    Matrix<Bool> sPok(solveParOK().nonDegenerate(1));
    Matrix<Float> sSNR(solveParSNR().nonDegenerate(1));
    
    for (size_t icor=0; icor != (size_t) nCorr; icor++) {
        const set<Int>& activeAntennas = drf.getActiveAntennasCorrelation(icor);
        for (Int iant=0; iant != nAnt(); iant++) {
            if (iant == refant()) {
                Double maxsnr = 999.0;
                sSNR(4*icor + 0, iant) = maxsnr;
                sSNR(4*icor + 1, iant) = maxsnr;
                sSNR(4*icor + 2, iant) = maxsnr;
                sSNR(4*icor + 3, iant) = maxsnr;
            }
            else if (activeAntennas.find(iant) != activeAntennas.end()) {
                Double delay = sRP(4*icor + 1, iant);
                Double rate = sRP(4*icor + 2, iant);
                // Note that DelayRateFFT::snr is also used to calculate SNRs for the least square values!
                Float snrval = drf.snr(icor, iant, delay, rate);
                sSNR(4*icor + 0, iant) = snrval;
                sSNR(4*icor + 1, iant) = snrval;
                sSNR(4*icor + 2, iant) = snrval;
                sSNR(4*icor + 3, iant) = snrval;
            } else {
                sPok(4*icor + 0, iant) = false;
                sPok(4*icor + 1, iant) = false;
                sPok(4*icor + 2, iant) = false;
                sPok(4*icor + 3, iant) = false;
            }
        }
    }
}




void
FringeJones::selfSolveOne(SDBList& sdbs) {
    solveRPar()=0.0;
    solveParOK()=false; 
    solveParErr()=1.0; // Does nothing?
    // Maybe we put refFreq, refTime stuff in here?
    Vector<Double> myRefFreqs;
    // Cannot assume we have a calibration table (ct_) in this method.
    // MSSpectralWindow msSpw(ct_->spectralWindow());
    /// MSSpWindowColumns spwCol(msSpw);
    // spwCol.refFrequency().getColumn(myRefFreqs, true);
    //Double ref_freq = myRefFreqs(currSpw());
    //Double ref_freq = sdbs.freqs()(0);
    if (DEVDEBUG) {
        cerr << "Solving for fringes for spw=" << currSpw() << endl
             << "Valid spws are " << freqMetaData_.validSpws() << endl;
    }
    Int spw = currSpw();
    Double reffreq = freqMetaData_.freq(spw)(0); 
    Double t0 = sdbs(0).time()(0);
    Double dt0 = refTime() - t0;
    //Double df0 = ref_freq - sdbs.freqs()(0);
    //Double df0 = 0; 
    Double df0 = reffreq - sdbs(0).freqs()(0);  // global center to global edge

    logSink() << "Solving for fringes for spw=" << currSpw() << " at t="
              << MVTime(refTime()/C::day).string(MVTime::YMD,7)  << LogIO::POST;

    std::map<Int, Double> aggregateTime;
    // Set the refant to the first choice that has data!
    refant() = findRefAntWithData(sdbs, refantlist(), prtlev());
    if (refant()<0) {
        logSink() << "Antennas " << refantlist() << LogIO::POST;
        logSink() << "dt0 " << dt0 
                  << " nSDB() " << sdbs.nSDB()
                  << LogIO::POST; 
        throw(AipsError("No valid reference antenna supplied."));
    }
    else {
        logSink() << "Using reference antenna " << refant() << LogIO::POST;
    }
    aggregateTimeCentroid(sdbs, refant(), aggregateTime);

    if (DEVDEBUG) {
        std::cerr << "Weighted time centroids" << endl; 
        for (auto it=aggregateTime.begin(); it!=aggregateTime.end(); ++it)
            std::cerr << it->first << " => " << it->second - t0 << std::endl;
    }

    // We arrange that we can use either the DelayRateFFT based on concatenation
    // or the new one that combines spectral windows after the FFT
    if (DEVDEBUG) {
        std::cerr << "Making a DelayRateFFT" << endl;
    }
    DelayRateFFT *drfp = DelayRateFFT::makeAChild(concatSPWs(), sdbs, refant(), delayWindow(), rateWindow());
    // DelayRateFFT *drfp = new DelayRateFFTConcat(sdbs, refant(), delayWindow(), rateWindow());
    // DelayRateFFT *drfp = new DelayRateFFTConcat(sdbs, refant(), delayWindow(), rateWindow());
    if (DEVDEBUG) {
        cerr << "Made a DelayRateFFT of some kind!" << endl;
    }
    drfp->FFT();
    drfp->searchPeak();
    Matrix<Float> sRP(solveRPar().nonDegenerate(1));
    Matrix<Bool> sPok(solveParOK().nonDegenerate(1));
    Matrix<Float> sSNR(solveParSNR().nonDegenerate(1));
    logSink() << "sPok " << sPok.shape() << LogIO::POST;
    
    // Map from MS antenna number to index
    // transcribe fft results to sRP
    Int ncol = drfp->param().ncolumn();
    Int nrow = drfp->param().nrow();
    if (DEVDEBUG) {
        std::cerr << "nrow " << nrow << ", ncol " << ncol << endl; 
        std::cerr << "drfp->flag() " << drfp->flag() << endl; 
    }

    for (Int i=0; i!=ncol; i++) { // i==iant
        for (Int j=0; j!=nrow; j++) { // j is parameter number
            Int oj = (j>=3) ? j+1 : j;
            sRP(IPosition(2, oj, i)) = drfp->param()(IPosition(2, j, i));
            sPok(IPosition(2, oj, i)) = !(drfp->flag()(IPosition(2, j, i)));
        }
        // Our estimate for dispersion is zero, unconditionally, and we stand by it.
        sPok(IPosition(2, 3, i)) = true;
        if (nrow > 3)
            sPok(IPosition(2, 7, i)) = true;
    }
    size_t nCorrOrig(sdbs(0).nCorrelations());
    size_t nCorr = (nCorrOrig> 1 ? 2 : 1); // number of p-hands
    calculateSNR(nCorr, *drfp);
    set<Int> belowThreshold;
    Float threshold = minSNR();
    
    for (size_t icor=0; icor != nCorr; icor++) {
        const set<Int>& activeAntennas = drfp->getActiveAntennasCorrelation(icor);
        for (Int iant=0; iant != nAnt(); iant++) {
            if (iant != refant() && (activeAntennas.find(iant) != activeAntennas.end())) {
                Float s = sSNR(4*icor + 0, iant);
		// Start the log message; finished below
		logSink() << "Antenna " << iant << " correlation " << icor << " has (FFT) SNR of " << s;
                if (s < threshold) {
                    belowThreshold.insert(iant);
                    logSink() << " below threshold (" << threshold << ")";
                    // Don't assume these will be flagged later; do it right away.
                    // (The least squares routine will eventually become optional.)
                    sPok(4*icor + 0, iant) = false;
                    sPok(4*icor + 1, iant) = false;
                    sPok(4*icor + 2, iant) = false;
                    sPok(4*icor + 3, iant) = false;
                }
		// Finish the log message
		logSink() << "." << LogIO::POST;
            }
        }
        // We currently remove the antennas below SNR threshold from
        // the object used to handle the FFT fringe search.
        drfp->removeAntennasCorrelation(icor, belowThreshold);
        if (DEVDEBUG) {
            drfp->printActive();
        }
    }
    if (globalSolve()) {
        logSink() << "Starting least squares optimization." << LogIO::POST;
        // Note that least_squares_driver is *not* a method of
        // FringeJones so we pass everything in, including the logSink
        // reference.  Note also that sRP is passed by reference and
        // altered in place.
        least_squares_driver(sdbs, sRP, sPok, sSNR, refant(), drfp->getActiveAntennas(), maxits(),
                             paramActive(), logSink());
    }
    else {
        logSink() << "Skipping least squares optimisation." << LogIO::POST;
    }

    if (DEVDEBUG) {
        cerr << "Ref time " << MVTime(refTime()/C::day).string(MVTime::YMD,7) << endl;
        cerr << "df0 " << df0 << " dt0 " << dt0 << " reffreq*dt0 " << reffreq*dt0 << endl;
        cerr << "reffreq " << reffreq << endl;
        cerr << "df0 " << df0 << " dt0 " << dt0 << " reffreq*dt0 " << reffreq*dt0 << endl;
        cerr << "sRP " << sRP << endl; 
    }

    for (Int iant=0; iant != nAnt(); iant++) {
        for (size_t icor=0; icor != nCorr; icor++) {
            const set<Int>& activeAntennas = drfp->getActiveAntennasCorrelation(icor);
            if (activeAntennas.find(iant) == activeAntennas.end()) {
                continue;
            }
            Double phi0 = sRP(4*icor + 0, iant);
            Double delay = sRP(4*icor + 1, iant);
            Double rate = sRP(4*icor + 2, iant);
            Double k_disp = sRP(4*icor + 3, iant);
            aggregateTime.find(iant);
            // We assume the reference frequency for fringe fitting
            // (which is NOT the one stored in the SPECTRAL_WINDOW
            // table) is the left-hand edge of the frequency grid.
            Double delta1 = df0*delay/1e9;
            Double delta2 = reffreq*dt0*rate;
            Double delta3 = C::_2pi*(delta1+delta2);
            Double dt;
            auto p = aggregateTime.find(iant);
            if (zeroRates() && p != aggregateTime.end()) {
                dt = p->second - t0;
            } else {
                dt = refTime() - t0;
            }
            if (DEVDEBUG) {
                cerr << "Antenna " << iant << ": phi0 " << phi0 << " delay " << delay << " rate " << rate << " dt " << dt << endl
                     << "k_disp " << k_disp << endl
                     << "reffreq "<< reffreq << " Adding corrections for frequency (" << 360*delta1 << ")" 
                     << " and time (" << 360*delta2 << ") degrees." << endl;
            }
            sRP(4*icor + 0, iant) += delta3;
         }
    }
    
    // We can zero the rates here (if needed) whether we did least squares or not.
    if (zeroRates()) {
        logSink() << "Zeroing delay rates in calibration table." << LogIO::POST;
        
        for (size_t icor=0; icor != nCorr; icor++) {
            for (Int iant=0; iant != nAnt(); iant++) {
                sRP(4*icor + 2, iant) = 0.0;
            }
        }
    }
    // Copy the results to the second polarisation for combined pols
    if (corrcomb()=="all") { 
        logSink() << "Correlations combined: Copying results to other correlation" << LogIO::POST;
        for (Int iant=0; iant != nAnt(); iant++) {
            for (Int i=0; i !=4; i++) {
                sRP(4+i, iant) = sRP(i, iant);
                sPok(4+i, iant) = sPok(i, iant);            
                sSNR(4+i, iant) = sSNR(i, iant);
            }
        }
    } 
    if (DEVDEBUG) {
        std::cerr << "sPok " << sPok << endl;
    }
    delete drfp;
}

void
FringeJones::solveOneVB(const VisBuffer&) {
    throw(AipsError("VisBuffer interface not supported!"));
}


void FringeJones::globalPostSolveTinker() {

  // Re-reference the phase, if requested
  if (refantlist()(0)>-1) applyRefAnt();
}

void FringeJones::applyRefAnt() {

  // TBD:
  // 1. Synchronize refant changes on par axis
  // 2. Implement minimum mean deviation algorithm

  if (refantlist()(0)<0) 
    throw(AipsError("No refant specified."));

  Int nUserRefant=refantlist().nelements();

  // Get the preferred refant names from the MS
  String refantName(msmc().antennaName(refantlist()(0)));
  if (nUserRefant>1) {
    refantName+=" (";
    for (Int i=1;i<nUserRefant;++i) {
      refantName+=msmc().antennaName(refantlist()(i));
      if (i<nUserRefant-1) refantName+=",";
    }
    refantName+=")";
  }

  logSink() << "Applying refant: " << refantName
	    << " refantmode = " << refantmode();
  if (refantmode()=="flex")
    logSink() << " (hold alternate refants' phase constant) when refant flagged";
  if (refantmode()=="strict")
    logSink() << " (flag all antennas when refant flagged)";
  logSink() << LogIO::POST;

  // Generate a prioritized refant choice list
  //  The first entry in this list is the user's primary refant,
  //   the second entry is the refant used on the previous interval,
  //   and the rest is a prioritized list of alternate refants,
  //   starting with the user's secondary (if provided) refants,
  //   followed by the rest of the array, in distance order.   This
  //   makes the priorities correct all the time, and prevents
  //   a semi-stochastic alternation (by preferring the last-used
  //   alternate, even if nominally higher-priority refants become
  //   available)


  // Extract antenna positions
  Matrix<Double> xyz;
  if (msName()!="<noms>") {
    MeasurementSet ms(msName());
    MSAntennaColumns msant(ms.antenna());
    msant.position().getColumn(xyz);
  }
  else {
    // TBD RO*
    CTColumns ctcol(*ct_);
    CTAntennaColumns& antcol(ctcol.antenna());
    antcol.position().getColumn(xyz);
  }

  // Calculate (squared) antenna distances, relative
  //  to last preferred antenna
  Vector<Double> dist2(xyz.ncolumn(),0.0);
  for (Int i=0;i<3;++i) {
    Vector<Double> row=xyz.row(i);
    row-=row(refantlist()(nUserRefant-1));
    dist2+=square(row);
  }
  // Move preferred antennas to a large distance
  for (Int i=0;i<nUserRefant;++i)
    dist2(refantlist()(i))=DBL_MAX;

  // Generated sorted index
  Vector<uInt> ord;
  genSort(ord,dist2);

  // Assemble the whole choices list
  Int nchoices=nUserRefant+1+ord.nelements();
  Vector<Int> refantchoices(nchoices,0);
  Vector<Int> r(refantchoices(IPosition(1,nUserRefant+1),IPosition(1,refantchoices.nelements()-1)));
  convertArray(r,ord);

  // set first two to primary preferred refant
  refantchoices(0)=refantchoices(1)=refantlist()(0);

  // set user's secondary refants (if any)
  if (nUserRefant>1) 
    refantchoices(IPosition(1,2),IPosition(1,nUserRefant))=
      refantlist()(IPosition(1,1),IPosition(1,nUserRefant-1));

  //cout << "refantchoices = " << refantchoices << endl;


  if (refantmode()=="strict") {
    nchoices=1;
    refantchoices.resize(1,True);
  }

  Vector<Int> nPol(nSpw(),2);  // TBD:or 1, if data was single pol

  if (nPar()==8) {
    // Verify that 2nd poln has unflagged solutions, PER SPW
    ROCTMainColumns ctmc(*ct_);

    Block<String> cols(1);
    cols[0]="SPECTRAL_WINDOW_ID";
    CTIter ctiter(*ct_,cols);
    Cube<Bool> fl;

    while (!ctiter.pastEnd()) {

      Int ispw=ctiter.thisSpw();
      fl.assign(ctiter.flag());

      IPosition blc(3,0,0,0), trc(fl.shape());
      trc-=1; blc(0)=nPar()/2;
      
      //      cout << "ispw = " << ispw << " nfalse(fl(1,:,:)) = " << nfalse(fl(blc,trc)) << endl;
      
      // If there are no unflagged solutions in 2nd pol, 
      //   avoid it in refant calculations
      if (nfalse(fl(blc,trc))==0)
	nPol(ispw)=1;

      ctiter.next();      
    }
  }
  //  cout << "nPol = " << nPol << endl;

  Bool usedaltrefant(false);
  Int currrefant(refantchoices(0)), lastrefant(-1);

  Block<String> cols(2);
  cols[0]="SPECTRAL_WINDOW_ID";
  cols[1]="TIME";
  CTIter ctiter(*ct_,cols);

  // Arrays to hold per-timestamp solutions
  Cube<Float> solA, solB;
  Cube<Bool> flA, flB;
  Vector<Int> ant1A, ant1B, ant2B;
  Matrix<Complex> refPhsr;  // the reference phasor [npol,nchan] 
  Int lastspw(-1);
  Bool first(true);
  while (!ctiter.pastEnd()) {
    Int ispw=ctiter.thisSpw();
    if (ispw!=lastspw) first=true;  // spw changed, start over

    // Read in the current sol, fl, ant1:
    solB.assign(ctiter.fparam());
    flB.assign(ctiter.flag());
    ant1B.assign(ctiter.antenna1());
    ant2B.assign(ctiter.antenna2()); 

    // First time thru, 'previous' solution same as 'current'
    if (first) {
      solA.reference(solB);
      flA.reference(flB);
      ant1A.reference(ant1B);
    }
    IPosition shB(solB.shape());
    IPosition shA(solA.shape());

    // Find a good refant at this time
    //  A good refant is one that is unflagged in all polarizations
    //     in the current(B) and previous(A) intervals (so they can be connected)
    Int irefA(0),irefB(0);  // index on 3rd axis of solution arrays
    Int ichoice(0);  // index in refantchoicelist
    Bool found(false);
    IPosition blcA(3,0,0,0),trcA(shA),blcB(3,0,0,0),trcB(shB);
    trcA-=1; trcA(0)=trcA(2)=0;
    trcB-=1; trcB(0)=trcB(2)=0;
    ichoice=0;
    while (!found && ichoice<nchoices) { 
      // Find index of current refant choice
      irefA=irefB=0;
      while (ant1A(irefA)!=refantchoices(ichoice) && irefA<shA(2)) ++irefA;
      while (ant1B(irefB)!=refantchoices(ichoice) && irefB<shB(2)) ++irefB;

      if (irefA<shA(2) && irefB<shB(2)) {

	//	cout << " Trial irefA,irefB: " << irefA << "," << irefB 
	//	     << "   Ants=" << ant1A(irefA) << "," << ant1B(irefB) << endl;

	blcA(2)=trcA(2)=irefA;
	blcB(2)=trcB(2)=irefB;
	found=true;  // maybe
	for (Int ipol=0;ipol<nPol(ispw);++ipol) {
	  blcA(0)=blcB(0)=ipol*nPar()/2;
	  trcA(0)=trcB(0)=(ipol+1)*nPar()/2-1;
	  found &= (nfalse(flA(blcA,trcA))>0);  // previous interval
	  found &= (nfalse(flB(blcB,trcB))>0);  // current interval
	} 
      }
      else
	// irefA or irefB out-of-range
	found=false;  // Just to be sure

      if (!found) ++ichoice;  // try next choice next round

    }

    if (found) {
      // at this point, irefA/irefB point to a good refant
      
      // Keep track
      usedaltrefant|=(ichoice>0);
      currrefant=refantchoices(ichoice);
      refantchoices(1)=currrefant;  // 2nd priorty next time

      //      cout << " currrefant = " << currrefant << " (" << ichoice << ")" << endl;

      //      cout << " Final irefA,irefB: " << irefA << "," << irefB 
      //	   << "   Ants=" << ant1A(irefA) << "," << ant1B(irefB) << endl;


      // Only report if using an alternate refant
      if (currrefant!=lastrefant && ichoice>0) {
	logSink() 
	  << "At " 
	  << MVTime(ctiter.thisTime()/C::day).string(MVTime::YMD,7) 
	  << " ("
	  << "Spw=" << ctiter.thisSpw()
	  << ", Fld=" << ctiter.thisField()
	  << ")"
	  << ", using refant " << msmc().antennaName(currrefant)
	  << " (id=" << currrefant 
	  << ")" << " (alternate)"
	  << LogIO::POST;
      }  

      // Form reference phasor [nPar,nChan]
      Matrix<Float> rA,rB;
      Matrix<Bool> rflA,rflB;
      rB.assign(solB.xyPlane(irefB));
      rflB.assign(flB.xyPlane(irefB));
      
      if (!first) {
	// Get and condition previous phasor for the current refant
	rA.assign(solA.xyPlane(irefA));
	rflA.assign(flA.xyPlane(irefA));
	rB-=rA;

	// Accumulate flags
	rflB&=rflA;
      }
      
      //      cout << " rB = " << rB << endl;
      //      cout << boolalpha << " rflB = " << rflB << endl;
      // TBD: fillChanGaps?
      
      // Now apply reference phasor to all antennas
      Matrix<Float> thissol;
      for (Int iant=0;iant<shB(2);++iant) {
	thissol.reference(solB.xyPlane(iant));
	thissol-=rB;
      }
      
      // Set refant, so we can put it back
      ant2B=currrefant;
      
      // put back referenced solutions
      ctiter.setfparam(solB);
      ctiter.setantenna2(ant2B);

      // Next time thru, solB is previous
      solA.reference(solB);
      flA.reference(flB);
      ant1A.reference(ant1B);
      solB.resize();  // (break references)
      flB.resize();
      ant1B.resize();
	
      lastrefant=currrefant;
      first=false;  // avoid first-pass stuff from now on
      
    } // found
    else {
      logSink() 
	<< "At " 
	<< MVTime(ctiter.thisTime()/C::day).string(MVTime::YMD,7) 
	<< " ("
	<< "Spw=" << ctiter.thisSpw()
	<< ", Fld=" << ctiter.thisField()
	<< ")"
	<< ", refant (id=" << currrefant 
	<< ") was flagged; flagging all antennas strictly."
	<< LogIO::POST;
      // Flag all solutions in this interval
      flB.set(True);
      ctiter.setflag(flB);
      if (ichoice == nchoices){
          logSink() << LogIO::WARN
          << "From time: "
          << MVTime(ctiter.thisTime()/C::day).string(MVTime::YMD,7)
          << "in Spw: "
          << ctiter.thisSpw()
          << " refant list exhausted, flagging all solutions"
          << LogIO::POST; 
      }
    }

    // advance to the next interval
    lastspw=ispw;
    ctiter.next();
  }

  if (usedaltrefant)
    logSink() << LogIO::NORMAL
	      << " NB: An alternate refant was used at least once to maintain" << endl
	      << "  phase continuity where the user's preferred refant drops out." << endl
	      << "  Alternate refants are held constant in phase (_not_ zeroed)" << endl
	      << "  during these periods, and the preferred refant may return at" << endl
	      << "  a non-zero phase.  This is generally harmless."
	      << LogIO::POST;

  return;

}

} //# NAMESPACE CASA - END


/*
Example of use in at Python level:

fringefit(vis="n14c2.ms", caltable="fail.fj", field="",spw="1",intent="",
          selectdata=True, timerange="", antenna="", scan="5", observation="",
          msselect="", solint="inf", refant="EF", minsnr=3.0, append=False,
          gaintable=['n14c2.gcal'], parang=False)
*/