amdyv.c

ADaM is a data mining and image processing toolkit

{
  int i;
  if ( d_1 -> size != d_2 -> size )
  {
    fprintf_dyv(stderr,"d_1",d_1,"\n");
    fprintf_dyv(stderr,"d_2",d_2,"\n");
    my_error("dyv_plus: dyvs (DYnamic Vectors) different shape");
  }

  for ( i = 0 ; i < r_d -> size ; i++ )
    r_d -> farr[i] = d_1->farr[i] + d_2 -> farr[i];
}

dyv *mk_dyv_plus(const dyv *a,const dyv *b)
{
  dyv *result = mk_dyv(a->size);
  dyv_plus(a,b,result);
  return(result);
}

void dyv_subtract(const dyv *d_1,const dyv *d_2,dyv *r_d)
{
  int i;
  if ( d_1 -> size != d_2 -> size )
  {
    fprintf_dyv(stderr,"d_1",d_1,"\n");
    fprintf_dyv(stderr,"d_2",d_2,"\n");
    my_error("dyv_subtract: dyvs (DYnamic Vectors) different shape");
  }

  for ( i = 0 ; i < r_d -> size ; i++ )
    r_d -> farr[i] = d_1->farr[i] - d_2 -> farr[i];
}

dyv *mk_dyv_subtract(const dyv *a,const dyv *b)
{
  dyv *result; 
  result = mk_dyv(a->size);
  dyv_subtract(a,b,result);
  return(result);
}

void dyv_abs( const dyv *dv, dyv *r_dv)
{
  int i;
  for (i=0; i<dyv_size( dv); ++i) dyv_set( r_dv, i, fabs(dyv_ref( dv, i)));
  return;
}

dyv *mk_dyv_abs( const dyv *dv)
{
  dyv *absdv;
  absdv = mk_dyv( dyv_size( dv));
  dyv_abs( dv, absdv);
  return absdv;
}

/***** More complex operations ******/

double dyv_scalar_product(const dyv *a,const dyv *b)
{
  int i,n= a->size;
  double result = 0.0;
  if ( b -> size != n)
  {
    fprintf(stderr,"dyv_scalar_product: sizes differ\n");
    wait_for_key();
    fprintf_dyv(stderr,"a",a,"\n");
    fprintf_dyv(stderr,"b",b,"\n");
    my_error("dyv_scalar_product: sizes differ");
  }

  for ( i = 0 ; i < n ; i++ )
    result += a->farr[i] * b->farr[i];
  return(result);
}

double dyv_magnitude(const dyv *a)
{
  double result = sqrt(dyv_scalar_product(a,a));
  return(result);
}

double dyv_mean(const dyv *dv)
{
  return(doubles_mean(dv->farr,dv->size));
}

double dyv_sdev(const dyv *dv)
{
  double sum_sq = 0.0;
  int i;
  double mean = dyv_mean(dv);
  double result;

  for ( i = 0 ; i < dv->size ; i++ )
    sum_sq += real_square(dv->farr[i] - mean);
    
  result = sqrt(sum_sq / int_max(dv->size-1,1));

  return(result);
}

double dyv_min(const dyv *dv)
{
  if ( dyv_size(dv) < 1 )
    my_error("dyv_min: empty dyv");
  return(doubles_min(dv->farr,dv->size));
}

double dyv_max(const dyv *dv)
{
  if ( dyv_size(dv) < 1 )
    my_error("dyv_max: empty dyv");
  return(doubles_max(dv->farr,dv->size));
}

int dyv_argmin(const dyv *dv)
{
  if ( dyv_size(dv) < 1 )
    my_error("dyv_argmin: empty dyv");
  return(doubles_argmin(dv->farr,dv->size));
}

int dyv_argmax(const dyv *dv)
{
  if ( dyv_size(dv) < 1 )
    my_error("dyv_argmax: empty dyv");
  return(doubles_argmax(dv->farr,dv->size));
}

/* Returns first index i such that dv[i] >= key.  If the last element of
   dv is strictly less than the key, then dyv_size(dv) is returned. */
int index_in_sorted_dyv( dyv *dv, double key)
{
  int idxlo, idxhi, idx;
  double val;

  /* idxhi is always >= desired idx.  idxlo is always <= desired idx. */

  idxlo = 0;
  idxhi = dyv_size( dv);
  if (idxhi == 0) return 0;

  while (1) {
    /* Get next value. */
    idx = (idxlo + idxhi) / 2;
    val = dyv_ref( dv, idx);

    /* Adjust bounds. */
    if (val >= key) idxhi = idx ;  /* idx is too high. */
    else idxlo = idx + 1;          /* idx is too low.  */

    /* If the boundaries meet, we're done. */
    if (idxhi == idxlo) break;
  }

  return idxlo;
}


/* We're not quite sure what this code computed.  On the bright
   side, since it isn't commented we can easily comment the
   entire function out. */
/* 
int index_in_sorted_dyv(dyv *d,double t){
  int i1 = dyv_size(d), i0 = 0;
  int i = (i1-i0)>>1;
  double n2 = dyv_ref(d,i), n1 = (i)? dyv_ref(d,i-1):-1;
  while((i1-i0)>1 && !(n2>=t && (!i || n1<t))){
    if(n2>t){
      i1 = i;
      i = i0+((i-i0)>>1);
    } else {
      i0 = i+1;
      i = i+((i1-i)>>1);
    }
  }
  return i;
}
*/

void dyv_sort(const dyv *dv,dyv *r_dv)
{
  int size;
  double *farr;

  size = dyv_size(dv);
  farr = mk_farr_from_dyv(dv);
  sort_realnums(farr,size,farr);
  copy_farr_to_dyv(farr,size,r_dv);
  am_free_realnums(farr,size);
}


dyv *mk_dyv_sort(const dyv *dv)
{
  dyv *result;
  result = mk_dyv(dyv_size(dv));
  dyv_sort(dv,result);
  return(result);
}

/*
 Creates a dyv of indices such that indices[i] is the origional
 location (in the unsorted dv) of the ith smallest value.
 Used when you want the location of the sorted values instead of
 the sorted vector itself.
*/
dyv *mk_sorted_dyv_indices(dyv *dv)
{
  dyv* indices;
  int size = dyv_size(dv);
  int i;
  int* iarr; 
  double* farr;

  farr = mk_farr_from_dyv(dv);
  iarr = am_malloc_ints(size); 

  indices_sort_realnums(farr,size,iarr);
  indices = mk_dyv(size);

  for(i=0;i<size;i++)
  {
    dyv_set(indices,i,iarr[i]);
  }

  am_free_realnums(farr,size);
  am_free_ints(iarr,size);

  return indices;
}

ivec *mk_ivec_sorted_dyv_indices(dyv *dv)
{
  ivec* indices;
  int size = dyv_size(dv);
  double* farr;

  farr = mk_farr_from_dyv(dv);

  indices = mk_ivec(size);
  indices_sort_realnums(farr,size,indices->iarr);

  am_free_realnums(farr,size);

  return indices;
}


/* Increases dv in length by 1 and shifts all elements
   with original index greater or equal to index one to the
   right and inserts val at index. */
void dyv_insert(dyv *dv,int index,double val)
{
  int i;
  add_to_dyv(dv,0);
  for ( i = dyv_size(dv)-1 ; i > index ; i-- )
    dyv_set(dv,i,dyv_ref(dv,i-1));
  dyv_set(dv,index,val);
}

void fprintf_oneline_dyv(FILE *s,const char *m1, const dyv *d, const char *m2)
{
  int i;
  fprintf(s,"%s ",m1);
  for ( i = 0 ; i < dyv_size(d) ; i++ )
    fprintf(s,"%8g%s",dyv_ref(d,i),(i==dyv_size(d)-1) ? "" : " ");
  fprintf(s,"%s",m2);
}






void indices_of_sorted_dyv(const dyv *dv,ivec *iv)
/*
   NOTE: ivec structure (integer vectors) defined in sortind.ch
   PRE: dv and iv must be same size. iv's contents irrelevant.
   POST: iv contains sorted indexes into dv, so that
         forall iv[j] is the j'th smallest element of dv.

         thus forall i,j, (i < j) => dyv_ref(dv,iv[i]) <= dyv_ref(dv,iv[j])
          and iv contains a permutation of [0 ... iv->size-1]
*/
{
  int *iarr = am_malloc_ints(dyv_size(dv));
  indices_sort_realnums(dv->farr,dv->size,iarr);
  copy_iarr_to_ivec(iarr,dv->size,iv);
  am_free_ints(iarr,dv->size);
}

ivec *mk_indices_of_sorted_dyv(const dyv *dv)
{
  ivec *iv = mk_ivec(dyv_size(dv));
  indices_of_sorted_dyv(dv,iv);
  return(iv);
}







double dyv_partial_sum( const dyv *dv, const ivec *indices)
{
  int imax, i, *iarr;
  double sum, *farr;
  imax = ivec_size(indices);
  sum = 0.0;
  iarr = indices->iarr;
  farr = dv->farr;
  for (i=0; i<imax; ++i) sum += farr[iarr[i]];
  return sum;
}

dyv *mk_dyv_from_ivec(const ivec *iv)
{
  dyv *d = mk_dyv(ivec_size(iv));
  int i;
  for ( i = 0 ; i < ivec_size(iv) ; i++ )
    dyv_set(d,i, ivec_ref(iv,i));

  return(d);
}

ivec *mk_ivec_from_dyv( const dyv *dv)
{
  int size, i;
  ivec *iv;
  size = dyv_size( dv);
  iv = mk_ivec( size);
  for (i=0; i<size; ++i) ivec_set( iv, i, (int) dyv_ref(dv,i));
  return iv;
}

/* Writes dv[a] to dv[b-1] into r_dyv.  If a or b is negative, then
   it is replaced by by len+a or len+b. Returns the number of elements
   written.  This is more-or-less equivalent to Python slice syntax.
*/
int dyv_slice( const dyv *dv, int a, int b, dyv *r_dv)
{
  int dvlen, lb, ub, r_dvlen, srcidx, dstidx;
  double dval;

  /* Set lower- and (strict) upper-bounds. */
  dvlen = dyv_size( dv);
  if (a < 0) lb = int_max(0, dvlen+a);
  else lb = int_min(dvlen, a);
  if (b < 0) ub = int_max(0, dvlen+b);
  else ub = int_min(dvlen, b);
  r_dvlen = ub - lb;

  /* Copy elements. */
  for (srcidx=lb,dstidx=0;  srcidx<ub;  ++srcidx,++dstidx) {
    dval = dyv_ref( dv, srcidx);
    dyv_set( r_dv, dstidx, dval);
  }

  return r_dvlen;
}

/* If b >= a >= 0, creates "slice" of dv on [a,b-1].  If a or b < 0, then
   it is replaced by len+a or len+b.  Indices which exceed bounds are
   replaced by appropriate indices.  Therefore this function always succeeds,
   sometimes returning an empty dyv.
*/
dyv *mk_dyv_slice( const dyv *dv, int a, int b)
{
  int dvlen, lb, ub, r_dvlen;
  dyv *r_dv;

  /* Set lower- and (strict) upper-bounds. */
  dvlen = dyv_size( dv);
  if (a < 0) lb = int_max(0, dvlen+a);
  else lb = int_min(dvlen, a);
  if (b < 0) ub = int_max(0, dvlen+b);
  else ub = int_min(dvlen, b);
  r_dvlen = ub - lb;

  r_dv = mk_dyv( r_dvlen);
  dyv_slice( dv, lb, ub, r_dv);
  return r_dv;
}