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));
}

lls_workspace *
lls_alloc(const size_t max_block, const size_t p)
{
  const gsl_multifit_robust_type *robust_t = gsl_multifit_robust_huber;
  const size_t nblock = GSL_MAX(max_block, p);
  lls_workspace *w;

  w = calloc(1, sizeof(lls_workspace));
  if (!w)
    return 0;

  /* A^T A matrix */
  w->ATA = gsl_matrix_alloc(p, p);
  w->work_A = gsl_matrix_alloc(p, p);
  if (!w->ATA || !w->work_A)
    {
      lls_free(w);
      return 0;
    }

  w->c = gsl_vector_alloc(p);
  w->ATb = gsl_vector_alloc(p);
  w->work_b = gsl_vector_alloc(p);
  w->AF = gsl_matrix_alloc(p, p);
  w->S = gsl_vector_alloc(p);
  if (!w->ATb || !w->work_b)
    {
      lls_free(w);
      return 0;
    }

  w->r = gsl_vector_alloc(nblock);
  w->w_robust = gsl_vector_alloc(nblock);
  w->robust_workspace_p = gsl_multifit_robust_alloc(robust_t, nblock, p);
  if (!w->w_robust || !w->r)
    {
      lls_free(w);
      return 0;
    }

  w->eigen_p = gsl_eigen_symm_alloc(p);
  w->eval = gsl_vector_alloc(p);

  w->p = p;
  w->max_block = nblock;

  w->residual = 0.0;
  w->chisq = 0.0;
  w->cond = -1.0;
  w->niter = 0;
  w->bTb = 0.0;

  /* initialize matrices/vectors to 0 */
  lls_reset(w);

  return w;
} /* lls_alloc() */

int
gsl_multilarge_nlinear_test (const double xtol, const double gtol,
                             const double ftol, int *info,
                             const gsl_multilarge_nlinear_workspace * w)
{
  int status;
  double gnorm, fnorm, phi;

  *info = 0;

  status = gsl_multifit_test_delta(w->dx, w->x, xtol*xtol, xtol);
  if (status == GSL_SUCCESS)
    {
      *info = 1;
      return GSL_SUCCESS;
    }

  /* compute gnorm = max_i( g_i * max(x_i, 1) ) */
  gnorm = scaled_infnorm(w->x, w->g);

  /* compute fnorm = ||f|| */
  fnorm = gsl_blas_dnrm2(w->f);
  phi = 0.5 * fnorm * fnorm;

#if 0
  fprintf(stderr, "gnorm = %.12e fnorm = %.12e gnorm/phi = %.12e\n", gnorm, fnorm, gnorm / phi);
#endif

  if (gnorm <= gtol * GSL_MAX(phi, 1.0))
    {
      *info = 2;
      return GSL_SUCCESS;
    }

#if 0
  if (dfnorm <= ftol * GSL_MAX(fnorm, 1.0))
    {
      *info = 3;
      return GSL_SUCCESS;
    }
#endif

  return GSL_CONTINUE;
}

static void
_nc_cluster_mass_benson_p_limits (NcClusterMass *clusterm, NcHICosmo *model, const gdouble *xi, const gdouble *xi_params, gdouble *lnM_lower, gdouble *lnM_upper)
{
  NcClusterMassBenson *msz = NC_CLUSTER_MASS_BENSON (clusterm);
  const gdouble xil = GSL_MAX (xi[0] - 7.0, msz->signif_obs_min);
  const gdouble zetal = _significance_to_zeta (clusterm, model, 2.0, xil) - 7.0 * D_SZ;
  const gdouble lnMl = GSL_MAX (_zeta_to_mass (clusterm, model, 2.0, zetal), log (NC_CLUSTER_MASS_BENSON_M_LOWER_BOUND));

  const gdouble xiu = xi[0] + 7.0;
  const gdouble zetau = _significance_to_zeta (clusterm, model, 0.0, xiu) + 7.0 * D_SZ;
  const gdouble lnMu = _zeta_to_mass (clusterm, model, 0.0, zetau);

  NCM_UNUSED (xi_params);
  
  *lnM_lower = lnMl;
  *lnM_upper = lnMu;

  return;
}

static double
cholesky_LDLT_norm1(const gsl_matrix * LDLT, const gsl_permutation * p, gsl_vector * work)
{
  const size_t N = LDLT->size1;
  gsl_vector_const_view D = gsl_matrix_const_diagonal(LDLT);
  gsl_vector_view diagA = gsl_vector_subvector(work, N, N);
  double max = 0.0;
  size_t i, j;

  /* reconstruct diagonal entries of original matrix A */
  for (j = 0; j < N; ++j)
    {
      double Ajj;

      /* compute diagonal (j,j) entry of A */
      Ajj = gsl_vector_get(&D.vector, j);
      for (i = 0; i < j; ++i)
        {
          double Di = gsl_vector_get(&D.vector, i);
          double Lji = gsl_matrix_get(LDLT, j, i);

          Ajj += Di * Lji * Lji;
        }

      gsl_vector_set(&diagA.vector, j, Ajj);
    }

  gsl_permute_vector_inverse(p, &diagA.vector);

  for (j = 0; j < N; ++j)
    {
      double sum = 0.0;
      double Ajj = gsl_vector_get(&diagA.vector, j);

      for (i = 0; i < j; ++i)
        {
          double *wi = gsl_vector_ptr(work, i);
          double Aij = gsl_matrix_get(LDLT, i, j);
          double absAij = fabs(Aij);

          sum += absAij;
          *wi += absAij;
        }

      gsl_vector_set(work, j, sum + fabs(Ajj));
    }

  for (i = 0; i < N; ++i)
    {
      double wi = gsl_vector_get(work, i);
      max = GSL_MAX(max, wi);
    }

  return max;
}

/** If there is an NaN anywhere in the row of data (including the matrix, the vector, the weights, and the text) then delete the row from the data set.

\li If every row has an NaN, then this returns \c NULL.
\li If \c apop_opts.db_nan is not \c NULL, then I will use that as a regular expression to check the text elements for bad data as well.
\li If \c inplace = 'y', then I'll free each element of the input data
    set and refill it with the pruned elements. I'll still take up (up to)
    twice the size of the data set in memory during the function. If
    every row has an NaN, then your \c apop_data set will end up with
    \c NULL vector, matrix, \dots. if \c inplace = 'n', then the original data set is left unmolested.
\li I only look at the first page of data (i.e. the \c more element is ignored).
\li This function uses the \ref designated syntax for inputs.

    \param d    The data, with NaNs
    \param inplace If \c 'y', clear out the pointer-to-\ref apop_data that
    you sent in and refill with the pruned data. If \c 'n', leave the
    set alone and return a new data set.
    \return     A (potentially shorter) copy of the data set, without
    NaNs. If <tt>inplace=='y'</tt>, redundant with the input. If the entire data set is
    cleared out, then this will be \c NULL.
*/
APOP_VAR_HEAD apop_data * apop_data_listwise_delete(apop_data *d, char inplace){
    apop_data * apop_varad_var(d, NULL);
    if (!d) return NULL;
    char apop_varad_var(inplace, 'n');
APOP_VAR_ENDHEAD
    Get_vmsizes(d) //defines firstcol, vsize, wsize, msize1, msize2.
    apop_assert_c(msize1 || vsize || d->textsize[0], NULL, 0, 
            "You sent to apop_data_listwise_delete a data set with NULL matrix, NULL vector, and no text. "
            "Confused, it is returning NULL.");
    //find out where the NaNs are
    int len = GSL_MAX(vsize ? vsize : msize1, d->textsize[0]); //still some size assumptions here.
    int not_empty = 0;
    int *marked = calloc(len, sizeof(int));
    for (int i=0; i< (vsize ? vsize: msize1); i++)
        for (int j=firstcol; j <msize2; j++){
            if (gsl_isnan(apop_data_get(d, i, j))){
                    marked[i] = 1;
                    break;
            }
        }
    for (int i=0; i< wsize; i++)
        if (gsl_isnan(gsl_vector_get(d->weights, i)))
            marked[i] = 1;
    if (d->textsize[0] && apop_opts.db_nan){
        regex_t    rex;
        int compiled_ok = !regcomp(&rex, apop_opts.db_nan, REG_EXTENDED +  REG_ICASE + REG_NOSUB);
        apop_assert(compiled_ok, "apop_opts.db_nan needs to be a regular expression that "
                                "I can use to check the text element of your data set for "
                                "NaNs, But compiling %s into a regex failed. Or, set "
                                "apop_opts.db_nan=NULL to bypass text checking.", apop_opts.db_nan);
        for(int i=0; i< d->textsize[0]; i++)
            if (!marked[i])
                for(int j=0; j< d->textsize[1]; j++)
                    if (!regexec(&rex, d->text[i][j], 0, 0, 0)){
                        marked[i] ++;
                        break;
                    }
        regfree(&rex);
    }

    //check that at least something isn't NULL.
    for (int i=0; i< len; i++)
        if (!marked[i]){
            not_empty ++;
            break;
        }
    if (!not_empty){
        free(marked);
        return NULL;
    }
    apop_data *out = (inplace=='y'|| inplace=='Y') ? d : apop_data_copy(d);
    apop_data_rm_rows(out, marked);
    free(marked);
    return out;
}