00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "array.h"
00022 #include "helper.h"
00023 #include <stdarg.h>
00024 #include <string.h>
00025 #include <gsl/gsl_rng.h>
00026 #include <time.h>
00027
00046 void* array_max( const Array *a ){
00047 ulong i;
00048 double maxel=DBL_MIN;
00049 double tmp;
00050 void *rmaxel=NULL, *mem=NULL;
00051 for( i=0; i<array_NUMEL(a); i++ ){
00052 mem = array_INDEXMEM1(a,i);
00053 array_dtype_to_double( &tmp, mem, a->dtype );
00054 if( tmp>maxel ){
00055 maxel=tmp;
00056 rmaxel=mem;
00057 }
00058 }
00059 return rmaxel;
00060 }
00061
00080 void* array_min( const Array *a ){
00081 ulong i;
00082 double minel=DBL_MAX;
00083 double tmp;
00084 void *rminel=NULL, *mem=NULL;
00085 for( i=0; i<array_NUMEL(a); i++ ){
00086 mem = array_INDEXMEM1(a,i);
00087 array_dtype_to_double( &tmp, mem, a->dtype );
00088 if( tmp<minel ){
00089 minel=tmp;
00090 rminel=mem;
00091 }
00092 }
00093 return rminel;
00094 }
00095
00108 void array_typecast( Array *a, DType target_type ){
00109 int srcsize=0,destsize=0;
00110 ulong nel;
00111 ulong i;
00112 array_SIZEOF_DTYPE( srcsize, a->dtype );
00113 array_SIZEOF_DTYPE( destsize, target_type );
00114 nel = array_NUMEL( a );
00115
00116 double tmpel;
00117 void *mem;
00118
00119 if( srcsize==destsize ){
00120 a->dtype=target_type;
00121 return;
00122 }
00123
00124
00125 void *newdata=malloc( destsize*nel );
00126
00127 for( i=0; i<nel; i++ ){
00128 mem=a->data+(i*srcsize);
00129 array_dtype_to_double( &tmpel, mem, a->dtype );
00130 array_MEMSET( newdata+(i*destsize), target_type, tmpel );
00131 }
00132
00133 free(a->data);
00134 a->data=newdata;
00135 a->dtype=target_type;
00136 a->dtype_size=destsize;
00137 a->nbytes=destsize*nel;
00138 }
00139
00140
00167 Array* array_concatenate( const Array *a, const Array *b, int dim ){
00168 Array *out;
00169 if( a==NULL && b==NULL ){
00170 warnprintf("concatenating two NULL-arrays\n");
00171 return NULL;
00172 }
00173 if( a==NULL ){
00174 out=array_copy( b, TRUE ); return out;
00175 }
00176 if( b==NULL ){
00177 out=array_copy( a, TRUE ); return out;
00178 }
00179 if( dim!=0 && dim!=1 ){
00180 errprintf("can only concatenate rows (0) or columns (1), got %i\n",dim);
00181 return NULL;
00182 }
00183 if( a->dtype!=b->dtype ){
00184 errprintf("Arrays must be of same data-type\n");
00185 return NULL;
00186 }
00187 if( a->ndim!=b->ndim ){
00188 errprintf("Arrays must be of same dimensionality\n");
00189 return NULL;
00190 }
00191 if( a->ndim>2 || b->ndim>2 ){
00192 errprintf("Arrays must be 1D or 2D\n");
00193 return NULL;
00194 }
00195 if( a->size[1-dim]!=b->size[1-dim] ){
00196 errprintf("Arrys must be of same dimension in dim %i, have %i vs. %i\n",
00197 1-dim, a->size[1-dim], b->size[1-dim] );
00198 return NULL;
00199 }
00200
00201 Array *sa,*sb;
00202 sa = array_fromptr2( a->dtype, 2, a->data, a->size[0], (a->ndim>1)?(a->size[1]):1 );
00203 sb = array_fromptr2( b->dtype, 2, b->data, b->size[0], (b->ndim>1)?(b->size[1]):1 );
00204
00205 int outr, outc;
00206 if( dim==0 ){
00207 outr=sa->size[0]+sb->size[0];
00208 outc=sa->size[1];
00209 } else {
00210 outr=sa->size[0];
00211 outc=sa->size[1]+sb->size[1];
00212 }
00213 out = array_new2( a->dtype, 2, outr, outc );
00214
00215 if( dim==0 ){
00216 memcpy( out->data, sa->data, sa->nbytes );
00217 memcpy( out->data+sa->nbytes, sb->data, sb->nbytes );
00218 } else {
00219 int i;
00220 for( i=0; i<a->size[0]; i++ ){
00221 memcpy( array_INDEXMEM2( out, i, 0), array_INDEXMEM2( a, i, 0 ), a->dtype_size*a->size[1] );
00222 memcpy( array_INDEXMEM2( out, i, a->size[1]),
00223 array_INDEXMEM2( b, i, 0 ), b->dtype_size*b->size[1] );
00224 }
00225 }
00226
00227 array_free( sa );
00228 array_free( sb );
00229
00230 array_dimred( out );
00231 return out;
00232 }
00233
00234
00244 void array_calc_rowindex( ulong offset, const uint *size, uint nsize, uint *index ){
00245 ulong prod=1;
00246 ulong accounted=0;
00247 int i;
00248 if( nsize<=0 ){
00249 errprintf("invalid dimensionality=%i\n", nsize);
00250 return;
00251 }
00252 for( i=0; i<nsize; i++ ) prod*=size[i];
00253
00254 for( i=0; i<nsize; i++ ){
00255 index[i] = (offset-accounted)/(prod/size[i]);
00256
00257 prod /= size[i];
00258 accounted += index[i]*prod;
00259 }
00260 }
00261
00271 void array_calc_colindex( ulong offset, const uint *size, uint nsize, uint *index ){
00272 ulong prod=1;
00273 ulong accounted=0;
00274 int i;
00275 if( nsize<=0 ){
00276 errprintf("invalid dimensionality=%i\n", nsize);
00277 return;
00278 }
00279 for( i=0; i<nsize; i++ ) prod*=size[i];
00280
00281 for( i=0; i<nsize; i++ ){
00282 index[nsize-i-1] = (offset-accounted)/(prod/size[nsize-i-1]);
00283
00284 prod /= size[nsize-i-1];
00285 accounted += index[nsize-i-1]*prod;
00286 }
00287 }
00288
00297 Array* array_convert_rowcolmajor( Array *a, bool alloc){
00298 ulong i;
00299 uint *idx=(uint*)malloc( a->ndim*sizeof(uint) );
00300 Array *b = array_copy( a, TRUE );
00301
00302 for( i=0; i<array_NUMEL( a ); i++ ){
00303 array_calc_colindex( i, a->size, a->ndim, idx );
00304 memcpy( array_INDEXMEM1( b, i ), array_index( a, idx ), a->dtype_size );
00305 }
00306
00307 if( !alloc ){
00308 memcpy( a->data, b->data, b->nbytes );
00309 array_free( b );
00310 b=a;
00311 }
00312
00313 return b;
00314 }
00315
00323 void array_reverse( Array *a ){
00324 int i;
00325 void *loc1,*loc2;
00326 int n = array_NUMEL(a);
00327 void *tmp;
00328 tmp = malloc( a->dtype_size );
00329 for( i=0; i<n/2; i++ ){
00330 loc1=array_INDEXMEM1( a, i );
00331 loc2=array_INDEXMEM1( a, n-i-1 );
00332 memcpy( tmp, loc1, a->dtype_size );
00333 memcpy( loc1, loc2, a->dtype_size );
00334 memcpy( loc2, tmp, a->dtype_size );
00335 }
00336 free( tmp );
00337 }
00338
00344 bool array_comparable( const Array *a, const Array *b ){
00345 int i;
00346 if( a->dtype!=b->dtype ){
00347 char *dts1="", *dts2="";
00348 array_DTYPESTRING( dts1, a->dtype );
00349 array_DTYPESTRING( dts2, b->dtype );
00350 errprintf("arrays do not have the same datatype: '%s' vs. '%s'\n", dts1, dts2 );
00351 return FALSE;
00352 }
00353 if( a->ndim != b->ndim ){
00354 errprintf("arrays do not have the same dimensionality, %i vs. %i\n", a->ndim, b->ndim );
00355 return FALSE;
00356 }
00357 for( i=0; i<a->ndim; i++ ){
00358 if( a->size[i] != b->size[i] ){
00359 errprintf("Arrays differ in dimension %i: %i vs. %i elements\n",
00360 i, a->size[i], b->size[i] );
00361 return FALSE;
00362 }
00363 }
00364 return TRUE;
00365 }
00366
00376 Array* array_copy( const Array *in, bool allocdata ){
00377 Array *out;
00378 if( allocdata ){
00379 out = array_new( in->dtype, in->ndim, in->size );
00380 memcpy( out->data, in->data, in->nbytes );
00381 } else {
00382 out = array_fromptr( in->dtype, in->ndim, in->data, in->size );
00383 }
00384 return out;
00385 }
00386
00395 int array_dimred( Array *a ){
00396 int i,j;
00397 int nd=0;
00398 uint *nsize;
00399 for( i=0; i<a->ndim; i++ ){
00400 if( a->size[i]!=1 )
00401 nd++;
00402 }
00403 if( nd==a->ndim ) return 0;
00404
00405 MALLOC( nsize, nd, uint );
00406 j=0;
00407 for( i=0; i<a->ndim; i++ ){
00408 if( a->size[i]!=1 )
00409 nsize[j++] = a->size[i];
00410 }
00411 free( a->size );
00412 a->size=nsize;
00413 a->ndim=nd;
00414
00415 return 0;
00416 }
00417
00422 int array_scale (Array * a, double x){
00423 long i;
00424 double r;
00425 for( i=0; i<a->nbytes/a->dtype_size; i++ ){
00426 array_dtype_to_double( &r, a->data+i*a->dtype_size, a->dtype );
00427 array_MEMSET( a->data+i*a->dtype_size, a->dtype, x*r );
00428 }
00429 return 0;
00430 }
00431
00433
00434 uint fill_buffer( void *data, void *buf, int nd, uint **ind, uint *size, uint dtsize, ulong *nbelow ){
00435 int i;
00436 int writ=0;
00437 dprintf("nd=%i,d=%p,b=%p\n",nd, data, buf);
00438 if( nd==1 ){
00439 for( i=0; i<*size; i++ ){
00440 memcpy( buf+i*dtsize, data+(*ind)[i]*dtsize, dtsize );
00441 writ += dtsize;
00442 }
00443 return writ;
00444 } else {
00445 for( i=0; i<*size; i++ ){
00446 writ += fill_buffer( data+((*ind)[i])*(*nbelow)*dtsize, buf+writ, nd-1, ind+1,
00447 size+1, dtsize, nbelow+1 );
00448 }
00449 return writ;
00450 }
00451 }
00456
00457 void parse_slicedesc( const Array *a, const char *slicedesc, uint *size, uint **ind ){
00458 int i,j;
00459 char *desc, *cptr1, *cptr2, *orgdesc;
00460
00461 desc = strdup( slicedesc );
00462
00463 orgdesc=desc;
00464 for( i=0; i<a->ndim; i++ ){
00465 cptr1=strchr( desc, ',' );
00466 if( cptr1 )
00467 *cptr1='\0';
00468 if( strchr( desc,':' ) ){
00469 size[i] = a->size[i];
00470 MALLOC( ind[i], size[i], uint );
00471 for( j=0; j<size[i]; j++ )
00472 ind[i][j]=j;
00473 dprintf("Keep original dimenson\n");
00474 } else if( (cptr2=strchr( desc,'[' ))!=NULL ){
00475 MALLOC( ind[i], a->size[i], uint );
00476 j=0;
00477 size[i]=0;
00478 cptr2++;
00479 while( *cptr2 != ']' ){
00480 ind[i][j]=atoi( cptr2 );
00481 dprintf("Found index '%i'\n", ind[i][j] );
00482 size[i]++;
00483 j++;
00484 cptr2=strchr( cptr2+1, ' ' );
00485 if(!cptr2) break;
00486 }
00487 } else if( (cptr2=strchr( desc, '-' ))!=NULL ){
00488 int a1,a2;
00489 a1=atoi( desc );
00490 a2=atoi( cptr2+1 );
00491 dprintf("range %i-%i\n", a1, a2 );
00492 size[i]=(a2-a1)+1;
00493 MALLOC( ind[i], size[i], uint );
00494 for( j=0; j<(a2-a1)+1; j++ ){
00495 dprintf("j=%i\n", j);
00496 ind[i][j] = a1+j;
00497 }
00498 } else {
00499 size[i] = 1;
00500 MALLOC( ind[i], 1, uint );
00501 ind[i][0]=atoi( desc );
00502 dprintf("Found index '%i'\n", ind[i][0] );
00503 }
00504 desc=cptr1+1;
00505 }
00506 free( orgdesc );
00507 }
00522 Array* array_slice( const Array *a, const char *slicedesc ){
00523 Array *b;
00524 uint *size, **ind;
00525 int i;
00526
00527 if( strcount( slicedesc, ',' )!=a->ndim-1 ){
00528 errprintf("Slice Description does not contain enough dimensions (need %i)\n",a->ndim );
00529 return NULL;
00530 }
00531 MALLOC( size, a->ndim, uint );
00532 MALLOC( ind, a->ndim, uint*);
00533
00534
00535 parse_slicedesc( a, slicedesc, size, ind );
00536
00537
00538 long bufn=0;
00539 ulong *nbelow;
00540 MALLOC( nbelow, a->ndim, ulong );
00541 nbelow[0]=a->nbytes/a->dtype_size;
00542 for( i=1; i<a->ndim; i++ ){
00543 nbelow[i] = nbelow[i-1]/a->size[i-1];
00544 }
00545 dprintf("a->ndim=%i\n", a->ndim );
00546 b = array_new( a->dtype, a->ndim, size );
00547 bufn=fill_buffer( a->data, b->data, a->ndim, ind, size, a->dtype_size, nbelow+1 );
00548 if( bufn!=b->nbytes ){
00549 errprintf("Somethings wrong, wrote %li bytes but needed to write %li\n",
00550 bufn, b->nbytes );
00551 }
00552
00553
00554 array_dimred( b );
00555
00556
00557 free( nbelow );
00558 free( size );
00559 for( i=0; i<a->ndim; i++ ){
00560 free( ind[i] );
00561 }
00562 free( ind );
00563
00564 return b;
00565 }
00566
00568
00569 void dump_data( FILE *out, DType dt, uint dtsize, uint nd, uint nel, void *data,
00570 uint *size, ulong *nbelow ){
00571 int i, len=nel;
00572 if( nel<=0 || nel>=*size )
00573 len=*size;
00574 if( nd==1 ){
00575 fprintf( out, "[ ");
00576 for( i=0; i<len; i++ ){
00577 array_DTYPEPRINT( out, dt, data+i*dtsize );
00578 fprintf( out, ", " );
00579 }
00580 if( len<*size )
00581 fprintf( out, "... " );
00582 fprintf( out, "]\n" );
00583 } else {
00584 fprintf( out, "[ " );
00585 for( i=0; i<len; i++ ){
00586 dump_data( out, dt, dtsize, nd-1, nel, data+i*(*nbelow)*dtsize, size+1, nbelow+1 );
00587 }
00588 if( len<*size )
00589 fprintf( out, "...\n" );
00590 fprintf( out, "]\n" );
00591 }
00592 }
00603 void array_print( Array *a, uint nel_per_dim, FILE *out ){
00604 int i;
00605 ulong *nbelow;
00606 MALLOC( nbelow, a->ndim, ulong );
00607 nbelow[0]=a->nbytes/a->dtype_size;
00608 for( i=1; i<a->ndim; i++ ){
00609 nbelow[i] = nbelow[i-1]/a->size[i-1];
00610 }
00611 char *dt=NULL;
00612 array_DTYPESTRING(dt, a->dtype);
00613 fprintf( out, "array(%s, ", dt );
00614 for( i=0; i<a->ndim; i++ ){
00615 fprintf( out, "%i, ", a->size[i] );
00616 }
00617 fprintf( out, "\b\b):\n" );
00618 dump_data( out, a->dtype, a->dtype_size, a->ndim, nel_per_dim, a->data, a->size, nbelow+1 );
00619 fprintf( out, "\n");
00620
00621 free( nbelow );
00622 }
00623
00638 void* array_index( const Array *a, uint *idx ){
00639 int i;
00640 uint index=0;
00641 uint cumdim;
00642
00643 if( a->nbytes % a->dtype_size!=0 ){
00644 errprintf("Corrupt array, nbytes not divisible by dtype_size\n");
00645 }
00646 cumdim = a->nbytes/a->dtype_size;
00647 for( i=0; i<a->ndim; i++ ){
00648
00649
00650
00651
00652
00653 cumdim /= a->size[i];
00654 index += idx[i]*cumdim*a->dtype_size;
00655 }
00656 return a->data+index;
00657 }
00658
00672 void* array_index2( const Array *a, ... ){
00673 va_list ap;
00674 uint *idx;
00675 int i;
00676 void *r;
00677 MALLOC( idx, a->ndim, uint );
00678
00679 va_start (ap, a );
00680 for( i=0; i<a->ndim; i++ ){
00681 idx[i] = (uint)va_arg( ap, uint );
00682 }
00683 va_end (ap);
00684 r = array_index( a, idx );
00685 free( idx );
00686
00687 return r;
00688 }
00689
00699 Array *array_new ( DType dtype, uint ndim, const uint *dims ){
00700 Array *a;
00701 int i;
00702 MALLOC( a, 1, Array );
00703 a->dtype=dtype;
00704 array_SIZEOF_DTYPE( a->dtype_size, dtype );
00705 a->ndim=ndim;
00706 a->free_data=TRUE;
00707 MALLOC( a->size, ndim, uint );
00708
00709
00710 a->nbytes=1;
00711 for( i=0; i<ndim; i++ ){
00712 a->size[i] = dims[i];
00713 a->nbytes *= a->size[i];
00714 }
00715 a->nbytes *= a->dtype_size;
00716
00717
00718 a->data = (void*)malloc( a->nbytes );
00719 if( !a->data ){
00720 errprintf("Failed to allocate memory!\n");
00721 } else {
00722 memset( a->data, 0, a->nbytes );
00723 }
00724
00725 return a;
00726 }
00727
00737 Array *array_new2( DType dtype, uint ndim, ... ){
00738 Array *a;
00739 va_list ap;
00740 uint *size;
00741 int i;
00742 MALLOC( size, ndim, uint );
00743
00744
00745 va_start (ap, ndim );
00746 for( i=0; i<ndim; i++ ){
00747 size[i] = (uint)va_arg( ap, uint );
00748 }
00749 va_end (ap);
00750
00751 a = array_new( dtype, ndim, size );
00752
00753 free( size );
00754 return a;
00755 }
00756
00767 Array *array_randunif( unsigned long seed, uint ndim, ... ){
00768 va_list ap;
00769 Array *a;
00770 long i;
00771 int n;
00772 uint *size;
00773 MALLOC( size, ndim, uint );
00774
00775
00776 va_start (ap, ndim );
00777 for( i=0; i<ndim; i++ ){
00778 size[i] = (uint)va_arg( ap, uint );
00779 }
00780 va_end (ap);
00781
00782 const gsl_rng_type * T;
00783 gsl_rng * r;
00784 gsl_rng_env_setup();
00785 T = gsl_rng_default;
00786 r = gsl_rng_alloc (T);
00787 if( seed==0 )
00788 seed=(unsigned long)time(NULL);
00789 gsl_rng_set (r, seed );
00790 a = array_new( DOUBLE, ndim, size );
00791 n = a->nbytes/a->dtype_size;
00792 for( i=0; i<n; i++ ){
00793 array_INDEX1( a, double, i )= gsl_rng_uniform (r);
00794 }
00795
00796 gsl_rng_free (r);
00797 free( size );
00798 return a;
00799 }
00800
00805 void array_shuffle( Array *a, unsigned long seed ){
00806 const gsl_rng_type * T;
00807 gsl_rng * r;
00808 gsl_rng_env_setup();
00809 T = gsl_rng_default;
00810 r = gsl_rng_alloc (T);
00811 if( seed==0 )
00812 seed=(unsigned long)time(NULL);
00813 gsl_rng_set (r, seed );
00814
00815 gsl_ran_shuffle( r, a->data, array_NUMEL(a), a->dtype_size );
00816
00817 gsl_rng_free (r);
00818 }
00819
00820
00830 Array *array_new_dummy( DType dtype, uint ndim, ... ){
00831 va_list ap;
00832 Array *a;
00833 long i;
00834 int n;
00835 uint *size;
00836 MALLOC( size, ndim, uint );
00837
00838
00839 va_start (ap, ndim );
00840 for( i=0; i<ndim; i++ ){
00841 size[i] = (uint)va_arg( ap, uint );
00842 }
00843 va_end (ap);
00844
00845 a = array_new( dtype, ndim, size );
00846 n = a->nbytes/a->dtype_size;
00847 for( i=0; i<n; i++ ){
00848 array_MEMSET( a->data+i*a->dtype_size, a->dtype, i );
00849 }
00850
00851 free( size );
00852 return a;
00853 }
00854
00869 Array *array_fromptr2( DType dtype, uint ndim, void *data, ... ){
00870 va_list ap;
00871 Array *a;
00872 int i;
00873 uint *size;
00874 MALLOC( size, ndim, uint );
00875
00876
00877 va_start (ap, data );
00878 for( i=0; i<ndim; i++ ){
00879 size[i] = (uint)va_arg( ap, uint );
00880 }
00881 va_end (ap);
00882
00883 a = array_fromptr( dtype, ndim, data, size );
00884
00885 free( size );
00886
00887 return a;
00888 }
00889
00904 Array *array_fromptr( DType dtype, uint ndim, void *data, const uint *size ){
00905 Array *a;
00906 int i;
00907 MALLOC( a, 1, Array );
00908 a->dtype=dtype;
00909 array_SIZEOF_DTYPE( a->dtype_size, dtype );
00910 a->ndim=ndim;
00911 a->free_data=FALSE;
00912 MALLOC( a->size, ndim, uint );
00913
00914
00915 a->nbytes=1;
00916 for( i=0; i<ndim; i++ ){
00917 a->size[i] = size[i];
00918 a->nbytes *= a->size[i];
00919 }
00920 a->nbytes *= a->dtype_size;
00921
00922 a->data = data;
00923
00924 return a;
00925 }
00926
00933 void array_dtype_to_double( double *out, void *mem, DType dt ){
00934 switch( dt ){
00935 case CHAR:
00936 *out = *((char*)mem);break;
00937 case UINT:
00938 *out = *((uint*)mem);break;
00939 case INT:
00940 *out = *((int*)mem);break;
00941 case LONG:
00942 *out = *((long*)mem);break;
00943 case ULONG:
00944 *out = *((ulong*)mem);break;
00945 case FLOAT:
00946 *out = *((float*)mem);break;
00947 case DOUBLE:
00948 *out = *((double*)mem);break;
00949 default:
00950 warnprintf("Do not know this datatype: %i\n", dt );
00951 }
00952 }
00953
00954
00955
00958 void array_free( Array *a ){
00959 if( !a ) return;
00960 if( a->free_data && a->data ) free( a->data );
00961 if( a->size ) free( a->size );
00962 free( a );
00963 }
