/* File: amdm.c Author: Andrew W. Moore Created: Thu Sep 15 21:01:13 EDT 1994 Description: Dynamically allocated and deallocated matrices Copyright (C) 1994, Andrew W. Moore */ #include #include #include "amdm.h" /* Dynamically allocated and deallocated matrices */ /*************** FIRST SOME STUFF COPIED OUT OF ***************/ /*************** NUMERICAL RECIPES IN `C'. ***************/ #define NEW_VERSION #ifdef NEW_VERSION #define MAX_SVD_ITERS 120 static float maxarg1,maxarg2; #define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\ (maxarg1) : (maxarg2)) #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) static double sqrarg; #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg) static int iminarg1,iminarg2; #define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\ (iminarg1) : (iminarg2)) float pythag(float a, float b) { float absa,absb; absa=fabs(a); absb=fabs(b); if (absa > absb) return absa*sqrt(1.0+SQR(absb/absa)); else return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+SQR(absa/absb))); } bool dym_copy_of_svdcmp(float **a, int m, int n, float w[], float **v) { float pythag(float a, float b); int flag,i,its,j,jj,k; int l = 0; /* To prevent "possible uninitialized variable" warning */ int nm = 0; /* To prevent "possible uninitialized variable" warning */ float anorm,c,f,g,h,s,scale,x,y,z,*rv1; rv1=am_malloc_floats(1 + n); /* change */ g=scale=anorm=0.0; for (i=1;i<=n;i++) { l=i+1; rv1[i]=scale*g; g=s=scale=0.0; if (i <= m) { for (k=i;k<=m;k++) scale += fabs(a[k][i]); if (scale) { for (k=i;k<=m;k++) { a[k][i] /= scale; s += a[k][i]*a[k][i]; } f=a[i][i]; g = -SIGN(sqrt(s),f); h=f*g-s; a[i][i]=f-g; for (j=l;j<=n;j++) { for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j]; f=s/h; for (k=i;k<=m;k++) a[k][j] += f*a[k][i]; } for (k=i;k<=m;k++) a[k][i] *= scale; } } w[i]=scale *g; g=s=scale=0.0; if (i <= m && i != n) { for (k=l;k<=n;k++) scale += fabs(a[i][k]); if (scale) { for (k=l;k<=n;k++) { a[i][k] /= scale; s += a[i][k]*a[i][k]; } f=a[i][l]; g = -SIGN(sqrt(s),f); h=f*g-s; a[i][l]=f-g; for (k=l;k<=n;k++) rv1[k]=a[i][k]/h; for (j=l;j<=m;j++) { for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k]; for (k=l;k<=n;k++) a[j][k] += s*rv1[k]; } for (k=l;k<=n;k++) a[i][k] *= scale; } } anorm=FMAX(anorm,(fabs(w[i])+fabs(rv1[i]))); } for (i=n;i>=1;i--) { if (i < n) { if (g) { for (j=l;j<=n;j++) v[j][i]=(a[i][j]/a[i][l])/g; for (j=l;j<=n;j++) { for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j]; for (k=l;k<=n;k++) v[k][j] += s*v[k][i]; } } for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0; } v[i][i]=1.0; g=rv1[i]; l=i; } for (i=IMIN(m,n);i>=1;i--) { l=i+1; g=w[i]; for (j=l;j<=n;j++) a[i][j]=0.0; if (g) { g=1.0/g; for (j=l;j<=n;j++) { for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j]; f=(s/a[i][i])*g; for (k=i;k<=m;k++) a[k][j] += f*a[k][i]; } for (j=i;j<=m;j++) a[j][i] *= g; } else for (j=i;j<=m;j++) a[j][i]=0.0; ++a[i][i]; } for (k=n;k>=1;k--) { for (its=1;its<=MAX_SVD_ITERS;its++) { /* change */ flag=1; for (l=k;l>=1;l--) { nm=l-1; if ((float)(fabs(rv1[l])+anorm) == anorm) { flag=0; break; } if ((float)(fabs(w[nm])+anorm) == anorm) break; } if (flag) { c=0.0; s=1.0; for (i=l;i<=k;i++) { f=s*rv1[i]; rv1[i]=c*rv1[i]; if ((float)(fabs(f)+anorm) == anorm) break; g=w[i]; h=pythag(f,g); w[i]=h; h=1.0/h; c=g*h; s = -f*h; for (j=1;j<=m;j++) { y=a[j][nm]; z=a[j][i]; a[j][nm]=y*c+z*s; a[j][i]=z*c-y*s; } } } z=w[k]; if (l == k) { if (z < 0.0) { w[k] = -z; for (j=1;j<=n;j++) v[j][k] = -v[j][k]; } break; } if (its == MAX_SVD_ITERS) /* change */ { am_free_floats(rv1,n+1); fprintf(stderr,"*** amdm.c, warning, SVD did not converge\n"); return(FALSE); } x=w[l]; nm=k-1; y=w[nm]; g=rv1[nm]; h=rv1[k]; f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y); g=pythag(f,1.0); f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x; c=s=1.0; for (j=l;j<=nm;j++) { i=j+1; g=rv1[i]; y=w[i]; h=s*g; g=c*g; z=pythag(f,h); rv1[j]=z; c=f/z; s=h/z; f=x*c+g*s; g = g*c-x*s; h=y*s; y *= c; for (jj=1;jj<=n;jj++) { x=v[jj][j]; z=v[jj][i]; v[jj][j]=x*c+z*s; v[jj][i]=z*c-x*s; } z=pythag(f,h); w[j]=z; if (z) { z=1.0/z; c=f*z; s=h*z; } f=c*g+s*y; x=c*y-s*g; for (jj=1;jj<=m;jj++) { y=a[jj][j]; z=a[jj][i]; a[jj][j]=y*c+z*s; a[jj][i]=z*c-y*s; } } rv1[l]=0.0; rv1[k]=f; w[k]=x; } } am_free_floats(rv1,n+1); return(TRUE); } #endif #ifdef OLD_VERSION static float at,bt,ct; #define PYTHAG(a,b) ((at=fabs(a)) > (bt=fabs(b)) ? \ (ct=bt/at,at*sqrt(1.0+ct*ct)) : (bt ? (ct=at/bt,bt*sqrt(1.0+ct*ct)): 0.0)) static float maxarg1,maxarg2; #define MAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\ (maxarg1) : (maxarg2)) #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) int dym_copy_of_svdcmp( float **a, int m, int n, float *w, float **v ) { int flag,i,its,j,jj,k,l,nm; float c,f,h,s,x,y,z; float anorm=0.0,g=0.0,scale=0.0; float *rv1 = am_malloc_floats(n+1); /* altered */ if ( Verbosity > 10000.0 ) { printf("Entered dym_copy_of_svdcmp().\n"); printf("m (number rows) = %d\n",m); printf("n (number cols) = %d\n",n); fprintf_2d_floats(stdout,"a",a,m+1,n+1,"\n"); wait_for_key(); } if (m < n) my_error("SVDCMP: You must augment A with extra zero rows"); for (i=1;i<=n;i++) { l=i+1; rv1[i]=scale*g; g=s=scale=0.0; if (i <= m) { for (k=i;k<=m;k++) scale += fabs(a[k][i]); if (scale) { for (k=i;k<=m;k++) { a[k][i] /= scale; s += a[k][i]*a[k][i]; } f=a[i][i]; g = -SIGN(sqrt(s),f); h=f*g-s; a[i][i]=f-g; if (i != n) { for (j=l;j<=n;j++) { for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j]; f=s/h; for (k=i;k<=m;k++) a[k][j] += f*a[k][i]; } } for (k=i;k<=m;k++) a[k][i] *= scale; } } w[i]=scale*g; g=s=scale=0.0; if (i <= m && i != n) { for (k=l;k<=n;k++) scale += fabs(a[i][k]); if (scale) { for (k=l;k<=n;k++) { a[i][k] /= scale; s += a[i][k]*a[i][k]; } f=a[i][l]; g = -SIGN(sqrt(s),f); h=f*g-s; a[i][l]=f-g; for (k=l;k<=n;k++) rv1[k]=a[i][k]/h; if (i != m) { for (j=l;j<=m;j++) { for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k]; for (k=l;k<=n;k++) a[j][k] += s*rv1[k]; } } for (k=l;k<=n;k++) a[i][k] *= scale; } } anorm=MAX(anorm,(fabs(w[i])+fabs(rv1[i]))); } for (i=n;i>=1;i--) { if (i < n) { if (g) { for (j=l;j<=n;j++) v[j][i]=(a[i][j]/a[i][l])/g; for (j=l;j<=n;j++) { for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j]; for (k=l;k<=n;k++) v[k][j] += s*v[k][i]; } } for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0; } v[i][i]=1.0; g=rv1[i]; l=i; } for (i=n;i>=1;i--) { l=i+1; g=w[i]; if (i < n) for (j=l;j<=n;j++) a[i][j]=0.0; if (g) { g=1.0/g; if (i != n) { for (j=l;j<=n;j++) { for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j]; f=(s/a[i][i])*g; for (k=i;k<=m;k++) a[k][j] += f*a[k][i]; } } for (j=i;j<=m;j++) a[j][i] *= g; } else { for (j=i;j<=m;j++) a[j][i]=0.0; } ++a[i][i]; } for (k=n;k>=1;k--) { for (its=1;its<=30;its++) { flag=1; for (l=k;l>=1;l--) { nm=l-1; if (fabs(rv1[l])+anorm == anorm) { flag=0; break; } if (fabs(w[nm])+anorm == anorm) break; } if (flag) { c=0.0; s=1.0; for (i=l;i<=k;i++) { f=s*rv1[i]; if (fabs(f)+anorm != anorm) { g=w[i]; h=PYTHAG(f,g); w[i]=h; h=1.0/h; c=g*h; s=(-f*h); for (j=1;j<=m;j++) { y=a[j][nm]; z=a[j][i]; a[j][nm]=y*c+z*s; a[j][i]=z*c-y*s; } } } } z=w[k]; if (l == k) { if (z < 0.0) { w[k] = -z; for (j=1;j<=n;j++) v[j][k]=(-v[j][k]); } break; } if (its == 30) { printf("WARNING: No convergence in 30 SVDCMP iterations. \nUsing mean output.\n"); am_free_floats(rv1,n+1); return 0; /* FAIL */ } x=w[l]; nm=k-1; y=w[nm]; g=rv1[nm]; h=rv1[k]; f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y); g=PYTHAG(f,1.0); f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x; c=s=1.0; for (j=l;j<=nm;j++) { i=j+1; g=rv1[i]; y=w[i]; h=s*g; g=c*g; z=PYTHAG(f,h); rv1[j]=z; c=f/z; s=h/z; f=x*c+g*s; g=g*c-x*s; h=y*s; y=y*c; for (jj=1;jj<=n;jj++) { x=v[jj][j]; z=v[jj][i]; v[jj][j]=x*c+z*s; v[jj][i]=z*c-x*s; } z=PYTHAG(f,h); w[j]=z; if (z) { z=1.0/z; c=f*z; s=h*z; } f=(c*g)+(s*y); x=(c*y)-(s*g); for (jj=1;jj<=m;jj++) { y=a[jj][j]; z=a[jj][i]; a[jj][j]=y*c+z*s; a[jj][i]=z*c-y*s; } } rv1[l]=0.0; rv1[k]=f; w[k]=x; } } am_free_floats(rv1,n+1); return 1; } #endif #undef SIGN #undef MAX #undef PYTHAG void dym_copy_of_svbksb( float **u, float w[], float **v, float b[], float x[], int m, int n ) { int jj,j,i; float s; float *tmp = am_malloc_floats(n+1); /* altered */ for (j=1;j<=n;j++) { s=0.0; if (w[j]) { for (i=1;i<=m;i++) s += u[i][j]*b[i]; s /= w[j]; } tmp[j]=s; } for (j=1;j<=n;j++) { s=0.0; for (jj=1;jj<=n;jj++) s += v[j][jj]*tmp[jj]; x[j]=s; } am_free_floats(tmp,n+1); } /*************** END OF STUFF COPIED OUT OF ***************/ /*************** NUMERICAL RECIPES IN `C'. ***************/ /**** Nicer interface to SVD ***/ bool float_is_nan(float x) /* Returns TRUE is x is NaN (which is different from Inf). This test is accomplished by seeing if x is a number for which <= 0 and >= 0 are both false. */ { bool non_neg = x >= 0.0; bool non_pos = x <= 0.0; return( !(non_neg || non_pos) ); } bool float_is_legal(float x) /* Returns TRUE if x represents a float number, and is not NaN or Inf or -Inf. */ { bool result = (x * 0.0) == 0.0; return(result); } bool tdarr_is_legal(float **tdarr,int rows,int cols) { bool result = TRUE; int i,j; for ( i = 0 ; result && i < rows ; i++ ) for ( j = 0 ; result && j < cols ; j++ ) result = (tdarr[i][j] * 0.0) == 0.0; return(result); } bool farr_is_legal(float *farr,int size) { bool result = TRUE; int i; for ( i = 0 ; result && i < size ; i++ ) result = (farr[i] * 0.0) == 0.0; return(result); } bool sing_val_decomp( float **tdarr, int rows, int cols, float **u, float *w, float **v ) /* Takes as input a matrix represented by a 2-d array "tdarr". The matrix has "rows" rows and "cols" columns, thus the legal entries in it are tdarr[i][j] where 0 <= i < rows and 0 <= j < cols. It computes the singular value decomposition and stores the result in u, v and w. (decribed below). PRECONDITIONS: tdarr is an array of vectors of floats. It contains "rows" vectors of floats. Each vector of floats contains "cols" elements. u is a similar 2-d array with exactly the same dimensions. v is a similar structure 2-d array, except its dimensions are: "cols" rows AND "cols" columns. w is a 1-d array (a vecoir) of floats with "cols" elements IMPORTANT: When this routine is called none of tdarr, u ,v or w must share any of the same memory. POST-CONDITIONS: Returns TRUE if succeeded in performing the computation. Returns FALSE, and sets everything to zero if failed. If successful: Let W = Diag(w) = matrix with W[ii] = w[i], W[ij] = 0 if i != j. Then on exit from this routine we will have tdarr = u * W * v^T where v^T is the transpose of v and * is matrix *'s furthermore v with be orthogonal and v will be orthogonal, meaning u^T u = Identity, thus, if square, U^-1 = U^T (ditto for V). This in turn means that if tdarr is a square matrix then tdarr^-1 = v * DW * u^T where DW[ii] = 1 / w[i] and DW[ij] = 0 if i != j See pages 59-70 of Numerical Recipes in C for further info. MODIFICATION: Feb 25th 1995, by Andrew Moore. When SVD gets called with an a matrix with all very small values (e.g. all O(1e-30)), it was returning some NaN values in the u and v vectors, as well as a Inf. This is due, I suspect, to magnitude errors in intermediate parts of the SVD calculation. The final result needn't have over or underflows. To fix this, the matrix passed to SVD is normalized so that its maximum magnitude component is 1. The resulting w vector is then unnormalized by the normalizing factor. This is okay, based on the fact (unproven but obviosish) that if the SVD of A is U Diag(W) V^T then the SVD of kA is U Diag(kW) V^T where k is a scalar. */ { int u_rows = rows; int u_cols = cols; int v_rows = cols; int v_cols = cols; int w_size = cols; bool result; float **nra = am_malloc_2d_floats(rows+1,cols+1); float **nrv = am_malloc_2d_floats(v_rows+1,v_cols+1); float *nrw = am_malloc_floats(w_size+1); int i,j; float max_mag_tdarr = 0.0; /* See Modification comments Feb 95 above */ char *problem; if ( cols > rows ) fprintf(stderr,"sing_val_decomp: warning: fewer rows than columns\n"); for ( i = 0 ; i < rows ; i++ ) for ( j = 0 ; j < cols ; j++ ) max_mag_tdarr = real_max(max_mag_tdarr,fabs(tdarr[i][j])); if ( max_mag_tdarr <= 0.0 ) max_mag_tdarr = 1.0; /* Hopeless case. The array is all zeros. Let's let numercial recipies SVD cope with this */ for ( i = 0 ; i < rows ; i++ ) for ( j = 0 ; j < cols ; j++ ) nra[i+1][j+1] = tdarr[i][j] / max_mag_tdarr; if ( Verbosity > 5000.0 ) { printf("SVD: max_mag_tdarr = %g\n",max_mag_tdarr); fprintf_2d_floats(stdout,"svd_tdarr",tdarr,rows,cols,"\n"); fprintf_2d_floats(stdout,"svd_nra",nra,rows+1,cols+1,"\n"); wait_for_key(); } if ( !dym_copy_of_svdcmp(nra,rows,cols,nrw,nrv) ) { /* If SVD fails */ if ( Verbosity > 60.0 ) { fprintf(stderr,"Singular Value Decomposition choked when it was"); fprintf(stderr,"asked to use the following matrix A:\n"); fprintf_2d_floats(stderr,"A",tdarr,rows,cols,"\n"); wait_for_key(); } set_2d_floats_constant(u,rows,cols,0.0); set_2d_floats_constant(v,v_rows,v_cols,0.0); set_floats_constant(w,w_size,0.0); result = FALSE; } else { for ( i = 0 ; i < u_rows ; i++ ) for ( j = 0 ; j < u_cols ; j++ ) u[i][j] = nra[i+1][j+1]; for ( i = 0 ; i < v_rows ; i++ ) for ( j = 0 ; j < v_cols ; j++ ) v[i][j] = nrv[i+1][j+1]; for ( i = 0 ; i < w_size ; i++ ) w[i] = nrw[i+1] * max_mag_tdarr; result = TRUE; } if ( Verbosity > 5000.0 ) { printf("SVD: result = %s\n",(result) ? "Okay" : "SVD Failed"); fprintf_2d_floats(stdout,"svd_nru",nra,rows+1,cols+1,"\n"); fprintf_2d_floats(stdout,"svd_u",u,rows,cols,"\n"); fprintf_2d_floats(stdout,"svd_nrv",nrv,v_rows+1,v_cols+1,"\n"); fprintf_2d_floats(stdout,"svd_v",v,v_rows,v_cols,"\n"); fprintf_floats(stdout,"svd_nrw",nrw,cols+1,"\n"); fprintf_floats(stdout,"svd_w",w,cols,"\n"); wait_for_key(); } problem = NULL; if ( !tdarr_is_legal(tdarr,rows,cols) ) problem = "bad-svd-input"; else if ( !tdarr_is_legal(u,rows,cols) || !tdarr_is_legal(v,rows,cols) || !farr_is_legal(w,w_size) ) problem = "bad-svd-output"; else if ( !result ) problem = "svd-failed"; if ( problem != NULL ) { int i,j; for ( i = 0 ; i < 6 ; i++ ) for ( j = 0 ; j <= i ; j++ ) fprintf(stderr,"*%s",(j==i)?"\n":""); fprintf(stderr,"******* Singular Value Decomp has a problem.\n"); for ( i = 0 ; i < 6 ; i++ ) for ( j = 0 ; j <= (6-i) ; j++ ) fprintf(stderr,"*%s",(j==(6-i))?"\n":""); fprintf(stderr,"\n\n"); if ( eq_string(problem,"bad-svd-input") ) { fprintf(stderr,"The problem is that SVD was passed a matrix with NaN\n"); fprintf(stderr,"and/or Inf entries (see tdarr parameter below)\n"); } else if ( eq_string(problem,"bad-svd-output") ) { fprintf(stderr,"The problem is that SVD was produced a matrix with NaN\n"); fprintf(stderr,"and/or Inf entries (see u,v,w parameter below)\n"); } else if ( eq_string(problem,"svd-failed") ) fprintf(stderr,"SVD failed to converge on anything.\n"); fprintf(stderr,"\n\n"); fprintf_2d_floats(stdout,"svd_tdarr",tdarr,rows,cols,"\n"); fprintf_2d_floats(stderr,"svd_u",u,rows,cols,"\n"); fprintf_2d_floats(stderr,"svd_v",v,v_rows,v_cols,"\n"); fprintf_floats(stderr,"svd_w",w,cols,"\n"); wait_for_key(); } free_2d_floats(nra,rows+1,cols+1); free_2d_floats(nrv,v_rows+1,v_cols+1); am_free_floats(nrw,w_size+1); return(result); } /****** End of nicer interface etc ***********/ #define DYM_CODE 4508 #define DYV_CODE 4509 void check_dym_code(dym *d,char *name) { if ( d == NULL ) { fprintf(stderr,"NULL dym passed in operation %s\n",name); my_error("dym data structure"); } if ( d->dym_code != DYM_CODE ) { fprintf(stderr,"Attempt to access a non-allocated DYnamic Matrix\n"); fprintf(stderr,"This is in the operation %s\n",name); my_error("dym data structure error"); } } void check_dyv_code(dyv *d,char *name) { if ( d == NULL ) { fprintf(stderr,"NULL dyv passed in operation %s\n",name); my_error("dyv data structure"); } if ( d->dyv_code != DYV_CODE ) { fprintf(stderr,"Attempt to access a non-allocated DYnamic Vector\n"); fprintf(stderr,"This is in the operation %s\n",name); my_error("dyv data structure error"); } } void check_dym_access(dym *d,int i, int j,char *name) { check_dym_code(d,name); if ( i < 0 || i >= d->rows || j < 0 || j >= d->cols ) { fprintf(stderr,"In operation \"%s\"\n",name); fprintf(stderr,"the dym (dynamic matrix) has rows = %d, cols = %d\n", d->rows,d->cols ); fprintf(stderr,"You tried to use indices i=%d j=%d\n",i,j); fprintf(stderr,"Here is the dym that was involved:\n"); fprintf_dym(stderr,"dm",d,"\n"); my_error("check_dym_access"); } } void check_dyv_access(dyv *d,int i, char *name) { check_dyv_code(d,name); if ( i < 0 || i >= d->size ) { fprintf(stderr,"In operation \"%s\"\n",name); fprintf(stderr,"the dyv (dynamic vector) has size = %d\n",d->size); fprintf(stderr,"You tried to use index i=%d\n",i); fprintf(stderr,"Here is the dyv that was involved:\n"); fprintf_dyv(stderr,"dv",d,"\n"); my_error("check_dyv_access"); } } void assert_dym_shape(dym *d,int rows, int cols,char *name) { check_dym_code(d,name); if ( rows != d->rows || cols != d->cols ) { fprintf(stderr,"In operation \"%s\"\n",name); fprintf(stderr,"the dym (dynamic matrix) has rows = %d, cols = %d\n", d->rows,d->cols ); fprintf(stderr,"But should have been predefined with the shape:\n"); fprintf(stderr,"rows = %d, cols = %d\n",rows,cols); my_error("assert_dym_shape"); } } void assert_dyv_shape(dyv *d,int size,char *name) { check_dyv_code(d,name); if ( size != d->size ) { fprintf(stderr,"In operation \"%s\"\n",name); fprintf(stderr,"the dyv (dynamic vector) has size = %d\n", d->size); fprintf(stderr,"But should have been predefined with the shape:\n"); fprintf(stderr,"size = %d\n",size); my_error("assert_dyv_shape"); } } int Dyms_mallocked = 0; int Dyms_freed = 0; int Dyvs_mallocked = 0; int Dyvs_freed = 0; dym *mk_dym(int rows,int cols) { dym *result = AM_MALLOC(dym); result -> dym_code = DYM_CODE; result -> rows = rows; result -> cols = cols; result -> rows_allocated = rows; result -> tdarr = am_malloc_2d_floats(rows,cols); Dyms_mallocked += 1; return(result); } dyv *mk_dyv(int size) { dyv *result = AM_MALLOC(dyv); result -> dyv_code = DYV_CODE; result -> size = size; result -> farr = am_malloc_floats(size); Dyvs_mallocked += 1; return(result); } void free_dym(dym *d) { int i; check_dym_code(d,"free_dym"); d -> dym_code = 7777; for ( i = 0 ; i < d -> rows ; i++ ) { am_free_floats(d->tdarr[i],d->cols); d->tdarr[i] = NULL; } am_free((char *) (d->tdarr),sizeof(float_ptr) * d->rows_allocated); am_free((char *)d,sizeof(dym)); Dyms_freed += 1; } void free_dyv(dyv *d) { check_dyv_code(d,"free_dyv"); d -> dyv_code = 7777; am_free_floats(d->farr,d->size); am_free((char *)d,sizeof(dyv)); Dyvs_freed += 1; } void dym_malloc_report(void) { printf("Dynamic Matrices (datatype dym) currently allocated: %d\n", Dyms_mallocked - Dyms_freed ); printf(" Number of dym allocations since program start: %d\n", Dyms_mallocked ); printf(" Number of dym frees since program start: %d\n\n", Dyms_freed ); printf("Dynamic Vectors (datatype dyv) currently allocated: %d\n", Dyvs_mallocked - Dyvs_freed ); printf(" Number of dyv allocations since program start: %d\n", Dyvs_mallocked ); printf(" Number of dyv frees since program start: %d\n", Dyvs_freed ); } void add_row(dym *d) { check_dym_code(d,"add_row"); if ( d->rows_allocated < d->rows ) my_error("oujbcowlbucv"); if ( d->rows_allocated == d->rows ) { int new_size = int_max(100,(int) (2.5 * d->rows_allocated)); float **new = AM_MALLOC_ARRAY(float_ptr,new_size); int i; for ( i = 0 ; i < new_size ; i++ ) new[i] = NULL; for ( i = 0 ; i < d->rows ; i++ ) new[i] = d->tdarr[i]; am_free((char *)d->tdarr,d->rows_allocated * sizeof(float_ptr)); d->tdarr = new; d->rows_allocated = new_size; } d->tdarr[d->rows] = am_malloc_floats(d->cols); set_floats_constant(d->tdarr[d->rows],d->cols,0.0); d->rows += 1; } float dym_ref(dym *d, int i,int j) { check_dym_access(d,i,j,"dym_ref"); return(d->tdarr[i][j]); } float dyv_ref(dyv *d, int i) { check_dyv_access(d,i,"dyv_ref"); return(d->farr[i]); } void dym_set(dym *d,int i,int j,float value) { check_dym_access(d,i,j,"dym_set"); d->tdarr[i][j] = value; } void dyv_set(dyv *d,int i,float value) { check_dyv_access(d,i,"dyv_set"); d->farr[i] = value; } void dym_increment(dym *d,int i,int j,float value) { check_dym_access(d,i,j,"dym_increment"); d->tdarr[i][j] += value; } void dyv_increment(dyv *d,int i,float value) { check_dyv_access(d,i,"dyv_increment"); d->farr[i] += value; } void dyv_increase_length(dyv *d,int extra_size) { float *longer = am_malloc_floats(d->size + extra_size); copy_floats(d->farr,longer,d->size); am_free_floats(d->farr,d->size); d->farr = longer; d->size += extra_size; } void copy_dym_to_tdarr(dym *d,float **tdarr) { copy_2d_floats(d->tdarr,tdarr,d->rows,d->cols); } float **mk_tdarr_from_dym(dym *d) { float **result; check_dym_code(d,"make_copy_tdarr"); result = am_malloc_2d_floats(d->rows,d->cols); copy_dym_to_tdarr(d,result); return(result); } void copy_dyv_to_farr(dyv *d, float *farr) { copy_floats(d->farr,farr,d->size); } float *mk_farr_from_dyv(dyv *d) { float *result; check_dyv_code(d,"make_copy_farr"); result = am_malloc_floats(d->size); copy_dyv_to_farr(d,result); return(result); } void copy_tdarr_to_dym(float **tdarr,int rows,int cols,dym *r_d) { assert_dym_shape(r_d,rows,cols,"copy_tdarr_to_dym"); copy_2d_floats(tdarr,r_d->tdarr,rows,cols); } dym *mk_dym_from_tdarr(float **tdarr,int rows,int cols) { dym *result = mk_dym(rows,cols); copy_tdarr_to_dym(tdarr,rows,cols,result); return(result); } void copy_farr_to_dyv(float *farr,int size,dyv *r_d) { assert_dyv_shape(r_d,size,"copy_farr_to_dyv"); copy_floats(farr,r_d->farr,size); } dyv *mk_dyv_from_farr(float *farr,int size) { dyv *result = mk_dyv(size); copy_farr_to_dyv(farr,size,result); return(result); } /***** Copying dyvs to and from rows and columns of dyms. And Making dyvs from rows and columns of dyms too. *********/ void copy_dyv_to_dym_row(dyv *dv,dym *dm,int row) { int i; check_dym_code(dm,"copy_dyv_to_dym_row"); if ( row < 0 || row >= dm->rows ) my_error("copy_dyv_to_dym_row: illegal row"); assert_dyv_shape(dv,dm->cols,"copy_dyv_to_dym_row"); for ( i = 0 ; i < dm->cols ; i++ ) dm->tdarr[row][i] = dv->farr[i]; } void copy_dyv_to_dym_col(dyv *dv,dym *dm,int col) { int i; check_dym_code(dm,"copy_dyv_to_dym_col"); if ( col < 0 || col >= dm->cols ) my_error("copy_dyv_to_dym_col: illegal col"); assert_dyv_shape(dv,dm->rows,"copy_dyv_to_dym_col"); for ( i = 0 ; i < dm->rows ; i++ ) dm->tdarr[i][col] = dv->farr[i]; } void copy_dym_row_to_dyv(dym *dm,dyv *dv,int row) { int i; check_dym_code(dm,"copy_dym_row_to_dyv"); if ( row < 0 || row >= dm->rows ) my_error("copy_dyv_to_dym_row: illegal row"); assert_dyv_shape(dv,dm->cols,"copy_dyv_to_dym_row"); for ( i = 0 ; i < dm->cols ; i++ ) dv->farr[i] = dm->tdarr[row][i]; } dyv *mk_dyv_from_dym_row(dym *dm,int row) { dyv *result = mk_dyv(dm->cols); copy_dym_row_to_dyv(dm,result,row); return(result); } void copy_dym_col_to_dyv(dym *dm,dyv *dv,int col) { int i; check_dym_code(dm,"copy_dym_col_to_dyv"); if ( col < 0 || col >= dm->cols ) my_error("copy_dym_col_to_dyv: illegal col"); assert_dyv_shape(dv,dm->rows,"copy_dym_col_to_dyv"); for ( i = 0 ; i < dm->rows ; i++ ) dv->farr[i] = dm->tdarr[i][col]; } dyv *mk_dyv_from_dym_col(dym *dm,int col) { dyv *result = mk_dyv(dm->rows); copy_dym_col_to_dyv(dm,result,col); return(result); } /***** Copying farrs to and from rows and columns of dyms. And Making farrs from rows and columns of dyms too. *********/ void copy_farr_to_dym_row(float *farr,dym *dm,int row) { int i; check_dym_code(dm,"copy_farr_to_dym_row"); if ( row < 0 || row >= dm->rows ) my_error("copy_farr_to_dym_row: illegal row"); for ( i = 0 ; i < dm->cols ; i++ ) dm->tdarr[row][i] = farr[i]; } void copy_farr_to_dym_col(float *farr,dym *dm,int col) { int i; check_dym_code(dm,"copy_farr_to_dym_col"); if ( col < 0 || col >= dm->cols ) my_error("copy_farr_to_dym_col: illegal col"); for ( i = 0 ; i < dm->rows ; i++ ) dm->tdarr[i][col] = farr[i]; } void copy_dym_row_to_farr(dym *dm,float *farr,int row) { int i; check_dym_code(dm,"copy_dym_row_to_farr"); if ( row < 0 || row >= dm->rows ) my_error("copy_farr_to_dym_row: illegal row"); for ( i = 0 ; i < dm->cols ; i++ ) farr[i] = dm->tdarr[row][i]; } float *mk_farr_from_dym_row(dym *dm,int row) { float *result = am_malloc_floats(dm->cols); copy_dym_row_to_farr(dm,result,row); return(result); } void copy_dym_col_to_farr(dym *dm,float *farr,int col) { int i; check_dym_code(dm,"copy_dym_col_to_farr"); if ( col < 0 || col >= dm->cols ) my_error("copy_dym_col_to_farr: illegal col"); for ( i = 0 ; i < dm->rows ; i++ ) farr[i] = dm->tdarr[i][col]; } float *mk_farr_from_dym_col(dym *dm,int col) { float *result = am_malloc_floats(dm->rows); copy_dym_col_to_farr(dm,result,col); return(result); } /***** Copying dyvs to and from rows and columns of tdarrs. And Making dyvs from rows and columns of tdarrs too. *********/ void copy_dyv_to_tdarr_row(dyv *dv,float **tdarr,int row) { int i; for ( i = 0 ; i < dv->size ; i++ ) tdarr[row][i] = dv->farr[i]; } void copy_dyv_to_tdarr_col(dyv *dv,float **tdarr,int col) { int i; for ( i = 0 ; i < dv->size ; i++ ) tdarr[i][col] = dv->farr[i]; } void copy_tdarr_row_to_dyv(float **tdarr,dyv *dv,int row) { int i; for ( i = 0 ; i < dv->size ; i++ ) dv->farr[i] = tdarr[row][i]; } dyv *mk_dyv_from_tdarr_row(float **tdarr,int row,int tdarr_cols) { dyv *result = mk_dyv(tdarr_cols); copy_tdarr_row_to_dyv(tdarr,result,row); return(result); } void copy_tdarr_col_to_dyv(float **tdarr,dyv *dv,int col) { int i; for ( i = 0 ; i < dv->size ; i++ ) dv->farr[i] = tdarr[i][col]; } dyv *mk_dyv_from_tdarr_col(float **tdarr,int col,int tdarr_rows) { dyv *result = mk_dyv(tdarr_rows); copy_tdarr_col_to_dyv(tdarr,result,col); return(result); } /**** Making whole dyms from dyvs ******/ dym *mk_col_dym_from_dyv(dyv *dv) { dym *result; check_dyv_code(dv,"mk_col_dym_from_dyv"); result = mk_dym(dv->size,1); copy_dyv_to_dym_col(dv,result,0); return(result); } dym *mk_row_dym_from_dyv(dyv *dv) { dym *result; check_dyv_code(dv,"mk_row_dym_from_dyv"); result = mk_dym(1,dv->size); copy_dyv_to_dym_row(dv,result,0); return(result); } dym *mk_diag_dym_from_dyv(dyv *dv) { dym *result; int i; check_dyv_code(dv,"mk_row_dym_from_dyv"); result = mk_dym(dv->size,dv->size); zero_dym(result); for ( i = 0 ; i < dv->size ; i++ ) result->tdarr[i][i] = dv->farr[i]; return(result); } /**** Making whole dyms from farrs ******/ dym *mk_col_dym_from_farr(float *farr,int farr_size) { dym *result; result = mk_dym(farr_size,1); copy_farr_to_dym_col(farr,result,0); return(result); } dym *mk_row_dym_from_farr(float *farr, int farr_size) { dym *result; result = mk_dym(1,farr_size); copy_farr_to_dym_row(farr,result,0); return(result); } dym *mk_diag_dym_from_farr(float *farr, int farr_size) { dym *result; int i; result = mk_dym(farr_size,farr_size); zero_dym(result); for ( i = 0 ; i < farr_size ; i++ ) result->tdarr[i][i] = farr[i]; return(result); } /***** Simple operations on dyms ******/ int dym_rows(dym *d) { check_dym_code(d,"dym_rows"); return(d->rows); } int dym_cols(dym *d) { check_dym_code(d,"dym_cols"); return(d->cols); } int dyv_size(dyv *d) { check_dyv_code(d,"dyv_size"); return(d->size); } int num_char_occurs(char *s,char c) { int result = 0; int i; for ( i = 0 ; s[i] != '\0' ; i++ ) if ( s[i] == c ) result += 1; return(result); } dym *old_formatted_dym_from_lex(lex *lx,char *fstring,char fcode,char *fname) { int cols = num_char_occurs(fstring,fcode); dym *d = mk_dym(0,cols); int i; for ( i = 0 ; i < lx -> num_lines ; i++ ) { int line_len = lex_line_length(lx,i); if ( line_len > 0 && lex_line_ref(lx,i,0)->type == NUMBER ) { int dym_col = 0; int lex_col; add_row(d); for ( lex_col = 0 ; lex_col < line_len ; lex_col++ ) { if ( fstring[lex_col] == fcode ) { lextok *lxt = lex_line_ref(lx,i,lex_col); if ( lxt->type != NUMBER ) { fprintf(stderr,"Syntax error attempting to read a number from\n"); fprintf(stderr,"file \"%s\". Line %d item %d should be number.\n", fname,i+1,lex_col+1); my_error("Non-number in datafile"); } dym_set(d,dym_rows(d)-1,dym_col,lxt->number); dym_col += 1; } } if ( dym_col != cols ) { int j; fprintf(stderr,"Syntax error attempting to read numbers from\n"); fprintf(stderr,"file \"%s\". Line %d should contain %d number%s.\n", fname,i+1,cols,(cols==1)?"":"s"); fprintf(stderr,"There should be numbers in following line items:\n"); for ( j = 0 ; fstring[j] != '\0' ; j++ ) if ( fstring[j] == fcode ) fprintf(stderr," %d ",j+1); fprintf(stderr,"\n"); my_error("Too short a line in datafile"); } } } return(d); } int items_on_first_numeric_row(lex *lx) { int result = 0; int i; for ( i = 0 ; i < lx->num_lines && result == 0 ; i++ ) if ( lex_line_length(lx,i) > 0 && lex_line_ref(lx,i,0)->type == NUMBER ) result = lex_line_length(lx,i); return(result); } char *make_fstring_for_everything(lex *lx) { int cols = items_on_first_numeric_row(lx); char *result = AM_MALLOC_ARRAY(char,1 + cols); int i; for ( i = 0 ; i < cols ; i++ ) result[i] = 'i'; result[i] = '\0'; return(result); } dym *old_read_dym(FILE *s,char *fname,char *format) { dym *d; lex lx; char *my_format; lex_read_stream(s,&lx); my_format = make_fstring_for_everything(&lx); if ( Verbosity > 1000.0 ) printf("my_format = %s\n",my_format); d = old_formatted_dym_from_lex(&lx,(format==NULL) ? my_format : format, 'i',fname ); lex_free_contents(&lx); am_free_string(my_format); return(d); } void old_read_io_dyms( FILE *s, char *fname, char *format, dym **r_in_dym, dym **r_out_dym ) { lex lx; char *my_format; lex_read_stream(s,&lx); my_format = make_fstring_for_everything(&lx); if ( strlen(my_format) > 0 ) my_format[strlen(my_format)-1] = 'o'; *r_in_dym = old_formatted_dym_from_lex(&lx,(format==NULL) ? my_format : format, 'i',fname ); *r_out_dym = old_formatted_dym_from_lex(&lx,(format==NULL) ? my_format : format, 'o',fname ); lex_free_contents(&lx); am_free_string(my_format); } void basic_message(char *fname,int lines) { printf(" (Read %6d lines from file %s)\n",lines,fname); } void basic_read_io_dyms( FILE *s, char *fname, char *format, bool last_col_is_output, dym **r_in_dym, dym **r_out_dym ) /* last_col_is_output is IGNORED except when format == NULL In this function format may be NULL or else a null-terminated string. If format is NULL and last_col_is_output is TRUE then you can treat the specification as though the function were called with format = iii..iio where the number of i's is one less than the number of numbers on the first numerical line of the file. If format is NULL and last_col_is_output is FALSE then you can treat the specification as though the function were called with format = iii..iii where the number of i's is the number of numbers on the first numerical line of the file. Let N be the number of characters in format. Then we read the file, assuming that every line which starts with a number as its first lexical item contains exactly N lexical items all numbers. Otherwise we'll signal a syntax error. What do we do with these numbers? The number on the i'th numerical row of the file, in the j'th numeric colum is either ignored, stored in *r_in_dym or stored in *r_out_dym. It is stored in *r_in_dym[i,k] if format[j] is the k'th 'i' character in format. It is stored in *r_out_dym[i,k] if format[j] is the k'th 'o' character in format. If format contains no i's , no dym is created, and *r_in_dym is set to NULL If format contains no o's , no dym is created, and *r_out_dym is set to NULL EXAMPLE: FILE STARTS HERE: .7 6 -9 2.1 # Line starts with non-numeric so ignored Actually, this ignored too -1 3 4 3 If called with format = ii-o would produce *r_in_dym = [ 0.7 6.0 ] *r_out_dym = [ 2.1 ] [ -1.0 3.0 ] [ 3.0 ] If called with format = --i- would produce *r_in_dym = [ -9.0 ] *r_out_dym = NULL [ 4.0 ] If called with format = o-io would produce *r_in_dym = [ -9.0 ] *r_out_dym = [ 0.7 2.1 ] [ 4.0 ] [ -1.0 3.0 ] If called with format = iiio would produce *r_in_dym = [ 0.7 6.0 -9.0 ] *r_out_dym = [ 2.1 ] [ -1.0 3.0 4.0 ] [ 3.0 ] If called with format = NULL and last_col_is_output == TRUE would also produce *r_in_dym = [ 0.7 6.0 -9.0 ] *r_out_dym = [ 2.1 ] [ -1.0 3.0 4.0 ] [ 3.0 ] If called with format = iiii would produce *r_in_dym = [ 0.7 6.0 -9.0 2.1 ] *r_out_dym = NULL [ -1.0 3.0 4.0 3.0 ] If called with format = NULL and last_col_is_output == FALSE would also produce *r_in_dym = [ 0.7 6.0 -9.0 2.1 ] *r_out_dym = NULL [ -1.0 3.0 4.0 3.0 ] If called with format = o-ioo would produce ERROR because numeric lines don't contain 5 numbers. */ { lex lx[1]; bool file_ended = FALSE; dym *din = NULL; dym *dout = NULL; int line_number = 0; int din_cols = (format==NULL) ? 0 : num_char_occurs(format,'i'); int dout_cols = (format==NULL) ? 0 : num_char_occurs(format,'o'); int total_cols = (format==NULL) ? 0 : strlen(format); while ( !file_ended ) { file_ended = !lex_can_read_line(s,lx); line_number += 1; if ( !file_ended ) { int line_len = lex_line_length(lx,0); if ( (line_number % 1000)==0 ) basic_message(fname,line_number); if ( line_len > 0 && lex_line_ref(lx,0,0)->type == NUMBER ) { int din_col = 0; int dout_col = 0; int lex_col; if ( format == NULL ) { dout_cols = (last_col_is_output) ? 1 : 0; din_cols = line_len - dout_cols; total_cols = din_cols + dout_cols; } if ( din == NULL && din_cols > 0 ) din = mk_dym(0,din_cols); if ( dout == NULL && dout_cols > 0 ) dout = mk_dym(0,dout_cols); if ( din_cols > 0 ) add_row(din); if ( dout_cols > 0 ) add_row(dout); for ( lex_col = 0 ; lex_col < int_min(line_len,total_cols) ; lex_col++ ) { bool is_input,is_output; if ( format != NULL ) is_input = format[lex_col] == 'i'; else is_input = lex_col < din_cols; if ( format != NULL ) is_output = format[lex_col] == 'o'; else is_output = lex_col == din_cols; if ( is_input || is_output ) { lextok *lxt = lex_line_ref(lx,0,lex_col); if ( lxt->type != NUMBER ) { fprintf(stderr,"Syntax error attempting to read number from\n"); fprintf(stderr,"file \"%s\". Line %d item %d should be number\n", fname,line_number,lex_col+1); my_error("Non-number in datafile"); } if ( is_input ) { dym_set(din,dym_rows(din)-1,din_col,lxt->number); din_col += 1; } else { dym_set(dout,dym_rows(dout)-1,dout_col,lxt->number); dout_col += 1; } } } if ( line_len != total_cols ) { int j; fprintf(stderr,"Syntax error attempting to read numbers from\n"); fprintf(stderr,"file \"%s\". Line %d should contain %d number%s.\n", fname,line_number,total_cols,(total_cols==1)?"":"s"); fprintf(stderr,"There should be numbers in following line items:\n"); for ( j = 0 ; j < total_cols ; j++ ) if ( format == NULL || format[j] == 'i' || format[j] =='o' ) fprintf(stderr," %d ",j+1); fprintf(stderr,"\n"); my_error("Too short a line in datafile"); } } } lex_free_contents(lx); } if ( din == NULL && dout == NULL ) { if ( format == NULL ) { fprintf(stderr,"Syntax error attempting to read numbers from\n"); fprintf(stderr,"file \"%s\". There are no lines containing numbers\n", fname ); fprintf(stderr,"and since you didn't give me any knowledge of the\n"); fprintf(stderr,"of the number of columns in the matrix/matrices\n"); fprintf(stderr,"required for loading, I can't just set them to\n"); fprintf(stderr,"matrices with zero rows without knowing the number\n"); fprintf(stderr,"of columns.\n"); my_error("In reading a dym (or dyms)"); } else { if ( din_cols > 0 ) din = mk_dym(0,din_cols); if ( dout_cols > 0 ) dout = mk_dym(0,dout_cols); } } if ( line_number > 1000 ) basic_message(fname,line_number); *r_in_dym = din; *r_out_dym = dout; } void read_io_dyms( FILE *s, char *fname, char *format, dym **r_in_dym, dym **r_out_dym ) /* In this function format may be NULL or else a null-terminated string. If format is NULL then you can treat the specification as though the function were called with format = iii..iio where the number of i's is one less than the number of numbers on the first numerical line of the file. Let N be the number of characters in format. Then we read the file, assuming that every line which starts with a number as its first lexical item contains exactly N lexical items all numbers. Otherwise we'll signal a syntax error. What do we do with these numbers? The number on the i'th numerical row of the file, in the j'th numeric colum is either ignored, stored in *r_in_dym or stored in *r_out_dym. It is stored in *r_in_dym[i,k] if format[j] is the k'th 'i' character in format. It is stored in *r_out_dym[i,k] if format[j] is the k'th 'o' character in format. If format contains no i's , no dym is created, and *r_in_dym is set to NULL If format contains no o's , no dym is created, and *r_out_dym is set to NULL EXAMPLE: FILE STARTS HERE: .7 6 -9 2.1 # Line starts with non-numeric so ignored Actually, this ignored too -1 3 4 3 If called with format = ii-o would produce *r_in_dym = [ 0.7 6.0 ] *r_out_dym = [ 2.1 ] [ -1.0 3.0 ] [ 3.0 ] If called with format = --i- would produce *r_in_dym = [ -9.0 ] *r_out_dym = NULL [ 4.0 ] If called with format = o-io would produce *r_in_dym = [ -9.0 ] *r_out_dym = [ 0.7 2.1 ] [ 4.0 ] [ -1.0 3.0 ] If called with format = iiio would produce *r_in_dym = [ 0.7 6.0 -9.0 ] *r_out_dym = [ 2.1 ] [ -1.0 3.0 4.0 ] [ 3.0 ] If called with format = NULL would also produce *r_in_dym = [ 0.7 6.0 -9.0 ] *r_out_dym = [ 2.1 ] [ -1.0 3.0 4.0 ] [ 3.0 ] If called with format = o-ioo would produce ERROR because numeric lines don't contain 5 numbers. */ { bool last_col_is_output = TRUE; basic_read_io_dyms(s,fname,format,last_col_is_output,r_in_dym,r_out_dym); } dym *read_dym(FILE *s,char *fname,char *format) /* In this function format may be NULL or else a null-terminated string. If format is NULL then you can treat the specification as though the function were called with format = iii..iii where the number of i's is the number of numbers on the first numerical line of the file. Let N be the number of characters in format. Then we read the file, assuming that every line which starts with a number as its first lexical item contains exactly N lexical items all numbers. Otherwise we'll signal a syntax error. What do we do with these numbers? The number on the i'th numerical row of the file, in the j'th numeric colum is either ignored, stored in dym *result (the dym we make and return) It is stored in result[i,k] if format[j] is the k'th 'i' character in format. If format contains no i's , no dym is created, and *r_in_dym is set to NULL EXAMPLE: FILE STARTS HERE: .7 6 -9 2.1 # Line starts with non-numeric so ignored Actually, this ignored too -1 3 4 3 If called with format = ii-- would produce result = [ 0.7 6.0 ] [ -1.0 3.0 ] If called with format = --i- would produce result = [ -9.0 ] [ 4.0 ] If called with format = --i- would produce result = [ -9.0 ] [ 4.0 ] If called with format = iii- would produce result = [ 0.7 6.0 -9.0 ] [ -1.0 3.0 4.0 ] If called with format = iiii would produce result = [ 0.7 6.0 -9.0 2.1 ] [ -1.0 3.0 4.0 3.0 ] If called with format = NULL would also produce result = [ 0.7 6.0 -9.0 2.1 ] [ -1.0 3.0 4.0 3.0 ] If called with format = o-ioo would produce ERROR because numeric lines don't contain 5 numbers. */ { dym *din,*dout; bool last_col_is_output = FALSE; basic_read_io_dyms(s,fname,format,last_col_is_output,&din,&dout); if ( dout != NULL ) free_dym(dout); /* Need to do this if someone gave us a format with 'o's in it. We are meant to ignore outputs in this function. */ return(din); } void save_dym(FILE *s,dym *d) { int i,j; for ( i = 0 ; i < dym_rows(d) ; i++ ) for ( j = 0 ; j < dym_cols(d) ; j++ ) fprintf(s,"%12g%s",dym_ref(d,i,j), (j==dym_cols(d)-1) ? "\n" : " " ); } void dym_main(int argc,char *argv[]) { char *fname = string_from_args("fname",argc,argv,"default.txt"); FILE *s = safe_fopen(fname,"r"); char *format = string_from_args("-fo",argc,argv,(char *)NULL); if ( index_of_arg("io",argc,argv) > 0 ) { dym *in; dym *out; read_io_dyms(s,fname,format,&in,&out); fprintf_dym(stdout,"in",in,"\n"); wait_for_key(); fprintf_dym(stdout,"out",out,"\n"); wait_for_key(); free_dym(in); free_dym(out); } else { dym *d = read_dym(s,fname,format); if ( input_bool("print out the dym>") ) { fprintf_dym(stdout,"d",d,"\n"); wait_for_key(); } free_dym(d); } am_malloc_report(); } #define DYM_SVD_NULL_THRESH (1e-18) /* NOTE: ***SEE DYM NOTES IN AMDM.H ***** */ void constant_dym(dym *r_d,float v) { int i,j; check_dym_code(r_d,"constant_dym"); for ( i = 0 ; i < r_d -> rows ; i++ ) for ( j = 0 ; j < r_d -> cols ; j++ ) r_d->tdarr[i][j] = v; } void constant_dyv(dyv *r_d,float v) { check_dyv_code(r_d,"constant_dyv"); set_floats_constant(r_d->farr,r_d->size,v); } void zero_dym(dym *r_d) { check_dym_code(r_d,"zero_dym"); constant_dym(r_d,0.0); } void zero_dyv(dyv *r_d) { check_dyv_code(r_d,"zero_dyv"); constant_dyv(r_d,0.0); } dym *mk_constant_dym(int rows,int cols,float v) { dym *result = mk_dym(rows,cols); constant_dym(result,v); return(result); } dyv *mk_constant_dyv(int size,float v) { dyv *result = mk_dyv(size); constant_dyv(result,v); return(result); } dym *mk_zero_dym(int rows,int cols) { dym *result = mk_dym(rows,cols); zero_dym(result); return(result); } dyv *mk_zero_dyv(int size) { dyv *result = mk_dyv(size); zero_dyv(result); return(result); } /**** Standard operations on dyms ****/ void dym_scalar_mult(dym *d, float alpha, dym *r_d) { int i,j; assert_dym_shape(r_d,d->rows,d->cols,"dym_scalar_mult"); for ( i = 0 ; i < r_d -> rows ; i++ ) for ( j = 0 ; j < r_d -> cols ; j++ ) r_d -> tdarr[i][j] = d->tdarr[i][j] * alpha; } dym *mk_dym_scalar_mult(dym *d,float alpha) { dym *result; check_dym_code(d,"mk_dym_scalar_mult"); result = mk_dym(d->rows,d->cols); dym_scalar_mult(d,alpha,result); return(result); } void copy_dym(dym *d, dym *r_d) { assert_dym_shape(r_d,d->rows,d->cols,"copy_dym"); dym_scalar_mult(d,1.0,r_d); } dym *mk_copy_dym(dym *d) { check_dym_code(d,"mk_copy_dym"); return(mk_dym_scalar_mult(d,1.0)); } void dym_plus(dym *d_1, dym *d_2, dym *r_d) { int i,j; if ( d_1 -> rows != d_2 -> rows || d_1 -> cols != d_2 -> cols ) { fprintf_dym(stderr,"d_1",d_1,"\n"); fprintf_dym(stderr,"d_2",d_2,"\n"); my_error("dym_plus: dyms (DYnamic Matrices) different shape"); } assert_dym_shape(r_d,d_1->rows,d_1->cols,"dym_plus"); for ( i = 0 ; i < r_d -> rows ; i++ ) for ( j = 0 ; j < r_d -> cols ; j++ ) r_d -> tdarr[i][j] = d_1->tdarr[i][j] + d_2 -> tdarr[i][j]; } dym *mk_dym_plus(dym *a,dym *b) { dym *result = mk_dym(a->rows,a->cols); dym_plus(a,b,result); return(result); } void dym_subtract(dym *d_1,dym *d_2,dym *r_d) { dym *a = mk_dym_scalar_mult(d_2,-1.0); dym_plus(d_1,a,r_d); free_dym(a); } dym *mk_dym_subtract(dym *a,dym *b) { dym *result = mk_dym(a->rows,a->cols); dym_subtract(a,b,result); return(result); } /********* Standard operations of dyvs *********/ void dyv_scalar_mult(dyv *d, float alpha, dyv *r_d) { int i; assert_dyv_shape(r_d,d->size,"dyv_scalar_mult"); for ( i = 0 ; i < r_d -> size ; i++ ) r_d -> farr[i] = d->farr[i] * alpha; } dyv *mk_dyv_scalar_mult(dyv *d,float alpha) { dyv *result; check_dyv_code(d,"mk_dyv_scalar_mult"); result = mk_dyv(d->size); dyv_scalar_mult(d,alpha,result); return(result); } void copy_dyv(dyv *d, dyv *r_d) { assert_dyv_shape(r_d,d->size,"copy_dyv"); dyv_scalar_mult(d,1.0,r_d); } dyv *mk_copy_dyv(dyv *d) { check_dyv_code(d,"mk_copy_dyv"); return(mk_dyv_scalar_mult(d,1.0)); } void dyv_plus(dyv *d_1, dyv *d_2, dyv *r_d) { int i; check_dyv_code(d_1,"dyv_plus (1st arg)"); check_dyv_code(d_2,"dyv_plus (2nd arg)"); if ( d_1 -> size != d_2 -> size ) { fprintf_dyv(stderr,"d_1",d_1,"\n"); fprintf_dyv(stderr,"d_2",d_2,"\n"); my_error("dyv_plus: dyvs (DYnamic Vectors) different shape"); } assert_dyv_shape(r_d,d_1->size,"dyv_plus"); for ( i = 0 ; i < r_d -> size ; i++ ) r_d -> farr[i] = d_1->farr[i] + d_2 -> farr[i]; } dyv *mk_dyv_plus(dyv *a,dyv *b) { dyv *result = mk_dyv(a->size); dyv_plus(a,b,result); return(result); } void dyv_subtract(dyv *d_1,dyv *d_2,dyv *r_d) { dyv *a = mk_dyv_scalar_mult(d_2,-1.0); check_dyv_code(d_1,"dyv_subtract (1st arg)"); check_dyv_code(d_2,"dyv_subtract (2nd arg)"); if ( d_1 -> size != d_2 -> size ) { fprintf_dyv(stderr,"d_1",d_1,"\n"); fprintf_dyv(stderr,"d_2",d_2,"\n"); my_error("dyv_subtract: dyvs (DYnamic Vectors) different shape"); } assert_dyv_shape(r_d,d_1->size,"dyv_subtract"); dyv_plus(d_1,a,r_d); free_dyv(a); } dyv *mk_dyv_subtract(dyv *a,dyv *b) { dyv *result = mk_dyv(a->size); dyv_subtract(a,b,result); return(result); } /***** More complex operations ******/ float dyv_scalar_product(dyv *a,dyv *b) { int i; float result = 0.0; if ( a->size != b -> size ) { fprintf(stderr,"dyv_scalar_product: sizes differ\n"); wait_for_key(); fprintf_dyv(stderr,"a",a,"\n"); fprintf_dyv(stderr,"b",b,"\n"); my_error("dyv_scalar_product: sizes differ"); } for ( i = 0 ; i < a -> size ; i++ ) result += a->farr[i] * b->farr[i]; return(result); } float dyv_dsqd(dyv *a,dyv *b) { dyv *d = mk_dyv_subtract(a,b); float result = dyv_scalar_product(d,d); free_dyv(d); return(result); } void dym_times_dyv(dym *a,dyv *b,dyv *result) { int i; dyv *temp = mk_dyv(a->rows); /* We need a copy in case b and result share memory */ if ( a->cols != b -> size ) my_error("dym_times_dyv: sizes wrong"); assert_dyv_shape(result,a->rows,"dym_times_dyv"); for ( i = 0 ; i < a->rows ; i++ ) { float sum = 0.0; int j; for ( j = 0 ; j < a->cols ; j++ ) sum += a->tdarr[i][j] * b->farr[j]; temp->farr[i] = sum; } copy_dyv(temp,result); free_dyv(temp); } dyv *mk_dym_times_dyv(dym *a,dyv *b) { dyv *result = mk_dyv(a->rows); dym_times_dyv(a,b,result); return(result); } void dym_mult(dym *d_1, dym *d_2, dym *r_d) { dym *a = mk_dym(d_1->rows,d_2->cols); /* Note we have to first do the multiplying to the result a, in case the routine was called with d_1's memory = r_d's memory or d_2's memory = r_d's memory */ int i,j; if ( d_1 -> cols != d_2 -> rows ) { fprintf_dym(stderr,"d_1",d_1,"\n"); fprintf_dym(stderr,"d_2",d_2,"\n"); my_error("dym_mult: dyms (DYnamic Matrices) wrong shape\n"); } assert_dym_shape(r_d,d_1 -> rows,d_2 -> cols,"dym_mult"); for ( i = 0 ; i < a -> rows ; i++ ) for ( j = 0 ; j < a -> cols ; j++ ) { float sum = 0.0; int k; for ( k = 0 ; k < d_1 -> cols ; k++ ) sum += d_1->tdarr[i][k] * d_2->tdarr[k][j]; a -> tdarr[i][j] = sum; } copy_dym(a,r_d); free_dym(a); } dym *mk_dym_mult(dym *a,dym *b) { dym *result = mk_dym(a->rows,b->cols); dym_mult(a,b,result); return(result); } void dym_transpose(dym *d, dym *r_d) { dym *a = mk_dym(d->cols,d->rows); /* Note we have to first do the transpose to the result a, in case the routine was called with d's memory = r_d's memory */ int i,j; assert_dym_shape(r_d,d->cols,d->rows,"dym_transpose"); for ( i = 0 ; i < d -> rows ; i++ ) for ( j = 0 ; j < d -> cols ; j++ ) a->tdarr[j][i] = d->tdarr[i][j]; copy_dym(a,r_d); free_dym(a); } dym *mk_dym_transpose(dym *a) { dym *result = mk_dym(a->cols,a->rows); dym_transpose(a,result); return(result); } void sing_val_w_check(dyv *w_diag,float thresh,dyv *winv_diag) /* Sets winv_diag[i] = 1 / w_diag[i] for each i, EXCEPT, let w_max = magnitude of w_diag[i] with largest magnitude. Forall w_diag[i] such that w_diag[i] <= thresh * w_max, winv_diag[i] = 0.0 */ { float max_mag = 0.0; float zero_below_me; bool be_verbose = Verbosity > 100.0; int i; assert_dyv_shape(winv_diag,w_diag->size,"sing_val_w_check"); for ( i = 0 ; i < w_diag->size ; i++ ) max_mag = real_max(max_mag,fabs(w_diag->farr[i])); zero_below_me = thresh * max_mag; if ( be_verbose ) printf("dym.c: SVD: thresh=%g, max_mag=%g\n",thresh,max_mag); for ( i = 0 ; i < w_diag->size ; i++ ) { if ( be_verbose ) printf("w_diag->farr[%d] = %g",i,w_diag->farr[i]); if ( w_diag->farr[i] == 0.0 || fabs(w_diag->farr[i]) <= zero_below_me ) { winv_diag->farr[i] = 0.0; if ( be_verbose ) printf("Singular value... zeroing!"); } else winv_diag->farr[i] = 1.0 / w_diag->farr[i]; if ( be_verbose ) printf("\n"); } } void dym_scale_rows(dym *d,dyv *w_diag,dym *r_d) /* Diag(w_diag) * d is copied to r_d */ { int i,j; assert_dym_shape(r_d,d->rows,d->cols,"dym_scale_rows::r_d"); assert_dyv_shape(w_diag,d->rows,"dym_scale_rows::w_diag"); for ( i = 0 ; i < d->rows ; i++ ) for ( j = 0 ; j < d->cols ; j++ ) r_d->tdarr[i][j] = d->tdarr[i][j] * w_diag->farr[i]; } dym *mk_dym_scale_rows(dym *d,dyv *w_diag) /* Returns Diag(w_diag) * d */ { dym *result = mk_dym(d->rows,d->cols); assert_dyv_shape(w_diag,d->rows,"mk_dym_scale_rows::w_diag"); dym_scale_rows(d,w_diag,result); return(result); } void dym_scale_cols(dym *d,dyv *w_diag,dym *r_d) /* d * Diag(w_diag) is copied to r_d */ { int i,j; assert_dym_shape(r_d,d->rows,d->cols,"dym_scale_cols::r_d"); assert_dyv_shape(w_diag,d->cols,"dym_scale_cols::w_diag"); for ( i = 0 ; i < d->rows ; i++ ) for ( j = 0 ; j < d->cols ; j++ ) r_d->tdarr[i][j] = d->tdarr[i][j] * w_diag->farr[j]; } dym *mk_dym_scale_cols(dym *d,dyv *w_diag) /* Returns d * Diag(w_diag) */ { dym *result = mk_dym(d->rows,d->cols); assert_dyv_shape(w_diag,d->cols,"mk_dym_scale_cols::w_diag"); dym_scale_cols(d,w_diag,result); return(result); } void copy_dym_to_farr(dym *d,float *farr) /* Copies either a column dym or a 1-row dym to an array of floats. It is an error to call with neither rows nor columns equal to 1 */ { if ( d->rows != 1 && d->cols != 1 ) { fprintf_dym(stderr,"d",d,"\n"); my_error("dym_to_floats(): should be 1-columns or 1-row"); } else if ( d->rows == 1 ) copy_floats(d->tdarr[0],farr,d->cols); else { int i; for ( i = 0 ; i < d->rows ; i++ ) farr[i] = d->tdarr[i][0]; } } void dym_svd(dym *d,dym *u,dyv *w_diagonal,dym *v) { bool okay; dym *u_temp = mk_dym(d->rows,d->cols); /* In case d and u share memory */ assert_dym_shape(u,d->rows,d->cols,"dym_svd::u"); assert_dyv_shape(w_diagonal,d->cols,"dym_svd::w"); assert_dym_shape(v,d->cols,d->cols,"dym_svd::v"); okay = sing_val_decomp(d->tdarr,d->rows,d->cols, u_temp->tdarr,w_diagonal->farr,v->tdarr); copy_dym(u_temp,u); free_dym(u_temp); if ( !okay ) fprintf(stderr,"amdm.c :: SVD failed\n"); } void make_svd_components(dym *d,dym **u,dyv **w_diagonal,dym **v) { *u = mk_dym(d->rows,d->cols); *w_diagonal = mk_dyv(d->cols); *v = mk_dym(d->cols,d->cols); dym_svd(d,*u,*w_diagonal,*v); } #ifdef NEVER float dym_determinant(dym *dm) { dym *u,*v; dyv *w_diag; float result = 1.0; int i; if ( dym_rows(dm) != dym_cols(dm) ) my_error("dym_determinant: you gave me a non-square dym"); make_svd_components(dm,&u,&w_diag,&v); fprintf_dym(stdout,"dm",dm,"\n"); fprintf_dym(stdout,"u",u,"\n"); fprintf_dyv(stdout,"w_diag",w_diag,"\n"); fprintf_dym(stdout,"v",v,"\n"); for ( i = 0 ; i < dyv_size(w_diag) ; i++ ) result *= dyv_ref(w_diag,i); free_dym(u); free_dym(v); free_dyv(w_diag); if ( dym_cols(dm) == 1 ) result = dym_ref(dm,0,0); else if ( dym_cols(dm) == 2 ) result = dym_ref(dm,0,0) * dym_ref(dm,1,1) - dym_ref(dm,1,0) * dym_ref(dm,0,1); else my_error("Sorry, dym_determinant not working for bigger than 2x2"); return(result); } #endif float **mk_nrecipes_matrix_from_dym(dym *d) { float **tdarr = am_malloc_2d_floats(d->rows+1,d->cols+1); int i,j; for ( i = 0 ; i < d->rows ; i++ ) for ( j = 0 ; j < d->cols ; j++ ) tdarr[i+1][j+1] = d->tdarr[i][j]; return(tdarr); } float *mk_nrecipes_vector_from_dyv(dyv *d) { float *farr = am_malloc_floats(d->size+1); int i; for ( i = 0 ; i < d->size ; i++ ) farr[i+1] = d->farr[i]; return(farr); } void copy_nrecipes_vector_to_dyv(float *nrfarr,dyv *d) { int i; for ( i = 0 ; i < d->size ; i++ ) d->farr[i] = nrfarr[i+1]; } void free_nrecipes_matrix(float **nrtdarr,dym *d) { free_2d_floats(nrtdarr,d->rows+1,d->cols+1); } void free_nrecipes_vector(float *nrfarr,dyv *d) { am_free_floats(nrfarr,d->size+1); } void dym_svd_backsub(dym *u,dyv *w_diag,dym *v,dyv *b,dyv *x) /* Solves A x = b where A = u * Diag(w_diag) * v^T x is the only thing updated. */ { float **nru = mk_nrecipes_matrix_from_dym(u); float *nrw = mk_nrecipes_vector_from_dyv(w_diag); float **nrv = mk_nrecipes_matrix_from_dym(v); float *nrb = mk_nrecipes_vector_from_dyv(b); float *nrx = mk_nrecipes_vector_from_dyv(x); assert_dyv_shape(w_diag,u->cols,"dym_svd_backsub::w_diag"); assert_dym_shape(v,u->cols,u->cols,"dym_svd_backsub::v"); assert_dyv_shape(b,u->rows,"dym_svd_backsub::b"); assert_dyv_shape(x,u->cols,"dym_svd_backsub::x"); if (Verbosity > 2000.0) { fprintf_2d_floats(stdout,"nru",nru,u->rows+1,u->cols+1,"\n"); fprintf_floats(stdout,"nrw",nrw,w_diag->size+1,"\n"); fprintf_2d_floats(stdout,"nrv",nrv,v->rows+1,v->cols+1,"\n"); fprintf_floats(stdout,"nrb",nrb,b->size+1,"\n"); fprintf_floats(stdout,"nrx B4",nrx,x->size+1,"\n"); wait_for_key(); } dym_copy_of_svbksb(nru,nrw,nrv,nrb,nrx,u->rows,u->cols); if (Verbosity > 2000.0) { fprintf_floats(stdout,"nrx AFTER",nrx,x->size+1,"\n"); wait_for_key(); } copy_nrecipes_vector_to_dyv(nrx,x); free_nrecipes_vector(nrw,w_diag); free_nrecipes_vector(nrb,b); free_nrecipes_vector(nrx,x); free_nrecipes_matrix(nru,u); free_nrecipes_matrix(nrv,v); } void dym_solve_vector(dym *a, dyv *b, dyv *r_x) /* Sets r_x so that (in the singular value decomp sense) it is as close as possible to a solution of a * r_x = b */ { dym *u,*v; dyv *w_diag; if ( a -> rows != b -> size ) { fprintf(stderr,"amdm.c :: dym_solve, solving a x = b\n"); fprintf(stderr,"a and b have incompatible numbers of rows\n"); fprintf_dym(stderr,"a",a,"\n"); fprintf_dyv(stderr,"b",b,"\n"); my_error("dym_solve(a,b,r_x)"); } assert_dyv_shape(r_x,a->cols,"dym_solve_vector"); make_svd_components(a,&u,&w_diag,&v); dym_svd_backsub(u,w_diag,v,b,r_x); free_dym(u); free_dym(v); free_dyv(w_diag); } dyv *mk_dym_solve_vector(dym *a, dyv *b) { dyv *result = mk_dyv(a->cols); dym_solve_vector(a,b,result); return(result); } void invert_dym(dym *d,dym *r_d) /* PRECONDITION: d must be a square dym. if d is singular or near singular returns the least squares inverse. */ { dym *u,*v; dyv *w_diag; dyv *winv_diag; check_dym_code(d,"invert_dym"); winv_diag = mk_dyv(d->cols); assert_dym_shape(r_d,d->rows,d->cols,"invert_dym"); if ( d->rows != d->cols ) { fprintf_dym(stderr,"d",d,"\n"); my_error("invert_dym(): the bove dym is not square"); } make_svd_components(d,&u,&w_diag,&v); if ( Verbosity > 2000.0 ) { fprintf_dym(stdout,"u",u,"\n"); wait_for_key(); fprintf_dyv(stdout,"w_diag",w_diag,"\n"); wait_for_key(); fprintf_dym(stdout,"v",v,"\n"); wait_for_key(); } sing_val_w_check(w_diag,DYM_SVD_NULL_THRESH,winv_diag); if ( Verbosity > 2000.0 ) { fprintf_dyv(stdout,"winv_diag",winv_diag,"\n"); wait_for_key(); } dym_transpose(u,r_d); if ( Verbosity > 2000.0 ) { fprintf_dym(stdout,"u^T",r_d,"\n"); wait_for_key(); } dym_scale_rows(r_d,winv_diag,r_d); if ( Verbosity > 2000.0 ) { fprintf_dym(stdout,"Diag(winv_diag) * u^T",r_d,"\n"); wait_for_key(); } dym_mult(v,r_d,r_d); if ( Verbosity > 2000.0 ) { fprintf_dym(stdout,"v * winv * u^T",r_d,"\n"); wait_for_key(); } free_dym(u); free_dyv(w_diag); free_dym(v); free_dyv(winv_diag); } dym *mk_invert_dym(dym *d) { dym *result = mk_dym(d->rows,d->cols); invert_dym(d,result); return(result); } bool is_same_float(float x,float y) { float eps = 1e-3; bool result = TRUE; if ( fabs(x) < eps && fabs(y) < eps ) result = TRUE; else if ( fabs(x) < eps && fabs(y) > 2.0 * eps ) result = FALSE; else if ( fabs(y) < eps && fabs(x) > 2.0 * eps ) result = FALSE; else if ( x < 0.0 && y > 0.0 ) result = FALSE; else if ( fabs(x) > fabs(y) ) result = fabs(x) < (1.0 + eps) * fabs(y); else result = fabs(y) < (1.0 + eps) * fabs(x); return(result); } bool is_dym_symmetric(dym *d) { int i,j; bool result = TRUE; if ( d->rows != d->cols ) { fprintf_dym(stderr,"d = ",d,"\n"); my_error("dym.c is_dym_symmetric(d). Requires a square dym.\n"); } for ( i = 0 ; result && i < d->rows ; i++ ) for ( j = 0 ; result && j < i ; j++ ) { result = is_same_float(d->tdarr[i][j],d->tdarr[j][i]); if ( Verbosity > 70.0 ) printf("i = %d, j = %d, result = %d\n",i,j,result); } if ( !result ) printf("Not Symmetric\n"); return(result); } bool attempt_cholesky_decomp(dym *a,dym *r_l) /* Returns FALSE if a is not symmetric or not positive definite. If it is symmetric and positive definite, returns TRUE and also returns a lower triangular dym in r_l such that r_l r_l^T = a [r_l]_ij = 0 for all i,j such that j > i As usual, it is fine if a and r_l share the same memory. */ { bool result = is_dym_symmetric(a); dym *l = mk_dym(a->rows,a->cols); float maxmat = 1e-6; int i,j; for ( i = 0 ; result && i < a->rows ; i++ ) for ( j = i ; result && j < a->rows ; j++ ) maxmat = real_max(fabs(a->tdarr[i][j]),maxmat); zero_dym(l); if ( Verbosity > 50.0 ) printf("is_dym_symmetric = %d\n",result); if ( result ) { /* COPIED from numerical recipes in C */ int k; for ( i = 0 ; result && i < a->rows ; i++ ) for ( j = i ; result && j < a->rows ; j++ ) { float sum = a->tdarr[i][j]; for ( k = i - 1 ; k >= 0 ; k-- ) sum -= l->tdarr[i][k] * l->tdarr[j][k]; if ( Verbosity > 50.0 ) printf("i = %d , j = %d , sum = %g\n",i,j,sum); if ( i == j ) { if ( sum < -1e-6 * maxmat ) result = FALSE; else if ( sum < 0.0 ) sum = 0.0; else l->tdarr[i][j] = sqrt(sum); } else l->tdarr[j][i] = sum / l->tdarr[i][i]; } } if ( result ) copy_dym(l,r_l); else zero_dym(r_l); free_dym(l); return(result); } bool is_dym_symmetric_positive_definite(dym *d) { dym *l = mk_dym(d->rows,d->cols); bool result = attempt_cholesky_decomp(d,l); free_dym(l); return(result); } void enforce_dym_symmetry(dym *d) /* Let d' = matrix after execution Let d = matrix before execution d' = 0.5 * ( d + d^T ) so d'[i][j] = 0.5 * (d[i][j] + d[j][i]) */ { int i,j; if ( dym_rows(d) != dym_cols(d) ) { fprintf_dym(stderr,"d",d,"\n"); my_error("enforce_dym_symmetry(): the bove dym is not square"); } for ( i = 0 ; i < d->rows ; i++ ) for ( j = 0 ; j < i ; j++ ) dym_set(d,i,j,(dym_ref(d,i,j) + dym_ref(d,j,i))/2.0); for ( i = 0 ; i < d->rows ; i++ ) for ( j = i+1 ; j < d->rows ; j++ ) dym_set(d,i,j,dym_ref(d,j,i)); } void invert_symmetric_dym(dym *d,dym *r_d) /* PRECONDITION: d must be a square symmetric dym. if d is singular or near singular returns the least squares inverse. */ { if ( !is_dym_symmetric(d) ) my_error("invert_symmetric_dym: was given a non-symmetric dym"); else { invert_dym(d,r_d); enforce_dym_symmetry(r_d); } } dym *mk_invert_symmetric_dym(dym *d) { dym *result = mk_dym(d->rows,d->cols); invert_symmetric_dym(d,result); return(result); } void test_dym(int argc,char *argv[]) { char *fname = (argc > 1) ? argv[1] : "default.dym"; FILE *s = safe_fopen(fname,"r"); float farr[200]; int i; dyv *dv = mk_dyv(8); dyv *dv2; dym *a = read_dym(s,fname,NULL); dym *b,*c,*d; fclose(s); for ( i = 0 ; i < 8 ; i++ ) dyv_set(dv,i,(float) i*i); fprintf_dyv(stdout,"dv",dv,"\n"); wait_for_key(); printf("dv . dv = %g\n",dyv_scalar_product(dv,dv)); wait_for_key(); free_dyv(dv); fprintf_dym(stdout,"a",a,"\n"); wait_for_key(); dv = mk_dyv_from_dym_row(a,0); fprintf_dyv(stdout,"mk_dyv_from_dym_row(a,0)",dv,"\n"); wait_for_key(); dv2 = mk_dym_times_dyv(a,dv); fprintf_dyv(stdout,"dv2 := a * dv",dv2,"\n"); wait_for_key(); free_dyv(dv); free_dyv(dv2); dv = mk_dyv(3); dyv_set(dv,0,1.0); dyv_set(dv,1,2.0); dyv_set(dv,2,0.0); dv2 = mk_dym_solve_vector(a,dv); fprintf_dyv(stdout,"dv",dv,"\n"); fprintf_dyv(stdout,"dv2, (where a dv2 = dv)",dv2,"\n"); wait_for_key(); dym_times_dyv(a,dv2,dv); fprintf_dyv(stdout,"dv := a dv2 ",dv,"\n"); wait_for_key(); free_dyv(dv); free_dyv(dv2); d = mk_invert_dym(a); fprintf_dym(stdout,"(d := a^-1)",d,"\n"); wait_for_key(); b = mk_dym_mult(a,d); fprintf_dym(stdout,"(b := a d)",b,"\n"); wait_for_key(); dym_scalar_mult(a,3.0,a); fprintf_dym(stdout,"(a := 3.0 * a)",a,"\n"); wait_for_key(); free_dym(d); d = mk_dym_subtract(a,b); fprintf_dym(stdout,"(d := a - b)",d,"\n"); wait_for_key(); dym_transpose(b,b); fprintf_dym(stdout,"(b := b^T)",b,"\n"); wait_for_key(); dym_mult(b,a,a); fprintf_dym(stdout,"(a := b a)",a,"\n"); wait_for_key(); for ( i = 0 ; i < 10 ; i++ ) farr[i] = sin((float) i); c = mk_col_dym_from_farr(farr,10); fprintf_dym(stdout,"(c := big column)",c,"\n"); wait_for_key(); am_malloc_report(); wait_for_key(); fprintf_dym(stdout,"[after wait_for_key();] a",a,"\n"); wait_for_key(); dym_transpose(a,b); fprintf_dym(stdout,"[after dym_transpose(a,b);] b",b,"\n"); wait_for_key(); dym_plus(a,b,b); fprintf_dym(stdout,"[after dym_plus(a,b,b);] b",b,"\n"); wait_for_key(); if ( !attempt_cholesky_decomp(b,a) ) my_error("cholesky decomp failed"); fprintf_dym(stdout,"[after attempt_cholesky_decomp(b,a)] a",a,"\n"); wait_for_key(); dym_transpose(a,b); fprintf_dym(stdout,"[after dym_transpose(a,b);] b",b,"\n"); wait_for_key(); dym_mult(a,b,b); fprintf_dym(stdout,"[after dym_mult(a,b,b);] b",b,"\n"); wait_for_key(); free_dym(a); free_dym(b); free_dym(c); free_dym(d); am_malloc_report(); } /******* printing dyms (DYnamic Matrices) *********/ typedef struct buftab_struct { int rows,cols; char ***strings; } buftab; char *bufstr(buftab *bt,int i,int j) { char *result; if ( i < 0 || i >= bt->rows || j < 0 || j >= bt->cols ) { result = NULL; my_error("bufstr()"); } else result = bt->strings[i][j]; if ( result == NULL ) result = "-"; return(result); } void print_buftab_debug(buftab *bt) { int i,j; for ( i = 0 ; i < bt->rows ; i++ ) for ( j = 0 ; j < bt->cols ; j++ ) printf("bt[%d][%d] = \"%s\"\n",i,j,bufstr(bt,i,j)); } void fprint_buftab( FILE *s, buftab *bt ) { int *widths = AM_MALLOC_ARRAY(int,bt->cols); int i,j; set_ints_constant(widths,bt->cols,0); for ( i = 0 ; i < bt->rows ; i++ ) for ( j = 0 ; j < bt->cols ; j++ ) widths[j] = int_max(widths[j],strlen(bufstr(bt,i,j))); #ifdef NEVER print_buftab_debug(bt); fprintf_ints(stdout,"widths = ",widths,bt->cols,"\n"); #endif for ( i = 0 ; i < bt->rows ; i++ ) for ( j = 0 ; j < bt->cols ; j++ ) { char ford[20]; sprintf(ford,"%%%ds%s",widths[j],(j==bt->cols-1) ? "\n" : " "); fprintf(s,ford,bufstr(bt,i,j)); } am_free((char *)widths,sizeof(int) * bt->cols); } void init_buftab( buftab *bt, int rows, int cols ) { if ( rows < 0 || cols < 0 ) my_error("init_buftab()"); else { int i,j; bt -> rows = rows; bt -> cols = cols; bt -> strings = AM_MALLOC_ARRAY(char_ptr_ptr,rows); for ( i = 0 ; i < rows ; i++ ) bt->strings[i] = AM_MALLOC_ARRAY(char_ptr,cols); for ( i = 0 ; i < rows ; i++ ) for ( j = 0 ; j < cols ; j++ ) bt->strings[i][j] = NULL; } } void free_buftab_contents(buftab *bt) { int i,j; for ( i = 0 ; i < bt->rows ; i++ ) for ( j = 0 ; j < bt->cols ; j++ ) if ( bt->strings[i][j] != NULL ) am_free((char *)(bt->strings[i][j]), sizeof(char) * (strlen(bt->strings[i][j]) + 1) ); for ( i = 0 ; i < bt->rows ; i++ ) am_free((char *)(bt->strings[i]),sizeof(char_ptr) * bt->cols); am_free((char *)bt->strings,sizeof(char_ptr_ptr) * bt->rows); } void set_buftab( buftab *bt, int i, int j, char *str ) { if ( i < 0 || i >= bt->rows || j < 0 || j >= bt->cols ) my_error("set_buftab()"); else if ( bt->strings[i][j] != NULL ) my_error("set_buftab: non null string"); else bt->strings[i][j] = make_copy_string(str); } void fprintf_dym(FILE *s,char *m1,dym *d,char *m2) { if ( d == NULL ) fprintf(s,"%s = (dym *)NULL%s",m1,m2); else if ( d->dym_code != DYM_CODE ) { fprintf(stderr,"fprintf_dym(s,\"%s\",d,\"\\n\"\n",m1); my_error("fprintf_dym called with a non-allocated dym (DYnamic Matrix)"); } else if ( d->rows <= 0 || d->cols <= 0 ) fprintf(s,"%s = %s", m1,d->rows,(d->rows==-1)?"":"s", d->cols,(d->cols==-1)?"":"s",m2 ); else { int i; buftab bt; init_buftab(&bt,d->rows,d->cols + 4); for ( i = 0 ; i < d->rows ; i++ ) { int j; set_buftab(&bt,i,0,(i == (d->rows-1)/2) ? m1 : ""); set_buftab(&bt,i,1,(i == (d->rows-1)/2) ? "=" : ""); set_buftab(&bt,i,2,"["); for ( j = 0 ; j < d -> cols ; j++ ) { char buff[100]; sprintf(buff," %g ",d->tdarr[i][j]); set_buftab(&bt,i,3+j,buff); } set_buftab(&bt,i,3+d->cols,"]"); } fprint_buftab(s,&bt); free_buftab_contents(&bt); } fprintf(s,"\n"); } void fprintf_dyv(FILE *s,char *m1,dyv *d,char *m2) { if ( d == NULL ) fprintf(s,"%s = (dyv *)NULL%s",m1,m2); else if ( d->dyv_code != DYV_CODE ) { fprintf(stderr,"fprintf_dyv(s,\"%s\",d,\"\\n\"\n",m1); my_error("fprintf_dyv called with a non-allocated dyv (DYnamic Vector)"); } else if ( d->size <= 0 ) fprintf(s,"%s = %s",m1,d->size,m2); else { int i; buftab bt; int cols = 1; init_buftab(&bt,d->size,cols + 4); for ( i = 0 ; i < d->size ; i++ ) { char buff[100]; set_buftab(&bt,i,0,(i == (d->size-1)/2) ? m1 : ""); set_buftab(&bt,i,1,(i == (d->size-1)/2) ? "=" : ""); set_buftab(&bt,i,2,"("); sprintf(buff," %g ",d->farr[i]); set_buftab(&bt,i,3,buff); set_buftab(&bt,i,3+cols,")"); } fprint_buftab(s,&bt); free_buftab_contents(&bt); } fprintf(s,"\n"); } void fprintf_dym_and_confidence( FILE *s, char *m1, dym *d, dym *conf, bool huge_uncertainty, char *m2 ) { if ( d == NULL ) fprintf(s,"%s = (dym *)NULL%s",m1,m2); else if ( d->rows <= 0 || d->cols <= 0 ) fprintf(s,"%s = %s", m1,d->rows,(d->rows==1)?"":"s", d->cols,(d->cols==1)?"":"s",m2 ); else if ( d->rows != conf->rows || d->cols != conf->cols ) my_error("fprintf_dym_and_confidence(). d and conf differ in shape"); else { int i; buftab bt; init_buftab(&bt,d->rows,3 * d->cols + 4); for ( i = 0 ; i < d->rows ; i++ ) { int j; set_buftab(&bt,i,0,(i == (d->rows-1)/2) ? m1 : ""); set_buftab(&bt,i,1,(i == (d->rows-1)/2) ? "=" : ""); set_buftab(&bt,i,2,"["); for ( j = 0 ; j < d -> cols ; j++ ) { char buff[100]; sprintf(buff," %g",d->tdarr[i][j]); set_buftab(&bt,i,3 + 3 * j,buff); set_buftab(&bt,i,4 + 3 * j,"+/-"); if ( huge_uncertainty ) sprintf(buff,"huge uncertainty "); else sprintf(buff,"%g ",conf->tdarr[i][j]); set_buftab(&bt,i,5 + 3 * j,buff); } set_buftab(&bt,i,3 + 3 * d->cols,"]"); } fprint_buftab(s,&bt); free_buftab_contents(&bt); } fprintf(s,"\n"); } void fprintf_dym_dym(FILE *s,char *m1,dym *d1,char *m2,dym *d2,char *m3) { int maxrows = int_max(dym_rows(d1),dym_rows(d2)); int crow = int_min(int_min(dym_rows(d1)-1,dym_rows(d2)-1),(maxrows-1)/2); buftab bt[1]; int i; init_buftab(bt,maxrows,6 + dym_cols(d1) + dym_cols(d2)); for ( i = 0 ; i < maxrows ; i++ ) { int j; if ( i < dym_rows(d1) ) { set_buftab(bt,i,0,(i == crow) ? m1 : ""); set_buftab(bt,i,1,"["); for ( j = 0 ; j < dym_cols(d1) ; j++ ) { char buff[100]; sprintf(buff," %g",dym_ref(d1,i,j)); set_buftab(bt,i,2 + j,buff); } set_buftab(bt,i,2 + dym_cols(d1),"]"); } if ( i < dym_rows(d2) ) { set_buftab(bt,i,3 + dym_cols(d1),(i == crow) ? m2 : ""); set_buftab(bt,i,4 + dym_cols(d1),"["); for ( j = 0 ; j < dym_cols(d2) ; j++ ) { char buff[100]; sprintf(buff," %g",dym_ref(d2,i,j)); set_buftab(bt,i,5 + dym_cols(d1) + j,buff); } set_buftab(bt,i,5 + dym_cols(d1) + dym_cols(d2),"]"); } } fprint_buftab(s,bt); free_buftab_contents(bt); fprintf(s,"\n"); } void fprintf_dym_dyv(FILE *s,char *m1,dym *d1,char *m2,dyv *d2,char *m3) { dym *dm = mk_col_dym_from_dyv(d2); fprintf_dym_dym(s,m1,d1,m2,dm,m3); free_dym(dm); } void boring_test(void) { float f[3]; dyv *v = mk_dyv(4); dym *m = mk_dym(2,3); float **t = am_malloc_2d_floats(4,2); f[0] = 2.0; f[1] = 1.0; f[2] = 3.0; dyv_set(v,0,5.0); dyv_set(v,1,0.0); dyv_set(v,2,6.0); dyv_set(v,3,3.0); dym_set(m,0,0,00.0); dym_set(m,0,1,01.0); dym_set(m,0,2,02.0); dym_set(m,1,0,10.0); dym_set(m,1,1,11.0); dym_set(m,1,2,12.0); t[0][0] = 00.0; t[0][1] = 01.0; t[1][0] = 10.0; t[1][1] = 11.0; t[2][0] = 20.0; t[2][1] = 21.0; t[3][0] = 30.0; t[3][1] = 31.0; fprintf_floats(stdout,"f = ",f,3,"\n"); wait_for_key(); fprintf_dyv(stdout,"v",v,"\n"); wait_for_key(); fprintf_dym(stdout,"m",m,"\n"); wait_for_key(); fprintf_2d_floats(stdout,"t",t,4,2,"\n"); wait_for_key(); fprintf_2d_floats(stdout,"mk_tdarr_from_dym(m)",mk_tdarr_from_dym(m),2,3,"\n"); wait_for_key(); fprintf_floats(stdout,"mk_farr_from_dyv(v) = ",mk_farr_from_dyv(v),4,"\n"); wait_for_key(); fprintf_dym(stdout,"mk_dym_from_tdarr(t,4,2)",mk_dym_from_tdarr(t,4,2),"\n"); wait_for_key(); fprintf_dyv(stdout,"mk_dyv_from_farr(f,3)",mk_dyv_from_farr(f,3),"\n"); wait_for_key(); fprintf_dyv(stdout,"mk_dyv_from_dym_row(m,1)",mk_dyv_from_dym_row(m,1),"\n"); wait_for_key(); fprintf_dyv(stdout,"mk_dyv_from_dym_col(m,2)",mk_dyv_from_dym_col(m,2),"\n"); wait_for_key(); fprintf_floats(stdout,"mk_farr_from_dym_row(m,0) = ",mk_farr_from_dym_row(m,0),3,"\n"); wait_for_key(); fprintf_floats(stdout,"mk_farr_from_dym_col(m,1) = ",mk_farr_from_dym_col(m,1),2,"\n"); wait_for_key(); fprintf_dyv(stdout,"mk_dyv_from_tdarr_row(t,3,2)",mk_dyv_from_tdarr_row(t,3,2),"\n"); wait_for_key(); fprintf_dyv(stdout,"mk_dyv_from_tdarr_col(t,0,4)",mk_dyv_from_tdarr_col(t,0,4),"\n"); wait_for_key(); fprintf_dym(stdout,"mk_col_dym_from_dyv(v)",mk_col_dym_from_dyv(v),"\n"); wait_for_key(); fprintf_dym(stdout,"mk_row_dym_from_dyv(v)",mk_row_dym_from_dyv(v),"\n"); wait_for_key(); fprintf_dym(stdout,"mk_diag_dym_from_dyv(v)",mk_diag_dym_from_dyv(v),"\n"); wait_for_key(); fprintf_dym(stdout,"mk_col_dym_from_farr(f,3)",mk_col_dym_from_farr(f,3),"\n"); wait_for_key(); fprintf_dym(stdout,"mk_row_dym_from_farr(f,3)",mk_row_dym_from_farr(f,3),"\n"); wait_for_key(); fprintf_dym(stdout,"mk_diag_dym_from_farr(f,3)",mk_diag_dym_from_farr(f,3),"\n"); wait_for_key(); } float dyv_mean(dyv *dv) { check_dyv_code(dv,"dyv_mean"); return(floats_mean(dv->farr,dv->size)); } float dyv_sdev(dyv *dv) { float sum_sq = 0.0; int i; float mean = dyv_mean(dv); float result; for ( i = 0 ; i < dv->size ; i++ ) sum_sq += real_square(dv->farr[i] - mean); result = sqrt(sum_sq / int_max(dv->size-1,1)); return(result); } float dyv_min(dyv *dv) { if ( dyv_size(dv) < 1 ) my_error("dyv_min: empty dyv"); return(floats_min(dv->farr,dv->size)); } float dyv_max(dyv *dv) { if ( dyv_size(dv) < 1 ) my_error("dyv_max: empty dyv"); return(floats_max(dv->farr,dv->size)); } int dyv_argmin(dyv *dv) { if ( dyv_size(dv) < 1 ) my_error("dyv_argmin: empty dyv"); return(floats_argmin(dv->farr,dv->size)); } int dyv_argmax(dyv *dv) { if ( dyv_size(dv) < 1 ) my_error("dyv_argmax: empty dyv"); return(floats_argmax(dv->farr,dv->size)); } dyv *mk_dyv_1(float x0) { dyv *result = mk_dyv(1); dyv_set(result,0,x0); return(result); } dyv *mk_dyv_2(float x0,float x1) { dyv *result = mk_dyv(2); dyv_set(result,0,x0); dyv_set(result,1,x1); return(result); } dyv *mk_dyv_3(float x0,float x1,float x2) { dyv *result = mk_dyv(3); dyv_set(result,0,x0); dyv_set(result,1,x1); dyv_set(result,2,x2); return(result); } dym *mk_dym_11(float x00) { dym *result = mk_dym(1,1); dym_set(result,0,0,x00); return(result); } dym *mk_dym_12(float x00, float x01) { dym *result = mk_dym(1,2); dym_set(result,0,0,x00); dym_set(result,0,1,x01); return(result); } dym *mk_dym_13(float x00, float x01, float x02) { dym *result = mk_dym(1,3); dym_set(result,0,0,x00); dym_set(result,0,1,x01); dym_set(result,0,2,x02); return(result); } dym *mk_dym_21(float x00, float x10) { dym *result = mk_dym(2,1); dym_set(result,0,0,x00); dym_set(result,1,0,x10); return(result); } dym *mk_dym_22(float x00, float x01, float x10, float x11) { dym *result = mk_dym(2,2); dym_set(result,0,0,x00); dym_set(result,0,1,x01); dym_set(result,1,0,x10); dym_set(result,1,1,x11); return(result); } dym *mk_dym_23(float x00, float x01, float x02, float x10, float x11, float x12) { dym *result = mk_dym(2,3); dym_set(result,0,0,x00); dym_set(result,0,1,x01); dym_set(result,0,2,x02); dym_set(result,1,0,x10); dym_set(result,1,1,x11); dym_set(result,1,2,x12); return(result); } dym *mk_dym_31(float x00, float x10, float x20) { dym *result = mk_dym(3,1); dym_set(result,0,0,x00); dym_set(result,1,0,x10); dym_set(result,2,0,x20); return(result); } dym *mk_dym_32(float x00, float x01, float x10, float x11, float x20, float x21) { dym *result = mk_dym(3,2); dym_set(result,0,0,x00); dym_set(result,0,1,x01); dym_set(result,1,0,x10); dym_set(result,1,1,x11); dym_set(result,2,0,x20); dym_set(result,2,1,x21); return(result); } dym *mk_dym_33(float x00, float x01, float x02, float x10, float x11, float x12, float x20, float x21, float x22) { dym *result = mk_dym(3,3); dym_set(result,0,0,x00); dym_set(result,0,1,x01); dym_set(result,0,2,x02); dym_set(result,1,0,x10); dym_set(result,1,1,x11); dym_set(result,1,2,x12); dym_set(result,2,0,x20); dym_set(result,2,1,x21); dym_set(result,2,2,x22); return(result); } float dym_xt_a_x_value(dyv *x,dym *a) /* Returns x^T a x */ { dyv *ax = mk_dym_times_dyv(a,x); float result = dyv_scalar_product(x,ax); free_dyv(ax); return(result); } dyv *mk_user_input_dyv(char *message,int dims) { dyv *res = mk_dyv(dims); int i = 0; char buff[100]; for ( i = 0 ; i < dims ; i++ ) { sprintf(buff,"%s (Dyv Component %d)> ",message,i); dyv_set(res,i,input_float(buff)); } return(res); } dym *mk_identity_dym(int dims) { dyv *dv = mk_constant_dyv(dims,1.0); dym *result = mk_diag_dym_from_dyv(dv); free_dyv(dv); return(result); } dyv *mk_basic_dyv_from_args(char *name,int argc,char *argv[],int size) { int index = index_of_arg(name,argc,argv); dyv *result; if ( index < 0 ) result = NULL; else { int i; bool ok = TRUE; result = mk_dyv(size); for ( i = 0 ; i < size && ok ; i++ ) { int j = index + i + 1; if ( j >= argc || !is_a_number(argv[j]) ) ok = FALSE; else dyv_set(result,i,atof(argv[j])); } if ( !ok ) { free_dyv(result); result = NULL; } } return(result); } dyv *mk_dyv_from_args(char *name,int argc,char *argv[],dyv *deflt) /* COPIES in deflt (if so required) */ { bool name_there = index_of_arg(name,argc,argv) >= 0; dyv *result; if ( !name_there ) result = mk_copy_dyv(deflt); else { result = mk_basic_dyv_from_args(name,argc,argv,dyv_size(deflt)); if ( result == NULL ) { fprintf(stderr,"COMMAND LINE USER ERROR (it's YOUR fault)\n"); fprintf(stderr,"...when attempting to read a dyv identified by\n"); fprintf(stderr,"the name \"%s\". Perhaps a non-number, or the\n",name); fprintf(stderr,"command line finished before all args found?\n"); fprintf_dyv(stderr,"deflt_dyv",deflt,"\n"); my_error("mk_dyv_from_args()"); } } return(result); } dym *mk_dym_from_args(char *name,int argc,char *argv[],dym *deflt) /* COPIES in deflt (if so required) */ { bool name_there = index_of_arg(name,argc,argv) >= 0; dym *result; if ( !name_there ) result = mk_copy_dym(deflt); else { int size = dym_rows(deflt) * dym_cols(deflt); dyv *dv = mk_basic_dyv_from_args(name,argc,argv,size); int i,j; if ( dv == NULL ) { fprintf(stderr,"COMMAND LINE USER ERROR (it's YOUR fault)\n"); fprintf(stderr,"...when attempting to read a dym identified by\n"); fprintf(stderr,"the name \"%s\". Perhaps a non-number, or the\n",name); fprintf(stderr,"command line finished before all args found?\n"); fprintf_dym(stderr,"deflt_dym",deflt,"\n"); my_error("mk_dym_from_args()"); } result = mk_dym(dym_rows(deflt),dym_cols(deflt)); for ( i = 0 ; i < dym_rows(deflt) ; i++ ) for ( j = 0 ; j < dym_cols(deflt) ; j++ ) dym_set(result,i,j,dyv_ref(dv,i * dym_cols(deflt) + j)); free_dyv(dv); } return(result); } #define TINY 1.0e-20; bool numrec_ludcmp(float **a, int n, int *indx, float *d) /* Returns TRUE if and only if it succeeded. Returns FALSE if the matrix is singular */ { int i,j,k; int imax = 0; /* Prevent "used uninitialized" warning */ float big,dum,sum,temp; float *vv; bool result = TRUE; vv=am_malloc_floats(n+1); *d=1.0; for (i=1;i<=n && result;i++) { big=0.0; for (j=1;j<=n;j++) if ((temp=fabs(a[i][j])) > big) big=temp; if (big == 0.0) result = FALSE; vv[i]=1.0/big; } for (j=1;j<=n && result;j++) { for (i=1;i= big) { big=dum; imax=i; } } if (j != imax) { for (k=1;k<=n;k++) { dum=a[imax][k]; a[imax][k]=a[j][k]; a[j][k]=dum; } *d = -(*d); vv[imax]=vv[j]; } indx[j]=imax; if (a[j][j] == 0.0) a[j][j]=TINY; if (j != n) { dum=1.0/(a[j][j]); for (i=j+1;i<=n;i++) a[i][j] *= dum; } } am_free_floats(vv,n+1); return(result); } #undef TINY bool is_determinant_positive(dym *a) { int size = dym_rows(a); float **nra = am_malloc_2d_floats(size+1,size+1); int *indx = am_malloc_ints(size+1); float d; int i,j; bool nonsing; bool result; if ( size != dym_cols(a) ) my_error("dym_determinant: a should be square"); for ( i = 0 ; i < size ; i++ ) for ( j = 0 ; j < size ; j++ ) nra[i+1][j+1] = dym_ref(a,i,j); nonsing = numrec_ludcmp(nra,size,indx,&d); if ( nonsing ) { for ( j = 1 ; j <= size ; j++ ) d *= nra[j][j]; result = d >= 0.0; } else result = TRUE; am_free_2d_floats(nra,size+1,size+1); am_free_ints(indx,size+1); return(result); } float dym_determinant(dym *dm) /* Note this is a pretty expensive way to get a determinant. I use SVD to get me the magnitude of the determinant. I then use LU decomp to get me the sign. Why go this way? Why not just use LU decomp for the whole thing? Because it's crrrrap , that's why. det ( 3333 3333 5555 5555 ) = 1.3562 according to numrec's ludcmp. I should probably be able to deduce the sign from the SVD decomposition. If *YOU* are reading this code and *YOU* can figure out how to do it, please let me know. (-awm). */ { dym *u,*v; dyv *w_diag; float result = 1.0; int i; if ( dym_rows(dm) != dym_cols(dm) ) my_error("dym_determinant: you gave me a non-square dym"); make_svd_components(dm,&u,&w_diag,&v); if ( Verbosity > 60.0 ) { fprintf_dym(stdout,"determinant dm",dm,"\n"); fprintf_dym(stdout,"determinant u",u,"\n"); fprintf_dyv(stdout,"determinant w_diag",w_diag,"\n"); fprintf_dym(stdout,"determinant v",v,"\n"); } for ( i = 0 ; i < dyv_size(w_diag) ; i++ ) result *= dyv_ref(w_diag,i); free_dym(u); free_dym(v); free_dyv(w_diag); result = (is_determinant_positive(dm) ? 1.0 : -1.0) * fabs(result); return(result); } void dyv_sort(dyv *dv,dyv *r_dv) { int size; float *farr; check_dyv_code(dv,"dyv_sort (1st arg)"); assert_dyv_shape(r_dv,dv->size,"dyv_sort"); size = dyv_size(dv); farr = mk_farr_from_dyv(dv); sort_floats(farr,size,farr); copy_farr_to_dyv(farr,size,r_dv); am_free_floats(farr,size); } dyv *mk_dyv_sort(dyv *dv) { dyv *result; check_dyv_code(dv,"mk_dyv_sort"); result = mk_dyv(dyv_size(dv)); dyv_sort(dv,result); return(result); } /**** WARNING THE FOLLOWING FUNCTION IS UNIMPLEMENTED AS YET. NUMERICAL RECIPES IS VERY UNCLEAR ABOUT PRECISELY WHAT ludcmp DOES, so our lu_decomp doesn't know what to do *****/ bool attempt_lu_decomp(dym *a,dym *l,dym *u) /* If result is true, the values of dym's l and u are set so that l is lower triangular, u is upper triangular, a = l u Note that a is unchanged. It is legal to pass l or u as a (but not both, of course) in which case the contents of a are overwritten. If result is FALSE then the matrix A was singular. PRECONDISTIONS: a is a square matrix, and has same dimensions as l and u */ { int size = dym_rows(a); float **nra = am_malloc_2d_floats(size+1,size+1); int *indx = am_malloc_ints(size+1); float d; int i,j; bool result; if ( size != dym_cols(a) ) my_error("lu_decomp: a should be square"); assert_dym_shape(l,size,size,"attempt_lu_decomp (l)"); assert_dym_shape(u,size,size,"attempt_lu_decomp (u)"); for ( i = 0 ; i < size ; i++ ) for ( j = 0 ; j < size ; j++ ) nra[i+1][j+1] = dym_ref(a,i,j); result = numrec_ludcmp(nra,size,indx,&d); zero_dym(l); zero_dym(u); if ( result ) { /* for ( i = 0 ; i < size ; i++ ) for ( j = 0 ; j < size ; j++ ) { int true_i = indx[i+1]; int true_j = j+1; if ( true_i > true_j ) */ fprintf_2d_floats(stdout,"nra",nra,size+1,size+1,"\n"); fprintf_ints(stdout,"indx",indx,size+1,"\n"); printf("d = %g\n",d); } else printf("Singular matrix\n"); my_error("not finished"); return(result); } void test_lu_decomp(int argc,char *argv[]) { int dims = int_from_args("dims",argc,argv,2); dym *default_l = mk_identity_dym(dims); dym *default_u = mk_identity_dym(dims); dym *l,*u; dym *a; int i,j; bool okay; for ( i = 0 ; i < dims ; i++ ) for ( j = 0 ; j < dims ; j++ ) { float v = (float) (i + i * j + 1); dym_set((i > j) ? default_l : default_u,i,j,v); } l = mk_dym_from_args("l",argc,argv,default_l); u = mk_dym_from_args("u",argc,argv,default_u); a = mk_dym_mult(l,u); fprintf_dym(stdout,"l",l,"\n"); wait_for_key(); fprintf_dym(stdout,"u",u,"\n"); wait_for_key(); fprintf_dym(stdout,"a",a,"\n"); wait_for_key(); okay = attempt_lu_decomp(a,l,u); fprintf_dym(stdout,"l (after)",l,"\n"); wait_for_key(); fprintf_dym(stdout,"u (after)",u,"\n"); wait_for_key(); dym_mult(l,u,a); fprintf_dym(stdout,"a (after)",a,"\n"); wait_for_key(); free_dym(default_l); free_dym(default_u); free_dym(l); free_dym(u); free_dym(a); am_malloc_report(); } void test_determinant(int argc,char *argv[]) { int dims = int_from_args("dims",argc,argv,2); dym *default_a = mk_identity_dym(dims); dym *a; int i,j; for ( i = 0 ; i < dims ; i++ ) for ( j = 0 ; j < dims ; j++ ) { float v = (int_random(201)-100) * 0.02; dym_set(default_a,i,j,v); } a = mk_dym_from_args("a",argc,argv,default_a); fprintf_dym(stdout,"a",a,"\n"); wait_for_key(); printf("Det(a) = %g\n",dym_determinant(a)); if ( dims == 2 ) printf("ad - bc = %g\n",dym_ref(a,0,0) * dym_ref(a,1,1) - dym_ref(a,1,0) * dym_ref(a,0,1) ); free_dym(default_a); free_dym(a); am_malloc_report(); }