text:
function ssmult_install (dotests)
%SSMULT_INSTALL compiles, installs, and tests ssmult.
% Note that the "lcc" compiler provided with MATLAB for Windows can generate
% slow code; use another compiler if possible.  Your current directory must be
% SSMULT for ssmult_install to work properly.  If you use Linux/Unix/Mac, I
% recommend that you use COPTIMFLAGS='-O3 -DNDEBUG' in your mexopts.sh file.
%
% Example:
%   ssmult_install          % compile and install
%   ssmult_install (0)      % just compile and install, do not test
%
% See also ssmult, ssmultsym, sstest, sstest2, mtimes.

% Copyright 2009, Timothy A. Davis, University of Florida

fprintf ('Compiling SSMULT:\n') ;

%-------------------------------------------------------------------------------
% compile ssmult and add it to the path
%-------------------------------------------------------------------------------

d = '' ;
if (~isempty (strfind (computer, '64')))
    % 64-bit MATLAB
    d = ' -largeArrayDims -DIS64' ;
end

v = getversion ;
if (v < 6.5)
    % mxIsDouble is false for a double sparse matrix in MATLAB 6.1 or earlier
    d = [d ' -DMATLAB_6p1_OR_EARLIER'] ;
end

cmd = sprintf ('mex -O%s ssmult.c ssmult_mex.c ssmult_saxpy.c ssmult_dot.c ssmult_transpose.c', d) ;
disp (cmd) ;
eval (cmd) ;

cmd = sprintf ('mex -O%s ssmultsym.c', d) ;
disp (cmd) ;
eval (cmd) ;

cmd = sprintf ('mex -O%s sptranspose.c ssmult_transpose.c', d) ;
disp (cmd) ;
eval (cmd) ;

addpath (pwd) ;
fprintf ('\nssmult has been compiled, and the following directory has been\n') ;
fprintf ('added to your MATLAB path.  Use pathtool to add it permanently:\n') ;
fprintf ('\n%s\n\n', pwd) ;
fprintf ('If you cannot save your path with pathtool, add the following\n') ;
fprintf ('to your MATLAB startup.m file (type "doc startup" for help):\n') ;
fprintf ('\naddpath (''%s'') ;\n\n', pwd) ;

%-------------------------------------------------------------------------------
% test ssmult and ssmultsym
%-------------------------------------------------------------------------------

if (nargin < 1)
    dotests = 1 ;
end
if (~dotests)
    return
end

fprintf ('Please wait while your new ssmult function is tested ...\n') ;

fprintf ('\nTesting large sparse column vectors (1e7-by-1)\n') ;
x = sprandn (1e7,1,1e-4) ;
y = sprandn (1e7,1,1e-4) ;
x (1) = pi ;
y (1) = exp (1) ;
tic ; a = x'*y ; t1 = toc ;
tic ; b = ssmult (x, y, 1) ; t2 = toc ;
fprintf ('s=x''*y in MATLAB: %8.3f seconds\n', t1) ;
fprintf ('s=ssmult(x,y,1):  %8.3f seconds; error %g\n', t2, abs (full(a-b))) ;
fprintf ('SSMULT speedup: %8.3g\n\n', t1/t2) ;

load west0479
A = west0479 ;
B = sprand (A) ;
C = A*B ;
D = ssmult (A,B) ;
err = norm (C-D,1) / norm (C,1) ;
fprintf ('west0479 error: %g\n', err) ;

fprintf ('\ntesting large matrices (may fail if you are low on memory):\n')
rand ('state', 0) ;

n = 10000 ;
A = sprand (n, n, 0.01) ;
B = sprand (n, n, 0.001) ;
test_large (A,B) ;

msg = { 'real', 'complex' } ;

% all of these calls to ssmult should fail:
fprintf ('\ntesting error handling (the errors below are expected):\n') ;
A = { 3, 'gunk', sparse(1), sparse(1), sparse(rand(3,2)) } ;
B = { 4,   0   , 5,         msg,       sparse(rand(3,4)) } ;
for k = 1:length(A)
    try
        % the following statement is supposed to fail
        C = ssmult (A {k}, B {k}) ;                                         %#ok
        error ('test failed\n') ;
    catch
        disp (lasterr) ;
    end
end
fprintf ('error handling tests: ok.\n') ;

% err should be zero:
rand ('state', 0)
for Acomplex = 0:1
    for Bcomplex = 0:1
        err = 0 ;
        fprintf ('\ntesting C = A*B where A is %s, B is %s\n', ...
            msg {Acomplex+1}, msg {Bcomplex+1}) ;
        for m = [ 0:30 100 ]
            fprintf ('.') ;
            for n = [ 0:30 100 ]
                for k = [ 0:30 100 ]
                    A = sprand (m,k,0.1) ;
                    if (Acomplex)
                        A = A + 1i*sprand (A) ;
                    end
                    B = sprand (k,n,0.1) ;
                    if (Bcomplex)
                        B = B + 1i*sprand (B) ;
                    end
                    C = A*B ;
                    D = ssmult (A,B) ;
                    s = ssmultsym (A,B) ;
                    err = max (err, norm (C-D,1)) ;
                    err = max (err, nnz (C-D)) ;
                    err = max (err, isreal (D) ~= (norm (imag (D), 1) == 0)) ;
                    err = max (err, s.nz > nnz (C)) ;
                    [i j x] = find (D) ;                                    %#ok
                    if (~isempty (x))
                        err = max (err, any (x == 0)) ;
                    end
                end
            end
        end
        fprintf (' maximum error: %g\n', err) ;
    end
end

sstest ;
fprintf ('\nSSMULT tests complete.\n') ;

%-------------------------------------------------------------------------------
function [v,pc] = getversion
% determine the MATLAB version, and return it as a double.
% only the primary and secondary version numbers are kept.
% MATLAB 7.0.4 becomes 7.0, version 6.5.2 becomes 6.5, etc.
v = version ;
t = find (v == '.') ;
if (length (t) > 1)
    v = v (1:(t(2)-1)) ;
end
v = str2double (v) ;
try
    % ispc does not appear in MATLAB 5.3
    pc = ispc ;
catch
    % if ispc fails, assume we are on a Windows PC if it's not unix
    pc = ~isunix ;
end

%-------------------------------------------------------------------------------
function test_large (A,B)
% test large matrices
n = size (A,1) ;
fprintf ('dimension %d   nnz(A): %d   nnz(B): %d\n', n, nnz (A), nnz (B)) ;
c = ssmultsym (A,B) ;
fprintf ('nnz(C): %d   flops: %g   memory: %g MB\n', ...
    c.nz, c.flops, c.memory/2^20) ;
try
    % warmup for accurate timings
    C = A*B ;                                                               %#ok
    D = ssmult (A,B) ;                                                      %#ok
    tic ;
    C = A*B ;
    t1 = toc ;
    tic ;
    D = ssmult (A,B) ;
    t2 = toc ;
    tic ;
    t3 = toc ;
    fprintf ('MATLAB time:          %g\n', t1) ;
    err = norm (C-D,1) ;
    fprintf ('SSMULT time:          %g err: %g\n', t2, err) ;
catch
    disp (lasterr)
    fprintf ('tests with large random matrices failed ...\n') ;
end
clear C D

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef grain2d < phaseList & dynProp
  % class representing two dimensional grains
  %
  % Syntax
  %   grains = grain2d(V, poly, ori, CSList, phaseId, phaseMap)
  %
  % Input
  %  V    - n x 2 list of vertices
  %  poly - cell array of the polyhedrons
  %  ori  - array of mean orientations
  %  CSList   - cell array of symmetries
  %  phaseId  - list of phaseId for each grain
  %  phaseMap -
  %
  % Example
  %
  %   V = [0 0; 1 0; 2 0; 0 1; 1 1; 2 1];
  %   poly = {[1 2 5 4 1];[2 3 6 5 2]};
  %   rot = rotation.rand(2,1);
  %   grains = grain2d(V,poly,rot,crystalSymmetry.load('quartz'))
  %   plot(grains,grains.meanOrientation)
  %
  % Class Properties
  %  phaseId - phase identifier of each grain
  %  id            - id of each grain
  %  poly          - cell list of the vertex ids of each grain (index to V)
  %  V             - list of verticies (x,y coordinates)
  %  boundary      - @grainBoundary
  %  innerBoundary - @grainBoundary
  %  triplePoints  - @triplePoints
  %  grainSize     - number if pixels belonging to the grain
  %  GOS           - grain orientation spread
  %  meanOrientation - average grain orientation (<GrainOrientationParameters.html only single phase>)
  %
  % See also
  % GrainReconstruction GrainSpatialPlots SelectingGrains ShapeParameter

  % properties with as many rows as data
  properties
    poly={}    % cell list of polygons forming the grains
    id=[]      % id of each grain
    grainSize = [] % number of measurements per grain
  end

  properties (Hidden = true)
    inclusionId = []; % number of elements in poly that model inclusions
  end

  % general properties
  properties
    boundary = grainBoundary % boundary of the grains
    innerBoundary = grainBoundary % inner grain boundary
  end

  properties (Dependent = true)
    meanOrientation  % mean orientation
    V                % vertices with x,y coordinates
    scanUnit         % unit of the vertice coordinates
    GOS              % intragranular average misorientation angle
    x                % x coordinates of the vertices of the grains
    y                % y coordinates of the vertices of the grains
    triplePoints     % triple points
  end

  properties (Dependent = true, Access = protected)
    idV % active vertices
  end

  methods

    function grains = grain2d(V, poly, ori, CSList, phaseId,  phaseMap, varargin)
      % constructor
      %
      % Input
      %  V    - n x 2 list of vertices
      %  poly - cell array of the polyhedrons
      %  ori  - array of mean orientations
      %  CSList   - cell array of symmetries
      %  phaseId  - list of phaseId for each grain
      %  phaseMap -

      if nargin == 0, return;end

      grains.poly = poly;
      grains.inclusionId = cellfun(@(p) length(p) - find(p(2:end)==p(1),1),poly)-1;

      if nargin>=3 && ~isempty(ori)
        grains.prop.meanRotation = ori;
      else
        grains.prop.meanRotation = rotation.nan(length(poly),1);
      end

      if nargin>=4
        grains.CSList = ensurecell(CSList);
      else
        grains.CSList = {'notIndexed'};
      end

      if nargin>=5
        grains.phaseId = phaseId;
      else
        grains.phaseId = ones(length(poly),1);
      end

      if nargin>=6
        grains.phaseMap = phaseMap;
      else
        grains.phaseMap = 1:length(grains.CSList);
      end

      grains.id = (1:numel(grains.phaseId)).';
      grains.grainSize = ones(size(poly));

      if isa(V,'grainBoundary') % grain boundary already given
        grains.boundary = V;
      else % otherwise compute grain boundary

        % compute boundary segments
        F = arrayfun(@(i) poly{i}([1:end-grains.inclusionId(i)-1;2:end-grains.inclusionId(i)]),...
          1:length(poly),'uniformOutput',false);
        F = [F{:}].';

        lBnd = cellfun(@(p) length(p),poly) - grains.inclusionId -1;

        grainIds = repelem(1:length(grains),lBnd);

        [~,iF,iG] = unique(sort(F,2),'rows','stable');
        F = F(iF,:);

        % compute boundary grainIds
        grainId = zeros(size(F));
        grainId(:,1) = grainIds(iF); % first column - first apperance

        % second column - second appearance
        col2 = true(size(iG,1),1);
        col2(iF) = false;
        grainId(iG(col2),2) = grainIds(col2);

        % compute misorientation from mean orientations of the grain
        mori = rotation.nan(size(F,1),1);
        isNotBoundary = all(grainId,2);
        mori(isNotBoundary) = ...
          inv(grains.prop.meanRotation(grainId(isNotBoundary,2))) ...
          .* grains.prop.meanRotation(grainId(isNotBoundary,1));

        grains.boundary = grainBoundary(V,F,grainId,1:max(grainId),...
          grains.phaseId,mori,grains.CSList,grains.phaseMap);

      end

    end

    function V = get.V(grains)
      V = grains.boundary.V;
    end

    function x = get.x(grains)
      x = grains.boundary.x;
    end

    function y = get.y(grains)
      y = grains.boundary.y;
    end

    function grains = set.V(grains,V)

      grains.boundary.V = V;
      grains.innerBoundary.V = V;

    end

    function idV = get.idV(grains)

      isCell = cellfun('isclass',grains.poly,'cell');
      polygons = grains.poly;
      polygons(isCell) = cellfun(@(x) [x{:}] ,grains.poly(isCell),'UniformOutput',false);
      idV = unique([polygons{:}]);

    end

    function varargout = size(grains,varargin)
      [varargout{1:nargout}] = size(grains.id,varargin{:});
    end

    function ori = get.meanOrientation(grains)
      if isempty(grains)
        ori = orientation;
      else
        ori = orientation(grains.prop.meanRotation,grains.CS);

        % set not indexed orientations to nan
        if ~all(grains.isIndexed), ori(~grains.isIndexed) = NaN; end
      end
    end

    function grains = set.meanOrientation(grains,ori)

      if ~isempty(grains)

        if isnumeric(ori) && all(isnan(ori(:)))
          grains.prop.meanRotation = rotation.nan(size(grains.prop.meanRotation));
        else
          % update rotation
          grains.prop.meanRotation = rotation(ori);

          % update phase
          grains.CS = ori.CS;
        end
      end

    end

    function gos = get.GOS(grains)
      gos = grains.prop.GOS;
    end

    function unit = get.scanUnit(grains)
      unit = grains.boundary.scanUnit;
    end

    function grains = set.scanUnit(grains,unit)
      grains.boundary.scanUnit = unit;
      grains.innerBoundary.scanUnit = unit;
    end

    function tP = get.triplePoints(grains)
      tP = grains.boundary.triplePoints;
    end

    function grains = set.triplePoints(grains,tP)
      grains.boundary.triplePoints = tP;
    end

    function grains = update(grains)

      grains.boundary = grains.boundary.update(grains);
      grains.innerBoundary = grains.innerBoundary.update(grains);
      grains.triplePoints = grains.triplePoints.update(grains);

    end

  end

  methods (Static = true)

    [grain2d] = load(fname)

  end

end

author: mtex-toolbox

repo: mtex

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


text:
function hP = plugin_board()
% Zynq Platform PCore
% Use Plugin API to create board plugin object

%   Copyright 2015 The MathWorks, Inc.

% Call the common board definition function
hP = AnalogDevices.fmcomms2.common.plugin_board('ZED');

author: analogdevicesinc


text:
function MimSave(file_path, filename_base, image_data, compression, reporting)
    % MimSave. Saves an image or a structure containing images, with compression support
    %
    %    MimSave allows you to save a structure as a .mat file, where every
    %    PTKImage in the structure has its raw data saved as a separate
    %    .raw file, with compression if supported. Use MimLoad to reload
    %    the structure.
    %
    %     Syntax
    %     ------
    %
    %         MimSave(file_path, filename_base, image_data, compression, reporting)
    %
    %             file_path       specify the location to save the files.
    %             filename_base   specify the file prefix. Suffixes will be
    %                             added automatically. A single header file
    %                             will be saved for the structure or image,
    %                             and a separate raw image file will be
    %                             created for each image in the structure
    %             image_data      is a PTKImage (or PTKDicomImage) class containing the image
    %                             to be saved, or a structure which could
    %                             contain one or more PTKImages
    %             compression     the compression to use when saving the
    %                             raw data in PTKImage files
    %             reporting       an object implementing CoreReportingInterface
    %                             for reporting progress and warnings
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    result = ConvertStructAndSaveRawImageData(image_data, file_path, filename_base, 0, compression, reporting);
    filename = [fullfile(file_path, filename_base) '.mat'];
    MimDiskUtilities.Save(filename, result);
end

function [new_structure, next_index] = ConvertStructAndSaveRawImageData(old_structure, file_path, filename_base, next_index, compression, reporting)
    if isstruct(old_structure)
        new_structure = struct;
        field_names = fieldnames(old_structure);
        for field = field_names'
            field_name = field{1};
            [new_structure.(field_name), next_index] = ConvertStructAndSaveRawImageData(old_structure.(field_name), file_path, filename_base, next_index, compression, reporting);
        end
    else
        if isa(old_structure, 'PTKImage')
            reporting.LogVerbose(['Saving raw image data for ' filename_base]);
            if next_index == 0
                file_suffix = '';
            else
                file_suffix = ['_' int2str(next_index)];
            end
            raw_image_file_name = [filename_base file_suffix '.raw'];

            if (length(size(old_structure.RawImage)) ~= 3) || (~MimDiskUtilities.CompressionSupported(compression, class(old_structure.RawImage), reporting))
                % Compression currently only supports 3D images
                compression = [];
            end

            % Create a header file if requested. The header is the image object
            % minus the raw image data, and contains the raw image filename
            header = old_structure.CreateHeader(raw_image_file_name, compression);

            % Save the pixel data with compression, if supported
            CoreSaveRawImage(old_structure.RawImage, file_path, raw_image_file_name, compression, reporting);

            next_index = next_index + 1;
            new_structure = header;
        else
            new_structure = old_structure;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [ z, rnorm ] = cg_it ( A, x )

%*****************************************************************************80
%
%% CG_IT carries out the conjugate gradient iteration.
%
%  Modified:
%
%    08 March 2010
%
%  Author:
%
%    Gaurav Sharma, Jos Martin
%
%  Reference:
%
%    Gaurav Sharma, Jos Martin,
%    MATLAB: A Language for Parallel Computing,
%    International Journal of Parallel Programming,
%    Volume 37, Number 1, pages 3-36, February 2009.
%
%  Parameters:
%
%    Input, real A(N,N), the matrix.
%
%    Input, real X(N,1), the right hand side.
%
%    Output, real Z(N,1), the estimated solution.
%
%    Output, real RNORM, the norm of the residual ( x - A * z );
%
  z = zeros ( size ( x ) );
  r = x;
  rho = r' * r;
  p = r;

  for i = 1 : 15
    q = A * p;
    alpha = rho / ( p' * q );
    z = z + alpha * p;
    rho0 = rho;
    r = r - alpha * q;
    rho = r' * r;
    beta = rho / rho0;
    p = r + beta * p;
  end

  rnorm =  norm ( x - A * z );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [F,feasible] = filter_normball(F_xw,Zmodel,x,w,allw,norm_p,ops,VariableType)

% Creates robustified version of the uncertain set of linear inequalities
% s.t A(w)*x <= b(w) for all norm(w,p)<r

feasible = 1;
if length(F_xw) == 0
    F = [];
    return
end

X = sdpvar(F_xw);

% Some pre-calc. Consider uncertainties not in this class as decision
% variables (i.e. other than the ones we are robustifying)
vars = [getvariables(x) depends(X) getvariables(allw)];
vars = setdiff(vars,getvariables(w));
x = recover(vars);

% Trying to cope with cases as x+norm(w,1)+norm(w,2)<1, s.t norm(w,2)<1
% In test_robust_5
%extended_in_w = intersect(VariableType.w_variables,yalmip('extvariables'));
%extended_in_F = intersect(vars,yalmip('extvariables'));
if ~(strcmp(ops.robust.auxreduce,'none') | strcmp(ops.robust.auxreduce,'affine'))
    if ~isempty(intersect(vars,intersect(VariableType.w_variables,yalmip('extvariables'))))
        F = F_xw;
        return
    end
end

% Create a bilinear decomposition of the constraints
try
xw = [x;w];
xind = find(ismembcYALMIP(getvariables(xw),getvariables(x)));
wind = find(ismembcYALMIP(getvariables(xw),getvariables(w)));
[Qs,cs,fs,dummy,nonquadratic] = vecquaddecomp(X,xw);
all_f = [];
all_c_w = [];
all_c_x = [];
all_Q_xw = [];
HigherOrder = [];
top = 1;
for i = 1:length(X)
    Q = Qs{i};
    c = cs{i};
    f = fs{i};
    if nonquadratic
        HigherOrder = [HigherOrder;i];
        %error('Constraints can be at most quadratic, with the linear term uncertain');
    else
        Q_ww = Q(wind,wind);
        Q_xw = Q(xind,wind);
        Q_xx{top} = Q(xind,xind);
        c_x = c(xind);
        c_w = c(wind);

        all_f = [all_f;f];
        all_c_w = [all_c_w;c_w'];
        all_c_x = [all_c_x;sparse(c_x)];
        all_Q_xw = [all_Q_xw Q_xw'];
        top = top + 1;
    end
end
% Linear uncertain constraint is (Bbetai*x + cdi) >= 0 for all w, or
% (bi' + (Bi*w)')*x + (ci'*w + di).
catch
    % Complicated epression, so go for polynomial decomposition below
    Q = [];
end

% Currently 3 special cases implemented
if ops.verbose
    disp([' - Eliminating uncertainty using explicit maximization of ' norm_p])
end

if length(Q)> 0
    switch norm_p
        case '1-norm'
            [F,feasible] = filter_norm_1(all_f,all_c_w,all_c_x,all_Q_xw,x,Zmodel,allw,X,Q_xx,VariableType);

        case '2-norm'
            [F,feasible] = filter_norm_2(all_f,all_c_w,all_c_x,all_Q_xw,x,Zmodel,allw,X,Q_xx,VariableType);

        case 'inf-norm'
            [F,feasible] = filter_norm_inf(all_f,all_c_w,all_c_x,all_Q_xw,x,Zmodel,allw,X,Q_xx,VariableType);

        otherwise
            error('The detected norm-ball has not been implemented yet')
    end
elseif isequal(norm_p,'inf-norm')
    % Quadratic decomposition failed. Let's try a general basis if simple
    % inf-norm
    F = [];
    feasible = 1;
    HigherOrder = [];
   % [c] = coefficients([1+sum(allw);X],allw);
    for i = 1:length(X)
       % try
            [c,v] = coefficients(X(i),allw,[1;allw]);
            % Great, affine in uncertainty
            if isa(c(2:end),'double') && nnz(c(2:end))==0
                % Trivial case, no uncertainty in model
                F = [F, X(i) >= 0];
            else
                % c(1) +  c(2:end)'*w >= 0
                % c(1) +  c(2:end)'*(center + r*wnormalized) >= 0
                % c(1) +  c(2:end)'*center + r*c(2:end)'*wnormalized >= 0
                center = (Zmodel.lb + Zmodel.ub)/2;
                r = (Zmodel.ub - Zmodel.lb)/2;
                t = sdpvar(length(allw),1);
                F = [F, c(1) + c(2:end)'*center - sum(t) >= 0, -t <= r.*c(2:end) <= t];
            end
      %  catch
      %      HigherOrder = [HigherOrder;i];
      %  end
    end
else
    F = [];
    feasible = 1;
end

if ~isempty(HigherOrder)
    F = [F,X(HigherOrder) >= 0];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
s = [1 1 1 1 1]';
t = [2 3 4 5 6]';
G = graph(s',t');
I = full(incidence(G))

n = length(s);

ind = (1:n)'
A = full(sparse([s;t],[ind;ind],[ones(n,1)*-1;ones(n,1)]));

A

I

A - I

s = [1 1 1 1 1]';
t = [6 3 4 5 2]';
G = graph(s',t');
I = full(incidence(G))

ind = (1:n)'
A = full(sparse([s;t],[ind;ind],[ones(n,1)*-1;ones(n,1)]));

A

I

A - I

s = [1 1 1 1 1 6]';
n = length(s);
t = [6 3 4 5 2 1]';
G = graph(s',t');%switches orientation of duplicates
I = full(incidence(G))

ind = (1:n)'
A = full(sparse([s;t],[ind;ind],[ones(n,1)*-1;ones(n,1)]));

A

I

A - I

s = [1 1 1 1 1 6]';
n = length(s);
t = [6 3 4 5 2 1]';
G = digraph(s',t'); %does not switch orientation of duplicates
I = full(incidence(G))

ind = (1:n)'
A = full(sparse([s;t],[ind;ind],[ones(n,1)*-1;ones(n,1)]));

A

I

A - I

author: opencobra

repo: cobratoolbox


text:
function chi2 = chi_squared(y,fit,P,eb)

% returns *reduced* chi^2 value for use in data modelling
% "y" is a vector of data, "fit" is a vector of model values (size(fit)=size(y)), P is the number of
% parameters fit in the model, and eb is a vector of error bars (1-to-1 correspondnce with y)
% Ref: John R. Taylor, "An Introduction to Error Analysis", (2nd ed., 1997)
% 11/11/01 Mike Scarpulla.  Please direct questions or comments to scarps@uclink.berkeley.edu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<3
    error('Wrong number of arguments passed to "chi_squared"')
end

% if error bars are not availible, evaluate chi^2 by normalizing deviation^2 by magnitude of data.
% This assumes that the STDEV of a value scales as SQRT(value).  USE WITH THIS CAVEAT IN MIND
if nargin==3
    N = max(size(y));
    terms = ((y-fit).^2)./abs(y);
    chi2 = 1/(N-P)*sum(terms);
end

%if error bars are availible, normalize the deviation to the expectred error
if nargin==4
    N = max(size(y));
    terms = ((y-fit)./eb).^2;
    chi2 = 1/(N-P)*sum(terms);
end

author: Sable

repo: mcbench-benchmarks


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  ETHZ - Computer Vision Lab
%
%  Jordi Pont-Tuset <jponttuset@vision.ee.ethz.ch>
%  September 2015
% ------------------------------------------------------------------------
% This file is part of the BOP package presented in:
%    Pont-Tuset J, Van Gool, Luc,
%    "Boosting Object Proposals: From Pascal to COCO"
%    International Conference on Computer Vision (ICCV) 2015.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
function db_root_dir = db_root_dir( database )
if strcmp(database,'Pascal') || strcmp(database,'COCO') ||  strcmp(database,'SBD') ||  strcmp(database,'ILSVRC') ||  strcmp(database,'BSDS500') || strcmp(database,'bsds_object_gt')
    db_root_dir = fullfile('/srv/glusterfs/jpont/gt_dbs/',database);
else
    error(['Unknown database: ' database]);
end

end

author: jponttuset

repo: mcg

save_path: github-repos/MATLAB/jponttuset-mcg


text:
%% example
function [Y, X, U] = example(U, varargin)

%% Syntax
%  function [Y, X, U] = example(U, varargin)

%% Description
% Function for simulating the nonlinear dynamic system:
%
%                y(k)
%   y(k+1) = ---------------  + [u(k)]^3
%            1 + y(k)*y(k)
%
% which is used to demonstrate the use of this toolbox.
% [1] K.S. Narendra and K. Parthasarathy. Identification
% and Control of Dynamical Systems Using Neural Networks,
% IEEE Transactions on NN, Vol.1 No. 1, 4-27, 1990.
%
% Usage:
% (1) [Y, X, U] = example(u, y0), u is vector, y0 is optional
% (2) [Y, X, U] = example(U, y0), U is Nx2 matrix
% Inputs:
% u .. vector of inputs
% U .. matrix U=[x u], where x is vector of past states nad u is vector of
%   inputs
% y0 .. intial state of the system, optional
% Outputs:
% Y .. vector of outputs
% X .. vector of past outputs: X(k) = Y(k-1)
% U .. vector of inputs
%
%

%% Examples
% demo_example_gp_data.m

%% See Also
% EXAMPLE_DERIVATIVE, EXAMPLE_LM_IDENT

if(nargin==1)
    y0 = 0;
else
    y0 = varargin{1};
end

if (size(U,2)>2 || size(U,2)<1)
    error('input U must have size 1 or 2');
elseif (size(U,2)==2)
    % first coloumn = y
    % second coloumn = u
    X = U(:,1);
    U = U(:,2);
    Y = X./(1+X.^2) + U.^3;
else % size(U,2)=1
    % step1
    X = y0;
    % steps
    for ii=1:length(U)
        Y(ii) = X(ii)/(1+X(ii).^2) + U(ii).^3;
        X(ii+1) = Y(ii);
    end

    X = X(1:end-1);
    X = X'; Y = Y';
end

author: Dynamic-Systems-and-GP

repo: GPdyn


text:
function livefg(ma)
% livefg(m) where m is a initialized handle to a framegrabber shows live image
% in a window.

% Vladimir Smutny 22.10.1999

global loop

f = figure;
set(f, 'Name', 'Camera Viewer',...
   'NumberTitle', 'off',...
   'Interruptible','on',...
   'KeyPressFcn','global loop; loop=0;' );
colormap gray(256)
set(gca,'Position',[0,0,1,1]);

%image(grabfg(ma)')

mov=[];
mov.cdata = grabfg(ma)';
mov.colormap=gray(256);
movsize=size(mov.cdata);
set(gcf,'Position',[50,50,movsize(2),movsize(1)]);
loop = 1;
while (loop ~= 0),
   mov.cdata=grabfg(ma)';
   movie(mov);
   pause(0.1);
end


text:
function y = scale_cols(x, s)
% SCALE_COLS       Scale each column of a matrix.
% SCALE_COLS(x,s) returns matrix y, same size as x, such that
% y(:,i) = x(:,i)*s(i)
% It is more efficient than x*diag(s).

%y = x.*repmat(s(:)', rows(x), 1);

for i = 1:size(x,2)
    y(:,i) = x(:,i)*s(1,i);
end

%y = x.*(ones(rows(x),1)*s(:)');

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
report_this_filefun(mfilename('fullpath'));

[lat,lon] = meshgrat(tmap,tmapleg);
[X , Y]  = meshgrid(gx,gy);

ren = interp2(X,Y,re3,lon,lat);
mi = 0;

figure_w_normalized_uicontrolunits('pos',[150 100 1000 700])

hold on; axis off
axesm('MapProjection','mercator',...
    'MapLatLimit',[s4_south s3_north],'MapLonLimit',[s2_west s1_east])

meshm(ren,tmapleg,size(tmap),tmap);

daspectm('m',5);
tightmap
view([0 90])
camlight; lighting phong
% set(gca,'projection','perspective');

plotm(coastline(:,2), coastline(:,1),'w','Linewidth',2);
zdatam(handlem('allline'),10000) % keep line on surface

j = jet(64);
%j = j(64:-1:1,:);
j = [ [ 0.85 0.9 0.9] ; j];

colormap(j); brighten(0.3);
%caxis([mi 1.4])

set(gcf,'color','w')

setm(gca,'ffacecolor','w')
setm(gca,'fedgecolor','w','flinewidth',1);

setm(gca,'mlabellocation',1)
setm(gca,'meridianlabel','on')
setm(gca,'plabellocation',1)
setm(gca,'parallellabel','on')
setm(gca,'Fontcolor','w','Fontweight','bold','FontSize',12)

h5 = colorbar;
set(h5,'position',[0.82 0.35 0.01 0.3])

set(gcf,'Inverthardcopy','off');

author: CelsoReyes


function [V_RF, V_U] = gevd_algorithm(V_RF, w, H1)

global Ns  Nrf Vn;
Nt = size(V_RF,1);
theta = 1/(Vn * w * Nt);  % for simplification
i = 0;
% trigger = 1;
% newtri = 100;
for m = 1: 2  %outer iteration
    for j = 1: Nrf
        V_m = V_RF;
        V_m(:,j) = [];
        Am = eye(Ns) + theta * H1' * V_m * V_m' * H1;
        Um = theta * H1 * Am^(-2) * H1';
        Wm = 1/Nt * eye(Nt) + theta * H1 * Am^(-1) * H1';
        [V,D] = eig(Um, Wm);
        % get the largest eigenvector
        [~,max_index] = max(diag(D));
        V_RF(:,j) = exp(1i * angle(V(:,max_index)));
        %v = V_RF(:,j);
    end
end

V_U = inv(V_RF'*H1 * H1'* V_RF+  Vn * w *(V_RF)'*V_RF)*V_RF'*H1;


text:
function yc = pwl_approx_1d ( nd, xd, yd, nc, xc )

%*****************************************************************************80
%
%% PWL_APPROX_1D determines the control values for a PWL approximant.
%
%  Discussion:
%
%    The piecewise linear approximant is defined by NC control pairs
%    (XC(I),YC(I)) and approximates ND data pairs (XD(I),YD(I)).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ND, the number of data points.
%    ND must be at least 1.
%
%    Input, real XD(ND,1), the data points.
%
%    Input, real YD(ND,1), the data values.
%
%    Input, integer NC, the number of control points.
%    NC must be at least 1.
%
%    Input, real XC(NC,1), the control points.
%
%    Output, real YC(NC,1), the control values.
%

%
%  Define the NDxNC linear system that determines the control values.
%
  a = pwl_approx_1d_matrix ( nd, xd, yd, nc, xc );
%
%  Solve the system.
%
  yc = a \ yd;

  return
end


text:
function [F, M, trpy, drpy] = controller(qd, t, qn, params)
% CONTROLLER quadrotor controller
% The current states are:
% qd{qn}.pos, qd{qn}.vel, qd{qn}.euler = [roll;pitch;yaw], qd{qn}.omega
% The desired states are:
% qd{qn}.pos_des, qd{qn}.vel_des, qd{qn}.acc_des, qd{qn}.yaw_des, qd{qn}.yawdot_des
% Using these current and desired states, you have to compute the desired controls

% =================== Your code goes here ===================

% ordinary linear
% % position controller params
% Kp = ones(3,1)*30;
% Kd = ones(3,1)*10;
%
% % attitude controller params
% KpM = ones(3,1)*10000;
% KdM = ones(3,1)*500;

% position controller params
Kp = [15;15;30];
% Kd = [15;15;10];
Kd = [12;12;10];

% attitude controller params
KpM = ones(3,1)*3000;
KdM = ones(3,1)*300;

% t = qd{qn}.vel_des/norm(qd{qn}.vel_des+eps);
% n = qd{qn}.acc_des/norm(qd{qn}.acc_des+eps);
% b = cross(t,n);
% ep = ((qd{qn}.pos_des - qd{qn}.pos).*n).*n + ((qd{qn}.pos_des - qd{qn}.pos).*b).*b;
% ev = qd{qn}.vel_des - qd{qn}.vel;

acc_des = qd{qn}.acc_des + Kd.*(qd{qn}.vel_des - qd{qn}.vel) + Kp.*(qd{qn}.pos_des - qd{qn}.pos);
% acc_des = qd{qn}.acc_des + Kd.*ev + Kp.*ep

% Desired roll, pitch and yaw
phi_des = 1/params.grav * (acc_des(1)*sin(qd{qn}.yaw_des) - acc_des(2)*cos(qd{qn}.yaw_des));
theta_des = 1/params.grav * (acc_des(1)*cos(qd{qn}.yaw_des) + acc_des(2)*sin(qd{qn}.yaw_des));
psi_des = qd{qn}.yaw_des;

euler_des = [phi_des;theta_des;psi_des];
pqr_des = [0;0; qd{qn}.yawdot_des];
% Thurst
qd{qn}.acc_des(3);
F  = params.mass*(params.grav + acc_des(3));
% Moment
M =  params.I*(KdM.*(pqr_des - qd{qn}.omega) + KpM.*(euler_des - qd{qn}.euler));
% =================== Your code ends here ===================

% Output trpy and drpy as in hardware
trpy = [F, phi_des, theta_des, psi_des];
drpy = [0, 0,       0,         0];

end

author: yrlu

repo: quadrotor

save_path: github-repos/MATLAB/yrlu-quadrotor


text:
function M = slpweval(X1, X2, f, varargin)
%SLPWEVAL Perform pairwise computation
%
% $ Syntax $
%   - M = slpweval(X1, X2, f)
%   - M = slpweval(X1, X2, f, ...)
%
% $ Arguments $
%   - X1:       the matrix of vectors serving as first argument of f
%   - X2:       the matrix of vectors serving as second argument of f
%   - f:        the function maps two vectors to a single scalar value
%   - M:        the matrix of pairwise evaluaton results
%
% $ Description $
%   - M = slpweval(X1, X2, f) takes vector arguments from X1 and X2, and
%     computes the pairwise evaluation result with f. The vectors in X1
%     and X2 are stored in a column-wise manner. Suppose X1 and X2 have
%     m and n columns respectively, then the resultant matrix M would be
%     of size m x n, with the M(i, j) = f(X1(:,i), X2(:,j)).
%
%   - M = slpweval(X1, X2, f, ...) conducts the computation with extra
%     parameters to f, i.e. M(i, j) = f(X1(:,i), X2(:,j), ...).
%
% $ Remarks $
%   - The vector length of the vectors in X1 and X2 are not necessarily
%     equal. The requirment on their dimensions depends on the callback
%     function f.
%   - For efficiency, the function would invoke f to evaluate in batch.
%     Thus f should support batch-evaluation. When the input arguments
%     to f have n columns, f should return an 1 x n row vector.
%
% $ History $
%   - Created by Dahua Lin on Apr 21, 2006
%   - Modified by Dahua Lin on Sep 10, 2006
%       - make some minor changes to suppress warnings
%

%% parse and verify input arguments

if nargin < 3
    raise_lackinput('slpweval', 3);
end
[d1, n1] = size(X1);
[d2, n2] = size(X2);
slignorevars(d1, d2);

%% compute

% prepare output matrix
M = zeros(n1, n2);

if n1 > n2      % expand each column in X2 to n1 copies

    inds_e = ones(1, n1);
    for i = 1 : n2
        x2 = X2(:, i);
        X2e = x2(:, inds_e);
        M(:, i) = feval(f, X1, X2e, varargin{:})';
    end

else            % expand each column in X1 to n2 copies

    inds_e = ones(1, n2);
    for i = 1 : n1
        x1 = X1(:, i);
        X1e = x1(:, inds_e);
        M(i, :) = feval(f, X1e, X2, varargin{:});
    end

end

author: lmthang


text:
function [H0,H1,H2] = spm_bilinear(A,B,C,D,x0,N,dt)
% Return global Volterra kernels for a MIMO Bilinear system
% FORMAT [H0,H1,H2] = spm_bilinear(A,B,C,D,x0,N,dt)
% A     - (n x n)     df(x(0),0)/dx                    - n states
% B     - (n x n x m) d2f(x(0),0)/dxdu                 - m inputs
% C     - (n x m)     df(x(0),0)/du - d2f(x(0),0)/dxdu*x(0)
% D     - (n x 1)     f(x(0).0) - df(x(0),0)/dx*x(0)
% x0    - (n x 1)     x(0)
% N     - kernel depth       {intervals}
% dt    - interval           {seconds}
%
% Volterra kernels:
%
% H0    - (n)                 = h0(t)         = y(t)
% H1    - (N x n x m)         = h1i(t,s1)     = dy(t)/dui(t - s1)
% H2    - (N x N x n x m x m) = h2ij(t,s1,s2) = d2y(t)/dui(t - s1)duj(t - s2)
%
% where n = p if modes are specified
%
%--------------------------------------------------------------------------
% Returns Volterra kernels for bilinear systems of the form
%
% dx/dt = f(x,u) = A*x + B1*x*u1 + ... Bm*x*um + C1u1 + ... Cmum + D
%  y(t) = x(t)
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_bilinear.m 5219 2013-01-29 17:07:07Z spm $

% Volterra kernels for bilinear systems
%==========================================================================

% parameters
%--------------------------------------------------------------------------
n     = size(A,1);                  % state variables
m     = size(C,2);                  % inputs
A     = full(A);
B     = full(B);
C     = full(C);
D     = full(D);

% eignvector solution {to reduce M0 to leading diagonal form}
%-------------------------------------------------------------------------
M0    = [0 zeros(1,n); D A];
[U,J] = eig(M0);
V     = pinv(U);

% Lie operator {M0}
%--------------------------------------------------------------------------
M0    = sparse(J);
X0    = V*[1; x0];

% 0th order kernel
%--------------------------------------------------------------------------
H0    = ex(N*dt*M0)*X0;

% 1st order kernel
%--------------------------------------------------------------------------
if nargout > 1

    % Lie operator {M1}
    %----------------------------------------------------------------------
    for i = 1:m
        M1(:,:,i)  = V*[0 zeros(1,n); C(:,i) B(:,:,i)]*U;
    end

    % 1st order kernel
    %----------------------------------------------------------------------
    H1    = zeros(N,n + 1,m);
    for p = 1:m
    for i = 1:N
        u1         = N - i + 1;
        H1(u1,:,p) = ex(u1*dt*M0)*M1(:,:,p)*ex(-u1*dt*M0)*H0;
    end
    end
end

% 2nd order kernels
%--------------------------------------------------------------------------
if nargout > 2
    H2    = zeros(N,N,n + 1,m,m);
    for p = 1:m
    for q = 1:m
    for j = 1:N
        u2         = N - j + 1;
        u1         = N - [1:j] + 1;
        H          = ex(u2*dt*M0)*M1(:,:,q)*ex(-u2*dt*M0)*H1(u1,:,p)';
        H2(u2,u1,:,q,p) = H';
        H2(u1,u2,:,p,q) = H';
    end
    end
    end
end

% project to state space and remove kernels associated with the constant
%--------------------------------------------------------------------------
if nargout > 0
    H0    = real(U*H0);
    H0    = H0([1:n] + 1);
end
if nargout > 1
    for p = 1:m
        H1(:,:,p) = real(H1(:,:,p)*U.');
    end
    H1    = H1(:,[1:n] + 1,:);
end
if nargout > 1
    for p = 1:m
    for q = 1:m
    for j = 1:N
        H2(j,:,:,p,q) = real(squeeze(H2(j,:,:,p,q))*U.');
    end
    end
    end
    H2    = H2(:,:,[1:n] + 1,:,:);
end

return

% matrix exponential function (for diagonal matrices)
%==========================================================================
function y = ex(x)
n          = length(x);
y          = spdiags(exp(diag(x)),0,n,n);
return

author: spm


text:
function [ h1 ] = plot_sample_trajectories( xurxps )
%PLOT_SAMPLE_TRAJECTORIES Summary of this function goes here
%   Detailed explanation goes here

h1 = figure;

subplot(1,2,1);
hold on; grid on;
for n=1:size(xurxps,1)

    xurxp = xurxps{n};

    x  = xurxp(:,1:2);
    u  = xurxp(:,3);
    r  = xurxp(:,4);
    xp = xurxp(:,5:6);

    % plot start
    scatter(x(1,1),x(1,2),15,[0 0 0],'filled');
    scatter(x(end,1),x(end,2),15,[1 0 0],'filled');
    scatter(x(:,1),x(:,2),2,r,'filled');
    plot(x(:,1),x(:,2));
end
title('recorded (x,r)');
xlabel('position');
ylabel('velocity');

subplot(1,2,2);
hold on; grid on;
for n=1:size(xurxps,1)

    xurxp = xurxps{n};

    r  = xurxp(:,4);
    xp = xurxp(:,5:6);

    scatter(xp(:,1),xp(:,2),5,r,'filled');
end
title('recorded (xp,r)');
xlabel('position');
ylabel('velocity');

end

author: epfl-lasa


text:
% This script was used once to create the demo dataset
% private/demo_waveforms.m. It has no further use and is kept only as
% sample. The originla dataset was California earthquake data in SAC files
% provided by Carl Tape. Original source unknown.
%
% To be run from within private directory. At the time demo data was
% created, a directory sample_data was in private/ as well.

% clear
% clear classes
% rmgismo
% addpath('/home/field/GISMO_BRANCH/GISMO/');
% startup_GISMO

% LOAD DATA
fileList = struct2cell(dir('sample_data/9155518*S*.BHZ.*'));
filesZ = strcat('sample_data/', fileList(1,:)');
fileList = struct2cell(dir('sample_data/9155518*S*.BHN.*'));
filesN = strcat('sample_data/', fileList(1,:)');
fileList = struct2cell(dir('sample_data/9155518*S*.BHE.*'));
filesE = strcat('sample_data/', fileList(1,:)');
files = [filesZ filesN filesE];
w = loadsacfile(files);
%save junk

% TRIM SIZE
w = w([36 19 20 4 21 9 18 10 17 11 8 5 14 26],:);

% ADJUST FIELDS
backAzimuth = get(w(:,1),'BAZ');
for n = 1:numel(w)
     w(n) = addfield( w(n) , 'stationLatitude' , get(w(n), 'STLA') );
     w(n) = addfield( w(n) , 'stationLongitude' , get(w(n), 'STLO') );
     w(n) = addfield( w(n) , 'originLatitude' , get(w(n), 'EVLA') );
     w(n) = addfield( w(n) , 'originLongitude' , get(w(n), 'EVLO') );
    channel = get(w(n), 'CHANNEL');
    w(n) = set( w(n) , 'CHANNEL' , channel(1:3) );
    network = get(w(n), 'NETWORK');
    w(n) = set( w(n) , 'NETWORK' , network(1:2) );
end
w = delfield(w,'AZ');
w = delfield(w,'BAZ');
w = delfield(w,'CMPAZ');
w = delfield(w,'CMPINC');
w = delfield(w,'DIST');
w = delfield(w,'EVDP');
w = delfield(w,'GCARC');
w = delfield(w,'HISTORY');
w = delfield(w,'IFTYPE');
w = delfield(w,'IZTYPE');
w = delfield(w,'KEVNM');
w = delfield(w,'KT1');
w = delfield(w, 'KT2');
w = delfield(w,'KUSER0');
w = delfield(w,'LCALDA');
w = delfield(w,'LOVROK');
w = delfield(w,'LPSPOL');
w = delfield(w,'NZHOUR');
w = delfield(w,'NZJDAY');
w = delfield(w,'NZMIN');
w = delfield(w,'NZMSEC');
w = delfield(w,'NZSEC');
w = delfield(w,'NZYEAR');
w = delfield(w,'O');
w = delfield(w,'SCALE');
w = delfield(w,'STDP');
w = delfield(w,'STEL');
w = delfield(w,'T1');
w = delfield(w,'T2');
w = delfield(w,'STLA');
w = delfield(w,'STLO');
w = delfield(w,'EVLA');
w = delfield(w,'EVLO');

% PLOT LOCATIONS
figure;
plot( get(w(:,1),'STATIONLONGITUDE') , get(w(:,1),'STATIONLATITUDE') , 'b+');
hold on;
plot( get(w(:,1),'ORIGINLONGITUDE') , get(w(:,1),'ORIGINLATITUDE') , 'ro');
set(gca,'DataAspectRatio',[1 cosd(35) 1]);
for n=1:length(backAzimuth)
   text( get(w(n,1),'STATIONLONGITUDE') , get(w(n,1),'STATIONLATITUDE') , num2str(n) );
end

% SAVE DEMO DATA
save('demo_waveforms.mat','w','backAzimuth');

author: geoscience-community-codes

repo: GISMO


text:
%MAT_INTO_HDR  The old versions of SPM (any version before SPM5) store
%	an affine matrix of the SPM Reoriented image into a matlab file
%	(.mat extension). The file name of this SPM matlab file is the
%	same as the SPM Reoriented image file (.img/.hdr extension).
%
%	This program will convert the ANALYZE 7.5 SPM Reoriented image
%	file into NIfTI format, and integrate the affine matrix in the
%	SPM matlab file into its header file (.hdr extension).
%
%	WARNING: Before you run this program, please save the header
%	file (.hdr extension) into another file name or into another
%	folder location, because all header files (.hdr extension)
%	will be overwritten after they are converted into NIfTI
%	format.
%
%  Usage: mat_into_hdr(filename);
%
%  filename:	file name(s) with .hdr or .mat file extension, like:
%		'*.hdr', or '*.mat', or a single .hdr or .mat file.
%	e.g.	mat_into_hdr('T1.hdr')
%		mat_into_hdr('*.mat')
%

%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
%-------------------------------------------------------------------------
function mat_into_hdr(files)

   pn = fileparts(files);
   file_lst = dir(files);
   file_lst = {file_lst.name};
   file1 = file_lst{1};
   [p n e]= fileparts(file1);

   for i=1:length(file_lst)
      [p n e]= fileparts(file_lst{i});
      disp(['working on file ', num2str(i) ,' of ', num2str(length(file_lst)), ': ', n,e]);
      process=1;

      if isequal(e,'.hdr')
         mat=fullfile(pn, [n,'.mat']);
         hdr=fullfile(pn, file_lst{i});

         if ~exist(mat,'file')
            warning(['Cannot find file "',mat  , '". File "', n, e, '" will not be processed.']);
            process=0;
         end
      elseif isequal(e,'.mat')
         hdr=fullfile(pn, [n,'.hdr']);
         mat=fullfile(pn, file_lst{i});

         if ~exist(hdr,'file')
            warning(['Can not find file "',hdr  , '". File "', n, e, '" will not be processed.']);
            process=0;
         end
      else
         warning(['Input file must have .mat or .hdr extension. File "', n, e, '" will not be processed.']);
         process=0;
      end

      if process
         load(mat);
         R=M(1:3,1:3);
         T=M(1:3,4);
         T=R*ones(3,1)+T;
         M(1:3,4)=T;

         [h filetype fileprefix machine]=load_nii_hdr(hdr);
         h.hist.qform_code=0;
         h.hist.sform_code=1;
         h.hist.srow_x=M(1,:);
         h.hist.srow_y=M(2,:);
         h.hist.srow_z=M(3,:);
         h.hist.magic='ni1';

         fid = fopen(hdr,'w',machine);
         save_nii_hdr(h,fid);
         fclose(fid);
      end
   end

   return;				% mat_into_hdr

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [C, sigma] = dataset3Params(X, y, Xval, yval)
%EX6PARAMS returns your choice of C and sigma for Part 3 of the exercise
%where you select the optimal (C, sigma) learning parameters to use for SVM
%with RBF kernel
%   [C, sigma] = EX6PARAMS(X, y, Xval, yval) returns your choice of C and
%   sigma. You should complete this function to return the optimal C and
%   sigma based on a cross-validation set.
%

% You need to return the following variables correctly.
C = 1;
sigma = 0.3;

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the optimal C and sigma
%               learning parameters found using the cross validation set.
%               You can use svmPredict to predict the labels on the cross
%               validation set. For example,
%                   predictions = svmPredict(model, Xval);
%               will return the predictions on the cross validation set.
%
%  Note: You can compute the prediction error using
%        mean(double(predictions ~= yval))
%

% C = best C
% sigma = best Sigma
values = [0.01 0.03 0.1 0.3 1, 3, 10 30]';

% all predictions
predictions = zeros(size(values));

for i=1:size(values,1), % C
    for j=1:size(values,1), % sigma

        model = svmTrain(X, y, values(i), @(x1, x2) gaussianKernel(x1, x2, values(j)));
        predict = svmPredict(model, Xval);
        predictions(i,j) = mean(double(predict ~= yval));

    end
end

[colmin, rowindex] = min(predictions);
[minerror, index] = min(colmin);

C = values(rowindex(index))
sigma = values(index)

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

end

author: rmarquis

repo: coursera-machinelearning


text:
function bnet_miss = gener_MAR_net(bnet_orig, base_proba)
% function bnet_miss = gener_MAR_net(bnet_orig, base_proba)
%
%   bnet_orig : a bnet
%   base_proba :  a probability for value to be missing
%
%   bnet_miss : a bnet that could be used in gener_data_from_bnet_miss function
%               to generate incomplete MAR dataset
%
% Francois.Olivier.C.H@gmail.com

%%%%%%%%%%%% INIT
if nargin<2, error('Not enougth arguments'); end

% création du réseau
dag = bnet_orig.dag;
N = size(dag,2);
ns = bnet_orig.node_sizes;

  ns_miss = zeros(1,3*N);
  ns_miss(1:N) = ns;
  ns_miss(N+1:2*N) = 2*ones(1,N); % 1= node i-N present, 2= node i-N missing
  ns_miss(2*N+1:3*N) = ns+1; % 1:ns, absent

  dag_miss = zeros(3*N,3*N);
  dag_miss(1:N,1:N) = dag;
%  dag_miss(2*N+1:3*N,N+1:2*N)=mk_rnd_dag(N,N-ceil(rand*N/2));
  lim = 1+(rand>.4)+(rand>.65)+(rand>.9);
  dag_miss(2*N+1:3*N,N+1:2*N)=mk_rnd_dag(N,lim);

  for i=1:N, dag_miss(i,2*N+i)=1; dag_miss(N+i,2*N+i)=1; dag_miss(2*N+i,i)=0; dag_miss(2*N+i,N+i)=0; end

  bnet_miss = mk_bnet(dag_miss, ns_miss);
  CPT = CPT_from_bnet(bnet_orig, 0);
  for i=1:N
    bnet_miss.CPD{i} = tabular_CPD (bnet_miss, i, CPT{i});
  end

  % CPD of nodes M
  for i=1:N,
     if find(bnet_miss.order==i)<find(bnet_miss.order==N+i),
        CPT_M=[];
        for j=1:ns_miss(i), for l=1:ns_miss(N+i), for k=1:ns_miss(2*N+i),
          CPT_M=[CPT_M (((j==k)&(l==1))|((k==ns_miss(2*N+i))&(l==2)))];
        end, end, end
     else
        CPT_M=[];
        for k=1:ns_miss(2*N+i), for l=1:ns_miss(N+i), for j=1:ns_miss(i),
          CPT_M=[CPT_M (((j==k)&(l==1))|((k==ns_miss(2*N+i))&(l==2)))];
        end, end, end
     end
     bnet_miss.CPD{2*N+i} = tabular_CPD (bnet_miss, 2*N+i, CPT_M);
  end

%%%%%%%%%%%% Base probability of missing value
p = base_proba;
 for i=1:N
  fam = find(dag_miss(:,N+i)==1)';
  semisize = prod(ns_miss(fam)); % as node N+i is binary to say i is present or missing
  CPT = zeros(1,2*semisize);
  CPT(1:semisize) = 1-p;
  CPT(semisize+1:2*semisize) = p;
  bnet_miss.CPD{N+i} = tabular_CPD (bnet_miss, N+i, CPT);
 end

   BETA = gener_discrete_dist(N, base_proba);

   order=[];
   missdagtmp = bnet_miss.dag; %(N+1:2*N,N+1:2*N);
   unprocessed = 1:N;
   while ~isempty(unprocessed)
     npar=[];
     for i=N+1:2*N, npar(end+1)=length(parents(missdagtmp,i));end,  % to be verifie from here
     [npar, ord] = sort(npar);
     while ~ismember(ord(1),unprocessed)
       ord=ord(2:end);
     end
     order = [order, ord(1)];
     missdagtmp(ord(1),:)=0;
     unprocessed = mysetdiff(unprocessed,ord(1));
   end

%%%%%%%%%%%% Update CPT with MCAR process
for i=1:length(BETA)
  fam_miss = find(dag_miss(:,N+order(i))==1)';
  p=BETA(i);
  semisize = prod(ns_miss(fam_miss));

  if isempty(fam_miss),
        CPT = zeros(1,2*semisize);
        CPT(1:semisize) = 1-p;
        CPT(semisize+1:2*semisize) = p;
        bnet_miss.CPD{N+order(i)} = tabular_CPD (bnet_miss, N+order(i), CPT);
  else

        MUi1k = gener_discrete_dist(semisize, p);
        CPT = zeros(1,2*semisize);
        CPT(1:semisize) = 1-MUi1k;
        CPT(semisize+1:2*semisize) = MUi1k;
        bnet_miss.CPD{N+order(i)} = tabular_CPD (bnet_miss, N+order(i), CPT);

  end

end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function pass = test_dlt(pref)

if ( nargin == 0 )
    pref = chebfunpref();
end
seedRNG(42)

n = 10;
x = legpts(n);
r = rand(n,1)./(1:n)';
rr = [r, r];
tol = 100*n*eps;
P = legpoly(0:n-1);

N = 5001;
R = rand(N,1)./(1:N)';
RR = [R, R];
Tol = 100*N*eps;

%% Basic test (small)
err = norm(P(x)*r - chebfun.dlt(r), inf);
pass(1) = err < tol;

%% Vector test (small)
err = norm(P(x)*rr - chebfun.dlt(rr), inf);
pass(2) = err < tol;

%% Test direct
err = norm(P(x)*r - dlt_direct(r), inf);
pass(3) = err < tol;

%% Basic test (large)
err = norm(dlt_direct(R) - chebfun.dlt(R), inf);
pass(4) = err < Tol;

%% Vector test (large)
err = norm(dlt_direct(RR) - chebfun.dlt(RR), inf);
pass(5) = err < Tol;

end

function v = dlt_direct(c)
N = size(c, 1);
if ( N == 0 ), v = 1 + 0*c; return, end
x = legpts(N);
v = repmat(c(1,:), length(x), 1) + bsxfun(@times, c(2,:), x);
Pm1 = 1 + 0*x; P = x;   % P_0 and P_1.
for n = 1:(N-2) % Recurrence relation.
    Pp1 = (2-1/(n+1))*(P.*x) - (1-1/(n+1))*Pm1;
    Pm1 = P;    P = Pp1;
    v = v + bsxfun(@times, c(n+2,:), P);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function modelOut = addExoMetToEFBA(model,exoMet,param)
% generates H and h to add a min (alpha/2)(v-h)'*H*(v-h) term to an EFBA problem
%
% USAGE:
%   modelOut = addExoMetToEFBA(model,exoMet,param)
%
% INPUTS:
%  model.S:
%  model.rxns:
%
%  exoMet.rxns:
%  exoMet.mean:
%  exoMet.SD:
%
% OPTIONAL INPUT
%  param.printLevel:
%  param.alpha:              alpha in (alpha/2)(v-h)'*H*(v-h), default = 10000;
%  param.metabolomicWeights: String indicating the type of weights to be applied to penalise the difference
%                            between of predicted and experimentally measured fluxes by, where
%                            'SD'   weights = 1/(1+exoMet.SD^2)
%                            'mean' weights = 1/(1+exoMet.mean^2)  (Default)
%                            'RSD'  weights = 1/((exoMet.SD./exoMet.mean)^2)
%  param.relaxBounds: True to relax bounds on reaction whose fluxes are to be fitted to exometabolomic data
%
% OUTPUTS:
%  modelOut:
%
% EXAMPLE:
%
% NOTE:
%
% Author(s):

if ~exist('param','var')
    param=struct;
end
if ~isfield(param,'printLevel')
    param.printLevel=0;
end
if ~isfield(param,'alpha')
    param.alpha=10000;
end
if ~isfield(param,'relaxBounds')
    param.relaxBounds=1;
end
if ~isfield(param,'metabolomicWeights')
    param.metabolomicWeights = 'mean';
end

[bool, locb] = ismember(exoMet.rxns, model.rxns);
if any(locb == 0)
    if printLevel > 0
        fprintf('%s\n','The following exometabolomic exchange reactions were not found in the model:')
        disp(exoMet.rxns(locb == 0))
    end
end

if length(unique(exoMet.rxns)) ~= length(exoMet.rxns) && printLevel > 0
    disp('There are duplicate rxnID entries in the metabolomic data! Only using data corresponding to first occurance')
end

% Assume mean model flux is equal to the mean experimental reaction flux.
[~,nRxn] = size(model.S);
vExp = NaN * ones(nRxn,1);
vExp(locb(bool)) = exoMet.mean(bool);

vSD = NaN * ones(nRxn,1);
vSD(locb(bool)) = exoMet.SD(bool);

% Set the weight on the Euclidean distance of the predicted steady state
% flux from the experimental steady state flux. Uniform at first.
weightExp = ones(nRxn,1);
% The weight penalty on relaxation from the experimental flux vector should be greater than that
% on the lower and upper bounds, if the experimental data is considered more reliable
% than the lower and upper bounds of the model.
% Assumes the lower and upper bound of standard deviation of
% experimental reaction flux are separated by two standard
% deviations.

% Penalise the relaxation from the experimental flux value by 1/(1+weights^2)
switch param.metabolomicWeights
    case 'SD'
        weightExp(locb(bool)) = 1 ./ (1 +  (vSD(locb(bool))).^2);
    case 'mean'
        weightExp(locb(bool)) = 1 ./ (1 + (vExp(locb(bool))).^2);
    case 'RSD'
        weightExp(locb(bool)) = 1 ./ ((vSD(locb(bool))./vExp(locb(bool))).^2);
    otherwise
        weightExp(locb(bool)) = 2;
end
% Weights are ignored on the reactions without experimental data, i.e. vExp(n) == NaN.
weightExp(~isfinite(vExp)) = 0;

%add
model.H=diag(sparse(weightExp*param.alpha));
model.h=vExp;
if param.printLevel>1
    figure;
    histogram(weightExp)
    xlabel('weights on the diagonal of model.H')
    ylabel('#reactions')
end

if param.relaxBounds
    bool=ismember(model.rxns,exoMet.rxns);
    model.lb(bool) = model.lb_preconstrainRxns(bool);
    model.ub(bool) = model.ub_preconstrainRxns(bool);
end

modelOut = model;
end


text:
%% ellipseCoord
% Below is a demonstration of the features of the |ellipseCoord| function

%%
clear; close all; clc;

%% Syntax
% |[V]=ellipseCoord(A,t);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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 3 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, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function fcnn = fcnnupdateweights(fcnn)
%FCNNUPDATEWEIGHTS Update the weights of fully-connected neural net.
%   FCNN = FCNUPDATEWEIGHTS(FCNN]) updates the weights of the
%   fully-connected net, FCNN. The following code implements Step 4 in
%   Table 14.9.
%
%   Type
%
%    >> help fcninfo
%
%   at the prompt for detailed explanations of the components of the
%   fully-connected neural net.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Number of layers.
L = numel([fcnn.NumNodes]);

% Learning rate constant.
alpha = fcnn(L).Alpha;

% Update the weights.
for k = 2:L
   fcnn(k).Weights = fcnn(k).Weights - alpha*(fcnn(k).D)*(fcnn(k - 1).A');
   fcnn(k).Biases = fcnn(k).Biases - alpha*sum(fcnn(k).D,2);
end

author: dipum

repo: dipum-toolbox

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


text:
function cleanInplane()
%
% function cleanInplane()
%
% Deletes:
%   Gray/coords.mat
%   Gray/dataType/corAnal.mat (for each dataType)
%   Gray/dataType/*.mat (all other parameter maps)
%
% If you change this function make parallel changes in:
%     cleanFlat, cleanDataType
%
% djh, 2/2001

global HOMEDIR
global dataTYPES

for typeNum = 1:length(dataTYPES)
    dataTypeName = dataTYPES(typeNum).name;
    % Delete corAnal and parameter map files
    datadir = fullfile(HOMEDIR,'Inplane',dataTypeName);
    delete(fullfile(datadir,'*.mat'));
    % Delete tSeries (if there are any)
    [nscans,scanDirList] = countDirs(fullfile(datadir,'TSeries','Scan*'));
    for s=1:nscans
        delete(fullfile(datadir,'TSeries',scanDirList{s},'*.mat'));
    end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function value = p08_f ( dim_num, point_num, x )

%*****************************************************************************80
%
%% P08_F evaluates the integrand for problem 08.
%
%  Dimension:
%
%    DIM_NUM arbitrary.
%
%  Region:
%
%    0 <= X(1:DIM_NUM) <= 1
%
%  Integrand:
%
%    ( sin ( (pi/4) * sum ( x(1:dim_num) ) ) )**2
%
%  Exact Integral:
%
%    1/2 - sqrt ( 2**(3*DIM_NUM) ) * cos ( DIM_NUM * pi / 4 ) ) / ( 2 * pi**DIM_NUM )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Richard Crandall,
%    Projects in Scientific Computing,
%    Springer, 2000,
%    ISBN: 0387950095,
%    LC: Q183.9.C733.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, real X(DIM_NUM,POINT_NUM), the evaluation points.
%
%    Output, real VALUE(POINT_NUM), the integrand values.
%
  value(1:point_num) = 0.0;

  for point = 1 : point_num
    value(point) = ( sin ( pi * sum ( x(1:dim_num,point) ) / 4.0 ) )^2;
  end

  p08_i4 ( 'I', '#', point_num );

  return
end

author: johannesgerer


text:
% IMAGE_SCAT_LAYER engine function of IMAGE_SCAT
%
% Usage
%     big_img = image_scat_layer(Scatt, renorm, dsp_legend)
%
% Input
%	Scatt (struct): a layer of  scattering (either U or S)
%	renorm (boolean): if 1 renormalize each path by its max. Default
%	value is set to 1.
%	dsp_legend (boolean): if set to 1, display legend. Default value is 1
%
% Output
%	big_img (numerical): a large image where all paths are concatenated

function big_img = image_scat_layer(Scatt, renorm, dsp_legend)
	if (numel(size(Scatt.signal{1})) == 3)
		p2 = 1;
		for p = 1:numel(Scatt.signal)
			for theta = 1:size(Scatt.signal{1},3)
				ScattBis.signal{p2} = Scatt.signal{p}(:,:,theta);
				ScattBis.meta.j(:,p2) = Scatt.meta.j(:,p);
				ScattBis.meta.theta2(:,p2) = Scatt.meta.theta2(:,p);
				ScattBis.meta.theta(:,p2) = theta;
				ScattBis.meta.k2(:,p2) = Scatt.meta.k2(:,p);
				p2 = p2 + 1;
			end
		end

		if ~exist('renorm','var');
			renorm  = 1;
		end
		if ~exist('dsp_legend','var');
			dsp_legend  = 1;
		end
		big_img = image_scat_layer(ScattBis, renorm, dsp_legend);
		return
	end
	if ~exist('renorm','var');
		renorm  = 1;
	end
	if ~exist('dsp_legend','var');
		dsp_legend  = 1;
	end

	numel_per_width= zeros(10000,1);
	numel_per_height = zeros(10000,1);
	if (~isfield(Scatt,'signal'))
		Scatt2 = Scatt;
		clear Scatt; % avoid warning
		Scatt.signal = Scatt2;
	end
	if (size(Scatt.signal{1},3)==3)
		for rgb = 1:3
			for p=1:numel(Scatt.signal)
				Scattrgb{rgb}.signal{p}=squeeze(Scatt.signal{p}(:,:,rgb));
			end
			big_img(:,:,rgb) = display_scatt_2d_all_layer(Scattrgb{rgb});
		end
	else
		for p = 1:numel(Scatt.signal)
			curr_sz = size(Scatt.signal{p});
			numel_per_width(curr_sz(2)) = numel_per_width(curr_sz(2)) + 1;
			numel_per_height(curr_sz(1)) = numel_per_height(curr_sz(1)) + 1;
		end

		% block style
		widths = find(numel_per_width);

		total_width = ceil(sqrt(numel_per_width(max(widths))))* max(widths);
		curr_y = 0;

		% first pass to compute size and allocate

		for j = numel(widths):-1:1
			w= widths(j);
			curr_x = 0;

			for p = 1:numel(Scatt.signal)
				curr_sig = Scatt.signal{p};
				if (size(curr_sig,2) == w)
					h = size(curr_sig,1);
					if (curr_x >= total_width)
						curr_x = 0;
						curr_y = curr_y + size(curr_sig,1);
					end
					curr_x = curr_x + w;
				end
			end
			curr_y = curr_y + h;

		end
		big_img = zeros(curr_y,total_width);
		% second pass to put in the matrix what has to be put in
		curr_y = 0;
		for j = numel(widths):-1:1
			w= widths(j);
			curr_x = 0;

			for p = 1:numel(Scatt.signal)
				curr_sig = Scatt.signal{p};
				if (size(curr_sig,2) == w)
					h = size(curr_sig,1);
					if (curr_x >= total_width)
						curr_x = 0;
						curr_y = curr_y + size(curr_sig,1);
					end
					K = max(max(abs(Scatt.signal{p})));
					if (renorm==0)
						ds = size(Scatt.signal{1},1)/size(Scatt.signal{p},1);
						K = ds;
					end
					if (renorm==3)
						K=1 ;% preserve_l2_norm is off;
					end
					if (renorm == 2) % input is in log2 : substract the constant !
						ds = size(Scatt.signal{1},1)/size(Scatt.signal{p},1);
						K = log2(ds);
						big_img( (1:size(Scatt.signal{p},1)) + curr_y , (1:size(Scatt.signal{p},2)) + curr_x ) = Scatt.signal{p} -K;
					else
						big_img( (1:size(Scatt.signal{p},1)) + curr_y , (1:size(Scatt.signal{p},2)) + curr_x ) = Scatt.signal{p}/K;
					end
					curr_x = curr_x + w;
				end
			end
			curr_y = curr_y + h;

		end

		% compute min and max for future renormalization of legend
		img_min = min(big_img(:));
		img_max = max(big_img(:));
		%big_img =  (big_img-m)/(M-m);

		% last pass for legend
		curr_y = 0;
		for j = numel(widths):-1:1
			w= widths(j);
			curr_x = 0;

			for p = 1:numel(Scatt.signal)
				curr_sig = Scatt.signal{p};
				if (size(curr_sig,2) == w)
					h = size(curr_sig,1);
					if (curr_x >= total_width)
						curr_x = 0;
						curr_y = curr_y + size(curr_sig,1);
					end

					if (isstruct(Scatt) && isfield(Scatt,'meta'))
						str_legend = meta2str(Scatt.meta,p);
						if (numel(str_legend)>0)
							str_legend_split = textscan(str_legend,'%s','delimiter',' ');
						else
							str_legend_split{1} ={};
						end
					else % no meta to display
						str_legend_split{1}={};
					end
					if (dsp_legend ==0)
						str_legend_split{1}={};
					end
					str_legend_split = str_legend_split{1};
					for jl = 1:numel(str_legend_split)

						str = str_legend_split{jl};

						imstr = str2imfast(str);
						% renormalize imstr
						imstr = img_min + (img_max-img_min)*imstr;
						[n,m] = size(imstr);

						y_min =   n*(jl-1) + 1 +  curr_y;
						y_max =   n*(jl-1) + n + curr_y;
						x_min =   curr_x + 1;
						x_max =   curr_x + m;

						y_min = min(y_min,size(big_img,1));
						x_min = min(x_min,size(big_img,2));
						y_max = min(y_max,size(big_img,1));
						x_max = min(x_max,size(big_img,2));
						ny = y_max - y_min +1;
						nx = x_max - x_min +1;
						big_img(y_min:y_max,x_min:x_max) = max(imstr(1:ny,1:nx),0.5* big_img(y_min:y_max,x_min:x_max));
					end

					curr_x = curr_x + w;
				end
			end
			curr_y = curr_y + h;

		end

	end
	imagesc(big_img);
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function files = getFilesInDir(varargin)
% List all files in the supplied (git tracked) Directory with their absolute path name
% based on the git ls-file command. If the directory is not git controlled, the
% type is assumed to be 'all' and all files (except for .git files
% will be returned).
%
% USAGE:
%    files = getFilesInDir(varargin)
%
% OPTIONAL INPUTS:
%    varargin:     Options as 'ParameterName' value pairs with the following options:
%
%                   - `dirToList`: the directory to list the files for (default: the current working directory)
%                   - `type`: Git type of files to return
%
%                      - `tracked`: Only tracked files
%                      - `ignored`: Only `git` ignored files excluding files that are tracked but ignored (specified in .gitignore). If the folder is not controlled by git, `all` will be used.
%                      - `untracked`: anything that is not ignored and not tracked (new files)
%                      - `all`: all files except for the git specific files (e.g. .git, .gitignore etc).
%                      - `ignoredByCOBRA`: use the COBRA Toolbox .gitignore file. Only return those files that match patterns specified there. Slower than `ignored` since all files have to be manually checked against the expressions specified in the .gitignore file, but available on all folders (default: `all`)
%                   - `restrictToPattern` - give a regexp pattern to filter the files, this option is ignored if empty (default: '', i.e. ignored)
%                   - `checkSubFolders - check the subfolders of the current directory (default: `true`)
%                   - `printLevel`
%
%                      - `0`: No print out (default)
%                      - `1`: print Information
%
%
% OUTPUTS:
%    files:         A Cell Array of files with absolute file pathes.
%                   present in this folder matching the options choosen.
%
% EXAMPLES:
%
%    % get all m files in the source folder:
%    files = getFilesInDir('dirToList', [CBTDIR filesep 'src'], 'restrictToPattern', '\.m$');
%
%    % Get the git tracked files in the test Directory.
%    files = getFilesInDir('dirToList', [CBTDIR filesep 'test'], 'type', 'tracked');
%
%    % get all git tracked files  which start with "MyFile" in the current directory
%    files = getFilesInDir('type', 'tracked', 'restrictToPattern', '^MyFile');
%
%    % get only the gitIgnored files in the current folder
%    files = getFilesInDir('type', 'ignored');

persistent COBRAIgnored

if isempty(COBRAIgnored)
    COBRAIgnored = regexptranslate('wildcard',getIgnoredFiles());
end

gitFileTypes = {'tracked','all','untracked','ignored','ignoredByCOBRA'};
parser = inputParser();
parser.addParamValue('dirToList',pwd,@(x) exist(x,'file') == 7);
parser.addParamValue('type','all',@(x) ischar(x) && any(strcmpi(x,gitFileTypes)));
parser.addParamValue('restrictToPattern','',@(x) ischar(x));
parser.addParamValue('checkSubFolders','',@(x) islogical(x) || (isnumeric(x) && x == 0 || x == 1));
parser.addParamValue('printLevel',0,@(x) isnumeric(x));

parser.parse(varargin{:});

printLevel = parser.Results.printLevel;

%get the absolute path of the files that are listed
currentDir = cd(parser.Results.dirToList);

%get the absolute path to the folder to be listed.
absPath = pwd;

%get the type of files to obtain.
selectedType = lower(parser.Results.type);

%test, whether the folder is git controlled
[gitStatus,~] = system('git status');
if gitStatus ~= 0 && ~strcmpi(selectedType,'ignoredbycobra')
    selectedType = 'all';
end
if gitStatus == 0 && strcmpi(selectedType,'ignoredbycobra')
    [~,repos] = system('git remote -v');
    if any(cellfun(@(x) ~isempty(strfind(x,'cobratoolbox.git')),strsplit(repos,'\n')))
        %So, we are on a COBRA repo. lets just use the ignored option
        selectedType = 'ignored';
    end
end

if gitStatus == 0 && printLevel > 0 && strcmp(selectedType,'ignored')
    [~,folder] = system('git rev-parse --show-toplevel');
    fprintf('Using the .gitignore files as specified in the repository under:\n%s\n',folder);
end

switch selectedType
    case 'all'
        if gitStatus == 0
            [status, trackedfiles] = system('git ls-files');
            [status, untrackedfiles] = system('git ls-files -o');
            trackedfiles = strsplit(strtrim(trackedfiles), '\n');
            untrackedfiles = strsplit(strtrim(untrackedfiles), '\n');
            files = [trackedfiles,untrackedfiles];
        else
           if parser.Results.checkSubFolders
               rdircall = ['**' filesep '*'];
           else
               rdircall = ['*'];
           end
           files = rdir(rdircall);
           files = {files.name}'; %Need to transpose to give consistent results.
        end

    case 'tracked'
        [status, files] = system('git ls-files');
        files = strsplit(strtrim(files), '\n');
    case 'untracked'
        %Files, which are not tracked but are not ignored.
        [status, files] = system('git ls-files -o --exclude-standard');
        files = strsplit(strtrim(files), '\n');
    case 'ignored'
        [~, files] = system('git ls-files -o -i --exclude-standard');
        files = strsplit(strtrim(files), '\n');
    case 'ignoredbycobra'
        if gitStatus == 0
            [status, trackedfiles] = system('git ls-files');
            [status, untrackedfiles] = system('git ls-files -o');
            trackedfiles = strsplit(strtrim(trackedfiles), '\n');
            untrackedfiles = strsplit(strtrim(untrackedfiles), '\n');
            files = [trackedfiles,untrackedfiles];
        else
           files = rdir(['**' filesep '*']);
           files = {files.name}'; %Need to transpose to give consistent results.
        end
        matching = false(size(files));
        for i = 1:numel(COBRAIgnored)
            matching = matching | ~cellfun(@(x) isempty(regexp(x,COBRAIgnored{i},'ONCE')),files);
        end
        files = files(matching);

end
if ~parser.Results.checkSubFolders
    files = files(cellfun(@(x) isempty(regexp(x,regexptranslate('escape',filesep))),files));
end

[mlt,nlt]=size(files);
if mlt > nlt || mlt==0
    files=files';
end
if isempty(files)
    files = strcat(repmat({absPath}, 1, length(files)), repmat({filesep}, 1, length(files)));
else
    files = strcat(repmat({absPath}, 1, length(files)), repmat({filesep}, 1, length(files)),files);
end

%Filter according to the restriction pattern.
if ~isempty(parser.Results.restrictToPattern)
    files = files(cellfun(@(x) ~isempty(regexp(x,parser.Results.restrictToPattern)), files));
end

cd(currentDir);
end

author: opencobra


text:
function nseq = bal_seq_enum ( n )

%*****************************************************************************80
%
%% BAL_SEQ_ENUM enumerates the balanced sequences.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the number of 0's (and 1's) in the sequence.
%    N must be nonnegative.
%
%    Output, integer NSEQ, the number of balanced sequences.
%
  nseq = i4_choose ( 2*n, n ) / ( n + 1 );

  return
end

repo: jburkardt-m


text:
function gee_its_simple_check (A, name, b)
%GEE_ITS_SIMPLE_CHECK private function to check input arguments
% Ensures the matrix A is square, and that the right-hand-side b has the same
% number of rows as A (if present).  All matrices must be 2D, as well.
%
% Example:
%   gee_its_simple_check (A, 'A', b)
%
% See also: gee_its_simple

% Copyright 2007, Timothy A. Davis.
% http://www.cise.ufl.edu/research/sparse

[m n] = size (A) ;
if (m ~= n)
    error ('%s must be square', name) ;
end

if (ndims (A) ~= 2)
    error ('%s must be a 2D matrix', name) ;
end

if (nargin > 2)
    if (m ~= size (b,1))
        error ('%s and b must have the same number of rows', name) ;
    end
    if (ndims (b) ~= 2)
        error ('b must be a 2D matrix') ;
    end
end

author: fieldtrip

repo: fieldtrip


text:
function varargout = iszero(varargin)
%ISZERO   Check if a DISKFUN is identically zero on its domain.
%   OUT = ISZERO( F ) return 1 if the DISKFUN is exactly the zero function, and
%   0 otherwise.

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

[varargout{1:nargout}] = iszero@separableApprox(varargin{:});

end


text:
function snirfSegment()
% Description: This function divides the nirs data into segments predefined by
% user input. User should input the time period of the segments of the nirs
% data as eg [0 200; 400:900; 950:1000] in *sec*
% Meryem A Yucel, Nov 2016
%%%%%%%%%%%%%%%
% Code refined for Homer3 and fNIRS files

% Homer3 Note
% Segment tool does not work for a flat file structure

global maingui;
if isempty(maingui)
    [files, pathnm] = uigetfile( '*.snirf', 'Pick the .snirf file', 'multiselect','on');
else
    fpath = [maingui.dataTree.currElem.path maingui.dataTree.currElem.name];
    [files, pathnm] = uigetfile( '*.snirf', 'Pick the .snirf file', fpath, 'multiselect','on');
end

if files==0
    return
end

if ~iscell(files)
    if files==0
        return;
    end
end

[~,name,~] = fileparts(files);

fsn = inputdlg( 'Select the time range (in seconds) for the segment of data to save in a separate file. You can enter multiple time ranges, separated by a '';'', to save different segments to different files. For example, [0 100; 300 400] would give you two segments of the original file from 0 to 100 and 300 to 400', 'Segment SNIRF file', 1 );
if isempty(fsn)
    return
end
fsn = str2num(fsn{1});

wd = cd;
cd(pathnm)

if ~iscell(files)
    foo{1} = files;
    files = foo;
end

for iFile = 1:length(files)
    snirfData = SnirfClass(files{iFile});
    fs = 1/(snirfData.data.time(2) - snirfData.data.time(1));
    maxT = max(snirfData.data.time);

    if fsn(1,1) == 0 % if time starts from 0, take the first data sample!
        fsn(1,1) = 1/fs;
    end

    for P = 1:size(fsn,1) % loop over different parts
        snirfData = SnirfClass(files{iFile});
        if fsn(P,1)> maxT || fsn(P,2) > maxT
            errordlg('Time period (in sec) exceeds the maximum time. Please re-try.','Retry');
            return
        elseif fsn(P,1) < 0 || fsn(P,2) < 0
            errordlg('Time period (in sec) should consist of positive numbers. Please re-try.','Retry');
            return
        end

        snirfData.data.dataTimeSeries = snirfData.data.dataTimeSeries(round(fsn(P,1)*fs):round(fsn(P,2)*fs),:);
        snirfData.data.time = snirfData.data.time(round(fsn(P,1)*fs):round(fsn(P,2)*fs));
        for iStim = 1:length(snirfData.stim)
            if ~isempty(snirfData.stim(iStim).data)
                snirfData.stim(iStim).data = snirfData.stim(iStim).data(snirfData.stim(iStim).data(:,1)>fsn(P,1) & snirfData.stim(iStim).data(:,1)<fsn(P,2),:);
            end
        end
        for iAux = 1:length(snirfData.aux)
            snirfData.aux(iAux).dataTimeSeries = snirfData.aux(iAux).dataTimeSeries(round(fsn(P,1)*fs):round(fsn(P,2)*fs));
            snirfData.aux(iAux).time = snirfData.aux(iAux).time(round(fsn(P,1)*fs):round(fsn(P,2)*fs));
        end

        snirfName = sprintf([name '_seg_' num2str(P) '.snirf']);
        snirfData.Save(snirfName);
        msgbox(['File created with name' snirfName], 'Notification');
    end
end

new_name = [files{1} '.orig'];
movefile (files{1}, new_name);
% if isempty(maingui)
%     %Nothing
% else
%     for iFile = 1:length(files)
%
%         maingui.dataTree.
%     end
% end

cd(wd);

author: BUNPC


text:
function result = isVertexInRadius(vertex, origin, radius)

vt = vertex - origin;

result = ( ((vt(1)^2 + vt(2)^2) <= (radius^2)) & ...
            ((vt(2)^2 + vt(3)^2) <= (radius^2)) & ...
            ((vt(1)^2 + vt(3)^2) <= (radius^2)) );


text:
function C = train_RSLDAshrink(xTr, yTr, sublab, varargin)
% TRAIN_RSLDASHRINK - Relevance Subclass LDA with mean shrinkage to subclasses and automatic covariance shrinkage
%
%Synopsis:
%   C = train_RLSDAshrink(XTR, YTR, SUBLABELS)
%   C = train_RLSDAshrink(XTR, YTR, SUBLABELS, OPTS)
%
%Arguments:
%   XTR: DOUBLE [NxM] - Data matrix, with N feature dimensions, and M training points/examples.
%   YTR: INT [CxM] - Class membership labels of points in X_TR. C by M matrix of training
%                     labels, with C representing the number of classes and M the number of training examples/points.
%                     Y_TR(i,j)==1 if the point j belongs to class i.
%   SUBLABELS[Mx1] - Subclass membership for each data point
%   OPT: PROPLIST - Structure or property/value list of optional
%                   properties.
%     'Whitening' - BOOL (default 1): If true, mean shrinkage is performed
%       on whitened data. This is important for ERP data, as it has a
%       highly skewed eigenspectrum and the mean shrinkage is dominated by
%       directions of high variance.
%     'ReturnRegularizationProfile' - BOOL (default 1): If true, returns regularization
%       profile as matrix MxM for both classes, with M being the number of
%       unique subclasses.
%
%Returns:
%   C: STRUCT - Trained classifier structure, with the subclassifier
%     hyperplanes given by fields C.subC{k}.w and C.subC{k}.b for each
%     subclass k
%
% C includes the fields:
%    'sublab_unique' :      unique subclass labels
%    'pattern' :            pattern (i.e. mu2-mu1) of each subclassifier
%    'subC' :               LDA classifiers (w and b) for each subclass
%    'regularization_list': list of regularization weights for each class
%       and subclass, stacked for both classes
%    'regularization_profile' : (optional) regularization profile as matrix
%
%Description:
%   TRAIN_RLSDA trains a Relevance Subclass LDA classifier on data X with
%   class labels given in LABELS and Subclass labels given in SUBLABELS.
%   The mean shrinkage parameter is selected by sotm and the covariance
%   shrinkage parameter is selected by the function clsutil_shrinkage.
%
%
%Example 1:
%%perform standard ERP analysis (loads data, subsampling and divide into epochs)
% demo_analysis_ERPs
%
%
%%only the last digit codes for the stimulus identity
% subclass_lab = mod(epo.event.desc, 10);
%
%%extract features from the epochs by subsampling
% ivals = [140:30:260 300:50:600; 170:30:300 350:50:650]';
%%create features
% fv = proc_flaten(proc_jumpingMeans(epo, ivals));
%%train RSLDA
% C = train_RSLDAshrink(fv.x, fv.y, subclass_lab);
%
%%Example 2: compare RSLDAshrink with LDAshrink
% nfolds = 4;
% [divTr, divTe] = sample_KFold(epo.y, nfolds);
% loss_list_RSLDA = nan(1, nfolds);loss_list_RLDA = loss_list_RSLDA;
% for kk = 1:nfolds
%     %generate training and test features
%     fvTr = proc_flaten(proc_jumpingMeans(proc_selectEpochs(epo, divTr{1}{kk}), ivals));
%     fvTe = proc_flaten(proc_jumpingMeans(proc_selectEpochs(epo, divTe{1}{kk}), ivals));
%     %train cls on traing data
%     C_RSLDA = train_RSLDAshrink(fvTr.x, fvTr.y, subclass_lab(divTr{1}{kk}));
%     C_LDA = train_RLDAshrink(fvTr.x, fvTr.y);
%     %apply RSLDA on the test data
%     out_RSLDA = apply_RSLDAshrink(C_RSLDA, fvTe.x, subclass_lab(divTe{1}{kk}));
%     out_RLDA = apply_separatingHyperplane(C_LDA , fvTe.x);
%
%     loss_list_RSLDA(kk) = loss_rocArea(fvTe.y,out_RSLDA);
%     loss_list_RLDA(kk) = loss_rocArea(fvTe.y,out_RLDA);
% end
%%compare loss for each fold
% loss_list_RSLDA - loss_list_RLDA
%%plot regularization profile for the last fold
% figure, imagesc(C_RSLDA.regularization_profile{2})
%
%See also:
%   APPLY_SEPARATINGHYPERPLANE, CLSUTIL_MEANMTS, CLSUTIL_SHRINKAGE, TRAIN_RLDAshrink
%   TRAIN_LDA

% 02-2015 Johannes Hoehne
if size(yTr,1)==1, yTr= [yTr<0; yTr>0]; end
yTr = logical(yTr);
nClasses= size(yTr,1);
props= {'Whitening'        1                             'BOOL'
        'ReturnRegularizationProfile'  1                 'BOOL'
	    'Prior'            ones(nClasses, 1)/nClasses    'DOUBLE[- 1]'
       };

% get props list of the subfunction
props_shrinkage= clsutil_meanMTS;

if nargin==0,
  C= opt_catProps(props, props_shrinkage);
  return
end

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props, props_shrinkage);

% preprocessing of input
yTr = logical(yTr);
if size(sublab,2)==1 && size(sublab,1)>1
    sublab = sublab'; %sublab is expected to be a row-vector of
end

% empirical class priors as an option, if nan (I leave 1/nClasses as default)
if isnan(opt.Prior)
  opt.Prior = sum(yTr, 2)/sum(sum(yTr));
end

C = {};

%start func
C.sublab_unique = unique(sublab);

% get meanfree X (without considereing subclass labels)
% compute Cov matrix
if opt.Whitening
    M1 = mean(xTr(:, yTr(1,:)),2);
    M2 = mean(xTr(:,yTr(2,:)),2);
    M = repmat(M1, 1, size(yTr,2));
    M(:, yTr(2,:)) = repmat(M2, 1, sum(yTr(2,:)));
    Xmeanfree_dum = xTr - M;
    globalCov = clsutil_shrinkage(Xmeanfree_dum);
    [ Eigvec, Eigval  ] = eig(globalCov);
    A_feat2white =  Eigvec * (Eigval^-.5) * Eigvec;
    A_white2feat = inv(A_feat2white);
else %do not perform whitening
%     --> this leads to a poor performance for data with a skewed eigen spectrum
    A_feat2white = eye(size(xTr,2));
    A_white2feat = eye(size(xTr,2));
end
Xwhite = A_feat2white' * xTr;

kk = 0;
M = nan(size(xTr)); %saves the classwise-shrinked means
for my_sublab = C.sublab_unique
    kk = kk+1;

    for ix_c = 1:nClasses %for each class
        %% estimate class means
        shrink_dat = {}; %initialize shrindat container
        shrink_dat.X = {Xwhite(:,((sublab == my_sublab) & yTr(ix_c,:)))' };

        %there might be sublabels without epoch in calss ix_c -> skip them as sublabel
        sublab_missing = [];
        for sublab_this = setdiff(C.sublab_unique, my_sublab)
            if sum((((sublab == sublab_this)) & yTr(ix_c,:)))>0
                shrink_dat.X{end+1} = Xwhite(:,(((sublab == sublab_this)) & yTr(ix_c,:)))';
            else
                sublab_missing(end+1) = sublab_this; %(only important for special cases)
            end
        end

        %PERFORM MULTI-TARGET SHINKAGE FOR CLASS 1
        % sublab=k & T   sublab~=k & T
        [Mest(ix_c,:), gamma_tmp] = clsutil_meanMTS(shrink_dat, 'convex', 1, 'variablewise', 0, 'conservative', 1);

        % (only important for special cases) deal with the arrangement of missing subclass labels
        for xx = sublab_missing %put missing subclasses as zeros in gamma to maintain the correct size
            if xx>=length(gamma_tmp)+1, gamma_tmp = [gamma_tmp; 0]; %the last one was missing
            else  gamma_tmp = [gamma_tmp(1:xx-1); 0; gamma_tmp(xx:end)]; %another one was missing
            end
        end
        gamma(ix_c,:) = gamma_tmp;

        %transform the Means back into feature space
        Mest_origSpace(ix_c,:) = A_white2feat' * Mest(ix_c,:)';

        % for each datapoint, save the dcorresponding subclass mean
        % --> build up Mean-template M with Target and NonTarget template
        M(:, find(yTr(ix_c,:) & sublab == my_sublab)) = repmat(Mest_origSpace(ix_c,:)', 1, length(find(yTr(ix_c,:) & sublab == my_sublab)));

        %save the mean estimators as we need them later to compute the cls
        list_M{kk}(ix_c,:) = Mest_origSpace(ix_c,:);
    end
    C.regularization_list(:,kk) =  gamma(:);
end

%% estimate Cov with optimized means
Xmeanfree = xTr - M;
Cest = clsutil_shrinkage(Xmeanfree);
C_invcov = pinv(Cest);

%% finalize cls
kk = 0;
for my_sublab = C.sublab_unique
    kk = kk+1;
    if nClasses>2
        C.subC{kk} =   comp_reg_cls(C_invcov, list_M{kk}', opt.Prior) ;
    else
        % get the Target and NonTartget Means for each subclass
        M1est = list_M{kk}(1,:)';
        M2est = list_M{kk}(2,:)';
        C.pattern(:,kk) = M2est - M1est;
        C.subC{kk} = comp_reg_cls2(C_invcov, M1est, M2est);
    end
    if ~isfinite(C.subC{kk}.b)
        error('sth went wrong')
    end
end

if opt.ReturnRegularizationProfile
    %reorganize the regularization_list such that we can interpret them
%                the regul parameters are reordered such that (dat)
%                    l1  l2
%                    l3  l4
%                    l5  l6
%                becomes a square Matrix A, which is better to interpret:
%                    (1-l1-l2)  l        l
%                      l3   (1-l3-l4)    l4
%                      l5       l6    (1-l5-l6)
    for i_class = 1:nClasses %for both classes
        nPara = size(C.regularization_list,1)/nClasses;
        dat = C.regularization_list((i_class-1)*nPara+1 : (i_class)*nPara,:)';
        A = zeros(size(dat,1));
        A(logical(eye(size(A)))) = 1-sum(dat,2);
        d = dat';
        A(~logical(eye(size(A)))) = d(:);
        C.regularization_profile{i_class} = A';
    end
end

return

end

function C = comp_reg_cls2(C_invcov, M1, M2)
C = [];
C.w= C_invcov*(M2 - M1);
C.b= -0.5*C.w'*(M2 + M1);
end

function C = comp_reg_cls(C_invcov, C_mean, Prior)
C = [];
C.w= C_invcov*C_mean;
C.b= -0.5*sum(C_mean.*C.w,1)' + log(Prior);
end

author: bbci


text:
function i_rect = im_project(im,ps,rr,dim)

rr_min   = min(rr(1,:));
rr_max   = max(rr(2,:));
rr_width = uint32(rr_max-rr_min);
nr_lines = size(rr,2);

i_rect = zeros(dim',class(im));

skipped=0;
for line = 1:nr_lines
    for jj = 1:size(ps{line},2)
        x = ps{line}(1,jj);
        y = ps{line}(2,jj);
        [val,flag] = get_pixel(im,[x;y]);
        if flag == 1
            if dim(3) == 1
                i_rect(line,uint32(rr(1,line)-rr_min)+jj) = val;
            else
                i_rect(line,uint32(rr(1,line)-rr_min)+jj,:) = val;
            end
            skipped = skipped + 1;
        end
    end
end

fprintf('skipped %d of total pixels\n',skipped/(rr_width*nr_lines));
fprintf('original image: %g by %g = %d pixls\n',dim(1),dim(2),prod(dim));
fprintf('new image: %d by %d = %d pixels\n',rr_width,nr_lines,rr_width*nr_lines);
fprintf('orig_size/new_size = %f\n',double(prod(dim))/double((rr_width*nr_lines)));

author: Sable


text:
function [ result ] = CREST_tracking( opts, varargin, config, display)
%Test function
global objSize;

LocGt=config.gt;
num_channels=64;

% training options (SGD)
opts.train = struct([]) ;
[opts, varargin] = vl_argparse(opts, varargin) ;

output_sigma_factor = 0.1;

gpuDevice(1);
[net1,avgImg]=initVGG16Net();

nFrame=config.nFrames;
[Gt]=config.gt;

result = zeros(length(nFrame), 4); result(1,:) = Gt(1,:);

% I only have a Titan Black GPU on my side. Its memory is limited.
% For objects with large size, I have to perform the downsampling.
% You can remove the downsampling on your side. The results might be little different.
scale=1;
global resize;
objSize=Gt(1,3:4);
if objSize(1)*objSize(2)>resize*resize
    scale=resize/max(objSize);
    disp('resized');
end

im1=imread(config.imgList{1});
im=imresize(im1,scale);
cell_size=4;
if size(im,3)==1
    im = cat(3, im, im, im);
    im1 = cat(3, im1, im1, im1);
end

targetLoc=round(Gt(1,:)*scale);
target_sz=[targetLoc(4) targetLoc(3)];
im_sz=size(im);
window_sz = get_search_window(target_sz, im_sz);
l1_patch_num = ceil(window_sz/ cell_size);
l1_patch_num=l1_patch_num-mod(l1_patch_num,2)+1;
cos_window = hann(l1_patch_num(1)) * hann(l1_patch_num(2))';

sz_window=size(cos_window);
pos = [targetLoc(2), targetLoc(1)] + floor(target_sz/2);
patch = get_subwindow(im, pos, window_sz);
meanImg=zeros(size(patch));
meanImg(:,:,1)=avgImg(1);
meanImg(:,:,2)=avgImg(2);
meanImg(:,:,3)=avgImg(3);
patch1 = single(patch) - meanImg;
net1.eval({'input',gpuArray(patch1)});

index=[23];
feat=cell(length(index),1);
for i=1:length(index)
    feat1 = gather(net1.vars(index(i)).value);
    feat1 = imResample(feat1, sz_window(1:2));
    feat{i} = bsxfun(@times, feat1, cos_window);
end
feat=feat{1};

[hf,wf,cf]=size(feat);
matrix=reshape(feat,hf*wf,cf);
coeff = pca(matrix);
coeff=coeff(:,1:num_channels);

target_sz1=ceil(target_sz/cell_size);
output_sigma = target_sz1*output_sigma_factor;
label=gaussian_shaped_labels(output_sigma, l1_patch_num);

label1=imresize(label,[size(im1,1) size(im1,1)])*255;
patch1=imresize(patch,[size(im1,1) size(im1,1)]);
imd=[im1];
%-------------------Display First frame----------
if display
    figure(2);
    set(gcf,'Position',[200 300 480 320],'MenuBar','none','ToolBar','none');
    hd = imshow(imd,'initialmagnification','fit'); hold on;
    rectangle('Position', Gt(1,:), 'EdgeColor', [0 0 1], 'Linewidth', 1);
    set(gca,'position',[0 0 1 1]);
    text(10,10,'1','Color','y', 'HorizontalAlignment', 'left', 'FontWeight','bold', 'FontSize', 30);
    hold off;
    drawnow;
end

%-------------------first frame initialization-----------
trainOpts.numEpochs=4000;

feat_=reshape(feat,hf*wf,cf);
feat_=feat_*coeff;
featPCA=reshape(feat_,hf,wf,num_channels);

net_online=initNet(target_sz1);

trainOpts.batchSize = 1 ;
trainOpts.numSubBatches = 1 ;
trainOpts.continue = true ;
trainOpts.gpus = 1 ;
trainOpts.prefetch = true ;

trainOpts.expDir = opts.expDir ;
trainOpts.learningRate=5e-8;
trainOpts.weightDecay= 1;

train=1;
imdb=[];
input={featPCA label};
featPCA1st=featPCA;
opts.train.gpus=1;
bopts.useGpu = numel(opts.train.gpus) > 0 ;
info = cnn_train_dag(net_online, imdb, input,getBatchWrapper(bopts), ...
                     trainOpts, ...
                     'train', train, ...
                     opts.train) ;

net_online.move('gpu');

%----------------online prediction------------------
motion_sigma_factor=0.6;
cell_size=4;
global num_update;
num_update=2;
cur=1;
feat_update=cell(num_update,1);
label_update=cell(num_update,1);
target_szU=target_sz;
for i=2:nFrame
    im1=imread(config.imgList{i});
    im=imresize(im1,scale);
    if size(im1,3)==1
        im = cat(3, im, im, im);
        im1 = cat(3, im1, im1, im1);
    end

    patch = get_subwindow(im, pos, window_sz);
    patch1 = single(patch) - meanImg;
    net1.eval({'input',gpuArray(patch1)});

    feat=cell(length(index),1);
    for j=1:length(index)
        feat1 = gather(net1.vars(index(j)).value);
        feat1 = imResample(feat1, sz_window(1:2));
        feat{j} = bsxfun(@times, feat1, cos_window);
    end
    feat=feat{1};

    feat_=reshape(feat,hf*wf,cf);
    feat_=feat_*coeff;
    featPCA=reshape(feat_,hf,wf,num_channels);

    net_online.eval({'input1',gpuArray(featPCA),'input2',gpuArray(featPCA1st)});
    regression_map=gather(net_online.vars(10).value);

    motion_sigma = target_sz1*motion_sigma_factor;
    motion_map=gaussian_shaped_labels(motion_sigma, l1_patch_num);

    response=regression_map.*motion_map;

    [vert_delta, horiz_delta] = find(response == max(response(:)), 1);
    vert_delta  = vert_delta  - ceil(hf/2);
    horiz_delta = horiz_delta - ceil(wf/2);

    pos = pos + cell_size * [vert_delta, horiz_delta];

    target_szU=scale_estimation(im,pos,target_szU,window_sz,...
            net1,net_online,coeff,meanImg,featPCA1st);

    targetLoc=[pos([2,1]) - target_szU([2,1])/2, target_szU([2,1])];
    result(i,:)=round(targetLoc/scale);

    label1=imresize(regression_map,[size(im1,1) size(im1,1)])*255;
    patch1=imresize(patch,[size(im1,1) size(im1,1)]);
    imd=[im1];
%    -----------Display current frame-----------------
    if display
        hc = get(gca, 'Children'); delete(hc(1:end-1));
        set(hd,'cdata',imd); hold on;
        rectangle('Position', result(i,:), 'EdgeColor', [0 0 1], 'Linewidth', 1);
        set(gca,'position',[0 0 1 1]);
        text(10,10,num2str(i),'Color','y', 'HorizontalAlignment', 'left', 'FontWeight','bold', 'FontSize', 30);
        hold off;
        drawnow;
    end

     %-----------Model update-------------------------
     labelU=circshift(label,[vert_delta,horiz_delta]);
     feat_update{cur}=featPCA;
     label_update{cur}=labelU;

     if cur==num_update

        trainOpts.batchSize = 1 ;
        trainOpts.numSubBatches = 1 ;
        trainOpts.continue = true ;
        trainOpts.gpus = 1 ;
        trainOpts.prefetch = true ;

        trainOpts.expDir = 'exp/update/' ;
        trainOpts.learningRate = 2e-9;
        trainOpts.weightDecay= 1;
        trainOpts.numEpochs = 2;

        train=1;
        imdb=[];
        input={feat_update featPCA1st label_update};
        opts.train.gpus=1;
        bopts.useGpu = numel(opts.train.gpus) > 0 ;

        info = cnn_train_dag_update(net_online, imdb, input,getBatchWrapper(bopts), ...
                             trainOpts, ...
                             'train', train, ...
                             opts.train) ;
        net_online.move('gpu');

        cur=1;
    else
        cur=cur+1;
    end

end

end

function fn = getBatchWrapper(opts)
% -------------------------------------------------------------------------
fn = @(imdb,batch) getBatch(imdb,batch,false,opts,'prefetch',nargout==0) ;
end

repo: CREST-Release


text:
function [lines, spdata, edgeIm] = APPgetLargeConnectedEdges(grayIm, minLen, imsegs)
% [lines, spdata, edgeIm] = APPgetLargeConnectedEdges(grayIm, minLen,
%                                 imsegs)
%
% Uses the method of Video Compass [Kosecka, et al 2002] to get long,
% straight edges.
%
% Input:
%   grayIm: grayscale image to be analyzed
%   minLen: minimum length in pixels for edge (suggested 0.025*diagonal)
%   imsegs: superpixel image structure (used to store sp statistics)
% Output:
%   lines: parameters for long, straight lines
%   spdata: statistics for lines in each superpixel
%   edgeIm (optional): for displaying lines found
%
% Copyright(C) Derek Hoiem, Carnegie Mellon University, 2005
% Permission granted to non-commercial enterprises for
% modification/redistribution under GNU GPL.
% Current Version: 1.0  09/30/2005

[dX, dY] = gradient(conv2(grayIm, fspecial('gaussian', 7, 1.5), 'same'));

im_canny = edge(grayIm, 'canny');
% remove border edges
im_canny([1 2 end-1 end], :) = 0;
im_canny(:, [1 2 end-1 end]) = 0;
width = size(im_canny, 2);
height = size(im_canny, 1);

ind = find(im_canny > 0);

num_dir = 8;

dX_ind = dX(ind);
dY_ind = dY(ind);
a_ind = atan(dY_ind ./ (dX_ind+1E-10));

% a_ind ranges from 1 to num_dir with bin centered around pi/2
a_ind = ceil(mod(a_ind/pi*num_dir-0.5, num_dir));
%[g, gn] = grp2idx(a_ind);

% get the indices of edges in each direction
for i = 1:num_dir
    direction(i).ind = ind(find(a_ind==i));
end

% remove edges that are too small and give all edges that have the same
% direction a unique id
% edges(height, width, [angle id])
if nargout>2
    edgeIm = grayIm*0.75;
end
lines = zeros(2000, 6);

nspdata = length(imsegs.npixels);
spdata = repmat(struct('lines', zeros(5, 1), 'edge_count', 0), nspdata,1);
bcount = zeros(nspdata, 1);

used = zeros(size(im_canny));

line_count = 0;
for k = 1:num_dir

    num_ind = 0;
    for m = (k-1):k+1
        num_ind = num_ind + sum(~used(direction(mod(m-1, num_dir)+1).ind));
    end

    ind = zeros(num_ind, 1);
    dir_im = zeros(size(im_canny));

    count = 0;
    for m = (k-1):k+1
        m2 = mod(m-1, num_dir)+1;
        tind = direction(m2).ind(find(~used(direction(m2).ind)));
        tmpcount = length(tind);
        ind(count+1:count+tmpcount) = tind;
        count = count + tmpcount;
    end
    dir_im(ind) = 1;

    [tmpL, num_edges] = bwlabel(dir_im, 8);

    % get the number of pixels in each edge
    edge_size = zeros(num_edges, 1);
    edges = repmat(struct('ind', zeros(200, 1)), num_edges, 1);
    for i = 1:length(ind)
        id = tmpL(ind(i));
        edge_size(id) = edge_size(id) + 1;
        edges(id).ind(edge_size(id)) = ind(i);
    end
    for i = 1:num_edges
        edges(i).ind = edges(i).ind(1:edge_size(i));
    end

    % get the endpoints of the long edges and an image of the long edges
    for id = 1:num_edges
        if edge_size(id) > minLen

            y = mod(edges(id).ind-1, height)+1;
            x = floor((edges(id).ind-1) / height)+1;

            %by = min(floor(y/block_size), ny-1);
            %bx = min(floor(x/block_size), nx-1);

            mean_x = mean(x);
            mean_y = mean(y);
            zmx = (x-mean_x);
            zmy = (y-mean_y);
            D = [sum(zmx.^2) sum(zmx.*zmy); sum(zmx.*zmy) sum(zmy.^2)];
            [v, lambda] = eig(D);
            theta = atan2(v(2, 2) , v(1, 2));
            if lambda(1,1)>0
                conf = lambda(2,2)/lambda(1,1);
            else
                conf = 100000;
            end

            %disp(conf)

            if conf >= 400
                line_count = line_count+1;

                used(edges(id).ind) = 1;
                bi = double(imsegs.segimage(edges(id).ind));
                [g, gn] = grp2idx(bi);
                for k = 1:length(bi)
                    if bi(k) > 0
                        spdata(bi(k)).edge_count = spdata(bi(k)).edge_count + 1;
                    end
                end
                for k = 1:length(gn)
                    tmpbi = str2num(gn{k});
                    if tmpbi>0
                        bcount(tmpbi) = bcount(tmpbi)+1;
                        spdata(tmpbi).lines(bcount(tmpbi)) = line_count;
                    end
                end

                %disp(num2str([lambda(1,1) lambda(2,2)]))
                r = sqrt((max(x)-min(x))^2 + (max(y)-min(y))^2);
                x1 = mean_x - cos(theta)*r/2;
                x2 = mean_x + cos(theta)*r/2;
                y1 = mean_y - sin(theta)*r/2;
                y2 = mean_y + sin(theta)*r/2;

                r = mean_x*cos(theta)+mean_y*sin(theta);
                %tr = x1*cos(theta) + y1*sin(theta);
                %disp(num2str([r tr]))

                lines(line_count, 1:6) = [x1 x2 y1 y2 theta r];

            end
        end
    end

end

if nargout>2
    for i = 1:line_count
        edgeIm = draw_line_image2(edgeIm, lines(i, 1:4)', i);
    end
end

for k = 1:length(spdata)
    spdata(k).lines = spdata(k).lines(1:bcount(k))';
end
lines = lines(1:line_count, :);

imsize = size(grayIm);
lines = normalizeLines(lines, imsize(1:2));

author: zouchuhang

repo: LayoutNet