/* EGAD: moremath.h Navin Pokala and Tracy Handel Dept. of Molecular and Cell Biology University of California, Berkeley Copyright (C) 2003 Regents of the University of California GNU Public License Aug 12 2003 Absolutely no warranties are made or are implied with the use of this program or its parts. This file is the header for moremath.cpp. It also contains math-related macros and inlined functions for trigonometry, random-number generation, and 3D vector math. */ #ifndef moremath_header_flag #define moremath_header_flag #include #include #include #include #include #include #include #include "structure_types.h" #include "io.h" extern long IDUM; double average_array(int n, double *array); /* return the value of the greater number */ #define max(a,b) ((a)<(b) ? (b) : (a)) /* return the value of the lesser number */ #define min(a,b) ((a)>(b) ? (b) : (a)) /* This function assigns a 0 to each element of an array of type double The array must be ended w/ a ENDFLAG value */ void erase_double_array(double *array); #define cosd(n) (cos((n)*PI_over180)) /* cos of degree argument */ #define sind(n) (sin((n)*PI_over180)) /* sin of degree argument */ #define tand(n) (tan((n)*PI_over180)) /* tan of degree argument */ #define arccos(n) (PI_under180*(acos((n)))) /* arccos of double argument; returns degree value */ #define POW2(n) ((n)*(n)) /* returns n^2 */ /* set CARTESIAN a to (0,0,0) */ #define CARTESIAN_zero(a) { (a).x = 0; (a).y = 0; (a).z = 0; } /* return squared distance between a and b */ inline double Distance_sqrd(CARTESIAN a, CARTESIAN b) { double deltaX, deltaY, deltaZ; deltaX = (a.x-b.x); deltaY = (a.y-b.y); deltaZ = (a.z-b.z); return (POW2(deltaX) + POW2(deltaY) + POW2(deltaZ)); } /* return distance between a and b */ inline double Distance(const CARTESIAN& a, const CARTESIAN& b) { double deltaX, deltaY, deltaZ; deltaX = (a.x-b.x); deltaY = (a.y-b.y); deltaZ = (a.z-b.z); return sqrt(POW2(deltaX) + POW2(deltaY) + POW2(deltaZ)); } /* return length of vector a */ inline double vector_length(CARTESIAN a) { return(sqrt(POW2(a.x) + POW2(a.y) + POW2(a.z))); } /* return dotproduct for two vectors CARTESIAN a and CARTESIAN b */ inline double dotproduct(CARTESIAN a, CARTESIAN b) { return((a.x*b.x) + (a.y*b.y) + (a.z*b.z)); } /* multiply the dimensions of b with a; return the scaled vector */ CARTESIAN scalar_multiply_CARTESIAN(double a, CARTESIAN b); /* return CARTESIAN whose dimensions are the sum of the dimensions of a and b */ CARTESIAN add_CARTESIAN(CARTESIAN a, CARTESIAN b); /* return CARTESIAN whose dimensions are the difference of the dimensions of a and b */ CARTESIAN subtract_CARTESIAN(CARTESIAN a, CARTESIAN b); /* return the vector between a and b */ CARTESIAN vector_ab(CARTESIAN a, CARTESIAN b); /* return the midpoint between two CARTESIAN points a & b */ CARTESIAN midpoint(CARTESIAN a, CARTESIAN b); /* return unit vector for vector ab: a(x,y,z) --> b(x,y,z) */ inline CARTESIAN unitvector(CARTESIAN a, CARTESIAN b) { CARTESIAN q; double r; r = Distance(a,b); q.x = (a.x-b.x)/r; q.y = (a.y-b.y)/r; q.z = (a.z-b.z)/r; return(q); } /* return cross product between vectors a and b */ inline CARTESIAN cross_product(CARTESIAN a, CARTESIAN b) { CARTESIAN q; q.x = (a.y)*(b.z) - (a.z)*(b.y); q.y = (a.z)*(b.x) - (a.x)*(b.z); q.z = (a.x)*(b.y) - (a.y)*(b.x); return(q); } /* return unit crossproduct for two vectors a and b */ inline CARTESIAN unitcrossprod(CARTESIAN a, CARTESIAN b) { CARTESIAN q; double r; q = cross_product(a, b); r = vector_length(q); q.x=q.x/r; q.y=q.y/r; q.z=q.z/r; return(q); } /* return projection of vector u onto vector v */ inline CARTESIAN projection(CARTESIAN u, CARTESIAN v) { return( scalar_multiply_CARTESIAN( ( dotproduct(u,v)/pow(vector_length(v),2) ), v) ); } /* Given a transform *matrix (given as an array) and base coordinates *p, this function returns the transformed coordinates. p(x,y,z) -> (matrix) -> q(x',y',z') , where matrix is a 12 member array corresponding to: x' = matrix[1]*x + matrix[2]*y + matrix[3]*z + matrix[4] y' = matrix[5]*x + matrix[6]*y + matrix[7]*z + matrix[8] z' = matrix[9]*x + matrix[10]*y + matrix[11]*z + matrix[12] */ inline int check_transform_matrix(double *matrix) { /* check to make sure that this is a valid transform matrix */ if(fabs(matrix[1]*matrix[5] + matrix[2]*matrix[6] + matrix[3]*matrix[7]) > EPS || fabs(matrix[1]*matrix[9] + matrix[2]*matrix[10] + matrix[3]*matrix[11]) > EPS || fabs(matrix[5]*matrix[9] + matrix[6]*matrix[10] + matrix[7]*matrix[11]) > EPS || fabs(matrix[1]*matrix[2] + matrix[5]*matrix[6] + matrix[9]*matrix[10]) > EPS || fabs(matrix[2]*matrix[3] + matrix[6]*matrix[7] + matrix[10]*matrix[11]) > EPS || fabs(matrix[1]*matrix[3] + matrix[5]*matrix[7] + matrix[9]*matrix[11]) > EPS ) { fprintf(stderr,"ILLEGAL TRANSFORM MATRIX: non-orthogonal matrix\n"); fprintf(stderr,"| 1 5 + 2 6 + 3 7 | = %lf\n",fabs(matrix[1]*matrix[5] + matrix[2]*matrix[6] + matrix[3]*matrix[7])); fprintf(stderr,"| 1 9 + 2 10 + 3 11 | = %lf\n",fabs(matrix[1]*matrix[9] + matrix[2]*matrix[10] + matrix[3]*matrix[11])); fprintf(stderr,"| 5 9 + 6 10 + 7 11 | = %lf\n",fabs(matrix[5]*matrix[9] + matrix[6]*matrix[10] + matrix[7]*matrix[11])); fprintf(stderr,"| 1 2 + 5 6 + 9 10 | = %lf\n",fabs(matrix[1]*matrix[2] + matrix[5]*matrix[6] + matrix[9]*matrix[10])); fprintf(stderr,"| 2 3 + 6 7 + 10 11 | = %lf\n",fabs(matrix[2]*matrix[3] + matrix[6]*matrix[7] + matrix[10]*matrix[11])); fprintf(stderr,"| 1 3 + 5 7 + 9 11 | = %lf\n",fabs(matrix[1]*matrix[3] + matrix[5]*matrix[7] + matrix[9]*matrix[11])); fprintf(stderr,"x\t%lf\t%lf\t%lf\t|\t%lf\n",matrix[1],matrix[2],matrix[3],matrix[4]); fprintf(stderr,"y\t%lf\t%lf\t%lf\t|\t%lf\n",matrix[5],matrix[6],matrix[7],matrix[8]); fprintf(stderr,"z\t%lf\t%lf\t%lf\t|\t%lf\n",matrix[9],matrix[10],matrix[11],matrix[12]); exit(1); } /* check to make sure that this is a valid transform matrix */ if(fabs(matrix[1]*matrix[1] + matrix[2]*matrix[2] + matrix[3]*matrix[3] - 1.0) > EPS || fabs(matrix[5]*matrix[5] + matrix[6]*matrix[6] + matrix[7]*matrix[7] - 1.0) > EPS || fabs(matrix[9]*matrix[9] + matrix[10]*matrix[10] + matrix[11]*matrix[11] - 1.0) > EPS || fabs(matrix[1]*matrix[1] + matrix[5]*matrix[5] + matrix[9]*matrix[9] - 1.0) > EPS || fabs(matrix[2]*matrix[2] + matrix[6]*matrix[6] + matrix[10]*matrix[10] - 1.0) > EPS || fabs(matrix[3]*matrix[3] + matrix[7]*matrix[7] + matrix[11]*matrix[11] - 1.0) > EPS ) { fprintf(stderr,"ILLEGAL TRANSFORM MATRIX: un-normalized matrix\n"); fprintf(stderr,"| 1 1 + 2 2 + 3 3 | = %lf\n",fabs(matrix[1]*matrix[1] + matrix[2]*matrix[2] + matrix[3]*matrix[3])); fprintf(stderr,"| 5 5 + 6 6 + 7 7 | = %lf\n",fabs(matrix[5]*matrix[5] + matrix[6]*matrix[6] + matrix[7]*matrix[7])); fprintf(stderr,"| 9 9 + 10 10 + 11 11 | = %lf\n",fabs(matrix[9]*matrix[9] + matrix[10]*matrix[10] + matrix[11]*matrix[11])); fprintf(stderr,"| 1 1 + 5 5 + 9 9 | = %lf\n",fabs(matrix[1]*matrix[1] + matrix[5]*matrix[5] + matrix[9]*matrix[9])); fprintf(stderr,"| 2 2 + 6 6 + 10 10 | = %lf\n",fabs(matrix[2]*matrix[2] + matrix[6]*matrix[6] + matrix[10]*matrix[10])); fprintf(stderr,"| 3 3 + 7 7 + 11 11 | = %lf\n",fabs(matrix[3]*matrix[3] + matrix[7]*matrix[7] + matrix[11]*matrix[11])); fprintf(stderr,"x\t%lf\t%lf\t%lf\t|\t%lf\n",matrix[1],matrix[2],matrix[3],matrix[4]); fprintf(stderr,"y\t%lf\t%lf\t%lf\t|\t%lf\n",matrix[5],matrix[6],matrix[7],matrix[8]); fprintf(stderr,"z\t%lf\t%lf\t%lf\t|\t%lf\n",matrix[9],matrix[10],matrix[11],matrix[12]); exit(1); } return(1); } inline CARTESIAN tranform_coordinates(CARTESIAN p, double *matrix) { CARTESIAN q; if(matrix == NULL) { q = p; return(q); } /* perform the actual coordinate transformation */ q.x = matrix[1]*p.x + matrix[2]*p.y + matrix[3]*p.z + matrix[4]; q.y = matrix[5]*p.x + matrix[6]*p.y + matrix[7]*p.z + matrix[8]; q.z = matrix[9]*p.x + matrix[10]*p.y + matrix[11]*p.z + matrix[12]; return(q); } inline void translate_rotate_array_to_transform_matrix(double *translate_rotate_array, double *transform_matrix) { double t,s,c; double txy, txz, tyz, sx, sy, sz; CARTESIAN r_hat, zero_vector; r_hat.x = translate_rotate_array[4]; r_hat.y = translate_rotate_array[5]; r_hat.z = translate_rotate_array[6]; CARTESIAN_zero(zero_vector); r_hat = unitvector(r_hat, zero_vector); s = sind(translate_rotate_array[7]); c = cosd(translate_rotate_array[7]); t = 1.0 - c; txy = t*r_hat.x*r_hat.y; txz = t*r_hat.x*r_hat.z; tyz = t*r_hat.y*r_hat.z; sx = s*r_hat.x; sy = s*r_hat.y; sz = s*r_hat.z; transform_matrix[1] = t*r_hat.x*r_hat.x + c; transform_matrix[2] = txy + sz; transform_matrix[3] = txz - sy; transform_matrix[4] = translate_rotate_array[1]; // dx; transform_matrix[5] = txy - sz; transform_matrix[6] = t*r_hat.y*r_hat.y + c; transform_matrix[7] = tyz + sx; transform_matrix[8] = translate_rotate_array[2]; // dy; transform_matrix[9] = txz + sy; transform_matrix[10] = tyz - sx; transform_matrix[11] = t*r_hat.z*r_hat.z + c; transform_matrix[12] = translate_rotate_array[3]; // dz; } /* Random number generator ran2 from Numerical Recipies in C */ double rand2(); /* return a random double between 0 and 1 */ #define dice() ((rand2())) /* return a random double between 0 and n */ #define rand_double(n) (((double)(n))*dice()) /* return a random int between 0 and n, inclusive */ #define randint(n) ((int)rand_double((n)+1)) /* return random double between lower_bnd and upper_bnd */ #define rand_double_window(lower_bnd, upper_bnd) (rand_double(fabs((lower_bnd) - (upper_bnd))) + (lower_bnd)) /* Returns the linear correlation between two double arrays x and y of n elements Modified from the pearsn function from Numerical Recipies in C */ double correllation(double *x, double *y, int n); /* this function returns the r^2, rmsd, and rse between two arrays of n doubles */ STATISTICS statistics(double *actual, double *calc, int n); /* This function takes an array of probabalities It picks a random number between 0 and 1. If the random number is within the boundry for a given choice, it returns the index */ inline int roll_loaded_dice(double *probability, int array_size) { double doh; int n; doh = dice(); for(n=1;n<=array_size;++n) { if(probability[n]>doh) { if(probability[n-1]<=doh) { return(n); } } } return(array_size); } /* for two double arrays a and b, return the sum of a[i]*b[i], i=start->end */ double dot_prd(double *a, double *b, int start, int end); /* c[i] = a[i] + b[i], for i=start->end */ #define add_arrays(c, a, b, start, end) \ { \ int iiii; \ for(iiii=start;iiii<=end;++iiii) \ c[iiii] = a[iiii] + b[iiii]; \ } inline double array_sum(double *a, int start, int end) { int i; double sum; sum=0; for(i=start;i<=end;++i) sum += a[i]; return(sum); } /* c[i] = a[i] - b[i], for i=start->end */ #define array_diff(c, a, b, start, end) \ { \ int iiii; \ for(iiii=start;iiii<=end;++iiii) \ c[iiii] = a[iiii] - b[iiii]; \ } /* c[i] = 0, for i=start->end */ #define array_0(c, start, end) \ { \ int iiii; \ for(iiii=start;iiii<=end;++iiii) \ c[iiii] = 0; \ } /* for arrays x and y of ndata points, fit to a line; m = slope, b = y-intercept, sigma = fitting stat adapted from numerical recipies in C */ void linear_fit(double *x, double *y, int ndata, double *sigma, double *m, double *b); /* given a line y[i] = m*x[i] + b and and array of n x-values, calculate the y-values y must be allocated by the calling function */ void linear_transform(double *y, double *x, int n, double *m, double *b); /* returns determinant for an n x n matrix */ double laplace_determinant(int n, double **matrix); /* solves a system of n linear equations w/ n unknowns; form of AX=D where matrix_A contains coefficients; vector_X = variable vector; vector_D = constant term ie: 2x + 3y + 8z = 12 5x + 8y + 3z = 9 7x + 4y + z = 3 2 3 8 x 12 A = 5 8 3 X = y D = 9 7 4 1 z 3 */ void cramers_rule(int n, double **matrix_A, double *vector_X, double *vector_D); inline void falsify_bool_array(bool *bool_array, int start_index, int end_index) { int i; for(i=start_index;i<=end_index;++i) bool_array[i] = 0; } inline void copy_bool_array(bool *target_array, int target_start_index, bool *source_array, int source_start_index, int source_end_index) { int i,j; i=target_start_index; j=source_start_index; while(j<=source_end_index) { target_array[i] = source_array[j]; ++i; ++j; } } inline void union_bool_arrays(bool *target_array, int target_start_index, bool *source1_array, int source1_start_index, int source1_end_index, bool *source2_array, int source2_start_index, int source2_end_index) { int i,j,k; if(source1_end_index - source1_start_index != source2_end_index - source2_start_index) failure_report("ERROR in moremath.h: union_bool_arrays, source1 and source2 not same length","abort"); i=target_start_index; j=source1_start_index; k=source2_start_index; while(j<=source1_end_index) { target_array[i] = source1_array[j] || source2_array[k]; ++i; ++j; ++k; } } inline void intersection_bool_arrays(bool *target_array, int target_start_index, bool *source1_array, int source1_start_index, int source1_end_index, bool *source2_array, int source2_start_index, int source2_end_index) { int i,j,k; if(source1_end_index - source1_start_index != source2_end_index - source2_start_index) failure_report("ERROR in moremath.h: intersection_bool_arrays, source1 and source2 not same length","abort"); i=target_start_index; j=source1_start_index; k=source2_start_index; while(j<=source1_end_index) { target_array[i] = source1_array[j] && source2_array[k]; ++i; ++j; ++k; } } inline int any_false(bool *bool_array, int start_index, int stop_index) { int i; for(i=start_index;i<=stop_index;++i) if(bool_array[i] == 0) return(1); /* at least one false */ return(0); } inline int count_false(bool *bool_array, int start_index, int stop_index) { int i; int ctr; ctr=0; for(i=start_index;i<=stop_index;++i) if(bool_array[i] == 0) ++ctr; return(ctr); } inline int count_true(bool *bool_array, int start_index, int stop_index) { int i; int ctr; ctr=0; for(i=start_index;i<=stop_index;++i) if(bool_array[i] == 1) ++ctr; return(ctr); } void fprintf_bool_array(FILE *file_ptr, bool *bool_array, int start_index, int stop_index); void find_cube_corners_near_far_from_point(double min_x, double min_y, double min_z, double max_x, double max_y, double max_z, CARTESIAN *point, CARTESIAN *near_cube_corner, CARTESIAN *far_cube_corner); #endif