text:
function [p] = mesh_open(p)

% mesh_open - calls functions to read a triangulation file
%
% Usage: [p] = mesh_open(p)
%
% p is a parameter structure (see eeg_toolbox_defaults for
% more details). In this function, it should contain at least
% the following string fields:
%
%       p.mesh.path - the directory location of the file to load
%       p.mesh.file - the name of the file to load
%       p.mesh.type - the file format:
%
%      'emse'
%      'brainstorm'
%      'ascii' or 'freesurfer_ascii' for *.asc/.tri files
%      'freesurfer_surf' for freesurfer binary surface
%      'freesurfer_curv' for freesurfer binary curvature
%      'freesurfer_overlay' for freesurfer binary overlay (.w)
%
% The last two file types are associated with the current cortex,
% if it contains the same number of vertices (eg, load 'freesurfer_surf'
% first, then 'freesurfer_curv' or 'freesurfer_overlay').
%
% The file formats supported here are described in more detail
% at their respective websites:
%
%       FreeSurfer: http://surfer.nmr.mgh.harvard.edu
%       EMSE:       http://www.sourcesignal.com
%       BrainStorm: http://neuroimage.usc.edu/brainstorm
%
% The return structure creates or updates p.mesh.data, which
% contains cell arrays:
%
%       p.mesh.data.meshtype     type of surface, strings
%       p.mesh.data.vertices     Mx3 (x,y,z) vertices
%       p.mesh.data.faces        Mx3 vertex indices
%       p.mesh.data.Cdata        scalar overlay, M vert x N overlays
%
% For example, having loaded a scalp and an inner skull mesh,
% p.mesh.data.meshtype could be:
%
%       p.mesh.data.meshtype{1} = 'scalp'
%       p.mesh.data.meshtype{2} = 'inner skull'
%
% To plot the data returned, see mesh_plot or try
%
% Hpatch = patch('Vertices',p.mesh.data.vertices{1}',...
%                'Faces',p.mesh.data.faces{1},...
%                'Property','PropertyValue',...);
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:57 $

% Licence:  GNU GPL, no express or implied warranties
% History:  02/2002 Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf('\nMESH_OPEN...\n'); tic;

if ~exist('p','var'),
 [p] = eeg_toolbox_defaults;
  fprintf('...creating default p structure.\n');
end

[path,name,ext] = fileparts(strcat(p.mesh.path,filesep,p.mesh.file));
file = fullfile(path,[name ext]);

if exist(file) ~= 2,
  msg = sprintf('\n\n...file does not exist:\n\t%s\n',file);
  error(msg);
end

type = lower(p.mesh.type);

switch type,

  case 'emse',

    fprintf('...loading EMSE mesh from:\n\t%s\n',file);

    % Get EMSE .wfr data (in meters)
    options = {'vertex','face'};
    [vertices,faces,edges,meshtype] = mesh_emse2matlab(file,options);

    vertex_matrix = [vertices.x; vertices.y; vertices.z]';
    face_matrix = [faces.vertex1;faces.vertex2;faces.vertex3]';

    % Rotate 90 degrees around Z for EMSE data
    vertex_matrix = Rz(vertex_matrix,90,'degrees');

    % Is this a new or replacement mesh?
    [p.mesh.current,meshExists] = mesh_check(p,meshtype);

    p.mesh.data.meshtype{p.mesh.current}    = meshtype;
    p.mesh.data.vertices{p.mesh.current}    = vertex_matrix;
    p.mesh.data.faces{p.mesh.current}       = face_matrix;
    p.mesh.data.lapmat{p.mesh.current}      = [];
    p.mesh.data.lapint{p.mesh.current}      = [];
    p.mesh.data.timeseries{p.mesh.current}  = [];
    p.mesh.data.Cdata{p.mesh.current}       = [];

  case 'brainstorm',

    fprintf('...loading BrainStorm data from:\n\t%s\n',file);
    load(file);

    p.mesh.data = [];

    for i=1:size(Comment,2),
      if isempty(Comment{i}), continue; end
      fprintf('...converting tesselation: %s\n',Comment{i});
      p.mesh.data.meshtype{i}   = Comment{i};
      p.mesh.data.vertices{i}   = Vertices{i}';  % transpose Vertices
      p.mesh.data.faces{i}      = Faces{i};
      p.mesh.data.lapmat{i}     = [];
      p.mesh.data.lapint{i}     = [];
      p.mesh.data.timeseries{i} = [];
      p.mesh.data.Cdata{i}      = [];
      p.mesh.current = i;
    end

  case {'ascii','freesurfer_ascii'},

    fprintf('...loading ASCII or FreeSurfer data from:\n\t%s\n',file);

    % Get Freesurfer data
    if findstr(file,'.tri'),
      [vertex_matrix,face_matrix] = freesurfer_read_tri(file);
    else
      [vertex_matrix,face_matrix] = freesurfer_read_ascii(file);
    end

    fprintf('...converting surface coordinates (mm to meters)\n');
    vertex_matrix = vertex_matrix ./ 1000;

    meshtype = freesurfer_meshtype(file); % see function below

    [p.mesh.current,meshExists] = mesh_check(p,meshtype);

    p.mesh.data.meshtype{p.mesh.current}   = meshtype;
    p.mesh.data.vertices{p.mesh.current}   = vertex_matrix;
    p.mesh.data.faces{p.mesh.current}      = face_matrix;
    p.mesh.data.lapmat{p.mesh.current}     = [];
    p.mesh.data.lapint{p.mesh.current}     = [];
    p.mesh.data.timeseries{p.mesh.current} = [];
    p.mesh.data.Cdata{p.mesh.current}      = [];

  case 'freesurfer_surf',

    fprintf('...loading FreeSurfer binary surface from:\n\t%s\n',file);

    meshtype = freesurfer_meshtype(file); % see function below

    [p.mesh.current,meshExists] = mesh_check(p,meshtype);

    [vertex_matrix, face_matrix] = freesurfer_read_surf(file);

    fprintf('...converting surface coordinates (mm to meters)\n');
    vertex_matrix = vertex_matrix ./ 1000;

    p.mesh.data.meshtype{p.mesh.current}   = meshtype;
    p.mesh.data.vertices{p.mesh.current}   = vertex_matrix;
    p.mesh.data.faces{p.mesh.current}      = face_matrix;
    p.mesh.data.lapmat{p.mesh.current}     = [];
    p.mesh.data.lapint{p.mesh.current}     = [];
    p.mesh.data.timeseries{p.mesh.current} = [];
    p.mesh.data.Cdata{p.mesh.current}      = [];

  case 'freesurfer_curv',

    fprintf('...loading FreeSurfer binary curvature from:\n\t%s\n',file);

    [curv, Nfaces] = freesurfer_read_curv(file);

    allocated = 0;
    for meshN = 1:length(p.mesh.data.meshtype),

      type = p.mesh.data.meshtype{meshN};

      if size(p.mesh.data.vertices{meshN},1) == size(curv,1),
        fprintf('...allocating overlay to Cdata for ''%s'' surface.\n',type);
        p.mesh.data.Cdata{meshN} = curv;
        p.mesh.current = meshN;
        allocated = 1;
      end

    end
    if allocated < 1,
      fprintf('...failed to allocate curvature to any surface, incompatible Nfaces!\n');
    end

  case 'freesurfer_overlay',

    fprintf('...loading FreeSurfer binary overlay from:\n\t%s\n',file);

    [w,vert] = freesurfer_read_wfile(file);

    Nvert = length(vert); clear vert;

    allocated = 0;
    for meshN = 1:length(p.mesh.data.meshtype),

      type = p.mesh.data.meshtype{meshN};

      if size(p.mesh.data.vertices{meshN},1) == Nvert,
        fprintf('...allocating overlay to Cdata for ''%s'' surface.\n',type);
        p.mesh.data.Cdata{meshN} = w;
        p.mesh.current = meshN;
        allocated = 1;
      end

    end
    if allocated < 1,
      fprintf('...failed to allocate overlay to any surface, incompatible vertices!\n');
    end

  otherwise,
    fprintf('...mesh format: %s\n', p.mesh.type);
    fprintf('...sorry, cannot load this data format at present.\n');
    return;
end

t=toc; fprintf('...done (%5.2f sec).\n\n',t);

return

repo: OpenBMI


text:
function blkStruct = slblocks

blkStruct.Name = sprintf('DC Motor Library');
blkStruct.OpenFcn = 'dc_motor_lib';
blkStruct.MaskDisplay = 'disp(''DC Motor Library'')';

% Information for Simulink Library Browser
Browser(1).Library = 'dc_motor_lib';
Browser(1).Name    = 'DC Motor Library';
Browser(1).IsFlat  = 1;% Is this library "flat"

blkStruct.Browser = Browser;

author: Sable

repo: mcbench-benchmarks


text:
function   [dc,dk] = bspderiv(d,c,k)
%
% Function Name:
%
%   bspdeval - Evaluate the control points and knot sequence of the derivative
%              of a univariate B-Spline.
%
% Calling Sequence:
%
%   [dc,dk] = bspderiv(d,c,k)
%
% Parameters:
%
%   d	: Degree of the B-Spline.
%
%   c	: Control Points, matrix of size (dim,nc).
%
%   k	: Knot sequence, row vector of size nk.
%
%   dc	: Control points of the derivative
%
%   dk	: Knot sequence of the derivative
%
% Description:
%
%   Evaluate the derivative of univariate B-Spline, which is itself a B-Spline.
%   This function provides an interface to a toolbox 'C' routine.
[mc,nc] = size(c);
nk = numel(k);
                                                     %
                                                     % int bspderiv(int d, double *c, int mc, int nc, double *k, int nk, double *dc,
                                                     %              double *dk)
                                                     % {
                                                     %   int ierr = 0;
                                                     %   int i, j, tmp;
                                                     %
                                                     %   // control points
                                                     %   double **ctrl = vec2mat(c,mc,nc);
                                                     %
                                                     %   // control points of the derivative
dc = zeros(mc,nc-1);                                 %   double **dctrl = vec2mat(dc,mc,nc-1);
                                                     %
for i=0:nc-2                                         %   for (i = 0; i < nc-1; i++) {
   tmp = d / (k(i+d+2) - k(i+2));                    %     tmp = d / (k[i+d+1] - k[i+1]);
   for j=0:mc-1                                      %     for (j = 0; j < mc; j++) {
       dc(j+1,i+1) = tmp*(c(j+1,i+2) - c(j+1,i+1));  %       dctrl[i][j] = tmp * (ctrl[i+1][j] - ctrl[i][j]);
   end                                               %     }
end                                                  %   }
                                                     %
dk = zeros(1,nk-2);                                  %   j = 0;
for i=1:nk-2                                         %   for (i = 1; i < nk-1; i++)
   dk(i) = k(i+1);                                   %     dk[j++] = k[i];
end                                                  %
                                                     %   freevec2mat(dctrl);
                                                     %   freevec2mat(ctrl);
                                                     %
                                                     %   return ierr;
                                                     % }</pre>


text:
function rho12=rho12(h)
sldata;
limitpoints;
slopes;
rho12=rho1*(T12(h)/T1)^-(1+(g/(m12*R)));


text:
function [xRgb,yRgb,zRgb,xform,xAxes,yAxes,zAxes] = dtiGetCurSlices(handles)
% Determine current slices and return them (and potentially a lot of other info)
%
%  [xRgb,yRgb,zRgb,xform,xAxes,yAxes,zAxes] = dtiGetCurSlices(handles)
%
% The potential background and overlay images are stored in the background
% slots, dtiH.bg().  This is a structure in which the image is always
% between 0 and 1, but the values of the original are stored in the
% .bg(1).minVal and .bg(1).maxVal slots.
%
% Note: mmPerVox and xform are for the current space that the slices are
% in, NOT the original image space. This is typically 1mm isotropic, ac-pc
% aligned space.
%
% HISTORY:
%  ?????: Dougherty & Wandell wrote it.
%  2005.06.09 RFD: minor code optimizations to make refreshes closer to 'realtime'.
%  2006.11.08 RFD: default interpolation method is now nearest-neighbor
%    rather than trilinear. We should make this a use option.
%  2006.11.10 RFD: interpolation method is now an option.
%  2009.06.16 RFD: added 0.5 offset to xAxes,yAxes,zAxes to make the axes
%    tick marks line up with the voxel centers.
%
% Brian and Bob (c) Stanford VISTA Team 2005

%% Get key parameters from the GUI
curPosition   = dtiGet(handles,'acpcpos');
overlayThresh = dtiGet(handles,'cur overlay thresh');
overlayAlpha  = dtiGet(handles,'cur overlay alpha');
mmPerVox      = dtiGet(handles,'render mm');
bb            = dtiGet(handles,'default Bounding Box');
anat          = dtiGet(handles,'current anatomy data');      % anat      = handles.bg(curBgNum).img;
anatXform     = dtiGet(handles,'bg img2acpc xform');
dispRange     = dtiGet(handles,'display range');
% handles.bg(curBgNum).displayValueRange;
% dispRange     = dtiGet(handles,handles.bg(curBgNum).displayValueRange;

%% Here we get the anatomical slices transformed into image space
[zRgb,x,y,z] = dtiGetSlice(anatXform, anat, 3, curPosition(3), bb, handles.interpType, mmPerVox, dispRange);
zAxes = [x(1), x(end); y(1), y(end)]+0.5;
[yRgb,x,y,z] = dtiGetSlice(anatXform, anat, 2, curPosition(2), bb, handles.interpType, mmPerVox, dispRange);
yAxes = [x(1), x(end); z(1), z(end)]+0.5;
[xRgb,x,y,z] = dtiGetSlice(anatXform, anat, 1, curPosition(1), bb, handles.interpType, mmPerVox, dispRange);
xAxes = [y(1), y(end); z(1), z(end)]+0.5;
%  mrvNewGraphWin; imagesc(flipud(zRgb)); axis image; colormap(gray)

if(ndims(anat) == 3)
    % Turn the images into RGB
    xRgb = repmat(xRgb,[1 1 3]);
    yRgb = repmat(yRgb,[1 1 3]);
    zRgb = repmat(zRgb,[1 1 3]);
elseif(ndims(anat)==4)
    % 3d views can only handle grayscale anatomy, so we'll just use a
    % luminance map from the RGB data.
    % We are not sure this works all that well. - FP
    anat = mean(anat,4);
end
%  mrvNewGraphWin; image(zRgb); axis image; colormap(gray)

% Force the background image to be between 0 and 1
xRgb(xRgb<0) = 0; xRgb(xRgb>1) = 1;
yRgb(yRgb<0) = 0; yRgb(yRgb>1) = 1;
zRgb(zRgb<0) = 0; zRgb(zRgb>1) = 1;

if(overlayAlpha>0)
    %[overlayImg, oMmPerVoxel, oXform] = dtiGetCurAnat(handles,1);
    n          = dtiGet(handles,'overlay number');
    overlayImg = dtiGet(handles,'overlay image',n);
    oXform     = dtiGet(handles,'overlay xform',n);
    dispRange  = dtiGet(handles,'overlay display range',n);

    % Get the color map and add a black where we place NaNs.
    cmap = handles.cmaps(get(handles.popupOverlayCmap,'Value')).rgb;
    cmap = [cmap; 0 0 0];

    interpType = dtiGet(handles,'interp type');

    if(ndims(overlayImg)==3)
        % Get the transformed overlay slices and combine them with the
        % background slices.
        % The disp range we send in is used to return an image between 0
        % and 1.  We need to remember the display range, somehow. The
        % spatial interpolation type is set somewhere.  Not sure where. It
        % is a string, by default n, which I think means nearest neighbor.
        ox = dtiGetSlice(oXform, overlayImg, 1, curPosition(1), bb, interpType, mmPerVox, dispRange);
        oy = dtiGetSlice(oXform, overlayImg, 2, curPosition(2), bb, interpType, mmPerVox, dispRange);
        oz = dtiGetSlice(oXform, overlayImg, 3, curPosition(3), bb, interpType, mmPerVox, dispRange);

        % Find the locations above the threshold
        mx = ox > overlayThresh;
        my = oy > overlayThresh;
        mz = oz > overlayThresh;

        % Here we merge the overlay image (which is colorful) with the
        % background image created above, which is gray scale.
        % We do this by taking every value in the overlay image and
        % assigning it an index into the color map, cmap.
        %
        % The cmap rows are 1:256.  We require the overlay to be between 0
        % and 1. That way the entires of the overlay, e.g., oz, map into a
        % row of cmap.  So, oz*255 + 1 is a an integer between 1 and 256.

        % We set the NaNs to the 257th entry of the table, which is 0 0 0.
        oz(isnan(oz)) = (256 / 255);
        oy(isnan(oy)) = (256 / 255);
        ox(isnan(ox)) = (256 / 255);

        % This makes a list of cmap values that is a vector, and we reshape
        % it to an RGB image.
        oz = reshape(cmap(round(oz*255+1),:), [size(mz) 3]);
        oy = reshape(cmap(round(oy*255+1),:), [size(my) 3]);
        ox = reshape(cmap(round(ox*255+1),:), [size(mx) 3]);

        % Now, we mix the background gray image with the overlay image at
        % the locations specified by mx.
        mx = repmat(mx,[1,1,3]);
        my = repmat(my,[1,1,3]);
        mz = repmat(mz,[1,1,3]);
        xRgb(mx) = (1-overlayAlpha).*xRgb(mx) + overlayAlpha.*ox(mx);
        yRgb(my) = (1-overlayAlpha).*yRgb(my) + overlayAlpha.*oy(my);
        zRgb(mz) = (1-overlayAlpha).*zRgb(mz) + overlayAlpha.*oz(mz);

    elseif(ndims(overlayImg)==4)
        % If there is a 4th dim, then this overlay is an RGB map.
        % dtiGetSlice will return an XxYx3 RGB slice.
        % BROKEN! _ See Franco P
        ox = dtiGetSlice(oXform, overlayImg, 1, curPosition(1), bb, interpType, mmPerVox, dispRange);
        oy = dtiGetSlice(oXform, overlayImg, 2, curPosition(2), bb, interpType, mmPerVox, dispRange);
        oz = dtiGetSlice(oXform, overlayImg, 3, curPosition(3), bb, interpType, mmPerVox, dispRange);
        % We'll treat these differnt from the intensity image overlay case.
        % Rather than use a fixed transparency everywhere, we'll set the
        % transparency per-voxel by modulating it with the luminaince
        % component.
        if(overlayAlpha>1)
            % alpha>1 applies a gamma to the transparency map
            g = 1/overlayAlpha;
            overlayAlpha = 1;
        else
            g = 1;
        end
        lumz = repmat(mean(oz,3).^g,[1,1,3]);
        lumy = repmat(mean(oy,3).^g,[1,1,3]);
        lumx = repmat(mean(ox,3).^g,[1,1,3]);
        mz = lumz>overlayThresh;
        my = lumy>overlayThresh;
        mx = lumx>overlayThresh;
        xRgb(mx) = (1-lumx(mx).*overlayAlpha).*xRgb(mx) + lumx(mx).*overlayAlpha.*ox(mx);
        yRgb(my) = (1-lumy(my).*overlayAlpha).*yRgb(my) + lumy(my).*overlayAlpha.*oy(my);
        zRgb(mz) = (1-lumz(mz).*overlayAlpha).*zRgb(mz) + lumz(mz).*overlayAlpha.*oz(mz);
    end
end

%% Set up the return transform from acpc

% The xform that we return is from acpc to the space that we put the slices
% in, which is ac-pc aligned on the specific sample grid (mmPerVox,
% typically 1mm isotropic) and filling the bounding box. The xform is
% returned and used in say, dtiRefreshFigure.
xform = [[diag(mmPerVox) bb(1,:)'-1];[0 0 0 1]];

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function info = cnn_imagenet_evaluate(varargin)
% CNN_IMAGENET_EVALUATE   Evauate MatConvNet models on ImageNet

run(fullfile(fileparts(mfilename('fullpath')), ...
  '..', '..', 'matlab', 'vl_setupnn.m')) ;

opts.dataDir = fullfile('data', 'ILSVRC2012') ;
opts.expDir = fullfile('data', 'imagenet12-eval-vgg-f') ;
opts.modelPath = fullfile('data', 'models', 'imagenet-vgg-f.mat') ;
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.imdbPath = fullfile(opts.expDir, 'imdb.mat');
opts.networkType = [] ;
opts.lite = false ;
opts.numFetchThreads = 12 ;
opts.train.batchSize = 128 ;
opts.train.numEpochs = 1 ;
opts.train.gpus = [] ;
opts.train.prefetch = true ;
opts.train.expDir = opts.expDir ;

opts = vl_argparse(opts, varargin) ;
display(opts);

% -------------------------------------------------------------------------
%                                                   Database initialization
% -------------------------------------------------------------------------

if exist(opts.imdbPath)
  imdb = load(opts.imdbPath) ;
else
  imdb = cnn_imagenet_setup_data('dataDir', opts.dataDir, 'lite', opts.lite) ;
  mkdir(opts.expDir) ;
  save(opts.imdbPath, '-struct', 'imdb') ;
end

% -------------------------------------------------------------------------
%                                                    Network initialization
% -------------------------------------------------------------------------

net = load(opts.modelPath) ;
if isfield(net, 'net') ;
  net = net.net ;
end
% Cannot use isa('dagnn.DagNN') because it is not an object yet
isDag = isfield(net, 'params') ;

if isDag
  opts.networkType = 'dagnn' ;
  net = dagnn.DagNN.loadobj(net) ;
  trainfn = @cnn_train_dag ;

  % Drop existing loss layers
  drop = arrayfun(@(x) isa(x.block,'dagnn.Loss'), net.layers) ;
  for n = {net.layers(drop).name}
    net.removeLayer(n) ;
  end

  % Extract raw predictions from softmax
  sftmx = arrayfun(@(x) isa(x.block,'dagnn.SoftMax'), net.layers) ;
  predVar = 'prediction' ;
  for n = {net.layers(sftmx).name}
    % check if output
    l = net.getLayerIndex(n) ;
    v = net.getVarIndex(net.layers(l).outputs{1}) ;
    if net.vars(v).fanout == 0
      % remove this layer and update prediction variable
      predVar = net.layers(l).inputs{1} ;
      net.removeLayer(n) ;
    end
  end

  % Add custom objective and loss layers on top of raw predictions
  net.addLayer('objective', dagnn.Loss('loss', 'softmaxlog'), ...
               {predVar,'label'}, 'objective') ;
  net.addLayer('top1err', dagnn.Loss('loss', 'classerror'), ...
               {predVar,'label'}, 'top1err') ;
  net.addLayer('top5err', dagnn.Loss('loss', 'topkerror', ...
                                     'opts', {'topK',5}), ...
               {predVar,'label'}, 'top5err') ;

  % Make sure that the input is called 'input'
  v = net.getVarIndex('data') ;
  if ~isnan(v)
    net.renameVar('data', 'input') ;
  end

  % Swtich to test mode
  net.mode = 'test' ;
else
  opts.networkType = 'simplenn' ;
  net = vl_simplenn_tidy(net) ;
  trainfn = @cnn_train ;
  net.layers{end}.type = 'softmaxloss' ; % softmax -> softmaxloss
end

% Synchronize label indexes used in IMDB with the ones used in NET
imdb = cnn_imagenet_sync_labels(imdb, net);

% Run evaluation
[net, info] = trainfn(net, imdb, getBatchFn(opts, net.meta), ...
                      opts.train, ...
                      'train', NaN, ...
                      'val', find(imdb.images.set==2)) ;

% -------------------------------------------------------------------------
function fn = getBatchFn(opts, meta)
% -------------------------------------------------------------------------
useGpu = numel(opts.train.gpus) > 0 ;
bopts = meta.normalization ;
bopts.numThreads = opts.numFetchThreads ;

% Most networks are trained by resizing images to 256 pixels and then
% cropping a slightly smaller subarea. Reproduce this effect (center
% crop) for a more accurate evaluation (it also avoids resizing the
% images if these have been pre-processed to be of this size, which
% accelerates everything).

bopts.border = 256 - meta.normalization.imageSize ;

switch lower(opts.networkType)
  case 'simplenn'
    fn = @(x,y) getSimpleNNBatch(bopts,x,y) ;
  case 'dagnn'
    fn = @(x,y) getDagNNBatch(bopts,useGpu,x,y) ;
end

% -------------------------------------------------------------------------
function [im,labels] = getSimpleNNBatch(opts, imdb, batch)
% -------------------------------------------------------------------------
images = strcat([imdb.imageDir filesep], imdb.images.name(batch)) ;
im = cnn_imagenet_get_batch(images, opts, ...
                            'prefetch', nargout == 0) ;
labels = imdb.images.label(batch) ;

% -------------------------------------------------------------------------
function inputs = getDagNNBatch(opts, useGpu, imdb, batch)
% -------------------------------------------------------------------------
images = strcat([imdb.imageDir filesep], imdb.images.name(batch)) ;
im = cnn_imagenet_get_batch(images, opts, ...
                            'prefetch', nargout == 0) ;
if nargout > 0
  if useGpu
    im = gpuArray(im) ;
  end
  inputs = {'input', im, 'label', imdb.images.label(batch)} ;
end

author: jiangqy


text:
function x = subsasgn( x, S, y )

%   Disciplined convex/geometric programming information for SUBSASGN:
%      Subscripting can be used to change values of elements or slices
%      of CVX variables in the same manner as with numeric arrays. All
%      conventions are preserved, including the colon ':' and 'end'
%      notation, as well as the ability to "expand" a CVX variable by
%      assigning a value to a location outside of its current dimensions
%      (e.g., X(end+1)=0).
%
%      One notable exception is this: if the right-hand side is a CVX
%      expression, then the left-hand side must be as well. So, for
%      example, the following will fail:
%         variable x;
%         y = ones(3,1);
%         y(2) = x;
%      This is because MATLAB does not know how to automatically promote
%      'y' to a CVX variable so that it can accept 'x' as an element. If
%     you want
%      to accomplish something like this, you must manually convert y into a
%      CVX variable first, as follows:
%         variable x;
%         y = cvx(ones(3,1));
%         y(2) = x;

%
% Test subscripts
%

szx = size( x );
szy = size( y );
nlx = prod( szx );
try
    temp = reshape( 1 : nlx, szx );
    ndx_x = builtin( 'subsasgn', temp, S, zeros( szy ) );
catch errmsg
    throw( errmsg );
end
szx_n = size( ndx_x );
if length( szx_n ) ~= length( szx ),
   szx_n(end+1:length(szx)) = 1;
   szx(end+1:length(szx_n)) = 1;
end

%
% Assign data
%

bx = cvx_basis( x );
if any( szx_n < szx ),
    bx = bx( :, ndx_x );
else
    if any( szx_n > szx ),
        bx( :, end + 1 ) = 0;
        ndx_x( ndx_x == 0 ) = size( bx, 2 );
        bx = bx( :, ndx_x );
        temp = reshape( 1 : prod( szx_n ), szx_n );
    end
    ndx_x = builtin( 'subsref', temp, S );
    ndx_x = ndx_x( : );
    nlz = length( ndx_x );
    by = cvx_basis( y );
    nx = size( bx, 1 );
    [ ny, my ] = size( by );
    if nx < ny,
        if issparse( by ) && ~issparse( bx ), bx = sparse( bx ); end
        bx( ny, : ) = 0;
    elseif nx > ny,
        if issparse( bx ) && ~issparse( by ), by = sparse( by ); end
        by( nx, : ) = 0;
    end
    if my < nlz,
        by = by( :, ones( 1, nlz ) );
    end
    bx( :, ndx_x ) = by;
end

%
% Create the new object
%

x = cvx( szx_n, bx );

% 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:
%   Generic rank-r Newton's iteration routine solving    F(z) = 0  where
%   F is a smooth mapping between vector spaces, assuming the Jacobian
%   is of rank r at the solution that can be non-isolated
%   (c.f. http://homepages.neiu.edu/~zzeng/Papers/Rank-r_Newton.pdf)
%
%     Syntax
%         [z,res,fcond] = Newton({@F,domain,para}, @G, z0)
%         [z,res,fcond] = Newton({@F,domain,para}, @G, z0, trac)
%         [z,res,fcond] = Newton({@F,domain,para}, @G, z0, trac, tol)
%         [z,res,fcond] = Newton({@F,domain,para}, {@G,r}, z0)
%         [z,res,fcond] = Newton({@F,domain,para}, {@G,r}, z0, trac)
%         [z,res,fcond] = Newton({@F,domain,para}, {@G,r}, z0, trac, tol)
%
%     Input:   F -- (function) Matlab function name for calculating y = F(z)
%                    with syntax >> y = F(z1,...zk,p1,...,pm) where z1,...,zk
%                    are components of z and p1,...,pm are fixed parameters.
%                    (F must run on F(domain{:},para{:}).
%         domain -- domain of the homomorphism represented by
%                        (i) an mxn matrix of 0 'and 1's representing
%                            0 entries and variable entries respectively, or
%                       (ii) a polynomial with variable coefficients being
%                              nonzero, or
%                      (iii) a cell array of matrices in (i) and/or (ii)
%                      assuming the domain is a product of matrix spaces and
%                      polynomial spaces
%            para -- (cell) parameters needed for running F
%               G -- (function) Matlab function name for calculating the
%                     Jacobian u = F_z(z0)y F as a homomorphism with syntax
%                        >> u = G(y1,...,yk,z1,...zk,p1,...,pm)
%                     where z1,...,zk,p1,...,pk are the same for F but
%                     considered fixed, the variables are y1,...,yk in the
%                     same domain of F.
%               r -- (integer, optional) if provided, the rank-r projection
%                       of the Jacobian is used for the iteration
%              z0 -- (cell/matrix/polynomial) the initial iterate for the
%                    Newton's iteration in the domain
%            trac -- (0,1,or 2, optional) flag for tracking the iteration
%                     if tracking = 0 or missing, no tracking
%                     if tracking = 1, track the first component z1 of z
%                     if tracking = 2, track all components of z
%              tol -- (numeric, optional) error tolerance of the iterate
%                     if missing, iteration stops when residule stops
%                     decreasing
%
%   Output:         z  -- the least squares solution in the domain
%                 res  -- (optional) the residual ||F(z)||
%               fcond  -- (optional) the condition number

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% StackExchange Signal Processing Q79314
% https://dsp.stackexchange.com/questions/79314
% Image Segmentation Using Deep Learning
% References:
%   1.
% Remarks:
%   1.  B
% TODO:
% 	1.  C
% Release Notes
% - 1.0.000     27/11/2021
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

imageSetFolderName  = 'ImageSet';
maskSetFolderName   = 'MaskSet';

resizePostName = 'Resize';

numRows     = 64;
numCols     = 64;

%% Generate Data

sFiles = dir(strcat(imageSetFolderName, resizePostName, FILE_SEP, '*.jpg'));
for ii = 1:length(sFiles)
    mI = imread(strcat(sFiles(ii).folder, FILE_SEP, sFiles(ii).name));
    assert(size(mI, 1) == numRows);
    assert(size(mI, 2) == numCols);
    assert(size(mI, 3) == 3);
end

sFiles = dir(strcat(maskSetFolderName, resizePostName, FILE_SEP, '*.png'));
for ii = 1:length(sFiles)
    mI = imread(strcat(sFiles(ii).folder, FILE_SEP, sFiles(ii).name));
    assert(size(mI, 1) == numRows);
    assert(size(mI, 2) == numCols);
    assert(size(mI, 3) == 1);
end

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function [X,Y,Z,CAPS] = extrude(varargin)
% EXTRUDE  Extrude a 2D curve along a 3D path to create
% tubes/cylinders/etc (and, optionally, fly though it!)
%
%   [X,Y,Z] = EXTRUDE(base,traj)
%   [X,Y,Z,CAPS] = EXTRUDE(base,traj,cap)
%   [X,Y,Z,CAPS] = EXTRUDE(base,traj,cap,alg)
%   [X,Y,Z,CAPS] = EXTRUDE(base,traj,cap,alg,fly) returns the surface and
%   handle to the surface plot object created by the 2D base curve and 3D
%   trajectory. SURF(X,Y,Z) can be used to display the object. If cap == 2,
%   data for the caps will be sent separately.
%
%   base = [x;y]   A 2-row Matrix defining the 2d base
%   traj = [x;y;z] A 3-row Matrix defining the 3d path
%   cap = Close off the ends of the tube:
%            0 -> No (default)
%            1 -> Yes
%            2 -> Yes but output them separately into the variable CAPS
%
%   alg = Creation algorithm:
%            1 Creates a natural looking surface with little twist (default)
%            2 More twisty, but preserves orientation vs. direction (which
%              is better for periodicity)
%   fly =  N  Fly through the path generated N times
%
%
%   DEMO: EXTRUDE with no arguments flies through a random periodic path
%         with a star/road/tube cross section using either algorithm
%
%   EXAMPLE 1: Generate a torus:
%
%       q = linspace(0,2*pi,33);
%       base = [cos(q); sin(q)];   % Base curve is a circle (radius = 1)
%       q = linspace(0,2*pi,101);
%       traj = 5*[cos(q); sin(q); 0*q]; %Trajectory is a circle (radius = 5)
%
%       [X,Y,Z] = extrude(base,traj);
%       figure; surf(X,Y,Z); axis equal;
%
%   EXAMPLE 2: Generate half a torus, with nice caps plotted separately:
%
%       q = linspace(0,2*pi,33);
%       base = [cos(q); sin(q)];   % Base curve is a circle (radius = 1)
%       q = linspace(0,pi,101);
%       traj = 5*[cos(q); sin(q); 0*q]; %Trajectory is a semicircle (radius = 5)
%
%       [X,Y,Z,CAPS] = extrude(base,traj,2);    %cap = 2 for separate caps
%       figure; surf(X,Y,Z); axis equal; hold on;
%       surf(CAPS(1).X,CAPS(1).Y,CAPS(1).Z,'linestyle','none');
%       surf(CAPS(2).X,CAPS(2).Y,CAPS(2).Z,'linestyle','none');
%
%
%   EXAMPLE 3: Generate a 3D square tube and fly through it once:
%
%       q = linspace(0,2*pi,5) + pi/4;
%       base = 0.15*[cos(q); sin(q)];   % Base curve is a square
%       t = linspace(0,1,2001);
%       traj = [sin(4*pi*t+1); sin(5*pi*t+2); sin(6*pi*t)];
%
%       [X,Y,Z] = extrude(base,traj,0,[],1);
%       figure; surf(X,Y,Z); axis equal;
%
%   2009-12-25 v3.1     Fixed an error with argument parsing
%   2009-12-18 v3.0     Added an option to close off the ends, and changed
%                       the order of the input arguments.
%   2009-08-20 v2.0     Added a better algorithm that allow for much less
%                       twistiness, a new (simpler) example, and better
%                       error checking.
%   2009-08-20 v1.1     Fixed an error with passing in only 2 arguments and
%                       clearer documentation
%   2009-08-19 v1.0     Initial Creation

%% Process input parameters
if nargin < 2  % Default base = star,   Default path = random
    disp('EXTRUDE DEMO: Fly through a random periodic path.')
    sw =  input('Pick a cross-section: (1)Star (2)Road (3)Tube ? ');
    switch sw
        case 1
            bs = 5; sharp = 0.5; bsize = .1;   %5 pointed star...
            basex = [cos(linspace(0,2*pi,bs+1));
                    sharp*cos(2*pi/2/bs + linspace(0,2*pi,bs+1))];
            basey = [sin(linspace(0,2*pi,bs+1));
                    sharp*sin(2*pi/2/bs + linspace(0,2*pi,bs+1))];
            base = bsize*[basex(1:2*bs+1); basey(1:2*bs+1)];
        case 2
            base = [-1 -1 1 1; 0 -.2 -.2 0]/10;
        case 3
            base = .1*[cos(linspace(0,2*pi,12)); sin(linspace(0,2*pi,12))];
        otherwise
            disp('Pick 1, 2, or 3! Try again.'); return;
    end

    alg = input('Pick an algorithm: (1) Normal (2) Twisty but orientation preserving ? ');
    if alg ~= 1 && alg ~= 2
        disp('Pick 1 or 2! Try again.'); return;
    end

    C = interpft(randn(10,3),1001)'; %Generate a smooth periodic path
    C = [C(:,end) C]; %close the loop
    C = 1*C./max(abs(C(:))); %normalize it
	fly = Inf;   %Fly through it forever
    cap = 0;    % No caps

elseif nargin >= 2;
    fly = []; alg = 1; cap = 0;
    C = varargin{2};
    base = varargin{1};
    if nargin >= 3 && ~isempty(varargin{3}); cap = varargin{3}; end
    if nargin >= 4 && ~isempty(varargin{4}); alg = varargin{4}; end
    if nargin == 5 && ~isempty(varargin{5}); fly = varargin{5}; end
end

% Check size to make size(C) = 3xN
if size(C,2) == 3 && size(C,1) ~= 3
    C = C';
end

% Check size to make size(base) = 2xN
if size(base,2) == 2 && size(base,1) ~= 2
    base = base';
end

%% Calculate derivatives and allocate matrices
npt = size(base,2);
base = [base; zeros(1,npt)];

if size(C,2) >= 3 %Use a 2nd order approximation for the derivatives of the trajectory
    dC = [C(:,1:3)*[-3; 4; -1]/2  [C(:,3:end) - C(:,1:end-2)]/2 C(:,end-2:end)*[1; -4; 3]/2];
else
    dC = C(:,[2 2]) - C(:,[1 1]);
end

dC0 = find(sum(abs(dC),1) == 0,1);    %Check for stagnation points
if ~isempty(dC0)
    warning('Removing stagnation points found in trajectory');
    dCgood = find(sum(abs(dC),1) ~= 0);
    C = C(:,dCgood); %
    dC = dC(:,dCgood); %
end

K = size(dC,2);
SUR = nan(3,npt,K);
camtar = zeros(3,K);
camup = zeros(3,K);
camdata = [[0;0;1],[0;1;0]];
%% Generate and plot the surface

dCvec_prev = [0;0;1];

switch alg
    case 1
        %% Natural Looking Tube (default)
        for k = 1:K
            % We want to rotate [0;0;1] 180degrees around an axis 'z' to become dC

            dCvec = dC(:,k)/norm(dC(:,k));
            z = cross(dCvec_prev,dCvec);

            if norm(z) ~= 0
                z = z/norm(z);
                q = real(acos(dot(dCvec_prev,dCvec)/norm(dCvec_prev)/norm(dCvec)));

                Z = repmat(z,1,npt);
                base = base*cos(q) + cross(Z,base)*sin(q)+Z*(1-cos(q))*diag(dot(Z,base));
                camdata = camdata*cos(q) + cross([z z],camdata)*sin(q)+[z z]*(1-cos(q))*diag(dot([z z],camdata));

                dCvec_prev = dCvec;
            end

            % Data for the camera used in the flying routine
            camtar(:,k) = camdata(:,1);
            camup(:,k) = camdata(:,2);

            SUR(:,:,k) = base + repmat(C(:,k),1,npt);

        end
    case 2
        %% Orientation preserving tube (more twisty but allows for periodicity)
        for k = 1:K
            % We want to rotate [0;0;1] 180degrees around an axis 'z' to become dC
            dCvec = dC(:,k)/norm(dC(:,k));
            if isequal(dCvec,[0;0;-1]) %Prevents a 0/0 = nan
                z = [0;1;0];
            else
                z = ([0; 0; 1] + dCvec)/2; z = z/norm(z);
            end
            z = repmat(z,1,npt);

            % Data for the camera used in the flying routine
            camtar(:,k) = z(:,1)*z(3)*2 - [0;0;1];
            camup(:,k) = z(:,1)*z(2)*2 - [0;1;0];

            SUR(:,:,k) = z*diag(dot(z,base))*2 - base + repmat(C(:,k),1,npt);

        end
    otherwise, error('Algorthm 1 or 2?');
end

CAPS = [];
if cap == 1; % Add caps
    SUR = cat(3,repmat(C(:,1),1,npt),SUR,repmat(C(:,K),1,npt));
end

X = squeeze(SUR(1,:,:));
Y = squeeze(SUR(2,:,:));
Z = squeeze(SUR(3,:,:));

if cap == 2 % Add separate caps
    SURCAP1 = cat(3,SUR(:,:,1),repmat(C(:,1),1,npt));
    SURCAP2 = cat(3,SUR(:,:,K),repmat(C(:,K),1,npt));
    CAPS(1).X = squeeze(SURCAP1(1,:,:));
    CAPS(1).Y = squeeze(SURCAP1(2,:,:));
    CAPS(1).Z = squeeze(SURCAP1(3,:,:));
    CAPS(2).X = squeeze(SURCAP2(1,:,:));
    CAPS(2).Y = squeeze(SURCAP2(2,:,:));
    CAPS(2).Z = squeeze(SURCAP2(3,:,:));
end

if isempty(fly) || fly <= 0
    return
end

if exist('sw','var') %If in DEMO mode, make two figures
    surf(X,Y,Z);
    axis equal;
end

curfig = figure;
h = surf(X,Y,Z);
axis equal;

%% Fix any NaN's (This should not happen...)
xisnan = find(isnan(X(1,:)), 1);
xnotnan = find(~isnan(X(1,:)));
if ~isempty(xisnan)
    warning('NaN''s found');
    X = X(:,xnotnan);
    Y = Y(:,xnotnan);
    Z = Z(:,xnotnan);
    camup = camup(:,xnotnan);
    camtar = camtar(:,xnotnan);
    C = C(:,xnotnan);
end

%% Fly through it
set(h,'facealpha',0.5,'edgecolor','flat'); colormap(hsv); %Make it look nice
set(gca,'visible','off','projection','perspective');
set(gcf,'color','k');
set(gcf,'closerequestfcn','assignin(''caller'',''stopme'',1); delete(gcbf)'); % Stop the loop when closing the figure
h = gcf;
K = size(C,2);
k = 0;
laps = 1;
stopme = 0;

while laps <= ceil(fly) && gcf == curfig && stopme == 0
    k = 1 + mod(k,K);  %Increment k but loop from K back to 1

    set(gca,'cameraviewangle',150);
    set(gcf,'name',['Lap ' num2str(laps) '/' num2str(fly) '  ' num2str(k) '/' num2str(K) '   Close the figure to stop.']);
    set(gca,'cameraposition',C(:,k),'cameratarget',C(:,k)+camtar(:,k),'cameraupvector',camup(:,k));
    drawnow;

    if k == K; laps=laps+1; end;

end
try close(h); catch; end; % If the figure remains after leaving the while loop


text:
function msg = mhe_init()

s = 'http://byu.apmonitor.com';
b = 'mhe';

apm(s,b,'clear all');

% load model and data
apm_load(s,b,'model.apm');
csv_load(s,b,'mhe.csv');

% configure MV / CV
apm_info(s,b,'FV','Kp');
apm_info(s,b,'FV','tau');
apm_info(s,b,'FV','zeta');
apm_info(s,b,'FV','TC_ss');
apm_info(s,b,'MV','Q1');
apm_info(s,b,'CV','TC');

% dynamic estimation
apm_option(s,b,'nlc.imode',5);
apm_option(s,b,'nlc.solver',3);

% tune FV
apm_option(s,b,'Kp.dmax',0.05);
apm_option(s,b,'Kp.lower',0.3);
apm_option(s,b,'Kp.upper',0.6);
apm_option(s,b,'tau.dmax',2.0);
apm_option(s,b,'tau.lower',40.0);
apm_option(s,b,'tau.upper',60.0);
apm_option(s,b,'zeta.dmax',0.05);
apm_option(s,b,'zeta.lower',1.3);
apm_option(s,b,'zeta.upper',1.7);
apm_option(s,b,'TC_ss.dmax',1.0);
apm_option(s,b,'TC_ss.lower',22.0);
apm_option(s,b,'TC_ss.upper',25.0);
% turn on FVs as degrees of freedom
apm_option(s,b,'Kp.status',1);
apm_option(s,b,'tau.status',1);
apm_option(s,b,'zeta.status',1);
apm_option(s,b,'TC_ss.status',1);

apm_option(s,b,'Kp.fstatus',0);
apm_option(s,b,'tau.fstatus',0);
apm_option(s,b,'zeta.fstatus',0);
apm_option(s,b,'TC_ss.fstatus',0);

% read Q, don't let optimize use MV
apm_option(s,b,'Q1.status',0);
apm_option(s,b,'Q1.fstatus',1);
% include CV in objective function
apm_option(s,b,'TC.status',1);
apm_option(s,b,'TC.fstatus',1);

% web-viewer option, update every second
apm_option(s,b,'nlc.web_plot_freq',2);

msg = 'initialization complete';

end

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
% set default value for a field if not set yet.
function para = SetDefaultValue(para, fieldName, defaultValue)
    if ~isfield(para, fieldName)
        para.(fieldName) = defaultValue;
    end
end

author: singaxiong

repo: SignalGraph


text:
%gaussian mixture model
%>
%> @param FeatureMatrix: features for all train observations (dimension iNumFeatures x iNumObservations)
%> @param k: number of gaussians
%> @param numMaxIter: maximum number of iterations (stop if not converged before)
%> @param prevState: internal state that can be stored to continue clustering later
%>
%> @retval mu means
%> @retval sigma standard deviations
%> @retval state result containing internal state (if needed)
% ======================================================================
function [mu, sigma, state] = ToolGmm(V, k, numMaxIter, prevState)

    if (nargin < 3)
        numMaxIter  = 1000;
    end
    if (nargin == 4)
        state = prevState;
    else
        % initialize state
        state = initState_I(V, k);
    end

    for j = 1:numMaxIter
        prevState = state;

        % compute weighted gaussian
        p = computeProb_I(V, state);

        % update clusters
        state = updateGaussians_I(V, p, state);

        % if we have converged, break
        if (max(sum(abs(state.m-prevState.m))) <= 1e-20)
            break;
        end
    end

    mu = state.m;
    sigma = state.sigma;
end

function [state] = updateGaussians_I(FeatureMatrix, p, state)

    % number of clusters
    K = size(state.m, 2);

    % update priors
    state.prior = mean(p, 1)';

    for k = 1:K
        s = 0;

        % update means
        state.m(:, k) = FeatureMatrix * p(:, k) / sum(p(:, k));

        % subtract mean
        F = FeatureMatrix - repmat(state.m(:, k), 1, size(FeatureMatrix, 2));

        for n = 1:size(FeatureMatrix, 2)
            s = s + p(n, k) * (F(:, n) * F(:, n)');
        end
        state.sigma(:, :, k) = s / sum(p(:, k));
    end
end

function [p] = computeProb_I(FeatureMatrix, state)

    K = size(state.m, 2);
    p = zeros(size(FeatureMatrix, 2), K);

    % for each cluster
    for k = 1:K
        % subtract mean
        F = FeatureMatrix - repmat(state.m(:, k), 1, size(FeatureMatrix, 2));

        % weighted gaussian
        p(:, k) = 1 / sqrt((2*pi)^size(F, 1) * det(state.sigma(:, :, k))) *...
            exp(-1/2 * sum((F .* (inv(state.sigma(:, :, k)) * F)), 1)');
        p(:, k) = state.prior(k) * p(:, k);
    end

    % norm over clusters
    p = p ./ repmat(sum(p, 2), 1, K);
end

function [state] = initState_I(FeatureMatrix, K)

    %init
    m       = zeros(size(FeatureMatrix, 1), K);
    sigma   = zeros(size(FeatureMatrix, 1), size(FeatureMatrix, 1), K);
    prior   = zeros(1, K);

    % pick random points as cluster means
    mIdx    = round(rand(1, K) * (size(FeatureMatrix, 2)-1)) + 1;

    % assign means etc.
    m       = FeatureMatrix(:, mIdx);
    prior   = ones(1, K) / K;
    sigma   = repmat(cov(FeatureMatrix'), 1, 1, K);

    % write initial state
    state   = struct('m', m, 'sigma', sigma, 'prior', prior);
end

author: alexanderlerch

repo: ACA-Code


text:
% The COBRAToolbox: testSWIFTCORE.m
%
% Purpose:
%     - tests the basic functionality of swiftcore algorithm
%
% Authors:
%     - Original file: Mojtaba Tefagh, March 2019
%

global CBTDIR

% require the specified toolboxes and solvers
solvers = prepareTest('needsLP', true, 'requireOneSolverOf', {'gurobi','ibm_cplex'},'excludeSolvers', {'matlab', 'lp_solve','pdco'},'requiredToolboxes', {'statistics_toolbox'});

% save the current path
currentDir = pwd;

% initialize the test
testPass = fileparts(which('testSWIFTCORE.m'));
cd(testPass);

% set the cobra solver
solverLPOK = changeCobraSolver(solvers.LP{1}, 'LP', 0);

% load the model
model = getDistributedModel('ecoli_core_model.mat');
model.rev = double(model.lb < 0);
A = swiftcc(model.S, model.rev,solvers.LP{1});
model.S = model.S(:, A);
model.rev = model.rev(A);
model.ub = model.ub(A);
model.lb = model.lb(A);
model.rxns = model.rxns(A);
n = length(model.rev);
core = randsample(n, round(n/2));

fprintf(' -- Running swiftcore w/o reduction and using the %s solver...\n', solvers.LP{1});
[~, coreInd, ~] = swiftcore(model, core, ones(n, 1), 1e-10, false, solvers.LP{1});
assert(all(coreInd(core)));
A = swiftcc(model.S(:, coreInd), model.rev(coreInd), solvers.LP{1});
%A = swiftcc(model.S(:, coreInd), model.rev(coreInd));
assert(all(A.' == 1:length(A)));
fprintf(' -- Running swiftcore w/ reduction and using the %s solver...\n', solvers.LP{1});
[~, coreInd, ~] = swiftcore(model, core, ones(n, 1), 1e-10, true, solvers.LP{1});
assert(all(coreInd(core)));
A = swiftcc(model.S(:, coreInd), model.rev(coreInd),solvers.LP{1});
tmp=nnz(A.' == 1:length(A))/length(A)
bool=all(A.' == 1:length(A));
assert(bool);

% output a success message
fprintf('\nDone.\n');

% load the model
model = getDistributedModel('ecoli_core_model.mat');
model.rev = double(model.lb < 0);
A = fastcc(model, 1e-4, 0);

fprintf('\n -- Running swiftcc using the %s solver...\n', solvers.LP{1});
consistent = swiftcc(model.S, model.rev, solvers.LP{1});
assert(all(A == consistent));

[solverName, solverOK] = getCobraSolver('LP');
fprintf('\n -- Running swiftcc using the default LP solver, which is %s,....\n', solverName);
consistent = swiftcc(model.S, model.rev,solvers.LP{1});
assert(all(A == consistent));

% output a success message
fprintf('\nDone.\n');

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function ens = comp_eff_node_sizes(ns, cnodes, ev, domain)

dnodes = mysetdiff(1:length(ns), cnodes);
odom = domain(~isemptycell(evidence(domain)));
cdom = myintersect(cnodes, domain);
ddom = myintersect(dnodes, domain);
cobs = myintersect(cdom, odom);
dobs = myintersect(ddom, odom);
ens = ns;
ens(cobs) = 0;
ens(dobs) = 1;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function a = repmat(a,m,n)
%REPMAT       Implements  repmat(a)  for Hessians
%
%Functionality as in Matlab.
%

% written  05/03/09     S.M. Rump
%

%Code adapted from Matlab R2008a
%

if nargin < 2
    error('repmat requires at least 2 inputs')
end

if nargin == 2
    if prod(size(m))==1
        siz = [m m];
    else
        siz = m;
    end
else
    siz = [m n];
end

if prod(size(a))==1
    %VVVV  a = a(ones(siz));
    s.type = '()'; s.subs = {ones(siz)}; a = subsref(a,s);
    %AAAA  Matlab bug fix
elseif ndims(a) == 2 & length(siz) == 2
    [m,n] = size(a);

    if (m == 1 & siz(2) == 1)
        %VVVV  a = a(ones(siz(1), 1), :);
        s.type = '()'; s.subs = {ones(siz(1), 1), ':'}; a = subsref(a,s);
        %AAAA  Matlab bug fix
    elseif (n == 1 & siz(1) == 1)
        %VVVV  a = a(:, ones(siz(2), 1));
        s.type = '()'; s.subs = {':', ones(siz(2), 1)}; a = subsref(a,s);
        %AAAA  Matlab bug fix
    else
        mind = (1:m)';
        nind = (1:n)';
        mind = mind(:,ones(1,siz(1)));
        nind = nind(:,ones(1,siz(2)));
        %VVVV  a = a(mind,nind);
        s.type = '()'; s.subs = {mind,nind}; a = subsref(a,s);
        %AAAA  Matlab bug fix
    end
else
    asiz = size(a);
    asiz = [asiz ones(1,length(siz)-length(asiz))];
    siz = [siz ones(1,length(asiz)-length(siz))];
    for i=length(asiz):-1:1
        ind = (1:asiz(i))';
        subs{i} = ind(:,ones(1,siz(i)));
    end
    %VVVV  a = a(subs{:});
    s.type = '()'; s.subs = {subs{:}}; a = subsref(a,s);
    %AAAA  Matlab bug fix
end

repo: OpenMAS


text:
function [T,P]=fnipals(X,w,T)

%FNIPALS nipals algorithm for PCA
%
% function [T,P]=fnipals(X,w,T)
%
% 'fnipals.m'
%
% This algorithm requires the presence of:
% 'missmean.m'
%
% ----------------------------------------------------
%        Find eigenvectors according to NIPALS
% ----------------------------------------------------
%
% [T,P]=fnipals(X,w,T);
% [T,P]=fnipals(X,w);
%
% T is found so that X = T*P', s.t ||T||=1 and T'T=I
%
% X        : The matrix to be decomposed.
% w        : Number of factors to extract.
%            If w is high (perhaps>20) consider using SVD.
% T        : Initial guess of the solution, optional.
%            If T is not specified, a little time will
%            be used on finding orthogonal random
%            starting values.
%
% You may want to calculate P afterwards by typing 'P=X*T'.
% Note that the T returned is orthonormal.
% Calculation of P is left of this implementation to save FLOP's.
% It handles missing values NaNs (very dispersed, less than 15%)
% If the problem is small enough you would prefer the SVD rather
% than NIPALS for finding T. NIPALS may be inaccurate when
% extracting too many factors, i.e., many more than the rank
% of X.

%scalar ConvLim WarnLim ItMax a b i

% $ Version 1.01 $ Date 18. June 1998 $ Not compiled $

% Copyright (C) 1995-2006  Rasmus Bro & Claus Andersson
% Copenhagen University, DK-1958 Frederiksberg, Denmark, rb@life.ku.dk
%
% 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., 51 Franklin
% Street, Fifth Floor, Boston, MA  02110-1301, USA.

ConvLim=1e-12;
WarnLim=1e-4;
ConvLimMiss=100*ConvLim;
ItMax=100;

filename='fnipals.m';

[a b]=size(X);

if (w>a | w>b) | w<1,
    help(filename);
    error(['Error in ' filename ': Number of factors to extract is invalid!'])
end;

np=isnan(X);
MissingExist=any(np);

if ~exist('T'),
    T=orth(randn(a,w));
end;

if exist('P'),
    P=[];
end;

if ~MissingExist
    if (size(T) == [a w]),
        if a>b,
            P=X'*T;
            l2=Inf;
            Z=X'*X;
            for i=1:w,
                p=P(:,i);
                d=1;
                it=0;
                while (d>ConvLim) & (it<ItMax),
                    it=it+1;
                    p=Z*p;
                    l1=sqrt(p'*p);
                    p=p/l1;
                    d=(l1-l2)^2;
                    l2=l1;
                end;
                P(:,i)=sqrt(l1)*p;
                Z=Z-P(:,i)*P(:,i)';
                WarnLim=sqrt(l1)/1000;
                if it>=ItMax & d>WarnLim,
                    disp('FNIPALS, High-X: Iterated up to the ItMax limit!')
                    disp('FNIPALS, High-X: The solution has not converged!')
                end;
            end;
            T=X*P;
        else
            P=[];
            l2=Inf;
            Z=X*X';
            for i=1:w,
                t=T(:,i);
                d=1;
                it=0;
                while (d>ConvLim) & (it<ItMax),
                    it=it+1;
                    t=Z*t;
                    l1=sqrt(t'*t);
                    t=t/l1;
                    d=(l1-l2).^2;
                    l2=l1;
                end;
                T(:,i)=sqrt(l1)*t;
                Z=Z-T(:,i)*T(:,i)';
                WarnLim=sqrt(l1)/1000;
                if it>=ItMax & d>WarnLim,
                    disp('FNIPALS, Wide-X: Iterated up to the ItMax limit!')
                    disp('FNIPALS, Wide-X: The solution has not converged!')
                end;
            end;
        end;
        T=gsm(T);
    else
        error(['Error in ' filename ': Number of factors to extract is invalid!'])
    end;
else
    MissIdx=find(np);
    [i j]=find(np);
    mnx=missmean(X)/2;
    mny=missmean(X')/2;
    n=size(i,1);
    for k=1:n,
        i_i=i(k);
        j_j=j(k);
        X(i_i,j_j) = mny(i_i) + mnx(j_j);
    end;
    mnz=(missmean(mnx)+missmean(mny))/2;

    ssmisold=sum(sum( X(MissIdx).^2 ));
    sstotold=sum(sum( X.^2 ));
    ssrealold=sstotold-ssmisold;
    iterate=1;
    while iterate

        if (size(T) == [a w]),
            if a>b,
                P=X'*T;
                l2=Inf;
                Z=X'*X;
                for i=1:w,
                    p=P(:,i);
                    d=1;
                    it=0;
                    while (d>ConvLim) & (it<ItMax),
                        it=it+1;
                        p=Z*p;
                        l1=sqrt(p'*p);
                        p=p/l1;
                        d=(l1-l2)^2;
                        l2=l1;
                    end;
                    P(:,i)=sqrt(l1)*p;
                    Z=Z-P(:,i)*P(:,i)';
                    WarnLim=sqrt(l1)/1000;
                    if it>=ItMax & d>WarnLim,
                        disp('FNIPALS, High-X: Iterated up to the ItMax limit!')
                        disp('FNIPALS, High-X: The solution has not converged!')
                    end;
                end;
                T=X*P;
            else
                P=[];
                l2=Inf;
                Z=X*X';
                for i=1:w,
                    t=T(:,i);
                    d=1;
                    it=0;
                    while (d>ConvLim) & (it<ItMax),
                        it=it+1;
                        t=Z*t;
                        l1=sqrt(t'*t);
                        t=t/l1;
                        d=(l1-l2).^2;
                        l2=l1;
                    end;
                    T(:,i)=sqrt(l1)*t;
                    Z=Z-T(:,i)*T(:,i)';
                    WarnLim=sqrt(l1)/1000;
                    if it>=ItMax & d>WarnLim,
                        disp('FNIPALS, Wide-X: Iterated up to the ItMax limit!')
                        disp('FNIPALS, Wide-X: The solution has not converged!')
                    end;
                end;
            end;
            T=gsm(T);
        else
            error(['Error in ' filename ': Number of factors to extract is invalid!'])
        end;

        P=X'*T;
        Xm=T*P';
        X(MissIdx)=Xm(MissIdx);
        ssmis=sum(sum( Xm(MissIdx).^2 ));
        sstot=sum(sum( X.^2 ));
        ssreal=sstot-ssmis;
        if abs(ssreal-ssrealold)<ConvLim*ssrealold & abs(ssmis-ssmisold)<ConvLimMiss*ssmisold,
            iterate=0;
        end;
        ssrealold=ssreal;
        ssmisold=ssmis;
    end;
end;
T=gsm(T);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function Op=operatornew(otype,varargin)
%OPERATORNEW  Construct a new operator
%   Usage: F=operatornew(otype,...);
%
%   `Op=operatornew(otype,...)` constructs a new operator object *Op* of type
%   *otype*. Arguments following *otype* are specific to the type of operator
%   chosen.
%
%   Frame multipliers
%   -----------------
%
%   `operatornew('framemul',Fa,Fs,s)` constructs a frame multiplier with
%   analysis frame *Fa*, synthesis frame *Fs* and symbol *s*. See the help on
%   |framemul|.
%
%   Spreading operators
%   -------------------
%
%   `operatornew('spread',s)` constructs a spreading operator with symbol
%   *s*. See the help on |spreadop|.
%
%   See also: operator, ioperator

if nargin<2
  error('%s: Too few input parameters.',upper(mfilename));
end;

if ~ischar(otype)
  error(['%s: First agument must be a string denoting the type of ' ...
         'frame.'],upper(mfilename));
end;

otype=lower(otype);

switch(otype)
  case 'framemul'
    Op.Fa=varargin{1};
    Op.Fs=varargin{2};
    Op.s =varargin{3};
    Op.L =framelengthcoef(Op.Fs,size(Op.s,1));
  case 'spread'
    Op.s =varargin{1};
    Op.L =length(Op.s);
  otherwise
    error('%s: Unknows operator type: %s',upper(mfilename),otype);
end;

Op.type=otype;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
classdef DaeHandler < handle

  properties
    ode
    alg

    userdata_p
  end

  methods

    function self = DaeHandler(userdata)
      self.ode = struct;
      self.alg = [];
      self.userdata_p = userdata;
    end

    function r = userdata(self)
      r = self.userdata_p;
    end

     function setODE(self,id,eq)
      if isfield(self.ode, id)
        ocl.utils.exception(['Ode for var ', id, ' already defined']);
      end
      self.ode.(id) = ocl.Variable.getValueAsColumn(eq);
    end

    function setAlgEquation(self,eq)
      self.alg = [self.alg;ocl.Variable.getValueAsColumn(eq)];
    end

    function r = getOde(self, nx, statesOrder)

      r = cell(length(statesOrder),1);
      for k=1:length(statesOrder)
        id = statesOrder{k};
        if ~isfield(self.ode,id)
          ocl.utils.exception(['Ode for state ', id, ' not defined.']);
        end
        r{k} = self.ode.(id);
        self.ode = rmfield(self.ode, id);
      end
      r = vertcat(r{:});

      if length(r) ~= nx
        ocl.utils.exception(['Number of ode equations does not match ',...
                      'number of state variables.']);
      end

      if numel(fieldnames(self.ode)) > 0
        ocl.utils.exception('ODE for variables defined that do not exist.');
      end
    end

    function alg = getAlg(self, nz)
      alg = self.alg;
      if length(alg) ~= nz
        ocl.utils.exception(['Number of algebraic equations does not match ',...
                      'number of algebraic variables.']);
      end
    end
  end
end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function requiredMetsStatus = checkModelFunction(model, requiredMets)
% Checks that the model is able to produce the list of metabolites defined in
% `requiredMets`
%
% USAGE:
%    requiredMetsStatus = checkModelFunction(model, requiredMets)
%
% INPUTS:
%    model:                 model structure
%    requiredMets:          cell array with the list of metabolites that
%                           the model need to be able to produce
%
% OUTPUT:
%    requiredMetsStatus:    0 = test fail,
%                           1 = test succeed
%
% Authors: - This script is an adapted version of the implementation from
%            https://github.com/jaeddy/mcadre.
%          - Modified and commented by A. Richelle,May 2017

    if isempty(requiredMets) %If requiredMets is empty, skip the step and always say the test is passed
        requiredMetsStatus = true;
        warning('No metabolites defined to check the model function')
    else
        % Add demand reactions for required metabolites
        % function evalc used to remove display of reaction added
        [~, model, requiredRxns] = evalc(['addDemandReaction(model,requiredMets);']);

        % Check that the demand reaction can carry fluxes (default value of
        % flux active - 1e-8)
        inactiveRequired = checkRxnFlux(model, requiredRxns);
        requiredMetsStatus = ~numel(inactiveRequired);
    end

end


text:
%%                   runCoordinateDescent.m

% This test file implements the Coordinate Descent solver. The code builds
% a synthetic formulation of the Phase Retrieval problem, b = |Ax| and
% computes an estimate to x. The code finally plots a convergence curve
% and also makes a scatter plot of the true vs recovered solution.

% PAPER TITLE:
%              Coordinate Descent Algorithms for Phase Retrieval

% ARXIV LINK:
%              https://arxiv.org/abs/1706.03474

% 1) Each test script starts out by defining the length of the unknown
% signal, n and the number of measurements, m. These mesurements can be
% made complex by setting the isComplex flag to be true.
%
% 2) We then build the test problem by invoking the function
% 'buildTestProblem' which generates random gaussian measurements according
% to the user's choices in step(1). The function returns the measurement
% matrix 'A', the true signal 'xt' and the measurements 'b0'.
%
% 3) We set the options for the PR solver. For example, the maximum
% number of iterations, the tolerance value, the algorithm and initializer
% of choice. These options are controlled by setting the corresponding
% entries in the 'opts' struct.  Please see the user guide for a complete
% list of options.
%
% 4) We solve the phase retrieval problem by running the following line
% of code:
%   >>  [x, outs, opts] = solvePhaseRetrieval(A, A', b0, n, opts)
% This solves the problem using the algorithm and initialization scheme
% specified by the user in the struct 'opts'.
%
% 5) Determine the optimal phase rotation so that the recovered solution
% matches the true solution as well as possible.
%
% 6) Report the relative reconstruction error. Plot residuals (a measure
% of error) against the number of iterations and plot the real part of the
% recovered signal against the real part of the original signal.

% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017
%% -----------------------------START----------------------------------

clc
clear
close all

n = 100;          % Dimension of unknown vector
m = 8 * n;        % Number of measurements
isComplex = true; % If the signal and measurements are complex

%%  Build a random test problem
fprintf('Building test problem...\n');
[A, xt, b0] = buildTestProblem(m, n, isComplex);

% Options
opts = struct;
opts.initMethod = 'Truncatedspectral';
% opts.initMethod = 'optimalspectral';
opts.algorithm = 'CoordinateDescent';
opts.isComplex = isComplex;
opts.maxIters = 500000;
opts.tol = 1e-6;
opts.verbose = 2;

%% Try to recover x
fprintf('Running algorithm...\n');
[x, outs, opts] = solvePhaseRetrieval(A, A', b0, n, opts);

%% Determine the optimal phase rotation so that the recovered solution
%  matches the true solution as well as possible.
alpha = (x'*xt)/(x'*x);
x = alpha * x;

%% Determine the relative reconstruction error.  If the true signal was
%  recovered, the error should be very small - on the order of the numerical
%  accuracy of the solver.
reconError = norm(xt-x)/norm(xt);
fprintf('relative recon error = %d\n', reconError);

% Plot a graph of error(definition depends on if opts.xt is provided) versus
% the number of iterations.
plotErrorConvergence(outs, opts)

% Plot a graph of the recovered signal x against the true signal xt.
plotRecoveredVSOriginal(x,xt);

author: tomgoldstein

repo: phasepack-matlab


text:
%  INTERNAL FUNCTION: generate candidates for optimization
%
%  ::
%
%    [x,f,viol,funevals]=generate_candidates(objective,lb,ub,n,...
%        restrictions,penalty,varargin)
%
%  Args:
%     - **objective** [function_handle]: objective to minimize
%     - **lb** [vector]: lower bound of the search space
%     - **ub** [vector]: upper bound of the search space
%     - **n** [integer]: number of candidates to generate
%     - **max_trials** [integer]: number of trials after which the procedure
%       crashes
%     - **restrictions** [empty|function_handle]: function evaluating the
%       violations
%     - **opt** [struct]: structure with fields
%
%       - **restrictions_in_objective** [true|false]:
%       - **returns_retcode** [true|false]:
%       - **restrictions_same_weights** [true|{false}]:
%       - **allow_restrictions_violations** [true|{false}]:
%
%     - **penalty** [numeric]: value functions exceeding this value in
%       absolute value are assigned this value
%     - **varargin** []: additional input arguments for the objective function
%
%  Returns:
%     :
%
%     - **x** [d x n matrix]: parameter vectors
%     - **f** [row vector]: value function for each parameter vector
%     - **viol** [row vector]: violation penalty for each parameter vector
%     - **funevals** [integer]: number of function evaluations
%
%

author: jmaih


text:
% this function is dedicated for crossover operation used by MATLAB ga
% function.
function xoverKids  = SP_crossover(parents,options,NVARS, ...
    FitnessFcn,thisScore,thisPopulation)

nKids = length(parents)/2;
xoverKids = cell(nKids,1); % Normally zeros(nKids,NVARS);
index = 1;
for i=1:nKids
    parent = thisPopulation{parents(index)};
    index = index + 2;
    % Flip a section of parent1.
    p1 = ceil((length(parent) -1) * rand);
    p2 = p1 + ceil((length(parent) - p1- 1) * rand);
    child = parent;
    child(p1:p2) = fliplr(child(p1:p2));
    xoverKids{i} = child; % Normally, xoverKids(i,:);
end


text:
% Test manipulation of nested cell arrays

% Test nested extraction
A = {{1 {2 3}} {4 5} {6 7 8} {9 10}};
C_real = vec(1:10);
C = flatten_nested(A);
for i=length(C_real)
   if C{i}~=C_real(i)
      fprintf(1,'flatten_nested test failed\n');
      assert(false);
   end
end
fprintf(1,'flatten_nested test passed\n');

% Test nested element-wise modification
fnc = @(x)x + 1;
A = apply_to_nested(A,fnc);
C2 = flatten_nested(A);
for i=length(C_real)
   if C2{i}~=C_real(i)+1
      fprintf(1,'apply_to_nested test failed\n');
      assert(false);
   end
end
fprintf(1,'apply_to_nested test passed\n');

% MULTI-NESTED
% Test nested element-wise modification
A = { {{1} {2}} {{3} {4}} {{5} {6}} };
B = { {{-1} {-2}} {{-3} {-4}} {{-5} {-6}} };
C_real = zeros(6,1);
fnc = @(x,y) x + y;
A = multi_apply_to_nested(fnc,A,B);
C3 = flatten_nested(A);
for i=length(C_real)
   if C3{i}~=C_real(i)
      fprintf(1,'multi_apply_to_nested test failed\n');
      assert(false);
   end
end
fprintf(1,'multi_apply_to_nested test passed\n');

% MULTI-NESTED
% Test nested element-wise modification
A = { {{1} {2}} {{3} {4}} {{5} {6}} };
B = { {{-1} {-2}} {{-3} {-4}} {{-5} {-6}} };
C = { {{5} {10}} {{5} {10}} {{5} {12}} };
C_real = [5 10 5 10 5 12];
fnc = @(x,y,z) x + y + z;
A = multi_apply_to_nested(fnc,A,B,C);
C4 = flatten_nested(A);
for i=length(C_real)
   if C4{i}~=C_real(i)
      fprintf(1,'multi_apply_to_nested 3-arg test failed\n');
      assert(false);
   end
end
fprintf(1,'multi_apply_to_nested 3-arg test passed\n');

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
function [max_r,ind] = lp_s(p,set)
%
% Computation of the maximal magnitude of the rational ADI function over
% a discrete subset of the left complex half plane.
%
%   Calling sequence:
%
%     [max_r,ind] = lp_s(p,set)
%
%   Input:
%
%     p        vector of ADI parameters;
%     set      vector representing the discrete set.
%
%   Output:
%
%     max_r    maximal magnitude of the rational ADI function over set;
%     ind      index - maximum is attained for set(ind).
%
%
%   LYAPACK 1.0 (Thilo Penzl, Jan 1999)
%

max_r = -1;
ind = 0;

for i = 1:length(set)

  x = set(i);

  rr = 1;
  for j = 1:length(p)

    rr = rr*abs(p(j)-x)/abs(p(j)+x);

  end

  if rr > max_r

    max_r = rr;
    ind = i;

  end

end

author: Sable


text:
% tiny
%d = 10;
% medium
d = 100;
% big
%d = 500;
n = 3*d;
%n = 2*d;

if 1
  clear run
if 0
  % experiment 1
  w = randn(d,1);
  % this ensures about 16% training error
  w = sqrt(2)*w/norm(w);
  % training data
  x = randn(d,n);
  s = 1./(1+exp(-w'*x));
  y = (rand(1,n) < s);
  y = 2*y-1;
  if 0
    % test data
    xt = randn(d,1e4);
    s = 1./(1+exp(-w'*xt));
    yt = (rand(1,cols(xt)) < s);
    yt = 2*yt-1;
  end
elseif 1
  % experiment 2
  % positive data
  x = dirichlet_sample(ones(d,1),n);
  w = dirichlet_sample(ones(d,1),2);
  w = log(w(:,1)./w(:,2));

  s = 1./(1+exp(-w'*x));
  y = (rand(1,n) < s);
  y = 2*y-1;
  if 1
    % test data
    xt = dirichlet_sample(ones(d,1),n);
    s = 1./(1+exp(-w'*xt));
    yt = (rand(1,cols(xt)) < s);
    yt = 2*yt-1;
  end
else
  % Collins's generator
  x = randn(d,n);
  w = randn(d,1);
  y = sign(w'*x);
  x = x + randn(d,n)*sqrt(0.8);
end

if 0
  % shift the data
  c = 10;
  x = [x+c; ones(1,n)];
  w = [w; -c*sum(w)];
end
end
[d,n] = size(x);

i1 = find(y > 0);
i0 = find(y < 0);
figure(1)
plot(x(1,i1), x(2,i1), 'o', x(1,i0), x(2,i0), 'x')
if d == 2
  draw_line_clip(w(1),w(end),-w(2),'Color','k');
end

w0 = zeros(d,1);
xy = scale_cols(x,y);

if 0
  xyt = scale_cols(xt,yt);
  % find the best lambda
  lambdas = exp(linspace(-4,4,20));
  f = [];
  for i = 1:length(lambdas)
    w = train_newton(xy,w0,lambdas(i));
    f(i) = logProb(xyt,w);
  end
  figure(1)
  semilogx(lambdas,f)
  axis_pct;
  [dummy,i] = max(f);
  lambda = lambdas(i);
else
  lambda = 1e-2;
end
fprintf('lambda = %g\n',lambda)

if exist('run') ~= 1 | ~isfield(run,'Newton')
  disp('Newton')
  [w,run.Newton] = train_newton(xy,w0,lambda);
  wbest = w;
  ebest = run.Newton.e(end);
  % training errors
  s = 1./(1+exp(-wbest'*x));
  fprintf('%g training error\n', mean((s > 0.5) ~= (y > 0)))
end

if ~isfield(run,'Coord')
  disp('Coord')
  %[w,run.Coord] = train_newton2(xy,w0,lambda);
end
if ~isfield(run,'CG')
  disp('CG')
  [w,run.CG] = train_cg(xy,w0,lambda);
end
%disp(length(cg.run.e))
%[w,cg2.run] = train_cg2(xy,w0);
%disp(length(cg2.run.e))
if 0
  % compare by iteration
  figure(3)
  i = 1:length(cg2.run.e);
  plot(i, cg.run.e(i), i, cg2.run.e(i))
  legend('CG','CG2',4)
  return
end

if ~isfield(run,'BFGS')
  disp('BFGS')
  [w,run.BFGS] = train_bfgs(xy,w0,lambda);
end
%[w,run.lmBFGS] = train_lmbfgs(xy,w0,lambda);

%[w,run.sg] = train_sg(xy,w0);
if 0
  fprintf('max(sg.run.e) = %g\n',max(sg.run.e))
  fprintf('max(cg.run.e) = %g\n',max(cg.run.e))
  figure(1)
  plot(sg.run.step)
  figure(2)
  hold on, plot(cg.run.e,'g'), hold off
  figure(3)
  plot(sg.run.w(1,:),sg.run.w(2,:))
  hold on, plot(cg.run.w(1,:),cg.run.w(2,:),'g'), hold off
  return
end

if ~isfield(run,'FixedH')
  disp('FixedH')
  [w,run.FixedH] = train_bohning(xy,w0,lambda);
end

if ~isfield(run,'Dual') & lambda > 0
  disp('Dual')
  [w,run.Dual] = train_dual(xy,w0,lambda,1);
end
%[w,run.Dual2] = train_dual(xy,w0,lambda,1);
%[w,run.Dual2] = train_dual_cg(xy,w0,lambda);

if ~isfield(run,'MIS') & lambda == 0
  [w,run.MIS] = train_mis(xy,w0,lambda);
end
if ~isfield(run,'IS') & lambda == 0
  %[w,run.IS] = train_is(x,y,w0,lambda);
end

color.Newton = 'k';
color.CG = 'g';
color.BFGS = 'r';
color.lmBFGS = 'm';
color.cg2 = 'r';
color.sg = 'g';
color.Coord = 'm';
color.FixedH = 'c';
color.Dual = 'b';
color.Dual2 = 'm';
color.MIS = 'g';
color.IS = 'k';

linespec.Newton = 'b-.';
linespec.CG = 'g--';
linespec.Coord = 'g-';
linespec.FixedH = 'c-';
linespec.BFGS = 'r-.';
linespec.lmBFGS = 'm-.';
linespec.Dual = 'y-.';
linespec.Dual2 = 'm-.';
linespec.MIS = 'm--';
linespec.IS = 'k--';

% plot cost vs. accuracy
figure(2)
ebest = -Inf;
for f = fieldnames(run)'
  thisrun = getfield(run,char(f));
  ebest = max([ebest max(thisrun.e)]);
end
for f = fieldnames(run)'
  thisrun = getfield(run,char(f));
  %semilogx(thisrun.flops, thisrun.e, getfield(linespec,char(f)));
  thisrun.err = (ebest - thisrun.e)/n;
  run = setfield(run,char(f),thisrun);
  loglog(thisrun.flops, thisrun.err, getfield(color,char(f)));
  hold on
end
hold off
xlabel('FLOPS')
%ylabel('Log-likelihood')
ylabel('Difference from optimal log-likelihood')
axis_pct;
if 1
  ax = axis;
  ax(3) = 1e-10;
  axis(ax);
end
f = fieldnames(run);
legend(f,4)
if 0
  legend off
  f = fieldnames(run);
  h = mobile_text(f{:});
end
set(gcf,'paperpos',[0.25 2.5 8 6])
%print -dpsc expt1_d100_n300.ps
%print -dpsc expt2_d100_n300.ps
%print -dpsc expt25_d100_n300.ps
%print -dpsc expt25_d500_n1500.ps
%print -dpsc expt3_d100_n300.ps
%print -dpsc expt4_d100_n300.ps
% print -dpsc expt2_d300_n1500.ps
% save expt4.mat x y run

if 0
% cosine distance from the best solution
wbest = wbest/norm(wbest);
figure(3)
for f = fieldnames(run)'
  thisrun = getfield(run,char(f));
  thisrun.w = scale_cols(thisrun.w, 1./sqrt(sum(thisrun.w.^2,1)));
  thisrun.acc = sqrt(clip(sqdist(wbest, thisrun.w)));
  run = setfield(run,char(f),thisrun);
  loglog(thisrun.flops, thisrun.acc, getfield(color,char(f)));
  hold on
end
hold off
legend(fieldnames(run),4)
end

for f = fieldnames(run)'
  thisrun = getfield(run,char(f));
  %i = convergence(thisrun.e);
  i = min(find(thisrun.err < 1e-4));
  fprintf('%-9s %11d\n',char(f),thisrun.flops(i))
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function [C_offset] = updateTemporal_endoscope(obj, Y, smin)
%% run HALS by fixating all spatial components
% input:
%   Y:  d*T, fluorescence data
%   smin: scalar, threshold for detecting one spikes (>smin*sigma)
% output:
%   C_raw: K*T, temporal components without being deconvolved

% Author: Pengcheng Zhou, Carnegie Mellon University, adapted from Johannes

% options
maxIter = obj.options.maxIter;
if ~exist('smin', 'var') || isempty(smin)
    smin = 3;
end
%% initialization
A = obj.A;
K = size(A, 2);     % number of components
C = obj.C;
C_raw = zeros(size(C));
C_offset = zeros(K, 1);
S = zeros(size(C));
A = full(A);
U = A'*Y;
V = A'*A;
aa = diag(V);   % squares of l2 norm all all components
sn =  zeros(1, K);
kernel = obj.kernel;
kernel_pars = zeros(K, 2);

%% updating
ind_del = false(K, 1);
for miter=1:maxIter
    for k=1:K
        if ind_del
            continue;
        end
        temp = C(k, :) + (U(k, :)-V(k, :)*C)/aa(k);
        % estimate noise
        if miter==1
            sn(k) = get_noise_fft(temp);
        end

        % remove baseline
        [temp, C_offset(k)] = remove_baseline(temp, sn(k));

        % deconvolution
        if obj.options.deconv_flag
            if miter==1
                [ck, sk, kernel] = deconvCa(temp, kernel, smin, true, false, sn(k));
                kernel_pars(k, :) = kernel.pars;
            else
                kernel.pars = kernel_pars(k, :);
                [ck, sk, kernel] = deconvCa(temp, kernel, smin, false, false, sn(k));
            end
        else
            ck = max(0, temp);
        end

        % save convolution kernels and deconvolution results
        C(k, :) = ck;

        if sum(ck(2:end))==0
            ind_del(k) = true;
        end
        % save the spike count in the last iteration
        if miter==maxIter
            if obj.options.deconv_flag
                S(k, :) = sk;
            end
            C_raw(k, :) = temp;
        end
    end
end
obj.A = bsxfun(@times, A, sn);
obj.C = bsxfun(@times, C, 1./sn');
obj.C_raw = bsxfun(@times, C_raw, 1./sn');
obj.S = bsxfun(@times, S, 1./sn');
obj.P.kernel_pars = kernel_pars;
obj.delete(ind_del);

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function sR = byVertices(V,varargin)
% define a spherical region by its vertices
%
% Syntax
%   sR = sphericalRegion.byVertices(V)
%
% Input
%  V - @vector3d
%
% Output
%  sR - @sphericalRegion
%
% See also
% sphericalRegion/sphericalRegion

N = normalize(cross(V,V([2:end,1])));
alpha = zeros(size(N));

sR = sphericalRegion(N,alpha,varargin{:});

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
%try_circ_fit
%
% IB
%
% revival of a 13 years old code

  % Create data for a circle + noise

  th = linspace(0,2*pi,20)';
  R=1.1111111;
  sigma = R/10;
  x = R*cos(th)+randn(size(th))*sigma;
  y = R*sin(th)+randn(size(th))*sigma;

   plot(x,y,'o'), title(' measured points')
   pause(1)

   % reconstruct circle from data
   [xc,yc,Re,a] = circfit(x,y);
      xe = Re*cos(th)+xc; ye = Re*sin(th)+yc;

     plot(x,y,'o',[xe;xe(1)],[ye;ye(1)],'-.',R*cos(th),R*sin(th)),
     title(' measured fitted and true circles')
      legend('measured','fitted','true')
      text(xc-R*0.9,yc,sprintf('center (%g , %g );  R=%g',xc,yc,Re))
     xlabel x, ylabel y
     axis equal

repo: mcbench-benchmarks


text:
function prob_test032 ( )

%*****************************************************************************80
%
%% TEST032 tests CHI_SQUARE_CDF, CHI_SQUARE_CDF_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST032:\n' );
  fprintf ( 1, '  CHI_SQUARE_CDF evaluates the cumulative\n' );
  fprintf ( 1, '    distribution function for the chi-square central\n' );
  fprintf ( 1, '    probability density function.\n' );
  fprintf ( 1, '  CHI_SQUARE_CDF_VALUES returns some exact values.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A    X       Exact F     CHI_SQUARE_CDF(A,X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, a, x, fx ] = chi_square_cdf_values ( n_data );

    if ( n_data == 0 );
      break
    end

    fx2 = chi_square_cdf ( x, a );

    fprintf ( 1, '  %4d  %8f  %14f  %14f\n', a, x, fx, fx2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function data = reshapeData(~, data, ~)
%RESHAPEDATA   Nothing to do for SPINOP2.

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

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% EX_STOKES_BIFURCATION_RT_MP: data file for Stokes problem in a pipe with a bifurcation.

% 1) PHYSICAL DATA OF THE PROBLEM
problem_data  = struct ();

% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = 'geo_bifurcation_mp.txt';

% Type of boundary conditions for each side of the domain
problem_data.drchlt_sides = 1:3;
problem_data.nmnn_sides = [];

% Physical parameters
problem_data.viscosity = @(x, y) ones (size (x));

% Force term
problem_data.f  = @(x, y) zeros ([2, size(x)]);

% Boundary terms
problem_data.h  = @test_stokes_bifurcation_mp_h_drchlt;

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
method_data = struct ();

method_data.element_name = 'rt';     % Element type for discretization

method_data.degree       = [3 3];  % Degree of the splines
method_data.regularity   = [2 2];  % Regularity of the splines
method_data.nsub         = [5 5];  % Number of subdivisions
method_data.nquad        = [5 5];  % Points for the Gaussian quadrature rule
method_data.Cpen = 10 * (method_data.degree(1)+1);

% 3) CALL TO THE SOLVER
[geometry, msh, space_v, vel, space_p, press] = mp_solve_stokes_div_conforming (problem_data, method_data);

% 4) POST-PROCESSING
% 4.1) EXPORT TO PARAVIEW
output_file  = 'bifurcation_2d_rt_mp_deg3_reg2_sub5';
vtk_pts = {linspace(0, 1, 20), linspace(0, 1, 20)};

fprintf ('results being saved in: %s_vel.pvd and %s_press.pvd\n', output_file, output_file)
sp_to_vtk (vel, space_v, geometry, vtk_pts, sprintf ('%s_vel', output_file), {'velocity', 'divergence'}, {'value', 'divergence'})
sp_to_vtk (press, space_p, geometry, vtk_pts, sprintf ('%s_press', output_file), 'pressure')

%!test
%! problem_data  = struct ();
%! problem_data.geo_name = 'geo_bifurcation_mp.txt';
%! problem_data.drchlt_sides = 1:3;
%! problem_data.nmnn_sides = [];
%! problem_data.viscosity = @(x, y) ones (size (x));
%! problem_data.f  = @(x, y) zeros ([2, size(x)]);
%! problem_data.h  = @test_stokes_bifurcation_mp_h_drchlt;
%! method_data = struct ();
%! method_data.element_name = 'rt';     % Element type for discretization
%! method_data.degree       = [3 3];  % Degree of the splines
%! method_data.regularity   = [2 2];  % Regularity of the splines
%! method_data.nsub         = [5 5];  % Number of subdivisions
%! method_data.nquad        = [5 5];  % Points for the Gaussian quadrature rule
%! method_data.Cpen = 10 * (method_data.degree(1)+1);
%! [geometry, msh, space_v, vel, space_p, press] = mp_solve_stokes_div_conforming (problem_data, method_data);
%! assert (msh.nel, 100)
%! assert (space_v.ndof, 552)
%! assert (space_p.ndof, 256)
%! for iptc = 1:4
%!   div = sp_eval (vel(space_v.gnum{iptc}) .* space_v.dofs_ornt{iptc}', space_v.sp_patch{iptc}, geometry(iptc), [20 20], 'divergence');
%!   assert (max (abs (div(:))) < 1e-12)
%! end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [mat, dim] = spm_get_matdim(img, vx, bb)
% Voxel-to-world matrix and image dimensions from image or bbox and vox-dim
%
% FORMAT [mat, dim] = spm_get_matdim(img, vx, bb)
%
% img - filename of image to use as reference (defaults to SPM's TPM.nii)
% vx  - [1 x 3] vector of voxel dimensions (mm).
% bb  - [2 x 3] array of the min and max X, Y, and Z coordinates (mm),
%       i.e. bb = [minX minY minZ; maxX maxY maxZ].
%
% mat - [4 x 4] matrix mapping voxel coordinates to world (mm) coordinates
% dim - [1 x 3] vector of image dimensions (number of voxels)
%       (both as in output from spm_vol)
%
%       Note that the output mat will correspond to the same orientation
%       as SPM's canonical templates (transverse and vx(1) forced negative)
%       if either or both bb and vx are specified (finite), but otherwise
%       will keep the orientation of the reference image.
%__________________________________________________________________________
% Copyright (C) 2013 Wellcome Trust Centre for Neuroimaging

% Ged Ridgway
% $Id: spm_get_matdim.m 5374 2013-03-29 17:26:24Z ged $

if nargin < 3, vx  = nan(1, 3); end
if nargin < 2, bb  = nan(2, 3); end
if nargin < 1, img = '';        end

% Use MNI space by default, based on tissue priors
try
    vol = spm_vol(img);
    if isempty(vol), error('Failed to read volume %s', img), end
catch
    vol = spm_vol(fullfile(spm('dir'), 'tpm', 'TPM.nii'));
end
mat = vol(1).mat;
dim = vol(1).dim;

valid_bb = all(isfinite(bb(:)));
valid_vx = all(isfinite(vx));
if ~valid_bb && ~valid_vx
    return
end

% User has specified one or both of bb or vx, over-ride appropriately
[BB VX] = spm_get_bbox(vol(1));
if ~valid_bb, bb = BB; end
if ~valid_vx, vx = VX; end

% Determine mat and dim from bb and vx, assuming "canonical" orientation
% (i.e. transverse with negative first voxel dimension)
vx  = [-1 1 1] .* abs(vx);
mn  = vx .* min(bb ./ repmat(vx, 2, 1)); % "first" voxel's mm coordinates
mx  = vx .* round(max(bb ./ repmat(vx, 2, 1))); % "last voxel's mm coords
% matrix that maps voxel [1 1 1] to mn
mat = spm_matrix([mn 0 0 0 vx]) * spm_matrix([-1 -1 -1]);
% dim such that mat * [dim 1]' == [mx 1]'
dim = mat \ [mx 1]';
dim = round(dim(1:3)');

author: spm


text:
function proc= xvalutil_procSetDefault(proc)
% XVALUTIL_PROCSETDEFAULTS - Function used in crossvalidation

% 2014-02 Benjamin Blankertz

if isempty(proc),
  return;
end

if ~isfield(proc, 'train'),
  proc.train= {};
end

if ~isfield(proc, 'apply'),
  proc.apply= {};
end

proc.train= setDefaultProc(proc.train);
proc.apply= setDefaultProc(proc.apply);

return

function proc= setDefaultProc(proc)

for k= 1:length(proc),
  if ~iscell(proc{k}),
    proc{k}= {proc{k}};
  end
end

author: bbci


text:
function rotMat = matRad_getRotationMatrix(gantryAngle,couchAngle,system)
% matRad function to return the rotation / transformation matrix
% for gantry and/or couch rotation. The Rotation matrix stands for a (1)
% counter-clockwise, (2) active rotation in the patient coordinate system
% that is performed on a (4) column vector (by premultiplying the matrix).
% Per change of one of these directions a matrix transpose of the returned
% matrix is required.
%
%
% call
%  rotMat = matRad_getRotationMatrix(gantryAngle,couchAngle,type,system)
%
% input
%   gantryAngle:    beam/gantry angle
%   couchAngle:     couch angle
%
%   system:         optional coordinate system the transformation matrix is
%                   requested for. So far, only the default option 'LPS' is
%                   supported (right handed system).
%
% output
%   rotMat:         3x3 matrix that performs an active rotation around the
%                   patient system origin via rotMat * x
%
% References
%   https://en.wikipedia.org/wiki/Rotation_matrix (2017, Mar 1)
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Parse arguments
%We need at least two and max 3 input arguments
narginchk(2,3);

% Coordinate System (only LPS so far)
if nargin < 3
    system = 'LPS';
end

%% Set Up requested Rotation Matrix
switch system
    case 'LPS'
        %The LPS system is right-handed, gantry rotates counter-clockwise
        %around z and couch rotates counter-clockwise around y

        %active, counter-clockwise rotation Matrix for Gantry around z
        %with pre-multiplication of the matrix (R*x)
        %Note: Gantry rotation is physically an active rotation of a beam
        %vector around the target / isocenterin the patient coordinate
        %system
        R_Gantry = [cosd(gantryAngle)  -sind(gantryAngle)  0; ...
            sind(gantryAngle)    cosd(gantryAngle)  0; ...
            0                            0           1];

        %active, counter-clockwise rotation for couch around y
        %with pre-multiplication of the matrix (R*x)
        %Note: Couch rotation is physically a passive rotation of the
        %patient system around the beam target point / isocenter
        R_Couch = [cosd(couchAngle) 0 sind(couchAngle); ...
            0 1 0; ...
            -sind(couchAngle)    0    cosd(couchAngle)];
    otherwise
        error('matRad only supports LPS system so far');
end

rotMat = R_Couch*R_Gantry;

end

repo: matRad


text:
classdef NodeCoordinatesComputerTester < Tester

    properties (Access = private)
        data
    end

    properties (Access = public)
        testName
        corrValues
        calcValues
    end

    methods (Access = public)

        function obj = NodeCoordinatesComputerTester(cParams)
            obj.data = cParams;
            obj.testName = 'NodeCoordinatesComputer';
            obj.loadCorrectValues();
            obj.obtainCalculatedData();
        end

    end

    methods (Access = protected)

        function loadCorrectValues(obj)
            switch obj.data.nvert
                case 4
                    vC = load('vertCoordQuad.mat');
                    bC = load('boundCoordQuad.mat');
                    tC = load('totalCoordQuad.mat');
                case 6
                    vC = load('vertCoordHex.mat');
                    bC = load('boundCoordHex.mat');
                    tC = load('totalCoordHex.mat');
            end
            obj.corrValues(1).Matrix = vC.vertCoord;
            obj.corrValues(2).Matrix = bC.boundary;
            obj.corrValues(3).Matrix = tC.coord;
        end

        function obtainCalculatedData(obj)
            solution = NodeCoordinatesComputer(obj.data);
            solution.computeCoordinates();
            obj.calcValues(1).Matrix = solution.vertCoord;
            obj.calcValues(2).Matrix = solution.boundCoord;
            obj.calcValues(3).Matrix = solution.totalCoord;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function result = calcPreRecRadiusLabel(queryLabel, retrievalLabel, qB, rB)
%% Function: calcPreRecRadiusLabel
%   calculate precision and recall within different radius based on Label.
% Input:
%   queryLabel: 0-1 label matrix (numQuery * numLabel) for query set.
%   retrievalLabel: 0-1 label matrix (numQuery * numLabel) for retrieval set.
%   qB: compressed binary code for query set.
%   rB: compressed binary code for retrieval set.
% Output:
%   result.Pre: maxR-dims vector. Precision within different hamming radius.
%   result.Rec: maxR-dims vector. Recall within different hamming radius.
Wtrue = queryLabel * retrievalLabel' > 0;
Dhamm = hammingDist(qB, rB);

maxHamm = max(Dhamm(:));
totalGoodPairs = sum(Wtrue(:));

% find pairs with similar codes
precision = zeros(maxHamm, 1);
recall = zeros(maxHamm, 1);
for n = 1: length(precision)
    j = (Dhamm <= ((n-1) + 00.001));
    retrievalGoodPairs = sum(Wtrue(j));

    retrievalPairs = sum(j(:));
    precision(n) = retrievalGoodPairs / (retrievalPairs + eps);
    recall(n) = retrievalGoodPairs / totalGoodPairs;
end

result.Pre = precision;
result.Rec = recall;
end

author: jiangqy


text:
function Model = readCbTmodel(modelName, Tmodel, varargin)
% Creates COBRA format of a specific model from `Tmodel`.
% (composite database).
% Called by `mergeModelsBorg`, `driveModelBorgifier`, calls `TmodelFields`, `orderModelFieldsBorg`, `organizeModelCool`, `buildRxnEquations`, `ismatlab`, `parseBoolean`.
%
% USAGE:
%
%    Model = readCbTmodel(modelName, Tmodel)
%
% INPUTS:
%    modelName:           Name of an available model (i.e. iAF1260). Available names
%                         can be found by looking in `Tmodel.Models`
%    Tmodel:              Otherwise `Tmodel` will be loaded from path below.
%
% OPTIONAL INPUTS:
%    emptyArrayChoice:    "y" or "n" to remove empty arrays or not
%    'Verbose':           Print progress stamps.
%
% OUTPUTS:
%    Model:         Model in COBRA format.
%
% Please cite:
% `Sauls, J. T., & Buescher, J. M. (2014). Assimilating genome-scale
% metabolic reconstructions with modelBorgifier. Bioinformatics
% (Oxford, England), 30(7), 1036?8`. http://doi.org/10.1093/bioinformatics/btt747
%
% ..
%    Edit the above text to modify the response to help addMetInfo
%    Last Modified by GUIDE v2.5 06-Dec-2013 14:19:28
%    This file is published under Creative Commons BY-NC-SA.
%
%    Correspondance:
%    johntsauls@gmail.com
%
%    Developed at:
%    BRAIN Aktiengesellschaft
%    Microbial Production Technologies Unit
%    Quantitative Biology and Sequencing Platform
%    Darmstaeter Str. 34-36
%    64673 Zwingenberg, Germany
%    www.brain-biotech.de

verbose = false ;  % Varaging.
emptyArrayChoice = [] ;
if ~isempty(varargin)
    if sum(strcmp('Verbose', varargin))
        verbose = true ;
    end
    if sum(strcmp('y', varargin))
        emptyArrayChoice = 'y' ;
    end
    if sum(strcmp('n', varargin))
        emptyArrayChoice = 'n' ;
    end
end

%% Declare variables.
% Check to see if model exists.
TmodelNames = fieldnames(Tmodel.Models) ;
modelCheck = 0 ;
while ~modelCheck
    modelCheck = sum(strcmp(modelName, TmodelNames)) ;
    if modelCheck
        if verbose ; fprintf('Extracting %s from Tmodel\n', modelName) ; end
        modelCheck = 1 ;
    else
        fprintf('ERROR: Model %s not found in Tmodel\n', modelName)
        % Print available model names and allow user to rechoose.
        fprintf('The available models are :\n')
        for iName = 1:length(TmodelNames)
            fprintf([TmodelNames{iName} '\n'])
        end
        modelName = input('Input desired model: ', 's') ;
    end
end

% Set name.
Model.description = modelName ;

otherModels = find(cellfun(@isempty, (strfind(TmodelNames, modelName)))) ;

% Get field names.
fields = TmodelFields ;
% Reaction related cell field names.
rxnFields = fields{1} ;
% Reaction related double Fields (which are kept model specific).
rNumFields = fields{2} ;
% Metabolite related cell field names.
metFields = fields{3} ;
% Metabolite related double field names.
mNumFields = fields{4} ;
% All fields that will end up in the Model.
allFields = fields{5} ;

% Identity arrays for desired model.
rxnLog = Tmodel.Models.(modelName).rxns ;
metLog = Tmodel.Models.(modelName).mets ;

%% Extract reaction information.
for iF = 1:length(rxnFields)
    Model.(rxnFields{iF}) = Tmodel.(rxnFields{iF})(rxnLog) ;
end
for iF = 1:length(rNumFields)
    Model.(rNumFields{iF}) = Tmodel.(rNumFields{iF}).(modelName)(rxnLog);
end

% Remove information from other models for grRules.
% Find where the rull starts for the model.
% searchStringStart = ['(?<=(^|\|)' modelName ':).*'] ;
grRulesStart = regexp(Model.grRules, modelName ) ; %searchStringStart) ;
fixGrRulesStart = @(startpos, shift) not(isempty(startpos)) * (startpos + shift) ;
grRulesStart = cellfun(@(startpos, shift)fixGrRulesStart(startpos, shift), ...
                        grRulesStart, repmat({length(modelName) + 1}, ...
                        size(grRulesStart)), 'Uniformoutput',false) ;
ruleIndexes = find(~cellfun(@isempty,grRulesStart)) ;
% Find where information for other models starts, if it is after the
% information for this model.
searchStringEnd = ['(?<=' modelName '.*).('] ;
for iMo = 1:length(otherModels)
    if iMo ~= 1
        searchStringEnd = [searchStringEnd '|'] ;
    end
    searchStringEnd = [searchStringEnd '\|'...
                        TmodelNames{otherModels(iMo)}] ;
end
searchStringEnd = [searchStringEnd '):'] ;
grRulesEnd = regexpi(Model.grRules,searchStringEnd, 'once') ;
% Clear gene rules and rebuild.
grRulesHold = Model.grRules ;
Model.grRules(:) = {''} ;
% Only include information from between the index where information starts
% and where it end.
for iRxn = 1:length(ruleIndexes)
    % if there was no ending index, set it as the length of the string.
    if isempty(grRulesEnd{ruleIndexes(iRxn)})
        grRulesEnd{ruleIndexes(iRxn)} = ...
            length(grRulesHold{ruleIndexes(iRxn)}) ;
    end
    Model.grRules{ruleIndexes(iRxn)} = ...
        grRulesHold{ruleIndexes(iRxn)} ...
            (grRulesStart{ruleIndexes(iRxn)}: ...
            grRulesEnd{ruleIndexes(iRxn)}) ;
end

% Clean up rxnID also. This is really key, beacause we need the original
% IDs to compare the removed model to the model before it was merged.
if ismatlab
    searchString = ['(?<=(^|\|)' modelName ':)[^\|]*'] ;
    IDs = regexpi(Model.rxnID, searchString, 'match') ;
else
    for ir = 1:length(Model.rxnID)
        nowID = Model.rxnID{ir} ;
        startpos = min(strfind(nowID, modelName) + length(modelName) + 1) ;
        endpos =   strfind(nowID, '|') -1 ;
        endpos = min(endpos(endpos > startpos)) ;
        IDs{ir} = nowID(startpos:endpos) ;
    end
end
IDIndexes = find(~cellfun(@isempty, IDs)) ;
for iRxn = 1:length(IDIndexes)
    Model.rxnID{IDIndexes(iRxn)} = char(IDs{iRxn}) ;
end

%% Extract metabolite information.
Model.mets = Tmodel.mets(metLog) ;
for iF = 1:length(metFields)
    Model.(metFields{iF}) = Tmodel.(metFields{iF})(metLog) ;
end
for iF = 1:length(mNumFields)
    Model.(mNumFields{iF}) = Tmodel.(mNumFields{iF})(metLog) ;
end

% Clean up metID
if ismatlab
    searchString = ['(?<=(^|\|)' modelName ':)[^\|]*'] ;
    IDs = regexpi(Model.metID, searchString, 'match') ;
else
    for ir = 1:length(Model.metID)
        nowID = Model.metID{ir} ;
        startpos = min(strfind(nowID, modelName) + length(modelName) + 1) ;
        endpos =   strfind(nowID, '|') -1 ;
        endpos = min(endpos(endpos > startpos)) ;
        IDs{ir} = nowID(startpos:endpos) ;
    end
end
IDIndexes = find(~cellfun(@isempty, IDs)) ;
for iMet = 1:length(IDIndexes)
    Model.metID{IDIndexes(iMet)} = char(IDs{iMet}) ;
end

%% Construct S Matrix.
Model.S = Tmodel.S(metLog, rxnLog) ;

%% Rebuild reaction equations (as bounds are different per Model)
Model = buildRxnEquations(Model) ;

%% Extract gene information.
if isfield(Tmodel.Models.(modelName), 'genes')

    nRxns = length(Model.rxns) ;

    allGenes = Tmodel.Models.(modelName).genes ;
    Model.genes = allGenes ; % put in model for saving.
    genesByRxn = cell(nRxns, 1) ;

    % Construct rules cell array
    rules = cell(nRxns, 1) ;
    rules(:) = {''} ;

    % Convert grRules to rules format and pull out genes per reaction.
    for iRxn = 1:nRxns
       [genesByRxn{iRxn}, rules{iRxn}] = parseBoolean(Model.grRules{iRxn});
    end

    % Construct gene to rxn mapping.
    rxnGeneMat = sparse(nRxns,length(allGenes)) ;

    % but only make the rules if there are grRules.
    if sum(rxnGeneMat) > 0
        fprintf('Building rxnGeneMat.\n');
        for i = 1:nRxns
            if iscell(genesByRxn{i})
                [~, geneInd] = ismember(genesByRxn{i}, allGenes);
            else
                [~, geneInd] = ismember(num2cell(genesByRxn{i}), allGenes);
            end

            if geneInd > 0
                rxnGeneMat(i, geneInd) = 1;
                for j = 1:length(geneInd)
                    rules{i} = strrep(rules{i}, ['x(' num2str(j) ')'], ...
                        ['x(' num2str(geneInd(j)) '_TMP_)']);
                end
                rules{i} = strrep(rules{i}, '_TMP_', '');
            else
                rules{i} = '' ;
            end
        end

        % Put things in the model.
        Model.rules = rules ;
        Model.rxnGeneMat = rxnGeneMat ;
    end
end

%% Remove all empty cell arrays.
if isempty(emptyArrayChoice)
    choice = input('Remove empty cell arrays? (y/n): ', 's') ;
else
    choice = emptyArrayChoice ;
end
if strcmpi(choice,'y') || strcmpi(choice, 'yes')
    if verbose ; fprintf('Removing empty cell arrays:\n') ; end
    for iF = 1:length(allFields)
        % Only remove the cell array fields.
        if isfield(Model, allFields{iF})
            if iscell(Model.(allFields{iF}))
                if cellfun(@isempty, Model.(allFields{iF}))
                    Model = rmfield(Model, allFields{iF}) ;
                    if verbose ; fprintf([allFields{iF} '\n']) ; end
                end
            end
        end
    end
end

%% Order model.
Model = orderModelFieldsBorg(Model) ;
Model = organizeModelCool(Model) ;

author: opencobra

repo: cobratoolbox