#!/usr/bin/perl
package EmpiricalDFofShares;
use strict;
use warnings;
use PDL;
use PDL::NiceSlice;
use Carp;
use FixedSizeTable;
use SimpleObject;

our @ISA = qw(SimpleObject);

sub new {
  my $c  = $_[0];
  my $class = ref($c) || $c;
  my @p = @$_[1];
  my @v = @$_[2];
  my $options = $_[3]; 
  $options->{discrete} = 0 unless (exists $options->{discrete});
  $options->{interpol} = "linear" unless (exists $options->{interpol});
  my $x = pdl @p;; 
  my @y;
  my $notfirst = 0;
  my $num = 0;
  for my $i ( @v) {
	if (ref($i) ne "ARRAY") {
		print STDERR "Error in EmpiricalDFofShares->new, not an array reference\n";
		exit(1);
	}
	if (scalar) {
		print STDERR "Error in EmpiricalDFofShares->new, not an array reference\n";
		exit(1);
	}
	$notfirst = 1;
	push @y , pdl (@$i);
  }
   = map {pdl $_ } ;
  
  my $self = $class->SUPER::create(
	    x => $x,
	    y => $y,
	    interpolation => $interpol,
	    init => 0,
	    options => $options,
	  );
  return $self;
}
sub print_points {
    my $self = $_[0];
    foreach (0 .. $self->{x}->nelem-1) {
	  print $self->{x}->index($_). "\t". $self->{y}->index($_) . "\n" ;
    }

}



sub set_interpolation {
  my $self = shift;
  $self->{interpolation} = shift;
}
sub set_pmch_end_d2 {
  my $self = shift;
  $self->{init} = 0;
  $self->{pmch_end_d2} = shift;
}
sub qf {
  my $self = shift;
  my $xi = pdl (@_);
  my $yi = $self->_i_qf($xi);
  return wantarray ? list($yi) : sclr ($yi->slice('0')) ;
}
sub _i_qf {
  my ($self, $xi) = @_;
  my ($yi, $e2);
  if ($self->{interpolation} eq "pmch") {
      $self->initialize unless $self->is_initialized;
      ($yi, $e2) =  chfe($self->{x}, $self->{y}, $self->{pmch_deriv}, 0, $xi);
  }
  elsif ($self->{interpolation} eq "linear") {
     $yi = interpol($xi, $self->{x}, $self->{y});
  }
  elsif ($self->{interpolation} eq "step") {
     $yi = $self->{y}->(vsearch($xi, $self->{x}) );
  }
  rint(inplace $yi) if ($self->{options}->{discrete});
  return $yi;
}
sub qf_interval {
  my $self = shift;
  my $xi = pdl (@_);
  my $ix = vsearch($xi, $self->{x}) ;
  my $M = cat ($ix , $ix -1);
  my $Yi = $self->{y}($M);
  return ( [ list($Yi(,0))], [ list($Yi(,1))]);
}
sub qf_integral {
 my $self = shift;
 my $lb = shift;
 my $ub = shift;
 my $integral;
 my $e3;
 if ($self->{interpolation} eq "pmch") {
      $self->initialize unless $self->is_initialized;
      ($integral, $e3) =  chia($self->{x}, $self->{y}, $self->{pmch_deriv}, 0,  $lb , $ub);
 }
 elsif ($self->{interpolation} eq "linear") {
     my $x = $self->{x};
     my $xi = cat (pdl ($lb), $x->where(($x > $lb) & ($x < $ub)), pdl ($ub));
     my $A = $xi->lags(0,1,2);
     $integral = inner(($A->slice(':,0') - $A->slice(':,1')) , ( 0.5 * (interpol($A->slice(':,0'), $self->{x}, $self->{y})  + interpol($A->slice(':,1'), $self->{x}, $self->{y}) )));
  }
  elsif ($self->{interpolation} eq "step") {
     my $x = $self->{x};
     my $xi = cat (pdl ($lb), $x->where(($x > $lb) & ($x < $ub)), pdl ($ub));
     my $A = $xi->lags(0,1,2);
     my $d = $A(:,0) - $A(:,1);
      $integral = inner($d, $self->qf($xi(1:)) );
  }
  return sclr($integral);
}
sub qf_sum_for_origN {
  my $self = shift;
  my $n = shift;
  croak ("N has to be larger than 0 (passed: $n)\n") unless ($n  > 0);
  my $lb = shift;
  my $ub = shift;
  $lb = 0 unless defined ($lb);
  $ub = 1 unless defined ($ub);
  my $xi = sequence($n);
  $xi += 1;
  $xi *=  ($ub -$lb) /$n;
  $xi += $lb;
  #print $xi->slice('1:5')."\n";
  my $yi = $self->_i_qf($xi);
  #print $yi ."\n";;
  return sum ( $yi);
}

sub lowest_nonzero {
    my $self = shift;
    return 0 unless ($self->{options}->{discrete} or $self->{interpolation} eq "step");
    return 1 unless ($self->{interpolation} eq "step");
    return $self->{y}->where($self->{y} > 0)->slice('0');
}

1;
