text:
function [xmin, fmin] = golden(npd, f, ax, bx, cx, tol, varargin)
%
%GOLDEN   Minimize function of one variable using golden section search
%
%   [xmin, fmin] = golden(npd, f, ax, bx, cx, tol) computes a local minimum
%   of f. xmin is the computed local minimizer of f and fmin is
%   f(xmin). xmin is computed to an relative accuracy of TOL.
%
%   The parameters ax, bx and cx must satisfy the following conditions:
%   ax < bx < cx, f(bx) < f(ax) and f(bx) < f(cx).
%
%   xmin satisfies ax < xmin < cx. golden is guaranteed to succeed if f
%   is continuous between ax and cx
%
%   Roman Geus, ETH Zuerich, 9.12.97
%
%   ATI -- added "npd" argument & made private to KDE class
C = (3-sqrt(5))/2;
R = 1-C;

x0 = ax;
x3 = cx;
if (abs(cx-bx) > abs(bx-ax)),
  x1 = bx;
  x2 = bx + C*(cx-bx);
else
  x2 = bx;
  x1 = bx - C*(bx-ax);
end
f1 = feval(f,x1,npd,varargin{:});
f2 = feval(f,x2,npd,varargin{:});

k = 1;
while abs(x3-x0) > tol*(abs(x1)+abs(x2)),
%  fprintf(1,'k=%4d, |a-b|=%e\n', k, abs(x3-x0));
  if f2 < f1,
    x0 = x1;
    x1 = x2;
    x2 = R*x1 + C*x3;   % x2 = x1+c*(x3-x1)
    f1 = f2;
    f2 = feval(f,x2,npd,varargin{:});
  else
    x3 = x2;
    x2 = x1;
    x1 = R*x2 + C*x0;   % x1 = x2+c*(x0-x2)
    f2 = f1;
    f1 = feval(f,x1,npd,varargin{:});
  end
  k = k+1;

%  [x0,x1,x2,x3,f1,f2]
end

if f1 < f2,
  xmin = x1;
  fmin = f1;
else
  xmin = x2;
  fmin = f2;
end

author: ShapeNet

repo: RenderForCNN


text:
function dfdocontext(varargin)
%DFDOCONTEXT Perform context menu actions for distribution fitting tool

% Copyright 2001-2004 The MathWorks, Inc.
% $Revision: 1.1.6.7 $  $Date: 2004/03/09 16:17:04 $
import com.mathworks.toolbox.stats.*;

% Special action to create context menus
if isequal(varargin{1},'create')
   makecontextmenu(varargin{2});
   return
end

% Get information about what invoked this function
obj = gcbo;
action = get(obj,'Tag');
h = gco;
if isempty(h), return; end
dffig = gcbf;

% Set up variables that define some menu items
[sizes styles markers] = getmenuitems;
styles{end+1} = 'none';

changed = true;   % did a line property change?

switch action

 % This case is triggered when we display the menu
 case {'fitcontext' 'datacontext'}
   % Store a handle to the object that triggered this menu
   set(obj,'UserData',h);

   hObject = get(h,'UserData');
   c = findall(obj,'Type','uimenu');
   hBounds = findall(c,'flat','Tag','confbounds');
   ftype = dfgetset('ftype');

   % Enable or disable as appropriate
   hMarker = findall(c,'flat','Tag','marker');
   if isequal(action,'datacontext')
      hLineStyle = findall(c,'flat','Tag','linestyle');
      hLineWidth = findall(c,'flat','Tag','linewidth');
      hBinRules = findall(c,'flat','Tag','binrules');
      if isequal(ftype,'probplot')
         set(hMarker,'Enable','on');
         set(hLineStyle,'Enable','off');
         set(hLineWidth,'Enable','off');
         set(hBinRules,'Enable','off');
      elseif isequal(ftype,'pdf')
         set(hMarker,'Enable','off');
         set(hLineStyle,'Enable','on');
         set(hLineWidth,'Enable','on');
         set(hBinRules,'Enable','on');
      else
         set(hMarker,'Enable','off');
         set(hLineStyle,'Enable','on');
         set(hLineWidth,'Enable','on');
         set(hBinRules,'Enable','off');
      end
   else
      if ~hObject.iscontinuous && isequal(ftype,'pdf')
         set(hMarker,'Enable','on');
      else
         set(hMarker,'Enable','off');
      end
   end
   try
      hasconfbounds = hObject.distspec.hasconfbounds;
   catch
      hasconfbounds = false;
   end
   if isequal(ftype,'pdf') || isequal(ftype,'probplot') || ...
      (isequal(action,'fitcontext') && ~hasconfbounds) || ...
      (isequal(action,'datacontext') && isequal(ftype,'icdf'))
      set(hBounds,'Enable','off');
   else
      set(hBounds,'Enable','on');
   end

   set(c,'Checked','off');

   % Fix check mark for confidence bounds
   if hObject.showbounds
      set(hBounds,'Checked','on');
   end

   % Fix check marks on line width and line style cascading menus
   w = get(h,'LineWidth');
   u = findall(c,'flat','Tag',num2str(w));
   if ~isempty(u)
      set(u,'Checked','on');
   end
   w = get(h,'LineStyle');
   u = findall(c,'flat','Tag',w);
   if ~isempty(u)
      set(u,'Checked','on');
   end
   w = get(h,'Marker');
   u = findall(c,'flat','Tag',w);
   if ~isempty(u)
      set(u,'Checked','on');
   end
   return

 % Remaining cases are triggered by selecting menu items
 case 'confbounds'
   hObject = get(h,'UserData');
   hObject.showbounds = ~hObject.showbounds;
   nm = get(hObject,'name');
   htag = get(h,'Tag');
   if isequal(htag,'dfdata')
      DataSetsManager.getDataSetsManager.dataSetChanged(java(hObject),nm,nm);
   else
      FitsManager.getFitsManager.fitChanged(java(hObject),nm,nm);
   end

 case 'color'
   oldcolor = get(h,'Color');
   newcolor = uisetcolor(oldcolor);
   if ~isequal(oldcolor,newcolor)
      set(h,'Color',newcolor);
   end

 case styles
   set(h,'LineStyle',action);

 case markers
   if isequal(action,'point')
      msize = 12;
   else
      msize = 6;
   end
   set(h,'Marker',action,'MarkerSize',msize);

 % Either delete a fit, or a hide a fit or data set
 case {'hidecurve' 'deletefit'}
   htag = get(h,'Tag');
   if isequal(htag,'distfit') || isequal(htag,'dfdata')
      hndl = get(h,'UserData');
      if isequal(action,'hidecurve')
         hndl.plot = 0;
         nm = get(hndl,'name');
      else
         % The delete action appears on the fit menu only, not the data set menu
         FitsManager.getFitsManager.deleteFits(java(hndl));
      end
   end
   changed = false;

 % Edit a fit
 case 'editfit'
   htag = get(h,'Tag');
   if isequal(htag,'distfit')  % should always be true
      hndl = get(h,'UserData');
      FitsManager.getFitsManager.editFit(hndl.name);
   end
   changed = false;

 % Bring up the "Set Bin Width Rules" dialog for this data set
 case 'binrules'
   htag = get(h,'Tag');
   if isequal(htag,'dfdata')
      hndl = get(h,'UserData');
      nm = get(hndl,'name');
      bw = com.mathworks.toolbox.stats.BinWidth.getBinWidth; % get dialog
      bw.displayBinWidth(nm); % display dialog for this data set
   end
   changed = false;

 % If the menu item is a number, it is a line width
 otherwise
   j = str2num(action);
   if ~isempty(j)
      set(h,'LineWidth',j);
   end

end

if changed
   % Save plot info in the fit or data set object
   hObject = get(h,'UserData');
   savelineproperties(hObject);

   % Update legend
   dfupdatelegend(dffig);
end

% ---------------------- helper to make context menu
function makecontextmenu(dffig)
%MAKECONTEXTMENU Creates context menu for curve fitting figure

% Create context menus for fits, data curve, probability plot data curves
cFit = uicontextmenu('Parent',dffig,'Tag','fitcontext','Callback',@dfdocontext);
uimenu(cFit,'Label','Color...','Tag','color','Callback',@dfdocontext);

% Add menu items for line and marker control
uwidth = uimenu(cFit,'Label','Line &Width','Tag','linewidth');
ustyle = uimenu(cFit,'Label','Line &Style','Tag','linestyle');
umark = uimenu(cFit,'Label','Marker','Tag','marker','Position',2);

% Add menu items to control confidence bounds
uimenu(cFit,'Label','Confidence &Bounds','Callback',@dfdocontext,...
            'Tag','confbounds');

% Get menu item labels and tags
[sizes styles markers slabels mlabels] = getmenuitems;

for j=1:length(markers)
   uimenu(umark,'Label',mlabels{j},'Callback',@dfdocontext,'Tag',markers{j});
end

% Sub-menus for line widths
for i = 1:length(sizes)
   val = num2str(sizes(i));
   uimenu(uwidth,'Label',val,'Callback',@dfdocontext,'Tag',val);
end

% Sub-menus for line styles
for j=1:length(styles)
   uimenu(ustyle,'Label',slabels{j},'Callback',@dfdocontext,'Tag',styles{j});
end

% Copy the fit menu to create a data menu
cData = copyobj(cFit,dffig);
set(cData,'Tag','datacontext')

% Add items for fit menus only
uimenu(cFit,'Label','&Hide Fit','Tag','hidecurve','Callback',@dfdocontext,...
       'Separator','on');
uimenu(cFit,'Label','&Delete Fit','Tag','deletefit','Callback',@dfdocontext);
uimenu(cFit,'Label','&Edit Fit','Tag','editfit','Callback',@dfdocontext);

% Add items for data menus only
uimenu(cData,'Label','&Hide Data','Tag','hidecurve',...
       'Callback',@dfdocontext,'Separator','on');
uimenu(cData,'Label','Set Bin &Rules','Tag','binrules',...
       'Callback',@dfdocontext,'Separator','on');

% -------------- helper to get menu item labels
function [sizes,styles,markers,slabels,mlabels] = getmenuitems
%GETMENUITEMS Get items for curve fitting context menus
sizes = [0.5 1 2 3 4 5 6 7 8 9 10];
styles = {'-' '--' ':' '-.'};
markers = {'+' 'o' '*' '.' 'x' 'square' 'diamond' ...
        'v' '^' '<' '>' 'pentagram' 'hexagram'};
slabels = {'solid' 'dash' 'dot' 'dash-dot'};
mlabels = {'plus' 'circle' 'star' 'point' 'x-mark' 'square' 'diamond' ...
           'triangle (down)' 'triangle (up)' 'triangle (left)' ...
           'triangle (right)' 'pentagram' 'hexagram'};

author: zouchuhang

repo: LayoutNet


text:
classdef MimWSEditSlice < MimModel
    methods (Access = protected)
        function value = run(obj)
            value = struct();
        end

        function ValueHasChanged(obj, value)
            overlayImage = obj.Callback.getModelValue(obj.Parameters.segmentationVolumeId);
            imageSliceNumber = obj.Parameters.imageSliceNumber;
            axialDimension = obj.Parameters.axialDimension;
            disp('New value for seg slice:');
            disp(value);

            slice = overlayImage.GetSlice(imageSliceNumber, axialDimension);
            strokesArray = value.data;

            for stroke = strokesArray
                colour = MimWSEditSlice.stringToColour(stroke{1}.colour);
                pointsArray = stroke{1}.points;
                for point = pointsArray
                    x = point{1}.x;
                    y = point{1}.y;
                    slice(max(1, round(x)), max(1, round(y))) = colour;
                end
            end

            overlayImage.ReplaceImageSlice(slice, imageSliceNumber, axialDimension);
            obj.Callback.setModelValue(obj.Parameters.segmentationVolumeId, overlayImage.Copy());
        end
    end

    methods (Static, Access = private)
        function colourIndex = stringToColour(colourString)
            switch(colourString)
                case 'blue'
                    colourIndex = 1;
                case 'green'
                    colourIndex = 2;
                case 'red'
                    colourIndex = 3;
                case 'cyan'
                    colourIndex = 4;
                case 'magenta'
                    colourIndex = 5;
                case 'yellow'
                    colourIndex = 6;
                case 'grey'
                    colourIndex = 7;
                otherwise
                    colourIndex = 1;
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [Fnew, obj, M,k,x,u,n,deg,linears,nonlinears,vecConstraints,isinequality,ulong] = momentmodel(F,obj,k,keepnonlinears)

if nargin < 2
    obj = [];
end
if nargin < 3
    k = [];
end
if nargin < 4
    keepnonlinears = 0;
end

vecConstraints = [];
sdpConstraints = [];
isinequality = [];
binaries = [];
xvars = [];
Fnew = ([]);
for i = 1:length(F)
    if is(F(i),'elementwise')
        X = sdpvar(F(i));
        vecConstraints = [vecConstraints;X(:)];
        isinequality = [isinequality ones(1,prod(size(X)))];
        xvars = [xvars depends(X(:))];
    elseif is(F(i),'equality')
        X = sdpvar(F(i));
        if is(X,'symmetric')
            X = X(find(triu(ones(length(X)))));
        end
        vecConstraints = [vecConstraints;-X(:)];
        isinequality = [isinequality zeros(1,prod(size(X)))];
        xvars = [xvars depends(X(:))];
    elseif is(F(i),'sdp')
        sdpConstraints{end+1} = sdpvar(F(i));
        xvars = [xvars depends(F(i))];
    elseif is(F(i),'binary')
        binaries = [binaries getvariables(F(i))];
    else
        Fnew = Fnew+F(i); % Should only be SOCP constraints
    end
end

% Recover the involved variables
x = recover(unique([depends(obj) xvars]));
n = length(x);

% Check degrees of constraints
deg = [];
for i = 1:length(vecConstraints)
    deg(end+1) = degree(vecConstraints(i));
end
for i = 1:length(sdpConstraints)
    deg(end+1) = degree(sdpConstraints{i});
end
if isempty(deg)
    deg = 0;
end

% Create lowest possible relaxation if k=[]
d = ceil((max(degree(obj),max(deg)))/2);
k_min = d;
if isempty(k)
    k = k_min;
else
    if k<k_min
        error('Higher order relaxation needed')
    end
end

% Generate monomials of order k
u{k} = monolist(x,k);
ulong{k} = monolist(x,2*k);

% Largest moment matrix. NOTE SHIFT M{k+1} = M_k.
M{k+1}=u{k}*u{k}';
% Moment matrices easily generated with this trick
% The matrices will NOT be rank-1 since the products
% generate the relaxed variables

% ... and lower degree localization matrices
M{1} = 1;
for i = 1:1:k-1;
    n_i = round(factorial(n+k-i)/(factorial(n)*factorial(k-i)));
    M{k-i+1} = M{k+1}(1:n_i,1:n_i);
end

% Lasserres relaxation (Lasserre, SIAM J. OPTIM, 11(3) 796-817)
Fmoments = (M{k+1}>=0);
for i = 1:length(vecConstraints)
    if isinequality(i)
        v_k = floor((degree(vecConstraints(i))+1)/2);
        Localizer = vecConstraints(i)*M{k-v_k+1};
        if isa(vecConstraints(i),'double')
            if vecConstraints(i)<0
                error('Problem is trivially infeasible due to negative constant')
            else
                continue
            end
        end
        Fmoments = Fmoments+(Localizer>=0);
    else
        if isa(vecConstraints(i),'double')
            if vecConstraints(i)~=0
                error('Problem is trivially infeasible due to non-zero constant in equality constraints')
            else
                continue
            end
        end
        Localizer = vecConstraints(i)*monolist(x,2*k-degree(vecConstraints(i)));
        Fmoments = Fmoments+(Localizer==0);
    end
end
for i = 1:length(sdpConstraints)
    v_k = floor((degree(sdpConstraints{i})+1)/2);
    Fmoments = Fmoments+(kron(M{k-v_k+1},sdpConstraints{i})>=0);
end

% Add them all
Fnew = Fnew + Fmoments;

% Get all binary and reduce problem
binaries = union(binaries,yalmip('binvariables'));
if ~isempty(binaries)
    if isa(obj,'sdpvar')
        obj = eliminateBinary(obj,binaries);
    end
    for i = 1:length(Fmoments)
        Fnew(i) = eliminateBinary(Fnew(i),binaries);
    end
    for i = 2:1:k+1;
        M{i} = eliminateBinary(M{i},binaries);
    end
end

vars = getvariables(Fnew);
for i = 1:length(M)
    vars = [vars getvariables(M{i})];
end
vars = unique([vars getvariables(obj)]);
[mt,variabletype] = yalmip('monomtable');
nonlinears = vars(find(variabletype(vars)));
newLinear = sdpvar(length(nonlinears),1);

if isa(obj,'sdpvar')
    obj = variablereplace(obj,nonlinears,getvariables(newLinear));
end
for i = 1:length(M)
    if isa(M{i},'sdpvar')
        M{i} = variablereplace(M{i},nonlinears,getvariables(newLinear));
    end
end
linears = getvariables(newLinear);
Fnew = variablereplace(Fnew,nonlinears,getvariables(newLinear));
linears = recover(linears);
nonlinears = recover(nonlinears);

end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [J, grad] = lrCostFunction(theta, X, y, lambda)
%LRCOSTFUNCTION Compute cost and gradient for logistic regression with
%regularization
%   J = LRCOSTFUNCTION(theta, X, y, lambda) computes the cost of using
%   theta as the parameter for regularized logistic regression and the
%   gradient of the cost w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Hint: The computation of the cost function and gradients can be
%       efficiently vectorized. For example, consider the computation
%
%           sigmoid(X * theta)
%
%       Each row of the resulting matrix will contain the value of the
%       prediction for that example. You can make use of this to vectorize
%       the cost function and gradient computations.
%
% Hint: When computing the gradient of the regularized cost function,
%       there're many possible vectorized solutions, but one solution
%       looks like:
%           grad = (unregularized gradient for logistic regression)
%           temp = theta;
%           temp(1) = 0;   % because we don't add anything for j = 0
%           grad = grad + YOUR_CODE_HERE (using the temp variable)
%

sigmo = sigmoid(X * theta);

Jreg = lambda / (2 * m) * (theta' * theta - theta(1) .^2);

J = 1 / m * (-y' * log(sigmo) - (1 - y)' * log(1 - sigmo)) + Jreg;

temp = theta;
temp(1) = 0;
gradreg = lambda / m * temp;

grad = 1 / m * X' * (sigmo - y) + gradreg;

% =============================================================

grad = grad(:);

end

author: benoitvallon

repo: coursera-machine-learning


text:
%TESTR MSE for regression
%
%      E = TESTR(X,W,TYPE)
%      E = TESTR(X*W,TYPE)
%      E = X*W*TESTR([],TYPE)
%      E = X*W*TESTR(TYPE)
%
% INPUT
%   X    Regression dataset
%   W    Regression mapping
%   TYPE Type of error measure, default: mean squared error
%
% OUTPUT
%   E    Mean squared error
%
% DESCRIPTION
% Compute the error of regression W on dataset X. The following error
% measures have been defined for TYPE:
% 'mse'    mean squared error (default)
% 'mad'    mean absolute deviation
%
% SEE ALSO (PRTools Guide)
%  RSQUARED, TESTC

% Copyright: D.M.J. Tax, D.M.J.Tax@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function e = testr(varargin)

  argin = shiftargin(varargin,'char');
  argin = setdefaults(argin,[],[],'mse');

  if mapping_task(argin,'definition')
    e = define_mapping(argin,'fixed');

  else	% Evaluate.

    [x,w,type] = deal(argin{:});

    if (ismapping(w) & istrained(w))
      x = x*w;
    end
    if ischar(w)
      type = w;
    end
    switch type
      case 'mse'
        e = mean((+x(:,1) - gettargets(x)).^2);
      case 'mad'
        e = mean(abs(+x(:,1) - gettargets(x)));
      otherwise
        error('Error %s is not implemented.',type);
    end

    if nargout==0
      %display results on the screen:
      fprintf('Error on %d objects: %f.\n',...
        size(x,1), e);
      clear e;
    end

  end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [fmin,fconstr]=fun(A,p,q)
% Example of target function to be used by SIMPS minimizer

% Handling of optional problem-dependent parameters
if nargin>1
%
 if nargin>2
%
 end
end

n=prod(size(A));
B=zeros(1,n); B(:)=A;
B=(sin((B.*(1:n)).^2)).^20;
fmin=abs(det(A))/500+sum(B)/1+norm(fix(n*A))/200;
% All three components are scaled to about the same order of
% magnitude. Observe that constants 500, 1 and 200
% - as well as 10 and 5 used for the constraints (see below) -
% are obtained statistically and they are hard-coded
% exclusively for the case of three-dimensional matrices A.

fconstr=[trace(abs(A))-10,5-min(abs(eig(A)))];

return

author: Sable

repo: mcbench-benchmarks


text:
function fv_fscore= proc_fisherScore(fv, varargin)
%PROC_FISHERSCORE - Computes the Fisher Score for each feature
%
%Synopsis:
% FV_SCORE= proc_fisherScore(FV, <OPT>)
%
%Arguments:
% FV  - data structure of feature vectors
% OPT - struct or property/value list of optional properties
%    .TolerateNans: observations with NaN value are skipped
%       (nanmean/nanstd are used instead of mean/std)
%
%Returns:
% FV_SCORE - data structute of Fisher scores (one sample only)
%
%See also:
% proc_t_scaled, proc_r_values, proc_r_square, proc_rocAreaValues
%
%Comment:
% Only standard case is tested.

% Author(s): Benjamin Blankertz
props= { 'TolerateNans'   0    '!BOOL'};

if nargin==0,
  fv_fscore = props; return
end

fv = misc_history(fv);
misc_checkType(fv, 'STRUCT(x clab y)');
opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

if size(fv.y,1)>2,  %% multi-class: TODO! now doing it pairwise:
  warning('calculating pairwise Fisher scores');
  combs= nchoosek(1:size(fv.y,1), 2);
  for ic= 1:length(combs),
    ep= proc_selectClasses(fv, combs(ic,:));
    fv0= proc_fisherScore(ep);
    if ic==1,
      fv_fscore= fv0;
    else
      fv_fscore= proc_appendEpochs(fv_fscore, fv0);
    end
  end
  return;
elseif size(fv.y,1)==1,
  warning('1 class only: calculating Fisher score against flat-line of same var');
  fv2= fv;
  szx= size(fv.x);
  fv2.x= fv2.x - repmat(mean(fv2.x,3), [1 1 size(fv2.x,3)]);
  fv2.className= {'flat'};
  fv= proc_appendEpochs(fv, fv2);
end

sz= size(fv.x);
fv.x= reshape(fv.x, [prod(sz(1:end-1)), sz(end)]);
cl1= find(fv.y(1,:));
cl2= find(fv.y(2,:));
if opt.TolerateNans,
  me= [nanmean(fv.x(:,cl1),2) nanmean(fv.x(:,cl2),2)];
  va= [nanstd(fv.x(:,cl1),0,2).^2 nanstd(fv.x(:,cl2),0,2).^2];
else
  me= [mean(fv.x(:,cl1),2) mean(fv.x(:,cl2),2)];
  va= [var(fv.x(:,cl1)')' var(fv.x(:,cl2)')'];
end
fscore= ((me(:,1)-me(:,2)).^2)./(va(:,1)+va(:,2)+eps);
fscore= reshape(fscore, [sz(1:end-1) 1]);

fv_fscore= rmfield(fv, intersect(fieldnames(fv),{'x','y','className'},'legacy') );
fv_fscore.x= fscore;
if isfield(fv, 'className'),
  fv_fscore.className= {sprintf('Fs( %s , %s )', fv.className{1:2})};
end
fv_fscore.y= 1;
fv_fscore.yUnit= 'au';

author: bbci


text:
classdef prtClassFld < prtClass
 %prtClassFld Fisher linear discriminant classifier
 %
 %    CLASSIFIER = prtClassFld returns a Fisher linear discriminant classifier
 %
 %    CLASSIFIER = prtClassFld(PROPERTY1, VALUE1, ...) constructs a
 %    prtClassFld object CLASSIFIER with properties as specified by
 %    PROPERTY/VALUE pairs.
 %
 %    A prtClassFld object inherits all properties from the abstract class
 %    prtClass. In addition is has the following properties:
 %
 %    w                  - regression weights, estimated during training
 %    plotBasis          - Flag indicating whether to plot the basis
 %                         functions when the PLOT function is called
 %    plotProjections    - Flag indicating whether to plot the projection
 %                         of points to the basis when the PLOT function is
 %                         called
 %
 %    For information on the Fisher Linear Discriminant algorithm, please
 %    refer to the following URL:
 %
 %    http://en.wikipedia.org/wiki/Linear_discriminant_analysis#Fisher.27s_linear_discriminant
 %
 %    A prtClassFld object inherits the TRAIN, RUN, CROSSVALIDATE and
 %    KFOLDS methods from prtAction. It also inherits the PLOT method from
 %    prtClass.
 %
 %    Example:
 %
 %    ds1 = prtDataGenUnimodal;       % Create some test and
 %    ds2 = prtDataGenUnimodal;   % training data
 %    classifier = prtClassFld;           % Create a classifier
 %    classifier = classifier.train(ds1);    % Train
 %    classified = run(classifier, ds2);         % Test
 %    subplot(2,1,1);
 %    classifier.plot;
 %    subplot(2,1,2);
 %    [pf,pd] = prtScoreRoc(classified);
 %    h = plot(pf,pd,'linewidth',3);
 %    title('ROC'); xlabel('Pf'); ylabel('Pd');
 %
 %   See also prtClass, prtClassLogisticDiscriminant, prtClassBagging,
 %   prtClassMap, prtClassCap, prtClassBinaryToMaryOneVsAll, prtClassDlrt,
 %   prtClassPlsda, prtClassKnn, prtClassRvm, prtClassGlrt,  prtClassSvm,
 %   prtClassTreeBaggingCap, prtClassKmsd, prtClassKnn

% Copyright (c) 2013 New Folder Consulting
%
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the
% "Software"), to deal in the Software without restriction, including
% without limitation the rights to use, copy, modify, merge, publish,
% distribute, sublicense, and/or sell copies of the Software, and to permit
% persons to whom the Software is furnished to do so, subject to the
% following conditions:
%
% The above copyright notice and this permission notice shall be included
% in all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
% OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
% MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
% NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
% DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
% OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
% USE OR OTHER DEALINGS IN THE SOFTWARE.

    properties (SetAccess=private)

        name = 'Fisher Linear Discriminant' % Fisher Linear Discriminant
        nameAbbreviation = 'FLD'            % FLD
        isNativeMary = false;  % False
    end

    properties (SetAccess = protected)
        % w is a dataSet.nDimensions x 1 vector of projection weights
        % learned during Fld.train(dataSet)

        w = []; % The vector of weights, learned during training

        % plotting options
        plotBasis = false; % Flag indicating whether or not to plot the basis
        plotProjections = false; % Flag indicating whether or not to plot the projections
    end

    methods

               % Allow for string, value pairs
        function self = prtClassFld(varargin)
            self = prtUtilAssignStringValuePairs(self,varargin{:});
        end

        function self = set.plotProjections(self,value)
            if islogical(value) || (isnumeric(value) && (value == 1 || value == 0))
                self.plotProjections = value;
            else
                error('prt:prtClassFld:plotProjections','plotProjections can only take true or false (boolean or 0/1) values; user speficied value %d',value);
            end
        end
    end

    methods (Access=protected, Hidden = true)

        function self = trainAction(self,dataSet)

            n = dataSet.nObservations;
            p = dataSet.nFeatures;

            if p > n
                warning('prt:prtClassFld:train:illconditioned','dataSet has n (%d) < p (%d); prtClassFld may not be stable',n,p);
            end
            if ~dataSet.isBinary
                error('prtClassFld:nonBinaryTraining','Input dataSet for prtClassFld.train must be binary');
            end

            dataH0 = dataSet.getObservationsByClassInd(1);
            dataH1 = dataSet.getObservationsByClassInd(2);

            mean0 = mean(dataH0,1);
            mean1 = mean(dataH1,1);

            cov0 = cov(dataH0);
            cov1 = cov(dataH1);
            covW = cov1 + cov0;

            self.w = covW\(mean1-mean0)'; %w = covW^-1 * (mean1-mean0)'; But better
            self.w = self.w./norm(self.w);

        end

        function dataSet = runAction(self,dataSet)
            %dataSet = prtDataSetClass((self.w'*dataSet.getObservations()')');
            dataSet.X = (self.w'*dataSet.getObservations()')';
        end

        function imageHandle = plotGriddedEvaledClassifier(self, DS, linGrid, gridSize, cMap)

            % Call the original plot function
            imageHandle = plotGriddedEvaledClassifier@prtClass(self, DS, linGrid, gridSize, cMap);

            W = self.w;
            limits = axis;
            nDims = length(W);

            if self.plotBasis
                hold on
                switch nDims
                    case 1
                        % Nothing
                    case 2
                        distances = zeros(4,1);
                        distances(1) = sqrt(sum([limits(2); limits(4)].^2));
                        distances(2) = sqrt(sum([limits(1); limits(3)].^2));
                        distances(3) = sqrt(sum([limits(2); limits(3)].^2));
                        distances(4) = sqrt(sum([limits(1); limits(4)].^2));

                        highPoint =  max(distances).*W;
                        lowPoint =  -max(distances).*W;

                        h = plot([lowPoint(1),highPoint(1)],[lowPoint(2),highPoint(2)],'k');
                        set(h,'linewidth',3);
                    case 3
                        distances = zeros(8,1);
                        distances(1) = sqrt(sum([limits(1); limits(3); limits(5)].^2));
                        distances(2) = sqrt(sum([limits(1); limits(3); limits(6)].^2));
                        distances(3) = sqrt(sum([limits(1); limits(4); limits(5)].^2));
                        distances(4) = sqrt(sum([limits(1); limits(4); limits(6)].^2));
                        distances(5) = sqrt(sum([limits(2); limits(3); limits(5)].^2));
                        distances(6) = sqrt(sum([limits(2); limits(3); limits(6)].^2));
                        distances(7) = sqrt(sum([limits(2); limits(4); limits(5)].^2));
                        distances(8) = sqrt(sum([limits(2); limits(4); limits(6)].^2));

                        highPoint =  max(distances).*W;
                        lowPoint =  -max(distances).*W;

                        h = plot3([lowPoint(1),highPoint(1)],[lowPoint(2),highPoint(2)],[lowPoint(3), highPoint(3)],'k');
                        set(h,'linewidth',3);
                    otherwise
                        error('prt:prtClassFld:tooManyDimensions','Too many dimensions for plotting.')
                end
            end

            if self.plotProjections && ~isempty(self.dataSet)
                OutputDataSet = run(self, self.dataSet);
                hold on;
                switch nDims
                    case 2
                        for i = 1:double(self.plotProjections):self.dataSet.nObservations
                            cX = self.dataSet.X(i,:);
                            cYout = OutputDataSet.X(i,:);
                            plot([cX(1),cYout*W(1)],[cX(2),cYout*W(2)],'k');
                        end
                    case 3
                        for i = 1:double(self.plotProjections):self.dataSet.nObservations
                            cX = self.dataSet.X(i,:);
                            cYout = OutputDataSet.X(i,:);
                            plot3([cX(1),cYout*W(1)],[cX(2),cYout*W(2)],[cX(3),cYout*W(3)],'k');
                        end
                end
                axis(limits);
            end
            hold off;
        end

    end

end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function spm_eeg_inv_group(S)
% Source reconstruction for a group ERP or ERF study
% FORMAT spm_eeg_inv_group(S)
%
% S  - string array of names of M/EEG mat files for inversion (optional)
%__________________________________________________________________________
%
% spm_eeg_inv_group inverts forward models for a group of subjects or ERPs
% under the simple assumption that the [empirical prior] variance on each
% source can be factorised into source-specific and subject-specific terms.
% These covariance components are estimated using ReML (a form of Gaussian
% process modelling) to give empirical priors on sources.  Source-specific
% covariance parameters are estimated first using the sample covariance
% matrix in sensor space over subjects and trials using multiple sparse
% priors (and,  by default, a greedy search).  The subject-specific terms
% are then estimated by pooling over trials for each subject separately.
% All trials in D.events.types will be inverted in the order specified.
% The result is a contrast (saved in D.mat) and a 3-D volume of MAP or
% conditional estimates of source activity that are constrained to the
% same subset of voxels.  These would normally be passed to a second-level
% SPM for classical inference about between-trial effects, over subjects.
%__________________________________________________________________________
%
% References:
% Electromagnetic source reconstruction for group studies. V. Litvak and
% K.J. Friston. NeuroImage, 42:1490-1498, 2008.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_eeg_inv_group.m 6862 2016-08-25 14:42:19Z guillaume $

SVNrev = '$Rev: 6862 $';

%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename, SVNrev);

%-Check if to proceed
%--------------------------------------------------------------------------
str = questdlg({'This will overwrite previous source reconstructions.', ...
    'Do you wish to continue?'},'M/EEG Group Inversion','Yes','No','Yes');
if ~strcmp(str,'Yes'), return, end

% Load data
%==========================================================================

% Give file names
%--------------------------------------------------------------------------
if ~nargin
    [S, sts] = spm_select(Inf, 'mat','Select M/EEG mat files');
    if ~sts, return; end
end
Ns    = size(S,1);
swd   = pwd;

% Load data and set method
%==========================================================================
for i = 1:Ns

    fprintf('checking for previous inversions: subject %i\n',i);
    D{i}                 = spm_eeg_load(deblank(S(i,:)));
    D{i}.val             = 1;
    D{i}.inv{1}.method   = 'Imaging';

    % clear redundant models
    %----------------------------------------------------------------------
    D{i}.inv = D{i}.inv(1);

    % clear previous inversions
    %----------------------------------------------------------------------
    try, D{i}.inv{1} = rmfield(D{i}.inv{1},'inverse' ); end
    try, D{i}.inv{1} = rmfield(D{i}.inv{1},'contrast'); end

    % save forward model parameters
    %----------------------------------------------------------------------
    save(D{i});

end

% Check for existing forward models and consistent Gain matrices
%--------------------------------------------------------------------------
Nd = zeros(1,Ns);
for i = 1:Ns
    fprintf('checking for forward models: subject %i\n',i);
    try
        [L, D{i}] = spm_eeg_lgainmat(D{i});
        Nd(i) = size(L,2);               % number of dipoles
    catch
        Nd(i) = 0;
    end
end

% use template head model where necessary
%==========================================================================
if max(Nd > 1024)
    NS = find(Nd ~= max(Nd));            % subjects requiring forward model
else
    NS = 1:Ns;
end
for i = NS

    cd(D{i}.path);

    % specify cortical mesh size (1 to 4; 1 = 5125, 2 = 8196 dipoles)
    %----------------------------------------------------------------------
    Msize  = 2;

    % use a template head model and associated meshes
    %======================================================================
    D{i} = spm_eeg_inv_mesh_ui(D{i}, 1, 1, Msize);

    % save forward model parameters
    %----------------------------------------------------------------------
    save(D{i});

end

% Get inversion parameters
%==========================================================================
inverse = spm_eeg_inv_custom_ui(D{1});

% Select modality
%==========================================================================
% Modality
%------------------------------------------------------------------
[mod, list] = modality(D{1}, 1, 1);
if strcmp(mod, 'Multimodal')
    [selection, ok]= listdlg('ListString', list, 'SelectionMode', 'multiple' ,...
        'Name', 'Select modalities' , 'InitialValue', 1:numel(list),  'ListSize', [400 300]);
    if ~ok
        return;
    end

    inverse.modality  = list(selection);

    if numel(inverse.modality) == 1
        inverse.modality = inverse.modality{1};
    end
else
    inverse.modality = mod;
end

for i = 2:Ns
    [mod, list] = modality(D{i}, 1, 1);
    if ~all(ismember(inverse.modality, list))
        error([inverse.modality ' modality is missing from ' D{i}.fname]);
    end
end

% and save them (assume trials = types)
%--------------------------------------------------------------------------
for i = 1:Ns
    D{i}.inv{1}.inverse = inverse;
end

% specify time-frequency window contrast
%==========================================================================
tfwin = spm_input('Time-Frequency contrast?','+1','y/n',[1,0],1);
if tfwin

    % get time window
    %----------------------------------------------------------------------
    woi              = spm_input('Time window (ms)','+1','r',[100 200]);
    woi              = sort(woi);
    contrast.woi     = round([woi(1) woi(end)]);

    % get frequency window
    %----------------------------------------------------------------------
    fboi             = spm_input('Frequency [band] of interest (Hz)','+1','r',0);
    fboi             = sort(fboi);
    contrast.fboi    = round([fboi(1) fboi(end)]);
    contrast.display = 0;
    contrast.smooth  = 4;

    str  = {'evoked','induced'};
    contrast.type = spm_input('Power of the energy or mean energy','+1','b',str,[],1);
else
    contrast = [];
end

% Register and compute a forward model
%==========================================================================
for i = NS

    fprintf('Registering and computing forward model (subject: %i)\n',i);

    % Forward model
    %----------------------------------------------------------------------
    D{i} = spm_eeg_inv_datareg_ui(D{i}, 1);
    try
        D{i}.inv{1}.forward.voltype = vltype;
        D{i}   = spm_eeg_inv_forward(D{i});
    catch
        D{i}   = spm_eeg_inv_forward_ui(D{i});
        vltype = D{i}.inv{1}.forward.voltype;
    end

    % save forward model
    %----------------------------------------------------------------------
    save(D{i});

end

% Invert the forward model
%==========================================================================
D     = spm_eeg_invert(D);
if ~iscell(D), D = {D}; end

% Save
%==========================================================================
for i = 1:Ns
    save(D{i});
end
clear D

% Compute conditional expectation of contrast and produce image
%==========================================================================
if ~isempty(contrast)

    % evaluate contrast and write image
    %----------------------------------------------------------------------
    for i = 1:Ns
        D     = spm_eeg_load(deblank(S(i,:)));
        D.inv{1}.contrast = contrast;
        D     = spm_eeg_inv_results(D);
        D     = spm_eeg_inv_Mesh2Voxels(D);
        save(D);
    end
end

% Cleanup
%==========================================================================
cd(swd);

author: spm


text:
function [TorF, vstr, rdate] = have_feature_yalmip()
%HAVE_FEATURE_YALMIP  Detect availability/version info for YALMIP
%
%   Feature detection function implementing 'yalmip' tag for HAVE_FEATURE
%   to detect availability/version of YALMIP modeling platform
%   (https://yalmip.github.io).
%
%   See also HAVE_FEATURE, YALMIP.

%   MP-Opt-Model
%   Copyright (c) 2004-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

TorF = exist('yalmip','file') == 2;
vstr = '';
rdate = '';
if TorF
    vstr = yalmip('version');
    if length(vstr) == 8
        yr = str2num(vstr(1:4));
        mo = str2num(vstr(5:6));
        dy = str2num(vstr(7:8));
        rdate = datestr([yr mo dy 0 0 0], 'dd-mmm-yyyy');
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% DATA
%
% Files
%   cubeafemdata     - Data of an example of AFEM in a cube
%   eddycurrentdata1 - Copyright (C) Long Chen. See COPYRIGHT.txt for details.
%   elasticitydata   - data for elasticity problem
%   fourorderdate    - data for four order problem
%   jumpdata1        - data for interface problem
%   jumpmgdata1      - Data of JUMPMG1
%   jumpmgdata2      - Data of JUMPMG1
%   Maxwelldata1     - positive definite and homogenous Neumann boundary condition
%   Maxwelldata2     - non-homogenous Dirichlet/Neumann boundary condition
%   Maxwelldata3     - polynomial data
%   Maxwelldata4     - homogenous Dirichlet boundary condition
%   Maxwelldata5     - linear polynomial data
%   mixBCdata        - mix boundary condition data for Poisson equation
%   mixBCdata3       - mix boundary condition data for Poisson equation in 3-D
%   planewavedata    - plane wave solution of Maxwell equations: complex solution
%   planewavedata1   - plane wave solution of Maxwell equations: real solution.
%   planewavedataC   - plane wave solution of Maxwell equations: real solution and complex coefficients
%   planewavedataH   - plane wave solution of Maxwell equations in terms of H
%   polydata1        - polynomial data for Poisson equation
%   polydata3        - polynomial data for Poisson equation in 3-D
%   sincosdata       - trigonometric  data for Poisson equation
%   sincosdata3      - trigonometric data for Poisson equation in 3-D
%   tensordata1      - data for interface problem

repo: ifem


text:
% ReadSegyFast : Reads a SEG Y rev 1 formatted file, without header values (faster than ReadSegy)
%
% Call :
% [Data]=ReadSegyFast(filename);
% and equivalent to :
% [Data]=ReadSegy(filename);
%
%
% Read only the data of a SegFile - NOT Their headers.
% Much faster than ReadSegy
%
% 'minmax', 'skip'
%

% Implemented using the syntax of the SEG-Y revised format :
% SEGY-Y rev 0, SEG-Y rev 1 as described in
% http://seg.org/publications/tech-stand/
%
% Extended Textual Header is not yet tested
% If you would like it implemented, please send me an SEGY file with
% that sort of information, as well as a description of the segy file
%
%
% (C) 2001, Thomas Mejer Hansen, thomas@cultpenguin.com
%
%    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
%
%
%
%
function [Data,SegyHeader,SegyTraceHeadersBinary]=ReadSegyFast(filename,varargin);

  mfilename='ReadSegyFast';

  if nargin==0,
    SegymatVerbose([mfilename,' : ', filename,' does not exist !'])
    Data=[];SegyTraceHeaders=[];SegyHeader=[];HeaderInfo=[];
    return
  end

  % Initialize varargin values
  dsf=[];
  revision=[];
  endian_tight=[];
  tmin=[];tmax=[];
  headermin=[];headermax=[];header=[];
  jump=[];
  SegyHeader=[];

  segyid = fopen(filename,'r','b');   % ALL DISK FILES ARE IN BIG
                                    % ENDIAN FORMAT, ACCORDING TO SEG
                                    % Y rev 1

  fseek(segyid,0,'eof');
  BytesInFile = ftell(segyid);
  fseek(segyid,0,'bof');

ninput=nargin;
% NEXT TWO LINES TO ENUSRE THAT VARARGIN CAN BE PASSED TO FUNCTION
if ninput==2
    % CALL USING VARARGIN
    ninput=1+length(varargin{1});
    varargin=varargin{1};
else
    % DIRECT CALL
    ninput=length(varargin);
end

% TRANSFORM VARARGING INTO PARAMETERS
cargin=1;
while (cargin<ninput)

    if strcmp(varargin{cargin},'jump')
       cargin=cargin+1;
       eval(['jump=',num2str(varargin{cargin}),';']);
       SegymatVerbose(['JUMP : Read only every ',num2str(jump),'th trace'])
    end

    if strcmp(varargin{cargin},'minmax')
       cargin=cargin+1;
       eval(['header=''',varargin{cargin},''';']);
       cargin=cargin+1;
       eval(['headermin=',num2str(varargin{cargin}),';']);
       cargin=cargin+1;
       eval(['headermax=',num2str(varargin{cargin}),';']);
       SegymatVerbose(['MIN MAX : Using header ',header,' from ',num2str(headermin),' to ',num2str(headermax)])
    end

    if strcmp(varargin{cargin},'trange')
       cargin=cargin+1;
       eval(['tmin=',num2str(varargin{cargin}),';']);
       cargin=cargin+1;
       eval(['tmax=',num2str(varargin{cargin}),';']);
       SegymatVerbose(['TRANGE : tmin=',num2str(tmin),' tmax=',num2str(tmax)])
    end

    if strcmp(varargin{cargin},'revision')
       cargin=cargin+1;
       eval(['revision=',num2str(varargin{cargin}),';']);
       SegymatVerbose(['USING REVISION : rev=',num2str(revision)])
    end

    if strcmp(varargin{cargin},'dsf')
       cargin=cargin+1;
       eval(['dsf=',num2str(varargin{cargin}),';']);
       SegymatVerbose(['USING Data Sample Format : dsf=',num2str(dsf)])
    end

    if strcmp(varargin{cargin},'SegyHeader')
       cargin=cargin+1;
       SegyHeader=varargin{cargin};
       SegymatVerbose(['USING LOADED SEGYHEADER'])
    end

    cargin=cargin+1;

end

% SKIP DATA SEEMS OUTDATED, COMMENTING OUT 2003-10-24

%if exist('SkipData')==0,
%    SkipData=0; % [0] READ ONLY HEADER VALUES, [1] READ IN ALL DATA
%end%%
%
%if SkipData==1, SegymatVerbose(['Not reading data - headers only']), end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BINARY HEADERS
if isempty(SegyHeader)==1,
  SegyHeader=GetSegyHeader(segyid);
else
  SegymatVerbose([mfilename,' - Using supplied SegyHeader'])
end

% APPLY CHANGES TO SEGY HEADER IF NEEDE
if isempty(revision)==0,
  SegyHeader.SegyFormatRevisionNumber=revision;
  SegymatVerbose([mfilename,' - Manually set SEG-Y revision to ',num2str(revision)])
end
if isempty(dsf)==0,
  SegyHeader.DataSampleFormat=dsf;
end

% JUST SOME INFORMATION TO WRITE TO SCREEN :

Revision=SegyHeader.SegyFormatRevisionNumber;
if Revision>0, Revision=1; end
FormatName=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).name;
Format=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).format;
BPS=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).bps;
txt=['SegyRevision ',sprintf('%0.4g',Revision),', ',FormatName,'(',num2str(SegyHeader.DataSampleFormat),')'];

ns=SegyHeader.ns;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% READ TEXTURAL FILE HEADER EXTENSION IF NEEDED
%if SegyHeader.NumberOfExtTextualHeaders~=0
%  SegymatVerbose(['---------------------------------------------------'])
%  SegymatVerbose(['extendeD textual file headers are not supported yet'])
%  SegymatVerbose(['They are simply skipped'])
%  SegymatVerbose(['---------------------------------------------------'])
%end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% READ DATA
%Segy=fread(segyid,4000,'float32');
fseek(segyid,0,'eof'); DataEnd=ftell(segyid);

DataStart=3600+3200*SegyHeader.NumberOfExtTextualHeaders;
fseek(segyid,DataStart,'bof');       % Go to the beginning of the file

ntraces=(DataEnd-DataStart)./(240+(SegyHeader.ns)*(BPS/8));

SegymatVerbose(['Number of Samples Per Trace=',num2str(SegyHeader.ns)])
SegymatVerbose(['Number of Traces=',num2str(ntraces)])

% Set Format :
Revision=SegyHeader.SegyFormatRevisionNumber;
if Revision>0, Revision=1; end
Format=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).format;
bitpersample=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).bps;
bytepersample=bitpersample/8;
ntraceheader=240/bytepersample;

Data=fread(segyid,[ns+ntraceheader ntraces],Format);
SegyTraceHeadersBinary=Data(1:ntraceheader,:);
Data=Data([ntraceheader+1:1:ns+ntraceheader],:);
%SegymatVerbose(['Now=',num2str(ftell(segyid)),' END=',num2str(DataEnd)])

% THE FOLLOWING DOES NOT WORK FOR OCTAVE (uint32)
SegyHeader.DataFormat=Format;
if (strcmp(Format,'uint32')==1), % IBM FLOATING POINT
        % CONVERT FROM FLOATING POINT
        verbose=2;
        SegymatVerbose([mfilename,'Converting from IBM, DataFormat :',SegyHeader.DataFormat],verbose);
        Data=ibm2num(uint32(Data));
end;

existJump=1-isempty(jump);
existHeader=1-isempty(header);
existTmin=1-isempty(tmin);
existTmax=1-isempty(tmax);

if existJump
  usetraces=[jump:jump:ntraces];
  ntraces=length(usetraces);
  Data=Data(:,usetraces);
end

if (existTmin==1)&(existTmax==1)
  SegymatVerbose('TRANGE')
  % NEXT LINE SHOULD CONSIDER THAT ns in Trace and Segy Header could vary !!!
  origtrange=[1:1:SegyHeader.ns].*SegyHeader.dt.*1e-6;
  gooddata=find(origtrange>tmin & origtrange<tmax);
  Data=Data(gooddata,:);
end

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
function [ p ] = ml_explained_variance( L, Var )
%EXPLAINED_VARIANCE Function that returns the optimal p given a desired
%   explained variance. The student should convert the Eigenvalue matrix
%   to a vector and visualize the values as a 2D plot.
%   input -----------------------------------------------------------------
%
%       o L      : (N x N), Diagonal Matrix composed of lambda_i
%
%   output ----------------------------------------------------------------
%
%       o p      : optimal principal components wrt. explained variance

% ====================== Implement Eq. 8 Here ======================
eigs = diag(L);
expl_var = eigs./sum(eigs);

% ====================== Implement Eq. 9 Here ======================
cum_expl_var = cumsum(expl_var);

% ====================== Implement Eq. 10 Here ======================
p = sum(cum_expl_var <= Var) + 1;

% Visualize Explained Variance from Eigenvalues
figure;
plot(cum_expl_var, '--r', 'LineWidth', 2) ; hold on;
plot(p,cum_expl_var(p),'or')
title('Explained Variance from EigenValues')
set(gca,'XTick',[1:1:length(eigs)])
ylabel('% Cumulative Variance Explained')
xlabel('Eigenvector index')
grid on

end

author: epfl-lasa


text:
function [f,g] = tt_cp_fun(x,Z,Znormsqr)
%TT_CP_FUN Calculate function and gradient for CP fit function.
%
%  [F,G] = TT_CP_FUN(X,Z) where X is a vector containing the entries of the
%  components of the model and Z is the tensor to be fit.
%
%  See also TT_CP_VEC_TO_FAC, TT_FAC_TO_VEC, TT_CP_FG, CP_OPT
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Convert x to a cell array of matrices
A = tt_cp_vec_to_fac(x,Z);

%% Call cp_fit and cp_gradient using cp_fg
[f,G] = tt_cp_fg(Z,A,Znormsqr);

%% Convert a cell array to a vector
g = tt_fac_to_vec(G);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% LDA - MATLAB subroutine to perform linear discriminant analysis
% by Will Dwinnell and Deniz Sevis
%
% Use:
% W = LDA(Input,Target,Priors)
%
% W       = discovered linear coefficients (first column is the constants)
% Input   = predictor data (variables in columns, observations in rows)
% Target  = target variable (class labels)
% Priors  = vector of prior probabilities (optional)
%
% Note: discriminant coefficients are stored in W in the order of unique(Target)
%
% Example:
%
% % Generate example data: 2 groups, of 10 and 15, respectively
% X = [randn(10,2); randn(15,2) + 1.5];  Y = [zeros(10,1); ones(15,1)];
%
% % Calculate linear discriminant coefficients
% W = LDA(X,Y);
%
% % Calulcate linear scores for training data
% L = [ones(25,1) X] * W';
%
% % Calculate class probabilities
% P = exp(L) ./ repmat(sum(exp(L),2),[1 2]);
%
%
% Last modified: Dec-11-2010

function W = LDA(Input,Target,Priors)

% Determine size of input data
[n m] = size(Input);

% Discover and count unique class labels
ClassLabel = unique(Target);
k = length(ClassLabel);

% Initialize
nGroup     = NaN(k,1);     % Group counts
GroupMean  = NaN(k,m);     % Group sample means
PooledCov  = zeros(m,m);   % Pooled covariance
W          = NaN(k,m+1);   % model coefficients

if  (nargin >= 3)  PriorProb = Priors;  end

% Loop over classes to perform intermediate calculations
for i = 1:k,
    % Establish location and size of each class
    Group      = (Target == ClassLabel(i));
    nGroup(i)  = sum(double(Group));

    % Calculate group mean vectors
    GroupMean(i,:) = mean(Input(Group,:));

    % Accumulate pooled covariance information
    PooledCov = PooledCov + ((nGroup(i) - 1) / (n - k) ).* cov(Input(Group,:));
end

% Assign prior probabilities
if  (nargin >= 3)
    % Use the user-supplied priors
    PriorProb = Priors;
else
    % Use the sample probabilities
    PriorProb = nGroup / n;
end

% Loop over classes to calculate linear discriminant coefficients
for i = 1:k,
    % Intermediate calculation for efficiency
    % This replaces:  GroupMean(g,:) * inv(PooledCov)
    Temp = GroupMean(i,:) / PooledCov;

    % Constant
    W(i,1) = -0.5 * Temp * GroupMean(i,:)' + log(PriorProb(i));

    % Linear
    W(i,2:end) = Temp;
end

% Housekeeping
clear Temp

end

% EOF


text:
function stage = BOPointInPolygon( P, q )
%% BOPointInPolygon - find point in polygon
%
%   REFERENCE:
%       This code is described in "Computational Geometry in C", Chapter 7.
%       The Matlab version is based on C code written by Joseph O'Rourke,
%       contributions by Min Xu, June 1997.
%
%   INPUT:
%       P       - input vector Nx2
%       q       - input query vector 1x2
%
%   OUTPUT:
%       stage   - output stage:
%                   'v' - q is a vertex
%                   'i' - q is inside
%                   'o' - q is outside
%                   'e' - q on the edge
%   USAGE:
%       t = 0:0.6:6.28;
%       seedx = 100; seedy = 100; scalefactor = 20;
%       P(:,1) = (seedx(1) + scalefactor*cos(t))';
%       P(:,2) = (seedy(1) + scalefactor*sin(t))';
%       PP = [P; P(1,:)];
%       plot(PP(:,1),PP(:,2),'-bs'); hold on
%       q = PP(10,:);
%       plot(q(:,1),q(:,2),'r*');
%       stage = BOPointInPolygon(PP,q);
%
%   AUTHOR:
%       Boguslaw Obara, http://boguslawobara.net/
%
%   VERSION:
%       0.1 - 25/05/2008 First implementation
%

%%
% i, i1     - point index; i1 = i-1 mod n
% d         - dimension index
% x         - x intersection of e with ray
% Rcross    - number of right edge/ray crossings
% Lcross    - number of left edge/ray crossings
%%
Rcross = 0; Lcross = 0;
n = length(P);
%%
% Shift so that q is the origin. Note this destroys the polygon.
for i=1:n
    for d=1:2
      P(i,d) = P(i,d) - q(d);
    end
end
%%
% For each edge e=(i-1,i), see if crosses ray.
for i=1:n
    % First see if q=(0,0) is a vertex.
    if P(i,1)==0 && P(i,2)==0
        stage = 'v';
        return;
    end
    i1 = mod((i + n - 2), n) + 1;

    % if e "straddles" the x-axis...
    % The commented-out statement is logically equivalent to the one following.
    % if( ((P(i,2)>0) && (P(i1,2)<=0)) || ((P(i1,2)>0) && (P(i,2)<=0))

    if (P(i,2)>0) ~= (P(i1,2)>0)
      % e straddles ray, so compute intersection with ray.
        x = (P(i,1) * P(i1,2) - P(i1,1) * P(i,2)) / (P(i1,2) - P(i,2));

      % crosses ray if strictly positive intersection.
      if x>0
          Rcross = Rcross + 1;
      end
    end

    % if e straddles the x-axis when reversed...
    % if ((P(i,2)<0) && (P(i1,2)>=0)) || ((P(i1,2)<0) && (P(i,2)>=0)

    if (P(i,2) < 0) ~= (P(i1,2) < 0)
      % e straddles ray, so compute intersection with ray.
      x = (P(i,1) * P(i1,2) - P(i1,1) * P(i,2)) / (P(i1,2) - P(i,2));
      % crosses ray if strictly positive intersection.
      if (x < 0)
          Lcross = Lcross + 1;
      end
    end
end
%%
% q on the edge if left and right cross are not the same parity.
if mod(Rcross,2) ~= mod(Lcross,2)
    stage = 'e';
    return;
end
% q inside if an odd number of crossings.
if mod(Rcross,2) == 1
    stage = 'i';
    return;
else
    stage = 'o';
    return;
end
%%
end


text:
% //This software is licensed under the BSD 3 Clause license: http://opensource.org/licenses/BSD-3-Clause
%
%
% //Copyright (c) 2013, University of Oxford
% //All rights reserved.
%
% //Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
%
% //Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
% //Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
% //Neither the name of the University of Oxford nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
% //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
%   The method implemented in this file has been patented by their original
%   authors. Commercial use of this code is thus strongly not
%   recocomended.
%
% //Authors: 	Gari D Clifford -
% //            Roberta Colloca -
% //			Julien Oster	-

function [ MAD ] = comput_MAD(data)

%comput_MAD computes the MAD value for the vector  'data'

%Input
%data RR interval series (column vector)

%Output
%MAD median of median value for heart rate deviation

len=length(data);
J=3; %number of segments
N=len/J;  % number of RR intervals per segment

%heart rates segments: J-1=seg(;,1) J=seg(;,2) J+1=seg(;,3)
data=reshape(data,N,J);
seg=1./data ;

%remove mean from each segment: J-1, J, J+1
seg_m=detrend(seg,'constant');

%remove linear trend
seg_m_lt=detrend(seg_m,'linear');

%determine absolute deviation from mean, absdev_i for i=1..N
%within segments J-1, J, J+1
abs_dev=abs(detrend(seg_m_lt,'constant'));

%determine the Median of absolute deviation from mean for each RR interval, mABSdev_j within segments J-1, J, J+1
 mABSDEV= median(abs_dev);

 %determine median absolute deviation of segment J, compared to adjacent
 %segments , J-1 and J+1
  MAD= median(mABSDEV);
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
  function st = ct_geom_par(type, varargin)
%|function st = ct_geom_par(type, varargin)
%|
%| Same as ct_geom, but adds a rebinned variable
%| ct_geom should be changed to include a new type instead
%|
%| Create the "CT geometry" structure that describes the sampling
%| characteristics of a cone-beam CT system (axial or helical).
%| (Use sino_geom() for 2D fan-beam or parallel-beam systems.)
%|
%| in
%|	type	'fan'		multi-slice fan-beam - recommended
%|		'par'		(parallel-beam only weakly supported)
%|
%| options for all geometries
%|	'orbit_start'		default: 0
%|	'orbit'			[degrees] default: 180 for parallel / mojette
%|					or 360 for fan (negative for CW)
%|					can be 'short' for fan-beam short scan
%|	'down'			down-sampling factor, for testing
%|					can be [down_s down_t down_a]
%|	'units'			string to print distance units (default: 'mm')
%|
%| options for fan-beam
%|	'ns'			# of horizontal samples
%|	'nt'			# of vertical samples
%|	'na' | 'nbeta'		# of angular samples
%|	'ds'			horizontal sample spacing (default: 1)
%|	'dt'			vertical sample spacing (default: -ds)
%|				or {'dz', dz} to use dz * dsd / dso (usual CT)
%|	'offset_s'		unitless fraction of a channel (default: 0)
%|				(relative to line between two center channels).
%|				use 0.25 or 1.25 for "quarter-detector offset"
%|	'offset_t'		unitless (default: 0)
%|
%| options for partial scans (added by jang-hwan cho)
%|	'nframe'		# of frames to divide orbit into (default: 1)
%|	'frame'			which frame? (default: 1)
%|
%| options for helical (or step-and-shoot)
%|	'pitch'			bed_travel_per_rotation / axial_fov. default: 0.
%|				(unitless. can be negative. usually near 1.0)
%|	'source_z0'		z-location of source for first view. Default 0.
%|					It must have same units as dt and ds.
%|				Use 'center' to center the helix around z=0.
%|
%|	'user_source_zs' [na]	user-specified source z-locations for each view.
%|				usually this is empty (default) in which case
%|				source_zs is computed internally from "pitch"
%|				It must have same units as dt and ds.
%|
%|	fan beam distances:
%|	'dsd' | 'dis_src_det'	default: inf (parallel beam)
%|	'dso' | 'dis_src_iso'	default: inf (parallel beam)
%|	'dod' | 'dis_iso_det'	default: 0
%|	'dfs' | 'dis_foc_src'	default: 0 (3rd generation CT arc),
%|					use 'inf' for flat detector
%|
%| out
%|	st	(struct)	initialized structure
%|
%| methods
%|	st.shape(sino)		reshape sinograms that are columns into 3d array
%|	st.s			s sample locations
%|	st.t			t sample locations
%|	st.gamma		[nb] gamma sample values [radians]
%|	st.gamma_max		half of fan angle [radians]
%|	st.ws			(ns-1)/2 + st.offset_s
%|	st.wt			(nt-1)/2 + st.offset_t
%|	st.ad			[na] source angles in degrees
%|	st.ar			[na] source angles in radians
%|	st.dim			dimensions: [st.ns st.nt st.na]
%|	st.downsample(down)	reduce sampling by integer factor
%|	st.ones			ones(ns,nt,na, 'single')
%|	st.zeros		zeros(ns,nt,na, 'single')
%|	st.rmax			max radius within FOV
%|	st.footprint_size(ig)	max footprint width in 's'
%|	st.zfov			axial FOV
%|	st.source_zs		[na] z-locations of source for each view
%|	st.shape(sino(:))	reshape to [ns,nt,na,?]
%|	st.unitv(is,it,ia)	unit 'vector' with one nonzero element
%|	st.plot([ig])		show geometry
%|
%|	trick: you can make orbit=0 and orbit_start = column vector (length na)
%|	if you need nonuniformly spaced projection view angles.
%|
%| Copyright 2006-1-18, Jeff Fessler and Jang-Hwan Cho, University of Michigan
%|
%| 2009-12-04 modified source_zs definition to use source_z0, eliminate na/2
%| 2012-12-09 modified by Greg Handy to include the rebinned variable

if nargin == 1 && streq(type, 'test'), ct_geom_test, return, end
if nargin < 1, help(mfilename), error(mfilename), end

if streq(type, 'ge1') % special case: GE fan-beam
	st = ct_geom_ge1(type, varargin{:});
return
end

if streq(type, 'ge2') % special case: GE axial or helical
	st = ct_geom_ge2(type, varargin{:});
return
end

if streq(type, 'hd1') % special case: GE HD (UM only)
	st = ct_geom_hd1(type, varargin{:});
return
end

% defaults
st.type = type;
st.rebinned = false;
st.ns = [];
st.nt = [];
st.na = [];
st.down = 1;
st.nframe = 1; % entire scan as 1 "frame"
st.frame = 1;
st.orbit_start = 0;
st.pitch = 0; % default for axial
st.source_z0 = 0; % default for axial
st.units = 'mm';
st.user_source_zs = [];

if streq(type, 'fan')
	st = ct_geom_fan(st, varargin{:});
%elseif streq(type, 'par')
%	st = ct_geom_par(st, varargin{:});
%elseif streq(type, 'moj')
%	st = ct_geom_moj(st, varargin{:});
else
	fail('unknown sinotype %s', type)
end

if isempty(st.na), st.na = 2 * floor(st.ns * pi/2 / 2); end

ct_geom_checks(st)

meth = { ...
	's', @ct_geom_s, '()'; ...
	't', @ct_geom_t, '()'; ...
	'ws', @ct_geom_ws, '()'; ...
	'wt', @ct_geom_wt, '()'; ...
	'ad', @ct_geom_ad, '()'; ...
	'ar', @ct_geom_ar, '()'; ...
	'gamma', @ct_geom_gamma, '()'; ...
	'gamma_max', @ct_geom_gamma_max, '()'; ...
	'zfov', @ct_geom_zfov, '()'; ...
	'source_dz_per_view', @ct_geom_source_dz_per_view, '()'; ...
	'source_zs', @ct_geom_source_zs, '() -> [na]'; ...
	'orbit_short', @ct_geom_orbit_short, '()'; ...
	'xds', @ct_geom_xds, '()'; ...
	'yds', @ct_geom_yds, '()'; ...
	'downsample', @ct_geom_downsample, '()'; ...
	'dim', @ct_geom_dim, '()'; ...
	'ones', @ct_geom_ones, '()'; ...
	'rmax', @ct_geom_rmax, '()'; ...
	'footprint_size', @ct_geom_footprint_size, '(ig)'; ...
	'unitv', @ct_geom_unitv, '() | (is,it,ia)'; ...
	'zeros', @ct_geom_zeros, '()'; ...
	'shape', @ct_geom_shape, '()'; ...
	'plot', @ct_geom_plot, '() | (ig)';
	'plot3', @ct_geom_plot3, '() | (ig)';
	};

st = strum(st, meth);

if streq(st.source_z0, 'center')
	st.source_z0 = -st.na/2 * st.source_dz_per_view;
end

if any(st.down ~= 1)
	down = st.down; st.down = 1; % trick
	st = st.downsample(down);
end

if streq(type, 'fan') && streq(st.orbit, 'short')
	st.orbit = st.orbit_short / st.nframe; % jc
	st.orbit_start = st.orbit_start + (st.frame - 1) * st.orbit; % jc/jf
	if st.frame ~= 1, warn 'untested and possibly incorrect', end % jf
end

% ct_geom_orbit_short()
function os = ct_geom_orbit_short(st)
os = 180 + 2 * rad2deg(st.gamma_max);

% ct_geom_dim()
function dim = ct_geom_dim(st)
dim = [st.ns st.nt st.na];
if isempty(st.ns) || isempty(st.nt) || isempty(st.na)
	error 'dim requested without ns,nt,na'
end

% ct_geom_ones()
% sinogram of all ones
function out = ct_geom_ones(st)
out = ones(st.dim, 'single');

% ct_geom_zeros()
% sinogram of all zeros
function out = ct_geom_zeros(st)
out = zeros(st.dim, 'single');

% ct_geom_unitv()
% sinogram with just one ray
function out = ct_geom_unitv(st, is, it, ia)
out = st.zeros;
if ~isvar('is') || isempty(is)
	is = floor(st.ns/2 + 1);
	it = floor(st.nt/2 + 1);
	ia = 1;
end
out(is,it,ia) = 1;

% ct_geom_rmax()
% max radius within fov
function rmax = ct_geom_rmax(st)
smax = max(abs(st.s));
if streq(st.type, 'fan')
	if isinf(st.dso) % parallel
		rmax = smax;
	elseif st.dfs == 0 % arc
		rmax = st.dso * sin(smax / st.dsd);
	elseif isinf(st.dfs) % flat
		rmax = st.dso * sin(atan(smax / st.dsd));
	else
		error 'unknown case'
	end
end

% ct_geom_ws()
% 'middle' sample position
function ws = ct_geom_ws(st)
ws = (st.ns-1)/2 + st.offset_s;

% ct_geom_wt()
% 'middle' sample position
function wt = ct_geom_wt(st)
wt = (st.nt-1)/2 + st.offset_t;

% ct_geom_s()
% sample locations ('radial')
function s = ct_geom_s(st, varargin)
s = st.ds * ([0:st.ns-1]' - st.ws);
if length(varargin)
	s = s(varargin{:});
end

% ct_geom_t()
% sample locations ('radial')
function t = ct_geom_t(st, varargin)
t = st.dt * ([0:st.nt-1]' - st.wt);
if length(varargin)
	t = t(varargin{:});
end

% ct_geom_ad()
% angular sample locations (degrees)
function ang = ct_geom_ad(st, varargin)
ang = [0:st.na-1]'/st.na * st.orbit + st.orbit_start;
ang = ang(varargin{:});

% ct_geom_ar()
% angular sample locations (radians)
function ang = ct_geom_ar(st, varargin)
ang = deg2rad(ct_geom_ad(st));
ang = ang(varargin{:});

% ct_geom_shape()
% reshape into sinogram array
function sino = ct_geom_shape(st, sino)
sino = reshapee(sino, st.ns, st.nt, st.na, []);

% ct_geom_downsample()
% down-sample (for testing)
function st = ct_geom_downsample(si, down)
st = si;
st.down = st.down .* down;

switch numel(down)
case 1
	down_s = down;
	down_t = down;
	down_a = down;
case 3
	down_s = down(1);
	down_t = down(2);
	down_a = down(3);
otherwise
	fail('bad down %g', down)
end
clear down

st.ns = 4 * ceil(st.ns / down_s / 4); % multiple of 4 for simd
st.nt = 2 * ceil(st.nt / down_t / 2); % keep it even

user_specified = false; % did user specify zs or orbit_start vector?

if ~isempty(st.user_source_zs)
	st.user_source_zs = st.user_source_zs(1:down_a:st.na);
	user_specified = true;
end

if numel(st.orbit_start) > 1
	st.orbit_start = st.orbit_start(1:down_a:si.na);
	user_specified = true;
end

if streq(st.type, 'fan')
	st.ds = st.ds * down_s;
	st.dt = st.dt * down_t;
else
	fail('unknown sinotype "%s"', type)
end

if user_specified % if user-specified then just decimate
	st.na = length([1:down_a:st.na]);

elseif all(diff(si.source_zs) == 0) % axial
	st.na = max(1, round(si.na / down_a)); % at least one view

else % helical described by pitch and source_z0
	nturn = si.orbit / 360;
	na1 = si.na / nturn; % # of views in 1 turn, originally

	% if it is helical with equal view spacing, then adjust orbit slightly,
	% to preserve integer number of views per turn (if applicable)
	tol = 1e-4;
	if abs(round(na1) - na1) < tol % integer views per turn
		na2 = round(na1 / down_a); % new # of views in 1 turn (integer)
		tmp = nturn * na2; % new na value, roughly na/down
		if abs(round(tmp) - tmp) < tol
			st.na = round(tmp);
		else % adjust orbit
			st.na = floor(tmp);
			st.orbit = 360 / na2 * st.na;
		end
	else
		st.na = round(si.na / down_a);
	end
end

% ct_geom_zfov()
% axial FOV, considering magnification factor, as collimated at iso
function zfov = ct_geom_zfov(st)
if isinf(st.dso) || isinf(st.dsd) % parallel beam
	zfov = st.nt * st.dt;
else % cone
	zfov = st.dso / st.dsd * st.nt * st.dt;
end

% ct_geom_source_dz_per_view()
% source axial travel for each view
function out = ct_geom_source_dz_per_view(st)
if ~isempty(st.user_source_zs), fail 'undefined', end
if st.na == 1 % for tomosynthesis
	out = 0;
	return
end
if ~st.pitch
	out = 0;
	return
end
if length(st.orbit) ~= 1 || st.orbit == 0
	fail 'dz_per_view undefined for user angles'
end
na_per_360 = st.na * (360 / st.orbit); % # views per turn
out = st.pitch * st.zfov / na_per_360;

% ct_geom_checks
% check vectors user_source_zs and orbit_start
function ct_geom_checks(st)
if ~isempty(st.user_source_zs) && (st.pitch ~= 0)
	fail('only one of "pitch" (recommended) or user_source_zs can be given')
end
if ~isempty(st.user_source_zs) && length(st.user_source_zs) ~= st.na
	fail('user_source_zs size mismatch')
end
if numel(st.orbit_start) > 1 && numel(st.orbit_start) ~= st.na
	fail('orbit_start vector length')
end

% ct_geom_source_zs()
% source z locations
function source_zs = ct_geom_source_zs(st, varargin)

if ~isempty(st.user_source_zs) % user-provided
	source_zs = st.user_source_zs;
	source_zs = source_zs(varargin{:});
else
	source_dz = st.source_dz_per_view;
	source_zs = st.source_z0 + [0:st.na-1]' * source_dz;
	source_zs = source_zs(varargin{:});
end

% ct_geom_gamma()
% gamma sample values
function gamma = ct_geom_gamma(st, varargin)
switch st.dfs
case 0
	gamma = st.s / st.dsd; % 3rd gen: equiangular arc
case inf
	gamma = atan(st.s / st.dsd); % flat
otherwise
	error 'dfs not done'
end
gamma = gamma(varargin{:});

% ct_geom_gamma_max()
function gamma_max = ct_geom_gamma_max(st)
gamma_max = max(st.gamma);

% ct_geom_xds()
% center positions of detectors
function xds = ct_geom_xds(st, varargin)
switch st.type
case 'par'
	xds = st.s;
case 'fan'
	switch st.dfs
	case 0 % arc
		xds = st.dsd * sin(st.gamma);
	case inf % flat
		xds = st.s;
	otherwise
		error 'not done'
	end
otherwise
	error 'bug'
end

% ct_geom_yds()
% center positions of detectors
function yds = ct_geom_yds(st, varargin)
switch st.type
case 'par'
	yds = zeros(size(st.s));

case 'fan'
	switch st.dfs
	case 0 % arc
		yds = st.dso - st.dsd * cos(st.gamma);
	case inf % flat
		yds = -st.dod * ones(size(st.s));
	otherwise
		error 'not done'
	end
otherwise
	error 'bug'
end

% ct_geom_footprint_size()
% biggest possible footprint over image
function footprint_size = ct_geom_footprint_size(st, ig, varargin)
di = sqrt(ig.dx^2 + ig.dy^2);
smax = max(abs(st.s));
rfov = max(ig.nx * ig.dx, ig.ny * ig.dy) / 2;
dso = st.dso;
dsd = st.dsd;
if rfov > 0.99 * dso, error 'bad dso', end

if isinf(st.dso) % parallel
	footprint_size = di;
elseif st.dfs == 0 % arc3
	footprint_size = di * dsd / (dso - rfov);
elseif isinf(st.dfs) % flat
	footprint_size = di * sqrt(dsd^2 + smax^2) / (dso - rfov);
else
	error 'unknown case'
end
footprint_size = footprint_size / st.ds; % unitless

%
% ct_geom_fan()
%
function st = ct_geom_fan(st, varargin);

% defaults
st.orbit = 360; % [degrees]
st.ds		= 1;
st.dt		= [];
st.offset_s	= 0;
st.offset_t	= 0;

st.dsd = [];	% dis_src_det
st.dso = [];	% dis_src_iso
st.dod = [];	% dis_iso_det
st.dfs = 0;	% dis_foc_src (3rd gen CT)

subs = { ...
	'src_det_dis', 'dsd';
	'dis_src_det', 'dsd';
	'dis_src_iso', 'dso';
	'dis_iso_det', 'dod';
	'dis_foc_src', 'dfs';
	'nbeta', 'na';
	'source_zs', 'user_source_zs';
	};
st = vararg_pair(st, varargin, 'subs', subs);

% work out distances
if (~isempty(st.dsd) && isinf(st.dsd)) ...
|| (~isempty(st.dso) && isinf(st.dso)) % handle parallel-beam case gracefully
	st.dsd = inf; st.dso = inf; st.dod = 1;
end
if isempty(st.dsd) + isempty(st.dso) + isempty(st.dod) > 1
	error 'must provide at least two of dsd, dso, dod'
end
if isempty(st.dsd), st.dsd = st.dso + st.dod; end
if isempty(st.dso), st.dso = st.dsd - st.dod; end
if isempty(st.dod), st.dod = st.dsd - st.dso; end
if st.dso + st.dod ~= st.dsd
	error 'bad fan-beam distances'
end

if isempty(st.dt), st.dt = -st.ds; end
if iscell(st.dt)
	if length(st.dt) == 2 && streq(st.dt{1}, 'dz')
		dz = st.dt{2};
		st.dt = dz * st.dsd / st.dso;
		printm('dt = dz * dsd / dso = %g * %g / %g = %g', ...
			dz, st.dsd, st.dso, st.dt)
	else
		fail 'bad dt cell'
	end
end

% ct_geom_plot2()
% picture of 2D source position / detector geometry
function ct_geom_plot2(st, ig, varargin)
arg.tomosyn = false; % set to 1 for offset_x plotting trick
arg = vararg_pair(arg, varargin);
if arg.tomosyn, fail 'not done', end
if ~streq(st.type, 'fan'), error 'only fan done', end
x0 = 0;
y0 = st.dso;
t = linspace(0,2*pi,1001);

if isinf(st.dsd) % parallel beam
	rfov = max(abs(st.s));
	plot(	0, 0, '.', ...
		rfov * cos(t), rfov * sin(t), 'm:', ...
		st.s, -rfov, 'yo')

else % fan beam

	rot = deg2rad(st.orbit_start);
	rot = [cos(rot) sin(rot); -sin(rot) cos(rot)];
	p0 = rot * [x0; y0];
	pd = rot * [st.xds'; st.yds'];
	rfov = st.dso * sin(max(abs(st.gamma)));

	plot(	p0(1), p0(2), 'ys', ... % source
		st.dso * cos(t), st.dso * sin(t), 'c--', ... % source circle
		pd(1,:), pd(2,:), 'yo') % detector elements

	hold on
	if isvar('ig') && ~isempty(ig)
		xmin = min(ig.x); xmax = max(ig.x);
		ymin = min(ig.y); ymax = max(ig.y);
		im(ig.x, ig.y, ig.mask(:,:,1))

		plot([xmax xmin xmin xmax xmax], [ymax ymax ymin ymin ymax], 'g-')
	end

	plot(	0, 0, '.', ...
		[pd(1,1) p0(1) pd(1,end)], [pd(2,1) p0(2) pd(2,end)], 'r-', ...
		rfov * cos(t), rfov * sin(t), 'm:')

	axis equal, axis tight
	hold off

end

if isvar('ig') && ~isempty(ig)
	hold on
	xmin = min(ig.x); xmax = max(ig.x);
	ymin = min(ig.y); ymax = max(ig.y);
	plot([xmax xmin xmin xmax xmax], [ymax ymax ymin ymin ymax], 'g-')
	hold off
end
titlef('fov = %g', rfov)
axis square, zoom on

% ct_geom_plot3()
% picture of 3D helical CT geometry
function dummy = ct_geom_plot3(st, ig)
dummy = [];
if ~streq(st.type, 'fan'), error 'only fan done', end

if isinf(st.dso) || isinf(st.dsd)
	warn 'parallel-beam plot not done'
return
end

t1 = -st.dso * sin(st.ar); % source x pos
t2 = st.dso * cos(st.ar); % source y pos
t3 = st.source_zs; % source z pos
plot3(t1, t2, t3, 'y.') % source trajectory
view(-110, 22)
hold on
for ia = [1 1+floor(st.na/2) st.na]
	text(t1(ia), t2(ia), t3(ia), sprintf('%d', ia))
end

% axes
r100 = @(x) 100 * ceil(max(abs(x)) / 100);
r10 = @(x) 10 * ceil(max(abs(x)) / 10);
plot3([-1 1] * r100(t1), [0 0], [0 0])
text(r100(t1), 0, 'x')
plot3([0 0], [-1 1] * r100(t2), [0 0])
text(0, r100(t2), 'y')
plot3([0 0], [0 0], [-1 1] * r10(t3))
text(0, 0, r10(t3), 'z')

grid on
xlabelf('x (%s)', st.units)
ylabelf('y (%s)', st.units)
zlabelf('z (%s)', st.units)

if isvar('ig') && ~isempty(ig) && isfield(ig.meth,'z') % image box
	xmin = min(ig.x); xmax = max(ig.x);
	ymin = min(ig.y); ymax = max(ig.y);
	zmin = min(ig.z); zmax = max(ig.z);
	plot3(	[xmin xmax xmax xmin xmin ], ...
		[ymin ymin ymax ymax ymin ], ...
		[zmin zmin zmin zmin zmin ], ...
		'g-')
	plot3(	[xmin xmax xmax xmin xmin], ...
		[ymin ymin ymax ymax ymin], ...
		[zmax zmax zmax zmax zmax], ...
		'g-')
	plot3(xmin*[1 1], ymin*[1 1], [zmin zmax], 'g-')
	plot3(xmax*[1 1], ymin*[1 1], [zmin zmax], 'g-')
	plot3(xmin*[1 1], ymax*[1 1], [zmin zmax], 'g-')
	plot3(xmax*[1 1], ymax*[1 1], [zmin zmax], 'g-')
	plot3(xmin*ones(ig.nz,1), ymin*ones(ig.nz,1), ig.z, 'g.')
end

% detector array for source at each position
detcolor = {'r.', 'c.', 'm.'};
ia_list = [1 1+floor(st.na/2) st.na];
for ia = ia_list
	src = [t1(ia) t2(ia) t3(ia)];
	unit_vec = [-src(1) -src(2) 0] / sqrt(src(1)^2+src(2)^2);
	% diametrically opposite point of the source on the detector
	det_cen_loc = src + st.dsd * unit_vec;
	plot3(det_cen_loc(1), det_cen_loc(2), det_cen_loc(3), 'rs')

	rot = st.ar(ia);
	rot = [cos(rot) -sin(rot); sin(rot) cos(rot)];
	pd = rot * [st.xds'; st.yds'];

	for it=1:st.nt
		plot3(pd(1,:), pd(2,:), src(3) + st.t(it)*ones(1,st.ns), ...
			detcolor{find(ia == ia_list)})
	end

	if 1 && ia == 1 % line from source to middle of first detector row
		plot3([src(1) det_cen_loc(1)], [src(2) det_cen_loc(2)], ...
			[src(3) st.t(1)], 'm-')
	end
end
hold off

% ct_geom_plot()
function out = ct_geom_plot(st, varargin)
im clf
if all(st.source_zs == 0) % 2D or axial case
	ct_geom_plot2(st, varargin{:});
else
	ct_geom_plot3(st, varargin{:});
end
if nargout, out = []; end

% ct_geom_ge1()
% 'lightspeed';
% these numbers are published in IEEE T-MI Oct. 2006, p.1272-1283 wang:06:pwl
function st = ct_geom_ge1(type, varargin)

tmp.orbit = 360;
tmp = vararg_pair(tmp, varargin, 'allow_new', true);
if streq(tmp.orbit, 'short')
	na = 642; % trick: reduce na for short scans
	for ii=1:length(varargin)
		if (streq(varargin{ii}, 'orbit'))
			varargin{ii+1} = na/984*360;
			break
		end
	end
else
	na = 984;
end

st = ct_geom('fan', ...
	'ns', 888, ... % detector channels
	'nt', 1, ... % detector rows
	'na', na, ... % angular samples
	'orbit', 360, ...
	'offset_s', 1.25, ... % quarter-detector offset
	'dsd', 949.075, ...
	'dso', 541, ...
	'dfs', 0, ... % arc
	'ds', 1.0239, ... % detector pitch
	'dt', 1.0964, ... % detector row spacing for 0.625mm slices, 2009-12-06
	varargin{:});
%	'dod', 408.075, ...
% 'strip_width', [], ...

% ct_geom_ge2()
% helical CT
function st = ct_geom_ge2(type, varargin)
st = ct_geom_ge1(type, ...
	'nt', 64, ... % 64-slice
	'dt', 949.0750 / 541 * 0.625, ... % about 1.0964
	varargin{:});

% ct_geom_test()
function ct_geom_test
% axial cone-beam test
cg = ct_geom('fan', 'ns', 888, 'nt', 64, 'na', 984, ...
	'offset_s', 1.25, ...
	'dsd', 949, 'dod', 408);
cg.ad(2);
cg.downsample(2);
cg.rmax;
cg.ws;
cg.nframe; cg.frame; cg.gamma; cg.gamma_max;
cg.s(cg.ns/2+1);
if im
	cg.plot;
prompt
end

% test user_source_zs
cg = ct_geom('fan', 'dsd', 949, 'dod', 408, ...
	'ns', 888, 'ds', 1.0239, 'offset_s', 1.25, ...
	'nt', 8, 'dt', 1.0964, ...
	'na', 2*984, 'orbit', 2*360, ...
	'source_zs', []);
cg.source_zs(1:2:7);

% helical cone-beam test
cg = ct_geom('fan', 'dsd', 949, 'dod', 408, ...
	'ns', 888, 'ds', 1.0239, 'offset_s', 1.25, ...
	'nt', 8, 'dt', 1.0964, ...
	'na', 2*984, 'orbit', 2*360, ...
	'pitch', 1.0);

cg.source_zs(1:2:7);
cg.downsample(2);
cg.s(cg.ns/2+1);
if im
	cg.plot;
prompt
end

if 1 % footprint
	cg = ct_geom('ge2');
	ig = image_geom('nx', 512, 'fov', 500);
	cg.footprint_size(ig);
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function gX = mgaussianNoise3dPlot(noise, X)

% MGAUSSIANNOISE3DPLOT Draws a 3D or contour plot for the MGAUSSIAN noise model.
% FORMAT
% DESC draws a 3D or contour plot for the multiple output Gaussian noise model.
% ARG noise : the noise structure for which the plot is required.
% ARG plotType : string containing the name of the plotting function (for example mesh, contour).
% ARG X : the input X data in the form of a 'mesh' matrix.
% ARG Y : the input Y data in the form of a 'mesh' matrix.
% ARG mu : the input mean in the form of a 'mesh' matrix.
% ARG varSigma : the input variance in the form of a 'mesh' matrix.
% ARG P1, P2, P3 ... : optional additional arguments for the given plot type.
% RETURN h : the gradients of the diagonal with respect to each element
% of X. The returned matrix has the same dimensions as X.
%
% SEEALSO : mgaussianNoiseParamInit, noise3dPlot,
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% NOISE

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%TestRunMonitor Abstract base class for monitoring a running test suite
%   The abstract TestRunMonitor class defines an object that can observe and
%   record the results of running a test suite.  The run() method of a
%   TestComponent object takes a TestRunMonitor object as an input argument.
%
%   Different test suite logging or reporting functionality can be achieved by
%   subclassing TestRunMonitor.  For example, see the TestRunLogger and the
%   CommandWindowTestRunDisplay classes.
%
%   TestRunMonitor methods:
%       TestRunMonitor        - Constructor
%       testComponentStarted  - Called at beginning of test component run
%       testComponentFinished - Called when test component run finished
%       testCaseFailure       -   Called when a test case fails
%       testCaseError         - Called when a test case causes an error
%
%   See also CommandWindowTestRunDisplay, TestRunLogger, TestCase, TestSuite

%   Steven L. Eddins
%   Copyright 2008-2009 The MathWorks, Inc.

classdef TestRunMonitor < handle

    methods (Abstract)

        testComponentStarted(self, component)

        testComponentFinished(self, component, did_pass)

        testCaseFailure(self, test_case, failure_exception)

        testCaseError(self, test_case, error_exception)

    end
end

author: ricedsp


text:
function b = subsref(a,s)

Q = 'a'; Q2='a';

%% ---- preprocess adding children calls each time there isn't one!!! --
%%  to minimic f.child{1}.child{1}  by just doing f{1}{1} !!

t=s; j=1; isfile=isobject(a); lasttext=0;
for i=1:length(s)

  if ischar(s(i).subs)  %% check if last thing accessed is an object
    isfile=logical(myexist(s(i).subs,'file')==2);
   end

  if equal(s(i).type,'{}') | equal(s(i).type,'()')  %% time for child!
    before=i-1; if i==1 before=1; end; %% look at s(before)
    if ~equal(s(before).subs,'child') & isfile
      % if it exists as file is hopefully an object, isnt child add one!
      t(j)=s(1); t(j).subs='child'; t(j).type='.';
      j=j+1;
    end
  end

  if ischar(s(i).subs)
    lasttext=j;
  end

  t(j)=s(i); j=j+1;
end
s=t;
%%-------------------------------------------------------------

for i=1:length(s)
  if equal(s(i).type, '.')
    Q=[Q, '.', s(i).subs];

  elseif equal(s(i).type, '()')
    tmp = length(s(i).subs);
    Q=[Q, '('];
    for j=1:(tmp-1)
      Q=[Q, 's(', num2str(i), ').subs{', num2str(j), '}, '];
    end;
    Q=[Q, 's(', num2str(i), ').subs{', num2str(tmp), '})'];

  elseif equal(s(i).type, '{}')
    tmp = length(s(i).subs);
    Q=[Q, '{'];
    for j=1:(tmp-1)
      Q=[Q, 's(', num2str(i), ').subs{', num2str(j), '}, '];
    end;
    Q=[Q, 's(', num2str(i), ').subs{', num2str(tmp), '}}'];
  end ;

  if i<lasttext Q2=Q; end;
end;

%% ---- if member doesn't exist search in the children,algorithm.. --------
%% -- differs from subsasgn in that it only returns first one   ----------

global recursive_subsasgn_off;
if isempty(recursive_subsasgn_off) recursive_subsasgn_off=0;end;

if ~isfile &  (~recursive_subsasgn_off)
  name=s(lasttext).subs;  % this is the member we are searching for
  Qend=Q(length(Q2)+length(name)+2:length(Q));

  notfound=1;
  stack=[]; stack{1}=Q2;
  while notfound %% keep searching into children..

      cur=stack{1};
      b=eval(cur);   % eval top line
      if isfield(struct(b),name) notfound=0; end; % check for  member

      if notfound                      %% check for aliases
	if ~isfield(struct(b),'alias') % as long asn't currently in @algorithm
	  ali=b.algorithm.alias;
	  for i=1:2:length(ali)
	    if strcmp(name,ali{i})
	      if isfield(struct(b),ali{i+1}) %% found it!
		name=ali{i+1}; notfound=0;
	      end
	    end
	  end
	end
      end

      if notfound
	if isfield(struct(b),'algorithm') % check algorithm path
	  stack{length(stack)+1}=[ cur '.algorithm'];
	end
	if isfield(struct(b),'child') % check child path
	  if ~iscell(b.child)
	    stack{length(stack)+1}=[ cur '.child'];
	  else
	    for i=1:length(b.child)
	      stack{length(stack)+1}=[ cur '.child{' num2str(i) '}'];
	    end
	  end
	end

	stack=stack(2:length(stack));   % remove top of stack
      end
  end
  Q=[stack{1} '.' name Qend];

end

%%-------------------------------------------------------------

b=eval(Q) ;  %% evaluate final answer

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function varargout=skew3d(varargin)
%
% For calling details please see v_skew3d.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_skew3d.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: skew3d.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_skew3d(varargin{:});
else
    v_skew3d(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function v = cvxtuple( varargin )

if nargin == 1,
    v = varargin{1};
    switch class( v ),
        case 'struct',
            if numel( v ) ~= 1,
                cvx_throw( 'struct arrays not permitted in cvx tuple objects.' );
            end
        case 'cell',
            v = reshape( v, 1, numel( v ) );
        otherwise,
            return
    end
else
    v = varargin;
end
v = class( struct( 'value_', { v }, 'dual_', { [] } ), 'cvxtuple' );

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function handles = switchMode(handles, newMode)
%SWITCHMODE   Go from one mode of CHEBGUI to another.
% Calling sequence:
%   HANDLES = SWITCHMODE(HANDLES, NEWMODE)
% where
%   HANDLES:    A MATLAB handle object for the CHEBGUI figure.
%   NEWMODE:    Mode we want to switch to.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Do a lot of disabling/enabling and hiding/showing objects on the CHEBGUI
% figure.
if ( strcmp(newMode, 'bvp') ) % Going into BVP mode
    handles.guifile.type = 'bvp';

    set(handles.button_bvp, 'Value', 1)
    set(handles.button_ivp, 'Value', 0)
    set(handles.button_pde, 'Value', 0)
    set(handles.button_eig, 'Value', 0)

    set(handles.panel_DEs, 'Title', 'Differential equations')
    set(handles.toggle_useLatest, 'Visible', 'on')

    set(handles.panel_timedomain, 'Visible', 'off')

    set(handles.button_realplot, 'Visible', 'off')
    set(handles.button_imagplot, 'Visible', 'off')
    set(handles.popupmenu_bottomFig, 'Visible', 'on')
    set(handles.toggle_useLatest, 'Value', 0)

    if ( isempty(fieldnames(handles.latest)) )
        set(handles.toggle_useLatest, 'Enable', 'off')
    end

    % Change heading for BC/IC input field:
    set(handles.panel_BCs, 'Title', 'Boundary conditions')

    % Hide the iter_list box:
    handles = hideIterList(handles);

    % Show the initial guess panel, and ensure it has the correct title:
    set(handles.panel_initialGuess, 'Visible', 'on', 'Title', 'Initial guess');

    set(handles.panel_eigopts, 'Visible', 'Off')

    % Switch from LBC and RBC fields to BC
    handles = toggleBCinput(handles, 'off', 'on');

    % Change the list of available options.
    % Enable ODE menu options
    set(handles.menu_odedampednewton, 'Enable', 'On')
    set(handles.menu_odeplotting, 'Enable', 'On')
    % Disable IVP solver option
    set(handles.menu_ivpSolver, 'Enable', 'off');
    % Disable PDE menu options
    set(handles.menu_pdeSolver, 'Enable', 'Off')
    set(handles.menu_pdeplotting, 'Enable', 'Off')
    set(handles.menu_pdeholdplot, 'Enable', 'Off')
    set(handles.menu_pdefix, 'Enable', 'Off')
    set(handles.menu_fixN, 'Enable', 'Off')

    % Make the discretization visible
    set(handles.panel_discretization, 'Visible', 'on');

    % Hide the IVP solver panel visible
    set(handles.panel_IVPsolver, 'Visible', 'off')

    % Clear the top figure
    chebguiController.initialiseFigureTop(handles)

    % Always clear the bottom figure
    chebguiController.initialiseFigureBottom(handles)

    % Make the discretization panel visible
    set(handles.panel_discretization, 'Visible', 'on')

    % Set correct header to top figure
    set(handles.panel_figSol, 'Title', 'Solution')

    % Note: The line below is a nice way to disable to objects, rather than hide
    % them.
    % Enable all childrens of discretisation panel
    %     set(findall(handles.panel_discretization,  '-property',  'enable'),  ...
    %         'visible',  'on')
elseif ( strcmp(newMode, 'ivp') ) % Going into IVP mode
    % IVP mode is almost identical to BVP, so do the same adjustments, and then
    % minor corrections at the end:
    handles = chebguiController.switchMode(handles, 'bvp');

    % Set the type of the problem
    handles.guifile.type = 'ivp';

    % Change heading for BC/IC input field:
    set(handles.panel_BCs, 'Title', 'Initial/final conditions')

    % Hide the discretization panel
    set(handles.panel_discretization, 'Visible', 'off');

    % Make the IVP solver panel visible
    set(handles.panel_IVPsolver, 'Visible', 'on')

    set(handles.button_bvp, 'Value', 0)
    set(handles.button_ivp, 'Value', 1)
    % Enable IVP solver option
    set(handles.menu_ivpSolver, 'Enable', 'on');

    % If we're solving in time-stepping mode, hide the initial guess box:
    if ( ~isempty(strfind(handles.guifile.options.ivpSolver, 'ode')) )
        set(handles.panel_initialGuess,'Visible','off');
    end

elseif ( strcmp(newMode, 'pde') ) % Going into PDE mode
    handles.guifile.type = 'pde';

    set(handles.button_bvp, 'Value', 0)
    set(handles.button_ivp, 'Value', 0)
    set(handles.button_pde, 'Value', 1)
    set(handles.button_eig, 'Value', 0)

    % Hide the iter_list box:
    handles = hideIterList(handles);

    % Show the initial guess panel, and ensure it has the correct title:
    set(handles.panel_initialGuess,'Visible','on','Title','Initial condition');

    set(handles.toggle_useLatest, 'Visible', 'off')
    set(handles.panel_DEs, 'Title', 'Differential equations')

    set(handles.input_GUESS, 'Enable', 'On')
    set(handles.toggle_useLatest, 'Value', 0)
    set(handles.toggle_useLatest, 'Enable', 'off')

    set(handles.panel_timedomain, 'Visible', 'on')

    set(handles.button_realplot, 'Visible', 'off')
    set(handles.button_imagplot, 'Visible', 'off')
    set(handles.popupmenu_bottomFig, 'Visible', 'off')

    set(handles.panel_eigopts, 'Visible', 'Off')

    % Switch from LBC and RBC fields to BC
    handles = toggleBCinput(handles, 'on', 'off');

    % Change the list of available options
    % Disable ODE menu options
    set(handles.menu_odedampednewton, 'Enable', 'Off')
    set(handles.menu_odeplotting, 'Enable', 'Off')
    % Disable IVP solver option
    set(handles.menu_ivpSolver, 'Enable', 'off');
    % Enable PDE menuoptions
    set(handles.menu_pdeSolver, 'Enable', 'On')
    set(handles.menu_pdeplotting, 'Enable', 'On')
    set(handles.menu_pdeholdplot, 'Enable', 'On')
    set(handles.menu_pdefix, 'Enable', 'On')
    set(handles.menu_fixN, 'Enable', 'On')

    % Clear the top figure
    chebguiController.initialiseFigureTop(handles)

    % Always clear the bottom figure
    chebguiController.initialiseFigureBottom(handles)

    % Hide the discretization and IVP solver panels
    set(handles.panel_discretization, 'Visible', 'off')
    set(handles.panel_IVPsolver, 'Visible', 'off')

    % Set correct header to top figure
    set(handles.panel_figSol, 'Title', 'Initial condition')

else % Going into EIG mode
    handles.guifile.type = 'eig';

    set(handles.button_bvp, 'Value', 0)
    set(handles.button_ivp, 'Value', 0)
    set(handles.button_pde, 'Value', 0)
    set(handles.button_eig, 'Value', 1)

    set(handles.toggle_useLatest, 'Visible', 'off')
    set(handles.panel_DEs, 'Title', 'Differential operator')
    set(handles.text_initial, 'String', 'Look for')
    set(handles.toggle_useLatest, 'Enable', 'off')

    set(handles.panel_timedomain, 'Visible', 'off')

    set(handles.button_realplot, 'Visible', 'on')
    set(handles.button_realplot, 'Value', 1)
    set(handles.button_imagplot, 'Visible', 'on')
    set(handles.button_imagplot, 'Value', 0)
    set(handles.popupmenu_bottomFig, 'Visible', 'off')

    set(handles.panel_eigopts,'Visible','On')
    if ( ~isempty(handles.guifile.sigma) )
        switch ( handles.guifile.sigma )
            case ''
                set(handles.popupmenu_sigma, 'Value', 1);
            case 'lm'
                set(handles.popupmenu_sigma, 'Value', 2);
            case 'sm'
                set(handles.popupmenu_sigma, 'Value', 3);
            case 'lr'
                set(handles.popupmenu_sigma, 'Value', 4);
            case 'sr'
                set(handles.popupmenu_sigma, 'Value', 5);
            case 'li'
                set(handles.popupmenu_sigma, 'Value', 6);
            case 'si'
                set(handles.popupmenu_sigma, 'Value', 7);
        end
    else
        set(handles.popupmenu_sigma, 'Value', 1)
    end

    if ( isfield(handles.guifile.options,'numeigs') ...
        && ~isempty(handles.guifile.options.numeigs) )
        set(handles.edit_eigN, 'String', handles.guifile.options.numeigs);
    else
        set(handles.edit_eigN, 'String', 6);
    end

    % Hide the initial guess panel, and ensure it has the correct title:
    set(handles.panel_initialGuess,'Visible','off');

    % Change heading for BC/IC input field:
    set(handles.panel_BCs, 'Title', 'Boundary conditions')

    % Hide the iter_list box:
    handles = hideIterList(handles);

    % Switch from LBC and RBC fields to BC
    handles = toggleBCinput(handles, 'off', 'on');

    % Change the list of available options.

    % Disable ODE options
    set(handles.menu_odedampednewton, 'Enable', 'Off')
    set(handles.menu_odeplotting, 'Enable', 'Off')

    % Disable IVP solver option
    set(handles.menu_ivpSolver, 'Enable', 'off');

    % Disable PDE options
    set(handles.menu_pdeSolver, 'Enable', 'Off')
    set(handles.menu_pdeplotting, 'Enable', 'Off')
    set(handles.menu_pdeholdplot, 'Enable', 'Off')
    set(handles.menu_pdefix, 'Enable', 'Off')
    set(handles.menu_fixN, 'Enable', 'Off')

    % Make the discretization panel visible
    set(handles.panel_discretization, 'Visible', 'on')

    % Hide the IVP solver panel
    set(handles.panel_IVPsolver, 'Visible', 'off')

    % Set correct header to top figure
    set(handles.panel_figSol, 'Title', 'Eigenvalues (imag vs real)')

    % Clear the top figure
    chebguiController.initialiseFigureTop(handles)

    % Always clear the bottom figure
    chebguiController.initialiseFigureBottom(handles)

end

end

function handles = toggleBCinput(handles, onOff1, onOff2)
set(handles.input_LBC, 'Visible', onOff1)
set(handles.panel_leftBCs, 'Visible', onOff1)
set(handles.input_RBC, 'Visible', onOff1)
set(handles.panel_rightBCs, 'Visible', onOff1)
set(handles.input_BC,  'Visible', onOff2)
set(handles.panel_BCs,  'Visible', onOff2)
end

function handles = hideIterList(handles)
set(handles.iter_list, 'Visible', 'off')
set(handles.iter_list, 'String', '')
set(handles.iter_list, 'Value', 0)
set(handles.iter_text, 'Visible', 'off')
set(handles.iter_text, 'String', '')
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function v = norm( F )
%NORM   Frobenius norm of a CHEBFUN2V.
%   V = NORM(F) returns the Frobenius norm of the two/three components, i.e.
%       V = sqrt(norm(F1).^2 + norm(F2).^2),
%   or
%       V = sqrt(norm(F1).^2 + norm(F2).^2 + norm(F3).^2) .

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Empty check:
if ( isempty( F ) )
    v = [];
    return
end

nF = F.nComponents;
v = 0;
for jj = 1:nF
    v = v + sum( svd( F.components{jj} ).^2 );
%     v = v + sum2( power( F.components{jj}, 2 ) );
end
v = sqrt(v);

end


text:
% Use the VAD stream (the second stream) to select speech frames in the first stream
% Assume the input data is a matrix of DxT size.
function output = SynchronizeStreamsVAD(data_in, action)
data = data_in{1};
vad = data_in{2};
[D,nFr] = size(data);
precision = class(gather(data_in{1}(1,1,1)));

words = ExtractWordsFromString_v2(action);
action = words{1};

% first
switch lower(action)
    case 'concatenation'
        [n1,n2,n3] = size(data);
        if n3>1
            nFr = n3;
        else
            nFr = n2;   % this would be erronous if the input is a tensor with just one frame. But this should be almost impossible as we will not use vad in this case.
        end
        vad(nFr+1:end) = [];
        if n3>1
            output = data(:,:,vad==1);
        else
            output = data(:,vad==1);
        end

    case 'segmentation'
        seglen = str2num(words{2});
        segshift = str2num(words{3});
        vad_seg = label2seg(vad);
        output_tmp = {};
        for j=1:length(vad_seg.label)
            if vad_seg.label(j)==0; continue; end
            if vad_seg.stop(j)-vad_seg.start(j) < seglen; continue; end
            start = max(1,vad_seg.start(j));
            stop = min(vad_seg.stop(j), nFr);
            output_tmp{end+1} = DivideSent2Segments(data(:,start:stop), seglen, segshift, 0);
        end
        % move data into a tensor
        nSeg = [];
        for j=1:length(output_tmp)
            nSeg(j) = size(output_tmp{j}, 3);
        end
        if length(nSeg) == 0
            pause(0.1);
        end
        output = zeros(D, seglen, sum(nSeg), precision);
        for j=1:length(output_tmp)
            output(:,:, sum(nSeg(1:j-1))+1 : sum(nSeg(1:j)) ) = output_tmp{j};
        end

    otherwise
        fprintf('SynchronizeStreamsVAD: Error: unknown action %s\n', action);
end

end

author: singaxiong

repo: SignalGraph