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