#!/usr/bin/perl
package EmpiricalDF;
use strict;
use warnings;
use PDL;
use PDL::Graphics::PGPLOT;
use PGPLOT;
use PDL::Slatec;
use PDL::NiceSlice;
use Carp;
use SimpleObject;

our @ISA = qw(SimpleObject);

sub new {
  my $c  = shift;
  my $class = ref($c) || $c;
  my @p = @{shift @_ } ; #probabilities, frequencies, ...
  my @v = @{shift @_}; #values
  my $type = shift || "qf";
  my $interpol = shift || "linear"; # pmch = piecewise monotonic cubic Hermite, step = stepwise constant
  my $options = shift; 
  $options->{discrete} = 0 unless (exists $options->{discrete});
  my $x; 
  my $y = pdl @v;
  if ($type eq "freq") { 
	my $tmp = pdl @p;
	$x = cumusumover ( $tmp ) / sum($tmp);
	
  }
  elsif ($type eq "qf") {
	$x = pdl @p;
  }
  my $self = $class->SUPER::create(
	    x => $x,
	    y => $y,
	    interpolation => $interpol,
	    init => 0,
	    options => $options,
	    pmch_end_d2 => 0.1
	  );
  return $self;
}
sub print_points {
    my $self = $_[0];
    foreach (0 .. $self->{x}->nelem-1) {
	  print $self->{x}->index($_). "\t". $self->{y}->index($_) . "\n" ;
    }

}


sub is_initialized {
    my $self = shift;
    return $self->{init};
}
sub initialize {
    my $self = shift;
    if ($self->{interpolation} eq "pmch") {
	  my $ic = short( 0, -2 );
	  my $vc = float( 0, $self->{pmch_end_d2} ); 
	  my $wk = zeroes( 2 * $self->{x}->getdim(0) );
	  my ($d, $e) = chic( $ic, $vc, 0, $self->{x}, $self->{y}, $wk );
	  $self->{init} = 1;
	  $self->{pmch_deriv} = $d;
    }
}
sub set_interpolation {
  my $self = shift;
  $self->{init} = 0;
  $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 qf_plot {
    my $self = shift;
    my $fn = shift;
    my ($lx, $ux, $ly, $uy) =  @_;
    if (defined $fn) {
      dev($fn.'/png');
    }
    else {
      dev ('/xs');
    }
    env($lx, $ux, $ly, $uy)  if (defined $lx and defined $ux  and defined $ly and defined $uy);
    pgscr(0,1,1,1); # set index 0 to white
    pgscr(1,0,0,0); # set index 1 to black
    hold
    points $self->{x}, $self->{y};
    my $xi = sequence (10001);
    $xi *= 0.0001;
    line $xi, $self->_i_qf($xi);
    release;
}
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');
}
sub qf_compare_plot {
    
    my ($self,$fn, $y, $lx, $ux, $ly, $uy) =  @_;
    my $ny = pdl (@$y);
    qsort(inplace ($ny));
    my $nx = ones($ny);
    inplace($nx / $ny->getdim(0));
    cumusumover(inplace($nx));
    if (defined $fn) {
      dev($fn.'/png');
    }
    else {
      dev ('/xs');
    }
    env($lx, $ux, $ly, $uy)  if (defined $lx and defined $ux  and defined $ly and defined $uy);
    pgscr(0,1,1,1); # set index 0 to white
    pgscr(1,0,0,0); # set index 1 to black
    hold
    points $nx , $ny(-1:0);
    my $xi = sequence (10001);
    $xi *= 0.0001;
    line $xi, $self->_i_qf($xi);
    release;
}
1;