text:
function test_failed=test_frames
%TEST_FRAMES  Test the frames methods

test_failed=0;

disp(' ===============  TEST_FRAMES ================');
global LTFAT_TEST_TYPE;

tolchooser.double=1e-9;
tolchooser.single=2e-4;
tolerance = tolchooser.(LTFAT_TEST_TYPE);

% Iterative algorithms need a bigger tolerance
tolchooseriter.double=1e-1;
tolchooseriter.single=1e-1;
toleranceiter = tolchooseriter.(LTFAT_TEST_TYPE);

Fr=cell(1,32);

L=200;
Fr{1}  = frame('dgt','gauss',10,20);
Fr{2}  = frame('dgtreal','gauss',10,20);
Fr{3}  = frame('dwilt','gauss',20);
Fr{4}  = frame('wmdct','gauss',20);
Fr{5}  = frame('gen',tester_crand(200,300),20);

Fr{6}  = frametight(frame('dgt','gauss',10,20));
Fr{7}  = frametight(frame('dgtreal','gauss',10,20));
Fr{8}  = frametight(frame('dwilt','gauss',20));
Fr{9}  = frametight(frame('wmdct','gauss',20));
Fr{10} = frametight(frame('gen',tester_crand(200,300),20));

Fr{11} = frame('dft');
Fr{12} = frame('dcti');
Fr{13} = frame('dctii');
Fr{14} = frame('dctiii');
Fr{15} = frame('dctiv');
Fr{16} = frame('dsti');
Fr{17} = frame('dstii');
Fr{18} = frame('dstiii');
Fr{19} = frame('dstiv');

% Repeat generation of the filters until they have a nice condition number
condnum = 1e10;
while condnum > 1e3
   gfilt={tester_rand(30,1),...
          tester_rand(20,1),...
          tester_rand(15,1),...
          tester_rand(10,1)};

   gfilt=cellfun(@(gEl) cast(gEl,'double'),gfilt,'UniformOutput',0);

   % These two frames might be badly conditioned,
   Fr{20} = frame('ufilterbank',    gfilt,3,4);
   Fr{21} = frame('ufilterbankreal',gfilt,3,4);

   condnum = framebounds(Fr{20},128);
end

Fr{22} = frame('dgt','gauss',4,6,'lt',[1 2]);
Fr{23} = frame('identity');
Fr{24} = frame('fusion',[1 1],Fr{1},Fr{1});
Fr{25} = frametight(frame('dgt','hamming',10,20));
Fr{26} = frametight(frame('wmdct','hamming',20));

g={randn(30,1),randn(50,1),randn(70,1),randn(90,1)};
a=[20,40,60,80];
M=[30,50,70,100];

Fr{27} = frametight(frame('nsdgt',g,a,M));
Fr{28} = frametight(frame('unsdgt',g,a,100));
Fr{29} = frametight(frame('nsdgtreal',g,a,M));
Fr{30} = frametight(frame('unsdgtreal',g,a,100));

Fr{31} = frametight(frame('dftreal'));

Fr{32} = frame('fwt','ana:spline2:2',5);
Fr{33} = frame('wfbt',{'syn:spline2:2',5});
%Fr{34} = frame('wpfbt',{'db4',5});
%Fr{35} = frame('wpfbt',{'db4',5});
%Fr{36} = frame('wpfbt',{'db4',5});

Fr{37} = frame('ufwt','db4',4);
Fr{38} = frame('ufwt','db4',4,'scale');
Fr{39} = frame('ufwt','db4',4,'noscale');

Fr{40} = frame('uwfbt',{'db4',4});
Fr{41} = frame('uwfbt',{'db4',4},'scale');
Fr{42} = frame('uwfbt',{'db4',4},'noscale');

Fr{43} = frame('uwpfbt',{'db4',4});
Fr{44} = frame('uwpfbt',{'db4',4},'scale');
Fr{45} = frame('uwpfbt',{'db4',4},'noscale');

%Fr{36} = frame('uwfbt',{'db4',5});
%Fr{37} = frame('uwpfbt',{'db4',5});

% The tensor frame implementation is currenly broken
%Fr{33} = frame('tensor',Fr{11});

%Fr{31} = frame('filterbank',     gfilt,[4 3 2 2],4);
%Fr{32} = frame('filterbankreal', gfilt,[4 3 2 2],4);

Fr{60} = frame('erbletfb',44100,L,'real','regsampling');
Fr{61} = frame('erbletfb',44100,L,'complex','regsampling');

Fr{62} = frame('erbletfb',44100,L,'real','fractional');
Fr{63} = frame('erbletfb',44100,L,'complex','fractional');

Fr{64} = frametight(Fr{60});
Fr{65} = frametight(Fr{62});

Fr{66} = frame('cqtfb',44100,200,20000,20,L,'real','regsampling');
Fr{67} = frame('cqtfb',44100,200,20000,20,L,'complex','regsampling');

Fr{68} = frame('cqtfb',44100,200,20000,20,L,'real','fractional');
Fr{69} = frame('cqtfb',44100,200,20000,20,L,'complex','fractional');

Fr{70} = frametight(Fr{66});
Fr{71} = frametight(Fr{67});

for cmpx = {'real','complex'}

    if strcmp(cmpx{1},'real')
        f=tester_rand(L,1);
    else
        f=tester_crand(L,1);
    end

for ii=1:numel(Fr)

  F=Fr{ii};

  % To avoid holes in Fr
  if isempty(F)
    continue;
  end;

  % Do not test real-only frames with complex arrays
  if strcmp(cmpx{1},'complex') && F.realinput
      continue;
  end

  Fd=framedual(F);

  c=frana(F,f);
  r=frsyn(Fd,c);
  res=norm(r(1:L)-f);

  lendiff = size(c,1) - frameclength(F,L);
  [test_failed,fail]=ltfatdiditfail(lendiff ,test_failed);
  s=sprintf(['FRAMES CLENGTH        frameno:%3i %s %0.5g %s'],ii,F.type,lendiff,fail);
  disp(s);

  [test_failed,fail]=ltfatdiditfail(res,test_failed,tolerance);
  s=sprintf(['FRAMES DUAL REC       frameno:%3i %s %0.5g %s'],ii,F.type,res,fail);
  disp(s);

  % Checking equality F == framedual(framedual(F))
  Fdd = framedual(framedual(F));
  cdd = frana(Fdd,f);
  res_dd = norm(cdd-c);

  [test_failed,fail]=ltfatdiditfail(res_dd,test_failed,tolerance);
  s=sprintf(['FRAMES DUAL DUAL      frameno:%3i %s %0.5g %s'],ii,F.type,res_dd,fail);
  disp(s);

  F2=frameaccel(F,L);
  F2d=frameaccel(Fd,L);

  c=frana(F2,f);
  r=frsyn(F2d,c);
  res=norm(r(1:L)-f);

  [test_failed,fail]=ltfatdiditfail(res,test_failed,tolerance);
  s=sprintf(['FRAMES ACCEL DUAL REC frameno:%3i %s %0.5g %s'],ii,F.type,res,fail);
  disp(s);

  % Checking equality F == framedual(framedual(F)) after acceleration
  F2dd = framedual(framedual(F2));
  c2dd = frana(Fdd,f);
  res2_dd = norm(c2dd-c);

  [test_failed,fail]=ltfatdiditfail(res2_dd,test_failed,tolerance);
  s=sprintf(['FRAMES ACCEL DUAL DUAL      frameno:%3i %s %0.5g %s'],ii,F.type,res2_dd,fail);
  disp(s);

  % Test that framebounds are able to run, not actual resting is done on
  % the values.
  [A,B]=framebounds(F,L);

  %% Test iterative analysis and synthesis
  r=frsyniter(F,c);

  res=norm(r(1:L)-f);
  [test_failed,fail]=ltfatdiditfail(res,test_failed,toleranceiter);
  s=sprintf(['FRSYNITER             frameno:%3i %s %0.5g %s'],ii,F.type,res,fail);
  disp(s);

  c2=franaiter(Fd,f);
  res=norm(c2-c);
  [test_failed,fail]=ltfatdiditfail(res,test_failed,toleranceiter);
  s=sprintf(['FRANAITER             frameno:%3i %s %0.5g %s'],ii,F.type,res,fail);
  disp(s);

  %% Test matrix representations
  if (~F.realinput)
    LL=framelength(F,L);
    G=frsynmatrix(F,LL);
    res=norm(c-G'*postpad(f,LL));

    [test_failed,fail]=ltfatdiditfail(res,test_failed,tolerance);
    s=sprintf(['FRAMES ANA MATRIX     frameno:%3i %s %0.5g %s'],ii,F.type,res,fail);
    disp(s);

    % We create a different set of coefficients here.
    % The old code used c and failed the test for some filterbank frames
    ctmp = tester_crand(numel(c),1);
    res=norm(frsyn(F,ctmp)-G*ctmp);
    [test_failed,fail]=ltfatdiditfail(res,test_failed,tolerance);
    s=sprintf(['FRAMES SYN MATRIX     frameno:%3i %s %0.5g %s'],ii,F.type,res,fail);
    disp(s);

  end;

  %% Test the frame multipliers: test framemul, framemuladj and
  %% iframemul
  if F.realinput
      m=1+0.01*tester_rand(size(c,1),1);
  else
      m=1+1i+0.01*tester_crand(size(c,1),1);
  end;
  ff=framemul(f,F,Fd,m);
  fr=iframemul(ff,F,Fd,m);
  res=norm(f-fr(1:L))/norm(f);
  [test_failed,fail]=ltfatdiditfail(res,test_failed,tolerance);
  s=sprintf('IFRAMEMUL             frameno:%3i %s %0.5g %s',ii, ...
            F.type,res,fail);
  disp(s);

end;
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


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

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

% Last Modified by GUIDE v2.5 27-Mar-2004 20:18:45

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @histGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @histGUI_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 histGUI is made visible.
function histGUI_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to histGUI (see VARARGIN)

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

% Set the figure background color
set(hObject,'Color',[0.6 1 1]);

% Set the initial number of bins
handles.nbins = 11;

% Update handles structure
guidata(hObject, handles);

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

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

function NBins_Callback(hObject, eventdata, handles)
% hObject    handle to NBins (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get number of bins, round to integer, and update field
nbins = str2num(get(hObject,'String'));
nbins = round(nbins);
if nbins < 1
   nbins = 1;
end
set (handles.NBins,'String',int2str(nbins));

% Save in handles structure
handles.nbins = nbins;
guidata(gcbf, handles);

% Re-display data, if available
if handles.nbins > 0 & ~isempty(handles.x)

   % Create histogram
   hist(handles.x,handles.nbins);

end

% --- Executes during object creation, after setting all properties.
function NBins_CreateFcn(hObject, eventdata, handles)
% hObject    handle to NBins (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --------------------------------------------------------------------
function Open_Callback(hObject, eventdata, handles)
% hObject    handle to Open (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get file name
[filename,path] = uigetfile('*.dat','Load Data File');
if filename ~= 0

   % Read data
   x = textread([path filename],'%f');

   % Save in handles structure
   handles.x = x;
   guidata(gcbf, handles);

   % Create histogram
   hist(handles.x,handles.nbins);

   % Set axis labels
   xlabel('\bfValue');
   ylabel('\bfCount');

   % Calculate statistics
   ave = mean(x);
   med = median(x);
   sd  = std(x);
   n   = length(x);

   % Update fields
   set (handles.MeanData,'String',sprintf('%7.2f',ave));
   set (handles.MedianData,'String',sprintf('%7.2f',med));
   set (handles.StdDevData,'String',sprintf('%7.2f',sd));
   set (handles.TitleString,'String',['Histogram (N = ' int2str(n) ')']);

end

% --------------------------------------------------------------------
function Exit_Callback(hObject, eventdata, handles)
% hObject    handle to Exit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

close(gcbf);

% --------------------------------------------------------------------
function File_Callback(hObject, eventdata, handles)
% hObject    handle to File (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


text:
function [x,fval,exitflag,output,grad]=fminlbfgs(funfcn,x_init,optim)
%FMINLBFGS finds a local minimum of a function of several variables.
%   This optimizer is developed for image registration methods with large
%	amounts of unknown variables.
%
%   Optimization methods supported:
%	- Quasi Newton Broyden-Fletcher-Goldfarb-Shanno (BFGS)
%   - Limited memory BFGS (L-BFGS)
%   - Steepest Gradient Descent optimization.
%
%   [X,FVAL,EXITFLAG,OUTPUT,GRAD] = FMINLBFGS(FUN,X0,OPTIONS)
%
%   Inputs,
%		FUN: Function handle or string which is minimized, returning an
%				error value and optional the error gradient.
%		X0: Initial values of unknowns can be a scalar, vector or matrix
%	 (optional)
%		OPTIONS: Structure with optimizer options, made by a struct or
%				optimset. (optimset doesnot support all input options)
%
%   Outputs,
%		X : The found location (values) which minimize the function.
%		FVAL : The minimum found
%		EXITFLAG : Gives value, which explain why the minimizer stopt
%		OUTPUT : Structure with all important ouput values and parameters
%		GRAD : The gradient at this location
%
%   Extended description of input/ouput variables
%   OPTIONS,
%		OPTIONS.GoalsExactAchieve : If set to 0, a line search method is
%               used which uses a few function calls to do a good line
%               search. When set to 1 a normal line search method with Wolfe
%				conditions is used (default).
%		OPTIONS.GradConstr, Set this variable to true if gradient calls are
%				cpu-expensive (default). If false more gradient calls are
%				used and less function calls.
%	    OPTIONS.HessUpdate : If set to 'bfgs', Broyden�Fletcher�Goldfarb�Shanno
%				optimization is used (default), when the number of unknowns is
%				larger then 3000 the function will switch to Limited memory BFGS,
%				or if you set it to 'lbfgs'. When set to 'steepdesc', steepest
%				decent optimization is used.
%		OPTIONS.StoreN : Number of itterations used to approximate the Hessian,
%			 	in L-BFGS, 20 is default. A lower value may work better with
%				non smooth functions, because than the Hessian is only valid for
%				a specific position. A higher value is recommend with quadratic equations.
%		OPTIONS.GradObj : Set to 'on' if gradient available otherwise finited difference
%				is used.
%     	OPTIONS.Display : Level of display. 'off' displays no output; 'plot' displays
%				all linesearch results in figures. 'iter' displays output at  each
%               iteration; 'final' displays just the final output; 'notify'
%				displays output only if the function does not converge;
%	    OPTIONS.TolX : Termination tolerance on x, default 1e-6.
%	    OPTIONS.TolFun : Termination tolerance on the function value, default 1e-6.
%		OPTIONS.MaxIter : Maximum number of iterations allowed, default 400.
% 		OPTIONS.MaxFunEvals : Maximum number of function evaluations allowed,
%				default 100 times the amount of unknowns.
%		OPTIONS.DiffMaxChange : Maximum stepsize used for finite difference gradients.
%		OPTIONS.DiffMinChange : Minimum stepsize used for finite difference gradients.
%		OPTIONS.OutputFcn : User-defined function that an optimization function calls
%				at each iteration.
%		OPTIONS.rho : Wolfe condition on gradient (c1 on wikipedia), default 0.01.
%		OPTIONS.sigma : Wolfe condition on gradient (c2 on wikipedia), default 0.9.
%		OPTIONS.tau1 : Bracket expansion if stepsize becomes larger, default 3.
%		OPTIONS.tau2 : Left bracket reduction used in section phase,
%		default 0.1.
%		OPTIONS.tau3 : Right bracket reduction used in section phase, default 0.5.
%   FUN,
%		The speed of this optimizer can be improved by also providing
%   	the gradient at X. Write the FUN function as follows
%   	function [f,g]=FUN(X)
%       	f , value calculation at X;
%   	if ( nargout > 1 )
%       	g , gradient calculation at X;
%   	end
%	EXITFLAG,
%		Possible values of EXITFLAG, and the corresponding exit conditions
%		are
%  		1, 'Change in the objective function value was less than the specified tolerance TolFun.';
%  		2, 'Change in x was smaller than the specified tolerance TolX.';
%  		3, 'Magnitude of gradient smaller than the specified tolerance';
%  		4, 'Boundary fminimum reached.';
%  		0, 'Number of iterations exceeded options.MaxIter or number of function evaluations exceeded options.FunEvals.';
%  		-1, 'Algorithm was terminated by the output function.';
%  		-2, 'Line search cannot find an acceptable point along the current search';
%
%   Examples
%       options = optimset('GradObj','on');
%       X = fminlbfgs(@myfun,2,options)
%
%   	% where myfun is a MATLAB function such as:
%       function [f,g] = myfun(x)
%       f = sin(x) + 3;
%	    if ( nargout > 1 ), g = cos(x); end
%
%   See also OPTIMSET, FMINSEARCH, FMINBND, FMINCON, FMINUNC, @, INLINE.
%
%   Function is written by D.Kroon University of Twente (March 2009)

% Contributed to SOM Toolbox vs2, December 3rd, 2012 by Alexander Schulz
% See license file BSD_license.txt for license information.
% http://www.cis.hut.fi/projects/somtoolbox/

% Read Optimalisation Parameters
defaultopt = struct('Display','final','HessUpdate','bfgs','GoalsExactAchieve',1,'GradConstr',true,  ...
			'TolX',1e-6,'TolFun',1e-6,'GradObj','off','MaxIter',400,'MaxFunEvals',100*numel(x_init)-1,  ...
			'DiffMaxChange',1e-1,'DiffMinChange',1e-8,'OutputFcn',[], ...
			'rho',0.0100,'sigma',0.900,'tau1',3,'tau2', 0.1, 'tau3', 0.5,'StoreN',20);

if (~exist('optim','var'))
    optim=defaultopt;
else
    f = fieldnames(defaultopt);
    for i=1:length(f),
        if (~isfield(optim,f{i})||(isempty(optim.(f{i})))), optim.(f{i})=defaultopt.(f{i}); end
    end
end

% Initialize the data structure
data.fval=0;
data.gradient=0;
data.fOld=[];
data.xsizes=size(x_init);
data.numberOfVariables = numel(x_init);
data.xInitial = x_init(:);
data.alpha=1;
data.xOld=data.xInitial;
data.iteration=0;
data.funcCount=0;
data.gradCount=0;
data.exitflag=[];
data.nStored=0;
%%% data.timeTotal=tic;
data.timeExtern=0;
% Switch to L-BFGS in case of more than 3000 unknown variables
if(false && optim.HessUpdate(1)=='b')
    if(data.numberOfVariables<3000),
        optim.HessUpdate='bfgs';
    else
        optim.HessUpdate='lbfgs';
    end
end

if(optim.HessUpdate(1)=='l')
    data.deltaX=zeros(data.numberOfVariables,optim.StoreN);
    data.deltaG=zeros(data.numberOfVariables,optim.StoreN);
    data.saveD=zeros(data.numberOfVariables,optim.StoreN);
end

exitflag=[];

% Display column headers
if(strcmp(optim.Display,'iter'))
    disp('     Iteration  Func-count   Grad-count         f(x)         Step-size');
end

% Calculate the initial error and gradient
data.initialStepLength=1;
[data,fval,grad]=gradient_function(data.xInitial,funfcn, data, optim);
data.gradient=grad;
data.dir = -data.gradient;
data.gOld=grad;
data.fInitial = fval;
data.fPrimeInitial= data.gradient'*data.dir(:);

gNorm = norm(data.gradient,Inf);  % Norm of gradient
data.initialStepLength = min(1/gNorm,5);

% Show the current iteration
if(strcmp(optim.Display,'iter'))
        s=sprintf('     %5.0f       %5.0f       %5.0f       %13.6g    ',data.iteration,data.funcCount,data.gradCount,data.fInitial); disp(s);
end

% Hessian intialization
if(optim.HessUpdate(1)=='b')
	data.Hessian=eye(data.numberOfVariables);
end

% Call output function
if(call_output_function(data,optim,'init')), exitflag=-1; end

% Start Minimizing
while(true)
    % Update number of itterations
    data.iteration=data.iteration+1;

    % Set current lineSearch parameters
    data.TolFunLnS = eps(max(1,abs(data.fInitial )));
    data.fminimum = data.fInitial - 1e16*(1+abs(data.fInitial));

	% Make arrays to store linesearch results
    data.storefx=[]; data.storepx=[]; data.storex=[]; data.storegx=[];

    % If option display plot, than start new figure
    if(optim.Display(1)=='p'), figure, hold on; end

    % Find a good step size in the direction of the gradient: Linesearch
    if(optim.GoalsExactAchieve==1)
		data=linesearch(funfcn, data,optim);
    else
        data=linesearch_simple(funfcn, data, optim);
    end

	% Make linesearch plot
	if(optim.Display(1)=='p');
		plot(data.storex,data.storefx,'r*');
		plot(data.storex,data.storefx,'b');

		alpha_test= linspace(min(data.storex(:))/3, max(data.storex(:))*1.3, 10);
		falpha_test=zeros(1,length(alpha_test));
        for i=1:length(alpha_test)
			[data,falpha_test(i)]=gradient_function(data.xInitial(:)+alpha_test(i)*data.dir(:),funfcn, data, optim);
        end
		plot(alpha_test,falpha_test,'g');
        plot(data.alpha,data.f_alpha,'go','MarkerSize',8);
	end

    % Check if exitflag is set
    if(~isempty(data.exitflag)),
        exitflag=data.exitflag;
        data.xInitial=data.xOld;
        data.fInitial=data.fOld;
        data.gradient=data.gOld;
        break,
    end;

    % Update x with the alpha step
    data.xInitial = data.xInitial + data.alpha*data.dir;

    % Set the current error and gradient
    data.fInitial =  data.f_alpha;
	data.gradient = data.grad;

    % Set initial steplength to 1
    data.initialStepLength = 1;

    gNorm = norm(data.gradient,Inf);  % Norm of gradient

    % Set exit flags
    if(gNorm <optim.TolFun), exitflag=1; end
    if(max(abs(data.xOld-data.xInitial)) <optim.TolX), exitflag=2; end
    if(data.iteration>=optim.MaxIter), exitflag=0; end

    % Check if exitflag is set
    if(~isempty(exitflag)), break, end;

    % Update the inverse Hessian matrix
    if(optim.HessUpdate(1)~='s')
        % Do the Quasi-Neton Hessian update.
        data = updateQuasiNewtonMatrix_LBFGS(data,optim);
    else
        data.dir = -data.gradient;
    end

    % Derivative of direction
    data.fPrimeInitial= data.gradient'*data.dir(:);

    % Call output function
    if(call_output_function(data,optim,'iter')), exitflag=-1; end

    % Show the current iteration
    if(strcmp(optim.Display(1),'i')||strcmp(optim.Display(1),'p'))
        s=sprintf('     %5.0f       %5.0f       %5.0f       %13.6g   %13.6g',data.iteration,data.funcCount,data.gradCount,data.fInitial,data.alpha); disp(s);
    end

    % Keep the variables for next iteration
    data.fOld=data.fInitial;
    data.xOld=data.xInitial;
    data.gOld=data.gradient;
end
% Set output parameters
fval=data.fInitial;
grad=data.gradient;
x = data.xInitial;

% Reshape x to original shape
x=reshape(x,data.xsizes);

% Call output function
if(call_output_function(data,optim,'done')), exitflag=-1; end

% Make exist output structure
if(optim.HessUpdate(1)=='b'), output.algorithm='Broyden-Fletcher-Goldfarb-Shanno (BFGS)';
elseif(optim.HessUpdate(1)=='l'), output.algorithm='limited memory BFGS (L-BFGS)';
else output.algorithm='Steepest Gradient Descent';
end
output.message=getexitmessage(exitflag);
output.iteration = data.iteration;
output.funccount = data.funcCount;
output.fval = data.fInitial;
output.stepsize = data.alpha;
output.directionalderivative = data.fPrimeInitial;
output.gradient = reshape(data.gradient, data.xsizes);
output.searchdirection = data.dir;
%%% output.timeTotal=toc;
output.timeExtern=data.timeExtern;
%%%oupput.timeIntern=output.timeTotal-output.timeExtern;
% Display final results
if(~strcmp(optim.Display,'off'))
    disp('    Optimizer Results')
    disp(['        Algorithm Used: ' output.algorithm]);
    disp(['        Exit message : ' output.message]);
    disp(['        iterations : '  int2str(data.iteration)]);
    disp(['        Function Count : ' int2str(data.funcCount)]);
    disp(['        Minimum found : ' num2str(fval)]);
%%%    disp(['        Intern Time : ' num2str(oupput.timeIntern) ' seconds']);
%%%    disp(['        Total Time : ' num2str(output.timeTotal) ' seconds']);
end

function message=getexitmessage(exitflag)
    switch(exitflag)
        case 1, message='Change in the objective function value was less than the specified tolerance TolFun.';
        case 2, message='Change in x was smaller than the specified tolerance TolX.';
        case 3, message='Magnitude of gradient smaller than the specified tolerance';
        case 4, message='Boundary fminimum reached.';
        case 0, message='Number of iterations exceeded options.MaxIter or number of function evaluations exceeded options.FunEvals.';
        case -1, message='Algorithm was terminated by the output function.';
        case -2, message='Line search cannot find an acceptable point along the current search';
        otherwise, message='Undefined exit code';
    end

function stopt=call_output_function(data,optim,where)
stopt=false;
if(~isempty(optim.OutputFcn))
    output.iteration = data.iteration;
    output.funccount = data.funcCount;
    output.fval = data.fInitial;
    output.stepsize = data.alpha;
    output.directionalderivative = data.fPrimeInitial;
    output.gradient = reshape(data.gradient, data.xsizes);
    output.searchdirection = data.dir;
    stopt=feval(optim.OutputFcn,reshape(data.xInitial,data.xsizes),output,where);
end

function data=linesearch_simple(funfcn, data, optim)
% Find a bracket of acceptable points
data = bracketingPhase_simple(funfcn, data, optim);

if (data.bracket_exitflag  == 2)
  % BracketingPhase found a bracket containing acceptable points;
  % now find acceptable point within bracket
  data = sectioningPhase_simple(funfcn, data, optim);
  data.exitflag = data.section_exitflag;
else
  % Already acceptable point found or MaxFunEvals reached
  data.exitflag = data.bracket_exitflag;
end

function data = bracketingPhase_simple(funfcn, data,optim)
% Number of itterations
itw=0;

% Point with smaller value, initial
data.beta=0;
data.f_beta=data.fInitial;
data.fPrime_beta=data.fPrimeInitial;

% Initial step is equal to alpha of previous step.
alpha = data.initialStepLength;

% Going up hill
hill=false;

% Search for brackets
while(true)
    % Calculate the error registration gradient
    if(optim.GradConstr)
        [data,f_alpha]=gradient_function(data.xInitial(:)+alpha*data.dir(:),funfcn, data, optim);
        fPrime_alpha=nan;
        grad=nan;
    else
        [data,f_alpha, grad]=gradient_function(data.xInitial(:)+alpha*data.dir(:),funfcn, data,optim);
        fPrime_alpha = grad'*data.dir(:);
    end

	% Store values linesearch
	data.storefx=[data.storefx f_alpha];
    data.storepx=[data.storepx fPrime_alpha];
	data.storex=[data.storex alpha];
	data.storegx=[data.storegx grad(:)];

    % Update step value
    if(data.f_beta<f_alpha),
        % Go to smaller stepsize
        alpha=alpha*optim.tau3;

        % Set hill variable
        hill=true;
    else
        % Save current minium point
        data.beta=alpha; data.f_beta=f_alpha; data.fPrime_beta=fPrime_alpha; data.grad=grad;
        if(~hill)
            alpha=alpha*optim.tau1;
        end
    end

    % Update number of loop iterations
    itw=itw+1;

    if(itw>(log(optim.TolFun)/log(optim.tau3))),
      % No new optium found, linesearch failed.
      data.bracket_exitflag=-2; break;
    end

    if(data.beta>0&&hill)
            % Get the brackets around minimum point
            % Pick bracket A from stored trials
            [t,i]=sort(data.storex,'ascend');
            storefx=data.storefx(i);storepx=data.storepx(i); storex=data.storex(i);
            [t,i]=find(storex>data.beta,1);
            if(isempty(i)), [t,i]=find(storex==data.beta,1); end
            alpha=storex(i); f_alpha=storefx(i); fPrime_alpha=storepx(i);

            % Pick bracket B from stored trials
            [t,i]=sort(data.storex,'descend');
            storefx=data.storefx(i);storepx=data.storepx(i); storex=data.storex(i);
            [t,i]=find(storex<data.beta,1);
            if(isempty(i)), [t,i]=find(storex==data.beta,1); end
            beta=storex(i); f_beta=storefx(i); fPrime_beta=storepx(i);

            % Calculate derivatives if not already calculated
            if(optim.GradConstr)
                gstep=data.initialStepLength/1e6;
                if(gstep>optim.DiffMaxChange), gstep=optim.DiffMaxChange; end
                if(gstep<optim.DiffMinChange), gstep=optim.DiffMinChange; end
                [data,f_alpha2]=gradient_function(data.xInitial(:)+(alpha+gstep)*data.dir(:),funfcn, data, optim);
                [data,f_beta2]=gradient_function(data.xInitial(:)+(beta+gstep)*data.dir(:),funfcn, data, optim);
                fPrime_alpha=(f_alpha2-f_alpha)/gstep;
                fPrime_beta=(f_beta2-f_beta)/gstep;
            end

            % Set the brackets A and B
            data.a=alpha; data.f_a=f_alpha; data.fPrime_a=fPrime_alpha;
            data.b=beta; data.f_b=f_beta; data.fPrime_b=fPrime_beta;

            % Finished bracketing phase
            data.bracket_exitflag  = 2; return
    end

	% Reached max function evaluations
	if(data.funcCount>=optim.MaxFunEvals), data.bracket_exitflag=0; return; end
end

function data = sectioningPhase_simple(funfcn, data, optim)
% Get the brackets
brcktEndpntA=data.a; brcktEndpntB=data.b;

% Calculate minimum between brackets
[alpha,f_alpha_estimated] = pickAlphaWithinInterval(brcktEndpntA,brcktEndpntB,data.a,data.b,data.f_a,data.fPrime_a,data.f_b,data.fPrime_b,optim);
if(isfield(data,'beta')&&(data.f_beta<f_alpha_estimated)), alpha=data.beta; end

[t,i]=find(data.storex==alpha,1);
if((~isempty(i))&&(~isnan(data.storegx(i))))
    f_alpha=data.storefx(i); grad=data.storegx(:,i);
else
    % Calculate the error and gradient for the next minimizer itteration
    [data,f_alpha, grad]=gradient_function(data.xInitial(:)+alpha*data.dir(:),funfcn, data,optim);
    if(isfield(data,'beta')&&(data.f_beta<f_alpha)),
        alpha=data.beta;
        if((~isempty(i))&&(~isnan(data.storegx(i))))
            f_alpha=data.storefx(i); grad=data.storegx(:,i);
        else
            [data,f_alpha, grad]=gradient_function(data.xInitial(:)+alpha*data.dir(:),funfcn, data,optim);
        end
    end
end

% Store values linesearch
data.storefx=[data.storefx f_alpha]; data.storex=[data.storex alpha];

fPrime_alpha = grad'*data.dir(:);
data.alpha=alpha;
data.fPrime_alpha= fPrime_alpha;
data.f_alpha= f_alpha;
data.grad=grad;

% Set the exit flag to succes
data.section_exitflag=[];

function data=linesearch(funfcn, data, optim)

% Find a bracket of acceptable points
data = bracketingPhase(funfcn, data,optim);

if (data.bracket_exitflag  == 2)
  % BracketingPhase found a bracket containing acceptable points;
  % now find acceptable point within bracket
  data = sectioningPhase(funfcn, data, optim);
  data.exitflag = data.section_exitflag;
else
  % Already acceptable point found or MaxFunEvals reached
  data.exitflag = data.bracket_exitflag;
end

function data = sectioningPhase(funfcn, data, optim)
%
% sectioningPhase finds an acceptable point alpha within a given bracket [a,b]
% containing acceptable points. Notice that funcCount counts the total number of
% function evaluations including those of the bracketing phase.

while(true)

    % Pick alpha in reduced bracket
    brcktEndpntA = data.a + min(optim.tau2,optim.sigma)*(data.b - data.a);
    brcktEndpntB = data.b - optim.tau3*(data.b - data.a);

    % Find global minimizer in bracket [brcktEndpntA,brcktEndpntB] of 3rd-degree
    % polynomial that interpolates f() and f'() at "a" and at "b".
    alpha = pickAlphaWithinInterval(brcktEndpntA,brcktEndpntB,data.a,data.b,data.f_a,data.fPrime_a,data.f_b,data.fPrime_b,optim);

    % No acceptable point could be found
    if (abs( (alpha - data.a)*data.fPrime_a ) <= data.TolFunLnS), data.section_exitflag = -2; return; end

    % Calculate value (and gradient if no extra time cost) of current alpha
    if(~optim.GradConstr)
        [data,f_alpha, grad]=gradient_function(data.xInitial(:)+alpha*data.dir(:),funfcn, data, optim);
        fPrime_alpha = grad'*data.dir(:);
    else
        gstep=data.initialStepLength/1e6;
        if(gstep>optim.DiffMaxChange), gstep=optim.DiffMaxChange; end
        if(gstep<optim.DiffMinChange), gstep=optim.DiffMinChange; end
        [data,f_alpha]=gradient_function(data.xInitial(:)+alpha*data.dir(:),funfcn, data,optim);
        [data,f_alpha2]=gradient_function(data.xInitial(:)+(alpha+gstep)*data.dir(:),funfcn, data, optim);
        fPrime_alpha=(f_alpha2-f_alpha)/gstep;
    end

	% Store values linesearch
	data.storefx=[data.storefx f_alpha]; data.storex=[data.storex alpha];

    % Store current bracket position of A
    aPrev = data.a;
    f_aPrev = data.f_a;
    fPrime_aPrev = data.fPrime_a;

    % Update the current brackets
    if ((f_alpha > data.fInitial + alpha*optim.rho*data.fPrimeInitial) || (f_alpha >= data.f_a))
        % Update bracket B to current alpha
        data.b = alpha; data.f_b = f_alpha; data.fPrime_b = fPrime_alpha;
    else
        % Wolfe conditions, if true then acceptable point found
        if (abs(fPrime_alpha) <= -optim.sigma*data.fPrimeInitial),
            if(optim.GradConstr)
                % Gradient was not yet calculated because of time costs
                [data,f_alpha, grad]=gradient_function(data.xInitial(:)+alpha*data.dir(:),funfcn, data, optim);
                fPrime_alpha = grad'*data.dir(:);
            end
            % Store the found alpha values
            data.alpha=alpha; data.fPrime_alpha= fPrime_alpha; data.f_alpha= f_alpha;
            data.grad=grad;
            data.section_exitflag = []; return,
        end

        % Update bracket A
        data.a = alpha; data.f_a = f_alpha;  data.fPrime_a = fPrime_alpha;

        if (data.b - data.a)*fPrime_alpha >= 0
            % B becomes old bracket A;
            data.b = aPrev; data.f_b = f_aPrev;  data.fPrime_b = fPrime_aPrev;
        end
    end

    % No acceptable point could be found
    if (abs(data.b-data.a) < eps), data.section_exitflag = -2; return, end

    % maxFunEvals reached
    if(data.funcCount >optim.MaxFunEvals), data.section_exitflag = -1; return, end
end

function data = bracketingPhase(funfcn, data, optim)
% bracketingPhase finds a bracket [a,b] that contains acceptable points; a bracket
% is the same as a closed interval, except that a > b is allowed.
%
% The outputs f_a and fPrime_a are the values of the function and the derivative
% evaluated at the bracket endpoint 'a'. Similar notation applies to the endpoint
% 'b'.

% Parameters of bracket A
data.a = [];
data.f_a = [];
data.fPrime_a = [];

% Parameters of bracket B
data.b = [];
data.f_b = [];
data.fPrime_b = [];

% First trial alpha is user-supplied
% f_alpha will contain f(alpha) for all trial points alpha
% fPrime_alpha will contain f'(alpha) for all trial points alpha
alpha = data.initialStepLength;
f_alpha = data.fInitial;
fPrime_alpha = data.fPrimeInitial;

% Set maximum value of alpha (determined by fminimum)
alphaMax = (data.fminimum - data.fInitial)/(optim.rho*data.fPrimeInitial);
alphaPrev = 0;

while(true)
  % Evaluate f(alpha) and f'(alpha)
  fPrev = f_alpha;
  fPrimePrev = fPrime_alpha;

  % Calculate value (and gradient if no extra time cost) of current alpha
  if(~optim.GradConstr)
      [data,f_alpha, grad]=gradient_function(data.xInitial(:)+alpha*data.dir(:),funfcn, data, optim);
      fPrime_alpha = grad'*data.dir(:);
  else
      gstep=data.initialStepLength/1e6;
      if(gstep>optim.DiffMaxChange), gstep=optim.DiffMaxChange; end
      if(gstep<optim.DiffMinChange), gstep=optim.DiffMinChange; end
      [data,f_alpha]=gradient_function(data.xInitial(:)+alpha*data.dir(:),funfcn, data, optim);
      [data,f_alpha2]=gradient_function(data.xInitial(:)+(alpha+gstep)*data.dir(:),funfcn, data, optim);
      fPrime_alpha=(f_alpha2-f_alpha)/gstep;
  end

  % Store values linesearch
  data.storefx=[data.storefx f_alpha]; data.storex=[data.storex alpha];

  % Terminate if f < fminimum
  if (f_alpha <= data.fminimum), data.bracket_exitflag = 4; return; end

  % Bracket located - case 1 (Wolfe conditions)
  if (f_alpha > (data.fInitial + alpha*optim.rho*data.fPrimeInitial)) || (f_alpha >= fPrev)
    % Set the bracket values
    data.a = alphaPrev; data.f_a = fPrev;  data.fPrime_a = fPrimePrev;
    data.b = alpha; data.f_b = f_alpha;  data.fPrime_b = fPrime_alpha;
    % Finished bracketing phase
    data.bracket_exitflag  = 2; return
  end

  % Acceptable steplength found
  if (abs(fPrime_alpha) <= -optim.sigma*data.fPrimeInitial),
      if(optim.GradConstr)
          % Gradient was not yet calculated because of time costs
          [data,f_alpha, grad]=gradient_function(data.xInitial(:)+alpha*data.dir(:),funfcn, data, optim);
          fPrime_alpha = grad'*data.dir(:);
      end
      % Store the found alpha values
      data.alpha=alpha;
      data.fPrime_alpha= fPrime_alpha; data.f_alpha= f_alpha; data.grad=grad;
      % Finished bracketing phase, and no need to call sectioning phase
      data.bracket_exitflag = [];  return
  end

  % Bracket located - case 2
  if (fPrime_alpha >= 0)
    % Set the bracket values
    data.a = alpha; data.f_a = f_alpha;  data.fPrime_a = fPrime_alpha;
    data.b = alphaPrev; data.f_b = fPrev; data.fPrime_b = fPrimePrev;
    % Finished bracketing phase
    data.bracket_exitflag  = 2; return
  end

  % Update alpha
  if (2*alpha - alphaPrev < alphaMax )
      brcktEndpntA = 2*alpha-alphaPrev;
      brcktEndpntB = min(alphaMax,alpha+optim.tau1*(alpha-alphaPrev));
      % Find global minimizer in bracket [brcktEndpntA,brcktEndpntB] of 3rd-degree polynomial
      % that interpolates f() and f'() at alphaPrev and at alpha
      alphaNew = pickAlphaWithinInterval(brcktEndpntA,brcktEndpntB,alphaPrev,alpha,fPrev, ...
                                         fPrimePrev,f_alpha,fPrime_alpha,optim);
      alphaPrev = alpha;
      alpha = alphaNew;
  else
      alpha = alphaMax;
  end

  % maxFunEvals reached
  if(data.funcCount >optim.MaxFunEvals), data.bracket_exitflag = -1; return, end
end

function [alpha,f_alpha]= pickAlphaWithinInterval(brcktEndpntA,brcktEndpntB,alpha1,alpha2,f1,fPrime1,f2,fPrime2,optim)
% finds a global minimizer alpha within the bracket [brcktEndpntA,brcktEndpntB] of the cubic polynomial
% that interpolates f() and f'() at alpha1 and alpha2. Here f(alpha1) = f1, f'(alpha1) = fPrime1,
% f(alpha2) = f2, f'(alpha2) = fPrime2.

% determines the coefficients of the cubic polynomial with c(alpha1) = f1,
% c'(alpha1) = fPrime1, c(alpha2) = f2, c'(alpha2) = fPrime2.
coeff = [(fPrime1+fPrime2)*(alpha2-alpha1)-2*(f2-f1) ...
    3*(f2-f1)-(2*fPrime1+fPrime2)*(alpha2-alpha1) (alpha2-alpha1)*fPrime1 f1];

% Convert bounds to the z-space
lowerBound = (brcktEndpntA - alpha1)/(alpha2 - alpha1);
upperBound = (brcktEndpntB - alpha1)/(alpha2 - alpha1);

% Swap if lowerbound is higher than the upperbound
if (lowerBound  > upperBound), t=upperBound; upperBound=lowerBound; lowerBound=t; end

% Find minima and maxima from the roots of the derivative of the polynomial.
sPoints = roots([3*coeff(1) 2*coeff(2) coeff(3)]);

% Remove imaginaire and points outside range

sPoints(imag(sPoints)~=0)=[];
sPoints(sPoints<lowerBound)=[]; sPoints(sPoints>upperBound)=[];

% Make vector with all possible solutions
sPoints=[lowerBound sPoints(:)' upperBound];

% Select the global minimum point
[f_alpha,index]=min(polyval(coeff,sPoints)); z=sPoints(index);

% Add the offset and scale back from [0..1] to the alpha domain
alpha = alpha1 + z*(alpha2 - alpha1);

% Show polynomial search
if(optim.Display(1)=='p');
    vPoints=polyval(coeff,sPoints);
    plot(sPoints*(alpha2 - alpha1)+alpha1,vPoints,'co');
    plot([sPoints(1) sPoints(end)]*(alpha2 - alpha1)+alpha1,[vPoints(1) vPoints(end)],'c*');
    xPoints=linspace(lowerBound/3, upperBound*1.3, 50);
    vPoints=polyval(coeff,xPoints);
    plot(xPoints*(alpha2 - alpha1)+alpha1,vPoints,'c');
end

function [data,fval,grad]=gradient_function(x,funfcn, data, optim)
    % Call the error function for error (and gradient)
    if ( nargout <3 )
        %%% timem=tic;
        fval=funfcn(reshape(x,data.xsizes));
        %%% data.timeExtern=data.timeExtern+toc(timem);
        data.funcCount=data.funcCount+1;
    else
        if(strcmp(optim.GradObj,'on'))
            %%% timem=tic;
            [fval, grad]=feval(funfcn,reshape(x,data.xsizes));
            %%% data.timeExtern=data.timeExtern+toc(timem);
            data.funcCount=data.funcCount+1;
            data.gradCount=data.gradCount+1;
        else
            % Calculate gradient with forward difference if not provided by the function
            grad=zeros(length(x),1);
            fval=funfcn(reshape(x,data.xsizes));
            gstep=data.initialStepLength/1e6;
            if(gstep>optim.DiffMaxChange), gstep=optim.DiffMaxChange; end
            if(gstep<optim.DiffMinChange), gstep=optim.DiffMinChange; end
            for i=1:length(x),
                x_temp=x; x_temp(i)=x_temp(i)+gstep;
                %%% timem=tic;
                [fval_g]=feval(funfcn,reshape(x_temp,data.xsizes)); data.funcCount=data.funcCount+1;
                %%% data.timeExtern=data.timeExtern+toc(timem);
                grad(i)=(fval_g-fval)/gstep;
            end
        end
        grad=grad(:);
    end

function data = updateQuasiNewtonMatrix_LBFGS(data,optim)
% updates the quasi-Newton matrix that approximates the inverse to the Hessian.
% Two methods are support BFGS and L-BFGS, in L-BFGS the hessian is not
% constructed or stored.
% Calculate position, and gradient diference between the
% itterations
deltaX=data.alpha* data.dir;
deltaG=data.gradient-data.gOld;

if ((deltaX'*deltaG) >= sqrt(eps)*max( eps,norm(deltaX)*norm(deltaG) ))

    if(optim.HessUpdate(1)=='b')
        % Default BFGS as described by Nocedal
        p_k = 1 / (deltaG'*deltaX);
        Vk = eye(data.numberOfVariables) - p_k*deltaG*deltaX';
        % Set Hessian
        data.Hessian = Vk'*data.Hessian *Vk + p_k * deltaX*deltaX';
        % Set new Direction
        data.dir = -data.Hessian*data.gradient;
    else
        % L-BFGS with scaling as described by Nocedal

        % Update a list with the history of deltaX and deltaG
        data.deltaX(:,2:optim.StoreN)=data.deltaX(:,1:optim.StoreN-1); data.deltaX(:,1)=deltaX;
        data.deltaG(:,2:optim.StoreN)=data.deltaG(:,1:optim.StoreN-1); data.deltaG(:,1)=deltaG;

        data.nStored=data.nStored+1; if(data.nStored>optim.StoreN), data.nStored=optim.StoreN; end

        % Initialize variables
        a=zeros(1,data.nStored);
        p=zeros(1,data.nStored);

        q = data.gradient;
        for i=1:data.nStored
            p(i)= 1 / (data.deltaG(:,i)'*data.deltaX(:,i));
            a(i) = p(i)* data.deltaX(:,i)' * q;
            q = q - a(i) * data.deltaG(:,i);
        end
        % Scaling of initial Hessian (identity matrix)
        p_k = data.deltaG(:,1)'*data.deltaX(:,1) / sum(data.deltaG(:,1).^2);

        % Make r = - Hessian * gradient
        r = p_k * q;
        for i=data.nStored:-1:1,
            b = p(i) * data.deltaG(:,i)' * r;
            r = r + data.deltaX(:,i)*(a(i)-b);
        end

        % Set new direction
        data.dir = -r;
    end
end

function earlystoppedreturn = LVQ_progresser(variables, output, where)
%
% control actions during each iteration of the fminlbfgs optimization
%
% Conditions of GNU General Public License, version 2 and BSD License apply.
% See file 'license-gpl2.txt' and 'BSD_license.txt' enclosed in this package.
% Programs are not for use in critical applications!
% Kerstin Bunte based on the code from Marc Strickert
%
global earlystopped threshstop

if sum(where == 'iter') < 4
    earlystopped = false;
    earlystoppedreturn = earlystopped;
    return
end

earlystopped = exist ('stop','file') > 0;
earlystoppedreturn = earlystopped;

if earlystopped
    delete('stop');
    return
end

if isempty(threshstop)
    return
end

if threshstop > 0  % early stopping by simple thresholding: if too good: stop
    earlystopped = output.fval < threshstop;  % close to perfect
    earlystoppedreturn = earlystopped;
    if earlystopped
      return
    end
end

global useEarlyStopping

if isempty(useEarlyStopping)
    return
end

%%% early stopping

global training_data training_label prototypeLabel% n_vec bestvariables
nb_prototypes = numel(prototypeLabel);
model.w = variables(1:nb_prototypes,:);
model.c_w = prototypeLabel;
model.omega = variables(nb_prototypes+1:end,:);
estimatedLabels = GMLVQ_classify(training_data, model);
ref = mean( training_label ~= estimatedLabels );
% fprintf('error: %f',ref);
%   if isempty(n_vec) || n_vec < 1 % GRLVQ
%     lam = variables(1,:); % weights
%     protos = variables(2:end,:); % protos
%     [cls C] = applyprotoseuc(protos, protolabl, lam, datval, labval);
%   else % MRLVQ
%     lam = variables(1:n_vec,:).';
%     protos = variables((n_vec+1):end,:);
%     [cls C] = applyprotosmat(protos, protolabl, lam, datval, labval);
%   end
% ref = trace(C)/sum(C(:))-1; % negative classif error

persistent last penalty
% disp(last);
if isempty(last)
%     bestvariables = variables;
    penalty = -25;  % penalty counter, skip first 15 iterations (burn-in phase)
%     last = ref * ones(5,1);  % memory of 5 elements
    last = ref;
else
%     [srt idxsrt] = sort([last; ref],'descend');
%     fvalval = srt(1);
%     fpos = find(idxsrt == 6);  % index of added element
%     if fpos == 6  % last position?
      penalty = penalty + 1;
      if penalty > 10  % allow a maximum of 10 times to fail
        earlystopped = true;
        earlystoppedreturn = earlystopped;
        return
      end
%     else
% %       last = srt(1:5);
% %       if fpos == 1  % best position
% %         bestvariables = variables;
% %       end
%       if penalty < 0  % still in burn-in phase
%         penalty = penalty + 1;
%       end
%       if penalty > 0  % do not go below here
%         penalty = penalty - 1;
%       end
%     end
end
earlystoppedreturn = earlystopped;

author: ilarinieminen

repo: SOM-Toolbox


text:
function score = BIC_score_CPD(CPD, fam, data, ns, cnodes)
% BIC_score_CPD Compute the BIC score of a tabular CPD
% score = BIC_score_CPD(CPD, fam, data, ns, cnodes)

if iscell(data)
  local_data = cell2num(data(fam,:));
else
  local_data = data(fam, :);
end
counts = compute_counts(local_data, CPD.sizes);
CPT = mk_stochastic(counts); % MLE
tiny = exp(-700);
CPT = CPT + (CPT==0)*tiny;  % replace 0s by tiny
LL = sum(log(CPT(:)) .* counts(:));
N = size(data, 2);
score = LL - 0.5*CPD.nparams*log(N);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
function boxes = labels2boxes(superpixels, labels)

% Get boxes of the superpixels
tmp = regionprops(superpixels,'BoundingBox');
tmp = cat(1,tmp.BoundingBox);
sp_bboxes = [tmp(:,2)+0.5, tmp(:,1)+0.5, tmp(:,2)+tmp(:,4)-0.5, tmp(:,1)+tmp(:,3)-0.5];

% Compute the boxes of the labels from those of the superpixels
boxes = zeros(length(labels),4);
for ii=1:length(labels)
    curr_sp_bboxes = sp_bboxes(labels{ii},:);
    boxes(ii,1) = min(curr_sp_bboxes(:,1));
    boxes(ii,2) = min(curr_sp_bboxes(:,2));
    boxes(ii,3) = max(curr_sp_bboxes(:,3));
    boxes(ii,4) = max(curr_sp_bboxes(:,4));
end

author: jponttuset

repo: mcg

save_path: github-repos/MATLAB/jponttuset-mcg


text:
% eeg_point2lat() - convert latency in data points to latency in ms relative
%                   to the time locking. Used in eeglab().
% Usage:
%       >> [newlat ] = eeg_point2lat( lat_array, [], srate);
%       >> [newlat ] = eeg_point2lat( lat_array, epoch_array,...
%                                 srate, timelimits, timeunit);
% Inputs:
%   lat_array   - latency array in data points assuming concatenated
%                 data epochs (see eeglab() event structure)
%   epoch_array - epoch number corresponding to each latency value
%   srate       - data sampling rate in Hz
%   timelimits  - [min max] timelimits in 'timeunit' units (see below)
%   timeunit    - time unit in second. Default is 1 = seconds.
%
% Outputs:
%   newlat      - converted latency values (in 'timeunit' units) for each epoch
%
% Example:
%   tmpevent = EEG.event;
%   eeg_point2lat( [ tmpevent.latency ], [], EEG.srate);
%   % returns the latency of all events in second for a continuous
%   % dataset EEG
%
%   eeg_point2lat( [ tmpevent.latency ], [ tmpevent.epoch ],
%                 EEG.srate, [EEG.xmin EEG.xmax]*1000, 1E-3);
%   % returns the latency of all events in millisecond for a dataset
%   % containing data epochs.
%
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2 Mai 2002
%
% See also: eeg_lat2point(), eeglab(), pop_editieventvals(), pop_loaddat()

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

function newlat = eeg_point2lat( lat_array, epoch_array, srate, timewin, timeunit);

if nargin <3
    help eeg_point2lat;
    return;
end;
if isempty( epoch_array )
    epoch_array = ones( size(lat_array) );
end;
if nargin <4
    timewin = 0;
end;
if nargin <5
	timeunit = 1;
end;

if length(lat_array) ~= length(epoch_array)
	if length(epoch_array)~= 1
		disp('eeg_point2lat: latency and epoch arrays must have the same length'); return;
	else
		epoch_array = ones(1,length(lat_array))*epoch_array;
	end;
end;
if length(timewin) ~= 2
    disp('eeg_point2lat: timelimits array must have length 2'); return;
end;
if iscell(epoch_array)
	epoch_array = [ epoch_array{:} ];
end;
if iscell(lat_array)
	lat_array = [ lat_array{:} ];
end

timewin = timewin*timeunit;

if length(timewin) == 2
    pnts = (timewin(2)-timewin(1))*srate+1;
else
    pnts = 0;
end;
newlat  = ((lat_array - (epoch_array-1)*pnts-1)/srate+timewin(1))/timeunit;
newlat = round(newlat*1E9)*1E-9;

author: PatternRecognition

repo: OpenBMI


text:
function sys = setdualize(sys,state)

sys.dualized = state;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%% uniform grid
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],1);
showmesh3(node,elem);
findnode3(node,'all');
[node,hexelem] = tet2hex(node,elem,HB);

%% Adaptive grid
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],1/4);
[node,elem] = delmesh(node,elem,'x<0 & y<0');
bdFlag = setboundary3(node,elem,'Neumann');
for k = 1:2
    eta = abs(sign(f(node(elem(:,1),:))) + sign(f(node(elem(:,2),:)))...
            + sign(f(node(elem(:,3),:))) + sign(f(node(elem(:,4),:))));
    refineElem = find(eta < 4);
    [node,elem,bdFlag,HB] = bisect3(node,elem,refineElem,bdFlag,HB);
end
[tempvar,bdFace] = findboundary3(elem,bdFlag);
showmesh(node,bdFace);
[node,hexelem] = tet2hex(node,elem,HB);

%% sphere
function s = f(p)
    s = sum(p.^2,2) - (0.5)^2;
end

repo: ifem


text:
function X=solveSylvesterEq(A,B,C)
%%SOLVESYLVESTEREQ Solve the Sylester equation A*X+X*B=C for X, where A, B,
%           and C can be real or complex. This implementation only works
%           when a unique solution for X exists, which is the case when
%           none of the eigenvalues of A and B when summed together equals
%           0.
%
%INPUTS: A A pXp matrix.
%        B An rXr matrix.
%        C A pXr matrix.
%
%OUTPUTS: X The pXr matrix solving A*X+X*B=C.
%
%This function implements the Bartels-Stewart algorithm of [1]. The back
%substitution algorithm of Algorithm 7.6.2 in [2] is used for the second
%half of the algorithm, which also means that upper-triangular Schur
%decompositions are used for the first half of the algorithm.
%
%EXAMPLE:
%In this example, we solve a Sylvester equation and show that the result
%Tolerance is about 7e-14, which is good agreement.
% A=[3,  0,   0, 23;
%   -2,  5, -14, -6;
%   11, 11,  -7,  7;
%  -10, 15, -10, -1];
% B=[8,  4,   2,  6;
%   -7, -1,   1,  8;
%  -14, -1,  15, -2;
%  -14, 14,  -8,  2];
% C=[16,  3,  6,  8;
%     3, 11,  8, 11;
%     6,  8,  6, 13;
%     8, 11, 13,  1];
% X=solveSylvesterEq(A,B,C);
% AbsTol=max(max(abs(A*X+X*B-C)))
%
%REFERENCES
%[1] R. H. Bartels and G. W. Stewart, "Algorithm 432: Solution of the
%    matrix equation AX+XB=C [F4]," Communications of the ACM, pp. 820-826,
%    Sep. 1972.
%[2] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed.
%    Baltimore: The Johns Hopkins Press, 2013.
%
%December 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%We use the complex Schur decomposition. If we used the real one, then
%Ap and Bp are not guaranteed to be upper triangular if the matrices have
%complex eigenvalues.
[U,Ap]=schur(A,'complex');
[V,Bp]=schur(B,'complex');
Cp=U'*C*V;

%Solve the transformed system using Algorithm 6.7.2 of [1].
Bp=-Bp;

r=size(Ap,1);
p=size(Bp,1);

I=eye(r,r);

for k=1:r
    Cp(1:p,k)=Cp(1:p,k)+Cp(1:p,1:(k-1))*Bp(1:(k-1),k);
    z=(Ap-Bp(k,k)*I)\Cp(1:p,k);
    Cp(1:p,k)=z;
end

%Cp hold the transformed result. Now undo the transformation.
X=U*Cp*V';

if(isreal(A)&&isreal(B)&&isreal(C))
    %If all of the inputs are real, then Cp should be real, not accounting
    %for finite precision limits.
    X=real(X);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function audnoise(ns_file,outfile)

%
%  Implements the audible-noise suppression algorithm [1].
%
%  Usage:  audnoise(noisyFile, outputFile)
%
%         infile - noisy speech file in .wav format
%         outputFile - enhanced output file in .wav format
%
%  It runs 2 iterations, but one could change the number of iterations by
%  modifying accordingly the variable iter_num on line 33.
%
%  Example call:  audnoise('sp04_babble_sn10.wav','out_aud.wav');
%
%  References:
%   [1] Tsoukalas, D. E., Mourjopoulos, J. N., and Kokkinakis, G. (1997). Speech
%       enhancement based on audible noise suppression. IEEE Trans. on Speech and
%       Audio Processing, 5(6), 497-514.
%
% Authors: Yi Hu and Philipos C. Loizou
%
% Copyright (c) 2006 by Philipos C. Loizou
% $Revision: 0.0 $  $Date: 10/09/2006 $
%-------------------------------------------------------------------------

if nargin<2
   fprintf('Usage: audnoise(noisyfile.wav,outFile.wav) \n\n');
   return;
end

iter_num=2;  % number of iterations
NF_SABSENT= 6;
%this is the number of speech-absent frames to estimate the initial
%noise power spectrum

[nsdata, Fs, bits]= wavread( ns_file);	%nsdata is a column vector

aa=0.98;
mu=0.98;
eta=0.15;

nwind= floor( 20* Fs/ 1000);	%this corresponds to 20ms window
if rem( nwind, 2)~= 0 nwind= nwind+ 1; end	%made window length even
noverlap= nwind/ 2;
w= hamming( nwind);
rowindex= ( 1: nwind)';

%we assume the first NF_SABSENT frames are speech absent, we use them to estimate the noise power spectrum
noisedata= nsdata( 1: nwind* NF_SABSENT);	noise_colindex= 1+ ( 0: NF_SABSENT- 1)* nwind;
noisematrixdata = zeros( nwind, NF_SABSENT);
noisematrixdata( :)= noisedata( ...
    rowindex( :, ones(1, NF_SABSENT))+ noise_colindex( ones( nwind, 1), :)- 1);
noisematrixdata= noisematrixdata.* w( :, ones( 1, NF_SABSENT)) ;	%WINDOWING NOISE DATA
noise_ps= mean( (abs( fft( noisematrixdata))).^ 2, 2); %NOTE!!!! it is a column vector

% ----- estimate noise in CBs ------------------
%
noise_b=zeros(nwind/2+1,1);
[CB_FREQ_INDICES]=find_CB_FREQ_INDICES(Fs,nwind,16,nwind/2);

for i = 1:length(CB_FREQ_INDICES)
    noise_b(CB_FREQ_INDICES{i})=ones(size(CB_FREQ_INDICES{i},2),1)*mean(noise_ps(CB_FREQ_INDICES{i}));
end
noise_b1=[noise_b; fliplr(noise_b(2:nwind/2))];

nslide= nwind- noverlap;

x= nsdata;
nx= length( x);	ncol= fix(( nx- noverlap)/ nslide);
colindex = 1 + (0: (ncol- 1))* nslide;
if nx< (nwind + colindex(ncol) - 1)
    x(nx+ 1: nwind+ colindex(ncol) - 1) = ...
        rand( nwind+ colindex( ncol)- 1- nx, 1)* (2^ (-15));   % zero-padding
end

es_old= zeros( noverlap, 1);
%es_old is actually the second half of the previous enhanced speech frame,
%it is used for overlap-add

for k= 1: ncol

    y= x( colindex( k): colindex( k)+ nwind- 1);
    y= y.* w;	%WINDOWING NOISY SPEECH DATA

    y_spec= fft( y);	y_specmag= abs( y_spec);	y_specang= angle( y_spec);
    %they are the frequency spectrum, spectrum magnitude and spectrum phase, respectively

    y_ps= y_specmag.^ 2;	%power spectrum of noisy speech
    y_ps1=y_ps(1:nwind/2+1);

    % ====start of vad ===
    gammak=min(y_ps./noise_ps,40);  % post SNR
    if k==1
        ksi=aa+(1-aa)*max(gammak-1,0);
    else
        ksi=aa*Xk_prev./noise_ps + (1-aa)*max(gammak-1,0);     % a priori SNR
    end

    log_sigma_k= gammak.* ksi./ (1+ ksi)- log(1+ ksi);
    vad_decision= sum( log_sigma_k)/ nwind;
    if (vad_decision < eta)
        % noise only frame found
        noise_ps= mu* noise_ps+ (1- mu)* y_ps;
    end

    for i = 1:length(CB_FREQ_INDICES)
        noise_b(CB_FREQ_INDICES{i})=...
            ones(size(CB_FREQ_INDICES{i},2),1)*mean(noise_ps(CB_FREQ_INDICES{i}));
    end

    % ===end of vad===

    x_cons1=max(y_ps-noise_ps,0.001);
    % conservative estimate of x from power spectral subtraction
    x_cons = x_cons1(1:nwind/2+1);

    % --- Estimate masking thresholds iteratively (as per page 505) ----
    %
    Tk0=mask(x_cons,nwind,Fs,16);
    Xp=y_ps1;
    for j=1:iter_num
        ab = noise_b+(noise_b.^2)./Tk0;  % Eq. 41
        Xp=(Xp.^2)./(ab+Xp);             % Eq. 40
        Tk0=mask(Xp,nwind,Fs,16);
    end

    % --- Estimate alpha ------
    %
    alpha = (noise_b+Tk0).*(noise_b./Tk0);
    % eq. 26 for Threshold (T) method with ni(b)=1

    % ---- Apply suppression rule --------------
    %
    H0 = (Xp./(alpha+Xp));
    H=[H0(1:nwind/2+1); flipud(H0(2:nwind/2))];

    x_hat = H.*y_spec;
    Xk_prev= abs( x_hat).^ 2;

    es_tmp=real(ifft(x_hat));

    % ---- Overlap and add ---------------

    es_data( colindex( k): colindex( k)+ nwind- 1)= [es_tmp( 1: noverlap)+ es_old;...
        es_tmp( noverlap+ 1: nwind)];
    %overlap-add
    es_old= es_tmp( nwind- noverlap+ 1: nwind);
end

wavwrite( es_data, Fs, bits, outfile);

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

function [CB_FREQ_INDICES]=find_CB_FREQ_INDICES(Fs,dft_length,nbits,frame_overlap)
% This function is from Matlab STSA Toolbox for Audio Signal Noise Reduction
% Copyright (C) 2001  Patrick J. Wolfe

freq_val = (0:Fs/dft_length:Fs/2)';
freq=freq_val;
crit_band_ends = [0;100;200;300;400;510;630;770;920;1080;1270;1480;1720;2000;2320;2700;3150;3700;4400;5300;6400;7700;9500;12000;15500;Inf];
imax = max(find(crit_band_ends < freq(end)));
num_bins = length(freq);
LIN_TO_BARK = zeros(imax,num_bins);
i = 1;
for j = 1:num_bins
    while ~((freq(j) >= crit_band_ends(i)) & (freq(j) < crit_band_ends(i+1))),i = i+1;end
    LIN_TO_BARK(i,j) = 1;
end
% Calculation of critical band frequency indices--i.e., which bins are in which critical band for i = 1:imax
for i=1:imax,
    CB_FREQ_INDICES{i} = find(LIN_TO_BARK(i,:));
end

author: jtkim-kaist

repo: Speech-enhancement


text:
function [] = tapas_display_table(summaries, figname)
%% Display a figure with the summaries.
%
% Input
%       summaries       -- Table with the summaries.
%       figname         -- Optional. String with the name for the figure.
%
% Output
%

% aponteeduardo@gmail.com
% copyright (C) 2019
%

n = 1;

n = n + 1;
if nargin < n
    figname = 'Table';
end

figure('name', figname);
uitable('Data', summaries{:, :}, ...
    'ColumnName', summaries.Properties.VariableNames,...
    'RowName', summaries.Properties.RowNames, ...
    'Units', 'Normalized', ...
    'Position',[0, 0, 1, 1]);

end

author: translationalneuromodeling

repo: tapas


text:
% Firmament Autopilot
% MLog (Mat Log) Parser
% @brief: parse mlog.bin file
function [LogHeader, LogMsg] = mlog_parser(logfile)

%% BLog Data Type Definition
TYPE_CAST = ["int8=>int8", "uint8=>uint8", "int16=>int16", "uint16=>uint16",...
    "int32=>int32", "uint32=>uint32", "float=>float", "double=>double"];

% mlog message start/end flag
MLOG_BEGIN_MSG1 = hex2dec('92');
MLOG_BEGIN_MSG2 = hex2dec('05');
MLOG_END_MSG = hex2dec('26');

%% Open Log File
fileID = fopen(logfile, 'r');
fileDir = dir(logfile);
if ~fileID
    disp('log file open fail');
    return;
end

%% Read Log Header
% Read Log Information
LogHeader.version = fread(fileID, 1, 'uint16=>uint16');
LogHeader.timestamp = fread(fileID, 1, 'uint32=>uint32');
LogHeader.max_name_len = fread(fileID, 1, 'uint16=>uint16');
LogHeader.max_desc_len = fread(fileID, 1, 'uint16=>uint16');
LogHeader.max_model_info_len = fread(fileID, 1, 'uint16=>uint16');
LogHeader.description = fread(fileID, LogHeader.max_desc_len, 'char=>char');
LogHeader.model_info = fread(fileID, LogHeader.max_model_info_len, 'char=>char');

fprintf('MLog File:%s\n', logfile);
fprintf('Version:%d\n', LogHeader.version);
fprintf('Timestamp:%d(ms)\n', LogHeader.timestamp);
fprintf('Description:%s\n', LogHeader.description);
fprintf('Model Information:\n%s\n', LogHeader.model_info);

% Read bus number
LogHeader.num_bus = fread(fileID, 1, 'uint8=>uint8');
% Preallocate memory
MsgCount = zeros(1, LogHeader.num_bus);
LogMsg = cell(1, LogHeader.num_bus);
% Read bus data
for n = 1:LogHeader.num_bus
    LogHeader.bus(n).name = fread(fileID, [1 LogHeader.max_name_len], 'uint8=>char');
    LogHeader.bus(n).msg_id = fread(fileID, 1, 'uint8=>uint8');
    LogHeader.bus(n).num_elem = fread(fileID, 1, 'uint8=>uint8');
    % Read bus elements
    for k = 1:LogHeader.bus(n).num_elem
        LogHeader.bus(n).elem_list(k).name = fread(fileID, [1 LogHeader.max_name_len], 'uint8=>char');
        LogHeader.bus(n).elem_list(k).type = fread(fileID, 1, 'uint16=>uint16');
        LogHeader.bus(n).elem_list(k).number = fread(fileID, 1, 'uint16=>uint16');
    end
    % Preallocate memory
    LogMsg{n} = cell(1, LogHeader.bus(n).num_elem);
end

% Read Parameter
LogHeader.num_param_group = fread(fileID, 1, 'uint8=>uint8');

for n = 1:LogHeader.num_param_group
    LogHeader.param_group_list(n).name = fread(fileID, [1 LogHeader.max_name_len], 'uint8=>char');
    LogHeader.param_group_list(n).num_param = fread(fileID, 1, 'uint32=>uint32');
    LogHeader.param_group_list(n).name = LogHeader.param_group_list(n).name(~isspace(LogHeader.param_group_list(n).name));

    for k = 1:LogHeader.param_group_list(n).num_param
        LogHeader.param_group_list(n).param(k).name = fread(fileID, [1 LogHeader.max_name_len], 'uint8=>char');
        LogHeader.param_group_list(n).param(k).type = fread(fileID, 1, 'uint8=>uint8');
        index = LogHeader.param_group_list(n).param(k).type+1;
        LogHeader.param_group_list(n).param(k).val = fread(fileID, 1, TYPE_CAST(index));
        LogHeader.param_group_list(n).param(k).name = LogHeader.param_group_list(n).param(k).name(~isspace(LogHeader.param_group_list(n).param(k).name));
    end
end

%% Read Log Msg
while ~feof(fileID) && ftell(fileID)<fileDir.bytes

    %%% Read Msg Begin Flag %%%
    begin_state = 0;
    while begin_state < 2
        switch begin_state
            case 0
                % search for MLOG_BEGIN_MSG1
                msg_begin = fread(fileID, 1, 'uint8=>uint8');
                if msg_begin == MLOG_BEGIN_MSG1
                    begin_state = 1;
                end
            case 1
                % search for MLOG_BEGIN_MSG2
                msg_begin = fread(fileID, 1, 'uint8=>uint8');
                if msg_begin == MLOG_BEGIN_MSG2
                    begin_state = 2;
                elseif msg_begin == MLOG_BEGIN_MSG1
                    begin_state = 1;
                else
                    begin_state = 0;
                end
            otherwise
                break;
        end
    end

    %%% Read Msg ID %%%
    msg_id = fread(fileID, 1, 'uint8=>uint8');

    % msg_id start from 0 and index start from 1
    index = msg_id + 1;
    if index > LogHeader.num_bus
        index = -1; % invalid msg id
    end

    if index < 0
        fprintf('invalid msg id:%d\r\n', msg_id);
        % continue to search next msg
        continue;
    end

    %%% Read Msg Payload %%%

    % read bus elements
    for k = 1:LogHeader.bus(index).num_elem
        type = LogHeader.bus(index).elem_list(k).type+1;
        len = LogHeader.bus(index).elem_list(k).number;

        [elem_val, rb] = fread(fileID, [len, 1], TYPE_CAST(type));

        if rb < len
            fprintf('%s %s cnt %d read err, delete it\n', LogHeader.bus(index).name, LogHeader.bus(index).elem_list(k).name, MsgCount(index));
            % TODO: handle this error
            break;
        else
            % LogMsg{msg_id}{elem_index}(len:cnt)
            LogMsg{index}{k}(1:len, MsgCount(index)+1) = elem_val;
        end
    end

    %%% Read Msg End Flag %%%
    msg_end = fread(fileID, 1, 'uint8=>uint8');
    if msg_end == MLOG_END_MSG
        % valid msg received
        MsgCount(index) = MsgCount(index) + 1;
    else
        fprintf('invalid msg end flag:%d, msg id:%d\r\n', msg_end, msg_id);
        % delete invalid msg
        for k = 1:LogHeader.bus(index).num_elem
            try
                LogMsg{index}{k}(:, MsgCount(index)+1) = [];
            catch
                continue
            end
        end
    end
end
fclose(fileID);

%% Output Parsed Information
fprintf('\nLog parse finish!\n');
% fprintf('log elapsed time: %f(s)\n', )
for n = 1:LogHeader.num_bus
    BusName = strrep(LogHeader.bus(n).name, '"', '');
    BusName = BusName(~isspace(BusName));
	fprintf('%s: %d msg recorded\n', BusName, MsgCount(n));
end

%% Generate .MAT file
generate_MAT(LogHeader, LogMsg, fileDir.folder);
end

function generate_MAT(LogHeader, LogMsg, Path)
%% create Bus variable and save into .mat file
save(strcat(Path, '/LogHeader.mat'), 'LogHeader'); % save log header

for n = 1:LogHeader.num_bus
    msg_id = LogHeader.bus(n).msg_id;
    index = msg_id + 1;
    if isempty(LogMsg{index})
       continue;
    end

    BusName = strrep(LogHeader.bus(n).name, '"', '');
    BusName = BusName(~isspace(BusName));

    % find timestamp
    timestamp_id = 0;
    for k = 1:LogHeader.bus(n).num_elem
        ElemName = strrep(LogHeader.bus(n).elem_list(k).name, '"', '');
        ElemName = ElemName(~isspace(ElemName));
        if strcmp(ElemName, "timestamp_ms") || strcmp(ElemName, "timestamp")
            timestamp_id = k;
        end
    end

    if timestamp_id <= 0
       fprintf("can't find timestamp element in %s\n", LogHeader.bus(n).name);
       continue;
    else
       time_stamp = double(LogMsg{index}{timestamp_id}-LogMsg{index}{timestamp_id}(1)) * 0.001;   % milli second to second
    end

    % construct Bus variable
    for k = 1:LogHeader.bus(n).num_elem
        ElemName = strrep(LogHeader.bus(n).elem_list(k).name, '"', '');
        ElemName = ElemName(~isspace(ElemName));

        exp = sprintf('timeseries(LogMsg{index}{k}'', time_stamp);');
        if ~isempty(ElemName)
            eval([BusName, '.', ElemName, '=', exp]);
        else
            eval([BusName, '=', exp]);
        end

    end

    % save as .mat file
    out_file = strcat(Path, ['/', BusName, '.mat']);
    save(out_file, BusName);
    fprintf("Save to path:%s\n", out_file);
end
end

author: Firmament-Autopilot

repo: FMT-Model


text:
function im = color(input)
% Convert input image to color.
%   im = color(input)

if size(input, 3) == 1
  im(:,:,1) = input;
  im(:,:,2) = input;
  im(:,:,3) = input;
elseif size(input, 3) == 4 % CMYK
  h = size(input,1);
  w = size(input,2);
  im = reshape(input,[h*w,4]);
  im = rgb2cmyk(im);
  im = reshape(im,[h,w,3]);
else
  im = input;
end

end

author: aayushbansal

repo: PixelNet


text:
function str = showmydoc(item, indent)

% function str = showmydoc(item, indent)
% Display help text for a cfg_choice and all of its options, without
% recursive calls to child nodes.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: showmydoc.m 2085 2008-09-12 10:26:59Z volkmar $

rev = '$Rev: 2085 $'; %#ok
str1 = showdoc(item.cfg_item, indent);
citems = subsref(item, substruct('.','values'));
str2{1} = 'One of the following options must be selected:';
% Display short listing of cfg_choice value items first
str3 = cell(1,numel(citems));
for k = 1:numel(citems)
    str3{k} = sprintf('* %s', subsref(citems{k}, substruct('.','name')));
end;
valitem = subsref(item, substruct('.','val'));
if ~isempty(valitem)
    str4 = {'Currently selected option:', ...
            sprintf('* "%s"',subsref(valitem{1}, substruct('.','name')))};
else
    str4 = {};
end;
str = [str1(:); {''}; str2(:); str3(:); str4(:)];

author: spm


text:
%%******************************************************************
%% ops:
%%
%%   Z = ops(X,operand,Y,alpha);
%%
%%  INPUT:        X = a matrix or a scalar
%%                    or a CELL ARRAY consisting only of matrices
%%          operand = sym, transpose, triu, tril,
%%                    real, imag, sqrt, abs, max, min, nnz,
%%                    spdiags, ones, zeros, norm, sum, row-norm, blk-norm
%%                    rank1, rank1inv, inv
%%                    +,  -, *, .*,  ./, .^
%%     Y (optional) = a matrix or a scalar
%%                    or a CELL ARRAY consisting only of matrices
%% alpha (optional) = a scalar
%%                    or the variable blk.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%******************************************************************

function Z = ops(X,operand,Y,alpha)

spdensity = 0.4;

if (nargin == 2)
    if strcmp(operand,'sym');
        if ~iscell(X);
            [m,n] = size(X);
            if (m == n);
                Z = (X+X')/2;
            elseif (n == 1);
                Z = X;
            else
                error('X must be square matrix or a column vector');
            end;
        else
            Z = cell(size(X));
            for p = 1:length(X);
                [m,n] = size(X{p});
                if (m == n);
                    Z{p} = (X{p}+X{p}')/2;
                elseif (n == 1);
                    Z{p} = X{p};
                else
                    error('X{p} must be square matrix or a column vector');
                end;
            end;
        end;
    elseif strcmp(operand,'sqrt') || strcmp(operand,'abs') || ...
            strcmp(operand,'real') || strcmp(operand,'imag');
        if ~iscell(X);
            eval(['Z = ',operand,'(X);']);
        else
            Z = cell(size(X));
            for p = 1:length(X);
                eval(['Z{p} = ',operand,'(X{p});']);
            end;
        end;
    elseif strcmp(operand,'max') || strcmp(operand,'min') || ...
            strcmp(operand,'sum');
        if ~iscell(X);
            eval(['Z = ',operand,'(X);']);
        else
            Z = [];
            for p = 1:length(X);
                eval(['Z = [Z  ',operand,'(X{p})','  ];']);
            end;
        end;
        eval(['Z = ',operand,'(Z);']);
    elseif strcmp(operand,'transpose') || strcmp(operand,'triu') || ...
            strcmp(operand,'tril');
        if ~iscell(X);
            if (size(X,1) == size(X,2));
                Z = feval(operand,X);
            elseif (size(X,2) == 1);
                Z = X;
            else
                error('X must be square matrix or a column vector');
            end;
        else
            Z = cell(size(X));
            for p = 1:length(X);
                if (size(X{p},1) == size(X{p},2));
                    Z{p} = feval(operand,X{p});
                elseif (size(X{p},2) == 1);
                    Z{p} = X{p};
                else
                    error('X{p} must be square matrix or a column vector');
                end;
            end;
        end;
    elseif strcmp(operand,'norm');
        if ~iscell(X);
            Z = full(sqrt(sum(sum(X.*X))));
        else
            Z = 0;
            for p = 1:length(X); Z = Z + sum(sum(X{p}.*X{p})); end;
            Z = sqrt(Z);
        end;
    elseif strcmp(operand,'blk-norm');
        if ~iscell(X);
            Z = full(sqrt(sum(sum(X.*X))));
        else
            Z = zeros(length(X),1);
            for p = 1:length(X); Z(p) = sum(sum(X{p}.*X{p})); end;
            Z = sqrt(Z);
        end;
    elseif strcmp(operand,'inv');
        if ~iscell(X);
            [m,n] = size(X); n2 = n*n;
            if (m==n)
                Z = inv(X);
                if (nnz(Z) > spdensity*n2)
                    Z = full(Z);
                else
                    Z = sparse(Z);
                end
            elseif (m > 1 && n == 1);
                Z = 1./X;
                if (nnz(Z) > spdensity*n)
                    Z = full(Z);
                else
                    Z = sparse(Z);
                end
            end
        else
            Z = cell(size(X));
            for p = 1:length(X);
                [m,n] = size(X{p}); n2 = n*n;
                if (m==n)
                    Z{p} = inv(X{p});
                    if (nnz(Z{p}) > spdensity*n2)
                        Z{p} = full(Z{p});
                    else
                        Z{p} = sparse(Z{p});
                    end
                elseif (m > 1 && n == 1);
                    Z{p} = 1./X{p};
                    if (nnz(Z{p}) > spdensity*n)
                        Z{p} = full(Z{p});
                    else
                        Z{p} = sparse(Z{p});
                    end
                end
            end
        end
    elseif strcmp(operand,'getM');
        if ~iscell(X);
            Z = size(X,1);
        else
            for p = 1:length(X); Z(p) = size(X{p},1); end; %#ok
            Z = sum(Z);
        end;
    elseif strcmp(operand,'nnz');
        if ~iscell(X);
            Z = nnz(X);
        else
            for p = 1:length(X);
                Z(p) = nnz(X{p}); %#ok
            end;
            Z = sum(Z);
        end;
    elseif strcmp(operand,'ones');
        if ~iscell(X);
            Z = ones(size(X));
        else
            Z = cell(size(X));
            for p = 1:length(X);
                Z{p} = ones(size(X{p}));
            end
        end
    elseif strcmp(operand,'zeros');
        if ~iscell(X);
            [m,n] = size(X);
            Z = sparse(m,n);
        else
            Z = cell(size(X));
            for p = 1:length(X);
                [m,n] = size(X{p});
                Z{p} = sparse(m,n);
            end
        end
    elseif strcmp(operand,'identity');
        blk = X;
        Z = cell(size(blk,1),1);
        for p = 1:size(blk,1)
            pblk = blk(p,:); n = sum(pblk{2});
            if strcmp(pblk{1},'s')
                Z{p} = speye(n,n);
            elseif strcmp(pblk{1},'q')
                s = 1+[0, cumsum(pblk{2})];
                len = length(pblk{2});
                Z{p} = zeros(n,1);
                Z{p}(s(1:len)) = ones(len,1);
            elseif strcmp(pblk{1},'l')
                Z{p} = ones(n,1);
            elseif strcmp(pblk{1},'u')
                Z{p} = zeros(n,1);
            end
        end
    elseif strcmp(operand,'row-norm');
        if ~iscell(X);
            if (size(X,2) == size(X,1));
                Z = sqrt(sum((X.*conj(X))'))'; %#ok
            elseif (size(X,2) == 1);
                Z = abs(X);
            end
        else
            Z = cell(size(X));
            for p = 1:length(X);
                if (size(X{p},2) == size(X{p},1));
                    Z{p} = sqrt(sum((X{p}.*conj(X{p}))'))'; %#ok
                elseif (size(X{p},2) == 1);
                    Z{p} = abs(X{p});
                end
            end
        end
    end
end
%%
if (nargin == 3)
    if strcmp(operand,'spdiags');
        if ~iscell(Y);
            [m,n] = size(Y);
            if (m == n);
                Z = spdiags(X,0,m,n);
            else
                Z = X;
            end
        else
            Z = cell(size(Y));
            for p = 1:length(Y);
                [m,n] = size(Y{p});
                if (m == n);
                    Z{p} = spdiags(X{p},0,m,n);
                else
                    Z{p} = X{p};
                end
            end
        end
    elseif strcmp(operand,'inprod')
        if ~iscell(X) && ~iscell(Y)
            Z = (Y'*X)';
        elseif iscell(X) && iscell(Y)
            Z = zeros(size(X{1},2),1);
            for p=1:length(X)
                Z = Z + (Y{p}'*X{p})';
            end
        end
    elseif strcmp(operand,'+') || strcmp(operand,'-') || ...
            strcmp(operand,'/') || strcmp(operand,'./') || ...
            strcmp(operand,'*') || strcmp(operand,'.*') || ...
            strcmp(operand,'.^');
        if (~iscell(X) && ~iscell(Y));
            eval(['Z = X',operand,'Y;']);
        elseif (iscell(X) && iscell(Y))
            Z = cell(size(X));
            for p = 1:length(X);
                if (size(X{p},2) == 1) && (size(Y{p},2) == 1) && ...
                        (strcmp(operand,'*') || strcmp(operand,'/'));
                    eval(['Z{p} = X{p}.',operand,'Y{p};']);
                else
                    eval(['Z{p} = X{p} ',operand,'Y{p};']);
                end
            end
        elseif (iscell(X) && ~iscell(Y));
            if (length(Y) == 1); Y = Y*ones(length(X),1); end
            Z = cell(size(X));
            for p = 1:length(X);
                eval(['Z{p} = X{p}',operand,'Y(p);']);
            end
        elseif (~iscell(X) && iscell(Y));
            Z = cell(size(Y));
            if (length(X) == 1); X = X*ones(length(Y),1); end
            for p = 1:length(Y);
                eval(['Z{p} = X(p)',operand,'Y{p};']);
            end
        end
    else
        error([operand,' is not available, check input arguments']);
    end
end
%%
if (nargin == 4)
    if strcmp(operand,'rank1') || strcmp(operand,'rank1inv');
        Z = cell(size(alpha,1),1);
        for p = 1:size(alpha,1);
            if ~strcmp(alpha{p,1},'diag');
                blktmp = alpha{p,2};
                if (length(blktmp) == 1);
                    if strcmp(operand,'rank1');
                        Z{p} = (X{p}*Y{p}' + Y{p}*X{p}')/2;
                    else
                        Z{p} = 2./(X{p}*Y{p}' + Y{p}*X{p}');
                    end
                else
                    Xp = X{p};
                    Yp = Y{p};
                    n = sum(blktmp);
                    Zp = sparse(n,n);
                    s = [0 cumsum(blktmp)];
                    if strcmp(operand,'rank1');
                        for i = 1:length(blktmp)
                            pos = s(i)+1 : s(i+1);
                            x = Xp(pos);
                            y = Yp(pos);
                            Zp(pos,pos) = sparse((x*y' + y*x')/2); %#ok
                        end;
                        Z{p} = Zp;
                    else
                        for i = 1:length(blktmp)
                            pos = s(i)+1 : s(i+1);
                            x = Xp(pos);
                            y = Yp(pos);
                            Zp(pos,pos) = sparse(2./(x*y' + y*x')); %#ok
                        end
                        Z{p} = Zp;
                    end
                end
            elseif strcmp(alpha{p,1},'diag');
                if strcmp(operand,'rank1');
                    Z{p} = X{p}.*Y{p};
                else
                    Z{p} = 1./(X{p}.*Y{p});
                end
            end
        end
    elseif strcmp(operand,'+') || strcmp(operand,'-');
        if ~iscell(X) && ~iscell(Y);
            eval(['Z = X',operand,'alpha*Y;']);
        elseif (iscell(X) && iscell(Y));
            Z = cell(size(X));
            if (length(alpha) == 1);
                alpha = alpha*ones(length(X),1);
            end
            if strcmp(operand,'+'),
                for p = 1:length(X)
                    Z{p} = X{p} + alpha(p) * Y{p};
                end
            else
                for p = 1:length(X);
                    Z{p} = X{p} - alpha(p) * Y{p};
                end
            end
        else
            error('X, Y are different objects');
        end
    else
        error([operand,' is not available']);
    end
end
%%============================================================

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% VISUALIZEMCMCMARGINALS
%
% This function accepts a list of sample lists, each from a different MCMC run.  It then visualizes
% the estimated marginals for each variable in V over the lifetime of the MCMC run.
%
% samples_list - a list of sample lists; each sample list is a m-by-n matrix where m is the
% number of samples and n is the number of variables in the state of the Markov chain
%
% V - an array of variables
% D - the dimensions of the variables in V
% F - a list of factors (used in computing likelihoods of each sample)
% window_size - size of the window over which to aggregate samples to compute the estimated
%               marginal at a given time
% ExactMarginals - the exact marginals of V (optional)
%
% Copyright (C) Daphne Koller, Stanford University, 2012

function VisualizeMCMCMarginals(samples_list, V, D, F, window_size, ExactMarginals, tname)

for i = 1:length(V)
    figure;
    v = V(i);
    d = D(i);
    title(['Marginal for Variable ', num2str(v)]);
    if exist('ExactMarginals') == 1, M = ExactMarginals(i); end;
    for j = 1:length(samples_list)
        samples_v = samples_list{j}(:, v);
        indicators_over_time = zeros(length(samples_v), d);
        for k = 1:length(samples_v)
            indicators_over_time(k, samples_v(k)) = 1;
        end

        % estimated_marginal = cumsum(indicators_over_time, 1);
        estimated_marginal = [];
        for k = 1:size(indicators_over_time, 2)
            estimated_marginal = [estimated_marginal, smooth(indicators_over_time(:, k), window_size)];
        end
        % Prune ends
        estimated_marginal = estimated_marginal(window_size/2:end - window_size/2, :);

        estimated_marginal = estimated_marginal ./ ...
            repmat(sum(estimated_marginal, 2), 1, size(estimated_marginal, 2));
        hold on;
        plot(estimated_marginal, '-', 'LineWidth', 2);
        title(['Est marginals for entry ' num2str(i) ' of samples for ' tname])
        if exist('M') == 1
            plot([1; size(estimated_marginal, 1)], [M.val; M.val], '--', 'LineWidth', 3);
        end
        set(gcf,'DefaultAxesColorOrder', rand(d, 3));
    end
end

% Visualize likelihood of sample at each time step
all_likelihoods = [];
for i = 1:length(samples_list)
    samples = samples_list{i};
    likelihoods = [];
    for j = 1:size(samples, 1)
        likelihoods = [likelihoods; LogProbOfJointAssignment(F, samples(j, :))];
    end
    all_likelihoods = [all_likelihoods, likelihoods];
end
figure;
title('Likelihoods')
plot(all_likelihoods, '-', 'LineWidth', 2);

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function area = disk01_area ( )

%*****************************************************************************80
%
%% DISK01_AREA returns the area of the unit disk.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real AREA, the area of the unit disk.
%
  r = 1.0;
  area = pi * r * r;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%IBBOX Find bounding box
%
% BOX = IBBOX(P) is the minimal bounding box that contains the points
% described by the columns of P (2xN).
%
% BOX = IBBOX(IM) as above but the box minimally contains the non-zero
% pixels in the image IM.
%
% Notes::
% - The bounding box is a 2x2 matrix [XMIN XMAX; YMIN YMAX].

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.

function box = ibbox(I)

    if numrows(I) == 2
        % input is a set of points
        u = I(1,:);
        v = I(2,:);
    else
        % input is an image, find the non-zero elements
        [v,u] = find(I);
    end
    umin = min(u);
    umax = max(u);
    vmin = min(v);
    vmax = max(v);

    box = [umin umax; vmin vmax];

author: petercorke

repo: machinevision-toolbox-matlab


text:
function vis_debug_bbox_pred(imdb, rcnn_model, bbox_reg, cls)
% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

image_ids = imdb.image_ids;
min_overlap = 0.3;
feat_opts = bbox_reg.training_opts;

roidb = imdb.roidb_func(imdb);
class_id = imdb.class_to_id(cls);

bbox_regressor = bbox_reg.models{class_id};

perm1 = randperm(length(image_ids));
%for i = 1:length(image_ids)
for i = perm1
  fprintf('%s: debug (%s) %d/%d\n', procid(), imdb.name, i, length(image_ids));

  if exist('cls','var')
    if isempty(find(roidb.rois(i).class == class_id))
      continue;
    end
  end

  d = rcnn_load_cached_pool5_features(feat_opts.cache_name, ...
      imdb.name, image_ids{i});
  if isempty(d.feat)
    continue;
  end

  d.feat = rcnn_pool5_to_fcX(d.feat, feat_opts.layer, rcnn_model);
  d.feat = rcnn_scale_features(d.feat, feat_opts.feat_norm_mean);

  if feat_opts.binarize
    d.feat = (d.feat > 0);
  end

  im = imread(imdb.image_at(i));

  if exist('cls','var')
    sel_gt = find(d.class == class_id);
  else
    sel_gt = find(d.class > 0);
  end
  gt_boxes = d.boxes(sel_gt, :);
  gt_classes = d.class(sel_gt);

  if exist('cls','var')
    max_ov = d.overlap(:, class_id);
  else
    max_ov = max(d.overlap, [], 2);
  end
  sel_ex = find(max_ov >= min_overlap);
  ex_boxes = d.boxes(sel_ex, :);
  ex_feat = d.feat(sel_ex, :);

  %for j = 1:size(ex_boxes, 1)
  perm = randperm(size(ex_boxes,1), min(10, size(ex_boxes,1)));
  for j = perm
    ex_box = ex_boxes(j, :);
    ov = boxoverlap(gt_boxes, ex_box);
    [max_ov, assignment] = max(ov);
    gt_box = gt_boxes(assignment, :);
    cls = gt_classes(assignment);

    src_w = ex_box(3) - ex_box(1) + eps;
    src_h = ex_box(4) - ex_box(2) + eps;
    src_ctr_x = ex_box(1) + 0.5*src_w;
    src_ctr_y = ex_box(2) + 0.5*src_h;

    gt_w = gt_box(3) - gt_box(1) + eps;
    gt_h = gt_box(4) - gt_box(2) + eps;
    gt_ctr_x = gt_box(1) + 0.5*gt_w;
    gt_ctr_y = gt_box(2) + 0.5*gt_h;

%    dst_ctr_x = (gt_ctr_x - src_ctr_x) * 1/src_w;
%    dst_ctr_y = (gt_ctr_y - src_ctr_y) * 1/src_h;
%    dst_scl_w = log(gt_w / src_w);
%    dst_scl_h = log(gt_h / src_h);
%
%    target = [dst_ctr_x dst_ctr_y dst_scl_w dst_scl_h];

    pred_box = rcnn_predict_bbox_regressor(bbox_regressor, ex_feat(j,:), ex_box);

    pred_box(1) = max(pred_box(1), 1);
    pred_box(2) = max(pred_box(2), 1);
    pred_box(3) = min(pred_box(3), size(im,2));
    pred_box(4) = min(pred_box(4), size(im,1));

    ovs = boxoverlap(cat(1, pred_box, ex_box), gt_box);

    showboxesc(im, gt_box, 'g', '-');
    showboxesc([], ex_box, 'b', '-');
    showboxesc([], pred_box, 'r', '-');
    hold on;
    plot(gt_ctr_x, gt_ctr_y, 'gd');
    plot(src_ctr_x, src_ctr_y, 'rd');
    hold off;
    title(sprintf('green = GT box;  blue = original box;  red = predicted box;  orig %.3f  pred %.3f', ovs(2), ovs(1)));

    pause;
  end
end

author: rbgirshick

repo: rcnn

save_path: github-repos/MATLAB/rbgirshick-rcnn


text:
function [net, stats] = loadState_deploy(fileName)
% -------------------------------------------------------------------------
load(fileName, 'net', 'stats') ;
net = dagnn.DagNN.loadobj(net) ;

repo: DBCNN


text:
% [INPUT]
% data = A float t-by-n matrix containing the time series to be converted into rolling windows.
% bw = An integer [2,Inf) representing the dimension of each rolling window.
% truncate = A boolean that indicates whether to exclude all the rolling windows with a dimension less than the bandwidth (optional, default=false).
%
% [OUTPUT]
% windows = A column cell array of bw-by-n float matrices representing the rolling windows.
%
% [NOTES]
% If the number of observations is less than or equal to the specified bandwidth, a single rolling window containing all the observations is returned.

function windows = extract_rolling_windows(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('data',@(x)validateattributes(x,{'double'},{'2d' 'nonempty'}));
        ip.addRequired('bw',@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 2 'scalar'}));
        ip.addOptional('truncate',false,@(x)validateattributes(x,{'logical'},{'scalar'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    data = ipr.data;
    bw = ipr.bw;
    truncate = ipr.truncate;

    nargoutchk(1,1);

    windows = extract_rolling_windows_internal(data,bw,truncate);

end

function windows = extract_rolling_windows_internal(data,bw,truncate)

    t = size(data,1);

    if (bw >= t)
        windows = cell(1,1);
        windows{1} = data;
        return;
    end

    limit = t - bw + 1;

    if (truncate)
        windows = cell(limit,1);

        for i = 1:limit
            windows{i} = data(i:bw+i-1,:);
        end
    else
        windows = cell(t,1);

        k = max(round(nthroot(bw,1.81),0),5);

        for i = 1:(bw - 1)
            windows{i} = data(1:max(i,k),:);
        end

        for i = 1:limit
            windows{i+bw-1} = data(i:bw+i-1,:);
        end
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function [TKenergy]=TKEO(x)
%
%% Utility function to calculate measures based on the nonlinear energy operator
%
% Function to estimate the TKEO of input data, follows the classical rule
% of Teager and Kaiser
%
% Inputs:  x            -> any time-series signal (vector)
%
% =========================================================================
% Output:  TKenergy     -> (nonlinear) Teager-Kaiser Energy Operator (TKEO)
% =========================================================================
%
% Part of the "Speech Disorders" Toolbox
%
% -----------------------------------------------------------------------
% Useful references:
%
% 1) J. Kaiser: On a simple algorithm to calculate the 'energy' of a
%    signal, Proc. IEEE International Conference on Acoustics, Speech, and
%    Signal Processing (ICASSP '90), pp. 381-384, Albuquerque, NM, USA,
%    April 1990
%
% -----------------------------------------------------------------------
%
% Last modified on 24 August 2014
%
% Copyright (c) Athanasios Tsanas, 2014
%
% ********************************************************************
% If you use this program please cite:
%
% 1) A. Tsanas: "Accurate telemonitoring of Parkinson's disease symptom
%    severity using nonlinear speech signal processing and statistical
%    machine learning", D.Phil. thesis, University of Oxford, 2012
% ********************************************************************
%
% For any question, to report bugs, or just to say this was useful, email
% tsanasthanasis@gmail.com

%% Algorithm computation

% The algorithm is computed using the instantaneous value squared minus the
% previous step value times the next step value:[x_n]^2-[x_(n-1)]*(x_(n+1)]
%
%
% 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.

% 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.
%

data_length=length(x);
TKenergy=zeros(data_length,1);

TKenergy(1)=(x(1))^2; % first sample

for n=2:data_length-1
    TKenergy(n)=(x(n))^2-x(n-1)*x(n+1);
end

%  TKenergy(2:data_length-1) = x(2:data_length-1).^2 - x(1:data_length-2).*x(3:data_length); % alternative, vectorized version

TKenergy(data_length)=(x(data_length))^2; % last sample

author: fernandoandreotti