C++ GSL_MAX函数代码示例

int bbfind_compute(bbfind*bbf, BB2 bbox) {
	double ul[2], ur[2], ll[2], lr[2];
	
	if(1) {
		if(!getBoundingBox(bbf->buf, bbf->num, ul, ur, ll, lr)) {
			sm_error("Could not compute bounding box.\n");
			return 0;
		}
		bbox->pose[0] = ll[0];
		bbox->pose[1] = ll[1];
		bbox->pose[2] = atan2(lr[1]-ll[1], lr[0]-ll[0]);
		bbox->size[0] = distance_d(lr, ll);
		bbox->size[1] = distance_d(ll, ul);
	} else {
		double bb_min[2] = {bbf->buf[0].x,bbf->buf[0].y}, 
				bb_max[2] = {bbf->buf[0].x,bbf->buf[0].y};
		int i; for(i=0;i<bbf->num; i++) {
			bb_min[0] = GSL_MIN(bb_min[0], bbf->buf[i].x);
			bb_min[1] = GSL_MIN(bb_min[1], bbf->buf[i].y);
			bb_max[0] = GSL_MAX(bb_max[0], bbf->buf[i].x);
			bb_max[1] = GSL_MAX(bb_max[1], bbf->buf[i].y);
		}
		bbox->pose[0] = bb_min[0];
		bbox->pose[1] = bb_min[1];
		bbox->pose[2] = 0;
		bbox->size[0] = bb_max[0] - bb_min[0];
		bbox->size[1] = bb_max[1] - bb_min[1];
	}
	return 1;
}

/** Run the Kolmogorov-Smirnov test to determine whether two distributions are identical.

\param m1, m2  Two models, most likely of \ref apop_pmf type. I will ue the cdf method, so if your function doesn't have one, expect this to run the slow default. I run it for each row of each data set, so if your model has a \c NULL at the data, I won't know what to do. 
 
\return An \ref apop_data set including the \f$p\f$-value from the Kolmogorov test that the two distributions are equal.

\li I assume that the data sets are sorted.

\include ks_tests.c

\ingroup histograms
*/
apop_data *apop_test_kolmogorov(apop_model *m1, apop_model *m2){
    //version for not a pair of histograms
    Apop_assert(m1->data, "I will test the CDF at each point in the data set, but the first model has a NULL data set. "
                          "Maybe generate, then apop_data_sort, a few thousand random draws?");
    Apop_assert(m2->data, "I will test the CDF at each point in the data set, but the second model has a NULL data set. "
                          "Maybe generate, then apop_data_sort, a few thousand random draws?");
    int maxsize1, maxsize2;
    {Get_vmsizes(m1->data); maxsize1 = maxsize;}//copy one of the macro's variables 
    {Get_vmsizes(m2->data); maxsize2 = maxsize;}//  to the full function's scope.
    double largest_diff=GSL_NEGINF;
    for (size_t i=0; i< maxsize1; i++){
        Apop_data_row(m1->data, i, arow);
        largest_diff = GSL_MAX(largest_diff, fabs(apop_cdf(arow, m1)-apop_cdf(arow, m2)));
    }
    for (size_t i=0; i< maxsize2; i++){     //There should be matched data rows, so there is redundancy.
        Apop_data_row(m2->data, i, arow); // Feel free to submit a smarter version.
        largest_diff = GSL_MAX(largest_diff, fabs(apop_cdf(arow, m1)-apop_cdf(arow, m2)));
    }
    apop_data *out = apop_data_alloc();
    sprintf(out->names->title, "Kolmogorov-Smirnov test");
    apop_data_add_named_elmt(out, "max distance", largest_diff);
    double ps = psmirnov2x(largest_diff, maxsize1, maxsize2);
    apop_data_add_named_elmt(out, "p value, 2 tail", 1-ps);
    apop_data_add_named_elmt(out, "confidence, 2 tail", ps);
    return out;
}

inline
static double locMAX4(double x, double y, double z, double w)
{
  double xy  = GSL_MAX(x,  y);
  double xyz = GSL_MAX(xy, z);
  return GSL_MAX(xyz, w);
}

inline
static
int locMax3(const int a, const int b, const int c)
{
  int d = GSL_MAX(a, b);
  return GSL_MAX(d, c);
}

static int
magcal_init(const satdata_mag *data, magcal_workspace *w)
{
  int s = 0;
  size_t i;
  size_t n = 0;

  for (i = 0; i < data->n; ++i)
    {
      /* don't store flagged data */
      if (data->flags[i])
        continue;

      /* don't process high latitude data */
      if (fabs(data->latitude[i]) > MAGCAL_MAX_LATITUDE)
        continue;

      w->Ex[n] = SATDATA_VEC_X(data->B_VFM, i);
      w->Ey[n] = SATDATA_VEC_Y(data->B_VFM, i);
      w->Ez[n] = SATDATA_VEC_Z(data->B_VFM, i);

      w->F[n] = data->F[i];

      ++n;
    }

  if (n < 200)
    {
      fprintf(stderr, "magcal_init: insufficient data points for calibration: %zu\n",
              n);
      return -1;
    }

  if (n != w->n)
    {
      gsl_multifit_fdfsolver_free(w->fdf_s);
      gsl_multifit_fdfridge_free(w->fdf_ridge);
      w->fdf_s = gsl_multifit_fdfsolver_alloc(w->fdf_type, n, w->p);
      w->fdf_ridge = gsl_multifit_fdfridge_alloc(w->fdf_type, n, w->p);
      w->n = n;
    }

#if MAGCAL_SCALE
  w->B_s = GSL_MAX(gsl_stats_sd(w->Ex, 1, n),
                   GSL_MAX(gsl_stats_sd(w->Ey, 1, n),
                           gsl_stats_sd(w->Ez, 1, n)));
#endif

  /* center and scale data arrays */
  for (i = 0; i < n; ++i)
    {
      w->Ex[i] /= w->B_s;
      w->Ey[i] /= w->B_s;
      w->Ez[i] /= w->B_s;
      w->F[i] /= w->B_s;
    }

  return s;
} /* magcal_init() */

 apop_data * apop_bootstrap_cov_base(apop_data * data, apop_model *model, gsl_rng *rng, int iterations, char keep_boots, char ignore_nans, apop_data **boot_store){
#endif
    Get_vmsizes(data); //vsize, msize1, msize2
    apop_model *e = apop_model_copy(model);
    apop_data *subset = apop_data_copy(data);
    apop_data *array_of_boots = NULL,
              *summary;
    //prevent and infinite regression of covariance calculation.
    Apop_model_add_group(e, apop_parts_wanted); //default wants for nothing.
    size_t i, nan_draws=0;
    apop_name *tmpnames = (data && data->names) ? data->names : NULL; //save on some copying below.
    if (data && data->names) data->names = NULL;

    int height = GSL_MAX(msize1, GSL_MAX(vsize, (data?(*data->textsize):0)));
	for (i=0; i<iterations && nan_draws < iterations; i++){
		for (size_t j=0; j< height; j++){       //create the data set
			size_t randrow	= gsl_rng_uniform_int(rng, height);
            apop_data_memcpy(Apop_r(subset, j), Apop_r(data, randrow));
		}
		//get the parameter estimates.
		apop_model *est = apop_estimate(subset, e);
        gsl_vector *estp = apop_data_pack(est->parameters);
        if (!gsl_isnan(apop_sum(estp))){
            if (i==0){
                array_of_boots	      = apop_data_alloc(iterations, estp->size);
                apop_name_stack(array_of_boots->names, est->parameters->names, 'c', 'v');
                apop_name_stack(array_of_boots->names, est->parameters->names, 'c', 'c');
                apop_name_stack(array_of_boots->names, est->parameters->names, 'c', 'r');
            }
            gsl_matrix_set_row(array_of_boots->matrix, i, estp);
        } else if (ignore_nans=='y'){
            i--; 
            nan_draws++;
        }
        apop_model_free(est);
        gsl_vector_free(estp);
	}
    if(data) data->names = tmpnames;
    apop_data_free(subset);
    apop_model_free(e);
    int set_error=0;
    Apop_stopif(i == 0 && nan_draws == iterations, apop_return_data_error(N),
                1, "I ran into %i NaNs and no not-NaN estimations, and so stopped. "
                       , iterations);
    Apop_stopif(nan_draws == iterations,  set_error++;
            apop_matrix_realloc(array_of_boots->matrix, i, array_of_boots->matrix->size2),
                1, "I ran into %i NaNs, and so stopped. Returning results based "
                       "on %zu bootstrap iterations.", iterations, i);
	summary	= apop_data_covariance(array_of_boots);
    if (boot_store) *boot_store = array_of_boots;
    else            apop_data_free(array_of_boots);
    if (set_error) summary->error = 'N';
	return summary;
}

static int
secs2d_green_df(const double r, const double theta, const double phi,
                const double theta0, const double phi0,
                double B[3], secs2d_state_t *state)
{
  const double s = GSL_MIN(r, state->R) / GSL_MAX(r, state->R);
  double thetap, stp, ctp;
  size_t i;

  secs2d_transform(theta0, phi0, theta, phi, &thetap);
  stp = sin(thetap);
  ctp = cos(thetap);

  if (r >= state->R)
    {
      B[0] = ((1.0 - s * ctp) / sqrt(1 + s*s - 2.0*s*ctp) - 1.0) / stp;
      B[1] = 0.0;
      B[2] = -s * (1.0 / sqrt(1.0 + s*s - 2.0*s*ctp) - 1.0);
    }

  for (i = 0; i < 3; ++i)
    B[i] *= MAGFIT_MU_0 / (4.0 * M_PI * r);

  /* transform from SECS-centered system to geographic */
  secs2d_transform_vec(theta0, phi0, theta, phi, B, B);

  return GSL_SUCCESS;
}

long AT_n_bins_for_DSB_distribution(const long n_bins_f,
		const double f_d_Gy[],
		const double f_dd_Gy[],
		const double f[],
		const double enhancement_factor[],
		const double DSB_per_Gy_per_domain){


	  //double* avg_DSB = (double*)calloc(n_bins_f, sizeof(double));

	  double max_number_of_DSB = 0.0;

	  for(long i = 0; i < n_bins_f; i++){
		  max_number_of_DSB = GSL_MAX(max_number_of_DSB, f_d_Gy[i] * enhancement_factor[i] * DSB_per_Gy_per_domain);
	  }

	  if(isnan(max_number_of_DSB)){
	    return(0);
	  }

	  // Use max + 5 * st.dev as safety margin
	  max_number_of_DSB = floor(max_number_of_DSB) + 1.0;
	  max_number_of_DSB += floor(5 * sqrt(max_number_of_DSB)) + 1;

	  // Add one for zero bin
	  return((long)max_number_of_DSB + 1);
}

size_t
track_flag_incomplete(const double qd_min, const double qd_max, satdata_mag *data, track_workspace *w)
{
  size_t i;
  size_t nflagged = 0;        /* number of points flagged */
  size_t ntrack_flagged = 0;  /* number of entire tracks flagged for missing data */

  if (data->n == 0)
    return 0;

  for (i = 0; i < w->n; ++i)
    {
      track_data *tptr = &(w->tracks[i]);
      double qd0 = GSL_MIN(data->qdlat[tptr->start_idx], data->qdlat[tptr->end_idx]);
      double qd1 = GSL_MAX(data->qdlat[tptr->start_idx], data->qdlat[tptr->end_idx]);

      if (qd0 > qd_min || qd1 < qd_max)
        {
          nflagged += track_flag_track(i, TRACK_FLG_INCOMPLETE, data, w);
          ++ntrack_flagged;
        }
    }

  fprintf(stderr, "track_flag_incomplete: flagged %zu/%zu (%.1f%%) tracks due to missing data\n",
          ntrack_flagged, w->n, (double) ntrack_flagged / (double) w->n * 100.0);

  return nflagged;
} /* track_flag_incomplete() */

int
gsl_multifit_linear_lreg (const double smin, const double smax,
                          gsl_vector * reg_param)
{
  if (smax <= 0.0)
    {
      GSL_ERROR("smax must be positive", GSL_EINVAL);
    }
  else
    {
      const size_t N = reg_param->size;

      /* smallest regularization parameter */
      const double smin_ratio = 16.0 * GSL_DBL_EPSILON;
      const double new_smin = GSL_MAX(smin, smax*smin_ratio);
      double ratio;
      size_t i;

      gsl_vector_set(reg_param, N - 1, new_smin);

      /* ratio so that reg_param(1) = s(1) */
      ratio = pow(smax / new_smin, 1.0 / ((double)N - 1.0));

      /* calculate the regularization parameters */
      for (i = N - 1; i > 0 && i--; )
        {
          double rp1 = gsl_vector_get(reg_param, i + 1);
          gsl_vector_set(reg_param, i, ratio * rp1);
        }

      return GSL_SUCCESS;
    }
}

/* interpolate Ne data near time t to latitude qdlat */
double
interp_ne(const time_t t, const double qdlat, const satdata_mag *data)
{
  double t_cdf = satdata_timet2epoch(t);
  double Ne;

  if (t_cdf >= data->t[0] && t_cdf <= data->t[data->n - 1])
    {
      int idx = (int) bsearch_double(data->t, t_cdf, 0, data->n - 1);
      int half_window = 15; /* number of minutes in half time window */
      int idx1 = GSL_MAX(0, idx - half_window*60*2);
      int idx2 = GSL_MIN(data->n - 1, idx + half_window*60*2);
      int qidx;
      
      if (data->qdlat[idx1] < data->qdlat[idx2])
        qidx = (int) bsearch_double(data->qdlat, qdlat, idx1, idx2);
      else
        qidx = (int) bsearch_desc_double(data->qdlat, qdlat, idx1, idx2);

      Ne = interp1d(data->qdlat[qidx], data->qdlat[qidx + 1],
                    data->ne[qidx], data->ne[qidx + 1], qdlat);
    }
  else
    Ne = 0.0;

  return Ne;
}

int
gsl_sf_pochrel_e(const double a, const double x, gsl_sf_result * result)
{
  const double absx = fabs(x);
  const double absa = fabs(a);

  /* CHECK_POINTER(result) */

  if(absx > 0.1*absa || absx*log(GSL_MAX(absa,2.0)) > 0.1) {
    gsl_sf_result lnpoch;
    double sgn;
    int stat_poch = gsl_sf_lnpoch_sgn_e(a, x, &lnpoch, &sgn);
    if(lnpoch.val > GSL_LOG_DBL_MAX) {
      OVERFLOW_ERROR(result);
    }
    else {
      const double el = exp(lnpoch.val);
      result->val  = (sgn*el - 1.0)/x;
      result->err  = fabs(result->val) * (lnpoch.err + 2.0 * GSL_DBL_EPSILON);
      result->err += 2.0 * GSL_DBL_EPSILON * (fabs(sgn*el) + 1.0) / fabs(x);
      return stat_poch;
    }
  }
  else {
    return pochrel_smallx(a, x, result);
  }
}

int
gsl_sf_ellint_RC_e(double x, double y, gsl_mode_t mode, gsl_sf_result * result)
{
  const double lolim = 5.0 * GSL_DBL_MIN;
  const double uplim = 0.2 * GSL_DBL_MAX;
  const gsl_prec_t goal = GSL_MODE_PREC(mode);
  const double errtol = ( goal == GSL_PREC_DOUBLE ? 0.001 : 0.03 );
  const double prec   = gsl_prec_eps[goal];

  if(x < 0.0 || y < 0.0 || x + y < lolim) {
    DOMAIN_ERROR(result);
  }
  else if(GSL_MAX(x, y) < uplim) { 
    const double c1 = 1.0 / 7.0;
    const double c2 = 9.0 / 22.0;
    double xn = x;
    double yn = y;
    double mu, sn, lamda, s;
    while(1) {
      mu = (xn + yn + yn) / 3.0;
      sn = (yn + mu) / mu - 2.0;
      if (fabs(sn) < errtol) break;
      lamda = 2.0 * sqrt(xn) * sqrt(yn) + yn;
      xn = (xn + lamda) * 0.25;
      yn = (yn + lamda) * 0.25;
    }
    s = sn * sn * (0.3 + sn * (c1 + sn * (0.375 + sn * c2)));
    result->val = (1.0 + s) / sqrt(mu);
    result->err = prec * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    DOMAIN_ERROR(result);
  }
}

static int
lmniel_set(void *vstate, const gsl_vector *swts,
           gsl_multifit_function_fdf *fdf, gsl_vector *x,
           gsl_vector *f, gsl_vector *dx)
{
  int status;
  lmniel_state_t *state = (lmniel_state_t *) vstate;
  const size_t p = x->size;
  size_t i;

  /* initialize counters for function and Jacobian evaluations */
  fdf->nevalf = 0;
  fdf->nevaldf = 0;

  /* evaluate function and Jacobian at x and apply weight transform */
  status = gsl_multifit_eval_wf(fdf, x, swts, f);
  if (status)
   return status;

  if (fdf->df)
    status = gsl_multifit_eval_wdf(fdf, x, swts, state->J);
  else
    status = gsl_multifit_fdfsolver_dif_df(x, swts, fdf, f, state->J);
  if (status)
    return status;

  /* compute rhs = -J^T f */
  gsl_blas_dgemv(CblasTrans, -1.0, state->J, f, 0.0, state->rhs);

#if SCALE
  gsl_vector_set_zero(state->diag);
#else
  gsl_vector_set_all(state->diag, 1.0);
#endif

  /* set default parameters */
  state->nu = 2;

#if SCALE
  state->mu = state->tau;
#else
  /* compute mu_0 = tau * max(diag(J^T J)) */
  state->mu = -1.0;
  for (i = 0; i < p; ++i)
    {
      gsl_vector_view c = gsl_matrix_column(state->J, i);
      double result; /* (J^T J)_{ii} */

      gsl_blas_ddot(&c.vector, &c.vector, &result);
      state->mu = GSL_MAX(state->mu, result);
    }

  state->mu *= state->tau;
#endif

  return GSL_SUCCESS;
} /* lmniel_set() */

static void scaling(size_t const *elmts, size_t const n,  gsl_vector *weights, double const in_sum){
    double fit_sum = 0;
    for(size_t i=0; i < n; i ++)
        fit_sum += weights->data[elmts[i]];
    if (!fit_sum) return; //can happen if init table is very different from margins.
    for(size_t i=0; i < n; i ++)
        weights->data[elmts[i]] *= in_sum/fit_sum;
	overall_max_dev = GSL_MAX(overall_max_dev, fabs(in_sum-fit_sum));
}