#!/usr/bin/perl
use strict;
use warnings;
use diagnostics;

#Force the same set of random numbers for each run until it works...
#use Math::Random::MT qw(rand);
srand 0;

open (F1,'>>event_times_new2.mcm');
open (F2,'>>lambda_star_new2.mcm');

my ($lambda0, $alpha, $beta, $lambda_star);

my $time=8000000;
my $counterfails=0;
my $num=1;

$lambda0=1.2;
$alpha=0.8;
$beta=0.9;

#ALGO it self

my (@d, @s, @u, @event_times);
$lambda_star = $lambda0;
$d[0] = 0;
$u[0] = rand($num);
$s[0] = -log($u[0]) / $lambda_star;

die if $s[0] >= $time;
$event_times[0] = $s[0];

print F1 "$event_times[0]\n";
print F2 "$lambda_star\n";

my $valuelambda = lambda($event_times[0]) + $lambda0;
$lambda_star=$valuelambda+$alpha;

START:

$u[1]=rand($num);
$s[1] = $s[0] - log($u[1])/$lambda_star;
die if $s[1] > $time;

$d[1]=rand($num) ;
my $rejectionlambda=lambda($s[1]);
if ($d[1] < $rejectionlambda/$lambda_star) {
    $event_times[1]=$s[1];

    print F1 "$event_times[1]\n";
    print F2 "$lambda_star\n";
}else{
    $counterfails+=1;
    $lambda_star=$rejectionlambda;
    goto START;
}

print F1 "The number of event-times that failed the A/R procedure is $counterfails\n";

close F1;
close F2;

################################################

#subroutine that computes intensity function
sub lambda {
    my $t = $_[0];

    my $term;
    my $sum=0;
    for (my $i=0; $i<@event_times; $i++) {
        $term = $alpha*exp( -$beta * ($t - $event_times[$i]) );
        $sum += $term;
    }
    $sum -= $alpha if $term == $alpha;

    return $sum;
}