text:
function mov = aviread_rawy8(filename,varargin)
%AVIREAD Read AVI file.
%   AVIREAD will be removed in a future release. Use VIDEOREADER instead.
%
%   MOV = AVIREAD(FILENAME) reads the AVI movie FILENAME into the
%   MATLAB movie structure MOV. If FILENAME does not include an
%   extension, then '.avi' will be used. MOV has two fields, "cdata"
%   and "colormap". If the frames are truecolor images, MOV.cdata will
%   be Height-by-Width-by-3 and MOV.colormap will be empty. If the
%   frames are indexed images, then MOV.cdata field will be
%   Height-by-Width and MOV.colormap will be M-by-3. On UNIX, FILENAME
%   must be an uncompressed AVI file.
%
%   MOV = AVIREAD(...,INDEX) reads only the frame(s) specified by
%   INDEX.  INDEX can be a single index or an array of indices into the
%   video stream, where the first frame is number one.
%
%   The supported frame types are 8-bit (indexed or grayscale),
%   16-bit grayscale or 24-bit (TrueColor) frames.
%
%   See also VIDEOREADER, VIDEOWRITER, MOVIE.

%   Copyright 1984-2011 The MathWorks, Inc.
%   $Revision: 1.1.6.18.2.1 $  $Date: 2011/12/13 18:36:41 $

%warning(message('MATLAB:audiovideo:aviread:FunctionToBeRemoved'));

% Initialization
map = [];
index = -10000;

% Validate input/output.
error(nargoutchk(1,1,nargout, 'struct'));
error(nargchk(1,2,nargin, 'struct'));

if nargin == 2
    if isnumeric(varargin{1})
        index = varargin{1};
    else
        error(message('MATLAB:audiovideo:aviread:nonNumericIndex'));
    end
end

if ~ischar(filename)
    error(message('MATLAB:audiovideo:aviread:invalidFilename'));
end

[path,name,ext] = fileparts(filename);
if isempty(ext)
    filename = strcat(filename,'.avi');
end

info = getAviInfo( filename );

if ismember(lower(info.VideoFrameHeader.CompressionType),{'y800','y8'}),
  info.VideoFrameHeader.CompressionType = 'RAW ';
end

if (index(1) ~= -10000)
    if ( isfield(info.MainHeader,'HasIndex') == 0 )
        error(message('MATLAB:audiovideo:aviread:indexParameterNotSupported', filename));
    elseif any(index > info.MainHeader.TotalFrames)
        error(message('MATLAB:audiovideo:aviread:indexOutOfRange', info.MainHeader.TotalFrames));
    elseif any(index <= 0)
        error(message('MATLAB:audiovideo:aviread:invalidZeroIndex'));
    end
end

% if ispc
%     try
%         if index == -10000
%             X = readavi(info.Filename,-1);
%         else
%             X = readavi(info.Filename,index-1);
%         end
%     catch exception
%         if strcmpi(exception.identifier,'MATLAB:audiovideo:aviread:invalidcodec') == 1 && ...
%             strcmpi(info.VideoFrameHeader.CompressionType,'indeo5')
%             % if we get an invalid codec error and are asking for indeo5
%             % inform the user of the problem and link to mathworks.com
%             error(message('MATLAB:audiovideo:aviread:indeo5NotFound'));
%         else
%             throw(exception);
%         end
%     end
%
%     % The width of the frames may be padded to be on 4 byte boundaries because
%     % this is how bitmaps are typically written.  However, AVI files do not
%     % restrict the frames to lie on 4 byte boundries.  Because of
%     % inconsistencies with the header information in AVI files, the READAVI
%     % MEX-function determines if the frame is padded and reads the correct
%     % number of bytes for each frame.  If the returned data is more than
%     % expected (different for indexed and RGB) then the data is padded.  The
%     % padding is removed before frames are returned by AVIREAD.
%
%     width = info.VideoFrameHeader.Width;
%     if ((info.VideoFrameHeader.Width*info.VideoFrameHeader.Height ~= length(X(1).cdata)) && (info.VideoFrameHeader.BitDepth == 8)) || ((info.VideoFrameHeader.Width*info.VideoFrameHeader.Height*3 ~= length(X(1).cdata)) && (info.VideoFrameHeader.BitDepth == 24) )
%         paddedWidth = 4*ceil(width/4);
%     else
%         paddedWidth = width;
%     end
%     height = info.VideoFrameHeader.Height;
%
%     if info.VideoFrameHeader.BitDepth == 8 || info.VideoFrameHeader.BitDepth == 16
%         % Indexed frames
%         for i=1:length(X)
%             if height<0
%                 % Movie frames are stored top down (height is negative).
%                 mov(i).cdata = reshape(X(i).cdata, paddedWidth, abs(height))';
%             else
%                 mov(i).cdata = rot90(reshape(X(i).cdata, paddedWidth, height));
%             end
%             if paddedWidth ~= width
%                 mov(i).cdata = mov(i).cdata(:,1:width);
%             end
%             map = reshape(X(i).colormap, 4, ...
%                 info.VideoFrameHeader.NumColormapEntries);
%             map = double(flipud(map(1:3,:))')/255;
%             mov(i).colormap = map;
%         end
%     elseif info.VideoFrameHeader.BitDepth == 24
%         paddedWidth = 4*ceil(width*3/4);
%         % Truecolor frames
%         for i=1:length(X)
%             f = X(i).cdata;
%             if height<0
%                 % Movie frames are stored top down (height is negative).
%                 f = permute(reshape(f, paddedWidth, abs(height)),[2 1 3]);
%             else
%                 f = rot90(reshape(f, paddedWidth,height));
%             end
%             if paddedWidth ~= width
%                 f = f(:,1:width*3);
%             end
%             % Movie frames are stored top down (height is negative).
%             RGB(1:abs(height), 1:width,3) = f(:,1:3:end);
%             RGB(:, :, 2) = f(:,2:3:end);
%             RGB(:, :, 1) = f(:,3:3:end);
%             mov(i).cdata = RGB;
%             mov(i).colormap = [];
%         end
%     else
%         error(message('MATLAB:audiovideo:aviread:invalidColorBitDepth'));
%     end
% end %End of PC specific code

% if isunix

    if isempty(strmatch(lower(info.VideoFrameHeader.CompressionType),...
            {'dib ', 'raw ','none','raw ',char([0 0 0 0])}))
        error(message('MATLAB:audiovideo:aviread:compressedFileInvalid'));
    end

    if strcmpi(info.VideoFrameHeader.CompressionType,char([0 0 0 0]))
        info.VideoFrameHeader.CompressionType = 'none';
    end

    fid = fopen(filename,'r','l');
    if fid == -1
        error(message('MATLAB:audiovideo:aviread:unableToOpenFile', filename));
    end

    % Find RIFF chunk
    [chunk, msg,msgID] = findchunk(fid,'RIFF');
    errorWithFileClose(msgID,msg,fid);

    % Read AVI chunk
    [rifftype,msg,msgID] = readfourcc(fid);
    errorWithFileClose(msgID,msg,fid);
    if ( strcmp(rifftype,'AVI ') == 0 )
        errorWithFileClose('MATLAB:audiovideo:aviread:MissingAVIChunk',xlate('Not a valid AVI file. Missing ''AVI '' chunk.'), fid);
    end

    % Find hdrl LIST chunk
    [hdrlsize,msg,msgID] = findlist(fid,'hdrl');
    errorWithFileClose(msgID,msg,fid);

    % Find and skip avih chunk
    [chunk,msg,msgID] = findchunk(fid,'avih');
    errorWithFileClose(msgID,msg,fid);
    [msg,msgID] = skipchunk(fid,chunk);
    errorWithFileClose(msgID,msg,fid);

    % Find the video stream
    for  i = 1:info.MainHeader.NumStreams
        % Find strl LIST chunk
        [strlsize,msg,msgID] = findlist(fid,'strl');
        errorWithFileClose(msgID,msg,fid);
        % Read strh chunk
        [strhchunk,msg,msgID] = findchunk(fid,'strh');
        errorWithFileClose(msgID,msg,fid);
        % Determine stream type
        streamType = readfourcc(fid);
        % Break if it is a video stream
        if(strcmp(streamType,'vids'))
            found = 1;
            break;
        else
            found  = 0;
            % Seek to end of strl list minus the amount read
            if ( fseek(fid,listsize - 16,0) == -1 )
                errorWithFileClose('MATLAB:audiovideo:aviread:invalidChunkSize',getString(message('MATLAB:audiovideo:aviread:invalidChunkSizeAVI')), fid);
            end
        end
    end

    if (found == 0)
        errorWithFileClose('MATLAB:audiovideo:aviread:noVideoStreamFound',getString(message('MATLAB:audiovideo:aviread:noVideoStreamFound')), fid);
    end

    % Skip the strh chunk minus the fourcc (4 bytes) already read.
    strhchunk.cksize = strhchunk.cksize - 4;
    [msg,msgID] = skipchunk(fid,strhchunk);
    errorWithFileClose(msgID,msg,fid);

    % Read strf chunk
    [strfchunk,msg,msgID] = findchunk(fid,'strf');
    errorWithFileClose(msgID,msg,fid);

    if info.VideoFrameHeader.BitDepth == 24
        % For TrueColor images, skip the Stream Format chunk
        [msg,msgID] = skipchunk(fid,strfchunk);
        errorWithFileClose(msgID,msg,fid);
    elseif  info.VideoFrameHeader.BitDepth == 8
        % If bitmap has a palette Seek past the BITMAPINFOHEADER to put the
        % file pointer at the beginning of the colormap
        if  fseek(fid,info.VideoFrameHeader.BitmapHeaderSize,0) == -1
            errorWithFileClose('MATLAB:audiovideo:aviread:invalidBitmapInfoHeaderSize',getString(message('MATLAB:audiovideo:aviread:invalidBitmapInfoHeaderSize')),fid);
        end
        map = readColormap(fid,info.VideoFrameHeader.NumColorsUsed);
    else
        errorWithFileClose('MATLAB:audiovideo:aviread:invalidColorBitDepth',getString(message('MATLAB:audiovideo:aviread:invalidColorBitDepthBitmap')), fid);
    end

    % Search for the movi LIST
    [movisize,msg,msgID] = findlist(fid,'movi');
    errorWithFileClose(msgID,msg,fid);
    % movioffset will be used when using idx1. The offsets stored in idx1 are
    % with respect to just after the 'movi' LIST (not including 'movi')
    movioffset = ftell(fid) -4;

    %Determine method of reading movie
    if ( index ~= -10000 )
        method = 'UseIndex';
    elseif ( isfield(info,'IsInterleaved') == 1 )
        method = 'Interleaved';
    else
        method = 'Normal';
    end

    totalReadFrames = 1;
    switch method
        case 'Interleaved'
            % Read movies that are interleaved (contain rec lists)
            while(totalReadFrames < info.MainHeader.TotalFrames )
                [recsize, msg] = findlist(fid,'rec ');
                currentPos = ftell(fid);
                % The recListPos is the current position plus the rec list size minus
                % 4 bytes because we read the four character LIST name in findlist
                recListEndPos = currentPos + recsize - 4;
                chunk = readchunk(fid);
                while(~strcmpi(chunk.ckid,'00db') && ~strcmpi(chunk.ckid,'00dc'))
                    [msg,msgID] = skipchunk(fid,chunk);
                    errorWithFileClose(msgID,msg,fid);
                    if (ftell(fid)==recListEndPos)
                        [recsize,msg,msgID] = findlist(fid,'rec ');
                        errorWithFileClose(msgID,msg,fid);
                        currentPos = ftell(fid);
                        recListEndPos = currentPos + recsize - 4; % minus 4 because we read
                    end
                    chunk = readchunk(fid);
                end

                % Prepare input for readbmpdata
                tempinfo = info;
                tempinfo.ImageDataOffset = ftell(fid);
                tempinfo.CompressionType = info.VideoFrameHeader.CompressionType;
                tempinfo.BitDepth = info.VideoFrameHeader.BitDepth;
                % readbmpdata opens and closes the file so aviread must also to
                % have a valid fid.
                status = fclose(fid);

                RGB = readbmpdata(tempinfo);

                fid = fopen(filename,'r','l');
                fseek(fid,tempinfo.ImageDataOffset,'bof');
                % readbmpdata does not move the file pointer
                [msg,msgID] = skipchunk(fid,chunk);
                errorWithFileClose(msgID,msg,fid);
                % Assign to MATLAB movie structure
                mov(totalReadFrames).cdata = RGB;
                totalReadFrames = totalReadFrames+1;
            end
        case 'Normal'
            % Find each of the 00db or 00dc chunks and read the frames
            for i = 1:info.MainHeader.TotalFrames
                chunk = readchunk(fid);
                while (~strcmpi(chunk.ckid,'00db') && ~strcmpi(chunk.ckid,'00dc') );
                    [msg,msgID] = skipchunk(fid,chunk);
                    errorWithFileClose(msgID,msg,fid);
                    chunk = readchunk(fid);
                end

                % Prepare data to be sent to readbmpdata
                tempinfo.Filename = info.Filename;
                tempinfo.Width = info.VideoFrameHeader.Width;
                tempinfo.Height = info.VideoFrameHeader.Height;
                tempinfo.ImageDataOffset = ftell(fid);
                tempinfo.CompressionType = 'none';
                tempinfo.BitDepth = info.VideoFrameHeader.BitDepth;
                % readbmpdata opens and closes the file so aviread must also to
                % have a valid fid
                status = fclose(fid);
                % Read RGB frame
                RGB = readbmpdata(tempinfo);

                fid = fopen(filename,'r','l');
                fseek(fid,tempinfo.ImageDataOffset,'bof');
                [msg,msgID] = skipchunk(fid,chunk);
                errorWithFileClose(msgID,msg,fid);
                % Assign to MATLAB movie structure
                mov(i).cdata = RGB;
            end
        case 'UseIndex'
            % Skip the movi LIST (minus 4 because 'movi' was read) and use idx1.
            if ( fseek(fid,movisize-4,0) == -1 )
                errorWithFileClose('MATLAB:audiovideo:aviread:InvalidChunkSize',getString(message('MATLAB:audiovideo:aviread:InvalidChunkSizeWAV')), fid);
            end
            % Find idx1 chunk
            [idx1chunk,msg,msgID] = findchunk(fid,'idx1');
            errorWithFileClose(msgID,msg,fid);
            idx1ChunkPos = ftell(fid);

            for j = 1:length(index)
                fseek(fid,idx1ChunkPos,'bof');
                for i = 1:index(j)
                    found = 0;
                    while(found == 0)
                        id = readfourcc(fid);
                        if (strcmpi(id,'00db') || strcmpi(id,'00dc'))
                            found = 1;
                        end
                        [idx1data,msg,msgID] = readIDX1(fid);
                        errorWithFileClose(msgID,msg,fid);
                    end

                    % If the very first index inside the 'idx1' chunk contains a
                    % data offset that is larger than the 'movi' list offset, then
                    % it is an absolute offset.  Otherwise, we will calculate the
                    % offset relative to the beginning of the 'movi' list.
                    if (i == 1)
                        if (idx1data.offset > movioffset)
                            isIdx1OffsetAbsolute = true;
                        else
                            isIdx1OffsetAbsolute = false;
                        end
                    end
                end

                % Prepare data to be sent to readbmpdata
                tempinfo.Filename = info.Filename;
                tempinfo.Width = info.VideoFrameHeader.Width;
                tempinfo.Height = info.VideoFrameHeader.Height;
                tempinfo.BitDepth = info.VideoFrameHeader.BitDepth;
                % 8 is the riffheadersize

                % Add the 'movi' offset only if the data offset is not absolute.
                if isIdx1OffsetAbsolute
                    tempinfo.ImageDataOffset = idx1data.offset + 8;
                else
                    tempinfo.ImageDataOffset = movioffset + idx1data.offset + 8;
                end
                %tempinfo.CompressionType = info.VideoFrameHeader.CompressionType;
                tempinfo.CompressionType = 'none';
                % readbmpdata opens and closes the file so aviread must also to
                % have a valid fid
                currPos = ftell(fid);
                status = fclose(fid);

                frame = readbmpdata(tempinfo);

                fid = fopen(filename,'r','l');
                fseek(fid,currPos,'bof');
                mov(j).cdata = frame;
            end
    end
    fclose(fid);
    % Formulate outputs
    [mov.colormap] = deal(map);
    varargout{1} = mov;
%end
return

function map = readColormap(fid,numColors)
% Read colormap for 8-bit indexed images
map = fread(fid,numColors*4,'*uint8');
map = reshape(map,4,numColors);
map = double(flipud(map(1:3,:))')/255;
return;

function [idx1data,msg,msgID] = readIDX1(fid)
% Read the data in the idx1 chunk.
msg = '';
msgID='';
[idx1data.flags, count] = fread(fid,1,'uint32');
if ( count ~= 1 )
    msgID = 'MATLAB:audiovideo:aviread:IncorrectIDX1ChunkSize';
    msg = getString(message(msgID));
end
[idx1data.offset, count] = fread(fid,1,'uint32');
if ( count ~= 1 )
    msgID = 'MATLAB:audiovideo:aviread:IncorrectIDX1ChunkSize';
    msg = getString(message(msgID));
end
[idx1data.length, count] = fread(fid,1,'uint32');
if ( count ~= 1 )
    msgID = 'MATLAB:audiovideo:aviread:IncorrectIDX1ChunkSize';
    msg = getString(message(msgID));
end
return;

function errorWithFileClose(msgID,msg,fid)
%Close open file the error
if ~isempty(msgID)
    fclose(fid);
    error(msgID,msg);
end
return;

function info = getAviInfo( filename )
    aviinfoWarningID = 'MATLAB:audiovideo:aviinfo:FunctionToBeRemoved';
    S = warning('OFF', aviinfoWarningID);
    cleaner = onCleanup(@()warning(S));

    % Save the state of lastwarn
    [preWarningMsg, preWarningID] = lastwarn();

    info = aviinfo(filename,'Robust');

    [~, lastWarningID] = lastwarn();

    % If aviinfo warns about deprecation restore
    % the penultimate last warning
    if (strcmp(lastWarningID, aviinfoWarningID))
        lastwarn( preWarningMsg, preWarningID );
    end

author: kristinbranson

repo: JAABA


text:
% eeg_mergelocs() - merge channel structure while preserving channel
%                   order
%
%      >> mergedlocs = eeg_mergelocs(loc1, loc2, loc3, ...);
%
% Inputs:
%     loc1     - EEGLAB channel location structure
%     loc2     - second EEGLAB channel location structure
%
% Output:
%     mergedlocs - merged channel location structure
%     warning    - [0|1] dissimilar structures found (0=false, 1=true)
%
% Author: Arnaud Delorme, August 2006

% Copyright (C) Arnaud Delorme, CERCO, 2006, arno@salk.edu
%
% 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 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

function [alllocs warn] = eeg_mergelocs(varargin)

persistent warning_shown;
warn = 0;

try
    % sort by length
    % --------------
    len = cellfun(@length, varargin);
    [tmp so] = sort(len, 2, 'descend');
    varargin = varargin(so);

    alllocs = varargin{1};
    for index = 2:length(varargin)

        % fuse while preserving order (assumes the same channel order)
        % ------------------------------------------------------------
        tmplocs = varargin{index};
        newlocs = myunion(alllocs, tmplocs);

        if length(newlocs) > length(union({ alllocs.labels }, { tmplocs.labels }))

            warn = 1;
            if isempty(warning_shown)
                disp('Warning: different channel montage or electrode order for the different datasets');
                warning_shown = 1;
            end;

            % trying to preserve order of the longest array
            %----------------------------------------------
            if length(alllocs) < length(tmplocs)
                tmp     = alllocs;
                alllocs = tmplocs;
                tmplocs = tmp;
            end;
            allchans = { alllocs.labels tmplocs.labels };
            [uniquechan ord1 ord2 ]  = unique_bc( allchans );

            [tmp rminds] = intersect_bc( uniquechan, { alllocs.labels });
            ord1(rminds) = [];
            tmplocsind = ord1-length(alllocs);

            newlocs = [ alllocs tmplocs(tmplocsind) ];

        end;
        alllocs = newlocs;
    end;
catch,
    % temporary fix for dissimilar structures
    % should check channel structure consistency instead
    % using checkchan function
    disp('Channel merging warning: dissimilar fields in the two structures');
    [alllocs warn ] = eeg_mergelocs_diffstruct(varargin{:});
end;

% Checking consistency of chanlocs
alllocs = eeg_checkchanlocs(alllocs);

% union of two channel location structure
% without loosing the order information
% ---------------------------------------
function alllocs = myunion(locs1, locs2)

   labs1 = { locs1.labels };
   labs2 = { locs2.labels };

   count1 = 1;
   count2 = 1;
   count3 = 1;
   alllocs = locs1; alllocs(:) = [];
   while count1 <= length(locs1) || count2 <= length(locs2)

       if count1 > length(locs1)
           alllocs(count3) = locs2(count2);
           count2 = count2 + 1;
           count3 = count3 + 1;
       elseif count2 > length(locs2)
           alllocs(count3) = locs1(count1);
           count1 = count1 + 1;
           count3 = count3 + 1;
       elseif strcmpi(labs1{count1}, labs2{count2})
           alllocs(count3) = locs1(count1);
           count1 = count1 + 1;
           count2 = count2 + 1;
           count3 = count3 + 1;
       elseif isempty(strmatch(labs1{count1}, labs2, 'exact'))
           alllocs(count3) = locs1(count1);
           count1 = count1 + 1;
           count3 = count3 + 1;
       else
           alllocs(count3) = locs2(count2);
           count2 = count2 + 1;
           count3 = count3 + 1;
       end;

   end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function test_bug2444

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_selectdata ft_selectdata_new

%% make some data

comp = [];
comp.label = {'comp1', 'comp2'};
comp.topolabel = {'chan1', 'chan2'};
comp.time = {1:1000, 1:1000};
comp.trial = {randn(2,1000), randn(2,1000)};
comp.topo = eye(2);
comp.unmixing = eye(2);

%% select some data

cfg = [];
cfg.latency = [10 20];
selected = ft_selectdata(cfg, comp);

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [varargout]=tricontf(varargin);
% TRICONT  Filled contours data on a triangular mesh
%     [CS,h]=TRICONTF(X,Y,M,Z) takes the mesh specified by points
%     (X,Y), with heights Z, and the  Nx3 triangulation M (where
%     each row gives the indexes into X/Y vectors of the triangle
%     corners), and make filled contours.
%
%     TRICONTF(...,LEVELS) fills contours at the specified levels.
%
%
%     TRICONTF(...,LINESPEC) takes line style and colour
%     parameters for the lines.
%
%     [CS,h] are the inputs needed by CLABEL, but if you want
%     labelled contours you are better off doing something like
%
%
%     [CS,h]=TRICONTF(...);
%     set(h,'edgecolor','none');
%     hold on;
%     [CS,h]=TRICONT(...);
%     hold off
%     clabel(CS,h)
%
%     since the CS and h returned by TRICONTF contain boundary
%     information which is not usually required for labelling.
%
%     The advantage of using TRICONTF over the GRIDDATA method is the
%     the triangulation is already specified (and may be non-convex
%     as well as containing holes), and the contours are also exact
%     on the triangulation (no weird boundary jaggies). However,
%     LINEAR finite elements are assumed.
%
%     Note - unlike contourf, TRICONTF will not
%     work 'properly' with NaN values in Z - if you
%     have bad data remove those triangles from M!
%
%     See also TRICONT
%
%     Rich Pawlowicz (rpawlowicz@eos.ubc.ca)  March/2013
%

% Copy a bunch of lines from contourf
error(nargchk(4,9,nargin,'struct'));
[cax,args,nargs] = axescheck(varargin{:});

cax = newplot(cax);

% Check for empty arguments.
for i = 1:nargs
  if isempty(args{i})
    error('MATLAB:contourf:EmptyInput','Input argument is empty');
  end
end

% Trim off the last arg if it's a string (line_spec).
nin = nargs;
if ischar(args{end})
  [lin,col,mark,msg] = colstyle(args{end}); %#ok
  if ~isempty(msg), error(msg); end %#ok
  nin = nin - 1;
else
  lin = '';
  col = '';
end

% Closed contours are made up of a) contours through the
% interior, and b) curves around the grid boundaries (there
% can be interior "islands" in an FEM mesh too).  The
% general strategy is to get a bunch of line segments, and
% then join them into continuous curves.

% First, get all the interior contours (this is pretty fast)
% but don't draw anything yet; this also does some checking
% on M to make sure triangles are oriented.

[Xp,Yp,M,Zp,nv,CS,xx,yy,zz]=tricont(args{1:nin});
nCS=[find(isnan(CS(2,:))) size(CS,2)+1];
lCS=CS(1,nCS(1:end-1));

% Don't fill contours below the lowest level specified in nv.
% To fill all contours, specify a value of nv lower than the
% minimum of the surface.
i = find(isfinite(Zp));
minz = min(Zp(i));
maxz = max(Zp(i));
draw_min=0;
if any(nv <= minz),
  draw_min=1;
end

% Second step - get all the boundary curves.

BS=findboundary(M);
iBS=find(isnan(BS));

% Third step:
% Once we have interior and boundary curves, join *them* together.
% However, we have to have separate boundary curves for every
% level, because we need only the parts of the boundary where the
% interior is higher than that level.

fCS=NaN(2,size(CS,2)+size(BS,2));iCS=1;
ncurves = 0;
I = [];
Area=[];levs=[];

for l=1:length(nv),   % For each level

  lvlBS=NaN(2,length(CS)+length(BS)+50);jBS=1;

  iCSlvl=find(lCS==nv(l));
  if any(iCSlvl),

     % Create a structure with all the lines at level
     % nv(k). First add the interior contours

     for k=1:length(iCSlvl),
       nseg=nCS(iCSlvl(k)+1)-nCS(iCSlvl(k))-1;
       level=CS(1,nCS(iCSlvl(k)));
       xseg=CS(1,nCS(iCSlvl(k))+[1:nseg]);
       yseg=CS(2,nCS(iCSlvl(k))+[1:nseg]);
       lvlBS(:,jBS+[0:nseg])=[ nv(l) xseg ; NaN yseg ];
       jBS=jBS+1+nseg;
     end;
  end;

     % Now look through all the boundaries, and for each boundary segment
     % take only the parts that are higher than nv(l). In general we will
     % have to interpolate to get the first and last position of the curve,
     % except when it it as the beginning or end of a boundary curve.

     for k=1:length(iBS)-1,
       Bseg=BS(iBS(k)+1:iBS(k+1)-1);
       id=Zp(Bseg)>nv(l);

       segstart=find(diff(id)== 1)+1;
       if id(1)==1,segstart=[1;segstart]; end;
       segend=  find(diff(id)==-1);
       if id(end)==1,segend=[segend;length(Bseg)]; end;

       for m=1:length(segstart),
	 if segstart(m)==1,  % First point is above...
	   xstart=[];
	   ystart=[];
	 else                % ...otherwise interpolate it
	   xstart=(Xp(Bseg(segstart(m)))*(nv(l)-Zp(Bseg(segstart(m)-1))) + Xp(Bseg(segstart(m)-1))*(Zp(Bseg(segstart(m)))-nv(l)) )./(Zp(Bseg(segstart(m)))-Zp(Bseg(segstart(m)-1)));
	   ystart=(Yp(Bseg(segstart(m)))*(nv(l)-Zp(Bseg(segstart(m)-1))) + Yp(Bseg(segstart(m)-1))*(Zp(Bseg(segstart(m)))-nv(l)) )./(Zp(Bseg(segstart(m)))-Zp(Bseg(segstart(m)-1)));
	 end;
	 if segend(m)==length(Bseg),  % Last point is above....
	   xend=[];
	   yend=[];
	 else;               % ...otherwise interpolate it.
	   xend=(Xp(Bseg(segend(m)+1))*(nv(l)-Zp(Bseg(segend(m)    ))) + Xp(Bseg(segend(m)    ))*(Zp(Bseg(segend(m)+1))-nv(l)) )./(Zp(Bseg(segend(m)+1))-Zp(Bseg(segend(m)    ))) ;
	   yend=(Yp(Bseg(segend(m)+1))*(nv(l)-Zp(Bseg(segend(m)    ))) + Yp(Bseg(segend(m)    ))*(Zp(Bseg(segend(m)+1))-nv(l)) )./(Zp(Bseg(segend(m)+1))-Zp(Bseg(segend(m)    )));
	 end;
         xseg=[ xstart , Xp(Bseg(segstart(m):segend(m)))' , xend ];
         yseg=[ ystart , Yp(Bseg(segstart(m):segend(m)))' , yend ];
	 lvlBS(:,jBS+[0:length(xseg)])=[ nv(l) xseg ; NaN yseg ];
         jBS=jBS+1+length(xseg);
        end;
      end;
      lvlBS(:,jBS+1:end)=[];

      ii=find(isnan(lvlBS(2,:)));

      % Now, apply the line joing algorithm *again* to put the different curves together
      % into a bunch of closed curves.

      % Begin and end of all line segments
      x2=[lvlBS(1,ii(1:end-1)+1) ; lvlBS(1,ii(2:end)-1) ];
      y2=[lvlBS(2,ii(1:end-1)+1) ; lvlBS(2,ii(2:end)-1) ];

      iZ=ones(1,size(x2,2));
      iN=1;
      iZ(iN)=0;
      seg=iN;
      rev1=1;rev2=2;  % keeps track of which direction to add points
       while any(iN),

	% Any line segments with the same endpoint that aren't used?
	% Use == for real numbers because the end points are calculated using
	% the same formula and hence should be exactly the same
	iN=find( x2(rev1,:)==x2(rev2,iN) & y2(rev1,:)==y2(rev2,iN) & iZ,1);

	if any(iN),  % If yes...
	  seg=[seg iN];  % Add it to the segment list...
	  iZ(iN)=0;      % ...and mark it used
	else         % If none...
	  if rev1==1;  % If we haven't searched backwards yet
            iN=seg(1);  % Take the other end of the lines
            seg=seg(end:-1:1); % Reverse the segment list
            rev1=2;rev2=1;
	  else         % we *have* searched backwards, and it really is the end
	    xseg=[];yseg=[];
	    for k=1:length(seg),
	      xseg=[lvlBS(1,ii(seg(k))+1:ii(seg(k)+1)-1) xseg];
	      yseg=[lvlBS(2,ii(seg(k))+1:ii(seg(k)+1)-1) yseg];
	    end;

	    fCS(:,iCS)=[nv(l);length(xseg)];
	    fCS(:,iCS+[1:length(xseg)])=[xseg ;yseg];

            % need these stats later
	    ncurves = ncurves + 1;
            levs(ncurves)=nv(l);
            I(ncurves) = iCS;
            Area(ncurves)=sum( diff(xseg).*(yseg(1:end-1)+yseg(2:end))/2 );

	    iCS=iCS+1+length(xseg);

            iN=find(iZ,1); % ...and go find the next unused segment
	    seg=iN;
            iZ(iN)=0;
            rev1=1;rev2=2;   % ...and search forward.
	  end;
	end;
      end;

end;
fCS(:,iCS:end)=[];

% OK, now we have all the closed curves formed.

% Plot patches in order of decreasing size. This makes sure that
% all the levels get drawn, not matter if we are going up a hill or
% down into a hole. When going down we shift levels though, you can
% tell whether we are going up or down by checking the sign of the
% area (since curves are oriented so that the high side is always
% the same side). Lowest curve is largest and encloses higher data
% always.

% The areas of 'interior holes' will be the same at all
% levels, as will any levels that go all around the outside (which will
% happen if there are interior valleys, which causes a problem (this
% didn't happen in contourf because boundary-finding was
% done differently) so we use sortrows to make sure the highest level
% gets drawn on top for positive areas (i.e. hills), but the lowest
% level gets drawn on top for NEGATIVE areas (holes).

%[FA,IA]=sort(-abs(Area));
[FA,IA]=sortrows([-abs(Area)' (sign(Area).*levs)']);

% below here code is basically identical to contourf

if ~ishold(cax),
    view(cax,2);
    set(cax,'Box','on','Layer','top');
    grid(cax,'off')
end

fig = ancestor(cax,'figure');
H=[];

% This is the colour for holes
if ~ischar(get(cax,'color'))
  bg = get(cax,'color');
else
  bg = get(fig,'color');
end

if isempty(col)
  edgec = get(fig,'defaultsurfaceedgecolor');
else
  edgec = col;
end
if isempty(lin)
  edgestyle = get(fig,'defaultpatchlinestyle');
else
  edgestyle = lin;
end

for jj=IA',
  nl=fCS(2,I(jj));
  lev=fCS(1,I(jj));
  if (lev ~= minz || draw_min ),
    xp=fCS(1,I(jj)+(1:nl));
    yp=fCS(2,I(jj)+(1:nl));
    clev = lev;           % color for filled region above this level
    if (sign(Area(jj)) ~=sign(Area(IA(1))) ),
      kk=find(nv==lev);
      kk0 = 1 + sum(nv<=minz) * (~draw_min);
      if (kk > kk0)
        clev=nv(kk-1);    % in valley, use color for lower level
      elseif (kk == kk0)
        clev=NaN;
      else
        clev=NaN;         % missing data section
        lev=NaN;
      end
    end

    if (isfinite(clev)),
      H=[H;patch(xp,yp,clev,'facecolor','flat','edgecolor',edgec, ...
              'linestyle',edgestyle,'userdata',lev,'parent',cax)];
    else
      H=[H;patch(xp,yp,clev,'facecolor',bg,'edgecolor',edgec, ...
              'linestyle',edgestyle,'userdata',fCS(1,I(jj)),'parent',cax)];
    end
  end;
end;

% Contourf strips out bnoundary points but I can't be bothered - if you
% want labelled contours you should really follow a tricontf call with
% a tricont call.

varargout={fCS,H};

%-----------------------------------------------------------------
function BS=findboundary(M);
% Finds the curves that make up the
% boundary of a triangulation
%
% First, examining ALL 2 point line segments in M
% to see if they are contained in one or
% two triangles in the mesh! This is slow.
%
lenM=size(M,1);

%Bseg=NaN(2,lenM);iB=1;
%
%pat=[1 2;  % Edge indices in order when triangles
%     2 3;  % are arranged CW
%     3 1];
%for k=1:lenM,
% for l=1:3,
%    [id,jid]  =find( (      M==M(k,pat(l,1))) );
%    [id2,jid2]=find( (M(id,:)==M(k,pat(l,2))) );
%    if length(id2)==1,
%      Bseg(:,iB)=M(k,pat(l,:))';
%      iB=iB+1;
%    end;
% end;
%end;
%Bseg(:,iB:end)=[];

% Try another way - this works MUCH MUCH faster
% (can be many orders of magnitude faster for large
% problems!), but will fail if the mesh is screwy
% (i.e. edges do not appear ONLY one or two times
% in the whole list of triangles()

% Put all the triangle edges in a list
alledge=[ M(:,[1 2]) ; M(:,[2 3]) ; M(:,[3 1]) ];

% put rows in increasing numerical order, then sort
% so the rows will either be unique edges, or pairs
% of edges
[sortedge,I]=sortrows(sort(alledge,2));
% Now identify the boundaries between pairs and
% unique edges - the diff will == 0 between
% pairs of identical edges
chges=[1;any(diff(sortedge)~=0,2);1];
% so find places where we have 1s one after another,
% since this marks a transition through a unique
% edge
ibdy=find(diff(chges)==0);
% ...and now get the original rows, before the
% numerical sorting
Bseg=alledge(I(ibdy),:)'; % back to original list

% Now I join all the individual boundary segments
% into a smaller number of actual curves using the
% follow foward/follow backward algorithm, all curves
% should traverse with the interior to their right.

BS=NaN(1,length(Bseg)+50);iBS=1;  % preallocate
iZ=ones(1,size(Bseg,2));  % flag to tell me if a line segment has not been added
                        % to a line
iN=1;
iZ(iN)=0;
seg=iN;
rev1=1;rev2=2;  % keeps track of which direction to add points
                % (differs if I am going forward or backward from a point)
while any(iN),

  % Any line segments with the same endpoint that aren't used?
  iN=find( Bseg(rev1,:)==Bseg(rev2,iN) & iZ ,1);

  if any(iN),  % If yes...
    seg=[seg iN];  % Add it to the segment list...
    iZ(iN)=0;      % ...and mark it used
  else         % If none...
    if rev1==1;  % If we haven't searched backwards yet
      iN=seg(1);  % Take the other end of the lines
    %  iZ(iN)=0;
      seg=seg(end:-1:1); % Reverse the segment list
      rev1=2;rev2=1;
    else         % we *have* searched backwards, and it really is the end
      xseg=[Bseg(2,seg),Bseg(1,seg(end))];  % make up the line

      BS(iBS+[0:length(xseg)])=[NaN xseg(end:-1:1) ];
      iBS=iBS+1+length(xseg);

      iN=find(iZ,1); % ...and go find the next unused segment
      seg=iN;
      iZ(iN)=0;
      rev1=1;rev2=2;   % ...and search forward.
    end;
  end;
end;
BS(iBS+1:end)=[];

author: Sable

repo: mcbench-benchmarks


text:
function [f0v,vrv,dfv,nf,aav]=fixpF0VexMltpBG4(x,fs,f0floor,nvc,nvo,mu,imgi,shiftm,smp,minm,pc,nc)

%	Fixed point analysis to extract F0
%	[f0v,vrv,dfv,nf]=fixpF0VexMltpBG4(x,fs,f0floor,nvc,nvo,mu,imgi,shiftm,smp,minm,pc,nc)
%	x	: input signal
%	fs	: sampling frequency (Hz)
%	f0floor	: lowest frequency for F0 search
%	nvc	: total number of filter channels
%	nvo	: number of channels per octave
%	mu	: temporal stretching factor
%	imgi	: image display indicator (1: display image, default)
%	shiftm	: frame shift in ms
%	smp	: smoothing length relative to fc (ratio)
%	minm	: minimum smoothing length (ms)
%	pc	: exponent to represent nonlinear summation
%	nc	: number of harmonic component to use (1,2,3)

%	Designed and coded by Hideki Kawahara
%	28/March/1999
%	04/April/1999 revised to multi component version
%	07/April/1999 bi-reciprocal smoothing for multi component compensation
%	01/May/1999 first derivative of Amplitude is taken into account
%	17/Dec./2000 display bug fix
%	19/Sep./2002 bug fix (mu information was discarded.)
%   07/Dec./2002 waitbar was added
%	30/April/2005 modification for Matlab v7.0 compatibility
%	10/Aug./2005 modified by Takahashi on waitbar
%	10/Sept./2005 mofidied by Kawahara on waitbar
%   11/Sept./2005 fixed waitbar problem

%f0floor=40;
%nvo=12;
%nvc=52;
%mu=1.1;

x=cleaninglownoise(x,fs,f0floor);

fxx=f0floor*2.0.^((0:nvc-1)/nvo)';
fxh=max(fxx);

dn=max(1,floor(fs/(fxh*6.3)));

if nc>2
	pm3=multanalytFineCSPB(decimate(x,dn),fs/dn,f0floor,nvc,nvo,mu,3,imgi); % error crrect 2002.9.19 (mu was fixed 1.1)
	pif3=zwvlt2ifq(pm3,fs/dn);
	[~,mm]=size(pif3);
	pif3=pif3(:,1:3:mm);
	pm3=pm3(:,1:3:mm);
end;

if nc>1
	pm2=multanalytFineCSPB(decimate(x,dn),fs/dn,f0floor,nvc,nvo,mu,2,imgi);% error crrect 2002.9.19(mu was fixed 1.1)
	pif2=zwvlt2ifq(pm2,fs/dn);
	[~,mm]=size(pif2);
	pif2=pif2(:,1:3:mm);
	pm2=pm2(:,1:3:mm);
end;

pm1=multanalytFineCSPB(decimate(x,dn*3),fs/(dn*3),f0floor,nvc,nvo,mu,1,imgi);% error crrect 2002.9.19(mu was fixed 1.1)
%%%% safe guard added on 15/Jan./2003
mxpm1=max(max(abs(pm1)));
eeps=mxpm1/10000000;
pm1(pm1==0)=pm1(pm1==0)+eeps;
%%%% safe guard end
pif1=zwvlt2ifq(pm1,fs/(dn*3));
%keyboard;

[~,mm1]=size(pif1);
mm=mm1;
if nc>1
	[~,mm2]=size(pif2);
	mm=min(mm1,mm2);
end;

if nc>2
	[~,mm3]=size(pif3);
	mm=min([mm1 mm2 mm3]);
end;

if nc == 2
	for ii=1:mm
		pif2(:,ii)=(pif1(:,ii).*(abs(pm1(:,ii))).^pc ...
			+pif2(:,ii)/2.*(abs(pm2(:,ii))).^pc )...
			./((abs(pm1(:,ii))).^pc+(abs(pm2(:,ii))).^pc);
	end;
end;
if nc == 3
	for ii=1:mm
		pif2(:,ii)=(pif1(:,ii).*(abs(pm1(:,ii))).^pc ...
			+pif2(:,ii)/2.*(abs(pm2(:,ii))).^pc ...
			+pif3(:,ii)/3.*(abs(pm3(:,ii))).^pc )...
			./((abs(pm1(:,ii))).^pc+(abs(pm2(:,ii))).^pc+(abs(pm3(:,ii))).^pc);
	end;
end;
if nc == 1
	pif2=pif1;
end;

%pif2=zwvlt2ifq(pm,fs/dn)*2*pi;
pif2=pif2*2*pi;
dn=dn*3;

[slp,~]=zifq2gpm2(pif2,f0floor,nvo);
[nn,mm]=size(pif2);
dpif=(pif2(:,2:mm)-pif2(:,1:mm-1))*fs/dn;
dpif(:,mm)=dpif(:,mm-1);
[dslp,~]=zifq2gpm2(dpif,f0floor,nvo);

damp=(abs(pm1(:,2:mm))-abs(pm1(:,1:mm-1)))*fs/dn;
damp(:,mm)=damp(:,mm-1);
damp=damp./abs(pm1);

%[c1,c2]=znormwght(1000);
fxx=f0floor*2.0.^((0:nn-1)/nvo)'*2*pi;
mmp=0*dslp;
[c1,c2b]=znrmlcf2(1);
if imgi==1; hpg=waitbar(0,'P/N map calculation'); end; % 07/Dec./2002 by H.K.%10/Aug./2005
for ii=1:nn
%	[c1,c2]=znrmlcf2(fxx(ii)/2/pi); % This is OK, but the next Eq is much faster.
	c2=c2b*(fxx(ii)/2/pi)^2;
	cff=damp(ii,:)/fxx(ii)*2*pi*0;
	mmp(ii,:)=(dslp(ii,:)./(1+cff.^2)/sqrt(c2)).^2+(slp(ii,:)./sqrt(1+cff.^2)/sqrt(c1)).^2;
    if imgi==1; waitbar(ii/nn); end; %,hpg); % 07/Dec./2002 by H.K.%10/Aug./2005
end;
if imgi==1; close(hpg); end;%10/Aug./2005

if smp~=0
	smap=zsmoothmapB(mmp,fs/dn,f0floor,nvo,smp,minm,0.4);
else
	smap=mmp;
end;

fixpp=zeros(round(nn/3),mm);
fixvv=fixpp+100000000;
fixdf=fixpp+100000000;
fixav=fixpp+1000000000;
nf=zeros(1,mm);
if imgi==1; hpg=waitbar(0,'Fixed pints calculation'); end; % 07/Dec./2002 by H.K.%10/Aug./2005
for ii=1:mm
	[ff,vv,df,aa]=zfixpfreq3(fxx,pif2(:,ii),smap(:,ii),dpif(:,ii)/2/pi,pm1(:,ii));
	kk=length(ff);
	fixpp(1:kk,ii)=ff;
	fixvv(1:kk,ii)=vv;
	fixdf(1:kk,ii)=df;
	fixav(1:kk,ii)=aa;
	nf(ii)=kk;
    if imgi==1 && rem(ii,10)==0; waitbar(ii/mm); end;% 07/Dec./2002 by H.K.%10/Aug./2005
end;
if imgi==1; close(hpg); end; % 07/Dec./2002 by H.K.%10/Aug./2005
fixpp(fixpp==0)=fixpp(fixpp==0)+1000000;

%keyboard
%[vvm,ivv]=min(fixvv);
%
%for ii=1:mm
%	ff00(ii)=fixpp(ivv(ii),ii);
%	esgm(ii)=fixvv(ivv(ii),ii);
%end;
np=max(nf);
f0v=fixpp(1:np,round(1:shiftm/dn*fs/1000:mm))/2/pi;
vrv=fixvv(1:np,round(1:shiftm/dn*fs/1000:mm));
dfv=fixdf(1:np,round(1:shiftm/dn*fs/1000:mm));
aav=fixav(1:np,round(1:shiftm/dn*fs/1000:mm));
nf=nf(round(1:shiftm/dn*fs/1000:mm));

if imgi==1
	cnmap(fixpp,smap,fs,dn,nvo,f0floor,shiftm);
end;
%ff00=ff00(round(1:shiftm/dn*fs/1000:mm));
%esgm=sqrt(esgm(round(1:shiftm/dn*fs/1000:mm)));
%keyboard;

return;
%------------------------------------------------------------------
function okid=cnmap(fixpp,smap,fs,dn,nvo,f0floor,shiftm)

% 	This function had a bug in map axis.
%	17/Dec./2000 bug fix by Hideki Kawahara.

dt=dn/fs;
[nn,mm]=size(smap);
aa=figure;
set(aa,'PaperPosition',[0.3 0.25 8 10.9]);
set(aa,'Position',[30 130 520 680]);
subplot(211);
imagesc([0 (mm-1)*dt*1000],[1 nn],20*log10(smap(:,round(1:shiftm/dn*fs/1000:mm))));axis('xy')
hold on;
tx=((1:shiftm/dn*fs/1000:mm)-1)*dt*1000;
plot(tx,(nvo*log(fixpp(:,round(1:shiftm/dn*fs/1000:mm))/f0floor/2/pi)/log(2)+0.5)','ko');
plot(tx,(nvo*log(fixpp(:,round(1:shiftm/dn*fs/1000:mm))/f0floor/2/pi)/log(2)+0.5)','w.');
hold off
xlabel('time (ms)');
ylabel('channel #');
colormap(jet);

okid=1;
return;

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

%%function pm=zmultanalytFineCSPm(x,fs,f0floor,nvc,nvo,mu,mlt);

%       Dual waveleta analysis using cardinal spline manipulation
%               pm=multanalytFineCSP(x,fs,f0floor,nvc,nvo);
%       Input parameters
%
%               x       : input signal (2kHz sampling rate is sufficient.)
%               fs      : sampling frequency (Hz)
%               f0floor : lower bound for pitch search (60Hz suggested)
%               nvc     : number of total voices for wavelet analysis
%               nvo     : number of voices in an octave
%				mu		: temporal stretch factor
%       Outpur parameters
%               pm      : wavelet transform using iso-metric Gabor function
%
%       If you have any questions,  mailto:kawahara@hip.atr.co.jp
%
%       Copyright (c) ATR Human Information Processing Research Labs. 1996
%       Invented and coded by Hideki Kawahara
%       30/Oct./1996

%t0=1/f0floor;
%lmx=round(6*t0*fs*mu);
%wl=2^ceil(log(lmx)/log(2));
%x=x(:)';
%nx=length(x);
%tx=[x,zeros(1,wl)];
%gent=((1:wl)-wl/2)/fs;

%nvc=18;

%wd=zeros(nvc,wl);
%wd2=zeros(nvc,wl);
%ym=zeros(nvc,nx);
%pm=zeros(nvc,nx);
%mpv=1;
%mu=1.0;
%for ii=1:nvc
%  t=gent*mpv;
%  t=t(abs(t)<3.5*mu*t0);
%  wbias=round((length(t)-1)/2);
%  wd1=exp(-pi*(t/t0/mu).^2);%.*exp(i*2*pi*t/t0);
%  wd2=max(0,1-abs(t/t0/mu));
%  wd2=wd2(wd2>0);
%  wwd=conv(wd2,wd1);
%  wwd=wwd(abs(wwd)>0.0001);
%  wbias=round((length(wwd)-1)/2);
%  wwd=wwd.*exp(i*2*pi*mlt*t(round((1:length(wwd))-wbias+length(t)/2))/t0);
%  pmtmp1=fftfilt(wwd,tx);
%  pm(ii,:)=pmtmp1(wbias+1:wbias+nx)*sqrt(mpv);
%  mpv=mpv*(2.0^(1/nvo));
%  keyboard;
%end;
%[nn,mm]=size(pm);
%pm=pm(:,1:mlt:mm);

%----------------------------------------------------------------
function pif=zwvlt2ifq(pm,fs)
%	Wavelet to instantaneous frequency map
%	fqv=wvlt2ifq(pm,fs)

%	Coded by Hideki Kawahara
%	02/March/1999

[~,mm]=size(pm);
pm=pm./(abs(pm));
pif=abs(pm(:,:)-[pm(:,1),pm(:,1:mm-1)]);
pif=fs/pi*asin(pif/2);
pif(:,1)=pif(:,2);

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

function [slp,pbl]=zifq2gpm2(pif,f0floor,nvo)
%	Instantaneous frequency 2 geometric parameters
%	[slp,pbl]=ifq2gpm(pif,f0floor,nvo)
%	slp		: first order coefficient
%	pbl		: second order coefficient

%	Coded by Hideki Kawahara
%	02/March/1999

[nn,~]=size(pif);
fx=f0floor*2.0.^((0:nn-1)/nvo)*2*pi;

c=2.0^(1/nvo);
g=[1/c/c 1/c 1;1 1 1;c*c c 1];
h=inv(g);

%slp=pif(1:nn-2,:)*h(1,1)+pif(2:nn-1,:)*h(1,2)+pif(3:nn,:)*h(1,3);
slp=((pif(2:nn-1,:)-pif(1:nn-2,:))/(1-1/c) ...
    +(pif(3:nn,:)-pif(2:nn-1,:))/(c-1))/2;
slp=[slp(1,:);slp;slp(nn-2,:)];

pbl=pif(1:nn-2,:)*h(2,1)+pif(2:nn-1,:)*h(2,2)+pif(3:nn,:)*h(2,3);
pbl=[pbl(1,:);pbl;pbl(nn-2,:)];

for ii=1:nn
	slp(ii,:)=slp(ii,:)/fx(ii);
	pbl(ii,:)=pbl(ii,:)/fx(ii);
end;

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

%function [c1,c2]=znormwght(n)

%zz=0:1/n:3;
%hh=[diff(zGcBs(zz,0)) 0]*n;
%c1=sum((zz.*hh).^2)/n;
%c2=sum((2*pi*zz.^2.*hh).^2)/n;

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

function p=zGcBs(x,k)

tt=x+0.0000001;
p=tt.^k.*exp(-pi*tt.^2).*(sin(pi*tt+0.0001)./(pi*tt+0.0001)).^2;

%--------------------------------------------
function smap=zsmoothmapB(map,fs,f0floor,nvo,mu,mlim,pex)

[nvc,mm]=size(map);
%mu=0.4;
t0=1/f0floor;
lmx=round(6*t0*fs*mu);
wl=2^ceil(log(lmx)/log(2));
gent=((1:wl)-wl/2)/fs;

smap=map;
mpv=1;
zt=0*gent;
iiv=1:mm;
for ii=1:nvc
	t=gent*mpv; %t0*mu/mpv*1000
	t=t(abs(t)<3.5*mu*t0);
	wbias=round((length(t)-1)/2);
	wd1=exp(-pi*(t/(t0*(1-pex))/mu).^2);
	wd2=exp(-pi*(t/(t0*(1+pex))/mu).^2);
	wd1=wd1/sum(wd1);
	wd2=wd2/sum(wd2);
	tm=fftfilt(wd1,[map(ii,:) zt]);
	tm=fftfilt(wd2,[1.0./tm(iiv+wbias) zt]);
	smap(ii,:)=1.0./tm(iiv+wbias);
	if t0*mu/mpv*1000 > mlim
		mpv=mpv*(2.0^(1/nvo));
	end;
end;

%--------------------------------------------
%function [ff,vv,df]=zfixpfreq2(fxx,pif2,mmp,dfv)
%
%nn=length(fxx);
%iix=(1:nn)';
%cd1=pif2-fxx;
%cd2=[diff(cd1);cd1(nn)-cd1(nn-1)];
%cdd1=[cd1(2:nn);cd1(nn)];
%fp=(cd1.*cdd1<0).*(cd2<0);
%ixx=iix(fp>0);
%ff=pif2(ixx)+(pif2(ixx+1)-pif2(ixx)).*cd1(ixx)./(cd1(ixx)-cdd1(ixx));
%vv=mmp(ixx);
%vv=mmp(ixx)+(mmp(ixx+1)-mmp(ixx)).*(ff-fxx(ixx))./(fxx(ixx+1)-fxx(ixx));
%df=dfv(ixx)+(dfv(ixx+1)-dfv(ixx)).*(ff-fxx(ixx))./(fxx(ixx+1)-fxx(ixx));

%--------------------------------------------
function [ff,vv,df,aa]=zfixpfreq3(fxx,pif2,mmp,dfv,pm)

aav=abs(pm);
nn=length(fxx);
iix=(1:nn)';
cd1=pif2-fxx;
cd2=[diff(cd1);cd1(nn)-cd1(nn-1)];
cdd1=[cd1(2:nn);cd1(nn)];
fp=(cd1.*cdd1<0).*(cd2<0);
ixx=iix(fp>0);
ff=pif2(ixx)+(pif2(ixx+1)-pif2(ixx)).*cd1(ixx)./(cd1(ixx)-cdd1(ixx));
%vv=mmp(ixx);
vv=mmp(ixx)+(mmp(ixx+1)-mmp(ixx)).*(ff-fxx(ixx))./(fxx(ixx+1)-fxx(ixx));
df=dfv(ixx)+(dfv(ixx+1)-dfv(ixx)).*(ff-fxx(ixx))./(fxx(ixx+1)-fxx(ixx));
aa=aav(ixx)+(aav(ixx+1)-aav(ixx)).*(ff-fxx(ixx))./(fxx(ixx+1)-fxx(ixx));

%--------------------------------------------
function [c1,c2]=znrmlcf2(f)

n=100;
x=0:1/n:3;
g=zGcBs(x,0);
dg=[diff(g) 0]*n;
dgs=dg/2/pi/f;
xx=2*pi*f*x;
c1=sum((xx.*dgs).^2)/n*2;
c2=sum((xx.^2.*dgs).^2)/n*2;

%--------------------------------------------
function x=cleaninglownoise(x,fs,f0floor)

flm=50;
flp=round(fs*flm/1000);
nn=length(x);
wlp=fir1(flp*2,f0floor/(fs/2));
wlp(flp+1)=wlp(flp+1)-1;
wlp=-wlp;

tx=[x(:)' zeros(1,2*length(wlp))];
ttx=fftfilt(wlp,tx);
x=ttx((1:nn)+flp);

return;

author: HidekiKawahara


text:
classdef TimeDepthAnalysisWindow < AnalysisWindow
    % CUMRATEANALYSISWINDOW shows the cumulative event rate
    properties
    end
    methods
        function obj=TimeDepthAnalysisWindow(ax)
            obj@AnalysisWindow(ax);
        end

        function prepare_axes(obj)
            % prepare the cumulative event rate axes
            if isempty(obj.ax.Tag)
                obj.ax.Tag = 'dvTimeDepth';
            end
            obj.ax.YDir = 'reverse';
            obj.ax.Title.String='Depth through Time';
            obj.ax.XLabel.String='Time';
            obj.ax.YLabel.String='Depth';
        end

        function [x,y]=calculate(obj,catalog)
            % calculate the cumulative number of events (y) over time (x)
            x=sort(catalog.Date);
            y=catalog.Depth;
        end

    end
end

author: CelsoReyes


text:
function [psi] = bar2psi(bar)
% Convert pressure from bar of pounds per square inch.
% Chad Greene 2012
psi = bar*14.5038;


text:
function imageInGray_InSpatialDomain = IFFT_image(imageInGray_fft_amplitudeOnly, imageInGray_fft_phaseOnly)
%IFFT_IMAGE

% interleaved the 2 FFT2 results
image_amplitudeWith_phaseInFreqDomain = imageInGray_fft_amplitudeOnly.*imageInGray_fft_phaseOnly;
% take only the real parts
imageInGray_InSpatialDomain = real(ifft2(image_amplitudeWith_phaseInFreqDomain));

end


text:
function [inference] = tapas_vlinear_inference(inference, pars)
%%
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

if ~isfield(inference, 'estimate_method')
    inference.estimate_method = @tapas_mcmc_blocked_estimate;
end

if ~isfield(inference, 'initialize_states')
    inference.initialize_states = @tapas_linear_init_states;
end

if ~isfield(inference, 'initialize_state')
    inference.initialize_state = @tapas_vlinear_init_state;
end

if ~isfield(inference, 'sampling_methods')
    inference.sampling_methods = {
        @(d, m, i, s) tapas_mh_mc3_sample_node(d, m, i, s, 2), ...
        @(d, m, i, s) tapas_sampler_vlinear_gibbs_node(d, m, i, s, 3), ...
            @tapas_sampler_mc3};
end

if ~isfield(inference, 'metasampling_methods')
    inference.metasampling_methods = {@tapas_mcmc_meta_diagnostics, ...
        @tapas_mcmc_meta_adaptive};
end

if ~isfield(inference, 'get_stored_state')
    inference.get_stored_state = @tapas_vlinear_get_stored_state;
end

if ~isfield(inference, 'prepare_posterior')
    inference.prepare_posterior = @tapas_linear_prepare_posterior;
end

if ~isfield(inference, 'mh_sampler')
    inference.mh_sampler = cell(4, 1);
end

if ~isfield(inference.mh_sampler{2}, 'propose_sample')
    inference.mh_sampler{2}.propose_sample = ...
        @tapas_mh_mc3_propose_gaussian_sample;
    inference.mh_sampler{2}.ar_rule = ...
        @tapas_mh_mc3_arc;
end

if ~isfield(inference, 'niter')
    if isfield(pars, 'niter')
        inference.niter = pars.niter;
    else
        inference.niter = 5000;
    end
end

if ~isfield(inference, 'nburnin')
    if isfield(pars, 'nburnin')
        inference.nburnin = pars.nburnin;
    else
        inference.nburnin = 5000;
    end
end

if ~isfield(inference, 'mc3it')
    if isfield(pars, 'mc3it')
        inference.mc3it = pars.mc3it;
    else
        inference.mc3it = 10;
    end
end

if ~isfield(inference, 'thinning')
    if isfield(pars, 'thinning')
        inference.thinning = pars.thinning;
    else
        inference.thinning = 0;
    end
end

if ~isfield(inference, 'ndiag')
    if isfield(pars, 'ndiag')
        inference.ndiag = pars.ndiag;
    else
        inference.ndiag = 200;
    end
end

end

author: translationalneuromodeling

repo: tapas


text:
function [ gx, gy ] = g02_xy ( gn )

%*****************************************************************************80
%
%% G02_XY returns the grid points for grid 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer GN, the grid size.
%
%    Output, real GX(GN,1), GY(GN,1), the grid coordinates.
%
  gx = [ ...
    0.05,  0.00,  0.00,  0.00,  0.10, ...
    0.10,  0.15,  0.20,  0.25,  0.30, ...
    0.35,  0.50,  0.50,  0.55,  0.60, ...
    0.60,  0.60,  0.65,  0.70,  0.70, ...
    0.70,  0.75,  0.75,  0.75,  0.80, ...
    0.80,  0.85,  0.90,  0.90,  0.95, ...
    1.00,  1.00,  1.00 ]';
  gy = [ ...
    0.45,  0.50,  1.00,  0.00,  0.15, ...
    0.75,  0.30,  0.10,  0.20,  0.35, ...
    0.85,  0.00,  1.00,  0.95,  0.25, ...
    0.65,  0.85,  0.70,  0.20,  0.65, ...
    0.90,  0.10,  0.35,  0.85,  0.40, ...
    0.65,  0.25,  0.35,  0.80,  0.90, ...
    0.00,  0.50,  1.00 ]';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% this script is to provide complete reaction identifiers (#174) that will be used for Met Atlas web-portal.

% load reaction ids
rxnIDs = importTsvFile('reactions.tsv');

% replace underscore with `R`
rxnIDs.rxnMAID = strrep(rxnIDs.rxnMAID,'_','R');
indToFill = getNonEmptyList(rxnIDs.rxnMAID,false);

% prepare a list to fill
idsToAdd = cell(19999,1);
for i=1:19999
    idsToAdd{i} = strcat('MAR',sprintf('%05d', i));
end
idsToAdd = setdiff(idsToAdd, rxnIDs.rxnMAID);

% sequentially fill in blank MA ids
rxnIDs.rxnMAID(indToFill) = idsToAdd(1:length(indToFill));
if isequal(length(rxnIDs.rxnMAID), length(unique(rxnIDs.rxnMAID)))
    fprintf('the filling looks okay.\n');
end

% insert dash between MA and RDDDDD
rxnIDs.rxnMAID = strrep(rxnIDs.rxnMAID,'MA','MA-');

exportTsvFile(rxnIDs,'reactions.tsv');

author: SysBioChalmers

repo: Human-GEM


text:
function [maxCI minCI] = calculateConfidenceInterval(tau, sPeriod, rCount,alpha,stabDev)
%This function calculates max and min confidence intervals for the XDEV
%calculation. The method used is found in the NIST Handbook of Frequency
%Stability Analysis page 37 section 5.3.1. tau is the array of tau values
%used in the stability calculation. Sample Period is used to calculate the
%number of values used in a particular tau calculation. rCount is the
%number of total readings at the sample period. alpha is an array of values
%for determine power law noise type. stabDev is the array of XDEV
%calculations.

tCount = length(tau); %get total number of calculations
maxCI = zeros(1,tCount); %initiate arrays
minCI = zeros(1,tCount);

for i = 1:tCount %for number of calculations
    avgCount = tau(i) / sPeriod; %gets count for averaging
    vCount = rCount / avgCount; %get number of loop iterations needed
    vCount = floor(vCount); %convert loopCount to integer in case it is a non int

    %perform confidence interval calculation based on noise type
    if alpha(i) >= 1.5 %dominant noise type is white PM
        temp = .99*(stabDev(i)/sqrt(vCount));
    elseif alpha(i) < 1.5 && alpha(i) >= .5 %dominant noise type is flicker PM
        temp = .99*(stabDev(i)/sqrt(vCount));
    elseif alpha(i) < .5 && alpha(i) >= -.5 %dominant noise type is white FM
        temp = .87*(stabDev(i)/sqrt(vCount));
    elseif alpha(i) < -.5 && alpha(i) >= -1.5 %dominant nosie type is flicker FM
        temp = .77*(stabDev(i)/sqrt(vCount));
    else %dominant noise type is random walk FM
        temp = .75*(stabDev(i)/sqrt(vCount));
    end

    %set min and max confidence limits
    maxCI(i) = stabDev(i) + temp;
    minCI(i) = stabDev(i) - temp;

end


%------------------------------------------------------------------
% Programed by:
%   - Lucas Rath (lucasrm25@gmail.com, https://github.com/lucasrm25)
%   -
%   -

%   Generate fancy plots (fp)
%------------------------------------------------------------------

classdef fp
    properties
    end

    methods(Static)

        function savefig(fname,varargin)

            p = inputParser;

            defaultFormat  = 'epsc';
            expectedFormat = {'epsc','svg','png','jpg'};
            addParameter(p,'format',defaultFormat, @(x) any(validatestring(x,expectedFormat)));
            addParameter(p,'fighandle',gcf, @(x) isa(x,'matlab.ui.Figure'));
            addParameter(p,'folder',fullfile(pwd,'images'), @(x) isa(x,'char'));
            parse(p,varargin{:});

            if ~ exist(p.Results.folder,'dir'), mkdir(p.Results.folder); end
            % set(p.Results.fighandle.CurrentAxes,'LooseInset',p.Results.fighandle.CurrentAxes.TightInset)
            saveas(p.Results.fighandle, fullfile(p.Results.folder,fname), p.Results.format)
        end

        function fig = f()
            fig = figure('Color','white');
            hold on, grid on;
        end

        function latexcode = m2latex(matrix)
            if ~isa(matrix,'sym')
                matrix = sym(matrix);
            end
            latexcode = latex(vpa(simplify(matrix)));
            if numel(latexcode)>=6 && strcmp(latexcode(6),'(') && strcmp(latexcode(end),')')
                latexcode(6) = '[';
                latexcode(end) = ']';
            end
            clipboard('copy',latexcode);
        end

        function str = figpos()
            a = gcf;
            pos = a.Position;
            str = ['figure(''Color'',''white'',''Position'',[' num2str(pos) ']);'];
            clipboard('copy',str);
        end

        % varargin{1}: line opacity
        function cl = getColor(n, varargin)
            colrs = lines(max(n));
            if nargin >= 2
                opacity = varargin{1};
            else
                opacity = [];
            end
            cl = [colrs(n,:), repmat(opacity,numel(n),1)];
        end
    end
end


text:
function [ area, radin, side ] = polygon_outrad_data ( n, radout )

%*****************************************************************************80
%
%% POLYGON_OUTRAD_DATA determines polygonal data from its outer radius.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 September 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of sides of the polygon.
%    N must be at least 3.
%
%    Input, real RADOUT, the outer radius of the polygon, that is,
%    the radius of the smallest circle that can be described
%    around the polygon.
%
%    Output, real AREA, the area of the regular polygon.
%
%    Output, real RADIN, the inner radius of the polygon, that is,
%    the radius of the largest circle that can be inscribed
%    within the polygon.
%
%    Output, real SIDE, the length of one side of the polygon.
%
  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'POLYGON_OUTRAD_DATA - Fatal error!\n' );
    fprintf ( 1, '  Input value of N must be at least 3.\n' );
    fprintf ( 1, '  but your input value was N = %d\n', n );
    error ( 'POLYGON_OUTRAD_DATA - Fatal error!' );
  end

  angle = pi / n;
  area = 0.5 * n * radout * radout * sin ( 2.0 * angle );
  side = 2.0 * radout * sin ( angle );
  radin = 0.5 * side / tan ( angle );

  return
end


text:
clear,clc

addpath(pwd);
cd cvx;
addpath(genpath(pwd));
cd ..;

load('Ns=3.mat');

Ns = 3;

NRF = 3;

SNR_dB = -35:5:5;
SNR = 10.^(SNR_dB./10);
realization = size(H,3);
smax = length(SNR);% enable the parallel

parfor reali = 1:realization
    [ FRF, FBB ] = SDR_AltMin( Fopt(:,:,reali), NRF);
    [ WRF, WBB ] = Receiver( Wopt(:,:,reali), NRF);
    for s = 1:smax
        R(s,reali) = log2(det(eye(Ns) + SNR(s)/Ns * pinv(WRF * WBB) * H(:,:,reali) * FRF * FBB * FBB' * FRF' * H(:,:,reali)' * WRF * WBB));
    end
end
plot(SNR_dB,sum(R,2)/realization,'Marker','diamond','LineWidth',1.5,'Color',[0.87058824300766 0.490196079015732 0]);
grid on
hold on

author: yuxianghao


text:
classdef IntersectionCoordComputer < TotalCoordinatesCalculator

    properties (Access = public)
        c
        theta
        nodes
        vertCoord
        boundCoord
        div
        totalCoord
    end

    properties (Access = private)
        vA
        vB
    end

    methods (Access = public)

        function obj = IntersectionCoordComputer(cParams)
            obj.init(cParams);
            obj.initBoundary();
            obj.obtainPrincipalVectors();
            obj.computeIntersections();
        end

    end

    methods (Access = private)

        function obtainPrincipalVectors(obj)
            A = obj.vertCoord(1,:);
            B = obj.vertCoord(2,:);
            C = obj.vertCoord(3,:);
            obj.vA = IntersectionCoordComputer.computeUnitaryVector(A,B);
            obj.vB = IntersectionCoordComputer.computeUnitaryVector(B,C);
        end

        function computeIntersections(obj)
            nodesX = obj.div(1)-1;
            nodesY = obj.div(2)-1;
            intNode = obj.nodes.bound +1;
            for jNodes = 1:nodesY
                pB = obj.boundCoord(obj.nodes.bound+1-jNodes,:);
                for iNodes = 1:nodesX
                    pA = obj.boundCoord(obj.nodes.vert+iNodes,:);
                    if obj.vA(1) == 0
                        [x,y] = IntersectionCoordComputer.computeVerticalIntersection(pA,pB,obj.vB);
                    elseif obj.vB(1) == 0
                        [x,y] = IntersectionCoordComputer.computeVerticalIntersection(pB,pA,obj.vA);
                    else
                        [x,y] = IntersectionCoordComputer.computeGeneralIntersection(pA,pB,obj.vA,obj.vB);
                    end
                    obj.totalCoord(intNode,:) = [x y];
                    intNode = intNode+1;
                end
            end
        end

    end

    methods (Static)

        function v_norm = computeUnitaryVector(A,B)
            v = B-A;
            m = norm(v);
            v_norm = v/m;
        end

        function [x,y] = computeVerticalIntersection(p1,p2,v)
            x = p2(1);
            y = (x-p1(1))*v(2)/v(1)+p1(2);
        end

        function [x,y] = computeGeneralIntersection(p1,p2,v1,v2)
            x = (p2(2)-p1(2)+p1(1)*v2(2)/v2(1)-p2(1)*v1(2)/v1(1))/(v2(2)/v2(1)-v1(2)/v1(1));
            y = (x-p2(1))*v1(2)/v1(1)+p2(2);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [name mfcc] = readKaldiFeature(fileName)

endian = 'l';

% FID = fopen('D:\Data\MagorBN\MFCC_LDA_MLLT_MVN\Ark\cml_mfcc_train.10.ark', 'r');
FID = fopen(fileName, 'r');
name = {};
mfcc = {};
byte_cnt = 1;
while 1
    [tmp,byte_cnt] = readUttName(FID, byte_cnt);
    if length(tmp)==0
        break;
    end
    name{end+1} = tmp;
    % fprintf('name{%d} = %s\n', length(name), name{end});
    header = readHeader(FID, byte_cnt, endian);
    data = fread(FID, header.nframe * header.dim, 'float32', 0, endian)';
    mfcc{end+1} = reshape(data, header.dim, header.nframe);
    if 0
        imagesc(mfcc{end});
        title(regexprep(name{end}, '_', '\\_'));
        colorbar;
        pause
    end
end
fclose(FID);

%%
function [header,byte_cnt] = readHeader(FID,byte_cnt, endian)
header.format = [];
byte_cnt = skipGap(FID, byte_cnt);
while 1
    tmp = fread(FID, 1, 'char');
    byte_cnt = byte_cnt + 1;
    if tmp==' '
        break;
    else
        header.format(end+1) = tmp;
    end
end
header.format = char(header.format);
byte_cnt = skipGap(FID,byte_cnt);
header.nframe = fread(FID, 1, 'int32', 0, endian); byte_cnt = byte_cnt + 4;
byte_cnt = skipGap(FID,byte_cnt);
header.dim = fread(FID, 1, 'int32', 0, endian); byte_cnt = byte_cnt + 4;

function byte_cnt = skipGap(FID, byte_cnt)
while 1
    tmp = fread(FID, 1, 'int8');
    byte_cnt = byte_cnt + 1;
    if tmp==0 || tmp==4
        break;
    end
end

function [uttName,byte_cnt] = readUttName(FID,byte_cnt)
uttName = [];
while 1
    tmp = fread(FID, 1, 'char');
    if length(tmp) == 0 % read end of the file
        break;
    elseif tmp==' '
        break;
    else
        uttName(end+1) = tmp;
    end
    byte_cnt = byte_cnt + 1;
end
uttName = char(uttName);

author: singaxiong

repo: SignalGraph


text:
% moveaxes() - move, resize, or copy Matlab axes using the mouse
%
% Usage: >> moveaxes
%        >> moveaxes(fig)
%        >> moveaxes off
%
% Note: clicking the left mouse button selects an axis
%       dragging the left mouse button resizes a selected axis
%       dragging the right mouse button copies a selected axis
%       clicking the middle mouse button deselects a selected axis
%
% Authors: Colin Humphries & Scott Makeig, CNL / Salk Institute, La Jolla

% Copyright (C) Colin Humphries & Scott Makeig, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% 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 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

% 01-25-02 reformated help & license -ad

function moveaxes(fig)

if nargin<1
  fig = gcf;
end

if ~isstr(fig)
   tmp = findobj('tag','moveaxes');
   if ~isempty(tmp)   % turn off previous moveaxes Tags
      ax = get(tmp,'children');
      set(ax,'ButtonDownFcn','','Selected','Off');
      set(tmp,'Tag','','UserData',[]);
   end

   ax = findobj('parent',fig,'type','axes');
   for a=1:length(ax)                    % make all axes visible
     axvis = get(ax(a),'visible');
     set(ax(a),'visible','on','UserData',axvis);
   end

   set(ax,'ButtonDownFcn','selectmoveresize');
   set(fig,'UserData',ax,'Tag','moveaxes');

elseif strcmp(fig,'off')
   fig=findobj('Tag','moveaxes');
   ax = get(fig,'UserData');
   for a=1:length(ax)                    % reset original axis visibility
     axvis= get(ax(a),'UserData')
     set(ax(a),'visible',axvis);
   end
   set(ax,'ButtonDownFcn','','Selected','off');
   set(fig,'Tag','','UserData',[]);
end

author: PatternRecognition

repo: OpenBMI


text:
% SVMLOO - Computes the exact leave-one-out estimate of the error rate
%          for a two-class SVM using Cauwenbergh's method.
%
% Syntax: [loo_est,conf_matrix,g_loo] = svmloo;
%         (evaluates the current SVM in memory)
%
%         [loo_est,conf_matrix,g_loo] = svmloo(X,y,a,b,g,ind,type,scale,Rs,Q)
%         (evaluates the given SVM)
%
%     loo_est: number of leave-one-out error errors
% conf_matrix: confusion matrix
%       g_loo: the resulting g for each example after the example is unlearned
%           X: matrix of training vectors stored columnwise
%           y: column vector of class labels (-1/+1) for training vectors
%           a: alpha coefficients
%           b: bias
%           g: partial derivatives of cost function w.r.t. alpha coefficients
%         ind: cell array containing indices of margin, error and reserve vectors
%         		ind{1}: indices of margin vectors
%         		ind{2}: indices of error vectors
%         		ind{3}: indices of reserve vectors
%        type: kernel type
%                1: linear kernel        X'*Y
%              2-4: polynomial kernel    (scale*X'*Y + 1)^type
%                5: Gaussian kernel with variance 1/(2*scale)
%       scale: kernel scale
%          Rs: inverse of extended kernel matrix for margin vectors
%           Q: extended kernel matrix for all vectors
%      g_flag: flag indicating whether or not to compute g_loo for the error vectors with g < -1
%
% Version 3.22e -- Comments to diehl@alumni.cmu.edu
%

function [loo_est,conf_matrix,g_loo] = svmloo(varargin)

% flags for example state
MARGIN    = 1;
ERROR     = 2;
RESERVE   = 3;
UNLEARNED = 4;

if (nargin == 0)

   % define global variables
   global a;
   global b;
   global g;
   global ind;
   global Q;
   global Rs;
   global scale;
   global type;
   global X;
   global y;

else

   % define arguments
   X = varargin{1};
   y = varargin{2};
   a = varargin{3};
   b = varargin{4};
   g = varargin{5};
   ind = varargin{6};
   type = varargin{7};
   scale = varargin{8};
   Rs = varargin{9};
   Q = varargin{10};
   g_flag = varargin(11);

end;

% if the user wants g_loo, make sure to compute g_loo for the error vectors with initial g < -1.
% if we only care about the error rate, we don't need to unlearn these examples because they are
% guaranteed to be classified incorrectly.
if (nargout == 3)
   g_flag = 1;
else
   g_flag = 0;
end;

% initialize variables
num_MVs = length(ind{MARGIN});      % number of margin vectors
loo_est = 0;                        % number of leave-one-out errors
a_orig = a;                         % original value of a
b_orig = b;                         % original value of b
Rs_orig = Rs;                       % original value of Rs
Q_orig = Q;                         % original value of Q
g_orig = g;                         % original value of g
num_MVs_orig = num_MVs;             % original value of num_MVs
ind_orig = ind;                     % original value of ind
g_loo = g;

% initialize confusion matrix
conf_matrix = zeros(2,3);
if (length(ind{RESERVE}) > 0)
	conf_matrix(1,1) = sum(y(ind{RESERVE}) == 1);
	conf_matrix(2,2) = sum(y(ind{RESERVE}) == -1);
end;

% begin leave-one-out estimation
ind_loo = [ind{MARGIN} ind{ERROR}];
num_tested = 1;
disp('Beginning LOO error rate estimation.');
for i = 1:length(ind_loo)

   % select example to unlearn
   indc = ind_loo(i);

   % unlearn example
   if ((g(indc) >= -1) | (g_flag))
      unlearn(indc);
      g_loo(indc) = g(indc);
   end;

   if (mod(num_tested,50) == 0)
      s = sprintf('Unlearned and tested %d examples.',num_tested);
      disp(s);
   end;
   num_tested = num_tested + 1;

   % check to see if the example is now misclassified and record results
   loo_est = loo_est + (g(indc) < -1);
   if (y(indc) == 1)
      j = 1 + (g(indc) <= -1) + (g(indc) == -1);
      conf_matrix(1,j) = conf_matrix(1,j) + 1;
   else
      j = 1 + (g(indc) >= -1) + (g(indc) == -1);
      conf_matrix(2,j) = conf_matrix(2,j) + 1;
   end;

   % reset to original state prior to unlearning example
   a = a_orig;
   b = b_orig;
   Rs = Rs_orig;
   Q = Q_orig;
   g = g_orig;
   ind = ind_orig;
   num_MVs = num_MVs_orig;

end;
if (mod(num_tested-1,50) ~= 0)
   s = sprintf('Unlearned and tested %d examples.',num_tested-1);
   disp(s);
end;
s = sprintf('Process complete!\n');
disp(s);

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function complextest_AD1()
% Test AD for a complex optimization problem on a product manifold (struct)

    % Verify that Manopt was indeed added to the Matlab path.
    if isempty(which('spherecomplexfactory'))
        error(['You should first add Manopt to the Matlab path.\n' ...
		       'Please run importmanopt.']);
    end

    % Verify that the deep learning tool box was installed
    assert(exist('dlarray', 'file') == 2, ['Deep learning tool box is '...
    'needed for automatic differentiation.\n Please install the'...
    'latest version of the deep learning tool box and \nupgrade to Matlab'...
    ' R2021b if possible.'])

    % Generate the problem data.
    n = 100;
    A = randn(n) + 1i*randn(n);
    A = .5*(A+A');

    % Create the product manifold
    S = spherecomplexfactory(n);
    manifold.x = S;
    manifold.y = S;
    problem.M = productmanifold(manifold);  %struct

    % For Matlab R2021b or later, define the problem cost function as usual
    % problem.cost  = @(X) -real(X.x'*A*X.y);

    % For Matlab R2021a or earlier, translate the cost function into a
    % particular format with the basic functions in /functions_AD
    problem.cost  = @(X) -creal(cprod(cprod(ctransp(X.x), A), X.y));

    % Define the gradient and the hessian via automatic differentiation
    problem = manoptAD(problem);

    % Numerically check gradient and Hessian consistency.
    figure;
    checkgradient(problem);
    figure;
    checkhessian(problem);

    % Solve.
    [x, xcost, info] = trustregions(problem);          %#ok<ASGLU>

    % Test
    ground_truth = svd(A);
    distance = abs(ground_truth(1) - (-problem.cost(x)));
    fprintf('The distance between the ground truth and the solution is %e \n',distance);

end

author: NicolasBoumal

repo: manopt


text:
function [yd3] = l2yd3(l)
% Convert volume from liters to cubic yards.
% Chad Greene 2012
yd3 = l*0.0013079506193;


text:
function [] = ControlPanel()

% initialize control variables
Lx = 200*10^-3; % billet size
Ly = 200*10^-3;
dx = 01*10^-3;
dy = 01*10^-3;
dt = 0.001;

% Temperature definitions
initial_Temperature = 510+273;
left_Temperature = 65+273;
upper_Temperature = 65+273;
right_Temperature = 65+273;
bottom_Temperature = 65+273;
BC_temperature = left_Temperature;
new_Temperature = BC_temperature;
reference_Temperature = 200+273;

% Secondary variables
boolean = 0;
timeRequired = 0;
time = 26;

% Material properties
alpha = 6.584*10^-5; % Steel
rho = 2700;
Cp = 900;

courant_Number = alpha*dt/(dx^2);

% Define computational domain (Geometry)

% storage parameters
x_intervals = Lx/dx + 1;
y_intervals = Ly/dy + 1;

% Matrices
T_old =  zeros(x_intervals,y_intervals); % 2-D storage for previous time step
T_new = zeros(x_intervals,y_intervals); % 2-D storage for current time step
%T = zeros((time/dt),x_intervals,y_intervals); % Mega Temperature matrix to stroe in 2-D for each time step, hence 3-D

% Initialize domain
T_old = InitializeSolution(T_old,initial_Temperature,left_Temperature,upper_Temperature,right_Temperature,bottom_Temperature,x_intervals,y_intervals);
T_boundary = zeros(1,time/dt);

% formulation : EXPLICIT scheme
for time_index = 1:1:(time/dt)

    T_new = ExplicitScheme(T_old,alpha,dx,dy,dt,x_intervals,y_intervals);
    % Introduce the water bath dynaic BC factor
    new_Temperature = CalculatedBathAdjustedBoundaryTemperature(T_new,T_old,x_intervals,y_intervals,dt,rho,Cp,new_Temperature);
    boolean = ExitConditionEvaluation(T_new,reference_Temperature,x_intervals,y_intervals);

    if(boolean == 1 && timeRequired == 0)
        timeRequired = time_index;
    end

    T_boundary(time_index) = new_Temperature;

    T_old = ReInitializeSolution(T_new,new_Temperature,x_intervals,y_intervals);

end

T_plot = T_new';
T_plotting = CalculatePlottingMatrix(T_plot,x_intervals,y_intervals); % Plotting matrix

%Contour plot on the 2-D spatial domain
T_constant = initial_Temperature:.001:max([initial_Temperature,left_Temperature,upper_Temperature,right_Temperature,bottom_Temperature]);
[C,h] = contour(T_plotting,T_constant);

end


text:
function [ y, next ] = p12_equil ( neqn, next )

%*****************************************************************************80
%
%% P12_EQUIL returns equilibrium solutions of problem p12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, integer NEXT, the index of the previous
%    equilibrium, which should be 0 on first call.
%
%    Output, real Y(NEQN), the "next" equilibrium solution, if any.
%
%    Output, integer NEXT, the index of the current equilibrium,
%    or 0 if there are no more.
%
  if ( next == 0 )
    next = 1;
    y(1:neqn,1) = 0.0;
  else
    next = 0;
    y = [];
  end

  return
end


text:
function another_colormap
%ANOTHER_COLORMAP try another color map
%
% Example:
%   another_colormap
% See also: testall

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

j = jet (128) ;
j = j (48:112, :) ;

% jj = linspace (0,1,64)' ./ sum (jet,2) ;
% j (:,1) = j (:,1) .* jj ;
% j (:,2) = j (:,2) .* jj ;
% j (:,3) = j (:,3) .* jj ;

% white = [1 1 1] ;
% gray = [.5 .5 .5] ;                                                       %#ok

% j = [white ; purple ; j ] ;
disp ('j = ') ;
disp (j)

image (1:size(j,1)) ;
colormap (j) ;


text:
function [fx,dfdx,dfdP] = f_Id(x,P,u,in)
% this function applies the identity mapping
fx = x;
dfdx = eye(size(x,1));
dfdP = zeros(size(P,1),size(x,1));

author: MBB-team

repo: VBA-toolbox


text:
function pout = propagate_bounds_from_arbitrary_quadratics(p)

pout = p;
if p.bilinears~=0
    F_struc = p.F_struc;

    p.F_struc = [-p.F_struc(1:p.K.f,:);p.F_struc];
    p.K.f=2*p.K.f;

    InequalityConstraintState = [p.EqualityConstraintState;p.EqualityConstraintState;p.InequalityConstraintState];

    if p.K.l+p.K.f>0
        quadratic_variables = find(p.bilinears(:,2) == p.bilinears(:,3));
        if ~isempty(quadratic_variables)
            quadratic_variables = p.bilinears(quadratic_variables,1);
            for i = 1:length(quadratic_variables)
                k = quadratic_variables(i);

                if (p.lb(k) < p.ub(k)-1e-4)
                    x = p.bilinears(p.bilinears(:,1)==k,2);% x^2
                    candidates = find((InequalityConstraintState==1) & p.F_struc(1:p.K.f+p.K.l,1+k))';
                    for j = candidates
                        a = p.F_struc(j,2:end);
                        b = p.F_struc(j,1);
                        aij = a(k);

                        % rest + cij*x + (+-)x^2 >= 0
                        cij = a(x);
                        cij = cij/abs(aij);
                        b = b/abs(aij);
                        a = a/abs(aij);
                        %  aij = sign(aij);

                        a(k) = 0;
                        a(x) = 0;
                        indNEG = find(a < 0);
                        indPOS = find(a > 0);
                        LB = p.lb;
                        UB = p.ub;
                        LB(k) = 0;
                        UB(k) = 0;
                        LB(x) = 0;
                        UB(x) = 0;
                        a(k) = 0;
                        a(x) = 0;

                        if aij < 0
                            % rest + cij*x - x^2 >= 0

                            % Derive upper bound on rest
                            rest = (b+a([indPOS(:);indNEG(:)])*[UB(indPOS);LB(indNEG)]);
                            % Write as rest >= (x - center)^2-(center)^2
                            center = cij/2;
                            radii2 = rest + (cij/2)^2;
                            if radii2 > 0
                                newUB = center + sqrt(radii2);
                                newLB = center - sqrt(radii2);
                                p.lb(x) = max(p.lb(x),newLB);
                                p.ub(x) = min(p.ub(x),newUB);
                            end
                        else
                            % rest + cij*x + x^2 >= 0
                            % (x+cij/2)^2 - (cij/2)^2 + rest >= 0
                            % (x+cij/2)^2 >= radii
                            % Lower bound on rest
                            rest = (b+a([indPOS(:);indNEG(:)])*[UB(indPOS);LB(indNEG)]);
                            center = -cij/2;
                            radii2 = -rest + (cij/2)^2;
                            if radii2 > 0
                                left  = center-sqrt(radii2);
                                right = center+sqrt(radii2);
                                if p.ub(x) < right
                                    p.ub(x) = min(p.ub(x),left);
                                end
                                if p.lb(x) > left
                                    p.lb(x) = max(p.lb(x),right);
                                end
                            end
                        end
                    end

                elseif p.lb(k)<0

                end
            end
        end
    end

    if p.K.l+p.K.f>0
        bilinear_variables = find(p.bilinears(:,2) ~= p.bilinears(:,3));
        if ~isempty(bilinear_variables)
            bilinear_variables = p.bilinears(bilinear_variables,1);
            for i = 1:length(bilinear_variables)
                k = bilinear_variables(i);
                if p.lb(k) >= -5000000000 & (p.lb(k) < p.ub(k)-1e-4)
                    x = p.bilinears(p.bilinears(:,1)==k,2);% x^2
                    y = p.bilinears(p.bilinears(:,1)==k,3);% x^2
                    candidates = find((InequalityConstraintState==1) & p.F_struc(1:p.K.f+p.K.l,1+k))';
                    for j = candidates
                        a = p.F_struc(j,2:end);
                        aij = a(k);
                        if aij > 0
                            indNEG = find(a < 0);
                            indPOS = find(a > 0);
                            LB = p.lb;
                            UB = p.ub;
                            LB(k) = 0;
                            UB(k) = 0;
                            a(k) = 0;
                            newLB = (-p.F_struc(j,1)-a([indPOS(:);indNEG(:)])*[UB(indPOS);LB(indNEG)])/aij;
                            p.lb(k) = max(p.lb(k),newLB);

                        elseif aij < 0
                            indNEG = find(a < 0);
                            indPOS = find(a > 0);
                            LB = p.lb;
                            UB = p.ub;
                            LB(k) = 0;
                            UB(k) = 0;
                            a(k) = 0;
                            newUB = (p.F_struc(j,1)+a([indPOS(:);indNEG(:)])*[UB(indPOS);LB(indNEG)])/(-aij);
                            p.ub(k) = min(p.ub(k),newUB);
                        end
                    end

                    if p.lb(k)>0 & p.lb(x)>0 & p.lb(y)>0
                        p.ub(x) = min(p.ub(x), p.ub(k)/p.lb(y));
                        p.ub(y) = min(p.ub(y), p.ub(k)/p.lb(x));
                        p.lb(x) = max(p.lb(x), p.lb(k)/p.ub(y));
                        p.lb(y) = max(p.lb(y), p.lb(k)/p.ub(x));
                    end

                elseif p.lb(k)<0

                end
            end
        end
    end
end

if ~isequal([p.lb p.ub],[pout.lb pout.ub])
    quad_v = find(p.bilinears(:,2) == p.bilinears(:,3));
    quad_x = p.bilinears(quad_v,2);
    quad_v = p.bilinears(quad_v,1);
    if ~isempty(quad_v)
        % y = x^2, x>=0, y >= L means x >= sqrt(L)
        k = find(p.lb(quad_x)>=0 & p.lb(quad_v)>0);
        if ~isempty(k)
            p.lb(quad_x(k)) = max(p.lb(quad_x(k)),sqrt(p.lb(quad_v(k))));
        end
    end
end

if ~isequal([p.lb p.ub],[pout.lb pout.ub])
    pout.changedbounds = 1;
end
pout.lb = p.lb;
pout.ub = p.ub;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [f,lambda,delta] = fit(ori,varargin)
% determines the fibre that fits best a list of orientations
%
% Syntax
%   f = fibre.fit(ori) % fit fibre to a list of orientations
%
%   f = fibre.fit(ori,'local')
%
%   f = fibre.fit(odf) % fit fibre to a list of orientations
%
% Input
%  ori1, ori2, ori - @orientation
%  odf - @SO3Fun
%
% Output
%  f       - @fibre
%  lambda  - eigenvalues of the orientation tensor
%
% Options
%  local   - fast approach for locally concentrated orientations

if isa(ori,'SO3Fun')

  odf = ori;
  cs = odf.CS.properGroup;

  hGrid = Miller(equispacedS2Grid('resolution',10*degree,...
    cs.fundamentalSector),cs);
  rGrid = equispacedS2Grid('resolution',1.25*degree,odf.SS.fundamentalSector);
  lGrid = equispacedS2Grid('resolution',0.5*degree,'maxTheta',10*degree);

  % search through all pole figures
  v = -inf;
  for ih = 1:length(hGrid)

    % detect maximum in pole figure
    [~,ir] = max(calcPDF(odf,hGrid(ih),rGrid));

    % detect local maximum in corresponding inverse pole figure
    lhGrid = rotation.byAxisAngle(zvector + hGrid(ih),pi) * lGrid;
    [~,ih2] = max(calcPDF(odf,lhGrid,rGrid(ir)));

    % detect local maximum in corresponding pole figure - second iteration
    lrGrid = rotation.byAxisAngle(zvector + rGrid(ir),pi) * lGrid;
    [vTmp,ir2] = max(calcPDF(odf,lhGrid(ih2),lrGrid));

    % maybe we have found a new optimum
    if vTmp > v
      f = fibre(Miller(lhGrid(ih2),odf.CS),lrGrid(ir2));
      v = vTmp;
    end

  end
  return

elseif check_option(varargin,{'noSymmetry','local'}) || ori.CS.numProper==1

  [~,~,lambda,eigv] = mean(ori,varargin{:});

  ori12 = orientation(quaternion(eigv(:,4:-1:3)),ori.CS,ori.SS);
  f = fibre(ori12(1),ori12(2),'full',varargin{:});

  if nargout == 3
    delta = norm(angle(f,ori)) / length(ori);
  end

  return

end

% the global approach
odf = calcDensity(ori,'halfwidth',10*degree);
cs = ori.CS.properGroup;
rotSym = cs.rot;
fs = cs.fundamentalSector;

hGrid = Miller(equispacedS2Grid('resolution',10*degree,fs),cs);
rGrid = equispacedS2Grid('resolution',1.25*degree,ori.SS.fundamentalSector);
% a local grid
lGrid = equispacedS2Grid('resolution',0.5*degree,'maxTheta',10*degree);
%localhGrid = [localhGrid;-localhGrid];

% search through all pole figures
v = inf;
for ih = 1:length(hGrid)

  % detect maximum in pole figure
  [~,ir] = max(calcPDF(odf,hGrid(ih),rGrid));

  % detect local maximum in corresponding inverse pole figure
  lhGrid = rotation.byAxisAngle(zvector + hGrid(ih),pi) * lGrid;
  [~,ih2] = max(calcPDF(odf,lhGrid,rGrid(ir)));

  % project to fibre
  d = angle_outer(times(rotSym,lhGrid(ih2),1), times(inv(ori),rGrid(ir),0),'noSymmetry').';
  [m,idSym] = min(d,[],2);

  % robust estimate
  oriP = times(ori(m<quantile(m,0.8)), rotSym(idSym(m<quantile(m,0.8))),0);
  %oriP = times(ori, rotSym(idSym),0);

  % fit fibre
  fTmp = fibre.fit(oriP,'noSymmetry');
  %fTmp = fibre(Miller(rhGrid(id2),ori.CS),rGrid(id));

  % maybe we have found a new optimum
  vTmp = mean(angle(fTmp,ori));
  if vTmp < v, f = fTmp; v = vTmp; end

end

function test

  cs = crystalSymmetry('432');
  f = fibre.rand(cs);
  odf = fibreODF(f,'halfwidth',20*degree);

  ori = discreteSample(odf,100000);

  tic
  fRec = fibre.fit(ori,'refine')
  toc

  mean(angle(fRec,ori))./degree

  mean(angle(f,ori))./degree

  angle(fRec,f)./degree

end

end

author: mtex-toolbox

repo: mtex

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


text:
function g = cell2mat(f)
%CELL2MAT   Convert an array of CHEBTECH objects into an array-valued CHEBTECH.
%   G = CELL2MAT(F) converts the CHEBTECH array F into a single array-valued
%   CHEBTECH G. F should be a vector array (i.e., not a matrix).
%
% Example:
%   f = chebtech.constructor(@sin);
%   g = chebtech.constructor(@cos);
%   h = cell2mat([f g]);
%
% See also MAT2CELL.
%
% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% TODO: This function is probably not needed anymore.

% Return an empty result:
if ( isempty(f) || numel(f) == 1 )
    g = f;
    return
end

% Extract the coeffs from each of the CHEBTECH objects:
coeffs = { f.coeffs };

% These should all have the same length...
lengths = cellfun(@(x) size(x, 1), coeffs);
minlength = min(lengths);
maxlength = max(lengths);

% If not, then prolong and get the new values:
if ( minlength ~= maxlength )
    for k = 1:numel(f)
        f(k) = prolong(f(k), maxlength);
    end
    coeffs = { f.coeffs };
end

% Append new data to an empty CHEBTECH:
g = f.make(); % Make an empty CHEBTECH.
g.ishappy = min([f.ishappy]);
g.coeffs = cell2mat(coeffs);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function imgOut = bilateralNoiseRemoval(img, sigma_s, sigma_r)
%
%       imgOut = bilateralNoiseRemoval(img, sigma_s, sigma_r)
%
%
%        Input:
%           -img: input image to be denoised
%           -simga_s: spatial sigma; size of the neighborhood to be used
%           for denosing
%           -sigma_r: range sigma; what intesities to use for filtering
%
%        Output:
%           -imgOut: filtered denoised image
%
%     Copyright (C) 2013  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

col = size(img,3);

if(~exist('sigma_s', 'var'))
    sigma_r = 4;
end

if(sigma_s<4)
    sigma_s = 4;
end

if(~exist('sigma_r', 'var'))
    sigma_r = 0.01;
end

if(sigma_r <= 0.0)
    sigma_r = 0.01;
end

switch col
	case 1
    	minC = min(img(:));
        maxC = max(img(:));
        delta = (maxC - minC);
        img = (img - minC) / delta;
        imgOut = bilateralFilter(img, [], 0.0, 1.0, sigma_s, sigma_r) * delta + minC;

    case 3
        imgLab = ConvertXYZtoCIELab(ConvertRGBtoXYZ(img, 0), 0);

        for i=1:col
            tmp = imgLab(:,:,i);
            minC = min(tmp(:));
            maxC = max(tmp(:));
            delta = (maxC - minC);
            tmp = (tmp - minC) / delta;
            imgLab(:,:,i) = bilateralFilter(tmp,[], 0.0, 1.0,sigma_s,sigma_r) * delta + minC;
        end

        imgOut = ConvertRGBtoXYZ(ConvertXYZtoCIELab(imgLab,1),1);

    otherwise
        imgOut = zeros(size(img));

        for i=1:col
            tmp = img(:,:,i);
            minC = min(tmp(:));
            maxC = max(tmp(:));
            delta = (maxC - minC);
            tmp = (tmp - minC) / delta;
            imgOut(:,:,i) = bilateralFilter(tmp, [], 0.0, 1.0,sigma_s, sigma_r) * delta + minC;
        end
end

end

author: banterle


text:
% Computes log(besseli(0,x)) robustly.  Computing it directly causes
% numerical problems at high x, but the function has asymptotic linear
% behaviour, which we approximate here for high x.
%
% author: Daniel C Alexander (d.alexander@ucl.ac.uk)
%

function lb0 = logbesseli0(x)

% For very large arguments to besseli, we approximate log besseli using a
% linear model of the asymptotic behaviour.
% The linear parameters come from this command:
% app=regress(log(besseli(0,500:700))',[ones(201,1) (500:700)']);
app = [-3.61178295877576 0.99916157999904];

lb0 = zeros(length(x), 1);
exact = find(x<700);
approx = find(x>=700);
lb0(exact) = log(besseli(0, x(exact)));
%lb0(approx) = x(approx)*app(2) + app(1);

% This is a more standard approximation.  For large x, I_0(x) -> exp(x)/sqrt(2 pi x).
lb0(approx) = x(approx) - log(2*pi*x(approx))/2;

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [Population,FrontNo,d2] = EnvironmentalSelection(Population,RPSet,N)
% The environmental selection of RPD-NSGA-II

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Normalization
    PopObj = Population.objs;
    Fmin   = min(PopObj,[],1);
    Fmax   = max(PopObj,[],1);
    PopObj = (PopObj-repmat(Fmin,size(PopObj,1),1))./repmat(Fmax-Fmin,size(PopObj,1),1);

    %% Association
    normP   = sqrt(sum(PopObj.^2,2));
    Cosine  = 1 - pdist2(PopObj,RPSet,'cosine');
    d1      = repmat(normP,1,size(RPSet,1)).*Cosine;
    d2      = repmat(normP,1,size(RPSet,1)).*sqrt(1-Cosine.^2);
    [d2,RP] = min(d2,[],2);
    d1      = d1((1:length(RP))'+(RP-1)*length(RP));

    %% Favor extreme solutions
    ND              = find(NDSort(PopObj,1)==1);
    [~,Extreme]     = max(PopObj(ND,:),[],1);
    d1(ND(Extreme)) = 0;
    d2(ND(Extreme)) = 0;

    %% Non-RPD-dominated sorting
    [FrontNo,MaxFNo] = NRPDDSort(PopObj,d1,d2,RP,N);
    Next = FrontNo < MaxFNo;

    %% Select the solutions in the last front
    Last     = find(FrontNo==MaxFNo);
    [~,Rank] = sort(d2(Last));
    Next(Last(Rank(1:N-sum(Next)))) = true;

    %% Population for next generation
    Population = Population(Next);
    FrontNo    = FrontNo(Next);
    d2         = d2(Next);
end

function [FrontNo,MaxFNo] = NRPDDSort(PopObj,d1,d2,RP,nSort)
% Non-RPD-dominated sorting (based on deductive sort)

    [N,M]   = size(PopObj);
    FrontNo = inf(1,N);
    MaxFNo  = 0;
    while sum(FrontNo<inf) < min(nSort,N)
        MaxFNo = MaxFNo + 1;
        Sorted = FrontNo ~= inf;
        D      = Sorted;
        for i = 1 : N
            if ~D(i)
                for j = i+1 : N
                    if ~D(j)
                        domi = 0;
                        for m = 1 : M
                            if PopObj(i,m) < PopObj(j,m)
                                if domi == -1
                                    domi = 0;
                                    break;
                                else
                                    domi = 1;
                                end
                            elseif PopObj(i,m) > PopObj(j,m)
                                if domi == 1
                                    domi = 0;
                                    break;
                                else
                                    domi = -1;
                                end
                            end
                        end
                        % The definition of dominance is modified here
                        % since the one in the original paper is
                        % questionable
                        if domi == 0 && RP(i)==RP(j)
                            if d1(i)+5*d2(i) < d1(j)+5*d2(j)
                            	domi = 1;
                            elseif d1(i)+5*d2(i) > d1(j)+5*d2(j)
                            	domi = -1;
                            end
                        end
                        if domi == 1
                            D(j) = true;
                        elseif domi == -1
                            D(i) = true;
                            break;
                        end
                    end
                end
                if ~D(i)
                    FrontNo(i) = MaxFNo;
                end
            end
        end
    end
end


text:
function f = ivmKernelObjective(params, model)

% IVMKERNELOBJECTIVE Compute the negative of the IVM log likelihood approximation.
% FORMAT

% DESC computes the IVM negative log likelihood approximation at a
% given set of kernel parameters. Adiditonally if there is any
% regularisation present on the kernel parameter it computes the prior
% probability and adds it in.
% ARG params : the parameter values where the obective is to be
% evaluated.
% ARG model : the model structure for which the objective is to be
% evaluated.
%
% SEEALSO : ivmKernelGradient, ivmOptimiseKernel, kernExpandParam,
% ivmApproxLogLikelihood, kernPriorLogProb
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% IVM

%/~
if any(isnan(params))
  warning('Parameter is NaN')
end
%~/

model.kern = kernExpandParam(model.kern, params);
f = ivmApproxLogLikelihood(model);
f = f + kernPriorLogProb(model.kern);
f = -f;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [fc_requested, fc_programmed, fs_requested, fs_programmed] = read_header_from_bin_fid(fid)

fc_requested = inf;
fc_programmed = inf;
fs_requested = inf;
fs_programmed = inf;

% fid = fopen(bin_filename);
%
% if fid==-1
%     disp('read_header_from_bin: Can not open file for read!');
%     return;
% end

magic1 = fread(fid, 1, 'double');
tmp1 = fread(fid, 1, 'uint64');

magic2 = fread(fid, 1, 'double');
tmp2 = fread(fid, 1, 'uint64');

magic3 = fread(fid, 1, 'double');
tmp3 = fread(fid, 1, 'uint64');

magic4 = fread(fid, 1, 'double');
tmp4 = fread(fid, 1, 'uint64');

magic5 = fread(fid, 1, 'double');
fread(fid, 1, 'uint64');

magic6 = fread(fid, 1, 'double');
fread(fid, 1, 'uint64');

magic7 = fread(fid, 1, 'double');
fread(fid, 1, 'uint64');

magic8 = fread(fid, 1, 'double');
fread(fid, 1, 'uint64');

% fclose(fid);

fc_requested_magic = 73492.215;
fc_programmed_magic = -0.7923597;
fs_requested_magic = -189978508;
fs_programmed_magic = 93.126712;

reserve1_magic = -53243.129;
reserve2_magic = 0.0008123898;
reserve3_magic = -6.0098321;
reserve4_magic = 237.09983;

if magic1 == fc_requested_magic && ...
   magic2 == fc_programmed_magic && ...
   magic3 == fs_requested_magic && ...
   magic4 == fs_programmed_magic && ...
   magic5 == reserve1_magic && ...
   magic6 == reserve2_magic && ...
   magic7 == reserve3_magic && ...
   magic8 == reserve4_magic

    fc_requested = tmp1;
    fc_programmed = tmp2;
    fs_requested = tmp3;
    fs_programmed = tmp4;
end

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function [roiImg, imgXform, bb] = dtiRoiToImg(coords, imgXform, bb)
%
% [roiImg, imgXform, bb] = dtiRoiToImg(roi, [imgXform=eye(4)], [bb])
%
% If not the identity matrix, imgXform is typically set to xformToAcpc. This
% has the effect of forcing the roiImg to have the same voxel size as the
% image that xformToAcpc is based on.
%
% bb is the bounding box, defined in the xformed spaced. Defaults to
% [min(coords)-10; max(coords)+10].
%
% [roiImg, imgXform, bb] = dtiRoiToImg(roi);
% % Do some processing on the ROI
% % perimImg = bwperim(roiImg);
% perimRoi = dtiRoiFromImg(roiImg, imgXform, bb);
%
% HISTORY:
% 2009.08.19 RFD wrote it.

if(isstruct(coords))
    coords = coords.coords;
end

if(~exist('bb','var')||isempty(bb))
    bb = [min(coords)-10; max(coords)+10];
    if(~exist('imgXform','var')||isempty(imgXform))
        imgXform = eye(4);
        imgXform(1:3,4) = bb(1,:)'-1;
    end
end
if(~exist('imgXform','var')||isempty(imgXform))
    imgXform = eye(4);
end

sz = abs(diff(ceil(mrAnatXformCoords(inv(imgXform), bb))))+1;
roiImg = false((sz));

% Remove coords outside the bounding box
badCoords = coords(:,1)<bb(1,1) | coords(:,1)>bb(2,1) ...
          | coords(:,2)<bb(1,2) | coords(:,2)>bb(2,2) ...
          | coords(:,3)<bb(1,3) | coords(:,3)>bb(2,3);
coords = coords(~badCoords,:);

%coords(:,1) = coords(:,1) - bb(1,1) + 1;
%coords(:,2) = coords(:,2) - bb(1,2) + 1;
%coords(:,3) = coords(:,3) - bb(1,3) + 1;
coords = round(mrAnatXformCoords(inv(imgXform), coords));
roiImg(sub2ind(size(roiImg), coords(:,1), coords(:,2), coords(:,3))) = true;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function fileName = fgplvmWriteResult(model, dataSet, number)

% FGPLVMWRITERESULT Write a FGPLVM result.
% FORMAT
% DESC writess an FGPLVM result.
% ARG model : the model to write.
% ARG dataSet : the name of the data set to write.
% ARG number : the number of the FGPLVM run to write.
% RETURN fileName : the file name used to write.
%
% SEEALSO : fgplvmLoadResult
%
% COPYRIGHT : Neil D. Lawrence, 2009, 2011

% FGPLVM

dataSet(1) = upper(dataSet(1));
type = model.type;
type(1) = upper(type(1));
fileName = ['dem' dataSet type num2str(number)];

[kern, noise, fgplvmInfo, X] = fgplvmDeconstruct(model);

save(fileName, 'kern', 'noise', 'fgplvmInfo', 'X');