text:
function results = test_tc( dataname )

load( ['..\..\data\', dataname] );
addpath( '..\utils' );

% data.W = tenzeros( size(data.X) );
% data.W(data.Omega) = data.b;
data.W = tenzeros( size(data.X) );
N = length( size(data.W) );
data.Lamb = cell(1,N);
for i = 1:N
    data.Lamb{i} = tenzeros( size(data.X) );
end
data.Lambda = tenzeros( size(data.X) );

params.mu0 = 1e1;
params.sigma = 1e-3;
params.lambda = 1;
params.max_iter = 100;
params.opt_tol = 1e-3;
params.R = 50;

while params.lambda < 100
    results = TC_ADAL_Nuc( data, params );
%     results = TC_MSA( data, params );
    data.W = results.X;
    params.lambda = params.lambda * 5;
end

end

author: andrewssobral

repo: lrslibrary


text:
function randomfacegenerator;
% randomfacegenerator;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Written by Joseph Hollmann, February 2013
% This file may be copied, used, or modified for educational and research purposes provided
% that this header information is not removed or altered. Other distribution is prohibited
% without permission. We assume no liability for use of the code and no obligation to provide
% support, Nevertheless we would like to improve this package, and would be happy to receive comments
% Email - hollmann.j@husky.neu.edu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% As written, the function does not require any inputs
% generates a randomly created face and plots it. Nothing else

% face
dx = 0.01;
rx = .75+.1*(0.5-rand(1));
ry = 1+.1*(0.5-rand(1));
x = [-1:dx:1];
y1 = real(ry*sqrt(1-(x/rx).^2));
y2 = real(-ry*sqrt(1-(x/rx).^2));
z = find(y1 == 0);
[~,ctfnd] = max(diff(z));
ct = [1:z(ctfnd-1),z(ctfnd+1)+1:length(y1)];
y1(ct) = []; y2(ct) = []; x(ct) = [];

% smile
len = 0.5+.1*rand(1);
posy = -0.5;
amp = 0.1*(0.5-rand(1));
phs = round(rand(1)*4);
xsm = -len/2:dx:len/2;
ysm = posy+amp*cos(pi*xsm/xsm(end)+phs);

% nose
len = 0.3+0.1*rand(1);
wid = 0.3+0.1*rand(1);

xns = [-wid/3 -wid/2 0 wid/2 wid/3];
yns = [len/2 -len/4 -len/2 -len/4 len/2];

% eyes -- draw ovals
rxe = 0.1+0.1*(0.5-rand(1));
rye = 0.1+0.1*(0.5-rand(1));
xe = [-0.1:dx:0.1];
y1e = 0.3+real(rye*sqrt(1-(xe/rxe).^2));
y2e = 0.3+real(-rye*sqrt(1-(xe/rxe).^2));
xelft = -0.3+xe;
xert = 0.3+xe;
yce = 0.3;
xce = 0.3;

% hair
hrstrt = 0.25+0.75*rand(1);
hairind = find(y1 > hrstrt);

yhair = [y1(hairind);y1(hairind)+2*rand(1)];
xhair = [x(hairind);x(hairind)];

% % % %% ears
% % % rx = 0.1*rand(1);
% % % midx1 = x(1)-.1+.01*rand(1);
% % % midx2 = x(end)+.1+.01*rand(1);
% % % xer1 = midx1-rx:dx:midx1+rx;
% % % xer2 = midx2-rx:dx:midx2+rx;
% % % midy = 0.2*(0.5-rand(1));
% % % ry = 0.2*rand(1);
% % %
% % % yer11 = real(ry*sqrt(1-(xer1/rx).^2));
% % % yer12 = real(ry*sqrt(1-(xer2/rx).^2));
% % %
% % % yer21 = real(ry*sqrt(1-(xer1/rx).^2));
% % % yer22 = real(ry*sqrt(1-(xer2/rx).^2));

% define colors
p = [1 0.78 0.80];

hc = {'y','g','r','k'};
sk = {'k','y'};
e = {'b','g','k'};
et = {'*','o','.','x','+'};
lpp = 'r';

hcp = hc{randi(4,1)};
skp = sk{randi(2,1)};
ep = [e{randi(3,1)},et{randi(5,1)}];
% lpp = lp{randi(2,1)};

if strcmpi(get(gcf,'Name'),'Your New Face'),
    figure(gcf);
else,
%     figure;
    set(gcf,'name','Your New Face');
end;
plot(x,y1,skp,x,y2,skp,xsm,ysm,lpp,xns,yns,skp,xhair,yhair,hcp...
    ,xelft,y1e,skp,xelft,y2e,skp,xert,y1e,skp,xert,y2e,skp,-xce,yce,ep,xce,yce,ep,'linewidth',4)
axis tight
set(gca,'XTickLabel',[],'yticklabel',[],'xtick',0,'ytick',0)

author: Sable

repo: mcbench-benchmarks


text:
function test_ft_appendfreq

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_appendfreq

% make some dummy frequency structures
freq1.label = {'1';'2'};
freq1.freq  = 1:10;
freq1.time  = 1:5;
freq1.dimord = 'chan_freq_time';
freq1.powspctrm = randn(2,10,5);

cfg = [];
cfg.parameter = 'powspctrm';

freq2         = freq1;
cfg.appenddim = 'rpt';
freqrpt       = ft_appendfreq(cfg, freq1, freq2);
cfg.appenddim = 'auto';
freqrptauto   = ft_appendfreq(cfg, freq1, freq2);

freq2         = freq1;
freq2.label   = {'3';'4'};
cfg.appenddim = 'chan';
freqchan      = ft_appendfreq(cfg, freq1, freq2);
cfg.appenddim = 'auto';
freqchanauto  = ft_appendfreq(cfg, freq1, freq2);

freq2         = freq1;
freq2.freq    = 11:20;
cfg.appenddim = 'freq';
freqfreq      = ft_appendfreq(cfg, freq1, freq2);
if ~isfield(freqfreq,'freq')
  error('freq field is not appeneded: see bugs 1984 and 2187');
end

cfg.appenddim = 'auto';
freqfreqauto  = ft_appendfreq(cfg, freq1, freq2);
if ~isfield(freqfreqauto,'freq')
  error('freq field is not appeneded: see bugs 1984 and 2187');
end

freq2         = freq1;
freq2.time    = 6:10;
cfg.appenddim = 'time';
freqtime      = ft_appendfreq(cfg, freq1, freq2);
if ~isfield(freqfreqauto,'freq')
  error('freq field is not appeneded: see bugs 1984 and 2187');
end

cfg.appenddim = 'auto';
freqtimeauto  = ft_appendfreq(cfg, freq1, freq2);
if ~isfield(freqtimeauto,'freq')
  error('freq field is not appeneded: see bugs 1984 and 2187');
end

% now test for numerical inaccurracies, should concatenate across 'rpt'
freq2          = freq1;
freq2.time     = freq1.time+0.0000001;
cfg.appenddim  = 'auto';
freqrpt2       = ft_appendfreq(cfg, freq1, freq2);
if ~isfield(freqrpt2,'freq')
  error('freq field is not appeneded: see bugs 1984 and 2187');
end

%% test for data with labels shuffled around
freq1 = [];
freq1.label = {'1';'2'};
freq1.freq  = 1:10;
freq1.time  = 1:5;
freq1.dimord = 'chan_freq_time';
freq1.powspctrm = randn(2,10,5);

freq2 = freq1;
freq2.label = {'2';'1'};
freq2.powspctrm = randn(2,10,5);

freqshuffled = ft_appendfreq(cfg, freq1, freq2);

if ~strcmp(freqshuffled.dimord, 'rpt_chan_freq_time')
  error('unexpected dimord when appending freqs with differently permuted labels');
end
if ~isfield(freqshuffled,'freq')
  error('freq field is not appeneded: see bugs 1984 and 2187');
end

% now check whether channels were correctly appended
[a,b] = match_str(freq1.label, freqshuffled.label);
x1 = freq1.powspctrm(a(1),:,:);
y1 = freqshuffled.powspctrm(1,b(1),:,:);

[a,b] = match_str(freq2.label, freqshuffled.label);
x2 = freq2.powspctrm(a(1),:,:);
y2 = freqshuffled.powspctrm(2,b(1),:,:);

if ~all(x1(:) == y1(:)) || ~all(x2(:) == y2(:))
  error('data was wrongly appended when channel labels are differently ordered in input arguments');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


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

l = isnan(tmap);
tmap(l) = 1;

%l = tmap< 0.1;
%tmap(l) = nan;

[lat,lon] = meshgrat(tmap,tmapleg);
%[smap,smapleg] = country2mtx('switzerland',100);
%[lat0, lon0] = meshgrat(smap,smapleg);

% tmap = km2deg(tmap/1);
[X , Y]  = meshgrid(gx,gy);

%sw = interp2(lon0,lat0,smap,lon,lat);

ren = interp2(X,Y,re4,lon,lat);

mi = min(min(ren));
l =  isnan(ren);
ren(l) = mi-100;

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

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

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

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

load worldlo
%h = displaym(POline); set(h(1),'color',[0.9 0.9 0.9],'Linewidth',1.7)
h2 = displaym(PPpoint);
%h = displaym(PPtext); trimcart(h);
plotm(mainfault(:,2), mainfault(:,1),'m','Linewidth',4);

pl = plotm(ma(:,2),ma(:,1),'hw');
set(pl,'LineWidth',1.5,'MarkerSize',12,...
    'MarkerFaceColor','y','MarkerEdgeColor','k')

zdatam(handlem('allline'),10000) % keep line on surface
%zdatam(handlem('alltext'),10000) % keep line on surface

j = jet;
j = j(64:-1:1,:);
j = [ [ 0.85 0.9 0.9] ; j];
caxis([ 95 2000]);

colormap(j); brighten(0.1);

axis off; set(gcf,'color','k')

setm(gca,'ffacecolor','k')
setm(gca,'fedgecolor','w','flinewidth',3);

setm(gca,'mlabellocation',0.5)
setm(gca,'meridianlabel','on')
setm(gca,'plabellocation',0.5)
setm(gca,'parallellabel','on')
setm(gca,'Fontcolor','w','Fontweight','bold','FontSize',15,'Labelunits','dm')

h5 = colorbar;
set(h5,'position',[0.8 0.35 0.01 0.3],'TickDir','out','Ycolor','w','Xcolor','w',...
    'Fontweight','bold','FontSize',15);
set(gcf,'Inverthardcopy','off');

author: CelsoReyes


text:
%% FILESYNC Synchronize data using a USB key.
%
%   PLEASE MAKE SURE YOU UNDERSTAND THE CODE COMPLETELY BECAUSE IT DOES OVERWRITE DATA!
%
%   This file is designed to synchronize data between two computers using a
%   USB key. It only updates files that have been modified since the last
%   update. Files with the more recent modified date overwrite existing
%   files.
%
%   The first time this function is executed, it initializes the USB key.
%   However, since it doesn't know when the last backup was performed, no
%   data is copied to the USB key. Each time after the initialization,
%   this function will update the USB key to keep data between the two
%   computer synchronized.
%
%   Variables that need to be set are:
%      <basedir>    - (1 x 2) CELLL ARRAY of STRINGs of full paths to the base directory
%                     on each computer to be synchronized.
%      <staticname> - (1 x 2) CELL ARRAY of STRINGs that include the computer name
%                     of each computer retruned by "getenv('computername')"
%      <USBdir>     - STRING of directory path on the USB key used for synchronization
%      <USBtmp>     - STRING of directory path on the USB key used for temporary
%                     purposes during synchronization
%
% This function was tested on Windows XP, SP2, running Matlab 2006b (ver 7.3).
% There is a GETENV function call that expects 'computername' to exist as an
% environmet variable. This will need to be modified for other operating
% systems.
%
% Please send comments/suggestions to:
%
% Jesse Norris
% jnorris@wfubmc.edu
% http://jesse.musculoskeletaldynamics.com
%
% Last modified: 15 April 2007
%

%% filesync
function filesync

   % Start stopwatch, clear command window and store current directory
   %
   tic; clc; currdir = cd;

%% ~~~~~ USER SPECIFIC VARIABLES ~~~~~

   % Set the static computer names and base directories on those computers
   % that will be synchronized. These static computer names should match those
   % returned by the GETENV function with 'computername' argument.
   %
   staticname{1} = '3390-99RGXAZ';  % Computer 1 <computername>
   basedir{1}    = 'C:\Documents and Settings\jnorris\My Documents'; % Base directory on Computer 1
   staticname{2} = '1612-8QDLG11';   % Computer 2 <computername>
   basedir{2}    = 'C:\Documents and Settings\jnorris\My Documents'; % Base directory on Computer 2

   % Set the directories on the USB key used for synchronization.
   %
   USBdir   = 'E:\File Transfer';                  % Directories for storing updated files
   USBtmp   = 'E:\Temp';                           % Temporary directory

%% ~~~~~ END USER SPECIFIC VARIABLES ~~~~~

   % Retrieve current computer name. NOTE: The <getenv> function assumes
   % 'computername' is an enviroment variable.
   %
   computername  = getenv('computername');
   if isempty(computername); error('''computername'' does not exist in this computer''s environment variables.'); end;

   % Set files for keeping track of synchronization times.
   %
   USBfile  = ['\' '0 Files updated.txt'];         % File where updated files are stored.
   USBts{1} = ['\1 ' staticname{1} ' to USB.txt']; % Files on the USB key that are used solely for time stamps.
   USBts{2} = ['\1 ' staticname{2} ' to USB.txt'];

   % Check if directories & file exist. If not, attempt to create. If create fails, return error.
   %
   chkmkdir(USBdir);
   chkmkfile([USBdir USBfile]);
   chkmkfile([USBdir USBts{1}]);
   chkmkfile([USBdir USBts{2}]);

   % Ensure that user is ready to begin synchronization.
   %
   button = questdlg('Begin synchronization?', 'USB Synchronization', 'Yes', 'No', 'No');

   if strcmp(button,'Yes');

       % Determine which computer was most recently synchronized to the USB key
       %
       comp1date = filemodified([USBdir USBts{1}]);
       comp2date = filemodified([USBdir USBts{2}]);
       if (comp1date > comp2date); USBupdatedfrom = staticname{1}; else USBupdatedfrom = staticname{2}; end;

       % Collapse <basedir> and <USBts> to just STRING variables based on the current computer.
       %
       if strcmp(computername,staticname{1});
           basedir = basedir{1};
           USBts   = USBts{1};
       elseif strcmp(computername,staticname{2});
           basedir = basedir{2};
           USBts   = USBts{2};
       else
           error('filesync:unknowncomp','<computername> ''%s'' does not match either name in the variable <staticname>',computername);
       end;

       % If the current computer was the most recent to be synchronized with the USB key,
       % it is only necessary to update the USB key with any new data.
       %
       if strcmp(computername,USBupdatedfrom);
           updateUSB(basedir,USBdir,USBfile,USBts); % Update the USB key

       % Otherwise, it is necessary to
       % (1) Update any changes that have occured on the current computer to <USBtmp> on the USB key.
       % (2) Then update the current computer with USB data.
       % (3) Lastly, copy from <USBtmp> to <USBdir> for updates back to the other computer.
       %
       else
           chkmkdir(USBtmp);                          % Check if <USBtmp> exists. If not, attempt to create. If create fails, return error.
           copyfile([USBdir USBts],[USBtmp USBfile]); % Copy the time stamp file to the update file in the temporary directory
           updateUSB(basedir,USBtmp,USBfile,USBts);   % Update <USBtmp>  with data on computer
           updatecomputer(basedir,USBdir,USBfile);    % Update computer with data in <USBdir> on USB key,
           movefile([USBtmp '\*'],USBdir);            % Move all files from the <USBtmp> to <USBdir>, and
           rmdir(USBtmp);                             % Remove <USBtmp> from the USB key
       end;

       % Return to original directory
       %
       try
           cd(currdir);
       catch
           cd(basedir);
           warning('filesync:cderror','Unable to change directory to: %s',currdir);
       end

       toc; % Elapsed time

       fprintf(' ~~~ \n\n');
   else
       fprintf('Synchronization Cancelled.\n\n');
   end;
end

%% chkmkdir
function message = chkmkdir(dirname)
   % Check if <dirname> is a directory. If not, attempt to create.
   % If create fails return error.
   %
   if ~exist(dirname,'dir')
       try
           mkdir(dirname); % Create directory
           message = 'Directory created.';
       catch
           error('Could not create path: ''%s''.\n~Check initialization value of variable: <%s>.',dirname,inputname(1));
       end;
   else
       message = 'Directory exists.';
   end;
end

%% chkmkfile
function message = chkmkfile(fname)
   % Check if <fname> is a file. If not, attempt to create.
   % If create fails return error.
   %
   if ~exist(fname,'file')
       try
           fid = fopen(fname,'w'); fclose(fid); % Create file
           message = 'File created.';
       catch
           error('Could not create file: ''%s''.\n~Check initialization value of variable: <%s>.',fname,inputname(1));
       end;
   else
       message = 'File exists.';
   end;
end

%% removeemptyfolders
function removeemptyfolders(rmdirname)
   if isdir(rmdirname);
       cd(rmdirname);          % CD to directory
       while length(dir) == 2; % Recurrisvely remove until an unempty folder is encountered
           rmdirname = cd;     % Current directory
           cd ..;              % CD to parent directory
           rmdir(rmdirname);   % Remove the subdirectory
       end;
   end;
end

%% filemodified
function modified = filemodified(fullpath)
   if exist(fullpath,'file')
       fileinfo = dir(fullpath);   % Determine last modification time
       modified = datenum(fileinfo.date,'dd-mmm-yyyy HH:MM:SS');
   else
       modified = [];
   end;
end

%% updatecomputer
function updatecomputer(basedir,USBdir,USBfile)
   fprintf('\n... Updating the computer with files on USB key ...\n\n');

   USBfile = [USBdir USBfile];                 % Update to full path
   fid     = fopen(USBfile,'r');                 % Open file that includes full directory paths and file names
   files   = textscan(fid,'%s','delimiter',';'); % Load all the files to update
   files   = files{1};
   fclose(fid);                                % Close the file

   % Step through each file
   % -> check that file exists on USB key
   % -> check that the date on the USB key is more recent than that on the computer
   % -> copy file from USB key to computer
   for fileind = 1:length(files)
       [subdir name ext] = fileparts(files{fileind});
       fprintf('%s...',[subdir '\' name ext]); % Update workspace

       if exist([USBdir subdir '\' name ext],'file')
           USBdate = filemodified([USBdir  subdir '\' name ext]);
           Cdate   = filemodified([basedir subdir '\' name ext]);
           if ~isempty(Cdate)
               if  USBdate > Cdate
                   movefile([USBdir subdir '\' name ext],[basedir subdir '\' name ext],'f');
                   fprintf('%s\n','Updated.');
               else
                   switch questdlg(['Confirm Overwrite of ' name ext],'Overwrite','Yes','No - Delete from USB','No - Leave on USB','No - Leave on USB')
                       case 'Yes'
                           movefile([USBdir subdir '\' name ext],[basedir subdir '\' name ext],'f');
                           fprintf('%s\n','Updated.');
                       case 'No - Delete from USB'
                           delete([USBdir subdir '\' name ext]);
                           fprintf('%s\n','Not updated - deleted from USB.');
                       otherwise
                           fprintf('%s\n','Not updated - left on USB.');
                   end;
               end;
           else
               if ~isdir([basedir subdir]); mkdir([basedir subdir]); end; % Create new directory
               movefile([USBdir subdir '\' name ext],[basedir subdir '\' name ext],'f');
               fprintf('%s\n','Updated.');
           end;
       else
           fprintf('%s\n','File not found on USB key.');
       end;
       removeemptyfolders([USBdir subdir]); % Recurrisvely remove until an unempty folder is encountered
   end;
end

%% updateUSB
function updateUSB(basedir,USBdir,USBfile,USBts)
   fprintf('\n... Updating the USB key with files from the computer ...\n\n');

   USBfile = [USBdir USBfile];      % Update to full path
   USBts   = [USBdir USBts];        % Update to full path
   USBdate = filemodified(USBfile); % Check when it was modified
   fid     = fopen(USBfile,'a');    % Open file to append the full path strings to

   % Load all the paths to scan. This expression is based on a line in
   % Matlab function PARSEDIRS.
   %
   paths = regexp(genpath(basedir), sprintf('[^\\s%s][^%s]*', pathsep, pathsep), 'match')';

   % Step through all paths
   % -> examine if files have been modified since <USBdate>
   % -> copy updated files to the USB key
   % -> append <fid> with directory where file was originally from
   %
   for pathind = 1:length(paths)
       cd(paths{pathind}); % Change to directory
       files = dir;        % Load file information

       % Step through all files
       %
       for fileind = 3:length(files)
           % NOTE: THIS EXPRESSIONS STOPS THIS FUNCTION FROM COPYING OUTLOOK MAILBOX FILES (*.pst)
           % The one on my computer is too large to be copied with just a USB key.
           %
           if ~files(fileind).isdir && (datenum(files(fileind).date,'dd-mmm-yyyy HH:MM:SS') > USBdate) && ~strcmpi(files(fileind).name(end-2:end),'pst')
               if (files(fileind).bytes > 10e6)
                   okay2copy = questdlg(['File ' files(fileind).name ' is ' sprintf('%.0f',files(fileind).bytes/1e6) ' Mb. Copy file to the USB key?'],'Large file.','Yes','No','No');
                   okay2copy = strcmp(okay2copy,'Yes');
               else
                   okay2copy = true;
               end;
               if okay2copy
                   try
                       if ~isdir([USBdir paths{pathind}(47:end)]); mkdir([USBdir paths{pathind}(47:end)]); end; % Create directory on USB key
                       copyfile(files(fileind).name,[USBdir paths{pathind}(47:end) '\' files(fileind).name]);   % Copy file to USB key
                       fprintf(fid,'%s\r\n',[paths{pathind}(47:end) '\' files(fileind).name '; ']);             % Update the txt file
                       fprintf(      '%s\n',[paths{pathind}(47:end) '\' files(fileind).name]);                  % Update workspace
                   catch
                       warning('updateUSB:copyfile','Error caught in function <updateUSB>');
                       handleerror(lasterror);
                   end
               else
                   fidnoupdate = fopen([USBdir '\0 Files NOT updated'],'a');
                   fprintf(fidnoupdate,'%s\r\n',[paths{pathind}(47:end) '\' files(fileind).name '; ']);     % Update the txt file
                   fclose(fidnoupdate);
               end;
           end;
       end;
   end;
   fclose(fid); % Close the file

   fid = fopen(USBts,'w'); fclose(fid); % Update successful, replace appropriate time stamp file

   cd(basedir); % Return to base directory
end

%% genpath
% p = genpath(d)
% OUTPUT arguments
%      <p>: a path string that includes all subdirectories of <d>
%
% INPUT arguments
%      <d>: returns a path string starting in D, plus, recursively, all
%           the subdirectories of D, including empty subdirectories.
%
% This function is based on Matlab's genpath function. It was written because
% the genpath function skips method directories (i.e. those beginning with @),
% and 'private' directories.
%
function p = genpath(d)
  p = '';                % Initialize path string
  files = dir(d);        % Generate path based on given root directory
  if isempty(files)
      return;
  end;
  p     = [p d pathsep]; % Add d to the path even if it is empty.
  isdir = [files.isdir]; % Create logical vector for subdirectory entries in d
  dirs  = files(isdir);  % Select only directory entries from the current list

  for i = 1:length(dirs)
      dirname = dirs(i).name;
      if ~strcmp(dirname,'.') && ~strcmp(dirname,'..')
          p = [p genpath(fullfile(d,dirname))]; % recursive calling of this function.
      end;
  end;
end

%% handleerror
% handleerror(errorin)
% INPUT arguments
%      <errorin>: a structure containing at least the following two fields:
%
%        message    : the text of the error message
%        identifier : the message identifier of the error message
%
%      Additionally, <errorin> may contain field <stack>.
%
function handleerror(errorin)
  try
      fprintf('Message:    %s\n',errorin.message);
      fprintf('Identifier: %s\n',errorin.identifier);
      fprintf('Stack:\n');
      if(isfield(errorin,'stack'));
          while(length(errorin.stack) >=1)
              %fprintf('      %s\n',errorin.stack(1).file);
              fprintf('       <%s> line: %d\n',errorin.stack(1).name,...
                                               errorin.stack(1).line);
              errorin.stack = errorin.stack(2:end);
          end;
      end;
  catch
      warning('handleerror:errorcaught','Error in function <handleerror>');
  end
end


text:
function numgrad = computeNumericalGradient_SR(image, label, weight_idx)
    global config mem;
    estimatedGrad = mem.grads;
    epsilon = config.NEW_MEM(0.01);
    % Initialize numgrad with zeros
    numgrad = zeros(size(config.weights{1}{weight_idx}));

    for x = 1:size(config.weights{1}{weight_idx}, 1)
        for y = 1:size(config.weights{1}{weight_idx}, 2)
            config.weights{1}{weight_idx}(x, y) = config.weights{1}{weight_idx}(x, y) + epsilon;
            op_train_pipe(image, label);
            cost1 = config.cost;
            config.weights{1}{weight_idx}(x, y) = config.weights{1}{weight_idx}(x, y)  - (2*epsilon);
            op_train_pipe(image, label);
            cost2 = config.cost;
            config.weights{1}{weight_idx}(x, y) = config.weights{1}{weight_idx}(x, y) + epsilon;

            % compute the numerical
            temp = (cost1 - cost2) / (2 * epsilon);
            numgrad(x, y) = gather(temp);

            diff = norm(temp-estimatedGrad{1}{weight_idx}(x, y));
            threshold = 10^-10;
            %if(diff > threshold)
                fprintf('numerical: %i, ', temp);
                fprintf('estimated: %i. \n', estimatedGrad{1}{weight_idx}(x, y));
                fprintf('loop (%d, %d)\n', x, y);
                fprintf('diff too large! diff: %i.\n', diff);
            %end
        end
    end
end

author: jimmy-ren


text:
function D=imdir(dirname,filter)

% IMDIR lists image content in the directory specified by dirname.
%   Returns a structure similar to the one returned by MATLAB's dir
%   function. Additional filter on the image file name can be specified.
%
% See also DIR
% (c) Michael Rubinstein, MIT
%

if nargin<2
    filter='';
end

if ~exist(dirname,'dir')
    error([dirname ' not found.']);
end

D=dir(fullfile(dirname,sprintf('*%s*',filter)));
D=D(cellfun(@isImageFile,{D(:).name}));

function tf=isImageFile(filename)

% supported formats
formats = {'jpg','jpeg','png','bmp','pgm','gif','ppm','tif','tiff'};

[path,name,ext]=fileparts(filename);
ext = ext(2:end); % remove '.'

tf = any(strcmpi(ext,formats));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
 function ph = ir_unwrap(ph, varargin)
%function ph = ir_unwrap(ph, varargin)
%|
%| Unwrap phase along all dimensions.
%|
%| 2014-08-19, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if streq(ph, 'test'), ir_unwrap_test, return, end

for id = 1:ndims(ph)
	ph = unwrap(ph, [], id);
end

% ir_unwrap_test
function ir_unwrap_test
p0 = linspace(-pi,2*pi,98)' * linspace(1,2,99);
p0 = p0.'; % octave worked better this way!?
x = exp(1i * p0);
ph = angle(x);
pu = ir_unwrap(ph);
clim = [-1 2]*pi;
im plc 2 2
im(1, p0, 'true', clim), cbar
im(2, ph, 'angle', clim), cbar
im(3, pu, 'unwrapped', clim), cbar
im(4, angle(exp(1i*(pu -p0)))), cbar

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function ASAversionchk(ASAcontrol)
%ASAVERSIONCHK Function version check
%   ASAVERSIONCHK(ASACONTROL) checks the compatibility of two functions
%   according to the information provided in the structure variable
%   ASACONTROL, if called from within the body of a m-file. If
%   incompatibility is detected, it halts execution and prints an error
%   message.
%
%   To apply a version check procedure when for example CALLER_FUNCTION
%   calls CALLED_FUNCTION, the latter must invoke ASAVERSIONCHK.
%   CALLER_FUNCTION must supply ASACONTROL as the last input argument
%   when it calls CALLED_FUNCTION. ASACONTROL must contain a (nested)
%   field named 'req_version.CALLED_FUNCTION' containing a version
%   identifier that defines the minimum version of CALLED_FUNCTION, as
%   required by CALLER_FUNCTION to guarantee consistency. In addition,
%   CALLED_FUNCTION must add two fields to ASACONTROL, named 'is_version'
%   and 'comp_version' prior to invoking ASAVERSIONCHK. Both fields
%   containing version identifiers, are used respectively to identify
%   CALLED_FUNCTION itself, with its own actual version, and to define
%   its downward compatibility with former versions. Operable version
%   check implementations can be found in all ARMASA main functions.
%
%   See also: ASACONTROL, STRUCT

caller_stack = dbstack;
stack_length = length(caller_stack);
if stack_length == 1
   error(ASAerr(7))
else
   [caller_path chk_caller] = fileparts(caller_stack(2).name);
   if stack_length >= 3
      [caller_path funct_caller] = fileparts(caller_stack(3).name);
   end
   if ~any(strcmp(fieldnames(ASAcontrol),'is_version'))|~any(strcmp(fieldnames(ASAcontrol),'comp_version'))
      error(ASAerr(22))
   end
end

if ~any(strcmp(fieldnames(ASAcontrol),'req_version'))
   chk = 0;
elseif ~any(strcmp(fieldnames(ASAcontrol.req_version),chk_caller))
   chk = 0;
else
   chk = 1;
   v = ASAcontrol.is_version;
   formatchk(v,'ASAcontrol.is_version',chk_caller);
   is_version = datenum(v(1),v(2),v(3),v(4),v(5),v(6));
   v = ASAcontrol.comp_version;
   formatchk(v,'ASAcontrol.comp_version',chk_caller);
   comp_version = datenum(v(1),v(2),v(3),v(4),v(5),v(6));
   v = getfield(ASAcontrol.req_version,chk_caller);
   formatchk(v,['ASAcontrol.req_version' chk_caller],chk_caller);
   req_version = datenum(v(1),v(2),v(3),v(4),v(5),v(6));
end

if isequal(chk,0)
   if stack_length > 2
      show_warn = 0; % <----  User adjustable: show_warn = 1 will enable warning messages
   else
      show_warn = 0;
   end
   if show_warn
      v = ASAcontrol.is_version;
      formatchk(v,'ASAcontrol.is_version',chk_caller);
      string = sprintf('%0.0f %0.0f %0.0f %0.0f %0.0f %0.0f',v);
      warning(ASAwarn(8,...
         {funct_caller;chk_caller;funct_caller;chk_caller;...
            string;chk_caller;funct_caller},ASAcontrol))
   end
elseif req_version > is_version
   if stack_length == 2
      error(ASAerr(3,chk_caller))
   else
      error(ASAerr(4,{chk_caller;funct_caller}))
   end
elseif req_version < comp_version
   if stack_length == 2
      error(ASAerr(5,chk_caller))
   else
      error(ASAerr(6,{chk_caller;funct_caller}))
   end
end

function formatchk(version,name_version,chk_caller)

if ~isnumeric(version) | ~isequal(length(version),6)
   error(ASAerr(10,{name_version;chk_caller}))
end

%Program history
%======================================================================
%
% Version                Programmer(s)          E-mail address
% -------                -------------          --------------
% [2000 12 30 20 0 0]    W. Wunderink           wwunderink01@freeler.nl


text:
function zi = pwl_interp_2d ( nxd, nyd, xd, yd, zd, ni, xi, yi )

%*****************************************************************************80
%
%% PWL_INTERP_2D: piecewise linear interpolant to data defined on a 2D grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NXD, NYD, the number of X and Y data values.
%
%    Input, real XD(NXD), YD(NYD), the sorted X and Y data.
%
%    Input, real ZD(NXD,NYD), the Z data.
%
%    Input, integer NI, the number of interpolation points.
%
%    Input, real XI(NI), YI(NI), the coordinates of the
%    interpolation points.
%
%    Output, real ZI(NI), the value of the interpolant.
%
  zi = zeros ( ni, 1 );

  for k = 1 : ni

    i = r8vec_bracket5 ( nxd, xd, xi(k) );
    if ( i == -1 )
      zi(k) = r8_huge ( );
      continue
    end

    j = r8vec_bracket5 ( nyd, yd, yi(k) );
    if ( j == -1 )
      zi(k) = r8_huge ( );
      continue
    end

    if ( yi(k) < yd(j+1) ...
      + ( yd(j) - yd(j+1) ) * ( xi(i) - xd(i) ) / ( xd(i+1) - xd(i) ) )

      dxa = xd(i+1) - xd(i);
      dya = yd(j)   - yd(j);

      dxb = xd(i)   - xd(i);
      dyb = yd(j+1) - yd(j);

      dxi = xi(k)   - xd(i);
      dyi = yi(k)   - yd(j);

      det = dxa * dyb - dya * dxb;

      alpha = ( dxi * dyb - dyi * dxb ) / det;
      beta =  ( dxa * dyi - dya * dxi ) / det;
      gamma = 1.0 - alpha - beta;

      zi(k) = alpha * zd(i+1,j) + beta * zd(i,j+1) + gamma * zd(i,j);

    else

      dxa = xd(i)   - xd(i+1);
      dya = yd(j+1) - yd(j+1);

      dxb = xd(i+1) - xd(i+1);
      dyb = yd(j)   - yd(j+1);

      dxi = xi(k)   - xd(i+1);
      dyi = yi(k)   - yd(j+1);

      det = dxa * dyb - dya * dxb;

      alpha = ( dxi * dyb - dyi * dxb ) / det;
      beta =  ( dxa * dyi - dya * dxi ) / det;
      gamma = 1.0 - alpha - beta;

      zi(k) = alpha * zd(i,j+1) + beta * zd(i+1,j) + gamma * zd(i+1,j+1);

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function p = hen_polynomial_value ( m, n, x )

%*****************************************************************************80
%
%% HEN_POLYNOMIAL_VALUE evaluates Hen(i,x).
%
%  Discussion:
%
%    Hen(i,x) is the normalized probabilist's Hermite polynomial of degree I.
%
%    These polynomials satisfy the orthonormality condition:
%
%      Integral ( -oo < X < +oo ) exp ( - 0.5 * X^2 ) * Hen(M,X) Hen(N,X) dX
%      = delta ( N, M )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Frank Olver, Daniel Lozier, Ronald Boisvert, Charles Clark,
%    NIST Handbook of Mathematical Functions,
%    Cambridge University Press, 2010,
%    ISBN: 978-0521192255,
%    LC: QA331.N57.
%
%  Parameters:
%
%    Input, integer M, the number of evaluation points.
%
%    Input, integer N, the highest order polynomial to compute.
%    Note that polynomials 0 through N will be computed.
%
%    Input, real X(M,1), the evaluation points.
%
%    Output, real P(M,N+1), the values of the polynomials of index 0 through N.
%
  p = zeros ( m, n + 1 );

  p(1:m,1) = 1.0;

  if ( n == 0 )
    return
  end

  p(1:m,2) = x(1:m,1);

  for j = 2 : n
    p(1:m,j+1) = x(1:m,1) .* p(1:m,j) - ( j - 1 ) * p(1:m,j-1);
  end
%
%  Normalize.
%
  d = diag ( 1.0 ./ sqrt ( gamma ( 1:n+1) .* sqrt ( 2.0 * pi ) ) );

  p = p * d;

  return
end


text:
function [yarray0, dyarray0, tarray] = calctraj(via,dt,tacc)
ptlen = size(via,2);
y0 = via(1,:);
yarray0 = mstraj(via(2:end,:),ones(1,ptlen),[],y0,dt,tacc); % row vector for each xr
yarray0 = [y0;yarray0];
dyarray0 = [diff(yarray0)./dt;zeros(1,ptlen)];
tarray = 0:dt:(size(yarray0,1)-1)*dt;
end

repo: paper-simulation


text:
super = 2;
opt.filter_type = {'morlet_1d','spline_1d'};
opt.B = [8 1];
opt.Q = [super*8 1];
opt.filter_format = 'fourier_truncated';
M = 3;

Js = [9:2:15];

N = 2^17;

src = phone_src('/path/to/timit');
files = src.files;
rs = RandStream.create('mt19937ar','Seed',floor(pi*1e9));
files = files(rs.randperm(length(files)));

scat_opt.M = M;
scat_opt.oversampling = 2;

for l = 1:length(Js)
	opts{l} = opt;

	opts{l}.J = T_to_J(2^Js(l),opts{l});

	[Wop{l},filters{l}] = wavelet_factory_1d(N, opts{l}, scat_opt);
end

E = zeros(length(Js),M+1,length(files));

for k = 1:length(files)
	fprintf('%s\n',files{k});
	x = data_read(files{k});
	x = x(1:min(N,length(x)));
	x = [x; zeros(N-length(x),1)];
	x = x-mean(x);
	x = x/sqrt(sum(abs(x(:)).^2));
	for l = 1:length(Js)
		[t,meta] = format_scat(scat(x,Wop{l}));

		t = squeeze(t);
		for m = 0:max(meta.order)
			t1 = t(meta.order==m,:);
			E(l,m+1,k) = sum(abs(t1(:)).^2);
		end
	end
	mean(E(:,:,1:k),3)
end

% At the end, we should have
%	0.0000    0.9453    0.0484    0.0023
%	0.0000    0.6800    0.2900    0.0191
%	0.0000    0.3486    0.5325    0.1156
%	0.0000    0.2772    0.5607    0.2471

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function DEM_demo_ALAP
% This demonstration is essentially the same as DEM_demo_LAP - however
% here, we compare two generalised filtering schemes that are implemented
% very differently: the first integrates the generative process in
% parallel with the inversion, while the standard spm_LAP scheme inverts a
% model given pre-generated data. The advantage of generating and modelling
% data  contemporaneously is that it allows the inversion scheme to couple
% back to the generative process through action (see active inference
% schemes): spm_ALAP.
%__________________________________________________________________________
% Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: DEM_demo_ALAP.m 7679 2019-10-24 15:54:07Z spm $

% get basic convolution model
%==========================================================================
M       = spm_DEM_M('convolution model');

% gradient functions for speed (not implemented here)
%--------------------------------------------------------------------------
% M(1).fx = inline('P.f','x','v','P');
% M(1).fv = inline('P.h','x','v','P');
% M(1).gx = inline('P.g','x','v','P');
% M(1).gv = inline('sparse(4,1)','x','v','P');

M(1).E.nN = 8;                                 % number of E steps
M(1).E.nE = 8;                                 % number of E steps
M(1).E.nD = 1;                                 % number of time steps
M(1).E.s  = 1;                                 % smoothness
G(1).E.s  = 1;                                 % smoothness
M(1).E.d  = 2;                                 % order
M(1).E.n  = 6;                                 % order

% free parameters
%--------------------------------------------------------------------------
P       = M(1).pE;                             % true parameters
ip      = [1 2 5 9];                           % free parameters
pE      = spm_vec(P);
np      = length(pE);
pE(ip)  = 0;
pE      = spm_unvec(pE,P);
pC      = sparse(ip,ip,exp(4),np,np);
M(1).pE = pE;
M(1).pC = pC;

% free hyperparameters
%--------------------------------------------------------------------------
M(1).Q  = {speye(M(1).l,M(1).l)};
M(1).R  = {speye(M(1).n,M(1).n)};
M(1).hE = 8;
M(1).gE = 6;
M(1).hC = 1/4;
M(1).gC = 1/4;

% generative process
%==========================================================================
G(1).f  = M(1).f;
G(1).g  = M(1).g;
G(1).x  = M(1).x;
G(1).V  = exp(8);
G(1).W  = exp(6);
G(1).pE = P;

G(2).v  = 0;
G(2).V  = exp(16);

% hidden cause
%--------------------------------------------------------------------------
N      = 32;
U      = exp(-((1:N) - 12).^2/(2.^2));

% invert
%==========================================================================
DEM.M  = M;
DEM.G  = G;
DEM.C  = U;

% generate and filter responses
%--------------------------------------------------------------------------
LAP    = spm_ALAP(DEM);

% filter generated responses
%--------------------------------------------------------------------------
DEM.Y  = LAP.Y;
DEM.pU = LAP.pU;
DEM.pP = LAP.pP;
DEM    = spm_LAP(DEM);

% Show results for LAP (standard scheme)
%==========================================================================
spm_figure('GetWin','Figure 1: Generalised filtering - standard scheme');

% overlay true values
%--------------------------------------------------------------------------
spm_DEM_qU(DEM.qU,DEM.pU)

% parameters
%--------------------------------------------------------------------------
qP    = spm_vec(DEM.qP.P);
qP    = qP(ip);
tP    = spm_vec(DEM.pP.P);
tP    = tP(ip);

subplot(2,2,4)
bar([tP qP])
axis square
legend('true','GF - standard')
title('parameters','FontSize',16)

cq    = 1.64*sqrt(diag(DEM.qP.C(ip,ip)));
for i = 1:length(qP),hold on
    plot([i i] + 1/8,qP(i) + [-1 1]*cq(i),'LineWidth',4,'color','r')
end, hold off

% Show results for ALAP (parallel scheme)
%==========================================================================
spm_figure('GetWin','Figure 2: Generalised filtering - parallel scheme');

% overlay true values
%--------------------------------------------------------------------------
spm_DEM_qU(LAP.qU,LAP.pU)

% parameters
%--------------------------------------------------------------------------
qP    = spm_vec(LAP.qP.P);
qP    = qP(ip);
tP    = spm_vec(LAP.pP.P);
tP    = tP(ip);

subplot(2,2,4)
bar([tP qP])
axis square
legend('true','GF - parallel')
title('parameters','FontSize',16)

cq    = 1.64*sqrt(diag(LAP.qP.C(ip,ip)));
for i = 1:length(qP),hold on
    plot([i i] + 1/8,qP(i) + [-1 1]*cq(i),'LineWidth',4,'color','r')
end, hold off

% Compare
%==========================================================================
spm_figure('GetWin','Figure 3: Comparison of integration schemes');

% hyperparameters
%--------------------------------------------------------------------------
qL    = spm_vec({LAP.qH.h LAP.qH.g});
qD    = spm_vec({DEM.qH.h DEM.qH.g});
vL    = spm_vec({LAP.qH.V LAP.qH.W});
vD    = spm_vec({DEM.qH.V DEM.qH.W});
qh    = log([G(1).V; G(1).W]);

subplot(2,2,1)
bar([qh qL qD])
axis square
legend('true','parallel','standard')
title('log-precisions','FontSize',16)

cq    = 1.64*sqrt(vL);
for i = 1:length(qL),hold on
    plot([i i] + 0,qL(i) + [-1 1]*cq(i),'LineWidth',4,'color','r')
end, hold off

cq    = 1.64*sqrt(vD);
for i = 1:length(qD),hold on
    plot([i i] + 1/4,qD(i) + [-1 1]*cq(i),'LineWidth',4,'color','r')
end, hold off

% Log-evidence
%--------------------------------------------------------------------------
subplot(2,2,2)
nL   = length(LAP.F);
nD   = length(DEM.F);
plot(1:nL,LAP.F,1:nD,DEM.F)
axis square
legend('parallel (F)','standard (F)')
title('log-evidence ','FontSize',16)
xlabel('iteration','FontSize',12)

author: spm


text:
function hermite_polynomial_test16 ( )

%*****************************************************************************80
%
%% HERMITE_POLYNOMIAL_TEST16 tests Hermite projection.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 December 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_POLYNOMIAL_TEST16:\n' );
  fprintf ( 1, '  As a sanity check, make sure that the projection of\n' );
  fprintf ( 1, '  He(i,x) is 1 for the i-th component and zero for all others.\n' );

  n = 3;

  [ x, w ] = he_quadrature_rule ( n + 1 );

  phi = hen_polynomial_value ( n + 1, n, x );

  for j = 0 : n

    c = zeros ( n + 1, 1 );

    f_vec(1:n+1,1) = phi ( 1 : n + 1, j + 1 );

    for i = 1 : n + 1
      phiw(i,1:n+1) = phi(i,1:n+1) * w(i);
    end

    c = phiw' * f_vec;
    title = sprintf ( '  Coefficients for He(%d,x)', j );

    r8vec_print ( n + 1, c, title );

  end

  return
end


text:
function table2 = r8mat_border_add ( m, n, table )

%*****************************************************************************80
%
%% R8MAT_BORDER_ADD adds a "border" to an R8MAT.
%
%  Discussion:
%
%    We suppose the input data gives values of a quantity on nodes
%    in the interior of a 2D grid, and we wish to create a new table
%    with additional positions for the nodes that would be on the
%    border of the 2D grid.
%
%                  0 0 0 0 0 0
%      * * * *     0 * * * * 0
%      * * * * --> 0 * * * * 0
%      * * * *     0 * * * * 0
%                  0 0 0 0 0 0
%
%    The illustration suggests the situation in which a 3 by 4 array
%    is input, and a 5 by 6 array is to be output.
%
%    The old data is shifted to its correct positions in the new array.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the table data.
%
%    Output, real TABLE2(M+2,N+2), the augmented table data.
%
  table2 = zeros(m+2,n+2);

  table2(1,1:n+2) = 0.0;
  table2(m+2,1:n+2) = 0.0;
  table2(2:m+1,1) = 0.0;
  table2(2:m+1,n+2) = 0.0;

  table2(2:m+1,2:n+1) = table(1:m,1:n);

  return
end


text:
function [cm3] = ml2cm3(ml)
% Convert volume from milliliters to cubic centimeters.
% Chad Greene 2012
cm3 = ml;


text:
function [ n_data, mu, sigma, x, fx ] = cauchy_cdf_values ( n_data )

%*****************************************************************************80
%
%% CAUCHY_CDF_VALUES returns some values of the Cauchy CDF.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      Needs["Statistics`ContinuousDistributions`"]
%      dist = CauchyDistribution [ mu, sigma ]
%      CDF [ dist, x ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real MU, the mean of the distribution.
%
%    Output, real SIGMA, the standard deviation of the distribution.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 12;

  fx_vec = [ ...
     0.5000000000000000E+00, ...
     0.8524163823495667E+00, ...
     0.9220208696226307E+00, ...
     0.9474315432887466E+00, ...
     0.6475836176504333E+00, ...
     0.6024163823495667E+00, ...
     0.5779791303773693E+00, ...
     0.5628329581890012E+00, ...
     0.6475836176504333E+00, ...
     0.5000000000000000E+00, ...
     0.3524163823495667E+00, ...
     0.2500000000000000E+00 ];

  mu_vec = [ ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.3000000000000000E+01, ...
     0.4000000000000000E+01, ...
     0.5000000000000000E+01 ];

  sigma_vec = [ ...
     0.5000000000000000E+00, ...
     0.5000000000000000E+00, ...
     0.5000000000000000E+00, ...
     0.5000000000000000E+00, ...
     0.2000000000000000E+01, ...
     0.3000000000000000E+01, ...
     0.4000000000000000E+01, ...
     0.5000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.2000000000000000E+01 ];

  x_vec = [ ...
     0.1000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.3000000000000000E+01, ...
     0.4000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.3000000000000000E+01, ...
     0.3000000000000000E+01, ...
     0.3000000000000000E+01, ...
     0.3000000000000000E+01 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    mu = 0.0;
    sigma = 0.0;
    x = 0.0;
    fx = 0.0;
  else
    mu = mu_vec(n_data);
    sigma = sigma_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function table2 = r8mat_border_add ( m, n, table )

%*****************************************************************************80
%
%% R8MAT_BORDER_ADD adds a "border" to an R8MAT.
%
%  Discussion:
%
%    An R8MAT is an array of R8's.
%
%    We suppose the input data gives values of a quantity on nodes
%    in the interior of a 2D grid, and we wish to create a new table
%    with additional positions for the nodes that would be on the
%    border of the 2D grid.
%
%                  0 0 0 0 0 0
%      * * * *     0 * * * * 0
%      * * * * --> 0 * * * * 0
%      * * * *     0 * * * * 0
%                  0 0 0 0 0 0
%
%    The illustration suggests the situation in which a 3 by 4 array
%    is input, and a 5 by 6 array is to be output.
%
%    The old data is shifted to its correct positions in the new array.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the table data.
%
%    Output, real TABLE2(M+2,N+2), the augmented table data.
%
  table2 = zeros(m+2,n+2);

  table2(1,1:n+2) = 0.0;
  table2(m+2,1:n+2) = 0.0;
  table2(2:m+1,1) = 0.0;
  table2(2:m+1,n+2) = 0.0;

  table2(2:m+1,2:n+1) = table(1:m,1:n);

  return
end


function visible = visiblevertices(FV, R)

V2 = (R*FV.vertices')';

% Store the vertex depths for z-buffering
Z = V2(:, 3);

% Compute the projected vertices in the image plane
UV(:, 1)    = V2(:, 1) ;	% orthographic projection for x
UV(:, 2)    = V2(:, 2) ;	% orthographic projection for y

% Transform to the pixel plane (the axes remain switched)
UV(:, 1)    = UV(:, 1)-min(UV(:, 1));
UV(:, 2)    = UV(:, 2)-min(UV(:, 2));
UV          = UV + 1;

UV=(UV./max(UV(:))).*1000;
width=1000;
height=1000;

% Get the triangle vertices
v1      = FV.faces(:, 1);
v2      = FV.faces(:, 2);
v3      = FV.faces(:, 3);
Nfaces  = size(FV.faces, 1);

% Compute bounding boxes for the projected triangles
x       = [UV(v1, 1), UV(v2, 1), UV(v3, 1)];
y       = [UV(v1, 2), UV(v2, 2), UV(v3, 2)];
minx    = ceil (min(x, [], 2));
maxx    = floor(max(x, [], 2));
miny    = ceil (min(y, [], 2));
maxy    = floor(max(y, [], 2));

clear x y

% Frustum culling
minx    = max(1,        minx);
maxx    = min(width,    maxx);
miny    = max(1,        miny);
maxy    = min(height,   maxy);

% Construct the pixel grid (can speed up by precomputing if shared among the images)
[rows, cols] = meshgrid(1: width, 1: height);

% Initialize the depth-, face- and weight-buffers
zbuffer     = -inf(height, width);
fbuffer     = zeros(height, width);
wbuffer1    = NaN(height, width);
wbuffer2    = NaN(height, width);
wbuffer3    = NaN(height, width);

% For each triangle (can speed up by comparing the triangle depths to the z-buffer and priorly sorting the triangles by increasing depth)
for i = 1: Nfaces

    % If some pixels lie in the bounding box
    if minx(i) <= maxx(i) && miny(i) <= maxy(i)

        % Get the pixels lying in the bounding box
        px = rows(miny(i): maxy(i), minx(i): maxx(i));
        py = cols(miny(i): maxy(i), minx(i): maxx(i));
        px = px(:);
        py = py(:);

        % Compute the edge vectors
        e0 = UV(v1(i), :);
        e1 = UV(v2(i), :) - e0;
        e2 = UV(v3(i), :) - e0;

        % Compute the barycentric coordinates (can speed up by first computing and testing a solely)
        det     = e1(1) * e2(2) - e1(2) * e2(1);
        tmpx    = px - e0(1);
        tmpy    = py - e0(2);
        a       = (tmpx * e2(2) - tmpy * e2(1)) / det;
        b       = (tmpy * e1(1) - tmpx * e1(2)) / det;

        % Test whether the pixels lie in the triangle
        test = a >= 0 & b >= 0 & a + b <= 1;

        % If some pixels lie in the triangle
        if any(test)

            % Get the pixels lying in the triangle
            px = px(test);
            py = py(test);

            % Interpolate the triangle depth for each pixel
            w2 = a(test);
            w3 = b(test);
            w1 = 1 - w2 - w3;
            pz = Z(v1(i)) * w1 + Z(v2(i)) * w2 + Z(v3(i)) * w3;

            % For each pixel lying in the triangle
            for j = 1: length(pz)

                % Frustum culling
%                if pz(j) <= -near && pz(j) >= -far

                    % Update the depth-, face- and weight-buffers
                    if pz(j) > zbuffer(py(j), px(j))
                        zbuffer(py(j), px(j))   = pz(j);
                        fbuffer(py(j), px(j))   = i;
                        wbuffer1(py(j), px(j))  = w1(j);
                        wbuffer2(py(j), px(j))  = w2(j);
                        wbuffer3(py(j), px(j))  = w3(j);
                    end

%                end

            end

        end

    end

end

%figure; imshow(zbuffer,[])

clear UV Z zbuffer px py pz minx maxx miny maxy

% Get the vertices to render
test    = fbuffer ~= 0;
f       = unique(fbuffer(test));
v       = unique([v1(f); v2(f); v3(f)]);
f       = find(any(ismember(FV.faces, v), 2));
Nfaces  = length(f);

visible = v;

return

% Compute the edge vectors
e1s = V2(v2(f), :) - V2(v1(f), :);
e2s = V2(v3(f), :) - V2(v1(f), :);
e3s = V2(v2(f), :) - V2(v3(f), :);

clear V2

% Normalize the edge vectors
e1s_norm = e1s ./ repmat(sqrt(sum(e1s.^2, 2)), 1, 3);
e2s_norm = e2s ./ repmat(sqrt(sum(e2s.^2, 2)), 1, 3);
e3s_norm = e3s ./ repmat(sqrt(sum(e3s.^2, 2)), 1, 3);

% Compute the angles
angles(:, 1) = acos(sum(e1s_norm .* e2s_norm, 2));
angles(:, 2) = acos(sum(e3s_norm .* e1s_norm, 2));
angles(:, 3) = pi - (angles(:, 1) + angles(:, 2));

% Compute the triangle weighted normals
triangle_normals    = cross(e1s, e3s, 2);
w1_triangle_normals = triangle_normals .* repmat(angles(:, 1), 1, 3);
w2_triangle_normals = triangle_normals .* repmat(angles(:, 2), 1, 3);
w3_triangle_normals = triangle_normals .* repmat(angles(:, 3), 1, 3);

clear e1s e2s e3s e1s_norm e2s_norm e3s_norm angles triangle_normals

% Initialize the vertex normals
normals = zeros(Nvertices, 3);

% Update the vertex normals
for i = 1: Nfaces
    normals(v1(f(i)), :) = normals(v1(f(i)), :) + w1_triangle_normals(i, :);
    normals(v2(f(i)), :) = normals(v2(f(i)), :) + w2_triangle_normals(i, :);
    normals(v3(f(i)), :) = normals(v3(f(i)), :) + w3_triangle_normals(i, :);
end

clear w1_triangle_normals w2_triangle_normals w3_triangle_normals

% Normalize the vertex normals
normals = normals(v, :);
normals = normals ./ repmat(sqrt(sum(normals.^2, 2)), 1, 3);

% Self-occlusions
test = normals(:, 3) >= 0;

% Interpolate the z-buffer at the projected vertices
vZ          = inf(size(FV.vertices, 1), 1);
vZ(test)    = interp2(rows, cols, zbuffer, UV(test, 1), UV(test, 2), '*linear');

% Fix the border
tmp_test            = vZ(test) == -inf;
vZ(test(tmp_test))	= interp2(rows, cols, zbuffer, UV(test(tmp_test), 1), UV(test(tmp_test), 2), '*nearest');

% Compute the vertex visibility index in terms of relative depth compared to the z-buffer
test(test) = abs(vZ(test) - Z(test)) <= 0.05 * (max(Z) - min(Z));

clear ambient1 diffuse1 specular1 ambient2 diffuse2 specular2 ambient3 diffuse3 specular3 v

% Initialize the image
im1 = NaN(height, width);
im2 = NaN(height, width);
im3 = NaN(height, width);

% Rasterize the image
v1          = v1(fbuffer(test));
v2          = v2(fbuffer(test));
v3          = v3(fbuffer(test));
w1          = wbuffer1(test);
w2          = wbuffer2(test);
w3          = wbuffer3(test);
im1(test)   = w1 .* texture1(v1) + w2 .* texture1(v2) + w3 .* texture1(v3);
im2(test)   = w1 .* texture2(v1) + w2 .* texture2(v2) + w3 .* texture2(v3);
im3(test)   = w1 .* texture3(v1) + w2 .* texture3(v2) + w3 .* texture3(v3);

clear v1 v2 v3 w1 w2 w3

% Reshape the image
im = cat(3, im1, im2, im3);

clear im1 im2 im3 test

im = flipud(fliplr(im));

end


text:
%% Gauss Laguerre Quadrature
% solution for class exercise:

clear
clc
close all

%% Aproximate solution
% Tau(m) is a function of m.

m=5.5   % Exponent value of f(x).
n=9;    % number of points used to compute approximate solution.

[x w]=GaussLaguerre(n,0); % built in function to generate weight and points

% Initalizing row vectors:
l=length(x);
f=zeros(1,l);
t=zeros(1,l);

for i=1:l
    f(i)=x(i)^(m-1);
    t(i)=w(i)*f(i);
end

Gamma=sum(t)

%% Exact Solution
% Gamma=(m-1)!
% Carefull! m can only be an integer number!

Gamma2=factorial(floor(m)-1)
% I'm using a round to the floor in case m is a rational number.

%% Exact solution of the Gamma function
% this time using the Matlab's gamma function.

Gamma3=gamma(m)

repo: Aero-matlab


text:
% Local Linear Embedding
% This function implements the weight computation defined in
% Sam T. Roweis, Lawrence K. Saul, "Nonlinear Dimensionality
% Reduction by Local Linear Embedding", Science, 2000.
% 'w' is the weights for representing the row-vector 'pt' in terms
% of the dimensions x neighborCount matrix 'neighbors'.
% 'conditionerMult' is the multiplier of the identity matrix added
% to the neighborhood correlation matrix before inversion.

function w = localLinearEmbedding(pt, neighbors, conditionerMult)
    % each column of neighbors represent a neighbor, each row a dimension
    % pt is a row vector
    corr = neighbors' * neighbors + conditionerMult * eye(size(neighbors, 2));
    ptDotN = neighbors' * pt;
    alpha = 1 - sum(corr \ ptDotN);
    beta = sum(corr \ ones(size(corr, 1), 1)); % sum of elements of inv(corr)
    lagrangeMult = alpha / beta;
    w = corr \ (ptDotN + lagrangeMult);
end

author: yaksoy

repo: AffinityBasedMattingToolbox


text:
function determ = sweet3_determinant ( )

%*****************************************************************************80
%
%% SWEET3_DETERMINANT returns the determinant of the SWEET3 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real DETERM, the determinant.
%
  determ = - 5.4056067E+07;

  return
end


text:
function [ imdb, roidb ] = imdb_sample_frames(imdb, roidb,max_samples_from_cls, offset )
    if nargin < 4
      offset = 0;
    end
    valid = zeros(size(imdb.image_ids));
    if iscell(max_samples_from_cls)
      [C,ia,ib] = intersect(imdb.image_ids, max_samples_from_cls);
      valid = ia;
    else
      nFrames = imdb.num_frames;
      count=0;
      for i=1:numel(nFrames)
        % shuffled_frames{i} = randperm(nFrames(i),min(nFrames(i),max_samples_from_cls)) + count;
        shuffled_frames{i} = round(linspace(1+offset,nFrames(i)-offset,min(nFrames(i)-2*offset,max_samples_from_cls))) + count;
        count = count + nFrames(i);
      end
      valid = [shuffled_frames{:}];
    end

     % clear images without any object from target_imdb

    roidb.rois = roidb.rois(valid);
    imdb.is_blacklisted = imdb.is_blacklisted(valid);
    imdb.sizes = imdb.sizes(valid,:);
    if isfield(imdb,'flip_from'),  imdb.flip_from = imdb.flip_from(valid); end
    imdb.image_ids = imdb.image_ids(valid);
    imdb.image_at = @(i) ...
    fullfile(imdb.image_dir, [imdb.image_ids{i} '.' imdb.extension]);
    if isfield(imdb,'vid_id'),    imdb.vid_id = imdb.vid_id(valid); end

end

author: feichtenhofer

repo: Detect-Track


text:
function [hostName, databaseName, username, password, customPort] = SQL_ShowConnSettings(doDisplay)
% SQL_ShowConnSettings  Display the database connection settings,
% read in from the sql_settings.conf file.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
% Check Inputs:
% ------------------------------------------------------------------------------

if nargin < 1 || isempty(doDisplay)
    doDisplay = true;
end

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

theConfigFile = which('sql_settings.conf');
if isempty(theConfigFile)
    % no sql_settings.conf found
    error(['No sql_settings.conf file found.\n  ' ...
    'Please create sql_settings.conf file using SQL_CreateDatabase, SQL_ChangeDatabase,' ...
        ' or generate one yourself (cf. Documentation).'])
else
    fid = fopen(theConfigFile);
end

% Read the file:
s = textscan(fid,'%s%s%s%s%u','Delimiter',',','CommentStyle','%');
fclose(fid);

% Interpret the output:
hostName = s{1}{1};
databaseName = s{2}{1};
username = s{3}{1};
password = s{4}{1};
customPort = s{5};

if isempty(customPort) || customPort==0
    customPort = 3306; % default port
end

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

if doDisplay
    fprintf(1,'Connecting to %s at %s as %s (through port %u)\n',...
                    databaseName,hostName,username,customPort);
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function pick = nms(boxes, overlap)
% Non-maximum suppression.
%   pick = nms(boxes, overlap)
%
%   Greedily select high-scoring detections and skip detections that are
%   significantly covered by a previously selected detection.
%
% Return value
%   pick      Indices of locally maximal detections
%
% Arguments
%   boxes     Detection bounding boxes (see pascal_test.m)
%   overlap   Overlap threshold for suppression
%             For a selected box Bi, all boxes Bj that are covered by
%             more than overlap are suppressed. Note that 'covered' is
%             is |Bi \cap Bj| / |Bj|, not the PASCAL intersection over
%             union measure.

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
% Copyright (C) 2008, 2009, 2010 Pedro Felzenszwalb, Ross Girshick
% Copyright (C) 2007 Pedro Felzenszwalb, Deva Ramanan
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

if isempty(boxes)
  pick = [];
else
  x1 = boxes(:,1);
  y1 = boxes(:,2);
  x2 = boxes(:,3);
  y2 = boxes(:,4);
  s = boxes(:,end);
  area = (x2-x1+1) .* (y2-y1+1);

  [vals, I] = sort(s);
  pick = [];
  while ~isempty(I)
    last = length(I);
    i = I(last);
    pick = [pick; i];
    suppress = [last];
    for pos = 1:last-1
      j = I(pos);
      xx1 = max(x1(i), x1(j));
      yy1 = max(y1(i), y1(j));
      xx2 = min(x2(i), x2(j));
      yy2 = min(y2(i), y2(j));
      w = xx2-xx1+1;
      h = yy2-yy1+1;
      if w > 0 && h > 0
        % compute overlap
        o = w * h / area(j);
        if o > overlap
          suppress = [suppress; pos];
        end
      end
    end
    I(suppress) = [];
  end
end

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function varargout = HELP(varargin)
% HELP M-file for HELP.fig
%      HELP, by itself, creates a new HELP or raises the existing
%      singleton*.
%
%      H = HELP returns the handle to a new HELP or the handle to
%      the existing singleton*.
%
%      HELP('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in HELP.M with the given input arguments.
%
%      HELP('Property','Value',...) creates a new HELP or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before HELP_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to HELP_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help HELP

% Last Modified by GUIDE v2.5 06-Feb-2008 22:52:28

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @HELP_OpeningFcn, ...
                   'gui_OutputFcn',  @HELP_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before HELP is made visible.
function HELP_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to HELP (see VARARGIN)

% Choose default command line output for HELP
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes HELP wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = HELP_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


text:
%%  Introduction to *wavelet_layer_2d*
% *wavelet_layer_2d* is computing the wavelet transform coefficient from
% the coefficient of the previous layer. It should be applied iteratively
% to an input signal.
%
%% Usage
% [A, V] = *wavelet_layer_2d*(U, filters, options), documentation is given in
% <matlab:doc('wavelet_layer_2d') wavelet_layer_2d>
%
%% Description
% It is possible to create some wavelet filters with wavelet_factory_2d for
% instance. The filters size have to be adapted to the size of the input
% signal $x$.
clear; close all;

x = mandrill;

% Create $ U[\empty]x $
U{1}.signal{1} = x;
U{1}.meta.j = zeros(0,1);
U{1}.meta.q = zeros(0,1);
U{1}.meta.resolution=0;
filters = morlet_filter_bank_2d(size(x));

% A corresponds to the output scattering coefficients, and V to the wavelet
% coefficients.
[A, V] = wavelet_layer_2d(U{1}, filters);

colormap gray
subplot(121)
imagesc(real(V.signal{1}))
axis off
title('Real part of the first wavelet transform coefficient');
subplot(122)
imagesc(imag(V.signal{1}))
axis off
title('Imaginary part of the first wavelet transform coefficient');

%% Options
% The options are the same as in *wavelet_2d*.


text:
function A = projection_precon_mnls(X, B, A)

    s = svd(B'*B);
    L = max(s);
    %eigen_values = eig(B'*B);
    %L = max(eigen_values);

    grad = - X * B + A*(B' * B);
    A = A - 1/L * grad / (B'*B);
    A = max(A, 0);
end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function varargout = atan(varargin)
%ATAN (overloaded)

switch class(varargin{1})

    case 'sdpvar'
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'

        operator = struct('convexity','none','monotonicity','increasing','definiteness','none','model','callback');
        operator.convexhull = @convexhull;
        operator.bounds = @bounds;
        operator.derivative = @(x)((1+x.^2).^-1);
        operator.inverse = @(x)(tan(x));
        operator.range = [-pi/2 pi/2];

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error('SDPVAR/ATAN called with CHAR argument?');
end

function [L,U] = bounds(xL,xU)
L = atan(xL);
U = atan(xU);

function [Ax, Ay, b] = convexhull(xL,xU)
fL = atan(xL);
fU = atan(xU);
dfL = 1/(1+xL^2);
dfU = 1/(1+xU^2);
if xL >= 0
    % Concave region
    [Ax,Ay,b] = convexhullConcave(xL,xU,fL,fU,dfL,dfU);
elseif xU <= 0
    % Convex region
    [Ax,Ay,b] = convexhullConvex(xL,xU,fL,fU,dfL,dfU);
else
    % Changes convexity. We're lazy and let YALMIP sample instead
    Ax = [];
    Ay = [];
    b = [];
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function Y = utl_picktimes(X,wnd)
% Average the input array within the given index ranges, for each row and plane.
% Result = utl_picktimes(Data, Windows)
%
% In:
%   Data    : [Channels x Samples x Trials] data array
%   Windows : [Ranges x 2] array specifying the beginning and end of each range, in samples
%
% Out:
%   Result  : averaged sub-ranges of the original data array; sized [Channels x Ranges x Trials]
%
% Examples:
%   % for a given 3d data array, average intervals of 20:30 samples, 50:100 samples, and 100:200 samples
%   % for each epoch and each channel
%   averages = utl_picktimes(EEG.data,[20 30;50 100;100 200])
%
% See also:
%   set_picktimes
%
%                               Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                               2010-04-20

[C,S,T] = size(X);

if any(wnd(:,1) > wnd(:,2))
    error('Some of the windows are malformed: %s',hlp_tostring(wnd(wnd(:,1) > wnd(:,2),:))); end
if any(wnd(:,1)>S | wnd(:,2)<1)
    error('Some windows are entirely outside the data dimensions (%i): %s',S,hlp_tostring(wnd(wnd(1,:)>S|wnd(2,:)<1,:))); end
if any(wnd(:)>S | wnd(:)<1)
    disp_once('WARNING: some windows passed to utl_picktimes exceed the data dimensions; windows trimmed.'); end

% fix the windows, if necessary
wnd = min(max(wnd,1),S);
wnd(:,2) = wnd(:,2)+1;
W = size(wnd,1);

% calc offset & coverage values for first & last sample
fo = floor(wnd(:,1));
fc = min(fo+1,wnd(:,2)) - wnd(:,1);
lo = ceil(wnd(:,2)-1);
lc = wnd(:,2) - max(lo,wnd(:,1));
% calc fully overlapped sample range, inverse length
for r=1:W
    full{r} = fo(r)+1:lo(r)-1; end %#ok<AGROW>
ilen = 1./(fc+cellfun('length',full)'+lc);

Y = zeros(C,W,T);
% accumulate for every range...
for r=1:W
    Y(:,r,:) = (X(:,fo(r),:)*fc(r)+sum(X(:,full{r},:),2)+X(:,lo(r),:)*lc(r))*ilen(r); end

repo: REST


text:
function DPMpcStandAlongJingOneBeam3(leak, spectrum_File, planC_File, nhist, OutputError, whichBeam, PBMaxWidth, gradsense, scatterThreshold, batch)
% written by KZ, JC
% Calculate doseV dose3D directly from DPM, without compressing it.
% Use generateDPMdose.m instead of DPMInfluenceJing.m
% Calculate IM struct, call DPM, get
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

% Usage: IIMCalc = calculateBeamIM(0.018, 'DPM_10beams_planC', 1, 1, 5)
%   leak = 0.018
%   spectrum_File = file name of photon spectrum. eg. photon6MV.spectrum
%   planC_File = file name contained planC, eg. 'DPM_10beams_planC'
%   indexBeam = Beam Index: 1, 2, ...
%   imin = Beamlet index
%   imax = Beamlet index

%   eval(['load Beam',num2str(indexBeam),'_w']);
%   where, Beam1_w.mat (Beam2_w.mat) need to be in the current directory. It's the weight of
%   the beamlets for each beam.

% JC Feb 27 06
% Call in DPMInfluence.m; use sum(10*clock) as the random seeds input to DPM.

% JC. Aug 10, 2005, Need to make planC read from the disk, since in the
% stand-alone mode, no matlab run is expected.

% open planC_File
% In matlab release 14, "eval" can be compiled.
Error = []; %Initialize Error to NULL.
load(planC_File);
%eval(['load ',planC_File]);
% This planC_File should contains planC and stateS
%leak = 0.018; %leakage in per *100

if ischar(leak)
    leak = str2num(leak)
end

if ischar(nhist)
    nhist = str2num(nhist)
end

if ischar(OutputError)
    OutputError = str2num(OutputError)
end

if ischar(whichBeam)
    whichBeam = str2num(whichBeam)
end

if ischar(PBMaxWidth)
    PBMaxWidth = str2num(PBMaxWidth)
end

if ischar(gradsense)
    gradsense = str2num(gradsense)
end

if ischar(scatterThreshold)
    scatterThreshold = str2num(scatterThreshold)
end

if ischar(batch)
    batch = str2num(batch)
end

currentDir = cd;

%for indexBeam = 1 : planC{7}.FractionGroupSequence.Item_1.NumberOfBeams
for indexBeam = whichBeam
    bf = planC{7}.FractionGroupSequence.Item_1.ReferencedBeamSequence.(['Item_' num2str(indexBeam)]);
    %if (bf.BeamDose ~= 0 & indexBeam ~= 2)
    if (bf.BeamDose ~= 0)
        bs = planC{7}.BeamSequence.(['Item_' num2str(indexBeam)]);

        LS = getDICOMLeafPositions(bs)

        [inflMap, xV, yV, colDividerXCoord, rowDividerYCoord, rowLeafPositions] = getLSInfluenceMapFactor(LS,leak,bs.BeamNumber);

        gA = bs.ControlPointSequence.Item_1.GantryAngle;
        iC = bs.ControlPointSequence.Item_1.IsocenterPosition;
        IMBase.beams(1).gantryAngle = gA;

        if ~isfield(planC{7}.PatientSetupSequence,(['Item_' num2str(indexBeam)]))
            position = {planC{7}.PatientSetupSequence.(['Item_' num2str(1)]).PatientPosition};
        else
            position = {planC{7}.PatientSetupSequence.(['Item_' num2str(indexBeam)]).PatientPosition};
        end

        if strcmpi(position, 'HFP')
            IMBase.beams(1).isocenter.x = iC(1)/10;
            IMBase.beams(1).isocenter.y = iC(2)/10;
            IMBase.beams(1).isocenter.z = -iC(3)/10;
        else
            IMBase.beams(1).isocenter.x = iC(1)/10;
            IMBase.beams(1).isocenter.y = -iC(2)/10;
            IMBase.beams(1).isocenter.z = -iC(3)/10;
        end

        IMBase.beams(1).isodistance = bs.SourceAxisDistance/10;

        IMBase.beams(1).xRel = IMBase.beams(1).isodistance * sindeg(IMBase.beams(1).gantryAngle);

        IMBase.beams(1).yRel = IMBase.beams(1).isodistance * cosdeg(IMBase.beams(1).gantryAngle);

        IMBase.beams(1).beamEnergy = bs.ControlPointSequence.Item_1.NominalBeamEnergy;

        %[PBX] = getPBXCoorFromFluence(inflMap, rowLeafPositions, 2, 50, colDividerXCoord, rowDividerYCoord, 10);
        [PBX] = getPBXCoorFromFluence(inflMap, rowLeafPositions, 2, PBMaxWidth, colDividerXCoord, rowDividerYCoord, gradsense);

        [xPosV, yPosV, beamlet_delta_x, beamlet_delta_y, w_field] = getPBVectors(PBX,rowLeafPositions,colDividerXCoord,rowDividerYCoord,inflMap);

%        %Visualization
%                 figure;hAxis1 = axes;imagesc([xV(1)+.05 xV(end)-.05], [yV(1)+.05, yV(end)-.05], inflMap);
%                 set(gcf, 'renderer', 'zbuffer')
%                 figure;hAxis2 = axes;hold on;
%                 set(gcf, 'renderer', 'zbuffer')
%                 xL = get(hAxis1, 'xlim');
%                 yL = get(hAxis1, 'ylim');
%                 set(hAxis2, 'xlim', xL);
%                 set(hAxis2, 'ylim', yL);
%                 axis(hAxis2, 'manual');
%             %     w_colors = floor((w_field ./ max(w_field))*255)+1;
%                 set(gcf, 'doublebuffer', 'on');
%                 for i=1:length(xPosV)
%                     patch([xPosV(i) - beamlet_delta_x(i)/2 xPosV(i) - beamlet_delta_x(i)/2 xPosV(i) + beamlet_delta_x(i)/2 xPosV(i) + beamlet_delta_x(i)/2 xPosV(i) - beamlet_delta_x(i)/2], [yPosV(i) - beamlet_delta_y(i)/2 yPosV(i) + beamlet_delta_y(i)/2 yPosV(i) + beamlet_delta_y(i)/2 yPosV(i) - beamlet_delta_y(i)/2 yPosV(i) - beamlet_delta_y(i)/2], w_field(i));
%                 end
%                 axis([hAxis1 hAxis2], 'ij');
%                 kids = get(hAxis2, 'children');
%             %     set(kids, 'edgecolor', 'none');
%                 cMap = colormap('jet');
%                 set(hAxis2, 'color', cMap(1,:));

        xPosV = xPosV/10;
        yPosV = yPosV/10;
        beamlet_delta_x = beamlet_delta_x/10;
        beamlet_delta_y = beamlet_delta_y/10;

        IMBase.goals.PBMargin = 0.5;
        IMBase.goals.structNum = 2;
        IMBase.params.xyDownsampleIndex = 1;
        IMBase.goals.isTarget(1) = 'y';
        IMBase.goals.xySampleRate = 2;
        IMBase.params.numCTSamplePts = 300;

        IMBase.beams(1).zRel = 0;
        IMBase.beams(1).xRel =  IMBase.beams(1).isodistance * sindeg(IMBase.beams(1).gantryAngle);
        IMBase.beams(1).yRel =  IMBase.beams(1).isodistance * cosdeg(IMBase.beams(1).gantryAngle);

        %IMBase.params.algorithm = 'VMC++';
        %IMCalc = IMRTP_MapCheck(IMBase);

        IMBase.params.algorithm = 'DPM';
        IMCalc = IMRTP_MapCheck(IMBase, planC, stateS, xPosV, yPosV, beamlet_delta_x, beamlet_delta_y, gA);
        % JC. 11 Aug, 05 Added more arguments to IMRTP_MapCheck,

        %disp('OK_BIG1'); pause(5);
        % JC. Aug. 3, 2005
        % Add fields to IMCalc, necessary for generateDPMInfluence.m
        IMCalc.beams.couchAngle = 0;
        IMCalc.beams.collimatorAngle = 0;
        IMCalc.params.ScatterMethod = 'exponential';
        IMCalc.params.Scatter.Threshold = scatterThreshold;
        IMCalc.params.Scatter.RandomStep = 30;
        IMCalc.beamlets = struct;

        % JC Jun 06 2006
        % Still need to fill in the structure of IMCalc, but don't need put
        % it back(?), since it's not to be used to calculate dose.
        % Will also need to use 'DoseScale' inside generateDPMdose.m,
        % since DPM output dose on a per partical base.
        % Need to use w_field as a way of telling DPM how to add beamlets
        % up.

        [doseV indV] = generateDPMdose(IMCalc,spectrum_File,nhist, OutputError,planC,stateS, indexBeam, w_field);

       siz = getUniformizedSize(planC);
       dose3D = zeros(siz);

       if (OutputError == 1)
            doseV = doseV(:,1);
            Error_3D = zeros(siz);
            Error_3D(indV') = doseV(:,2);
        else
            disp('Do not output error')
            dose3D(indV') = doseV;
        end

        filename = ['doseV_',num2str(indexBeam),'_', num2str(nhist), '_', num2str(batch)];
        save(filename, 'doseV', 'indV');

        filename = ['dose3D_',num2str(indexBeam),'_', num2str(nhist), '_', num2str(batch)];
        save(filename, 'dose3D');

        % Only Output w_field for the first batch
        if(batch == 1)
            filename = ['w_field',num2str(indexBeam)];
            save(filename, 'w_field');%dose3D = (dose3D/max(dose3D(:)));
        end

        clear IMCalc dose3D w_field

    end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function[p] = gui_erp_plot(p)

% gui_erp_plot - eeg_toolbox plot for ERP data
%
%[p] = erp_plot(p)
%
% creates a GUI interface for ERP measurement and
% topographic mapping functions.
%

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

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

plotfig = figure('Name',p.volt.file,...
  'NumberTitle','off',...
  'UserData',p);

movegui(plotfig,'center');

plot(p.volt.timeArray,p.volt.data);
axis tight;
eeg_plot_metric;

[Xpoint,Ypoint] = eeg_crosshair('init',p);

return

author: PatternRecognition

repo: OpenBMI


text:
function train_id_net_vgg16(varargin)
% -------------------------------------------------------------------------
% Part 4.1: prepare the data
% -------------------------------------------------------------------------

imdb = load('./dataset/Flickr30k-prepare/url_data.mat');
imdb = imdb.imdb;
load('./dataset/Flickr30k-prepare/dense_feature_word2.1.mat');
imdb.charcnn = wordcnn;
%imdb.charmean = mean(imdb.charcnn(:,:,:,imdb.images.set==1),4);
% -------------------------------------------------------------------------
% Part 4.2: initialize a CNN architecture
% -------------------------------------------------------------------------
net = resnet52_new_hope_word_Rankloss();
net.conserveMemory = true;
im_mean = imdb.rgbMean;
net.meta.normalization.averageImage = im_mean;
%net.meta.normalization.charmean = imdb.charmean;
% -------------------------------------------------------------------------
% Part 4.3: train and evaluate the CNN

% -------------------------------------------------------------------------
opts.train.averageImage = net.meta.normalization.averageImage;
opts.train.batchSize = 32;
opts.train.continue = true;
opts.train.gpus = 3;
opts.train.prefetch = false ;
opts.train.nesterovUpdate = true ;
opts.train.expDir = './data/res52_batch32_Rankloss_2:1:0.1_margin0.5_img0.75_shift_hard_256_no_ranking';
opts.train.derOutputs = {'objective_f',0,'objective_img',1,'objective_txt',0.1} ;
opts.train.candidate = 256;
%opts.train.gamma = 0.9;
opts.train.momentum = 0.9;
%opts.train.constraint = 100;
opts.train.learningRate = [0.1*ones(1,40),0.01*ones(1,20)] ;
opts.train.weightDecay = 0.0001;
opts.train.numEpochs = numel(opts.train.learningRate) ;
[opts, ~] = vl_argparse(opts.train, varargin) ;
% Call training function in MatConvNet
[net,info] = cnn_train_dag_batchsize_net(net, imdb, @getBatch,opts) ;

% --------------------------------------------------------------------
function inputs = getBatch(imdb,batch,net,opts)
% --------------------------------------------------------------------
%-- img data
im_url = imdb.images.data(batch) ;
im = vl_imreadjpeg(im_url,'Pack','Resize',[224,224],'Flip',...
    'CropLocation','random','CropSize',[0.8,1],...
    'Interpolation', 'bicubic','NumThreads',16,... %'Brightness', double(0.1*imdb.rgbCovariance),...
    'SubtractAverage',imdb.rgbMean,...
    'CropAnisotropy',[3/4,4/3]);
oim = gpuArray(im{1});
label_img =  imdb.images.label(batch);

%----- first half img feature
net.mode = 'test' ;
%net.vars(net.getVarIndex('fc2bn_n')).value = [];
%net.vars(net.getVarIndex('data')).value = [];
%net.vars(net.getVarIndex('data2')).value = [];
%net.vars(net.getVarIndex('label_img')).value = [];
%net.vars(net.getVarIndex('label_txt')).value = [];
reset(net);
batchsize = numel(batch);
half = batchsize/2;
f_img = getFeature2(net,oim(:,:,:,1:half),[],'data','fc1bn_n');
f_img = reshape(f_img,[],half)';

%------ rand txt input
ind = imdb.images.set==1;  % training set
ind(batch) = 0;  % remove the same label
ind = find(ind==1);  % change 0-1 to real index
ind = (ind-1)*5 + randi(5,numel(ind),1); % img index to char index
rand_ind = ind(randi(numel(ind),opts.candidate,1)); % random select 256 sample

txt = single(imdb.charcnn(:,rand_ind));
txtinput_rand = zeros(1,32,20074,opts.candidate,'single');
for i=1:opts.candidate
    len = sum(txt(:,i)>0);
    location = randi(33-len);
    for j=1:len
        v = txt(j,i);
        txtinput_rand(1,location,v,i)=1;
        location = location+1;
    end
end

%---txt feature
net.vars(net.getVarIndex('fc1bn_n')).value = [];
f_txt = getFeature2(net,gpuArray(txtinput_rand),[],'data2','fc2bn_n');
size4 = size(f_txt,4);
f_txt = reshape(f_txt,[],size4)';

txt_branch = bsxfun(@times,(batch-1),5) + randi(5,batchsize,1); %correct sample
txt = single(imdb.charcnn(:,txt_branch));
txtinput = zeros(1,32,20074,batchsize,'single');
% correct half
for i=1:half
    len = sum(txt(:,i)>0);
    location = randi(33-len);
    for j=1:len
       v = txt(j,i);
       txtinput(1,location,v,i)=1;
       location = location+1;
    end
end
%--- select hard txt sample
for i= half+1:batchsize
    [~,max_ind] = max(f_img(i-half,:)*f_txt'); %1x2048 2048x256
    txt_branch(i) = rand_ind(max_ind);
    txtinput(:,:,:,i) = txtinput_rand(:,:,:,max_ind);
end
label_txt =  imdb.images.label(floor((txt_branch-1)/5)+1);

%----- reset net
net.vars(net.getVarIndex('fc1bn_n')).value = [];
net.vars(net.getVarIndex('fc2bn_n')).value = [];
net.mode = 'normal' ;

%--
inputs = {'data',oim,'data2',gpuArray(txtinput),'label_img',label_img,'label_txt',label_txt};

author: layumi

repo: Image-Text-Embedding