text: /****************************************************************************
   ****************************************************************************
   *                                                                          *
   *  Copyright (C) 2018  Genome Research Ltd.                                *
   *                                                                          *
   *  Author: Zemin Ning (zn1@sanger.ac.uk)                                   *
   *                                                                          *
   *  This file is part of ScaffHiC pipeline.                                 *
   *                                                                          *
   *  ScaffHiC is a free software: you can redistribute it and/or modify it   *
   *  under the terms of the GNU General Public License as published by the   *
   *  Free Software Foundation, either version 3 of the License, or (at your  *
   *  option) any later version.                                              *
   *                                                                          *
   *  This program is distributed in the hope that it will be useful, but     *
   *  WITHOUT ANY WARRANTY; without even the implied warranty of              *
   *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU        *
   *  General Public License for more details.                                *
   *                                                                          *
   *  You should have received a copy of the GNU General Public License along *
   *  with this program.  If not, see <http://www.gnu.org/licenses/>.         *
   *                                                                          *
   ****************************************************************************
   ****************************************************************************/
  /****************************************************************************/
  
   
  #include <math.h>
  #include <values.h>
  #include <stdio.h>
  #include <netinet/in.h>
  #include <stdlib.h>
  #include <dirent.h>
  #include <string.h>
  #include <ctype.h>
  #include <unistd.h>
  #include "fasta.h"
  
  #define GT '>'
  #define GT4 (((((GT<<8)+GT)<<8)+GT)<<8)+GT
  
  #define ENDS_EXTRA 0
  #define PADCHAR '-'
  #define MAX_N_BRG 50000 
  #define MAX_N_ROW 40000 
  #define Max_N_NameBase 400 
  #define Max_N_NameBase2 400 
  #define Max_N_Pair 100
  static char bindir[2000];
  static char tmpdir[2000];
  static char **S_Name,**R1_Name,**R2_Name,**R1_Name1,**R2_Name1,**R1_Name2,**R2_Name2;
  static int *insert_siz,*insert_dev,*core_list,*ctg_list,*ctg_head,*read2contig;
  static int *readIndex,*reads_tag1,*reads_tag2;
  
  /* SSAS default parameters   */
  static int n_group=0;
  static int num_reads=10;
  static int num_links=8;
  static int len_block=50000;
  static int len_edges=50000;
  static int len_matrx=2000;
  static int min_edge = 5;
  static int gap_len = 100;
  static int uplinks = 50;
  static int mscore = 20;
  static int n_longread = 1;
  static int n_files = 1000;
  static int run_align = 1;
  static int min_len = 3000;
  
  typedef struct
  {
         int foffset;
         int fsindex;
  } SIO;
  
  
  int main(int argc, char **argv)
  {
      int i,nSeq,args,n_files,n_r1,n_r2;
      char *st,*ed;
      char **cmatrix(B64_long nrl,B64_long nrh,B64_long ncl,B64_long nch);
      void ArraySort_String2(int n,char **Pair_Name,int *brr);
      fasta *seq;
      FILE *fp,*namef,*namef2;
      int size_file;
      int m_score,n_nodes,n_debug,num_sigma;
      void decodeReadpair(int nSeq);
      void HashFasta_Head(int i, int nSeq);
      void HashFasta_Table(int i, int nSeq);
      void Search_SM(fasta *seq,int nSeq);
      void Assemble_SM(int arr,int brr);
      void Read_Index(int nSeq,char *namefile);
      void Read_Group(fasta *seq,int nSeq,int nRead,int cindex);
      void File_Output(int aaa);
      void Memory_Allocate(int arr);
      char* substring(char*, int, int);
      char readname[2000],line[2000],tempa[2000],tempc[2000],syscmd[2000],file_tarseq[2000],file_scaff[2000],workdir[2000];
      char reads_file[2000],file_read1[2000],file_read2[2000],samname[500],toolname[500];
      int systemRet = system (syscmd);
      int systemChd = chdir(tmpdir);
      pid_t pid;
  
      seq=NULL;
      fflush(stdout);
      if(system("ps aux | grep scaffHiC; date") == -1)
      {
  //        printf("System command error:\n);
      }
  
      if(argc < 2)
      {
           printf("Usage: %s <Input_read_file> <Output_read1> <Output_read2>\n",argv[0]);
           printf("q1=/lustre/scratch117/sciops/team117/hpag/zn1/project/fish/reads/test/22047_4_AACCGTAA_S16_L004_R1_001.fastq.gz\n");
           printf("q2=/lustre/scratch117/sciops/team117/hpag/zn1/project/fish/reads/test/22047_4_AACCGTAA_S16_L004_R2_001.fastq.gz\n");
           printf("...\n");
           exit(1);
      }
  
      m_score = 50;
      n_nodes = 20;
      n_debug = 0;
  
      strcpy(toolname,"bwa");
      nSeq=0;
      args=1;
      for(i=1;i<argc;i++)
      {
         if(!strcmp(argv[i],"-gap"))
         {
           sscanf(argv[++i],"%d",&gap_len); 
           args=args+2;
         }
         else if(!strcmp(argv[i],"-help"))
         {
           printf("Usage: %s <Input_read_file> <Output_read1> <Output_read2>\n",argv[0]);
           printf("q1=/lustre/scratch117/sciops/team117/hpag/zn1/project/fish/reads/test/22047_4_AACCGTAA_S16_L004_R1_001.fastq.gz\n");
           printf("q2=/lustre/scratch117/sciops/team117/hpag/zn1/project/fish/reads/test/22047_4_AACCGTAA_S16_L004_R2_001.fastq.gz\n");
           printf("...\n");
           exit(1);
         }
         else if(!strcmp(argv[i],"-debug"))
         {
           sscanf(argv[++i],"%d",&n_debug);
           args=args+2;
         }
         else if(!strcmp(argv[i],"-link"))
         {
           sscanf(argv[++i],"%d",&num_links); 
           args=args+2;
         }
      }
  
      pid = getpid();
      memset(tempa,'\0',2000);
      if (!getcwd(tempa, sizeof(tempa)))
      {
        exit(1);
      } 
      memset(tmpdir,'\0',2000);
      memset(workdir,'\0',2000);
      sprintf(tmpdir,"%s/",tempa);
      sprintf(workdir,"%s/tmp_rununik_%d/",tempa,pid);
  
      memset(syscmd,'\0',2000);
      sprintf(syscmd,"mkdir %s",workdir);
      if(system(syscmd) == -1)
      {
        printf("System command error:\n");
      }
  //    system (syscmd);
  
  
      if((fp = fopen(argv[args],"r")) == NULL)
      {
        printf("ERROR main:: reads group file \n");
        exit(1);
      }
      n_files = 0;
      while(!feof(fp))
      {
        if(fgets(line,2000,namef) == NULL)
        {
  //        printf("fgets command error:\n);
        }
        if(feof(fp)) break;
        n_files++;
      }
      fclose(fp);
  
      if((reads_tag1 = (int *)calloc(n_files,sizeof(int))) == NULL)
      {
        printf("fmate: calloc - reads_tag\n");
        exit(1);
      }
      if((reads_tag2 = (int *)calloc(n_files,sizeof(int))) == NULL)
      {
        printf("fmate: calloc - reads_tag\n");
        exit(1);
      }
  
      R1_Name=cmatrix(0,n_files+10,0,2000);
      R2_Name=cmatrix(0,n_files+10,0,2000);
      R1_Name1=cmatrix(0,n_files+10,0,2000);
      R2_Name1=cmatrix(0,n_files+10,0,2000);
      R1_Name2=cmatrix(0,n_files+10,0,2000);
      R2_Name2=cmatrix(0,n_files+10,0,2000);
  
      if((fp = fopen(argv[args],"r")) == NULL)
      {
        printf("ERROR main:: reads group file \n");
        exit(1);
      }
      i = 0;
  
      n_r1 = 0;
      n_r2 = 0;
      while(fscanf(fp,"%s",readname)!=EOF)
      {
          int clen = strlen(readname);
          char *st,*ed;
          st = strchr(readname,'=');;
          ed = strrchr(readname,'/');
          if((strncmp(readname,"q1=",3))==0)
          {
            strcpy(R1_Name[n_r1],st+1);
            strcpy(R1_Name1[n_r1],ed+1);
            if((strncmp(substring(readname,clen-2,clen),".gz",3))==0)
              reads_tag1[n_r1] = 1;
            else
              reads_tag1[n_r1] = 2;
  //      printf("file: %d %d %s %s\n",clen,n_r1,R1_Name1[n_r1],substring(readname,clen-2,clen));
  
            memset(syscmd,'\0',2000);
            sprintf(syscmd,"cp %s %s",R1_Name[n_r1],workdir);
            if(system(syscmd) == -1)
            {
  //           printf("System command error:\n);
            }
  
            n_r1++;
          }
          else if((strncmp(readname,"q2=",3))==0) 
          {
            strcpy(R2_Name[n_r2],st+1);
            strcpy(R2_Name1[n_r2],ed+1);
            if((strncmp(substring(readname,clen-2,clen),".gz",3))==0)
              reads_tag2[n_r2] = 1;
            else
              reads_tag2[n_r2] = 2;
  //      printf("file: %d %s \n",n_r2,R2_Name[n_r2]);
            memset(syscmd,'\0',2000);
            sprintf(syscmd,"cp %s %s",R2_Name[n_r2],workdir);
            if(system(syscmd) == -1)
            {
  //           printf("System command error:\n);
            }
            n_r2++;
          }
          i++;
      }
      fclose(fp);
  
      if(n_r1 != n_r2)
      {
        printf("Number of read1 files: %d NOT equal to read2 files: %d\n",n_r1,n_r2);
        memset(syscmd,'\0',2000);
        sprintf(syscmd,"rm -rf %s > /dev/null",workdir);
        if(system(syscmd) == -1)
        {
  //      printf("System command error:\n);
        }
        exit(1);
      }
      
      if(chdir(workdir) == -1)
      {
        printf("System command error: chdir\n");
      }
      
      st = argv[0];
      ed = strrchr(argv[0],'/');
      memset(tempc,'\0',2000);
      strncpy(tempc,argv[0],ed-st);
      memset(bindir,'\0',2000);
      sprintf(bindir,"%s/scaff-bin",tempc);
  
      memset(file_read1,'\0',2000);
      memset(file_read2,'\0',2000);
      sprintf(file_read1,"%s/%s",tempa,argv[args+1]);
      sprintf(file_read2,"%s/%s",tempa,argv[args+2]);
  
      memset(syscmd,'\0',2000);
      sprintf(syscmd,"pwd");
      if(system(syscmd) == -1)
      {
  //     printf("System command error:\n);
      }
  
      for(i=0;i<n_r1;i++)
      {
         if(reads_tag1[i]==1)
         {
           int clen = strlen(R1_Name1[i]);
           memset(syscmd,'\0',2000);
           sprintf(syscmd,"gunzip %s",R1_Name1[i]);
           if(system(syscmd) == -1)
           {
  //          printf("System command error:\n);
           }
           strncpy(R1_Name2[i],R1_Name1[i],clen-3);
         }
         else
           strcpy(R1_Name2[i],R1_Name1[i]);
  
      }
   
      for(i=0;i<n_r2;i++)
      {
         if(reads_tag2[i]==1)
         {
           int clen = strlen(R2_Name1[i]);
           memset(syscmd,'\0',2000);
           sprintf(syscmd,"gunzip %s",R2_Name1[i]);
           if(system(syscmd) == -1)
           {
  //          printf("System command error:\n);
           }
           strncpy(R2_Name2[i],R2_Name1[i],clen-3);
         }
         else
           strcpy(R2_Name2[i],R2_Name1[i]);
         memset(syscmd,'\0',2000);
         sprintf(syscmd,"%s/scaff_BC-reads-1 %s %s.RC1 %s.name > try.out",bindir,R1_Name2[i],R1_Name2[i],R1_Name2[i]);
         if(system(syscmd) == -1)
         {
  //       printf("System command error:\n);
         }
         memset(syscmd,'\0',2000);
         sprintf(syscmd,"%s/scaff_BC-reads-2 %s.name %s %s.RC2 > try.out",bindir,R1_Name2[i],R2_Name2[i],R2_Name2[i]);
         if(system(syscmd) == -1)
         {
  //       printf("System command error:\n);
         }
      } 
      memset(readname,'\0',2000);
      memset(syscmd,'\0',2000);
      for(i=0;i<n_r1;i++)
      {
         int clen = strlen(R1_Name2[i]);
         strncat(readname,R1_Name2[i],clen);
         strncat(readname,".RC1 ",5);  
      }
  
      printf("cat read files %s\n",readname);
      sprintf(syscmd,"cat %s > 10X-reads_1.fastq",readname);
      if(system(syscmd) == -1)
      {
  //    printf("System command error:\n);
      }
  
      memset(readname,'\0',2000);
      memset(syscmd,'\0',2000);
      for(i=0;i<n_r2;i++)
      {
         int clen = strlen(R2_Name2[i]);
         strncat(readname,R2_Name2[i],clen);
         strncat(readname,".RC2 ",5);  
      }
  
      printf("cat read files %s\n",readname);
      sprintf(syscmd,"cat %s > 10X-reads_2.fastq",readname);
      if(system(syscmd) == -1)
      {
  //    printf("System command error:\n);
      }
  
      memset(syscmd,'\0',2000);
      sprintf(syscmd,"mv 10X-reads_1.fastq %s",file_read1);
      if(system(syscmd) == -1)
      {
  //      printf("System command error:\n);
      }
  
      memset(syscmd,'\0',2000);
      sprintf(syscmd,"mv 10X-reads_2.fastq %s",file_read2);
      if(system(syscmd) == -1)
      {
  //      printf("System command error:\n);
      }
  
      if(n_debug == 0)
      {
        memset(syscmd,'\0',2000);
        sprintf(syscmd,"rm -rf * > /dev/null");
        if(system(syscmd) == -1)
        {
  //      printf("System command error:\n);
        }
  
        if(chdir(tmpdir) == -1)
        {
  //      printf("System command error:\n);
        }
  
        memset(syscmd,'\0',2000);
        sprintf(syscmd,"rm -rf %s > /dev/null",workdir);
        if(system(syscmd) == -1)
        {
  //      printf("System command error:\n);
        }
      }
      
      return EXIT_SUCCESS;
  
  }
  /* end of the main */
  
  
  
  /*   subroutine to sort out read pairs    */
  /* =============================== */
  void Read_Index(int nSeq, char *namefile)
  /* =============================== */
  {
       int i,j,nseq;
       int i_reads,n_reads,c_reads,insertsize=0;
       FILE *namef;
       char *ptr;
       char line[500];
       char **cmatrix(B64_long nrl,B64_long nrh,B64_long ncl,B64_long nch);
  
  
       if((namef = fopen(namefile,"r")) == NULL)
       {
         printf("ERROR main:: reads group file \n");
         exit(1);
       }
       nseq = 0;
       i_reads = 0;
       n_reads = 0;
       c_reads = 0;
       while(!feof(namef))
       {
         if(fgets(line,500,namef) == NULL)
              printf("Data input file problem!\n");
         if(feof(namef)) break;
         nseq++;
       }
       fclose(namef);
  
       nseq = 2*nseq;
       if((ctg_head = (int *)calloc(nseq,sizeof(int))) == NULL)
       {
         printf("fmate: calloc - ctg_head\n");
         exit(1);
       }
       if((ctg_list = (int *)calloc(nseq,sizeof(int))) == NULL)
       {
         printf("fmate: calloc - ctg_list\n");
         exit(1);
       }
       if((core_list = (int *)calloc(nseq,sizeof(int))) == NULL)
       {
         printf("fmate: calloc - core_list\n");
         exit(1);
       }
       if((read2contig = (int *)calloc(nseq,sizeof(int))) == NULL)
       {
         printf("fmate: calloc - read2contig\n");
         exit(1);
       }
  
       nseq = nseq*3;
       S_Name=cmatrix(0,nseq,0,Max_N_NameBase);
       if((insert_siz = (int *)calloc(nseq,sizeof(int))) == NULL)
       {
         printf("fmate: calloc - insert\n");
         exit(1);
       }
       if((insert_dev = (int *)calloc(nseq,sizeof(int))) == NULL)
       {
         printf("fmate: calloc - insert\n");
         exit(1);
       }
   
       if((namef = fopen(namefile,"r")) == NULL)
       {
         printf("ERROR main:: reads group file \n");
         exit(1);
       }
  
       j = 0;
       insertsize = 0;
       while(!feof(namef))
       {
         int nPair=0,len=0;
         char line2[500],line3[500],base[500];
  
         if(fgets(line,500,namef) == NULL)
              printf("Data input file problem!\n");
         if(feof(namef)) break;
         strcpy(line2,line);
         strcpy(line3,line);
         insertsize = 0;
         if((strncmp(line,"readnames",9))==0)
         {
           i=0;
           c_reads = 0;
           for(ptr=strtok(line2," ");ptr!=NULL;ptr=strtok((char *)NULL," "),i++)
           {
              if(i==5)
              {
                  memset(base,'\0',500);
  //                len=strlen(ptr);
  //                strncpy(base,ptr,len-1);
                  strcat(base,ptr);
                  c_reads = atoi(base);
              }
           }
  //       printf("creads: %d %d\n",c_reads,n_reads);
           if(n_group>0)
             ctg_list[n_group-1]=n_reads;
           n_group++;
           n_reads = 0;
         }
         else
         {      
           for(ptr=strtok(line," ");ptr!=NULL;ptr=strtok((char *)NULL," "),nPair++)
           {
           }
           i=0;
           for(ptr=strtok(line3," ");ptr!=NULL;ptr=strtok((char *)NULL," "),i++)
           {
              if(nPair>1)
              {
                if(i==(nPair-2))
                {
                  memset(base,'\0',500);
                  strcat(base,ptr);
                  insertsize = atoi(base);
                }
              }
           }
           i=0;
           j=0;
           for(ptr=strtok(line2," ");ptr!=NULL;ptr=strtok((char *)NULL," "),i++)
           {
              if(i==0)
              {
  
                len=strlen(ptr);
                if(nPair==1)
                  strncpy(S_Name[i_reads],ptr,len-1);
                else
                  strncpy(S_Name[i_reads],ptr,len);
  //       printf("reads: %d %d %d %d %s\n",j,i_reads,insertsize,c_reads,S_Name[i_reads]);
                i_reads++;
                j++;
              }
  //            else if(insertsize<50000)
              else if((insertsize<50000)&&(c_reads>4))
              {
                if(i==1)
                {
                  len=strlen(ptr);
                  strncpy(S_Name[i_reads],ptr,len);
                  i_reads++;
                  j++;
                }
                else if((i==2)&&(i<(nPair-2)))
                {
                  len=strlen(ptr);
                  strncpy(S_Name[i_reads],ptr,len);
                  i_reads++;
                  j++;
                }
  /*              else if(i==(nPair-2))
                {
                  memset(base,'\0',500);
                  strcat(base,ptr);
                  for(k=0;k<j;k++)
                     insert_siz[i_reads-k-1] = atoi(base);
                }
                else if(i==(nPair-1))
                {
                  memset(base,'\0',500);
                  strcat(base,ptr);
                  for(k=0;k<j;k++)
                     insert_dev[i_reads-k-1] = atoi(base);
                }   */
              }
           }
           n_reads = n_reads+j;
         }
         c_reads++;
       }
       fclose(namef);
       printf("contig: %d %d\n",n_reads,i_reads);
       ctg_list[n_group-1]=n_reads;
       ctg_head[0]=0;
       for(i=1;i<n_group;i++)
          ctg_head[i]=ctg_head[i-1]+ctg_list[i-1];
        
  }
  
  #define SWAP(a,b) temp=(a);(a)=b;(b)=temp;
  
  /*   Subroutine to sort an array arr[0,...,n-1] into ascending order while
       making the corresponding reaarangement of the array brr[0,...,n-1]
       by the use of Quicksort (Sedgwick, R. 1978, Communications o fthe ACM,
       vol. 21, pp. 847-857) also see Numerical Recipes in C                  */  
  
  /* =============================== */
  void ArraySort_Long(int n, B64_long *arr)
  /* =============================== */
  {
       int i,ir=n-1,j,k,m=0,jstack=0,NSTACK=50,istack[NSTACK];
       B64_long a,temp,MIN=7;
  
       for(;;)
       {
  /*      Insertion sort when subarray is small enough    */
          if(ir-m<MIN)
          {
            for(j=m+1;j<=ir;j++)
            {
               a=arr[j];
               for(i=j-1;i>=m;i--)
               {
                  if(arr[i]<=a) break;
                  arr[i+1]=arr[i];
               }
               arr[i+1]=a;
            }
            if(!jstack) return;
            ir=istack[jstack--];
            m=istack[jstack--];
          }
          else
          {
            k=(m+ir)>>1;
            SWAP(arr[k],arr[m+1]);
  
            if(arr[m]>arr[ir])
            {
              SWAP(arr[m],arr[ir]);
            }
  
            if(arr[m+1]>arr[ir])
            {
              SWAP(arr[m+1],arr[ir]);
            }
  
            if(arr[m]>arr[m+1])
            {
              SWAP(arr[m],arr[m+1]);
            }
  
            i=m+1;
            j=ir;
            a=arr[m+1];
            for(;;)
            {
               do i++; while (arr[i]<a);
               do j--; while (arr[j]>a);
               if(j<i) break;
               SWAP(arr[i],arr[j]);
            }
            arr[m+1]=arr[j];
            arr[j]=a;
            jstack+=2;
  
  /*        Push pointers to larger subarray on stack      */
  /*        process smaller subarray immediately           */
            if(jstack>NSTACK)
            {
               printf("Stack error: NSTACK too small\n");
               exit(0);
            }
            if(ir-i+1>=j-m)
            {
              istack[jstack]=ir;
              istack[jstack-1]=i;
              ir=j-1;
            }
            else
            {
              istack[jstack]=j-1;
              istack[jstack-1]=m;
              m=i;
            }
          }
       }
  }
  
  /*C substring function: It returns a pointer to the substring */
   
  char *substring(char *string, int position, int length) 
  {
     char *pointer;
     int c;
   
     pointer = malloc(length+1);
   
     if (pointer == NULL)
     {
        printf("Unable to allocate memory.\n");
        exit(1);
     }
   
     for (c = 0 ; c < length ; c++)
     {
        *(pointer+c) = *(string+position-1);      
        string++;   
     }
   
     *(pointer+c) = '\0';
   
     return pointer;
  }
  
  /* =============================== */
  void ArraySort_Int(int n, int *arr)
  /* =============================== */
  {
       int i,ir=n-1,j,k,m=0,jstack=0,NSTACK=50,istack[NSTACK];
       int a,temp,MIN=7;
  
       for(;;)
       {
  /*      Insertion sort when subarray is small enough    */
          if(ir-m<MIN)
          {
            for(j=m+1;j<=ir;j++)
            {
               a=arr[j];
               for(i=j-1;i>=m;i--)
               {
                  if(arr[i]<=a) break;
                  arr[i+1]=arr[i];
               }
               arr[i+1]=a;
            }
            if(!jstack) return;
            ir=istack[jstack--];
            m=istack[jstack--];
          }
          else
          {
            k=(m+ir)>>1;
            SWAP(arr[k],arr[m+1]);
  
            if(arr[m]>arr[ir])
            {
              SWAP(arr[m],arr[ir]);
            }
  
            if(arr[m+1]>arr[ir])
            {
              SWAP(arr[m+1],arr[ir]);
            }
  
            if(arr[m]>arr[m+1])
            {
              SWAP(arr[m],arr[m+1]);
            }
  
            i=m+1;
            j=ir;
            a=arr[m+1];
            for(;;)
            {
               do i++; while (arr[i]<a);
               do j--; while (arr[j]>a);
               if(j<i) break;
               SWAP(arr[i],arr[j]);
            }
            arr[m+1]=arr[j];
            arr[j]=a;
            jstack+=2;
  
  /*        Push pointers to larger subarray on stack      */
  /*        process smaller subarray immediately           */
            if(jstack>NSTACK)
            {
               printf("Stack error: NSTACK too small\n");
               exit(0);
            }
            if(ir-i+1>=j-m)
            {
              istack[jstack]=ir;
              istack[jstack-1]=i;
              ir=j-1;
            }
            else
            {
              istack[jstack]=j-1;
              istack[jstack-1]=m;
              m=i;
            }
          }
       }
  }
  
  
  /* =============================== */
  void ArraySort_Mix(int n, B64_long *arr, int *brr)
  /* =============================== */
  {
       int i,ir=n-1,j,k,m=0,jstack=0,b,NSTACK=50,istack[NSTACK];
       B64_long a,temp,MIN=7;
  
       for(;;)
       {
  /*      Insertion sort when subarray is small enough    */
          if(ir-m<MIN)
          {
            for(j=m+1;j<=ir;j++)
            {
               a=arr[j];
               b=brr[j];
               for(i=j-1;i>=m;i--)
               {
                  if(arr[i]<=a) break;
                  arr[i+1]=arr[i];
                  brr[i+1]=brr[i];
               }
               arr[i+1]=a;
               brr[i+1]=b;
            }
            if(!jstack) return;
            ir=istack[jstack--];
            m=istack[jstack--];
          }
          else
          {
            k=(m+ir)>>1;
            SWAP(arr[k],arr[m+1]);
            SWAP(brr[k],brr[m+1]);
  
            if(arr[m]>arr[ir])
            {
              SWAP(arr[m],arr[ir]);
              SWAP(brr[m],brr[ir]);
            }
  
            if(arr[m+1]>arr[ir])
            {
              SWAP(arr[m+1],arr[ir]);
              SWAP(brr[m+1],brr[ir]);
            }
  
            if(arr[m]>arr[m+1])
            {
              SWAP(arr[m],arr[m+1]);
              SWAP(brr[m],brr[m+1]);
            }
  
            i=m+1;
            j=ir;
            a=arr[m+1];
            b=brr[m+1];
            for(;;)
            {
               do i++; while (arr[i]<a);
               do j--; while (arr[j]>a);
               if(j<i) break;
               SWAP(arr[i],arr[j]);
               SWAP(brr[i],brr[j]);
            }
            arr[m+1]=arr[j];
            arr[j]=a;
            brr[m+1]=brr[j];
            brr[j]=b;
            jstack+=2;
  
  /*        Push pointers to larger subarray on stack      */
  /*        process smaller subarray immediately           */
            if(jstack>NSTACK)
            {
               printf("Stack error: NSTACK too small\n");
               exit(0);
            }
            if(ir-i+1>=j-m)
            {
              istack[jstack]=ir;
              istack[jstack-1]=i;
              ir=j-1;
            }
            else
            {
              istack[jstack]=j-1;
              istack[jstack-1]=m;
              m=i;
            }
          }
       }
  }
  
  /* =============================== */
  void ArraySort_Int2(int n, int *arr, int *brr)
  /* =============================== */
  {
       int i,ir=n-1,j,k,m=0,jstack=0,b,NSTACK=50,istack[NSTACK];
       int a,temp,MIN=7;
  
       for(;;)
       {
  /*      Insertion sort when subarray is small enough    */
          if(ir-m<MIN)
          {
            for(j=m+1;j<=ir;j++)
            {
               a=arr[j];
               b=brr[j];
               for(i=j-1;i>=m;i--)
               {
                  if(arr[i]<=a) break;
                  arr[i+1]=arr[i];
                  brr[i+1]=brr[i];
               }
               arr[i+1]=a;
               brr[i+1]=b;
            }
            if(!jstack) return;
            ir=istack[jstack--];
            m=istack[jstack--];
          }
          else
          {
            k=(m+ir)>>1;
            SWAP(arr[k],arr[m+1]);
            SWAP(brr[k],brr[m+1]);
  
            if(arr[m]>arr[ir])
            {
              SWAP(arr[m],arr[ir]);
              SWAP(brr[m],brr[ir]);
            }
  
            if(arr[m+1]>arr[ir])
            {
              SWAP(arr[m+1],arr[ir]);
              SWAP(brr[m+1],brr[ir]);
            }
  
            if(arr[m]>arr[m+1])
            {
              SWAP(arr[m],arr[m+1]);
              SWAP(brr[m],brr[m+1]);
            }
  
            i=m+1;
            j=ir;
            a=arr[m+1];
            b=brr[m+1];
            for(;;)
            {
               do i++; while (arr[i]<a);
               do j--; while (arr[j]>a);
               if(j<i) break;
               SWAP(arr[i],arr[j]);
               SWAP(brr[i],brr[j]);
            }
            arr[m+1]=arr[j];
            arr[j]=a;
            brr[m+1]=brr[j];
            brr[j]=b;
            jstack+=2;
  
  /*        Push pointers to larger subarray on stack      */
  /*        process smaller subarray immediately           */
            if(jstack>NSTACK)
            {
               printf("Stack error: NSTACK too small\n");
               exit(0);
            }
            if(ir-i+1>=j-m)
            {
              istack[jstack]=ir;
              istack[jstack-1]=i;
              ir=j-1;
            }
            else
            {
              istack[jstack]=j-1;
              istack[jstack-1]=m;
              m=i;
            }
          }
       }
  }
  
  /*   function to sort an array into a decreasing order:  a>b>c>....    */  
  /* =============================== */
  void ArraySort2_Int2(int n, int *arr, int *brr)
  /* =============================== */
  {
       int i,ir=n-1,j,k,m=0,jstack=0,b,NSTACK=50,istack[NSTACK];
       int a,temp,MIN=7;
  
       for(;;)
       {
  /*      Insertion sort when subarray is small enough    */
          if(ir-m<MIN)
          {
            for(j=m+1;j<=ir;j++)
            {
               a=arr[j];
               b=brr[j];
               for(i=j-1;i>=m;i--)
               {
                  if(arr[i]>=a) break;
                  arr[i+1]=arr[i];
                  brr[i+1]=brr[i];
               }
               arr[i+1]=a;
               brr[i+1]=b;
            }
            if(!jstack) return;
            ir=istack[jstack--];
            m=istack[jstack--];
          }
          else
          {
            k=(m+ir)>>1;
            SWAP(arr[k],arr[m+1]);
            SWAP(brr[k],brr[m+1]);
  
            if(arr[m]<arr[ir])
            {
              SWAP(arr[m],arr[ir]);
              SWAP(brr[m],brr[ir]);
            }
  
            if(arr[m+1]<arr[ir])
            {
              SWAP(arr[m+1],arr[ir]);
              SWAP(brr[m+1],brr[ir]);
            }
  
            if(arr[m]<arr[m+1])
            {
              SWAP(arr[m],arr[m+1]);
              SWAP(brr[m],brr[m+1]);
            }
  
            i=m+1;
            j=ir;
            a=arr[m+1];
            b=brr[m+1];
            for(;;)
            {
               do i++; while (arr[i]>a);
               do j--; while (arr[j]<a);
               if(j<i) break;
               SWAP(arr[i],arr[j]);
               SWAP(brr[i],brr[j]);
            }
            arr[m+1]=arr[j];
            arr[j]=a;
            brr[m+1]=brr[j];
            brr[j]=b;
            jstack+=2;
  
  /*        Push pointers to larger subarray on stack      */
  /*        process smaller subarray immediately           */
            if(jstack>NSTACK)
            {
               printf("Stack error: NSTACK too small\n");
               exit(0);
            }
            if(ir-i+1>=j-m)
            {
              istack[jstack]=ir;
              istack[jstack-1]=i;
              ir=j-1;
            }
            else
            {
              istack[jstack]=j-1;
              istack[jstack-1]=m;
              m=i;
            }
          }
       }
  }
  
  /* =============================== */
  void ArraySort_Mix3(int n, B64_long *arr, int *brr, int *crr)
  /* =============================== */
  {
       int i,ir=n-1,j,k,m=0,jstack=0,b,c,NSTACK=50,istack[NSTACK];
       B64_long a,temp,MIN=7;
  
       for(;;)
       {
  /*      Insertion sort when subarray is small enough    */
          if(ir-m<MIN)
          {
            for(j=m+1;j<=ir;j++)
            {
               a=arr[j];
               b=brr[j];
               c=crr[j];
               for(i=j-1;i>=m;i--)
               {
                  if(arr[i]<=a) break;
                  arr[i+1]=arr[i];
                  brr[i+1]=brr[i];
                  crr[i+1]=crr[i];
               }
               arr[i+1]=a;
               brr[i+1]=b;
               crr[i+1]=c;
            }
            if(!jstack) return;
            ir=istack[jstack--];
            m=istack[jstack--];
          }
          else
          {
            k=(m+ir)>>1;
            SWAP(arr[k],arr[m+1]);
            SWAP(brr[k],brr[m+1]);
            SWAP(crr[k],crr[m+1]);
  
            if(arr[m]>arr[ir])
            {
              SWAP(arr[m],arr[ir]);
              SWAP(brr[m],brr[ir]);
              SWAP(crr[m],crr[ir]);
            }
  
            if(arr[m+1]>arr[ir])
            {
              SWAP(arr[m+1],arr[ir]);
              SWAP(brr[m+1],brr[ir]);
              SWAP(crr[m+1],crr[ir]);
            }
  
            if(arr[m]>arr[m+1])
            {
              SWAP(arr[m],arr[m+1]);
              SWAP(brr[m],brr[m+1]);
              SWAP(crr[m],crr[m+1]);
            }
  
            i=m+1;
            j=ir;
            a=arr[m+1];
            b=brr[m+1];
            c=crr[m+1];
            for(;;)
            {
               do i++; while (arr[i]<a);
               do j--; while (arr[j]>a);
               if(j<i) break;
               SWAP(arr[i],arr[j]);
               SWAP(brr[i],brr[j]);
               SWAP(crr[i],crr[j]);
            }
            arr[m+1]=arr[j];
            arr[j]=a;
            brr[m+1]=brr[j];
            brr[j]=b;
            crr[m+1]=crr[j];
            crr[j]=c;
            jstack+=2;
  
  /*        Push pointers to larger subarray on stack      */
  /*        process smaller subarray immediately           */
            if(jstack>NSTACK)
            {
               printf("Stack error: NSTACK too small\n");
               exit(0);
            }
            if(ir-i+1>=j-m)
            {
              istack[jstack]=ir;
              istack[jstack-1]=i;
              ir=j-1;
            }
            else
            {
              istack[jstack]=j-1;
              istack[jstack-1]=m;
              m=i;
            }
          }
       }
  }
  
  
  /*   to swap the string arrays           */
  /* ============================================= */
  void s_swap(char **Pair_Name, int i, int j)
  /* ============================================= */
  {
       char temp[Max_N_NameBase];
  
       strcpy(temp,Pair_Name[i]);
       strcpy(Pair_Name[i],Pair_Name[j]);
       strcpy(Pair_Name[j],temp);
  }
  
  
  /*   to sort the string array in order          */
  /* ============================================= */
  void ArraySort_String(int n, char **Pair_Name, int *brr)
  /* ============================================= */
  {
       int i,ir=n-1,j,k,m=0,jstack=0,b,NSTACK=50,istack[NSTACK];
       int temp,MIN=7;
       char p[Max_N_NameBase];
  
       for(;;)
       {
  /*      Insertion sort when subarray is small enough    */
          if(ir-m<MIN)
          {
            for(j=m+1;j<=ir;j++)
            {
               strcpy(p,Pair_Name[j]);
               b=brr[j];
               for(i=j-1;i>=m;i--)
               {
                  if(strcmp(Pair_Name[i],p)<=0) break;
                  strcpy(Pair_Name[i+1],Pair_Name[i]);
                  brr[i+1]=brr[i];
               }
               strcpy(Pair_Name[i+1],p);
               brr[i+1]=b;
            }
            if(!jstack) return;
            ir=istack[jstack--];
            m=istack[jstack--];
          }
          else
          {
            k=(m+ir)>>1;
            s_swap(Pair_Name,k,m+1);
            SWAP(brr[k],brr[m+1]);
  
            if(strcmp(Pair_Name[m],Pair_Name[ir])>0)
            {
              s_swap(Pair_Name,m,ir);
              SWAP(brr[m],brr[ir]);
            }
  
            if(strcmp(Pair_Name[m+1],Pair_Name[ir])>0)
            {
              s_swap(Pair_Name,m+1,ir);
              SWAP(brr[m+1],brr[ir]);
            }
  
            if(strcmp(Pair_Name[m],Pair_Name[m+1])>0)
            {
              s_swap(Pair_Name,m,m+1);
              SWAP(brr[m],brr[m+1]);
            }
  
            i=m+1;
            j=ir;
            strcpy(p,Pair_Name[m+1]);
            b=brr[m+1];
            for(;;)
            {
               do i++; while (strcmp(Pair_Name[i],p)<0);
               do j--; while (strcmp(Pair_Name[j],p)>0);
               if(j<i) break;
               s_swap(Pair_Name,i,j);
               SWAP(brr[i],brr[j]);
            }
            strcpy(Pair_Name[m+1],Pair_Name[j]);
            strcpy(Pair_Name[j],p);
            brr[m+1]=brr[j];
            brr[j]=b;
            jstack+=2;
  
  /*        Push pointers to larger subarray on stack      */
  /*        process smaller subarray immediately           */
            if(jstack>NSTACK)
            {
               printf("Stack error: NSTACK too small\n");
               exit(0);
            }
            if(ir-i+1>=j-m)
            {
              istack[jstack]=ir;
              istack[jstack-1]=i;
              ir=j-1;
            }
            else
            {
              istack[jstack]=j-1;
              istack[jstack-1]=m;
              m=i;
            }
          }
       }
  }
  
  /*   to swap the string arrays           */
  /* ============================================= */
  void s_swap2(char **Pair_Name, int i, int j)
  /* ============================================= */
  {
       char temp[Max_N_NameBase2];
  
       strcpy(temp,Pair_Name[i]);
       strcpy(Pair_Name[i],Pair_Name[j]);
       strcpy(Pair_Name[j],temp);
  }
  
  
  /*   to sort the string array in order          */
  /* ============================================= */
  void ArraySort_String2(int n, char **Pair_Name, int *brr)
  /* ============================================= */
  {
       int i,ir=n-1,j,k,m=0,jstack=0,b,NSTACK=50,istack[NSTACK];
       int temp,MIN=7;
       char p[Max_N_NameBase2];
  
       for(;;)
       {
  /*      Insertion sort when subarray is small enough    */
          if(ir-m<MIN)
          {
            for(j=m+1;j<=ir;j++)
            {
               strcpy(p,Pair_Name[j]);
               b=brr[j];
               for(i=j-1;i>=m;i--)
               {
                  if(strcmp(Pair_Name[i],p)<=0) break;
                  strcpy(Pair_Name[i+1],Pair_Name[i]);
                  brr[i+1]=brr[i];
               }
               strcpy(Pair_Name[i+1],p);
               brr[i+1]=b;
            }
            if(!jstack) return;
            ir=istack[jstack--];
            m=istack[jstack--];
          }
          else
          {
            k=(m+ir)>>1;
            s_swap2(Pair_Name,k,m+1);
            SWAP(brr[k],brr[m+1]);
  
            if(strcmp(Pair_Name[m],Pair_Name[ir])>0)
            {
              s_swap2(Pair_Name,m,ir);
              SWAP(brr[m],brr[ir]);
            }
  
            if(strcmp(Pair_Name[m+1],Pair_Name[ir])>0)
            {
              s_swap2(Pair_Name,m+1,ir);
              SWAP(brr[m+1],brr[ir]);
            }
  
            if(strcmp(Pair_Name[m],Pair_Name[m+1])>0)
            {
              s_swap2(Pair_Name,m,m+1);
              SWAP(brr[m],brr[m+1]);
            }
  
            i=m+1;
            j=ir;
            strcpy(p,Pair_Name[m+1]);
            b=brr[m+1];
            for(;;)
            {
               do i++; while (strcmp(Pair_Name[i],p)<0);
               do j--; while (strcmp(Pair_Name[j],p)>0);
               if(j<i) break;
               s_swap2(Pair_Name,i,j);
               SWAP(brr[i],brr[j]);
            }
            strcpy(Pair_Name[m+1],Pair_Name[j]);
            strcpy(Pair_Name[j],p);
            brr[m+1]=brr[j];
            brr[j]=b;
            jstack+=2;
  
  /*        Push pointers to larger subarray on stack      */
  /*        process smaller subarray immediately           */
            if(jstack>NSTACK)
            {
               printf("Stack error: NSTACK too small\n");
               exit(0);
            }
            if(ir-i+1>=j-m)
            {
              istack[jstack]=ir;
              istack[jstack-1]=i;
              ir=j-1;
            }
            else
            {
              istack[jstack]=j-1;
              istack[jstack-1]=m;
              m=i;
            }
          }
       }
  }
  
  
  /* creat an int matrix with subscript ange m[nrl...nrh][ncl...nch]  */
  int     **imatrix(B64_long nrl,B64_long nrh,B64_long ncl,B64_long nch)
  {
          B64_long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
          int  **m;
  
          /* allocate pointers to rows        */
          if((m=(int **)calloc(nrow,sizeof(int*)))==NULL)
          {
             printf("error imatrix: calloc error No. 1 \n");
             return(NULL);
          }
          m+=0;
          m-=nrl;
  
          /* allocate rows and set pointers to them        */
          if((m[nrl]=(int *)calloc(nrow*ncol,sizeof(int)))==NULL)
          {
             printf("error imatrix: calloc error No. 2 \n");
             return(NULL);
          }
          m[nrl]+=0;
          m[nrl]-=nrl;
  
          for(i=nrl+1;i<=nrh;i++)
             m[i]=m[i-1]+ncol;
          /* return pointer to array of pointers to rows   */
          return m;
  }
  
  /* creat char matrix with subscript ange cm[nrl...nrh][ncl...nch]  */
  char    **cmatrix(B64_long nrl,B64_long nrh,B64_long ncl,B64_long nch)
  {
          B64_long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
          char **cm;
  
          /* allocate pointers to rows        */
          if((cm=(char **)calloc(nrow,sizeof(char*)))==NULL)
          {
             printf("error cmatrix: calloc error No. 1 \n");
             return(NULL);
          }
          cm+=0;
          cm-=nrl;
  
          /* allocate rows and set pointers to them        */
          if((cm[nrl]=(char *)calloc(nrow*ncol,sizeof(char)))==NULL)
          {
             printf("error cmatrix: calloc error No. 2 \n");
             return(NULL);
          }
          cm[nrl]+=0;
          cm[nrl]-=nrl;
  
          for(i=nrl+1;i<=nrh;i++)
             cm[i]=cm[i-1]+ncol;
          /* return pointer to array of pointers to rows   */
          return cm;
  }
  
  