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 <<EOF;
Usage for some drsutils.pm functions:

  makeFileExtExists : splits (or adds) file name into parts 

  isNumeric : checks if input is a number  
  isInt : checks if input is an integer
  isPosInt : checks if input is a positve integer
  isNegInt : checks if input is a negative integer
  isNonNegInt : checks if input is 0 or a positive integer

  stripComment : strips input line of comments
  isComment : checks if input line is wholly a comment

  arrayEq : checks if two arrays of numbers are identical
  arrayEqSTR : checks if two arrays of strings are identical
  indexFirstOf : searches for occurrence of number in array 
  indexFirstOfSTR : searches for occurrence of string in array 
  diffArray : finds differences of successive elements of array of numbers
  copyArrayRef : makes copy of array given a reference to it

  numdp : finds minimum number of decimal places to represent input number
  num2string : creates a formatted string to represent a number
  numarray2string : returns string of formatted input array of numbers
  array2string : returns string with elements input array

  array2FH : outputs input array to file handle
  array2stdout : outputs input array to screen
  printArray : outputs array with many options    
  printHash : prints Hash of elements to screen
  printHashOfArrays : prints Hash of arrays of elements to screen

  nearZero : checks if input is close to zero  

  makeLog : returns array of real numbers from 0 to 1 that increase logarithmically
  makeExp : returns array of real numbers from 0 to 1 that increase exponentially

  mag : returns magnitude of input vector
  normalize : returns normalized version of input vector
  gramschmidt : finds orthonormal basis of n input d-dim vectors (n<=d)
EOF
}



sub makeFileExtExists{
  # used to split a file into its prefix and its extension e.g. blah.jpg -> (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) <fnameS> EXISTS (i.e. is a file in the current directory)
  #      (b) <fnameS>.<ext1S> or ... or <fnameS>.<extnS> EXISTS
  #      (c) Neither (a) nor (b)
  #   If (a) or (b) 
  #     returns (okfname,xS,eS, eS was one of extnS ,"") where okfname = <xS>.<eS> EXISTS
  #     If (a) 
  #       If <fnameS> has a period then <fnameS> equals <xS>.<eS> where <xS> may have periods but <eS> does not.
  #       If <fnameS> 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 <fnameS>.<extnS> 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<scalar(@{$aref}); $i++)
   { return $i if ($aref->[$i] == $x);}
  return -1;
}

sub indexFirstOfSTR{
  # as for indexFirstOf but for arrays of strings

  my ($x,$aref) = @_;
  for (my $i=0; $i<scalar(@{$aref}); $i++)
   { return $i if ($aref->[$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<j} <ej,vi>*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;