text:
%ADABOOSTC
%
% [W,V,ALF] =  ADABOOSTC(A,CLASSF,N,RULE,VERBOSE);
%
% INPUT
%   A       Dataset
%   CLASSF  Untrained weak classifier
%   N       Number of classifiers to be trained
%   RULE    Combining rule (default: weighted voting)
%   VERBOSE Suppress progress report if 0 (default)
%
% OUTPUT
%   W       Combined trained classifier
%   V       Cell array of all classifiers
%           Use VC = stacked(V) for combining
%   ALF     Weights
%
% DESCRIPTION
%
% Computation of a combined classifier according to adaboost.
%
% In total N weighted versions of the training set A are generated
% iteratevely and used for the training of the specified classifier.
% Weights, to be used for the probabilities of the objects in the training
% set to be selected, are updated according to the Adaboost rule.
%
% The entire set of generated classifiers is given in V.
% The set of classifier weigths, according to Adaboost is returned in ALF
%
% Various aggregating possibilities can be given in
% the final parameter rule:
% []:      WVOTEC, weighted voting.
% VOTEC    voting
% MEANC    sum rule
% AVERAGEC averaging of coeffients (for linear combiners)
% PRODC    product rule
% MAXC     maximum rule
% MINC     minimum rule
% MEDIANC  median rule
%
% REFERENCE
% Ji Zhu, Saharon Rosset, Hui Zhou and Trevor Hastie,
% Multiclass Adaboost. A multiclass generalization of the Adaboost
% algorithm, based on a generalization of the exponential loss.
% http://www-stat.stanford.edu/~hastie/Papers/samme.pdf
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands
% (Multiclass correction by Marcin Budka, Bournemouth Univ., UK)

%function [W,V,alf] = adaboostc(a,clasf,n,rule,verbose)
function [out,V,alf] = adaboostc(varargin)

%% INITIALISATION
argin = setdefaults(varargin,[],nmc,100,[],0);
if mapping_task(argin,'definition')

  out = define_mapping(argin,'untrained','Adaboost');

%% TRAINING
elseif mapping_task(argin,'training')

  [a,clasf,n,rule,verbose] = deal(argin{:});
  [m,k,c] = getsize(a);
  V = [];
  laba = getlab(a);
  u = ones(m,1)/m;			% initialise object weights
  alf = zeros(1,n);			% space for classifier weights
  isseparable = 0;          % check if we can make 0 error
  if verbose && k == 2
    figure(verbose);
    scatterd(a);
  end

  %% generate n classifiers
  for i = 1:n
    b = gendatw(a,u,m);             % sample training set
    b = setprior(b,getprior(a));	% use original priors
    w = b*clasf;                    % train weak classifier
    ra = a*w;                       % test weak classifier

    if verbose && k == 2
      plotc(w,1); drawnow
    end

    labc = labeld(ra);
    diff = sum(labc~=laba,2)~=0;	% objects erroneously classified
    erra = sum((diff).*u);          % weighted error on original dataset

    if (erra==0)
        isseparable = 1;
        V = w;
        break;
    end
    if (erra < (1-1/c))        % if classifier better then random guessing...
      alf(i) = 0.5*(log((1-erra)/erra) + log(c-1));
      correct = find(diff==0); % find correctly classified objects
      wrong = find(diff==1);   % find incorrectly classified objects
      u(correct) = u(correct)*exp(-alf(i));	% give them the ...
      u(wrong) = u(wrong)*exp(alf(i));	  	% proper weights
      u = u./sum(u);                        % normalize weights
    else
      alf(i) = 0;
    end

    if verbose
      disp([erra alf(i) sum(alf)])
    end
    V = [V w];                       % store all classifiers

  end

  %% combine and return
  if isseparable
      W = V;
      W = setname(W,['Boosted ',getname(V)]);
  else
    if isempty(rule)
        W = wvotec(V,alf);             % default is weighted combiner
    else
        W = traincc(a,V,rule);         % otherwise, use user supplied combiner
    end
  end

  if verbose > 0 && k == 2
    plotc(W,'r',3)
    ee = a*W*testc;
    title(['Error: ', num2str(ee)]);
  end

  out = W;

else
  error('Illegal call')
end

return

author: marianux

repo: ecg-kit

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


text:
function accuracies = hkl_test(model,outputs,Xtest,varargin);
%%%%%%%%%%%%%%%%%%%%%%%
% HKL
%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% required parameters
%
% Xtest					test data (input)
%
% optional parameters
%
% Ytest					test data (output), to compute predictive performance

Ytest = [];                         % test data (response)

% READ OPTIONAL PARAMETERS
args = varargin;
nargs = length(args);
for i=1:2:nargs
	switch args{i},
		case 'Ytest',           Ytest = args{i+1};
	end
end

data = model.data;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NORMALIZE DATA (in input space)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~strcmp(model.kernel,'base kernels') || ~strcmp(model.kernel,'base kernels-mkl') || ~strcmp(model.kernel,'base kernels-bimkl')
    switch model.data_normalization
        case  'scale'
            [ntest , p ] = size(Xtest);
            Xtest = Xtest - repmat(model.mean,ntest,1);
            Xtest = Xtest ./ repmat(model.std,ntest,1);
        case  'center'
            [ntest , p ] = size(Xtest);
            Xtest = Xtest - repmat(model.mean,ntest,1);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PREPROCESS DATA
% the structure data store all data necessary to compute the kernel for a given node
%
% for directed grids where each kernel in the graph is the product of base kernels
% each base kernels is represented through a Cholesky decomposition
%
% data.n				number of observations
% data.p				dimension of the grid
% data.q				maximal depth of the grid
% data.qs				number of kernels in each direction of the grid
% data.Xs				storing all dimensions in the same vector
% data.ind_Xs			indices corresponding to each base kernel
% data.d_Xs				rank of each base kernel
%
% for directed grids, the (unique root) is always the constant kernel (and thus zeroed out by centering)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch model.kernel

	case {'polynomial','polynomial-mkl','polynomial-bimkl'}

		datatest = get_kernel_data(Xtest,model.kernel,model.kernel_params);
		datatest.dag_type =  'grid + input_space';

	case {'hermite', 'hermite-mkl', 'hermite-bimkl'}

		datatest = get_kernel_data(Xtest,model.kernel,model.kernel_params);
		datatest.dag_type =  'grid + input_space';

	case {'anova','anova-mkl','anova-bimkl'}

		datatest = get_kernel_data(Xtest,model.kernel,model.kernel_params,model.X,model.data);
		datatest.dag_type =  'grid + input_space';

	case {'spline','spline-mkl','spline-bimkl'}

		datatest = get_kernel_data(Xtest,model.kernel,model.kernel_params,model.X,model.data);
		datatest.dag_type =  'grid + input_space';

    case {'base kernels','base kernels-mkl','base kernels-bimkl'}

		datatest = get_kernel_data(Xtest,model.kernel,model.kernel_params,model.X,model.data);
		datatest.dag_type =  'grid + input_space';

	case {'gauss-hermite','gauss-hermite-mkl','gauss-hermite-bimkl'}

		datatest = get_kernel_data(Xtest,model.kernel,model.kernel_params);
		datatest.dag_type =  'grid + input_space';

	case {'gauss-hermite-full','gauss-hermite-full-mkl','gauss-hermite-full-bimkl'}

		datatest = get_kernel_data(Xtest,model.kernel,model.kernel_params,model.X,model.data);
		datatest.dag_type =  'grid + input_space';

end

switch model.kernel

	case {'hermite-mkl', 'gauss-hermite-mkl', 'gauss-hermite-full-mkl', 'polynomial-mkl', 'anova-mkl',  'spline-mkl' }
		dag_type = 'mkl + input_space';
		data.dag_type = dag_type;
		datatest.dag_type = dag_type;
	case {'hermite-bimkl', 'gauss-hermite-bimkl', 'gauss-hermite-full-bimkl', 'polynomial-bimkl', 'anova-bimkl',  'spline-bimkl' }
		dag_type = 'bimkl + input_space';
		data.dag_type = dag_type;
		datatest.dag_type = dag_type;
end

% NORMALIZATION OF BASE KERNELS
switch model.kernel_normalization
	case 'center'  % center before taking product of base kernels
		switch datatest.dag_type
			case { 'grid + input_space', 'mkl + input_space', 'bimkl + input_space'}
				for i=1:p
					for j=2:data.q+1
						datatest.Xs(:,data.ind_Xs{i}{j}) = datatest.Xs(:,data.ind_Xs{i}{j}) - repmat( model.meanXs{i}{j},n,1);
					end
				end

		end

	case 'scale'  % scale (and center) before taking product of base kernels
		switch datatest.dag_type
			case { 'grid + input_space', 'mkl + input_space', 'bimkl + input_space'}
				for i=1:p
					for j=2:data.q+1
						datatest.Xs(:,data.ind_Xs{i}{j}) = datatest.Xs(:,data.ind_Xs{i}{j}) - repmat( model.meanXs{i}{j},ntest,1);
						datatest.Xs(:,data.ind_Xs{i}{j}) = datatest.Xs(:,data.ind_Xs{i}{j}) / model.stdXs{i}{j};
					end
				end

		end

end

accuracies.predtest = NaN*ones(length(Ytest),length(outputs));
accuracies.predtest_unreg = NaN*ones(length(Ytest),length(outputs));

if ~isempty(Ytest)
	accuracies.testing_error = NaN*ones(1,length(outputs));
	accuracies.testing_error_unreg = NaN*ones(1,length(outputs));
	switch model.loss
		case 'logistic'
			accuracies.testing_error_class = NaN*ones(1,length(outputs));
			accuracies.testing_error_class_unreg = NaN*ones(1,length(outputs));
	end
end

for ilambda = 1:length(outputs)

	[predtrain,predtest] = predict_train_test(datatest,model.data,outputs{ilambda}.hull,outputs{ilambda}.alpha,outputs{ilambda}.b,outputs{ilambda}.zeta);
	[predtrain_unreg,predtest_unreg] = predict_train_test_unreg(datatest,model.data,model.Y,Ytest,outputs{ilambda}.hull);
	accuracies.predtest(:,ilambda) = predtest;
	accuracies.predtrain(:,ilambda) = predtrain;
	accuracies.predtest_unreg(:,ilambda) = predtest_unreg;
	accuracies.predtrain_unreg(:,ilambda) = predtrain_unreg;
	if ~isempty(Ytest)
		switch model.loss
			case 'square'
				accuracies.training_error(ilambda) =  sum( ( model.Y - predtrain ).^2 ) /length(model.Y);
				accuracies.testing_error(ilambda) =  sum( ( Ytest - predtest ).^2 ) /length(Ytest);
				accuracies.training_error_unreg(ilambda) =  sum( ( model.Y - predtrain_unreg ).^2 ) /length(model.Y);
				accuracies.testing_error_unreg(ilambda) =  sum( ( Ytest - predtest_unreg ).^2 ) /length(Ytest);

			case 'logistic'
				accuracies.training_error_class(ilambda) =  sum( abs( model.Y - (sign(predtrain-.5)+1)/2 ) ) /length(model.Y);
				accuracies.testing_error_class(ilambda) =  sum( abs( Ytest - (sign(predtest-.5)+1)/2 ) ) /length(Ytest);
				accuracies.training_error_class_unreg(ilambda) =  sum( abs( model.Y - (sign(predtrain_unreg-.5)+1)/2 ) ) /length(model.Y);
				accuracies.testing_error_class_unreg(ilambda) =  sum( abs( Ytest - (sign(predtest_unreg-.5)+1)/2 ) ) /length(Ytest);
				accuracies.training_error(ilambda) =  sum(  model.Y .* log( 1 + exp( -predtrain) ) + ( 1 - model.Y ) .* log( 1 + exp(predtrain) )  ) /length(model.Y);
				accuracies.testing_error(ilambda) =  sum(  Ytest .* log( 1 + exp( -predtest) ) + ( 1 - Ytest ) .* log( 1 + exp(predtest) )  ) /length(Ytest);
				accuracies.training_error_unreg(ilambda) =  sum(  model.Y .* log( 1 + exp( -predtrain_unreg) ) + ( 1 - model.Y ) .* log( 1 + exp(predtrain_unreg) )  ) /length(model.Y);
				accuracies.testing_error_unreg(ilambda) =  sum(  Ytest .* log( 1 + exp( -predtest_unreg) ) + ( 1 - Ytest ) .* log( 1 + exp(predtest_unreg) )  ) /length(Ytest);
		end

	end

end

repo: REST


text:
classdef SettingsMeshPlotter < AbstractSettings

    properties (Access = protected)
        defaultParamsName = 'paramsMeshPlotter.json'
    end

    properties (Access = public)
        mesh
        isBackground
        faceColor
        faceAlpha
        edgeAlpha
    end

    methods (Access = public)

        function obj = SettingsMeshPlotter(varargin)
            switch nargin
                case 1
                    obj.loadParams(varargin{1});
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function im = crop_qim(im, bbx, imgsz)
% CROP_QIM crops query image with defined bounding box.
%
%   IM = crop_qim(IM, BBX)  Crop IM with BBX.
%
%   IM = crop_qim(IM, BBX, IMG_SIZE) Resize image so that longer edge is IMG_SIZE.
%     After that crop IM with BBX.

	allow_increase = 0;

	if size(im, 1) == 1, im = imread(im); end
	if ~exist('imgsz'),
		im = im(bbx(2):min(bbx(4),size(im,1)), bbx(1):min(bbx(3),size(im,2)), :);
		return;
	end

	bbx = uint32(max(imgsz * (bbx + 1) / max(size(im,1), size(im,2)), 1));
	im = imresizemaxd(im, imgsz, allow_increase);
	im = im(bbx(2):min(bbx(4),size(im,1)), bbx(1):min(bbx(3),size(im,2)), :);

author: filipradenovic

repo: cnnimageretrieval


text:
clear; clc;
disp('compile native libraries')
mex -O -largeArrayDims sparse_inp.c
mex -O -largeArrayDims sparse_update.c
disp('done compiling.')
test_consist1();
test_consist2();

author: andrewssobral

repo: lrslibrary


text:
% @KTENSOR
%
% Files
%   arrange      - Arranges the rank-1 components of a ktensor.
%   datadisp     - Special display of a ktensor.
%   disp         - Command window display for a ktensor.
%   display      - Command window display for a ktensor.
%   double       - Convert a ktensor to a double array.
%   end          - Last index of indexing expression for ktensor.
%   extract      - Creates a new ktensor with only the specified components.
%   fixsigns     - Fix sign ambiguity of a ktensor.
%   full         - Convert a ktensor to a (dense) tensor.
%   innerprod    - Efficient inner product with a ktensor.
%   isequal      - True if each component of two ktensor's is numerically equal.
%   ktensor      - Tensor stored as a Kruskal operator (decomposed).
%   minus        - Binary subtraction for ktensor.
%   mtimes       - Implement A*B (scalar multiply) for ktensor.
%   mttkrp       - Matricized tensor times Khatri-Rao product for ktensor.
%   ncomponents  - Number of components for a ktensor.
%   ndims        - Number of dimensions for a ktensor.
%   norm         - Frobenius norm of a ktensor.
%   normalize    - Normalizes the columns of the factor matrices.
%   nvecs        - Compute the leading mode-n vectors for a ktensor.
%   permute      - Permute dimensions of a ktensor.
%   plus         - Binary addition for ktensor.
%   redistribute - Distribute lambda values to a specified mode.
%   score        - Checks if two ktensors match except for permutation.
%   size         - Size of ktensor.
%   subsasgn     - Subscripted assignement for ktensor.
%   subsref      - Subscripted reference for a ktensor.
%   times        - Element-wise multiplication for ktensor.
%   tocell       - Convert X to a cell array.
%   ttm          - Tensor times matrix for ktensor.
%   ttv          - Tensor times vector for ktensor.
%   uminus       - Unary minus for ktensor.
%   uplus        - Unary plus for a ktensor.

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [dset] = cell2dataset(inputData)
    dset = dataset();
    if (~isfield(inputData, 'data'))
        error('[cell2dataset] input dataset structure has no ".data" field');
    end
    if (~isfield(inputData, 'columnnames'))
        error('[cell2dataset] input dataset structure has no ".columnnames" field');
    end

    for i = 1:size(inputData.columnnames,2)
        if(isnumeric(inputData.data{1,i}))
            if(isempty(dset))
                dset = dataset({cell2mat(inputData.data(:,i)),cell2mat(inputData.columnnames(i))});
            else
                dset = [dset dataset({cell2mat(inputData.data(:,i)),cell2mat(inputData.columnnames(i))})];
            end
            disp('numeric')
        else
            if(isempty(dset))
                dset = dataset({inputData.data(:,i),cell2mat(inputData.columnnames(i))});
            else
                dset = [dset dataset({inputData.data(:,i),cell2mat(inputData.columnnames(i))})];
            end
            disp('other')
        end
    end
end

author: Sable

repo: mcbench-benchmarks


text:
%JSINGU Show the linearly dependent joints in a Jacobian matrix
%
% JSINGU(J) displays the linear dependency of joints in a Jacobian matrix.
% This dependency indicates joint axes that are aligned and causes singularity.
%
% See also SerialLink.jacobn.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB 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 RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function jsingu(J)

    % convert to row-echelon form
    [R, jb] = rref(J);

    depcols = setdiff( 1:numcols(J), jb);

    fprintf('%d linearly dependent joints:\n', length(depcols));
    for d=depcols
        fprintf('  q%d depends on: ', d)
        for k=find(R(:,d))
            fprintf('q%d ', k);
        end
        fprintf('\n');
    end

author: Allopart

repo: rbpf-gmapping


text:
function [ o, x, w ] = cn_jac_00_1 ( n, alpha, beta )

%*****************************************************************************80
%
%% CN_JAC_00_1 implements the midpoint rule for region CN_JAC.
%
%  Discussion:
%
%    The rule has order O = 1.
%
%    The rule has precision P = 0.
%
%    CN is the cube [-1,+1]^N with the Jacobi (beta) weight function
%
%      w(alpha,beta;x) = product ( 1 <= i <= n ) (1-x(i))^beta (1+x(i))^alpha.
%
%    with -1 < alpha, -1 < beta.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Input, real ALPHA, BETA, the parameters.
%    -1.0 < ALPHA, -1.0 < BETA.
%
%    Output, integer O, the order.
%
%    Output, real X(N,O), the abscissas.
%
%    Output, real W(O), the weights.
%
  if ( alpha <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CN_JAC_00_1 - Fatal error!\n' );
    fprintf ( 1, '  ALPHA <= -1.0\n' );
    error ( 'CN_JAC_00_1 - Fatal error!\n' );
  end

  if ( beta <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CN_JAC_00_1 - Fatal error!\n' );
    fprintf ( 1, '  BETA <= -1.0\n' );
    error ( 'CN_JAC_00_1 - Fatal error!\n' );
  end

  o = 1;

  expon = 0;
  volume = c1_jac_monomial_integral ( alpha, beta, expon );
  volume = volume ^ n;

  x = zeros ( n, o );
  w = zeros ( o, 1 );

  k = 0;
%
%  1 point.
%
  k = k + 1;
  w(k) = volume;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [J, grad] = linearRegCostFunction(X, y, theta, lambda)
%LINEARREGCOSTFUNCTION Compute cost and gradient for regularized linear
%regression with multiple variables
%   [J, grad] = LINEARREGCOSTFUNCTION(X, y, theta, lambda) computes the
%   cost of using theta as the parameter for linear regression to fit the
%   data points in X and y. Returns the cost in J and the gradient in grad

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost and gradient of regularized linear
%               regression for a particular choice of theta.
%
%               You should set J to the cost and grad to the gradient.
%

h_theta = X*theta;
J = (1/(2*m))*sum((h_theta - y).^2) + (lambda/(2*m))*(theta(2:end)'*theta(2:end));
grad = (1/m)*((h_theta - y)'*X)'
grad = grad + (lambda/m)*[0;theta(2:end,:)]

% =========================================================================

grad = grad(:);

end

author: anirudhjayaraman

repo: Machine-Learning


text:
function [err,v, typ] = BrikInfo_SectionValue (BRIKinfo, sSection)
%
%   [err,s, typ] = BrikInfo_SectionValue (sHead, sSection)
%
%Purpose:
%   Get the values of a section in the .HEAD string
%   This function is dedicated for BrikInfo function
%
%Input Parameters:
%   sHead, a string vector containing the .HEAD info
%   sSection, the section name such as: 'DATASET_DIMENSIONS'
%
%
%Output Parameters:
%   err : 0 No Problem
%       : 1 Mucho Problems
%   v a vector (or a string) containing the values in that section
%     v is empty (and err = 1) if the section can't be found
%   typ : is a string to indicate the type of parameter (string or number)
%      that v is
%
%Key Terms:
%
%More Info :
%  see also BrikInfo
%
%
%
%     Author : Ziad Saad
%     Date : Mon Oct 18 13:46:28 CDT 1999

%Define the function name for easy referencing
FuncName = 'BrikInfo_SectionValue';

%Debug Flag
DBG = 1;

%initailize return variables
err = 1;

N_BRIKinfo = length(BRIKinfo);

%itmp = findstr (BRIKinfo, sSection); This fails if sSection appears in places
%like History note and history note appears first.
%Nick OOsterhof suggested the following fix:
htmp = regexp (BRIKinfo, ['name\s*=\s*' sSection]);
if isempty(htmp) | length(htmp) > 1, v = []; return; end
itmp = findstr(BRIKinfo(htmp:N_BRIKinfo),sSection)+htmp-1;
if isempty(itmp), v = []; return; end

%findout where this section ends
inxt = strfind(BRIKinfo(itmp:N_BRIKinfo),'type');
inxt = setdiff(inxt, strfind(BRIKinfo(itmp:N_BRIKinfo),'_type')+1); %JM added to avoid getting stuck on data_type
if isempty(inxt), inxt = N_BRIKinfo; else inxt = inxt(1) + itmp -2; end

%Get the count value
[err,sn, jnk, strt,stp] = GetNextLine(BRIKinfo(itmp:inxt),2);
ic = findstr(sn,'count'); %watch it, spacing sensitive section
if ~isempty(ic)
  %skip the = sign
  sn = sn(ic+5:length(sn));
  ic = findstr(sn, '=');
  sn = sn(ic+1:length(sn));
else
  err = ErrEval(FuncName,'Err_could not find count field');
  return;
end
%get the count value
n = str2num(sn);
if ~n %empty field
  v = [];
  typ = '';
  return;
end

%advance to the next line, after count
itmp = itmp + stp + 1;

%read all values
Svals = BRIKinfo(itmp:inxt); %that's where the values are
inl = find(char(Svals) == 10 | char(Svals) == 13); %find any new lines left in Svals (13 was added to work in DOS (CR and LF))

if (~isempty(inl))
  Svals(inl) = ' '; %replace them by space characters
end

%if (sum(isletter(Svals))), %isletter fails when you have numbers such as 1.4e-3
v = zdeblank(Svals);
if isempty(v)
  typ = 'string';
elseif strcmp(v(1),'''') && strcmp(v(length(v)),'~')
  %This parameter is a string, do not change to numbers
  %remove first and last chars
  v = v(2:length(v)-1);
  typ = 'string';
else
  v = str2num(Svals);
  typ = 'number';
end

%pause

err = 0;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function planDBGUI(action, varargin)
% Gui to access and browse data in planDB.
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

if(nargin == 0)
    action = 'init';
end

switch upper(action)
    case 'INIT'
		screenSize = get(0,'ScreenSize');
		y = 700; %Initial size of figure in pixels.
		x = 700;

		h = figure('position', [(screenSize(3)-x)/2 (screenSize(4)-y)/2 x y], 'units', 'normalized');
        set(h, 'resizefcn', 'myResize');

        ud = []; %init UD.
		%Make bottom Axis, table it
		ud.handles.botAxis = axes('position', [.25 .05 .7 .4]);
   		table(ud.handles.botAxis);
        table(ud.handles.botAxis, 'SETDBCALLBACK', 'planDBGUI(''planclicked'')');

        %Make top Axis, table it
		ud.handles.topAxis = axes('position', [.25 .67 .7 .28]);
  		table(ud.handles.topAxis);
        table(ud.handles.topAxis, 'SETDBCALLBACK', 'planDBGUI(''fieldclicked'')');

        ud.handles.midAxis = axes('position', [.25 .55 .7 .1]);
  		table(ud.handles.midAxis);
        table(ud.handles.midAxis,'SETDBCALLBACK',  'planDBGUI(''filterclicked'')');

        %Init upper planDB controls
        uicontrol('style', 'frame', 'units', 'normalized', 'position', [.025 .55 .2 .4]);
        ud.handles.searchBox = uicontrol('style', 'edit', 'units', 'normalized', 'position', [.05 .875 .15 .03], 'callback', 'planDBGUI(''SEARCHENTERED'');', 'backgroundcolor', [1 1 1]);
        uicontrol('style', 'text', 'units', 'normalized', 'position', [.05 .84 .15 .03], 'string', 'Search Fieldnames');
        ud.handles.clearFieldFilter = uicontrol('style', 'pushbutton', 'string', 'Clear Search', 'units', 'normalized', 'position', [.05 .78 .15 .05], 'callback', 'planDBGUI(''clearFieldFilter'');');

        ud.handles.planNameSearchBox = uicontrol('style', 'edit', 'units', 'normalized', 'position', [.05 .675 .15 .03], 'callback', 'planDBGUI(''PLANNAMESEARCH'');', 'backgroundcolor', [1 1 1]);
        uicontrol('style', 'text', 'units', 'normalized', 'position', [.05 .64 .15 .03], 'string', 'Search plan names');
        ud.handles.clearPlanSearchFilter = uicontrol('style', 'pushbutton', 'string', 'Clear Search', 'units', 'normalized', 'position', [.05 .58 .15 .05], 'callback', 'planDBGUI(''clearPlanSearchFilter'');');

        %Initialize lower planDB controls
        uicontrol('style', 'frame', 'units', 'normalized', 'position', [.025 .05 .2 .4]);
        uicontrol('style', 'text', 'units', 'normalized', 'position', [.05 .375 .1 .05], 'string', 'Database:', 'horizontalalignment', 'left');

        uicontrol('style', 'text', 'units', 'normalized', 'position', [.05 .265 .1 .05], 'string', '.mat files:', 'horizontalalignment', 'left');
        uicontrol('style', 'text', 'units', 'normalized', 'position', [.05 .235 .1 .05], 'string', 'plan files:', 'horizontalalignment', 'left');
        uicontrol('style', 'text', 'units', 'normalized', 'position', [.05 .205 .1 .05], 'string', 'present plans:', 'horizontalalignment', 'left');
        uicontrol('style', 'text', 'units', 'normalized', 'position', [.05 .175 .1 .05], 'string', 'num fields:', 'horizontalalignment', 'left');

        ud.handles.planNameTx = uicontrol('style', 'text', 'units', 'normalized', 'position', [.05 .355 .1 .05], 'string', '<None>', 'horizontalalignment', 'left');
        ud.handles.matFileTx = uicontrol('style', 'text', 'units', 'normalized', 'position', [.13 .265 .07 .05], 'string', '<N/A>', 'horizontalalignment', 'right');
        ud.handles.planFileTx = uicontrol('style', 'text', 'units', 'normalized', 'position', [.13 .235 .07 .05], 'string', '<N/A>', 'horizontalalignment', 'right');
        ud.handles.presPlanTx = uicontrol('style', 'text', 'units', 'normalized', 'position', [.13 .205 .07 .05], 'string', '<N/A>', 'horizontalalignment', 'right');
        ud.handles.numFieldTx = uicontrol('style', 'text', 'units', 'normalized', 'position', [.13 .175 .07 .05], 'string', '<N/A>', 'horizontalalignment', 'right');

        ud.handles.loadDBButton = uicontrol('style', 'pushbutton', 'string', 'Load planDB', 'units', 'normalized', 'position', [.05 .075 .15 .05], 'callback', 'planDBGUI(''loadplanDB'');');

        if nargin == 2
            planDB = varargin{1};
            setappdata(h, 'planDB', planDB);
        end
        set(h, 'userdata', ud);

    case 'LOADPLANDB'
        planDB = [];
        ud = get(gcbf, 'userdata');
        [filename, pathname] = uigetfile('*.mat', 'Pick a .mat planDB file');
        if filename == 0
            return;
        else
            load(fullfile(pathname, filename));
        end

        if exist('planDB')
            %planDB = updateplandb(planDB);
            setappdata(gcbf, 'planDB', planDB);
            ud.DB.filename = filename;
            ud.DB.path = pathname;
            set(gcbf, 'userdata', ud);
            planDBGUI('UPDATEPLANSTATS');
            planDBGUI('SHOWALLPLANS');
            planDBGUI('SHOWALLFIELDS');
            planDBGUI('SHOWFILTERS');
        else
            return;
        end

    case 'UPDATEPLANSTATS'
        planDB = getappdata(gcbf, 'planDB');
        ud = get(gcbf, 'userdata');

        set(ud.handles.planNameTx, 'String', ud.DB.filename);
        set(ud.handles.matFileTx, 'String', num2str(length(planDB.matFiles)));
        set(ud.handles.planFileTx, 'String', num2str(length(find([planDB.matFiles.isPlanC]))));
        set(ud.handles.presPlanTx, 'String', num2str(length(find([planDB.matFiles.isPlanC] & [planDB.matFiles.isPresent]))));
        set(ud.handles.numFieldTx, 'String', num2str(length(planDB.fieldIndex)));
        set(gcbf, 'userdata', ud);

    case 'CLEARFIELDFILTER'
        ud = get(gcbf, 'userdata');
        set(ud.handles.searchBox, 'string', '');
        planDBGUI('SHOWALLFIELDS');

    case 'SHOWALLPLANS'
        planDB = getappdata(gcbf, 'planDB');
        ud = get(gcbf, 'userdata');
        ud.planList = [planDB.matFiles([planDB.matFiles.isPlanC]).info]; %set planlist to file info.
        set(gcbf, 'userdata', ud);
        updatePlanList(gcbf);
        numPlans = length(ud.planList);
        title(['All Plans (' num2str(numPlans) ')'], 'interpreter', 'none');

    case 'SHOWALLFIELDS'
        planDB = getappdata(gcbf, 'planDB');
        ud = get(gcbf, 'userdata');
        ud.fieldStruct = searchFieldNames('.*', planDB.fieldIndex);
        set(gcbf, 'userdata', ud);
        updateFieldList(gcbf);
        numFields = length(ud.fieldStruct);
        title(['All Fields (' num2str(numFields) ')'], 'interpreter', 'none')

    case 'SHOWFILTERS'
        ud = get(gcbf, 'userdata');
        ud.filters = struct('Filters', {}, 'fieldname', {}, 'action', {}, 'regexp', {}, 'bool', {}, 'invert', {}, 'indices', {});
        set(gcbf, 'userdata', ud);

    case 'SEARCHENTERED'
        planDB = getappdata(gcbf, 'planDB');
        ud = get(gcbf, 'userdata');
        searchString = get(gcbo, 'String');
        if isempty(searchString)
            planDBGUI('showallfields');
            return;
        end
        fieldS = searchFieldNames(searchString, planDB.fieldIndex);
        ud.fieldStruct = fieldS;
        set(gcbf, 'userdata', ud);

        updateFieldList(gcbf);

        title(['Fields matching ''' searchString ''''], 'interpreter', 'none');

        set(gcbf, 'userdata', ud);

    case 'FIELDCLICKED'
        planDBGUI('newfilter', 'AND', 'exists', '0');

    case 'PLANCLICKED'
        ud = get(gcbf, 'userdata');
        element = table(ud.handles.botAxis, 'getcurrentelement');
        planDB = getappdata(gcbf, 'planDB');

        planList = ud.planList;
        name = planList(element).name;
        path = planList(element).path;

        file = fullfile(path, name);

        CERRSliceViewer('load', file);

    case 'FILTERCLICKED'
        ud = get(gcbf, 'userdata');
        index = table(gca, 'getcurrentelement');
        ud.filters(index) = [];
        set(gcbf, 'userdata', ud);
        updateFilterList(gcbf);

    case 'NEWFILTER'
        bool = varargin{1};
        action = varargin{2};
        invert = str2num(varargin{3});

        ud = get(gcbf, 'userdata');
        index = table(ud.handles.topAxis, 'getcurrentelement');
        field = ud.fieldStruct(index).fieldname;

        ud.filters(end+1).fieldname = field;
        ud.filters(end).bool = bool;
        ud.filters(end).action = action;
        ud.filters(end).indices = [];
        if strcmpi(action, 'contains') | strcmpi(action, 'does not contain')
            regexp = inputdlg('Regexp:','Enter a regular expression');
        else
            regexp = {''};
        end

        ud.filters(end).regexp = regexp;
        ud.filters(end).Filters = [bool ' ' field ' ' action ' ' regexp{:} '.'];

        set(gcbf, 'userdata', ud);
        updateFilterList(gcbf);

end

function hMenu = rtMenu()
%Create right click menu to be passed to table.
	h = uicontextmenu;

    hAnd = 	uimenu(h, 'label', 'AND this field...');
    hOr  = 	uimenu(h, 'label', 'OR this field...');

	uimenu(hAnd, 'label', 'contains <expression>.', 'separator', 'off', 'callback', 'planDBGUI(''newfilter'', ''AND'', ''contains'', ''0'')');
	uimenu(hAnd, 'label', 'does not contain <expression>.', 'callback', 'planDBGUI(''newfilter'', ''AND'', ''contains'', ''1'')');
	uimenu(hAnd, 'label', 'exists.', 'callback', 'planDBGUI(''newfilter'', ''AND'', ''exists'', ''0'')');
	uimenu(hAnd, 'label', 'does not exist.', 'callback', 'planDBGUI(''newfilter'', ''AND'', ''exists'', ''1'')');
	uimenu(hAnd, 'label', 'is empty.', 'callback', 'planDBGUI(''newfilter'', ''AND'', ''is empty'', ''0'')');
	uimenu(hAnd, 'label', 'is not empty.', 'callback', 'planDBGUI(''newfilter'', ''AND'', ''is empty'', ''1'')');

	uimenu(hOr, 'label', 'contains <expression>.', 'separator', 'off', 'callback', 'planDBGUI(''newfilter'', ''OR'', ''contains'', ''0'')');
	uimenu(hOr, 'label', 'does not contain <expression>.', 'callback', 'planDBGUI(''newfilter'', ''OR'', ''contains'', ''1'')');
	uimenu(hOr, 'label', 'exists.', 'callback', 'planDBGUI(''newfilter'', ''OR'', ''exists'', ''0'')');
	uimenu(hOr, 'label', 'does not exist.', 'callback', 'planDBGUI(''newfilter'', ''OR'', ''exists'', ''1'')');
	uimenu(hOr, 'label', 'is empty.', 'callback', 'planDBGUI(''newfilter'', ''OR'', ''is empty'', ''0'')');
	uimenu(hOr, 'label', 'is not empty.', 'callback',  'planDBGUI(''newfilter'', ''OR'', ''is empty'', ''1'')');
    hMenu = h;

function updateFieldList(hFigure)
    ud = get(hFigure, 'userdata');
    table(ud.handles.topAxis, 'data', rmfield(ud.fieldStruct, 'planIndices'));
    table(ud.handles.topAxis, 'SETRIGHTCLICKMENU', rtMenu);

function updateFilterList(hFigure)
    ud = get(hFigure, 'userdata');
    table(ud.handles.midAxis, 'data', rmfield(ud.filters, {'fieldname', 'bool', 'action', 'regexp', 'indices', 'invert'}));
    applyFilters(hFigure);

function updatePlanList(hFigure)
    ud = get(hFigure, 'userdata');
	table(ud.handles.botAxis);
    table( ud.handles.botAxis, 'SETDBCALLBACK', 'planDBGUI(''planclicked'')');
    table(ud.handles.botAxis, 'data', rmfield(ud.planList, {'path', 'lastMod'}));
    numPlans = length(ud.planList);
    title(['Filtered Plans (' num2str(numPlans) ')'], 'interpreter', 'none');

function applyFilters(hFigure)
    ud = get(hFigure, 'userdata');
    filters = ud.filters;
    if isempty(filters)
        planDBGUI('showallplans');
    else
        ud.planList = [];
        planDB = getappdata(gcbf, 'planDB');
        allFieldNames = {planDB.fieldIndex.fieldname};

        plans = [];
        for i=1:length(filters)
            if ~isempty(ud.filters(i).indices)
                indices = ud.filters(i).indices;
            else
                fIndex = find(strcmpi(allFieldNames, filters(i).fieldname));
                fieldS = planDB.fieldIndex(fIndex);
                filters(i) = runFilter(filters(i), planDB);
            end
            ud.filters(i).indices = filters(i).indices;
            if isempty(plans), plans = filters(i).indices;, end;
            if strcmpi(filters(i).bool, 'and')
                plans = intersect(plans, filters(i).indices);
            elseif strcmpi(filters(i).bool, 'or')
                plans = union(plans, filters(i).indices);
            end
        end

        fieldsToDisplay = unique({filters.fieldname});
        ud.planList = [planDB.matFiles(plans).info];
        for i=1:length(plans)
            for j=1:length(fieldsToDisplay)
                data = getFieldContents(fieldsToDisplay{j}, planDB.matFiles(plans(i)).extract);
                ud.planList(i).(['data' num2str(j)]) = data{1};
            end
        end
        set(hFigure, 'userdata', ud);
        updatePlanList(hFigure);
    end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
classdef RWMOP17 < PROBLEM
% <multi> <real> <constrained>
% Bulk carriers design problem

%------------------------------- Reference --------------------------------
% A. Kumar, G. Wu, M. Ali, Q. Luo, R. Mallipeddi, P. Suganthan, and S. Das,
% A benchmark-suite of real-world constrained multi-objective optimization
% problems and some baseline results, Swarm and Evolutionary Computation,
% 2021, 67: 100961.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Initialization
        function Setting(obj)
            obj.M        = 3;
            obj.D        = 6;
            obj.lower    = [150.0 20.0 13.0 10.0 14.0 0.63];
            obj.upper    = [274.32 32.31 25.0 11.71 18.0 0.75];
            obj.encoding = ones(1,obj.D);
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            x   = varargin{1};
            L   = x(:,1);
            B   = x(:,2);
            D   = x(:,3);
            T   = x(:,4);
            V_k = x(:,5);
            C_B = x(:,6);

            a   = 4977.06.*C_B.^2 - 8105.61.*C_B + 4456.51;
            b   = -10847.2.*C_B.^2 + 12817.*C_B - 6960.32;
            F_n = 0.5144./(9.8065 .* L).^0.5;
            P   = ((1.025.*L.*B.*T.*C_B).^(2/3).*V_k.^3)./(a + b.*F_n);

            W_s = 0.034.*L.^1.7.*B.^0.6.*D.^0.4.*C_B.^0.5;
            W_o = L.^0.8.*B.^0.6.*D.^0.3.*C_B.^0.1;
            W_m = 0.17.*P.^0.9;
            ls  = W_s+W_o+W_m;

            D_wt  = 1.025.*L.*B.*T.*C_B-ls;
            F_c   = 0.19.*24.*P./1000 + 0.2;
            D_cwt = D_wt - F_c.*((5000.*V_k)./24 + 5)-2.*D_wt.^0.5;
            R_trp = 350./((5000.*V_k)./24 + 2.*(D_cwt./8000 + 0.5));
            ac    = D_cwt.*R_trp;
            S_d   = 5000.*V_k./24;

            C_c = 0.2.*1.3 .* (2000.*W_s.^0.85 + 3500.*W_o + 2400.*P.^0.8);
            C_r = 40000.*D_wt.^0.3;
            C_v = (1.05.*100.*F_c.*S_d + 6.3.*D_wt.^0.8).*R_trp;

            % Objectives
            f(:,1) = (C_c + C_r + C_v)./ac;
            f(:,2) = ls;
            f(:,3) = -ac;

            % Constraints
            g(:,1) = L./B - 6;
            g(:,2) = 15 - L./D;
            g(:,3) = 19 - L./T;
            g(:,4) = 0.45.*D_wt.^0.31 - T;
            g(:,5) = 0.7.*D + 0.7 - T;
            g(:,6) = 0.32 - F_n;
            g(:,7) = 0.53.*T + ((0.085.*C_B - 0.002).*B.^2)./(T.*C_B)-(1 + 0.52.*D) - 0.07.*B;
            g(:,8) = D_wt - 3000;
            g(:,9) = 500000 - D_wt;
            g      = -g;
            Population = SOLUTION(varargin{1},f,g,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
        %% Generate a point for hypervolume calculation
        function R = GetOptimum(obj,~)
            R = [-3.1514157e+03   8.2606298e+03   8.1260004e+02];
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function U = phi2u(phi)

m = length(phi);
k = ceil(sqrt(2*m));
U = eye(k);
c = cos(phi);
s = sin(phi);
count = 1;
for i=1:k
    for j=i+1:k
        R = eye(k);
        R(i,i) = c(count);
        R(i,j) = -s(count);
        R(j,i) = s(count);
        R(j,j) = c(count);
        U = U*R;
        count = count + 1;
    end
end

author: bashtage

repo: mfe-toolbox


text:
% -------------------------------------------------------------------------------------------------
function net = make_siameseFC(opts)
%MAKE_SIAMESEFC
%   Creates Siamese Fully-Convolutional network,
%   made by duplicating a vanilla AlexNet in two branches
%   and joining the branches with a cross-correlation layer
%
%   Luca Bertinetto, Jack Valmadre, Joao F. Henriques, 2016
% -------------------------------------------------------------------------------------------------------
    net_opts.siamese = true;
    net_opts.batchNormalization = true;
    net_opts.strides = [2, 2, 1, 2];
    net_opts = vl_argparse(net_opts, {opts.net.conf});

    if opts.pretrain
        error('not yet implemented for alexnet');
    end

    % create siamese branch of the network (5 conv layers)
    fprintf('construct network\n');
    % four different net sizes to tradeoff accuracy/speed: vid_create_net, _small1, _small2 and _small3
    branch = vid_create_net(...
        'exemplarSize',       opts.exemplarSize * [1 1], ...
        'instanceSize',       opts.instanceSize * [1 1], ...
        'batchNormalization', net_opts.batchNormalization, ...
        'networkType',        'simplenn', ...
        'weightInitMethod',   opts.init.weightInitMethod, ...
        'scale',              opts.init.scale, ...
        'initBias',           opts.init.initBias, ...
        'strides',            net_opts.strides);

    branch = dagnn.DagNN.fromSimpleNN(branch);

    % Add common final stream of the network
    net = dagnn.DagNN();
    rename_io_vars(branch, 'in', 'out');
    add_pair_of_streams(net, branch, ...
                        {'exemplar', 'instance'}, ...
                        {'a_feat', 'b_feat'}, ...
                        net_opts.siamese);

    net.addLayer('xcorr', XCorr(), ...
                 {'a_feat', 'b_feat'}, ...
                 {'xcorr_out'}, ...
                 {});

    add_adjust_layer(net, 'adjust', 'xcorr_out', 'score', ...
                 {'adjust_f', 'adjust_b'}, 1e-3, 0, 0, 1);

end

author: bertinetto

repo: siamese-fc


text:
function [TorF, vstr, rdate] = have_feature_gurobi()
%HAVE_FEATURE_GUROBI  Detect availability/version info for Gurobi
%
%   Feature detection function implementing 'gurobi' tag for HAVE_FEATURE
%   to detect availability/version of Gurobi optimizer
%   (https://www.gurobi.com).
%
%   See also HAVE_FEATURE, QPS_MASTER, MIQPS_MASTER, GUROBI.

%   MP-Opt-Model
%   Copyright (c) 2004-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

TorF = exist('gurobi', 'file') == 3;
vstr = '';
rdate = '';
if TorF
    try
        model = struct( ...
            'A', sparse(1), ...
            'rhs', 1, ...
            'sense', '=', ...
            'vtype', 'C', ...
            'obj', 1, ...
            'modelsense', 'min' ...
        );
        params = struct( ...
            'outputflag', 0 ...
        );
        result = gurobi(model, params);
        vstr = sprintf('%d.%d.%d', result.versioninfo.major, result.versioninfo.minor, result.versioninfo.technical);
    catch % gurobiError
        TorF = 0;
        fprintf('Gurobi Error!\n');
%         disp(gurobiError.message);
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%% Estimate the quality of QRS complex detections
% Calculate the quality of QRS detections based on the the likelihood of
% measurements from RR interval series with respect to a trained model. See
% for further details. See 'A Pattern-Recognition Approach for
% Lead-Selection in Heartbeat Detection' CINC 2014.
%
% Example
%
%   [ ratio, estimated_labs ] = CalcRRserieRatio(time_serie, ECG_header, start_end, ECG_type)
%
% Arguments:
%      +time_serie: [cell] REQUIRED
%
%           Cell array of size [nsig 1] with the detections
%           performed for each lead.
%
%     +ECG_header: [struct] OPTIONAL.
%
%             Description of the ECG typically available in the
%             ECG_header. Structure with fields:
%
%               -freq: Sampling rate in Hz. (1)
%
%               -nsig: Number of ECG leads. (size(ECG,2))
%
%               -nsamp: Number of ECG samples. (size(ECG,1))
%
%     +start_end: [numeric] OPTIONAL. Vector of size [2 1] with the start
%                           and end samples.
%
%     +ECG_type: [numeric] OPTIONAL. Char indicating if you know something
%                          about the ECG recording type a priori. Possible
%                          values are: 'stress' 'arrhythmia' 'longterm' 'sinus'
%
% Output:
%     + ratio: a measure of quality of the QRS detections
%
%     + estimated_labs: a label for each heartbeat: {TP, FP and FN}
%
% See also ECGtask_QRS_detection
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 23/4/2013
% Copyright 2008-2015
function [ ratio, estimated_labs ] = CalcRRserieRatio(time_serie, ECG_header, start_end, ECG_type)

    %% constants

    min_pattern_separation = 350;  % ms
    max_pattern_separation = 1500; % ms
    cECGtypes = {'stress' 'arrhythmia' 'longterm' 'sinus'};

    %% start

    if( nargin < 2 || isempty(time_serie) || isempty(ECG_header) )
        ratio = 0;
        return
    end

    if( nargin < 3 || isempty(start_end) )
        start_end = [1 ECG_header.nsamp];
    end

    if( nargin < 4 || ~any(strcmpi(ECG_type, cECGtypes)) )
        ECG_type = 'stress';
    end

    start_sample = start_end(1);
    end_sample = start_end(2);

    pb = progress_bar( 'Calculating ratios', 'Start' );

    % co_ocurrences respect other leads
    co_ocurrence = calc_co_ocurrences(time_serie);

    if( isempty(co_ocurrence) )
        return
    end

    dummy = prdataset([]);
    dummy = prmapping([]);

    % load the trained classifier.
    aux_load = load([ 'QRS_q_' ECG_type '.mat']);

    w_lablist = cellstr(getlabels(aux_load.wTrained_Classifier));
    FN_lab = find(strcmpi(w_lablist, 'FN'));
    FP_lab = find(strcmpi(w_lablist, 'FP'));
    TP_lab = find(strcmpi(w_lablist, 'TP'));

    lreferences = length(time_serie);

    pb.Loops2Do = lreferences;
    pb.checkpoint('Calculating ratios.');

    ratio = nan(lreferences,1);

    estimated_labs = cell(lreferences,1);

    for ii = 1:lreferences

        pb.start_loop();

        this_ts = colvec(time_serie{ii});

        if( length(this_ts) < (ECG_header.nsamp / max_pattern_separation) )
%             warning( [mfilename ':few_anns'], 'Few detections in  %d\n', ii );
            ratio(ii) = 0;
        else

            pb.checkpoint([]);

            k_gaps = calc_k_gaps(this_ts);

            % build the feature matrix [ RR_i-1 RR_i RR_10 RR_60 co_ocurrences]
            % RR_i-1
            this_ts = round(this_ts * 1000 / ECG_header.freq);

            pb.checkpoint([]);

            aux_fm = [ calculate_RR_features(this_ts) colvec(co_ocurrence{ii}) ];

            pb.checkpoint([]);

            ds_result = prdataset(aux_fm) * aux_load.wTrained_Classifier;

            estimated_labs{ii} = renumlab(ds_result * labeld, char(w_lablist) );

            aux_val = sum(estimated_labs{ii} == TP_lab | estimated_labs{ii} == FN_lab);
            if( aux_val == 0 )
                this_se = 0;
            else
                this_se = sum(estimated_labs{ii} == TP_lab) / aux_val;
            end

            aux_val = sum(estimated_labs{ii} == TP_lab | estimated_labs{ii} == FP_lab);
            if( aux_val == 0 )
                this_pp = 0;
            else
                this_pp = sum(estimated_labs{ii} == TP_lab) / aux_val;
            end

            % metric used for the paper experiment
            %this_q = (2*this_se + this_pp)/3;
            % F1 score
            this_q = 2/(1/this_se + 1/(eps+this_pp));
            ratio(ii) = k_gaps * this_q;

        end

        pb.end_loop();
    end

%     clear pb

    function k = calc_k_gaps( time_serie )

        time_serie = time_serie( time_serie >= start_sample & time_serie <= end_sample );

        % add some samples in order to calculate the gaps among heartbeats (FN)
        if( time_serie(1) ~= start_sample )
            time_serie = [start_sample; time_serie];
        end

        if( time_serie(end) ~= end_sample )
            time_serie = [time_serie; end_sample];
        end

        time_serie = round(time_serie * 1000 / ECG_header.freq);

        RRserie = colvec(diff(time_serie));
        RRserie = [0;RRserie];

        gap_idx = find(RRserie > max_pattern_separation);
%         gap_end_idx = gap_idx;
%         gap_start_idx = gap_start_idx - 1;

        % Percent of the time we have gaps.
        k = 1 - (sum(RRserie(gap_idx)) / time_serie(end));

    end

    function fm = calculate_RR_features( time_serie )

        RRserie = diff(time_serie);
        RRserie = [RRserie(1); RRserie];

        % build the feature matrix [ RR_i-1 RR_i RR_10 RR_60 ]
        % RR_i-1
        fm = [ RRserie(1); RRserie(1:end-1) ];
        fm = [fm RRserie CalcFeatureRRx(time_serie, RRserie, 10000) CalcFeatureRRx(time_serie, RRserie, 60000)];

    end

    function start_end_aux = findStartEnd( bAux )

        start_aux = find(  bAux, 1, 'first' );
        end_aux = find(  bAux, 1, 'last' );
        start_end_aux = [start_aux end_aux];

    end

    function aux_mean = CalcFeatureRRx(time_serie, RRserie, win_size)
%         win_size in milliseconds

        aux_seq = 1:length(RRserie);
        aux_idx = arrayfun(@(a)( findStartEnd(  time_serie >= (time_serie(a) - win_size) & ...
                                                time_serie <= time_serie(a) )), ...
                           aux_seq, 'UniformOutput', false);

        aux_mean = colvec(cellfun(@(a)(round(median(RRserie(a(1):a(2))))), aux_idx));

    end

end


text:
function Y = expvarMeanOut(model, X);

% EXPVARMEANOUT Output of an EXPVARMEAN model.
% FORMAT
% DESC gives the output of the mean of a variational approximation
% to a given kernel. This is for use with the exponentiated kernel
% when making a variational approximation to an exponentiated
% Gaussian process.
% ARG model : the model for which the output is required.
% ARG X : the input data for which the output is required.
% RETURN Y : the output.
%
%
% SEEALSO :  modelOut, expvarMeanCreate
%
% COPYRIGHT : Neil D. Lawrence, 2006

% SHEFFIELDML

Y = exp(0.5*kernDiagCompute(model.kern.argument, X));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function e = legendre_monomial_quadrature ( n, x, w, p )

%*****************************************************************************80
%
%% LEGENDRE_MONOMIAL_QUADRATURE applies a quadrature rule to a monomial.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points in the rule.
%
%    Input, real X(N), the quadrature points.
%
%    Input, real W(N), the quadrature weights.
%
%    Input, integer P, the exponent.
%
%    Output, real E, the quadrature error.
%

%
%  Get the exact value of the integral.
%
  t = legendre_integral ( p );
%
%  Evaluate the monomial at the quadrature points.
%
  v(1:n,1) = x(1:n).^p;
%
%  Compute the weighted sum.
%
  q = w' * v;
%
%  Error:
%
  if ( t == 0.0 )
    e = abs ( q - t );
  else
    e = abs ( ( q - t ) / t );
  end

  return
end


text:
function [F,h,x] = sedumi2yalmip(At,b,c,K)

nvars = length(b);
x = sdpvar(nvars,1);

if size(At,2)~=length(b)
    At = At';
end

F = ([]);
top = 1;

if isvalidfield(K,'f')
    X = c(top:top+K.f-1)-At(top:top+K.f-1,:)*x;
    F = F + (X(:) == 0);
    top = top + K.f;
end

if isvalidfield(K,'l')
    X = c(top:top+K.l-1)-At(top:top+K.l-1,:)*x;
    F = F + (X(:)>=0);
    top = top + K.l;
end

if isvalidfield(K,'q')
    for i = 1:length(K.q)
        X = c(top:top+K.q(i)-1)-At(top:top+K.q(i)-1,:)*x;
        F = F + (cone(X(2:end),X(1)));
        top = top + K.q(i);
    end
end

if isvalidfield(K,'r')
    for i = 1:length(K.r)
        X = c(top:top+K.r(i)-1)-At(top:top+K.r(i)-1,:)*x;
        F = F + (rcone(X(3:end),X(2),X(1)));
        top = top + K.r(i);
    end
end

if isvalidfield(K,'s')
    for i = 1:length(K.s)
        [ix,iy,iv] = find([c(top:top+K.s(i)^2-1) At(top:top+K.s(i)^2-1,:)]);
        off = (ix-1)/(K.s(i)+1);
        if all(off == round(off))
            X = c(top:top+K.s(i)^2-1)-At(top:top+K.s(i)^2-1,:)*x;
            if isa(X,'sdpvar')
                F = F + (diag(reshape(X,K.s(i),K.s(i))) >= 0);
            else
                X
                i
                'silly data!'
            end
            top = top + K.s(i)^2;
        else
            X = c(top:top+K.s(i)^2-1)-At(top:top+K.s(i)^2-1,:)*x;
            X = reshape(X,K.s(i),K.s(i));
            X = (X+X')/2;
            F = F + (X >= 0);
            top = top + K.s(i)^2;
        end
    end
end

h = -b'*x;

function ok = isvalidfield(K,fld)
ok = 0;
if isfield(K,fld)
    s = getfield(K,fld);
    if prod(size(s))>0
        if s(1)>0
            ok = 1;
        end
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
 function [ob,m,mz] = Gocrf(k,fov,N,mdomswitch,dt,varargin)
%function [ob,m,mz] = Gocrf(k,fov,N,m0,mdomswitch,dt,varargin)
%
% Construct Gocrf object for fast optimal control large-tip-angle
% multidimensional (and parallel) MR RF pulse design. Used in
% conjunction with an optimization routine, this object can be
% used to design small pulse updates that, when summed with an
% underlying pulse, result in an accurate excitation.
% The vectors m and mz are the results of initial Bloch sim at
% object creation.
%
% The object is created by calling:
%    Gocrf = Gocrf( ... );
% and it is then used by commands like:
%    y = Gocrf * x;
% which will evaluate perturbations to the magnetization
% produced by pulse perturbations x. It also passes x with
% unit weighting, so that if G is the sub-object relating
% x to its perturbations, then y = [G;I]*x
%
% in:
%    k            [Nt,d]   excitation k-space trajectory
%                          (= -gam int_t^T g(t') dt')
%                          (inverse fov units)
%    fov          [d,1]    field of view in each dimension
%    N            [d,1]    size of design grid in each
%                          dimension
%    mdomswitch   [1]      switch to evaluate magnetization
%                          or spinors (spinors not tested yet)
%    dt           [1]      RF sampling period
%
% options:
%    m0           [3,1]    initial magnetization state
%                          (e.g., if equilibrium, [0 0 1] (default))
%    'sens'       [prod(N),Ncoils]   transmit (B1+) sensitivities
%    'd'          [[N]]    desired excitation pattern, in
%                          arbitrary flip angle units. this is used
%                          to determine sub-space over which to
%                          simulate
%    'indmask'    [[N]]    mask to choose sub-space points from
%    'baseB1'     [Ncoils*Nt]    baseline RF pulse(s), used for
%                                initial Bloch simulation during
%                                object creation
%    'g'          [Nt,d]   gradient waveforms, for Bloch simulation
%                          (gauss/fov units)
%    'Nsvd'       [1]      # of baseline excitation expansion terms
%                          (default is 4)
%    'a0B','b0B'  [Nt,Nsvd]      baseline temporal basis function
%                                matrices
%    'a0Ct','b0Ct' [prod(N),Nsvd]   baseline spatial coefficient
%                                   matrices
%    'a0','b0'    [prod(N),1]    final baseline spinors
%    'nufft'      cell     nufft arguments (passed to Gmri_SENSE)
%    'RFbasis'    [Nt,Nbasis]    RF pulse basis functions for
%                                parameterized pulse design
%
% options for field-corrected pulse design (untested in this object!):
%    'ti'         [Nt,1]   sample times (0:T in seconds)
%    'we'         [[N]]    field_map (Hz)
%    'L'          [1]      # of approximation terms (default is 4)
%
% After building this object, the user can Bloch-simulate
% a new pulse and update the object using:
%     [ob,m,mz] = feval(ob.arg.update, ob, baseB1);
%
% or a Bloch simulation (without object update) can be performed
% using:
%     [m,mz,a0,b0] = feval(ob.arg.blochsim,baseB1,arg,pos,sens,we);
%
% OR an update can be made (without Bloch sim) using:
%     ob = feval(ob.arg.new_B_Ct,ob,a0B,a0Ct,b0B,b0Ct,a0,b0)
%
% Copyright 2008-10-8, Will Grissom, Stanford University

if nargin < 6, help(mfilename), error(mfilename), end

%
% default object
%
% Gmri objects
arg.a0G = [];arg.b0G = [];
% svd matrices
arg.a0B = [];arg.a0Ct = [];arg.b0B = [];arg.b0Ct = [];
% underlying final alpha0, beta0
arg.a0 = [];arg.b0 = [];
% initial mag (global)
arg.bsq = 0; %
arg.exc = 0; % just return mxy
arg.rfb = 0; % no rf basis
arg.m0 = [0 0 1];
arg.mdomswitch = mdomswitch;
arg.we = [];
arg.nufft_args = {N,6*ones(size(N)),2*N,N/2,'table',2^10,'minmax:kb'};
arg.new_B_Ct = @Gocrf_new_B_Ct;
arg.blochsim = @Gocrf_blochsim;
arg.subindcalc = @Gocrf_subindcalc;
arg.update = @Gocrf_update;
arg.subind = []; % set of points to Bloch sim and do SVD on
arg.nsubbins = 10; % default # of histogram bins to choose sub-points from
arg.simpos = []; % spatial locations for Bloch sim
arg.sens = []; % need sens for Bloch sims
arg.d = []; % desired pattern (for subind calc)
arg.indmask = []; % mask (for subind calc)
arg.simmask = []; % mask for Bloch simulation
arg.dt = dt; % sampling period
arg.ti = []; % time points (for field-map corrected design)
arg.fov = fov;
arg.k = k; % excitation k-space traj
arg.N = N; % spatial dimensions
arg.L = 4;
arg.Nsvd = 4;
arg.subindlen = 15*arg.Nsvd; % default # of sub-points for svd sim
% field map bases and coeffs
arg.B = [];arg.Ct = [];
arg.baseB1 = []; % base B1 pulses (for Bloch sim)
arg.g = []; % gradient waveforms (for Bloch sim)
arg.fmod = []; % modulation freq (for eg dual band designs)
arg.aux1 = 1; % vector to multiply a or mxy by
arg.aux2 = 1; % vector to multiply b or b^2 or mz by
arg.RFbasis = [];
ob.Nbasis = 0;
% Bloch-simulated patterns (will be empty if basis and coeff
% matrices are provided)
m = [];mz = [];

gambar = 4257;             % gamma/2pi in Hz/g
gam = gambar*2*pi;         % gamma in radians/g

arg = vararg_pair(arg, varargin, 'subs', ...
	{'basis', 'basis_args'; 'nufft', 'nufft_args'});

% if not provided, get the subindices for Bloch sim
if isempty(arg.subind) & ~isempty(arg.d) & ~isempty(arg.indmask)
  arg = Gocrf_subindcalc(arg,arg.d,arg.indmask);
elseif isempty(arg.subind)
  disp ['No sub indices or means to calculate them provided. update() ' ...
        'will not work'];
end

ktmp = -[arg.k - repmat(arg.k(1,:)/2,[size(arg.k,1) 1])];
if isempty(arg.we) || sum(abs(arg.we(:))) == 0  % no field map

	if ~isempty(arg.sens)
		arg.a0G = Gmri_SENSE(ktmp,true(N),'fov',fov, ...
			'basis',{'dirac'},'nufft',arg.nufft_args, ...
			'exact',0,'sens',arg.sens*(-1i*gam*arg.dt/2))';
		arg.b0G = Gmri_SENSE(ktmp,true(N),'fov',fov, ...
			'basis',{'dirac'},'nufft',arg.nufft_args, ...
			'exact',0,'sens',arg.sens*(1i*gam*arg.dt/2))';
	else
%{
		arg.a0G = Gmri_SENSE(ktmp,true(N),'fov',fov,'basis',{'dirac'},'nufft',arg.nufft_args, ...
                         'exact',0,'sens',ones(prod(N),1)*(-1i*gam*arg.dt/2))';
		arg.b0G = Gmri_SENSE(ktmp,true(N),'fov',fov,'basis',{'dirac'},'nufft',arg.nufft_args, ...
                         'exact',0,'sens',ones(prod(N),1)*(1i*gam*arg.dt/2))';
%}
		% jf version:
		arg.a0G = (-1i*gam*arg.dt/2)' * ...
			Gmri(ktmp, true(N), 'fov', fov, 'basis', {'dirac'}, ...
				'nufft',arg.nufft_args)';
		arg.b0G = (1i*gam*arg.dt/2)' * ...
			Gmri(ktmp, true(N), 'fov', fov, 'basis', {'dirac'}, ...
				'nufft', arg.nufft_args)';

	end

else % non-zero field map

	if ~isempty(arg.sens)
		arg.a0G = Gmri_SENSE(ktmp,true(N),'fov',fov,...
			'basis',{'dirac'},'nufft',arg.nufft_args, ...
			'exact',0,'sens',arg.sens*(-1i*gam*arg.dt/2), ...
			'ti',-arg.ti+arg.ti(end)/2,...
			'zmap',1i*2*pi/2*arg.we,'L',arg.L)';
		arg.b0G = Gmri_SENSE(ktmp,true(N),'fov',fov,...
			'basis',{'dirac'},'nufft',arg.nufft_args, ...
			'exact',0,'sens',arg.sens*(1i*gam*arg.dt/2), ...
			'ti',-arg.ti+arg.ti(end)/2, ...
			'zmap',1i*2*pi/2*arg.we,'L',arg.L)';
	else
%{
		arg.a0G = Gmri_SENSE(ktmp,true(N),'fov',fov,...
			'basis',{'dirac'},'nufft',arg.nufft_args, ...
			'exact',0,'sens',ones(prod(N),1)*(-1i*gam*arg.dt/2), ...
			'ti',-arg.ti+arg.ti(end)/2, ...
			'zmap',1i*2*pi/2*arg.we,'L',arg.L)';
		arg.b0G = Gmri_SENSE(ktmp,true(N),'fov',fov, ...
			'basis',{'dirac'},'nufft',arg.nufft_args, ...
			'exact',0,'sens',ones(prod(N),1)*(1i*gam*arg.dt/2), ...
 			'ti',-arg.ti+arg.ti(end)/2, ...
			'zmap',1i*2*pi/2*arg.we,'L',arg.L)';
%}
		% jf version:
		arg.a0G =(-1i*gam*arg.dt/2)' * ...
			Gmri(ktmp,true(N),'fov',fov,...
				'basis',{'dirac'},'nufft',arg.nufft_args, ...
				'ti',-arg.ti+arg.ti(end)/2, ...
				'zmap',1i*2*pi/2*arg.we,'L',arg.L)';
		arg.b0G = (1i*gam*arg.dt/2)' * ...
			Gmri(ktmp,true(N),'fov',fov, ...
				'basis',{'dirac'},'nufft',arg.nufft_args, ...
 				'ti',-arg.ti+arg.ti(end)/2, ...
				'zmap',1i*2*pi/2*arg.we,'L',arg.L)';
	end

	arg.B = arg.a0G.arg.B;
	arg.Ct = arg.a0G.arg.Ct;
end

tmp = arg.a0G;
[dim1,dim2] = size(tmp);
if arg.exc
  arg.dim = [dim1 + 2*dim2, 2*dim2];
else
  arg.dim = [2*dim1+2*dim2, 2*dim2];
end

if ~isempty(arg.RFbasis)

  arg.rfb = 1;
  arg.Nbasis = size(arg.RFbasis,2);

  % assuming that there is only one coil when we parameterize the
  % pulse, adjust object dims
  if arg.exc
    arg.dim = [dim1 + 2*arg.Nbasis, 2*arg.Nbasis];
  else
    arg.dim = [2*dim1 + 2*arg.Nbasis, 2*arg.Nbasis];
  end

  % allocate space for basis perturbations
  arg.Gat = zeros(dim1,arg.Nbasis);
  arg.Gbt = zeros(dim1,arg.Nbasis);

end

arg.is.empty = false;

% now build object
%ob = Fatrix(arg.dim, arg, 'caller', mfilename, ...
%	'forw', @Gocrf_forw, 'back', @Gocrf_back);
ob = fatrix2('odim', arg.dim(1), 'idim', arg.dim(2), ...
	'arg', arg, ...
	'forw', @Gocrf_forw, 'back', @Gocrf_back); % jf

% if SVD matrices not supplied, do the necessary Bloch sims
if isempty(ob.arg.a0B)
  % get Bloch simulation grid
  ndgridarg = ['-fov(1)/2:fov(1)/N(1):fov(1)/2-fov(1)/N(1)'];
  posarg = ['x1'];
  for ii = 2:length(N)
    ndgridarg = [ndgridarg sprintf(',-fov(%d)/2:fov(%d)/N(%d):fov(%d)/2-fov(%d)/N(%d)',ii,ii,ii,ii,ii,ii)];
    posarg = [posarg sprintf(',x%d',ii)];
  end
  eval(sprintf('[%s] = ndgrid(%s);',posarg,ndgridarg));
  ob.arg.simpos = [x1(:)];
  for ii = 2:length(N)
    ob.arg.simpos = [ob.arg.simpos eval(sprintf('x%d(:)',ii))];
  end
  % do the Bloch sims and expansion calc
  [ob,m,mz] = Gocrf_update(ob,ob.arg.baseB1);
  % we don't need base B1 anymore
  ob.arg.baseB1 = [];
end

%
% Gocrf_new_B_Ct()
% new values of (svd) bases and coefficients
%
function G = Gocrf_new_B_Ct(G, a0B, a0Ct, b0B, b0Ct, a0, b0)

% apply frequency modulation to basis matrices
if ~isempty(G.arg.fmod)
  t = [0:G.arg.dt:(size(G.arg.k,1)-1)*G.arg.dt] - size(G.arg.k,1)*G.arg.dt/2;
  Amod = repmat(exp(1i*2*pi*G.arg.fmod*t(:)),[1 G.arg.Nsvd]);
  a0B = a0B .* Amod;
  b0B = b0B .* Amod;
end

if ~isempty(G.arg.we) & sum(abs(G.arg.we(:))) ~= 0
  % multiply out new expansion with field map expansion
  newa0B = [];newa0Ct = [];
  newb0B = [];newb0Ct = [];
  for ll = 1:G.arg.L
    newa0B = [newa0B spdiag(G.arg.B(:,ll))*a0B];
    newa0Ct = [newa0Ct spdiag(G.arg.Ct(:,ll))*a0Ct];
    newb0B = [newb0B spdiag(G.arg.B(:,ll))*b0B];
    newb0Ct = [newb0Ct spdiag(G.arg.Ct(:,ll))*b0Ct];
  end
  % update G with provided B and C.' matrices
  G.arg.a0G = feval(G.arg.a0G.arg.new_B_Ct,G.arg.a0G,newa0B,newa0Ct);
  G.arg.b0G = feval(G.arg.b0G.arg.new_B_Ct,G.arg.b0G,newb0B,newb0Ct);
else
  % set the expansion directly
  G.arg.a0G = feval(G.arg.a0G.arg.new_B_Ct,G.arg.a0G,a0B,a0Ct);
  G.arg.b0G = feval(G.arg.b0G.arg.new_B_Ct,G.arg.b0G,b0B,b0Ct);
end

G.arg.a0 = a0;
G.arg.b0 = b0;

%
% Gocrf_forw(): y = G * x
%
function v = Gocrf_forw(arg, u)

% assume the RF vector is [real(db1); imag(db1)]
db1 = u(1:length(u)/2) + 1i*u(length(u)/2+1:end);

if ~arg.rfb
  if any(db1)
    at = arg.b0G * conj(db1);
    bt = arg.a0G * conj(db1);
  else
    at = zeros(size(arg.b0G,1),1);
    bt = zeros(size(arg.a0G,1),1);
  end
else % pulse is parameterized
  if any(db1)
    at = arg.Gat * conj(db1);
    bt = arg.Gbt * conj(db1);
  else
    at = zeros(size(arg.Gat,1),1);
    bt = zeros(size(arg.Gbt,1),1);
  end
end

if arg.mdomswitch
  % dm terms, equilibrium
  dm = 0;dmz = 0;
  if arg.m0(3) ~= 0
    % dm terms, equilibrium
    dm = -2*arg.m0(3)*conj(arg.a0(:).*bt + at.*arg.b0(:));
    % dmz, equilibrium
    dmz = 2*arg.m0(3)*real(arg.a0(:).*conj(at)) - 2*arg.m0(3)*real(arg.b0(:).*conj(bt));
  end
  if any(arg.m0(1:2) ~= 0)
    % dm terms, non-equilibrium
    dm = dm - (arg.m0(1)-1i*arg.m0(2))*conj(2*arg.b0(:).*bt);
    dm = dm + (arg.m0(1)+1i*arg.m0(2))*conj(2*arg.a0(:).*at);
    % dmz, non-equilibrium
    dmz = dmz + 2*real((arg.m0(1)+1i*arg.m0(2))*(conj(arg.a0(:)).*bt + conj(at).*arg.b0(:)));
  end
  if ~arg.exc
    v = [arg.aux1.*dm;arg.aux2.*dmz;u];
  else
    v = [arg.aux1.*dm;u];
  end
else
  if ~arg.bsq
    v = [arg.aux1.*at;arg.aux2.*bt;u];
  else
    v = [2*arg.aux1.*arg.a0(:).*at;2*arg.aux2.*arg.b0(:).*bt;u];
  end
end

%
% Gocrf_back(): x = G' * y
%
function v = Gocrf_back(arg,u)

v = zeros(size(arg.a0G,2),1);
Ns = size(arg.a0G,1);
tmp1 = 0; tmp2 = 0;

% tmp1 contains terms that will multiply by system matrix for at
if arg.mdomswitch
  if arg.m0(3) ~= 0
    if ~arg.exc
      tmp1 = arg.m0(3)*2*conj(-arg.b0(:).*u(1:Ns)+conj(arg.a0(:)).*u(Ns+1:2*Ns));
    else
      tmp1 = arg.m0(3)*2*conj(-arg.b0(:).*u(1:Ns));
    end
  end
  if any(arg.m0(1:2) ~= 0)
    % additional alphat terms from Mxy, if non-equilibrium ic
    tmp1 = tmp1 + 2*conj((arg.m0(1)-1i*arg.m0(2))*arg.a0(:).*u(1:Ns));
    % additional alphat terms from Mz, if non-equilibrium ic
    if ~arg.exc
      tmp1 = tmp1 + 2*(arg.m0(1)+1i*arg.m0(2))*arg.b0(:).*conj(u(Ns+1:2*Ns));
    end
  end
  % tmp2 contains terms that will multiply by system matrix for betat
  if arg.m0(3) ~= 0
    if ~arg.exc
      tmp2 = -arg.m0(3)*2*conj(arg.a0(:).*u(1:Ns)+conj(arg.b0(:)).*u(Ns+1:2*Ns));
    else
      tmp2 = -arg.m0(3)*2*conj(arg.a0(:).*u(1:Ns));
    end
  end
  if any(arg.m0(1:2) ~= 0)
    % additional betat terms from Mxy, if non-equilibrium ic
    tmp2 = tmp2 - 2*conj((arg.m0(1)+1i*arg.m0(2))*arg.b0(:).*u(1:Ns));
    % additional betat terms from Mz, if non-equilibrium ic
    if ~arg.exc
      tmp2 = tmp2 + 2*conj((arg.m0(1)+1i*arg.m0(2))*conj(arg.a0(:)).*u(Ns+1:2*Ns));
    end
  end
else

  if ~arg.bsq
    tmp1 = u(1:Ns);
    tmp2 = u(Ns+1:2*Ns);
  else
    tmp1 = 2*conj(arg.a0(:)).*u(1:Ns);
    tmp2 = 2*conj(arg.b0(:)).*u(Ns+1:2*Ns);
  end
end

if ~arg.rfb
  v = conj(arg.a0G' * (conj(arg.aux2).*tmp2));
  v = v + conj(arg.b0G' * (conj(arg.aux1).*tmp1));
else
  v = conj(arg.Gbt' * (conj(arg.aux2).*tmp2));
  v = v + conj(arg.Gat' * (conj(arg.aux1).*tmp1));
end

if ~arg.exc
  v = [real(v);imag(v)] + u(2*Ns+1:end);
else
  v = [real(v);imag(v)] + u(Ns+1:end);
end

%
% Gocrf_subindcalc(): Determine subset of points to simulate for
%                     SVD evaluation
function arg = Gocrf_subindcalc(arg,d,indmask)

% get bin centers
if min(d(:)) == max(d(:))
  fa = min(d(:))*ones(arg.nsubbins+1,1);
else
  fa = min(d(:)):(max(d(:))-min(d(:)))/arg.nsubbins:max(d(:));
end

dind = 1:prod(size(d));

% arrange indices into bins
indices{1} = dind(col(logical(indmask.*(d >= fa(1) & d <= fa(2)))));
for ii = 2:length(fa)-1
  indices{ii} = dind(col(logical(indmask.*(d > fa(ii) & d <= fa(ii+1)))));
end

% loop through bins again and again, pulling one from each bin
% each time around at random (use rand*size(bin)), until
% index vector is full. This works much better than my
% previous method, which failed for filter-based desired patterns
indvec = [];
indind = 0;
rand('twister',1000);
while length(indvec) < arg.subindlen
  if ~isempty(indices{indind+1})
    binind = ceil(length(indices{indind+1})*rand(1,1));
    indvec = [indvec indices{indind+1}(binind)];
    indices{indind+1} = [indices{indind+1}(1:binind-1) ...
                        indices{indind+1}(binind+1:end)];
  end
  indind = mod(indind + 1,arg.nsubbins);
end

arg.subind = indvec;

%
% Gocrf_update(): Perform Bloch sims to update SVD expansions
%
function [G,m,mz] = Gocrf_update(G,baseB1)

% perform reduced Bloch sim
if isfield(G.arg, 'sens') && ~isempty(G.arg.sens) % jf

  if ~isempty(G.arg.we) & sum(abs(G.arg.we(:))) ~= 0
    [foo,fooz,alpha0,beta0] = Gocrf_blochsim(baseB1,G.arg,G.arg.simpos(G.arg.subind,:), ...
                                             G.arg.sens(G.arg.subind,:),G.arg.we(G.arg.subind));
  else
    [foo,fooz,alpha0,beta0] = Gocrf_blochsim(baseB1,G.arg,G.arg.simpos(G.arg.subind,:), ...
                                             G.arg.sens(G.arg.subind,:),G.arg.we);
  end

else

  if ~isempty(G.arg.we) & sum(abs(G.arg.we(:))) ~= 0
    [foo,fooz,alpha0,beta0] = Gocrf_blochsim(baseB1,G.arg,G.arg.simpos(G.arg.subind,:), ...
                                             ones(size(G.arg.subind)),G.arg.we(G.arg.subind));
  else
    [foo,fooz,alpha0,beta0] = Gocrf_blochsim(baseB1,G.arg,G.arg.simpos(G.arg.subind,:), ...
                                             ones(size(G.arg.subind)),G.arg.we);
  end

end

% calculate rotating frame transformation
ktmp = -[G.arg.k - repmat(G.arg.k(1,:),[size(G.arg.k,1) 1])];
Asub = exp(1i*2*pi/2*G.arg.simpos(G.arg.subind,:)*ktmp');
if ~isempty(G.arg.fmod)
  t = 0:G.arg.dt:G.arg.dt*(size(G.arg.k,1)-1);
  Asub = Asub.*exp(-1i*2*pi/2*repmat(G.arg.fmod*t,[length(G.arg.subind) 1]));
end

if ~isempty(G.arg.we) & sum(abs(G.arg.we(:))) ~= 0
  t = 0:G.arg.dt:G.arg.dt*(size(G.arg.k,1)-1);
  Asub = Asub.*exp(-1i*2*pi/2*G.arg.we(G.arg.subind).'*t);
end

% get temporal basis funcs via svd
[foo1,foo2,G.arg.a0B] = svd(conj(alpha0)./Asub,'econ');
[foo1,foo2,G.arg.b0B] = svd(conj(beta0)./Asub,'econ');

% truncate basis
G.arg.a0B = G.arg.a0B(:,1:G.arg.Nsvd);
G.arg.b0B = G.arg.b0B(:,1:G.arg.Nsvd);

% re-simulate to get coeffs at all spatial locs
if ~isempty(G.arg.sens)
  [m,mz,G.arg.a0,G.arg.b0,G.arg.a0Ct,G.arg.b0Ct] = Gocrf_blochsim(baseB1,G.arg,G.arg.simpos,G.arg.sens,G.arg.we,G.arg.simmask);
else
  [m,mz,G.arg.a0,G.arg.b0,G.arg.a0Ct,G.arg.b0Ct] = Gocrf_blochsim(baseB1,G.arg,G.arg.simpos,ones(size(G.arg.simpos,1),1),G.arg.we,G.arg.simmask);
end

m = reshape(m,G.arg.N);
mz = reshape(mz,G.arg.N);
G.arg.a0 = reshape(G.arg.a0,G.arg.N);
G.arg.b0 = reshape(G.arg.b0,G.arg.N);

% recalculate subindices based on Bloch sim results, in
% anticipation of the next update
%G.arg = Gocrf_subindcalc(G.arg,real(asin(abs(m))),G.arg.indmask);
jf_tmp = Gocrf_subindcalc(G.arg,real(asin(abs(m))),G.arg.indmask);
G = subsasgn_trick(G, 'arg', jf_tmp); % jf kludge for now

% then stick them into the Gmri objects
G = Gocrf_new_B_Ct(G,G.arg.a0B,G.arg.a0Ct,G.arg.b0B,G.arg.b0Ct,G.arg.a0,G.arg.b0);

if ~isempty(G.arg.RFbasis)
  % evaluate perturbations for basis functions
  for ll = 1:G.arg.Nbasis
    G.arg.Gat(:,ll) = G.arg.b0G * conj(G.arg.RFbasis(:,ll));
    G.arg.Gbt(:,ll) = G.arg.a0G * conj(G.arg.RFbasis(:,ll));
  end
end

%
% Gocrf_blochsim(): Perform a Bloch simulation
%
function [m,mz,a,b,a0Ct,b0Ct] = Gocrf_blochsim(B1,arg,pos,sens,we,simmask)

gambar = 42570000;               % gamma/2pi in Hz/T
gam = gambar*2*pi/10000;         % gamma in radians/g

if nargin == 6
  if ~isempty(simmask)
    % filter out un-sim'd spatial locs
    pos = pos(logical(simmask(:)),:);
    sens = sens(logical(simmask(:)),:);
    if ~isempty(we)
      we = we(logical(simmask));
    end
  end
end

Ns = size(pos,1);
Nt = size(arg.g,1);

statea = ones(Ns,1);
stateb = zeros(Ns,1);

% sum up RF over coils
B1 = reshape(B1,[size(arg.g,1) size(sens,2)]).';
bxy = sens * B1;

% sum up gradient over channels
bz = pos * arg.g';

% add off-resonance
if ~isempty(we) & sum(abs(we(:))) ~= 0
  bz = bz + repmat(we(:)/gam*2*pi,1,Nt);
end

if ~isempty(arg.fmod)
  bz = bz + arg.fmod/gam*2*pi;
end

tmp = zeros(2*Ns,1);

if nargout == 4
  % if we are only calling it with four output args,
  % assume we are doing the subspace sim and return all
  % states. Otherwise, we only return the final state.
  returnallstate = 1;
  a = zeros(Ns,Nt);
  b = zeros(Ns,Nt);
else; returnallstate = 0; end

% check if we need to do the running inner product to get spatial coeffs
if nargout == 6;
  inprod = 1;
  a0Ct = 0;b0Ct = 0;
  tvec = 0:arg.dt:(Nt-1)*arg.dt;
else; inprod = 0; end

ktmp = -[arg.k - repmat(arg.k(1,:),[size(arg.k,1) 1])];

for tt = 1:Nt

  phi = arg.dt*gam*(abs(bxy(:,tt)).^2+bz(:,tt).^2).^0.5;
  normfact = arg.dt*gam*(phi.^(-1));normfact(~isfinite(normfact)) = 0;
  nxy = normfact.*bxy(:,tt);nxy(~isfinite(nxy)) = 0;
  nz = normfact.*bz(:,tt);nz(~isfinite(nz)) = 0;
  cp = cos(phi/2);
  sp = sin(phi/2);
  alpha = cp+1i*nz.*sp;
  beta = 1i*conj(nxy).*sp;

  tmpa = alpha.*statea + beta.*stateb;
  tmpb = -conj(beta).*statea + conj(alpha).*stateb;

  statea = tmpa;stateb = tmpb;
  if any(~isfinite(statea)); keyboard; end;
  if any(~isfinite(stateb)); keyboard; end;
  if returnallstate
    a(:,tt) = statea;
    b(:,tt) = -conj(stateb);
  end

  if inprod
    % build drf transform vector
    Asub = exp(-1i*2*pi/2*pos*ktmp(tt,:).');
    if ~isempty(we) & sum(abs(we(:))) ~= 0
      Asub = Asub.*exp(1i*2*pi/2*we(:)*tvec(tt));
    end
    if ~isempty(arg.fmod)
      Asub = Asub.*exp(1i*2*pi/2*arg.fmod*ones(size(pos,1),1)*tvec(tt));
    end
    % get B matrix coeffs
    a0Ct = a0Ct + conj((conj(statea).*Asub)*arg.a0B(tt,:));
    b0Ct = b0Ct - conj((stateb.*Asub)*arg.b0B(tt,:));
  end

end

% return final alpha, beta if not returning the whole
% state progression
if ~returnallstate
  a = statea;
  b = -conj(stateb);
end

% calculate final magnetization state
mxy0 = arg.m0(1)+1i*arg.m0(2);
mz0 = arg.m0(3);
m = mz0*2*conj(statea).*stateb;
m = m + mxy0*conj(statea).^2;
m = m - conj(mxy0)*stateb.^2;
mz = mz0*(statea.*conj(statea) - stateb.*conj(stateb));
mz = mz + 2*real(mxy0*conj(statea).*-conj(stateb));

% embed results, if we did a masked sim
if nargin == 6
  if ~isempty(simmask)
    tmp = zeros(size(simmask(:)));tmp(logical(simmask)) = a;a = tmp;
    tmp = zeros(size(simmask(:)));tmp(logical(simmask)) = b;b = tmp;
    tmp = zeros(size(simmask(:)));tmp(logical(simmask)) = m;m = tmp;
    tmp = zeros(size(simmask(:)));tmp(logical(simmask)) = mz;mz = tmp;
    tmpa0Ct = zeros(length(simmask(:)),size(arg.a0B,2));
    tmpb0Ct = zeros(length(simmask(:)),size(arg.b0B,2));
    for ii = 1:size(arg.a0B,2)
      tmpa0Ct(logical(simmask),ii) = a0Ct(:,ii);
      tmpb0Ct(logical(simmask),ii) = b0Ct(:,ii);
    end
    a0Ct = tmpa0Ct;b0Ct = tmpb0Ct;
  end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [] = testLSTM(modelFiles, beamSize, stackSize, batchSize, outputFile,varargin)
% Test a trained LSTM model by generating translations.
% Arguments:
%   modelFiles: single or multiple models to decode. Multiple models are
%     separated by commas.
%   beamSize: number of hypotheses kept at each time step.
%   stackSize: number of translations retrieved.
%   batchSize: number of sentences decoded simultaneously. We only ensure
%     accuracy of batchSize = 1 for now.
%   outputFile: output translation file.
%   varargin: other optional arguments.
%
% Thang Luong @ 2015, <lmthang@stanford.edu>
% Hieu Pham @ 2015, <hyhieu@cs.stanford.edu>

  addpath(genpath(sprintf('%s/..', pwd)));

  %% Argument Parser
  p = inputParser;
  % required
  addRequired(p,'modelFiles',@ischar);
  addRequired(p,'beamSize',@isnumeric);
  addRequired(p,'stackSize',@isnumeric);
  addRequired(p,'batchSize',@isnumeric);
  addRequired(p,'outputFile',@ischar);

  % optional
  addOptional(p,'gpuDevice', 0, @isnumeric); % choose the gpuDevice to use: 0 -- no GPU
  addOptional(p,'align', 0, @isnumeric); % 1 -- output aignment from attention model
  addOptional(p,'assert', 0, @isnumeric); % 1 -- assert
  addOptional(p,'debug', 0, @isnumeric); % 1 -- debug
  addOptional(p,'minLenRatio', 0.5, @isnumeric); % decodeLen >= minLenRatio * srcMaxLen
  addOptional(p,'maxLenRatio', 1.5, @isnumeric); % decodeLen <= maxLenRatio * srcMaxLen
  addOptional(p,'testPrefix', '', @ischar); % to specify a different file for decoding
  addOptional(p,'hasTgt', 1, @isnumeric); % 0 -- no ref translations (groundtruth)
  addOptional(p,'continueId', 0, @isnumeric); % > 0: start decoding from this continueId (base 1) sent and append the results
  addOptional(p,'unkDiscount', 0, @isnumeric); % > 0: to prevent long unk sequences

  % force decoding: always feed the correct words (groundtruth)
  addOptional(p,'forceDecoder', 0, @isnumeric);

  p.KeepUnmatched = true;
  parse(p,modelFiles,beamSize,stackSize,batchSize,outputFile,varargin{:})
  decodeParams = p.Results;
  if decodeParams.batchSize==-1 % decode sents one by one
    decodeParams.batchSize = 1;
  end

  % GPU settings
  decodeParams.isGPU = 0;
  if decodeParams.gpuDevice
    n = gpuDeviceCount;
    if n>0 % GPU exists
      fprintf(2, '# %d GPUs exist. So, we will use GPUs.\n', n);
      decodeParams.isGPU = 1;
      gpuDevice(decodeParams.gpuDevice)
      decodeParams.dataType = 'single';
    else
      decodeParams.dataType = 'double';
    end
  else
    decodeParams.dataType = 'double';
  end
  printParams(2, decodeParams);

  %% load multiple models
  tokens = strsplit(decodeParams.modelFiles, ',');
  numModels = length(tokens);
  models = cell(numModels, 1);
  for mm=1:numModels
    modelFile = tokens{mm};
    [savedData] = load(modelFile);
    models{mm} = savedData.model;
    models{mm}
    models{mm}.params = savedData.params;

    % backward compatible
    if (isfield(models{mm}.params,'softmaxFeedInput')) && (~isfield(models{mm}.params,'feedInput'))
      models{mm}.params.feedInput = models{mm}.params.softmaxFeedInput;
    end
    [models{mm}.params] = backwardCompatible(models{mm}.params, {'charOpt', 'charSrcRep', 'charTgtGen', 'normLocalAttn'}, 0);

    % convert absolute paths to local paths
    fieldNames = fields(models{mm}.params);
    for ii=1:length(fieldNames)
     field = fieldNames{ii};
     if ischar(models{mm}.params.(field))
       if strfind(models{mm}.params.(field), '/afs/ir/users/l/m/lmthang') ==1
         models{mm}.params.(field) = strrep(models{mm}.params.(field), '/afs/ir/users/l/m/lmthang', '~');
       end
       if strfind(models{mm}.params.(field), '/afs/cs.stanford.edu/u/lmthang') ==1
         models{mm}.params.(field) = strrep(models{mm}.params.(field), '/afs/cs.stanford.edu/u/lmthang', '~');
       end
       if strfind(models{mm}.params.(field), '/home/lmthang') ==1
         models{mm}.params.(field) = strrep(models{mm}.params.(field), '/home/lmthang', '~');
       end
     end
    end

    % load vocabs
    [models{mm}.params] = prepareVocabs(models{mm}.params);

    % make sure all models have the same vocab, and the number of layers
    if mm>1
      for ii=1:models{mm}.params.srcVocabSize
        assert(strcmp(models{mm}.params.srcVocab{ii}, models{1}.params.srcVocab{ii}), '! model %d, mismatch src word %d: %s vs. %s\n', mm, ii, models{mm}.params.srcVocab{ii}, models{1}.params.srcVocab{ii});
      end
      for ii=1:models{mm}.params.tgtVocabSize
        assert(strcmp(models{mm}.params.tgtVocab{ii}, models{1}.params.tgtVocab{ii}), '! model %d, mismatch tgt word %d: %s vs. %s\n', mm, ii, models{mm}.params.tgtVocab{ii}, models{1}.params.tgtVocab{ii});
      end
      models{mm}.params = rmfield(models{mm}.params, {'srcVocab', 'tgtVocab'});
    end

    % copy fields
    fieldNames = fields(decodeParams);
    for ii=1:length(fieldNames)
      field = fieldNames{ii};
      if strcmp(field, 'testPrefix')==1 && strcmp(decodeParams.(field), '')==1 % skip empty testPrefix
        continue;
      elseif strcmp(field, 'testPrefix')==1
        fprintf(2, '# Decode a different test file %s\n', decodeParams.(field));
      end
      models{mm}.params.(field) = decodeParams.(field);
    end
  end

  params = models{1}.params;

  % force decode
  if decodeParams.forceDecoder
    params.stackSize = 1;
    params.beamSize = 1;
  end

  if params.continueId > 0 % appending
    fileOpt = 'a';
  else
    fileOpt = 'w';
  end

  params.fid = fopen(params.outputFile, fileOpt);
  if params.charTgtGen
    params.charFid = fopen([params.outputFile '.char'], fileOpt);
  end

  params.logId = fopen([outputFile '.log'], fileOpt);
  % align
  if params.align
    params.alignId = fopen([params.outputFile '.align'], fileOpt);
  end
  printParams(2, params);

  % load test data
  [srcSents, tgtSents, numSents]  = loadBiData(params, params.testPrefix, params.srcVocab, params.tgtVocab, -1, params.hasTgt);

  %%%%%%%%%%%%
  %% decode %%
  %%%%%%%%%%%%
  numBatches = floor((numSents-1)/batchSize) + 1;

  fprintf(2, '# Decoding %d sents, %s\n', numSents, datestr(now));
  fprintf(params.logId, '# Decoding %d sents, %s\n', numSents, datestr(now));
  startTime = clock;
  for batchId = 1 : numBatches
    % prepare batch data
    startId = (batchId-1)*batchSize+1;
    endId = batchId*batchSize;
    if endId > numSents
      endId = numSents;
    end

    % continue training
    if params.continueId > startId
      continue;
    end

    [decodeData] = prepareData(srcSents(startId:endId), tgtSents(startId:endId), 1, params);
    decodeData.startId = startId;

    % call lstmDecoder
    [candidates, candScores, alignInfo, otherInfo] = lstmDecoder(models, decodeData, params);

    % print results
    printDecodeResults(decodeData, candidates, candScores, alignInfo, otherInfo, params, 1);
  end

  endTime = clock;
  timeElapsed = etime(endTime, startTime);
  fprintf(2, '# Complete decoding %d sents, time %.0fs, %s\n', numSents, timeElapsed, datestr(now));
  fprintf(params.logId, '# Complete decoding %d sents, time %.0fs, %s\n', numSents, timeElapsed, datestr(now));

  fclose(params.fid);
  if params.charTgtGen
    fclose(params.charFid);
  end
  fclose(params.logId);
end

%     % convert local paths to absolute paths
%     fieldNames = fields(models{mm}.params);
%     for ii=1:length(fieldNames)
%       field = fieldNames{ii};
%       if ischar(models{mm}.params.(field))
%         if strfind(models{mm}.params.(field), '~lmthang/') ==1
%           models{mm}.params.(field) = strrep(models{mm}.params.(field), '~lmthang/', '/afs/ir/users/l/m/lmthang/');
%         end
%         if strfind(models{mm}.params.(field), '~lmthang/') ==1
%           models{mm}.params.(field) = strrep(models{mm}.params.(field), '~lmthang/', '/afs/cs.stanford.edu/u/lmthang/');
%         end
%         if strfind(models{mm}.params.(field), '~lmthang/') ==1
%           models{mm}.params.(field) = strrep(models{mm}.params.(field), '~lmthang/', '/home/lmthang/');
%         end
%       end
%     end

author: lmthang


text:
function this = MrDimInfo_constructor(this, dimInfoVariants)
% Unit test for MrDimInfo Constructor evoking all 6 variants
%
%   Y = MrUnitTest()
%   run(Y, 'MrDimInfo_constructor');
%
% This is a method of class MrUnitTest.
%
% IN
%
% OUT
%
% EXAMPLE
%   MrDimInfo_constructor
%
%   See also MrUnitTest

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

% Unit test for MrDimInfo Constructor evoking all 6 variants via
% dimInfoVariants

% construct MrDimInfo object from sampling points - this is the actual
% solution for variant (1) and the expected solution for all other cases
% using make_dimInfo_reference
dimInfo = this.make_dimInfo_reference(0);

switch dimInfoVariants
    %% (1) explicit setting of sampling points for dimension(s)
    case '1'

        % define actual solution
        actSolution = dimInfo;
        % load expected solution
        classesPath = tapas_uniqc_get_path('classes');
        solutionFileName = fullfile(classesPath, '@MrUnitTest' , 'dimInfo.mat');
        expSolution = load(solutionFileName);
        expSolution = expSolution.dimInfo;

        %% (2) nSamples + ranges
    case '2'
        % define expected solution
        expSolution = dimInfo;
        % define actual solution
        % Initialize MrDimInfo via nSamples + ranges (variant (2))
        actSolution = MrDimInfo(...
            'dimLabels', expSolution.dimLabels, ...
            'units', expSolution.units, ...
            'nSamples', expSolution.nSamples, ...
            'ranges', expSolution.ranges);
        % overwrite samplingPoints of coil due to non-uniform sampling
        actSolution.coil.samplingPoints = expSolution.coil.samplingPoints;
        % overwrite samplingWidths since they are retained in case of later
        % changes, but not set because of the construction of the reference
        % MrDimInfo object (= expected Solution)
        actSolution.coil.samplingWidths = nan;

        %% (3) nSamples + resolutions + samplingPoint + arrayIndex
    case '3'
        % define expected solution
        expSolution = dimInfo;
        % define actual solution
        arrayIndex = [81 54 14 321 1];

        actSolution = MrDimInfo(...
            'dimLabels', expSolution.dimLabels, ...
            'units', expSolution.units, ...
            'nSamples', expSolution.nSamples, ...
            'resolutions', expSolution.resolutions, ...
            'samplingPoint', expSolution.index2sample(arrayIndex), ....
            'arrayIndex', arrayIndex);

        % overwrite samplingPoints of coil due to non-uniform sampling
        actSolution.coil.samplingPoints = expSolution.coil.samplingPoints;

        % overwrite samplingWidths since they are retained in case of later
        % changes, but not set because of the construction of the reference
        % MrDimInfo object (= expected Solution)
        actSolution.coil.samplingWidths = nan;

        %% (4) nSamples + resolutions + firstSamplingPoint
    case '4'
        % define expected solution
        expSolution = dimInfo;
        % define actual solution
        arrayIndex = [1 1 1 1 1];

        actSolution = MrDimInfo(...
            'dimLabels', expSolution.dimLabels, ...
            'units', expSolution.units, ...
            'nSamples', expSolution.nSamples, ...
            'resolutions', expSolution.resolutions, ...
            'firstSamplingPoint', expSolution.index2sample(arrayIndex));

        % overwrite samplingPoints of coil due to non-uniform sampling
        actSolution.coil.samplingPoints = expSolution.coil.samplingPoints;

        % overwrite samplingWidths since they are retained in case of later
        % changes, but not set because of the construction of the reference
        % MrDimInfo object (= expected Solution)
        actSolution.coil.samplingWidths = nan;

        %% (5) nSamples + resolutions + lastSamplingPoint
    case '5'
        % define expected solution
        expSolution = dimInfo;
        % define actual solution
        arrayIndex = expSolution.nSamples;

        actSolution = MrDimInfo(...
            'dimLabels', expSolution.dimLabels, ...
            'units', expSolution.units, ...
            'nSamples', expSolution.nSamples, ...
            'resolutions', expSolution.resolutions, ...
            'lastSamplingPoint', expSolution.index2sample(arrayIndex));

        % overwrite samplingPoints of coil due to non-uniform sampling
        actSolution.coil.samplingPoints = expSolution.coil.samplingPoints;

        % overwrite samplingWidths since they are retained in case of later
        % changes, but not set because of the construction of the reference
        % MrDimInfo object (= expected Solution)
        actSolution.coil.samplingWidths = nan;

    case '6'
        %% (6) nSamples only
        % get nSamples from reference object
        nSamples = dimInfo.nSamples;
        nDims = numel(nSamples);

        % make actual solution
        actSolution = MrDimInfo('nSamples', nSamples);

        % make expected solution from ranges and sampling points

        for n = 1:nDims
            if n < 4
                % x y z
                samplingPoints{n} = -(nSamples(n)-1)/2:(nSamples(n)-1)/2;
            else
                samplingPoints{n} = 1:(nSamples(n));
            end
        end
        expSolution = MrDimInfo('samplingPoints', samplingPoints);

    case '7'
        %% (7) resolutions only
        % get resolutions from reference object
        resolutions = dimInfo.resolutions;
        % use only valid (i.e. finite)
        validResolutions = isfinite(resolutions);
        nDims = sum(validResolutions);

        % make actual solution
        actSolution = MrDimInfo('resolutions', resolutions(validResolutions));

        % make expected solution
        for n = 1:nDims
            if n < 4
                % x y z
                samplingPoints{n} = -(resolutions(n))/2:resolutions(n):(resolutions(n))/2;
            else
                samplingPoints{n} = resolutions(n):resolutions(n):2*resolutions(n);
            end
        end

        expSolution = MrDimInfo('samplingPoints', samplingPoints);

    case '8'
        %% (8) ranges only

        % get ranges from reference object
        ranges = dimInfo.ranges;

        % make actual solution
        actSolution = MrDimInfo('ranges', ranges);

        % make expected solution from sampling points
        samplingPoints = mat2cell(ranges', [1 1 1 1 1], 2)';
        expSolution = MrDimInfo('samplingPoints', samplingPoints);

    case '9'
        %% (9) dimLabels only

        % get dimLabels from reference object
        dimLabels = dimInfo.dimLabels;
        nDims = numel(dimLabels);

        % make actual solution
        actSolution = MrDimInfo('dimLabels', dimLabels);

        % make expected solution with standard sampling points, i.e [-0.5
        % 0.5] for dims x,y,z and [1 2] for all other dims

        expSolution.nDims = nDims;
        expSolution.nSamples = 2*ones(1, nDims);
        expSolution.resolutions = ones(1, nDims);
        expSolution.ranges = [-0.5 -0.5 -0.5 1 1; 0.5 0.5 0.5 2 2];
        for n = 1:nDims
            expSolution.dimLabels{n} = dimInfo.get_default_dim_labels(n);
            expSolution.units{n} = dimInfo.get_default_dim_units(n);
        end
        expSolution.samplingPoints = {[-0.5 0.5] [-0.5 0.5] [-0.5 0.5] [1 2] [1 2]};
        expSolution.samplingWidths = expSolution.resolutions;

    case '10'
        %% (10) units only

        % get units from reference object
        units = dimInfo.units;
        nDims = numel(units);

        % make actual solution
        actSolution = MrDimInfo('units', units);

        % make expected solution with standard sampling points, i.e [-0.5
        % 0.5] for dims x,y,z and [1 2] for all other dims
        expSolution.nDims = nDims;
        expSolution.nSamples = 2*ones(1, nDims);
        expSolution.resolutions = ones(1, nDims);
        expSolution.ranges = [-0.5 -0.5 -0.5 1 1; 0.5 0.5 0.5 2 2];
        for n = 1:nDims
            expSolution.dimLabels{n} = dimInfo.get_default_dim_labels(n);
            expSolution.units{n} = dimInfo.get_default_dim_units(n);
        end
        expSolution.samplingPoints = {[-0.5 0.5] [-0.5 0.5] [-0.5 0.5] [1 2] [1 2]};
        expSolution.samplingWidths = expSolution.resolutions;

    case '11'
        %% (11) samplingWidths only
        % note that samplingWidths does not impact resolutions and, thus,
        % sampling points and can therefore also be nan

        % get samplingWidths from reference object
        samplingWidths = dimInfo.samplingWidths;
        nDims = numel(samplingWidths);

        % make actual solution
        actSolution = MrDimInfo('samplingWidths', samplingWidths);

        % make expected solution with standard sampling points, i.e [-0.5
        % 0.5] for dims x,y,z and [1 2] for all other dims
        expSolution.nDims = nDims;
        expSolution.nSamples = 2*ones(1, nDims);
        expSolution.resolutions = ones(1, nDims);
        expSolution.ranges = [-0.5 -0.5 -0.5 1 1; 0.5 0.5 0.5 2 2];
        for n = 1:nDims
            expSolution.dimLabels{n} = dimInfo.get_default_dim_labels(n);
            expSolution.units{n} = dimInfo.get_default_dim_units(n);
        end
        expSolution.samplingPoints = {[-0.5 0.5] [-0.5 0.5] [-0.5 0.5] [1 2] [1 2]};
        expSolution.samplingWidths = samplingWidths;

end

% verify whether expected and actual solution are identical
% Note: convert to struct, since the PublicPropertyComparator (to allow
% nans to be treated as equal) does not compare properties of objects that
% overload subsref

warning('off', 'MATLAB:structOnObject');
this.verifyEqual(struct(actSolution), struct(expSolution), 'absTol', 10e-7);
warning('on', 'MATLAB:structOnObject');
end

author: translationalneuromodeling

repo: tapas


text:
b = ones(size(A,1),1);
tol = 1e-6; maxit = 100;
% tic; L1 = chol(A,'lower'); toc;
% fprintf('\n Incomplete chol decomposition');
% tic;
% L1 = ichol(A);
% toc;
% fprintf('\n ichol as Preconditioner');
% tic;
% [x1,fl1,rr1,it1,rv1] = pcg(A,b,tol,maxit,@(x)icholpre(x,A,L1,L1'));
% toc;
% fprintf('#dof: %8.0u,  iter: %2.0u\n',size(A,1), it1)
% semilogy(0:it1,rv1./norm(b),'r.');
%
fprintf('\n Approximate chol decomposition');
tic;
[L2,p,Ac] = achol(A);
toc;
fprintf('\n Achol as Preconditioner');
tic;
Ap = A(p,p);
[x2,fl2,rr2,it2,rv2] = pcg(A,b,tol,maxit,@(r)acholpre(r,Ap,L2,L2',p,Ac));
toc;
fprintf('#dof: %8.0u,  iter: %2.0u\n',size(A,1), it2)
semilogy(0:it2,rv2./norm(b),'b.');
hold on

repo: Aero-matlab


text:
function dat=func_ar(dat, order, varargin)
%Calculating the auto regression parameter / aar
%In dat      -     input the data structure of OpenBMI segementation original data
%   order    -     order of AR setting
%   varadgin -     model selection for obtatining AR parameter
%                  deafualt is 'aryule';
%                  model : 'arburg', 'arcov', 'armcov'

%out dat     -     data structure of otanined ar parameter in OpenBMI sturcture

% Example code func_AR(dat, 7, {'method','arburg'})
% Example code

opt=opt_cellToStruct(varargin{:});
opt=struct('method',opt.method);

if isempty(dat)
    warning('[OpenBMI] Warning! data is empty.');
end

if isempty(order)
    warning('[OpenBMI] Order is not exist.');
end

if isempty(opt.method) %method selection
   opt.method='aryule';
end

[T, nEvents , nChans]= size(dat.x);

temp_ar= [];
for i= 1:nChans*nEvents,
  ar= feval(opt.method, dat.x(:,i), order);
  temp_ar(:,i)= ar(2:end)';
end

dat.x= temp_ar;

author: PatternRecognition

repo: OpenBMI


text:
function debug(msg)
disp(msg);

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function [PC,CC,DC,IC] = trim_with_spline(PA,CA,PB,CB,tol)
  % TRIM_WITH_SPLINE Trim a given spline (PA,CA) with a "solid" spline (PB,CB)
  %
  % [PC,CC,DC,IC] = trim_with_spline(PA,CA,PB,CB,tol)
  %
  % Inputs:
  %   PA  #PA by dim list of control point locations
  %   CA  #CA by 4 list of indices into PA of cubic Bezier curves
  %   PB  #PB by dim list of control point locations
  %   CB  #CB by 4 list of indices into PB of cubic Bezier curves
  %   tol  tolerance for intersection {1e-7}
  % Outputs:
  %   PC  #PC by dim list of control point locations
  %   CC  #CC by 4 list of indices into PC of cubic Bezier curves
  %   DC  #CC list of flags whether inside B
  %   IC  #CC list of indices into CA
  %

  function [PA,CA,DA,IA] = trim_with_spline_helper(PA,CA,PB,CB)

    [A1,A2] = box_each_element(PA,CA);
    [B1,B2] = box_each_element(PB,CB);
    I = box_intersect(A1,A2,B1,B2);
    T = cell(size(CA,1),1);
    % consider each intersection, gather splits
    for ii = 1:size(I,1)
      ca = I(ii,1);
      cb = I(ii,2);
      Tii = cubic_cubic_intersect(PA(CA(ca,:),:),PB(CB(cb,:),:),tol);
      if ~isempty(Tii)
        T{ca} = [T{ca};Tii(:,1)];
      end
    end
    % conduct all splits
    IA = (1:size(CA,1))';
    for ca = 1:size(CA,1)
      if ~isempty(T{ca})
        Tca = sort(T{ca});
        [PAa,CAa] = cubic_subdivide(PA(CA(ca,:),:),Tca);
        CAa = reshape([CA(ca,1) CAa(2:end-1)+size(PA,1)-2 CA(ca,end)],size(CAa));
        IA([ca size(CA,1)+(1:size(CAa,1)-1)],:) = IA(ca);
        CA([ca size(CA,1)+(1:size(CAa,1)-1)],:) = CAa;
        PA = [PA;PAa(3:end,:)];
      end
    end

    %%fprintf('trim_with_spline...\n');
    %% loop order is imporant
    %IA = (1:size(CA,1))';
    %for cb = 1:size(CB,1)
    %  %fprintf('  %04d/%04d:\n',cb,size(CB,1));
    %  for ca = 1:size(CA,1)
    %    %progressbar(ca,size(CA,1));
    %    T = cubic_cubic_intersect(PA(CA(ca,:),:),PB(CB(cb,:),:),tol);
    %    if ~isempty(T)
    %      [PAa,CAa] = cubic_subdivide(PA(CA(ca,:),:),T(:,1));
    %      CAa = reshape([CA(ca,1) CAa(2:end-1)+size(PA,1)-2 CA(ca,end)],size(CAa));
    %      IA([ca size(CA,1)+(1:size(CAa,1)-1)],:) = IA(ca);
    %      CA([ca size(CA,1)+(1:size(CAa,1)-1)],:) = CAa;
    %      PA = [PA;PAa(3:end,:)];
    %    end
    %  end
    %end

    % midpoints of each cubic
    %fprintf('midpoints...\n');
    M = zeros(size(CA,1),2);
    for ca = 1:size(CA,1)
      M(ca,:) = cubic_eval(PA(CA(ca,:),:),0.5);
    end
    %fprintf('spline_winding_number...\n');
    DA = abs(spline_winding_number(PB,CB,M))>0.5;
  end

  m = size(CA,1);
  embed = @(P,C) reshape(P(C,:),size(C,1),8);
  CB8 = embed(PB,CB);
  % In either direction
  RB8 = embed(PB,fliplr(CB));
  CA8 = embed(PA,CA);
  % find perfect matches (only considering that every control point is exactly
  % the same; not finding partial co-incidence )
  [idx,dist] = rangesearch([CB8;RB8],CA8,0);
  perfect = cellfun(@(i) ~isempty(i),idx);
  JN = find(~perfect);
  [PA,CN,DN,IN] = trim_with_spline_helper(PA,CA(JN,:),PB,CB);
  CA = [CN;CA(perfect,:)];
  JP = find(perfect);
  IA = [JN(IN);JP];
  DA = [DN;true(numel(JP),1)];
  assert(size(CA,1) == numel(IA));
  assert(max(IA) <= m);
  assert(size(CA,1) == numel(DA));

  % rename
  PC = PA;
  CC = CA;
  DC = DA;
  IC = IA;
end

author: alecjacobson

repo: gptoolbox