--- /dev/null
+#include <math.h>
+#include <stdio.h>
+#include "scconst.h"
+#include "scio.h"
+
+
+/** This computes 1-erf(x) with fractional error everywhere less than **/
+/** 1.2 * 10^-7. ***/
+/** I changed things from "float" to "double" from the numerical methods **/
+/** in C code. **/
+
+double erfcc(x)
+double x;
+{
+ double t,z,ans;
+
+ z=fabs(x);
+ t=1.0/(1.0+0.5*z);
+ ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
+ t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
+ t*(-0.82215223+t*0.17087277)))))))));
+ return x >= 0.0 ? ans : 2.0-ans;
+}
+
+
+
+
+double erf(double x)
+{
+ return(1-erfcc(x));
+}
+
+
+void printweight(long pfreqs[AANUM][POSNUM],
+ long pfreqp[AANUM][AANUM][POSNUM][POSNUM],
+ double gprobs[NUM_RES_TYPE],
+ double gprobp[NUM_RES_TYPE][NUM_RES_TYPE][POSNUM],
+ double weights[AANUM][POSNUM],
+ double weightp[AANUM][AANUM][POSNUM][POSNUM],
+ long totals[POSNUM],
+ long totalp[POSNUM][POSNUM], FILE *weight_file,
+ double pprobs[NUM_RES_TYPE][POSNUM],
+ double pprobp[NUM_RES_TYPE][NUM_RES_TYPE][POSNUM][POSNUM],
+ double bonnies_pprobs[NUM_RES_TYPE][POSNUM],
+ double bonnies_pprobp[NUM_RES_TYPE][NUM_RES_TYPE][POSNUM][POSNUM])
+{
+ int reg1, reg2, res1, res2;
+ double expect_prob, quantity,value;
+ long Num_samples;
+ int alert=0;
+
+ fprintf(weight_file,"\n\n SINGLE WEIGHTS \n\n");
+ fprintf(weight_file," ");
+ /* 9 spaces to get register to line up in columns*/
+ for(reg1=0; reg1<POSNUM; reg1++) {
+ fprintf(weight_file, "\t%c", 'a' + reg1);
+ }
+
+ fprintf(weight_file,"\n");
+
+ for(res1=0; res1<AANUM; res1++) {
+/* First print out the single weights. */
+ fprintf(weight_file,"\n\n%c, weight",numaa(res1));
+ for(reg1=0; reg1<POSNUM; reg1++)
+ fprintf(weight_file,"\t%.2lf", weights[res1][reg1]);
+
+/* Now print out the value that led to that weight. **/
+ fprintf(weight_file,"\n%c, values",numaa(res1));
+ for(reg1=0; reg1<POSNUM; reg1++) {
+ Num_samples = totals[reg1];
+ expect_prob = exp(gprobs[res1]);
+ quantity = ((double)pfreqs[res1][reg1]- (double)Num_samples*expect_prob);
+ value= fabs(quantity)/sqrt(Num_samples*expect_prob*(1-expect_prob));
+ fprintf(weight_file,"\t%.2lf", value);
+ }
+
+
+/* Now print out the number of samples. */
+ fprintf(weight_file,"\n%c, N ",numaa(res1));
+ for(reg1=0; reg1<POSNUM; reg1++) {
+ Num_samples = totals[reg1];
+ fprintf(weight_file,"\t%.2ld", Num_samples);
+ }
+/* Now print out the gprob values. */
+ fprintf(weight_file,"\n%c, gprobs",numaa(res1));
+ for(reg1=0; reg1<POSNUM; reg1++) {
+ expect_prob = exp(gprobs[res1]);
+ fprintf(weight_file,"\t%.2lf", expect_prob);
+ }
+/* Now print out the pprob values. */
+ fprintf(weight_file,"\n%c, pprobs",numaa(res1));
+ for(reg1=0; reg1<POSNUM; reg1++) {
+ fprintf(weight_file,"\t%.2lf", exp(pprobs[res1][reg1]));
+ }
+/* Now print out bonnie's pprob values. */
+ fprintf(weight_file,"\n%c, bprobs",numaa(res1));
+ for(reg1=0; reg1<POSNUM; reg1++) {
+ fprintf(weight_file,"\t%.2lf", exp(bonnies_pprobs[res1][reg1]));
+ }
+/* Now print out the num of pos samples values. */
+ fprintf(weight_file,"\n%c, NumPos",numaa(res1));
+ for(reg1=0; reg1<POSNUM; reg1++) {
+ fprintf(weight_file,"\t%.2ld", pfreqs[res1][reg1]);
+ }
+
+ }
+
+
+ fprintf(weight_file,"\n\n\nPAIR WEIGHTS");
+
+ for(reg1=0; reg1<POSNUM; reg1++) {
+ fprintf(weight_file,"\n\n");
+ for(reg2=0; reg2<POSNUM; reg2++)
+ fprintf(weight_file, "\t%c,%c", 'a' + reg1,'a' + reg2);
+ for(res1=0; res1<AANUM; res1++) {
+ for(res2=0; res2<AANUM; res2++) {
+ fprintf(weight_file,"\n\n%c,%c ",numaa(res1),numaa(res2));
+ for(reg2=0; reg2<POSNUM; reg2++)
+ fprintf(weight_file,"\t%.2lf",
+ weightp[res1][res2][reg1][reg2]);
+
+
+/* Now print out the pprob values. */
+ fprintf(weight_file,"\n%c,%c log_pprob",numaa(res1),numaa(res2));
+ for(reg2=0; reg2<POSNUM; reg2++) {
+ fprintf(weight_file,"\t%.2lf", pprobp[res1][res2][reg1][reg2]);
+ }
+/* Now print out bonnie's pprob values. */
+ fprintf(weight_file,"\n%c,%c log_bprob",numaa(res1),numaa(res2));
+ for(reg2=0; reg2<POSNUM; reg2++) {
+ fprintf(weight_file,"\t%.2lf",
+ bonnies_pprobp[res1][res2][reg1][reg2]);
+ }
+
+/* Now print out gprob values. */
+ fprintf(weight_file,"\n%c,%c log_gprob",numaa(res1),numaa(res2));
+ for(reg2=0; reg2<POSNUM; reg2++) {
+ fprintf(weight_file,"\t%.2lf",
+ gprobp[res1][res2][(reg2-reg1)%POSNUM]);
+ }
+
+/* Now print out the count on the amino acids. */
+ alert=0;
+ fprintf(weight_file,"\n%c,%c",numaa(res1),numaa(res2));
+ for(reg2=0; reg2<POSNUM; reg2++) {
+ fprintf(weight_file,"\t%.2ld",
+ pfreqp[res1][res2][reg1][reg2]);
+ if (pfreqp[res1][res2][reg1][reg2]<=1)
+ alert=1;
+ }
+ if (alert) fprintf(weight_file, "\nALERT: 0-1 frequency");
+
+
+ }
+ }
+ }
+}
+
+
+
+
+
+void printfreq_prob(long pfreqs[AANUM][POSNUM],
+ long pfreqp[AANUM][AANUM][POSNUM][POSNUM],
+ FILE *weight_file,
+ double pprobs[NUM_RES_TYPE][POSNUM],
+ double pprobp[NUM_RES_TYPE][NUM_RES_TYPE][POSNUM][POSNUM])
+{
+ int reg1, reg2, res1, res2;
+ double expect_prob, quantity,value;
+ long Num_samples;
+ int alert=0;
+
+ fprintf(weight_file,"\n\n SINGLE WEIGHTS \n\n");
+ fprintf(weight_file," ");
+ /* 9 spaces to get register to line up in columns*/
+ for(reg1=0; reg1<POSNUM; reg1++) {
+ fprintf(weight_file, "\t%c", 'a' + reg1);
+ }
+
+ fprintf(weight_file,"\n");
+
+ for(res1=0; res1<AANUM; res1++) {
+ /* Now print out the pprob values. */
+ fprintf(weight_file,"\n%c, pprobs",numaa(res1));
+ for(reg1=0; reg1<POSNUM; reg1++) {
+ fprintf(weight_file,"\t%.2lf", exp(pprobs[res1][reg1]));
+ }
+
+/* Now print out the num of pos samples values. */
+ fprintf(weight_file,"\n%c, NumPos",numaa(res1));
+ for(reg1=0; reg1<POSNUM; reg1++) {
+ fprintf(weight_file,"\t%.2ld", pfreqs[res1][reg1]);
+ }
+
+}
+
+
+ fprintf(weight_file,"\n\n\nPAIR WEIGHTS");
+
+ for(reg1=0; reg1<POSNUM; reg1++) {
+ fprintf(weight_file,"\n\n");
+ for(reg2=0; reg2<POSNUM; reg2++)
+ fprintf(weight_file, "\t%c,%c", 'a' + reg1,'a' + reg2);
+ for(res1=0; res1<AANUM; res1++) {
+ for(res2=0; res2<AANUM; res2++) {
+ fprintf(weight_file,"\n\n%c,%c ",numaa(res1),numaa(res2));
+
+/* Now print out the pprob values. */
+ fprintf(weight_file,"\n%c,%c log_pprob",numaa(res1),numaa(res2));
+ for(reg2=0; reg2<POSNUM; reg2++) {
+ fprintf(weight_file,"\t%.2lf", pprobp[res1][res2][reg1][reg2]);
+ }
+/* Now print out the count on the amino acids. */
+ alert=0;
+ fprintf(weight_file,"\n%c,%c",numaa(res1),numaa(res2));
+ for(reg2=0; reg2<POSNUM; reg2++) {
+ fprintf(weight_file,"\t%.2ld",
+ pfreqp[res1][res2][reg1][reg2]);
+ if (pfreqp[res1][res2][reg1][reg2]<=1)
+ alert=1;
+ }
+ if (alert) fprintf(weight_file, "\nALERT: 0-1 frequency");
+
+
+
+ }
+ }
+ }
+}
+
+
+/* Let x_i be 1 if kth sample at given positions give current amino acids.*/
+/* Then by central limit thm (sum(x_i) *Np)/sqrt(Np(1-p)) is normal. */
+int calcweight(long pfreqs[AANUM][POSNUM],
+ long pfreqp[AANUM][AANUM][POSNUM][POSNUM],
+ double gprobs[NUM_RES_TYPE],
+ double gprobp[NUM_RES_TYPE][NUM_RES_TYPE][POSNUM],
+ double weights[AANUM][POSNUM],
+ double weightp[AANUM][AANUM][POSNUM][POSNUM],
+ long totals[POSNUM],
+ long totalp[POSNUM][POSNUM], FILE *weight_file)
+{
+
+ int res1, res2, pos1, pos2;
+ long Num_samples;
+ double value;
+ double expect_prob;
+ double quantity;
+
+
+ for (pos1=0; pos1<POSNUM; pos1++) {
+ Num_samples = totals[pos1];
+ for (res1= 0; res1 < AANUM; res1++) {
+ expect_prob = exp(gprobs[res1]);
+ quantity = ((double)pfreqs[res1][pos1]- (double)Num_samples*expect_prob);
+ value= fabs(quantity)/sqrt(Num_samples*expect_prob*(1-expect_prob));
+ weights[res1][pos1] = erf(value/sqrt(2));
+ if (weights[res1][pos1]>1) weights[res1][pos1]=1;
+ }
+ }
+
+
+ for (pos1=0; pos1<POSNUM; pos1++)
+ for (pos2=0; pos2<POSNUM; pos2++) {
+ Num_samples = totalp[pos1][pos2];
+ for (res1= 0; res1 < AANUM; res1++)
+ for (res2= 0; res2 < AANUM; res2++) {
+ expect_prob = exp(gprobp[res1][res2][(res2 -res1)%POSNUM]);
+ quantity = ((double)pfreqp[res1][res2][pos1][pos2] -
+ (double)Num_samples*expect_prob);
+ value= fabs(quantity)/
+ sqrt(Num_samples*expect_prob*(1-expect_prob));
+ weightp[res1][res2][pos1][pos2] = erf(value/sqrt(2));
+ if (weightp[res1][res2][pos1][pos2] >1)
+ weightp[res1][res2][pos1][pos2]=1;
+
+ }
+ }
+
+/******* Now done in PairCoilData in sc2seq_interface.c
+ printweight(pfreqs,pfreqp,gprobs,gprobp,weights, weightp, totals,totalp,weight_file);
+*********/
+}
+
+
+
+/** Give the combination table a wt. that is the max of the wt. in the **/
+/** other two tables. i.e. if it is significant in one table, it should **/
+/** be significant in the combination??? Could just average??? **/
+
+
+void combine_many_weight(int current_table)
+{
+ int reg1, reg2, res1, res2;
+/*********
+ extern double many_weights[MAX_TABLE_NUMBER][AANUM][POSNUM];
+ extern double many_weightp[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM];
+**********/
+ double many_weights[MAX_TABLE_NUMBER][AANUM][POSNUM];
+ double many_weightp[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM];
+
+ for (res1=0; res1<AANUM; res1++)
+ for (reg1=0; reg1<POSNUM; reg1++) {
+ if (many_weights[0][res1][reg1] > many_weights[1][res1][reg1])
+ many_weights[current_table][res1][reg1] =
+ many_weights[0][res1][reg1];
+ else many_weights[current_table][res1][reg1] =
+ many_weights[1][res1][reg1];
+
+ for (res2=0; res2<AANUM; res2++)
+ for (reg2=0; reg2<POSNUM; reg2++)
+ if (many_weightp[0][res1][res2][reg1][reg2] >
+ many_weightp[1][res1][res2][reg1][reg2])
+ many_weightp[current_table][res1][res2][reg1][reg2] =
+ many_weightp[0][res1][res2][reg1][reg2];
+ else many_weightp[current_table][res1][res2][reg1][reg2] =
+ many_weightp[1][res1][res2][reg1][reg2];
+ }
+}
+
+
+
+/** Note that these weights should all be <= 1....Maybe even sum to 1 if ***/
+/** want to use them is a weighted avg. But in "geometric averaging" don't */
+/** have to sum to 1. Note that the bigger the sum though, the smaller **/
+/** the output prob. **/
+void calc_distance_weights(double dist_weight[POSNUM])
+{
+ int i;
+
+ dist_weight[0]= .3;
+ dist_weight[1]= .1;
+ dist_weight[2]= .3;
+ dist_weight[3]= .3;
+ dist_weight[4]= .2;
+ dist_weight[5]= .05;
+ dist_weight[6]= .1;
+
+ for(i=0; i<POSNUM; i++)
+ dist_weight[i]= 1;
+}
+
+
+
+double hard_cutoff(double x, double cutoff)
+{
+ double result;
+
+ if (x>cutoff) return(1);
+ else return(0);
+}
+
+/*** Want a sigma that maps 0 to 0 and 1 to 1. ***/
+/*** middle is the value that things should have things above it be mapped */
+/*** towards 1, and below it mapped towards 0. **/
+/*** Use 1+middle since want 0 mapped to 0 and middle to middle. **/
+/*** The larger power is, the more values are mapped towards the **/
+/*** extremes. **/
+double sigma(double x, int power, double middle)
+{
+ double result;
+
+ if (power ==-1) /** signal to do hardcutoff **/
+ return(hard_cutoff(x,middle));
+
+ if (power == 0) /** Just do nothing. **/
+ return(x);
+
+ if (x<=middle)
+ result= pow(1+middle,pow(x/middle,power)) - 1;
+ else
+ result= 2-pow(2-middle,pow( (1 - x)/(1-middle),power));
+
+ return(result);
+}
+
git://source.jalview.org / jpred.git / blobdiff