text:
function varargout = buildRun1GUI(varargin)
% function varargout = buildRun1GUI(varargin)
% Purpose: There are three methods for which the GUI can be used to generate
%          the user-defined GUI:
%          1) User provide menu definition thru a GUI interactive input
%             gui4gui.fig. By selecting the GUI generation button, the
%             collected menu definition data will be passed to this
%             function to generate the user GUI.
%          2) User enters menu definition directly into a file. Then, the
%             user manually run this program to generate the user GUI.
%          3) Combination of the above two methods.
% varargin: Caller provides menu definition data as a structure GID.
%           Only varargin{1} is relevant, any additional input will be
%           ignored. The caller is usually gui4gui. It collects GID
%           interactively (or semi-interactively) and then calls genGUI to
%           generate the application GUI.
% If nargin == 0, then this GUI assumes that the menu definition data has
% already been generated and saved in gui4gui_input.mat. It uses this
% file to generate the GUI without going thru gui4gui
% Date created: May 21, 2009
% Kadin Tseng, SCV, Boston University
version = '1.0';

if nargin > 0 & isstruct(varargin{1})
  GID = varargin{1};
else
% loads the data structure that defines the menus
% for the MODE code, there are two situations where this GUI generator might
% be used: Real time task or Fixed duration task. They are both very similar
% in the input requirement but nevertheless requires distinctive m-files
%  varargin{1}
  load run1_inputfile;
%  eval('load varargin{1}');
end

% additional GUI settings:
% fontsize and fontweight for the RUN button
% The title fontsize is 2 sizes larger than the basic user setting
fontsize_title = GID.fontsize + 2;
fontweight_title = GID.fontweight;
fontsize_run = GID.fontsize + 2;
fontweight_run = GID.fontweight;

fid = fopen([GID.guiName '.m'],'w+');

fprintf(fid,'function %s(varargin)\n',GID.guiName);
fprintf(fid,'fh = figure(''Visible'',''off'', ...\n');
fprintf(fid,'            ''HandleVisibility'',''callback'', ...\n');
fprintf(fid,'            ''NumberTitle'',''off'', ...  %% turns off MATLAB GUI window title\n');
fprintf(fid,'            ''Name'',''%s'', ...    %% now, define my own\n',GID.guiName);
fprintf(fid,'            ''Units'',''normalized'', ...\n');
fprintf(fid,'            ''Position'',[%d %d %d %d], ...\n',GID.position);
fprintf(fid,'            ''Color'',''%s'', ...   %% match bg of MODE_new.jpg\n',GID.bgcolor);
fprintf(fid,'            ''Resize'',''%s'');\n',GID.resize);

fprintf(fid,'set(fh,''MenuBar'',''none'');  %% removes MATLAB default menu bar\n');

fprintf(fid,'x0 = 0.1; y0 = 0.85;\n');
fprintf(fid,'x0run = 0.825; y0run = 0.075; %% lower-left corner coord. of the run button\n');
fprintf(fid,'dxst1 = 0.075; dyst1 = 0.05;  %% for static text 1\n');
fprintf(fid,'dxet1 = 0.100; dyet1 = 0.05;  %% for edit text 1\n');
fprintf(fid,'dxst2 = 0.400; dyst2 = 0.10;  %% for static text 2\n');
fprintf(fid,'dxpb1 = 0.050; dypb1 = 0.05;  %% for PDF logo pushbuttons\n');
fprintf(fid,'dxrun = 0.100; dyrun = 0.075; %% for "Run" pushbutton\n');

fprintf(fid,'%% First, create a Heading for the secondary GUI with a static text box\n');
fprintf(fid,'dx = 0.6; dy = 0.125;\n');

fprintf(fid,'st0h = uicontrol(fh, ...\n');
fprintf(fid,'            ''style'',''text'', ...\n');
brkpt = regexp(GID.title,'!');  % check for "!", definde here as line break

%%%% to do: regexprep would be more direct; implementation problem with fprintf

numLines = length(brkpt) + 1;
switch numLines
  case 1    % one-line title
    fprintf(fid,'            ''String'',[''%s''], ...\n',GID.title);
  case 2    % two-line title
    fprintf(fid,'            ''String'',[''%s'' char(10) ''%s''], ...\n',GID.title(1:brkpt-1),GID.title(brkpt+1:end));
  case 3    % three-line title
    fprintf(fid,'            ''String'',[''%s'' char(10) ''%s'' char(10) ''%s''], ...\n', ...
      GID.title(1:brkpt(1)-1),GID.title(brkpt(1)+1:brkpt(2)-1),GID.title(brkpt(2)+1:end));
  otherwise
    disp('Title with more than 3 lines is not anticipated.')
    disp(' "!" is reserved by the GUI builder to force a new line.')
    disp('Make sure not to use it as part of your GUI title text.')
    return
end
fprintf(fid,'            ''Units'',''normalized'', ...\n');
fprintf(fid,'            ''Position'',[x0,y0,dx,dy], ...\n');
fprintf(fid,'            ''FontSize'',%d, ...\n',fontsize_run);
fprintf(fid,'            ''FontWeight'',''%s'');\n',fontweight_run);
fprintf(fid,'\n');
fprintf(fid,'%% Create static box, edit box, static box, logo (that points to a doc; optional)\n');
fprintf(fid,'y0 = y0 - dy*1.05; space = dxst1*0.20;\n');
nlines = size(GID.string,1);
for line=1:nlines
str = GID.string(line,:);
fprintf(fid,'data.str{%d} = ''%s'';\n',line,str{1});
fprintf(fid,'data.default{%d} = %f;\n',line,str2double(str{2}));
fprintf(fid,'  st1h(%d) = uicontrol(fh,''style'',''text'', ...\n',line);
fprintf(fid,'      ''String'',''%s'', ...\n',str{1});
fprintf(fid,'      ''Units'',''normalized'', ...\n');
fprintf(fid,'      ''Position'',[x0,y0,dxst1,dyst1], ...\n');
fprintf(fid,'      ''FontSize'',%d, ...\n',GID.fontsize);
fprintf(fid,'      ''FontWeight'',''%s'');\n',GID.fontweight);
fprintf(fid,'  et1h(%d) = uicontrol(fh,''Style'',''edit'', ...\n',line);
fprintf(fid,'      ''String'',''%s'', ...\n',str{2});
fprintf(fid,'      ''Units'',''normalized'', ...\n');
fprintf(fid,'      ''Position'',[x0+dxst1+space,y0,dxet1,dyet1], ...\n');
fprintf(fid,'      ''FontSize'',%d, ...\n',GID.fontsize);
fprintf(fid,'      ''FontWeight'',''%s'', ...\n','normal');
fprintf(fid,'      ''callback'',{@edit1_Callback ''%s''});\n',str{1});
fprintf(fid,'  st2h(%d) = uicontrol(fh,''style'',''text'', ...\n',line);
fprintf(fid,'      ''String'',''%s'', ...\n',str{3});
fprintf(fid,'      ''Units'',''normalized'', ...\n');
fprintf(fid,'      ''Position'',[x0+dxst1+dxet1+space*2,y0-0.02,dxst2,dyst2], ...\n');
fprintf(fid,'      ''FontSize'',%d, ...\n',GID.fontsize);
fprintf(fid,'      ''FontWeight'',''%s'');\n',GID.fontweight);
fprintf(fid,'  pb1h(%d) = uicontrol(fh,''Style'',''pushbutton'', ...\n',line);
fprintf(fid,'      ''String'',''PDF'', ...\n');
fprintf(fid,'      ''Units'',''normalized'', ...\n');
fprintf(fid,'      ''Position'',[x0+dxst1+dxet1+dxst2+space*3,y0,dxpb1,dypb1], ...\n');
fprintf(fid,'      ''FontSize'',%d, ...\n',GID.fontsize);
fprintf(fid,'      ''FontWeight'',''%s'', ...\n',GID.fontweight);
fprintf(fid,'      ''callback'',{@opendoc1_Callback ''%s''});\n',str{4});
fprintf(fid,'  y0 = y0 - max([dyst1 dyst2 dyet1 dypb1])*1.1;\n');

end   % end of "line" loop

fprintf(fid, ...
    ['%% save the edit text handles to data\n' ...
     'data.et1h = et1h(1:end);\n']);
fprintf(fid,'for i=1:length(data.str)\n');
fprintf(fid,'   eval([''data.'' data.str{i} ''= data.default{i};'']);\n');
fprintf(fid,'end\n');

fprintf(fid,'data.guiName = ''%s'';\n',GID.guiName);
fprintf(fid,'setappdata(fh, ''myData'', data);\n');

% Foonotes and RUN button
fprintf(fid,'%%align(st0h,''HorizontalAlignment'',''center'');\n');
fprintf(fid,'st3h = uicontrol(fh,''style'',''text'', ...\n');
fprintf(fid,'      ''String'',''%s'', ...\n',GID.fn1);
fprintf(fid,'      ''Units'',''normalized'', ...\n');
fprintf(fid,'      ''Position'',[0.05,0.08,0.60,0.05], ...\n');
fprintf(fid,'      ''HorizontalAlignment'',''Left'', ...\n');
fprintf(fid,'      ''foregroundColor'',[1 0 0], ...\n');
fprintf(fid,'      ''FontSize'',%d, ...\n',GID.fontsize);
fprintf(fid,'      ''FontWeight'',''%s'');\n', GID.fontweight);
fprintf(fid,'st4h = uicontrol(fh,''style'',''text'', ...\n');
fprintf(fid,'      ''String'',''%s'', ...\n',GID.fn2);
fprintf(fid,'      ''Units'',''normalized'', ...\n');
fprintf(fid,'      ''Position'',[0.05,0.02,0.60,0.05], ...\n');
fprintf(fid,'      ''HorizontalAlignment'',''Left'', ...\n');
fprintf(fid,'      ''foregroundColor'',[1 0 0], ...\n');
fprintf(fid,'      ''FontSize'',%d, ...\n',GID.fontsize);
fprintf(fid,'      ''FontWeight'',''%s'');\n',GID.fontweight);
fprintf(fid,'\n');

% this is for either RT or FD (by virtue of GID value)advanced run

fprintf(fid,'runh = uicontrol(fh,''Style'',''pushbutton'', ...\n');
fprintf(fid,'      ''String'',''RUN'', ...\n');
fprintf(fid,'      ''Units'',''normalized'', ...\n');
fprintf(fid,'      ''Position'',[x0run,y0run,dxrun,dyrun], ...\n');
fprintf(fid,'      ''FontSize'',%d, ...\n',fontsize_run);
fprintf(fid,'      ''FontWeight'',''%s'', ...\n',fontweight_run);
fprintf(fid,'      ''callback'',{@run1_Callback ''%s''});\n',GID.run);
fprintf(fid,'PDFbuttons(pb1h);\n');
fprintf(fid,'\n\n');

% use the bgcolor of uicontrol for the GUI figure
fprintf(fid,'etbgcolor = get(et1h(1),''BackgroundColor'');\n');
fprintf(fid,'set(fh,''Color'',etbgcolor);\n');
% Now, turns on the GUI figure
fprintf(fid,'set(fh,''Visible'',''on'');\n');
fprintf(fid,'\n\n');

% Callback functions needed fot this GUI are generated below
% callback function of the "doc1" type
fprintf(fid, ...
    ['function opendoc1_Callback(hObject,eventdata,varargin)\n' ...
     '%% hObject    handle to Exit_menu_item (see GCBO)\n' ...
     '%% eventdata  reserved - to be defined in a future version of MATLAB\n' ...
     'open(varargin{1});\n' ...
     '\n\n']);

% edit text callback
fprintf(fid, ...
  ['function edit1_Callback(hObject, eventdata, varargin)\n' ...
   '%% hObject    handle to edit1 (see GCBO)\n' ...
   '%% eventdata  reserved - to be defined in a future version of MATLAB\n' ...
   '\n' ...
   '%% Hints: get(hObject,''String'') returns contents of edit1 as text\n' ...
   '%%        str2double(get(hObject,''String'')) returns contents of edit1 as a double\n' ...
   '\n' ...
   'fh = get(hObject,''parent'');\n' ...
   'data = getappdata(fh, ''myData'');\n' ...
   'value = str2num(get(hObject,''String''));\n' ...
   'eval([''data.'' varargin{1} ''= value;'']);\n' ...
   'setappdata(fh,''myData'',data);\n' ...
   '\n\n']);

% RUN callback
fprintf(fid, ...
  ['function run1_Callback(hObject, eventdata, varargin)\n' ...
   '%% hObject    handle to run1_Callback \n' ...
   '%% eventdata  reserved - to be defined in a future version of MATLAB\n' ...
   '%% Hints: get(hObject,''String'') returns contents of edit1 as text\n' ...
   '%%        str2double(get(hObject,''String'')) returns contents of edit1 as a double\n' ...
   'fh = get(hObject,''parent'');\n' ...
   'data = getappdata(fh, ''myData'');\n']);
   fprintf(fid,'for i=1:length(data.str)\n');
   fprintf(fid,'   eval([data.str{i} ''= data.'' data.str{i} '';'']);  %% e.g., A9 = data.A9;\n');
   fprintf(fid,'end\n');

fprintf(fid, ...
  ['%% Prevent user from changing the contents of edit boxes\n' ...
   'editBoxEnable2(data,''Off'');\n' ...
   'old_color=get(hObject,''BackgroundColor''); %% save original button color\n' ...
   'set(hObject, ''BackgroundColor'', [1 0 0]); %% Change button color\n' ...
   'set(hObject, ''String'', ''running ...'');\n' ...
   '%% open(''RTcaption.pdf'');\n' ...
   'drawnow    %% to make sure that the "Run" button change color right away\n' ...
   'eval(varargin{1})   %% run the m-file specified by varargin{1}\n']);
fprintf(fid, ...
  ['set(hObject, ''String'', ''Run'');\n' ...
   'set(hObject, ''BackgroundColor'', old_color);\n' ...
   'editBoxEnable2(data,''On'');\n']);
fprintf(fid,'\n\n');

fclose(fid);
disp([char(10) '***** The ' GID.guiName ' GUI is built successfully *****' char(10)]);


text:
function [F,L,Qc,H,Pinf,dF,dQc,dPinf,params] = cf_matern32_to_ss(magnSigma2, lengthScale)
% CF_MATERN32_TO_SS - Matern covariance functions, nu = 3/2, to state space
%
% Syntax:
%   [F,L,Qc,H,Pinf,dF,dQc,dPinf,params] = cf_matern32_to_ss(magnSigma2, lengthScale)
%
% In:
%   magnSigma2  - Matern magnitude scale parameter (default: 1)
%   lengthScale - Matern distance scale parameter (default: 1)
%
% Out:
%   F           - Feedback matrix
%   L           - Noise effect matrix
%   Qc          - Spectral density of white noise process w(t)
%   H           - Observation model matrix
%   Pinf        - Covariance of the stationary process
%   dF          - Derivatives of F w.r.t. parameters
%   dQc         - Derivatives of Qc w.r.t. parameters
%   dPinf       - Derivatives of Pinf w.r.t. parameters
%   params      - Input and output parameter information
%
% Description:
%   This function converts one-dimensional covariance functions of
%   the Matern class to state space models. The covariance function
%   parametrization is as follows
%
%     k(tau) = magnSigma2 (1+sqrt(3) |tau|/lengthScale) exp(-sqrt(3) |tau|/lengthScale),
%
%   where magnSigma2 is the magnitude scale parameter, lengthScale the
%   distance scale parameter, and tau the time difference between states,
%   tau = t-t'.
%     This function takes the covariance function parameters as inputs and
%   outputs the corresponding state space model matrices. The state space
%   model is given as follows in terms of a stochastic differential
%   equation
%
%      df(t)/dt = F f(t) + L w(t),
%
%   where w(t) is a white noise process with spectral denisty Qc. The
%   observation model for discrete observation y_k of f(t_k) at step k,
%   is as follows
%
%      y_k = H f(t_k) + r_k, r_k ~ N(0, R),
%
%   where r_k is the Gaussian measurement noise with covariance R.
%     Pinf is the stationary covariance, where the value of Pinf(i,j),
%   is defined as follows
%
%      Pinf(i,j) = E[(f_i(t)-E[f_i(t)])(f_j(t)-E[f_j(t)])],
%
%   where f_i(t) is component i of state vector f(t).
%     Derivatives: All have same form. For example, dF has the following
%   form:
%
%       dF(:,:,1) = dF/d(magnSigma2 = input parameter_1),
%       dF(:,:,i) = dF/d(input parameter_i).
%
% References:
%
%   [1] Simo Sarkka, Arno Solin, Jouni Hartikainen (2013).
%       Spatiotemporal learning via infinite-dimensional Bayesian
%       filtering and smoothing. IEEE Signal Processing Magazine,
%       30(4):51-61.
%
%   [2] Jouni Hartikainen and Simo Sarkka (2010). Kalman filtering and
%       smoothing solutions to temporal Gaussian process regression
%       models. Proceedings of IEEE International Workshop on Machine
%       Learning for Signal Processing (MLSP).
%
% See also:
%   COV_MATERN32, SPEC_MATERN32
%
% Copyright:
%   2012-2014   Arno Solin
%   2013-2014   Jukka Koskenranta
%
%  This software is distributed under the GNU General Public
%  License (version 3 or later); please refer to the file
%  License.txt, included with the software, for details.

%% Apply defaults

  % Check if magnSigma2 is given
  if nargin < 1 || isempty(magnSigma2), magnSigma2 = 1; end

  % Check if lengthScale is given
  if nargin < 2 || isempty(lengthScale), lengthScale = 1; end

%% Form state space model

  % Derived constants
  lambda = sqrt(3)/lengthScale;

  % Feedback matrix
  F = [0,          1;
       -lambda^2,  -2*lambda];

  % Noise effect matrix
  L = [0;   1];

  % Spectral density
  Qc = 12*sqrt(3)/lengthScale^3*magnSigma2;

  % Observation model
  H = [1,   0];

%% Stationary covariance

  % Calculate Pinf only if requested
  if nargout > 4,

    Pinf = [magnSigma2, 0;
            0,          3*magnSigma2/lengthScale^2];

  end

%% Calculate derivatives

  % Calculate derivatives only if requested
  if nargout > 5

    % Derivative of F w.r.t. parameter magnSigma2
    dFmagnSigma2    =  [0,  0;
                        0,  0];

    % Derivative of F w.r.t parameter lengthScale
    dFlengthScale   =  [0,                  0;
                        6/lengthScale^3,    2*sqrt(3)/lengthScale^2];
    % Derivative of Qc w.r.t. parameter magnSigma2
    dQcmagnSigma2   =   12*sqrt(3)/lengthScale^3;

    % Derivative of Qc w.r.t. parameter lengthScale
    dQclengthScale  =   -3*12*sqrt(3)/lengthScale^4*magnSigma2;

    % Derivative of Pinf w.r.t. parameter magnSigma2
    dPinfmagnSigma2 = [1,0;
                       0,3/lengthScale^2];

    % Derivative of Pinf w.r.t. parameter lengthScale
    dPinflengthScale = [0,0;
                        0,-6*magnSigma2/lengthScale^3];

    % Stack all derivatives
    dF = zeros(2,2,2);
    dQc = zeros(1,1,2);
    dPinf = zeros(2,2,2);

    dF(:,:,1) = dFmagnSigma2;
    dF(:,:,2) = dFlengthScale;
    dQc(:,:,1) = dQcmagnSigma2;
    dQc(:,:,2) = dQclengthScale;
    dPinf(:,:,1) = dPinfmagnSigma2;
    dPinf(:,:,2) = dPinflengthScale;

  end

%% Return parameter names

  % Only return if requested
  if nargout > 8

    % Stationarity
    p.stationary = true;

    % Input parameter information
    p.in{1}.name = 'magnSigma2';  p.in{1}.default = 1; p.in{1}.opt = true;
    p.in{2}.name = 'lengthScale'; p.in{2}.default = 1; p.in{2}.opt = true;

    % Return parameter setup
    params = p;

  end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function [ pivot, lu, info ] = r8mat_to_r8plu ( n, a )

%*****************************************************************************80
%
%% R8MAT_TO_R8PLU factors a general matrix.
%
%  Discussion:
%
%    This routine is a simplified version of the LINPACK routine DGEFA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, real A(N,N), the matrix to be factored.
%
%    Output, integer PIVOT(N), a vector of pivot indices.
%
%    Output, real LU(N,N), an upper triangular matrix U and
%    the multipliers L which were used to obtain it.  The factorization
%    can be written A = L * U, where L is a product of permutation and
%    unit lower triangular matrices and U is upper triangular.
%
%    Output, integer INFO, singularity flag.
%    0, no singularity detected.
%    nonzero, the factorization failed on the INFO-th step.
%
  lu(1:n,1:n) = a(1:n,1:n);

  info = 0;

  for k = 1 : n-1
%
%  Find L, the index of the pivot row.
%
    l = k;
    for i = k+1 : n
      if ( abs ( lu(l,k) ) < abs ( lu(i,k) ) )
        l = i;
      end
    end

    pivot(k) = l;
%
%  If the pivot index is zero, the algorithm has failed.
%
    if ( lu(l,k) == 0.0 )
      info = k;
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8MAT_TO_R8PLU - Fatal error!\n' );
      fprintf ( 1, '  Zero pivot on step %d\n', info );
      return
    end
%
%  Interchange rows L and K if necessary.
%
    if ( l ~= k )
      temp    = lu(l,k);
      lu(l,k) = lu(k,k);
      lu(k,k) = temp;
    end
%
%  Normalize the values that lie below the pivot entry A(K,K).
%
    lu(k+1:n,k) = -lu(k+1:n,k) / lu(k,k);
%
%  Row elimination with column indexing.
%
    for j = k+1 : n

      if ( l ~= k )
        temp    = lu(l,j);
        lu(l,j) = lu(k,j);
        lu(k,j) = temp;
      end

      lu(k+1:n,j) = lu(k+1:n,j) + lu(k+1:n,k) * lu(k,j);

    end

  end

  pivot(n) = n;

  if ( lu(n,n) == 0.0 )
    info = n;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_TO_R8PLU - Fatal error!\n' );
    fprintf ( 1, '  Zero pivot on step %d\n', info );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_rereference

% WALLTIME 00:10:00
% MEM 3gb
% DEPENDENCY ft_preprocessing ft_prepare_montage ft_apply_montage preproc

%% avg

data = [];
data.label = {'1', '2', '3'};
data.time{1} = (1:1000)/1000;
data.trial{1} = rand(length(data.label), length(data.time{1}));

cfg = [];
cfg.implicitref = [];
cfg.reref = 'yes';
cfg.refmethod = 'avg';
cfg.refchannel = 'all';
data1 = ft_preprocessing(cfg, data);

assert(~isequal(data.trial{1}, data1.trial{1}))
assert(length(data1.label)==3);
assert(all(abs(sum(data1.trial{1},1)) < 10*eps)); % the mean should be zero

cfg = [];
cfg.implicitref = '4';
cfg.reref = 'yes';
cfg.refmethod = 'avg';
cfg.refchannel = 'all';
data2 = ft_preprocessing(cfg, data);

assert(length(data2.label)==4);
assert(all(abs(sum(data2.trial{1},1)) < 10*eps)); % the mean should be zero

cfg = [];
cfg.implicitref = [];
cfg.reref = 'yes';
cfg.refmethod = 'avg';
cfg.refchannel = {'1', '2'};
data3 = ft_preprocessing(cfg, data);
assert(all(abs(sum(data3.trial{1}(1:2,:),1)) < 10*eps)); % the mean over the first 2 channels should be zero
assert(~all(abs(mean(data3.trial{1},1)) < 10*eps)); % the mean over all channels should NOT be zero

%% median

data = [];
data.label = {'1', '2', '3'};
data.time{1} = (1:1000)/1000;
data.trial{1} = rand(length(data.label), length(data.time{1}));

cfg = [];
cfg.implicitref = [];
cfg.reref = 'yes';
cfg.refmethod = 'median';
cfg.refchannel = 'all';
data1 = ft_preprocessing(cfg, data);

assert(~isequal(data.trial{1}, data1.trial{1}))
assert(length(data1.label)==3);
assert(all(abs(median(data1.trial{1},1)) < 10*eps)); % the median should be zero
assert(~all(abs(mean(data1.trial{1},1)) < 10*eps)); % the mean should NOT be zero

cfg = [];
cfg.implicitref = '4';
cfg.reref = 'yes';
cfg.refmethod = 'median';
cfg.refchannel = 'all';
data2 = ft_preprocessing(cfg, data);

assert(length(data2.label)==4);
assert(all(abs(median(data2.trial{1},1)) < 10*eps)); % the median should be zero
assert(~all(abs(mean(data2.trial{1},1)) < 10*eps)); % the mean should NOT be zero

cfg = [];
cfg.implicitref = [];
cfg.reref = 'yes';
cfg.refmethod = 'median';
cfg.refchannel = {'1', '2'};
data3 = ft_preprocessing(cfg, data);
assert(all(abs(median(data3.trial{1}(1:2,:),1)) < 10*eps)); % the median over the first 2 channels should be zero
assert(~all(abs(median(data3.trial{1},1)) < 10*eps)); % the median over all channels should NOT be zero
assert(~all(abs(mean(data3.trial{1},1)) < 10*eps)); % the mean over all channels should NOT be zero

%% rest

elec = [];
elec.label = arrayfun(@num2str, 1:162, 'UniformOutput', false)';
elec.elecpos = mesh_sphere(162);

headmodel = [];
headmodel.type = 'singlesphere';
headmodel.r = 1;
headmodel.o = [0 0 0];
headmodel.cond = 1;

cfg = [];
cfg.headmodel = headmodel;
cfg.elec = elec;
cfg.method = 'basedonvol';
cfg.inwardshift = 0.1;
sourcemodel = ft_prepare_sourcemodel(cfg);

if false
  figure
  ft_plot_headmodel(headmodel);
  alpha 0.5
  ft_plot_mesh(sourcemodel);
  ft_plot_sens(elec, 'label', 'label', 'elecshape', 'disc');
end

cfg = [];
cfg.headmodel = headmodel;
cfg.elec = elec;
cfg.sourcemodel = sourcemodel;
leadfield = ft_prepare_leadfield(cfg);

lfmat = cat(2, leadfield.leadfield{:});
assert(isequal(size(lfmat), [162, 642*3]));
assert(~any(isnan(lfmat(:))));
assert(all(abs(mean(lfmat,1)) < 10*eps)); % the leadfield should be average referenced

data = [];
data.label = arrayfun(@num2str, 1:162, 'UniformOutput', false)';
data.time{1} = (1:1000)/1000;
data.trial{1} = rand(length(data.label), length(data.time{1}));

cfg = [];
cfg.implicitref = [];
cfg.reref = 'yes';
cfg.refmethod = 'rest';
cfg.refchannel = 'all';
cfg.leadfield = leadfield;
data1 = ft_preprocessing(cfg, data);

% the REST rereferenced data should have a mean that is close to zero
% since the electrode and dipole distribution are very homogenous

assert(all(abs(mean(data1.trial{1},1)) < 10*eps));

%% bipolar

nelec = 10; % per shaft
nshaft = 5;

shaft = {'A', 'B', 'C', 'D', 'E'};

data = [];
data.label = {};
for i=1:nshaft
  for j=1:nelec
    data.label{end+1} = sprintf('%s%d', shaft{i}, j);
  end
end
data.time{1} = (1:1000)/1000;
data.trial{1} = rand(length(data.label), length(data.time{1}));

% note that the ordering should be 1, 2, 3, ... 10, and not 1, 10, 2, 3, 4, ...

cfg = [];
cfg.reref = 'yes';
cfg.refmethod = 'bipolar';
cfg.groupchans = 'no';
data1 = ft_preprocessing(cfg, data);
assert(length(data1.label)==49); % 50 channels, hence 49 differences

a1_min_a2 = data.trial{1}(1,:) - data.trial{1}(2,:);
assert(all(abs(data1.trial{1}(1,:)-a1_min_a2)<10*eps));

cfg = [];
cfg.reref = 'yes';
cfg.refmethod = 'bipolar';
cfg.groupchans = 'yes';
data2 = ft_preprocessing(cfg, data);
assert(length(data2.label)==45); % 5x10 channels, hence 5x9 differences

a1_min_a2 = data.trial{1}(1,:) - data.trial{1}(2,:);
assert(all(abs(data2.trial{1}(1,:)-a1_min_a2)<10*eps));

%% laplace

nelec = 10; % per shaft
nshaft = 5;

shaft = {'A', 'B', 'C', 'D', 'E'};

data = [];
data.label = {};
for i=1:nshaft
  for j=1:nelec
    data.label{end+1} = sprintf('%s%d', shaft{i}, j);
  end
end
data.time{1} = (1:1000)/1000;
data.trial{1} = rand(length(data.label), length(data.time{1}));

% note that the ordering should be 1, 2, 3, ... 10, and not 1, 10, 2, 3, 4, ...

cfg = [];
cfg.reref = 'yes';
cfg.refmethod = 'laplace';
cfg.groupchans = 'no';
data1 = ft_preprocessing(cfg, data);

assert(length(data1.label)==50);

cfg = [];
cfg.reref = 'yes';
cfg.refmethod = 'laplace';
cfg.groupchans = 'yes';
data2 = ft_preprocessing(cfg, data);

assert(length(data1.label)==50);

% channel 1 is in both cases on the edge, hence the same
% channel 10 and 11 are in one case in the middle, in the other at two separate edges of the shafts
% channel 2-9 are always embedded in the shaft, so the same

assert( isequal(data1.trial{1}(1,:), data2.trial{1}(1,:)));
assert(~isequal(data1.trial{1}(10,:), data2.trial{1}(10,:)));
assert(~isequal(data1.trial{1}(11,:), data2.trial{1}(11,:)));
assert( isequal(data1.trial{1}(2:9,:), data2.trial{1}(2:9,:)));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Data = loadsamples(imgpathlistfile, exc_setlabel)
%LOADSAMPLES Summary of this function goes here
%   Function: load samples from dbname database
%   Detailed explanation goes here
%   Input:
%        dbname: the name of one database
%        exc_setlabel: excluded set label
%   Output:
%        Data: loaded data from the database
imgpathlist = textread(imgpathlistfile, '%s', 'delimiter', '\n');

Data = cell(length(imgpathlist), 1);

setnames = {'train' 'test'};

% Create a cascade detector object.
% faceDetector = vision.CascadeObjectDetector();
% bboxes_facedet = zeros(length(imgpathlist), 4);
% bboxes_gt      = zeros(length(imgpathlist), 4);
% isdetected = zeros(length(imgpathlist), 1);

parfor i = 1:length(imgpathlist)
    img = im2uint8(imread(imgpathlist{i}));
    Data{i}.width_orig    = size(img, 2);
    Data{i}.height_orig   = size(img, 1);

    % Data{i}.img      = img
    % shapepath = strrep(imgpathlist{i}, 'png', 'pts');
    shapepath = strcat(imgpathlist{i}(1:end-3), 'pts');
    Data{i}.shape_gt = double(loadshape(shapepath));
    % Data{i}.shape_gt = Data{i}.shape_gt(params.ind_usedpts, :);
    % bbox     =  bounding_boxes_allsamples{i}.bb_detector; %
    Data{i}.bbox_gt = getbbox(Data{i}.shape_gt); % [bbox(1) bbox(2) bbox(3)-bbox(1) bbox(4)-bbox(2)];

    % cut original image to a region which is a bit larger than the face
    % bounding box
    region = enlargingbbox(Data{i}.bbox_gt, 2.0);

    region(2) = double(max(region(2), 1));
    region(1) = double(max(region(1), 1));

    bottom_y = double(min(region(2) + region(4) - 1, Data{i}.height_orig));
    right_x = double(min(region(1) + region(3) - 1, Data{i}.width_orig));

    img_region = img(region(2):bottom_y, region(1):right_x, :);
    Data{i}.shape_gt = bsxfun(@minus, Data{i}.shape_gt, double([region(1) region(2)]));

    % to save memory cost during training
    if exc_setlabel == 2
        ratio = min(1, sqrt(single(150 * 150) / single(size(img_region, 1) * size(img_region, 2))));
        img_region = imresize(img_region, ratio);
        Data{i}.shape_gt = Data{i}.shape_gt .* ratio;
    end
    Data{i}.bbox_gt = getbbox(Data{i}.shape_gt);

    Data{i}.bbox_facedet = getbbox(Data{i}.shape_gt);
    % perform face detection using matlab face detector
    %{
    bbox = step(faceDetector, img_region);
    if isempty(bbox)
        % if face detection is failed
        isdetected(i) = 1;
        Data{i}.bbox_facedet = getbbox(Data{i}.shape_gt);
    else
        int_ratios = zeros(1, size(bbox, 1));
        for b = 1:size(bbox, 1)
            area = rectint(Data{i}.bbox_gt, bbox(b, :));
            int_ratios(b) = (area)/(bbox(b, 3)*bbox(b, 4) + Data{i}.bbox_gt(3)*Data{i}.bbox_gt(4) - area);
        end
        [max_ratio, max_ind] = max(int_ratios);

        if max_ratio < 0.4  % detection fail
            isdetected(i) = 0;
        else
            Data{i}.bbox_facedet = bbox(max_ind, 1:4);
            isdetected(i) = 1;
            % imgOut = insertObjectAnnotation(img_region,'rectangle',Data{i}.bbox_facedet,'Face');
            % imshow(imgOut);
        end
    end
    %}
    % recalculate the location of groundtruth shape and bounding box
    % Data{i}.shape_gt = bsxfun(@minus, Data{i}.shape_gt, double([region(1) region(2)]));
    % Data{i}.bbox_gt = getbbox(Data{i}.shape_gt);

    if size(img_region, 3) == 1
        Data{i}.img_gray = img_region;
    else
        % hsv = rgb2hsv(img_region);
        Data{i}.img_gray = rgb2gray(img_region);
    end

    Data{i}.width    = size(img_region, 2);
    Data{i}.height   = size(img_region, 1);
end

ind_valid = ones(1, length(imgpathlist));
parfor i = 1:length(imgpathlist)
    if ~isempty(exc_setlabel)
        ind = strfind(imgpathlist{i}, setnames{exc_setlabel});
        if ~isempty(ind) % | ~isdetected(i)
            ind_valid(i) = 0;
        end
    end
end

% learn the linear transformation from detected bboxes to groundtruth bboxes
% bboxes = [bboxes_gt bboxes_facedet];
% bboxes = bboxes(ind_valid == 1, :);

Data = Data(ind_valid == 1);

end

function shape = loadshape(path)
% function: load shape from pts file
file = fopen(path);

if ~isempty(strfind(path, 'COFW'))
    shape = textscan(file, '%d16 %d16 %d8', 'HeaderLines', 3, 'CollectOutput', 3);
else
    shape = textscan(file, '%d16 %d16', 'HeaderLines', 3, 'CollectOutput', 2);
end
fclose(file);

shape = shape{1};
end

function region = enlargingbbox(bbox, scale)

region(1) = floor(bbox(1) - (scale - 1)/2*bbox(3));
region(2) = floor(bbox(2) - (scale - 1)/2*bbox(4));

region(3) = floor(scale*bbox(3));
region(4) = floor(scale*bbox(4));

% region.right_x = floor(region.left_x + region.width - 1);
% region.bottom_y = floor(region.top_y + region.height - 1);

end

author: jwyang

repo: face-alignment


text:
function [max_min_diff, total_max_min_diff] = minimax(M)
    max_rows = max(M');
    min_rows = min(M');
    max_min_diff = max_rows - min_rows;
    total_max_min_diff = max(max_rows) - min(min_rows);
end

author: anishLearnsToCode

repo: introduction-to-programming-with-matlab


text:
function opts = getGlobalDefaults()
% Returns a structure with default option values for Manopt.
%
% function opts = getGlobalDefaults()
%
% Returns a structure opts containing the global default options such as
% verbosity level etc. Typically, global defaults are overwritten by solver
% defaults, which are in turn overwritten by user-specified options.
% See the online Manopt documentation for details on options.
%
% See also: mergeOptions

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:

    % There should be no reason to modify this file.
    % For better compatibility with future Manopt versions,
    % use the options structure of solvers.
    %
    % Really: don't modify it.

    % Verbosity level: 0 is no output at all. The higher the verbosity, the
    % more info is printed / displayed during solver execution.
    opts.verbosity = 3;

    % If debug is set to true, additional computations may be performed and
    % debugging information is outputed during solver execution.
    opts.debug = false;

    % Maximum number of store structures to store. If set to 0, caching
    % capabilities are not disabled, but the cache will be emptied at each
    % iteration of iterative solvers (more specifically: every time the
    % solver calls to purge the storedb).
    opts.storedepth = 20;

    % Maximum amount of time a solver may execute, in seconds.
    opts.maxtime = inf;

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function [M,G]=triSurf2ImSpec(F,V,voxelSize)

% function [M,G,bwLabels]=triSurf2Im(F,V,voxelSize)
% -----------------------------------------------------------------------
% This function converts the input triangulated surface, specified by the
% faces F and the vertices V into an image based on the voxel size
% specified in the input voxelSize. If the latter is empty then the
% voxelSize is based on the maximum edge length. If required the surface is
% resampled to subvoxel resolution using the |subtri| function. Then
% surface vertices are simply mapped to an image coordinate system and
% since they are densely sampled with respect to the voxel size form an
% enclosing boundary of voxels. The exterior, boundary and interior are
% then formulated using the bwlabeln function. The exterior, boundary and
% intertior voxels are labelled using 0's, 1's and 2's in the output image
% M. The second output G is a structure containing the fields G.voxelSize
% and G.origin which form the geometric information for the image.
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 28/08/2013
%------------------------------------------------------------------------

%%
%Checking edge lenghts of surface
TR = triangulation(F,V);
E=edges(TR);
Ed=0;
for q=1:1:size(V,2)
    Ed=Ed+(V(E(:,1),q)-V(E(:,2),q)).^2;
end
edgeLengths=sqrt(Ed);
maxEdgeLength=max(edgeLengths(:));

if isempty(voxelSize)
    voxelSize=maxEdgeLength;
end

%Resample surface if voxelsize is small with respect to edgelenghts
n=maxEdgeLength/voxelSize;
if n>(1-eps(1))
    n=floor(n);
    [~,V]=subtri(F,V,n);
end

%Determine surface set coordinate minima
minV=min(V,[],1);

%Determine shift so all coordinates are positive
imOrigin=-(minV-voxelSize);

%Shift points using origin
V=V+imOrigin(ones(size(V,1),1),:);

%Convert to image coordinates
V_IJK=V;
[V_IJK(:,1),V_IJK(:,2),V_IJK(:,3)]=cart2im(V(:,1),V(:,2),V(:,3),voxelSize*ones(1,3));

%Rounding image coordinates to snap to voxel
V_IJK=round(V_IJK);

%Determin image size
siz=max(V_IJK,[],1)+1;

%Get linear indices of points
indV=sub2ind(siz,V_IJK(:,1),V_IJK(:,2),V_IJK(:,3));

%Create surface boundary image
L=false(siz);
L(indV)=1;

%Create boundary, interior and exterior image
LL = bwlabeln(~L,6); %Get labels for ~L image which will segment interior, exterior and boundary
uniqueLabels=unique(LL(:));

labelsBoundary=LL(L); %The label numbers for the boundary

indExteriorVoxel=1; %First is outside since image is at least a voxel too big on all sides
labelExterior=LL(indExteriorVoxel); %The label number for the exterior

labelsInterior=uniqueLabels(~ismember(uniqueLabels,[labelsBoundary(:); labelExterior])); %Labels for the interior (possibly multiple)

M=zeros(size(L)); %The exterior is set to 0
M(L)=1; %The boundary is set to 1
M(ismember(LL,labelsInterior))=2; %Interior is set to 2

%Storing image geometry metrics
G.voxelSize=voxelSize;
G.origin=imOrigin;

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function x = r8row_to_r8vec ( m, n, a )

%*****************************************************************************80
%
%% R8ROW_TO_R8VEC converts an R8ROW into an R8VEC.
%
%  Example:
%
%    M = 3, N = 4
%
%    A =
%      11 12 13 14
%      21 22 23 24
%      31 32 33 34
%
%    X = ( 11, 12, 13, 14, 21, 22, 23, 24, 31, 32, 33, 34 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), the R8ROW.
%
%    Output, real X(M*N), a vector containing the M rows of A.
%
  j = 1;
  for i = 1 : m
    x(j:j+n-1) = a(i,1:n);
    j = j + n;
  end

  return
end


text:
function varargout = sure2(varargin)
% SURE2 M-file for sure2.fig
%      SURE2, by itself, creates a new SURE2 or raises the existing
%      singleton*.
%
%      H = SURE2 returns the handle to a new SURE2 or the handle to
%      the existing singleton*.
%
%      SURE2('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in SURE2.M with the given input arguments.
%
%      SURE2('Property','Value',...) creates a new SURE2 or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before sure2_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to sure2_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 sure2

% Last Modified by GUIDE v2.5 21-Jan-2008 02:39:36

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @sure2_OpeningFcn, ...
                   'gui_OutputFcn',  @sure2_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 sure2 is made visible.
function sure2_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 sure2 (see VARARGIN)

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

% Update handles structure
guidata(hObject, handles);

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

% --- Outputs from this function are returned to the command line.
function varargout = sure2_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;

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

close('VScope');
%close('VAudio')
close('sure2');
%close('VAudio');
close('sure');

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

close('sure2');

author: Sable

repo: mcbench-benchmarks


text:
%% load
load cnt_ssvep;
CNT = cnt_ssvep;

% cnt variables
% cnt.t  : time information
% cnt.fs : sampling frequency
% cnt.y_dec : class information (e.g., up = 1, left = 2, right = 3, down = 4)
% cnt.y_logic : logical format of class inforamtion
% cnt_y_class : class name (e.g., up, left, right, and down)
% cnt.class : number of class
% cnt.chan : number of electrodes
% cnt. x : raw eeg signals

% revised 2017.11.11 - Oyeon Kwon (oy_kwon@korea.ac.kr)
%% Initialization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
time = 4;
freq = [5, 7, 9, 11];
fs = CNT.fs;
interval = [0 4000];  %
marker  = {'1','up';'2', 'left';'3', 'right';'4', 'down'};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CCA - Analysis
SMT = prep_segmentation(CNT, {'interval', interval});

numTrials = size(SMT.x, 2);
count= numTrials;

for i = 1: size(SMT.x, 2)
    cca_result = ssvep_cca_analysis(squeeze(SMT.x(:,i,:)),{'marker',marker;'freq', freq;'fs', fs;'time',time});
    [~, ind] = max(cca_result);
    if SMT.y_dec(i) ~= ind
        count = count -1;
    end
end
fprintf('Accuracy: %.2f%%\n', count/numTrials);
Accuracy = count/numTrials;

author: PatternRecognition

repo: OpenBMI


text:
function [Jul1,Jul2]=UTC2TT(Jul1,Jul2)
%%UTC2TT Convert from universal coordinated time (UTC) given as a two-part
%        pseudo-Julian date to terrestrial time (TT), represented as a two-
%        part Julian date.
%
%INPUTS: Jul1, Jul2 Matrices of two parts of a pseudo-Julian date given in
%                   UTC. The units of the date are days. The full date is
%                   the sum of both terms. The date is broken into two
%                   parts to provide more bits of precision. It does not
%                   matter how the date is split. Corresponding elements in
%                   each matrix are times that are converted.
%
%OUTPUTS: Jul1, Jul2 The time as a Julian date in TT with the same
%                    dimensionalities as the input sets of dates.
%
%The UTC date is only pseudo-Julian, because there is not a fixed number
%of seconds in a Julian day. The convention used in the IAU standard is
%that the Julian day matches the UTC day regardless of whether the UTC day
%is 86399, 86400 or 86401 SI seconds (depending on the presence of leap
%seconds).
%
%UTC began at 1960 January 1.0 (JD 2436934.5) and this function should not
%be called with an earlier date.
%
%This just calls a number of intermediate conversion functions out of the
%International Astronomical Union's (IAU) Standard's of Fundamental
%Astronomy library.
%
%Many temporal coordinate systems standards are compared in [1].
%
%REFERENCES:
%[1] D. F. Crouse, "An Overview of Major Terrestrial, Celestial, and
%    Temporal Coordinate Systems for Target Tracking," Formal Report, Naval
%    Research Laboratory, no. NRL/FR/5344--16-10,279, 10 Aug. 2016, 173
%    pages.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

[Jul1,Jul2]=UTC2TAI(Jul1,Jul2);
[Jul1,Jul2]=TAI2TT(Jul1,Jul2);
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:
function r = ldivide(a,b)
%LDIVIDE      Gradient elementwise left division a .\ b
%

% written  11/02/05     S.M. Rump
%

  r = b ./ a;

repo: OpenMAS


text:
function varargout = triangleGrid(bounds, origin, size, varargin)
%TRIANGLEGRID Generate triangular grid of points in the plane.
%
%   usage
%   PTS = triangleGrid(BOUNDS, ORIGIN, SIZE)
%   generate points, lying in the window defined by BOUNDS, given in form
%   [xmin ymin xmax ymax], starting from origin with a constant step equal
%   to size.
%   SIZE is constant and is equals to the length of the sides of each
%   triangles.
%
%   TODO: add possibility to use rotated grid
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 06/08/2005.
%

dx = size(1);
dy = size(1)*sqrt(3);

% consider two square grids with different centers
pts1 = squareGrid(bounds, origin, [dx dy], varargin{:});
pts2 = squareGrid(bounds, origin + [dx dy]/2, [dx dy], varargin{:});

% gather points
pts = [pts1;pts2];

% process output
if nargout>0
    varargout{1} = pts;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
% Convert seconds to milliseconds
function ms = sec2ms(sec)
    ms = sec * 1000;
end

author: kavli-ntnu


text:
function mat_gram = GetKernelMat(data1,data2,kernel_type,kernel_param)
% Ccomputes the Gram matrix of a specified kernel function.
%Input: data1, data2-Input data matix
%            kernel_type: 'linear' | 'poly' | 'rbf'|'intersect'|'sigmoid'
% 2016-10-16 by jlfeng

n1=size(data1,1);
n2=size(data2,1);

switch kernel_type
    case 'linear'
        mat_gram=data2*data1';
    case 'poly'
        mat_gram=(data2*data1'+1).^kernel_param;
    case 'polyhomog'
        mat_gram=(data2*data1').^kernel_param;
    case 'rbf'
        mat_gram = exp(-(repmat(sum(data1.*data1,2)',n2,1) + ...
                repmat(sum(data2.*data2,2),1,n1) - 2*data2*data1') ...
                /(2*kernel_param^2));
    case 'intersect'
        mat_gram=zeros(n2,n1);
        for kk=1:n1
            temp=min(data2,repmat(data1(kk,:),[n2 1]));
            mat_gram(:,kk)=sum(temp,2);
        end
    case 'sigmod'
        mat_gram=tanh(data2*data1'*kernel_param(1)+kernel_param(2));
    otherwise
        error('Unknown kernel function.');
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [u,v,I]=localcross(Y,tol)
% Full-pivoted cross for truncating one ket TT block instead of SVD

if (exist('localcross_mex','file')>0)&&(isreal(Y))
    % Check for a faster mex
    if (nargout>2)
        [u,v,I]=localcross_mex(Y,tol);
    else
        [u,v]=localcross_mex(Y,tol);
    end;
    return;
end;

[n,m,b]=size(Y);

minsz = min(n, m*b);
u = zeros(n, minsz);
v = zeros(minsz, m*b);
res = reshape(Y, n, m*b);
% Return also the indices
I = zeros(1, minsz);
val_max = max(abs(Y(:)));
for r=1:minsz
    res = reshape(res, [], 1);
    [val,piv]=max(abs(res));
    piv = tt_ind2sub([n, m*b], piv);
    if (val<=tol*val_max)
        break;
    end;
    res = reshape(res, n, m*b);
    u(:,r) = res(:, piv(2));
    v(r,:) = res(piv(1), :)/res(piv(1),piv(2));
    res = res - u(:,r)*v(r,:);
    I(r)=piv(1);
end;
% Return indices
r = find(I==0,1);
if (isempty(r))
    r=minsz;
else
    r=r-1;
end;
I = I(1:r);
u = u(:, 1:r);
v = v(1:r, :);
if (r==0)
    u = zeros(n,1);
    v = zeros(1,m);
    I = 1;
end;
% qr u, in case we don't have enrichment
[u,rv]=qr(u,0);
v = rv*v;
end

author: oseledets

repo: TT-Toolbox


text:
%% GSA_FAST_GetSi: calculate the FAST sensitivity indices
% Ref: Cukier, R.I., C.M. Fortuin, K.E. Shuler, A.G. Petschek and J.H.
% Schaibly (1973). Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients. I Theory. Journal of Chemical Physics
%
% Max number of input variables: 50
%
% Usage:
%   Si = GSA_FAST_GetSi(pro)
%
% Inputs:
%    pro                project structure
%
% Output:
%    Si                 vecotr of sensitivity coefficients
%
% ------------------------------------------------------------------------
% See also
%
% Author : Flavio Cannavo'
% e-mail: flavio(dot)cannavo(at)gmail(dot)com
% Release: 1.0
% Date   : 01-05-2011
%
% History:
% 1.0  01-05-2011  First release.
%%

function Si = GSA_FAST_GetSi(pro)

k = length(pro.Inputs.pdfs);

M = 4;

W = fnc_FAST_getFreqs(k);

Wmax = W(k);
N = 2*M*Wmax+1;
q = (N-1)/2;

S = pi/2*(2*(1:N)-N-1)/N;
alpha = W'*S;

NormedX = 0.5 + asin(sin(alpha'))/pi;

X = fnc_FAST_getInputs(pro, NormedX);

Y = nan(N,1);

for j=1:N
    Y(j) = pro.Model.handle(X(j,:));
end

A = zeros(N,1);
B = zeros(N,1);
N0 = q+1;

for j=2:2:N
    A(j) = 1/N*(Y(N0)+(Y(N0+(1:q))+Y(N0-(1:q)))'* ...
        cos(pi*j*(1:q)/N)');
end

for j=1:2:N
    B(j) = 1/N*(Y(N0+(1:q))-Y(N0-(1:q)))'* ...
        sin(pi*j*(1:q)/N)';
end

V = 2*(A'*A+B'*B);

for i=1:k
    Vi=0;
    for j=1:M
        Vi = Vi+A(j*W(i))^2+B(j*W(i))^2;
    end
    Vi = 2*Vi;
    Si(i) = Vi/V;
end


text:
function showmeshdensity(node,elem)
%% SHOWMESHDENSITY show the density of the mesh
%
% showmeshdensity(node,elem) display the size of each triangle using different
% color.
%
% Example
%   load lakemesh
%   showmeshdensity(node,elem)
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

ve2 = node(elem(:,1),:)-node(elem(:,3),:);
ve3 = node(elem(:,2),:)-node(elem(:,1),:);
area = 0.5*abs(-ve3(:,1).*ve2(:,2) + ve3(:,2).*ve2(:,1));
h = sqrt(area);
showsolution(node,elem,h)
colorbar
axis equal; axis tight; axis off;

repo: ifem


text:
function a = sym2polys(p,x)

%Sym2Polys  Extract the coefficients of a symbolic polynomial.
%           This function is an extension of the Matlab SYM2POLY and COEFFS
%           functions in that it allows the coefficients to be symbolic and
%           returns the full coefficient vector including the zero coefficients.
%
%Usage: c = sym2polys(p,x)
%       where p is the (multi) symbolic polynomial and x is the
%       independent variable. If x is not specified then the variable
%       alphabetically closest to x is used as the independent variable.
%
%Example:    If p = a*b*x^3 + b*c*x + c*d
%            then sym2polys(p) returns [a*b, 0, b*c, c*d]
%            whereas sym2polys(p,'b') returns [a*x^3+c*x, c*d]
%            Note that coeffs(p,x) returns [c*d, b*c, a*b]

%see also: sym2poly, coeffs

% Mukhtar Ullah
% mukhtar.ullah@informatik.uni-rostock.de
% September 2, 2004

if nargin == 1, x = findsym(p,1); end

[c,t] = coeffs(p,x);
i = sym2poly(sort(sum(t*(1:numel(t)).')));
a = sym(i);
a(i>0) = c(i(i>0));

if isempty(findsym(a)), a = double(a); end


text:
function varargout=bss_decomp_filt(varargin)

% decompose an estimated source into target/interference/noise/artefacts components, assuming the admissible distortion is a pure time-invariant filter.
% components, assuming the admissible distortion is a pure time-invariant
% filter.
%
% Usage:
%
% [s_target,e_interf[,e_noise],e_artif]=bss_decomp_filt(se,index,S[,N],L)
%
% Input:
%   - se: row vector of length T containing the estimated source,
%   - index: points which component of S se has to be compared to,
%   - S: n x T matrix containing the original sources,
%   - N: m x T matrix containing the noise on the observations (if any).
%   - L: the number of lags of the allowed filter
%
% Output:
%   - s_target: row vector of length T containing the target source(s)
%   contribution,
%   - e_interf: row vector of length T containing the interferences
%   contribution,
%   - e_noise: row vector of length T containing the noise contribution (if
%   any),
%   - e_artif: row vector of length T containing the artifacts
%   contribution.
%
% Developers:  - Cedric Fevotte (fevotte@tsi.enst.fr) - Emmanuel Vincent
% (emmanuel.vincent@irisa.fr) - Remi Gribonval (remi.gribonval@irisa.fr)

se=varargin{1}; index=varargin{2}; S=varargin{3};

switch nargin
    case 4
    N=[];
    L=varargin{4};
    case 5
    N=varargin{4};
    L=varargin{5};
    otherwise
    disp('Wrong number of arguments.')
end

[ne,Te]=size(se);
[n,T]=size(S);

%%%%%%%%%% WARNINGS %%%%%%%%%%%%%
switch isempty(N)
    case 1
        if n>T | ne>Te, disp('Watch out: signals must be in rows.'), return; end
        if ne~=1, disp('Watch out: se must contain only one row.'), return; end
        if T~=Te, disp('Watch out: se and S have different lengths.'), return; end
    case 0
        [m,Tm]=size(N);
        if n>T | ne>Te | m>Tm, disp('Watch out: signals must be in rows.'), return; end
        if ne~=1, disp('Watch out: se must contain only one row.'), return; end
        if T~=Te, disp('Watch out: S and Se have different lengths.'), return; end
        if T~=Tm, disp('Watch out: N, S and Se have different lengths.'), return; end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Create the space of target source(s)
target_space=bss_make_lags(S(index,:),L);
% Create the space of sources
sources_space=bss_make_lags(S,L);
% Create the noise space
noise_space=bss_make_lags(N,L);

s_target=zeros(1,T);
e_interf=zeros(1,T);
e_artif=zeros(1,T);
if isempty(noise_space)==0, e_noise=zeros(1,T); end

%%% Target source(s) contribution %%%
s_target = bss_proj(se,target_space);

%%% Interferences contribution %%%
P_S_se = bss_proj(se,sources_space);
e_interf = P_S_se - s_target;

switch isempty(noise_space)
    case 1 % No noise
        %%% Artifacts contribution %%%
        e_artif= se - P_S_se;

        %%% Output %%%
        varargout{1}=s_target;
        varargout{2}=e_interf;
        varargout{3}=e_artif;

    case 0 % Noise
        %%% Noise contribution %%%
        P_SN_se= bss_proj(se,[sources_space;noise_space]);
        e_noise=P_SN_se-P_S_se;

        %%% Artifacts contribution %%%
        e_artif=se-P_SN_se;

        %%% Output %%%
        varargout{1}=s_target;
        varargout{2}=e_interf;
        varargout{3}=e_noise;
        varargout{4}=e_artif;
end

author: posenhuang

repo: deeplearningsourceseparation


text:
function [TrainData, TestData, TrainLabel, TestLabel] = createData_TriTL()
    % source domain (train)
    A = load(strcat('data/TrainSelect_8.data'));
    B = spconvert(A);
    n1 = size(B,2);
    for i = 1:n1
        TrainY1(1,i) = 1;
    end
    A1 = load(strcat('data/TrainSelect_12.data'));
    B1 = spconvert(A1);
    n2 = size(B1,2);
    for i = 1:n2
        TrainY1(1,i+n1) = 2;
    end
    TrainX1 = [B,B1];
    clear A B A1 B1;

    A = load(strcat('data/TrainSelect_9.data'));
    B = spconvert(A);
    n1 = size(B,2);
    for i = 1:n1
        TrainY2(1,i) = 1;
    end
    A1 = load(strcat('data/TrainSelect_13.data'));
    B1 = spconvert(A1);
    n2 = size(B1,2);
    for i = 1:n2
        TrainY2(1,i+n1) = 2;
    end
    TrainX2 = [B,B1];
    clear A B A1 B1;

    A = load(strcat('data/TrainSelect_10.data'));
    B = spconvert(A);
    n1 = size(B,2);
    for i = 1:n1
        TrainY3(1,i) = 1;
    end
    A1 = load(strcat('data/TrainSelect_14.data'));
    B1 = spconvert(A1);
    n2 = size(B1,2);
    for i = 1:n2
        TrainY3(1,i+n1) = 2;
    end
    TrainX3 = [B,B1];
    clear A B A1 B1;

    % target domain (test)
    A = load(strcat('data/TrainSelect_11.data'));
    B = spconvert(A);
    n1 = size(B,2);
    for i = 1:n1
        TestY(1,i) = 1;
    end
    A1 = load(strcat('data/TrainSelect_15.data'));
    B1 = spconvert(A1);
    n2 = size(B1,2);
    for i = 1:n2
        TestY(1,i+n1) = 2;
    end
    TestX = [B,B1];
    clear A B A1 B1;

    TrainData{1,1} = TrainX1;
    TrainLabel{1,1} = TrainY1;
    TrainData{1,2} = TrainX2;
    TrainLabel{1,2} = TrainY2;
    TrainData{1,3} = TrainX3;
    TrainLabel{1,3} = TrainY3;
    TestData{1,1} = TestX;
    TestLabel{1,1} = TestY;
    save inputData.mat TrainData TrainLabel TestData TestLabel
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
classdef MrCompare < MrCopyData
    % This class loads and compares different objects, e.g. images/rois
    %
    % EXAMPLE
    %   MrCompare(imageArray, extraDimInfo)
    %
    %   See also MrDimInfo MrImage

    % Author:   Saskia Klein & Lars Kasper
    % Created:  2019-02-25
    % Copyright (C) 2019 Institute for Biomedical Engineering
    %                    University of Zurich and ETH Zurich
    %
    % This file is part of the TAPAS UniQC Toolbox, which is released
    % under the terms of the GNU General Public Licence (GPL), version 3.
    % You can redistribute it and/or modify it under the terms of the GPL
    % (either version 3 or, at your option, any later version).
    % For further details, see the file COPYING or
    %  <http://www.gnu.org/licenses/>.
    %

    properties

        % cell(nRows,1) of strings with detailed image information, e.g.
        % previous processing steps to arrive at that image
        % (masking/thresholding...)
        % 'detailed image information'; 'given in cell strings'
        info    = {};

        % A short string identifier of the image, used e.g. as plot-title
        name    = 'MrCompare';

        % cell(nExtraDim1, ..., nExtraDimN) holding handles to all objects
        % that shall be compared in different aspects
        % The dimensions are described by the dimInfo
        data = {};

        % MrDimInfo for data, describing which items are held in data
        dimInfo

    end % properties

    methods

        function this = MrCompare(objectArray, extraDimInfo)
            % Constructor of class
            % IN
            %   objectArray cell(nExtraDim1, ..., nExtraDimN)
            %                   OR
            %               cell(nObjects,1)
            %               holding handles to all objects of the same class
            %               (MrImage, MrRoi, MrSeries, or MrDataNd,
            %               MrImageSpm4D) to be included in comparison
            %   extraDimInfo
            %               MrDimInfo, describing meta-data of objectArray,
            %               e.g., objects from different subjects, sessions,
            %               protocols

            if nargin >= 1 % needed to allow empty constructor for copyobj etc
                % select only non-singleton-dims for dimInfo
                nSamples = size(objectArray);
                iDims = nSamples > 1;
                nDims = max(sum(iDims),1);
                nSamples = nSamples(iDims);
                if isempty(nSamples), nSamples = 1; end;

                if nargin >= 2
                    this.dimInfo = extraDimInfo.copyobj();
                else

                    for iDim = 1:nDims
                        dimLabels{iDim} = sprintf('dim%d', iDim);
                        units{iDim}      = '';
                    end
                    this.dimInfo = MrDimInfo('dimLabels', dimLabels, 'units', units, ...
                        'nSamples', nSamples);
                end

                % needed for reshape...
                if nDims == 1
                    nSamples(2) = 1;
                end

                this.data = reshape(squeeze(objectArray),nSamples);

                for iObject = 1:numel(this.data);
                    this.name = sprintf('%s %s', this.name, objectArray{iObject}.name);
                end
                this.info{end+1,1} = sprintf('Constructed from %s', this.name);
            end
        end

        % NOTE: Most of the methods are saved in separate function.m-files in this folder;
        %       except: constructor, delete, set/get methods for properties.

    end % methods

end

author: translationalneuromodeling

repo: tapas


text:
%% FUNCTION Least_SRMTL
%   Sparse Structure-Regularized Learning with Least Squares Loss.
%
%% OBJECTIVE
%   argmin_W { sum_i^t (0.5 * norm (Y{i} - X{i}' * W(:, i))^2)
%            + rho1 * norm(W*R, 'fro')^2 + rho2 * \|W\|_1}
%
%% R encodes structure relationship
%   1)Structure order is given by using [1 -1 0 ...; 0 1 -1 ...; ...]
%    e.g.: R=zeros(t,t-1);R(1:(t+1):end)=1;R(2:(t+1):end)=-1;
%   2)Ridge penalty term by setting: R = eye(t)
%   3)All related regularized: R = eye (t) - ones (t) / t
%
%% INPUT
%   X: {n * d} * t - input matrix
%   Y: {n * 1} * t - output matrix
%   R: regularization structure
%   rho1: structure regularization parameter
%   rho2: sparsity controlling parameter
%
%% OUTPUT
%   W: model: d * t
%   funcVal: function value vector.
%
%% LICENSE
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%   Copyright (C) 2011 - 2012 Jiayu Zhou and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 3, 2012.
%
%% RELATED PAPERS
%
% [1] Evgeniou, T. and Pontil, M. Regularized multi-task learning, KDD 2004
% [2] Zhou, J. Technical Report. http://www.public.asu.edu/~jzhou29/Software/SRMTL/CrisisEventProjectReport.pdf
%
%% RELATED FUNCTIONS
%  Logistic_SRMTL, init_opts

%% Code starts here
function [W, funcVal] = Least_SRMTL(X, Y, R, rho1, rho2, opts)

if nargin <5
    error('\n Inputs: X, Y, R, rho1, and rho2 should be specified!\n');
end
X = multi_transpose(X);

if nargin <6
    opts = [];
end

% initialize options.
opts=init_opts(opts);

if isfield(opts, 'rho_L2')
    rho_L2 = opts.rho_L2;
else
    rho_L2 = 0;
end

task_num  = length (X);
dimension = size(X{1}, 1);
funcVal = [];

% precomputation.
RRt = R * R';
XY = cell(task_num, 1);
W0_prep = [];
for t_idx = 1: task_num
    XY{t_idx} = X{t_idx}*Y{t_idx};
    W0_prep = cat(2, W0_prep, XY{t_idx});
end

% initialize a starting point
if opts.init==2
    W0 = zeros(dimension, task_num);
elseif opts.init == 0
    W0 = W0_prep;
else
    if isfield(opts,'W0')
        W0=opts.W0;
        if (nnz(size(W0)-[dimension, task_num]))
            error('\n Check the input .W0');
        end
    else
        W0=W0_prep;
    end
end

bFlag=0; % this flag tests whether the gradient step only changes a little

Wz= W0;
Wz_old = W0;

t = 1;
t_old = 0;

iter = 0;
gamma = 1;
gamma_inc = 2;

while iter < opts.maxIter
    alpha = (t_old - 1) /t;

    Ws = (1 + alpha) * Wz - alpha * Wz_old;

    % compute function value and gradients of the search point
    gWs  = gradVal_eval(Ws, rho1);
    Fs   = funVal_eval  (Ws, rho1);

    while true
        [Wzp l1c_wzp] = l1_projection(Ws - gWs/gamma, 2 * rho2 / gamma);
        Fzp = funVal_eval  (Wzp, rho1);

        delta_Wzp = Wzp - Ws;
        r_sum = norm(delta_Wzp, 'fro')^2;
        Fzp_gamma = Fs + trace(delta_Wzp' * gWs) + gamma/2 * norm(delta_Wzp, 'fro')^2;

        if (r_sum <=1e-20)
            bFlag=1; % this shows that, the gradient step makes little improvement
            break;
        end

        if (Fzp <= Fzp_gamma)
            break;
        else
            gamma = gamma * gamma_inc;
        end
    end

    Wz_old = Wz;
    Wz = Wzp;

    funcVal = cat(1, funcVal, Fzp + rho2 * l1c_wzp);

    if (bFlag)
        % fprintf('\n The program terminates as the gradient step changes the solution very small.');
        break;
    end

    % test stop condition.
    switch(opts.tFlag)
        case 0
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <= opts.tol)
                    break;
                end
            end
        case 1
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <=...
                        opts.tol* funcVal(end-1))
                    break;
                end
            end
        case 2
            if ( funcVal(end)<= opts.tol)
                break;
            end
        case 3
            if iter>=opts.maxIter
                break;
            end
    end

    iter = iter + 1;
    t_old = t;
    t = 0.5 * (1 + (1+ 4 * t^2)^0.5);

end

W = Wzp;

% private functions

    function [z l1_comp_val] = l1_projection (v, beta)
        % this projection calculates
        % argmin_z = \|z-v\|_2^2 + beta \|z\|_1
        % z: solution
        % l1_comp_val: value of l1 component (\|z\|_1)

        z = zeros(size(v));
        vp = v - beta/2;
        z (v> beta/2)  = vp(v> beta/2);
        vn = v + beta/2;
        z (v< -beta/2) = vn(v< -beta/2);

        l1_comp_val = sum(sum(abs(z)));
    end

    function [grad_W] = gradVal_eval(W, rho1)

        if opts.pFlag
            grad_W = zeros(size(W));
            parfor t_ii = 1:task_num
                XWi = X{t_ii}' * W(:,t_ii);
                XTXWi = X{t_ii}* XWi;
                grad_W(:, t_ii) = XTXWi - XY{t_ii};
                %grad_W = cat(2, grad_W, X{t_ii}*(X{t_ii}' * W(:,t_ii)-Y{t_ii}) );
            end
        else
            grad_W = [];
            for t_ii = 1:task_num
                XWi = X{t_ii}' * W(:,t_ii);
                XTXWi = X{t_ii}* XWi;
                grad_W = cat(2, grad_W, XTXWi - XY{t_ii});
                %grad_W = cat(2, grad_W, X{t_ii}*(X{t_ii}' * W(:,t_ii)-Y{t_ii}) );
            end
        end
        grad_W = grad_W + rho1 * 2 *  W * RRt...
            + rho_L2 * 2 * W;
    end

    function [funcVal] = funVal_eval (W, rho1)
        funcVal = 0;
        if opts.pFlag
            parfor i = 1: task_num
                funcVal = funcVal + 0.5 * norm (Y{i} - X{i}' * W(:, i))^2;
            end
        else
            for i = 1: task_num
                funcVal = funcVal + 0.5 * norm (Y{i} - X{i}' * W(:, i))^2;
            end
        end
        funcVal = funcVal + rho1 * norm(W*R, 'fro')^2 ...
            + rho_L2 * norm(W, 'fro')^2;
    end

end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function [B, model] = circulant_learning(X, para)
% CBE-opt
% X training data
% model.r optimized circulant vector
% model.bernoulli random bernoulli vector
% B CBE code for X

%N = size(X,1);
d = size(X,2);
rr = randn(1,d);
rr(rr > 0) = 1;
rr(rr <= 0 ) = -1;
% randomly flipping the sign
for i = 1:size(X,1)
    X(i,:) = X(i,:).*rr;
end
model.bernoulli = rr;

if ~isfield(para, 'bit')
   para.bit = size(X,2);
end

if ~isfield(para, 'iter')
   para.iter = 500;
end

if ~isfield(para, 'lambda')
   para.lambda = 1;
end

if ~isfield(para, 'verbose')
   para.verbose = 0;
end

% pre-compute the fft of X
fft_X = fft(X,[],2);

% pre-compute m
m = sum(real(fft_X).^2,1)+ sum(imag(fft_X).^2,1);
m = m'/d;

% initilization
if (~isfield(para, 'init_r'))
    r = randn(d,1);
else
    r = para.init_r;
end

fft_r = fft(r);
B = optimize_B(fft_r, fft_X, para);
ops = 100;
obj = inf;
iter = 0;

while(ops > 0.1 && iter <= para.iter)
    fprintf('iteration %d, obj = %f\n', iter, obj);
    if (para.verbose)
        fprintf('obj = %f \n', obj);
        fprintf('obj = %f \n', compute_obj_time(B, X, real(ifft(fft_r)), para.lambda));
        fft_B = fft(B, [], 2);
        fprintf('obj = %f \n', compute_obj_freq(fft_B, fft_X, fft_r, para.lambda));
    end

    [fft_r, obj_new] = optimize_R(B, fft_X, fft_r, m, para);
    B = optimize_B(fft_r, fft_X, para);
    ops = obj - obj_new;
    obj = obj_new;
    iter = iter + 1;
end

r = ifft(fft_r);
r = real(r);
model.r = r;

end

function obj = compute_obj_time(B, X, r, lambda)
   R = circulant(r, 1);
   obj = sum(sum((B - X*R').^2)) + lambda * sum(sum((R'*R - diag(ones(size(R,1), 1))).^2));
end

function obj = compute_obj_freq(fft_B, fft_X, fft_r, lambda)
    d = length(fft_r);
    Z = fft_B - repmat(fft_r.', size(fft_X,1) ,1).*fft_X;
    obj =  1/d * sum(sum(real(Z).^2 + imag(Z).^2));
    obj = obj + lambda*sum((real(fft_r).^2 + imag(fft_r).^2 - 1).^2);
end

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
function yhat=lfsmooth(varargin)
%
% a simple interface to locfit.
% output is a vector of smoothed values, at each data point.
% all locfit options, except evaluation structures, are valid.
%
% Example, to smooth a time series of observations,
%
% t = (1:100)';
% y = 2*sin(t/10) + normrnd(0,1,100,1);
% plot(t,y,'.');
% hold on;
% plot(t,lfsmooth(t,y,'nn',0.5));
% hold off;
%

% Minimal input validation
if nargin < 1
   error( 'At least one input argument required' );
end

fit = locfit(x,varargin{:},'module','simple');
yhat = fit.fit_points.fitted_values;

return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%SOLVEP3P  Finds an object pose from 3 3D-2D point correspondences
%
%     [rvecs, tvecs, solutions] = cv.solveP3P(objectPoints, imagePoints, cameraMatrix)
%     [...] = cv.solveP3P(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __objectPoints__ Array of object points in the object coordinate space,
%   1xNx3/Nx1x3 or Nx3 array, where `N=3` is the number of points, or cell
%   array of length `N=3` of 3-element vectors can be also passed here
%   `{[x1,y1,z1], [x2,y2,z2], [x3,y3,z3]}`.
% * __imagePoints__ Array of corresponding image points, 1xNx2/Nx1x2 or Nx2
%   array, where `N=3` is the number of points, or cell array of length `N=3`
%   of 2-element vectors can be also passed here `{[x1,y1], [x2,y2], [x3,y3]}`.
% * __cameraMatrix__ Input camera matrix `A = [fx 0 cx; 0 fy cy; 0 0 1]`.
%
% ## Output
% * __rvecs__ Output rotation vectors (see cv.Rodrigues) that, together with
%   `tvecs`, brings points from the model coordinate system to the camera
%   coordinate system. A P3P problem has up to 4 solutions.
% * __tvecs__ Output translation vectors.
% * __solutions__ number of solutions.
%
% ## Options
% * __DistCoeffs__ Input vector of distortion coefficients
%   `[k1,k2,p1,p2,k3,k4,k5,k6,s1,s2,s3,s4,taux,tauy]` of 4, 5, 8, 12 or 14
%   elements. If the vector is empty, the zero distortion coefficients are
%   assumed. default empty.
% * __Method__ Method for solving the P3P problem. One of the following:
%   * __P3P__ (default) Method is based on the paper [gao2003complete].
%   * __AP3P__ Method is based on the paper [Ke17].
%
% The function estimates the object pose given 3 object points, their
% corresponding image projections, as well as the camera matrix and the
% distortion coefficients.
%
% ## References
% [gao2003complete]:
% > X.S. Gao, X.R. Hou, J. Tang, H.F. Chang; "Complete Solution
% > Classification for the Perspective-Three-Point Problem",
% > IEEE Trans. on PAMI, vol. 25, No. 8, p. 930-943, August 2003.
%
% [Ke17]:
% > T. Ke, S. Roumeliotis; "An Efficient Algebraic Solution to the
% > Perspective-Three-Point Problem", IEEE Conference on Computer Vision and
% > Pattern Recognition (CVPR), 2017
% > [PDF](https://arxiv.org/pdf/1701.08237.pdf)
%
% See also: cv.solvePnP
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function hs = plotGaussEllipses( mus, Cs, rad )
% Plots 2D ellipses derived from 2D Gaussians specified by mus & Cs.
%
% USAGE
%  hs = plotGaussEllipses( mus, Cs, [rad] )
%
% INPUTS
%  mus     - k x 2 matrix of means
%  Cs      - 2 x 2 x k  covariance matricies
%  rad     - [2] Number of std to create the ellipse to
%
% OUTPUTS
%  hs      - handles to ellipses
%
% EXAMPLE
%  plotGaussEllipses( [ 10 10; 10 10 ], cat(3,eye(2),eye(2)*2) );
%
% See also PLOTELLIPSE, GAUSS2ELLIPSE
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if (nargin<3 || isempty(rad) ); rad=2; end
colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k']; nc = length(colors);

washeld = ishold; if (~washeld); hold('on'); end
hs = zeros( size(mus,1),1 );
for i=1:size( mus,1)
  [ cRow, ccol, ra, rb, phi ] = gauss2ellipse( mus(i,:), Cs(:,:,i), rad );
  hs(i)=plotEllipse( cRow, ccol, ra, rb, phi, colors( mod(i-1,nc)+1) );
end
if (~washeld); hold('off'); end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function lik = lik_binomial(varargin)
%LIK_BINOMIAL  Create a Binomial likelihood structure
%
%  Description
%    LIK = LIK_BINOMIAL creates Binomial likelihood structure.
%
%    The likelihood is defined as follows:
%                  __ n
%      p(y|f, z) = || i=1 [ p_i^(y_i)*(1-p_i)^(z_i-y_i)) *
%                           gamma(z_i+1)/(gamma(y_i+1)*gamma(z_i-y_i+1))]
%    where p_i = exp(f_i)/ (1+exp(f_i)) is the succes probability,
%    which is a function of the latent variable f_i and z is a
%    vector of numbers of trials.
%
%    When using Binomial likelihood you need to give the vector z
%    as an extra parameter to each function that requires y also.
%    For example, you should call gp_optim as follows
%      gp_optim(gp, x, y, 'z', z)
%
%  See also
%    GP_SET, LIK_*
%
% Copyright (c) 2009-2010 Jaakko Riihimäki & Jarno Vanhatalo
% Copyright (c) 2010-2011 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'LIK_BINOMIAL';
  ip.addOptional('lik', [], @isstruct);
  ip.parse(varargin{:});
  lik=ip.Results.lik;

  if isempty(lik)
    init=true;
    lik.type = 'Binomial';
  else
    if ~isfield(lik,'type') || ~isequal(lik.type,'Binomial')
      error('First argument does not seem to be a valid likelihood function structure')
    end
    init=false;
  end

  if init
    % Set the function handles to the subfunctions
    lik.fh.pak = @lik_binomial_pak;
    lik.fh.unpak = @lik_binomial_unpak;
    lik.fh.ll = @lik_binomial_ll;
    lik.fh.llg = @lik_binomial_llg;
    lik.fh.llg2 = @lik_binomial_llg2;
    lik.fh.llg3 = @lik_binomial_llg3;
    lik.fh.tiltedMoments = @lik_binomial_tiltedMoments;
    lik.fh.predy = @lik_binomial_predy;
    lik.fh.predprcty = @lik_binomial_predprcty;
    lik.fh.invlink = @lik_binomial_invlink;
    lik.fh.recappend = @lik_binomial_recappend;
  end

end

function [w,s,h] = lik_binomial_pak(lik)
%LIK_BINOMIAL_PAK  Combine likelihood parameters into one vector.
%
%  Description
%    W = LIK_BINOMIAL_PAK(LIK) takes a likelihood structure LIK
%    and returns an empty verctor W. If Binomial likelihood had
%    parameters this would combine them into a single row vector
%    W (see e.g. likelih_negbin). This is a mandatory subfunction
%    used for example in energy and gradient computations.
%
%  See also
%    LIK_NEGBIN_UNPAK, GP_PAK

  w = []; s = {}; h=[];
end

function [lik, w] = lik_binomial_unpak(lik, w)
%LIK_BINOMIAL_UNPAK  Extract likelihood parameters from the vector.
%
%  Description
%    W = LIK_BINOMIAL_UNPAK(W, LIK) Doesn't do anything.
%
%    If Binomial likelihood had parameters this would extracts
%    them parameters from the vector W to the LIK structure.
%    This is a mandatory subfunction used for example in energy
%    and gradient computations.
%
%  See also
%    LIK_BINOMIAL_PAK, GP_UNPAK

  lik=lik;
  w=w;

end

function ll = lik_binomial_ll(lik, y, f, z)
%LIK_BINOMIAL_LL  Log likelihood
%
%  Description
%    LL = LIK_BINOMIAL_LL(LIK, Y, F, Z) takes a likelihood
%    structure LIK, succes counts Y, numbers of trials Z, and
%    latent values F. Returns the log likelihood, log p(y|f,z).
%    This subfunction is needed when using Laplace approximation
%    or MCMC for inference with non-Gaussian likelihoods. This
%    subfunction is also used in information criteria (DIC, WAIC)
%    computations.
%
%  See also
%    LIK_BINOMIAL_LLG, LIK_BINOMIAL_LLG3, LIK_BINOMIAL_LLG2, GPLA_E

  if isempty(z)
    error(['lik_binomial -> lik_binomial_ll: missing z!'...
           'Binomial likelihood needs the expected number of   '...
           'occurrences as an extra input z. See, for         '...
           'example, lik_binomial and gp_optim.             ']);
  end

  expf = exp(f);
  p = expf ./ (1+expf);
  N = z;
  ll =  sum(gammaln(N+1)-gammaln(y+1)-gammaln(N-y+1)+y.*log(p)+(N-y).*log(1-p));
end

function llg = lik_binomial_llg(lik, y, f, param, z)
%LIK_BINOMIAL_LLG    Gradient of the log likelihood
%
%  Description
%    LLG = LIK_BINOMIAL_LLG(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, succes counts Y, numbers of trials Z and
%    latent values F. Returns the gradient of the log likelihood
%    with respect to PARAM. At the moment PARAM can be 'param' or
%    'latent'. This subfunction is needed when using Laplace
%    approximation or MCMC for inference with non-Gaussian
%    likelihoods.
%
%  See also
%    LIK_BINOMIAL_LL, LIK_BINOMIAL_LLG2, LIK_BINOMIAL_LLG3, GPLA_E

  if isempty(z)
    error(['lik_binomial -> lik_binomial_llg: missing z!'...
           'Binomial likelihood needs the expected number of   '...
           'occurrences as an extra input z. See, for         '...
           'example, lik_binomial and gp_optim.             ']);
  end

  switch param
    case 'latent'
      expf = exp(f);
      N = z;

      llg = y./(1+expf) - (N-y).*expf./(1+expf);
  end
end

function llg2 = lik_binomial_llg2(lik, y, f, param, z)
%LIK_BINOMIAL_LLG2  Second gradients of the log likelihood
%
%  Description
%    LLG2 = LIK_BINOMIAL_LLG2(LIK, Y, F, PARAM) takes a
%    likelihood structure LIK, succes counts Y, numbers of trials
%    Z, and latent values F. Returns the Hessian of the log
%    likelihood with respect to PARAM. At the moment PARAM can be
%    only 'latent'. G2 is a vector with diagonal elements of the
%    Hessian matrix (off diagonals are zero). This subfunction
%    is needed when using Laplace approximation or EP for inference
%    with non-Gaussian likelihoods.
%
%  See also
%    LIK_BINOMIAL_LL, LIK_BINOMIAL_LLG, LIK_BINOMIAL_LLG3, GPLA_E

  if isempty(z)
    error(['lik_binomial -> lik_binomial_llg2: missing z!'...
           'Binomial likelihood needs the expected number of    '...
           'occurrences as an extra input z. See, for          '...
           'example, lik_binomial and gp_optim.              ']);
  end

  switch param
    case 'latent'
      expf = exp(f);
      N = z;

      llg2 = -N.*expf./(1+expf).^2;
  end
end

function llg3 = lik_binomial_llg3(lik, y, f, param, z)
%LIK_BINOMIAL_LLG3  Third gradients of the log likelihood
%
%  Description
%    LLG3 = LIK_BINOMIAL_LLG3(LIK, Y, F, PARAM) takes a
%    likelihood structure LIK, succes counts Y, numbers of trials
%    Z and latent values F and returns the third gradients of the
%    log likelihood with respect to PARAM. At the moment PARAM
%    can be only 'latent'. G3 is a vector with third gradients.
%    This subfunction is needed when using Laplace appoximation
%    for inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_BINOMIAL_LL, LIK_BINOMIAL_LLG, LIK_BINOMIAL_LLG2, GPLA_E, GPLA_G

  if isempty(z)
    error(['lik_binomial -> lik_binomial_llg3: missing z!'...
           'Binomial likelihood needs the expected number of    '...
           'occurrences as an extra input z. See, for          '...
           'example, lik_binomial and gp_optim.              ']);
  end

  switch param
    case 'latent'
      expf = exp(f);
      N = z;
      llg3 = N.*(expf.*(expf-1))./(1+expf).^3;
  end
end

function [logM_0, m_1, sigm2hati1] = lik_binomial_tiltedMoments(lik, y, i1, sigma2_i, myy_i, z)
%LIK_BINOMIAL_TILTEDMOMENTS  Returns the marginal moments for EP algorithm
%
%  Description
%    [M_0, M_1, M2] = LIK_BINOMIAL_TILTEDMOMENTS(LIK, Y, I, S2,
%    MYY, Z) takes a likelihood structure LIK, succes counts Y,
%    numbers of trials Z, index I and cavity variance S2 and mean
%    MYY. Returns the zeroth moment M_0, mean M_1 and variance
%    M_2 of the posterior marginal (see Rasmussen and Williams
%    (2006): Gaussian processes for Machine Learning, page 55).
%    This subfunction is needed when using EP for inference with
%    non-Gaussian likelihoods.
%
%  See also
%    GPEP_E

%  if isempty(z)
%    error(['lik_binomial -> lik_binomial_tiltedMoments: missing z!'...
%           'Binomial likelihood needs the expected number of               '...
%           'occurrences as an extra input z. See, for                     '...
%           'example, lik_binomial and gp_optim.                         ']);
%  end

  yy = y(i1);
  N = z(i1);
  logM_0=zeros(size(yy));
  m_1=zeros(size(yy));
  sigm2hati1=zeros(size(yy));

  for i=1:length(i1)
    if isscalar(sigma2_i)
      sigma2ii = sigma2_i;
    else
      sigma2ii = sigma2_i(i);
    end

    % Create function handle for the function to be integrated
    % (likelihood * cavity) and useful integration limits
    [tf,minf,maxf]=init_binomial_norm(yy(i),myy_i(i),sigma2ii,N(i));

    % Integrate with quadrature
    RTOL = 1.e-6;
    ATOL = 1.e-10;
    [m_0, m_1(i), m_2] = quad_moments(tf,minf, maxf, RTOL, ATOL);
    if isnan(m_0)
      logM_0=NaN;
      return
    end
    sigm2hati1(i) = m_2 - m_1(i).^2;

    % If the second central moment is less than cavity variance
    % integrate more precisely. Theoretically for log-concave
    % likelihood should be sigm2hati1 < sigm2_i.
    if sigm2hati1(i) >= sigma2ii
      ATOL = ATOL.^2;
      RTOL = RTOL.^2;
      [m_0, m_1(i), m_2] = quad_moments(tf, minf, maxf, RTOL, ATOL);
      sigm2hati1(i) = m_2 - m_1(i).^2;
      %    if sigm2hati1 >= sigm2_i
      %      error('lik_binomial_tilted_moments: sigm2hati1 >= sigm2_i');
      %    end
    end
    logM_0(i) = log(m_0);
  end
end

function [lpy, Ey, Vary] = lik_binomial_predy(lik, Ef, Varf, yt, zt)
%LIK_BINOMIAL_PREDY  Returns the predictive mean, variance and density of y
%
%  Description
%    [LPY] = LIK_BINOMIAL_PREDY(LIK, EF, VARF YT, ZT)
%    Returns logarithm of the predictive density PY of YT, that is
%        p(yt | y, zt) = \int p(yt | f, zt) p(f|y) df.
%    This requires also the succes counts YT, numbers of trials ZT.
%    This subfunction is needed when computing posterior predictive
%    distributions for future observations.
%
%    [LPY, EY, VARY] = LIK_BINOMIAL_PREDY(LIK, EF, VARF) takes a
%    likelihood structure LIK, posterior mean EF and posterior
%    Variance VARF of the latent variable and returns the
%    posterior predictive mean EY and variance VARY of the
%    observations related to the latent variables. This subfunction
%    is needed when computing posterior predictive distributions for
%    future observations.
%
%
%  See also
%    GPEP_PRED, GPLA_PRED, GPMC_PRED

  if isempty(zt)
    error(['lik_binomial -> lik_binomial_predy: missing z!'...
           'Binomial likelihood needs the expected number of       '...
           'occurrences as an extra input z. See, for             '...
           'example, lik_binomial and gp_optim.                 ']);
  end

  if nargout > 1
  % Here we approximate the inverse logistic function by a linear combination of inverse probit function
  % logitinv(f) = a probit_inv(f/c1S) + (1 - a) probit_inv(f/c2s)
  % see Demidenko (2004). Mixed models: Theory and applications
      p1 = 0.4353;   p2 = 0.5647; % these are a and (1 - a)
      c1S = 2.2967^2;  c2S = 1.3017^2;

      z1 = Ef ./ sqrt(c1S + Varf);
      z2 = Ef ./ sqrt(c2S + Varf);

      % unconditional expectation
      p = (p1 .* normcdf(z1) + p2 .* normcdf(z2));
      Ey = zt .* p;

      % unconditional variance
      if any(zt ~= 1)
          nt = length(Ef);
          zz11 = zeros(nt, 1);
          zz12 = zeros(nt, 1);
          zz22 = zeros(nt, 1);

          for i = 1:length(Ef)
              zz11(i) = mvncdf([Ef(i) Ef(i)], [0 0], ...
                  [Varf(i) + c1S Varf(i); Varf(i) Varf(i) + c1S]);

              zz12(i) = mvncdf([Ef(i) Ef(i)], [0 0], ...
                  [Varf(i) + c1S Varf(i); Varf(i) Varf(i) + c2S]);

              zz22(i) = mvncdf([Ef(i) Ef(i)], [0 0], ...
                  [Varf(i) + c2S Varf(i); Varf(i) Varf(i) + c2S]);

          end
          Ey12 = p1^2 * zz11 + 2*p1*p2 * zz12 + p2^2 * zz22;

      else
          Ey12 = 0;

      end

      % unconditional variance
      Vary = zt .* (p - Ey12) + zt.^2 .* (Ey12 - p.^2);

  %     nt=length(Ef);
  %     Ey=zeros(nt,1);
  %     EVary = zeros(nt,1);
  %     VarEy = zeros(nt,1);
  %     for i1=1:nt
  %        ci = sqrt(Varf(i1));
  %       F  = @(x)zt(i1)./(1+exp(-x)).*norm_pdf(x,Ef(i1),sqrt(Varf(i1)));
  %        Ey(i1) = quadgk(F,Ef(i1)-6*ci,Ef(i1)+6*ci);
  %
  %       F2  = @(x)zt(i1)./(1+exp(-x)).*(1-1./(1+exp(-x))).*norm_pdf(x,Ef(i1),sqrt(Varf(i1)));
  %       EVary(i1) = quadgk(F2,Ef(i1)-6*ci,Ef(i1)+6*ci);
  %
  %       F3  = @(x)(zt(i1)./(1+exp(-x))).^2.*norm_pdf(x,Ef(i1),sqrt(Varf(i1)));
  %       VarEy(i1) = quadgk(F3,Ef(i1)-6*ci,Ef(i1)+6*ci) - Ey(i1).^2;
  %     end
  %     Vary = EVary+VarEy;
  end

  nt=length(yt);
  lpy=zeros(nt,1);
  if (size(Ef,2) > 1) && (size(Ef,2) > 1) && size(yt,2) == 1
    % Approximate integral with sum of grid points when using corrected
    % marginal posterior pf
    for i1=1:length(yt)
      py = arrayfun(@(f) exp(lik.fh.ll(lik, yt(i1), f, zt(i1))), Ef(i1,:));
      pf = Varf(i1,:)./sum(Varf(i1,:));
      lpy(i1) = log(sum(py.*pf));
    end
  else
    for i1=1:nt
      ci = sqrt(Varf(i1));
      F  = @(x)exp(gammaln(zt(i1)+1)-gammaln(yt(i1)+1)-gammaln(zt(i1)-yt(i1)+1) + yt(i1).*log(1./(1+exp(-x))) + (zt(i1)-yt(i1)).*log(1-(1./(1+exp(-x))))).*norm_pdf(x,Ef(i1),sqrt(Varf(i1)));
      lpy(i1) = log(quadgk(F,Ef(i1)-6*ci,Ef(i1)+6*ci));
    end
  end

end

function prctys = lik_binomial_predprcty(lik, Ef, Varf, zt, prcty)
%LIK_BINOMIAL_PREDPRCTY  Returns the percentiled of predictive density of y
%
%  Description
%    PRCTY = LIK_BINOMIAL_PREDPRCTY(LIK, EF, VARF YT, ZT)
%    Returns percentiles of the predictive density PY of YT, that is
%    This requires also the succes counts YT, numbers of trials ZT. This
%    subfunction is needed when using function gp_predprcty.
%
%  See also
%    GP_PREDPCTY

  if isempty(zt)
    error(['lik_binomial -> lik_binomial_predprcty: missing z!'...
           'Binomial likelihood needs the expected number of       '...
           'occurrences as an extra input z. See, for             '...
           'example, lik_binomial and gp_optim.                 ']);
  end

  opt=optimset('TolX',.5,'Display','off');
  nt=size(Ef,1);
  prctys = zeros(nt,numel(prcty));
  prcty=prcty/100;
  for i1=1:nt
    ci = sqrt(Varf(i1));
    for i2=1:numel(prcty)
      a=floor(fminbnd(@(a) (quadgk(@(f) binocdf(a,zt(i1),logitinv(f)).*norm_pdf(f,Ef(i1),ci),Ef(i1)-6*ci,Ef(i1)+6*ci,'AbsTol',1e-4)-prcty(i2)).^2,binoinv(prcty(i2),zt(i1),logitinv(Ef(i1)-1.96*ci)),binoinv(prcty(i2),zt(i1),logitinv(Ef(i1)+1.96*ci)),opt));
      if quadgk(@(f) binocdf(a,zt(i1),logitinv(f)).*norm_pdf(f,Ef(i1),ci),Ef(i1)-6*ci,Ef(i1)+6*ci,'AbsTol',1e-4)<prcty(i2)
        a=a+1;
      end
      prctys(i1,i2)=a;
    end
  end
end

function [df,minf,maxf] = init_binomial_norm(yy,myy_i,sigm2_i,N)
%INIT_BINOMIAL_NORM
%
%  Description
%    Return function handle to a function evaluating Binomial *
%    Gaussian which is used for evaluating (likelihood * cavity)
%    or (likelihood * posterior) Return also useful limits for
%    integration. This is private function for lik_binomial. This
%    subfunction is needed by subfunctions tiltedMoments and predy.
%
% See also
%   LIK_BINOMIAL_TILTEDMOMENTS, LIK_BINOMIAL_PREDY

% avoid repetitive evaluation of constant part
  ldconst = gammaln(N+1)-gammaln(yy+1)-gammaln(N-yy+1) - log(sigm2_i)/2 - log(2*pi)/2;
%   ldconst = log(factorial(N)/(factorial(yy)*factorial(N-yy))-log(sigm2_i)/2 -log(2*pi)/2;

 % Create function handle for the function to be integrated
  df = @binomial_norm;
 % use log to avoid underflow, and derivates for faster search
  ld = @log_binomial_norm;
  ldg = @log_binomial_norm_g;
  ldg2 = @log_binomial_norm_g2;

  % Set the limits for integration
  % Binomial likelihood is log-concave so the binomial_norm
  % function is unimodal, which makes things easier
  if yy==0 || yy==N
    % with yy==0 or yy==N the mode of the likelihood is not defined
    % use the mode of the Gaussian (cavity or posterior) as a first guess
    modef = myy_i;
  else
    % use precision weighted mean of the Gaussian approximation of the
    % binomial likelihood and Gaussian
    mean_app = log(yy./(N-yy));
    ld0=1/(1+exp(-mean_app));
    ld1=(1-ld0)*ld0;
    ld2=ld0-3*ld0^2+2*ld0^3;
    var_app=inv(-( yy*(ld2*ld0-ld1^2)/ld0^2 + (N-yy)*(ld2*(ld0-1)-ld1^2)/(ld0-1)^2 ));

    modef = (myy_i/sigm2_i + mean_app/var_app)/(1/sigm2_i + 1/var_app);
%     sigm_app = sqrt((1/sigm2_i + 1/var_app)^-1);
  end
  % find the mode of the integrand using Newton iterations
  % few iterations is enough, since the first guess in the right direction
  niter=3;       % number of Newton iterations
  mindelta=1e-6; % tolerance in stopping Newton iterations
  for ni=1:niter
      g = ldg(modef);
      h = ldg2(modef);
      delta=-g/h;
      modef=modef+delta;
      if abs(delta)<mindelta
          break
      end
  end
  if abs(delta)>1 || isinf(delta) || isnan(delta)
    % Newton algorithm didn't work properly so do binary search
    modef=myy_i;
    a=modef-5.*sqrt(sigm2_i); b=modef+5.*sqrt(sigm2_i); delta=1;
    while ldg(a)<0
      a=a-5.*sqrt(sigm2_i);
    end
    while ldg(b)>0
      b=b+5.*sqrt(sigm2_i);
    end
    while delta > 0.1
      modef=(a+b)/2;
      if ldg(modef) > 0
        a=modef;
      else
        b=modef;
      end
      delta=b-a;
    end
    h=ldg2(modef);
  end
  % integrand limits based on Gaussian approximation at mode
  modes=sqrt(-1/h);
  minf=modef-4*modes;
  maxf=modef+4*modes;
  modeld=ld(modef);
  iter=0;
  % check that density at end points is low enough
  lddiff=12; % min difference in log-density between mode and end-points
  minld=ld(minf);
  step=1;
  while minld>(modeld-lddiff)
    minf=minf-step*modes;
    minld=ld(minf);
    iter=iter+1;
    step=step*2;
    if iter>100
      error(['lik_binomial -> init_binomial_norm: ' ...
             'integration interval minimun not found ' ...
             'even after looking hard!'])
    end
  end
  maxld=ld(maxf);
  iter=0;
  step=1;
  while maxld>(modeld-lddiff)
    maxf=maxf+step*modes;
    maxld=ld(maxf);
    iter=iter+1;
    step=step*2;
    if iter>100
      error(['lik_binomial -> init_binomial_norm: ' ...
             'integration interval maximun not found ' ...
             'even after looking hard!'])
    end
  end

  function integrand = binomial_norm(f)
  % Logit * Gaussian
    integrand = exp(ldconst + yy*log(1./(1.+exp(-f)))+(N-yy)*log(1-1./(1.+exp(-f)))...
                   - 0.5 * (f-myy_i).^2./sigm2_i);
%     integrand = exp(ldconst ...
%                     +yy*log(x)+(N-yy)*log(1-x) ...
%                     -0.5*(f-myy_i).^2./sigm2_i);
    integrand(isnan(integrand)|isinf(integrand))=0;
  end

  function log_int = log_binomial_norm(f)
  % log(Binomial * Gaussian)
  % log_binomial_norm is used to avoid underflow when searching
  % integration interval

    log_int = ldconst + yy*log(1./(1.+exp(-f)))+(N-yy)*log(1-1./(1.+exp(-f)))...
                   - 0.5 * (f-myy_i).^2./sigm2_i;
%     log_int = ldconst ...
%               -log(1+exp(-yy.*f)) ...
%               -0.5*(f-myy_i).^2./sigm2_i;
    log_int(isnan(log_int)|isinf(log_int))=-Inf;
  end

  function g = log_binomial_norm_g(f)
  % d/df log(Binomial * Gaussian)
  % derivative of log_binomial_norm
    g = -(f-myy_i)./sigm2_i - exp(-f).*(N-yy)./((1+exp(-f)).^2.*(1-1./(1+exp(-f)))) ...
        + exp(-f).*yy./(1+exp(-f));
%     g = yy./(exp(f*yy)+1)...
%         + (myy_i - f)./sigm2_i;
  end

  function g2 = log_binomial_norm_g2(f)
  % d^2/df^2 log(Binomial * Gaussian)
  % second derivate of log_binomial_norm
    g2 = - (1+exp(2.*f)+exp(f).*(2+N*sigm2_i)./((1+exp(f))^2*sigm2_i));
%     a=exp(f*yy);
%     g2 = -a*(yy./(a+1)).^2 ...
%          -1/sigm2_i;
  end

end

function p = lik_binomial_invlink(lik, f, z)
%LIK_BINOMIAL_INVLINK  Returns values of inverse link function
%
%  Description
%    P = LIK_BINOMIAL_INVLINK(LIK, F) takes a likelihood structure LIK and
%    latent values F and returns the values of inverse link function P.
%    This subfunction is needed when using gp_predprctmu.
%
%     See also
%     LIK_BINOMIAL_LL, LIK_BINOMIAL_PREDY

  p = logitinv(f);
end

function reclik = lik_binomial_recappend(reclik, ri, lik)
%RECAPPEND  Append the parameters to the record
%
%  Description
%    RECLIK = GPCF_BINOMIAL_RECAPPEND(RECLIK, RI, LIK) takes a
%    likelihood record structure RECLIK, record index RI and
%    likelihood structure LIK with the current MCMC samples of
%    the parameters. Returns RECLIK which contains all the old
%    samples and the current samples from LIK. This subfunction
%    is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC

  if nargin == 2
    reclik.type = 'Binomial';

    % Set the function handles
    reclik.fh.pak = @lik_binomial_pak;
    reclik.fh.unpak = @lik_binomial_unpak;
    reclik.fh.ll = @lik_binomial_ll;
    reclik.fh.llg = @lik_binomial_llg;
    reclik.fh.llg2 = @lik_binomial_llg2;
    reclik.fh.llg3 = @lik_binomial_llg3;
    reclik.fh.tiltedMoments = @lik_binomial_tiltedMoments;
    reclik.fh.invlink = @lik_binomial_invlink;
    reclik.fh.predprcty = @lik_binomial_predprcty;
    reclik.fh.predy = @lik_binomial_predy;
    reclik.fh.recappend = @likelih_binomial_recappend;
    return
  end

end

author: gpstuff-dev

repo: gpstuff


text:
function X = flushmidfactors(X)
for i = 1:length(X.midfactors)
    if isa(X.midfactors{i},'sdpvar')
        X.midfactors{i} = flush(X.midfactors{i});
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
classdef TestLUCID
    %TestLUCID

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','tsukuba.png');
    end

    methods (Static)
        function test_compute_img
            obj = cv.LUCID('LucidKernel',1, 'BlurKernel',2);
            typename = obj.typeid();
            ntype = obj.defaultNorm();

            img = imread(TestLUCID.im);
            kpts = cv.FAST(img, 'Threshold',20);
            [desc, kpts2] = obj.compute(img, kpts);
            validateattributes(kpts2, {'struct'}, {'vector'});
            assert(all(ismember(fieldnames(kpts), fieldnames(kpts2))));
            validateattributes(desc, {obj.descriptorType()}, ...
                {'size',[numel(kpts2) obj.descriptorSize()]});
        end

        function test_compute_imgset
            img = imread(TestLUCID.im);
            imgs = {img, img};
            kpts = cv.FAST(img, 'Threshold',20);
            kpts = {kpts, kpts};

            obj = cv.LUCID();
            [descs, kpts] = obj.compute(imgs, kpts);
            validateattributes(kpts, {'cell'}, {'vector'});
            validateattributes(descs, {'cell'}, {'vector', 'numel',numel(imgs)});
            cellfun(@(d,k) validateattributes(d, {obj.descriptorType()}, ...
                {'size',[numel(k) obj.descriptorSize()]}), descs, kpts);
        end

        function test_error_1
            try
                cv.LUCID('foobar');
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end


text:
function TT_Tools_demo(project_file)
% TT_Tools_demo(project_file)
%
% This function demonstrates basic functionality of the Tracking Tools API
% from Natural Point (to be used with Optitrack). Before use it is
% essential to have calibrated cameras and created the desirable
% trackable/s and save these in a project file. The project file is then
% passed to the function and the library loaded if need be before
% attempting to plot any data
%
% Input - project_file - a string containing path and filename for the
%               project file that is to be used (default value assigned if
%               not assigned)
%
% Written by Glen Lichtwark, University of Queensland, Australia
% Last updated: 22nd Jan 2010
% Please acknowledge in any academic papers which may utilise this code

if nargin < 1
    project_file = 'C:\Optitrack\LabSpace.ttp';
end

% load the NPTrackingTools library if it is not already loaded
if ~libisloaded('NPTrackingTools')

addpath('C:\Program Files\NaturalPoint\TrackingTools\lib'); % change if necessary
addpath('C:\Program Files\NaturalPoint\TrackingTools\inc'); % change if necessary

[notfound,warnings]=loadlibrary('NPTrackingTools','NPTrackingTools.h');
end

% libfunctionsview NPTrackingTools --> use this to see available functions

% initialise cameras
calllib('NPTrackingTools', 'TT_Initialize');

% load the project file which sets up cameras correctly
calllib('NPTrackingTools', 'TT_LoadProject', project_file);

%define the outputs types from TT_TrackableLocation function
X = 0;Y = 0;Z = 0;
qx = 0;qy = 0;qz = 0;qw = 0;
yaw = 0;pitch = 0;roll = 0;

figure(1)
clf
set(gcf,'Position',[127 73 933 602])

TC = calllib('NPTrackingTools','TT_TrackableCount');

TrackableNum = TC-1; % change this value to view different trackable object (starts at 0)

%loop through and plot the marker positions using frame results
for i = 1:500

    M = [];

    %update frame and get time stamp
    calllib('NPTrackingTools', 'TT_UpdateSingleFrame');
    D.T(i) = calllib('NPTrackingTools', 'TT_FrameTimeStamp');

    %find out how many markers are visible store data for X Y Z coordinates
    %of each
    marker_count = calllib('NPTrackingTools', 'TT_FrameMarkerCount');
    for j = 1:marker_count
        M(j,1) = calllib('NPTrackingTools', 'TT_FrameMarkerX',j-1);
        M(j,2) = calllib('NPTrackingTools', 'TT_FrameMarkerY',j-1);
        M(j,3) = calllib('NPTrackingTools', 'TT_FrameMarkerZ',j-1);
    end
    D.dat(i) = {M};

    % find the location of any trackable and plot the XYZ position on one
    % plot and euler angles on another

    [X,Y,Z,qx,qy,qz,qw,yaw,pitch,roll] = calllib('NPTrackingTools', 'TT_TrackableLocation',TrackableNum,X,Y,Z,qx,qy,qz,qw,yaw,pitch,roll);
    D.trans_dat(i,:) = [X Y Z];
    D.rot_dat(i,:) = [yaw pitch roll];
    % plot data --> note that this slows the frame rate considerably so
    % only do it every now and again
    if rem(i,4) == 0
        subplot(2,2,2), plot(D.T-D.T(1), D.trans_dat)
        xlabel('Time (s)?')
        ylabel('Object Position (m)?')
        subplot(2,2,4), plot(D.T-D.T(1), D.rot_dat)
        xlabel('Time (s)?')
        ylabel('Object Orientation (deg)?')

        % make 3D plot of marker positions and the trackable position
        if ~isempty(M)
            subplot(1,2,1), plot3(M(:,1),M(:,2),M(:,3),'ko',X,Y,Z,'ro');axis equal
            axis([-0.4 0.4 -0.4 0.4 -0.4 0.4])
            xlabel('X')
            ylabel('Y')
            zlabel('Z')
        end
        drawnow
    end

end

calllib('NPTrackingTools', 'TT_Shutdown')

author: Sable

repo: mcbench-benchmarks