text:
% pop_epoch() - Convert a continuous EEG dataset to epoched data by extracting
%               data epochs time locked to specified event types or event indices.
%               May also sub-epoch an already epoched dataset (if sub-epochs are
%               same size or smaller). This pop_function calls epoch().
% Usage:
%   >> OUTEEG = pop_epoch( EEG); % pop-up a data entry window
%   >> OUTEEG = pop_epoch( EEG, events, timelimits);
%   >> [OUTEEG, indices] = pop_epoch( EEG, typerange, timelimits,'key1', value1 ...);
%
% Graphic interface:
%   "Time-locking event type(s)" - [edit box] Select 'Edit > Event values'
%                to see a list of event.type values; else use the push button.
%                To use event types containing spaces, enter in single-quotes.
%                epoch() function command line equivalent: 'typerange'
%   "..."      - [push button] scroll event types.
%   "Epoch limits" - [edit box] epoch latency range [start, end] in seconds relative
%                to the time-locking events. epoch() function equivalent: 'timelim'
%   "Name for the new dataset" - [edit box]
%                epoch() function equivalent: 'newname'
%   "Out-of-bounds EEG ..." - [edit box] Rejection limits ([min max], []=none).
%                epoch() function equivalent: 'valuelim'
% Inputs:
%   EEG        - Input dataset. Data may already be epoched; in this case,
%                extract (shorter) subepochs time locked to epoch events.
%   typerange  - Cell array of event types to time lock to. 'eventindices'
%                {default {} --> time lock epochs to any type of event}
%                (Note: An event field called 'type' must be defined in the
%                'EEG.event' structure. The command line argument is
%                'eventindices' below).
%   timelim    - Epoch latency limits [start end] in seconds relative to
%                the time-locking event {default: [-1 2]}
%
% Optional inputs:
%   'eventindices'- [integer vector] Extract data epochs time locked to the
%                indexed event numbers.
%   'valuelim' - [min max] or [max]. Lower and upper bound latencies for
%                trial data. Else if one positive value is given, use its
%                negative as the lower bound. The given values are also
%                considered outliers (min max) {default: none}
%   'verbose'  - ['yes'|'no'] {default: 'yes'}
%   'newname'  - [string] New dataset name {default: "[old_dataset] epochs"}
%   'epochinfo'- ['yes'|'no'] Propagate event information into the new
%                epoch structure {default: 'yes'}
%
% Outputs:
%   OUTEEG     - output dataset
%   indices    - indices of accepted events
%
% Authors: Arnaud Delorme and Hilit Serby, SCCN, INC, UCSD, 2001
%
% See also: eeglab, epoch

% deprecated
%   'timeunit' - Time unit ['seconds'|'points'] If 'seconds,' consider events
%                times to be in seconds. If 'points,' consider events as
%                indices into the data array. {default: 'points'}

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% 01-25-02 reformated help & license -ad
% 02-13-02 introduction of 'key', val arguments -ad
% 02-13-02 rereferencing of events -ad
% 03-18-02 interface and debugging -ad
% 03-27-02 interface and debugging -ad & sm

function [EEG, indices, com] = pop_epoch( EEG, events, lim, varargin );

if nargin < 1
   help pop_epoch;
	return;
end;
com = '';
indices = [];

if isempty(EEG.event)
    if EEG.trials > 1 & EEG.xmin <= 0 & EEG.xmax >=0
        disp('No EEG.event structure found: creating events of type ''TLE'' (Time-Locking Event) at time 0');
        EEG.event(EEG.trials).epoch = EEG.trials;
        for trial = 1:EEG.trials
            EEG.event(trial).epoch   = trial;
            EEG.event(trial).type    = 'TLE';
            EEG.event(trial).latency = -EEG.xmin*EEG.srate+1+(trial-1)*EEG.pnts;
        end;
    else
        disp('Cannot epoch data with no events'); beep;
        return;
    end;
end;
if ~isfield(EEG.event, 'latency'),
    disp( 'Absent latency field in event array/structure: must name one of the fields ''latency''');
     beep; return;
end;
if size(EEG.data,3) > 1
    epochlim = [num2str( round(EEG.xmin*1000,1)) '  '  num2str(round(EEG.xmax*1000,1))];
else
    epochlim = '-1 2';
end
OLDEEG = EEG;

if nargin < 3
	% popup window parameters
	% -----------------------
   promptstr    = { strvcat('Time-locking event type(s) ([]=all):', ...
                    'Select ''Edit > Event values'' to see type values.'), ...
                    'Epoch limits [start, end] in seconds:', ...
                    'Name for the new dataset:', ...
					'Out-of-bounds EEG rejection limits ([min max], []=none):'  };

   cbevent = ['if ~isfield(EEG.event, ''type'')' ...
				   '   errordlg2(''No type field'');' ...
				   'else' ...
                   '   tmpevent = EEG.event;' ...
                   '   if isnumeric(EEG.event(1).type),' ...
				   '        [tmps,tmpstr] = pop_chansel(unique([ tmpevent.type ]));' ...
				   '   else,' ...
                   '        [tmps,tmpstr] = pop_chansel(unique({ tmpevent.type }));' ...
                   '   end;' ...
				   '   if ~isempty(tmps)' ...
				   '       set(findobj(''parent'', gcbf, ''tag'', ''events''), ''string'', tmpstr);' ...
				   '   end;' ...
				   'end;' ...
				   'clear tmps tmpevent tmpv tmpstr tmpfieldnames;' ];

   geometry = { [2 1 0.5] [2 1 0.5] [2 1.5] [2 1 0.5] };
   uilist = { { 'style' 'text'       'string' 'Time-locking event type(s) ([]=all)' } ...
              { 'style' 'edit'       'string' '' 'tag' 'events' } ...
              { 'style' 'pushbutton' 'string' '...' 'callback' cbevent } ...
              { 'style' 'text'       'string' 'Epoch limits [start, end] in seconds' } ...
              { 'style' 'edit'       'string' epochlim } ...
              { } ...
              { 'style' 'text'       'string' 'Name for the new dataset' } ...
              { 'style' 'edit'       'string'  fastif(isempty(EEG.setname), '', [ EEG.setname ' epochs' ]) } ...
              { 'style' 'text'       'string' 'Out-of-bounds EEG limits if any [min max]' } ...
              { 'style' 'edit'       'string' '' } { } };

   result = inputgui( geometry, uilist, 'pophelp(''pop_epoch'')', 'Extract data epochs - pop_epoch()');
   if length(result) == 0 return; end;

   if strcmpi(result{1}, '[]'), result{1} = ''; end;
   if ~isempty(result{1})
       if strcmpi(result{1}(1),'''')   % If event type appears to be in single-quotes, use comma
                                       % and single-quote as delimiter between event types. toby 2.24.2006
                                       % fixed Arnaud May 2006
            events = eval( [ '{' result{1} '}' ] );
       else events = parsetxt( result{1});
       end;
   else events = {};
   end
   lim = eval( [ '[' result{2} ']' ] );

   args = {};
   if ~isempty( result{3} ),  args = { args{:}, 'newname', result{3} }; end;
   if ~isempty( result{4} ),  args = { args{:}, 'valuelim', eval( [ '[' result{4} ']' ] ) }; end;
   args = { args{:}, 'epochinfo', 'yes' };

else % no interactive inputs
    args = varargin;
end;

% create structure
% ----------------
if ~isempty(args)
   try, g = struct(args{:});
   catch, disp('pop_epoch(): wrong syntax in function arguments'); return; end;
else
    g = [];
end;

% test the presence of variables
% ------------------------------
try, g.epochfield; 	 	  catch, g.epochfield = 'type'; end; % obsolete
try, g.timeunit; 	 	  catch, g.timeunit = 'points'; end;
try, g.verbose; 	      catch, g.verbose = 'on'; end;
try, g.newname; 	      catch, g.newname = fastif(isempty(EEG.setname), '', [EEG.setname ' epochs' ]); end;
try, g.eventindices;      catch, g.eventindices = 1:length(EEG.event); end;
try, g.epochinfo;         catch, g.epochinfo = 'yes'; end;
try, if isempty(g.valuelim), g.valuelim = [-Inf Inf]; end; catch, g.valuelim = [-Inf Inf]; end;

% transform string events into a int array of column indices
% ----------------------------------------------------------
tmpevent = EEG.event;
tmpeventlatency = [ tmpevent(:).latency ];
[tmpeventlatency Itmp] = sort(tmpeventlatency);
EEG.event = EEG.event(Itmp);  % sort by ascending time
Ievent = g.eventindices;

if ~isempty( events )
    % select the events for epoching
    % ------------------------------
    Ieventtmp = [];
    tmpevent = EEG.event;
    tmpeventtype  = { tmpevent.type };
    if iscell(events)
		if isstr(EEG.event(1).type)
			for index2 = 1:length( events )
				tmpevent = events{index2};
				if ~isstr( tmpevent ), tmpevent = num2str( tmpevent ); end;
				Ieventtmp = [ Ieventtmp ; strmatch(tmpevent, tmpeventtype, 'exact') ];
			end;
		else
			for index2 = 1:length( events )
				tmpevent = events{index2};
				if isstr( tmpevent ),tmpevent = str2num( tmpevent ); end;
				if isempty( tmpevent ), error('pop_epoch(): string entered in a numeric field'); end;
				Ieventtmp = [ Ieventtmp find(tmpevent == [ tmpeventtype{:} ]) ];
			end;
		end;
    else
        error('pop_epoch(): multiple event types must be entered as {''a'', ''cell'', ''array''}'); return;
    end;
    Ievent = sort(intersect(Ievent, Ieventtmp));
end;

% select event latencies for epoching
%------------------------------------
Ievent = sort(Ievent);
alllatencies = tmpeventlatency(Ievent);

if isempty(alllatencies)
   error('pop_epoch(): empty epoch range (no epochs were found).'); return;
end;
fprintf('pop_epoch():%d epochs selected\n', length(alllatencies));

try
    % ----------------------------------------------------
    % For AMICA probabilities...Temporarily add model probabilities as channels
    %-----------------------------------------------------
    if isfield(EEG.etc, 'amica') && ~isempty(EEG.etc.amica) && isfield(EEG.etc.amica, 'v_smooth') && ~isempty(EEG.etc.amica.v_smooth) && ~isfield(EEG.etc.amica,'prob_added')
        if isfield(EEG.etc.amica, 'num_models') && ~isempty(EEG.etc.amica.num_models)
            if size(EEG.data,2) == size(EEG.etc.amica.v_smooth,2) && size(EEG.data,3) == size(EEG.etc.amica.v_smooth,3) && size(EEG.etc.amica.v_smooth,1) == EEG.etc.amica.num_models

                EEG = eeg_formatamica(EEG);
                %--------------------
                [EEG indices com] = pop_epoch(EEG,events,lim,args{:});
                %---------------------------------

                EEG = eeg_reformatamica(EEG);
                EEG = eeg_checkamica(EEG);
                return;
            else
                disp('AMICA probabilities not compatible with size of data, model probabilities cannot be epoched...')

            end
        end
    end
    % ----------------------------------------------------
catch
    warnmsg = strcat('your dataset contains amica information, but the amica plugin is not installed.  Continuing and ignoring amica information.');
    warning(warnmsg)
end

% change boundaries in rare cases when limits do not include time-locking events
% ------------------------------------------------------------------------------
tmpevents = EEG.event;
if lim(1) > 0 && ischar(EEG.event(1).type)
   % go through all onset latencies
   for Z1 = length(alllatencies):-1:1
      % if there is any event in between trigger and epoch onset which are boundary events
      selEvt = find([tmpevents.latency] > alllatencies(Z1) & [tmpevents.latency] < alllatencies(Z1) + lim(1) * EEG.srate);
      selEvt = selEvt(strcmp({tmpevents(selEvt).type}, 'boundary'));
      if any(selEvt)
          if sum([tmpevents(selEvt).duration]) > lim(1) * EEG.srate
              alllatencies(Z1) = [];
          else
              % correct the latencies by the duration of the data that were cutout
              alllatencies(Z1) = alllatencies(Z1) - sum([tmpevents(selEvt).duration]);
          end;
      end
   end
end
if lim(2) < 0 && ischar(EEG.event(1).type)
   % go through all onset latencies
   for Z1 = length(alllatencies):-1:1
      % if there is any event in between trigger and epoch onset which are boundary events
      selEvt = find([tmpevents.latency] < alllatencies(Z1) & [tmpevents.latency] > alllatencies(Z1) + lim(2) * EEG.srate);
      selEvt = selEvt(strcmp({tmpevents(selEvt).type}, 'boundary'));
      if any(selEvt)
          if sum([tmpevents(selEvt).duration]) > -lim(2) * EEG.srate
              alllatencies(Z1) = [];
          else
              % correct the latencies by the duration of the data that were cutout
              alllatencies(Z1) = alllatencies(Z1) + sum([tmpevents(selEvt).duration]);
          end;
      end
   end
end

% select event time format and epoch
% ----------------------------------
switch lower( g.timeunit )
 case 'points',	[EEG.data tmptime indices epochevent]= epoch(EEG.data, alllatencies, [lim(1) lim(2)]*EEG.srate, ...
                                                    'valuelim', g.valuelim, 'allevents', tmpeventlatency);
  tmptime = tmptime/EEG.srate;
 case 'seconds',	[EEG.data tmptime indices epochevent]= epoch(EEG.data, alllatencies, lim, 'valuelim', g.valuelim, ...
                                                    'srate', EEG.srate, 'allevents', tmpeventlatency);
	otherwise, disp('pop_epoch(): invalid event time format'); beep; return;
end;
alllatencies = alllatencies(indices);
fprintf('pop_epoch():%d epochs generated\n', length(indices));

% update other fields
% -------------------
if lim(1) ~= tmptime(1) & lim(2)-1/EEG.srate ~= tmptime(2)
	fprintf('pop_epoch(): time limits have been adjusted to [%3.3f %3.3f] to fit data points limits\n', ...
		tmptime(1), tmptime(2)+1/EEG.srate);
end;
EEG.xmin = tmptime(1);
EEG.xmax = tmptime(2);
EEG.pnts = size(EEG.data,2);
EEG.trials = size(EEG.data,3);
EEG.icaact = [];
if ~isempty(EEG.setname)
	if ~isempty(EEG.comments)
		EEG.comments = strvcat(['Parent dataset "' EEG.setname '": ----------'], EEG.comments);
	end;
	EEG.comments = strvcat(['Parent dataset: ' EEG.setname ], ' ', EEG.comments);
end;
EEG.setname = g.newname;

% count the number of events to duplicate and duplicate them
% ----------------------------------------------------------
totlen = 0;
for index=1:EEG.trials, totlen = totlen + length(epochevent{index}); end;
EEG.event(1).epoch = 0;          % create the epoch field (for assignment consistency afterwards)
if totlen ~= 0
    newevent(totlen) = EEG.event(1); % reserve array
else
    newevent = [];
end;

% modify the event structure accordingly (latencies and add epoch field)
% ----------------------------------------------------------------------
allevents = [];
count = 1;
for index=1:EEG.trials
    for indexevent = epochevent{index}
		newevent(count)         = EEG.event(indexevent);
        newevent(count).epoch   = index;
	    newevent(count).latency = newevent(count).latency ...
			- alllatencies(index) - tmptime(1)*EEG.srate + 1 + EEG.pnts*(index-1);
		count = count + 1;
    end;
end;
EEG.event = newevent;
EEG.epoch = [];
EEG = eeg_checkset(EEG, 'eventconsistency');

% check for boundary events
% -------------------------
disp('pop_epoch(): checking epochs for data discontinuity');
if ~isempty(EEG.event) & isstr(EEG.event(1).type)
    tmpevent = EEG.event;
	boundaryindex = strmatch('boundary', { tmpevent.type });
	if ~isempty(boundaryindex)
		indexepoch = [];
		for tmpindex = boundaryindex
            if isfield(tmpevent, 'epoch')
    			indexepoch = [indexepoch tmpevent(tmpindex).epoch ];
            else
                indexepoch = 1; % only one epoch
            end;
		end;
		EEG = pop_select(EEG, 'notrial', indexepoch);
        % update the "indices of accepted events", too
        indices = indices(setdiff(1:length(indices),indexepoch));
	end;
end;

% generate text command
% ---------------------
com = sprintf('%s = pop_epoch( %s, { ', inputname(1), inputname(1));
for j=1:length(events);
    if isstr( events{j} )   com = sprintf('%s ''%s'' ', com, events{j} );
    else                    com = sprintf('%s [%s] ',   com, num2str(events{j}) );
    end;
end;
com = sprintf('%s }, [%s]', com, num2str(lim));
for i=1:2:length(args)
    if ~isempty( args{i+1} )
        if isstr( args{i+1} )   com = sprintf('%s, ''%s'', ''%s''', com, args{i}, args{i+1} );
        else                    com = sprintf('%s, ''%s'', [%s]', com, args{i}, num2str(args{i+1}) );
        end;
    end;
end;
com = [com ');'];
return; % text command

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


pose = icp0301(7300:end,:);
plot(pose(:,4), pose(:,8), 'k.')


text:
function outstr = midrad(c,name,restricted)
%MIDRAD       Display of intervals by midpoint/radius (rigorous)
%
%The call
%
%  midrad(c)        displays interval c in < mid , rad > representation
%  str = midrad(c)  puts output into string str
%
%Output in string str columnwise; can be used for input by intval(str).
%

%for internal use:
%  name                name of output variable
%  restricted == 1     no header, no extra lines output
%
%Call only with 1 or 3 input arguments
%
%Special call:  outstr = midrad(x,[],[])
%  for column vector x output in outstr.exp and outstr.str
%

% written  11/30/98     S.M. Rump
% modified 06/10/98     S.M. Rump  multi-dimensional arrays, exceptions
% modified 06/23/99     S.M. Rump  dble2str -> dble2str_rnd in \private,
%                                  sparse arrays
% modified 08/29/00     S.M. Rump  special output in outstr.exp and .str added
%                                  rounding unchanged after use
% modified 10/03/02     S.M. Rump  improved for Matlab 6.1f and sparse output
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    take care of Matlab sparse Inf/NaN bug
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 11/20/05     S.M. Rump  fast check for rounding to nearest
% modified 02/11/06     S.M. Rump  SparseInfNanFlag removed
%                                    linear indices for huge matrices
% modified 08/25/07     S.M. Rump  huge indices for sparse matrices
% modified 10/14/08     S.M. Rump  huge indices for sparse matrices
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/04/12     S.M. Rump  take care of sparse input
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if nargin<=2
    name = inputname(1);
    restricted = 0;
  end

  loose = strcmp(get(0,'FormatSpacing'),'loose');

  siz = size(c);
  if ~restricted & ~nargout & ( length(siz)<3 )
    if loose, disp(' '); end
    disp([ 'intval ' name ' = ' ])
    if loose, disp(' '); end
  end

  if isempty(c)          % empty interval
    if c.complex
      c.mid
    else
      c.inf
    end
    if loose & ~restricted & ~nargout
      disp(' ')
    end
    setround(rndold)
    return
  end

  if length(siz)>2 & ~nargout
    siz = siz(3:end);
    len = length(siz);
    prodsiz = cumprod([1 siz(1:len-1)]);

    for k=0:prod(siz)-1
      y = mod(floor(k./prodsiz),siz);
      str = sprintf(',%d',y+1);
%VVVV eval(['cs = c(:,:' str ');']);
      s.type = '()';
      eval(['s.subs = {'':'','':''' str '};']);
      cs = subsref(c,s);
%AAAA Matlab V5.2 bug fix
      midrad(cs,[name '(:,:' str ')'],0);
    end
    setround(rndold)
    return
  end

  cmid = mid(c);
  crad = rad(c);
  [m n] = size(cmid);

  format = get(0,'format');
  if ~isequal(format,'short') & ~isequal(format,'long') & ...
     ~isequal(format,'shortE') & ~isequal(format,'longE')
    format = 'short';
  end

  INTLAB_INTVAL_POWER10 = getappdata(0,'INTLAB_INTVAL_POWER10');
  commonexp = 0;
  if ~isequal(format(end),'E')
    % calculate exponent range and common factor
    if c.complex
      cre = abs(real(cmid));
      cim = abs(imag(cmid));
      cre = nonzeros(cre);
      cim = nonzeros(cim);
      if isempty(cre)                   % careful: max(1,[]) = []
        cre = 0;
      end
      if isempty(cim)
        cim = 0;
      end
      cre(isinf(cre)) = 0;
      cim(isinf(cim)) = 0;
      cremax = max(cre(:));
      cimmax = max(cim(:));
      cmax = max(cremax,cimmax);
    else
      cmid_ = abs(cmid);
      cmid_ = nonzeros(cmid_);
      if isempty(cmid_)                 % careful: max(1,[]) = []
        cmid_ = 0;
      end
      cmid_(isinf(cmid_)) = 0;
      cmax = max(cmid_);
      cmax = max(cmax(:));
    end
    crad_ = abs(crad);
    crad_ = nonzeros(crad_);
      if isempty(crad_)                 % careful: max(1,[]) = []
        crad_ = 0;
      end
    crad_(isinf(crad_)) = 0;
    cmax = max( cmax , max(crad_(:)) );
    % care for sparse matrices
    cmax = full(cmax);
    if isempty(cmax)
      cmax = 0;
    end

    emax = floor(log10( cmax + (cmax==0) ));
    if c.complex
      if emax >= 2
        commonexp = emax;
      end
      if emax <= -3
        commonexp = emax+1;
      end
    else
      if isequal(format,'short') & emax >= 3
        commonexp = emax;
      elseif isequal(format,'long') & emax >= 2
        commonexp = emax;
      end
      if emax <= -4
        commonexp = emax+1;
      end
    end
  end
  factor = INTLAB_INTVAL_POWER10.sup(1,commonexp+341);    % 10^commonexp

  if commonexp~=0                      % common factor
    if nargout                         % output to str: be sure with exponent
      if ~( isempty(name) & isempty(restricted) )   % if not special output
        commonexp = 0;
        if ~isequal(format(end),'E')
          format = [ format 'E' ];
        end
      end
    else
      strexp = sprintf('%04d',commonexp);
      if commonexp >= 0
        strexp(1) = '+';
      end
      disp([ '  1.0e' strexp ' *' ])
      if loose, disp(' '); end
    end
  end

  if c.complex
    switch format
      case 'short',  len = 9;  prec = 4;  expon = 0;
      case 'long',   len = 19; prec = 14; expon = 0;
      case 'shortE', len = 13; prec = 4;  expon = 1;
      case 'longE',  len = 24; prec = 15; expon = 1;
    end
    len1 = 3*len+5;    % length of one element in current format
  else
    switch format
      case 'short',  len = 10; prec = 4;  expon = 0;
      case 'long',   len = 19; prec = 14; expon = 0;
      case 'shortE', len = 13; prec = 4;  expon = 1;
      case 'longE',  len = 24; prec = 15; expon = 1;
    end
    len1 = 2*len+2;    % length of one element in current format
  end
  INTLAB_DISPLAY_WIDTH = getappdata(0,'INTLAB_DISPLAY_WIDTH');    % minimum 110, so columns>=1
  columns = floor((INTLAB_DISPLAY_WIDTH+1)/(len1+1));

  formatstr = [ '%' int2str(len) '.' int2str(prec) 'f' ];
  if expon
    formatstr(end) = 'e';
  end

  if nargout & isempty(name) & isempty(restricted)
    % special case: input c column vector, output in outstr.exp and outstr.str

    if commonexp~=0                      % common factor
      strexp = sprintf('%04d',commonexp);
      if commonexp >= 0
        strexp(1) = '+';
      end
      outstr.exp = [ '  1.0e' strexp ' *' ];
    else
      outstr.exp = '';
    end

    if c.complex           % complex intervals
      outstr.str = char(zeros(m,3*len+5));
      for i=1:m
        s = '';
        cmidi = full(cmid(i));
        strrealmid = sprintf(formatstr,real(cmidi)/factor);
        signimagmid = sign(imag(cmidi));
        strimagmid = sprintf(formatstr,abs(imag(cmidi))/factor);
        strimagmid = strimagmid(2:end);     % omit leading space
        [cmidrealinf,cmidrealsup] = str2intval(strrealmid,commonexp);
        if signimagmid==1
          [cmidimaginf,cmidimagsup] = str2intval(strimagmid,commonexp);
        else
          [cmidimagsup,cmidimaginf] = str2intval(strimagmid,commonexp);
          cmidimagsup = - cmidimagsup ;
          cmidimaginf = - cmidimaginf ;
        end
        setround(1)
        re = max( real(cmidi)-cmidrealinf,cmidrealsup-real(cmidi) );
        im = max( imag(cmidi)-cmidimaginf,cmidimagsup-imag(cmidi) );
        if isequal(crad,0)
          cradij = abs( re + sqrt(-1)*im );
        else
          cradij = crad(i) + abs( re + sqrt(-1)*im );
        end
        cradij = full(cradij);
        chsign = '+';
        if signimagmid==-1
          chsign = '-';
        end
        strrad = dble2str_rnd(cradij,commonexp,len-1,prec,expon,+1);
        s = [ s '<' strrealmid ' ' chsign strimagmid 'i,' ...
              strrad '> ' ];
        outstr.str(i,:) = s;
      end
    else                   % real intervals
      outstr.str = char(zeros(m,2*len+3));
      for i=1:m
        s = '';
        cmidi = full(cmid(i));
        strmid = sprintf(formatstr,cmidi/factor);
        [cmidinf,cmidsup] = str2intval(strmid,commonexp);
        setround(1)
        cradij = full( crad(i) + max(cmidi-cmidinf,cmidsup-cmidi) );
        s = [ s '<' strmid ',' ...
              dble2str_rnd(cradij,commonexp,len-1,prec,expon,+1) '> ' ];
        outstr.str(i,:) = s;
      end
    end
    setround(rndold)
    return
  end

  if nargout
    outstr = [];
    if c.complex
      for i=1:prod(size(cmid))
        cmidi = full(cmid(i));
        strrealmid = sprintf(formatstr,real(cmidi));
        strimagmid = sprintf(formatstr,imag(cmidi));
        strimagmid(isspace(strimagmid)) = '';
        if strimagmid(1)~='-'
          strimagmid = [ '+' strimagmid ];
        end
        [cmidrealinf,cmidrealsup] = str2intval(strrealmid,0);
        [cmidimaginf,cmidimagsup] = str2intval(strimagmid,0);
        setround(1)
        re = max( real(cmidi)-cmidrealinf,cmidrealsup-real(cmidi) );
        im = max( imag(cmidi)-cmidimaginf,cmidimagsup-imag(cmidi) );
        if isequal(crad,0)
          cradi = abs( re + sqrt(-1)*im );
        else
          cradi = crad(i) + abs( re + sqrt(-1)*im );
        end
        cradi = full(cradi);
        strrad = dble2str_rnd(cradi,commonexp,len-1,prec,expon,+1);
        outstr = [ outstr '<' strrealmid strimagmid 'i,' strrad '> ' ];
      end
    else
      for i=1:prod(size(cmid))
        cmidi = full(cmid(i));
        strmid = sprintf(formatstr,cmidi);
        [cmidinf,cmidsup] = str2intval(strmid,0);
        setround(1)
        cradi = full( crad(i) + max(cmidi-cmidinf,cmidsup-cmidi) );
        outstr = [ outstr '<' strmid ',' ...
              dble2str_rnd(cradi,commonexp,len-1,prec,expon,+1) '> ' ];
      end
    end
    setround(rndold)
    return
  end

  if issparse(c)
    if isequal(crad,0)                            % mid or rad maybe zero
      [I,J] = find(spones(cmid));
    else
      [I,J] = find(spones(cmid)+spones(crad));
    end
    if length(I)==0
      mid(c)
      setround(rndold)
      return
    end
    if c.complex
      for i=1:length(I)
        str = sprintf('  (%d,%d)',I(i),J(i));
        str = [ str blanks(20-length(str)) ];
        cmidij = full(cmid(I(i),J(i)));
        strrealmid = sprintf(formatstr,real(cmidij)/factor);
        signimagmid = sign(imag(cmidij));
        strimagmid = sprintf(formatstr,abs(imag(cmidij))/factor);
        strimagmid = strimagmid(2:end);     % omit leading space
        [cmidrealinf,cmidrealsup] = str2intval(strrealmid,commonexp);
        if signimagmid==1
          [cmidimaginf,cmidimagsup] = str2intval(strimagmid,commonexp);
        else
          [cmidimagsup,cmidimaginf] = str2intval(strimagmid,commonexp);
          cmidimagsup = - cmidimagsup ;
          cmidimaginf = - cmidimaginf ;
        end
        setround(1)
        re = max( real(cmidij)-cmidrealinf,cmidrealsup-real(cmidij) );
        im = max( imag(cmidij)-cmidimaginf,cmidimagsup-imag(cmidij) );
        if isequal(crad,0)
          cradij = abs( re + sqrt(-1)*im );
        else
          cradij = crad(I(i),J(i)) + abs( re + sqrt(-1)*im );
        end
        cradij = full(cradij);
        chsign = '+';
        if signimagmid==-1
          chsign = '-';
        end
        strrad = dble2str_rnd(cradij,commonexp,len-1,prec,expon,+1);
        str = [ str '<' strrealmid ' ' chsign strimagmid 'i,' ...
                strrad '> ' ];
        disp( str )
      end
    else
      for i=1:length(I)
        str = sprintf('  (%d,%d)',I(i),J(i));
        str = [ str blanks(20-length(str)) ];
        cmidij = full(cmid(I(i),J(i)));
        strmid = sprintf(formatstr,cmidij/factor);
        [cmidinf,cmidsup] = str2intval(strmid,commonexp);
        setround(1)
        cradij = full(crad(I(i),J(i)) + max(cmidij-cmidinf,cmidsup-cmidij));
        str = [ str '<' strmid ',' ...
                dble2str_rnd(cradij,commonexp,len-1,prec,expon,+1) '> ' ];
        disp( str )
      end
    end
    setround(rndold)
    return
  end

  for jj=1:ceil(n/columns)
    j1 = (jj-1)*columns+1;
    if jj*columns<n
      j2 = jj*columns;
    else
      j2 = n;
    end
    if n>columns
      if j1~=j2
        disp(['  Columns ' sprintf('%d',j1) ' through ' sprintf('%d',j2)]);
      else
        disp(['  Column ' sprintf('%d',j1)]);
      end
      if loose, disp(' '); end
    end
    if c.complex           % complex intervals
      for i=1:m
        s = '';
        for j = j1:j2
          cmidij = full(cmid(i,j));
          strrealmid = sprintf(formatstr,real(cmidij)/factor);
          signimagmid = sign(imag(cmidij));
          strimagmid = sprintf(formatstr,abs(imag(cmidij))/factor);
          strimagmid = strimagmid(2:end);     % omit leading space
          [cmidrealinf,cmidrealsup] = str2intval(strrealmid,commonexp);
          if signimagmid==1
            [cmidimaginf,cmidimagsup] = str2intval(strimagmid,commonexp);
          else
            [cmidimagsup,cmidimaginf] = str2intval(strimagmid,commonexp);
            cmidimagsup = - cmidimagsup ;
            cmidimaginf = - cmidimaginf ;
          end
          setround(1)
          re = max( real(cmidij)-cmidrealinf,cmidrealsup-real(cmidij) );
          im = max( imag(cmidij)-cmidimaginf,cmidimagsup-imag(cmidij) );
          if isequal(crad,0)
            cradij = abs( re + sqrt(-1)*im );
          else
            cradij = crad(i,j) + abs( re + sqrt(-1)*im );
          end
          cradij = full(cradij);
          chsign = '+';
          if signimagmid==-1
            chsign = '-';
          end
          strrad = dble2str_rnd(cradij,commonexp,len-1,prec,expon,+1);
          s = [ s '<' strrealmid ' ' chsign strimagmid 'i,' ...
                strrad '> ' ];
        end
        disp(s)
      end
    else                   % real intervals
      for i=1:m
        s = '';
        for j = j1:j2
          cmidij = full(cmid(i,j));
          strmid = sprintf(formatstr,cmidij/factor);
          [cmidinf,cmidsup] = str2intval(strmid,commonexp);
          setround(1)
          cradij = full( crad(i,j) + max(cmidij-cmidinf,cmidsup-cmidij) );
          s = [ s '<' strmid ',' ...
                dble2str_rnd(cradij,commonexp,len-1,prec,expon,+1) '> ' ];
        end
        disp(s)
      end
    end
    if loose & ~restricted, disp(' '); end
  end

  setround(rndold)

function [xinf,xsup] = str2intval(str,commonexp)
%STR2INTVAL   internal routine for output routines midrad
%

%Rigorous conversion of string to lower/upper bound with directed rounding
%  and common exponent; only for one real input.
%Input string in e- or f-format, possibly with exponent
%Input number is assumed to be correct.
%

  % exception handling
  if any( str=='a' ) | any( str=='f' )    % NaN or Inf
    xinf = NaN;
    xsup = NaN;
    return
  end

  INTLAB_INTVAL_POWER10 = getappdata(0,'INTLAB_INTVAL_POWER10');
  [m i] = max(str~=' ');
  str(1:i-1) = '';             % omit leading blanks

  sign = 1;                    % calculate sign
  if str(1)=='+'
    str = str(2:end);
  elseif str(1)=='-'
    sign = -1;
    str = str(2:end);
  end

  indexdp = find(str=='.');    % index of decimal point
  str(indexdp) = '';

  indexexp = find(str=='e');   % check exponent
  if ~isempty(indexexp)
    Exp = ( 2*(str(indexexp+1)=='+') - 1 ) * ...
          (str(indexexp+2:end)-'0')*(10.^((length(str)-indexexp-2):-1:0)');
    indexend = indexexp-1;
  else
    Exp = 0;
    indexend = length(str);
  end

  m = str(indexend:-1:1) - '0';

  % convert back to double
  offset = 9*(indexdp-2+Exp+commonexp+341);
  mant = 9*(indexend:-1:1);
  index = (m~=0).*(offset - mant + m) + (m==0);
  if sign==1
    setround(-1)
    xinf = sum( INTLAB_INTVAL_POWER10.inf(index) );
    setround(1)
    xsup = sum( INTLAB_INTVAL_POWER10.sup(index) );
  else
    setround(1)
    xinf = - sum( INTLAB_INTVAL_POWER10.sup(index) );
    setround(-1)
    xsup = - sum( INTLAB_INTVAL_POWER10.inf(index) );
  end
  setround(0)

repo: OpenMAS


text:
%quad_diff   solve Poisson problem in quadrilateral domain
%   IFISS scriptfile: DJS; 4 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
clear
%% define geometry
pde=1; domain=8;
global viscosity
viscosity=1;
quad_domain
load quad_grid
%
%% set up matrices
qmethod=default('Q1/Q2 approximation 1/2? (default Q1)',1);
if qmethod ==2,
   [x,y,xy] = q2grid(x,y,xy,mv,bound);
   [A,M,f] = femq2_diff(xy,mv);
else
   [ev,ebound] = q1grid(xy,mv,bound,mbound);
   [A,M,f] = femq1_diff(xy,ev);
end
%
%% boundary conditions
[Agal,fgal] = nonzerobc(A,f,xy,bound);
%
%% save resulting system
fprintf('system saved in quad_diff.mat ...\n')
gohome
cd datafiles
save quad_diff qmethod Agal M  fgal  xy x y
%% compute solution
fprintf('solving linear system ...  ')
x_gal=Agal\fgal;
fprintf('done\n')
save quad_diff x_gal  -append
% plot solution
if qmethod > 0
   solplot(x_gal,xy,x,y,10);
   title(['Q',int2str(qmethod),' solution'])
   drawnow
end

author: johannesgerer

repo: jburkardt-m


text:
% Program for Newton-Raphson Load Flow with STATCOM

Y = ybusppg();              % Calling ybusppg.m..
busdata = busdata30();       % Calling Busdata..
statdata = statdata30();     % Statcom Data..
baseMVA = 100;              % Base MVA..
bus = busdata(:,1);         % Bus Number..
type = busdata(:,2);        % Type of Bus 1-Slack, 2-PV, 3-PQ..
V = busdata(:,3);           % Specified Voltage..
del = busdata(:,4);         % Voltage Angle..
Pg = busdata(:,5);          % PGi..
Qg = busdata(:,6);          % QGi..
Pl = busdata(:,7);          % PLi..
Ql = busdata(:,8);          % QLi..
Qmin = busdata(:,9);        % Minimum Reactive Power Limit..
Qmax = busdata(:,10);       % Maximum Reactive Power Limit..
nbus = max(bus);            % To get no. of buses..
P = Pg - Pl;                % Pi = PGi - PLi..
Q = Qg - Ql;                % Qi = QGi - QLi..
P = P/baseMVA;              % Converting to p.u..
Q = Q/baseMVA;
Qmin = Qmin/baseMVA;
Qmax = Qmax/baseMVA;
Tol = 1;
Iter = 1;
Psp = P;
Qsp = Q;
G = real(Y);    % Conductance..
B = imag(Y);    % Susceptance..
Vsp = V;

% Details of STATCOM
statb = statdata(:,1); % Buses at which statcoms are placed..
Vsh = statdata(:,2);
Thst = statdata(:,3);
Qsmx = statdata(:,4);
Qsmn = statdata(:,5);
gsh = 0.9901;
bsh = -9.901;
Vshmx = 1.1; Vshmn = 0.9;
Thstmx = pi; Thstmn = -pi;
np = length(statb); % Number of STATCOMs..

pv = find(type == 2 | type == 1); % Index of PV Buses..
pq = find(type == 3); % Index of PQ Buses..

npv = length(pv); % Number of PV buses..
npq = length(pq); % Number of PQ buses..

while (Tol > 1e-5 && Iter <= 50)   % Iteration starting..

    P = zeros(nbus,1);
    Q = zeros(nbus,1);
    % Calculate P and Q
    for i = 1:nbus
        for k = 1:nbus
            P(i) = P(i) + V(i)* V(k)*(G(i,k)*cos(del(i)-del(k)) + B(i,k)*sin(del(i)-del(k)));
            Q(i) = Q(i) + V(i)* V(k)*(G(i,k)*sin(del(i)-del(k)) - B(i,k)*cos(del(i)-del(k)));
        end
        m = find(statb == i);
        if ~isempty(m)
            P(i) = P(i) + V(i)^2*gsh - V(i)*Vsh(m)*(gsh*cos(del(i)-Thst(m)) + bsh*sin(del(i)-Thst(m)));
            Q(i) = Q(i) - V(i)^2*bsh + V(i)*Vsh(m)*(bsh*cos(del(i)-Thst(m)) - gsh*sin(del(i)-Thst(m)));
        end
    end

    % Checking Q-limit violations..
    if Iter >= 2
        for n = 1:nbus
            if type(n) == 2
                if Q(n) < Qmin(n)
                    V(n) = V(n) + 0.01;
                elseif Q(n) > Qmax(n)
                    V(n) = V(n) - 0.01;
                end
            end
         end
    end

    % Calculating PE,
    PE = zeros(np,1);
    for i = 1:np
        m = statb(i);
        PE(i) = Vsh(i)^2*gsh - V(m)*Vsh(i)*(gsh*cos(del(m)-Thst(i)) - bsh*sin(del(m)-Thst(i)));
    end

    % Calculate F, Control Variables
    F = zeros(np,1);
    for i = 1:np
        m = statb(i);
        F(i) = V(m) - Vsp(m);
    end

    % Calculate change from specified value
    dPa = Psp-P;
    dQa = Qsp-Q;
    dQ = zeros(npq,1);
    k = 1;
    for i = 1:nbus
        if type(i) == 3
            dQ(k,1) = dQa(i);
            k = k+1;
        end
    end
    dP = dPa(2:nbus);
    dPE = -PE;
    dF = -F;
    M = [dP; dQ; dPE; dF];       % Mismatch Vector

    % Jacobian
    % J11 - Derivative of Real Power Injections with Angles..
    J11 = zeros(nbus-1,nbus-1);
    for i = 1:(nbus-1)
        m = i+1;
        p = find(statb == m);
        for k = 1:(nbus-1)
            q = k+1;
            if q == m
                for n = 1:nbus
                    J11(i,k) = J11(i,k) + V(m)* V(n)*(-G(m,n)*sin(del(m)-del(n)) + B(m,n)*cos(del(m)-del(n)));
                end
                J11(i,k) = J11(i,k) - V(m)^2*B(m,m);
                if ~isempty(p)
                    J11(i,k) = J11(i,k) + V(m)*Vsh(p)*(gsh*sin(del(m)-Thst(p)) - bsh*cos(del(m)-Thst(p)));
                end
            else
                J11(i,k) = V(m)* V(q)*(G(m,q)*sin(del(m)-del(q)) - B(m,q)*cos(del(m)-del(q)));
            end
        end
    end

    % J12 - Derivative of Real Power Injections with V..
    J12 = zeros(nbus-1,npq);
    for i = 1:(nbus-1)
        m = i+1;
        p = find(statb == m);
        for k = 1:npq
            q = pq(k);
            if q == m
                for n = 1:nbus
                    J12(i,k) = J12(i,k) + V(n)*(G(m,n)*cos(del(m)-del(n)) + B(m,n)*sin(del(m)-del(n)));
                end
                J12(i,k) = J12(i,k) + V(m)*G(m,m);
                if ~isempty(p)
                    J12(i,k) = J12(i,k) + 2*V(m)*gsh - Vsh(p)*(gsh*cos(del(m)-Thst(p)) + bsh*sin(del(m)-Thst(p)));
                end
            else
                J12(i,k) = V(m)*(G(m,q)*cos(del(m)-del(q)) + B(m,q)*sin(del(m)-del(q)));
            end
        end
    end

    % J13 - Derivative of Real Power Injections with Vsh..
    % J14 - Derivative of Real Power Injections with Thsh..
    J13 = zeros(nbus-1,np);
    J14 = zeros(nbus-1,np);
    for i = 1:(nbus-1)
        m = i+1;
        for k = 1:np
            p = statb(k);
            if m == p
                J13(i,k) = -V(m)*(gsh*cos(del(m)-Thst(k)) + bsh*sin(del(m)-Thst(k)));
                J14(i,k) = -V(m)*Vsh(k)*(gsh*sin(del(m)-Thst(k)) - bsh*cos(del(m)-Thst(k)));
            end
        end
    end

    % J21 - Derivative of Reactive Power Injections with Angles..
    J21 = zeros(npq,nbus-1);
    for i = 1:npq
        m = pq(i);
        p = find(statb == m);
        for k = 1:(nbus-1)
            q = k+1;
            if q == m
                for n = 1:nbus
                    J21(i,k) = J21(i,k) + V(m)* V(n)*(G(m,n)*cos(del(m)-del(n)) + B(m,n)*sin(del(m)-del(n)));
                end
                J21(i,k) = J21(i,k) - V(m)^2*G(m,m);
                if ~isempty(p)
                    J21(i,k) = J21(i,k) - V(m)*Vsh(p)*(gsh*cos(del(m)-Thst(p)) + bsh*sin(del(m)-Thst(p)));
                end
            else
                J21(i,k) = -V(m)* V(q)*(G(m,q)*cos(del(m)-del(q)) + B(m,q)*sin(del(m)-del(q)));
            end
        end
    end

    % J22 - Derivative of Reactive Power Injections with V..
    J22 = zeros(npq,npq);
    for i = 1:npq
        m = pq(i);
        p = find(statb == m);
        for k = 1:npq
            q = pq(k);
            if q == m
                for n = 1:nbus
                    J22(i,k) = J22(i,k) + V(n)*(G(m,n)*sin(del(m)-del(n)) - B(m,n)*cos(del(m)-del(n)));
                end
                J22(i,k) = J22(i,k) - V(m)*B(m,m);
                if ~isempty(p)
                    J22(i,k) = J22(i,k) - 2*V(m)*bsh - Vsh(p)*(gsh*sin(del(m)-Thst(p)) - bsh*cos(del(m)-Thst(p)));
                end
            else
                J22(i,k) = V(m)*(G(m,q)*sin(del(m)-del(q)) - B(m,q)*cos(del(m)-del(q)));
            end
        end
    end

    % J23 - Derivative of Reactive Power Injections with Vsh..
    % J24 - Derivative of Reactive Power Injections with Thsh..
    J23 = zeros(npq,np);
    J24 = zeros(npq,np);
    for i = 1:npq
        q = pq(i);
        m = i+1;
        for k = 1:np
            p = statb(k);
            if q == p
                J23(i,k) = -V(m)*(gsh*sin(del(m)-Thst(k)) - bsh*cos(del(m)-Thst(k)));
                J24(i,k) = V(m)*Vsh(k)*(gsh*cos(del(m)-Thst(k)) + bsh*sin(del(m)-Thst(k)));
            end
        end
    end

    % J31 - Derivative of PE with Angles..
    % J41 - Derivative of F with Angles..
    J31 = zeros(np,nbus-1);
    J41 = zeros(np,nbus-1);
    for i = 1:np
        m = statb(i);
        for k = 1:(nbus-1)
            if m == k+1
                J31(i,k) = V(m)*Vsh(i)*(gsh*sin(del(m)-Thst(i)) + bsh*cos(del(m)-Thst(i)));
            end
        end
    end

    % J32 - Derivative of PE with V..
    % J42 - Derivative of F with V..
    J32 = zeros(np,npq);
    J42 = zeros(np,npq);
    for i = 1:np
        m = statb(i);
        for k = 1:npq
            if m == pq(k)
                J32(i,k) = -Vsh(i)*(gsh*cos(del(m)-Thst(i)) - bsh*sin(del(m)-Thst(i)));
            end
            if m == pq(k)
                J42(i,k) = 1;
            end
        end
    end

    % J33 - Derivative of PE with Vsh..
    % J34 - Derivative of PE with Thsh..
    % J43 - Derivative of F with Vsh..
    % J44 - Derivative of F with Thsh..
    J33 = zeros(np,np);
    J34 = zeros(np,np);
    J43 = zeros(np,np);
    J44 = zeros(np,np);
    for i = 1:np
        m = statb(i);
        for k = 1:np
            p = statb(k);
            if m == p
                J33(i,k) = 2*Vsh(k)*gsh - V(m)*(gsh*cos(del(m)-Thst(k)) - bsh*sin(del(m)-Thst(k)));
                J34(i,k) = -V(m)*Vsh(k)*(gsh*sin(del(m)-Thst(k)) + bsh*cos(del(m)-Thst(k)));
            end
        end
    end

    J = [J11 J12 J13 J14; J21 J22 J23 J24; J31 J32 J33 J34; J41 J42 J43 J44];     % Jacobian
    clear J11 J12 J13 J14 J21 J22 J23 J24 J31 J32 J33 J34 J41 J42 J43 J44

    X = inv(J)*M;           % Correction Vector
    dTh = X(1:nbus-1);
    dV = X(nbus:nbus+npq-1);
    dVsh = X(nbus+npq:nbus+npq+np-1);
    dThst = X(nbus+npq+np:nbus+npq+2*np-1);
    del(2:nbus) = dTh + del(2:nbus);
    k = 1;
    for i = 2:nbus
        if type(i) == 3
            V(i) = dV(k) + V(i);
            k = k+1;
        end
    end
    Vsh = Vsh + dVsh;
    Thst = Thst + dThst;

    % Calculate Qsh..
    Qsh = zeros(np,1);
    for m = 1:np
        i = statb(m);
        Qsh(m) = -V(i)^2*bsh + V(i)*Vsh(m)*(bsh*cos(del(i)-Thst(m)) - gsh*sin(del(i)-Thst(m)));
    end
    Iter = Iter + 1;
    Tol = max(abs(M));
end

Iter = Iter - 1; % Number of Iterations took..
V;
Del = 180/pi*del;
Thst = 180/pi*Thst;
E2 = [V Del]; % Bus Voltages and angles..
disp('------------------------------');
disp('|  Bus  |    V    |  Angle   | ');
disp('|  No   |   pu    |  Degree  | ');
disp('------------------------------');
for m = 1:nbus
    fprintf('%4g', m), fprintf('    %8.4f', V(m)), fprintf('    %8.4f', Del(m)),fprintf('\n');
end
disp('-----------------------------');
disp('----------------------------------------');
disp('| STATCOM |  Vsh   |  Thst    |   Qsh  |');
disp('|   Bus   |   pu   |  Degree  |    pu  |');
disp('----------------------------------------');
for m = 1:np
    fprintf('  %4g',statb(m)), fprintf('   %8.4f', Vsh(m)), fprintf('    %8.4f', Thst(m)),fprintf('  %8.4f', Qsh(m)), fprintf('\n');
end
disp('----------------------------------------');

author: Sable

repo: mcbench-benchmarks


text:
%  Imbed Zoom, Interp, and Info menu to view_nii window.
%
%  Usage: view_nii_menu(fig);
%

%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
%--------------------------------------------------------------------
function menu_hdl = view_nii_menu(fig, varargin)

   if isnumeric(fig)
      menu_hdl = init(fig);
      return;
   end

   menu_hdl = [];

   switch fig
   case 'interp'
      if nargin > 1
         fig = varargin{1};
      else
         fig = gcbf;
      end

      nii_menu = getappdata(fig, 'nii_menu');
      interp_on_state = get(nii_menu.Minterp,'Userdata');

      if (interp_on_state == 1)
         opt.useinterp = 1;
         view_nii(fig,opt);
         set(nii_menu.Minterp,'Userdata',0,'Label','Interp off');
         reset_zoom(fig);
      else
         opt.useinterp = 0;
         view_nii(fig,opt);
         set(nii_menu.Minterp,'Userdata',1,'Label','Interp on');
         reset_zoom(fig);
      end
   case 'reset_zoom'
      if nargin > 1
         fig = varargin{1};
      else
         fig = gcbf;
      end

      reset_zoom(fig);
   case 'orient'
      orient;
   case 'editvox'
      editvox;
   case 'img_info'
      img_info;
   case 'save_disp'
      save_disp;
   end

   return					% view_nii_menu

%--------------------------------------------------------------------
function menu_hdl = init(fig)

   %  search for edit, view menu
   %
   nii_menu.Mfile = [];
   nii_menu.Medit = [];
   nii_menu.Mview = [];
   menuitems = findobj(fig, 'type', 'uimenu');

   for i=1:length(menuitems)
      filelabel = get(menuitems(i),'label');

      if strcmpi(strrep(filelabel, '&', ''), 'file')
         nii_menu.Mfile = menuitems(i);
      end

      editlabel = get(menuitems(i),'label');

      if strcmpi(strrep(editlabel, '&', ''), 'edit')
         nii_menu.Medit = menuitems(i);
      end

      viewlabel = get(menuitems(i),'label');

      if strcmpi(strrep(viewlabel, '&', ''), 'view')
         nii_menu.Mview = menuitems(i);
      end
   end

   set(fig, 'menubar', 'none');

   if isempty(nii_menu.Mfile)
      nii_menu.Mfile = uimenu('Parent',fig, ...
   	   'Label','File');

      nii_menu.Mfile_save = uimenu('Parent',nii_menu.Mfile, ...
   	   'Label','Save displayed image as ...', ...
           'Callback','view_nii_menu(''save_disp'');');
   else
      nii_menu.Mfile_save = uimenu('Parent',nii_menu.Mfile, ...
   	   'Label','Save displayed image as ...', ...
           'separator','on', ...
           'Callback','view_nii_menu(''save_disp'');');
   end

   if isempty(nii_menu.Medit)
      nii_menu.Medit = uimenu('Parent',fig, ...
   	   'Label','Edit');

      nii_menu.Medit_orient = uimenu('Parent',nii_menu.Medit, ...
   	   'Label','Convert to RAS orientation', ...
           'Callback','view_nii_menu(''orient'');');

      nii_menu.Medit_editvox = uimenu('Parent',nii_menu.Medit, ...
   	   'Label','Edit voxel value at crosshair', ...
           'Callback','view_nii_menu(''editvox'');');
   else
      nii_menu.Medit_orient = uimenu('Parent',nii_menu.Medit, ...
   	   'Label','Convert to RAS orientation', ...
           'separator','on', ...
           'Callback','view_nii_menu(''orient'');');

      nii_menu.Medit_editvox = uimenu('Parent',nii_menu.Medit, ...
   	   'Label','Edit voxel value at crosshair', ...
           'Callback','view_nii_menu(''editvox'');');
   end

   if isempty(nii_menu.Mview)
      nii_menu.Mview = uimenu('Parent',fig, ...
   	   'Label','View');

      nii_menu.Mview_info = uimenu('Parent',nii_menu.Mview, ...
   	   'Label','Image Information', ...
           'Callback','view_nii_menu(''img_info'');');
   else
      nii_menu.Mview_info = uimenu('Parent',nii_menu.Mview, ...
   	   'Label','Image Information', ...
           'separator','on', ...
           'Callback','view_nii_menu(''img_info'');');
   end

   nii_menu.Mzoom = rri_zoom_menu(fig);

   nii_menu.Minterp = uimenu('Parent',fig, ...
   	   'Label','Interp on', ...
	   'Userdata', 1, ...
           'Callback','view_nii_menu(''interp'');');

   setappdata(fig,'nii_menu',nii_menu);
   menu_hdl = nii_menu.Minterp;

   return					% init

%----------------------------------------------------------------
function reset_zoom(fig)

   old_handle_vis = get(fig, 'HandleVisibility');
   set(fig, 'HandleVisibility', 'on');

   nii_view = getappdata(fig, 'nii_view');
   nii_menu = getappdata(fig, 'nii_menu');

   set(nii_menu.Mzoom,'Userdata',1,'Label','Zoom on');
   set(fig,'pointer','arrow');
   zoom off;

   axes(nii_view.handles.axial_axes);
   setappdata(get(gca,'zlabel'), 'ZOOMAxesData', ...
			[get(gca, 'xlim') get(gca, 'ylim')])
%   zoom reset;
 %  zoom getlimits;
   zoom out;

   axes(nii_view.handles.coronal_axes);
   setappdata(get(gca,'zlabel'), 'ZOOMAxesData', ...
			[get(gca, 'xlim') get(gca, 'ylim')])
%   zoom reset;
 %  zoom getlimits;
   zoom out;

   axes(nii_view.handles.sagittal_axes);
   setappdata(get(gca,'zlabel'), 'ZOOMAxesData', ...
			[get(gca, 'xlim') get(gca, 'ylim')])
%   zoom reset;
 %  zoom getlimits;
   zoom out;

   set(fig, 'HandleVisibility', old_handle_vis);

   return;					% reset_zoom

%----------------------------------------------------------------
function img_info

   nii_view = getappdata(gcbf, 'nii_view');
   hdr = nii_view.nii.hdr;

   max_value = num2str(double(max(nii_view.nii.img(:))));
   min_value = num2str(double(min(nii_view.nii.img(:))));

   dim = sprintf('%d  %d  %d', double(hdr.dime.dim(2:4)));
   vox = sprintf('%.3f  %.3f  %.3f', double(hdr.dime.pixdim(2:4)));

   if double(hdr.dime.datatype) == 1
      type = '1-bit binary';
   elseif double(hdr.dime.datatype) == 2
      type = '8-bit unsigned integer';
   elseif double(hdr.dime.datatype) == 4
      type = '16-bit signed integer';
   elseif double(hdr.dime.datatype) == 8
      type = '32-bit signed integer';
   elseif double(hdr.dime.datatype) == 16
      type = '32-bit single float';
   elseif double(hdr.dime.datatype) == 64
      type = '64-bit double precision';
   elseif double(hdr.dime.datatype) == 128
      type = '24-bit RGB true color';
   elseif double(hdr.dime.datatype) == 256
      type = '8-bit signed integer';
   elseif double(hdr.dime.datatype) == 511
      type = '96-bit RGB true color';
   elseif double(hdr.dime.datatype) == 512
      type = '16-bit unsigned integer';
   elseif double(hdr.dime.datatype) == 768
      type = '32-bit unsigned integer';
   elseif double(hdr.dime.datatype) == 1024
      type = '64-bit signed integer';
   elseif double(hdr.dime.datatype) == 1280
      type = '64-bit unsigned integer';
   end

   msg = {};
   msg = [msg {''}];
   msg = [msg {['Dimension:  [', dim, ']']}];
   msg = [msg {''}];
   msg = [msg {['Voxel Size:  [', vox, ']']}];
   msg = [msg {''}];
   msg = [msg {['Data Type:  [', type, ']']}];
   msg = [msg {''}];
   msg = [msg {['Max Value:  [', max_value, ']']}];
   msg = [msg {''}];
   msg = [msg {['Min Value:  [', min_value, ']']}];
   msg = [msg {''}];

   if isfield(nii_view.nii, 'fileprefix')
      if isfield(nii_view.nii, 'filetype') & nii_view.nii.filetype == 2
         msg = [msg {['File Name:  [', nii_view.nii.fileprefix, '.nii]']}];
         msg = [msg {''}];
      elseif isfield(nii_view.nii, 'filetype')
         msg = [msg {['File Name:  [', nii_view.nii.fileprefix, '.img]']}];
         msg = [msg {''}];
      else
         msg = [msg {['File Prefix:  [', nii_view.nii.fileprefix, ']']}];
         msg = [msg {''}];
      end
   end

   h = msgbox(msg, 'Image Information', 'modal');
   set(h,'color',[1 1 1]);

   return;					% img_info

%----------------------------------------------------------------
function orient

   fig = gcbf;
   nii_view = getappdata(fig, 'nii_view');
   nii = nii_view.nii;

   if ~isempty(nii_view.bgimg)
      msg = 'You can not modify an overlay image';
      h = msgbox(msg, 'Error', 'modal');
      return;
   end

   old_pointer = get(fig,'Pointer');
   set(fig,'Pointer','watch');

   [nii orient] = rri_orient(nii);

   if isequal(orient, [1 2 3])		% do nothing
      set(fig,'Pointer',old_pointer);
      return;
   end

   oldopt = view_nii(fig);
   opt.command = 'updatenii';
   opt.usecolorbar = oldopt.usecolorbar;
   opt.usepanel = oldopt.usepanel;
   opt.usecrosshair = oldopt.usecrosshair;
   opt.usestretch = oldopt.usestretch;
   opt.useimagesc = oldopt.useimagesc;
   opt.useinterp = oldopt.useinterp;
   opt.setarea = oldopt.area;
   opt.setunit = oldopt.unit;
   opt.setviewpoint = oldopt.viewpoint;
   opt.setscanid = oldopt.scanid;
   opt.setcbarminmax = oldopt.cbarminmax;
   opt.setcolorindex = oldopt.colorindex;
   opt.setcolormap = oldopt.colormap;
   opt.setcolorlevel = oldopt.colorlevel;

   if isfield(oldopt,'highcolor')
      opt.sethighcolor = oldopt.highcolor;
   end

   view_nii(fig, nii, opt);
   set(fig,'Pointer',old_pointer);
   reset_zoom(fig);

   return;					% orient

%----------------------------------------------------------------
function editvox

   fig = gcbf;
   nii_view = getappdata(fig, 'nii_view');

   if ~isempty(nii_view.bgimg)
      msg = 'You can not modify an overlay image';
      h = msgbox(msg, 'Error', 'modal');
      return;
   end

   nii = nii_view.nii;
   oldopt = view_nii(fig);
   sag = nii_view.imgXYZ.vox(1);
   cor = nii_view.imgXYZ.vox(2);
   axi = nii_view.imgXYZ.vox(3);

   if nii_view.nii.hdr.dime.datatype == 128
      imgvalue = [double(nii.img(sag,cor,axi,1,nii_view.scanid)) double(nii.img(sag,cor,axi,2,nii_view.scanid)) double(nii.img(sag,cor,axi,3,nii_view.scanid))];
      init_val = sprintf('%7.4g %7.4g %7.4g',imgvalue);
   elseif nii_view.nii.hdr.dime.datatype == 511
      R = double(nii.img(sag,cor,axi,1,nii_view.scanid)) * (nii_view.nii.hdr.dime.glmax - ...
		nii_view.nii.hdr.dime.glmin) + nii_view.nii.hdr.dime.glmin;
      G = double(nii.img(sag,cor,axi,2,nii_view.scanid)) * (nii_view.nii.hdr.dime.glmax - ...
		nii_view.nii.hdr.dime.glmin) + nii_view.nii.hdr.dime.glmin;
      B = double(nii.img(sag,cor,axi,3,nii_view.scanid)) * (nii_view.nii.hdr.dime.glmax - ...
		nii_view.nii.hdr.dime.glmin) + nii_view.nii.hdr.dime.glmin;
      imgvalue = [R G B];
      init_val = sprintf('%7.4g %7.4g %7.4g',imgvalue);
   else
      imgvalue = double(nii.img(sag,cor,axi,nii_view.scanid));
      init_val = sprintf('%.6g',imgvalue);
   end

   old_pointer = get(fig,'Pointer');
   set(fig,'Pointer','watch');

   repeat = 1;
   while repeat
      if nii_view.nii.hdr.dime.datatype == 128 | nii_view.nii.hdr.dime.datatype == 511
         init_val = inputdlg({'Replace the current voxel values with 3 new numbers:'}, ...
		'Edit voxel value at crosshair', 1, {num2str(init_val)});
      else
         init_val = inputdlg({'Replace the current voxel value with 1 new number:'}, ...
		'Edit voxel value at crosshair', 1, {num2str(init_val)});
      end

      if isempty(init_val)
         set(fig,'Pointer',old_pointer);
         return
      end

      imgvalue = str2num(init_val{1});

      if ( (nii_view.nii.hdr.dime.datatype == 128 | nii_view.nii.hdr.dime.datatype == 511) ...
		& length(imgvalue) ~= 3 ) | ...
         ( (nii_view.nii.hdr.dime.datatype ~= 128 & nii_view.nii.hdr.dime.datatype ~= 511) ...
		& length(imgvalue) ~= 1 )
	% do nothing
      else
         repeat = 0;
      end
   end

   if nii_view.nii.hdr.dime.datatype == 128
      nii.img(sag,cor,axi,1,nii_view.scanid) = imgvalue(1);
      nii.img(sag,cor,axi,2,nii_view.scanid) = imgvalue(2);
      nii.img(sag,cor,axi,3,nii_view.scanid) = imgvalue(3);
   elseif nii_view.nii.hdr.dime.datatype == 511
      nii.img(sag,cor,axi,1,nii_view.scanid) = (imgvalue(1) - nii_view.nii.hdr.dime.glmin) ...
		/ (nii_view.nii.hdr.dime.glmax - nii_view.nii.hdr.dime.glmin);
      nii.img(sag,cor,axi,2,nii_view.scanid) = (imgvalue(2) - nii_view.nii.hdr.dime.glmin) ...
		/ (nii_view.nii.hdr.dime.glmax - nii_view.nii.hdr.dime.glmin);
      nii.img(sag,cor,axi,3,nii_view.scanid) = (imgvalue(3) - nii_view.nii.hdr.dime.glmin) ...
		/ (nii_view.nii.hdr.dime.glmax - nii_view.nii.hdr.dime.glmin);
   else
      nii.img(sag,cor,axi,nii_view.scanid) = imgvalue;
   end

   opt.command = 'updatenii';
   opt.usecolorbar = oldopt.usecolorbar;
   opt.usepanel = oldopt.usepanel;
   opt.usecrosshair = oldopt.usecrosshair;
   opt.usestretch = oldopt.usestretch;
   opt.useimagesc = oldopt.useimagesc;
   opt.useinterp = oldopt.useinterp;
   opt.setarea = oldopt.area;
   opt.setunit = oldopt.unit;
   opt.setviewpoint = oldopt.viewpoint;
   opt.setscanid = oldopt.scanid;
   opt.setcbarminmax = oldopt.cbarminmax;
   opt.setcolorindex = oldopt.colorindex;
   opt.setcolormap = oldopt.colormap;
   opt.setcolorlevel = oldopt.colorlevel;

   if isfield(oldopt,'highcolor')
      opt.sethighcolor = oldopt.highcolor;
   end

   view_nii(fig, nii, opt);
   set(fig,'Pointer',old_pointer);
   reset_zoom(fig);

   return;					% editvox

%----------------------------------------------------------------
function save_disp

   [filename pathname] = uiputfile('*.*', 'Save displayed image as (*.nii or *.img)');

   if isequal(filename,0) | isequal(pathname,0)
      return;
   else
      out_imgfile = fullfile(pathname, filename);	% original image file
   end

   old_pointer = get(gcbf,'Pointer');
   set(gcbf,'Pointer','watch');

   nii_view = getappdata(gcbf, 'nii_view');
   nii = nii_view.nii;

   try
      save_nii(nii, out_imgfile);
   catch
      msg = 'File can not be saved.';
      msgbox(msg, 'File write error', 'modal');
   end

   set(gcbf,'Pointer',old_pointer);

   return;					% save_disp

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function u = unicycle_unrank ( n, rank )

%*****************************************************************************80
%
%% UNICYCLE_UNRANK "unranks" a unicycle.
%
%  Discussion:
%
%    That is, given a rank, it computes the corresponding unicycle.
%
%    The value of the rank should be between 0 and (N-1)%-1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2012
%
%  Author:
%
%    John Burkardt.
%
%  Reference:
%
%    Dennis Stanton, Dennis White,
%    Constructive Combinatorics,
%    Springer, 1986,
%    ISBN: 0387963472,
%    LC: QA164.S79.
%
%  Parameters:
%
%    Input, integer N, the number of elements in the set.
%
%    Input, integer RANK, the desired rank of the permutation.
%
%    Output, integer U(N), the unicycle.
%
  p = perm_lex_unrank ( n - 1, rank );

  u(1) = 1;
  u(2:n) = p(1:n-1) + 1;

  return
end


text:
function userhelp(subject,page)
% USERHELP On-Line Help facility for The Controls Tutor

% Author: Craig Borghesani
% Date: 11/26/94
% Revised:
% Copyright (c) 1999, Prentice-Hall

if subject == 0,
 f = gcf;
else
 f2 = gcf;
 f = get(f2,'userdata');
end
ui_data = get(f,'userdata');
ui_han = ui_data{1};
ui_obj = ui_data{2};
win_str = get(ui_han(48),'userdata');
ui_han2 = get(ui_han(53),'userdata');

if ~figflag(['Controls Tutor Help',win_str]) | subject == 0, % initial page/setup

 if ~figflag(['Controls Tutor Help',win_str]),
  scrn_size = get(0,'screensize');
  fig_h = 378;
  fig_w = 498;
  fig_lft = (scrn_size(3) - fig_w)/2;
  fig_btm = (scrn_size(4) - (fig_h+20))/2;

  f2 = figure('pos',[fig_lft,fig_btm,fig_w,fig_h],'numbertitle','off','menubar','none',...
              'name',['Controls Tutor Help',win_str],'userdata',f);

  ui_han2(1) = uicontrol('style','text','pos',[5,331,488,17],...
              'horizontalalignment','center','foregroundcolor','k',...
              'backgroundcolor',[0.8,0.8,0.8]);
  ui_han2(2) = uicontrol('style','edit','pos',[5,5,488,323],...
              'horizontalalignment','left','enable','inactive','max',2,...
              'background','w');

  ui_han2(3) = uicontrol('style','push','pos',[20,353,60,20],...
              'string','Close','callback','set(gcf,''vis'',''off'')');
  ui_han2(4) = uicontrol('style','push','pos',[149,353,60,20],...
              'string','Next >>');
  ui_han2(5) = uicontrol('style','push','pos',[86,353,60,20],...
              'string','<< Prev');

  popup_str = ['Introduction|Program Conventions|Input Center|File Menu|',...
               'Options Menu|Format Menu|View Menu|Frequency Response|',...
               'Root Locus|Gain Plots|Time Response|Getting HELP'];

  ui_han2(6) = uicontrol('style','popup','pos',[214,353,160,20],'string',...
                     popup_str,'callback','userhelp(1)','backgroundcolor','w');

  set(ui_han2,'units','norm');

  set(ui_han(53),'userdata',ui_han2);

 end

 if subject == 0,  subject = 2; end

end

lab_han = ui_han2(1);
txt_han = ui_han2(2);
nxt_han = ui_han2(4);
prv_han = ui_han2(5);
pop_han = ui_han2(6);

cr = sprintf('\n');
tab = [32,32,32];

if subject == 1,
 subject = get(pop_han,'value')+1;
 page = 1;
else
 set(pop_han,'value',subject-1);
end

if subject == 2, % Introduction

 set(lab_han,'string','Introduction - page 1 of 1');
 str = [cr,'Welcome to the Controls Tutor On-line Help',cr,...
        cr,...
        'The Controls Tutor On-line Help facility provides you with a quick',cr,...
        'and easy reference to the many options provided by the Controls Tutor.',cr,...
        cr,...
        'Use the << Prev and Next >> buttons to turn the pages of this Help',cr,...
        'facility.',cr,...
        cr,...
        'Use the popup menu to select a specific area that you need help with.',cr,...
        cr,...
        'Words that are capitalized and separated with a -> represent pulldown',cr,...
        'menu options.  For example, FILE->OPEN means to select the Open option',cr,...
        'under the File pulldown menu.',cr,...
        cr,...
        'It is recommended that if you do not have a manual that you should',cr,...
        'quickly browse this Help facility.'];

 set(txt_han,'string',str);
 set(prv_han,'enable','off');
 set(nxt_han,'callback','userhelp(3,1)','enable','on');

elseif subject == 3, % Program Conventions

 if page == 1,
  set(prv_han,'enable','on','callback','userhelp(2,1)');
  set(nxt_han,'enable','on','callback','userhelp(3,2)');

  set(lab_han,'string','Program Conventions - page 1 of 5');

  str = [cr,...
         'The conventions for the Controls Tutor can be broken up into 5',cr,...
         'subject areas:',cr,...
         cr,...
         tab,'Screen Layout',cr,...
         tab,'Current Items',cr,...
         tab,'System',cr,...
         tab,'Element',cr,...
         tab,'Mouse',cr,...
         cr,...
         'Screen Layout',cr,...
         '------------------------',cr,...
         cr,...
         'The screen layout consists of the following:',cr,...
         cr,...
         tab,'Menus - Pull down menus defining all available options.',cr,...
         cr,...
         tab,'Input Center - Tools provided for manipulating transfer function',cr,...
         tab,'elements.  Located in the upper-left hand corner.'];

  set(txt_han,'string',str);

 elseif page == 2,
  set(prv_han,'enable','on','callback','userhelp(3,1)');
  set(nxt_han,'enable','on','callback','userhelp(3,3)');

  set(lab_han,'string','Program Conventions - page 2 of 5');

  str = [cr,...
         tab,'Pointer Location - Display of mouse location.  Located below',cr,...
         tab,'current plot button.',cr,...
         cr,...
         tab,'System Information - Various information regarding the current',cr,...
         tab,'system.  Located below current system info button.',cr,...
         cr,...
         tab,'Transfer Function - Display of accepted elements.  Located at',cr,...
         tab,'the top of the page.',cr,...
         cr,...
         tab,'Page - Display area that contains all plotting environments.',cr,...
         cr,...
         tab,'Status Bar - Display of program status.',cr,...
         cr,...
         'Current Items',cr,...
         '-----------------------',cr,...
         cr,...
         'The Controls Tutor works with a number of ''current'' items: the',cr,...
         'current system, the current transfer function, and the current',cr,...
         'plot'];

  set(txt_han,'string',str);

 elseif page == 3,
  set(prv_han,'enable','on','callback','userhelp(3,2)');
  set(nxt_han,'enable','on','callback','userhelp(3,4)');

  set(lab_han,'string','Program Conventions - page 3 of 5');

  str = [tab,'The current system is one of the three systems that are provided',cr,...
         tab,'by the Controls Tutor.  The System Information section tells you what',cr,...
         tab,'system you are currently modifying.  If you display all the systems,',cr,...
         tab,'the current system will be green and the non-current red.',cr,...
         cr,...
         tab,'The current transfer function is one of the three transfer functions',cr,...
         tab,'provided by the Controls Tutor: plant, forward controller, or feed-',cr,...
         tab,'back controller.  The Input Center section indicates what transfer',cr,...
         tab,'function you are currently modifying.',cr,...
         cr,...
         tab,'The current plot is one of the plots that is displayed on your page.',cr,...
         tab,'The Pointer Location section indicates what plot you are currently',cr,...
         tab,'allowed to modify.'];

  set(txt_han,'string',str);

 elseif page == 4,
  set(prv_han,'enable','on','callback','userhelp(3,3)');
  set(nxt_han,'enable','on','callback','userhelp(3,5)');

  set(lab_han,'string','Program Conventions - page 4 of 5');

  str = [cr,...
         'System',cr,...
         '------------',cr,...
         cr,...
         tab,'The Controls Tutor defines a system as containing a forward',cr,...
         tab,'controller, G(s), a plant, P(s), and a feedback controller, H(s).',cr,...
         tab,'Hence, the open-loop equation would be G(s)P(s)H(s) and the',cr,...
         tab,'closed-loop equation would be G(s)P(s)/(1 + G(s)P(s)H(s)).',cr,...
         cr,...
         tab,'The Controls Tutor allows you to work with up to 3 separate systems.',cr...
         cr,...
         'Element',cr,...
         '---------------',cr,...
         cr,...
         tab,'The Controls Tutor defines an element as a component of a transfer',cr,...
         tab,'function.  For example, integrators, real poles, and complex zeros',cr,...
         tab,'are all elements.'];

  set(txt_han,'string',str);

 elseif page == 5,
  set(prv_han,'enable','on','callback','userhelp(3,4)');
  set(nxt_han,'enable','on','callback','userhelp(4,1)');

  set(lab_han,'string','Program Conventions - page 5 of 5');

  str = [cr,...
         'Mouse',cr,...
         '-----------',cr,...
         cr,...
         tab,'The Controls Tutor makes extensive use of the mouse.  The mouse',cr,...
         tab,'can be used to identify parts of a plot and also select elements',cr,...
         tab,'from a displayed transfer function.  It can also be used to switch',cr,...
         tab,'between systems by double-clicking on a non-current system.',cr,...
         cr,...
         tab,'The results of a mouse selection are displayed both in the Status',cr,...
         tab,'Bar and Pointer Location.'];

  set(txt_han,'string',str);
 end

elseif subject == 4, % Input Center

 if page == 1,
  set(prv_han,'enable','on','callback','userhelp(3,5)');
  set(nxt_han,'enable','on','callback','userhelp(4,2)');

  set(lab_han,'string','The Input Center - page 1 of 2');

  str = [cr,...
         'The Input Center gives you the interface to enter various',cr,...
         'transfer function elements.  This is accomplished with the following',cr,...
         'interface objects:',cr,...
         cr,...
         tab,'Popup Menu - used to choose your transfer function element',cr,...
         tab,'Edit Boxes - used to enter values related to your chosen element',cr,...
         tab,'Sliders - used to iterate over values of selected elements',cr,...
         tab,'Add/Update Button - used to accept values in edit boxes',cr,...
         tab,'Delete Button - used to remove selected element',cr,...
         cr,...
         'The Input Center is placed into Add mode when an element is selected',cr,...
         'using the Popup Menu.',cr,...
         cr,...
         'The Input Center is placed into Update mode when an element is',cr,...
         'selected using the mouse pointer.  This is accomplished by depressing',cr,...
         'the mouse button while over the desired element in the displayed',cr,...
         'transfer function.  However, in order for you to select the element,',cr,...
         'the transfer function must be displayed with FORMAT->INDIVIDUAL and',cr,...
         'FORMAT->FACTORED selected.'];

  set(txt_han,'string',str);

 elseif page == 2,
  set(prv_han,'enable','on','callback','userhelp(4,1)');
  set(nxt_han,'enable','on','callback','userhelp(5,1)');

  set(lab_han,'string','The Input Center - page 2 of 2');

  str = [cr,...
         'Transfer function elements provided by the popup menu:',cr,...
         cr,...
         'Gain - update the gain of the current transfer function.',cr,...
         'Integrator - enter the number (n) of integrators for 1/s^n.',cr,...
         'Differentiator - enter the number (n) of differentiators for s^n.',cr,...
         'Real Pole - enter the pole value (p) for 1/(s + p).',cr,...
         'Real Zero - enter the zero value (z) for (s + z).',cr,...
         'Complex Pole - enter zeta and wn for 1/(s^2 + 2*zeta*wn*s + wn^2).',cr,...
         'Complex Zero - enter zeta and wn for (s^2 + 2*zeta*wn*s + wn^2).',cr,...
         'Lead/Lag - enter the phase change and corner frequency.',cr,...
         'Num/Den - enter the numerator and denominator in descending powers',cr,...
         tab,tab,'of s and enclosed in brackets. i.e. [1, 2, 3]',cr,...
         'Delay Time - enter the time value (T) for exp(-s*T).',cr,...
         'Design Point - enter the zeta value and wn value.  These values will',cr,...
         tab,tab,'then be applied to the root locus plot in the form of a zeta',cr,...
         tab,tab,'line and wn circle.'];

  set(txt_han,'string',str);
 end

elseif subject == 5, % File Menu

 if page == 1,
  set(prv_han,'enable','on','callback','userhelp(4,2)');
  set(nxt_han,'enable','on','callback','userhelp(5,2)');

  set(lab_han,'string','The File Menu - page 1 of 3');

  str = [cr,...
         'The File menu controls all the operations involving the disk,',cr,...
         'MATLAB workspace, printer, and exiting the Controls Tutor.',cr,...
         cr,...
         'NEW clears The Controls Tutor environment.',cr,...
         cr,...
         'OPEN... displays a dialog box that allows you to select a',cr,...
         'Tutor file which was saved with the SAVE option.',cr,...
         cr,...
         'SAVE saves the current transfer function settings.  If',cr,...
         'it is a [Untitled], then the SAVE AS... option is invoked.',cr,...
         cr,...
         'SAVE AS... displays a dialog box that allows you to specify a file',cr,...
         'name to store your Controls Tutor transfer function settings. Tutor',cr,...
         'files are designated with the *.tut extension.'];

  set(txt_han,'string',str);

 elseif page == 2,
  set(prv_han,'enable','on','callback','userhelp(5,1)');
  set(nxt_han,'enable','on','callback','userhelp(5,3)');

  set(lab_han,'string','The File Menu - page 2 of 3');

  str = [cr,...
         'IN FROM WORKSPACE... allows you to take variables from the',cr,...
         'MATLAB workspace and pass them into the Controls Tutor to define',cr,...
         'any of the three transfer functions.',cr,...
         cr,...
         'The formats and transfer functions are selected using the provided',cr,...
         'radio buttons and pressing OK brings the data into the transfer function.',cr,...
         cr,...
         'There are three different formats which the Controls Tutor can accept:',cr,...
         'Num/Den, Zero/Pole/Gain, and State-Space.  Each of these formats can',cr,...
         'then be passed into the selected transfer function: plant, forward',cr,...
         'controller, or feedback controller.  To specify two different transfer',cr,...
         'functions, you need to choose this option twice.',cr,...
         cr,...
         'Note: the Controls Tutor can only handle SISO systems.  If you try to',cr,...
         'pass in a MIMO system, the Tutor will only take the transfer function',cr,...
         'related to the first input.'];

  set(txt_han,'string',str);

 elseif page == 3,
  set(prv_han,'enable','on','callback','userhelp(5,2)');
  set(nxt_han,'enable','on','callback','userhelp(6,1)');

  set(lab_han,'string','The File Menu - page 3 of 3');

  str = [cr,...
         'OUT TO WORKSPACE... works the same as IN FROM WORKSPACE... except',cr,...
         'that it places its results into the MATLAB workspace with your specified',cr,...
         'variable names.',cr,...
         cr,...
         'PRINT prints all graphical information currently displayed:',cr,...
         'transfer function, plots, pointer location, and system information.',cr,...
         cr,...
         'EXIT ends the Controls Tutor session.'];

  set(txt_han,'string',str);

 end

elseif subject == 6, % Options Menu

 if page == 1,
  set(prv_han,'enable','on','callback','userhelp(5,3)');
  set(nxt_han,'enable','on','callback','userhelp(6,2)');

  set(lab_han,'string','The Options Menu - page 1 of 3');

  str = [cr,...
         'The Options menu determines the current system, current transfer',cr,...
         'function, and type of controller.',cr,...
         cr,...
         'SYSTEM allows you to switch between the 3 systems provided',cr,...
         'by the Controls Tutor.  The SHOW OTHERS option allows you view all',cr,...
         '3 systems at one time.  The current system will be green and the non-',cr,...
         'current systems will be red.',cr,...
         cr,...
         'Moving between systems can be done by using this menu option, double-',cr,...
         'clicking on a non-current system with the mouse, or pressing the',cr,...
         'System # Info button at the left.',cr,...
         cr,...
         'PLANT, P(S) allows you to set the plant of the current system',cr,...
         'as the displayed transfer function.'];

  set(txt_han,'string',str);

 elseif page == 2,
  set(prv_han,'enable','on','callback','userhelp(6,1)');
  set(nxt_han,'enable','on','callback','userhelp(6,3)');

  set(lab_han,'string','The Options Menu - page 2 of 3');

  str = [cr,...
         'FORWARD CONTROLLER, G(S) allows you to set the forward',cr,...
         'controller of the current system as the displayed transfer function.',cr,...
         cr,...
         'FEEDBACK CONTROLLER, H(S) allows you to set the feedback',cr,...
         'controller of the current system as the displayed transfer function.',cr,...
         cr,...
         'SET G(S) = 1 allows you to set the forward controller to',cr,...
         'one.  This option changes depending upon which transfer function is',cr,...
         'current.  If the forward controller were current, then this option',cr,...
         'would be SET P(S) = 1.',cr,...
         cr,...
         'SET H(S) = 1 allows you to set the feedback controller to',cr,...
         'one.  This option changes depending upon which transfer function is',cr,...
         'current.  If the feedbacl controller were current, then this option',cr,...
         'would be SET G(S) = 1.'];

  set(txt_han,'string',str);

 elseif page == 3,
  set(prv_han,'enable','on','callback','userhelp(6,2)');
  set(nxt_han,'enable','on','callback','userhelp(7,1)');

  set(lab_han,'string','The Options Menu - page 3 of 3');

  str = [cr,...
         'CONTROLLER TYPE allows you to change the format of your',cr,...
         'controller from a num/den format to one of two PID (proportional-',cr,...
         'integral-derivative) formats.',cr,...
         cr,...
         'The first PID format is of the basic equation:',cr,...
         cr,...
         'G(s) = Kp + Ki/s + Kds',cr,...
         cr,...
         'The second PID format is a manipulation of the first and is made up',cr,...
         'of a PD component and a PI component and has the following equation:',cr,...
         cr,...
         'G(s) = (1 + KDs)(KP + KI/s)',cr,...
         cr,...
         'where',cr,...
         cr,...
         'Kp = KP + KD*KI',cr,...
         'Ki = KI',cr,...
         'Kd = KD*KP'];

  set(txt_han,'string',str);

 end
elseif subject == 7, % Format Menu

 if page == 1,
  set(prv_han,'enable','on','callback','userhelp(6,3)');
  set(nxt_han,'enable','on','callback','userhelp(8,1)');

  set(lab_han,'string','The Format Menu - page 1 of 1');

  str = [cr,...
         'The Format menu controls how you want to view the current',cr,...
         'transfer function.  The first three options and the second three',cr,...
         'options are separated so that you can mix and match to determine',cr,...
         'your desired format.',cr,...
         cr,...
         'INDIVIDUAL allows you to display the current transfer',cr,...
         'function.',cr,...
         cr,...
         'OPEN-LOOP allows you to display the current system''s',cr,...
         'open-loop transfer function.',cr,...
         cr,...
         'CLOSED-LOOP allows you to display the current system''s',cr,...
         'closed-loop transfer function.',cr,...
         cr,...
         'FACTORED allows you to display the transfer function in',cr,...
         'elemental form.',cr,...
         cr,...
         'POLYNOMIAL allows you to display the transfer function in',cr,...
         'descending powers of s.',cr,...
         cr,...
         'INVERSE LAPLACE allows you to display the transfer function',cr,...
         'in its inverse Laplace form.'];

  set(txt_han,'string',str);
 end

elseif subject == 8, % View Menu

 if page == 1,
  set(prv_han,'enable','on','callback','userhelp(7,1)');
  set(nxt_han,'enable','on','callback','userhelp(8,2)');

  set(lab_han,'string','The View Menu - page 1 of 2');

  str = [cr,...
         'The View Menu allows you to control the axis limits, turn the grid',cr,...
         'on or off, and change The Control''s Tutor page layout.',cr,...
         cr,...
         'FULL allows you to set the axis limits so that all visible',cr,...
         'objects can be viewed.',cr,...
         cr,...
         'IN allows you to zoom in by 20%.',cr,...
         cr,...
         'OUT allows you to zoom out by 20%.',cr,...
         cr,...
         'ZOOM allows you to define a specific region to be your new',cr,...
         'axis limits.  When selected, the mouse pointer will change into a',cr,...
         'crosshair.  Press and drag the mouse to change the bounding box that',cr,...
         'appears; releasing the mouse button will change the axis limits.'];

  set(txt_han,'string',str);

 elseif page == 2,
  set(prv_han,'enable','on','callback','userhelp(8,1)');
  set(nxt_han,'enable','on','callback','userhelp(9,1)');

  set(lab_han,'string','The View Menu - page 2 of 2');

  str = [cr,...
         'LAST returns you to the previously defined axis limits.',cr,...
         cr,...
         'GRID allows you to toggle the grids on or off for the current',cr,...
         'axis.',cr,...
         cr,...
         'PAGE LAYOUT... allows you to change which plots you want, and',cr,...
         'on what pages you want them.',cr,...
         cr,...
         'The Controls Tutor has 8 pages that plots can be placed on.  Each page',cr,...
         'can hold up to 4 plots.  You move from page to page with the popup',cr,...
         'menu on the left.  Plot selection is performed by selecting from the',cr,...
         'list of 8 plots.  A plot is selected when its background color changes',cr,...
         'from white to its unique color.  At the same time, you will see the same',cr,...
         'color patch appear on the sample page to the right.  Page configurations',cr,...
         'are changed by selecting and dragging the patch objects.',cr,...
         cr,...
         'Clear - clears all plots from all pages.',cr,...
         'Default - returns all pages to their default layouts.'];

  set(txt_han,'string',str);

 end
elseif subject == 9, % Frequency Response

 if page == 1,
  set(prv_han,'enable','on','callback','userhelp(8,2)');
  set(nxt_han,'enable','on','callback','userhelp(9,2)');

  set(lab_han,'string','Frequency Response - page 1 of 2');

  str = [cr,...
         'The Tools Menu for Frequency Response are independent choices',cr,...
         'and toggle between checked and unchecked.  When an option is',cr,...
         'selected it becomes checked.  The Controls Tutor will always',cr,...
         'display that option until it becomes deselected, or unchecked.',cr,...
         cr,...
         'OPEN-LOOP allows you to display the open-loop frequency response.',cr,...
         cr,...
         'CLOSED-LOOP allows you to display the closed-loop frequency response.',cr,...
         cr,...
         'EXACT allows you to display the exact frequency response.',cr,...
         cr,...
         'STRAIGHT-LINE allows you to display the straight line approximation',cr,...
         'of the frequency response.  This option is only applicable to the',cr,...
         'Bode plots.'];

  set(txt_han,'string',str);

 elseif page == 2,
  set(prv_han,'enable','on','callback','userhelp(9,1)');
  set(nxt_han,'enable','on','callback','userhelp(10,1)');

  set(lab_han,'string','Frequency Response - page 2 of 2');

  str = [cr,...
         'MARGINS allows you to display the gain and phase margins of the',cr,...
         'frequency response.',cr,...
         cr,...
         'FREQUENCY... allows you to change the current working frequency',cr,...
         'vector.  The frequency vector is defined in the same manner as MATLAB''s',cr,...
         'logspace function.  The first and second edit boxes define the beginning',cr,...
         'and ending frequencies (log base10) and the third defines the number',cr,...
         'of frequencies to span this range.  For example, we could enter a -2',cr,...
         'in the first edit box, a -3 in the second, and a 100 in the third.',cr,...
         'This would define a frequency vector from 10^-2 (0.01) to 10^3 (1000)',cr,...
         'with a length of 100.'];

  set(txt_han,'string',str);

 end
elseif subject == 10, % Root Locus

 if page == 1,
  set(prv_han,'enable','on','callback','userhelp(9,2)');
  set(nxt_han,'enable','on','callback','userhelp(10,2)');

  set(lab_han,'string','Root Locus - page 1 of 3');

  str = [cr,...
         'The Tools Menu for Root Locus are independent choices and',cr,...
         'toggle between checked and unchecked.  When an option is',cr,...
         'selected it becomes checked.  The Controls Tutor will always',cr,...
         'display that option until it becomes deselected, or unchecked.',cr,...
         cr,...
         'POSITIVE LOCUS allows you to set all computations to be for the',cr,...
         'positive root locus, meaning the gain values range from 0 to +inf.',cr,...
         cr,...
         'NEGATIVE LOCUS allows you to set all computations to be for the',cr,...
         'negative root locus, meaning the gain values range from 0 to -inf.',cr,...
         cr,...
         'REAL AXIS allows you to display the real axis root loci.',cr,...
         cr,...
         'ASYMPTOTES allows you to display the asymptotes of the root locus.',cr,...
         cr,...
         'BREAKAWAY-REENTRY allows you to display where the root loci either',cr,...
         'breakaway or reenter the real axis.'];

  set(txt_han,'string',str);

 elseif page == 2,
  set(prv_han,'enable','on','callback','userhelp(10,1)');
  set(nxt_han,'enable','on','callback','userhelp(10,3)');

  set(lab_han,'string','Root Locus - page 2 of 3');

  str = [cr,...
         'IMAGINARY CROSSINGS allows you to display where the root loci cross',cr,...
         'the imaginary axis.',cr,...
         cr,...
         'ARRIVAL-DEPARTURE allows you to display the arrival and departure',cr,...
         'angles in degrees from complex zeros and poles, respectively.',cr,...
         cr,...
         'ROOT LOCUS allows you to compute the entire root locus of the current',cr,...
         'system',cr,...
         cr,...
         'CLOSED-LOOP POLES allows you to display the location of the closed-',cr,...
         'loop poles for the current system.'];

  set(txt_han,'string',str);

 elseif page == 3,
  set(prv_han,'enable','on','callback','userhelp(10,2)');
  set(nxt_han,'enable','on','callback','userhelp(11,1)');

  set(lab_han,'string','Root Locus - page 3 of 3');

  str = [cr,...
         'GAIN... allows you to manually set the gains used to compute the root',cr,...
         'locus.  The gain vector is defined in the same manner as MATLAB''s',cr,...
         'logspace function.  The first and second edit boxes define the beginning',cr,...
         'and ending gains (log base10) and the third defines the number of gains',cr,...
         'to span this range.  For example, we could enter a -2 in the first',cr,...
         'edit box, a -3 in the second, and a 100 in the third.  This would',cr,...
         'define a gain vector from 10^-2 (0.01) to 10^3 (1000) with a length',cr,...
         'of 100.',cr,...
         cr,...
         'The Auto-Select option tells the Controls Tutor to iteratively select',cr,...
         'its own gain values.  If this option is selected, all values in the',cr,...
         'edit boxes are cleared.'];

  set(txt_han,'string',str);
 end

elseif subject == 11, % Gain Plots

 if page == 1,
  set(prv_han,'enable','on','callback','userhelp(10,2)');
  set(nxt_han,'enable','on','callback','userhelp(12,1)');

  set(lab_han,'string','Gain Plots - page 1 of 1');

  str = [cr,...
         'The gain plots are related to root locus in the same way that bode',cr,...
         'plots are related to nyquist.  For nyquist and root locus, the vary-',cr,...
         'ing parameters, frequency and gain, respectively, are implicit; you',cr,...
         'cannot readily see which value of the result corresponds to the value',cr,...
         'of the parameter.  However, in bode and gain plots, this varying',cr,...
         'parameter is displayed along the x-axis and the corresponding result',cr,...
         'is split into a magnitude and phase component and displayed along the',cr,...
         'y-axis.  Hence, gain plots display the magnitude and phase component',cr,...
         'of the closed-loop roots of the characteristic equation.',cr,...
         cr,...
         'The only tools supplied for the gain plots is the ROOT LOCUS->GAIN...',cr,...
         'option.  Press the << Prev button for a description.'];

  set(txt_han,'string',str);

 end
elseif subject == 12, % Time Response

 if page == 1,
  set(prv_han,'enable','on','callback','userhelp(11,1)');
  set(nxt_han,'enable','on','callback','userhelp(12,2)');

  set(lab_han,'string','Time Response - page 1 of 3');

  str = [cr,...
         'The Tools Menu for Time Response are independent choices and',cr,...
         'toggle between checked and unchecked.  When an option is',cr,...
         'selected it becomes checked.  The Controls Tutor will always',cr,...
         'display that option until it becomes deselected, or unchecked.',cr,...
         cr,...
         'OPEN-LOOP allows you to display all the selected inputs with the',cr,...
         'open-loop transfer function.',cr,...
         cr,...
         'CLOSED-LOOP allows you to display all the selected inputs with the',cr,...
         'closed-loop transfer function.',cr,...
         cr,...
         'UNIT IMPULSE allows you to display the current system''s response to',cr,...
         'a unit impulse input.',cr,...
         cr,...
         'UNIT STEP allows you to display the current system''s response to a',cr,...
         'unit step input.',cr,...
         cr,...
         'UNIT RAMP allows you to display the current system''s response to a',cr,...
         'unit ramp input.'];

  set(txt_han,'string',str);

 elseif page == 2,
  set(prv_han,'enable','on','callback','userhelp(12,1)');
  set(nxt_han,'enable','on','callback','userhelp(12,3)');

  set(lab_han,'string','Time Response - page 2 of 3');
  str = [cr,...
         'RISE TIME allows you to display the rise time for the open- or',cr,...
         'closed-loop step response.  The rise time is the time required for',cr,...
         'the underdamped response to rise from 10 to 90% of its final value;',cr,...
         'or the overdamped response to rise from 0 to 100% of its final value',cr,...
         cr,...
         'DELAY TIME allows you to display the delay time for the open- or',cr,...
         'closed-loop step response.  The delay time is the time required for',cr,...
         'the response to reach half the final value the very first time.',cr,...
         cr,...
         'PEAK TIME allows you to display the peak time for the open- or',cr,...
         'closed-loop step response.  The peak time is the time required for',cr,...
         'the response to reach the first peak.  If the system is over-damped,',cr,...
         'peak time can not be computed.',cr,...
         cr,...
         '% OVERSHOOT allows you to display the percent overshoot for the open-',cr,...
         'or closed-loop step response.  The percent overshoot is the precentage',cr,...
         'amount that the response proceeds beyond the desired response.'];

  set(txt_han,'string',str);

 elseif page == 3,
  set(prv_han,'enable','on','callback','userhelp(12,2)');
  set(nxt_han,'enable','on','callback','userhelp(13,1)');

  set(lab_han,'string','Time Response - page 3 of 3');
  str = [cr,...
         '2% SETTLING allows you to display the 2% settling time for the open-',cr,...
         'or closed-loop step response.  The 2% settling time is the time',cr,...
         'required for the response curve to reach and stay within a 2% range',cr,...
         'about the final value.',cr,...
         cr,...
         '5% SETTLING allows you to display the 5% settling time for the open-',cr,...
         'or closed-loop step response.  The 5% settling time is the time',cr,...
         'required for the response curve to reach and stay within a 5% range',cr,...
         'about the final value.',cr,...
         cr,...
         'TIME... allows you to set the values at which to compute the time',cr,...
         'response.  The time vector is defined in the same manner as MATLAB''s',cr,...
         'linspace function.  The first and second edit boxes define the be-',cr,...
         'ginning and ending time values and the third defines the number of',cr,...
         'time values to span this range.'];

  set(txt_han,'string',str);

 end

elseif subject == 13, % Getting HELP

 if page == 1,
  set(prv_han,'enable','on','callback','userhelp(12,3)');
  set(nxt_han,'enable','off');

  set(lab_han,'string','Getting HELP - page 1 of 1');

  str = [cr,...
         'If you encounter problems either related to software or a general',cr,...
         'misunderstanding, send email to the appropriate address:',cr,...
         cr,...
         tab,'info@terasoft.com',cr...
         tab,'support@terasoft.com',cr,...
         tab,'suggest@terasoft.com',cr,...
         tab,'bugs@terasoft.com',cr,...
         cr,...
         'If this is an illegal copy, please consider purchasing the program',cr...
         'through proper channels so that you can become a registered user.'];

  set(txt_han,'string',str);

 end
end


text:
function [ cvx_optval, success ] = quad_form( x, Q, tol )

%QUAD_FORM   Internal cvx version.

%
% Check sizes and types
%

error( nargchk( 2, 3, nargin ) );
if nargin < 3, tol = 4 * eps; end
sx = size( x );
if length( sx ) ~= 2 || all( sx > 1 ),
    error( 'The first argument must be a row or column.' );
else
    sx = prod( sx );
end

sQ = size( Q );
if length( sQ ) ~= 2 || sQ( 1 ) ~= sQ( 2 ),
    error( 'The second argument must be a scalar or a square matrix.' );
elseif sQ( 1 ) ~= sx && sQ( 1 ) ~= 1,
    error( 'Sizes are incompatible.' );
end

x = vec( x );
success = true;
if cvx_isconstant( x ),

    if isreal( Q ) || isreal( x ),

        %
        % Constant x, affine Q, real case
        %

        x = real( x );
        Q = real( Q );
        cvx_optval = x' * ( Q * x );

    else

        %
        % Constant x, affine Q, complex case
        %

        xR = real( x );
        xI = imag( x );
        cvx_optval = xR' * ( real( Q ) * xR ) + xI' * ( imag( Q ) * xI );

    end

elseif ~cvx_isaffine( x ),

    error( 'First argument must be affine.' );

elseif sQ( 1 ) == 1,

    %
    % Constant scalar Q, affine x
    %

    cvx_optval = real( Q ) * sum_square_abs( x );

else

    %
    % Constant matrix Q, affine x
    %

    cvx_optval = [];
    while true,
        Q = cvx_constant( Q );

        %
        % Quick exit for a zero Q
        %

        nnzQ = nnz( Q );
        if nnzQ == 0,
            cvx_optval = 0;
            break
        end

        %
        % Remove zero rows and columns from Q. If a diagonal element of Q is
        % zero but there are elements on that row or column that are not,
        % then we know that neither Q nor -Q is PSD.
        %

        dQ = diag( Q );
        if ~all( dQ ),
            tt = dQ ~= 0;
            Q = Q( tt, tt );
            if nnz( Q ) ~= nnzQ,
                break
            end
            dQ = dQ( tt );
            x = cvx_subsref( x, tt, ':' );
        end

        %
        % Determine the sign of the elements of Q. If they are not all of
        % the same sign, then neither Q nor -Q is PSD.
        %

        dQ = dQ > 0;
        if all( dQ ),
            alpha = +1;
        elseif any( dQ ),
            break
        else
            alpha = -1;
            Q = -Q;
        end

        %
        % Now perform a Cholesky factorization. If it succeeds then we
        % know that alpha * Q is PSD.
        %

        Q = 0.5 * ( Q + Q' );
        if cvx_use_sparse( Q ),
            Q = sparse( Q );
            prm = symamd( Q );
            R = cholinc( Q( prm, prm ), 'inf' );
            R( :, prm ) = R;
            tt = any( isinf( R ), 2 );
            valid = ~any( tt );
            if ~valid,
                R( tt, : ) = [];
            end
        else
            [ R, p ] = chol( full( Q ) );
            valid = p == 0;
            if ~valid,
                R = [ R, R' \ Q(1:p-1,p:end) ];
            end
        end
        if ~valid,
            valid = false; % normest( Q - R' * R ) < tol * normest( Q );
        end

        %
        % If more accuracy is needed, perform a Schur decomposition.
        %

        if valid,
            Dmax = max(R(:));
            R = R / Dmax;
            alpha = alpha * Dmax * Dmax;
        else
            [ V, D ] = eig( full( Q ) );
            if cvx_use_sparse( V ),
                V = sparse( V );
            end
            D = diag( D );
            Dmax = max( D );
            Derr = tol * Dmax;
            if min( D ) < - Derr,
                break
            end
            tt = find( D > Derr );
            alpha = alpha * Dmax;
            R = diag(sparse(sqrt(D(tt)/Dmax))) * V( :, tt )';
        end

        cvx_optval = alpha * sum_square_abs( R * x );
        success = true;
        break;

    end

    if isempty( cvx_optval ),
        if nargout > 1,
            success = false;
        else
            error( 'The second argument must be positive or negative semidefinite.' );
        end
    end

end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function [ fvP ] = primalObjD( Xdiag, yvect, Th_vecIdx, lambda1, lambda2, W )
% primal objective (diagonalized)
%  P(W)  sum_i^m ||Xi wi - yi|| + lambda1 ||W||_{1,2} + lambda2/2 ||W||_F^2
fvP = lambda1 * sum(sqrt(sum(W.^2, 2))) + lambda2 /2 * sum(sum(W.^2))...
    + segL2 (Xdiag * W(:) - yvect, Th_vecIdx);

end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function [clusLengths,biggestEvent,mbg,bg,clustnumbers] = funBuildclu(mycat, biggestEvent, clus, mbg)
% FUNBUILDCLU builds cluster out out of information stored in clus
%
%  FUNBUILDCLU(mycat, biggestEvent, clus, mbg)
%
% mycat : catalog
% clus: size of catalog, containing cluster index numbers
% biggestEvent: size nClusters
% mbg: max mag in cluster
% bg : size nClusters
%
    %  originally, "mycat" was "newcat"
    % A.Allmann

% modified by C Reyes 2017

%global mycat clus mbg k1 clust clustnumbers
%global bgevent bg

 % count # events in each cluster
clusLengths = histcounts(clus,'BinLimits',[1 max(clus)], 'BinMethod','integers');

% cull any clusters that don't have events
empty_clusters = clusLengths == 0;      %numbers of clusters that are empty  eg. [1:9 2:0 3: 8 4:0] ->  [0 1 0 1]
clusLengths(empty_clusters)=[];
biggestEvent(empty_clusters)=[];
mbg(empty_clusters)=[];

bg = biggestEvent;
biggestEvent = mycat.subset(bg); %biggest event in a cluster(if more than one,take first)

clustnumbers=(1:numel(clusLengths));    % stores numbers of clusters

author: CelsoReyes


text:
function label_vec  = convert_matlabel_to_vec(label_mat, num)
    label_vec = zeros(1, num);
    for i=1:num
        nonzero_idx = find(label_mat(:,i)>0);
        if length(nonzero_idx) > 1
            nonzero_idx = min(nonzero_idx);
        end
        label_vec(1, i) = nonzero_idx;
    end
end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
%% Copyright (C) 2014, 2016, 2018-2019 Colin B. Macdonald
%% Copyright (C) 2016 Lagu
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym real (@var{z})
%% Real part of a symbolic expression.
%%
%% Examples:
%% @example
%% @group
%% syms z
%% real(z)
%%   @result{} ans = (sym) re(z)
%% @end group
%%
%% @group
%% syms x real
%% real(x)
%%   @result{} ans = (sym) x
%%
%% real([x  sym(pi) + 6i  7  3i])
%%   @result{} ans = (sym) [x  π  7  0]  (1×4 matrix)
%% @end group
%% @end example
%%
%% @seealso{@@sym/imag, @@sym/conj, @@sym/ctranspose}
%% @end defmethod

function y = real(z)
  if (nargin ~= 1)
    print_usage ();
  end

  y = elementwise_op ('re', z);

end

%!assert (isequal (real (sym (4) + 3i),4))

%!test
%! syms x y real
%! z = x + 1i*y;
%! assert (isequal (real (z),x))

%!test
%! syms x y real
%! Z = [4  x + 1i*y; x 4 + 3i];
%! assert (isequal (real (Z),[4 x; x 4]))

%!test
%! syms x real
%! d = exp (x*i);
%! assert (isequal (real (d), cos (x)))

%!test
%! % round trip
%! syms x
%! d = 3 - 5i;
%! f = real (x);
%! A = real (d);
%! h = function_handle (f);
%! B = h (d);
%! assert (A, B)

repo: octsympy


text:
%  PARSE_ARGUMENTS -- arguments parser
%
%  More About
%  ------------
%
%  - This function is an alias for utils.miscellaneous.parse_arguments.
%  Check it out for more explanations
%
%  Examples
%  ---------
%
%  See also: UTILS.MISCELLANEOUS.PARSE_ARGUMENTS
%

author: jmaih