C code

AgingOMApp/C/assymheavi.c

@@ -0,0 +1,575 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <gsl/gsl_integration.h>
+
+
+double * assymheavi_pdf(double * x, int n, double k1, double k2, double d ) {
+  double * result;
+  double q;
+  
+  result = (double *) malloc(n * sizeof(double));
+  
+  for (int i = 0; i < n; ++i) {
+    q = x[i] - d;
+    if (q < 0)
+      result[i] = k2 * exp(-k1 * q) / pow(1 + exp(-k1 * q), 2);
+    else
+      result[i] = k2 * exp(- k2 * q) / pow(1 + exp(-k2 * q), 2);
+  }
+  
+  return result;
+}
+
+//for a single x avoid annoying malloc (used for integration)
+double assymheavi_pdf_singleval(double x, double k1, double k2, double d ) {
+  double result;
+  double q;
+  
+  q = x - d;
+  if (q < 0)
+    result = k2 * exp(-k1 * q) / pow(1 + exp(-k1 * q), 2);
+  else
+    result = k2 * exp(- k2 * q) / pow(1 + exp(-k2 * q), 2);
+  
+  return result;
+}
+
+double assymheavi_cdf(double a, double b, double k1, double k2, double d) {
+  double result, alpha, inv1, inv2;
+  
+  //shift distribution
+  a = a - d;
+  b = b - d;
+  
+  //handling some exceptions
+  if (a == b) {
+    return(0.0);
+  } else if (a > b) {
+    printf("Error in assymheavi_cdf: a < b.");
+    exit(1);
+  }
+  
+  //normal cases
+  alpha = k2/k1;
+  if (a < 0 && b <= 0) {
+    inv1 = 1 / (1+exp(-k1*b));
+    inv2 = 1 / (1+exp(-k1*a));
+                         
+    result = alpha * (inv1 - inv2);
+  } else if (a <= 0 && b > 0) {
+    inv1 = 1 / (1+exp(-k1*a));
+    inv2 = 1 / (1+exp(-k2*b));
+    
+    result = alpha * (0.5 - inv1) + inv2 - 0.5;
+  } else {
+    inv1 = 1 / (1+exp(-k2*b));
+    inv2 = 1 / (1+exp(-k2*a));
+    
+    result = inv1 - inv2;
+  }
+  
+  return result;
+}
+
+double * assymheavi_cdf_vec(double * x, int n, double k1, double k2, double d) {
+  double * result;
+  double a;
+  
+  if (n == 1) {
+    result = (double *) malloc(sizeof(double));
+    result[0] = 0.0;
+    return(result);
+  } else {
+    a = x[0];
+    result = (double *) malloc( (n - 1) * sizeof(double));
+    
+    for (int i = 1; i < n; ++i) {
+      result[i-1]  = assymheavi_cdf(a, x[i], k1, k2, d);
+    }
+    
+    return(result);
+  }
+}
+
+double assymheavi_intcdf(double a, double b, double k1, double k2, double d)  {
+  double result, alpha;
+  
+  //shift distribution
+  a = a - d;
+  b = b - d;
+  
+  //some exceptions handling weird input
+  if (a == b) {
+    return(0.0);
+  } else if (a > b) {
+    printf("Error in assymheavi_intcdf: a < b.");
+    exit(1);
+  }
+  
+  alpha = k2/k1;
+  
+  if (a < 0 && b <= 0) {
+    double ea, eb, C1;
+    
+    ea = exp(k1 * a);
+    eb = exp(k1 * b);
+    
+    C1 = -alpha * (log(ea+1)/k1 - (a * ea)/(ea +1));
+    
+    result = alpha * (log(eb +1)/k1 - (b*ea)/(ea+1)) + C1;
+  } else if (a < 0 && b >= 0) {
+    double e1a, e1ai, e2a, C1, C2, end_left, begin_right;
+    const double log2 = 0.69314718056;
+    
+    //some recurrent expressions
+    e1a = exp(k1 * a);
+    e1ai = 1/e1a;
+    e2a = exp(k2 * a);
+    
+    //integral constants
+    C1 = -alpha * (log(e1a +1)/k1 - (a * e1a)/(e1a+1));
+    end_left = alpha * (log2/k1) + C1;
+    begin_right = -(alpha * (0.5 - 1/(1+e1ai)) - 0.5) * a +
+      (1/k2) * log(2/(e2a+1));
+    C2 = end_left - begin_right;
+    
+    result = (alpha * (0.5 - 1/(1+e1ai)) - 0.5) * (b-a) +
+      (1/k2) * log((exp(k2*b)+1)/(e2a+1)) + C2;
+  } else {
+    double ea, eb, C3;
+    
+    ea = exp(k2 * a);
+    eb = exp(k2 * b);
+    
+    C3 = - (log(ea + 1)/k2 - (a * ea)/(ea + 1));
+    
+    result = log(eb + 1)/k2 - (b * ea)/(ea + 1) + C3;
+  }
+  
+  return result;
+}
+
+double * assymheavi_intcdf_vec(double * x, int n, double k1, double k2, double d) {
+  double * result;
+  double a;
+  
+  if (n == 1) {
+    result = (double *) malloc(sizeof(double));
+    result[0] = 0.0;
+    return(result);
+  } else {
+    a = x[0];
+    result = (double *) malloc( (n - 1) * sizeof(double));
+    
+    for (int i = 1; i < n; ++i) {
+      result[i-1]  = assymheavi_intcdf(a, x[i], k1, k2, d);
+    }
+    
+    return(result);
+  }
+}
+
+double assymheavi_Cmu(double a, double b, double k1, double k2, double d) {
+  double p1, p2;
+  p1 = assymheavi_cdf(a, b, k1, k2, d) * b;
+  p2 = - assymheavi_intcdf(a, b, k1, k2, d);
+  
+  return (p1 + p2);
+}
+
+double * find_bounds(double k1, double k2, double d, double rel_threshold, 
+                     double min, double max) {
+  const int nsteps = 1000;
+  double step, maxval;
+  int maxindex;
+  double x[nsteps];
+  double * f, * bounds;
+  
+  int success_flag;
+  
+  
+  step = (max - min) / (nsteps - 1);
+  
+  //define x sequence
+  for (int i = 0; i < nsteps; ++i) {
+    x[i] = min + step * i;
+  }
+  
+  //calculate the pdf values
+  f = assymheavi_pdf(x, nsteps, k1, k2, d);
+  
+  
+  //determine maximum value
+  maxval = f[0];
+  maxindex = 0;
+  for (int i = 1; i < nsteps; ++i) {
+    if (f[i] > maxval) {
+      maxval = f[i];
+      maxindex = i;
+    }
+  }
+  
+  bounds = (double *) malloc(2 * sizeof(double));
+  
+  //determine the left bound
+  success_flag = 0;
+  for (int i = maxindex-1; i > -1; --i) {
+    if ( (f[i]/maxval) < rel_threshold) {
+      bounds[0] = x[i];
+      success_flag = 1;
+      break;
+    }
+  }
+  if (success_flag == 0)
+    bounds[0] = min;
+  
+  //determine the right bound
+  success_flag = 0;
+  for (int i = maxindex+1; i < nsteps; ++i) {
+    if ( (f[i]/maxval) < rel_threshold) {
+      bounds[1] = x[i];
+      success_flag = 1;
+      break;
+    }
+  }
+  if (success_flag == 0)
+    bounds[1] = max;
+  
+  //free memory allocated in assymheavi_pdf
+  free(f);
+  
+  return bounds;
+}
+
+
+double mom_func(double x, void * parm) {
+  double * p;
+  double k1, k2, d, center, power, pdf;
+  
+  //read out parms
+  p = (double *) parm;
+  k1 = p[0];
+  k2 = p[1];
+  d = p[2];
+  center = p[3];
+  power = p[4];
+  
+  //compute function to be integrated
+  pdf = assymheavi_pdf_singleval(x, k1, k2, d);
+  return(pdf * pow( (x-center), power));
+}
+
+//calculate concentration, mean, and centralized higher moments
+void determine_moments(double * k1, double * k2, double * d, 
+                       double * xmin, double * xmax, 
+                       int * use_noncentral_moment, double * center,
+                       double * moments) {
+
+  double * bounds;
+  
+  //determine bounds for integrals
+  bounds = find_bounds(*k1, *k2, *d, 1e-9, *xmin, *xmax);
+  
+  //determine the concentration
+  moments[0] = assymheavi_cdf(*xmin, *xmax, *k1, *k2, *d);
+  
+  //determine the mean (this is not a central moment, as then it would be 0)
+  moments[1] = assymheavi_Cmu(*xmin, *xmax, *k1, *k2, *d) / moments[0];
+  
+  //determine the variance and skewness
+  for (int i_moment = 2; i_moment < 4; ++i_moment) {
+    gsl_integration_workspace * w;
+    gsl_function F;
+    double result, error;
+    
+    double * parm;
+    parm = (double *) malloc(sizeof(double)*5);
+    parm[0] = *k1;
+    parm[1] = *k2;
+    parm[2] = *d;
+
+    if (use_noncentral_moment == 0)
+      parm[3] = *center;
+    else
+      parm[3] = moments[1];
+    
+    parm[4] = (double) i_moment;
+    
+    w = gsl_integration_workspace_alloc (1000);
+    F.function = &mom_func;
+    F.params = parm;
+    
+    gsl_integration_qag(&F, //const gsl_function *f
+                        bounds[0], //begin integration interval
+                        bounds[1], //end integration interval
+                        0, //double epsabs
+                        1e-6, //double epsrel
+                        1000, //size_t limit (see allocation w)
+                        GSL_INTEG_GAUSS15, //int key
+                        w, //gsl_integration_workspace *workspace
+                        &result, //double *result
+                        &error); //double *abserr
+    
+    moments[i_moment] = result / moments[0];
+    
+    gsl_integration_workspace_free(w);
+    free(parm);
+  }
+  
+  free(bounds);
+}
+
+
+//central_moms [0] is mean, [1] is central var, [2] is central skew 
+void transform_moments(double * central_moms, 
+                         double * center, double * new_moms) {
+  
+  double diff_center;
+  
+  //mu0 is always 1
+  
+  //mu1
+  diff_center = (central_moms[0] - center[0]);
+  new_moms[0] = diff_center;
+  
+  //mu2
+  new_moms[1] = pow(diff_center, 2) + central_moms[1];
+  
+  //mu3
+  new_moms[2] = pow(diff_center, 3) + 3 * central_moms[1] * diff_center +
+    central_moms[2];
+}
+
+//k1, k2, d: are parms to be tested
+//xmin and xmax: are the domain of the pdf
+//use_noncentral_moms: boolean to indicate if transformations are requested
+//new_center: vector with centers for transformations
+//ncenter: the length of vector new_center
+//include_conc: boolean indication if moment vector contains the concentration
+//moms: the moments vector, starting from mu0 = concentration
+//total_error: the result, length must be same as ncenter
+void cost_function(double * k1, double * k2, double * d,
+                     double * xmin, double * xmax,
+                     int * use_noncentral_moms, double * new_center, int * ncenter,
+                     int * include_conc, double * moms, double * total_error) {
+  
+  double * calc_moms, * trans_moms, * which_moms; 
+  
+  if (*use_noncentral_moms == 1) {
+    trans_moms = (double *) malloc(sizeof(double) * 4);
+    trans_moms[0] = moms[0];
+  }
+  
+  calc_moms = (double *) malloc(sizeof(double) * 4);
+  
+  for (int i_center = 0; i_center < *ncenter; ++i_center) {
+    //calc_moms: [0] conc, [1] mean, [2] var, [3]skew
+    determine_moments(k1, k2, d, 
+                      xmin, xmax, 
+                      use_noncentral_moms, &new_center[i_center],
+                      calc_moms);
+    
+    if (*use_noncentral_moms == 1) {
+      //transform given moments to be fitted
+      transform_moments(&moms[1], &new_center[i_center], &trans_moms[1]);
+      //direct pointer for comparison to these transformed given moments
+      which_moms = trans_moms;
+    } else
+      which_moms = moms; //use standard given moments
+    
+    //calculate error of the moments
+    double scale_mom1, scale_mom2, scale_mom3;//, stdev;
+    //stdev = sqrt(which_moms[2]);
+    scale_mom1 = sqrt(which_moms[1]);
+    scale_mom2 = which_moms[2];
+    scale_mom3 = pow(which_moms[3], 1.5);
+    
+    if (*include_conc == 1) {
+      total_error[i_center] = fabs(calc_moms[0] - which_moms[0]) +
+        fabs((calc_moms[1] - which_moms[1]) / scale_mom1) +
+        fabs((calc_moms[2] - which_moms[2]) / scale_mom2) +
+        fabs((calc_moms[3] - which_moms[3]) / scale_mom3);
+    } else {
+      total_error[i_center] = 
+        fabs((calc_moms[1] - which_moms[1]) / scale_mom1) +
+        fabs((calc_moms[2] - which_moms[2]) / scale_mom2) +
+        fabs((calc_moms[3] - which_moms[3]) / scale_mom3);
+    }
+  }
+  
+  //free up memory
+  free(calc_moms);
+  if (*use_noncentral_moms == 1)
+    free(trans_moms);
+}
+
+double cost_function_sum(double * k1, double * k2, double * d,
+                         double * xmin, double * xmax,
+                         int * use_noncentral_moms, double * new_center, 
+                         int * ncenter, int * include_conc, double * moms) {
+  
+  double * total_error;
+  total_error = (double *) malloc(sizeof(double) * ncenter[0]);
+
+  
+  cost_function(k1, k2, d,
+                xmin, xmax,
+                use_noncentral_moms, new_center, ncenter,
+                include_conc, moms, total_error);
+
+  double sum_error;
+  if (* ncenter == 1) {
+    sum_error = total_error[0];
+  } else {
+    sum_error = 0;
+    for (int i = 0; i < * ncenter; ++i) {
+       sum_error = sum_error + total_error[i];
+    }
+  }
+  
+  free(total_error);
+    
+  return(sum_error);
+}
+
+
+void multiopt_boxed(double * start, double * xmin, double * xmax,
+                    double * lower, double * upper, 
+                    double * moms, int * use_noncentral_moms, double * center, 
+                    int * ncenter, int * niter, int * include_conc, 
+                    double * result) {
+  double k1, k2, d, mult, perturb, f_old;
+  int npar;
+  
+  npar = 3 + include_conc[0];
+  
+  for (int i = 0; i < npar; ++i) {
+    result[i] = start[i];
+  }
+  
+  perturb = 0.05;
+  
+  f_old = cost_function_sum(&result[0], &result[1], &result[2],
+                            xmin, xmax,
+                            use_noncentral_moms, center, ncenter,
+                            include_conc, moms);
+
+  for (int i = 0; i < niter[0]; ++i) {
+    perturb = perturb * 0.5;
+    for(int j_par = 0; j_par < npar; ++j_par)  {
+
+      int max_search_iter;
+      max_search_iter = ceil( (upper[j_par] - lower[j_par])/perturb);
+      
+      for (int k_search = 0; k_search < max_search_iter; ++k_search) {
+        double f_pert_plus, f_pert_min;
+        double errorsdx_plus, errorsdx_min;
+        
+        result[j_par] = result[j_par] + perturb;
+
+        f_pert_plus = cost_function_sum(&result[0], &result[1], &result[2],
+                                        xmin, xmax,
+                                        use_noncentral_moms, center, ncenter,
+                                        include_conc, moms);
+        
+        result[j_par] = result[j_par] - 2 * perturb;
+        
+        f_pert_min = cost_function_sum(&result[0], &result[1], &result[2],
+                                        xmin, xmax,
+                                        use_noncentral_moms, center, ncenter,
+                                        include_conc, moms);
+        //reset to original value
+        result[j_par] = result[j_par] + perturb;
+        
+        //errors up and down
+        errorsdx_plus = fabs(f_pert_plus) - fabs(f_old);  //scalar
+        errorsdx_min = fabs(f_pert_min) - fabs(f_old);  //scalar
+        
+        if (errorsdx_plus < 0 && ((result[j_par]+perturb) < upper[j_par]) ) { 
+          //error becomes smaller for greater x
+          result[j_par] = (result[j_par]+perturb);
+          f_old = f_pert_plus;
+        } else if (errorsdx_min < 0 && ((result[j_par]-perturb) > lower[j_par])) {
+          result[j_par] = result[j_par] - perturb;
+          f_old = f_pert_min;
+        } else {
+          break; //exit inner loop
+        }
+      }
+    }
+  }
+}
+
+/*
+	Main function is for testing purposes
+*/
+int main() {
+  double k1, k2, d;
+  double * res;
+  
+  double x[3] = {22, 24, 26};
+  
+  k1 = 4;
+  k2 = 4;
+  d = 24.7;
+  
+  res = assymheavi_pdf(x, 3, k1, k2, d );
+  
+  printf("The results are: \n");
+  for (int i = 0; i < 3; ++i) {
+    printf("x = %e gives %e\n", x[i], res[i]);
+  }
+  
+  printf("Test Cmu %e\n",
+         assymheavi_Cmu(25, 26, k1, k2, d)
+         );
+  
+  double * bounds;
+  
+  bounds = find_bounds(k1, k2, d, 1e-8, 0, 100);
+  
+  printf("Left bound: %f, right bound %f\n", bounds[0], bounds[1]);
+  
+  
+  // test calculation of moments 
+  double xmin, xmax;
+  double * moments, * new_moms;
+  int use_noncentral_mom;
+  
+  k1 = 3.9075678;
+  k2 = 0.9210129;
+  d = 1.2927411;
+  xmin = 0;
+  xmax = 100;
+  moments = (double *) malloc(sizeof(double) * 4);
+  
+  use_noncentral_mom = 0;
+  
+  
+  determine_moments(&k1, &k2, &d, 
+                    &xmin, &xmax, 
+                    &use_noncentral_mom, NULL,
+                    moments);
+  
+  printf("conc %e ... mean %e ... var %e ... skew %e \n",
+         moments[0], moments[1], moments[2], moments[3]);
+  
+  printf("Transform moments to center = 13\n");
+  new_moms = (double *) malloc(sizeof(double)*3);
+  double center = 13.0;
+  
+  transform_moments(&moments[1], &center, new_moms);
+  
+  printf("The moments centered around 13 are: %e ... %e ... %e \n",
+         new_moms[0], new_moms[1], new_moms[2]);
+  
+  free(bounds);
+  free(res);
+  free(moments);
+  free(new_moms);
+
+  
+  return 0;
+}
+

AgingOMApp/C/comp.sh

@@ -0,0 +1,6 @@
+filename="assymheavi"
+
+rm ${filename}.o
+rm ${filename}.so
+#Using R for compilation. GNU scientific library must have been installed.
+R CMD SHLIB ${filename}.c -lm -lgsl