CyclicPoisson.m

classdef CyclicPoisson < ContinuousDistribution
    % CyclicPoisson Cyclic Poisson arrival process with deterministic phase switching.
    %
    % In phase i, arrivals follow Poisson(rates(i)) for exactly durations(i)
    % time units, then the rate switches deterministically to phase mod(i,n)+1.
    % The schedule repeats cyclically.
    %
    % Unlike a MAP, which uses exponential holding times to approximate the
    % phase schedule stochastically, this process fires rate-change events at
    % exact phase boundaries in the LDES simulation engine.
    %
    % The process representation stores:
    %   process{1} : 1-by-n row vector of arrival rates
    %   process{2} : 1-by-n row vector of phase durations
    %
    % Solver support. Only the LDES simulation engine honours the exact
    % phase schedule. CyclicPoisson is declared solely in the LDES feature
    % set, so every analytical solver (MVA, NC, CTMC, MAM, FLD, JMT, SSA)
    % rejects a model using it via the standard unsupported-feature check.
    %
    % Copyright (c) 2012-2026, Imperial College London
    % All rights reserved.
 
    properties
        rates;      % 1-by-n vector of arrival rates per phase
        durations;  % 1-by-n vector of phase durations
        process;    % {rates, durations} cell, for serialization
    end
 
    methods
        function self = CyclicPoisson(rates, durations)
            % SELF = CYCLICPOISSON(RATES, DURATIONS)
            self@ContinuousDistribution('CyclicPoisson', 0, [0, Inf]);
            if nargin < 2 || isempty(rates) || isempty(durations) ...
                    || numel(rates) ~= numel(durations)
                line_error(mfilename, 'CyclicPoisson: rates and durations must be non-empty and equal length');
            end
            self.rates = rates(:).';
            self.durations = durations(:).';
            self.process = {self.rates, self.durations};
            self.mean = 1.0 / self.getTimeAverageRate();
            self.immediate = false;
        end
 
        function rate = getTimeAverageRate(self)
            % RATE = GETTIMEAVERAGERATE()
            % Time-average arrival rate: sum(rates.*durations) / sum(durations).
            totalArrivals = sum(self.rates .* self.durations);
            totalDuration = sum(self.durations);
            if totalDuration > 0
                rate = totalArrivals / totalDuration;
            else
                rate = 0.0;
            end
        end
 
        function r = getRates(self)
            r = self.rates;
        end
 
        function d = getDurations(self)
            d = self.durations;
        end
 
        function n = getNumPhases(self)
            n = numel(self.rates);
        end
 
        function mean = getMean(self)
            % MEAN = GETMEAN()
            mean = 1.0 / self.getTimeAverageRate();
        end
 
        function rate = getRate(self)
            rate = self.getTimeAverageRate();
        end
 
        function scv = getSCV(self)
            % SCV = GETSCV() — placeholder (used only if Erlang conversion is attempted).
            scv = 1.0;
        end
 
        function skew = getSkewness(self)
            skew = 2.0;
        end
 
        function F = evalCDF(self, t)
            rate = self.getTimeAverageRate();
            F = 1.0 - exp(-rate * t);
        end
 
        function L = evalLST(self, s)
            rate = self.getTimeAverageRate();
            L = rate ./ (rate + s);
        end
 
        function proc = getProcess(self)
            proc = self.process;
        end
 
        function X = sample(self, n)
            % X = SAMPLE(N) — samples Exp(time-average rate). The exact phase
            % schedule is honoured only by the LDES simulator, not by sample().
            if nargin < 2
                n = 1;
            end
            rate = self.getTimeAverageRate();
            X = -log(1 - rand(n, 1)) / rate;
        end
    end
end