################################################################################
#
#  kPerl Logics/Mathematics Library Group
#  Object Library for N-dimensional Kartesian Point or Point Vector
#
#  copyright (c)
#    Karol Szafranski, 2005
#  author
#    Karol Szafranski, karol.szafranski@leibniz-fli.de
#
################################################################################
#
#  DESCRIPTION
#
# - purpose and background
#   Object-based computation with N-dimensional Kartesian point or point vector
#
# - switches may be set in two ways:
#   - library-wide, i.e. all active objects are affected. Here, the way to
#     assign switches is just an assignment to the library switch hash, like:
#       $Math::GeometNdPoint::LibGlob{switch}{-debug} = 1;
#   - object-associated: assign switches via method AddSwitch().
#   Note that lirary-wide switch assignments take effect in the moment of
#   object construction. Later changes to the library properties will not
#   take effect on existing objects.
#
# - individual description of functions at the beginning of the code blocks
#
################################################################################
#
#  OBJECT OPERATORS
#
# - data type conversion
#   bool           always TRUE
#   @{}            [ $x, $y, ... ]
#
# - calculation
#   neg            multiply with -1
#   +=  +          vector addition
#   -=  -          vector addition
#   *=  *          stretch vector
#   /=  /          stretch vector
#
#
#  OBJECT METHODS  for external access
#
# - housekeeping
#   new            return object, possible initializations (via ini()):
#                  - with 1 numerical arg, yielding zero-vector of chosen
#                    dimension
#                  - with x args: yielding explicit vector ($arg1,$arg2,...)
#                  - with 1 arg (object), yielding copy of object
#   ini            object initialization, args of new()
#   AddSwitch      assign object switches
#   clone          return copy of object
#   copy           return copy of object, alias to clone()
#   string         represent object in compact string format, matrix values in
#                  nested parentheses
#
# - attributes and arithmetics
#   dim            return/set dimension of point vector
#   coord          return/set coordinate array
#   amount         determine amount of point vector
#   norm_direct    return normalized direction vector
#   stretch        stretch vector, i.e. multiply with float
#   add            add another vector
#   move           move point along vector, alias to add()
#
# - arithmetics with a group of points/vectors
#   points2vectors turn a list of points to a list of connecting vectors
#                  exported as function &geometnd_pts2vects
#  (dependent      determine dependency for a set of vectors
#                  exported as function &geometnd_dependent
#   scalar_prod    return scalar product of two vectors
#                  exported as function &geometnd_scalarprod
#   strip_vector   remove the vector content along a given axis
#  (cross_prod     return cross product of a set of vectors
#                  exported as function &geometnd_crossprod
#
# - interface to Math::GeometNdTrform
#   transform      transform vector/point by matrix
#
#
#  OBJECT DATA STRUCTURE
#  (hash)
#
#   dim            dimension
#   coo            reference to coordinate array
#
################################################################################
#
#  FUNCTIONS, DATA
#
#   @EXPORT
#   @ISA
#   overload
#   %LibGlob
#
# - housekeeping
#   $LibGlob{switch}
#   &new  see MainLib::DefaultObjHash.pm
#   &ini  called by &new
#   &AddSwitch
#   &clone
#   &copy
#   &string
#
# - attributes and arithmetics
#   &dim
#   &coord
#   &amount
#   &stretch
#   &add
#   &move
#
# - arithmetics with a group of points/vectors
#   &points2vectors
#   &dependent
#   &scalar_prod
#   &strip_vector
#   &cross_prod
#
# - interface to Math::GeometNdTrform
#   &transform
#
#
#  STD OPTIONS
#
#  -debug      print debug protocol to STDERR
#
################################################################################
#
#  DEBUG, CHANGES, ADDITIONS
#
#  - There is no linear formula to compute the determinant beyond 3D space.
#
################################################################################

package Math::GeometNdPoint;

# includes
#use strict; use warnings;  # OK 20050623
use MainLib::DefaultObjHash;
use Math::Analysis;
### beware: ###
### use Math::kCalc qw(&stretch);
  # stretch() is defined as an object-specific method here

# symbol export
use Exporter qw(import);
our @EXPORT = qw (
  &geometnd_pts2vects &geometnd_dependent
    &geometnd_crossprod &geometnd_scalarprod
  );

# inheritance
our @ISA = qw(MainLib::DefaultObjHash);

# overload
use overload (
'""'   => # make the stringify operator work like for non-overloaded object references
          sub { \%{$_[0]} },
'bool' => sub { 1 },
'@{}'  => sub { $_[0]->coord() },
'neg'  => sub { $_[0]->clone()->stretch(-1) },
'+='   => sub { $_[0]->add($_[1]) },
'+'    => sub { $_[0]->clone()->add($_[1]) },
'-='   => sub { $_[0]->add(-$_[1]) },
'-'    => sub { $_[0]->clone()->add(-$_[1]) },
'*='   => sub { $_[0]->stretch($_[1]) },
'*'    => sub { $_[0]->clone()->stretch($_[1]) },
'/='   => sub { $_[0]->stretch(1/$_[1]) },
'/'    => sub { $_[0]->clone()->stretch(1/$_[1]) },
);

# package-wide constants and variables
our (%LibGlob);


################################################################################
# housekeeping
################################################################################


$LibGlob{switch} = {};


# parametric initialization
#
# INTERFACE
# - argument 0:  object reference, from method new()
# - argument 1:  dimension
# - argument 2*: initialization argument(s), either:
#                - existing object (function as a copy constructor)
#                - array of Kartesian coordinates
#                - reference to array of Kartesian coordinates
#                - none => point for Kartesian origin
#
# - library options:
#   -debug      [STD]
#
# - return val: object reference
#
# DESCRIPTION
# - this function will be called by method new(), as a constructor sub-function.
#
sub ini {
  my ($this,$dim,@arg) = @_;

  # switches
  $this->{switch} = { %{$LibGlob{switch}} };
  my $debug = $this->{switch}{-debug};

  # Kartesian dimensions
  if (($this->{dim}=$dim) < 2) {
    die sprintf "%s. ERROR: bad dimension: %d\n", (caller(0))[3], $dim;
  }

  # construct from submitted data argument(s)
  if (@arg) {

    # convert refs to array, likewise for object of this same class
    if (int(@arg) == 1) {
      if (ref($arg[0]) eq 'ARRAY' and int(@{$arg[0]})==$dim) {
        return &ini ($this, $dim, @{$arg[0]});
      } elsif (ref($arg[0]) eq __PACKAGE__) {
        $debug and printf STDERR "%s. called as a copy constructor\n", (caller(0))[3];
        $this = $arg[0]->clone();
      } else {
        die sprintf "%s. ERROR: unable to interpret arguments: %s\n", (caller(0))[3],
          join(',',@arg)||"''";
      }
    }

    # assign from array of coordinate values
    elsif (int(@arg) == $dim) {
      $this->{coo} = [ @arg ];
      $debug and printf STDERR "%s. assigning from arguments -> %s\n", (caller(0))[3], $this->string();
    }

    # arg number error
    else {
      die sprintf "%s. ERROR: too many arguments: %d\n", (caller(0))[3], int(@arg);
    }
  }

  # assign null point
  else {
    $this->{coo} = [ (0) x $dim ];
    $debug and printf STDERR "%s. assigning default: %s\n", (caller(0))[3], $this->string();
  }

  return $this;
}


# enter object switches
#
# INTERFACE
# - return val: success status (boolean)
#
sub AddSwitch {
  my ($this, %oopt) = @_;

  # loop over switches
  while (my($key,$val) = each(%oopt)) {
    if (0) { }

    # we don't enter -TabType into object switches
    elsif ($key eq '-TabType') { next }

    # options that we just have to enter
    else {
      if (defined $val) { $this->{switch}{$key} = $val; }
      else       { delete $this->{switch}{$key}; }
    }
  }

  # return success status
  return 1;
}


# copy object
#
# INTERFACE
# - return val: reference to copy of object
#
# DEVELOPER'S NOTES
# - perform limited recursion copy, since overload data records may contain
#   large-tree or circular references
#
sub clone {
  my ($this) = @_;
  my $pCopy = { %$this };
  $pCopy->{coo} = [ @{$this->{coo}} ];
  bless ($pCopy, ref($this));
    # this way of blessing also works for derived classes
}

sub copy { &clone(@_) }


# compact string representation of object
#
# INTERFACE
# - options:
#   -fmt        full sprintf format string. Note: the sprintf expression has
#               as much numerical arguments as dimensions. Use this with care!
#   -fmtfloat   sprintf format substring for float conversion, default %s
# - return val: string representation
#
sub string {
  my ($this,%opt) = @_;
  my $fmtfloat = $opt{-fmtfloat} || '%s';
  my $lbl = $opt{-label} || 'point';
  my $fmt = $opt{-fmt} || sprintf("$lbl(%s)",join(',',($fmtfloat)x$this->{dim}));
  eval 'sprintf($fmt,$this->coord())' || '* fmt error *';
}


################################################################################
# attributes and arithmetics
################################################################################


# return/set dimension of point vector
#
# INTERFACE
# - argument 1*: dimension (assignment mode)
# - return val:  dimension
#
sub dim {
  my ($this,$arg) = @_;
  if (defined ($arg)) {
    if ($arg < 2) {
      die sprintf "%s. ERROR: bad argument: %s\n", (caller(0))[3], $arg||"''";
    }
    $this->{dim} = int($arg);
    $#{$this->{coo}} = int($arg)-1;
  }
  return $this->{dim};
}


# return/set coordinate array
#
# INTERFACE
# - argument 1*: coordinate array (assignment mode)
# - return val:  coordinate array
#
sub coord {
  my ($this,@arg) = @_;
  if (@arg) {
    if (int(@arg) < 2) {
      die sprintf "%s. ERROR: bad argument number %d\n", (caller(0))[3], int(@arg);
    }
    $this->{dim} = int(@arg);
    $this->{coo} = [ @arg ];
  }
  return @{$this->{coo}};
}


# vector's amount
#
# INTERFACE
# - return val: amount
#
sub amount {
  my ($this) = @_;
  my $sum=0;
  for (my $i=0; $i<$this->dim(); ++$i) {
    $sum += $this->{coo}[$i]**2;
  }
  return sqrt($sum);
}


# return normalized direction vector
#
# INTERFACE
# - argument 0: object reference
# - options:
#   -amount     provide pre-calculated vector amount
# - return val: vector with normalized length
#
# DESCRIPTION
# - This method works as a constructor.
# - A null vector evaluates to a null vector.
#
sub norm_direct {
  my ($this,%opt) = @_;
  my $amt = $opt{-amount} || $this->amount();
  if ($amt) {
    return $this->copy()->stretch(1/$amt);
  } else {
    return $this->new($this->dim());
  }
}


# stretch vector (multiply with floating point), overload on operator "*="
#
# INTERFACE
# - argument 1: multiplier
# - return val: passed object reference
#
# DEVELOPER'S NOTES
# - multiplication makes sense only for vectors, actually not for points. This
#   method will move to derived class "GeometNdVector" if it's once created.
#
sub stretch {
  my ($this,$multiplier) = @_;
  $multiplier||=0;
  if (ref ($multiplier)) {
    die sprintf "%s. ERROR: unable to interpret argument: %s\n", (caller(0))[3], $multiplier||"''";
  }

  # apply multiplier
  my $dim = $this->dim();
  for (my $i=0; $i<$dim; ++$i) {
    $this->{coo}[$i] *= $multiplier;
  }

  return $this;
}


# add another vector, overload on operator "+="
#
# INTERFACE
# - argument 1: reference to vector
# - return val: passed object reference
#
# DEVELOPER'S NOTES
# - the rather long sequence to obtain the addend vector is to allow addend
#   operands other than this same class, but process operands of this same class
#   as fast as possible.
#
sub add {
  my ($this,@ArgVec) = @_;
  my ($pAddend);
  if (ref($pAddend=$ArgVec[0]) ne __PACKAGE__) {
    unless ($pAddend = eval("\$this->new(\@ArgVec)")) {
      die sprintf "%s. ERROR: unable to interpret arguments: %s\n", (caller(0))[3],
        join (',', @ArgVec);
    }
  }

  # recalculate coordinates
  my @addcoo = $pAddend->coord();
  my $dim = $this->dim();
  for (my $i=0; $i<$dim; ++$i) {
    $this->{coo}[$i] += $addcoo[$i];
  }

  return $this;
}

# move point by vector, alias to add()
sub move { &add(@_) }


################################################################################
# arithmetics with a group of points/vectors
################################################################################


# turn a list of points to a list of connecting vectors
#
# INTERFACE
# - argument 1+: list of points
# - return val:  list of vectors (n-1)
#
sub points2vectors {
  my ($this,@ArgVec) = @_;
  if (int(@ArgVec) < 1) {
    die sprintf "%s. ERROR: bad number of arguments\n", (caller(0))[3], int(@ArgVec);
  }

  my @ConnectVec;
  while (my $pVec = shift @ArgVec) {
    push @ConnectVec, $pVec-$this;
    $this = $pVec;
  }
  return @ConnectVec;
}

sub geometnd_pts2vects { &points2vectors(@_) }


# determine dependency for a set of vectors
#
# INTERFACE
# - argument 1:  reference to another vector
# - argument 2*: reference to more vectors
# - return val:  dependency (boolean)
#
# DEBUG, CHANGES, ADDITIONS
# - There is no linear formula to compute the determinant beyond 3D space.
#
sub dependent {
  my ($this,@ArgVec) = @_;
  my $debug = $this->{switch}{-debug};
  if (grep{ ref($_) ne __PACKAGE__ } @ArgVec) {
    die sprintf "%s. ERROR: unable to interpret one of arguments\n", (caller(0))[3];
  }
  # we could allow mixed dimensionality and convert all vectors to maximum
  # dimensionality
  #require Math::kCalc; Math::kCalc->import('&Max');
  #my $dim = &Max(map{ $_->dim() } $this,@ArgVec);
  if (grep{ $_->dim() != $this->dim() } @ArgVec) {
    warn sprintf "%s. WARNING: varying dimensions in arguments\n", (caller(0))[3];
  }

  # construct matrix for determinant analysis
  my @WorkVec;
  foreach my $pVec ($this,@ArgVec) {
    # for mixed dimensionality, clone and redimension here
    #my $pVecCp = $pVec->clone();
    #$pVecCp->dim($dim);
    push @WorkVec, [ $pVec->coord() ];
  }
  if ($debug) {
    require MainLib::Data;
    printf STDERR "%s. matrix given to &determinant\n", (caller(0))[3];
    &MainLib::Data::DataPrint(\@WorkVec,-handle=>\*STDERR,-space=>2);
  }

  return &determinant(\@WorkVec,-debug=>$debug) < 1e-5;
}

sub geometnd_dependent { &dependent(@_) }


# return scalar product of two vectors
#
# INTERFACE
# - argument 1: reference to another vector
# - return val: scalar product
#
sub scalar_prod {
  my ($this,$pVec) = @_;
  if (ref($pVec) ne __PACKAGE__) {
    die sprintf "%s. ERROR: unable to interpret argument %s\n", (caller(0))[3], $pVec;
  }
  # find maximum dimensionality
  if ($pVec->dim() != $this->dim()) {
    if ($pVec->dim() > $this->dim()) {
      ($this=$this->clone())->dim($pVec->dim());
    } {
      ($pVec=$pVec->clone())->dim($this->dim());
    }
  }

  # calculate sum of products
  my $sum=0;
  for (my $i=0; $i<$this->dim(); ++$i) {
    $sum += $this->{coo}[$i] * $pVec->{coo}[$i];
  }

  return $sum;
}

sub geometnd_scalarprod { &scalar_prod(@_) }


# remove the vector content along a given axis
#
# INTERFACE
# - argument 1: reference to another vector
# - return val: passed object reference
#
sub strip_vector {
  my ($this,$pVec) = @_;
  if (ref($pVec) ne __PACKAGE__) {
    die sprintf "%s. ERROR: unable to interpret argument %s\n", (caller(0))[3], $pVec;
  }
  my $aa=$this->amount();
  my $an=$this->norm_direct(-amount=>$aa);
  my $bn=$pVec->norm_direct();

  my $co = $an->scalar_prod($bn);
  $bn->stretch(-$co*$aa);
  return $this->add($bn);
}


# return cross product of a set of vectors
#
# INTERFACE
# - argument 1:  reference to another vector
# - argument 2*: reference to more vectors
# - return val:  cross product (orthogonal vector)
#
# DESCRIPTION
# - The number of arguments directs the dimensionality of the calculation. For
#   n vectors (including the base object), the result will be (n+1)-dimensional.
#   Accordingly, the dimensionality of the vector arguments will be set to n.
#
# DEBUG, CHANGES, ADDITIONS
# - There is no linear formula to compute the determinant beyond 3D space.
#
sub cross_prod {
  my ($this,@ArgVec) = @_;
  my $debug = $this->{switch}{-debug};
  if (grep{ ref($_) ne __PACKAGE__ } @ArgVec) {
    die sprintf "%s. ERROR: unable to interpret one of arguments\n", (caller(0))[3];
  }
  if (grep{ $_->dim() != $this->dim() } @ArgVec) {
    warn sprintf "%s. WARNING: varying dimensions in arguments\n", (caller(0))[3];
  }

  # construct matrix for determinant analysis
  my $dim = int(@_) + 1;
  my @WorkVec;
  foreach my $pVec ($this,@ArgVec) {
    my $pVecCp = $pVec->clone();
    $pVecCp->dim($dim);
    push @WorkVec, [ $pVecCp->coord() ];
  }

  # construct vector component in matrix for determinant analysis
  my $pVecComp = [ 1, (0) x ($dim-1) ];
  my @VecComp;
  for (my $i=0; $i<$dim; ++$i) {
    push @VecComp, $this->new($dim,$pVecComp);
    unshift @$pVecComp, pop @$pVecComp;
  }
  unshift @WorkVec, [ @VecComp ];
  if ($debug) {
    require MainLib::Data;
    printf STDERR "%s. matrix given to &determinant\n", (caller(0))[3];
    &MainLib::Data::DataPrint(\@WorkVec,-handle=>\*STDERR,-space=>2);
  }

  return &determinant(\@WorkVec,-init=>$this->new($dim),-debug=>$debug);
}

sub geometnd_crossprod { &cross_prod(@_) }


################################################################################
# interface to Math::GeometNdTrform
################################################################################


# transform by matrix object
#
# INTERFACE
# - argument 1: reference to matrix object
# - return val: transformed point (passed object reference)
#
# DEBUG, CHANGES, ADDITIONS
# - adjust the transformation matrix (copy) to the requested dimension?
#   But, how to tell what is the requested dimensionality?
#
sub transform {
  my ($this,$pMat) = @_;
  if (ref($pMat) ne 'Math::GeometNdTrform') {
    die sprintf "%s. ERROR: unable to interpret argument: %s\n", (caller(0))[3], $pMat||"''";
  }
  if ($pMat->dim() != $this->dim()) {
    warn sprintf "%s. WARNING: dimensions differ between point/transformation\n", (caller(0))[3];
  }

  # apply matrix
  my $dim = $this->dim();
  my $pPt0 = $this->clone(); $pPt0->{coo}[$dim] = 1;
  for (my $t=0; $t<$dim; ++$t) {
    my $cell=0;
    for (my $s=0; $s<=$dim; ++$s) {
      $cell += $pMat->{mat}[$s][$t]*$pPt0->{coo}[$s];
    }
    $this->{coo}[$t] = $cell;
  }

  return $this;
}


1;
# $Id: GeometNdPoint.pm,v 1.13 2018/06/05 18:02:56 szafrans Exp $
