binaries/src/ViennaRNA/lib/PS_dot.c

   1 /*
   2         PostScript and GML output for RNA secondary structures
   3                     and pair probability matrices
   4
   5                  c  Ivo Hofacker and Peter F Stadler
   6                           Vienna RNA package
   7 */
   8 #include <config.h>
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <math.h>
  12 #include <string.h>
  13 #include <ctype.h>
  14 #include "utils.h"
  15 #include "fold_vars.h"
  16 #include "PS_dot.h"
  17 #include "pair_mat.h"
  18 #include "aln_util.h"
  19 #include "gquad.h"
  20 #include "plot_layouts.h"
  21
  22 static char UNUSED rcsid[] = "$Id: PS_dot.c,v 1.38 2007/02/02 15:18:13 ivo Exp $";
  23
  24 #ifndef PI
  25 #define  PI       3.141592654
  26 #endif
  27 #define  PIHALF       PI/2.
  28 #define SIZE 452.
  29 #define PMIN 0.00001
  30
  31 /* local functions */
  32 PRIVATE FILE  *PS_dot_common(char *seq, char *wastlfile, char *comment,
  33                              int winsize);
  34 PRIVATE char **annote(const char *structure, const char *AS[]);
  35
  36
  37 /*---------------------------------------------------------------------------*/
  38
  39 /* options for gml output:
  40    uppercase letters: print sequence labels
  41    lowercase letters: no sequence lables
  42    graphics information:
  43    x X  simple xy plot
  44    (nothing else implemented at present)
  45    default:           no graphics data at all
  46 */
  47
  48 PUBLIC int gmlRNA(char *string, char *structure, char *ssfile, char option)
  49 {
  50   FILE *gmlfile;
  51   int i;
  52   int length;
  53   int labels=0;
  54   short *pair_table;
  55   float *X, *Y;
  56
  57   if (isupper(option)) labels = 1;
  58
  59   gmlfile = fopen(ssfile, "w");
  60   if (gmlfile == NULL) {
  61      fprintf(stderr, "can't open file %s - not doing xy_plot\n", ssfile);
  62      return 0;
  63   }
  64
  65   length = strlen(string);
  66
  67   pair_table = make_pair_table(structure);
  68
  69   switch(option){
  70   case 'X' :
  71   case 'x' :
  72     /* Simple XY Plot */
  73     X = (float *) space((length+1)*sizeof(float));
  74     Y = (float *) space((length+1)*sizeof(float));
  75     if (rna_plot_type == 0)
  76       i = simple_xy_coordinates(pair_table, X, Y);
  77     else
  78       i = naview_xy_coordinates(pair_table, X, Y);
  79
  80     if(i!=length) fprintf(stderr,"strange things happening in gmlRNA ...\n");
  81     break;
  82   default:
  83     /* No Graphics Information */
  84     X = NULL;
  85     Y = NULL;
  86   }
  87
  88   fprintf(gmlfile,
  89           "# Vienna RNA Package %s\n"
  90           "# GML Output\n"
  91           "# CreationDate: %s\n"
  92           "# Name: %s\n"
  93           "# Options: %s\n", VERSION, time_stamp(), ssfile, option_string());
  94   fprintf(gmlfile,
  95           "graph [\n"
  96           " directed 0\n");
  97   for (i=1; i<=length; i++){
  98      fprintf(gmlfile,
  99           " node [ id %d ", i);
 100      if (option) fprintf(gmlfile,
 101           "label \"%c\"",string[i-1]);
 102      if ((option == 'X')||(option=='x'))
 103        fprintf(gmlfile,
 104                "\n  graphics [ x %9.4f y %9.4f ]\n", X[i-1], Y[i-1]);
 105      fprintf(gmlfile," ]\n");
 106   }
 107   for (i=1; i<length; i++)
 108     fprintf(gmlfile,
 109             "edge [ source %d target %d ]\n", i, i+1);
 110   for (i=1; i<=length; i++) {
 111      if (pair_table[i]>i)
 112         fprintf(gmlfile,
 113                 "edge [ source %d target %d ]\n", i, pair_table[i]);
 114   }
 115   fprintf(gmlfile, "]\n");
 116   fclose(gmlfile);
 117
 118   free(pair_table);
 119   free(X); free(Y);
 120   return 1; /* success */
 121 }
 122
 123 /*---------------------------------------------------------------------------*/
 124 int PS_rna_plot(char *string, char *structure, char *ssfile) {
 125   return PS_rna_plot_a(string, structure, ssfile, NULL, NULL);
 126 }
 127
 128 static const char *RNAss_head =
 129 "%%BeginProlog\n"
 130 "/RNAplot 100 dict def\n"
 131 "RNAplot begin\n"
 132 "/fsize  14 def\n"
 133 "/outlinecolor {0.2 setgray} bind def\n"
 134 "/paircolor    {0.2 setgray} bind def\n"
 135 "/seqcolor     {0   setgray} bind def\n"
 136 "/cshow  { dup stringwidth pop -2 div fsize -3 div rmoveto show} bind def\n"
 137 "/min { 2 copy gt { exch } if pop } bind def\n"
 138 "/max { 2 copy lt { exch } if pop } bind def\n"
 139 "/arccoords { % i j arccoords\n"
 140 "  % puts optimal x1 y1 x2 y2 coordinates used in bezier curves from i to j\n"
 141 "  % onto the stack\n"
 142 "  dup 3 -1 roll dup 4 -1 roll lt dup dup 5 2 roll {exch} if\n"
 143 "  dup 3 -1 roll dup 3 -1 roll exch sub 1 sub dup\n"
 144 "  4 -2 roll 5 -1 roll {exch} if 4 2 roll\n"
 145 "  sequence length dup 2 div exch 3 1 roll lt \n"
 146 "  {exch 5 -1 roll pop 4 -2 roll exch 4 2 roll}\n"
 147 "  { 4 2 roll 5 -1 roll dup 6 1 roll {exch} if\n"
 148 "    4 -2 roll exch pop dup 3 -1 roll dup 4 1 roll\n"
 149 "    exch add 4 -1 roll dup 5 1 roll sub 1 sub\n"
 150 "    5 -1 roll not {4 -2 roll exch 4 2 roll} if\n"
 151 "  }ifelse\n"
 152 "   % compute the scalingfactor and prepare (1-sf) and sf*r\n"
 153 "  2 mul exch cpr 3 1 roll div dup\n"
 154 "  3 -1 roll mul exch 1 exch sub exch\n"
 155 "   % compute the coordinates\n"
 156 "  3 -1 roll 1 sub coor exch get aload pop % get coord for i\n"
 157 "  4 -1 roll dup 5 1 roll mul 3 -1 roll dup 4 1 roll add exch % calculate y1\n"
 158 "  4 -1 roll dup 5 1 roll mul 3 -1 roll dup 4 1 roll add exch % calculate x1\n"
 159 "  5 -1 roll 1 sub coor exch get aload pop % get coord for j\n"
 160 "  % duplicate j coord\n"
 161 "  dup 3 -1 roll dup 4 1 roll exch 8 2 roll\n"
 162 "  6 -1 roll dup 7 1 roll mul 5 -1 roll dup 6 1 roll add exch % calculate y2\n"
 163 "  6 -1 roll mul 5 -1 roll add exch % calculate x2\n"
 164 "  6 -2 roll % reorder\n"
 165 "} bind def\n"
 166 "/drawoutline {\n"
 167 "  gsave outlinecolor newpath\n"
 168 "  coor 0 get aload pop 0.8 0 360 arc % draw 5' circle of 1st sequence\n"
 169 "  currentdict /cutpoint known        % check if cutpoint is defined\n"
 170 "  {coor 0 cutpoint getinterval\n"
 171 "   {aload pop lineto} forall         % draw outline of 1st sequence\n"
 172 "   coor cutpoint 1 add get aload pop\n"
 173 "   2 copy moveto 0.8 0 360 arc       % draw 5' circle of 2nd sequence\n"
 174 "   coor cutpoint 1 add coor length cutpoint 1 add sub getinterval\n"
 175 "   {aload pop lineto} forall}        % draw outline of 2nd sequence\n"
 176 "  {coor {aload pop lineto} forall}   % draw outline as a whole\n"
 177 "  ifelse\n"
 178 "  stroke grestore\n"
 179 "} bind def\n"
 180 "/drawpairs {\n"
 181 "  paircolor\n"
 182 "  0.7 setlinewidth\n"
 183 "  [9 3.01] 9 setdash\n"
 184 "  newpath\n"
 185 "  pairs {aload pop\n"
 186 "      currentdict (cpr) known\n"
 187 "      { exch dup\n"
 188 "        coor  exch 1 sub get aload pop moveto\n"
 189 "        exch arccoords curveto\n"
 190 "      }\n"
 191 "      { coor exch 1 sub get aload pop moveto\n"
 192 "        coor exch 1 sub get aload pop lineto\n"
 193 "      }ifelse\n"
 194 "  } forall\n"
 195 "  stroke\n"
 196 "} bind def\n"
 197 "% draw bases\n"
 198 "/drawbases {\n"
 199 "  [] 0 setdash\n"
 200 "  seqcolor\n"
 201 "  0\n"
 202 "  coor {\n"
 203 "    aload pop moveto\n"
 204 "    dup sequence exch 1 getinterval cshow\n"
 205 "    1 add\n"
 206 "  } forall\n"
 207 "  pop\n"
 208 "} bind def\n\n"
 209 "/init {\n"
 210 "  /Helvetica findfont fsize scalefont setfont\n"
 211 "  1 setlinejoin\n"
 212 "  1 setlinecap\n"
 213 "  0.8 setlinewidth\n"
 214 "  72 216 translate\n"
 215 "  % find the coordinate range\n"
 216 "  /xmax -1000 def /xmin 10000 def\n"
 217 "  /ymax -1000 def /ymin 10000 def\n"
 218 "  coor {\n"
 219 "      aload pop\n"
 220 "      dup ymin lt {dup /ymin exch def} if\n"
 221 "      dup ymax gt {/ymax exch def} {pop} ifelse\n"
 222 "      dup xmin lt {dup /xmin exch def} if\n"
 223 "      dup xmax gt {/xmax exch def} {pop} ifelse\n"
 224 "  } forall\n"
 225 "  /size {xmax xmin sub ymax ymin sub max} bind def\n"
 226 "  72 6 mul size div dup scale\n"
 227 "  size xmin sub xmax sub 2 div size ymin sub ymax sub 2 div\n"
 228 "  translate\n"
 229 "} bind def\n"
 230 "end\n";
 231
 232 static const char *anote_macros =
 233 "RNAplot begin\n"
 234 "% extra definitions for standard anotations\n"
 235 "/min { 2 copy gt { exch } if pop } bind def\n"
 236 "/BLACK { 0 0 0 } def\n"
 237 "/RED   { 1 0 0 } def\n"
 238 "/GREEN { 0 1 0 } def\n"
 239 "/BLUE  { 0 0 1 } def\n"
 240 "/WHITE { 1 1 1 } def\n"
 241 "/LabelFont { % font size LabelFont\n"
 242 "  exch findfont exch fsize mul scalefont setfont\n"
 243 "} bind def\n"
 244 "/Label { % i dx dy (text) Label\n"
 245 "  % write text at base i plus offset dx, dy\n"
 246 "  4 3 roll 1 sub coor exch get aload pop moveto\n"
 247 "  3 1 roll fsize mul exch fsize mul exch rmoveto\n"
 248 "  show\n"
 249 "} bind def\n"
 250 "/cmark { % i cmark   draw circle around base i\n"
 251 "  newpath 1 sub coor exch get aload pop\n"
 252 "  fsize 2 div 0 360 arc stroke\n"
 253 "} bind def\n"
 254 "/gmark { % i j c gmark\n"
 255 "  % draw basepair i,j with c counter examples in gray\n"
 256 "  gsave\n"
 257 "  3 min [0 0.33 0.66 0.9] exch get setgray\n"
 258 "  1 sub dup coor exch get aload pop moveto\n"
 259 "  sequence exch 1 getinterval cshow\n"
 260 "  1 sub dup coor exch get aload pop moveto\n"
 261 "  sequence exch 1 getinterval cshow\n"
 262 "  grestore\n"
 263 "} bind def\n"
 264 "/segmark { % f i j lw r g b segmark\n"
 265 "  % mark segment [i,j] with outline width lw and color rgb\n"
 266 "  % use omark and Fomark instead\n"
 267 "  gsave\n"
 268 "  setrgbcolor setlinewidth\n"
 269 "  newpath\n"
 270 "  1 sub exch 1 sub dup\n"
 271 "  coor exch get aload pop moveto\n"
 272 "  currentdict (cpr) known\n"
 273 "  {\n"
 274 "    3 -1 roll dup 4 1 roll dup\n"
 275 "    {\n"
 276 "      3 1 roll dup 3 -1 roll dup\n"
 277 "      4 1 roll exch 5 2 roll exch\n"
 278 "    }\n"
 279 "    {\n"
 280 "      3 1 roll exch\n"
 281 "    } ifelse\n"
 282 "    1 exch { coor exch get aload pop lineto } for\n"
 283 "    {\n"
 284 "      dup 3 1 roll 1 add exch 1 add arccoords pop pop\n"
 285 "      4 2 roll 5 -1 roll coor exch get aload pop curveto\n"
 286 "    } if\n"
 287 "  }\n"
 288 "  {\n"
 289 "    exch 1 exch {\n"
 290 "      coor exch get aload pop lineto\n"
 291 "    } for\n"
 292 "  } ifelse\n"
 293 "  { closepath fill } if  stroke\n"
 294 "  grestore\n"
 295 "} bind def\n"
 296 "/omark { % i j lw r g b omark\n"
 297 "  % stroke segment [i..j] with linewidth lw, color rgb\n"
 298 "  false 7 1 roll segmark\n"
 299 "} bind def\n"
 300 "/Fomark { % i j r g b Fomark\n"
 301 "  % fill segment [i..j] with color rgb\n"
 302 "  % should precede drawbases\n"
 303 "  1 4 1 roll true 7 1 roll segmark\n"
 304 "} bind def\n"
 305 "/BFmark{ % i j k l r g b BFmark\n"
 306 "  % fill block between pairs (i,j) and (k,l) with color rgb\n"
 307 "  % should precede drawbases\n"
 308 "  gsave\n"
 309 "  setrgbcolor\n"
 310 "  newpath\n"
 311 "  currentdict (cpr) known\n"
 312 "  {\n"
 313 "    dup 1 sub coor exch get aload pop moveto % move to l\n"
 314 "    dup 1 sub 4 -1 roll dup 5 1 roll 1 sub 1 exch\n"
 315 "    { coor exch get aload pop lineto } for % lines from l to j\n"
 316 "    3 -1 roll 4 -1 roll dup 5 1 roll arccoords curveto % curve from j to i\n"
 317 "    exch dup 4 -1 roll 1 sub exch 1 sub 1 exch\n"
 318 "    { coor exch get aload pop lineto } for % lines from i to k\n"
 319 "    exch arccoords curveto% curve from k to l\n"
 320 "  }\n"
 321 "  {  exch 4 3 roll exch 1 sub exch 1 sub dup\n"
 322 "     coor exch get aload pop moveto\n"
 323 "     exch 1 exch { coor exch get aload pop lineto } for\n"
 324 "     exch 1 sub exch 1 sub dup\n"
 325 "     coor exch get aload pop lineto\n"
 326 "     exch 1 exch { coor exch get aload pop lineto } for\n"
 327 "  } ifelse\n"
 328 "    closepath fill stroke\n"
 329 "   grestore\n"
 330 "} bind def\n"
 331 "/hsb {\n"
 332 "  dup 0.3 mul 1 exch sub sethsbcolor\n"
 333 "} bind def\n"
 334 "/colorpair { % i j hue sat colorpair\n"
 335 "  % draw basepair i,j in color\n"
 336 "  % 1 index 0.00 ne {\n"
 337 "  gsave\n"
 338 "  newpath\n"
 339 "  hsb\n"
 340 "  fsize setlinewidth\n"
 341 "  currentdict (cpr) known\n"
 342 "  {\n"
 343 "    exch dup\n"
 344 "    coor  exch 1 sub get aload pop moveto\n"
 345 "    exch arccoords curveto\n"
 346 "  }\n"
 347 "  { 1 sub coor exch get aload pop moveto\n"
 348 "    1 sub coor exch get aload pop lineto\n"
 349 "  } ifelse\n"
 350 "   stroke\n"
 351 "   grestore\n"
 352 "   % } if\n"
 353 "} bind def\n"
 354  "end\n\n";
 355
 356 int PS_rna_plot_a(char *string, char *structure, char *ssfile, char *pre, char *post)
 357 {
 358   float  xmin, xmax, ymin, ymax, size;
 359   int    i, length;
 360   float *X, *Y;
 361   FILE  *xyplot;
 362   short *pair_table;
 363   char *c;
 364
 365   length = strlen(string);
 366
 367   xyplot = fopen(ssfile, "w");
 368   if (xyplot == NULL) {
 369     fprintf(stderr, "can't open file %s - not doing xy_plot\n", ssfile);
 370     return 0;
 371   }
 372
 373   pair_table = make_pair_table(structure);
 374
 375   X = (float *) space((length+1)*sizeof(float));
 376   Y = (float *) space((length+1)*sizeof(float));
 377   switch(rna_plot_type){
 378     case VRNA_PLOT_TYPE_SIMPLE:   i = simple_xy_coordinates(pair_table, X, Y);
 379                                   break;
 380     case VRNA_PLOT_TYPE_CIRCULAR: {
 381                                     int radius = 3*length;
 382                                     i = simple_circplot_coordinates(pair_table, X, Y);
 383                                     for (i = 0; i < length; i++) {
 384                                       X[i] *= radius;
 385                                       X[i] += radius;
 386                                       Y[i] *= radius;
 387                                       Y[i] += radius;
 388                                     }
 389                                   }
 390                                   break;
 391     default:                      i = naview_xy_coordinates(pair_table, X, Y);
 392                                   break;
 393   }
 394   if(i!=length) fprintf(stderr,"strange things happening in PS_rna_plot...\n");
 395
 396   xmin = xmax = X[0];
 397   ymin = ymax = Y[0];
 398   for (i = 1; i < length; i++) {
 399      xmin = X[i] < xmin ? X[i] : xmin;
 400      xmax = X[i] > xmax ? X[i] : xmax;
 401      ymin = Y[i] < ymin ? Y[i] : ymin;
 402      ymax = Y[i] > ymax ? Y[i] : ymax;
 403   }
 404   size = MAX2((xmax-xmin),(ymax-ymin));
 405
 406   fprintf(xyplot,
 407           "%%!PS-Adobe-3.0 EPSF-3.0\n"
 408           "%%%%Creator: %s, ViennaRNA-%s\n"
 409           "%%%%CreationDate: %s"
 410           "%%%%Title: RNA Secondary Structure Plot\n"
 411           "%%%%BoundingBox: 66 210 518 662\n"
 412           "%%%%DocumentFonts: Helvetica\n"
 413           "%%%%Pages: 1\n"
 414           "%%%%EndComments\n\n"
 415           "%%Options: %s\n", rcsid+5, VERSION, time_stamp(), option_string());
 416   fprintf(xyplot, "%% to switch off outline pairs of sequence comment or\n"
 417           "%% delete the appropriate line near the end of the file\n\n");
 418   fprintf(xyplot, "%s", RNAss_head);
 419
 420   if (pre || post) {
 421     fprintf(xyplot, "%s", anote_macros);
 422   }
 423   fprintf(xyplot, "%%%%EndProlog\n");
 424
 425   fprintf(xyplot, "RNAplot begin\n"
 426           "%% data start here\n");
 427
 428   /* cut_point */
 429   if ((c = strchr(structure, '&'))) {
 430     int cutpoint;
 431     cutpoint = c - structure;
 432     string[cutpoint] = ' '; /* replace & with space */
 433     fprintf(xyplot, "/cutpoint %d def\n", cutpoint);
 434   }
 435
 436   /* sequence */
 437   fprintf(xyplot,"/sequence (\\\n");
 438   i=0;
 439   while (i<length) {
 440     fprintf(xyplot, "%.255s\\\n", string+i);  /* no lines longer than 255 */
 441     i+=255;
 442   }
 443   fprintf(xyplot,") def\n");
 444   /* coordinates */
 445   fprintf(xyplot, "/coor [\n");
 446   for (i = 0; i < length; i++)
 447     fprintf(xyplot, "[%3.8f %3.8f]\n", X[i], Y[i]);
 448   fprintf(xyplot, "] def\n");
 449   /* correction coordinates for quadratic beziers in case we produce a circplot */
 450   if(rna_plot_type == VRNA_PLOT_TYPE_CIRCULAR)
 451     fprintf(xyplot, "/cpr %6.2f def\n", (float)3*length);
 452   /* base pairs */
 453   fprintf(xyplot, "/pairs [\n");
 454   for (i = 1; i <= length; i++)
 455     if (pair_table[i]>i)
 456       fprintf(xyplot, "[%d %d]\n", i, pair_table[i]);
 457   fprintf(xyplot, "] def\n\n");
 458
 459   fprintf(xyplot, "init\n\n");
 460   /* draw the data */
 461   if (pre) {
 462     fprintf(xyplot, "%% Start Annotations\n");
 463     fprintf(xyplot, "%s\n", pre);
 464     fprintf(xyplot, "%% End Annotations\n");
 465   }
 466   fprintf(xyplot,
 467           "%% switch off outline pairs or bases by removing these lines\n"
 468           "drawoutline\n"
 469           "drawpairs\n"
 470           "drawbases\n");
 471
 472   if (post) {
 473     fprintf(xyplot, "%% Start Annotations\n");
 474     fprintf(xyplot, "%s\n", post);
 475     fprintf(xyplot, "%% End Annotations\n");
 476   }
 477   fprintf(xyplot, "%% show it\nshowpage\n");
 478   fprintf(xyplot, "end\n");
 479   fprintf(xyplot, "%%%%EOF\n");
 480
 481   fclose(xyplot);
 482
 483   free(pair_table);
 484   free(X); free(Y);
 485   return 1; /* success */
 486 }
 487
 488 int PS_rna_plot_a_gquad(char *string,
 489                         char *structure,
 490                         char *ssfile,
 491                         char *pre,
 492                         char *post){
 493   float  xmin, xmax, ymin, ymax, size;
 494   int    i, length;
 495   int    ee, gb, ge, Lg, l[3];
 496   float *X, *Y;
 497   FILE  *xyplot;
 498   short *pair_table, *pair_table_g;;
 499   char *c;
 500
 501   length = strlen(string);
 502
 503   xyplot = fopen(ssfile, "w");
 504   if (xyplot == NULL) {
 505     fprintf(stderr, "can't open file %s - not doing xy_plot\n", ssfile);
 506     return 0;
 507   }
 508
 509   pair_table = make_pair_table(structure);
 510   pair_table_g = make_pair_table(structure);
 511
 512   ge=0;
 513   while ( (ee=parse_gquad(structure+ge, &Lg, l)) >0 ) {
 514     ge += ee;
 515     gb=ge-Lg*4-l[0]-l[1]-l[2]+1;
 516     /* add pseudo-base pair encloding gquad */
 517     for (i=0; i<Lg; i++) {
 518       pair_table_g[ge-i]=gb+i;
 519       pair_table_g[gb+i]=ge-i;
 520     }
 521   }
 522
 523   X = (float *) space((length+1)*sizeof(float));
 524   Y = (float *) space((length+1)*sizeof(float));
 525   switch(rna_plot_type){
 526     case VRNA_PLOT_TYPE_SIMPLE:   i = simple_xy_coordinates(pair_table_g, X, Y);
 527                                   break;
 528     case VRNA_PLOT_TYPE_CIRCULAR: {
 529                                     int radius = 3*length;
 530                                     i = simple_circplot_coordinates(pair_table_g, X, Y);
 531                                     for (i = 0; i < length; i++) {
 532                                       X[i] *= radius;
 533                                       X[i] += radius;
 534                                       Y[i] *= radius;
 535                                       Y[i] += radius;
 536                                     }
 537                                   }
 538                                   break;
 539     default:                      i = naview_xy_coordinates(pair_table_g, X, Y);
 540                                   break;
 541   }
 542   if(i!=length) fprintf(stderr,"strange things happening in PS_rna_plot...\n");
 543
 544   xmin = xmax = X[0];
 545   ymin = ymax = Y[0];
 546   for (i = 1; i < length; i++) {
 547      xmin = X[i] < xmin ? X[i] : xmin;
 548      xmax = X[i] > xmax ? X[i] : xmax;
 549      ymin = Y[i] < ymin ? Y[i] : ymin;
 550      ymax = Y[i] > ymax ? Y[i] : ymax;
 551   }
 552   size = MAX2((xmax-xmin),(ymax-ymin));
 553
 554   fprintf(xyplot,
 555           "%%!PS-Adobe-3.0 EPSF-3.0\n"
 556           "%%%%Creator: %s, ViennaRNA-%s\n"
 557           "%%%%CreationDate: %s"
 558           "%%%%Title: RNA Secondary Structure Plot\n"
 559           "%%%%BoundingBox: 66 210 518 662\n"
 560           "%%%%DocumentFonts: Helvetica\n"
 561           "%%%%Pages: 1\n"
 562           "%%%%EndComments\n\n"
 563           "%%Options: %s\n", rcsid+5, VERSION, time_stamp(), option_string());
 564   fprintf(xyplot, "%% to switch off outline pairs of sequence comment or\n"
 565           "%% delete the appropriate line near the end of the file\n\n");
 566   fprintf(xyplot, "%s", RNAss_head);
 567
 568   if (pre || post) {
 569     fprintf(xyplot, "%s", anote_macros);
 570   }
 571   fprintf(xyplot, "%%%%EndProlog\n");
 572
 573   fprintf(xyplot, "RNAplot begin\n"
 574           "%% data start here\n");
 575
 576   /* cut_point */
 577   if ((c = strchr(structure, '&'))) {
 578     int cutpoint;
 579     cutpoint = c - structure;
 580     string[cutpoint] = ' '; /* replace & with space */
 581     fprintf(xyplot, "/cutpoint %d def\n", cutpoint);
 582   }
 583
 584   /* sequence */
 585   fprintf(xyplot,"/sequence (\\\n");
 586   i=0;
 587   while (i<length) {
 588     fprintf(xyplot, "%.255s\\\n", string+i);  /* no lines longer than 255 */
 589     i+=255;
 590   }
 591   fprintf(xyplot,") def\n");
 592   /* coordinates */
 593   fprintf(xyplot, "/coor [\n");
 594   for (i = 0; i < length; i++)
 595     fprintf(xyplot, "[%3.8f %3.8f]\n", X[i], Y[i]);
 596   fprintf(xyplot, "] def\n");
 597   /* correction coordinates for quadratic beziers in case we produce a circplot */
 598   if(rna_plot_type == VRNA_PLOT_TYPE_CIRCULAR)
 599     fprintf(xyplot, "/cpr %6.2f def\n", (float)3*length);
 600   /* base pairs */
 601   fprintf(xyplot, "/pairs [\n");
 602   for (i = 1; i <= length; i++)
 603     if (pair_table[i]>i)
 604       fprintf(xyplot, "[%d %d]\n", i, pair_table[i]);
 605   /* add gquad pairs */
 606   ge=0;
 607   while ( (ee=parse_gquad(structure+ge, &Lg, l)) >0 ) {
 608     int k;
 609     fprintf(xyplot, "%% gquad\n");
 610     ge += ee;
 611     gb=ge-Lg*4-l[0]-l[1]-l[2]+1; /* add pseudo-base pair encloding gquad */
 612     for (k=0; k<Lg; k++) {
 613       int ii, jj, il;
 614       for (il=0, ii=gb+k; il<3; il++) {
 615         jj = ii+l[il]+Lg;
 616         fprintf(xyplot, "[%d %d]\n", ii, jj);
 617         ii = jj;
 618       }
 619       jj = gb+k;
 620       fprintf(xyplot, "[%d %d]\n", jj, ii);
 621     }
 622   }
 623
 624   fprintf(xyplot, "] def\n\n");
 625
 626   fprintf(xyplot, "init\n\n");
 627   /* draw the data */
 628   if (pre) {
 629     fprintf(xyplot, "%% Start Annotations\n");
 630     fprintf(xyplot, "%s\n", pre);
 631     fprintf(xyplot, "%% End Annotations\n");
 632   }
 633   fprintf(xyplot,
 634           "%% switch off outline pairs or bases by removing these lines\n"
 635           "drawoutline\n"
 636           "drawpairs\n"
 637           "drawbases\n");
 638
 639   if (post) {
 640     fprintf(xyplot, "%% Start Annotations\n");
 641     fprintf(xyplot, "%s\n", post);
 642     fprintf(xyplot, "%% End Annotations\n");
 643   }
 644   fprintf(xyplot, "%% show it\nshowpage\n");
 645   fprintf(xyplot, "end\n");
 646   fprintf(xyplot, "%%%%EOF\n");
 647
 648   fclose(xyplot);
 649
 650   free(pair_table);
 651   free(pair_table_g);
 652   free(X); free(Y);
 653   return 1; /* success */
 654 }
 655
 656
 657 int PS_rna_plot_snoop_a(char *string, char *structure, char *ssfile, int *relative_access, const char *seqs[])
 658 {
 659   float  xmin, xmax, ymin, ymax, size;
 660   int    i, length;
 661   float *X, *Y;
 662   FILE  *xyplot;
 663   short *pair_table;
 664   short *pair_table_snoop;
 665
 666   length = strlen(string);
 667
 668   xyplot = fopen(ssfile, "w");
 669   if (xyplot == NULL) {
 670     fprintf(stderr, "can't open file %s - not doing xy_plot\n", ssfile);
 671     return 0;
 672   }
 673
 674   pair_table = make_pair_table(structure);
 675   pair_table_snoop = make_pair_table_snoop(structure);
 676
 677   X = (float *) space((length+1)*sizeof(float));
 678   Y = (float *) space((length+1)*sizeof(float));
 679   if (rna_plot_type == 0)
 680     i = simple_xy_coordinates(pair_table, X, Y);
 681   else
 682     i = naview_xy_coordinates(pair_table, X, Y);
 683   if(i!=length) fprintf(stderr,"strange things happening in PS_rna_plot...\n");
 684 /*   printf("cut_point %d\n", cut_point); */
 685   xmin = xmax = X[0];
 686   ymin = ymax = Y[0];
 687 /*   for (i = 1; i < length; i++) { */
 688 /*     printf("%d X %f Y %f \n", i, X[i], Y[i]); */
 689 /*     xmin = X[i] < xmin ? X[i] : xmin; */
 690 /*     xmax = X[i] > xmax ? X[i] : xmax; */
 691 /*     ymin = Y[i] < ymin ? Y[i] : ymin; */
 692 /*     ymax = Y[i] > ymax ? Y[i] : ymax; */
 693 /*   } */
 694   /* localize centre of the interaction bucket. Geometry */
 695
 696   for (i = 1; i < cut_point; i++) {  /* interior loop of size 0 */
 697     if(pair_table_snoop[i] != 0){
 698       X[i-1]=X[pair_table_snoop[i]-1];
 699       Y[i-1]=Y[pair_table_snoop[i]-1];
 700     }
 701     else if(pair_table_snoop[i-1] && pair_table_snoop[i+1]){ /* interior loop of size 1 */
 702       X[i-1]=X[pair_table_snoop[i-1] -1-1];
 703       Y[i-1]=Y[pair_table_snoop[i-1] -1-1];
 704     }
 705     else if(pair_table_snoop[i-1] && pair_table_snoop[i+2]){ /* interior loop of size 2 */
 706       if(pair_table_snoop[i-1] - pair_table_snoop[i+2] ==2){
 707         X[i-1]=X[pair_table_snoop[i-1]-2];
 708         Y[i-1]=Y[pair_table_snoop[i-1]-2];
 709         X[i]=X[pair_table_snoop[i+2]];
 710         Y[i]=Y[pair_table_snoop[i+2]];
 711         i++;
 712       }
 713       else if(pair_table[pair_table_snoop[i-1]-1]){
 714         X[i-1]=X[pair_table_snoop[i-1]-2];
 715         Y[i-1]=Y[pair_table_snoop[i-1]-2];
 716         X[i]=X[pair_table[pair_table_snoop[i-1]-1]-1];
 717         Y[i]=Y[pair_table[pair_table_snoop[i-1]-1]-1];
 718         i++;
 719       }
 720       else if(pair_table[pair_table_snoop[i-1]-2]){
 721         X[i-1]=X[pair_table_snoop[i-1]-3];
 722         Y[i-1]=Y[pair_table_snoop[i-1]-3];
 723         X[i]=X[pair_table[pair_table_snoop[i-1]-2]-1];
 724         Y[i]=Y[pair_table[pair_table_snoop[i-1]-2]-1];
 725         i++;
 726       }
 727       else if(pair_table[pair_table_snoop[i-1]-3]){
 728         X[i-1]=X[pair_table_snoop[i-1]-4];
 729         Y[i-1]=Y[pair_table_snoop[i-1]-4];
 730         X[i]=X[pair_table[pair_table_snoop[i-1]-3]-1];
 731         Y[i]=Y[pair_table[pair_table_snoop[i-1]-3]-1];
 732         i++;
 733       }
 734       else{
 735         X[i-1]=X[pair_table_snoop[i-1]-2];
 736         Y[i-1]=Y[pair_table_snoop[i-1]-2];
 737         X[i]=X[pair_table_snoop[i+2]];
 738         Y[i]=Y[pair_table_snoop[i+2]];
 739         i++;
 740       }
 741     }
 742     else if(pair_table_snoop[i-1] && pair_table_snoop[i+3]){ /* interior loop of size 2 */
 743       if(pair_table[pair_table_snoop[i-1]-1]){
 744         X[i-1]=0.5*(X[pair_table_snoop[i-1]-1]+X[pair_table_snoop[i-1]-2]);
 745         Y[i-1]=0.5*(Y[pair_table_snoop[i-1]-1]+Y[pair_table_snoop[i-1]-2]);
 746         X[i]=  0.5*(X[pair_table[pair_table_snoop[i-1]-1]-1]+X[pair_table_snoop[i-1]-2]);
 747         Y[i]=  0.5*(Y[pair_table[pair_table_snoop[i-1]-1]-1]+Y[pair_table_snoop[i-1]-2]);
 748         X[i+1]=0.5*(X[pair_table[pair_table_snoop[i-1]-1]-2]+X[pair_table[pair_table_snoop[i-1]-1]-1]);
 749         Y[i+1]=0.5*(Y[pair_table[pair_table_snoop[i-1]-1]-2]+Y[pair_table[pair_table_snoop[i-1]-1]-1]);
 750         i++;i++;
 751
 752       }
 753       else if(pair_table[pair_table_snoop[i-1]-2]){
 754         X[i-1]=0.5*(X[pair_table_snoop[i-1]-2]+X[pair_table_snoop[i-1]-3]);
 755         Y[i-1]=0.5*(Y[pair_table_snoop[i-1]-2]+Y[pair_table_snoop[i-1]-3]);
 756         X[i]=  0.5*(X[pair_table[pair_table_snoop[i-1]-2]-1]+X[pair_table_snoop[i-1]-3]);
 757         Y[i]=  0.5*(Y[pair_table[pair_table_snoop[i-1]-2]-1]+Y[pair_table_snoop[i-1]-3]);
 758         X[i+1]=0.5*(X[pair_table[pair_table_snoop[i-1]-2]-2]+X[pair_table[pair_table_snoop[i-1]-2]-1]);
 759         Y[i+1]=0.5*(Y[pair_table[pair_table_snoop[i-1]-2]-2]+Y[pair_table[pair_table_snoop[i-1]-2]-1]);
 760         i++;i++;
 761       }
 762       else if(pair_table[pair_table_snoop[i-1]-3]){
 763         X[i-1]=0.5*(X[pair_table_snoop[i-1]-3]+X[pair_table_snoop[i-1]-4]);
 764         Y[i-1]=0.5*(Y[pair_table_snoop[i-1]-3]+Y[pair_table_snoop[i-1]-4]);
 765         X[i]=  0.5*(X[pair_table[pair_table_snoop[i-1]-3]-1]+X[pair_table_snoop[i-1]-4]);
 766         Y[i]=  0.5*(Y[pair_table[pair_table_snoop[i-1]-3]-1]+Y[pair_table_snoop[i-1]-4]);
 767         X[i+1]=0.5*(X[pair_table[pair_table_snoop[i-1]-3]-2]+X[pair_table[pair_table_snoop[i-1]-3]-1]);
 768         Y[i+1]=0.5*(Y[pair_table[pair_table_snoop[i-1]-3]-2]+Y[pair_table[pair_table_snoop[i-1]-3]-1]);
 769         i++;i++;
 770       }
 771       else{
 772         X[i-1]=X[pair_table_snoop[i-1]-2];
 773         Y[i-1]=Y[pair_table_snoop[i-1]-2];
 774         X[i]=X[pair_table_snoop[i-1]-2];
 775         Y[i]=Y[pair_table_snoop[i-1]-2];
 776         X[i+1]=X[pair_table_snoop[i-1]-2];
 777         Y[i+1]=Y[pair_table_snoop[i-1]-2];
 778         i++;i++;
 779       }
 780     }
 781   }
 782   double xC;
 783   double yC;
 784   double R ;
 785   double d ;
 786   float X0=-1,Y0=-1,X1=-1,Y1=-1,X2=-1,Y2=-1;
 787 /*   int c1,c2,c3; */
 788   for(i=1;i<cut_point; i++){
 789     if(pair_table_snoop[i]){
 790       X0=X[pair_table_snoop[i]-1];Y0=Y[pair_table_snoop[i]-1];
 791   /*     c1=pair_table_snoop[i]; */
 792       i++;
 793       break;
 794     }
 795   }
 796   for(;i<cut_point; i++){
 797     if(pair_table_snoop[i]){
 798       X1=X[pair_table_snoop[i]-1];Y1=Y[pair_table_snoop[i]-1];
 799     /*   c2=pair_table_snoop[i]; */
 800       i++;
 801       break;
 802     }
 803   }
 804   for(;i<cut_point; i++){
 805     if(pair_table_snoop[i]){
 806       X2=X[pair_table_snoop[i]-1];Y2=Y[pair_table_snoop[i]-1];
 807     /*   c3=pair_table_snoop[i]; */
 808       i++;
 809       break;
 810     }
 811   }
 812 /*   for(i=cut_point-2;i>pair_table_snoop[c1]; i--){ */
 813 /*     if(pair_table_snoop[i]){ */
 814 /*       X1=X[pair_table_snoop[i]-1];Y1=Y[pair_table_snoop[i]-1]; */
 815 /*       c2=pair_table_snoop[i]; */
 816 /*       i++; */
 817 /*       break; */
 818 /*     } */
 819 /*   } */
 820 /*   for(i=pair_table_snoop[c1]+1;i<pair_table_snoop[c2]; i++){ */
 821 /*     if(pair_table_snoop[i]){ */
 822 /*       X2=X[pair_table_snoop[i]-1];Y2=Y[pair_table_snoop[i]-1]; */
 823 /*       c3=pair_table_snoop[i]; */
 824 /*       i++; */
 825 /*       break; */
 826 /*     } */
 827 /*   } */
 828  if(X0 < 0 || X1 < 0 || X2 < 0){
 829    printf("Could not get the center of the binding bucket. No ps file will be produced!\n");
 830    fclose(xyplot);
 831    free(pair_table);
 832    free(pair_table_snoop);
 833    free(X);free(Y);
 834    pair_table=NULL;pair_table_snoop=NULL;X=NULL;Y=NULL;
 835    return 0;
 836  }
 837   double alpha   =   (X0 -X1)/(Y1-Y0);
 838   double alpha_p =   (X1 -X2)/(Y2-Y1);
 839   double b =         (Y0+Y1 -alpha*(X0+X1))*0.5;
 840   double b_p =       (Y1+Y2 -alpha_p*(X1+X2))*0.5;
 841   /*    if(abs(alpha -alpha_p) > 0.0000001){ */
 842   xC  =  (b_p - b) / (alpha - alpha_p);
 843   yC  =  alpha * xC + b;
 844   R   =  sqrt((xC-X0)*(xC-X0) + (yC-Y0)*(yC-Y0));
 845   d   = sqrt((X0-X1)*(X0-X1) + (Y0-Y1)*(Y0-Y1));
 846    for (i = 1; i < cut_point; i++) {
 847      X[i-1] = X[i-1] + 0.25*(xC-X[i-1]);
 848      Y[i-1] = Y[i-1] + 0.25*(yC-Y[i-1]);
 849    }
 850   size = MAX2((xmax-xmin),(ymax-ymin));
 851   fprintf(xyplot,
 852           "%%!PS-Adobe-3.0 EPSF-3.0\n"
 853           "%%%%Creator: %s, ViennaRNA-%s\n"
 854           "%%%%CreationDate: %s"
 855           "%%%%Title: RNA Secondary Structure Plot\n"
 856           "%%%%BoundingBox: 66 210 518 662\n"
 857           "%%%%DocumentFonts: Helvetica\n"
 858           "%%%%Pages: 1\n"
 859           "%%%%EndComments\n\n"
 860           "%%Options: %s\n", rcsid+5, VERSION, time_stamp(), option_string());
 861   fprintf(xyplot, "%% to switch off outline pairs of sequence comment or\n"
 862           "%% delete the appropriate line near the end of the file\n\n");
 863   fprintf(xyplot, "%s", RNAss_head);
 864   char **A;
 865   fprintf(xyplot, "%s", anote_macros);
 866   if(seqs){
 867     fprintf(xyplot, "%s", anote_macros);
 868     A = annote(structure, (const char**) seqs);
 869   }
 870   fprintf(xyplot, "%%%%EndProlog\n");
 871
 872   fprintf(xyplot, "RNAplot begin\n"
 873           "%% data start here\n");
 874   /* cut_point */
 875   if (cut_point > 0 && cut_point <= strlen(string))
 876     fprintf(xyplot, "/cutpoint %d def\n", cut_point-1);
 877   /* sequence */
 878   fprintf(xyplot,"/sequence (\\\n");
 879   i=0;
 880   while (i<length) {
 881     fprintf(xyplot, "%.255s\\\n", string+i);  /* no lines longer than 255 */
 882     i+=255;
 883   }
 884   fprintf(xyplot,") def\n");
 885   /* coordinates */
 886   fprintf(xyplot, "/coor [\n");
 887   for (i = 0; i < length; i++)
 888     fprintf(xyplot, "[%3.3f %3.3f]\n", X[i], Y[i]);
 889   fprintf(xyplot, "] def\n");
 890   /* base pairs */
 891   fprintf(xyplot, "/pairs [\n");
 892   for (i = 1; i <= length; i++)
 893     if (pair_table[i]>i)
 894       fprintf(xyplot, "[%d %d]\n", i, pair_table[i]);
 895   for (i = 1; i <= length; i++)
 896     if (pair_table_snoop[i]>i)
 897       fprintf(xyplot, "[%d %d]\n", i, pair_table_snoop[i]);
 898   fprintf(xyplot, "] def\n\n");
 899   if(relative_access){
 900     fprintf(xyplot,"/S [\n");
 901     for(i=0;i<cut_point-1; i++){
 902       fprintf(xyplot, " %f\n", (float)relative_access[i]/100);
 903     }
 904     fprintf(xyplot,"]\n bind def\n");
 905     fprintf(xyplot,"/invert false def\n");
 906     fprintf(xyplot,"/range 0.8 def\n");
 907     fprintf(xyplot,"/drawreliability {\n"
 908                    "/Smax 2.6 def\n"
 909                    "  0        \n"
 910                    "  coor 0 cutpoint getinterval {\n"
 911                    "    aload pop\n"
 912                    "    S 3 index get\n"
 913                    "    Smax div range mul\n"
 914                    "    invert {range exch sub} if\n"
 915                    "    1 1 sethsbcolor\n"
 916                    "    newpath\n"
 917                    "    fsize 2.5 div 0 360 arc\n"
 918                    "    fill\n"
 919                    "    1 add\n"
 920                    "  } forall\n"
 921                    "\n"
 922                    "} bind def\n");
 923   }
 924   fprintf(xyplot, "init\n\n");
 925   /*raw the data */
 926   if (seqs) {
 927      fprintf(xyplot, "%% Start Annotations\n");
 928      fprintf(xyplot, "%s\n", A[0]);
 929      fprintf(xyplot, "%% End Annotations\n");
 930    }
 931
 932
 933   fprintf(xyplot,"%%switch off outline pairs or bases by removing these lines\n");
 934   if(relative_access){
 935     fprintf(xyplot,"drawreliability\n");
 936   }
 937   fprintf(xyplot,
 938           "drawoutline\n"
 939           "drawpairs\n"
 940           "drawbases\n");
 941   /* fprintf(xyplot, "%d cmark\n",c1); */
 942   /* fprintf(xyplot, "%d cmark\n",c2); */
 943   /* fprintf(xyplot, "%d cmark\n",c3); */
 944   if (seqs) {
 945      fprintf(xyplot, "%% Start Annotations\n");
 946      fprintf(xyplot, "%s\n", A[1]);
 947      fprintf(xyplot, "%% End Annotations\n");
 948    }
 949   fprintf(xyplot, "%% show it\nshowpage\n");
 950   fprintf(xyplot, "end\n");
 951   fprintf(xyplot, "%%%%EOF\n");
 952
 953   fclose(xyplot);
 954   if(seqs){free(A[0]);free(A[1]);free(A);}
 955   free(pair_table);free(pair_table_snoop);
 956   free(X); free(Y);
 957   return 1; /* success */
 958 }
 959
 960
 961 PRIVATE char **annote(const char *structure, const char *AS[]) {
 962   char *ps, *colorps, **A;
 963   int i, n, s, pairings, maxl;
 964   short *ptable;
 965   char * colorMatrix[6][3] = {
 966     {"0.0 1", "0.0 0.6",  "0.0 0.2"},  /* red    */
 967     {"0.16 1","0.16 0.6", "0.16 0.2"}, /* ochre  */
 968     {"0.32 1","0.32 0.6", "0.32 0.2"}, /* turquoise */
 969     {"0.48 1","0.48 0.6", "0.48 0.2"}, /* green  */
 970     {"0.65 1","0.65 0.6", "0.65 0.2"}, /* blue   */
 971     {"0.81 1","0.81 0.6", "0.81 0.2"} /* violet */
 972   };
 973
 974   make_pair_matrix();
 975   n = strlen(AS[0]);
 976   maxl = 1024;
 977
 978   A = (char **) space(sizeof(char *)*2);
 979   ps = (char *) space(maxl);
 980   colorps = (char *) space(maxl);
 981   ptable = alimake_pair_table(structure);
 982   for (i=1; i<=n; i++) {
 983     char pps[64], ci='\0', cj='\0';
 984     int j, type, pfreq[8] = {0,0,0,0,0,0,0,0}, vi=0, vj=0;
 985     if ((j=ptable[i])<i) continue;
 986     for (s=0; AS[s]!=NULL; s++) {
 987       type = pair[encode_char(AS[s][i-1])][encode_char(AS[s][j-1])];
 988       pfreq[type]++;
 989       if (type) {
 990         if (AS[s][i-1] != ci) { ci = AS[s][i-1]; vi++;}
 991         if (AS[s][j-1] != cj) { cj = AS[s][j-1]; vj++;}
 992       }
 993     }
 994     for (pairings=0,s=1; s<=7; s++) {
 995       if (pfreq[s]) pairings++;
 996     }
 997
 998     if ((maxl - strlen(ps) < 192) || ((maxl - strlen(colorps)) < 64)) {
 999       maxl *= 2;
1000       ps = realloc(ps, maxl);
1001       colorps = realloc(colorps, maxl);
1002       if ((ps==NULL) || (colorps == NULL))
1003           nrerror("out of memory in realloc");
1004     }
1005
1006     if (pfreq[0]<=2) {
1007       snprintf(pps, 64, "%d %d %s colorpair\n",
1008                i,j, colorMatrix[pairings-1][pfreq[0]]);
1009       strcat(colorps, pps);
1010     }
1011
1012     if (pfreq[0]>0) {
1013       snprintf(pps, 64, "%d %d %d gmark\n", i, j, pfreq[0]);
1014       strcat(ps, pps);
1015     }
1016     if (vi>1) {
1017       snprintf(pps, 64, "%d cmark\n", i);
1018       strcat(ps, pps);
1019     }
1020     if (vj>1) {
1021       snprintf(pps, 64, "%d cmark\n", j);
1022       strcat(ps, pps);
1023     }
1024   }
1025   free(ptable);
1026   A[0]=colorps;
1027   A[1]=ps;
1028   return A;
1029 }
1030
1031 /*--------------------------------------------------------------------------*/
1032
1033
1034 int svg_rna_plot(char *string, char *structure, char *ssfile)
1035 {
1036   float  xmin, xmax, ymin, ymax, size;
1037   int    i, length;
1038   float *X, *Y, *R = NULL, *CX = NULL, *CY = NULL;
1039   FILE  *xyplot;
1040   short *pair_table;
1041
1042   length = strlen(string);
1043
1044   xyplot = fopen(ssfile, "w");
1045   if (xyplot == NULL) {
1046     fprintf(stderr, "can't open file %s - not doing xy_plot\n", ssfile);
1047     return 0;
1048   }
1049
1050   pair_table = make_pair_table(structure);
1051
1052   X = (float *) space((length+1)*sizeof(float));
1053   Y = (float *) space((length+1)*sizeof(float));
1054
1055   switch(rna_plot_type){
1056     case VRNA_PLOT_TYPE_SIMPLE:   i = simple_xy_coordinates(pair_table, X, Y);
1057                                   break;
1058     case VRNA_PLOT_TYPE_CIRCULAR: {
1059                                     int radius = 3*length;
1060                                     int dr = 0;
1061                                     R = (float *) space((length+1)*sizeof(float));
1062                                     CX = (float *) space((length+1)*sizeof(float));
1063                                     CY = (float *) space((length+1)*sizeof(float));
1064                                     i = simple_circplot_coordinates(pair_table, X, Y);
1065                                     for (i = 0; i < length; i++) {
1066                                       if(i+1 < pair_table[i+1]){
1067                                         dr = (pair_table[i+1]-i+1 <= (length/2 + 1)) ? pair_table[i+1]-i : i + length - pair_table[i+1];
1068                                         R[i] = 1. - (2.*dr/(float)length);
1069                                       }
1070                                       else if(pair_table[i+1]){
1071                                         R[i] = R[pair_table[i+1]-1];
1072                                       }
1073                                       else{
1074                                         R[i] = 1.0;
1075                                       }
1076                                       CX[i] = X[i] * radius * R[i] + radius;
1077                                       CY[i] = Y[i] * radius * R[i] + radius;
1078                                       X[i] *= radius;
1079                                       X[i] += radius;
1080                                       Y[i] *= radius;
1081                                       Y[i] += radius;
1082                                     }
1083                                   }
1084                                   break;
1085     default:                      i = naview_xy_coordinates(pair_table, X, Y);
1086                                   break;
1087   }
1088
1089   if(i!=length) fprintf(stderr,"strange things happening in PS_rna_plot...\n");
1090
1091
1092   xmin = xmax = X[0];
1093   ymin = ymax = Y[0];
1094   for (i = 1; i < length; i++) {
1095      xmin = X[i] < xmin ? X[i] : xmin;
1096      xmax = X[i] > xmax ? X[i] : xmax;
1097      ymin = Y[i] < ymin ? Y[i] : ymin;
1098      ymax = Y[i] > ymax ? Y[i] : ymax;
1099   }
1100   for (i = 0; i < length; i++)
1101     Y[i] = ymin+ymax - Y[i]; /* mirror coordinates so they look as in PS */
1102
1103   if(rna_plot_type == VRNA_PLOT_TYPE_CIRCULAR)
1104     for (i = 0; i < length; i++){
1105       CY[i] = ymin+ymax - CY[i]; /* mirror coordinates so they look as in PS */
1106     }
1107
1108   size = MAX2((xmax-xmin),(ymax-ymin));
1109   size += 15; /* add some so the bounding box isn't too tight */
1110
1111   fprintf(xyplot,
1112           "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"yes\"?>\n"
1113           "<svg xmlns=\"http://www.w3.org/2000/svg\" height=\"452\" width=\"452\">\n");
1114   fprintf(xyplot,
1115           "<script type=\"text/ecmascript\">\n"
1116           "      <![CDATA[\n"
1117           "        var shown = 1;\n"
1118           "        function click() {\n"
1119           "             var seq = document.getElementById(\"seq\");\n"
1120           "             if (shown==1) {\n"
1121           "               seq.setAttribute(\"style\", \"visibility: hidden\");\n"
1122           "               shown = 0;\n"
1123           "             } else {\n"
1124           "               seq.setAttribute(\"style\", \"visibility: visible\");\n"
1125           "               shown = 1;\n"
1126           "             }\n"
1127           "         }\n"
1128           "        ]]>\n"
1129           "</script>\n");
1130   fprintf(xyplot,
1131           "  <rect style=\"stroke: white; fill: white\" height=\"452\" x=\"0\" y=\"0\" width=\"452\" onclick=\"click(evt)\" />\n"
1132           "  <g transform=\"scale(%7f,%7f) translate(%7f,%7f)\">\n",
1133           SIZE/size, SIZE/size, (size-xmin-xmax)/2, (size-ymin-ymax)/2);
1134
1135   fprintf(xyplot,
1136           "    <polyline style=\"stroke: black; fill: none; stroke-width: 1.5\" id=\"outline\" points=\"\n");
1137   for (i = 0; i < length; i++)
1138     fprintf(xyplot, "      %3.3f,%3.3f\n", X[i], Y[i]);
1139   fprintf(xyplot,"    \" />\n");
1140
1141   fprintf(xyplot,"    <g style=\"stroke: black; stroke-width: 1; fill: none;\" id=\"pairs\">\n");
1142   for (i = 1; i <= length; i++) {
1143     int j;
1144     if ((j=pair_table[i])>i){
1145       if(rna_plot_type == VRNA_PLOT_TYPE_CIRCULAR)
1146         fprintf(xyplot,
1147                 "      <path id=\"%d,%d\" d=\"M %6.15f %6.15f C %6.15f,%6.15f %6.15f,%6.15f %6.15f %6.15f\" />\n",
1148                 i,j, X[i-1], Y[i-1], CX[i-1], CY[i-1], CX[j-1], CY[j-1], X[j-1], Y[j-1]);
1149       else
1150         fprintf(xyplot,
1151                 "      <line id=\"%d,%d\" x1=\"%6.5f\" y1=\"%6.5f\" x2=\"%6.5f\" y2=\"%6.5f\" />\n",
1152                 i,j, X[i-1], Y[i-1], X[j-1], Y[j-1]);
1153     }
1154   }
1155   fprintf(xyplot, "    </g>\n");
1156   fprintf(xyplot, "    <g style=\"font-family: SansSerif\" transform=\"translate(-4.6, 4)\" id=\"seq\">\n");
1157   for (i = 0; i < length; i++)
1158     fprintf(xyplot, "      <text x=\"%.3f\" y=\"%.3f\">%c</text>\n", X[i], Y[i], string[i]);
1159   fprintf(xyplot, "    </g>\n");
1160   fprintf(xyplot, "  </g>\n");
1161   fprintf(xyplot, "</svg>\n");
1162
1163   fclose(xyplot);
1164
1165   free(pair_table);
1166   free(X); free(Y);
1167   if(R) free(R);
1168   if(CX) free(CX);
1169   if(CY) free(CY);
1170   return 1; /* success */
1171 }
1172
1173 /*--------------------------------------------------------------------------*/
1174
1175 PUBLIC int ssv_rna_plot(char *string, char *structure, char *ssfile)
1176 {           /* produce input for the SStructView java applet */
1177   FILE *ssvfile;
1178   int i, bp;
1179   int length;
1180   short *pair_table;
1181   float *X, *Y;
1182   float xmin, xmax, ymin, ymax;
1183
1184   ssvfile = fopen(ssfile, "w");
1185   if (ssvfile == NULL) {
1186      fprintf(stderr, "can't open file %s - not doing xy_plot\n", ssfile);
1187      return 0;
1188   }
1189   length = strlen(string);
1190   pair_table = make_pair_table(structure);
1191
1192   /* make coordinates */
1193   X = (float *) space((length+1)*sizeof(float));
1194   Y = (float *) space((length+1)*sizeof(float));
1195
1196   if (rna_plot_type == 0)
1197     i = simple_xy_coordinates(pair_table, X, Y);
1198   else
1199     i = naview_xy_coordinates(pair_table, X, Y);
1200   if (i!=length)
1201     fprintf(stderr,"strange things happening in ssv_rna_plot...\n");
1202
1203   /* make coords nonegative */
1204   xmin = xmax = X[0];
1205   ymin = ymax = Y[0];
1206   for (i = 1; i < length; i++) {
1207      xmin = X[i] < xmin ? X[i] : xmin;
1208      xmax = X[i] > xmax ? X[i] : xmax;
1209      ymin = Y[i] < ymin ? Y[i] : ymin;
1210      ymax = Y[i] > ymax ? Y[i] : ymax;
1211   }
1212   if (xmin<1) {
1213     for (i = 0; i <= length; i++)
1214       X[i] -= xmin-1;
1215     xmin = 1;
1216   }
1217   if (ymin<1) {
1218     for (i = 0; i <= length; i++)
1219       Y[i] -= ymin-1;
1220     ymin = 1;
1221   }
1222 #if 0
1223   {
1224     float size, xoff, yoff;
1225     float JSIZE = 500; /* size of the java applet window */
1226     /* rescale coordinates, center on square of size HSIZE */
1227     size = MAX2((xmax-xmin),(ymax-ymin));
1228     xoff = (size - xmax + xmin)/2;
1229     yoff = (size - ymax + ymin)/2;
1230     for (i = 0; i <= length; i++) {
1231       X[i] = (X[i]-xmin+xoff)*(JSIZE-10)/size + 5;
1232       Y[i] = (Y[i]-ymin+yoff)*(JSIZE-10)/size + 5;
1233     }
1234   }
1235 #endif
1236   /* */
1237
1238   fprintf(ssvfile,
1239           "# Vienna RNA Package %s\n"
1240           "# SStructView Output\n"
1241           "# CreationDate: %s\n"
1242           "# Name: %s\n"
1243           "# Options: %s\n", VERSION, time_stamp(), ssfile, option_string());
1244   for (i=1; i<=length; i++)
1245     fprintf(ssvfile, "BASE\t%d\t%c\t%d\t%d\n",
1246             i, string[i-1], (int) (X[i-1]+0.5), (int) (Y[i-1]+0.5));
1247   for (bp=1, i=1; i<=length; i++)
1248     if (pair_table[i]>i)
1249       fprintf(ssvfile, "BASE-PAIR\tbp%d\t%d\t%d\n", bp++, i, pair_table[i]);
1250   fclose(ssvfile);
1251
1252   free(pair_table);
1253   free(X); free(Y);
1254   return 1; /* success */
1255 }
1256
1257 /*---------------------------------------------------------------------------*/
1258 PUBLIC int xrna_plot(char *string, char *structure, char *ssfile)
1259 {           /* produce input for XRNA RNA drawing program */
1260   FILE *ss_file;
1261   int i;
1262   int length;
1263   short *pair_table;
1264   float *X, *Y;
1265
1266   ss_file = fopen(ssfile, "w");
1267   if (ss_file == NULL) {
1268     fprintf(stderr, "can't open file %s - not doing xy_plot\n", ssfile);
1269     return 0;
1270   }
1271
1272   length = strlen(string);
1273   pair_table = make_pair_table(structure);
1274
1275   /* make coordinates */
1276   X = (float *) space((length+1)*sizeof(float));
1277   Y = (float *) space((length+1)*sizeof(float));
1278
1279   if (rna_plot_type == 0)
1280     i = simple_xy_coordinates(pair_table, X, Y);
1281   else
1282     i = naview_xy_coordinates(pair_table, X, Y);
1283   if (i!=length)
1284     fprintf(stderr,"strange things happening in xrna_plot...\n");
1285
1286   fprintf(ss_file,
1287           "# Vienna RNA Package %s, XRNA output\n"
1288           "# CreationDate: %s\n"
1289           "# Options: %s\n", VERSION, time_stamp(), option_string());
1290   for (i=1; i<=length; i++)
1291     /* XRNA likes to have coordinate mirrored, so we use (-X, Y) */
1292     fprintf(ss_file, "%d %c %6.2f %6.2f %d %d\n", i, string[i-1],
1293             -X[i-1], Y[i-1], (pair_table[i]?1:0), pair_table[i]);
1294   fclose(ss_file);
1295
1296   free(pair_table);
1297   free(X); free(Y);
1298   return 1; /* success */
1299 }
1300
1301
1302 /*---------------------------------------------------------------------------*/
1303 const char *RNAdp_prolog =
1304 "%This file contains the square roots of the base pair probabilities in the form\n"
1305 "% i  j  sqrt(p(i,j)) ubox\n\n"
1306 "%%BeginProlog\n"
1307 "/DPdict 100 dict def\n"
1308 "DPdict begin\n"
1309 "/logscale false def\n"
1310 "/lpmin 1e-05 log def\n\n"
1311 "/box { %size x y box - draws box centered on x,y\n"
1312 "   2 index 0.5 mul sub            % x -= 0.5\n"
1313 "   exch 2 index 0.5 mul sub exch  % y -= 0.5\n"
1314 "   3 -1 roll dup rectfill\n"
1315 "} bind def\n\n"
1316 "/ubox {\n"
1317 "   logscale {\n"
1318 "      log dup add lpmin div 1 exch sub dup 0 lt { pop 0 } if\n"
1319 "   } if\n"
1320 "   3 1 roll\n"
1321 "   exch len exch sub 1 add box\n"
1322 "} bind def\n\n"
1323 "/lbox {\n"
1324 "   3 1 roll\n"
1325 "   len exch sub 1 add box\n"
1326 "} bind def\n\n"
1327 "/drawseq {\n"
1328 "% print sequence along all 4 sides\n"
1329 "[ [0.7 -0.3 0 ]\n"
1330 "  [0.7 0.7 len add 0]\n"
1331 "  [-0.3 len sub -0.4 -90]\n"
1332 "  [-0.3 len sub 0.7 len add -90]\n"
1333 "] {\n"
1334 "   gsave\n"
1335 "    aload pop rotate translate\n"
1336 "    0 1 len 1 sub {\n"
1337 "     dup 0 moveto\n"
1338 "     sequence exch 1 getinterval\n"
1339 "     show\n"
1340 "    } for\n"
1341 "   grestore\n"
1342 "  } forall\n"
1343 "} bind def\n\n"
1344 "/drawgrid{\n"
1345 "  0.01 setlinewidth\n"
1346 "  len log 0.9 sub cvi 10 exch exp  % grid spacing\n"
1347 "  dup 1 gt {\n"
1348 "     dup dup 20 div dup 2 array astore exch 40 div setdash\n"
1349 "  } { [0.3 0.7] 0.1 setdash } ifelse\n"
1350 "  0 exch len {\n"
1351 "     dup dup\n"
1352 "     0 moveto\n"
1353 "     len lineto \n"
1354 "     dup\n"
1355 "     len exch sub 0 exch moveto\n"
1356 "     len exch len exch sub lineto\n"
1357 "     stroke\n"
1358 "  } for\n"
1359 "  [] 0 setdash\n"
1360 "  0.04 setlinewidth \n"
1361 "  currentdict /cutpoint known {\n"
1362 "    cutpoint 1 sub\n"
1363 "    dup dup -1 moveto len 1 add lineto\n"
1364 "    len exch sub dup\n"
1365 "    -1 exch moveto len 1 add exch lineto\n"
1366 "    stroke\n"
1367 "  } if\n"
1368 "  0.5 neg dup translate\n"
1369 "} bind def\n\n"
1370 "end\n"
1371 "%%EndProlog\n";
1372
1373
1374 /*---------------------------------------------------------------------------*/
1375 int PS_color_dot_plot(char *seq, cpair *pi, char *wastlfile) {
1376   /* produce color PostScript dot plot from cpair */
1377
1378   FILE *wastl;
1379   int i, length;
1380
1381   length= strlen(seq);
1382   wastl = PS_dot_common(seq, wastlfile, NULL, 0);
1383   if (wastl==NULL)  return 0; /* return 0 for failure */
1384
1385   fprintf(wastl, "/hsb {\n"
1386           "dup 0.3 mul 1 exch sub sethsbcolor\n"
1387           "} bind def\n\n");
1388
1389   fprintf(wastl, "\n%%draw the grid\ndrawgrid\n\n");
1390   fprintf(wastl,"%%start of base pair probability data\n");
1391
1392   /* print boxes */
1393    i=0;
1394    while (pi[i].j>0) {
1395      fprintf(wastl,"%1.2f %1.2f hsb %d %d %1.6f ubox\n",
1396              pi[i].hue, pi[i].sat, pi[i].i, pi[i].j, sqrt(pi[i].p));
1397
1398      if (pi[i].mfe)
1399        fprintf(wastl,"%1.2f %1.2f hsb %d %d %1.4f lbox\n",
1400                pi[i].hue, pi[i].sat, pi[i].i, pi[i].j, pi[i].p);
1401      i++;
1402    }
1403
1404    fprintf(wastl,"showpage\n"
1405            "end\n"
1406            "%%%%EOF\n");
1407    fclose(wastl);
1408    return 1; /* success */
1409 }
1410
1411 /*---------------------------------------------------------------------------*/
1412
1413 const char *RNAdp_gquad_triangle =
1414 "/min { 2 copy gt { exch } if pop } bind def\n\n"
1415 "/utri{ % i j prob utri\n"
1416 "  gsave\n"
1417 "  1 min 2 div\n"
1418 "  0.85 mul 0.15 add 0.95  0.33\n"
1419 "  3 1 roll % prepare hsb color\n"
1420 "  sethsbcolor\n"
1421 "  % now produce the coordinates for lines\n"
1422 "  exch 1 sub dup len exch sub dup 4 -1 roll dup 3 1 roll dup len exch sub\n"
1423 "  moveto lineto lineto closepath fill\n"
1424 "  grestore\n"
1425 "} bind def\n";
1426
1427
1428 static int sort_plist_by_type_desc(const void *p1, const void *p2){
1429   if(((plist*)p1)->type > ((plist*)p2)->type) return -1;
1430   if(((plist*)p1)->type < ((plist*)p2)->type) return 1;
1431   return 0;
1432 }
1433
1434 static int sort_plist_by_prob_asc(const void *p1, const void *p2){
1435   if(((plist*)p1)->p > ((plist*)p2)->p) return 1;
1436   if(((plist*)p1)->p < ((plist*)p2)->p) return -1;
1437   return 0;
1438 }
1439
1440
1441
1442 PUBLIC int PS_dot_plot_list(char *seq,
1443                             char *wastlfile,
1444                             plist *pl,
1445                             plist *mf,
1446                             char *comment){
1447
1448   FILE *wastl;
1449   int length, pl_size, gq_num;
1450   double tmp;
1451   plist *pl1, *pl_tmp;
1452
1453   length= strlen(seq);
1454   wastl = PS_dot_common(seq, wastlfile, comment, 0);
1455   if (wastl==NULL) return 0; /* return 0 for failure */
1456
1457   fprintf(wastl, "%s\n", RNAdp_gquad_triangle);
1458
1459   fprintf(wastl,"%%data starts here\n");
1460
1461   /* sort the plist to bring all gquad triangles to the front */
1462   for(gq_num = pl_size = 0, pl1 = pl; pl1->i > 0; pl1++, pl_size++)
1463     if(pl1->type == 1) gq_num++;
1464   qsort(pl, pl_size, sizeof(plist), sort_plist_by_type_desc);
1465   /* sort all gquad triangles by probability to bring lower probs to the front */
1466   qsort(pl, gq_num, sizeof(plist), sort_plist_by_prob_asc);
1467
1468   /* print triangles for g-quadruplexes in upper half */
1469   fprintf(wastl,"\n%%start of quadruplex data\n");
1470   for (pl1=pl; pl1->type == 1; pl1++) {
1471     tmp = sqrt(pl1->p);
1472     fprintf(wastl, "%d %d %1.9f utri\n", pl1->i, pl1->j, tmp);
1473   }
1474
1475   fprintf(wastl, "\n%%draw the grid\ndrawgrid\n\n");
1476   fprintf(wastl,"%%start of base pair probability data\n");
1477
1478   /* print boxes in upper right half*/
1479   for (; pl1->i>0; pl1++) {
1480     tmp = sqrt(pl1->p);
1481     if(pl1->type == 0)
1482         fprintf(wastl,"%d %d %1.9f ubox\n", pl1->i, pl1->j, tmp);
1483   }
1484
1485
1486   /* print boxes in lower left half (mfe) */
1487   for (pl1=mf; pl1->i>0; pl1++) {
1488     tmp = sqrt(pl1->p);
1489     fprintf(wastl,"%d %d %1.7f lbox\n", pl1->i, pl1->j, tmp);
1490   }
1491
1492   fprintf(wastl,"showpage\n"
1493           "end\n"
1494           "%%%%EOF\n");
1495   fclose(wastl);
1496   return 1; /* success */
1497 }
1498
1499 const char *RNAdp_prolog_turn =
1500 "/drawseq_turn {"
1501 "% print sequence at bottom\n"
1502 "   gsave\n"
1503 "   len 2 sqrt div dup neg 0.28 add exch 0.78 sub translate\n"
1504 "    0 1 len 1 sub {\n"
1505 "     dup dup 2 sqrt mul 0 moveto\n"
1506 "     sequence exch 1 getinterval\n"
1507 "     show\n"
1508 "    } for\n"
1509 "   grestore\n"
1510 "} bind def\n"
1511 "/drawgrid_turn{\n"
1512 "  0.01 setlinewidth\n"
1513 "  len log 0.9 sub cvi 10 exch exp  % grid spacing\n"
1514 "  dup 1 gt {\n"
1515 "     dup dup 20 div dup 2 array astore exch 40 div setdash\n"
1516 "  } { [0.3 0.7] 0.1 setdash } ifelse\n"
1517 "  0 exch len {    %for (0, gridspacing, len) \n"
1518 "     dup dup      %duplicate what - gridspacing??\n"
1519 "     dup len exch sub moveto     %moveto diagonal?\n"
1520 "     dup winSize gt\n"
1521 "     {dup dup len exch sub winSize add lineto}\n"
1522 "     {dup len lineto}ifelse\n"
1523 "     dup len exch sub moveto  %moveto diagonal?\n"
1524 "     dup len winSize sub le\n"
1525 "     {dup dup len exch sub dup winSize exch sub len add exch lineto}\n"
1526 "     {dup dup len exch sub len exch lineto}ifelse"
1527 "     stroke pop pop\n"
1528 "  } for\n"
1529 "  len log 0.9 sub cvi 10 exch exp  % grid spacing\n"
1530 "      dup 1 gt {\n"
1531 "          dup dup 20 div dup 2 array astore exch 40 div setdash\n"
1532 "      } { [0.3 0.7] 0.1 setdash } ifelse\n"
1533 "      0 exch len {    %for (0, gridspacing, len) \n"
1534 "     dup dup      %duplicate what - gridspacing??\n"
1535 "     dup len exch sub moveto     %moveto diagonal?\n"
1536 "     len exch sub 0.7 sub exch 0.7 sub exch lineto\n"
1537 "     stroke\n"
1538 "   }for\n"
1539 " winSize len moveto  len winSize  lineto stroke\n"
1540 "  [] 0 setdash\n"
1541 "  0.04 setlinewidth \n"
1542 "  currentdict /cutpoint known {\n"
1543 "    cutpoint 1 sub\n"
1544 "    dup dup -1 moveto len 1 add lineto\n"
1545 "    len exch sub dup\n"
1546 "    -1 exch moveto len 1 add exch lineto\n"
1547 "   stroke\n"
1548 "  } if\n"
1549 "  0.5 neg dup translate\n"
1550   "} bind def \n\n";
1551
1552 int PS_color_dot_plot_turn(char *seq, cpair *pi, char *wastlfile, int winSize) {
1553   /* produce color PostScript dot plot from cpair */
1554
1555   FILE *wastl;
1556   int i, length;
1557
1558   length= strlen(seq);
1559   wastl = PS_dot_common(seq, wastlfile, NULL, winSize);
1560   if (wastl==NULL)
1561     return 0; /* return 0 for failure */
1562
1563   fprintf(wastl, "/hsb {\n"
1564           "dup 0.3 mul 1 exch sub sethsbcolor\n"
1565           "} bind def\n\n"
1566           "%%BEGIN DATA\n");
1567
1568   if(winSize > 0)
1569     fprintf(wastl, "\n%%draw the grid\ndrawgrid_turn\n\n");
1570   else
1571     fprintf(wastl, "\n%%draw the grid\ndrawgrid\n\n");
1572   fprintf(wastl,"%%start of base pair probability data\n");
1573
1574   /* print boxes */
1575    i=0;
1576    while (pi[i].j>0) {
1577      fprintf(wastl,"%1.2f %1.2f hsb %d %d %1.6f ubox\n",
1578              pi[i].hue, pi[i].sat, pi[i].i, pi[i].j, sqrt(pi[i].p));/*sqrt??*/
1579
1580      if (pi[i].mfe)
1581        fprintf(wastl,"%1.2f %1.2f hsb %d %d %1.4f lbox\n",
1582                pi[i].hue, pi[i].sat, pi[i].i, pi[i].j, pi[i].p);
1583      i++;
1584    }
1585
1586    fprintf(wastl,"showpage\n"
1587            "end\n"
1588            "%%%%EOF\n");
1589    fclose(wastl);
1590    return 1; /* success */
1591 }
1592
1593 int PS_dot_plot_turn(char *seq, struct plist *pl, char *wastlfile, int winSize) {
1594   /* produce color PostScript dot plot from cpair */
1595
1596   FILE *wastl;
1597   int i, length;
1598
1599   length= strlen(seq);
1600   wastl = PS_dot_common(seq, wastlfile, NULL, winSize);
1601   if (wastl==NULL)
1602     return 0; /* return 0 for failure */
1603
1604   if(winSize > 0)
1605     fprintf(wastl, "\n%%draw the grid\ndrawgrid_turn\n\n");
1606   else
1607     fprintf(wastl, "\n%%draw the grid\ndrawgrid\n\n");
1608   fprintf(wastl,"%%start of base pair probability data\n");
1609   /* print boxes */
1610   i=0;
1611   while (pl[i].j>0) {
1612     fprintf(wastl,"%d %d %1.4f ubox\n",
1613             pl[i].i, pl[i].j, sqrt(pl[i].p));
1614     i++;
1615   }
1616
1617   fprintf(wastl,"showpage\n"
1618                 "end\n"
1619                 "%%%%EOF\n");
1620   fclose(wastl);
1621   return 1; /* success */
1622 }
1623
1624 static FILE * PS_dot_common(char *seq, char *wastlfile,
1625                             char *comment, int winsize) {
1626   /* write PS header etc for all dot plot variants */
1627   FILE *wastl;
1628   char name[31], *c;
1629   int i, length;
1630
1631   length= strlen(seq);
1632   wastl = fopen(wastlfile,"w");
1633   if (wastl==NULL) {
1634     fprintf(stderr, "can't open %s for dot plot\n", wastlfile);
1635     return NULL; /* return 0 for failure */
1636   }
1637   strncpy(name, wastlfile, 30);
1638   if ((c=strrchr(name, '_'))!=0) *c='\0';
1639
1640   fprintf(wastl,
1641           "%%!PS-Adobe-3.0 EPSF-3.0\n"
1642           "%%%%Title: RNA Dot Plot\n"
1643           "%%%%Creator: %s, ViennaRNA-%s\n"
1644           "%%%%CreationDate: %s", rcsid+5, VERSION, time_stamp());
1645   if (winsize>0)
1646     fprintf(wastl, "%%%%BoundingBox: 66 530 520 650\n");
1647   else
1648     fprintf(wastl, "%%%%BoundingBox: 66 211 518 662\n");
1649   fprintf(wastl,
1650           "%%%%DocumentFonts: Helvetica\n"
1651           "%%%%Pages: 1\n"
1652           "%%%%EndComments\n\n"
1653           "%%Options: %s\n", option_string());
1654
1655   if (comment) fprintf(wastl,"%% %s\n",comment);
1656
1657   fprintf(wastl,"%s", RNAdp_prolog);
1658
1659   fprintf(wastl,"DPdict begin\n"
1660           "%%delete next line to get rid of title\n"
1661           "270 665 moveto /Helvetica findfont 14 scalefont setfont "
1662           "(%s) show\n\n", name);
1663
1664   fprintf(wastl,"/sequence { (\\\n");
1665   for (i=0; i<strlen(seq); i+=255)
1666     fprintf(wastl, "%.255s\\\n", seq+i);
1667   fprintf(wastl,") } def\n");
1668   if (winsize>0)
1669     fprintf(wastl,"/winSize %d def\n",winsize);
1670   fprintf(wastl,"/len { sequence length } bind def\n\n");
1671   if (cut_point>0) fprintf(wastl,"/cutpoint %d def\n\n", cut_point);
1672
1673
1674   if (winsize>0)
1675   fprintf(wastl,"292 416 translate\n"
1676           "72 6 mul len 1 add winSize add 2 sqrt mul div dup scale\n");
1677   else
1678     fprintf(wastl,"72 216 translate\n"
1679           "72 6 mul len 1 add div dup scale\n");
1680   fprintf(wastl, "/Helvetica findfont 0.95 scalefont setfont\n\n");
1681
1682   if (winsize>0) {
1683     fprintf(wastl, "%s", RNAdp_prolog_turn);
1684     fprintf(wastl,"0.5 dup translate\n"
1685           "drawseq_turn\n"
1686           "45 rotate\n\n");
1687   }
1688   else
1689     fprintf(wastl,"drawseq\n"
1690             "0.5 dup translate\n"
1691             "%% draw diagonal\n"
1692             "0.04 setlinewidth\n"
1693             "0 len moveto len 0 lineto stroke \n\n");
1694   return(wastl);
1695 }
1696
1697 int PS_color_aln(const char *structure, const char *filename,
1698                  const char *seqs[], const char *names[]) {
1699   /* produce PS sequence alignment color-annotated by consensus structure */
1700
1701   int N,i,j,k,x,y,tmp,columnWidth;
1702   char *tmpBuffer,*ssEscaped,*ruler, *cons;
1703   char c;
1704   float fontWidth, fontHeight, imageHeight, imageWidth,tmpColumns;
1705   int length, maxName, maxNum, currPos;
1706   float lineStep,blockStep,consStep,ssStep,rulerStep,nameStep,numberStep;
1707   float maxConsBar,startY,namesX,seqsX, currY;
1708   float score,barHeight,xx,yy;
1709   int match,block;
1710   FILE *outfile;
1711   short *pair_table;
1712   char * colorMatrix[6][3] = {
1713     {"0.0 1", "0.0 0.6",  "0.0 0.2"},  /* red    */
1714     {"0.16 1","0.16 0.6", "0.16 0.2"}, /* ochre  */
1715     {"0.32 1","0.32 0.6", "0.32 0.2"}, /* turquoise */
1716     {"0.48 1","0.48 0.6", "0.48 0.2"}, /* green  */
1717     {"0.65 1","0.65 0.6", "0.65 0.2"}, /* blue   */
1718     {"0.81 1","0.81 0.6", "0.81 0.2"} /* violet */
1719   };
1720
1721   const char *alnPlotHeader =
1722         "%%!PS-Adobe-3.0 EPSF-3.0\n"
1723         "%%%%BoundingBox: %i %i %i %i\n"
1724         "%%%%EndComments\n"
1725         "%% draws Vienna RNA like colored boxes\n"
1726         "/box { %% x1 y1 x2 y2 hue saturation\n"
1727         "  gsave\n"
1728         "  dup 0.3 mul 1 exch sub sethsbcolor\n"
1729         "  exch 3 index sub exch 2 index sub rectfill\n"
1730         "  grestore\n"
1731         "} def\n"
1732         "%% draws a box in current color\n"
1733         "/box2 { %% x1 y1 x2 y2\n"
1734         "  exch 3 index sub exch 2 index sub rectfill\n"
1735         "} def\n"
1736         "/string { %% (Text) x y\n"
1737         " 6 add\n"
1738         " moveto\n"
1739         "  show\n"
1740         "} def\n"
1741         "0 %i translate\n"
1742         "1 -1 scale\n"
1743         "/Courier findfont\n"
1744         "[10 0 0 -10 0 0] makefont setfont\n";
1745
1746
1747   outfile = fopen(filename, "w");
1748
1749   if (outfile == NULL) {
1750     fprintf(stderr, "can't open file %s - not doing alignment plot\n",
1751             filename);
1752     return 0;
1753   }
1754
1755   columnWidth=60;            /* Display long alignments in blocks of this size */
1756   fontWidth=6;               /* Font metrics */
1757   fontHeight=6.5;
1758   lineStep=fontHeight+2;     /* distance between lines */
1759   blockStep=3.5*fontHeight;  /* distance between blocks */
1760   consStep=fontHeight*0.5;   /* distance between alignment and conservation curve */
1761   ssStep=2;                  /* distance between secondary structure line and sequences */
1762   rulerStep=2;               /* distance between sequences and ruler */
1763   nameStep=3*fontWidth;             /* distance between names and sequences */
1764   numberStep=fontWidth;      /* distance between sequeces and numbers */
1765   maxConsBar=2.5*fontHeight; /* Height of conservation curve */
1766   startY=2;                     /* "y origin" */
1767   namesX=fontWidth;             /* "x origin" */
1768
1769   /* Number of columns of the alignment */
1770   length=strlen(seqs[0]);
1771
1772   /* Allocate memory for various strings, length*2 is (more than)
1773          enough for all of them */
1774   tmpBuffer = (char *) space((unsigned) MAX2(length*2,columnWidth)+1);
1775   ssEscaped=(char *) space((unsigned) length*2);
1776   ruler=(char *) space((unsigned) length*2);
1777
1778   pair_table=make_pair_table(structure);
1779   /* Get length of longest name and count sequences in alignment*/
1780
1781   for (i=maxName=N=0; names[i] != NULL; i++) {
1782     N++;
1783     tmp=strlen(names[i]);
1784     if (tmp>maxName)  maxName=tmp;
1785   }
1786
1787
1788   /* x-coord. where sequences start */
1789   seqsX=namesX+maxName*fontWidth+nameStep;
1790
1791   /* calculate number of digits of the alignment length */
1792   snprintf(tmpBuffer,length, "%i",length);
1793   maxNum=strlen(tmpBuffer);
1794
1795
1796   /* Calculate bounding box */
1797   tmpColumns=columnWidth;
1798   if (length<columnWidth){
1799         tmpColumns=length;
1800   }
1801   imageWidth=ceil(namesX+(maxName+tmpColumns+maxNum)*fontWidth+2*nameStep+fontWidth+numberStep);
1802   imageHeight=startY+ceil((float)length/columnWidth)*((N+2)*lineStep+blockStep+consStep+ssStep+rulerStep);
1803
1804   /* Write postscript header including correct bounding box */
1805   fprintf(outfile,alnPlotHeader,0,0,(int)imageWidth,(int)imageHeight,(int)imageHeight);
1806
1807   /* Create ruler and secondary structure lines */
1808   i=0;
1809   /* Init all with dots */
1810   for (i=0;i<(length);i++){
1811         ruler[i]='.';
1812   }
1813   i=0;
1814   for (i=0;i<length;i++){
1815         /* Write number every 10th position, leave out block breaks */
1816         if ((i+1)%10==0 && (i+1)%columnWidth!=0){
1817           snprintf(tmpBuffer,length,"%i",i+1);
1818           strncpy(ruler+i,tmpBuffer,strlen(tmpBuffer));
1819         }
1820   }
1821   ruler[length]='\0';
1822
1823   /* Draw color annotation first */
1824   /* Repeat for all pairs */
1825   for (i=1; i<=length; i++) {
1826     if ((j=pair_table[i])>i) {
1827       /* Repeat for open and closing position */
1828       for (k=0;k<2;k++){
1829         int pairings, nonpair, s, col;
1830         int ptype[8] = {0,0,0,0,0,0,0,0};
1831         char *color;
1832         col = (k==0)?i-1:j-1;
1833         block=ceil((float)(col+1)/columnWidth);
1834         xx=seqsX+(col-(block-1)*columnWidth)*fontWidth;
1835         /* Repeat for each sequence */
1836         for (s=pairings=nonpair=0; s<N; s++) {
1837           ptype[BP_pair[ENCODE(seqs[s][i-1])][ENCODE(seqs[s][j-1])]]++;
1838         }
1839         for (pairings=0,s=1; s<=7; s++) {
1840           if (ptype[s]) pairings++;
1841         }
1842         nonpair=ptype[0];
1843         if (nonpair <=2) {
1844           color = colorMatrix[pairings-1][nonpair];
1845           for (s=0; s<N; s++) {
1846             yy=startY+(block-1)*(lineStep*(N+2)+blockStep+consStep+rulerStep)+ssStep*(block)+(s+1)*lineStep;
1847
1848             /* Color according due color information in pi-array, only if base pair is possible */
1849             if (BP_pair[ENCODE(seqs[s][i-1])][ENCODE(seqs[s][j-1])]) {
1850
1851               fprintf(outfile, "%.1f %.1f %.1f %.1f %s box\n",
1852                       xx,yy-1,xx+fontWidth,yy+fontHeight+1,color);
1853             }
1854           }
1855         }
1856       }
1857     }
1858   }
1859   free(pair_table);
1860
1861   /* Process rest of the output in blocks of columnWidth */
1862
1863   currY=startY;
1864   currPos=0;
1865
1866   cons =  consensus(seqs);
1867
1868   while (currPos<length) {
1869
1870     /* Display secondary structure line */
1871     fprintf(outfile,"0 setgray\n");
1872     strncpy(tmpBuffer,structure+currPos,columnWidth);
1873     tmpBuffer[columnWidth]='\0';
1874
1875     x=0;y=0;
1876     while ((c=tmpBuffer[x])){
1877       if (c=='.'){
1878         ssEscaped[y++]='.';
1879       } else {
1880         ssEscaped[y++]='\\';
1881         ssEscaped[y++]=c;
1882       }
1883       x++;
1884     }
1885     ssEscaped[y]='\0';
1886
1887     fprintf(outfile, "(%s) %.1f %.1f string\n", ssEscaped,seqsX,currY);
1888     currY+=ssStep+lineStep;
1889
1890     /* Display names, sequences and numbers */
1891
1892     for (i=0; i<N; i++) {
1893
1894       strncpy(tmpBuffer,seqs[i]+currPos,columnWidth);
1895       tmpBuffer[columnWidth]='\0';
1896
1897       match=0;
1898       for (j=0;j<(currPos+strlen(tmpBuffer));j++){
1899         if (seqs[i][j] != '-') match++;
1900       }
1901
1902       fprintf(outfile, "(%s) %.1f %.1f string\n", names[i],namesX,currY);
1903       fprintf(outfile, "(%s) %.1f %.1f string\n", tmpBuffer,seqsX,currY);
1904       fprintf(outfile, "(%i) %.1f %.1f string\n", match,seqsX+fontWidth*(strlen(tmpBuffer))+numberStep,currY);
1905       currY+=lineStep;
1906     }
1907     currY+=rulerStep;
1908     strncpy(tmpBuffer,ruler+currPos,columnWidth);
1909     tmpBuffer[columnWidth]='\0';
1910     fprintf(outfile, "(%s) %.1f %.1f string\n", tmpBuffer,seqsX,currY);
1911
1912     currY+=lineStep;
1913     currY+=consStep;
1914
1915     /*Display conservation bar*/
1916
1917     fprintf(outfile,"0.6 setgray\n");
1918     for (i=currPos;(i<currPos+columnWidth && i<length);i++){
1919       match=0;
1920       for (j=0;j<N;j++){
1921         if (cons[i] == seqs[j][i]) match++;
1922         if (cons[i]=='U' && seqs[j][i]=='T') match++;
1923         if (cons[i]=='T' && seqs[j][i]=='U') match++;
1924       }
1925       score=(float)(match-1)/(N-1);
1926
1927       if (cons[i] == '-' ||
1928           cons[i] == '_' ||
1929           cons[i] == '.'){
1930         score=0;
1931       }
1932
1933       barHeight=maxConsBar*score;
1934       if (barHeight==0){
1935         barHeight=1;
1936       }
1937
1938       xx=seqsX+(i-(columnWidth*currPos/columnWidth))*fontWidth;
1939
1940       fprintf(outfile,"%.1f %.1f %.1f %.1f box2\n",
1941               xx,
1942               currY+maxConsBar-barHeight,
1943               xx+fontWidth,
1944               currY+maxConsBar);
1945     }
1946
1947     currY+=blockStep;
1948     currPos+=columnWidth;
1949   }
1950   free(cons);
1951
1952   fprintf(outfile,"showpage\n");
1953   fclose(outfile);
1954
1955   free(tmpBuffer);
1956   free(ssEscaped);free(ruler);
1957
1958   return 0;
1959
1960 }
1961
1962 int aliPS_color_aln(const char *structure, const char *filename,
1963                     const char *seqs[], const char *names[]) {
1964   /* produce PS sequence alignment color-annotated by consensus structure */
1965
1966   int N,i,j,k,x,y,tmp,columnWidth;
1967   char *tmpBuffer,*ssEscaped,*ruler, *cons;
1968   char c;
1969   float fontWidth, fontHeight, imageHeight, imageWidth,tmpColumns;
1970   int length, maxName, maxNum, currPos;
1971   float lineStep,blockStep,consStep,ssStep,rulerStep,nameStep,numberStep;
1972   float maxConsBar,startY,namesX,seqsX, currY;
1973   float score,barHeight,xx,yy;
1974   int match,block;
1975   FILE *outfile;
1976   short *pair_table;
1977   char * colorMatrix[6][3] = {
1978     {"0.0 1", "0.0 0.6",  "0.0 0.2"},  /* red    */
1979     {"0.16 1","0.16 0.6", "0.16 0.2"}, /* ochre  */
1980     {"0.32 1","0.32 0.6", "0.32 0.2"}, /* turquoise */
1981     {"0.48 1","0.48 0.6", "0.48 0.2"}, /* green  */
1982     {"0.65 1","0.65 0.6", "0.65 0.2"}, /* blue   */
1983     {"0.81 1","0.81 0.6", "0.81 0.2"} /* violet */
1984   };
1985
1986   const char *alnPlotHeader =
1987         "%%!PS-Adobe-3.0 EPSF-3.0\n"
1988         "%%%%BoundingBox: %i %i %i %i\n"
1989         "%%%%EndComments\n"
1990         "%% draws Vienna RNA like colored boxes\n"
1991         "/box { %% x1 y1 x2 y2 hue saturation\n"
1992         "  gsave\n"
1993         "  dup 0.3 mul 1 exch sub sethsbcolor\n"
1994         "  exch 3 index sub exch 2 index sub rectfill\n"
1995         "  grestore\n"
1996         "} def\n"
1997         "%% draws a box in current color\n"
1998         "/box2 { %% x1 y1 x2 y2\n"
1999         "  exch 3 index sub exch 2 index sub rectfill\n"
2000         "} def\n"
2001         "/string { %% (Text) x y\n"
2002         " 6 add\n"
2003         " moveto\n"
2004         "  show\n"
2005         "} def\n"
2006         "0 %i translate\n"
2007         "1 -1 scale\n"
2008         "/Courier findfont\n"
2009         "[10 0 0 -10 0 0] makefont setfont\n";
2010
2011
2012   outfile = fopen(filename, "w");
2013   if (outfile == NULL) {
2014     fprintf(stderr, "can't open file %s - not doing alignment plot\n",
2015             filename);
2016     return 0;
2017   }
2018
2019   columnWidth=100;            /* Display long alignments in blocks of this size */
2020   fontWidth=6;               /* Font metrics */
2021   fontHeight=6.5;
2022   lineStep=fontHeight+2;     /* distance between lines */
2023   blockStep=3.5*fontHeight;  /* distance between blocks */
2024   consStep=fontHeight*0.5;   /* distance between alignment and conservation curve */
2025   ssStep=2;                  /* distance between secondary structure line and sequences */
2026   rulerStep=2;               /* distance between sequences and ruler */
2027   nameStep=3*fontWidth;             /* distance between names and sequences */
2028   numberStep=fontWidth;      /* distance between sequeces and numbers */
2029   maxConsBar=2.5*fontHeight; /* Height of conservation curve */
2030   startY=2;                     /* "y origin" */
2031   namesX=fontWidth;             /* "x origin" */
2032
2033   /* Number of columns of the alignment */
2034   length=strlen(seqs[0]);
2035
2036   /* Allocate memory for various strings, length*2 is (more than)
2037          enough for all of them */
2038   tmpBuffer = (char *) space((unsigned) columnWidth + length*2 );
2039   ssEscaped=(char *) space((unsigned) length*2 );
2040   ruler=(char *) space((unsigned) length*2  );
2041 /*   char * structur; */
2042 /*   structur = (char*) space((length+1)*sizeof(char)); */
2043 /*   structur = strdup(structure); */
2044 /*   for(i=0; i<length;i++){ */
2045 /*     if(structur[i] == '<') structur[i]='('; */
2046 /*     if(structur[i] == '>') structur[i]=')'; */
2047 /*   } */
2048 /*   structur[length]='\0';    */
2049 /*   printf("%s \n", structur); */
2050    pair_table=alimake_pair_table(structure);
2051   /* Get length of longest name and count sequences in alignment*/
2052
2053   for (i=maxName=N=0; names[i] != NULL; i++) {
2054     N++;
2055     tmp=strlen(names[i]);
2056     if (tmp>maxName)  maxName=tmp;
2057   }
2058
2059
2060   /* x-coord. where sequences start */
2061   seqsX=namesX+maxName*fontWidth+nameStep;
2062
2063   /* calculate number of digits of the alignment length */
2064   snprintf(tmpBuffer,length, "%i",length);
2065   maxNum=strlen(tmpBuffer);
2066
2067
2068   /* Calculate bounding box */
2069   tmpColumns=columnWidth;
2070   if (length<columnWidth){
2071         tmpColumns=length;
2072   }
2073   imageWidth=ceil(namesX+(maxName+tmpColumns+maxNum)*fontWidth+2*nameStep+fontWidth+numberStep);
2074   imageHeight=startY+ceil((float)length/columnWidth)*((N+2)*lineStep+blockStep+consStep+ssStep+rulerStep);
2075
2076   /* Write postscript header including correct bounding box */
2077   fprintf(outfile,alnPlotHeader,0,0,(int)imageWidth,(int)imageHeight,(int)imageHeight);
2078
2079   /* Create ruler and secondary structure lines */
2080   i=0;
2081   /* Init all with dots */
2082   for (i=0;i<(length);i++){
2083         ruler[i]='.';
2084   }
2085   i=0;
2086   for (i=0;i<length;i++){
2087         /* Write number every 10th position, leave out block breaks */
2088         if ((i+1)%10==0 && (i+1)%columnWidth!=0){
2089           snprintf(tmpBuffer,length,"%i",i+1);
2090           strncpy(ruler+i,tmpBuffer,strlen(tmpBuffer));
2091         }
2092   }
2093   ruler[length]='\0';
2094
2095   /* Draw color annotation first */
2096   /* Repeat for all pairs */
2097   for (i=1; i<=length; i++) {
2098     if ((j=pair_table[i])>i) {
2099       /* Repeat for open and closing position */
2100       for (k=0;k<2;k++){
2101         int pairings, nonpair, s, col;
2102         int ptype[8] = {0,0,0,0,0,0,0,0};
2103         char *color;
2104         col = (k==0)?i-1:j-1;
2105         block=ceil((float)(col+1)/columnWidth);
2106         xx=seqsX+(col-(block-1)*columnWidth)*fontWidth;
2107         /* Repeat for each sequence */
2108         for (s=pairings=nonpair=0; s<N; s++) {
2109           ptype[BP_pair[ENCODE(seqs[s][i-1])][ENCODE(seqs[s][j-1])]]++;
2110         }
2111         for (pairings=0,s=1; s<=7; s++) {
2112           if (ptype[s]) pairings++;
2113         }
2114         nonpair=ptype[0];
2115         if (nonpair <=2) {
2116           color = colorMatrix[pairings-1][nonpair];
2117           for (s=0; s<N; s++) {
2118             yy=startY+(block-1)*(lineStep*(N+2)+blockStep+consStep+rulerStep)+ssStep*(block)+(s+1)*lineStep;
2119
2120             /* Color according due color information in pi-array, only if base pair is possible */
2121             if (BP_pair[ENCODE(seqs[s][i-1])][ENCODE(seqs[s][j-1])]) {
2122
2123               fprintf(outfile, "%.1f %.1f %.1f %.1f %s box\n",
2124                       xx,yy-1,xx+fontWidth,yy+fontHeight+1,color);
2125             }
2126           }
2127         }
2128       }
2129     }
2130   }
2131   free(pair_table);
2132
2133   /* Process rest of the output in blocks of columnWidth */
2134
2135   currY=startY;
2136   currPos=0;
2137
2138   cons =  consensus(seqs);
2139
2140   while (currPos<length) {
2141
2142     /* Display secondary structure line */
2143     fprintf(outfile,"0 setgray\n");
2144     strncpy(tmpBuffer,structure+currPos,columnWidth);
2145     tmpBuffer[columnWidth]='\0';
2146
2147     x=0;y=0;
2148     while ((c=tmpBuffer[x])){
2149       if (c=='.'){
2150         ssEscaped[y++]='.';
2151       } else {
2152         ssEscaped[y++]='\\';
2153         ssEscaped[y++]=c;
2154       }
2155       x++;
2156     }
2157     ssEscaped[y]='\0';
2158
2159     fprintf(outfile, "(%s) %.1f %.1f string\n", ssEscaped,seqsX,currY);
2160     currY+=ssStep+lineStep;
2161
2162     /* Display names, sequences and numbers */
2163
2164     for (i=0; i<N; i++) {
2165
2166       strncpy(tmpBuffer,seqs[i]+currPos,columnWidth);
2167       tmpBuffer[columnWidth]='\0';
2168
2169       match=0;
2170       for (j=0;j<(currPos+strlen(tmpBuffer));j++){
2171         if (seqs[i][j] != '-') match++;
2172       }
2173
2174       fprintf(outfile, "(%s) %.1f %.1f string\n", names[i],namesX,currY);
2175       fprintf(outfile, "(%s) %.1f %.1f string\n", tmpBuffer,seqsX,currY);
2176       fprintf(outfile, "(%i) %.1f %.1f string\n", match,seqsX+fontWidth*(strlen(tmpBuffer))+numberStep,currY);
2177       currY+=lineStep;
2178     }
2179     currY+=rulerStep;
2180     strncpy(tmpBuffer,ruler+currPos,columnWidth);
2181     tmpBuffer[columnWidth]='\0';
2182     fprintf(outfile, "(%s) %.1f %.1f string\n", tmpBuffer,seqsX,currY);
2183
2184     currY+=lineStep;
2185     currY+=consStep;
2186
2187     /*Display conservation bar*/
2188
2189     fprintf(outfile,"0.6 setgray\n");
2190     for (i=currPos;(i<currPos+columnWidth && i<length);i++){
2191       match=0;
2192       for (j=0;j<N;j++){
2193         if (cons[i] == seqs[j][i]) match++;
2194         if (cons[i]=='U' && seqs[j][i]=='T') match++;
2195         if (cons[i]=='T' && seqs[j][i]=='U') match++;
2196       }
2197       score=(float)(match-1)/(N-1);
2198
2199       if (cons[i] == '-' ||
2200           cons[i] == '_' ||
2201           cons[i] == '.'){
2202         score=0;
2203       }
2204
2205       barHeight=maxConsBar*score;
2206       if (barHeight==0){
2207         barHeight=1;
2208       }
2209
2210       xx=seqsX+(i-(columnWidth*currPos/columnWidth))*fontWidth;
2211
2212       fprintf(outfile,"%.1f %.1f %.1f %.1f box2\n",
2213               xx,
2214               currY+maxConsBar-barHeight,
2215               xx+fontWidth,
2216               currY+maxConsBar);
2217     }
2218
2219     currY+=blockStep;
2220     currPos+=columnWidth;
2221   }
2222   free(cons);
2223   fprintf(outfile,"showpage\n");
2224   fclose(outfile);
2225   free(tmpBuffer);
2226   free(ssEscaped);free(ruler);
2227
2228   return 0;
2229
2230 }
2231
2232 /*###########################################*/
2233 /*# deprecated functions below              #*/
2234 /*###########################################*/
2235
2236 int PS_dot_plot(char *string, char *wastlfile) {
2237   /* this is just a wrapper to call PS_dot_plot_list */
2238   int i, j, k, length, maxl, mf_num;
2239   struct plist *pl;
2240   struct plist *mf;
2241
2242   length = strlen(string);
2243   maxl = 2*length;
2244   pl = (struct plist *)space(maxl*sizeof(struct plist));
2245   k=0;
2246   /*make plist out of pr array*/
2247   for (i=1; i<length; i++)
2248     for (j=i+1; j<=length; j++) {
2249       if (pr[iindx[i]-j]<PMIN) continue;
2250       if (k>=maxl-1) {
2251         maxl *= 2;
2252         pl = (struct plist *)xrealloc(pl,maxl*sizeof(struct plist));
2253       }
2254       pl[k].i = i;
2255       pl[k].j = j;
2256       pl[k++].p = pr[iindx[i]-j];
2257     }
2258   pl[k].i=0;
2259   pl[k].j=0;
2260   pl[k++].p=0.;
2261   /*make plist out of base_pair array*/
2262   mf_num = base_pair ? base_pair[0].i : 0;
2263   mf = (struct plist *)space((mf_num+1)*sizeof(struct plist));
2264   for (k=0; k<mf_num; k++) {
2265     mf[k].i = base_pair[k+1].i;
2266     mf[k].j = base_pair[k+1].j;
2267     mf[k].p = 0.95*0.95;
2268   }
2269   mf[k].i=0;
2270   mf[k].j=0;
2271   mf[k].p=0.;
2272   i = PS_dot_plot_list(string, wastlfile, pl, mf, "");
2273   free(mf);
2274   free(pl);
2275   return (i);
2276 }
2277