--- /dev/null
+/* Ethan Wolf 1995 */
+
+#include <stdio.h>
+#include <math.h>
+#include "scio.h"
+#include "likelihood.h"
+#include "methods.h"
+#include "postscript.h"
+#include "scatter_ps.h"
+#include "compute_like.h"
+
+
+/* ax+b is the likelihood line approximation computed in likelihood region */
+/* from low % to hight % */
+void likelihood(long n[2], int mode, int filenum, double ratio, int minplace,
+ int maxplace, double *means[2], long *counts[2],
+ double gauss_mean[2], double sd[2], double down[2],
+ FILE *fout, FILE *foutps, FILE *scatter, FILE *scores,
+ double a[2], double b[2],
+ double low_percent, double high_percent, double error_allowed)
+{
+ int i, f;
+ double x, hist_height;
+ double _likelihood[2][201-2*MINSCORE];
+ double *likelihood[2];
+ double area1=0, area2=0;
+
+ int low[2]={0,0}; /* How many consectutive low % found. */
+ int high[2]={0,0}; /* How many consecutive high % found. */
+ int low_box[2], high_box[2]; /* These hold the first */
+ /* and last boxes that are low, high % likely, where these */
+ /* determined by finding the first box above the gaussian mean such */
+ /* that it has likelihood above low percent and the next two boxes */
+ /* are also above low-error percent. Similarly, high% is the first */
+ /* following box such that it has likelihood < %high and the two */
+ /* following have likelihood at least high-error */
+ double error=0;
+
+
+ likelihood[0] = &(_likelihood[0][-2*MINSCORE]);
+ likelihood[1] = &(_likelihood[1][-2*MINSCORE]);
+
+
+ if (mode & STOCK_METHOD) {
+ for(i=minplace; i<=maxplace; i++)
+ likelihood[0][i]= stock_like(ratio, (double)i/2, gauss_mean[0],sd[0],
+ down[0],gauss_mean[1], sd[1], down[1]);
+ filenum=1; /*** For printout purposes. ***/
+ }
+
+ else
+ for(f=0; f<filenum; f++) {
+
+ if (mode & RATIO_BOXES_GT_METHOD) {
+ area1= total_area_above_gauss(minplace, maxplace, counts[f],
+ gauss_mean[f],sd[f],down[f], n[f]);
+ area2= total_area(minplace,maxplace, counts[f], n[f]);
+ }
+ if (mode & FRACTION_POS_METHOD) {
+ area1=total_area_above_gauss(minplace, maxplace, counts[f],
+ gauss_mean[f],sd[f],down[f], n[f]);
+ area2=0;
+ }
+
+
+ for(i=minplace; i<=maxplace; i++) {
+ if (mode & MEAN_MODE) x= means[f][i];
+ else x= (double) i/2;
+
+/***** We need to convert counts of residue scores into heights that **/
+/***** are compatable with the gaussian (where area is assumed to be 1). **/
+/***** Hence, since each box is .5 wide, we take 2*count/n. **/
+ hist_height= 2*(double)counts[f][i]/n[f];
+
+ if (mode & RATIO_BOX_METHOD)
+ likelihood[f][i]= like_gauss_by_box(x, gauss_mean[f],sd[f],
+ down[f], hist_height);
+ if (mode & RATIO_TRAP_METHOD)
+ likelihood[f][i]= like_gauss_by_line(x, gauss_mean[f],sd[f],
+ down[f], hist_height);
+ if (mode & RATIO_BOXES_GT_METHOD)
+ likelihood[f][i]=ratio_PosToNeg_at_least_x(x, gauss_mean[f],
+ sd[f],down[f], &area1, &area2, hist_height);
+ if (mode & FRACTION_POS_METHOD)
+ likelihood[f][i]=fraction_of_positives_lower_than_x(x, gauss_mean[f],
+ sd[f], down[f], &area2,area1,hist_height);
+
+ if (low[f]<3) {
+ if ((likelihood[f][i] > low_percent - error) &&
+ ( (double)i/2 > gauss_mean[f]) ) { /* Guessing low percentage
+ is higher than mean. */
+ low[f]++;
+ error= error_allowed; /* For 2nd and 3d boxes.*/
+ if (low[f]==3) {
+ low_box[f]=i-2; /* First of 3 consecutive low-error boxes */
+ error=0; /* Now find first high percentage.*/
+ }
+ }
+ else {
+ low[f]=0;
+ error=0; /* Back to finding first low box. */
+ }
+ }
+ else if (high[f]<3)
+ if (likelihood[f][i] >= high_percent - error) {
+ high[f]++;
+ error=error_allowed;
+ if (high[f]==3)
+ high_box[f]=i-3; /* Last box < high_percent was three before.*/
+ }
+ else {
+ high[f]=0;
+ error=0; /* Back to finding first high box. */
+ }
+ }
+
+ if (high[f] < 3) /* didn't find a high_box, so make it last box. */
+ high_box[f]=maxplace;
+
+ if (low[f]< 3) /* didn't find a low_box so make it the mean. */
+ low_box[f] = (int)(2*gauss_mean[f]);
+
+
+
+ fprintf(fout, "For txt file %d:", f);
+ least_sq_fit(&a[f], &b[f],(double)low_box[f]/2,
+ (double)high_box[f]/2, means[f],likelihood[f],fout);
+ }
+
+ for(f=0;f<filenum; f++) {
+ fprintf(fout, "\n\nFor txt file %d:\n", f);
+ for(i=minplace; i<=maxplace; i++)
+ if (mode & MEAN_MODE)
+ fprintf(fout,"\nBox %lf (mean score %lf) has likelihood %lf",
+ (double)i/2 , means[f][i], likelihood[f][i]);
+ else
+ fprintf(fout,"\nBox %lf has likelihood %lf",
+ (double)i/2 , likelihood[f][i]);
+
+
+
+ /***Genetate postscript likelihood plot***/
+ post_script(foutps, minplace, maxplace, likelihood[f],
+ (mode & NO_HEADER_MODE), a[f], b[f],
+ (double)low_box[f]/2,(double)high_box[f]/2 );
+
+ /****Generate Scatter Plot****/
+ if (mode & SCATTER_MODE) {
+ scatter_header(scatter,a[f],b[f],10,90);
+ scatter_points(scatter,scores);
+ scatter_trailer(scatter);
+ }
+
+ }
+
+
+}
+
+
+/*** Likelihood line is l(x)=ax+b from low_score to high_score. ***/
+void post_script(FILE *foutps, int minplace, int maxplace, double *likelihood,
+ int no_header, double a, double b, double low_score,
+ double high_score)
+{
+ int i;
+
+ if (!no_header)
+ post_script_header(foutps);
+
+
+ for (i=minplace; i<=maxplace; i++)
+ fprintf(foutps,"%lf %lf bar\n",(double)i/2 , likelihood[i]);
+
+ if (!no_header)
+ post_script_like_line(foutps, a, b, low_score, high_score);
+
+ if (!no_header)
+ post_script_trailer(foutps);
+
+}
+
+
+
+/*** In doing the least square fit to line mx+b we get two equations: */
+/*** c1m + c2b = a1
+ c2m + c3b = a2
+ where c1 = # of points
+ c2= /sum x_value of points
+ c3= /sum x^2 for points
+ a1= /sum y_value of points
+ a2= /sum x*y for all the points
+*******/
+/******
+ Solving for m and b gives:
+ b= a2c2 -a1c3/denom and m= a1c2 - a2c1/denom
+ where denom = c2^2 - c1c3
+*********/
+
+void least_sq_fit(double *m, double *b, double leftend, double rightend,
+ double *means, double *likelihood, FILE *fout)
+{
+ int i;
+ double c1=0,c2=0,c3=0,a1=0,a2=0;
+
+ c1= 2*(rightend - leftend) +1;
+ for(i=2*leftend; i<=2*rightend; i++) {
+ c2 += means[i];
+ c3 += means[i]*means[i];
+ a1 += likelihood[i];
+ a2 += likelihood[i]*means[i];
+ }
+ *b= (a2*c2 - a1*c3)/ (c2*c2 - c1*c3);
+ *m= (a1*c2 - a2*c1)/ (c2*c2 - c1*c3);
+
+ fprintf(fout, "\nThe likelihood curve from %.2lf to %.2lf using means of the scores for x values is approximated by the line:\n", leftend, rightend);
+ fprintf(fout, " %.4lf * x + %.4lf\n\n",*m,*b);
+}
+
+
+
+/***Find an approximation to the score such that there are n1 residues with */
+/***lower scores and 2*n1/n=(1-percent). i.e. this mean gives a left side */
+/***which contains exactly half of the expected negatives (given by 1-%). */
+double percent_mean(long *counts, double *means, double percent,
+ int minplace,int maxplace, long number_residues)
+{
+ int box; /* This box will contain the score we want... */
+ int i;
+ long less_than; /** Number of scores less_than our current box. */
+ long greater_than;
+ double number_want_less_than, number_short, fraction_of_box;
+ double left_side, right_side;
+ double score_estimate;
+ double f; /* The fraction of the box we want. **/
+
+ number_want_less_than = (double)number_residues * (1-percent)/2;
+
+ i=minplace;
+ less_than=0;
+ while ( (double)less_than < number_want_less_than) {
+ less_than += counts[i];
+ i++;
+ }
+
+ box=i-1;
+ left_side= (double)box/2 - .25;
+ right_side= (double)box/2 +.25;
+ greater_than = number_residues - less_than;
+ less_than -= counts[box];
+
+ number_short= number_want_less_than - (double) less_than;
+/*** We want to find the "number_short"th score in box.***/
+ f= (double) number_short/counts[box];
+
+/*******OPTION 1*/
+ score_estimate= left_side +
+ f* (right_side -left_side);
+/*********/
+
+/********OPTION 2: quadratic through (0,left_side) (1,right) (1/2,mean)**
+ score_estimate = f*f*(2*(right_side+left_side) - 4*means[box])
+ + f*(4*means[box] - 3*left_side-right_side)
+ + left_side;
+***************/
+
+ return(score_estimate);
+}
+
+
+
+
+/***** half_box_mean compute the mean of the input file by computing **/
+/***** the mean of all points in boxes that have count at least half **/
+/***** of the max count **/
+
+
+double half_box_mean(int minplace, int maxplace, long *counts, double *means)
+{
+ int max_count=0, first_half_box, last_half_box, half_count=0, i;
+ double new_gauss_mean=0;
+
+ for (i=minplace; i<=maxplace; ++i)
+ if (counts[i]>max_count) max_count=counts[i];
+ for (first_half_box=minplace; counts[first_half_box] < (double)max_count/2;
+ ++first_half_box);
+ for (last_half_box=maxplace; counts[last_half_box] < (double)max_count/2;
+ --last_half_box);
+
+ for (i=first_half_box; i<=last_half_box; i++) {
+ new_gauss_mean += means[i]* (double)counts[i];
+ half_count += counts[i];
+ }
+ new_gauss_mean /= half_count;
+
+ fprintf(stderr, "\nGaussian mean based on half_max boxes = %lf\n",
+ new_gauss_mean);
+
+ return(new_gauss_mean);
+}
+
+
git://source.jalview.org / jpred.git / blobdiff