################################################################################
#
#  kPerl Logics/Mathematics Library Group
#  Object Library for 2D Data Plot
#
#  copyright (c)
#    Institute of Molecular Biotechnology Jena, Dept. Genome Analysis, 1999-2001, 2004
#  author
#    Karol Szafranski, karol.szafranski@leibniz-fli.de
#
################################################################################
#
#  DESCRIPTION
#
# - purpose
#   provide an OO representation of plot data, i.e. 2D value vectors.
#
# - individual description of functions at the beginning of the code blocks
#
################################################################################
#
#  OBJECT OPERATORS
#
#   bool           validity status
#   @{}            2D table data
#   -=             subtract scalar or another plot, method Subtract()
#
#
#  OBJECT METHODS  for external access
#
# - housekeeping
#   new            return object, initialize via ini()
#                  optionally, load plot data
#   ini            initialize object
#
# - data I/O
#   Clear          initialize object, keep switch settings
#   Data           $arg: load data from supplied table data reference
#                  no $arg: return plot data (formatted as table data)
#
# - evaluate data
#   Size           pixel number in plot
#   Xmin           minimum x value
#   Xmax           maximum x value
#   Ymin           minimum y value
#   Ymax           maximum y value
#   Ymean          mean of y values
#   Bound          return x boundaries of defined data pixels that frame an
#                  x value
#   Extrapolate    extrapolate for given x value that's out of defined x range
#   Interpolate    do linear interpolation of y for given x value
#                  if x gets out of defined x range Extrapolate() is called
#                  automatically
#   Integral       calculate integral of x value interval
#
# - recalculate data
#   FlipWing       flip wing of plot and add to opposite wing
#                  used for smoothening with fixed boundary values
#   SmoothVal      calculate smoothened y value
#   SmoothPlot     calculate smoothened plot
#   Subtract       recalculate by subtracting scalar or another plot
#
#
#  OBJECT DATA STRUCTURE
#
#   switch         hash reference for object switches, cf. ini()
#     -debug         print debug protocol to STDERR
#     -extrapolate   method of extrapolation:
#                    border  default: assign y of border x to all x in undefined
#                            range
#                    mean    assign mean y to all x in undefined range
#                    mirror  not implemented! mirror plot at end of defined
#                            range
#                    val=Y   assign defined y value (arg statment) to undefined
#                            x range
#     -SmoothProfile profile for smoothening
#     -SmoothWindow  window size for smoothening
#   pdata          reference to plot data. It's stored as table data structure,
#                  type AA.
#   SecData        reference to secondary data. There may be:
#     Ymin         minimum of y values
#     Ymax         maximum of y values
#     Ymean        mean of y values
#
################################################################################
#
#  FUNCTIONS, DATA
#
#   @ISA
#   %_LibParam
#
# - housekeeping
#   description of object parameters
#   $_LibParam{default}
#   &new  see MainLib::DefaultObjHash.pm
#   &ini
#   &AddSwitch
#   &_LocalSwitch  see MainLib::DefaultObjHash.pm
#
# - data I/O
#   &Clear
#   &_LoadData
#   &_LoadFile
#   &_ReturnData
#   &Data
#
# - evaluate plot
#   &Size
#   &Xmin
#   &Xmax
#   &Ymin
#   &Ymax
#   &Ymean
#   &Bound
#   &Interpolate
#   &Extrapolate
#   &Integral
#
# - recalculate plot
#   &FlipWing
#   description of smoothening profile function
#   %SmoothProfile
#   &SmoothVal
#   &SmoothPlot
#   &Subtract
#
################################################################################

package Math::Plot2D;

# includes
#use strict; use warnings;  # OK 20040507
use MainLib::DefaultObjHash;
use MainLib::Data qw(&DataClone &DataPrint);
use MainLib::Misc qw(&MySub);
use Math::kCalc;
use Math::Sort;
use database::Table;
use database::DbPlain;

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

# overload
use overload (
bool  => sub { 1; },
#'0+'  => "Numify",
'@{}' => "_ReturnData",
'-='  => "Subtract",
#'x'   => sub { print "%s::x. arguments:\n", __PACKAGE__;
#                map { printf "  %s\n", $_ } @_ },
);

# package-wide constants and variables
my %_LibParam;


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


# default object data structure
#
$_LibParam{default} = {
  TabType => 'AA',
  switch  => {
    -debug         => undef,
    -extrapolate   => 'border',
    -SmoothProfile => 'rectangle',
    -SmoothWindow  => 5.0,
    }
  };


# parametric initialization
#
# INTERFACE
# - argument 1: reference to data or name of table source file
#
# - options:    options may be used in conjunction with a data argument
#   -TabType    table data structure type if different from standard type
#
# - return val: - object reference
#               - undef if an error occurs
#
# DESCRIPTION
# - this function will be called by method new()
#
sub ini {
  my ($this,$arg,%opt) = @_;
  my $debug = $opt{-debug};

  # initialize data structure
  # Clear() is done only in conjunction Load...() functions
  $this->{switch} = $_LibParam{default}{switch};

  # enter data
  if (ref($arg) and ref($arg) !~ m/\b(FileHandle|GLOB)/) {
    $debug and printf STDERR "%s. entering data %s\n", &MySub, $arg;
    $this->_LoadData($arg,%opt) or return undef;
  } elsif ((!ref($arg) and (-e $arg or $arg eq '-'))
    or ref($arg) =~ m/\b(FileHandle|GLOB)/) {
    $debug and printf STDERR "%s. loading data from file %s\n", &MySub, $arg;
    $this->_LoadFile($arg,%opt) or return undef;
  } else {
    $debug and printf STDERR "%s. ERROR: unable to interpret data argument %s\n", &MySub, $arg||"''";
    return undef;
  }
  delete $opt{-TabType};
  $debug and printf STDERR "%s. %d data lines\n", &MySub, int(@{$this->{pdata}});

  # enter object switches
  $this->AddSwitch(%opt) or return undef;

  # return object
  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;
}


################################################################################
# data I/O
################################################################################


# clear table data
#
# INTERFACE
# - argument 1: reference to table data
#
sub Clear {
  my ($this) = @_;

  # drop data structures
  delete $this->{pdata};
  delete $this->{SecData};
}


# enter table data
#
# INTERFACE
# - argument 1: reference to table data
#
# - local options, overriding object options:
#   -debug      [STD]
#   -TabType    table data structure type if different from standard type
#               ($_LibParam{default}{TabType})
#
# - return val: success status (boolean)
#
sub _LoadData {
  my ($this,$pTable,%opt) = @_;
  my %lopt = $this->_LocalSwitch(%opt);
  my $debug = $lopt{-debug};
  my $dbg2  = $debug ? $debug-1 : undef;

  # clear any existing table data
  $this->Clear();

  # convert table data if needed
  if ($lopt{-TabType} and $lopt{-TabType} ne $_LibParam{default}{TabType}) {
    unless ($pTable = &TableConvert ($lopt{-TabType}, $_LibParam{default}{TabType}, $pTable, -debug=>$dbg2)) {
      $debug and printf STDERR "%s. ERROR: unable to convert table data type (%s -> %s)\n", &MySub,
        $lopt{-TabType}||"''", $_LibParam{default}{TabType}||"''";
      return undef;
    }
  }

  # enter sorted table data into object structure
  $this->{pdata} = [ sort { $$a[0]<=>$$b[0] } @$pTable ];

  return 1;
}


# load table data from file
#
# INTERFACE
# - argument 1: file path, process call, or file handle reference
#
# - local options, overriding object options:
#   ...         all options affecting &PlainToTable
#   -debug      [STD]
#
# - return val: success status (boolean)
#
sub _LoadFile {
  my ($this,$ArgTable,%opt) = @_;
  my %lopt = $this->_LocalSwitch(%opt);
  my $debug = $lopt{-debug};
  my $dbg2  = $debug ? $debug-1 : undef;

  # clear table data
  $this->Clear();

  # load table
  unless ($this->{pdata} = &PlainToTable ($ArgTable, -TabType=>$_LibParam{default}{TabType}, %opt, -debug=>$dbg2)) {
    $debug and printf STDERR "%s. ERROR: unable to load table data from file/handle %s\n", &MySub, $ArgTable;
    return undef;
  }

  # success
  return 1;
}


# return table data
#
# INTERFACE
# - local options, overriding object options:
#   -debug      [STD]
#   -TabType    table data structure type to be returned
#
# - return val: success status (boolean)
#
sub _ReturnData {
  my ($this,%opt) = @_;
  my %lopt = $this->_LocalSwitch(%opt);
  my $debug = $lopt{-debug};
  my $dbg2  = $debug ? $debug-1 : undef;

  # grab data from data anchor in object
  my $pTable = $this->{pdata};

  # convert table data?
  if ($lopt{TabType} and $lopt{-TabType} ne $_LibParam{default}{TabType}) {
    unless ($pTable = &TableConvert ($lopt{-TabType}, $_LibParam{default}{TabType}, $pTable, -debug=>$dbg2)) {
      $debug and printf STDERR "%s. ERROR: unable to convert table data type (%s -> %s)\n", &MySub,
        $_LibParam{default}{TabType}||"''", $lopt{-TabType}||"''";
      return undef;
    }
  }
  
  # return data
  return $pTable;
}


# object interface to table data
#
# INTERFACE
# - argument 1: (optional) reference to table data. If this argument is
#               given the function switches to the data enter mode (means:
#               _LoadData()).
#
# - options:
#   everything is handed over to either _LoadData() or _ReturnData()
#
# - return val: that of _LoadData() or _ReturnData()
#
sub Data {
  my ($this, @ArgPlus) = @_;
  my ($arg, %opt);

  # splice true argument from argument list
  while (defined($_=$ArgPlus[0]) and !m/^-/) {
    $arg ||= shift @ArgPlus;
  }
  %opt = @ArgPlus;
  
  # enter data
  # clear all existing data
  if ($arg) {
    return $this->_LoadData($arg,%opt);
  }
  
  # return data
  else {
    return $this->_ReturnData(%opt);
  }
}


################################################################################
# evaluate data
################################################################################


# total plot pixel number
#
# INTERFACE
# - return val: plot pixel number
#
sub Size {
  my ($this) = @_;
  return int (@{$this->{pdata}});
}


# left border of x value range
#
# INTERFACE
# - return val: - left border of x value range
#               - undef if an error occurred
#
sub Xmin {
  my ($this) = @_;
  int (@{$this->{pdata}}) or return undef;
  return $this->{pdata}[0][0];
}


# right border of x value range
#
# INTERFACE
# - return val: - right border of x value range
#               - undef if an error occurred
#
sub Xmax {
  my ($this) = @_;
  int (@{$this->{pdata}}) or return undef;
  return $this->{pdata}[$#{$this->{pdata}}][0];
}


# minimum of y values
#
# INTERFACE
# - return val: - minimum of y values
#               - undef if an error occurred
#
sub Ymin {
  my ($this,%opt) = @_;

  # get and store maximum of y values
  unless (defined $this->{SecData}{Ymin}) {
    int (@{$this->{pdata}}) or return undef;
    $this->{SecData}{Ymin} = &Min (map { $_->[1] } @{$this->{pdata}});
  }

  return $this->{SecData}{Ymin};
}


# maximum of y values
#
# INTERFACE
# - return val: - maximum of y values
#               - undef if an error occurred
#
sub Ymax {
  my ($this,%opt) = @_;

  # get and store maximum of y values
  unless (defined $this->{SecData}{Ymax}) {
    int (@{$this->{pdata}}) or return undef;
    $this->{SecData}{Ymax} = &Max (map { $_->[1] } @{$this->{pdata}});
  }

  return $this->{SecData}{Ymax};
}


# mean of y values
#
# INTERFACE
# - return val: - mean of y values
#               - undef if an error occurred
#
sub Ymean {
  my ($this,%opt) = @_;
  # get and store meanimum of y values
  $this->{SecData}{Ymean} ||= &Mean (map { $_->[1] } @{$this->{pdata}});

  return $this->{SecData}{Ymean};
}


# return defined x boundaries framing an x value
#
# INTERFACE
# - argument 1: x value
#
# - options:
#   -debug      this may change debug behaviour for this particular function
#               call
#
# - return val: - data structure reference
#                 [ [$BoundLeftNum,$BoundLeftX,$BoundLeftY], [$BoundRightNum,$BoundRightX,$BoundRightY] ]
#               - undef if an error occurred
#
sub Bound {
  my ($this,$ValX,%opt) = @_;
  my %lopt = $this->_LocalSwitch(%opt);
  my $debug = $lopt{-debug};
  my $dbg2  = $debug ? $debug-1 : undef;

  # look for plot data interval
  # - &ListLocateSorted yields right neighbor of $ValX
  # - if &ListLocateSorted finds an element that equals $ValX, that element will become the
  #   right neighbor (see sorting code)
  my $iBound = &ListLocateSorted ($this->{pdata}, $ValX,
    -sort  => sub { my $XNew=$_[0]; my $XList=${$_[1]}[0]; $XNew<=>$XList or -1; },
    -debug => $dbg2);
  if (not defined($iBound)) {
    die sprintf "%s. ERROR in ListLocateSorted() for x value %f\n", &MySub, $ValX;
  }

  # outside data range?
  if ($iBound == int(@{$this->{pdata}})) {
    return [ [ $#{$this->{pdata}}, @{$this->{pdata}[-1]} ], undef ];
  }
  if ($iBound==0 and $this->{pdata}[0][0]!=$ValX) {
    return [ undef, [ 0, @{$this->{pdata}[0]} ] ];
  }

  # exact match?
  if ($this->{pdata}[$iBound][0] == $ValX) {
    my @b = ($iBound, @{$this->{pdata}[$iBound]});
    return [ [ @b ], [ @b ] ];
  }

  # two independent boundaries
  return [ [ $iBound-1, @{$this->{pdata}[$iBound-1]} ], [ $iBound, @{$this->{pdata}[$iBound]} ] ];
}


# interpolate given x value to y value
#
# INTERFACE
# - argument 1: x value to be interpolated to y value
#
# - options:
#   -debug      this may change debug behaviour for this particular function
#               call
#
# - return val: - interpolated y value
#               - undef if an error occurred
#
# DESCRIPTION
# - on data range escape, $this->Extrapolate is called automatically
#
sub Interpolate {
  my ($this,$ValX,%opt) = @_;
  my %lopt = $this->_LocalSwitch(%opt);
  my $debug = $lopt{-debug};

  # look for plot data interval
  my $pBound = $this->Bound ($ValX);
  $debug and printf STDERR "%s. boundaries: [ %s, %s ]\n", &MySub,
    $$pBound[0] ? sprintf ("[%d, %f, %f]", @{$$pBound[0]}) : 'undef',
    $$pBound[1] ? sprintf ("[%d, %f, %f]", @{$$pBound[1]}) : 'undef';

  # outside of data range? => chain to &Extrapolate
  my $ValInterpol;
  if (grep { ! defined $_ } @$pBound) {
    $ValInterpol = $this->Extrapolate ($ValX, %opt);
    $debug and printf STDERR "%s. got value %f from \&Extrapolate\n", &MySub, $ValInterpol;
    return $ValInterpol;
  }

  # exact hit on plot pixel?
  if ($$pBound[0][0] == $$pBound[1][0]) { return $this->{pdata}[$$pBound[0][0]][1]; }

  # safety check
  unless ($this->{pdata}[$$pBound[1][0]][0] - $this->{pdata}[$$pBound[0][0]][0]) { return undef }

  # interpolate linearly
  $ValInterpol = $this->{pdata}[$$pBound[0][0]][1] +
    ($ValX                                - $this->{pdata}[$$pBound[0][0]][0]) *
    ($this->{pdata}[$$pBound[1][0]][1] - $this->{pdata}[$$pBound[0][0]][1]) /
    ($this->{pdata}[$$pBound[1][0]][0] - $this->{pdata}[$$pBound[0][0]][0]);

  # exit SUBs
  return $ValInterpol;
}


# extrapolate given x value to y value
#
# INTERFACE
# - argument 1: x value to be extrapolated to y value
#
# - options:
#   -debug      this may change debug behaviour for this particular function
#               call
#
# - return val: - extrapolated y value
#               - undef if an error occurred
#
sub Extrapolate {
  my ($this,$ValX,%opt) = @_;
  my %lopt = $this->_LocalSwitch(%opt);
  my $debug = $lopt{-debug};

  # extrapolation method 'border'
  if ($lopt{-extrapolate} eq 'border') {
    if ($ValX < $this->Xmin()) {
      return $this->{pdata}[0][1];
    }
    if ($ValX > $this->Xmax()) {
      return $this->{pdata}[$#{$this->{pdata}}][1];
    }
  }

  # extrapolation method 'mean'
  if ($lopt{-extrapolate} eq 'mean') {
    return $this->Ymean();
  }

  # extrapolation method 'val'
  if ($lopt{-extrapolate} =~ m/^val=(.+)/) {
    return $1;
  }

  # unknown extrapolation method
  return undef;
}


# calculate integral for x value interval
#
# INTERFACE
# - argument 1: left x interval border, default: Xmin
# - argument 2: right x interval border, default: Xmax
# - return val: - integral
#               - undef if an error occurred
#
sub Integral {
  my ($this,$BordLeft,$BordRight,%opt) = @_;
  my %lopt = $this->_LocalSwitch(%opt);
  my $debug = $lopt{-debug};
  my $dbg2  = $debug ? $debug-1 : undef;

  # calculate interval borders
  my %border;
  $border{left}{x}  = defined($BordLeft) ?  $BordLeft  : $this->{pdata}[0][0];
  $border{right}{x} = defined($BordRight) ? $BordRight : $this->{pdata}[$#{$this->{pdata}}][0];
  $border{left}{y}  = $this->Interpolate($border{left}{x},-debug=>$dbg2);
  $border{right}{y} = $this->Interpolate($border{right}{x},-debug=>$dbg2);
  $debug and printf STDERR "%s. interval: %.3f (y=%.3f) to %.3f (y=%.3f)\n", &MySub,
    (map{($_->{x},$_->{y})} $border{left}, $border{right});
  if ($border{right}{x} < $border{left}{x}) {
    die sprintf "%s. ERROR: inconsistent interval borders, x values [%s,%s]", &MySub,
      map{$border{$_}{x}} qw(left right);
  }

  # find leftmost (lower) intervening border
  $border{middle}{num} = &ListLocateSorted ($this->{pdata}, $border{left}{x},
    -sort=>sub {my $EntryNew=$_[0]; my $EntryList=${$_[1]}[0]; $EntryNew<=>$EntryList;}, -debug=>$dbg2);
  if ($border{middle}{num} > $#{$this->{pdata}}) {
    %{$border{middle}} = %{$border{right}};
  } else {
    $border{middle}{x} = $this->{pdata}[$border{middle}{num}][0];
    $border{middle}{y} = $this->{pdata}[$border{middle}{num}][1];
  }
  $debug and printf STDERR "%s. first (lowest) intervening border: %.3f (y=%.3f)\n", &MySub,
    $border{middle}{x}, $border{middle}{y};

  # go through first and intervening partial intervals
  my $integral;
  while ($border{middle}{x} < $border{right}{x}) {
    $debug and printf STDERR "%s. interval border at head of iteration: (%.3f,%.3f)\n", &MySub,
      $border{middle}{x}, $border{middle}{y};
    $integral += ($border{middle}{x} - $border{left}{x}) *
      &Mean ($border{middle}{y}, $border{left}{y});
    %{$border{left}} = %{$border{middle}};
    unless ($border{middle}{num} < $#{$this->{pdata}}) { last }
    $border{middle}{num} ++;
    $border{middle}{x} = $this->{pdata}[$border{middle}{num}][0];
    $border{middle}{y} = $this->{pdata}[$border{middle}{num}][1];
  }

  # last partial interval
  $integral += ($border{right}{x} - $border{left}{x}) *
    &Mean ($border{right}{y}, $border{left}{y});

  # exit SUB
  return $integral;
}


################################################################################
# recalculate
################################################################################


# flip wing of plot and add to opposite wing
#
# INTERFACE
# - argument 1: x value
# - argument 2: wing to flip, -1:=left, 1:=right
#
# - return val: - smoothened data value
#               - undef if an error occurred
#
sub FlipWing {
  my ($this,$ValX,$wing,%opt) = @_;
  my %lopt = $this->_LocalSwitch(%opt);
  my $debug = $lopt{-debug};

  # add outer wing to plot corpus
  my ($i,@v);
  my $minorx=0.0000000001;
  if ($wing<0) {
    for ($i=0; $i<int(@{$this->{pdata}}) and $this->{pdata}[$i][0]<$ValX; ++$i) {
      my $x=$ValX+$ValX-$this->{pdata}[$i][0];
      push @v, [$x, $this->{pdata}[$i][1]+$this->Interpolate($x)]
    }
    if ($this->{pdata}[$i][0]==$ValX) { ++$i }
    push @v, [$ValX+$minorx, $this->Interpolate($ValX-$minorx)+$this->Interpolate($ValX+$minorx)];
    for (; $i<int(@{$this->{pdata}}) and $this->{pdata}[$i][0]<$ValX+$ValX-$this->Xmin(); ++$i) {
      my $x=$ValX+$ValX-$this->{pdata}[$i][0];
      push @v, [$this->{pdata}[$i][0], $this->{pdata}[$i][1]+$this->Interpolate($x)]
    }
    @v = sort{ $a->[0]<=>$b->[0] } @v;
    for (my $j=0; $j<int(@v)-1; ++$j) {
      if ($v[$j][0]==$v[$j+1][0]) {
        $v[$j][1]=($v[$j][1]+$v[$j+1][1])/2;
        splice @v, $j+1, 1;
      }
    }
    splice @{$this->{pdata}}, 0, $i, @v;
  }
  else {
    for ($i=$#{$this->{pdata}}; $i>=0 and $this->{pdata}[$i][0]>$ValX; --$i) {
      my $x=$ValX+$ValX-$this->{pdata}[$i][0];
      push @v, [$x, $this->{pdata}[$i][1]+$this->Interpolate($x)]
    }
    if ($this->{pdata}[$i][0]==$ValX) { --$i }
    push @v, [$ValX+$minorx, $this->Interpolate($ValX-$minorx)+$this->Interpolate($ValX+$minorx)];
    for (; $i>=0 and $this->{pdata}[$i][0]>$ValX+$ValX-$this->Xmax(); --$i) {
      my $x=$ValX+$ValX-$this->{pdata}[$i][0];
      push @v, [$this->{pdata}[$i][0], $this->{pdata}[$i][1]+$this->Interpolate($x)]
    }
    @v = sort{ $a->[0]<=>$b->[0] } @v;
    for (my $j=0; $j<int(@v)-1; ++$j) {
      if ($v[$j][0]==$v[$j+1][0]) {
        $v[$j][1]=($v[$j][1]+$v[$j+1][1])/2;
        splice @v, $j+1, 1;
      }
    }
    splice @{$this->{pdata}}, $i+1, int(@{$this->{pdata}})-$i-1, @v;
  }

  return $this;
}


# description of smoothening profile function
#
# INTERFACE
# - argument 1: object reference
# - argument 2: x value
# - argument 3: window size
#               SmoothVal() ensures that this is > 0
# - return val: - smoothened data value
#               - undef if an error occurred
# 


my %SmoothProfile;

# smoothening profile 'rectangle'
#
$SmoothProfile{rectangle} = sub {
  my ($this,$ValX,$WinSize,%opt) = @_;
  my %lopt = $this->_LocalSwitch(%opt);
  my $debug = $lopt{-debug};
  my $dbg2  = $debug ? $debug-1 : undef;
  $WinSize or return undef;

  # calculate mean via integral
  $debug and printf STDERR "$%s::SmoothProfile\{'rectangle'\}. Smoothening at $ValX, window size $WinSize\n", __PACKAGE__;
  my $SmoothVal = $this->Integral ($ValX-0.5*$WinSize, $ValX+0.5*$WinSize, -debug=>$dbg2)
    / $WinSize;

  # exit SUB
  return $SmoothVal;
};


# calculate smoothened value
#
# INTERFACE
# - argument 1: x value
#
# - function options (beside the object options): 
#   -profile    type of smoothening profile to be used. Choose one of the
#               built-in profiles (see above) or provide a code reference for a
#               profile function (argument behaviour as described above).
#   -window     profile window size
#
# - return val: - smoothened data value
#               - undef if an error occurred
#
sub SmoothVal {
  my ($this,$ValX,%opt) = @_;
  my %lopt = $this->_LocalSwitch(%opt);
  my $debug = $lopt{-debug};
  my $profile = $lopt{-profile} || $lopt{-SmoothProfile};
  my $WinSize = $lopt{-window} || $lopt{-SmoothWindow} or return undef;

  # chain to profile function
  if (ref $profile eq 'CODE') {
    return &$profile ($this, $ValX, $WinSize);
  } elsif (ref($SmoothProfile{$profile}) eq 'CODE') {
    return &{$SmoothProfile{$profile}} ($this, $ValX, $WinSize);
  } else {
    $debug and printf STDERR "%s. ERRONEOUS profile argument %s\n", &MySub, $profile||"''";
    return undef;
  }

}


# calculate smoothened data plot
#
# INTERFACE
# - argument 1: re-plot step size
#
# - function options (beside the object options): 
#   -IntervOff  x offset of smooth re-plot, default: $this->Xmax()-$WinSize*0.95
#   -IntervEnd  x end of smooth re-plot, default: $this->Xmax()+$WinSize*0.95
#   -profile    type of smoothening profile to be used. Choose one of the
#               built-in profiles (see above) or provide a code reference for a
#               profile function (argument behaviour as described above).
#   -TabType    table data structure type of output if different from standard
#               type ($_LibParam{default}{TabType})
#   -window     profile window size, default: re-plot step size
#
# - return val: reference to smoothened data plot (table data structure,
#               table type "AA")
#
sub SmoothPlot {
  my ($this,$StepSize,%opt) = @_;
  my %lopt = $this->_LocalSwitch(%opt);
  my $debug = $lopt{-debug};
  $StepSize or die sprintf "%s. ERROR: essentially need a step size!\n", &MySub;
  my $WinSize = $lopt{-window} || $lopt{-SmoothWindow} || $StepSize;

  # re-plot offset, end and step size
  my $IntervEnd = defined($lopt{-IntervEnd}) ? $lopt{-IntervEnd} : $this->Xmax()+$WinSize*0.95;
  my $IntervOff = defined($lopt{-IntervOff}) ? $lopt{-IntervOff} : $this->Xmin()-$WinSize*0.95;
  unless ($IntervOff < $IntervEnd) {
    die sprintf "%s. ERROR: bad interval parameter: %s..%s\n", &MySub, $IntervOff, $IntervEnd
  }
  # ensure that the last step neatly fits to the rescanning interval
  # Or, could it be a problem to slightly modify the step size?
  $StepSize *= (int(($IntervEnd-$IntervOff)/$StepSize) + 0.95) / int(($IntervEnd-$IntervOff)/$StepSize);
  $debug and printf STDERR "%s. x value range: [%s,%s]\n", &MySub, $this->Xmin(), $this->Xmax();
  $debug and printf STDERR "%s. x interval parameters: [%s,%s], step %s, window size %s\n", &MySub,
    $IntervOff, $IntervEnd, $StepSize, $WinSize;

  # profile function
  my $profile = $lopt{-profile} || $lopt{-SmoothProfile};
  my $pAction;
  if (ref($profile) eq 'CODE') {
    $pAction = $profile;
  } elsif (ref($SmoothProfile{$profile}) eq 'CODE') {
    $pAction = $SmoothProfile{$profile};
  } else { die sprintf "%s. ERROR: bad profile argument %s\n", &MySub, $profile||"''" }

  # smoothening loop
  my @plot;
  for (my $CtVal=$IntervOff; $CtVal<=$IntervEnd; $CtVal+=$StepSize) {
    push @plot, [ $CtVal, &$pAction($this,$CtVal,$WinSize) ];
  }

  return \@plot;
}


# subtract operand
#
# INTERFACE
# - argument 1: subtraction operand
# - options:
#   -copy       return new object, keep object untouched
# - return val: - integral
#               - undef if an error occurred
#
sub Subtract {
  my ($this,$operand,%opt) = @_;
  my %lopt = $this->_LocalSwitch(%opt);
  my $debug = $lopt{-debug};
  my $dbg2  = $debug ? $debug-1 : undef;

  # create object copy
  if ($lopt{-copy}) {
    $this = &DataClone ( {%{$this}}, -debug=>$dbg2);
  }

  # recalculate by subtracting scalar operand
  if (! ref($operand)) {
    if (int ($operand)) {
      foreach (@{$this->{pdata}}) { $_->[1]-=$operand }
      delete $this->{SecData};
    }
  }

  # find first intervening border
  elsif (ref($operand) eq __PACKAGE__) {
    my (%ValX,@ValX,@PlotNew);
    foreach (sort{ $a<=>$b } map{ $_->[0] } map{ @{$_||[]} } $this->{pdata},$operand->Data()) {
      if ($ValX{$_}) { next }
      $ValX{$_} = 1;
      push @ValX, $_;
    }
    $debug and printf STDERR "%s. unique x values in both plots: %d\n", &MySub, int(@ValX);
    foreach my $itX (@ValX) {
      push @PlotNew, [ $itX,
        $this->Interpolate($itX,-debug=>$dbg2)
          - $operand->Interpolate($itX,-debug=>$dbg2) ];
    }
    $this->{pdata} = [ @PlotNew ];
    delete $this->{SecData};
  }

  # bad operand
  else {
    die printf "%s. ERROR: bad operand type %s\n", &MySub, ref($operand);
  }

  return $this;
}


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