Index: trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c =================================================================== --- trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c (revision 39584) +++ trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c (revision 39585) @@ -45,5 +45,12 @@ remove_argument (N, &argc, argv); } - + + + double binning_step = 0.0; + if ((N = get_argument (argc, argv, "-binning"))) { + remove_argument (N, &argc, argv); + binning_step = atof(argv[N]); + remove_argument (N, &argc, argv); + } if (argc != 6) { @@ -51,4 +58,5 @@ gprint (GP_ERR, " -max-iterations : maximum number of IRLS iterations to run (default 10)\n"); gprint (GP_ERR, " -outlier-limit : fraction of average weight to reject on (default 0.1)\n"); + gprint (GP_ERR, " -binning : fraction of a year to use in a bin for initial pass (default 0.0/no binning)\n"); return (FALSE); } @@ -111,5 +119,5 @@ fitdata.t, fitdata.pX, fitdata.pY, fitdata.Wx, fitdata.Wy, - fitdata.Npts, max_iterations, outlier_limit, VERBOSE)) { + fitdata.Npts, max_iterations, outlier_limit, binning_step, VERBOSE)) { return FALSE; } @@ -242,5 +250,5 @@ 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 Npts, int max_iterations, double outlier_limit, double binning_step, int VERBOSE) { int i,j; @@ -294,8 +302,41 @@ // Solve OLS equation - if (!weighted_LS_PLX(T,pR,pD,X,Wx,Y,Wy,Npts,A,B,VERBOSE)) { - // Handle fail cases gracefully. - return(FALSE); - } + + // Solve OLS equation + if (binning_step == 0.0) { + if (!weighted_LS_PLX(T,pR,pD,X,Wx,Y,Wy,Npts, + A,B,VERBOSE)) { + // Handle fail case + return(FALSE); + } + } + else { + int Nbins; + double *Tbin = NULL; + double *pRbin = NULL; + double *pDbin = NULL; + double *Xbin = NULL; + double *Ybin = NULL; + double *WXbin = NULL; + double *WYbin = NULL; + + if (!bin_points_PLX(T, pR, pD, X, Wx, Y, Wy, Npts, + &Tbin, &pRbin, &pDbin, &Xbin, &WXbin, &Ybin, &WYbin, &Nbins, binning_step) ) { + return(FALSE); + } + if (!weighted_LS_PLX(Tbin,pRbin,pDbin,Xbin,WXbin,Ybin,WYbin, Nbins, + A,B,VERBOSE)) { + return(FALSE); + } + + FREE(Tbin); + FREE(pRbin); + FREE(pDbin); + FREE(Xbin); + FREE(Ybin); + FREE(WXbin); + FREE(WYbin); + } + // Calculate r vector of residuals and least squares sigma @@ -576,2 +617,183 @@ return TRUE; } + +int bin_points_PLX (double *T, double *pR, double *pD, double *X, double *WX, double *Y, double *WY, int Npts, + double **Tbin, double **pRbin, double **pDbin, double **Xbin, double **WXbin, double **Ybin, double **WYbin, int *Nbins, double binning_step) { + double *T_tmp = NULL; + int i,j,k,l; + + int Nbins_test; + int *bin_test = NULL; + + double *T_int = NULL, *pD_int = NULL, *pR_int = NULL, *X_int = NULL, *Y_int = NULL, *WX_int = NULL, *WY_int = NULL; + double T_min, T_max; + + // Allocate more bins than we need + ALLOCATE(T_tmp,double,Npts); + for (i = 0; i < Npts; i++) { + T_tmp[i] = T[i]; + } + dsort(T_tmp,Npts); + + Nbins_test = floor((T_tmp[Npts - 1] - T_tmp[0]) / binning_step) + 1; + ALLOCATE(bin_test,int,Nbins_test); + + T_min = T_tmp[0]; + T_max = T_tmp[Npts - 1]; + + // Check how many bins are filled. + for (j = 0; j < Nbins_test; j++) { + bin_test[j] = 0; + } + for (i = 0; i < Npts; i++) { + j = floor((T_tmp[i] - T_tmp[0]) / binning_step); + bin_test[j] ++; + } + + *Nbins = 0; + for (j = 0; j < Nbins_test; j++) { + if (bin_test[j] != 0) { + *Nbins = *Nbins + 1; + } + } + + double *Tbin_in = NULL; + double *pRbin_in = NULL; + double *pDbin_in = NULL; + double *Xbin_in = NULL; + double *Ybin_in = NULL; + double *WXbin_in = NULL; + double *WYbin_in = NULL; + + ALLOCATE(Tbin_in,double,*Nbins); + ALLOCATE(pRbin_in, double, *Nbins); + ALLOCATE(pDbin_in, double, *Nbins); + ALLOCATE(Xbin_in,double,*Nbins); + ALLOCATE(WXbin_in,double,*Nbins); + ALLOCATE(Ybin_in,double,*Nbins); + ALLOCATE(WYbin_in,double,*Nbins); + + k = 0; + for (j = 0; j < Nbins_test; j++) { + if (bin_test[j] == 0) { // No data for this bin + continue; + } + else { + // Allocate internal arrays to hold the data that goes into this bin. + ALLOCATE(T_int,double, bin_test[j]); + ALLOCATE(pR_int,double, bin_test[j]); + ALLOCATE(pD_int,double, bin_test[j]); + ALLOCATE(X_int,double, bin_test[j]); + ALLOCATE(Y_int,double, bin_test[j]); + ALLOCATE(WX_int,double, bin_test[j]); + ALLOCATE(WY_int,double, bin_test[j]); + + // Fill those arrays. + l = 0; + for (i = 0; i < Npts; i++) { + if ((T[i] >= T_min + j * binning_step)&& + (T[i] < T_min + (j + 1) * binning_step)) { + T_int[l] = T[i]; + pR_int[l] = pR[i]; + pD_int[l] = pD[i]; + X_int[l] = X[i]; + Y_int[l] = Y[i]; + WX_int[l] = WX[i]; + WY_int[l] = WY[i]; + l++; + // printf("bin: %d Points: %d Point: %d %d %f %f %f : %f %f\n",j,bin_test[j],l,i,T[i],X[i],Y[i], + // T_min + j * binning_step,T_min + (j+1) * binning_step); + } + } // End loop over input points + + // Make a decision what the binned values should be. + if (l == 1) { + Tbin_in[k] = T_int[0]; + Xbin_in[k] = X_int[0]; + Ybin_in[k] = Y_int[0]; + pRbin_in[k] = pR_int[0]; + pDbin_in[k] = pD_int[0]; + WXbin_in[k] = WX_int[0]; + WYbin_in[k] = WY_int[0]; + } + else { + // I think I'm going with medians. + dsort(T_int,l); + dsort(pR_int,l); + dsort(pD_int,l); + dsort(X_int,l); + dsort(Y_int,l); + + // The median gives the midpoint of all the measurements + if ((l % 2) == 0) { + Tbin_in[k] = 0.5*(T_int[(int)(0.5*l)] + T_int[(int)(0.5*l) - 1]); + pRbin_in[k] = 0.5*(pR_int[(int)(0.5*l)] + pR_int[(int)(0.5*l) - 1]); + pDbin_in[k] = 0.5*(pD_int[(int)(0.5*l)] + pD_int[(int)(0.5*l) - 1]); + Xbin_in[k] = 0.5*(X_int[(int)(0.5*l)] + X_int[(int)(0.5*l) - 1]); + Ybin_in[k] = 0.5*(Y_int[(int)(0.5*l)] + Y_int[(int)(0.5*l) - 1]); + } else { + Tbin_in[k] = T_int[(int)(0.5*l)]; + pDbin_in[k] = pD_int[(int)(0.5*l)]; + pRbin_in[k] = pR_int[(int)(0.5*l)]; + Xbin_in[k] = X_int[(int)(0.5*l)]; + Ybin_in[k] = Y_int[(int)(0.5*l)]; + } + + // The scatter between points is probably the most useful measurement for the error. + for (i = 0; i < l; i++) { + X_int[i] = fabs(X_int[i] - Xbin_in[k]); + Y_int[i] = fabs(Y_int[i] - Ybin_in[k]); + } + dsort(X_int,l); + dsort(Y_int,l); + + if ((l % 2) == 0) { + WXbin_in[k] = 1.0 / (0.5*(X_int[(int)(0.5*l)] + X_int[(int)(0.5*l) - 1])); + WYbin_in[k] = 1.0 / (0.5*(Y_int[(int)(0.5*l)] + Y_int[(int)(0.5*l) - 1])); + } else { + WXbin_in[k] = 1.0 / (X_int[(int)(0.5*l)]); + WYbin_in[k] = 1.0 / (Y_int[(int)(0.5*l)]); + } + } + + if (WXbin_in[k] == 0.0) { + WXbin_in[k] = WX_int[0]; + } + if (WYbin_in[k] == 0.0) { + WYbin_in[k] = WY_int[0]; + } + if (!isfinite(WXbin_in[k])) { + WXbin_in[k] = 1.0; + } + if (!isfinite(WYbin_in[k])) { + WYbin_in[k] = 1.0; + } + + // Increment bin index. + k++; + + // Clean up + FREE(T_int); + FREE(pR_int); + FREE(pD_int); + FREE(X_int); + FREE(Y_int); + FREE(WX_int); + FREE(WY_int); + } + } // End loop over initial bin set. + + *Tbin = Tbin_in; + *pRbin = pRbin_in; + *pDbin = pDbin_in; + *Xbin = Xbin_in; + *Ybin = Ybin_in; + *WXbin = WXbin_in; + *WYbin = WYbin_in; + + FREE(T_tmp); + FREE(bin_test); + + return TRUE; +} + Index: trunk/Ohana/src/opihi/cmd.astro/fitpm_irls.c =================================================================== --- trunk/Ohana/src/opihi/cmd.astro/fitpm_irls.c (revision 39584) +++ trunk/Ohana/src/opihi/cmd.astro/fitpm_irls.c (revision 39585) @@ -25,6 +25,14 @@ } + double binning_step = 0.0; + if ((N = get_argument (argc, argv, "-binning"))) { + remove_argument (N, &argc, argv); + binning_step = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + if (argc != 6) { - gprint (GP_ERR, "USAGE: fitplx_irls (ra) (dR) (dec) (dD) (mjd) [-mask mask]\n"); + gprint (GP_ERR, "USAGE: fitplx_irls (ra) (dR) (dec) (dD) (mjd) [-mask mask] [-binning step_size]\n"); // what about the errors? return (FALSE); @@ -105,5 +113,5 @@ PlxFit fit; - if (!FitPMonly_IRLS (&fit, X, dX, Y, dY, t, n, VERBOSE)) { + if (!FitPMonly_IRLS (&fit, X, dX, Y, dY, t, n, binning_step, VERBOSE)) { return FALSE; } @@ -140,5 +148,5 @@ /* do we want an init function which does the alloc and a clear function to free? */ -int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) { +int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, double binning_step, int VERBOSE) { int i,j; @@ -196,9 +204,33 @@ } - // Solve OLS equation - if (!weighted_LS_PM(T,X,Wx,Y,Wy,Npts, - A,B,VERBOSE)) { - // Handle fail case - return(FALSE); + // Solve OLS equation + if (binning_step == 0.0) { + if (!weighted_LS_PM(T,X,Wx,Y,Wy,Npts, + A,B,VERBOSE)) { + // Handle fail case + return(FALSE); + } + } + else { + int Nbins; + double *Tbin = NULL; + double *Xbin = NULL; + double *Ybin = NULL; + double *WXbin = NULL; + double *WYbin = NULL; + + if (!bin_points(T, X, Wx, Y, Wy, Npts, + &Tbin, &Xbin, &WXbin, &Ybin, &WYbin, &Nbins, binning_step) ) { + return(FALSE); + } + if (!weighted_LS_PM(Tbin,Xbin,WXbin,Ybin,WYbin, Nbins, + A,B,VERBOSE)) { + return(FALSE); + } + FREE(Tbin); + FREE(Xbin); + FREE(Ybin); + FREE(WXbin); + FREE(WYbin); } @@ -424,4 +456,173 @@ } +int bin_points (double *T, double *X, double *WX, double *Y, double *WY, int Npts, + double **Tbin, double **Xbin, double **WXbin, double **Ybin, double **WYbin, int *Nbins, double binning_step) { + double *T_tmp = NULL; + int i,j,k,l; + + int Nbins_test; + int *bin_test = NULL; + + double *T_int = NULL, *X_int = NULL, *Y_int = NULL, *WX_int = NULL, *WY_int = NULL; + double T_min, T_max; + + + // printf("In binning code\n"); + + // Allocate more bins than we need + ALLOCATE(T_tmp,double,Npts); + for (i = 0; i < Npts; i++) { + T_tmp[i] = T[i]; + } + dsort(T_tmp,Npts); + + Nbins_test = floor((T_tmp[Npts - 1] - T_tmp[0]) / binning_step) + 1; + ALLOCATE(bin_test,int,Nbins_test); + + T_min = T_tmp[0]; + T_max = T_tmp[Npts - 1]; + + // printf("Found %d points, putting into %d bins between %f and %f\n",Npts,Nbins_test,T_min,T_max); + // Check how many bins are filled. + for (j = 0; j < Nbins_test; j++) { + bin_test[j] = 0; + } + for (i = 0; i < Npts; i++) { + j = floor((T_tmp[i] - T_tmp[0]) / binning_step); + bin_test[j] ++; + } + + *Nbins = 0; + for (j = 0; j < Nbins_test; j++) { + if (bin_test[j] != 0) { + *Nbins = *Nbins + 1; + } + } + + // printf("Using %d actual bins\n",*Nbins); + + double *Tbin_in = NULL; + double *Xbin_in = NULL; + double *Ybin_in = NULL; + double *WXbin_in = NULL; + double *WYbin_in = NULL; + + ALLOCATE(Tbin_in,double,*Nbins); + ALLOCATE(Xbin_in,double,*Nbins); + ALLOCATE(WXbin_in,double,*Nbins); + ALLOCATE(Ybin_in,double,*Nbins); + ALLOCATE(WYbin_in,double,*Nbins); + + k = 0; + for (j = 0; j < Nbins_test; j++) { + if (bin_test[j] == 0) { // No data for this bin + continue; + } + else { + // printf("%d bin %d, N = %d\n",k,j,bin_test[j]); + // Allocate internal arrays to hold the data that goes into this bin. + ALLOCATE(T_int,double, bin_test[j]); + ALLOCATE(X_int,double, bin_test[j]); + ALLOCATE(Y_int,double, bin_test[j]); + ALLOCATE(WX_int,double, bin_test[j]); + ALLOCATE(WY_int,double, bin_test[j]); + + // Fill those arrays. + l = 0; + for (i = 0; i < Npts; i++) { + if ((T[i] >= T_min + j * binning_step)&& + (T[i] < T_min + (j + 1) * binning_step)) { + T_int[l] = T[i]; + X_int[l] = X[i]; + Y_int[l] = Y[i]; + WX_int[l] = WX[i]; + WY_int[l] = WY[i]; + // printf("%d %d %f %f %f %f %f\n",i,l,T[i],X[i],Y[i],WX[i],WY[i]); + l++; + } + } // End loop over input points + // printf("%d bin %d, N = %d, l = %d\n",k,j,bin_test[j],l); + + // Make a decision what the binned values should be. + if (l == 1) { + Tbin_in[k] = T_int[0]; + Xbin_in[k] = X_int[0]; + Ybin_in[k] = Y_int[0]; + WXbin_in[k] = WX_int[0]; + WYbin_in[k] = WY_int[0]; + } + else { + // I think I'm going with medians. + dsort(T_int,l); + dsort(X_int,l); + dsort(Y_int,l); + + // The median gives the midpoint of all the measurements + if ((l % 2) == 0) { + Tbin_in[k] = 0.5*(T_int[(int)(0.5*l)] + T_int[(int)(0.5*l) - 1]); + Xbin_in[k] = 0.5*(X_int[(int)(0.5*l)] + X_int[(int)(0.5*l) - 1]); + Ybin_in[k] = 0.5*(Y_int[(int)(0.5*l)] + Y_int[(int)(0.5*l) - 1]); + } else { + Tbin_in[k] = T_int[(int)(0.5*l)]; + Xbin_in[k] = X_int[(int)(0.5*l)]; + Ybin_in[k] = Y_int[(int)(0.5*l)]; + } + + // The scatter between points is probably the most useful measurement for the error. + for (i = 0; i < l; i++) { + X_int[i] = fabs(X_int[i] - Xbin_in[k]); + Y_int[i] = fabs(Y_int[i] - Ybin_in[k]); + } + dsort(X_int,l); + dsort(Y_int,l); + + if ((l % 2) == 0) { + WXbin_in[k] = 1.0 / (0.5*(X_int[(int)(0.5*l)] + X_int[(int)(0.5*l) - 1])); + WYbin_in[k] = 1.0 / (0.5*(Y_int[(int)(0.5*l)] + Y_int[(int)(0.5*l) - 1])); + } else { + WXbin_in[k] = 1.0 / (X_int[(int)(0.5*l)]); + WYbin_in[k] = 1.0 / (Y_int[(int)(0.5*l)]); + } + } + + if (WXbin_in[k] == 0.0) { + WXbin_in[k] = WX_int[0]; + } + if (WYbin_in[k] == 0.0) { + WYbin_in[k] = WY_int[0]; + } + if (!isfinite(WXbin_in[k])) { + WXbin_in[k] = 1.0; + } + if (!isfinite(WYbin_in[k])) { + WYbin_in[k] = 1.0; + } + + // printf("%f %f %f %f %f\n",Tbin_in[k],Xbin_in[k],Ybin_in[k],WXbin_in[k],WYbin_in[k]); + // Increment bin index. + k++; + + // Clean up + FREE(T_int); + FREE(X_int); + FREE(Y_int); + FREE(WX_int); + FREE(WY_int); + } + } // End loop over initial bin set. + + *Tbin = Tbin_in; + *Xbin = Xbin_in; + *Ybin = Ybin_in; + *WXbin = WXbin_in; + *WYbin = WYbin_in; + + FREE(T_tmp); + FREE(bin_test); + + return TRUE; +} + + double weight_cauchy (double x) { double r = x / 2.385; @@ -450,5 +651,5 @@ dsort(x,N); - if (N % 2) { + if ((N % 2) == 0) { median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]); } else { @@ -462,5 +663,5 @@ dsort(x,N); - if (N % 2) { + if ((N % 2) == 0) { median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]); } else { @@ -468,4 +669,5 @@ } + FREE(x); return(median); } Index: trunk/Ohana/src/opihi/include/astro.h =================================================================== --- trunk/Ohana/src/opihi/include/astro.h (revision 39584) +++ trunk/Ohana/src/opihi/include/astro.h (revision 39585) @@ -58,9 +58,14 @@ /***** */ -int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, 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 FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, double binning_step, 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, double binning_step, int VERBOSE); int weighted_LS_PLX (double *T, double *pR, double *pD, double *X, double *WX, double *Y, double *WY, int Npts, double **A, double **B, int VERBOSE); int weighted_LS_PM (double *T, double *X, double *WX, double *Y, double *WY, int Npts, double **A, double **B, int VERBOSE); +int bin_points (double *T, double *X, double *WX, double *Y, double *WY, int Npts, + double **Tbin, double **Xbin, double **WXbin, double **Ybin, double **WYbin, int *Nbins, double binning_step); +int bin_points_PLX (double *T, double *pR, double *pD, double *X, double *WX, double *Y, double *WY, int Npts, + double **Tbin, double **pRbin, double **pDbin, double **Xbin, double **WXbin, double **Ybin, double **WYbin, int *Nbins, double binning_step); + double weight_cauchy (double x);