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