gs_pm_model.c

Go to the documentation of this file.
00001 #include <stdlib.h>
00002 #include <stdio.h>
00003 #include "gs_pm_model.h"
00004 #include <sys/types.h>
00005 #include <sys/uio.h>
00006 #include <unistd.h>
00007 #include <string.h>
00008 
00009 int
00010 gs_pm_save(gs_pm_model_t *md, int fd)
00011 {
00012   int i;
00013   ssize_t res;
00014   res = write(fd, &(md->nb_params), sizeof(md->nb_params));
00015   if(res == -1)
00016     return -1;
00017   res = write(fd, &(md->nb_categories), sizeof(md->nb_categories));
00018   if(res == -1)
00019     return -1;
00020   res = write(fd, &(md->max_run), sizeof(md->max_run));
00021   if(res == -1)
00022     return -1;
00023   res = write(fd, &(md->run_index), sizeof(md->run_index));
00024   if(res == -1)
00025     return -1;
00026   res = write(fd, &(md->full), sizeof(md->full));
00027   if(res == -1)
00028     return -1;
00029   res = write(fd, md->params, sizeof(md->params[0]) * md->nb_params);
00030   if(res == -1)
00031     return -1;
00032   res =
00033       write(fd, md->categories, sizeof(md->categories[0]) * md->nb_categories);
00034   if(res == -1)
00035     return -1;
00036   for(i = 0; i < md->max_run; i++) {
00037     res = write(fd, md->param_mat[i],
00038            sizeof(md->param_mat[i][0]) * (md->nb_params + md->nb_categories));
00039     if(res == -1)
00040       return -1;
00041   }
00042   res = write(fd, md->time_vec, sizeof(md->time_vec[0]) * md->max_run);
00043   if(res == -1)
00044     return -1;
00045   return 0;
00046 }
00047 
00048 int
00049 gs_pm_dump_model_terse(gs_pm_model_t *md, FILE *pf)
00050 {
00051   int i;
00052 
00053   fprintf(pf, "Num params:      %d\n", md->nb_params);
00054   fprintf(pf, "Num categories:  %d\n", md->nb_categories);
00055   fprintf(pf, "Max runs:        %d\n", md->max_run);
00056   fprintf(pf, "Run index:       %d\n", md->run_index);
00057   fprintf(pf, "Is full:         %d\n", md->full);
00058   fprintf(pf, "Model Ok:        %d\n", md->model_ok);
00059   fprintf(pf, "Params: ");
00060   for(i = 0; i < md->nb_params; i++)
00061     fprintf(pf, "%g ", md->params[i]);
00062   fprintf(pf, "\n");
00063   fprintf(pf, "Categories: ");
00064   for(i = 0; i < md->nb_categories; i++)
00065     fprintf(pf, "%.3g ", md->categories[i]);
00066   fprintf(pf, "\n");
00067   fprintf(pf, "Coefficients: ");
00068   for(i = 0; i < md->nb_params; i++)
00069     fprintf(pf, "%.3g ", md->coef_vec[i]);
00070   fprintf(pf, "\n");
00071 
00072   return 0;
00073 }
00074 
00075 int
00076 gs_pm_dump_model(gs_pm_model_t *md, FILE *pf)
00077 {
00078   int i, j;
00079 
00080   gs_pm_dump_model_terse(md, pf);
00081 
00082   for(i = 0; i < md->max_run; i++) {
00083     for(j = 0; j < md->nb_params + md->nb_categories; j++)
00084       fprintf(pf, "%.52e ", md->param_mat[i][j]);
00085     fprintf(pf, "\n");
00086   }
00087   for(i = 0; i < md->max_run; i++)
00088     fprintf(pf, "%.52e ", md->time_vec[i]);
00089   fprintf(pf, "\n");
00090 
00091   return 0;
00092 }
00093 
00094 int
00095 gs_pm_save_model(gs_pm_model_t *md, char *filename)
00096 {
00097   FILE *pf;
00098   int i, j;
00099   pf = fopen(filename, "w");
00100   if(!pf)
00101     return -1;
00102 
00103   fprintf(pf, "%d\n", md->nb_params);
00104   fprintf(pf, "%d\n", md->nb_categories);
00105   fprintf(pf, "%d\n", md->max_run);
00106   fprintf(pf, "%d\n", md->run_index);
00107   fprintf(pf, "%d\n", md->full);
00108   for(i = 0; i < md->nb_params; i++)
00109     fprintf(pf, "%.52e ", md->params[i]);
00110   fprintf(pf, "\n");
00111   for(i = 0; i < md->nb_categories; i++)
00112     fprintf(pf, "%.52e ", md->categories[i]);
00113   fprintf(pf, "\n");
00114   for(i = 0; i < md->max_run; i++) {
00115     for(j = 0; j < md->nb_params + md->nb_categories; j++)
00116       fprintf(pf, "%.52e ", md->param_mat[i][j]);
00117     fprintf(pf, "\n");
00118   }
00119   for(i = 0; i < md->max_run; i++)
00120     fprintf(pf, "%.52e ", md->time_vec[i]);
00121   fprintf(pf, "\n");
00122   fclose(pf);
00123   return 0;
00124 }
00125 
00126 gs_pm_model_t *
00127 gs_pm_load(int fd)
00128 {
00129   int i, res;
00130   gs_pm_model_t *md, dub;
00131   res = read(fd, &(dub.nb_params), sizeof(dub.nb_params));
00132   if(res == -1)
00133     return NULL;
00134   res = read(fd, &(dub.nb_categories), sizeof(dub.nb_categories));
00135   if(res == -1)
00136     return NULL;
00137   res = read(fd, &(dub.max_run), sizeof(dub.max_run));
00138   if(res == -1)
00139     return NULL;
00140   md = gs_pm_init_model(dub.nb_categories, dub.nb_params, dub.max_run);
00141   res = read(fd, &(md->run_index), sizeof(md->run_index));
00142   if(res == -1) {
00143     gs_pm_free_model(md);
00144     return NULL;
00145   }
00146   res = read(fd, &(md->full), sizeof(md->full));
00147   if(res == -1) {
00148     gs_pm_free_model(md);
00149     return NULL;
00150   }
00151   res = read(fd, md->params, sizeof(md->params[0]) * md->nb_params);
00152   if(res == -1) {
00153     gs_pm_free_model(md);
00154     return NULL;
00155   }
00156   res = read(fd, md->categories, sizeof(md->categories[0]) * md->nb_categories);
00157   if(res == -1) {
00158     gs_pm_free_model(md);
00159     return NULL;
00160   }
00161   for(i = 0; i < md->max_run; i++) {
00162     res =
00163         read(fd, md->param_mat[i],
00164              sizeof(md->param_mat[i][0]) * (md->nb_params + md->nb_categories));
00165     if(res == -1) {
00166       gs_pm_free_model(md);
00167       return NULL;
00168     }
00169   }
00170   res = read(fd, md->time_vec, sizeof(md->time_vec[0]) * md->max_run);
00171   if(res == -1) {
00172     gs_pm_free_model(md);
00173     return NULL;
00174   }
00175   md->model_ok = gs_pm_compute_model(md);
00176   return md;
00177 }
00178 
00179 gs_pm_model_t *
00180 gs_pm_load_model(char *filename)
00181 {
00182   FILE *pf;
00183   int i, j;
00184   gs_pm_model_t *md;
00185   int nb_params, nb_categories, max_run;
00186   pf = fopen(filename, "r");
00187   if(!pf)
00188     return NULL;
00189   fscanf(pf, "%d\n", &nb_params);
00190   fscanf(pf, "%d\n", &nb_categories);
00191   fscanf(pf, "%d\n", &max_run);
00192   md = gs_pm_init_model(nb_categories, nb_params, max_run);
00193   fscanf(pf, "%d\n", &(md->run_index));
00194   fscanf(pf, "%d\n", &(md->full));
00195   for(i = 0; i < md->nb_params; i++)
00196     fscanf(pf, "%le ", &(md->params[i]));
00197   fscanf(pf, "\n");
00198   for(i = 0; i < md->nb_categories; i++)
00199     fscanf(pf, "%le ", &(md->categories[i]));
00200   fscanf(pf, "\n");
00201   for(i = 0; i < md->max_run; i++) {
00202     for(j = 0; j < md->nb_params + md->nb_categories; j++)
00203       fscanf(pf, "%le ", &(md->param_mat[i][j]));
00204     fscanf(pf, "\n");
00205   }
00206   for(i = 0; i < md->max_run; i++)
00207     fscanf(pf, "%le ", &(md->time_vec[i]));
00208   fscanf(pf, "\n");
00209   md->model_ok = gs_pm_compute_model(md);
00210   fclose(pf);
00211   return md;
00212 }
00213 
00214 void
00215 gs_pm_free_model(gs_pm_model_t *md)
00216 {
00217   int i;
00218   if(md->nb_params)
00219     free(md->params);
00220   if(md->nb_categories)
00221     free(md->categories);
00222   if(md->max_run) {
00223     free(md->time_vec);
00224     for(i = 0; i < md->max_run; i++)
00225       if(md->nb_categories + md->nb_params)
00226         free(md->param_mat[i]);
00227     free(md->param_mat);
00228   }
00229   if(md->coef_vec)
00230     free(md->coef_vec);
00231   free(md);
00232 }
00233 
00234 gs_pm_model_t *
00235 gs_pm_init_model(int nb_categories, int nb_params, int nb_lines)
00236 {
00237   int i, j;
00238   gs_pm_model_t *res;
00239   res = (gs_pm_model_t *) malloc(sizeof(gs_pm_model_t));
00240   res->params = (double *) calloc(nb_params, sizeof(double));
00241   res->categories = (double *) calloc(nb_categories, sizeof(double));
00242   res->nb_categories = nb_categories;
00243   res->nb_params = nb_params;
00244   res->run_index = 0;
00245   res->max_run = nb_lines;
00246   res->param_mat = (double **) malloc(sizeof(double *) * nb_lines);
00247   for(i = 0; i < nb_lines; i++) {
00248     res->param_mat[i] =
00249         (double *) calloc(nb_categories + nb_params, sizeof(double));
00250     for(j = 0; j < nb_categories; j++)
00251       res->param_mat[i][j] = -1;
00252   }
00253   res->time_vec = (double *) calloc(nb_lines, sizeof(double));
00254   res->coef_vec = (double *) calloc(nb_lines, sizeof(double));
00255   res->full = 0;
00256   res->model_ok = 0;
00257   return res;
00258 }
00259 
00260 void
00261 gs_pm_store_timing(double time, gs_pm_model_t *md)
00262 {
00263   int i;
00264   int j;
00265   i = md->run_index;
00266   for(j = 0; j < md->nb_categories; j++)
00267     md->param_mat[i][j] = md->categories[j];
00268   for(j = 0; j < md->nb_params; j++)
00269     md->param_mat[i][j + md->nb_categories] = md->params[j];
00270   md->time_vec[i] = time;
00271   if(md->run_index == md->max_run - 1) {
00272     md->full = 1;
00273     md->run_index = 0;
00274   }
00275   else
00276     md->run_index++;
00277 
00278   md->model_ok = gs_pm_compute_model(md);
00279 }
00280 
00281 gs_pm_model_t *gs_pm_clone_model(gs_pm_model_t *md1){
00282   gs_pm_model_t *md;
00283   int i;
00284 
00285   md=gs_pm_init_model(md1->nb_categories,md1->nb_params,md1->max_run);
00286   md->run_index=md1->run_index;
00287   md->full=md1->full;;
00288   memcpy(md->params,md1->params,md1->nb_params*sizeof(double));
00289   memcpy(md->categories,md1->categories,md1->nb_categories*sizeof(double));
00290   for(i=0;i<md->max_run;i++)
00291     memcpy(md->param_mat[i],md1->param_mat[i],sizeof(double)*(md1->nb_params+md1->nb_categories));
00292   memcpy(md->time_vec,md1->time_vec,md1->max_run*sizeof(double));
00293   md->model_ok=gs_pm_compute_model(md);
00294   return md;
00295 }
00296 
00297 /* in :a model 
00298   out : table of categories and coefficients
00299   categories = NULL if no categories.
00300   Returns: number of lines of each table, 0 if no models*/
00301 int gs_pm_all_models(gs_pm_model_t *md,double ***categories, double ***coefficients){
00302   gs_pm_model_t *md_copy;
00303   int n,res,i,j;
00304   double *coef;
00305 
00306   if(md->nb_categories==0){
00307     *categories=NULL;
00308     if(gs_pm_compute_model(md)==0){
00309       coef=NULL;
00310       res=0;
00311     }else{
00312       res=1;
00313       coef=(double*)malloc(md->nb_params*sizeof(double));
00314       memcpy(coef,md->coef_vec,md->nb_params*sizeof(double));
00315     }
00316     *coefficients=(double**)malloc(sizeof(double*));
00317     *coefficients[0]=coef;
00318     return res;
00319   }
00320 
00321   
00322   md_copy=gs_pm_clone_model(md);
00323   
00324 
00325   n=md_copy->full?md_copy->max_run:md_copy->run_index-1;
00326   
00327 
00328   /* compute res the number of lines of categories and coefficients */
00329   i=0;
00330   res=0;
00331   do{
00332     memcpy(md->categories,md_copy->param_mat[i],md_copy->nb_categories*sizeof(double));
00333     if(gs_pm_compute_model(md)){
00334        res++;
00335     }
00336     for(j=i;j<n;j++){
00337       int k,ok=1;
00338       /*find lines of the matrix that correspond to the current run*/
00339       for(k=0;k<md->nb_categories;k++){
00340         if(md_copy->param_mat[j][k]!=md->categories[k]){
00341           ok=0;
00342           break;
00343         }
00344       }
00345       if(ok){
00346         md_copy->time_vec[j]=-1;
00347       }
00348     }
00349     while((md_copy->time_vec[i]==-1)&&(i<n))
00350       i++;
00351   }while(i!=n);
00352 
00353 
00354   
00355   memcpy(md_copy->time_vec,md->time_vec,n*sizeof(double));
00356 
00357   *coefficients=(double**)malloc(sizeof(double*)*res);
00358   *categories=(double**)malloc(sizeof(double*)*res);
00359 
00360   
00361 
00362   i=0;
00363   res=0;
00364   do{
00365     memcpy(md->categories,md_copy->param_mat[i],md_copy->nb_categories*sizeof(double));
00366     if(gs_pm_compute_model(md)){
00367       (*coefficients)[res]=(double*)malloc(sizeof(double)*md->nb_params);
00368       (*categories)[res]=(double*)malloc(sizeof(double)*md->nb_categories);
00369       memcpy((*coefficients)[res],md->coef_vec,md_copy->nb_params*sizeof(double));
00370       memcpy((*categories)[res],md->categories,md_copy->nb_categories*sizeof(double));
00371       res++;
00372     }
00373     for(j=i;j<n;j++){
00374       int k,ok=1;
00375       /*find lines of the matrix that correspond to the current run*/
00376       for(k=0;k<md->nb_categories;k++){
00377         if(md_copy->param_mat[j][k]!=md->categories[k]){
00378           ok=0;
00379           break;
00380         }
00381       }
00382       if(ok){
00383         md_copy->time_vec[j]=-1;
00384       }
00385     }
00386     while((md_copy->time_vec[i]==-1)&&(i<n))
00387       i++;
00388   }while(i!=n);
00389 
00390   gs_pm_free_model(md_copy);
00391 
00392   return res;
00393 }
00394 
00395     /*
00396      * retruns 0 if no model 1 if able to compute a model
00397      */
00398 int
00399 gs_pm_compute_model(gs_pm_model_t *md)
00400 {
00401   double *mat, *vec;
00402   int i, j, M, N, info, l, n;
00403   N = md->nb_params;
00404   mat = (double *) calloc(md->max_run * N, sizeof(double));
00405   vec = md->coef_vec;
00406   n = md->full ? md->max_run : md->run_index;
00407   for(i = 0, l = 0; i < n; i++) {
00408     int k, ok = 1;
00409 
00410     /*
00411      * find lines of the matrix that correspond to the current run 
00412      */
00413     for(k = 0; k < md->nb_categories; k++) {
00414       if(md->param_mat[i][k] != md->categories[k]) {
00415         ok = 0;
00416         break;
00417       }
00418     }
00419 
00420     /*
00421      * if ok, add this line to the input matrix of the solver 
00422      */
00423     if(ok) {
00424       vec[l] = md->time_vec[i];
00425       for(j = 0; j < N; j++)
00426         mat[j*md->max_run+l]=md->param_mat[i][j+md->nb_categories];
00427       l++;
00428     }
00429   }
00430   M = l;
00431 
00432   /*
00433    * printf("M=%d, N=%d n=%d\n",M,N,n);
00434    * 
00435    * for(i=0;i<M;i++){ for(j=0;j<N;j++) printf("%.2e
00436    * ",mat[j*md->max_run+i]); printf("\n"); }
00437    * 
00438    * for(i=0;i<M;i++) printf("%.2e ",vec[i]);
00439    * 
00440    * printf("\n"); 
00441    */
00442   if(M < 1) {
00443     free(mat);
00444     return 0;
00445   }
00446 
00447 #ifdef DGELS
00448   {
00449     double *work;
00450     int nrhs = 1, lwork = 1;
00451     if((!(md->full)) && (M < md->nb_params + md->nb_categories)) {
00452       free(mat);
00453       return 0;
00454     }
00455 
00456     /*
00457      * Workspace query 
00458      */
00459     work = (double *) malloc(N * sizeof(double));
00460     if(work == NULL) {
00461       printf("error of memory allocation\n");
00462       exit(0);
00463     }
00464     lwork = -1;
00465     dgels_("N", &M, &N, &nrhs, mat, &(md->max_run), vec, &(md->max_run), work,
00466            &lwork, &info);
00467     lwork = (int) work[0];
00468     free(work);
00469 
00470     /*
00471      * allocation of the workspace 
00472      */
00473     work = (double *) malloc(lwork * sizeof(double));
00474     if(work == NULL) {
00475       printf("error of memory allocation\n");
00476       exit(0);
00477     }
00478     dgels_("N", &M, &N, &nrhs, mat, &(md->max_run), vec, &(md->max_run),
00479            work, &lwork, &info);
00480     if(info != 0) {
00481       printf("\nError on parameter %d!\n", -info);
00482       exit(0);
00483     }
00484     free(work);
00485     free(mat);
00486     return 1;
00487   }
00488 
00489 #else                           /* */
00490   {
00491     double *x, rnorm, *w, *zz;
00492     int *index;
00493     x = (double *) malloc(N * sizeof(double));
00494     if(x == NULL) {
00495       printf("error of memory allocation\n");
00496       exit(0);
00497     }
00498     w = (double *) malloc(N * sizeof(double));
00499     if(w == NULL) {
00500       printf("error of memory allocation\n");
00501       exit(0);
00502     }
00503     zz = (double *) malloc(M * sizeof(double));
00504     if(zz == NULL) {
00505       printf("error of memory allocation\n");
00506       exit(0);
00507     }
00508     index = (int *) malloc(N * sizeof(int));
00509     if(index == NULL) {
00510       printf("error of memory allocation\n");
00511       exit(0);
00512     }
00513     nnls_c(mat, &(md->max_run), &M, &N, vec, x, &rnorm, w, zz, index, &info);
00514 
00515     /*
00516      * printf("info : %d, norm : %f\n",info,rnorm); 
00517      */
00518 
00519     free(w);
00520     free(zz);
00521     free(index);
00522     free(mat);
00523     for(i = 0; i < N; i++)
00524       vec[i] = x[i];
00525     free(x);
00526     return 1;
00527   }
00528 
00529 #endif                          /* */
00530 }
00531 
00532 double
00533 gs_pm_guess_time(gs_pm_model_t *md)
00534 {
00535   double *vec;
00536   double time = 0;
00537   int i;
00538   if(!md->model_ok)
00539     return -1;
00540   vec = md->coef_vec;
00541 
00542   /*
00543    * printf("\n"); for(i=0;i<nb_params;i++) printf("%e*%.0f
00544    * ",vec[i],param_vec[i]); printf("\n"); 
00545    */
00546   for(i = 0; i < md->nb_params; i++) {
00547     time+=vec[i]*md->params[i];
00548 
00549     /*
00550      * printf("%e*%f ",vec[i],param_vec[i]); 
00551      */
00552   }
00553 
00554   /*
00555    * printf("\n"); 
00556    */
00557   return time;
00558 }
00559 
00560 void
00561 gs_pm_display_model(gs_pm_model_t *md)
00562 {
00563   double *vec;
00564   int i;
00565   if(md == NULL)
00566     return;
00567   if(!md->model_ok)
00568     return;
00569   vec = md->coef_vec;
00570   for(i = 0; i < md->nb_params; i++)
00571     printf("%.1e ",vec[i]);
00572   printf("\n");
00573 }