text:
function A=signedPolygonArea(vertices)
%%SIGNEDPOLYGONAREA Calculate the signed area of a non-self-intersecting
%                   2D polygon given its vertices. The magnitude of the
%                   area is the typical definition of area that one might
%                   expect. The sign of the area is positive if the points
%                   are arranged in a counterclockwise order around the
%                   polygon and negative if the points are in a clockwise
%                   order.
%
%INPUTS: vertices A 2XN list of the N vertices of the polygon in order
%                 around the polygon. Duplicate vertices are allowed. It
%                 does not matter whether or not the first vertex is
%                 repeated at the end.
%
%OUTPUTS: A The signed area of the polygon.
%
%The formula for computing the signed area of a non-self-intersecting
%polygon is taken from [1].
%
%REFERENCES:
%[1] Weisstein, Eric W. "Polygon Area." From MathWorld--A Wolfram Web
%    Resource. http://mathworld.wolfram.com/PolygonArea.html
%
%December 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numVertices=size(vertices,2);

A=0;
prevVertex=vertices(:,end);
for curIdx=1:numVertices
    curVertex=vertices(:,curIdx);
    A=A+det([prevVertex,curVertex]);
    prevVertex=curVertex;
end
A=A/2;

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
function [config,fn]=cosmo_config(fn, config)
% return a struc with configuration settings, or store such settings
%
% Usages:
% - get the configuation (either default, or what is in '.cosmomvpa.cfg'):
%   >> config=cosmo_config();
%
% - read configuration from a specified file
%   >> config=cosmo_config(fn)
%
% - store configuration in a file
%   >> cosmo_config(fn, to_store)
%
% - related function (in CoSMoMVPA's 'examples/' directory):
%   >> cosmo_wizard_set_config
%
% Inputs:
%   fn          optional filename of a configuration file.
%               This can either be a path of a file, or (on Unix platforms)
%               a filename in the user's home directory.
%               If fn is omitted, then it defaults to '.cosmomvpa.cfg'
%               and will read the configuration from that file if it
%               exists in one of the aforementioned locations; see the
%               notes below for details.
%   to_store    optional struct with configuration to store.
%               If omitted then the configuration is not stored.
%
% Returns:
%   config      Struct with configurations
%
% Notes:
%  - the rationale for this function is to keep the example code fixed
%    (that is, without any paths hard-coded) and still allow each
%    user to store the example data in a directory of their choice.
%
%  - the format for a configuration file is of the form <key>=<value>,
%    where <key> cannot contain the '=' character and <value> cannot start
%    or end with a white-space character.

%  - if the .cosmomvpa.cfg file is in the same directory as this file
%   (presently, 'mvpa/'), it will be found there as well.
%
%  - an example configuration file (of three lines):
%
% # path for runnable examples
% tutorial_data_path=/Users/nick/organized/_datasets/CoSMoMVPA/tutorial_data
% output_data_path=/Users/nick/organized/_datasets/CoSMoMVPA/tutorial_data
%
%  - If a configuration file with the name '.cosmomvpa.cfg' is stored
%    in a directory that is in the matlab path, the user path, or (on Unix
%    platforms) the user's home directory, then calling this function with
%    no arguments will read that file and return the configuration stored
%    in it.
%
%  - The configuration can also be changed as follows:
%    >> % get default configuration
%    >> % If the following command gives an error, do: config=struct()
%    >> config=cosmo_config();
%    >>
%    >> % update settings
%    >> config.tutorial_data_path='/path/to/some/data';
%    >> config.output_data_path='/where/i/want/output';
%    >>
%    >> % get first directory in user path
%    >> matlab_path=cosmo_strsplit (userpath,':',1);
%    >>
%    >> % set configuration filename
%    >> config_fn=fullfile(matlab_path,'.cosmomvpa.cfg');
%    >> cosmo_config(config_fn, config);
%    >> fprintf('Configuration stored in %s\n', config_fn);
%    >>
%    >> % now check they are the same
%    >> loaded_config=cosmo_config();
%    >> assert(isequal(loaded_config,config));
%
%  - The configuration can be set using cosmo_wizard_set_config in
%    CoSMoMVPA's 'examples/' directory
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    default_config_fn='.cosmomvpa.cfg';

    if nargin==1 && ischar(fn)
        fn=find_config_file(fn);
        config=read_config(fn);
    elseif nargin==0
        % see if the configuration file can be found
        fn=find_config_file(default_config_fn);
        if isempty(fn)
            config=struct();
        else
            config=read_config(fn);
        end
    elseif nargin==2 && ischar(fn) && isstruct(config)
        write_config(fn,config);
    else
        error('Illegal input');
    end

    expected_fields={'tutorial_data_path','output_data_path'};
    missing_fields=setdiff(expected_fields,fieldnames(config));

    if numel(missing_fields)>0
        missing_fields_str=cosmo_strjoin(missing_fields,''', ''');
        cosmo_warning(['Using %s, fields ''%s'' are missing. This makes '...
                    'it more complicated to run CoSMoMVPA''s '...
                    'exercises and examples.\n'...
                    'To set the configuration, consider running\n\n'...
                    '    cosmo_wizard_set_config\n\n'...
                    'in the CoSMoMVPA ''examples/'' directory'],...
                    mfilename(),missing_fields_str);
    end

    validate_config(config);

function validate_config(config)
    % simple validation of config

    % poor-man version of OO
    path_exists=struct();
    path_exists.match=@(x)isempty(cosmo_strsplit(x,'_path',-1));
    path_exists.test=@(p) exist(p,'file');
    path_exists.msg=@(key, p) sprintf('%s: path "%s" not found. ',key,p);

    checks={path_exists};
    add_msg=sprintf('To set the configuration, run: help %s', mfilename());

    % perform checks on fieldnames present in 'checks'.

    fns=fieldnames(config);
    for k=1:numel(fns)
        fn=fns{k};
        for j=1:numel(checks)
            check=checks{j};
            if check.match(fn)
                test_func=check.test;
                value=config.(fn);

                if ~test_func(value)
                    msg_func=check.msg;
                    cosmo_warning('%s\n%s',msg_func(fn, value),add_msg);
                end
            end
        end
    end

function path_fn=find_config_file(fn, raise_)
% tries to find a configuration file by looking:
% - for the path of the file
% - in the matlab path
% - in the user's home directory (on Unix)
    if nargin<2, raise_=false; end

    exist_=@(fn_) exist(fn_,'file');

    path_fn=[];

    % simulate 'go-to' statement using a while loop with break at the end
    while true
        % does the file exist 'as is'?
        if exist_(fn)
            path_fn=fn;
            break
        end

        % is it in the mvpa directory?
        parent_dir=fileparts(mfilename('fullpath'));
        m_fn=fullfile(parent_dir,fn);
        if exist_(m_fn)
            path_fn=m_fn;
            break;
        end

        % is it in the user path?
        % (not supported on octave)
        if cosmo_wtf('is_matlab')
            upaths=cosmo_strsplit(userpath(),pathsep());

            for k=1:numel(upaths)
                u_fn=fullfile(upaths{k},fn);
                if exist_(u_fn)
                    path_fn=u_fn;
                    break;
                end
            end
            if ~isempty(path_fn)
                break;
            end
        end

        if isunix()
            % is it in the home directory?
            u_fn=fullfile(getenv('HOME'),fn);
            if exist_(u_fn)
                path_fn=u_fn;
                break;
            end
        end

        if raise_
            error('Cannot find config file "%s"', fn);
        end

        break;
    end

function config=read_config(fn)
% reads configuration from a file fn

    config=struct(); % space for output
    if isempty(fn)
        error('Input file not specified - cannot open')
    end

    fid=fopen(fn);
    if fid==-1
        error('Unable to open file %s',fn);
    end

    line_number=0;
    while true
        % read each line
        line=fgetl(fid);
        if ~ischar(line)
            % end of file
            break;
        end

        line_number=line_number+1;

        % ignore empty lines or lines starting with '#'
        if isempty(line) || line(1)=='#'
            continue
        end

        % look for lines of form '<key>=<value>'.
        % white spaces around key or value are ignored.
        m=regexp(line,'(?<key>[^=]+)\s*=\s*(?<value>.*)\s*','names');

        if isempty(m)
            cosmo_warning('Skipping non-recognized line "%s"', line);
            continue;
        end

        % get key and value
        key=m.key;
        value=m.value;

        verify_no_illegal_char(fn,line_number,[key value])

        % see if it can be converted to numeric
        value_num=str2double(value);
        if ~isnan(value_num)
            value=value_num;
        end

        config.(key)=value;
    end

    fclose(fid);

function verify_no_illegal_char(fn,line_number,s)
    illegal_chars=sprintf('"''');
    matching=bsxfun(@eq,s,illegal_chars(:));
    if any(matching(:))
        error(['File %s, line %d: found illegal character. '...
            'note that quote characters are not allowed '...
            'for keys or values'],fn,line_number)
    end

function write_config(fn, config)
% writes the config to a file fn
% no support for comments or empty lines

    fid=fopen(fn,'w');

    fns=fieldnames(config);
    for k=1:numel(fns)
        fn=fns{k};

        v=config.(fn);
        if isnumeric(v)
            % convert numeric to string
            v=sprintf('%d ',v);
        elseif ischar(v)
            % no converstion
        else
            cosmo_warning('Skipping unsupported data type for key "%s"',...
                                                fn);
        end
        fprintf(fid,'%s=%s\n',fn,v);
    end

    fclose(fid);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function z = bsxfun( func, x, y, varargin )

global cvx___
oflag = cvx___.broadcast;
try
    cvx___.broadcast = true;
    z = func( x, y, varargin{:} );
    cvx___.broadcast = oflag;
catch exc
    cvx___.broadcast = oflag;
    cvx_throw( exc );
end

author: yu-jiang

repo: radpbook

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


text:
%% Workspace initialization
%  The workspace is cleaned and some parameters and deffined
clear all ; % this removes all variables stored in your current workspace
close all force ; % this closes all plotted figures
clc ; % this clears the command window log

% PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Set up the next parameters as your convenience.
%
% (2) <PCA_level> This parameter must be an integer between 0 and 100.
%     It will be use to compute how many components you need explain a
%     given level of variance in the data.
%     For example:
%     PCA_level = 95 ;
      PCA_level = 95 ;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp( ['================'] ) ;
disp( ['= PCA analysis ='] ) ;
disp( ['================'] ) ; disp( char(10) ) ;

%% Loading and examinating the data
%  Data are read from an Excel or CSV file (genes at columns and samples at rows)
%  which must have the next structure:
%  - First column: the name of the class (cell-type, treatment, stage, etc)
%  - Second column: the name of the sample
%  - First row: a header with the name of the genes at each column. Notice
%  that the first two cells of the header corresponds to the sample's
%  class and samples's name.
tic ;  disp( [' - Loading data...'] ) ; disp( char(10) ) ;
[ input_file,input_path ] = uigetfile( {'*.csv','Comma Sparated Values (*.csv)';'*.xls','Excel 97-2003 (*.xls)';'*.xlsx','Excel 2010 (*.xlsx)'},'MultiSelect','off' ) ;
[ vector ] = strsplit( input_file,'.' ) ;
file_name = vector(1) ;
extension = vector(2) ;
if ( strcmp(extension,'csv') )
  T = readtable( [input_path,input_file],'Delimiter',',','ReadVariableNames',false ) ;
  if ( sum(strcmp(T{:,end},'')) > 1 )
    T(:,end) = [] ;
  end%if
  classes = T{2:end,1} ;
  names = T{2:end,2} ;
  if ( sum(sum(isnan(char(T{:,end})))) > 0 )
    genes = T{1,3:end-1} ;
    data = str2double(T{2:end,3:end-1}) ;
  else
    genes = T{1,3:end} ;
    data = str2double(T{2:end,3:end}) ;
  end%if
  clear T ;
else
  if ( strcmp(extension,'xls') ||  strcmp(extension,'xlsx') )
    [ data,txt ] = xlsread( [input_path,input_file] ) ;
    classes = txt(2:end,1) ;
    names = txt(2:end,2) ;
    genes = txt(1,3:end) ;
    clear txt ;
  else
    disp('   file-format error.')
    return ;
  end%if
end%if
whos('data','names','classes','genes') ;
%  "data" is the matrix containing the normalized expression profiles from
%  single-cell qPCR data, with genes as columns and samples as rows. The
%  name of the genes at each column are in "genes", a row-vector.
%  Similarly, the name of each sample at each row of the "data" matrix are
%  stored in "cell_names", a column-vector. The names on these vectors must
%  be unique. The name of the classes at each row are stored in
%  "classes_names". This column vector has non-unique classes names.
N = size(data,1) ;
G = size(data,2) ;
classes_unique = unique( sort(classes) ) ;
C = size(classes_unique,1) ;
classes_keys = zeros(N,1) ;
K = size(classes_unique(:),1) ;
for k = 1:C
    classes_keys(strcmp(classes_unique(k),classes)) = k ;
  end%for

% plotting parameters
  colors = hsv( C ) ;
  colors_soft = hsv2rgb( rgb2hsv(colors).*[ones(C,1) 0.5*ones(C,1) ones(C,1)] ) ;
  markers = {'+' 'o' '*' 's' 'x' 'd' '.' 'v' 'p' 'h'}' ;
  for k=1:ceil( C/size(markers,1)-1 )
    markers = [ markers(:) ; markers(:) ] ;
  end%for

%% PCA
%  PCA is calculated using the corresponding function from the MATLAB's
%  statistical toolbox.
   disp( [char(10),'Performing PCA analysis...'] ) ;
   [coefficient,score,latent,tsquare,explained] = pca( data,'Algorithm','svd' ) ;
 % score*coefficient' = data - repmat(mean(data),size(data,1),1)
   PCA2 = rotatefactors(score(:,1:2),'Method','varimax') ;
   PCA3 = rotatefactors(score(:,1:3),'Method','varimax') ;
   IC = cumsum(explained) ;
   disp( ['  Information content on the two first component: ',num2str(round(IC(2))),'%.'] ) ;
   disp( ['  Information content on the three first component: ',num2str(round(IC(3))),'%.'] ) ;
% writing output
  output_file = strjoin( [ input_path,file_name,'__PCA.csv' ],'' ) ;
  fid = fopen(output_file,'w') ;
  for k = 1:N
    fprintf( fid,'%s,%s,%f,%f,%f\r\n',names{k},classes{k},PCA3(k,:) ) ;
  end%for
  fclose(fid) ;
% plotting the 2D PCA
    d = 0.1*( max(PCA2(:,1)) - min(PCA2(:,1)) ) ;
    x_min = min(PCA2(:,1))-d ;
    x_max = max(PCA2(:,1))+d ;
    d = 0.1*( max(PCA2(:,2)) - min(PCA2(:,2)) ) ;
    y_min = min(PCA2(:,2))-d ;
    y_max = max(PCA2(:,2))+d ;
    [ x y ] = meshgrid( x_min:(x_max-x_min)/100:x_max,y_min:(y_max-y_min)/100:y_max ) ;
    x = x(:) ; y = y(:) ;
    classification_result = classify( [x y],PCA2(:,1:2),classes,'quadratic' ) ;
    figure(1) ; set(1,'WindowStyle','docked') ;
  % set(1,'WindowStyle','normal','units','normalized','outerposition',[0 0 1 1]) ;
    hold on ; cla ; title(['2D PCA plot and discriminant analysis classification.']) ;
    if C > 1
    for k = 1:C
      index = strcmp( classification_result,classes_unique(k) ) ;
      scatter( x(index),y(index),100,colors_soft(k,:),'s','fill' ) ;
    end%for
    end%if
    h=[] ; text = [] ;
    for k = 1:C
      index = [classes_keys == k] ;
      text = [ text ; unique(classes(index)) ] ;
      h(k) = plot( PCA2(index,1),PCA2(index,2),markers{k},'Color',colors(k,:) ) ;
    end%for
    hold off ; grid on ;
    legend( h(:),text,'EdgeColor',[1 1 1],'Location','EO' ) ;
    xlabel('PC1 (first component)') ; xlim([x_min x_max]) ;
    ylabel('PC2 (second component)') ; ylim([y_min y_max]) ;
% plotting the 3D PCA
    d = 0.1*( max(PCA3(:,3)) - min(PCA3(:,3)) ) ;
    z_min = min(PCA3(:,3))-d ;
    z_max = max(PCA3(:,3))+d ;
    figure(2) ; title(['3D PCA analysis.']) ;
    set(2,'WindowStyle','docked') ;
    cla ; hold on ;
    for k = 1:C
      index = [classes_keys == k] ;
      plot3( PCA3(index,1),PCA3(index,2),PCA3(index,3),markers{k},'Color',colors(k,:),'MarkerFaceColor','none','LineWidth',0.1 ) ;
    end%for
    hold off ; grid on ; view(135,25) ;
    legend( text,'EdgeColor',[1 1 1],'Location','EO' ) ;
    xlabel('PC1 (first component)') ; xlim([x_min x_max]) ;
    ylabel('PC2 (second component)') ; ylim([y_min y_max]) ;
    zlabel('PC3 (third component)') ; zlim( [z_min z_max] ) ;
% plotting PCA information
  figure(3) ; set(3,'WindowStyle','docked') ; cla ;
  noc = min( sum(IC <= PCA_level)+1,G )  ;
  title(['PCA analysis. The ',num2str(PCA_level),'% is explained by the first ',num2str(noc),' copmonents.']) ;
  subplot(2,2,[1,2]) ; title(['Percentage of explained variance accumulated at each PCA component.']) ;
  hold on ; cla ;
   bar( IC(1:noc),1,'r' ) ;
   bar( IC.*([zeros(noc,1);ones(numel(IC)-noc,1)]),1,'b' ) ;
  hold off ;
  grid on ;
  legend( {['Less or equal than ',num2str(PCA_level),'%.'],['More than the ',num2str(PCA_level),'%.']},'EdgeColor',[1 1 1],'Location','NW' ) ;
  xlabel('component') ; xlim([0.5,size(IC,1)+0.5]) ;
  ylabel('% of variance') ; ylim( [0 100] ) ;
  subplot(2,2,3) ; title('PCA reconstruction, using all the components.') ;
  hold on ; cla ;
   imagesc( score*coefficient' + repmat(mean(data),size(data,1),1),[-1 1]*max(abs(data(:))) ) ;
   colormap( redgreencmap(100) ) ;
  hold off ;
  axis off ;
  subplot(2,2,4) ; title(['PCA approximation using ',num2str(noc),' of components.']) ;
  hold on ; cla ;
   imagesc( score*(coefficient.*[ones(G,noc) zeros(G,numel(IC)-noc)])' + repmat(mean(data),size(data,1),1),[-1 1]*max(abs(data(:))) ) ;
   colormap( redgreencmap(100) ) ;
  hold off ;
  axis off ;

%% Picking data-points from the 2D plot
  choice = questdlg(['Do you want to identify a super class in you dataset?'],'Select superclass','Yes','No','Yes');
  switch choice
  case 'Yes'
    choice='No, I want to do it again' ;
    msgbox({'Select the data you are interested in by:','',' - doing click on a single point ',' - clicking and dragging to select all points in an area','','Press RETURN when you were ready.',''}) ;
    while strcmp(choice,'No, I want to do it again')
    figure(5) ; clf ; set(5,'WindowStyle','docked') ;
    hold on ;
    title(['2D PCA']) ;
    gscatter( PCA2(:,1),PCA2(:,2),classes,colors,strjoin([markers(1:C)'],'') ) ;
    hold off ; grid on ;
    xlim( [x_min x_max] ) ; xlabel('PCA1') ;
    ylim( [y_min y_max] ) ; ylabel('PCA2') ;
    legend( 'Location','EO' ) ;
    h=gname(names,5) ;
    superclass=get(h,'String') ; superclass=unique([superclass{:}]) ;
    fileID = fopen( strjoin( [ input_path,file_name,'__PCA_superclass.txt' ],'' ),'w' ) ;
    fprintf(fileID,'%s\r\n',superclass) ;
    fclose(fileID) ;
    choice = questdlg(['You have selected ',num2str(size(h,1)),' data points.',char(10),'Do you want to keep this selection or do you want to try it again'],'Are you Happy with your selection?','Yes, it is fine','No, I want to do it again','No, I want to do it again');
    end%while
  end%switch

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function c = tmat_mxm ( a, b )

%*****************************************************************************80
%
%% TMAT_MXM multiplies two geometric transformation matrices.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Foley, van Dam, Feiner, Hughes,
%    Computer Graphics, Principles and Practice,
%    Addison Wesley, Second Edition, 1990.
%
%  Parameters:
%
%    Input, real A(4,4), the first geometric transformation matrix.
%
%    Input, real B(4,4), the second geometric transformation
%    matrix.
%
%    Output, real C(4,4), the product A * B.
%
  c(1:4,1:4) = a(1:4,1:4) * b(1:4,1:4);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Fig. 8.4   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%

clear all;
close all;

zgrid('new')
axis('square')
plot([-1 1],[0 0])
plot([0 0],[-1 1])
title('Fig. 8.4 Natural frequency and damping loci in z-plane')
xlabel('Re(z)')
ylabel('Im(z)')
hold off

author: Sable

repo: mcbench-benchmarks


text:
function fit_gls_brain(imgs, X, arorder, conditionnames, maskimg, varargin)
% Run Generalized Least Squares model with AR-model autoregression (optional)
% at each voxel in a set of images.
%
% :Usage:
% ::
%
%     fit_gls_brain(imgs, X, arorder, conditionnames, maskimg, ['contrasts', contrast_mtx, 'contrastnames', contrastnames], ['weights', weights])
%
% Single trial model to get trial amplitude, width, AUC, etc.
% Take one of those (e.g., AUC) and get a list of trial images
% Pass that into this function.
%
% fit_gls_brain(trial_amp_imgs, eventdesign{1}, contrasts
%
% You must add the intercept to X yourself if you want one!
%
% contrasts can be empty, and if it is, this function will write images
% for statistics on individual predictors.  If you enter contrast values,
% then it will write images for contrasts instead.
%
% conditionnames = column or contrast image names, in cell array
%
% e.g. eventnames{1} = {'high'    'medium'    'low'    'warm'}
%
% :Examples:
%
% This one looks @ significance for betas in X model:
% ::
%
%    load ../Multilev_mediation-try4(resliced)_10k/mediation_SETUP.mat
%    imgs = SETUP.data.M{1};
%    X = eventdesign{1};
%    arorder = 1;
%    contrasts = [];
%    conditionnames = eventnames{1};
%    maskimg = spm_select(1); % try gray matter mask...
%    fit_gls_brain(imgs, X, arorder, contrasts, conditionnames, maskimg)
%
% Now define contrasts and re-run on contrast values:
% ::
%
%    contrasts = [3 1 -1 -3; .25 .25 .25 .25; 1 0 0 -1]'
%    conditionnames = {'Linearpain' 'Average_resp' 'High-Low'};
%    fit_gls_brain(imgs, X, arorder, contrasts, conditionnames, maskimg)

    % ..
    %    Optional inputs
    % ..
    nobs = size(imgs, 1);
    weights = ones(nobs, 1);        % default weights
    contrasts = [];
    contrastnames = {};

    for i = 1:length(varargin)
        if ischar(varargin{i})
            switch varargin{i}
                % functional commands
                case 'contrasts', contrasts = varargin{i+1};

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

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

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

    % Check weights
    wtol = max(sqrt(weights))*eps^(1/3);
    if any(weights == 0) || any(weights < wtol)
        fprintf('Some weights are too low or badly scaled. Setting zero weights to %3.4f\n.', wtol);
        weights(weights < wtol) = wtol;
    end

    W = diag(weights);

    % [t, df, beta/contrast, Phi, sigma,stebeta, F] =
    % fit_gls(y,X,c,p,[PX, equal to pinv(X), for speed])

    % Map mask image into space of imgs, if not already, and write
    % 'mask.img'
    scn_map_image(maskimg, imgs(1, :), 'write', 'mask.img');
    maskimg = 'mask.img';

    % Setup inputs and function handle
    % -------------------------------------------------------------
    wh = intercept(X, 'which');
    if isempty(wh), disp('Warning! No intercept found in model.'); end

    %PX = pinv(X);
    PX = inv(X' * W * X) * X' * W;  % General weighted

    % Check if X is OK
    condx = cond(X);
    fprintf('Condition number of design matrix: %3.2f\n', condx);
    if condx > 10000, error('Design matrix is not estimable!'); end

    fhandle = @(y) fit_gls(y, X, contrasts, arorder, PX, weights);

    % Setup output image names
    % -------------------------------------------------------------
    % [beta, t, pvals, convals, con_t, con_pvals, sigma, Phi, df, stebeta, conste, F]

    %if isempty(contrasts), betaconname = 'beta'; else betaconname = 'contrast'; end

    names = cell(1, 12);
    names{7} = 'sigma.img';
    names{8} = 'phi.img';
    names{9} = 'df.img';
    names{12} = 'omnibus_F.img';

    for i = [1:6 10:11]  % Potential Multi-image outputs
        for j = 1:length(conditionnames)
            switch i
                case 1
                    names{i} = char(names{i}, ['beta_' conditionnames{j} '.img']);
                case 2
                    names{i} = char(names{i}, ['t_' conditionnames{j} '.img']);
                case 3
                    names{i} = char(names{i}, ['p_' conditionnames{j} '.img']);

                case 10
                    names{i} = char(names{i}, ['ste_beta_' conditionnames{j} '.img']);

            end
        end

        for j = 1:length(contrastnames)
            switch i
                case 4
                    names{i} = char(names{i}, ['con_beta_' contrastnames{j} '.img']);
                case 5
                    names{i} = char(names{i}, ['con_t_' contrastnames{j} '.img']);
                case 6
                    names{i} = char(names{i}, ['con_p_' contrastnames{j} '.img']);
                case 11
                    names{i} = char(names{i}, ['ste_con_' contrastnames{j} '.img']);
            end
        end

    end

    for i = [1:6 10:11]
        names{i} = names{i}(2:end, :);
    end

    save fit_gls_brain_SETUP imgs X arorder contrasts names conditionnames contrastnames maskimg

    % Run
    % -------------------------------------------------------------
    [beta, t, pvals, convals, con_t, con_pvals, sigma, Phi, df, stebeta, conste, F] = ...
        image_eval_function(imgs, fhandle, 'mask', maskimg, 'outimagelabels' , names);

    % other optional args: 'outimagelabels' , 'connames'

end

% % %
% % %
% % %             [t, df, beta, Phi, sigma, stebeta] = fit_gls(h, trial2ndX, [], 1, trial2ndpx);
% % %
% % %             conditiont(i(k), j(k), :) = t;
% % %             conditionste(i(k), j(k), :) = stebeta;
% % %             conditionmean(i(k), j(k), :) = beta;
% % %             trialphi(i(k), j(k), :) = Phi;
% % %
% % %             conditionp(i(k), j(k), :) = 2 .* (1 - tcdf(abs(t), df));
% % %
% % %             if ~isempty(c)
% % %                 stecons = diag(c' * diag(stebeta).^2 * c).^.5;
% % %                 con_est = c' * beta;
% % %                 tcons = con_est ./ stecons;
% % %
% % %                 contrastest(i(k), j(k), :) = con_est;
% % %                 contrastt(i(k), j(k), :) = tcons;
% % %                 contrastste(i(k), j(k), :) = stecons;
% % %                 contrastp(i(k), j(k), :) = 2 .* (1 - tcdf(abs(tcons), df));
% % %             end
% % %
% % %
% % %             % 2nd-level estimation of ratings model with AR(1) model
% % %             % use first param, ignore intercept
% % %             % -------------------------------------------------------------
% % %             if exist('ratings','var')
% % %                 [t, df, beta, Phi, sigma, stebeta] = fit_gls(h, ratings, [], 1, ratingspx);
% % %
% % %                 % 1st col. is rating effect
% % %                 ratingsest(i(k), j(k), 1) = beta(1);
% % %                 ratingst(i(k), j(k), 1) = t(1);
% % %                 ratingsste(i(k), j(k), 1) = stebeta(1);
% % %                 ratingsp(i(k), j(k), 1) = 2 .* (1 - tcdf(abs(t(1)), df));
% % %
% % %                 % 2nd column is intercept
% % %                 intcptest(i(k), j(k), 1) = beta(2);
% % %                 intcptt(i(k), j(k), 1) = t(2);
% % %                 intcptste(i(k), j(k), 1) = stebeta(2);
% % %                 intcptp(i(k), j(k), 1) = 2 .* (1 - tcdf(abs(t(2)), df));
% % %
% % %             end
% % %         end
% % %
% % %         if k == 1000, fprintf(1, '%3.0f s per 1000 vox.', etime(clock, et)), end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function b = uplus(a)
% function B=uplus(A)
%
% DESCRIPTION
%   Unary plus
%
% INPUTS
%   A: polynomial
%
% OUTPUTS
%   B: unary plus of input, +A
%
% SYNTAX
%   B = +A;
%     Unary plus.
%   B = uplus(A);
%     Function-call form for unary plus.

% 6/8/2002: PJS  Initial Coding

b = a;


text:
function xp = spm_beta_compare(alpha1,alpha2,Nsamp)
% Compute probability that r1 > r2
% FORMAT xp = spm_beta_compare(alpha1,alpha2,Nsamp)
%
% Input:
% alpha1    - Beta parameters for first density
% alpha2    - Beta parameters for second density
% Nsamp     - number of samples used to compute xp [default = 1e4]
%
% Output:
% xp        - exceedance probability
%
% Compute probability that r1 > r2 where p(r1)=Beta(r1|alpha1),
% p(r2)=Beta(r2|alpha2). Uses sampling.
% Useful for comparing groups in RFX model inference
%__________________________________________________________________________
% Copyright (C) 2009 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_beta_compare.m 3458 2009-10-13 10:05:35Z maria $

if nargin < 3
    Nsamp = 1e4;
end

Nk = length(alpha1);

% Perform sampling in blocks
%--------------------------------------------------------------------------
blk = ceil(Nsamp*Nk*8 / 2^28);
blk = floor(Nsamp/blk * ones(1,blk));
blk(end) = Nsamp - sum(blk(1:end-1));

xp = 0;
for i=1:length(blk)

    % Sample from univariate gamma densities then normalise
    % (see Dirichlet entry in Wikipedia or Ferguson (1973) Ann. Stat. 1,
    % 209-230)
    %----------------------------------------------------------------------
    r1 = zeros(blk(i),Nk);
    for k = 1:Nk
        r1(:,k) = spm_gamrnd(alpha1(k),1,blk(i),1);
    end
    sr1 = sum(r1,2);
    for k = 1:Nk
        r1(:,k) = r1(:,k)./sr1;
    end

    r2 = zeros(blk(i),Nk);
    for k = 1:Nk
        r2(:,k) = spm_gamrnd(alpha2(k),1,blk(i),1);
    end
    sr2 = sum(r2,2);
    for k = 1:Nk
        r2(:,k) = r2(:,k)./sr2;
    end

    % Exceedance probabilities:

    xp = xp + length(find (r1(:,1)>r2(:,1)));

end
xp = xp / Nsamp;

author: spm


text:
function editParamsROE(hObj,hData,hFig,planC)
% Edit model parameters
% AI 05/12/21

ud = guidata(hFig);
tag = get(hObj,'Tag');
indexS = planC{end};

%Get input data
idx = hData.Indices(1);
val = hData.EditData;
val2num = str2num(val);
if isempty(val2num) %Convert from string if numerical
    val2num = val;
end
prtcNum = ud.PrtcNum;
modelsC = ud.Protocols(prtcNum).model;
if isfield(ud,'ModelNum')
    modelNum = ud.ModelNum;
end

%Update parameter
switch(tag)
    case 'strSel'
        if hData.Indices(2)==1
            parameterS = modelsC{modelNum}.parameters;
            inputStructC = fieldnames(parameterS.structures);
            inputStructC = strcat('Select structure',{' '},inputStructC);
            matchIdx = find(strcmp(inputStructC,val));
            modelsC{modelNum}.inputStrNum = matchIdx;
            if ~isfield(modelsC{modelNum},'strNum') || ...
                    modelsC{modelNum}.strNum(matchIdx)==0
                hObj.Data{2} = 'Select from list';
            else
                allStrC = {planC{indexS.structures}.structureName};
                strIdx = modelsC{modelNum}.strNum(matchIdx);
                hObj.Data{2} = allStrC{strIdx};
            end
        else
            strListC = {'Select structure',planC{indexS.structures}.structureName};
            matchIdx = find(strcmp(strListC,val));
            inputStrNum = modelsC{modelNum}.inputStrNum;
            modelsC{modelNum}.strNum(inputStrNum) = matchIdx - 1;
            if isfield(ud.Protocols(prtcNum),'constraints')
                criteriaS = ud.Protocols(prtcNum).constraints;
                expectedStrName = strrep(hObj.Data{1},'Select structure ','');
                selectedStrName = strListC{matchIdx};
                %Update expected str name in criteria data stucture
                if isfield(criteriaS.structures,expectedStrName)
                    expS = criteriaS.structures.(expectedStrName);
                    selectedStrName = strrep(selectedStrName,' ','_');
                    criteriaS.structures.(selectedStrName) = expS;
                    criteriaS.structures = rmfield(criteriaS.structures,expectedStrName);
                    ud.Protocols(prtcNum).constraints = criteriaS;
                end
            end
        end
    case 'doseSel'
        if hData.Indices(2)==1
            return
        else
            dosListC = {'Select Plan',planC{indexS.dose}.fractionGroupID};
            matchIdx = find(strcmp(dosListC,val));
            %modelsC{modelNum}.planNum = matchIdx - 1;
            ud.planNum = matchIdx - 1;
            %Clear pre-computed dvh
            for p = 1:length(ud.Protocols)
                modelC = ud.Protocols(p).model;
                for modNum = 1:length(modelC)
                    if isfield(modelC{modNum},'dv')
                        modelC{modNum} = rmfield(modelC{modNum},'dv');
                    end
                end
                ud.Protocols(p).model = modelC;
            end
            %Auto-populate precribed dose if available
            RxField = ud.handle.tab1H(12);
            if ~isfield(planC{indexS.dose}(matchIdx - 1),'prescribedDose')
                try
                    prescribedDose = getPrescribedDose(matchIdx-1,planC);
                    set(RxField,'String',num2str(prescribedDose));
                    ud.handle.tab1H(12) = RxField;
                end
            else
                prescribedDose = planC{indexS.dose}(matchIdx - 1).prescribedDose;
                set(RxField,'String',num2str(prescribedDose));
                ud.handle.tab1H(12) = RxField;
            end
        end
    case 'fieldEdit'
        modelsC{modelNum} = modelsC{modelNum};
        parName = hObj.Data{idx,1};
        modelsC{modelNum}.(parName) = val2num;
        modelsC{modelNum} = modelsC{modelNum};
        set(ud.handle.tab1H(8),'Enable','On');  %Enable save
    case 'paramEdit'
        %Update modelC
        parName = hObj.Data{idx,1};
        strParam = 0;
        if isfield(modelsC{modelNum}.parameters,'structures')
            structS = modelsC{modelNum}.parameters.structures;
            if isstruct(structS)
                stC = fieldnames(structS);
                found = 0;
                t = 0;
                while ~found & t<numel(stC)
                    found = isfield(structS.(stC{t+1}),parName);
                    t = t+1;
                end
                if found
                    strParam = 1;
                    type = structS.(stC{t}).(parName).type;
                end
            end
        end
        if ~strParam
            type = modelsC{modelNum}.parameters.(parName).type;
        end
        if strcmpi(type{1},'bin')
            desc = modelsC{modelNum}.parameters.(parName).desc;
            ctgIdx = strcmp(desc,val2num);
            value = modelsC{modelNum}.parameters.(parName).cteg(ctgIdx);
            modelsC{modelNum}.parameters.(parName).val = value;
        else
            modelsC{modelNum}.parameters.(parName).val = val2num;
        end
        set(ud.handle.tab1H(8),'Enable','On');  %Enable save
end
ud.Protocols(prtcNum).model = modelsC;
guidata(hFig,ud);

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function I=montec2vvar(g,cx,dx,a,b,c,d,n)

%I=montec2vvar('g',cx,dx,a,b,c,d,n)
%Calculates integral on no rectangular space using monte-carlo method
%a,b, fixed boundaries
%cx,dx, variable boundaries (of x)
%c,d fixed boundaries (of y) that form a rectangular space with a and b
%boundaries.
t=rand(1,n);
u=rand(1,n);
x=a+t.*(b-a);			%extremos a,b fijos
y=c+u.*(d-c);			%extremos c,d variabales

cont=0;
s=0;

for i=1:n
   if (y(i)>feval(cx,x(i))&&y(i)<feval(dx,x(i)))
      cont=cont+1;
      fx=feval(g,x(i),y(i));
      s=s+fx;
   end
end

area=(simpson(dx,a,b,20)-simpson(cx,a,b,20))
I=((area./cont).*s);


text:
function failed_eeg_leadfield_units

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_convert_units ft_datatype_sens ft_convert_vol_sens ft_compute_leadfield current_dipole

[pnt, tri] = mesh_sphere(162);
sel = find(pnt(:,3)>0);

elec = [];
elec.pnt = pnt(sel,:) .* 10; % distributed on a 10 cm sphere
% elec.tra = eye(length(sel));
for i=1:length(sel)
  elec.label{i} = sprintf('electrode%d', i);
end
elec.unit = 'cm';

mesh = [];
mesh(1).pnt = pnt * 10.0; % in cm
mesh(1).tri = tri;
mesh(2).pnt = pnt *  9.5;
mesh(2).tri = tri;
mesh(3).pnt = pnt *  9.0;
mesh(3).tri = tri;

elec = ft_datatype_sens(elec);
elec = ft_convert_units(elec, 'm');
mesh = ft_convert_units(mesh, 'm');

vol0 = [];
vol0.type = 'infinite_currentdipole';
vol0.unit = 'm';

cfg = [];
cfg.method = 'singlesphere';
vol1 = ft_prepare_headmodel(cfg, mesh(1)); % only pass the outermost surface from the mesh

cfg = [];
cfg.method = 'concentricspheres';
cfg.conductivity = [1 1 1]; % functionally identical to a single sphere
vol2 = ft_prepare_headmodel(cfg, mesh);

cfg = [];
cfg.method = 'bemcp';
cfg.conductivity = [1 1 1]; % functionally identical to a single sphere
vol3 = ft_prepare_headmodel(cfg, mesh);

cfg = [];
cfg.method = 'openmeeg';
cfg.conductivity = [1 1 1]; % functionally identical to a single sphere
vol4 = ft_prepare_headmodel(cfg, mesh);

dip = [0 0 0.08]; % in meter

% this is to make a selection of the MEG channels
[vol0, elec] = ft_prepare_vol_sens(vol0, elec);
[vol1, elec] = ft_prepare_vol_sens(vol1, elec);
[vol2, elec] = ft_prepare_vol_sens(vol2, elec);
[vol3, elec] = ft_prepare_vol_sens(vol3, elec);
[vol4, elec] = ft_prepare_vol_sens(vol4, elec);

lf0 = ft_compute_leadfield(dip, elec, vol0);
lf1 = ft_compute_leadfield(dip, elec, vol1);
lf2 = ft_compute_leadfield(dip, elec, vol2);
lf3 = ft_compute_leadfield(dip, elec, vol3);
lf4 = ft_compute_leadfield(dip, elec, vol4);

n0 = norm(lf0);
n1 = norm(lf1);
n2 = norm(lf2);
n3 = norm(lf3);
n4 = norm(lf4);
% these should be relatively close to two, due to the mirror sources for the boundary condition
assert(abs(n1/n0-2)<0.3);
assert(abs(n2/n0-2)<0.3);
assert(abs(n3/n0-2)<0.3);
assert(abs(n3/n0-2)<0.3);
assert(abs(n4/n0-2)<0.3);

figure
ft_plot_dipole(dip, [1 0 0], 'unit', 'm');
ft_plot_headmodel(vol3, 'edgecolor', 'none', 'facealpha', 0.2);
ft_plot_sens(elec);
ft_plot_topo3d(elec.chanpos, lf3(:,1), 'facealpha', 0.6);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [cj,ier]=nufft2d2(nj,xj,yj,iflag,eps,ms,mt,fk)
%NUFFT2D2: Nonuniform FFT in R^2 - Type 2.
%
%  [CJ,IER] = NUFFT2D2(NJ,XJ,YJ,IFLAG,EPS,MS,MT,FK);
%
%     cj(j) = SUM   fk(k1,k2) exp(+/-i k1 xj(j)) exp(+/-i k2 yj(j))
%             k1,k2
%                            for j = 1,...,nj
%
%     where -ms/2 <= k1 <= (ms-1)/2, -mt/2 <= k2 <= (mt-1)/2
%
%     If (iflag .ge.0) the + sign is used in the exponential.
%     If (iflag .lt.0) the - sign is used in the exponential.
%
%  Input parameters:
%
%     nj     number of output values   (integer)
%     xj,yj  location of output values (real *8 array)
%     iflag  determines sign of FFT (see above)
%     eps    precision request  (between 1.0e-15 and 1.0e-1)
%     ms     number of Fourier modes given  [ -ms/2: (ms-1)/2 ]
%     mt     number of Fourier modes given  [ -mt/2: (mt-1)/2 ]
%     fk     Fourier coefficient values (complex *16 array)
%
%  Output parameters:
%
%     cj     output values (complex *16 array)
%     ier    error return code
%            ier = 0  => normal execution.
%            ier = 1  => precision eps requested is out of range.
%
%

cj=zeros(nj,1)+1i*zeros(nj,1);
ier=0;

mex_id_ = 'nufft2d2f90(i int[x], i double[], i double[], io dcomplex[], i int[x], i double[x], i int[x], i int[x], i dcomplex[], io int[x])';
[cj, ier] = nufft2d(mex_id_, nj, xj, yj, cj, iflag, eps, ms, mt, fk, ier, 1, 1, 1, 1, 1, 1);

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%SETFEATSIZE Fixed mapping, (re)set the feature size of a dataset
%
%   A = SETFEATSIZE(A,FEATSIZE)
%   A = A*SETFEATSIZE([],FEATSIZE)
%   A = A*SETFEATSIZE(FEATSIZE)
%
% INPUT
%   A         Dataset
%   FEATSIZE  Feature size vector, default: scalar with actual featsize
%
% OUTPUT
%   A         Dataset
%
% DESCRIPTION
% By default the feature size of a dataset is its number of features, i.e.
% the number of columns in the DATA field of A. If the features are samples
% of a multi-dimensional data item, e.g. the pixels of an image, the
% original size of this data item may be stored in FEATSIZE. The product of
% all elements in FEATSIZE has to be equal to the number of columns in the
% DATA field of A.
%
% This routine is particulary useful as A = A*SETFEATSIZE to remove image
% size settings for the feature size.
%
% There is a low level version of this routine with the same name in the
% PRDATASET directory.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com

function out = setfeatsize(varargin)

  varargin = shiftargin(varargin,'double');
  argin = setdefaults(varargin,[],[]);
  if mapping_task(argin,'definition')
    out = define_mapping(argin,'fixed')
  else
    % this will not happen. prmap directly jumps to @prdataset/setfeatsize
    [a,featsize] = deal(argin{:});
    if isempty(featsize)
      featsize = prod(getfeatsize(a));
    end
    out = setfeatsize(a,featsize);
  end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function chebyshev_polynomial_test17 ( )

%*****************************************************************************80
%
%% TEST17 tests T_TRIPLE_PRODUCT_INTEGRAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 April 2012
%
%  Author:
%
%    John Burkardt
%
  test_num = 20;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHEBYSHEV_POLYNOMIAL_TEST17:\n' );
  fprintf ( 1, '  T_TRIPLE_PRODUCT_INTEGRAL computes the triple integral\n' );
  fprintf ( 1, '  Tijk = integral ( -1 <= x <= 1 ) T(i,x) T(j,x) T(k,x) / sqrt ( 1-x^2) dx\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   I   J   K     Tijk           Tijk\n' );
  fprintf ( 1, '                 computed       exact\n' );
  fprintf ( 1, '\n' );

  n = 15;
  [ x, w ] = t_quadrature_rule ( n );

  seed = 123456789;

  for test = 1 : test_num
    [ i, seed ] = i4_uniform ( 2, 6, seed );
    [ j, seed ] = i4_uniform ( 1, 3, seed );
    [ k, seed ] = i4_uniform ( 0, 4, seed );
    fx1 = t_triple_product_integral ( i, j, k );
    fx2 = 0.0;
    for l = 1 : n
      ti = t_polynomial_value ( i, x(l) );
      tj = t_polynomial_value ( j, x(l) );
      tk = t_polynomial_value ( k, x(l) );
      fx2 = fx2 + w(l) * ti * tj * tk;
    end
    fprintf ( 1, '  %2d  %2d  %2d  %14.6g  %14.6g\n', i, j, k, fx1, fx2 );
  end

  return
end


text:
%Function to place dose into the running CERR plan.
%Just edit the strings below.
%JOD.

% [CTUniform3D, CTUniformInfoS] = getUniformizedCTScan;

register = 'UniformCT';  %Currently only option supported.  Dose has the same shape as the uniformized CT scan.
doseError = [];
fractionGroupID = 'CERR test';
doseEdition = 'CERR test';
description = 'Test PB distribution.'
overWrite = 'no';  %Overwrite the last CERR dose?
dose2CERR(dose3D,doseError,fractionGroupID,doseEdition,description,register,[],overWrite);

global planC;
global stateS;
indexS = planC{end};

stateS.doseToggle = 1;

stateS.doseSetChanged = 1;
stateS.CTDisplayChanged = 1;
stateS.structsChanged = 1;

stateS.doseSet = length(planC{indexS.dose});
%plancheckCallback('refresh');


text:
function [img, newInds] = rotateStack90(img, axis, varargin)
%ROTATESTACK90 Rotate a 3D image by 90 degrees around one image axis
%
%   RES = rotateStack90(IMG, AXIS);
%   IMG is a 3D image (either gray scale or color), and AXIS is the axis
%   number, in matrix convention: 1-> y axis, 2->x-axis, 3->z-axis.
%   AXIS can also be specified as letter: 'x', 'y' or 'z'.
%
%   RES = rotateStack90(IMG, AXIS, NUMBER);
%   Apply NUMBER rotation around the axis. NUMBER is the number of
%   rotations to apply, between 1 and 3. NUMER can also be negative, in
%   this case the rotation is performed in reverse direction.
%
%   [RES INDS] = rotateStack90(...);
%   Returns also informations about rotation.
%
%
%   Example
%   img = discreteBall(1:90, 1:100, 1:110, [50 50 50 30 20 10]);
%   img2 = rotateStack90(img, 'x');
%   figure(1); clf;
%   subplot(121); imshow(rot90(img(:,:,50)));
%   title('rotated slice');
%   subplot(122); imshow(squeeze(img2(50,:,:)));
%   title('slice of rotated image');
%
%   See also
%   imStacks, stackRotate90
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-05-18,    using Matlab 7.9.0.529 (R2009b)
% http://www.pfl-cepia.inra.fr/index.php?page=slicer
% Copyright 2010 INRA - Cepia Software Platform.

% parse axis, and check bounds
axis = parseAxisIndex_ijk(axis);

% positive or negative rotation
n = 1;
if ~isempty(varargin)
    n = varargin{1};
end

% ensure n is between 0 (no rotation) and 3 (rotation in inverse direction)
n = mod(mod(n, 4) + 4, 4);

% convert between 1 and 4
if n==0
    n = 4;
end

% all permutations: row for axis (ijk ordering), column for number of
% rotations, from 1 to 4
% permDimParams = {...
%     [2 3], [], [2 3]; ...
%     [1 3], [], [1 3]; ...
%     [1 2], [], [1 2] };
permDimParams = {...
    [1 3 2], [1 2 3], [1 3 2], [1 2 3]; ...
    [3 2 1], [1 2 3], [3 2 1], [1 2 3]; ...
    [2 1 3], [1 2 3], [2 1 3], [1 2 3] };

% all axis flipping: row for axis (ijk ordering), column for number of
% rotations, from 1 to 4
flipDimParams = {...
    3, [2 3], 2, []; ...
    1, [1 3], 3, []; ...
    2, [1 2], 1, [] } ;

newIndParams = {...
    [3 2 1], [1 2 3], [3 2 1], [1 2 3]; ...
    [1 3 2], [1 2 3], [1 3 2], [1 2 3]; ...
    [2 1 3], [1 2 3], [2 1 3], [1 2 3]};

% compute rotation params in ijk ordering
permInds = permDimParams{axis, n};
flipInds = flipDimParams{axis, n};
newInds  = newIndParams{axis, n};

% in case of a color image, need to adapt indices
if length(size(img))>3
    % convert index 3 to index 4
    permInds(permInds==3) = 4;
    % insert color index into indices array
    permInds = [permInds(1:2) 3 permInds(3)];

    % convert indices of dimensions to permute
    flipInds(flipInds==3) = 4;
end

% apply matrix dimension permutation
img = permute(img, permInds);

% depending on rotation, some dimensions must be fliped
for i=1:length(flipInds)
    img = flip(img, flipInds(i));
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function sF = rotate_outer(sF, rot)
% rotate a function by a rotation
%
% Syntax
%   sF = sF.rotate_outer(rot)
%
% Input
%  sF - @S2FunHarmonic
%  rot - @rotation
%
% Output
%  sF - @S2FunHarmonic
%

if sF.bandwidth ~= 0
  f = @(v) sF.eval(rotate(v, inv(rot)));
  sF = S2FunHarmonic.quadrature(f, 'bandwidth', sF.bandwidth);
end

end

author: mtex-toolbox

repo: mtex

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


text:
close all
clear all

% r1 = importdata('C:\Users\hakki\Downloads\camera-pipeline-UI-master\camera-pipeline-UI-master\image\r.txt');
% g1 = importdata('C:\Users\hakki\Downloads\camera-pipeline-UI-master\camera-pipeline-UI-master\image\g.txt');
% b1 = importdata('C:\Users\hakki\Downloads\camera-pipeline-UI-master\camera-pipeline-UI-master\image\b.txt');
% [m,n] = size(r1);
% image1 = zeros(m,n,3);
% image1(:,:,1) = r1;
% image1(:,:,2) = g1;
% image1(:,:,3) = b1;

% save image as binary file
% fd = fopen('x.txt','w');
%     fwrite(fd,image1(:,:,1)','double');
%     fclose(fd);
%
% fd2 = fopen('y.txt','w');
%     fwrite(fd2,image1(:,:,2)','double');
%     fclose(fd2);
% fd3 = fopen('z.txt','w');
%     fwrite(fd3,image1(:,:,3)','double');
%     fclose(fd3);

% image = im2double(image)/65536;
% figure, imshow(image1)

    fd1 = fopen('r.txt','r');
    fd2 = fopen('g.txt','r');
    fd3 = fopen('b.txt','r');
    r = fread(fd1,[ 3008 2000], 'double');
    fclose(fd1);
    g = fread(fd2,[ 3008 2000], 'double');
    fclose(fd2);
    b = fread(fd3,[ 3008 2000], 'double');
    fclose(fd3);

image = zeros(2000, 3008, 3);
image(:,:,1) = r';
image(:,:,2) = g';
image(:,:,3) = b';
image = im2double(image);
figure, imshow(image)

%% proof that txt files and stage 4 size is not equal.
% prophoto2 = im2double(imread('..\current_result.tif'));
%center aligned
eucliudian_error1 = sqrt((image(:,:,1)-...
    image1(:,:,1)).^2 + ...
    (image(:,:,2)-...
    image1(:,:,2)).^2 + ...
    (image(:,:,3)-...
    image1(:,:,3)).^2);
sum(eucliudian_error1(:))
figure,
imagesc(eucliudian_error1, [0 0.01]);

% abs()

author: karaimer

repo: camera-pipeline-UI


text:
function [A,B] = getLinSys(x,y,u)
%GETLINSYS
%    [A,B] = GETLINSYS(X,Y,U)

%    This function was generated by the Symbolic Math Toolbox version 6.0.
%    19-Nov-2014 20:59:12

A = reshape([0.0,cos(x),cos(y),0.0],[2, 2]);
if nargout > 1
    B = [sin(u);cos(u)];
end

author: MatthewPeterKelly

repo: dscTutorials


text:
%% hexMeshCylinder
% Below is a demonstration of the features of the |hexMeshCylinder| function

%%

clear; close all; clc;

%%
% Plot settings
fontSize=15;
faceAlpha1=1;
edgeColor=0.25*ones(1,3);
edgeWidth=2;

%% Creating a solid hexahedral mesh sphere
% Creating a solid hexahedral mesh sphere

pointSpacing=0.5;
cylRadius=1;
cylLength=2;

[meshStruct]=hexMeshCylinder(cylRadius,cylLength,pointSpacing);

% Access model element and patch data
Fb=meshStruct.facesBoundary;
Cb=meshStruct.boundaryMarker;
V=meshStruct.nodes;
E=meshStruct.elements;

%%
% Visualize mesh

hFig=cFigure;
subplot(1,2,1); hold on;
title('Boundary surfaces','FontSize',fontSize);
gpatch(Fb,V,Cb,'k',0.5,2);
% patchNormPlot(Fb,V);
axisGeom(gca,fontSize);
colormap(gjet); icolorbar;
camlight headlight;

hs=subplot(1,2,2); hold on;
title('Cut view of solid mesh','FontSize',fontSize);
optionStruct.hFig=[hFig hs];
gpatch(Fb,V,'kw','none',0.25);
meshView(meshStruct,optionStruct);
axisGeom(gca,fontSize);
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function test_suite = test_type()
% tests for cosmo_type
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_type_basics()
    fn=cosmo_make_temp_filename();
    cleaner=onCleanup(@()delete(fn));

    data=write_data(fn);

    s=cosmo_type(fn);
    assertEqual(s,data);

function test_dim_type_fprintf()
    if cosmo_skip_test_if_no_external('!evalc')
        return;
    end

    fn=cosmo_make_temp_filename();
    cleaner=onCleanup(@()delete(fn));

    data=write_data(fn);

    expr=sprintf('cosmo_type(''%s'')',fn);

    s=evalc(expr);
    s_fixed=regexprep(s,'ans\s*=\s*','');
    s_fixed=regexprep(s_fixed,'\s*$','');

    assertEqual(s_fixed,data);

function data=write_data(fn)
    fid=fopen(fn,'w');
    data=sprintf('foo\nbar%s',cosmo_make_temp_filename());
    fprintf(fid,data);
    fclose(fid);


text:
%--- help for abstvar/autocov ---
%
%  Compute the autocovariances (and the auto-correlation) of endogenous variables given the parameter values
%
%  ::
%
%     [C,R] = autocov(self);
%     [C,R] = autocov(self, params);
%     [C,R] = autocov(self, params, max_periods);
%
%  Args:
%     self (var object): var object
%     params (cell of struct): struct containing var model related parameters (default: [])
%     max_periods (integer): maximum number of period to calculate auto-covariance (default: 5)
%
%  Returns:
%     :
%
%     - **C** [4-dimensional array]: auto-covariance where dimensions correspond to
%
%        - 1,2: covariance
%        - 3: time lags
%        - 4: Number of parameters
%
%     - **R** [4-dimensional array]: auto-correlation with same dimensions
%
%

author: jmaih


text:
classdef winstonsource < spatiotemporal_database
   %UNTITLED2 Summary of this class goes here
   %   Detailed explanation goes here

   properties
   end

   methods
      function data = retrieve(where_, from_ , until_)
         error('unimplemented base class');
      end

   end

end

author: geoscience-community-codes

repo: GISMO


text:
% Initialize trial with default values
%
function trial = initTrial()
    % Recording system. Possible values are 'Axona' and 'NeuraLynx'.
    trial.system = nan;

    % Session base name. I.e. for single session '2103201101'. For combined sessions
    % it will be '2103201101+02'. More than two combined sessions are not supported.
    trial.basename = '';

    % Session path, normally the folder of the first session
    trial.path = '';

    trial.positions = []; % Matrix Nx5, N - number of entries.
                                % columns:
                                % 1 timestamps
                                % 2 LED1 x values
                                % 3 LED1 y values
                                % 4 LED2 x values, also known as x2
                                % 5 LED2 y values, also knows as y2

    trial.spikes = [];    % Matrix Nx3, N - number of entries. 1 column stores timestamps, spike data.
                                % 2 column contains tetrode number (integer), 3 column - unit
                                % number (integer).

    % Extra information, structure. Could contain fields
    % room - string with information about recording room;
    % shape - struct with information about recording box shape, fields:
    %         descr - raw string from file
    %         type - integer. 1 = box, 2 = cylinder, 3 = linear track
    %         value - double. Side length or diameter of the arena. If dimensions are not equal,
    %                 then value will have several values. For example for box arena value(1)
    %                 defines x-dimension, and value(2) defines y-dimension.
    % startTime - string
    % calibration - structure with following fields:
    %         file - Path to file with camera calibration information
    %         data - calibration data
    %
    % innerSize - structure with following fields:
    %         rec - Path to recording with LEDs that define arena inner size
    %         data - extracted inner size data
    % outerSize - structure with following fields
    %         rec - Path to recording with LEDs that define arena outer size
    %         data - extracted outer size data
    % subjectId - string, animal/subject ID
    %
    trial.extraInfo = struct();

    % Sampling time. Default values are from Axona system.
    trial.sampleTime = 0.02; % 50 Hz
    trial.videoSamplingRate = 50;

    trial.sessions = {}; % session names
    trial.cuts = {}; % cut filenames. Matrix of size NxM. Each row contains cut filenames
                     % for a particular tetrode (tetrodes are obtained from trial.units).
                     % Values in columns allow to load multiple cut files for cells,
                     % or have multiple cuts per tetrode in case of combined sessions.
    % Information about cut files that have been used to define units.
    % The key is an integer that represents the cell: each tetrode and cell number is
    % mapped in a single value by function helpers.tetrodeCellLinearIndex.
    % Each value is a Nx3 cell array. N is the number of cut files that have been used
    % to define this cell. First column contains path to cut file, second column
    % contains corresponding MD5 hash value. Both columns are strings.
    % Third column contains session index of a session to which this cut file belongs.
    % This is needed to adjust spike offset in case of multiple session loading.
    % Used to determine if cut files have been changed since last load of data in Matlab.
    trial.cutHash = containers.Map('KeyType', 'int64', 'ValueType', 'any');

    trial.tetrode = [];
    trial.units = []; % Matrix of units, Nx2. First column contains
                      % tetrode numbers, second - cell numbers.

    trial.startIndices = []; % When multiple sessions are combined, it can be
                             % useful to track position of each session individually.
                             % This vector stores starting index of each session. Index
                             % is applied to position matrix.

     % List of LFP channels to load/use. -1 has a special meaning being
     % 'all available channels'.
     trial.lfpChannels = [];

     % Flag that determines if this trial inherits units of a previous
     % trial or not. Default behaviour is that units are not inherited.
     trial.inheritUnits = false;

     % Identifier of cutting software. See also bntConstants.CutApp.
     trial.cutApp = [];

     % Fields that are going to be initialized later if at all:
     % lfpInfo - instance of io.LfpInfoInterface class
end

author: kavli-ntnu


text:
% Observation models and vision library.
%
% Sensor management
%   composeRobSen                   - Compose robot and sensor frames and parameters.
%
% Pin-hole model functions
%   distort                         - Distort projected point with radial distortion.
%   undistort                       - Undistorts projected point with radial distortion.
%   pixellise                       - Metric to pixellic conversion
%   depixellise                     - Pixellic to metric conversion
%   intrinsic                       - Build intrinsic matrix
%   invIntrinsic                    - Build inverse intrinsic matrix
%   essential                       - Essential matrix from frame specification.
%   fundamental                     - Fundamental matrix between 2 cameras.
%   pinHole                         - Pin-hole camera model, with optional radial distortion.
%   invPinHole                      - Inverse pin-hole camera model, with radial distortion correction.
%   pinHoleIdp                      - Pin-hole camera model for Inverse depth points, with radial distortion.
%   invPinHoleIdp                   - Inverse pin-hole camera model for IDP, with radial distortion correction.
%   pinHolePlucker                  - Projects plucker line.
%   invPinHolePlucker               - Retro-projects plucker line
%   aInvPinHolePlucker              - Inverse pin hole model for Plucker lines
%   pinHoleSegment                  - Pin hole projection of a segment.
%   pluckerInvCamera                - Inverse Plucker projection matrix
%   isVisible                       - Points visible from pinHole camera.
%   invDistortion                   - Radial distortion correction calibration.
%   invPinHoleHmg                   - Retro-project anchored homogeneous point AHP.
%   visibleSegment                  - Visible segment.
%   invPinHoleAPlucker              - Retro-projects anchored plucker line
%   invPinHoleIdpLin                - IDP line retro projection.
%   invPinHoleAhm                   - Retro-project anchored homogeneous point AHP.
%   persp_project                   - Project point into plane using pin-hole camera model
%   persp_retro                     - Retroproject pixel into 3D space.
%   pinHoleHmg                      - Pin-hole camera model for HMG points, with optional radial distortion.
%   projAhmPntIntoPinHole           - Project AHM pnt into pinhole camera.
%   projIdpPntIntoPinHole           - Project Idp pnt into pinhole.
%   invPinHoleDepth                 -
%   pinHoleDepth                    - Pin hole projection with distance measurement
%
% Functions related to Omnicam camera model
%   depixelliseOmniCam              - TODO: can inverse A up front - will save calculation esp in Jac part
%   invOmniCam                      - Inverse omnidirectional camera model
%   invOmniCamAhm                   - Retro-project anchored homogeneous point AHP.
%   invPinHoleAhmLin                - AHM line retro projection.
%   invPinHoleHmgLin                - IDP line retro projection.
%   omniCam                         - Gives projected pixel u of 3D point p for Omnidirectional Camera model
%   omni_project                    - Gives projected pixel u of 3D point p for Omnidirectional Camera model
%   omni_retro                      - Retroproject pixel into 3D space.
%   pixelliseOmniCam                - Omnicam Affine correction step
%
% Observation models, with 2 frame transforms
%   projEucPntIntoPinHoleOnRob      - Project Euc pnt into pinhole on robot.
%   projIdpPntIntoPinHoleOnRob      - Project Idp pnt into pinhole on robot.
%   projHmgPntIntoPinHoleOnRob      - Project Hmg pnt into pinhole on robot.
%   projSegLinIntoPinHoleOnRob      - Project segment line into pinhole on robot.
%   projPlkLinIntoPinHoleOnRob      - Project Plucker line into pinhole on robot.
%   projAplLinIntoPinHoleOnRob      - Project anchored Plucker line into pinhole on robot.
%   projIdpLinIntoPinHoleOnRob      - Project Idp line into pinhole on robot.
%   retroProjIdpPntFromPinHoleOnRob - Retro-project idp from pinhole on robot.
%   retroProjHmgPntFromPinHoleOnRob - Retro-proj. Hmg pnt from pinhole on rob.
%   retroProjPlkLinFromPinHoleOnRob - Retro-project Plucker line from pinhole on robot.
%   retroProjPlkEndPnts             - Retro project Plucker endpoints.
%   retroProjAplLinFromPinHoleOnRob - Retro-project anchored Plucker line from pinhole on robot.
%   retroProjAplEndPnts             - Retro project anchored Plucker endpoints.
%   retroProjIdpLinFromPinHoleOnRob - retroprj Idp Line from pinhole on robot.
%   projAhmPntIntoPinHoleOnRob      - Project Ahm pnt into pinhole on robot.
%   retroProjAhmPntFromPinHoleOnRob - Retro-project ahm from pinhole on robot.
%   projAhmLinIntoPinHoleOnRob      - Project Ahm line into pinhole on robot.
%   projHmgLinIntoPinHoleOnRob      - Project Hmg line into pinhole on robot.
%   retroProjAhmLinFromPinHoleOnRob - retroprj Ahm Line from pinhole on robot.
%   retroProjHmgLinFromPinHoleOnRob - retroprj Hmg Line from pinhole on robot.
%   projAhmPntIntoOmniCamOnRob      - Project Ahm pnt into omnidirectional camera on robot.
%   projEucPntIntoOmniCamOnRob      - Project Euc pnt into omnidirec cam on robot.
%   retroProjAhmPntFromOmniCamOnRob - Retro-project ahm from omnicam on robot.
%   projEucPntIntoPhdOnRob          - Project Eucliden point into Pinhole-depth in Robot
%   retroProjEucPntFromPhdOnRob     - Retro-proj Euc. pnt. from Pinhole-depth in Rob

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function event = read_brainvision_vmrk(filename)

% READ_BRAINVISION_VMRK reads the markers and latencies
% it returns the stimulus/response code and latency in ms.
%
% Use as
%   event = read_brainvision_vmrk(filename)
%
% See also READ_BRAINVISION_VHDR, READ_BRAINVISION_EEG

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

fid = fopen_or_error(filename, 'rt');

event = [];
line  = [];

readTime = ft_platform_supports('datetime');

while ischar(line) || isempty(line)
  line = fgetl(fid);
  if ~isempty(line) && ~(isnumeric(line) && line==-1)
    if startsWith(line, 'Mk')
      % this line starts with "Mk", so it probably contains a marker
      tok = strsplit(line, '=');
      if length(tok)~=2
        ft_warning('skipping unexpected formatted line in BrainVision marker file');
      else
        % the line looks like "MkXXX=YYY", which is ok
        % the interesting part now is in the YYY, i.e. the second token
        %tok = strsplit(tok{2}, ','); % this is not robust if there are
        %empty spaces between the ',', see github issue 2026
        tok = tokenize(tok{2}, ',');

        if isempty(tok{1})
          tok{1}  = [];
        end
        if isempty(tok{2})
          tok{2}  = [];
        end
        event(end+1).type     = tok{1};
        event(end  ).value    = tok{2};
        event(end  ).sample   = str2double(tok{3});
        event(end  ).duration = str2double(tok{4});
        if numel(tok)>5 && readTime
          try
            event(end).timestamp = datetime(tok{6}, 'InputFormat', 'yyyyMMddHHmmssSSSSSS');
          catch
            ft_warning('skipping invalid datetime in BrainVision marker file');
            readTime = false;
          end
        end
      end
    end
  end
end

fclose(fid);


text:
function centroid = circle_sector_centroid_2d ( r, center, theta1, theta2 )

%*****************************************************************************80
%
%% CIRCLE_SECTOR_CENTROID_2D returns the centroid of a circular sector in 2D.
%
%  Discussion:
%
%    A circular sector is formed by a circular arc, and the two straight line
%    segments that join its ends to the center of the circle.
%
%    A circular sector is defined by
%
%      ( X - CENTER(1) )^2 + ( Y - CENTER(2) )^2 = R^2
%
%    and
%
%      Theta1 <= Theta <= Theta2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Adrian Bowyer and John Woodwark,
%    A Programmer's Geometry,
%    Butterworths, 1983.
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%
%    Input, real CENTER(2,1), the center of the circle.
%
%    Input, real THETA1, THETA2, the angles defining the arc,
%    in radians.  Normally, THETA1 < THETA2.
%
%    Output, real CENTROID(2,1), the coordinates of the centroid
%    of the sector.
%
  theta = theta2 - theta1;

  if ( theta == 0.0 )
    d = 2.0 * r / 3.0;
  else
    d = 4.0 * r * sin ( 0.5 * theta ) / ( 3.0 * theta );
  end

  centroid(1,1) = center(1,1) + d * cos ( theta );
  centroid(2,1) = center(2,1) + d * sin ( theta );

  return
end

author: johannesgerer


text:
function pass = test_mean(pref)

if ( nargin == 0 )
    pref = chebfunpref();
end

%% Two arguments:
f = chebfun(@sin, pref);
g = chebfun(@cos, pref);
h = chebfun(@(x) .5*(sin(x) + cos(x)), pref);
pass(1) = normest(mean(f, g) - h) < 10*eps;

f = chebfun(@(x) [sin(x), cos(x)], pref);
g = chebfun(@(x) [cos(x), 1i*exp(x)], pref);
h = .5*(f+g);
pass(2) = normest(mean(f, g) - h) < 10*eps;

%% One argument:
f = chebfun(@sin, pref);
pass(3) = abs(mean(f)) < eps;
pass(4) = abs(mean(f.')) < eps;

f = chebfun(@(x) [sin(x), x], pref);
pass(5) = norm(mean(f), inf) < eps;
pass(6) = norm(mean(f.'), inf) < eps;

f = chebfun(@(x) [sin(x), x], [0, 6], pref);
pass(7) = norm(mean(f) - sum(f)/6, inf) < vscale(f)*eps;
pass(8) = norm(mean(f.') - sum(f).'/6, inf) < vscale(f)*eps;

%% singular function: a finite case

% define the domain:
dom = [-2 7];

op = @(x) sin(100*x)./((x-dom(1)).^0.5.*(x-dom(2)).^0.5);
f = chebfun(op, dom, 'exps', [-0.5 -0.5], 'splitting', 'on');
m = mean(f);
m_exact = -0.01273522016443600i;
pass(9) = abs(m - m_exact) < 1e5*eps*abs(m_exact);

%% singular function: an infinite case

% define the domain:
dom = [-2 7];

op = @(x) sin(100*x)./((x-dom(1)).^1.5.*(dom(2)-x).^0.5);
f = chebfun(op, dom, 'exps', [-1.5 -0.5], 'splitting', 'on');
m = mean(f);
pass(10) = ( isinf(m) );

%% singular function: a NaN case

% define the domain:
dom = [-2 7];

op = @(x) sin(98*x)./((x-dom(1)).^1.5.*(dom(2)-x).^1.5);
f = chebfun(op, dom, 'exps', [-1.5 -1.5], 'splitting', 'on');
m = mean(f);
pass(11) = ( isnan(m) );

%% Test for functions defined on unbounded domain:

% Functions on [-inf inf]:

% Set the domain:
dom = [-Inf 2 Inf];

op1 = @(x) x.^2.*exp(-x.^2);
op2 = @(x) (1-exp(-x.^2))./x.^2 + 2;
f = chebfun({op1 op2}, dom);
M = mean(f);
pass(12) = isnan(M);

% Function defined on [0 Inf]:

% Specify the domain:
dom = [0 Inf];

op = @(x) 0.75+sin(10*x)./exp(x);
f = chebfun(op, dom, 'splitting', 'on');
M = mean(f);
pass(13) = abs(M - 0.75) < 1e2*eps*vscale(f);

x = chebfun('x');
f = mean([x 3*x], 2);
pass(14) = norm(f-2*x) < 1e2*eps;

f = mean([x 3*x], 1);
pass(15) = norm( f - mean([x 3*x]) ) == 0;

pass(16) = (norm(mean(x,1) - mean(x.',2)) == 0);

pass(17) = (norm(mean(x') - 0) < 1e2*eps);

pass(18) = (norm(mean(x',1) - x') == 0);

end

function out = normest(f, dom)

% Generate a few random points to use as test values.
seedRNG(6178);
if ( nargin == 1 )
    x = 2 * rand(100, 1) - 1;
else
    x = sum(dom) * rand(10, 1) - dom(1);
end

out = norm(feval(f, x), inf);
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun