~i-martividal/uvmultifit/trunk-1 : revision 30

65

#include <sys/types.h>

66

#include <new>

67

#include <complex>

68

#include <iostream>

69

#include <fstream>

70

//#include <sstream>

71

//#include <stdlib.h>

72

//#include <string.h>

68

73

69

74

70

75

#if QUINN_FITTER == 0

135

140

static int Nspw, NCPU, nui, cIF, ncomp, npar, HankelOrder=1, maxnchan=0;

136

141

static int mode, Nants, NphClos, NampClos, Nbas;

137

142

static double cosDecRef, sinDecRef;

138

static bool compFixed, isModel, MixedG, isTime, doAmpClos, doPhClos, doClos, onlyClos;

143

static bool compFixed, isModel, MixedG, isTime, doAmpClos, doPhClos, doClos, onlyClos, doDump;

139

144

static double phClosWgt, ampClosWgt;

140

145

static int NparMax = 7; // Degrees of freedom for components (i.e., RA, Dec, Flux, Major, Ratio, PA, Inner).

141

146

398

403

399

404

int k,i,t,j,m,p, ii, currow[NparMax+HankelOrder];

400

405

401

double phase, uu, vv, ww, UVRad, Ampli, Ampli0, DerivRe[npar], DerivIm[npar];

406

double phase, uu, vv, ww, UVRad, Ampli, DerivRe[npar], DerivIm[npar];

402

407

double SPA, CPA, tA, tB, tempChi, deltat, Ellip, Ellip0;

403

Ellip0 = 1.; Ellip = 1.; Ampli0 = 1.;

408

Ellip0 = 1.; Ellip = 1.;

404

409

double wgtcorrA, tempRe[npar+1], tempIm[npar+1], tempAmp;

405

410

tempChi = 0.0;

406

411

441

446

};

442

447

443

448

449

450

FILE* phClosFile;

451

FILE* apClosFile;

452

453

if(doDump){

454

if(doPhClos){ phClosFile = fopen("UVMultiFit_Phase_Closures.dat", "wb");};

455

if(doAmpClos){apClosFile = fopen("UVMultiFit_Amplitude_Closures.dat", "wb");};

456

};

457

458

444

459

//bool write2mod = (mode == -3 || mode == -4 || mode == -5 || mode >= 0);

445

460

//bool writeDer = (mode==-1);

461

446

462

bool TriInv;

447

463

double tempD, tempR, tempI, cosphase, sinphase, cosphase0, sinphase0, PA;

448

464

cplx64 *totGain;

546

562

547

563

548

564

549

if (writeDer) {

565

// if (writeDer) {

550

566

551

567

552

568

553

569

// Closure Phases:

554

570

for(ii=0;ii<NphClos;ii++){

555

571

556

if(Observed[mod.phaseClosIdx[0][ii]] && Observed[mod.phaseClosIdx[1][ii]] && Observed[mod.phaseClosIdx[2][ii]]){

572

if(Observed[mod.phaseClosIdx[0][ii]] && Observed[mod.phaseClosIdx[1][ii]] && Observed[mod.phaseClosIdx[2][ii]] && mod.closBufferWgt[mod.phaseClosIdx[0][ii]]>0.0 && mod.closBufferWgt[mod.phaseClosIdx[1][ii]] && mod.closBufferWgt[mod.phaseClosIdx[2][ii]]>0.0){

557

573

558

574

// Error propagation:

559

575

tempD = 1./(mod.closBufferWgt[mod.phaseClosIdx[0][ii]]);

590

606

tempres0 = (TriSpec.real() - TriSpec0.real());

591

607

tempres1 = (TriSpec.imag() - TriSpec0.imag());

592

608

609

tempChi += (tempres0*tempres0 + tempres1*tempres1)*tempD;

610

611

if(doDump){

612

fwrite(&mod.BasIdx[0][mod.phaseClosIdx[0][ii]],sizeof(int),1,phClosFile);

613

fwrite(&mod.BasIdx[1][mod.phaseClosIdx[0][ii]],sizeof(int),1,phClosFile);

614

fwrite(&mod.BasIdx[1][mod.phaseClosIdx[1][ii]],sizeof(int),1,phClosFile);

615

fwrite(&deltat,sizeof(double),1,phClosFile);

616

tempR = std::arg(TriSpec);

617

fwrite(&tempR,sizeof(double),1,phClosFile);

618

tempR = std::arg(TriSpec0);

619

fwrite(&tempR,sizeof(double),1,phClosFile);

620

};

621

593

622

594

623

// Compute all gradients of the closure model w.r.t. fitting parameters:

624

if(writeDer){

595

625

for(p=1;p<npar+1;p++){

596

626

// New module of the model trispectrum:

597

627

tempR = mod.closBufferAbsMod[p][mod.phaseClosIdx[0][ii]]*mod.closBufferAbsMod[p][mod.phaseClosIdx[1][ii]]/mod.closBufferAbsMod[p][mod.phaseClosIdx[2][ii]];

607

637

mod.WorkGrad[widx][p-1] += tempres1*tempD*DerivIm[p-1];

608

638

};

609

639

610

tempChi += (tempres0*tempres0 + tempres1*tempres1)*tempD;

611

640

612

641

// Add to the Hessian:

613

642

for(p=0;p<npar;p++){

616

645

};

617

646

};

618

647

648

};

649

619

650

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

620

651

// FOR TESTING:

621

652

// tempD *= 1000.;

626

657

627

658

for(p=0;p<Nants;p++){

628

659

if(mod.triSpec[p]==1){

629

tri0 = 0; tri1 = 0; tri2 = 0;

660

tri0 = 0; tri1 = 0; tri2 = 0; other0 = -1; other1 = -1; other2=-1;

630

661

if(mod.BasIdx[0][mod.phaseClosIdx[0][ii]]==p){tri0=1;other0=mod.BasIdx[1][mod.phaseClosIdx[0][ii]];

631

662

} else if(mod.BasIdx[1][mod.phaseClosIdx[0][ii]]==p){tri0=2;other0=mod.BasIdx[0][mod.phaseClosIdx[0][ii]];};

632

663

if(mod.BasIdx[0][mod.phaseClosIdx[1][ii]]==p){tri1=1;other1=mod.BasIdx[1][mod.phaseClosIdx[1][ii]];

714

745

715

746

// printf("TriSpec: %i-%i / %i-%i / %i-%i: %.4e/%.4e ; %.4e/%.4e\n",mod.BasIdx[0][mod.phaseClosIdx[0][ii]],mod.BasIdx[1][mod.phaseClosIdx[0][ii]],mod.BasIdx[0][mod.phaseClosIdx[1][ii]],mod.BasIdx[1][mod.phaseClosIdx[1][ii]],mod.BasIdx[0][mod.phaseClosIdx[2][ii]],mod.BasIdx[1][mod.phaseClosIdx[2][ii]],TriSpec.real(),TriSpec0.real(), TriSpec.imag(),TriSpec0.imag());

716

747

748

tempChi += (tempres0*tempres0 + tempres1*tempres1)*tempD;

749

750

751

if(writeDer){

717

752

718

753

for(j=1;j<npar+1;j++){

719

754

TriSpec1 = BL1m[j]*BL2m[j]/BL3m[j];

724

759

mod.WorkGrad[widx][j-1] += tempres1*tempD*DerivIm[j-1];

725

760

};

726

761

727

tempChi += (tempres0*tempres0 + tempres1*tempres1)*tempD;

728

762

729

763

// Add to the Hessian:

730

764

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

731

for(m=j;m<npar;m++){

732

mod.WorkHess[widx][npar*j+m] += tempD*(DerivRe[j]*DerivRe[m]+DerivIm[j]*DerivIm[m]);

765

for(m=j;m<npar;m++){

766

mod.WorkHess[widx][npar*j+m] += tempD*(DerivRe[j]*DerivRe[m]+DerivIm[j]*DerivIm[m]);

767

};

733

768

};

734

};

769

770

};

735

771

736

772

};

737

773

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

740

776

};

741

777

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

742

778

743

744

745

779

};

746

780

};

747

781

748

782

749

783

750

751

752

753

754

755

756

757

758

784

// Closure Amplitudes:

759

785

for(ii=0;ii<NampClos;ii++){

760

786

if(Observed[mod.ampClosIdx[0][ii]] && Observed[mod.ampClosIdx[1][ii]] && Observed[mod.ampClosIdx[2][ii]] && Observed[mod.ampClosIdx[3][ii]]){

761

auxAmp1 = mod.closBufferAbs[mod.ampClosIdx[0][ii]];

787

// auxAmp1 = mod.closBufferAbs[mod.ampClosIdx[0][ii]];

762

788

tempD = 1./(mod.closBufferWgt[mod.ampClosIdx[0][ii]]);// *auxAmp1);//*auxAmp1);

763

auxAmp1 = mod.closBufferAbs[mod.ampClosIdx[1][ii]];

789

// auxAmp1 = mod.closBufferAbs[mod.ampClosIdx[1][ii]];

764

790

tempD += 1./(mod.closBufferWgt[mod.ampClosIdx[1][ii]]);// *auxAmp1);//*auxAmp1);

765

auxAmp1 = mod.closBufferAbs[mod.ampClosIdx[2][ii]];

791

// auxAmp1 = mod.closBufferAbs[mod.ampClosIdx[2][ii]];

766

792

tempD += 1./(mod.closBufferWgt[mod.ampClosIdx[2][ii]]);// *auxAmp1);//*auxAmp1);

767

auxAmp1 = mod.closBufferAbs[mod.ampClosIdx[3][ii]];

793

// auxAmp1 = mod.closBufferAbs[mod.ampClosIdx[3][ii]];

768

794

tempD += 1./(mod.closBufferWgt[mod.ampClosIdx[3][ii]]);// *auxAmp1);//*auxAmp1);

769

795

tempD = ampClosWgt/tempD;

770

796

auxAmp = mod.closBufferLog[mod.ampClosIdx[0][ii]]+mod.closBufferLog[mod.ampClosIdx[1][ii]]-(mod.closBufferLog[mod.ampClosIdx[2][ii]]+mod.closBufferLog[mod.ampClosIdx[3][ii]]);

772

798

773

799

774

800

tempres0 = (auxAmp-auxAmp0);

801

tempChi += tempres0*tempres0*tempD;

802

803

if(doDump){

804

fwrite(&mod.BasIdx[0][mod.ampClosIdx[0][ii]],sizeof(int),1,apClosFile);

805

fwrite(&mod.BasIdx[1][mod.ampClosIdx[0][ii]],sizeof(int),1,apClosFile);

806

fwrite(&mod.BasIdx[0][mod.ampClosIdx[1][ii]],sizeof(int),1,apClosFile);

807

fwrite(&mod.BasIdx[1][mod.ampClosIdx[1][ii]],sizeof(int),1,apClosFile);

808

fwrite(&deltat,sizeof(double),1,apClosFile);

809

fwrite(&auxAmp,sizeof(double),2,apClosFile);

810

fwrite(&auxAmp0,sizeof(double),2,apClosFile);

811

};

812

813

814

if(writeDer){

775

815

776

816

for(p=1;p<npar+1;p++){

777

817

auxAmp1 = mod.closBufferLogMod[p][mod.ampClosIdx[0][ii]]+mod.closBufferLogMod[p][mod.ampClosIdx[1][ii]]-(mod.closBufferLogMod[p][mod.ampClosIdx[2][ii]]+mod.closBufferLogMod[p][mod.ampClosIdx[3][ii]]);

779

819

mod.WorkGrad[widx][p-1] += tempres0*tempD*DerivRe[p-1];

780

820

};

781

821

782

tempChi += tempres0*tempres0*tempD;

783

822

784

823

785

824

for(p=0;p<npar;p++){

787

826

mod.WorkHess[widx][npar*p+m] += tempD*(DerivRe[p]*DerivRe[m]);

788

827

};

789

828

};

829

830

};

831

790

832

};

791

833

};

792

834

793

835

794

};

836

// };

795

837

796

838

797

839

1016

1058

} else {Ampli = 1.0;};

1017

1059

1018

1060

Ampli *= wgtcorrA;

1019

// if(p==0){Ampli0=Ampli;};

1020

1061

1021

// } else if (!PBlimit){Ampli = Ampli0;} else {Ampli = 0.0;};

1022

1062

} else {Ampli = 0.0;};

1023

1063

1024

1064

1142

1182

};

1143

1183

1144

1184

1185

if(doDump){

1186

if(doPhClos){fclose(phClosFile);};

1187

if(doAmpClos){fclose(apClosFile);};

1188

};

1189

1190

1191

1145

1192

1146

1193

if (Iam != -1){mod.Chi2[Iam] = tempChi; pthread_exit((void*) 0);}

1147

1194

else {mod.Chi2[0] = tempChi; return (void*) 0;};

1729

1776

int i,j;

1730

1777

int nparsq = npar*npar;

1731

1778

1732

int auxClos, auxClos2;

1779

int auxClos, auxClos2,Dump;

1733

1780

1734

1781

/* Parse the input tuple */

1735

if (!PyArg_ParseTuple(args, "iiiii",&cIF,&nui,&mode, &auxClos,&auxClos2)){printf("FAILED modelcomp!\n"); fflush(stdout); PyObject *ret = Py_BuildValue("i",-1); return ret;};

1782

if (!PyArg_ParseTuple(args, "iiiiii",&cIF,&nui,&mode, &auxClos,&auxClos2,&Dump)){printf("FAILED modelcomp!\n"); fflush(stdout); PyObject *ret = Py_BuildValue("i",-1); return ret;};

1736

1783

1737

1784

1738

1785

doClos = auxClos == 1;

1739

1786

onlyClos = auxClos2 == 1;

1787

doDump = Dump == 1;

1740

1788

1741

1789

// Zero the workers memory:

1742

1790

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

1856

1904

1857

1905

npy_intp Dims[1]; Dims[0] = Nants;

1858

1906

1859

printf("\nNants: %i DIMS: %i/%.3e/%.3e \n",Nants,Dims[0], Bins[0],Bins[1]);

1907

//printf("\nNants: %i DIMS: %i/%.3e/%.3e \n",Nants,Dims[0], Bins[0],Bins[1]);

1860

1908

1861

1909

PyRate = PyArray_SimpleNewFromData(1, Dims, NPY_FLOAT64, (void *) Rates);

1862

1910

PyDelay = PyArray_SimpleNewFromData(1, Dims, NPY_FLOAT64, (void *) Delays);