text:
function submit(partId, webSubmit)
%SUBMIT Submit your code and output to the pgm-class servers
%   SUBMIT() will connect to the pgm-class server and submit your solution
%   There is no penalty for submitting, so go ahead and try this!
%
%   If this function does not work for you, use the web-submission mechanism.
%   Call the submitWeb function, which will produce a file for the part you
%   wish to submit. Then, submit the file to the class servers using the
%   "Web Submission" button on the Programming Assignments page on the course
%   website.
%
%   webSubmit is a boolean variable that specifies whether to prepare
%   a file for web-submission (if webSubmit = 1) or to submit directly
%   to the server (if webSubmit = 0, or is not specified); you should call
%   submitWeb if you want to do a web-submission.
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

fprintf('======================================================\n');
fprintf('                 Honor Code Statement:\n');
fprintf('------------------------------------------------------\n');
fprintf('The Honor Code is an undertaking of the students,\n');
fprintf('individually and collectively, that they will not\n');
fprintf('give or receive unpermitted aid in class work and\n');
fprintf('will do their share and take an active part in seeing\n');
fprintf('to it that others as well as themselves uphold the\n');
fprintf('spirit and letter of the Honor Code.\n\n');
fprintf('Sharing, collaboration, or looking at any code related\n');
fprintf('to Programming Assignments that is not your own is\n');
fprintf('considered a violation of the Honor Code.\n');
fprintf('------------------------------------------------------\n');
fprintf('By submitting the Programming Assignment, I\n');
fprintf('acknowledge and accept the Honor Code.\n');
fprintf('======================================================\n\n');

  fprintf('==\n== [pgm-class] Submitting Solutions | Programming Assignment %s\n==\n', ...
	  homework_id());

  if ~exist('partId', 'var') || isempty(partId)
    partId = promptPart();
  end

  if ~exist('webSubmit', 'var') || isempty(webSubmit)
    webSubmit = 0; % submit directly by default
  end

  % Check valid partId
  partNames = validParts();
  if ~isValidPartId(partId)
    fprintf('!! Invalid assignment part selected.\n');
    fprintf('!! Expected an integer from 1 to %d.\n', numel(partNames) + 1);
    fprintf('!! Submission Cancelled\n');
    return
  end

  if ~exist('pgm_login_data.mat','file')
    [login password] = loginPrompt();
    save('pgm_login_data.mat','login','password');
  else
    load('pgm_login_data.mat');
    [login password] = quickLogin(login,password);
    save('pgm_login_data.mat','login','password');
  end

  if isempty(login)
    fprintf('!! Submission Cancelled\n');
    return
  end

  fprintf('\n== Connecting to pgm-class ... ');
  if exist('OCTAVE_VERSION')
    fflush(stdout);
  end

  % Setup submit list
  if partId == numel(partNames) + 1
    submitParts = 1:numel(partNames);
  else
    submitParts = [partId];
  end

  for s = 1:numel(submitParts)
    % Submit this part
    partId = submitParts(s);

    for thisPartId = subParts(partId)
      if (~webSubmit) % submit directly to server

        [login, ch, signature, auxstring] = getChallenge(login, thisPartId);
        if isempty(login) || isempty(ch) || isempty(signature)
          % Some error occured, error string in first return element.
          fprintf('\n!! Error: %s\n\n', login);
          return
        end

        % Attempt Submission with Challenge
        ch_resp = challengeResponse(login, password, ch);

        [result, str] = submitSolution(login, ch_resp, thisPartId, ...
               output(thisPartId, auxstring), source(partId), signature);

        if (~isTest(thisPartId))
          partName = partNames{partId};
        else
          partName = [partNames{partId} ' (test)'];
        end

        fprintf('\n== [pgm-class] Submitted Assignment %s - Part %d - %s\n', ...
          homework_id(), partId, partName);
        fprintf('== %s\n', strtrim(str));

        if exist('OCTAVE_VERSION')
          fflush(stdout);
        end

      else % make web submission files

        [result] = submitSolutionWeb(login, thisPartId, output(thisPartId), ...
                            source(partId));
        result = base64encode(result);

        if (~isTest(thisPartId))
          partType = 'sample';
        else
          partType = 'test';
        end

        fprintf('\nSave as submission file [submit_pa%s_part%d_%s.txt (enter to accept default)]:', ...
          homework_id(), partId, partType);
        saveAsFile = input('', 's');
        if (isempty(saveAsFile))
          saveAsFile = sprintf('submit_pa%s_part%d_%s.txt', homework_id(), partId, partType);
        end

        fid = fopen(saveAsFile, 'w');
        if (fid)
          fwrite(fid, result);
          fclose(fid);
          fprintf('\nSaved your solutions to %s.\n\n', saveAsFile);
          fprintf(['You can now submit your solutions using the ' ...
             'Web Submission button\non the Assignments page\n']);
        else
          fprintf('Unable to save to %s\n\n', saveAsFile);
          fprintf(['You can create a submission file by saving the \n' ...
             'following text in a file: (press enter to continue)\n\n']);
          pause;
          fprintf(result);
        end

      end
    end
  end
end

% ====================== SERVER CONFIGURATION ===========================

% ***************** REMOVE -staging WHEN YOU DEPLOY *********************
function url = site_url()
  url = 'http://www.coursera.org/pgm';
end

function url = challenge_url()

  url = [site_url() '/assignment/challenge'];
end

function url = submit_url()
  url = [site_url() '/assignment/submit'];
end

% ================== CONFIGURABLES FOR EACH HOMEWORK ==================

function id = homework_id()
  id = '7';
end

function [partNames] = validParts()
  partNames = { 'StochasticGradientDescent', ...
                'LRSearchLambdaSGD', ...
                'LogZ', ...
                'CRFNegLogLikelihood', ...
                'CRFGradient'
              };
end

function srcs = sources()
  % Separated by part
  srcs = { { 'StochasticGradientDescent.m' }, ...
           { 'LRSearchLambdaSGD.m' }, ...
           { 'CliqueTreeCalibrate.m' }, ...
           { 'InstanceNegLogLikelihood.m'}, ...
           { 'InstanceNegLogLikelihood.m'}
         };
end

% defines the shown part to back-end part mappings
function parts = subParts(partId)
  first = 2 * (partId - 1) + 1;
  parts = [first, first + 1];
end

% specifies which parts are test parts
function result = isTest(partId)
  if (mod(partId, 2) == 0)
      result = true;
  else
      result = false;
  end
end

function out = output(partId, auxstring)
if partId == 1
  load 'Train1X.mat';
  load 'Train1Y.mat';
  out = SerializeVector(LRTrainSGD(Train1X, Train1Y, 1));
elseif partId == 2
  load 'Train2X.mat';
  load 'Train2Y.mat';
  out = SerializeVector(LRTrainSGD(Train2X, Train2Y, 1));
elseif partId == 3
  load 'Train1X.mat';
  load 'Train1Y.mat';
  load 'Validation1X.mat';
  load 'Validation1Y.mat';
  Test3Lambdas = [2 8];
  out = SerializeVector(LRSearchLambdaSGD(Train1X, Train1Y, Validation1X, Validation1Y, Test3Lambdas));
elseif partId == 4
  load 'Train2X.mat';
  load 'Train2Y.mat';
  load 'Validation2X.mat';
  load 'Validation2Y.mat';
  Test4Lambdas = [2 8];
  out = SerializeVector(LRSearchLambdaSGD(Train2X, Train2Y, Validation2X, Validation2Y, Test4Lambdas));
elseif partId == 5
    load Part2Sample.mat;
    [dummy, logZ] = CliqueTreeCalibrate(sampleUncalibratedTree, false);
    out = num2str(logZ);
elseif partId == 6
    load Part2LogZTest.mat;
    [dummy, logZ] = CliqueTreeCalibrate(logZTestCliqueTree, false);
    out = num2str(logZ);
elseif partId == 7
    load Part2Sample.mat;
    [nll, grad] = InstanceNegLogLikelihood(sampleX, sampleY, sampleTheta, sampleModelParams);
    out = num2str(nll);
elseif partId == 8
    load Part2Test.mat;
    [nll, grad] = InstanceNegLogLikelihood(testX, testY, testTheta, testModelParams);
    out = num2str(nll);
elseif partId == 9
    load Part2Sample.mat;
    [nll, grad] = InstanceNegLogLikelihood(sampleX, sampleY, sampleTheta, sampleModelParams);
    out = SerializeVector(grad);
elseif partId == 10
    load Part2Test.mat;
    [nll, grad] = InstanceNegLogLikelihood(testX, testY, testTheta, testModelParams);
    out = SerializeVector(grad);
end

end

function out = SerializeVector(x)
  % Serializes a numeric (floating point) vector.
  numLines = length(x);
  lines = cell(numLines,1);
  for i=1:numLines
    lines{i} = sprintf('%.4f\n', x(i));
  end
  out = sprintf('%s', lines{:});
end

function out = SerializeIntVector(x)
  % Serializes an integer vector.
  numLines = length(x);
  lines = cell(numLines,1);
  for i=1:numLines
    lines{i} = sprintf('%d\n', x(i));
  end
  out = sprintf('%s', lines{:});
end

% ========================= CHALLENGE HELPERS =========================

function src = source(partId)
  src = '';
  src_files = sources();
  if partId <= numel(src_files)
      flist = src_files{partId};
      for i = 1:numel(flist)
          fid = fopen(flist{i});
          if (fid == -1)
            error('Error opening %s (is it missing?)', flist{i});
          end
          line = fgets(fid);
          while ischar(line)
            src = [src line];
            line = fgets(fid);
          end
          fclose(fid);
          src = [src '||||||||'];
      end
  end
end

function ret = isValidPartId(partId)
  partNames = validParts();
  ret = (~isempty(partId)) && (partId >= 1) && (partId <= numel(partNames) + 1);
end

function partId = promptPart()
  fprintf('== Select which part(s) to submit:\n', ...
          homework_id());
  partNames = validParts();
  srcFiles = sources();
  for i = 1:numel(partNames)
    fprintf('==   %d) %s [', i, partNames{i});
    fprintf(' %s ', srcFiles{i}{:});
    fprintf(']\n');
  end
  fprintf('==   %d) All of the above \n==\nEnter your choice [1-%d]: ', ...
          numel(partNames) + 1, numel(partNames) + 1);
  selPart = input('', 's');
  partId = str2num(selPart);
  if ~isValidPartId(partId)
    partId = -1;
  end
end

function [email,ch,signature,auxstring] = getChallenge(email, part)
  str = urlread(challenge_url(), 'post', {'email_address', email, 'assignment_part_sid', [homework_id() '-' num2str(part)], 'response_encoding', 'delim'});

  str = strtrim(str);
  r = struct;
  while(numel(str) > 0)
    [f, str] = strtok (str, '|');
    [v, str] = strtok (str, '|');
    r = setfield(r, f, v);
  end

  email = getfield(r, 'email_address');
  ch = getfield(r, 'challenge_key');
  signature = getfield(r, 'state');
  auxstring = getfield(r, 'challenge_aux_data');
end

function [result, str] = submitSolutionWeb(email, part, output, source)

  result = ['{"assignment_part_sid":"' base64encode([homework_id() '-' num2str(part)], '') '",' ...
            '"email_address":"' base64encode(email, '') '",' ...
            '"submission":"' base64encode(output, '') '",' ...
            '"submission_aux":"' base64encode(source, '') '"' ...
            '}'];
  str = 'Web-submission';
end

function [result, str] = submitSolution(email, ch_resp, part, output, ...
                                        source, signature)

  params = {'assignment_part_sid', [homework_id() '-' num2str(part)], ...
            'email_address', email, ...
            'submission', base64encode(output, ''), ...
            'submission_aux', base64encode(source, ''), ...
            'challenge_response', ch_resp, ...
            'state', signature};

  str = urlread(submit_url(), 'post', params);

  % Parse str to read for success / failure
  result = 0;

end

% =========================== LOGIN HELPERS ===========================

function [login password] = loginPrompt()
  % Prompt for password
  [login password] = basicPrompt();

  if isempty(login) || isempty(password)
    login = []; password = [];
  end
end

function [login password] = basicPrompt()
  login = input('Login (Email address): ', 's');
  password = input('Submission Password (from Assignments page): ', 's');
end

function [login password] = quickLogin(login,password)
  cont_token = input(['You are currently logged in as ' login '.\nIs this you? (y/n - type n to reenter password)'],'s');
  if(isempty(cont_token) || cont_token(1)=='Y'||cont_token(1)=='y')
    return;
  else
    [login password] = loginPrompt();
  end
end

function [str] = challengeResponse(email, passwd, challenge)
  str = sha1([challenge passwd]);
end

% =============================== SHA-1 ================================

function hash = sha1(str)

  % Initialize variables
  h0 = uint32(1732584193);
  h1 = uint32(4023233417);
  h2 = uint32(2562383102);
  h3 = uint32(271733878);
  h4 = uint32(3285377520);

  % Convert to word array
  strlen = numel(str);

  % Break string into chars and append the bit 1 to the message
  mC = [double(str) 128];
  mC = [mC zeros(1, 4-mod(numel(mC), 4), 'uint8')];

  numB = strlen * 8;
  if exist('idivide')
    numC = idivide(uint32(numB + 65), 512, 'ceil');
  else
    numC = ceil(double(numB + 65)/512);
  end
  numW = numC * 16;
  mW = zeros(numW, 1, 'uint32');

  idx = 1;
  for i = 1:4:strlen + 1
    mW(idx) = bitor(bitor(bitor( ...
                  bitshift(uint32(mC(i)), 24), ...
                  bitshift(uint32(mC(i+1)), 16)), ...
                  bitshift(uint32(mC(i+2)), 8)), ...
                  uint32(mC(i+3)));
    idx = idx + 1;
  end

  % Append length of message
  mW(numW - 1) = uint32(bitshift(uint64(numB), -32));
  mW(numW) = uint32(bitshift(bitshift(uint64(numB), 32), -32));

  % Process the message in successive 512-bit chs
  for cId = 1 : double(numC)
    cSt = (cId - 1) * 16 + 1;
    cEnd = cId * 16;
    ch = mW(cSt : cEnd);

    % Extend the sixteen 32-bit words into eighty 32-bit words
    for j = 17 : 80
      ch(j) = ch(j - 3);
      ch(j) = bitxor(ch(j), ch(j - 8));
      ch(j) = bitxor(ch(j), ch(j - 14));
      ch(j) = bitxor(ch(j), ch(j - 16));
      ch(j) = bitrotate(ch(j), 1);
    end

    % Initialize hash value for this ch
    a = h0;
    b = h1;
    c = h2;
    d = h3;
    e = h4;

    % Main loop
    for i = 1 : 80
      if(i >= 1 && i <= 20)
        f = bitor(bitand(b, c), bitand(bitcmp(b), d));
        k = uint32(1518500249);
      elseif(i >= 21 && i <= 40)
        f = bitxor(bitxor(b, c), d);
        k = uint32(1859775393);
      elseif(i >= 41 && i <= 60)
        f = bitor(bitor(bitand(b, c), bitand(b, d)), bitand(c, d));
        k = uint32(2400959708);
      elseif(i >= 61 && i <= 80)
        f = bitxor(bitxor(b, c), d);
        k = uint32(3395469782);
      end

      t = bitrotate(a, 5);
      t = bitadd(t, f);
      t = bitadd(t, e);
      t = bitadd(t, k);
      t = bitadd(t, ch(i));
      e = d;
      d = c;
      c = bitrotate(b, 30);
      b = a;
      a = t;

    end
    h0 = bitadd(h0, a);
    h1 = bitadd(h1, b);
    h2 = bitadd(h2, c);
    h3 = bitadd(h3, d);
    h4 = bitadd(h4, e);

  end

  hash = reshape(dec2hex(double([h0 h1 h2 h3 h4]), 8)', [1 40]);

  hash = lower(hash);

end

function ret = bitadd(iA, iB)
  ret = double(iA) + double(iB);
  ret = bitset(ret, 33, 0);
  ret = uint32(ret);
end

function ret = bitrotate(iA, places)
  t = bitshift(iA, places - 32);
  ret = bitshift(iA, places);
  ret = bitor(ret, t);
end

% =========================== Base64 Encoder ============================
% Thanks to Peter John Acklam
%

function y = base64encode(x, eol)
%BASE64ENCODE Perform base64 encoding on a string.
%
%   BASE64ENCODE(STR, EOL) encode the given string STR.  EOL is the line ending
%   sequence to use; it is optional and defaults to '\n' (ASCII decimal 10).
%   The returned encoded string is broken into lines of no more than 76
%   characters each, and each line will end with EOL unless it is empty.  Let
%   EOL be empty if you do not want the encoded string broken into lines.
%
%   STR and EOL don't have to be strings (i.e., char arrays).  The only
%   requirement is that they are vectors containing values in the range 0-255.
%
%   This function may be used to encode strings into the Base64 encoding
%   specified in RFC 2045 - MIME (Multipurpose Internet Mail Extensions).  The
%   Base64 encoding is designed to represent arbitrary sequences of octets in a
%   form that need not be humanly readable.  A 65-character subset
%   ([A-Za-z0-9+/=]) of US-ASCII is used, enabling 6 bits to be represented per
%   printable character.
%
%   Examples
%   --------
%
%   If you want to encode a large file, you should encode it in chunks that are
%   a multiple of 57 bytes.  This ensures that the base64 lines line up and
%   that you do not end up with padding in the middle.  57 bytes of data fills
%   one complete base64 line (76 == 57*4/3):
%
%   If ifid and ofid are two file identifiers opened for reading and writing,
%   respectively, then you can base64 encode the data with
%
%      while ~feof(ifid)
%         fwrite(ofid, base64encode(fread(ifid, 60*57)));
%      end
%
%   or, if you have enough memory,
%
%      fwrite(ofid, base64encode(fread(ifid)));
%
%   See also BASE64DECODE.

%   Author:      Peter John Acklam
%   Time-stamp:  2004-02-03 21:36:56 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   if isnumeric(x)
      x = num2str(x);
   end

   % make sure we have the EOL value
   if nargin < 2
      eol = sprintf('\n');
   else
      if sum(size(eol) > 1) > 1
         error('EOL must be a vector.');
      end
      if any(eol(:) > 255)
         error('EOL can not contain values larger than 255.');
      end
   end

   if sum(size(x) > 1) > 1
      error('STR must be a vector.');
   end

   x   = uint8(x);
   eol = uint8(eol);

   ndbytes = length(x);                 % number of decoded bytes
   nchunks = ceil(ndbytes / 3);         % number of chunks/groups
   nebytes = 4 * nchunks;               % number of encoded bytes

   % add padding if necessary, to make the length of x a multiple of 3
   if rem(ndbytes, 3)
      x(end+1 : 3*nchunks) = 0;
   end

   x = reshape(x, [3, nchunks]);        % reshape the data
   y = repmat(uint8(0), 4, nchunks);    % for the encoded data

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Split up every 3 bytes into 4 pieces
   %
   %    aaaaaabb bbbbcccc ccdddddd
   %
   % to form
   %
   %    00aaaaaa 00bbbbbb 00cccccc 00dddddd
   %
   y(1,:) = bitshift(x(1,:), -2);                  % 6 highest bits of x(1,:)

   y(2,:) = bitshift(bitand(x(1,:), 3), 4);        % 2 lowest bits of x(1,:)
   y(2,:) = bitor(y(2,:), bitshift(x(2,:), -4));   % 4 highest bits of x(2,:)

   y(3,:) = bitshift(bitand(x(2,:), 15), 2);       % 4 lowest bits of x(2,:)
   y(3,:) = bitor(y(3,:), bitshift(x(3,:), -6));   % 2 highest bits of x(3,:)

   y(4,:) = bitand(x(3,:), 63);                    % 6 lowest bits of x(3,:)

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Now perform the following mapping
   %
   %   0  - 25  ->  A-Z
   %   26 - 51  ->  a-z
   %   52 - 61  ->  0-9
   %   62       ->  +
   %   63       ->  /
   %
   % We could use a mapping vector like
   %
   %   ['A':'Z', 'a':'z', '0':'9', '+/']
   %
   % but that would require an index vector of class double.
   %
   z = repmat(uint8(0), size(y));
   i =           y <= 25;  z(i) = 'A'      + double(y(i));
   i = 26 <= y & y <= 51;  z(i) = 'a' - 26 + double(y(i));
   i = 52 <= y & y <= 61;  z(i) = '0' - 52 + double(y(i));
   i =           y == 62;  z(i) = '+';
   i =           y == 63;  z(i) = '/';
   y = z;

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Add padding if necessary.
   %
   npbytes = 3 * nchunks - ndbytes;     % number of padding bytes
   if npbytes
      y(end-npbytes+1 : end) = '=';     % '=' is used for padding
   end

   if isempty(eol)

      % reshape to a row vector
      y = reshape(y, [1, nebytes]);

   else

      nlines = ceil(nebytes / 76);      % number of lines
      neolbytes = length(eol);          % number of bytes in eol string

      % pad data so it becomes a multiple of 76 elements
      y = [y(:) ; zeros(76 * nlines - numel(y), 1)];
      y(nebytes + 1 : 76 * nlines) = 0;
      y = reshape(y, 76, nlines);

      % insert eol strings
      eol = eol(:);
      y(end + 1 : end + neolbytes, :) = eol(:, ones(1, nlines));

      % remove padding, but keep the last eol string
      m = nebytes + neolbytes * (nlines - 1);
      n = (76+neolbytes)*nlines - neolbytes;
      y(m+1 : n) = '';

      % extract and reshape to row vector
      y = reshape(y, 1, m+neolbytes);

   end

   % output is a character array
   y = char(y);

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
% PLOT_ELLIPSE
% h=plot_ellipse(x,y,theta,a,b)
%
% This routine plots an ellipse with centre (x,y), axis lengths a,b
% with major axis at an angle of theta radians from the horizontal.

%
% Author: P. Fieguth
%         Jan. 98
%
%http://ocho.uwaterloo.ca/~pfieguth/Teaching/372/plot_ellipse.m

function h=plot_ellipse(x,y,theta,a,b)

np = 100;
ang = [0:np]*2*pi/np;
R = [cos(theta) -sin(theta); sin(theta) cos(theta)];
pts = [x;y]*ones(size(ang)) + R*[cos(ang)*a; sin(ang)*b];
h=plot( pts(1,:), pts(2,:) );

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% lsconvex;
% check convexity
%

if nmin>1,
  convex=0;
else
  convex=1;
  for i=2:s-1,
    f12=(flist(i)-flist(i-1))/(alist(i)-alist(i-1));
    f13=(flist(i)-flist(i+1))/(alist(i)-alist(i+1));
    f123=(f13-f12)/(alist(i+1)-alist(i-1));
    if f123<0,
      if prt>1, disp('not convex'); end;
      convex=0;
      break;
    end;
  end;
  if prt>1 & convex, disp('convex'); end;
end;

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function [varargout] = nanmin(varargin)
%NANMIN Minimum value, ignoring NaNs.
%   M = NANMIN(A) returns the minimum of A with NaNs treated as missing.
%   For vectors, M is the smallest non-NaN element in A.  For matrices, M
%   is a row vector containing the minimum non-NaN element from each
%   column.  For N-D arrays, NANMIN operates along the first non-singleton
%   dimension.
%
%   [M,NDX] = NANMIN(A) returns the indices of the minimum values in A.  If
%   the values along the first non-singleton dimension contain more than
%   one minimal element, the index of the first one is returned.
%
%   M = NANMIN(A,B) returns an array the same size as A and B with the
%   smallest elements taken from A or B.  Either one can be a scalar.
%
%   [M,NDX] = NANMIN(A,[],DIM) operates along the dimension DIM.
%
%   See also MIN, NANMAX, NANMEAN, NANMEDIAN, NANVAR, NANSTD.

%   Copyright 1993-2004 The MathWorks, Inc.
%   $Revision: 1.1.8.1 $  $Date: 2010/03/16 00:15:52 $

% Call [m,ndx] = min(a,b) with as many inputs and outputs as needed
[varargout{1:nargout}] = min(varargin{:});

author: marianux

repo: ecg-kit

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


text:
filename='Gripping_quad_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'nonadjoint_compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
constraint_case = 'INEQUALITY';
optimizer = 'MMA'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.75;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%% Example using ESS for the BCIT to resample
ess2Dir = 'D:\TestData\LargData\NCTU\NCTU_DAS\Level2';
levelDerivedDir = 'D:\TestData\LargData\NCTU\NCTU_DAS\Level2_256Hz';

% ess2Dir = 'D:\TestData\BCITV2\Data\STDL2\ARL_BCIT_CalibrationDriving_Dataset_level_2';
% levelDerivedDir = 'D:\TestData\BCITV2\Data\STDL2_256Hz\ARL_BCIT_CalibrationDriving_Dataset_level_2';

% ess2Dir = 'D:\TestData\BCITV2\Data\STDL2\ARL_BCIT_TrafficComplexity_Dataset_level_2';
% levelDerivedDir = 'D:\TestData\BCITV2\Data\STDL2_256Hz\ARL_BCIT_TrafficComplexity_Dataset_level_2';

% ess2Dir = 'D:\TestData\LargData\VEP\ARL_VEP_v1.1.0_Level2';
% levelDerivedDir = 'D:\TestData\LargData\VEP\ARL_VEP_v1.1.0_Level2_256Hz';

% ess2Dir = 'D:\TestData\LargData\UCSD_RSVP\Level_2';
% levelDerivedDir = 'D:\TestData\LargData\UCSD_RSVP\Level_2_256Hz';

%% Check to make sure level 2 study validates
ess2File = [ess2Dir filesep 'studyLevel2_description.xml'];
obj1 = level2Study('level2XmlFilePath', ess2File);
obj1.validate();

%% Create a level 2 derived study
obj = levelDerivedStudy('parentStudyXmlFilePath', ess2File);
callbackAndParameters = {@resampleAndDealias, {'resampleOff', false, ...
                        'resampleFrequency', 256, 'lowPassFrequency', 100}};
obj = obj.createLevelDerivedStudy(callbackAndParameters, ...
     'filterDescription', 'Downsample and then low pass to remove alias', ...
     'filterLabel', 'resample', 'levelDerivedFolder', levelDerivedDir);

author: VisLab

repo: EEG-Clean-Tools


text:
% Fig. 5.27  Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
% script for Figure 5.27
n=[1 5.4];
d=conv([1 20],[1 1 0]);
rlocus(n,d)
hold on
r=roots([1 7 49]);
plot(r,'*')
axis([-20 4 -9 9])
title('Figure 5.27 Root locus for (s+5.4)/s(s+1)(s+20)')
z=0:.1:.9;
 wn=2:2:19;
 sgrid(z, wn)
 hold off

author: Sable

repo: mcbench-benchmarks


text:
function test_2D
%TEST_2D Test 2D Plot with DRAGZOOM

x = -pi*4:0.1:pi*4;
y1 = sin(x);
y1n = y1 + 0.5 * randn(1, length(x));
y2 = cos(x);
y2n = y2 + 0.5 * randn(1, length(x));

figure();
hax = axes();
plot(hax, x, y1, '.-b', x, y1n, '*g', x, y2, '.-r', x, y2n, 'om');
legend('plot 1', 'plot 2', 'plot 3', 'plot 4')

dragzoom();


text:
function M = Vec(M)
M = M(:);
end

author: KuangJuiHsu


text:
function [no,fc]=gridsurf(x,y,z, varargin)
%
% [no,fc]=gridsurf(x,y,z)
%    or
% [no,fc]=gridsurf(x,y,z, opt)
% [no,fc]=gridsurf(x,y,z, 'param1', value1, 'param2', value2, ...)
%
% convert a grid-shaped surface (used as input for surf) to a quad or triangular mesh
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%   x,y,z: parameter used as the input for surf()/mesh()
%   : a surface mesh triangle list (ne x 3)
%   opt: a list of optional parameters, currently surfacenorm supports:
%        'Type': [3|4|int] if set to 3 (default), output triangular mesh; if set to 4,
%                output a quad mesh where fc is an Nx4 array; otherwise,
%                output a quad mesh where fc is a cell array
%        'Nodup': [0|1] if set to 0 (default), no duplicated nodes in x/y/z are removed;
%                if set to 1, duplicated nodes are removed
%
% output:
%   no: output surface node coordinates
%   fc: output surface connections - if 'type' set to 3 or 4, fc is a numerical array,
%       otherwise, fc is a cell array with each element containing 4 integers
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

opt=varargin2struct(varargin{:});

s=size(x);
if(~all(s==size(y) & s==size(z)))
    error('x/y/z must be a 2D array of the same size');
end

no=[x(:),y(:),z(:)];
nodelen=size(no,1);
fc=zeros((s(1)-1)*s(2),4);
row=[(1:s(1)-1)', (2:s(1))', (s(1)+2:2*s(1))', (s(1)+1:2*s(1)-1)'];

for i=0:s(2)-1
    fc(i*(s(1)-1)+1:(i+1)*(s(1)-1),:)=row+i*s(1);
end

fc(fc>nodelen)=fc(fc>nodelen)-nodelen;

if(jsonopt('nodup',0,opt))
    [no,fc]=removedupnodes(no,fc);
end

outputtype=jsonopt('type',3,opt);
if(outputtype==3)
    fc=[fc(:,[1 2 3]); fc(:,[1 3 4])];
elseif(outputtype==4)
    fc=num2cell(fc,2);
end

author: fangq


text:
function volCoords = functional2volXformCoords(inplane, gray, preserveExactValues)
% Return coords from Volume view that correspond to coords in
% mrVista functional inplane coords
%
%  volCoords = functional2volXformCoords(inplane, gray, preserveExactValues)
%
% We first get a 3xn matrix of functional view coordinates, then find the
% corresponding inplane anatomical coordinates, and finally find the
% corresponding volume  coordinates. We use this transform when
% we want to convert functional data (e.g., a parameter map, coranal, or
% time series) from the volume view to the inplane view.
%
% INPUTS
%   inplane: mrVista view structure (must be an inplane view)
%   gray: mrVista view structure (must be a gray view)
%   preserveExactValues: boolean. If false, return integer coordinates. If
%                   true, return the calculated (non-integer values). If
%                   non-integer values are returned, then the parent
%                   function will have to deal with these, e.g., via
%                   interpolation.
% OUTPUTS
%   volCoords: 3xn matrix of coordinates in Volume space
%                   corresponding to 3xn matrix of inplane functional  coords
%
% Example:
%   volCoords = ip2volXformCoords(inplane, gray)
%
%
% CO & JW 2016.01.14

% Don't do this unless inplane is really an inplane and volume is really a volume
if ~strcmp(viewGet(inplane, 'viewType'),'Inplane')
    myErrorDlg('ip2volParMap can only be used to transform from inplane to volume/gray.');
end
if ~strcmp(viewGet(gray, 'viewType'),'Volume') &&~strcmp(viewGet(gray, 'viewType'),'Gray')
    myErrorDlg('ip2volParMap can only be used to transform from inplane to volume/gray.');
end

% check inputs
if ~exist('preserveExactValues', 'var'), preserveExactValues = false; end

% we need mrSESSION for the alignment matrix
mrGlobals;

% The gray coords are the integer-valued (y,x,z) volume
% coordinates that correspond to the inplanes.  Convert to
% homogeneous form by adding a row of ones.
preserveCoords = true;
functionalCoords  = ip2functionalCoords(inplane,    viewGet(inplane, 'coords'), ...
    [], preserveCoords, preserveExactValues);
nVoxels = size(functionalCoords, 2);
functionalCoords  = double([functionalCoords; ones(1,nVoxels)]);

% inplane2VolXform is the 4x4 homogeneous transform matrix that
% takes inplane (y',x',z',1) coordinates into Volume (y,x,z,1)
% coordinates.
inplane2VolXform = sessionGet(mrSESSION,'alignment');

% We don't care about the last coordinate in (y,x,z,1), so we
% toss the fourth row of Xform.  Then our outputs will be (y,x,z).
%
inplane2VolXform = inplane2VolXform(1:3,:);

% Transform coord positions to the volume.  Hence, grayCoords
% contains the volume position of each of the inplane  voxels.  These
% will generally not fall on integer-valued coordinates, rather they will
% fall between voxels.  Other functions that rely on the output of this
% function will require interpolation to get the data at these
% between-voxel positions.
%
volCoords = inplane2VolXform*functionalCoords;

% % Convert coords from inplane anatomical space to inplane functional space.
% % We do this because the anatomical inplane is often higher resolution than
% % the functional data.  We preserve the number of coords so that the number
% % of output functional voxels is identical to the number of gray or volume
% % voxels. If requested, we preserve the exact (non-integer) values, which
% % means that the functional coordinates will lie in locations between the
% % actual functional data points.
% preserveCoords = true;
% volCoords   = ip2functionalCoords(inplane, ipAnatomicalCoords, ...
%     [], preserveCoords, preserveExactValues);

% Some confusion about zeros in the output coordinates. It seems best to
% leave the calculations as they are unless there is a problem that someone
% understands. See below.

% HH, 6/16/2012: If we use the code below, we elimintate some voxels
% and then nVoxels is different from the size of ipCoords, leading to a
% mismatch in the size of the expected time series (tSeries) and the
% size of interpolated tSeries, causing an error in  the line,
%  tSeries(frame,:) = interp3(subData, ...
%                                   ipFuncCoords(2,:), ...
%                                   ipFuncCoords(1,:), ...
%                                   ipFuncCoords(3,:), ...
%                                   method);
%
% Hence we comment out the lines below. Leaving in the values with 0
% did not produce an error, at least for the function we tested it on,
% ip2volTSeries.
%
% % ras 12/20/04:
% % occasionally the xformed grayCoords will include a 0
% % coordinate. While it seems this should not happen,
% % I'm applying this band-aid for the time being:
% [badRows badCols] = find(ipFuncCoords==0); %#ok<ASGLU>
% goodCols = setdiff(1:nVoxels, badCols);
% ipFuncCoords = ipFuncCoords(:,goodCols);
% if length(badCols)>1
%     fprintf('%i voxels mapped to slice 0...\n',length(badCols));
% end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [tt]=tt_qutrtoepl(x)

% returns the multilevel upper-triangular Toeplitz matrix tt generated by the multi-dimensional input vector x
% in the QTT format
%
% If the size of the input vector is N,
% then the size of the output matrix is N x N.
%
% April 20, 2011
% Vladimir Kazeev
% vladimir.kazeev@gmail.com
% INM RAS
% Moscow, Russia
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For details please see the preprint
% http://www.mis.mpg.de/publications/preprints/2011/prepr2011-36.html
% Vladimir A. Kazeev, Boris N. Khoromskij and Eugene E. Tyrtyshnikov
% Multilevel Toeplitz matrices generated by QTT tensor-structured vectors and convolution with logarithmic complexity
% January 12, 2012
% Vladimir Kazeev,
% Seminar for Applied Mathematics, ETH Zurich
% vladimir.kazeev@sam.math.ethz.ch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

d=size(x,1);
tt=tt_qshiftstack_r(d);
tt=tt_qreshape(tt,3,[4*ones(d,1),2*ones(d,1)]);
tt=tt_mv(tt,x);
tt=tt_qreshape(tt,1,2*ones(d,2));

return
end

author: oseledets

repo: TT-Toolbox


text:
function [midpoint,minDist,t]=closestPointBetween2Lines(a1,b1,a2,b2)
%%CLOSESTPOINTBETWEEN2LINES Given two non-parallel lines in numDim
%              dimensions that do not necessarily intersect, find the point
%              closest to both lines. This is done by finding the point on
%              each line that is closest to the other line, and then
%              averaging the points (get the point halfway between the
%              lines).
%
%INPUTS: a1,b1 numDimXN matrix of parameters for the first line in N pairs
%              of lines, the closes point between each being desired. The
%              ith line is represented in parametic form as
%              r=a1(:,i)*t+b1(:,i), where t is a scalar parameter and p is
%              a point on the line.
%        a2,b2 The numDimXN set of parameters for the second lines.
%
%OUTPUTS: midpoint The numDimXN set of closest points between each of the N
%                  pairs of lines.
%          minDist The NX1 set of minimum distances between the pairs of
%                  lines.
%                t The 2XN set of values of t used in the linear equations
%                  to get the minimum points. These values can be positive
%                  or negative.
%
%Line 1 is r1=a1*t1+b1 and line 2 is r2=a2*t2+b2. We want to find the
%values of t1 and t2 such that norm(r1-r2)^2 is minimized.
%norm(r1-r2)^2=norm(a1*t1+b1-a2*t2-b2)^2
%             =norm(A*T+c)^2 where A=[a1,-a2] and c=b1-b2 and T=[t1;t2]
%             =T'*A'*A*T+2*c'*A*T+c'*c
%Setting the derivative with respect to the vector T equal to zero, we get
%             2*A'*A*T+2*A'*c=0
%The solution is T=-inv(A'*A)*A'*c, which (being the zero deivative point)
%is an extrema and being the only extrema (unless the lines are parallel)
%is the closest point between the lines. The pinv function is used to
%evaluate (A'*A)\A' as it is more stable.
%
%EXAMPLE 1:
%This is a simple 2D case, so the lines intersect. The starting points are
%drawn in blue and lines with positive t are drawn in blue. The
%intersection point is marked in red.
% x1=[4;3];
% u1=[3;2];
% u1=u1/norm(u1);
% x2=[6;7];
% u2=[6;-1];
% u2=u2/norm(u2);
% r=10;
% p1=[x1,x1+r*u1];
% p2=[x2,x2+r*u2];
% figure(1)
% clf
% hold on
% scatter(x1(1),x1(2),400,'.b')
% scatter(x2(1),x2(2),400,'.b')
% plot(p1(1,:),p1(2,:),'-b')
% plot(p2(1,:),p2(2,:),'-b')
% pt=closestPointBetween2Lines(u1,x1,u2,x2);
% scatter(pt(1),pt(2),300,'.r')
%
%EXAMPLE 2:
%This is the same as example 1, but is in 3D. This time, the lines do not
%intersect in 3D, so we use t to compute the two nearest points on either
%line and we draw a red line between them.
% x1=[4;3;8];
% u1=[3;2;0];
% u1=u1/norm(u1);
% x2=[6;7;-2];
% u2=[6;-1;4];
% u2=u2/norm(u2);
%
% r=20;
% p1=[x1,x1+r*u1];
% p2=[x2,x2+r*u2];
% figure(1)
% clf
% hold on
% scatter3(x1(1),x1(2),x1(3),400,'.b')
% scatter3(x2(1),x2(2),x2(3),400,'.b')
% plot3(p1(1,:),p1(2,:),p1(3,:),'-b')
% plot3(p2(1,:),p2(2,:),p2(3,:),'-b')
% [p,~,t]=closestPointBetween2Lines(u1,x1,u2,x2);
% scatter3(p(1),p(2),p(3),200,'.r')
% pt1=x1+t(1)*u1;
% pt2=x2+t(2)*u2;
% plot3([pt1(1);pt2(1)],[pt1(2);pt2(2)],[pt1(3);pt2(3)],'-r')
% view(8,90)
%
%March 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

N=size(a1,2);
numDim=size(a1,1);

minDist=zeros(N,1);
midpoint=zeros(numDim,N);
t=zeros(2,N);
for i=1:N
    A=[a1(:,i),-a2(:,i)];
    c=b1(:,i)-b2(:,i);

    %pinv(A)=(A'*A)\A', we are doing pinv(A)*c using lsqminnorm.
    TMin=-lsqminnorm(A,c);

    t(:,i)=TMin;

    r1Min=a1*TMin(1)+b1;
    r2Min=a2*TMin(2)+b2;

    minDist(i)=norm(r1Min-r2Min);
    midpoint(:,i)=(r1Min+r2Min)/2;
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "Hill+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: updates the spring attributes!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function springs_info = update_Springs(dt,current_time,xLag,yLag,springs_info)

%springs_info: col 1: starting spring pt (by lag. discretization)
%              col 2: ending spring pt. (by lag. discretization)
%              col 3: spring stiffness
%              col 4: spring resting lengths

%RL = springs_info(:,4); % resting-length vector

dist = 1.250077e-01; %coming from Jellyfish.m (distance between ends of jellyfish bell)

% CHANGE RESTING LENGTH BTWN SIDES OF JELLYFISH BELL
springs_info(end,4) = dist*( 1 - 0.5*abs( sin(10*pi*current_time) ) );

author: nickabattista


text:
function [x,P]=onePointCartInit(zCart,SRCart,higherDerivStdDev,matType)
%%ONEPOINTCARTINIT This function implements single-point initialization for
%              target states that consist of components of position,
%              velocity, acceleration, etc. This function initializes
%              tracks from a single measurement by setting the position
%              components to Cartesian measurement value with its
%              associated covariance and then setting the diagonal elements
%              of the rest of the components based on fixed standard
%              deviations. For example, in [1] it is suggested to make the
%              standard deviation for velocity vMax/3 and in Chapter 3.2.2
%              of [2], vMax/2. Similar ad-hoc values could be used for
%              higher moments. This function does not use Doppler/ range
%              rate.
%
%INPUTS: zCart A zDimXnumMeas set of Cartesian measurements for which
%              single-point differencing should be used to start tracks.
%       SRCart If matType is omitted or is 0, then this is a
%              zDimXzDimXnumMeas set of lower-triangular square root
%              covariance matrices associated with the measurements in
%              zCart. If all of the matrices are the same, then a single
%              zDimXzDim matrix can be passed. If matType=1, then this is a
%              set of covariance matrices.
% higherDerivStdDev A numMomentsX1 or 1XnumMoments vector containing the
%              standard deviations to use for each of the moments
%              (position, velocity, etc) that cannot be estimated from the
%              data. As mentioned in [1] and in and in Chapter 3.2.2
%              of [2], for velocity, this might be vMax/sqrt(2) or
%              vMax/sqrt(3).
%      matType An optional input specifying whether SRCart is a set of
%              lower-triangular square roots of the covariance matrix, or
%              whether it is the set of covariance matrices. Possible
%              values are:
%              0 (The default if omitted or an empty matrix is passed)
%                SRCart holds lower-triangular square root covariance
%                matrices.
%              1 SRCart holds covariance matrices.
%
%OUTPUTS: x The xDimXnumMeas set of target state estimates. All
%           non-position components are zero. xDim=zDim*(numMoments+1). The
%           components are arranged position, velocity, acceleration, etc.
%           For example, [x;y;z;xDot;yDot;zDot].
%         P The xDimXxDimXnumMeas set of initial target state covariance
%           matrices associated with x.
%
%One-point differencing is discussed in [1] and Chapter 3.2.2 of [2].
%
%REFERENCES:
%[1] M. Mallick and B. La Scala, "Comparison of single-point and two-point
%    difference track initiation algorithms using position measurements".
%    Acta Automatica Sinica, vol.34, no. 3, pp 258-265, Mar. 2008.
%[2] Y. Bar-Shalom, P. K. Willett, and X. Tian, Tracking and Data Fusion.
%    Storrs, CT: YBS Publishing, 2011.
%
%November 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(matType))
    matType=0;
end

zDim=size(zCart,1);
numMeas=size(zCart,2);

if(size(SRCart,3)==1)
    SRCart=repmat(SRCart,1,1,numMeas);
end

numMoments=length(higherDerivStdDev);
xDim=zDim*(numMoments+1);

x=zeros(xDim,numMeas);
P=zeros(xDim,xDim,numMeas);

x(1:zDim,:)=zCart;
switch(matType)
    case 0
        for curMeas=1:numMeas
            P(1:zDim,1:zDim,curMeas)=SRCart(:,:,curMeas)*SRCart(:,:,curMeas)';
        end
    case 1
        for curMeas=1:numMeas
            P(1:zDim,1:zDim,curMeas)=SRCart(:,:,curMeas);
        end
    otherwise
        error('Unknown matrix type specified.')
end

sel=(zDim+1):xDim;
P(sel,sel,:)=repmat(kron(diag(higherDerivStdDev(:).^2),eye(zDim,zDim)),1,1,numMeas);
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
%  Examplescript for usage of omcompile, omrun, omparameter, omimport
%  and omcleanup
%  For more help, look into the help of the functions and the functions
%  themselves. E.g. it's important to adapt the path to omc.exe in
%  omcompile.m. You could maybe improve performance by changing some value in omimport.m
%
%  Feedback/problems: Christian Schaad, ingenieurbuero@christian-schaad.de

clear all; close all;hold on;

modelname='HelloWorld.test';   % Set Modelname

omcompile('hello.mo',modelname);  % Load hello.mo, compile and run HelloWorld.test, do NOT load Modelica library (0), no debug (0)

omcleanup(modelname)  % Delete the compiler logs, makefile, etc.

%omparameter(modelname,'stop value',2.5);  % Change "stop value" in parameter file

%omparameter(modelname,'step value',1e-02);% Change "step value" in parameter file

j=1;
for value=2:1.5:10;   % Loop over values of parameter a

omparameter(modelname,'Hello1.a',num2str(value));  % Change parameter Hello1.a in parameter file to value

omrun(modelname)  % Run the compiled model

omimport(modelname);  % Read results into Matlab workspace

plot(time,Hello1.x)   % plot, legend, save or whatever...
hold all;
legende(j).name=['Hello1.a=',num2str(value)];
j=j+1;
end;

legend(legende.name)

clear legende i j
grid on;

% Use omparvar.m for automatical creation of plots with different
% parameters...


text:
% Image construction from overlapping patches
function [result] = overlap_add(patches, img_size, grid)

result = zeros(img_size);
weight = zeros(img_size);

for i = 1:size(grid, 3)
    patch = reshape(patches(:, i), size(grid, 1), size(grid, 2));
    result(grid(:, :, i)) = result(grid(:, :, i)) + patch;
    weight(grid(:, :, i)) = weight(grid(:, :, i)) + 1;
end

I = logical(weight);
result(I) = result(I) ./ weight(I);

author: thomas-koehler

repo: SupER


text:
% ============================================================= %
% Files of the Matlab programs included in the book:            %
%     Xin-She Yang, Nature-Inspired Metaheuristic Algorithms,   %
%     Second Edition, Luniver Press, (2010).   www.luniver.com  %
% ============================================================= %
% The Accelerated Particle Swarm Optimization:                  %
% For detail, section Chapter 5 (section 5.3) in the book by    %
% Xin-She Yang, Nature-Inspired Metaheuristic Algorithms,       %
% First Edition, Luniver Press, (2008). Also in Chapter 8       %
% (section 8.3, page 65) of the second edition (2010).          %
% ------------------------------------------------------------- %

% =============================================================
% Usage: apso(number_of_particles,timesteps)
%  eg:   best=apso(20,10);
%  where best=[xbest ybest zbest]  %a n by 3 matrix
% -------------------------------------------------------------
function [best]=apso(n,timesteps)
% n=number of particles
% timesteps=total number of time steps
if nargin<2,   timesteps=15;  end
if nargin<1,   n=20;         end
% -------------------------------------------------------------
% Michalewicz Function f*=-1.8013 at [2.20319, 1.57049]
funstr='-sin(x)*(sin(x^2/pi))^20-sin(y)*(sin(2*y^2/pi))^20';
f=inline(vectorize(funstr));
% Lb=[xmin  ymin] and Ub=[xmax ymax];
Lb=[0 0]; Ub=[4 4];
% -------------------------------------------------------------
% Setting  parameters alpha, beta
% Randomization amplitude of roaming particles alpha=[0,1]
% alpha=gamma^t=0.7^t;
% Speed of convergence (0->1)=(slow->fast); % beta=0.5
beta=0.5;
% -------------------------------------------------------------
% Grid values for show the topology of the function
% to be optimized. For display only, not used for optimization.
Ngrid=100; % Ngrid=100;
dx=(Ub(1)-Lb(1))/Ngrid;
dy=(Ub(2)-Lb(2))/Ngrid;

[x,y]=meshgrid(Lb(1):dx:Ub(1),Lb(2):dy:Ub(2));
% get the shape of the function
% you can define any function to be optimized in 'fun(x,y)'
% funflag=1 or 2 or 3 choose a function from a list

z=f(x,y);

% Display the shape of the function to be optimized
figure(1);
surfc(x,y,z);
% -----------------------------------------------------------
best=zeros(timesteps,3);

% ===========================================================
% ------------- Start Particle Swarm Optimization -----------
% generating the initial locations of n particles
[xn,yn]=init_pso(n,Lb,Ub);
% Display the paths of particles in a figure with
% contours of the function to be optimized
   figure(2);

% Time-marching or iterations
for i=1:timesteps,     %%%%% start time-stepping
   % Show the contour of the function
  contour(x,y,z,15); hold on;
  colormap cool;

% Find the current best location (xo,yo)
zn=f(xn,yn);
zn_min=min(zn);
xo=min(xn(zn==zn_min));
yo=min(yn(zn==zn_min));
zo=min(zn(zn==zn_min));

% Trace the paths of all roaming particles
% Display these roaming particles
plot(xn,yn,'.',xo,yo,'d','markersize',10,'markerfacecolor','g');

% The accelerated PSO alpha=alpha_0* gamma^t with alpha_0=1
  gamma=0.7; alpha=gamma.^i;
% Move all particle to new locations
  [xn,yn]=pso_moving(xn,yn,xo,yo,alpha,beta,Lb,Ub);
  drawnow;

% If you keep "hold on", then the paths of particles are displayed
	hold off;
% History
   best(i,1)=xo; best(i,2)=yo; best(i,3)=zo;
% Output the results (xo,yo) to screen
   str=strcat('Best estimates: =',num2str([xo yo]));
   str=strcat(str,'  fmin='); str=strcat(str,num2str(best(i,3)));
   str=strcat(str,'  Iteration='); str=strcat(str,num2str(i));
   disp(str);
end   %%%%% end of time-stepping

% Postprocessing & display best
bestsolution=best(end,1:2)
fmin=best(end,3)
% -----------------------------------------------------------------
% All subfunctions are listed here
% -----------------------------------------------------------------
% Intial locations of particles
function [xn,yn]=init_pso(n,Lb,Ub)
xn=rand(1,n)*(Ub(1)-Lb(1))+Lb(1);
yn=rand(1,n)*(Ub(2)-Lb(1))+Lb(2);

% Move all the particles toward (xo,yo)
function [xn,yn]=pso_moving(xn,yn,xo,yo,alpha,beta,Lb,Ub)
% Get the number of particles
nn=size(yn,2);
xn=xn.*(1-beta)+xo.*beta+alpha.*randn(1,nn);
yn=yn.*(1-beta)+yo.*beta+alpha.*randn(1,nn);
[xn,yn]=findrange(xn,yn,Lb,Ub);

% Make sure the particles are not outside of the range
function [xn,yn]=findrange(xn,yn,Lb,Ub)
nn=length(yn);
for i=1:nn,
   if xn(i)<=Lb(1), xn(i)=Lb(1); end
   if xn(i)>=Ub(1), xn(i)=Ub(1); end
   if yn(i)<=Lb(2), yn(i)=Lb(2); end
   if yn(i)>=Ub(2), yn(i)=Ub(2); end
end
% End of the APSO -------------------------------------------------


text:
function params = retinoCheckParams(view, dt, scans);
% Check if visual field mapping ("retinotopy", though it doesn't
% need to be retinotopic) parameters have been set for the specified
% scans, giving the user an option to set them if they haven't been
% assigned yet.
%
% params = retinoCheckParams(view, <dt, scans>);
%
% Returns the parameters as a struct array of length(nScans).
% Will prompt the user to set any scans that haven't been set yet
% in the scan list.
%
% ras, 01/06: testing the waters if this code is needed. I see many
% other places where similar parameters are set, but none of them
% seem immediately useable to me.
if notDefined('view'), view = getCurView;                   end
if notDefined('dt'), dt = viewGet(view, 'curdt');           end
if notDefined('scans'), scans = viewGet(view, 'curscan');   end

mrGlobals;
params = [];

if isnumeric(dt)
    dtNum = dt;
    dt = dataTYPES(dtNum).name;
else
    dtNum = existDataType(dt);
end

for i = 1:length(scans)
    check = retinoGetParams(view, dt, scans(i));

    if isempty(check)
        % not assigned: ask user if he/she wants to assign now
        annotation = dataTYPES(dtNum).scanParams(scans(i)).annotation;
        q = [sprintf('%s, scan %i (%s): ', dt, scans(i), annotation) ...
            'This scan doesn''t have any visual field-mapping ' ...
            'parameters assigned yet. Do you want to assign them now?'];
        resp = questdlg(q, mfilename, 'Yes', 'No', 'Yes');

        if isequal(resp, 'Yes')
            check = retinoSetParams(view, dt, scans(i));
        else
            error('No visual field map parameters set -- user aborted.');
        end
    end

    % at this point, check should be assigned, or we would have aborted
    if isempty(params)  % first scan
        params = check;
    else
        % add the current params to the struct array: copy a field
        % at a time, since the params for different scans may have
        % different fields (different design types)
        for f = fieldnames(check)'
            params(i).(f{1}) = check.(f{1});
        end
    end
end

return


text:
function idx = getNearest(data, centers,trans)
% idx = getNearest(data, centers)
%
% Gets indices of centers closest (in Euclidean space) to each data point
%
% data(ndata, nvars)
% centers(ncenters, nvars)
% idx(ndata, 1)

% dist(a, b) = ||a-b|| = sum(a.^2) - 2*a*b' + sum(b.^2) for 1xnvars vectors
% a and b.

if(~exist('trans','var'))
    centers = centers';
end
centerssq = sum(centers.^2, 1);

distmat = repmat(centerssq, [size(data,1) 1]);
distmat = distmat - 2*data * centers;

[dist, idx] = min(distmat, [], 2);

author: Cloud-CV

repo: object-proposals


text:
function bpab_test ( )

%*****************************************************************************80
%
%% BPAB_TEST tests BPAB.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BPAB_TEST\n' );
  fprintf ( 1, '  BPAB evaluates Bernstein polynomials.\n' );
  fprintf ( 1, '\n' );

  x = 0.3;
  a = 0.0;
  b = 1.0;
  bern = bpab ( n, x, a, b );

  fprintf ( 1, '  The Bernstein polynomials of degree %d\n', n );
  fprintf ( 1, '  based on the interval from %f\n', a );
  fprintf ( 1, '  to %f\n', b );
  fprintf ( 1, '  evaluated at X = %f\n', x );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   I        Bern(I,X)\n' );
  fprintf ( 1, '\n' );

  for i = 0 : n
    fprintf ( 1, '  %2d  %12f\n', i, bern(i+1) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_failed = test_libltfat_normalize(varargin)
test_failed = 0;

fprintf(' ===============  %s ================ \n',upper(mfilename));

definput.flags.complexity={'double','single'};
[flags]=ltfatarghelper({},definput,varargin);
dataPtr = [flags.complexity, 'Ptr'];

[~,~,enuminfo]=libltfatprotofile;
Cenumnorms = enuminfo.ltfat_norm_t;

Larr = [1,9,11,110];
normpairs = {{'null',Cenumnorms.LTFAT_NORM_NULL},...
             {'1', Cenumnorms.LTFAT_NORM_1},...
             {'2', Cenumnorms.LTFAT_NORM_2}};

for do_complex = 0:1
    complexstring = '';
    if do_complex, complexstring = 'complex'; end

    funname = makelibraryname('normalize',flags.complexity,do_complex);

    for L = Larr
        for npId = 1:numel(normpairs)
            normpair = normpairs{npId};

            if do_complex
                z = cast(randn(L,1)+1i*randn(L,1),flags.complexity);
                zi = complex2interleaved(z);
                zout = randn(size(zi),flags.complexity);

                ziPtr = libpointer(dataPtr,zi);
                zoutPtr = libpointer(dataPtr,zout);
            else
                z = cast(randn(L,1),flags.complexity);
                zi = z;
                zout = randn(size(zi),flags.complexity);

                ziPtr = libpointer(dataPtr,zi);
                zoutPtr = libpointer(dataPtr,zout);
            end

            trueres = normalize(z,normpair{1});

            status = calllib('libltfat',funname, ziPtr,L,normpair{2},zoutPtr);

            if do_complex
                res = norm(trueres - interleaved2complex(zoutPtr.Value));
            else
                res = norm(trueres - zoutPtr.Value);
            end

            [test_failed,fail]=ltfatdiditfail(res-status,test_failed);
            fprintf(['NORMALIZE OP L:%3i, norm: %s %s %s %s %s\n'],L,normpair{1},flags.complexity,complexstring,ltfatstatusstring(status),fail);

            status = calllib('libltfat',funname, ziPtr,L,normpair{2},ziPtr);

            if do_complex
                res = norm(trueres - interleaved2complex(ziPtr.Value));
            else
                res = norm(trueres - ziPtr.Value);
            end

            [test_failed,fail]=ltfatdiditfail(res-status,test_failed);
            fprintf(['NORMALIZE IP L:%3i, norm: %s %s %s %s %s\n'],L,normpair{1},flags.complexity,complexstring,ltfatstatusstring(status),fail);
        end
    end
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%########################################################################
%
%	- PPGI Toolbox -
%   A MATLAB toolbox for Photoplethysmography Imaging (PPGI)
%
% Author   : Christian S. Pilz
% Company  : The Nature of Space of Time
% Date     : 07.05.2019
%
% Contact  : cpi@partofthestars.com
% Web Page : www.partofthestars.com
%
% Version  : beta0.1
%
%########################################################################
%
%	test_stochastic_resonator.m:
%
% Description:
%
% 	The frequency trajectories are drawn randomly such that the
% 	random walk (Wiener process) defines a trajectory that is
% 	transformed into a frequency time series f (in Hertz). The oscillator state space
% 	method is applied to the data after the simulated signal has been
% 	transformed to downsampled observations (defined by dts). The mean
% 	squared error results are then captured and visualized.
%

%
close all; %% Simulate M draws

% Number of random draws
M = 1;

% Different time discretizations (TR) to consider
dts = [0.01 0.05 0.1:0.1:1 1.2:.2:2.4];
% Allocate space for results
MSE = zeros(M,numel(dts));
C = zeros(M,numel(dts));

for i=1:M
	[MSE(i,:),C(i,:),x,f] = simulate_and_estimate(dts,25);
	figure;
	subplot(2,1,1)
	plot(x(1,:),'black');
	title('Stochastic Oscillator');
	ylabel('Amplitude');
	xlabel('Time in seconds')
	x_ticks = 0:0.01:25;
	set(gca,'XTickLabel',x_ticks(1:2500).*500 );
	subplot(2,1,2)
	plot(f,'black')
	ylabel('Frequency in Hz');
	xlabel('Time in frames')
	title('Oscillator Frequency')
end

author: partofthestars

repo: PPGI-Toolbox


text:
function matrixPSD = enforcePSD(matrix)
    if size(matrix,1) ~= size(matrix,2)
        error('Input matrix is not symmetric.');
    else
        matrixPSD = matrix;
        for r = 1:size(matrixPSD,1)
            for c = 1:size(matrixPSD,2)
                if r == c
                    matrixPSD(r,c) = abs(matrixPSD(r,c));
                else
                     offDiagElement = mean([matrixPSD(r,c),matrixPSD(c,r)]);
%                      offDiagElement = matrixPSD(r,c);
                     matrixPSD(c,r) = offDiagElement;
                     matrixPSD(r,c) = offDiagElement;
                end
            end
        end
    end
end

author: utiasSTARS

repo: msckf-swf-comparison


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A heuristic way to further refine the output windows
%
% For each small output window, we run our method on the
% this ROI again and extract the output that has the
% largest IOU with the orignal window for replacement.
%
% NMS is further applied to remove duplicate windows.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function res = refineWin(I, res, net, param)

imsz = [size(I,1) size(I,2)];
param.lambda = 0.05;
for i = 1:size(res,2)
    bb = res(:,i);
    bbArea = (bb(3)-bb(1))*(bb(4)-bb(2));
    % only refine small windows
    if bbArea < 0.125*imsz(1)*imsz(2)
        margin = (bb(3)-bb(1)+bb(4)-bb(2))*0.2;
        bb = round(expandROI(bb, imsz, margin));
        Itmp = I(bb(2):bb(4), bb(1):bb(3) , :);
        [Ptmp, Stmp] = getProposals(Itmp, net, param);
        restmp = propOpt(Ptmp, Stmp, param);
        if ~isempty(restmp)
            restmp = getROIBBox(restmp, bb);
            [~,ii] = max(getIOUFloat(restmp', res(:,i)));
            res(:,i) = restmp(:,ii);
        end
    end
end
res = doNMS(res, 0.5);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function value = r8mat_norm_fro ( m, n, a )

%*****************************************************************************80
%
%% R8MAT_NORM_FRO returns the Frobenius norm of an R8MAT.
%
%  Discussion:
%
%    The Frobenius norm is defined as
%
%      value = sqrt ( sum ( 1 <= I <= M ) sum ( 1 <= j <= N ) A(I,J)**2 )
%
%    The matrix Frobenius norm is not derived from a vector norm, but
%    is compatible with the vector L2 norm, so that:
%
%      vec_norm_l2 ( A * x ) <= mat_norm_fro ( A ) * vec_norm_l2 ( x ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in A.
%
%    Input, integer N, the number of columns in A.
%
%    Input, real A(M,N), the matrix whose Frobenius
%    norm is desired.
%
%    Output, real VALUE, the Frobenius norm of A.
%
  value = sqrt ( sum ( sum ( a(1:m,1:n).^2 ) ) );

  return
end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function [RCM1, RCM2, RCM3, RCM4, RCM5] = region_covariance(img, options)
% Region covariance matrix (GCM) algorithm
%
% Inputs:
%       img                 image data
%       options             options
% Output:
%       RCM1                region covariance matrix of size 7x7 (RCM1)
%                           pixel location, intensity, 1st&2nd-order gradient
%       RCM2                region covariance matrix of size 7x7 (RCM2)
%                           pixel location, 1st&2nd-order gradient, and edge orientation
%       RCM3                region covariance matrix of size 6x6 (RCM3)
%                           pixel location, and 1st&2nd-order gradient
%       RCM4                region covariance matrix of size 8x8 (RCM4) with
%                           pixel location, intensity, 1st&2nd-order gradient, and
%                           edge orientation
%       RCM5                region covariance matrix of size 5x5 (RCM5) with
%                           intensity, 1st&2nd-order gradient (not appear in the paper)
%
%
% References:
%       Yanwei Pang, Yuan Yuan, and Xuelong Li
%       'Gabor-Based Region Covariance Matrices for Face Recognition'
%       IEEE Transactions on Circuits and Systems for Video Technology vol.18, no.7, 2008.
%
%
% Created by H.Kasai on June 23, 2017

    symm = @(X) .5*(X+X');

    % check correct number of arguments
    if (nargin ~= 2)
        error('Please use the correct number of input arguments!')
    end

    % extract options
    if ~isfield(options, 'spd_projection')
        spd_projection = true;
    else
        spd_projection = options.spd_projection;
    end

    % check if the input image is grayscale
    if size(img,3) == 3
        warning('The input RGB image is converted to grayscale!')
        img = rgb2gray(img);
    end

    % obtain image size
    [ysize, xsize] = size(img);
    num_of_pixels = ysize * xsize;

    % prepare array
    Feature_Mat = zeros(8, num_of_pixels);

    % calculate first-order and second-order derivatives
    img = double(img);
    [Ix, Iy] = gradient(img);   % first order partials
    [Ixx, Ixy] = gradient(Ix);  % second order partials
    [Iyx, Iyy] = gradient(Iy);  % second order partials

    % calculate 9-dimensional feature vectior for every pixels
    pixel_cnt = 0;
    for y = 1 : ysize
        for x = 1 : xsize
            %fprintf('[%d %d] %d)\n', y, x, pixel_cnt);
            pixel_cnt = pixel_cnt + 1;

            Feature_Mat(1, pixel_cnt)   = x;
            Feature_Mat(2, pixel_cnt)   = y;
            Feature_Mat(3, pixel_cnt)   = img(y,x);

            Feature_Mat(4, pixel_cnt)   = abs(Ix(y,x));
            Feature_Mat(5, pixel_cnt)   = abs(Iy(y,x));

            Feature_Mat(6, pixel_cnt)   = abs(Ixx(y,x));
            Feature_Mat(7, pixel_cnt)   = abs(Iyy(y,x));

            theta = atan(abs(Iy(y,x))/abs(Ix(y,x)));
            Feature_Mat(8, pixel_cnt)   = theta;

        end
    end

    % calculate covariance matrix
    Feature_Mat_type1 = Feature_Mat;
    Feature_Mat_type1(8,:) = [];            % remove orientation (theta)
    RCM1 = cov(Feature_Mat_type1');         % case for Eq.(5)

    Feature_Mat_type2 = Feature_Mat;
    Feature_Mat_type2(3,:) = [];            % remove intensity
    RCM2 = cov(Feature_Mat_type2');         % case for Eq.(6)

    Feature_Mat_type3 = Feature_Mat;
    Feature_Mat_type3(8,:) = [];            % remove edge orientation (theta)
    Feature_Mat_type3(3,:) = [];            % remove intensity
    RCM3 = cov(Feature_Mat_type3');         % case for Eq.(8)

    RCM4 = cov(Feature_Mat');               % case for Eq.(9)

    Feature_Mat_type5 = Feature_Mat;
    Feature_Mat_type5(8,:) = [];            % remove edge orientation (theta)
    Feature_Mat_type5(1:2,:) = [];          % remove location
    RCM5 = cov(Feature_Mat_type5');         % additional case

    % project onto SPD if needed
    if spd_projection
        RCM1 = symm(spd_project(RCM1));
        RCM2 = symm(spd_project(RCM2));
        RCM3 = symm(spd_project(RCM3));
        RCM4 = symm(spd_project(RCM4));
        RCM5 = symm(spd_project(RCM5));
    end

end


text:
% uiputfile2() - same as uigputfile but remember folder location.
%
% Usage: >> uiputfile2(...)
%
% Inputs: Same as uiputfile
%
% Author: Arnaud Delorme & Hilit Serby, Scott Makeig, SCCN, UCSD, 2004
%         Thanks to input from Bas Kortmann
%
% Copyright (C) Arnaud Delorme & Hilit Serby, Scott Makeig, SCCN, UCSD, 2004

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

function varargout = uiputfile2(varargin);

    if nargin < 1
        help uiputfile2;
        return;
    end;

    % remember old folder
    %% search for the (mat) file which contains the latest used directory
    % -------------------
    olddir = pwd;
    try,
        eeglab_options;
        if option_rememberfolder
            tmp_fld = getenv('TEMP');
            if isempty(tmp_fld) & isunix
                if exist('/tmp') == 7
                    tmp_fld = '/tmp';
                end;
            end;
            if exist(fullfile(tmp_fld,'eeglab.cfg'))
                load(fullfile(tmp_fld,'eeglab.cfg'),'Path','-mat');
                s = ['cd([''' Path '''])'];
                if exist(Path) == 7, eval(s); end;
            end;
        end;
    catch, end;

    %% Show the open dialog and save the latest directory to the file
    % ---------------------------------------------------------------
    [varargout{1} varargout{2}] = uiputfile(varargin{:});
    try,
        if option_rememberfolder
            if varargout{1} ~= 0
                Path = varargout{2};
                try, save(fullfile(tmp_fld,'eeglab.cfg'),'Path','-mat','-V6'); % Matlab 7
                catch,
                    try,  save(fullfile(tmp_fld,'eeglab.cfg'),'Path','-mat');
                    catch, error('uigetfile2: save error, out of space or file permission problem');
                    end
                end
                if isunix
                    eval(['cd ' tmp_fld]);
                    system('chmod 777 eeglab.cfg');
                end
            end;
        end;
    catch, end;
    cd(olddir)

author: PatternRecognition

repo: OpenBMI


text:
function [ h ] = tapas_rdcm_get_convolution_bm(options)
% [ h ] = tapas_rdcm_get_convolution_bm(options)
%
% Creates a fixed hemodynamic response function (HRF) by convolving a single
% event (impulse) with the standard Balloon model from DCM.
%
%   Input:
%       options     - options structure with relevant information
%
%   Output:
%       h           - hemodynamic response function (HRF)
%

% ----------------------------------------------------------------------
%
% Authors: Stefan Fraessle (stefanf@biomed.ee.ethz.ch), Ekaterina I. Lomakina
%
% Copyright (C) 2016-2022 Translational Neuromodeling Unit
%                         Institute for Biomedical Engineering
%                         University of Zurich & ETH Zurich
%
% This file is part of the TAPAS rDCM Toolbox, which is released under the
% terms of the GNU General Public License (GPL), version 3.0 or later. You
% can redistribute and/or modify the code under the terms of the GPL. For
% further see COPYING or <http://www.gnu.org/licenses/>.
%
% Please note that this toolbox is in an early stage of development. Changes
% are likely to occur in future releases.
%
% ----------------------------------------------------------------------

% compile source code of integrator
tapas_rdcm_compile()

% get the DCM
DCM = options.DCM;

% number of input data points
N = size(DCM.U.u,1);

% check if full convolution should be generated
if isfield(options,'conv_full') & options.conv_full
    r_dt = 1;
else
    r_dt = N/options.conv_length;
end

% construct a dummy DCM
DCM.a    = -1;
DCM.b    = 0;
DCM.c    = 1;
DCM.d    = zeros(1,1,0);
DCM.Ep   = tapas_rdcm_empty_par(DCM);
DCM.Ep.A = DCM.a;
DCM.Ep.C = DCM.c*16;

% setting input of the fake DCM
DCM.U.u = zeros(size(DCM.U.u,1),1);
DCM.U.u(1:r_dt,:) = 1;

% setting additional parameters of the fake DCM
DCM.Y.dt = DCM.U.dt;
DCM.n    = 1;
DCM.v    = N;
DCM.ns   = N;

% create the HRF from the dummy DCM
DCM = tapas_dcm_euler_make_indices(DCM);
y   = tapas_dcm_euler_gen(DCM, DCM.Ep);

% sample the HRF at the sampling rate of the data
h = y(1:r_dt:end,1);
h = h;

end

author: translationalneuromodeling

repo: tapas


text:
% shift elements from specified dimension
% e.g. if dimIdx = 2, shift = 1, the class shifts the
% elements in dimension 2 by 1 to the right.
%
classdef ElementShiftNode < GraphNode
    properties
        dimIdx = 1;         % which dimension to select from
        shift;   % which elements to select
        circular = 0;       % whether to use circular shift
        shiftIdx;
    end
    methods
        function obj = ElementShiftNode(shift, dimIdx, circular)
            obj = obj@GraphNode('ElementShift',length(shift));
            obj.shift = shift;
            if nargin>=2
                obj.dimIdx = dimIdx;
            end
            if nargin>=3
                obj.circular = circular;
            end
        end

        function obj = forward(obj,prev_layers)
            obj = obj.preprocessingForward(prev_layers);
            input = prev_layers{1}.a;
            [D(1), D(2), D(3), D(4)] = size(input);
            nElement = D(obj.dimIdx);
            assert(nElement >= obj.shift, ...
                sprintf('%s:forward, Error: shift (%d) is larger than data length (%d), exit. \n', obj.name, obj.shift, nElement));
            obj.shiftIdx = 1:nElement;
            obj.shiftIdx = obj.shiftIdx + obj.shift;
            if obj.circular
                tmp = obj.shiftIdx>nElement;
                obj.shiftIdx(tmp) = obj.shiftIdx(tmp) - nElement;
                tmp = obj.shiftIdx<1;
                obj.shiftIdx(tmp) = obj.shiftIdx(tmp) + nElement;
            else
                obj.shiftIdx(obj.shiftIdx>nElement) = [];
                obj.shiftIdx(obj.shiftIdx<1) = [];
            end

            if obj.dimIdx==1
                obj.a = input(obj.shiftIdx,:,:,:);
            elseif obj.dimIdx==2
                obj.a = input(:,obj.shiftIdx,:,:);
            elseif obj.dimIdx==3
                obj.a = input(:,:,obj.shiftIdx,:);
            elseif obj.dimIdx==4
                obj.a = input(:,:,:,obj.shiftIdx);
            end
            obj = forward@GraphNode(obj, prev_layers);
        end

        function obj = backward(obj,future_layers, prev_layers)
            future_grad = GetFutureGrad(future_layers, obj);

            input = prev_layers{1}.a;
            [D1,D2,T,N] = size(input);

            obj.grad{1} = obj.AllocateMemoryLike([D1 D2 T N], future_grad);
            if obj.dimIdx==1
                obj.grad{1}(obj.shiftIdx,:,:,:) = future_grad;
            elseif obj.dimIdx==2
                obj.grad{1}(:,obj.shiftIdx,:,:) = future_grad;
            elseif obj.dimIdx==3
                obj.grad{1}(:,:,obj.shiftIdx,:) = future_grad;
            elseif obj.dimIdx==4
                obj.grad{1}(:,:,:,obj.shiftIdx) = future_grad;
            end
        end

    end

end

author: singaxiong

repo: SignalGraph


text:
function [ f, g, H ] = opt15_fgh ( x, flag )

%*****************************************************************************80
%
%% OPT15_FGH evaluates F, G and H for test case #15.
%
%  Discussion:
%
%    This example, if started at X = (-1,0), seems to get stuck.
%
%  Modified:
%
%    09 January 2008
%
%  Author:
%
%    Jeff Borggaard,
%    Gene Cliff,
%    Virginia Tech.
%
%  Reference:
%
%    John Dennis, Robert Schnabel,
%    Numerical Methods for Unconstrained Optimization
%    and Nonlinear Equations,
%    SIAM, 1996,
%    ISBN13: 978-0-898713-64-0,
%    LC: QA402.5.D44.
%
%  Parameters:
%
%    Input, real X(2), the evaluation point.
%
%    Input, string FLAG, indicates what must be computed.
%    'f' means only the value of F is needed,
%    'g' means only the value of G is needed,
%    'all' means F, G and H (if appropriate) are needed.
%    It is acceptable to behave as though FLAG was 'all'
%    on every call.
%
%    Output, real F, the optimization function.
%
%    Output, real G(2,1), the gradient column vector.
%
%    Output, real H(2,2), the Hessian matrix.
%
  n = length ( x );

  if ( n ~= 2 )
    fprintf ( '\n' );
    fprintf ( 'OPT15_FGH - Fatal error!\n' );
    fprintf ( '  The input vector X should have length 2.\n'),
    fprintf ( '  Instead, it has length = %d.\n', n );
    keyboard
  end

  f = - ( x(1) + x(2) ) ...
      + 0.5 * ( 1 - ( x(1)^2 + x(2)^2 ) ) ...
      + 5.0 * ( 1 - ( x(1)^2 + x(2)^2 ) )^2;

  g(1,1) = -1 - 21 * x(1) + 20 * x(1) * ( x(1)^2 + x(2)^2 );
  g(2,1) = -1 - 21 * x(2) + 20 * x(2) * ( x(1)^2 + x(2)^2 );

  H = zeros(n,n);

  H(1,1) = - 21 + 60 * x(1)^2 + 20 * x(2)^2;
  H(1,2) = 40 * x(1) * x(2);

  H(2,1) = 40 * x(1) * x(2);
  H(2,2) =  - 21 + 20 * x(1)^2 + 60 * x(2)^2;


text:
function data_out = model2data(model,fordb)
% Load complete reconstructions and output format for rBioNets ReconstructionCreator.
%
% USAGE:
%
%    data_out = model2data(model, fordb)
%
% INPUTS:
%    model:    model
%    fordb:    1 if metabolite table is wanted, default 0 if not
%
% OUTPUT:
%    data_out:
%
% .. Author: - Stefan G. Thorleifsson 2011
%
% .. rBioNet is published under GNU GENERAL PUBLIC LICENSE 3.0+
% .. Thorleifsson, S. G., Thiele, I., rBioNet: A COBRA toolbox extension for
% .. reconstructing high-quality biochemical networks, Bioinformatics, Accepted.
% .. rbionet@systemsbiology.is

if nargin < 2
    fordb = 0;
end
%   rxn table ->data
%   description
%   genes
%   model with setup change.

handles.model = model;

names = fieldnames(handles.model);

name = {'genes','data','description'};
match = zeros(size(names));
for i = 1:size(name,2)
    match(i) = strcmp(name{i},names{i});
end

if sum(match) == length(name)
    msgbox(['This file is not an complete reconstruction. Try opening'...
        ' it as a mat file.'], 'Incorrect file type','Help');
    return
end

%This field is necessary for the output.
if ~isfield(model,'grRules')
    model = creategrRulesField(model);
end

Mandatory = {'rxns','rxnNames','grRules', 'lb', 'ub','S','genes'};
%Opt numbers: CS = 1, subsystem = 2, citations = 3, comments = 4, ecNumbers = 5,
% rxnKEGGID = 6, description = 7

% reaction data lineup:
% 1.rxns
% 2.rxnNames
% 3.formula
% 4.reversible
% 5.grRules,
% 6.lb
% 7.ub
% 8.confidenceScores
% 9.subSystems
% 10.citations
% 11.comments
% 12.ecNumbers
% 13.KEGGID.

%enable is set infront afterwards.

%Check Mandatory

for i = 1:length(Mandatory)
    k = ismember(Mandatory(i),names);
    if k == 0
        msgbox(['Group ' Mandatory{i} ' is missing from structure of model.']...
            ,'Structure incorrect.','error');
        return
    end
end

S1 = size(handles.model.rxns); %Get size of model
data = cell(S1(1),13); %create cell matrix for data
data(:,1) = handles.model.rxns;
data(:,2) = handles.model.rxnNames;
data(:,3) = printRxnFormula(handles.model, handles.model.rxns,false);
data(:,4) = num2cell(double(model.lb < 0));
data(:,5) = handles.model.grRules;
data(:,6) = num2cell(handles.model.lb);
data(:,7) = num2cell(handles.model.ub);

Optional = {'confidenceScores', 'subSystems','citations', 'comments', 'ecNumbers',...
    'rxnKEGGID','description'};
%Check Optional
opt = [];
for i = 1:length(Optional)

    % Recon had different names of fields, this is a quick fix.
    if strcmp(Optional{i},'ecNumbers')
        if ismember('rxnECNumbers',names)
            Optional{i} = 'rxnECNumbers';
            continue;
        end

    elseif strcmp(Optional{i},'subSystems')
        if ismember('rxnSubSystems',names)
            Optional{i} = 'rxnSubSystems';
            continue;
        end
        if ismember('subSystems',names)
            Optional{i} = 'subSystems';
            if ~ischar([handles.model.subSystems{:}])
                handles.model.subSystems = cellfun(@(x) strjoin(x,';'),handles.model.subSystems,'UniformOutput',0);
            end
            continue;
        end
    elseif strcmp(Optional{i},'confidenceScores')
        if ismember('rxnConfidenceScores',names)
            Optional{i} = 'rxnConfidenceScores';
            continue;
        end

    elseif strcmp(Optional{i},'citations')
        if ismember('rxnReferences',names)
            Optional{i} = 'rxnReferences';
            continue;
        end

    elseif strcmp(Optional{i},'comments')
        if ismember('rxnNotes',names)
            Optional{i} = 'rxnNotes';
            continue;
        end
    end
    k = ismember(Optional(i),names);
    if k == 0
        opt = [opt,i]; %data missing
    end
end
%Opt numbers: CS = 1, subsystem = 2, citations = 3, comments = 4,
%ecNumbers = 5, rxnKEGGID = 6, description = 7
for i = 1:length(Optional)-1
    if isfield(model,Optional{i})
        if isnumeric(handles.model.(Optional{i}))
            %Confidence Scores are numeric values....
            data(:,i+7) = num2cell(handles.model.(Optional{i}));
        else
            data(:,i+7) = handles.model.(Optional{i});
        end
    else %data missing
        data(:,i+7) = cell(S1(1),1);
    end
end

%description
if ~any(ismember(opt,7)) %data exist
    if isa(handles.model.description,'char')

        s_b = [NaN NaN];
        description{1} = handles.model.description;
        description{2} = '';
        description{3} = '';
        description{4} = '';
        description{5} = '';
        description{6} = '';
        description{7} = '';
        handles.description = description;
    elseif isa(handles.model.description,'struct')
        b = fieldnames(handles.model.description);
        s_b = size(b);
    else
        s_b = [NaN NaN];
    end

    if s_b(1) == 7 %Date is not loaded into description.
        description{1} = handles.model.description.name;
        description{2} = handles.model.description.organism;
        description{3} = handles.model.description.author;
        description{4} = handles.model.description.geneindex;
        description{5} = handles.model.description.genedate;
        description{6} = handles.model.description.genesource;
        description{7} = handles.model.description.notes;
        handles.description = description;
    end

end

%This is a quick fix after adding the logical term in the front of the
%reaction table: 4. Mar. 2011
data(:,2:end+1) = data;
%Matlab does not allow data(:,1) = true;
for k=1:size(data,1)
    data{k,1} = true;
end
disabled = {};
if isfield(handles.model,'disabled')
    for i = 1:size(handles.model.disabled,1)
        disabled(i,:) = handles.model.disabled(i,:);
    end
else
    handles.model.disabled = {};
end

% Note disabled data does not carry metabolite data.

if fordb == 0
    data_out = {[data; disabled], description, handles.model.genes, handles.model};
elseif fordb == 1
    % From the metabolite database the only thing missing is the Neutral
    % formula.
     met_fields = {'mets','metNames','metFormulas','metCharges','metKEGGID',...
         'metPubChemID','metChEBIID','metInChIString','metSmiles','metHMDB',...
         'metHepatoNetID','metEHMNID'};

     names = fieldnames(model);
     %Fix incorrect cases compared to structure
     for i=1:length(met_fields)
         if isempty(find(strcmp(names,met_fields{i}),1)) %Fannst ekki
             match = find(strcmpi(names,met_fields{i}),1);
             if isempty(match) % Field not here
                 continue;
             else % Field named incorrectly regarding case.
                 met_fields{i} = names{match};
             end
         end
     end

     metabolites={};
     for i = 1:length(met_fields)
         if isfield(model,met_fields{i})
             if isempty(model.(met_fields{i})) %fields can exist but be empty.
                 metabolites(:,i) = cell(size(model.mets,1),1);
             elseif strcmp(met_fields{i}, 'metCharges') %numeric fix
                 metabolites(:,i) = num2cell(model.(met_fields{i}));
             else
                 metabolites(:,i) = model.(met_fields{i});
             end
         else %no data
             metabolites(:,i) = cell(size(model.mets,1),1);
         end
     end
     data_out = {[data; disabled], description, handles.model.genes, handles.model, metabolites};

else %input incorrect
    error('fordb is supposed to be 0 or 1.');
end

author: opencobra

repo: cobratoolbox


text:
function [x, counter] = conjugate_gradient(K, P, tol)
    % 共轭梯度法对方程进行求解

    % 默认精度
    if nargin == 2
        tol = 1e-5;
    end

    [m, n] = size(K);
    [p, q] = size(P);

    % 初始化
    x = sparse(zeros(m,1));

    % 根据P初始化 r0 and d
    r0 = P;
    d  = P;
    % 计数器
    counter = 0;

    % 如果成立，直接返回
    if sqrt(r0'*r0) < tol;
        return;
    end

    v  = K * d;
    % 计算步长
    w  = r0' * r0 / (d' * v);
    % 更新 x
    x  = x + w * d;
    % 计算剩余的
    r1 = r0 - w * v;
    counter = counter + 1;

    while sqrt(r1'*r1) > tol
        P  = r1' * r1 / (r0' * r0);
        % 更新搜索方向
        d  = r1 + P * d;
        r0 = r1;
        v  = K * d;
        % 计算步长
        w  = r0' * r0 / (d' * v);
        x  = x + w * d;
        r1 = r0 - w * v;
        counter = counter + 1;
    end

author: Meelfy

repo: FEM

save_path: github-repos/MATLAB/Meelfy-FEM


text:
function p09_default ( dim_num )

%*****************************************************************************80
%
%% P09_DEFAULT sets default values for problem 09.
%
%  Discussion:
%
%    If a problem uses vector parameters, then the spatial dimension
%    DIM_NUM is needed as input, so that the vector parameters can be
%    properly dimensioned.
%
%    Every problem keeps a count of the number of function calls.  Calling
%    this routine causes that count to be reset to 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the problem.
%
  i4 = 0;
  p09_i4 ( 'D', '*', i4 );

  r8vec = [];
  p09_r8vec ( 'D', '*', dim_num, r8vec );

  return
end


text:
function s = y2s(y)
% S = y2s(Y)
%
% Admittance to Scattering transformation
% for square matrices at multiple frequencies
%
% s = inv(I+y) * (I-y)
% ver 0.0 original	31.03.1998
% ver 0.1 +freq		27.09.2002
% ver 0.2 +octave supp.	01.06.2005

  if size(size(y),2) > 2
      nF = size(y,3);
    else nF = 1;
    end;

I = diag(ones(1, size(y,2)));

  for i=1:nF
    s(:,:,i) = inv(I+y(:,:,i)) * (I-y(:,:,i));
  end;

author: Sable

repo: mcbench-benchmarks


text:
function lse = mylogsumexp(b)
% does logsumexp across columns
B = max(b,[],2);
lse = log(sum(exp(b-repmat(B,[1 size(b,2)])),2))+B;

% Old version that used repmatC
%lse = log(sum(exp(b-repmatC(B,[1 size(b,2)])),2))+B;
end

author: zouchuhang

repo: LayoutNet


text:
function ind = subsind(pf,varargin)
% subindexing of PoleFigure data
%

while iscell(varargin{1}), varargin = [varargin{:}];end

% pf('111') or pf({'100','110'})
if iscellstr(varargin)

  h = cellfun(@(x) string2Miller(x,pf.CS),varargin,'uniformOutput',false);

elseif isa(varargin{1},'Miller') % pf(Miller(1,0,0,cs))

  h = cellfun(@(x) pf.CS.ensureCS(x),varargin);

else % pf({1,2,3})

  ind = [varargin{:}];
  return

end

ind = [];

%
for i = 1:length(h)
  for j = 1:length(pf.allH)
    if any(isnull(angle_outer(h{i},pf.allH{j})))
      ind = [ind,j]; %#ok<AGROW>
    end
  end
end

author: mtex-toolbox

repo: mtex

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


text:
% pop_writeeeg - write EEGLAB dataset to disk in EDF/GDF or BDF format
%
% pop_writeeeg( EEG ) % pops up a window
% pop_writeeeg( EEG, filename, 'key', 'val' )
%
% Inputs:
%  EEG        - EEGLAB dataset
%  filename   - [string] filename
%
% Optional keys (same as writeeeg):
%  'TYPE'         - ['GDF'|'EDF'|'BDF'|'CFWB'|'CNT'] file format for writing
%                   default is 'EDF'.
%  See writeeeg for more information
%
% Author: Arnaud Delorme, SCCN, UCSD/CERCO, 2009
%         Based on BIOSIG, sopen and swrite

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

function [command] = pop_writeeeg(EEG, filename, varargin);
command = '';

if nargin < 2
    if EEG.trials > 1
        res = questdlg2( [ 'This dataset contains data epochs.' 10 'Do you want to export the concatenated' 10 'data epochs?' ], '', 'No', 'Yes', 'Yes');
        if strcmpi(res, 'No')
            return;
        end;
    end;

	% ask user
	[filename, filepath] = uiputfile('*.*', 'Enter a file name -- pop_writeeeg()');
	if filename == 0 return; end;
	filename = fullfile(filepath,filename);

    % file format
    % -----------
    fileformats = { 'EDF' 'GDF' 'BDF' };
    uilist = { { 'style' 'text' 'String' 'File format' } ...
               { 'style' 'listbox' 'string' strvcat(fileformats) 'value' 1 } };
    geom = [1 1];
    result = inputgui( 'geometry', geom, 'uilist', uilist, 'helpcom', 'pophelp(''pop_writeeeg'')', ...
                     'title', 'Write data using BIOSIG -- pop_writeeeg()', 'geomvert', [1 2.5]);
    if length(result) == 0 return; end;

    options = { 'TYPE' fileformats{result{1}} };
else
    options = varargin;
end;

warning('off', 'MATLAB:intConvertNonIntVal');
if ~isempty(EEG.chanlocs)
    tmpchanlocs = EEG.chanlocs;
    writeeeg(filename, EEG.data(:,:), EEG.srate, 'label', { tmpchanlocs.labels }, 'EVENT', EEG.event, options{:});
else
    writeeeg(filename, EEG.data(:,:), EEG.srate, 'EVENT', EEG.event, options{:});
end;
warning('on', 'MATLAB:intConvertNonIntVal');

command = sprintf('pop_writeeeg(EEG, ''%s'', %s);', filename, vararg2str(options));

return;

repo: REST


text:
function [SwarmRBF,DeltaRBF] = RBFOperator(net,Demons,SwarmRBF,DeltaRBF,Gbest,Problem)
% OperatorPSO - The operator of particle swarm optimization.

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Parameter setting
    factor = 0;
    [N,D]  = size(SwarmRBF);
    D      = D - 1;
    BU     = Problem.upper;
    BD     = Problem.lower;

    %% SL-PSO
    for i = 1 : N
        % Choose an individual to learn one dimension to follow
        BetterInd = find(Demons(:,D+1)<SwarmRBF(i,D+1));
        ChoseId   = zeros(1,D);
        if length(BetterInd) == 1
            ChoseId = Demons(BetterInd(1),1:D);
        elseif length(BetterInd) == 2
            Mix = rand(1,D) > 0.5;
            ChoseId(Mix) = Demons(BetterInd(1),Mix);
            ChoseId(~Mix) = Demons(BetterInd(2),~Mix);
        elseif length(BetterInd) > 2
            for t = 1 : D
                ChoseId(t) = Demons(BetterInd(randperm(length(BetterInd),1)),t);
            end
        end
        if ~isempty(BetterInd)
            DeltaRBF(i,:) = rand(1,D).*DeltaRBF(i,:) + rand(1,D).*(ChoseId-SwarmRBF(i,1:D));
            % Update
            DeltaRBF(i,1:D) = max(DeltaRBF(i,1:D),BD);
            DeltaRBF(i,1:D) = min(DeltaRBF(i,1:D),BU);
            SwarmRBF(i,1:D) = DeltaRBF(i,1:D);
        end
    end
    %% Repair
    SwarmRBF(:,1:D) = max(SwarmRBF(:,1:D),repmat(Problem.lower,N,1));
    SwarmRBF(:,1:D) = min(SwarmRBF(:,1:D),repmat(Problem.upper,N,1));
    % Approximate the fitness of each particle
    SwarmRBF(:,1+D) = sim(net,SwarmRBF(:,1:D)')';
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function new_padding_value = MimRemovePaddingValues(image_wrapper, metadata, reporting)

    new_padding_value = [];
    padding_indices = [];

    % Replace padding value with zero
    if (isfield(metadata, 'PixelPaddingValue'))
        padding_value = metadata.PixelPaddingValue;

        padding_indices = find(image_wrapper.RawImage == padding_value);
    end

    % Check for unspecified padding value in GE images
    if isfield(metadata, 'Manufacturer')
        if strcmp(metadata.Manufacturer, 'GE MEDICAL SYSTEMS')
            padding_indices = find(image_wrapper.RawImage == -2000);
        end
    end

    if ~isempty(padding_indices)
        image_wrapper.RawImage(padding_indices) = max(image_wrapper.RawImage(:));
        new_padding_value = min(image_wrapper.RawImage(:)) - 1;
        image_wrapper.RawImage(padding_indices) = new_padding_value;
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
% --------------------------------------------------------------------------
function per_frame_annotation(root, folder_index)
% Read per-frame XML annotations and write bbox and track info on txt files
% Argument is folder_index [0-4]
% e.g. per_frame_annotations('/path/to/ILSVRC2015/', 0)
% --------------------------------------------------------------------------

    addpath(genpath('../..'));
    i = folder_index + 3;

    CLASS_IDS = {'n02691156','n02419796','n02131653','n02834778','n01503061','n02924116','n02958343','n02402425','n02084071','n02121808', ...
                     'n02503517','n02118333','n02510455','n02342885','n02374451','n02129165','n01674464','n02484322','n03790512','n02324045', ...
                     'n02509815','n02411705','n01726692','n02355227','n02129604','n04468005','n01662784','n04530566','n02062744','n02391049'};

%     CLASS_NAMES = {'airplane','antelope','bear','bicycle','bird','bus','car','cattle','dog','domestic_cat', ...
%                    'elephant','fox','giant_panda','hamster','horse','lion', 'lizard','monkey','motorcycle','rabbit', ...
%                    'red_panda','sheep','snake','squirrel','tiger','train','turtle','watercraft','whale','zebra'};

    OLD_FOLDERS = {'ILSVRC2015_VID_train_0000','ILSVRC2015_VID_train_0001','ILSVRC2015_VID_train_0002','ILSVRC2015_VID_train_0003','val'};
    NEW_FOLDERS = {'a','b','c','d','e'};

    root_anno = [root 'Annotations/VID/train'];
    root_imgs = [root 'Data/VID/train'];
    slash = '/';

    dir_root_anno = dir(root_anno);
    dir_root_anno = {dir_root_anno.name};
    % we always start from 3 because at 1 we have '.' and at 2 '..'
    level1 = [root_anno slash dir_root_anno{i}];
    dir_level1 = dir(level1);
    dir_level1 = {dir_level1.name};
    % j iterates across videos
    for j=3:numel(dir_level1)
        level2 = [root_anno slash dir_root_anno{i} slash dir_level1{j}];
        dir_level2 = dir(level2);
        dir_level2 = {dir_level2.name};
        xml_path = [root_anno slash dir_root_anno{i} slash dir_level1{j} slash dir_level2{3}];
        struct_from_xml = xml2struct_custom(xml_path);
        % get frame size from first frame of the video
        frame_sz = [str2double(struct_from_xml.annotation.size.height.Text), str2double(struct_from_xml.annotation.size.width.Text)];

        % k iterates across frames
        for k=3:numel(dir_level2)
            fprintf('Processing i: %d j: %d k: %d\n', i-3, j-3, k-3);
            xml_path = [root_anno slash dir_root_anno{i} slash dir_level1{j} slash dir_level2{k}];
            struct_from_xml = xml2struct_custom(xml_path);
            dir_level2_img = strrep(dir_level2{k}, '.xml', '.JPEG');
            im_path = [root_imgs slash dir_root_anno{i} slash dir_level1{j} slash dir_level2_img];
            im_anno_path = strrep(im_path, '.JPEG', '.txt');
            fid = fopen(im_anno_path,'w');
            % check if the frame contains objects at all
            if isfield(struct_from_xml.annotation, 'object')
                n_objects = numel(struct_from_xml.annotation.object);
                % k iterates across objects
                for o=1:n_objects
                    if n_objects==1
                        o_class = struct_from_xml.annotation.object(o).name.Text;
                        o_xmax = str2double(struct_from_xml.annotation.object(o).bndbox.xmax.Text);
                        o_xmin = str2double(struct_from_xml.annotation.object(o).bndbox.xmin.Text);
                        o_ymax = str2double(struct_from_xml.annotation.object(o).bndbox.ymax.Text);
                        o_ymin = str2double(struct_from_xml.annotation.object(o).bndbox.ymin.Text);
                        trackid = str2double(struct_from_xml.annotation.object(o).trackid.Text);
                    else
                        o_class = struct_from_xml.annotation.object{o}.name.Text;
                        o_xmax = str2double(struct_from_xml.annotation.object{o}.bndbox.xmax.Text);
                        o_xmin = str2double(struct_from_xml.annotation.object{o}.bndbox.xmin.Text);
                        o_ymax = str2double(struct_from_xml.annotation.object{o}.bndbox.ymax.Text);
                        o_ymin = str2double(struct_from_xml.annotation.object{o}.bndbox.ymin.Text);
                        trackid = str2double(struct_from_xml.annotation.object{o}.trackid.Text);
                    end
                    % convert ImageNet WNID of the class to integer [1-30]
                    [~, o_class] = ismember(o_class, CLASS_IDS);
                    o_sz = [o_ymax - o_ymin + 1 , o_xmax - o_xmin + 1];

                    fprintf(fid, '%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%s\n', ...
                        i, j, k, o, trackid, o_class, ...
                        frame_sz(2), frame_sz(1),  ...
                        o_xmin, o_ymin, o_sz(2), o_sz(1), im_path);
                end
            else
                    fprintf(fid, '%d,%d,%d,%d,%d,%s\n', ...
                        i, j, k, ...
                        frame_sz(2), frame_sz(1), im_path);
            end
            fclose(fid);
        end
    end
end

author: bertinetto

repo: siamese-fc