text:
function v = kspace3d(v,M)
% 3D rigid body transformation performed as shears in 1D Fourier space.
% FORMAT v1 = kspace3d(v,M)
% Inputs:
% v - the image stored as a 3D array.
% M - the rigid body transformation matrix.
% Output:
% v - the transformed image.
%
% The routine is based on the excellent papers:
% R. W. Cox and A. Jesmanowicz (1999)
% Real-Time 3D Image Registration for Functional MRI
% Submitted to MRM (April 1999) and avaliable from:
% http://varda.biophysics.mcw.edu/~cox/index.html.
% and:
% W. F. Eddy, M. Fitzgerald and D. C. Noll (1996)
% Improved Image Registration by Using Fourier Interpolation
% Magnetic Resonance in Medicine 36(6):923-931
%_______________________________________________________________________

% taken from SPM8

[S0,S1,S2,S3] = shear_decomp(M);

d  = [size(v) 1 1 1];
g = 2.^ceil(log2(d));
if any(g~=d),
    tmp = v;
    v   = zeros(g);
    v(1:d(1),1:d(2),1:d(3)) = tmp;
    clear tmp;
end

% XY-shear
tmp1 = -sqrt(-1)*2*pi*([0:((g(3)-1)/2) 0 (-g(3)/2+1):-1])/g(3);
for j=1:g(2),
    t        = reshape( exp((j*S3(3,2) + S3(3,1)*(1:g(1)) + S3(3,4)).'*tmp1) ,[g(1) 1 g(3)]);
    v(:,j,:) = real(ifft(fft(v(:,j,:),[],3).*t,[],3));
end

% XZ-shear
tmp1 = -sqrt(-1)*2*pi*([0:((g(2)-1)/2) 0 (-g(2)/2+1):-1])/g(2);
for k=1:g(3),
    t        = exp( (k*S2(2,3) + S2(2,1)*(1:g(1)) + S2(2,4)).'*tmp1);
    v(:,:,k) = real(ifft(fft(v(:,:,k),[],2).*t,[],2));
end

% YZ-shear
tmp1 = -sqrt(-1)*2*pi*([0:((g(1)-1)/2) 0 (-g(1)/2+1):-1])/g(1);
for k=1:g(3),
    t        = exp( tmp1.'*(k*S1(1,3) + S1(1,2)*(1:g(2)) + S1(1,4)));
    v(:,:,k) = real(ifft(fft(v(:,:,k),[],1).*t,[],1));
end

% XY-shear
tmp1 = -sqrt(-1)*2*pi*([0:((g(3)-1)/2) 0 (-g(3)/2+1):-1])/g(3);
for j=1:g(2),
    t        = reshape( exp( (j*S0(3,2) + S0(3,1)*(1:g(1)) + S0(3,4)).'*tmp1) ,[g(1) 1 g(3)]);
    v(:,j,:) = real(ifft(fft(v(:,j,:),[],3).*t,[],3));
end

if any(g~=d), v = v(1:d(1),1:d(2),1:d(3)); end
return;
%_______________________________________________________________________

%_______________________________________________________________________

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function w = ymdf_to_weekday_islamic_a ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_WEEKDAY_ISLAMIC_A returns the weekday of an Islamic A YMDF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, D, real F, the YMDF date.
%
%    Output, integer W, is the week day number of the date, with
%    1 for Sunday, through 7 for Saturday.
%
  jed = ymdf_to_jed_islamic_a ( y, m, d, f );

  [ w, f2 ] = jed_to_weekday ( jed );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function value = p20_r8 ( action, name, value )

%*****************************************************************************80
%
%% P20_R8 sets or gets R8 parameters for problem 20.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION,
%    'D' sets the parameter to its default value;
%    'G' gets a parameter.
%    'R' sets the parameter to a random value.
%    'S' sets a parameter,
%
%    Input, string NAME, the name of the parameter.
%    'A' is the value of all lower integral bounds.
%    'B' is the value of all upper integral bounds.
%    'P' is the value of the exponent in the integrand.
%
%    Input/output, real VALUE.
%    * If ACTION = 'S', then VALUE is an input quantity, and is the
%      new value to be assigned to NAME.
%    * If ACTION = 'G' or 'R', then VALUE is an output quantity,
%      and is the current value of NAME.
%
  persistent a;
  persistent b;
  persistent p;

  if ( action == 'D' | action == 'd' )

    if ( name == 'A' | name == 'a' | name == '*' )
      a = 0.0;
    end

    if ( name == 'B' | name == 'b' | name == '*' )
      b = 1.0;
    end

    if ( name == 'P' | name == 'p' | name == '*' )
      p = 2.0;
    end

  elseif ( action == 'G' | action == 'g' )

    if ( name == 'A' | name == 'a' )
      value = a;
    elseif ( name == 'B' | name == 'b' )
      value = b;
    elseif ( name == 'P' | name == 'p' )
      value = p;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P20_R8 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P20_R8 - Fatal error!' );
    end

  elseif ( action == 'R' | action == 'r' )

    if ( name == 'A' | name == 'a' )
      a = rand ( 1, 1 );
      value = a;
    elseif ( name == 'B' | name == 'b' )
      b = rand ( 1, 1 );
      value = b;
    elseif ( name == 'P' | name == 'p' )
      p = rand ( 1, 1 );
      value = p;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P20_R8 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P20_R8 - Fatal error!' );
    end

  elseif ( action == 'S' | action == 's' )

    if ( name == 'A' | name == 'a' )
      a = value;
    elseif ( name == 'B' | name == 'b' )
      b = value;
    elseif ( name == 'P' | name == 'p' )
      p = value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P20_R8 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P20_R8 - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P20_R8 - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action = "%s".\n', action );
    error ( 'P20_R8 - Fatal error!' );

  end

  return
end


text:
clear, clc;

% This is an example for running the function mcLeastC
%
%  Problem:
%
%  min  1/2 || A x - y||^2
%  s.t. sum_j ||x^j||_q <=z
%
%  The current program only implements q=2
%
%  x is grouped into k groups according to opts.ind
%  The indices of x_j in x is (ind(j)+1):ind(j+1)
%
% For detailed description of the function, please refer to the Manual.
%
%% ------------   History --------------------
% First version on August 10, 2009.
%
% September 5, 2009: adaptive line search is added
%
% For any problem, please contact Jun Liu (j.liu@asu.edu)

cd ..
cd ..

root=cd;
addpath(genpath([root '/SLEP']));
                     % add the functions in the folder SLEP to the path

% change to the original folder
cd Examples/L1Lq;

m=100;  n=100;       % The data matrix is of size m x n
k=10;                % the number of classes (tasks)
q=2;                 % the value of q in the L1/Lq regularization
z=50;             % the regularization parameter
randNum=1;           % a random number

% ---------------------- generate random data ----------------------
randn('state',(randNum-1)*3+1);
A=randn(m,n);        % the data matrix

randn('state',(randNum-1)*3+2);
xOrin=randn(n,k);

randn('state',(randNum-1)*3+3);
noise=randn(m,k);
y=A*xOrin +...
    noise*0.01;      % the response

%----------------------- Set optional items -----------------------
opts=[];

% Starting point
opts.init=2;        % starting from a zero point

% Termination
opts.tFlag=5;       % run .maxIter iterations
opts.maxIter=100;   % maximum number of iterations

% Normalization
opts.nFlag=0;       % without normalization

% Group Property
opts.q=q;           % set the value for q

%----------------------- Run the code mcLeastC -----------------------
fprintf('\n lFlag=0 \n');
opts.mFlag=0;       % treating it as compositive function
opts.lFlag=0;       % Nemirovski's line search
tic;
[x1, funVal1, ValueL1]= mcLeastC(A, y, z, opts);
toc;

opts.maxIter=200;

fprintf('\n lFlag=1 \n');
opts.mFlag=1;       % smooth reformulation
opts.lFlag=1;       % adaptive line search
opts.tFlag=2; opts.tol= funVal1(end);
tic;
[x2, funVal2, ValueL2]= mcLeastC(A, y, z, opts);
toc;

figure;
plot(funVal1,'-r');
hold on;
plot(funVal2,'--b');
legend('lFlag=0', 'lFlag=1');
xlabel('Iteration (i)');
ylabel('The objective function value');

% % --------------------- compute the pathwise solutions ----------------
% opts.fName='mcLeastC';    % set the function name to 'mcLeastC'
% Z=[0.9, 0.8, 0.5, 0.3];   % set the parameters
%
% % run the function pathSolutionLeast
% fprintf('\n Compute the pathwise solutions, please wait...');
% X=pathSolutionLeast(A, y, Z, opts);

author: zzstefan

repo: BrainNetClass


text:
function [data,p]=create_data(p)
if nargin==0
    p=struct();
end

gpr=load('gpr.dat');

p.zss = exp(mean(gpr(:,1)));
p.paiss = exp(mean(gpr(:,2)));
p.beta = (p.zss*p.paiss)/(exp(mean(gpr(:,3))));

start_date='1983Q1';
gt = gpr(:,1) - log(p.zss);
pit = gpr(:,2) - log(p.paiss);
rt = gpr(:,3) - log(p.zss) + log(p.beta) - log(p.paiss);

data=ts(start_date,[gt,pit,rt],{'GHAT','PAIHAT','RHAT'});
data=pages2struct(data);

author: jmaih


text:
classdef VARtype < uint8
    %VARTYPE enmieration class to define the differnt type of VAR problems
    %that can be solved in BEAR toolbox

    enumeration
        OLS   (1) % OLS VAR
        BVAR  (2) % Bayesian VAR
        PANEL (4) % Panel BVAR
        SV    (5) % Stochastic Volatility BVAR
        TVP   (6) % Time-Varying Parameter BVAR
        MFVAR (7) % Time-Varying Parameter BVAR
    end

end

author: european-central-bank

repo: BEAR-toolbox


text:
% Copyright (C) Daphne Koller, Stanford University, 2012

function [MEU OptimalDecisionRule] = OptimizeMEU( I )

  % Inputs: An influence diagram I with a single decision node and a single utility node.
  %         I.RandomFactors = list of factors for each random variable.  These are CPDs, with
  %              the child variable = D.var(1)
  %         I.DecisionFactors = factor for the decision node.
  %         I.UtilityFactors = list of factors representing conditional utilities.
  % Return value: the maximum expected utility of I and an optimal decision rule
  % (represented again as a factor) that yields that expected utility.

  % We assume I has a single decision node.
  % You may assume that there is a unique optimal decision.
  D = I.DecisionFactors(1);
  DF = CalculateExpectedUtilityFactor(I);
  MEU = 0;
  D.val = D.val*0;
  ca=D.card(1);
  l = length(D.val)/ca;
  for i = 1:l
	  start = i*ca-ca+1;
  endd = i*ca;
  [x,y] = max(DF.val(start:endd));
	MEU+=x;
	D.val(start+y-1) = 1;
end
OptimalDecisionRule = D;
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  % YOUR CODE HERE...
  %
  % Some other information that might be useful for some implementations
  % (note that there are multiple ways to implement this):
  % 1.  It is probably easiest to think of two cases - D has parents and D
  %     has no parents.
  % 2.  You may find the Matlab/Octave function setdiff useful.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function x=v_pcma2lin(p,m,s)
%V_PCMU2LIN Convert A-law PCM to linear X=(P,M,S)
%	lin = v_pcma2lin(pcma,m,s) where pcma contains a vector or matrix
%	of A-law values in the range 0 to 255.
%	No checking is performed to see that numbers are in this range.
%
%	Input values are exclusive ored with m (default=85)
%
%	Output values are divided by the scale factor s:
%
%		   s		Output Range
%
%		   1		+-4032	(integer values)
%		2017.396342	+-1.998616 (default)
%		4032		+-1
%		4096		+-0.984375 (+-1 nominal full scale)
%
%	The default value of s is 2017.396342 which equals
%	sqrt((1120^2 + 2624^2)/2). This factor follows ITU standard G.711 and
%	the sine wave with PCM-A values [225 244 244 225 97 116 116 97]
%	has a mean square value of unity corresponding to 0 dBm0.

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: v_pcma2lin.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<3
  t=4.95688418E-4;
  if nargin<2 m=85; end
else
  t=1/s;
end

if m q=bitxor(p,m); else q=p; end;
k=rem(q,16);
g=floor(q/128);
e=(q-k-128*g)/16;
f=(abs(e-1)-e+1)/2;
x=(2*g-1).*(pow2(k+16.5,e)+f.*(k-15.5))*t;

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function snew = structprefixfieldnames(s,pfix)
assert(isstruct(s));
assert(ischar(pfix));

fnames = fieldnames(s);
pfixfnames = cellfun(@(x)[pfix x],fnames,'uni',0);

c = struct2cell(s);
snew = cell2struct(c,pfixfnames,1);

assert(isequal(size(s),size(snew)));

author: kristinbranson

repo: JAABA


text:
function varargout = histosettingGUI(varargin)
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @histosettingGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @histosettingGUI_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 histosettingGUI is made visible.
function histosettingGUI_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = hObject;

try
        def       = varargin{1};
        binvalue  = def{1};
        normhisto = def{2};
        chisto    = def{3};
        fitnormd  = def{4};
        cfitnorm  = def{5};
        nvalue   = def{6};
catch
        binvalue  = 'auto';
        normhisto = 0;
        chisto    = [1 0.5 0.2];
        fitnormd  = 0;
        cfitnorm  = [1 0 0];
        nvalue   = 100;
end
if isempty(binvalue)
        binvalue = 'auto';
end
if isempty(normhisto)
        normhisto = 0;
end
if isempty(chisto)
        chisto = [1 0.5 0.2];
end
if isempty(cfitnorm)
        cfitnorm = [1 0 0];
end
if isempty(fitnormd)
        fitnormd = 0;
end
if isempty(nvalue)
        nvalue = 100;
end

handles.binvalue  = binvalue;
handles.normhisto = normhisto;
handles.chisto    = chisto;
handles.cfitnorm  = cfitnorm;
handles.fitnormd  = fitnormd;
handles.nvalue    = nvalue;

%
% Name & version
%
version = geterplabversion;
set(handles.gui_chassis,'Name', ['ERPLAB ' version '   -   Histogram Settings GUI'],'WindowStyle','modal')

%
% Color GUI
%
handles = painterplab(handles);

%
% Set font size
%
handles = setfonterplab(handles);

% Update handles structure
guidata(hObject, handles);

if ischar(binvalue)
        set(handles.edit_nbinh, 'String', binvalue)
else
        set(handles.edit_nbinh, 'String', vect2colon(binvalue,'Delimiter', 'off'))
end
set(handles.checkbox_normhist, 'Value', normhisto)
set(handles.checkbox_fitnormal,'Value', fitnormd)

% if normhisto
%         set(handles.checkbox_fitnormal, 'Enable', 'on')
%         set(handles.checkbox_fitnormal,'Value', fitnormd)
%         set(handles.pushbutton_colornormdist, 'Enable', 'on')
% else
%         set(handles.checkbox_fitnormal, 'Value', 0)
%         set(handles.checkbox_fitnormal, 'Enable', 'off')
%         set(handles.pushbutton_colornormdist, 'Enable', 'off')
% end

% UIWAIT makes histosettingGUI wait for user response (see UIRESUME)
uiwait(handles.gui_chassis);

% --- Outputs from this function are returned to the command line.
function varargout = histosettingGUI_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;

% The figure can be deleted now
delete(handles.gui_chassis);
pause(0.1)

%----------------------------------------------------------------
function edit_nbinh_Callback(hObject, eventdata, handles)

%----------------------------------------------------------------
function edit_nbinh_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%----------------------------------------------------------------
function pushbutton_cancel_Callback(hObject, eventdata, handles)
handles.output= [];
% Update handles structure
guidata(hObject, handles);
uiresume(handles.gui_chassis);

%----------------------------------------------------------------
function pushbutton_OK_Callback(hObject, eventdata, handles)

binvalue  = get(handles.edit_nbinh, 'String');
if isempty(binvalue)
        msgboxText =  'You have to specify one value at list!';
        title = 'Error: histosettingGUI';
        errorfound(sprintf(msgboxText), title);
        return
end
binvalue = strrep(binvalue,'''','');
if ~strcmpi(binvalue,'auto') && isempty(str2num(binvalue))
        msgboxText =  ['%s is not a valid input.\n\n'...
                'Enter a single value, a monotonically non-decreasing vector, or ''auto''.\n'...
                'Or just click the "suggest" button.'];
        title = 'Error: histosettingGUI';
        errorfound(sprintf(msgboxText, binvalue), title);
        return

elseif ~strcmpi(binvalue,'auto') && ~isempty(str2num(binvalue))
        binvalue = str2num(binvalue);
end

normhisto = get(handles.checkbox_normhist, 'Value');
chisto    = handles.chisto;
fitnormd  = get(handles.checkbox_fitnormal, 'Value');
cfitnorm  = handles.cfitnorm;

erpworkingmemory('BinHisto', binvalue);
erpworkingmemory('NormHisto', normhisto);
erpworkingmemory('FitNormd', fitnormd);

handles.output = {binvalue, normhisto, chisto, fitnormd, cfitnorm};
% Update handles structure
guidata(hObject, handles);
uiresume(handles.gui_chassis);

%----------------------------------------------------------------
function pushbutton_colorhist_Callback(hObject, eventdata, handles)
chisto = handles.chisto;
chisto = uisetcolor(chisto,'Histogram color') ;
erpworkingmemory('HistoColor', chisto);

handles.chisto = chisto;
% Update handles structure
guidata(hObject, handles);

%----------------------------------------------------------------
function checkbox_normhist_Callback(hObject, eventdata, handles)
% if get(hObject,'Value')
%         set(handles.checkbox_fitnormal, 'Enable', 'on')
%         set(handles.pushbutton_colornormdist, 'Enable', 'on')
% else
%         set(handles.checkbox_fitnormal, 'Value', 0)
%         set(handles.checkbox_fitnormal, 'Enable', 'off')
%         set(handles.pushbutton_colornormdist, 'Enable', 'off')
% end

%----------------------------------------------------------------
function checkbox_fitnormal_Callback(hObject, eventdata, handles)

%----------------------------------------------------------------
function pushbutton_colornormdist_Callback(hObject, eventdata, handles)
cfitnorm = handles.cfitnorm;
cfitnorm = uisetcolor(cfitnorm,'Curve color') ;
erpworkingmemory('FnormColor', cfitnorm);

handles.cfitnorm = cfitnorm;
% Update handles structure
guidata(hObject, handles);

%----------------------------------------------------------------
function pushbutton_autobinhisto_Callback(hObject, eventdata, handles)
nvalue = handles.nvalue;
binvalue = round(sqrt(nvalue));
handles.binvalue = binvalue;
% Update handles structure
guidata(hObject, handles);

set(handles.edit_nbinh, 'String', num2str(binvalue));

%--------------------------------------------------------------------------
function gui_chassis_CloseRequestFcn(hObject, eventdata, handles)
if isequal(get(handles.gui_chassis, 'waitstatus'), 'waiting')
        handles.output= [];
        % Update handles structure
        guidata(hObject, handles);
        uiresume(handles.gui_chassis);
else
        % The GUI is no longer waiting, just close it
        delete(handles.gui_chassis);
end

% chisto   = erpworkingmemory('HistoColor'); % histogram color
% cfitnorm = erpworkingmemory('FnormColor'); % line color for fitted normal distribution
% fitnormd = erpworkingmemory('FitNormd'); % select window measurement by mouse option

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function [W] = hpe2W(hpe)
% Convert power from electrical horsepower to watts.
% Chad A. Greene 2012
W = hpe*746;

author: Sable

repo: mcbench-benchmarks


text:
function title = p34_title ( )

%*****************************************************************************80
%
%% P34_TITLE returns the title of problem P34.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'Problem 34, The Linearized Damped Pendulum';

  return
end


text:
function percentCorrect = prtEvalPercentCorrect(classifier,dataSet,nFolds)
% prtEvalPercentCorrect    Calculate percent correct of a classification operation
%
%   percentCorrect = prtEvalPercentCorrect(prtClassifier, prtDataSet)
%   returns the percentage of correctly classified elements of prtDataSet
%   when classifier by prtClassifier. prtDataSet must be a labeled,
%   prtDataSetStandard object. prtClassifier must be a prtClass object.
%
%   percentCorrect = prtEvalPercentCorrect(prtClassifier, prtDataSet,
%   nFolds)  returns the percentage of correctly classified elements of
%   prtDataSet when classifier by prtClassifier with K-fold
%   cross-validation. prtDataSet must be a labeled, prtDataSetStandard
%   object. prtClassifier must be a prtClass object. nFolds is the number
%   of folds in the K-fold cross-validation.
%
%   percentCorrect = prtEvalPercentCorrect(prtClassifier, prtDataSet,
%   xValInds) same as above, but use crossValidation with specified indices
%   instead of random folds.
%
%   Example:
%   dataSet = prtDataGenSpiral;
%   classifier = prtClassDlrt;
%   pc =  prtEvalPercentCorrect(classifier, dataSet)
%
%   See Also: prtEvalPdAtPf, prtEvalPfAtPd, prtEvalAuc,
%   prtEvalMinCost

if nargin < 3 || isempty(nFolds)
    nFolds = 1;
end
results = prtUtilEvalParseAndRun(classifier,dataSet,nFolds);

%(note: can't check results.nFeatures here any more...)
if dataSet.nClasses == 2 %binary classifier
    [pf,pd] = prtScoreRoc(results.getObservations,dataSet.getTargets);
    pe = prtUtilPfPd2Pe(pf,pd);
    minPe = min(pe);
    percentCorrect = 1-minPe;
else
    %Note, this is a hack; we need to fix this.
    if isa(classifier,'prtAlgorithm')
        guess = results.getObservations;
    else
        if classifier.includesDecision
            guess = results.getObservations;
        else
           %Naive MAP decision:
           % [twiddle,guess] = max(results.getObservations,[],2); %#ok<ASGLU>
           %Naive MAP decision:
           prtMap = prtDecisionMap;
           prtMap = train(prtMap,results);
           results = prtMap.run(results);
           guess = results.getObservations;
        end
    end
    percentCorrect = prtScorePercentCorrect(guess,dataSet.getTargets);
    %confusionMatrix = prtScoreConfusionMatrix(guess,dataSet.getTargets);
    %percentCorrect = prtUtilConfusion2PercentCorrect(confusionMatrix);
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [t, storedb] = getLinesearch(problem, x, d, storedb)
% Returns a hint for line-search algorithms.
%
% function [t, storedb] = getLinesearch(problem, x, d, storedb)
%
% For a line-search problem at x along the tangent direction d, computes
% and returns t such that retracting t*d at x yields a good point around
% where to look for a line-search solution. That is: t is a hint as to "how
% far to look" along the line.
%
% The cache database storedb is passed along, possibly modified and
% returned in the process.
%
% See also: canGetLinesearch

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, July 17, 2014.
% Contributors:
% Change log:

    if isfield(problem, 'linesearch')
    %% Compute the line-search hint function using linesearch.

		is_octave = exist('OCTAVE_VERSION', 'builtin');
		if ~is_octave
			narg = nargin(problem.linesearch);
		else
			narg = 3;
		end

        % Check whether the linesearch function wants to deal with the
        % store structure or not.
        switch narg
            case 2
                t = problem.linesearch(x, d);
            case 3
                % Obtain, pass along, and save the store structure
                % associated to this point.
                store = getStore(problem, x, storedb);
                [t, store] = problem.linesearch(x, d, store);
                storedb = setStore(problem, x, storedb, store);
            otherwise
                up = MException('manopt:getLinesearch:badfun', ...
                    'linesearch should accept 2 or 3 inputs.');
                throw(up);
        end

    else
    %% Abandon computing the line-search function.

        up = MException('manopt:getLinesearch:fail', ...
            ['The problem description is not explicit enough to ' ...
             'compute a line-search hint.']);
        throw(up);

    end

end

author: andrewssobral

repo: lrslibrary


text:
function view = writeAnalyzeInTSeriesForm(view,scanList);
%
% view = writeAnalyzeInTSeriesForm(view,[scanList])
%
% Could be for either INPLANE or VOLUME
% Transform the analyze format images into time series format.
% By definition, each analyze image corresponds to one time
% frame and is a 3D matrix.
%
% scanList:
%   0 - do all scans
%   number or list of numbers - do only those scans
%   default - prompt user via selectScans dialog
%
% The user will identify, by GUI, where the source img files are, and where
% the tseries will be saved to.
%
% saveDir: if not given, then
%   ask a user to select a DataType or create a new DataType
%   save to 'Inplane/DataTypeName/TSeries_fromImgFiles/Scan*/TSeries*.mat'
%
% SPM Matlab toolbox required.
%
% MA, 10/18/2004. Based on ana2ts.m
% JL, 11/8/2004. Reorganize a little bit

mrGlobals;

% Choose the img file source
sourceChoice = questdlg('Where are the Analyze format Files in ?','choose source file','Current DataType','Other DataType','I Choose!','Current DataType');
switch sourceChoice;
    case 'Current DataType';
        curDataDir = dataDir(view);
    case 'Other DataType';
        %%%% Pop out menu to let user choose datatype.
        dtCount = size(dataTYPES,2);
        for i=1:dtCount; dtList{i} = dataTYPES(i).name; end;
        [selDataType,ok] = listdlg('ListString',dtList, 'SelectionMode','single', ...
            'InitialValue',curDataType, 'Name','Please, select a DataType:','ListSize', [300,200]);
        if ok;
            curDataDir = fullfile(viewGet(view,'subdir'),dtList{selDataType});
        end
    otherwise
        curDataDir = pwd;
end
% If no dir found, ask user to identify.
sourceDir = fullfile(curDataDir,'TSeries_imgFiles');
if ~exist(sourceDir,'dir') | strcmp(curDataDir,pwd);
    if ~strcmp(curDataDir,pwd); % which means the user did not choose 'I Choose!'
        disp('Warning: cannot find TSeries_imgFiles/ in the dataType.');
    end
    sourceDir = 0;
    while sourceDir == 0;     % Must choose one
        sourceDir = uigetdir(curDataDir,'Please choose the folder containing Scan* subfolders that have Analyze format files');
    end
end

% In the sourceDir, there must be folders Scan1, Scan2, ... Scan*
% (Re-)set scanList
if ~exist('scanList','var'); % User choose #scans
    scanList = selectScans(view,'Choose scans with IMG-files to Xform',sourceDir);
elseif scanList == 0 % all scans available
    [nScans, fileNames] = countDirs('Scan*',sourceDir);
    scanList = 1:nScans;
elseif ~isnumeric(scanList) | isempty(scanList)
    myErrorDlg('Xform aborted');
end

% Now, choose target save location
targetChoice = questdlg('Where to save the Xform tSeries?','choose save location','Current DataType','I Choose!','I Choose!');
switch targetChoice;
    case 'Current DataType';
        targetDir = tSeriesDir(view);
        disp('Warning: The old tSeries of this datatype will be overwritten.');
    otherwise
        targetDir = 0;
        while targetDir == 0;     % Must choose one
            targetDir = uigetdir(pwd,'Please choose the folder where Xform-ed time series is saved to');
        end
        disp('Warning: The Xform-ed time series may not match the scan params of current dataType.');
end

% current DataType params:
curDataType = viewGet(view,'currentDataType');
dtName = dataTYPES(curDataType).name;
scanParams = dataTYPES(curDataType).scanParams;

%select analyze files:
nScans = length(scanList);
for scanIndex=1:nScans
    scanNum = scanList(scanIndex);
    scanDir = ['Scan',int2str(scanNum)];
	sourcePath = fullfile(sourceDir,scanDir);
	targetPath = fullfile(targetDir,scanDir);

    %get files:
    files = dir(fullfile(sourcePath,'*.img'));
    for iList = 1:size(files);
        fName = files(iList).name;
        pat = '^\d+.img'; % must start with number and end with .img
        [s f tag] = regexp(lower(fName),pat);
        if ~isempty(s);
            seqNo = str2num(fName(1:end-4));
            F(seqNo+1,:) = fullfile(sourcePath,files(iList).name);
        end
    end
%    dpos=findstr(F(1,:),FILESEP);
%	sourcePath=F(1,1:dpos(end));

    if ~exist('F','var');
        myErrorDlg(['writeAnalyzeInTSeriesForm did not find any ', sourcePath, ' *.img files']);
    end;

    if scanIndex == 1
        msg = ['Saving TSeries Files for ',scanDir,' in ',targetPath];
        disp(msg);
        waitHandle = mrvWaitbar(0,[msg, '. Please wait...']);
    end

	%load volumes
	V = spm_vol(F);

    %check sizes. If mismatch, if user chose current datatype to save to, then error
    %and quit, but if user chose any location to save to, then just warning.
    errorString = {[]};
    if scanNum > length(scanParams);
        errorString{1} = ['Scan # exceeds the total number of scans in current dataType.']
    end
    params = scanParams(scanNum);
    if size(V) ~= params.nFrames;
        errorString{2} = ['Number of img-files does not match to Number of Frames for current dataType.'];
    end
    if V(1).dim(1:2) ~= params.cropSize;
        errorString{3} = ['cropSize in img-files does not match for current dataType.'];
    end
    if V(1).dim(3) ~= size(params.slices,2)
        errorString{4} = ['Number of Slices in img-files does not match for current dataType.'];
    end
    for ierror = 1:length(errorString);
        if ~isempty(errorString{ierror});
            if strcmp(targetChoice, 'Current DataType');
                myErrorDlg(errorString{ierror});
            else
                disp(['Warning!: ',errorString{ierror}]);
            end
        end
    end

    if ~exist(targetPath,'dir'); mkdir(targetDir,scanDir); end;

    %go through slices
    nSlices = V(1).dim(3);
	for slice=1:nSlices
        %get slice signal
        tSeries = [];
        C=[1 0 0 0;0 1 0 0;0 0 1 -slice;0 0 0 1];
        for it=1:size(V,1)
            temp = squeeze( spm_slice_vol(V(it),inv(C),V(it).dim(1:2),0) );
            tSeries(it,:) = temp(:)';
        end
        fileName = ['tSeries',num2str(slice),'.mat'];
        fprintf('Save slice %d in %s\n', slice, fullfile(targetPath,fileName));
        save(fullfile(targetPath,fileName),'tSeries');
        msg = [scanDir, ': saved slice in ', targetPath];
        mrvWaitbar((scanIndex-1+slice/nSlices)/nScans, waitHandle, msg);
	end
end

disp(['TSeries Files saved to ',targetDir]);
close(waitHandle);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [VIn,MIn] = partition_distance(Cx,Cy)
%PARTITION_DISTANCE     Distance between community partitions
%
%   This function quantifies the distance between pairs of community
%   partitions with information theoretic measures.
%
%   VIn       = partition_distance(Cx,Cy)
%   [VIn MIn] = partition_distance(Cx,Cy)
%
%   Inputs:         Cx,         Community affiliation vector X
%                   Cy,         Community affiliation vector Y
%
%   Outputs:        VIn,        Normalized variation of information
%                   MIn,        Normalized mutual information
%
%   (Definitions:
%       VIn = [H(X) + H(Y) - 2MI(X,Y)]/log(n)
%       MIn = 2MI(X,Y)/[H(X)+H(Y)]
%    where H is entropy, MI is mutual information and n is number of nodes)
%
%   Reference: Meila M (2007) J Multivar Anal 98, 873-895.
%
%   2011, Mika Rubinov, UNSW

%   Modification History:
%   Mar 2011: Original

%#ok<*ASGLU>

n = numel(Cx);                                  %n

[dum,dum,Cx]  = unique(Cx);
[dum,dum,Cy]  = unique(Cy);
[dum,dum,Cxy] = unique([Cx(:) Cy(:)],'rows');

Px  = hist(Cx,1:max(Cx))/n;                     %P(x)
Py  = hist(Cy,1:max(Cy))/n;                     %P(y)
Pxy = hist(Cxy,1:max(Cxy))/n;                   %P(x,y)

Hx  = -sum(Px.*log(Px));                        %H(x)
Hy  = -sum(Py.*log(Py));                        %H(y)
Hxy = -sum(Pxy.*log(Pxy));                      %H(x,y)

VIn = (2*Hxy-Hx-Hy)/log(n);                     %VIn
MIn = 2*(Hx+Hy-Hxy)/(Hx+Hy);                    %MIn


text:
function Q = vbpcamv(Y, d, varargin)
% Modified RPPCA by Jaakko Luttinen
% Uses separate outlier models for each observed signal.
% p(y|W,x,mu,u,tau) = N(y | W*x+mu, diag(1/u*1/tau))
% Uses hierarchical models and variational learning.
%
% If you have troubles with bad minima, you may want to try changing the
% time when the algorithm starts to update the hyperparameters
% (startupdatehyper) and/or to rotate (startrotate).

% Last updated 13th February 2009, Jaakko Luttinen

opts = struct( ...
    'init',          [],...
    'prior', [], ...
    'hiermean', false, ...
    'rotate', true, ...
    'commonnoise', true, ...
    'startrotate', 1, ...
    'startupdatehyper', 1, ...
    'autosavetime', 0,...
    'autosavefile', 'vbpcamv_autosave',...
    'testset', [], ...
    'fixw', false, ...
    'maxiters', 1000, ...
    'convergence', eps);

[ opts, errmsg, wrnmsg ] = argschk( opts, varargin{:} );
if ~isempty(errmsg), error( errmsg ), end
if ~isempty(wrnmsg), warning( wrnmsg ), end

[m,n] = size(Y);

% Missing values
if issparse(Y)
  Obs = Y~=0; %spones(Y);
else
  Obs = ~isnan(Y);
end
[ObsI, ObsJ] = find(Obs);

%%%%%%%%%%%%%%%%%%%
%% INITIALIZATION

nmv = sum(Obs,2);

% NAMING CONVENTION:
% variables are named without '_': e.g. a->tau->mu ==> ataumu
% but the values of posterior distribution parameters are named with one
% '_': e.g. a->tau->mu ==> a_taumu
% So, priors (or hyperparameter variables) have no '_'. Understand the
% difference!
% But for clarity mu_W is just W and mu_X is just X and so on. That is,
% plain variable name (W,X,mu,tau,logtau,U,logU) refers to posterior
% expectation if there is no naming conflict.

% Mean (mu)
% $$$ mumumu = 0;     % prior for hierarchical stuff
% $$$ taumumu = 1e-5; % prior for hierarchical stuff
v_mumu = 0;
ataumu = 1e-3;  % prior for hierarchical stuff
btaumu = 1e-3;  % prior for hierarchical stuff
mumu = 0;       % init posterior or prior if no hierarchical mean
taumu = 1e-6;   % init posterior or prior if no hierarchical mean
v_mu = 1e0*ones(m,1); % init posterior
mu = zeros(m,1);      % init posterior

% Loading matrix (W)
%aw = 1e-16; % prior
%bw = 1e-13; % prior
aw = 1e-10; % prior
bw = 1e-5; % prior
a_w = aw * ones(1,d); % init size
b_w = bw * ones(1,d); % init size
w = 1e-6 * ones(1,d);    % init posterior
logw = nan * ones(1,d);  % init size
W = orth(randn(m,d));          % init posterior
CovW = 1e-0*repmat(eye(d),[1,1,m]); % init posterior
varW = nan*ones(m,d);              % init size

% Observation noise
atau = 1e-4; % prior
btau = 1e-4; % prior
a_tau = nan * ones(m,1);
b_tau = nan * ones(m,1);
tau = 1e2 * ones(m,1); % init posterior (THIS SEEMS TO BE IMPORTANT!!!)
%tau = 1./varmv(Y,2) * (m^2)/((m-d)^2)   % adhoc init posterior
logtau = nan * ones(m,1); % init size

% Outlier parameters
% $$$ nu = 10 * ones(m,1); % variable init
% $$$ a_U = nan * ones(m,n); % init size
% $$$ b_U = nan * ones(m,n); % init size
U = ones(m,n);     % init posterior
% $$$ logU = nan * zeros(m,n); % init size

% Initialize sizes
X = nan * zeros(d,n); % init size
Sv = cell(n,1);
Sv(:) = {nan*zeros(d,d)};
CovX = nan*zeros(d,d,n);
XX = nan*zeros(d,d,n);
UXX = nan*zeros(d,d,m);

% Prior information
if isstruct(opts.prior)
  if isfield(opts.prior, 'mumumu')
    mumumu(1) = opts.prior.mumumu
  end
  if isfield(opts.prior, 'taumumu')
    taumumu(1) = opts.prior.taumumu
  end
  if isfield(opts.prior, 'ataumu')
    ataumu(1) = opts.prior.ataumu
  end
  if isfield(opts.prior, 'btaumu')
    btaumu(1) = opts.prior.btaumu
  end
  if isfield(opts.prior, 'aw')
    aw(1) = opts.prior.aw
  end
  if isfield(opts.prior, 'bw')
    bw(1) = opts.prior.bw
  end
  if isfield(opts.prior, 'atau')
    atau(1) = opts.prior.atau
  end
  if isfield(opts.prior, 'btau')
    btau(1) = opts.prior.btau
  end
end

% Use given initialization (size matching is checked by using ':' )
if isstruct(opts.init)
  if isfield(opts.init, 'W')
    W(:,:) = opts.init.W;
    fprintf('Using old W.\n');
  end
  if isfield(opts.init, 'CovW')
    CovW(:,:,:) = opts.init.CovW;
    fprintf('Using old CovW.\n');
  end
  if isfield(opts.init, 'w')
    w(:) = opts.init.w
    fprintf('Using old w.\n');
  end
  if isfield(opts.init, 'X')
    X(:,:) = opts.init.X;
    fprintf('Using old X.\n');
  end
  if isfield(opts.init, 'CovX')
    CovX(:,:,:) = opts.init.CovX;
    fprintf('Using old CovX.\n');
  end
  if isfield(opts.init, 'Sv')
    Sv(:) = opts.init.Sv;
    fprintf('Using old Sv.\n');
  end
  if isfield(opts.init, 'mu')
    mu(:) = opts.init.mu;
    fprintf('Using old mu.\n');
  end
  if isfield(opts.init, 'v_mu')
    v_mu(:) = opts.init.v_mu;
    fprintf('Using old v_mu.\n');
  end
  if isfield(opts.init, 'mumu')
    mumu(1) = opts.init.mumu;
    fprintf('Using old mumu.\n');
  end
  if isfield(opts.init, 'taumu')
    taumu(1) = opts.init.taumu;
    fprintf('Using old taumu.\n');
  end
  if isfield(opts.init, 'tau')
    tau(:) = opts.init.tau;
    fprintf('Using old tau.\n');
  end
  if isfield(opts.init, 'nu')
    nu(:) = opts.init.nu;
    fprintf('Using old nu.\n');
  end
  if isfield(opts.init, 'U')
    U(:,:) = opts.init.U;
    fprintf('Using old U.\n');
  end
end

WW = zeros(d,d,m);
for i=1:m
  WW(:,:,i) = W(i,:)'*W(i,:) + CovW(:,:,i);
end

% Log-likelihood
logP = -Inf;

Im = eye(m);
Id = eye(d);

lastsave = now;

N = 1:n;
vM = 1:m;
m_mv = sum(Obs,1);
n_mv = sum(Obs,2);
nm_mv = sum(Obs(:)); % number of observations

log2pi = log(2*pi);

oldcosts = -inf;

% Monitor the overall time of the iteration
starttime = cputime;
duration = nan*zeros(opts.maxiters,1);
rmse = nan*zeros(opts.maxiters,1);
rmse_test = nan*zeros(opts.maxiters,1);

cost = nan*zeros(opts.maxiters,1);
for k=1:opts.maxiters

  itertime = cputime;

  % This is used to monitor the rotation of the subspace after each step
  W_old = W;

  % Update X
  for j=1:n
    imv = Obs(:,j);
    UTau = diag(U(imv,j).*tau(imv));
    UWWTau = 0;
    for i=vM(imv)
      UWWTau = UWWTau + U(i,j)*tau(i) * WW(:,:,i);
    end
    CovX(:,:,j) = inv(UWWTau + Id);
    X(:,j) = CovX(:,:,j) * W(imv,:)' * UTau * (Y(imv,j)-mu(imv));
    XX(:,:,j) = CovX(:,:,j) + X(:,j)*X(:,j)';
    Sv{j} = CovX(:,:,j); % only for backward compatibility..
  end

  for i=1:m
    jmv = Obs(i,:);

    % Update sum_j(<u><xx>)
    UXX(:,:,i) = 0;
    for j=N(jmv) % sum over observed time instances
      UXX(:,:,i) = UXX(:,:,i) + U(i,j)*XX(:,:,j);%(CovX(:,:,j) + X(:,j)*X(:,j)');
    end

    % Update W
    CovW(:,:,i) = inv(UXX(:,:,i)*tau(i) + diag(w));
    W(i,:) = (U(i,jmv).*(Y(i,jmv)-mu(i))*X(:,jmv)') * CovW(:,:,i) * tau(i);
    WW(:,:,i) = W(i,:)'*W(i,:) + CovW(:,:,i);

    % Update mu
    v_mu(i) = 1 / (taumu + sum(U(i,jmv))*tau(i));
    mu(i) = v_mu(i) * (taumu*mumu + U(i,jmv)*(Y(i,jmv)-W(i,:)*X(:,jmv))'*tau(i));

  end

  % Rotate
  % THINK/TODO: WHAT WOULD BE A GOOD TIME TO START ROTATING??
  % (rotating too early may lead to bad local minima and not rotating
  % makes converging very slow) should one measure how "converged" the
  % loglikelihood is and then decide whether to start rotate?
  if opts.rotate && k>=opts.startrotate
% $$$     before = trace( sum(CovX,3) * sum(CovW,3) )
    [W,CovW,X,Sv,CovX,mu] = orthogonalize(W,CovW,X,Sv,CovX,mu,opts.fixw,tau,Obs);
% $$$     after = trace( sum(CovX,3) * sum(CovW,3) )
%    [W,CovW,X,Sv,CovX,mu] = orthogonalize_TEST(W,CovW,X,Sv,CovX,mu,tau(1)); false
%    [W,CovW,X,CovX,mu] = RotateToPCA(W,CovW,X,CovX,mu);
  end

  for i=1:m
    jmv = Obs(i,:); % observed

    WW(:,:,i) = W(i,:)'*W(i,:) + CovW(:,:,i);
    varW(i,:) = diag(CovW(:,:,i));

    % Update sum_j(<u><xx>)
    UXX(:,:,i) = 0;
    for j=N(jmv) % sum over observed time instances
      UXX(:,:,i) = UXX(:,:,i) + U(i,j)*(CovX(:,:,j) + X(:,j)*X(:,j)');
    end

    % Update tau
    UWXXW = sum(sum( WW(:,:,i) .* UXX(:,:,i) ));
    a_tau(i) = atau + 0.5*sum(jmv);
    b_tau(i) = btau ...
        + 0.5 * (UWXXW + ...
                 sum(-2*W(i,:)*X(:,jmv)*(U(i,jmv).*(Y(i,jmv)-mu(i)))') + ...
                 sum(U(i,jmv).*(Y(i,jmv)-mu(i)).^2) + ...
                 sum(U(i,jmv)*v_mu(i)));

  end

  % "AD HOC"(?): Use common variance
  if opts.commonnoise
    a_tau(:) = sum(a_tau-atau) + atau;
    b_tau(:) = sum(b_tau-btau) + btau;
  end
  tau = a_tau./b_tau;
  logtau = psi(a_tau) - log(b_tau);
  % TODO: DEBUGGING: USE ML ESTIMATE
% $$$   logtau = log(tau);

  % Update hyperparameters for mu
  if ~opts.hiermean
    logtaumu = log(taumu);
  else
    a_taumu = ataumu + 0.5*m;
    b_taumu = btaumu + 0.5*sum( (mu-mumu).^2 + v_mu + v_mumu);
    taumu = a_taumu ./ b_taumu; % posterior mean of accuracy hyperparameter
    logtaumu = psi(a_taumu) - log(b_taumu);
  end

  % Update hyperparameter for W (you can try to start updating these
  % hyperparameters only after some iteration if some problems seem to
  % appear.)
  if ~opts.fixw
    if k >= opts.startupdatehyper
      a_w(:) = aw + 0.5*m;
      b_w(:) = bw + 0.5*sum(W.^2 + varW, 1);
    end
    w = a_w ./ b_w; % posterior mean of precision hyperparameter
    logw = psi(a_w) - log(b_w);
  else
    logw = log(w);
  end

  % Change of the principal subspace (in radians)
  angle = subspace(W,W_old);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Calculate lower bound of the log-likelihood

  % Cost from Y
  cost_Y = 0;
  for j=1:n
    imv = Obs(:,j);
    for i=vM(imv)
      mean_muY = W(i,:)*X(:,j)+mu(i);
      var_muY = W(i,:)*CovX(:,:,j)*W(i,:)' + X(:,j)'*CovW(:,:,i)*X(:,j) ...
                + sum(sum( CovX(:,:,j).*CovW(:,:,i) )) + v_mu(i);
      tauY = tau(i);
      logtauY = logtau(i);
      cost_Y = cost_Y + cost_gaussian(Y(i,j),[], mean_muY,var_muY, ...
                                      tauY,logtauY);
    end
  end

  % Cost from W
  cost_W = 0;
  for i=1:m
    cost_W = cost_W + cost_gaussian(W(i,:)',CovW(:,:,i), 0,0, w',logw');
  end

  % Cost from X
  cost_X = 0;
  for j=1:n
    cost_X = cost_X + cost_gaussian(X(:,j),CovX(:,:,j), 0,0, 1,0);
  end

  % Cost from mu
  cost_mu = cost_gaussian(mu,diag(v_mu), mumu,v_mumu, taumu,logtaumu);

  % Cost from tau
  if opts.commonnoise
    cost_tau = cost_gamma(tau(1),logtau(1),a_tau(1),b_tau(1),atau,btau);
  else
    cost_tau = sum(cost_gamma(tau,logtau,a_tau,b_tau,atau,btau));
  end

  % Cost from w
  if ~opts.fixw
    cost_w = sum(cost_gamma(w,logw,a_w,b_w,aw,bw));
  else
    cost_w = 0;
  end

  if opts.hiermean
    cost_taumu = cost_gamma(taumu,logtaumu,a_taumu,b_taumu,ataumu,btaumu);
  end

  oldlogP = logP;

  % Lower bound of loglikelihood
  logP = cost_Y + cost_W + cost_X + cost_mu + cost_tau + cost_w;
  %logP = cost_Y;
  costs = [cost_Y; cost_W; cost_X; cost_mu; cost_tau; cost_w];

  cost(k) = logP;

  % Show progress
  time = cputime - itertime;
  if k > 1
    duration(k) = duration(k-1) + time;
  else
    duration(1) = time;
  end
  fprintf('Step %d: loglike=%e, angle=%e  (%f seconds)\n', ...
          k, logP, angle, time);

  % Debugging..
  if logP < oldlogP
%    diff_costs = costs - oldcosts
%    error('Log-likelihood bound not improved! Bug in code or likelihood?');
    warning('Log-likelihood bound decreased! Bug in code or likelihood?');
  end
  oldcosts = costs;

  % RMSE
  Yh = W*X + repmat(mu, 1,n);
  err = Y-Yh;
  err = err(~isnan(err));
  rmse(k) = sqrt( mean(err(:).^2) );
  if ~isempty(opts.testset)
    err = opts.testset-Yh;
    err = err(~isnan(err));
    rmse_test(k) = sqrt( mean(err(:).^2) );
  end

  % Check whether to save the results
  if (now-lastsave)*3600*24 >= opts.autosavetime && opts.autosavetime > 0
    fprintf('Saving to %s...', opts.autosavefile);
    loglikelihood = logP;
    if opts.hiermean
      save(opts.autosavefile, 'W','CovW','w','a_w','b_w','X','CovX','Sv', ...
           'mu','v_mu','mumu','v_mumu','taumu','a_taumu','b_taumu','tau', ...
           'a_tau','b_tau','U','a_U','b_U','nu','loglikelihood');
    else
      save(opts.autosavefile, 'W','CovW','w','a_w','b_w','X','CovX','Sv', ...
           'mu','v_mu','tau', 'a_tau','b_tau', 'loglikelihood', 'cost', ...
           'duration', 'rmse', 'rmse_test');
    end
    lastsave = now;
    fprintf(' done.\n');
  end

% $$$   if angle < 1e-14
% $$$     fprintf('Stopping iteration: only minor change in subspace\n');
% $$$     break
% $$$   end
  improvement = abs((logP - oldlogP) / logP);
  if improvement < opts.convergence && k > opts.startupdatehyper
    fprintf('Stopping iteration: only minor improvement in loglikelihood\n');
    break
  end

end

fprintf('Saving to %s...', opts.autosavefile);
loglikelihood = logP;
if opts.hiermean
  save(opts.autosavefile, 'W','CovW','w','a_w','b_w','X','CovX','Sv', ...
       'mu','v_mu','mumu','v_mumu','taumu','a_taumu','b_taumu','tau', ...
       'a_tau','b_tau','U','a_U','b_U','nu','loglikelihood');
else
  save(opts.autosavefile, 'W','CovW','w','a_w','b_w','X','CovX','Sv', ...
       'mu','v_mu','tau', 'a_tau','b_tau', 'loglikelihood', 'cost', ...
       'duration', 'rmse', 'rmse_test');
end
fprintf(' done.\n');

% Results as a struct if desired
if true %nargout == 1
  results.W = W;
  results.CovW = CovW;
  results.w = w;
  results.a_w = a_w;
  results.b_w = b_w;
  results.X = X;
  results.CovX = CovX;
  results.Sv = Sv;
  results.mu = mu;
  results.v_mu = v_mu;
  if opts.hiermean
    results.mumu = mumu;
    results.v_mumu = v_mumu;
    results.taumu = taumu;
    results.a_taumu = a_taumu;
    results.b_taumu = b_taumu;
  end
  results.tau = tau;
  results.a_tau = a_tau;
  results.b_tau = b_tau;
% $$$   results.U = U;
% $$$   results.a_U = a_U;
% $$$   results.b_U = b_U;
% $$$   results.nu = nu;
  results.loglikelihood = logP;
  results.cost = cost;
  results.time = duration;
  Q = results;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c = cost_gaussian(x,Covx,mu,Covmu,tau,logtau)
% c = cost_gaussian(x,Covx,mu,Covmu,tau,logtau)
%
% Calculates the difference between <log p(X)> - <log q(X)>
% where expectation is over q(X).
%
% A lower bound for a Gaussian:
% p(X) = N(MU,1/TAU)
% q(X) = N(x,Covx)
%
% Rest of the parameters are defined as:
% q(MU) = N(mu,Covmu)
% <TAU> = tau
% <log TAU> = logtau
%
% If Covx==[], then the term <log q(X)> is not calculated.
% This is useful when X is, e.g., observations.

c = 0;

% Cost from q-posterior
if ~isempty(Covx)
  entropy = entropy_gaussian(Covx);
  c = entropy;
else
  Covx = 0;
end

% Cost from prior
err2 = ((x-mu).^2 + diag(Covx) + diag(Covmu));
c = c + 0.5 * sum( -log(2*pi) + logtau - tau .* err2 );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c = cost_gamma(x,logx,apost,bpost,aprior,bprior)
% A lower bound for a Gamma distributed X:
% p(X) = G(aprior,bprior)
% q(X) = G(apost,bpost)
% <X> = x
% <log X> = logx

% Cost from prior
c = aprior*log(bprior) - gammaln(aprior) + (aprior-1)*logx - bprior*x;
% Cost from q-posterior
c = c + entropy_gamma(apost,bpost);

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

function H = entropy_gaussian(Cov)
% Cov is covariance matrix
d = size(Cov,1);

log2pi = log(2*pi);
%H = 0;

%for j=1:size(Cov,3)
L = chol(Cov);
H = d/2*log2pi + 0.5*2*logdettri(L) + d/2;
%end

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

function H = entropy_gamma(a,b)
% a is shape, b is inverse scale

H = a - log(b) + gammaln(a) - (a-1).*psi(a);
%H = sum(H(:));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ A, Av, S, Sv, mu ] = ...
    RotateToPCA( A, Av, S, Sv, mu);

n1 = size(A,1);
n2 = size(S,2);

% TODO: Take into account the prior for Mu, A???
mS = mean(S,2);
dMu = A*mS;
S = S - repmat(mS,1,n2);
mu = mu + dMu;

covS = S*S';
for j = 1:n2
  covS = covS + Sv(:,:,j);
end

R = 1;

% w.r.t. S
covS = covS / n2;
[VS,D] = eig(covS);
RA = VS*sqrt(D);
A = A*RA;
covA = A'*A;
for i = 1:n1
  Av(:,:,i) = RA'*Av(:,:,i)*RA;
  covA = covA + Av(:,:,i);
end
R = diag(1./sqrt(diag(D)))*VS';

% $$$ % w.r.t. A
% $$$ covA = covA / n1;
% $$$ [VA,DA] = eig(covA);
% $$$ [DA,I] = sort( -diag(DA) );
% $$$ DA = -DA;
% $$$ VA = VA(:,I);
% $$$ A = A*VA;
% $$$ for i = 1:n1
% $$$   Av(:,:,i) = VA'*Av(:,:,i)*VA;
% $$$ end
% $$$ R = VA'*R;

S = R*S;
for j = 1:length(Sv)
    Sv(:,:,j) = R*Sv(:,:,j)*R';
end

return

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

function [W,CovW,X,Sv,CovX,mu] = orthogonalize(W,CovW,X,Sv,CovX,mu,fixw,tau,Obs)
[m,c] = size(W);
n = size(X,2);

WW = W'*W + sum(CovW,3);

% $$$ old_recon = W*X + repmat(mu,1,n);

% Use ML zero mean (mu should be updated after this function!!) ..
if false
  disp('Using complicated but exact mean translation');
  tau = reshape(tau,[1,1,m]);
  nom = 0;
  denom = 0;
  I = eye(c);
  for j=1:n
    imv = Obs(:,j);
    tmp = I + sum(bsxfun(@times,tau(1,1,imv),CovW(:,:,imv)),3);
    nom = nom + tmp * X(:,j);
    denom = denom + tmp;
  end
  dmu = denom \ nom;
else
  dmu = mean(X,2);
end
% .. or my analytic zero mean with proper fix to mu:
% $$$  dmu = inv(n*eye(c)+taumu*WW) * (sum(X,2) + taumu*W'*(mumu-mu));

% Move bias
X = X - repmat(dmu,1,n);
mu = mu + W*dmu;
% $$$ nobiasmove = true

Qx = 1;

% Whiten w.r.t. X
XX = X*X' + sum(CovX,3);
if fixw
  [Vx,Dx,tmp] = svd(XX/(n-m)); % USE THIS IF FIXED w ??
else
  [Vx,Dx,tmp] = svd(XX/n);
end
Qx = diag(1./sqrt(diag(Dx))) * Vx';
Qw = Vx*sqrt(Dx);
W = W * Qw;
for i=1:size(CovW,3)
  CovW(:,:,i) = Qw'*CovW(:,:,i)*Qw;
end
X = Qx * X;
for j = 1:size(CovX,3)
  Sv{j} = Qx*Sv{j}*Qx';
  CovX(:,:,j) = Qx*CovX(:,:,j)*Qx';
end

% Check that XX is really whitened! (because of numerical issues!!)
XX = X*X' + sum(CovX,3);
if fixw
  [Vx,Dx,tmp] = svd(XX/(n-m)); % USE THIS IF FIXED w ??
else
  [Vx,Dx,tmp] = svd(XX/n);
end
if Dx(1) > 1.1
  warning('Needs to whiten X again. See vbpcamv->orthogonalize');
  % Whiten w.r.t. X AGAIN
  Qx = diag(1./sqrt(diag(Dx))) * Vx';
  Qw = Vx*sqrt(Dx);
  W = W * Qw;
  for i=1:size(CovW,3)
    CovW(:,:,i) = Qw'*CovW(:,:,i)*Qw;
  end
  X = Qx * X;
  for j = 1:size(CovX,3)
    Sv{j} = Qx*Sv{j}*Qx';
    CovX(:,:,j) = Qx*CovX(:,:,j)*Qx';
  end
end

Qx = 1;
% Diagonalize w.r.t. W
WW = W'*W + sum(CovW,3);
[Vw,Dw,tmp] = svd(WW);
%[Dw,I] = sort(diag(Dw), 'descend');
%Vw = Vw(:,I);
Qx = Vw' * Qx;
Qw = Vw;
W = W * Qw;
for i=1:size(CovW,3)
  CovW(:,:,i) = Qw'*CovW(:,:,i)*Qw;
end
X = Qx * X;
for j = 1:size(CovX,3)
  Sv{j} = Qx*Sv{j}*Qx';
  CovX(:,:,j) = Qx*CovX(:,:,j)*Qx';
end

% $$$ xx = X*X' + sum(CovX,3)
% $$$ ww = W'*W + sum(CovW,3)

% $$$ new_recon = W*X + repmat(mu,1,n);
% $$$
% $$$ mse = mean((old_recon(:) - new_recon(:)).^2)

return
% $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ function [W,CovW,X,Sv,CovX,mu] = orthogonalize_TEST(W,CovW,X,Sv,CovX,mu,tau)
% $$$ [m,c] = size(W);
% $$$ n = size(X,2);
% $$$
% $$$ WW = W'*W + sum(CovW,3);
% $$$
% $$$ old_recon = W*X + repmat(mu,1,n);
% $$$
% $$$ % Use ML zero mean (mu should be updated after this function!!) ..
% $$$ % $$$ dmu = mean(X,2);
% $$$ % $$$ tau*sum(CovW,3)
% $$$ dmu = inv(n*eye(c)+tau*n*sum(CovW,3)) * (sum((eye(c)+tau*sum(CovW,3))*X, 2));
% $$$ %(dmu2-dmu)./dmu
% $$$ % .. or my analytic zero mean with proper fix to mu:
% $$$ % $$$  dmu = inv(n*eye(c)+taumu*WW) * (sum(X,2) + taumu*W'*(mumu-mu));
% $$$
% $$$ % Move bias
% $$$ X = X - repmat(dmu,1,n);
% $$$ mu = mu + W*dmu;
% $$$
% $$$ Qx = 1;
% $$$
% $$$ % Whiten w.r.t. X
% $$$ XX = X*X' + sum(CovX,3);
% $$$ [Vx,Dx] = eig(XX/n);
% $$$ %[Vx,Dx] = eig(XX/(n+m));
% $$$ %[Vx,Dx] = eig(XX/(n-m));
% $$$ Qx = diag(1./sqrt(diag(Dx))) * Vx';
% $$$ Qw = Vx*sqrt(Dx);
% $$$ W = W * Qw;
% $$$ for i=1:size(CovW,3)
% $$$   CovW(:,:,i) = Qw'*CovW(:,:,i)*Qw;
% $$$ end
% $$$
% $$$ % Whiten w.r.t. W
% $$$ WW = W'*W + sum(CovW,3);
% $$$ [Vw,Dw] = eig(WW);
% $$$ [Dw,I] = sort(diag(Dw), 'descend');
% $$$ Vw = Vw(:,I);
% $$$ Qx = Vw' * Qx;
% $$$ Qw = Vw;
% $$$ W = W * Qw;
% $$$ for i=1:size(CovW,3)
% $$$   CovW(:,:,i) = Qw'*CovW(:,:,i)*Qw;
% $$$ end
% $$$
% $$$ X = Qx * X;
% $$$ for j = 1:size(CovX,3)
% $$$   Sv{j} = Qx*Sv{j}*Qx';
% $$$   CovX(:,:,j) = Qx*CovX(:,:,j)*Qx';
% $$$ end
% $$$
% $$$ new_recon = W*X + repmat(mu,1,n);
% $$$
% $$$ mse = mean((old_recon(:) - new_recon(:)).^2);
% $$$
% $$$ return
% $$$

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
% Spectral factorization using Kolmogorov 1939 approach
% (code follows pp. 232-233, Signal Analysis, by A. Papoulis)
%
% Computes the minimum-phase impulse response which satisfies
% given auto-correlation.
%
% Input:
%   r: top-half of the auto-correlation coefficients
%      starts from 0th element to end of the auto-corelation
%      should be passed in as a column vector
% Output
%   h: impulse response that gives the desired auto-correlation

function h = spectral_fact(r)

% length of the impulse response sequence
nr = length(r);
n  = (nr+1)/2;

% over-sampling factor
mult_factor = 30;        % should have mult_factor*(n) >> n
m = mult_factor*n;

% computation method:
% H(exp(jTw)) = alpha(w) + j*phi(w)
% where alpha(w) = 1/2*ln(R(w)) and phi(w) = Hilbert_trans(alpha(w))

% compute 1/2*ln(R(w))
w = 2*pi*[0:m-1]/m;
R = exp( -j*kron(w',[-(n-1):n-1]) )*r;
R = abs(real(R)); % remove numerical noise from the imaginary part
figure; plot(20*log10(R));
alpha = 1/2*log(R);

% find the Hilbert transform
alphatmp = fft(alpha);
alphatmp(floor(m/2)+1:m) = -alphatmp(floor(m/2)+1:m);
alphatmp(1) = 0;
alphatmp(floor(m/2)+1) = 0;
phi = real(ifft(j*alphatmp));

% now retrieve the original sampling
index  = find(rem([0:m-1],mult_factor)==0);
alpha1 = alpha(index);
phi1   = phi(index);

% compute the impulse response (inverse Fourier transform)
h = ifft(exp(alpha1+j*phi1),n);

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [nll,g,G] = LogisticLossGN(w,X,y,weights)
% Negative log likelihood for binary logistic regression
% This code returns the Gauss Newton approximation to the Hessian
% w: d*1
% X: n*d
% y: n*1, should be -1 or 1

% This file is from pmtk3.googlecode.com

[n, d] = size(X);
if nargin < 4, weights = ones(n,1); end
y01 = (y+1)/2;
mu = sigmoid(X*w);
mu = max(mu, eps); % bound away from 0
mu = min(1-eps, mu); % bound away from 1
nll = -sum(weights .* (y01 .* log(mu) + (1-y01) .* log(1-mu)));
Xw = X .* repmat(colvec(weights), 1, d);
if nargout > 1
  g = Xw'*(mu-y01);
end

if nargout > 2
  g = Xw.*(mu-y01);
  G = g'*g;
  %H = Xw'*diag(mu.*(1-mu))*X;
end

end

author: emtiyaz

repo: vadam

save_path: github-repos/MATLAB/emtiyaz-vadam


text:
function obsNames = prtUtilGenerateDefaultObservationNames(indices2)
% obsNames = prtUtilGenerateDefaultObservationNames(indices2)

% Previously static in prtDataSetBase
%     function obsNames = generateDefaultObservationNames(indices2)
%         obsNames = prtUtilCellPrintf('Observation %d',num2cell(indices2));
%         obsNames = obsNames(:);
%     end
obsNames = prtUtilCellPrintf('Observation %d',num2cell(indices2));
obsNames = obsNames(:);


text:
function v = std3(f)
%STD3   Standard deviation of a CHEBFUN3.
%   V = STD3(F) computes the standard deviation of the CHEBFUN3 object F,
%   i.e.,
%             STD3(F)^2 = 1 / V*sum3(|F(x,y,z) - M|^2)
%
%   where V is the volume of the domain of F and M is the mean of F.
%
% See also CHEBFUN3/MEAN, CHEBFUN3/MEAN2, CHEBFUN3/MEAN3, CHEBFUN3/STD and
% CHEBFUN3/STD2.

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

h = f - mean3(f);
v = sqrt(mean3(h .* conj(h)));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% Test sharedmatrix mex file.
%
% Copyright (c) 2010,2011 Joshua V Dillon
% All rights reserved. (See file header for details.)

%C = sparse(rand(10000,10000).*(rand(10000,10000)>.1)); % big matrix
C = sparse(rand(100,100).*(rand(100,100)>.1)); % small matrix

shmkey=12345;

try,sharedmatrix('free',shmkey);catch,end;

shmsiz=sharedmatrix('clone',shmkey,C);clear C;
%fprintf('[shmkey shmsiz]=deal(%d,%d);\n',shmkey,shmsiz);

if exist('x','var')==1
	sharedmatrix('detach',shmkey,x);
	clear x;
end

for i=1:100
	x=sharedmatrix('attach',shmkey);
	I=randperm(size(x,2));
	a=sum(x,1);
%	pause(.5);
	sharedmatrix('detach',shmkey,x);
%	clear x;
end


text:
function varargout = svmtrain(varargin)
% VL_SVMTRAIN   Train a Support Vector Machine
%   [W B] = VL_SVMTRAIN(X, Y, LAMBDA) trains a linear Support Vector
%   Machine (SVM) from the data vectors X and the labels Y. X is a D
%   by N matrix, with one column per example and D feature dimensions
%   (SINGLE or DOUBLE). Y is a DOUBLE vector with N elements with a
%   binary (-1 or +1) label for each training point. To a first order
%   approximation, the function computes a weight vector W and offset
%   B such that the score W'*X(:,i)+B has the same sign of LABELS(i)
%   for all i.
%
%   VL_SVMTRAIN(DATASET, LABELS, LAMBDA) takes as input a DATASET
%   structure, which allows more sophisticated input formats to be
%   supported (see VL_SVMDATASET()).
%
%   [W, B, INFO] = VL_SVMTRAIN(...) additionally returns a structure
%   INFO with the following fields:
%
%   iteration::
%     Number of iterations performed.
%
%   epoch::
%     Number of iterations over number of training data points.
%
%   elapsedTime::
%     Time elapsed since the start of training.
%
%   objective::
%     SVM objective value.
%
%   regularizer::
%     Regularizer value.
%
%   loss::
%     Loss value.
%
%   scoreVariation:: [SGD only]
%     Mean square root of the difference between the last two
%     values of the SVM scores for each point.
%
%   dualObjective:: [SDCA only]
%     Dual objective value.
%
%   dualLoss:: [SDCA only]
%     Dual loss value::
%
%   dualityGap:: [SDCA only]
%     Difference between the objective and the dual objective.
%
%   [W, B, INFO, SCORES] = VL_SVMTRAIN(X, Y, LABMDA) returns a row
%   vector of the SVM score for each training point. This can be used
%   in combination with the options SOLVER, MODEL, and BIAS to
%   evaluate an existing SVM on new data points. Furthermore INFO will
%   contain the corresponding SVM loss, regularizer, and objective
%   function value. If this information is not of interest, it is
%   possible to pass a null vector Y instead of the actual labels as
%   well as a null regularizer.
%
%   VL_SVMTRAIN() accepts the following options:
%
%   Verbose::
%     Specify one or multiple times to increase the verbosity level.
%     Given only once, produces messages at the beginning and end of
%     the learning. Verbosity of at least 2 prints information at
%     every diagnostic step.
%
%   Epsilon:: 1e-3
%     Tolerance for the stopping criterion.
%
%   MaxNumIterations:: 10/LAMBDA
%     Maximum number of iterations.
%
%   BiasMultiplier:: 1
%     Value of the constant B0 used as bias term (see below).
%
%   BiasLearningRate:: 0.5
%     Learning rate for the bias (SGD solver only).
%
%   DiagnosticFunction:: []
%     Diagnostic function callback. The callback takes the INFO
%     structure as only argument. To trace energies and plot graphs,
%     the callback can update a global variable or, preferably, be
%     defined as a nested function and update a local variable in the
%     parent function.
%
%   DiagnosticFrequency:: Number of data points
%     After how many iteration the diagnostic is run. This step check
%     for convergence, and is done rarely, typically after each epoch
%     (pass over the data). It also calls the DiangosticFunction,
%     if any is specified.
%
%   Loss:: HINGE
%     Loss function. One of HINGE, HINGE2, L1, L2, LOGISTIC.
%
%   Solver:: SDCA
%     One of SGD (stochastic gradient descent [1]), SDCA (stochastic
%     dual coordinate ascent [2,3]), or NONE (no training). The
%     last option can be used in combination with the options MODEL
%     and BIAS to evaluate an existing SVM.
%
%   Model:: null vector
%     Specifies the initial value for the weight vector W (SGD only).
%
%   Bias:: 0
%     Specifies the initial value of the bias term (SGD only).
%
%   Weights:: []
%     Specifies a weight vector to assign a different non-negative
%     weight to each data point. An application is to rebalance
%     unbalanced datasets.
%
%   FORMULATION
%
%   VL_SVMTRAIN() minimizes the objective function of the form:
%
%     LAMBDA/2 |W|^2 + 1/N SUM_i LOSS(W' X(:,i), Y(i))
%
%   where LOSS(W' Xi,Yi) is the loss (hinge by default) for i-th
%   data point. The bias is incorporated by extending each data
%   point X with a feature of constant value B0, such that the
%   objective becomes
%
%    LAMBDA/2 (|W|^2 + WB^2) 1/N SUM_i LOSS(W' X(:,i) + WB B0, Y(i))
%
%   Note that this causes the learned bias B = WB B0 to shrink
%   towards the origin.
%
%   Example::
%     Learn a linear SVM from data X and labels Y using 0.1
%     as regularization coefficient:
%
%       [w, b] = vl_svmtrain(x, y, 0.1) ;
%
%     The SVM can be evaluated on new data XTEST with:
%
%       scores = w'*xtest + b ;
%
%     Alternatively, VL_SVMTRAIN() can be used for evaluation too:
%
%       [~,~,~, scores] = vl_svmtrain(xtest, y, 0, 'model', w, 'bias', b, 'solver', 'none') ;
%
%     The latter form is particularly useful when X is a DATASET structure.
%
%   See also: SVM fundamentals,
%   VL_SVMDATASET(), VL_HELP().
[varargout{1:nargout}] = vl_svmtrain(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
clear; close all; clc;

XLIMs = [-2*pi, 2*pi];
YLIMs = [-2 2];
x = linspace(XLIMs(1), XLIMs(2) ,100);
y = sin(x);

plot(x,y,'linewidth',2)
grid on;
hold on;
line([0 0],[YLIMs(1) YLIMs(2)],'color','k')
line([XLIMs(1) XLIMs(2)],[0 0],'color','k')
xlim(XLIMs)
ylim(YLIMs)

ExpansionPoint = 1;

for n= 1:10

    syms x
    f = sin(x);
    T = taylor(f, 'ExpansionPoint', ExpansionPoint, 'Order', n);
    h = fplot(T,XLIMs,'linewidth',2,'color','r');

    title('e^x and its Taylor series');
    tp1 = 0.8;
    tp2 = -1.5;
    if n==1
        t = text(tp1,tp2,['added upto 1st term']);
    elseif n==2
        t = text(tp1,tp2,['added upto 2nd term']);
    elseif n==3
        t = text(tp1,tp2,['added upto 3rd term']);
    else
        t = text(tp1,tp2,['added upto ',num2str(n),'th term']);
    end
    t.FontSize=12;
    pause;
    if n<10
        delete(h);
        delete(t);
    end
end

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function [ fmin, a, b ] = p00_fmin ( a, b, problem, tol )

%*****************************************************************************80
%
%% P00_FMIN seeks a minimizer of a scalar function of a scalar variable.
%
%  Discussion:
%
%    FMIN seeks an approximation to the point where F attains a minimum on
%    the interval (A,B).
%
%    The method used is a combination of golden section search and
%    successive parabolic interpolation.  Convergence is never much
%    slower than that for a Fibonacci search.  If F has a continuous
%    second derivative which is positive at the minimum (which is not
%    at A or B), then convergence is superlinear, and usually of the
%    order of about 1.324....
%
%    The function F is never evaluated at two points closer together
%    than EPS * ABS ( FMIN ) + (TOL/3), where EPS is approximately the
%    square root of the relative machine precision.  If F is a unimodal
%    function and the computed values of F are always unimodal when
%    separated by at least EPS * ABS ( XSTAR ) + (TOL/3), then FMIN
%    approximates the abcissa of the global minimum of F on the
%    interval [A, B] with an error less than 3 * EPS * ABS ( FMIN ) + TOL.
%    If F is not unimodal, then FMIN may approximate a local, but
%    perhaps non-global, minimum to the same accuracy.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Richard Brent,
%    Algorithms for Minimization without Derivatives,
%    Prentice Hall, 1973.
%
%    David Kahaner, Cleve Moler, Steven Nash,
%    Numerical Methods and Software,
%    Prentice Hall, 1988.
%
%  Parameters
%
%    Input, real A, B, the left and right endpoints of the initial interval.
%
%    Input, integer PROBLEM, the index of a problem.
%
%    Input, real TOL, the desired length of the interval of
%    uncertainty of the final result.  TOL must not be negative.
%
%    Output, real FMIN, the abcissa approximating the
%    minimizer of f.
%
%    Output, real A, B, the lower and upper bounds for the minimizer.
%
  c = 0.5 * ( 3.0 - sqrt ( 5.0 ) );
%
%  C is the squared inverse of the golden ratio.
%
%  EPSILON is the square root of the relative machine precision.
%
  epsilon = sqrt ( eps );
%
%  Initialization.
%
  v = a + c * ( b - a );
  w = v;
  x = v;
  e = 0.0;
  fx = p00_f ( problem, x );
  fv = fx;
  fw = fx;
%
%  The main loop starts here.
%
  while ( 1 )

    midpoint = 0.5 * ( a + b );
    tol1 = epsilon * abs ( x ) + tol / 3.0;
    tol2 = 2.0 * tol1;
%
%  Check the stopping criterion.
%
    if ( abs ( x - midpoint ) <= ( tol2 - 0.5 * ( b - a ) ) )
      break
    end
%
%  Is golden-section necessary?
%
    if ( abs ( e ) <= tol1 )
      if ( midpoint <= x )
        e = a - x;
      else
        e = b - x;
      end

      d = c * e;
%
%  Consider fitting a parabola.
%
    else

      r = ( x - w ) * ( fx - fv );
      q = ( x - v ) * ( fx - fw );
      p = ( x - v ) * q - ( x - w ) * r;
      q = 2.0 * ( q - r );
      if ( 0.0 < q )
        p = -p;
      end
      q = abs ( q );
      r = e;
      e = d;
%
%  Choose a golden-section step if the parabola is not advised.
%
      if ( ...
        ( abs ( 0.5 * q * r ) <= abs ( p ) ) | ...
        ( p <= q * ( a - x ) ) | ...
        ( q * ( b - x ) <= p ) )

        if ( midpoint <= x )
          e = a - x;
        else
          e = b - x;
        end

        d = c * e;
%
%  Choose a parabolic interpolation step.
%
      else

        d = p / q;
        u = x + d;

        if ( ( u - a ) < tol2 )
          d = abs ( tol1 ) * r8_sign ( midpoint - x );
        end

        if ( ( b - u ) < tol2 )
          d = abs ( tol1 ) * r8_sign ( midpoint - x );
        end

     end

   end
%
%  F must not be evaluated too close to X.
%
    if ( tol1 <= abs ( d ) )
      u = x + d;
    end

    if ( abs ( d ) < tol1 )
      u = x + abs ( tol1 ) * r8_sign ( d );
    end

    fu = p00_f ( problem, u );
%
%  Update the data.
%
    if ( fu <= fx )

      if ( x <= u )
        a = x;
      else
        b = x;
      end

      v = w;
      fv = fw;
      w = x;
      fw = fx;
      x = u;
      fx = fu;
      continue

    end

    if ( u < x )
      a = u;
    else
      b = u;
    end

    if ( fu <= fw | w == x )
      v = w;
      fv = fw;
      w = u;
      fw = fu;
    elseif ( fu <= fv | v == x | v == w )
      v = u;
      fv = fu;
    end

  end

  fmin = x;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% test omp vs bp

n = 100;
m = 200;
s = 20;
p = 100; % nbr tests

err_bp = [];
err_omp = [];

D = randn(n,m);
D = D./repmat( sqrt(sum(D.^2)), [n 1] );

X0 = randn(m,p);
t = sort(abs(X0));
t = t(end-s,:);
t = repmat(t, [m 1]);
X0 = X0 .* (abs(X0)>t);
Y = D*X0;

options.use_mex = 1;
for k=1:2*s
    progressbar(k,2*s);
    options.nbr_max_atoms = k;
    options.sparse_coder = 'mp';
    X = perform_omp(D,Y,options);
    err_bp(end+1) = norm(Y-D*X,'fro');

    options.sparse_coder = 'omp';
    X = perform_omp(D,Y,options);
    err_omp(end+1) = norm(Y-D*X,'fro');
end

clf;
plot(1:2*s, err_omp, 1:2*s, err_bp);
axis tight;
legend('OMP', 'BP');

author: gpeyre

repo: matlab-toolboxes


text:
function removeDataType(dataTypeName, queryFlag);
%
%     removeDataType(dataTypeName,[queryFlag=1]);
%
% Remove a dataType from mrSESSION and updates open views appropriately.
% This routine depends on dataTYPES being a global variable.
%
% See also: removeDataTypeByIndex
%
% djh, 3/2001
% dar, 3/2007 - added queryflag (default = 1) to pass to saveSession &
% cleanSession - default loads a dialoge asking if you want to overwrite /
% delete.
mrGlobals
if strcmp(dataTypeName,'Original')
    error('removeDataType: Cannot remove original dataType.\n');
end
if ~exist('queryFlag','var')
    queryFlag=1;
end
dataTypeNum = existDataType(dataTypeName);

% if query flag is 1, ask the user: and DON'T DELETE THE DATA
% TYPE if the user says 'no'.
if queryFlag==1
	q = sprintf('Completely delete data type %s?', dataTypeName);
    confirm = questdlg(q, mfilename);
	if ~isequal(confirm, 'Yes')
		fprintf('[%s]: Aborted removal.', mfilename);
		return
	end
end

% If the dataType exists, (1) remove it from mrSESSION, (2) remove it from the
% popup menu options, (3) adjust the value of the curDataType
% and (4) check whether or not to delete the corresponding subdirectories
if dataTypeNum
    % Adjust the mrSESSION structure
    fprintf('removeDataType: removing %s dataType\n',dataTypeName);
    other = setdiff(1:length(dataTYPES),dataTypeNum);
    dataTYPES = dataTYPES(other);

    saveSession;

    if dataTypeNum > 1 % should always happen, except unusual circumstances
        dataTypeNum = dataTypeNum - 1;
    end

    % Delete all the data files for this dataType
    cleanDataType(dataTypeName, 0);

    % Loop through the open views, switch their curDataType appropriately,
    % and update the dataType popups
    INPLANE = resetDataTypes(INPLANE, dataTypeNum);
    VOLUME = resetDataTypes(VOLUME, dataTypeNum);
    FLAT = resetDataTypes(FLAT, dataTypeNum);

	% ras 02/07: check if a mrVista2 session is open (this shouldn't
	% break anything if the toolbox is not installed)
	global GUI;
	if ~isempty(GUI), try, sessionGUI_selectDataType(dataTypeNum); end; end
else
    fprintf('removeDataType: %s dataType not found.\n',dataTypeName);
end
return;


text:
function flag = is_lv_on_the_left(scimat, m, rot)
% IS_LV_ON_THE_LEFT  Check whether the left ventricle is on the left hand
% side of the image.
%
% FLAG = is_lv_on_the_left(SCIMAT, M, ROT)
%
%   SCIMAT is a struct with the whole tissue segmentation (see "help
%   scimat" for details on SCIMAT structs).
%
%   M, ROT are a 3-vector and a (3, 3)-matrix that represent the centroid
%   and main axes of the heart, respectively.
%
%   FLAG is a boolean. It is true of the left ventricle is on the left hand
%   side of the image with respect to the coordinate system defined by ROT.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2010-2014 University of Oxford
% Version: 0.1.2
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(3, 3);
nargoutchk(0, 1);

% cut a whole tissue plane orthogonal to the vertical axis
im = scimat_intersect_plane(scimat, m, rot(:, 3));

% the computed "X axis" of the orientation can be pointing more or less in
% the direction of of X or -X. We compute the angle
alpha = acosd(dot(rot(:, 1), [1 0 0]'));

% now we rotate the intersected plane to get an idea of whether the LV is
% going to be on the left or right hand side of the rotated volume
im = imrotate(im, alpha);

% label connected components
[lab, nlabs] = bwlabel(~im);

% get number of voxels in each label
nvox = zeros( nlabs, 1 );
for L = 1:nlabs
    nvox(L) = sum( lab(:) == L );
end

% sort labels by decreasing size
[nvox, idx] = sort(nvox, 1, 'descend');

% the largest component is the background; the 2nd and 3rd largest should
% be the RV and LV, in either order

% compute the eccentricity of each label (the opposite of the roundness)
stat2 = regionprops(double(lab == idx(2)), 'Eccentricity');
stat3 = regionprops(double(lab == idx(3)), 'Eccentricity');

% we vote whether the 2nd largest component is actually the RV: it is more
% eccentric
if (stat2.Eccentricity > stat3.Eccentricity)
    idxrv = idx(2);
    idxlv = idx(3);
else
    idxrv = idx(3);
    idxlv = idx(2);
end

% get the centroids of both components
crv = regionprops(double(lab == idxrv), 'Centroid');
clv = regionprops(double(lab == idxlv), 'Centroid');

% is the LV placed on the left hand side?
flag = clv.Centroid(1) < crv.Centroid(1);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [RHS, RHSx, RHSy] = ...
    convectionTvdRHS2D(u, phi, FL)
% This function uses the upwind scheme to discretize a 1D
% convection term in the form \grad (u \phi) where u is a face vactor
% It also returns the x and y parts of the matrix of coefficient.
%
% SYNOPSIS:
%   [M, RHS, Mx, My, RHSx, RHSy] = ...
%    convectionTvdTerm2D(u, phi, FL)
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% extract data from the mesh structure
Nx = u.domain.dims(1);
Ny = u.domain.dims(2);
G=reshape(1:(Nx+2)*(Ny+2), Nx+2, Ny+2);
DXp = repmat(u.domain.cellsize.x(2:end-1), 1, Ny);
DYp = repmat(u.domain.cellsize.y(2:end-1)', Nx, 1);
dx=repmat(0.5*(u.domain.cellsize.x(1:end-1)+u.domain.cellsize.x(2:end)), 1, Ny);
dy=repmat(0.5*(u.domain.cellsize.y(1:end-1)+u.domain.cellsize.y(2:end))', Nx, 1);
psiX_p = zeros(Nx+1,Ny);
psiX_m = zeros(Nx+1,Ny);
psiY_p = zeros(Nx,Ny+1);
psiY_m = zeros(Nx,Ny+1);

% define the vectors to stores the sparse matrix data
mnx = Nx*Ny;	mny = Nx*Ny;

% extract the velocity data
% note: size(ux) = [1:m+1, 1:n] and size(uy) = [1:m, 1:n+1]
ux = u.xvalue;
uy = u.yvalue;

% calculate the upstream to downstream gradient ratios for u>0 (+ ratio)
% x direction
dphiX_p = (phi.value(2:Nx+2, 2:Ny+1)-phi.value(1:Nx+1, 2:Ny+1))./dx;
rX_p = dphiX_p(1:end-1,:)./fsign(dphiX_p(2:end,:));
psiX_p(2:Nx+1,:) = 0.5*FL(rX_p).*(phi.value(3:Nx+2,2:Ny+1)-phi.value(2:Nx+1, 2:Ny+1));
psiX_p(1, :) = 0; % left boundary will be handled in the main matrix
% y direction
dphiY_p = (phi.value(2:Nx+1, 2:Ny+2)-phi.value(2:Nx+1, 1:Ny+1))./dy;
rY_p = dphiY_p(:,1:end-1)./fsign(dphiY_p(:,2:end));
psiY_p(:,2:Ny+1) = 0.5*FL(rY_p).*(phi.value(2:Nx+1,3:Ny+2)-phi.value(2:Nx+1, 2:Ny+1));
psiY_p(:,1) = 0; % Bottom boundary will be handled in the main matrix

% calculate the upstream to downstream gradient ratios for u<0 (- ratio)
% x direction
rX_m = dphiX_p(2:end,:)./fsign(dphiX_p(1:end-1,:));
psiX_m(1:Nx,:) = 0.5*FL(rX_m).*(phi.value(1:Nx, 2:Ny+1)-phi.value(2:Nx+1, 2:Ny+1));
psiX_m(Nx+1,:) = 0; % right boundary
% y direction
rY_m = dphiY_p(:,2:end)./fsign(dphiY_p(:,1:end-1));
psiY_m(:,1:Ny) = 0.5*FL(rY_m).*(phi.value(2:Nx+1, 1:Ny)-phi.value(2:Nx+1, 2:Ny+1));
psiY_m(:, Ny+1) = 0; % top boundary will be handled in the main matrix

% reassign the east, west, north, and south velocity vectors for the
% code readability
ue = ux(2:Nx+1,:);		uw = ux(1:Nx,:);
vn = uy(:,2:Ny+1);       vs = uy(:,1:Ny);

% find the velocity direction for the upwind scheme
ue_min = min(ue,0);	ue_max = max(ue,0);
uw_min = min(uw,0);	uw_max = max(uw,0);
vn_min = min(vn,0);	vn_max = max(vn,0);
vs_min = min(vs,0);	vs_max = max(vs,0);

% build the sparse matrix based on the numbering system
rowx_index = reshape(G(2:Nx+1,2:Ny+1),mnx,1); % main diagonal x
rowy_index = reshape(G(2:Nx+1,2:Ny+1),mny,1); % main diagonal y

% calculate the TVD correction term
div_x = -(1./DXp).*((ue_max.*psiX_p(2:Nx+1,:)+ue_min.*psiX_m(2:Nx+1,:))- ...
              (uw_max.*psiX_p(1:Nx,:)+uw_min.*psiX_m(1:Nx,:)));
div_y = -(1./DYp).*((vn_max.*psiY_p(:,2:Ny+1)+vn_min.*psiY_m(:,2:Ny+1))- ...
              (vs_max.*psiY_p(:,1:Ny)+vs_min.*psiY_m(:,1:Ny)));
% define the RHS Vector
RHS = zeros((Nx+2)*(Ny+2),1);
RHSx = zeros((Nx+2)*(Ny+2),1);
RHSy = zeros((Nx+2)*(Ny+2),1);

% assign the values of the RHS vector
RHS(rowx_index) = reshape(div_x+div_y,Nx*Ny,1);
RHSx(rowx_index) = reshape(div_x,Nx*Ny,1);
RHSy(rowy_index) = reshape(div_y,Nx*Ny,1);

end

function phi_out = fsign(phi_in)
% This function checks the value of phi_in and assigns an eps value to the
% elements that are less than or equal to zero, while keeping the signs of
% the nonzero elements
    phi_out = (abs(phi_in)>=eps).*phi_in+eps*(phi_in==0)+eps*(abs(phi_in)<eps).*sign(phi_in);
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function [mappedX, mapping] = em_pca(X, no_dims, max_iter)
%EMPCA Run an EM-based implementation of (probabilistic) PCA
%
%   [mappedX, mapping] = em_pca(X, no_dims)
%
% Performs probabilistic PCA on dataset X in order to reduce its
% dimensionality to no_dims. The dimensionality reduction is performed by
% means of an EM algorithm. The resulting low-dimensional counterpart of X
% is returned in mappedX. Information on the applied mapping (allowing for,
% e.g., out-of-sample extension) is returned in mapping.
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    if ~exist('max_iter', 'var')
        max_iter = 200;
    end

    % Initialize some variables
    [n D] = size(X);                        % data dimensions
    Ez = zeros(no_dims, n);                 % expectation of latent vars
    Ezz = zeros(no_dims, no_dims, n);       % expectation of cov(z)
    Q = Inf;                                % log-likelihood
    mapping = struct;

    % Randomly initialize W and sigma
    W = rand(D, no_dims) * 2;               % factor loadings
    sigma2 = rand(1) * 2;                   % variance ^ 2
    % The covariance of the Gaussian is: C = W * W' + sigma2 * eye(D);

    % Make data zero-mean (possible because data mean is ML estimate for mu)
    mapping.mean = mean(X, 1);
    X = bsxfun(@minus, X, mapping.mean);

    % Compute data covariance and transpose data
    S = cov(X);
    X = X';

    % Perform EM iterations
    converged = 0;
    iter = 0;
    inW = W' * W;
    while ~converged && iter <= max_iter

        % Update iteration number
        iter = iter + 1;
        if rem(iter, 5) == 0
            fprintf('.');
        end

        % Perform E-step
        invM = inv(inW + sigma2 * eye(no_dims));
        for i=1:n
            Ez(:,i)    = invM * W' * X(:,i);
            Ezz(:,:,i) = sigma2 * invM + Ez(:,i) * Ez(:,i)';
        end

        % Perform M-step (maximize mapping W)
        Wp1 = zeros(D, no_dims);
        Wp2 = zeros(no_dims, no_dims);
        for i=1:n
            Wp1 = Wp1 + X(:,i) * Ez(:,i)';
            Wp2 = Wp2 + Ezz(:,:,i);
        end
        W = Wp1 / Wp2;
        inW = W' * W;

        % Perform M-step (maximize discarded variance sigma)
        normX = sum(X .^ 2, 1);
        sigma2 = 0;
        for i=1:n
            sigma2 = sigma2 + (normX(i) - 2 * Ez(:,i)' * W' * X(:,i) + trace(Ezz(:,:,i) * inW));
        end
        sigma2 = (1 / (n * D)) * sigma2;

        % Compute likelihood of new model
        oldQ = Q;
        if iter > 1
            invC = ((1 / sigma2) * eye(D)) - ((1 / sigma2) * W * invM * W');
            detC = det(sigma2 * eye(D)) * det(eye(no_dims) + W' * ((sigma2 .^ -1) * eye(D)) * W);
            Q = (-n / 2) * (D * log(2 * pi) + log(detC) + trace(invC * S));
        end

        % Stop condition to detect convergence
        if abs(oldQ - Q) < 1e-3
            converged = 1;
        end
    end

    % Compute mapped data
    disp(' ');
    mapping.M = (inW \ W')';
    mapping.sigma2 = sigma2;
    mappedX = X' * mapping.M;

repo: cluster


text:
function [cfg, artifact] = ft_artifact_clip(cfg, data)

% FT_ARTIFACT_CLIP scans the data segments of interest for channels that clip. These
% artifacts are detected by the signal being completely flat for a given amount of
% time.
%
% Use as
%   [cfg, artifact] = ft_artifact_clip(cfg)
% with the configuration options
%   cfg.dataset     = string with the filename
% or
%   cfg.headerfile  = string with the filename
%   cfg.datafile    = string with the filename
% and optionally
%   cfg.headerformat
%   cfg.dataformat
%
% Alternatively you can use it as
%   [cfg, artifact] = ft_artifact_clip(cfg, data)
% where the input data is a structure as obtained from FT_PREPROCESSING.
%
% In both cases the configuration should also contain
%   cfg.artfctdef.clip.channel       = Nx1 cell-array with selection of channels, see FT_CHANNELSELECTION for details
%   cfg.artfctdef.clip.pretim        = 0.000;  pre-artifact rejection-interval in seconds
%   cfg.artfctdef.clip.psttim        = 0.000;  post-artifact rejection-interval in seconds
%   cfg.artfctdef.clip.timethreshold = number, minimum duration in seconds of a datasegment with consecutive identical samples to be considered as 'clipped'
%   cfg.artfctdef.clip.amplthreshold = number, minimum amplitude difference in consecutive samples to be considered as 'clipped' (default = 0)
%                                      string, percent of the amplitude range considered as 'clipped' (i.e. '1%')
%   cfg.continuous                   = 'yes' or 'no' whether the file contains continuous data
%
% The output argument "artifact" is a Nx2 matrix comparable to the
% "trl" matrix of FT_DEFINETRIAL. The first column of which specifying the
% beginsamples of an artifact period, the second column contains the
% endsamples of the artifactperiods.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% If you specify this option the input data will be read from a *.mat
% file on disk. This mat files should contain only a single variable named 'data',
% corresponding to the input structure.
%
% See also FT_REJECTARTIFACT, FT_ARTIFACT_CLIP, FT_ARTIFACT_ECG, FT_ARTIFACT_EOG,
% FT_ARTIFACT_JUMP, FT_ARTIFACT_MUSCLE, FT_ARTIFACT_THRESHOLD, FT_ARTIFACT_ZVALUE

% Copyright (C) 2005-2011, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble provenance
ft_preamble loadvar data

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed',    {'datatype', 'continuous'});
cfg = ft_checkconfig(cfg, 'renamed',    {'artfctdef.clip.thresh', 'artfctdef.clip.timethreshold'});
cfg = ft_checkconfig(cfg, 'renamedval', {'continuous', 'continuous', 'yes'});

% set default rejection parameters for clip artifacts if necessary.
if ~isfield(cfg, 'artfctdef'),                    cfg.artfctdef                    = [];    end
if ~isfield(cfg.artfctdef, 'clip'),               cfg.artfctdef.clip               = [];    end
if ~isfield(cfg.artfctdef.clip, 'channel'),       cfg.artfctdef.clip.channel       = 'all'; end
if ~isfield(cfg.artfctdef.clip, 'timethreshold'), cfg.artfctdef.clip.timethreshold = 0.010; end
if ~isfield(cfg.artfctdef.clip, 'amplthreshold'), cfg.artfctdef.clip.amplthreshold = 0.000; end
if ~isfield(cfg.artfctdef.clip, 'pretim'),        cfg.artfctdef.clip.pretim        = 0.000; end
if ~isfield(cfg.artfctdef.clip, 'psttim'),        cfg.artfctdef.clip.psttim        = 0.000; end
if ~isfield(cfg, 'headerformat'),                cfg.headerformat                 = [];    end
if ~isfield(cfg, 'dataformat'),                  cfg.dataformat                   = [];    end

% for backward compatibility
if isfield(cfg.artfctdef.clip, 'sgn')
  cfg.artfctdef.clip.channel = cfg.artfctdef.clip.sgn;
  cfg.artfctdef.clip         = rmfield(cfg.artfctdef.clip, 'sgn');
end

% start with an empty artifact list
artifact = [];

% the data is either passed into the function by the user or read from file with cfg.inputfile
hasdata = exist('data', 'var');

if ~hasdata
  cfg = ft_checkconfig(cfg, 'dataset2files', 'yes');
  cfg = ft_checkconfig(cfg, 'required', {'headerfile', 'datafile'});
  hdr = ft_read_header(cfg.headerfile, 'headerformat', cfg.headerformat);
  trl = cfg.trl;
else
  data = ft_checkdata(data, 'hassampleinfo', 'yes');
  cfg  = ft_checkconfig(cfg, 'forbidden', {'dataset', 'headerfile', 'datafile'});
  hdr  = ft_fetch_header(data);
  if isfield(data, 'sampleinfo')
    trl = data.sampleinfo;
    for k = 1:numel(data.trial)
      trl(k,3) = time2offset(data.time{k}, data.fsample);
    end
  else
    ft_error('the input data does not contain a valid description of the sampleinfo');
  end
end

% set default cfg.continuous
if ~isfield(cfg, 'continuous')
  if hdr.nTrials==1
    cfg.continuous = 'yes';
  else
    cfg.continuous = 'no';
  end
end

% find the channel labels present in the data and their indices
label = ft_channelselection(cfg.artfctdef.clip.channel, hdr.label);
sgnindx = match_str(hdr.label, label);

% make a local copy for convenience
artfctdef = cfg.artfctdef.clip;

ntrl = size(trl,1);
nsgn = length(sgnindx);
for trlop=1:ntrl
  fprintf('searching for clipping artifacts in trial %d\n', trlop);
  % read the data of this trial
  if hasdata
    dat = ft_fetch_data(data,        'header', hdr, 'begsample', trl(trlop,1), 'endsample', trl(trlop,2), 'chanindx', sgnindx);
  else
    dat = ft_read_data(cfg.datafile, 'header', hdr, 'begsample', trl(trlop,1), 'endsample', trl(trlop,2), 'chanindx', sgnindx, 'checkboundary', strcmp(cfg.continuous, 'no'), 'dataformat', cfg.dataformat);
  end
  % get time
  if size(trl,2)>=3
    time = offset2time(trl(trlop,3), hdr.Fs, size(dat,2));
  elseif hasdata
    time = data.time{trlop};
  end
  datflt = preproc(dat, label, time, artfctdef);

  %check if cfg.artfctdef.clip.amplthreshold is an string indicating percentage (e.g. '10%')
  if ~isempty(cfg.artfctdef.clip.amplthreshold) && ischar(cfg.artfctdef.clip.amplthreshold) && cfg.artfctdef.clip.amplthreshold(end)=='%'
    ratio = sscanf(cfg.artfctdef.clip.amplthreshold, '%f%%');
    ratio = ratio/100;
    identical = abs(datflt(:,1:(end-1))-datflt(:,2:end));
    r = range(identical,2);
    for sgnlop=1:length(sgnindx);
      identical(sgnlop,:) = (identical(sgnlop,:)/r(sgnlop))*100;
    end
    identical = identical <= ratio;
  else
    % detect all samples that have the same value as the previous sample
    identical = abs(datflt(:,1:(end-1))-datflt(:,2:end))<=cfg.artfctdef.clip.amplthreshold;
  end

  % ensure that the number of samples does not change
  identical = [identical zeros(nsgn,1)];

  % determine the number of consecutively identical samples
  clip = zeros(size(dat));
  for sgnlop=1:length(sgnindx)
    up = find(diff([0 identical(sgnlop,:)], 1, 2)== 1);
    dw = find(diff([identical(sgnlop,:) 0], 1, 2)==-1);
    for k=1:length(up)
      clip(sgnlop,up(k):dw(k)) = dw(k)-up(k);
    end
  end
  % collapse over cannels
  clip = max(clip,[],1);

  % detect whether there are intervals in which the number of consecutive
  % identical samples is larger than the threshold
  thresh = (clip>=artfctdef.timethreshold*hdr.Fs);

  % remember the thresholded parts as artifacts
  artup = find(diff([0 thresh])== 1) + trl(trlop,1) - 1;
  artdw = find(diff([thresh 0])==-1) + trl(trlop,1) - 1;
  for k=1:length(artup)
    artifact(end+1,:) = [artup(k) artdw(k)];
  end
end

if ~isempty(artifact)
  % add the pretim and psttim to the detected artifacts
  artifact(:,1) = artifact(:,1) - artfctdef.pretim * hdr.Fs;
  artifact(:,2) = artifact(:,2) + artfctdef.psttim * hdr.Fs;
end

% remember the details that were used here
cfg.artfctdef.clip          = artfctdef;
cfg.artfctdef.clip.label    = label;
cfg.artfctdef.clip.trl      = trl;
cfg.artfctdef.clip.artifact = artifact;

% do the general cleanup and bookkeeping at the end of the function
ft_postamble provenance
ft_postamble previous data
ft_postamble savevar

author: spm


text:
function B = cleanpoly(A,tol,deg)
% function B = cleanpoly(A,tol,deg)
%
% DESCRIPTION
%   Cleans up the input polynomial.  The output polynomial includes only
%   terms whose coefficients have magnitude greater than or equal to TOL
%   and whose monomial degree is specified by DEG.
%
% INPUTS
%   A: polynomial
%   tol: scalar double specifying the coefficient tolerance
%   deg: vector of non-negative integers specifying the degrees of
%        mononmials to retain. Alternatively deg can be an N-by-2
%        cell array with deg{i,1} specifying a variable and
%        deg{i,2} specifying a vector of non-negative integers.
%        This will retain only monomials whose degree in variable
%        deg{i,1} is specified in deg{i,2}.
%
% OUTPUTS
%   B: polynomial which only contains the terms of A whose coefficients
%      have magnitude greater than or equal to tol and whose monomial
%      degree is listed in deg.
%
% SYNTAX
%   B=cleanpoly(A,tol);
%   B=cleanpoly(A,[],deg);
%   B=cleanpoly(A,tol,deg);
%
% EXAMPLE
%   pvar x1 x2 u;
%   p = 9*u^3 + u*x1^2 + 1e-6*u^2*x1*x2 + 1e-5*u*x2^2 + 2*x1^3 ...
%        - x1*x2 + 3*u + x1 + 2*x2;
%
%   % Remove terms whose coefficients has magnitude < tol
%   tol = 1e-4;
%   p1 = cleanpoly(p,tol)
%
%   % Retain linear and quadratic terms
%   p2 = cleanpoly(p,[],1:2)
%
%   % Retain terms linear in u but of degree 0,1,2,3 in x1 and x2
%   p3 = cleanpoly(p,[],{x1, 0:3; x2, 0:3; u 1})

% 1/28/2008: PJS  Added option to grab terms of a specified degree
% 4/22/2009: PJS  Added functionality to specify deg as a cell array

if nargin<3
    deg = [];
end

if ~isempty(tol)
    if isa(tol,'double') && isscalar(tol) && tol>=0
        Acoef = get(A,'Coefficient');
        idx = find(abs(Acoef) < tol);
        if ~isempty(idx)
            Acoef(idx) = 0;
            chkval = 0; % skip validity check
            A=polynomial(Acoef,A.degmat,A.varname,size(A),chkval);
        end
    else
        error('tol must be a non-negative double');
    end
end

if ~isempty(deg)

    if isa(deg,'double') && ndims(deg)==2 && ...
            all(floor(deg)==ceil(deg)) && all(deg>=0)

        % deg is a vector of non-negative integers

        Adeg = sum(A.degmat,2);
        idx = [];
        for i1=1:length(deg)
            %idx=find(Adeg~=deg(i1));
            idx=[idx; find(Adeg==deg(i1))];
        end
        idx = setdiff(1:length(Adeg),idx);

        Acoef = A.coefficient;
        Acoef(idx,:) = 0;
        chkval = 0; % skip validity check
        A=polynomial(Acoef,A.degmat,A.varname,size(A),chkval);

    elseif iscell(deg) && ndims(deg)==2 && size(deg,2)==2
        % deg is a cell array of pvars and vectors of non-neg integers

        Acoef = A.coefficient;
        for i1=1:size(deg,1)
            var = deg{i1,1};
            if ispvar(var) && all(size(var)==[1 1])
                var = char(var);
            elseif ~ischar(var)
                error('First column of cell array deg must contain strings or pvars');
            end
            vidx = find(strcmp(A.varname,var));

            if ~isempty(vidx)
                degi = deg{i1,2};
                if ~( isa(degi,'double') && ndims(degi)==2 && ...
                        all(floor(degi)==ceil(degi)) && all(degi>=0) )
                    error(['Second column of cell array deg must '...
                        'contain a vector of non-negative integers.']);
                end

                Adeg = A.degmat(:,vidx);
                cidx = [];
                for i2=1:length(degi)
                    cidx=[cidx; find(Adeg==degi(i2))];
                end
                cidx = setdiff(1:length(Adeg),cidx);
                Acoef(cidx,:) = 0;
            end
        end
        A=polynomial(Acoef,A.degmat,A.varname,size(A));
    else
        error('deg must be a vector of non-negative integers or an Nx2 cell array');
    end

end

% Combine to remove terms
B = combine(A);

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [dag,best_score] = learn_struct_hc(data, nodesizes, seeddag, varargin)
%
% LEARN_STRUCT_HC(data,seeddag) learns a structure of Bayesian net by Hill Climbing.
% dag = learn_struct_hc(data, nodesizes, seeddag)
%
% dag: the final structurre matrix
% Data : training data, data(i,m) is the m obsevation of node i
% Nodesizes: the size array of different nodes
% seeddag: given seed Dag for hill climbing, optional
%
% by Gang Li @ Deakin University (gli73@hotmail.com)

[N ncases] = size(data);
if (nargin < 3 )
    seeddag = zeros(N,N); % mk_rnd_dag(N); %call BNT function
elseif ~acyclic(seeddag)
    seeddag = mk_rnd_dag(N); %zeros(N,N);
end;

% set default params
scoring_fn = 'bic';
verbose  = 'yes';

% get params
args = varargin;
nargs = length(args);
if length(args) > 0
    if isstr(args{1})
    	for i = 1:2:nargs
    		switch args{i}
    		case 'scoring_fn', scoring_fn = args{i+1};
    		case 'verbose',  verbose  = strcmp(args{i+1},'yes');
    		end;
    	end;
    end;
end;

done = 0;
best_score = score_dags(data,nodesizes, {seeddag},'scoring_fn',scoring_fn);
while ~done
    [dags,op,nodes] = mk_nbrs_of_dag(seeddag);
    nbrs = length(dags);
    scores = score_dags(data, nodesizes, dags,'scoring_fn',scoring_fn);
    max_score = max(scores);
    new = find(scores == max_score );
    if ~isempty(new) & (max_score > best_score)
        p = sample_discrete(normalise(ones(1, length(new))));
        best_score = max_score;
        seeddag = dags{new(p)};
    else
        done = 1;
    end;
end;

dag = seeddag;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [out1, out2] = templateFunction(arg1, arg2, arg3)
% This function does this and that
%
% USAGE:
%
%   [out1, out2] = templateFunction(arg1, arg2, arg3)
%
% INPUTS:
%    arg1: a short description for this argument
%    arg2: another argument
%
% OPTIONAL INPUT:
%    arg3: optional argument
%
% OUTPUT:
%    out1: the product of arg1 and arg2
%
% OPTIONAL OUTPUT:
%    out2: the product of arg1, arg2 and arg3
%
% EXAMPLE:
%
%   [out1, out2] = templateFunction(1, 2, 3);
%   >> out1 = 2
%   >> out2 = 6
%
% .. Authors: - myself
%             - laurent
%
if nargin < 3
    arg3 = 1;
end

out1 = arg1 * arg2;

if nargout > 1
    out2 = out1 * arg3;
end

author: opencobra

repo: cobratoolbox


text:
function a = hartley ( n )

%*****************************************************************************80
%
%% HARTLEY returns the HARTLEY matrix.
%
%  Formula:
%
%    A(I,J) = SIN ( 2*PI*(i-1)*(j-1)/N ) + COS( 2*PI*(i-1)*(j-1)/N )
%
%  Example:
%
%    N = 5
%
%    1.0000    1.0000    1.0000    1.0000    1.0000
%    1.0000    1.2601   -0.2212   -1.3968   -0.6420
%    1.0000   -0.2212   -0.6420    1.2601   -1.3968
%    1.0000   -1.3968    1.2601   -0.6420   -0.2212
%    1.0000   -0.6420   -1.3968   -0.2212    1.2601
%
%  Properties:
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A arises in the Hartley transform.
%
%    A * A = N * I, in other words, A is "almost" involutional,
%    and inverse ( A ) = ( 1 / N ) * A.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    D Bini, P Favati,
%    On a matrix algebra related to the discrete Hartley transform,
%    SIAM Journal on Matrix Analysis and Applications,
%    Volume 14, 1993, pages 500-507.
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n

      angle = 2.0 * pi * ( i - 1 ) * ( j - 1 ) / n;

      a(i,j) = sin ( angle ) + cos ( angle );

    end
  end

  return
end


text:
function pde = elasticitydata(para)
%% ELASTICITYDATA data for elasticity problem
%
% - mu \Delta u - (mu + lambda) grad(div u) = f    in \Omega
%                                         u = g_D  on \partial \Omega
%
% Created by Huayi Wei Monday, 27 June 2011.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if nargin == 0
    lambda = 1;
    mu = 1;
else
    if ~isstruct(para)
        exit('we need a struct data');
    end
    if ~isfield(para,'lambda') || isempty(para.lambda)
        lambda = 1;
    else
        lambda = para.lambda;
    end
    if ~isfield(para,'mu') || isempty(para.mu)
        mu = 1;
    else
        mu = para.mu;
    end
end

pde = struct('lambda',lambda,'mu',mu, 'f', @f, 'exactu',@exactu,'g_D',@g_D);
%%%%%% subfunctions %%%%%%
    function z = f(p)
       z = mu*2*pi^2.*exactu(p);
    end
    function z = exactu(p)
       x = p(:,1); y = p(:,2);
       z = [cos(pi*x).*cos(pi*y), sin(pi*x).*sin(pi*y)];
    end
    function z = g_D(p)
       z = exactu(p);
    end
end

repo: Aero-matlab


text:
clc;
close all;
clearvars;
rng('default');
n = 10;
% Let's get a symmetric positive definite matrix
A = gallery('gcdmat', n);
% Compute its Cholesky decomposition
U = chol(A);
b = randn(n, 1);
% solve the problem
x = U \b
x = spx.la.tris.back_col(U, b)

author: indigits

repo: sparse-plex