text:
% REJSTATEPOCH - reject bad eeg trials based a statistical measure. Can
%                 be applied either to the raw eeg data or the ICA
%                 component activity. This is an interactive function.
%
% Usage:
%   >> [ Irej, Irejdetails, n, threshold, thresholdg] = ...
%                         rejstatepoch( signal, rej, 'key1', value1...);
%
% Inputs:
%   signal     - 3 dimensional data array channel x points x trials
%                (instead of channels, one might also use independent
%                components).
%   rej        - rejection array, one value per trial and per channel or
%                component (size channel x trials).
%                By default, values are normalized by this
%                function and trehshold is expressed in term of standard
%                deviation of the mean.
%
% Optional inputs:
%   'plot'       - ['on'|'off'] interactive mode or just rejection.
%                  In the interactive mode, it plots the normalized
%                  entropy of original signal (blue) and the limits
%                  (red) (default:'on')
%   'threshold'  - percentage error threshold (default 1-0.25/nb_trials)
%                  for individual trials of individual channel/component.
%                  This threshold is expressed in term of standard
%                  deviation from the mean (default 5).
%   'global'     - ['on'|'off'], also perform threshold on the global
%                  measure (by default, the mean over all channel or
%                  electrodes).
%   'rejglob'    - rejection array, one value per trials. Use this
%                  argument when the global measure for all channels or
%                  does not correspond to their mean.
%   'thresholdg' - global threshold for the reunion of all channels
%                  or components. By default, it is equal to 'threshold'.
%   'normalize'  - ['on'|'off'], normalize values before applying the
%                  threshold. Default is 'on'.
%   'plotcom'    - string command to plot single trials. Default
%                  is none.
%   'title'      - title of the graph. Default is none.
%   'labels'     - labels for electrodes (array not cell).
%
% Outputs:
%   Irej        - indexes of trials to be rejected
%   Irejdetails - array for rejected components or channel (nb_rejected x
%                 nb_channel or nb_rejected x nb_components)
%   n           - number of trials rejected
%   threshold   - percentage error threshold
%   thresholdg  - percentage error threshold for global rejection
%
% See also: EEGLAB

% Algorithm:
%   normalise the measure given as input and reject trials based on
%   an uniform distribution of the data

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% userdata
% gcf : plotsig - signal to plot in the pop_out window
%       pnts - number of points per epoch
%       Irej - rejection of trials
%       Irejelec - rejection of trials and electrodes
%       labels - labels for curves
% plotwin : rej1 - rejection array (with all electrodes x trials)
%           rej2 - global rejection array (size trials)
%           thr1 - threshold (electrodes)
%           thr2 - threshold global

function [ Irej, Irejdetails, n, threshold, thresholdg] = rejstatepoch( signal, ...
		rej, varargin); % pnts, th_E, th_rejg, command, commandplot, typetitle, E, rejg);

if nargin < 1
	help rejstatepoch;
	return;
end

if ~ischar( signal )

	if nargin < 2
		help rejstatepoch;
		return;
	end

	if ~isempty( varargin ), g=struct(varargin{:});
	else g= []; end
	try, g.plot; 			catch, g.plot='on'; end
	try, g.threshold; 		catch, g.threshold=5; end
	try, g.thresholdg; 		catch, g.thresholdg=5; end
	try, g.global; 			catch, g.global='on'; end
	try, g.rejglob; 		catch, g.rejglob=[]; end
	try, g.normalize; 		catch, g.normalize='on'; end
	try, g.plotcom; 		catch, g.plotcom=''; end
	try, g.title;	 		catch, g.title=''; end
	try, g.labels;	 		catch, g.labels=''; end

	g.rej = rej;
	clear rej
	switch lower(g.plot)
		case {'on', 'off'} ;
		otherwise disp('Error: Plot must be either ''on'' or ''off'''); return;
	end;
	switch lower(g.global)
		case {'on', 'off'} ;
		otherwise disp('Error: Global must be either ''on'' or ''off'''); return;
	end;
	switch lower(g.normalize)
		case {'on', 'off'} ;
		otherwise disp('Error: Normalize must be either ''on'' or ''off'''); return;
	end;
	if ~ischar(g.plotcom)
		disp('Error: Plotcom must be a string to evaluate'); return;
	end;
	if ~ischar(g.title)
		disp('Error: Title must be a string'); return;
	end;
	try, g.threshold*2;
		catch, disp('Error: Threhsold must be a number'); return;
	end;
	try, g.thresholdg*2;
		catch, disp('Error: Threhsoldg must be a number'); return;
	end;
	if length(g.threshold(:)) > 1
		disp('Error: Threhsold must be a single number'); return;
	end;
	if length(g.thresholdg(:)) > 1
		disp('Error: Threhsoldg must be a single number'); return;
	end;
	if ~isempty(g.rejglob)
		if length(g.rejglob) ~= size(g.rej,2)
			disp('Error: Rejglob must be have the same length as rej columns'); return;
		end
	else
		switch lower(g.global), case 'on', g.rejglob = sum(g.rej,1); end
	end;
	if size(signal,3) ~= size(g.rej,2)
		disp('Error: Signal must be have the same number of element in 3rd dimension as rej have columns'); return;
	end
	if isempty(g.labels)
		for index = 1:size(g.rej,1)
			g.labels(index,:) = sprintf('%3d', index);
		end
		if ~isempty(g.rejglob)
			g.labels(index+2,:) = 'g. ';
		end
	end;

	switch lower(g.normalize),
		case 'on',
			g.rej = (g.rej-mean(g.rej,2)*ones(1, size(g.rej,2)))./ (std(g.rej, 0, 2)*ones(1, size(g.rej,2)));
			switch lower(g.global), case 'on', g.rejglob = (g.rejglob(:)-mean(g.rejglob(:)))./ std(g.rejglob(:)); end
	end;
	switch lower(g.global), case 'off',g.rejglob = []; end

	% plot the buttons
	% ----------------
	try, icadefs; catch, GUIBUTTONCOLOR = [0.8 0.8 0.8]; BACKCOLOR = [0.8 0.8 0.8]; end;
	figure('color', BACKCOLOR);
	set(gcf, 'name', 'Rejectrials');
	pos = get(gca,'position'); % plot relative to current axes
	set( gca, 'tag', 'mainaxis');
	q = [pos(1) pos(2) 0 0];
	s = [pos(3) pos(4) pos(3) pos(4)]./100; % allow to use normalized position [0 100] for x and y
	axis('off');

	plotsig = sum(abs(signal(:,:)),1);
	set(gcf, 'userdata', { plotsig, size(signal, 2), [], [] }); % the two last arguments are the rejection

	% Create axis
	% -----------
	h6 = axes('Units','Normalized', 'tag',  'Plotwin', 'Position',[-10 12 120 84].*s+q);
	title(g.title);
	set( h6, 'userdata', { g.rej g.rejglob g.threshold g.thresholdg g.labels }); % g.rej was put twice because it is used to compute the global entropy

	% CANCEL button
	% -------------
	h  = uicontrol(gcf, 'Style', 'pushbutton', 'string', 'Cancel', 'Units','Normalized','Position',...
			[-15 -6 15 6].*s+q, 'callback', 'close(gcf);', 'backgroundcolor', GUIBUTTONCOLOR);

	% Entropy component text and threshold
	% ------------------------------------
    makebutton( 'Single-channel', 'Plotwin'  ,  [3 -6 27 6].*s+q, [30 -6 15 6].*s+q, 3, g.threshold, GUIBUTTONCOLOR);
    makebutton( 'All-channel', 'Plotwin'     ,  [50 -6  27 6].*s+q, [77 -6  15 6].*s+q, 4, g.thresholdg, GUIBUTTONCOLOR);

	% EEGPLOT button
	% --------------
	if ~isempty(g.plotcom)
			h  = uicontrol(gcf, 'backgroundcolor', GUIBUTTONCOLOR, 'Style', 'pushbutton', 'string', 'EEGPLOT', 'Units','Normalized','Position',[95 -2 15 6].*s+q, ...
					'callback',['TMPEEG = get(gcbf, ''userdata'');' ...
								'TMPREJ = TMPEEG{3};' ...
								'TMPREJELEC = TMPEEG{4};' ...
								g.plotcom ] );
			posaccept = [95 -10 15 6];
	else	posaccept = [95 -6 15 6];
	end;

	% ACCEPT button
	% -------------
	command = 'fprintf(''Rejection indices has been put in the matrix TMPREJ\n'')';
	haccept  = uicontrol(gcf, 'backgroundcolor', GUIBUTTONCOLOR, 'Style', 'pushbutton', 'string', 'UPDATE', 'Units','Normalized','Position', posaccept.*s+q, ...
					'callback', [	'set(gcbo, ''userdata'', 1);' ... %signal to signify termination
									'TMPEEG = get(gcbf, ''userdata'');' ...
									'TMPREJ = TMPEEG{3};' ...
									'TMPREJELEC = TMPEEG{4};' ...
									command ] );

	command = [ 'entwin = get(gcbf, ''currentobject'');' ...
				'tmptype = get(entwin,''type'');' ...
				'if tmptype(1) == ''l'' entwin = get(entwin, ''parent''); end;' ...
				'tagwin = get(entwin,''tag'');' ... % either entropy or kurtosis
				'switch tagwin,' ...
				' case ''Plotwin'',' ... % check if the user clicked on the right window
				'   alldata = get(gcf, ''userdata'');' ...
				'   plotsig = alldata{1};' ...
				'   pnts = alldata{2};' ...
				'   fig = figure(''position'', [100 300 600 400],''color'', [1 1 1]);' ...
				'   I = alldata{3};' ...
				'   sweeps = size(plotsig,2) /pnts;' ...
				'   h1 = axes(''parent'', fig, ''Units'',''Normalized'', ''Position'',[0.6 0.11 0.27 0.815]);' ...
				'   pos = get(entwin, ''currentpoint'');' ...
				'   component = round(pos(1) / sweeps + 0.5);' ... % determine the component number
				'   alldata = get(entwin, ''userdata'');' ...
		    	'   rej = alldata{1};' ...
				'   if component <= size(rej,1)' ... % component
				'        rej_threshold = alldata{3};' ...
				'        component = max(1,min(component, size(rej,1))); ' ...
		    	'   	 rej = rej(component, :);' ...
				'        titlegraph = sprintf(''' g.title ' #%d'', component);' ...
				'        colorgraph = ''b'';' ...
				'   else' ...                        % global
		    	'   	 rej = alldata{2};' ...
				'        rej_threshold = alldata{4};' ...
				'        titlegraph = sprintf(''' g.title ' global'');' ...
				'        colorgraph = ''g'';' ...
				'   end;' ...
				'   plot([1:length(rej)], rej, ''color'', colorgraph);' ...
				'   ss = get(h1, ''xlim'');' ...
				'   set(h1, ''view'', [90 90]);' ...
				'   set(h1, ''xdir'', ''reverse'');' ...
				'   set(h1, ''XLim'', ss);' ...
				'   hold on;' ...  % plot component
				'   yl = get(h1, ''ylim'');' ...
				'   set(h1, ''xtickmode'', ''manual'', ''xtick'', sweeps/2, ''xticklabel'', component, ''xlim'', [ 0 sweeps ]);' ...
				'   title( titlegraph );' ...
				'   plot( get(h1, ''xlim''), [rej_threshold rej_threshold], ''r'');' ... % plot limit
				'   plot( get(h1, ''xlim''), [-rej_threshold -rej_threshold], ''r'');' ... % plot limit
				'   set(h1, ''xticklabel'', []);' ...
				'   hold on;' ...
				'   h2 = axes(''parent'', fig,''Units'',''Normalized'', ''Position'',[0.13 0.11 0.27 0.815]);' ...
				'   erpimage( plotsig, ones(1,size(plotsig,2)/pnts), [0:1/(pnts-1):1], '''', 3, 1, ''nosort'', ''noplot'');' ...
				'   title(''Currentset all chans''); xlabel(''''); ylabel(''''); ' ...
				'   set(gca, ''xticklabel'', []);' ...
				'   hold on;' ...
				'   h3 = axes(''parent'', fig,''Units'',''Normalized'', ''Position'',[0.45 0.11 0.1 0.815]);' ...
				'   if any(I ~= 0)' ...
				'      rejImage = (I'' * ones(1, 10))'';' ...
				'      imagesc( rejImage'' );' ...
				'      set(gca, ''ydir'', ''normal'');' ...
				'   end;' ...
				'   title(''Rejected (all)''); xlabel(''''); ylabel('''');' ...
				'   set(gca, ''xticklabel'', [], ''yticklabel'', []);' ...
				'end;' ...
				'clear fig tmptype tagwin alldata rej rejImage plotsig sweeps pnts rej_threshold ss q s h1 h2 h3 pos component yl;' ];

%				'      erpimage( rejImage(:)'', ones(1,size(I,2)), [0:1/(10-1):1], '''', 1, 0, ''nosort'', ''noplot'');' ...
	set(gcf, 'WindowButtonDownFcn', command);

	rejstatepoch('draw');
	switch g.plot,
		case 'on', waitfor( haccept, 'userdata'); drawnow;
	end

	threshold  = g.threshold;
	thresholdg = g.thresholdg;
	Irej = [];
	Irejdetails = [];
	n = 0;
	try
		TMPEEG = get(gcf, 'userdata');
		Irej = TMPEEG{3};
		Irejdetails = TMPEEG{4};
		n = length(find(Irej == 1));

		plothandler = findobj( 'parent', gcf, 'tag', 'Plotwin');
		TMPEEG = get(plothandler, 'userdata');
		threshold = TMPEEG{3};
		thresholdg = TMPEEG{4};
		close(gcf);
	catch, end
else %if signal is a string draw everything

	% retrieve data
	% -------------
	gcfdata = get(gcf, 'userdata');
	plotsig = gcfdata {1};
	pnts    = gcfdata {2};
    sweeps  = size(plotsig,2)/pnts;

	h6 = findobj('parent', gcf, 'tag', 'Plotwin');
	alldata = get(h6, 'userdata');
	g.rej       = alldata {1};
	g.rejg      = alldata {2};
	g.threshold   = alldata {3};
	g.thresholdg  = alldata {4};
	set(h6, 'userdata', alldata);

	nbchans = size(g.rej,1);

	% reject trials
	% -------------
	rejelec = abs(g.rej) > g.threshold;
	rej  = max(rejelec,[],1);
	n1 = sum(rej(:));
	if ~isempty( g.rejg )
		rej2 = abs(g.rejg) > g.thresholdg;
		n2 = sum(rej2(:));
		rej = rej | rej2(:)';
	end
	fprintf('%d trials rejected (single:%d, all:%d)\n', sum(rej(:)), n1, n2);
	gcfdata {3} = rej;
	gcfdata {4} = rejelec;
	set(gcf, 'userdata', gcfdata);

	% plot the sorted entropy curve
	% -----------------------------
	plotstat( 'Plotwin');

end
return;

function plotstat( id );

	% plot the sorted entropy curve
	% -----------------------------
	h6 = findobj('parent', gcf, 'tag', id);
	axes(h6); cla;
	ttmp = get(gca, 'title');
	oldtitle = get(ttmp, 'string');

	% get datas
    % ---------
   	alldata = get(gca, 'userdata');
	g.rej       = alldata {1};
	g.rejg      = alldata {2};
	g.threshold   = alldata {3};
	g.thresholdg  = alldata {4};
	g.labels      = alldata {5};
	nbchans = size(g.rej,1);
	sweeps  = size(g.rej,2);

	% plot datas
    % ----------
	g.rej = g.rej'; plot(g.rej(:)); g.rej = g.rej';
	hold on;
	yl = get(gca, 'ylim');

	% plot vertival bars to separate components and the trehsold
	% ----------------------------------------------------------
	set( gca, 'tag',  id, 'ylimmode', 'manual');
	set(gca, 'xtickmode', 'manual', 'xtick', [0:sweeps:(size(g.rej(:),1)-1+2*sweeps)] + sweeps/2, ...
			 'xticklabel', g.labels, 'xlim', [ 0 (size(g.rej(:),1)-1+2*sweeps)]);
	plot( [1 size(g.rej(:),1)], [-g.threshold -g.threshold], 'r');	% plot threshold
	plot( [1 size(g.rej(:),1)], [g.threshold g.threshold], 'r');	% plot threshold

	if ~isempty(g.rejg) % plot global ?
		plot([size(g.rej(:),1)+sweeps:size(g.rej(:),1)+2*sweeps-1],  g.rejg(:), 'g');
		pp = patch([size(g.rej(:),1) size(g.rej(:),1) size(g.rej(:),1)+sweeps size(g.rej(:),1)+sweeps], [yl(1)-1 yl(2)+1 yl(2)+1 yl(1)-1], get(gcf, 'color'), 'clipping', 'off');
		set(pp, 'EdgeColor',  get(gcf, 'color'));
		plot( [size(g.rej(:),1)+sweeps length(g.rejg)+size(g.rej(:),1)+sweeps], [-g.thresholdg -g.thresholdg], 'r');	% plot threshold
		plot( [size(g.rej(:),1)+sweeps length(g.rejg)+size(g.rej(:),1)+sweeps], [g.thresholdg g.thresholdg], 'r');	% plot threshold
		plot([size(g.rej(:),1)+sweeps size(g.rej(:),1)+sweeps], yl, 'k');
	else
		pp = patch([size(g.rej(:),1) size(g.rej(:),1) size(g.rej(:),1)+2*sweeps size(g.rej(:),1)+2*sweeps], [yl(1)-1 yl(2)+1 yl(2)+1 yl(1)-1], get(gcf, 'color'), 'clipping', 'off');
		set(pp, 'EdgeColor',  get(gcf, 'color'));
	end
	for index = 0:sweeps:size(g.rej(:),1);
		plot([index index], yl, 'k');
	end

	% restore properties
	title(oldtitle);
	set(h6, 'userdata', alldata);

return;

function makebutton( string, ax, pos1, pos2, userindex, init, GUIBUTTONCOLOR );
	h  = uicontrol(gcf , 'backgroundcolor', GUIBUTTONCOLOR, 'Style', 'radiobutton', 'string', string, 'value', fastif(init == 0, 0, 1), 'Units','Normalized', 'Position', pos1, ...
					'callback', [ 'textresh = findobj(''parent'', gcbf, ''tag'', ''' string ''');' ...
								  'checkstatus = get(gcbo, ''value'');' ...
								  'ax = findobj(''parent'', gcbf, ''tag'', ''' ax ''');' ...
								  'tmpdat = get(ax, ''userdata'');' ...
								  'if checkstatus' ... % change status of the textbox
								  '   set(textresh, ''enable'', ''on'');' ...
								  '   tmpdat{' int2str(userindex) '} = str2num(get(textresh, ''string''));' ...
								  'else' ...
								  '   set(textresh, ''enable'', ''off'');' ...
								  '   tmpdat{' int2str(userindex) '} = 0;' ...
								  'end;' ...
								  'set(ax, ''userdata'' , tmpdat);' ...
								  'rejstatepoch(''draw'');' ...
								  'clear tmpdat ax checkstatus textresh;'  ] );
	h  = uicontrol(gcf, 'Style', 'edit', 'backgroundcolor', [1 1 1], 'tag', string, 'string', num2str(init), 'enable', fastif(init == 0, 'off', 'on'), 'Units','Normalized', 'Position', pos2, ...
					'callback', [ 	'ax = findobj(''parent'', gcbf, ''tag'', ''' ax ''');' ...
								    'tmpdat = get(ax, ''userdata'');' ...
									'tmpdat{' int2str(userindex) '} = str2num(get(gcbo, ''string''));' ...
									'set(ax, ''userdata'' , tmpdat);' ...
								    'rejstatepoch(''draw'');' ...
									'clear tmpdat ax;' ] );
return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [outputSignal, inputImages] = computeSignalsFromImages(inputImages,inputMovie, varargin)
	% CIAtah package wrapper for applyImagesToMovie.
	% Applies images to a 3D movie matrix in order to get a signal based on a thresholded version of the image.
	% Biafra Ahanonu
	% started: 2013.10.11 [2020.10.27 [12:56:27] for wrapper]
	% inputs
		% inputImages - [x y signalNo] of images, signals will be calculated for each image from the movie.
		% inputMovie - [x y frame] or char string path to the movie.
	% outputs
		% outputSignal - [signalNo frame] matrix of each signal's activity trace extracted directly from the movie.
		% inputImages - [x y signalNo], same as input.

	% changelog
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

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

	% ========================
	% DESCRIPTION
	% OPTIONS ARE THE SAME AS loadMovieList.
	% options.exampleOption = '';
	% get options
	% options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	% ========================

	try
		[outputSignal, inputImages] = applyImagesToMovie(inputImages,inputMovie, 'passArgs',varargin);
	catch err
		disp(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		disp(repmat('@',1,7))
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [AC] = ACbuilder(X, Spec)
%   Creates an aircraft from decision variables contained in X and the base
%   AC, which contains project-fixed parameters and function handles.

AC = Spec.baseAC;

%% Add decision variables
AC.Vars = cell2struct(num2cell(X),Spec.MyDesignXlabels);

%% Brute force loop method
%initial weight guess
AUMestimate = AC.Vars.Wpayload*1000; %can't be less than this!

changefactor = 0; %factor used for convergence
tol = 1e-7; %fraction equivalent to .1kg tolerance for a 10000kg helicopter
err = 1;
ii = 0;
while abs(err)>tol;
    ii = ii+1;
    AC.W.AUM = AUMestimate;
    AC = add_sized_parts(AC);
    [AUMestimatenew,AC] = AC.Ass.weightestimate(AC);
    err = AUMestimate - AUMestimatenew;
    AUMestimate = AUMestimatenew + err*changefactor;
    err = err/AUMestimatenew;
end

AC.W.fuel = AC.W.AUM*AC.Vars.FuelFrac;
AC.W.W = AC.W.AUM*g0;
end

function AC = add_sized_parts(AC)

AC.W.fuelfrac = AC.Vars.FuelFrac;
AC.W.W = AC.W.AUM*g0;

%% Build Rotor
AC.Rotor.R =    sqrt(AC.W.AUM/(pi*AC.Vars.DiscLoading));
AC.Rotor.TS =   AC.Vars.TipSpeed;
AC.Rotor.b =    round(AC.Vars.Blades);
AC.Rotor.c =    AC.Vars.Solidity*pi*AC.Rotor.R/AC.Rotor.b;
AC.Rotor.tw =   AC.Vars.Twist*pi/180;

AC.Rotor.D = AC.Rotor.R*2;
AC.Rotor.Omega = AC.Rotor.TS./AC.Rotor.R;
AC.Rotor.DA = pi*AC.Rotor.R.^2;
AC.Rotor.AR = AC.Rotor.R/AC.Rotor.c;
AC.Rotor.BA = AC.Rotor.b*AC.Rotor.c*AC.Rotor.R;
AC.Rotor.sigma = AC.Rotor.BA/AC.Rotor.DA;

%% Build Wing
AC.Wing.S =     AC.Vars.SpecificWingArea*AC.W.AUM/1000;
AC.Wing.AR =    AC.Vars.WingAR;

AC.Wing.b = sqrt(AC.Wing.S*AC.Wing.AR);

%% Build Power
AC.Power.P =    AC.Vars.P_to_W*AC.W.AUM*1000; %X in kW/kg
AC.Power.T =    AC.Vars.T_to_W*AC.W.W; %X in kN

%% Build Tail Rotor
AC = simpleTRsizing(AC);
% gives us R, sigma, c, Omega, TS, and b for tail rotor

end

author: Sable

repo: mcbench-benchmarks


text:
function hermite_product_display ( name, i, j )

%*****************************************************************************80
%
%% HERMITE_PRODUCT_DISPLAY displays a 2D Hermite product such as H(i,x)*H(j,y).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string NAME, the name of the polynomial:
%    'h', physicist's Hermite;
%    'hn', normalized physicist's Hermite;
%    'he', probabilist's Hermite;
%    'hen', normalized probabilist's Hermite;
%    'hf', physicist's Hermite * exp(-x^2).
%
%    Input, integer I, J, the indices of the two polynomials.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_PRODUCT_DISPLAY:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Display the function defined by a 2d product of two \n' );
  fprintf ( 1, '  Hermite polynomials like H(i,x) * H(j,y).\n' );
%
%  Get missing input.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The polynomial name must be one of:\n' );
    fprintf ( 1, '  ''h'', physicist''s Hermite\n' );
    fprintf ( 1, '  ''hn'', normalized physicist''s Hermite\n' );
    fprintf ( 1, '  ''he'', probabilist''s Hermite\n' );
    fprintf ( 1, '  ''hen'', normalized probabilist''s Hermite\n' );
    fprintf ( 1, '  ''hf'', physicist''s Hermite * exp(-x^2)\n' );
    name = input ( '  Enter the polynomial name:  ' );
  end

  if ( nargin < 2 )
    i = input ( '  Enter the polynomial index I:  ' );
  end

  if ( nargin < 3 )
    j = input ( '  Enter the polynomial index J:  ' );
  end
%
%  Set up the graph.
%
  x_min = - 2.0;
  x_max = + 2.0;
  y_min = - 2.0;
  y_max = + 2.0;
%
%  Clear the graphics page.
%
  clf
  nx = 101;
  x = r8vec_linspace ( nx, x_min, x_max );
  ny = 101;
  y = r8vec_linspace ( ny, y_min, y_max );
  [ X, Y ] = meshgrid ( x, y );

  nx = 101;
  if ( s_eqi ( name, 'h' ) )
    hx = h_polynomial ( nx, i, x );
  elseif ( s_eqi ( name, 'hn' ) )
    hx = hn_polynomial ( nx, i, x );
  elseif ( s_eqi ( name, 'he' ) )
    hx = he_polynomial ( nx, i, x );
  elseif ( s_eqi ( name, 'hen' ) )
    hx = hen_polynomial ( nx, i, x );
  elseif ( s_eqi ( name, 'hf' ) )
    hx = hf_function ( nx, i, x );
  end
  hx = hx(1:nx,i+1);
  ny = 101;
  if ( s_eqi ( name, 'h' ) )
    hy = h_polynomial ( ny, j, y );
  elseif ( s_eqi ( name, 'hn' ) )
    hy = hn_polynomial ( ny, j, y );
  elseif ( s_eqi ( name, 'he' ) )
    hy = he_polynomial ( ny, j, y );
  elseif ( s_eqi ( name, 'hen' ) )
    hy = hen_polynomial ( ny, j, y );
  elseif ( s_eqi ( name, 'hf' ) )
    hy = hf_function ( ny, j, y );
  end
  hy = hy(1:ny,j+1);

  HXY = hy * hx';
  z_min = min ( min ( HXY ) );
  z_max = max ( max ( HXY ) );

  if ( z_min == z_max )
    z_min = z_min - 0.5;
    z_max = z_max + 0.5;
  end

  caxis ( [ -0.4, 1.2 ] )
  surface ( X, Y, HXY, 'FaceColor', 'interp' )

  shading interp

  axis ( [ x_min, x_max, y_min, y_max, z_min, z_max ] );
  axis equal

  xlabel ( '--X axis--' );
  ylabel ( '--Y axis--' );
  zlabel ( '--Z axis--' );

  title_string = sprintf ( '%s(%d,x)*%s(%d,y)', name, i, name, j );
  title ( title_string );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_PRODUCT_DISPLAY:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function c = ch_cap ( c )

%*****************************************************************************80
%
%% CH_CAP capitalizes a single character.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, the character to capitalize.
%
%    Output, character C, the capitalized character.
%
  if ( 'a' <= c && c <= 'z' )
    c = c + 'A' - 'a';
  end

  c = char ( c );

  return
end
function p = h_polynomial ( m, n, x )

%*****************************************************************************80
%
%% H_POLYNOMIAL evaluates the physicist's Hermite polynomials H(n,x).
%
%  Differential equation:
%
%    Y'' - 2 X Y' + 2 N Y = 0
%
%  First terms:
%
%      1
%      2 X
%      4 X^2     -  2
%      8 X^3     - 12 X
%     16 X^4     - 48 X^2     + 12
%     32 X^5    - 160 X^3    + 120 X
%     64 X^6    - 480 X^4    + 720 X^2    - 120
%    128 X^7   - 1344 X^5   + 3360 X^3   - 1680 X
%    256 X^8   - 3584 X^6  + 13440 X^4  - 13440 X^2   + 1680
%    512 X^9   - 9216 X^7  + 48384 X^5  - 80640 X^3  + 30240 X
%   1024 X^10 - 23040 X^8 + 161280 X^6 - 403200 X^4 + 302400 X^2 - 30240
%
%  Recursion:
%
%    H(0,X) = 1,
%    H(1,X) = 2*X,
%    H(N,X) = 2*X * H(N-1,X) - 2*(N-1) * H(N-2,X)
%
%  Orthogonality:
%
%    Integral ( -oo < X < +oo ) exp ( - X^2 ) * H(M,X) H(N,X) dX
%    = sqrt ( pi ) * 2^N * N! * delta ( N, M )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Frank Olver, Daniel Lozier, Ronald Boisvert, Charles Clark,
%    NIST Handbook of Mathematical Functions,
%    Cambridge University Press, 2010,
%    ISBN: 978-0521192255,
%    LC: QA331.N57.
%
%  Parameters:
%
%    Input, integer M, the number of evaluation points.
%
%    Input, integer N, the highest order polynomial to compute.
%    Note that polynomials 0 through N will be computed.
%
%    Input, real X(M,1), the evaluation points.
%
%    Output, real P(M,N+1), the values of the physicist's Hermite polynomials
%    of index 0 through N.
%
  p = zeros ( m, n + 1 );

  p(1:m,1) = 1.0;

  if ( n == 0 )
    return
  end

  p(1:m,2) = 2.0 * x(1:m,1);

  for j = 2 : n
    p(1:m,j+1) = 2.0 * x(1:m,1) .* p(1:m,j) - 2.0 * ( j - 1 ) * p(1:m,j-1);
  end

  return
end
function p = he_polynomial ( m, n, x )

%*****************************************************************************80
%
%% HE_POLYNOMIAL evaluates the probabilist's Hermite polynomials at X.
%
%  Differential equation:
%
%    ( exp ( - 0.5 * x^2 ) * He(n,x)' )' + n * exp ( - 0.5 * x^2 ) * He(n,x) = 0
%
%  First terms:
%
%   1
%   X
%   X^2  -  1
%   X^3  -  3 X
%   X^4  -  6 X^2 +   3
%   X^5  - 10 X^3 +  15 X
%   X^6  - 15 X^4 +  45 X^2 -   15
%   X^7  - 21 X^5 + 105 X^3 -  105 X
%   X^8  - 28 X^6 + 210 X^4 -  420 X^2 +  105
%   X^9  - 36 X^7 + 378 X^5 - 1260 X^3 +  945 X
%   X^10 - 45 X^8 + 630 X^6 - 3150 X^4 + 4725 X^2 - 945
%
%  Recursion:
%
%    He(0,X) = 1,
%    He(1,X) = X,
%    He(N,X) = X * He(N-1,X) - (N-1) * He(N-2,X)
%
%  Orthogonality:
%
%    Integral ( -oo < X < +oo ) exp ( - 0.5 * X^2 ) * He(M,X) He(N,X) dX
%    = sqrt ( 2 * pi ) * N! * delta ( N, M )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Frank Olver, Daniel Lozier, Ronald Boisvert, Charles Clark,
%    NIST Handbook of Mathematical Functions,
%    Cambridge University Press, 2010,
%    ISBN: 978-0521192255,
%    LC: QA331.N57.
%
%  Parameters:
%
%    Input, integer M, the number of evaluation points.
%
%    Input, integer N, the highest order polynomial to compute.
%    Note that polynomials 0 through N will be computed.
%
%    Input, real X(M,1), the evaluation points.
%
%    Output, real P(M,N+1), the values of the probabilist's Hermite polynomials
%    of index 0 through N.
%
  p = zeros ( m, n + 1 );

  p(1:m,1) = 1.0;

  if ( n == 0 )
    return
  end

  p(1:m,2) = x(1:m,1);

  for j = 2 : n
    p(1:m,j+1) = x(1:m,1) .* p(1:m,j) - ( j - 1 ) * p(1:m,j-1);
  end

  return
end
function p = hen_polynomial ( m, n, x )

%*****************************************************************************80
%
%% HEN_POLYNOMIAL evaluates the normalized probabilist's Hermite polynomials.
%
%  Discussion:
%
%    These polynomials satisfy the orthonormality condition:
%
%      Integral ( -oo < X < +oo ) exp ( - 0.5 * X^2 ) * Hen(M,X) Hen(N,X) dX
%      = delta ( N, M )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Frank Olver, Daniel Lozier, Ronald Boisvert, Charles Clark,
%    NIST Handbook of Mathematical Functions,
%    Cambridge University Press, 2010,
%    ISBN: 978-0521192255,
%    LC: QA331.N57.
%
%  Parameters:
%
%    Input, integer M, the number of evaluation points.
%
%    Input, integer N, the highest order polynomial to compute.
%    Note that polynomials 0 through N will be computed.
%
%    Input, real X(M,1), the evaluation points.
%
%    Output, real P(M,N+1), the values of the polynomials of index 0 through N.
%
  p = zeros ( m, n + 1 );

  p(1:m,1) = 1.0;

  if ( n == 0 )
    return
  end

  p(1:m,2) = x(1:m,1);

  for j = 2 : n
    p(1:m,j+1) = x(1:m,1) .* p(1:m,j) - ( j - 1 ) * p(1:m,j-1);
  end
%
%  Normalize.
%
  d = diag ( 1.0 ./ sqrt ( gamma ( 1:n+1) .* sqrt ( 2.0 * pi ) ) );

  p = p * d;

  return
end
function f = hf_function ( m, n, x )

%*****************************************************************************80
%
%% HF_FUNCTION evaluates the Hermite function.
%
%  Discussion:
%
%    The Hermite function of degree n is related to the physicist's
%    Hermite polynomial H(n,x):
%
%      Hf(n,x) = H(n,x) * exp ( - 0.5 * x^2 ) / sqrt ( 2^n n! sqrt ( pi ) )
%
%    The Hermite functions are orthonormal:
%
%      Integral ( -oo < x < +oo ) Hf(m,x) Hf(n,x) dx = delta ( m, n )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Frank Olver, Daniel Lozier, Ronald Boisvert, Charles Clark,
%    NIST Handbook of Mathematical Functions,
%    Cambridge University Press, 2010,
%    ISBN: 978-0521192255,
%    LC: QA331.N57.
%
%  Parameters:
%
%    Input, integer M, the number of evaluation points.
%
%    Input, integer N, the highest order polynomial to compute.
%    Note that polynomials 0 through N will be computed.
%
%    Input, real X(M,1), the evaluation points.
%
%    Output, real F(M,N+1), the values of the Hermite functions
%    of index 0 through N at the evaluation points.
%
  f = zeros ( m, n + 1 );

  f(1:m,1) =       exp ( - 0.5 * x(1:m,1).^2 )             / sqrt (       sqrt ( pi ) );

  if ( n == 0 )
    return
  end

  f(1:m,2) = 2.0 * exp ( - 0.5 * x(1:m,1).^2 ) .* x(1:m,1) / sqrt ( 2.0 * sqrt ( pi ) );

  for j = 2 : n
    f(1:m,j+1) = ( sqrt ( 2.0 ) * x(1:m,1) .* f(1:m,j) - sqrt ( j - 1 ) * f(1:m,j-1) ) ...
      / sqrt ( j );
  end

  return
end
function p = hn_polynomial ( m, n, x )

%*****************************************************************************80
%
%% HN_POLYNOMIAL evaluates the normalized physicist's Hermite polynomials.
%
%  Discussion:
%
%    These polynomials satisfy the orthonormality condition:
%
%      Integral ( -oo < X < +oo ) exp ( - X^2 ) * Hn(M,X) Hn(N,X) dX = delta ( N, M )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Frank Olver, Daniel Lozier, Ronald Boisvert, Charles Clark,
%    NIST Handbook of Mathematical Functions,
%    Cambridge University Press, 2010,
%    ISBN: 978-0521192255,
%    LC: QA331.N57.
%
%  Parameters:
%
%    Input, integer M, the number of evaluation points.
%
%    Input, integer N, the highest order polynomial to compute.
%    Note that polynomials 0 through N will be computed.
%
%    Input, real X(M,1), the evaluation points.
%
%    Output, real P(M,N+1), the values of the polynomials of index 0 through N.
%
  p = zeros ( m, n + 1 );

  p(1:m,1) = 1.0;

  if ( n == 0 )
    return
  end

  p(1:m,2) = 2.0 * x(1:m,1);

  for j = 2 : n
    p(1:m,j+1) = 2.0 * x(1:m,1) .* p(1:m,j) - 2.0 * ( j - 1 ) * p(1:m,j-1);
  end
%
%  Normalize.
%
  d = diag ( 1.0 ./ sqrt ( gamma ( 1:n+1) .* 2.^(0:n) * sqrt ( pi ) ) );

  p = p * d;

  return
end
function a = r8vec_linspace ( n, a_first, a_last )

%*****************************************************************************80
%
%% R8VEC_LINSPACE creates a vector of linearly spaced values.
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%    While MATLAB has the built in command
%
%      a = linspace ( a_first, a_last, n )
%
%    that command has the distinct disadvantage of returning a ROW vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, real A_FIRST, A_LAST, the first and last entries.
%
%    Output, real A(N), a vector of linearly spaced data.
%
  a = zeros ( n, 1 );

  if ( n == 1 )
    a(1,1) = ( a_first + a_last ) / 2.0;
  else
    a = ( linspace ( a_first, a_last, n ) )';
  end

  return
end
function value = s_eqi ( s1, s2 )

%*****************************************************************************80
%
%% S_EQI is a case insensitive comparison of two strings for equality.
%
%  Example:
%
%    S_EQI ( 'Anjana', 'ANJANA' ) is 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S1, S2, the strings to compare.
%
%    Output, logical VALUE, is TRUE if the strings are equal.
%
  len1 = length ( s1 );
  len2 = length ( s2 );
  lenc = min ( len1, len2 );

  for i = 1 : lenc

    c1 = ch_cap ( s1(i) );
    c2 = ch_cap ( s2(i) );

    if ( c1 ~= c2 )
      value = 0;
      return
    end

  end

  for i = lenc + 1 : len1
    if ( s1(i) ~= ' ' )
      value = 0;
      return
    end
  end

  for i = lenc + 1 : len2
    if ( s2(i) ~= ' ' )
      value = 0;
      return
    end
  end

  value = 1;

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef SPX_SteepestDescent < handle
% Steepest descent algorithm implementation

    properties
        % Settable and gettable properties

        % Maximum number of iterations for which the algorithm can run
        MaxIterations
        % Threshold of norm (in terms of percentage)
        NormThreshold
    end

    properties(SetAccess=private)
        % Gettable properties

        % The sparse symmetric positive definite matrix operator
        A
        % The vectors to be solved
        B
        % Problem dimension
        N
        % Number of equations to be solved
        S
        % Number of iterations taken for solving the problem
        Iterations
        % Residual norms at the end
        ResidualNorms
        % The solution vectors
        X
        % Residual vectors
        Residuals
        % Indicates if the problems converged
        % This would be false only if the algorithm crossed max iterations
        Converged
    end

    methods
        % Public methods
        function self = SPX_SteepestDescent(A, B)
            % Constructor
            if isa(A, 'spx.dict.Operator')
                self.A = A;
            elseif ismatrix(A)
                self.A = spx.dict.MatrixOperator(A);
            else
                error('Unsupported operator.');
            end
            self.B = B;
            [self.N, self.S] = size(B);
            self.MaxIterations = self.N * 20;
            self.NormThreshold = 1e-6;
        end
        function result = solve(self)
            aa = self.A;
            bb = self.B;
            % Initial estimate vectors are all zeros
            xx = zeros(self.N, self.S);
            result = xx;
            self.Iterations = zeros(self.S, 1);
            self.ResidualNorms = zeros(self.S, 1);
            self.Converged = false(self.S, 1);
            % Initial residual vectors
            rr = bb  - aa * xx;
            self.Residuals  = rr;
            % Norm squared of initial residual vectors
            deltas = spx.norm.inner_product_cw(rr, rr);
            % The factor with which the norm needs to be reduced
            epsilon = self.NormThreshold;
            % Target limits on norm squared of residuals
            limits = epsilon^2 * deltas;
            % Maximum number of iterations for which the algorithm
            % is allowed to run
            imax = self.MaxIterations;
            % Number of problems being solved
            ns = self.S;
            for s=1:ns
                % Initialize iteration counter
                i = 0;
                % The quantities for this problem
                limit = limits(s);
                delta = deltas(s);
                r = rr(:, s);
                x = xx(:, s);
                b = bb(:, s);
                while i < imax && delta > limit
                    % Compute the intermediate variable
                    q = aa * r;
                    % the line search scale factor in the direction of steepest descent
                    alpha = delta / (r' * q);
                    % Update estimate in the direction of steepest descent
                    x = x + alpha * r;
                    if mod(i , 50) == 0
                        % In order to avoid propagation of floating point
                        % errors, we will recompute the value of residual
                        r = b - aa * x;
                    else
                        % Otherwise we use a shortcut
                        r = r  - alpha * q;
                    end
                    % Update residual norm squared
                    delta = r' * r;
                    % Increase iteration counter
                    i = i + 1;
                end
                % The problem has been solved
                result(:, s) = x;
                % Number of iterations taken to solve this problem
                self.Iterations(s) = i;
                self.ResidualNorms(s) = sqrt(delta);
                self.Residuals(:, s) = r;
                self.Converged(s) = (delta <= limit);
            end
            % Maintain the result for reference
            self.X = result;
        end

        function result = hasConverged(self)
            % Returns if all the solutions have converged
            result = all(self.Converged);
        end

        function printResults(self)
            ns = self.S;
            nn = self.N;
            for s = 1:ns
                fprintf('Problem: %d\n', s);
                fprintf('Iterations: %d\n', self.Iterations(s));
                fprintf('Residual norm: %.2f, Converged: %d\n', ...
                    self.ResidualNorms(s), self.Converged(s));
                if nn < 10
                    % We will print the solutions too
                    fprintf('Solution vector: ');
                    fprintf('%.4f ', self.X(:, s));
                    fprintf('\n');
                    fprintf('Residual vector: ');
                    fprintf('%.4f ', self.Residuals(:, s));
                    fprintf('\n');
                end
            end
        end
    end

    methods(Access=private)
        % Private methods

    end

    methods(Static)
        % Public static methods

    end

end

author: indigits

repo: sparse-plex


text:
% Originally Coded By: Niall Mangan
% Copyright 2017, All Rights Reserved
% Code by Niall Mangan for paper "Inferring biological networks by sparse
% identification of nonlinear dynamics"
% by N. M. Mangan S. L. Brunton, J. L. Proctor, and J. N. Kutz
% Original Code: https://github.com/niallmm/iSINDy
%%
% This file is modified to compare the performance of i-SINDy on
% insufficient data.
% Modified by:K, 2019/07/16
% Find equations for yeast glycolysis state variable 2
%%
clc;close all;clear all;

addpath('./utils');
addpath('./bioutils');

% define libarary parameters
laurentorder = 0;
polyorder = 6;
usesine = 0;
dyorder = 1;

% search for equation for the 2nd state-variable
% the pareto front and details for the 2nd variable were plotted in the paper
% clear results for other state variables
clear Theta Thetastring Xi indTheta lambdavec numterms errorv indopt
clear indTheta1 Xi1 numterms1 nT

% Here we load the previously simulated data
load('TrainingData.mat')

% Define the new data length
percent=0.01;new_length=round(percent*length(xt));

% Shuffel the original data
Sequence=randperm(size(xt,1));
xt=xt(Sequence,:);
dxt=dxt(Sequence,:);

% Assign the value to the new variables
Data=xt(1:new_length,:);dData=dxt(1:new_length,:);

% Define the number of states
n=size(Data,2);

% pool Data  (i.e., build library of nonlinear time series)
[Theta, Thetastring] = poolDatady(Data,n,polyorder,usesine, laurentorder, dData(:,2), dyorder);

% compute Sparse regression using ADM
tol = 2e-3;

lambda = 2e-3;

jj = 1; % counter
num= 1; % initialize the number of nonzero terms found for the lambda
errorvec= 0;

MaxIter = 1e3;

% for now calculate null space using null function
nT = null(Theta);

% Define the parameters for plooting
plottag=2;

[indTheta1, Xi1, numterms1] = ADMinitvary(nT,lambda,MaxIter,tol, plottag);

Thetastring(indTheta1)'
n0Xi = Xi1(Xi1~=0); % terms need to be rearranged to recover coefficients
n0Xi/n0Xi(end)

save('Results/2nd_state_variable_1.mat')

author: dynamicslab

repo: SINDy-PI


text:
function [As,A,eil,Ei] = indexed_sparse(i,j,v,m,n,varargin)
% INDEXED_SPARSE Create a sparse matrix with indexed edges.
%
% [As,A,eil,Ei] = indexed_sparse(i,j,v,m,n) creates a sparse matrix A just
% like A = sparse(i,j,v,m,n).  However, indexed_sparse returns additional
% information.  The matrix As is a structural matrix for A which
% corresponds to As = sparse(i,j,1,m,n).  Thus, As(i,j) != 0 for all edges.
% The vector eil is a permutation for the vector v, such that v(eil) is the
% correct input for the edge_weight parameter.  The matrix Ei lists the
% index of each edge in the vector v, so that
% A = sparse(j,i,v(nonzeros(Ei)),m,n)' unless options.istrans = 0.
%
% This function handles the case when v(k) == 0.  For v(k) = 0,
% A(i(k),j(k)) = 0, but As(i(k),j(k)) = 1, and the vector v(eil) provides
% an appropriate input to the edge_weight parameter for all the algorithms.
%
% See the examples reweighted_edges for more information.
%
% ... = indexed_sparse(A,...) takes a set of
% key-value pairs or an options structure.  See set_matlab_bgl_options
% for the standard options.
%    options.undirected: output edge indices for an undirected graph [{0} | 1]
%      See the note about undirected inputs.
%
% Note (Undirected inputs): If options.undirected = 1, the input to the
% graph still must contain both undirected edges and the corresponding
% weight.
%
% Example:
%   % see example/reweighted_edges

% David Gleich
% Copyright, Stanford University, 2007-2008

%% History
%  2007-07-13: Changed input options to use undirected as the option name.
%  2008-10-07: Changed options parsing
%%

[trans check]  = get_matlab_bgl_options(varargin{:});

options = struct('undirected', 0);
options = merge_options(options, varargin{:});

symmetric = options.undirected;

Ei = accumarray([i j], 1:length(v), [m n], @min, 0, true);
if symmetric, Ei = min(Ei,Ei'); end

A = sparse(i, j, v, m, n);
As = spones(Ei);

if trans, eil = nonzeros(Ei');
else eil = nonzeros(Ei);
end

if check
    % make sure there were no elements in i,j repeated
    % that would cause problems
    if length(v) ~= nnz(Ei)
        warning('matlab_bgl:indexed_sparse', ...
            'duplicated elements detected, using the first occurrence as the index');
    end

    if symmetric
        [ei ej] = find(As');
        if ~isequal(A,sparse(ei,ej,v(eil),m,n))
            warning('matlab_bgl:indexed_sparse', ...
                'the output matrix failed a symmetry test.  This indicates a non-symmetric value array v.');
        end
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function layerstruct = riverplot_layer_names(layerstruct, layernames, varargin)
%
% This is a subfunction of the riverplot toolbox.
%  - Adds or removes names from a layer struct, and corresponding text on plot
%  - Toggles between adding and removing, depending on whether names exist
%  in layerstruct.  Layerstruct is created with riverplot_draw_layer.
%
% :Usage:
% ::
%
%     layerstruct = riverplot_layer_names(layerstruct, layernames, [optional inputs])
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2016  Tor Wager
%
%     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/>.
% ..
%
% :Inputs:
%
%   **layerstruct:**
%        Structure with handles and decription of imrect class objects for
%        a layer. Layerstruct is created with riverplot_draw_layer.
%
%   **layernames:**
%        A cell array of strings for names
%
% :Optional Inputs:
%
%
% :Outputs:
%
%   **layerstruct:**
%        structure with names and name info added
%
%
% :Examples:
% ::
%
%
% :References:
%
%
% :See also:
%   - list other functions related to this one, and alternatives*
%

% ..
%    Programmers' notes:
%    Created by Tor Wager, July 2, 2016
%
%    8/21/2017 Stephan Geuter
%       - updated text alignment
% ..

% ..
%    DEFAULTS AND INPUTS
% ..

ax = gca;
colors = custom_colors([.2 .2 .2], [.2 .2 .2], length(layerstruct)); % series of cells with colors for each layer
leftright_string = 'left';

if ~iscell(layernames), layernames = cellstr(layernames); end

% optional inputs with default values
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case 'ax', ax = varargin{i+1}; varargin{i+1} = [];
            %case 'position', myposition = varargin{i+1}; varargin{i+1} = [];
            case {'color', 'colors'}, colors = varargin{i+1}; varargin{i+1} = [];

            case 'left', leftright_string = 'left';

            case 'right', leftright_string = 'right';

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

axes(ax)

n_rects = length(layerstruct);

if n_rects > length(layernames)
    error('Not enough names entered to riverplot_layer_names.m');
end

for i = 1:n_rects

    % Get y position: midpoint of layer
    yval = layerstruct{i}.topleft(2) - .5 * (layerstruct{i}.topleft(2) - layerstruct{i}.bottomleft(2));

    % Get x position:
    switch leftright_string
        case 'left'

            xval = layerstruct{i}.topleft(1) - .3;

            hh = text(xval, yval, layernames{i}, 'Color', colors{i}, 'FontSize', 16,...
                      'HorizontalAlignment','right');

        case 'right'

            xval = layerstruct{i}.topright(1) + .3;

            hh = text(xval, yval, layernames{i}, 'Color', colors{i}, 'FontSize', 16,...
                      'HorizontalAlignment','left');
    end

% removed because text property 'HorizontalAlignment' is better and adjusts
% after re-sizing the figure. SG 8/21/2017
%     hh = text(xval, yval, layernames{i}, 'Color', colors{i}, 'FontSize', 16);
%
%     switch leftright_string
%         case 'left'
%             hh = shift_text(hh);
%     end

    layerstruct{i}.name_handles(i) = hh;
    layerstruct{i}.names(i) = layernames(i);

end

end % function

% function text_handle = shift_text(text_handle)
% % Center text - Adjust for width of text string
% % Given text handle, subtract 1/2 the width of text string from x and y positions
% % Borrowed/adapted from John Barber, http://www.mathworks.com/matlabcentral/fileexchange/30671-calcticks
%
% textExt = get(text_handle, 'Extent');
% textHeight = textExt(4);
% textWidth = textExt(3);
%
% % If using a proportional font, shrink text width by a fudge factor to
% % account for kerning.
% % ax = gca;
% %
% % if ~strcmpi(ax.FontName,'FixedWidth')
% %     textWidth = textWidth*0.8;
% % end
%
% mypos = get(text_handle, 'Position');
% mypos(1) = mypos(1) - textWidth;
% set(text_handle, 'Position', mypos);
%
% end


text:
function x = fejer2 ( n )

%*****************************************************************************80
%
%% FEJER2 returns the Type 2 Fejer points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2018
%
%  Author:
%
%    John Burkardt.
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real X(N), the points.
%
  theta(1:n) = ( n : -1 : 1 ) * pi / ( n + 1 );
  x(1:n) = cos ( theta(1:n) );

  return
end


text:
function [evaluation_info time]=evaluate(B_train, B_test, gnd_train, gnd_test)
start_t=clock;
Ntrain = length(gnd_train);
Ntest = length(gnd_test);
% Ntest = 100;

ap_all = zeros(Ntest,1);
ph2_all = zeros(Ntest,1);
for i = 1:Ntest
    % compute your distance
    D_code = hammingDist(B_test(i,:) ,B_train);

    % evaluation
    [ap ph2] = phmap(D_code, gnd_train, gnd_test(i));
    ap_all(i)=ap;
    ph2_all(i)=ph2;
end

evaluation_info.AP=mean(ap_all);
evaluation_info.PH2=mean(ph2_all);
end_t=clock;
time=etime(end_t,start_t);

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
function wh_holdout_final = xval_select_holdout_set(Y, covs, nfolds, holdout_proportion, varargin)
    %
    % wh_holdout_final = xval_select_holdout_set(Y, covs, nfolds, holdout_proportion, [doplot])
    %
    % Selects one or more holdout sets balancing on an outcome and set of
    % covariates.  The idea is inspired by Rubin's propensity score method, and
    % the logic is that selecting holdout sets that are balanced on the outcome
    % and other variables of no interest will reduce the variance in the
    % estimated predictive accuracy in a cross-validated parameter optimization
    % scheme.  The goal of this algorithm is thus to choose holdout sets where
    % the probability of holdout assignment is balanced on (independent of) the covariates.
    %
    % For assessing accuracy with little bias, a leave-one-out holdout method
    % can be a good choice, as it has minimum bias.  However, the accuracy
    % estimates have high variance, which often motivates k-fold schemes with
    % larger holdout sets.  For model/parameter selection (e.g., in inner
    % cross-validation schemes, or when model selection is the primary goal),
    % larger holdout sets and small training sets may be desirable.  See Shao,
    % 1998, JASA.
    %
    % It is also advantageous to have training/holdout sets that are independent
    % of one another, but this is not optimized here.
    %
    % Input variables:
    % Y = outcome variable
    % covs = covariates to balance on
    % nfolds = number of holdout sets desired
    % holdout_proportion = how much of the sample is reserved for holdout, [0 1]
    %
    % Output variables:
    % wh_holdout = a matrix of [obs x nfolds], containing optimized holdout sets
    %
    % Example inputs:
    % Y = randn(30, 1);
    % covs = randn(N, 1);
    % nfolds = 10;
    % holdout_proportion = .5;
    %
    % Example function call:
    % wh_holdout = xval_select_holdout_set(Y, covs, 20, .5);
    %
    % Notes:
    % There is a plot option hard-coded in the script for demonstration
    % purposes.

    doplot = 0;
    if ~isempty(varargin), doplot = varargin{1}; end

    anyzero = 1;
    wh_holdout_final = [];

    wh_holdout_final = generate_balanced_sets(Y, covs, nfolds, holdout_proportion, doplot);

    % find any missing observations
    % Make sure all observations are selected equally often, to the
    % degree possible
    % Find the holdout sets that select the most frequently used observations
    obssum = sum(wh_holdout_final, 2);
    ismissing = double(obssum == 0);

    while any(ismissing)

        wh_holdout_new = generate_balanced_sets(Y, covs, nfolds, holdout_proportion, doplot); % obs x sets

        new_to_save = wh_holdout_new' * ismissing; % these holdout sets have the missing obs from before

        %frequent_obs_score = wh_holdout_new' * obssum; % higher values mean the holdout set uses more frequently selected observations (bad)

        % Replace the ones with the most frequent scores with the new ones
        % to keep
        frequent_obs_score = wh_holdout_final' * obssum; % higher values mean the holdout set uses more frequently selected observations (bad)
        [f2, indx] = sort(frequent_obs_score, 1, 'descend');

        if any(new_to_save)

            new_to_save = find(new_to_save);
            wh_to_replace = indx(1:length(new_to_save));
            wh_holdout_final(:, wh_to_replace) = wh_holdout_new(:, new_to_save);

            % it is possible that we've replaced one that causes another obs to
            % be missing, so we need the while loop
            obssum = sum(wh_holdout_final, 2);
            ismissing = double(obssum == 0);
        else

            %% do nothing   ismissing = 1;

        end

    end

end % function

function wh_holdout_final = generate_balanced_sets(Y, covs, nfolds, holdout_proportion, doplot)

    N = length(Y);

    niter = nfolds * 10;
    [betavals, tvals] = deal(zeros(1, niter));
    wh_holdout = false(N, niter);

    warning off % iteration limit...

    for i = 1:niter

        wh_holdout(:, i) = zeros(N, 1);
        nh = round(N .* holdout_proportion);

        wh = randperm(N);
        wh = wh(1:nh);
        wh_holdout(wh, i) = 1;

        % For testing purposes
        %covs = mvnrnd([1 3], [1 .4; .4 1], N); % random outcome and covs

        % p-value is linear with respect to betas near the origin (50%)
        % so we can select subsets based on beta alone.
        %
        b = glmfit([Y covs], wh_holdout(:, i), 'binomial', 'link', 'logit');

        % This is used for testing purposes
        if doplot
            [b, dev, stat] = glmfit([Y covs], wh_holdout(:, i), 'binomial', 'link', 'logit');

            %
            %         % plot partial fit wrt Y
            %         yfit = [min(Y):.1:max(Y)];
            %         xfit = glmval(b(1:2), yfit, 'logit');
            %
            %         create_figure; plot(Y, wh_holdout(:, i), 'ko','MarkerFaceColor', 'k');
            %         set(gca, 'XLim', [-.5 1.5]);
            %         plot(yfit, xfit, 'r');
            %         ylabel('P(holdout)'); xlabel('outcome');

            tvals(i) = stat.t(2);
        end

        betavals(i) = b(2);
        %

    end

    warning on

    % Save the most balanced holdout sets

    cutoff = prctile(abs(betavals), 10);  % save top 10% -- that ensures us the right size holdout set
    wh_to_keep = abs(betavals) <= cutoff;
    wh_holdout_final = wh_holdout(:, wh_to_keep);

    if doplot
        create_figure('plot'); plot_correlation_samefig(betavals, tvals);
        plot(betavals(wh_to_keep), tvals(wh_to_keep), 'go');
        xlabel('Logistic reg beta'); ylabel('T-value');
        title('Selection set estimates for logistic reg, outcome')

        % plot partial fit wrt Y for one example holdout set
        [b, dev, stat] = glmfit([Y covs], wh_holdout(:, 1), 'binomial', 'link', 'logit');

        yfit = [min(Y):.1:max(Y)];
        xfit = glmval(b(1:2), yfit, 'logit');

        create_figure('Example holdout set'); plot(Y, wh_holdout(:, 1), 'ko','MarkerFaceColor', 'k');
        set(gca, 'XLim', [-.5 1.5]);
        plot(yfit, xfit, 'r');
        ylabel('P(holdout)'); xlabel('Outcome');
        axis auto
        title('Example holdout set: P(holdout|Y)')
        drawnow
    end

end


text:
function [ UH ] = uh_5_half( d_base, delta_t )
%uh_5_half Unit Hydrograph [days] with half a triangle (exponential decay)

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

%   Inputs
%   d_base  - time base of routing delay [d]
%   delta_t - time step size [d]
%
%   Output
%   UH      - unit hydrograph [nx2]
%               uh's first row contains coeficients to splut flow at each
%               of n timesteps forward, the second row contains zeros now,
%               these are the still-to-flow values.
%
%   Unit hydrograph spreads the input volume over a time period delay.
%   Percentage of input returned only decreases.
%   I.e. d_base = 3.8 [days], delta_t = 1:
%   UH(1) = 0.841  [% of inflow]
%   UH(2) = 0.133
%   UH(3) = 0.021
%   UH(4) = 0.004

%%TIME STEP SIZE
delay = d_base/delta_t;
if delay == 0; delay = 1; end       % any value below t = 1 means no delay,
                                    % but zero leads to problems
tt = 1:ceil(delay);                 % time series for which we need UH
                                    % ordinates [days]

%%UNIT HYDROGRAPH
% The Unit Hydrograph follows exponential decay y=exp(-x). For a given
% delay time, the fraction of flow per time step is thus the integral of
% t-1 to t of the exponential decay curve. The curve has range [0,Inf>.
% We impose the arbitrary boundary of [0,7] here (exp(-7) = 9e-4) as the
% point where the decay curve 'ends'. This allows to divide the range [0,7]
% in n delay steps, and so calculate the UH.

%%Find integral limits
stepsize        = (7-0)/delay;      % Range over which the decay curve is
                                    % calculated, divided by required
                                    % number of delay steps
limits          = 0:stepsize:7;
limits(end+1)   = 7;

%%EMPTIES
UH = zeros(1,length(tt));

%%UNIT HYDROGRAPH
for t = 1:length(tt)
    UH(t) = integral(@(x) exp(-x),limits(t),limits(t+1));
end

%%ACCOUNT FOR <7,Inf> PART OF THE CURVE (i.e. add the missing tail end of
% the curve to the last delay step, to ensure that 100% of flow is routed).
UH(end) = UH(end)+(1-sum(UH));

UH(2,:) = zeros(size(UH));

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
classdef MotionCorrection < handle

    % MotionCorrection object for the NoRMCorre algorithm
    properties
        file ='';            % filename
        dims;                % dimensionality of the FOV
        nd;                  % number of dimensions (2d or 3d)
        T;                   % number of timesteps
        options;             % motion correction options structure
        shifts;              % inferred shifts
        shifts_x;            % shifts in x direction (matrix format)
        shifts_y;            % shifts in y direction (matrix format)
        shifts_z;            % shifts in z direction (matrix format)
        template;            % inferred template
        M;                   % motion corrected file
        meanY;               % mean image raw data
        meanM;               % mean image corrected data
        corrY                % correlation with mean raw data
        corrM;               % correlation with mean corrected data
        crispY;              % crispness of mean image (raw)
        crispM;              % crispness of mean image (corrected)
        col_shift = 0;       % relative shift due to bidirectional scanning
        gSig = 7;            % std of gaussian kernel for high pass filtering
        gSiz = 21;           % size of high pass filter
        file_orig;           % original file (used for 1p correction)
    end

    methods
        %% construct object, lazily load file and set dimensionality variables
        function obj = MotionCorrection(filename)
            obj.file = filename;
            if ischar(filename)
                Ytemp = read_file(filename,1,1);
                obj.dims = size(Ytemp);
            else    % filename is an array already loaded in memory
                sizY = size(filename);
                obj.T = sizY(end);          % assuming more than 1 frames
                obj.dims = sizY(1:end-1);
            end
            obj.nd = length(obj.dims);
            obj.options.d1 = obj.dims(1);
            obj.options.d2 = obj.dims(2);
            if obj.nd == 3; obj.options.d3 = obj.dims(3); end
        end

        %% load and update options
        function optionsSet(obj, varargin)
            obj.options = NoRMCorreSetParms(obj.options,varargin{:});
        end

        %% motion correct serial
        function motionCorrectSerial(obj,options)
            if exist('options','var')
                obj.optionsSet(obj.options,options);
            end
            obj.options.d1 = obj.dims(1);
            obj.options.d2 = obj.dims(2);
            if obj.nd == 3; obj.options.d3 = obj.dims(3); end
            [obj.M,obj.shifts,obj.template,obj.options,obj.col_shift] = normcorre(obj.file,obj.options);
            obj.T = length(obj.shifts);
            shifts_res = cat(ndims(obj.shifts(1).shifts)+1,obj.shifts(:).shifts);
            shifts_res = reshape(shifts_res,[],obj.nd,obj.T);
            obj.shifts_x = squeeze(shifts_res(:,1,:))';
            obj.shifts_y = squeeze(shifts_res(:,2,:))';
            if obj.nd == 3; obj.shifts_z = squeeze(shifts_res(:,3,:))'; end
        end

        %% motion correct Parallel(obj)
        function motionCorrectParallel(obj,options)
            if exist('options','var')
                obj.optionsSet(obj.options,options);
            end
            obj.options.d1 = obj.dims(1);
            obj.options.d2 = obj.dims(2);
            if obj.nd == 3; obj.options.d3 = obj.dims(3); end
            [obj.M,obj.shifts,obj.template,obj.options,obj.col_shift] = normcorre_batch(obj.file,obj.options);
            obj.T = length(obj.shifts);
            shifts_res = cat(ndims(obj.shifts(1).shifts)+1,obj.shifts(:).shifts);
            shifts_res = reshape(shifts_res,[],obj.nd,obj.T);
            obj.shifts_x = squeeze(shifts_res(:,1,:))';
            obj.shifts_y = squeeze(shifts_res(:,2,:))';
            if obj.nd == 3; obj.shifts_z = squeeze(shifts_res(:,3,:))'; end
        end

        %% apply shifts to a different file
        function M = applyShifts(obj,newfile)
            M = apply_shifts(newfile,obj.shifts,obj.options,0,0,0,obj.col_shift);
        end

        %% high pass filter
        function HPF(obj,gSig,gSiz)
            if exist('gSig','var'); obj.gSig = gSig; end
            if exist('gSiz','var'); obj.gSiz = gSiz; end
            psf = fspecial('gaussian', round(gSiz), gSig);
            ind_nonzero = (psf(:)>=max(psf(:,1)));
            psf = psf-mean(psf(ind_nonzero));
            psf(~ind_nonzero) = 0;
            if ischar(obj.file)
                obj.file_orig = read_file(obj.file);
            else
                obj.file_orig = obj.file;
            end
            obj.file = imfilter(single(obj.file_orig),psf,'symmetric');
        end

        %% compute mean
        function computeMean(obj)
            if obj.nd == 2; batch = 1000; elseif obj.nd == 3; batch = 10; end
            if isempty(obj.meanY)
                if ischar(obj.file)
                    obj.meanY = zeros(obj.dims,'single');
                    for t = 1:batch:obj.T
                        Y_temp = single(read_file(obj.file,t,batch));
                        obj.meanY = ((t-1)*obj.meanY+size(Y_temp,obj.nd+1)*mean(Y_temp,obj.nd+1))...
                                        /(t-1+size(Y_temp,obj.nd+1));
                    end
                else
                    obj.meanY = mean(obj.Y,obj.nd+1);
                end
            end
            if isempty(obj.meanM)
                if ischar(obj.M)
                    obj.meanM = zeros(obj.dims,'single');
                    for t = 1:batch:obj.T
                        M_temp = single(read_file(obj.M,t,batch));
                        obj.meanM = ((t-1)*obj.meanM+size(M_temp,obj.nd+1)*mean(M_temp,obj.nd+1))...
                                        /(t-1+size(M_temp,obj.nd+1));
                    end
                else
                    obj.meanM = mean(obj.M,obj.nd+1);
                end
            end
        end

        %% crispness (norm of gradient of mean)
        function crispness(obj,bnd)
            if ~exist('bnd','var'); bnd = zeros(6,1); end
            if isscalar(bnd); bnd = bnd*ones(6,1); end
            computeMean(obj);
            if obj.nd == 2
                mean_Y = obj.meanY(bnd(1)+1:end-bnd(2),bnd(3)+1:end-bnd(4));
                mean_M = obj.meanM(bnd(1)+1:end-bnd(2),bnd(3)+1:end-bnd(4));
            elseif obj.nd == 3
                mean_Y = obj.meanY(bnd(1)+1:end-bnd(2),bnd(3)+1:end-bnd(4),bnd(5)+1:end-bnd(6));
                mean_M = obj.meanM(bnd(1)+1:end-bnd(2),bnd(3)+1:end-bnd(4),bnd(5)+1:end-bnd(6));
            end
            if ismatrix(obj.meanY)
                [gx,gy] = gradient(mean_Y);
                ng = norm(sqrt(gx.^2+gy.^2),'fro');
            elseif obj.nd == 3
                [gx,gy,gz] = gradient(mean_Y);
                ng = sqrt(sum(gx(:).^2+gy(:).^2+gz(:).^2));
            end
            obj.crispY = ng;
            if ismatrix(obj.meanM)
                [gx,gy] = gradient(mean_M);
                ng = norm(sqrt(gx.^2+gy.^2),'fro');
            elseif obj.nd == 3
                [gx,gy,gz] = gradient(mean_M);
                ng = sqrt(sum(gx(:).^2+gy(:).^2+gz(:).^2));
            end
            obj.crispM = ng;
        end

        %% correlation with mean
        function correlationMean(obj,bnd)
            if ~exist('bnd','var'); bnd = zeros(6,1); end
            if isscalar(bnd); bnd = bnd*ones(6,1); end
            computeMean(obj);
            if obj.nd == 2
                mean_Y = obj.meanY(bnd(1)+1:end-bnd(2),bnd(3)+1:end-bnd(4));
                mean_M = obj.meanM(bnd(1)+1:end-bnd(2),bnd(3)+1:end-bnd(4));
            elseif obj.nd == 3
                mean_Y = obj.meanY(bnd(1)+1:end-bnd(2),bnd(3)+1:end-bnd(4),bnd(5)+1:end-bnd(6));
                mean_M = obj.meanM(bnd(1)+1:end-bnd(2),bnd(3)+1:end-bnd(4),bnd(5)+1:end-bnd(6));
            end
            if obj.nd == 2; batch = 1000; elseif obj.nd == 3; batch = 10; end
            if ~isempty(obj.file)
                if ischar(obj.file)
                    obj.corrY = zeros(obj.T,1,'single');
                    for t = 1:batch:obj.T
                        Y_temp = single(read_file(obj.file,t,batch));
                        if obj.nd == 2; Y_temp = Y_temp(bnd(1)+1:end-bnd(2),bnd(3)+1:end-bnd(4),:); end
                        if obj.nd == 3; Y_temp = Y_temp(bnd(1)+1:end-bnd(2),bnd(3)+1:end-bnd(4),bnd(5)+1:end-bnd(6),:); end
                        obj.corrY(t:t+size(Y_temp,obj.nd+1)-1) = corr(reshape(Y_temp,[],size(Y_temp,obj.nd+1)),mean_Y(:));
                    end
                else
                    if obj.nd == 2; obj.corrY = corr(mean_Y(:),single(reshape(obj.Y(bnd(1)+1:end-bnd(2),bnd(3)+1:end-bnd(4),:),[],obj.T))); end
                    if obj.nd == 3; obj.corrY = corr(mean_Y(:),single(reshape(obj.Y(bnd(1)+1:end-bnd(2),bnd(3)+1:end-bnd(4),bnd(5)+1:end-bnd(6),:),[],obj.T))); end
                end
            end
            if ~isempty(obj.M)
                if ischar(obj.M)
                    obj.meanM = zeros(obj.T,1,'single');
                    for t = 1:batch:obj.T
                        M_temp = single(read_file(obj.M,t,batch));
                        if obj.nd == 2; M_temp = M_temp(bnd(1)+1:end-bnd(2),bnd(3)+1:end-bnd(4),:); end
                        if obj.nd == 3; M_temp = M_temp(bnd(1)+1:end-bnd(2),bnd(3)+1:end-bnd(4),bnd(5)+1:end-bnd(6),:); end
                        obj.corrM(t:t+size(M_temp,obj.nd+1)-1) = corr(reshape(M_temp,[],size(M_temp,obj.nd+1)),mean_M(:));
                    end
                else
                    if obj.nd == 2; obj.corrM = corr(mean_M(:),single(reshape(obj.M(bnd(1)+1:end-bnd(2),bnd(3)+1:end-bnd(4),:),[],obj.T))); end
                    if obj.nd == 3; obj.corrM = corr(mean_M(:),single(reshape(obj.M(bnd(1)+1:end-bnd(2),bnd(3)+1:end-bnd(4),bnd(5)+1:end-bnd(6),:),[],obj.T))); end
                end
            end
        end

    end
end

author: flatironinstitute

repo: NoRMCorre


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni

% Relationship between DFT and DTFT.
% The Transforms are plotted in -pi<w<pi

n=0:11;
x=0.9.^n;
syms w
Xdtft=sum(x.*exp(-j*w*n));
Xdft=fft(x);
Xdft.'

N=12;
k=-N/2: N/2-1;
wk=2*pi*k/N;
wk'

Xshift= fftshift(Xdft)
Xshift.'
ezplot(abs(Xdtft),[-pi,pi]);
hold ;
plot(wk,abs(Xshift),'o')
hold on
ylim([0 8])
legend('DTFT','DFT')
title('|DTFT| & |DFT|')
hold off

figure
w1=-pi:.01:pi;
Xdtft=subs(Xdtft,w,w1);
plot(w1,angle(Xdtft));
hold on
plot(wk,angle(Xshift),'o')
legend('DTFT','DFT')
xlim([-pi pi])
title('\angle DTFT &\angle DFT')
hold off


text:
function p17_demo ( iteration_max, h )

%*****************************************************************************80
%
%% P17_DEMO runs the 2D demo problem #17, with mesh size H.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 July 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer ITERATION_MAX, the maximum number of iterations that DISTMESH
%    should take.  (The program might take fewer iterations if it detects convergence.)
%
%    Input, real H, the mesh spacing parameter.
%
  if ( nargin < 1 )
    iteration_max = 200;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P17_DEMO - Note:\n' );
    fprintf ( 1, '  No value of ITERATION_MAX was supplied.\n' );
    fprintf ( 1, '  The default value ITERATION_MAX = %d will be used.\n', ...
      iteration_max );
  end

  if ( nargin < 2 )
    h = 0.10;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P17_DEMO - Note:\n' );
    fprintf ( 1, '  No value of H was supplied.\n' );
    fprintf ( 1, '  The default value H = %f will be used.\n', h );
  end
%
%  Put the random number generator into a fixed initial state.
%
  rand ( 'state', 111 );
%
%  Set the rendering method for the current figure to Z-buffering.
%
  set ( gcf, 'rend', 'z' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Problem 17:\n' );
  fprintf ( 1, '  Rectangle with a Reuleuax triangle hole, h = %f\n', h )

  fd = @p17_fd;
  fh = @p17_fh;
  box = [ 0.0, -1.0; 5.0, 1.0 ];
  fixed = [];

  [ p, t ] = distmesh_2d ( fd, fh, h, box, iteration_max, fixed );

  post_2d ( p, t, fh )
%
%  Write a PostScript image of the triangulation.
%
  [ node_num, junk ] = size ( p );
  [ tri_num, junk ] = size ( t );
  p = p';
  t = t';
  node_show = 0;
  triangle_show = 1;

  triangulation_order3_plot ( 'p17_mesh.eps', node_num, p, tri_num, ...
    t, node_show, triangle_show );
%
%  Write a text file containing the nodes.
%
  r8mat_write ( 'p17_nodes.txt', 2, node_num, p );
%
%  Write a text file containing the triangles.
%
  i4mat_write ( 'p17_elements.txt', 3, tri_num, t );

  return
end


text:
function plottingMomentumParameter
n = 20000;
lambda(1) = 0;
gamma(1) = 0;
for i = 2:n
    lambda(i) = 0.5*(1 + sqrt(1+4*lambda(i-1)^2));
    gamma(i) = (lambda(i-1) - 1)/lambda(i);
end

plot(1:n,[lambda;gamma],'+-')
legend('\lambda','\gamma','Location','Best')
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function y = rosen(x)
%ROSEN	Calculate Rosenbrock's function.
%
%	Description
%	Y = ROSEN(X) computes the value of Rosenbrock's function at each row
%	of X, which should have two columns.
%
%	See also
%	DEMOPT1, ROSEGRAD
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Calculate value of Rosenbrock's function: x should be nrows by 2 columns

y = 100 * ((x(:,2) - x(:,1).^2).^2) + (1.0 - x(:,1)).^2;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% Demonstrates how to use the matrix editor.
%
% See also: MatrixEditor

% Copyright 2010 Levente Hunyadi
function example_matrixeditor

fig = figure( ...
    'MenuBar', 'none', ...
    'Name', 'Matrix editor demo - Copyright 2010 Levente Hunyadi', ...
    'NumberTitle', 'off', ...
    'Toolbar', 'none');
editor = MatrixEditor(fig, ...
    'Item', [1,2,3,4;5,6,7,8;9,10,11,12], ...
    'Type', PropertyType('denserealdouble','matrix'));
uiwait(fig);
disp(editor.Item);

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [spectrum,f,s] = MTPointSpectrum(times,varargin)

%MTPointSpectrum - Compute point process spectrum by multi-taper estimation.
%
%  The spectrum is computed using the chronux toolbox.
%
%  USAGE
%
%    [spectrum,f,s] = MTPointSpectrum(times,<options>)
%
%    times          event times, e.g. spike times
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'frequency'   sampling rate (in Hz, default = 20,000 Hz)
%     'range'       frequency range (in Hz) (default = all)
%     'tapers'      relative resolution and order of the tapers [NW K]
%                   (default = [3 5])
%     'pad'         FFT padding (see help for mtspectrumpt) (default = 0)
%     'show'        plot spectrum (default = 'off')
%    =========================================================================
%
%  OUTPUT
%
%    spectrum       spectrum
%    s              error
%    f              frequency bins
%
%  DEPENDENCIES
%
%    This function requires the chronux toolbox.
%
%  SEE
%
%    See also MTSpectrum, MTSpectrogram, SpectrogramBands.

% Copyright (C) 2010-2014 by Michaël Zugaro
%
% 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.

% Make sure chronux is installed and functional
CheckChronux('mtspectrumpt');

% Defaults
frequency = 20000;
range = [];
show = 'off';
tapers = [3 5];
pad = 0;
err = [1 0.95];

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

% Check parameter sizes
if size(times,2) ~= 1 && size(times,2) ~= 2,
	error('Parameter ''times'' is not a vector or a Nx2 matrix (type ''help MTPointSpectrum'' for details).');
end

% Parse parameter list
v = {};
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+2) ' is not a property (type ''help MTPointSpectrum'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'frequency',
			frequency = varargin{i+1};
			if ~isdscalar(frequency,'>0'),
				error('Incorrect value for property ''frequency'' (type ''help MTPointSpectrum'' for details).');
			end
		case 'range',
			range = varargin{i+1};
			if ~isdvector(range,'#2','<','>=0'),
				error('Incorrect value for property ''range'' (type ''help MTPointSpectrum'' for details).');
			end
		case 'tapers',
			tapers = varargin{i+1};
			if ~isivector(tapers,'#2'),
				error('Incorrect value for property ''tapers'' (type ''help MTPointSpectrum'' for details).');
			end
		case 'pad',
			pad = varargin{i+1};
			if ~isiscalar(pad,'>-1'),
				error('Incorrect value for property ''pad'' (type ''help MTPointSpectrum'' for details).');
			end
		case 'show',
			show = varargin{i+1};
			if ~isstring_FMAT(show,'on','off'),
				error('Incorrect value for property ''show'' (type ''help MTPointSpectrum'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help MTPointSpectrum'' for details).']);
	end
	if ~strcmp(varargin{i},'show'), v = {v{:},varargin{i:i+1}}; end
end

% Compute point spectrogram and moments
[spectrogram,~,f] = MTPointSpectrogram(lfp,v{:});
spectrogram = spectrogram';
spectrum = mean(spectrogram);
s = var(spectrogram);

% Plot log, i.e. mean E[log(spectrogram)] and stdev sqrt(Var[log(spectrogram)])
% (see http://en.wikipedia.org/wiki/Taylor_expansions_for_the_moments_of_functions_of_random_variables)
if strcmp(lower(show),'on'),
	figure;
	%logSpectrum = log(spectrum);         % classic value
	logSpectrum = log(mu)-v./(2*mu.*mu);  % corrected value, valid for mu/s > 1.5
	logS = sqrt(v./mu.^2);               % valid for mu/s > 2.5
	PlotMean(f,logSpectrum,logSpectrum-logS,logSpectrum+logS,':');
	xlabel('Frequency (Hz)');
	ylabel('Power');
	title('Power Spectrum');
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%      Copyright (C) 2009-2013 Biomedizinische NMR Forschungs GmbH
%                     http://www.biomednmr.mpg.de
%            Author: Tilman Johannes Sumpf <tsumpf@gwdg.de>
%
%       Distributed under the Boost Software License, Version 1.2.
%          (See accompanying file LICENSE_1_0.txt or copy at
%                 http://www.boost.org/LICENSE_1_0.txt)

classdef asIconClass

    properties (SetAccess = private, GetAccess = public)
        asBrowse
        colorbar
        dontSend
        download
        lineup
        magnify
        pause
        play
        refresh
        rotLeft
        rotRight
        send
        squeeze
        upload
        wsObj       % assign workspace object
    end

    methods (Access = public)
        function obj = asIconClass(iconPath)

            % determine property names of this class
            iconNames = properties(mfilename);
            iconPaths = cellfun(@(x)[iconPath,filesep,x,'.png'],iconNames,...
                'UniformOutput',false);

            % load pngs
            for i = 1 : length(iconNames)
                obj.(iconNames{i}) = iconRead(iconPaths{i});
            end
        end
    end
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [nudos]=sap2000nudos(filenam)
%programa registrado por MECAVIL ingenieria-ciencia
%bresler_lin@hotmail.com
%carga los nudos del archivo sap2000
%genera cuatros archivos txt
%[filenam pathnam]=recibetodo
%%%%%%%%%%%%%%%%%%           SAP2000 v14           %%%%%%%%%%%%%%%%%%%%%
%pathnam
%filenam
arc=fopen(filenam,'r');
arc;
encuentra =0;
while encuentra==0;
    linea=fgets(arc);
    if not(isempty(findstr('JOINT COORDINATES',linea)));
        encuentra=1;
   end
end
 nudos=[0;0;0;0];
 linea=fgets(arc);
 f=1;
 encuentra=0;
 while encuentra==0;
     [s1 s2 s3 s4 s5 s6 s7] = strread(linea,'%s %s %s %s %s %s %s');
      s1=char(s1);s4=char(s4);s5=char(s5);s6=char(s6);
      v=[s1,s4,s5,s6];
 [e1 e2]=strread(v,'%s %f ', 'delimiter',' =')  ;
 nudos=[nudos e2(:,1)];
     f=f+1;
     linea=fgets(arc);
     linea;
     if not (isempty(findstr('JOINT PATTERN DEFINITIONS',linea)));
         encuentra=1;
     end
 end
 fclose(arc);
  nudos(:,1)=[];
  nudos=nudos';
xx=nudos(:,2);
yy=nudos(:,3);
zz=nudos(:,4);
%    class(nudos)
save nudos.txt nudos -ASCII
%save quimicaXXX.txt xx -ASCII
%save quimicaYYY.txt yy -ASCII
%save quimicaZZZ.txt zz -ASCII
%  save (nudos, 'nudos', '-ascii')

%%%%%%%%%%%%%%%%%%           SAP2000 v14           %%%%%%%%%%%%%%%%%%%%%


text:
function radians = dms_to_radians ( degrees, minutes, seconds )

%*****************************************************************************80
%
%% DMS_TO_RADIANS converts an angle from degrees/minutes/seconds to radians.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DEGREES, MINUTES, SECONDS, an angle in degrees, minutes,
%    and seconds.
%
%    Output, real RADIANS, the equivalent angle in radians.
%
  angle =   degrees + ( minutes + ( seconds / 60.0 ) ) / 60.0;

  radians = ( angle / 180.0 ) * pi;

  return
end


text:
function DEMO_SLR
% Demo of sparse logistic regression. This demonstration routines
% illustrates the use of sparse acoustic regression (as implemented in
% spm_sparse_regression.m) to recover a small number of causes that best
% explain some dependent variable. For example, imagine we thought that a
% small number (e.g., four) of SNPs or copy number variants a generic study
% had sufficiently large effect sizes on some phenotypic measure to be
% worth pursuing. We might treat these (rare variants) as being expressed
% in the context of random effects (e.g., common variants and phenotypic
% measurement noise); however, we have many more potential causes than
% observations. This problem is addressed using (logistic) regression
% under sparsity constraints specified in terms of hyper priors over the
% precision (i.e. inverse variance) of model parameters. This provides the
% Bayesian shrinkage estimators of the regression coefficients that,
% crucially, can then be subject to Bayesian model reduction. Bayesian
% model reduction effectively eliminates redundant parameters that provided
% the optimal balance between accuracy and complexity.
%
% in the example below, we assume that we have 32 subjects with 128
% independent variables (e.g., following some initial dimension reduction).
% The simulated data is generated with just four of the independent to see
% whether these can be identified using sparse logistic regression and
% Bayesian model reduction.
%
% If the dependent variables are classes or probabilities a logistic
% transform is automatically applied. However, one can also use this
% routine for continuous (i.e., parametric phenotypes). the graphics
% produced by this demo report the results of sparse logistic regression
% using variational Laplace (i.e., approximate Bayesian inference and hyper
% priors). In addition, it reports the results and summary of the
% subsequent Bayesian model reduction.
%
% see also: spm_sparse_regression.m
%__________________________________________________________________________
% Copyright (C) 2014 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: DEMO_SLR.m 7454 2018-10-19 19:38:50Z karl $

% Genomic data (G)
%--------------------------------------------------------------------------
Ny   = 32;                                % number of observations
Nx   = 128;                               % number of explanatory variables

% Design matrix (U)
%--------------------------------------------------------------------------
X    = randn(Ny,Nx);
X0   = ones(Ny,1);

% Relevant explanatory variables (j)
%--------------------------------------------------------------------------
j    = 1:4;
P    = zeros(Nx,1);
P(j) = randn(length(j),1);

% Simulate (classification probability) data with a softmax
%--------------------------------------------------------------------------
y    = X*P + randn(Ny,1)/8;
y    = spm_softmax(y);

% Bayesian model reduction
%==========================================================================
RCM  = spm_sparse_regression(y,X,X0);

% show results
%--------------------------------------------------------------------------
spm_figure('Getwin','Model posterior (over families)'); clf

subplot(2,2,1)
spm_plot_ci(RCM.Ep(j),diag(RCM.Cp(j,j)));
xlabel('parameter')
title('Posterior density','FontSize',16)
axis square

subplot(2,2,2)
bar(P(j));
xlabel('parameter')
title('True (nonzero) parameters','FontSize',16)
axis square

subplot(2,2,3)
bar(RCM.Ep)
title('Recovered (sparse) profile','FontSize',16)
xlabel('parameter')
axis square, a = axis;

subplot(2,2,4)
bar(P)
title('True profile','FontSize',16)
xlabel('parameter')
axis square, axis(a)

author: spm


text:
function headmodel = triangle4pt(headmodel)

% TRIANGLE4PNT takes the volume model and estimates the 4th point of each
% triangle of each mesh.
%
% Use as
%   headmodel = triangle4pt(headmodel)
%
% In each headmodel.bnd sub-structure, a field '.pnt4' is added. The '.pnt4'
% field is a Ntri*3 matrix, with the coordinates of a point for each
% triangle in the meshed surface.
%
% Explanations:
% The point is that for some BEM, specifically 'solid angle', calculation
% it is necessary to estimate the local curvature of the true surface which
% is approximated by the flat triangle. One way to proceed is to use
% "close by" vertices to estimate the overall area's curvature.
% A more elegant(?) way uses a 4th point for each triangle: the "centroid"
% of the triangle is simply pusehd away from the triangle surface to fix
% the local surface curvature (assuming the surface is smooth enough).
% This 4th point is thus hovering above/under the triangle and can be used
% to fit a sphere on the triangle in a realistic way.
%
% Method:
% - The 4th point can/could be defined at the tessalation stage, based on
%   the anatomical images directly.
% - With any model, the curvature can be estimated/approximated by looking
%   at the vertices around the triangle considered and fit a sphere on
%   those few vertices, assuming the surface is smooth enough
% The latter option is the one followed here.
% The extra-vertices considered here are those 3 which are linked to the
% triangle by 2 edges.
%__________________________________________________________________________
%
% written by Christophe Phillips, 2009/01/19
% Cyclotron Research Centre, University of li?ge, belgium
%
% $Id$

Ns = length(headmodel.bnd);
for ii=1:Ns % treat each mesh one at a time
  tri = headmodel.bnd(ii).tri;
  pos = headmodel.bnd(ii).pos;
  Nt = size(tri,1);
  pnt4 = zeros(Nt,3);
  for jj=1:Nt % treat each triangle on a t a time
    lt1 = find(tri(:,1)==tri(jj,1) | tri(:,2)==tri(jj,1) | ...
      tri(:,3)==tri(jj,1));
    lt2 = find(tri(:,1)==tri(jj,2) | tri(:,2)==tri(jj,2) | ...
      tri(:,3)==tri(jj,2));
    lt3 = find(tri(:,1)==tri(jj,3) | tri(:,2)==tri(jj,3) | ...
      tri(:,3)==tri(jj,3));
    lt = [intersect(lt1,lt2); intersect(lt2,lt3); intersect(lt3,lt1)];
    lt(lt==jj) = [];
    % list of 3 directly surrounding triangles
    lv = tri(lt,:);
    lv = setxor(lv(:)',tri(jj,:));
    % list of 3 voxels connected by 2 edges to the jj_th triangle.
    sph_pnt = pos([tri(jj,:) lv],:);
    [center,radius] = fitsphere(sph_pnt);
    % best fitting sphere radius & centre, for the 6 points chosen
    pos_c = sum(pos(tri(jj,:),:))/3;
    % centroid of the triangle treated
    if isfinite(radius)
      tmp = pos_c-center;
      vd = tmp/norm(tmp);
      % unit vector from sphere center in direction of middle triangle
      pnt4(jj,:) = center+vd*radius;
      % projection of the centroid, at 'radius' from the center
    else
      % radius at infinite, i.e. flat surface
      pnt4(jj,:) = pos_c;
    end
  end
  headmodel.bnd(ii).pnt4 = pnt4;
end

return

% figure, plot3(pos(:,1),pos(:,2),pos(:,3),'.')
% axis equal, hold on
% plot3(pos(tri(jj,:),1),pos(tri(jj,:),2),pos(tri(jj,:),3),'*')
% plot3(pos(lv,1),pos(lv,2),pos(lv,3),'o')
% plot3(pnt4(jj,1),pnt4(jj,2),pnt4(jj,3),'rs')

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [j, dist] = LMgistquery(gistQuery, gist, IsCorrelation)

if exist('IsCorrelation','var')~=1
    IsCorrelation = 0;
end

if IsCorrelation
    % Normalized correlation
    gistQuery = normalizeGist(gistQuery);
    gist = normalizeGist(gist);
    dist = 2-2*gist*gistQuery';
    [dist,j] = sort(dist);
else
    % L2
    dist = sum((gist - repmat(gistQuery, [size(gist,1) 1])).^2,2);
    [dist,j] = sort(dist);
end

author: CSAILVision

repo: LabelMeToolbox


text:
function [T,totalError,exitCode]=LotkaVolterraPeriod(x0,y0,a,b,c,d,intParams)
%%LOTKAVOLTERRAPERIOD Determine the period of a given set of Lotka-Volterra
%           equations. These are also known as predator-prey equations.
%           These are the following set of two coupled differential
%           equations:
%           dx=a*x-b*x*y
%           dy=c*x*v-d*y
%           with a,b,c,d>0 and x,y with initial conditions >0. Such
%           equations are useful for analyzing the performance of
%           techniques for solving initial values problems.
%
%INPUTS: x0, y0 The real initial conditions of the two coupled processes.
%               x0,y0 >0.
%       a,b,c,d The real coefficients of the model, all >0.
%     intParams Optionally, one can pass a structure affecting how the
%               numerical integration in this function is performed. This
%               can be a structure with fields RelTol, AbsTol, n, and
%               maxSearchReg having respective default values of 1e-8,
%               1e-11, 8, and 1000. These values correspond to the
%               same-named parameters of the integral1DAdaptive function.
%
%OUTPUTS: T The period of the Lotka-Volterra system.
% totalError An estimate of the total error of the solution. This is
%           derived from the same outpuf of the integral1DAdaptive
%           function.
%  exitCode The exit code returned by the integral1DAdaptive function. 0
%           means that the absolute ore relative error condition was
%           fulfilled and thus, the algorithm converged.
%
%This function implements the integral given by Theorem 6.2 of [1].
%
%EXAMPLE:
%Here, we determine the period of s system and then we perform numerical
%integration over that period to verify that the solution looks like a
%closed loop.
% a=2/3;
% b=4/3;
% c=1;
% d=1;
% x0=1/5;
% y0=1/5;
% [T,totalError,exitCode]=LotkaVolterraPeriod(x0,y0,a,b,c,d)
% %The Lotka-Volterra differential system.
% f=@(x,t)([a;-d].*x+[-b;c]*prod(x));
% numSteps=1e4;
% xStart=[x0;y0];
% theTimes=linspace(0,T,numSteps);
% xList=RungeKAtTimes(xStart,theTimes,f,Inf);
% figure(1)
% clf
% plot(xList(1,:),xList(2,:),'linewidth',2)
% h1=xlabel('x');
% h2=ylabel('y');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%One will see that the plot is a closedl oop and that reducing the
%integration time a little bit opens the loop, indicating that the solution
%was correct.
%
%REFERENCES:
%[1] S.-D. Shih, "The period of a Lotka-Volterra system," Taiwanese Journal
%    of Mathematics, vol. 1, no. 4, pp. 451-470, Dec. 1997.
%
%November 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

RelTol=1e-8;
AbsTol=1e-11;
n=8;
maxSearchReg=1000;

if(nargin>6&&~isempty(intParams))
    if(isfield(intParams,'RelTol'))
        RelTol=intParams.RelTol;
    end

    if(isfield(intParams,'AbsTol'))
        AbsTol=intParams.AbsTol;
    end

    if(isfield(intParams,'n'))
        n=intParams.n;
    end

    if(isfield(intParams,'maxSearchReg'))
        n=intParams.maxSearchReg;
    end
end

%The energy of the system, Equation 1.5. E>=0.
E=max(0,c*x0-d+b*y0-a-a*log((b/a)*y0)-d*log((c/d)*x0));

func=@(s)intArg(s,E,a,d);

[intVal,totalError,exitCode]=integral1DAdaptive(func,[0;E],n,2,[],RelTol,AbsTol,maxSearchReg);
T=(1/(a*d))*intVal;
totalError=(1/(a*d))*totalError;
end

function val=intArg(s,E,a,d)
%%INTARG The argument fo the integral in Equation 6.6 of [1].
%
%REFERENCES:
%[1] S.-D. Shih, "The period of a Lotka-Volterra system," Taiwanese Journal
%    of Mathematics, vol. 1, no. 4, pp. 451-470, Dec. 1997.
%
%November 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
val=Phi(s/d).*Phi((E-s)/a);
val(~isfinite(val))=0;%If the endpoints were evaluated.

end

function val=Phi(s)
%%PHI The Phi function of Equation 6.7 of [1].
%
%REFERENCES:
%[1] S.-D. Shih, "The period of a Lotka-Volterra system," Taiwanese Journal
%    of Mathematics, vol. 1, no. 4, pp. 451-470, Dec. 1997.
%
%November 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

z=-exp(-1-s);
val=1./(1+LambW(z,0))-1./(1+LambW(z,-1));
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 cmap = hsvCmap(numGrays,numColors)
%
% cmap = hsvCmap(numGrays,numColors)
%
% Makes colormap array with:
%   gray scale - 1:numGrays
%   hsv colors - numGrays+1:numGrays+numColors
%
% djh 1/98

if ~exist('numGrays','var')
  numGrays=128;
end
if ~exist('numColors','var')
  numColors=96;
end

cmap = zeros(numGrays+numColors,3);
cmap(1:numGrays+numColors,:) = [gray(numGrays);hsv(numColors)];

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% PLANEWAVEMAXWELL plane wave solutions to Maxwell equations in a cube.
%
% Test Maxwell function can solve problems with complex solution and
% indefinite case.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

close all;

%% Generate an initial mesh
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],1);
bdFlag = setboundary3(node,elem,'Neumann');

%% Parameters
maxIt = 4;
N = zeros(maxIt,1);
h = zeros(maxIt,1);
energyErr = zeros(maxIt,1);
L2Err = zeros(maxIt,1);
energyErrImag = zeros(maxIt,1);
L2ErrImag = zeros(maxIt,1);

%% Get the data of the pde
global d P omega
d = [1, pi/2, pi/2];  % in sphereical coordinate
P = [1, 0, 0];
omega = 1;
r = d(1); theta = d(2); phi = d(3);
d = [r*sin(phi)*cos(theta), r*sin(phi)*sin(theta), r*cos(phi)];
% pde = planewavedataC; % plane wave with complex coefficients
pde = planewavedata; % plane wave with real coefficients and complex solution

%% Finite Element Method
for k = 1:maxIt
    % refine grid
    [node,elem,bdFlag] = uniformrefine3(node,elem,bdFlag);
    % solve the equation
    [u,edge,eqn] = Maxwell(node,elem,bdFlag,pde);
    % compute the error
    uI = edgeinterpolate(pde.exactu,node,edge);
    L2Err(k) = sqrt(abs(real(u-uI)'*eqn.M*real(u-uI)));
    energyErr(k) = sqrt(abs(real(u-uI)'*eqn.A*real(u-uI)) + L2Err(k)^2);
    L2ErrImag(k) = sqrt(abs(imag(u-uI)'*eqn.M*imag(u-uI)));
    energyErrImag(k) = sqrt(abs(imag(u-uI)'*eqn.A*imag(u-uI)) + L2ErrImag(k)^2);
%     energyErr(k) = getHcurlerror3ND(node,elem,pde.curlu,u);
%     L2Err(k) = getL2error3ND(node,elem,pde.exactu,u);
    N(k) = length(u);
    h(k) = 1./(size(node,1)^(1/3)-1);
end

%% Plot convergence rates
figure(1); clf;
showrateh2(h,energyErr,1,'r-+','|| Re(u-u_h)||_A',...
           h,L2Err,1,'b-+','|| Re(u-u_h)||');
figure(2); clf;
showrateh2(h,energyErrImag,1,'r-+','|| Img(u-u_h)||_A',...
           h,L2ErrImag,1,'b-+','|| Img(u-u_h)||');

repo: ifem