text:
function varargout = plot(varargin)
%PLOT   Basic linear plot for CHEBFUN objects.
%   PLOT(F) plots the CHEBFUN object F in the interval where it is defined. If F
%   is complex valued, PLOT(F) is equivalent to PLOT(real(F), imag(F)).
%
%   PLOT(F, G) plots the CHEBFUN G versus the CHEBFUN F. Quasimatrices and
%   array-valued CHEBUFUN objects are also supported in the natural way.
%
%   PLOT(F, S) allows various line types, plot symbols, and colors to be used
%   when S is a character string made from one element from any or all the
%   following 3 columns:
%
%            b     blue          .     point              -     solid
%            g     green         o     circle             :     dotted
%            r     red           x     x-mark             -.    dashdot
%            c     cyan          +     plus               --    dashed
%            m     magenta       *     star             (none)  no line
%            y     yellow        s     square
%            k     black         d     diamond
%            w     white         v     triangle (down)
%                                ^     triangle (up)
%                                <     triangle (left)
%                                >     triangle (right)
%                                p     pentagram
%                                h     hexagram
%
%   The entries from the centre columns are plotted at the grid being used to
%   represent F (typically Chebyshev).
%
%   The X,Y pairs, or X,Y,S triples, can be followed by parameter/value pairs to
%   specify additional properties of the lines. For example,
%            f = chebfun(@sin);
%            plot(f, 'LineWidth', 2, 'Color', [.6 0 0])
%   will create a plot with a dark red line width of 2 points.
%
%   PLOT(F1, S1, F2, S2, F3, S3, ...) or PLOT(F1, G1, S1, F2, G2, S2, ...)
%   combines the plots defined by the (F,G,S) triples or (F,S) doubles, where
%   the F's and G's are CHEBFUN object and the S's are strings.
%
%   [HLINE, HPOINT, HJUMP] = PLOT(F) returns column vectors of handles to
%   lineseries objects (one for each column in the case of array-valued CHEBFUN
%   objects), corresponding to the line, point, and jump plots, respectively.
%
%   PLOT(F, 'interval', [A, B]) restricts the plot to the interval [A, B], which
%   can be useful when the domain of F is infinite, or for 'zooming in' on, say,
%   oscillatory CHEBFUN objects. If plotting an array-valued CHEBFUN or more
%   than one CHEBFUN in a call like PLOT(F, 'b', G, '--r', 'interval', [A, B])
%   this property is applied globally. Markers, such as 'o', or '.', are ignored
%   if the interval flag is used.
%
%   Besides the usual parameters that control the specifications of lines (see
%   linespec), the parameter JumpLine and DeltaLine determines the linestyle for
%   the discontinuities and the delta functions of the CHEBFUN F, respectively.
%   For example, PLOT(F, 'JumpLine', '-r') will plot discontinuities as solid
%   red lines and PLOT(F, 'deltaline', '-r') will plot the delta functions as
%   solid red lines. By default the plotting style for JumpLine is ':', and '-'
%   for delta functions and colours are chosen to match the lines they
%   correspond to. It is possible to modify other properties of JumpLines syntax
%   like PLOT(F, 'JumpLine', {'color', 'r', 'LineWidth', 5}). JumpLines and
%   deltaLines can be suppressed with the argument 'none'.
%
%   Note that the PLOT(F, 'numpts', N) option for V4 is deprecated, and this
%   call now has no effect.
%
% See also PLOTDATA, PLOT3.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEVELOPER NOTE:
%  We actually plot a fifth 'dummy' plot which contains both the line and marker
%  styles. This is the only plot which has 'handlevis', 'on', and so will be
%  included in subsequent calls to LEGEND(). The handle to this plot is included
%  as the fifth output, and the plot itself is set to 'visible', 'off'. Note
%  that if a user modifies the plot by adjusting HLINE or HPOINT, subsequent
%  calls the LEGEND will not have the updated style. There's no way to get
%  around this without a complicated callback procedure.
%
%  Note, to clarify, the reason we need all this is that we must plot the line
%  and points as different plots. MATLAB/PLOT() doesn't have this problem.
%
%  Also note that the handles are now output as individually, rather than
%  forced in to an array.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% TODO: Document plotting of DELTFAFUN objects. (In particular, handle output)

%% Deal with an empty input:
if ( isempty(varargin{1}) )
    if ( nargout == 1 )
        varargout{1} = plot([]);
    end
    return
end

%% Initialization:

% Suppress inevitable warning for growing these arrays:
%#ok<*AGROW>

% Store the hold state of the current axis:
holdState = ishold;

% Store the current X and Y-limits, and whether ylim is in manual or auto mode.
if ( holdState )
    xLimCurrent = get(gca, 'xlim');
    yLimCurrent = get(gca, 'ylim');
    yLimModeCurrent = get(gca, 'ylimmode');
end

% Initialize flags:
isComplex = false;
intervalIsSet = false;

% Initialize XLIM and YLIM. Note that the first entries are initialized to be
% UPPER limits on the LOWER parts of the axes, while the second entries
% correspond to LOWER bounds on the UPPER parts of the axes. Hence, this
% somewhat strange convention.
xLim = [inf, -inf];
yLim = [inf, -inf];
defaultXLim = 1;
defaultYLim = 1;

% Check to see if the 'interval' flag has been set:
interval = [];
loc = find(strcmpi(varargin, 'interval'));
if ( any(loc) )
    interval = varargin{loc+1};
    varargin(loc:loc+1) = [];
    intervalIsSet = true;
else
    % TODO: Do we want to support this?
    % Allow plot(f, [a, b]) as shorthand for plot(f, 'interval', [a, b]):
    loc = find(cellfun(@(f) isa(f, 'chebfun'), varargin));
    for k = 1:numel(loc)
        if ( loc(k) < nargin && isnumeric(varargin{loc(k)+1}) )
            interval = varargin{loc(k)+1};
            varargin(loc(k)+1) = [];
            intervalIsSet = true;
            break
        end
    end
end

% Support 'interval' by evaluating on a fixed grid size (2000 points). See #602.
if ( intervalIsSet )
    for k = 1:numel(varargin)
        if ( isa(varargin{k}, 'chebfun') )
            dom = union(domain(varargin{k}), interval);
            dom(dom < interval(1) | dom > interval(end)) = [];
            varargin{k} = chebfun(@(x) feval(varargin{k}, x), dom, 2000);
        end
    end
end

% Initialise storage:
lineData = {};
pointData = {};
jumpData = {};
deltaData = {};
jumpLineIsSet = any(cellfun(@(v) ischar(v) && strcmpi(v, 'JumpLine'), varargin));

% Remove global plotting options from input arguments.
[lineStyle, pointStyle, jumpStyle, deltaStyle, varargin] = ...
    chebfun.parsePlotStyle(varargin{:});

%% Preparation of the data:

% Get the data for plotting from PLOTDATA():
while ( ~isempty(varargin) )

    % Acquire plotting data for each CHEBFUN / pair of CHEBFUNs:
    if ( (numel(varargin) > 1) && ... % PLOT(f, g).
            (isa(varargin{2}, 'chebfun') || isnumeric(varargin{2})) )
        % Remove CHEBFUN objects from array input:
        f = varargin{1};
        g = varargin{2};
        varargin(1:2) = [];

        % We can only plot real against real:
        isComplex = false;
        if ( ~isreal(f) || ~isreal(g) )
            warning('CHEBFUN:CHEBFUN:plot:complex', ...
                'Imaginary parts of complex X and/or Y arguments ignored.');
            f = real(f);
            g = real(g);
        end

        % Call PLOTDATA():
        if ( isnumeric(f) )
            % For plotting numeric data in a CHEBFUN/PLOT() call.
            newData = plotData(chebfun());
            newData.xLine = f;
            newData.yLine = g;
            newData.xPoints = f;
            newData.yPoints = g;
            newData.xJumps = NaN;
            newData.yJumps = NaN;
            newData.xDeltas = NaN;
            % Do nothing
        elseif ( numel(f) == 1 && numel(g) == 1 )
            % Array-valued CHEBFUN case:
            newData = plotData(f, g);
        else
            % QUASIMATRIX case:
            f = num2cell(f);
            g = num2cell(g);
            if ( numel(f) > 1 && numel(g) > 1 )
                if ( numel(f) ~= numel(g) )
                    error('CHEBFUN:CHEBFUN:plot:dim', ...
                    'CHEBFUN objects must have the same number of columns.');
                end
                for k = 1:numel(f)
                    newData(k) = plotData(f{k}, g{k});
                end
            elseif ( numel(f) > 1 && numel(g) == 1 )
                for k = 1:numel(f)
                    newData(k) = plotData(f{k}, g{1});
                end
            elseif ( numel(f) == 1 && numel(g) > 1 )
                for k = 1:numel(f)
                    newData(k) = plotData(f{1}, g{k});
                end
            end
        end

    else  % PLOT(f).

        % Remove CHEBFUN from array input:
        f = varargin{1};
        varargin(1) = [];
        isComplex = ~isreal(f);

        % Loop over the columns:
        for k = 1:numel(f)
            if ( isComplex )
                newData(k) = plotData(real(f(k)), imag(f(k)));
            else
                newData(k) = plotData(f(k));
            end
        end

    end

    % Style data.
    pos = 0;
    styleData = [];

    lv = length(varargin);
    % Find the location of the next CHEBFUN in the input array:
    while ( (pos < lv) && ~isa(varargin{pos+1}, 'chebfun') )
        if ( pos+1 < lv && isnumeric(varargin{pos+1}) && ...
                isnumeric(varargin{pos+2}) )
            break
        end
        pos = pos + 1;
    end

    if ( pos > 0 )
        styleData = varargin(1:pos);
        varargin(1:pos) = [];
        % Remove deprecated 'numpts' option:
        idx = find(strcmp(styleData, 'numpts'), 1);
        if ( any(idx) )
            styleData(idx:(idx+1)) = [];
        end
    end

    % Loop over the columns:
    for k = 1:numel(newData)

        % Update axis limits:
        xLim = [min(newData(k).xLim(1), xLim(1)), ...
            max(newData(k).xLim(2), xLim(2))];
        yLim = [min(newData(k).yLim(1), yLim(1)), ...
            max(newData(k).yLim(2), yLim(2))];

        % Append new data:
        lineData = [lineData, newData(k).xLine, newData(k).yLine, styleData];
        pointData = [pointData, newData(k).xPoints, newData(k).yPoints, ...
            styleData];
        jumpData = [jumpData, newData(k).xJumps, newData(k).yJumps, styleData];
        deltaData{k} = {newData(k).xDeltas, newData(k).yDeltas, ...
            newData(k).yDeltaBase, styleData};

        defaultXLim = defaultXLim & newData(k).defaultXLim;
        defaultYLim = defaultYLim & newData(k).defaultYLim;
    end

    % If xLim(1) == xLim(2), let Matlab figure out a proper xLim:
    if ( ~diff(xLim) )
        defaultXLim = 1;
    end

    % If yLim(1) == yLim(2), let Matlab figure out a proper yLim:
    if ( ~diff(yLim) )
        defaultYLim = 1;
    end

end

%% Plotting starts here:

% Acquire initial color cycle if running R2014b+.
if ( ~verLessThan('matlab', '8.4') )
    if ( ~holdState )
        set(gca, 'ColorOrderIndex', 1);
    end
    originalColorOrder = get(gca, 'ColorOrderIndex');
end

% Plot the lines:
h1 = plot(lineData{:});
set(h1, 'Marker', 'none', lineStyle{:})

% Ensure the plot is held:
hold on

% Get color cycle prior to point plot if running R2014b.
if ( ~verLessThan('matlab', '8.4') )
    newColorOrder = get(gca, 'ColorOrderIndex');
    set(gca, 'ColorOrderIndex', originalColorOrder)
end

% Plot the points:
h2 = plot(pointData{:});
% Change the style accordingly:
set(h2, 'LineStyle', 'none', pointStyle{:})
if ( intervalIsSet )
    % Markers are meaningless if the 'interval' flag is used.
    set(h2, 'Marker', 'none', pointStyle{:})
end

% Reset color cycle prior to jump plot if running R2014b.
if ( ~verLessThan('matlab', '8.4') )
    set(gca, 'ColorOrderIndex', originalColorOrder);
end

% Plot the jumps:
if ( isempty(jumpData) || ischar(jumpData{1}) )
    jumpData = {[]};
end
h3 = plot(jumpData{:});
% Change the style accordingly:
if ( ~isempty(jumpStyle) )
    set(h3, jumpStyle{:});
end
if ( ~jumpLineIsSet )
    if ( isComplex )
        set(h3, 'LineStyle', 'none')
    else
        set(h3, 'LineStyle', ':')
    end
    set(h3, 'Marker', 'none')
end

% Plot the Delta functions:
if ( isempty(deltaData) || ~isnumeric(deltaData{1}{1}) )
    h4 = plot([]);
else
    h4 = plotDeltas(deltaData);
end
if ( ~isempty(deltaStyle) )
    set(h4, deltaStyle{:});
end

% Reset colors prior to legend data plot if running R2014b.
if ( ~verLessThan('matlab', '8.4') )
    set(gca, 'ColorOrderIndex', originalColorOrder);
end

% Generate data for the dummy plot for the legend.  dummyData is the same as
% lineData but with the numeric entries replaced by NaNs of suitable dimension
% so that they don't show up on the plot but are still "active" for the legend.
dummyData = lineData;
for i = 1:1:length(dummyData)
    if ( isnumeric(lineData{i}) )
        dummyData{i} = NaN(1, size(dummyData{i}, 2));
    end
end

% Plot the dummy data, which includes both line and point style:
hDummy = plot(dummyData{:});
if ( ~isempty(lineStyle) || ~isempty(pointStyle) )
    set(hDummy, lineStyle{:}, pointStyle{:});
end

% Reset colors prior to legend data plot if running R2014b.
if ( ~verLessThan('matlab', '8.4') )
    set(gca, 'ColorOrderIndex', newColorOrder);
end

%% Setting xLim and yLim:

% If holding, then make sure not to shrink the X-limits.
if ( holdState )
    xLim = [min(xLimCurrent(1), xLim(1)), max(xLimCurrent(2), xLim(2))];
    yLim = [min(yLimCurrent(1), yLim(1)), max(yLimCurrent(2), yLim(2))];
end

% We always want to set the x-limits. Otherwise, plots like
%   plot(chebfun(@(x) sin(x), [0 pi])
% will have extra white space around the ends of the domain, and look ugly.

if ( diff(xLim) )
    set(gca, 'xLim', xLim)
end

% Set the Y-limits if appropriate values have been suggested, or if we were
% holding on when we entered this method:
if ( ~defaultYLim || (holdState && strcmp(yLimModeCurrent, 'manual')) ) && ...
        ( diff(yLim) )
    set(gca, 'yLim', yLim)
end

%% Misc:

% Return hold state to what it was before:
if ( ~holdState )
    hold off
end

% We don't want these guys to be included in LEGEND(), so we turn them off.
set(h1, 'handlevis', 'off');
set(h2, 'handlevis', 'off');
set(h3, 'handlevis', 'off');
set(h4, 'handlevis', 'off');

% Ensure that the both the handle for the dummy plot and the dummy plot itself
% are visible for LEGEND().  The dummy plot could in theory have its 'Visible'
% property set to 'off', but the data is NaN, so it doesn't draw anything
% anyway.  Moreover, if 'Visible' is 'off', the entries in the legend will be
% grayed out in R2016b and later.
set(hDummy, 'handlevis', 'on', 'visible', 'on');     %  ¯\_(o.O)_/¯

% Give an output to the plot handles if requested:
if ( nargout > 0 )
    varargout = {h1 ; h2 ; h3 ; h4 ; hDummy};
end

end

function h = plotDeltas(deltaData)
%PLOTDELTAS   Plots delta functions.
    h = [];

    % Get and save the current ColorOrder if running on R2014a or earlier.
    if ( verLessThan('matlab', '8.4') )
        originalColorOrder = get(gca, 'ColorOrder');
        colorOrder = circshift(originalColorOrder, 1);
    end

    for k = 1:1:numel(deltaData)
        % Set color for the next delta function plot.
        if ( verLessThan('matlab', '8.4') )
            % Manually manipulate the ColorOrder for R2014a or earlier.
            colorOrder = circshift(colorOrder, -1);
            set(gca, 'ColorOrder', colorOrder);
        else
            % Use ColorOrderIndex for R2014b and later.
            set(gca, 'ColorOrderIndex', k);
        end

        h = [h, plotDeltasForOneFunction(deltaData{k}{:})];
    end

    % Restore the ColorOrder if running on R2014a or earlier.
    if ( verLessThan('matlab', '8.4') )
        set(gca, 'ColorOrder', originalColorOrder);
    end
end

function h = plotDeltasForOneFunction(varargin)
%PLOTDELTASFORONEFUNCTION   Plots delta functions for one function.

h = [];
j = 1;
while ( ~isempty(varargin) )
    % Extract data which is in triplets:
    xData = varargin{1};  % locations of delta functions
    yData = varargin{2};  % magnitude of delta functions
    yBase = varargin{3};  % starting height for delta functions

    % Delete these arguments, since they have been copied:
    varargin(1:3) = [];

    % Handle the delta style argument:
    style = '';
    if ( ~isempty(varargin{1}) && iscell(varargin{1}) )
        style = varargin{1}{1};
    end
    varargin(1) = [];

    % Ignore complete NaN data:
    if ( all(isnan(xData)) )
        continue
    end

    if ( isnumeric(xData) )
        % Remove mixed NaN data:
        nanIdx = isnan(xData);
        xData(nanIdx) = [];
        yData(nanIdx) = [];
        yBase(nanIdx) = [];

        % merge duplicate delta functions.
        [yData, xData, dupIdx] = deltafun.mergeColumns(yData.', xData.');
        xData = xData.';
        yData = yData.';

        % For delta functions at discontinuities, set the base value to the
        % average of the function values:
        yBase(dupIdx-1) = 1/2*(yBase(dupIdx-1)+yBase(dupIdx));
        yBase(dupIdx) = [];
        yFinish = yBase + yData;

        % Add NaNs in order to plot isolated vertical lines:

        % NaNs for xData:
        xFullData = zeros(3*size(xData,1), 1);
        xFullData(1:3:end) = xData;
        xFullData(2:3:end) = xData;
        xFullData(3:3:end) = NaN;

        % NaNs for yData
        yFullData = zeros(3*size(xData,1), 1);
        yFullData(1:3:end) = yBase;
        yFullData(2:3:end) = yFinish;
        yFullData(3:3:end) = NaN;

        fullData{1} = xFullData;
        fullData{2} = yFullData;
        fullData{3} = style;
    end

    % Plot the vertical lines:
    h(j) = plot(fullData{:});

    % Plot the markers for delta functions:
    for jj = 1:length(xData)
        % For positive delta functions, we plot an arrow:
        if ( yData(jj) >= 0 )
            marker = '^';
        end
        % For negative delta functions, we plot a v:
        if ( yData(jj) < 0 )
            marker = 'v';
        end

        % Plot the markers and match the color with the delta lines:
        hjj = plot(xData(jj), yFinish(jj), marker, ...
            'markersize', 6, 'linestyle', 'none', 'handlevis', 'off');
        set(hjj, 'color', get(h(j), 'color'));
        set(hjj, 'markerfacecolor', get(h(j), 'color'))
    end

    j = j + 1;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% Copyright 2018 Marc René Schädler
%
% This file is part of the mobile hearing aid prototype project
% The the mobile hearing aid prototype project 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.
%
% The mobile hearing aid prototype project 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 the mobile hearing aid prototype project. If not, see http://www.gnu.org/licenses/.

function [feedback, range] = estimate_feedback(signal, fs, method)

if nargin < 3 || isempty(method)
  method = 'median';
end

signal = signal(1000:end-1000,:);
num_samples = size(signal,1);
max_samples = round(20.*fs./1000); % ~20 ms
max_range = round(3.*fs./1000); % ~3 ms
min_freq = fs./max_samples;

step_samples = round(0.1.*fs);
analysis_samples = round(0.5.*fs);
num_frames = floor((num_samples-analysis_samples)/step_samples);
fft_samples = 2.^nextpow2(2.*analysis_samples);

% Use only single precision
signal = single(signal);

% Function to calculate the ifft of a reduced spectrum which belongs to a real function
real_ifft = @(x) real(ifft([x;conj(x(end-1:-1:2))]));

if strcmp(method,'deconv')
  % Unsegmented reference implementation
  signal_fft = fft(signal);
  signal_fft(1:1+round(min_freq./fs.*size(signal,1)),3:4) = 0;
  signal_fft(1+end-round(min_freq./fs.*size(signal,1)):end,3:4) = 0;
  feedback = real(ifft(signal_fft(:,[3,4])./signal_fft(:,[1,2])));
  feedback_short = feedback_ref(1:max_samples,:);
else
  % Calculate feedback paths for each segment
  feedback_left_fft = zeros(num_frames,fft_samples./2+1,'single');
  feedback_right_fft = zeros(num_frames,fft_samples./2+1,'single');
  for i=1:num_frames
    start = 1+(i-1).*step_samples;
    stop = analysis_samples+(i-1).*step_samples;
    signal_fft = fft(signal(start:stop,:),fft_samples);
    feedback_left_fft(i,:) = signal_fft(1:end/2+1,3)./signal_fft(1:end/2+1,1);
    feedback_right_fft(i,:) = signal_fft(1:end/2+1,4)./signal_fft(1:end/2+1,2);
  end
end

switch(method)
  case 'median'
    % Function to calculate median of real and imaginary part independently
    cmedian = @(x) median(real(x))+1i.*median(imag(x));
    % Remove outliers and frequencies which dont "fit" into "max_range"
    feedback_left_fft_median = cmedian(feedback_left_fft);
    feedback_right_fft_median = cmedian(feedback_right_fft);
    feedback_left_fft_median(1:1+round(min_freq./fs.*fft_samples)) = 0;
    feedback_right_fft_median(1:1+round(min_freq./fs.*fft_samples)) = 0;
    feedback_left_median = real_ifft(feedback_left_fft_median(:));
    feedback_right_median = real_ifft(feedback_right_fft_median(:));
    feedback_left_median_short = feedback_left_median(1:max_samples);
    feedback_right_median_short = feedback_right_median(1:max_samples);
    feedback1 = single(feedback_left_median_short);
    feedback2 = single(feedback_right_median_short);
  case 'mean'
    % A variant using the mean instead of the median values
    feedback_left_fft_mean = mean(feedback_left_fft);
    feedback_right_fft_mean = mean(feedback_right_fft);
    feedback_left_fft_mean(1:1+round(min_freq./fs.*fft_samples)) = 0;
    feedback_right_fft_mean(1:1+round(min_freq./fs.*fft_samples)) = 0;
    feedback_left_mean = real_ifft(feedback_left_fft_mean(:));
    feedback_right_mean = real_ifft(feedback_right_fft_mean(:));
    feedback_left_mean_short = feedback_left_mean(1:max_samples);
    feedback_right_mean_short = feedback_right_mean(1:max_samples);
    feedback1 = single(feedback_left_mean_short);
    feedback2 = single(feedback_right_mean_short);
end

% Find level for first path
level1 = 0;
while diff(selectmain(feedback1,level1)) <= max_range
  level1 = level1 - 1;
end
level1 = level1 + 1;

% Find level for second path
level2 = 0;
while diff(selectmain(feedback2,level2)) <= max_range
  level2 = level2 - 1;
end
level2 = level2 + 1;

% Get corrsponding ranges ranges
range1 = selectmain(feedback1,level1)-1;
range2 = selectmain(feedback2,level2)-1;

feedback = [feedback1,feedback2];
range = [range1;range2];
end

function range = selectmain(in,level)
  % Dynamically select the window that covers the most energy

  min_level = 10.^(level/20);
  [~, maxidx] =  max(abs(in));
  start = find(abs(in)>min_level,1,'first');
  stop = find(abs(in)>min_level,1,'last');
  if isempty(start)
    start = maxidx;
  end
  if isempty(stop)
    stop = maxidx;
  end
  range = [start stop];
end

author: m-r-s

repo: hearingaid-prototype


text:
  function fit = de_ftab_fit_exp3(xrs, mac, varargin)
%|function fit = de_ftab_fit_exp3(xrs, mac, [options])
%|
%| Fit to DE table fm(), suitable for subsequent interpolation / extrapolation.
%| Uses a special experimental exponential model:
% todo -log(sum_k p_k exp(-m_k . s))
% fit a 3-term exponential model that has the proper derivatives
% at 0, and at +/- infinity
%|
%| in
%|	xrs	strum		X-ray spectra; see xray_read_spectra.m
%|	mac	[Ne,L]		mass attenuation coefficients
%|
%| option
%|	'show'	1|0		plot?
%|
%| out
%|	fit	strum		strum object for fitted f_m(s_1,...,s_L)
%|
%|	methods:
%|	.fmfun(sll)		fm function evaluation, for stacked array sll
%|	.fgrad(sll)		fm gradient evaluation, for stacked array sll
%|	.show_sp(en, sp)	plot true spectrum vs fitted spectrum
%|	.show_fm(sl, fm)	mesh plot of fm and its fit
%|	.show_err(sl, fm)	mesh plot of fit error
%|	.mac_eff		effective mass atten coef based on fit
%|				(valid for 'exp' only)
%|
%| Copyright 2008-8-10, Jeff Fessler, University of Michigan

%if nargin == 1 && streq(xrs, 'test'), de_ftab_fit_exp3_test, return, end
if nargin < 2, ir_usage, end

arg.show = false;
arg = vararg_pair(arg, varargin);

LL = ncol(mac);
if LL ~= 1, fail 'only 1 component done', end

MM = ncol(xrs.Ide);

fit.MM = MM;
fit.type = 'exp3';

massbar = (xrs.Ide' * mac) ./ sum(xrs.Ide)'; % [M,1]

%if ~isempty(wt), warn 'weighting ignored', end
% if isempty(wt), wt = num2cell(ones(MM,1)); end

fit.kev = cell(1,MM);
fit.mac = cell(1,MM);
fit.coef = cell(1,MM);
for mm=1:MM
	good = xrs.sp(:,mm) > 0;
	sp_m = xrs.sp(good,mm); % nonzero spectrum
	sp_m = sp_m / sum(sp_m);
	[mac_m ii] = sort(mac(good));
	if any(mac_m(2:end)) == mac_m(1), error 'non unique min', end
	if any(mac_m(1:end-1)) == mac_m(end), error 'non unique max', end
	kev_m = xrs.en(good);
	kev_m = kev_m(ii);
	q1 = sp_m(1);
	qN = sp_m(end);
	b1 = mac_m(1);
	bN = mac_m(end);
	mac_mid = (massbar(mm) - q1 * b1 - qN * bN) / (1 - q1 - qN);
	fit.mac{mm} = [b1 mac_mid bN]';
	fit.coef{mm} = [q1 1-q1-qN qN]';

	kev_mid = interp1(mac_m, kev_m, mac_mid);
	fit.kev{mm} = [kev_m(1) kev_mid kev_m(end)]';
end

meth = {'fmfun', @de_ftab_fit_exp_eval, '(sll)'; ...
	'fgrad', @de_ftab_fit_exp_grad, '(sll)'; ...
	'show_err', @de_ftab_fit_show_err, '(sl, fm)'; ...
	'show_fm', @de_ftab_fit_show_fm, '(sl, fm)'; ...
	'show_sp', @de_ftab_fit_show_sp, '(en, sp)'; ...
	};
fit = strum(fit, meth);

if arg.show
	fit.show_fm(sl, fm);
end

end % de_ftab_fit_exp3()

%
% de_ftab_fit_exp()
% todo: more documentation
%
function [fit, ffun, fgrad] = de_ftab_fit_exp(sl, fm, MM, wt, mtype, mac, kev)

sll = ndgrid_jf('mat', sl{:});

if isempty(mac)
	if isempty(mtype), error 'mac or mtype required', end
	mac = xray_read_atten(mtype, kev);
end

LL = length(sl);

if isempty(wt), wt = num2cell(ones(MM,1)); end

Ab = 1;
for ll=1:LL
	sl = col(stackpick(sll, ll));
	Ab = Ab .* exp(-sl * mac(:,ll)'); % [#s*, #E]
end

fit.kev = cell(1,MM);
fit.mac = cell(1,MM);
fit.coef = cell(1,MM);
for mm=1:MM
	if LL == 1
		dat = fm(:,mm);
	elseif LL == 2
		dat = fm(:,:,mm);
	else
		fail 'not done'
	end
	dat = stackpick(fm,mm);
	y = exp(-dat);
	Wh = spdiag(sqrt(wt{mm}(:)), 'nowarn');
	x = wls_simplex(Ab, y(:), Wh); % coefficients for each energy

	ie = x > 1e-6; % find key energies
	fit.kev{mm} = kev(ie);
	fit.mac{mm} = mac(ie,:); % [#E,L]
	A = 1;
	for ll=1:LL
		sl = col(stackpick(sll,ll));
		A = A .* exp(-sl * fit.mac{mm}(:,ll)'); % [#s*, #E]
	end
	fit.coef{mm} = wls_simplex(A, y(:), Wh); % refit with key energies
end
ffun = @de_ftab_fit_exp_eval;
fgrad = @de_ftab_fit_exp_grad;

end % de_ftab_fit_exp()

%
% de_ftab_fit_exp_eval()
% evaluate
% in
%	sll	[(Nd),L]	stackup of s1,s2,...,s_L
% out
%	f	[(Nd),M]	stackup of f1,f2,...,f_M
%
function f = de_ftab_fit_exp_eval(fit, sll)
Nd = size(sll); LL = Nd(end); Nd = Nd(1:end-1);
sll = reshape(sll, [], LL); % [*Nd,L]
MM = fit.MM;
f = zeros(prod(Nd),MM);
for mm=1:MM
	A = 1;
	mac = fit.mac{mm};
	for ll=1:LL
		sl = sll(:,ll);
		A = A .* exp(-sl * mac(:,ll)'); % [*Nd,ne]
	end
	tmp = -log(A * fit.coef{mm}); % [*Nd,1]
	f(:,mm) = tmp;
end
f = reshape(f, [Nd MM]);

end % de_ftab_fit_exp_eval()

%
% de_ftab_fit_exp_grad()
% evaluate gradient of f for each of the given s vectors.
% in
%	sll	[(Nd),L]	stackup of s1,s2,...,s_L
% out
%	g	[(Nd),L,M]	stackup of gradients of f(s)
%
function g = de_ftab_fit_exp_grad(fit, sll)
Nd = size(sll); LL = Nd(end); Nd = Nd(1:end-1);
sll = reshape(sll, [], LL); % [*Nd,L]
MM = fit.MM;
g = zeros(prod(Nd), LL, MM);
for mm=1:fit.MM
	A = 1;
	mac = fit.mac{mm}; % [ne,L]
	alf = fit.coef{mm}; % [ne,1]
	for ll=1:LL
		sl = sll(:,ll);
		A = A .* exp(-sl * mac(:,ll)'); % [*Nd,ne]
	end
	vm = A * alf; % [*Nd,1]
	tmp = A * (mac .* repmat(alf, [1 LL])); % [*Nd,L]
	g(:,:,mm) = tmp ./ repmat(vm, [1 LL]);
end
g = reshape(g, [Nd LL MM]);

end % de_ftab_fit_exp_grad()

%
% de_ftab_fit_show_sp()
% compare true spectra to fitted spectra
%
function out = de_ftab_fit_show_sp(fit, en, sp)
if nargin < 3, fail 'de_ftab_fit_show_sp(fit, en, sp)', end

if ~streq(fit.type, 'exp3'), printm 'show_sp only done for exp3', return, end
if im
	clf, pl = (fit.MM+1)*100 + 10 + 1;
	subplot(pl)
	plot(en, sp * diag(1 ./ max(sp)))
	if isfield(fit, 'kev')
		for mm=1:fit.MM
			subplot(pl+mm)
			bar(fit.kev{mm}, fit.coef{mm})
			axis tight, axisx(minmax(en))
		end
	else
		warn 'kev unknown'
	end
end

if nargout, out = []; end

end % de_ftab_fit_show_sp()

%
% de_ftab_fit_show_fm()
% compare fit to sampled fm
%
function out = de_ftab_fit_show_fm(fit, sl, fm)
if nargin < 3, error 'de_ftab_fit_show_fm(fit, sl, fm)', end

sll = ndgrid_jf('mat', sl{:});

if fit.MM ~= 2, error 'only M=2 done', end

fh = fit.fmfun(sll);

switch length(sl) % LL
case 1
	s1 = sl{1};
	im clf
	plot(	s1, fm(:,1), 'c.', s1, fm(:,2), 'y.', ...
		s1, fh(:,1), 'c-', s1, fh(:,2), 'y-')
	xlabel 's1'
	ylabel 'fm(s1)'
	legend('true m=1', 'true m=2', 'fit m=1', 'fit m=2', 4)

case 2
	s1 = sl{1};
	s2 = sl{2};
	smax(1) = max(s1);
	smax(2) = max(s2);
	fmax = max(fm(:));

	ax = [0 smax(1) 0 smax(2) 0 fmax];
	cax = [0 fmax];

	im clf
	im pl 2 2

	show(1, s1, s2, fm(:,:,1), ax, cax, 'f_1(s)')
	% text(-60, 10, '[cm^2/g]')
	show(2, s1, s2, fm(:,:,2), ax, cax, 'f_2(s)')

	show(3, s1, s2, fh(:,:,1), ax, cax, 'f_1 approx')
	show(4, s1, s2, fh(:,:,2), ax, cax, 'f_2 approx')
otherwise
	fail 'not done'
end

if nargout, out = []; end

end % de_ftab_fit_show_fm()

%
% de_ftab_fit_show_err()
% show fit errors, and relate to HU
% f = mac s, mac(70kev,H2O) = 0.2 cm^2 / g = 1000 HU
% mh = mac * 1000 HU / (0.2 g/cm^2)
% so Df/Dmh = Df/Dmac * Dmac/Dmh = (50 g/cm^2) (0.2 cm^2/g) / 1000 HU = 1/100 HU
% so Dmh = 100 HU * Df for 50 cm of H2O.
%
function out = de_ftab_fit_show_err(fit, sl, fm)
if nargin < 3, error 'de_ftab_fit_show_err(fit, sl, fm)', end

sll = ndgrid_jf('mat', sl{:});

fh = fit.fmfun(sll);
err = fh - fm;
printm('worst model error = %g of %g', max(abs(err(:))), max(fm(:)))
printm('worst model error = %g HU over 50cm H2O', 100*max(abs(err(:))))
for mm=1:fit.MM
	ee = stackpick(err,mm);
	printm('worst error (mm=%d): %g', mm, max(abs(ee(:))))
end

if max(abs(err(:))) == 0, return, end
if fit.MM > 2, error 'only M=2 done', end
if 0
	e1 = err(:,:,1);
	e2 = err(:,:,2);
	disp([minmax(e1); minmax(e2)]')
	printm('worst error1 %g', max(col(abs(err(:,1,:)))))
	printm('worst error2 %g', max(col(abs(err(:,2,:)))))
end
%err = abs(err);

switch length(sl) % LL
case 1
	s1 = sl{1};
	plot(s1, err(:,1), 'c-', s1, err(:,2), 'y-')
	xlabel 's1'
	ylabel 'error'
	legend('m=1', 'm=2', 4)

case 2
	s1 = sl{1};
	s2 = sl{2};
	smax(1) = max(s1);
	smax(2) = max(s2);

	elim = minmax(err(:))';
	%elim = [-1 1] * 0.01; % +/- 1 HU
	ax = [0 smax(1) 0 smax(2) elim];
	im plc 1 2
	for mm=1:fit.MM
		show(mm, s1, s2, err(:,:,mm), ax, elim, sprintf('f_%d error', mm))
	end
otherwise
	fail 'not done'
end

if nargout, out = []; end

end % de_ftab_fit_show_err()

%
% show()
%
function show(pl, x, y, f, ax, cax, ti)
if ~im, return, end
im('subplot', pl)
if 1
	mesh(x,y,f')
	colormap hsv, caxis(cax), cbar
	axis(ax)
	xtick, ytick, ztick, zwhite
else
	im(x,y,f), cbar
	xtick, ytick
end
xlabel 's_1', ylabel 's_2', title(ti)

end % show()

%
% de_ftab_fit_exp3_test()
%
function de_ftab_fit_exp3_test

%stype = 'mono,70';
stype = 'ps1';
xrs = xray_read_spectra(stype);

if 1 % L=1 component
	sl{1} = linspace(0, 50, 26);
	mtype = {'water'};
	mas = xray_read_mac(mtype);
	mac = mas.mac(xrs.en);
	if im
		clf, semilogy(xrs.en, mac), legend(mtype{:})
	end
	sll = ndgrid_jf('mat', sl{:});
	fm = de_ftab_fm(sll, mac, xrs.Ide);
	fit = de_ftab_fit(sl, fm, 'type', 'exp', 'mtype', mtype);
	g = fit.fgrad(sll);
%	fit = de_ftab_fit(sl, fm, 'type', 'poly')

	if 1
		fit.show_sp(xrs.en, xrs.sp);
		prompt
		fit.show_fm(sl, fm);
		prompt
		fit.show_err(sl, fm);
		prompt
	end
end

if 1 % L=2
	sl{1} = linspace(0, 50, 26);
	sl{2} = linspace(0, 30, 31);
	mtype = {'water', 'bone'};
	mas = xray_read_mac(mtype);
	mac = mas.mac(xrs.en);
	if im
		clf, semilogy(xrs.en, mac), legend(mtype{:})
	end
	sll = ndgrid_jf('mat', sl{:});
	fm = de_ftab_fm(sll, mac, xrs.Ide);
	fit = de_ftab_fit(sl, fm, 'type', 'exp', 'mtype', mtype);
	g = fit.fgrad(sll);
	%fit = de_ftab_fit(sl, fm, 'type', 'poly')

	if 1
		fit.show_sp(xrs.en, xrs.sp);
		prompt
		fit.show_fm(sl, fm);
		prompt
		fit.show_err(sl, fm);
		prompt
	end
end

end % de_ftab_fit_exp3_test()

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function criterion = gbc_variants(ori,CS,Dl,Dr,delta,varargin)
errorcheck = all(length(ori)==length(delta) & all(floor(delta) == delta));
assert(errorcheck,'Provide a list of valid variant IDs to compute grains from variant Ids');
% variant grain criterion
% check whether the two ebsd measurements are seperated by a grain boundary
if size(delta,2)>1
    criterion = delta(Dl,1)-delta(Dr,1) + delta(Dl,2)-delta(Dr,2) == 0;
else
    criterion = delta(Dl,1)-delta(Dr,1) == 0;
end

author: mtex-toolbox

repo: mtex

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


text:
function G = myLinearize(BlockData)
%Copyright 2013 The MathWorks, Inc.

% assignin('base','BlockData',BlockData)
[x,u] = findop(BlockData.Parameters.Value,'snapshot',2,ones(1000,1));
G = linearize(BlockData.Parameters.Value,u,x);

end

author: Sable

repo: mcbench-benchmarks


text:
function s=p2s(p)
% function s=p2s(p)
%
% DESCRIPTION
%   Converts from a multipoly polynomial to a symbolic toolbox polynomial.
%
% INPUTS
%   p: Polynomial created using the multipoly toolbox
%
% OUTPUTS
%   s: Polynomial created using the symbolic toolbox
%
% SYNTAX
%   s = p2s(p)
%
% See also s2p

% 1/30/2003  PJS  Initial Coding
% 11/23/2010 PJS  Updated for matrix polynomials

% Get polynomial info
[nr nc] = size(p);
[nt,nv] = size(p.degmat);
coef = p.coefficient;
deg = p.degmat;
var = p.varname;

% Create symbolic toolbox variables
for i1 = 1:nv
    eval(['syms ' var{i1} ' real']);
end

% Construct monomial by monomial
s = zeros(nr,nc);
for i1 = 1:nt
    termexp = deg(i1,:);
    idx = find(termexp);
    term = 1;
    for i2 = 1:length(idx)
        eval(['term = term*(' var{idx(i2)} ')^termexp(idx(i2));']);
    end
    s = s + reshape(coef(i1,:),[nr,nc])*term;
end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function pass = test_std3(pref)
% Test the chebfun3/max command.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e3*pref.cheb3Prefs.chebfun3eps;

f = chebfun3(10); % (x,y,z)-std of any fixed value is zero.
exact = 0;
h = std3(f);
pass(1) = norm(h - exact) < tol;

end


text:
function varargout = guiDatePicker(varargin)
% GUIDATEPICKER MATLAB code for guiDatePicker.fig
%      GUIDATEPICKER, by itself, creates a new GUIDATEPICKER or raises the existing
%      singleton*.
%
%      H = GUIDATEPICKER returns the handle to a new GUIDATEPICKER or the handle to
%      the existing singleton*.
%
%      GUIDATEPICKER('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in GUIDATEPICKER.M with the given input arguments.
%
%      GUIDATEPICKER('Property','Value',...) creates a new GUIDATEPICKER or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before guiDatePicker_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to guiDatePicker_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help guiDatePicker

% Last Modified by GUIDE v2.5 21-Apr-2011 12:54:26

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @guiDatePicker_OpeningFcn, ...
                   'gui_OutputFcn',  @guiDatePicker_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before guiDatePicker is made visible.
function guiDatePicker_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to guiDatePicker (see VARARGIN)
year = 2000:1:2020;
dt = varargin{1}; %
if ~isempty(dt)
    [YYi, MMi, DDi, hhI, mmI, ssI] = datevec(dt);
    wo = find(year==YYi);
    set(handles.Year,'Value',wo);
    set(handles.Month,'Value',MMi);
    set(handles.Hours,'String',num2str(hhI));
    set(handles.Minutes,'String',num2str(mmI));
end

y = get(handles.Year,'String');
yi = get(handles.Year,'Value');
m = get(handles.Month,'String');
mi = get(handles.Month,'Value');
yy = y(yi);

YY = str2double(yy);
MM = mi;
ldom = eomday(YY,MM);
for i = 1:ldom
    dtn(i) = datenum([YY MM i 0 0 0]);
    dts{i} = datestr(dtn);
    %disp(char(dts(i)));
end
    [n,s] = weekday(datenum([YY MM 1 0 0 0]));
% clear calendar
for i = 1:7
    for j = 1:6
        D{j,i} = '';
    end
end
% update calendar
num = 0;
for i = n:7
    num = num + 1;
    D{1,i} = num;
end
for j = 2:6
    for i = 1:7
        num = num + 1;
        if num > ldom
            break;
        end
        D{j,i} = num;
    end
end

set(handles.Calendar,'Data',D);

% Choose default command line output for guiDatePicker
handles.output = hObject;

% Close handles structure
guidata(hObject, handles);

% UIWAIT makes guiDatePicker wait for user response (see UIRESUME)
uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = guiDatePicker_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;
close(gcf);

% --- Executes on selection change in Year.
function Year_Callback(hObject, eventdata, handles)
% hObject    handle to Year (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns Year contents as cell array
%        contents{get(hObject,'Value')} returns selected item from Year

y = get(handles.Year,'String');
yi = get(handles.Year,'Value');
m = get(handles.Month,'String');
mi = get(handles.Month,'Value');
yy = y(yi);

YY = str2double(yy);
MM = mi;
ldom = eomday(YY,MM);
for i = 1:ldom
    dtn(i) = datenum([YY MM i 0 0 0]);
    dts{i} = datestr(dtn);
    %disp(char(dts(i)));
end
    [n,s] = weekday(datenum([YY MM 1 0 0 0]));
% clear calendar
for i = 1:7
    for j = 1:6
        D{j,i} = '';
    end
end
% update calendar
num = 0;
for i = n:7
    num = num + 1;
    D{1,i} = num;
end
for j = 2:6
    for i = 1:7
        num = num + 1;
        if num > ldom
            break;
        end
        D{j,i} = num;
    end
end

set(handles.Calendar,'Data',D);

% --- Executes during object creation, after setting all properties.
function Year_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Year (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in Month.
function Month_Callback(hObject, eventdata, handles)
% hObject    handle to Month (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns Month contents as cell array
%        contents{get(hObject,'Value')} returns selected item from Month
y = get(handles.Year,'String');
yi = get(handles.Year,'Value');
m = get(handles.Month,'String');
mi = get(handles.Month,'Value');
yy = y(yi);

YY = str2double(yy);
MM = mi;
ldom = eomday(YY,MM);
for i = 1:ldom
    dtn(i) = datenum([YY MM i 0 0 0]);
    dts{i} = datestr(dtn);
    %disp(char(dts(i)));
end
    [n,s] = weekday(datenum([YY MM 1 0 0 0]));
% clear calendar
for i = 1:7
    for j = 1:6
        D{j,i} = '';
    end
end
% update calendar
num = 0;
for i = n:7
    num = num + 1;
    D{1,i} = num;
end
for j = 2:6
    for i = 1:7
        num = num + 1;
        if num > ldom
            break;
        end
        D{j,i} = num;
    end
end

set(handles.Calendar,'Data',D);

% --- Executes during object creation, after setting all properties.
function Month_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Month (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in Close.
function Close_Callback(hObject, eventdata, handles)
% hObject    handle to Close (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

y = get(handles.Year,'String');
yi = get(handles.Year,'Value');
m = get(handles.Month,'String');
mi = get(handles.Month,'Value');
yy = y(yi);

YY = str2double(yy);
MM = mi;
if isfield(handles,'Day')
    DD = handles.Day;
    hh = str2double(get(handles.Hours,'String'));
    mm = str2double(get(handles.Minutes,'String'));

    % x = datenum([YY MM DD hh mm 0])
    handles.output = datenum([YY MM DD hh mm 0]);
    % Close handles structure
    guidata(hObject, handles);

    uiresume;
else
   errordlg('Day was not selected','DatePicker Error');
end

% --- Executes during object creation, after setting all properties.
function Calendar_CreateFcn(hObject, eventdata, handles)
% hObject    handle to CTDtab (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% --- Executes when selected cell(s) is changed in Calendar.
function Calendar_CellSelectionCallback(hObject, eventdata, handles)
% hObject    handle to Calendar (see GCBO)
% eventdata  structure with the following fields (see UITABLE)
%	Indices: row and column indices of the cell(s) currently selecteds
% handles    structure with handles and user data (see GUIDATA)
        sel = eventdata.Indices;      % Get selection indices (row, col)
                                     % Noncontiguous selections are ok
        selcols = unique(sel(:,2));   % Get all selected data col IDs
        selrows = unique(sel(:,1));
        num = get(hObject,'Data');  % Get copy of uitable data

x = num(selrows,selcols);
handles.Day = x{1};

% Close handles structure
guidata(hObject, handles);

function Hours_Callback(hObject, eventdata, handles)
% hObject    handle to Hours (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of Hours as text
%        str2double(get(hObject,'String')) returns contents of Hours as a double

% --- Executes during object creation, after setting all properties.
function Hours_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Hours (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function Minutes_Callback(hObject, eventdata, handles)
% hObject    handle to Minutes (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of Minutes as text
%        str2double(get(hObject,'String')) returns contents of Minutes as a double

% --- Executes during object creation, after setting all properties.
function Minutes_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Minutes (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


text:
function [njm] = tapas_ti_zeromat(jm)
%% Zeros all but the first entry in each column of jm.
%
% Input
%   jm      -- Matrix.
%
% Output
%   njm     -- Matrix.
%

% aponteeduardo@gmail.com
% copyright (C) 2015
%

njm = zeros(size(jm));

j = 0;
i = 0;

while j < size(jm, 2)
    j = j + 1;
    i = 0;
    while i < size(jm, 1)
        i = i + 1;
        if jm(i, j);
            njm(i, j) = 1;
            break
        end
    end
end

end

author: translationalneuromodeling

repo: tapas


text:
function [confusion_matrix, classes]=cosmo_confusion_matrix(ds, varargin)
% Returns a confusion matrix
%
% Usage 1: mx=cosmo_confusion_matrix(ds)
% Usage 2: mx=cosmo_confusion_matrix(targets, predicted)
%
%
% Inputs:
%   targets     Nx1 targets for N samples, or a dataset struct with
%               .sa.targets
%   predicted   NxM predicted labels (from a classifier), for N samples and
%               M predictions per set of samples
%
% Returns:
%   mx          PxPxM matrix assuming there are P unique targets.
%               mx(i,j,k)==c means that the i-th target class was classified
%               as the j-th target class c times for the k-th set of
%               samples.
%   classes     Px1 class labels.
%
% Example:
%     ds=cosmo_synthetic_dataset('ntargets',3,'nchunks',4);
%     args=struct();
%     args.partitions=cosmo_nchoosek_partitioner(ds,1);
%     args.output='winner_predictions';
%     args.classifier=@cosmo_classify_lda;
%     pred_ds=cosmo_crossvalidation_measure(ds,args);
%     confusion=cosmo_confusion_matrix(pred_ds.sa.targets,pred_ds.samples);
%     cosmo_disp(confusion)
%     %|| [ 3         0         1
%     %||   0         3         1
%     %||   1         0         3 ]
%     confusion_alt=cosmo_confusion_matrix(pred_ds);
%     isequal(confusion,confusion_alt)
%     %|| true
%     %
%     % run a searchlight with tiny radius of 1 voxel (3 is more common)
%     nbrhood=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);
%     measure=@cosmo_crossvalidation_measure;
%     sl_ds=cosmo_searchlight(ds,nbrhood,measure,args,'progress',false);
%     %
%     % the confusion matrix is 3x3x6, that is 6 3x3 confusion
%     % matrices. Here the dataset is passed directly
%     sl_confusion=cosmo_confusion_matrix(sl_ds);
%     cosmo_disp(sl_confusion)
%     %|| <double>@3x3x6
%     %||    (:,:,1) = [ 4         0         0
%     %||                0         4         0
%     %||                0         1         3 ]
%     %||    (:,:,2) = [ 4         0         0
%     %||                0         4         0
%     %||                0         1         3 ]
%     %||    (:,:,3) = [ 2         1         1
%     %||                0         4         0
%     %||                1         0         3 ]
%     %||    (:,:,4) = [ 4         0         0
%     %||                0         3         1
%     %||                0         1         3 ]
%     %||    (:,:,5) = [ 3         0         1
%     %||                0         4         0
%     %||                1         1         2 ]
%     %||    (:,:,6) = [ 3         0         1
%     %||                0         4         0
%     %||                1         1         2 ]
%
%     % using samples that are not predictions gives an error
%     ds=cosmo_synthetic_dataset('ntargets',3,'nchunks',4);
%     confusion=cosmo_confusion_matrix(ds)
%     %|| error('72 predictions mismatch targets, first is (1,1)=2.211999e+00')
%
% Notes:
%   - this function counts the number of times each sample was classified
%     as any target
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    [targets,predicted]=get_data(ds,varargin{:});

    % see which classes there are
    [class_indices,classes]=cosmo_index_unique(targets);
    nclasses=numel(class_indices);

    % allocate space for output
    nfeatures=size(predicted,2);
    confusion_matrix=zeros([nclasses,nclasses,nfeatures]);

    % keep track which predicted samples were in targets
    visited=false(size(predicted));
    % >@@>
    for k=1:nclasses
        % rows for k-th class
        idxs=class_indices{k};
        for j=1:nclasses
            match_msk=bsxfun(@eq,classes(j),predicted(idxs,:));
            confusion_matrix(k,j,:)=sum(match_msk,1);
            visited(idxs,:)=visited(idxs,:) | match_msk;
        end
    end

    % <@@<

    missing=~(visited | isnan(predicted));
    if any(missing(:))
        n=sum(missing(:));
        [i,j]=find(missing,1);
        error(['%d predictions mismatch targets, '...
                'first is (%d,%d)=%d'],n,i,j,predicted(i,j));
    end

function [targets,predicted]=get_data(ds, predicted)
    has_predicted=nargin>=2;
    is_ds=isstruct(ds);
    if is_ds
        if has_predicted
            error('Need exactly one argument when input is struct');
        end
        % input is a dataset
        cosmo_isfield(ds,'sa.targets',true);
        cosmo_isfield(ds,'samples',true);

        predicted=ds.samples;
        targets=ds.sa.targets;
    elseif isnumeric(ds)
        if ~has_predicted
            error('Need two arguments when first argument is numeric');
        end
        targets=ds;
    else
        error('Illegal input: need struct or numeric vector');
    end

    if ~isvector(targets) || size(targets,2)~=1
        error('targets must be column vector');
    end

    if numel(size(predicted))~=2
        error('predictions must be matrix');
    end

    nsamples=numel(targets);
    if size(predicted,1)~=nsamples
        error(['Size mismatch: predictions has %d values on first '...
                'dimension, but targets has %d values'],...
                    size(predicted,1),nsamples);
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [m,c]=cep2pow(u,v,mode)
%CEP2POW convert cepstral means and variances to the power domain
% Inputs:
%    u: vector giving the cepstral means with u(1) the 0'th cepstral coefficient
%    v: cepstral covariance matrix or else a vector containing the diagonal elements
% mode: 'c'  pow=exp(irdct(cep))    [default]
%       'f'  pow=exp(rsfft(cep)/n)  [fft length even]
%       'fo' pow=exp(rsfft(cep)/n)  [fft length odd]
%       'i'  pow=exp(cep)           [ no transformation ]
%
% Outputs:
%    m: row vector giving means in the power domain
%    c: covariance matrix in the power domain

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

if nargin<3 mode='c'; end
if min(size(v))==1
   v=diag(v);
end
u=u(:)';    % force u to be a row vector
if any(mode=='f')
   n=2*length(u)-2;
   if any(mode=='o')
      n=n+1;
   end
   p=rsfft(u',n)/n;
   q=rsfft(rsfft(v,n)',n)/n^2;
elseif any(mode=='i')
    p=u';
    q=v';
else
   p=irdct(u');
   q=irdct(irdct(v)');
end
m=exp(p+0.5*diag(q))';
c=(m'*m).*(exp(q)-1);

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
%DEMOEXPANDPOLYGON  Expand a polygon by a given distance
%   demoExpandPolygon
%
%   Example
%   demoExpandPolygon
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-06-17,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

%% Base polygon

% create a non-convex polygon
poly = [200 100;200 150;150 150;150 200;50 200;150 100];

% set up display
figure(1); clf; hold on;
axis equal; axis([20 230 70 230]);

% show the polygon
drawPolygon(poly, 'linewidth', 2);

print(gcf, 'expandPolygon_initial.png', '-dpng');

%% Positive expansion

% expand the polygon by a positive distance (outside of the polygon)
polyOut = expandPolygon(poly, 10);

% draw the expanded polygon
drawPolygon(polyOut, 'k');

print(gcf, 'expandPolygon_outer.png', '-dpng');

%% Negative expansion

% expand the polygon by a negative distance (inside the polygon)
polyIn = expandPolygon(poly, -20, 'cleanuploops', true);

% draw the expanded polygon
drawPolygon(polyIn, 'g');

% decorate
legend('Original', 'Outside', 'Inside')

print(gcf, 'expandPolygon_OuterAndInner.png', '-dpng');

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function X = mtraceprod(A,B)
% X = MTRACEPROD(A, B)
%
% Evaluates trace(A*B') efficiently.
%
% X = MTRACEPROD(A)
%
% Evaluates trace(A*A') efficiently.
%
% If A and/or B has third dimension, the third dimension of A results in
% rows to X and the third dimension of B columns to X:
%
% X(I,J) = trace(A(:,:,I)*B(:,:,J)')

%warning('Deprecated.');

if nargin < 2
  [M,N,D] = size(A);
  X = reshape(A,[M*N,D])' * reshape(A,[M*N,D]);
else
  [M,N,Da] = size(A);
  [M,N,Db] = size(B);
  X = reshape(A,[M*N,Da])' * reshape(B,[M*N,Db]);
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function SO3F = rotate_outer(SO3F,rot,varargin)
% rotate a function on SO(3)
%
% Syntax
%   SO3F = SO3F.rotate_outer(rot)
%
% Input
%  SO3F - @SO3FunHarmonic
%  rot  - @rotation
%
% Output
%  SO3F - @SO3FunHarmonic
%

L = SO3F.bandwidth;
D = WignerD(rot,'bandwidth',L);

for l = 0:L
  SO3F.fhat(deg2dim(l)+1:deg2dim(l+1)) = ...
    reshape(SO3F.fhat(deg2dim(l)+1:deg2dim(l+1)),2*l+1,2*l+1) * ...
    reshape(D(deg2dim(l)+1:deg2dim(l+1)),2*l+1,2*l+1).' ;
end

end


text:
function DCM = spm_dcm_ind(DCM)
% Estimate parameters of a (bilinear) DCM of induced spectral responses
% FORMAT DCM = spm_dcm_ind(DCM)
%
% DCM
%    name: name string
%       M:  Forward model
%              M.dipfit - leadfield specification
%       xY: data   [1x1 struct]
%       xU: design [1x1 struct]
%
%   Sname: cell of source name strings
%       A: {[nr x nr double]  [nr x nr double]  [nr x nr double]}
%       B: {[nr x nr double], ...}   Connection constraints
%       C: [nr x 1 double]
%
%   options.Nmodes       - number of frequency modes
%   options.Tdcm         - [start end] time window in ms
%   options.D            - time bin decimation       (usually 1 or 2)
%   options.h            - number of DCT drift terms (usually 1 or 2)
%   options.type         - 'ECD' (1) or 'Imaging' (2) (see spm_erp_L)
%   options.onset        - stimulus onset (ms)
%   options.dur          - and dispersion (sd)
%______________________________________________________________________
% This routine inverts dynamic causal models (DCM) of induced or spectral
% responses as measured with the electroencephalogram (EEG) or the
% magnetoencephalogram (MEG). It models the time-varying power, over a
% range of frequencies, as the response of a distributed system of coupled
% electromagnetic sources to a spectral perturbation. The model parameters
% encode the frequency response to exogenous input and coupling among
% sources and different frequencies. Bayesian inversion of this model,
% given data enables inferences about the parameters of a particular model
% and allows one to compare different models, or hypotheses. One key aspect
% of the model is that it differentiates between linear and non-linear
% coupling; which correspond to within and between frequency coupling
% respectively.
%
% The number of nodes can be optimised using Bayesian model selection. The
% data are reduced to a fixed number of principal components that capture
% the greatest variation inspection responses never peristimulus time. The
% number of nodes specified by the user tries to reconstruct the response
% in the space of the principle components or eigenmodes using a reduced
% set of eigenmodes. The number of modes corresponding to data features can
% be changed (from Nf = 8) by editing spm_dcm_ind_data.m
%
% see also: spm_dcm_ind_data; spm_gen_ind; spm_fx_ind and spm_lx_ind
%
% See: Chen CC, Kiebel SJ, Friston KJ.
% Dynamic causal modelling of induced responses.
% Neuroimage. 2008 Jul 15;41(4):1293-312.
% ______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_ind.m 7143 2017-07-29 18:50:38Z karl $

% check options
%==========================================================================
clear spm_erp_L
DCM.options.analysis  = 'IND';

% Filename and options
%--------------------------------------------------------------------------
try, DCM.name;                  catch, DCM.name           = 'DCM_IND'; end
try, DCM.options.Nmodes;        catch, DCM.options.Nmodes = 4;         end
try, onset = DCM.options.onset; catch, onset              = 80;        end
try, dur   = DCM.options.dur;   catch, dur                = 32;        end
try, DATA  = DCM.options.DATA;  catch, DATA               = 1;         end

% Data and spatial model
%==========================================================================
if DATA
    DCM      = spm_dcm_erp_dipfit(DCM, 1);
    if ~isfield(DCM.xY,'source')
        DCM  = spm_dcm_ind_data(DCM);
    end
end
xY     = DCM.xY;
xU     = DCM.xU;
xU.dt  = xY.dt;

% dimensions
%--------------------------------------------------------------------------
Nm     = DCM.options.Nmodes;          % number of frequency modes modelled
Nf     = size(xY.U,2);                % number of frequency modes explained
Nt     = length(xY.y);                % number of trials
Nr     = size(DCM.C,1);               % number of sources
nu     = size(DCM.C,2);               % number of neuronal inputs
Ns     = size(xY.y{1},1);             % number of samples
Nu     = size(xU.X,2);                % number of trial-specific effects
nx     = Nr*Nf;                       % number of states

% assume noise precision is the same over modes:
%--------------------------------------------------------------------------
xY.Q   = {spm_Q(1/2,Ns,1)};
xY.X0  = sparse(Ns,0);

% Inputs
%==========================================================================

% trial-specific effects
%--------------------------------------------------------------------------
try
    if length(DCM.B) ~= Nu;
        warndlg({'please ensure number of trial specific effects', ...
                 'encoded by DCM.xU.X & DCM.B are the same'})
    end
catch
    DCM.B = {};
end

% model specification and nonlinear system identification
%==========================================================================
M      = DCM.M;
try, M = rmfield(M,'g');  end
try, M = rmfield(M,'FS'); end

% prior moments
%--------------------------------------------------------------------------
[pE,gE,pC,gC] = spm_ind_priors(DCM.A,DCM.B,DCM.C,Nm,Nf);

% hyperpriors (assuming about 99% signal to noise)
%--------------------------------------------------------------------------
hE    = 8 - log(var(spm_vec(xY.y)));
hC    = exp(-4);

% likelihood model
%--------------------------------------------------------------------------
M.IS  = 'spm_gen_ind';
M.f   = 'spm_fx_ind';
M.G   = 'spm_lx_ind';
M.x   = sparse(nx,1);
M.pE  = pE;
M.pC  = pC;
M.gE  = gE;
M.gC  = gC;
M.hE  = hE;
M.hC  = hC;
M.m   = nu;
M.n   = nx;
M.l   = Nr*Nf;
M.ns  = Ns;
M.ons = onset - xY.pst(1);
M.dur = dur;

% EM: inversion
%--------------------------------------------------------------------------
[Qp,Qg,Cp,Cg,Ce,F] = spm_nlsi_N(M,xU,xY);

% Data ID
%==========================================================================
if isfield(M,'FS') && DATA
    try
        ID  = spm_data_id(feval(M.FS,xY.y,M));
    catch
        ID  = spm_data_id(feval(M.FS,xY.y));
    end
else
    ID  = spm_data_id(xY.y);
end

% Bayesian inference {threshold = prior} NB Prior on A,B  and C = exp(0) = 1
%==========================================================================
warning('off','SPM:negativeVariance');
dp  = spm_vec(Qp) - spm_vec(pE);
Pp  = spm_unvec(1 - spm_Ncdf(0,abs(dp),diag(Cp)),Qp);
warning('on','SPM:negativeVariance');

% neuronal and sensor responses (x and y)
%--------------------------------------------------------------------------
L   = feval(M.G, Qg,M);        % get gain matrix
x   = feval(M.IS,Qp,M,xU);     % prediction (source space)

% trial-specific responses (in mode, channel and source space)
%--------------------------------------------------------------------------
for i = 1:Nt

    s  = x{i};                 % prediction (source space)
    y  = s*L';                 % prediction (sensor space)
    r  = xY.y{i} - y;          % residuals  (sensor space)

    % parse frequency modes
    %----------------------------------------------------------------------
    for j = 1:Nm
        f      = (1:Nr) + (j - 1)*Nr;
        H{i,j} = y(:,f);
        R{i,j} = r(:,f);
        K{i,j} = s(:,f);
    end
end

% store estimates in DCM
%--------------------------------------------------------------------------
DCM.M  = M;                    % model specification
DCM.xY = xY;                   % data structure
DCM.xU = xU;                   % input structure
DCM.Ep = Qp;                   % conditional expectation f(x,u,p)
DCM.Cp = Cp;                   % conditional covariances G(g)
DCM.Eg = Qg;                   % conditional expectation
DCM.Cg = Cg;                   % conditional covariances
DCM.Pp = Pp;                   % conditional probability
DCM.H  = H;                    % conditional responses (y), channel space
DCM.K  = K;                    % conditional responses (x)
DCM.R  = R;                    % conditional residuals (y), channel space
DCM.Ce = Ce;                   % ReML error covariance
DCM.F  = F;                    % Laplace log evidence
DCM.ID = ID;                   % data ID
DCM.L  = [];

% and save
%--------------------------------------------------------------------------
save(DCM.name, 'DCM', spm_get_defaults('mat.format'));
assignin('base','DCM',DCM)

author: spm


text:
function [sys,x0,str,ts]=ESO2_NL(t,x,u,flag,a1,a2,d,Beta,b0)
switch flag,
    case 0
      sys=[3,0,3,2,0,0,1];
      x0=[0;0;0];
      str=[];
      ts=[0 0];
    case 1
      e=x(1)-u(1);
      sys(1)=x(2)-Beta(1)*e;
      sys(2)=x(3)-Beta(2)*fal(e,a1,d)+b0*u(2);
      sys(3)=-Beta(3)*fal(e,a2,d);
    case 3
        sys=x;
    case {2,4,9}
        sys=[];
    otherwise
        error(['Unhandled flag=',num2str(flag)]);
end
function f=fal(e,a,d)
    if abs(e)<d
        f=e*d^(a-1);
    else f=(abs(e))^a*sign(e);
    end

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
function l = GridLength(G)
% return number of points

l = cellfun('prodofsize',{G.points});


text:
function display(v,varargin)
% standard output

if check_option(varargin,'onlyShowVectors')
  d = [v.x(:),v.y(:),v.z(:)];
  d(abs(d) < 1e-10) = 0;

  cprintf(d,'-L','  ','-Lc',{'x' 'y' 'z'});
  return
end

displayClass(v,inputname(1),varargin{:});%...
  %'moreInfo',char(v.refSystem,'compact'),varargin{:});

if length(v) ~= 1, disp([' size: ' size2str(v)]);end

if v.antipodal, disp(' antipodal: true'); end

% display resolution
if isOption(v,'resolution')
  disp([' resolution: ',xnum2str(getOption(v,'resolution')/degree),mtexdegchar]);
  v.opt = rmfield(v.opt,'resolution');
end

% display all other options
disp(char(dynOption(v)));

% display coordinates
if check_option(varargin,'skipCoordinates') || isempty(v) || getMTEXpref('generatingHelpMode')

elseif check_option(varargin,'all') || (length(v) < 20)

  display(v,'onlyShowVectors')

else

  disp(' ')
  s = setappdata(0,'data2beDisplayed',v);
  disp(['  show vectors'])
  disp(' ')

end


text:
% Test routine for mtimesx, multi-dimensional speed and equality to MATLAB
%******************************************************************************
%
%  MATLAB (R) is a trademark of The Mathworks (R) Corporation
%
%  Function:    mtimesx_test_nd
%  Filename:    mtimesx_test_nd.m
%  Programmer:  James Tursa
%  Version:     1.40
%  Date:        October 4, 2010
%  Copyright:   (c) 2009,2010 by James Tursa, All Rights Reserved
%
%  This code uses the BSD License:
%
%  Redistribution and use in source and binary forms, with or without
%  modification, are permitted provided that the following conditions are
%  met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
%  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%  POSSIBILITY OF SUCH DAMAGE.
%
%  Syntax:
%
%    A = mtimesx_test_nd    % default n=4 is used
%    A = mtimesx_test_nd(n)
%
%    where n = number of repetitions (should be 4 <= n <= 100)
%
%  Output:
%
%    Prints out speed and equality test results.
%    A = cell array with tabled results.
%
%  2010/Oct/04 --> 1.40, Added OpenMP support for custom code
%                        Expanded sparse * single and sparse * nD support
%
%--------------------------------------------------------------------------

function Cr = mtimesx_test_nd(n)
mtimesx; % load the mex routine into memory
if( nargin == 0 )
    n = 4;
else
    n = floor(n);
    if( ~(n >= 4 && n <= 100) )
        n = 4;
    end
end
cn = sprintf('%g',n);

disp(' ');
disp('MTIMESX multi-dimensional equality and speed tests');
disp('--------------------------------------------------');
disp(' ');
disp('(M x K) * ( K x N) equality tests, SPEED mode, M,K,N <= 4');
trans = 'NGTC';
cmpx = {'real ','cmpx '};
mtimesx('speed');
smallok = true;
for m=1:4
    for k=1:4
        for n=1:4
            for transa=1:4
                if( transa <= 2 )
                    ma = m;
                    ka = k;
                else
                    ma = k;
                    ka = m;
                end
                for transb=1:4
                    if( transb <= 2 )
                        kb = k;
                        nb = n;
                    else
                        kb = n;
                        nb = k;
                    end
                    for cmplxa=1:2
                        if( cmplxa == 1 )
                            A = floor(rand(ma,ka)*100+1);
                        else
                            A = floor(rand(ma,ka)*100+1) + floor(rand(ma,ka)*100+1)*1i;
                        end
                        for cmplxb=1:2
                            if( cmplxb == 1 )
                                B = floor(rand(kb,nb)*100+1);
                            else
                                B = floor(rand(kb,nb)*100+1) + floor(rand(kb,nb)*100+1)*1i;
                            end
                            Cm = mtimesx_sparse(A,trans(transa),B,trans(transb));
                            Cx = mtimesx(A,trans(transa),B,trans(transb));
                            if( isequal(Cm,Cx) )
                                disp(['(' cmpx{cmplxa} num2str(m) ' x ' num2str(k) ')' trans(transa) ...
                                        ' * (' cmpx{cmplxb} num2str(k) ' x ' num2str(n) ')' trans(transb) '  EQUAL']);
                            else
                                disp(['(' cmpx{cmplxa} num2str(m) ' x ' num2str(k) ')' trans(transa) ...
                                        ' * (' cmpx{cmplxb} num2str(k) ' x ' num2str(n) ')' trans(transb) '  NOT EQUAL']);
                                smallok = false;
                            end
                        end
                    end
                end
            end
        end
    end
end

if( mtimesx('openmp') )
disp(' ');
disp('(M x K) * ( K x N) equality tests, SPEEDOMP mode, M,K,N <= 4');
mtimesx('speedomp');
smallokomp = true;
for m=1:4
    for k=1:4
        for n=1:4
            for transa=1:4
                if( transa <= 2 )
                    ma = m;
                    ka = k;
                else
                    ma = k;
                    ka = m;
                end
                for transb=1:4
                    if( transb <= 2 )
                        kb = k;
                        nb = n;
                    else
                        kb = n;
                        nb = k;
                    end
                    for cmplxa=1:2
                        if( cmplxa == 1 )
                            A = floor(rand(ma,ka)*100+1);
                        else
                            A = floor(rand(ma,ka)*100+1) + floor(rand(ma,ka)*100+1)*1i;
                        end
                        A = reshape(repmat(A,1000,1),ma,ka,1000);
                        for cmplxb=1:2
                            if( cmplxb == 1 )
                                B = floor(rand(kb,nb)*100+1);
                            else
                                B = floor(rand(kb,nb)*100+1) + floor(rand(kb,nb)*100+1)*1i;
                            end
                            B = reshape(repmat(B,1000,1),kb,nb,1000);
                            Cm = mtimesx_sparse(A(:,:,1),trans(transa),B(:,:,1),trans(transb));
                            Cx = mtimesx(A,trans(transa),B,trans(transb));
                            if( isequal(Cm,Cx(:,:,1)) )
                                disp(['(' cmpx{cmplxa} num2str(m) ' x ' num2str(k) ')' trans(transa) ...
                                        ' * (' cmpx{cmplxb} num2str(k) ' x ' num2str(n) ')' trans(transb) '  EQUAL']);
                            else
                                disp(['(' cmpx{cmplxa} num2str(m) ' x ' num2str(k) ')' trans(transa) ...
                                        ' * (' cmpx{cmplxb} num2str(k) ' x ' num2str(n) ')' trans(transb) '  NOT EQUAL']);
                                smallokomp = false;
                            end
                        end
                    end
                end
            end
        end
    end
end
end

disp(' ');
if( smallok )
    disp('All small matrix multiplies are OK in SPEED mode');
else
    disp('ERROR --> One or more of the small matrix multiplies was not equal in SPEED mode');
end
if( mtimesx('openmp') )
if( smallokomp )
    disp('All small matrix multiplies are OK in SPEEDOMP mode');
else
    disp('ERROR --> One or more of the small matrix multiplies was not equal in SPEEDOMP mode');
end
end

disp(' ');
disp(['mtimesx multi-dimensional test routine using ' cn ' repetitions']);

if( mtimesx('OPENMP') )
    topm = 6;
else
    topm = 4;
end
Cr = cell(6,topm+1);
Cr{1,1} = 'All operands real';

for m=2:topm+1
if( m == 2 )
    mtimesx('BLAS');
elseif( m == 3 )
    mtimesx('LOOPS');
elseif( m == 4 )
    mtimesx('MATLAB');
elseif( m == 5 )
    mtimesx('SPEED');
elseif( m == 6 )
    mtimesx('LOOPSOMP');
else
    mtimesx('SPEEDOMP');
end
Cr{1,m} = mtimesx;

disp(' ');
disp('--------------------------------------------------------------');
disp('--------------------------------------------------------------');
disp(' ');
disp(['MTIMESX mode: ' mtimesx]);
disp(' ');
disp('(real 3x5x1x4x3x2x1x8) * (real 5x7x3x1x3x2x5) example');
Cr{2,1} = '(3x5xND) *(5x7xND)';
A = rand(3,5,1,4,3,2,1,8);
B = rand(5,7,3,1,3,2,5);
% mtimes
tm = zeros(1,n);
for k=1:n
clear Cm
A(1) = 2*A(1);
B(1) = 2*B(1);
tic
Cm = zeros(3,7,3,4,3,2,5,8);
for k1=1:3
    for k2=1:4
        for k3=1:3
            for k4=1:2
                for k5=1:5
                    for k6=1:8
                        Cm(:,:,k1,k2,k3,k4,k5,k6) = A(:,:,1,k2,k3,k4,1,k6) * B(:,:,k1,1,k3,k4,k5);
                    end
                end
            end
        end
    end
end
tm(k) = toc;
end
% mtimesx
tx = zeros(1,n);
for k=1:n
clear Cx
tic
Cx = mtimesx(A,B);
tx(k) = toc;
end
% results
tm = median(tm);
tx = median(tx);
if( tx < tm )
    faster = sprintf('%7.1f',100*(tm)/tx-100);
    slower = '';
else
    faster = sprintf('%7.1f',-(100*(tx)/tm-100));
    slower = ' (i.e., slower)';
end
Cr{2,m} = faster;
disp(' ');
disp(['mtimes  Elapsed time ' num2str(tm) ' seconds.']);
disp(['MTIMESX Elapsed time ' num2str(tx) ' seconds.']);
disp(['MTIMESX ' mtimesx ' mode is ' faster '% faster than MATLAB mtimes' slower])
if( isequal(Cx,Cm) )
    disp(['MTIMESX ' mtimesx ' mode result matches mtimes:  EQUAL'])
else
    dx = max(abs(Cx(:)-Cm(:)));
    disp(['MTIMESX ' mtimesx ' mode result does not match mtimes:  NOT EQUAL , max diff = ' num2str(dx)])
end

disp(' ');
disp('--------------------------------------------------------------');
disp('(real 3x3x1000000) * (real 3x3x1000000) example');
Cr{3,1} = '(3x3xN) *(3x3xN)';
A = rand(3,3,1000000);
B = rand(3,3,1000000);
% mtimes
tm = zeros(1,n);
for k=1:n
clear Cm
A(1) = 2*A(1);
B(1) = 2*B(1);
tic
Cm = zeros(3,3,1000000);
for k1=1:1000000
    Cm(:,:,k1) = A(:,:,k1) * B(:,:,k1);
end
tm(k) = toc;
end
% mtimesx
tx = zeros(1,n);
for k=1:n
clear Cx
tic
Cx = mtimesx(A,B);
tx(k) = toc;
end
% results
tm = median(tm);
tx = median(tx);
if( tx < tm )
    faster = sprintf('%7.1f',100*(tm)/tx-100);
    slower = '';
else
    faster = sprintf('%7.1f',-(100*(tx)/tm-100));
    slower = ' (i.e., slower)';
end
Cr{3,m} = faster;
disp(' ');
disp(['mtimes  Elapsed time ' num2str(tm) ' seconds.']);
disp(['MTIMESX Elapsed time ' num2str(tx) ' seconds.']);
disp(['MTIMESX ' mtimesx ' mode is ' faster '% faster than MATLAB mtimes' slower])
if( isequal(Cx,Cm) )
    disp(['MTIMESX ' mtimesx ' mode result matches mtimes:  EQUAL'])
else
    dx = max(abs(Cx(:)-Cm(:)));
    disp(['MTIMESX ' mtimesx ' mode result does not match mtimes:  NOT EQUAL , max diff = ' num2str(dx)])
end

disp(' ');
disp('--------------------------------------------------------------');
disp('(real 2x2x2000000) * (real 2x2x2000000) example');
Cr{4,1} = '(2x2xN) *(2x2xN)';
A = rand(2,2,2000000);
B = rand(2,2,2000000);
% mtimes
tm = zeros(1,n);
for k=1:n
clear Cm
A(1) = 2*A(1);
B(1) = 2*B(1);
tic
Cm = zeros(2,2,2000000);
for k1=1:2000000
    Cm(:,:,k1) = A(:,:,k1) * B(:,:,k1);
end
tm(k) = toc;
end
% mtimesx
tx = zeros(1,n);
for k=1:n
clear Cx
tic
Cx = mtimesx(A,B);
tx(k) = toc;
end
% results
tm = median(tm);
tx = median(tx);
if( tx < tm )
    faster = sprintf('%7.1f',100*(tm)/tx-100);
    slower = '';
else
    faster = sprintf('%7.1f',-(100*(tx)/tm-100));
    slower = ' (i.e., slower)';
end
Cr{4,m} = faster;
disp(' ');
disp(['mtimes  Elapsed time ' num2str(tm) ' seconds.']);
disp(['MTIMESX Elapsed time ' num2str(tx) ' seconds.']);
disp(['MTIMESX ' mtimesx ' mode is ' faster '% faster than MATLAB mtimes' slower])
if( isequal(Cx,Cm) )
    disp(['MTIMESX ' mtimesx ' mode result matches mtimes:  EQUAL'])
else
    dx = max(abs(Cx(:)-Cm(:)));
    disp(['MTIMESX ' mtimesx ' mode result does not match mtimes:  NOT EQUAL , max diff = ' num2str(dx)])
end

disp(' ');
disp('--------------------------------------------------------------');
disp('(real 2x2x2000000) * (real 1x1x2000000) example');
Cr{5,1} = '(2x2xN) *(1x1xN)';
A = rand(2,2,2000000);
B = rand(1,1,2000000);
% mtimes
tm = zeros(1,n);
for k=1:n
clear Cm
A(1) = 2*A(1);
B(1) = 2*B(1);
tic
Cm = zeros(2,2,2000000);
for k1=1:2000000
    Cm(:,:,k1) = A(:,:,k1) * B(:,:,k1);
end
tm(k) = toc;
end
% mtimesx
tx = zeros(1,n);
for k=1:n
clear Cx
tic
Cx = mtimesx(A,B);
tx(k) = toc;
end
% results
tm = median(tm);
tx = median(tx);
if( tx < tm )
    faster = sprintf('%7.1f',100*(tm)/tx-100);
    slower = '';
else
    faster = sprintf('%7.1f',-(100*(tx)/tm-100));
    slower = ' (i.e., slower)';
end
Cr{5,m} = faster;
disp(' ');
disp(['mtimes  Elapsed time ' num2str(tm) ' seconds.']);
disp(['MTIMESX Elapsed time ' num2str(tx) ' seconds.']);
disp(['MTIMESX ' mtimesx ' mode is ' faster '% faster than MATLAB mtimes' slower])
if( isequal(Cx,Cm) )
    disp(['MTIMESX ' mtimesx ' mode result matches mtimes:  EQUAL'])
else
    dx = max(abs(Cx(:)-Cm(:)));
    disp(['MTIMESX ' mtimesx ' mode result does not match mtimes:  NOT EQUAL , max diff = ' num2str(dx)])
end

try
bsxfun(@times,1,1);
Cr{6,1} = 'above vs bsxfun';
A = rand(2,2,2000000);
B = rand(1,1,2000000);
% bsxfun
tm = zeros(1,n);
for k=1:n
clear Cm
A(1) = 2*A(1);
B(1) = 2*B(1);
tic
Cm = bsxfun(@times,A,B);
tm(k) = toc;
end
% mtimesx
tx = zeros(1,n);
for k=1:n
clear Cx
tic
Cx = mtimesx(A,B);
tx(k) = toc;
end
% results
tm = median(tm);
tx = median(tx);
if( tx < tm )
    faster = sprintf('%7.1f',100*(tm)/tx-100);
    slower = '';
else
    faster = sprintf('%7.1f',-(100*(tx)/tm-100));
    slower = ' (i.e., slower)';
end
Cr{6,m} = faster;
disp(' ');
disp(['bsxfun  Elapsed time ' num2str(tm) ' seconds.']);
disp(['MTIMESX Elapsed time ' num2str(tx) ' seconds.']);
disp(['MTIMESX ' mtimesx ' mode is ' faster '% faster than MATLAB bsxfun with @times' slower])
if( isequal(Cx,Cm) )
    disp(['MTIMESX ' mtimesx ' mode result matches bsxfun with @times:  EQUAL'])
else
    dx = max(abs(Cx(:)-Cm(:)));
    disp(['MTIMESX ' mtimesx ' mode result does not match bsxfun with @times:  NOT EQUAL , max diff = ' num2str(dx)])
end
catch
    disp('Could not perform comparison with bsxfun, possibly because your version of');
    disp('MATLAB does not have it. You can download a substitute for bsxfun from the');
    disp('FEX here: http://www.mathworks.com/matlabcentral/fileexchange/23005-bsxfun-substitute');
end

end

disp(' ');
disp('Percent Faster Results Table');
disp(' ');
disp(Cr);

disp(' ');
disp('Done');
disp(' ');

end

author: martin-danelljan

repo: ECO


text:
clear all;
fprintf('Compiling perform_arithmetic_coding_mex and perform_arithmetic_coding_escape ... ');

% compile mex files
mex mex/perform_arithmetic_coding_mex.cpp mex/ac.cpp
mex mex/perform_arithmetic_coding_escape.cpp mex/coder/Arith.cpp mex/coder/BitIO.cpp mex/coder/IntCoding.cpp mex/coder/coder.cpp mex/coder/entropy.cpp mex/coder/global.cpp mex/coder/iHisto.cpp
mex mex/perform_arithmetic_coding_fixed.cpp mex/nr/arcmak.cpp mex/nr/arcode.cpp mex/nr/arcsum.cpp mex/nr/nrutil.cpp

disp('done.');

author: gpeyre

repo: matlab-toolboxes


text:
function [data g] = pre_diffData(varargin)
%
% Apply a difference filter to data. Differencing is a standard operation
% to improve stationarity of a time series. A first-order difference filter
% for input X is given by Y(t) = X(t) - X(t-1). This operation can be
% applied repeatedly to obtain an nth-order difference filter [1]
%
% Inputs:
%
%   EEG:        EEG data structure
%
% Optional:     <'Name',Value> pairs
%
%     VerbosityLevel:    Verbosity level. 0 = no output, 1 = text, 2 = graphical
%                        Possible values: 0,1,2
%                        Default value  : 1
%                        Input Data Type: real number (double)
%
%     DifferencingOrder: Differencing order
%                        Number of times to difference data
%                        Input Range  : [0  10]
%                        Default value: 1
%                        Input Data Type: real number (double)
% Outputs:
%
%   EEG:        processed EEG structure
%   g:          argument specification structure
%
% See Also: pop_pre_prepData(), pre_prepData(), diff()
%
% References:
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual. Section 6.5.1
%   Available at: http://www.sccn.ucsd.edu/wiki/Sift
%
% Author: Tim Mullen 2010, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

g = arg_define([0 1], varargin, ...
        arg_norep('data',mandatory), ...
        arg({'verb','VerbosityLevel'},int32(1),{int32(0) int32(1) int32(2)},'Verbosity level. 0 = no output, 1 = text, 2 = graphical'), ...
        arg({'difforder','DifferencingOrder'},1,[0 10],'Differencing order. Number of times to difference data'));

% commit data to workspace
data = g.data;
g = rmfield(g,'data');

[nchs npnt ntr] = size(data);

if g.difforder==0
    return;
end

if g.verb
    disp(['Differencing ' num2str(g.difforder) ' times...']);
end

if g.verb==2
    multiWaitbar('Differencing','Reset','Color',hlp_getNextUniqueColor);
end

datmp = zeros(nchs,npnt-g.difforder,ntr);

% difference each trial
for tr=1:ntr
    if g.verb==2
        multiWaitbar('Differencing',tr/ntr);
    end
    datmp(:,:,tr) = diff(data(:,:,tr),g.difforder,2);
end

% differencing reduces the number of datapoints by difforder
% so here we insert difforder random samples at the beginning
% of each time-series
noiseVar = 0.01*var(datmp(:));
data      = cat(2,noiseVar*randn(nchs,g.difforder,ntr),datmp);

if g.verb==2
    multiWaitbar('Differencing','Close');
end

repo: REST


text:
close all; clear all;

%% Setting of the problem
global s
pde = fracLapdata8;
pde.L = 1;
option.theta = 0.3;
option.estType = 'star';
option.maxIt = 17;
option.maxN = 5e4;
option.solver = 'mg';
option.tol = 1e-6;
[node,elem] = squaremesh([-1,1,-1,1],0.5);
bdFlag = setboundary(node,elem,'Dirichlet');

% %% s = 0.2
% s = 0.2;
% afemfracLap(node,elem,pde,bdFlag,option);
%
% %% s = 0.4
% s = 0.4;
% afemfracLap(node,elem,pde,bdFlag,option);

%% s = 0.6
s = 0.6;
afemfracLap(node,elem,pde,bdFlag,option);

%% s = 0.8
s = 0.8;
afemfracLap(node,elem,pde,bdFlag,option);

repo: ifem


text:
function hRD = plugin_rd_rx
% Reference design definition

% Call the common reference design definition function
hRD = AnalogDevices.adrv9361z7035.common.plugin_rd('ccbob_lvds','Rx & Tx');
AnalogDevices.adrv9361z7035.common.add_io(hRD, 'Rx & Tx', 'ccbob_lvds');

author: analogdevicesinc


text:
function model = rmSearchFit_twoGaussiansDoG(model,data,params,wProcess,t)
% rmSearchFit_twoGaussiansDoG - wrapper for 'fine' two DoG Gaussian fit
%
% model = rmSearchFit_twoGaussiansDoG(model,prediction,data,params);
%
% Second gaussian negative only.
%
% 2008/01 SOD: split of from rmSearchFit.

% add upper and lower limit:
expandRange    = params.analysis.fmins.expandRange;

% fminsearch options
searchOptions = params.analysis.fmins.options;
vethresh      = params.analysis.fmins.vethresh;

% convert to double just in case
params.analysis.X = double(params.analysis.X);
params.analysis.Y = double(params.analysis.Y);
params.analysis.allstimimages = double(params.analysis.allstimimages);

% amount of negative fits
nNegFit  = 0;
trends   = t.trends;
t_id     = t.dcid+2;

% get starting upper and lower range and reset TolFun
% (raw rss computation (similar to norm) and TolFun adjustments)
[range TolFun] = rmSearchFit_range(params,model,data);

% initialize
if ~isfield(model,'rss2')
    model.rss2 = zeros(size(model.rss));
end

if ~isfield(model,'rssPos')
    model.rsspos = zeros(size(model.rss));
end

if ~isfield(model,'rssNeg')
    model.rssneg = zeros(size(model.rss));
end

%-----------------------------------
% Go for each voxel
%-----------------------------------
progress = 0;tic;
for ii = 1:numel(wProcess),

    % progress monitor (10 dots)
    if floor(ii./numel(wProcess)*10)>progress,
        % print out estimated time left
        if progress==0,
            esttime = toc.*10;
            if floor(esttime./3600)>0,
                fprintf(1,'[%s]:Estimated processing time: %d voxels: %d hours.\n',...
                    mfilename,numel(wProcess),ceil(esttime./3600));
            else
                fprintf(1,'[%s]:Estimated processing time: %d voxels: %d minutes.\n',...
                    mfilename,numel(wProcess),ceil(esttime./60));
            end;
            fprintf(1,'[%s]:Nonlinear optimization (x,y,sigma):',mfilename);
        end;
        fprintf(1,'.');drawnow;
        progress = progress + 1;
    end;

    % volume index
    vi = wProcess(ii);
    vData = double(data(:,ii));

    % reset tolFun: Precision of evaluation function.
    % We define RMS improvement relative to the initial raw 'no-fit' data
    % RMS. So, 1 means stop if there is less than 1% improvement on the fit:
    searchOptions.TolFun = TolFun(ii);

    % actual fitting routine
    if searchOptions.MaxIter>0
        outParams = ...
            fmincon(@(x) rmModelSearchFit_twoGaussiansDoG(x,vData,...
            params.analysis.X,...
            params.analysis.Y,...
            params.analysis.allstimimages,trends),...
            range.start(:,vi),[],[],[],[],range.lower(:,vi),range.upper(:,vi),...
            @(x) distanceCon(x,range.start(:,vi),range.step(:,vi).*expandRange,params.analysis.minSigmaRatio),searchOptions);
    else
        outParams = range.start(:,vi);
    end

    % make predictions
    Xv = params.analysis.X - outParams(1);   % positive x0 moves center right
    Yv = params.analysis.Y - outParams(2);   % positive y0 moves center up
    rf(:,1) = exp( (Yv.*Yv + Xv.*Xv) ./ (-2.*(outParams(3).^2)) );
    rf(:,2) = exp( (Yv.*Yv + Xv.*Xv) ./ (-2.*(outParams(4).^2)) );

    % test 1 gaussian alone (must be done prior to with second gaussian)
    if searchOptions.MaxIter==0
        % without second gaussian
        X = [params.analysis.allstimimages*rf(:,1) trends];
        b = pinv(X)*vData;
        % force positive fit (is a constraint in the rmModelSearchFit_twoGaussian, so it should return here to get the right values)
        b(1) = abs(b(1));
        rss2  = norm(vData-X*b).^2;
        % make only the positive and the negative rf
        % find out where the rf is positive and where it is negative
        % for the positive rf put all the negative rf-values to zero, for
        % the negative rf put all the positive rf-values to zero.
        X = [params.analysis.allstimimages*rf trends];
        b    = pinv(X)*vData;

        % force positive fit
        b(1) = abs(b(1));
        %force negative b2 fit
        b(2) = -(abs(b(2)));

        rfBeta = b(1).*rf(:,1) + b(2).*rf(:,2);
        posInd = rfBeta > 0;
        negInd = rfBeta < 0;
        rfPos = rf;
        rfNeg = rf;
        rfPos(negInd,1) = 0;
        rfPos(negInd,2) = 0;
        rfNeg(posInd,1) = 0;
        rfNeg(posInd,2) = 0;
        XPos = [params.analysis.allstimimages*rfPos trends];
        XNeg = [params.analysis.allstimimages*rfNeg trends];
        rssPos = norm(vData-XPos*b).^2;
        rssNeg = norm(vData-XNeg*b).^2;

        %rfBeta = rf*b(1:2);
        %rssPos = norm(vData-max(rfBeta,0)).^2;
        %rssNeg = norm(vData-min(rfBeta,0)).^2;

%         if length(rssPos)==0
%             disp('rssPos is empty!');
%             rssPos = 1;
%         end
    end

    % with second gaussian
    X = [params.analysis.allstimimages*rf trends];
    b    = pinv(X)*vData;
    % force positive fit
    b(1) = abs(b(1));
    %force negative b2 fit
    b(2) = -(abs(b(2)));
    rss  = norm(vData-X*b).^2;

    % store results only if the first beta is positive, somehow fmincon
    % outputs negative fits. If the fit is negative keep old (grid) fit.
    if b(1)>0 && b(1)>-b(2) && b(2) <=0,
        model.x0(vi)   = outParams(1);
        model.y0(vi)   = outParams(2);
        model.s(vi)    = outParams(3);
        model.s_major(vi)    = outParams(3);
        model.s_minor(vi)    = outParams(3);
        model.s_theta(vi)    = 0;
        model.s2(vi)   = outParams(4);
        model.rss(vi)  = rss;
        model.b([1 2 t_id],vi) = b;
        if searchOptions.MaxIter==0
            model.rss2(vi) = rss2;
            model.rsspos(vi) = rssPos;
            model.rssneg(vi) = rssNeg;
        end
    else
        % change the percent variance explained to be just under the
        % current vethresh. So it counts as a 'coarse'-fit but can still be
        % included in later 'fine'-fits
        model.rss(vi)  = (1-max((vethresh-0.01),0)).*model.rawrss(vi);
        nNegFit = nNegFit + 1;
        if searchOptions.MaxIter==0
            model.rss2(vi)  = (1-max((vethresh-0.01),0)).*model.rawrss(vi);
            model.rsspos(vi)  = (1-max((vethresh-0.01),0)).*model.rawrss(vi);
            model.rssneg(vi)  = (1-max((vethresh-0.01),0)).*model.rawrss(vi);
        end
    end;
end;

% end time monitor
et  = toc;
if floor(et/3600)>0,
    fprintf(1,'Done [%d hours].\n',ceil(et/3600));
else
    fprintf(1,'Done [%d minutes].\n',ceil(et/60));
end;
fprintf(1,'[%s]:Removed negative fits: %d (%.1f%%).\n',...
    mfilename,nNegFit,nNegFit./numel(wProcess).*100);

return;

%-----------------------------------
% make sure that the pRF can only be moved "step" away from original
% poisiton "startParams"
% For the two Gaussian model we add the additional constraint that the
% second Gaussian is at least twice as large as the first.
function [C, Ceq]=distanceCon(x,startParams,step,minRatio)
Ceq = [];
dist = x([1 2])-startParams([1 2]);
C(1) = sqrt(dist(1).^2+dist(2).^2) - step;
C(2) = minRatio - 0.001 - x(4)./x(3);
return;
%-----------------------------------

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [bestEpsilon, bestF1] = selectThreshold(yval, pval)
    %% SELECTTHRESHOLD Find the best threshold (epsilon) to use for selecting
    %   outliers
    %   [bestEpsilon bestF1] = SELECTTHRESHOLD(yval, pval) finds the best
    %   threshold to use for selecting outliers based on the results from a
    %   validation set (pval) and the ground truth (yval).
    %

    bestEpsilon = 0;
    bestF1 = 0;
    F1 = 0;

    stepsize = (max(pval) - min(pval)) / 1000;
    for epsilon = min(pval) : stepsize : max(pval)
        % ====================== YOUR CODE HERE ======================
        % Instructions: Compute the F1 score of choosing epsilon as the
        %               threshold and place the value in F1. The code at the
        %               end of the loop will compare the F1 score for this
        %               choice of epsilon and set it to be the best epsilon if
        %               it is better than the current choice of epsilon.
        %
        % Note: You can use predictions = (pval < epsilon) to get a binary vector
        %       of 0's and 1's of the outlier predictions
        predictions = (pval < epsilon);  % entry = 1 if anomalous

        % calculate how well we did
        tp = sum(predictions == 1 & yval == 1);
        fp = sum(predictions == 1 & yval == 0);
        fn = sum(predictions == 0 & yval == 1);

        % calculate precision and recall
        precision = tp / (tp + fp);
        recall = tp / (tp + fn);
        F1 = (2 * precision * recall) / (precision + recall);

        % =============================================================

        if F1 > bestF1
           bestF1 = F1;
           bestEpsilon = epsilon;
        end

        %fprintf('tp=%d, fp=%d, fn=%d, precision=%f, recall=%f, best=%f was F1=%f\n', ...
        %    tp, fp, fn, precision, recall, bestEpsilon, bestF1);
    end
end

author: worldveil

repo: coursera-ml


text:
function [x,fval,exitflag,info] = opti_lmder(fun,grad,x0,ydata,opts)
%OPTI_LMDER Solve a NLS using LM_DER (Levenberg-Marquardt MINPACK Routine)
%
%   min sum[ (F(x) - ydata)^2 ]
%    x
%
%   x = opti_lmder(fun,grad,x0,ydata) solves a NLS where fun is the fitting
%   function. grad is an optional gradient of the fitting function and x0
%   is a starting guess. ydata is the data to fit the function to. Use
%   ydata = zeros() for solving a SNLE.
%
%   x = opti_lmder(fun,grad,x0,ydata,opts) uses opts to pass optiset
%   options to the solver.
%
%   [x,fval,exitflag,info] = opti_lmder(...) returns the objective value at
%   the solution, together with the solver exitflag, and an information
%   structure.
%
%   THIS IS A WRAPPER FOR LM_DER + LM_DIF
%   See supplied License

%   Copyright (C) 2012 Jonathan Currie (I2C2)

if(nargin < 5), opts = optiset; end
if(nargin < 4), error('LM_DER requires at least 4 arguments'); end

%Setup display level
opts.display = dispLevel(opts.display);

%Check we have a valid x0
if(isempty(x0) || any(isnan(x0)))
    error('LMDER requires an initial guess, x0!');
end

t = tic;
% Run LM_DER
[x, fval, exitflag, feval] = lmder(fun,grad,x0,ydata,opts);

%Collect Results
info.FuncEvals = feval;
info.Time = toc(t);
info.Algorithm = 'LM_DER: MINPACK Levenberg-Marquardt';

switch(exitflag)
    case 1
        info.Status = 'Optimal';
    case 0
        info.Status = 'Exceeded Function Evaluations';
    case -1
        info.Status = 'Infeasible / Could not Converge';
    case -2
        info.Status = 'LM_DER Error';
    case -5
        info.Status = 'User Exited';
    otherwise
        info.Status = 'LM_DER Error';
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus