function Stable=pmap_calc(varargin)
% pmap_calc
%
% Version 1.0
% 3-18-11
%
% Iteratively bisects and tests a parameter space region for stability
% using vertex-polynomials (now companion matrices) which are generated
% from Sign-Definite Decomposition.
%
% [1] L.H. Keel and S.P. Bhattacharyya, "Fixed order multivariable
% controller synthesis: A new algorithm," Proceedings of the 47th IEEE
% Conference on Decision and Control, Cancun, Mexico, December, Cancun,
% Mexico: 2008.
%
% [2] M.J. Knap, L.H. Keel, and S.P. Bhattacharyya, "Robust Stability of
% Complex Systems with Applications to Performance Attainment
% Problems," Proceedings of the 2010 American Control Conference,
% Baltimore, MD, USA: 2010, pp. 3907-3912
%
%
% Copyright (C) 2011 Michael Knap
% FIXME: Throuroughly introduce the different parameter value pairs
%% First, parse input
p=inputParser;
p.addParamValue('nChunks',128,@(x) x>1 && mod(x,1)==0);
p.addParamValue('bisect1method','longest',@(str) any(strcmpi(str,{'longest','modulus','order'})));
p.addParamValue('bisect2method','longest',@(str) any(strcmpi(str,{'longest','modulus','order'})));
p.addParamValue('threshold',[],@(x) x>4 && x<30 && mod(x,1)==0);
p.addParamValue('prune1method','none',@(str) any(strcmpi(str,{'none','coef','routh'})));
p.addParamValue('prune2method','none',@(str) any(strcmpi(str,{'none','coef','routh'})));
p.addParamValue('order1',[],@(x) true) % TODO - validator for parsing order bisection
p.addParamValue('order2',[],@(x) true) % TODO - validator for parsing order bisection
p.addParamValue('debug',0, @(d) isscalar(d));
p.KeepUnmatched=true;
p.StructExpand=true;
p.parse(varargin{:});
% First we want to "import" our unmatched entries ... these are from the
% Test structure generated with pmap_generate
fs=fieldnames(p.Unmatched);
for k=1:length(fs)
Stable.(fs{k})=p.Unmatched.(fs{k});
end
% Then we want to get our resuls... these come from parsing the input
fs=fieldnames(p.Results);
for k=1:length(fs)
Stable.(fs{k})=p.Results.(fs{k});
end
if Stable.debug>3
disp('-------------------------------------------------------------------')
disp('pmap_calc')
disp('Sign-Definite decomposition to find stability in parameter space')
fprintf('\n');
disp('Parse Results (p.Rresults) :')
disp(p.Results)
fprintf('\n');
disp('Unmatched Parse Results (p.Unmatched) :')
disp(p.Unmatched)
disp('-------------------------------------------------------------------')
%keyboard
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% get our poolsize and Announce what we are doing for debugging
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Stable.n=matlabpool('size');
if Stable.debug>1 && Stable.n >1
fprintf('Using %d parallel workers.\n',Stable.n);
elseif Stable.n==0
warning('pmap:matlabpool','NOT RUNNING. pmap is meant to be ran with a matlabpool.')
%return;
end
% fs=fieldnames(Test);
% for k=1:length(fs)
% Stable.(fs{k})=Test.(fs{k});
% end
if Stable.debug>0
Stable.threshold=p.Results.threshold;
disp(['Using threshold : ' num2str(Stable.threshold)]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Now, here we want to get our phase interval for the following test :
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Note that I moved this code to here from the generate function.
% the idea is that I only wanted to generate the test structure once.
% In that way, I can pass the phasemargin to this function for different
% testing;
PhaseIntervals=GetPhaseIntervals(Stable.debug,Stable.phasemargin);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Chunking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Chunking is an initial division of the parameter space to facilitate
% paralellization and future divisional testing
% bisectmethod used for chunking
if Stable.debug>0
fprintf('Chunking... ');
end
tchunk=tic;
switch Stable.bisect1method
case 'modulus'
bisect=@bisectmod;
case 'longest'
bisect=@bisectlongest;
case 'order'
bisect=@bisectorder;
end
switch Stable.prune1method
case 'none'
condCheck=@(~) true;
case 'routh'
condCheck=@RouthCondCheck;
case 'coef'
condCheck=@CoefCondCheck;
end
chunks=cell(1,length(Stable.orthants)); % chunks are our orthants
parfor k=1:length(Stable.orthants)
% we give each chunk their own properties
chunks{k}.pts=Stable.orthants{k};
chunks{k}.threshold=Stable.threshold;
chunks{k}.Pcompan={Stable.Pcompan{PhaseIntervals,k}};
if isfield(Stable,'coefCond')
chunks{k}.cond=Stable.coefCond{k};
end
if isfield(Stable,'routhCond')
chunks{k}.rcond=Stable.routhCond{k};
end
if strcmpi(Stable.bisect1method,'order')
chunks{k}.order=Stable.order1;
end;
if strcmpi(Stable.bisect2method,'order')
chunks{k}.order2=Stable.order2;
end;
end
prune1=[];
% chunk loop
while length(chunks) < Stable.nChunks
if condCheck(chunks{1})
[l u]=bisect(chunks{1});
if ~isempty(l)
chunks=[chunks {l u}];
end
else
tmp.pts=chunks{1}.pts;
tmp.threshold=chunks{1}.threshold;
prune1=[prune1 tmp];
end
chunks(1)=[];
end
%Stable.chunks=chunks;
Stable.t.chunk=toc(tchunk);
Stable.prune1=prune1;
if Stable.debug>0
fprintf('Created %d chunks (%4.2f s).\n',length(chunks),Stable.t.chunk);
end;
% display([' Created ' num2str(length(chunks)) ' chunks in ' ...
% num2str(Stable.t.chunk) ' seconds with ' ...
% num2str(length(Stable.prune1)) ' pruned.']);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The Parallel Test
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tcalc=tic;
%disp(['Performing ' Stable.testmethod ' test ...']);
fprintf('Performing %s test... ',Stable.testmethod);
bm=Stable.bisect2method;
pm=Stable.prune2method;
switch(Stable.testmethod)
case 'vertex'
parfor k=1:length(chunks)
[s{k} prune2{k}]=HurwitzTest(chunks{k},...
bm,pm,Stable.phasemargin,Stable.gainmargin);
% if ~isempty(s{k})
% display('Got something here!');
% else
% display('Empty')
% end
end
case 'vertex_compan'
parfor k=1:length(chunks)
[s{k} prune2{k}]=HurwitzTest_compan(chunks{k},...
bm,pm,Stable.phasemargin,Stable.gainmargin);
end
otherwise
error(['Unknown test method: ' Stable.testmethod]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Finish up
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Stable.t.calc=toc(tcalc);
Stable.s=horzcat(s{:});
Stable.prune2=horzcat(prune2{:});
Stable.version.pmap=1.0;
Stable.date.pmap=datestr(now);
fprintf('Calculated %d stable regions (%4.2f s).\n',length(Stable.s),Stable.t.calc);
if Stable.debug>=5
disp(Stable);
end
% if Stable.debug==5
% testplot(Stable);
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Supporting Functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% validateTest
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TODO : write validateTest fcn
function bool = validateTest(~)
bool=true;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% bisectorder(region)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [lower upper]=bisectorder(r)
midpoint=@(a,b) (a+b)/2 ;
%% Only bisect if threshold >0
if r.threshold>0 % Only bisect if threshold >0
%% Bisect according to threshold modulus and either provided order or order of dimensions(default)
try
n=length(r.order);
e=mod(r.threshold,n);
bisectdim=r.order(e+1);
catch emsg
disp(emsg);
warning('bkk2:info','Order method was selected, but no order field provided. defaulting to mod')
e=mod(r.threshold, length(r.pts));
bisectdim=length(r.pts)-e;
end
%% The bisection redefines upper and lower pts from midpoint
lower=r;
upper=r;
%
% lower.pts=r.pts;
% upper.pts=r.pts;
lower.threshold=r.threshold-1;
upper.threshold=r.threshold-1;
% lower.order=r.order;
% upper.order=r.order;
%
lower.pts(1,bisectdim)=midpoint(r.pts(1,bisectdim),r.pts(2,bisectdim));
upper.pts(2,bisectdim)=midpoint(r.pts(1,bisectdim),r.pts(2,bisectdim));
else
lower=[];
upper=[];
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% bisectlongest
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [lower upper]=bisectlongest(r)
midpoint=@(a,b) (a+b)/2 ;
% Only bisect if threshold >0
if r.threshold>0 % Only bisect if threshold >0
% Bisect longest edge
[~,bisectdim]=max(abs(r.pts(1,:)-r.pts(2,:)));
% The bisection redefines upper and lower pts from midpoint
lower=r;
upper=r;
%
% lower.pts=r.pts;
% upper.pts=r.pts;
lower.threshold=r.threshold-1;
upper.threshold=r.threshold-1;
% lower.order=r.order;
% upper.order=r.order;
%
lower.pts(1,bisectdim)=midpoint(r.pts(1,bisectdim),r.pts(2,bisectdim));
upper.pts(2,bisectdim)=midpoint(r.pts(1,bisectdim),r.pts(2,bisectdim));
else
lower=[];
upper=[];
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% bisectmod(region)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [lower upper]=bisectmod(r)
midpoint=@(a,b) (a+b)/2 ;
% Only bisect if threshold >0
if r.threshold>0 % Only bisect if threshold >0
% Bisect according to threshold modulus and either provided order or order of dimensions(default)
e=mod(r.threshold, length(r.pts));
bisectdim=length(r.pts)-e;
% The bisection redefines upper and lower pts from midpoint
lower=r;
upper=r;
%
% lower.pts=r.pts;
% upper.pts=r.pts;
lower.threshold=r.threshold-1;
upper.threshold=r.threshold-1;
% lower.order=r.order;
% upper.order=r.order;
%
lower.pts(1,bisectdim)=midpoint(r.pts(1,bisectdim),r.pts(2,bisectdim));
upper.pts(2,bisectdim)=midpoint(r.pts(1,bisectdim),r.pts(2,bisectdim));
else
lower=[];
upper=[];
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PhaseIntervals = GetPhaseIntervals(phasemargin)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIXME : This has got to be more elegant !
% Right now, this is a very crude hack just to get some testing
% done, it is not full-featured or complete.
function PhaseIntervals = GetPhaseIntervals(debug,margin)
if debug>3
fprintf('Getting phaseintervals...');
end
I{1}=[-180 -90];
I{2}=[-90 0];
I{3}=[0 90];
I{4}=[90 180];
if margin(1)==margin(2)
Imat=vertcat(I{:});
bools = margin(1)>=Imat(:,1) & margin(1)<=Imat(:,2);
PhaseIntervals=bools;
PhaseIntervals=[false false true false];
else
fprintf('Checking lower interval ')
for i=1:4
fprintf('%d',i);
if margin(1)>=I{i}(1) && margin(1)<I{i}(2)
fprintf('\n');
I1=i;
break;
end
fprintf('...');
end
fprintf('Checking upper interval ');
for i=1:4
fprintf('%d',i);
if margin(2)<=I{i}(2) && margin(2)>I{i}(1)
fprintf('\n');
I2=i;
break;
end
fprintf('...');
end
PhaseIntervals=I1:I2
end
if debug>3
fprintf(' Done.\n');
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% HurwitzTest_compan
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is my latest addition to these algorithms and programs. Though I
% haven't tested for numerical instability, this is now the default. The
% use of companion matrices has made this code about 3 times faster.
%
function [StableRegions PruneRegions]=HurwitzTest_compan(orthant,...
bisectmethod,prunemethod,pmarg,gmarg) %,bisectmethod)
% We want to pass a smaller amount of data around, so we define a region
region.pts=orthant.pts;
region.threshold=orthant.threshold;
% Assign bisect method
switch bisectmethod
case 'modulus'
bisect=@bisectmod;
case 'longest'
bisect=@bisectlongest;
case 'order'
bisect=@bisectorder;
region.order=orthant.order2;
end
% Assign prune method
switch prunemethod
case 'none'
condCheck=@(varargin) true;
case 'routh'
condCheck=@(varargin) orthant.rcond(varargin{:}) > 0;
case 'coef'
condCheck=@(varargin) orthant.cond(varargin{:}) > 0;
end
% Initialize vars
to_test=region;
StableRegions=[];
PruneRegions=[];
% The loop:
while isempty(to_test)~=true
region=to_test(1);
to_test(1)=[];
% this block of code is used frequently, so ...
% TODO - vectorize the getting of min and max magnitudes from pts
[~, minI]=min(abs(region.pts));
[~, maxI]=max(abs(region.pts));
for j=1:length(minI)
minsandmaxs{j}=region.pts(minI(j),j);
minsandmaxs{j+length(minI)}=region.pts(maxI(j),j);
end
% This is where the phasemargin and gainmargin comes in %
args=[minsandmaxs, ...
{pmarg(1)*pi/180, pmarg(2)*pi/180, gmarg(1), gmarg(2)} ];
% 1 for min and max gain, zero for degree will change
if condCheck(args{:})
% old way:
% parfor k=1:8
% r(:,k)=roots(orthant.polies{k}(args{:}));
% end
%r=cellfun(@eig,orthant.Pcompan(args{:}))
%r=eig(orthant.Pcompan(args{:}));
for j=1:length(orthant.Pcompan)
rp{j}=eig(orthant.Pcompan{j}(args{:}));
end
r=vertcat(rp{:});
if real(r) < 0 % for continuous case
StableRegions=[StableRegions region];
else
[l u]=bisect(region);
to_test=[to_test; l ; u];
end
else
PruneRegions=[PruneRegions region];
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% %% HurwitzTest
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % TODO - can the diff. tests be consolidated into one ?
% function [StableRegions PruneRegions]=HurwitzTest(orthant,...
% bisectmethod,prunemethod,pmarg,gmarg) %,bisectmethod)
% % We want to pass a smaller amount of data around, so we define a region
% region.pts=orthant.pts;
% region.threshold=orthant.threshold;
%
% % Assign bisect method
% switch bisectmethod
% case 'modulus'
% bisect=@bisectmod;
% case 'longest'
% bisect=@bisectlongest;
% case 'order'
% bisect=@bisectorder;
% region.order=orthant.order2;
% end
%
% % Assign prune method
% switch prunemethod
% case 'none'
% condCheck=@(varargin) true;
% case 'routh'
% condCheck=@(varargin) orthant.rcond(varargin{:}) > 0;
% case 'coef'
% condCheck=@(varargin) orthant.cond(varargin{:}) > 0;
% end
%
% % Initialize vars
% to_test=region;
% StableRegions=[];
% PruneRegions=[];
%
% % The loop:
% while isempty(to_test)~=true
% region=to_test(1);
% to_test(1)=[];
%
% % TODO - vectorize the getting of min and max magnitudes from pts
% [~, minI]=min(abs(region.pts));
% [~, maxI]=max(abs(region.pts));
% for j=1:length(minI)
% minsandmaxs{j}=region.pts(minI(j),j);
% minsandmaxs{j+length(minI)}=region.pts(maxI(j),j);
% end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % FIXME: This is where the phasemargin and gainmargin comes in %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% args=[minsandmaxs, ...
% {pmarg(1)*pi/180, pmarg(2)*pi/180, gmarg(1), gmarg(2)} ];
%
%
% % 1 for min and max gain, zero for degree will change
%
% if condCheck(args{:})
% parfor k=1:8
% r(:,k)=roots(orthant.polies{k}(args{:}));
% end
%
% if real(r) <0 % for continuous case
% %if abs(r) <1 % for discrete case (needs a whole new approach)
% StableRegions=[StableRegions region];
% else
% [l u]=bisect(region);
% to_test=[to_test; l ; u];
% end
% else
% PruneRegions=[PruneRegions region];
% end
% end
% end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% %% GainTest
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [StableRegions]=GainTest(orthant,bisectmethod,prunemethod,gainmargin) %,bisectmethod)
% % We want to pass a smaller amount of data around, so we define a region
% region.pts=orthant.pts;
% region.threshold=orthant.threshold;
%
% % Assign bisect method
% switch bisectmethod
% case 'modulus'
% bisect=@bisectmod;
% case 'longest'
% bisect=@bisectlongest;
% case 'order'
% bisect=@bisectorder;
% end
%
% % Initialize vars
% to_test=region;
% StableRegions=[];
% %PrunedRegions=[];
%
% % The loop:
% while isempty(to_test)~=true
% region=to_test(1);
% to_test(1)=[];
% [~, minI]=min(abs(region.pts));
% [~, maxI]=max(abs(region.pts));
% for j=1:length(minI)
% minsandmaxs{j}=region.pts(minI(j),j);
% minsandmaxs{j+length(minI)}=region.pts(maxI(j),j);
% end
% args=[minsandmaxs {gainmargin(1),gainmargin(2) ,0}]; % 1 for min and max gain, zero for degree will change
%
%
%
% % p1=orthant.p1Fcn(minsandmaxs{:});
% % p2=orthant.p2Fcn(minsandmaxs{:});
% % p3=orthant.p3Fcn(minsandmaxs{:});
% % p4=orthant.p4Fcn(minsandmaxs{:});
% %
% % q1=orthant.q1Fcn(minsandmaxs{:});
% % q2=orthant.q2Fcn(minsandmaxs{:});
% % q3=orthant.q3Fcn(minsandmaxs{:});
% % q4=orthant.q4Fcn(minsandmaxs{:});
% %
% % Phase 2 pruning will be here, but right now, it is a little buggy
% % Prune the dataset by checking the coefficients, if all vertex
% % below line, then we can remove the
% % if ~conditiontest(orthant,region) % if the conditiontest fails for the region, then prune it
% % PrunedRegions=[PrunedRegions region];
% % else
%
% parfor k=1:8
% r(:,k)=roots(orthant.polies{k}(args{:}));
% end
%
% if real(r) <0
% StableRegions=[StableRegions region];
% else
% [l u]=bisect(region);
% to_test=[to_test; l ; u];
% end
% end
% end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% %% GainPhaseTest
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [StableRegions]=GainPhaseTest(orthant,bisectmethod,prunemethod,gainmargin,phasemargin) %,bisectmethod)
% % We want to pass a smaller amount of data around, so we define a region
% region.pts=orthant.pts;
% region.threshold=orthant.threshold;
%
% % Assign bisect method
% switch bisectmethod
% case 'modulus'
% bisect=@bisectmod;
% case 'longest'
% bisect=@bisectlongest;
% case 'order'
% bisect=@bisectorder;
% end
%
% % Initialize vars
% to_test=region;
% StableRegions=[];
% %PrunedRegions=[];
%
% % The loop:
% while isempty(to_test)~=true
% region=to_test(1);
% to_test(1)=[];
%
% [~, minI]=min(abs(region.pts));
% [~, maxI]=max(abs(region.pts));
% for j=1:length(minI)
% minsandmaxs{j}=region.pts(minI(j),j);
% minsandmaxs{j+length(minI)}=region.pts(maxI(j),j);
% end
% for d=1:(phasemargin(2)-phasemargin(1))
% args=[minsandmaxs {gainmargin(1),gainmargin(2) ,d-1}]; % 1 for min and max gain, zero for degree will change
% parfor k=1:8
% r(:,d,k)=roots(orthant.polies{k}(args{:}));
% end
% end
% if real(r) <0
% StableRegions=[StableRegions region];
% else
% [l u]=bisect(region);
% to_test=[to_test; l ; u];
% end
% end
% end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% %% PhaseTest
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [StableRegions PruneRegions]=PhaseTest(orthant,bisectmethod,prunemethod,phasemargin) %,bisectmethod)
% % We want to pass a smaller amount of data around, so we define a region
% region.pts=orthant.pts;
% region.threshold=orthant.threshold;
%
% % Assign bisect method
% switch bisectmethod
% case 'modulus'
% bisect=@bisectmod;
% case 'longest'
% bisect=@bisectlongest;
% case 'order'
% bisect=@bisectorder;
% end
%
% % Initialize vars
% to_test=region;
% StableRegions=[];
% PruneRegions=[];
%
% % The loop:
% while isempty(to_test)~=true
% region=to_test(1);
% to_test(1)=[];
% [~, minI]=min(abs(region.pts));
% [~, maxI]=max(abs(region.pts));
% for j=1:length(minI)
% minsandmaxs{j}=region.pts(minI(j),j);
% minsandmaxs{j+length(minI)}=region.pts(maxI(j),j);
% end
% for d=1:(phasemargin(2)-phasemargin(1))+1
% args=[minsandmaxs {1,1 ,d-1}]; % 1 for min and max gain, zero for degree will change
%
% parfor k=1:8
% r(:,d,k)=roots(orthant.polies{k}(args{:}));
% end
% end
% if real(r) <0
% StableRegions=[StableRegions region];
% else
% [l u]=bisect(region);
% to_test=[to_test; l ; u];
% end
% end
% end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% %% RouthCondCheck - for chunking
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function bool = RouthCondCheck(orth)
% % get our args for the condition check:
% [~, minI]=min(abs(orth.pts));
% [~, maxI]=max(abs(orth.pts));
% for j=1:length(minI)
% minsandmaxs{j}=orth.pts(minI(j),j);
% minsandmaxs{j+length(minI)}=orth.pts(maxI(j),j);
% end
% args=[minsandmaxs {1,1,0}]; % 1 for min and max gain, zero for degree, will change
% bool=orth.rcond(args{:}) > 0;
% end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% %% CoefCondCheck - for chunking
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function bool=CoefCondCheck(orth)
% % Get our args for the condition check:
% % TODO : Can a vectorized function for this be created ( This snippet is used elsewhere)
% [~, minI]=min(abs(orth.pts));
% [~, maxI]=max(abs(orth.pts));
% for j=1:length(minI)
% minsandmaxs{j}=orth.pts(minI(j),j);
% minsandmaxs{j+length(minI)}=orth.pts(maxI(j),j);
% end
% args=[minsandmaxs {1,1,0}]; % 1 for min and max gain, zero for degree, will change
%
% bool= orth.cond(args{:}) > 0;
% end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% %% DiscreteTest - nonsense
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % TODO - can the diff. tests be consolidated into one ?
% % Added Feb 10 2011 just because I wanted something to try. Will compare
% % with Tim. THIS DOESN'T EVEN MAKE SENSE. A test like this would require a
% % different method. I think Tim was doing Cebyshev deconstruction.
%
% % function [StableRegions PruneRegions]=DiscreteTest(orthant,bisectmethod,prunemethod) %,bisectmethod)
% % % We want to pass a smaller amount of data around, so we define a region
% % region.pts=orthant.pts;
% % region.threshold=orthant.threshold;
% %
% % % Assign bisect method
% % switch bisectmethod
% % case 'modulus'
% % bisect=@bisectmod;
% % case 'longest'
% % bisect=@bisectlongest;
% % case 'order'
% % bisect=@bisectorder;
% % region.order=orthant.order2;
% % end
% %
% % % Assign prune method
% % switch prunemethod
% % case 'none'
% % condCheck=@(varargin) true;
% % case 'routh'
% % condCheck=@(varargin) orthant.rcond(varargin{:}) > 0;
% % case 'coef'
% % condCheck=@(varargin) orthant.cond(varargin{:}) > 0;
% % end
% %
% % % Initialize vars
% % to_test=region;
% % StableRegions=[];
% % PruneRegions=[];
% %
% % % The loop:
% % while isempty(to_test)~=true
% % region=to_test(1);
% % to_test(1)=[];
% %
% % % TODO - vectorize the getting of min and max magnitudes from pts
% % [~, minI]=min(abs(region.pts));
% % [~, maxI]=max(abs(region.pts));
% % for j=1:length(minI)
% % minsandmaxs{j}=region.pts(minI(j),j);
% % minsandmaxs{j+length(minI)}=region.pts(maxI(j),j);
% % end
% % args=[minsandmaxs {1,1,0}]; % 1 for min and max gain, zero for degree will change
% %
% % if condCheck(args{:})
% % parfor k=1:8
% % r(:,k)=roots(orthant.polies{k}(args{:}));
% % end
% %
% % if abs(r) <1
% % StableRegions=[StableRegions region];
% % else
% % [l u]=bisect(region);
% % to_test=[to_test; l ; u];
% % end
% % else
% % PruneRegions=[PruneRegions region];
% % end
% % end
% % end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
