sources/multicoil/interface_precuts.c

   1 /* Code by Ethan Wolf 1995.  Graphics code by David Wilson. */
   2 /*** See dimer_trimer_like.c, function calculate_likelihood to see how */
   3 /*** to compute likelihood line number before call DimerTrimerScore()  */
   4
   5 #define XCOORD_SCORED_BY 430
   6 #define XCOORD_SEQSCORES 0
   7 #define YCOORD_SEQSCORES 55 /* Each char takes "15" so put us below printing */
   8 #define WIDTH_ADD_ON 30  /** Empty space after seq in main window. **/
   9 #define xmin  120  /*  Want this far enough over so have room for seqscoreb */
  10                    /** seqscoreb is 108 wide.  xmin used to be 60 **/
  11
  12 #include "sc.h"
  13 #include <stdio.h>
  14 #include "sc2seq.h"   /* This also includes scconst.h, sc.h, interface.h */
  15 #include <X11/Xlib.h>
  16 #include <X11/keysym.h>
  17 #include <math.h>
  18 #include "scio.h"
  19 #include "options.h"
  20
  21 Display *display;
  22 Window window, zoom, printb, quit, next, prev, autob, zoomb, methodb, paramb, helpb, help, seqscoreb;
  23 GC gc, gcgray, gcdash;
  24
  25
  26 #include <stdlib.h>
  27 #include <sys/types.h>
  28 #include <sys/time.h>
  29 #include <string.h>
  30 #include "config.h"
  31 #include "switches.h"
  32
  33 #include "hmm.h"
  34
  35
  36 /* FILE *debug_out; */
  37
  38
  39 /*************************** MULTICOIL LIKELIHOOD STUFF. *********************/
  40 double both_determinants[2][MAX_CLASS_NUMBER];
  41 double *determinants= &both_determinants[0][0];
  42
  43 double both_inverse_covars[2][MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
  44 double *inverse_covars= &both_inverse_covars[0][0][0][0];
  45
  46 double both_covars_submatrix[2][MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
  47
  48 double both_means_submatrix[2][MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM];
  49 double *means_submatrix= &both_means_submatrix[0][0][0];
  50
  51 double init_class_prob[MAX_CLASS_NUMBER];
  52
  53 int multi_functnum[MAX_TABLE_NUMBER];
  54 int number_score_dim=0;
  55 int number_classes=0;
  56 int number_multi_lib[MAX_TABLE_NUMBER];
  57                       /* Used to detrmine the number of different library */
  58                       /* distance values to do gaussian param over.       */
  59
  60 #ifdef TEST_VERSION
  61  FILE *ftotal_like=NULL;   /*  Extern for now.  ***/
  62  double total_gauss_like[GRID_SIZE_DIMER][GRID_SIZE_TRIMER];
  63                     /*  When ftotal_like is read from get_defaults(), the */
  64                     /*  i,j entry will hold the total gaussian like when  */
  65                     /*  assuming dimers have init_prob =                  */
  66                     /*            MIN_DIMER_PROB + i*DIMER_GRID_STEP */
  67                     /*  and trimers have init_prob =                      */
  68                     /*            MIN_TRIMER_PROB + j*TRIMER_GRID_STEP    */
  69                     /*  These constants are defined in scconst.h.          */
  70                     /*  The value is updated during "output_seq()".        */
  71
  72  /*************************** PAIRDIFF LIKELIHOOD STUFF. *********************/
  73  double differ_determinants[MAX_CLASS_NUMBER];
  74  double differ_inverse_covars[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
  75  double differ_covars_submatrix[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
  76  double differ_means_submatrix[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM];
  77
  78  int number_differ_score_dim =0;
  79  int number_differ_lib=0;
  80  int differ_window_length = 28;
  81  int which_differentiator=1;   /** Default is just the score for the window **/
  82                            /** 2 means to give each residue a pair of scores**/
  83                            /** corresponding to the max and min of the      **/
  84                            /** windows containing it.                       **/
  85  int differ_functnum=0;
  86 #endif
  87
  88 /************************** General Scoring parameters. ********************/
  89 int functnum=0, pair_functnum=0;
  90 int preproc_table=0;
  91 double log_bound = 0; /*In PRN_MODE only scores > bound output to logfile.*/
  92
  93 int WINDOW = PAIRWINDOW;
  94 int window_length[MAX_TABLE_NUMBER];
  95
  96 double SCALE0 = DEFAULT_SCALE0;
  97 double scale0s[MAX_TABLE_NUMBER];  /** Both used as bonnie's scale0 and **/
  98 double scale0p[MAX_TABLE_NUMBER];  /** in computing posterior prob. from **/
  99                                   /** prior as sigma^2(prior)=scale0*sigma^2 */
 100                                  /** where these are the standard deviation *
 101                                   /** of the n*(prior_prob) and the distib */
 102                                  /** of number of hits when sampling       */
 103
 104 int log_offset_to_use;
 105 int offset_to_use=-1;    /*  Default offset in DimerTrimerScore is all.    */
 106 int avg_max=1;   /* Default means use max coil scores in log file.  */
 107                  /** Value of 0 means do average, 1 means do max.   **/
 108                  /** 2 means do both.  */
 109 int number_tables =0;
 110
 111 int method = MultiCoil;
 112 int main_method = MultiCoil;
 113 int main_preprocessor_method = MultiCoil;
 114 int preprocessor_method = MultiCoil;
 115
 116 int main_table = 0;
 117 int table = 0;
 118
 119 /*********************** Parameters involved in Coil Scores ***********/
 120 int by_coil_or_seq =0;   /*** To determe what type of coil_out file to do. **/
 121 int weighted_avg = 0;  /*  Determines which method for determining coil */
 122                          /*  score will be used in DimerTrimerScore.      */
 123 int start_coil_at_reg_shift =0;  /* For averaging scores over coil length. */
 124 int Coil_Score =0;   /*  0 means for MultiCoil do residues score, 1 means */
 125                      /*  average scores over coils, 2 means take max coil */
 126                      /*  score.                                           */
 127
 128
 129 Sequence sequence;
 130
 131 char *methodname[] = {"MultiCoil", "Coil","NEWCOILS", "HMMCoil","ActualCoil",
 132                          "PairCoilDiff",
 133                          "PairCoilDiffAvg"};  /* "Coil" is PairCoil or */
 134
 135
 136
 137 char table_diff_stat_filename[MAXLINE]="/dev/null";
 138 char weightfile_outname[MAXLINE]="/dev/null";
 139
 140
 141
 142
 143
 144 /***/
 145  void *hmm;
 146 /***/
 147
 148 double seq_scores[MAX_CLASS_NUMBER];
 149 double *all_scores;
 150 double *scores;    /** The space these pointers points to is alocated by **/
 151                    /** using static variable arrays in ScoreSeq */
 152             /** all_scores is scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM +1]*/
 153              /** scores is scores[current_score_algor][MAXSEQLEN][POSNUM+1] **/
 154
 155 double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
 156                   /** Since preprocess_likelihood is used everywhere, just */
 157                   /** make it an external variable.                      */
 158                   /** This is given a value in ScoreSeq and sometime in  */
 159                   /** ActualCoil().                                      */
 160 double *preprocess_like;
 161                           /* This will point to the current table version */
 162                           /* of all_preprocess_like */
 163 char all_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN];
 164
 165
 166
 167 /************************ XWindow stuff.  ***********************************/
 168 int zoomptr=0;
 169 int ps_res_per_line=0;    /* 0 is a flag that it wasn't set in get_defaults */
 170                           /* so should just print it using nominal_length   */
 171                           /* as done on the screen.  If it is set to -1     */
 172                           /* then should print it all on one line.          */
 173
 174
 175 char *motd[] = {
 176
 177   "           *********************************",
 178   "           **                             **",
 179   "           **  MultiCoil (version 1.0)    **",
 180   "           **                             **",
 181   "           **                             **",
 182   "           *********************************",
 183   "",
 184   "MultiCoil is described in:",
 185   "      Ethan Wolf, Peter Kim, Bonnie Berger,",
 186   "      `MultiCoil:  A Program for Predicting Two- and
 187   "       Three-stranded Coiled Coils'.",
 188   " ",
 189   "Email inquiries to:  multicoil-help@theory.lcs.mit.edu.",
 190   ""
 191
 192 };
 193
 194 char *hmotd[] = {
 195   "help:    Clicking on help displays this message.",
 196   " ",
 197   "next:    Displays the the next sequence.",
 198   "         Same as pressing 'n' or down-arrow.",
 199   " ",
 200   "prev:    Displays the previous sequence.",
 201   "         Same as pressing 'p' or up-arrow.",
 202   " ",
 203   "auto:    Automatically steps through the sequences",
 204   "         of a file.",
 205   " ",
 206   "zoom:    Gives more information on the highest",
 207   "         scoring residue.  Same as pressing 'z'.",
 208   "         Clicking on a specific region will give",
 209   "         more information on that region.",
 210   " ",
 211   "print:   Prints the sequence.",
 212   " ",
 213 /**********
 214   "table:  Switches the table used by PairCoil and",
 215   "        by TableDiff methods, same as 't'.  Note that",
 216   "        key 'i' (input_table) changes the preprocessor's table.",
 217   " ",
 218 **********/
 219   "method:  Switches the scoring method, same as 'm'.",
 220   "         'w' swithces Which preprocessor is used. ",
 221   " ",
 222   "key 'Up' or 'Down' moves the bound up or down by .01",
 223   "key 'Return' rescores the sequence after changing the bound",
 224   " ",
 225   "key 's' or 'l':  switch between singles and pair",
 226   "                 probabilities at lib distances.",
 227   " ",
 228   "key 'o':   Calc paircoil class prob. with respect only to",
 229   "           positive oligimerization classes (ignore pdb- class)",
 230   "key '7':   MultiCoil max score is obtained by taking tuple",
 231   "           of each dimension's max score (combo register).",
 232   "key 'x':   MultiCoil computes scores for each offset, taking",
 233   "           max of those scores to be the max score (max reg).",
 234   "key 'r':   Cycle through offsets 0-7 for register (7=combo reg)",
 235   "",
 236   "key 'c':   Switch between MultiCoil register score, average coil",
 237   "           score, and max coil score",
 238   "key 'b':   For Coil score, toggle whether start new coils at",
 239   "             a register shift (break in structure).",
 240   "key 'a':   For Coil score, toggle between unweighted and weighted avg",
 241   "",
 242   "quit:    Quits the program and completes log file.",
 243   "         'q' quits without completing the log.",
 244   "         Pressing 'q' in any window also closes   ",
 245   "         that window.",
 246   ""
 247
 248 };
 249
 250
 251 ShowTitle (Window win)
 252 {
 253   int i;
 254
 255   for (i = 0; motd[i][0]; ++i)
 256     XDrawString(display, win, gc, 5, 17*i+20, motd[i], strlen(motd[i]));
 257   XDrawString(display,win,gcgray,100,17*i+20,
 258               "Press Any Key To Continue",25);
 259 }
 260
 261 ShowHelp (Window win)
 262 {
 263   int i;
 264
 265   /* Places hmotd text in win */
 266   for (i=0; hmotd[i][0]; ++i)
 267     XDrawString(display,win,gc,5,17*i+20, hmotd[i],strlen(hmotd[i]));
 268
 269 }
 270
 271
 272 static int seqnum=0, sequp=0, seqlow=1;
 273
 274 char autoadvance=0;
 275 main (int argc, char *argv[])
 276 {
 277   struct itimerval interval;
 278   fd_set fds, fdnone;
 279   FILE *fin, *ps;
 280   XEvent report;
 281   char buff[MAXLINE+MAXSEQLEN];
 282
 283
 284
 285 /* Note that fout_raw[MAX_TABLE_NUMBER+1] and fout_coils[MAX_TABLE_NUMBER+1] */
 286 /* will be files of tuples of each residue's scores according to all tables. */
 287 /* Never mind... **/
 288   char *command_line = NULL;
 289
 290   FILE *fgin=NULL,  *fout=NULL,
 291        *flog=NULL, *flog_coil_conflicts=NULL, *fout_coils=NULL;
 292   FILE *fpin[MAX_TABLE_NUMBER];
 293   FILE *pir=NULL;
 294   char pir_name[MAXLINE];
 295   char likelihoods[MAX_TABLE_NUMBER][MAXLINE];
 296            /* Files of likelihood lines.  */
 297   char printfile[500];   /*  The filename to output postscript print to. */
 298   char *print= printfile;
 299
 300   extern int number_classes;
 301
 302   char class_sc_filenames[2][MAX_CLASS_NUMBER][MAXLINE];
 303   double both_class_covars[2][MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM]
 304   double both_class_means[2][MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM];
 305                                                [MAX_NUM_SCORE_DIM];
 306   char gauss_param[2][MAXLINE];
 307                    /** 0 entry is for combo offset, 1 entry is for **/
 308                    /** max of other registers.   ***/
 309   extern int functnum;   /* What is used in scscore.c.  Takes value from */
 310                          /* pair_functnum and multi_functnum.            */
 311
 312   char lib[MAXFUNCTNUM];   /*  Use MAXFUNCTUM=7 here since are 7 distances  */
 313                            /*  in coil, so we know AT MOST 7 lib functions. */
 314   extern int pair_functnum;   /*  The number of distances in lib[].*/
 315
 316   char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM];
 317   extern int  multi_functnum[MAX_TABLE_NUMBER];
 318
 319 #ifdef TEST_VERSION
 320   char differ_class_sc_filenames[MAX_CLASS_NUMBER][MAXLINE];
 321   double differ_means[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM];
 322   double differ_covars[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
 323   char differ_gauss_param[MAXLINE];
 324   char differ_lib[MAXFUNCTNUM];   /*  Distance values for PairCoilDiff */
 325   extern int differ_functnum;
 326 #ifdef TEST_VERSION
 327
 328
 329
 330   int acnt = 1;
 331   int mode = 0;     /*  If the i-th option is set, mode gets 1 in i-th bit. */
 332   double bound =0;
 333
 334   double upper_bound;   /* In DEBUG_MODE (for printing out the predicted    */
 335                         /* heptad repeat), scores in  [bound,upper_bound]   */
 336                         /* are printed in capital letters.                  */
 337
 338   double m[MAX_TABLE_NUMBER],
 339          b[MAX_TABLE_NUMBER];    /*The likelihood line is mx+b for paircoil. */
 340   double m_single[MAX_TABLE_NUMBER],
 341          b_single[MAX_TABLE_NUMBER];  /*  Like line for singlecoil. */
 342   int i, class;
 343
 344   int table_to_remove_from = -1;
 345                               /*    Used if want to remove ALL input  seq */
 346                               /* from that table at outset.    **/
 347
 348   int old_num_dim_table[2];
 349   int new_num_dim_table[2];
 350   int old_num_dim_differ =0;
 351                           /* The gauss param file for multiple distances is */
 352                          /* contains the parameters for all distances.     */
 353
 354
 355 #ifdef TEST_VERSION
 356   double prior_freq_single[MAX_TABLE_NUMBER],
 357          prior_freq_pair[MAX_TABLE_NUMBER];
 358   int which_priors =0, which_priorp =0;
 359   int good_turing_fixed_disc =0;    /** If set to 1 do the good_turing_fixed_disc estimate of probs*/
 360   int structural_pos[POSNUM+1];  /** Holds list of which positions are **/
 361                                /** structurally important.  All others **/
 362                                /** get genbnk probabilities....       **/
 363                                /** Value -1 indicates no more entries. **/
 364 #endif
 365
 366
 367 /*** INITIALIZATION ***/
 368
 369   fprintf(stderr,"\nAbout to initialize");
 370   initialize(&number_classes, &number_multi_lib[0], window_length, scale0s,
 371              scale0p,fpin[0],
 372              &likelihoods[0][0], &pir_name[0], print, gauss_param,
 373              differ_gauss_param,
 374              &mode, prior_freq_single, prior_freq_pair,
 375              structural_pos);
 376
 377   fprintf(stderr,"\nInitialized");
 378   read_command_line(argc, argv, &command_line, multi_functnum,
 379                           &mode, number_multi_lib, multi_lib);
 380   fprintf(stderr,"\n\nRead commandline, command_line=%s",command_line);
 381
 382 /*  Get the default options in the .paircoil file!!! */
 383   get_defaults(argv[1], &mode, &log_bound, &fgin, fpin,
 384                &ftotal_like, &fout,
 385                &flog, &flog_coil_conflicts,
 386                &fout_coils, &by_coil_or_seq, likelihoods,
 387                pir_name, lib, differ_lib, multi_lib,
 388                &pair_functnum, &differ_functnum, multi_functnum,
 389                print, &main_method, &main_preprocessor_method,
 390                &main_table, &number_tables,  &offset_to_use,
 391                &avg_max, &weighted_avg,   &Coil_Score,
 392                &start_coil_at_reg_shift,
 393                table_diff_stat_filename, weightfile_outname,
 394                &ps_res_per_line, &number_classes,
 395                class_sc_filenames, differ_class_sc_filenames,
 396                gauss_param, differ_gauss_param, number_multi_lib,
 397                &number_differ_lib,
 398                init_class_prob,
 399                &table_to_remove_from,
 400                command_line, window_length, scale0s, scale0p,
 401                prior_freq_single, prior_freq_pair, &which_priors,
 402                &which_priorp, &good_turing_fixed_disc,
 403                structural_pos, &which_differentiator, &differ_window_length);
 404
 405   fprintf(stderr,"\nOut of get_def");
 406
 407
 408   bound = log_bound;
 409
 410 #ifdef TEST_VERSION
 411   compute_multivariate_num_differ_dim(&number_differ_score_dim,
 412                                number_differ_lib,
 413                                differ_functnum, &old_num_dim_differ,
 414                                which_differentiator);
 415 #endif
 416
 417
 418   compute_multivariate_num_dim(mode, &number_score_dim,
 419                                number_multi_lib,multi_functnum, number_tables,
 420                                new_num_dim_table, old_num_dim_table);
 421
 422   if (number_classes >0)  {        /*** Get gaussian param for the classes. */
 423     if (ftotal_like) {      /** Want to output total gauss like over classes*/
 424       for (i=0; i<number_classes; i++)  /* for various init_prob settings.  */
 425         init_class_prob[i]=1;   /** All 1 values is flag to just output     */
 426       init_totals(total_gauss_like); /** gauss values, not like values in   */
 427       }                              /** convert_raw_scores_to_gauss_like.   */
 428
 429
 430     for (i=0; i<2; i++)
 431       if (gauss_param[i][0] != ',') {
 432         get_gauss_param_all_classes2(class_sc_filenames[i],
 433                                      old_num_dim_table[0]+old_num_dim_table[1],
 434                                     both_class_means[i], both_class_covars[i],
 435                                    number_classes, gauss_param[i]);
 436
 437
 438 /******Note we store the submatrices in the original matrices....****/
 439         compute_submatrices2(both_class_means[i],both_means_submatrix[i],
 440                                    both_class_covars[i],number_classes,
 441                                    both_inverse_covars[i],
 442                                    both_determinants[i], multi_lib,
 443                                    old_num_dim_table[0], old_num_dim_table[1],
 444                                    new_num_dim_table[0], new_num_dim_table[1],
 445                                    mode & MULTI_TRIMER_PAIRS);
 446 /************/
 447           for (class=0; class<number_classes; class++) {
 448             fprintf(stderr,"\n\nMean Submatrix %d:\n",class);
 449             print_matrix(both_means_submatrix[i][class],1,
 450                        new_num_dim_table[0] + new_num_dim_table[1],stderr);
 451             fprintf(stderr,"\n\nInver Covar Submatrix %d:\n", class);
 452             print_matrix(both_inverse_covars[i][class],
 453                          new_num_dim_table[0] + new_num_dim_table[1],
 454                        new_num_dim_table[0] + new_num_dim_table[1],stderr);
 455
 456           }
 457 /**********/
 458         fprintf(stderr,"\nDet0 =%lf, det1= %lf, det2= %lf\n",
 459                   both_determinants[i][0],both_determinants[i][1],
 460                   both_determinants[i][2]);
 461
 462         fprintf(stderr,"\nnew_num_dim_tab0=%d, nuw_num_dim_tab1=%d, old_num_dim_tab0=%d, old_num_dim_tab1=%d,   number_score_dim=%d", new_num_dim_table[0], new_num_dim_table[1], old_num_dim_table[0], old_num_dim_table[1], number_score_dim);
 463
 464
 465
 466       }
 467
 468   }
 469   switch_gauss_param(offset_to_use, main_method);
 470
 471   fprintf(stderr, "\n Out of get_gauss_param");
 472
 473
 474
 475 #ifdef TEST_VERSION
 476 /*************************************************************************/
 477 /******************Now do gauss parameters for PairCoilDiff method.  *****/
 478 /*************************************************************************/
 479 /**** Note the "2" parameter represents that are 2 classes (dimers, trimers)***/
 480 /**** not three, since just differentiating so don't need pdb-.    ************/
 481
 482   if (differ_gauss_param[0] != ',') {
 483     char differ_lib_hack[2][MAXFUNCTNUM];
 484                      /** If are computing max and min window scores, **/
 485                      /** then need to extract both dimenion i and   **/
 486                      /** i + old_num_dim_differ/2. Treat them like they **/
 487                      /** were from another table.   ***/
 488
 489
 490     get_gauss_param_all_classes2(differ_class_sc_filenames,
 491                                  old_num_dim_differ,
 492                                  differ_means, differ_covars,
 493                                  2, differ_gauss_param);
 494
 495     if (which_differentiator ==2) {
 496       for (i=0; i<number_differ_score_dim/2; i++) {
 497         differ_lib_hack[0][i] = differ_lib[i];
 498         differ_lib_hack[1][i] = differ_lib[i];
 499       }
 500       compute_submatrices2(differ_means,differ_means_submatrix,
 501                          differ_covars,2,
 502                          differ_inverse_covars,
 503                          differ_determinants, differ_lib_hack,
 504                          old_num_dim_differ/2, old_num_dim_differ/2,
 505                          number_differ_score_dim/2, number_differ_score_dim/2,
 506                          1);
 507     }
 508     else
 509       compute_submatrices2(differ_means,differ_means_submatrix,
 510                          differ_covars,2,
 511                          differ_inverse_covars,
 512                          differ_determinants, differ_lib,
 513                          old_num_dim_differ, 0,
 514                          number_differ_score_dim, 0,
 515                          1);
 516
 517 /************/
 518           for (class=0; class<2; class++) {
 519             fprintf(stderr,"\n\nDiffer Mean Submatrix %d:\n",class);
 520             print_matrix(differ_means_submatrix[class],1,
 521                         number_differ_score_dim,stderr);
 522             fprintf(stderr,"\n\nDiffer Inver Covar Submatrix %d:\n", class);
 523             print_matrix(differ_inverse_covars[class],
 524                          number_differ_score_dim,number_differ_score_dim,
 525                          stderr);
 526
 527           }
 528 /**********/
 529
 530     fprintf(stderr,"\nDet0 =%lf, det1= %lf\n",
 531             differ_determinants[0],differ_determinants[1]);
 532
 533     fprintf(stderr,"\nnew_num_dim_differ=%d, old_num_dim_differ=%d,", number_differ_score_dim, old_num_dim_differ);
 534
 535
 536
 537   }
 538 #endif
 539
 540
 541   fprintf(stderr, "\n Out of get_gauss_param");
 542
 543
 544
 545   log_offset_to_use = offset_to_use;
 546   method = main_method;
 547   preprocessor_method = main_preprocessor_method;
 548   table = main_table;
 549
 550
 551   if ( (mode & PRN_MODE) && (flog || flog_coil_conflicts)){
 552     log_file_header(flog, flog_coil_conflicts,  mode,
 553                     argc, argv, avg_max, log_bound,
 554                     preprocessor_method, preproc_table, log_offset_to_use,
 555                     start_coil_at_reg_shift, weighted_avg, main_method,
 556                     main_table);
 557   }
 558
 559   if (mode & POS_MODE) ActualCoil_Method = available;
 560
 561
 562 /****** Now if didn't get locations from .paircoil, get from ENV variable */
 563   if (!fgin) {                                   /* Didn't get fgin in */
 564     fprintf(stderr,"\n\n **********************");
 565     fprintf(stderr,"\nDid not get genbnk location from config file.");
 566     fprintf(stderr,"\n     Will instead use uniform distrib. as underlying distrib.");
 567     fprintf(stderr,"\n**********************\n\n");
 568
 569 /****    exit(-1); ****/                               /*  get_defaults.   */
 570   }
 571   if (!fpin[0]) {                                /* Didn't get fpin in */
 572     fprintf(stderr,"\nDid not get positive data set from config file.");
 573     exit(-1);                                   /*  get_defaults.   */
 574   }
 575 /******************************************/
 576
 577
 578   PairCoilData(fgin,fpin,&number_tables, mode & PROLINE_FREE_MODE,
 579                table_to_remove_from, argv[1], mode,
 580                prior_freq_single, prior_freq_pair, which_priors,
 581                which_priorp, good_turing_fixed_disc, structural_pos);
 582                                   /* If fpin[1] != NULL this will call    */
 583                                   /* TrimerDimerDiff to set up table diff */
 584                                   /* statistics.                          */
 585                                   /* If it is NULL then table is average  */
 586                                   /* of probabilities from 1 and 2.       */
 587
 588
 589 /***/
 590   hmm = CCHMM(fgin,fpin);
 591 /***/
 592
 593
 594   if (mode & USE_LIKE_LINE) {
 595                            /** Don't want  to  compute  likelihood line if  */
 596                            /**  not using.   */
 597     functnum = pair_functnum;
 598     get_likelihood_line(likelihoods, m, b, m_single, b_single,
 599                       lib, pair_functnum,pir_name,mode,number_tables);
 600   }
 601   for (i=0; fpin[i]; i++)
 602     fprintf(stderr, "\n m[%d]= %lf, b[%d]= %lf\n",i,m[i],i,b[i]);
 603
 604
 605
 606   i=0;
 607   while (fpin[i]) {
 608     sclose(fpin[i]);
 609     i++;
 610   }
 611   if (fgin) sclose(fgin);
 612   NEWCOILSinit();
 613
 614   sequence.seqlen = 0;
 615
 616
 617   fin = sopen(argc>1?argv[1]:"-","r");
 618   fprintf(stderr, "\nfin is opened to %s",argc>1?argv[1]:"-");
 619
 620   /***  Note: If can't open xdisplay then OpenScreen will write log and */
 621  /***  out files if those options are set.                              */
 622   fprintf(stderr,"\nInitializing\n");
 623   OpenScreen(mode,log_bound,fin,flog,flog_coil_conflicts, fout_coils,
 624              fout,  m,b,m_single, b_single,
 625              lib,differ_lib,multi_lib,
 626              main_method,main_preprocessor_method, &seqnum, which_priors,
 627              which_priorp, prior_freq_single, prior_freq_pair,
 628              good_turing_fixed_disc, structural_pos);
 629   painter();
 630   /* Let user read message while initializing */
 631
 632
 633
 634   FD_ZERO(&fds);
 635   FD_ZERO(&fdnone);
 636   interval.it_value.tv_sec = 2;
 637   interval.it_value.tv_usec = 0;
 638   interval.it_interval.tv_sec = 2;
 639   interval.it_interval.tv_usec = 0;
 640
 641   nopainter();
 642   XSelectInput(display, window, KeyPressMask | ButtonPressMask | ExposureMask);
 643   while (1) {
 644     if (autoadvance) do {
 645       FD_SET(ConnectionNumber(display),&fds);
 646       if (!select(FD_SETSIZE,&fds,&fdnone,&fdnone,&interval))
 647         NextSeq(fin,lib,differ_lib,multi_lib,
 648                 m,b,m_single,b_single,
 649                 mode,bound,flog,flog_coil_conflicts,fout_coils,
 650                 fout, which_priors, which_priorp, prior_freq_single,
 651                 prior_freq_pair, good_turing_fixed_disc, structural_pos);
 652     }
 653     while ( (!XPending(display))  && (autoadvance));
 654          /* Autoadvance until button is hit or end of seq file.  */
 655     XNextEvent(display,&report);
 656     switch (report.type) {
 657     case Expose:
 658       if (report.xany.window == zoom)
 659         Zoomer();
 660       else if (report.xany.window == help)
 661         ShowHelp(help);
 662       else if (sequence.seqlen)
 663         ShowScores(scores,&sequence,NULL,bound, mode);
 664       else
 665         ShowTitle(window);
 666       break;
 667     case ButtonPress:
 668 /*      if (!sequence.seqlen) {NextSeq(fin,lib,differ_lib,multi_lib,
 669                                      m,b,m_single, b_single,
 670                                      bound,flog,flog_coil_conflicts,
 671                                      fout_coils,fout, which_priors,
 672                                      which_priorp, prior_freq_single,
 673                                      prior_freq_pair, good_turing_fixed_disc,
 674                                      structural_pos);
 675                                      break;}
 676 */
 677       if (report.xany.window == help) {XUnmapWindow(display,help); break;}
 678       if (report.xany.window == zoom) {XUnmapWindow (display, zoom); break;}
 679       if (report.xany.window == window) {
 680         /* Zoom in and show register */
 681         int tmp;
 682         tmp = XYtoPos(report.xbutton.x,report.xbutton.y);
 683         if (tmp) {
 684           ZoomPnt(); zoomptr=tmp; Zoomer(); ZoomPnt();
 685         }
 686       }
 687       if (report.xany.window == helpb)
 688         XMapWindow(display,help);
 689       if (report.xany.window == next)
 690         NextSeq(fin,lib,differ_lib,multi_lib,
 691                 m,b,m_single, b_single,
 692                 mode,bound,flog,flog_coil_conflicts,fout_coils,fout,
 693                 which_priors,which_priorp, prior_freq_single, prior_freq_pair,
 694                 good_turing_fixed_disc, structural_pos);
 695       if (report.xany.window == prev)
 696         PrevSeq(fin,lib,differ_lib,multi_lib,
 697                 m,b,m_single,b_single,
 698                 mode,bound, which_priors,which_priorp, prior_freq_single,
 699                 prior_freq_pair, good_turing_fixed_disc, structural_pos);
 700       if (report.xany.window == autob)
 701         /* toggle the automatic sequencer */
 702         { if (autoadvance ^= 1) NextSeq(fin,lib,differ_lib,multi_lib,
 703                                         m,b,m_single,b_single,
 704                                         mode,bound,flog,
 705                                         flog_coil_conflicts,fout_coils,
 706                                         fout, which_priors,which_priorp,
 707                                         prior_freq_single, prior_freq_pair,
 708                                         good_turing_fixed_disc,
 709                                         structural_pos); }
 710       else
 711         autoadvance = 0;
 712       if (report.xany.window == zoomb) {
 713         if (zoomptr) ZoomPnt();
 714         zoomptr = GoodZoom(); Zoomer(); ZoomPnt();
 715       }
 716       if (report.xany.window == printb) {
 717         XBell(display,-10);
 718         ps = sopen(print?print:defprint,"a");
 719         PSDefs(ps, sequence.seqlen);
 720         ShowScores(scores,&sequence,ps,bound,mode);
 721         sclose(ps);
 722       }
 723 /***
 724       if (report.xany.window == tableb) {
 725         table = (table +1) % number_tables;
 726         ShowSeq(lib,differ_lib,multi_lib,
 727                 m,b,m_single,b_single,mode,0,0, bound);
 728         if (zoomptr) Zoomer();
 729       }
 730 *****/
 731       if (report.xany.window == methodb) {
 732         method_cycle(&method,NUMBER_METHODS);
 733         ShowSeq(lib,differ_lib,multi_lib,
 734                 m,b,m_single,b_single,mode,1,0, bound);
 735         if (zoomptr) Zoomer();
 736       }
 737       if (report.xany.window == quit) {
 738         if (flog  || flog_coil_conflicts  || fout || ftotal_like || fout_coils)
 739           {
 740             seqnum=sequp;
 741             finish_log(mode,log_bound,fin,flog,flog_coil_conflicts,fout_coils,
 742                      fout,m,b,m_single,b_single,
 743                      lib,differ_lib,multi_lib,main_method,main_table,
 744                      main_preprocessor_method,
 745                      &seqnum, avg_max, which_priors,which_priorp,
 746                        prior_freq_single, prior_freq_pair,
 747                        good_turing_fixed_disc, structural_pos);
 748           }
 749
 750         sclose(fin);
 751         if (flog) sclose(flog); if (fout) sclose(fout);
 752         if (fout_coils) sclose(fout_coils);
 753         if (flog_coil_conflicts) sclose(flog_coil_conflicts);
 754         if (ftotal_like)  sclose(ftotal_like);
 755         exit(0); }
 756
 757       break;
 758     case KeyPress:
 759       autoadvance = 0;
 760 /*      if (!sequence.seqlen) {NextSeq(fin,lib,differ_lib,multi_lib,
 761                                      m,b,m_single, b_single,
 762                                      mode,bound,flog,flog_coil_conflicts,
 763                                      fout_coils,fout, which_priors,
 764                                      which_priorp, prior_freq_single,
 765                                      prior_freq_pair, good_turing_fixed_disc,
 766                                      structural_pos);
 767                                      break;}
 768 */
 769       if (report.xany.window == help) {XUnmapWindow(display,help); break;}
 770       switch (XLookupKeysym(&(report.xkey),0)) {
 771       case XK_Up:
 772         bound = bound + .01;
 773         if (bound > 1) bound =1;
 774         draw_bound_box(bound, paramb);
 775         break;
 776       case XK_Down:
 777         bound = bound - .01;
 778         if (bound < 0) bound =0;
 779         draw_bound_box(bound, paramb);
 780         break;
 781
 782       case XK_Return:   /** Signal that bound has been set, should rescore **/
 783   /** -1 then means just need to recompute MultiCoil coil and seq scores.    */
 784         ShowSeq(lib,differ_lib,multi_lib,
 785                 m,b,m_single,b_single,mode,-1,-1, bound);
 786         if (zoomptr) Zoomer();
 787         break;
 788
 789
 790       case XK_C:         /** Switch between register, avg, and max **/
 791       case XK_c:         /** coil score for MultiCoil.             **/
 792         Coil_Score = (Coil_Score +1) % 3;
 793         if (Coil_Score == 0)
 794           ShowSeq(lib,differ_lib,multi_lib,  /** Change to res score **/
 795                 m,b,m_single,b_single,mode,1,1, -2,-2, bound);
 796         else ShowSeq(lib,differ_lib,multi_lib, /** Change to coil score with */
 797                 m,b,m_single,b_single,mode,-1,-1, bound); /* out recomputing */
 798         if (zoomptr) Zoomer();                            /* res score.      */
 799         break;
 800
 801       case XK_O:         /** Oligomerization ratio. ***/
 802       case XK_o:
 803         mode ^= ONLY_COILED_CLASSES;  /* To XOR to complement that bit. */
 804         if (number_classes >0) {
 805           if (method == MultiCoil)
 806             type_class_score_convert(sequence,all_scores,bound,
 807                                     mode & ONLY_COILED_CLASSES,number_tables,
 808                                     number_classes);
 809           if (preprocessor_method == MultiCoil)
 810             type_class_score_convert(sequence,all_preprocess_like,bound,
 811                                     mode & ONLY_COILED_CLASSES,number_tables,
 812                                     number_classes);
 813         }
 814         ShowSeq(lib,differ_lib,multi_lib,
 815                 m,b,m_single,b_single,mode,0,2, bound);
 816         if (zoomptr) Zoomer();
 817         break;
 818
 819       case XK_Print:
 820         XBell(display,-10);
 821         ps = sopen(print?print:defprint,"w");
 822         PSDefs(ps, sequence.seqlen);
 823         ShowScores(scores,&sequence,ps,bound,mode);
 824         sclose(ps);
 825         XBell(display,-10);
 826         break;
 827       case XK_Left:
 828         if (zoomptr>1)
 829           {ZoomPnt(); --zoomptr; Zoomer(); ZoomPnt();}
 830         break;
 831       case XK_Right:
 832         if (zoomptr && zoomptr<sequence.seqlen)
 833           {ZoomPnt(); ++zoomptr; Zoomer(); ZoomPnt();}
 834         break;
 835       case XK_Z:
 836       case XK_z:
 837         if (zoomptr) ZoomPnt();
 838         zoomptr = GoodZoom(); Zoomer(); ZoomPnt();
 839         break;
 840
 841 #ifdef TEST_VERSION
 842       case XK_A:         /** Switch between unweighted avg and  **/
 843       case XK_a:         /** weighted avg coil score.             **/
 844         weighted_avg = !weighted_avg;
 845         if (Coil_Score !=0) { /** CoilScore ==0 means not doing coil score. **/
 846           ShowSeq(lib,differ_lib,multi_lib,
 847                   m,b,m_single,b_single,mode,-1,-1, bound);
 848           if (zoomptr) Zoomer();
 849         }
 850         break;
 851
 852       case XK_b:
 853       case XK_B:         /** For CoilScore, change how deal with reg shift. **/
 854         start_coil_at_reg_shift = !start_coil_at_reg_shift;
 855         if (Coil_Score !=0) { /** CoilScore ==0 means not doing coil score. **/
 856           ShowSeq(lib,differ_lib,multi_lib,
 857                   m,b,m_single,b_single,mode,-1,-1, bound);
 858           if (zoomptr) Zoomer();
 859         }
 860         break;
 861
 862
 863       case XK_R:
 864       case XK_r:
 865         offset_to_use = (offset_to_use +1) % 8;   /* Cycle through offsets */
 866 /***    switch_gauss_param(offset_to_use);  **/
 867         ShowSeq(lib,differ_lib,multi_lib,
 868                 m,b,m_single,b_single,mode,1,1, bound);
 869         if (zoomptr) Zoomer();
 870         break;
 871       case XK_x:
 872         offset_to_use =-1;   /* Do max scoring register */
 873         ShowSeq(lib,differ_lib,multi_lib,
 874                 m,b,m_single,b_single,mode,1,1, bound);
 875         if (zoomptr) Zoomer();
 876         break;
 877       case XK_7:
 878         offset_to_use =7;   /* Do MultiCoil max score as combo of max score */
 879         ShowSeq(lib,differ_lib,multi_lib,         /* in each dimension.     */
 880                 m,b,m_single,b_single,mode,1,1, bound);
 881         if (zoomptr) Zoomer();
 882         break;
 883
 884       case XK_s:        /*  Use single probabilities. */
 885       case XK_S:
 886         if (mode & PAIRCOIL_PAIRS) {
 887           mode -= PAIRCOIL_PAIRS; /** Turn it off **/
 888           ShowSeq(lib,differ_lib,multi_lib,
 889                 m,b,m_single,b_single,mode,1,1, bound);
 890           if (zoomptr) Zoomer();
 891         }
 892         break;
 893
 894       case XK_l:        /*  Use pair probabilities at lib distances */
 895       case XK_L:
 896         if (!(mode & PAIRCOIL_PAIRS)) {
 897           mode |= PAIRCOIL_PAIRS;  /* Turn it on. */
 898           ShowSeq(lib,differ_lib,multi_lib,
 899                   m,b,m_single,b_single,mode,1,1, bound);
 900           if (zoomptr) Zoomer();
 901         }
 902         break;
 903 #endif
 904
 905       case XK_h:
 906       case XK_H:
 907       case XK_question:
 908       case XK_Help:
 909         XMapWindow(display,help);
 910         break;
 911       case XK_Q:
 912       case XK_q:
 913       case XK_Cancel:
 914         if (report.xany.window == window)
 915           {sclose(fin); if (flog) sclose(flog);
 916            if (flog_coil_conflicts) sclose(flog_coil_conflicts);
 917            if (ftotal_like)  sclose(ftotal_like);
 918            if (fout_coils)
 919              sclose(fout_coils);
 920
 921            if (fout) sclose(fout);
 922            exit(0);}
 923         if (report.xany.window == zoom) {
 924           ZoomPnt();
 925           XUnmapWindow(display, zoom);
 926           zoomptr=0;
 927         }
 928         break;
 929       case XK_M:
 930       case XK_m:
 931         method_cycle(&method,NUMBER_METHODS);
 932
 933         if (method == MultiCoil)
 934           table = table % number_classes;
 935         else
 936           table = table % number_tables;
 937
 938         ShowSeq(lib,differ_lib,multi_lib,
 939                 m,b,m_single,b_single,mode,1,0, bound);
 940         if (zoomptr) Zoomer();
 941         break;
 942
 943       case XK_W:
 944       case XK_w:
 945         method_cycle(&preprocessor_method,NUMBER_PREPROCESSORS);
 946
 947         if (preprocessor_method == MultiCoil)
 948           preproc_table = preproc_table% number_classes;
 949         else
 950           preproc_table = preproc_table%number_tables;
 951
 952         ShowSeq(lib,differ_lib,multi_lib,
 953                 m,b,m_single,b_single,mode,0,1, bound);
 954         if (zoomptr) Zoomer();
 955         break;
 956
 957       case XK_i:
 958         if (preprocessor_method == MultiCoil)
 959           preproc_table = (preproc_table +1) % number_classes;
 960         else
 961           preproc_table = (preproc_table +1) %number_tables;
 962         ShowSeq(lib,differ_lib,multi_lib,
 963                 m,b,m_single,b_single,mode,0,2,  bound);
 964         if (zoomptr) Zoomer();
 965         break;
 966       case XK_T:
 967       case XK_t:
 968         if (method == MultiCoil)
 969           table = (table +1) % number_classes;
 970         else
 971           table = (table +1) % number_tables;
 972         ShowSeq(lib,differ_lib,multi_lib,
 973                 m,b,m_single,b_single,mode,0,0, bound);
 974         if (zoomptr) Zoomer();
 975         break;
 976       case XK_P:
 977       case XK_p:
 978       case XK_Prior:
 979         if (report.xany.window == window) PrevSeq(fin, lib,differ_lib,
 980                      multi_lib,m,b,m_single,b_single,mode,bound, which_priors,
 981                             which_priorp, prior_freq_single,prior_freq_pair,
 982                             good_turing_fixed_disc, structural_pos);
 983         break;
 984
 985       case XK_N:
 986       case XK_n:
 987       case XK_Next:
 988       case XK_space:
 989       case XK_Linefeed:
 990         if (report.xany.window == window)
 991           NextSeq(fin,lib,differ_lib,multi_lib,
 992                   m,b,m_single,b_single,
 993                   mode,bound,flog,flog_coil_conflicts,fout_coils,
 994                   fout, which_priors, which_priorp,
 995                   prior_freq_single,prior_freq_pair, good_turing_fixed_disc,
 996                   structural_pos);
 997         break;
 998       }
 999       break;
1000     }
1001   }
1002 }
1003
1004 static void Alarm ()
1005 {
1006   XEvent report;
1007   if (XCheckTypedEvent(display,Expose,&report)) {
1008     ShowTitle(window);
1009     XDrawString(display,window,gcgray,500,17,"(initializing)",14);
1010   }
1011 }
1012 #include <signal.h>
1013 painter ()
1014 {
1015   struct itimerval value;
1016   signal(SIGALRM,Alarm);
1017   getitimer(ITIMER_REAL,&value);
1018   value.it_value.tv_sec = 0;
1019   value.it_value.tv_usec = 500;
1020   value.it_interval.tv_sec = 0;
1021   value.it_interval.tv_usec = 500;
1022   setitimer(ITIMER_REAL,&value,NULL);
1023   Alarm();
1024 }
1025 nopainter ()
1026 {
1027   struct itimerval value;
1028   getitimer(ITIMER_REAL,&value);
1029   value.it_value.tv_sec = 0;
1030   value.it_value.tv_usec = 0;
1031   setitimer(ITIMER_REAL,&value,NULL);
1032   XDrawImageString(display,window,gcgray,500,17,"              ",14);
1033   XDrawImageString(display,window,gcgray,225,729,"press any key to continue",25);
1034 }
1035 ZoomPnt ()
1036 {
1037   int x,y;
1038   XPoint pts[3];
1039   if (!zoomptr) return;
1040   PostoXY(zoomptr,&x,&y);
1041   pts[0].x = x;  pts[0].y = y;
1042   pts[1].x = -5; pts[1].y = 5;
1043   pts[2].x = 10; pts[2].y = 0;
1044   XFillPolygon(display,window,gcdash,pts,3,Convex,CoordModePrevious);
1045   XDrawLine(display,window,gcdash,x,y+5,x,y+8);
1046 }
1047
1048
1049 /* Find most likely residue, break ties arbitrarily */
1050 GoodZoom ()
1051 {
1052   int i, zoom = -1;
1053   double like=-1;
1054   for (i=0; i<sequence.seqlen; ++i)
1055     if (scores[8*i+7] > like) {
1056       like = scores[8*i+7];
1057       zoom = i;
1058     }
1059   return zoom+1;
1060 }
1061
1062
1063 Zoomer ()
1064 { int i, offset;
1065   double *sc;
1066   char buf[32];
1067   XMapWindow(display, zoom);
1068   /* Indicate where in the sequence we are */
1069   sprintf(buf,"Position %d    ",zoomptr);
1070   XDrawImageString(display,zoom,gc,7,20,
1071                    buf,strlen(buf));
1072   for (i = -3; i<=3; ++i) {
1073     if (zoomptr+i>0 && zoomptr+i<=sequence.seqlen)
1074       buf[0] = numaa(sequence.seq[zoomptr-1+i]);
1075     else
1076       buf[0] = ' ';
1077     XDrawImageString(display,zoom,gc,159+9*i,20, buf,1);
1078   }
1079   XDrawLine(display,zoom,gc,159,22,168,22);
1080
1081   /* Indicate the probability of being in any given state */
1082   sc = &(scores[8*(zoomptr-1)]);
1083
1084   for (i=0; i<=7; ++i)
1085     if (sc[i] == sc[7]) {offset=i; break;}  /*  Used to print out maxscoring */
1086                                             /*  register as 'a' + offset.  */
1087     /* Offset 7 will be a flag that no particular offset scored max.       */
1088
1089   /* Give likelihood */
1090   sprintf(buf,"%5.1f%% chance",100*sc[7]);
1091   XDrawImageString(display,zoom,gc,220,20, buf,strlen(buf));
1092
1093 /*****
1094   fprintf(stderr,"\noffset= %d, reg = %c",all_offsets[table][zoomptr-1],
1095                     'a' + (zoomptr-1 + all_offsets[table][zoomptr-1])%POSNUM);
1096 *******/
1097
1098   if ( sc[7] != 0)
1099     if (offset != 7)
1100       sprintf(buf,"register %c",'a'+offset);
1101     else
1102       sprintf(buf,"combo reg ");
1103   else sprintf(buf,"           ");
1104   XDrawImageString(display,zoom,gc,230,40, buf,strlen(buf));
1105
1106   for (i=0; i<7; ++i) {
1107     sprintf(buf,"%c:%4.1f%%  ",'a'+i,(100*sc[i]));
1108     XDrawImageString(display,zoom,gc,15+92*(i&3),60+20*(i/4),buf,strlen(buf));
1109   }
1110   XDrawString(display,zoom,gc,10,110,"Click on viewgram or use arrow",30);
1111   XDrawString(display,zoom,gc,10,127,"keys to view other positions.",29);
1112   XDrawString(display,zoom,gc,10,144,"You may need to move this",25);
1113   XDrawString(display,zoom,gc,10,161,"window to see the viewgram.",27);
1114 }
1115
1116
1117
1118
1119
1120
1121
1122
1123 /*********The following is for storing sequences.   ***/
1124 static char seqbuff[40000];
1125
1126 static Sequence seqs[30];
1127 static int seqbufpos[30];
1128 static int seqsizes[30];
1129
1130 ReadSeq()
1131 { int size, pos, limsize;
1132 /*
1133   size = seqsizes[sequp-seqnum];
1134   pos = seqbufpos[sequp-seqnum];
1135 */
1136   sequence = seqs[sequp-seqnum];   /** Retrieve seq from memory. **/
1137
1138 /*
1139   limsize = (pos+size<=40000) ? size : 40000-pos;
1140   bcopy(&(seqbuff[pos]),sequence.code,limsize);
1141   bcopy(seqbuff,&(sequence.code[limsize]),size-limsize);
1142 */
1143 }
1144
1145
1146 StoreSeq()
1147 { int size, pos, limsize,i;
1148   char *temp, *current;
1149
1150   bcopy(seqs,&(seqs[1]),29*sizeof(Sequence));
1151   bcopy(seqbufpos,&(seqbufpos[1]),29*sizeof(int));
1152   bcopy(seqsizes,&(seqsizes[1]),29*sizeof(int));
1153   if (sequp-seqlow>=30) seqlow = sequp-29;
1154
1155 /* The 2 accounts for the end of string characters in .code and .title. */
1156   size = 2 + strlen(sequence.code) + strlen(sequence.title) +
1157     sequence.seqlen * (sequence.reg ? 2 : 1);
1158
1159   pos = seqbufpos[1] + seqsizes[1];  /* Number of characters in array */
1160
1161 /*  if (pos>=40000) pos -= 40000;*/
1162 /*  limsize = (pos+size<=40000) ? size : 40000-pos; */
1163
1164   if (pos+size >=40000) pos=0;  /* Array will be full so start at beginning. */
1165   limsize = (pos+size<=40000) ? size : 40000-pos;
1166
1167
1168 /* Save the "code+title+seq+reg" in something permanant, and wrap around
1169 at end of array. */
1170
1171 /*  David uses bcopy because wants to copy past end of string characters */
1172 /*
1173     bcopy(sequence.code,&(seqbuff[pos]),limsize);
1174     bcopy(&(sequence.code[limsize]),seqbuff,size-limsize);
1175 */
1176
1177   sequence.code=strcpy(&(seqbuff[pos]),sequence.code);
1178   sequence.title=
1179     strcpy(&(seqbuff[pos+strlen(sequence.code)+1]),sequence.title);
1180
1181   temp=sequence.seq;
1182   sequence.seq=
1183     &(seqbuff[pos+strlen(sequence.code) + strlen(sequence.title) + 2]);
1184   for(i=0; i<sequence.seqlen; i++)
1185     sequence.seq[i] = temp[i];
1186
1187
1188   if (sequence.reg) {
1189     temp=sequence.reg;
1190     sequence.reg=&(sequence.seq[sequence.seqlen]);
1191     for(i=0; i<sequence.seqlen; i++)
1192       sequence.reg[i] = temp[i];
1193   }
1194
1195
1196   while ((seqlow<sequp) && (seqbufpos[sequp-seqlow] >= pos && seqbufpos[sequp-seqlow] < pos+limsize || seqbufpos[sequp-seqlow] < size-limsize))
1197     ++seqlow;
1198   seqs[0] = sequence;
1199   seqs[0].seqlen= sequence.seqlen;
1200   seqsizes[0] = size;
1201   seqbufpos[0] = pos;
1202 }
1203
1204
1205
1206
1207
1208
1209 method_cycle(int *Method, int Number_of_Methods)
1210 {
1211   do {
1212     *Method = (*Method+1)%Number_of_Methods;
1213   }
1214   while ((*Method == MultiCoil && MultiCoil_Method == unavailable) ||
1215          (*Method == PairCoil && PairCoil_Method == unavailable) ||
1216          (*Method == NEWCOILS && NEWCOILS_Method == unavailable) ||
1217          (*Method == HMMCoil && HMMCoil_Method == unavailable) ||
1218 /****    (*Method == TableDiffRes &&
1219           TableDiffRes_Method == unavailable)  ||
1220          (*Method == TableDiffCoil &&
1221           TableDiffCoil_Method == unavailable) ||
1222 ****/
1223          (*Method == ActualCoil && ActualCoil_Method == unavailable) ||
1224          (*Method == PairCoilDiff && PairCoilDiff_Method == unavailable)||
1225          (*Method == PairCoilDiffAvg && PairCoilDiffAvg_Method == unavailable) ||
1226
1227          (*Method == ActualCoil) && (!sequence.reg) /* Don't do actual coil */
1228           );                                        /* if didn't input register */
1229
1230 }
1231
1232
1233
1234 PrevSeq (FILE *fin,   char lib[MAXFUNCTNUM], char differ_lib[MAXFUNCTNUM],
1235          char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1236          double *m, double *b,
1237          double *m_single, double *b_single,
1238          int mode, double bound, int which_priors, int which_priorp,
1239          double prior_freq_single[MAX_TABLE_NUMBER],
1240          double prior_freq_pair[MAX_TABLE_NUMBER], int good_turing_fixed_disc,
1241          int structural_pos[POSNUM+1])
1242 {
1243   if (seqnum<=seqlow) {XBell(display,0); return;}
1244   --seqnum;
1245   ReadSeq();
1246   NewSeq(mode,sequence, which_priors, which_priorp, prior_freq_single,
1247          prior_freq_pair, good_turing_fixed_disc, structural_pos);
1248   ShowSeq(lib, differ_lib, multi_lib,
1249           m, b, m_single, b_single, mode,1,1, bound);
1250 }
1251
1252
1253 NextSeq (FILE *fin,char lib[MAXFUNCTNUM], char differ_lib[MAXFUNCTNUM],
1254          char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1255          double *m, double *b, double *m_single, double *b_single,
1256          int mode,
1257          double bound, FILE *flog, FILE* flog_coil_conflicts,
1258          FILE *fout_coils,
1259          FILE *fout, int which_priors, int which_priorp,
1260          double prior_freq_single[MAX_TABLE_NUMBER],
1261          double prior_freq_pair[MAX_TABLE_NUMBER], int good_turing_fixed_disc,
1262          int structural_pos[POSNUM+1])
1263 {
1264   char title[MAXLINE],code[MAXLINE];
1265   char seq[MAXSEQLEN], reg[MAXSEQLEN];
1266   Sequence temp_seq;
1267   int first_time_seq=0;
1268
1269 /********
1270   fprintf(stderr,"\nIn next seq");
1271 *********/
1272
1273   if (seqnum < sequp) {
1274     ++seqnum;
1275     ReadSeq();
1276   }
1277   else {
1278     temp_seq= sequence;
1279     sequence.seqlen = 0;
1280     sequence.seq = seq;
1281     sequence.reg = reg;
1282     sequence.title = title;
1283     sequence.code = code;
1284
1285
1286 /*    if (  ( (mode & POS_MODE) ? !getpos(fin,&sequence) :
1287                                 !getseq2(fin,&sequence) )  )  */
1288 /*** getseq2 now gets register info too... ***/
1289
1290     if (!getseq2(fin,&sequence))
1291       {XBell(display,0); sequence= temp_seq; autoadvance=0; return;}
1292
1293
1294     ++seqnum; ++sequp;
1295     StoreSeq();
1296     first_time_seq =1;
1297   }
1298
1299   NewSeq(mode,sequence, which_priors, which_priorp,
1300          prior_freq_single, prior_freq_pair, good_turing_fixed_disc,
1301          structural_pos);
1302   if (first_time_seq) {
1303
1304     output_seq(lib,differ_lib, multi_lib,
1305                m,b,m_single, b_single,
1306                mode,log_bound,flog,flog_coil_conflicts,fout_coils,fout,
1307                avg_max, main_method, main_preprocessor_method, main_table);
1308
1309   }
1310
1311
1312   ShowSeq(lib, differ_lib,multi_lib,
1313           m, b, m_single, b_single,mode,1,1, bound);
1314
1315 }
1316
1317
1318 /** This is stuff need to always do before rescore a sequence. **/
1319 NewSeq_nonX (int mode, Sequence sequence, int which_priors, int which_priorp,
1320              double prior_freq_single[MAX_TABLE_NUMBER],
1321              double prior_freq_pair[MAX_TABLE_NUMBER],
1322              int good_turing_fixed_disc, int structural_pos[POSNUM+1])
1323 {
1324   if (mode & TST_MODE0)
1325     recalc_prob(sequence, 0, mode, which_priors, which_priorp,
1326                 prior_freq_single, prior_freq_pair, good_turing_fixed_disc,
1327                 structural_pos);
1328   if (mode & TST_MODE1)
1329     recalc_prob(sequence,1, mode, which_priors, which_priorp,
1330                 prior_freq_single, prior_freq_pair, good_turing_fixed_disc,
1331                 structural_pos);
1332 }
1333
1334 NewSeq (int mode, Sequence sequence, int which_priors, int which_priorp,
1335         double prior_freq_single[MAX_TABLE_NUMBER],
1336         double prior_freq_pair[MAX_TABLE_NUMBER], int good_turing_fixed_disc,
1337         int structural_pos[POSNUM+1])
1338 {
1339   NewSeq_nonX (mode, sequence, which_priors, which_priorp,
1340                prior_freq_single, prior_freq_pair, good_turing_fixed_disc,
1341                structural_pos);
1342
1343   XUnmapWindow(display, zoom);
1344   zoomptr=0;
1345 }
1346
1347
1348 /* The likelihood line is mx +b */
1349 void output_seq(char lib[MAXFUNCTNUM],char differ_lib[MAXFUNCTNUM],
1350                 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1351                 double *m,double *b,
1352                 double *m_single, double *b_single,
1353                 int mode,
1354                 double log_bound,FILE *flog, FILE *flog_coil_conflicts,
1355                 FILE *fout_coils,
1356                 FILE *fout,
1357                 int avg_max,
1358                 int main_method, int main_preprocessor_method,
1359                 int main_table)
1360 {
1361
1362   double maxother=0; double maxscore;
1363   extern double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
1364
1365   int i,j;
1366   int do_preproc;
1367
1368   double junk_scores[MAXSEQLEN][POSNUM+1];
1369
1370 /***
1371   fprintf(stderr,"\nIn outputseq");
1372 ****/
1373
1374   if (is_not_differentiator(main_method))
1375     do_preproc = 0;
1376   else do_preproc =1;
1377
1378   all_scores = ScoreSeq (lib, differ_lib, multi_lib,
1379                      m, b, m_single, b_single, mode,
1380                      main_table,
1381                      main_method, main_preprocessor_method,
1382                      log_offset_to_use,  &maxscore,
1383                      1, do_preproc, log_bound);
1384
1385
1386   if (ftotal_like) {
1387     add_to_total_likes(all_scores, total_gauss_like, sequence.seqlen,
1388                        ftotal_like);
1389   }
1390 /******
1391   fprintf(stderr,"\nGot all_scores");
1392 *******/
1393   scores = &all_scores[main_table* MAXSEQLEN*(POSNUM+1) ];
1394
1395      /* This is all_scores[main_table + 0*number_tables][0][0] */
1396      /* so is score on main_table using library # 0 **/
1397      /* Similarly for preproces_like. **/
1398
1399   preprocess_like = &all_preprocess_like[main_table][0][0];
1400
1401
1402   for (i=0; i<sequence.seqlen; i++)
1403     if (preprocess_like[i*(POSNUM+1) + 7] > maxother)
1404       maxother=preprocess_like[i*(POSNUM +1) + 7];
1405                    /* This is preprocess_like[i][7] */
1406
1407
1408   switch (main_method) {
1409   case STOCK_PAIR:
1410     if (mode & VER_MODE) {
1411       NEWCOILSScore(mode, sequence, &maxother,junk_scores,
1412                     log_offset_to_use);
1413                                         /* just dump */
1414                                         /*  scores in junk_scores */
1415       log_output_ver(SC2SEQ | SCSTOCK, maxscore, maxother,
1416                      maxscore,seqnum, sequence.title,
1417                      sequence.code, flog);
1418     }
1419     break;
1420
1421   case MultiCoil:
1422 /****/
1423     if (fout_coils)
1424       /*** For pos files coils, does a log file version  of coil/seq score. **/
1425       output_pos_scores(fout_coils, sequence, all_scores, log_bound, mode,
1426                         by_coil_or_seq, weighted_avg);
1427
1428
1429     if (mode & PRN_MODE)
1430
1431       if (Coil_Score && flog) {
1432 /***  Were having  trouble with this in printing out coils, since offset
1433  ***  would sometime change at ends of coils between residue offset and
1434  ***  coil score offset.  To make log file print out nice, use res. offset.
1435         char all_offs[MAX_NUM_SCORE_DIM][MAXSEQLEN];
1436         get_offsets(all_scores, sequence,all_offs,3);
1437 ***********/
1438         log_coil_output(all_scores, all_offsets,sequence, log_bound, flog,
1439                         seqnum,mode);
1440       }
1441       else if (flog)
1442         log_output2_prn(mode,maxscore, maxscore, scores, scores,
1443                         sequence.seqlen,
1444                       seqnum, sequence, log_bound, log_bound, log_bound,
1445                         flog, avg_max);
1446  /****/
1447     break;
1448
1449   case PairCoil:
1450   case NEWCOILS:
1451  /****
1452    if (mode & VER_MODE)
1453         log_output_ver(SC2SEQ , maxscore, 0,
1454                   maxscore,seqnum, sequence.title,
1455                      sequence.code, flog);
1456 ***/
1457
1458     if (mode & PRN_MODE)
1459       log_output2_prn(mode,maxscore, maxscore, scores, scores,sequence.seqlen,
1460                       seqnum, sequence, log_bound, log_bound, log_bound,
1461                       flog, avg_max);
1462
1463
1464     if (flog_coil_conflicts)
1465       coil_conflicts(mode,all_scores,all_scores,sequence, seqnum,log_bound,
1466                      log_bound,number_tables, flog_coil_conflicts,avg_max);
1467    break;
1468
1469 /*  case PairCoilDiff: */
1470   case PairCoilDiffAvg:
1471     if (mode & PRN_MODE) {
1472       log_output2_prn(mode,maxscore, maxother, scores,preprocess_like,
1473                       sequence.seqlen,
1474                       seqnum, sequence, .5, .5, log_bound,
1475                       flog, avg_max);
1476     }
1477
1478     if (flog_coil_conflicts) {
1479
1480       coil_conflicts(mode,all_scores,all_preprocess_like,sequence, seqnum,
1481                      .5, log_bound, number_tables, flog_coil_conflicts,
1482                      avg_max);
1483
1484     }
1485     break;
1486
1487   case ActualCoil:
1488     break;
1489   }
1490
1491
1492 /***************
1493   fprintf(stderr,"\n About to go into fout.");
1494 ***************/
1495
1496
1497   if ( fout)
1498     if ( (!(mode & RAW_OUT)) ||
1499         ( (main_method != MultiCoil) && (main_method != PairCoilDiff))) {
1500       if (main_method== MultiCoil)
1501         tuple_output2(sequence, mode, fout, all_scores,3, 1,log_bound,3);
1502                                   /* Output based on 3 tables, meaning */
1503                                   /* output the dimer and trimer likes */
1504                                   /* and total like. */
1505       else txt_output2(sequence, mode, fout, scores, log_bound);
1506      }
1507     else {   /* Want raw tuple_output for multi_coil and PairCoilDiff*/
1508
1509       if (main_method == MultiCoil)
1510         if (log_offset_to_use==7)   /* Combo offset just output max score  */
1511           tuple_output2(sequence, mode, fout, /** for each dimension.      */
1512                         all_scores, number_tables,
1513                         number_multi_lib,log_bound, number_score_dim);
1514         else  /** Max reg for raw score can't really be determined. So instead */
1515               /* for coils, output scores for the known correct register, and */
1516               /* for pdb- output ALL possible 7 register scores (all negs).  */
1517           tuple_output_for_max(sequence,mode,fout, all_scores, number_tables,
1518                                number_multi_lib, log_bound, number_score_dim);
1519
1520       else if (main_method == PairCoilDiff)
1521                               /* Output max score for each dimension. */
1522                    /** Had to do a hack since tuple_out2 expects number_lib  */
1523                   /** to be array of size num_tab, so pass in 1 for num_tab */
1524         tuple_output2(sequence, mode, fout, all_scores, 1,
1525                       &number_differ_lib,log_bound, number_differ_score_dim);
1526
1527     }
1528
1529 }
1530
1531
1532 int copy_offsets(char preproc_offsets[MAXSEQLEN],
1533                  Sequence sequence,
1534                  char all_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN],
1535                  int number_dimens,
1536                  double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
1537                  int offset_to_use)
1538 {
1539   int tablenum,i;
1540
1541   for (tablenum =0; tablenum < number_dimens; tablenum++)
1542     for(i=0; i<sequence.seqlen; i++)  {
1543       all_offsets[tablenum][i]=preproc_offsets[i];
1544       /*** Make sure the output score corresponds to this offset... **/
1545       /*** Also do it for backup residue scores in tablenum +3 ***/
1546       if ((offset_to_use ==7) || (offset_to_use == -1)) {
1547         all_scores[tablenum][i][7] =
1548           all_scores[tablenum][i][(i+all_offsets[tablenum][i]) %POSNUM];
1549         all_scores[tablenum+3][i][7] =
1550           all_scores[tablenum+3][i][(i+all_offsets[tablenum][i]) %POSNUM];    }
1551     }
1552 }
1553
1554 int get_offsets(double all_likes[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
1555                  Sequence sequence, char offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN],
1556                  int number_dimens)
1557 {
1558   int tablenum;
1559   int got_offset, max_offset;
1560   double max_off_score;
1561   int i,j;
1562
1563   for (tablenum =0; tablenum < number_dimens; tablenum++)
1564     for(i=0; i<sequence.seqlen; i++) {
1565       offsets[tablenum][i]=-1;
1566       got_offset=0;
1567       max_off_score=-HUGE_VAL;
1568       for (j=0; j<POSNUM; j++) {
1569         if (all_likes[tablenum][i][j] == all_likes[tablenum][i][7]) {
1570           offsets[tablenum][i] = (j-i) % POSNUM;
1571           if (offsets[tablenum][i]<0)
1572             offsets[tablenum][i] += POSNUM;
1573           got_offset=1;
1574         }
1575         if (all_likes[tablenum][i][j] > max_off_score) {
1576           max_off_score=all_likes[tablenum][i][j];
1577           max_offset=(j-i) % POSNUM;
1578           if (max_offset <0) max_offset +=POSNUM;
1579         }
1580       }
1581       if (!got_offset) offsets[tablenum][i] = max_offset;
1582
1583       /* This is a new thing so that if an offset change is not real */
1584       /* we keep the old offset. **********/
1585       if ( (i>0) && (all_likes[tablenum][i][(i+offsets[tablenum][i-1]) %POSNUM]
1586                  == all_likes[tablenum][i][(i+offsets[tablenum][i])%POSNUM]) )
1587         offsets[tablenum][i] = offsets[tablenum][i-1];
1588
1589     }
1590 }
1591
1592
1593
1594
1595
1596
1597
1598
1599 void  preprocessor_score(char lib[MAXFUNCTNUM],
1600                          char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1601                          double *m, double *b,
1602                        double *m_single, double *b_single, int mode,
1603                          int Preprocessor_Method, int Offset_to_Use,
1604             double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
1605                          Sequence sequence,
1606                         char all_preproc_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN],
1607                          double bound)
1608 {
1609   double maxscore;
1610   char offsets[MAXSEQLEN];
1611   int tablenum;
1612   char  *local_lib;
1613   int dist;
1614   int use_pairs=1;
1615   int dimension, num_multilib;
1616
1617
1618
1619   for (tablenum =0; tablenum < number_tables; tablenum++) {
1620     switch_tables(tablenum);  /* Set prob. to be for that table. */
1621     if (tablenum ==0) use_pairs =1;  /*  Use pairs for cctb */
1622     else if (tablenum ==1) use_pairs = (mode & MULTI_TRIMER_PAIRS);
1623                /* Do what config file said to do for trimers. **/
1624
1625     switch(Preprocessor_Method) {
1626     case MultiCoil:
1627       functnum = multi_functnum[tablenum];   /*  For use in scscore.c **/
1628       num_multilib= compute_num_multilib(tablenum, number_multi_lib,
1629                                          mode & MULTI_TRIMER_PAIRS);
1630       for (dist=0; dist<num_multilib; dist++) {
1631         if (number_multi_lib[tablenum]>1) {
1632           local_lib = &multi_lib[tablenum][dist];
1633         }
1634         else local_lib = multi_lib[tablenum];
1635         dimension = compute_dimension(tablenum, number_tables, dist,
1636                                       number_score_dim, number_multi_lib,
1637                                       mode & MULTI_TRIMER_PAIRS);
1638
1639         MultiCoilDimensionScore(mode, sequence, local_lib,
1640                       &maxscore,tablenum, offsets,
1641                       all_preprocess_like[dimension],
1642                       Offset_to_Use,use_pairs);
1643       }
1644       break;
1645
1646     case PairCoil:
1647       if (mode & PAIRCOIL_PAIRS) {
1648         functnum = pair_functnum;   /* For use in scscore.c **/
1649         PairCoilScore(mode, sequence, lib, m[tablenum], b[tablenum],
1650                       &maxscore,tablenum, offsets,
1651                       all_preprocess_like[tablenum],
1652                       Offset_to_Use, 0);
1653       }
1654       else    /*  Use singles probabilities... for likelihood use line 0 */
1655         SingleCoilScore(mode,sequence,m_single[tablenum],
1656                                b_single[tablenum],
1657                                &maxscore, tablenum,
1658                                all_preprocess_like[tablenum], Offset_to_Use,
1659                                mode & RAW_OUT);
1660       break;
1661
1662     case NEWCOILS:
1663       NEWCOILSScore(mode, sequence, &maxscore, all_preprocess_like[tablenum], Offset_to_Use);
1664       break;
1665
1666     case HMMCoil:
1667       HMMScore(sequence.seq, sequence.seqlen, hmm, all_preprocess_like[tablenum],
1668                tablenum, Offset_to_Use);
1669       break;
1670
1671     case ActualCoil:
1672   /* if (mode & POS_MODE)  */    /* Score the coils in the posfile based on */
1673                                /* their having 100% likelihood. */
1674       ActualCoils(sequence, all_preprocess_like[tablenum], Offset_to_Use,1);
1675       break;
1676
1677     default:
1678       fprintf(stderr,"\nBad preprocessor, using PairCoil.\n");
1679       functnum = pair_functnum;   /* For use in scscore.c **/
1680       PairCoilScore(mode, sequence, lib, m[tablenum], b[tablenum],
1681                     &maxscore,tablenum, offsets,
1682                     all_preprocess_like[tablenum],
1683                     Offset_to_Use,  0);
1684       break;
1685     }
1686   }
1687
1688   if (Preprocessor_Method !=MultiCoil)
1689     get_offsets(all_preprocess_like, sequence,all_preproc_offsets,
1690                 number_tables);
1691   else {
1692 /***** Get copy of res scores in tablenum and tablenum+3 **/
1693     switch_gauss_param(Offset_to_Use, Preprocessor_Method);
1694
1695     convert_raw_scores_to_gauss_prob_like2(all_preprocess_like,sequence.seqlen,
1696                                  number_tables,means_submatrix,inverse_covars,
1697                                  determinants, number_classes, init_class_prob,
1698                                  number_score_dim,
1699                                      mode & ONLY_COILED_CLASSES,
1700                                      bound,Offset_to_Use,  1);
1701
1702     get_offsets(all_preprocess_like, sequence,all_preproc_offsets,3);
1703
1704     if (Coil_Score)  {  /* Value of 1 means do average, 2 means do max. */
1705              /*** Need to save residue scores, since the coil scores **/
1706              /*** will replace them, but old residue scores might be needed **/
1707              /*** again in tablenum+3. ****/
1708       for (tablenum =0;  tablenum <3; tablenum++) {
1709         average_score_over_coils3(sequence, all_preprocess_like[tablenum+3],
1710                                   all_preprocess_like[tablenum],
1711                              all_preprocess_like[2+3],/** Total Coil Like. **/
1712                                 weighted_avg, offset_to_use,
1713                                 all_preproc_offsets[tablenum],
1714                                 start_coil_at_reg_shift, bound,
1715                                 Coil_Score -1);
1716       }
1717     }
1718   }
1719
1720
1721
1722
1723 }
1724
1725
1726
1727
1728 double *ScoreSeq (char lib[MAXFUNCTNUM], char differ_lib[MAXFUNCTNUM],
1729                   char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
1730                   double *m_paircoil, double *b_paircoil,
1731                   double *m_single, double *b_single,
1732                   int mode, int Table, int Method,
1733                   int Preprocessor_Method, int Offset_to_Use,
1734                   double *maxscore,
1735                   int rescore_seq, int rescore_preproc, double bound)
1736 {
1737   /** rescore_preproc == 1 means need to rescore.  If it or rescore_seq are  */
1738   /** -1 then means just need to recompute MultiCoil coil and seq scores.    */
1739   /** -2 means need to change from MultiCoil coil score back to res score.   */
1740   /**  If it is another non-zero */
1741   /** value, don't need to rescore it, but do need to rescore anything that  */
1742   /** uses it as a preproc.                                                  */
1743
1744   double max_paircoil_score;
1745   static char all_preproc_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN];
1746   extern char all_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN];
1747   extern double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
1748   static double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
1749   int i,j;
1750   int tablenum;
1751   double maxsc;
1752   char *local_lib;
1753   int dist;
1754
1755   int use_pairs=1;
1756   int num_multilib, dimension;
1757   int local_number_tables;
1758
1759   double *preproc_res_like;
1760   int scored_already =0;
1761
1762 /**************************************************************************/
1763 /********* PreprocLike is used in the 2-3 stranded differentiators.   ****/
1764
1765
1766   if ( (rescore_preproc==1) &&
1767                       (!ftotal_like)) {  /* A hack to save computation time. */
1768     preprocessor_score(lib, multi_lib,
1769                        m_paircoil, b_paircoil, m_single, b_single,
1770                        mode,  Preprocessor_Method,
1771                        Offset_to_Use,
1772                        all_preprocess_like, sequence, all_preproc_offsets,
1773                        bound);
1774   }
1775
1776 /*****************************************************************************/
1777   if (Preprocessor_Method == MultiCoil) {
1778     if ( ( rescore_preproc==-1) && Coil_Score) {
1779              /* Coil_Score value of 1 means do average, 2 means do max. */
1780              /*** Need to save residue scores, since the coil scores **/
1781              /*** will replace them, but old residue scores might be needed **/
1782              /*** again in tablenum+3. ****/
1783       for (tablenum =0;  tablenum <3; tablenum++)
1784         average_score_over_coils3(sequence, all_preprocess_like[tablenum+3],
1785                                   all_preprocess_like[tablenum],
1786                              all_preprocess_like[2+3],/** Total Coil Like. **/
1787                                 weighted_avg, offset_to_use,
1788                                 all_preproc_offsets[tablenum],
1789                                 start_coil_at_reg_shift, bound,
1790                                 Coil_Score -1);
1791       if (mode & ONLY_COILED_CLASSES)  /** Convert to olig state. **/
1792         type_class_score_convert(sequence,all_preprocess_like,bound,
1793                                   mode & ONLY_COILED_CLASSES,number_tables,
1794                                   number_classes);
1795     }
1796     else if (rescore_preproc == -2)           /** Switch back to res score.*/
1797       for (tablenum =0; tablenum <3 ; tablenum++)
1798         for (j =0; j <=POSNUM; j++)
1799           for (i=0; i<sequence.seqlen; i++)  {
1800             all_preprocess_like[tablenum][i][j] =
1801               all_preprocess_like[tablenum+3][i][j];
1802             if ((mode & ONLY_COILED_CLASSES) && (tablenum <2)&&
1803                 (all_preprocess_like[2+3][i][j]>bound))
1804               all_preprocess_like[tablenum][i][j]/=
1805                 all_preprocess_like[2+3][i][j];
1806           }
1807
1808   }
1809 /****************************************************************************/
1810
1811
1812
1813
1814
1815   if (Method == PairCoilDiff)
1816     local_number_tables = 1;
1817   else local_number_tables= number_tables;
1818
1819
1820 /*** Clean up this conditional, since really the only thing that might ***/
1821 /*** change for paircoildiffer when change preproc is the offset to use **/
1822 /*** (i.e. all_scores[table][i][7]. ***/
1823
1824   if ((rescore_seq>0) ||
1825     (!is_not_differentiator(Method) && rescore_preproc) )  {
1826     scored_already=1;
1827     for (tablenum=0; tablenum<local_number_tables; tablenum++) {
1828       switch_tables(tablenum);    /* Switch prob to right table. */
1829
1830       if (tablenum ==0) use_pairs =1;  /*  Use pairs for cctb */
1831       else if (tablenum ==1) use_pairs = (mode & MULTI_TRIMER_PAIRS);
1832                /* Do what config file said to do for trimers. **/
1833
1834       maxsc=0;  /* Initialization. */
1835
1836
1837       switch (Method) {
1838
1839       case STOCK_PAIR:
1840         functnum = pair_functnum;   /* For use in scscore.c **/
1841         PairCoilScore(mode, sequence, lib, m_paircoil[tablenum],
1842                    b_paircoil[tablenum],
1843                    &maxsc,tablenum, all_offsets[tablenum],all_scores[tablenum],
1844                    Offset_to_Use, 0);
1845         break;
1846
1847
1848       case MultiCoil:
1849         functnum = multi_functnum[tablenum];   /* For use in scscore.c **/
1850         num_multilib = compute_num_multilib(tablenum,number_multi_lib,
1851                                             mode & MULTI_TRIMER_PAIRS);
1852         for (dist=0; dist<num_multilib; dist++) {
1853           if (number_multi_lib[tablenum] >1)
1854             local_lib = &multi_lib[tablenum][dist];
1855           else local_lib = multi_lib[tablenum];
1856           dimension = compute_dimension(tablenum, number_tables, dist,
1857                                         number_score_dim, number_multi_lib,
1858                                         mode & MULTI_TRIMER_PAIRS);
1859
1860           MultiCoilDimensionScore(mode, sequence, local_lib, &maxsc,
1861                         tablenum,all_offsets[dimension],
1862                         all_scores[dimension],
1863                         Offset_to_Use,use_pairs);
1864         }
1865
1866         break;
1867
1868       case PairCoil:
1869         if (mode & PAIRCOIL_PAIRS) {
1870           functnum = pair_functnum;   /* For use in scscore.c **/
1871           PairCoilScore(mode, sequence, lib, m_paircoil[tablenum],
1872                         b_paircoil[tablenum], &maxsc,
1873                         tablenum,all_offsets[tablenum],
1874                         all_scores[tablenum],
1875                         Offset_to_Use, (mode & RAW_OUT));
1876         }
1877         else /*Use singles probabilities... for likelihood use line number 0 */
1878           SingleCoilScore(mode,sequence,m_single[tablenum],
1879                           b_single[tablenum],
1880                           &maxsc, tablenum,
1881                           all_scores[tablenum], Offset_to_Use,
1882                           (mode & RAW_OUT));
1883         break;
1884       case NEWCOILS:
1885         NEWCOILSScore(mode,sequence,&maxsc, all_scores[tablenum],
1886                       Offset_to_Use);
1887         break;
1888
1889       case HMMCoil:
1890         HMMScore(sequence.seq, sequence.seqlen, hmm, all_scores[tablenum],
1891                  tablenum, Offset_to_Use);
1892         break;
1893
1894       case PairCoilDiff:
1895         functnum = differ_functnum;   /* For use in scscore.c **/
1896         for (dimension=0; dimension<number_differ_lib; dimension++) {
1897           if (number_differ_lib >1) {
1898             local_lib = &differ_lib[dimension];
1899           }
1900           else local_lib = differ_lib;
1901
1902
1903           PairCoilDifferDimension(mode, sequence, local_lib, &maxsc,
1904                                   all_preproc_offsets[preproc_table],
1905                                   all_preprocess_like[preproc_table],
1906                                   all_scores,
1907                                   all_offsets, dimension,
1908                                   Offset_to_Use, number_differ_score_dim,
1909                                   bound);
1910
1911         }
1912 /***********************
1913           PairCoilDiffer(mode, sequence, local_lib, &maxsc,
1914                              all_preproc_offsets[tablenum],
1915                              all_preprocess_like[preproc_table],
1916                              Offset_to_Use, tablenum,
1917                              all_scores[tablenum],0,
1918                              weighted_avg,start_coil_at_reg_shift);
1919
1920 *************************/
1921         break;
1922
1923 /*************
1924       case PairCoilDiffAvg:
1925         functnum = differ_functnum;
1926         PairCoilDiffer(mode, sequence, differ_lib, &maxsc,
1927                        all_preproc_offsets[preproc_table],
1928                        all_preprocess_like[preproc_table],
1929                        Offset_to_Use, tablenum,
1930                        all_scores[tablenum],
1931                        1, weighted_avg,start_coil_at_reg_shift);
1932         break;
1933 ****************/
1934
1935       case ActualCoil:
1936         ActualCoils(sequence, all_scores[tablenum], Offset_to_Use,0);
1937         break;
1938       }
1939       if (tablenum == Table) *maxscore = maxsc;
1940
1941     }
1942
1943
1944
1945
1946
1947 /*********************************************************************/
1948 /*********************************************************************/
1949 /********* This ends the scoring along all and all dimensions ********/
1950 /********* Now convert the raw scores that haven't been converted ****/
1951 /********* yet into likelihoods (i.e. the multidimensional scores.****/
1952
1953
1954 /***** Note:  Can use same "init_class_prob" here, because since set number **/
1955 /***  of classes to 2, the 3rd class (Pdb-) will be ignored, so just get  ****/
1956 /***  ratio of two to three stranded likelihoods, BOTH weighted by their  ****/
1957 /***  init probaabiities.       ********/
1958
1959     if ( (Method != MultiCoil)  && (Method !=PairCoilDiff))
1960       get_offsets(all_scores, sequence,all_offsets,number_tables);
1961     else if (!(mode & RAW_OUT)) {
1962
1963       switch_gauss_param(Offset_to_Use, Method);
1964
1965
1966 /***********Do conversions on PairCoilDiff scores, using preprocessor. **/
1967 /*********** Note for differentiator there are only 2 possible classes **/
1968 /*********** (dimer and trimer). **********/
1969       if (Method == PairCoilDiff) {
1970         convert_raw_scores_to_gauss_prob_like2(all_scores,sequence.seqlen,
1971                           number_tables, differ_means_submatrix,
1972                           differ_inverse_covars,differ_determinants, 2,
1973                           init_class_prob, number_differ_score_dim,
1974                            0, bound, Offset_to_Use,0);
1975         if (Preprocessor_Method == MultiCoil)
1976           preproc_res_like = &all_preprocess_like[preproc_table+3][0][0];
1977         else preproc_res_like = &all_preprocess_like[preproc_table][0][0];
1978
1979
1980 /*  Only use window scores that lie completely in coil */
1981         if ( (which_differentiator == 3) || (which_differentiator ==4))
1982           zero_out_non_coils(all_scores, all_preprocess_like[preproc_table],
1983                              sequence.seqlen, bound);
1984
1985         if (which_differentiator == 1)
1986           zero_out_bad_windows(differ_window_length,
1987                                all_scores,all_preprocess_like[preproc_table],
1988                                sequence.seqlen,bound);
1989
1990         copy_offsets(all_preproc_offsets[preproc_table],sequence,
1991                       all_offsets,2, all_scores, Offset_to_Use);
1992
1993
1994         if (Coil_Score) /* Value of 1 means do average, 2 means do max. */
1995           for (tablenum=0; tablenum<2; tablenum++) {
1996             average_score_over_coils3(sequence, all_scores[tablenum+3],
1997                                     all_scores[tablenum],
1998                                     preproc_res_like,
1999                                     /** Total Coil Like. **/
2000                                 weighted_avg, offset_to_use,
2001                                   all_preproc_offsets[preproc_table],
2002                                 start_coil_at_reg_shift, bound,
2003                                 Coil_Score -1);
2004         }
2005       }
2006
2007 /**********Do conversions on MultiCoil dimension scores. ******/
2008
2009       else if (Method == MultiCoil) {
2010         convert_raw_scores_to_gauss_prob_like2(all_scores,sequence.seqlen,
2011                              number_tables, means_submatrix,inverse_covars,
2012                              determinants, number_classes,init_class_prob,
2013                              number_score_dim,
2014                              mode & ONLY_COILED_CLASSES, bound,
2015                              Offset_to_Use,1);
2016         get_offsets(all_scores, sequence,all_offsets,3);
2017
2018
2019
2020         if (Coil_Score) /* Value of 1 means do average, 2 means do max. */
2021           for (tablenum=0; tablenum<3; tablenum++)
2022             average_score_over_coils3(sequence, all_scores[tablenum+3],
2023                                     all_scores[tablenum],
2024                                 all_scores[2+3],  /** Total Coil Like. **/
2025                                 weighted_avg, offset_to_use,
2026                                   all_offsets[tablenum],
2027                                 start_coil_at_reg_shift, bound,
2028                                 Coil_Score -1);
2029       }
2030     }
2031
2032
2033   }   /* Ends if rescore_seq */
2034
2035
2036 /*****************************************************************************/
2037   if ( (Method == MultiCoil) || (Method == PairCoilDiff)){
2038
2039     if (Method == MultiCoil)
2040       if ( (rescore_seq ==1 ) ||  /** Signals rescored MultiCoil **/
2041           (rescore_seq==-1)) /*Signals bound changed, need to redo seq_score */
2042         get_seq_scores(seq_scores, sequence, all_scores[0+3],
2043                                         all_scores[1+3],bound, weighted_avg);
2044
2045     if (rescore_seq==-1 && !scored_already && (Method==PairCoilDiff)) {
2046       if ( (which_differentiator == 3) || (which_differentiator ==4))
2047         zero_out_non_coils(all_scores, all_preprocess_like[preproc_table],
2048                              sequence.seqlen, bound);
2049       if (which_differentiator == 1) /*  Only use window scores that lie */
2050         zero_out_bad_windows(differ_window_length, /* completely in coil */
2051                                all_scores, all_preprocess_like[preproc_table],
2052                                sequence.seqlen,bound);
2053     }
2054
2055     if ( (rescore_seq==-1) && Coil_Score && !scored_already) {
2056       if (Method == MultiCoil) {
2057         for (tablenum=0; tablenum<3; tablenum++) {
2058           average_score_over_coils3(sequence, all_scores[tablenum+3],
2059                                     all_scores[tablenum],
2060                                 all_scores[2+3],  /** Total Coil Like. **/
2061                                 weighted_avg, offset_to_use,
2062                                   all_offsets[tablenum],
2063                                 start_coil_at_reg_shift, bound,
2064                                 Coil_Score -1);
2065           if (mode & ONLY_COILED_CLASSES)   /** Convert to olig. state. **/
2066             type_class_score_convert(sequence,all_scores,bound,
2067                                      mode & ONLY_COILED_CLASSES, number_tables,
2068                                      number_classes);
2069         }
2070       }
2071       else if (Method == PairCoilDiff)  {
2072         if (Preprocessor_Method == MultiCoil)
2073           preproc_res_like = &all_preprocess_like[preproc_table+3][0][0];
2074         else preproc_res_like = &all_preprocess_like[preproc_table][0][0];
2075
2076
2077         for (tablenum=0; tablenum<2; tablenum++)
2078           average_score_over_coils3(sequence, all_scores[tablenum+3],
2079                                     all_scores[tablenum],
2080                                     preproc_res_like,
2081                                     /** Total Coil Like. **/
2082                                     weighted_avg, offset_to_use,
2083                                     all_offsets[tablenum],
2084                                     start_coil_at_reg_shift, bound,
2085                                     Coil_Score -1);
2086       }
2087     }
2088
2089
2090     else if (rescore_seq == -2)           /** Switch back to reg score.*/
2091       for (tablenum =0; tablenum <3 ; tablenum++)
2092         for (j =0; j <=POSNUM; j++)
2093           for (i=0; i<sequence.seqlen; i++) {
2094             all_scores[tablenum][i][j] =
2095               all_scores[tablenum+3][i][j];
2096             if ((Method == MultiCoil) &&
2097                 (mode & ONLY_COILED_CLASSES) && (tablenum <2) &&
2098                 (all_scores[2+3][i][j]>bound))
2099               all_scores[tablenum][i][j]/= all_scores[2+3][i][j];
2100           }
2101   }
2102 /****************************************************************************/
2103
2104
2105
2106
2107   return (double*)all_scores;
2108
2109 }
2110
2111
2112
2113
2114 ShowSeq(char lib[MAXFUNCTNUM], char differ_lib[MAXFUNCTNUM],
2115         char  multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
2116         double *m, double *b,
2117          double *m_single, double *b_single,
2118          int mode, int rescore_sequence, int rescore_preproc, double bound)
2119 {
2120   double maxscore;
2121
2122 /*  fprintf(stderr,"\nIn show seq"); */
2123
2124
2125   if ( (rescore_sequence) || (rescore_preproc))
2126     all_scores = ScoreSeq(lib, differ_lib, multi_lib, m, b,m_single, b_single,
2127                     mode, table,
2128                     method, preprocessor_method,
2129                     offset_to_use, &maxscore,
2130                     rescore_sequence, rescore_preproc, bound);
2131
2132
2133 /*  fprintf(stderr,"\n Got all_scores in showseq"); */
2134
2135   scores = &all_scores[table * MAXSEQLEN*(POSNUM+1)];
2136      /* This is all_scores[main_table][0][0] */
2137   preprocess_like = &all_preprocess_like[preproc_table][0][0];
2138
2139   ShowScores(scores,&sequence,NULL, bound,mode);
2140
2141 }
2142 /*---------------------------------------------------------------------------*/
2143
2144
2145 Font font;
2146
2147 OpenScreen(int mode, double log_bound, FILE *fin, FILE *flog,
2148                 FILE* flog_coil_conflicts,FILE *fout_coils,
2149                 FILE *fout,
2150                 double *m, double *b, double *m_single, double *b_single,
2151                 char lib[MAXFUNCTNUM],
2152                 char differ_lib[MAXFUNCTNUM],
2153                 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
2154                 int main_method,
2155                 int main_preprocessor_method,
2156                 int *seqcnt, int which_priors, int which_priorp,
2157                 double prior_freq_single[MAX_TABLE_NUMBER],
2158                 double prior_freq_pair[MAX_TABLE_NUMBER],
2159                 int good_turing_fixed_disc, int structural_pos[POSNUM+1])
2160 {
2161   int WIDTH, HEIGHT;
2162   XGCValues values;
2163   XSetWindowAttributes attr;
2164   Pixmap icon;
2165   static unsigned char dashl[2] = {1,1};
2166   static unsigned char dashd[2] = {2,2};
2167   int i;
2168
2169   WIDTH = 597;
2170   HEIGHT = 729;
2171   display = XOpenDisplay("");
2172
2173     /*  Check mode to see if want GUI. */
2174   if (!display || (mode & NO_GUI)) {
2175                      /** COmplete logfile, etc. if can't open xdisplay. */
2176     if (!(mode & NO_GUI))
2177       fprintf(stderr,"Sorry, bad display.\nMaybe your DISPLAY environment variable isn't set properly,\nor you haven't xhosted this machine.\n");
2178     if ( fout || ( (mode & PRN_MODE) || (mode & VER_MODE)) ) {
2179       fprintf(stderr,"\nThe log file will still be written by:\n\n");
2180       for (i = 0; motd[i][0]; ++i)   /* Print out title. */
2181         fprintf(stderr,"    %s\n",motd[i]);
2182       seqnum=sequp;
2183       finish_log(mode,log_bound,fin,flog,flog_coil_conflicts,fout_coils,
2184                  fout,
2185                  m,b, m_single, b_single,
2186                      lib,differ_lib,multi_lib, main_method, main_table,
2187                      main_preprocessor_method, seqcnt, avg_max, which_priors,
2188                      which_priorp, prior_freq_single, prior_freq_pair,
2189                      good_turing_fixed_disc, structural_pos);
2190     }
2191
2192     sclose(fin);
2193     if (flog) sclose(flog);
2194     if (flog_coil_conflicts) sclose(flog_coil_conflicts);
2195     if (fout) sclose(fout);
2196     if (ftotal_like)  sclose(ftotal_like);
2197     if (fout_coils)
2198       sclose(fout_coils);
2199     exit(0);
2200   }
2201
2202   font = XLoadFont(display,"9x15");
2203   window = XCreateSimpleWindow(display, RootWindow(display, 0),
2204                                0, 0,
2205                                WIDTH, HEIGHT/4 - 20,
2206                                0,
2207                                BlackPixel(display, 0),
2208                                WhitePixel(display, 0));
2209
2210   help   = XCreateSimpleWindow(display, RootWindow(display, 0),
2211                                0, 0,
2212                                WIDTH / 5 * 4 +60, HEIGHT / 5 * 3 + 200,
2213                                0,
2214                                BlackPixel(display, 0),
2215                                WhitePixel(display, 0));
2216   helpb  = XCreateSimpleWindow(display, window,
2217                               90, 0, 30, 15, 1,
2218                               BlackPixel(display, 0), WhitePixel(display, 0));
2219   next  = XCreateSimpleWindow(display, window,
2220                               125, 0, 30, 15, 1,
2221                               BlackPixel(display, 0), WhitePixel(display, 0));
2222   prev  = XCreateSimpleWindow(display, window,
2223                               160, 0, 30, 15, 1,
2224                               BlackPixel(display, 0), WhitePixel(display, 0));
2225   autob  = XCreateSimpleWindow(display, window,
2226                               195, 0, 30, 15, 1,
2227                               BlackPixel(display, 0), WhitePixel(display, 0));
2228   zoomb = XCreateSimpleWindow(display, window,
2229                               230, 0, 30, 15, 1,
2230                               BlackPixel(display, 0), WhitePixel(display, 0));
2231   printb = XCreateSimpleWindow(display, window,
2232                               265, 0, 36, 15, 1,
2233                               BlackPixel(display, 0), WhitePixel(display, 0));
2234   methodb=XCreateSimpleWindow(display, window,
2235                               306, 0, 42, 15, 1,
2236                               BlackPixel(display, 0), WhitePixel(display, 0));
2237 /***
2238   tableb=XCreateSimpleWindow(display, window,
2239                               353, 0, 36, 15, 1,
2240                               BlackPixel(display, 0), WhitePixel(display, 0));
2241 ***/
2242   paramb=XCreateSimpleWindow(display, window,
2243                               353, 0, 36, 30, 1,
2244                               BlackPixel(display, 0), WhitePixel(display, 0));
2245   quit  = XCreateSimpleWindow(display, window,
2246                               393, 0, 30, 15, 1,
2247                               BlackPixel(display, 0), WhitePixel(display, 0));
2248
2249 /** To get size:  **/
2250 /** Printing " trimer x.xx" takes 12 char = 12*7 =84 **/
2251 /** Each line takes 15 of y coord. giving 3*15 =45 + 10 for blank line **/
2252   seqscoreb=XCreateSimpleWindow(display, window,
2253                               XCOORD_SEQSCORES, YCOORD_SEQSCORES,
2254                               84, 55, 1,
2255                               BlackPixel(display, 0), WhitePixel(display, 0));
2256
2257   zoom   = XCreateSimpleWindow(display, RootWindow(display, 0),
2258                                0, 0,
2259                                400, 169,
2260                                0,
2261                                BlackPixel(display, 0),
2262                                WhitePixel(display, 0));
2263   XSelectInput(display, helpb, ButtonPressMask);
2264   XSelectInput(display, next, ButtonPressMask);
2265   XSelectInput(display, prev, ButtonPressMask);
2266   XSelectInput(display, autob, ButtonPressMask);
2267   XSelectInput(display, zoomb, ButtonPressMask);
2268   XSelectInput(display, printb, ButtonPressMask);
2269   XSelectInput(display, methodb, ButtonPressMask);
2270 /***
2271   XSelectInput(display, tableb, ButtonPressMask);
2272 ****/
2273   XSelectInput(display, paramb, ButtonPressMask);
2274
2275   XSelectInput(display, quit, ButtonPressMask);
2276   XSelectInput(display, window, ExposureMask);
2277   XSelectInput(display, zoom, ExposureMask | KeyPressMask | ButtonPressMask);
2278   XSelectInput(display, help, ExposureMask | KeyPressMask | ButtonPressMask);
2279   attr.bit_gravity = NorthWestGravity;
2280   XChangeWindowAttributes(display,window,CWBitGravity,&attr);
2281   XMapWindow(display, window);
2282   XFlush(display);
2283
2284   values.foreground = BlackPixel(display, 0);
2285   values.background = WhitePixel(display, 0);
2286   values.font = font;
2287   gc = XCreateGC(display, window, GCForeground | GCBackground | GCFont, &values);
2288   values.line_style = LineDoubleDash;
2289   values.font = XLoadFont(display,"6x13");
2290   gcgray = XCreateGC(display, window, GCForeground | GCBackground
2291                      | GCLineStyle | GCFont, &values);
2292   XSetDashes(display, gcgray, 1, dashl, 2);
2293   values.function = GXxor;
2294   values.foreground = 1;
2295   gcdash = XCreateGC(display, window, GCForeground
2296                      | GCLineStyle | GCFunction, &values);
2297   XSetDashes(display, gcdash, 0, dashd, 2);
2298
2299
2300 /*
2301   icon = XCreateBitmapFromData(display,window,bart_bits,bart_width,bart_height);*/
2302  /* XSetStandardProperties(display,window,"Coiled Coil Finder","xcoil",icon,NULL,0,NULL);*/
2303 /*  XSetStandardProperties(display,zoom,"XCoil Zoomer","xcoil",icon,NULL,0,NULL);*/
2304 /*  XSetStandardProperties(display,help,"XCoil Help","xcoil",icon,NULL,0,NULL);*/
2305     XFlush(display);
2306 }
2307
2308 PSDefs (FILE *ps, int seqlen)
2309 {
2310   double yscaler=1;
2311   double xscaler;
2312   int number_lines;
2313   int nominallength= 300;  /* Default value for number residues per line. */
2314
2315   if ( (seqlen%nominallength <= 50) && (seqlen >= nominallength)) {
2316     nominallength = 350;
2317   }
2318
2319 /*************************/
2320   if (ps_res_per_line != 0)  /*  0 signifies do the default as above. */
2321     if (ps_res_per_line == -1)  /* -1 signifies put the seq all on 1 line. */
2322       nominallength = seqlen;
2323     else
2324       nominallength = ps_res_per_line;  /* How many res per line.  */
2325 /************************/
2326
2327   if (seqlen ==0) number_lines=1;
2328   else
2329     number_lines = (int)(seqlen-1)/nominallength  +1;
2330                /* nominallength amino acids per line. */
2331               /*  C rounds down and we want 0-nominallength to go to 1 */
2332               /*  nominallength+1 to 2*noiminallength to go to 2.. */
2333
2334
2335   fputs("%!\n",ps);
2336
2337 /*****Scale and translate to fit the page if want it to rotate page.  */
2338 /* unscaled there are 9 lines per page. */
2339 /*                      .89 scaling if 300 res per line.
2340   xscaler = .89 * 300/nominallength;
2341   fputs("90 rotate\n",ps);
2342   fputs("78 -60 translate\n",ps);
2343   if (number_lines > 9) yscaler = (double) 9/number_lines;
2344   fprintf(ps,"%lf %lf scale\n", xscaler, yscaler);
2345 */
2346
2347 /*****Scale and translate to fit the page if don't want to rotate page.  */
2348 /* unscaled there can fit 12 lines per page, so 3600. */
2349   xscaler = .75 * 300/nominallength;   /* .75 scaling if 300 res per line. */
2350   fputs("25 740 translate\n",ps);
2351   if (number_lines > 12) yscaler = (double) 12/number_lines;
2352   fprintf(ps,"%lf %lf scale\n", xscaler,yscaler);
2353
2354
2355
2356   fputs("/black {0 setgray} def\n",ps);
2357   fputs("/gray {0.6 setgray} def\n",ps);
2358   fputs("/tt /Courier findfont 15 scalefont def\n",ps);
2359   fputs("/rm /Times-Roman findfont 12 scalefont def\n",ps);
2360   fputs("/cshow {dup stringwidth pop -2 div 0 rmoveto show}def\n",ps);
2361   fputs("/rshow {dup stringwidth pop neg 0 rmoveto show}def\n",ps);
2362   fputs("/m {neg moveto} def\n",ps);
2363   fputs("/l {neg lineto stroke} def\n",ps);
2364   fputs("/rl {neg rlineto} def\n",ps);
2365   fputs("/rl_stroke {neg rlineto stroke} def\n",ps);
2366   fputs("/bar {dup 0.0019 lt {pop} {0 exch rlineto stroke} ifelse} def\n",ps);
2367 }
2368
2369 /*** static int xmin = 60; *****  Made it global variable *****/
2370                           /*Originally was 29 */
2371                           /** xmin determines where in window start seq. **/
2372 static int nominallength;
2373 static int length;
2374 static int ybase;
2375 static int ydelta=57;
2376 /* If click at x,y, determine which residue position it corresponds to.
2377    Return 0 if none.
2378    */
2379 XYtoPos (int x, int y)
2380 {
2381   int i;
2382   if (x<xmin) return 0;
2383   i = (x-xmin)/2;
2384   if (i>= nominallength) return 0;
2385   for (y -= ybase+ydelta; y>0; y-=ydelta, i+=nominallength);
2386   if (y<-30) return 0;
2387   if (i>=length) return 0;
2388   return i+1;
2389 }
2390 PostoXY (int i, int *x, int *y)
2391 {
2392   --i;
2393   *x = xmin + 2*(i%nominallength);
2394   *y = ybase+ydelta + ydelta*(i/nominallength);
2395 }
2396 psstr (char *str, int len, FILE *ps)
2397 {
2398   int i;
2399   for (i=0; i<len; ++i)
2400     switch (str[i]) {
2401     case '(':
2402     case ')':
2403     case '\\':
2404       fputc('\\',ps);
2405     default:
2406       fputc(str[i],ps);
2407     }
2408 }
2409
2410
2411 /*** ybase determines where top of the seq score box should be so don't ***/
2412 /*** print over the sequenc title.  **/
2413 draw_seqscore_box(double bound, double seq_scores[MAX_CLASS_NUMBER],
2414                   Window *seqscoreb, FILE *ps, int ybase)
2415 {
2416   char buf[14];
2417 /*  XWindowChanges values; */
2418
2419
2420
2421   XClearWindow(display,*seqscoreb);
2422
2423  /** Relocate seqscoreb base on ybase**/
2424 /**This is complex way to move it.****
2425   values.y= YCOORD_SEQSCORES+ybase-30;
2426   XConfigureWindow(display,*seqscoreb,CWY,&values);
2427 ***/
2428
2429   XMoveWindow(display,*seqscoreb,XCOORD_SEQSCORES,  YCOORD_SEQSCORES+ybase-30);
2430
2431
2432
2433 /** 3 means print 3 into the box horiz, 10 means print 10 into box vert? */
2434
2435   XMapWindow(display, *seqscoreb);
2436   XDrawString(display,*seqscoreb,gcgray,3,10,"Seq Prob:",9);
2437
2438   if (ps) {     /* Draw the box. **/
2439     fprintf(ps,"%d %d m\n",10, 50+ybase-30);
2440     fprintf(ps,"0 40 rl\n");
2441     fprintf(ps,"79 0 rl\n");
2442     fprintf(ps,"0 -40 rl\n");
2443     fprintf(ps,"-79 0 rl_stroke\n");
2444   }
2445
2446   if (seq_scores[0] + seq_scores[1] + seq_scores[2] > 0) { /* Signals above */
2447     sprintf(buf," Dimer  %5.2lf", seq_scores[0]);         /*  bound.      */
2448     XDrawImageString(display,*seqscoreb,gcgray,3,35,buf,strlen(buf));
2449     sprintf(buf," Trimer %5.2lf", seq_scores[1]);
2450     XDrawImageString(display,*seqscoreb,gcgray,3,50,buf,strlen(buf));
2451
2452     if (ps) {
2453       fprintf(ps,"rm setfont ");
2454       fprintf(ps,"%d %d m (Seq Prob for) show\n", 15,60 +ybase-30);
2455
2456       fprintf(ps,"%d %d m (bound %5.2lf:) show\n", 15,70 +ybase-30,bound);
2457
2458       fprintf(ps,"%d %d m (   Dimer  %5.2lf) show\n",15,80 +ybase-30,
2459                                                 seq_scores[0]);
2460       fprintf(ps,"%d %d m (   Trimer %5.2lf) show\n",15,90 +ybase-30,
2461                                                 seq_scores[1]);
2462     }
2463
2464
2465   }
2466   else {
2467     sprintf(buf," No Residue  ");
2468     XDrawImageString(display,*seqscoreb,gcgray,3,35,buf,strlen(buf));
2469     sprintf(buf," above %5.2lf ", bound);
2470     XDrawImageString(display,*seqscoreb,gcgray,3,50,buf,strlen(buf));
2471
2472     if (ps) {
2473       fprintf(ps,"rm setfont ");
2474       fprintf(ps,"%d %d m (Seq Prob:) show\n", 15,60 +ybase-30);
2475
2476       fprintf(ps,"%d %d m ( No Residue) show\n", 15, 75 +ybase-30);
2477       fprintf(ps,"%d %d m ( above %5.2lf) show\n",15,85 +ybase-30, bound);
2478     }
2479
2480   }
2481
2482
2483 }
2484
2485 draw_bound_box(double bound, Window paramb)
2486 {
2487   char buf[6];
2488
2489   XMapWindow(display, paramb);
2490   XDrawString(display,paramb,gcgray,3,10,"bound",5);
2491   sprintf(buf,"%5.2lf", bound);
2492   XDrawImageString(display,paramb,gcgray,3,25,buf,strlen(buf));
2493 }
2494
2495 ShowScores (double *scores, Sequence *sequence, FILE *ps, double bound,
2496             int mode)
2497 {
2498   int i,q, x, y0, xmax;
2499   int width;
2500   int rows;
2501   int charperline;
2502   int xcoord;
2503   char MethodTitleArray[500];  /* For printing the methodname to ps file. */
2504   char *MethodTitle = MethodTitleArray;
2505   char PreprocTitleArray[500]; /* For printing the preproc.name to ps file. */
2506   char *PreprocTitle = PreprocTitleArray;
2507
2508
2509   MethodTitleArray[0] = '\0';  /* Initialization. */
2510   PreprocTitleArray[0]='\0';
2511
2512   nominallength = 300;
2513   length = sequence->seqlen;
2514
2515 /*  This makes display look nice (so don't have a line with < 50 residues. */
2516   if ((length%nominallength <= 50) && (length>= nominallength)){
2517     nominallength = 350;
2518 /*    if (ps) fputs("-50 0 translate\n",ps);   */
2519   }
2520
2521
2522 /*************************/
2523   if (ps)         /* Do this stuff only for postscript printout. */
2524     if (ps_res_per_line != 0)  /*  0 signifies do the default as above. */
2525       if (ps_res_per_line == -1)  /* -1 signifies put the seq all on 1 line. */
2526         nominallength = sequence->seqlen;
2527       else
2528         nominallength = ps_res_per_line;  /* How many res per line.  */
2529 /************************/
2530
2531
2532   charperline = (xmin+2*nominallength-5)/9;  /*Each char is "9" units witdth */
2533   rows = (length-1)/nominallength+1;    /* Number of rows to display. */
2534   if (rows==1)
2535     xmax = 2*length+xmin;              /* The max x position in current row. */
2536   else
2537     xmax = 2*nominallength+xmin;
2538   width = xmin+2*nominallength+ WIDTH_ADD_ON; /* Width of window to display. */
2539                                   /* The +WIDTH_ADD_ON used to be +11, but  */
2540                                   /* was cutting of some numbers. */
2541   x = strlen(sequence->title);
2542   XClearWindow(display,window);
2543
2544 /**********This is where window size is adjusted to fit the sequence.  ***/
2545 /********** using width and number of rows.     **************************/
2546
2547   XResizeWindow(display,window,width,
2548                 38+15*((x+charperline-1)/charperline) +ydelta*rows);
2549
2550 /*************************************************************************/
2551   XMapWindow(display, helpb);
2552   XMapWindow(display, next);
2553   XMapWindow(display, prev);
2554   XMapWindow(display, autob);
2555   XMapWindow(display, zoomb);
2556   XMapWindow(display, printb);
2557   XMapWindow(display, methodb);
2558 /****
2559   XMapWindow(display, tableb);
2560   XMapWindow(display, paramb);            Now done in draw_bound_box()
2561 ******/
2562   XMapWindow(display, quit);
2563   XSync(display,0);
2564   { XEvent report;
2565     while (XCheckTypedEvent(display,Expose,&report));
2566   }
2567   XDrawString(display,helpb,gcgray,3,10,"help",4);
2568   XDrawString(display,next,gcgray,3,10,"next",4);
2569   XDrawString(display,prev,gcgray,3,10,"prev",4);
2570   XDrawString(display,autob,gcgray,3,10,"auto",4);
2571   XDrawString(display,zoomb,gcgray,3,10,"zoom",4);
2572   XDrawString(display,printb,gcgray,3,10,"print",5);
2573   XDrawString(display,methodb,gcgray,3,10,"method",6);
2574
2575 /*  XDrawString(display,tableb,gcgray,3,10,"table",5); */
2576
2577   draw_bound_box(bound, paramb);
2578 /**  draw_seqscore_box(bound, seq_scores,?? &seqscoreb, ps, ybase);
2579   MOVED TO LATER (AFTER ybase gets adjusted at end of this function). ****/
2580
2581
2582   XDrawString(display,quit,gcgray,3,10,"quit",4);
2583
2584
2585 /************This section prints out the method name..... ********************/
2586 /****  Add 9 to xcoord for each character.                               *****/
2587   XDrawString(display,window,gc, XCOORD_SCORED_BY,15,"Scored by ",10);
2588
2589   make_methodtitle(MethodTitle, method, table, mode);
2590   XDrawString(display,window,gc,XCOORD_SCORED_BY +90,15,MethodTitle,
2591               strlen(MethodTitle));
2592
2593   XDrawString(display,window,gc, XCOORD_SCORED_BY,30,"Filter by ",10);
2594   make_methodtitle(PreprocTitle, preprocessor_method, preproc_table,
2595                    mode);
2596   XDrawString(display,window,gc,XCOORD_SCORED_BY + 90,30,
2597               PreprocTitle,strlen(PreprocTitle));
2598
2599 /***************************************************************************/
2600
2601   XDrawString(display,window,gc, 5,30,sequence->code,strlen(sequence->code));
2602   for (i=0,ybase=30; x > i + charperline; i += charperline)
2603     XDrawString(display,window,gc, 5,ybase+=15,&(sequence->title[i]),charperline);
2604   XDrawString(display,window,gc, 5,ybase+15,&(sequence->title[i]),x-i);
2605
2606   if (ps) {
2607     fprintf(ps,"tt setfont 5 15 m (%s) show\n",sequence->code);
2608 /*    fprintf(ps,"175 15 m (Viewgram by XCoil) show\n");  */
2609 /*
2610     fprintf(ps,"430 15 m (Scored by %s) show\n",methodname[method]);
2611 */
2612     fprintf(ps,"180 15 m (Scored by %s) show\n",MethodTitle);
2613     fprintf(ps,"180 26 m (Filter by %s) show\n",PreprocTitle);
2614
2615
2616     for (i=0,ybase=30; x > i + charperline; i += charperline) {
2617       fprintf(ps,"5 %d m (",ybase+=15);
2618       psstr(&(sequence->title[i]),charperline,ps);
2619       fprintf(ps,") show\n");
2620     }
2621     fprintf(ps,"5 %d m (",ybase+15);
2622     psstr(&(sequence->title[i]),x-i,ps);
2623     fprintf(ps,") show\n");
2624     fputs("rm setfont 2 setlinewidth\n",ps);
2625   }
2626   for (i=0,y0=ybase; i<length; ++i,x+=2) {
2627     double y;
2628     if (!(i%nominallength)) {x=xmin; y0 += ydelta;}
2629     y = y0-30*scores[8*i+7];
2630     XDrawLine(display, window, gc, x, y0, x, (int)(y+0.5));
2631     XDrawLine(display, window, gc, x+1, y0, x+1, (int)(y+0.5));
2632     if (ps) {
2633       fprintf(ps,"%d %d m %lf bar\n",x+1, y0, y0-y);
2634     }
2635   }
2636 /********************* If in pos mode draws in the correct likelihoods every
2637 now and then? ******************
2638   if (sequence->reg && method==HMMCoil) {
2639     if (ps)
2640       fputs("gray\n",ps);
2641     for (i=0,y0=ybase; i<length; ++i,x+=2) {
2642       if (!(i%nominallength)) {x=xmin; y0 += ydelta;}
2643       if (sequence->reg[i] != '.') {
2644         double y = y0-30*scores[8*i+sequence->reg[i]];
2645         XDrawLine(display, window, gcgray, x, y0, x, (int)(y+0.5));
2646         XDrawLine(display, window, gcgray, x+1, y0-1, x+1, (int)(y+0.5));
2647         if (ps)
2648         fprintf(ps,"%d %d m %lf bar\n",x+1, y0, y0-y);
2649       }
2650     }
2651     if (ps) fputs("black\n",ps);
2652   }
2653 **********************/
2654
2655 /**********Now draw in the PreProcScore in gray every other line.   ****/
2656
2657   if (ps)
2658     fputs("gray\n",ps);
2659   for (i=0,y0=ybase; i<length; i+=2,x+=4) { /* Normally is ++i,x+=4, but */
2660      double y;                              /* we do every other line here. */
2661     if (!(i%nominallength)) {x=xmin; y0 += ydelta;}
2662
2663     y = y0-30*preprocess_like[i*(POSNUM +1) + 7];
2664                            /* This is preprocess_like[i][7]. */
2665     XDrawLine(display, window, gcgray, x, y0, x, (int)(y+0.5));
2666     XDrawLine(display, window, gcgray, x+1, y0-1, x+1, (int)(y+0.5));
2667     if (ps)
2668       fprintf(ps,"%d %d m %lf bar\n",x+1, y0, y0-y);
2669   }
2670
2671   if (ps) fputs("black\n",ps);
2672
2673 /*************************************************************************/
2674
2675   if (ps) fputs("0 setlinewidth\n",ps);
2676   for (i=9,y0=ybase; i<length; i+=10) {
2677     if (!((i-9)%nominallength)) {x=xmin-1; y0 += ydelta;}
2678     x += 20;
2679     XDrawLine(display, window, gc, x,y0,x,y0+2);         /* Draw ticks every */
2680     if (ps) fprintf(ps,"%d %d m %d %d l\n",x, y0, x, y0+2);  /* 10 residues  */
2681
2682   }
2683   for (i=49,y0=ybase; i<length; i+=50) {
2684     char buff[6];
2685     if (!((i-49)%nominallength)) {x=xmin-1; y0 += ydelta;}
2686     x += 100;
2687     XDrawLine(display, window, gc, x,y0,x,y0+5);
2688     if (ps) fprintf(ps,"%d %d m %d %d l\n",x, y0, x, y0+5);
2689     sprintf(buff,"%d",i+1);
2690     XDrawString(display, window, gcgray, x-7, y0+18, buff, strlen(buff));
2691     if (ps) fprintf(ps,"%d %d m (%d) cshow\n",x,y0+18,i+1);
2692   }
2693   for (i=1,y0=ybase+ydelta; i<=rows; ++i,y0+= ydelta) {
2694     if (i==rows)
2695       xmax = 2*(length-(rows-1)*nominallength)+xmin;  /* Max x position in */
2696                                                       /* last row.         */
2697     XDrawString(display,window,gcgray,xmin-25,y0-26,"100%",4);
2698     XDrawString(display,window,gcgray,xmin-19,y0-11,"50%",3);
2699     XDrawString(display,window,gcgray,xmin-13,y0+4,"0%",2);
2700     XDrawLine(display,window,gcgray,xmin,y0-30,xmax-1,y0-30);
2701     XDrawLine(display,window,gcgray,xmin,y0-15,xmax-1,y0-15);
2702     if (ps) {
2703       fprintf(ps,"%d %d m (100%%) rshow\n",xmin-1,y0-26);
2704       fprintf(ps,"%d %d m (50%%) rshow\n",xmin-1,y0-11);
2705       fprintf(ps,"%d %d m (0%%) rshow\n",xmin-1,y0+4);
2706       fprintf(ps,"%d %d m %d %d l\n",xmin,y0-30,xmax,y0-30);
2707       fprintf(ps,"%d %d m %d %d l\n",xmin,y0-15,xmax,y0-15);
2708       fprintf(ps,"%d %d m %d %d l\n",xmin,y0,xmax,y0);
2709     }
2710   }
2711
2712
2713   draw_seqscore_box(bound, seq_scores, &seqscoreb, ps, ybase);
2714
2715   if (zoomptr) ZoomPnt();
2716   XFlush(display);
2717   if (ps) {
2718     fputs("showpage\n",ps);
2719     fflush(ps);
2720   }
2721 }
2722
2723
2724
2725
2726
2727 make_methodtitle(char *MethodTitle, int Method, int which_table,
2728                  int mode)
2729 {
2730   char buf[4];
2731   int xcoord;
2732   int like2or3;
2733
2734
2735   if (which_table ==0) like2or3 =2;
2736   else like2or3=3;
2737
2738
2739 /************This section prints out the method name..... ********************/
2740 /****  Add 9 to xcoord for each character.                               *****/
2741
2742
2743   xcoord =520;
2744   if ((Method == PairCoilDiff) || (Method == PairCoilDiffAvg)) {
2745                                /* Does Differen. use dimer or trimer like. */
2746     if (like2or3==2) {
2747 /*      XDrawString(display,window,gc,xcoord,15,"Dimer",5);  */
2748       MethodTitle= strcat(MethodTitle, "Dimer");
2749       xcoord += 45;
2750     }
2751     else if (like2or3==3) {
2752 /*      XDrawString(display,window,gc,xcoord,15,"Trimer",6); */
2753       MethodTitle= strcat(MethodTitle, "Trimer");
2754       xcoord += 54;
2755     }
2756   }
2757
2758
2759   if (Method == PairCoil) {
2760     if (mode & PAIRCOIL_PAIRS) { /*  Are we using pair or single prob */
2761 /*      XDrawString(display,window,gc,xcoord,15,"Pair",4); */
2762       MethodTitle= strcat(MethodTitle, "Pair");
2763       xcoord += 36;
2764     }
2765     else {
2766 /*      XDrawString(display,window,gc,xcoord,15,"Single",6);  */
2767       MethodTitle= strcat(MethodTitle, "Single");
2768       xcoord += 54;
2769     }
2770
2771   }   /*  Ends the stuff for dealing with PairCoil, TDRes, TDCoil.   */
2772
2773
2774 /*****************************************************************************/
2775 /*************************** Now write the methodname.  *********************/
2776
2777 /*  XDrawString(display,window,gc,xcoord,15,methodname[method],strlen(methodname[method]));   */
2778   xcoord += 9*strlen(methodname[Method]);
2779   MethodTitle= strcat(MethodTitle, methodname[Method]);
2780
2781   if (Method == MultiCoil) {
2782     if (mode & ONLY_COILED_CLASSES)
2783       strcat(MethodTitle,"Olig");
2784     if (Coil_Score) {
2785       if (Coil_Score ==1)
2786         if (weighted_avg)
2787           strcat(MethodTitle,"WtdAvg");
2788         else strcat(MethodTitle,"Avg");
2789
2790       else if (Coil_Score ==2)
2791         strcat(MethodTitle, "MAX");
2792
2793       if (start_coil_at_reg_shift)
2794         strcat(MethodTitle, "Shift");
2795       else
2796         strcat(MethodTitle, "NoShifts");
2797     }
2798   }
2799
2800    /* Now print Table are using.  */
2801    /* And if are using a particular offset, print that.     */
2802   if ( (Method==MultiCoil) || (Method==PairCoil) || (Method == HMMCoil)
2803        || (Method==PairCoilDiff) || (Method==PairCoilDiffAvg)) {
2804
2805
2806     if (number_tables >1) {
2807          sprintf(buf,"%d",which_table+1);
2808 /*         XDrawString(display,window,gc, xcoord,15,buf,strlen(buf)); */
2809          xcoord += 9*strlen(buf);
2810          MethodTitle= strcat(MethodTitle, buf);
2811        }
2812   }
2813
2814   sprintf(buf,"");  /* Clear out buf. */
2815
2816
2817   if ((offset_to_use != -1) && (offset_to_use != 7))
2818     sprintf(buf,"%c",'a' + offset_to_use);
2819   else if (Method == MultiCoil)
2820     if (offset_to_use == -1) {
2821       strcat(MethodTitle,"Max");
2822       xcoord += 9*3;
2823     }
2824     else if (offset_to_use == 7) {
2825       strcat(MethodTitle,"Combo");
2826       xcoord += 9*5;
2827     }
2828
2829   /*    XDrawString(display,window,gc, xcoord,15,buf,strlen(buf));  */
2830   xcoord += 9*strlen(buf);
2831   MethodTitle= strcat(MethodTitle, buf);
2832
2833 }
2834
2835 /***************************************************************************/
2836
2837
2838 int is_not_differentiator(int Method)
2839 {
2840   if ( (Method == MultiCoil) || (Method == PairCoil) || (Method==HMMCoil) ||
2841       (Method==NEWCOILS) || (Method == ActualCoil) )
2842     return (1);
2843   else return(0);
2844 }
2845
2846
2847
2848
2849
2850 log_file_header(FILE *flog, FILE *flog_coil_conflicts,
2851                      int mode, int argc, char *argv[],
2852                      int avg_max,  double bound,
2853                      int preprocessor_method, int preproc_table,
2854                      int log_offset_to_use, int start_coil_at_reg_shift,
2855                      int weighted_avg, int main_method, int main_table)
2856 {
2857   int i;
2858   char MethodTitleArray[500];  /* For printing the methodname to ps file. */
2859   char *MethodTitle = MethodTitleArray;
2860
2861   MethodTitleArray[0] = '\0';  /* Initialization. */
2862   make_methodtitle(MethodTitle, method,main_table,mode);
2863
2864   if (flog) {
2865     fprintf(flog,"\nInput file is %s", argc>1?argv[1]:"-");
2866     fprintf(flog, "\n\nMethod Used is %s", MethodTitle);
2867     fprintf(flog, "\nThe bound for the coil locater is %.2lf",bound);
2868   }
2869   if (flog_coil_conflicts) {
2870     fprintf(flog_coil_conflicts,"\nInput file is %s", argc>1?argv[1]:"-");
2871     if (mode & TST_MODE0)
2872       fprintf(flog_coil_conflicts,"\nSubtracting sequence off from cctb before scoring it.");
2873     if (mode & TST_MODE1)
2874       fprintf(flog_coil_conflicts,"\nSubtracting sequence off from trimer table before scoring it.");
2875   }
2876
2877
2878   if (mode & WEIGHTED_PROBS) {
2879     if (flog_coil_conflicts)
2880       fprintf(flog_coil_conflicts,"\nStatistically weighted PairCoil score is ");
2881     if (flog)
2882       fprintf(flog,"\nStatistically weighted PairCoil score is ");
2883   }
2884   else {
2885     if (flog_coil_conflicts)
2886       fprintf(flog_coil_conflicts,"\nPairCoil score is ");
2887     if (flog)
2888       fprintf(flog,"\nPairCoil score is ");
2889   }
2890   if (avg_max ==0) {
2891     if (flog_coil_conflicts)
2892       fprintf(flog_coil_conflicts,"AVERAGE of residue scores.");
2893     if (flog)
2894       fprintf(flog,"AVERAGE of residue scores.");
2895   }
2896   else {
2897     if (flog_coil_conflicts)
2898       fprintf(flog_coil_conflicts,"MAX of residue scores.");
2899     if (flog)
2900       fprintf(flog,"MAX of residue scores.");
2901   }
2902
2903   if (mode & USE_LIKE_LINE) {
2904     if (flog_coil_conflicts)
2905       fprintf(flog_coil_conflicts,"\nPaircoil used the likelihood line.");
2906     if (flog)
2907       fprintf(flog,"\nPaircoil used the likelihood line.");
2908   }
2909   else {
2910     if (flog_coil_conflicts)
2911       fprintf(flog_coil_conflicts,"\nPaircoil used direct prob. formula instead of likelihood line.");
2912     if (flog)
2913       fprintf(flog,"\nPaircoil used direct prob. formula instead of likelihood line.");
2914   }
2915
2916   if (flog_coil_conflicts)
2917     fprintf(flog_coil_conflicts, "\nThe distances used by PairCoil are: ");
2918   if (flog)
2919     fprintf(flog, "\nThe distances used by PairCoil are: ");
2920
2921     for (i=0; i<pair_functnum; i++) {
2922       if (flog_coil_conflicts)
2923         fprintf(flog_coil_conflicts,"%d ",i);
2924       if (flog)
2925         fprintf(flog,"%d ",i);
2926     }
2927
2928     if (flog_coil_conflicts) {
2929       fprintf(flog_coil_conflicts, "\n\nMethod Used is %s", MethodTitle);
2930       fprintf(flog_coil_conflicts, "\nThe bound for the coil locater is %.2lf",bound);
2931     }
2932
2933   if ( (main_method==PairCoilDiff) || (main_method==PairCoilDiffAvg)) {
2934     char MethodTitleArray[500];  /* For printing the methodname. */
2935     char *MethodTitle = MethodTitleArray;
2936     int oligimerization;
2937
2938     MethodTitleArray[0] = '\0';  /* Initialization. */
2939     make_methodtitle(MethodTitle, preprocessor_method, preproc_table,mode);
2940
2941     if (main_table ==0)
2942       oligimerization == 2;
2943     else oligimerization ==3;
2944
2945     if (flog)
2946       fprintf(flog, "\nLikelihood of Coils being %d-stranded using offset_method %d and preprocessor %s", oligimerization, log_offset_to_use,MethodTitle);
2947     if (flog_coil_conflicts)
2948       fprintf(flog_coil_conflicts, "\nLikelihood of Coils being %d-stranded using offset_method %d and preprocessor %s", oligimerization, log_offset_to_use,MethodTitle);
2949                                          /* Whether logfile is for 2 or */
2950                                          /* 3 stranded predictions.     */
2951
2952     if (main_method==PairCoilDiffAvg) {
2953       if (weighted_avg == 0) {
2954         if (flog)
2955           fprintf(flog, "\n Averaging over coils was done by coil length.");
2956         if (flog_coil_conflicts)
2957           fprintf(flog_coil_conflicts, "\n Averaging over coils was done by coil length.");
2958       }
2959       else  { /* weighted_avg ==1 */
2960         if (flog)
2961           fprintf(flog, "\n Averaging over coils was done using coil length weighted at each position by the preprocessor coil likelihood.");
2962         if (flog_coil_conflicts)
2963           fprintf(flog_coil_conflicts, "\n Averaging over coils was done using coil length weighted at each position by the preprocessor coil likelihood.");
2964       }
2965
2966       if (start_coil_at_reg_shift) {
2967         if (flog)
2968           fprintf(flog,"\n New coils were considered to start at register shifts.");
2969         if (flog_coil_conflicts)
2970           fprintf(flog_coil_conflicts,"\n New coils were considered to start at register shifts.");
2971       }
2972       else {
2973         if (flog)
2974           fprintf(flog,"\n Register shifts did not start a new coil, the new register scores were just averaged in.");
2975         if (flog_coil_conflicts)
2976           fprintf(flog_coil_conflicts,"\n Register shifts did not start a new coil, the new register scores were just averaged in.");
2977       }
2978     }
2979   }
2980
2981   if (flog) {
2982     fprintf(flog,"\nMain Table is Table %d\n\n",main_table);
2983
2984     if (mode & POS_MODE)
2985       fprintf(flog,"\nIf the start or end of a coil is listed inside (), then that position does not occur in the pos file.\n");
2986   }
2987 }
2988
2989
2990
2991
2992
2993 coil_tuple_output(int number_tables, char *lib, int mode, int avg_max,
2994                        double bound, Sequence sequence, FILE *fout_coils,
2995                        int number_multi_lib[MAX_TABLE_NUMBER],
2996                        int number_score_dim)
2997 {
2998   char offsets[MAXSEQLEN];
2999   int dist;
3000   int  tab;
3001   int i,j;
3002   extern double all_preprocess_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1];
3003   extern double *all_scores;
3004
3005
3006   for (tab=0; tab< number_tables; tab++) {
3007     for (dist=0; dist<number_multi_lib[tab]; dist++)  {
3008
3009 /*************************************************************************/
3010 /**********This computes offset of that table/dist method of scoring. ***/
3011
3012       for (i=0; i<sequence.seqlen; i++)  {
3013         for (j=0; j<POSNUM; j++)
3014           if (all_preprocess_like[tab + number_tables*(int)dist][i][j] ==
3015               all_preprocess_like[tab + number_tables*(int)dist][i][7]) {
3016             offsets[i] = (j-i) % POSNUM;
3017             if (offsets[i]<0) offsets[i] += POSNUM;
3018           }
3019         /* This is a new thing so that if an offset change is not real */
3020             /* we keep the old offset. **********/
3021         if ( (i>0) && (offsets[i-1] != -1) &&
3022             (offsets[i] != offsets[i-1]) &&
3023             (all_preprocess_like[tab + number_tables*(int)dist]
3024              [i][(i+offsets[i-1]) %POSNUM]
3025              == all_preprocess_like[tab+number_tables*(int)dist][i][7]) )
3026           offsets[i] = offsets[i-1];
3027       }
3028
3029 /**************************************************************************/
3030 /****** Now compute coil scores for that method (table,dist) pair.     ***/
3031       get_raw_coil_score(sequence,
3032                          all_preprocess_like[tab+number_tables*(int)dist],
3033                          all_scores[tab+number_tables*(int)dist], mode,
3034                          avg_max,bound,offsets);
3035     }
3036   }
3037
3038 /************* Now output the tuples of coil scores.   ****/
3039   tuple_output2(sequence, mode, fout_coils,
3040                 all_scores, number_tables,
3041                 number_multi_lib,bound, number_score_dim);
3042 }
3043
3044
3045
3046 switch_gauss_param(int Offset_to_Use, int Method)
3047 {
3048   if (Method == MultiCoil) {
3049     if (Offset_to_Use == 7) {  /* Combo offset */
3050       determinants= &both_determinants[0][0];
3051       inverse_covars= &both_inverse_covars[0][0][0][0];
3052       means_submatrix= &both_means_submatrix[0][0][0];
3053     }
3054     else   {      /* max offset or some fixed offset */
3055       determinants= &both_determinants[1][0];
3056       inverse_covars= &both_inverse_covars[1][0][0][0];
3057       means_submatrix= &both_means_submatrix[1][0][0];
3058     }
3059   }
3060   else if (Method == PairCoilDiff) {
3061     determinants= &differ_determinants[0];
3062     inverse_covars= &differ_inverse_covars[0][0][0];
3063     means_submatrix= &differ_means_submatrix[0][0];
3064   }
3065
3066
3067 }
3068
3069
3070 void usage(char *prog_command)
3071 {
3072   int i;
3073   char *prog_name;
3074
3075 /** Strip off directories to get just the program name.  **/
3076   i=strlen(prog_command);
3077   while ((i>=0) && (prog_command[i] != '/') )
3078     i--;
3079   prog_name= &prog_command[i+1];
3080
3081   fprintf(stderr,"\nUsage: %s <file-name> [dimer, trimer, dimer-1, trimer-1, negatives, config_file] [multi_lib_dim1] [multi_lib_dim2] \n",prog_name);
3082
3083 }
3084
3085
3086
3087
3088
3089
3090
3091
3092
3093
3094
3095
3096
3097