% This script reads the coefficient files in the current directory.
% It checks the order condition and whether the function for 
% the positivity coefficient changes its sign near "C".
% The convention for the filename is 
%  'm' s ns p'.mat'
% where s is the number of stages, ns the number of shifted stages 
% and p the order.

% compatible to octave and matlab.

files = dir('m*.mat')

for fn = {files.name}
    fprintf('------- %s -------',fn{1});
    mf = load(fn{1})
    c = mf.c; A = mf.A; B = mf.B; R = mf.R; s = mf.s; ropt = mf.ropt;
    defect = [sum(B,2)-ones(s,1)];
    for j = 1:mf.p
        jf = factorial(j); 
        j1 = j-1;
        j1f = factorial(j1);
        defect = [defect,c.^j/jf-B*((c-1).^j/jf)-A*((c-1).^j1/j1f)-R*(c.^j1/j1f)];
    end
    defectInOrderConditions = norm(defect)
    G=@(r)min(min((eye(s)+r*R)\[A,R,B-r*A]));
    positivityNearC = arrayfun(G,ropt*[0.999,1,1.001])
end

