text:
function Callbacks_p_Tube_VT_GUI25(f,C)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=C{1,1};
y=C{1,2};
a=C{1,3};
b=C{1,4};
u=C{1,5};
v=C{1,6};
m=C{1,7};
n=C{1,8};
lengthbutton=C{1,9};
widthbutton=C{1,10};
enterType=C{1,11};
enterString=C{1,12};
enterLabel=C{1,13};
noPanels=C{1,14};
noGraphicPanels=C{1,15};
noButtons=C{1,16};
labelDist=C{1,17};%distance that the label is below the button
noTitles=C{1,18};
buttonTextSize=C{1,19};
labelTextSize=C{1,20};
textboxFont=C{1,21};
textboxString=C{1,22};
textboxWeight=C{1,23};
textboxAngle=C{1,24};
labelHeight=C{1,25};
fileName=C{1,26};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%PANELS
for j=0:noPanels-1
uipanel('Parent',f,...
'Units','Normalized',...
'Position',[x(1+4*j) y(1+4*j) x(2+4*j)-x(1+4*j) y(3+4*j)-y(2+4*j)]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GRAPHIC PANELS
for i=0:noGraphicPanels-1
switch (i+1)
case 1
graphicPanel1 = axes('parent',f,...
'Units','Normalized',...
'Position',[a(1+4*i) b(1+4*i) a(2+4*i)-a(1+4*i) b(3+4*i)-b(2+4*i)],...
'GridLineStyle','--');
case 2
graphicPanel2 = axes('parent',f,...
'Units','Normalized',...
'Position',[a(1+4*i) b(1+4*i) a(2+4*i)-a(1+4*i) b(3+4*i)-b(2+4*i)],...
'GridLineStyle','--');
case 3
graphicPanel3 = axes('parent',f,...
'Units','Normalized',...
'Position',[a(1+4*i) b(1+4*i) a(2+4*i)-a(1+4*i) b(3+4*i)-b(2+4*i)],...
'GridLineStyle','--');
case 4
graphicPanel4 = axes('parent',f,...
'Units','Normalized',...
'Position',[a(1+4*i) b(1+4*i) a(2+4*i)-a(1+4*i) b(3+4*i)-b(2+4*i)],...
'GridLineStyle','--');
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%TITLE BOXES
for k=0:noTitles-1
switch (k+1)
case 1
titleBox1 = uicontrol('parent',f,...
'Units','Normalized',...
'Position',[u(1+4*k) v(1+4*k) u(2+4*k)-u(1+4*k) v(3+4*k)-v(2+4*k)],...
'Style','text',...
'FontSize',textboxFont{k+1},...
'String',textboxString(k+1),...
'FontWeight',textboxWeight{k+1},...
'FontAngle',textboxAngle{k+1});
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%BUTTONS
for i=0:(noButtons-1)
enterColor='w';
if strcmp(enterType{i+1},'pushbutton')==1 ||strcmp(enterType{i+1},'text')==1
enterColor='default';
end
if (strcmp(enterLabel{1,(i+1)},'')==0 &&...
        strcmp(enterLabel{1,(i+1)},'...')==0) %i.e. there is a label
%creating a label for some buttons
uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i)-labelDist-labelHeight(i+1) ...
(m(2+2*i)-m(1+2*i)) labelHeight(i+1)],...
'Style','text',...
'String',enterLabel{i+1},...
'FontSize', labelTextSize(i+1),...
'HorizontalAlignment','center');
end
switch (i+1)
case 1
button1=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button1Callback);
case 2
button2=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button2Callback);
case 3
button3=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button3Callback);
case 4
button4=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button4Callback);
case 5
button5=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button5Callback);
case 6
button6=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button6Callback);
case 7
button7=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button7Callback);
case 8
button8=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button8Callback);
case 9
button9=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button9Callback);
case 10
button10=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button10Callback);
case 11
button11=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button11Callback);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%USER CODE FOR THE VARIABLES AND CALLBACKS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize Variables
rG=0.7; rL=0.7;
fs=10000; ls=1.75;
fno=1;
% area=ones(1,10);
c=35000;
nfft=4000;
p=10;
ipd=100;
fileout='out_p_tube';
sliderpos1=0.7;
sliderpos2=0.7;
filename='area_1';
source='figure5.40';
seq_length=4001;

% Name the GUI
    set(f,'Name','p_Tube_VT_GUI');

% CALLBACKS
% Callback for button1 -- mat file with areas and other parameters

set_slider_pos()

 function button1Callback(h,eventdata)
     fno=get(button1,'val');
     set(button6,'val',(rG+1)/2);
     set(button7,'val',(rL+1)/2);
 end

% Callback for button2 -- read in sampling rate, fs, of vocal tract
% synthesis
 function button2Callback(h,eventdata)
     fs=str2num(get(button2,'string'));
 end

% Callback for button3 -- report length of each uniform tube, ls
 function button3Callback(h,eventdata)
     ls=str2num(get(button3,'string'));
 end

% Callback for button4 -- glottal reflection coefficient, rG, edit button
 function button4Callback(h,eventdata)
     set(button4,'string',num2str(sliderpos1));
     rG=sliderpos1;
 end

% Callback for button5 -- lips reflection coefficient, rL, edit button
 function button5Callback(h,eventdata)
     set(button5,'string',num2str(sliderpos2));
     rL=sliderpos2;
 end

% Callback for button6 -- rG slider (range of -1 to +1) edit button
 function button6Callback(h,eventdata)
     sliderpos1=get(button6,'val')*2-1;
     button4Callback(h,eventdata);
 end

% Callback for button7 -- rL slider  (range of -1 to +1) edit button
 function button7Callback(h,eventdata)
     sliderpos2=get(button7,'val')*2-1;
     button5Callback(h,eventdata);
 end

% Callback for button8 -- name of output text file, fileout, for saving
% p-tube vocal tract resonances
 function button8Callback(h,eventdata)
     fileout=get(button8,'string');
 end

% Callback for button11 -- ipd:  pitch period in samples
    function button11Callback(h,eventdata)
        ipd=num2str(get(button11,'string'));
    end

% Callback for button9 -- run the p-tube VT model code
 function button9Callback(h,eventdata)

% send message to change either rG or rL to be less than 1 for realistic
% vowel sound
    if (rG == 1 & rL == 1)
    uiwait(msgbox('change rG or rL to be less than 1 for correct vowel sound',...
        'set rG, rL','modal'));
    end

% read in source file for areas and other parameters
     clear area;
     filename=['area_',num2str(fno),'.mat'];
     areas=load(filename);

% extract p-tube vocal tract information from areas structure, including the
% p-tube areas, areas.area. the source of the area data, areas.source,
% the number of tubes, areas.p, the size of fft for synthesis,
% areas.nfft, the sampling rate of the simulation, areas.fs, the
% speed of sound in cm/sec, areas.c, and the length of each of the
% p-uniform tubes, areas.ls
     area=areas.area;
     source=areas.source;
     p=areas.p;
     nfft=areas.nfft;
     fs=areas.fs;
     c=areas.c;
     ls=areas.ls;

% update buttons for fs and ls by writing the values obtained above into
% the appropriate edit buttons (button2 and button3)
    set(button2,'string',num2str(fs));
    set(button3,'string',num2str(ls));

% open text file for saving vocal tract resonance information
    fid=fopen(fileout,'wt');

% create input for plotting vocal tract area function across p sections
    yc=[area(1) area(1) area(1)];
    xc=[0 1 1 1 1];
    for index=2:p-1
        yc=[yc area(index) area(index) area(index) area(index)];
        xc=[xc index index index index];
    end
    yc=[yc area(p) area(p) area(p)];
    xc=[xc p];

% plot vocal tract area function in Graphics Panel 1
    grid off;
    reset(graphicPanel1);
    axes(graphicPanel1);

% plot p-tube area function
    amax=max(yc);
    plot(xc,yc,'b','LineWidth',2);
    hold on,plot(xc,-yc,'b','LineWidth',2);
    plot([0 p],[0 0],'r--','LineWidth',2);
    axis ([0 p -amax*1.1 amax*1.1]);
    xlabel('Distance from Glottis');
    ylabel('Area(square cm)');
    stitle=sprintf(' source: %s, p: %d, ls: %6.2f',source,p,ls);

% set title in titleBox1; adjust font size to 20
    stitle1=strcat('p-Tube Vocal Tract -- ',stitle);
    set(titleBox1,'String',stitle1);
    set(titleBox1,'FontSize',20);

% transform from vocal tract areas areas(1:p-1) to reflection coefficients,
% k(1:p-1)
    clear k;
    k(1:p-1)=(-area(2:p)+area(1:p-1))./(area(2:p)+area(1:p-1));
    r(1:p-1)=-k(1:p-1);
    k=[rG -r rL];

% plot set of reflection coefficients on Graphics Panel 2
    reset(graphicPanel2);
    axes(graphicPanel2);
    stem(0:p,k,'r','LineWidth',2);
    axis([0 p -1 1]);
    xlabel('Distance from Glottis'); ylabel('Reflection Coefficients');

% print out (to ascii text file) the set of reflection coefficients and areas
    fprintf(fid,'areas-to-spectrum using p: %d, nfft: %d,',p,nfft);
    fprintf(fid,' fs: %d \n',fs);
    for index=1:p-1
        fprintf(fid,'index: %d, k:%6.4f,',index,k(index));
        fprintf(fid,' areas:%6.3f, %6.3f \n',area(index+1),area(index));
    end

%
% solve nodal equations for the p tubes
%
% matrix1 columns are uG, u1+, u2+, ..., up+, u(p+1)+
% matrix1 rows are n=0,1,2,...,NL
% matrix2 columns are u1-, u2-, ...,up-
% matrix2 rows are n=1,2,...,NL+1
%
    NL=3000; % number of time slots for solution

% initialize nodal equations
    uplus=zeros(p+2,NL+1);
    uminus=zeros(p,NL+1);
    uplus(1:2,1)=[1,(1+rG)/2];

% solve nodal equations for NL iterations
    for n=1:NL
        uminus(p,n+1)=-rL*uplus(p+1,n);
        for node=p-1:-1:1
            uminus(node,n+1)=-r(node)*uplus(node+1,n)+...
                (1-r(node))*uminus(node+1,n+1);
        end
        uplus(2,n+1)=uminus(1,n+1)*rG;
        for node=3:p+1
            uplus(node,n+1)=uplus(node-1,n)*(1+r(node-2))+...
                uminus(node-1,n+1)*r(node-2);
        end
        uplus(p+2,n+1)=(1+rL)*uplus(p+1,n);
    end

% save output as vocal tract impulse response
    h(1:NL+1)=uplus(p+2,1:NL+1);

% truncate impulse response to nl=500 samples
    nl=500;
    nir=0:nl-1;

% plot resulting vocal tract impulse response in Figure 3
    reset(graphicPanel3);
    axes(graphicPanel3);

    plot(nir,h(1:nl),'b','LineWidth',2);
    xpp=['Time in Samples;  fs=',num2str(fs),' samples/second'];
    xlabel(xpp),ylabel('Amplitude');
    grid on, axis tight;
    s1=sprintf('source: %s, p: %d,',source,p);
    s2=sprintf(' rG,rL: %6.3f, %6.3f, fs:%6.1f',rG,rL,fs);
    stitle=strcat(s1,s2);
    legend('Vocal Tract Impulse Response');

% transform to spectral domain and plot log magnitude response
    h=[h, zeros(1,nfft-NL-1)];

% convert to log magnitude spectrum
    hs=fft(h,nfft);
    spectrum(1:nfft/2+1)=20*log10(1.e-1+abs(hs(1:nfft/2+1)));
    vmin=min(spectrum);
    vmax=max(spectrum);

% plot log magnitude spectrum in graphics Panel 4
    reset(graphicPanel4);
    axes(graphicPanel4);

    freq=0:fs/nfft:fs/2;
    npts=nfft/2+1;
    plot(freq(1:npts),spectrum(1:npts),'g','LineWidth',2); grid;
    xlabel('Frequency(kHz)'); ylabel('Log Magnitude(dB)');
    axis([freq(1) freq(npts) vmin-2 vmax+2]);
    hold on;

% find formant locations and include on plot
    fmt=[];
    loc=[];
    nfmt=0;
    for index=10:npts-2
        if (spectrum(index+1) > spectrum(index) && ...
                spectrum(index+1) > spectrum(index+2))
            fmt=[fmt freq(index+1)];
            loc=[loc index];
            nfmt=nfmt+1;
        end
    end

% next estimate formant bandwidths
    bwlevel=6; % (db) range of bandwidth around formant peak
    if (nfmt <= 10)
        for nf=1:nfmt
            f=fmt(nf);
            index=loc(nf);
            peak=spectrum(index);
            [peakl,peakh]=findpeak(spectrum,index,peak,bwlevel);
            fr(nf)=f;
            bw(nf)=freq(peakh)-freq(peakl);
        end
    else
        peakspectrum=max(spectrum);
        nfe=1;
        for nf=1:nfmt
            f=fmt(nf);
            index=loc(nf);
            peak=spectrum(index);
            if (peak >= peakspectrum-20)
                fr(nfe)=f;
                bw(nfe)=0;
                fmt(nfe)=f;
                nfe=nfe+1;
            end
        end
        nfmt=nfe-1;
    end

% print out results to output text file
    fprintf(fid,'p-tube model results: \n');
    fprintf(fid,'rL,rG: %4.1f %4.1f \n',rL,rG);
    for formant=1:nfmt
        fprintf(fid,'formant frequency: %6.1f,',fr(formant));
        fprintf(fid,' formant bandwidth: %6.1f \n',bw(formant));
    end
    fprintf(fid,'\n \n');

% print formant locations on plots
    for nf=1:nfmt
        frfmt=[fmt(nf) fmt(nf)];
        vrfmt=[vmin vmax];
        plot(frfmt,vrfmt,'r','LineWidth',2);
        string=strcat('fmt',int2str(nf),':',int2str(fmt(nf)));
        if (mod(nf,2) == 0)
            text(fmt(nf)-200,vmax*.9,string);
        else
            text(fmt(nf)-200,vmin*1.1,string);
        end
    end

% convolve p-tube impulse response with periodic pulse train of period ipd
% samples (default is 100 Hz fundamental at fs=10000 Hz), and play out resulting
% speech waveform -- note wierd sound for rG=rL=1 since impulse response
% never dies out with no loss
        ipd=str2num(get(button11,'string'));
                   if (ipd < 0)
            waitfor(errordlg('IPD must be positive'));
        end
        seq_length=40*ipd;  % sequence length in samples
        xin=zeros(1,seq_length);
        xin(1:ipd:seq_length)=1;
        yout=conv(xin,h);
        soundsc(yout,fs);

% close up printing file
    fclose('all');
 end

% Callback for button10 -- close the GUI
 function button10Callback(h,eventdata)
     close(gcf);
 end

function set_slider_pos()
    set(button6,'val',rG+0.15); %when set at exactly rG, slider bar is at 0.4
    set(button7,'val',rL+0.15); %when set at exactly rL, slider bar is at 0.4
end
end


text:
function mrUtilMakeColorbar(cmap, labels, desc, fileName, figNum, axisHandle)
% Simple utility to make and save a stand-alone colorbar.
% mrUtilMakeColorbar(cmap, labels, desc, [fileName = no file saved], [figNum=figure], [axisHandle=gca])
%
% labels can be a simple cell-array of stings. You'll get one tick mark per
% label and they will be evenly distributed along the colorbar.
%
% If you specify a file name, the figure will be 'print'ed to that file.
% Valid file name extensions are '.tif' (300 dpi), '.png' (120 dpi), and
% '.eps'. If you don't provide an extension (eg. just specify
% '/path/to/some/file'), then both a png (good for powerpoint) and an eps
% will be saved.
%
% EXAMPLE:
% mrUtilMakeColorbar(autumn(256),{'0','0.5','1.0'},'Data','DataLegend')
%
% To change the font and line color to white:
% mrUtilMakeColorbar(autumn(256),{'0','0.5','1.0'},'Data')
% set(gcf,'InvertHardCopy','off','color','k');
% set(gca,'ycolor','w','xcolor','w','color', 'k');
% set(get(gca,'XLabel'),'color','w');
% mrUtilPrintFigure(gcf, 'DataLegend');
%
% HISTORY:
% 2005.08.20 RFD: wrote it.

fontSize = 12;
fontName = 'Helvetica';

if(~exist('desc','var'))
    desc = '';
end
if(~exist('fileName','var'))
  fileName = '';
  figName = desc;
else
  [p,f,e] = fileparts(fileName);
  figName = f;
end
if(~exist('figNum','var')||isempty(figNum))
    figNum = figure;
    set(figNum,'Name',figName);
    figExists = false;
else
    figExists = true;
end
if(~exist('axisHandle','var')||isempty(axisHandle))
    figure(figNum);
    axisHandle = gca;
else
  axes(axisHandle);
end

if(size(labels,1)>size(labels,2))
    vert = true;
else
    vert = false;
end

nColors = size(cmap,1);
nLabels = length(labels);

figure(figNum); axes(axisHandle); cla(axisHandle);
height = ceil(nColors.*0.04);
if(vert)
    image(repmat([1:nColors]', 1, height));
    axis xy;
    set(axisHandle,'ytick',linspace(1,nColors,nLabels),'xtick',[],'FontSize',fontSize,'FontName',fontName);
    set(axisHandle,'yticklabel',labels,'xticklabel',[],'YAxisLocation','right');
    sz = [70 300];
    apos = [-0.2 0.055 0.945 0.88];
else
    image(repmat([1:nColors], height, 1));
    set(axisHandle,'xtick',linspace(1,nColors,nLabels),'ytick',[],'FontSize',fontSize,'FontName',fontName);
    set(axisHandle,'xticklabel',labels,'yticklabel',[]);
    sz = [300 70];
    apos = [0.055 0.35 0.88 0.945];
end
axis equal tight on;

if(~figExists)
   pos = get(gcf,'Position');
   set(figNum,'Position',[pos(1) pos(2) sz], 'PaperPositionMode', 'auto');
   set(axisHandle,'Position',apos);
   colormap(cmap);
end
if(vert)
    ylabel(desc,'FontSize',fontSize,'FontName',fontName);
else
    xlabel(desc,'FontSize',fontSize,'FontName',fontName);
end
if(~isempty(fileName))
    [p,f,e] = fileparts(fileName);
    if(isempty(e))
      print(figNum, '-depsc', '-tiff', '-cmyk', [fileName '.eps']);
      print(figNum, '-dpng', '-r120', [fileName '.png']);
    elseif(strcmp(e,'.eps'))
        print(figNum, '-depsc', '-tiff', '-cmyk', fileName);
    elseif(strcmp(e,'.tif'))
        print(figNum, '-dtiff', '-r300', fileName);
    else
        print(figNum, '-dpng', '-r120', fileName);
    end
end
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function Q = mg_ns_smooth(A,sweeps,smooth,stype)
%mg_ns_smooth   smoothers for GMG iteration (Navier-Stokes)
%   Q = mg_ns_smooth(A,sweeps,smooth,stype)
%   input
%          A         coefficient matrix
%          sweeps    number of directions for multidirectional Gauss-Seidel
%          smooth    smoother (1, Jacobi, 2, Gauss-Seidel, 3 ILU)
%          stype     type of smoother (1 point, 2 line)
%   output
%          Q         structure, smoothing operators in factored form
%
%   IFISS function: HCE; 18 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage

N=size(A,1);
n = sqrt(size(A,1));
if stype==2
   % line Gauss-Seidel
   Q1 = tril(A,1);
   [L1,U1] = lu(Q1);
   if sweeps>=2
      Q2 = diag(diag(A,0)) + diag(diag(A,-n),-n) + diag(diag(A,n),n) ...
         + diag(diag(A,-n-1),-n-1) + diag( diag(A,-1),-1)        ...
         + diag(diag(A,n-1),n-1);
      [L2,U2] = lu(Q2);
      if sweeps>=3
         Q3 = triu(A,-1);
         [L3,U3] = lu(Q3);
      if sweeps==4
         Q4 = diag(diag(A,0)) + diag(diag(A,-n),-n) + diag(diag(A,n),n) ...
            + diag(diag(A,n+1),n+1) + diag( diag(A,1),1)            ...
            + diag(diag(A,-n+1),-n+1);
         [L4,U4] = lu(Q4);
         else
            L4=sparse(N,N); U4=L4;
         end
      else
         L3=sparse(N,N); U3=L3; L4=L3; U4=L3;
      end
   else
      L2=sparse(N,N); U2=L2; L3=L2; U3=L2; L4=L2; U4=L2;
   end
else
   % point smoothers
   if smooth==3,
      % ILU
      [L1,U1]=ilu0(A);
      L2=sparse(size(L1)); U2=L2; L3=L2; U3=L2; L4=L2; U4=L2;
   elseif smooth==2,
      % point Gauss-Seidel
      Q1 = tril(A,0); [L1,U1] = lu(Q1);
      L2=sparse(size(L1)); U2=L2; L3=L2; U3=L2; L4=L2; U4=L2;
   else
      % point damped Jacobi
      omega=8/9; % relaxation factor for damped Jacobi
      Q1 = (1/omega)*spdiags(diag(A),0,n*n,n*n);[L1,U1] = lu(Q1);
      L2=sparse(size(L1)); U2=L2; L3=L2; U3=L2; L4=L2; U4=L2;
   end
end
Q = struct('L1',L1,'L2',L2,'L3',L3,'L4',L4, ...
           'U1',U1,'U2',U2,'U3',U3,'U4',U4);

author: johannesgerer

repo: jburkardt-m


text:
function X = tenrand(varargin)
%TENRAND Uniformly distributed pseudo-random tensor.
%
%   X = TENRAND(SZ) forms a tensor of size SZ with pseudo-random
%   values drawn from a uniform distribution on the unit interval.
%
%   TENRAND(SZ) is equivalent to TENSOR(RAND(SZ(1),SZ(2),...),SZ).
%
%   See also TENSOR, SPTENRAND, RAND.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

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

if nargin == 1
    sz = varargin{1};
else
    sz = cell2mat(varargin);
end

data = rand([sz 1 1]);
X = tensor(data,sz);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function x = exponential_01_cdf_inv ( cdf )

%*****************************************************************************80
%
%% EXPONENTIAL_01_CDF_INV inverts the Exponential 01 CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%    0.0 <= CDF <= 1.0.
%
%    Output, real X, the corresponding argument.
%
  if ( cdf < 0.0 | 1.0 < cdf )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'EXPONENTIAL_01_CDF_INV - Fatal error!\n' );
    fprintf ( 1, '  CDF < 0 or 1 < CDF.\n' );
    error ( 'EXPONENTIAL_01_CDF_INV - Fatal error!' );
  end

  x = - log ( 1.0 - cdf );

  return
end


text:
% convert the trainable parameters in the graph into a vector such that we
% can call standard optimization packages to optimize the network
% parameters in batch mode.
%
function [cost, grad] = DNN_cost_wrapper(W, layer, data, para, mode)

% retrieve the weights from W and assign it to the correct layers
layer = NetWeights_vec2layer(W, layer, 0);

% compute the forward and backward passes and get the gradients
[cost_func, layer] = DNN_Cost10(layer, data, para, mode);
cost = cost_func.cost;

% retrieve the gradient from layer and store it into vector format
[grad] = NetWeights_layer2vec(layer, 1, para.useGPU);

cost = gather(cost);
grad = gather(grad);

end

author: singaxiong

repo: SignalGraph


text:
function [map, dist] = imPointsInfluenceZones(varargin)
%IMPOINTSINFLUENCEZONES Maps influence zones of a set of 2D/3D points
%
%   MAP = imPointsInfluenceZones(DIM, POINTS)
%   DIM is a 1-by-2 or 1-by-3 row vector containing dimensions of ouput
%   label map. POINTS is a N-by-2 or N-by-3 array of coordinates.
%
%   MAP = imPointsInfluenceZones(LX, LY, POINTS)
%   MAP = imPointsInfluenceZones(LX, LY, LZ, POINTS)
%   LX and LY and optionnaly LZ are row vectors containing the values of
%   XData, YData, and ZData. The size of MAP if given by:
%   * length(LY)-by-length(LX) in case of 2D points
%   * length(LY)-by-length(LX)-by-length(LZ) in case of 3D points
%
%   [MAP, DIST] = imPointsInfluenceZones(...)
%   Also returns the distance map, containing for each pixel or voxel the
%   distance to the closest point in the input point set.
%
%   Example
%     % Planar example
%     points = rand(20, 2) * 200;
%     map = imPointsInfluenceZones([200 200], points);
%     rgb = label2rgb(map, jet(20), 'w', 'shuffle');
%     figure; imshow(rgb);
%
%     % 3D example
%     dim = [100 100 100];
%     np = 200;
%     points = rand(np, 3) * dim(1);
%     tic; [map, dist] = imPointsInfluenceZones(dim, points); toc
%     rgb = label2rgb3d(map, jet(np+1), 'w', 'shuffle');
%     figure; imshow(rgb(:,:,:,50));
%
%   See also
%     imDistance, imvoronoi2d, imvoronoi3d

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2019-05-10,    using Matlab 9.6.0.1072779 (R2019a)
% Copyright 2019 INRA - Cepia Software Platform.

%% Extract input arguments

% checkup on input argument number
if nargin < 2
    error('Requires at least two input arguments');
end

% extraction of image dimensions
var1 = varargin{1};
if size(var1, 1) == 1 && (size(var1, 2) == 2 || size(var1, 2) == 3)
    % first argument contains the size of the output image
    lx = 1:var1(2);
    ly = 1:var1(1);
    if length(var1) == 3
        lz = 1:var1(3);
    end
    varargin(1) = [];

elseif size(var1, 1) == 1 && size(varargin{2}, 1) == 1
    % first and second arguments contain vector for each coordinate
    % respectively
    lx = var1;
    ly = varargin{2};
    if size(varargin{3}, 1) == 1
        lz = varargin{3};
        varargin(1:3) = [];
    else
        varargin(1:2) = [];
    end

else
    error(['wrong input arguments in ' mfilename]);
end

% extraction of points
points = varargin{1};

%% Initialisations

% number of points
nPoints = size(points, 1);
nd = size(points, 2);

% size of output image
Nx = length(lx);
Ny = length(ly);

%% Generation of distance function

if nd == 2
    % allocate memory for label map
    map = zeros([Ny Nx]);
    % initialize distance map with arbitrarily large distance
    dist = inf * ones([Ny Nx]);

    % pixels coordinates
    [x, y] = meshgrid(lx, ly);

    % update distance for each point
    for i = 1:nPoints
        % squared distance from each pixel to current point
        di = (x - points(i,1)).^2 + (y - points(i,2)).^2;

        % update arrays for current influence zone
        inds = di < dist;
        map(inds) = i;
        dist(inds) = di(inds);
    end

    % convert squared distance to Euclidean distance
    dist = sqrt(dist);

elseif nd == 3
    % size in third dimension
    Nz = length(lz);

    % allocate memory for label map
    map = zeros([Ny Nx Nz]);

    % initialize distance map with arbitrarily large distance
    dist = inf * ones([Ny Nx Nz]);

    % pixels coordinates
    [x, y, z] = meshgrid(lx, ly, lz);

    % update distance for each point
    for i = 1:nPoints
        % distance from each pixel to current point
        di = (x - points(i,1)).^2 + (y - points(i,2)).^2 + (z - points(i,3)).^2;

        % update arrays for current influence zone
        inds = di < dist;
        map(inds) = i;
        dist(inds) = di(inds);
    end

    % convert squared distance to Euclidean distance
    dist = sqrt(dist);
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [Pro,Res] = interpolationAMGa(A,isC)
%% INTERPOLATIONAMGA construct prolongation and restriction matrices
%
% [Pro,Res] = INTERPOLATIONAMGA(A,isC) construct prolongation and
% restriction matrices use matrix-dependent interpolation. Each fine nodes
% use only one coarse node.
%
% In the input, A is a SPD matrix and isC is a logical array to indicate
% nodes in coarse matrix. In the output Pro and Res are prolongation and
% restriction matrices satisfying Res = Pro'.
%
% The submatrix A_{cf} is used to construct the interpolation of values on
% fine nodes from that of coarse nodes. The weight is normalized to
% preserve the constant.
%
% Example
%   load lakemesh
%   A = assemblematrix(node,elem);
%   [isC,As] = coarsenAMGc(A);
%   [Pro,Res] = interpolation(As,isC);
%
% See also: coarsenAMGc, amg
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

N = size(A,1);
%% Index map between coarse grid and fine grid
allNode = (1:N)';
fineNode = allNode(~isC);
% Nf = length(fineNode);
Nc = N-length(fineNode);
coarseNode = (1:Nc)';    % coarse node index
coarseNodeFineIdx = find(isC); % coarse node index in the fine grid

%% Construct prolongation and restriction operator
Afc = A(fineNode,coarseNodeFineIdx);     % matrix-dependent interpolation
[Dsum,j] = max(abs(Afc),[],2);
idx = (Dsum ~= 0);
ip = [coarseNodeFineIdx; fineNode(idx)];   % fine node index
jp = [coarseNode; j(idx)];                 % coarse node index
sp = [ones(Nc,1); ones(length(j(idx)),1)]; % weight = 1;
Pro = sparse(ip,jp,sp,N,Nc);
Res = Pro';

repo: ifem


text:
function fx = bary(x, fvals)
%BARY  Barycentric interpolation on a 2nd-kind Chebyshev grid.
%   BARY(X, FVALS) evaluates G(X) using the barycentric interpolation formula,
%   where F is the polynomial interpolant on a 2nd-kind Chebyshev grid to the
%   values stored in the columns of FVALS. X should be a column vector.
%
%   If size(FVALS, 2) > 1 then BARY returns values in the form [F_1(X), F_2(X),
%   ...], where size(F_k(X)) = size(X).
%
%   Example:
%     xcheb = chebtech2.chebpts(14);
%     fx = 1./( 1 + 25*x.^2 );
%     xx = linspace(-1, 1, 1000);
%     [xx, yy] = meshgrid(xx, xx);
%     ff = bary(xx + 1i*yy, fx);
%     h = surf(xx, yy, 0*xx, angle(-ff));
%     set(h, 'edgealpha', 0)
%     view(0,90), shg
%
% See also CHEBTECH.BARY, CHEBPTS, BARYWTS, FEVAL.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This method is basically a wrapper for BARY.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Parse inputs:
n = size(fvals, 1);

% Chebyshev nodes and barycentric weights:
xk = chebtech2.chebpts(n);
vk = chebtech2.barywts(n);

% Call BARY:
fx = bary(x, fvals, xk, vk);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
clear; close all;

%% 10MHz incident signal;
fs = 1000e6;
f = 10e6;
wavelength = fs/f;
sig = wavemaker(3.5, f, fs);
figure; subplot(3,1,1);
plot(sig, 'LineWidth', 1); title(strcat('input signal, wavelength =',num2str(wavelength),' data pts'));
axis([0,3300,-2,2])

%% Sparse signal;
N = 3000;                    % N : length of signal
s = zeros(N,1);
k = [50:(50+wavelength*0.1) 500:(500+wavelength*0.6) 1200:(1200+wavelength*1.6) 2200:(2200+wavelength*4.6)];
s(k) = 1;
subplot(3,1,2);
plot(s, 'LineWidth', 1); title('distribution: objective function');
axis([0,3300,0,2])

%% convoluted signal
y = conv(sig,s);
subplot(3,1,3);plot(y, 'LineWidth', 1);hold on;
plot(abs(hilbert(y)), 'LineWidth', 1);
title('convoluted signal between input and districution');
xlim([0,3300])

% My addition
vX = s(:);
vH = sig(:);
vY = y(:);

save('UserData', 'vX', 'vY', 'vH');
% End of my addition

function x = wavemaker(nCycles, fc, fs)
% function to generate wave packet;
nSample = round(fs / fc * nCycles);
ts      = 1 / fs;
T       = ts * nSample;
t_max   = ts * (nSample-1);
t       = 0: ts: t_max;

x       = sin( 2 * pi * fc .* t);

x = x.*hanning(nSample)';
end

author: RoyiAvital

repo: StackExchangeCodes


text:
function out = DT_IsSeasonal(y)
% DT_IsSeasonal     A simple test of seasonality.
%
% Fits a 'sin1' model to the time series using fit function from the Curve Fitting
% Toolbox. The output is binary: 1 if the goodness of fit, R^2, exceeds 0.3 and
% the amplitude of the fitted periodic component exceeds 0.5, and 0 otherwise.
%
%---INPUTS:
% y, the input time series
%
%---OUTPUT: Binary: 1 (= seasonal), 0 (= non-seasonal)

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
%% Preliminaries
%-------------------------------------------------------------------------------
% Check a curve-fitting toolbox license is available:
BF_CheckToolbox('curve_fitting_toolbox');

% Make sure the input time series, y, is a column vector
if size(y,2) > size(y,1)
    y = y';
end
N = length(y); % length of input time series
r = (1:N)'; % range over which to fit

%-------------------------------------------------------------------------------
%% Fit a sinusoidal model using the Curve-Fitting Toolbox
%-------------------------------------------------------------------------------
[cfun,gof] = fit(r,y,'sin1'); % fits the following form: a1*sin(b1*x+c1)

%-------------------------------------------------------------------------------
%% Two conditions for determining whether time series contains periodicities:
%-------------------------------------------------------------------------------
% Condition 1: fit is ok
th_fit = 0.3; % r2 > th_fit

% Condition 2: amplitude is not too small
th_ampl = 0.5; % a1 > th_ampl

if gof.rsquare > th_fit && abs(cfun.a1 > th_ampl)
    out = 1; % test thinks the time series has strong periodicities
else
    out = 0; % test thinks the time series doesn't have any strong periodicities
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function value = p14_exact ( dim_num )

%*****************************************************************************80
%
%% P14_EXACT returns the exact integral for problem 14.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Output, real VALUE, the exact value of the integral.
%
  value = ( - 1.0 / 3.0 ) * ( 1.0 - ( -1.0 / 2.0 )^dim_num );

  return
end


text:
function [k]=fetruss(el,leng,area,c,s)
%--------------------------------------------------------------
%  Purpose:
%     Compute stiffness matrices for the 2-d truss element
%     nodal dof {u_1 v_1 u_2 v_2}
%
%  Synopsis:
%     [k]=fetruss(el,leng,area,c,s)
%
%  Variable Description:
%     k - element stiffness matrix (size of 4x4)
%     el - elastic modulus ( E )
%     leng - element length
%     area - area of truss cross-section
%----------------------------------------------------------------
% syms X1 X2 X3
k=(area*el/leng)*[c*c c*s -c*c -c*s;...
                  c*s s*s -c*s -s*s;...
                  -c*c -c*s c*c c*s;...
                  -c*s -s*s c*s s*s];

author: Sable

repo: mcbench-benchmarks

path: github-repos/MATLAB/Sable-mcbench-benchmarks/mcbench-benchmarks-ba13b2f0296ef49491b95e3f984c7c41fccdb6d8/12401-optimize-truss-by-fsd-and-slp/Optimize Truss by FSD and SLP/fetruss.m


text:
function [retY] = get_y(d,ind,fInd)

%   get_y(DATA)         returns Y matrix of the data object
%   get_y(DATA,INDEXES) returns Y matrix of the data object for given indexes

if (nargin==1)|isempty(ind),
    ind=[1:size(d.Y,1)];
end;    % return all data
if nargin < 3 | isempty(fInd),
    fInd = [1:size(d.Y,2)];
end;
if ~isempty(d.Y)
  retY=d.Y(ind,fInd);
else
  retY=[];
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% TRIDEMO  Demonstration of triangular
% mesh contouring
%
% R. Pawlowicz (rpawlowicz@eos.ubc.ca) Mar/2013

% Example
res=1000;    % Number o random points

% Create a grid and a triangulation.
x=rand(res,1)*6-3;
y=rand(res,1)*6-3;
M=delaunay(x,y);
z=peaks(x,y);

% Now, remove triangles in the middle, and a bite
% out of the upper right corner

ii=find(sqrt(x.^2+y.^2)<.6 | sqrt((x-3).^2+(y-2).^2)<1.5 );
jj=zeros(size(M,1),1);
for l=1:size(M,1);
  jj(l)=any( M(l,1)==ii | M(l,2)==ii | M(l,3)==ii );
end;
M=M(~jj,:);

clf; orient landscape;
set(gcf,'defaultaxestickdir','out','defaultaxesfontsize',16);
subplot(121);
xx=[ x(M(:,[1 2 3 1])');NaN(1,size(M,1)) ];
yy=[ y(M(:,[1 2 3 1])');NaN(1,size(M,1)) ];
plot(xx(:),yy(:));
title({'Randomly generated non-convex','triangular mesh'});
axis([-3 3 -3 3]);

subplot(122);
[CS,h]=tricontf(x,y,M,z);
set(h,'edgecolor','none');
hold on;
[CS,h]=tricont(x,y,M,z,'-k');
clabel(CS,h,'fontsize',14);
hold off;
title({'...and the ''peaks'' function','contoured over that mesh'});
axis([-3 3 -3 3]);


text:
function est = refine_grid_estimates(f_obj, grid, est_idx)
%REFINE_GRID_ESTIMATES Refines DOA estimates obtained from a grid.
%Inputs:
%   f_obj - Objective function. Its local minimums identify the DOAs.
%   grid - Grid used for the original estimation.
%   est_idx - Indices (corresponding to the grid) of the original
%             estimates.
%Output:
%   est - Refined estimates.
est = zeros(size(est_idx));
if size(grid, 1) == 1
    % 1d
    n_iter = 10;
    subgrid_size = 10;
    for kk = 1:length(est_idx)
        % k-th DOA
        % init bounds
        if est_idx(kk) > 1
            lb = grid(est_idx(kk) - 1);
        else
            lb = grid(1);
        end
        if est_idx(kk) < length(grid)
            ub = grid(est_idx(kk) + 1);
        else
            ub = grid(end);
        end
        % refine
        for ii = 1:n_iter
            % find minimum over the subgrid
            subgrid = linspace(lb, ub, subgrid_size);
            obj_vals = f_obj(subgrid);
            [~, min_idx] = min(obj_vals);
            % update bounds
            if min_idx > 1
                lb = subgrid(min_idx - 1);
            else
                lb = subgrid(1);
            end
            if min_idx < subgrid_size
                ub = subgrid(min_idx + 1);
            else
                ub = subgrid(end);
            end
        end
        est(kk) = subgrid(min_idx);
    end
else
    % 2d
    error('Not implemented.');
end
end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
function exact = p40_exact ( )

%*****************************************************************************80
%
%% P40_EXACT returns the exact integral for problem 40.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  alpha = p40_param_get ( );

  exact = ( ( 1.0 - 0.25 * pi )^( alpha + 1.0 ) ...
          + (     + 0.25 * pi )^( alpha + 1.0 ) ) ...
          / ( alpha + 1.0 );

  return
end


text:
function trl = ft_trialfun_realtime(cfg)

% FT_TRIALFUN_REALTIME can be used to segment a continuous stream of
% data in real-time. Trials are defined as [begsample endsample offset
% condition]
%
% The configuration structure can contain the following specifications
%   cfg.minsample  = the last sample number that was already considered (passed from rt_process)
%   cfg.blocksize  = in seconds. In case of events, offset is with respect to the trigger.
%   cfg.offset     = the offset wrt the 0 point. In case of no events, offset is wrt
%                    prevSample. E.g., [-0.9 1] will read 1 second blocks with
%                    0.9 second overlap
%   cfg.bufferdata = {'first' 'last'}. If 'last' then only the last block of
%                    interest is read. Otherwise, all well-defined blocks are read (default = 'first')

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

if ~isfield(cfg,'minsample'),   cfg.minsample = 0;        end
if ~isfield(cfg,'blocksize'),   cfg.blocksize = 0.1;  end
if ~isfield(cfg,'offset'),      cfg.offset = 0; end
if ~isfield(cfg,'bufferdata'),  cfg.bufferdata = 'first'; end
if ~isfield(cfg,'triggers'),    cfg.triggers = [];        end

% blocksize and offset in terms of samples
cfg.blocksize = round(cfg.blocksize * cfg.hdr.Fs);
cfg.offset = round(cfg.offset * cfg.hdr.Fs);

% retrieve trials of interest
if isempty(cfg.event) % asynchronous mode
  trl = trialfun_asynchronous(cfg);
else % synchronous mode
  trl = trialfun_synchronous(cfg);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function trl = trialfun_asynchronous(cfg)

trl = [];

prevSample = cfg.minsample;

if strcmp(cfg.bufferdata, 'last') % only get last block

  % begsample starts blocksize samples before the end
  begsample  = cfg.hdr.nSamples*cfg.hdr.nTrials - cfg.blocksize;

  % begsample should be offset samples away from the previous read
  if begsample >= (prevSample + cfg.offset)

    endsample  = cfg.hdr.nSamples*cfg.hdr.nTrials;

    if begsample < endsample && begsample > 0
      trl = [begsample endsample 0 nan];
    end
  end

else % get all blocks

  while true

    % see whether new samples are available
    newsamples = (cfg.hdr.nSamples*cfg.hdr.nTrials-prevSample);

    % if newsamples exceeds the offset plus length specified in blocksize
    if newsamples >= (cfg.offset+cfg.blocksize)

      % we do not consider samples < 1
      begsample  = max(1,prevSample+cfg.offset);
      endsample  = max(1,prevSample+cfg.offset+cfg.blocksize);

      if begsample < endsample && endsample <= cfg.hdr.nSamples*cfg.hdr.nTrials
        trl = [trl; [begsample endsample 0 nan]];
      end
      prevSample = endsample;

    else
      break;
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function trl = trialfun_synchronous(cfg)

trl = [];

% process all events
for j=1:length(cfg.event)

  if isempty(cfg.triggers)
    curtrig = cfg.event(j).value;
  else
    [m1,curtrig] = ismember(cfg.event(j).value,cfg.triggers);
  end

  if isempty(curtrig), curtrig = nan; end

  if isempty(cfg.triggers) || (~isempty(m1) && m1)
    % catched a trigger of interest

    % we do not consider samples < 1
    begsample = max(1,cfg.event(j).sample + cfg.offset);
    endsample = max(1,begsample + cfg.blocksize);

    trl = [trl; [begsample endsample cfg.offset curtrig]];

  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function sdn = utc2sdn(utc)
% UTC2SDN   - Convert UTC seconds to matlab date format
%
% Use as: sdn = utc2sdn(utc);
%
% sdn = Matlab datenumber (se DATENUM and DATESTR)
% utc = Seconds since 01 Jan 1970 00:00:00

% Brian Schlining
% 12 Apr 2000

%datenum('01 Jan 1970 00:00:00') = 719529
sdn = utc/60/60/24 + 719529;

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function [r, pol, res, zer, zj, fj, wj, errvec, wt] = aaa(F, varargin)
%AAA   AAA and AAA-Lawson (near-minimax) real or complex rational approximation.
%   R = AAA(F, Z) computes the AAA rational approximant R (function handle) to
%   data F on the set of sample points Z.  F may be given by its values at Z,
%   or as a function handle or chebfun.  R = AAA(F, Z, 'degree', N) attempts to
%   compute the minimax approximation of degree N (i.e., rational type (N,N)).
%
%   [R, POL, RES, ZER] = AAA(F, Z) returns vectors of poles, residues, and zeros
%   of R.
%
%   [R, POL, RES, ZER, ZJ, FJ, WJ] = AAA(F, Z) also returns the vectors of
%   support points ZJ, approximation values FJ = r(ZJ), and weights WJ of the
%   barycentric representation of R.
%
%   [R, POL, RES, ZER, ZJ, FJ, WJ, ERRVEC] = AAA(F, Z) also returns the vector
%   of errors ||f-r||_infty in successive iterative steps of AAA.  Note that the
%   rational degrees are not 1,...,length(ERRVEC) but 0,...,length(ERRVEC)-1.
%
%   R = AAA(F, Z, NAME, VALUE) sets the following parameters:
%   - 'tol', TOL: relative tolerance (default TOL = 1e-13),
%   - 'degree', N: maximal degree (default N = 99).
%       The output rational approximant will be of degree at most N.  Like
%       'mmax', N+1, except that Lawson is turned on by default: see below.
%   - 'mmax', MMAX: maximal number of terms in the barycentric representation
%       (default MMAX = 100). R will be of degree at most MMAX-1.  Like
%       'degree', MMAX-1, except that Lawson is not turned on by default.
%   - 'dom', DOM: domain (default DOM = [-1, 1]).  No effect if Z is provided.
%   - 'cleanup', 'off' or 0: turns off automatic removal of Froissart doublets.
%   - 'cleanuptol', CLEANUPTOL: cleanup tolerance (default CLEANUPTOL = TOL).
%       Poles with residues less than this number times the geometric mean size
%       of F times the minimum distance to Z are deemed spurious by the cleanup
%       procedure.  If TOL = 0, then CLEANUPTOL defaults to 1e-13.
%   - 'lawson', NLAWSON: take NLAWSON iteratively reweighted least-squares steps
%       to bring approximation closer to minimax.  Specifying NLAWSON = 0
%       ensures there is no Lawson iteration.  See next paragraph.
%
%   If 'degree' is specified and 'lawson' is not, AAA attempts to find a minimax
%   approximant of degree N by AAA-Lawson iteration.  This will generally be
%   successful only if the minimax error is well above machine precision, and
%   is more reliable for complex problems than real ones.  If 'degree' and
%   'lawson' are both specified, then exactly NLAWSON Lawson steps are taken
%   (so NLAWSON = 0 corresponds to AAA approximation with no Lawson iteration).
%   The final weight vector WT of the Lawson iteration is available with
%   [R, POL, RES, ZER, ZJ, FJ, WJ, ERRVEC, WT] = AAA(F, Z).
%
%   Note that R may have fewer than N poles and zeros.  This may happen, for
%   example, if N is too large, or if F is even and N is odd, or if F is odd
%   and N is even.
%
%   One can also execute R = AAA(F), with no specification of a set Z.  If F is
%   a vector, this is equivalent to R = AAA(F, Z) with
%   Z = LINSPACE(-1, 1, LENGTH(F)).  If F is a function handle, AAA attempts
%   to resolve F on [-1,1] by default.
%
%   This standalone code works in GNU Octave as well as MATLAB.
%
% Examples:
%   r = aaa(@exp); xx = linspace(-1,1); plot(xx,r(xx)-exp(xx))
%
%   r = aaa(@exp,'degree',4); xx = linspace(-1,1); plot(xx,r(xx)-exp(xx))
%
%   X = linspace(-1,1,30); r = aaa(gamma(X),X);
%   fplot(r,[-5,5]), axis([-5 5 -15 15]), grid on
%
%   Z = exp(2i*pi*linspace(0,1,500));
%   [r,pol,res] = aaa(@tan,Z); disp([pol res])
%
%   X = linspace(-1,1,1000); F = tanh(20*X);
%   subplot(1,2,1)
%   r = aaa(F,X,'degree',15,'lawson',0); plot(X,F-r(X)), hold on
%   r = aaa(F,X,'degree',15); plot(X,F-r(X)), hold off
%
%   Z = exp(1i*pi*linspace(-1,1,1000)); G = exp(Z);
%   subplot(1,2,2)
%   r = aaa(G,Z,'degree',3,'lawson',0); plot(G-r(Z)), axis equal, hold on
%   r = aaa(G,Z,'degree',3); plot(G-r(Z)), axis equal, hold off
%
%   References on AAA and AAA-Lawson, respectively:
%
%   [1] Y. Nakatsukasa, O. Sete, and L. N. Trefethen, "The AAA algorithm
%   for rational approximation", SIAM J. Sci. Comp. 40 (2018), A1494-A1522.
%
%   [2] Y. Nakatsukasa and L. N. Trefethen, An algorithm for real and
%   complex rational minimax approximation, SIAM J. Sci. Comp. 42 (2020),
%   A3157-A3179.
%
% See also AAATRIG, CF, CHEBPADE, MINIMAX, PADEAPPROX, RATINTERP.

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

% Parse inputs:
[F, Z, M, dom, tol, mmax, cleanup_flag, cleanup_tol, needZ, mmax_flag, ...
    nlawson, degree_flag, degree] = parseInputs(F, varargin{:});

if ( needZ )
    % Z was not provided.  Try to resolve F on its domain.
    [r, pol, res, zer, zj, fj, wj, errvec] = ...
        aaa_autoZ(F, dom, tol, mmax, cleanup_flag, cleanup_tol, mmax_flag, ...
            nlawson, degree_flag, degree);
    return
end

% Remove infinite or NaN function values and repeated entries:
toKeep = ~isinf(F) & ~isnan(F);
F = F(toKeep); Z = Z(toKeep);
[Z, uni] = unique(Z,'stable'); F = F(uni);

% Initialization for AAA iteration:
M = length(Z);
abstol = tol*norm(F, inf);                 % Absolute tolerance
J = (1:M)';
zj = []; fj = []; C = []; A = [];
errvec = [];
R = mean(F)*ones(size(J));

% AAA iteration:
for m = 1:mmax
    % Introduce next support point:
    [~, jj] = max(abs(F(J) - R(J)));       % Select next support point
    zj = [zj; Z(J(jj))];                   % Update support points
    fj = [fj; F(J(jj))];                   % Update data values
    C = [C 1./(Z - Z(J(jj)))];             % Next column of Cauchy matrix
    J(jj) = [];                            % Update index vector
    A = [A, (F-fj(end)).*C(:,end)];        % Update Loewner matrix

    % Compute weights:
    if ( length(J) >= m )                  % The usual tall-skinny case
        [~, S, V] = svd(A(J,:), 0);        % Reduced SVD
        s = diag(S);
        mm = find( s == min(s) );          % Treat case of multiple min sing val
        nm = length(mm);
        wj = V(:,mm)*ones(nm,1)/sqrt(nm);  % Aim for non-sparse wt vector
    elseif ( length(J) >= 1 )
        V = null(A(J,:));                  % Fewer rows than columns
        nm = size(V,2);
        wj = V*ones(nm,1)/sqrt(nm);        % Aim for non-sparse wt vector
    else
        wj = ones(m,1)/sqrt(m);            % No rows at all (needed for Octave)
    end

    % Compute rational approximant:
    i0 = find(wj~=0);                      % Omit columns with wj = 0
    N = C(:,i0)*(wj(i0).*fj(i0));          % Numerator
    D = C(:,i0)*wj(i0);                    % Denominator
    R = N./D;
    Dinf = isinf(D);
    R(Dinf) = F(Dinf);                     % Interpolate at supp pts with wj~=0

    % Check if converged:
    maxerr = norm(F - R, inf);
    errvec = [errvec; maxerr];
    if ( maxerr <= abstol )
        break
    end
end
maxerrAAA = maxerr;                        % Error at end of AAA

% We now enter Lawson iteration: barycentric IRLS = iteratively reweighted
% least-squares if 'lawson' is specified with NLAWSON > 0 or 'mmax' is
% specified and 'lawson' is not.  In the latter case the number of steps
% is chosen adaptively.  Note that the Lawson iteration is unlikely to be
% successful when the errors are close to machine precision.

wj0 = wj; fj0 = fj;                        % Save params in case Lawson fails
wt = NaN(M,1); wt_new = ones(M,1);
if ( nlawson > 0 )                         % Lawson iteration

    maxerrold = maxerrAAA;
    maxerr = maxerrold;
    nj = length(zj);
    A = [];
    for j = 1:nj                           % Cauchy/Loewner matrix
        A = [A 1./(Z-zj(j)) F./(Z-zj(j))];
    end
    for j = 1:nj
        [i,~] = find(Z==zj(j));            % support pt rows are special
        A(i,:) = 0;
        A(i,2*j-1) = 1;
        A(i,2*j) = F(i);
    end
    stepno = 0;
    while ( (nlawson < inf) && (stepno < nlawson) ) || ...
          ( (nlawson == inf) && (stepno < 20) ) || ...
          ( (nlawson == inf) && (maxerr/maxerrold < .999) && (stepno < 1000) )
        stepno = stepno + 1;
        wt = wt_new;
        W = spdiags(sqrt(wt),0,M,M);
        [~,~,V] = svd(W*A,0);
        c = V(:,end);
        denom = zeros(M,1); num = zeros(M,1);
        for j = 1:nj
            denom = denom + c(2*j)./(Z-zj(j));
            num = num - c(2*j-1)./(Z-zj(j));
        end
        R = num./denom;
        for j = 1:nj
            [i,~] = find(Z==zj(j));        % support pt rows are special
            R(i) = -c(2*j-1)/c(2*j);
        end
        err = F - R; abserr = abs(err);
        wt_new = wt.*abserr; wt_new = wt_new/norm(wt_new,inf);
        maxerrold = maxerr;
        maxerr = max(abserr);
    end
    wj = c(2:2:end);
    fj = -c(1:2:end)./wj;
    % If Lawson has not reduced the error, return to pre-Lawson values.
    if ( (maxerr > maxerrAAA) && (nlawson == Inf) )
        wj = wj0; fj = fj0;
    end
end

% Remove support points with zero weight:
I = find(wj == 0);
zj(I) = []; wj(I) = []; fj(I) = [];

% Construct function handle and compute poles, residues and zeros:
r = @(zz) reval(zz, zj, fj, wj);
[pol, res, zer] = prz(zj, fj, wj);

if ( cleanup_flag == 1 && nlawson == 0 )      % Remove Froissart doublets
    [r, pol, res, zer, zj, fj, wj] = ...
        cleanup(r, pol, res, zer, zj, fj, wj, Z, F, cleanup_tol);
elseif ( cleanup_flag == 2 && nlawson == 0 )  % Alternative cleanup.  Currently
    a.zj = zj; a.fj = fj; a.wj = wj;          % an undocumented feature,
    a.Z = Z; a.F = F;                         % pending further investigation.
    a.cleanup_tol = max(cleanup_tol, eps);
    c = cleanup2(a);
    zj = c.zj; fj = c.fj; wj = c.wj;
    r = @(zz) reval(zz, zj, fj, wj);
    [pol, res, zer] = prz(zj, fj, wj);
end

end % of AAA()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   PARSEINPUTS   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [F, Z, M, dom, tol, mmax, cleanup_flag, cleanup_tol, ...
    needZ, mmax_flag, nlawson, degree_flag, degree] = parseInputs(F, varargin)

% Check if F is empty:
if ( isempty(F) )
    error('AAA:emptyF', 'No function given.')
elseif ( isa(F, 'chebfun') )
    if ( size(F, 2) ~= 1 )
        error('AAA:nColF', 'Input chebfun must have one column.')
    end
end

% Sample points:
if ( ~isempty(varargin) && isfloat(varargin{1}) )
    % Z is given.
    Z = varargin{1};
    if ( isempty(Z) )
        error('AAA:emptyZ', ...
            'If sample set is provided, it must be nonempty.')
    end
    varargin(1) = [];
end

% Set defaults for other parameters:
tol = 1e-13;                   % Relative tolerance
mmax = 100;                    % Maximum number of terms
degree = NaN;                  % Specified degree
cleanup_tol = 1e-13;           % Cleanup tolerance
nlawson = Inf;                 % Number of Lawson steps (Inf means adaptive)
% Domain:
if ( isa(F, 'chebfun') )
    dom = F.domain([1, end]);
else
    dom = [-1, 1];
end
cleanup_flag = 1;              % Cleanup on
mmax_flag = 0;                 % Checks if mmax manually specified
degree_flag = 0;               % Checks if degree specified
cleanup_set = 0;               % Checks if cleanup_tol manually specified
while ( ~isempty(varargin) )   % Check if parameters have been provided
    if ( strncmpi(varargin{1}, 'tol', 3) )
        if ( isfloat(varargin{2}) && isequal(size(varargin{2}), [1, 1]) )
            tol = varargin{2};
            if ( ~cleanup_set && tol > 0 ) % If not set, set cleanup_tol to tol
              cleanup_tol = tol;
            end
        end
        varargin([1, 2]) = [];

    elseif ( strncmpi(varargin{1}, 'degree', 6) )
        if ( isfloat(varargin{2}) && isequal(size(varargin{2}), [1, 1]) )
            if ( mmax_flag == 1 ) && ( mmax ~= varargin{2}+1 )
                error('AAA:degmmaxmismatch', ' mmax must equal degree+1.')
            end
            degree = varargin{2};
            mmax = degree + 1;
            mmax_flag = 1;
            degree_flag = 1;
        end
        varargin([1, 2]) = [];

    elseif ( strncmpi(varargin{1}, 'mmax', 4) )
        if ( isfloat(varargin{2}) && isequal(size(varargin{2}), [1, 1]) )
            if ( mmax_flag == 1 ) && ( mmax ~= varargin{2})
                error('AAA:degmmaxmismatch', ' mmax must equal degree+1.')
            end
            mmax = varargin{2};
            mmax_flag = 1;
        end
        varargin([1, 2]) = [];

    elseif ( strncmpi(varargin{1}, 'lawson', 6) )
        if ( isfloat(varargin{2}) && isequal(size(varargin{2}), [1, 1]) )
            nlawson = varargin{2};
        end
        varargin([1, 2]) = [];

    elseif ( strncmpi(varargin{1}, 'dom', 3) )
        if ( isfloat(varargin{2}) && isequal(size(varargin{2}), [1, 2]) )
            dom = varargin{2};
        end
        varargin([1, 2]) = [];
        if ( isa(F, 'chebfun') )
            if ( ~isequal(dom, F.domain([1, end])) )
                warning('AAA:dom', ...
                    ['Given domain does not match that of the chebfun.\n', ...
                    'Results may be inaccurate.'])
            end
        end

    elseif ( strncmpi(varargin{1}, 'cleanuptol', 10) )
        if ( isfloat(varargin{2}) && isequal(size(varargin{2}), [1, 1]) )
          cleanup_tol = varargin{2};
          cleanup_set = 1;
        end
        varargin([1, 2]) = [];

    elseif ( strncmpi(varargin{1}, 'cleanup', 7) )
        if ( strncmpi(varargin{2}, 'off', 3) || ( varargin{2} == 0 ) )
            cleanup_flag = 0;
        elseif ( varargin{2} == 2 )     % Alternative cleanup
            cleanup_flag = 2;
        end
        varargin([1, 2]) = [];

    else
        error('AAA:UnknownArg', 'Argument unknown.')
    end
end

% Deal with Z and F:
if ( ~exist('Z', 'var') && isfloat(F) )
    % F is given as data values, pick same number of sample points:
    Z = linspace(dom(1), dom(2), length(F)).';
end

if ( exist('Z', 'var') )
    % Z is given:
    needZ = 0;

    % Work with column vector:
    Z = Z(:);
    M = length(Z);

    % Function values:
    if ( isa(F, 'function_handle') || isa(F, 'chebfun') )
        % Sample F on Z:
        F = F(Z);
    elseif ( isnumeric(F) )
        % Work with column vector and check that it has correct length.
        F = F(:);
        if ( length(F) ~= M )
            error('AAA:lengthFZ', ...
                'Inputs F and Z must have the same length.')
        end
    elseif ( ischar(F) )
        % F is given as a string input. Convert it to a function handle.
        F = inline(vectorize(F));
        F = F(Z);
    else
        error('AAA:UnknownF', 'Input for F not recognized.')
    end

else
    % Z was not given.  Set flag that Z needs to be determined.
    % Also set Z and M since they are needed as output.
    needZ = 1;
    Z = [];
    M = length(Z);
end

if ( ~degree_flag && (nlawson == Inf) )
    nlawson = 0;
end

end % End of PARSEINPUTS.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   CLEANUP   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [r, pol, res, zer, z, f, w] = ...
    cleanup(r, pol, res, zer, z, f, w, Z, F, cleanup_tol)
% Remove spurious pole-zero pairs.

% Find negligible residues:
if any(F)
    geometric_mean_of_absF = exp(mean(log(abs(F(F~=0)))));
else
    geometric_mean_of_absF = 0;
end
Zdistances = NaN(size(pol));
for j = 1:length(Zdistances)
    Zdistances(j) = min(abs(pol(j)-Z));
end
ii = find(abs(res)./Zdistances < cleanup_tol * geometric_mean_of_absF);
ni = length(ii);
if ( ni == 0 )
    return
elseif ( ni == 1 )
    warning('AAA:Froissart','1 Froissart doublet');
else
    warning('AAA:Froissart',[int2str(ni) ' Froissart doublets']);
end

% For each spurious pole find and remove closest support point:
for j = 1:ni
    azp = abs(z-pol(ii(j)));
    jj = find(azp == min(azp),1);

    % Remove support point(s):
    z(jj) = []; f(jj) = [];
end

% Remove support points z from sample set:
for jj = 1:length(z)
    F(Z == z(jj)) = [];
    Z(Z == z(jj)) = [];
end
m = length(z);
M = length(Z);

% Build Loewner matrix:
SF = spdiags(F, 0, M, M);
Sf = diag(f);
C = 1./(Z-z.');               % Cauchy matrix.
A = SF*C - C*Sf;              % Loewner matrix.

% Solve least-squares problem to obtain weights:
[~, ~, V] = svd(A, 0);
w = V(:,m);

% Build function handle and compute poles, residues and zeros:
r = @(zz) reval(zz, z, f, w);
[pol, res, zer] = prz(z, f, w);

end % End of CLEANUP.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   CLEANUP2   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c = cleanup2(a)
%% Alternative cleanup procedure to remove spurious pole-zero pairs.
% This considers pole-zero distances.  Stefano Costa, August 2022.

z = a.zj; f = a.fj; w = a.wj;
[pol, res, zer] = prz(z, f, w);

cleanup_tol = a.cleanup_tol;

niter = 0;
while(true)
    niter = niter+1;
    Z = a.Z; F = a.F;
    ii = [];
    for jj = 1:length(pol)
        dz = min(abs(zer-pol(jj))); if isempty(dz), dz = 1e100; end
        dS = abs(Z-pol(jj));
        ds = min(dS);
        if any(F)
            q = 4*pi*abs(F).*dS;
            Q = mean(q);                % Arithmetic mean
        else
            Q = 0;
        end
        R = 8*cleanup_tol*Q/(4*pi);    % Equivalent residue value

        % Conditions to expunge poles
        % Expunge if either minimum distance is zero
        if (ds==0) || (dz==0)
            ii = [ii; jj];
        % Expunge if Z is a real interval
        elseif isreal(Z) && (abs(imag(pol(jj)))<eps) && ...
            (real(pol(jj))>=min(Z)) && (real(pol(jj))<=max(Z))
            ii = [ii; jj];
        % Expunge if Z is the unit disk
        elseif all(abs(Z)==1) && (abs(abs(pol(jj))-1)<eps)
            ii = [ii; jj];
        % Expunge if distance to closest zero is undetectable
        elseif (dz/ds<1) && (dz<max(cleanup_tol^2,eps))
            ii = [ii; jj];
        % Expunge if a nearby zero exists and residue is below the
        % equivalent value R. Two choices for real and complex F
        elseif ((dz/ds)<sqrt(cleanup_tol))
            if ( ~any(imag(F)) && (abs(real(res(jj))) < R) )
                ii = [ii; jj];
            elseif (abs(res(jj)) < R)
                ii = [ii; jj];
            end
        end
    end
    ii = unique(ii);

    ni = length(ii);
    if ( ni == 0 )
        % Nothing to do.
        break;
    elseif ( ni == 1 )
        warning('AAA:Froissart',...
            ['1 Froissart doublet, niter = ', int2str(niter)]);
    else
        warning('AAA:Froissart',...
            [int2str(ni) ' Froissart doublets, niter = ' int2str(niter)]);
    end

    % For each spurious pole find and remove closest support point:
    for j = 1:ni
        azp = abs(z-pol(ii(j)));
        jj = find(azp == min(azp),1);

        % Remove support point(s):
        z(jj) = [];
        f(jj) = [];
    end

    % Remove support points z from sample set:
    for jj = 1:length(z)
        F(Z == z(jj)) = [];
        Z(Z == z(jj)) = [];
    end
    m = length(z);
    M = length(Z);

    % Build Loewner matrix:
    SF = spdiags(F, 0, M, M);
    Sf = diag(f);
    C = 1./(Z-z.');             % Cauchy matrix.
    A = SF*C - C*Sf;            % Loewner matrix.

    % Solve least-squares problem to obtain weights:
    [~, ~, V] = svd(A, 0);
    w = V(:,m);

    % Compute poles, residues and zeros:
    [pol, res, zer] = prz(z, f, w);
end % End of while loop

c.zj = z; c.fj = f; c.wj = w;

end  % End of CLEANUP2.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   AAA_AUTOZ   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [r, pol, res, zer, zj, fj, wj, errvec] = ...
    aaa_autoZ(F, dom, tol, mmax, cleanup_flag, cleanup_tol, mmax_flag, ...
               nlawson, degree_flag, degree)
% Automated choice of sample set

% Flag if function has been resolved:
isResolved = 0;

% Main loop:
for n = 5:14
    % Sample points:
    % Next line enables us to do pretty well near poles
    Z = linspace(dom(1)+1.37e-8*diff(dom), dom(2)-3.08e-9*diff(dom), 1 + 2^n).';
    if degree_flag
       [r, pol, res, zer, zj, fj, wj, errvec] = aaa(F, Z, 'tol', tol, ...
          'mmax', mmax, 'cleanup', cleanup_flag, 'cleanuptol', cleanup_tol, ...
          'lawson', nlawson, 'degree', degree);
    else
       [r, pol, res, zer, zj, fj, wj, errvec] = aaa(F, Z, 'tol', tol, ...
          'mmax', mmax, 'cleanup', cleanup_flag, 'cleanuptol', cleanup_tol, ...
          'lawson', nlawson);
    end
    % Test if rational approximant is accurate:
    abstol = tol * norm(F(Z), inf);

    % On Z(n):
    err(1,1) = norm(F(Z) - r(Z), inf);

    Zrefined = linspace(dom(1)+1.37e-8*diff(dom), dom(2)-3.08e-9*diff(dom), ...
        round(1.5 * (1 + 2^(n+1)))).';
    err(2,1) = norm(F(Zrefined) - r(Zrefined), inf);
    if ( all(err < abstol) )
        % Final check that the function is resolved, inspired by sampleTest().
        % Pseudo random sample points in [-1, 1]:
        xeval = [-0.357998918959666; 0.036785641195074];
        % Scale to dom:
        xeval = (dom(2) - dom(1))/2 * xeval + (dom(2) + dom(1))/2;

        if ( norm(F(xeval) - r(xeval), inf) < abstol )
            isResolved = 1;
            break
        end
    end
end

if ( ( isResolved == 0 ) && ~mmax_flag )
    warning('AAA:notResolved', ...
        'Function not resolved using %d pts.', length(Z))
end

end % End of AAA_AUTOZ.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   PRZ   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [pol, res, zer] = prz(zj, fj, wj)
%   Compute poles, residues, and zeros of rational fun in barycentric form.

% Compute poles via generalized eigenvalue problem:
m = length(wj);
B = eye(m+1);
B(1,1) = 0;
E = [0 wj.'; ones(m, 1) diag(zj)];
pol = eig(E, B);
pol = pol(~isinf(pol));

% Compute residues via formula for res of quotient of analytic functions:
N = @(t) (1./(t-zj.')) * (fj.*wj);
Ddiff = @(t) -((1./(t-zj.')).^2) * wj;
res = N(pol)./Ddiff(pol);

% Compute zeros via generalized eigenvalue problem:
E = [0 (wj.*fj).'; ones(m, 1) diag(zj)];
zer = eig(E, B);
zer = zer(~isinf(zer));

end % End of PRZ.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   REVAL   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function r = reval(zz, zj, fj, wj)
%   Construct function handle to evaluate rational function in barycentric form.

zv = zz(:);                         % vectorize zz if necessary
CC = 1./(zv-zj.');                  % Cauchy matrix
r = (CC*(wj.*fj))./(CC*wj);         % vector of values

% Deal with input inf: r(inf) = lim r(zz) = sum(w.*f) / sum(w):
r(isinf(zv)) = sum(wj.*fj)./sum(wj);

% Deal with NaN:
ii = find(isnan(r));
for jj = 1:length(ii)
    if ( isnan(zv(ii(jj))) || ~any(zv(ii(jj)) == zj) )
        % r(NaN) = NaN is fine.
        % The second case may happen if r(zv(ii)) = 0/0 at some point.
    else
        % Clean up values NaN = inf/inf at support points.
        % Find the corresponding node and set entry to correct value:
        r(ii(jj)) = fj(zv(ii(jj)) == zj);
    end
end

% Reshape to input format:
r = reshape(r, size(zz));

end % End of REVAL.


text:
function [A,B,flag] = oprACAv(opr,k,X,Xnrm,Y,Ynrm,tol)
%+========================================================================+
%|                                                                        |
%|            OPENOPR - LIBRARY FOR SPECIFIC OPERATORS IN BEM             |
%|           openOpr is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2019.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : oprACAv.m                                     |
%|    #    |   VERSION    : 0.61                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 05.09.2019                                    |
%| ( === ) |   SYNOPSIS   : Adaptative Cross Approximation, partial       |
%|  `---'  |                & total pivoting                              |
%+========================================================================+

% Matrix dimensions
Nx = size(X,1);
Ny = size(Y,1);

% Initialize indices
Ir = (1:Nx)';
Ic = (1:Ny)';

% Partial pivoting with handle function
row = @(i) oprGreenKernel(opr,k,X(i,:),Xnrm(i,:),Y,Ynrm);
col = @(j) oprGreenKernel(opr,k,X,Xnrm,Y(j,:),Ynrm(j,:));

% First row (row index of maximum reference value)
i   = 1;
B   = row(i);
dim = size(B,3);

% First pivot (maximum of the first row)
[~,j] = max(max(abs(B),[],3));
delta = B(j,:,:);

% Securitys
delta(abs(delta)<1e-12) = 1e-12;

% First column
A = col(j)./delta;

% Unused indices1
Ir = [Ir(1:i-1);Ir(i+1:end)];
Ic = [Ic(1:j-1);Ic(j+1:end)];

% Frobenius norm for the initial tensor product
Rn2 = zeros(dim,1);
err = 0;
for l = 1:dim
    An2    = A(:,1,l)'*A(:,1,l);
    Bn2    = B(:,1,l)'*B(:,1,l);
    Rn2(l) = An2*Bn2;
    err    = max(err,sqrt(An2)*sqrt(Bn2)/sqrt(Rn2(l)));
end

% Iterative construction using frobenius norm
n = 1;
while (err > tol)
    % New pivot by maximum
    [~,i] = max(max(abs(A(Ir,n,:)),[],3));
    new   = row(Ir(i));
    for l = 1:dim
        new(:,:,l) = new(:,:,l) - B(:,:,l)*A(Ir(i),1:n,l).';
    end
    [~,j] = max(max(abs(new(Ic,:,:)),[],3));
    delta = new(Ic(j),:,:);

    % New pivot by minimum
    [~,iMin] = min(min(abs(A(Ir,n,:)),[],3));
    newMin   = row(Ir(iMin));
    for l = 1:dim
        newMin(:,:,l) = newMin(:,:,l) - B(:,:,l)*A(Ir(iMin),:,l).';
    end
    [~,jMin] = max(max(abs(newMin(Ic,:,:)),[],3));
    deltaMin = newMin(Ic(jMin),:,:);

    % Best pivot
    if (norm(deltaMin(:)) > norm(delta(:)))
        i     = iMin;
        new   = newMin;
        j     = jMin;
        delta = deltaMin;
    end

    % Security
    delta(abs(delta)<1e-12) = 1e-12;

    % Update row
    B(:,n+1,:) = new;

    % New column
    new = col(Ic(j));
    for l = 1:dim
        new(:,:,l) = (new(:,:,l) - A(:,:,l)*B(Ic(j),1:n,l).')./delta(l);    %B(:,:,l)*A(Ir(i),:,l).';
    end
    A(:,n+1,:) = new;

    % Update row and column indices
    Ir = [Ir(1:i-1);Ir(i+1:end)];
    Ic = [Ic(1:j-1);Ic(j+1:end)];

    % Incrementation
    n = n + 1;

    % Relative Frobenius error by block
    err = 0;
    for l = 1:dim
        An2    = A(:,n,l)'*A(:,n,l);
        Bn2    = B(:,n,l)'*B(:,n,l);
        u      = B(:,n,l)'*B(:,1:n-1,l);
        v      = A(:,n,l)'*A(:,1:n-1,l);
        AB     = v*u.';
        Rn2(l) = Rn2(l) + 2*real(AB) + An2*Bn2;
        err    = max(err,sqrt(An2)*sqrt(Bn2)/sqrt(Rn2(l)));
    end

    % Compression failed
    if (n*(Nx+Ny) > Nx*Ny)
        A    = [];
        B    = [];
        flag = 0;
        return
    end
%     sqrt(An2)*sqrt(Bn2)/sqrt(Rn2)
%     norm(A*B.' - A(:,1:n-1)*B(:,1:n-1).','fro')/norm(A*B.','fro')
%     norm(A(:,n)*B(:,n).','fro')/norm(A*B.','fro')
%     norm(ref,'inf')/nrm
%     '================================'
end

% B transposition lead to  A * B
B    = permute(B,[2 1,3]);
flag = 1;
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function WHVLoss = CalWHVLoss(PopObj,FrontNo,wz,AA,RA)
% Calculate the weighted hypervolume (WHV) loss of each solution

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

    %% Calculate the weight of each solution for WHV calculation
    Weight = ones(1,size(PopObj,1));
    if nargin > 2 && ~isempty(wz)
        for i = 1 : size(PopObj,1)
            if any(all(repmat(PopObj(i,:),size(AA,1),1)<=AA,2))
                % The solution belongs to the preferred region (Pr)
                Weight(i) = wz(3);
            elseif any(all(repmat(PopObj(i,:),size(RA,1),1)>RA,2))
                % The solution belongs to the dominated region (Do)
                Weight(i) = wz(1);
            else
                % The solution belongs to the no preference information
                % region (In)
                Weight(i) = wz(2);
            end
        end
    end

    %% Calculate the WHV loss of each solution front by front
    WHVLoss  = zeros(1,size(PopObj,1));
    RefPoint = max(PopObj,[],1) + 0.1;
    for f = setdiff(unique(FrontNo),inf)
        current  = find(FrontNo==f);
        totalWHV = CalWHV(PopObj(current,:),RefPoint,Weight(current));
        for i = 1 : length(current)
            drawnow('limitrate');
            currenti           = current([1:i-1,i+1:end]);
            WHVLoss(current(i))= totalWHV - CalWHV(PopObj(currenti,:),RefPoint,Weight(currenti));
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function H = calc_Hazard(H0,betaV,dkV)
%Hazard calculated using Cox model
%AI 04/20/22

H = H0*exp(sum(betaV.*dkV));

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% 股票预测问题 by 曲线拟合
%% 历史数据处理
clear;clc;
x = [2,3,4,5,8,9,10,11,12,15,16,17,18,19,22,23,24,25,26,29,30]; % 时间
y = [7.74,7.84,7.82,7.78,7.91, 7.97,7.9,7.76,7.9,8.04,8.06, 8.11,8.08,8.13,8.03,8.01,8.06,8.0,8.3,8.41,8.28];  % 当天收盘价

p=polyfit(x,y,3);  % 多项式拟合，返回次数为 3 的多项式 p(x) 的系数，该阶数是 y 中数据的最佳拟合（在最小二乘方式中）。p 中的系数按降幂排列，p 的长度为 n+1
plot(x,y,'*',x,polyval(p,x));

%% 开始预测
x1=[31,32,33];   % 待预测时间
xi=[x,x1];
y1=[8.27,8.17,9.54];  %  该股票后三个交易日的收盘价分别为8.27,8.17,9.54
plot(x,y,'*',xi,polyval(p,xi),x1,y1,'rp');
% 结论：预测结果仅供参考

repo: NumericalAnalysis


text:
function dcor = distcorr(x,y)

% This function calculates the distance correlation between x and y.
% Reference: http://en.wikipedia.org/wiki/Distance_correlation
% Date: 18 Jan, 2013
% Author: Shen Liu (shen.liu@hotmail.com.au)

% Check if the sizes of the inputs match
if size(x,1) ~= size(y,1);
    error('Inputs must have the same number of rows')
end

% Delete rows containing unobserved values
N = any([isnan(x) isnan(y)],2);
x(N,:) = [];
y(N,:) = [];

% Calculate doubly centered distance matrices for x and y
a = pdist2(x,x);
mcol = mean(a);
mrow = mean(a,2);
ajbar = ones(size(mrow))*mcol;
akbar = mrow*ones(size(mcol));
abar = mean(mean(a))*ones(size(a));
A = a - ajbar - akbar + abar;

b = pdist2(y,y);
mcol = mean(b);
mrow = mean(b,2);
bjbar = ones(size(mrow))*mcol;
bkbar = mrow*ones(size(mcol));
bbar = mean(mean(b))*ones(size(b));
B = b - bjbar - bkbar + bbar;

% Calculate squared sample distance covariance and variances
dcov = sum(sum(A.*B))/(size(mrow,1)^2);

dvarx = sum(sum(A.*A))/(size(mrow,1)^2);
dvary = sum(sum(B.*B))/(size(mrow,1)^2);

% Calculate the distance correlation
dcor = sqrt(dcov/sqrt(dvarx*dvary));

repo: REST


text:
function x = normalize(x, mu, sigma)
    x=bsxfun(@minus,x,mu);
	x=bsxfun(@rdivide,x,sigma);
end

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function g = sigmoidGradient(z)
%SIGMOIDGRADIENT returns the gradient of the sigmoid function
%evaluated at z
%   g = SIGMOIDGRADIENT(z) computes the gradient of the sigmoid function
%   evaluated at z. This should work regardless if z is a matrix or a
%   vector. In particular, if z is a vector or matrix, you should return
%   the gradient for each element.

g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the gradient of the sigmoid function evaluated at
%               each value of z (z can be a matrix, vector or scalar).

g = sigmoid(z) .* (1 - sigmoid(z));

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

end

author: zsiciarz

repo: ml-coursera


text:
function net = init_MSLapSRN_model(opts, mode)
% -------------------------------------------------------------------------
%   Description:
%       initialize MS-LapSRN model
%
%   Input:
%       - opts  : options generated from init_MSLapSRN_opts()
%
%   Output:
%       - net   : dagnn model
%
%   Citation:
%       Fast and Accurate Image Super-Resolution with Deep Laplacian Pyramid Networks
%       Wei-Sheng Lai, Jia-Bin Huang, Narendra Ahuja, and Ming-Hsuan Yang
%       arXiv, 2017
%
%   Contact:
%       Wei-Sheng Lai
%       wlai24@ucmerced.edu
%       University of California, Merced
% -------------------------------------------------------------------------

    %% parameters
    rng('default');
    rng(0) ;

    f   = opts.conv_f;
    n   = opts.conv_n;
    pad = floor(f/2);

    if( f == 3 )
        crop = [0, 1, 0, 1];
    elseif( f == 5 )
        crop = [1, 2, 1, 2];
    else
        error('Need to specify crop in deconvolution for f = %d\n', f);
    end

    %% initialize model
    net = dagnn.DagNN;

    %% multiscale training
    for s = 1:length(opts.scales)

        scale = opts.scales(s);
        level = ceil(log(scale) / log(2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Feature extraction branch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        %% input conv
        inputs  = { sprintf('x%dSR_LR', scale) };
        outputs = { sprintf('x%dSR_input_conv', scale) };
        params  = { 'input_conv_f', 'input_conv_b' };

        net.addLayer(outputs{1}, ...
                     dagnn.Conv('size', [f, f, 1, n], ...
                                'pad', pad, ...
                                'stride', 1), ...
                     inputs, outputs, params);

        level_input = outputs{1};

        %% feature embedding sub-network
        for l = 1 : level

            block_input = {};
            block_output = {};

            for r = 1:opts.recursive

                if r == 1
                    block_input{r} = level_input;
                else
                    block_input{r} = block_output{r-1};
                end

                %% recursive block
                for d = 1:opts.depth

                    if d == 1
                        next_input = block_input{r};
                    end

                    % ReLU
                    inputs  = { next_input };
                    outputs = { sprintf('x%dSR_level%d_R%d_relu%d', scale, l, r, d) };

                    net.addLayer(outputs{1}, ...
                                 dagnn.ReLU('leak', 0.2), ...
                             inputs, outputs);

                    next_input = outputs{1};

                    % conv
                    inputs  = { next_input };
                    outputs = { sprintf('x%dSR_level%d_R%d_conv%d', scale, l, r, d) };
                    params  = { sprintf('conv%d_f', d), ...
                                sprintf('conv%d_b', d)};

                    net.addLayer(outputs{1}, ...
                                 dagnn.Conv('size', [f, f, n, n], ...
                                            'pad', pad, ...
                                            'stride', 1), ...
                                 inputs, outputs, params);

                    next_input = outputs{1};

                end %% end of recursive block

                %% local skip connection
                if strcmp(opts.LRL, 'NS')

                    % no skip connection
                    block_output{r} = next_input;

                elseif strcmp(opts.LRL, 'DS')

                    % next_input + block_input
                    inputs  = { next_input, block_input{r} };
                    outputs = { sprintf('x%dSR_level%d_R%d_output', scale, l, r) };
                    net.addLayer(outputs{1}, ...
                        dagnn.Sum(), ...
                        inputs, outputs);

                    block_output{r} = outputs{1};

                elseif strcmp(opts.LRL, 'SS')

                    % next_input + level_input
                    inputs  = { next_input, level_input };
                    outputs = { sprintf('x%dSR_level%d_R%d_output', scale, l, r) };
                    net.addLayer(outputs{1}, ...
                        dagnn.Sum(), ...
                        inputs, outputs);

                    block_output{r} = outputs{1};

                else

                    error('Unknown local skip connection %s.', opts.LRL);

                end %% end of local skip connection

            end %% end of recursive

            %% features upsample layers
            % ReLU
            inputs  = { block_output{opts.recursive} };
            outputs = { sprintf('x%dSR_level%d_uprelu', scale, l) };

            net.addLayer(outputs{1}, ...
                         dagnn.ReLU('leak', 0.2), ...
                         inputs, outputs);

            next_input = outputs{1};

            % conv
            inputs  = { next_input };
            outputs = { sprintf('x%dSR_level%d_upconv', scale, l) };
            params  = { 'upconv_f', 'upconv_b' };

            net.addLayer(outputs{1}, ...
                         dagnn.ConvTranspose(...
                             'size', [f, f, n, n], ...
                             'upsample', 2, ...
                             'crop', crop, ...
                             'numGroups', 1, ...
                             'hasBias', true), ...
                         inputs, outputs, params) ;

            level_input = outputs{1};

        end %% end of level (feature extraction)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Image reconstruction branch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        for l = 1 : level

            if l == 1
                next_input = sprintf('x%dSR_LR', scale);
            else
                next_input = sprintf('x%dSR_%dx_output', scale, 2^(l-1));
            end

            %% image upsample layer
            inputs  = { next_input };
            outputs = { sprintf('x%dSR_level%d_img_up', scale, l) };
            params  = { 'img_up_f' };

            net.addLayer(outputs{1}, ...
                dagnn.ConvTranspose(...
                    'size', [4, 4, 1, 1], ...
                    'upsample', 2, ...
                    'crop', 1, ...
                    'numGroups', 1, ...
                    'hasBias', false), ...
                inputs, outputs, params) ;

            %% residual prediction layer (f x f x n x 1)
            inputs  = { sprintf('x%dSR_level%d_upconv', scale, l) };
            outputs = { sprintf('x%dSR_level%d_residual', scale, l) };
            params  = { 'residual_conv_f', 'residual_conv_b' };

            net.addLayer(outputs{1}, ...
                dagnn.Conv('size', [f, f, n, 1], ...
                           'pad', pad, ...
                           'stride', 1), ...
                inputs, outputs, params);

            %% addition layer
            inputs  = { sprintf('x%dSR_level%d_img_up', scale, l), ...
                        sprintf('x%dSR_level%d_residual', scale, l) };
            outputs = { sprintf('x%dSR_%dx_output', scale, 2^l) };
            net.addLayer(outputs{1}, ...
                dagnn.Sum(), ...
                inputs, outputs);

            next_input = outputs{1};

            %% Loss layer
            inputs  = { next_input, ...
                        sprintf('x%dSR_%dx_HR', scale, 2^l) };
            outputs = { sprintf('x%dSR_%dx_%s_loss', scale, 2^l, opts.loss) };

            net.addLayer(outputs{1}, ...
                     dagnn.vllab_dag_loss(...
                        'loss_type', opts.loss), ...
                     inputs, outputs);

        end %% end of level (image reconstruction)

    end %% end of multiscale

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% initialize parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if strcmp(mode, 'train')

        %% input conv
        params  = { 'input_conv_f', 'input_conv_b' };

        sigma   = opts.init_sigma;
        filters = sigma * randn(f, f, 1, n, 'single');
        biases  = zeros(1, n, 'single');

        idx = net.getParamIndex(params{1});
        net.params(idx).value         = filters;
        net.params(idx).learningRate  = 1;
        net.params(idx).weightDecay   = 1;

        idx = net.getParamIndex(params{2});
        net.params(idx).value         = biases;
        net.params(idx).learningRate  = 0.1;
        net.params(idx).weightDecay   = 1;

        %% deep conv
        for d = 1:opts.depth

            params  = { sprintf('conv%d_f', d), ...
                        sprintf('conv%d_b', d)};

            sigma   = sqrt( 2 / (f * f * n) );
            filters = sigma * randn(f, f, n, n, 'single');
            biases  = zeros(1, n, 'single');

            idx = net.getParamIndex(params{1});
            net.params(idx).value         = filters;
            net.params(idx).learningRate  = 1;
            net.params(idx).weightDecay   = 1;

            idx = net.getParamIndex(params{2});
            net.params(idx).value         = biases;
            net.params(idx).learningRate  = 0.1;
            net.params(idx).weightDecay   = 1;

        end

        %% feature upsample
        params  = { 'upconv_f', 'upconv_b' };

        sigma   = sqrt( 2 / (f * f * n) );
        filters = sigma * randn(f, f, n, n, 'single');
        biases  = zeros(1, n, 'single');

        idx = net.getParamIndex(params{1});
        net.params(idx).value         = filters;
        net.params(idx).learningRate  = 1;
        net.params(idx).weightDecay   = 1;

        idx = net.getParamIndex(params{2});
        net.params(idx).value         = biases;
        net.params(idx).learningRate  = 0.1;
        net.params(idx).weightDecay   = 1;

        %% image upsample
        params  = { 'img_up_f' };

        filters = single(bilinear_kernel(4, 1, 1));

        idx = net.getParamIndex(params{1});
        net.params(idx).value         = filters;
        net.params(idx).learningRate  = 1;
        net.params(idx).weightDecay   = 1;

        %% residual prediction
        params  = { 'residual_conv_f', 'residual_conv_b' };

        sigma   = sqrt(2 / (f * f * n));
        filters = sigma * randn(f, f, n, 1, 'single');
        biases  = zeros(1, 1, 'single');

        idx = net.getParamIndex(params{1});
        net.params(idx).value         = filters;
        net.params(idx).learningRate  = 1;
        net.params(idx).weightDecay   = 1;

        idx = net.getParamIndex(params{2});
        net.params(idx).value         = biases;
        net.params(idx).learningRate  = 0.1;
        net.params(idx).weightDecay   = 1;

    end

end

repo: LapSRN


text:
function h = p12_h ( n, x )

%*****************************************************************************80
%
%% P12_H evaluates the Hessian for problem 12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real H(N,N), the N by N Hessian matrix.
%
  h = zeros ( n, n );

  d1 = 2.0 / 3.0;

  for i = 1 : 99

    arg = i / 100.0;
    r = ( - 50.0 * log ( arg ) )^d1 + 25.0 - x(2);
    t1 = abs ( r )^x(3) / x(1);
    t2 = exp ( - t1 );
    t3 = t1 * t2 * ( t1 * t2 + ( t1 - 1.0 ) * ( t2 - arg ) );
    t = t1 * t2 * ( t2 - arg );
    logr = log ( abs ( r ) );

    h(1,1) = h(1,1) + 2.0 * t3 - 2.0 * t;
    h(1,2) = h(1,2) + 2.0 * t3 / r;
    h(1,3) = h(1,3) - 2.0 * t3 * logr;

    h(2,1) = h(2,1) + 2.0 * t3 / r;
    h(2,2) = h(2,2) + 2.0 * ( t + x(3) * t3 ) / r / r;
    h(2,3) = h(2,3) + 2.0 * ( t - x(3) * t3 * logr ) / r;

    h(3,1) = h(3,1) - 2.0 * t3 * logr;
    h(3,2) = h(3,2) + 2.0 * ( t - x(3) * t3 * logr ) / r;
    h(3,3) = h(3,3) + 2.0 * t3 * logr * logr;

  end

  h(1,1) = ( h(1,1) / x(1) ) / x(1);
  h(1,2) = h(1,2) * x(3) / x(1);
  h(1,3) = h(1,3) / x(1);

  h(2,1) = h(2,1) * x(3) / x(1);
  h(2,2) = h(2,2) * x(3);

  h(3,1) = h(3,1) / x(1);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function result = torus_6s2 ( func, r1, r2, result )

%*****************************************************************************80
%
%% TORUS_6S2 approximates an integral inside a torus in 3D.
%
%  Integration region:
%
%    Points (X,Y,Z) such that:
%
%      ( SQRT ( X**2 + Y**2 ) - R1 )**2 + Z**2 <= R2**2.
%
%  Discussion:
%
%    An 84 point 6-th degree formula is used, Stroud number TOR3-S2:6-1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971.
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied
%    function of three variables which is to be integrated,
%    of the form:
%      function value = func ( x, y, z )
%
%    Input, real R1, R2, the two radii that define the torus.
%
%    Output, real RESULT, the approximate integral of the function.
%
  norder = 2;
  s = [ 0.322914992E+00, 0.644171310E+00 ];
  weight = [ 0.387077796E+00, 0.165609800E+00 ];

  w = 1.0E+00 / ( 7.0E+00 * r1 * pi );

  quad = 0.0E+00;

  for n = 1 : 7

    u = 0.5E+00 * sqrt ( 3.0E+00 ) * r2;
    cth = cos ( 2.0E+00 * pi * n / 7.0E+00 );
    sth = sin ( 2.0E+00 * pi * n / 7.0E+00 );

    for i = 1 : 2

      u = -u;

      x = ( r1 + u ) * cth;
      y = ( r1 + u ) * sth;
      z = 0.0E+00;
      quad = quad + 0.232710567E+00 * w * ( r1 + u ) * feval ( func, x, y, z );

      x = r1 * cth;
      y = r1 * sth;
      z = u;
      quad = quad + 0.232710567E+00 * w * r1 * feval ( func, x, y, z );

    end

    for k = 1 : norder

      u = s(k) * r2;
      v = u;

      for i = 1 : 2

        u = -u;

        for j = 1 : 2

          v = -v;

          x = ( r1 + u ) * cth;
          y = ( r1 + u ) * sth;
          z = v;
          quad = quad + weight(k) * w * ( r1 + u ) * feval ( func, x, y, z );

        end
      end
    end
  end

  volume = torus_volume_3d ( r1, r2 );
  result = quad * volume;

  return
end


text:
function splabels = msLabelMap2Sp(lmap, smap)

nocell = false;
if ~iscell(lmap)
    nocell = true;
    lmap = {lmap};
    smap = {smap};
end

% if ~exist('maxlab', 'var') || isempty(maxlab)
%     maxlab = 0;
%     for f = 1:numel(lmap)
%         maxlab = max(maxlab, max(lmap{f}(:)));
%     end
% end

for f = 1:numel(lmap)

    nseg = max(smap{f}(:));
    stats = regionprops(smap{f}, 'PixelIdxList');
    idx = {stats.PixelIdxList};

    splabels{f} = zeros(nseg, 1);
    for s = 1:nseg
        splabels{f}(s) = mode(double(lmap{f}(idx{s})));
    end
end

if nocell
    splabels = splabels{1};
end

author: zouchuhang

repo: LayoutNet


text:
function T = RpToTrans(R, p)
% *** CHAPTER 3: RIGID-BODY MOTIONS ***
% Takes rotation matrix R and position p.
% Returns the corresponding homogeneous transformation matrix T in SE(3).
% Example Input:
%
% clear; clc;
% R = [[1, 0, 0]; [0, 0, -1]; [0, 1, 0]];
% p = [1; 2; 5];
% T = RpToTrans(R, p)
%
% Output:
% T =
%     1     0     0     1
%     0     0    -1     2
%     0     1     0     5
%     0     0     0     1

T = [R, p; 0, 0, 0, 1];
end

author: ShuoYangRobotics

repo: QuadrupedSim


text:
function [ clmParams, pdm_right, pdm_left ] = Load_CLM_params_eye()
%LOAD_CLM_PARAMS_WILD Summary of this function goes here
%   Detailed explanation goes here
    clmParams.window_size = [15,15; 13,13;];
    clmParams.numPatchIters = size(clmParams.window_size,1);

    % the PDM created from in the wild data
    pdmLoc = ['../models/hierarch_pdm/pdm_6_r_eye.mat'];

    load(pdmLoc);

    pdm_right = struct;
    pdm_right.M = double(M);
    pdm_right.E = double(E);
    pdm_right.V = double(V);

    pdmLoc = ['../models/hierarch_pdm/pdm_6_l_eye.mat'];

    load(pdmLoc);

    pdm_left = struct;
    pdm_left.M = double(M);
    pdm_left.E = double(E);
    pdm_left.V = double(V);

    % the default model parameters to use
    clmParams.regFactor = 0.1;
    clmParams.sigmaMeanShift = 2;
    clmParams.tikhonov_factor = 0;

    clmParams.startScale = 1;
    clmParams.num_RLMS_iter = 5;
    clmParams.fTol = 0.01;
    clmParams.useMultiScale = true;
    clmParams.use_multi_modal = 1;
    clmParams.tikhonov_factor = 0;
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function seed = pm_seed(angvar,mask,pxs)
% Find a suitable (hopefully) seed point from which
% to start watershed-based unwrapping.
% FORMAT: seed = pm_seed(angvar,mask,pxs)
%
% Input:
% angvar  : Map of variance of (voxelwise) estimates
%           of phase angle.
% mask    : Tells us which part of angvar to consider.
% pxs     : Array of voxel sizes, used to ensure
%           isotropic smoothing.
%
% Output:
% seed    : Coordinates of suitable seed point.
%
% In order to find a seed point we first threshold the
% variance map at a quarter of the variance of a U(-pi,pi)
% distribution. This gives us a binary image with ones only
% for low variance regions. This is then smoothed with a
% very wide gaussian kernel (50mm). The maximum of
% the smoothed map is then pretty much a centre-of-mass
% of the "low-variance volume". It could however in
% principle be a relatively high variance voxel
% surrounded by low-variance voxels. Therefore we pick
% a percentage of the highest voxels in the smooth map
% (i.e. we pick a neighbourhood) and then pick the location
% of those that has the lowest variance in the original
% variance map.
%___________________________________________________________
% Jesper Andersson 1/10-03
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Jesper Andersson
% $Id: pm_seed.m 4842 2012-08-15 18:02:30Z guillaume $

if nargin < 3
   mask = ones(size(angvar));
end

%
% First let us create a volume where a high
% value indicates a "high density of relatively
% low variance voxels".
%
dim = size(angvar);
if length(dim) == 2 dim(3) = 1; end

M = eye(4)*diag([pxs(1) pxs(2) pxs(3) 1]);
M = M - [zeros(4,3) M*[mean(1:dim(1)) mean(1:dim(2)) mean(1:dim(3)) 0]'];

P = struct('dim',     [dim 64],...
           'pinfo',   [1 0]',...
           'mat',     M);
P.dat = double(angvar<(pi^2)/12);
svol = zeros(size(angvar));
spm_smooth(P.dat,svol,50);

%
% A high value in svol "probably" indicates a
% voxel with a low variance, surrounded by a
% lot of other voxels with low variance (i.e.
% a good place to start unwrapping from).
% However, it COULD also be a voxel with
% "not so low variance" surrounded by low
% variance voxels. To avoid that trap we pick
% the voxel with the lowest variance in the
% unsmoothed variance map, out of the 5% of
% voxels in the svol with the highest values.
% It's all very heuristic and ugly.
%

[N,X] = hist(svol(logical(mask(:))),100);
indx = find(cumsum(N)>0.95*length(find(mask)));
thres = X(indx(1)-1);

indx = find(svol(:)>thres);
[mv,mi] = min(angvar(indx));
seed = zeros(1,3);
[seed(1),seed(2),seed(3)] = ind2sub(size(angvar),indx(mi));

return

author: spm


text:
tracker_label = 'Demo_Cpp';

% Example is stored here: <path-to-toolkit-dir>\tracker\examples\native
tracker_command = '<path-to-build-directory>\static_cpp_rgbd.exe';

% tracker_interpreter = []; % Set the interpreter used here as a lower case string. E.g. if you are using Matlab, write 'matlab'. (optional)

tracker_linkpath = {'<path-to-trax-build-dir>\bin', ...
    '<path-to-trax-build-dir>\lib'};

author: votchallenge

repo: toolkit-legacy