text:
function jdate = julian (month, day, year)

% Julian date

% Input

%  month = calendar month [1 - 12]
%  day   = calendar day [1 - 31]
%  year  = calendar year [yyyy]

% Output

%  jdate = Julian date

% special notes

%  (1) calendar year must include all digits

%  (2) will report October 5, 1582 to October 14, 1582
%      as invalid calendar dates and stop

% Orbital Mechanics with MATLAB

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

y = year;
m = month;
b = 0;
c = 0;

if (m <= 2)
   y = y - 1;
   m = m + 12;
end

if (y < 0)
   c = -.75;
end

% check for valid calendar date

if (year < 1582)
   % null
elseif (year > 1582)
   a = fix(y / 100);
   b = 2 - a + floor(a / 4);
elseif (month < 10)
   % null
elseif (month > 10)
   a = fix(y / 100);
   b = 2 - a + floor(a / 4);
elseif (day <= 4)
   % null
elseif (day > 14)
   a = fix(y / 100);
   b = 2 - a + floor(a / 4);
else
   clc; home;

   fprintf('\n\n  this is an invalid calendar date!!\n');

   keycheck;

   return;
end

jd = fix(365.25 * y + c) + fix(30.6001 * (m + 1));

jdate = jd + day + b + 1720994.5;

author: goGPS-Project


text:
%
% Author: Javier Lopez-Calderon
%
function cleanerplab

p = which('eegplugin_erplab');
p = p(1:findstr(p,'eegplugin_erplab.m')-1);
direrp = {'functions', 'pop_functions', 'Guis'};
recycle on;

for i=1:length(direrp)
        tempname1 = fullfile(p,direrp{i},'*.m~');
        tempname2 = fullfile(p,direrp{i},'*.asv');
        delete(tempname1)
        delete(tempname2)
end
recycle off;
disp('Clear!')

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function x = loop(x, continue_fcn, f, cleanup_fcn)

% x = loop(x, continue_fcn, f, cleanup_fcn)
%
% Looping function. Provides the ability to treating "looping" as a
% function instead of a command.
%
% See functional_programming_examples.m for a discussion of this function.
%
% Inputs:
% x            - Initial state (can be cell array of arguments to f)
% continue_fcn - Continue function, returns true iff the loop should go on
% f            - Function of the state (x) to run every iteration
% cleanup_fcn  - Function to select what to return from the state when the
%                looping is complete (optional)
%
% Outputs:
% x        - The updated state. If a cleanup function is supplied, this
%            will be cleanup(x).
%
% Fibonacci example:
%
% fib = @(n) loop({1, 0, 1}, ...          % {k, fib(k-1), fib(k)}
%                 @(x) x{1} < n, ...      % While k < n
%                 @(x) {x{1} + 1, ...     %   Increase k
%                       x{3}, ...         %   Store fib(k-1)
%                       x{3} + x{2}}, ... %   Store fib(k-1) + fib(k-2)
%                 @(x) x{3})              % End, returning only fib(k)
%
% arrayfun(fib, 1:10)
%
% Tucker McClure
% Copyright 2013 The MathWorks, Inc.

    while tern(iscell(x), @()continue_fcn(x{:}), @()continue_fcn(x))
        if iscell(x)
            x = f(x{:});
        else
            x = f(x);
        end
    end

    if nargin == 4
        if iscell(x)
            x = cleanup_fcn(x{:});
        else
            x = cleanup_fcn(x);
        end
    end

end

author: Sable

repo: mcbench-benchmarks


text:
% DEFINE_BIFF script to define the known BIFF chunk types
%
% This function should not be called separately, it is intended to
% be used from withing read_biff and write_biff.
%
% See also READ_BIFF, WRITE_BIFF

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

CAT_.desc = {'CAT ', 'cat', 'unknown'};
LIST.desc = {'LIST', 'list',    'unknown'};
BIFF.desc = {'BIFF', 'biff',    'group'};

BIFF.SEMG.desc = {'SEMG', 'semg',   'group'};

BIFF.SEMG.VERS.desc = {'VERS', 'version',       'group'};
BIFF.SEMG.DINF.desc = {'DINF', 'data_info',     'group'};
BIFF.SEMG.EXPI.desc = {'EXPI', 'experiment',        'group'};
BIFF.SEMG.MEAS.desc = {'MEAS', 'measurement',       'group'};
BIFF.SEMG.SGRD.desc = {'SGRD', 'semggrid',      'group'};
BIFF.SEMG.PATI.desc = {'PATI', 'patient',       'group'};

BIFF.SEMG.VERS.ID__.desc = {'ID  ', 'id',       'string'};
BIFF.SEMG.VERS.LIC_.desc = {'LIC ', 'license',      'string'};

BIFF.SEMG.DATI.CHAN.desc = {'CHAN', 'unknown'       'unknown'};
BIFF.SEMG.DATI.MODE.desc = {'MODE', 'unknown'       'unknown'};
BIFF.SEMG.DATI.TYPE.desc = {'TYPE', 'unknown'       'unknown'};
BIFF.SEMG.DATI.FS__.desc = {'FS  ', 'unknown'       'unknown'};
BIFF.SEMG.DATI.LABL.desc = {'LABL', 'unknown'       'unknown'};
BIFF.SEMG.DATI.UNIT.desc = {'UNIT', 'unknown'       'unknown'};
BIFF.SEMG.DATI.FILT.desc = {'FILT', 'unknown'       'unknown'};

BIFF.SEMG.EXPI.NAME.desc = {'NAME', 'name',     'string'};
BIFF.SEMG.EXPI.EXAM.desc = {'EXAM', 'examination',  'string'};
BIFF.SEMG.EXPI.EXRC.desc = {'EXRC', 'exercise',     'string'};
BIFF.SEMG.EXPI.RPHS.desc = {'RPHS', 'ref_physician',    'string'};
BIFF.SEMG.EXPI.DATE.desc = {'DATE', 'date',     'unknown'};
BIFF.SEMG.EXPI.TIME.desc = {'TIME', 'time',     'unknown'};

BIFF.SEMG.MEAS.TIME.desc = {'TIME', 'time',     'string'};
BIFF.SEMG.MEAS.ORFN.desc = {'ORFN', 'filename',     'string'};
BIFF.SEMG.MEAS.MUSC.desc = {'MUSC', 'muscle',       'string'};
BIFF.SEMG.MEAS.SIDE.desc = {'SIDE', 'side',     'unknown'};
BIFF.SEMG.MEAS.FNUM.desc = {'FNUM', 'filenumber',   'unknown'};
BIFF.SEMG.MEAS.TEMP.desc = {'TEMP', 'temperature',  'unknown'};
BIFF.SEMG.MEAS.ANOT.desc = {'ANOT', 'annotation',   'string'};
BIFF.SEMG.MEAS.FORC.desc = {'FORC', 'force',        'string'};

BIFF.SEMG.SGRD.NPD_.desc = {'NPD ', 'npd',      'unknown'};
BIFF.SEMG.SGRD.NML_.desc = {'NML ', 'nml',      'unknown'};
BIFF.SEMG.SGRD.IED_.desc = {'IED ', 'ied',      'unknown'};
BIFF.SEMG.SGRD.FIBD.desc = {'FIBD', 'fibre_dir',    'unknown'};
BIFF.SEMG.SGRD.POS_.desc = {'POS ', 'tle_pos',      'unknown'};
BIFF.SEMG.SGRD.CABL.desc = {'CABL', 'cable_pos',    'unknown'};
BIFF.SEMG.SGRD.GRID.desc = {'GRID', 'grid',     'int16vec'};

BIFF.SEMG.PATI.UPID.desc = {'UPID', 'unique_id',    'unknown'};
BIFF.SEMG.PATI.GEND.desc = {'GEND', 'gendre',       'unknown'};
BIFF.SEMG.PATI.BDAT.desc = {'BDAT', 'birthdate',    'unknown'};
BIFF.SEMG.PATI.HPID.desc = {'HPID', 'hospital_id',  'unknown'};
BIFF.SEMG.PATI.ANAM.desc = {'ANAM', 'anamnesis',    'unknown'};

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function obj = modelModifyRegionAnalysis(obj,varargin)
	% Remove PCAs in a particular region or exclude from preprocessing, etc.
	% Biafra Ahanonu
	% started: 2015.09.14 14:46:02
	% inputs
		%
	% outputs
		%

	% changelog
		% 2017.01.14 [20:06:04] - support switched from [nSignals x y] to [x y nSignals]
		% 2021.06.18 [21:41:07] - added modelVarsFromFilesCheck() to check and load signals if user hasn't already.
		% 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
		% 2022.02.28 [13:01:20] - Fix modelVarsFromFilesCheck and modelVarsFromFiles recursion edge case.
		% 2022.08.26 [15:00:28] - Misc code standard updates.
	% TODO
		%

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	options.defaultOption = 'loadPreviousSelections';

	% get options
	options = getOptions(options,varargin);
	%========================

	if obj.guiEnabled==1
		scnsize = get(0,'ScreenSize');
		[fileIdxArray, idNumIdxArray, nFilesToAnalyze, nFiles] = obj.getAnalysisSubsetsToAnalyze();

		% [fileIdxArray, ok] = listdlg('ListString',obj.fileIDNameArray,'ListSize',[scnsize(3)*0.2 scnsize(4)*0.25],'Name','which folders to analyze?');
		scnsize = get(0,'ScreenSize');
		dlgSize = [scnsize(3)*0.7 scnsize(4)*0.8];
		usrIdxChoiceStr = {'selectRegions','modifyExistingRegions','runAlreadySelectedRegions','loadPreviousSelections'};
		[sel, ok] = listdlg('ListString',usrIdxChoiceStr,'ListSize',dlgSize);
		analysisToRun = usrIdxChoiceStr{sel};
	else
		if isempty(obj.foldersToAnalyze)
			fileIdxArray = 1:length(obj.fileIDNameArray);
		else
			fileIdxArray = obj.foldersToAnalyze;
		end
		analysisToRun = options.defaultOption;
	end
	nFolders = length(fileIdxArray);

	framesPerSecond = obj.FRAMES_PER_SECOND;
	for thisFileNumIdx = 1:length(fileIdxArray)
		try
			fileNum = fileIdxArray(thisFileNumIdx);
			obj.fileNum = fileNum;
			display(repmat('=',1,21))
			display([num2str(fileNum) '/' num2str(nFolders) ': ' obj.fileIDNameArray{obj.fileNum}]);

			% obj.modelVarsFromFilesCheck(fileNum);

			if strcmp(analysisToRun,'loadPreviousSelections')
				regionFile = getFileList(obj.inputFolders{obj.fileNum},obj.regionModSaveStr);
				% if no file, skip
				if isempty(regionFile)
					disp('No region analysis to load!')
					continue;
				end
				regionFile = regionFile{1};
				fprintf('loading: %s',regionFile)
				loadFile = load(regionFile,'roipolyRegion','roipolyCoords');
				obj.analysisROIArray{obj.fileNum} = loadFile.roipolyRegion;
				obj.validRegionModPoly{obj.fileNum} = loadFile.roipolyCoords;
				% NOT DONE!!!!!!!!
			end
			if strcmp(analysisToRun,'modifyExistingRegions')
				regionFile = getFileList(obj.inputFolders{obj.fileNum},obj.regionModSaveStr);
				% if no file, skip
				if isempty(regionFile)
					disp('No region analysis to load!')
					continue;
				end
				loadFile = load(regionFile,'roipolyRegion','roipolyCoords');
				obj.analysisROIArray{obj.fileNum} = loadFile.roipolyRegion;
				obj.validRegionModPoly{obj.fileNum} = loadFile.roipolyCoords;
			end
			[inputSignals, inputImages, signalPeaks, signalPeaksArray] = modelGetSignalsImages(obj,'returnType','raw');
			% inputImages = thresholdImages(inputImages,'binary',0)/3;
			% [inputSignals inputImages signalPeaks signalPeaksArray] = modelGetSignalsImages(obj);

			% display cellmap and ask user to select a region
			[~, ~] = openFigure(obj.figNoAll, '');
			clf
			colormap gray
			subX = 3; subY = 3;
			if sum(strcmp(analysisToRun,{'selectRegions','modifyExistingRegions','runAlreadySelectedRegions','loadPreviousSelections'}))
				subplot(subX,subY,3);
					numPeakEvents = sum(signalPeaks,2);
					numPeakEvents = numPeakEvents/size(signalPeaks,2)*framesPerSecond;
					size(numPeakEvents)
					size(signalPeaks)
					[groupedImagesRates] = groupImagesByColor(inputImages,numPeakEvents);
					thisCellmap = createObjMap(groupedImagesRates);
					imagesc(thisCellmap);
					title('firing rate cell map')
					% colormap(obj.colormap);
				subplot(subX,subY,6);
					movieList = getFileList(obj.inputFolders{obj.fileNum}, 'concat');
					if ~isempty(movieList)
						movieFrame = loadMovieList(movieList{1},'convertToDouble',0,'frameList',1:2);
						movieFrame = squeeze(movieFrame(:,:,1));
						% movieFrameMean = mean(movieFrame(:));
						% imagesc(imadjust(movieFrame+cast(thisCellmap*10000,class(movieFrame))));
						E = normalizeVector(double(movieFrame),'normRange','zeroToOne')/2;
						% Comb = E;
						I = cast(normalizeVector(thisCellmap,'normRange','zeroToOne')*mean(movieFrame(:)),class(movieFrame));
						I = normalizeVector(double(I),'normRange','zeroToOne')/3;
						Comb(:,:,1) = E; % red
						Comb(:,:,2) = E+I; % green
						Comb(:,:,3) = E; % blue
						imagesc(Comb)
						% inputImagesThresholded = thresholdImages(inputImages,'binary',0)/3;
						% zeroMap = zeros(size(movieFrame));
						% oneMap = ones(size(movieFrame));
						% warning off
						% colorOverlay = imshow(cat(3, zeroMap, oneMap, zeroMap));
						% warning on
						% thisCellmap = createObjMap(inputImagesThresholded);
						% set(colorOverlay, 'AlphaData', thisCellmap);
					end
					title('movie frame + cell map')
				subplot(subX,subY,[1 2 4 5 7 8]);
					try obj.valid{obj.fileNum}.(obj.signalExtractionMethod).manual;check.manual=1; catch check.manual=0; end
					% if isfield(obj.valid,obj.signalExtractionMethod)&isfield(obj.valid.(obj.signalExtractionMethod),'manual')
					if check.manual==1
						[groupedImagesRates] = groupImagesByColor(inputImages,obj.valid{obj.fileNum}.(obj.signalExtractionMethod).manual);
					else
						[groupedImagesRates] = groupImagesByColor(inputImages,obj.validAuto{obj.fileNum});
					end
					thisCellmap = createObjMap(groupedImagesRates);
					% obj.rawImagesFiltered{obj.fileNum} = thisCellmap;
					imagesc(thisCellmap+0.1);
					% imagesc(obj.rawImagesFiltered{obj.fileNum}+0.1);
					title('draw selection on this image')
					% colormap(obj.colormap)
					ciapkg.overloaded.suptitle(strrep(obj.folderBaseSaveStr{obj.fileNum},'_',' | '))
					if ~sum(strcmp(analysisToRun,{'runAlreadySelectedRegions','loadPreviousSelections'}))
						% [polyCoords] = obj.validRegionModPoly{obj.fileNum};
						if strcmp(analysisToRun,'modifyExistingRegions')
							[polyCoords] = obj.validRegionModPoly{obj.fileNum};
							h = impoly(gca,polyCoords);
							h.wait
							polyCoords = h.getPosition;
							[obj.analysisROIArray{obj.fileNum}, xpoly, ypoly] = roipoly(thisCellmap+0.1,polyCoords(:,1),polyCoords(:,2));
						else
							[obj.analysisROIArray{obj.fileNum}, xpoly, ypoly] = roipoly;
						end
						obj.validRegionModPoly{obj.fileNum} = [xpoly ypoly];
					end
					% imagesc(obj.rawImagesFiltered{obj.fileNum}.*obj.analysisROIArray{obj.fileNum})
					imagesc(thisCellmap.*obj.analysisROIArray{obj.fileNum})
					drawnow;
			end
			if strcmp(analysisToRun,'runAlreadySelectedRegions')
				disp('Using previous ROI')
			end
			inputImagesROI = obj.analysisROIArray{obj.fileNum};
			inputImagesThres = thresholdImages(inputImages,'waitbarOn',1,'binary',1);
			% signalInROI = squeeze(nansum(nansum(bsxfun(@times,inputImages,permute(inputImages,[2 3 1])),1),2));
			disp('finding ROIs inside region...')
			signalInROI = squeeze(sum(sum(bsxfun(@times,inputImagesThres,inputImagesROI),1,'omitnan'),2,'omitnan'));

			% signalInROI = applyImagesToMovie(inputImages,permute(inputImages,[2 3 1]), 'alreadyThreshold',1);
			signalsToKeep = signalInROI~=0;

			if isempty(obj.validManual{obj.fileNum})
				validRegionModHere = obj.valid{obj.fileNum}.(obj.signalExtractionMethod).auto(:)&signalsToKeep(:);
			else
				validRegionModHere = obj.valid{obj.fileNum}.(obj.signalExtractionMethod).manual(:)&signalsToKeep(:);
				% if ~isempty(obj.validManual{obj.fileNum})
				% 	obj.validRegionMod{obj.fileNum} = obj.validAuto{obj.fileNum}(:)&signalsToKeep(:);
				% end
			end
			obj.validRegionMod{obj.fileNum} = validRegionModHere;
			obj.valid{obj.fileNum}.(obj.signalExtractionMethod).regionMod = validRegionModHere;

			[~, ~] = openFigure(obj.figNoAll, '');
				subplot(subX,subY,9);
				[filterImageGroups] = createObjMap(groupImagesByColor(inputImages,obj.validRegionMod{obj.fileNum}+1));
				% size(filterImageGroups)
				imagesc(filterImageGroups)
				title(sprintf('new cell map | %d cells',sum(validRegionModHere)))
				% colormap(obj.colormap)

			roipolyRegion = obj.analysisROIArray{obj.fileNum};
			roipolyCoords = obj.validRegionModPoly{obj.fileNum};

			thisDirSaveStr = [obj.inputFolders{obj.fileNum} filesep obj.folderBaseSaveStr{obj.fileNum}];
			savestring = [thisDirSaveStr obj.regionModSaveStr];
			display(['saving: ' savestring])
			save(savestring,'roipolyRegion','-v7.3','roipolyCoords');

			if sum(strcmp(analysisToRun,{'runAlreadySelectedRegions','loadPreviousSelections'}))
			else
				uiwait(ciapkg.overloaded.msgbox('press OK to move onto next folder','Success','modal'));
			end

		catch err
			display(repmat('@',1,7))
			disp(getReport(err,'extended','hyperlinks','on'));
			display(repmat('@',1,7))
		end
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
% this is tha main screen of the program, which includes every aspect of our project you can see the help for further information

function menu()
fig=figure('name','Cellular Project','menubar','none','numbertitle','off','units','normalized','color','w',...
    'position',[0.25 0.25 0.7 0.65])
h1=axes('units','normalized','position',[.01 0.3 0.3 .47],'xcolor',[1 1 1],'ycolor',[1 1 1])
t1=uicontrol('style','text','units','normalized','position',[0.4 0.8 0.3 0.1],'fontname','arial','fontsize',10,'string','Please Click On One Of The Following')
a1=uicontrol('style','listbox','string',[{'Traffic Cacluation'},{'Site Configuration'},{'Report'},{'Close'}],'units','normalized','position',[0.4 0.4 0.3 0.3],'fontname','comicsansms','fontsize',10,'callback',@push);
m1=uimenu('label','extra');
m2=uimenu('label','about');

uimenu(m2,'label','about me','callback',@me);
uimenu(m1,'label','maximize','callback',@men);
uimenu(m1,'label','traffic','callback',@tra,'separator','on');
uimenu(m1,'label','cellular','callback',@cel,'separator','off');
uimenu(m1,'label','close','callback',@cls,'separator','on');
m3=uimenu('label','help','callback',@help)
   axis(h1);
   ss=imread('ju1.jpg');
   imshow(ss)
	function help(varargin)
		open('help.doc')
	end
	function me(varargin)
str=sprintf(['Sanad Amin \n\n',...
    'Reg No:0057035 \n',...
    'Sanad_amin@hotmail.com'])
msgbox(str,'About Me','help')
    end
    function men(varargin)

       set(fig,'units','normalized','position',[0.01 0.01 1 0.95]) ;
    end
    function tra(varargin)
       run traffic
    end
    function cel(varargin)
       run cellular
    end
    function cls(varargin)
        close(gcf)

    end
    function push(varargin)
       val=get(a1,'value');
       if val==1;
         run traffic
       else if val==2;
               run cellular
		   else if val==3
				   open('report.doc')
			   end
		    if val==4;
                  close(gcf) ;
               end

           end
       end
    end

end


text:
% EX_MAXWELL_SRC_CAMEMBERT: solve Maxwell source problem in three quarters of the cylinder, where the exact solution is a singular function.

% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = 'geo_camembert.txt';

% Type of boundary conditions
problem_data.nmnn_sides   = [3 4 5 6];
problem_data.drchlt_sides = [1 2];

% Physical parameters
problem_data.c_mass  = @(x, y, z) ones(size(x));
problem_data.c_stiff = @(x, y, z) ones(size(x));

% Source and boundary terms
k = 1; % Constant that characterizes the singularity
problem_data.f = @(x, y, z) cat (1, ...
                          singular_function_maxwell (x, y, k), ...
                          zeros ([1, size(x)]));
problem_data.g = @(x, y, z, ind) zeros ([3, size(x)]);
problem_data.h = @(x, y, z, ind) cat (1, ...
                          singular_function_maxwell (x, y, k), ...
                          zeros ([1, size(x)]));

% Exact solution (optional)
problem_data.uex     = @(x, y, z) cat (1, ...
                          singular_function_maxwell (x, y, k), ...
                          zeros ([1, size(x)]));
problem_data.curluex = @(x, y, z) zeros ([3, size(x)]);

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.degree     = [2 2 2];     % Degree of the bsplines
method_data.regularity = [1 1 1];     % Regularity of the splines
method_data.nsub       = [3 3 3];     % Number of subdivisions
method_data.nquad      = [3 3 3];     % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
% Remove warnings for the degenerated elements
warning ('off', 'geopdes:jacdet_zero_at_quad_node')
warning ('off', 'geopdes:zero_measure_element')
[geometry, msh, space, u] = solve_maxwell_src (problem_data, method_data);

% 4) POST-PROCESSING
% 4.1) EXPORT TO PARAVIEW
output_file = 'maxwell_camembert_Deg2_Reg1_Sub3';

vtk_pts = {linspace(0, 1, 15), linspace(0, 0.99, 15), linspace(0, 1, 15)};
fprintf ('The result is saved in the file %s \n \n', output_file);
sp_to_vtk (u, space, geometry, vtk_pts, output_file, 'u')

% 4.2) Comparison with the exact solution
[error_hcurl, error_l2] = ...
    sp_hcurl_error (space, msh, u, problem_data.uex, problem_data.curluex)

%!demo
%! ex_maxwell_src_camembert

%!test
%! problem_data.geo_name = 'geo_camembert.txt';
%! problem_data.nmnn_sides   = [3 4 5 6];
%! problem_data.drchlt_sides = [1 2];
%! problem_data.c_mass  = @(x, y, z) ones(size(x));
%! problem_data.c_stiff = @(x, y, z) ones(size(x));
%! k = 1; % Constant that characterizes the singularity
%! problem_data.f = @(x, y, z) cat (1, ...
%!                           singular_function_maxwell (x, y, k), ...
%!                           zeros ([1, size(x)]));
%! problem_data.g = @(x, y, z, ind) zeros ([3, size(x)]);
%! problem_data.h = @(x, y, z, ind) cat (1, ...
%!                           singular_function_maxwell (x, y, k), ...
%!                           zeros ([1, size(x)]));
%! problem_data.uex     = @(x, y, z) cat (1, ...
%!                           singular_function_maxwell (x, y, k), ...
%!                           zeros ([1, size(x)]));
%! problem_data.curluex = @(x, y, z) zeros ([3, size(x)]);
%! method_data.degree     = [2 2 2];     % Degree of the bsplines
%! method_data.regularity = [1 1 1];     % Regularity of the splines
%! method_data.nsub       = [3 3 3];     % Number of subdivisions
%! method_data.nquad      = [3 3 3];     % Points for the Gaussian quadrature rule
%! warning ('off', 'geopdes:jacdet_zero_at_quad_node')
%! warning ('off', 'geopdes:zero_measure_element')
%! [geometry, msh, space, u] = solve_maxwell_src (problem_data, method_data);
%! [error_hcurl, error_l2] = sp_hcurl_error (space, msh, u, problem_data.uex, problem_data.curluex);
%! assert (msh.nel, 81)
%! assert (space.ndof, 820)
%! assert (error_l2, 0.0436842163001183, 1e-14)
%! assert (error_hcurl, 0.0436861291289525, 1e-14)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function expressions( varargin )

%EXPRESSIONS Declares one or more CVX expression holders.
%   EXPRESSIONS x1 x2 x3 ..., where x1, x2, x3, etc. are valid
%   variable names, declares multiple cvx expression holders. It is
%   exactly equivalent to issuing a separate EXPRESSION command
%   for each x1, x2, x3, ...
%
%   EXPRESSIONS allows the declaration of vector, matrix, and
%   array variables.
%
%   For more information about expression holders, see the help for
%   EXPRESSION or the CVX user guide.
%
%   Examples:
%      expressions x y z
%
%   See also EXPRESSION.

if nargin < 1,
    error( 'Incorrect syntax for EXPRESSIONS. Type HELP EXPRESSIONS for details.' );
elseif ~iscellstr( varargin ),
    error( 'All arguments must be strings.' );
end
for k = 1 : nargin,
    evalin( 'caller', [ 'expression ', varargin{k} ] );
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function filename = env_translatepath(filename)
% Translates platform-independent directories into a system-specific directories.
% SystemPath = env_translatepath(IndependentPath)
%
% BCILAB supports platform-independent paths for all its scripts and IO functions, which allows for
% script portability from system to system. It is especially important when data paths are mounted
% at different locations, depending on access mode and operating system. A side effect of portable
% path names is that there is one unique expression which computes any given data set, such as
% "flt_iir(flt_reref((flt_resample(io_loadset('data:/Projects/Test/test1.vhdr'),[],200))),[],[4 6 25
% 30])", and this in turn allows to share the same data set caches (which are indexed by expression)
% across machines and operating systems, minimizing redundant computations.
%
% Some of the path references, such as data:/, temp:/, store:/, private:/ can be set in the startup
% arguments of bcilab.m or bcilab_config.m (see documentation in those files). It is permitted to
% use other relative paths, such as home:/ when specifying those arguments (but circular references
% must be avoided).
%
% In:
%   IndependentPath : platform-independent path; may contain forward and/or backward slashes
%                     (forward slashes generally preferred), and may refer to locations such as
%                     store:/ (the store path) or data:/ (one of the data paths); can also be a
%                     relative path
%
% Out:
%   SystemPath : system-specific path (with slashes corrected and locations resolved)
%                if multiple data paths are present, the one where the minimum number of
%                directories (and files) would have to be created to write to the given file is
%                selected.
%
% Examples:
%   % translate a platform-independent reference to a subdirectory of the data path to one that is
%   % recognized by the operating system (output might be, e.g., 'C:\Projects\mydata\test.mat')
%   env_translatepath('data:/projects/test.mat');
%
%   % resolve a reference to the storage directory
%   env_translatepath('store:/studyXY/result.mat')
%
%   % resolve a reference to the BCILAB root directory
%   env_translatepath('bcilab:/userscripts/myscript.m')
%
%   % resolve a reference to the user's home directory
%   env_translatepath('home:/myconfig.m')
%
%   % resolve a reference to the current temp directory (as specified in the startup options)
%   env_translatepath('temp:/output/001.mat')
%
%   % resolve a reference to the resources directory
%   env_translatepath('resources:/workspaces/testing.mat')
%
% See also:
%   env_startup, io_loadset
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-06-29

global tracking;

if ~ischar(filename)
    error('The file name must be given as a string.'); end
if size(filename,1) ~= 1
    error('The file name must be a non-empty row vector of characters.'); end

% function to sanitize file names (replace file separators by system-dependent version)
sanitize = @(filename) strrep(strrep(strrep(filename,'\',filesep),'/',filesep),[filesep filesep],filesep);

% loop until all path references have been resolved
while true
    % resolve location references
    if strncmp('store:',filename,6)
        filename = [tracking.paths.store_path filename(1+length('store:'):end)];
    elseif strncmp('resources:',filename,10)
        filename = [tracking.paths.resource_path filename(1+length('resources:'):end)];
    elseif strncmp('temp:',filename,5)
        filename = [tracking.paths.temp_path filename(1+length('temp:'):end)];
    elseif strncmp('functions:',filename,10)
        filename = [tracking.paths.function_path filename(1+length('functions:'):end)];
    elseif strncmp('bcilab:',filename,7)
        filename = [tracking.paths.bcilab_path filename(1+length('bcilab:'):end)];
    elseif strncmp('dependencies:',filename,13)
        filename = [tracking.paths.dependency_path filename(1+length('dependencies:'):end)];
    elseif strncmp('private:',filename,8)
        filename = [tracking.paths.private_path filename(1+length('private:'):end)];
    elseif strncmp('home:',filename,5)
        filename = [hlp_homedir filename(1+length('home:'):end)];
    elseif strncmp('data:',filename,5)
        rest = filename(1+length('data:'):end);
        bestpath = 1; bestlen = -1;
        if length(tracking.paths.data_paths) > 1
            % find the data path that contains the longest prefix of the filename
            fpieces = hlp_split(sanitize(rest),filesep);
            for pidx=1:length(tracking.paths.data_paths)
                p = env_translatepath(tracking.paths.data_paths{pidx});
                % for each prefix of the filename (starting with the longest one)
                for k=length(fpieces):-1:0
                    % check if the data path plus the first k pieces of the filename exists
                    if exist([p filesep fpieces{1:k}],'file')
                        % found a match - check if it is a new length record among all our data paths...
                        if k>bestlen
                            bestlen = k;
                            bestpath = pidx;
                        end
                        break;
                    end
                end
            end
        end
        % resolve the reference using that data path which matches most of the filename,
        % where, if multiple data paths are equally well suited, the first one of them is taken
        filename = [tracking.paths.data_paths{bestpath} rest];
    else
        break;
    end
end

filename = sanitize(filename);

repo: REST


text:
function nz_num = wathen_st_size ( nx, ny )

%*****************************************************************************80
%
%% WATHEN_ST_SIZE: Size of Wathen matrix stored in sparse triplet format.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 June 2014
%
%  Author:
%
%    John Burkardt.
%
%  Reference:
%
%    Nicholas Higham,
%    Algorithm 694: A Collection of Test Matrices in MATLAB,
%    ACM Transactions on Mathematical Software,
%    Volume 17, Number 3, September 1991, pages 289-305.
%
%    Andrew Wathen,
%    Realistic eigenvalue bounds for the Galerkin mass matrix,
%    IMA Journal of Numerical Analysis,
%    Volume 7, Number 4, October 1987, pages 449-457.
%
%  Parameters:
%
%    Input, integer NX, NY, values which determine the size of the matrix.
%
%    Output, integer NZ_NUM, the number of items of data used to describe
%    the matrix.
%
  nz_num = nx * ny * 64;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [tg,theta] = tgmo(tmap,ntex,radius,norient,varargin)
% function [tg,theta] = tgmo(tmap,ntex,radius,norient,...)
%
% Compute the texture gradient at a single scale and multiple
% orientations.
%
% INPUT
%	tmap		Texton map, values in [1,ntex].
%	ntex		Number of textons.
%	radius		Radius of disc for texture gradient.
%	norient		Number of orientation at which to compute
%			the texture gradient.
%	'smooth'	Smoothing method, one of
%			{'gaussian','savgol','none'}, default 'none'.
%	'sigma'		Sigma for smoothing, default to radius.
%	'tsim'		Texton similarity matrix.  If not
%			provided, then use chi-squared.
%
% OUTPUT
%	tg		Size [h w norient] array of tg images.
%	theta		Vector of disc orientations (which are
%			orthogonal to the texture gradient).
%
% David R. Martin <dmartin@eecs.berkeley.edu>
% March 2003

% process options
smooth = 'none';
sigma = radius;
usechi2 = true;
for i = 1:2:numel(varargin),
  opt = varargin{i};
  if ~ischar(opt), error('option names not a string'); end
  if i==numel(varargin), error(sprintf('option ''%s'' has no value',opt)); end
  val = varargin{i+1};
  switch opt,
   case 'smooth',
    switch val,
     case {'none','gaussian','savgol'}, smooth=val;
     otherwise, error(sprintf('invalid option smooth=''%s''',val));
    end
   case 'sigma', sigma=val;
   case 'tsim', tsim=val; usechi2=false;
   otherwise, error(sprintf('invalid option ''%s''',opt));
  end
end

radius = max(1,radius);
norient = max(1,norient);
theta = (0:norient-1)/norient*pi;

[h,w] = size(tmap);
tg = zeros(h,w,norient);
fwrite(2,'[');
for i = 1:norient,
  fwrite(2,'.');
  if usechi2,
    tg(:,:,i) = tgso(tmap,ntex,radius,theta(i),...
                     'smooth',smooth,'sigma',sigma);
  else
    tg(:,:,i) = tgso(tmap,ntex,radius,theta(i),...
                     'smooth',smooth,'sigma',sigma,'tsim',tsim);
  end
end
fwrite(2,sprintf(']\n'));

author: Cloud-CV

repo: object-proposals


text:
function tricost(varargin)
%TRICOST draw quality-metrics for a 2-simplex triangulation
%embedded in the two-dimensional plane.
%   TRICOST(VERT,EDGE,TRIA,TNUM) draws histograms of quality
%   metrics for the triangulation.
%   VERT is a V-by-2 array of XY coordinates in the triangu-
%   lation, EDGE is an array of constrained edges, TRIA is a
%   T-by-3 array of triangles, and TNUM is a T-by-1 array of
%   part indices. Each row of TRIA and EDGE define an eleme-
%   nt. VERT(TRIA(II,1),:), VERT(TRIA(II,2),:) and VERT(TRIA
%   (II,3),:) are the coordinates of the II-TH triangle. The
%   edges in EDGE are defined in a similar manner. NUM is an
%   array of part indexing, such that TNUM(II) is the index
%   of the part in which the II-TH triangle resides.
%
%   TRICOST(...,HVRT) additionally draws histograms of rela-
%   tive edge-length, indicating conformance to the spacing
%   constraints. HVRT is a V-by-1 array of spacing informat-
%   ion, per an evaluation of the mesh-size function at the
%   mesh vertices VERT.
%
%   See also REFINE2, SMOOTH2

%-----------------------------------------------------------
%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 09/07/2018
%-----------------------------------------------------------

    vert = [] ; conn = [] ; tria = [] ;
    tnum = [] ; hvrt = [] ;

%---------------------------------------------- extract args
    if (nargin>=+1), vert = varargin{1}; end
    if (nargin>=+2), conn = varargin{2}; end
    if (nargin>=+3), tria = varargin{3}; end
    if (nargin>=+4), tnum = varargin{4}; end
    if (nargin>=+5), hvrt = varargin{5}; end

%---------------------------------------------- basic checks
    if ( ~isnumeric(vert) || ...
         ~isnumeric(conn) || ...
         ~isnumeric(tria) || ...
         ~isnumeric(tnum) || ...
         ~isnumeric(hvrt) )
        error('tricost:incorrectInputClass' , ...
            'Incorrect input class.') ;
    end

%---------------------------------------------- basic checks
    if (ndims(vert) ~= +2 || ...
        ndims(conn) ~= +2 || ...
        ndims(tria) ~= +2 )
        error('tricost:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end
    if (size(vert,2)~= +2 || ...
        size(conn,2) < +2 || ...
        size(tria,2) < +3 )
        error('tricost:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

    nvrt = size(vert,1) ;
    ntri = size(tria,1) ;

%---------------------------------------------- basic checks
    if (min(min(conn(:,1:2))) < +1 || ...
            max(max(conn(:,1:2))) > nvrt )
        error('tricost:invalidInputs', ...
            'Invalid EDGE input array.') ;
    end

    if (min(min(tria(:,1:3))) < +1 || ...
            max(max(tria(:,1:3))) > nvrt )
        error('tricost:invalidInputs', ...
            'Invalid TRIA input array.') ;
    end

%-- borrowed from the JIGSAW library!

%-- draw sub-axes directly -- sub-plot gives
%-- silly inconsistent ax spacing...!

    axpos31 = [.125,.750,.800,.150] ;
    axpos32 = [.125,.450,.800,.150] ;
    axpos33 = [.125,.150,.800,.150] ;

    axpos41 = [.125,.835,.800,.135] ;
    axpos42 = [.125,.590,.800,.135] ;
    axpos43 = [.125,.345,.800,.135] ;
    axpos44 = [.125,.100,.800,.135] ;

%-- draw cost histograms for 2-tria elements
    figure;
    set(gcf,'color','w','units','normalized', ...
        'position',[.05,.10,.30,.30]);
    if (~isempty(hvrt))

%-- have size-func data
    axes('position',axpos41); hold on;
    scrhist(triscr2(vert,tria),'tria3');
    axes('position',axpos42); hold on;
    anghist(triang2(vert,tria),'tria3');
    axes('position',axpos43); hold on;
    hfnhist(relhfn2(vert, ...
                    tria,hvrt),'tria3');
    axes('position',axpos44); hold on;
    deghist(trideg2(vert,tria),'tria3');

    else

%-- null size-func data
    axes('position',axpos31); hold on;
    scrhist(triscr2(vert,tria),'tria3');
    axes('position',axpos32); hold on;
    anghist(triang2(vert,tria),'tria3');
    axes('position',axpos33); hold on;
    deghist(trideg2(vert,tria),'tria3');

    end

end

function [mf] = mad(ff)
%MAD return mean absolute deviation (from the mean).

    mf = mean(abs(ff-mean(ff))) ;

end

function deghist(dd,ty)
%DEGHIST draw histogram for "degree" quality-metric.

    dd = dd(:);
    be = 1:max(dd);
    hc = histc(dd,be);

    r = [.85,.00,.00] ; y = [1.0,.95,.00] ;
    g = [.00,.90,.00] ; k = [.60,.60,.60] ;

    bar(be,hc,1.05,'facecolor',k,'edgecolor',k);

    axis tight;
    set(gca,'ycolor', get(gca,'color'),'ytick',[],...
        'xtick',2:2:12,'layer','top','fontsize',...
            14,'linewidth',2.,'ticklength',[.025,.025],...
                'box','off','xlim',[0,12]);

    switch (ty)
    case 'tria4'

    if ( ~(exist('OCTAVE_VERSION','builtin') > +0) )
        text(-.225,0,'$|d|_{\tau}$',...
            'horizontalalignment','right',...
                'fontsize',22,'interpreter','latex') ;
    else
        text(-.225,0, '|d|_{\tau}' ,...
            'horizontalalignment','right',...
                'fontsize',22,'interpreter',  'tex') ;
    end

    case 'tria3'

    if ( ~(exist('OCTAVE_VERSION','builtin') > +0) )
        text(-.225,0,'$|d|_{f}$',...
            'horizontalalignment','right',...
                'fontsize',22,'interpreter','latex') ;
    else
        text(-.225,0, '|d|_{\tau}' ,...
            'horizontalalignment','right',...
                'fontsize',22,'interpreter',  'tex') ;
    end

    end

end

function anghist(ad,ty)
%ANGHIST draw histogram for "angle" quality-metric.

    ad = ad(:);
    be = linspace(0.,180.,91);
    bm =(be(1:end-1)+be(2:end))/2.;
    hc = histc(ad,be);

    switch (ty)
    case 'tria4'
        poor = bm <  10.  | bm >= 160. ;
        okay =(bm >= 10.  & bm <  20. )| ...
              (bm >= 140. & bm <  160.);
        good =(bm >= 20.  & bm <  30. )| ...
              (bm >= 120. & bm <  140.);
        best = bm >= 30.  & bm <  120. ;

    case 'tria3'
        poor = bm <  15.  | bm >= 150. ;
        okay =(bm >= 15.  & bm <  30. )| ...
              (bm >= 120. & bm <  150.);
        good =(bm >= 30.  & bm <  45. )| ...
              (bm >= 90.  & bm <  120.);
        best = bm >= 45.  & bm <  90.  ;

    end

    r = [.85,.00,.00] ; y = [1.0,.95,.00] ;
    g = [.00,.90,.00] ; k = [.60,.60,.60] ;

    bar(bm(poor),hc(poor),1.05,...
        'facecolor',r,'edgecolor',r) ;
    bar(bm(okay),hc(okay),1.05,...
        'facecolor',y,'edgecolor',y) ;
    bar(bm(good),hc(good),1.05,...
        'facecolor',g,'edgecolor',g) ;
    bar(bm(best),hc(best),1.05,...
        'facecolor',k,'edgecolor',k) ;

    axis tight;
    set(gca,'ycolor', get(gca,'color'),'ytick',[],...
        'xtick',0:30:180,'layer','top','fontsize',...
            14,'linewidth',2.,'ticklength',[.025,.025],...
                'box','off','xlim',[0.,180.]) ;

    mina = max(1.000,min(ad)); %%!! so that axes don't obscure!
    maxa = min(179.0,max(ad));

    bara = mean(ad(:));
    mada = mad (ad(:));

    line([ mina, mina],...
        [0,max(hc)],'color','r','linewidth',1.5);
    line([ maxa, maxa],...
        [0,max(hc)],'color','r','linewidth',1.5);

    if ( mina > 25.0)
        text(mina-1.8,.90*max(hc),num2str(min(ad),'%16.1f'),...
            'horizontalalignment',...
                'right','fontsize',15) ;
    else
        text(mina+1.8,.90*max(hc),num2str(min(ad),'%16.1f'),...
            'horizontalalignment',...
                'left' ,'fontsize',15) ;
    end

    if ( maxa < 140.)
        text(maxa+1.8,.90*max(hc),num2str(max(ad),'%16.1f'),...
            'horizontalalignment',...
                'left' ,'fontsize',15) ;
    else
        text(maxa-1.8,.90*max(hc),num2str(max(ad),'%16.1f'),...
            'horizontalalignment',...
                'right','fontsize',15) ;
    end

    if ( maxa < 100.)

    if ( ~(exist('OCTAVE_VERSION','builtin') > +0) )
        text(maxa-16.,.45*max(hc),...
        '$\bar{\sigma}_{\theta}\!= $',...
            'horizontalalignment', 'left',...
                'fontsize',16,'interpreter','latex') ;

        text(maxa+1.8,.45*max(hc),num2str(mad(ad),'%16.2f'),...
            'horizontalalignment',...
                'left' ,'fontsize',15) ;
    end

    else

    if ( ~(exist('OCTAVE_VERSION','builtin') > +0) )
        text(maxa-16.,.45*max(hc),...
        '$\bar{\sigma}_{\theta}\!= $',...
            'horizontalalignment', 'left',...
                'fontsize',16,'interpreter','latex') ;

        text(maxa+1.8,.45*max(hc),num2str(mad(ad),'%16.3f'),...
            'horizontalalignment',...
                'left' ,'fontsize',15) ;
    end

    end

    switch (ty)
    case 'tria4'

    if ( ~(exist('OCTAVE_VERSION','builtin') > +0) )
        text(-9.0,0.0,'$\theta_{\tau}$',...
            'horizontalalignment','right',...
                'fontsize',22,'interpreter','latex') ;
    else
        text(-9.0,0.0, '\theta_{\tau}' ,...
            'horizontalalignment','right',...
                'fontsize',22,'interpreter',  'tex') ;
    end

    case 'tria3'

    if ( ~(exist('OCTAVE_VERSION','builtin') > +0) )
        text(-9.0,0.0,'$\theta_{f}$',...
            'horizontalalignment','right',...
                'fontsize',22,'interpreter','latex') ;
    else
        text(-9.0,0.0, '\theta_{f}' ,...
            'horizontalalignment','right',...
                'fontsize',22,'interpreter',  'tex') ;
    end

    end

end

function scrhist(sc,ty)
%SCRHIST draw histogram for "score" quality-metric.

    sc = sc(:);
    be = linspace(0.,1.,101);
    bm = (be(1:end-1)+be(2:end)) / 2.;
    hc = histc(sc,be);

    switch (ty)
    case{'tria4','dual4'}
        poor = bm <  .25 ;
        okay = bm >= .25 & bm <  .50 ;
        good = bm >= .50 & bm <  .75 ;
        best = bm >= .75 ;

    case{'tria3','dual3'}
        poor = bm <  .30 ;
        okay = bm >= .30 & bm <  .60 ;
        good = bm >= .60 & bm <  .90 ;
        best = bm >= .90 ;

    end

    r = [.85,.00,.00] ; y = [1.0,.95,.00] ;
    g = [.00,.90,.00] ; k = [.60,.60,.60] ;

    bar(bm(poor),hc(poor),1.05,...
        'facecolor',r,'edgecolor',r) ;
    bar(bm(okay),hc(okay),1.05,...
        'facecolor',y,'edgecolor',y) ;
    bar(bm(good),hc(good),1.05,...
        'facecolor',g,'edgecolor',g) ;
    bar(bm(best),hc(best),1.05,...
        'facecolor',k,'edgecolor',k) ;

    axis tight;
    set(gca,'ycolor', get(gca,'color'),'ytick',[],...
        'xtick',.0:.2:1.,'layer','top','fontsize',...
            14,'linewidth',2.,'ticklength',[.025,.025],...
                'box','off','xlim',[0.,1.]) ;

    mins = max(0.010,min(sc)); %%!! so that axes don't obscure!
    maxs = min(0.990,max(sc));

    line([ mins, mins],...
        [0,max(hc)],'color','r','linewidth',1.5);
    line([mean(sc),mean(sc)],...
        [0,max(hc)],'color','r','linewidth',1.5);

    if ( mins > .4)
        text(mins-.01,.9*max(hc),num2str(min(sc),'%16.3f'),...
            'horizontalalignment',...
                'right','fontsize',15) ;
    else
        text(mins+.01,.9*max(hc),num2str(min(sc),'%16.3f'),...
            'horizontalalignment',...
                'left' ,'fontsize',15) ;
    end

    if ( mean(sc) > mins + .150)
    text(mean(sc)-.01,.9*max(hc),num2str(mean(sc),'%16.3f'),...
        'horizontalalignment','right','fontsize',15) ;
    end

    switch (ty)
    case 'tria4'

    if ( ~(exist('OCTAVE_VERSION','builtin') > +0) )
        text(-.04,0.0, ...
        '$\mathcal{Q}^{\mathcal{T}}_{\tau}$',...
            'horizontalalignment','right',...
                'fontsize',22,'interpreter','latex') ;
    else
        text(-.04,0.0,'Q^{t}_{\tau}',...
            'horizontalalignment','right',...
                'fontsize',22,'interpreter',  'tex') ;
    end

    case 'tria3'

    if ( ~(exist('OCTAVE_VERSION','builtin') > +0) )
        text(-.04,0.0, ...
        '$\mathcal{Q}^{\mathcal{T}}_{f}$',...
            'horizontalalignment','right',...
                'fontsize',22,'interpreter','latex') ;
    else
        text(-.04,0.0,'Q^{t}_{f}',...
            'horizontalalignment','right',...
                'fontsize',22,'interpreter',  'tex') ;
    end

    case 'dual4'

    if ( ~(exist('OCTAVE_VERSION','builtin') > +0) )
        text(-.04,0.0, ...
        '$\mathcal{Q}^{\mathcal{D}}_{\tau}$',...
            'horizontalalignment','right',...
                'fontsize',22,'interpreter','latex') ;
    else
        text(-.04,0.0,'Q^{d}_{\tau}',...
            'horizontalalignment','right',...
                'fontsize',22,'interpreter',  'tex') ;
    end

    case 'dual3'

    if ( ~(exist('OCTAVE_VERSION','builtin') > +0) )
        text(-.04,0.0, ...
        '$\mathcal{Q}^{\mathcal{D}}_{f}$',...
            'horizontalalignment','right',...
                'fontsize',22,'interpreter','latex') ;
    else
        text(-.04,0.0,'Q^{d}_{f}',...
            'horizontalalignment','right',...
                'fontsize',22,'interpreter',  'tex') ;
    end

    end

end

function hfnhist(hf,ty)
%HFNHIST draw histogram for "hfunc" quality-metric.

    be = linspace(0.,2.,101);
    bm = (be(1:end-1)+be(2:end)) / 2.;
    hc = histc(hf,be);

    poor = bm <  .40 | bm >= 1.6  ;
    okay =(bm >= .40 & bm <  .60 )| ...
          (bm >= 1.4 & bm <  1.6 );
    good =(bm >= .60 & bm <  .80 )| ...
          (bm >= 1.2 & bm <  1.4 );
    best = bm >= .80 & bm <  1.2 ;

    r = [.85,.00,.00] ; y = [1.0,.95,.00] ;
    g = [.00,.90,.00] ; k = [.60,.60,.60] ;

    bar(bm(poor),hc(poor),1.05,...
        'facecolor',r,'edgecolor',r) ;
    bar(bm(okay),hc(okay),1.05,...
        'facecolor',y,'edgecolor',y) ;
    bar(bm(good),hc(good),1.05,...
        'facecolor',g,'edgecolor',g) ;
    bar(bm(best),hc(best),1.05,...
        'facecolor',k,'edgecolor',k) ;

    axis tight;
    set(gca,'ycolor', get(gca,'color'),'ytick',[],...
        'xtick',.0:.5:2.,'layer','top','fontsize',...
            14,'linewidth',2.,'ticklength',[.025,.025],...
                'box','off','xlim',[0.,2.]);

    line([max(hf),max(hf)],...
        [0,max(hc)],'color','r','linewidth',1.5);

    text(max(hf)+.02,.90*max(hc),num2str(max(hf),'%16.2f'),...
        'horizontalalignment','left','fontsize',15) ;

    if ( ~(exist('OCTAVE_VERSION','builtin') > +0) )

    text(max(hf)-.18,.45*max(hc),'$\bar{\sigma}_{h}\! = $',...
        'horizontalalignment','left',...
            'fontsize',16,'interpreter','latex') ;

    text(max(hf)+.02,.45*max(hc),num2str(mad(hf),'%16.2f'),...
        'horizontalalignment','left','fontsize',15) ;

    end

    if ( ~(exist('OCTAVE_VERSION','builtin') > +0) )

    text(-0.10,0.0,'$h_{r}$','horizontalalignment','right',...
        'fontsize',22,'interpreter','latex') ;

    else

    text(-0.10,0.0, 'h_{r}' ,'horizontalalignment','right',...
        'fontsize',22,'interpreter',  'tex') ;

    end

end

author: dengwirda


text:
function gdif = p00_gdif ( problem, n, x )

%*****************************************************************************80
%
%% P00_GDIF approximates the gradient via finite differences.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the problem number.
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the point where the gradient
%    is to be approximated.
%
%    Output, real GDIF(N), the approximated gradient vector.
%
  tol = eps^0.33;

  for i = 1 : n

    if ( 0.0 <= x(i) )
      dx = eps * ( x(i) + 1.0 );
    else
      dx = eps * ( x(i) - 1.0 );
    end

    xi = x(i);
    x(i) = xi + dx;
    fplus = p00_f ( problem, n, x );

    x(i) = xi - dx;
    fminus = p00_f ( problem, n, x );

    gdif(i) = ( fplus - fminus ) / ( 2.0 * dx );

    x(i) = xi;

  end

  return
end


text:
%------------------------------ PolyMesher -------------------------------%
% Ref: C Talischi, GH Paulino, A Pereira, IFM Menezes, "PolyMesher: A     %
%      general-purpose mesh generator for polygonal elements written in   %
%      Matlab," Struct Multidisc Optim, DOI 10.1007/s00158-011-0706-z     %
%-------------------------------------------------------------------------%
function d = dLine(P,x1,y1,x2,y2)
% By convention, a point located at the left hand side of the line
% is inside the region and it is assigned a negative distance value.
a = [x2-x1,y2-y1]; a = a/norm(a);
b = [P(:,1)-x1,P(:,2)-y1];
d = b(:,1)*a(2) - b(:,2)*a(1);
d = [d,d];
%-------------------------------------------------------------------------%

repo: ifem


text:
classdef PitchTermCondition < AbstractEventTerminationCondition
    %BankAngleTermCondition Summary of this class goes here
    %   Detailed explanation goes here

    properties
        steeringModel(1,1) AbstractSteeringModel = RollPitchYawPolySteeringModel.getDefaultSteeringModel();
        targetPitchAngle(1,1) double = 0;
        bodyInfo KSPTOT_BodyInfo
    end

    methods
        function obj = PitchTermCondition(targetPitchAngle)
            obj.targetPitchAngle = targetPitchAngle;
        end

        function evtTermCondFcnHndl = getEventTermCondFuncHandle(obj)
            evtTermCondFcnHndl = @(t,y) obj.eventTermCond(t,y, obj.targetPitchAngle, obj.steeringModel, obj.bodyInfo);
        end

        function initTermCondition(obj, initialStateLogEntry)
            obj.steeringModel = initialStateLogEntry.steeringModel;
            obj.bodyInfo = initialStateLogEntry.centralBody;
        end

        function name = getName(obj)
            name = sprintf('Pitch Angle (%.3f deg)', rad2deg(obj.targetPitchAngle));
        end

        function tf = shouldBeReinitOnRestart(obj)
            tf = false;
        end

        function params = getTermCondUiStruct(obj)
            params = struct();

            params.paramName = 'Target Pitch Angle';
            params.paramUnit = 'deg';
            params.useParam = 'on';
            params.useStages = 'off';
            params.useTanks = 'off';
            params.useEngines = 'off';
            params.useStopwatches = 'off';

            params.value = rad2deg(obj.targetPitchAngle);
            params.refStage = LaunchVehicleStage.empty(1,0);
            params.refTank = LaunchVehicleEngine.empty(1,0);
            params.refEngine = LaunchVehicleEngine.empty(1,0);
            params.refStopwatch = LaunchVehicleStopwatch.empty(1,0);
        end

        function optVar = getNewOptVar(obj)
            optVar = PitchAngleTermCondOptimVar(obj);
        end

        function optVar = getExistingOptVar(obj)
            optVar = obj.optVar;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end
    end

    methods(Static)
        function termCond = getTermCondForParams(paramValue, stage, tank, engine, stopwatch)
            termCond = PitchTermCondition((paramValue));
        end
    end

    methods(Static, Access=private)
        function [value,isterminal,direction] = eventTermCond(t,y, targetPitchAngle, steeringModel, bodyInfo)
            ut = t;
            rVect = y(1:3);
            vVect = y(4:6);

            dcm = steeringModel.getBody2InertialDcmAtTime(ut, rVect, vVect, bodyInfo);
            [~,pitchAngle,~] = computeEulerAnglesFromInertialBodyAxes(ut, rVect, vVect, bodyInfo, dcm(:,1), dcm(:,2), dcm(:,3));

            value = pitchAngle - targetPitchAngle;
            isterminal = 1;
            direction = 0;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%EUL2TR Convert Euler angles to homogeneous transform
%
% T = EUL2TR(PHI, THETA, PSI, OPTIONS) is a SE(3) homogeneous
% transformation matrix (4x4) equivalent to the specified Euler angles.
% These correspond to rotations about the Z, Y, Z axes respectively. If
% PHI, THETA, PSI are column vectors (Nx1) then they are assumed to
% represent a trajectory and R is a three-dimensional matrix (4x4xN), where
% the last index corresponds to rows of PHI, THETA, PSI.
%
% T = EUL2TR(EUL, OPTIONS) as above but the Euler angles are taken from
% consecutive columns of the passed matrix EUL = [PHI THETA PSI].  If EUL
% is a matrix (Nx3) then they are assumed to represent a trajectory and T
% is a three-dimensional matrix (4x4xN), where the last index corresponds
% to rows of EUL which are assumed to be [PHI, THETA, PSI].
%
% Options::
%  'deg'      Compute angles in degrees (radians default)
%
% Note::
% - The vectors PHI, THETA, PSI must be of the same length.
% - The translational part is zero.
%
% See also EUL2R, RPY2TR, TR2EUL.

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

function T = eul2tr(phi, varargin)

    R = eul2r(phi, varargin{:});
    T = r2t(R);

author: Allopart

repo: rbpf-gmapping


text:
function [ops, stat, res] = fast_clustering_with_neuropil(ops, U, Sv)

U =  reshape(U, [], size(U,ndims(U)));
iplane = ops.iplane;

for i = 1:size(U,2)
  U(:,i) = U(:,i)  * Sv(i).^.5;
end
U = U';
[nSVD, Npix] = size(U);
%% clustering options
Ly = numel(ops.yrange);
Lx = numel(ops.xrange);

ops.Nk0 = floor(sqrt(ops.Nk0))^2;

Nk      = ops.Nk0;
nsqrt = round(sqrt(Nk));

xs = repmat(round(linspace(1, nsqrt, Lx)), Ly, 1);
ys = repmat(round(linspace(1, nsqrt, Ly))', 1, Lx);
iclust = xs + (ys-1) * nsqrt;

clear xs ys

niter   = ops.niterclustering;

% xs = repmat(1:Lx, Ly, 1);
% ys = repmat((1:Ly)', 1, Lx);
%
% randx = rand(1, Nk) * Lx;
% randy = rand(1, Nk) * Ly;
%
% dx = repmat(xs(:), 1, Nk) - repmat(randx, numel(xs(:)), 1);
% dy = repmat(ys(:), 1, Nk) - repmat(randy, numel(ys(:)), 1);
%
% dxy = dx.^2 + dy.^2;
% [~, iclust] = min(dxy, [], 2);
%%
% clear dx dy

if ops.ShowCellMap
    figure( 'Units', 'pixels', 'position', [100 100 900 900])
    colormap('hsv')
    axes('position', [.05 .05 .925 .925])
    set(gcf, 'Color', 'w')
end

r   = rand(1, Nk);
Sat = ones(Ly, Lx);

err = zeros(niter,1);
ops.meanV = gather(sum(Sv)/(Ly*Lx));

Nk = ops.Nk0;
Nkiter = round(linspace(ops.Nk0, ops.Nk, niter-2));
Nkiter(end+1:(niter+1)) = ops.Nk;

%%
M = .0001 * ones(1, Npix, 'single');

ison = true(Nk,1);
%
nBasis = 10;

xc = linspace(1, Lx, nBasis);
yc = linspace(1, Ly, nBasis);
yc = yc';
xs = 1:Lx;
ys = 1:Ly;

sigx = 4*(Lx - 1)/nBasis;
sigy = 4*(Ly - 1)/nBasis;

S = zeros(Ly, Lx, nBasis, nBasis, 'single');
for kx = 1:nBasis
    for ky = 1:nBasis
        cosx = 1+cos(2*pi*(xs - xc(kx))/sigx);
        cosy = 1+cos(2*pi*(ys - yc(ky))/sigy);
        cosx(abs(xs-xc(kx))>sigx/2) = 0;
        cosy(abs(ys-yc(ky))>sigy/2) = 0;

        S(:, :,ky, kx) = cosy' * cosx;
    end
end
S = reshape(S, [], nBasis^2);
S = normc(S);

%%
% nBasis = 0;
% S = zeros(Npix, nBasis^2);

StS = S' * S;
StU = S' * U';
LtL = zeros(Nk, Nk, 'single');
LtU = zeros(Nk, nSVD, 'single');
LtS = zeros(Nk, nBasis^2, 'single');
Ireg = diag([ones(Nk,1); zeros(nBasis^2,1)]);

tic
for k = 1:niter
    for i = 1:Nk
        ix = find(iclust==i);
        if numel(ix)==0
            LtU(i,:) = 0;
            LtL(i,i) = 0;
            LtS(i,:) = 0;
        else
            LtU(i,:) = M(ix) * U(:, ix)';
            LtL(i,i) = sum(M(ix).^2);
            LtS(i,:) = M(ix) * S(ix, :);
        end
    end
    %
    covL = [LtL LtS; LtS' StS];
    LtXS = [LtU; StU];

    Lam = (covL + 1e-4 * Ireg) \ LtXS;

    vs = normc(Lam(1:Nk, :)');
    neuropil = Lam(Nk + [1:nBasis^2], :)' * S';

    xs = vs' * (U - neuropil);

    [M, iclust] = max(xs,[],1);
    if 1
        %---------------------------------------------%
        xs(iclust + (0:Nk:numel(xs)-1)) = 0;
        [M2, iclust2] = max(abs(xs),[],1);

        dM = M - M2;
        indrem = Nkiter(k) - Nkiter(k+1);
        dMk = 0*ones(Nk,1);
        icl = cell(Nk,1);
        for j = 1:Nk
            icl{j} = iclust==j;
            if ~isempty(icl{j})
                dMk(j) = sum(dM(icl{j}));
            end
        end

        vlk = zeros(Nk, 1);
        vlk(~ison) = Inf;
        while indrem>0
            [Xmin, imin]          = min(dMk + vlk);
            if isinf(Xmin)
                Nkiter(k+1) = sum(ison);
                break;
            end
            newi               = iclust2(icl{imin});
            iclust(icl{imin})  = newi;
            M(icl{imin})       = M2(icl{imin});
            dMk(unique(newi))  = Inf;
            dMk(imin)          = Inf;
            dMk(unique(iclust(iclust2==imin))) ...
                = Inf;

            ison(imin) = 0;
            indrem             = indrem - 1;
        end
    end
    %---------------------------------------------%

    err(k) = sum(M(:));

    if (rem(k,10)==1 || k==niter) && ops.ShowCellMap
        %%
        lam = M;
        for i = 1:Nk
            ix = find(iclust==i);
            nT0 = numel(ix);
            if nT0>0
                vM = lam(ix);
%                 vM = vM/sum(vM.^2)^.5;
                lam(ix) = vM;
            end
        end
%         V = max(0, min(10 * reshape(lam, Ly, Lx), 1));
        V = max(0, min(.5 * reshape(lam, Ly, Lx)/mean(lam(:)), 1));
        H = reshape(r(iclust), Ly, Lx);
        rgb_image = hsv2rgb(cat(3, H, Sat, V));
        imagesc(rgb_image)
        axis off
        drawnow
        fprintf('explained variance is %2.6f time %2.2f \n', err(k), toc)
    end

end

lam = M;
for i = 1:Nk
    ix = find(iclust==i);

    nT0 = numel(ix);
    if nT0>0
        vM = lam(ix);
        lam(ix) = vM/sum(vM.^2)^.5;
    end
end

%%
newindx = cumsum(ison);
iclust  = newindx(iclust);
Nk      = numel(unique(iclust));
%
clear res

res.iclust  = iclust;
res.M       = M;
res.S       = S;
res.covL    = covL + 1e-4 * Ireg;
res.lambda  = lam;

%
res.Ly  = Ly;
res.Lx  = Lx;
stat    = get_stat(res);

if ~exist(ops.ResultsSavePath, 'dir')
    mkdir(ops.ResultsSavePath)
end
save(sprintf('%s/F_%s_%s_plane%d_Nk%d.mat', ops.ResultsSavePath, ...
    ops.mouse_name, ops.date, iplane, Nk),  'ops', 'res', 'stat')

%%
% sk = skewness(F,[],2);
% [~, isk] = sort(sk, 'descend');
% clf
% for i = 1:20
%    plot(5*i + zscore(F(isk(i), :)))
%    hold all
% end
% axis tight

author: cortex-lab


text:
function [ha] = yd22ha(yd2)
% Convert area from square yards to hectares.
% Chad A. Greene 2012
ha = yd2*0.000083612736;


text:
function Y = pdist(X,s,varargin)
%PDIST Pairwise distance between observations.
%   Y = PDIST(X) returns a vector Y containing the Euclidean distances
%   between each pair of observations in the M-by-N data matrix X.  Rows of
%   X correspond to observations, columns correspond to variables.  Y is an
%   (M*(M-1)/2)-by-1 vector, corresponding to the M*(M-1)/2 pairs of
%   observations in X.
%
%   Y = PDIST(X, DISTANCE) computes Y using DISTANCE.  Choices are:
%
%       'euclidean'   - Euclidean distance
%       'seuclidean'  - Standardized Euclidean distance, each coordinate
%                       in the sum of squares is inverse weighted by the
%                       sample variance of that coordinate
%       'cityblock'   - City Block distance
%       'mahalanobis' - Mahalanobis distance
%       'minkowski'   - Minkowski distance with exponent 2
%       'cosine'      - One minus the cosine of the included angle
%                       between observations (treated as vectors)
%       'correlation' - One minus the sample correlation between
%                       observatons (treated as sequences of values).
%       'hamming'     - Hamming distance, percentage of coordinates
%                       that differ
%       'jaccard'     - One minus the Jaccard coefficient, the
%                       percentage of nonzero coordinates that differ
%       function      - A distance function specified using @, for
%                       example @DISTFUN
%
%   A distance function must be of the form
%
%         function D = DISTFUN(XI, XJ, P1, P2, ...),
%
%   taking as arguments two L-by-N matrices XI and XJ each of which
%   contains rows of X, plus zero or more additional problem-dependent
%   arguments P1, P2, ..., and returning an L-by-1 vector of distances D,
%   whose Kth element is the distance between the observations XI(K,:)
%   and XJ(K,:).
%
%   Y = PDIST(X, DISTFUN, P1, P2, ...) passes the arguments P1, P2, ...
%   directly to the function DISTFUN.
%
%   Y = PDIST(X, 'minkowski', P) computes Minkowski distance using the
%   positive scalar exponent P.
%
%   The output Y is arranged in the order of ((1,2),(1,3),..., (1,M),
%   (2,3),...(2,M),.....(M-1,M)), i.e. the upper right triangle of the full
%   M-by-M distance matrix.  To get the distance between the Ith and Jth
%   observations (I < J), either use the formula Y((I-1)*(M-I/2)+J-I), or
%   use the helper function Z = SQUAREFORM(Y), which returns an M-by-M
%   square symmetric matrix, with the (I,J) entry equal to distance between
%   observation I and observation J.
%
%   Example:
%
%      X = randn(100, 5);                 % some random points
%      Y = pdist(X, 'euclidean');         % unweighted distance
%      Wgts = [.1 .3 .3 .2 .1];           % coordinate weights
%      Ywgt = pdist(X, @weucldist, Wgts); % weighted distance
%
%      function d = weucldist(XI, XJ, W) % weighted euclidean distance
%      d = sqrt((XI-XJ).^2 * W');
%
%   See also SQUAREFORM, LINKAGE, SILHOUETTE.

%   An example of distance for data with missing elements:
%
%      X = randn(100, 5);     % some random points
%      X(unidrnd(prod(size(X)),1,20)) = NaN; % scatter in some NaNs
%      D = pdist(X, @naneucdist);
%
%      function d = naneucdist(XI, XJ) % euclidean distance, ignoring NaNs
%      sqdx = (XI-XJ).^2;
%      nk = sum(~isnan(sqdx),2); nk(nk == 0) = NaN;
%      d = sqrt(nansum(sqdx')' .* size(XI,2) ./ nk); %correct for missing coords

%   Copyright 1993-2002 The MathWorks, Inc.
%   $Revision: 1.15 $

if nargin < 2
    s = 'euc';
    distfun = @distcalc;
    distargs = {s};
else
    if ischar(s)
        methods = strvcat('euclidean','seuclidean','cityblock','mahalanobis','minkowski','cosine','correlation','hamming','jaccard');
        i = strmatch(lower(s), methods);
        if length(i) > 1
            error(sprintf('Ambiguous ''DISTANCE'' argument:  %s.', s));
        elseif isempty(i)
            % error(sprintf('Unknown ''DISTANCE'' argument:  %s.', s));
            distfun = str2func(s);
            distargs = varargin;
            s = 'usr';
        else
            s = lower(methods(i,1:3));
            distfun = @distcalc;
            distargs = {s};
        end
    elseif isa(s, 'function_handle') |  isa(s, 'inline')
        distfun = s;
        distargs = varargin;
        s = 'usr';
    else
        error('The ''DISTANCE'' argument must be a string or a function.');
    end
end

[m, n] = size(X);
if any(imag(X(:))) & ~isequal(s,'usr')
   error('PDIST does not accept complex inputs.');
end

switch s
case 'seu' % Standardized Euclidean weights by coordinate variance
   distargs{end+1} = 1 ./ var(X)';
case 'mah' % Mahalanobis
   distargs{end+1} = inv(cov(X));
case 'min' % Minkowski distance needs a third argument
   if nargin < 3  % use default value for exponent
      distargs{end+1} = 2;
   elseif varargin{1} > 0
      distargs{end+1} = varargin{1}; % get exponent from input args
   else
      error('The exponent for the Minkowski metric must be positive.');
   end
case 'cos' % Cosine
   Xnorm = sqrt(sum(X.^2, 2));
   if min(Xnorm) <= eps * max(Xnorm)
       error(['Some points have small relative magnitudes, making them ', ...
              'effectively zero.\nEither remove those points, or choose a ', ...
              'distance other than cosine.'], []);
   end
   X = X ./ Xnorm(:,ones(1,n));
case 'cor' % Correlation
   X = X - repmat(mean(X,2),1,n);
   Xnorm = sqrt(sum(X.^2, 2));
   if min(Xnorm) <= eps * max(Xnorm)
       error(['Some points have small relative standard deviations, making ', ...
              'them effectively constant.\nEither remove those points, or ', ...
              'choose a distance other than correlation.'], []);
   end
   X = X ./ Xnorm(:,ones(1,n));
end

if m < 2
   % Degenerate case, just return an empty of the proper size.
   Y = zeros(1,0);
   return;
end

% Create (I,J) defining all pairs of points
p = (m-1):-1:2;
I = zeros(m*(m-1)/2,1);
I(cumsum([1 p])) = 1;
I = cumsum(I);
J = ones(m*(m-1)/2,1);
J(cumsum(p)+1) = 2-p;
J(1)=2;
J = cumsum(J);

% For large matrices, process blocks of rows as a group
n = length(I);
ncols = size(X,2);
blocksize = 1e4;                     % # of doubles to process as a group
M = max(1,ceil(blocksize/ncols));    % # of rows to process as a group
nrem = rem(n,M);
if nrem==0, nrem = min(M,n); end

Y = zeros(1,n);
ii = 1:nrem;
try
    Y(ii) = feval(distfun,X(I(ii),:),X(J(ii),:),distargs{:})';
catch
    if isa(distfun, 'inline')
        error(['The inline distance function generated the following ', ...
               'error:\n%s'], lasterr);
    elseif strfind(lasterr, ...
                   sprintf('Undefined function ''%s''', func2str(distfun)))
        error('The distance function ''%s'' was not found.', func2str(distfun));
    else
        error(['The distance function ''%s'' generated the following ', ...
               'error:\n%s'], func2str(distfun),lasterr);
    end
end;
for j=nrem+1:M:n
    ii = j:j+M-1;
    try
        Y(ii) = feval(distfun,X(I(ii),:),X(J(ii),:),distargs{:})';
    catch
        if isa(distfun, 'inline')
            error(['The inline distance function generated the following ', ...
                    'error:\n%s'], lasterr);
        else
            error(['The distance function ''%s'' generated the following', ...
                    'error:\n%s'], func2str(distfun),lasterr);
        end;
    end;
end

% ----------------------------------------------
function d = distcalc(XI,XJ,s,arg)
%DISTCALC Perform distance calculation for PDIST.
switch s
case 'euc',   d = sqrt(sum((XI-XJ).^2,2));            % Euclidean
case 'seu',   d = sqrt(((XI-XJ).^2) * arg);           % Standardized Euclidean
case 'cit',   d = sum(abs((XI-XJ)),2);                % City Block
case 'mah',   Y = XI - XJ;
              d = sqrt(sum((Y*arg).*Y,2));            % Mahalanobis
case 'min',   d = sum(abs((XI-XJ)).^arg,2).^(1/arg);  % Minkowski
case 'cos',   d = 1 - sum(XI.*XJ,2);                  % Cosine
case 'cor',   d = 1 - sum(XI.*XJ,2);                  % Correlation
case 'ham',   d = sum(XI ~= XJ,2) / size(XI,2);       % Hamming
case 'jac',   nz = XI ~= 0 | XJ ~= 0;
              ne = XI ~= XJ;
              d = sum(ne&nz,2) ./ sum(nz,2);          % Jaccard
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% seconds of the day.
function sod = mydatesod (epoch, epoch0)
    if (nargin < 2),  epoch0 = [];  end
    if isempty(epoch)
        sod = epoch;
        return;
    end
    num0 = mydatesod_aux (epoch, epoch0);

    if (size(epoch,2) == 1)
        % epoch is in mydatenum format
        num = epoch;
        %no need to calculate vec.
        %vec = mydatevec(num);
    else
        % epoch is in mydatevec format
        vec = epoch;
        num = mydatenum(vec);
    end

    %%
    sod = (num - num0);
    %num, num0  % DEBUG
    %vec, vec0  % DEBUG
end

%!test
%! d = [2000 1 1 0 5 0];
%! sod_correct = 300;
%! sod_answer = mydatesod(mydatenum(d));
%! myassert (sod_answer, sod_correct, -sqrt(eps));

%!test
%! d  = [2000 1 2 0 5 0];
%! d0 = [2000 1 1 0 0 0];
%! sod_correct = 300 + 24*60^2;
%! sod_answer = mydatesod(mydatenum(d), d0);
%! myassert (sod_answer, sod_correct, -sqrt(eps));
%! sod_answer = mydatesod(mydatenum(d), mydatenum(d0));
%! myassert (sod_answer, sod_correct, -sqrt(eps));

%!test
%! in = NaN(2,1);
%! out = mydatesod(in);
%! myassert(out, in)


text:
% SP_PLOT_SOLUTION: Plot the computed solution, given the degrees of freedom.
%
%   [eu, F] = sp_plot_solution (u, space, geometry, pts, [ncuts=2]);
%   [eu, F] = sp_plot_solution (u, space, geometry, [npts], [ncuts=2]);
%
% INPUT:
%
%     u:           vector of dof weights
%     space:       object defining the discrete space (see sp_multipatch)
%     geometry:    an array of geometry structures (see mp_geo_load)
%     pts:         cell array with coordinates of points along each parametric direction
%     npts:        number of points along each parametric direction
%     ncuts:       only for volumetric domains, number of internal "cuts" in each parametric direction.
%                    The zero value will plot the solution on the boundary.
%
%    This function only plots the value of the solution. To plot other
%     quantities, such as the gradient, compute them with sp_eval.
%
% Copyright (C) 2016 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function sp_plot_solution (u, space, geometry, npts, ncuts)

if (nargin < 4)
  npts = [];
end
if (nargin < 5)
  ncuts = [];
end

if (isa (space.sp_patch{1}, 'sp_vector'))
  disp ('Warning: a different scaling is used for each patch')
end

hold_flag = ishold ;
for iptc = 1:space.npatch
  if (isempty (space.dofs_ornt))
    sp_plot_solution (u(space.gnum{iptc}), space.sp_patch{iptc}, geometry(iptc), npts, ncuts);
  else
    sp_plot_solution (u(space.gnum{iptc}) .* space.dofs_ornt{iptc}.', space.sp_patch{iptc}, geometry(iptc), npts, ncuts);
  end
  hold on
end

if (~hold_flag)
  hold off
end

end


text:
function I = displayDictionaryElementsAsImage(coefplots,D, numRows, numCols,X,Y,sortVarFlag)
% function I = displayDictionaryElementsAsImage(D, numRows, numCols, X,Y)
% displays the dictionary atoms as blocks. For activation, the dictionary D
% should be given, as also the number of rows (numRows) and columns
% (numCols) for the atoms to be displayed. X and Y are the dimensions of
% each atom.

borderSize = 1;
columnScanFlag = 1;
strechEachVecFlag = 1;
showImFlag = 1;

if (length(who('X'))==0)
    X = 8;
    Y = 8;
end
if (length(who('sortVarFlag'))==0)
    sortVarFlag = 0;
end

numElems = size(D,2);
if (length(who('numRows'))==0)
    numRows = floor(sqrt(numElems));
    numCols = numRows;
end
if (length(who('strechEachVecFlag'))==0)
    strechEachVecFlag = 0;
end
if (length(who('showImFlag'))==0)
    showImFlag = 1;
end

%%% sort the elements, if necessary.
%%% construct the image to display (I)
sizeForEachImage = sqrt(size(D,1))+borderSize;
I = zeros(sizeForEachImage*numRows+borderSize,sizeForEachImage*numCols+borderSize,3);
%%% fill all this image in blue
I(:,:,1) = 0;%min(min(D));
I(:,:,2) = 0; %min(min(D));
I(:,:,3) = 1; %max(max(D));

%%% now fill the image squares with the elements (in row scan or column
%%% scan).
if (strechEachVecFlag)
    for counter = 1:size(D,2)
        D(:,counter) = D(:,counter)-min(D(:,counter));
        if (max(D(:,counter)))
            D(:,counter) = D(:,counter)./max(D(:,counter));
        end
    end
end

if (sortVarFlag)
    vars = var(D);
    [V,indices] = sort(vars');
    indices = fliplr(indices);
    D = [D(:,1:sortVarFlag-1),D(:,indices+sortVarFlag-1)];
    signs = sign(D(1,:));
    signs(find(signs==0)) = 1;
    D = D.*repmat(signs,size(D,1),1);
    D = D(:,1:numRows*numCols);
end

zigzag256 = zigzag(16)+1;
    % work on columns
    zigzag256 = zigzag256(:);

%     [coefplots] = coefplot(coefs);
    [num ind] = sort(coefplots,1,'descend');
    ind = ind';
    newD = zeros(64,256);
    for i0 = 1:256
        newD(:,i0) = D(:,ind(i0));
    end

% for j = 1:numRows
%     for i = 1:numCols
%         if (strechEachVecFlag)
%             D(:,counter) = D(:,counter)-min(D(:,counter));
%             D(:,counter) = D(:,counter)./max(D(:,counter));
%         end
%         if (columnScanFlag==1)
for counter = 1:256
    zignum = zigzag256(counter);
    [i,j] = convertind(zignum,16,16);

            I(borderSize+(i-1)*sizeForEachImage+1:i*sizeForEachImage,borderSize+(j-1)*sizeForEachImage+1:j*sizeForEachImage,1)=reshape(newD(:,counter),8,8);
            I(borderSize+(i-1)*sizeForEachImage+1:i*sizeForEachImage,borderSize+(j-1)*sizeForEachImage+1:j*sizeForEachImage,2)=reshape(newD(:,counter),8,8);
            I(borderSize+(i-1)*sizeForEachImage+1:i*sizeForEachImage,borderSize+(j-1)*sizeForEachImage+1:j*sizeForEachImage,3)=reshape(newD(:,counter),8,8);
%         else
            % Go in Column Scan:
%             I(borderSize+(j-1)*sizeForEachImage+1:j*sizeForEachImage,borderSize+(i-1)*sizeForEachImage+1:i*sizeForEachImage,1)=reshape(D(:,counter),X,Y);
%             I(borderSize+(j-1)*sizeForEachImage+1:j*sizeForEachImage,borderSize+(i-1)*sizeForEachImage+1:i*sizeForEachImage,2)=reshape(D(:,counter),X,Y);
%             I(borderSize+(j-1)*sizeForEachImage+1:j*sizeForEachImage,borderSize+(i-1)*sizeForEachImage+1:i*sizeForEachImage,3)=reshape(D(:,counter),X,Y);
%         end
end
%     end
% end

if (showImFlag)
    I = I-min(min(min(I)));
    I = I./max(max(max(I)));
    imshow(I,[]);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [ x, w ] = legendre_ek_compute ( n )

%*****************************************************************************80
%
%% LEGENDRE_EK_COMPUTE: Legendre quadrature rule by the Elhay-Kautsky method.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2011
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer N, the order.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%

%
%  Define the zero-th moment.
%
  zemu = 2.0;
%
%  Define the Jacobi matrix.
%
  bj = zeros ( n, 1 );

  for i = 1 : n
    bj(i) = ( i * i ) / ( 4 * i * i - 1 );
  end
  bj(1:n) = sqrt ( bj(1:n) );

  x = zeros ( n, 1 );

  w = zeros ( n, 1 );
  w(1) = sqrt ( zemu );
%
%  Diagonalize the Jacobi matrix.
%
  [ x, w ] = imtqlx ( n, x, bj, w );

  w(1:n) = w(1:n).^2;

  return
end


text:
function y = prod(varargin)
% Fixes special case prod([],1), prod([],1,...), prod('',1), prod('',1,...)
% Fixes special case prod(sym([])), prod(sym([]),...) (without specifying dimension)
if nargin>1 && isequal(size(varargin{1}), [0 0]) && isequal(varargin{2}, 1)
    y = ones(1,0);
elseif ( nargin==1 && isa(varargin{1}, 'sym') && isequal(size(varargin{1}), [0 0]) ) || ...
       ( nargin>1  && isa(varargin{1}, 'sym') && isequal(size(varargin{1}), [0 0]) && ~isnumeric(varargin{2}) )
    y = 1;
else
    if ~isa(varargin{1},'sym')
        y = builtin('prod', varargin{:});
    else
        y = builtin('@sym/prod', varargin{:});
    end
end
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL