text:
classdef(SharedTestFixtures = {DownloadGPT2Fixture}) tGPT2Tokenizer < matlab.unittest.TestCase
    % tGPT2Tokenizer   Tests for the text encoder class

    % Copyright 2020 The MathWorks, Inc.

    properties
        % Setup a default encoder for testing
        Encoder
    end

    properties(TestParameter)
        % UnicodeTextAndEncodedText
        UnicodeTextAndEncodedText = iGetUnicodeTextAndEncodedText()

        % Strings to be tokenized by the regex and expected results
        StringToTokenize = iStringToTokenize()

        % Test strings for encoding
        StringToEncode = iStringToEncode()
    end

    methods(TestClassSetup)
        function setup(test)
            test.setupEncoder();
        end
    end

    methods(Test)
        function encodeGivesCorrectResults(test, UnicodeTextAndEncodedText)
            % Unpack the test parameter
            inputText = UnicodeTextAndEncodedText.UnicodeText;
            expectedOutputText = UnicodeTextAndEncodedText.EncodedText;

            % Create the encoder
            enc = gpt2.tokenizer.GPT2Tokenizer(iModelName(), iModelPath());

            % Get the output text
            actualOutputText = enc.encode(inputText);

            % Verify the output is correct
            test.verifyThat( ...
                actualOutputText, iIsEqualTo(expectedOutputText) );
        end

        function decodeGivesCorrectResults(test, UnicodeTextAndEncodedText)
            % Unpack the test parameter
            inputText = UnicodeTextAndEncodedText.EncodedText;
            expectedOutputText = UnicodeTextAndEncodedText.UnicodeText;

            % Create the encoder
            enc = gpt2.tokenizer.GPT2Tokenizer(iModelName(), iModelPath());

            % Get the output text
            actualOutputText = enc.decode(inputText);

            % Verify the output is correct
            test.verifyThat( ...
                actualOutputText, iMatches(expectedOutputText) );
        end

        function testTokenization(test, StringToTokenize)
            % Test the regex tokenization used by the Encoder.
            enc = test.Encoder;
            toks = iTokenize(StringToTokenize.string,enc);
            test.verifyEqual(toks,StringToTokenize.exp);
        end

        function byteEncoderSize(test)
            % Verify the expected size of the byte encoder.
            enc = test.Encoder;
            test.verifyEqual(numel(unique(enc.ByteEncoder)),256);
        end

        function independenceOfInputClass(test, StringToEncode)
            % Verify that it doesn't matter if input text is string or char
            enc = test.Encoder;
            str = convertCharsToStrings(StringToEncode);
            ch = convertStringsToChars(StringToEncode);
            test.verifyEqual(enc.encode(str),enc.encode(ch));
        end

        function commentNotInBPE(test)
            % The vocab.bpe includes a comment that we must manually strip
            % out. Ensure we do that.
            vocabFile = gpt2.internal.getSupportFilePath("gpt2_vocab.bpe");
            fid = fopen(vocabFile,'r','n','utf-8');
            s = textscan(fid,'%s', 'Delimiter', '\n');
            fclose(fid);
            s = s{1};
            % first line of vocab.bpe is the comment
            firstLine = s(1);
            % split as in Encoder and check we don't accidentally add these
            % values to BPERanks.
            splitfirstLine = split(convertCharsToStrings(firstLine));
            enc = test.Encoder;
            actBPE = enc.BPERanks;
            % splitfirstLine is a string, ensure actBPE is
            test.assertClass(actBPE,'string');
            actBPEMatchesComment = any(actBPE==splitfirstLine(1),'all') || any(actBPE==splitfirstLine(2),'all');
            test.verifyFalse(actBPEMatchesComment, "A match was found between the encoder's BPERanks and the comment in vocab.bpe");
        end

        function decodeInvertsEncode(test,StringToEncode)
            % Ensure decoding inverts encoding
            enc = test.Encoder;
            encoded = enc.encode(StringToEncode);
            decoded = enc.decode(encoded);
            test.verifyMatches(decoded,StringToEncode);
        end
    end

    methods(Access=private)
        function setupEncoder(test)
            % Setup an encoder to be used by multiple tests.
            % Note this is a handle.
            test.Encoder = gpt2.tokenizer.GPT2Tokenizer(iModelName(),iModelPath());
        end
    end
end

function modelName = iModelName()
modelName = 'gpt2-355M';
end

function modelPath = iModelPath()
modelPath = fullfile(getRepoRoot());
end

function constraint = iIsEqualTo(varargin)
constraint = matlab.unittest.constraints.IsEqualTo(varargin{:});
end

function constraint = iMatches(varargin)
constraint = matlab.unittest.constraints.Matches(varargin{:});
end

function parameter = iGetUnicodeTextAndEncodedText()

parameter = struct;

parameter.NormalSentenceFragment1 = struct( ...
    'UnicodeText', 'In this tutorial we will see', ...
    'EncodedText', [818 428 11808 356 481 766] + 1 );

parameter.NonsenseFragment1 = struct( ...
    'UnicodeText', 'esgarghr', ...
    'EncodedText', [274 4563 456 81] + 1 );

% Weird edge case I found when I wasn't doing the Unicode decoding
% correctly. If you set "k=1" and input the sentence "In this tutorial ",
% into GPT-2 the first output is a special Unicode character known as a
% "Non-breaking space". After thinking about this, this weird behaviour is
% likely due to the fact that the tokenizer for GPT-2 tokenizes words so
% that they are preceded by spaces.
parameter.NonBreakingSpace = struct( ...
    'UnicodeText', native2unicode([194 160], 'UTF-8'), ...
    'EncodedText', 1849 + 1 );

end

function c = iStringToEncode()
% A test parameter of strings to be passed to the encoder.
c = {
    'In this tutorial we will see'
    'esgarghr'
    'foo bar baz'
    'A word and numbers 123 or 456word and decimal 7.89'
    'particles like don''t do that'
    'multiple white    space characters'
    };
end

function s = iStringToTokenize()
% A test parameter of strings to be tokenized using the regular expression
% attached to tokenizer.GPT2Tokenizer, and the expected output tokens.
s = struct( ...
    'oneWord', iStringToTokenizeStruct("foo","foo"), ...
    'twoWords', iStringToTokenizeStruct("foo bar",["foo"," bar"]), ...
    'possessive', iStringToTokenizeStruct("foo's",["foo","'s"]), ...
    'not', iStringToTokenizeStruct("foo't", ["foo","'t"]), ...
    'are', iStringToTokenizeStruct("foo're", ["foo","'re"]), ...
    'have', iStringToTokenizeStruct("foo've", ["foo","'ve"]), ...
    'm', iStringToTokenizeStruct("foo'm",["foo","'m"]), ...
    'will', iStringToTokenizeStruct("foo'll",["foo","'ll"]), ...
    'would', iStringToTokenizeStruct("foo'd",["foo","'d"]), ...
    'numbers', iStringToTokenizeStruct("foo 123bar",["foo"," 123","bar"]), ...
    'multiWS', iStringToTokenizeStruct("foo   bar",["foo", "  ", " bar"]));
end

function s = iStringToTokenizeStruct(str,expTokens)
% A struct for a case of the StringToTokenize test parameter.
s = struct('string',str,'exp',expTokens);
end

function toks = iTokenize(str,enc)
% Tokenize equivalently to tokenizer.GPT2Tokenizer
[toks,~] = regexp(str,enc.TokenizationExpression,'match','split');
end

author: matlab-deep-learning

repo: transformer-models


text:
%STUMPC Decision stump classifier
%
%   W = STUMPC(A,CRIT,N)
%   W = A*STUMPC([],CRIT,N)
%   W = A*STUMPC(CRIT,N)
%
% Computation of a decision tree classifier out of a dataset A using
% a binary splitting criterion CRIT:
%   INFCRIT  -  information gain
%   MAXCRIT  -  purity (default)
%   FISHCRIT -  Fisher criterion
% Just N (default N=1) nodes are computed.
%
% see also DATASETS, MAPPINGS, TREEC, TREE_MAP

% 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

% $Id: stumpc.m,v 1.2 2009/07/10 11:19:20 duin Exp $

function w = stumpc(varargin)

	mapname = 'Decision Stump';
  argin = shiftargin(varargin,{'char','integer'});
  argin = shiftargin(argin,'integer',2);
  argin = setdefaults(argin,[],'maxcrit',1);

  if mapping_task(argin,'definition')
    w = define_mapping(argin,'untrained',mapname);

  elseif mapping_task(argin,'training')			% Train a mapping.

    [a,crit,n] = deal(argin{:});
    islabtype(a,'crisp');
    isvaldset(a,1,2); % at least 1 object per class, 2 classes

    % First get some useful parameters:
    [m,k,c] = getsize(a);
    nlab = getnlab(a);
    tree = maketree(+a,nlab,c,crit,n);

    % Store the results:
    w = prmapping('tree_map','trained',{tree,1},getlablist(a),k,c);
    w = setname(w,mapname);
    w = setcost(w,a);

  end

	return

%MAKETREE General tree building algorithm
%
% 	tree = maketree(A,nlab,c,crit,stop)
%
% Constructs a binary decision tree using the criterion function
% specified in the string crit ('maxcrit', 'fishcrit' or 'infcrit'
% (default)) for a set of objects A. stop is a counter for the number of
% branches that are allowed below the present.
%
% Definition of the resulting tree:
%
% 	tree(n,1) - feature number to be used in node n
% 	tree(n,2) - threshold t to be used
% 	tree(n,3) - node to be processed if value <= t
% 	tree(n,4) - node to be processed if value > t
% 	tree(n,5:4+c) - aposteriori probabilities for all classes in
% 			node n
%
% If tree(n,3) == 0, stop, class in tree(n,1)
%
% This is a low-level routine called by treec.
%
% See also infstop, infcrit, maxcrit, fishcrit and mapt.

% Authors: Guido te Brake, TWI/SSOR, Delft University of Technology
%     R.P.W. Duin, TN/PH, Delft University of Technology
% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

function tree = maketree(a,nlab,c,crit,n)
		[m,k] = size(a);

	% Construct the tree:

	% When all objects have the same label, create an end-node:
	if all([nlab == nlab(1)])
		% Avoid giving 0-1 probabilities, but 'regularize' them a bit using
		% a 'uniform' Bayesian prior:
		p = ones(1,c)/(m+c); p(nlab(1)) = (m+1)/(m+c);
		tree = [nlab(1),0,0,0,p];
	else
		% now the tree is recursively constructed further:
		[f,j,t] = feval(crit,+a,nlab); % use desired split criterion

		p = sum(expandd(nlab),1);
		if length(p) < c, p = [p,zeros(1,c-length(p))]; end
		% When the stop criterion is not reached yet, we recursively split
		% further:

		if n >= 1
			% Make the left branch:
			J = find(a(:,j) <= t);
      tl = maketree(+a(J,:),nlab(J),c,crit,n-1);

      % Make the right branch:
			K = find(a(:,j) > t);
			tr = maketree(+a(K,:),nlab(K),c,crit,n-1);

      % Fix the node labelings before the branches can be 'glued'
			% together to a big tree:
			[t1,t2] = size(tl);
      tl = tl + [zeros(t1,2) tl(:,[3 4])>0 zeros(t1,c)];
			[t3,t4] = size(tr);
			tr = tr + (t1+1)*[zeros(t3,2) tr(:,[3 4])>0 zeros(t3,c)];

      % Make the complete tree: the split-node and the branches:
			tree= [[j,t,2,t1+2,(p+1)/(m+c)]; tl; tr];
		else
			% We reached the stop criterion, so make an end-node:
			[mt,cmax] = max(p);
			tree = [cmax,0,0,0,(p+1)/(m+c)];
		end

	end
	return

%MAXCRIT Maximum entropy criterion for best feature split.
%
% 	[f,j,t] = maxcrit(A,nlabels)
%
% Computes the value of the maximum purity f for all features over
% the data set A given its numeric labels. j is the optimum feature,
% t its threshold. This is a low level routine called for constructing
% decision trees.
%
% [1] L. Breiman, J.H. Friedman, R.A. Olshen, and C.J. Stone,
% Classification and regression trees, Wadsworth, California, 1984.

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

function [f,j,t] = maxcrit(a,nlab)
		[m,k] = size(a);
	c = max(nlab);
	% -variable T is an (2c)x k matrix containing:
	%      minimum feature values class 1
	%      maximum feature values class 1
	%      minimum feature values class 2
	%      maximum feature values class 2
	%            etc.
	% -variable R (same size) contains:
	%      fraction of objects which is < min. class 1.
	%      fraction of objects which is > max. class 1.
	%      fraction of objects which is < min. class 2.
	%      fraction of objects which is > max. class 2.
	%            etc.
	% These values are collected and computed in the next loop:
	T = zeros(2*c,k); R = zeros(2*c,k);
	for j = 1:c
		L = (nlab == j);
		if sum(L) == 0
			T([2*j-1:2*j],:) = zeros(2,k);
			R([2*j-1:2*j],:) = zeros(2,k);
		else
			T(2*j-1,:) = min(a(L,:),[],1);
			R(2*j-1,:) = sum(a < ones(m,1)*T(2*j-1,:),1);
			T(2*j,:) = max(a(L,:),[],1);
			R(2*j,:) = sum(a > ones(m,1)*T(2*j,:),1);
		end
	end
	% From R the purity index for all features is computed:
	G = R .* (m-R);
	% and the best feature is found:
	[gmax,tmax] = max(G,[],1);
	[f,j] = max(gmax);
	Tmax = tmax(j);
	if Tmax ~= 2*floor(Tmax/2)
		t = (T(Tmax,j) + max(a(find(a(:,j) < T(Tmax,j)),j)))/2;
	else
		t = (T(Tmax,j) + min(a(find(a(:,j) > T(Tmax,j)),j)))/2;
	end
	return

%INFCRIT The information gain and its the best feature split.
%
% 	[f,j,t] = infcrit(A,nlabels)
%
% Computes over all features the information gain f for its best
% threshold from the dataset A and its numeric labels. For f=1:
% perfect discrimination, f=0: complete mixture. j is the optimum
% feature, t its threshold. This is a lowlevel routine called for
% constructing decision trees.

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

function [g,j,t] = infcrit(a,nlab)
		[m,k] = size(a);
	c = max(nlab);
	mininfo = ones(k,2);
	% determine feature domains of interest
	[sn,ln] = min(a,[],1);
	[sx,lx] = max(a,[],1);
	JN = (nlab(:,ones(1,k)) == ones(m,1)*nlab(ln)') * realmax;
	JX = -(nlab(:,ones(1,k)) == ones(m,1)*nlab(lx)') * realmax;
	S = sort([sn; min(a+JN,[],1); max(a+JX,[],1); sx]);
	% S(2,:) to S(3,:) are interesting feature domains
	P = sort(a);
	Q = (P >= ones(m,1)*S(2,:)) & (P <= ones(m,1)*S(3,:));
	% these are the feature values in those domains
	for f=1:k,		% repeat for all features
		af = a(:,f);
		JQ = find(Q(:,f));
		SET = P(JQ,f)';
		if JQ(1) ~= 1
			SET = [P(JQ(1)-1,f), SET];
		end
		n = length(JQ);
		if JQ(n) ~= m
			SET = [SET, P(JQ(n)+1,f)];
		end
		n = length(SET) -1;
		T = (SET(1:n) + SET(2:n+1))/2; % all possible thresholds
		L = zeros(c,n); R = L;     % left and right node object counts per class
		for j = 1:c
			J = find(nlab==j); mj = length(J);
			if mj == 0
				L(j,:) = realmin*ones(1,n); R(j,:) = L(j,:);
			else
				L(j,:) = sum(repmat(af(J),1,n) <= repmat(T,mj,1)) + realmin;
				R(j,:) = sum(repmat(af(J),1,n) > repmat(T,mj,1)) + realmin;
			end
		end
		infomeas =  - (sum(L .* log10(L./(ones(c,1)*sum(L)))) ...
			       + sum(R .* log10(R./(ones(c,1)*sum(R))))) ...
		    ./ (log10(2)*(sum(L)+sum(R))); % criterion value for all thresholds
		[mininfo(f,1),j] = min(infomeas);     % finds the best
		mininfo(f,2) = T(j);     % and its threshold
	end
	g = 1-mininfo(:,1)';
	[finfo,j] = min(mininfo(:,1));		% best over all features
	t = mininfo(j,2);			% and its threshold
	return

%FISHCRIT Fisher's Criterion and its best feature split
%
% 	[f,j,t] = fishcrit(A,nlabels)
%
% Computes the value of the Fisher's criterion f for all features
% over the dataset A with given numeric labels. Two classes only. j
% is the optimum feature, t its threshold. This is a lowlevel
% routine called for constructing decision trees.

% Copyright R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

function [f,j,t] = fishcrit(a,nlab)
		[m,k] = size(a);
	c = max(nlab);
	if c > 2
		error('Not more than 2 classes allowed for Fisher Criterion')
	end
	% Get the mean and variances of both the classes:
	J1 = find(nlab==1);
	J2 = find(nlab==2);
	u = (mean(a(J1,:),1) - mean(a(J2,:),1)).^2;
	s = std(a(J1,:),0,1).^2 + std(a(J2,:),0,1).^2 + realmin;
	% The Fisher ratio becomes:
	f = u ./ s;
	% Find then the best feature:
	[ff,j] = max(f);
	% Given the feature, compute the threshold:
	m1 = mean(a(J1,j),1);
	m2 = mean(a(J2,j),1);
	w1 = m1 - m2; w2 = (m1*m1-m2*m2)/2;
	if abs(w1) < eps % the means are equal, so the Fisher
			 % criterion (should) become 0. Let us set the thresold
			 % halfway the domain
			 t = (max(a(J1,j),[],1) + minc(a(J2,j),[],1)) / 2;
	else
		t = w2/w1;
	end
	return

%INFSTOP Quinlan's Chi-square test for early stopping
%
% 	crt = infstop(A,nlabels,j,t)
%
% Computes the Chi-square test described by Quinlan [1] to be used
% in maketree for forward pruning (early stopping) using dataset A
% and its numeric labels. j is the feature used for splitting and t
% the threshold.
%
% [1] J.R. Quinlan, Simplifying Decision Trees,
% Int. J. Man - Machine Studies, vol. 27, 1987, pp. 221-234.
%
% See maketree, treec, classt, prune

% Guido te Brake, TWI/SSOR, TU Delft.
% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

function crt = infstop(a,nlab,j,t)
		[m,k] = size(a);
	c = max(nlab);
	aj = a(:,j);
	ELAB = expandd(nlab);
	L = sum(ELAB(aj <= t,:),1) + 0.001;
	R = sum(ELAB(aj > t,:),1) + 0.001;
	LL = (L+R) * sum(L) / m;
	RR = (L+R) * sum(R) / m;
	crt = sum(((L-LL).^2)./LL + ((R-RR).^2)./RR);
	return

%PRUNEP Pessimistic pruning of a decision tree
%
% 	tree = prunep(tree,a,nlab,num)
%
% Must be called by giving a tree and the training set a. num is the
% starting node, if omitted pruning starts at the root. Pessimistic
% pruning is defined by Quinlan.
%
% See also maketree, treec, mapt

% Guido te Brake, TWI/SSOR, TU Delft.
% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

function tree = prunep(tree,a,nlab,num)
		if nargin < 4, num = 1; end;
	[N,k] = size(a);
	c = size(tree,2)-4;
	if tree(num,3) == 0, return, end;
	w = prmapping('treec','trained',{tree,num},[1:c]',k,c);
	ttt=tree_map(prdataset(a,nlab),w);
	J = testc(ttt)*N;
	EA = J + nleaves(tree,num)./2;   % expected number of errors in tree
	P = sum(expandd(nlab,c),1);     % distribution of classes
					%disp([length(P) c])
					[pm,cm] = max(P);     % most frequent class
					E = N - pm;     % errors if substituted by leave
					SD = sqrt((EA * (N - EA))/N);
					if (E + 0.5) < (EA + SD)	     % clean tree while removing nodes
						[mt,kt] = size(tree);
						nodes = zeros(mt,1); nodes(num) = 1; n = 0;
						while sum(nodes) > n;	     % find all nodes to be removed
							n = sum(nodes);
							J = find(tree(:,3)>0 & nodes==1);
							nodes(tree(J,3)) = ones(length(J),1);
							nodes(tree(J,4)) = ones(length(J),1);
						end
						tree(num,:) = [cm 0 0 0 P/N];
						nodes(num) = 0; nc = cumsum(nodes);
						J = find(tree(:,3)>0);% update internal references
						tree(J,[3 4]) = tree(J,[3 4]) - reshape(nc(tree(J,[3 4])),length(J),2);
						tree = tree(~nodes,:);% remove obsolete nodes
					else
						K1 = find(a(:,tree(num,1)) <= tree(num,2));
						K2 = find(a(:,tree(num,1)) >  tree(num,2));

						tree = prunep(tree,a(K1,:),nlab(K1),tree(num,3));
						tree = prunep(tree,a(K2,:),nlab(K2),tree(num,4));
					end
					return

%PRUNET Prune tree by testset
%
% 	tree = prunet(tree,a)
%
% The test set a is used to prune a decision tree.

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

function tree = prunet(tree,a)
		[m,k] = size(a);
	[n,s] = size(tree);
	c = s-4;
	erre = zeros(1,n);
	deln = zeros(1,n);
	w = prmapping('treec','trained',{tree,1},[1:c]',k,c);
	[f,lab,nn] = tree_map(a,w);  % bug, this works only if a is dataset, labels ???
	[fmax,cmax] = max(tree(:,[5:4+c]),[],2);
	nngood = nn([1:n]'+(cmax-1)*n);
	errn = sum(nn,2) - nngood;% errors in each node
	sd = 1;
	while sd > 0
		erre = zeros(n,1);
		deln = zeros(1,n);
		endn = find(tree(:,3) == 0)';	% endnodes
		pendl = max(tree(:,3*ones(1,length(endn)))' == endn(ones(n,1),:)');
		pendr = max(tree(:,4*ones(1,length(endn)))' == endn(ones(n,1),:)');
		pend = find(pendl & pendr);		% parents of two endnodes
		erre(pend) = errn(tree(pend,3)) + errn(tree(pend,4));
		deln = pend(find(erre(pend) >= errn(pend))); % nodes to be leaved
		sd = length(deln);
		if sd > 0
			tree(tree(deln,3),:) = -1*ones(sd,s);
			tree(tree(deln,4),:) = -1*ones(sd,s);
			tree(deln,[1,2,3,4]) = [cmax(deln),zeros(sd,3)];
		end
	end
	return

%NLEAVES Computes the number of leaves in a decision tree
%
% 	number = nleaves(tree,num)
%
% This procedure counts the number of leaves in a (sub)tree of the
% tree by using num. If num is omitted, the root is taken (num = 1).
%
% This is a utility used by maketree.

% Guido te Brake, TWI/SSOR, TU Delft
% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

function number = nleaves(tree,num)
		if nargin < 2, num = 1; end
	if tree(num,3) == 0
		number = 1 ;
	else
		number = nleaves(tree,tree(num,3)) + nleaves(tree,tree(num,4));
	end
	return

author: marianux

repo: ecg-kit

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


text:
function reproducescript(reproducescript_dir, filename, highest_ft, cfg, outputcfg)

% This is a helper function to create a script that reproduces the analysis. It
% appends the configuration and the function call to a MATLAB script.

% Copyright (C) 2019, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if isfield(cfg, 'callinfo') && isfield(cfg.callinfo, 'usercfg')
  tmpcfg = removefields(cfg.callinfo.usercfg, ignorefields('reproducescript'));
else
  tmpcfg = removefields(cfg, ignorefields('reproducescript'));
end

tmpcfg = copyfields(cfg, tmpcfg, {'inputfile', 'outputfile'});
tmpcfg = printstruct('cfg', tmpcfg);

% We're going to explicitly mention any input files going into the pipeline
% at this point in a comment in the script file. We should mention any
% files that (1) have not been produced as an output in the same pipeline
% and that (2) have not been used as an input before. Do this by simply
% checking any .mat file references containing "input" in the stringified cfg against any .mat
% file references already in the script file.
re = ['(?<=' regexptranslate('escape', reproducescript_dir) '[/\\]{1})(\w+_input_\w+(\.mat){1})'];
if exist(filename, 'file')
  script = fileread(filename);
  existing_infiles = regexp(script, re, 'match');
else
  existing_infiles = {};
end
new_infiles = regexp(tmpcfg, re, 'match');
new_infiles = unique(setdiff(new_infiles, existing_infiles));

ft_info('writing script to file ''%s''\n', filename);
fid = fopen(filename, 'a+');
fprintf(fid, '%%%%\n\n');

for k = 1:numel(new_infiles)
  fprintf(fid, '%% a new input variable is entering the pipeline here: %s\n', new_infiles{k});
end
if numel(new_infiles) > 0
  fprintf(fid, '\n');
end

fprintf(fid, 'cfg = [];\n');
fprintf(fid, '%s\n', tmpcfg);

if outputcfg
  % this is for ft_definetrial, ft_artifact_zvalue, etc.
  fprintf(fid, 'cfg = %s(cfg);\n\n', highest_ft);
else
  fprintf(fid, '%s(cfg);\n\n', highest_ft);
end

fclose(fid);

author: spm


text:
if eps == 0
    DERIV_NOISE = 0;
end

%% Compute Derivative
Nt = size(x,1);
if any(strcmp(InputSignalType,{'sine2', 'chirp','prbs', 'sphs','mixed', 'noise','unforced'})==1) && DERIV_NOISE == 0%eps==0 % eps~=0
    % compute derivative using fourth order central difference
    % use TVRegDiff if more error
    dx = zeros(Nt-5,Nvar+1,Nic);
    for iIC = 1:Nic
    dx_tmp = zeros(Nt-5,3);
    for i=3:Nt-3
        for k=1:size(x,2)
            dx_tmp(i-2,k) = (1/(12*dt))*(-x(i+2,k,iIC)+8*x(i+1,k,iIC)-8*x(i-1,k,iIC)+x(i-2,k,iIC));
        end
    end
    dx(:,1:5,iIC) = dx_tmp;
    end

    % concatenate
    xnew = [x(3:end-3,:,:)];
    unew = [u(3:end-3,:,:)];

    disp(['DONE: computed derivative'])
else %excitation using Gaussian white noise

end

%% Reshape
% uensemble = [u;uv];
% uensemble = u;
X =[];
DX = [];
uensemble = [];
for iIC = 1:Nic
    X = [X; xnew(:,:,iIC)];
    DX = [DX; dx(:,:,iIC)];
    uensemble = [uensemble; unew(:,:,iIC)];
end
% xaug = [X repmat(uensemble(3:end-3,:),[Nic 1])];
xaug = [X uensemble];
DX(:,Nvar+1:size(xaug,2)) = repmat(0*DX(:,Nvar),[1 size(u,2)]);

M = size(xaug,1);

n = size(DX,2)-1;
%% Sparse regression
clear Theta Xi
Theta = poolData(xaug,n,polyorder,usesine);
Theta_norm = zeros(size(Theta,2),1);
for i = 1:size(Theta,2)
   Theta_norm(i) = norm(Theta(:,i));
   Theta(:,i) = Theta(:,i)./Theta_norm(i);
end
m = size(Theta,2);

if exist('lambda_vec') == 1
    Xi = sparsifyDynamicsIndependent(Theta,DX,lambda_vec,n);
else
    Xi = sparsifyDynamics(Theta,DX,lambda,n);
end

if n == 3
    str_vars = {'x','y','u'};
elseif n == 4
    str_vars = {'x','y','z','u'};
elseif n == 5
    str_vars = {'x1','x2','x3','x4','x5'};
elseif n == 6
    str_vars = {'x1','x2','x3','x4','x5','u'};
end

for i = 1:size(Theta,2)
   Xi(i,:) = Xi(i,:)./Theta_norm(i);
end

yout = poolDataLIST(str_vars,Xi,n,polyorder,usesine);

author: eurika-kaiser

repo: SINDY-MPC


text:
function rev = reversePolygon(poly)
%REVERSEPOLYGON Reverse a polygon, by iterating vertices from the end
%
%   POLY2 = reversePolygon(POLY)
%   POLY2 has same vertices as POLY, but in different order. The first
%   vertex of the polygon is still the same.
%
%   Example
%   reversePolygon
%
%   See also
%   polygons2d, reversePolyline
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-06-30,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

rev = poly([1 end:-1:2], :);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function view_stressmap(bvg) % autogenerated function wrapper
    % Script: view_stressmap.m
    % Author: S. Wiemer
    % last update: 19.05.2005, j.woessner@sed.ethz.ch
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    report_this_filefun(mfilename('fullpath'));

    SA = 1;
    SA2 = 3;
    % Matrix bvg contains:
    % bvg : [S1Trend S1Plunge S2Trend S2Plunge S3Trend S3Plunge Variance Radius b-value]
    % ste : [S1Plunge S1Trend+180 S2Plunge S2Trend+180 S3Plunge S3Trend+180 Variance];
    ste = [bvg(:,2) bvg(:,1)+180  bvg(:,4) bvg(:,3)+180 bvg(:,6) bvg(:,5)+180 bvg(:,7) ];
    sor = ste;
    ste2 = s;
    % sor : [S1Plunge S1Trend+270 S2Plunge S2Trend+180 S3Plunge S3Trend+180 Variance];
    sor(:,SA*2) = sor(:,SA*2)+90;

    % Create matrices
    normlap2=NaN(length(tmpgri(:,1)),1);

    re3=reshape(normlap2,length(yvect),length(xvect));
    s11 = re3;

    normlap2(ll)= bvg(:,8);
    rama=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,7);
    r=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,1);
    s11=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,4);
    s31=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,2);
    s12=reshape(normlap2,length(yvect),length(xvect));

    %old1 = re3;

    % Create figure
    figure_w_normalized_uicontrolunits('visible','off')
    l_normal =  ste(:,1) > 52 &   ste(:,5) < 35 ;
    l_notnormal = l_normal < 1;

    plq = quiver(newgri(l_notnormal,1),newgri(l_notnormal,2),-cos(sor(l_notnormal,SA*2)*pi/180),sin(sor(l_notnormal,SA*2)*pi/180),0.6,'.');
    set(plq,'LineWidth',0.5,'Color','k')
    px = get(plq,'Xdata');
    py = get(plq,'Ydata');

    close

    figure_w_normalized_uicontrolunits('visible','off')

    plq_n = quiver(newgri(l_normal,1),newgri(l_normal,2),-cos(sor(l_normal,SA2*2)*pi/180),sin(sor(l_normal,SA2*2)*pi/180),0.6,'.');
    set(plq_n,'LineWidth',0.5,'Color','r')

    drawnow
    px_n = get(plq_n,'Xdata');
    py_n = get(plq_n,'Ydata');
    close

    figure_w_normalized_uicontrolunits('Name','Faulting style map','pos',[100 100 860 600])
    watchon;
    %whitebg(gcf);
    set(gcf,'color','w');
    axes('pos',[0.12 0.12 0.8 0.8]);
    hold on
    n = 0;
    l0 = []; l1 = []; l2 = []; l3 = []; l4 = []; l5 = [];

    ste = sor(l_notnormal,:);
    for i = 1:3:length(px)-1
        n = n+1;j = jet;
        col = floor(ste(n,SA*2-1)/60*62)+1;
        if col > 64 ; col = 64; end
        pl = plot(px(i:i+1),py(i:i+1),'k','Linewidth',1.,'Markersize',1,'color',[ 0 0 0  ] );
        hold on

        dx = px(i)-px(i+1);
        dy = py(i) - py(i+1);
        pl2 = plot(px(i),py(i),'ko','Markersize',0.1,'Linewidth',0.5,'color',[0 0 0] );
        l0 = pl2;
        pl3 = plot([px(i) px(i)+dx],[py(i) py(i)+dy],'k','Linewidth',1.,'color',[0 0 0] );
        % Select faulting style according to Zoback(1992)
        if ste(n,1) < 40  && ste(n,3) > 45  && ste(n,5) < 20 ; set([pl pl3],'color',[0.2 0.8 0.2]); set(pl2,'color',[0.2 0.8 0.2]); l3 = pl; end
        if ste(n,1) < 20  && ste(n,3) > 45  && ste(n,5) < 40 ;  set([pl pl3],'color',[0.2 0.8 0.2]); set(pl2,'color',[0.2 0.8 0.2]);l3 = pl; end
        if ste(n,1) > 40  &&                 ste(n,1) < 52  && ste(n,5) < 20 ;  set([pl pl3],'color','m'); set(pl2,'color','m'); l2 = pl; end
        if ste(n,1) < 20  &&                 ste(n,5) > 40  &&  ste(n,5) <  52 ; set([pl pl3],'color','c'); set(pl2,'color','c');l4 = pl; end
        if ste(n,1) < 37  && ste(n,5) > 47  ; set([pl pl3],'color','b'); set(pl2,'color','b');l5 = pl;  end

    end
    %drawnow
    ste = sor(l_normal,:);
    n = 0;

    for i = 1:3:length(px_n)-1
        n = n+1;j = jet;
        col = floor(ste(n,SA*2-1)/60*62)+1;
        if col > 64 ; col = 64; end
        dx_n = px_n(i)-px_n(i+1);
        dy_n= py_n(i) - py_n(i+1);
        pl_n = plot(px_n(i:i+1),py_n(i:i+1),'k','Linewidth',1.,'Markersize',1,'color',[ 0 0 0  ] );
        hold on
        dx = px_n(i)- px_n(i+1);
        dy = py_n(i) - py_n(i+1);
        pl2_n = plot(px_n(i),py_n(i),'ko','Markersize',0.1,'Linewidth',0.5,'color',[0 0 0] );
        l0 = pl2;
        pl3_n = plot([px_n(i) px_n(i)+dx_n],[py_n(i) py_n(i)+dy_n],'k','Linewidth',1.,'color',[0 0 0] );

        if ste(n,1) > 52  &&                 ste(n,5) < 35 ;                 set([pl_n pl3_n],'color','r'); set(pl2_n,'color','r'); l1 = pl_n;
        end
    end

    if isempty(l1); pl2 = plot(px,py,'kx','Linewidth',1.,'color','r'); l1 = pl2; set(l1,'visible','off'); end
    if isempty(l2); pl2 = plot(px,py,'kx','Linewidth',1.,'color','m'); l2 = pl2; set(l2,'visible','off'); end
    if isempty(l3); pl2 = plot(px,py,'kx','Linewidth',1.,'color',[0.2 0.8 0.2] ); l3 = pl2; set(l3,'visible','off'); end
    if isempty(l4); pl2 = plot(px,py,'kx','Linewidth',1.,'color','c' ); l4 = pl2; set(l4,'visible','off'); end
    if isempty(l5); pl2 = plot(px,py,'kx','Linewidth',1.,'color','b' ); l5 = pl2; set(l5,'visible','off'); end
    if isempty(l0); l0 = plot(px,py,'kx','Linewidth',1.,'color',[0 0 0 ] );  set(l0,'visible','off'); end

    try
        legend([l1 l2 l3 l4 l5 l0],'NF','NS','SS','TS','TF','U');
    catch
        disp('Legend could not be drawn')
    end

    % Figure settings
    hold on
    axis('equal')
    update(mainmap())
    set(gca,'aspectratio',[0.827 1])
    axis([ s2_west s1_east s4_south s3_north])
    title([name ';  '   num2str(t0b) ' to ' num2str(teb) ],'FontSize',ZmapGlobal.Data.fontsz.s,...
        'Color','k','FontWeight','normal');
    mygca = gca;
    xlabel('Longitude ','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s)
    ylabel('Latitude ','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s)
    set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','normal',...
        'FontWeight','normal','LineWidth',1.,...
        'Box','on','TickDir','out','Ticklength',[0.01 0.01])

    watchoff;

    % View the variance map
    re3 = r;sha = 'in';
    view_varmap([],re3);

end

author: CelsoReyes


text:
function [ scaled_regions ] = map_regions_to_mult_scales(region_params, regions, image_size, scales )
% map_regions_to_mult_scales: given a set of regions that comes from an
% image of size image_size, it maps each of them to the appropriate image
% scale from set of scales given from the parameter scales.
%
% INPUTS:
% 1) region_params: (type struct) the region pooling parameters. Some of
% its fields are:
%   a) sz_conv_standard: (scalar value) the last convolution size
%   b) step_standard: (scalar value) is the stride in pixels of the output
%   of the last convolutional layer of the network where the regions are
%   going to be projected (for VGG16 is equal to 16).
% 2) regions: is a N x 8 array with the N input regions. Each region is
% represented by 8 values [xo0, yo0, xo1, yo1, xi0, yi0, xi1, yi1] that
% correspond to its outer rectangle [xo0, yo0, xo1, yo1] and its inner
% rectangle [xi0, yi0, xi1, yi1].
% 3) image_size: 2 x 1 or 1 x 2 array with the size of the original image
% 4) scales: NS x 1 or 1 x NS array with the image scales that are used. NS
% is the number of images. The i-th value is the size in pixels of the
% smallest dimension of the image in the i-th scale.
%
% OUTPUT:
% 1) scaled_regions: is a N x 9 array with the N output regions. Each region is
% represented by 9 values [scale_id, xo0, yo0, xo1, yo1, xi0, yi0, xi1, yi1]
% that correspond to its outer rectangle [xo0, yo0, xo1, yo1] and its inner
% rectangle [xi0, yi0, xi1, yi1]. scale_id corresponds to the scale id to
% which each region is mapped.
%
% This file is part of the code that implements the following paper:
% Title      : "LocNet: Improving Localization Accuracy for Object Detection"
% Authors    : Spyros Gidaris, Nikos Komodakis
% Institution: Universite Paris Est, Ecole des Ponts ParisTech
% ArXiv link : http://arxiv.org/abs/1511.07763
% code       : https://github.com/gidariss/LocNet
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2016 Spyros Gidaris
%
% Title     : "LocNet: Improving Localization Accuracy for Object Detection"
% ArXiv link: http://arxiv.org/abs/1511.07763
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

boxes_out = regions(:,1:4);
boxes_in  = regions(:,5:8);

min_img_sz = min(image_size(1:2));

if length(scales) > 1
    box_areas      = (boxes_out(:,3) - boxes_out(:, 1) + 1) .* (boxes_out(:,4) - boxes_out(:,2) + 1);
    expected_scale = region_params.sz_conv_standard * region_params.step_standard * min_img_sz ./ sqrt(box_areas);
    expected_scale = round(expected_scale(:));
    [~, best_scale_ids] = min(abs(bsxfun(@minus, scales, expected_scale(:))), [], 2);
else
    best_scale_ids = ones(size(boxes_out, 1), 1);
end

boxes_scales     = scales(best_scale_ids(:));
scaled_boxes_out = bsxfun(@times, (boxes_out - 1), (boxes_scales(:) - 1)) / (min_img_sz - 1) + 1;
scaled_boxes_in  = bsxfun(@times, (boxes_in  - 1), (boxes_scales(:) - 1)) / (min_img_sz - 1) + 1;
scaled_regions   = [best_scale_ids(:), scaled_boxes_out, scaled_boxes_in];
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
close all;
clear all;
clc;
rng('default');

png_export = true;
pdf_export = false;

mf = spx.graphics.Figures();

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

mf.new_figure('Recovery probability with S for different MMV recovery algorithms');

load ('bin/figure_2_success_with_s_comparison.mat');

hold all;
legends = cell(1, 4);

plot(Ss, success_with_s.ra_ormp, '-+');
legends{1} = 'RA-ORMP';
plot(Ss, success_with_s.ra_omp, '-o');
legends{2} = 'RA-OMP';
plot(Ss, success_with_s.somp, '-s');
legends{3} = 'SOMP';
plot(Ss, success_with_s.ra_thresholding, '-d');
legends{4} = 'RA-Thresholding';

grid on;
xlabel('Number of signals');
ylabel('Recovery Probability');
legend(legends, 'Location', 'southeast');
title('Comparison of recovery performance for MMV algorithms');

author: indigits

repo: sparse-plex


text:
function burnInfoText = printDVManeuversFMSToTextbox(hDvManInfoText, departBody, flybyBody, dVDepartVectNTW, ePreDepartOrbit, eTA, departBodyInfo, flyByDVVectNTW,  form, paddLen)
%printDVManeuversFMSToTextbox Summary of this function goes here
%   Detailed explanation goes here
        hRule = getHRule();

        ePeriod = computePeriod(ePreDepartOrbit(1), departBodyInfo.gm);

        burnInfoText{1} = ['Burn Information to Depart ', cap1stLetter(departBody)];
        burnInfoText{end+1} = hRule;
        burnInfoText{end+1} = [paddStr('Total Delta-V = ',paddLen), num2str(norm(dVDepartVectNTW), form), ' km/s'];
        burnInfoText{end+1} = [paddStr('Prograde Delta-V = ',paddLen), num2str(1000*dVDepartVectNTW(1), form), ' m/s'];
        burnInfoText{end+1} = [paddStr('Orbit Normal Delta-V = ',paddLen), num2str(1000*dVDepartVectNTW(2), form), ' m/s'];
        burnInfoText{end+1} = [paddStr('Radial Delta-V = ',paddLen), num2str(1000*dVDepartVectNTW(3), form), ' m/s'];
        burnInfoText{end+1} = '---------------------';
        burnInfoText{end+1} = [paddStr('Departure True Anomaly = ',paddLen), num2str(rad2deg(eTA), form), ' deg'];
        MA = AngleZero2Pi(computeMeanFromTrueAnom(eTA, ePreDepartOrbit(2)));
        eN = computeMeanMotion(ePreDepartOrbit(1), departBodyInfo.gm);
        secFromPeri = MA/eN;
        burnInfoText{end+1} = [paddStr('Departure Time Past Peri. = ',paddLen), num2str(secFromPeri, form), ' sec'];
        burnInfoText{end+1} = [paddStr('Departure Time Before Peri. = ',paddLen), num2str(ePeriod-secFromPeri, form), ' sec'];
        burnInfoText{end+1} = hRule;
        burnInfoText{end+1} = ['Burn Information to Depart ', cap1stLetter(flybyBody)];
        burnInfoText{end+1} = hRule;
        burnInfoText{end+1} = [paddStr('Total Delta-V = ',paddLen), num2str(norm(flyByDVVectNTW), form), ' km/s'];
        burnInfoText{end+1} = [paddStr('Prograde Delta-V = ',paddLen), num2str(1000*flyByDVVectNTW(1), form), ' m/s'];
        burnInfoText{end+1} = [paddStr('Orbit Normal Delta-V = ',paddLen), num2str(1000*flyByDVVectNTW(2), form), ' m/s'];
        burnInfoText{end+1} = [paddStr('Radial Delta-V = ',paddLen), num2str(1000*flyByDVVectNTW(3), form), ' m/s'];
        %flyby burn occurs at periapse per the algorithm - this is hardcoded below (next two lines)
        burnInfoText{end+1} = '---------------------';
        burnInfoText{end+1} = [paddStr('Departure True Anomaly = ',paddLen), num2str(rad2deg(0.0), form), ' deg'];
        burnInfoText{end+1} = [paddStr('Departure Time Past Peri. = ',paddLen), num2str(0.0, form), ' sec'];
        burnInfoText{end+1} = [paddStr('Departure Time Before Peri. = ',paddLen), num2str(0.0, form), ' sec'];
        set(hDvManInfoText,'String',burnInfoText);

        %first number is type of time in the second element: 0 is UT, 1 is
        %time past periapse;
        %Note: Code after this method is called may modify UserData
        %(specifically to change the type to strict UT time)
        set(hDvManInfoText, 'UserData', [1, secFromPeri, 1000*dVDepartVectNTW(1), 1000*dVDepartVectNTW(2), 1000*dVDepartVectNTW(3)]);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ IDX,timestamps ] = bz_INTtoIDX(INT,varargin)
%[IDX] = bz_INTtoIDX(INT) Converts state on/offsets to vector of indices
%
%INPUT
%   INT:    {nstates} cell array of [nintervals x 2] start and end times.
%                       (optional) can be TSObject intervalSet
%           -or-
%           structure with fields: INT.NAMEstate
%   (optional)
%   'statenames'    cell array of state names (position correponds to state number)
%   'sf'            desired sampling frequency of the output vector (default: 1s)
%   'length'        desired length of the index vector (default: max end time)
%
%OUTPUT
%   IDX:    [len x 1] vector of state indices, where states are identified by
%           integers starting from 1, 0 are unmarked.
%           -or-
%           structure with fields:  IDX.statenames
%           (see buzcode wiki)      IDX.states
%                                   IDX.timestamps
%   timestamps
%
%DLevenstein 2015
%Updated 2018 for buzcode
%% DEV
%INT = SleepState.ints;
%statenames = {'WAKE','','NREM','','REM'};
%%
p = inputParser;
addParameter(p,'statenames',[],@iscell)
addParameter(p,'sf',1)
addParameter(p,'length',Inf)
parse(p,varargin{:})
statenames = p.Results.statenames;
sf = p.Results.sf;
len = p.Results.length;
%% Deal with Input Types
%For Buzcode ints structure
if isstruct(INT)
    STRUCTIN = true;
    %Get the NAMEs of the states from the structure (INT.NAMEstate)
    fields = fieldnames(INT);
    endState = cellfun(@(X) contains(X,'state'),fields)';
    fieldstates = cellfun(@(X) char(extractBefore(X,'state')),fields(endState),'uniformoutput',false)';

    if ~isempty(statenames) %Check if there are any states the user didn't put in
        samestates = ismember(fieldstates,statenames);
        %Here: if you find any states in struct not in input, prompt user
        for ss = 1:length(samestates); if samestates(ss)==0
                statenum = inputdlg(['What number would you like state ''',fieldstates{ss},''' to be']);
                statenum = str2num(statenum{1});
                statenames{statenum} = fieldstates{ss};
        end; end
    else
        statenames = fieldstates;
    end

    %Convert to the cell array format needed for this function
    for ss = 1:length(statenames)
        if isempty(statenames{ss})
            continue;
        elseif  ~ismember(statenames{ss},fieldstates)
            INTtemp{ss} = [];
            continue;
        end
        INTtemp{ss} = INT.([statenames{ss},'state']);
    end
    INT = INTtemp;
else
    STRUCTIN = false;
end

%TSToolbox
if isa(INT,'intervalSet')
    INT = {[Start(INT,'s'), End(INT,'s')]};
end

%%

%Convert from seconds to dt = 1/sf
INT = cellfun(@(X) X*sf,INT,'UniformOutput',false);

%Getting the length of the output
if isinf(len)
    allints = cat(1,INT{:});
    len = ceil(max(allints(:,2)));
end

%%

IDX = zeros(len,1);

numstates = length(INT);
statenums = 1:numstates; %for possible later implementation of non 1:numstates nums
for ss = statenums
    if isempty(INT{ss}); continue; end
    stateints = round(INT{ss});

    stateints(stateints==0)=1;
    stateints(isinf(stateints))=len;

    numints = length(stateints(:,1));
    for ii = 1:numints
        IDX(stateints(ii,1):stateints(ii,2))=ss;
    end
end

switch numstates
    case 1
        IDX = logical(IDX);
    otherwise
end

timestamps = [1:length(IDX)]'./sf;

%% Convert to structure
if STRUCTIN
    IDXstruct.states = IDX;
    IDXstruct.timestamps = timestamps;
    IDXstruct.statenames = statenames;
    IDX = IDXstruct;
end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function []=sb_make_closure_loops()
%SB_MAKE_CLOSURE_LOOPS
%   Produce a matrix containign 1,-1 or 0 elements definining which
%   interferograms are used in a closed loop test.
%   Ekbal Hussain @ Uni Leeds
%   25/10/2013
%
%   ======================================================================
%   06/2014 AH: Modify for StaMPS
%   08/2014 EH: Fixed bug for determining higher order loops
%   08/2014 EH: Ignores interferograms not connected in a closed loop, bug fixes
%   01/2015 EH: User input required for individual ints (not in a loop) connecting two looped networks
%   01/2015 EH: Bug fixes
%   04/2015 EH: Checks if you are running small baselines and not single master. Speed optimisations
%   04/2015 EH: Removed user input setion for ints connecting two loops
%   04/2015 EH: Allows for up to 12 interferograms in a loop
%   ======================================================================

logit;
fprintf('Making closure loops...\n')

load psver
psname=['ps',num2str(psver)];
loopname=['phuw_loops',num2str(psver)];

ps=load(psname);

% Check you are running small baselines and not single master
small_baseline_flag = getparm('small_baseline_flag');
if strcmpi(small_baseline_flag,'n')
    error('Error: This method does not work for single master. Use small baselines instead.')
end

%Interferogram indices
x = 1:ps.n_ifg; x = x';
intfg = [ps.ifgday_ix,x];          %where: col1=date 1, col2=dat2, col3=int_no.

%Remove interferograms dropped in previous stamps steps
drop_ifg = getparm('drop_ifg_index',1);
if ~isempty(drop_ifg)
    for n = 1:size(drop_ifg,2)
       intfg = intfg(intfg(:,3)~=drop_ifg(n),:);
    end
end

% loading the single master bperp
small_baseline_flag = getparm('small_baseline_flag');
if strcmpi(small_baseline_flag,'y')
    temp = load('../psver');
    psname_ps = ['..' filesep 'ps' num2str(temp.psver) '.mat'];
    load(psname_ps,'bperp')
else
   load(psname,'bperp')
end

clear x n

%% Create interferogram triangulalation matrix

%Acquisition date (node) indexes
iz = unique(intfg(:,1));

%Interferograms used
ints_used = [];
closed_loops = [];
for m=1:size(iz,1)

    %set node ID
    node_id = iz(m);

    iy = intfg(intfg(:,1)==node_id,:);      %Interferograms that connect to this node [nID nID intID]
    a = unique(iy(:,1:2));                  %The nodes used in the combinations

    if length(a) >= 3                       %Can be used in one or more loops

        %Create all possible combinations using these nodes and output the subset that
        %only begins with the selected start node.
        combinations = combnk(a,3);
        usecom = combinations((combinations(:,1)==node_id),:);

        int_list1 = [usecom(:,1) usecom(:,2)];  %Positive interfeograms in the loop
        int_list2 = [usecom(:,2) usecom(:,3)];  %Positive interfeograms in the loop
        int_list3 = [usecom(:,1) usecom(:,3)];  %Negative interfeograms in the loop

        %Find the interferogram index connecting the nodes
        %First positive interferogram
        temp_column = nan(size(int_list1,1),1);
        for k=1:size(int_list1,1)
           [~,~,ix] = intersect(int_list1(k,:),intfg(:,1:2),'rows');
           if isempty(ix)
               temp_column(k,1)=NaN;
           else
               temp_column(k,1)=intfg(ix,3);
           end
        end
        int_list1 = [int_list1 temp_column];
        inan1 = ~isnan(int_list1(:,3));
        clear temp_column temp

        %Second positive interferogram
        temp_column = nan(size(int_list2,1),1);
        for k=1:size(int_list2,1)
           [~,~,ix] = intersect(int_list2(k,:),intfg(:,1:2),'rows');
           if isempty(ix)
               temp_column(k,1)=NaN;
           else
               temp_column(k,1)=intfg(ix,3);
           end
        end
        int_list2 = [int_list2 temp_column];
        inan2 = ~isnan(int_list2(:,3));
        clear temp_column temp

        %The negative interferogram
        temp_column = nan(size(int_list3,1),1);
        for k=1:size(int_list3,1)
           [~,~,ix] = intersect(int_list3(k,:),intfg(:,1:2),'rows');
           if isempty(ix)
               temp_column(k,1)=NaN;
           else
               temp_column(k,1)=intfg(ix,3);
           end
        end
        int_list3 = [int_list3 temp_column];
        inan3 = ~isnan(int_list3(:,3));
        clear temp_column temp

        %Remove any rows that contain NaNs
        igood = [inan1, inan2, inan3];
        for k=1:size(igood,1)
            if igood(k,:) == [1 1 1];
                ixgood(k,1) = 1;
            else
                ixgood(k,1) = NaN;
            end
        end
        int_list1 = int_list1(~isnan(ixgood),:);
        int_list2 = int_list2(~isnan(ixgood),:);
        int_list3 = int_list3(~isnan(ixgood),:);
        clear inan1 inan2 inan3 ixgood

        %Create a vector of the interferograms used
        ints_used = unique([ints_used;int_list1(:,3);int_list2(:,3);int_list3(:,3);]);

        %Add these to a loop closure matrix
        closed_loops = [closed_loops; int_list1(:,3) int_list2(:,3) -int_list3(:,3);];      %This array contains only int IDs
        clear int_list1 int_list2 int_list3
    end
end

clear combinations usecom iy a iz k igood ix node_id small_baseline_flag

%% Higher order loops for those not fitted by a triangle

%Find the ints that were not fitted by a triangle
A = [closed_loops(:,1:2),-closed_loops(:,3)];
missing_loops = intfg(~ismember(intfg(:,3),A),:);

%Ignore unconnected single interferograms
missing = [];
if ~isempty(missing_loops)
    for i=1:size(missing_loops,1)
        occurences_node1  = sum(sum(intfg(:,1:2)==missing_loops(i,1)),2);
        occurences_node2  = sum(sum(intfg(:,1:2)==missing_loops(i,2)),2);
        if occurences_node1 ~= 1 && occurences_node2 ~= 1
            missing = vertcat(missing,missing_loops(i,3));
        end
    end
end

remove_from_missing = [];
if ~isempty(missing)

    % loading the single master perpendicular baselines
    temp = load('../psver');
    psname_ps = ['..' filesep 'ps' num2str(temp.psver) '.mat'];
    sm_bp = load(psname_ps,'bperp');
    bp = sm_bp.bperp;

    Nodes = [(1:size(bp,1))', ps.day, bp];
    segments = [intfg(:,3), intfg(:,1:2)];

    % Allow for loops containing up to 8 interferograms
    higherO_loops = nan(size(missing,1),12);
    for m=1:size(missing,1)

        A = intfg(intfg(:,3)==missing(m,1),1:2);

        %Remove the shortest link between the two nodes
        Segments = segments;
        Segments(Segments(:,1) == missing(m,1),:) = [];

        %Calculate the shortest path that links the two nodes using
        %dijkstra's algorithm
        [~,path] = dijkstra(Nodes,Segments,A(:,1),A(:,2));

        %Remove individual ints (not in a looped network) connecting two looped networks
        if isnan(path)
            disp(['Interferogram ',num2str(missing(m,:)),' is not connected to a loop.']);

            %remove this interferogram from the missing list
            remove_from_missing = [remove_from_missing;m];
        else
            %close the loop
            path = [path, path(1,1)];

            %Find the interferogram index for this path
            for k=1:(size(path,2)-1)
                if path(1,k+1) > path(1,k);
                    temp_vector = [path(1,k),path(1,k+1)];
                    [~,~,ix] = intersect(temp_vector,intfg(:,1:2),'rows');
                    higherO_loops(m,k)=intfg(ix,3);
                else
                    temp_vector = [path(1,k+1),path(1,k)];
                    [~,~,ix] = intersect(temp_vector,intfg(:,1:2),'rows');
                    higherO_loops(m,k)=-intfg(ix,3);
                end
            end
        end
    end

    %Remove rows and colums containing just nans
    higherO_loops = higherO_loops(~all(isnan(higherO_loops),2),~all(isnan(higherO_loops),1));
    %replace nans with zero (makes the next bit easier)
    higherO_loops(isnan(higherO_loops)) = 0;

    %Remove duplicated loops
    x = sort(abs(higherO_loops),2);
    y =  unique(x,'rows');
    [~,ix] = ismember(y,x,'rows');

    higherO_loops = higherO_loops(ix,:);
    %replace any zeros with nans
    higherO_loops(higherO_loops==0) = nan;
end

clear A i occurences_node1 occurences_node2 missing_loops psname_ps temp sm_bp bp m path k temp_vector ix Nodes segments Segments x y

%% The final closed network

%Remove the single interferograms connecting two loops
missing(remove_from_missing,:) = [];

if ~isempty(missing)
    intfg_network = nan((size(closed_loops,1)+size(higherO_loops,1)),size(higherO_loops,2));

    %Create the new matrix here. Add the triangular loops
    intfg_network(1:size(closed_loops,1),1:size(closed_loops,2)) = closed_loops;
    %Add the higher order loops
    intfg_network((size(closed_loops,1)+1):end,:) = higherO_loops;
else
    intfg_network = closed_loops;
end

%Turn this into a full matrix including flags, where 1,-1 = use and 0 = don't use
%intfg_loops: rows = number of closed loops, columns = total no. of ints in network
intfg_loops = zeros(size(intfg_network,1),size(intfg,1));

for m=1:size(intfg_network,1);
    for n=1:nnz(~isnan(intfg_network(m,:)));
        if intfg_network(m,n)<0;
            intfg_loops(m,(abs(intfg_network(m,n))))=-1;
        else
            intfg_loops(m,intfg_network(m,n))=1;
        end
    end
end
logit(1);
save(loopname,'intfg_loops')

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
% wavelet_factory_1d: Create wavelet cascade from filters
% Usage
%    [Wop, filters] = wavelet_factory_1d(N, filter_options, scat_options)
% Input
%    N: The size of the signals to be transformed.
%    filter_options: The filter options, same as for filter_bank.
%    scat_options: General options to be passed to wavelet_1d and
%        wavelet_layer_1d. Contains scat_options.M, which determines the max-
%        imal order of the scattering transform.
% Output
%    Wop: A cell array of wavelet transforms needed for the scattering trans-
%       form.
%    filters: A cell array of the filters used in defining the wavelets.

function [Wop, filters] = renorm_wavelet_layer_factory_1d(N, filter_options, scat_options)
	filters = filter_bank(N, filter_options);

	for m = 0:scat_options.M
		filt_ind = min(numel(filters), m+1);
		Wop{m+1} = @(X)(renorm_layers_wavelet_1d(X, filters{filt_ind}, ...
			scat_options));
	end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function varargout = plot(varargin)
%PLOT  Plots the feasible region of a set of constraints
%
% p = plot(C,x,c,n,options)
%
% Note that only convex sets are allowed, or union of convex sets
% represented using binary variables (either defined explictly or
% introduced by YALMIP when modelling, e.g., mixed integer linear
% programming representable operators)
%
% C:  Constraint object
% x:  Plot variables [At most three variables]
% c:  color [double] ([r g b] format) or char from 'rymcgbk'
% n:  #vertices [double ] (default 100 in 2D and 300 otherwise)
% options: options structure from sdpsettings

% Get the onstraints
if nargin<1
    return
end

F = varargin{1};

if length(F)==0
    return;
end

if nargin < 5
    opts = sdpsettings('verbose',0);
else
    opts = varargin{5};
    if isempty(opts)
        opts = sdpsettings('verbose',0);
    end
end
opts.verbose = max(opts.verbose-1,0);

if any(is(F,'uncertain'))
    F = robustify(F,[],opts);
end

if nargin < 3
    color=['rymcgbk']';
else
    color = varargin{3};
    if isa(color,'sdpvar')
        error('The variables should be specified in the second argument.');
    end
    color = color(:)';
    if isempty(color)
        color=['rymcgbk']';
    end
end

% Plot onto this projection (at most in 3D)
if nargin < 2
    x = [];
else
    x = varargin{2};
    if ~isempty(x)
        if ~isa(x,'sdpvar')
            error('Second argument should be empty or an SDPVAR');
        end
        x = x(:);
        x = x(1:min(3,length(x)));
    end
end

if isempty(F)
    return
end

if any(is(F,'sos'))
    % Use image representation, feels safer (I'm not sure about the logic
    % in the code at the moment. Can the dualization mess up something
    % otherwise...)
    if ~(opts.sos.model == 1)
        opts.sos.model = 2;
    end
    % Assume the variables we are plotting are the parametric
    F = compilesos(F,[],opts,x);
end

% Create a model in YALMIPs low level format
% All we change later is the cost vector
%sol = solvesdp(F,sum(x),opts);
[model,recoverdata,diagnostic,internalmodel] = export(F,[],opts,[],[],0);
if isempty(internalmodel) | (~isempty(diagnostic) && diagnostic.problem)
    error('Could not create model. Can you actually solve problems with this model?')
end
internalmodel.options.saveduals = 0;
internalmodel.getsolvertime = 0;
internalmodel.options.dimacs = 0;

% Try to find a suitable set to plot
if isempty(x)
    if isempty(internalmodel.extended_variables) & isempty(internalmodel.aux_variables)
        x = depends(F);
        x = x(1:min(3,length(x)));
        localindex = 1;
        localindex = find(ismember(recoverdata.used_variables,x));
    else
        % not extended variables
        candidates = setdiff(1:length(internalmodel.c),[ internalmodel.aux_variables(:)' internalmodel.extended_variables(:)']);
        % Not nonlinear variables
        candidates = candidates(find(internalmodel.variabletype(candidates)==0));
        % Settle with this guess
        localindex = candidates(1:min(3,length(candidates)));
        x = localindex;
    end
else
    localindex = [];
    for i = 1:length(x)
        localindex = [localindex find(ismember(recoverdata.used_variables,getvariables(x(i))))];
    end
end

if nargin < 4
    if length(x)==3
        n = 300;
    else
        n = 100;
    end
else
    n = varargin{4};
    if isempty(n)
        if length(x)==3
            n = 300;
        else
            n = 100;
        end
    end
    if ~isa(n,'double')
        error('4th argument should be an integer>0');
    end
end

if ~isempty(internalmodel.integer_variables)
    error('PLOT can currently not display sets involving integer variables');
end

% Clean up in case we are plotting a relaxed set and solver has support for
% evaluation, but we don't intend to use that
if internalmodel.solver.evaluation==0
    internalmodel.evalMap = [];
end

% Use common function for lmi/plot and optimizer/plot
if nargout == 0
    plotInternalModel(internalmodel,x,n,localindex,color,opts);
else
    varargout{1} = plotInternalModel(internalmodel,x,n,localindex,color,opts);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function y = TpFerwerda(x)
%
%        y = TpFerwerda(x)
%
%
%       The gamma function used in Ferwerda TMO for Photopic levels
%
%       Input:
%           -x: a value
%
%       Output:
%           -y: application of the gamma function
%
%
%     Copyright (C) 2010-15 Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

t = log10(x);

if(t <= -2.6)
    y = -0.72;
else
    if(t >= 1.9)
        y = t - 1.255;
    else
        y = (0.249 * t + 0.65).^2.7 - 0.72;
    end
end

y = 10.^y;

end

author: banterle


text:
function [yUpdate, PInvUpdate]=infoFilterUpdate(yPred,PInvPred,z,R,H)
%%INFOFILTERUPDATE Perform the measurement update step in the standard
%                  linear information filter. Using an information filter
%                  means that instead of propagating a state x and its
%                  covariance matrix P, one propagates inv(P) and
%                  y=inv(P)*x.
%
%INPUTS: yPred The xDimX1 predicted information state. The information
%              state is the inverse covariance matrix times the target
%              state.
%     PInvPred The xDimXxDim inverse of the predicted state covariance
%              matrix.
%            z The zDimX1 vector measurement.
%            R The zDimXzDim measurement covariance matrix. This must be
%              positive definite.
%            H The zDimXxDim measurement matrix for a linear measurement
%              model. That is z=H*x+w, where w is measurement noise having
%              covariance matrix R.
%
%OUTPUTS: yUpdate The xDimX1 updated (posterior) information state vector.
%      PInvUpdate The updated xDimXxDim inverse state covariance matrix.
%
%The information filter is algebraically equivalent to the standard linear
%Kalman filter, but allows for track propagation with very uncertain
%states, such as when starting tracks. The implementation of the update
%given here is from the flow chart given in Appendix H of [1].
%
%The standard information filter has the implied linear measurement
%equation
%z=H*x+w
%where z is the measurement, w is the zero-mean Gaussian measurement noise
%with covariance matrix R and x is the true state. However, instead of
%propagating the true state x with its covariance matrix P, the information
%filter propagates
%y=inv(P)*x
%and instead of propagating the true covariance matrix P, it propagates
%PInv=inv(P)
%which means that the filter can be used even when PInv is singular.
%
%More information on information filtering is given in Chapter 7.2 of [2].
%
%EXAMPLE:
%In this example, we feed two measurements to the information filter and
%show that the result has a consistent NEES and the RMSE is the same as the
%second measurement.
% T=1;%Sample period.
% numRuns=1000;
% %Parameters for the noise process dynamics and measurement.
% H=[eye(2,2),zeros(2,2)];
% zDim=size(H,1);
% xDim=size(H,2);
% R=diag([40;40]);
% SR=chol(R,'lower');
%
% %Statistics for the initial state.
% x0Mean=[1e3;0;75;-50];
% x0Cov=diag([1e3^2;100^2;25^2;25^2]);
% x0S=chol(x0Cov,'lower');
%
% %Parameters for the state dynamics.
% q=processNoiseSuggest('PolyKal-ROT',9.8,1);
% F=FPolyKal(T,xDim,1);
% Q=QPolyKal(T,xDim,1,q);
% SQ=chol(Q,'lower');
%
% RMSEMeas=0;
% RMSE=0;
% NEES=0;
% for curRun=1:numRuns
%     %Draw the initial state.
%     x1True=x0Mean+x0S*randn(xDim,1);
%     %Get the first measurement.
%     z1=H*x1True+SR*randn(zDim,1);
%
%     x2True=F*x1True+SQ*randn(xDim,1);
%     z2=H*x2True+SR*randn(zDim,1);
%
%     %Two-point initialization.
%     y0=zeros(xDim,1);
%     PInv0=zeros(xDim,xDim);
%     [yUpdate,PInvUpdate]=infoFilterUpdate(y0,PInv0,z1,R,H);
%     [yPred, PInvPred]=infoFilterDiscPred(yUpdate,PInvUpdate,F,Q);
%     [yUpdate,PInvUpdate]=infoFilterUpdate(yPred,PInvPred,z2,R,H);
%     xEst2=PInvUpdate\yUpdate;
%
%     diff=xEst2-x2True;
%     NEES=NEES+diff'*PInvUpdate*diff;
%     RMSE=RMSE+sum(diff(1:2).^2);
%     diff=z2-x2True(1:2);
%     RMSEMeas=RMSEMeas+sum(diff(1:2).^2);
% end
% RMSEMeas=sqrt(RMSEMeas/numRuns)
% RMSE=sqrt(RMSE/numRuns)
% NEES=NEES/(xDim*numRuns)
%One will see that RMSEMeas equal RMSE, because with just two measurements,
%one cannot smmooth the position estimate syet. Additionally, NEES will be
%close to 1 indicating covariance consistency.
%
%REFERENCES:
%[1] David F. Crouse , "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%[2] Y. Bar-Shalom, X. R. Li, and T. Kirubarajan, Estimation with
%    Applications to Tracking and Navigation. New York: John Wiley and
%    Sons, Inc, 2001.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    yUpdate=yPred+H'/R*z;
    PInvUpdate=PInvPred+H'/R*H;
    %Ensure symmetry
    PInvUpdate=(PInvUpdate+PInvUpdate')/2;
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function metaImageWrite(img, fileName, varargin)
% Save image data using MetaImage file format (MHD).
%
%   metaImageWrite(IMG, FILENAME);
%   IMG is a matlab array, and FILENAME is the generic name (without
%   extension) of the metaimage file.
%   The functions tries to determine which parameters are appropriate for
%   the given image.
%
%   metaImageWrite(IMG, FILENAME, INFO);
%   Gives an additional structure as argument, containing additional
%   properties.
%
%   metaImageWrite(..., 'colorImage', true);
%   Forces writing as a RGB color array.
%
%   metaImageWrite(..., PARAM, VALUE);
%   Specifies additional properties as parameter key and value pairs. Case
%   is relevant. See http://www.itk.org/Wiki/MetaIO/Documentation for
%   details.
%
%   Example
%   [x y z] = meshgrid(1:80, 1:80, 1:80);
%   img = ((x-40).^2 + (y-40).^2 + (z-40).^2)<1000;
%   img = uint8(img*255);
%   metaImageWrite(img, 'ball');
%   metaImageWrite(img, 'ball.mhd');    % works also, same result
%   metaImageWrite(img, 'ball', 'ElementType', 'MET_USHORT', ...
%       'ElementByteOrderMSB', 'True')
%
%   See also
%     metaImageWrite, metaImageInfo
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2010-02-03,    using Matlab 7.9.0.529 (R2009b)
% http://www.pfl-cepia.inra.fr/index.php?page=slicer
% Copyright 2010 INRA - Cepia Software Platform.

%% Initialisations

% extract meta-information
info = struct;
if ~isempty(varargin)
    var = varargin{1};
    if isstruct(var)
        info = var;
        varargin(1) = [];
    end
end

% permute image dimensions to use x as first index
img = permute(img, [2 1 3:ndims(img)]);

% extract image dimension
dims = size(img);
nd = length(dims);

% flag for color image (false by default)
isColor = false;

%% Process file names

% extract file name parts
[path, name, ext] = fileparts(fileName);

% if extension is different from '.mhd', keep it as part of the name
if ~strcmp(ext, '.mhd')
    name = [name ext];
end

% header and binary file names
headerFileName = [name '.mhd'];
binaryFileName = [name '.raw'];

%% Create file info structure

% check main information are present
if ~isfield(info, 'ObjectType')
    info.ObjectType = 'Image';
end
if ~isfield(info, 'NDims')
    info.NDims = nd;
end
if ~isfield(info, 'DimSize')
    info.DimSize = dims;
end
if ~isfield(info, 'ElementType')
    info.ElementType = metaTypeToString(class(img));
end
if ~isfield(info, 'ElementDataFile')
    info.ElementDataFile = binaryFileName;
end

% parse additional input arguments and populate the "info" structure
while length(varargin) > 1
    key = varargin{1};
    if ~ischar(key)
        error('following parameter must be a string: %s', key);
    end

    % check function-specific options.
    if strcmpi(key, 'colorImage')
        isColor = varargin{2};
    else
        info.(key) = varargin{2};
    end

    varargin(1:2) = [];
end

% in case of color image, update the info structure
if isColor
    info.NDims = info.NDims - 1;
    info.DimSize = info.DimSize([1 2 4:end]);
    info.ElementNumberOfChannels = 3;
end

% ensure MSB info is written for data other than 8-bits
elementByteOrderDefined = isfield(info, 'ElementByteOrderMSB');
binaryDataByteOrderDefined = isfield(info, 'BinaryDataByteOrderMSB');
if ~ischar(img) && ~(elementByteOrderDefined || binaryDataByteOrderDefined)
    info.BinaryData = true;
    if ispc
        info.ElementByteOrderMSB = false;
    else
        info.ElementByteOrderMSB = true;
    end
end

% open file
f = fopen(fullfile(path, headerFileName), 'wt');

%% write minimal information tags

fprintf(f, '%s = %s\n', 'ObjectType', 'Image');
fprintf(f, '%s = %d\n', 'NDims', info.NDims);
fprintf(f, '%s =%s\n', 'DimSize', ...
    sprintf(repmat(' %d', 1, info.NDims), info.DimSize));
fprintf(f, '%s = %s\n', 'ElementType', info.ElementType);

%% write each additional tag

names = fieldnames(info);
for i = 1:length(names)
    name = names{i};

    % some tags are either in the very beginning or at the end, so they are
    % not processed in the loop
    if ismember(name, {'ObjectType', 'NDims', 'DimSize', 'ElementType', ...
            'ElementDataFile'})
        continue;
    end

    fprintf(f, '%s = %s\n', name, convertToString(info.(name)));
end

%% write binary data info into header

fprintf(f, '%s = %s\n', 'ElementDataFile', binaryFileName);

fclose(f);

%% Write binary data

% Specify little- or big-endian ordering
byteOrder = determineByteOrder(info);

% extract data precision
precision = parseMetaType(info.ElementType);

% open file for binary writing
f = fopen(fullfile(path, binaryFileName), 'wb');

% write image data in correct order
if ~isColor
    % write grayscale data
    fwrite(f, img(:), precision, 0, byteOrder);
else
    % first permute data to write RGB of each element, then x, y, and z
    data = permute(img, [3 1 2 4:ndims(img)]);
    fwrite(f, data(:), precision, 0, byteOrder);
end

% close binary file
fclose(f);

function type = parseMetaType(string)
% Convert MetaImage type string into a Matlab data type name.

switch string
    case 'MET_UCHAR'
        type = 'uint8';
    case 'MET_CHAR'
        type = 'int8';
    case 'MET_USHORT'
        type = 'uint16';
    case 'MET_SHORT'
        type = 'int16';
    case 'MET_UINT'
        type = 'uint32';
    case 'MET_INT'
        type = 'int32';
    case 'MET_FLOAT'
        type = 'single';
    case 'MET_DOUBLE'
        type = 'double';
    otherwise
        error('unknown element type in metaimage header: %s', string);
end

function string = metaTypeToString(type)
% Convert a Matlab data type name into MetaImage Type string.

switch type
    case 'int8'
        string = 'MET_CHAR';
    case 'uint8'
        string = 'MET_UCHAR';
    case 'int16'
        string = 'MET_SHORT';
    case 'uint16'
        string = 'MET_USHORT';
    case 'int32'
        string = 'MET_INT';
    case 'uint32'
        string = 'MET_UINT';
    case 'single'
        string = 'MET_FLOAT';
    case 'double'
        string = 'MET_DOUBLE';
    otherwise
        error('unknown pixel type: %s', type);
end

function string = convertToString(data)
% Utility function that converts logical or numeric data to char array.

N = length(data(:));

if ischar(data)
    string = data;

elseif islogical(data)
    strings = {'false', 'true'};
    pattern = ['%s' repmat(' %s', 1, N-1)];
    string = sprintf(pattern, strings{data(:)+1});

elseif isinteger(data)
    pattern = ['%d' repmat(' %d', 1, N-1)];
    string = sprintf(pattern, data(:));

elseif isfloat(data)
    pattern = ['%g' repmat(' %g', 1, N-1)];
    string = sprintf(pattern, data(:));

else
    error('Unknwon type of data');
end

function byteOrder = determineByteOrder(info)
% Identify the byte-order character to use for low-level fwrite function.

% default byte order given by system (native)
byteOrder = 'n';

% first check the ElementByteOrderMSB field
if isfield(info, 'ElementByteOrderMSB')
    if info.ElementByteOrderMSB
        byteOrder = 'b';
    else
        byteOrder = 'l';
    end
end

% also check the BinaryDataByteOrderMSB field
if isfield(info, 'BinaryDataByteOrderMSB')
    if info.BinaryDataByteOrderMSB
        byteOrder = 'b';
    else
        byteOrder = 'l';
    end
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function vrot=givapp(c,s,vin,k)
%givapp   apply a sequence of Givens rotations
%   input
%          c, s      vectors of length k-1 defining rotations
%          vin       vector of length k to which rotations are applied
%          k         k-1 = number of rotations
%   output
%          tranformed vector after rotations are applied
% called by gmres_r
%   IFISS function: HCE; 15 March 2005.

%
%  C. T. Kelley, July 10, 1994
%Copyright 1994 C. T. Kelley.
%Reproduced and distributed with permission of the copyright holder.
%
% This code comes with no guarantee or warranty of any kind.
%
%  function vrot=givapp(c, s, vin, k)
%
vrot=vin;
for i=1:k
    w1=c(i)*vrot(i)-s(i)*vrot(i+1);        % Change on next line, 6/3/97
    w2=s(i)*vrot(i)+conj(c(i))*vrot(i+1);  % w2=s(i)*vrot(i)+c(i)*vrot(i+1);
    vrot(i:i+1)=[w1,w2];
end

author: johannesgerer

repo: jburkardt-m


text:
function w = m_ilut_2dir(v,aparams,mparams)
%m_ilut_2dir     two-directional ILU preconditioning
%   w = m_ilut_2dir(v,aparams,mparams)
%   input
%          v            operand for preconditioning operator
%          aparams      structure defining coefficient matrix
%          mparams      structure defining preconditioning matrix
%                          factored splitting operators L*U, Lv*Uv
%                          permutation vectors p, ip
%   output
%          w            result of preconditioning operation
%
%   IFISS function: HCE; 28 February 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage

w = mparams.U\(mparams.L\v);
r = v - aparams.A*w;
rp = r(mparams.ip);
rp = mparams.Uv\(mparams.Lv\rp);
rpp = rp(mparams.p);
w = w + rpp;