#!/usr/bin/perl -w

#
# randomly generate process observations that simulate a
# binomially-distributed process in the state of statistical control
# (p_setup.csv) and simulate the same process experiencing a drift of
# magnitude $drift starting two hours into the $shift shift
# (p_monitoring.csv)
#

use strict;
use Math::Random;

my %shiftSchedule = (
    "first" =>  { "start" =>  6.00, "end" => 14.00 },
    "second" => { "start" => 14.00, "end" => 22.00 },
    "third" =>  { "start" => 22.00, "end" =>  6.00 }
);
my $shift = "third";         # shift to monitor
my $inspectionRate = 1 / 2;  # every 1/2 hour
my $drift = 1.5;             # sigma drift to simulate
my $m = 25;                  # samples in control chart setup
my $target = 0.10;           # fraction nonconforming target

my $p = $target;
my $n =                      # observations per sample
    int((9 * $p * (1 - $p)) / ($drift * $p * $drift * $p) + 0.5);
my $hour;
my $i;
my $minute;
my $observation;
my $setupM = $m;

print "n = $n\n";

#
# simulate control chart setup
#
open(SETUPCSV, ">p_setup.csv") || die "! can't open \"p_setup.csv\" ($!)\n";
for ($i = 1; $i <= $m; $i++) {
    $observation = Math::Random::random_binomial(1, $n, $p);
    print SETUPCSV $observation . "\r\n";
}
close(SETUPCSV);

#
# simulate control chart monitoring
#
open(MONITORINGCSV, ">p_monitoring.csv") || die "! can't open \"p_monitoring.csv\" ($!)\n";
$m = $shiftSchedule{$shift}{"end"} - $shiftSchedule{$shift}{"start"};
if ($m < 0) {
    $m += 24;
}
$m /= $inspectionRate;
for ($i = 1; $i <= $m; $i++) {
    $hour = int($i * $inspectionRate + $shiftSchedule{$shift}{"start"});
    if ($hour >= 24) {
        $hour -= 24;
    }
    $minute = ($i & 0x1) ? (60 * $inspectionRate) : 0;
    if ($i >= (0.25 * $m)) {
        if ($i < (0.75 * $m)) {
            $p = $target + ($drift * $target / (0.5 * $m)) * ($i - (0.25 * $m));
        } else {
            $p = $target + $drift * $target;
        }
    }
    $observation = Math::Random::random_binomial(1, $n, $p);
    printf MONITORINGCSV "'%d:%02d',%d\r\n", $hour, $minute, $observation;
}
close(MONITORINGCSV);