text:
classdef ThreeDimKspWindTunnelDragCoeffModel < AbstractDragCoefficientModel
    %ThreeDimKspWindTunnelDragCoeffModel Summary of this class goes here
    %   Detailed explanation goes here

    properties
        slices(1,:) ThreeDimKWTDragCoefficientModelSlice

        speed(1,:) double
        altitude(:,1) double
        aoa(1,:) double
        data double

        gi
    end

    properties(Constant)
        enum = DragCoefficientModelEnum.ThreeDimWindTunnel
    end

    methods
        function obj = ThreeDimKspWindTunnelDragCoeffModel(dataSlices)
            obj.slices = dataSlices;
            obj.sortSlicesByAoA();
        end

        function CdA = getDragCoeff(obj, ut, rVect, vVect, bodyInfo, mass, altitude, pressure, density, vVectECEFMag, totalAoA, aoa, sideslip)
            arguments
                obj(1,1) ThreeDimKspWindTunnelDragCoeffModel
                ut(1,1) double
                rVect(3,1) double
                vVect(3,1) double
                bodyInfo(1,1) KSPTOT_BodyInfo
                mass(1,1) double
                altitude(1,1) double
                pressure(1,1) double
                density(1,1) double
                vVectECEFMag(1,1) double
                totalAoA(1,1) double = 0;
                aoa(1,1) double = 0;
                sideslip(1,1) double = 0;
            end

            CdA = obj.gi(altitude, vVectECEFMag, totalAoA);
            CdA = max([CdA, 0]);
        end

        function addSlice(obj, newSlice)
            obj.slices(end+1) = newSlice;
            obj.sortSlicesByAoA();
        end

        function removeSlice(obj, slice)
            obj.slices(obj.slices == slice) = [];
            obj.sortSlicesByAoA();
        end

        function tf = usesTotalAoA(obj)
            tf = true;
        end

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

        function createGriddedInterpFromSlices(obj)
            obj.sortSlicesByAoA();

            obj.speed = obj.slices(1).twoDSlice.speed;
            obj.altitude = obj.slices(1).twoDSlice.altitude;
            obj.aoa = [obj.slices.aoa];

            obj.data = [];
            for(i=1:length(obj.slices)) %#ok<*NO4LP>
                slice = obj.slices(i);
                obj.data(:,:,i) = slice.twoDSlice.data;
            end

            gridVecs = {obj.altitude, obj.speed, obj.aoa};
            obj.gi = griddedInterpolant(gridVecs, obj.data, 'makima', 'nearest');
        end

        function sortSlicesByAoA(obj)
            AoAs = [obj.slices.aoa];
            [~,I] = sort(AoAs);
            obj.slices = obj.slices(I);
        end

        function [sliceGoodTf, messages] = getStatusOfSlices(obj)
            sliceGoodTf = true(1, length(obj.slices));
            messages = cell(1, length(obj.slices));

            AoAs = [obj.slices.aoa];

            for(i=1:length(obj.slices))
                slice = obj.slices(i);

                [thisSliceGoodTf, sliceMessages] = slice.getSliceStatus();
                sliceGoodTf(i) = sliceGoodTf(i) && thisSliceGoodTf;
                messages{i} = horzcat(messages{i}, sliceMessages);

                if(length(obj.slices) == 1)
                    sliceGoodTf(i) = false;
                    messages{i}(end+1) = 'There must be at least two AoA data sets.';
                end

                if(i >= 2)
                    subAoAs = AoAs;
                    subAoAs(i) = [];
                    if(any(AoAs(i) == subAoAs))
                        sliceGoodTf(i) = false;
                        messages{i}(end+1) = 'Duplicate angle of attack data sets are not allowed.';
                    end

                    slice1 = obj.slices(1);
                    twoDModelSlice1 = slice1.twoDSlice;

                    altitudeSlice1 = twoDModelSlice1.altitude;
                    speedSlice1 = twoDModelSlice1.speed;

                    twoDModel = slice.twoDSlice;
                    altitudeSliceI = twoDModel.altitude;
                    speedSliceI = twoDModel.speed;

                    if(numel(altitudeSliceI) == numel(altitudeSlice1))
                        bool = altitudeSliceI(:) - altitudeSlice1(:) ~= 0;
                        if(any(bool))
                            sliceGoodTf(i) = false;
                            messages{i}(end+1) = "The altitude axis points do not match those of the first AoA data set.  All axis points must be identical in all data sets.";
                        end
                    else
                        sliceGoodTf(i) = false;
                        messages{i}(end+1) = "There are not the same number of altitude axis points of those of the first AoA data set.  All axis points must be identical in all data sets.";
                    end

                    if(numel(speedSliceI) == numel(speedSlice1))
                        bool = speedSliceI(:) - speedSlice1(:) ~= 0;
                        if(any(bool))
                            sliceGoodTf(i) = false;
                            messages{i}(end+1) = "The speed axis points do not match those of the first AoA data set.  All axis points must be identical in all data sets.";
                        end
                    else
                        sliceGoodTf(i) = false;
                        messages{i}(end+1) = "There are not the same number of speed axis points of those of the first AoA data set.  All axis points must be identical in all data sets.";
                    end
                end

                if(sliceGoodTf(i) == true)
                    messages{i} = "No warnings or errors found.";
                end
            end
        end

        function [listBoxStr, slices] = getListBoxStrs(obj)
            for(i=1:length(obj.slices))
                listBoxStr(i) = string(obj.slices(i).getListboxStr()); %#ok<AGROW>
            end

            slices = obj.slices;
        end

        function useTf = openEditDialog(obj, lvdData)
            out = AppDesignerGUIOutput({false});
            lvd_Edit3DKerbalWindTunnelDragPropertiesGUI_App(obj, lvdData, out);
            useTf = out.output{1};
        end
    end

    methods(Static)
        function obj = loadobj(obj)
            if(isempty(obj.slices))
                slice1 = ThreeDimKspWindTunnelDragCoeffModel.produceDataSlice(0, '');
                obj.addSlice(slice1);

                slice2 = ThreeDimKspWindTunnelDragCoeffModel.produceDataSlice(deg2rad(1), '');
                obj.addSlice(slice2);
            end

            try
                obj.createGriddedInterpFromSlices();
            catch ME
                % nothing to do here
            end
        end

        function slice = produceDataSlice(aoa, sliceDataFile)
            slice = ThreeDimKWTDragCoefficientModelSlice(aoa, sliceDataFile);
        end

        function dragCoeffModel = createDefault()
            slice1 = ThreeDimKspWindTunnelDragCoeffModel.produceDataSlice(0, '');
            slice2 = ThreeDimKspWindTunnelDragCoeffModel.produceDataSlice(deg2rad(1), '');
            slices = [slice1, slice2];

            dragCoeffModel = ThreeDimKspWindTunnelDragCoeffModel(slices);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
 function tester_tomo2(A1, mask, varargin)
%function tester_tomo2(A1, mask, [options])
%|
%| Test suite for a Fatrix-type 2D or 3D system object, including Gblock tests.
%|
%| in
%|	A1	Fatrix or fatrix2
%|	mask	logical [nx ny] or [nx ny nz]
%|
%| option
%|	'A2'	Fatrix	optional 2nd object for testing and comparisons
%|			(it is also put through all the same tests)
%|	'multi'	0|1	see Fatrix_test_basic()
%|	'halt'	0|1	see Fatrix_test_basic()
%|	'equiv_thresh'	instead of jf_equal, use equivs() with this threshold
%|			default: [], which means use jf_equal
%|	'nblock'	# of blocks for testing block version (default 2)
%|			set to 0 to prevent block test
%|
%| Copyright 2005-8-2, Jeff Fessler, University of Michigan

if nargin < 2, ir_usage, end

arg.multi = true;
arg.A2 = [];
arg.halt = true;
arg.equiv_thresh = [];
arg.nblock = 2; % for block test
arg = vararg_pair(arg, varargin, 'subs', {'G2', 'A2'});

Fatrix_test_basic(A1, mask, 'multi', arg.multi, ...
	'halt', arg.halt, ...
	'name', inputname(1), 'caller', caller_name)

switch ndims(mask)
case 2
	[nx, ny] = size(mask);
	ig = image_geom('nx', nx, 'ny', ny, 'dx', 1);
	x = ellipse_im(ig, []);

case 3
	[nx, ny, nz] = size(mask);
	ig = image_geom('nx', nx, 'ny', ny, 'nz', nz, 'dx', 1);
	x = ellipsoid_im(ig, '');
	x = x .* mask;
end

if arg.nblock
	tester_tomo2_block(A1, mask, x, arg.nblock, arg.multi, arg.equiv_thresh)
end

if ~isempty(arg.A2)
	Fatrix_test_basic(arg.A2, mask, 'multi', arg.multi, ...
		'name', inputname(4), 'caller', caller_name)
	tester_tomo2_compare(A1, arg.A2, x)
end

%
% tester_tomo2_compare()
%
function tester_tomo2_compare(A1, A2, x)

y1 = A1 * x;
y2 = A2 * x;
my_compare(y1, y2, 'A*x')

x1 = A1' * y1;
x2 = A2' * y1;
equivs(x1, x2)

j = round(size(A1,2) / 2); % roughly a middle pixel
y1 = A1(:,[j j+1]);
y2 = A2(:,[j j+1]);
y1 = [y1(:,1) y1(:,2)]; % for fatrix2
y2 = [y2(:,1) y2(:,2)]; % for fatrix2
my_compare(y1, y2, 'A(:,j)')

% check A(:,:)
if 0 && size(x,1) < 100
	t1 = A1(:,:);
	t2 = A2(:,:);
	my_compare(t1, t2, '(:,:)');
%	mpd = max_percent_diff(t1,t2);
%	printf('A(:,:)	mpd %g', mpd)
%	if mpd/100 > 1e-6, error 'A(:,:)', end
end

printm('passed %s', A1.caller)

%
% tester_tomo2_block()
% now block version
%
function tester_tomo2_block(A1, mask, x, nblock, multi, equiv_thresh)

B1 = Gblock(A1, nblock);

% B * x
y1 = A1 * x;
y2 = B1 * x;
my_compare(y1, y2, 'B*x')

y0 = y1;
odim = A1.arg.odim;
na = odim(end);
switch length(odim)
case 3
	getp = @(x,i) x(:,:,i);
case 2
	getp = @(x,i) x(:,i);
otherwise
	fail 'not done'
end
% caution: handle single view 3D projectors [ns nt 1 nrep]
catp = @(a,b) cat(length(odim)+1, a, b);

% B' * y
x1 = A1' * y0;
x2 = B1' * y0;
my_compare(x1, x2, 'B''*y')

comp = @(t1, t2, str) my_compare2(t1, t2, str, equiv_thresh);

%
% block operations
%
for k=1:nblock
	ia = k:nblock:na;
	str = sprintf('B{%d}', k);

	% check B{k} * x
	t1 = A1 * x;
	t1 = getp(t1,ia);
	t2 = B1{k} * x;
	comp(t1, t2, 'B{k} * x')

	% B{k} * [x x]
	if multi
		t2 = B1{k} * stackup(x,x);
		comp(catp(t1, t1), t2, [str ' * [x x]'])
	end

	% check B{k} * x(mask)
	t2 = B1{k} * x(mask);
	t2 = reshape(t2, size(t1));
	comp(t1, t2, [str ' * x()'])

	% check B{k} * [x(mask) x(mask)]
	if multi
		t2 = B1{k} * [x(mask) x(mask)];
		comp([t1(:) t1(:)], t2, [str ' * [x() x()]'])
	end

	% check B{k}' * y()
	tmp = block_insert(ia, odim, getp(y0,ia));
	t1 = A1' * tmp(:);
	t2 = B1{k}' * col(getp(y0,ia));
	if isempty(equiv_thresh)
%		my_compare(t1, t2, [str ''' * y'])
		mpd = max_percent_diff(t1,t2);
		if mpd
			printf([str '''*y mpd %g'], mpd)
			if mpd/100 > 1e-6
				fail 'B{k}'' * y()'
%				keyboard
			end
		end
	else
		comp(t1, t2, [str ''' * y'])
	end
end

% my_compare()
function my_compare(t1, t2, arg)
%max_percent_diff(t1, t2, arg)
%try
	jf_equal(t1, t2)
%catch
%	warn('%s failed', arg)
%	keyboard
%end

function my_compare2(t1, t2, arg, equiv_thresh)
if ~isempty(equiv_thresh)
	equivs(t1, t2, 'thresh', equiv_thresh)
else
	try
		my_compare(t1, t2, arg)
	catch
	%	keyboard
		fail('testing %s', arg)
	end
end

% block_insert()
% for block or OS methods, make data array that is all zeros but at "ia"
function out = block_insert(ia, dims, data)
out = zeros(dims);
switch length(dims)
case 2
	out(:,ia) = data;
case 3
	out(:,:,ia) = data;
otherwise
	fail 'not done'
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function u = nvecs(t,n,r,opts)
%NVECS Compute the leading mode-n vectors for a sparse tensor.
%
%   U = NVECS(X,n,r) computes the r leading eigenvalues of Xn*Xn'
%   (where Xn is the mode-n matricization of X), which provides
%   information about the mode-n fibers. In two-dimensions, the r
%   leading mode-1 vectors are the same as the r left singular vectors
%   and the r leading mode-2 vectors are the same as the r right
%   singular vectors.
%
%   U = NVECS(X,n,r,OPTS) specifies options:
%   OPTS.eigsopts: options passed to the EIGS routine [struct('disp',0)]
%   OPTS.flipsign: make each column's largest element positive [true]
%
%   See also SPTENSOR, SPTENMAT, EIGS.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if ~exist('opts','var')
    opts = struct;
end

if isfield(opts,'eigsopts')
    eigsopts = opts.eigsopts;
else
    eigsopts.disp = 0;
end

tnt = double(sptenmat(t,n,'t'));
y = tnt' * tnt;
opts.disp = 0;
[u,d] = eigs(y,r,'LM',eigsopts);

%tn = sptenmat(t,n);
%[u,d] = eigs(@(x)aatx(tn,x), size(t,n), r, 'LM', eigsopts);

if isfield(opts,'flipsign')
    flipsign = opts.flipsign;
else
    flipsign = true;
end

if flipsign
    % Make the largest magnitude element be positive
    [val,loc] = max(abs(u));
    for i = 1:r
        if u(loc(i),i) < 0
            u(:,i) = u(:,i) * -1;
        end
    end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function out = CO_PartialAutoCorr(y,maxTau,whatMethod)
% CO_PartialAutoCorr   Compute the partial autocorrelation of an input time series
%
%---INPUTS:
% y, a scalar time series column vector.
%
% maxTau, the maximum time-delay. Returns for lags up to this maximum.
%
% whatMethod, the method used to compute: 'ols' or 'yule_walker'
%
%---OUTPUT: the partial autocorrelations across the set of time lags.
%

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function 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/>.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%% Check inputs and set defaults:
% ------------------------------------------------------------------------------
if nargin < 2
    % Use a maximum lag of 10 by default
    maxTau = 10;
end

if nargin < 3 || isempty(whatMethod)
    % ordinary least square by default
    whatMethod = 'ols';
end

%-------------------------------------------------------------------------------
%% Initial checks on maxTau
%-------------------------------------------------------------------------------
N = length(y); % time-series length

assert(maxTau > 0)

if maxTau < 0
    error('Negative time lags not applicable')
end

% ------------------------------------------------------------------------------
%% Do the computation
% ------------------------------------------------------------------------------

pacf = parcorr(y,'NumLags',maxTau,'Method',whatMethod);

% Zero lag is the first entry in the PACF (and should always be 1)

for i = 1:maxTau
    out.(sprintf('pac_%u',i)) = pacf(i+1);
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [ eqtc, varargout ] = import_events(dbname, chanfile, starttime, endtime, phase, delay, enddelay, varargin)
%IMPORT_EVENTS Import waveforms associated with earthquake events.
% EQTC = IMPORT_EVENTS(DBNAME, CHANFILE, STARTTIME, ENDTIME, ...
%                      PHASE, DELAY, ENDDELAY) returns a cell array of
% threecomp objects, waveform objects, and empty cells, containing
% waveforms organized column-wise by station and row-wise by event.  Input
% arguments are as follows:
%
%    DBNAME: Filename of an Antelope database.
%    CHANFILE: Filename of a comma-delimited list of station-channel
%    associations (see IE_READCHANS helpfile for the format).
%    STARTTIME: Start date/time in Antelope-readable string format.
%    ENDTIME: End date/time in Antelope-readable string format.
%    PHASE: Arrival phase to use (e.g. 'P', 'S')
%    DELAY: Seconds preceding arrival time to retrieve into waveforms.
%    ENDDELAY: Seconds following arrival time to retrieve into waveforms.
%
% The intent of this function is to import waveform data from multiple
% events and stations, and to store it in an organized cell array which
% preserves associations between arrivals.  Three-component data is
% retrieved if it is available, and converted to threecomp objects.
% Single-component data remains as waveform objects.  Values ETYPE, ORID,
% OTIME, SEAZ, and STYPE are also retrieved from the database and stored as
% custom fields within each waveform object.  If there is no arrival for an
% event at a station, that cell is left empty.  If three-component data is
% available but cannot be converted into a threecomp object (usually
% indicating data problems), the cell is also left empty.
%
% [EQTC, EQWF] = IMPORT_EVENTS(...) returns the aforementioned cell array
% in EQTC as well as a cell array of unadulterated waveform objects in
% EQWF.  These will be 1x1 waveforms for single-component data and 1x3
% waveforms for three-component data (including those that could not be
% converted to threecomp objects and are not present in EQTC).
%
% ... = IMPORT_EVENTS(..., DS) uses specified GISMO datasource DS to
% retrieve the waveform data.  The earthquake catalogue and phase arrivals
% are always taken from the Antelope database DBNAME, but actual waveform
% data can be in any GISMO-readable format.  Examples: DS =
% datasource('uaf_continuous') or DS =
% datasource('winston','myserver.uaf.edu',12345).
%
% RETURN VALUE EXAMPLES
%
% In both EQTC and EQWF, waveforms are organized into columns by station
% and rows by event.  Here is a hypothetical three-station five-earthquake
% example; stations EG1 and EG2 are three-component, station EG3 is
% single-component.
%
% EQTC =  % EG1                EG2                EG3
%                  []    [1x1 threecomp]    [1x1 waveform]   % Event 1
%     [1x1 threecomp]    [1x1 threecomp]                []   % Event 2
%                  []                 []    [1x1 waveform]   % Event 3
%     [1x1 threecomp]    [1x1 threecomp]    [1x1 waveform]   % Event 4
%                  []                 []    [1x1 waveform]   % Event 5
%
% EQWF =  % EG1                EG2                EG3
%                  []     [1x3 waveform]    [1x1 waveform]   % Event 1
%      [1x3 waveform]     [1x3 waveform]                []   % Event 2
%                  []                 []    [1x1 waveform]   % Event 3
%      [1x3 waveform]     [1x3 waveform]    [1x1 waveform]   % Event 4
%      [1x3 waveform]                 []    [1x1 waveform]   % Event 5
%
% Notice that EQTC contains both threecomp objects and waveform objects.
% EQWF contains waveform objects only, but with size 1x3 for
% three-component stations and 1x1 for single-component stations.  When
% events are not recorded at a station, the cell is left empty [].
% IMPORTANT: The cell in EQTC is also left empty when a 1x3 waveform cannot
% be converted into a threecomp object, such as with event 5 at station EG1 in
% the above example.
%
% CAVEATS AND LIMITATIONS:
%
% 1. The earthquake list is taken from the origin table without
% consideration of the event table or the value of prefor.  If more than
% one origin is present for an event, they will all be retrieved.
%
% 2. The value of STYPE is taken from the chosen phase arrival, and is assumed
% to be the same across all components of three-component data.
%
% 3. Unassociated arrivals will not be retrieved.  The returned data
% structure is versatile enough that unassociated arrivals _could_ be
% stored in an intuitive way (one per row) but so far this has not been
% implemented.
%
% Author: Christopher Bruton, Geophysical Institute, University of Alaska Fairbanks
% $Date$
% $Revision$

warning('This function likely broken. It needs to be rebuilt on top of core GISMO classes for Catalog and Arrival objects')
%% Check for valid input

% dbname is checked by dbopen()

% chanfile is checked by ie_readchans()

% Convert dates.  We could pass the user input directly to Antelope in the
% query, but malformed dates could cause unexpected results.  There doesn't
% seem to be a simple way to otherwise check for valid date formats or to
% escape special characters. This way we at least ensure a valid epoch number.
starttime = str2epoch(starttime);
endtime = str2epoch(endtime);

% Set flag if dates are not given; this will prevent subset on starttime
% and endtime and instead return the whole catalogue.  Note that
% str2epoch()'s return value does not distinguish between invalid input
% (e.g. 'asdfjd') and true zero values (e.g. 'Jan 1 1970'), so this flag
% could get set inadvertently.
if (starttime == 0) && (endtime == 0)
	all_times = true;
else
	all_times = false;
	% Check date range
	if (endtime < starttime)
		error('Argument endtime must be after argument starttime.');
	end
end

% Check delay times.  These should be numbers
if ~isscalar(delay) || ~isnumeric(delay)
	error('Argument delay must be a scalar number.');
end
if ~isscalar(enddelay) || ~isnumeric(enddelay)
	error('Argument enddelay must be a scalar number.');
end

% Check that waveform time range is valid
if (enddelay + delay) <= 0
	error('Arguments delay and enddelay must denote a positive non-zero time range.');
end

% Check that phase is an alphanumeric string (perhaps this cannot be
% expected to be universally true, but the alternative would require
% checking that special characters are safe).
if ~all(isstrprop(phase,'alphanum'))
	error('Argument phase must be an alphanumeric string.');
end

% Check if optional datasource object was provided
switch size(varargin,2)
	case 0
		ds = datasource('antelope',dbname); % Standard antelope datasource
	case 1
		ds = varargin{1}; % User-supplied datasource
	otherwise
		error('Incorrect number of input arguments.');
end

%% Check output arguments

switch nargout
	case 1
		ret_eqwf = false;
	case 2
		ret_eqwf = true;
	otherwise
		error('Incorrect number of output arguments.');
end

%% Read and verify list of stations and channels

fprintf(2,'Loading and verifying station/channel list...\n');

% Open list
chanlist = ie_readchans(chanfile);

% The above function does some basic error checking but we still need to
% confirm things with the database (sitechan table)

% Open DB, read sitechan table, and close
db = dbopen(dbname,'r');
db = dblookup_table(db,'sitechan');
db = dbsort(db,'sta','chan');
[tsta,tchan] = dbgetv(db,'sta','chan');
dbclose(db);

% Now check that every station-channel pair in chanlist exists
stations = chanlist(:,1);
for n = 1:length(stations)
	channels = chanlist(n,2:4);
	channels(strcmp('',channels)) = []; % Delete empty elements

	tchans = tchan(strcmp(stations(n),tsta)); % Channels from tchan corresponding to station n
	for m = 1:length(channels)
		if ~any(strcmp(channels(m),tchans))
			error(['Channel ' channels{m} ' at station ' stations{n} ' does not exist in sitechan table.']);
		end
	end
end

%% Load the earthquake catalogue

fprintf(2,'Loading earthquake catalogue...');

% Open database
db = dbopen(dbname,'r');

% Get all arrivals (of desired phase) for origins within time period
db = dblookup_table(db,'origin');
if all_times == false
	db = dbsubset(db,['time >= ' num2str(starttime,'%f') ' && time < ' num2str(endtime,'%f')]);
end
db1 = dblookup_table(db,'assoc');
db1 = dbsubset(db1,['phase == "' phase '"']);
db = dbjoin(db,db1);

% Create station subset string
stastring = ['sta=="' stations{1} '"'];
if length(stations) > 1
	for n = 1:length(stations)
		stastring = [stastring ' || sta=="' stations{n} '"'];
	end
end

% Subset to only chosen stations
db = dbsubset(db,stastring);

% Join to arrival table
db1 = dblookup_table(db,'arrival');
db = dbjoin(db,db1);

% Sort by origin time and id
db = dbsort(db,'origin.time','orid');

% Get the values
[sta,atime,otime,orid,etype,stype,seaz] = dbgetv(db,'sta','arrival.time','origin.time','orid','etype','stype','seaz');
atime = epoch2datenum(atime);
otime = epoch2datenum(otime);

% We also want a list of unique orids.  Have to use DB operations for this
% (rather than unique()) so that sorting (by time) is preserved
db = dbgroup(db,{'origin.time','orid'});
unique_orids = dbgetv(db,'orid');

% Close database link
dbclose(db);

% Display counts
fprintf(2,'%d origins, %d arrivals\n',numel(unique_orids),numel(atime));

%% Arrange info from earthquake catalogue into a more useful format

% We want a column for each station and a row for each earthquake (origin).
% Store arrival times in this format for now, so we can later replace them
% with loaded waveforms.

% Number of stations.  Use user-specified list and order (rather than
% unique values from arrivals table) for greater flexibility.
nstations = length(stations);

% Number of earthquakes
norigins = length(unique_orids);

% Allocate an array to store arrival times
eqtimes = zeros(norigins,nstations);

% Allocate an additional array to store the indices, so we can retrieve
% other data (e.g. event type) later
eqindices = zeros(norigins,nstations);

% Now go through each arrival and place it in the array
for n = 1:length(orid)
	sta_index = find(strcmp(sta{n},stations),1);
	orid_index = find(unique_orids == orid(n),1); % Could do something more efficient here; orid is sorted in ascending order so a find is not really necessary
	eqtimes(orid_index,sta_index) = atime(n);
	eqindices(orid_index,sta_index) = n;
end

%% Load waveforms for all arrivals

% We now have an array full of arrival times at stations.  We want to load
% those waveforms into a cell array of the same shape.

% Allocate the array
eqwf = cell(norigins,nstations);

fprintf(2,'Loading waveforms:\n');
% We'll go column by column
for ns = 1:nstations
	% Get channel names and create scnlobject
	channels = chanlist(ns,2:4);
	channels(strcmp('',channels)) = [];
	scnl = scnlobject(stations{ns},channels,'','');

	% Is this one or three components?
	if length(channels) == 3
		is_tc = true;
	else
		is_tc = false;
	end
	fprintf(2,'%s (%d of %d)...\n',stations{ns},ns,nstations);

	% Create a list of arrival times for this station
	satime_indices = find(eqtimes(:,ns) > 0);
	satimes = eqtimes(satime_indices,ns); % Usually better practice to use logical indexing directly instead of a find(), but we need the indexes themselves for later

	% Calculate actual start and end times
	end_times = satimes + enddelay/86400;
	start_times = satimes - delay/86400;

	% Load waveforms for this station and channels
	wtmp = waveform(ds,scnl,start_times,end_times);

	% Reshape to a Nx1 vector; sometimes waveform() returns unusual
	% dimensions (???)
	wtmp = wtmp(:);

	% If 3-component data, rearrange to put channels together
	if is_tc
		wtmp = reshape(wtmp,size(wtmp,1)/3,3);
	end

	% Store event arrival times in waveform objects, then store the
	% waveforms in our big cell array eqwf.  There doesn't seem to be an
	% obvious way to vectorize these steps.
	% Also store other fields retrieved from database
	skip_count = 0;
	for na = 1:numel(satimes);

		% Check that the waveform is the one we expect; otherwise, skip to
		% the next origin
		wtime = get(wtmp(na-skip_count,1),'START');
		if abs(wtime - start_times(na)) > 0.1/86400
			skip_count = skip_count + 1;
			continue;
		end

		wtmp(na-skip_count,:) = addfield(wtmp(na-skip_count,:),'EVENT_START',satimes(na));
		wtmp(na-skip_count,:) = addfield(wtmp(na-skip_count,:),'ETYPE',cell2mat(etype(eqindices(satime_indices(na),ns))));
 		wtmp(na-skip_count,:) = addfield(wtmp(na-skip_count,:),'ORID',orid(eqindices(satime_indices(na),ns)));
		wtmp(na-skip_count,:) = addfield(wtmp(na-skip_count,:),'OTIME',otime(eqindices(satime_indices(na),ns)));
		wtmp(na-skip_count,:) = addfield(wtmp(na-skip_count,:),'SEAZ',seaz(eqindices(satime_indices(na),ns)));

		% Technically, there can be a different signal type for each picked
		% arrival.  Here we make the assumption that the value of stype for
		% the chosen iphase (typically P) applies to all components.  This
		% will need to be kept in mind if developing a classification
		% routine that makes use of stype, and training catalogues designed
		% accordingly.
		wtmp(na-skip_count,:) = addfield(wtmp(na-skip_count,:),'STYPE',stype(eqindices(satime_indices(na),ns)));

		% Store waveforms in cell array.
		eqwf{satime_indices(na),ns} = wtmp(na-skip_count,:);
	end
end

clear('wtmp'); % No longer needed so free up some space

%% Convert waveforms to threecomp objects where appropriate

% We've loaded all the waveforms.  Now convert the three-component sets
% into threecomp objects

fprintf(2,'Converting to threecomp objects if applicable...\n');

eqtc = cell(size(eqwf));
for no = 1:norigins
	for ns = 1:nstations
		if length(eqwf{no,ns}) == 3 % If it's three-compoenent
			try
				eqtc{no,ns} = threecomp(eqwf{no,ns},eqtimes(no,ns));
			catch me
				warning(['Error converting earthquake (' num2str(no) ',' num2str(ns) ') to threecomp object.  Leaving empty.']);
				eqtc{no,ns} = [];
			end
		else % If it's 1-component (or empty)...
			eqtc{no,ns} = eqwf{no,ns}; % ...store the original value
		end
	end
end

%% Return eqwf if it was requested

if ret_eqwf
	varargout{1} = eqwf;
end

fprintf(2,'Done.\n');

end

author: geoscience-community-codes

repo: GISMO


text:
classdef TestGroupRectanglesMeanShift
    %TestGroupRectanglesMeanShift

    methods (Static)
        function test_cell
            rcts = {[ 0, 1,10,10], ...
                    [ 1, 1,11,11], ...
                    [10,10,20,20], ...
                    [12,12,21,21], ...
                    [30,40,10,20]};
            weights = ones(size(rcts));
            scales = ones(size(rcts));
            [rslts,w] = cv.groupRectangles_meanshift(rcts, weights, scales);
            if ~isempty(rslts)
                validateattributes(rslts, {'cell'}, {'vector'});
                cellfun(@(r) validateattributes(r, {'numeric'}, ...
                    {'vector', 'integer', 'numel',4}), rslts);
                validateattributes(w, {'numeric'}, ...
                    {'vector', 'numel',numel(rslts)});
            end
        end

        function test_error_argnum
            try
                cv.groupRectangles_meanshift();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [obj] = ft_convert_coordsys(obj, target, method, template)

% FT_CONVERT_COORDSYS changes the coordinate system of the input object to
% the specified coordinate system. The coordinate system of the input
% object is determined from the structure field object.coordsys, or needs to
% be determined and specified interactively by the user.
%
% Use as
%   [object] = ft_convert_coordsys(object)
% to only determine the coordinate system, or
%   [object] = ft_convert_coordsys(object, target)
%   [object] = ft_convert_coordsys(object, target, method)
%   [object] = ft_convert_coordsys(object, target, method, template)
% to determine and convert the coordinate system.
%
% With the optional method input argument you can determine whether to use
% SPM for an affine or non-linear transformation. This option is passed on
% to the private functions: align_ctf2acpc and align_neuromag2acpc.
%
% The following input objects are supported
%   anatomical mri, see FT_READ_MRI
%   anatomical or functional atlas, see FT_READ_ATLAS
%   (not yet) electrode definition
%   (not yet) gradiometer array definition
%   (not yet) volume conductor definition
%   (not yet) dipole grid definition
%
% Possible input coordinate systems are 'ctf', 'bti', '4d', 'neuromag' and 'itab'.
% Possible target coordinate systems are 'acpc'.
%
% Note that the conversion will be an automatic one, which means that it
% will be an approximate conversion, not taking into account differences in
% individual anatomies/differences in conventions where to put the
% fiducials.
%
% See also FT_DETERMINE_COORDSYS, ALIGN_CTF2ACPC, ALIGN_NEUROMAG2ACPC, ALIGN_FSAVERAGE2MNI

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

if ~isfield(obj, 'coordsys') || isempty(obj.coordsys)
  % determine the coordinate system of the input object
  obj = ft_determine_coordsys(obj, 'interactive', 'yes');
end

if ~isfield(obj, 'coordsys') || isempty(obj.coordsys)
  % the call to ft_determine_coordsys should have taken care of this, but
  % it is possible that the user aborted the coordinate system
  % determination. See http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2526
  ft_error('the coordinate system of the geometrical object is not specified');
end

if any(strcmp(target, {'spm', 'mni', 'tal'}))
  % see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=3304
  ft_warning('Not applying any scaling, using ''acpc'' instead of ''%s''. See http://bit.ly/2sw7eC4', target);
  target = 'acpc';
end

% set default behavior to use an approximate alignment, followed by a call
% to spm_normalise for a better quality alignment
if nargin<3
  method = 2;
end

if isdeployed && (method==1 || method==2)
  needtemplate = true;
else
  needtemplate = false;
end

if needtemplate && nargin<4
  ft_error('you need to specify a template filename for the coregistration');
end

hastemplate = nargin>3;

if nargin>1 && ~strcmpi(target, obj.coordsys)
  % convert to the desired coordinate system
  switch target
    case 'acpc'
      switch obj.coordsys
        case {'ctf' 'bti' '4d'}
          fprintf('converting the coordinate system from %s to %s\n', obj.coordsys, target);
          if hastemplate
            obj = align_ctf2acpc(obj, method, template);
          else
            obj = align_ctf2acpc(obj, method);
          end
        case {'itab' 'neuromag'}
          fprintf('converting the coordinate system from %s to %s\n', obj.coordsys, target);
          if hastemplate
            obj = align_neuromag2acpc(obj, method, template);
          else
            obj = align_neuromag2acpc(obj, method);
          end
        case {'mni', 'fsaverage'}
          fprintf('not converting the coordinate system from %s to %s, these are similar enough\n', obj.coordsys, target);
          % these are close enough to ACPC, so nothing needs to be done
        otherwise
          ft_warning('conversion from %s to %s is not supported', obj.coordsys, target);
      end % switch obj.coordsys
    otherwise
      ft_error('conversion from %s to %s is not supported', obj.coordsys, target);
  end % switch target
end

author: spm


text:
function y = lt( x, y )

if isa(x,'cvxin')||~isa(y,'cvxin')||y.active,
    error( 'CVX error: improper use of the <in> pseudo-operator.' );
end
y.active = true;
y.value  = x;

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
load('shot.mat');                                                   % loading image
colormap('gray');                                                   % choosing colourmap
image(shot);                                                        % showing image
[py, px] = findcentroidlaserspotusingfourier(shot)                  % calculating centroid
if (isnan(py) == 0) && (isnan(px) == 0)                             % if calculated values are not NaN
    line([px-10; px+10], [py, py], 'color', 'g', 'LineWidth', 1);   % showing haircross
    line([px, px], [py-10, py+10], 'color', 'g', 'LineWidth', 1);
end

author: Sable

repo: mcbench-benchmarks


text:
function m = moreland(n)
%MORELAND   Divergent red-white-blue color map
%
%   Usage: m = moreland([n])
%
%   Input parameters:
%       n - optional length of colormap (default uses the figure default length)
%
%   Output parameters:
%       m - colormap [n 3]
%
%   MORELAND(N) returns an N-by-3 matrix containing a divergent colormap.
%   Without a given N the same length as the current figure's colormap is used.
%   For details on the colormap have a look at:
%   http://www.kennethmoreland.com/color-advice/
%
%   To change the colormap current figure run: colormap(moreland)
%
%   Original written by Lewis Marshall.
%
%   See also: generate_colormap, yellowred, plot_sound_field

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking input parameters =======================================
nargmin = 0;
nargmax = 1;
narginchk(nargmin,nargmax);
if nargin < nargmax
    n = size(get(gcf,'colormap'),1);
end

%% ===== Computation =====================================================
table = [ 59,76,192
60,78,194
61,80,195
62,81,197
63,83,198
64,85,200
66,87,201
67,88,203
68,90,204
69,92,206
70,93,207
71,95,209
73,97,210
74,99,211
75,100,213
76,102,214
77,104,215
79,105,217
80,107,218
81,109,219
82,110,221
84,112,222
85,114,223
86,115,224
87,117,225
89,119,226
90,120,228
91,122,229
93,123,230
94,125,231
95,127,232
96,128,233
98,130,234
99,131,235
100,133,236
102,135,237
103,136,238
104,138,239
106,139,239
107,141,240
108,142,241
110,144,242
111,145,243
112,147,243
114,148,244
115,150,245
116,151,246
118,153,246
119,154,247
120,156,247
122,157,248
123,158,249
124,160,249
126,161,250
127,163,250
129,164,251
130,165,251
131,167,252
133,168,252
134,169,252
135,171,253
137,172,253
138,173,253
140,174,254
141,176,254
142,177,254
144,178,254
145,179,254
147,181,255
148,182,255
149,183,255
151,184,255
152,185,255
153,186,255
155,187,255
156,188,255
158,190,255
159,191,255
160,192,255
162,193,255
163,194,255
164,195,254
166,196,254
167,197,254
168,198,254
170,199,253
171,199,253
172,200,253
174,201,253
175,202,252
176,203,252
178,204,251
179,205,251
180,205,251
182,206,250
183,207,250
184,208,249
185,208,248
187,209,248
188,210,247
189,210,247
190,211,246
192,212,245
193,212,245
194,213,244
195,213,243
197,214,243
198,214,242
199,215,241
200,215,240
201,216,239
203,216,238
204,217,238
205,217,237
206,217,236
207,218,235
208,218,234
209,219,233
210,219,232
211,219,231
213,219,230
214,220,229
215,220,228
216,220,227
217,220,225
218,220,224
219,220,223
220,221,222
221,221,221
222,220,219
223,220,218
224,219,216
225,219,215
226,218,214
227,218,212
228,217,211
229,216,209
230,216,208
231,215,206
232,215,205
232,214,203
233,213,202
234,212,200
235,212,199
236,211,197
236,210,196
237,209,194
238,209,193
238,208,191
239,207,190
240,206,188
240,205,187
241,204,185
241,203,184
242,202,182
242,201,181
243,200,179
243,199,178
244,198,176
244,197,174
245,196,173
245,195,171
245,194,170
245,193,168
246,192,167
246,191,165
246,190,163
246,188,162
247,187,160
247,186,159
247,185,157
247,184,156
247,182,154
247,181,152
247,180,151
247,178,149
247,177,148
247,176,146
247,174,145
247,173,143
247,172,141
247,170,140
247,169,138
247,167,137
247,166,135
246,164,134
246,163,132
246,161,131
246,160,129
245,158,127
245,157,126
245,155,124
244,154,123
244,152,121
244,151,120
243,149,118
243,147,117
242,146,115
242,144,114
241,142,112
241,141,111
240,139,109
240,137,108
239,136,106
238,134,105
238,132,103
237,130,102
236,129,100
236,127,99
235,125,97
234,123,96
233,121,95
233,120,93
232,118,92
231,116,90
230,114,89
229,112,88
228,110,86
227,108,85
227,106,83
226,104,82
225,102,81
224,100,79
223,98,78
222,96,77
221,94,75
220,92,74
218,90,73
217,88,71
216,86,70
215,84,69
214,82,67
213,80,66
212,78,65
210,75,64
209,73,62
208,71,61
207,69,60
205,66,59
204,64,57
203,62,56
202,59,55
200,57,54
199,54,53
198,51,52
196,49,50
195,46,49
193,43,48
192,40,47
190,37,46
189,34,45
188,30,44
186,26,43
185,22,41
183,17,40
181,11,39
180,4,38];
m = generate_colormap(table,n);

author: sfstoolbox

repo: sfs-matlab


text:
function r = isfinite(a)
%ISFINITE     Array of 1's for finite components
%
%   r = isfinite(a)
%

% written  05/21/09     S.M. Rump
%

  r = reshape(all(isfinite(a.t)),a.size);

repo: OpenMAS


text:
function [lmval,indd]=lmax_pw(xx, dx)
%   Find  piece-wise  local maxima in vector XX,where
%	LMVAL is the output vector with maxima values, INDD  is the
%	corresponding indexes, DX is length of piece where maxima is searched,
%   IMPORTANT:      FIRST and LAST point in vector are excluded
%   IMPORTANT:      XX must be single column vector
%   IMPORTANT:      Length of DX must be very carefully selected
%	For example compare dx=30; and dx=1000;
%
%   dx=150; xx=[0:0.01:35]'; y=sin(xx .* cos(xx /4.5)) + cos(xx);
%    plot(xx,y); grid; hold on;
%   %   Excluding first and last points
%   [b,a]=lmax_pw(y,dx); plot(xx(a),y(a),'r+')
%   % Way to include first and last points can be as:
%   y(1)=1.5; yy=[0; y; -1;];   % padd with smaller values
%   [b,a]=lmax_pw(yy,dx); a=a-1; plot(xx(a),y(a),'go')
%
%	see also LMIN, LMAX,  LMIN_PW, MATCH

% 	Sergei Koptenko, Applied Acoustic Technologies, Toronto, Canada
%   sergei.koptenko@sympatico.ca,  March/11/2003

if nargin <2,
	disp('Not enough arguments'); return
end

len_x = length(xx);
xx = [xx; xx(len_x); xx(len_x)];
nn=floor(len_x/dx);
ncount=1; lmval=[]; indd=[];
	for ii=1:nn,
        [lm,ind] = max(xx(ncount: ii*dx+2)) ;
        ind=ind+(ii-1)*dx;
                 if (ind ~=ncount) & (ind~=ii*dx+2),
                    lmval=[lmval, lm]; indd=[indd, ind];
                end
        ncount=ncount +dx;
	end
[lm,ind] = max(xx(ii*dx:len_x));
        if (ind ~=len_x) & (ind~=ii*dx),
            lmval=[lmval, lm]; indd=[indd, (ind+ii*dx-1)];
        end

       if indd(end)==len_x,
           indd=  indd(1:end-1);
           lmval=lmval(1:end-1);
       end
return


text:
function [m3] = in32m3(in3)
% Convert volume from cubic inches to cubic meters.
% Chad Greene 2012
m3 = in3*0.000016387064;


text:
clear
if(isunix)
    features_exe = '"../../../build/bin/FeatureExtraction"';
else
    features_exe = '"../../../x64/Release/FeatureExtraction.exe"';
end
ck_loc = 'E:\Datasets\ck+\cohn-kanade-images\';

out_loc = 'E:\datasets\face_datasets_processed/ck+';

% Go two levels deep
ck_dirs = dir(ck_loc);
ck_dirs = ck_dirs(3:end);

parfor f1=1:numel(ck_dirs)

    ck_dirs_level_2 = dir([ck_loc, ck_dirs(f1).name]);
    ck_dirs_level_2 = ck_dirs_level_2(3:end);

    for f2=1:numel(ck_dirs_level_2)

        if(~isdir([ck_loc, ck_dirs(f1).name, '/', ck_dirs_level_2(f2).name]))
           continue;
        end

        curr_vid = [ck_loc, ck_dirs(f1).name, '/', ck_dirs_level_2(f2).name];

        name = [ck_dirs(f1).name, '_',  ck_dirs_level_2(f2).name];

        command = sprintf('%s -fdir "%s" -out_dir "%s" -of %s -hogalign -pdmparams', features_exe, curr_vid, out_loc, name);

        dos(command);

    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function apply_derivative_boost(imgtype, varargin)
% Allows for recalculation of amplitude images from the fitted responses
% for each trial type. Necessary for calculating group statistics when
% using multiple basis functions
%
% :Usage:
% ::
%
%     function apply_derivative_boost(varargin)
%
% :NOTES:
%    DB estimation only works for 2 specific basis sets:
%    timeonly, with 2 parameters (canonical + time derivative)
%    timedispersion,' with 3 parameters (canonical hrf + temporal and spatial dispersion)
%
% This function loads the SPM.mat file in the current directory and uses
% the basis set specified in the loaded SPM structure.
%
% : Required Inputs:
% 	indicate whether images are'.img' or '.nii'
%
% :Optional Inputs:
%
%   **'amplitudes'**
%        will only create amping images (combination of betas
%        across basis functions)
%
%   **'contrasts'**
%        assumes that amping images are already created; will only
%        create contrast images
%
%   **'all'**
%        will run both the amplitudes and contrasts sections
%
%        In addition, 'amplitudes' now has two separate parts:
%          - The first uses Vince Calhoun's derivative boost (Calhoun, 2004) to
%             estimate amplitudes.  NOTE: *We have not worked out the scaling yet, so
%             I'm not sure this is working right*
%             To turn this OFF, enter 'nodb' as an optional argument
%
%          - The second way uses our HTW code to estimate height, time to peak,
%            width, and area under the curve (see Lindquist & Wager 2007 for
%            simulations using a version of this code).
%            It requires SCANlab specific functions, in SCN_Core_Support
%            (unlike the deriv. boost).
%            To turn this OFF, enter 'nohtw' as an optional argument
%
%   **'startend'**
%        followed by starting and ending values in seconds for amplitude
%        estimation window (for HTW estimation only).
%        If you do not enter this, it will show you a plot and ask you to pick
%        these values.
%        If you enter them here as inputs, you can skip the interactive step and
%        loop over subjects.
%
%   **'condition_numbers'**
%        followed by which index numbers in your list should
%        be used to calculate h, t, w from.  You should use this if you
%        are entering regressors of no interest, besides the intercepts.
%
% :Important for Contrasts:
%    disp('Using contrasts entered as F-contrasts. Assuming the first contrast vector in each F-contrast '
%
%    disp('is a contrast of interest across the CANONICAL basis function regressors.')
%
% :Examples:
% ::
%
%    % RUN THIS IN COMMAND WINDOW TO BATCH
%    subj = dir('06*')
%    for i = 1:length(subj), cd(subj(i).name), apply_derivative_boost, cd('..'); end
%
%    % ANOTHER BATCH EXAMPLE:
%    d = dir('remi*'); d = d(cat(2, d.isdir)); [mydirs{1:length(d)}] = deal(d.name)
%    for i = 1:length(mydirs), cd(mydirs{i}), apply_derivative_boost('all', 'nodb', 'startend', [4 15]), cd('..'); end
%
%    %An example for an event-related design, specifying condition numbers to get HTW from:
%    apply_derivative_boost('all', 'nodb', 'contrasts', 'condition_numbers', 1:14, 'startend', [4 10]);
%
%    % CALCULATE CONTRASTS ONLY ON ALREADY-ESTIMATED HTW IMAGES
%    apply_derivative_boost('contrasts','condition_numbers',1:14);
%

    imgtype = imgtype;
    spmname = fullfile(pwd, 'SPM.mat');     % SETUP INPUTS
    if ~exist(spmname, 'file'), error('You must be in an SPM 1st-level results directory with SPM5 SPM.mat file.'); end

    load(spmname);

    nbf = size(SPM.xBF.bf, 2);

    switch nbf
        case 2
            derivative_case = 'timeonly';
        case 3
            derivative_case = 'timedispersion';

        otherwise
            warning('Deriv. Boost only works for SPM canonical hrf with time or time + dispersion derivatives.  This SPM.mat doesn''t match those specs.');
    end

    % Not used; only for image eval function
    % boost = @(b) sign(b(1)) .* sqrt(sum(b .^ 2));

    docontrasts = 0;
    doamps = 0;
    nodb = 0;
    nohtw = 0;
    condition_numbers = [];
    do_downsample = [];          % downsample; default = 1 sec if units are seconds (1/dt)

    for i = 1:length(varargin)
        if ischar(varargin{i})
            switch varargin{i}
                % reserved keywords
                case 'all', docontrasts = 1; doamps = 1;

                case 'amplitudes', doamps = 1;

                case 'contrasts', docontrasts = 1;

                case 'nodb', nodb = 1; % skip DB estimation

                case 'nohtw', nohtw = 1; % skip HTW estimation

                case 'startend', startend = varargin{i + 1};  % starting and ending values in seconds for amplitude estimation window (for HTW estimation only).

                case 'condition_numbers',  condition_numbers = varargin{i + 1};

                case 'nodownsample', do_downsample = 0;

                case 'downsample', do_downsample = varargin{i + 1};

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

    if isempty(do_downsample)
        % default downsampling
        do_downsample = round(1 ./ SPM.xBF.dt);
    end

    if ~(docontrasts || doamps)
        disp('Nothing to do! Enter ''contrasts'' ''amplitudes'' or ''all'' as input argument.');
        return
    end

    if doamps
        % ---------------------------------------------
        % ---------------------------------------------

        % ESTIMATE AMPLITUDES

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

        % ---------------------------------------------
        % FILE NAMES
        % ---------------------------------------------

        imgs = dir(sprintf(['beta*', imgtype])); imgs = char(imgs.name);
        n = size(imgs, 1);

        fprintf('Found %3.0f images', n); fprintf('\n');

        %load(spmname);
        nsess = length(SPM.Sess);
        fprintf('I think there are %3.0f sessions (runs)', nsess); fprintf('\n');

        if ~isempty(condition_numbers)
            wh_intercept = true(1, size(imgs, 1));  % exclude these
            wh_intercept(condition_numbers) = 0;
            wh_intercept = find(wh_intercept);

            fprintf('\nIncluding only these images: ');
            fprintf('%3.0f ', condition_numbers);
            fprintf('\n');

        else
            wh_intercept = SPM.xX.iB;
            fprintf('\nI think these images are intercepts, and am not using them: ');
            fprintf('%3.0f ', wh_intercept);
            fprintf('\n');
        end

        imgs(wh_intercept, :) = [];

        %not used; only when using image_eval_function
        %mask_img = './mask.img';

        n = size(imgs, 1);

        if n / nbf ~= round(n / nbf), error('Error!  Wrong number of images for the specified number of basis functions.'); end

        if nodb
            % skip Derivative Boost estimation and go straight to HTW
        else

            % ---------------------------------------------
            % CALCULATE
            % ---------------------------------------------
            cond_indx = 1;

            for i = 1: nbf : (n - nbf + 1)

                imgs_cond = imgs(i : i+nbf - 1, :);

                disp('Working on :')
                disp(imgs_cond)

                out_name = sprintf(['db_amplitude_%03d',imgtype ], cond_indx);

                % This code uses SCN lab tools to create images
                % ---------------------------------------------
                % %     y = image_eval_function(imgs_cond, boost, 'mask', mask_img, ...
                % %         'outimagelabels', {out_name});
                % ---------------------------------------------

                % This code uses SPM instead
                % ---------------------------------------------
                switch derivative_case
                    case 'timeonly'
                        spm_imcalc(imgs_cond, out_name, 'sign(i1) .* sqrt(i1.^2 + i2.^2)');
                    case 'timedispersion'
                        spm_imcalc(imgs_cond, out_name, 'sign(i1) .* sqrt(i1.^2 + i2.^2 + i3.^2)');
                    otherwise
                        error('Basis set is incompatible with DB estimation!');
                end

                fprintf('Created %s\n', out_name);
                % ---------------------------------------------

                cond_indx = cond_indx + 1;
            end

            % Get and save names
            db_amp_names = [];
            for i = 1:nsess
                sessnames = char(SPM.Sess(i).Fc.name);
                sessnames = [repmat(sprintf('Sess%02d_', i), size(sessnames, 1), 1) sessnames];
                db_amp_names = strvcat(db_amp_names, sessnames);
            end

            save db_amplitude_names db_amp_names
            disp(db_amp_names);
            disp(' ')
            disp('Saved DB amplitude condition names for each image in db_amplitude_names.mat');

            fprintf('\n*-----------------------------*\nApplied DB successfully\n*-----------------------------*\n')

        end

        if nohtw
            % skip this

        else

            % ---------------------------------------------
            % Estimated amplitude from fit: HTW
            % ---------------------------------------------
            % This code uses SCN lab tools to create images

            disp(' ')
            disp('Next: Estimating amplitude, time to peak, width, and area-under-curve images from fitted response using SCN lab code.')
            disp(' ')

            % downsample bf, if requested
            if do_downsample
                mytimeres = SPM.xBF.dt * do_downsample;
                SPM.xBF.bf = downsample(SPM.xBF.bf, do_downsample);

            else
                mytimeres = SPM.xBF.dt;
            end

            if exist('startend', 'var')
                % just check, and use input values
                if length(startend) ~= 2, error('Startend input must have two values, a starting and ending value in seconds for the amp. estimate window'), end
            else
                % Set range in sec
                htw_from_fit(SPM.xBF.bf, ones(size(SPM.xBF.bf, 2), 1), mytimeres, 'plot', 'verbose');

                disp(' ')
                disp('Enter the range in seconds within which to estimate peak amplitude.')
                disp('Example: type [2.5 9] and press return for a typical event-related setup.');
                disp('More sustained responses, like pain responses, may require a longer window.');
                disp('This estimates the amplitude of the IMPULSE RESPONSE, before convolution with the stimulus function')
                disp('so if you have an epoch design, a typical window of [2.5 9] sec is still appropriate.');
                disp('Also note: AUC images are calculated as the area under the curve within the window you specify.')
                disp(' ')

                startend = input('Enter your choice in [ ] and press return: ');
            end

            % Test your choice by showing you a plot
            htwfunction = @(b) htw_from_fit(SPM.xBF.bf, b, mytimeres, 'startval', startend(1), 'endval', startend(2), 'plot', 'verbose');
            htwfunction(ones(size(SPM.xBF.bf, 2), 1))
            drawnow

            % Create without plot option for loop through brain.
            htwfunction = @(b) htw_from_fit(SPM.xBF.bf, b, mytimeres, 'startval', startend(1), 'endval', startend(2));

            disp('Check the screen for a plot of your choice of window.')
            disp(' ')

            % ---------------------------------------------
            % CALCULATE
            % ---------------------------------------------
            cond_indx = 1;

            for i = 1: nbf : (n - nbf + 1)

                imgs_cond = imgs(i : i+nbf - 1, :);

                disp('Working on :')
                disp(imgs_cond)

                clear out_name
                out_name{1} = sprintf(['htw_amplitude_%03d',imgtype], cond_indx);
                out_name{2} = sprintf(['htw_time_to_peak_%03d',imgtype], cond_indx);
                out_name{3} = sprintf(['htw_width_%03d',imgtype], cond_indx);
                out_name{4} = sprintf(['htw_area_under_curve_%03d',imgtype], cond_indx);

                % ---------------------------------------------
                [h, t, w, auc] = image_eval_function(imgs_cond, htwfunction, 'mask', ['mask', imgtype], ...
                    'outimagelabels', out_name);

                h;t;w;auc;  % we need the outputs above to tell it to write 4 images.
                % ---------------------------------------------

                cond_indx = cond_indx + 1;
            end

            % Get and save names
            htw_amp_names = [];
            for i = 1:nsess
                sessnames = char(SPM.Sess(i).Fc.name);
                sessnames = [repmat(sprintf('Sess%02d_', i), size(sessnames, 1), 1) sessnames];
                htw_amp_names = strvcat(htw_amp_names, sessnames);
            end

            if ~exist(fullfile(pwd, 'db_amplitude_names.mat'), 'file')
                save db_amplitude_names htw_amp_names
            else
                save db_amplitude_names -append htw_amp_names
            end
            disp(htw_amp_names);
            disp(' ')
            disp('Saved HTW amplitude condition names for each image in db_amplitude_names.mat');

            fprintf('\n*-----------------------------*\nApplied HTW estimation successfully\n*-----------------------------*\n')

        end

    end  % amplitudes

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

    % CREATE CONTRAST FOR THIS SUBJECT

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

    if docontrasts

        % Load contrast vectors

        % -----------------------------------------
        %load(spmname);
        % nsess = length(SPM.Sess);

        % Define which indices to exclude from contrasts
        if ~isempty(condition_numbers)
            nconvals = length(SPM.xCon(1).c(:, 1));  % test contrast
            wh_intercept = true(1, nconvals);  % exclude these
            wh_intercept(condition_numbers) = 0;
            wh_intercept = find(wh_intercept);

            fprintf('\nIncluding only these images: ');
            fprintf('%3.0f ', condition_numbers);
            fprintf('\n');

        else
            wh_intercept = SPM.xX.iB;
            fprintf('\nI think these images are intercepts, and am not using them: ');
            fprintf('%3.0f ', wh_intercept);
            fprintf('\n');
        end

        if ~isfield(SPM, 'xCon')
            error('Enter F-contrasts first, with the first contrast vector in each F-contrast the contrast across the CANONICAL basis function.');
        else
            disp('Using contrasts entered as F-contrasts. Assuming the first contrast vector in each F-contrast ')
            disp('is a contrast of interest across the CANONICAL basis function regressors.')
        end

        wh_F = strmatch('F', char(SPM.xCon.STAT), 'exact');

        if isempty(wh_F), error('No F-contrasts entered yet.'); end

        % All sets of contrast images
        % ------------------------------------------
        ampimgs_name = sprintf(['db_amplitude_*',imgtype]);

        ampimgs = dir(ampimgs_name); ampimgs = char(ampimgs.name);

        if ~isempty(ampimgs)
            contrast_image_names_dbamp = calc_contrasts(ampimgs_name, ampimgs, wh_F, SPM, wh_intercept, nbf);
        else
            disp(['Checked for but did not find: ']);
        end

        % HTW amplitude
        % ------------------------------------------
        ampimgs_name = sprintf(['htw_amplitude_*',imgtype]);

        ampimgs = dir(ampimgs_name); ampimgs = char(ampimgs.name);

        if ~isempty(ampimgs)
            contrast_image_names_htwamp = calc_contrasts(ampimgs_name, ampimgs, wh_F, SPM, wh_intercept, nbf);
        else
            disp(['Checked for but did not find: ' ampimgs_name]);
        end

        % HTW time
        % ------------------------------------------
        ampimgs_name = sprintf(['htw_time_to_peak_*',imgtype]);

        ampimgs = dir(ampimgs_name); ampimgs = char(ampimgs.name);

        if ~isempty(ampimgs)
            contrast_image_names_htwtime = calc_contrasts(ampimgs_name, ampimgs, wh_F, SPM, wh_intercept, nbf);
        else
            disp(['Checked for but did not find: ' ampimgs_name]);
        end

        % HTW width
        % ------------------------------------------
        ampimgs_name = sprintf(['htw_width_*',imgtype]);

        ampimgs = dir(ampimgs_name); ampimgs = char(ampimgs.name);

        if ~isempty(ampimgs)
            contrast_image_names_htwwid = calc_contrasts(ampimgs_name, ampimgs, wh_F, SPM, wh_intercept, nbf);
        else
            disp(['Checked for but did not find: ' ampimgs_name]);
        end

        % HTW area
        % ------------------------------------------
        ampimgs_name = sprintf(['htw_area_under_curve_*',imgtype]);

        ampimgs = dir(ampimgs_name); ampimgs = char(ampimgs.name);

        if ~isempty(ampimgs)
            contrast_image_names_htwarea = calc_contrasts(ampimgs_name, ampimgs, wh_F, SPM, wh_intercept, nbf);
        else
            disp(['Checked for but did not find: ' ampimgs_name]);
        end

        check_it = whos('contrast_image_names*');

        if isempty(check_it)
            disp('No valid images found to create contrasts on');

        else
            save('db_amplitude_names', '-append', 'contrast_image_names*');
            disp('Saved lists of contrast image names in db_amplitude_names.mat');
        end

    end

    %% INLINE

end  % main function

function contrast_image_names = calc_contrasts(ampimgs_name, ampimgs, wh_F, SPM, wh_intercept, nbf)

    n = size(ampimgs, 1);

    fprintf('Found %3.0f images:\n', n);
    disp(ampimgs)

    disp('Reading image data.')
    V = spm_vol(ampimgs);
    vols = spm_read_vols(V);

    % spm_check_registration(ampimgs);
    % colormap jet

    contrast_image_names = [];

    for i = 1:length(wh_F)

        name = SPM.xCon(wh_F(i)).name;
        disp(['Calculating contrast on: ' name])
        original_name = name;

        name = deblank(name);
        wh_bad = (name == ' ' | name == ',' | name == '.' | name == '^' | name == '~' | name == '''' | name == ':' | name == '*' | name == '%' | name == ';' | name == '@' | name == '&');
        name(wh_bad) = [];

        name = ['con_', ampimgs_name(1:8), '_', name, imgtype];

        c = SPM.xCon(wh_F(i)).c(:, 1);
        c(wh_intercept) = [];
        c = c(1 : nbf : end);

        if length(c) ~= n, error('Contrast is wrong length for some reason!  Coding error in this function?  Or wrong number of db_amplitude images.'); end

        fprintf('Contrast values: ')
        fprintf('%01d ', c)
        fprintf('\n')

        % calculate and save

        contrast_image_calc(name, c, vols, V, original_name)

        disp(['Written: ' name]);

        contrast_image_names = strvcat(contrast_image_names, name);

        disp(' ');

    end

    fprintf('\n*-----------------------------*\nContrasts Done successfully!\n*-----------------------------*\n')
    spm_check_registration(contrast_image_names);

end

function contrast_image_calc(Q, myc, vols, V, original_name)

    % FROM:
    % function contrast_image(Q, myc)
    %
    % Tor Wager
    %
    % Creates a contrast image called Q (do not include path)
    % Given a list of img files P (spm format, with path)
    % and a contrast vector myc
    % In the directory of 1st image in P.

    if ~(length(myc) == size(vols,4))
        error('Contrast vector length is not equal to number of image files.')
    end

    myc2 = zeros(size(vols));

    for i = 1:length(myc)
        myc2(:,:,:,i) = myc(i);
    end

    cvol = vols .* myc2;
    cvol = sum(cvol,4);

    % -------------------------
    % write
    % -------------------------
    dd = fileparts(V(1).fname);
    Q = fullfile(dd, Q);
    Vo = V(1);
    Vo.fname = Q;
    Vo.descrip = ['Contrast ' original_name];

    spm_write_vol(Vo,cvol);

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function printSampleStats(sampledModel, commonModel, sampleNames, fileName)
% Prints out sample statistics for multiple samples
%
% USAGE:
%
%    printSampleStats(samples, commonModel, sampleNames, fileName)
%
% INPUTS:
%    sampledModel:   Samples to plot
%    commonModel:    COBRA model structure
%    sampleNames:    Names of the models
%
% OPTIONAL INPUT:
%    fileName:       Name of tab delimited CSV file to generate
%                    (Default = print to command window)
%
% .. Author: - Markus Herrgard

if nargin > 3
    fid = fopen(fileName, 'w');
else
    fid = 1;
end

sampleStats = calcSampleStats(sampledModel);

fprintf(fid, 'Rxn\t');
if (isfield(commonModel, 'subSystems'))
    fprintf(fid, 'Subsystem\t');
end

for i = 1:length(sampleNames)
    fprintf(fid, '%s-mode\t', sampleNames{i});
end
for i = 1:length(sampleNames)
    fprintf(fid, '%s-mean\t', sampleNames{i});
end
for i = 1:length(sampleNames)
    fprintf(fid, '%s-median\t', sampleNames{i});
end
for i = 1:length(sampleNames)
    fprintf(fid, '%s-std\t', sampleNames{i});
end
fprintf(fid, '\n');

for i = 1:length(commonModel.rxns)
    fprintf(fid, '%s\t', commonModel.rxns{i});
    if (isfield(commonModel, 'subSystems'))
        fprintf(fid, '%s\t', strjoin(commonModel.subSystems{i},';'));
    end
    %for j = 1:length(samples)
    fprintf(fid, '%8.6f\t', sampleStats.mode(i));
    %end
    %for j = 1:length(samples)
    fprintf(fid, '%8.6f\t', sampleStats.mean(i, :));
    %end
    %for j = 1:length(samples)
    fprintf(fid, '%8.6f\t', sampleStats.median(i, :));
    %end
    %for j = 1:length(samples)
    fprintf(fid, '%8.6f\t', sampleStats.std(i, :));
    %end
    fprintf(fid, '\n');
end
if (fid > 1)
    fclose(fid);
end

author: opencobra

repo: cobratoolbox


text:
function [ft3] = m32ft3(m3)
% Convert volume from cubic meters to cubic feet.
% Chad Greene 2012
ft3 = m3*35.314666721;


text:
function test_failed = test_libltfat_fftcircshift(varargin)
test_failed = 0;

fprintf(' ===============  %s ================ \n',upper(mfilename));

definput.flags.complexity={'double','single'};
[flags]=ltfatarghelper({},definput,varargin);
dataPtr = [flags.complexity, 'Ptr'];

Larr = [1,9,11,110,111];
shiftarr = [-10, 100, 0, -1, 1, 3, 52];

for L = Larr
    for shift = shiftarr

        z = cast((1:L)' + 1i*(L:-1:1)',flags.complexity);
        zi = complex2interleaved(fft(z));
        zout = randn(size(zi),flags.complexity);

        ziPtr = libpointer(dataPtr,zi);
        zoutPtr = libpointer(dataPtr,zout);

        trueres = circshift(z,shift);

        status = calllib('libltfat',makelibraryname('fftcircshift',flags.complexity,1),...
            ziPtr,L,shift,zoutPtr);

        res = norm(trueres - ifft(interleaved2complex(zoutPtr.Value)));

        [test_failed,fail]=ltfatdiditfail(res+status,test_failed);
        fprintf(['FFTCIRCSHIFT OP L:%3i, shift:%3i %s %s %s\n'],L,shift,flags.complexity,ltfatstatusstring(status),fail);

        status = calllib('libltfat',makelibraryname('fftcircshift',flags.complexity,1),...
            ziPtr,L,shift,ziPtr);

        res = norm(trueres - ifft(interleaved2complex(ziPtr.Value)));

        [test_failed,fail]=ltfatdiditfail(res+status,test_failed);
        fprintf(['FFTCIRCSHIFT IP L:%3i, shift:%3i %s %s %s\n'],L,shift,flags.complexity,ltfatstatusstring(status),fail);

        z = cast((1:L)',flags.complexity);
        zi = complex2interleaved(fftreal(z));
        zout = randn(size(zi),flags.complexity);

        ziPtr = libpointer(dataPtr,zi);
        zoutPtr = libpointer(dataPtr,zout);

        trueres = circshift(z,shift);

        status = calllib('libltfat',makelibraryname('fftrealcircshift',flags.complexity,1),...
            ziPtr,L,shift,zoutPtr);

        res = norm(trueres - real(ifftreal(interleaved2complex(zoutPtr.Value),L)));

        [test_failed,fail]=ltfatdiditfail(res+status,test_failed);
        fprintf(['FFTREALCIRCSHIFT OP L:%3i, shift:%3i %s %s %s\n'],L,shift,flags.complexity,ltfatstatusstring(status),fail);

        status = calllib('libltfat',makelibraryname('fftrealcircshift',flags.complexity,1),...
            ziPtr,L,shift,ziPtr);

        res = norm(trueres - real(ifftreal(interleaved2complex(ziPtr.Value),L)));

        [test_failed,fail]=ltfatdiditfail(res+status,test_failed);
        fprintf(['FFTREALCIRCSHIFT IP L:%3i, shift:%3i %s %s %s\n'],L,shift,flags.complexity,ltfatstatusstring(status),fail);
    end
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [StopFlag, Status] = StopCriterion(grad, nhxk, Niter, Nmap, Ngmap, MaxNumIter, MaxNumMapEval, MaxNumGmapEval, T, TimeLimit, epsilon, nhx0, ngradx0, Stopping_Crit)
% `StopCriterion` is a function checking that one of the stopping criteria
% holds to terminate LLM and GLM. It perepare the status determining why
% the algorithm is stopped.
%
% USAGE:
%
%    [StopFlag, Status] = StopCriterion(grad, nhxk, Niter, Nmap, Ngmap, MaxNumIter, MaxNumMapEval, MaxNumGmapEval, T, TimeLimit, epsilon, nhx0, ngradx0, Stopping_Crit)
%
% INPUTS:
%    grad:              gradient of the merit funcrion
%    nhxk:              the norm 2 of `h(xk)`
%    MaxNumIter:        maximum number of iterations
%    MaxNumMapEval:     maximum number of function evaluations
%    MaxNumGmapEval:    maximum number of subgradient evaluations
%    TimeLimit:         maximum running time
%    epsilon:           accuracy parameter
%    Stopping_Crit:     stopping criterion
%
%                         1. stop if :math:`||grad|| \leq \epsilon`
%                         2. stop if :math:`||nhxk|| \leq \epsilon`
%                         3. stop if `MaxNumIter` is reached
%                         4. stop if `MaxNumMapEval` is reached
%                         5. stop if `MaxNumGmapEval` is reached
%                         6. stop if `TimeLimit` is reached
%                         7. stop if :math:`||grad|| \leq \textrm{max}(\epsilon, \epsilon^2 * ngradx0)`
%                         8. stop if :math:`||nhxk|| \leq \textrm{max}(\epsilon, \epsilon^2 * nhx0)`
%                         9. stop if (default) :math:`||hxk|| \leq \epsilon` or `MaxNumIter` is reached
%
% OUTPUTS:
%    StopFlag:          1: if one of the stopping criteria holds, 0: if none of the stopping criteria holds
%    Status:            the reason of the scheme termination

switch Stopping_Crit
  case 1
    if norm(grad) <= epsilon
      StopFlag = 1;
      Status   = 'Local (global) solution of merit function is found.';
    else
      StopFlag = 0;
      Status   = [];
    end
  case 2
    if nhxk <= epsilon
      StopFlag = 1;
      Status   = 'A solution of nonlinear system is found.';
    else
      StopFlag = 0;
      Status   = [];
    end
  case 3
    if Niter >= MaxNumIter
      StopFlag = 1;
      Status   = 'Maximum number of iterations is reached.';
    else
      StopFlag = 0;
      Status   = [];
    end
  case 4
    if Nmap >= MaxNumMapEval
      StopFlag = 1;
      Status   = 'Maximum number of mapping evaluations is reached.';
    else
      StopFlag = 0;
      Status   = [];
    end
  case 5
    if Ngmap >= MaxNumGmapEval
      StopFlag = 1;
      Status   = 'Maximum number of gradient evaluations is reached.';
    else
      StopFlag = 0;
      Status   = [];
    end
  case 6
    if T >= TimeLimit
      StopFlag = 1;
      Status   = 'Time limit is reached.';
    else
      StopFlag = 0;
      Status   = [];
    end

  case 7
    if norm(grad) <= max(epsilon,epsilon^2*ngradx0)
      StopFlag = 1;
      Status   = 'A possible stationary point is found.';
    else
      StopFlag = 0;
      Status   = [];
    end

  case 8
    if nhxk <= max(epsilon,epsilon^2*nhx0)
      StopFlag = 1;
      Status   = 'A possible solution is found.';
    else
      StopFlag = 0;
      Status   = [];
    end

  case 9
    if (nhxk <= epsilon || Niter >= MaxNumIter)
      StopFlag = 1;
      if Niter < MaxNumIter
          Status = 'a solution is found.';
      else
          Status = 'Maximum number of iterations is reached.';
      end
    else
      StopFlag = 0;
      Status   = [];
    end

end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% End of StopCriterion.m %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


text:
%% Machine Learning Online Class
%  Exercise 5 | Regularized Linear Regression and Bias-Variance
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  exercise. You will need to complete the following functions:
%
%     linearRegCostFunction.m
%     learningCurve.m
%     validationCurve.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% =========== Part 1: Loading and Visualizing Data =============
%  We start the exercise by first loading and visualizing the dataset.
%  The following code will load the dataset into your environment and plot
%  the data.
%

% Load Training Data
fprintf('Loading and Visualizing Data ...\n')

% Load from ex5data1:
% You will have X, y, Xval, yval, Xtest, ytest in your environment
load ('ex5data1.mat');

% m = Number of examples
m = size(X, 1);

% Plot training data
plot(X, y, 'rx', 'MarkerSize', 10, 'LineWidth', 1.5);
xlabel('Change in water level (x)');
ylabel('Water flowing out of the dam (y)');

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =========== Part 2: Regularized Linear Regression Cost =============
%  You should now implement the cost function for regularized linear
%  regression.
%

theta = [1 ; 1];
J = linearRegCostFunction([ones(m, 1) X], y, theta, 1);

fprintf(['Cost at theta = [1 ; 1]: %f '...
         '\n(this value should be about 303.993192)\n'], J);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =========== Part 3: Regularized Linear Regression Gradient =============
%  You should now implement the gradient for regularized linear
%  regression.
%

theta = [1 ; 1];
[J, grad] = linearRegCostFunction([ones(m, 1) X], y, theta, 1);

fprintf(['Gradient at theta = [1 ; 1]:  [%f; %f] '...
         '\n(this value should be about [-15.303016; 598.250744])\n'], ...
         grad(1), grad(2));

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =========== Part 4: Train Linear Regression =============
%  Once you have implemented the cost and gradient correctly, the
%  trainLinearReg function will use your cost function to train
%  regularized linear regression.
%
%  Write Up Note: The data is non-linear, so this will not give a great
%                 fit.
%

%  Train linear regression with lambda = 0
lambda = 1;
[theta] = trainLinearReg([ones(m, 1) X], y, lambda);

%  Plot fit over the data
plot(X, y, 'rx', 'MarkerSize', 10, 'LineWidth', 1.5);
xlabel('Change in water level (x)');
ylabel('Water flowing out of the dam (y)');
hold on;
plot(X, [ones(m, 1) X]*theta, '--', 'LineWidth', 2)
hold off;

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =========== Part 5: Learning Curve for Linear Regression =============
%  Next, you should implement the learningCurve function.
%
%  Write Up Note: Since the model is underfitting the data, we expect to
%                 see a graph with "high bias" -- slide 8 in ML-advice.pdf
%

lambda = 1;
[error_train, error_val] = ...
    learningCurve([ones(m, 1) X], y, ...
                  [ones(size(Xval, 1), 1) Xval], yval, ...
                  lambda);

plot(1:m, error_train, 1:m, error_val);
title('Learning curve for linear regression')
legend('Train', 'Cross Validation')
xlabel('Number of training examples')
ylabel('Error')
axis([0 13 0 150])

fprintf('# Training Examples\tTrain Error\tCross Validation Error\n');
for i = 1:m
    fprintf('  \t%d\t\t%f\t%f\n', i, error_train(i), error_val(i));
end

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =========== Part 6: Feature Mapping for Polynomial Regression =============
%  One solution to this is to use polynomial regression. You should now
%  complete polyFeatures to map each example into its powers
%

p = 8;

% Map X onto Polynomial Features and Normalize
X_poly = polyFeatures(X, p);
[X_poly, mu, sigma] = featureNormalize(X_poly);  % Normalize
X_poly = [ones(m, 1), X_poly];                   % Add Ones

% Map X_poly_test and normalize (using mu and sigma)
X_poly_test = polyFeatures(Xtest, p);
X_poly_test = bsxfun(@minus, X_poly_test, mu);
X_poly_test = bsxfun(@rdivide, X_poly_test, sigma);
X_poly_test = [ones(size(X_poly_test, 1), 1), X_poly_test];         % Add Ones

% Map X_poly_val and normalize (using mu and sigma)
X_poly_val = polyFeatures(Xval, p);
X_poly_val = bsxfun(@minus, X_poly_val, mu);
X_poly_val = bsxfun(@rdivide, X_poly_val, sigma);
X_poly_val = [ones(size(X_poly_val, 1), 1), X_poly_val];           % Add Ones

fprintf('Normalized Training Example 1:\n');
fprintf('  %f  \n', X_poly(1, :));

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% =========== Part 7: Learning Curve for Polynomial Regression =============
%  Now, you will get to experiment with polynomial regression with multiple
%  values of lambda. The code below runs polynomial regression with
%  lambda = 0. You should try running the code with different values of
%  lambda to see how the fit and learning curve change.
%

lambda = 1;
[theta] = trainLinearReg(X_poly, y, lambda);

% Plot training data and fit
figure(1);
plot(X, y, 'rx', 'MarkerSize', 10, 'LineWidth', 1.5);
plotFit(min(X), max(X), mu, sigma, theta, p);
xlabel('Change in water level (x)');
ylabel('Water flowing out of the dam (y)');
title (sprintf('Polynomial Regression Fit (lambda = %f)', lambda));

figure(2);
[error_train, error_val] = ...
    learningCurve(X_poly, y, X_poly_val, yval, lambda);
plot(1:m, error_train, 1:m, error_val);

title(sprintf('Polynomial Regression Learning Curve (lambda = %f)', lambda));
xlabel('Number of training examples')
ylabel('Error')
axis([0 13 0 100])
legend('Train', 'Cross Validation')

fprintf('Polynomial Regression (lambda = %f)\n\n', lambda);
fprintf('# Training Examples\tTrain Error\tCross Validation Error\n');
for i = 1:m
    fprintf('  \t%d\t\t%f\t%f\n', i, error_train(i), error_val(i));
end

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =========== Part 8: Validation for Selecting Lambda =============
%  You will now implement validationCurve to test various values of
%  lambda on a validation set. You will then use this to select the
%  "best" lambda value.
%

[lambda_vec, error_train, error_val] = ...
    validationCurve(X_poly, y, X_poly_val, yval);

close all;
plot(lambda_vec, error_train, lambda_vec, error_val);
legend('Train', 'Cross Validation');
xlabel('lambda');
ylabel('Error');

fprintf('lambda\t\tTrain Error\tValidation Error\n');
for i = 1:length(lambda_vec)
	fprintf(' %f\t%f\t%f\n', ...
            lambda_vec(i), error_train(i), error_val(i));
end

fprintf('Program paused. Press enter to continue.\n');
pause;

author: vugsus

repo: coursera-machine-learning


text:
function [m2] = ft22m2(ft2)
% Convert area from square feet to square meters.
% Chad A. Greene 2012
m2 = ft2*0.09290304;


text:
function L = co_reggui(action,varargin)
% GUI for Manual co-registration of different image modalities to go with
% CERR
% - Detailed Description: this is a manual co-registration routine for two
% set of images (2d/3d) assuming that a possible resmapling and affine
% transformation are sufficient for alignment. The quality of registration is measured by
% using the mutual information criterion (MI)
%
% Written By: Issam El Naqa    Date: 08/28/03
% Revised by:                  Date:
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

if nargin<1,
    action='Initializeco_reggui';
end;

switch action

    case 'Initializeco_reggui'   %Initialization
        Initializeco_reggui;

    case 'RegisterMode'
        mode = get(findobj('Tag', 'dispregpopupmenu'),'Value');
        switch mode
            case 1 %manual
            case 2 %ctrl points
            case 3 %auto control points
        end

    case 'LoadRefImage'
        set(findobj(gcbf, 'Tag', 'StatusText'), 'String', 'Loading images...');
        h = get(gcbf,'Userdata');
        pathname 	= h.pathname;
        [filename, pathname] = uigetfile('*.*', 'Load reference image');
        load([pathname,filename]);
        x = cropImage(x);
        h.refimage=min_maxnorm(x, 0:255);
        set(gcbf,'Userdata',h);
        DispImage('axref',x);
        drawnow

    case 'LoadPlanCUnregImage'

    case 'LoadPlanCRefImage'

    case 'LoadUnregImage'
        set(findobj(gcbf, 'Tag', 'StatusText'), 'String', 'Loading images...');
        h = get(gcbf,'Userdata');
        pathname 	= h.pathname;
        [filename, pathname] = uigetfile('*.*', 'Load unregistered image');
        load([pathname,filename]);
        x = cropImage(x);
        h.unregimage=min_maxnorm(x, 0:255);
        h.regimage=h.unregimage; % initailize registered image to unregistered one!
        set(gcbf,'Userdata',h);
        DispImage('axunreg',x);
        drawnow
        set(findobj(gcbf, 'Tag', 'StatusText'), 'String', 'Set registration parameters...');
    case 'EditRegParam'
        h = get(gcbf,'Userdata');
        h.Tx=str2num(get(findobj(gcbf,'Tag','TXed'),'String'));
        h.Ty=str2num(get(findobj(gcbf,'Tag','TYed'),'String'));
        h.Tz=str2num(get(findobj(gcbf,'Tag','TZed'),'String'));
        h.rot=str2num(get(findobj(gcbf,'Tag','Roted'),'String'));
        h.Sx=str2num(get(findobj(gcbf,'Tag','SXed'),'String'));
        h.Sy=str2num(get(findobj(gcbf,'Tag','SYed'),'String'));
        h.Sz=str2num(get(findobj(gcbf,'Tag','SZed'),'String'));
        h.samp=str2num(get(findobj(gcbf,'Tag','Resamped'),'String'));
        set(findobj(gcbf, 'Tag', 'regpushbutton'), 'Enable', 'on');
        set(findobj(gcbf, 'Tag', 'dispmodepopupmenu'), 'Enable', 'on');
        set(findobj(gcbf, 'Tag', 'StatusText'), 'String', 'Press register button for processing...');
        set(gcbf,'UserData',h);
    case 'DisplayMode'
        setDisplayMode;

    case 'enlargeOverlay'
        h = get(gcbf,'Userdata');
        displayOverlays(h.refimage, h.regimage);

    case 'ApplyRegistration'
        h = get(gcbf,'Userdata');

        mode = get(findobj('Tag', 'dispregpopupmenu'),'Value');
        switch mode
            case 1 %manual
                tempimage=imtranslate2d(h.unregimage,h.Tx,h.Ty);
                tempimage=imrotate2d(tempimage,h.rot);
                tempimage=imscale2d(tempimage,h.Sx,h.Sy);
                tempimage=imresample2d(tempimage,h.samp);
                h.regimage=tempimage;
            case 2 %ctrl points
                getControlPoints('init');
                getControlPoints('load', h.refimage, h.unregimage);
                waitfor(findobj('Tag', 'RegGui'), 'Tag', 'RegGuiDone');
                [ref_pts, target_pts] = getControlPoints('getpoints', h.refimage, h.unregimage);
                delete(findobj('Tag', 'RegGuiDone'));
                [tempimage,A]=compute_aff_transform(h.unregimage, ref_pts, target_pts);
                h.regimage=tempimage;
            case 3
                getControlPoints('init');
                getControlPoints('load', h.refimage, h.unregimage);
                waitfor(findobj('Tag', 'RegGui'), 'Tag', 'RegGuiDone');
                [ref_pts, target_pts] = getControlPoints('getpoints', h.refimage, h.unregimage);
                delete(findobj('Tag', 'RegGuiDone'));
                [tempimage,A]=compute_perspect_transform(h.unregimage, ref_pts, target_pts);
                h.regimage=tempimage;
            case 4 %auto control points
                numPoints = 20;
                sigma = 1;
                get_control_points(h.unregimage, h.refimage, numPoints, sigma);

        end

        % make joint histogram plot
        [h.mi, h.jhist]=get_mutualinfo(h.refimage,h.regimage);
        %set(findobj(gcbf, 'Tag', 'mitext'), 'Visible', 'on');
        %set(findobj(gcbf, 'Tag', 'mivalue'), 'Visible', 'on');
        % make a cross-correlation image
        h.xcorr=fxcorr(h.refimage,h.regimage);
        % save current data
        set(gcbf,'UserData',h);
        % display images
        setDisplayMode;
        DisplayJHist(h.jhist);
        set(findobj(gcbf, 'Tag', 'mivalue'), 'String', num2str(h.mi));
        DispImage('axcorr',h.xcorr);
        drawnow
    case 'info'
        helpwin co_reggui;
    case 'close'
        close(gcbf);
end

return

% supplmentary routines

function setDisplayMode() % set display mode for registration
h = get(gcbf,'Userdata');
h.displaymode = get(findobj(gcbf,'Tag','dispmodepopupmenu'),'Value');
switch h.displaymode
    case 1  % Single registered
        DispImage('axreg',h.regimage);
    case 2  % Overlayed images
        alpha=0.8; % transperancy factor
        alternate_time=0.5; % refresh time
        AlternateOverlayed(h.refimage,h.regimage,'axreg',alpha,alternate_time);
    case 3  % Sliceomatic
end
set(gcbf,'UserData',h);
return

% image display function
function DispImage(imTag,x)
set(gcbf,'CurrentAxes',findobj(gcbf,'Tag',imTag));
cla, h=imagesc(x), axis image, colormap('hot'), axis ij, axis off
set(h, 'ButtonDownFcn', ['co_reggui(''Ctrl_' imTag ''');']);
return

% image translation function
function fh=imtranslate2d(f,xoff,yoff)
[h,w]=size(f);
[x,y]=meshgrid([1:1:h],[1:1:w]);
xd=x(:)+xoff; yd=y(:)+yoff;
fh=reshape(bilinear_interpolation(f,xd,yd),w,h)';
return

% image rotation function around the middle of the image
function fh=imrotate2d(f,ang)
[h,w]=size(f);
phi = ang*pi/180; % Convert to radians
vx=[-floor(h/2):ceil(h/2)-1]; % center around the middle
vy=[-floor(w/2):ceil(w/2)-1];

[x,y]=meshgrid(vx,vy);
x=x(:); y=y(:);
xd=x*cos(phi)+y*sin(phi)+floor(h/2)+1;
yd=-x*sin(phi)+y*cos(phi)+floor(w/2)+1;
fh=reshape(bilinear_interpolation(f,xd,yd),w,h)';
return

% image scaling function (exapnsion/shrinking)
function fh=imscale2d(f,Sx,Sy)
[h,w]=size(f);
[x,y]=meshgrid([1:1:h],[1:1:w]);
xd=x(:)*Sx; yd=y(:)*Sy;
fh=reshape(bilinear_interpolation(f,xd,yd),w,h)';
return

% image resampling function
function fh=imresample2d(f,q)
[h,w]=size(f);
% use a binomial filter of order 5 for smoothing
d=[1 4 6 4 1]/16; B=d'*d;
fh=f;
if q==1
    return
elseif q>1 % upsample
    N=round(q);
    y = zeros(N*h,N*w);
    y(1:N:end,1:N:end)=fh;
    fh=conv2(y,B,'same');
elseif (q<1 & q>0) % downsample
    D=round(1/q);
    fh=conv2(fh,B,'same');
    fh=fh(1:D:end,1:D:end);
else
    errordlg('The resampling factor should a positive number!', 'co_reggui Error', 'replace');
end

% perform fast cross-correlation in frequency domain
function c=fxcorr(x,y)
Fsize=size(x)+size(y)-1;
Fx = fft2(rot90(x,2),Fsize(1),Fsize(2));
Fy = fft2(y,Fsize(1),Fsize(2));
c = real(ifft2(Fx .* Fy));
return

% compute the mutual information using the joint histogram
function [mi, histxy]=get_mutualinfo(x,y)
% x, y : the two images,
siz=min([size(x);size(y)]); % if sizes are different
nbits=8; ngray=2^nbits; % assume 256 levels is sufficient approximation!
x=double(uint8(double(x)+1)); y=double(uint8(double(y)+1)); % convert to 8 bits
histxy=zeros(ngray,ngray);

[iM,jM] = meshgrid(1:siz(1),1:siz(2));

indV = (jM(:) - 1) * siz(1) + iM(:);

xV =double(x(indV));
yV = double(y(indV));

ind2V = (yV - 1) * ngray + xV;

for i=1:length(ind2V)
    histxy(ind2V(i)) = histxy(ind2V(i)) + 1;
end

%for i=1:siz(1)
%    for j=1:siz(2)
%        histxy(x(i,j),y(i,j))= histxy(x(i,j),y(i,j))+1;
%    end
%end

histxy=histxy/sum(histxy(:)); % normalize
% compute marginal distributions
histx=sum(histxy,2); histy=sum(histxy,1); % by integrating out the joint
mi=sum(sum(histxy.*log2(histxy./(histx*histy+eps)+eps)));
return

% joint histogram display
function DisplayJHist(jhist)
set(gcbf,'CurrentAxes',findobj(gcbf,'Tag','axjhist'));
set(findobj(gcbf,'Tag','axjhist'),'Box','off');
cla, view(-37.5,30), colormap('hot'), mesh(jhist);
return

% display transparent in a cyclic fashion
function DisplayOverlayed(x,y,imTag,alpha)
% try alphas, or linear combinations...
%mask = ones(size(y));
%mask(find(y<1)) = alpha;
set(gcbf,'CurrentAxes',findobj(gcbf,'Tag',imTag));
sx=size(x); sy=size(y); % linear combination
siz=max([sx;sy]);
xa=zeros(siz); xa(1:sx(1),1:sx(2))=x;
ya=zeros(siz); ya(1:sy(1),1:sy(2))=y;
cla, imagesc(double(xa)+0.75*double(ya), 'ButtonDownFcn', 'co_reggui(''enlargeOverlay'');'), axis image, colormap('hot'), axis ij, axis off
% hold on
% hi=imagesc();
% axis image, colormap('hot'), axis ij, axis off;
%set(hi,'AlphaData',mask);
return

function AlternateOverlayed(x,y,imTag,alpha, atime)
DisplayOverlayed(x,y,imTag,alpha);
% pause(atime);
% DisplayOverlayed(y,x,imTag,alpha);
return

function x = cropImage(x)
    minCol = 1, minRow = 1;
    [maxCol, maxRow] = size(x);
    rows = find(~max(x'));
    cols = find(~max(x));
    for i=1:length(rows)
        if rows(i) ~= i
            minRow = rows(i-1)+1;
            maxRow = rows(i)-1;
            break;
        end
    end
    for i=1:length(cols)
        if cols(i) ~= i
            minCol = cols(i-1)+1;
            maxCol = cols(i)-1;
            break;
        end
    end
    x = imcrop(x,[minCol, minRow,  maxCol-minCol, maxRow-minRow]);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [out, W]= proc_localAverageReference(dat, mnt, varargin)
%PROC_LOCALAVERAGEREFERENCE - Rereference all channels to local average
%
%Synopsis:
%  [DAT, W]= proc_localAverageReference(DAT, MNT, RADIUS)
%  [DAT, W]= proc_localAverageReference(DAT, MNT, <OPT>)
%
%Arguments:
%  DAT: Data structure of continuous or epoched data
%  MNT: Electrode montage, see getElectrodePositions
%  RADIUS: see OPT.radius
%  OPT: struct or property/value list of optional properties:
%    .radius: For a radius of 0.8, the neighborhood of Cz extends from
%      C3 to C4 and from Fz to Pz. For a radius of 0.4 it extends from
%      C1 to C2 and from FCz to CPz. Default 0.8.
%    .clab: Labels of channels which are to be retrieved.
%    .median: If true, reference is calculated as median across
%      neighboring channels, instead of mean. Default 0.
%    .verbose: If true, text output shows the rereferencing
%
%Returns:
%  DAT: updated data structure
%  W:   filter matrix that can be used, e.g. in proc_linearDerivation
%       (before you need to exclude the same channels
%
%Description:
% Rereference signals to local average, i.e. to the average of all
% electrodes within a given radius.
%
%SEE:
%  proc_commonAverageReference, proc_laplacian, getElectrodePositions

% Author: Benjamin Blankertz
dat = misc_history(dat);

props= { 'radius'       0.8     'DOUBLE[1]'
         'clab'         '*'     'CHAR|CELL{CHAR}'
         'ignore_clab'  {'E*'}  'CHAR|CELL{CHAR}'
         'median'       0       'BOOL'
         'verbose'      0       'BOOL'};

if nargin==0,
  out = props; return
end

misc_checkType(dat, 'STRUCT(x clab)');
if length(varargin)==1 & isnumeric(varargin{1}),
  opt= struct('radius', varargin{1});
else
  opt= opt_proplistToStruct(varargin{:});
end
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

if ~iscell(opt.ignore_clab),
  opt.ignore_clab= {opt.ignore_clab};
end

mnt= mnt_adaptMontage(mnt, dat);
if ~isequal(mnt.clab, dat.clab),
  error('channel mismatch');
end
rc= util_chanind(dat, {'not', opt.ignore_clab{:}});
idx_tbf= util_chanind(dat, opt.clab);
out= proc_selectChannels(dat, opt.clab);
W= zeros(length(dat.clab), length(idx_tbf));
for ci= 1:length(idx_tbf),
  cc= idx_tbf(ci);
  pos= repmat(mnt.pos_3d(:,cc), [1 length(rc)]);
  dist= sqrt(sum( (mnt.pos_3d(:,rc)-pos).^2) );
  iRef= find(dist>0 & dist<opt.radius);
  if opt.verbose,
    fprintf('%s: ref''ed to: %s\n', dat.clab{cc}, str_vec2str(dat.clab(rc(iRef))));
  end
  W(cc,ci)= 1;
  if ~isempty(iRef),
    if opt.median,
      lar= median(dat.x(:,rc(iRef),:), 2);
    else
      lar= mean(dat.x(:,rc(iRef),:), 2);
    end
    out.x(:,ci,:)= dat.x(:,cc,:) - lar;
    out.clab{ci}= [dat.clab{cc} ' lar'];
    W(rc(iRef),ci)= -1/length(iRef);
  end
end

author: bbci


text:
function [ ur, vr, pr ] = resid_stokes1 ( n, x, y )

%*****************************************************************************80
%
%% RESID_STOKES1 returns residuals of the exact Stokes solution #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Junping Wang, Yanqiu Wang, Xiu Ye,
%    A robust numerical method for Stokes equations based on divergence-free
%    H(div) finite element methods,
%    SIAM Journal on Scientific Computing,
%    Volume 31, Number 4, 2009, pages 2784-2802.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), Y(N), the coordinates of the points.
%
%    Output, real UR(N), VR(N), PR(N), the residuals in the U,
%    V and P equations.
%
  ur = zeros ( n, 1 );
  vr = zeros ( n, 1 );
  pr = zeros ( n, 1 );
%
%  Get the right hand sides.
%
  [ f, g, h ] = rhs_stokes1 ( n, x, y );
%
%  Form the functions and derivatives.
%
  for i = 1 : n

    u = - 2.0 ...
          * x(i) ^ 2 * ( x(i) - 1.0 ) ^ 2  ...
          * y(i) * ( y(i) - 1.0 ) * ( 2.0 * y(i) - 1.0 );

    ux = - 2.0  ...
          * ( 4.0 * x(i) ^ 3 - 6.0 * x(i) ^ 2  ...
          + 2.0 * x(i) ) ...
          * y(i) * ( y(i) - 1.0 ) * ( 2.0 * y(i) - 1.0 );

    uxx = - 2.0  ...
          * ( 12.0 * x(i) ^ 2 - 12.0 * x(i) + 2.0 ) ...
          * ( 2.0 * y(i) ^ 3 - 3.0 * y(i) ^ 2 + y(i) );

    uy = - 2.0  ...
          * x(i) ^ 2 * ( x(i) - 1.0 ) ^ 2  ...
          * ( 6.0 * y(i) ^ 2 - 3.0 * y(i) + 1.0 );

    uyy = - 2.0  ...
          * ( x(i) ^ 4 - 2.0 * x(i) ^ 3 + x(i) ^ 2 ) ...
          * ( 12.0 * y(i) - 6.0 );

    v =   2.0  ...
          * x(i) * ( x(i) - 1.0 ) * ( 2.0 * x(i) - 1.0 ) ...
          * y(i) ^ 2 * ( y(i) - 1.0 ) ^ 2;

    vx =   2.0  ...
          * ( 6.0 * x(i) ^ 2 - 6.0 * x(i) + 1.0 ) ...
          * y(i) ^ 2 * ( y(i) - 1.0 ) ^ 2;

    vxx =   2.0  ...
          * ( 12.0 * x(i) - 6.0 ) ...
          * y(i) ^ 2 * ( y(i) - 1.0 ) ^ 2;

    vy =   2.0  ...
          * x(i) * ( x(i) - 1.0 ) * ( 2.0 * x(i) - 1.0 ) ...
          * ( 4.0 * y(i) ^ 3 - 6.0 * y(i) ^ 2  ...
          + 2.0 * y(i) );

    vyy =   2.0  ...
          * x(i) * ( x(i) - 1.0 ) * ( 2.0 * x(i) - 1.0 ) ...
          * ( 12.0 * y(i) ^ 2 - 12.0 * y(i) + 2.0 );

    p = 0.0;
    px = 0.0;
    py = 0.0;

    ur(i) = px - ( uxx + uyy ) - f(i);
    vr(i) = py - ( vxx + vyy ) - g(i);
    pr(i) = ux + vy - h(i);

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function exactness_test02 ( )

%*****************************************************************************80
%
%% EXACTNESS_TEST02 tests Fejer Type 2 rules for Legendre integrals.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EXACTNESS_TEST02\n' );
  fprintf ( 1, '  Test Fejer Type 2 rules on Legendre integrals.\n' );
  fprintf ( 1, '  Density function rho(x) = 1.\n' );
  fprintf ( 1, '  Region: -1 <= x <= +1.\n' );
  fprintf ( 1, '  Exactness: N   for N odd,\n' );
  fprintf ( 1, '             N-1 for N even.\n' );

  for n = 1 : 5

    [ x, w ] = fejer2_set ( n );
    if ( mod ( n, 2 ) == 1 )
      p_max = n + 1;
    else
      p_max = n;
    end
    legendre_exactness ( n, x, w, p_max );

  end

  return
end