text:
function [ha,hf,lr] = unwrapVBvolatileOTO(posterior,out)
% OTO: recovers posterior sufficient statistics from VB inversion

if isempty(posterior)
    ha =[];
    return
end

[n,n_t] = size(posterior.muX);
if n==5

    ka = out.options.inF.lev2* VBA_sigmoid(posterior.muTheta(1), 'scale', out.options.inF.kaub);
    om = posterior.muTheta(2);
    mux = posterior.muX([2,4],:);
    sx = exp(posterior.muX([3,5],:));

    hf = figure('color',[1 1 1]);
    ha = subplot(2,2,1,'parent',hf,'ygrid','on','xlim',[1,n_t]);
    plotUncertainTimeSeries(mux,sx,1:n_t,ha);
    title(ha,'Learner''s belief estimate')
    legend(ha,{'x2 = outcome probability [LOG space]','x3 = expected volatility [LOG space]'})
    box(ha,'off')

    ha(2) = subplot(2,2,2,'parent',hf,'nextplot','add','ygrid','on','xlim',[1,n_t]);
    plot(ha(2),VBA_sigmoid(mux(1,:)));
    mv = exp(ka*mux(2,:)+om);
    plot(ha(2),mv,'g');

    % plot belief about outcome probability
    er1 = mux(1,:) - sqrt(sx(1,:));
    er2 = mux(1,:) + sqrt(sx(1,:));
    yp = [VBA_sigmoid(er2),fliplr(VBA_sigmoid(er1))];
    xp = [1:size(mux,2),fliplr(1:size(mux,2))];
    hfi = fill(xp,yp,'r','facecolor','b','edgealpha',0,'facealpha',0.25,'parent',ha(2));
    ev1 = exp(ka*(mux(2,:)-sqrt(sx(2,:)))+om);
    ev2 = exp(ka*(mux(2,:)+sqrt(sx(2,:)))+om);
    yp = [ev2,fliplr(ev1)];
    hfi = fill(xp,yp,'r','facecolor','g','edgealpha',0,'facealpha',0.25,'parent',ha(2));
    legend(ha(2),{'outcome probability','volatility'})
    box(ha(2),'off')

    x = posterior.muX;
    x(3,:) = exp(x(3,:));
    x(5,:) = exp(x(5,:));
    s1h = VBA_sigmoid(x(2,:)).*(1-VBA_sigmoid(x(2,:))); % likelihood precision
    s2h = x(3,:) + exp(ka*x(4,:)+om); % 2nd-level prediction variance
    lr = 1./(s2h.^-1 + s1h); % posterior variance
    ha(3) = subplot(2,2,3,'parent',hf,'nextplot','add','ygrid','on','xlim',[1,n_t]);
    plot(ha(3),lr);
    title(ha(3),'LEARNING RATE')
    box(ha(3),'off')

else

    ka = out.options.inF.lev2*VBA_sigmoid(posterior.muTheta(1),'scale',out.options.inF.kaub);
    om = posterior.muTheta(2);
    mux = posterior.muX([2,4,7,9],:);
    sx = exp(posterior.muX([3,5,8,10],:));

    hf = figure('color',[1 1 1]);
    ha = subplot(2,2,1,'parent',hf,'ygrid','on','xlim',[1,n_t]);
    plotUncertainTimeSeries(mux,sx,1:n_t,ha);
    title(ha,'Learner''s belief estimate (log-space)')
    legend(ha,{'1st cue: outcome probability','1st cue: expected volatility','2nd cue: outcome probability','2nd cue: expected volatility'})
    box(ha,'off')

    ha(2) = subplot(2,2,2,'parent',hf,'nextplot','add','ygrid','on','xlim',[1,n_t]);
    mr = mux([1,3],:);
    mv = exp(ka*mux([2,4],:)+om);
    col = getColors(4);
    for i=1:2
        plot(ha(2),VBA_sigmoid(mr(i,:)),'color',col(2*(i-1)+1,:));
        plot(ha(2),mv(i,:),'color',col(2*(i-1)+2,:));
    end
    er1 = mux([1,3],:) - sqrt(sx([1,3],:));
    er2 = mux([1,3],:) + sqrt(sx([1,3],:));
    ev1 = exp(ka*(mux(2,:)-sqrt(sx(2,:)))+om);
    ev2 = exp(ka*(mux(2,:)+sqrt(sx(2,:)))+om);
    for i=1:2
        yp = [VBA_sigmoid(er2(i,:)),fliplr(VBA_sigmoid(er1(i,:)))];
        xp = [1:size(mux,2),fliplr(1:size(mux,2))];
        hfi = fill(xp,yp,'r','facecolor',col(2*(i-1)+1,:),'edgealpha',0,'facealpha',0.25,'parent',ha(2));
        yp = [ev2,fliplr(ev1)];
        hfi = fill(xp,yp,'r','facecolor',col(2*(i-1)+2,:),'edgealpha',0,'facealpha',0.25,'parent',ha(2));
    end
    legend(ha(2),{'1st cue: outcome probability','1st cue: expected volatility','2nd cue: outcome probability','2nd cue: expected volatility'})
    box(ha(2),'off')
    title(ha(2),'Learner''s belief estimate')

    x = posterior.muX;
    x([3,5,8,10],:) = exp(x([3,5,8,10],:));
    s1h = VBA_sigmoid(x([2,7],:)).*(1-VBA_sigmoid(x([2,7],:))); % likelihood precision
    s2h = x([3,8],:) + exp(ka*x([4,9],:)+om); % 2nd-level prediction variance
    lr = 1./(s2h.^-1 + s1h); % posterior variance
    ha(3) = subplot(2,2,3,'parent',hf,'nextplot','add','ygrid','on','xlim',[1,n_t]);
    plot(ha(3),lr');
    title(ha(3),'LEARNING RATES')
    box(ha(3),'off')

end

X = [posterior.muX];
C = corrcoef(X');
ha(4) = subplot(2,2,4,'parent',hf);
imagesc(C,'parent',ha(4))
axis(ha(4),'square')
colorbar('peer',ha(4))
set(ha(4),'clim',[-1,1])
title(ha(4),'STATES'' CORRELATION')

VBA_getSubplots ();

author: MBB-team

repo: VBA-toolbox


text:
function onl_sendLSL(data,outlet)
% send a data structure over an LSL outlet

str = hlp_serialize(data);
outlet.push_sample({str});

repo: REST


text:
function signal = flt_clean_flatlines(varargin)
% Remove channels with abnormal data from a continuous data set.
% Signal = flt_clean_flatlines(Signal,MaxFlatlineDuration,MaxAllowedJitter)
%
% This is an automated artifact rejection function which ensures that
% the data contains no flat-lined channels.
%
% In:
%   Signal : continuous data set, assumed to be appropriately high-passed (e.g. >0.5Hz or
%            with a 0.5Hz - 2.0Hz transition band)
%
%   MaxFlatlineDuration : Maximum tolerated flatline duration. In seconds. If a channel has a longer
%                         flatline than this, it will be considered abnormal. Default: 5
%
%   MaxAllowedJitter : Maximum tolerated jitter during flatlines. As a multiple of epsilon.
%                      Default: 20
%
%   InitializeOn : Initialize on time range. If a time range is given as [start,end], either in
%                  seconds or as fractions of the whole data (if both <= 1), then where applicable
%                  the filter will be initialized only on that subset of the whole data. As a
%                  result, the filter will not have to be retrained in each cross-validation
%                  fold. (default: [])
% Out:
%   Signal : data set with flat channels removed
%
% Examples:
%   % use with defaults
%   eeg = flt_clean_flatlines(eeg);
%
% See also:
%   flt_clean_settings
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-07-06

% flt_clean_flatlines_version<1.0> -- for the cache

if ~exp_beginfun('filter') return; end;

declare_properties('name','FlatlineCleaning', 'follows','flt_selchans', 'precedes','flt_laplace', 'independent_channels',true, 'independent_trials','initialize_on');

arg_define(varargin, ...
    arg_norep({'signal','Signal'}), ...
    arg({'max_flatline_duration','MaxFlatlineDuration'}, 5, [0 Inf], 'Maximum tolerated flatline duration. In seconds. If a channel has a longer flatline than this, it will be considered abnormal.'), ...
    arg({'max_allowed_jitter','MaxAllowedJitter'}, 20, [0 Inf], 'Maximum tolerated jitter during flatlines. As a multiple of epsilon.'), ...
    arg({'initialize_on','InitializeOn'},[],[0 0 600 Inf],'Initialize on time range. If a time range is given as [start,end], either in seconds or as fractions of the whole data (if both <= 1), then where applicable the filter will be initialized only on that subset of the whole data. As a result, it will not have to be retrained in each cross-validation fold.','shape','row'),...
    arg_deprecated({'min_stddev','MinimumStdDev'}, 0.0001, [0 Inf], 'This parameter has no function any more.'), ...
    arg_norep('removed_channels',unassigned));

% flag channels
if ~exist('removed_channels','var')
    if ~isempty(initialize_on)
        ref_section = exp_eval(set_selinterval(signal,initialize_on,quickif(all(initialize_on<=1),'fraction','seconds')));
    else
        ref_section = signal;
    end
    removed_channels = [];
    for c=1:ref_section.nbchan
        zero_intervals = reshape(find(diff([false abs(diff(ref_section.data(c,:)))<(max_allowed_jitter*eps) false])),2,[])';
        if max(zero_intervals(:,2) - zero_intervals(:,1)) > max_flatline_duration*ref_section.srate
            removed_channels(end+1) = c; end
    end
end

% execute
if ~isempty(removed_channels)
    retain_channels = true(1,size(signal.data,1));
    retain_channels(removed_channels) = false;
    signal.data = signal.data(retain_channels,:,:,:,:,:,:,:);
    signal.chanlocs = signal.chanlocs(retain_channels);
    signal.nbchan = size(signal.data,1);
end

exp_endfun('append_online',{'removed_channels',removed_channels});


text:
%% QUAD_FSU uses the fsuClusterMatlab command to run the QUAD computation.
%
%  Discussion:
%
%    The arguments to fsuClusterMatlab mean:
%
%    [] allows us to specify an output directory;
%    [] allows us to specify switches to the MOAB queue handler, such as more time.
%    's' indicates this is a "simple" (task) job;
%    'w' indicates that the MATLAB session should wait for completion;
%    task_num is the number of tasks;
%    @quad_task indicates the function we want to execute;
%    args contains input arguments to each of the tasks.
%
%    The function "quad_task" must correspond to a MATLAB M-file "quad_task.m"
%    and that file must be in MATLAB's path.
%
%    The easiest way to ensure that the file is in MATLAB's path is to
%    create a subdirectory called "matlab" immediately below your home directory,
%    and place the function file there (as well as any other files needed by the job).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 April 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUAD_FSU:\n' );
  fprintf ( 1, '  Use MATLAB''s task computing on the FSU HPC cluster.\n' );
  fprintf ( 1, '  Here, we estimate an integral over an interval by\n' );
  fprintf ( 1, '  dividing the interval into subintervals.\n' );

  n = 100001;
  a = 0.0;
  b = 1.0;

  task_num = 4;
  n_task = 1 + floor ( ( n - 1 ) / task_num );

  args = {};
  for task = 1 : task_num
    a_task = ( task - 1 ) / task_num;
    b_task =   task       / task_num;
    args{task} = { n_task, a_task, b_task };
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Call fsuClusterMatlab:\n' );

  results = fsuClusterMatlab ( [], [], 's', 'w', task_num, @quad_task, args )
%
%  Convert RESULTS, which is a TASK_NUMx1 cell array, to
%  an ordinary MATLAB numeric array.
%
  results = cell2mat ( results );
%
%  Now we can use the SUM command to get our result.
%
  q = sum ( results );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Quadrature result is %f\n', q );

author: johannesgerer

repo: jburkardt-m


text:
function y = ExtendArray(x,s,fill)

%ExtendArray - Extend existing array to a given (larger) size.
%
%  USAGE
%
%    y = ExtendArray(x,s)
%
%    x              input array
%    s              output size, must have the same number of dimensions as x
%    fill           optional value for new elements (default: nan)
%

% Copyright (C) 2013 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% Check inputs
if nargin < 2,
	error('Incorrect number of parameters (type ''help ExtendArray'' for details).');
end
currentSize = size(x);
if length(currentSize) ~= length(s),
	error('Array and size do not have the same number of dimensions (type ''help ExtendArray'' for details).');
end

% Create extended array, fill with requested value
if nargin < 3,
	y = nan(s);
else
	if ~isscalar(fill),
		error('Incorrect fill value (type ''help ExtendArray'' for details).');
	end
	switch(fill),
		case 0,
			y = zeros(s);
		case 1,
			y = ones(s);
		case inf,
			y = inf(s);
		otherwise,
			y = fill*ones(s);
	end
end

% Copy old array into new array
indices = arrayfun(@(x) 1:x,currentSize,'uniformoutput',0);
y(indices{:}) = x;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function ExecFlag=DoUpdRateChk

global VSeq
global VCtl

ExecFlag = 1;
URTolerance = -abs(VCtl.MinUpdRate * 0.01);
% check update rate
% check rf
TxCoilID=unique(VSeq.rfCoilLine);
for i = 1: length(TxCoilID)
    Flags = VSeq.flagsLine(1,:);
    t = VSeq.tsLine(Flags~=0);
    t = t(VSeq.rfCoilLine == TxCoilID(i));
    dt =  diff(t);
    dt (dt < 0.5*eps) =[]; % avoid unique function rounding threshold (0.5*eps)
    if ~isempty(find((abs(dt)-VCtl.MinUpdRate)<URTolerance, 1))
        errordlg(['rf' num2str(TxCoilID(i)) ' sequence line exceeds minimum update rate!']);
        ExecFlag = 0;
    end
end

% check Gz
Flags = VSeq.flagsLine(2,:);
t = VSeq.tsLine(Flags~=0);
dt =  diff(t);
dt (dt < 0.5*eps) =[];
if ~isempty(find((abs(dt)-VCtl.MinUpdRate)<URTolerance, 1))
    errordlg('GzSS sequence line exceeds minimum update rate!');
    ExecFlag = 0;
end

% check Gy
Flags = VSeq.flagsLine(3,:);
t = VSeq.tsLine(Flags~=0);
dt =  diff(t);
dt (dt < 0.5*eps) =[];
if ~isempty(find((abs(dt)-VCtl.MinUpdRate)<URTolerance, 1))
    errordlg('GyP sequence line exceeds minimum update rate!');
    ExecFlag = 0;
end

% check Gx
Flags = VSeq.flagsLine(4,:);
t = VSeq.tsLine(Flags~=0);
dt =  diff(t);
dt (dt < 0.5*eps) =[];
if ~isempty(find((abs(dt)-VCtl.MinUpdRate)<URTolerance, 1))
    errordlg('GxR sequence line exceeds minimum update rate!');
    ExecFlag = 0;
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function esi = esiFeature(config, x)
%ESIFEATURE     ESI feature extraction.
%   Format: esi = esiFeature(config, x)
%   Inputs:
%       config: Configuration for ESI features.
%       x:      Input waveform or salience map.
%   Output:
%       esi:    ESI features, one column per frame.

    if (any(size(x) == 1)) x = salience(config, x); end
    esi = config.G * x;
    for f = 1:size(esi,2)
        esi(:,f) = esi(:,f) / sum(esi(:,f));
    end
end

author: MaigoAkisame


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         Demonstration of Basic Utilities of JSONlab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rngstate = rand ('state');
randseed=hex2dec('623F9A9E');
clear data2json json2data

if(exist('isequaln')==0)
    isequaln=@isequal;
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a simple scalar value \n')
fprintf(1,'%%=================================================\n\n')

data2json=pi
savejson('',data2json)
json2data=loadjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  an empty array \n')
fprintf(1,'%%=================================================\n\n')

data2json=[]
savejson('',data2json)
json2data=loadjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  an ampty string \n')
fprintf(1,'%%=================================================\n\n')

data2json=''
savejson('emptystr',data2json)
json2data=loadjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a simple row vector \n')
fprintf(1,'%%=================================================\n\n')

data2json=1:3
savejson('',data2json)
json2data=loadjson(ans)
if(~isequaln(json2data,data2json))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a simple column vector \n')
fprintf(1,'%%=================================================\n\n')

data2json=(1:3)'
savejson('',data2json)
json2data=loadjson(ans)
if(~isequaln(json2data,data2json))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a string array \n')
fprintf(1,'%%=================================================\n\n')

data2json=['AC';'EG']
savejson('',data2json)
json2data=loadjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a string with escape symbols \n')
fprintf(1,'%%=================================================\n\n')

data2json=sprintf('AB\tCD\none"two')
savejson('str',data2json)
json2data=loadjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a mix-typed cell \n')
fprintf(1,'%%=================================================\n\n')

data2json={'a',true,[2;3]}
savejson('',data2json)
json2data=loadjson(ans)
if(~isequaln(json2data,data2json))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a 3-D array in nested array form\n')
fprintf(1,'%%=================================================\n\n')

data2json=reshape(1:(2*4*6),[2,4,6]);
savejson('',data2json,'NestArray',1)
json2data=loadjson(ans)
if(~isequaln(json2data,data2json))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a 3-D array in annotated array form\n')
fprintf(1,'%%=================================================\n\n')

data2json=reshape(1:(2*4*6),[2,4,6]);
savejson('',data2json,'NestArray',0)
json2data=loadjson(ans)
if(~isequaln(json2data,data2json))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a 4-D array in annotated array form\n')
fprintf(1,'%%=================================================\n\n')

data2json=reshape(1:(2*4*3*2),[2,4,3,2]);
savejson('',data2json,'NestArray',0)  % nestarray for 4-D or above is not working
json2data=loadjson(ans)
if(~isequaln(json2data,data2json))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a 3-D array in nested array form (JSONLab 1.9)\n')
fprintf(1,'%%=================================================\n\n')

data2json=reshape(1:(2*4*6),[2,4,6]);
savejson('',data2json,'NestArray',1,'FormatVersion',1.8)
json2data=loadjson(ans,'FormatVersion',1.8)
if(~isequaln(json2data,data2json))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a 3-D array in annotated array form (JSONLab 1.9 or earlier)\n')
fprintf(1,'%%=================================================\n\n')

data2json=reshape(1:(2*4*6),[2,4,6]);
savejson('',data2json,'NestArray',0,'FormatVersion',1.8)
json2data=loadjson(ans,'FormatVersion',1.8)
if(~isequaln(json2data,data2json))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a complex number\n')
fprintf(1,'%%=================================================\n\n')

data2json=1+2i
savejson('',data2json)
json2data=loadjson(ans)
if(~isequaln(json2data,data2json))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a complex matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=magic(6);
data2json=data2json(:,1:3)+data2json(:,4:6)*1i
savejson('',data2json)
json2data=loadjson(ans)
if(~isequaln(json2data,data2json))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  MATLAB special constants\n')
fprintf(1,'%%=================================================\n\n')

data2json=[NaN Inf -Inf]
savejson('specials',data2json)
json2data=loadjson(ans)
if(~isequaln(json2data.specials,data2json))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a real sparse matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=sprand(10,10,0.1)
savejson('sparse',data2json,'FloatFormat','%.18g')
json2data=loadjson(ans)
if(~isequaln(json2data.sparse,data2json))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a complex sparse matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=sprand(10,10,0.1);
data2json=data2json-data2json*1i
savejson('complex_sparse',data2json,'FloatFormat','%.18g')
json2data=loadjson(ans)
if(~isequaln(json2data.complex_sparse,data2json))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  an all-zero sparse matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=sparse(2,3);
savejson('all_zero_sparse',data2json)
json2data=loadjson(ans)
if(~isequaln(json2data.all_zero_sparse,data2json))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  an empty sparse matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=sparse([]);
savejson('empty_sparse',data2json)
json2data=loadjson(ans)
if(~isequaln(json2data.empty_sparse,data2json))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  an empty 0-by-0 real matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=[];
savejson('empty_0by0_real',data2json)
json2data=loadjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  an empty 0-by-3 real matrix\n')
fprintf(1,'%%=================================================\n\n')

data2json=zeros(0,3);
savejson('empty_0by3_real',data2json)
json2data=loadjson(ans)
if(~isequaln(json2data.empty_0by3_real,data2json))
    warning('conversion does not preserve original data');
end

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a sparse real column vector\n')
fprintf(1,'%%=================================================\n\n')

data2json=sparse([0,3,0,1,4]');
savejson('sparse_column_vector',data2json)
json2data=loadjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a sparse complex column vector\n')
fprintf(1,'%%=================================================\n\n')

data2json=data2json-1i*data2json;
savejson('complex_sparse_column_vector',data2json)
json2data=loadjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a sparse real row vector\n')
fprintf(1,'%%=================================================\n\n')

data2json=sparse([0,3,0,1,4]);
savejson('sparse_row_vector',data2json)
json2data=loadjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a sparse complex row vector\n')
fprintf(1,'%%=================================================\n\n')

data2json=data2json-1i*data2json;
savejson('complex_sparse_row_vector',data2json)
json2data=loadjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a structure\n')
fprintf(1,'%%=================================================\n\n')

data2json=struct('name','Think Different','year',1997,'magic',magic(3),...
                 'misfits',[Inf,NaN],'embedded',struct('left',true,'right',false))
savejson('astruct',data2json,struct('ParseLogical',1))
json2data=loadjson(ans)
class(json2data.astruct.embedded.left)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a structure array\n')
fprintf(1,'%%=================================================\n\n')

data2json=struct('name','Nexus Prime','rank',9);
data2json(2)=struct('name','Sentinel Prime','rank',9);
data2json(3)=struct('name','Optimus Prime','rank',9);
savejson('Supreme Commander',data2json)
json2data=loadjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a cell array\n')
fprintf(1,'%%=================================================\n\n')

data2json=cell(3,1);
data2json{1}=struct('buzz',1.1,'rex',1.2,'bo',1.3,'hamm',2.0,'slink',2.1,'potato',2.2,...
              'woody',3.0,'sarge',3.1,'etch',4.0,'lenny',5.0,'squeeze',6.0,'wheezy',7.0);
data2json{2}=struct('Ubuntu',['Kubuntu';'Xubuntu';'Lubuntu']);
data2json{3}=[10.04,10.10,11.04,11.10]
savejson('debian',data2json,struct('FloatFormat','%.2f'))
json2data=loadjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  invalid field-name handling\n')
fprintf(1,'%%=================================================\n\n')

json2data=loadjson('{"ValidName":1, "_InvalidName":2, ":Field:":3, "项目":"绝密"}')

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a function handle\n')
fprintf(1,'%%=================================================\n\n')

data2json=@(x) x+1
savejson('handle',data2json)
json2data=loadjson(ans)

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a 2D cell array\n')
fprintf(1,'%%=================================================\n\n')

data2json={{1,{2,3}},{4,5},{6};{7},{8,9},{10}};
savejson('data2json',data2json)
json2data=loadjson(ans)  % only savejson works for cell arrays, loadjson has issues

fprintf(1,'\n%%=================================================\n')
fprintf(1,'%%  a 2D struct array\n')
fprintf(1,'%%=================================================\n\n')

data2json=repmat(struct('idx',0,'data','structs'),[2,3])
for i=1:6
    data2json(i).idx=i;
end
savejson('data2json',data2json)
json2data=loadjson(ans)

if(exist('datetime'))
    fprintf(1,'\n%%=================================================\n')
    fprintf(1,'%%  datetime object \n')
    fprintf(1,'%%=================================================\n\n')

    data2json=datetime({'8 April 2015','9 May 2015'}, 'InputFormat','d MMMM yyyy')
    savejson('',data2json)
    json2data=loadjson(ans)
end

if(exist('containers.Map'))
    fprintf(1,'\n%%=================================================\n')
    fprintf(1,'%%  a container.Maps object \n')
    fprintf(1,'%%=================================================\n\n')

    data2json=containers.Map({'Andy','William','Om'},{21,21,22})
    savejson('',data2json)
    json2data=loadjson(ans)
end

if(exist('istable'))
    fprintf(1,'\n%%=================================================\n')
    fprintf(1,'%%  a table object \n')
    fprintf(1,'%%=================================================\n\n')

    Names={'Andy','William','Om'}';
    Age=[21,21,22]';
    data2json=table(Names,Age)
    savejson('table',table(Names,Age))
    json2data=loadjson(ans)
end

if(exist('bandwidth'))
    fprintf(1,'\n%%=================================================\n')
    fprintf(1,'%%  use _ArrayShape_ \n')
    fprintf(1,'%%=================================================\n\n')

    lband=2;
    uband=3;
    data2json=spdiags(true(8,lband+uband+1),-uband:lband,5,8);
    data2json=full(double(data2json));
    data2json(data2json~=0)=find(data2json)

    savejson('',data2json,'usearrayshape',1)
    json2data=loadjson(ans,'fullarrayshape',1)

    savejson('',tril(data2json),'usearrayshape',1)
    json2data=loadjson(ans,'fullarrayshape',1)

    savejson('',triu(data2json+1i*data2json),'usearrayshape',1)
    json2data=loadjson(ans,'fullarrayshape',1)

    savejson('',tril(triu(int8(data2json))),'usearrayshape',1)
    json2data=loadjson(ans,'fullarrayshape',1)

    savejson('',data2json(:,1:5)+data2json(:,1:5)','usearrayshape',1)
    json2data=loadjson(ans,'fullarrayshape',1)
end

try
    val=zlibencode('test');
    fprintf(1,'\n%%=================================================\n')
    fprintf(1,'%%  a 2-D array in compressed array format\n')
    fprintf(1,'%%=================================================\n\n')

    data2json=eye(10);
    data2json(20,1)=1;
    savejson('',data2json,'Compression','zlib','CompressArraySize',0)  % nestarray for 4-D or above is not working
    json2data=loadjson(ans)
    if(~isequaln(json2data,data2json))
        warning('conversion does not preserve original data');
    end
catch
end

rand ('state',rngstate);

author: fangq

repo: jsonlab

save_path: github-repos/MATLAB/fangq-jsonlab


text:
function r = gamrand(a, b, varargin);
%GAMRAND Random matrices from gamma distribution.
%
%   R = GAMRAND(A,B) returns a matrix of random numbers chosen
%   from the gamma distribution with parameters A and B.
%   The size of R is the common size of A and B if both are matrices.
%   If either parameter is a scalar, the size of R is the size of the other
%   parameter. Alternatively, R = GAMRAND(A,B,M,N) returns an M by N matrix.
%
%   Note: Parameterization as in (Neal, 1996).
%      A is mean of the distribution
%      B is degrees of freedom
%
%	See also INVGAMRAND
%
% Copyright (c) 1999 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

error('No mex-file for this architecture. See Matlab help and convert.m in ./linuxCsource or ./winCsource for help.')

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
    if inpr3 == 1
        FILNAME='private/gshhs_c.b';
        m_gshhs_c('save','coastl.mat');

    elseif inpr3 == 2
        FILNAME='private/gshhs_l.b';
        m_gshhs_l('save','coastl.mat');

    elseif inpr3 == 3
        FILNAME='private/gshhs_i.b';
        m_gshhs_i('save','coastl.mat');

    elseif inpr3 == 4
        FILNAME='private/gshhs_h.b';
        m_gshhs_h('save','coastl.mat');
    end
    load coastl.mat

    coastline = ncst;
    update(mainmap())
    clear  ncst coastl

author: CelsoReyes


text:
classdef testShowingError < test

    properties (Access = protected)
        error
    end

    properties (Abstract, Access = protected)
        tol
    end

    methods (Access = protected)
        function printTestPassed(obj)
            cprintf('green',obj.FileName);
            cprintf('green',' PASSED.');
            cprintf('black',['Error: ',num2str(obj.error),'\n']);
        end

        function printTestNotPassed(obj)
            cprintf('red',obj.FileName);
            cprintf('red',' FAILED.');
            cprintf('red',['Error: ',num2str(obj.error),'\n']);
        end

        function hasPassed = hasPassed(obj)
            obj.computeError()
            hasPassed = obj.error < obj.tol();
        end
    end

    methods (Abstract, Access = protected)
        computeError(obj)
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function lambda = combin_eigenvalues ( alpha, beta, n )

%*****************************************************************************80
%
%% COMBIN_EIGENVALUES returns the eigenvalues of the COMBIN matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, BETA, scalars that define A.
%
%    Input, integer N, the order of A.
%
%    Output, real LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );

  lambda(1:n-1,1) = alpha;
  lambda(n,1) = alpha + n * beta;

  return
end


text:
% circpatch() - draw a dashed arc between two points
%
% Usage:
%   >> [handles] = circpatch( X, Y, curv, color, linewidth, segments, offset );
%
% Inputs:
%   X          - abscissae of the points (2 abscices for the 2 points)
%   Y          - ordinates of the points (2 ordinates for the 2 points)
%   curv       - curvature (0 = no curvature, 1 = round, -1 = round in
%                the other direction)
%   color      - line color (default black)
%   linewidth  - line thickness (same as line 'linewidth' property, default:1)
%   segments   - number of segments composing each continuous part of the arc (default:50)
%   offset     - change the offset of the dash (range: [0,1]) (default:0)
%   middle     - put a small line at the middle of the arc
%
% Outputs:
%   handles    - handles of all the objects composing the arc

% arno@salk.edu, Arnaud Delorme, CNL / Salk Institute, 2001

% This program is free software; you can redistribute it and/or modify it.
% 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.

function h = circpatch( X, Y, circfactor, color, thickness, segments, offset, middle );

PADSTRIPE  = 0.8; % percentage
SIZESTRIPE = 0.2; % percentage

if nargin < 3
	error('Not enough arguments. Type help circpatch');
	return;
end;
if nargin < 4
	color = 'k';
end;
if nargin < 5
	thickness = 1;
end;
if nargin < 6
	segments = 50;
end;
if nargin < 7
	offset = 0;
end;
if nargin < 8
	middle = 0;
end;

% middle of X and Y
% -----------------
middleXY = [(X(1)+X(2))/2 (Y(1)+Y(2))/2];

% normal to this line
% -------------------
diffXY = [ (X(2)-X(1))/2 (Y(2)-Y(1))/2];
normdiffXY = [diffXY(2) -diffXY(1)];

% third coordinate
% ----------------
fisrtPoint  = [ X(1) Y(1) ];
secondPoint = [ X(2) Y(2) ];
thirdPoint  = middleXY + normdiffXY*circfactor;

% compute the two lines equation
% ------------------------------
a1 = normdiffXY(2)/normdiffXY(1);
b1 = thirdPoint(2)-thirdPoint(1)*a1;

middle13 = [(fisrtPoint(1)+thirdPoint(1))/2 (fisrtPoint(2)+thirdPoint(2))/2];
a2 = [1 (fisrtPoint(2)-thirdPoint(2))/(fisrtPoint(1)-thirdPoint(1))];
a2 = [a2(2) -a2(1)]; % take the normal
a2 = a2(2)/a2(1);
b2 = middle13(2)-middle13(1)*a2;

% solve the equation (the center of the circle is the solution)
% -------------------------------------------------------------
centreX = (b1-b2)/(a2-a1);
centreY = a1*centreX + b1;
centreY = a2*centreX + b2;
radius  = sqrt((fisrtPoint(1)-centreX)^2+(fisrtPoint(2)-centreY)^2);

%circle( fisrtPoint(1), fisrtPoint(2), 0.05);
%circle( secondPoint(1), secondPoint(2), 0.05);
%circle( thirdPoint(1), thirdPoint(2), 0.05);

% compute the angle of the points
% -------------------------------
angle1 = atan( (fisrtPoint(2)-centreY)/(fisrtPoint(1)-centreX) )/pi*180;
angle2 = atan( (secondPoint(2)-centreY)/(secondPoint(1)-centreX) )/pi*180;
if ((fisrtPoint(1)-centreX) < 0) % negative cosinus
	angle1 = 180+angle1;
end;
if ((secondPoint(1)-centreX) < 0)
	angle2 = angle2+180;
end;
if angle1 < 0 angle1 = angle1+360; end;
if angle2 < 0 angle2 = angle2+360; end;

if circfactor<0 % which side
	tmp = angle1;
	angle1 = angle2;
	angle2 = tmp;
end;
if angle2 < angle1 angle2 = angle2+360; end;

% if offset exist, drawstripes
% ----------------------------
if offset ~= -1
	offsetshifted = mod(0.5 + offset + SIZESTRIPE/2,1);
	vectpad    = PADSTRIPE*abs(angle2-angle1);
	vectstripe = SIZESTRIPE*abs(angle2-angle1);
	nblines    = 1/(PADSTRIPE+SIZESTRIPE);
	currentangle = angle1;
	cont =1;
	h = [];
	firstpass = 1;
	while cont
		if firstpass == 1
			if offsetshifted/nblines < SIZESTRIPE
				currentangle = currentangle + vectstripe*(offsetshifted/nblines)/SIZESTRIPE;
				targetangle = currentangle + vectpad;
			else
				targetangle = currentangle + vectstripe*(offsetshifted/nblines-SIZESTRIPE)/SIZESTRIPE; %consider the offset from 0 to 1
			end;
			firstpass = 0;
		else
			targetangle  = currentangle + vectpad;
		end;
		if targetangle >= angle2
			targetangle = angle2;
			cont = 0;
		end;
		hh = circle( centreX, centreY, radius, 0, color, currentangle, targetangle, 0, thickness, ceil((targetangle-currentangle)/4*radius) );
		h = [h hh];

		currentangle = targetangle+vectstripe;
		if currentangle >= angle2
			cont = 0;
		end;
	end;
else
	h = circle( centreX, centreY, radius, 0, color, angle1, angle2, 0, thickness, segments );
end;

if middle
	x = centreX + cos((angle1+angle2)/2/180*pi)*radius;
	y = centreY + sin((angle1+angle2)/2/180*pi)*radius;
	hold on; h = plot( x, y, '*k', 'markersize', thickness/2);
end;

return;

clf;
for i=9:31
	fprintf('******************************* %d \n', mod(i/10, 1));
	circpatch( [ 0+i 1+i ] , [ 10 20 ], 0.5, 'b', 5, 10, mod(i/10, 1),1);
end;

author: PatternRecognition

repo: OpenBMI


text:
function Mliimg = Image_DeSpeckle(mlistack,SHP)

%   Inputs:
%   - mlistack: A height by width by page (real) matrix,e.g., SAR single-look
%               intensity series
%   - SHP:      See script "SHP_SelPoint.m" for details
%   Outputs:
%   - Mliimg:   filtered intensity images
%

if nargin < 1
    help DeSpeckling
    return
end

tic;
[nlines,nwidths,npages]=size(mlistack);
Mliimg = mlistack;

CalWin =SHP.CalWin;
RadiusRow=(CalWin(1)-1)/2;
RadiusCol=(CalWin(2)-1)/2;
mlistack = padarray(mlistack,[RadiusRow RadiusCol],'symmetric');

%Despeckling
for ii=1:npages
    temp = mlistack(:,:,ii);
    num=1;
    for jj = 1:nwidths
        for kk= 1:nlines
            x_global  = jj+RadiusCol;
            y_global  = kk+RadiusRow;
            MliValue  = temp(y_global-RadiusRow:y_global+RadiusRow,x_global-RadiusCol:x_global+RadiusCol);
            MliValue  = MliValue(SHP.PixelInd(:,num));
            Mliimg(kk,jj,ii) = mean(MliValue);
            num=num+1;
        end
    end
    fprintf(' ADP. DESPECKLING: %d / %d is finished...\n',ii,npages);
end

t=toc;
disp(['DeSpeckling operation completed in ',num2str(t/60),' minute(s).']);
disp('Done!');

author: DinhHoTongMinh

repo: TomoSAR


text:
function d=v_distchpf(pf1,pf2,mode)
%V_DISTCHPF calculates the cosh spectral distance between power spectra D=(PF1,PF2,MODE)
%
% Inputs: PF1,PF2     Power spectra to be compared. Each row represents a power spectrum: the first
%                     and last columns represent the DC and Nyquist terms respectively.
%                     PF1 and PF2 must have the same number of columns.
%
%         MODE        Character string selecting the following options:
%                         'x'  Calculate the full distance matrix from every row of PF1 to every row of PF2
%                         'd'  Calculate only the distance between corresponding rows of PF1 and PF2
%                              The default is 'd' if PF1 and PF2 have the same number of rows otherwise 'x'.
%
% Output: D           If MODE='d' then D is a column vector with the same number of rows as the shorter of PF1 and PF2.
%                     If MODE='x' then D is a matrix with the same number of rows as PF1 and the same number of columns as PF2'.
%
% The COSH spectral distance is the average over +ve and -ve frequency of
%
%                     cosh(log(p1/p2))-1   =   (p1-p2)^2/(2p1*p2)   =   (p1/p2 + p2/p1)/2 - 1
%
% The COSH distance is a symmetrical version of the Itakura-Saito distance: v_distchpf(x,y)=(v_distispf(x,y)+v_distispf(y,x))/2

% The Cosh distance can also be calculated directly from AR coefficients; providing np is large
% enough, the values of d0 and d1 in the following will be very similar:
%
%         np=255; d0=v_distchar(ar1,ar2); d1=v_distchpf(v_lpcar2pf(ar1,np),v_lpcar2pf(ar2,np))
%

% Ref: A.H.Gray Jr and J.D.Markel, "Distance measures for speech processing", IEEE ASSP-24(5): 380-391, Oct 1976
%      L. Rabiner abd B-H Juang, "Fundamentals of Speech Recognition", Section 4.5, Prentice-Hall 1993, ISBN 0-13-015157-2

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

[nf1,p2]=size(pf1);
p1=p2-1;
nf2=size(pf2,1);
if nargin<3 | isempty(mode) mode='0'; end
if any(mode=='d') | (mode~='x' & nf1==nf2)
   nx=min(nf1,nf2);
   r=pf1(1:nx,:)./pf2(1:nx,:);
   q=r+r.^(-1);
   d=(2*sum(q(:,2:p1),2)+q(:,1)+q(:,p2))/(4*p1)-1;
else
   r=permute(pf1(:,:,ones(1,nf2)),[1 3 2])./permute(pf2(:,:,ones(1,nf1)),[3 1 2]);
   q=r+r.^(-1);
   d=(2*sum(q(:,:,2:p1),3)+q(:,:,1)+q(:,:,p2))/(4*p1)-1;
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function tests = test_ft_preproc_derivative

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_preproc_derivative

if nargout
  % assume that this is called by RUNTESTS
  tests = functiontests(localfunctions);
else
  % assume that this is called from the command line
  fn = localfunctions;
  for i=1:numel(fn)
    feval(fn{i});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function testOptions(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nchan   = 8;
nsample = 1000;
dat     = randn(nchan, nsample) + 1;

result = {};
result{end+1} = ft_preproc_derivative(dat, 1);
result{end+1} = ft_preproc_derivative(dat, 2);
result{end+1} = ft_preproc_derivative(dat, 4);
result{end+1} = ft_preproc_derivative(dat, 8);

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequal(result{i}, result{j}), 'the results %d and %d should not be equal', i, j);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function I = integral(f, varargin)
%INTEGRAL   Line integral of a CHEBFUN3 over a parametric curve.
%   I = INTEGRAL(F, G), returns the integral of the CHEBFUN3 object F
%   along a parametric curve defined by the Inf x 3 quasimatrix G. Columns
%   of G represent parametrization of the 3D curve.
%
%   I = INTEGRAL(F), returns the triple definite integral of the CHEBFUN3
%   object F over its domain of definition.
%
% See also CHEBFUN3/INTEGRAL2, CHEBFUN3/INTEGRAL3, CHEBFUN3/SUM,
% CHEBFUN3/SUM2, and CHEBFUN3/SUM3.

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

% Empty check:
if ( isempty(f) )
    I = [];
    return
end

if ( nargin == 1 )                       % Another way to do sum3(f)
    % Triple definite integral:
    I = sum3(f);

else
    if ( (isa(varargin{1}, 'chebfun')) && size(varargin{1}, 2) == 3 )
        % Line integral over an Inf x 3 quasimatrix.
        % Get curve:
        curve = varargin{1};
        xCurve = curve(:, 1);
        yCurve = curve(:, 2);
        zCurve = curve(:, 3);
        diffC = diff(curve);
        ds_squared = diffC(:, 1).^2 + diffC(:, 2).^2 + diffC(:, 3).^2;
        I = sum(feval(f, xCurve, yCurve, zCurve) .* sqrt(ds_squared), ...
            curve.domain);
    else
        error('CHEBFUN:CHEBFUN3:integral:badInputs',...
        'Unrecognised input arguments.');
    end
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function y = m_fp_amg(x_it,aparams,qparams)
%m_fp_amg     AMG iterated pressure convection-diffusion preconditioner
%   y = m_fp_amg(x_it,aparams,qparams);
%   input
%          x_it         operand for preconditioning operator
%          aparams      structure defining coefficient matrix
%          mparams      structure defining preconditioning matrix
%   output
%          y            result of preconditioning operation
%
%   calls FEMLAB function amgsol
%   Global array variables AMGLOBA and AMGLOBF define the
%   amg data for amgsol for Ap and F respectively
%
%   IFISS function: DJS; 23 April 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage

global AMGLOBA
global AMGLOBF

nv = length(aparams.F);
nu = nv/2;
np = size(aparams.B,1);

rv=x_it(1:nv); rp=x_it(nv+1:nv+np);

%% pressure solve
nv = length(aparams.F);
nu = nv/2;
np = size(aparams.B,1);

rv=x_it(1:nv); rp=x_it(nv+1:nv+np);
%% pressure solve
if qparams.domain==1,
   n_null = qparams.n_null;
   minor = [1:n_null-1,n_null+1:np]';
   yp = zeros(np,1);
   yp(minor) =  amgsol(rp(minor),AMGLOBA);
   zp = -qparams.Mp\(qparams.Fp*yp);
else
   zp = -(qparams.Mp)\((qparams.Fp)*((qparams.Ap)\rp));
end

%% velocity solve
rv = rv-(aparams.B')*zp;
zv = amgsol(rv,AMGLOBF);
y = [zv;zp];


text:
function out = mirror_extend(in,bx,by)
% Pad array with mirror reflections of itself.
% bx and by specify the amount of padding to add.
% ***********************************************

[h,w] = size(in);

%First flip up and down
u = flipud(in(2:1+bx,:));
d = flipud(in(h-bx:h-1,:));

in2 = [u' in' d']';

%Next flip left and right
l = fliplr(in2(:, 2:1+by));
r = fliplr(in2(:,w-by:w-1));

%set the 'mirrored' image to out.
out = [l in2 r];

return

%test
% A = [1 2 3 4;5 6 7 8;9 10 11 12]
% B = mirror_extend(A,2,2)

author: thfylsty


text:
%% Copyright (C) 2014-2016, 2019 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym isequaln (@var{f}, @var{g})
%% @defmethodx @@sym isequaln (@var{f}, @var{g}, @dots{})
%% Test if contents of arrays are equal, even with nan.
%%
%% Here NaN's are considered equal:
%% @example
%% @group
%% syms x
%% snan = sym(nan);
%% isequaln([1 snan x], [1 snan x])
%%   @result{} 1
%% @end group
%% @end example
%% To get the usual NaN convention, @pxref{isequal}.
%%
%% @seealso{@@sym/isequal, @@sym/logical, @@sym/isAlways, @@sym/eq}
%% @end defmethod

function t = isequaln(x, y, varargin)

  if (nargin < 2)
    print_usage ();
  end

  % isequal does not care about type, but if you wanted it to...
  %if ( ~ ( isa (x, 'sym') && isa (y, 'sym')))
  %  t = false;
  %  return
  %end

  %% some special cases
  if ~(is_same_shape(x, y))
    t = false;
    return
  end

  % In symy, nan == nan is true by structural (not mathematical)
  % equivalence, so we don't need to detect it ourselves.
  % Sympy's == returns a scalar for arrays, no special case.

  cmd = 'return (_ins[0] == _ins[1],)';

  t = pycall_sympy__ (cmd, sym(x), sym(y));

  if (~ islogical(t))
    error('nonboolean return from python');
  end

  if (nargin >= 3)
    t = t && isequaln(x, varargin{:});
  end

end

%!test
%! a = sym([1 2]);
%! b = a;
%! assert (isequaln (a, b))
%! b(1) = 42;
%! assert (~isequaln (a, b))

%!test
%! a = sym([1 2; 3 4]);
%! b = a;
%! assert (isequaln (a, b))
%! b(1) = 42;
%! assert (~isequaln (a, b))

%!test
%! a = sym([nan; 2]);
%! b = a;
%! assert (isequaln (a, b))

%!test
%! a = sym([nan 2; 3 4]);
%! b = a;
%! assert (isequaln (a, b))

%!test
%! % more than two arrays
%! a = sym([nan 2 3]);
%! b = a;
%! c = a;
%! assert (isequaln (a, b, c))
%! c(1) = 42;
%! assert (~isequaln (a, b, c))

repo: octsympy


text:
function [cl, cp, cs] = dtiComputeWestinShapes(eigVal, denominator)
%
% [cl, cp, cs] = dtiComputeWestinShapes(eigVal, [denominator='lsum'])
% Computes Westin's shape indices
%'denominator' is a parameter allowing to choose an appropriate definition
%for Westin shapes.
%
% denominator='lsum' (default):
%    cl = (lambda_1 - lambda_2) / (lambda_1 + lambda_2 + lambda_3)
%    cp = 2 * (lambda_2 - lambda_3) / (lambda_1 + lambda_2 + lambda_3)
%    cs = 3 * lambda_3 / (lambda_1 + lambda_2 + lambda_3)
%
% This was the originial formulation described in:
%
% C-F. Westin, S. Peled, H. Gubbjartsson, R. Kikinis, and F.A. Jolesz.
% Geometrical diffusion measures for MRI from tensor basis analysis.
% In Proceedings 5th Annual ISMRM, 1997.
%
% denominator='l1':
%
% In later work, Westin et. al. adopted a simpler normalization
% formulation (e.g., see Westin et. al. 2002 Med. Image Anal.; PMID:
% 12044998) where the constants are dropped and the denominator is simply
% lambda_1. This definition produces cl+cp+cs=1 which is convenient for
% tensor shape representation in barycentric coordinates.
% In practice, the two methods produce very similar maps.
% cl = (lambda_1 - lambda_2) / (lambda_1)
% cp = 2 * (lambda_2 - lambda_3) / (lambda_1)
% cs = 3 * lambda_3 / (lambda_1)
%
% eigVal: XxYxZx3 or nx3xN array of tensor eigenvalues. Or, you can pass
% the dt6 array (XxYxZx6) and we'll compute the eigenvalues for you.
%
%
% HISTORY:
% 2004.07.22 RFD (bobd@stanford.edu) & ASH wrote it.
% 2007.01.02 SHC made changes to allow nx3xN format.
% 2009.01.05 EIR added option of using a more recent definition of Westin shapes

if ~exist('denominator', 'var')
    denominator='lsum';
end

if size(eigVal,4)==6
    [eigVec, eigVal] = dtiEig(eigVal);
    clear eigVec;
end

% Check inputs
switch ndims(eigVal)
    case 2,
        Ind    = 1; % Data in indexed nx6 format
        eigVal = shiftdim(eigVal, -2);
    case 3,
        if size(eigVal,2)==6 && size(eigVal,3)~=6
            Ind    = 1; % Data in indexed nx6xN format
            eigVal = shiftdim(eigVal, -2);
        else
            Ind    = 2; % Data in XxYx6 format
            eigVal = shiftdim(eigVal, -1);
        end
    otherwise,
        Ind = 0; % Data in XxYxZx6xN format
end

epsilon = 1e-10;
switch denominator
    case 'lsum'
        denum=sum(eigVal,4);

    case 'l1'
        denum=eigVal(:, :, :, 1);
    otherwise fprintf('Wrong denominator option');
end

nz = denum>epsilon;

% Avoid divide-by-zero (we'll replace these values with zeros below).
denum(~nz) = 1;

cl = (eigVal(:,:,:,1)-eigVal(:,:,:,2))./denum;
cp = 2*(eigVal(:,:,:,2)-eigVal(:,:,:,3))./denum;
cs = 3*eigVal(:,:,:,3)./denum;

if strcmp(denominator,'l1')
    cp=cp./2;
    cs=cs./3;
end

cl(~nz) = 0; cp(~nz) = 0; cs(~nz) = 0;
cl(cl>1.0) = 1.0;
cp(cp>1.0) = 1.0;
cs(cs>1.0) = 1.0;
cl(cl<0.0) = 0.0;
cp(cp<0.0) = 0.0;
cs(cs<0.0) = 0.0;

% Adjust output
switch Ind
    case 1,
        cl = shiftdim(cl, 2);
        cp = shiftdim(cp, 2);
        cs = shiftdim(cs, 2);
    case 2,
        cl = shiftdim(cl, 1);
        cp = shiftdim(cp, 1);
        cs = shiftdim(cs, 1);
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function title = p05_title ( )

%*****************************************************************************80
%
%% P05_TITLE returns a title for problem 5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'The piecewise linear function';

  return
end


text:
function ao=lpcrf2ao(rf)
%LPCRF2AO Convert reflection coefficients to area ratios AO=(RF)
% ao(k) is the ratio between the k'th and the (k+1)'th segment

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

ao =(1-rf)./(1+rf);

author: decouples


text:
function [A,B] = vrep2hrep(pts);
%
% Given a set of points, computes an H-representation of its convex hull
%
%  A x <= B
%
% Calls Komei Fukuda's CDD (floating point version)

% This file is part of SOSTOOLS - Sum of Squares Toolbox ver 3.00.
%
% Copyright (C)2002, 2004, 2013  A. Papachristodoulou (1), J. Anderson (1),
%                                G. Valmorbida (1), S. Prajna (2),
%                                P. Seiler (3), P. A. Parrilo (4)
% (1) Department of Engineering Science, University of Oxford, Oxford, U.K.
% (2) Control and Dynamical Systems - California Institute of Technology,
%     Pasadena, CA 91125, USA.
% (3) Aerospace and Engineering Mechanics Department, University of
%     Minnesota, Minneapolis, MN 55455-0153, USA.
% (4) Laboratory for Information and Decision Systems, M.I.T.,
%     Massachusetts, MA 02139-4307
%
% Send bug reports and feedback to: sostools@cds.caltech.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, see <http://www.gnu.org/licenses/>.

% PP 22/10/2002

[npts,dim]=size(pts);

cdd = cddpath;

filename = 'temp.ext' ;

fid = fopen(filename,'w');

fprintf(fid,'*\n');
fprintf(fid,'* Generated automatically by vrep2hrep.m\n');
fprintf(fid,'*\n');
fprintf(fid,'V-representation\n');
fprintf(fid,'begin\n');
fprintf(fid,'%d %d real\n',npts,dim+1);

fmt = ['1 ' repmat('%.15E ',1,dim) '\n'];

% Print the data
fprintf(fid,fmt,pts');

fprintf(fid,'end\n');
fprintf(fid,'stdout_off\n');
fclose(fid);

% Now, run CDD
cmd = [cdd ' ' filename];

vv = version;
if str2num(vv(1))>5
  [dummy1,dummy2] = system(cmd);
else
  if ~isunix
    [dummy1,dummy2] = dos(cmd);
  else
    [dummy1,dummy2] = unix(cmd);
  end;
end;

% Read the output

filename = 'temp.ine' ;

fid = fopen(filename,'r');

% Skip everything, until 'begin'
tline = [];
while strcmp(tline,'begin') == 0;
 tline = fgetl(fid) ;
end

nineqs = fscanf(fid,'%d',1);
dims   = fscanf(fid,'%d',1);
dummy2 = fscanf(fid,'%s',1);

H = fscanf(fid,'%E',[dims,nineqs]);
fclose(fid);

% Get rid of the files
delete('temp.ext');
delete('temp.ine');

H = H' ;

B = H(:,1);
A = -H(:,2:dims);

author: yu-jiang

repo: radpbook

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


text:
function [band_offsets, quantif] = ReadS2RadiometricOffset(FilesInfo)
% This funciton is to read the Sentinel-2's radiometric offset for the
% version with Basline 4.00
%
% Return Empty means old version data
%
% Created on 16/12/2021

    filepath_meta = dir(fullfile(FilesInfo.InspireXML, 'M*L1C.xml')); % which is same location as inspire xml
    S = xml2struct(fullfile(FilesInfo.InspireXML, filepath_meta.name));
    quantif = str2num(S.n1_colon_Level_dash_1C_User_Product.n1_colon_General_Info.Product_Image_Characteristics.QUANTIFICATION_VALUE.Text);
    try
        Radiometric_Offset_List = S.n1_colon_Level_dash_1C_User_Product.n1_colon_General_Info.Product_Image_Characteristics.Radiometric_Offset_List.RADIO_ADD_OFFSET;

        band_offsets = zeros(length(Radiometric_Offset_List), 2);
        for i = 1: length(Radiometric_Offset_List)
            offset_value = str2num(Radiometric_Offset_List{1, i}.Text);
            band_id = str2num(Radiometric_Offset_List{1, i}.Attributes.band_id);
            band_offsets(i,1) = band_id;
            band_offsets(i,2) = offset_value;
        end
    catch
        band_offsets = [];
    end
end

author: GERSL

repo: Fmask

save_path: github-repos/MATLAB/GERSL-Fmask


% pivot_direction.m // Jon Lee
% view the direction zbar corresponding to increasing x_eta(j).
% syntax is pivot_direction(j)

function [zbar] = pivot_direction(j)
global beta eta m n Abar_eta
if (j > n-m)
    display('error: not so many directions')
    return
end
zbar = zeros(n,1);
zbar(beta) = -Abar_eta(:,j);
zbar(eta(j))=1;
zbar;
end


text:
function [requ_clab, W, neighbor_clab]= procutil_getClabForLaplacian(dat, varargin)
%GETCLABFORLAPLACIAN - Get channels which are required for Laplacian filtering
%
%Synopsis:
%   CLAB= procutil_getClabForLaplacian(DAT, <OPT>)
%   CLAB= procutil_getClabForLaplacian(DAT, CLAB, <OPT>)
%
%Arguments:
%   DAT:    STRUCT   - Continuous or epoched data, or directly the full
%                      channel set from which the clabs for laplacian should
%                      be extracted
%   CLAB:   CELL     - Names of channels that should be filtered by
%                      proc_laplacian. Accepts wildcards.
%   OPT: 	PROPLIST - struct or property/value list of optional properties:
%     'FilterType': CHAR (default 'small') - Name of the type of Laplace
%                      filter that should be used. Simply add your filter
%                      to the subfunction 'getLaplaceFilter'.
%     'CLab': CELL     - Equal to CLAB as a second parameter.
%     'IgnoreCLab': CELL (default {'E*'}) - Names of the channels that are
%                      to be left out of the calculation.
%     'GridFcn': FUNC
%     'RequireCompleteNeighborhood: BOOL (default 1) - Enforce that a
%                      Laplace channel can only be returned if all the
%                      necessary neighbors are available.
%
%Returns:
%   REQU_CLAB:         cell array of channel labels, necessary for the
%                      given Laplace channels
%   W:                 filter matrix that can be used, e.g. in
%                      proc_linearDerivation, to calculate the Laplace
%                      channels
%   NEIGHBOR_CLAB:     Channels around the Laplace channels that are used
%                      for referencing
%
%Examples:
%   procutil_getClabForLaplacian(dat, {'C3', 'C4'});
%   procutil_getClabForLaplacian(dat, 'FilterType', 'large', 'CLab', {'C3', 'C4'});

%See:
% proc_laplacian
%
% Benjamin Blankertz
% Martijn Schreuder 06/12 - Updated the help documentation

props= {'CLab'                         '*'                    'CHAR|CELL{CHAR}'
        'IgnoreCLab'                   {'E*'}                 'CHAR|CELL{CHAR}'
        'GridFcn'                      @util_gridForLaplacian 'FUNC'
        'FilterType'                   'small'                'CHAR(small large horizontal vertical bip_to_anterior bip_to_posterior bip_to_left bip_to_right diagonal diagonal_small)'
        'RequireCompleteNeighborhood'  1                      'BOOL'};

if mod(nargin,2)==0,
  misc_checkType(varargin{1}, 'CHAR|CELL{CHAR}', 'CLAB');
  opt= opt_proplistToStruct(varargin{2:end});
  opt.CLab= varargin{1};
else
  opt= opt_proplistToStruct(varargin{:});
end

[opt,isdefault] = opt_setDefaults(opt, props);
opt_checkProplist(opt, props);
misc_checkType(dat, 'CELL{CHAR}|STRUCT(clab)');
tcstate= bbci_typechecking('off');

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

if isstruct(dat)
  clab= dat.clab;
else
  clab= dat;
end

laplace= [];
laplace.grid= opt.GridFcn();
laplace.filter= procutil_lapGetLaplacianFilter(opt.FilterType);

rc= util_chanind(clab, cat(2, {'not'}, opt.IgnoreCLab));
nOrigChans= length(clab);
pos= zeros(2, nOrigChans);
for ic= 1:nOrigChans,
  pos(:,ic)= procutil_lapGetCoordinates(clab{ic}, laplace.grid);
end
pos(:,setdiff(1:nOrigChans,rc,'legacy'))= inf;

idx_tbf= util_chanind(clab, opt.CLab);
W= zeros(length(clab), length(idx_tbf));
lc= 0;
requ_clab= {};
neighbor_clab= cell(length(idx_tbf), 1);
for ci= 1:length(idx_tbf),
  cc= idx_tbf(ci);
  refChans= [];
  nRefs= size(laplace.filter,2);
  for ir= 1:nRefs,
    ri= find( pos(1,:)==pos(1,cc)+laplace.filter(1,ir) & ...
      pos(2,:)==pos(2,cc)+laplace.filter(2,ir) );
    refChans= [refChans ri];
  end
  if length(refChans)==nRefs | ~opt.RequireCompleteNeighborhood,
    lc= lc+1;
    W(cc,lc)= 1;
    if ~isempty(refChans),
      W(refChans,lc)= -1/length(refChans);
    end
    requ_clab= unique(cat(2, requ_clab, clab([cc refChans])),'legacy');
    neighbor_clab{ci}= clab(refChans);
  end
end
W= W(util_chanind(clab, requ_clab),:);

bbci_typechecking(tcstate);

author: bbci


text:
function hermite_cubic_test02 ( )

%*****************************************************************************80
%
%% HERMITE_CUBIC_TEST02 tests HERMITE_CUBIC_VALUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_CUBIC_TEST02:\n' );
  fprintf ( 1, '  HERMITE_CUBIC_VALUE evaluates a Hermite cubic polynomial.\n' );
  fprintf ( 1, '  Try out data from a cubic function:\n' );
  fprintf ( 1, '  on [0,10] and [-1.0,1.0] and [0.5,0.75]\n' );

  for x_interval = 1 : 3

    if ( x_interval == 1 )
      x1 = 0.0;
      x2 = 10.0;
    elseif ( x_interval == 2 )
      x1 = -1.0;
      x2 = +1.0;
    elseif ( x_interval == 3 )
      x1 = 0.5;
      x2 = 0.75;
    end

    [ f1, d1, s1, t1 ] = cubic_value ( x1 );
    [ f2, d2, s2, t2 ] = cubic_value ( x2 );

    fprintf ( 1, '\n' );
    fprintf ( 1, '    J      X           F           D           S           T\n' );

    for j = - 3 : 12
      x = ( ( 10 - j ) * x1   ...
          +        j   * x2 ) ...
            / 10.0;

      [ f, d, s, t ] = hermite_cubic_value ( x1, f1, d1, x2, f2, d2, 1, x );
      [ fc, dc, sc, tc ] = cubic_value ( x );

      fprintf ( 1, '\n' );
      if ( j == 0 )
        fprintf ( 1, '*Data  %10f  %10f  %10f\n', x1, f1, d1 );
      end
      fprintf ( 1, 'Exact  %10f  %10f  %10f  %10f  %10f\n',    x,  fc, dc, sc, tc );
      fprintf ( 1, '  %3d  %10f  %10f  %10f  %10f  %10f\n', j, x,  f,  d,  s,  t );
      if ( j == 10 )
        fprintf ( 1, '*Data  %10f  %10f  %10f\n', x2, f2, d2 );
      end

    end

  end

  return
end

author: johannesgerer


text:
function u = forward ( e, r, sigma, t1, Nx, Nt, L )

%*****************************************************************************80
%
%% FORWARD uses the forward difference method to value a European call option.
%
%  Modified:
%
%    28 March 2005
%
%  Author:
%
%    Original MATLAB version by Desmond Higham
%
%  Reference:
%
%    Desmond Higham,
%    Black-Scholes for Scientific Computing Students,
%    Computing in Science and Engineering,
%    November/December 2004, Volume 6, Number 6, pages 72-79.
%
%  Parameters:
%
%    Input, real E, the exercise price.
%
%    Input, real R, the interest rate.
%
%    Input, real SIGMA, the volatility of the asset.
%
%    Input, real T1, the expiry date.
%
%    Input, integer NX, the number of "space" steps used to divide the
%    interval [0,L].
%
%    Input, integer NT, the number of time steps.
%
%    Input, real L, the maximum value of S to consider.
%
%    Output, real U(NX-1,NT+1), the value of the European call option.
%
  k = t1 / Nt;
  h = L / Nx;
  TR1 = diag ( ones ( Nx-2,1),1) - diag ( ones(Nx-2,1), -1 );
  TR2 = -2 * eye(Nx-1,Nx-1) + diag(ones(Nx-2,1),1) + diag(ones(Nx-2,1),-1);
  mvec = [1:Nx-1];
  D1 = diag ( mvec );
  D2 = diag ( mvec.^2 );
  Aftcs = ( 1.0 - r * k ) * eye ( Nx-1,Nx-1 ) ...
        + 0.5 * k * sigma^2 * D2 * TR2 ...
        + 0.5 * k * r * D1 * TR1;

  u = zeros ( Nx-1,Nt+1);
  uzero = max ( [h:h:L-h]'-e, 0.0 );
  u(:,1) = uzero;
  p = zeros(Nx-1,1);

  for i = 1 : Nt
    tau = ( i - 1 ) * k;
    p(end) = 0.5 * k * (Nx-1) * ((sigma^2)*(Nx-1)+r) * ( L - e * exp ( - r * tau ) );
    u(:,i+1) = Aftcs * u(:,i) + p;
  end

  return
end


text:
 function h = ir_str2func(string)
%function h = ir_str2func(string)
%|
%| wrapper around str2func() for octave / matlab compatability
%|
%| Background: matlab is phasing out inline
%| http://www.mathworks.com/help/matlab/ref/inline.html
%|
%| So given a string like 'x^2' the way to make it function is
%| str2func('@(x) x^2')
%|
%| Unfortunately, octave 3.8.2 does not support that syntax,
%| so this is the work around.

if nargin < 1, ir_usage, end
if streq(string, 'test') ir_str2func_test, return, end

if ir_is_octave
	h = eval(string);
else
	h = str2func(string);
end

% ir_str2func_test()
% test routine (that is longer than the actual function)
function ir_str2func_test
s = '@(x) x.^2';
f = ir_str2func(s);
g = @(x) x.^2;
jf_equal(f, g)
x = 1:5;
y = f(x);
jf_equal(y, x.^2)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


function f = draw_3Dskeleton_MSRA(d2img,coord_world)

    cubicSz = 150;
    coord_world = squeeze(coord_world);

    x = coord_world(1,:);
    y = coord_world(2,:);
    z = coord_world(3,:);
    coord_world(1,:) = -z;
    coord_world(2,:) = x;
    coord_world(3,:) = y;
    line_width = 8;

    f = figure('Position',[100 100 600 450]);
    set(f, 'visible', 'off');
    hold on;
    grid on;

    plot3([coord_world(1,1),coord_world(1,18)],[coord_world(2,1),coord_world(2,18)],[coord_world(3,1),coord_world(3,18)],'Color',[255/255,153/255,153/255],'LineWidth',line_width) %wrist to thumb mcp
    plot3([coord_world(1,18),coord_world(1,19)],[coord_world(2,18),coord_world(2,19)],[coord_world(3,18),coord_world(3,19)],'Color',[255/255,102/255,102/255],'LineWidth',line_width) %thumb mcp to thumb pip
    plot3([coord_world(1,19),coord_world(1,20)],[coord_world(2,19),coord_world(2,20)],[coord_world(3,19),coord_world(3,20)],'Color',[255/255,51/255,51/255],'LineWidth',line_width) %thumb pip to thumb dip
    plot3([coord_world(1,20),coord_world(1,21)],[coord_world(2,20),coord_world(2,21)],[coord_world(3,20),coord_world(3,21)],'Color',[255/255,0/255,0/255],'LineWidth',line_width) %thumb dip to thumb tip

    plot3([coord_world(1,1),coord_world(1,2)],[coord_world(2,1),coord_world(2,2)],[coord_world(3,1),coord_world(3,2)],'Color',[153/255,255/255,153/255],'LineWidth',line_width) %wrist to index mcp
    plot3([coord_world(1,2),coord_world(1,3)],[coord_world(2,2),coord_world(2,3)],[coord_world(3,2),coord_world(3,3)],'Color',[102/255,255/255,102/255],'LineWidth',line_width) %index mcp to index pip
    plot3([coord_world(1,3),coord_world(1,4)],[coord_world(2,3),coord_world(2,4)],[coord_world(3,3),coord_world(3,4)],'Color',[51/255,255/255,51/255],'LineWidth',line_width) %index pip to index dip
    plot3([coord_world(1,4),coord_world(1,5)],[coord_world(2,4),coord_world(2,5)],[coord_world(3,4),coord_world(3,5)],'Color',[0/255,255/255,0/255],'LineWidth',line_width) %index dip to index tip

    plot3([coord_world(1,1),coord_world(1,6)],[coord_world(2,1),coord_world(2,6)],[coord_world(3,1),coord_world(3,6)],'Color',[255/255,204/255,153/255],'LineWidth',line_width) %wrist to middle mcp
    plot3([coord_world(1,6),coord_world(1,7)],[coord_world(2,6),coord_world(2,7)],[coord_world(3,6),coord_world(3,7)],'Color',[255/255,178/255,102/255],'LineWidth',line_width) %middle mcp to middle pip
    plot3([coord_world(1,7),coord_world(1,8)],[coord_world(2,7),coord_world(2,8)],[coord_world(3,7),coord_world(3,8)],'Color',[255/255,153/255,51/255],'LineWidth',line_width) %middle pip to middle dip
    plot3([coord_world(1,8),coord_world(1,9)],[coord_world(2,8),coord_world(2,9)],[coord_world(3,8),coord_world(3,9)],'Color',[255/255,128/255,0/255],'LineWidth',line_width) %middle dip to middle tip

    plot3([coord_world(1,1),coord_world(1,10)],[coord_world(2,1),coord_world(2,10)],[coord_world(3,1),coord_world(3,10)],'Color',[153/255,204/255,255/255],'LineWidth',line_width) %wrist to ring mcp
    plot3([coord_world(1,10),coord_world(1,11)],[coord_world(2,10),coord_world(2,11)],[coord_world(3,10),coord_world(3,11)],'Color',[102/255,178/255,255/255],'LineWidth',line_width) %ring mcp to ring pip
    plot3([coord_world(1,11),coord_world(1,12)],[coord_world(2,11),coord_world(2,12)],[coord_world(3,11),coord_world(3,12)],'Color',[51/255,153/255,255/255],'LineWidth',line_width) %ring pip to ring dip
    plot3([coord_world(1,12),coord_world(1,13)],[coord_world(2,12),coord_world(2,13)],[coord_world(3,12),coord_world(3,13)],'Color',[0/255,128/255,255/255],'LineWidth',line_width) %ring dip to ring tip

    plot3([coord_world(1,1),coord_world(1,14)],[coord_world(2,1),coord_world(2,14)],[coord_world(3,1),coord_world(3,14)],'Color',[255/255,153/255,255/255],'LineWidth',line_width) %wrist to little mcp
    plot3([coord_world(1,14),coord_world(1,15)],[coord_world(2,14),coord_world(2,15)],[coord_world(3,14),coord_world(3,15)],'Color',[255/255,102/255,255/255],'LineWidth',line_width) %little mcp to little pip
    plot3([coord_world(1,15),coord_world(1,16)],[coord_world(2,15),coord_world(2,16)],[coord_world(3,15),coord_world(3,16)],'Color',[255/255,51/255,255/255],'LineWidth',line_width) %little pip to little dip
    plot3([coord_world(1,16),coord_world(1,17)],[coord_world(2,16),coord_world(2,17)],[coord_world(3,16),coord_world(3,17)],'Color',[255/255,0/255,255/255],'LineWidth',line_width) %little dip to little tip

    colorList = [
    230/255 230/255 0/255;

    153/255 255/255 153/255;
    102/255 255/255 102/255;
    51/255 255/255 51/255;
    0/255 255/255 0/255;

    255/255 204/255 153/255;
    255/255 178/255 102/255;
    255/255 153/255 51/255;
    255/255 128/255 0/255;

    153/255 204/255 255/255;
    102/255 178/255 255/255;
    51/255 153/255 255/255;
    0/255 128/255 255/255;

    255/255 153/255 255/255;
    255/255 102/255 255/255;
    255/255 51/255 255/255;
    255/255 0/255 255/255;

    255/255 153/255 153/255;
    255/255 102/255 102/255;
    255/255 51/255 51/255;
    255/255 0/255 0/255;
    ];
    scatter3(coord_world(1,:),coord_world(2,:),coord_world(3,:),150,colorList,'filled');

    set(gca, 'color', [255/255 255/255 255/255])

    set(gca,'XTickLabel',[]);
    set(gca,'YTickLabel',[]);
    set(gca,'ZTickLabel',[]);

    xcenter = (min(coord_world(1,:)) + max(coord_world(1,:)))/2;
    ycenter = (min(coord_world(2,:)) + max(coord_world(2,:)))/2;
    zcenter = (min(coord_world(3,:)) + max(coord_world(3,:)))/2;

    xmin = xcenter - 3*cubicSz;
    xmax = max(coord_world(1,:));
    ymin = min(ycenter-cubicSz/2, min(coord_world(2,:)));
    ymax = max(coord_world(2,:))+30;
    zmin = min(zcenter-cubicSz/2, min(coord_world(3,:)));
    zmax = max(zcenter+cubicSz/2, max(coord_world(3,:)));

    xlim([xmin xmax]);
    ylim([ymin ymax]);
    zlim([zmin zmax]);

    h_img = surf([xmin;xmin],[ymin ymax;ymin ymax],[zmax zmax;zmin zmin],'CData',d2img,'FaceColor','texturemap');
    set(h_img,'edgecolor',[96/255 96/255 96/255]);

    view(62,7);

end


text:
function [metabolite_structure,IDsAdded] = parseCHOmineWebpage(metabolite_structure,startSearch,endSearch)
% try to guess chomine abbreviation based on VMH ID

annotationSource = 'CHOmine website matching';
annotationType = 'automatic';

Mets = fieldnames(metabolite_structure);
if ~exist('startSearch','var')
    startSearch = 1;
end
if ~exist('endSearch','var')
    endSearch = length(Mets);
end

a = 1;
IDsAdded = '';

for i = startSearch : endSearch
    if ~isempty(metabolite_structure.(Mets{i}).VMHId) && isempty(find(isnan(metabolite_structure.(Mets{i}).VMHId),1))

        if (isempty(metabolite_structure.(Mets{i}).chodb_id) || ~isempty(find(isnan(metabolite_structure.(Mets{i}).chodb_id),1)))
            try
                url = strcat('https://chomine.boku.ac.at/chomine/portal.do?externalid=',metabolite_structure.(Mets{i}).VMHId,'&class=Species');
                syst = urlread(url);
                if ~contains(syst,'No matches found')
                    metabolite_structure.(Mets{i}).chodb_id = Mets{i};
                    metabolite_structure.(Mets{i}).chodb_id_source =   [annotationSource,':',annotationType,':',datestr(now)];
                    IDsAdded{a,1} = Mets{i};
                    IDsAdded{a,2} = 'chodb_id';
                    IDsAdded{a,3} = char(idNew);
                    a = a + 1;
                end
            end
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function Vo = camview(varargin)
%CAMVIEW  Records or sets the viewpoint of the current axes
%
%   Vo = camview([hAx], [Vi])
%
% Records or sets the viewpoint of the current axes, including projection
% type. This is useful for giving multiple axes (with the same coordinate
% frame) the same viewpoint, or resetting the viewpoint to a known point.
% A data aspect ratio of [1 1 1] is assumed.
%
% The view is recorded and set using a 4x4 matrix. The upper 3x4 quadrant
% of this matrix is a projection matrix from the axes coordinate frame to
% the camera coordinate frame (the view frustrum coordinates range from -1
% to 1 in x and y directions). Padding the last row with [0 0 0 1] makes it
% possible to use projection matrices computed externally.
%
% IN:
%    hAx - Handle to the axes in question. Default: gca.
%    Vi - 4x4 matrix defining the viewpoint to set the current axes to. The
%         matrix should defined by calling camview on previous axes.
%
% OUT:
%    Vo - 4x4 matrix specifying the viewpoint of the current axes just
%         prior to the function being called.

% Set default inputs
hAx = gca;
Vi = [];
% Parse the inputs
for a = 1:nargin
    if ishandle(varargin{a})
        hAx = varargin{a};
    else
        Vi = varargin{a};
    end
end

if nargout > 0
    % Get the current viewpoint
    t = get(hAx, 'CameraPosition');
    d = get(hAx, 'CameraTarget') - t;
    K = eye(3);
    K([1 5]) = 1 / tan(get(hAx, 'CameraViewAngle') * pi / 360);
    R(:,3) = d / norm(d);
    R(:,2) = get(hAx, 'CameraUpVector');
    R(:,1) = cross(R(:,3), R(:,2));
    Vo = K * R' * [eye(3) -t'];
    Vo(4,:) = [norm(d) 0 0 strcmp(get(hAx, 'Projection'), 'perspective')];
end
if ~isempty(Vi)
    % Decompose the projection matrix
    st = @(M) M(end:-1:1,end:-1:1)';
    [R, K] = qr(st(Vi(1:3,1:3)));
    K = st(K);
    I = diag(K) < 0;
    K(:,I) = -K(:,I);
    R = st(R);
    R(I,:) = -R(I,:);
    t = (K * R) \ -Vi(1:3,4);
    K = K / K(3,3);

    % Set the current viewpoint
    projection = {'perspective', 'orthographic'};
    set(hAx, 'CameraTarget', t'+R(3,:)*(Vi(4)+(Vi(4)==0)), ...
             'CameraPosition', t, ...
             'CameraUpVector', R(2,:), ...
             'CameraViewAngle', atan(1/K(5))*360/pi, ...
             'Projection', projection{(Vi(16)==0)+1});
end
return

author: mtex-toolbox

repo: mtex

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


text:
function [ n_data, a, b, x, fx ] = von_mises_cdf_values ( n_data )

%*****************************************************************************80
%
%% VON_MISES_CDF_VALUES returns some values of the von Mises CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Kanti Mardia and Peter Jupp,
%    Directional Statistics,
%    Wiley, 2000, QA276.M335
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real A, B, the parameters of the function.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 23;

  a_vec = [ ...
     0.0E+00, ...
     0.0E+00, ...
     0.0E+00, ...
     0.0E+00, ...
     0.0E+00, ...
     0.1E+01, ...
     0.1E+01, ...
     0.1E+01, ...
     0.1E+01, ...
     0.1E+01, ...
     0.1E+01, ...
    -0.2E+01, ...
    -0.1E+01, ...
     0.0E+01, ...
     0.1E+01, ...
     0.2E+01, ...
     0.3E+01, ...
     0.0E+00, ...
     0.0E+00, ...
     0.0E+00, ...
     0.0E+00, ...
     0.0E+00, ...
     0.0E+00 ];

  b_vec = [ ...
     0.1E+01, ...
     0.1E+01, ...
     0.1E+01, ...
     0.1E+01, ...
     0.1E+01, ...
     0.2E+01, ...
     0.2E+01, ...
     0.2E+01, ...
     0.2E+01, ...
     0.2E+01, ...
     0.2E+01, ...
     0.3E+01, ...
     0.3E+01, ...
     0.3E+01, ...
     0.3E+01, ...
     0.3E+01, ...
     0.3E+01, ...
     0.0E+00, ...
     0.1E+01, ...
     0.2E+01, ...
     0.3E+01, ...
     0.4E+01, ...
     0.5E+01 ];

  fx_vec = [ ......
    0.2535089956281180E-01, ...
    0.1097539041177346E+00, ...
    0.5000000000000000E+00, ...
    0.8043381312498558E+00, ...
    0.9417460124555197E+00, ...
    0.5000000000000000E+00, ...
    0.6018204118446155E+00, ...
    0.6959356933122230E+00, ...
    0.7765935901304593E+00, ...
    0.8410725934916615E+00, ...
    0.8895777369550366E+00, ...
    0.9960322705517925E+00, ...
    0.9404336090170247E+00, ...
    0.5000000000000000E+00, ...
    0.5956639098297530E-01, ...
    0.3967729448207649E-02, ...
    0.2321953958111930E-03, ...
    0.6250000000000000E+00, ...
    0.7438406999109122E+00, ...
    0.8369224904294019E+00, ...
    0.8941711407897124E+00, ...
    0.9291058600568743E+00, ...
    0.9514289900655436E+00 ];

  x_vec = [ ......
    -0.2617993977991494E+01, ...
    -0.1570796326794897E+01, ...
     0.0000000000000000E+00, ...
     0.1047197551196598E+01, ...
     0.2094395102393195E+01, ...
     0.1000000000000000E+01, ...
     0.1200000000000000E+01, ...
     0.1400000000000000E+01, ...
     0.1600000000000000E+01, ...
     0.1800000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.0000000000000000E+00, ...
     0.0000000000000000E+00, ...
     0.0000000000000000E+00, ...
     0.0000000000000000E+00, ...
     0.0000000000000000E+00, ...
     0.0000000000000000E+00, ...
     0.7853981633974483E+00, ...
     0.7853981633974483E+00, ...
     0.7853981633974483E+00, ...
     0.7853981633974483E+00, ...
     0.7853981633974483E+00, ...
     0.7853981633974483E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    a = 0.0;
    b = 0.0;
    x = 0.0;
    fx = 0.0;
  else
    a = a_vec(n_data);
    b = b_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

repo: jburkardt-m


text:
function [SE_MR,SE_RZF,SE_MMMSE,SE_ZF,SE_SMMSE] = functionComputeSE_UL(Hhat,C,R,tau_c,tau_p,nbrOfRealizations,M,K,L,p)
%Compute UL SE for different receive combining schemes using Theorem 4.1.
%
%INPUT:
%Hhat              = M x nbrOfRealizations x K x L x L matrix with the MMSE
%                    channel estimates
%C                 = M x M x K x L x L matrix with estimation error
%                    correlation matrices when using MMSE estimation
%R                 = M x M x K x L x L matrix with spatial correlation
%                    matrices
%tau_c             = Length of coherence block
%tau_p             = Length of pilot sequences
%nbrOfRealizations = Number of channel realizations
%M                 = Number of antennas per BS
%K                 = Number of UEs per cell
%L                 = Number of BSs and cells
%p                 = Uplink transmit power per UE (same for everyone)
%
%OUTPUT:
%SE_MR    = K x L matrix where element (k,l) is the uplink SE of UE k in
%           cell l achieved with MR combining
%SE_RZF   = Same as SE_MR but with RZF combining
%SE_MMMSE = Same as SE_MR but with M-MMSE combining
%SE_ZF    = Same as SE_MR but with ZF combining
%SE_SMMSE = Same as SE_MR but with S-MMSE combining
%
%
%This Matlab function was developed to generate simulation results to:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.01 (Last edited: 2020-05-15)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Store identity matrices of different sizes
eyeK = eye(K);
eyeM = eye(M);

%Compute the pre-log factor (normalized with number of channel realizations)
%assuming only uplink transmission
prelogFactor = (tau_c-tau_p)/(tau_c*nbrOfRealizations);

%Prepare to store simulation results
SE_MR = zeros(K,L);

if nargout > 1
    SE_RZF = zeros(K,L);
end

if nargout > 2
    SE_MMMSE = zeros(K,L);

    %Compute sum of all estimation error correlation matrices at every BS
    C_totM = reshape(p*sum(sum(C,3),4),[M M L]);

end

if nargout > 3
    SE_ZF = zeros(K,L);
end

if nargout > 4
    SE_SMMSE = zeros(K,L);

    %Compute sum of intra-cell estimation error correlation matrices at every BS
    CR_totS = zeros(M,M,L);

    for j = 1:L
        CR_totS(:,:,j) = p*(sum(C(:,:,:,j,j),3)+sum(sum(R(:,:,:,[1:j-1 j+1:end],j),3),4));
    end
end

%% Go through all channel realizations
for n = 1:nbrOfRealizations

    %Go through all cells
    for j = 1:L

        %Extract channel estimate realizations from all UEs to BS j
        Hhatallj = reshape(Hhat(:,n,:,:,j),[M K*L]);

        %Compute MR combining in (4.11)
        V_MR = Hhatallj(:,K*(j-1)+1:K*j);

        if nargout > 1 %Compute RZF combining in (4.9)
            V_RZF = (p*V_MR)/(p*(V_MR'*V_MR)+eyeK);
        end

        if nargout > 2 %Compute M-MMSE combining in (4.7)
            V_MMMSE = (p*(Hhatallj*Hhatallj')+C_totM(:,:,j)+eyeM)\(p*V_MR);
        end

        if nargout > 3 %Compute ZF combining in (4.10), with the small regularization term 1e-12 for numerical stability
            V_ZF = V_MR/(V_MR'*V_MR+1e-12*eyeK);
        end

        if nargout > 4 %Compute S-MMSE combining in (4.8)
            V_SMMSE = (p*(V_MR*V_MR')+CR_totS(:,:,j)+eyeM)\(p*V_MR);
        end

        %Go through all UEs in cell j
        for k = 1:K

            %%MR combining
            v = V_MR(:,k); %Extract combining vector

            %Compute numerator and denominator of instantaneous SINR in (4.3)
            numerator = p*abs(v'*Hhat(:,n,k,j,j))^2;
            denominator = p*sum(abs(v'*Hhatallj).^2) + v'*(C_totM(:,:,j)+eyeM)*v - numerator;

            %Compute instantaneous SE for one channel realization
            SE_MR(k,j) = SE_MR(k,j) + prelogFactor*real(log2(1+numerator/denominator));

            %%ZF combining
            if nargout > 3

                v = V_ZF(:,k); %Extract combining vector

                %Compute numerator and denominator of instantaneous SINR in (4.3)
                numerator = p*abs(v'*Hhat(:,n,k,j,j))^2;
                denominator = p*sum(abs(v'*Hhatallj).^2) + v'*(C_totM(:,:,j)+eyeM)*v - numerator;

                %Compute instantaneous SE for one channel realization
                SE_ZF(k,j) = SE_ZF(k,j) + prelogFactor*real(log2(1+numerator/denominator));

            end

            %%RZF combining
            if nargout > 1

                v = V_RZF(:,k); %Extract combining vector

                %Compute numerator and denominator of instantaneous SINR in (4.3)
                numerator = p*abs(v'*Hhat(:,n,k,j,j))^2;
                denominator = p*sum(abs(v'*Hhatallj).^2) + v'*(C_totM(:,:,j)+eyeM)*v - numerator;

                %Compute instantaneous SE for one channel realization
                SE_RZF(k,j) = SE_RZF(k,j) + prelogFactor*real(log2(1+numerator/denominator));

            end

            %%S-MMSE combining
            if nargout > 4

                v = V_SMMSE(:,k); %Extract combining vector

                %Compute numerator and denominator of instantaneous SINR in (4.3)
                numerator = p*abs(v'*Hhat(:,n,k,j,j))^2;
                denominator = p*sum(abs(v'*Hhatallj).^2) + v'*(C_totM(:,:,j)+eyeM)*v - numerator;

                %Compute instantaneous SE for one channel realization
                SE_SMMSE(k,j) = SE_SMMSE(k,j) + prelogFactor*real(log2(1+numerator/denominator));

            end

            %%M-MMSE combining
            if nargout > 2

                v = V_MMMSE(:,k); %Extract combining vector

                %Compute numerator and denominator of instantaneous SINR in (4.3)
                numerator = p*abs(v'*Hhat(:,n,k,j,j))^2;
                denominator = p*sum(abs(v'*Hhatallj).^2) + v'*(C_totM(:,:,j)+eyeM)*v - numerator;

                %Compute instantaneous SE for one channel realization
                SE_MMMSE(k,j) = SE_MMMSE(k,j) + prelogFactor*real(log2(1+numerator/denominator));

            end

        end

    end

end

author: emilbjornson

repo: massivemimobook


text:
function [P,E,logw,acc] = spm_mci_ais_single_vl (mcmc,M,U,Y,vl)
% Produce a single independent sample from posterior using AIS
% FORMAT [P,E,logw,acc] = spm_mci_ais_single_vl (mcmc,M,U,Y,vl)
%
% mcmc      Sampling settings
% M         Model structure
% U         Input structure
% Y         Data
% vl        Variational Laplace solution
%               .Ep                 Posterior Mean
%               .Cp                 Posterior Covariance
%
% P         [Np x 1] sample
% E         Negative log joint
% logw      Contribution to model evidence
% acc       acc(j) is acceptance rate at temperature j
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_mci_ais_single_vl.m 6548 2015-09-11 12:39:47Z will $

beta=mcmc.beta;
nprop=mcmc.nprop;
prop=mcmc.prop;
scale=mcmc.scale;

J=length(beta);

% Sample from VL posterior (in subspace)
x(:,1) = spm_normrnd (vl.mr,vl.Cr,1);
[L(1),L2] = spm_mci_joint (x(:,1),M,U,Y);

% Sample in original space
xorg = M.V*x(:,1)+M.vpE;
e = xorg - vl.Ep;
log_pvl = vl.const - 0.5*e'*vl.Lambdap*e;
Lsum = (beta(1)-beta(2))*log_pvl+(beta(2)-beta(1))*L(1);

acc(1) = 1;
for j=2:J,
    xs=x(:,j-1); Ls=L(:,j-1);acc(j)=0;

    % Generate sample at next temperature using
    % Langevin Monte Carlo
    lgv_mcmc.init=xs;
    lgv_mcmc.maxits=nprop;

    [tmp,stats] = spm_mci_lgv_vl (lgv_mcmc,M,U,Y,vl,beta(j));

    xs = stats.P(:,end);
    Ls = -stats.E(:,end);
    L2 = stats.L2(end);
    acc (j) = any(stats.acc(2:end));

    x(:,j)=xs;
    L(j)=Ls;
    if j < J
        xorg = M.V*x(:,j)+M.vpE;
        e = xorg - vl.Ep;
        log_pvl = vl.const - 0.5*e'*vl.Lambdap*e;
        Lsum = Lsum + (beta(j)-beta(j+1))*log_pvl+(beta(j+1)-beta(j))*L(j);
    end
end
logw=Lsum;
P=x(:,J);
E=-L(J);

if mcmc.rec_traj
    nj=size(x,2);
    traj=spm_vec(M.pE)*ones(1,nj)+V*x;
else
    traj=[];
end

author: spm


text:
function [truth, normal_axis] = is2dlike(N)

chkarg(istypesizeof(N, 'int', [1, Axis.count]), ...
	'"N" should be length-%d row vector with integer elements.');

[~, imin] = min(N);
n = Axis.elems(imin);
[h, v] = cycle(n);

truth =  N(n)/N(h) < 0.1 && N(n)/N(v) < 0.1;
normal_axis = n;

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function X = spinv(A,tol,method)
%PINV  sparse Pseudoinverse.
%   Copyright 1984-2013 The MathWorks, Inc.
if nargin < 3 || strcmpi(method,'svd')
    [U,S,V] = svds(A,size(A,1));
    s = diag(S);
    if nargin < 2
        tol = max(size(A)) * eps(norm(s,inf));
    end
    r1 = sum(s > tol)+1;
    V(:,r1:end) = [];
    U(:,r1:end) = [];
    s(r1:end) = [];
    s = 1./s(:);
    X = bsxfun(@times,V,s.')*U';
elseif strcmpi(method,'qr')
    [H,R,E] = qr(A);
    dr = diag(R);
    if isempty(tol)
    tol = 1e4*eps(max(abs(dr)));
    end
    k = abs(dr)>tol;
    X = E*[inv(R(k,k))*H(:,k)'; sparse(sum(~k),size(H,2))];
end

author: goGPS-Project


text:
classdef DISOpticalFlow < handle
    %DISOPTICALFLOW  DIS optical flow algorithm
    %
    % This class implements the Dense Inverse Search (DIS) optical flow
    % algorithm. More details about the algorithm can be found at
    % [Kroeger2016]. Includes three presets with preselected parameters to
    % provide reasonable trade-off between speed and quality. However, even
    % the slowest preset is still relatively fast, use cv.calcOpticalFlowDF
    % if you need better quality and don't care about speed.
    %
    % This implementation includes several additional features compared to the
    % algorithm described in the paper, including spatial propagation of flow
    % vectors (cv.DISOpticalFlow.UseSpatialPropagation), as well as an option
    % to utilize an initial flow approximation passed to cv.DISOpticalFlow.calc
    % (which is, essentially, temporal propagation, if the previous frame's
    % flow field is passed).
    %
    % ## References
    % [Kroeger2016]:
    % > Till Kroeger, Radu Timofte, Dengxin Dai, and Luc Van Gool.
    % > "Fast optical flow using dense inverse search". In Proceedings of the
    % > European Conference on Computer Vision (ECCV), 2016.
    %
    % See also: cv.DISOpticalFlow.calc
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % Finest level of the Gaussian pyramid on which the flow is computed
        % (zero level corresponds to the original image resolution). The final
        % flow is obtained by bilinear upscaling.
        FinestScale
        % Size of an image patch for matching (in pixels). Normally, default
        % 8x8 patches work well enough in most cases.
        PatchSize
        % Stride between neighbor patches. Must be less than patch size. Lower
        % values correspond to higher flow quality.
        PatchStride
        % Maximum number of gradient descent iterations in the patch inverse
        % search stage. Higher values may improve quality in some cases.
        GradientDescentIterations
        % Number of fixed point iterations of variational refinement per scale.
        % Set to zero to disable variational refinement completely. Higher
        % values will typically result in more smooth and high-quality flow.
        VariationalRefinementIterations
        % Weight of the smoothness term.
        VariationalRefinementAlpha
        % Weight of the color constancy term.
        VariationalRefinementDelta
        % Weight of the gradient constancy term.
        VariationalRefinementGamma
        % Whether to use mean-normalization of patches when computing patch
        % distance. It is turned on by default as it typically provides a
        % noticeable quality boost because of increased robustness to
        % illumination variations. Turn it off if you are certain that your
        % sequence doesn't contain any changes in illumination.
        UseMeanNormalization
        % Whether to use spatial propagation of good optical flow vectors.
        % This option is turned on by default, as it tends to work better on
        % average and can sometimes help recover from major errors introduced
        % by the coarse-to-fine scheme employed by the DIS optical flow
        % algorithm. Turning this option off can make the output flow field a
        % bit smoother, however.
        UseSpatialPropagation
    end

    %% Constructor/destructor
    methods
        function this = DISOpticalFlow(varargin)
            %DISOPTICALFLOW  Creates an instance of DISOpticalFlow
            %
            %     obj = cv.DISOpticalFlow()
            %     obj = cv.DISOpticalFlow('OptionName',optionValue, ...)
            %
            % ## Options
            % * __Preset__ preset one of:
            %   * __UltraFast__
            %   * __Fast__ (default)
            %   * __Medium__
            %
            % See also: cv.DISOpticalFlow.calc
            %
            this.id = DISOpticalFlow_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.DISOpticalFlow
            %
            if isempty(this.id), return; end
            DISOpticalFlow_(this.id, 'delete');
        end
    end

    %% DenseOpticalFlow
    methods
        function flow = calc(this, I0, I1, varargin)
            %CALC  Calculates an optical flow
            %
            %     flow = obj.calc(I0, I1)
            %     flow = obj.calc(I0, I1, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __I0__ first 8-bit single-channel input image.
            % * __I1__ second input image of the same size and the same type
            %   as `I0`.
            %
            % ## Output
            % * __flow__ computed flow image that has the same size as `I0`
            %   and type `single` (2-channels).
            %
            % ## Options
            % * __InitialFlow__ specify the initial flow. Not set by default.
            %
            % See also: cv.DISOpticalFlow
            %
            flow = DISOpticalFlow_(this.id, 'calc', I0, I1, varargin{:});
        end

        function collectGarbage(this)
            %COLLECTGARBAGE  Releases all inner buffers
            %
            %     obj.collectGarbage()
            %
            DISOpticalFlow_(this.id, 'collectGarbage');
        end
    end

    %% Algorithm
    methods (Hidden)
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.DISOpticalFlow.empty
            %
            DISOpticalFlow_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Returns true if the algorithm is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the algorithm is empty (e.g. in the very
            %   beginning or after unsuccessful read).
            %
            % See also: cv.DISOpticalFlow.clear
            %
            b = DISOpticalFlow_(this.id, 'empty');
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.DISOpticalFlow.save, cv.DISOpticalFlow.load
            %
            name = DISOpticalFlow_(this.id, 'getDefaultName');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in a file storage.
            %
            % See also: cv.DISOpticalFlow.load
            %
            DISOpticalFlow_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from a file storage.
            % The previous model state is discarded.
            %
            % See also: cv.DISOpticalFlow.save
            %
            DISOpticalFlow_(this.id, 'load', fname_or_str, varargin{:});
        end
    end

    %% Getters/Setters
    methods
        function value = get.FinestScale(this)
            value = DISOpticalFlow_(this.id, 'get', 'FinestScale');
        end
        function set.FinestScale(this, value)
            DISOpticalFlow_(this.id, 'set', 'FinestScale', value);
        end

        function value = get.PatchSize(this)
            value = DISOpticalFlow_(this.id, 'get', 'PatchSize');
        end
        function set.PatchSize(this, value)
            DISOpticalFlow_(this.id, 'set', 'PatchSize', value);
        end

        function value = get.PatchStride(this)
            value = DISOpticalFlow_(this.id, 'get', 'PatchStride');
        end
        function set.PatchStride(this, value)
            DISOpticalFlow_(this.id, 'set', 'PatchStride', value);
        end

        function value = get.GradientDescentIterations(this)
            value = DISOpticalFlow_(this.id, 'get', 'GradientDescentIterations');
        end
        function set.GradientDescentIterations(this, value)
            DISOpticalFlow_(this.id, 'set', 'GradientDescentIterations', value);
        end

        function value = get.VariationalRefinementIterations(this)
            value = DISOpticalFlow_(this.id, 'get', 'VariationalRefinementIterations');
        end
        function set.VariationalRefinementIterations(this, value)
            DISOpticalFlow_(this.id, 'set', 'VariationalRefinementIterations', value);
        end

        function value = get.VariationalRefinementAlpha(this)
            value = DISOpticalFlow_(this.id, 'get', 'VariationalRefinementAlpha');
        end
        function set.VariationalRefinementAlpha(this, value)
            DISOpticalFlow_(this.id, 'set', 'VariationalRefinementAlpha', value);
        end

        function value = get.VariationalRefinementDelta(this)
            value = DISOpticalFlow_(this.id, 'get', 'VariationalRefinementDelta');
        end
        function set.VariationalRefinementDelta(this, value)
            DISOpticalFlow_(this.id, 'set', 'VariationalRefinementDelta', value);
        end

        function value = get.VariationalRefinementGamma(this)
            value = DISOpticalFlow_(this.id, 'get', 'VariationalRefinementGamma');
        end
        function set.VariationalRefinementGamma(this, value)
            DISOpticalFlow_(this.id, 'set', 'VariationalRefinementGamma', value);
        end

        function value = get.UseMeanNormalization(this)
            value = DISOpticalFlow_(this.id, 'get', 'UseMeanNormalization');
        end
        function set.UseMeanNormalization(this, value)
            DISOpticalFlow_(this.id, 'set', 'UseMeanNormalization', value);
        end

        function value = get.UseSpatialPropagation(this)
            value = DISOpticalFlow_(this.id, 'get', 'UseSpatialPropagation');
        end
        function set.UseSpatialPropagation(this, value)
            DISOpticalFlow_(this.id, 'set', 'UseSpatialPropagation', value);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function X = recon(D,P)

X = D.*cos(P) + i*D.*sin(P);

author: posenhuang

repo: deeplearningsourceseparation


text:
function H = psd_project_rows(H,tol)
  % H = psd_project_rows(H)
  if nargin < 2
    tol = 0;
  end
  m = round(sqrt(size(H,2)));
  assert(m*m==size(H,2));
  for i = 1:size(H,1)
    A = reshape(H(i,:),m,m);
    [V,D] = eig(triu(A)+triu(A,1)','vector');
    B = V*(max(D,tol).*V');
    H(i,:) = B(:);
    if any(~isreal(B(:)))
      keyboard
    end
  end

end

author: alecjacobson

repo: gptoolbox