package drsutils; use POSIX qw(ceil floor abs); # Dinoj Surendran dinoj@cs.uchicago.edu / dinojs@gmail.com # Jan 2007 use constant PI => 4 * atan2(1, 1); # following are taken from http://hea-www.harvard.edu/~alexey/calc-src.txt # by Alexey Vikhlinin sub ln {log(\$_[0]);} sub lg {log(\$_[0])/\$log10;} sub r2d {\$_[0]*180.0/PI} sub d2r {\$_[0]*PI/180.0} sub atan {atan2(\$_[0],1)} sub tan {my \$z = \$_[0]; sin(\$z)/cos(\$z)} sub acos {my \$z = \$_[0]; atan2(sqrt(1-\$z*\$z), \$z)} #sub acos {my \$z = \$_[0]; my \$zsq = \$z*\$z; atan2((\$zsq <= 1 ? sqrt(1-\$z*\$z) : 0), \$z)} #sub asin {my \$z = \$_[0]; my \$zsq = \$z*\$z; atan2(\$z, (\$zsq <= 1 ? sqrt(1-\$z*\$z) : 0))} sub asin {my \$z = \$_[0]; atan2(\$z, sqrt(1-\$z*\$z))} sub sind {sin(\$_[0]*PI/180)} sub cosd {cos(\$_[0]*PI/180)} sub tand {tan(\$_[0]*PI/180)} sub asind {asin(\$_[0])*180/PI} sub acosd {acos(\$_[0])*180/PI} sub atand {atan(\$_[0])*180/PI} # end of Alexey's functions sub help { print STDERR < (blah,jpg) # Notational Notes: anything ending with S is a string variable # the value of a variable v is [v] # input (fnameS,[ext1S,...,extNS]) ... best to provide extnS in CAPITALS # fnameS can include directory stuff e.g. "../../yo/bar/image.jpg" # Checks if either # (a) EXISTS (i.e. is a file in the current directory) # (b) . or ... or . EXISTS # (c) Neither (a) nor (b) # If (a) or (b) # returns (okfname,xS,eS, eS was one of extnS ,"") where okfname = . EXISTS # If (a) # If has a period then equals . where may have periods but does not. # If has no period then fnameS == xS and eS = "" # If (b) then xS == fnameS and eS = extnS where n is the minimum n(=1..N) such that . EXISTS # If (c) # returns ("","","",0,\$errmsg); my \$f=shift(@_); my @exts=(); foreach \$x(@_) { push(@exts,\$x); if (lc(\$x) ne \$x){ push(@exts,lc(\$x)); } } my \$okfname=""; my \$pre=""; my \$post=""; # filename \$f is [\$pre].[\$post] where \$post is one of @exts (upper or lower case) my \$postInExts=0; my \$errmsg="cannot find any of the following files:\n"; if (-e \$f){ \$okfname = \$f; if (\$f =~ m/(.*)\.(\S+)/){ \$pre = \$1; \$post = \$2; foreach \$ext(@exts){ \$postInExts=1 if (\$ext eq \$post); } }else{ \$pre = \$f; } }else{ \$errmsg = "\$errmsg \$f\n"; foreach \$ext(@exts){ if (-e (\$f.".\$ext")){ \$pre = \$f; \$post = \$ext; \$postInExts=1; \$okfname = \$f.".\$ext"; }else{ \$errmsg.= " \$f".".\$ext\n"; } } } if (! length(\$pre)){ return ("","","",0,\$errmsg); }else{ return (\$okfname,\$pre,\$post,\$postInExts,""); } } sub isNumeric{ # usage 1 : # input : \$x # output : true iff \$x is a number # # usage 2: # input : @X = (\$x0,\$x1,\$x2,...) # output : 0 if @X has at least one non-numeric entry that isnt "" # ("" is considered neither numeric nor non-numeric) # or n if @X has n numeric entries and no non-numeric entries (other than "") # # e.g. ("",3,"",0) or (3,0) results in 2 being returned # e.g. ("",3,"a",0) results in 0 being returned my \$numNumbers = 0; my \$x; if (\$#_ > 0){ my @tmp=@_; foreach \$x(@tmp){ if (\$x ne ""){ if (isNumeric(\$x)){ \$numNumbers++; }else{ return 0; } } } }elsif (\$#_ == 0){ \$x = \$_[0]; if (\$x =~ m/^[-+]?\d+\$/) {\$numNumbers = 1;} elsif (\$x =~ m/^[-+]?\d*\.\d+\$/) {\$numNumbers = 1;} elsif (\$x =~ m/(^[\d+-\.]+)[e|E]([\d-+\.]+)\$/){ \$numNumbers = (((\$1 !~ m/eE/) && (\$2 !~ m/eE/)) && (isNumeric(\$1) && isNumeric(\$2))); } } return \$numNumbers; } sub arrayEq{ # input is (n1,n2,\@a1,\@a2) where n1 = numelements in @a1, ditto n2,@a2 # returns true iff @a1 has same elements as @a2 my (\$n1,\$n2,\$aref1,\$aref2) = @_; return 0 if (\$n1 != \$n2); for (my \$i=0; \$i<\$n1; \$i++){ return 0 if (\$aref1->[\$i] != \$aref2->[\$i]); } return 1; } sub arrayEqSTR{ # as for arrayEq but for arrays of strings my (\$n1,\$n2,\$aref1,\$aref2) = @_; return 0 if (\$n1 != \$n2); for (my \$i=0; \$i<\$n1; \$i++){ return 0 if (\$aref1->[\$i] ne \$aref2->[\$i]); } return 1; } sub isInt{ # returns true iff input is an integer # return ((isNumeric(\$_[0])) && ((POSIX::floor(\$_[0]) * 1.0 == \$_[0]))); # print "isInt : checking <",\$_[0],"> = ",\$_[0]-POSIX::floor(\$_[0]),"\n"; return (isNumeric(\$_[0]) && nearZero(\$_[0]-POSIX::floor(\$_[0]),1e-16)); } sub isNonNegInt{ # returns true iff input is a positive number or zero return ((isInt(\$_[0])) && (\$_[0] >= 0)); } sub isNegInt{# returns true iff input is a negative number return ((isInt(\$_[0])) && (\$_[0] < 0)); } sub isPosInt{ # returns true iff input is a positive number return ((isInt(\$_[0])) && (\$_[0] > 0)); } sub stripComment{ # strips lines of comments # # input : (\$line,@commentchars) # if any element of @commentchars is found at i-th position of \$line # then line[0..i-1] is returned # else line is returned my (\$line,@cmt)=@_; foreach \$c(@cmt){ if ((\$line =~ m/^(.*)[\$c].*\$/)){ return \$1; } } return \$line; } sub lineIsComment{ #input : (\$x,@commentchars) # returns 1 iff if \$x starts with an element of @commentchars (after possibly some whitespace) my (\$x,@cmt)=@_; foreach \$c(@cmt){ return 1 if ((\$_[0] =~ m/^\s*\$c.*/)); } } sub searchArrayIgnoreCase{ # input : (\$x,@a) for some numeric or char array @a # output: -1 if @a doesnt contain \$x # i>=0 if \$a[\$i] == \$x for the first possible \$i for \$i(1..\$#_){ return \$i-1 if (lc(\$_[\$i]) eq lc(\$_[0])); } return -1; } sub searchArray{ if (isNumeric(\$_[0])){ return searchArrayNumeric(@_); }else{ return searchArraySTR(@_); } } sub searchArraySTR{ # input : (\$x,@a) for some numeric or char array @a # output: -1 if @a doesnt contain \$x # i>=0 if \$a[\$i] == \$x for the first possible \$i for \$i (1 .. \$#_){ return \$i-1 if (\$_[\$i] eq \$_[0]); } return -1; } sub searchArrayNumeric{ for \$i (1 .. \$#_){ return \$i-1 if (\$_[\$i] == \$_[0]); } return -1; } sub indexFirstOf { # input : (\$x,\@a) for some numeric array @a # output: -1 if @a doesnt contain \$x # i>=0 if \$a[\$i] == \$x for the first possible \$i my (\$x,\$aref) = @_; for (my \$i=0; \$i[\$i] == \$x);} return -1; } sub indexFirstOfSTR{ # as for indexFirstOf but for arrays of strings my (\$x,\$aref) = @_; for (my \$i=0; \$i[\$i] eq \$x);} return -1; } sub copyArrayRef{ # input : \@b for some array @b # output: \@c for new array @c that is copy of @b my (\$aref) = @_; my @a = @{\$aref}; my @b = @a; return \@b; } sub numdp{ # finds number of decimal places required to represent a number # input is (x,[dpmax] # dpmax (15 by default) is the value returned by numdp if x requires at least maxdp decimal places. # output is the minimum number (or dpmax) of decimal places needed to represent x my \$n = (\$_[0] >= 0 ? \$_[0] : -\$_[0]); my \$dpmax = (\$#_ >= 1 ? \$_[1] : 15); my \$eps = 0.1 ** (\$dpmax-1); my \$dp = 0; my \$remnant = POSIX::abs(\$n - POSIX::floor(.5+\$n)); while (! nearZero(\$remnant, \$eps) && \$dp <= \$dpmax){ \$remnant *= 10; \$remnant = POSIX::abs(\$remnant - POSIX::floor(.5+\$remnant)); \$dp++; \$eps *= 10; } return \$dp; } sub array2string{ join(" ",@_); } sub numarray2string{ # prints formatted numeric array # input : (k,@a) # output : string with elements of @a on a single line. # If k==-1 then the numbers in @a are printed with as few decimal places as possible # If k==0 then the numbers in @a are printed as integers # If k>0 then the numbers in @a are printed with k decimal places # # e.g. (-1,0.134,13, 6,-41.3) -> "[ 0.134 13 6 -41.3 ]" # e.g. (0 ,0.134,13, 6,-41.3) -> "[ 0 13 6 -41 ]" # e.g. (2 ,0.134,13, 6,-41.3) -> "[ 0.13 13.00 6.00 -41.30 ]" my \$numdec = shift(@_); my \$s = ""; foreach \$x(@_) {\$s = sprintf "%s %s ",\$s,num2string(\$x,\$numdec);} \$s; } sub num2string{ # converts number to string with user-specified (or as few as possible) number of decimal places # input : (n,[dp]) # dp is taken to be numdp(n) by default or if dp is <0 # output : n with dp decimal places (dp is 0 for integer) return \$_[0] if (! isNumeric(\$_[0])); my \$dp = (\$#_>=1 && \$_[1] >= 0 ? \$_[1] : numdp(\$_[0])); my \$a = (\$dp == 0 ? sprintf( "%d " , \$_[0]) : sprintf( "%.\$dp"."f " , \$_[0])); return \$a; } sub diffArray{ # compute difference of array # input is array @d with N elements # output is array @v with N-1 elements and v[t]=d[t+1]-d[t] my @v = (); for (my \$i=1; \$i<=\$#_; \$i++){ \$v[\$i-1] = \$_[\$i] - \$_[\$i-1]; } return @v; } sub array2FH { # prints array to file handle # input is (\$OFH,\$numdec,@a) # array @a is printed to output file handle OFH printArray(\$#_-1,\$@[2..\$#],1,\$_[1],\$_[0]); } sub array2stdout { # prints array to file handle # input is (\$numdec,@a) # array @a is printed to output file handle OFH # numdec should be -1 if you want to use as few decimal places as possible printArray(\$#_,@_[1..\$#_],1,\$_[0]); } sub printArray{ # prints array to stdout or file handle # input is (n,@a,[printNewLine],[numdecplaces],[outputFileHandle]) # where n = #elements in @a # printNewLine is 1 if a final newline should be output and 0 (default) otherwise # outputFileHandle should be a valid output file handle if specified. If not specified, output is to stdout my \$n = \$_[0]; my \$newline = ( (\$#_ >= \$n+1) ? \$_[\$n+1] : 0); my \$numdecplaces = ( (\$#_ >= \$n+2) ? \$_[\$n+2] : -1); my \$t; if (\$#_ >= \$n+3){ my \$ofh = \$_[\$n+3]; print \$ofh "[ "; foreach \$t(@_[1..\$n]) { printf \$ofh "%s ", (isNumeric(\$t) ? num2string(\$t,\$numdecplaces) : \$t); } print \$ofh "]"; print \$ofh "\n" if (\$newline); }else{ print "[ "; foreach \$t(@_[1..\$n]) { printf "%s ", (isNumeric(\$t) ? num2string(\$t,\$numdecplaces) : \$t); } print "]"; print "\n" if (\$newline); } } sub printHash{ my (\$href,\$hname) = @_; %h = %{\$href}; my @k = keys %h; printf "\$hname has ", \$#k+1, " entries\n"; foreach \$x (@k){ print ("\$hname [",\$x,"] = [",\$h{\$x},"]\n"); } } sub printHashOfArrays{ my (\$href,\$hname) = @_; %h = %{\$href}; my @k = keys %h; print ("\$hname has ", \$#k+1, " entries\n"); foreach \$x (@k){ print ("\$hname [",\$x,"] = "); printArray ( @{\$h{\$x}}); } } sub nearZero{ # returns true if input is near zero # input : (\$x,[\$eps]) # (default value of \$eps is 10^-6) # returns true iff -\$eps <= \$x <= \$eps my \$eps = (\$#_ >= 1 ? \$_[1] : 1e-6); return ((\$_[0] > (-1*\$eps)) && (\$_[0]<\$eps)); } sub makeLog{ # creates array of numbers increasing from 0 to 1 # input (T,a) where # T = final timestep (first timestep is indexed 0) # a = acceleration i.e. (frac[t]-frac[t-1])/(frac[t-1]-frac[t-2]) # (by default a=1 i.e. constant acceleration) # output : array @frac with \$T+1 elements that increase from 0 to 1 # \$frac[0] = 0, # \$frac[t] = fraction of total distance covered by at time t, 0<=t<=T # (frac[t] < frac[t+1]) # \$frac[T] = 1 # # e.g. makeLog (4,1) or makeLog (4) returns (0,.25,.5,.75,1) # makeLog (5,1.5) returns (0, 0.032, 0.097, 0.226, 0.484, 1) my \$T = shift(@_); my \$a = (\$#_ >= 0 ? \$_[0] : 1.0); die if (\$T<=0); die if (\$a<=0); die if (! isInt(\$T)); my \$t, @frac = (0); if (\$a == 1){ @frac = map {\$_/\$T} (0..\$T); }else{ \$frac[1] = (1-\$a)/(1-(\$a ** \$T)); foreach \$t (2..\$T){ \$frac[\$t] = \$frac[\$t-1] + \$a*(\$frac[\$t-1]-\$frac[\$t-2]); } @frac = map {\$_/\$frac[\$T]} @frac; # to be sure } return @frac; } sub makeExp{ # creates array of numbers increasing from 0 to 1 # input (T,[v0]) where # T = final timestep (first timestep is indexed 0) # v0 = optional initial value for frac[1] # = initial velocity as fraction of total distance # (v0 only used if <1/T and T>2) # # output : array @frac with \$T+1 elements that increase from 0 to 1 # \$frac[0] = 0, # \$frac[t] = fraction of total distance covered by at time t, 0<=t<=T # (frac[t] < frac[t+1]) # \$frac[T] = 1 # my \$T = shift(@_); my @frac=(); if (\$#_>=0 && \$_[0]*\$T<1 && \$T > 2){ # \$v0 = \$_[0]; my @tmp = map{(1-(\$_[0]*(\$T-1)))*\$_} makeExp(\$T-1); @frac = (0, \$_[0], map{\$_[0]*\$_ + \$tmp[\$_]} (1 .. \$T-1)); }else{ @frac = (0,map {exp(\$_)/exp(\$T)} (1..\$T)); } return @frac; } sub mag { # returns magnitude of input vector my \$sumsq=0; map {\$sumsq+=(\$_**2) } @_; sqrt(\$sumsq); } sub normalize { # input is numeric array representing vector, # output is normalized vector my \$m = mag(@_); map { \$_/(\$m ? \$m : 1) } @_ ; } sub gramschmidt{ # does Gram-Schmidt orthonomalization of any n d-dimensional vectors that # form a n-dimensional basis (n<=d) # input : d @v0, @v1, ..., @vn # output : @e0, @e1, ..., @en : orthonormal basis my(\$d) = shift(@_); \$n = (@_ / \$d) - 1; my @e=(); for (my \$i=0; \$i<=\$n; \$i++){ my @inds = (\$i*\$d .. \$i*\$d + \$d-1); splice(@_,\$i*\$d,\$d,normalize(@_[@inds])); my @sum=(0,0,0); for (my \$j=0; \$j<\$i; \$j++){ # ei = normalize( vi - sum_{i*ej ) my @jnds = (\$j*\$d .. \$j*\$d+\$d-1); my \$dotprod=0; for \$k(0..\$#jnds) {\$dotprod += (\$e[\$jnds[\$k]] * \$_[\$inds[\$k]]);} for \$k(0..\$#sum) {\$sum[\$k] += \$dotprod*\$e[\$jnds[\$k]];} } for \$k(0..2) { \$e[\$i*\$d+\$k] = \$_[\$i*\$d+\$k] - \$sum[\$k];} splice(@e,\$i*\$d,\$d,normalize(@e[@inds])); } return @e; } 1;