← Back to branch summary

~i-martividal/casa-poltools/trunk

« back to all changes in this revision

Viewing changes to PolSolver.cpp

  • Committer: I. Marti-Vidal
  • Date: 2019-11-19 14:43:44 UTC
  • Revision ID: i.marti-vidal@uv.es-20191119144344-z4wzo47zlyvqtgdh
DAR version

Show diffs side-by-side

added added

removed removed

Lines of Context:
PolSolver.cpp
59
59
int **A1, **A2;
60
60
double *VARS, *Hessian, *DerVec;
61
61
double ***Wgt; // = new double**[1];
 
62
double **WgtCorr;
62
63
 
63
64
double *TotFlux;
64
65
cplx64d **DATA, **COMPS;
91
92
{
92
93
 
93
94
  PyObject *DATAPy, *A1Py, *A2Py, *COMPSPy, *EPAsPy, *EMAsPy, *PSouPy, *VARSPy, *FluxPy; 
94
 
  PyObject *PAntRPy, *PAntLPy, *WgtPy, *VSouPy, *VAntRPy, *VAntLPy, *HessianPy, *ResVecPy;
 
95
  PyObject *PAntRPy, *PAntLPy, *WgtPy, *WgtCorrPy, *VSouPy, *VAntRPy, *VAntLPy, *HessianPy, *ResVecPy;
95
96
 
96
97
  PyObject *Err;
97
98
  int i,j;
98
99
 
99
100
 
100
 
  if (!PyArg_ParseTuple(args,"OOOOOOOOOOOOOOOOObb", &DATAPy, &WgtPy, &A1Py, &A2Py, &COMPSPy, &EPAsPy, &EMAsPy, &PSouPy, &PAntRPy, &PAntLPy, &VSouPy, &VAntRPy, &VAntLPy, &VARSPy, &HessianPy, &ResVecPy,&FluxPy, &doCirc, &doOrd1))
 
101
  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))
101
102
    {printf("Failed setData! Wrong arguments!\n"); fflush(stdout); Err = Py_BuildValue("i",-1); return Err;};
102
103
 
103
104
//  delete NSou, Nvis; 
112
113
  A1 = new int*[NCalSour];
113
114
  A2 = new int*[NCalSour];
114
115
  Wgt = new double**[NCalSour];
 
116
  WgtCorr = new double *[NCalSour];
115
117
 
116
118
  Hessian = (double *)PyArray_DATA(HessianPy);
117
119
  DerVec = (double *)PyArray_DATA(ResVecPy);
143
145
 
144
146
    if(DEBUG){printf("   %i visibilities; %i subcomponents. \n",Nvis[i],NSou[i]); fflush(stdout);};
145
147
 
 
148
 
146
149
    DATA[i] = (cplx64d *)PyArray_DATA(PyList_GetItem(DATAPy,i));
147
150
    A1[i] = (int *)PyArray_DATA(PyList_GetItem(A1Py,i));
148
151
    A2[i] = (int *)PyArray_DATA(PyList_GetItem(A2Py,i));
159
162
 
160
163
    };
161
164
 
 
165
    WgtCorr[i] = (double *)PyArray_DATA(PyList_GetItem(WgtCorrPy,i));
 
166
 
162
167
  };
163
168
 
164
169
  if(DEBUG){printf("Getting pointers to parameter indices\n"); fflush(stdout);};
253
258
 
254
259
  double Ifac, ChiSq, res, F[4];
255
260
  double *IRatio, *Iwt;
256
 
  bool doDeriv;
 
261
  bool doDeriv, doResid, doWgt;
257
262
  int i,j,k,l, s;
258
263
  long Ndata;
259
264
 
268
273
  }; // Relative weights for RR, RL, LR, LL / XX, XY, YX, YY
269
274
  
270
275
  
271
 
  double currWgt[4], ord2;
 
276
  double currWgt[4], ord2, Pf2Abs, PfAbs, ItotAbs;
272
277
  
273
278
  
274
279
  // Turn on/off 2nd order corrections:
281
286
  int iaux1, iaux2, iaux3, iaux4;
282
287
  int printedGood = 0;
283
288
 
284
 
  if (!PyArg_ParseTuple(args,"b", &doDeriv))
 
289
  if (!PyArg_ParseTuple(args,"bbb", &doDeriv, &doResid, &doWgt))
285
290
    {printf("Failed getHessian! Wrong arguments!\n"); 
286
291
     Err = Py_BuildValue("i",-1); fflush(stdout); return Err;};
287
292
 
303
308
  cplx64d RRc, RLc, LRc, LLc;
304
309
  cplx64d Im = cplx64d(0.,1.);
305
310
 
 
311
  cplx64d DTEff;
 
312
 
306
313
  cplx64d RR, RL, LR, LL, auxC;
307
314
  cplx64d resid[4];
308
315
  cplx64d AllDer[NPar][4];
441
448
 
442
449
// Apply leakage to model visibilities:
443
450
      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]]);
 
451
 
444
452
      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]]);
 
453
      PfAbs = std::abs(RLc);
 
454
 
445
455
      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]]);
 
456
      Pf2Abs = std::abs(LRc);
 
457
 
 
458
      ItotAbs = std::abs(Itot);
 
459
 
 
460
      if (doWgt && j==0 && PfAbs>0.0 && Pf2Abs > 0.0 && ItotAbs > 0.0){
 
461
        WgtCorr[s][i] = std::max(PfAbs/ItotAbs,Pf2Abs/ItotAbs);
 
462
      };
 
463
 
446
464
      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]]);
447
465
 
448
 
 
449
466
      if(DEBUG && NITER<2 && printedGood<10 && j==Nchan/2){
450
467
        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());
451
468
        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());
617
634
         if(k!=l){Hessian[VisPar[l]*NPar + VisPar[k]] = Hessian[VisPar[k]*NPar + VisPar[l]];};
618
635
       };
619
636
     };
620
 
    } else { // Residuals for each correlation product:
 
637
    }; 
 
638
 
 
639
    if (doResid){ // Residuals for each correlation product:
621
640
      DATA[s][iaux2  ] = resid[0];
622
641
      DATA[s][iaux2+1] = resid[1]; 
623
642
      DATA[s][iaux2+2] = resid[2];
624
643
      DATA[s][iaux2+3] = resid[3];
625
 
   };
 
644
    };
626
645
 
627
646
     };  // Comes from if(Wgt...)
628
647