IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39585


Ignore:
Timestamp:
May 27, 2016, 5:17:08 PM (10 years ago)
Author:
watersc1
Message:

Added binned initial step method for irls pm/plx fitting.

Location:
trunk/Ohana/src/opihi
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c

    r39333 r39585  
    4545    remove_argument (N, &argc, argv);
    4646  }
    47        
     47
     48 
     49  double binning_step = 0.0;
     50  if ((N = get_argument (argc, argv, "-binning"))) {
     51    remove_argument (N, &argc, argv);
     52    binning_step = atof(argv[N]);
     53    remove_argument (N, &argc, argv);
     54  }
    4855 
    4956  if (argc != 6) {
    … …  
    5158    gprint (GP_ERR, "  -max-iterations : maximum number of IRLS iterations to run (default 10)\n");
    5259    gprint (GP_ERR, "  -outlier-limit : fraction of average weight to reject on (default 0.1)\n");
     60    gprint (GP_ERR, "  -binning <step_size> : fraction of a year to use in a bin for initial pass (default 0.0/no binning)\n");
    5361    return (FALSE);
    5462  }
    … …  
    111119                         fitdata.t, fitdata.pX, fitdata.pY,
    112120                         fitdata.Wx, fitdata.Wy,
    113                          fitdata.Npts, max_iterations, outlier_limit, VERBOSE)) {
     121                         fitdata.Npts, max_iterations, outlier_limit, binning_step, VERBOSE)) {
    114122    return FALSE;
    115123  }
    … …  
    242250int FitPMandPar_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD,
    243251                      double *Wx, double *Wy,
    244                       int Npts, int max_iterations, double outlier_limit, int VERBOSE) {
     252                      int Npts, int max_iterations, double outlier_limit, double binning_step, int VERBOSE) {
    245253
    246254  int i,j;
    … …  
    294302 
    295303  // Solve OLS equation
    296   if (!weighted_LS_PLX(T,pR,pD,X,Wx,Y,Wy,Npts,A,B,VERBOSE)) {
    297     // Handle fail cases gracefully.
    298     return(FALSE);
    299   }
     304 
     305  // Solve OLS equation
     306  if (binning_step == 0.0) {
     307    if (!weighted_LS_PLX(T,pR,pD,X,Wx,Y,Wy,Npts,
     308                         A,B,VERBOSE)) {
     309      // Handle fail case
     310      return(FALSE);
     311    }
     312  }
     313  else {
     314    int Nbins;
     315    double *Tbin = NULL;
     316    double *pRbin = NULL;
     317    double *pDbin = NULL;
     318    double *Xbin = NULL;
     319    double *Ybin = NULL;
     320    double *WXbin = NULL;
     321    double *WYbin = NULL;
     322
     323    if (!bin_points_PLX(T, pR, pD, X, Wx, Y, Wy, Npts,
     324                        &Tbin, &pRbin, &pDbin, &Xbin, &WXbin, &Ybin, &WYbin, &Nbins,  binning_step) ) {
     325      return(FALSE);
     326    }
     327    if (!weighted_LS_PLX(Tbin,pRbin,pDbin,Xbin,WXbin,Ybin,WYbin, Nbins,
     328                         A,B,VERBOSE)) {
     329      return(FALSE);
     330    }
     331
     332    FREE(Tbin);
     333    FREE(pRbin);
     334    FREE(pDbin);
     335    FREE(Xbin);
     336    FREE(Ybin);
     337    FREE(WXbin);
     338    FREE(WYbin);
     339  }
     340
    300341
    301342  // Calculate r vector of residuals and least squares sigma
    … …  
    576617  return TRUE;
    577618}
     619
     620int bin_points_PLX (double *T, double *pR, double *pD, double *X, double *WX, double *Y, double *WY, int Npts,
     621                    double **Tbin, double **pRbin, double **pDbin, double **Xbin, double **WXbin, double **Ybin, double **WYbin, int *Nbins, double binning_step) {
     622  double *T_tmp = NULL;
     623  int i,j,k,l;
     624
     625  int Nbins_test;
     626  int *bin_test = NULL;
     627
     628  double *T_int = NULL, *pD_int = NULL, *pR_int = NULL, *X_int = NULL, *Y_int = NULL, *WX_int = NULL, *WY_int = NULL;
     629  double T_min, T_max;
     630 
     631  // Allocate more bins than we need
     632  ALLOCATE(T_tmp,double,Npts);
     633  for (i = 0; i < Npts; i++) {
     634    T_tmp[i] = T[i];
     635  }
     636  dsort(T_tmp,Npts);
     637
     638  Nbins_test = floor((T_tmp[Npts - 1] - T_tmp[0]) / binning_step) + 1;
     639  ALLOCATE(bin_test,int,Nbins_test);
     640
     641  T_min = T_tmp[0];
     642  T_max = T_tmp[Npts - 1];
     643 
     644  // Check how many bins are filled.
     645  for (j = 0; j < Nbins_test; j++) {
     646    bin_test[j] = 0;
     647  }
     648  for (i = 0; i < Npts; i++) {
     649    j = floor((T_tmp[i] - T_tmp[0]) / binning_step);
     650    bin_test[j] ++;
     651  }
     652
     653  *Nbins = 0;
     654  for (j = 0; j < Nbins_test; j++) {
     655    if (bin_test[j] != 0) {
     656      *Nbins = *Nbins + 1;
     657    }
     658  }
     659
     660  double *Tbin_in = NULL;
     661  double *pRbin_in = NULL;
     662  double *pDbin_in = NULL;
     663  double *Xbin_in = NULL;
     664  double *Ybin_in = NULL;
     665  double *WXbin_in = NULL;
     666  double *WYbin_in = NULL;
     667 
     668  ALLOCATE(Tbin_in,double,*Nbins);
     669  ALLOCATE(pRbin_in, double, *Nbins);
     670  ALLOCATE(pDbin_in, double, *Nbins);
     671  ALLOCATE(Xbin_in,double,*Nbins);
     672  ALLOCATE(WXbin_in,double,*Nbins);
     673  ALLOCATE(Ybin_in,double,*Nbins);
     674  ALLOCATE(WYbin_in,double,*Nbins);
     675
     676  k = 0;
     677  for (j = 0; j < Nbins_test; j++) {
     678    if (bin_test[j] == 0) { // No data for this bin
     679      continue;
     680    }
     681    else {
     682      // Allocate internal arrays to hold the data that goes into this bin.
     683      ALLOCATE(T_int,double, bin_test[j]);
     684      ALLOCATE(pR_int,double, bin_test[j]);
     685      ALLOCATE(pD_int,double, bin_test[j]);
     686      ALLOCATE(X_int,double, bin_test[j]);
     687      ALLOCATE(Y_int,double, bin_test[j]);
     688      ALLOCATE(WX_int,double, bin_test[j]);
     689      ALLOCATE(WY_int,double, bin_test[j]);
     690
     691      // Fill those arrays.
     692      l = 0;
     693      for (i = 0; i < Npts; i++) {
     694        if ((T[i] >= T_min + j * binning_step)&&
     695            (T[i] <  T_min + (j + 1) * binning_step)) {
     696          T_int[l] = T[i];
     697          pR_int[l] = pR[i];
     698          pD_int[l] = pD[i];
     699          X_int[l] = X[i];
     700          Y_int[l] = Y[i];
     701          WX_int[l] = WX[i];
     702          WY_int[l] = WY[i];
     703          l++;
     704          //      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],
     705          //     T_min + j * binning_step,T_min + (j+1) * binning_step);
     706        }
     707      } // End loop over input points
     708
     709      // Make a decision what the binned values should be.
     710      if (l == 1) {
     711        Tbin_in[k] = T_int[0];
     712        Xbin_in[k] = X_int[0];
     713        Ybin_in[k] = Y_int[0];
     714        pRbin_in[k] = pR_int[0];
     715        pDbin_in[k] = pD_int[0];
     716        WXbin_in[k] = WX_int[0];
     717        WYbin_in[k] = WY_int[0];
     718      }
     719      else {
     720        // I think I'm going with medians.
     721        dsort(T_int,l);
     722        dsort(pR_int,l);
     723        dsort(pD_int,l);
     724        dsort(X_int,l);
     725        dsort(Y_int,l);
     726       
     727        // The median gives the midpoint of all the measurements
     728        if ((l % 2) == 0) {
     729          Tbin_in[k] = 0.5*(T_int[(int)(0.5*l)] + T_int[(int)(0.5*l) - 1]);
     730          pRbin_in[k] = 0.5*(pR_int[(int)(0.5*l)] + pR_int[(int)(0.5*l) - 1]);
     731          pDbin_in[k] = 0.5*(pD_int[(int)(0.5*l)] + pD_int[(int)(0.5*l) - 1]);
     732          Xbin_in[k] = 0.5*(X_int[(int)(0.5*l)] + X_int[(int)(0.5*l) - 1]);
     733          Ybin_in[k] = 0.5*(Y_int[(int)(0.5*l)] + Y_int[(int)(0.5*l) - 1]);
     734        } else {
     735          Tbin_in[k] = T_int[(int)(0.5*l)];
     736          pDbin_in[k] = pD_int[(int)(0.5*l)];
     737          pRbin_in[k] = pR_int[(int)(0.5*l)];
     738          Xbin_in[k] = X_int[(int)(0.5*l)];
     739          Ybin_in[k] = Y_int[(int)(0.5*l)];
     740        }
     741       
     742        // The scatter between points is probably the most useful measurement for the error.
     743        for (i = 0; i < l; i++) {
     744          X_int[i] = fabs(X_int[i] - Xbin_in[k]);
     745          Y_int[i] = fabs(Y_int[i] - Ybin_in[k]);
     746        }
     747        dsort(X_int,l);
     748        dsort(Y_int,l);
     749       
     750        if ((l % 2) == 0) {
     751          WXbin_in[k] = 1.0 / (0.5*(X_int[(int)(0.5*l)] + X_int[(int)(0.5*l) - 1]));
     752          WYbin_in[k] = 1.0 / (0.5*(Y_int[(int)(0.5*l)] + Y_int[(int)(0.5*l) - 1]));
     753        } else {
     754          WXbin_in[k] = 1.0 / (X_int[(int)(0.5*l)]);
     755          WYbin_in[k] = 1.0 / (Y_int[(int)(0.5*l)]);
     756        }
     757      }
     758     
     759      if (WXbin_in[k] == 0.0) {
     760        WXbin_in[k] = WX_int[0];
     761      }
     762      if (WYbin_in[k] == 0.0) {
     763        WYbin_in[k] = WY_int[0];
     764      }
     765      if (!isfinite(WXbin_in[k])) {
     766        WXbin_in[k] = 1.0;
     767      }
     768      if (!isfinite(WYbin_in[k])) {
     769        WYbin_in[k] = 1.0;
     770      }
     771     
     772      // Increment bin index.
     773      k++;
     774     
     775      // Clean up
     776      FREE(T_int);
     777      FREE(pR_int);
     778      FREE(pD_int);
     779      FREE(X_int);
     780      FREE(Y_int);
     781      FREE(WX_int);
     782      FREE(WY_int);         
     783    }
     784  } // End loop over initial bin set.
     785
     786  *Tbin = Tbin_in;
     787  *pRbin = pRbin_in;
     788  *pDbin = pDbin_in;
     789  *Xbin = Xbin_in;
     790  *Ybin = Ybin_in;
     791  *WXbin = WXbin_in;
     792  *WYbin = WYbin_in;
     793
     794  FREE(T_tmp);
     795  FREE(bin_test);
     796 
     797  return TRUE;
     798}
     799

  • trunk/Ohana/src/opihi/cmd.astro/fitpm_irls.c

    r39228 r39585  
    2525  }
    2626
     27  double binning_step = 0.0;
     28  if ((N = get_argument (argc, argv, "-binning"))) {
     29    remove_argument (N, &argc, argv);
     30    binning_step = atof(argv[N]);
     31    remove_argument (N, &argc, argv);
     32  }
     33 
     34 
    2735  if (argc != 6) {
    28     gprint (GP_ERR, "USAGE: fitplx_irls (ra) (dR) (dec) (dD) (mjd) [-mask mask]\n");
     36    gprint (GP_ERR, "USAGE: fitplx_irls (ra) (dR) (dec) (dD) (mjd) [-mask mask] [-binning step_size]\n");
    2937    // what about the errors?
    3038    return (FALSE);
    … …  
    105113
    106114  PlxFit fit;
    107   if (!FitPMonly_IRLS (&fit, X, dX, Y, dY, t, n, VERBOSE)) {
     115  if (!FitPMonly_IRLS (&fit, X, dX, Y, dY, t, n, binning_step, VERBOSE)) {
    108116    return FALSE;
    109117  }
    … …  
    140148
    141149/* do we want an init function which does the alloc and a clear function to free? */
    142 int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) {
     150int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, double binning_step, int VERBOSE) {
    143151
    144152  int i,j;
    … …  
    196204  }
    197205 
    198   // Solve OLS equation 
    199   if (!weighted_LS_PM(T,X,Wx,Y,Wy,Npts,
    200                    A,B,VERBOSE)) {
    201     // Handle fail case
    202     return(FALSE);
     206  // Solve OLS equation
     207  if (binning_step == 0.0) {
     208    if (!weighted_LS_PM(T,X,Wx,Y,Wy,Npts,
     209                        A,B,VERBOSE)) {
     210      // Handle fail case
     211      return(FALSE);
     212    }
     213  }
     214  else {
     215    int Nbins;
     216    double *Tbin = NULL;
     217    double *Xbin = NULL;
     218    double *Ybin = NULL;
     219    double *WXbin = NULL;
     220    double *WYbin = NULL;
     221   
     222    if (!bin_points(T, X, Wx, Y, Wy, Npts,
     223                    &Tbin, &Xbin, &WXbin, &Ybin, &WYbin, &Nbins,  binning_step) ) {
     224      return(FALSE);
     225    }
     226    if (!weighted_LS_PM(Tbin,Xbin,WXbin,Ybin,WYbin, Nbins,
     227                        A,B,VERBOSE)) {
     228      return(FALSE);
     229    }
     230    FREE(Tbin);
     231    FREE(Xbin);
     232    FREE(Ybin);
     233    FREE(WXbin);
     234    FREE(WYbin);
    203235  }
    204236
    … …  
    424456}
    425457
     458int bin_points (double *T, double *X, double *WX, double *Y, double *WY, int Npts,
     459                double **Tbin, double **Xbin, double **WXbin, double **Ybin, double **WYbin, int *Nbins, double binning_step) {
     460  double *T_tmp = NULL;
     461  int i,j,k,l;
     462
     463  int Nbins_test;
     464  int *bin_test = NULL;
     465
     466  double *T_int = NULL, *X_int = NULL, *Y_int = NULL, *WX_int = NULL, *WY_int = NULL;
     467  double T_min, T_max;
     468
     469
     470  //  printf("In binning code\n");
     471
     472  // Allocate more bins than we need
     473  ALLOCATE(T_tmp,double,Npts);
     474  for (i = 0; i < Npts; i++) {
     475    T_tmp[i] = T[i];
     476  }
     477  dsort(T_tmp,Npts);
     478
     479  Nbins_test = floor((T_tmp[Npts - 1] - T_tmp[0]) / binning_step) + 1;
     480  ALLOCATE(bin_test,int,Nbins_test);
     481
     482  T_min = T_tmp[0];
     483  T_max = T_tmp[Npts - 1];
     484
     485  // printf("Found %d points, putting into %d bins between %f and %f\n",Npts,Nbins_test,T_min,T_max);
     486  // Check how many bins are filled.
     487  for (j = 0; j < Nbins_test; j++) {
     488    bin_test[j] = 0;
     489  }
     490  for (i = 0; i < Npts; i++) {
     491    j = floor((T_tmp[i] - T_tmp[0]) / binning_step);
     492    bin_test[j] ++;
     493  }
     494
     495  *Nbins = 0;
     496  for (j = 0; j < Nbins_test; j++) {
     497    if (bin_test[j] != 0) {
     498      *Nbins = *Nbins + 1;
     499    }
     500  }
     501
     502  // printf("Using %d actual bins\n",*Nbins);
     503
     504  double *Tbin_in = NULL;
     505  double *Xbin_in = NULL;
     506  double *Ybin_in = NULL;
     507  double *WXbin_in = NULL;
     508  double *WYbin_in = NULL;
     509 
     510  ALLOCATE(Tbin_in,double,*Nbins);
     511  ALLOCATE(Xbin_in,double,*Nbins);
     512  ALLOCATE(WXbin_in,double,*Nbins);
     513  ALLOCATE(Ybin_in,double,*Nbins);
     514  ALLOCATE(WYbin_in,double,*Nbins);
     515
     516  k = 0;
     517  for (j = 0; j < Nbins_test; j++) {
     518    if (bin_test[j] == 0) { // No data for this bin
     519      continue;
     520    }
     521    else {
     522      // printf("%d bin %d, N = %d\n",k,j,bin_test[j]);
     523      // Allocate internal arrays to hold the data that goes into this bin.
     524      ALLOCATE(T_int,double, bin_test[j]);
     525      ALLOCATE(X_int,double, bin_test[j]);
     526      ALLOCATE(Y_int,double, bin_test[j]);
     527      ALLOCATE(WX_int,double, bin_test[j]);
     528      ALLOCATE(WY_int,double, bin_test[j]);
     529
     530      // Fill those arrays.
     531      l = 0;
     532      for (i = 0; i < Npts; i++) {
     533        if ((T[i] >= T_min + j * binning_step)&&
     534            (T[i] <  T_min + (j + 1) * binning_step)) {
     535          T_int[l] = T[i];
     536          X_int[l] = X[i];
     537          Y_int[l] = Y[i];
     538          WX_int[l] = WX[i];
     539          WY_int[l] = WY[i];
     540          // printf("%d %d %f %f %f %f %f\n",i,l,T[i],X[i],Y[i],WX[i],WY[i]);
     541          l++;
     542        }
     543      } // End loop over input points
     544      // printf("%d bin %d, N = %d, l = %d\n",k,j,bin_test[j],l);
     545
     546      // Make a decision what the binned values should be.
     547      if (l == 1) {
     548        Tbin_in[k] = T_int[0];
     549        Xbin_in[k] = X_int[0];
     550        Ybin_in[k] = Y_int[0];
     551        WXbin_in[k] = WX_int[0];
     552        WYbin_in[k] = WY_int[0];
     553      }
     554      else {
     555        // I think I'm going with medians.
     556        dsort(T_int,l);
     557        dsort(X_int,l);
     558        dsort(Y_int,l);
     559       
     560        // The median gives the midpoint of all the measurements
     561        if ((l % 2) == 0) {
     562          Tbin_in[k] = 0.5*(T_int[(int)(0.5*l)] + T_int[(int)(0.5*l) - 1]);
     563          Xbin_in[k] = 0.5*(X_int[(int)(0.5*l)] + X_int[(int)(0.5*l) - 1]);
     564          Ybin_in[k] = 0.5*(Y_int[(int)(0.5*l)] + Y_int[(int)(0.5*l) - 1]);
     565        } else {
     566          Tbin_in[k] = T_int[(int)(0.5*l)];
     567          Xbin_in[k] = X_int[(int)(0.5*l)];
     568          Ybin_in[k] = Y_int[(int)(0.5*l)];
     569        }
     570       
     571        // The scatter between points is probably the most useful measurement for the error.
     572        for (i = 0; i < l; i++) {
     573          X_int[i] = fabs(X_int[i] - Xbin_in[k]);
     574          Y_int[i] = fabs(Y_int[i] - Ybin_in[k]);
     575        }
     576        dsort(X_int,l);
     577        dsort(Y_int,l);
     578       
     579        if ((l % 2) == 0) {
     580          WXbin_in[k] = 1.0 / (0.5*(X_int[(int)(0.5*l)] + X_int[(int)(0.5*l) - 1]));
     581          WYbin_in[k] = 1.0 / (0.5*(Y_int[(int)(0.5*l)] + Y_int[(int)(0.5*l) - 1]));
     582        } else {
     583          WXbin_in[k] = 1.0 / (X_int[(int)(0.5*l)]);
     584          WYbin_in[k] = 1.0 / (Y_int[(int)(0.5*l)]);
     585        }
     586      }
     587
     588      if (WXbin_in[k] == 0.0) {
     589        WXbin_in[k] = WX_int[0];
     590      }
     591      if (WYbin_in[k] == 0.0) {
     592        WYbin_in[k] = WY_int[0];
     593      }
     594      if (!isfinite(WXbin_in[k])) {
     595        WXbin_in[k] = 1.0;
     596      }
     597      if (!isfinite(WYbin_in[k])) {
     598        WYbin_in[k] = 1.0;
     599      }
     600     
     601      //      printf("%f %f %f %f %f\n",Tbin_in[k],Xbin_in[k],Ybin_in[k],WXbin_in[k],WYbin_in[k]);
     602      // Increment bin index.
     603      k++;
     604     
     605      // Clean up
     606      FREE(T_int);
     607      FREE(X_int);
     608      FREE(Y_int);
     609      FREE(WX_int);
     610      FREE(WY_int);
     611    }
     612  } // End loop over initial bin set.
     613
     614  *Tbin = Tbin_in;
     615  *Xbin = Xbin_in;
     616  *Ybin = Ybin_in;
     617  *WXbin = WXbin_in;
     618  *WYbin = WYbin_in;
     619
     620  FREE(T_tmp);
     621  FREE(bin_test);
     622 
     623  return TRUE;
     624}
     625
     626
    426627double weight_cauchy (double x) {
    427628  double r = x / 2.385;
    … …  
    450651  dsort(x,N);
    451652
    452   if (N % 2) {
     653  if ((N % 2) == 0) {
    453654    median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
    454655  } else {
    … …  
    462663  dsort(x,N);
    463664
    464   if (N % 2) {
     665  if ((N % 2) == 0) {
    465666    median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
    466667  } else {
    … …  
    468669  }
    469670
     671  FREE(x);
    470672  return(median);
    471673}

  • trunk/Ohana/src/opihi/include/astro.h

    r39229 r39585  
    5858/***** */
    5959
    60 int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE);
    61 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);
     60int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, double binning_step, int VERBOSE);
     61int 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);
    6262
    6363int 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);
    6464int weighted_LS_PM (double *T, double *X, double *WX, double *Y, double *WY, int Npts, double **A, double **B, int VERBOSE);
     65int bin_points (double *T, double *X, double *WX, double *Y, double *WY, int Npts,
     66                double **Tbin, double **Xbin, double **WXbin, double **Ybin, double **WYbin, int *Nbins, double binning_step);
     67int bin_points_PLX (double *T, double *pR, double *pD, double *X, double *WX, double *Y, double *WY, int Npts,
     68                    double **Tbin, double **pRbin, double **pDbin, double **Xbin, double **WXbin, double **Ybin, double **WYbin, int *Nbins, double binning_step);
     69
    6570
    6671double weight_cauchy (double x);
Note: See TracChangeset for help on using the changeset viewer.