image thumbnail

pmap - Parameter Space Stability Mapping Suite

by

 

30 Mar 2011 (Updated )

Finds the Hurwitz stability region in a parameter space of a continuous time system.

Stable=pmap_calc(varargin)
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
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Contact us