text:
From smalone@engin.umich.edu Wed Jul 15 13:37:52 1998
To: "Jeffrey A. Fessler" <fessler@eecs.umich.edu>
Subject: Re: pwls_sor

Prof. Fessler,

I'm sending along two of my pwls_sor routines.  Both compute one iteration
of the routine.  One uses precomputed diag[G'*W*G].  The other does not.
I could combine the two, if you want, using the nargin command to
determine which one to use.  But I'm assuming your mostly interested in
the routine that uses the precomputed vector.  Let me know if you need any
changes to the code.

Brendhan

----pwls_sor with precomputed diag[G'*W*G]---

function [x, r] = pwls_sor(x, r, G, W, R, y, beta, GTWG)

%    --	PWLS+SOR Iterative Routine --
%	Copyright Brendhan Givens, University of Michigan, July 9, 1998
%
%	Details of this algorithm are described in
%	"Penalized Weighted Least-Squares Image Reconstruction
%	for Positron Emission Tomography", Jeffrey A. Fessler,
%	IEEE Transactions on Medical Imaging, Vol. 13, No. 2,
%	June 1994.
%
%	[x, r] = pwls_sor(x, r, G, W, R, y, beta, GTWG)
%
%	pwls_sor.m computes one iteration of the overdetermined symmetric
%	positive definite linear system (G'*W*G + beta*R)x = G'*W*y using
%	a sequential minimization procedure with a non-negativity
%	constraint.  Input matrices are assumed to be in sparse format.
%	Input residual vector must be initialized as r = y - G*x.  The
%	diag[G'*W*G] is precomputed for faster iterations.
%
%	input	x	REAL approximation vector
%		r	REAL residual vector
%		G	REAL FULL RANK matrix (sparse)
%		W	REAL diagonal weight matrix (sparse)
%		R	REAL SYMMETRIC penalty matrix (sparse)
%		y	REAL projection data
%		beta	REAL scalar
%		GTWG	REAL SYMMETRIC POSITIVE INDEFINITE vector
%			where elements gtwg(j) = G(:,j)'*W*G(:,j)
%
%	output	x	REAL NON-NEGATIVE approximation vector
%		r	REAL residual vector

[m,n] = size(G);			% initialization
w = diag(W);

for j = 1:n				% begin iteration
   d = ((G(:,j).*w)'*r - beta*R(:,j)'*x)/(GTWG(j) + beta*R(j,j));
   x_old = x(j);
   x(j) = max(0, x(j) + d);		% update with non-negativity
constraint
   r = r + (x_old - x(j))*G(:,j);	% update residual
end

----pwls_sor without precompute----

function [x, r] = pwls_sor(x, r, G, W, R, y, beta)

%    --	PWLS+SOR Iterative Routine --
%	Brendhan Givens, University of Michigan, July 9, 1998
%
%	Details of this algorithm are described in
%	"Penalized Weighted Least-Squares Image Reconstruction
%	for Positron Emission Tomography", Jeffrey A. Fessler,
%	IEEE Transactions on Medical Imaging, Vol. 13, No. 2,
%	June 1994.
%
%	[x, r] = pwls_sor(x, r, G, W, R, y, beta)
%
%	pwls_sor.m computes one iteration of the overdetermined symmetric
%	positive definite linear system (G'*W*G + beta*R)x = G'*W*y using
%	a sequential minimization procedure with a non-negativity
%	constraint.  Input matrices are assumed to be in sparse format.
%	Input residual vector must be initialized as r = y - G*x.
%
%	input	x	REAL approximation vector
%		r	REAL residual vector
%		G	REAL FULL RANK matrix (sparse)
%		W	REAL diagonal weight matrix (sparse)
%		R	REAL SYMMETRIC penalty matrix (sparse)
%		y	REAL projection data
%		beta	REAL scalar
%
%	output	x	REAL NON-NEGATIVE approximation vector
%		r	REAL residual vector

[m,n] = size(G);			% initialization
w = diag(W);

for j = 1:n				% begin iteration
   p = (G(:,j).*w);
   d = (p'*r - beta*R(:,j)'*x)/(p'*G(:,j) + beta*R(j,j));
   x_old = x(j);
   x(j) = max(0, x(j) + d);		% update with non-negativity
constraint
   r = r + (x_old - x(j))*G(:,j);	% update residual
end

On Wed, 15 Jul 1998, Jeffrey A. Fessler wrote:

>
> please email me your pwls_sor.m routine - the one that is just
> the update without the loop wrapper liked we discussed.  there is
> a nuc. eng. student who wants to use it.
> in fact, i think i'll put a collection of such routines on my web page!
> (with your name on this one of course).
> jf
>

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [wrongTable, absentModelTable, absentMapTable, duplicateTable] = compareModelMapFormulas(model, map, excelName)
% Checks the errors in a given map using a given base model by
% comparing the reactions formulas. As different errors can exist, the
% output is separated in 4 different tables that can later be exported
% into Excel spreadsheets (see commented lines at the end).
%
% USAGE:
%
%   [wrongTable, absentModelTable, absentMapTable, duplicateTable] = compareModelMapFormulas(model, map, excelName)
%
% INPUTS:
%   model:          COBRA structure of a model
%   map:            MATLAB structure of the map obtained from the
%                   function `transformXML2Map`.
%
% OPTIONAL INPUT:
%   excelName:      Name of the excel file in which to export the info
%
% OUTPUTS:
%   wrongTable:             Table containing the information on wrong
%                           reactions. The fields are arranged as
%                           followed:
%
%                            * rxnName - Name of the reaction in the map
%                            * rxnID - ID of the reaction in the map
%                            * modelFormula - Formula of the reaction from the model
%                            * mapFormula - Formula of the reaction from the map
%
%   absentModelTable:       Table containing the information on
%                           reactions present in the map but absent
%                           from the model. The fields are arranged as
%                           followed:
%
%                            * rxnName - Name of the reaction in the map
%                            * rxnID - ID of the reaction in the map
%                            * mapFormula - Formula of the reaction from the map
%
%   absentMapTable:         Table containing the information on
%                           reactions present in the model but absent
%                           from the map. The fields are arranged as
%                           followed:
%
%                            * rxnName - Name of the reaction in the model
%                            * modelFormula - Formula of the reaction from the model
%
%   duplicateTable:         Table containing the information on
%                           duplicated reactions in the map. The fields
%                           are arranged as followed:
%
%                            * rxnName - Name of the reaction in the model
%                            * rxnID - ID of the reaction in the map
%                            * modelFormula - Formula of the reaction from the model
%                            * mapFormula - Formula of the reaction from the map
%
% .. Author: - N.Sompairac - Institut Curie, Paris, 25/07/2017.

modelReactionNameList = model.rxns;  % Getting names from model and map
mapReactionNameList = map.rxnName;

% Getting the formulas from the model and the map for further comparison
modelFormulasList = printRxnFormula(model, modelReactionNameList,0);
[mapFormulasList, mapReactionNameList] = mapFormula(map, mapReactionNameList);

% Deleting the stoechiometric numbers from the model formulas
modelFormulasList = regexprep(modelFormulasList, '[0-9.]+ ', '');

% Initialising lists that will contain corresponding info on reactions
wrong = 1;
wrongList = [];
dupl = 1;
duplicateList = [];
abs = 1;
absentMapList = [];

% Looping over the model's reaction names
for rxn = 1:length(modelReactionNameList)
    % Test if the reaction name is contained in the map
    if any(strcmp(modelReactionNameList{rxn}, mapReactionNameList))
        % Getting the index of the model reaction name in the map list
        index = find(strcmp(modelReactionNameList{rxn}, mapReactionNameList));
        % Test if there is only one reaction with this name in the map
        if length(index) == 1
            % Deleting the stoechiometric numbers from the model formula
            % Splitting the model formula for further comparison
            modelFormulaSplit = strsplit(modelFormulasList{rxn}, {'<=>', '->'});
            leftModel = strtrim(strsplit(modelFormulaSplit{1}, '+'));
            rightModel = strtrim(strsplit(modelFormulaSplit{2}, '+'));
            % Splitting the map formula for further comparison
            mapFormulaSplit = strsplit(mapFormulasList{index}, {'<=>', '->'});
            leftMap = strtrim(strsplit(mapFormulaSplit{1}, '+'));
            rightMap = strtrim(strsplit(mapFormulaSplit{2}, '+'));
            % Testing if the formulas are different and storing the info
            leftTest = setxor(leftModel, leftMap);
            rightTest = setxor(rightModel, rightMap);
            if ~isempty(leftTest) || ~isempty(rightTest)
                wrongList.name{wrong} = mapReactionNameList{index};
                wrongList.ID{wrong} = map.rxnID{strcmp(modelReactionNameList{rxn}, map.rxnName)};
                wrongList.modelFormula{wrong} = modelFormulasList{rxn};
                wrongList.mapFormula{wrong} = mapFormulasList{index};
                wrong = wrong + 1;
            end
            % Case where a reaction name is duplicated in the map
        else
            % Finding the IDs of the duplicated reactions in the map
            duplicateIDs = map.rxnID(strcmp(modelReactionNameList{rxn}, map.rxnName));
            % Looping over the duplicates to get the relevant info
            for d = 1:length(index)
                duplicateList.name{dupl} = modelReactionNameList{rxn};
                duplicateList.modelFormula{dupl} = modelFormulasList{rxn};
                duplicateList.mapFormula{dupl} = mapFormulasList{index(d)};
                duplicateList.ID{dupl} = duplicateIDs{d};
                dupl = dupl + 1;
            end
        end
        % Case where reactions are absent in the map and present in the model
    else
        absentMapList.name{abs} = modelReactionNameList{rxn};
        absentMapList.modelFormula{abs} = modelFormulasList{rxn};
        abs = abs + 1;
    end
end

if ~isempty(wrongList)

    wrongTable = table(wrongList.name', wrongList.ID', wrongList.modelFormula', wrongList.mapFormula');
    wrongTable.Properties.VariableNames = {'rxnName', 'rxnID', 'modelFormula', 'mapFormula'};

else

    wrongTable = [];

end

if ~isempty(duplicateList)
    duplicateTable = table(duplicateList.name', duplicateList.ID', duplicateList.modelFormula', duplicateList.mapFormula');
    duplicateTable.Properties.VariableNames = {'rxnName', 'rxnID', 'modelFormula', 'mapFormula'};
else
    duplicateTable = [];
end

if ~isempty(absentMapList)
    absentMapTable = table(absentMapList.name', absentMapList.modelFormula');
    absentMapTable.Properties.VariableNames = {'rxnName', 'modelFormula'};
else
    absentMapTable = [];
end

% Finding reaction names in the map that are not present in the model
differentMapRxnNamesList = setdiff(mapReactionNameList, modelReactionNameList);

if ~isempty(differentMapRxnNamesList)
    % Finding reaction names in the map in case of multiple presence
    differentMapRxnNamesList = map.rxnName(ismember(map.rxnName, differentMapRxnNamesList));
    % Finding reaction ID in the map that are not present in the model
    differentMapRxnIdList = map.rxnID(ismember(map.rxnName, differentMapRxnNamesList));
    % Finding reaction formulas in the map that are not present in the model
    differentMapRxnFormulasList = mapFormulasList(ismember(mapReactionNameList, differentMapRxnNamesList));

    % Storing the relevant info on missing reactions in the model from the map
    absentModelTable = table(differentMapRxnNamesList, differentMapRxnIdList, differentMapRxnFormulasList);
    absentModelTable.Properties.VariableNames = {'rxnName', 'rxnID', 'mapFormula'};
else
    absentModelTable = [];
end

if nargin == 3
    % Commented part to use a possible Excel output.
    fileName = excelName;
    warning('off', 'MATLAB:xlswrite:AddSheet');
    if ~isempty(wrongTable)
        writetable(wrongTable, fileName, 'Sheet', 'wrongReactions')
    end
    if ~isempty(absentMapTable)
        writetable(absentMapTable, fileName, 'Sheet', 'absentFromMapReactions')
    end
    if ~isempty(absentModelTable)
        writetable(absentModelTable, fileName, 'Sheet', 'absentFromModelReactions')
    end
    if ~isempty(duplicateTable)
        writetable(duplicateTable, fileName, 'Sheet', 'duplicatedReactions')
    end
end

author: opencobra

repo: cobratoolbox


text:
function [ X2 Y2 ] = LKTrackWrapper( imgseq,X1,Y1 )
% [ X2 Y2 ] = LKTrackWrapper( imgseq,X1,Y1 )
%	calls LKTrackPyr for each 2 neighboring frames in imgseq.
%	imgseq(M-by-N-by-F) are F grayscale images, it can be either from
%	a image sequence or a video(like from VideoReader).
%	if X1 and Y1 are not specified, function calls corner_ST to find
%	maxPtsNum corners. if X2 and Y2 are not required to output, function
%	will show these points in a figure. You should press a key to
%	start the show.

if nargin > 0 % image sequence or video
	[M N imgNum] = size(imgseq);
	maxPtsNum = 50;
	if ~exist('Y1','var')
		[Y1 X1] = corner_ST(imgseq(:,:,1),maxPtsNum);
		borderTh = 10;
		discard = Y1<borderTh | Y1>M-borderTh |...
			X1<borderTh | X1>N-borderTh; % discard points on the image border
		Y1 = Y1(~discard);
		X1 = X1(~discard);
	end

	X2 = zeros(length(X1),imgNum);
	Y2 = zeros(length(X1),imgNum);
	X2(:,1) = X1;
	Y2(:,1) = Y1;

	for p = 2:imgNum
		[X2(:,p) Y2(:,p)] = LKTrackPyr(imgseq(:,:,p-1),imgseq(:,:,p),...
			X2(:,p-1),Y2(:,p-1));
	end

	if nargout == 0
		figure,pause % show the points
		sc = min(M,N)/30; % scale of the line indicating the speed
		for p = 1:imgNum
			imshow(imgseq(:,:,p)),hold on
			X2p = X2(:,p); Y2p = Y2(:,p);
			plot(X2p,Y2p,'go');
			if p > 1 % draw speed lines
				plot([X2p,X2p+(X2p-X2pl)*sc]',...
					[Y2p,Y2p+(Y2p-Y2pl)*sc]','m-');
			end
			pause;
			X2pl = X2p; Y2pl = Y2p;
		end
	end

else % webcam, not finished

	vid = videoinput('winvideo',1); % specify your own name and index
	vid.TriggerRepeat = Inf;
	vid.FrameGrabInterval = 3;
	vid.ReturnedColorSpace='grayscale'; % 'rgb'
	figure
	while vid.FramesAcquired<=100
		data = getdata(vid,1);
		imshow(data);
		drawnow
	end
	stop(vid);
	delete(vid);

end

author: Sable

repo: mcbench-benchmarks


text:
function readJV(file)
[nv,nf]=nnu(file);
fprintf('number of point:%d\n',nv);
fprintf('number of face:%d\n',nf);
fid=fopen(file);
fid2=fopen('point.off','w');
while 1
    tline = fgetl(fid);
    if ~ischar(tline),   break,   end
    %disp(tline)
    n=length(tline);
    if n > 0
        header=findstr(tline,'#JV');
        J=isempty(header);
        if J==1
        disp(tline)
        fprintf(fid2,'%s \n',tline);
        end
    end
end
fclose(fid);
fclose(fid2);
function [nv,nf]=nnu(file)
fid = fopen(file);
C=0;
while 1
    tline = fgetl(fid);
    C=C+1;
    if ~ischar(tline),   break,   end
    if C==3
        CC = textscan(tline, '%f %f %f');
        break;
    end
    %disp(tline)
end
fclose(fid);
nv=CC{1,1};
nf=CC{1,2};


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

% get calibrated data
load 'MODEL_Hist';

% Discretization parameters
NTime = 12;                    % number of time steps
NSim = 10000;                    % number of simulations
NBatches = 1;                  % number of batches

NHist = 251;                     % number of historic scenarios

NModels = 10;                  % number of models

% Options
%K = 6800;                                % Strike
%C = 1;                                  % call (1) or put (0)
LF = 0; LC = 0.04; GF = 0; GC = 0.15;   % cliquet params
L = 80; U = 120;
T= 1;                                   % maturity of option

%option1 = 'Arithmetic Asian Call';
%option2 = 'Arithmetic Asian Put';
%exotic_price_c = @(x,y,r) Price_ArithmeticAsian(x,y,1,r,T);
%exotic_price_p = @(x,y,r) Price_ArithmeticAsian(x,y,0,r,T);

% exotic_price = @(x,r) Price_ArithmeticCliquet(x,LF,LC,GF,GC,r,T);
option1 = 'FixedStrike Lookback Call';
option2 = 'FixedStrike Lookback Put';
exotic_price_c = @(x,y,r) Price_LookBackFixedStrike(x,y,1,r,T);
exotic_price_p = @(x,y,r) Price_LookBackFixedStrike(x,y,0,r,T);

%option1 = 'FloatingStrike Lookback Call';
%option2 = 'FloatingStrike Lookback Put';
%exotic_price_c = @(x,y,r) Price_LookBackFloatingStrike(x,y,1,r,T);
%exotic_price_p = @(x,y,r) Price_LookBackFloatingStrike(x,y,0,r,T);

%exotic_price = @(x,r) Price_KnockOut(x,x,L,U,C,r,T);
vanilla_price_c = @(x,y,r) Price_CallPut(x,y,1,r,T);
vanilla_price_p = @(x,y,r) Price_CallPut(x,y,0,r,T);

% output
ep1c = zeros(NHist,NModels); %ep2c = ep1c; ep3c = ep1c;
ep1p = ep1c; %ep2p = ep1c; ep3p = ep1c;

vp_c = ep1c; vp_p = ep1c;
qp_c = ep1c; qp_p = ep1c;

vc = zeros(NHist,1); vp = vc;

% Black Scholes
for l = 1: NHist
    pathS = MC_B(Spot(l),r1y(l),0,T,BS_Par_Hist(l),NTime,NSim,NBatches);
    ep1c(l,1) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,1) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,1) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,1) = vanilla_price_p(pathS,Spot(l),r1y(l));
    [vc(l), vp(l)] = blsprice(Spot(l), Spot(l), r1y(l), T, BS_Par_Hist(l),0);
    qp_c(l,1) = ep1c(l,1)/vp_c(l,1);
    qp_p(l,1) = ep1p(l,1)/vp_p(l,1);
end

% Merton
for l = 1: NHist
    pathS = MC_M(Spot(l),r1y(l),0,T,...
        MJ_Par_Hist(l,1),MJ_Par_Hist(l,2),MJ_Par_Hist(l,3),MJ_Par_Hist(l,4),...
        NTime,NSim,NBatches);
    ep1c(l,2) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,2) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,2) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,2) = vanilla_price_p(pathS,Spot(l),r1y(l));
    qp_c(l,2) = ep1c(l,2)/vp_c(l,2);
    qp_p(l,2) = ep1p(l,2)/vp_p(l,2);
end

% Heston
for l = 1: NHist
    [pathS, pathV] = MC_QE(Spot(l),r1y(l),0,T,...
        HESTON_Par_Hist(l,1),HESTON_Par_Hist(l,2),HESTON_Par_Hist(l,3),...
        HESTON_Par_Hist(l,4),HESTON_Par_Hist(l,5),...
        NTime,NSim,NBatches);
    ep1c(l,3) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,3) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,3) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,3) = vanilla_price_p(pathS,Spot(l),r1y(l));
    qp_c(l,3) = ep1c(l,3)/vp_c(l,3);
    qp_p(l,3) = ep1p(l,3)/vp_p(l,3);
end

% Bates
for l = 1: NHist
    [pathS, pathV] = MC_QE_j(Spot(l),r1y(l),0,T,...
        BATES_Par_Hist(l,1),BATES_Par_Hist(l,2),BATES_Par_Hist(l,4),BATES_Par_Hist(l,5),...
        BATES_Par_Hist(l,3),BATES_Par_Hist(l,7),BATES_Par_Hist(l,7),BATES_Par_Hist(l,8),...
        NTime,NSim,NBatches);
    ep1c(l,4) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,4) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,4) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,4) = vanilla_price_p(pathS,Spot(l),r1y(l));
    qp_c(l,4) = ep1c(l,4)/vp_c(l,4);
    qp_p(l,4) = ep1p(l,4)/vp_p(l,4);
end

% % Variance Gamma CGM
% for l = 1:NHist
%     pathS  = MC_VG_CGM(Spot(l),r1y(l),0,T,...
%         VG_Par_Hist(l,1),VG_Par_Hist(l,2),VG_Par_Hist(l,3),...
%         NTime,NSim,NBatches);
%     ep1c(l,5) = exotic_price_c(pathS,Spot(l),r1y(l));
%     ep1p(l,5) = exotic_price_p(pathS,Spot(l),r1y(l));
%     vp_c(l,5) = vanilla_price_c(pathS,Spot(l),r1y(l));
%     vp_p(l,5) = vanilla_price_p(pathS,Spot(l),r1y(l));
%     qp_c(l,5) = ep1c(l,5)/vp_c(l,1);
%     qp_p(l,5) = ep1p(l,5)/vp_p(l,1);
% end
%
% % Normal Inverse Gaussian
% for l = 1:NHist
%     pathS  = MC_NIG(Spot(l),r1y(l),0,T,...
%         NIG_Par_Hist(l,1),NIG_Par_Hist(l,2),NIG_Par_Hist(l,3),...
%         NTime,NSim,NBatches);
%     ep1c(l,6) = exotic_price_c(pathS,Spot(l),r1y(l));
%     ep1p(l,6) = exotic_price_p(pathS,Spot(l),r1y(l));
%     vp_c(l,6) = vanilla_price_c(pathS,Spot(l),r1y(l));
%     vp_p(l,6) = vanilla_price_p(pathS,Spot(l),r1y(l));
%     qp_c(l,6) = ep1c(l,1)/vp_c(l,1);
%     qp_p(l,6) = ep1p(l,6)/vp_p(l,1);
% end
%
% % Variance Gamma - GOU
% for l = 1:NHist
%     pathS  = MC_VGGOU(Spot(l),r1y(l),0,T,...
%         VGGOU_Par_Hist(l,1), VGGOU_Par_Hist(l,2), VGGOU_Par_Hist(l,3),...
%         VGGOU_Par_Hist(l,6), VGGOU_Par_Hist(l,4), VGGOU_Par_Hist(l,5),...
%         NTime,NSim,NBatches);
%     ep1c(l,7) = exotic_price_c(pathS,Spot(l),r1y(l));
%     ep1p(l,7) = exotic_price_p(pathS,Spot(l),r1y(l));
%     vp_c(l,1) = vanilla_price_c(pathS,Spot(l),r1y(l));
%     vp_p(l,7) = vanilla_price_p(pathS,Spot(l),r1y(l));
%     qp_c(l,7) = ep1c(l,7)/vp_c(l,1);
%     qp_p(l,7) = ep1p(l,7)/vp_p(l,1);
% end
%
% % % Variance Gamma - CIR
% for l = 1: NHist
%     pathS =  MC_VGCIR(Spot(l),r1y(l),0,T,...
%         VGCIR_Par_Hist(l,1),VGCIR_Par_Hist(l,2),VGCIR_Par_Hist(l,3), ...
%         VGCIR_Par_Hist(l,4),VGCIR_Par_Hist(l,5),VGCIR_Par_Hist(l,6),...
%         NTime,NSim,NBatches);
%     ep1c(l,8) = exotic_price_c(pathS,Spot(l),r1y(l));
%     ep1p(l,8) = exotic_price_p(pathS,Spot(l),r1y(l));
%     vp_c(l,8) = vanilla_price_c(pathS,Spot(l),r1y(l));
%     vp_p(l,8) = vanilla_price_p(pathS,Spot(l),r1y(l));
%     qp_c(l,8) = ep1c(l,8)/vp_c(l,1);
%     qp_p(l,8) = ep1p(l,8)/vp_p(l,1);
% end
%
% % NIG - GOU
% for l = 1: NHist
%     pathS =  MC_NIGGOU(Spot(l),r1y(l),0,T,...
%         NIGGOU_Par_Hist(l,1),NIGGOU_Par_Hist(l,2),NIGGOU_Par_Hist(l,3), ...
%         NIGGOU_Par_Hist(l,4),NIGGOU_Par_Hist(l,5),NIGGOU_Par_Hist(l,6),NTime,NSim,NBatches);
%     ep1c(l,9) = exotic_price_c(pathS,Spot(l),r1y(l));
%     ep1p(l,9) = exotic_price_p(pathS,Spot(l),r1y(l));
%     vp_c(l,9) = vanilla_price_c(pathS,Spot(l),r1y(l));
%     vp_p(l,9) = vanilla_price_p(pathS,Spot(l),r1y(l));
%     qp_c(l,9) = ep1c(l,9)/vp_c(l,1);
%     qp_p(l,9) = ep1p(l,9)/vp_p(l,1);
% end
%
% % NIG - CIR
% for l = 1: NHist
%     pathS =  MC_NIGCIR(Spot(l),r1y(l),0,T,...
%         NIGCIR_Par_Hist(l,1),NIGCIR_Par_Hist(l,2),NIGCIR_Par_Hist(l,3), ...
%         NIGCIR_Par_Hist(l,5),NIGCIR_Par_Hist(l,6),NIGCIR_Par_Hist(l,4),NTime,NSim,NBatches);
%     ep1c(l,10) = exotic_price_c(pathS,Spot(l),r1y(l));
%     ep1p(l,10) = exotic_price_p(pathS,Spot(l),r1y(l));
%     vp_c(l,10) = vanilla_price_c(pathS,Spot(l),r1y(l));
%     vp_p(l,10) = vanilla_price_p(pathS,Spot(l),r1y(l));
%     qp_c(l,10) = ep1c(l,10)/vp_c(l,1);
%     qp_p(l,10) = ep1p(l,10)/vp_p(l,1);
% end
m1 = 'BS'; m2 = 'MJ'; m3 = 'Heston'; m4 = 'Bates';

plotit3(ep1c,option1,m1, m2, m3, m4);
plotit3(ep1p,option2,m1, m2, m3, m4);
plotit3(qp_c,strcat('Normalized ',option1),m1, m2, m3, m4);
plotit3(qp_p,strcat('Normalized ',option2),m1, m2, m3, m4);

modelriskc = repmat(max(ep1c,[],2),1,NModels)-ep1c;
modelriskp = -repmat(max(ep1p,[],2),1,NModels)+ep1p;


text:
function dist = sphere_distance_xyz ( xyz1, xyz2 )

%*****************************************************************************80
%
%% SPHERE_DISTANCE_XYZ computes great circle distances on a sphere.
%
%  Discussion:
%
%    XYZ coordinates are used.
%
%    We assume the points XYZ1 and XYZ2 lie on the same sphere.
%
%    This computation is a special form of the Vincenty formula.
%    It should be less sensitive to errors associated with very small
%    or very large angular separations.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    "Great-circle distance",
%    Wikipedia.
%
%  Parameters:
%
%    Input, real XYZ1(3), the coordinates of the first point.
%
%    Input, real XYZ2(3), the coordinates of the second point.
%
%    Output, real DIST, the great circle distance between the points.
%
  r = norm ( xyz1 );

  lat1 = r8_asin ( xyz1(3) );
  lon1 = r8_atan ( xyz1(2), xyz1(1) );

  lat2 = r8_asin ( xyz2(3) );
  lon2 = r8_atan ( xyz2(2), xyz2(1) );

  top = ( cos ( lat2 ) * sin ( lon1 - lon2 ) ).^2 ...
      + ( cos ( lat1 ) * sin ( lat2 ) ...
      -   sin ( lat1 ) * cos ( lat2 ) * cos ( lon1 - lon2 ) ).^2;

  top = sqrt ( top );

  bot = sin ( lat1 ) * sin ( lat2 ) ...
      + cos ( lat1 ) * cos ( lat2 ) * cos ( lon1 - lon2 );

  dist = r * atan2 ( top, bot );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function as_l_d
%
%  Deletes the global data that has been generated by 'as_l_i'.
%
%  Calling sequence:
%
%    as_l_d
%
%
%  LYAPACK 1.0 (Thilo Penzl, May 1999)

global LP_U

clear global LP_U


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function y = asian_bintree(S, df, p, K, type)
% pricing an american asian option using a binomial tree
n = size(S,2)-1;
v = cell(1, n+1);               % option value

tmp = repmat(1:n, length(S), 1);
average_val = cumsum(S(:,2:end), 2) ./ tmp; % average value

if type == 0
   cp = 1;
else
    cp = -1;
end

v{n+1} = max(cp*(average_val(:,end)-K), 0); % payoff at t_N

iVec = 1:length(S);

for i = n:-1:2
    Jset = 1:2:2^i;
    expected_val = p * v{i+1}(Jset) * df ...
        + (1-p) * v{i+1}(1+Jset) * df;

    iVec = iVec(1:length(iVec)/2) * 2;

    v{i} = max(cp*(average_val(iVec,i-1)-K),0);   % payoff at t_i
    v{i}=max(v{i}, expected_val);                 % option value t_i
end

y = p * v{2}(1) * df + (1-p) * v{2}(2) * df;      % option value t_0


text:
function sphere_triangle_quad_test04 ( )

%*****************************************************************************80
%
%% SPHERE_TRIANGLE_TEST04 tests SPHERE01_TRIANGLE_QUAD_ICOS1M.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2014
%
%  Author:
%
%    John Burkardt
%
  e_test = [ ...
    0, 0, 0; ...
    1, 0, 0; ...
    0, 1, 0; ...
    0, 0, 1; ...
    2, 0, 0; ...
    0, 2, 2; ...
    2, 2, 2; ...
    0, 2, 4; ...
    0, 0, 6; ...
    1, 2, 4; ...
    2, 4, 2; ...
    6, 2, 0; ...
    0, 0, 8; ...
    6, 0, 4; ...
    4, 6, 2; ...
    2, 4, 8; ...
   16, 0, 0 ]';
  n_mc1 = 1000;
  n_mc2 = 10000;
  n_mc3 = 100000;

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_TRIANGLE_QUAD_TEST04\n' );
  fprintf ( 1, '  SPHERE01_TRIANGLE_QUAD_ICOS1M approximates the\n' );
  fprintf ( 1, '  integral of a function over a spherical triangle on\n' );
  fprintf ( 1, '  the surface of the unit sphere using a midpoint rule.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We do not have an exact result, so we compare each\n' );
  fprintf ( 1, '  estimate to the final one.\n' );
%
%  Choose three points at random to define a spherical triangle.
%
  [ v1, seed ] = sphere01_sample ( 1, seed );
  [ v2, seed ] = sphere01_sample ( 1, seed );
  [ v3, seed ] = sphere01_sample ( 1, seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Vertices of random spherical triangle:\n' );
  fprintf ( 1, '\n' );
  r8vec_transpose_print ( 3, v1, '  V1:' );
  r8vec_transpose_print ( 3, v2, '  V2:' );
  r8vec_transpose_print ( 3, v3, '  V3:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  FACTOR   N   RESULT\n' );

  for j = 1 : 17

    e(1:3) = e_test(1:3,j);

    e = polyterm_exponent ( 'SET', e );

    polyterm_exponent ( 'PRINT', e );

%   factor = 2 ^ 11;
    factor = 2 ^ 7;
    [ best, node_num ] = sphere01_triangle_quad_icos1m ( v1, v2, v3, factor, ...
      @polyterm_value_3d );

    factor = 1;
    for factor_log = 0 : 5

      [ result, node_num ] = sphere01_triangle_quad_icos1m ( v1, v2, v3, ...
        factor, @polyterm_value_3d );

      error = abs ( result - best );

      fprintf ( 1, '  %4d  %8d  %16.8g  %10.2e\n', ...
        factor, node_num, result, error );

      factor = factor * 2;

    end

  end

  return
end


text:
function data = LoadBinary(filename,varargin)

%LoadBinary - Load data from a multiplexed binary file.
%
%  Reading a subset of the data can be done in two different manners: either
%  by specifying start time and duration (more intuitive), or by indicating
%  the position and size of the subset in terms of number of samples per
%  channel (more accurate).
%
%  USAGE
%
%    data = LoadBinary(filename,<options>)
%
%    filename       file to read
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'frequency'   sampling rate (in Hz, default = 20kHz)
%     'start'       position to start reading (in s, default = 0)
%     'duration'    duration to read (in s, default = Inf)
%     'offset'      position to start reading (in samples per channel,
%                   default = 0)
%     'samples'     number of samples (per channel) to read (default = Inf)
%     'nChannels'   number of data channels in the file (default = 1)
%     'channels'    channels to read (default = all)
%     'precision'   sample precision (default = 'int16')
%     'skip'        number of bytes to skip after each value is read
%                   (default = 0)
%     'downsample'  factor by which to downample by (default = 1)
%    =========================================================================

% Copyright (C) 2004-2011 by Michaël Zugaro
%Modified by DLevenstein 2016 to include downsampling
%
% 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.

% Default values
nChannels = 1;
precision = 'int16';
%error('this function is now deprecated, use bz_LoadBinary.m')
%warning(['this function assumes int16 precision, if your file is not int16 use load/resample.m'])
skip = 0;
frequency = 20000;
channels = [];
start = 0;
duration = Inf;
offset = 0;
nSamplesPerChannel = Inf;
time = false;
samples = false;
downsamplefactor = 1;

if nargin < 1 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help LoadBinary'' for details).');
end

% Parse options
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+3) ' is not a property (type ''help LoadBinary'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'frequency',
			frequency = varargin{i+1};
			if ~isdscalar(frequency,'>0'),
				error('Incorrect value for property ''frequency'' (type ''help LoadBinary'' for details).');
			end
		case 'start',
			start = varargin{i+1};
			if ~isdscalar(start),
				error('Incorrect value for property ''start'' (type ''help LoadBinary'' for details).');
			end
			if start < 0, start = 0; end
			time = true;
		case 'duration',
			duration = varargin{i+1};
			if ~isdscalar(duration,'>=0'),
				error('Incorrect value for property ''duration'' (type ''help LoadBinary'' for details).');
			end
			time = true;
		case 'offset',
			offset = varargin{i+1};
			if ~isiscalar(offset),
				error('Incorrect value for property ''offset'' (type ''help LoadBinary'' for details).');
			end
			if offset < 0, offset = 0; end
			samples = true;
		case 'samples',
			nSamplesPerChannel = varargin{i+1};
			if ~isdscalar(nSamplesPerChannel,'>=0'),
				error('Incorrect value for property ''samples'' (type ''help LoadBinary'' for details).');
			end
			samples = true;
		case 'nchannels',
			nChannels = varargin{i+1};
			if ~isiscalar(nChannels,'>0'),
				error('Incorrect value for property ''nChannels'' (type ''help LoadBinary'' for details).');
			end
		case 'channels',
			channels = varargin{i+1};
			if ~isivector(channels,'>0'),
				error('Incorrect value for property ''channels'' (type ''help LoadBinary'' for details).');
			end
		case 'precision',
			precision = varargin{i+1};
			if ~isa(precision,'char'),
				error('Incorrect value for property ''precision'' (type ''help LoadBinary'' for details).');
			end
		case 'skip',
			skip = varargin{i+1};
			if ~isiscalar(skip,'>=0'),
				error('Incorrect value for property ''skip'' (type ''help LoadBinary'' for details).');
            end
		case 'downsample',
			downsamplefactor = varargin{i+1};
			if ~isiscalar(downsamplefactor,'>=0'),
				error('Incorrect value for property ''downsample'' (type ''help LoadBinary'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help LoadBinary'' for details).']);
	end
end

% Either start+duration, or offset+size
if time && samples,
	error(['Data subset can be specified either in time or in samples, but not both (type ''help LoadBinary'' for details).']);
end

% By default, load all channels
if isempty(channels),
	channels = 1:nChannels;
end

% Check consistency between channel IDs and number of channels
if any(channels>nChannels),
	error('Cannot load specified channels (listed channel IDs inconsistent with total number of channels).');
end

% Open file
if ~exist(filename),
	error(['File ''' filename ''' not found.']);
end
f = fopen(filename,'r');
if f == -1,
	error(['Cannot read ' filename ' (insufficient access rights?).']);
end

% Size of one data point (in bytes)
sampleSize = 0;
switch precision,
	case {'uchar','unsigned char','schar','signed char','int8','integer*1','uint8','integer*1'},
		sampleSize = 1;
	case {'int16','integer*2','uint16','integer*2'},
		sampleSize = 2;
	case {'int32','integer*4','uint32','integer*4','single','real*4','float32','real*4'},
		sampleSize = 4;
	case {'int64','integer*8','uint64','integer*8','double','real*8','float64','real*8'},
		sampleSize = 8;
end

% Position and number of samples (per channel) of the data subset
if time,
	dataOffset = floor(start*frequency)*nChannels*sampleSize;
	nSamplesPerChannel = floor((duration*frequency));
else
	dataOffset = offset*nChannels*sampleSize;
end

% Position file index for reading
status = fseek(f,dataOffset,'bof');
if status ~= 0,
	fclose(f);
	error('Could not start reading (possible reasons include trying to read past the end of the file).');
end

% Determine total number of samples in file
fileStart = ftell(f);
status = fseek(f,0,'eof');
if status ~= 0,
	fclose(f);
	error('Error reading the data file (possible reasons include trying to read past the end of the file).');
end
fileStop = ftell(f);
% (floor in case all channels do not have the same number of samples)
maxNSamplesPerChannel = floor(((fileStop-fileStart)/nChannels/sampleSize));
frewind(f);
status = fseek(f,dataOffset,'bof');
if status ~= 0,
	fclose(f);
	error('Could not start reading (possible reasons include trying to read past the end of the file).');
end

if isinf(nSamplesPerChannel) || nSamplesPerChannel > maxNSamplesPerChannel,
	nSamplesPerChannel = maxNSamplesPerChannel;
end

if downsamplefactor>1
%     precision = [num2str(nChannels),'*',precision]; % this line is
%     incorrect, the precision variable is a string that does not depend on
%     the number of channels
    skip = nChannels*(downsamplefactor-1)*sampleSize;
    nSamplesPerChannel = floor(nSamplesPerChannel./downsamplefactor);
end

% For large amounts of data, read chunk by chunk
maxSamplesPerChunk = 10000;
nSamples = nSamplesPerChannel*nChannels;
if nSamples <= maxSamplesPerChunk,
	data = LoadChunk(f,nChannels,channels,nSamples/nChannels,precision,skip);
else
	% Determine chunk duration and number of chunks
	nSamplesPerChunk = floor(maxSamplesPerChunk/nChannels)*nChannels;
	nChunks = floor(nSamples/nSamplesPerChunk);
	% Preallocate memory
	data = zeros(nSamplesPerChannel,length(channels),precision);
	% Read all chunks
	i = 1;
	for j = 1:nChunks,
		d = LoadChunk(f,nChannels,channels,nSamplesPerChunk/nChannels,precision,skip);
		[m,n] = size(d);
		if m == 0, break; end
		data(i:i+m-1,:) = d;
		i = i+m;
	end
	% If the data size is not a multiple of the chunk size, read the remainder
	remainder = nSamples - nChunks*nSamplesPerChunk;
	if remainder ~= 0,
		d = LoadChunk(f,nChannels,channels,remainder/nChannels,precision,skip);
		[m,n] = size(d);
		if m ~= 0,
			data(i:i+m-1,:) = d;
		end
	end
end

fclose(f);

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

function data = LoadChunk(fid,nChannels,channels,nSamples,precision,skip)

if skip ~= 0,
	data = fread(fid,[nChannels nSamples],[num2str(nChannels),'*',precision '=>' precision],skip);
    %data = fread(fid,[nChannels nSamples],[num2str(nChannels),'*',precision],nChannels*skip);
else
	data = fread(fid,[nChannels nSamples],[precision '=>' precision]);
end
data=data';

if isempty(data),
	warning('No data read (trying to read past file end?)');
elseif ~isempty(channels),
	data = data(:,channels);
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function y = Fault_decision_wr2_t1(u)
global st6

% Using model defined by the structure st6 to make decision about type1 wr2 fault

sker=st6.x2sup+(abs(u))'*abs(u)*ones(st6.Nlsup,1)-2*st6.xsup*abs(u);
y=(st6.w)'*exp(-sker./(2*(st6.sigma).^2))+st6.b;

if y>=0
    y=1;
else
    y=0;
end
end


text:
function channels = smoothAngleChannels(channels, skel);

% SMOOTHANGLECHANNELS Try and remove artificial discontinuities associated with angles.

% MOCAP

for i=1:length(skel.tree)
  for j=1:length(skel.tree(i).rotInd)
    col = skel.tree(i).rotInd(j);
    if col
      for k=2:size(channels, 1)
        diff=channels(k, col)-channels(k-1, col);
        if abs(diff+360)<abs(diff)
          channels(k:end, col)=channels(k:end, col)+360;
        elseif abs(diff-360)<abs(diff)
          channels(k:end, col)=channels(k:end, col)-360;
        end
      end
    end
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [Hr,w,theta] = AmpResp(h,FreqPoint,FreqRange,ShiftFreq)

% AmpResp           Calculate Amplitude Response of Symmetric FIR Filter
%
% Arguments:
% h                 Symmetric impulse response:
%                     Type I linear-phase FIR filter, L is odd
%                     Type II lienar-phase FIR fitler, L is even
% FreqPoint         Range of amplitude response from 0 to 2pi or -pi to pi
% FreqRange         FreqRange = 1, frequency range is from 0 to 2pi
%                             = 2, frequency range -pi to pi
% ShiftFreq         Freq shift of amplitude response with specified amount
% Hr                Amplitude (zero-phase) response (column vectors)
%
% by Kong A. Lee, Feb 2005

% Check input argument values

if nargin <4; ShiftFreq = 0; end
if nargin <3; FreqRange = 1; end
if nargin <2; FreqPoint = 4096; end

h = h(:);                       % Make sure h is a column vector
L = length(h);                  % Length of filter

% Check selected frequency range

if FreqRange ==1
    w = (0:FreqPoint-1)*2*pi/FreqPoint;
elseif FreqRange ==2
    FreqPoint = FreqPoint/2;
    w = (-FreqPoint+1:1:FreqPoint-1)*pi/FreqPoint;
end

if mod(L,2)==1                  % L is an odd number

    alpha = (L-1)/2;
	a = [h(alpha + 1); 2*h(alpha:-1:1)];
	n = 0:alpha;
	Hr = cos((w + ShiftFreq)'*n)*a;

end

if mod(L,2)==0                  % L is an even number

    alpha = L/2;
    b = 2*h(alpha:-1:1);
    n = (1:alpha) - 0.5;
    Hr = cos((w + ShiftFreq)'*n)*b;

end

theta = -((L-1)/2)*w;

author: CharlesThaCat

repo: acoustic-interference-cancellation


text:
function fg = dtiXformFibersToMrVista(handles, fgNum, scale)
%
% fg = dtiXformFibersToMrVista(handles, [fiberGroupNum], [scale])
%
% Uses the xformVAnatToAcpc (see dtiXformVanatCompute) to convert the specified
% fiber group (ac-pc coords) to mrVista vAnatomy coords.
%
% HISTORY:
% 2004.05.03 RFD (bob@white.stanford.edu) wrote it.

if(~exist('fgNum','var') | isempty(fgNum))
    fgNum = handles.curFiberGroup;
end
if(~exist('scale','var') | isempty(scale))
    scale = [1 1 1];
end
if(length(scale(:))<3)
    scale = repmat(scale(1),1,3);
end

fg = handles.fiberGroups(fgNum);
xform = diag([scale 1])*inv(handles.xformVAnatToAcpc);
fg = dtiXformFiberCoords(fg, xform);
if(~isempty(fg.seeds))
  fg.seeds = mrAnatXformCoords(xform, fg.seeds);
end
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [h,ax] = scatter3d(v,varargin)
% plot spherical data
%
% Syntax
%   scatter3d(v,data)
%
% Input
%
% See also
% savefigure

% where to plot
if check_option(varargin,'parent')
  ax = get_option(varargin,'parent');
else
  ax = gca;
end

% plot a inner sphere that is not transluent
plotEmptySphere(ax);

% normalize vectors
v = reshape(v,[],1);
v = 1.02 .* v ./ norm(v);

if nargin > 1 && isnumeric(varargin{1})
  data = varargin{1};
  data = reshape(data,length(v),[]);
  varargin{1} = [];
else
  data = {};
end

if v.antipodal   %#ok<BDSCI,BDLGI>
  v = [v;-v];
  data = [data;data];
end

% markerSize
if ~check_option(varargin,{'scatter_resolution','MarkerSize'},'double')
  res = max(v.resolution,0.5*degree);
else
  res = get_option(varargin,'scatter_resolution',1*degree);
end
MarkerSize  = get_option(varargin,'MarkerSize',max(1,min(getMTEXpref('markerSize'),50*res)));

% plot
data = ensurecell(data);
if isempty(data), data = {}; end
h = optiondraw(scatter3(v.x(:),v.y(:),v.z(:),MarkerSize.^2,data{:},'filled','parent',ax),varargin{:});

% add transperency if required
if check_option(varargin,{'MarkerAlpha','MarkerFaceAlpha','MarkerEdgeAlpha'})

  faceAlpha = round(255*get_option(varargin,{'MarkerAlpha','MarkerFaceAlpha'},1));
  edgeAlpha = round(255*get_option(varargin,{'MarkerAlpha','MarkerEdgeAlpha'},1));

  % we have to wait until the markes have been drawn
  mh = [];
  while isempty(mh)
    pause(0.01);
    hh = handle(h);
    mh = [hh.MarkerHandle];
  end

  for j = 1:length(mh)
    mh(j).FaceColorData(4,:) = faceAlpha;
    mh(j).FaceColorType = 'truecoloralpha';

    mh(j).EdgeColorData(4,:) = edgeAlpha;
    mh(j).EdgeColorType = 'truecoloralpha';
  end

end

axis(ax,'equal','vis3d','off');

set(ax,'XDir','rev','YDir','rev',...
'XLim',[-1.02,1.02],'YLim',[-1.02,1.02],'ZLim',[-1.02,1.02]);

hold(ax,'off')

if nargout == 0, clear h;end

end

% since the legend entry for patch object is not nice we draw an
% invisible scatter dot just for legend
%if check_option(varargin,'DisplayName')
%  holdState = get(ax,'nextPlot');
%  set(ax,'nextPlot','add');
  %optiondraw(scatter([],[],'parent',ax,'MarkerFaceColor',mfc,...
  %  'MarkerEdgeColor',mec),varargin{:});%
  %set(ax,'nextPlot',holdState);
%end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function out = deriv(f, xx, varargin)
%DERIV   Evaluate a derivative of a CHEBFUN.
%
%   DERIV(F, X) evaluates the first derivative of the CHEBFUN F at the points in
%   X. If F is a quasimatrix with columns F1, ..., FN, then the result will be
%   [F1(X), ..., FN(X)], the horizontal concatenation of the results of
%   evaluating each column at the points in X.
%
%   DERIV(F, X, M) is the same as above, but returns the values of the Mth
%   derivative of F.
%
%   DERIV(F, X, S) or DERIV(F, X, S, M) where S is one of the strings 'left',
%   'right', '+', or '-', evaluates the left- or right-sided limit as described
%   in chebfun/feval.
%
%   Example 1:
%     u = chebfun(@sin);
%     deriv(u, 1) % u'(1)
%     deriv(u, .5, 2) % u''(.5)
%     deriv(u, 0.1:0.01:0.2) % u'(0.1:0.01:0.2)
%
%   Example 2:
%     x = chebfun('x')
%     u = abs(x);
%     deriv(u, 0, 'left')  % ans = -1
%     deriv(u, 0, 'right') % ans = +1
%
% See also chebfun/diff, chebfun/feval.

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

% The most trivial case:
if ( isempty(f) )
    return
end

% By default, compute first derivative:
m = 1;

% Parse the inputs:
if ( nargin == 3 )
    if ( isnumeric(varargin{1}) ) % DERIV(F, X, M)
        m = varargin{1};
        varargin = {};
    else                          % DERIV(F, X, S)
        m = 1;  % By default, compute first derivative
    end
elseif ( nargin == 4 )            % DERIV(F, X, S, M)
    m = varargin{2};
    varargin{2} = [];
end

if ( m == 0 )
    % Trivial case
    out = feval(f, xx, varargin{:});
else
    out = feval(diff(f, m), xx, varargin{:});
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [fMshift, fProbability, fBic, mProblikelihood] = calc_loglikelihood_dM_Mc(mCat1, mCat2)
% function [fMshift, fProbability, fBic, mProblikelihood] = calc_loglikelihood_dM_Mc(mCat1, mCat2);
% ----------------------------------------------------------------------------------------------
% Calculate log-likelihood estimation of magnitude shift dM between to periods
%
% Incoming variable
% mCat1 : EQ catalog period 1 (Catalog to be modified)
% mCat2 : EQ catalog period 2 (Observed catalog)
%
% Outgoing variable
% fProbability : log-likelihood probabilty
% fMshift      : Magnitude shift with the lowest  max. lieklihood score
% fBic         : Bayesian Information Criterion value
% mProblikelihood : Solution matrix shift and likelihood score
%
% Author: J. Woessner, woessner@seismo.ifg,.ethz.ch
% last update: 28.10.02

% Initialize
vProbability = [];
vMc = [];
vABValue =[];
mFitRes = [];
vX_res = [];
vNCumTmp = [];
mDataPred = [];
mProblikelihood = [];
vfProbability = [];
vMshift = [];
fBinning = 0.1;

% Cut catalog for magnitudes M >= 0
mCatalog = [mCat1; mCat2];
vSel0 = (mCatalog(:,6) >= 0);
mCatalog = mCatalog(vSel0,:);
mCat = mCatalog;

% Determine exact time period
fPeriod1 = max(mCatalog(:,3)) - min(mCatalog(:,3));

% Determine max. and min. magnitude
fMinMag = floor(min(mCatalog(:,6)));
% if fMinMag > 0
%     fMinMag = 0;
% end
fMaxMag = ceil(10 * max(mCatalog(:,6))) / 10;

% Calculate FMD for original catalog
[vFMDorg, vNonCFMDorg] = calc_FMD(mCatalog);
for fMc =fMinMag:0.1:fMaxMag-1
    vFMD = vFMDorg;
    vNonCFMD = vNonCFMDorg;
    vNonCFMD = fliplr(vNonCFMD);
    mCatalog = mCat;
    % Calculate a and b-value for GR-law and distribution vNCum
    [nIndexLo, fMagHi, vSel, vMagnitudes] = fMagToFitBValue(mCatalog, vFMD, fMc);
    if (length(mCatalog(vSel,1)) >= 20)
        [fMeanMag, fBValue, fStdDev, fAValue] =  calc_bmemag(mCatalog(vSel,:), fBinning);
        % Normalize to time period
        vFMD(2,:) = ceil(vFMD(2,:)./fPeriod1);
        vNonCFMD(2,:) = ceil(vNonCFMD(2,:)./fPeriod1);
        % Compute quantity of earthquakes by power law
        vMstep = [fMinMag:0.1:fMaxMag];
        vNCum = 10.^(fAValue-fBValue.*vMstep); % Cumulative number

        % Compute non-cumulative numbers vN
        fNCumTmp = 10^(fAValue-fBValue*(fMaxMag+0.1));
        vNCumTmp  = [vNCum fNCumTmp];
        vN = abs(diff(vNCumTmp));
        %     % Normalizea-value
        fAValue = fAValue./fPeriod1;
        % Normlize vN
        vN = vN./fPeriod1;
        % Data selection
        % Adjust magnitude ranges of vN and vNonCFMD
        vSelM = (vNonCFMD(1,:) >= fMinMag);
        vNonCFMD = vNonCFMD(:,vSelM);
        mData = [vN' vNonCFMD']; % Non cumulative
        vSel = (mData(:,2) >= fMc);
        mDataTest = mData(~vSel,:);
        fNmax = max(mDataTest(:,3));
        % For cutted catalog move lowest magnitude to 0
        mDataTest(:,2) = mDataTest(:,2)-fMinMag;
        if (~isempty(fNmax) & ~isnan(fNmax) & fNmax ~= 0 & length(mDataTest(:,1) > 4))
            mDataTest(:,3) = mDataTest(:,3)/fNmax; % Normalize datavalues for fitting with CDF
            mDataTmp = mData.subset(vSel);
            % Curve fitting: Non cumulative part below Mc
            options = optimset;
            options = optimset('Display','iter','Tolfun',1e-5,'TolX',0.001,'MaxFunEvals', 100000);
            [vX, resnorm, resid, exitflag, output, lambda, jacobian]=lsqcurvefit(@calc_lognormal,[0.1 0.1], mDataTest(:,2), mDataTest(:,3));
            mDataTest(:,1) = logncdf(mDataTest(:,2), vX(1), vX(2))*fNmax;
            if (length(mDataTest(:,2)) > length(vX(1,:)))
                %% Confidence interval
                [vPred,delta] = nlpredci(@calc_lognormal,mDataTest(:,2),vX, resid, jacobian);
            else
                vPred = nan;
                delta = nan;
            end; % END: This section is due for errors produced with datasets less long than amount of parameters in vX
            % For cutted catalog move lowest magnitude back to what it is
            mDataTest(:,2) = mDataTest(:,2)+fMinMag;
            % Results of fitting procedure
            mFitRes = [mFitRes; vX resnorm exitflag];
            %% Set data together
            mDataTest(:,3) = mDataTest(:,3)*fNmax;
            mDataPred = [mDataTest; mDataTmp];
            %%%% Now search dM for best fitting model
            mCat1Mod = mDataPred;
            for fMshift = -0.5:0.1:0.5
                % Apply shift
                mCat1Mod(:,6) = mCat1Mod(:,6)+fMshift;
                % Initialize values
                fMinMag = min([min(mCat1Mod(:,6)) min(mCat2(:,6))]);
                fMaxMag = max([max(mCat1Mod(:,6)) max(mCat2(:,6))]);

                [vPredFMD,vBin1] = hist(mCat1Mod(:,6),0:0.1:fMaxMag);
                [vObsFMD,vBin2] = hist(mCat2(:,6),0:0.1:fMaxMag);
                % Time normalization and round due to Poisson distribution calculation in calc_log10poisspdf
                vPredFMD = ceil(vPredFMD./fPeriod1);
                vObsFMD = ceil(vObsFMD./fPeriod2);

                % Calculate the likelihoods for both models
                vProb_ = calc_log10poisspdf(vObsFMD', vPredFMD');
                % Sum the probabilities
                fProbability = (-1) * sum(vProb_);
                vfProbability = [vfProbability; fProbability];
                vMshift = [vMshift; fMshift];
                mCat1Mod = mCat1;
            end
            vProb_ = calc_log10poisspdf(mDataPred(:,3), ceil(mDataPred(:,1))); % Non-cumulative
            % Sum the probabilities
            fProbability = (-1) * sum(vProb_);
            vProbability = [vProbability; fProbability];
            vMc = [vMc; fMc];
            vABValue = [vABValue; fAValue fBValue];
            % Keep best fitting model
            if (fProbability == min(vProbability))
                fProbMin = fProbability;
                vPredBest = vPred*fNmax;
                vDeltaBest = delta;
                mDat = mDataTest;
                vNBest = vN;
                fMcBest = fMc;
                mDatPredBest = mDataPred;
                vX_res = [vX resnorm exitflag];
                fNmaxBest = fNmax;
            end
        else
            fProbMin = nan;
            fMcBes = nan;
        end; % END of IF fNmax
    end; % END of IF length(mCatalog(vSel,1))

    % Clear variables
    vNCumTmp = [];
    mModelDat = [];
    vNCum = [];
    vSel = [];
    mDataTest = [];
    mDataPred = [];
end; % END of FOR fMc
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Initialize
% mProblikelihood = [];
% vfProbability = [];
% vMshift = [];
% fBinning = 0.1;
%
% % Determine exact time period
% fPeriod1 = max(mCat1(:,3)) - min(mCat1(:,3));
% fPeriod2 = max(mCat2(:,3)) - min(mCat2(:,3));
%
% % Calculate Mc for the two periods
% [fProbMin1, fMcBest1] = calc_McCdf(mCat1, fBinning);
% [fProbMin2, fMcBest2] = calc_McCdf(mCat2, fBinning);
%
% mCat1Mod = mCat1;
%
% for fMshift = -0.5:0.1:0.5
%     % Apply shift
%     mCat1Mod(:,6) = mCat1Mod(:,6)+fMshift;
%     % Initialize values
%     fMinMag = min([min(mCat1Mod(:,6)) min(mCat2(:,6))]);
%     fMaxMag = max([max(mCat1Mod(:,6)) max(mCat2(:,6))]);
%
%     [vPredFMD,vBin1] = hist(mCat1Mod(:,6),0:0.1:fMaxMag);
%     [vObsFMD,vBin2] = hist(mCat2(:,6),0:0.1:fMaxMag);
%     % Time normalization and round due to Poisson distribution calculation in calc_log10poisspdf
%     vPredFMD = ceil(vPredFMD./fPeriod1);
%     vObsFMD = ceil(vObsFMD./fPeriod2);
%
%     % Calculate the likelihoods for both models
%     vProb_ = calc_log10poisspdf(vObsFMD', vPredFMD');
%     % Sum the probabilities
%     fProbability = (-1) * sum(vProb_);
%     vfProbability = [vfProbability; fProbability];
%     vMshift = [vMshift; fMshift];
%     mCat1Mod = mCat1;
% end
%
% %%% Find the minimum loglikelihodd score: if the minimum score is obtained several times, calculate MEAN
% %%% of the magnitude shift
% vdMloglikeli = [vfProbability vMshift];
% vSel = (vdMloglikeli == min(vdMloglikeli(:,1)));
% vdMloglikeli = vdMloglikeli(vSel,:);
% if length(vdMloglikeli(:,1)) > 1
%     fProbability = min(vdMloglikeli(:,1));
%     fMshift = mean(vdMloglikeli(:,2));
% else
%     fProbability = vdMloglikeli(:,1);
%     fMshift = vdMloglikeli(:,2);
% end
% % Solution matrix
% mProblikelihood = [vMshift vfProbability];
%
% %% Bayesian Information Criterion (BIC)
% nDegFree = 1; % Magnitude shift is the degree of freedom
% n_samples = length(mCat1(:,6));
% fBic = 2*fProbability + 2*log(n_samples)*nDegFree;

author: CelsoReyes


text:
function sor_test01 ( w )

%*****************************************************************************80
%
%% SOR_TEST01 tests SOR1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 July 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SOR_TEST01:\n' );

  it_num = 400;
  n = 20;

  x_exact = ( 1 : n )';
  a = dif2 ( n );
  b = a * x_exact;

  x = zeros ( n, 1 );
  x_plot(1:n,1) = x;

  step = 1 : it_num + 1;
  e = nan ( it_num+1, 1 );
  xm = nan ( it_num+1, 1 );

  e(1,1) = ( norm ( a * x - b ) ).^2;

  for it = 1 : it_num

    x_new = sor1 ( n, a, b, x, w );

    e(it+1,1) = ( norm ( a * x_new - b ) ).^2;
    x_plot(1:n,it+1) = x_new;
    xm(it,1) = sum ( ( x_new(:) - x(:) ).^2 ) / n;
%
%  Update the solution
%
    x = x_new;

  end
%
%  Display the error.
%
  figure ( 1 )
  plot ( step, log ( e ), 'm-*' )
  title ( 'Log (Error^2)' )
  xlabel ( 'Step' )
  ylabel ( 'Error' )
  grid
%
%  Display the motion.
%
  figure ( 2 )
  plot ( step, log ( xm ), 'm-*' )
  title ( 'Log (Average generator motion)' )
  xlabel ( 'Step' )
  ylabel ( 'Energy' )
  grid
%
%  Plot the evolution of the locations of the generators.
%
  figure ( 3 )

  y = ( 0 : it_num );
  for k = 1 : n
    plot ( x_plot(k,1:it_num+1), y )
    hold on;
  end
  grid on
  hold off;

  title ( 'Generator evolution.' );
  xlabel ( 'Generator positions' );
  ylabel ( 'Iterations' );

  return
end


text:
function out = double(fa)
% Convert to double precision
% FORMAT double(fa)
% fa - a file_array
%__________________________________________________________________________
% Copyright (C) 2005-2018 Wellcome Trust Centre for Neuroimaging

%
% $Id: double.m 7501 2018-11-30 12:16:58Z guillaume $

out = double(full(fa));

author: spm


text:
function output_noise=addnoise(sig,input_noise,snr)
noise=input_noise;
noise_std_var=sqrt(10^(-snr/10)*(sig(:)'*sig(:))/(noise(:)'*noise(:)));
output_noise=noise_std_var*noise;
end

author: LimingShi


text:
function [ readerobj ] = open_ctiff_reader_noxml( filename, symmetry )
%OPEN_CTIFF_READER Read a TIFF that uses a complex SampleFormat
%
% Unfortunately, we have to write our own code here since the TIFF library
% that MATLAB uses does not support a complex SampleFormat.  Actually, this
% reader works for other TIFFs as well, but for those we can use the
% builtin MATLAB functions.
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

if nargin<2
    symmetry=[0 0 0];
end

%% Get TIFF metadata
tags = read_tiff_tags(filename);

%% Error check
% This function is intentionally only limited capability, purely to
% overcome the fact that the MATLAB TIFF reader doesn't handle complex data
% types.  None of TIFF files with complex data that we have ever seen have
% fallen into one of these cases checked for below.  If we find one, we may
% have to extend the reader to handle that case.  These checks will at
% least warn a user to explain why the reader is not working.
fixed_rowlength = all(tags{1}.StripByteCounts==tags{1}.StripByteCounts(1));
if isequal(class(tags{1}.StripOffsets), class(tags{1}.StripByteCounts))
    % Don't change if possible
    StripOffsets = tags{1}.StripOffsets;
    StripByteCounts = tags{1}.StripByteCounts;
else  % Sometimes come in different integer classes, so cast to common
    StripOffsets = uint64(tags{1}.StripOffsets);
    StripByteCounts = uint64(tags{1}.StripByteCounts);
end
contiguous = all((StripOffsets(1:(end-1)) + ...
    StripByteCounts(1:(end-1))) == StripOffsets(2:end));
if ~fixed_rowlength || ~contiguous
    error('OPEN_CTIFF_READER_NOXML:NONCONTIGUOUS_DATA','This TIFF reader only handles contiguous data.');
end
if isfield(tags{1},'Compression') && tags{1}.Compression~=1
    error('OPEN_CTIFF_READER_NOXML:COMPRESSED_DATA','This TIFF reader does not handle compressed data.');
end
if isfield(tags{1},'FillOrder') && tags{1}.FillOrder~=1
    error('OPEN_CTIFF_READER_NOXML:COMPRESSED_DATA','This TIFF reader does not handle bit-reversed data.');
end
% Reader currently only handle pixel-interleaved ("chunky") data
if isfield(tags{1},'PlanarConfiguration') && tags{1}.PlanarConfiguration~=1
    error('OPEN_CTIFF_READER_NOXML:COMPRESSED_DATA','This TIFF reader does not handle "planar" format (band-interleaved) data.');
end
if isfield(tags{1},'ExtraSamples') && ~isempty(tags{1}.ExtraSamples)
    error('OPEN_CTIFF_READER_NOXML:EXTRA_SAMPLES','This TIFF reader does not handle extra samples per pixel.');
end
if isfield(tags{1},'TileOffsets')
    error('OPEN_CTIFF_READER_NOXML:COMPRESSED_DATA','This TIFF reader does not handle tiled data.');
end

%% Testing pixel reading here
% For now, we just use first image.  Later if we run into TIFF files with
% more than one image, we will have to include all of them and return a
% cell array of reader objects.
image_n = 1;
% Setup SICD metadata stub
meta.native.tiff=tags{image_n};
if symmetry(3)
    meta.ImageData.NumCols=uint32(tags{image_n}.ImageLength);
    meta.ImageData.NumRows=uint32(tags{image_n}.ImageWidth);
else
    meta.ImageData.NumCols=uint32(tags{image_n}.ImageWidth);
    meta.ImageData.NumRows=uint32(tags{image_n}.ImageLength);
end
datasize=double([tags{image_n}.ImageWidth tags{image_n}.ImageLength]);
% Compute parameters necessary for reading pixel data
% Determine TIFF endianness
fid = fopen(filename,'r','b','US-ASCII');
endian = fread(fid,2,'uint8=>char').';
fclose(fid);
switch endian
    case 'II'
        endian = 'l';
    case 'MM'
        endian = 'b';
    % No need for otherwise since this would have already been checked in
    % the read_tiff_tags() call above.
end
if isfield(tags{image_n},'SampleFormat')
    switch tags{image_n}.SampleFormat(1) % Assumes same format for all bands
        case 1
            datatype = 'uint';
            complextype = 0;
        case 2
            datatype = 'int';
            complextype = 0;
        case 3
            datatype = 'float';
            complextype = 0;
        case 5
            datatype = 'int';
            complextype = 1;
        case 6
            datatype = 'float';
            complextype = 1;
        otherwise
            error('Unrecognized datatype.');
    end
else
    datatype = 'uint';
end
bits_per_sample = tags{image_n}.BitsPerSample(1) / (complextype+1); % Assumes same format for all bands
datatype = [datatype num2str(bits_per_sample)];
bands = double(tags{image_n}.SamplesPerPixel);
% Initialize reader object
readerobj=open_generic_reader(filename, datasize,...
    datatype, complextype, double(tags{image_n}.StripOffsets(1)), endian, symmetry, bands);
readerobj.get_meta=@() meta;

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function  outputData = depatchify(patchData, DParam)
% *************************************************************************
% Video Super-Resolution with Convolutional Neural Networks
% patchify
%
% creates image from patches
%
% DParam:     training parameters
%
%
% Version 1.0
%
% Created by:   Armin Kappeler
% Date:         07/21/2015
%
%
% *************************************************************************
USE_PADDING = 1; %make it divisible by "stride" by extending the image

bordercropping = 0;
overlap = DParam.patchSize - DParam.stride;
offset = (DParam.patchSize - DParam.outputPatchSize)/2; % because we only have center part

outputData = zeros([DParam.inputImgSize(2),DParam.inputImgSize(1)]);
mask = zeros([DParam.inputImgSize(2),DParam.inputImgSize(1)]);

if USE_PADDING == 1
    paddingSize = mod(DParam.stride-(DParam.inputImgSize - DParam.patchSize),DParam.stride);
    [ xPos, yPos ] = getPatchPos( [DParam.inputImgSize(2)+paddingSize(2),DParam.inputImgSize(1)+paddingSize(1)],DParam.patchSize,overlap,bordercropping);
else
    [ xPos, yPos ] = getPatchPos( [DParam.inputImgSize(2),DParam.inputImgSize(1)],DParam.patchSize,overlap,bordercropping);
end
xPos = xPos + offset;
yPos = yPos + offset;

if size(patchData,4) ~= length(xPos)
    error('Error the number of patches is not consistent with the image size / patch size / overlap configuration')
end

if DParam.zeropadding==1
    offset = DParam.outputPatchSize-1;
    for i = 1:size(patchData,4)

        offsety = min(offset,DParam.inputImgSize(2)-yPos(i));
        offsetx = min(offset,DParam.inputImgSize(1)-xPos(i));
        patchStartIdx = (DParam.patchSize-DParam.outputPatchSize)/2 + 1;
        %if offset~=offsetx || offsety ~=offset
        %   display('error')
        %end

    %    outputData(yPos(i):yPos(i)+DParam.outputPatchSize-1,xPos(i):xPos(i)+DParam.outputPatchSize-1,:) = patchData(:,:,:,i);
    %    outputData(yPos(i):yPos(i)+DParam.outputPatchSize-1,xPos(i):xPos(i)+DParam.outputPatchSize-1,:) = patchData(13:24,13:24,:,i);
        outputData(yPos(i):yPos(i)+offsety,xPos(i):xPos(i)+offsetx,:) = patchData(patchStartIdx:patchStartIdx+offsety,patchStartIdx:patchStartIdx+offsetx,:,i);
        mask(yPos(i):yPos(i)+offsety,xPos(i):xPos(i)+offsetx) = mask(yPos(i):yPos(i)+offsety,xPos(i):xPos(i)+offsetx) + 1;
    end
else
    for i = 1:size(patchData,4)
        outputData(yPos(i):yPos(i)+DParam.outputPatchSize-1,xPos(i):xPos(i)+DParam.outputPatchSize-1,:) = patchData(:,:,:,i);
        mask(yPos(i):yPos(i)+DParam.outputPatchSize-1,xPos(i):xPos(i)+DParam.outputPatchSize-1) = mask(yPos(i):yPos(i)+DParam.outputPatchSize-1,xPos(i):xPos(i)+DParam.outputPatchSize-1) + 1;
    end
end

mask(mask==0) = 1; %avoid division by zero
outputData = outputData./mask;

end

author: thomas-koehler

repo: SupER


text:
% MLRC Muli-response Linear Regression Combiner
%
%   W = A*(WU*MLRC)
%   W = WT*MLRC(B*WT)
%   D = C*W
%
% INPUT
%   A   Dataset used for training base classifiers as well as combiner
%   B   Dataset used for training combiner of trained base classifiers
%   C   Dataset used for testing (executing) the combiner
%   WU  Set of untrained base classifiers, see STACKED
%   WT  Set of trained base classifiers, see STACKED
%
% OUTPUT
%   W   Trained Muli-response Linear Regression Combiner
%   D   Dataset with prob. products (over base classifiers) per class
%
% DESCRIPTION
% Using dataset A that contains the posterior probabilities of each instance
% belonging to each class predicted by the base classifiers to train a
% multi-response linear regression combiner.
% If the original classification problem has K classes, it is converted
% into K seperate regression problems, where the problem for class c has
% instances with responses equal to 1 when they have label c and zero
% otherwise. Put in another way, this function establish a multi-response
% linear regression model for each class and utilize these models to estimate
% the probability that the instances belong to each class.
% Note that in the model for class c, only the probabilities of class c
% predicted by the set of base classifiers are used.
%
% REFERENCE
% 1. Ting, KM, Witten IH. Issues in stacked generalization, Journal of
% Artificial Intelligent Research, 1999, 10: 271-289.
% 2. Dzeroski S, Zenko B. Is combining classifiers with stacking better
% than selecting the best one? Machine Learning, 2004, 54(3): 255-273.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, STACKED, CLASSC, TESTD, LABELD

% Copyright: Chunxia Zhang, 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

function W = mlrc(A)

        name = 'MLR combiner';

    % If there are no inputs, return an untrained mapping.
    % Handle untrained calls like MLRC([])
    if nargin < 1 || isempty(A)
        W = prmapping(mfilename);
        W = setname(W,name);
        return
    end

    islabtype(A,'crisp');       % allow crisp labels only
    isvaldfile(A,1,2);          % at least one object per class and 2 classes

    A = testdatasize(A,'features'); % test whether they fit
    A = setprior(A,getprior(A));    % avoid many warnings
    [m,k,c] = getsize(A);           % size of training set; (m objects; k features; c classes)
    L = k/c;                        % compute the number of classifiers
    A = setfeatlab(A,repmat([1:c]',L,1)); % reset the feature labels of dataset A such that the first c features correspond to
                                          % the first classifier, the next c features correspond to the second classifier, ect.
    C = zeros(c*L,c);                     % register the coefficients of each model
    options = optimset('ToLX',1e-4);
    for i = 1:c                         % run over all classes
        Res = zeros(m,1);
        Index = find(A.featlab == i);   % find the indices correspond to the jth class
        B = seldat(A,[],Index);         % select the data corresponding to the jth class(m x L matrix)
        I = A.nlab == i;
        Res(I) = 1;
        [x,resnorm,residual,exitflag] = lsqnonneg(B.data,Res,[],options); % compute the nonnegative coefficients
        if exitflag == 0
            resnorm
        end

        for j = 1:L
            C((j-1)*c+i,i) = x(j);
        end
    end

    W = affine(C,[],A,getlablist(A),k,c);
    W = setname(W,name);

return

author: marianux

repo: ecg-kit

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


text:
function x=realCepstrum(x,N)
%%REALCEPSTRUM Compute the real cepstrum. The real cepstrum is the inverse
%              Fourier transform of the real logarithm of the magnitude of
%              the Fourier transform of a real sequence. The cepstrum is
%              often used in signal identification. The real cepstrum plays
%              a role in the design of certain filters.
%
%INPUTS: x An xDimX1 or 1XxDim real sequence.
%        N The length of the discrete Fourier transforms to use in the
%          computation of the cepstrum. If this parameter is omitted or an
%          empty matrix is passed, then the default of N=length(x) is used.
%
%OUTPUTS: x The real cepstrum of the input signal x. This will be length-N.
%
%The computation of the real cepstrum is described in [1]. However, the
%implementation here is more straighforward.
%
%REFERENCES:
%[1] T. F. Quatieri and J. M. Tribolet, "Computation of the real cepstrum
%    and minimum-phase reconstruction," in Programs for Digital Signal
%    Processing, Digital Signal Processing Committee, IEEE Acoustics,
%    Speech, and Signal Processing Society, Ed. New York: IEEE Press, 1979,
%    ch. 7.2.
%
%January 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(N))
    N=length(x);
end

if(~isreal(x))
   error('The input must be real')
end

%The ifft should be real; this just deals with any possible finite
%precision issues.
x=real(ifft(log(abs(fft(x,N))),N));

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 [m, maxnum] = mymode(x)
% computes the mode of x (the most commonly occuring value)
% x must be discrete and numeric
% if several values tie for frequency of occurance, each will be returned
% m - the most common value
% maxnum - the number of times that x==m

[g, gn] = grp2idx(x);
maxnum = 0;
m = [];
for i = 1:length(gn)
    count = sum(g==i);
    if count > maxnum
        maxnum = count;
        m = str2num(gn{i});
        c = 1;
    elseif count == maxnum
        m(c+1) = str2num(gn{i});
        c = c + 1;
    end
end

m = sort(m);

author: zouchuhang

repo: LayoutNet


text:
addpath('./scripts');
%% figures
gam = .8;
T = 25;
noise = .2;
seed = 3;
[y, ~, trueSpikes] = gen_data(gam, noise, T, 1, 0.1, [], 1, seed);

fig2_demo_deconvolveAR1(y, gam, 0.4, false, trueSpikes);

%% video
gam = .8;
T = 130;
noise = .2;
seed = 3;
[y, truth, trueSpikes] = gen_data(gam, noise, T, 1, 0.1, [], 1, seed);
fig2_demo_deconvolveAR1(y, gam, 0.4, true, trueSpikes);

author: zhoupc


text:
function[N,B,E] = isi(data,T,err,Nbins,plt)
% Calculate the inter-spike-interval histogram
%    Usage: [N,B,E] = isi(data,T,err,Nbins,plt)
%
% Input:
% Note that all times have to be consistent.
%
% data   - structure array of spike times  (required)
% T      - time interval of interest (default all)
% err    - 0 for no error bars, 1 for jackknife errors
%
% Nbins  - number of bins in the isi
%
% Output:
%
% N      - count in bins
% B      - bin centres
% E      - errorbar (this is 2 sig deviation
%          calculated using a jackknife over trials)

if nargin < 1; error('I need data!'); end
data=padNaN(data); % create a zero padded data matrix from input structural array
data=data'; % transposes data to get it in a form compatible with Murray's routine
if nargin < 2; T = [min(data(:,1)) max(max(data))]; end
if nargin < 3; err = 0;end
if nargin < 4; Nbins = -1; end
if nargin < 5; plt = 'r'; end

if isempty(T); T = [min(min(data)) max(max(data))]; end
if isempty(err); err = 0;end
if isempty(Nbins); Nbins = -1; end
if isempty(plt); plt = 'r'; end

%  get the number of intervals in each trial and the indices of spike times
%  that are kept

NT = length(data(1,:)); % number of trials
NI=zeros(1,NT);
index(1:NT)=struct('keep',[]);
for n=1:NT
  indx = find(data(:,n) >=  T(1) & data(:,n) <=  T(2) ...
                                 & ~isnan(data(:,n)));
  if isempty(indx)
    NI(n) = 0;
  else
    NI(n) = length(indx)-1;
    index(n).keep=indx;
  end
end

% calculate intervals...

I = zeros(NT,max(NI));
IT = [];
for n=1:NT
  I(n,1:NI(n)) = diff(data(index(n).keep,n));
  IT = [IT I(n,1:NI(n))];
end

Mx = max(IT);
if Nbins == -1
  Nbins = floor(sum(NI)/30);
  Med = median(IT);
  Nbins = max(floor(Nbins*Mx/Med),10);
end

B = linspace(0,Mx,Nbins);

N = zeros(NT,Nbins);
for n=1:NT
  N(n,:) = hist(I(n,1:NI(n)),B);
end

% answer...

if NT > 1;Ns = sum(N)/NT;else Ns = N;end
if ~strcmp(plt,'n')
  bar(B,NT*Ns);
end

% Jackknife iver trials to estimate std...

if NT > 4 && err == 1
  MN = 0;
  SN = 0;
  for n=1:NT
    JK = (NT*Ns - N(n,:))/(NT-1);
    MN = MN + JK;
    SN = SN + JK.^2;
  end
  MN = MN/NT;
  SN = SN/NT;
  E = sqrt((NT-1)*(SN - MN.^2));
  if ~strcmp(plt,'n')
    hold on
    errorbar(B,NT*Ns,NT*2*E,'r-')
    hold off
  end
end
N = NT*Ns;


text:
function [bestLabels, bestEpsilon] = find_best_segmentation(result, sequenceData,  groupCount, epsilon)

% find_best_segmentation.m
%
%   Given a known group count, choose the segmentation in a list of segmentations
%   that best codes the data. The approach uses a modified voting scheme to
%   pick a segmentation that has the correct number of groups. If no such
%   segmentation exists, than the segmentation with the minimum penalized
%   coding length is returned, along with a warning.
%
% Inputs:
%   result              - a structure containing the results of applying
%                           greedy algorithm to a particular motion
%                           sequence for a range of epsilons.
%   sequenceData        - a matrix where each column is a vector obtained
%                           from a tracked feature in the motion sequence.
%   groupCount          - the known number of groups.
%   epsilon             - a list containing the range of epsilons for which these results
%                           were obtained.
%
% Outputs:
%   bestLabels          - a labeling of the motion data that is the "best"
%                           segmentation.
%   bestEpsilon         - a choice of epsilon that will cause ALC to
%                         produce this "best" segmentation.
% Dependencies:
%   distinct_labels, total_coding_length
%
% Mar. '08  Shankar Rao -- srrao@uiuc.edu
%
% Copyright 2008, University of Illinois. All rights reserved.

OUTLIER_GROUP_COUNT = 4;

[sampleCount, trialCount] = size(result.labels);
properGroupings = [];
isSingleGrouping = false(1,trialCount);

for trialIndex = 1:trialCount
    currentLabels = result.labels(:, trialIndex);
    distribution = histc(currentLabels, 1:max(currentLabels));
    if (sum(distribution > OUTLIER_GROUP_COUNT) == groupCount)
        properGroupings = [properGroupings trialIndex];
        % Label all outliers as one group, relabel inliers
        % modified 10-24-07
        outlierLabels = find(distribution <= OUTLIER_GROUP_COUNT);
        inlierLabels = find(distribution > OUTLIER_GROUP_COUNT);
        outlierIndices = find(ismember(currentLabels, outlierLabels));
        currentLabels(outlierIndices) = -1;
        for labelIndex = 1:length(inlierLabels)
            currentLabels(result.labels(:, trialIndex) == inlierLabels(labelIndex)) = labelIndex;
        end
        currentLabels(outlierIndices) = groupCount+1;
        result.labels(:, trialIndex) = currentLabels;
    elseif (max(currentLabels) == 1)
        isSingleGrouping(trialIndex) = true;
    end
end

if isempty(properGroupings)
    warning('No segmentation has the correct number of groups. Returning the one with minimal penalized coding length.');
    [ignore winnerIndex] = min(result.penalty);
    bestLabels = result.labels(:, winnerIndex);
    bestEpsilon = epsilon(winnerIndex);
else
    % ignore any epsilon that results in all of the samples being grouped
    % together.
    nonSingleGroupings = find(~isSingleGrouping);

    % only epsilons that result in the proper number of group can be
    % considered as candidate groupings.
    properGroupingsCount = length(properGroupings);

    [distinctLabels labelIndices] = distinct_labels(result.labels(:, properGroupings));
    candidateCount = size(distinctLabels, 2);
    codingLengths = zeros(candidateCount, length(nonSingleGroupings));

    % however, all non single-group segmentations get a vote
    for trialIndex = 1:length(nonSingleGroupings)
        for candidateIndex = 1:candidateCount
            codingLengths(candidateIndex, trialIndex) = ...
                total_coding_length(sequenceData, distinctLabels(:, candidateIndex), epsilon(nonSingleGroupings(trialIndex)), false);
        end
    end

    [ignore votes] = min(codingLengths, [] , 1);
    voteCounts = histc(votes, 1:candidateCount);
    [ignore winnerIndex] = max(voteCounts);
    bestLabels = distinctLabels(:, winnerIndex);
    winnerLabelIndices = find(labelIndices == winnerIndex);
    bestEpsilon = epsilon(uint32(floor(median(winnerLabelIndices))));
end

author: SuTanTank

repo: VideoStitchingViaShakinessRemoving


text:
function [thr, q, obj, im] = gmthr_seg(im, nobj, nsubs)
% GMTHR_SEG  Segment an image estimating threshold as intersection of two
% Gaussians from Gaussian mixture model
%
% THR = gmthr_seg(IM)
%
%   IM is an input n-dim image. This function assumes that IM contains a
%   darker object over a brighter background.
%
%   THR is a scalar with the estimated threshold value between dark
%   voxels (object) and lighter voxels (background). A Gaussian mixture
%   model is fitted to the image intensities, and the intersection point
%   between the Gaussian maxima is computed. The object in the image is
%   segmented using this intersection value as the segmentation threshold.
%
%   If the object and the background are too similar compared to the number
%   of samples in the image (i.e. the Gaussians intersect outside of the
%   interval between the Gaussian maxima), then this method cannot provide
%   a threshold to separate object and background. In that case, THR is
%   returned as NaN. This is the case, for example, if the image only
%   contains background, or only object voxels.
%
% [THR, Q, OBJ, BW] = gmthr_seg(IM, NOBJ, NSUBS)
%
%   NOBJ is a scalar. Only the largest NOBJ are kept in the segmentation.
%   This last step is useful to remove segmentation noise. By default,
%   NOBJ=1.
%
%   NSUBS is a scalar. For large images, the variance estimate will be too
%   small for the Gaussian mixture fitting function, which will return an
%   error. This problem can be solved doing a random subsampling of the
%   image to estimate the Gaussian mixture model. NSUBS is the subsampling
%   factor. E.g. NSUBS=100 will randomly sample numel(IM)/NSUBS voxels in
%   the image. By default, NSUBS=1 and no subsampling is performed.
%
%   Q is a quality measure of the threshold. Q takes values in [0, 1].
%   Values close to 0 mean that both Gaussians have a lot of overlap, so
%   the threshold between object and background cannot be trusted very
%   much. Values close to 1 mean that both Gaussians are well separated,
%   and the threshold value can be trusted to provide a good segmentation.
%
%   OBJ is the Gaussian mixture object. See help('gmdistribution.fit') for
%   details. The mean and variance of the Gaussians can be extracted as
%   obj.mu and obj.Sigma, respectively.
%
%   BW is an output segmentation mask, where voxels == true correspond to
%   the darker object.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2012 University of Oxford
% Version: 0.5.2
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

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

% defaults
if (nargin < 2 || isempty(nobj))
    nobj = 1;
end
if (nargin < 3 || isempty(nsubs))
    nsubs = 1;
end

% % DEBUG: approximate pdf of whole image
% [ftot, xout] = hist(im(:), 100);
% inc = xout(2) - xout(1);
% ftot = ftot / numel(im) / inc;

% compute gaussian mixture model. Note that we need to randomly subsample
% the image so that the variance of the Gaussians is not too small
% (otherwise, gmdistribution.fit() gives an error)
if (nsubs > 1)
    idx = randi(numel(im), round(numel(im)/nsubs), 1);
    obj = gmdistribution.fit(im(idx), 2);
else
    obj = gmdistribution.fit(im(:), 2);
end

% get mixture of Gaussians parameters
[mutis, idx] = min(obj.mu);
vartis = obj.Sigma(idx);
[mubak, idx] = max(obj.mu);
varbak = obj.Sigma(idx);

% % DEBUG: create Gaussian curves for display purposes
% ftis = normpdf(xout, mutis, sqrt(vartis));
% fbak = normpdf(xout, mubak, sqrt(varbak));

% compute intersection points between two gaussians
thr = intersect_gaussians(mutis, mubak, sqrt(vartis), sqrt(varbak));

% keep the one that is between both mean values
thr = thr(thr > mutis & thr < mubak);

% if there's no intersection point between the maxima, then returning a
% threshold is meaningless, and instead we return NaN. This is the case,
% e.g. if there's only background and no object
if isempty(thr)
    thr = nan;
end

% quality of the clustering measure. Integral under the tissue Gaussian
% in [thr, Inf] and integral under the background Gaussian in [-Inf, thr]:
% the sum represents the Gaussian overlap area. This overlap has a value in
% [0, 1], with 0 for a lot of overlap, and 1 for no overlap. The quality
% measure is then 1-overlap
if (nargout < 2)
    return
end
if (isnan(thr))
    q = nan;
else
    q = 1 - normcdf(2*mutis-thr, mutis, sqrt(vartis))...
        - normcdf(thr, mubak, sqrt(varbak));
end

% % DEBUG: plot histogram curves
% hold off
% plot(xout, ftot)
% hold on
% plot(xout, ftis, 'r')
% plot(xout, fbak, 'g')
% plot([thr, thr], [0 max([ftis(:); fbak(:)])], 'k')
% legend('all', 'tissue', 'background')
% xlabel('intensity')

% no need to waste time segmenting the image if the user doesn't ask for
% the output segmentation
if (nargout < 4)
    return
end

% threshold segmentation
im = (im <= thr);

%% remove segmentation noise
% we could use function bwrmsmallcomp() here, but we don't want having to
% replicate the segmentation data too many times. That could create memory
% problems for very large volumes. Instead, we have copied the code in
% bwrmsmallcomp() directly here:

% get connected components
cc = bwconncomp(im);

% keep only the largest components to remove background noise
%
% note: it's better to clear the whole image, and then add the largest
% components, than trying to delete the smaller components. The latter
% doesn't remove all the noise, for some reason.
len = cellfun(@length, cc.PixelIdxList);
[~, idx] = sort(len, 2, 'descend');
im = zeros(size(im), 'uint8');
if ~isempty(idx)
    idx = cat(1, cc.PixelIdxList{idx(1:nobj)});
    im(idx) = 1;
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function B = rowFirstNonZero(A)
%DESCRIPTION: get first non zero element of each column of a matrix
[~, c] = max( A ~=0, [], 2 );
d = (c-1)*size(A,2) +[1:size(A,1)]';
B=A(d);
end

author: goGPS-Project


text:
function [am,em]=v_lpcrr2am(rr);
%V_LPCRR2AM Convert autocorrelation coefs to ar coef matrix [AM,EM]=(RR)
%AM is a 3-dimensional matrix of size (p+1,p+1,nf) where p is the lpc order
%and nf the number of frames.
%The matrix AM(:,:,*) is upper triangular with 1's on the main diagonal
%and contains the lpc coefficients for all orders from p down to 0.
%
%For lpc order p+1-r, AM(r,r:p+1,*), AM(p+1:-1:r,p+1,*) and EM(*,r) contain
%the lpc coefficients, reflection coefficients and the residual energy respectively.
%
%If A=am(:,:,*), R=toeplitz(rr(*,:)) and E=diag(em(*,:)), then
% A*R*A'=E; inv(R)=A'*(1/E)*A; A*R is lower triangular with the same diagonal as E
%
% This routine is equivalent to: c=chol(inv(toeplitz(rr))); d=diag(c).^-1; em=d.^2; am=diag(d)*c

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

[nf,p1]=size(rr);
p=p1-1;
p2=p1+1;
am=zeros(nf,p1,p1);
em=zeros(nf,p1);
am(:,p1,p1)=1;
em(:,p1)=rr(:,1);
ar=ones(nf,p1);
ar(:,2) = -rr(:,2)./rr(:,1);
e = rr(:,1).*(ar(:,2).^2-1);
for n = 2:p
   q=p2-n;
   em(:,q)=-e;
   am(:,q:p1,q)=ar(:,1:n);
   k = (rr(:,n+1)+sum(rr(:,n:-1:2).*ar(:,2:n),2)) ./ e;
   ar(:,2:n) = ar(:,2:n)+k(:,ones(1,n-1)).*ar(:,n:-1:2);
   ar(:,n+1) = k;
   e = e.*(1-k.^2);
end
em(:,1)=-e;
am(:,:,1)=ar;
am=permute(am,[3 2 1]);

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function pass = example_graze_pair_filt(sicdmeta1, sicdmeta2)
%EXAMPLE_GRAZE_PAIR_FILT An example image pair filter used by 'image_pair_search'
% EXAMPLE_GRAZE_PAIR_FILT(sicdmeta1, sicdmeta) takes two SICD-compatible
% structures and returns a boolean indicating pass/fail.  Any image/SICD
% structure pair which yields a "pass" (boolean true) will be included in
% the list returned by image_pair_search.
%
% Note: No this isn't a useful filter.  At all.  You should never use it.
% We (i.e. you) should, of course, provide a filter specific to your
% particular use.  This is meant to be a stub...
%
% Just FYI...  This filter will return images whos grazing angles are
% within 5 degrees and have lat,lon locations within 0.001.  No.  There's
% no real point to it.
%
% See Also:  example_graze_filt, write_SICD_KML_placemark,
%            write_KML_header, write_KML_footer
%
% Author: Tom Krauss, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

  pass = false;

  if isfield(sicdmeta1,'SCPCOA')  &&  ...
     isfield(sicdmeta1.SCPCOA,'GrazeAng')  &&  ...
     isfield(sicdmeta1, 'GeoData')  &&  ...
     isfield(sicdmeta1.GeoData, 'SCP')  &&  ...
     isfield(sicdmeta1.GeoData.SCP, 'LLH')  &&  ...
     isfield(sicdmeta2,'SCPCOA')  &&  ...
     isfield(sicdmeta2.SCPCOA,'GrazeAng')  &&  ...
     isfield(sicdmeta2, 'GeoData')  &&  ...
     isfield(sicdmeta1.GeoData, 'SCP')  &&  ...
     isfield(sicdmeta2.GeoData.SCP, 'LLH')

    if abs(sicdmeta1.SCPCOA.GrazeAng-sicdmeta2.SCPCOA.GrazeAng) < 5 && ...
       abs(sicdmeta1.GeoData.SCP.LLH.Lat-sicdmeta2.GeoData.SCP.LLH.Lat) < 0.001 && ...
       abs(sicdmeta1.GeoData.SCP.LLH.Lon-sicdmeta2.GeoData.SCP.LLH.Lon) < 0.001
      pass = true;
    end

  end
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////


text:
function linplus_test50 ( )

%*****************************************************************************80
%
%% TEST50 tests R8PO_FA, R8PO_SL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST50\n' );
  fprintf ( 1, '  R8PO_FA factors a positive definite symmetric\n' );
  fprintf ( 1, '    linear system,\n' );
  fprintf ( 1, '  R8PO_SL solves a factored system.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );

  for i = 1 : n
    for j = 1 : n
      a(i,j) = min ( i, j );
    end
  end
%
%  Set the desired solution.
%
  x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
  b = r8po_mxv ( n, a, x );
%
%  Factor the matrix.
%
  [ a_lu, info ] = r8po_fa ( n, a );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Fatal error!\n' );
    fprintf ( 1, '  R8PO_FA declares the matrix is singular!\n' );
    fprintf ( 1, '  The value of INFO is %d\n', info );
    return
  end
%
%  Solve the linear system.
%
  x = r8po_sl ( n, a_lu, b );

  r8vec_print ( n, x, '  Solution:' );
%
%  Set the desired solution.
%
  x(1:n) = 1;
%
%  Compute the corresponding right hand side, using the factored matrix.
%
  b = r8po_ml ( n, a, x );
%
%  Solve the linear system.
%
  x = r8po_sl ( n, a, b );

  r8vec_print ( n, x, '  Solution:' );

  return
end

author: johannesgerer


text:
function varargout = plot(v,varargin)
% plot vectors as two dimensional projections on the sphere
%
% Syntax
%   plot(v)
%   plot(v,value)
%   plot(v,rgb)
%   plot(v,'MarkerSize',10)
%   plot(v,'contourf')
%   plot(v,'contour')
%
% Input
%   v - @vector3d
%   value - values to be displayed
%   rgb -
%
% Options
%  Marker           - 'square', 'triangle', 'o','diamond'
%  MarkerSize       -
%  MarkerFaceColor  -
%  MarkerEdgeColor  -
%
% Flags
%  smooth   - plot point cloud as colored density
%  contourf - plot point cloud as filled contours
%  contour  - plot point cloud as contours
%

% maybe we should add this do all subplots
if check_option(varargin,'add2all')
  mtexFig = gcm;
  if isempty(gcm)
    ax = gca;
  else
    ax = mtexFig.children;
  end
  ax = get_option(varargin,'parent',ax);
  varargin = delete_option(varargin,'parent',1);
  varargin = delete_option(varargin,'add2all');

  h = [];
  for i = 1:length(ax)
    h = [h,plot(v,varargin{:},'parent',ax(i),'doNotDraw')]; %#ok<AGROW>
  end

  if nargout >=1, varargout{1} = h; end
  if nargout >=2, varargout{2} = ax; end

  return
end

% extract plot type
plotTypes = {'contour','contourf','smooth','scatter','text','quiver',...
  'line','plane','circle','surf','pcolor','custom','3d','scatter3d'};
plotType = get_flag(varargin,plotTypes);

% default plot type
if isempty(plotType)
  if isOption(v,'plot') && v.opt.plot

    if ~isempty(varargin) && isnumeric(varargin{1}) && all(size(varargin{1}) == [length(v),3])
      plotType = 'surf';
    else
      plotType = 'smooth';
    end

  else
    plotType = 'scatter';
  end
end
varargin = delete_option(varargin,plotTypes(3:end));

% if data is vector3d type is quiver
if ~isempty(varargin) && isa(varargin{1},'vector3d')
  plotType = 'quiver';
end

if any(strcmpi(plotType,{'smooth','contourf','contour','pcolor'}))
  varargin = ensureData(varargin);
end

% call plotting routine according to type
switch lower(plotType)

  case '3d'
    if v.isOption('plot')
      [varargout{1:nargout}] = v.plot3d(varargin{:});
    else
      [varargout{1:nargout}] = v.scatter3d(varargin{:});
    end

  case 'scatter3d'

    [varargout{1:nargout}] = v.plot3d(varargin{:});

  case 'scatter'

    [varargout{1:nargout}] = v.scatter(varargin{:});

  case 'smooth'

    [varargout{1:nargout}] = v.smooth(varargin{:});

  case 'surf'

    [varargout{1:nargout}] = v.surf(varargin{:});

  case 'contourf'

    [varargout{1:nargout}] = v.contourf(varargin{:});

  case 'contour'

    [varargout{1:nargout}] = v.contour(varargin{:});

  case 'pcolor'

    [varargout{1:nargout}] = v.pcolor(varargin{:});

  case 'quiver'

    [varargout{1:nargout}] = v.quiver(varargin{:});

  case 'line'

    [varargout{1:nargout}] = v.line(varargin{:});

  case 'circle'

    [varargout{1:nargout}] = v.circle(varargin{:});

  case 'plane'

    [varargout{1:nargout}] = v.circle(90*degree,varargin{:});

  case 'text'

    [varargout{1:nargout}] = v.text(varargin{:});

  case 'custom'

    [varargout{1:nargout}] = v.plotCustom(varargin{:});

end
end

function v = ensureData(v)
  if ~isempty(v) && ~(islogical(v{1}) ||  isnumeric(v{1}))
    v = [{[]},v];
  end
end


text:
function emap = snakeMap(f,varargin)
%SNAKEMAP Computes an edge map for use in the snake iterative algorithm.
%   EMAP = SNAKEMAP(F,T,SIG,NSIG,ORDER) computes the edge map, EMAP, of
%	 input image F. If only F is provided in the input, EMAP will equal
%	 the magnitude of the gradient of F (see Chpapter 3) without
%	 thresholding. If T is also provided, EMAP is thresholded so that, at
%	 any point, EMAP = 1 if EMAP > T; it is zero otherwise. If T is the
%	 string 'auto', thresholding is done automatically by function edge.
%	 Otherwise, T is expected to be a number in the range [0,1]. No
%	 smoothing filter is applied. If all inputs are provided, the
%	 magnitude of the gradient is thresholded and filtered using a
%	 Gaussian kernel of size NSIG*SIG-by-NSIG*SIG where SIG is the
%	 standard deviation and NISG is an integer giving the number of
%	 standard deviations we want the kernel size to be. The function
%	 adjusts the kernel size to be the odd size closest to the specified
%	 size. ORDER is a character string with the following possible
%	 values: If order = 'before' then image F is filtered before the map
%	 (magnitude of the gradient) is computed. If order = 'after'
%	 filtering is performed on the edge map after it is computed. If
%	 ORDER = 'both' filtering is performed on the image first, and then
%	 after the edge map is computed. If ORDER = 'none', no filtering is
%	 performed. This is the default. Gradients are computed using the
%	 Sobel kernels.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% PRELIMINARIES.
% Scale image to the floating point intensity range [0,1].
f = im2double(f);
% Set defaults.
thresholding = true;
filtering = true;
if nargin == 1
	% Only one input. No thresholding.
	thresholding = false;
	% No filtering will be done.
	filtering = false;
elseif numel(varargin) == 1
	% Only two inputs.
	T = varargin{1};
	filtering = false;
elseif numel(varargin) == 4
	T = varargin{1};
	sig = varargin{2};
	nsig = varargin{3};
	order = varargin{4};
else
	% Filtering is expected, but not enough inputs were provided.
	error('No enough inputs to define the filtering operation.')
end

% FILTER THE INPUT IMAGE.
% If filtering is called for, form the Gaussian lowpass kernel.
if filtering
	% Form Gaussian kernel using function fspecial.
	kernelSize = floor(nsig*sig);
	% Find the closest odd number of kernelSize.
	if iseven(kernelSize)
      kernelSize = kernelSize + 1;
	end
	% Gaussian kernel.
	w = fspecial('gaussian',kernelSize,sig);
end
% Check to see if filtering is to be done here.
if filtering
	if isequal(order,'before') || isequal(order,'both')
      f = im2double(imfilter(f,w,'conv','symmetric'));
	end
end

% COMPUTE THE EDGE MAP. We use function edge to compute the gradient
% image. Because thresholding is done as part of the function, we have
% to determine at this point if thresholding is to be done and, if so,
% whether 'auto' was specified or not. Either way, the resulting map is
% thinned automatically by function edge.
if thresholding
	if isequal(T,'auto')
      emap = double(edge(f,'sobel'));
   else
      emap = double(edge(f,'sobel',T));
	end
else
   % No thresholding; use edge with a 0 threshold.
	emap = double(edge(f,'sobel',0));
end
if min(emap(:)) == max(emap(:))
   error('Threshold is too low/high. All values of emap are equal')
end

% FILTER EDGE MAP.
% Check to see if filtering was specified.
if filtering
	% If order is 'after', or 'both' then filter the edge map now.
	if isequal(order,'after') || isequal(order,'both')
      emap = imfilter(im2double(intensityScaling(emap)),...
         w,'conv','symmetric');
	end
end

% SCALE THE EDGE MAP TO THE FULL RANGE [0,1].
emap = intensityScaling(emap);

author: dipum

repo: dipum-toolbox

save_path: github-repos/MATLAB/dipum-dipum-toolbox


text:
function [AtA,A] = corrMatrix3D(obj)

nCha = size(obj.kCalib,4);

if(isreal(obj.kCalib))
    A = zeros(prod(obj.calibSize - obj.kernelSize + 1), prod(obj.kernelSize)*nCha,obj.measPara.precision);
else
    A = complex(zeros(prod(obj.calibSize - obj.kernelSize + 1), prod(obj.kernelSize)*nCha,obj.measPara.precision),zeros(prod(obj.calibSize - obj.kernelSize + 1), prod(obj.kernelSize)*nCha,obj.measPara.precision));
end
% A = [];
counter = 1;
for n=1:nCha
    if(isreal(obj.kCalib(:,:,:,n)))
        if(strcmp(obj.measPara.precision,'single'))
            A(:,counter:counter+prod(obj.kernelSize)-1) = im3colRSingle(obj.kCalib(:,:,:,n),obj.kernelSize).';
        else
            A(:,counter:counter+prod(obj.kernelSize)-1) = im3colR(obj.kCalib(:,:,:,n),obj.kernelSize).'; % before: tmp =
        end
    else
        if(strcmp(obj.measPara.precision,'single'))
            A(:,counter:counter+prod(obj.kernelSize)-1) = im3colCSingle(obj.kCalib(:,:,:,n),obj.kernelSize).';
        else
            A(:,counter:counter+prod(obj.kernelSize)-1) = im3colC(obj.kCalib(:,:,:,n),obj.kernelSize).';
        end
    end
    counter = counter + prod(obj.kernelSize);
% 	A = [A, tmp];
end

AtA = A'*A; %%%%%% mtimesx faster better???? %%%%%%%

end

author: thomaskuestner