gs_pm_model.c

Go to the documentation of this file.

#include <stdlib.h>
#include <stdio.h>
#include "gs_pm_model.h"
#include <sys/types.h>
#include <sys/uio.h>
#include <unistd.h>
#include <string.h>

int
gs_pm_save(gs_pm_model_t *md, int fd)
{
  int i;
  ssize_t res;
  res = write(fd, &(md->nb_params), sizeof(md->nb_params));
  if(res == -1)
    return -1;
  res = write(fd, &(md->nb_categories), sizeof(md->nb_categories));
  if(res == -1)
    return -1;
  res = write(fd, &(md->max_run), sizeof(md->max_run));
  if(res == -1)
    return -1;
  res = write(fd, &(md->run_index), sizeof(md->run_index));
  if(res == -1)
    return -1;
  res = write(fd, &(md->full), sizeof(md->full));
  if(res == -1)
    return -1;
  res = write(fd, md->params, sizeof(md->params[0]) * md->nb_params);
  if(res == -1)
    return -1;
  res =
      write(fd, md->categories, sizeof(md->categories[0]) * md->nb_categories);
  if(res == -1)
    return -1;
  for(i = 0; i < md->max_run; i++) {
    res = write(fd, md->param_mat[i],
           sizeof(md->param_mat[i][0]) * (md->nb_params + md->nb_categories));
    if(res == -1)
      return -1;
  }
  res = write(fd, md->time_vec, sizeof(md->time_vec[0]) * md->max_run);
  if(res == -1)
    return -1;
  return 0;
}

int
gs_pm_dump_model_terse(gs_pm_model_t *md, FILE *pf)
{
  int i;

  fprintf(pf, "Num params:      %d\n", md->nb_params);
  fprintf(pf, "Num categories:  %d\n", md->nb_categories);
  fprintf(pf, "Max runs:        %d\n", md->max_run);
  fprintf(pf, "Run index:       %d\n", md->run_index);
  fprintf(pf, "Is full:         %d\n", md->full);
  fprintf(pf, "Model Ok:        %d\n", md->model_ok);
  fprintf(pf, "Params: ");
  for(i = 0; i < md->nb_params; i++)
    fprintf(pf, "%g ", md->params[i]);
  fprintf(pf, "\n");
  fprintf(pf, "Categories: ");
  for(i = 0; i < md->nb_categories; i++)
    fprintf(pf, "%.3g ", md->categories[i]);
  fprintf(pf, "\n");
  fprintf(pf, "Coefficients: ");
  for(i = 0; i < md->nb_params; i++)
    fprintf(pf, "%.3g ", md->coef_vec[i]);
  fprintf(pf, "\n");

  return 0;
}

int
gs_pm_dump_model(gs_pm_model_t *md, FILE *pf)
{
  int i, j;

  gs_pm_dump_model_terse(md, pf);

  for(i = 0; i < md->max_run; i++) {
    for(j = 0; j < md->nb_params + md->nb_categories; j++)
      fprintf(pf, "%.52e ", md->param_mat[i][j]);
    fprintf(pf, "\n");
  }
  for(i = 0; i < md->max_run; i++)
    fprintf(pf, "%.52e ", md->time_vec[i]);
  fprintf(pf, "\n");

  return 0;
}

int
gs_pm_save_model(gs_pm_model_t *md, char *filename)
{
  FILE *pf;
  int i, j;
  pf = fopen(filename, "w");
  if(!pf)
    return -1;

  fprintf(pf, "%d\n", md->nb_params);
  fprintf(pf, "%d\n", md->nb_categories);
  fprintf(pf, "%d\n", md->max_run);
  fprintf(pf, "%d\n", md->run_index);
  fprintf(pf, "%d\n", md->full);
  for(i = 0; i < md->nb_params; i++)
    fprintf(pf, "%.52e ", md->params[i]);
  fprintf(pf, "\n");
  for(i = 0; i < md->nb_categories; i++)
    fprintf(pf, "%.52e ", md->categories[i]);
  fprintf(pf, "\n");
  for(i = 0; i < md->max_run; i++) {
    for(j = 0; j < md->nb_params + md->nb_categories; j++)
      fprintf(pf, "%.52e ", md->param_mat[i][j]);
    fprintf(pf, "\n");
  }
  for(i = 0; i < md->max_run; i++)
    fprintf(pf, "%.52e ", md->time_vec[i]);
  fprintf(pf, "\n");
  fclose(pf);
  return 0;
}

gs_pm_model_t *
gs_pm_load(int fd)
{
  int i, res;
  gs_pm_model_t *md, dub;
  res = read(fd, &(dub.nb_params), sizeof(dub.nb_params));
  if(res == -1)
    return NULL;
  res = read(fd, &(dub.nb_categories), sizeof(dub.nb_categories));
  if(res == -1)
    return NULL;
  res = read(fd, &(dub.max_run), sizeof(dub.max_run));
  if(res == -1)
    return NULL;
  md = gs_pm_init_model(dub.nb_categories, dub.nb_params, dub.max_run);
  res = read(fd, &(md->run_index), sizeof(md->run_index));
  if(res == -1) {
    gs_pm_free_model(md);
    return NULL;
  }
  res = read(fd, &(md->full), sizeof(md->full));
  if(res == -1) {
    gs_pm_free_model(md);
    return NULL;
  }
  res = read(fd, md->params, sizeof(md->params[0]) * md->nb_params);
  if(res == -1) {
    gs_pm_free_model(md);
    return NULL;
  }
  res = read(fd, md->categories, sizeof(md->categories[0]) * md->nb_categories);
  if(res == -1) {
    gs_pm_free_model(md);
    return NULL;
  }
  for(i = 0; i < md->max_run; i++) {
    res =
        read(fd, md->param_mat[i],
             sizeof(md->param_mat[i][0]) * (md->nb_params + md->nb_categories));
    if(res == -1) {
      gs_pm_free_model(md);
      return NULL;
    }
  }
  res = read(fd, md->time_vec, sizeof(md->time_vec[0]) * md->max_run);
  if(res == -1) {
    gs_pm_free_model(md);
    return NULL;
  }
  md->model_ok = gs_pm_compute_model(md);
  return md;
}

gs_pm_model_t *
gs_pm_load_model(char *filename)
{
  FILE *pf;
  int i, j;
  gs_pm_model_t *md;
  int nb_params, nb_categories, max_run;
  pf = fopen(filename, "r");
  if(!pf)
    return NULL;
  fscanf(pf, "%d\n", &nb_params);
  fscanf(pf, "%d\n", &nb_categories);
  fscanf(pf, "%d\n", &max_run);
  md = gs_pm_init_model(nb_categories, nb_params, max_run);
  fscanf(pf, "%d\n", &(md->run_index));
  fscanf(pf, "%d\n", &(md->full));
  for(i = 0; i < md->nb_params; i++)
    fscanf(pf, "%le ", &(md->params[i]));
  fscanf(pf, "\n");
  for(i = 0; i < md->nb_categories; i++)
    fscanf(pf, "%le ", &(md->categories[i]));
  fscanf(pf, "\n");
  for(i = 0; i < md->max_run; i++) {
    for(j = 0; j < md->nb_params + md->nb_categories; j++)
      fscanf(pf, "%le ", &(md->param_mat[i][j]));
    fscanf(pf, "\n");
  }
  for(i = 0; i < md->max_run; i++)
    fscanf(pf, "%le ", &(md->time_vec[i]));
  fscanf(pf, "\n");
  md->model_ok = gs_pm_compute_model(md);
  fclose(pf);
  return md;
}

void
gs_pm_free_model(gs_pm_model_t *md)
{
  int i;
  if(md->nb_params)
    free(md->params);
  if(md->nb_categories)
    free(md->categories);
  if(md->max_run) {
    free(md->time_vec);
    for(i = 0; i < md->max_run; i++)
      if(md->nb_categories + md->nb_params)
        free(md->param_mat[i]);
    free(md->param_mat);
  }
  if(md->coef_vec)
    free(md->coef_vec);
  free(md);
}

gs_pm_model_t *
gs_pm_init_model(int nb_categories, int nb_params, int nb_lines)
{
  int i, j;
  gs_pm_model_t *res;
  res = (gs_pm_model_t *) malloc(sizeof(gs_pm_model_t));
  res->params = (double *) calloc(nb_params, sizeof(double));
  res->categories = (double *) calloc(nb_categories, sizeof(double));
  res->nb_categories = nb_categories;
  res->nb_params = nb_params;
  res->run_index = 0;
  res->max_run = nb_lines;
  res->param_mat = (double **) malloc(sizeof(double *) * nb_lines);
  for(i = 0; i < nb_lines; i++) {
    res->param_mat[i] =
        (double *) calloc(nb_categories + nb_params, sizeof(double));
    for(j = 0; j < nb_categories; j++)
      res->param_mat[i][j] = -1;
  }
  res->time_vec = (double *) calloc(nb_lines, sizeof(double));
  res->coef_vec = (double *) calloc(nb_lines, sizeof(double));
  res->full = 0;
  res->model_ok = 0;
  return res;
}

void
gs_pm_store_timing(double time, gs_pm_model_t *md)
{
  int i;
  int j;
  i = md->run_index;
  for(j = 0; j < md->nb_categories; j++)
    md->param_mat[i][j] = md->categories[j];
  for(j = 0; j < md->nb_params; j++)
    md->param_mat[i][j + md->nb_categories] = md->params[j];
  md->time_vec[i] = time;
  if(md->run_index == md->max_run - 1) {
    md->full = 1;
    md->run_index = 0;
  }
  else
    md->run_index++;

  md->model_ok = gs_pm_compute_model(md);
}

gs_pm_model_t *gs_pm_clone_model(gs_pm_model_t *md1){
  gs_pm_model_t *md;
  int i;

  md=gs_pm_init_model(md1->nb_categories,md1->nb_params,md1->max_run);
  md->run_index=md1->run_index;
  md->full=md1->full;;
  memcpy(md->params,md1->params,md1->nb_params*sizeof(double));
  memcpy(md->categories,md1->categories,md1->nb_categories*sizeof(double));
  for(i=0;i<md->max_run;i++)
    memcpy(md->param_mat[i],md1->param_mat[i],sizeof(double)*(md1->nb_params+md1->nb_categories));
  memcpy(md->time_vec,md1->time_vec,md1->max_run*sizeof(double));
  md->model_ok=gs_pm_compute_model(md);
  return md;
}

/* in :a model 
  out : table of categories and coefficients
  categories = NULL if no categories.
  Returns: number of lines of each table, 0 if no models*/
int gs_pm_all_models(gs_pm_model_t *md,double ***categories, double ***coefficients){
  gs_pm_model_t *md_copy;
  int n,res,i,j;
  double *coef;

  if(md->nb_categories==0){
    *categories=NULL;
    if(gs_pm_compute_model(md)==0){
      coef=NULL;
      res=0;
    }else{
      res=1;
      coef=(double*)malloc(md->nb_params*sizeof(double));
      memcpy(coef,md->coef_vec,md->nb_params*sizeof(double));
    }
    *coefficients=(double**)malloc(sizeof(double*));
    *coefficients[0]=coef;
    return res;
  }

  
  md_copy=gs_pm_clone_model(md);
  

  n=md_copy->full?md_copy->max_run:md_copy->run_index-1;
  

  /* compute res the number of lines of categories and coefficients */
  i=0;
  res=0;
  do{
    memcpy(md->categories,md_copy->param_mat[i],md_copy->nb_categories*sizeof(double));
    if(gs_pm_compute_model(md)){
       res++;
    }
    for(j=i;j<n;j++){
      int k,ok=1;
      /*find lines of the matrix that correspond to the current run*/
      for(k=0;k<md->nb_categories;k++){
        if(md_copy->param_mat[j][k]!=md->categories[k]){
          ok=0;
          break;
        }
      }
      if(ok){
        md_copy->time_vec[j]=-1;
      }
    }
    while((md_copy->time_vec[i]==-1)&&(i<n))
      i++;
  }while(i!=n);


  
  memcpy(md_copy->time_vec,md->time_vec,n*sizeof(double));

  *coefficients=(double**)malloc(sizeof(double*)*res);
  *categories=(double**)malloc(sizeof(double*)*res);

  

  i=0;
  res=0;
  do{
    memcpy(md->categories,md_copy->param_mat[i],md_copy->nb_categories*sizeof(double));
    if(gs_pm_compute_model(md)){
      (*coefficients)[res]=(double*)malloc(sizeof(double)*md->nb_params);
      (*categories)[res]=(double*)malloc(sizeof(double)*md->nb_categories);
      memcpy((*coefficients)[res],md->coef_vec,md_copy->nb_params*sizeof(double));
      memcpy((*categories)[res],md->categories,md_copy->nb_categories*sizeof(double));
      res++;
    }
    for(j=i;j<n;j++){
      int k,ok=1;
      /*find lines of the matrix that correspond to the current run*/
      for(k=0;k<md->nb_categories;k++){
        if(md_copy->param_mat[j][k]!=md->categories[k]){
          ok=0;
          break;
        }
      }
      if(ok){
        md_copy->time_vec[j]=-1;
      }
    }
    while((md_copy->time_vec[i]==-1)&&(i<n))
      i++;
  }while(i!=n);

  gs_pm_free_model(md_copy);

  return res;
}

    /*
     * retruns 0 if no model 1 if able to compute a model
     */
int
gs_pm_compute_model(gs_pm_model_t *md)
{
  double *mat, *vec;
  int i, j, M, N, info, l, n;
  N = md->nb_params;
  mat = (double *) calloc(md->max_run * N, sizeof(double));
  vec = md->coef_vec;
  n = md->full ? md->max_run : md->run_index;
  for(i = 0, l = 0; i < n; i++) {
    int k, ok = 1;

    /*
     * find lines of the matrix that correspond to the current run 
     */
    for(k = 0; k < md->nb_categories; k++) {
      if(md->param_mat[i][k] != md->categories[k]) {
        ok = 0;
        break;
      }
    }

    /*
     * if ok, add this line to the input matrix of the solver 
     */
    if(ok) {
      vec[l] = md->time_vec[i];
      for(j = 0; j < N; j++)
        mat[j*md->max_run+l]=md->param_mat[i][j+md->nb_categories];
      l++;
    }
  }
  M = l;

  /*
   * printf("M=%d, N=%d n=%d\n",M,N,n);
   * 
   * for(i=0;i<M;i++){ for(j=0;j<N;j++) printf("%.2e
   * ",mat[j*md->max_run+i]); printf("\n"); }
   * 
   * for(i=0;i<M;i++) printf("%.2e ",vec[i]);
   * 
   * printf("\n"); 
   */
  if(M < 1) {
    free(mat);
    return 0;
  }

#ifdef DGELS
  {
    double *work;
    int nrhs = 1, lwork = 1;
    if((!(md->full)) && (M < md->nb_params + md->nb_categories)) {
      free(mat);
      return 0;
    }

    /*
     * Workspace query 
     */
    work = (double *) malloc(N * sizeof(double));
    if(work == NULL) {
      printf("error of memory allocation\n");
      exit(0);
    }
    lwork = -1;
    dgels_("N", &M, &N, &nrhs, mat, &(md->max_run), vec, &(md->max_run), work,
           &lwork, &info);
    lwork = (int) work[0];
    free(work);

    /*
     * allocation of the workspace 
     */
    work = (double *) malloc(lwork * sizeof(double));
    if(work == NULL) {
      printf("error of memory allocation\n");
      exit(0);
    }
    dgels_("N", &M, &N, &nrhs, mat, &(md->max_run), vec, &(md->max_run),
           work, &lwork, &info);
    if(info != 0) {
      printf("\nError on parameter %d!\n", -info);
      exit(0);
    }
    free(work);
    free(mat);
    return 1;
  }

#else                           /* */
  {
    double *x, rnorm, *w, *zz;
    int *index;
    x = (double *) malloc(N * sizeof(double));
    if(x == NULL) {
      printf("error of memory allocation\n");
      exit(0);
    }
    w = (double *) malloc(N * sizeof(double));
    if(w == NULL) {
      printf("error of memory allocation\n");
      exit(0);
    }
    zz = (double *) malloc(M * sizeof(double));
    if(zz == NULL) {
      printf("error of memory allocation\n");
      exit(0);
    }
    index = (int *) malloc(N * sizeof(int));
    if(index == NULL) {
      printf("error of memory allocation\n");
      exit(0);
    }
    nnls_c(mat, &(md->max_run), &M, &N, vec, x, &rnorm, w, zz, index, &info);

    /*
     * printf("info : %d, norm : %f\n",info,rnorm); 
     */

    free(w);
    free(zz);
    free(index);
    free(mat);
    for(i = 0; i < N; i++)
      vec[i] = x[i];
    free(x);
    return 1;
  }

#endif                          /* */
}

double
gs_pm_guess_time(gs_pm_model_t *md)
{
  double *vec;
  double time = 0;
  int i;
  if(!md->model_ok)
    return -1;
  vec = md->coef_vec;

  /*
   * printf("\n"); for(i=0;i<nb_params;i++) printf("%e*%.0f
   * ",vec[i],param_vec[i]); printf("\n"); 
   */
  for(i = 0; i < md->nb_params; i++) {
    time+=vec[i]*md->params[i];

    /*
     * printf("%e*%f ",vec[i],param_vec[i]); 
     */
  }

  /*
   * printf("\n"); 
   */
  return time;
}

void
gs_pm_display_model(gs_pm_model_t *md)
{
  double *vec;
  int i;
  if(md == NULL)
    return;
  if(!md->model_ok)
    return;
  vec = md->coef_vec;
  for(i = 0; i < md->nb_params; i++)
    printf("%.1e ",vec[i]);
  printf("\n");
}