/* Phase Function Approach (Version 0.8) */ /* --- a software for reconstruction of rate */ /* constant distribution in complex decay data */ /* */ /* Zhuang Lab, Harvard University, 2006 */ /* */ /* Tested on Turbo C 2.0 */ #include #include #include #include #include #include #include #include #define PI 3.14159265358979323846264338 #define LC 0.88137358701954302523260932 /*Lattice constant h=log(1+sqrt(2))*/ #define LOG10 1 #define LINEAR 0 #define MAXOCTAVE 2048 /*Partition number for fast Fourier transform*/ #define MAXDATASIZE 8000 /*Maximum bins in the t-domain*/ #define RESOLUTION 1 /*Time resolution limit of t-domain data*/ #define LOGKBINS 416 /*Maximum bins in the lnk-domain (activation energy-domain)*/ /* Kernel routines for the "Phase Function Approach" */ /* */ /* Reference: Y. Zhou & X. Zhuang, Robust Reconstruction */ /* of the Rate Constant Distribution using the Phase */ /* Function Method, Biophys. J. 91:4045-4053 (2006) */ void ClearMat(int n, float *a); /* Sets matrix elements to zero */ void Survival(int datasize, float binsize, float *data, float *bins); /* Computes survival function from decay data */ void mesh(int left, int top, int right, int bottom, int scale); /* Draws meshes for data output */ float Max(float a,float b); /* ... */ float Min(float a,float b); /* ... */ float phi(float s, float c); /* Some float point functions */ float rate(int i); /* ... */ float Abs(float x); /* ... */ float FJ(float x, float y); /* ... */ float* ObtainPhase(int datasize, float *data); /* Computes phase function by Fourier transform and extrapolation */ float* SmoothOut(float *rugged, int option); /* Smoothes out noise in the extrapolated phase function beyond the diffraction limit */ float* finephi(float *f, float *phase, float meanlife); /* Refines phase function by compensating remainder of finite-difference */ /* T. Ooura's fast Fourier transform routine */ /* */ /* Reference: T. Ooura, Ooura's mathematical */ /* software packages */ /* http://momonga.t.u-tokyo.ac.jp/~ooura/ */ void rdft(int n, float wr, float wi, float *a); void cdft(int n, float wr, float wi, float *a); /* Main body of the program */ void main(){ int i,j,k,datasize; float d[MAXDATASIZE],hist; float meanlife; float cphase[LOGKBINS],Cumul[LOGKBINS+1]; float *func,f,g,KK,kk; FILE *Data; FILE *Anal; char file1[15], file2[15]; int GraphDriver=VGA; int GraphMode=VGAHI; printf("\n\n\n\n"); printf("\t\t Phase Function Approach (Version 0.8)\n\n"); printf("\t\t --- a software for reconstruction of rate\n\t\t constant distribution in complex decay data\n\n\t\t Released 2006\n\t\t Zhuang Lab, Harvard University\n"); printf("________________________________________________________________________________"); printf("\n\n\t\t (Please enter in `FILENAME.EXT' format) \n\t\t Input data file name: "); scanf("%s",file1); printf("\t\t Output data file name: "); scanf("%s",file2); /* Read t-domain data from the input file */ Data=fopen(file1,"r"); datasize=0; meanlife=0; while (fscanf(Data,"%f\n",&hist)!=EOF){ d[datasize]=hist; datasize++;} for(i=datasize;iKK/2) cphase[k]=1; } setcolor(YELLOW); Cumul[0]=0;Cumul[1]=0; cphase[0]=0; cphase[1]=0; /* Writing k-domain information to the output file */ setcolor(CYAN); for(i=0;ii) f*=pow(Abs(1.0-(rate(i)+rate(i+1))/(rate(j)+rate(j+1))),(cphase[j]-cphase[j+1])); else f*=1; } f*=pow((1+sqrt(1-0.25*(1-rate(i)/rate(i+1))*(1-rate(i)/rate(i+1))))/(1-rate(i)/rate(i+1))/0.5,cphase[i+1]-cphase[i]); f*=(rate(i)+rate(i+1))*meanlife*FJ(cphase[i],cphase[i+1]); Cumul[i+1]=Min(1,Cumul[i]+f*(1-2*rate(i)/(rate(i)+rate(i+1)))); *(func+i+48)=Cumul[i+1]-Cumul[i]; } finephi(&Cumul[0], &cphase[0], meanlife); for(i=0;i<48;i++) *(func+i)=0; for(i=512-48;i<512;i++) *(func+i)=0; func=SmoothOut(&func[0],1); for(i=1;i<512;i++) *(func+i)=*(func+i-1)+*(func+i); setcolor(MAGENTA); k=1; do while (k0.01){ if (/*(Cumul[min(LOGKBINS,k+j)]-Cumul[max(0,k-j)]>=0.1) ||*/ ((Cumul[min(LOGKBINS,k+j)]-Cumul[max(0,k-j)]<=0.005))) g=(Cumul[min(LOGKBINS,k+j)]-Cumul[max(0,k-j)])*0.5/((float)j); else g=Max(0,*(func+48+k+j)-(*(func+48+k-j)))*0.5/((float)j); fprintf(Anal,"%f\t%f\t%f\t%f\t%f\n",(rate(k)+rate(k+1))/RESOLUTION,Cumul[k],0.037*sqrt(Cumul[k]),j*0.025,g); circle((log10(rate(k))+4.1)*(int)(390.0/4.0),(10.0-500.0*g+0.5)*(int)(460.0/11.0),0); line((log10(rate(k-j))+4.1)*(int)(390.0/4.0),(10.0-500.0*g+0.5)*(int)(460.0/11.0),(log10(rate(k+j))+4.1)*(int)(390.0/4.0),(10.0-500.0*g+0.5)*(int)(460.0/11.0)); j=60; } if (j==36) { if (Cumul[min(LOGKBINS,k+j)]-Cumul[max(0,k-j)]<=0.005) g=(Cumul[min(LOGKBINS,k+j)]-Cumul[max(0,k-j)])*0.5/((float)j); else g=Max(0,*(func+48+k+j)-(*(func+48+k-j)))*0.5/((float)j); fprintf(Anal,"%f\t%f\t%f\t%f\t%f\n",(rate(k)+rate(k+1))/RESOLUTION,Cumul[k],0.037*sqrt(Cumul[k]),j*0.025,g); circle((log10(rate(k))+4.1)*(int)(390.0/4.0),(10.0-500.0*g+0.5)*(int)(460.0/11.0),0); line((log10(rate(k-j))+4.1)*(int)(390.0/4.0),(10.0-500.0*g+0.5)*(int)(460.0/11.0),(log10(rate(k+j))+4.1)*(int)(390.0/4.0),(10.0-500.0*g+0.5)*(int)(460.0/11.0)); } } k++; } while (k0)?x:(-x);} float rate(int i){ return((0.5E-4)*exp(i*LC/36.0)); } float phi(float s, float c){ if (Abs(c)<1E-20) return 0.5; else return (atan(s/c)/PI+((s*c)<0)); } float Max(float a, float b){ return (a>b)?a:b;} float Min(float a, float b){ return (a> 2; m2 = m << 1; n2 = n - 2; k = 0; for (j = 0; j <= m2 - 4; j += 4) { if (j < k) { xr = a[j]; xi = a[j + 1]; a[j] = a[k]; a[j + 1] = a[k + 1]; a[k] = xr; a[k + 1] = xi; } else if (j > k) { j1 = n2 - j; k1 = n2 - k; xr = a[j1]; xi = a[j1 + 1]; a[j1] = a[k1]; a[j1 + 1] = a[k1 + 1]; a[k1] = xr; a[k1 + 1] = xi; } k1 = m2 + k; xr = a[j + 2]; xi = a[j + 3]; a[j + 2] = a[k1]; a[j + 3] = a[k1 + 1]; a[k1] = xr; a[k1 + 1] = xi; l = m; while (k >= l) { k -= l; l >>= 1; } k += l; } } void cdft(int n, float wr, float wi, float *a) { void bitrv2(int n, float *a); int i, j, k, l, m; float wkr, wki, wdr, wdi, ss, xr, xi; m = n; while (m > 4) { l = m >> 1; wkr = 1; wki = 0; wdr = 1 - 2 * wi * wi; wdi = 2 * wi * wr; ss = 2 * wdi; wr = wdr; wi = wdi; for (j = 0; j <= n - m; j += m) { i = j + l; xr = a[j] - a[i]; xi = a[j + 1] - a[i + 1]; a[j] += a[i]; a[j + 1] += a[i + 1]; a[i] = xr; a[i + 1] = xi; xr = a[j + 2] - a[i + 2]; xi = a[j + 3] - a[i + 3]; a[j + 2] += a[i + 2]; a[j + 3] += a[i + 3]; a[i + 2] = wdr * xr - wdi * xi; a[i + 3] = wdr * xi + wdi * xr; } for (k = 4; k <= l - 4; k += 4) { wkr -= ss * wdi; wki += ss * wdr; wdr -= ss * wki; wdi += ss * wkr; for (j = k; j <= n - m + k; j += m) { i = j + l; xr = a[j] - a[i]; xi = a[j + 1] - a[i + 1]; a[j] += a[i]; a[j + 1] += a[i + 1]; a[i] = wkr * xr - wki * xi; a[i + 1] = wkr * xi + wki * xr; xr = a[j + 2] - a[i + 2]; xi = a[j + 3] - a[i + 3]; a[j + 2] += a[i + 2]; a[j + 3] += a[i + 3]; a[i + 2] = wdr * xr - wdi * xi; a[i + 3] = wdr * xi + wdi * xr; } } m = l; } if (m > 2) { for (j = 0; j <= n - 4; j += 4) { xr = a[j] - a[j + 2]; xi = a[j + 1] - a[j + 3]; a[j] += a[j + 2]; a[j + 1] += a[j + 3]; a[j + 2] = xr; a[j + 3] = xi; } } if (n > 4) { bitrv2(n, a); } } void rdft(int n, float wr, float wi, float *a) { void cdft(int n, float wr, float wi, float *a); int j, k; float wkr, wki, wdr, wdi, ss, xr, xi, yr, yi; if (n > 4) { wkr = 0; wki = 0; wdr = wi * wi; wdi = wi * wr; ss = 4 * wdi; wr = 1 - 2 * wdr; wi = 2 * wdi; if (wi >= 0) { cdft(n, wr, wi, a); xi = a[0] - a[1]; a[0] += a[1]; a[1] = xi; } for (k = (n >> 1) - 4; k >= 4; k -= 4) { j = n - k; xr = a[k + 2] - a[j - 2]; xi = a[k + 3] + a[j - 1]; yr = wdr * xr - wdi * xi; yi = wdr * xi + wdi * xr; a[k + 2] -= yr; a[k + 3] -= yi; a[j - 2] += yr; a[j - 1] -= yi; wkr += ss * wdi; wki += ss * (0.5 - wdr); xr = a[k] - a[j]; xi = a[k + 1] + a[j + 1]; yr = wkr * xr - wki * xi; yi = wkr * xi + wki * xr; a[k] -= yr; a[k + 1] -= yi; a[j] += yr; a[j + 1] -= yi; wdr += ss * wki; wdi += ss * (0.5 - wkr); } j = n - 2; xr = a[2] - a[j]; xi = a[3] + a[j + 1]; yr = wdr * xr - wdi * xi; yi = wdr * xi + wdi * xr; a[2] -= yr; a[3] -= yi; a[j] += yr; a[j + 1] -= yi; if (wi < 0) { a[1] = 0.5 * (a[0] - a[1]); a[0] -= a[1]; cdft(n, wr, wi, a); } } else { if (wi < 0) { a[1] = 0.5 * (a[0] - a[1]); a[0] -= a[1]; } if (n > 2) { xr = a[0] - a[2]; xi = a[1] - a[3]; a[0] += a[2]; a[1] += a[3]; a[2] = xr; a[3] = xi; } if (wi >= 0) { xi = a[0] - a[1]; a[0] += a[1]; a[1] = xi; } } } void mesh(int left, int top, int right, int bottom, int scale) { int i,j,x=(right-left)/4,y=(bottom-top)/11.0; char str[3]; setcolor(CYAN); setviewport(left, top, right, bottom, 1); clearviewport(); rectangle(0,0,right-left,bottom-top); for(j=0;j<4;j++){ for(i=1;i<=10;i++){ if (scale) { setcolor(BLUE); line((j+log10((float)(i))+0.1)*x,1,(j+log10((float)(i))+0.1)*x,bottom-top-1); } else { setcolor(MAGENTA); line((j*10+i+0.1)*x/10,1,(j*10+i+0.1)*x/10,bottom-top-1); } } } setviewport(left, top, right, bottom, 1); for(j=0;j<=10;j++) line(1,((float)j+0.5)*y,right-left-1,((float)j+0.5)*y); } float* SmoothOut(float *rugged, int option){ float c,a,b,f; int i,j; rdft(512,cos(PI/512.0),sin(PI/512.0),&rugged[0]); if (option==1){ for(i=2;i<512;i++) (*(rugged+i))*=sin(PI*(float)(i/2)/14.22)/(PI*(float)(i/2)/14.22); } if (option==3){ for(i=2;i<512;i++) (*(rugged+i))*=sin(PI*(float)(i/2)/21.33)/(PI*(float)(i/2)/21.33); } rdft(512,cos(PI/512.0),-sin(PI/512.0),&rugged[0]); for(i=0;i<512;i++) *(rugged+i)/=256.0; return(&rugged[0]); } float* ObtainPhase(int datasize, float *d){ float p[MAXOCTAVE],a[3][512],y,C,S,theta=PI/(float)MAXOCTAVE,big,small,yStep=12*LC/18.0/PI; int i,j,k; ClearMat(512*3,a); setviewport(400,10,630,470,1); for(k=0;k<=1;k++){ y=(1-2*LC/PI-5*yStep)+k*yStep; for(i=12-12*k;i<512-12+12*k;i++){ ClearMat(MAXOCTAVE,&p[0]); Survival(datasize, 2*PI/cos(y*PI/2)/(float)(MAXOCTAVE)/rate(i-48), d, p); S = 1*(datasize-p[0])*MAXOCTAVE*pow(cos(y*PI/2.0),2)/(2.0*PI); C = -S*tan(y*PI/2.0); for(j=0;ji) f*=pow(Abs(1.0-(rate(i)+rate(i+1))/(rate(j)+rate(j+1))),(a[2][24+j]-a[2][24+j+1])); else f*=1; } f*=pow((1+sqrt(1-0.25*(1-rate(i)/rate(i+1))*(1-rate(i)/rate(i+1))))/(1-rate(i)/rate(i+1))/0.5,a[2][24+i+1]-a[2][24+i]); f*=(rate(i)+rate(i+1))*meanlife*FJ(a[2][i+24],a[2][i+25]); Cumul[i+1]=Min(1,Cumul[i]+f*(1-2*rate(i)/(rate(i)+rate(i+1)))); } for (k=0;k<2;k++){ y=1-4*LC/PI/3.0/12.0+k*2*LC/PI/3.0/12.0; for (i=0;i<512;i++){ C = 0 ; S = 0; for (j=0;ji) f*=pow(Abs(1.0-(rate(i)+rate(i+1))/(rate(j)+rate(j+1))),(a[2][24+j]-a[2][24+j+1])); else f*=1; } f*=pow((1+sqrt(1-0.25*(1-rate(i)/rate(i+1))*(1-rate(i)/rate(i+1))))/(1-rate(i)/rate(i+1))/0.5,a[2][24+i+1]-a[2][24+i]); f*=FJ(a[2][i+24],a[2][i+25]); Cumul[i+1]=Min(1,Cumul[i]+f*(1-2*rate(i)/(rate(i)+rate(i+1)))); line((log10(rate(i))+4.1)*(int)(390.0/4.0),(10.0-10.0*Cumul[i]+0.5)*(int)(460.0/11.0), (log10(rate(i+1))+4.1)*(int)(390.0/4.0),(10.0-10.0*Cumul[i+1]+0.5)*(int)(460.0/11.0)); } for(i=0;i