sources/multicoil/scatter_ps.c

   1 /*  Ethan Wolf  1995 */
   2
   3 #include <stdio.h>
   4 #include <math.h>
   5 #include "scatter_ps.h"
   6
   7 /*** To print a point do     ". NEWCOIL_SCORE PAIRCOIL_SCORE\n" **/
   8 /*** Then end list of all points with "showpage"**/
   9 /*** The parameters are:
  10   PairCoil likelihood linear region has equation like(x)=ax+b
  11          where the linear region runs from likelihood lower_like to upper_like.
  12
  13   For computing NewCoil likelihoods we use gaussian approximations for their
  14   scores on the negatives in genbnk with parameter mean-g and sd-g
  15   and for the positives with parameter mean-cc and sd-cc.
  16   Then their likelihood is   gauss-cc(x)/[gauss-cc(x) + 30*gauss-g(x)]
  17   (assuming that negatives are 30 times more prevalent than coils.
  18 ******/
  19
  20 void scatter_trailer(FILE *scatter)
  21 {
  22   fprintf(scatter,"showpage\n");
  23 }
  24
  25
  26
  27
  28 /** Read in scores output by "scboth -v"  which are in format:
  29       "seq# stock_score paircoil_score seq_code seq_name".
  30     Then output them as points in scatter plot as "stock_score paircoil_score ."
  31 **/
  32 void scatter_points(FILE *scatter, FILE *scores)
  33 {
  34   char line[300];
  35   double stock_score, paircoil_score;
  36
  37   while (fgets(line,300,scores)) {
  38     sscanf(line,"%*d %lf %lf", &stock_score, &paircoil_score);
  39     fprintf(scatter,"%lf %lf .\n", stock_score, paircoil_score);
  40   }
  41
  42 }
  43
  44
  45 void scatter_header(FILE *scatter, double a, double b, double lower_like, double upper_like)
  46 {
  47  fprintf(scatter,"%%!PS-Adobe-\n");
  48  fprintf(scatter,"%%%%Center and scale the figure\n");
  49  fprintf(scatter,"8.5 36 mul 11 36 mul translate\n");
  50  fprintf(scatter,"0.8 0.8 scale\n");
  51  fprintf(scatter,"-8.5 36 mul -11 36 mul translate\n");
  52
  53  fprintf(scatter,"\n%%%%The Figure Title.\n");
  54  fprintf(scatter,"/name () def\n");
  55  fprintf(scatter,"\n%%For printing on vertical and horizontal axis.\n");
  56  fprintf(scatter,"/rshow {dup stringwidth pop neg 0 rmoveto show} def\n");
  57  fprintf(scatter,"/cshow {dup stringwidth pop -2 div 0 rmoveto show} def\n");
  58  fprintf(scatter,"/label {/Times-Roman findfont 20 scalefont setfont cshow} def\n");
  59  fprintf(scatter,"/ttlabel {/Courier findfont 20 scalefont setfont cshow} def\n");
  60  fprintf(scatter,"%%%%Center and print the Title.\n");
  61  fprintf(scatter,"4.5 72 mul 9.5 72 mul moveto name label\n");
  62
  63  fprintf(scatter,"\n/adjust {ymin sub 10 mul exch xmin sub 10 mul exch} def\n");
  64
  65
  66  fprintf(scatter,"\n%%%%The following is used for printing out the scatter plot points.\n");
  67 /** Unnecessary Stuff??
  68  fprintf(scatter,"/choose { pop true } def   %% choose will be redefined for selective plotting\n");
  69  fprintf(scatter,"/. {choose {adjust .25 0 360 arc fill} {pop pop true} ifelse} def\n");
  70 **/
  71  fprintf(scatter,"/. {1 index 1 index ymin gt exch xmin gt and {adjust .9 0 360 arc fill} {pop pop}ifelse} def\n\n");
  72
  73
  74  fprintf(scatter,"/xmin -2 def\n");
  75  fprintf(scatter,"/ymin -26 def\n");
  76  fprintf(scatter,"/xmax 30 def\n");
  77  fprintf(scatter,"/ymax 28 def\n");
  78  fprintf(scatter,"/vline 1.35 ln 28 mul def\n");
  79
  80  fprintf(scatter,"/e 2.7182818284590452353 def\n\n");
  81
  82  fprintf(scatter,"%%%%For computing down*gaussian(x) with input x mean sd\n");
  83  fprintf(scatter,"/tgauss {  \n");
  84  fprintf(scatter,"        /stddev exch def sub stddev div dup mul -2 div  %% -[(x-mean)/sd]^2/2\n");
  85  fprintf(scatter,"        e exch exp                  %% e^<above>\n");
  86  fprintf(scatter,"        2.5066282746310002 div stddev div down mul\n");
  87  fprintf(scatter,"} def\n\n");
  88
  89
  90  stock_likelihood_setup(scatter);
  91
  92
  93  define_paircoil_likelihoodline(scatter, a, b, lower_like, upper_like);
  94
  95
  96
  97  fprintf(scatter,"\n/axes {             %% Do the axes\n");
  98  fprintf(scatter,"  90 100 translate\n");
  99  fprintf(scatter,"%%%% Horizontal\n");
 100  fprintf(scatter,"  0.5 setlinewidth\n");
 101  fprintf(scatter,"  xmin ymin adjust moveto xmax ymin adjust lineto stroke\n");
 102  fprintf(scatter,"    xmax ymin adjust moveto -8  4 rlineto stroke\n");
 103  fprintf(scatter,"    xmax ymin adjust moveto -8 -4 rlineto stroke\n");
 104  fprintf(scatter,"    vline ymin adjust moveto 0 -66 rmoveto\n");
 105  fprintf(scatter,"    (NEWCOILS) ttlabel\n");
 106  fprintf(scatter,"%%%% Vertical\n");
 107  fprintf(scatter,"  xmin ymin adjust moveto xmin ymax adjust lineto stroke\n");
 108  fprintf(scatter,"    xmin ymax adjust moveto  4 -8 rlineto stroke\n");
 109  fprintf(scatter,"    xmin ymax adjust moveto -4 -8 rlineto stroke\n");
 110  fprintf(scatter,"    xmin hline adjust moveto -80 0 rmoveto\n");
 111  fprintf(scatter,"    gsave 90 rotate (PairCoil) ttlabel grestore\n");
 112  fprintf(scatter,"%%%% Vertical Ticks\n");
 113  fprintf(scatter,"  -24 2 26 {xmin exch adjust moveto -3 0 rmoveto 6 0 rlineto stroke} for\n");
 114  fprintf(scatter,"%%%% Horizontal Ticks\n");
 115  fprintf(scatter,"  [10 11 12 13 14 15 16 18 20 22 24 26 28]\n");
 116  fprintf(scatter,"   {10 div ln 28 mul ymin adjust moveto 0 -3 rmoveto 0 6 rlineto stroke} forall\n");
 117  fprintf(scatter,"%%%% Horizontal Numbers\n");
 118  fprintf(scatter,"  /Times-Roman findfont 15 scalefont setfont\n");
 119  fprintf(scatter,"  [10 11 12 13 14 16 18 20 24 28]\n");
 120  fprintf(scatter,"   {10 div dup ln 28 mul ymin adjust -22 add moveto\n");
 121  fprintf(scatter,"                  5 string cvs cshow} forall\n");
 122
 123
 124  stock_percentages(scatter);
 125
 126
 127  fprintf(scatter,"%%%% Vertical Numbers (paircoil scores)\n");
 128  fprintf(scatter,"  -24 2 26 {dup xmin exch adjust moveto -10 -4 rmoveto\n");
 129  fprintf(scatter,"               4 string cvs rshow} for\n");
 130  fprintf(scatter,"%%%% Vertical Likelihoods (percentages)\n");
 131  fprintf(scatter,"  -10 2 18 {dup xmin exch adjust moveto -10 -4 rmoveto -60 0 rmoveto sfit show (%%) show} for\n");
 132  fprintf(scatter,"  newpath\n");
 133  fprintf(scatter,"} def          %%Ends the Definition of the Axes\n\n\n");
 134
 135
 136
 137  fprintf(scatter,"0 setgray\n");
 138  fprintf(scatter,"axes\n");
 139  fprintf(scatter,"newpath\n");
 140  fprintf(scatter,"vline dup\n");
 141  fprintf(scatter,"ymin adjust moveto ymax adjust lineto\n");
 142  fprintf(scatter,"stroke\n");
 143  fprintf(scatter,"hline dup\n");
 144  fprintf(scatter,"xmin exch adjust moveto xmax exch adjust lineto\n");
 145  fprintf(scatter,"stroke\n");
 146  fprintf(scatter,"newpath\n");
 147
 148 }
 149
 150
 151
 152
 153
 154 /***This prints out post script stuff for the least square line ***/
 155 /***between 20% and 80% determined by the rest of the likelihood program ***/
 156
 157 void define_paircoil_likelihoodline(FILE *scatter, double a, double b, double lower_like, double upper_like)
 158 {
 159   /** x sfit should ouput the likelihood score for x in the linear **/
 160   /** region of the likelihood graph (between lower_like and upper_like) **/
 161   /** where the curve is like=ax+b.                                      **/
 162
 163   fprintf(scatter,"\n%%%%x sfit should ouput the likelihood score for x in the linear region of the\n");
 164   fprintf(scatter,"%%%%likelihood graph (between lower_like and upper_like) where the curve is like=ax+b.\n");
 165
 166   fprintf(scatter,"/a %lf def\n", a);
 167   fprintf(scatter,"/b %lf def\n", b);
 168   fprintf(scatter,"%%%%hline is used to draw the 50 percent likelihood line for PairCoil\n");
 169   fprintf(scatter,"/hline .5 b sub a div def\n");
 170   fprintf(scatter,"/lower_like %.0lf def\n",lower_like);
 171   fprintf(scatter,"/upper_like %.0lf def\n",upper_like);
 172   fprintf(scatter,"/sfit {\n");
 173   fprintf(scatter,"  a mul b add 100 mul\n");
 174   fprintf(scatter,"  .5 add cvi               %%Rounding to nearest percentage.\n");
 175   fprintf(scatter,"  dup upper_like gt        %%If likelihood > upper_like\n");
 176   fprintf(scatter,"    {(>) show upper_like 3 string cvs}\n");
 177   fprintf(scatter,"                           %%then print > and make string upper_like.\n");
 178   fprintf(scatter,"    {dup lower_like lt     %%else if < lower_like \n");
 179   fprintf(scatter,"      {(<) show  lower_like 3 string cvs}\n");
 180   fprintf(scatter,"                           %%then print < and make string lower_like.\n");
 181   fprintf(scatter,"      {( ) show 3 string cvs} ifelse  %%else print a space and\n");
 182   fprintf(scatter,"                                      %%make likelihood a string and convert to text\n");
 183   fprintf(scatter,"    } ifelse               %%ends else for when upper_like>likelihood(x)\n");
 184   fprintf(scatter,"} def\n\n");
 185
 186 }
 187
 188
 189 void stock_likelihood_setup(FILE *scatter)
 190 {
 191  fprintf(scatter,"%%%%To compute the likelihood of STOCK use:\n");
 192  fprintf(scatter,"%%%%      gauss_positives(x)/[gauss_pos(x)+30*gauss_genbnk(x)]\n");
 193
 194  fprintf(scatter,"%%%%NewCoil parameters for computing likelihoods.\n");
 195  fprintf(scatter,"/mean-g 0.77 def\n");
 196  fprintf(scatter,"/sd-g 0.20 def\n");
 197  fprintf(scatter,"/mean-cc 1.63 def\n");
 198  fprintf(scatter,"/sd-cc 0.24 def\n");
 199  fprintf(scatter,"/down 1 def\n\n");
 200
 201  fprintf(scatter,"/prob { dup mean-g sd-g tgauss 30 mul exch\n");
 202  fprintf(scatter,"        mean-cc sd-cc tgauss dup 3 1 roll add div\n");
 203  fprintf(scatter,"} def\n");
 204  fprintf(scatter,"/prob-log {28 div 2.7182818284590452353 exch exp prob} def\n\n");
 205 }
 206
 207
 208 void stock_percentages(FILE *scatter)
 209 { fprintf(scatter,"%%%% Stock Percentages at 1.2, 1.3, 1.4.\n");
 210  fprintf(scatter,"  1.1 ln 28 mul ymin adjust -44 add moveto (< 1%) cshow\n");
 211  fprintf(scatter,"  1.6 ln 28 mul ymin adjust -44 add moveto (>99%) cshow\n");
 212  fprintf(scatter,"  [12 13 14]\n");
 213  fprintf(scatter,"   {10 div ln 28 mul dup ymin adjust -44 add moveto\n");
 214  fprintf(scatter,"    prob-log 100 mul cvi 12 string cvs cshow} forall\n\n");
 215 }