util.cpp

Go to the documentation of this file.

#include <math.h>
#include <iostream>
#include <fstream>   // file I/O
#include <float.h>
#include <limits.h>
#include <algorithm>
#include <string>
#include <cctype>
#ifdef PARALLEL
   #include <mpi.h>
#elif defined(WIN)
  #include <time.h>
#else
   #include <sys/time.h>
#endif

#include "util.h"

using namespace std;

bool util::strcmp_ignore_case(string first, string second) {

  // Convert both strings to upper case by transfrom() before compare.
  std::transform(first.begin(), first.end(), first.begin(), static_cast < int(*)(int) > (toupper));
  std::transform(second.begin(), second.end(), second.begin(), static_cast < int(*)(int) > (toupper));

  if (first == second) 
        return true; 
  else 
        return false;
}

int util::current_time(double *t) {
    
    
    #ifdef PARALLEL

        *t=MPI_Wtime();                
    #elif defined(__GNUG__)

        struct timeval tv;
        struct timezone tz;

        gettimeofday(&tv,&tz);
        *t=(double)tv.tv_sec+(double)tv.tv_usec/1000000.0;      
    #elif defined(_WIN32)
                *t=clock() / CLOCKS_PER_SEC;
    #endif 


    return 0;    
}


int util::get_proc_id_nump(int &proc_id,int &nump)   {
    proc_id=0; nump=1;    
    #ifdef PARALLEL
        MPI_Comm_size(MPI_COMM_WORLD,&nump); 
        MPI_Comm_rank(MPI_COMM_WORLD,&proc_id);  
    #endif 
    return 0;
}

string util::get_time_str() {

        time_t rawtime;
        struct tm * timeinfo;
        time ( &rawtime );
        timeinfo = localtime ( &rawtime );
        string str=asctime(timeinfo);
        str.erase(std::remove(str.begin(), str.end(), '\n'), str.end());

        return str;
}
void util::fast_gcd(double *a, t_UnsignedInt len, double *gcd) {
  double g=a[0],y,r;
  t_UnsignedInt k;
  
  g=a[0];

  for (k = 0;k<len-1;k++) {

    y = a[k+1];
    while (y > 0) {
      r = g - y *  floor(g / y);
      g = y;
      y = r;
    }
    if (g == (double)1.0) {
      *gcd=g;
      return;
      }
    
  }
   *gcd=g;

} 
double util::infnormvec(vector< vector<double> > A) {
       t_UnsignedInt i, j;

       double s,max_s;
       if (A.size()==0)
                return 0;
        
       t_UnsignedInt m, n;
       m=A.size(); n=A[0].size();

       for (j=0,max_s=0;j<n;j++)
           max_s+=fabs(A[0][j]);
           
       for (i=1;i<m;i++) {
           for (j=0,s=0;j<n;j++)
               s+=fabs(A[i][j]);
               
           if (s>max_s)
              max_s=s;
       }
       
       return max_s;
}
int util::right_nullspace_integer(vector< vector<double> > &A, vector< vector<double> > &K) {

    t_UnsignedInt m, n;
    m=A.size(); n=A[0].size();

    double tol,eps=2.2204e-16;
    vector< vector<double> > tab((m+n), vector<double>((n),0));  
    vector<double> mc;
    vector<t_UnsignedInt> ic,rws, re_rows, re_cols,unused;

    t_UnsignedInt  i,j,k ;
    double *a_gcd, piv_el,min_el;
    t_UnsignedInt piv_c_ind,piv_r_ind,piv_r,piv_c;
    t_UnsignedInt col;

    vector<t_UnsignedInt> search_cols; //, search_cols_size;

    t_UnsignedInt tmp_el,kn_cols=0,all_zeros;
    double g, fp, fo;
        
    if (m==0 || n==0) {
        return 0;
    }

    for (i=0;i<m+n;i++) {
            if (i<m) {
                 for (j=0;j<n;j++)
                        tab[i][j]=A[i][j];           
            } else {
                 for (j=0;j<n;j++)
                        tab[i][j]=0;           
                 tab[i][i-m]=1;
            }
    }
        


    tol=0;

    re_rows.resize(m);
    re_cols.resize(n);

    for (i=0;i<re_rows.size();i++)
            re_rows[i]=i;
    for (i=0;i<re_cols.size();i++)
            re_cols[i]=i;
    for (i=0;i<n;i++)
            unused.push_back(1);

        //-----------------
    while (re_rows.size()>0 && re_cols.size()>0) {
                
        //if (prior_cols.size()==0) {
                        search_cols.resize(re_cols.size());
                        for (i=0;i<re_cols.size();i++)
                                search_cols[i]=re_cols[i];
        /*
        } else {
                //search_cols=prior_cols;
                search_cols.resize(prior_cols.size());
                for (i=0;i<prior_cols.size();i++)
                        search_cols[i]=prior_cols[i];
        }
        */              
                
                       
                        //tmp=fabs(tab(re_rows, search_cols));
                        vector< vector<double> > tmp(re_rows.size(),vector<double>(search_cols.size()));

                        for (i=0;i<re_rows.size();i++) {
                            for (j=0;j<search_cols.size();j++) {
                                tmp[i][j]=fabs(tab[re_rows[i]][search_cols[j]]);
                                    if (tmp[i][j]==0) 
                                          tmp[i][j]=FLT_MAX;   //tmp(tmp==0)=+inf;                              
                            }
                        }                                       

                        //[mc,ic]=min(tmp);  //mc is min value in a column, ic is the index of the row of the max value in the column                   
                        mc.resize(search_cols.size());
                        ic.resize(search_cols.size());

                        for (j=0;j<search_cols.size();j++)
                             for (i=1,mc[j]=tmp[0][j],ic[j]=0;i<re_rows.size();i++) {
                                if (tmp[i][j]<mc[j]) {
                                    mc[j]=tmp[i][j];
                                    ic[j]=i;   
                                 }
                             }

                        if (re_rows.size()==1) {
                           for (i=1,min_el=mc[0];i<search_cols.size();i++)
                              if (mc[i]<min_el) {
                                    ic[0]=i;               
                                    min_el=mc[i];
                              }
                           //mc_size=1;
                           mc.resize(1);
                           mc[0]=min_el;    
                        }

                        //[piv_el,piv_c_ind]=min(mc);                            
                        for (i=1,piv_el=mc[0],piv_c_ind=0;i<mc.size();i++)
                            if (mc[i]<piv_el) {
                               piv_el=mc[i];
                               piv_c_ind=i;                  
                            }
                        //piv_r_ind=ic(piv_c_ind);
                        piv_r_ind=ic[piv_c_ind];
                        
                        if (piv_el==FLT_MAX)
                                piv_el=0;

                        //for (i=0;i<re_rows.size();i++) 
                        //    custom_free(tmp[i]); 
                        //custom_free(tmp); 
                
                if (re_rows.size()==1) {
                        piv_r_ind=0;
                        piv_c_ind=ic[0];
                }
                ic.clear();

                piv_c=search_cols[piv_c_ind];
                piv_r=re_rows[piv_r_ind];
                
                /*
                if (prior_cols.size()!=0) {
                    //prior_cols[piv_c_ind]=[];
                    for (i=piv_c_ind;i<prior_cols.size()-1;i++)
                        prior_cols[i]=prior_cols[i+1];
                    prior_cols.resize(prior_cols.size()-1);
                    
                }
                */
                
                //j=find(re_cols==piv_c)  POSSIBLY INCORECT
                for (i=0;i<re_cols.size();i++)
                   if (re_cols[i]==piv_c) {
                      j=i;
                      break;
                   }
                //re_cols[j]=[]
                for (i=j;i<re_cols.size()-1;i++)
                    re_cols[i]=re_cols[i+1];
                re_cols.resize(re_cols.size()-1);
                
                if (piv_el<=tol) {
                   for (i=0;i<re_rows.size();i++)
                        tab[re_rows[i]][piv_c]=0;
                   continue;                 
                }    
        
                unused[piv_c]=0;
                //unused[piv_c]=1;
                piv_el=tab[piv_r][piv_c];
       
                //rws=custom_malloc((re_rows_size+tab_m-m),t_UnsignedInt,__FILE__,__LINE__);
                //rws_size=re_rows_size+tab_m-m;
                rws.resize(re_rows.size()+n); //tab_m-m);
                //rws= [re_rows, (m+1):tab_m];
                for (i=0;i<rws.size();i++) {
                        if (i<re_rows.size())
                                rws[i]=re_rows[i];
                        else
                                rws[i]=i-re_rows.size()+m;
        

                }
                
                for (i=piv_r_ind;i<re_rows.size()-1;i++)
                    re_rows[i]=re_rows[i+1];
                re_rows.resize(re_rows.size()-1);
                      
                unused[piv_c]=0;
                //unused[piv_c]=1;
                
                for (j=0;j<re_cols.size();j++) {
                    col=re_cols[j];
                   
                        
                       //a_gcd=custom_malloc(10,double,__FILE__,__LINE__);
                        a_gcd=new double[10];
                       //a_gcd[0]=1; a_gcd[1]=2;
                       a_gcd[0]=fabs(tab[piv_r][col]);                  
                       a_gcd[1]=fabs(piv_el);
                       fast_gcd(a_gcd,2,&g); //modified                 
                        delete [] a_gcd; //  delete [] a_gcd; //free(a_gcd);
                       fp=tab[piv_r][col]/g;
                       fo=piv_el/g;
                           //printf tab

                       for (i=0;i<rws.size();i++) 
                           tab[rws[i]][col]=fo*tab[rws[i]][col]-fp*tab[rws[i]][piv_c];
                        
                       //a_gcd=custom_malloc(rws.size(),double,__FILE__,__LINE__);
                       a_gcd=new double[rws.size()];
                       for (i=0;i<rws.size();i++)
                           a_gcd[i]=fabs(tab[rws[i]][col]);
                       fast_gcd(a_gcd,rws.size(),&g); //modified

                       for (i=0;i<rws.size();i++)
                           tab[rws[i]][col]/=g;                           
                        
                        delete [] a_gcd; //delete [] a_gcd; // free(a_gcd);                                                                
                
                }
                rws.clear();
                //custom_free(rws); //delete [] rws;//    free(rws);  //FREED
        }    

        if (m == 1) {
           // K= tab(2:tab_m, ~tab(1, :));
           
           // tab(1..m+n-1; 
           // TO DO
           cout << endl << "here";
          
        } else {
            //K= tab((m+1):tab_m, ~any(tab(1:m, :)));
            for (i=0,kn_cols=0;i<n;i++) {
                for (j=0,all_zeros=1;j<m;j++)
                        if (tab[j][i]) {
                                all_zeros=0;
                                break;
                        }
                if (all_zeros)
                        kn_cols++;
            }

            vector<t_UnsignedInt> kn_cols_ind;
            //t_UnsignedInt *kn_cols_ind=custom_malloc(kn_cols,t_UnsignedInt,__FILE__,__LINE__);
            kn_cols_ind.resize(kn_cols);


            for (i=0,kn_cols=0;i<n;i++) {
                for (j=0,all_zeros=1;j<m;j++)
                        if (tab[j][i]) {
                                all_zeros=0;
                                break;
                        }
                if (all_zeros)
                        kn_cols_ind[kn_cols++]=i;
            }

            K.resize(kn_cols, vector<double>(n));
            for (i=0;i<kn_cols;i++) 
                    for (j=0;j<n;j++)
                        K[i][j]=tab[j+m][kn_cols_ind[i]];       
        }

        t_UnsignedInt ii=0;
        for (i=0, ii=0;i<n;i++)
                if (unused[i])
                        unused[ii++]=i;
        unused.resize(ii);
        tab.clear();
        re_rows.clear();
        re_cols.clear();

        double elem;
        for (i=0;i<unused.size();i++) {
                //inspect row K[unused[i]]
                for (j=0;j<K.size();j++) {
                        if (K[j][unused[i]]!=0) {
                                col=j;
                                elem=K[j][unused[i]];
                                break;
                        }
                }
                for (j=0;j<n;j++)
                        K[col][j]=K[col][j]/elem;
        }
        /*
        double big_sum=0;
        for (i=0;i<A.size();i++)
                for (j=0;j<K.size();j++) {
                        double part_sum;
                        for ( k=0,  part_sum=0;k<A[0].size();k++)
                                part_sum+=A[i][k]*K[j][k];
                        big_sum+=fabs(part_sum);
                }
        */
        
        //cout << endl <<"(right_nullspace_integer) big_sum: "<<big_sum<<endl;


        return 0;       
}
int util::right_nullspace_nonint(vector< vector<double> > &A,  vector< vector<double> > &K) {

    t_UnsignedInt m, n;
    m=A.size(); n=A[0].size();

    double tol,eps=1e-15; //2.2204e-16; //, **tmp;
    vector< vector<double> > tab((m+n), vector<double>((n),0));  
    vector<double> mc;
    vector<t_UnsignedInt> ic,rws, re_rows, re_cols;

    t_UnsignedInt  i,j,k ;
    double  piv_el,min_el;
    t_UnsignedInt piv_c_ind,piv_r_ind,piv_r,piv_c;
    t_UnsignedInt col;

    vector<t_UnsignedInt> search_cols; //, search_cols_size;

    t_UnsignedInt kn_cols=0,all_zeros;
    double tmp_el;

    if (m==0 || n==0) 
        return 0;
    

    for (i=0;i<m+n;i++) {
            if (i<m) {
                 for (j=0;j<n;j++)
                        tab[i][j]=A[i][j];           
            } else {
                 for (j=0;j<n;j++)
                        tab[i][j]=0;           
                 tab[i][i-m]=1;
            }
    }



    tol=eps*max(m,n)*infnormvec(A); //inf The infinity norm, or largest row sum of A, max(sum(abs(A'))).

    re_rows.resize(m);
    re_cols.resize(n);

    for (i=0;i<re_rows.size();i++)
            re_rows[i]=i;
    for (i=0;i<re_cols.size();i++)
            re_cols[i]=i;

    while (re_rows.size()>0 && re_cols.size()>0) {
                
                search_cols.resize(re_cols.size());
                for (i=0;i<re_cols.size();i++)
                        search_cols[i]=re_cols[i];


        //[mc,ic]=max(tmp);  //mc is min value in a column, ic is the index of the row of the max value in the column                   
        mc.resize(search_cols.size());                   
        ic.resize(search_cols.size());
        
        for (j=0;j<search_cols.size();j++) {
             //mc[j]=tmp[re_rows[0]][search_cols[j]];
             mc[j]=fabs(tab[re_rows[0]][search_cols[j]]);
             ic[j]=0;
             for (i=1;i<re_rows.size();i++) {
                 double elem=fabs(tab[re_rows[i]][search_cols[j]]);     
                 if (elem>mc[j]) {
                
                        mc[j]=elem;
                        ic[j]=i;
                
                 }
             }
        }

         //[piv_el,piv_c_ind]=max(mc);                           
         for (i=1,piv_el=mc[0],piv_c_ind=0;i<search_cols.size();i++)
            if (mc[i]>piv_el) {
                   piv_el=mc[i];
                   piv_c_ind=i;                  
            }
        

         //piv_r_ind=ic(piv_c_ind);
         piv_r_ind=ic[piv_c_ind];

         mc.clear();
        
         if (re_rows.size()==1) {
                piv_r_ind=0;
                piv_c_ind=ic[0];
         }
         ic.clear();

         piv_c=search_cols[piv_c_ind];
         piv_r=re_rows[piv_r_ind];
                
         //j=find(re_cols==piv_c)  POSSIBLY INCORECT
         for (i=0;i<re_cols.size();i++)
               if (re_cols[i]==piv_c) {
                      j=i;
                      break;
               }
         //re_cols[j]=[]
         for (i=j;i<re_cols.size()-1;i++)
              re_cols[i]=re_cols[i+1];
         re_cols.resize(re_cols.size()-1);
                
         if (piv_el<=tol) {
             for (i=0;i<re_rows.size();i++)
                tab[re_rows[i]][piv_c]=0;
               continue;                 
         }    
        
         piv_el=tab[piv_r][piv_c];
         
         //rws_size=re_rows_size+tab_m-m;
         rws.resize(re_rows.size()+n); //tab_m-m);
         //rws= [re_rows, (m+1):tab_m];

         for (i=0;i<rws.size();i++) {
                if (i<re_rows.size())
                        rws[i]=re_rows[i];
                else
                        rws[i]=i-re_rows.size()+m;

         }

         
         for (i=0;i<rws.size();i++) 
             tab[rws[i]][piv_c]=tab[rws[i]][piv_c]/piv_el;         

         for (i=piv_r_ind;i<re_rows.size()-1;i++)
            re_rows[i]=re_rows[i+1];
         re_rows.resize(re_rows.size()-1);                   

         for (j=0;j<re_cols.size();j++) {
               col=re_cols[j];          
               tmp_el=tab[piv_r][col];
               for (i=0;i<rws.size();i++) 
                   tab[rws[i]][col]-=tmp_el*tab[rws[i]][piv_c];                
               for (i=0;i<rws.size();i++) 
                    if (fabs(tab[rws[i]][col])<=tol)
                              tab[rws[i]][col]=0;                                   
         }
        
         rws.clear();

        }    

        if (m == 1) {
           // K= tab(2:tab_m, ~tab(1, :));
           
           // tab(1..m+n-1; 
           // TO DO
          
        } else {
            //K= tab((m+1):tab_m, ~any(tab(1:m, :)));
                /*
            for (i=0,kn_cols=0;i<n;i++) {
                for (j=0,all_zeros=1;j<m;j++)
                        if (tab[j][i]) {
                                all_zeros=0;
                                break;
                        }
                if (all_zeros)
                        kn_cols++;
            }
        */

            vector<t_UnsignedInt> kn_cols_ind;
            //t_UnsignedInt *kn_cols_ind=custom_malloc(kn_cols,t_UnsignedInt,__FILE__,__LINE__);
            kn_cols_ind.resize(kn_cols);


            for (i=0,kn_cols=0;i<n;i++) {
                for (j=0,all_zeros=1;j<m;j++)
                        if (tab[j][i]) {
                                all_zeros=0;
                                break;
                        }
                if (all_zeros)
                        kn_cols_ind.push_back(i);
                        //kn_cols_ind[kn_cols++]=i;
            }

            //K.resize(kn_cols, vector<double>(n));
            K.resize(kn_cols_ind.size(), vector<double>(n));
            for (i=0;i<kn_cols_ind.size();i++) 
                    for (j=0;j<n;j++)
                        K[i][j]=tab[j+m][kn_cols_ind[i]];       
        }

        return 0;       
}

t_UnsignedInt util::row_echelon(vector< vector<double> > mat, t_UnsignedInt reduced, vector<t_UnsignedInt> rowmap, vector<t_UnsignedInt> colmap) {

        t_UnsignedInt m,n;
        m=mat.size();n=mat[0].size();
        t_UnsignedInt rows=m, cols=n;
        t_UnsignedInt prows=(rowmap.size()==0)? m : minval(m,rowmap.size());
        t_UnsignedInt pcols=(colmap.size()==0)? n : minval(n,colmap.size());
        t_UnsignedInt pivs=minval(prows,pcols);
        
        double mTolerance=1e-15;
        t_UnsignedInt tmp;
        
        t_UnsignedInt prow=-1, pcol=-1;
        //int prow=-1, pcol=-1;
        double pval=-DBL_MAX;
        for (t_UnsignedInt row=0;row<prows;row++) {
                for (t_UnsignedInt col=0;col<pcols;col++) {
                        double val=fabs(mat[row][col]);
                        if (val>pval) {
                                pval=val;
                                prow=row;
                                pcol=col;
                        }
                }
        }

        
        for (t_UnsignedInt pivot=0;pivot<pivs;pivot++) {
                if (pval<=mTolerance)
                        return pivot;
                if (prow!=pivot) {
                        swap_rows(mat,prow,pivot);
                        if (rowmap.size()>0) {
                                tmp=rowmap[prow];
                                rowmap[prow]=rowmap[pivot];
                                rowmap[pivot]=tmp;
                        }
                }
                if (pcol!=pivot) {
                        swap_columns(mat,pcol,pivot);
                        if (colmap.size()>0) {
                                tmp=colmap[pcol];
                                colmap[pcol]=colmap[pivot];
                                colmap[pivot]=tmp;
                        }
                }
                pval=mat[pivot][pivot];
                for (t_UnsignedInt col=pivot+1;col<cols;col++) 
                        mat[pivot][col]*=1.0/pval;
                mat[pivot][pivot]=1.0;
                pval=-DBL_MAX;
                for (t_UnsignedInt row=pivot+1;row<rows;row++) {
                        double rpiv=mat[row][pivot];
                        mat[row][pivot]=0;
                        for (t_UnsignedInt col=pivot+1;col<cols;col++) {
                                double val=mat[row][col];
                                double sub=mat[pivot][col];
                                val-=sub*rpiv;
                                mat[row][col]=val;
                                if (row<prows && col<pcols) {
                                        if (val<0) val=-val;
                                        if (val>pval) {
                                                pval=val;
                                                prow=row;
                                                pcol=col;
                                        }
                                }
                        }
                }
                if (reduced) {
                        for (t_UnsignedInt row=0;row<pivot;row++) {
                                double rpiv=mat[row][pivot];
                                mat[row][pivot]=0;
                                for (t_UnsignedInt col=pivot+1;col<cols;col++) {
                                        double val=mat[row][col];
                                        double sub=mat[pivot][col];
                                        mat[row][col]=val-sub*rpiv;
                                }
                        }
                }
        }
        
        return pivs;
        
}

int util::swap_rows(vector< vector<double> > &mat, t_UnsignedInt row1, t_UnsignedInt row2) {

        double tmp;
        t_UnsignedInt m,n;

        m=mat.size();
        if (mat.size()==0)
                return -1;
        n=mat[0].size();

        if ((row1>=m) || (row2>=m)) {
                printf("\nswap_rows: Incorrect row values.");
                return -1;
        }
        
        for (t_UnsignedInt i=0;i<n;i++) {
                tmp=mat[row1][i];
                mat[row1][i]=mat[row2][i];
                mat[row2][i]=tmp;
        }       
        return 0;
}
int util::swap_columns(vector< vector<double> > &mat, t_UnsignedInt col1, t_UnsignedInt col2) {

        double tmp;
        t_UnsignedInt i,m,n;

        m=mat.size();
        if (mat.size()==0)
                return -1;
        n=mat[0].size();


        if ((col1>=n) || (col2>=n)) {
                printf("\nswap_columns: Incorrect column index values.");
                return -1;
        }
        
        for (i=0;i<m;i++) {
                tmp=mat[i][col1];
                mat[i][col1]=mat[i][col2];
                mat[i][col2]=tmp;
        }
        return 0;
}

int util::nullspace(vector< vector<double> > mat,  vector< vector<double> > &kernel) {

         t_UnsignedInt m, n,ndim;
         //find dimensions of matrix (m,n). CONTROL missing!
         m=mat.size(); n=mat[0].size();


        t_UnsignedInt i,j,rank;
        vector< vector<double> > rref;
        vector<t_UnsignedInt> rowmap,colmap; //*rowmap, rowmap_len, *colmap, colmap_len;
        
        //rowmap_len=0;rowmap=NULL;
        //colmap_len=n;
        rowmap.resize(0);
        colmap.resize(n);

        //colmap=custom_malloc(colmap_len,t_UnsignedInt,__FILE__,__LINE__);
        for (i=0;i<colmap.size();i++)
                colmap[i]=i;
        
        //rref=custom_malloc(m,double*,__FILE__,__LINE__);
        rref.resize(m,vector<double>(n,0));

        for (i=0;i<m;i++) {
                //rref[i]=custom_malloc(n,double,__FILE__,__LINE__);
                for (j=0;j<n;j++)
                        rref[i][j]=mat[i][j];
        }
        
        //find row echelon form of the matrix mat.
        rank=row_echelon(rref,1,rowmap,colmap);
        //store nullity of the matrix mat in ndim
        ndim=n-rank;

        //allocate space for nullspace matrix
        kernel.resize(n,vector<double>(ndim,0));
        for (i=0;i<n;i++)  
                for (j=0;j<ndim;j++)
                        kernel[i][j]=0;



        //initialize values in the nullspace matrix 'kernel'
        for (t_UnsignedInt row=0;row<rank;row++) {
                for (t_UnsignedInt col=0;col<ndim;col++) {
                        kernel[colmap[row]][col]=-rref[row][col+rank];
                }
        }
        for (t_UnsignedInt row=0;row<ndim;row++)
                kernel[colmap[row+rank]][row]=1.0;

        
        return 0;

}



bool util::file_exists(const string &fname ) {   

        return ifstream( fname.c_str() ).is_open();

}