~i-martividal/casa-poltools/trunk

#include <Python.h>
#include <numpy/arrayobject.h>
#include <stdio.h>  
#include <sys/types.h>
#include <new>
#include <complex>
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION


/* Docstrings */
static char module_docstring[] =
    "Solver of leakage terms and source polarization.";
static char setData_docstring[] =
    "Reads data pointers and arranges data.";
static char getHessian_docstring[] =
    "Computes the Hessian matrix and residuals vector, given the parameter values";



/* Available functions */
static PyObject *setData(PyObject *self, PyObject *args);
static PyObject *getHessian(PyObject *self, PyObject *args);



/* Module specification */
static PyMethodDef module_methods[] = {
    {"setData", setData, METH_VARARGS, setData_docstring},
    {"getHessian", getHessian, METH_VARARGS, getHessian_docstring},
    {NULL, NULL, 0, NULL}
};


/* Initialize the module */
PyMODINIT_FUNC initPolSolver(void)
{
    PyObject *m = Py_InitModule3("PolSolver", module_methods, module_docstring);
    if (m == NULL)
        return;

/* Load `numpy` functionality. */
    import_array();


}




// Global variables:
static bool DEBUG = false;
bool doCirc, doOrd1;  
typedef std::complex<double> cplx64d;

int Nchan, NCalSour; 
int *Nvis, *NSou; 
int NAnt, NPar;
int NITER = 0;
int **A1, **A2;
double *VARS, *Hessian, *DerVec;
double ***Wgt; // = new double**[1];
double **WgtCorr;

double *TotFlux;
cplx64d **DATA, **COMPS;
cplx64d *DR, *DL, **EPA, **EMA;

/*
TotFlux = new double[1];
NSou = new int[1];
Nvis = new int[1];
DATA = new cplx64d*[1];
COMPS = new cplx64d*[1];
EPA = new cplx64d*[1];
EMA = new cplx64d*[1];
A1 = new int*[1];
A2 = new int*[1];
Wgt = new double**[1];
Wgt[0] = new double*[4];
*/

int *PSou, *PAntR, *PAntL, *VSou, *VAntR, *VAntL;






//////////////////////////////////
// Gets the pointers to data, metadata, Hessian and residuals vector:
static PyObject *setData(PyObject *self, PyObject *args)
{

  PyObject *DATAPy, *A1Py, *A2Py, *COMPSPy, *EPAsPy, *EMAsPy, *PSouPy, *VARSPy, *FluxPy; 
  PyObject *PAntRPy, *PAntLPy, *WgtPy, *WgtCorrPy, *VSouPy, *VAntRPy, *VAntLPy, *HessianPy, *ResVecPy;

  PyObject *Err;
  int i,j;


  if (!PyArg_ParseTuple(args,"OOOOOOOOOOOOOOOOOObb", &DATAPy, &WgtPy, &WgtCorrPy, &A1Py, &A2Py, &COMPSPy, &EPAsPy, &EMAsPy, &PSouPy, &PAntRPy, &PAntLPy, &VSouPy, &VAntRPy, &VAntLPy, &VARSPy, &HessianPy, &ResVecPy,&FluxPy, &doCirc, &doOrd1))
    {printf("Failed setData! Wrong arguments!\n"); fflush(stdout); Err = Py_BuildValue("i",-1); return Err;};

//  delete NSou, Nvis; 
//  delete[] DATA, COMPS, EPA, EMA, A1, A2, Wgt, DR, DL;
  NCalSour = PyList_Size(DATAPy);
  NSou = new int[NCalSour];
  Nvis = new int[NCalSour];
  DATA = new cplx64d*[NCalSour];
  COMPS = new cplx64d*[NCalSour];
  EPA = new cplx64d*[NCalSour];
  EMA = new cplx64d*[NCalSour];
  A1 = new int*[NCalSour];
  A2 = new int*[NCalSour];
  Wgt = new double**[NCalSour];
  WgtCorr = new double *[NCalSour];

  Hessian = (double *)PyArray_DATA(HessianPy);
  DerVec = (double *)PyArray_DATA(ResVecPy);
  NPar =  PyArray_SHAPE(reinterpret_cast<PyArrayObject*>(ResVecPy))[0];

  TotFlux = (double *)PyArray_DATA(FluxPy);


  if(DEBUG){printf("setData variables read. %i sources; %i parameters\n",NCalSour,NPar); fflush(stdout);};

// get pointers to the data:


// Array must be of proper full-polarization data:
  int Ndims = PyArray_NDIM(PyList_GetItem(DATAPy,0));
  if (Ndims != 3){printf("Array must have three dimenstions!\n"); fflush(stdout); Err = Py_BuildValue("i",-2); return Err;};
  NAnt = PyArray_SHAPE(reinterpret_cast<PyArrayObject*>(PAntRPy))[0];
  Nchan = PyArray_SHAPE(reinterpret_cast<PyArrayObject*>(PyList_GetItem(DATAPy,0)))[1];

  if(DEBUG){printf("%i antennas; %i frequency channels\n",NAnt,Nchan); fflush(stdout);};

  for(i=0;i<NCalSour;i++){

    if(DEBUG){printf("Setting data for field %i:\n",i); fflush(stdout);};

    Nvis[i] = PyArray_SHAPE(reinterpret_cast<PyArrayObject*>(PyList_GetItem(DATAPy,i)))[0];
//    NSou[i] = PyArray_SHAPE(reinterpret_cast<PyArrayObject*>(PyList_GetItem(PSouPy,i)))[0];
    NSou[i] = PyArray_SHAPE(reinterpret_cast<PyArrayObject*>(PyList_GetItem(COMPSPy,i)))[2] - 1;

    if(DEBUG){printf("   %i visibilities; %i subcomponents. \n",Nvis[i],NSou[i]); fflush(stdout);};


    DATA[i] = (cplx64d *)PyArray_DATA(PyList_GetItem(DATAPy,i));
    A1[i] = (int *)PyArray_DATA(PyList_GetItem(A1Py,i));
    A2[i] = (int *)PyArray_DATA(PyList_GetItem(A2Py,i));
    COMPS[i] = (cplx64d *)PyArray_DATA(PyList_GetItem(COMPSPy,i));
    EPA[i] = (cplx64d *)PyArray_DATA(PyList_GetItem(EPAsPy,i));
    EMA[i] = (cplx64d *)PyArray_DATA(PyList_GetItem(EMAsPy,i));

    if(DEBUG){printf("Setting weights for field %i\n",i); fflush(stdout);};
    Wgt[i] = new double*[4];
    for(j=0;j<4;j++){
      if(DEBUG){printf("PolProduct %i ",j); fflush(stdout);};
      Wgt[i][j] = (double *)PyArray_DATA(PyList_GetItem(PyList_GetItem(WgtPy,i),j));
      if(DEBUG){printf("Done.\n"); fflush(stdout);};

    };

    WgtCorr[i] = (double *)PyArray_DATA(PyList_GetItem(WgtCorrPy,i));

  };

  if(DEBUG){printf("Getting pointers to parameter indices\n"); fflush(stdout);};

  PSou = (int *)PyArray_DATA(PSouPy);
  PAntR = (int *)PyArray_DATA(PAntRPy);
  PAntL = (int *)PyArray_DATA(PAntLPy);
  VSou = (int *) PyArray_DATA(VSouPy);
  VAntR = (int *)PyArray_DATA(VAntRPy);
  VAntL = (int *)PyArray_DATA(VAntLPy);
  VARS = (double *)PyArray_DATA(VARSPy);


  if(DEBUG){printf("Array pointers set.\n"); fflush(stdout);};


  DR = reinterpret_cast<cplx64d*>(&VARS[VAntR[0]]);
  DL = reinterpret_cast<cplx64d*>(&VARS[VAntL[0]]);

  if(DEBUG){printf("Dterm pointers re-casted and set.\n"); fflush(stdout);};
  

//  DR = new cplx64d[NAnt]; //reinterpret_cast<cplx64d*>(&VARS[VAntR[0]]);
//  DL = new cplx64d[NAnt]; //reinterpret_cast<cplx64d*>(&VARS[VAntL[0]]);
//  for (i=0;i<NAnt;i++){
//    DR[i] = reinterpret_cast<cplx64d>(VARS[VAntR[i]]);
//    DR[i] = reinterpret_cast<cplx64d>(VARS[VAntL[i]]);
//  };

if (DEBUG){
  int k=0;
  for (j=0; j<NCalSour;j++){
    printf("Field id %i has %i visibs with %i channels\n",j,Nvis[j],Nchan);

    printf("Source subcomponents: \n");
    for (i=0; i< NSou[j]; i++){
      printf("#%i / %i:  PSou = %i,  VSou = %i\n",i,k,PSou[k],VSou[k]);
      printf("  Starting values: %.3e, %.3e\n\n",VARS[VSou[k]],VARS[VSou[k]+1]);
      k+=1;
    };
  };

  printf("ANTENNAS: \n");
  for (i=0; i< NAnt; i++){
    printf("#%i:  PAntR = %i,  VAntR = %i\n",i,PAntR[i],VAntR[i]);
    printf("#%i:  PAntL = %i,  VAntL = %i\n",i,PAntL[i],VAntL[i]);
    printf("  Starting values: %.3e, %.3e | %.3e, %.3e\n\n",DR[i].real(),DR[i].imag(),DL[i].real(),DL[i].imag());
  };

  int chi=Nchan/2;
 // j = 0;
  cplx64d RR, RL, LR, LL;
  double wgtaux;
  for (j=0; j<NCalSour;j++){
  printf("FIRST %i VISIBILITIES (FIELD %i; channel %i): \n",2*NAnt,j,chi);

  for (i=0; i < 2*NAnt; i++){
     RR = DATA[j][chi*4 + i*4*Nchan];
     RL = DATA[j][chi*4 + i*4*Nchan+1];
     LR = DATA[j][chi*4 + i*4*Nchan+2];
     LL = DATA[j][chi*4 + i*4*Nchan+3];
     wgtaux = 0.5*(Wgt[j][0][i*Nchan + chi] + Wgt[j][3][i*Nchan + chi]);
     printf("#%i (wgt %.2e): RR = (%-.2e,%-.2e) | RL = (%-.2e,%-.2e) | LR = (%-.2e,%-.2e) | LL = (%-.2e,%-.2e)\n",i,wgtaux,RR.real(),RR.imag(), RL.real(),RL.imag(), LR.real(),LR.imag(), LL.real(),LL.imag());

    };
  };
};



// Return success:
  NITER = 0;
  Err = Py_BuildValue("i",0);

  return Err;


};





/////////////////////////////////////////////////////////////
// Computes Hessian, residuals vector, and Chi Square.

static PyObject *getHessian(PyObject *self, PyObject *args)
{


  PyObject *Err;

  double Ifac, ChiSq, res, F[4];
  double *IRatio, *Iwt;
  bool doDeriv, doResid, doWgt;
  int i,j,k,l, s;
  long Ndata;

  IRatio = new double[NCalSour];
  Iwt = new double[NCalSour];


  if(doCirc){ 
      F[0] = 0.; F[1] = 1.; F[2] = 1.; F[3] = 0.;
  } else {
      F[0] = 1.; F[1] = 1.; F[2] = 1.; F[3] = 1.;
  }; // Relative weights for RR, RL, LR, LL / XX, XY, YX, YY
  
  
  double currWgt[4], ord2, Pf2Abs, PfAbs, ItotAbs;
  
  
  // Turn on/off 2nd order corrections:
  if(doOrd1){ord2=0.0;} else {ord2 = 1.0;};
  

  int VisPar[NPar];
  int currPar;

  int iaux1, iaux2, iaux3, iaux4;
  int printedGood = 0;

  if (!PyArg_ParseTuple(args,"bbb", &doDeriv, &doResid, &doWgt))
    {printf("Failed getHessian! Wrong arguments!\n"); 
     Err = Py_BuildValue("i",-1); fflush(stdout); return Err;};

  if (DEBUG){
     printf("\n Called getHessian. Num. sources: %i, Num. Params: %i\n",NCalSour,NPar);
     Ndata = 0;
     for (s=0; s<NCalSour; s++){
       for (i=0; i<Nvis[s]; i++){
         for(j=0; j<Nchan; j++){
           iaux3 = i*Nchan + j;
           if(Wgt[s][0][iaux3]>0.0 || Wgt[s][3][iaux3]>0.0){Ndata += 1;};
         };
       };
     };
     printf("There are %d good visibilities.\n",Ndata);
  };

  cplx64d Itot, Qtot, Utot;
  cplx64d RRc, RLc, LRc, LLc;
  cplx64d Im = cplx64d(0.,1.);

  cplx64d DTEff;

  cplx64d RR, RL, LR, LL, auxC;
  cplx64d resid[4];
  cplx64d AllDer[NPar][4];

  ChiSq = 0.0;


// set matrix and vector to zero:
  for(i=0;i<NPar;i++){
    DerVec[i] = 0.;
    for(j=0;j<NPar;j++){
      Hessian[i*NPar + j] = 0.;
    };
  };

  
  Ifac = 0.;
  Ndata = 0;
  



  s=0;
  int NtotSou = 0; // Total number of fittable source subcomponents.
  for(i=0; i<NCalSour; i++){

    IRatio[i] = 0.0; Iwt[i] = 0.0; 
    if(DEBUG){printf("Source %i has %i components.\n",i+1,NSou[i]);};

    for(j=0; j<NSou[i]; j++){
      if(PSou[j+s]>=0){NtotSou += 1;};
    };
    s += NSou[i];
  };


  if(DEBUG){printf("There is a total of %i components to fit.\n",NtotSou);};

// Loop over calibrator sources:
  for (s=0; s<NCalSour; s++){


  int Ns = 0;  // Number of subcomponents already computed.
  int Nf = 0;  // Number of fittable subcomponents already computed.

//////////
  for(i=0; i<s; i++){
    for(j=0; j<NSou[i]; j++){
      if(PSou[j+Ns]>=0){Nf += 1;};
    };
    Ns += NSou[i];
  };
/////////



// Loop over all data in proper order (time, channel, polariz):
  for(i=0; i<Nvis[s]; i++){

    for(j=0; j<Nchan; j++){

// Indices in 1D for the multi-dimensional arrays:
     iaux1 = (i*Nchan + j)*(NSou[s]+1);
     iaux2 = 4*(j + i*Nchan);
     iaux3 = i*Nchan + j;

// Index of current parameter in Hessian matrix:
   //  currPar = 2*Nf;
     currPar = 0;

// Reset vector of derivatives to zero:
     for(k=0;k<NPar;k++){for(iaux4=0;iaux4<4;iaux4++){AllDer[k][iaux4] = cplx64d(0.,0.);};};

// Proceed only if data are good:
     if(Wgt[s][0][iaux3]>0.0 || Wgt[s][3][iaux3]>0.0){

      Ndata += 1;
      
// Add Fourier transforms of all source components:
      Itot = COMPS[s][iaux1+NSou[s]]/TotFlux[s]; // cplx64d(0.,0.);
      Qtot = cplx64d(0.,0.);
      Utot = cplx64d(0.,0.);

      for(k=0; k<NSou[s]; k++){
    //    Itot += COMPS[iaux1 + k];
        Qtot += VARS[VSou[k+Ns]]*COMPS[s][iaux1 + k]/TotFlux[s];
        Utot += VARS[VSou[k+Ns]+1]*COMPS[s][iaux1 + k]/TotFlux[s];
        iaux4 = PSou[k+Ns];
        if(doDeriv && iaux4>=0){
    //      printf("%i %i %i\n",iaux4,currPar,NSou[s]);
          VisPar[currPar] = iaux4;
          VisPar[currPar+1] = iaux4+1;
          currPar += 2;
// Derivative w.r.t. Q for all corr products (RR, RL, LR, LL):
          AllDer[iaux4][0] = COMPS[s][iaux1+k]*(std::conj(DR[A2[s][i]])*EPA[s][i] + DR[A1[s][i]]/EPA[s][i]);
          AllDer[iaux4][1] = COMPS[s][iaux1+k]*(EPA[s][i] + ord2*std::conj(DL[A2[s][i]])*DR[A1[s][i]]/EPA[s][i]);
          AllDer[iaux4][2] = COMPS[s][iaux1+k]*(1./EPA[s][i] + ord2*std::conj(DR[A2[s][i]])*DL[A1[s][i]]*EPA[s][i]);
          AllDer[iaux4][3] = COMPS[s][iaux1+k]*(std::conj(DL[A2[s][i]])/EPA[s][i] + DL[A1[s][i]]*EPA[s][i]);
// Derivative w.r.t. U:
          AllDer[iaux4+1][0] = Im*COMPS[s][iaux1+k]*(std::conj(DR[A2[s][i]])*EPA[s][i] - DR[A1[s][i]]/EPA[s][i]);
          AllDer[iaux4+1][1] = Im*COMPS[s][iaux1+k]*(EPA[s][i] - ord2*std::conj(DL[A2[s][i]])*DR[A1[s][i]]/EPA[s][i]);
          AllDer[iaux4+1][2] = Im*COMPS[s][iaux1+k]*(-1./EPA[s][i] + ord2*std::conj(DR[A2[s][i]])*DL[A1[s][i]]*EPA[s][i]);
          AllDer[iaux4+1][3] = Im*COMPS[s][iaux1+k]*(-std::conj(DL[A2[s][i]])/EPA[s][i] + DL[A1[s][i]]*EPA[s][i]);
          
     ////////////////////////////    
     // FOR TESTING:     
     //     AllDer[iaux4][1] = cplx64d((COMPS[iaux1+k]*EPA[i]).real(),0.);
     //     AllDer[iaux4][2] = cplx64d((COMPS[iaux1+k]/EPA[i]).real(),0.);
     //     AllDer[iaux4+1][1] = cplx64d(0.,(COMPS[iaux1+k]*EPA[i]).imag());
     //     AllDer[iaux4+1][2] = cplx64d(0.,-(COMPS[iaux1+k]/EPA[i]).imag());
     ////////////////////////////    

          
          
        };
      };

      Ifac = std::abs(Itot);
      currWgt[0] = Wgt[s][0][iaux3]*Ifac;
      currWgt[1] = Wgt[s][1][iaux3]*Ifac;
      currWgt[2] = Wgt[s][2][iaux3]*Ifac;
      currWgt[3] = Wgt[s][3][iaux3]*Ifac;

      // Stokes I (i.e., (RR + LL)/2 ) from the data, corrected by inverse of Dt matrices:
      Iwt[s] += 0.5*std::abs((DATA[s][iaux2]*(1. + ord2*std::conj(DL[A2[s][i]])*DL[A1[s][i]])/EMA[s][i] + DATA[s][iaux2+3]*EMA[s][i]*(1. + ord2*std::conj(DR[A2[s][i]])*DR[A1[s][i]]) - DATA[s][iaux2+1]*(std::conj(DR[A2[s][i]])/EMA[s][i] + DL[A1[s][i]]*EMA[s][i]) - DATA[s][iaux2+2]*(DR[A1[s][i]]/EMA[s][i] + std::conj(DL[A2[s][i]])*EMA[s][i]))/(1.-ord2*DR[A1[s][i]]*DL[A1[s][i]])/(1.-ord2*std::conj(DR[A2[s][i]])*std::conj(DL[A2[s][i]])));
    
      // Stokes I from the model:
      IRatio[s] += Ifac;
      
      
// Correlation products in antenna frame with no leakage:
      RRc = Itot*EMA[s][i];
      RLc = (Qtot + Im*Utot)*EPA[s][i];
      LRc = (Qtot - Im*Utot)/EPA[s][i];
      LLc = Itot/EMA[s][i];

// Apply leakage to model visibilities:
      RR = (RRc + RLc*std::conj(DR[A2[s][i]]) + LRc*DR[A1[s][i]] + ord2*LLc*std::conj(DR[A2[s][i]])*DR[A1[s][i]]);

      RL = (RLc + RRc*std::conj(DL[A2[s][i]]) + LLc*DR[A1[s][i]] + ord2*LRc*std::conj(DL[A2[s][i]])*DR[A1[s][i]]);
      PfAbs = std::abs(RLc);

      LR = (LRc + RRc*DL[A1[s][i]] + LLc*std::conj(DR[A2[s][i]]) + ord2*RLc*std::conj(DR[A2[s][i]])*DL[A1[s][i]]);
      Pf2Abs = std::abs(LRc);

      ItotAbs = std::abs(Itot);

      if (doWgt && j==0 && PfAbs>0.0 && Pf2Abs > 0.0 && ItotAbs > 0.0){
        WgtCorr[s][i] = std::max(PfAbs/ItotAbs,Pf2Abs/ItotAbs);
      };

      LL = (LLc + LRc*std::conj(DL[A2[s][i]]) + RLc*DL[A1[s][i]] + ord2*RRc*std::conj(DL[A2[s][i]])*DL[A1[s][i]]);

      if(DEBUG && NITER<2 && printedGood<10 && j==Nchan/2){
        printf("Model Noleak %i: RR = (%.2e, %.2e); RL = (%.2e, %.2e); LR = (%.2e, %.2e); LL = (%.2e, %.2e)\n",i,RRc.real(),RRc.imag(),RLc.real(),RLc.imag(),LRc.real(),LRc.imag(),LLc.real(),LLc.imag());
        printf("Model Leak %i: RR = (%.2e, %.2e); RL = (%.2e, %.2e); LR = (%.2e, %.2e); LL = (%.2e, %.2e)\n",i,RR.real(),RR.imag(),RL.real(),RL.imag(),LR.real(),LR.imag(),LL.real(),LL.imag());
        printedGood += 1;
      };

// Residuals for each correlation product:
      resid[0] =  (DATA[s][iaux2  ] - RR);
      resid[1] =  (DATA[s][iaux2+1] - RL); 
      resid[2] =  (DATA[s][iaux2+2] - LR);
      resid[3] =  (DATA[s][iaux2+3] - LL);

// Vector of derivative*residuals (first, for source components):
      for(k=0; k<NSou[s]; k++){
        iaux4 = PSou[k+Ns];
        if(DEBUG){printf("Sou pid (%i): %i  ",s,iaux4);};
        if(doDeriv && iaux4>=0){
          res = 0.0;
          for(l=0;l<4;l++){
            res += AllDer[iaux4][l].real()*resid[l].real()*F[l]*currWgt[l];
            res += AllDer[iaux4][l].imag()*resid[l].imag()*F[l]*currWgt[l];
          }; 
          DerVec[iaux4] += res; res = 0.0;
          for(l=0;l<4;l++){
            res += AllDer[iaux4+1][l].real()*resid[l].real()*F[l]*currWgt[l];
            res += AllDer[iaux4+1][l].imag()*resid[l].imag()*F[l]*currWgt[l];
          }; 
          DerVec[iaux4+1] += res;
        };
      };

// Derivatives w.r.t. Dterms:
      iaux4 = PAntR[A1[s][i]];
   //   currPar = 2*NtotSou;

      if(doDeriv && iaux4>=0){
        VisPar[currPar] = iaux4;
        VisPar[currPar+1] = iaux4+1;
        currPar += 2;
        if(DEBUG){printf("DR1 pid: %i/%i  ",iaux4,currPar);};
        
        AllDer[iaux4][1] = LLc + ord2*LRc*std::conj(DL[A2[s][i]]);
        AllDer[iaux4+1][1] = Im*AllDer[iaux4][1];

        AllDer[iaux4][0] = LRc + ord2*LLc*std::conj(DR[A2[s][i]]);
        AllDer[iaux4+1][0] = Im*AllDer[iaux4][0];

        
        res  = AllDer[iaux4][0].real()*resid[0].real()*F[0]*currWgt[0];
        res += AllDer[iaux4][0].imag()*resid[0].imag()*F[0]*currWgt[0];
        res += AllDer[iaux4][1].real()*resid[1].real()*F[1]*currWgt[1];
        res += AllDer[iaux4][1].imag()*resid[1].imag()*F[1]*currWgt[1];
        DerVec[iaux4] += res;

        res  = AllDer[iaux4+1][0].real()*resid[0].real()*F[0]*currWgt[0];
        res += AllDer[iaux4+1][0].imag()*resid[0].imag()*F[0]*currWgt[0];
        res += AllDer[iaux4+1][1].real()*resid[1].real()*F[1]*currWgt[1];
        res += AllDer[iaux4+1][1].imag()*resid[1].imag()*F[1]*currWgt[1];
        DerVec[iaux4+1] += res;
      };

      iaux4 = PAntL[A1[s][i]];

      if(doDeriv && iaux4>=0){
        VisPar[currPar] = iaux4;
        VisPar[currPar+1] = iaux4+1;
        currPar += 2;
        if(DEBUG){printf("DL1 pid: %i/%i  ",iaux4,currPar);};
        
        AllDer[iaux4][2] = RRc + ord2*RLc*std::conj(DR[A2[s][i]]);
        AllDer[iaux4+1][2] = Im*AllDer[iaux4][2];

        AllDer[iaux4][3] = RLc + ord2*RRc*std::conj(DL[A2[s][i]]);
        AllDer[iaux4+1][3] = Im*AllDer[iaux4][3];

        res  = AllDer[iaux4][2].real()*resid[2].real()*F[2]*currWgt[2];
        res += AllDer[iaux4][2].imag()*resid[2].imag()*F[2]*currWgt[2];
        res += AllDer[iaux4][3].real()*resid[3].real()*F[3]*currWgt[3];
        res += AllDer[iaux4][3].imag()*resid[3].imag()*F[3]*currWgt[3];
        DerVec[iaux4] += res;

        res  = AllDer[iaux4+1][2].real()*resid[2].real()*F[2]*currWgt[2];
        res += AllDer[iaux4+1][2].imag()*resid[2].imag()*F[2]*currWgt[2];
        res += AllDer[iaux4+1][3].real()*resid[3].real()*F[3]*currWgt[3];
        res += AllDer[iaux4+1][3].imag()*resid[3].imag()*F[3]*currWgt[3];
        DerVec[iaux4+1] += res;
      };

      iaux4 = PAntR[A2[s][i]];
      if(doDeriv && iaux4>=0){
        VisPar[currPar] = iaux4;
        VisPar[currPar+1] = iaux4+1;
        currPar += 2;
        if(DEBUG){printf("DR2 pid: %i/%i  ",iaux4,currPar);};
      
        AllDer[iaux4][2] = LLc + ord2*RLc*DL[A1[s][i]];
        AllDer[iaux4+1][2] = -Im*AllDer[iaux4][2];

        
        AllDer[iaux4][0] = RLc + ord2*LLc*DR[A1[s][i]];
        AllDer[iaux4+1][0] = -Im*AllDer[iaux4][0];

        res  = AllDer[iaux4][0].real()*resid[0].real()*F[0]*currWgt[0];
        res += AllDer[iaux4][0].imag()*resid[0].imag()*F[0]*currWgt[0];
        res += AllDer[iaux4][2].real()*resid[2].real()*F[2]*currWgt[2];
        res += AllDer[iaux4][2].imag()*resid[2].imag()*F[2]*currWgt[2];
        DerVec[iaux4] += res;

        res  = AllDer[iaux4+1][0].real()*resid[0].real()*F[0]*currWgt[0];
        res += AllDer[iaux4+1][0].imag()*resid[0].imag()*F[0]*currWgt[0];
        res += AllDer[iaux4+1][2].real()*resid[2].real()*F[2]*currWgt[2];
        res += AllDer[iaux4+1][2].imag()*resid[2].imag()*F[2]*currWgt[2];
        DerVec[iaux4+1] += res;
      };

      iaux4 = PAntL[A2[s][i]];
      if(doDeriv && iaux4>=0){
        VisPar[currPar] = iaux4;
        VisPar[currPar+1] = iaux4+1;
        currPar += 2;
        if(DEBUG){printf("DL2 pid: %i/%i  \n",iaux4,currPar);};

        AllDer[iaux4][1] = RRc + ord2*LRc*DR[A1[s][i]];
        AllDer[iaux4+1][1] = -Im*AllDer[iaux4][1];

        AllDer[iaux4][3] = LRc + ord2*RRc*DL[A1[s][i]];
        AllDer[iaux4+1][3] = -Im*AllDer[iaux4][3];

        res  = AllDer[iaux4][1].real()*resid[1].real()*F[1]*currWgt[1];
        res += AllDer[iaux4][1].imag()*resid[1].imag()*F[1]*currWgt[1];
        res += AllDer[iaux4][3].real()*resid[3].real()*F[3]*currWgt[3];
        res += AllDer[iaux4][3].imag()*resid[3].imag()*F[3]*currWgt[3];
        DerVec[iaux4] += res;

        res  = AllDer[iaux4+1][1].real()*resid[1].real()*F[1]*currWgt[1];
        res += AllDer[iaux4+1][1].imag()*resid[1].imag()*F[1]*currWgt[1];
        res += AllDer[iaux4+1][3].real()*resid[3].real()*F[3]*currWgt[3];
        res += AllDer[iaux4+1][3].imag()*resid[3].imag()*F[3]*currWgt[3];
        DerVec[iaux4+1] += res;
      };

      // Add up to the Chi Square:
      for(iaux4=0;iaux4<4;iaux4++){
        ChiSq += resid[iaux4].real()*resid[iaux4].real()*currWgt[iaux4]*F[iaux4];
        ChiSq += resid[iaux4].imag()*resid[iaux4].imag()*currWgt[iaux4]*F[iaux4];
      };

      
// Add up the product of derivatives to the Hessian matrix:
 //   if(DEBUG){
 //   for (k=0;k<4;k++){  printf(" POL %i:\n",k);
 //     for (l=0;l<NPar;l++){ printf(" (% .3e, % .3ej) ",AllDer[l][k].real());  };
 //   };
 //  };

    if(doDeriv){
     for(k=0;k<currPar;k++){
       for(l=0;l<=k;l++){
         for(iaux4=0;iaux4<4;iaux4++){
           res  = AllDer[VisPar[k]][iaux4].real()*AllDer[VisPar[l]][iaux4].real()*F[iaux4]*currWgt[iaux4];
           res += AllDer[VisPar[k]][iaux4].imag()*AllDer[VisPar[l]][iaux4].imag()*F[iaux4]*currWgt[iaux4];
           Hessian[VisPar[k]*NPar + VisPar[l]] += res;

           if(DEBUG && VisPar[k]>39 && VisPar[l]>39 ){
      printf(" POL %i; DER %i / %i: (% .3e, % .3ej) | DER %i / %i: (% .3e, % .3ej) | H: % .3e \n",iaux4,k,VisPar[k],AllDer[VisPar[k]][iaux4].real(),AllDer[VisPar[k]][iaux4].imag(),l,VisPar[l],AllDer[VisPar[l]][iaux4].real(),AllDer[VisPar[l]][iaux4].imag(),Hessian[VisPar[k]*NPar + VisPar[l]]);
           };

         };
         if(k!=l){Hessian[VisPar[l]*NPar + VisPar[k]] = Hessian[VisPar[k]*NPar + VisPar[l]];};
       };
     };
    }; 

    if (doResid){ // Residuals for each correlation product:
      DATA[s][iaux2  ] = resid[0];
      DATA[s][iaux2+1] = resid[1]; 
      DATA[s][iaux2+2] = resid[2];
      DATA[s][iaux2+3] = resid[3];
    };

     };  // Comes from if(Wgt...)

    }; // Comes from loop over channels.
  };  // Comes from loop over visibilities

  }; // Comes from loop over calibrator sources.

// Reduced Chi Square:
  ChiSq /= (double)(Ndata);

// Ratio of Stokes I between model and data: 
  for(s=0; s<NCalSour; s++){  
    IRatio[s] /= Iwt[s];  
// Update ratio (only if it is worth, i.e. changes above 0.1%):  
    if (std::abs(IRatio[s] - 1.0) > 1.e-3){TotFlux[s] *= IRatio[s];}; 
  };
  
/*
  if (DEBUG){
    for (i=0;i<NPar;i++){
      for (j=0; j<NPar;j++){
        printf(" % .4e",Hessian[j*NPar+i]);
      };
      printf("\n");
    };
  };
*/

// Return the reduced ChiSquare (and the extreme I ratios):
  Err = Py_BuildValue("[d,d]",ChiSq,TotFlux);

  NITER += 1;
  printf("\r ITER %i. ChiSq %.3e ; Model Factors: ",NITER,ChiSq);
  for (s=0; s<NCalSour;s++){printf(" %i -> %.3e ",s,IRatio[s]);};
  fflush(stdout);
  return Err;

};