text:
function [RFcov, figHandle, all_models, weight, data] = rmPlotCoverage(vw, varargin)
% [RFcov, figHandle, all_models, weight, data] = rmPlotCoverage(vw, varargin)
% rmPlotCoverage - calulate the visual field coverage within an ROI
%
%
% Before you run this script, you have to load 'variance explained', 'eccentricity',
% 'polar-angle' and 'prf size' into 'co', 'map', 'ph' and 'amp' fields, respectively
%
% OUTPUT
%  RFcov
% INPUT
%  prf_size:        0 = plot pRF center; 1 = use pRF size
%  fieldRange:      maximum eccentricity to plot (deg)
%  method:          'sum','max', 'clipped average', 'signed profile'
%  newfig:          make a new figure (1) or not (0). (-1 indicates don't plot
%                       anything, just return the coverage map.)
%  nboot:           the number of bootstrapping (0 means no bootstrapping)
%  normalizeRange:  if true, scale z axis to [0 1]
%  smoothSigma:     median smoothing default: 2 nearest values
%  cothresh:        threshold by variance explained in model
%  eccthresh:       2-vector ecc limits (default = [0 1.5*fieldRange])
%  nsamples:        num samples in square grid (default = 128)
%  weight:          any of {'fixed', 'parameter map', 'variance explained'} (default = fixed)
%  weightBeta:      use beta values from rmModel to weight pRFs (default = false)
%  threshByCoh:     if true, threshold by values in coherence map, not variance explained in model
%                       (note: these are often the same, but don't have to be)  %   addcenters
%  addcenters:      1 = superimpose dots for pRF centers; 0 = do not show centers
%  dualVEthresh:    use both pRF model VE and GLM fit VE as threshold
%

% 08/02 KA wrote
% 08/04 KA added bootstrapping
% 08/04 SD various mods
% 09/02 SD large rearrangments
% 09/08 JW allow superposition of pRF centers; various minor debugs
% 02/10 MB added method 'betasum'
% 09/16 RL edited code to allow for the edge case of ROIs that are 1 or 2 voxels large
if notDefined('vw'),       error('View must be defined.'); end
if notDefined('varargin'), varargin{1} = 'dialog'; end

%% default parameters
vfc.prf_size = true;
vfc.fieldRange = min(30, vw.rm.retinotopyParams.analysis.maxRF);
vfc.method = 'max';
vfc.newfig = true;
vfc.nboot = 50;
vfc.normalizeRange = true;
vfc.smoothSigma = true;
vfc.cothresh = viewGet(vw, 'co thresh');
vfc.eccthresh = [0 1.5*vfc.fieldRange];
vfc.nSamples = 128;
vfc.meanThresh = 0;
vfc.weight = 'fixed';
vfc.weightBeta = 0;
vfc.cmap = 'jet';
vfc.clipn = 'fixed';
vfc.threshByCoh = false;
vfc.addCenters = true;
vfc.verbose = prefsVerboseCheck;
vfc.dualVEthresh = 0;

compVolume = false;

%% parse options
if strcmpi(varargin{1},'dialog')
    % get parameters from a dialog
    vfc = rmPlotCoverageDialog(vfc);
else
    for ii = 1:2:length(varargin)
        switch lower(varargin{ii})
            case 'prf_size',        vfc.prf_size        = varargin{ii+1}; % deg
            case 'fieldrange',      vfc.fieldRange      = varargin{ii+1}; % deg
            case 'method',          vfc.method          = varargin{ii+1}; % string
            case 'newfig',          vfc.newfig          = varargin{ii+1}; % boolean
            case 'nboot',           vfc.nboot           = varargin{ii+1}; % integer
            case 'normalizerange',  vfc.normalizeRange  = varargin{ii+1}; % boolean
            case 'smoothsigma',     vfc.smoothSigma     = varargin{ii+1}; % boolean
            case 'cothresh',        vfc.cothresh        = varargin{ii+1};
            case 'eccthresh',       vfc.eccthresh       = varargin{ii+1};
            case 'nsamples',        vfc.nSamples        = varargin{ii+1};
			case 'minmeanmap',      vfc.meanThresh		= varargin{ii+1};
            case 'weight',          vfc.weight          = varargin{ii+1};
            case 'weightbeta',      vfc.weightBeta      = varargin{ii+1};
			case 'cmap',			vfc.cmap			= varargin{ii+1};
            case 'threshbycoh',     vfc.threshByCoh     = varargin{ii+1};
            case 'addcenters',      vfc.addCenters      = varargin{ii+1}; % boolean
            case 'vfc.verbose',     vfc.verbose         = varargin{ii+1}; % boolean
            case 'dualvethresh',    vfc.dualVEthresh    = varargin{ii+1}; % boolean
        end
    end
end

%% load different pRF parameters
try
    rmModel   = viewGet(vw,'rmSelectedModel');
catch %#ok<CTCH>
    error('Need retModel information. Try using rmSelect. ');
end

% Get coordinates for current ROI
roi.coords   = viewGet(vw, 'roiCoords');
roi.indices  = viewGet(vw, 'roiIndices');
roi.name     = viewGet(vw, 'roiName');
curScan      = viewGet(vw, 'curScan');

% Get co and ph (vectors) for the current scan, within the
% current ROI.
vt      = vw.viewType;
co      = rmCoordsGet(vt, rmModel,'varexp',     roi.indices);
sigma1  = rmCoordsGet(vt, rmModel,'sigmamajor', roi.indices);
sigma2  = rmCoordsGet(vt, rmModel,'sigmaminor', roi.indices);
theta   = rmCoordsGet(vt, rmModel,'sigmatheta', roi.indices);
beta    = rmCoordsGet(vt, rmModel,'beta',       roi.indices);
x0      = rmCoordsGet(vt, rmModel,'x0',         roi.indices);
y0      = rmCoordsGet(vt, rmModel,'y0',         roi.indices);
clear rmModel

%%%%%%%%%%%%%%
% y flip note (r.a.s., 10/2009):
% ------------------------------
% I believe the stimulus generation code has had problems (from around
% 2006-2009) where all stimuli were up/down flipped with respect to the pRF
% sampling grid [X, Y]. As a consequence, all models solved in this time
% seem to have a Y flip. The values saved on disk are off; you can
% manually test this with an ROI like V2d which covers a quarterfield.
%
% I'm now trying to fix this issue. All accessor functions seem to have
% implicitly corrected the flip, but in hard to trace ways. (This involves
% a lot of post-hoc calls to functions like 'flipud' which made things
% confusing.) I'm trying to (1) fix the core problems in the stim code; and
% (2) remove the post-hoc corrections to the accessor code.
%
% But for now, there are a lot of models saved on disk, and it will take a
% long time to fix them. So, I'm keeping the y-flip correction, but making
% it explicit here. When the code is fixed and most models saved on disk
% are correct, we can remove this.
% y0 = -y0;

% ok. I think it is time to remove. I suggest putting in a flag to flip the
% y-dimension if requested, but otherwise not to.
if getpref('VISTA', 'verbose')
    warning('Negative y values plotted as lower visual field. This may be incorrect for old pRF models, as old code used to treat neagitve y as upper field.'); %#ok<*WNTAG>
end
%%%%%%%%%%%%%%

% grabbing both (x0, y0) and (pol, ecc) are redundant, and allow for the
% two specifications to get separated (because of the y-flip issue). So,
% re-express ecc and pol using x0 and y0.
[ph ecc] = cart2pol(x0, y0);

% If 'threshByCoh' is set (i.e., true), thresholding will be set based on the
% view struct's coherence field, instead of from the variance explained in
% the model.
if vfc.threshByCoh, co = vw.co{curScan}(roi.indices); end
if ~any(co), co = []; end

% If 'dualVEthresh' is set, thresholding will be based on both the VE of
% the model AND the VE of the GLM fit
if vfc.dualVEthresh == 1
    fprintf('[%s] Using dual VE thresholding \n',mfilename);
    if isempty(vw.dualVE)
        error('Need VE from GLM fit');
    else
        co = (co + vw.dualVE(roi.indices)) ./ 2;
    end
end

% Remove NaNs from subCo and subAmp that may be there if ROI
% includes volume voxels where there is no data.
NaNs = sum(isnan(co));
if NaNs
    fprintf('[%s]:WARNING:ROI includes voxels that have no data. These voxels are being ignored.',mfilename);
    notNaNs = ~isnan(co);
    co      = co(notNaNs);
    ph      = ph(notNaNs);
    ecc     = ecc(notNaNs);
end

% Find voxels which satisfy cothresh and eccthresh.
coIndices = co>vfc.cothresh & ...
    ecc>=vfc.eccthresh(1) & ecc<=vfc.eccthresh(2);
if ~any(coIndices)
   fprintf(1,'[%s]:No values above threshold.\n',mfilename);
   RFcov = zeros(vfc.nSamples);
   figHandle    = [];
   all_models   = [];
   weight       = [];
   data         = [];
   return
end

% also select by the mean map, if that's selected
if vfc.meanThresh > 0
	meanMapFile = fullfile(dataDir(vw), 'meanMap.mat');
	if ~exist(meanMapFile, 'file')
		warning(['A mean map threshold is specified, but no mean map ' ...
				 'exists. This threshold will be ignored for now.'])
	elseif isempty(map{vw.curScan})
		warning(['A mean map threshold is specified, but no mean map ' ...
				 'is computed for the current scan. ' ...
				 'This threshold will be ignored for now.'])
	else
		load(meanMapFile, 'map');
		meanVals = map{vw.curScan};
		if NaNs
			meanVals = meanVals(notNaNs);
		end
		coIndices = coIndices & (meanVals > vfc.meanThresh);
	end
end

% check
if vfc.verbose
    fprintf(1,'[%s]:co-thresh:%.2f.\n',mfilename,vfc.cothresh);
    fprintf(1,'[%s]:ecc-thresh:[%.2f %.2f].\n',mfilename,vfc.eccthresh(1),vfc.eccthresh(2));
    fprintf(1,'[%s]:Number of voxels above thresh in ROI: %d (total=%d).\n',...
        mfilename,sum(coIndices),numel(coIndices));
end

% Pull out co and ph for desired pixels
subCo    = co(coIndices);
subPh    = ph(coIndices);
subEcc   = ecc(coIndices);
subSize1 = single(sigma1(coIndices));
subSize2 = single(sigma2(coIndices));
subTheta = single(theta(coIndices));
subx0    = x0(coIndices);
suby0    = y0(coIndices);

% smooth sigma
if vfc.smoothSigma

    % Cannot do sigma smoothing if the ROI has less than 3 voxels.
    if size(subx0,2) < 3
        error('Cannot perform sigma smoothing when the ROI has less than 3 voxels. ')
    end

    if vfc.smoothSigma == 1
        vfc.smoothSigma = 3; %default
    end
    n = vfc.smoothSigma;

    %check sigma1==sigma2
    if subSize1 == subSize2
        for ii = 1:length(subSize1)
            %compute nearest coords
            dev = sqrt(abs(subx0(ii) - subx0).^2 + abs(suby0(ii) - suby0).^2);
            [dev, ix] = sort(dev); %#ok<*ASGLU>
            subSize1(ii) = median(subSize1(ix(1:n)));
        end
        subSize2 = subSize1;
    else
        for ii = 1:length(subSize1)
            %compute nearest coords
            dev = sqrt(abs(subx0(ii) - subx0).^2 + abs(suby0(ii) - suby0).^2);
            [dev, ix] = sort(dev);
            subSize1(ii) = median(subSize1(ix(1:n)));
            subSize2(ii) = median(subSize2(ix(1:n)));
        end
    end
end

% polar plot
subX = single(subEcc .* cos(subPh));
subY = single(subEcc .* sin(subPh));

% visual field
x = single( linspace(-vfc.fieldRange, vfc.fieldRange, vfc.nSamples) );
[X,Y] = meshgrid(x,x);

% gather this data to make accessible in the plot
if vfc.newfig  > -1   % -1 is a flag that we shouldn't plot the results
	data.figHandle = gcf;
	data.co        = co;
	data.ph        = ph;
	data.subCo     = subCo;
	data.subPh     = subPh;
	data.subEcc    = subEcc;
	data.subx0     = subx0;
	data.suby0     = suby0;
    data.subSize1  = subSize1;
    data.subSize2  = subSize2;
    data.X         = X;
	data.Y         = Y;

end

% For the pRF center plot, use a small constant pRF size
if vfc.prf_size==0
   subSize1 = ones(size(subSize1)) * 0.1;
   subSize2 = ones(size(subSize2)) * 0.1;
   subTheta = zeros(size(subTheta));
end

switch lower(vfc.weight)
    case 'fixed'
        weight = ones(size(subCo));

    case 'parameter map'
        weight = getCurDataROI(vw,'map',curScan,roi.coords);
        weight = weight(coIndices);

    case {'variance explained', 'varexp', 've'}
        weight = subCo;

    otherwise
        error('Unknown weight parameter: %s',vfc.weight);
end

if vfc.weightBeta==1
    weight = weight .* beta(coIndices);
end
weight = single(weight);

%% special case: for the 'density' coverage option, we don't need to
%% do a lot of memory-hungry steps like making all pRFs. So, I've set those
%% computations aside in their own subroutine. (ras)
if isequal( lower(vfc.method), 'density' )
	RFcov = prfCoverageDensityMap(vw, subx0, suby0, subSize1, X, Y);

	all_models = []; % not created for this option
	if vfc.newfig==-1
		figHandle = [];
	else
		[figHandle, data]  = createCoveragePlot(vw, RFcov, vfc, roi, data);
	end

	return
end

%% make all pRFs:
% make in small steps so we don't go into swap space for large ROIs
n = numel(subX);
s = [(1:ceil(n./1000):n-2) n+1];

% For the line above (which assumes that we have at least 3 voxels,
% probably for median smoothing),  s is an empty vector when n < 3
% -- and an empty rfcov is  returned when we try to plot the
% coverage. So modify s accordingly for these edge cases:
% The definition of s is not very intuitive
if n < 3
    s = [1:n+1];
end

all_models = zeros( numel(X), n, 'single' );
fprintf(1,'[%s]:Making %d pRFs:...', mfilename, n);
drawnow;
for n=1:numel(s)-1,
    % make rfs
    rf   = rfGaussian2d(X(:), Y(:),...
						subSize1(s(n):s(n+1)-1), ...
						subSize2(s(n):s(n+1)-1), ...
						subTheta(s(n):s(n+1)-1), ...
						subX(s(n):s(n+1)-1), ...
						subY(s(n):s(n+1)-1));
    all_models(:,s(n):s(n+1)-1) = rf;
end;
clear n s rf pred;
fprintf(1, 'Done.\n');
drawnow;

% Correct volume
if compVolume
    tmp = ones(size(all_models, 1), 1, 'single');

    vol = sigma1(coIndices).^2;
    vol = vol * (2 * pi);

    all_models = all_models ./ (tmp * vol);
end

% For the pRF center plot, put a constant value (1) within each Gaussian
if vfc.prf_size==0
    all_models(all_models>0.1)=1;
end

% weight all models
if isequal( lower(vfc.weight), 'fixed' )
	% if the weights are even, we avoid the redundant, memory-hungry
	% multiplication step that would otherwise be done.
	all_models_weighted = all_models;
else
	tmp = ones(size(all_models, 1), 1, 'single');
	all_models_weighted = all_models .* (tmp * weight);
	clear tmp
end

%% Different ways of combining them:
% 1) bootstrap (yes, no) 2) which statistic (sum, max, etc),
% bootstrap

% If we are only working with 1 voxel, bootstrapping does not do anything.
% We turn it off because the bootstp function does not handle this case well.
if size(subX,2) == 1
    vfc.nboot = 0;
end

if vfc.nboot>0
    if isempty(which('bootstrp'))
        warndlg('Bootstrap requires statistics toolbox');
        RFcov = [];
        return;
    end
    all_models(isnan(all_models))=0;

    switch lower(vfc.method)
        case {'sum','add','avg','average everything', 'average'}
            m = bootstrp(vfc.nboot, @mean, all_models');

        case {'max','profile','maximum profile' 'maximum'}
            m = bootstrp(vfc.nboot, @max, all_models');

        otherwise
            error('Unknown method %s',vfc.method)
    end
    RFcov=median(m,1)';

% no bootstrap
else
    switch lower(vfc.method)

        % coverage = sum(pRF(i)*w(i)) / (sum(pRF(i))
        case {'beta-sum','betasum','weight average'}
            RFcov = sum(all_models_weighted, 2) ./ sum(all_models,2);

        % coverage = sum(pRF(i)*w(i)) / (sum(pRF(i)) + clipping
        case {'clipped beta-sum','clippedbeta','clipped weight average'}
            % set all pRF beyond 2 sigmas to zero
            clipval = exp( -.5 *((2./1).^2));
            all_models(all_models<clipval) = 0;
            n = all_models > 0;

            % recompute all_models_weighted
			tmp = ones( size(all_models,1), 1, 'single' );
            all_models_weighted = all_models .* (tmp*weight);

            % compute weighted clipped sum/average
            sumn = sum(n,2);
            mask = sumn==0;
            sumn(mask) = 1; % prevent dividing by 0
            RFcov = sum(all_models_weighted,2) ./ sum(all_models,2);
            RFcov(mask) = 0;

            %clip to zero if n<clipn
            if isnumeric(vfc.clipn)
                RFcov(sumn<=vfc.clipn) = 0;
            end

        % coverage = sum(pRF(i)*w(i)) / (sum(w(i))
        case {'sum','add','avg','average','prf average'}
            RFcov = sum(all_models_weighted, 2) ./ sum(weight);

        % coverage = sum(pRF(i)*w(i)) / (sum(w(i)) + clipping
        case {'clipped average','clipped','clipped prf average'}
            % set all pRF beyond 2 sigmas to zero
            clipval = exp( -.5 *((2./1).^2));
            all_models(all_models<clipval) = 0;
            n = all_models > 0;

            % recompute all_models_weighted
			tmp = ones( size(all_models,1), 1, 'single' );
            all_models_weighted = all_models .* (tmp*weight);

            % compute weighted clipped mean
            sumn = sum(weight.*n);
            mask = sumn==0;
            sumn(mask) = 1; % prevent dividing by 0
            RFcov = sum(all_models_weighted,2) ./ sumn;
            RFcov(mask) = 0;

            %clip to zero if n<clipn
            if isnumeric(vfc.clipn)
                RFcov(sumn<=vfc.clipn) = 0;
            end

        % coverage = max(pRF(i))
        case {'maximum profile', 'max', 'maximum'}
            RFcov = max(all_models_weighted,[],2);

        case {'signed profile'}
            RFcov  = max(all_models_weighted,[],2);
            covmin = min(all_models_weighted,[],2);
            ii = RFcov<abs(covmin);
            RFcov(ii)=covmin(ii);

        case {'p','probability','weighted statistic corrected for upsampling'}
            RFcov = zeros(vfc.nSamples);

			% I guess this upsample factor assumes your functional data are
			% 2.5 x 2.5 x 3 mm?
            upsamplefactor = 2.5*2.5*3; % sigh.....
            for ii = 1:size(all_models,1)
                s = wstat(all_models(ii,:),weight,upsamplefactor);
                if isfinite(s.tval)
                    RFcov(ii) = 1 - t2p(s.tval,1,s.df);
                end
            end

        otherwise
            error('Unknown method %s',vfc.method)
    end
end

% convert 1D to 2D
RFcov = reshape( RFcov, [1 1] .* sqrt(numel(RFcov)) );

% When no voxels exceed threshold, return nan matrix rather than empty
% matrix
if sum(size(RFcov))==0
    RFcov=nan(nSamples,nSamples);
end

% if the newfig flag is set to -1, just return the image
if vfc.newfig==-1,
    figHandle = [];
else
	[figHandle, data] = createCoveragePlot(vw, RFcov, vfc, roi, data);
end

return
% /--------------------------------------------------------------------/ %

% /--------------------------------------------------------------------/ %
function [figHandle, data] = createCoveragePlot(vw, RFcov, vfc, roi, data)
% plotting subroutine for rmPlotCoverage. Broken off by ras 10/2009.
if vfc.newfig
    figHandle = figure('Color', 'w');
else
	figHandle = selectGraphWin;
end

headerStr = sprintf('Visual field coverage, ROI %s, scan %i', ...
					roi.name, vw.curScan);
set(gcf, 'Name', headerStr);

% normalize the color plots to 1
if vfc.normalizeRange,
	rfMax = max(RFcov(:));
else
	rfMax = 1;
end

img = RFcov ./ rfMax;
mask = makecircle(length(img));
img = img .* mask;
imagesc(data.X(1,:), data.Y(:,1), img);
set(gca, 'YDir', 'normal');
grid on

data.img = img;

colormap(vfc.cmap);
colorbar;

% start plotting
hold on;

% t = 0:.01:2*pi;
%
% % rings every 5 deg
% for n=(1:3)/3*vfc.fieldRange
%     polar(t,ones(size(t))*n,'w');
% end
% plot([0 0],[-vfc.fieldRange vfc.fieldRange],'w')
% plot([-sqrt(vfc.fieldRange^2/2) sqrt(vfc.fieldRange^2/2)],[-sqrt(vfc.fieldRange^2/2) sqrt(vfc.fieldRange^2/2)],'w')
% plot([-vfc.fieldRange vfc.fieldRange],[0 0],'w')
% plot([-sqrt(vfc.fieldRange^2/2) sqrt(vfc.fieldRange^2/2)],[sqrt(vfc.fieldRange^2/2) -sqrt(vfc.fieldRange^2/2)],'w')

% add polar grid on top
p.ringTicks = (1:3)/3*vfc.fieldRange;
p.color = 'w';
polarPlot([], p);

% add pRF centers if requested
if vfc.addCenters,
    inds = data.subEcc < vfc.fieldRange;
    plot(data.subx0(inds), data.suby0(inds), '.', ...
		'Color', [.5 .5 .5], 'MarkerSize', 4);
end

% scale z-axis
if vfc.normalizeRange
	if isequal( lower(vfc.method), 'maximum profile' )
		caxis([.5 1]);
	else
	    caxis([0 1]);
	end
else
    if min(RFcov(:))>=0
        caxis([0 ceil(max(RFcov(:)))]);
    else
        caxis([-1 1] * ceil(max(abs(RFcov(:)))));
    end
end
axis image;   % axis square;
xlim([-vfc.fieldRange vfc.fieldRange])
ylim([-vfc.fieldRange vfc.fieldRange])

title(roi.name, 'FontSize', 24, 'Interpreter', 'none');

% Save the data in gca('UserData')
set(gca, 'UserData', data);

return;
% /------------------------------------------------------------------/ %

% /------------------------------------------------------------------/ %
function RFcov = prfCoverageDensityMap(vw, x0, y0, sigma, X, Y) %#ok<INUSL>
% for each point (x, y) in visual space, this returns
% the proportion of voxels in the ROI for which (x, y) is
% within one standard deviation of the pRF center.
mask = NaN( size(X, 1), size(X, 2), length(x0) );

for v = 1:length(x0)
	% make a binary mask within one sigma of the center
	R = sqrt( (X - x0(v)) .^ 2 + (Y - y0(v)) .^ 2 );
	mask(:,:,v) = ( R < 2*sigma(v) );
end

% average (sum?) across all masks
RFcov = nansum(mask, 3);

return
% /------------------------------------------------------------------/ %

% /------------------------------------------------------------------/ %
function vfc = rmPlotCoverageDialog(vfc)
%% dialog to get parameters for rmPlotCoverage.
dlg(1).fieldName = 'method';
dlg(end).style = 'popup';
dlg(end).list = {'maximum profile' 'average' 'clipped average' ...
				 'density' 'signed profile' 'probability'};
dlg(end).string = 'Method for combining pRFs?';
dlg(end).value = vfc.(dlg(end).fieldName);

dlg(end+1).fieldName = 'weight';
dlg(end).style = 'popup';
dlg(end).list = {'fixed' 'variance explained' 'parameter map'};
dlg(end).string = 'Method for weighting pRFs?';
dlg(end).value = vfc.(dlg(end).fieldName);

dlg(end+1).fieldName = 'fieldRange';
dlg(end).style = 'number';
dlg(end).string = 'Visual Field Range (deg)?';
dlg(end).value = vfc.(dlg(end).fieldName);

dlg(end+1).fieldName = 'nboot';
dlg(end).style = 'number';
dlg(end).string = 'Number of bootstrapping steps?';
dlg(end).value = vfc.(dlg(end).fieldName);

dlg(end+1).fieldName = 'cmap';
dlg(end).style = 'popup';
dlg(end).list = mrvColorMaps;
dlg(end).string = 'If plotting, color map for coverage?';
dlg(end).value = vfc.(dlg(end).fieldName);

dlg(end+1).fieldName = 'normalizeRange';
dlg(end).style = 'checkbox';
dlg(end).list = {};
dlg(end).string = 'Normalize data range to [0 1]';
dlg(end).value = vfc.(dlg(end).fieldName);

dlg(end+1).fieldName = 'smoothSigma';
dlg(end).style = 'checkbox';
dlg(end).list = {};
dlg(end).string = 'Smooth sigma (medianfilter)';
dlg(end).value = vfc.(dlg(end).fieldName);

dlg(end+1).fieldName = 'prf_size';
dlg(end).style = 'checkbox';
dlg(end).list = {};
dlg(end).string = 'Use pRF sizes from model';
dlg(end).value = vfc.(dlg(end).fieldName);

dlg(end+1).fieldName = 'newfig';
dlg(end).style = 'checkbox';
dlg(end).string = 'Show results in new figure';
dlg(end).value = vfc.(dlg(end).fieldName);

dlg(end+1).fieldName = 'addCenters';
dlg(end).style = 'checkbox';
dlg(end).string = 'Add dots to show pRF centers';
dlg(end).value = vfc.(dlg(end).fieldName);

[resp ok] = generalDialog(dlg, mfilename);
if ~ok
	error('User Aborted.')
end
drawnow;

vfc = mergeStructures(vfc, resp);

return


text:
function moment = truncated_normal_ab_moment ( order, mu, s, a, b )

%*****************************************************************************80
%
%% TRUNCATED_NORMAL_AB_MOMENT returns a moment of the truncated Normal PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Phoebus Dhrymes,
%    Moments of Truncated Normal Distributions,
%    May 2005.
%
%  Parameters:
%
%    Input, integer ORDER, the order of the moment.
%    0 <= ORDER.
%
%    Input, real MU, S, the mean and standard deviation of the
%    parent Normal distribution.
%    0 < S.
%
%    Input, real A, B, the lower and upper truncation limits.
%    A < B.
%
%    Output, real MOMENT, the moment of the PDF.
%
  if ( order < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n' );
    fprintf ( 1, '  ORDER < 0.\n' );
    error ( 'TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n' );
  end

  if ( s <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n' );
    fprintf ( 1, '  S <= 0.0.\n' );
    error ( 'TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n' );
  end

  if ( b <= a )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n' );
    fprintf ( 1, '  B <= A.\n' );
    error ( 'TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n' );
  end

  a_h = ( a - mu ) / s;
  a_pdf = normal_01_pdf ( a_h );
  a_cdf = normal_01_cdf ( a_h );

  if ( a_cdf == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n' );
    fprintf ( 1, '  PDF/CDF ratio fails, because A_CDF is too small.\n' );
    fprintf ( 1, '  A_PDF = %g\n', a_pdf );
    fprintf ( 1, '  A_CDF = %g\n', a_cdf );
    error ( 'TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n' );
  end

  b_h = ( b - mu ) / s;
  b_pdf = normal_01_pdf ( b_h );
  b_cdf = normal_01_cdf ( b_h );

  if ( b_cdf == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n' );
    fprintf ( 1, '  PDF/CDF ratio fails, because B_CDF too small.\n' );
    fprintf ( 1, '  B_PDF = %g\n', b_pdf );
    fprintf ( 1, '  B_CDF = %g\n', b_cdf );
    error ( 'TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n' );
  end

  moment = 0.0;
  irm2 = 0.0;
  irm1 = 0.0;

  for r = 0 : order

    if ( r == 0 )
      ir = 1.0;
    elseif ( r == 1 )
      ir = - ( b_pdf - a_pdf ) / ( b_cdf - a_cdf );
    else
      ir = ( r - 1 ) * irm2 ...
        - ( b_h ^ ( r - 1 ) * b_pdf - a_h ^ ( r - 1 ) * a_pdf ) ...
        / ( b_cdf - a_cdf );
    end

    moment = moment + r8_choose ( order, r ) ...
      * mu ^ ( order - r ) ...
      * s ^ r * ir;

    irm2 = irm1;
    irm1 = ir;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Matlab Script: Cum_timeplot.m
% Modified timeplot-script of ZMAP
% ---------------------------------
% Removed buttons "Reset"  and " "Keep as new catalog"
% Author: J. Woessner
%
% This .m file "timeplot" plots the events select by "circle"
% or by other selection button as a cummultive number versus
% time plot in window 2.
% Time of events with a Magnitude greater than minmag will
% be shown on the curve.  Operates on newt2, resets  b  to newt2,
%     newcat is reset to:
%                       - "a" if either "Back" button or "Close" button is         %                          pressed.
%                       - newt2 if "Save as Newcat" button is pressed.
%Last modification 11/95

report_this_filefun(mfilename('fullpath'));

global tmvar                      %for P-Value
global par1 pplot tmp1 tmp2 tmp3 tmp4 difp loopcheck Info_p
global cplot mess til plo2 cum newt2 statime
global magco selt hndl2 wls_button ml_button;;

zmap_message_center.set_info(' ','Plotting cumulative number plot...');

if ~exist('nosort', 'var') ; nosort = 'of'  ; end

if nosort == 'of'
    [s,is] = sort(newt2.Date);
    newt2 = newt2(is(:,1),:) ;
else  % f
    if t3>t2
        l = min(find(newt2.Date > t2))-1;
        newt2 = [newt2(1:l,:) ; newt2(1,:)*nan ; newt2(l+1:newt2.Count,:) ];
    end

end

if ~exist('xt','var')
xt=[]; % time series that will be used
end
if ~exist('as','var')
    as=[]; % z values, maybe? used by the save callback.
end

think
report_this_filefun(mfilename('fullpath'));

% This is the info window text
%
ttlStr='The Cumulative Number Window                  ';
hlpStr1= ...
    ['                                                     '
    ' This window displays the seismicity in the sel-     '
    ' ected area as a cumulative number plot.             '
    ' Options from the Tools menu:                        '
    ' Cuts in magnitude and  depth: Opens input para-     '
    '    meter window                                     '
    ' Decluster the catalog: Will ask for declustering    '
    '     input parameter and decluster the catalog.      '
    ' AS(t): Evaluates significance of seismicity rate    '
    '      changes using the AS(t) function. See the      '
    '      Users Guide for details                        '
    ' LTA(t), Rubberband: dito                            '
    ' Overlay another curve (hold): Allows you to plot    '
    '       one or several more curves in the same plot.  '
    '       select "Overlay..." and then selext a new     '
    '       subset of data in the map window              '
    ' Compare two rates: start a comparison and moddeling '
    '       of two seimicity rates based on the assumption'
    '       of a constant b-value. Will calculate         '
    '       Magnitude Signature. Will ask you for four    '
    '       times.                                        '
    '                                                     '];
hlpStr2= ...
    ['                                                      '
    ' b-value estimation:    just that                     '
    ' p-value plot: Lets you estimate the p-value of an    '
    ' aftershock sequence.                                 '
    ' Save cumulative number cure: Will save the curve in  '
    '        an ASCII file                                 '
    '                                                      '
    ' The "Keep as newcat" button in the lower right corner'
    ' will make the currently selected subset of eartquakes'
    ' in space, magnitude and depth the current one. This  '
    ' will also redraw the Map window!                     '
    '                                                      '
    ' The "Back" button will plot the original cumulative  '
    ' number curve without statistics again.               '
    '                                                      '];

% Find out of figure already exists
%
[existFlag,figNumber]=figure_exists('Cumulative Number',1);
newCumWindowFlag=~existFlag;
cum = figNumber;

% Set up the Cumulative Number window

if newCumWindowFlag
    cum = figure_w_normalized_uicontrolunits( ...
        'Name','Cumulative Number',...
        'NumberTitle','off', ...
        'NextPlot','replace', ...
        'backingstore','on',...
        'Visible','off', ...
        'Position',[ 100 100 (ZmapGlobal.Data.map_len - [100 20]) ]);

    matdraw

    selt='in';
    add_menu_divider();
    options = uimenu('Label','ZTools');

    uimenu(options,'Label','Cuts in time, magnitude and depth', 'Callback','inpu2')
    uimenu(options,'Label','Cut in Time (cursor) ', 'Callback','timesel(4);timeplot;');
    uimenu(options,'Label','Date Ticks in different format', 'Callback','newtimetick');

    uimenu (options,'Label','Decluster the catalog', 'Callback','inpudenew;')
    iwl = ZG.compare_window_yrs*365/par1;
    uimenu(options,'Label','Overlay another curve (hold)', 'Callback','ho2=true; ')
    uimenu(options,'Label','Compare two rates (fit)', 'Callback','dispma2')
    uimenu(options,'Label','Compare two rates ( No fit)', 'Callback','ic=0;dispma3')
    %uimenu(options,'Label','Day/Night split ', 'Callback','daynigt')

    op3D  =   uimenu(options,'Label','Time series ');
    uimenu(op3D,'Label','Time-depth plot ',...
        'Callback',@(~,~)TimeDepthPlotter.plot(newt2));
    uimenu(op3D,'Label','Time magnitude plot ',...
        'Callback',@(~,~)TimeMagnitudePlotter.plot(newt2));

    op4B = uimenu(options,'Label','Rate changes (beta and z-values) ');

    uimenu(op4B, 'Label', 'beta values: LTA(t) function',...
         'Callback', 'sta = ''bet'',newsta')
    uimenu(op4B, 'Label', 'beta values: "Triangle" Plot',...
         'Callback', ';betatriangle')
    uimenu(op4B,'Label','z-values: AS(t)function',...
         'Callback','set(gcf,''Pointer'',''watch'');sta = ''ast'';newsta')
    uimenu(op4B,'Label','z-values: Rubberband function',...
         'Callback','set(gcf,''Pointer'',''watch'');sta = ''rub'';newsta')
    uimenu(op4B,'Label','z-values: LTA(t) function ',...
         'Callback','set(gcf,''Pointer'',''watch'');sta = ''lta'';newsta')

    op4 = uimenu(options,'Label','Mc and b-value estimation');
    uimenu(op4,'Label','automatic', 'Callback','ho=false,selt = ''in'',; bdiff2')
    uimenu(op4,'Label','automatic - overlay existing plot', 'Callback','ho=true,selt = ''in'',;,bdiff2')
 %   uimenu(op4,'Label','manual', 'Callback','bfitnew(newt2)')
    uimenu(op4,'Label','Estimate Mc', 'Callback','mcperc')
    uimenu(op4,'Label','b with depth', 'Callback','bwithde2')
    uimenu(op4,'Label','b with time', 'Callback','selt=''in'';bwithti')
    uimenu(op4,'label','b with magnitude', 'Callback','bwithmag');
    uimenu(op4,'label','Magnitude of completeness (Mc) with time', 'Callback','mcwti');
    uimenu(op4,'label','Create synthetic catalog', 'Callback','synthb');

    pstring=['global freq_field1 freq_field2 freq_field3 freq_field4 freq_field5 tmp1 tmp2 tmp3 tmp4 tmm magn hpndl1 ctiplo mtpl ttcat;ttcat=newt2;'];
    ptstring=[pstring ' cltipval(2);'];
    pmstring=[pstring ' cltipval(1);'];

    op5 = uimenu(options,'Label','p-value estimation');

    %uimenu(op5,'Label','automatic', 'Callback','global hndl1;ttcat =newt2; clpval(3)')
    %The following instruction calls a program for the computation of the parameters in Omori formula, for the catalog of which the cumulative number graph" is
    %displayed (the catalog newt2).
    uimenu(op5,'Label','Completeness in days after mainshock', 'Callback','mcwtidays')
    uimenu(op5,'Label','Define mainshock and estimate p', 'Callback','ho=false;inpu_main')
    uimenu(op5,'Label','Compute p and overlay existing plot', 'Callback','ho=true; pvalcat')
    %In the following instruction the program pvalcat2.m is called. This program computes a map of p in function of the chosen values for the minimum magnitude and
    %initial time.
    uimenu(op5,'Label','p as a function of time and magnitude', 'Callback','pvalcat2')
    %uimenu(op5,'Label','Run aspar3x', 'Callback',' runasp')
    %The previous option is available only under UNIX. The file must contain the mainshock!!
    % uimenu(op5,'Label',' Help/Info on  aspar3x', 'Callback','  do = [ ''web '' hodi ''/help/aspar.htm ;'' ];eval(do) ')
    uimenu(op5,'Label','Cut catalog at mainshock time', 'Callback','l = min(find( newt2.Magnitude == max(newt2.Magnitude) ));newt2 = newt2(l+1:newt2.Count,:);timeplot ')

    op6 = uimenu(options,'Label','Fractal dimension estimation');
    uimenu(op6,'Label','Compute the fractal dimension D', 'Callback',' E = newt2; org = [2]; startfd');
    uimenu(op6,'Label','Compute D for random catalog', 'Callback',' org = [5]; startfd;');
    uimenu(op6,'Label','Compute D with time', 'Callback',' org = [6]; startfd;');
    uimenu(op6,'Label',' Help/Info on  fractal dimension', 'Callback',' showweb(''fractal''); ')

    uimenu(options,'Label','get coordinates with cursor ',...
         'Callback','gi = ginput(1),plot(gi(1),gi(2),''+'');')
    %uimenu(options,'Label','CUSUM ',...
    % 'Callback','cusum(newt2);')
    %uimenu(options,'Label','CUFIT ',...
    % 'Callback','cufit(newt2);')
    uimenu(options,'Label','Cumlative Moment Release ', 'Callback','morel')
    % uimenu(options,'Label','Time to failure  ', 'Callback','savebufe; bufestefan')
    %uimenu(options,'Label','Compute maximum possibe magnitude (Kijko method)', 'Callback','maxmagkijko')

    op7 = uimenu(options,'Label','Stress Tensor Inversion Tools');
    uimenu(op7,'Label','Invert for stress-tensor - Michael''s Method ', 'Callback','doinvers_michael')
    uimenu(op7,'Label','Invert for stress-tensor - Gephart''s Method ', 'Callback','doinversgep_pc')
    uimenu(op7,'Label','Stress tensor with time', 'Callback','stresswtime')
    uimenu(op7,'Label','Stress tensor with depth', 'Callback','stresswdepth')
    uimenu(op7,'Label',' Help/Info on  stress tensor inversions', 'Callback','  showweb(''stress'') ')
    op5C = uimenu(options,'Label','Histograms');

    uimenu(op5C,'Label','Magnitude',...
        'Callback','global histo;hisgra(newt2.Magnitude,stt1);');
    uimenu(op5C,'Label','Depth',...
        'Callback','global histo;hisgra(newt2.Depth,stt2);');
    uimenu(op5C,'Label','Time',...
        'Callback','global histo;hisgra(newt2.Date,''Time '');');
    uimenu(op5C,'Label','Hr of the day',...
        'Callback','global histo;hisgra(newt2.Date.Hour,''Hr '');');

    %uimenu(options,'Label',' Magnitude signature', 'Callback','dispma0')
    uimenu(options,'Label','Save cumulative number curve', 'Callback',{@calSave1, xt, cumu2})

    uimenu(options,'Label','Save cum #  and z value', 'Callback',{@calSave7, xt, cumu2, as})

    %

%     uicontrol('Units','normal','Position',[.0 .0 .1 .05],'String','Reset', 'Callback','nosort = ''of'';newcat = newcat; newt2 = newcat; stri = ['' '']; stri1 = ['' '']; close(cum); timeplot','tooltip','Resets the catalog to the original selection')
%     uicontrol('Units','normal','Position',[.70 .0 .3 .05],'String','Keep as newcat', 'Callback','newcat = newt2;a=newt2;update(mainmap())','tooltip','Plots this subset in the map window')

    ho2=false;

end
%end;    if figure exist

if ho2
    cumu = 0:1:(tdiff*365/par1)+2;
    cumu2 = 0:1:(tdiff*365/par1)-1;
    cumu = cumu * 0;
    cumu2 = cumu2 * 0;
    n = newt2.Count;
    [cumu, xt] = hist(newt2.Date,(t0b:par1/365:teb));
    cumu2 = cumsum(cumu);

    hold on
    axes(ht)
    tiplot2 = plot(newt2.Date,(1:newt2.Count),'r','era','xor');
    set(tiplot2,'LineWidth',2.0)

    ho2=false
    return
end

figure_w_normalized_uicontrolunits(cum)
delete(gca)
delete(gca)
reset(gca)
dele = 'delete(sicum)';er = 'disp('' '')'; eval(dele,er);
cla
hold off
watchon;

set(gca,'visible','off','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','normal',...
    'LineWidth',1.5,...
    'Box','on','SortMethod','childorder')

if isempty(newcat), newcat =a; end

% select big events ( > minmag)
%
l = newt2.Magnitude > minmag;
big = newt2(l,:);
%big=[];
%calculate start -end time of overall catalog
%R
statime=[];
par2=par1;
t0b = min(a.Date);
n = newt2.Count;
teb = max(a.Date);
ttdif=(teb - t0b)*365;
if ttdif>10                 %select bin length respective to time in catalog
    %par1 = ceil(ttdif/300);
elseif ttdif<=10  &&  ttdif>1
    %par1 = 0.1;
elseif ttdif<=1
    %par1 = 0.01;
end

if par1>=1
    tdiff = round((teb - t0b)*365/par1);
    %tdiff = round(teb - t0b);
else
    tdiff = (teb-t0b)*365/par1;
end
% set arrays to zero
%
%if par1>=1
cumu = 0:1:((teb-t0b)*365/par1)+2;
cumu2 = 0:1:((teb-t0b)*365/par1)-1;
%else
%  cumu = 0:par1:tdiff+2*par1;
%  cumu2 =  0:par1:tdiff-1;
%end
% cumu = cumu * 0;
% cumu2 = cumu2 * 0;

%
% calculate cumulative number versus time and bin it
%
n = newt2.Count;
if par1 >=1
    [cumu, xt] = hist(newt2.Date,(t0b:par1/365:teb));
else
    [cumu, xt] = hist((newt2.Date-newt2(1,3)+par1/365)*365,(0:par1:(tdiff+2*par1)));
end
cumu2=cumsum(cumu);
% plot time series
%
%orient tall
set(gcf,'PaperPosition',[0.5 0.5 5.5 8.5])
rect = [0.25,  0.18, 0.60, 0.70];
axes('position',rect)
hold on
%tiplo = plot(xt,cumu2,'ob');
set(gca,'visible','off')
%tiplo2 = plot(xt,cumu2,'b');
%set(tiplo2,'LineWidth',2.5)

%d = datenum(ceil(a.Date)+1900,a.Date.Month,a.Date.Day,a.Date.Hour,a.Date.Minute,a.Date.Minute*0);
%tiplo2 = plot(d,(1:length(d)),'r-.');
%datetick('x',2)

nu = (1:newt2.Count); nu(newt2.Count) = newt2.Count;

tiplot2 = plot([newt2.Date ],nu,'b');
set(tiplot2,'LineWidth',2.0)

% plot end of data
% pl = plot(teb,newt2.Count,'rs');
%set(pl,'LineWidth',1.0,'MarkerSize',6,...
%  'MarkerFaceColor','r','MarkerEdgeColor','g');

% plot big events on curve
%
if par1>=1
    if ~isempty(big)
        %if ceil(big(:,3) -t0b) > 0
        %f = cumu2(ceil((big(:,3) -t0b)*365/par1));
        l = newt2.Magnitude > minmag;
        f = find( l  == 1);
        bigplo = plot(big(:,3),f,'hm');
        set(bigplo,'LineWidth',1.0,'MarkerSize',10,...
            'MarkerFaceColor','y','MarkerEdgeColor','k')
        stri4 = [];
        [le1,le2] = size(big);
        for i = 1:le1
            s = sprintf('  M=%3.1f',big(i,6));
            stri4 = [stri4 ; s];
        end   % for i

        %te1 = text(big(:,3),f,stri4);
        %set(te1,'FontWeight','normal','Color','k','FontSize',8)
        %end

        %option to plot the location of big events in the map
        %
        % figure_w_normalized_uicontrolunits(map)
        % plog = plot(big(:,1),big(:,2),'or','EraseMode','xor');
        %set(plog,'MarkerSize',ms10,'LineWidth',2.0)
        %figure_w_normalized_uicontrolunits(cum)

    end
end %if big

if exist('stri', 'var')
    %v = axis;
    %if par1>=1
    % axis([ v(1) ceil(teb) v(3) v(4)+0.05*v(4)]);
    %end
    %tea = text(v(1)+0.5,v(4)*0.9,stri) ;
    % set(tea,'FontSize',ZmapGlobal.Data.fontsz.s,'Color','k')
else
    strib = [file1];
end %% if stri

strib = [name];

title2(strib,'FontWeight','normal',...
    'FontSize',ZmapGlobal.Data.fontsz.s,...
    'Color','k')

if par1>=1
    xlabel('Time in years ','FontSize',ZmapGlobal.Data.fontsz.s)
else
    statime=newt2(1,3)-par1/365;
    xlabel(['Time in days relative to ',num2str(statime)],'FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
end
ylabel('Cumulative Number ','FontSize',ZmapGlobal.Data.fontsz.s)
ht = gca;
set(gca,'Color',color_bg);

%clear strib stri4 s l f bigplo plog tea v
% Make the figure visible
%
set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.s,...
    'LineWidth',1.0,'TickDir','out','Ticklength',[0.02 0.02],...
    'Box','on')
figure_w_normalized_uicontrolunits(cum);
%sicum = signatur('ZMAP','',[0.65 0.98 .04]);
%set(sicum,'Color','b')
axes(ht);
set(cum,'Visible','on');
watchoff(cum)
watchoff(map)
zmap_message_center.clear_message();
%par1=par2;
done

author: CelsoReyes


text:
function [cfg] = ft_definetrial(cfg)

% FT_DEFINETRIAL defines the trials or segments of data that will be used for further
% processing and analysis, i.e. the pieces of data that will be read in by
% FT_PREPROCESSING. Trials are defined by their begin and end sample in the data file
% and each trial has an offset that defines where the relative t=0 point (usually the
% sample at which the trigger is detected) is for that trial or segment.
%
% Use as
%   [cfg] = ft_definetrial(cfg)
% where the configuration structure should contain
%   cfg.trialdef       = structure with the details of trial definition, see below
%   cfg.trialfun       = string with the function name, see below (default = 'ft_trialfun_general')
%   cfg.representation = 'numeric' or 'table', determines how the trial definition is returned (default is automatic)
% and furthermore
%   cfg.dataset        = string with the filename
% or
%   cfg.headerfile     = string with the filename
%   cfg.datafile       = string with the filename
% and optionally
%   cfg.headerformat   = string, see FT_FILETYPE (default is automatic)
%   cfg.dataformat     = string, see FT_FILETYPE (default is automatic)
%   cfg.eventformat    = string, see FT_FILETYPE (default is automatic)
%
% In general, a call to FT_DEFINETRIAL results in the trial definition "trl" being
% added to the output configuration. The trials are defined on the basis of events or
% triggers by a user-specified MATLAB function that is subsequently referred to as
% the trial function. The user can specify their own custom function tailored to the
% experimental paradigm, or use one of the default trial functions (see below).
%
% Simple trial definitions (for example based on a single trigger) are supported by
% FT_TRIALFUN_GENERAL, which is specified as the default. It supports the following
% options
%   cfg.trialdef.eventtype  = string, or cell-array with strings
%   cfg.trialdef.eventvalue = number, string, or list with numbers or strings
%   cfg.trialdef.prestim    = number, latency in seconds (optional)
%   cfg.trialdef.poststim   = number, latency in seconds (optional)
%
% To read all data from a continuous file in a single or in multiple segments,
% FT_TRIALFUN_GENERAL understands the following options
%    cfg.trialdef.triallength = duration in seconds (can also be 1 or Inf)
%    cfg.trialdef.ntrials     = number of trials (can also be 1 or Inf)
%    cfg.trialdef.overlap     = number between 0 and 1 (exclusive) specifying the fraction of overlap between snippets (0 = no overlap)
%
% To display a list with the events in your data on screen, you can use
% FT_TRIALFUN_SHOW. This is useful for diagnostics; no actual trials will be defined.
%
% To display a graphical user interface dialog that allows you to select events of
% interest, you can use FT_TRIALFUN_GUI.
%
% The trial definition "trl" is an Nx3 matrix or table, where N is the number of
% trials. The first column contains the sample-indices of the start of each trial
% relative to the start of the raw data, the second column contains the sample
% indices of the end of each trial, and the third column contains the offset of the
% trigger with respect to the trial. An offset of 0 means that the first sample of
% the trial corresponds to the trigger. A positive offset indicates that the first
% sample is later than the trigger, and a negative offset indicates that the trial
% begins before the trigger.
%
% Besides the required three columns in the trial definition "trl" that represent
% start, end and offset, it can have contain additional columns . These additional
% columns can be used by a custom trialfun to provide  information about each trial,
% such as trigger codes, response latencies, trial type and response correctness.
% After FT_PREPROCESSING these additional columns of the "trl" matrix will be
% represented in the "trialinfo" field.
%
% If FT_TRIALFUN_GENERAL or FT_TRIALFUN_GUI has been used to generate the "trl"
% matrix or table, the function may return a fourth column that refers to the
% event-code for the corresponding trial. Whether or not this column is returned
% depends on the acquisition system. In general, this fourth column is generated by
% default if the event codes are represented numerically, or as a string starting
% with 'S' or 'R' (for BrainVision data).
%
% If you need to define the segments of interest on the basis of a conditional
% sequence of events (e.g. stimulus trigger followed by a correct response) or on
% basis of some signal feature that needs to be detected in the data, you should
% supply in cfg.trialfun the name of a function that you wrote yourself and that
% FT_DEFINETRIAL will call. The function receives the cfg structure as input and
% should return a NxM matrix in the same format as "trl" as the output. You can add
% extra custom fields to cfg.trialdef to pass to your own trialfun. See below for
% pointers to some examples.
%
% See also FT_PREPROCESSING, FT_READ_HEADER, FT_READ_EVENT, FT_TRIALFUN_GENERAL,
% FT_TRIALFUN_GUI, FT_TRIALFUN_SHOW, FT_TRIALFUN_BIDS, FT_TRIALFUN_EXAMPLE1,
% FT_TRIALFUN_EXAMPLE2

% Undocumented local options:
% cfg.datafile
% cfg.dataset
% cfg.event
% cfg.trl
% cfg.version

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

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

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble provenance

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

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'dataset2files', 'yes');

% this is deprecated as of 26 Sept 2021, it can be removed in a few years
if issubfield(cfg, 'trialdef.eventtype')
  if isequal(cfg.trialdef.eventtype, '?')
    ft_warning('cfg.trialdef.eventtype=''?'' is deprecated, please specify cfg.trialfun=''ft_trialfun_show''');
    cfg.trialfun = 'ft_trialfun_show';
    cfg.trialdef = rmfield(cfg.trialdef, 'eventtype');
  elseif isequal(cfg.trialdef.eventtype, 'gui')
    ft_warning('cfg.trialdef.eventtype=''gui'' is deprecated, please specify cfg.trialfun=''ft_trialfun_gui''');
    cfg.trialfun = 'ft_trialfun_gui';
    cfg.trialdef = rmfield(cfg.trialdef, 'eventtype');
  end
end % default based on eventtype

if issubfield(cfg, 'trialdef.eventvalue')
  if isequal(cfg.trialdef.eventvalue, '?')
    ft_warning('cfg.trialdef.eventvalue=''?'' is deprecated, please specify cfg.trialfun=''ft_trialfun_show''');
    cfg.trialfun = 'ft_trialfun_show';
    cfg.trialdef = rmfield(cfg.trialdef, 'eventvalue');
  elseif isequal(cfg.trialdef.eventvalue, 'gui')
    ft_warning('cfg.trialdef.eventvalue=''gui'' is deprecated, please specify cfg.trialfun=''ft_trialfun_gui''');
    cfg.trialfun = 'ft_trialfun_gui';
    cfg.trialdef = rmfield(cfg.trialdef, 'eventvalue');
  end
end % default based on eventvalue

% set the default options
cfg.trialfun        = ft_getopt(cfg, 'trialfun', 'ft_trialfun_general'); % the default might already be set above
cfg.representation  = ft_getopt(cfg, 'representation'); % this can 'numeric' or 'table'

if isfield(cfg, 'trl')
  % the trial definition is already part of the configuration
  ft_info('retaining existing trial definition\n');
  trl = cfg.trl;
  if isfield(cfg, 'event')
    ft_info('retaining existing event information\n');
    event = cfg.event;
  else
    event = [];
  end

else
  % try to locate the trialfun
  trialfunSpecified = cfg.trialfun;
  cfg.trialfun = ft_getuserfun(cfg.trialfun, 'trialfun');

  if isempty(cfg.trialfun)
    ft_error('the specified trialfun ''%s'' was not found', trialfunSpecified);
  else
    ft_info('evaluating trial function ''%s''\n', func2str(cfg.trialfun));
  end

  % determine the number of outpout arguments of the user-supplied trial function
  try
    % the nargout function in MATLAB 6.5 and older does not work on function handles
    num = nargout(cfg.trialfun);
  catch
    num = 1;
  end
  if num==1
    % the user-defined function only gives back the trial definition
    trl   = feval(cfg.trialfun, cfg);
    event = [];
  else
    % the user-defined function also gives back detailed information about
    % conditions, reaction time or any other information
    [trl, event] = feval(cfg.trialfun, cfg);
  end

end % if trl is present or trialfun specified

if size(trl,1)<1
  if isequal(cfg.trialfun, 'ft_trialfun_show') || isequal(cfg.trialfun, str2func('ft_trialfun_show'))
    % give a gentle message instead of an error
    ft_notice('no trials have been defined yet, see FT_DEFINETRIAL for further help\n');
  else
    % give an error in all other cases
    ft_error('no trials were defined, see FT_DEFINETRIAL for help');
  end
elseif size(trl,2) < 3
  ft_error('trl must have at least 3 columns, see FT_DEFINETRIAL for help');
end

% set trl to requested type
if isempty(cfg.representation)
  % keep the format as is
elseif strcmp(cfg.representation, 'numeric') && istable(trl)
  % convert from table to numeric
  trl = table2array(trl);
elseif strcmp(cfg.representation, 'table') && isnumeric(trl)
  % convert the numeric array to a table
  trl = array2table(trl);
  % the 3 first columns are named begsample, endsample and offset, and the rest get default names
  trl.Properties.VariableNames(1:3) = {'begsample', 'endsample', 'offset'};
end

% add the events to the output configuration
if isstruct(event)
  ft_info('found %d events\n', length(event));
  cfg.event = event;
elseif istable(event)
  ft_info('found %d events\n', size(event, 1));
  cfg.event = event;
end

% add the trial definition to the output configuration
if size(trl,1)>0
  ft_info('created %d trials\n', size(trl,1));
  cfg.trl = trl;
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble provenance
ft_postamble debug

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% written by Jianxiong Xiao http://mit.edu/jxiao/

function svm_prune(alwaysKeep)

global n;
global X;
global Y;
global D;
global SV;

SV = unique([SV; alwaysKeep(:)]);
SV = SV(SV<=n);
n = numel(SV);

X(1:n,:) = X(SV,:);
Y(1:n) = Y(SV,:);
D(1:n) = D(SV);

author: jianxiongxiao

repo: ProfXkit


text:
function test_suite = test_monotonic2

%   Run specific demo, save values and compare the results to the expected.
%   Works for both xUnit Test Framework package by Steve Eddins and for
%   the built-in Unit Testing Framework (as of Matlab version 2013b).
%
%   See also
%     TEST_ALL, DEMO_MONOTONIC2
%
% Copyright (c) 2011-2012 Ville Tolvanen
% Copyright (c) 2014 Tuomas Sivula

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  % Check which package to use
  if exist('initTestSuite', 'file')
    initTestSuite_path = which('initTestSuite');
  else
    initTestSuite_path = '';
  end
  if ~isempty(initTestSuite_path) && ...
     exist([initTestSuite_path(1:end-15) 'runtests'], 'file')
    % xUnit package
    initTestSuite;
  else
    % Built-in package
    % Use all functions except the @setup
    tests = localfunctions;
    tests = tests(~cellfun(@(x)strcmp(func2str(x),'setup'),tests));
    test_suite = functiontests(tests);
  end
end

% -------------
%     Tests
% -------------

function testRunDemo(testCase)
  % Run the correspondin demo and save the values. Note this test has to
  % be run at lest once before the other test may succeed.
  rundemo(getName())
end

function testPredictiveMeanAndVariance_normal(testCase)
  verifyVarsEqual(testCase, getName(), {'Eft','Varft'}, ...
    'RelTolElement', 0.05, 'RelTolRange', 0.01)
end

function testPredictiveMeanAndVariance_monotonic(testCase)
  verifyVarsEqual(testCase, getName(), {'Eftm','Varftm'}, ...
    'RelTolElement', 0.05, 'RelTolRange', 0.01)
end

% ------------------------
%     Helper functions
% ------------------------

function testCase = setup
  % Helper function to suply empty array into variable testCase as an
  % argument for each test function, if using xUnit package. Not to be
  % used with built-in test framework.
  testCase = [];
end

function name = getName
  % Helperfunction that returns the name of the demo, e.g. 'binomial1'.
  name = mfilename;
  name = name(6:end);
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function [] = tet_mesh_boundary(mesh_root)
%%
%  tet_mesh_boundary - a function that reads in a tet_mesh along with neighbor
%                      information and produces lists of boundary nodes and
%                      boundary edges.  The tet_mesh can correspond to either
%                      linear (4 vertices with right-hand orientation) or
%                      quadratic (4 vertices + 6 midside nodes).  In the
%                      latter case, we assume the following ordering:
%                        1  2  3  4 (1+2) (2+3) (1+3) (2+4) (3+4) (1+4).
%
%%

  e_conn = load(strcat(mesh_root,'.elem'));
  tet_adj = load(strcat(mesh_root,'.neighbors'));

  [n_elem,n_dof] = size(e_conn);

  [ boundary_faces ] = define_boundary_faces( e_conn, tet_adj );

  %  preallocate storage for boundary_nodes array
  [ n_faces, npf ] = size( boundary_faces );

  boundary_nodes = zeros( 1,n_faces*npf );
  boundary_nodes_counter = 0;

  for n_f=1:n_faces
    boundary_nodes( boundary_nodes_counter+1:boundary_nodes_counter+npf ) =...
      boundary_faces( n_f, 1:npf );
    boundary_nodes_counter = boundary_nodes_counter+npf;
  end

  [ boundary_nodes ] = unique( boundary_nodes );

  filename = strcat(mesh_root,'.boundary');
  [fid] = fopen(filename,'w');

  print_str = ' %10d\n ';

  fprintf(fid,print_str,boundary_nodes);

  fclose(fid);

end % function  tet_mesh_boundary

function [ boundary_faces ] = define_boundary_faces( e_conn, tet_adj )
%%
%  define_boundary_faces - a function that loops over elements, finds
%                          boundary faces from -1 entries to the
%                          adjacency array and stores information
%--------------------------------------------------------------------------
  [n_elem,n_dof] = size(e_conn);

  %  Estimate the number of boundary faces for preallocation
  %  (these will always be unique)
  boundary_faces = zeros(n_elem,1+n_dof/2);
  boundary_faces_counter = 0;

  if ( n_dof==4 )
    %----------------------------------------------------------------------
    %  Case: Linear Tetrahedron (4 vertices)
    %----------------------------------------------------------------------
    for n_el=1:n_elem
      % A 1-FACE
      if ( tet_adj(n_el,1) < 0 )
        boundary_faces_counter = boundary_faces_counter+1;
        boundary_faces(boundary_faces_counter,:) = ...
         [ e_conn(n_el, 2)  e_conn(n_el, 3)  e_conn(n_el, 4)     ];
      end

      % A 2-FACE
      if ( tet_adj(n_el,2) < 0 )
        boundary_faces_counter = boundary_faces_counter+1;
        boundary_faces(boundary_faces_counter,:) = ...
         [ e_conn(n_el, 3)  e_conn(n_el, 1)  e_conn(n_el, 4)     ];
      end

      % A 3-FACE
      if ( tet_adj(n_el,3) < 0 )
        boundary_faces_counter = boundary_faces_counter+1;
        boundary_faces(boundary_faces_counter,:) = ...
         [ e_conn(n_el, 1)  e_conn(n_el, 2)  e_conn(n_el, 4)     ];
      end

      % A 4-FACE
      if ( tet_adj(n_el,4) < 0 )
        boundary_face_counter = boundary_face_counter+1;
        boundary_face(boundary_face_counter,:) = ...
         [ e_conn(n_el, 1)  e_conn(n_el, 3)  e_conn(n_el, 2)     ];
      end

    end % element loop

    % end linear tetrahedron case

  elseif ( n_dof==10 )
    %----------------------------------------------------------------------
    %  Case: Quadratic Tetrahedron (10 noded)
    %----------------------------------------------------------------------
    for n_el=1:n_elem
      % A 1-FACE
      if ( tet_adj(n_el,1) < 0 )
        boundary_faces_counter = boundary_faces_counter+1;
        boundary_faces(boundary_faces_counter,:) = ...
         [ e_conn(n_el, 2)  e_conn(n_el, 3)  e_conn(n_el, 4)  ...
           e_conn(n_el, 6)  e_conn(n_el, 9)  e_conn(n_el, 8)     ];
      end

      % A 2-FACE
      if ( tet_adj(n_el,2) < 0 )
        boundary_faces_counter = boundary_faces_counter+1;
        boundary_faces(boundary_faces_counter,:) = ...
         [ e_conn(n_el, 3)  e_conn(n_el, 1)  e_conn(n_el, 4)  ...
           e_conn(n_el, 7)  e_conn(n_el,10)  e_conn(n_el, 9)     ];
      end

      % A 3-FACE
      if ( tet_adj(n_el,3) < 0 )
        boundary_faces_counter = boundary_faces_counter+1;
        boundary_faces(boundary_faces_counter,:) = ...
         [ e_conn(n_el, 1)  e_conn(n_el, 2)  e_conn(n_el, 4)  ...
           e_conn(n_el, 5)  e_conn(n_el, 8)  e_conn(n_el,10)     ];
      end

      % A 4-FACE
      if ( tet_adj(n_el,4) < 0 )
        boundary_faces_counter = boundary_faces_counter+1;
        boundary_faces(boundary_faces_counter,:) = ...
         [ e_conn(n_el, 1)  e_conn(n_el, 3)  e_conn(n_el, 2)  ...
           e_conn(n_el, 7)  e_conn(n_el, 6)  e_conn(n_el, 5)     ];
      end

    end % element loop

    % end quadratic tetrahedron case

  else

    error('tet_mesh_boundary: the routine only accounts for 4 noded\n')
    error('of 10 noded tetrahedra\n')

  end

  boundary_faces = boundary_faces(1:boundary_faces_counter,:);

end % function define_boundary_faces


text:
% resize image to 256*256
% and calculate the image mean
mkdir imgs_256
p = dir('./imgs/');
m1 = [];
for i = 3:numel(p)  %remove . & ..
    str = [ './imgs/',p(i).name,'/*.jpg'];
    pp = dir(str);
    if(numel(pp)==0)
        str = [ './imgs/',p(i).name,'/*.bmp'];
        pp = dir(str);
    end
    if(numel(pp)==0)
        str = [ './imgs/',p(i).name,'/*.png'];
        pp = dir(str);
    end
    for j = 1:numel(pp)
        img_str = [ './imgs/',p(i).name,'/',pp(j).name];
        wimg_str = [ './imgs_256/',p(i).name,'/',pp(j).name(1:end-3),'jpg'];
        wdir_str = ['./imgs_256/',p(i).name];
        mkdir(wdir_str);
        im = imresize(imread(img_str),[256,256]);
        imwrite(im,wimg_str);
        m = mean(mean(im,1),2);
        m1 = cat(4,m1,m);
    end
end

mm = mean(m1,4);

author: layumi

repo: Image-Text-Embedding


text:
function prob_test144 ( )

%*****************************************************************************80
%
%% TEST144 tests TRIANGLE_MEAN, TRIANGLE_SAMPLE, TRIANGLE_VARIANCE;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST144\n' );
  fprintf ( 1, '  For the Triangle PDF:\n' );
  fprintf ( 1, '  TRIANGLE_MEAN returns the mean;\n' );
  fprintf ( 1, '  TRIANGLE_SAMPLE samples;\n' );
  fprintf ( 1, '  TRIANGLE_VARIANCE returns the variance;\n' );

  a = 1.0;
  b = 3.0;
  c = 10.0;

  check = triangle_check ( a, b, c );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST144 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =             %14f\n', a );
  fprintf ( 1, '  PDF parameter B =             %14f\n', b );
  fprintf ( 1, '  PDF parameter C =             %14f\n', c );

  mean = triangle_mean ( a, b, c );
  variance = triangle_variance ( a, b, c );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter MEAN =          %14f\n', mean );
  fprintf ( 1, '  PDF parameter VARIANCE =      %14f\n', variance );

  for i = 1 : nsample
    [ x(i), seed ] = triangle_sample ( a, b, c, seed );
  end

  mean = r8vec_mean ( nsample, x );
  variance = r8vec_variance ( nsample, x );
  xmax = max ( x(1:nsample) );
  xmin = min ( x(1:nsample) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %6d\n', nsample );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %14f\n', xmax );
  fprintf ( 1, '  Sample minimum =  %14f\n', xmin );

  return
end


text:
function [con_vec] = spm_mlm_makecon (mlm,w)
% Make contrast to test if the subset of coefficients indexed by w = 0 ?
% FORMAT [con_vec] = spm_mlm_makecon (mlm,w)
%
% mlm           MLM data structure containing
%               [p x d] matrix of regression coefficients mlm.wmean
% w             [p x d] matrix of comprising 1's and 0's with
%               1s selecting the coefficients of interest
%
% con_vec       Vectorised contrast matrix that can be passed
%               to spm_mlm_posthoc.m
%
%___________________________________________________________________________
% Copyright (C) 2011 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_mlm_makecon.m 4651 2012-02-09 16:03:39Z will $

[p,d]=size(mlm.wmean);

con_vec=[];
for i=1:p,
    for j=1:d,
        if w(i,j)==1
            con=zeros(p,d);
            con(i,j)=1;
            con_vec=[con_vec;con(:)'];
        end
    end
end

author: spm


text:
function test_failed=test_thresh
%TEST_THRESH  Compare sparse and full thesholding

test_failed=0;
disp(' ===============  TEST_THRESH ================');
global LTFAT_TEST_TYPE;
if ~strcmpi(LTFAT_TEST_TYPE,'double')
   disp(sprintf('Skipping. Cannot work with sparse matrices of type %s.',LTFAT_TEST_TYPE));
   return;
end

lambda=0.1;

ttypes={'hard','soft','wiener'};

for ii=1:2

  if ii==1
    g=tester_rand(3,4);
    field='REAL  ';
    g(2,2)=lambda;
  else
    g=tester_crand(3,4);
    field='CMPLX ';
    g(2,2)=lambda;
  end;

  for jj=1:3
    ttype=ttypes{jj};

    [xo_full, Nfull] = thresh(g,lambda,ttype,'full');
    [xo_sparse, Nsp] = thresh(g,lambda,ttype,'sparse');

    res = xo_full-xo_sparse;
    res = norm(res(:));

    res2 = Nfull-Nsp;

    [test_failed,fail]=ltfatdiditfail(res,test_failed);
    s=sprintf(['THRESH   %s %s %0.5g %s'],field,ttype,res,fail);
    disp(s);

    [test_failed,fail]=ltfatdiditfail(res2,test_failed);
    s=sprintf(['THRESH N %s %s %0.5g %s'],field,ttype,res2,fail);
    disp(s);

    % Extend lambda to:
    % a) Vector
    lambdavec = lambda*ones(numel(g),1);

    [xo_full2, Nfull2] = thresh(g,lambdavec,ttype,'full');
    [xo_sparse2, Nsp2] = thresh(g,lambdavec,ttype,'sparse');

    res_full = xo_full2-xo_full;
    res_sparse = xo_sparse2-xo_sparse;
    res_full = norm(res_full(:));
    res_sparse = norm(res_sparse(:));

    res_nfull = Nfull2-Nfull;
    res_nsparse = Nfull2-Nfull;

    [test_failed,fail]=ltfatdiditfail(res,test_failed);
    s=sprintf(['THRESH VEC FULL  %s %s %0.5g %s'],field,ttype,res_full,fail);
    disp(s);

    [test_failed,fail]=ltfatdiditfail(res,test_failed);
    s=sprintf(['THRESH VEC SPARSE  %s %s %0.5g %s'],field,ttype,res_sparse,fail);
    disp(s);

    [test_failed,fail]=ltfatdiditfail(res2,test_failed);
    s=sprintf(['THRESH VEC N FULL %s %s %0.5g %s'],field,ttype,res_nfull,fail);
    disp(s);

    [test_failed,fail]=ltfatdiditfail(res2,test_failed);
    s=sprintf(['THRESH VEC N SPARSE %s %s %0.5g %s'],field,ttype,res_nsparse,fail);
    disp(s);

    % b) Same shape as g
    lambdamat = lambda*ones(size(g));

    [xo_full3, Nfull3] = thresh(g,lambdamat,ttype,'full');
    [xo_sparse3, Nsp3] = thresh(g,lambdamat,ttype,'sparse');

    res_full = xo_full3-xo_full;
    res_sparse = xo_sparse3-xo_sparse;
    res_full = norm(res_full(:));
    res_sparse = norm(res_sparse(:));

    res_nfull = Nfull3-Nfull;
    res_nsparse = Nfull3-Nfull;

    [test_failed,fail]=ltfatdiditfail(res,test_failed);
    s=sprintf(['THRESH MAT FULL  %s %s %0.5g %s'],field,ttype,res_full,fail);
    disp(s);

    [test_failed,fail]=ltfatdiditfail(res,test_failed);
    s=sprintf(['THRESH MAT SPARSE %s %s %0.5g %s'],field,ttype,res_sparse,fail);
    disp(s);

    [test_failed,fail]=ltfatdiditfail(res2,test_failed);
    s=sprintf(['THRESH MAT N FULL %s %s %0.5g %s'],field,ttype,res_nfull,fail);
    disp(s);

     [test_failed,fail]=ltfatdiditfail(res2,test_failed);
    s=sprintf(['THRESH MAT N SPARSE %s %s %0.5g %s'],field,ttype,res_nsparse,fail);
    disp(s);

  end;

end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [estimate] = ft_inverse_harmony(sourcemodel, sens, headmodel, dat, varargin)

% FT_INVERSE_HARMONY computes a linear estimate of the current in a distributed
% source model using a mesh harmonic based low-pass filter.
%
% Use as
%   [estimate] = ft_inverse_harmony(sourcemodel, sens, headmodel, dat, ...)
% where
%   sourcemodel is the input source model, see FT_PREPARE_SOURCEMODEL
%   sens        is the gradiometer or electrode definition, see FT_DATATYPE_SENS
%   headmodel   is the volume conductor definition, see FT_PREPARE_HEADMODEL
%   dat         is the data matrix with the ERP or ERF
% and
%   estimate    contains the estimated source parameters
%
% Additional input arguments should be specified as key-value pairs and can include
%   'noisecov'             = Nchan x Nchan matrix with noise covariance
%   'noiselambda'          = scalar value, regularisation parameter for the noise covariance matrix (default=0)
%   'filter_order'         = scalar, order of the mesh Butterwirth filter
%   'filter_bs'            = scalar, stop-band of the mesh Butterworth filter
%   'number_harmonics'     = Integer, number of mesh harmonics used (can be empty, the default will then be identity)
%   'lambda'               = scalar, regularisation parameter (can be empty, it will then be estimated from snr)
%   'snr'                  = scalar, signal to noise ratio
%   'scalesourcecov'       = 'no' or 'yes', scale the source covariance matrix R such that trace(leadfield*R*leadfield')/trace(C)=1
%   'connected_components' = number of connected components of the source mesh (1 or 2)
%   'prewhiten'            = 'no' or 'yes', prewhiten the leadfield matrix with the noise covariance matrix C
%
% These options influence the forward computation of the leadfield
%   'reducerank'      = 'no' or number  (default = 3 for EEG, 2 for MEG)
%   'backproject'     = 'yes' or 'no', in the case of a rank reduction this parameter determines whether the result will be backprojected onto the original subspace (default = 'yes')
%   'normalize'       = 'no', 'yes' or 'column' (default = 'no')
%   'normalizeparam'  = parameter for depth normalization (default = 0.5)
%   'weight'          = number or Nx1 vector, weight for each dipole position to compensate for the size of the corresponding patch (default = 1)
%
% This implements
% - Petrov Y (2012) Harmony: EEG/MEG Linear Inverse Source Reconstruction in the
%   Anatomical Basis of Spherical Harmonics. PLoS ONE 7(10): e44439.
%   doi:10.1371/journal.pone.0044439
%
% See also FT_SOURCEANALYSIS, FT_PREPARE_HEADMODEL, FT_PREPARE_SOURCEMODEL

% TODO implement the following options
% - keepleadfield

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

if mod(nargin-4,2)
  % the first 4 arguments are fixed, the other arguments should come in pairs
  ft_error('invalid number of optional arguments');
end

% get the optional input arguments, or use defaults
noisecov       = ft_getopt(varargin, 'noisecov');
lambda         = ft_getopt(varargin, 'lambda');  % can be empty, it will then be estimated based on SNR
noiselambda    = ft_getopt(varargin, 'noiselambda', []);
snr            = ft_getopt(varargin, 'snr');     % is used to estimate lambda if lambda is not specified
keepfilter     = istrue(ft_getopt(varargin, 'keepfilter', false));
dowhiten       = istrue(ft_getopt(varargin, 'prewhiten',  false));
doscale        = istrue(ft_getopt(varargin, 'scalesourcecov', false));

% construct the low-level options for the leadfield computation as key-value pairs, these are passed to FT_COMPUTE_LEADFIELD
leadfieldopt = {};
leadfieldopt = ft_setopt(leadfieldopt, 'reducerank',     ft_getopt(varargin, 'reducerank'));
leadfieldopt = ft_setopt(leadfieldopt, 'backproject',    ft_getopt(varargin, 'backproject'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalize',      ft_getopt(varargin, 'normalize'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalizeparam', ft_getopt(varargin, 'normalizeparam'));
leadfieldopt = ft_setopt(leadfieldopt, 'weight',         ft_getopt(varargin, 'weight'));

% flags to avoid calling isfield repeatedly in the loop over grid positions (saves a lot of time)
hasmom        = isfield(sourcemodel, 'mom');
hasleadfield  = isfield(sourcemodel, 'leadfield');
hasfilter     = isfield(sourcemodel, 'filter');

if isempty(lambda) && isempty(snr) && ~hasfilter
  ft_error('either lambda or snr should be specified');
elseif ~isempty(lambda) && ~isempty(snr)
  ft_error('either lambda or snr should be specified, not both');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find the dipole positions that are inside/outside the brain
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isfield(sourcemodel, 'inside')
  if hasfilter
    sourcemodel.inside = ~cellfun(@isempty, sourcemodel.filter);
  elseif hasleadfield
    sourcemodel.inside = ~cellfun(@isempty, sourcemodel.leadfield);
  else
    sourcemodel.inside = ft_inside_headmodel(sourcemodel.pos, headmodel);
  end
end

% convert to logical representation
sourcemodel = fixinside(sourcemodel);

% keep the original details on inside and outside positions
originside = sourcemodel.inside;
origpos    = sourcemodel.pos;

% select only the dipole positions inside the brain for scanning
sourcemodel.pos    = sourcemodel.pos(originside,:);
sourcemodel.inside = true(size(sourcemodel.pos,1),1);

if hasmom
  sourcemodel.mom = sourcemodel.mom(:,originside);
end

if hasfilter
  ft_info('using precomputed filters\n');
  sourcemodel.filter = sourcemodel.filter(originside);
elseif hasleadfield
  ft_info('using precomputed leadfields\n');
  sourcemodel.leadfield = sourcemodel.leadfield(originside);
else
  ft_info('computing forward model on the fly\n');
  if hasmom
    for i=size(sourcemodel.pos,1)
      % compute the leadfield for a fixed dipole orientation
      sourcemodel.leadfield{i} = ft_compute_leadfield(sourcemodel.pos(i,:), sens, headmodel, leadfieldopt{:}) * sourcemodel.mom(:,i);
    end
  else
    for i=1:size(sourcemodel.pos,1)
      % compute the leadfield
      sourcemodel.leadfield{i} = ft_compute_leadfield(sourcemodel.pos(i,:), sens, headmodel, leadfieldopt{:});
    end
  end
end

% Harmony parameters
num_dip      = size(sourcemodel.leadfield{1}, 2);
num_pos      = size(sourcemodel.leadfield, 2);
num_comp     = ft_getopt(varargin, 'connected_components', 1);
num_harm     = ft_getopt(varargin, 'number_harmonics',     5);
filt_ord     = ft_getopt(varargin, 'filter_order',         5);
lc           = ft_getopt(varargin, 'filter_bs',            5);

% Create the leadfield matrix

% count the number of channels and leadfield components
Nchan   = size(sourcemodel.leadfield{1},1);
Nsource = 0;
for i=1:size(sourcemodel.pos,1)
  Nsource = Nsource + size(sourcemodel.leadfield{i}, 2);
end

% concatenate the leadfield components of all sources into one large matrix
lf = zeros(Nchan, Nsource);
n = 1;
for i=1:size(sourcemodel.pos,1)
  cbeg = n;
  cend = n + num_dip - 1;
  lf(:,cbeg:cend) = sourcemodel.leadfield{i};
  n = n + size(sourcemodel.leadfield{i}, 2);
end

% Create brain harmonics matrix

% Mesh harmonics
[dum,H,d] = mesh_spectrum(sourcemodel,num_harm,num_comp);

% Elimination of pathological zeroth frequencies and concatenating
% harmonics numbers
l = [];
for n = 1:num_comp
  d{n}(1) = 0;
  H{n}(:,1) = 1/sqrt(numel(H{n}(:,1)));
  d_cat = [];
  for m = 1:num_dip
    d_cat = [d_cat; d{n}];
  end
  l = [l; -d_cat];
end

% Total harmonics matrix
if num_comp == 1
  H_blk = H{1};
elseif num_comp == 2
  H_blk = blkdiag(H{1},H{2});
else
  ft_error('this function only supports meshes with 1 or 2 connected components')
end

if num_dip == 1
  H_tot = H_blk;
elseif num_dip == 2
  H_tot = blkdiag(H_blk,H_blk);
elseif num_dip == 3
  H_tot = blkdiag(H_blk,H_blk,H_blk);
else
  ft_error('this function only supports source model with 1, 2 or 3 dimensional source values')
end

% Create harmonic leadfield
lf_h = zeros(size(lf,1),num_harm,num_comp,num_dip); %leadifield in the harmonic domain

for k = 1:num_comp
  for j = 1:num_dip
    ind = (j-1) + 1:num_dip:num_dip*num_pos;
    lf_h(:,:,k,j) = lf(:,ind)*H{k};
  end
end

lf_h = reshape(lf_h,size(lf,1),num_comp*num_harm*num_dip);

% Create harmonic source covariance matrix
q = 1./sqrt(1 + (l/lc).^(2*filt_ord));
sourcecov = diag(q);

% Compute the Harmony spatial filters

% count the number of channels and leadfield components
Nchan   = size(sourcemodel.leadfield{1},1);
Nsource = 0;
for i=1:size(sourcemodel.pos,1)
  Nsource = Nsource + size(sourcemodel.leadfield{i}, 2);
end

% Take the rieal part of the noise cross-spectral density matrix
if isreal(noisecov)
  ft_info('taking the real part of the noise cross-spectral density matrix\n');
  noisecov = real(noisecov);
end

% compute the inverse of the forward model, this is where prior information
% on source and noise covariance would be useful
if isempty(noisecov)
  % use an unregularised minimum norm solution, i.e. using the Moore-Penrose pseudoinverse
  ft_warning('computing a unregularised minimum norm solution. This typically does not work due to numerical accuracy problems');
  w = pinv(lf_h);
elseif ~isempty(noisecov)
  ft_info('using the noise covariance for regularisation\n');
  % the noise covariance has been given and can be used to regularise the solution
  if isempty(sourcecov)
    sourcecov = speye(Nsource);
  end
  % rename some variables for consistency with the publications
  A = lf_h;
  R = sourcecov;
  C = noisecov;

  if dowhiten
    ft_info('prewhitening the leadfields using the noise covariance\n');

    % compute the prewhitening matrix
    if ~isempty(noiselambda)
      ft_info('using a regularized noise covariance matrix\n');
      % note: if different channel types are present, one should probably load the diagonal with channel-type specific stuff
      [U,S,V] = svd(C+eye(size(C))*noiselambda);
    else
      [U,S,V] = svd(C);
    end

    Tol     = 1e-12;
    diagS   = diag(S);
    sel     = find(diagS>Tol.*diagS(1));
    P       = diag(1./sqrt(diag(S(sel,sel))))*U(:,sel)'; % prewhitening matrix
    A       = P*A; % prewhitened leadfields
    C       = eye(size(P,1)); % prewhitened noise covariance matrix
  end

  if doscale
    % estimate sourcecov such that trace(ARA')/trace(C) = 1 (see
    % http://martinos.org/mne/manual/mne.html. In the case of prewhitening
    % C reduces to I (and then lambda^2 ~ 1/SNR); note that in mixed
    % channel type covariance matrices prewhitening should be applied in
    % order for this to make sense (otherwise the diagonal elements of C
    % have different units)
    ft_info('scaling the source covariance\n');
    scale = trace(A*(R*A'))/trace(C);
    R     = R./scale;
  end

  if ~isempty(snr)
    % the regularisation parameter can be estimated from the noise covariance,
    % see equation 6 in Lin et al. 2004
    lambda = trace(A * R * A')/(trace(C)*snr^2);
  end

  if dowhiten
    % as documented on MNE website, this is replacing the part of the code below, it gives
    % more stable results numerically.
    Rc      = chol(R, 'lower');
    [U,S,V] = svd(A * Rc, 'econ');
    s  = diag(S);
    ss = s ./ (s.^2 + lambda);
    w  = Rc * V * diag(ss) * U';

    % unwhiten the filters to bring them back into signal subspace
    w = w*P;

  else
    % equation 5 from Lin et al 2004 (this implements Dale et al 2000, and Liu et al. 2002)
    denom = (A*R*A'+(lambda^2)*C);
    if cond(denom)<1e12
      w = R * A' / denom;
    else
      ft_info('taking pseudo-inverse due to large condition number\n');
      w = R * A' * pinv(denom);
    end
  end

end % if empty noisecov

% for each of the timebins, estimate the source strength
w_x = H_tot*w; % The Harm matrix project the Harmony level filters back to the source locations

if isreal(dat)
  ft_info('the input consists of time-series: computing the dipole moments\n')
  mom = w_x * dat;
  mom_ind = 1;
elseif size(dat,1)==size(dat,2)&&sum(sum(dat-dat'))<10^-5*sum(diag(dat))
  ft_info('the input consists of a cross-spectral density: computing source-level power\n')
  pow_tot = real(sum((w_x*dat).*w_x,2));
  pow = 0;
  for j = 1:num_dip
    pow = pow + pow_tot((1 + (j-1)*end/num_dip):(j*end/num_dip));
    mom_ind = 0;
  end
else
  ft_info('the input consists of Fourier coefficients: computing source-level Fourier coefficients\n')
  mom = w_x*dat; % The Harm matrix project the Harmonic activity back to the source locations
  mom_ind = 1;
end

% assign the estimated source strength to each dipole
if mom_ind == 1
  n = 1;
  for i=1:size(sourcemodel.pos,1)
    cbeg = n;
    cend = n + size(sourcemodel.leadfield{i}, 2) - 1;
    estimate.mom{i} = mom(cbeg:cend,:);
    n = n + size(sourcemodel.leadfield{i}, 2);
  end
end

% compute power (over the three orientations) at each location and for each time

if mom_ind == 1
  estimate.pow = nan(size(sourcemodel.pos,1), size(dat,2));
  for i=1:size(sourcemodel.pos,1)
    estimate.pow(i,:) = sum(abs(estimate.mom{i}).^2, 1);
  end
else
  estimate.pow = pow;
end

% deal with keepfilter option
if keepfilter && ~hasfilter
  % spatial filters have been computed, store them in the output
  % re-assign spatial filter to conventional 1 cell per dipole location
  n = 1;
  for i=1:size(sourcemodel.pos,1)
    cbeg = n;
    cend = n + size(sourcemodel.leadfield{i}, 2) - 1;
    estimate.filter{i} = w(cbeg:cend,:);
    n    = n + size(sourcemodel.leadfield{i}, 2);
  end
elseif keepfilter
  estimate.filter = sourcemodel.filter;
end

% deal with noise covariance
% if ~isempty(noisecov) && ~hasfilter
%   % compute estimate of the projected noise
%   n = 1;
%   for i=1:size(sourcemodel.pos,1)
%     cbeg = n;
%     cend = n + size(sourcemodel.leadfield{i}, 2) - 1;
%     estimate.noisecov{i} = w(cbeg:cend,:)*noisecov*w(cbeg:cend,:)';
%     n    = n + size(sourcemodel.leadfield{i}, 2);
%   end
% elseif ~isempty(noisecov)
%   % compute estimate of the projected noise
%   for i=1:size(sourcemodel.pos,1)
%     estimate.noisecov{i} = sourcemodel.filter{i}*noisecov*sourcemodel.filter{i}';
%   end
%end

% reassign the estimated values over the inside and outside grid positions
estimate.inside  = originside;
estimate.pos     = origpos;
if isfield(estimate, 'pow')
  estimate.pow( originside,:) = estimate.pow;
  estimate.pow(~originside,:) = nan;
end
if isfield(estimate, 'mom')
  estimate.mom( originside) = estimate.mom;
  estimate.mom(~originside) = {[]};
end
if isfield(estimate, 'noisecov')
  estimate.noisecov( originside) = estimate.noisecov;
  estimate.noisecov(~originside) = {[]};
end
if isfield(estimate, 'filter')
  estimate.filter( originside) = estimate.filter;
  estimate.filter(~originside) = {[]};
end


text:
function r8vec_permute_test ( )

%*****************************************************************************80
%
%% RR8VEC_PERMUTE_TEST tests R8VEC_PERMUTE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 October 2014
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  p = [ 2, 4, 5, 1, 3 ];
  x = [ 1.0, 2.0, 3.0, 4.0, 5.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_PERMUTE_TEST\n' );
  fprintf ( 1, '  R8VEC_PERMUTE permutes an R8VEC.\n' );

  r8vec_print ( n, x, '  Original array X[]:' )

  i4vec_print ( n, p, '  Permutation vector P[]:' )

  x = r8vec_permute ( n, x, p );

  r8vec_print ( n, x, '  Permuted array X[P[]]:' );

  return
end


text:
function Trans = PUMA_TransMatrices(alpha, a, d, theta)
% This creates the Robot Manipulator Link(s) Transformation Matrix.
% The entries of the DH Table are passed in as the parameter.

%% Handle errors (if any) for the parameters passed into the function

% Check for the dimensions of the parameters passed
[~, c(1)] = size(alpha);
[~, c(2)] = size(a);
[~, c(3)] = size(d);
[~, c(4)] = size(theta);

if (c(1) ~= c(2) || c(1) ~= c(3) || c(1) ~= c(4))
    error('Invalid arguments. Size of all the arguments should be same')
    return
end

%% Generate the Transformation Matrices from the DH-parameters passed

% Find the Transformation matrix from the first row of DH-Table
% This will give us the position and orientation of the end of first link
% of the Robot Manipulator
Trans(:,:,1) = PUMA_Transformation(alpha(1), a(1), d(1), theta(1));

% Find the Transformation matrix from the rest of the rows of DH-Table
% This will give us the position and orientation of the rest of the links
% of the Robot Manipulator
for i = 2:c(1)
    Trans(:,:,i)= PUMA_Transformation(alpha(i), a(i), d(i), theta(i), Trans(:,:, i-1));
end

end

author: YashBansod

repo: Robotics-Planning-Dynamics-and-Control


text:
function value = f1 ( x )

%*****************************************************************************80
%
%% F1 evaluates a function that can be used for Chebyshev interpolation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 September 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(), one or more points where the function is to be evaluated.
%
%    Output, real VALUE(), the function values.
%
  value = cos ( 2.0 * pi * x ) .* sin ( 3.0 * pi * x );

  return
end


text:
%DEMO_SURVIVAL_COMPARISON  Survival model comparison
%
%  Description:
%
%    By using kfc-validation and Bayesian bootstrap we compare the
%    predictive ability of different models by estimating various
%    assessment statistics .
%
%    We will compare two Cox proportional hazars model, the first
%    model will have less covariates than the second model.
%
%    The censoring indicator ye is
%
%      ye = 0 for uncensored event
%      ye = 1 for right censored event.
%
%    Example data set is leukemia survival data in Northwest
%    England presented in (Henderson, R., Shimakura, S., and Gorst,
%    D. (2002). Modeling spatial variation in leukemia survival
%    data. Journal of the American Statistical Association,
%    97:965–972). Data set was downloaded from
%    http://www.math.ntnu.no/%7Ehrue/r-inla.org/examples/leukemia/leuk.dat
%
%  See also  DEMO_SURVIVAL_COMPARISON2, DEMO_SURVIVAL_COXPH

% Copyright (C) 2012 Ernesto Ulloa, Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

%% First load data
S = which('demo_survival_weibull');
L = strrep(S,'demo_survival_weibull.m','demodata/leukemia.txt');
leukemiadata=load(L);

% leukemiadata consists of:
% 'time', 'cens', 'xcoord', 'ycoord', 'age', 'sex', 'wbc', 'tpi', 'district'

% survival times
y=leukemiadata(:,1);
% scale survival times
y=y/max(y);

ye=1-leukemiadata(:,2); % event indicator, ye = 0 for uncensored event
                        %                        ye = 1 for right censored event

%  we choose for the first model: 'age' and 'sex'covariates
x01=leukemiadata(:,5:6);
x1=x01;

%  we choose for the second model: 'age', 'sex', 'wbc', and 'tpi' covariates
x02=leukemiadata(:,5:8);
x2=x02;

% normalize continuous covariates

x1(:,1)=normdata(x01(:,1));
x2(:,[1 3:4])=normdata(x02(:,[1 3:4]));

[n1, nin1]=size(x1);
[n2, nin2]=size(x2);

% number of time intervals
ntime=50;
% create finite partition of time axis
S=linspace(0,max(y)+0.001,ntime+1);

%% obtain predictions

% Create the covariance functions
pl = prior_t('s2',1, 'nu', 4);
pm = prior_t('s2',1, 'nu', 4);

% covariance for hazard function
gpcfh1 = gpcf_sexp('lengthScale', 1, 'magnSigma2', 1.1, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);
gpcfh2 = gpcf_sexp('lengthScale', 1, 'magnSigma2', 1.1, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);

% covariance for proportional part
gpcf1 = gpcf_sexp('lengthScale', ones(1,size(x1,2)), 'magnSigma2', 1.2, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);
gpcf2 = gpcf_sexp('lengthScale', ones(1,size(x2,2)), 'magnSigma2', 1.2, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);

% Create the likelihood structure
lik = lik_coxph('S', S);

gp1 = gp_set('lik', lik, 'cf', {gpcfh1 gpcf1}, 'jitterSigma2', 1e-6, 'comp_cf', {[1] [2]});
gp2 = gp_set('lik', lik, 'cf', {gpcfh2 gpcf2}, 'jitterSigma2', 1e-6, 'comp_cf', {[1] [2]});

% Set the approximate inference method to Laplace
gp1 = gp_set(gp1, 'latent_method', 'Laplace');
gp2 = gp_set(gp2, 'latent_method', 'Laplace');

opt=optimset('TolFun',1e-2,'TolX',1e-4,'Display','iter','Derivativecheck','off');

% obtain predictions for both models using kfc-validation

%* first we set tau
tt=0.1:.1:1;

% set D event indicator vector for each time in tt (Di=0 if i experienced
% the event before tau and Di=1 otherwise)
% Also we set YY, the observed time vector for each time value in tt
for i=1:size(tt,2)
  for i2=1:size(ye,1)
    if y(i2)>tt(i)
      yytemp(i2)=tt(i);
      Dtemp(i2)=1;
    else
      if ye(i2)==1
        Dtemp(i2)=1;
      else
        Dtemp(i2)=0;
      end
      yytemp(i2)=y(i2);
    end
  end
  yyi{i}=yytemp';
  Di{i}=Dtemp';
end
for i=1:size(Di,2)
  D(:,i)=Di{i};
end
for i=1:size(yyi,2)
  yy(:,i)=yyi{i};
end

% set time vector to make predictions
yt=bsxfun(@times,ones(size(y)),tt);

% Obtain predictions
% (This takes several minutes)
crit1=gp_kfcv_cdf(gp1,x1,y,'z',D,'yt',yt,'opt',opt);
crit2=gp_kfcv_cdf(gp2,x2,y,'z',D,'yt',yt,'opt',opt);

%% Calculate statics and compare models
% MODEL COMPARISON

%% AUC
% AUC for Binary outcomes P(Pi>Pj | Di=1,Dj=0)
[auc1,fps1,tps1]=aucs(crit1(:,length(tt)),D(:,length(tt)));
[auc2,fps2,tps2]=aucs(crit2(:,length(tt)),D(:,length(tt)));

fprintf('AUC at end of study for model 1: %.3f \n', auc1);
fprintf('AUC at end of study for model 2: %.3f \n', auc2);
hold on
plot(fps1,tps1,'b')
plot(fps2,tps2,'r')
title('ROC curve')
legend('model 1', 'model 2',4)
xlabel('False positives')
ylabel('True positives')
hold off

%% Harrell's C
% Obtain for both models Binary AUC(t) = P(Pi>Pj | Di(t)=1,Dj(t)=0) and
% Harrell's C(t) = P(Pi>Pj | Di(ti)=1, ti<tj, ti<tt) for every element of tt
ct1=hct(crit1,yy,D,tt);
ct2=hct(crit2,yy,D,tt);
auct1=auct(crit1,yy,D,tt);
auct2=auct(crit2,yy,D,tt);
c=[ct1 ct2];

% Plot for both models Harrells C in function of time
plot(tt,c(:,1),'r');
hold on;
plot(tt,c(:,2),'g');
legend('Old model','New model')
title('Harrolds C in function of time ');
xlabel('Time');
ylabel('Harrell''s C');
hold off;

%% Estimated density
% Use bayesian bootsrap to obtain Harrells (C1-C2) statistic density at tt=1
[c1,bb1]=hcs(crit1(:,end),y,ye,1,'rsubstream',1);
[c2,bb2]=hcs(crit2(:,end),y,ye,1,'rsubstream',1);
title('Estimated density of C2-C1')
lgpdens(bb2-bb1)
xlabel('Difference in Harrell''s C statistics (C2-C1)');

% We integrate the (C1-C2) estimated density in the (0,inf) interval
zc=lgpdens_cum(bb2-bb1,0,inf);
fprintf('Estimated c statistics for model 1 and 2 respectively:  %.3f, %.3f \n', c1, c2);
fprintf('cumulative probability in the (0,inf) interval:  %.2f \n', zc);

%% IDI
%Estimate R^2 for both models, idi, its density and the cumulative
%probability in the (0,inf) interval, al at time 1

[idi,r1,r2,bbid] = idis(crit1(:,end),crit2(:,end),'rsubstream',1);
zidi=lgpdens_cum(bbid,0,inf);
title('IDI estimated density')
lgpdens(bbid)

fprintf('R^2 statistic for model 1: %.3f \n', r1);
fprintf('R^2 statistic for model 2: %.3f \n', r2);

fprintf('Estimated idi: %.3f ', idi);
fprintf('cumulative probability in the (0,inf) interval: %.2f\n', zidi);

%% EXT AUC

% Ext_AUC for different subsets of tt
Indxtmp{1}=1:1:size(tt,2);
Indx{1}=1:1:size(tt,2);
j=2;
k=round(size(tt,2)/2);
for i=2:k
  Indxtmp{i}=1:i:size(tt,2);
  if length(Indxtmp{i})~=length(Indxtmp{i-1})
    Indx{j}=Indxtmp{i};
    j=j+1;
  end
end

for i=1:size(Indx,2)
  l(i)=length(Indx{i});
end

for i=1:size(Indx,2)
  ea1(i) = ext_auc(crit1(:,Indx{i}),tt(:,Indx{i}),tt(:,Indx{i}(size(Indx{i},2))));
  ea2(i) = ext_auc(crit2(:,Indx{i}),tt(:,Indx{i}),tt(:,Indx{i}(size(Indx{i},2))));
end

extauc1 = ext_auc(crit1,tt,tt(:,end));
extauc2 = ext_auc(crit2,tt,tt(:,end));

fprintf('ExtAUC at end of study for model 1: %.3f \n', extauc1);
fprintf('ExtAUC at end of study for model 2: %.3f \n', extauc2);


text:
function Area=unitVectors2SpherTriangArea(u,r)
%%UNITVECTORS2SPHERTRIANGAREA Given unit vectors pointing from the center
%     of a 3D sphere to three points on the sphere, as well as the radius
%     of the sphere, determine the surface area of a spherical triangle.
%
%INPUTS: u A 3X3 matrix such that u(:,i) is the ith unit vector pointing
%          from the center of the sphere to the surface of the sphere.
%        r The scalar radius of the sphere. If this is omitted or an empty
%          matrix is passed, the default of 1 is used.
%
%OUTPUTS: Area The area of the specified triangle on the surface of a
%              sphere of radius r.
%
%Equation 1 in [1] is implemented. The absolute value of the atan function
%is taken so that the result is always positive.
%
%EXAMPLE:
%We find the surface area of a quadrant of a unit sphere. Since the surface
%area of an entire sphere is 4*pi*r^2, the area of a quadrant will be pi/2.
%The relative error is zero.
% u=[[1;0;0],[0;1;0],[0;0;1]];
% r=1;
% Area=unitVectors2SpherTriangArea(u,r);
% trueArea=pi/2;
% RelErr=(Area-trueArea)./trueArea
%
%REFERENCES:
% [1] F. Eriksson, "On the measure of solid angles," Mathematics Magazine,
%     vol. 63, no. 3, pp. 184-187, 1990.
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(r))
    r=1;
end

a=u(:,1);
b=u(:,2);
c=u(:,3);

Area=2*r^2*abs(atan(dot(a,cross(b,c))/(1+dot(b,c)+dot(c,a)+dot(a,b))));

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 [A, B] = reflexKronApprox( K, m, n )
%
%       [A, B] = reflexKronApprox( K, m, n );
%
% computes a kronecker sum approximation to the blurring matrix K
% that arises from the input PSF under reflexive boundary conditions.
% This approximation is done a-la Nagy-Ng-Perrone (see paper for details).
%
%  Input:
%         K - psfMatrix object
%         m - size of matrix A (assumed square)
%         n - size of matrix B (assumed square)
%
%  Output:
%      Matrices A and B such that K \aprrox A \otimes B.
%

%  L. Perrone, 4/28/02

%  Modifications:
%  5/25/02, J. Nagy
%           Cosmetic changes to incorporate into RestoreTools
%
%
%  11/17/02, J. Nagy
%            This was designed for image processing problems, where
%            it is common to use lexicographical (row) ordering.
%            But @kronMatrix functions where designed using vec (column)
%            ordering.  This inconsistency has been fixed.

P1 = K.psf;
P2 = P1.image;
PSF = P2{1};
c1 = P1.center;
center = c1{1};

[mp, np] = size(PSF);

if ( mp ~= np )
  error('For now, we expect PSF to be square')
end

%
% Compute weighted PSF.
%
c = zeros(mp,1);
c(1) = mp;
c(2:2:end) = 1;
R = chol( toeplitz(c) );

Phat = R*PSF*R';

%
% Compute SVD of weighted PSF, which is then used to construct
% the separable approximation.
%
[U,S,V] = svd( Phat );

%
% check to make sure first column looks like
% a Gaussian, and is not inverted.
%
minU = abs(min(min(U(:,1))));
maxU = max(max(abs(U(:,1))));
if minU == maxU
  U = -U;
  V = -V;
end

%
% Construct approximation.
%
a = R \ ( U(:,1) * sqrt(S(1,1)) );
b = R \ ( V(:,1) * sqrt(S(1,1)) );

%
%  This construction corresponds to lexicographical ordering,
%  but kronMatrix does everything corresponding to vec ordering.
%  Thus, these two do not work together.
%%
%%A = build_toep(a, center(1), n) + buildHank(a, center(1), n);
%%B = build_toep(b, center(2), m) + buildHank(b, center(2), m);

A = build_toep(b, center(2), m) + buildHank(b, center(2), m);
B = build_toep(a, center(1), n) + buildHank(a, center(1), n);

repo: IRtools


text:
function [K, Kbase, Prinv, Pqrinv, P, n2] = ggKernCompute(kern, x, x2)

% GGKERNCOMPUTE Compute the GG kernel given the parameters and X.
% FORMAT
% DESC computes the kernel parameters for
%	the gaussian gaussian kernel given inputs associated with rows and
%	columns.
% RETURN K : the kernel matrix computed at the given points.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : the input matrix associated with the rows of the kernel.
% ARG x2 the input matrix associated with the columns of the kernel.
%
% FORMAT
% DESC computes the kernel matrix for the
%	gaussian kernel given a design matrix of inputs.
%	 Returns:
% RETURN k : the kernel matrix computed at the given points.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG X : input data matrix in the form of a design matrix.
%
% SEEALSO : ggKernParamInit, kernCompute, kernCreate, ggKernDiagCompute
%
% COPYRIGHT : Mauricio A. Alvarez and Neil D. Lawrence, 2008
%
% MODIFICATIONS : Mauricio A. Alvarez, 2009.

% KERN

if nargin < 3
  x2 = x;
end

Pr  = kern.precisionU;
Pqr = kern.precisionG;
Prinv = 1./Pr;
Pqrinv = 1./Pqr;
Pinv = Prinv + 2*Pqrinv;
P = 1./Pinv;
if kern.isArd,
    sqrtP = sparseDiag(sqrt(P));
    sqrtPx = x*sqrtP;
    sqrtPx2 = x2*sqrtP;
    n2 = dist2(sqrtPx, sqrtPx2);
    Kbase = exp(-0.5*n2);
else
    n2 = dist2(x, x2);
    Kbase = exp(-0.5*P*n2);
end
K = kern.sigma2Latent*kern.sensitivity^2*...
    Kbase;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat