Index: trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c =================================================================== --- trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c (revision 39180) +++ trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c (revision 39181) @@ -15,5 +15,7 @@ double *dY; - double *IRLS_u; + double *Wx; + double *Wy; + int *index; int Npts; @@ -32,5 +34,8 @@ int Nfit; int getChisq; -} PlxFit_IRLS; +} PlxFit; + + +int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE); int VectorRobustStats_IRLS (Vector *vector, double *median, double *sigma); @@ -45,9 +50,9 @@ int PlxBootstrapResample_IRLS (PlxFitData *src, PlxFitData *tgt); -int FitPMandPar_IRLS (PlxFit_IRLS *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, double *IRLS_u, int Npts, int max_iterations, int VERBOSE); +int FitPMandPar_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, double *Wx, double *Wy, int Npts, int max_iterations, double outlier_limit, int VERBOSE); int sun_ecliptic_IRLS (double mjd, double *lambda, double *beta, double *epsilon, double *Radius); int ParFactor_IRLS (double *pR, double *pD, double RA, double DEC, double Time); -int IRLS_converged (PlxFit_IRLS *fit); +int IRLS_converged (PlxFit *fit); int Plx_weighted_LS (double *T, double *pR, double *pD, double *X, double *WX, double *Y, double *WY, int Npts, double **A, double **B, int VERBOSE); @@ -114,4 +119,19 @@ remove_argument (N, &argc, argv); } + + Vector *outMask = NULL; + if ((N = get_argument (argc, argv, "-out-mask"))) { + remove_argument (N, &argc, argv); + if (!(outMask = SelectVector(argv[N], ANYVECTOR, TRUE))) return FALSE; + remove_argument (N, &argc, argv); + } + + double outlier_limit = 0.1; + if ((N = get_argument (argc, argv, "-outlier-limit"))) { + remove_argument (N, &argc, argv); + outlier_limit = atof(argv[N]); + remove_argument (N, &argc, argv); + } + if (argc != 6) { @@ -119,5 +139,7 @@ gprint (GP_ERR, " -outlier-tests Nsamples dPsigMax : run Nsample bootstrap-resamples to define the path deviations and reject based on dPsigMax\n"); gprint (GP_ERR, " -dPsig vec : save path deviations in vec\n"); - gprint (GP_ERR, " -max_iterations : maximum number of IRLS iterations to run\n"); + gprint (GP_ERR, " -max-iterations : maximum number of IRLS iterations to run (default 10)\n"); + gprint (GP_ERR, " -outMask vect : save outlier mask in vec\n"); + gprint (GP_ERR, " -outlier-limit : fraction of average weight to reject on (default 0.1)\n"); return (FALSE); } @@ -151,5 +173,6 @@ int Ntotal = tvec->Nelements; // XXX check other lengths if (dPvec) ResetVector (dPvec, OPIHI_FLT, Ntotal); - + if (outMask) ResetVector (outMask, OPIHI_INT, Ntotal); + double Rmean, Dmean, Tmean; PlxSetMeanEpoch_IRLS (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal); @@ -166,5 +189,5 @@ PlxSetEpochPosition_IRLS (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean); - PlxFit_IRLS fit; memset (&fit, 0, sizeof(PlxFit_IRLS)); + PlxFit fit; memset (&fit, 0, sizeof(PlxFit)); if (0) { @@ -195,15 +218,41 @@ fitdata.Y, fitdata.dY, fitdata.t, fitdata.pX, fitdata.pY, - fitdata.IRLS_u, - fitdata.Npts, max_iterations, VERBOSE)) { + fitdata.Wx, fitdata.Wy, + fitdata.Npts, max_iterations, outlier_limit, VERBOSE)) { return FALSE; } - + // Set the output mask based on the input mask and the outlier limits. + if (outMask) { + double Sum_Wx = 0; + double Sum_Wy = 0; + int *outMask_V = outMask->elements.Int; + for (i = 0; i < fitdata.Npts; i++) { + Sum_Wx += fitdata.Wx[i]; + Sum_Wy += fitdata.Wy[i]; + } + for (i = 0; i < Ntotal; i++) { + outMask_V[i] = mask ? mask[i] : 0; + } + + for (i = 0; i < fitdata.Npts; i++) { + if ((fitdata.Wx[i] > outlier_limit * Sum_Wx / (1.0 * fitdata.Npts))|| + (fitdata.Wy[i] > outlier_limit * Sum_Wy / (1.0 * fitdata.Npts))) { + int n = fitdata.index[i]; + outMask_V[n] = 1; + + if (VERBOSE == 2) { + fprintf (stderr, "%f %f : %f %d : %f %f %f : %f %f %f %f\n", R[n], D[n], T[n], outMask_V[n], fitdata.t[i], fitdata.X[i], fitdata.Y[i], fitdata.Wx[i], fitdata.Wy[i], Sum_Wx, Sum_Wy); + } + } + } + } + + if (Nresample){ PlxFitData sample; PlxFitDataAlloc_IRLS (&sample, fitdata.Npts); - PlxFit_IRLS *testfit = NULL; - ALLOCATE (testfit, PlxFit_IRLS, Nresample); + PlxFit *testfit = NULL; + ALLOCATE (testfit, PlxFit, Nresample); int Ngood = 0; @@ -215,9 +264,9 @@ // fit the sample testfit[Ngood].getChisq = FALSE; - if (!FitPMandPar_IRLS (&testfit[Ngood], - sample.X, sample.dX, - sample.Y, sample.dY, sample.t, - sample.pX, sample.pY, - sample.IRLS_u, sample.Npts, 1, VERBOSE)) continue; + if (!FitPMandPar (&testfit[Ngood], + sample.X, sample.dX, + sample.Y, sample.dY, sample.t, + sample.pX, sample.pY, + sample.Npts, VERBOSE)) continue; Ngood ++; } @@ -292,11 +341,13 @@ set_variable ("Nfit", fit.Nfit); + + return (TRUE); } /* do we want an init function which does the alloc and a clear function to free? */ -int FitPMandPar_IRLS (PlxFit_IRLS *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, - double *IRLS_u, - int Npts, int max_iterations, int VERBOSE) { +int FitPMandPar_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, + double *Wx, double *Wy, + int Npts, int max_iterations, double outlier_limit, int VERBOSE) { int i,j; @@ -310,7 +361,7 @@ double sigma_ols, sigma_hat; // Sigma estimates. - double *Wx, *Wy; // Weight vectors. Not errors. + // double *Wx, *Wy; // Weight vectors. Not errors. double *rx, *ry; // Deviation from model - // double *u; // Deviation magnitude + double *u; // Deviation magnitude int converged; int iterations; @@ -338,8 +389,9 @@ ALLOCATE(Beta, double, 5); ALLOCATE(Beta_prev, double, 5); - ALLOCATE(Wx, double, Npts); - ALLOCATE(Wy, double, Npts); + // ALLOCATE(Wx, double, Npts); + // ALLOCATE(Wy, double, Npts); ALLOCATE(rx, double, Npts); ALLOCATE(ry, double, Npts); + ALLOCATE(u, double, Npts); // Convert the measurement errors into the initial weights. @@ -416,9 +468,9 @@ rx[i] = X[i] - (T[i] * B[1][0] + B[0][0] + B[4][0] * pR[i]); ry[i] = Y[i] - (T[i] * B[3][0] + B[2][0] + B[4][0] * pD[i]); - IRLS_u[i] = sqrt(SQ(rx[i] / dX[i]) + SQ(ry[i] / dY[i])); + u[i] = sqrt(SQ(rx[i] / dX[i]) + SQ(ry[i] / dY[i])); } // Calculate sigma_hat from distribution of residual magnitudes - sigma_hat = Plx_MAD(IRLS_u,Npts) / 0.6745; + sigma_hat = Plx_MAD(u,Npts) / 0.6745; // Check convergence @@ -511,6 +563,6 @@ double wx,wy; for (i = 0; i < Npts; i++) { - if ((Wx[i] > 0.1 * Sum_Wx / (1.0 * Npts))|| - (Wy[i] > 0.1 * Sum_Wy / (1.0 * Npts))) { + if ((Wx[i] > outlier_limit * Sum_Wx / (1.0 * Npts))|| + (Wy[i] > outlier_limit * Sum_Wy / (1.0 * Npts))) { Xf = fit[0].Ro + fit[0].uR*T[i] + fit[0].p*pR[i]; @@ -584,5 +636,6 @@ ALLOCATE (data->pX, double, N); ALLOCATE (data->pY, double, N); - ALLOCATE (data->IRLS_u, double, N); + ALLOCATE (data->Wx, double, N); + ALLOCATE (data->Wy, double, N); ALLOCATE (data->index, int, N); return TRUE; @@ -597,5 +650,6 @@ FREE (data->pX); FREE (data->pY); - FREE (data->IRLS_u); + FREE (data->Wx); + FREE (data->Wy); FREE (data->index); } @@ -614,5 +668,6 @@ tgt->pX[i] = src->pX[N]; tgt->pY[i] = src->pY[N]; - tgt->IRLS_u[i] = src->IRLS_u[N]; + tgt->Wx[i] = src->Wx[N]; + tgt->Wy[i] = src->Wy[N]; } return TRUE; @@ -674,5 +729,6 @@ pYmax = MAX (pYmax, fitdata->pY[Nsubset]); - fitdata->IRLS_u[Nsubset] = 0.0; + fitdata->Wx[Nsubset] = 1.0; + fitdata->Wy[Nsubset] = 1.0; fitdata->index[Nsubset] = i; Nsubset++; @@ -706,5 +762,5 @@ int i, n; - PlxFit_IRLS testfit; + PlxFit testfit; testfit.getChisq = FALSE; @@ -728,8 +784,8 @@ // fit the sample - if (!FitPMandPar_IRLS (&testfit, - sample.X, sample.dX, - sample.Y, sample.dY, sample.t, - sample.pX, sample.pY, sample.IRLS_u, sample.Npts, 1, VERBOSE)) continue; + if (!FitPMandPar (&testfit, + sample.X, sample.dX, + sample.Y, sample.dY, sample.t, + sample.pX, sample.pY, sample.Npts, VERBOSE)) continue; // fprintf (stderr, "%f +/- %f | %f %f\n", testfit.p, testfit.dp, testfit.uR, testfit.uD);