--- /dev/null
+/* Quadruple Statistics on an arbitrary alphabet
+
+ Vienna RNA Package --- Peter F Stadler 1993
+
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <strings.h>
+#include <string.h>
+#include <ctype.h>
+#include "utils.h"
+#include "PS3D.h"
+#include "distance_matrix.h"
+
+#define PUBLIC
+#define PRIVATE static
+
+#define MIN2(A, B) ((A) < (B) ? (A) : (B))
+#define MIN4(A, B, C, D) MIN2( (MIN2((A),(B))), (MIN2((C),(D))) )
+#define MAX2(A, B) ((A) > (B) ? (A) : (B))
+#define MAX4(A, B, C, D) MAX2( (MAX2((A),(B))), (MAX2((C),(D))) )
+
+PRIVATE int IBox[16];
+
+
+PUBLIC float *statgeom(char **seqs, int n_of_seqs);
+PUBLIC float *statgeom4(char **ss[4], int nn[4]);
+PUBLIC void printf_stg(float *B);
+
+PRIVATE void SingleBox(char *x1, char *x2, char *x3, char *x4);
+PRIVATE void SortSingleBox(void);
+
+/* ----------------------------------------------------------------------- */
+
+PUBLIC float *statgeom(char **seqs, int n_of_seqs)
+{
+ int i,j,k,l;
+ int i1;
+
+ float *B;
+ float temp;
+
+ if(n_of_seqs < 4) {
+ fprintf(stderr,"Less than 4 sequences for statistical geometry.\n");
+ return NULL;
+ }
+
+ B = (float *) space(16*sizeof(float));
+
+ for(i=3; i<n_of_seqs; i++) {
+ for(j=2; j<i; j++) {
+ for(k=1; k<j; k++) {
+ for(l=0; l<k; l++) {
+ SingleBox(seqs[i], seqs[j], seqs[k], seqs[l]);
+ SortSingleBox();
+ B[0] = (float) IBox[0]; /* transfer length */
+ for(i1=1;i1<=15;i1++) B[i1] += ( (float) IBox[i1] );
+ }
+ }
+ }
+ }
+
+ temp = (float) n_of_seqs;
+ temp = temp*(temp-1.)*(temp-2.)*(temp-3.)/24.;
+ for(i1=1;i1<=15;i1++) B[i1] /= (temp*B[0]);
+
+ return B;
+}
+
+/* ----------------------------------------------------------------------- */
+
+PUBLIC float *statgeom4(char **ss[4], int nn[4])
+{
+ int i,j,k,l;
+ int i1;
+
+ float *B;
+ float temp;
+
+ B = (float *) space(16*sizeof(float));
+
+ for(i=0; i<nn[0]; i++) {
+ for(j=0; j<nn[1]; j++) {
+ for(k=0; k<nn[2]; k++) {
+ for(l=0; l<nn[3]; l++) {
+ SingleBox(ss[0][i], ss[1][j], ss[2][k], ss[3][l]);
+ B[0] = (float) IBox[0]; /* transfer length */
+ for(i1=1;i1<=15;i1++) B[i1] += ( (float) IBox[i1] );
+ }
+ }
+ }
+ }
+
+ for (temp = 1, i=0;i<4;i++) temp *= ((float) nn[i]) ;
+ for(i1=1;i1<=15;i1++) B[i1] /= (temp*B[0]);
+
+ return B;
+}
+/* ------------------------------------------------------------------------- */
+
+PUBLIC void printf_stg(float *B)
+{
+ printf("> Statistical Geometry.\n");
+ printf("> %d (sequence length)\n", (int) B[0]);
+ printf("> AAAA\n");
+ printf(" %7.5f\n", B[1]);
+ printf("> BAAA ABAA AABA AAAB\n");
+ printf(" %7.5f %7.5f %7.5f %7.5f\n", B[2],B[3],B[4],B[5]);
+ printf("> AABB ABAB ABBA\n");
+ printf(" %7.5f %7.5f %7.5f\n", B[6],B[7],B[8]);
+ printf("> AABC ABAC ABCA BAAC BACA BCAA\n");
+ printf(" %7.5f %7.5f %7.5f %7.5f %7.5f %7.5f\n"
+ ,B[9],B[10],B[11],B[12],B[13],B[14]);
+ printf("> ABCD\n");
+ printf(" %7.5f\n",B[15]);
+}
+
+/* ------------------------------------------------------------------------- */
+
+
+PUBLIC void SimplifiedBox(float *B, char *filename)
+{
+ char *tmp, temp[50];
+ float X,Y,Z;
+ float T,P;
+ float t1, t2, t3, t4;
+ float p1, p2, p3, p4;
+ float x1, x2, y1, y2,z1,z2;
+ FILE *fp;
+ float view[3] = {0.3, 1.5, 0.1};
+ float axis[3] = {1.0, 0.0, 0.0};
+
+ t1 = B[9]+B[10]+B[11]+B[15];
+ t2 = B[9]+B[12]+B[13]+B[15];
+ t3 = B[10]+B[12]+B[14]+B[15];
+ t4 = B[11]+B[13]+B[14]+B[15];
+
+ p1 = B[2]; p2 = B[3]; p3 = B[4]; p4 = B[5];
+
+ x1 = B[6] + B[14];
+ x2 = B[6] + B[9];
+ y1 = B[7] + B[13];
+ y2 = B[7] + B[10];
+ z1 = B[8] + B[12];
+ z2 = B[8] + B[11];
+
+ X = (x1+x2)/2.; Y = (y1+y2)/2.; Z = (z1+z2)/2;
+ T = (t1+t2+t3+t4)/4.;
+ P = (p1+p2+p3+p4)/4.;
+
+ printf("> X = %7.5f \n",X);
+ printf("> Y = %7.5f \n",Y);
+ printf("> Z = %7.5f \n",Z);
+ printf("> T = %7.5f \n",T);
+ printf("> P = %7.5f \n",P);
+
+ tmp = get_taxon_label(-1); /* retrieve the dataset identifier */
+
+ temp[0]='\0';
+ if(tmp) { strcat(temp,tmp);strcat(temp,"_"); free(tmp); }
+ strcat(temp,filename);
+
+ fp = fopen(temp,"w");
+ if (fp!=NULL) {
+ ps3d_Preambel(fp, view, axis, "N");
+ PS_DrawSimplifiedBox(X,Y,Z,T,P, fp);
+ fclose(fp);
+ } else fprintf(stderr,"couldn't open %s -- not drawing box\n", temp);
+}
+
+
+
+/* ------------------------------------------------------------------------- */
+
+PRIVATE void SingleBox(char *x1, char *x2, char *x3, char *x4)
+{
+ int len,i1,j1,k1,i,M,m;
+ int d[4];
+ char t[4];
+
+ len=strlen(x1);
+ if(strlen(x2)!=len) nrerror("Sequences of unequal length in 'SingleBox'");
+ if(strlen(x3)!=len) nrerror("Sequences of unequal length in 'SingleBox'");
+ if(strlen(x4)!=len) nrerror("Sequences of unequal length in 'SingleBox'");
+
+ IBox[0] = len;
+ for(i=1; i<=15; i++) IBox[i] = 0;
+
+ for(i1=0;i1<len;i1++) {
+ t[0] = x1[i1];
+ t[1] = x2[i1];
+ t[2] = x3[i1];
+ t[3] = x4[i1];
+ for(j1=0;j1<4;j1++) {
+ d[j1]=4;
+ for(k1=0;k1<4;k1++) d[j1] -=(t[j1]!=t[k1]);
+ }
+ M=MAX4(d[0],d[1],d[2],d[3]);
+ switch(M) {
+ case 4 : /* Four of a kind */
+ IBox[1]++; /* A A A A */
+ break;
+ case 3 : /* Three of a kind */
+ if(d[0] == 1) IBox[2]++; /* B A A A */
+ else if(d[1] == 1) IBox[3]++; /* A B A A */
+ else if(d[2] == 1) IBox[4]++; /* A A B A */
+ else IBox[5]++; /* A A A B */
+ break;
+ case 2 :
+ m=MIN4(d[0],d[1],d[2],d[3]);
+ if(m==2){ /* Two Pairs */
+ if(t[1]==t[0]) IBox[6]++; /* A A B B */
+ else if(t[2]==t[0]) IBox[7]++; /* A B A B */
+ else IBox[8]++; /* A B B A */
+ }
+ else { /* One Pair */
+ if(d[0]==2){ /* 0 is in the pair */
+ if(t[1]==t[0]) IBox[9]++; /* A A B C */
+ else if(t[2]==t[0]) IBox[10]++; /* A B A C */
+ else IBox[11]++ ;/* A B C A */
+ }
+ else if(d[1]==2){ /* 1 is in the pair */
+ if(t[2]==t[1]) IBox[12]++; /* B A A C */
+ else IBox[13]++; /* B A C A */
+ }
+ else IBox[14]++; /* B C A A */
+ }
+ break;
+ case 1 : /* No Pair */
+ IBox[15]++; /* A B C D */
+ break;
+ default:
+ nrerror("This can't happen.");
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------- */
+
+PRIVATE void SortSingleBox(void)
+{
+ int i;
+ int M;
+ int IBB[16];
+ int s[4];
+
+ M = MAX2(MAX2( IBox[6], IBox[7] ), IBox[8] );
+
+ if( M== IBox[6] ) { /* 12|34 */
+ IBB[6] = IBox[6];
+ if( IBox[9] >= IBox[14] ) { /* 1,2 > 3,4 */
+ IBB[9] = IBox[9];
+ IBB[14] = IBox[14];
+ if( IBox[7] >= IBox[8] ) { /* 13|24 */
+ IBB[7] = IBox[7];
+ IBB[8] = IBox[8];
+ if( IBox[10] >= IBox[13] ) { /* 1>2>3>4 */
+ IBB[10] = IBox[10];
+ IBB[13] = IBox[13];
+ IBB[11] = IBox[11];
+ IBB[12] = IBox[12];
+ s[0]=0; s[1]=1; s[2]=2; s[3]=3; /* 1 2 3 4 */
+ }
+ else { /* 2>1>4>3 */
+ IBB[10] = IBox[13];
+ IBB[13] = IBox[10];
+ IBB[11] = IBox[12];
+ IBB[12] = IBox[11];
+ s[0]=1; s[1]=0; s[2]=3; s[3]=2; /* 2 1 4 3 */
+ }
+ }
+ else { /* 14|23 */
+ IBB[7] = IBox[8];
+ IBB[8] = IBox[7];
+ if( IBox[11] >= IBox[12]) { /* 1>2>4>3 */
+ IBB[10] = IBox[11];
+ IBB[13] = IBox[12];
+ IBB[11] = IBox[10];
+ IBB[12] = IBox[13];
+ s[0]=0; s[1]=1; s[2]=3; s[3]=2; /* 1 2 4 3 */
+ }
+ else { /* 2>1>3>4 */
+ IBB[10] = IBox[12];
+ IBB[13] = IBox[11];
+ IBB[11] = IBox[13];
+ IBB[12] = IBox[10];
+ s[0]=1; s[1]=0; s[2]=2; s[3]=3; /* 2 1 3 4 */
+ }
+ }
+ }
+ else { /* 3,4 > 1,2 */
+ IBB[9] = IBox[14];
+ IBB[14]= IBox[9];
+ if( IBox[7] >= IBox[8] ) { /* 31|42 */
+ IBB[7] = IBox[7];
+ IBB[8] = IBox[8];
+ if(IBox[10] >= IBox[13]) { /* 3>4>1>2 */
+ IBB[10] = IBox[10];
+ IBB[13] = IBox[13];
+ IBB[11] = IBox[12];
+ IBB[12] = IBox[11];
+ s[0]=2; s[1]=3; s[2]=0; s[3]=1; /* 3 4 1 2 */
+ }
+ else { /* */
+ IBB[10] = IBox[13];
+ IBB[13] = IBox[10];
+ IBB[11] = IBox[11];
+ IBB[12] = IBox[12];
+ s[0]=3; s[1]=2; s[2]=1; s[3]=0; /* 4 3 2 1 */
+ }
+ }
+ else { /* 32|14 */
+ IBB[7] = IBox[8];
+ IBB[8] = IBox[7];
+ if( IBox[11] >= IBox[12]) { /* 3>4>2>1 */
+ IBB[10] = IBox[12];
+ IBB[13] = IBox[11];
+ IBB[11] = IBox[10];
+ IBB[12] = IBox[13];
+ s[0]=2; s[1]=3; s[2]=1; s[3]=0; /* 3 4 2 1 */
+ }
+ else { /* 2>1>3>4 */
+ IBB[10] = IBox[10];
+ IBB[13] = IBox[13];
+ IBB[11] = IBox[12];
+ IBB[12] = IBox[11];
+ s[0]=2; s[1]=3; s[2]=0; s[3]=1; /* 3 4 1 2 */
+ }
+ }
+ }
+ }
+ else if (M == IBox[7] ) { /* 13|24 */
+ IBB[6] = IBox[7];
+ if( IBox[10] >= IBox[13] ) { /* 1,3 > 2,4 */
+ IBB[9] = IBox[10];
+ IBB[14] = IBox[13];
+ if( IBox[6] >= IBox[8]) { /* 12|34 */
+ IBB[7] = IBox[6];
+ IBB[8] = IBox[8];
+ if( IBox[9] >= IBox[14] ) { /* 1,2>3,4 */
+ IBB[10] = IBox[9];
+ IBB[13] = IBox[14];
+ IBB[11] = IBox[11];
+ IBB[12] = IBox[12];
+ s[0]=0; s[1]=2; s[2]=1; s[3]=3; /* 1 3 2 4 */
+ }
+ else {
+ IBB[10] = IBox[14];
+ IBB[13] = IBox[9];
+ IBB[11] = IBox[12];
+ IBB[12] = IBox[11];
+ s[0]=2; s[1]=0; s[2]=3; s[3]=1; /* 3 1 4 2 */
+ }
+ }
+ else { /* 14|23 */
+ IBB[7] = IBox[8];
+ IBB[8] = IBox[6];
+ if( IBox[11] >= IBox[12] ){ /* 1,4 > 2,3 */
+ IBB[10] = IBox[11];
+ IBB[13] = IBox[12];
+ IBB[11] = IBox[9];
+ IBB[12] = IBox[14];
+ s[0]=0; s[1]=2; s[2]=3; s[3]=1; /* 1 3 4 2 */
+ }
+ else {
+ IBB[10] = IBox[12];
+ IBB[13] = IBox[11];
+ IBB[11] = IBox[14];
+ IBB[12] = IBox[9];
+ s[0]=2; s[1]=0; s[2]=1; s[3]=3; /* 3 1 2 4 */
+ }
+ }
+ }
+ else { /* 2,4 > 1,3 */
+ IBB[9] = IBox[13];
+ IBB[14] = IBox[10];
+ if( IBox[6] >= IBox[8]) { /* 21|43 */
+ IBB[7] = IBox[6];
+ IBB[8] = IBox[8];
+ if( IBox[9] >= IBox[14] ) { /* 2,1>4,3 */
+ IBB[10] = IBox[9];
+ IBB[13] = IBox[14];
+ IBB[11] = IBox[12];
+ IBB[12] = IBox[11];
+ s[0]=1; s[1]=3; s[2]=0; s[3]=2; /* 2 4 1 3 */
+ }
+ else {
+ IBB[10] = IBox[14];
+ IBB[13] = IBox[9];
+ IBB[11] = IBox[11];
+ IBB[12] = IBox[12];
+ s[0]=3; s[1]=1; s[2]=2; s[3]=0; /* 4 2 3 1 */
+ }
+ }
+ else { /* 14|23 */
+ IBB[7] = IBox[8];
+ IBB[8] = IBox[6];
+ if( IBox[11] >= IBox[12] ){ /* 1,4 > 2,3 */
+ IBB[10] = IBox[11];
+ IBB[13] = IBox[12];
+ IBB[11] = IBox[14];
+ IBB[12] = IBox[9];
+ s[0]=3; s[1]=1; s[2]=0; s[3]=2; /* 4 2 1 3 */
+ }
+ else {
+ IBB[10] = IBox[12];
+ IBB[13] = IBox[11];
+ IBB[11] = IBox[9];
+ IBB[12] = IBox[14];
+ s[0]=1; s[1]=3; s[2]=2; s[3]=0; /* 2 4 3 1 */
+ }
+ }
+ }
+ }
+ else { /* 14 | 23 */
+ IBB[6] = IBox[8];
+ if( IBox[11] >= IBox[12] ) { /* 1,4 > 2,3 */
+ IBB[9] = IBox[11];
+ IBB[14] = IBox[12];
+ if(IBox[6] >= IBox[7]) { /* 12|34 */
+ IBB[7] = IBox[6];
+ IBB[8] = IBox[7];
+ if(IBox[9] >= IBox[14]) { /* 1,2>3,4 */
+ IBB[10] = IBox[9];
+ IBB[13] = IBox[14];
+ IBB[11] = IBox[10];
+ IBB[12] = IBox[13];
+ s[0]=0; s[1]=3; s[2]=1; s[3]=2; /* 1 4 2 3 */
+ }
+ else { /* 4,3>2,1 */
+ IBB[10] = IBox[14];
+ IBB[13] = IBox[9];
+ IBB[11] = IBox[13];
+ IBB[12] = IBox[10];
+ s[0]=3; s[1]=0; s[2]=2; s[3]=1; /* 4 1 3 2 */
+ }
+ }
+ else { /* 13|24 */
+ IBB[7] = IBox[7];
+ IBB[8] = IBox[6];
+ if( IBox[10] >= IBox[13]) { /* 1,3 > 2,4 */
+ IBB[10] = IBox[10];
+ IBB[13] = IBox[13];
+ IBB[11] = IBox[9];
+ IBB[12] = IBox[14];
+ s[0]=0; s[1]=3; s[2]=2; s[3]=1; /* 1 4 3 2 */
+ }
+ else {
+ IBB[10] = IBox[13];
+ IBB[13] = IBox[10];
+ IBB[11] = IBox[14];
+ IBB[12] = IBox[9];
+ s[0]=3; s[1]=0; s[2]=1; s[3]=2; /* 4 1 2 3 */
+ }
+ }
+ }
+ else { /* 2,3 > 1,4 */
+ IBB[9] = IBox[12];
+ IBB[14] = IBox[11];
+ if(IBox[6] >= IBox[7]) { /* 34|12 */
+ IBB[7] = IBox[6];
+ IBB[8] = IBox[7];
+ if(IBox[9] >= IBox[14]) { /* 2,1>4,3 */
+ IBB[10] = IBox[9];
+ IBB[13] = IBox[14];
+ IBB[11] = IBox[13];
+ IBB[12] = IBox[10];
+ s[0]=1; s[1]=2; s[2]=0; s[3]=3; /* 2 3 1 4 */
+ }
+ else { /* 4,3>2,1 */
+ IBB[10] = IBox[14];
+ IBB[13] = IBox[9];
+ IBB[11] = IBox[10];
+ IBB[12] = IBox[13];
+ s[0]=2; s[1]=1; s[2]=3; s[3]=0; /* 3 2 4 1 */
+ }
+ }
+ else { /* 31|42 */
+ IBB[7] = IBox[7];
+ IBB[8] = IBox[6];
+ if( IBox[10] >= IBox[13]) { /* 1,3 > 2,4 */
+ IBB[10] = IBox[10];
+ IBB[13] = IBox[13];
+ IBB[11] = IBox[14];
+ IBB[12] = IBox[9];
+ s[0]=2; s[1]=1; s[2]=0; s[3]=4; /* 3 2 1 4 */
+ }
+ else {
+ IBB[10] = IBox[13];
+ IBB[13] = IBox[10];
+ IBB[11] = IBox[9];
+ IBB[12] = IBox[14];
+ s[0]=1; s[1]=2; s[2]=3; s[3]=0; /* 2 3 4 1 */
+ }
+ }
+ }
+ }
+ /* HEUREKA */
+ IBB[0] = IBox[0];
+ IBB[1] = IBox[1];
+ IBB[2] = IBox[2+s[0]];
+ IBB[3] = IBox[2+s[1]];
+ IBB[4] = IBox[2+s[2]];
+ IBB[5] = IBox[2+s[3]];
+ /*
+ IBB[6...14] see above :-)
+ */
+ IBB[15] = IBox[15];
+
+ for(i=0;i<=15;i++) IBox[i] = IBB[i];
+}
+
+/* ----------------------------------------------------------------------- */
+
git://source.jalview.org / jabaws.git / blobdiff