unit mcobj;
 {UTILITIES FOR MONTE CARLO METHODS - (RANDOM NUMBER GENERATORS) }
 { ***************************************** }
{Modified from standard Pascal toolbox in
          Wolfgang Kreutzer,
          System Simulation, programming styles and Languages, 1986,
          ISBN 0-201-12914-0, p 245 ff.}
{with OBJECTS For Turbo Pascal 6.0, 7.0 and Borland Pascal 7.0}
{Modified by:
          Lin Jensen
          Bishop's University
          September 1996     }

interface
uses clock,   {ONLY FOR REPORTS to print the simulation clock time
               see procedure RNG.Report, if you would rather not use this.}
     SIMobj;  {for global reporting}

 TYPE               {Data Type Declarations}

{***** RANDOM GENERATORS to sample from probability distributions}

RNG = object(aSimObject)             {"random number generator" - the common features}
     CONSTRUCTOR Init ( seedling : longint);
 {    procedure   SpecificLabel (CallMe : STRING);   }
          procedure reset;       VIRTUAL;{all the same except Description}
          procedure report;      VIRTUAL;
   private
                 seed          : longint;
                 samplesTaken  : 0..MAXINT;
                 integral      : REAL;
          function  U01 : real;         {uniform dist. in [0..1) }
          function  MeanOfQ : REAL;
end;

RandomDraw = OBJECT (rng)
     CONSTRUCTOR Init ( seedling : longint;
                        percentP : real);
     function  Draw : BOOLEAN;
   private
      percentTRUE   : REAL;
end;

RandomInteger = OBJECT (RNG)
     CONSTRUCTOR Init ( seedling : longint;
                        lowP, highP: INTEGER);
     function RandInt  : INTEGER;
  private
                    low, high     : INTEGER ;
                   end;

 RandomUniform = OBJECT(rng)
     CONSTRUCTOR Init ( seedling : longint;
                        lowP, highP: REAL);
     function Uniform  : REAL;
  private
         low, high     : REAL;
  end;

RandomNegExp = OBJECT      (rng)
     CONSTRUCTOR Init ( seedling : longint;
                        meanRateP: REAL);
     function NegExp : REAL;
  private
                   lambda : REAL ;
  end;

RandomNormal = OBJECT (rng)
     CONSTRUCTOR Init ( seedling : longint;
                        meanP, sigmaP: REAL);
     function Normal  : REAL;
  private
                  mean, sigma : REAL;
  end;

 { * * * * * EXECUTION MONITOR for conducting sampling experiments     }
 {this data type serves no purpose here. It is only included to allow  }
 {declarations of statistics collection classes to compile as one unit }

 aStateSamplingMonitorType = OBJECT
                              samplesTaken, toBeTaken : 0..MAXINT ;
                             end;

 { *******************************************}
 { * TheMonitor is assumed to be a GLOBAL of *}
 { * aStateSamplingMonitorType               *}
 { *******************************************}

 {>>>>>  NOT SURE THIS IS REQUIRED ANYWHERE - Lin Jensen, Sept 1996 <<<<<<}

{Utility Functions and Procedures for Monte Carlo Methods }
{ = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = }

{a" static" clock must be instantiated to allow the declarations}
{of some statistics gathering devices to compile!
   - This is done in unit CLOCK}

IMPLEMENTATION

{ * * * * * Messages understood by a state SAMPLING MONITO R * * * * * * * * * }

procedure InitStateSamplingMonitor (var sampler : aStateSamplingMonitorType);
var
 samplesToTake: INTEGER;

begin
  if (samplesToTake > 0) then
  sampler.toBeTaken := samplesToTake
  else
    WRITELN ('only positive number of samples please !');
    sampler.samplesTaken := 0;
end;

function AnotherSampleP (var sampler: aStateSamplingMonitorType): BOOLEAN;

{NOTE: this function has a side effect on the monitor's <samplesTaken> }
{ field !!! Each time it is used it increments it by one. }

begin
  sampler.samplesTaken := sampler.samplesTaken + 1;
  AnotherSampleP := sampler.samplesTaken <= sampler.toBeTaken;
end;


function ReportOnSamplesTakenQ (sampler : aStateSamplingMonitorType) :INTEGER;

begin
  ReportOnSamplesTakenQ := sampler. samplesTaken - 1;
end;

{***************** Methods for RANDOM SAMPLING ***************************}
{these methods are based on descriptions in Fishman (1978) p.345... where }
{generators for other distributions can also be found }


CONSTRUCTOR RNG.Init ( seedling : longint);
  begin
    aSimObject.Init('generator');
    seed := seedling;
    samplesTaken:= 0;
    integral := 0.0;
  end;

{"private" method used by other functions; Kreutzer used a
 0-1 generator adopted from Sedgewick ( 1982).
 we are using Turbo Pascal's Random}

function RNG.U01 : REAL;
begin
  Randseed := seed;    {our seed to Pascal's global}
  U01 := Random;       {Pascal's generator}
  Seed := Randseed;
  inc(samplesTaken);
end; { GenerateZeroToOneNumber function }

procedure RNG.Reset;
begin
 aSimObject.reset;
 samplesTaken := 0;
 integral := 0.0;
end;

function RNG.MeanOfQ : REAL;
begin
  if samplestaken > 0 then MeanOfQ := integral / samplesTaken;
end;

procedure RNG.Report;
begin
  aSimObject.report;
  WRITELN ( samplesTaken: 5,' samples were taken.');
  WRITELN ('The sample mean was', MeanOfQ :7:2);
end;

constructor RandomDraw.Init ( seedling : longint;
                             percentP: REAL);
begin
  if(percentP > 1.0) or (percentP < 0) then
    begin
      WRITELN ('Percentage parameter out of range ! (0 to 1 only)');
      halt;
    end
  else
    begin
      RNG.init (seedling);
      percentTrue := percentP;
      SpecificLabel ( 'Binomial (T-F)');
    end;
end; {of procedure InitDrawGen}


constructor RandomInteger.Init ( seedling : longint;
                             lowP, highP: INTEGER);
begin
 if lowP > highP then
 WRITELN ('bounds parameters out of range ! (low must be < = high)')
 else
   begin
     rng.init (seedling);
     SpecificLabel ( 'uniform integer');
     low := lowP;
     high := highP;
   end;
 end; {of procedure InitRandIntGen }

constructor RandomUniform.Init ( seedling : longint;
                             lowP, highP: REAL);
begin
 if lowP > highP then
 WRITELN ('bounds parameters out of range ! (low must be < = high)')
 else
   begin
     rng.init(seedling);
     SpecificLabel ( 'uniform real');
     low := lowP;
     high := highP;
   end;
end; {of procedure InitUniformGen}

constructor RandomNegExp.Init ( seedling : longint;
                              meanRateP: REAL);

begin
  if(meanRateP <= 0) then
    begin
      WRITELN ('Arrival rate for Neg.Exp. out of range ! ( > 0 only)');
      halt;
    end
  else
    begin
      rng.init(seedling);
      SpecificLabel ( 'negative exponential');
      lambda := meanRateP;
    end;
end; {of constructor InitNegExpGen}

constructor RandomNormal.Init ( seedling     : longint;
                              meanP, sigmaP: REAL);
begin
  if sigmaP < 0 then
  WRITELN ('the sigma value of a normal distribution must be positive !')
  else
   begin
     rng.init(seedling);
     SpecificLabel ( 'normal');
     mean := meanP;
     sigma := sigmaP;
   end;
end; {of constructor InitNormalGen}

function RandomDraw.Draw  : BOOLEAN;

{NOTE: this function produces a side effect! It changes the <seed> and }
 { statistics fields of the generator }
var
 temp: BOOLEAN;
begin
 if U01 <= percentTRUE then
 temp:= TRUE
 else
 temp:= FALSE;
 if temp then
 integral := integral + 1;
 Draw:= temp;
end; {of DRAW function}

function RandomInteger.RandInt  : INTEGER;
    {NOTE: this function produces a side effect. It changes the <seed> and }
    { statistics fields of the generator  }
var
 temp :INTEGER;

begin
  temp := TRUNC( (high - low + 1) * U01 + low );
  integral := integral + temp;
  RandInt := temp;
end; {of RandInt function }


function RandomUniform.Uniform  : REAL;
     {NOTE:  this function produces a side effect. It changes the <seed> and    }
     {           statistics fields of the generator                                     }
var
 temp: REAL;

begin
  temp:= ((high - low) * U01 + low);
  integral := integral + temp;
  Uniform := temp;
end; {of Uniform function}

function RandomNegExp.NegExp  : REAL;
     {NOTE:  this function produces a side effect. It changes the <seed> and    }
     {           statistics fields of the generator                                       }
var
 temp: REAL;

begin
  temp             := - LN (U01) /
                        lambda;
   integral        :=  integral + temp;
  NegExp          := temp;
end; { oNegExp fnction}

function RandomNormal.Normal : REAL;
      {NOTE: this function produces a side effect. It changes the <seed> and } { statistics fields of the generator }
 var
   temp : REAL;
   firstRN, secondRN: REAL; { the function uses a number of RNs to generate }
                            {one normally distributed RV (see Fisan (1978)}
   savesamples : integer;
 begin
   savesamples := samplesTaken;     {to only count one sample in RNG}
   repeat
     firstRN := - LN (U01);
     secondRN:= - LN (U01);
   until (2 * firstRN) >= (SQR (secondRN - 1 ));
   samplesTaken  := savesamples;     {to only count one sample in RNG}
   if (U01 < 0.5) then             {count one here}
   secondRN:= - secondRN;
   temp := sigma * secondRN + mean;
    integral :=  integral + temp;
  Normal := temp;
end; { of ;Normal function}

begin
end.