text:
% SPM5 UPDATE 14/11/06
% UPDATE 27/01/05
% Sets the default values for the FieldMap toolbox
%
% FORMAT pm_defaults_Sonata
%_______________________________________________________________________
%
% This file is intended for use with the Siemens fieldmap sequence
% on the Sonata scanner at the FIL with the being phased out sequence
% ralf_epi_tbr and nw_mepi_silent_v1a.
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Chloe Hutton and Jesper Andersson
% $Id: pm_defaults_Sonata.m 5015 2012-10-24 13:40:07Z guillaume $

global pm_def

% Defaults for creating field map. (See pm_make_fieldmap.m and
%                                   FieldMap.man for more info.)
%=======================================================================
pm_def.INPUT_DATA_FORMAT = 'PM';      % 'RI' = load two real and
                                      % imaginary image pairs
                                      % 'PM' = load one or two
                                      % phase and magnitude image
                                      % pairs.
pm_def.SHORT_ECHO_TIME = 10.0;        % Short echo time in ms for Sonata
pm_def.LONG_ECHO_TIME = 14.76;        % Long echo time in ms for Sonata
pm_def.MASKBRAIN = 1;                 % Do brain masking (1 or 0,
                      % 0 for EPI fieldmaps)

% Defaults for unwrapping options. (See pm_make_fieldmap.m and
%                                   FieldMap.man for more info.)
%=======================================================================
pm_def.UNWRAPPING_METHOD = 'Mark3D';  % Unwrapping options are:
                                      % 'Huttonish', 'Mark3D' or 'Mark2D'
pm_def.FWHM = 10;                     % FWHM of Gaussian filter used to
                                      % implement weighted smoothing of
                                      % unwrapped maps.
pm_def.PAD = 0;                       % Size of padding kernel if required.
pm_def.WS = 1;                        % Weighted or normal smoothing.

% Flags for brain extraction
%=======================================================================
pm_def.MFLAGS.TEMPLATE = fullfile(spm('Dir'),'toolbox','FieldMap','T1.nii');
pm_def.MFLAGS.FWHM = 5; % In mm
pm_def.MFLAGS.NERODE = 2;% In voxels
pm_def.MFLAGS.NDILATE = 4; % In voxels
pm_def.MFLAGS.THRESH = 0.5;
pm_def.MFLAGS.REG = 0.02; % A larger value helps segmentation to converge
pm_def.MFLAGS.GRAPHICS = 0; % A larger value helps segmentation to converge

% Defaults for converting field map to voxel displacement map.
%=======================================================================
pm_def.EPI_BASED_FIELDMAPS = 0;         % EPI=1, other=0.
pm_def.K_SPACE_TRAVERSAL_BLIP_DIR = -1; % +ve k-space = 1, -ve = -1.
pm_def.TOTAL_EPI_READOUT_TIME = 32;     % Sonata EPI RO time (500E-6*64)

% Defaults for Unwarping.
%=======================================================================
pm_def.DO_JACOBIAN_MODULATION = 0;    % Do jacobian modulation to adjust
                                      % for compression or stretching
                                      % No = 0, Yes = 1

author: spm


text:
classdef ChRotatorForFiberHomogenizer < handle

    methods (Access = public)

        function Ch = rotate(obj, dir,Ch)
            angle  = obj.computeRotationAngle(dir);
            dir    = obj.computeZdirection();
            C      = obj.makeChTensor(Ch);
            Ch = Rotator.rotate(C,angle,dir);
        end

    end

    methods (Access = private, Static)

         function angle = computeRotationAngle(dir)
            fiberDirection = dir.getValue;
            angle = -acos(dot(fiberDirection,[1 0 0]));
        end

        function C = makeChTensor(Ch)
            C = SymmetricFourthOrderPlaneStressVoigtTensor();
            C.setValue(Ch);
        end

        function dir = computeZdirection()
            d = [ 0 0 1];
            dir = Vector3D;
            dir.setValue(d);
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function y_line=runline(y,n,dn)
% Running line fit (local linear regression)
%
% Usage: y_line=runline(y,n,dn);
%
% Inputs:
% y: input 1-d time series (real)
% n: length of running window in samples
% dn: stepsize of window in samples
%
% Outputs:
% y_line: local line fit to data

y=y(:);
nt=length(y);
y_line=zeros(nt,1);
norm=y_line;
nwin=ceil((nt-n)/dn);
yfit=zeros(nwin,n);
xwt=((1:n)-n/2)/(n/2);
wt=(1-abs(xwt).^3).^3;
for j=1:nwin,
	tseg=y(dn*(j-1)+1:dn*(j-1)+n);
	y1=mean(tseg);
	y2=mean((1:n)'.*tseg)*2/(n+1);
	a=(y2-y1)*6/(n-1); b=y1-a*(n+1)/2;
	yfit(j,:)=(1:n)*a+b;
	y_line((j-1)*dn+(1:n))=y_line((j-1)*dn+(1:n))+(yfit(j,:).*wt)';
	norm((j-1)*dn+(1:n))=norm((j-1)*dn+(1:n))+wt';
end
mask=find(norm>0); y_line(mask)=y_line(mask)./norm(mask);
indx=(nwin-1)*dn+n-1;
npts=length(y)-indx+1;
y_line(indx:end)=(n+1:n+npts)'*a+b;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function f = eval(SO3F,rot,varargin)
% pointwise evaluation of SO3F at rotations rot
%
% sum_i \sum_j \sum_k  weights * psi(g*s_j h_i . s_k r_i)

% if isa(rot,'orientation')
%   ensureCompatibleSymmetries(SO3F,rot)
% end

s = size(rot);
rot = rot(:);

h = symmetrise(SO3F.h.normalize,'unqiue');

% rot x SS x CS x h
gh = reshape((SO3F.SS.properGroup * rot).' * h, length(rot),[],length(SO3F.h));

f = mean(SO3F.psi.eval(dot(gh,reshape(normalize(SO3F.r),1,1,[]),'noSymmetry')),2);

f = reshape(reshape(f,length(rot),[]) * SO3F.weights,s);

if isalmostreal(f)
  f = real(f);
end

author: mtex-toolbox

repo: mtex

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


text:
function triangle_plots ( )

%*****************************************************************************80
%
%% TRIANGLE_PLOTS plots the surface and the R(Theta) function for a triangle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 May 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGLE_PLOTS:\n' );
  fprintf ( 1, '  For a 2D triangle defined by a 0/1 characteristic function,\n' );
  fprintf ( 1, '  plot the surface, and R(Theta),\n' );
  fprintf ( 1, '  using a centered point, and then an offcentered point.\n' );
%
%  1): Plot the surface, using a centered base point.
%
  n = 61;
  theta = linspace ( 0.0, 2.0 * pi, n );

  f = zeros(n,1);
  x0 = [ 2.0 / 3.0, 1.0 / 3.0 ];
  for i = 1 : n
    r(i) = bisect_characteristic ( x0, theta(i), @triangle_characteristic );
  end
  x = x0(1) + r(1:n) .* cos ( theta(1:n) );
  y = x0(2) + r(1:n) .* sin ( theta(1:n) );

  xc = x0(1) + 0.5 .* cos ( theta(1:n) );
  yc = x0(2) + 0.5 .* sin ( theta(1:n) );

  clf
  hold on
  plot ( x0(1), x0(2), 'b.', 'MarkerSize', 50 );
  plot ( xc, yc, 'b-', 'LineWidth', 2 );
  plot ( x, y, 'ro', 'LineWidth', 3 );
  axis equal
  grid on
  xlabel ( '<---X--->', 'FontSize', 24 );
  ylabel ( '<---Y--->', 'FontSize', 24 );
  title ( 'Triangle transition surface', 'FontSize', 24)
  hold off
  filename =  'triangle_centered_surface.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Created plotfile "%s".\n', filename );
%
%  2): Plot R as a function of Theta for a centered base point.
%
  plot ( theta, r, 'LineWidth', 3 );
  xlabel ( '<---Theta--->', 'FontSize', 24 )
  ylabel ( '<---R(Theta)--->', 'FontSize', 24 )
  title ( 'R(Theta) from base point to surface', 'FontSize', 24 )
  grid on
  filename = 'triangle_centered_plot.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Created plotfile "%s".\n', filename );
%
%  3): Plot the surface, using an offcentered base point.
%
  n = 61;
  theta = linspace ( 0.0 , 2.0 * pi, n );

  f = zeros(n,1);
  x0 = [ 1.0, 1.50 ];
  for i = 1 : n
    r(i) = bisect_characteristic ( x0, theta(i), @triangle_characteristic );
  end

  x = x0(1) + r(1:n) .* cos ( theta(1:n) );
  y = x0(2) + r(1:n) .* sin ( theta(1:n) );
%
%  Coordinates of nominal triangle around (x0,y0), to
%  suggest coordinate system.
%
  xc = x0(1) + 0.5 .* cos ( theta(1:n) );
  yc = x0(2) + 0.5 .* sin ( theta(1:n) );

  clf
  hold on
  plot ( x0(1), x0(2), 'b.', 'MarkerSize', 50 );
  plot ( xc, yc, 'b-', 'LineWidth', 2 );
  plot ( x, y, 'ro', 'LineWidth', 3 );
  axis equal
  grid on
  xlabel ( '<---X--->', 'FontSize', 24 );
  ylabel ( '<---Y--->', 'FontSize', 24 );
  title ( 'Triangle transition surface', 'FontSize', 24 )
  hold off
  filename = 'triangle_offcentered_surface.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Created plotfile "%s".\n', filename );
%
%  4): Plot R as a function of Theta for an offcentered base point.
%
  plot ( theta, r, 'LineWidth', 3 );
  xlabel ( '<---Theta--->', 'FontSize', 24 )
  ylabel ( '<---R(Theta)--->', 'FontSize', 24 )
  title ( 'R(Theta) from base point to surface', 'FontSize', 24 )
  grid on
  filename = 'triangle_offcentered_plot.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Created plotfile "%s".\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function val=intCotPow(u,n)
%%INTCOTPOWER Evaluate the integral of cot(u)^n du.  A definite integral
%             can be evaluated, or an indefinite integral (with a
%             particular additive constant).
%
%INPUTS: u A 2XN (for definite integral) or a 1XN (for indefinite
%          integrals) set of N points. For definite integrals, u(1,:) are
%          the real lower bounds and u(2,:) are the real upper bounds. For
%          indefinite integrals, the integral is evaluated at the points in
%          u. The values in u should be real.
%        n The positive integer exponent of the cotangent term.
%
%OUTPUTS: val The 1XN set of values of the integral of cot(u)^n.
%
%This function simply implements the recursion that arises from integration
%by parts from basic calculus, as are given in the tables in the back of
%[1].
%
%REFERENCES:
%[1] J. Stewart, Calculus, 7th ed. Belmont, CA: Brooks/Cole, 2012.
%
%October 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numDim=size(u,1);

if(isempty(u))
   val=[];
   return;
end

if(numDim==1)%An indefinite integral
    val=indefIntCotPow(u,n);
else%A definite integral
    val=indefIntCotPow(u(2,:),n)-indefIntCotPow(u(1,:),n);
end

end

function val=indefIntCotPow(u,n)

if(n==0)
    val=u;
    return;
elseif(n==1)
    val=log(abs(sin(u)));
    return;
end

%If here, n>=1
cotVal=cot(u);

val=0;
if(mod(n,2)==0)%If n is even
    endVal=2;
else
    endVal=3;
end

coeffProd=1;
for k=n:-2:endVal
    val=val-coeffProd*(1/(k-1))*cotVal.^(k-1);
    coeffProd=-coeffProd;
end

if(endVal==2)
    %The final integral is over cot^0
    val=val+coeffProd*u;
else
    %The final integral is over cot^1
    val=val+coeffProd*log(abs(sin(u)));
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function check_NFSOFT
%
% checks for Y_l(g h) = T_l(g) Y_l(h)
%

%% input data
L = 16;

%q = axis2quat(xvector+yvector,25*degree);
%q = axis2quat(zvector,90*degree);
%h = [xvector,-xvector,yvector];

qq = quaternion(SO3Grid(10));
h = equispacedS2Grid('points',20,'antipodal');

progress(0,length(qq));

for iq = 1:length(qq)

  progress(iq,length(qq));

  q = qq(iq);

  %% convert to export parameters
  g = Euler(q,'nfft');
%   alpha = fft_rho(alpha); %-->z
%   beta  = fft_theta(beta);
%   gamma = fft_rho(gamma); %-->z
%   g = 2*pi*[alpha;beta;gamma];
  %g = [0;pi/2;0];

  %% set parameters
  c = 1;
  A = ones(1,L+1);

  %% run NFSOFT
  T = call_extern('odf2fc','EXTERN',g,c,A); % conjugate(D)

  % extract result
  T = complex(T(1:2:end),T(2:2:end));

  %% check result

  for l = 0:L

    Y = sphericalY(l,h).'; % -> Y
    gY = sphericalY(l,q*h).'; % -> Y
    Tl = reshape(T(deg2dim(l)+1:deg2dim(l+1)),2*l+1,2*l+1);
    %    TY = flipud(Tl * flipud(Y));
    %    TY = flipud(fliplr(Tl) * Y);
    TY = Tl * Y;
    er(l+1,iq) = sqrt(norm(TY(:) - gY(:)));
    %TY = conj(Tl) * Y;norm(TY(:) - gY(:))
    %TY = Tl' * Y;norm(TY(:) - gY(:))
    %TY = Tl.' * Y;norm(TY(:) - gY(:))

  end

end

%plot(mean(er,2));
if mean(er(:)) > 0.001
  error('Error in NFSOFT');
else
  disp('checking NFSOFT: ok')
  disp(mean(er(:)))
end

pcolor(er)

function d = deg2dim(l)
% dimension of the harmonic space up to order l

d = l*(2*l-1)*(2*l+1)/3;


text:
function  bdiff_bdepth(mycat)
    %  This routine etsimates the b-value of a curve automatically
    %  The b-valkue curve is differenciated and the point
    %  of maximum curvature marked. The b-value will be calculated
    %  using this point and the point half way toward the high
    %  magnitude end of the b-value curve.

    %  originally, "mycat" was "newcat"
    %  Stefan Wiemer 1/95
    %
    global cluscat mess bfig backcat magsteps_desc bvalsum3  bval aw bw t1 t2 t3 t4 dloop leg1 leg2
    global teb t0b cua  ew onesigma mrt bvalsumhold
    global gBdiff % contains b1, n1, b2, n2
    global mxlkbt lsbt ni
    ZG=ZmapGlobal.Data;

    report_this_filefun();

    bfig=findobj('Type','Figure','-and','Name','frequency-magnitude distribution');
    if ~isempty(bfig)
        if dloop == 2
            ZG.hold_state=true;
        end
    else
        bfig=figure_w_normalized_uicontrolunits(...                  %build figure for plot
            'Units','normalized','NumberTitle','off',...
            'Name','frequency-magnitude distribution',...
            'visible','off',...
            'pos',[ 0.300  0.3 0.4 0.6]);
        ZG.hold_state=false;

        uicontrol('Units','normal',...
            'Position',[.0 .85 .08 .06],'String','Info ',...
            'callback',@(~,~)infoz(1));
        uicontrol('Units','normal',...
            'Position',[.0 .45 .10 .06],'String','Manual ',...
            'callback',@(~,~)bfitnew(mycat));

        uicontrol('Units','normal',...
            'Position',[.0 .35 .10 .06],'String','RecTime ',...
            'callback',@callbackfun_003);

        uicontrol('Units','normal',...
            'Position',[.0 .25 .10 .06],'String','TimePlot ',...
            'callback',@cb_timeplot);

        uicontrol('Units','normal',...
            'Position',[.0 .65 .08 .06],'String','Save ',...
            'Callback',{@calSave9,magsteps_desc, bvalsum3});
    end

    maxmag = ceil(10*max(mycat.Magnitude))/10;
    mima = min(mycat.Magnitude);
    mima = min(mima, 0);

    % number of mag units
    nmagu = (maxmag*10)+1;

    [bval,xt2] = hist(mycat.Magnitude,(mima:0.1:maxmag));
    bvalsum = cumsum(bval); % N for M <=
    bval2 = bval(end:-1:1);
    bvalsum3 = cumsum(bval(end:-1:1));    % N for M >= (counted backwards)
    magsteps_desc = (maxmag:-0.1:mima);

    backg_ab = log10(bvalsum3);
    orient tall

    if ZG.hold_state
        axes(cua)
        disp("set(gca,'NextPlot','add')")
        set(gca,'NextPlot','add')
    else
        figure(bfig);delete(findobj(bfig,'Type','axes'));
        rect = [0.2,  0.3, 0.70, 0.6];           % plot Freq-Mag curves
        axes('position',rect);
    end

    pldepth =semilogy(magsteps_desc,bvalsum3,'sb');
    set(pldepth,'LineWidth',1.0,'MarkerSize',6,...
        'MarkerFaceColor','r','MarkerEdgeColor','b');
    set(gca,'NextPlot','add')
    difb = [0 diff(bvalsum3) ];

    % Marks the point of maximum curvature
    %
    i = find(difb == max(difb));
    i = max(i);
    te = semilogy(magsteps_desc(i),bvalsum3(i),'xk');
    set(te,'LineWidth',1.5,'MarkerSize',ms10)

    % Estimate the b-value
    %
    i2 = 1 ;
    te = semilogy(magsteps_desc(i2),difb(i2),'xk');
    set(te,'LineWidth',1.5,'MarkerSize',ms10)
    te = semilogy(magsteps_desc(i2),bvalsum3(i2),'xk');
    set(te,'LineWidth',1.5,'MarkerSize',ms10)

    xlabel('Magnitude','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s)
    ylabel('Cumulative Number','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s)
    %set(gca,'Color',color_bg)
    set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','normal',...
        'FontWeight','normal','LineWidth',1.0,...
        'Box','on','Tag','cufi')

    cua = gca;

    M1b = [];
    M1b = [magsteps_desc(i) bvalsum3(i)];
    tt3=num2str(fix(100*M1b(1))/100);

    M2b = [];
    M2b =  [magsteps_desc(i2) bvalsum3(i2)];
    tt4=num2str(fix(100*M2b(1))/100);

    ll = magsteps_desc >= M1b(1)-0.05  & magsteps_desc <= M2b(1) +0.05;
    x = magsteps_desc(ll);

    l2 = mycat.Magnitude >= M1b(1)- 0.05  & mycat.Magnitude <= M2b(1)+ 0.05;
    [ bv, onesigma, av] = calc_bmemag(mycat.Magnitude(l2), 0.1) ;

    bv = -bv;

    pause(0.1)

    y = backg_ab(ll);
    [aw bw,  S, ew] = wls(x',y');
    p = [bw aw];
    %[p,S] = polyfit(x,y,1)                    % fit a line to background
    p2 = [bw+onesigma aw];
    p3 = [bw-onesigma aw];
    x2 = 1:0.1:6;
    f = polyval(p,x);
    f2 = polyval(p2,x);
    f3 = polyval(p3,x);
    [f4,delta] = polyval(p,x,S);
    Tr  = (teb-t0b)/(10.^ polyval(p,mrt));
    fprintf('Recurrence time Tr(M%g) = %g years\n', mrt, Tr);
    f = 10.^f;
    f2 = 10.^f2;
    f3 = 10.^f3;
    f4 = 10.^f4;
    delta = 10.^delta;
    set(gca,'NextPlot','add')
    ttm= semilogy(x,f,'r');                         % plot linear fit to backg
    set(ttm,'LineWidth',1)
    ttm= semilogy(x,f2,'k');                         % plot linear fit to backg
    set(ttm,'LineWidth',1)
    ttm= semilogy(x,f3,'k');                         % plot linear fit to backg
    set(ttm,'LineWidth',1)
    ttm= semilogy(x,f4-delta,'k-.');                         % plot linear fit to backg
    set(ttm,'LineWidth',1)
    ttm= semilogy(x,f4+delta,'k-.');
    set(ttm,'LineWidth',1)
    set(gca,'XLim',[min(mycat.Magnitude)-0.5  max(mycat.Magnitude)+0.5])
    set(gca,'YLim',[1 (mycat.Count+20)*1.4]);

    r = corrcoef(x,y);
    r = r(1,2);
    %std_backg = std(y - polyval(p,x));      % standard deviation of fit
    std_backg = ew;      % standard deviation of fit

    p=-p(1,1);
    p=fix(100*p)/100;
    std_backg=fix(100*std_backg)/100;
    tt1=num2str(bw,3);
    tt2=num2str(std_backg);
    tt4=num2str(bv,3);
    tt5=num2str(onesigma,2);

    rect=[0 0 1 1];
    h2=axes('position',rect);
    set(h2,'visible','off');

    %if ZG.hold_state
    ge_symbol = char(8805);
    if dloop == 2
        set(pldepth,'LineWidth',1.0,'MarkerSize',6,...
            'MarkerFaceColor','y','MarkerEdgeColor','g','Marker','o');
        set(cua,'Ylim',[ 1   ni ] );

        txt1=text(.10, .08,['Bottom Zone b-value (w LS, M  ', ge_symbol, ' ', num2str(M1b(1)) '): ',tt1, ' +/- ', tt2 ',a-value = ' , num2str(aw) ]);
        set(txt1,'FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s,'Color','r')
        txt1=text(.10, .04,['Bottom Zone b-value (max lik, M ', ge_symbol, ' ', num2str(min(mycat.Magnitude)) '): ',tt4, ' +/- ', tt5,',a-value = ' , num2str(av)]);
        set(txt1,'FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s, 'Color', 'r')
        lsbb = bw; mxlkbb = bv;

        lsb = lsbt/lsbb; mxlkb = mxlkbt/mxlkbb;
        slsb = num2str(lsb);
        smxlkb = num2str(mxlkb);
        txt3 = text(.25, .94,['LS b ratio = ', slsb,'           max lik b ratio = ',smxlkb]);
        set(txt3,'FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s,'Color','b')
        %  leg(2)=pldepth
    else
        txt1=text(.10, .18,['Top Zone b-value (w LS, M  ', ge_symbol, ' ', num2str(M1b(1)) '): ',tt1, ' +/- ', tt2, ',a-value = ' , num2str(aw) ]);
        set(txt1,'FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s)
        txt1=text(.10, .14,['Top Zone b-value (max lik, M ', ge_symbol, ' ', num2str(min(mycat.Magnitude)) '): ',tt4, ' +/- ', tt5,',a-value = ' , num2str(av)]);
        set(txt1,'FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s)
        set(gcf,'PaperPosition',[0.5 0.5 4.0 5.5])
        lsbt = bw; mxlkbt = bv;
    end

    if dloop == 1
        leg1 = pldepth;
    end
    if dloop == 2
        leg2 = pldepth;
        legend([leg1,leg2],'Top depth zone','Bottom depth zone');
    end
    set(gcf,'visible','on');

    if ZG.hold_state
        % calculate the probability that the two distributions are different
        l = mycat.Magnitude >=  M1b(1);
        gBdiff.b2 = str2double(tt1); gBdiff.n2 = M1b(2);
        n = gBdiff.n1+gBdiff.n2;
        da = -2*n*log(n) + 2*gBdiff.n1*log(gBdiff.n1+gBdiff.n2*gBdiff.b1/gBdiff.b2) + 2*gBdiff.n2*log(gBdiff.n1*gBdiff.b2/gBdiff.b1+gBdiff.n2) -2;
        pr = exp(-da/2-2);
        fprintf('Probability: %g\n', pr);
        txt1=text(.60, .75,['p=  ', num2str(pr,2)],'Units','normalized');
        set(txt1,'FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s)
        txt1=text(.60, .70, sprintf('gBdiff.n1: %g, gBdiff.n2: %g, gBdiff.b1: %g, gBdiff.b2: %g', ...
            gBdiff.n1, gBdiff.n2, gBdiff.b1, gBdiff.b2));
        set(txt1,'FontSize',8,'Units','normalized')
    else
        gBdiff.b1 = str2double(tt1); gBdiff.n1 = M1b(2);
    end

    bvalsumhold = bvalsum3;
    da = 10^(aw+bw*6.5);
    db = 10^(aw+bw*6.5)*(-6.5);
    dp = sqrt(da^2*ew^2+db^2*0.05^2);
    dr = 1/dp;

    set(gca,'NextPlot','replace');

    function callbackfun_003(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        est_recurrence_time_prob(mhcat, onesigma, aw, bw);
    end

    function cb_timeplot(mysrc,myevt)
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ZG=ZmapGlobal.Data;
        ZG.newt2 = mycat;
        ctp=CumTimePlot(mycat);
        ctp.plot();
    end
end

author: CelsoReyes


text:
function [b0,dt6] = dtiInterpDt6(mmPerVoxOut, dt6PathName, outPathName)
%
% [b0,dt6] = dtiInterpDt6(mmPerVoxOut, b0PathName, outPathName)
%
% Takes dt6 file and upsamples b0 and dt6 data to desired resolution
% Uses C versions of pajevic and interp3 algorithms.
% The input file should be already preped using dtiMakeDt6 to
% some resolution
%
% HISTORY:
% 2004.02.27 GSM (gmulye@white.stanford.edu) wrote it (based on
% dtiMakeDt6)

if ~exist('mmPerVoxOut','var') | isempty(mmPerVoxOut)
    mmPerVoxOut = [1 1 1];
end

if ~exist('dt6PathName','var') | isempty(dt6PathName)
   [f, p] = uigetfile({'*.mat','dt6 files (*.mat)';'*.*','All files'}, 'Select one of the dt6 files...');
   if(isnumeric(f)) error('Need a dt6 file to continue...'); end
   dt6PathName = fullfile(p, f);
   disp(dt6PathName);
end

% Upsampling tensors
disp('Interpolating tensors...');
oldDt6File = load(dt6PathName);
dt6 = oldDt6File.dt6;
xform = oldDt6File.xformToAcPc;
mmPerVoxIn = oldDt6File.mmPerVox;
sz = size(dt6);
boundingBox = [0 0 0; sz(1:3).*mmPerVoxIn-1]; %with this config, A(i) = A(2i-2)
%===================================
%Ripped from dtiResliceTensorAffine
%===================================
%We leave out all the transforms in dtiTensorAffine
%Tensors do not have to be rotated or brought to image space
x = (boundingBox(1,1):mmPerVoxOut(1):boundingBox(2,1));
y = (boundingBox(1,2):mmPerVoxOut(2):boundingBox(2,2));
z = (boundingBox(1,3):mmPerVoxOut(3):boundingBox(2,3));

[X,Y,Z] = meshgrid(x, y, z);
clear x y z;
coords = [X(:) Y(:) Z(:)]; %coords is Nx3
newSize = size(X);
clear X Y Z;
% HACK TO FIND CORRECT PAJEVIC LIBRARY
old = pwd;
cd(fileparts(which('dtiTensorInterp_Pajevic')));
dt6 = dtiTensorInterp_Pajevic(dt6, coords, mmPerVoxIn, 1, mmPerVoxIn./2);
cd(old);
dt6 = reshape(dt6, [newSize, 6]);

% Upsampling B0 grid
disp('Interpolating B0...');
% mmPerVox = b0.cannonical_mmPerVox;
% notes = b0.notes;
% xformToAnat = b0.anatXform;
% xformToTal = b0.talXform;
% xformToAcPc = b0.acpcXform;
% b0 = int16(round(b0.cannonical_img));
% anat.img = int16(t1.cannonical_img);
% anat.mmPerVox = t1.cannonical_mmPerVox;
% anat.xformToTal = t1.talXform;
% anat.xformToAcPc = t1.acpcXform;
% save(fullfile(outPathName,'dt6'), 'dt6', 'mmPerVox', 'notes', 'xformToAnat', 'xformToTal', 'xformToAcPc', 'b0', 'anat');

%============================
%Ripped from dtiResliceVolume
%============================
% myCinterp3 likes [rows,columns,slices], so we need a permute here
% myCinterp3 also wants pure image space, so we remove the mmPerVox.
coords(:,1) = coords(:,1)./mmPerVoxIn(1);
coords(:,2) = coords(:,2)./mmPerVoxIn(2);
coords(:,3) = coords(:,3)./mmPerVoxIn(3);
b0 = oldDt6File.b0;
b0 = double(permute(b0,[2,1,3]));
interpedb0 = myCinterp3(b0, [size(b0,1) size(b0,2)], size(b0,3), coords, 0.0);
b0 = squeeze(reshape(interpedb0,newSize));
clear interpedb0;

if(~exist('outPathName','var') | isempty(outPathName))
    [fullpath fName] = fileparts(dt6PathName);
    outPathName = fullfile(fullpath, [fName sprintf('%dx%dx%dmm',mmPerVoxIn) '_upto_' sprintf('%dx%dx%dmm',mmPerVoxOut)]);
end

%Find new transform to anatomy image
oldXformToAnat = oldDt6File.xformToAnat;
anat = oldDt6File.anat;
scale = mmPerVoxOut./anat.mmPerVox;
scale = diag([scale 1]);
xformToAnat = scale*oldXformToAnat;

mmPerVox = mmPerVoxOut;
notes = oldDt6File.notes;
b0 = int16(round(b0));
b0 = permute(b0,[2,1,3]); %Undo earlier permute
dt6 = permute(dt6,[2 1 3 4]); %Undo permute from earlier

disp(['Saving to ' outPathName '...']);
save(outPathName, 'b0', 'xformToAnat','anat','notes','mmPerVox','dt6');
save([outPathName '_anon'], 'b0', 'xformToAnat','anat','notes','mmPerVox','dt6');

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Upsampling procedure.
%
% Argments:
%   'I': greyscale image
%   'odd': 2-vector of binary values, indicates whether the upsampled image
%   should have odd size for the respective dimensions
%   'filter': upsampling filter
%
% If image width W is odd, then the resulting image will have width (W-1)/2+1,
% Same for height.
%
% tom.mertens@gmail.com, August 2007
%

function R = upsample(I,odd,filter)

% increase resolution
I = padarray(I,[1 1 0],'replicate'); % pad the image with a 1-pixel border
r = 2*size(I,1);
c = 2*size(I,2);
k = size(I,3);
R = zeros(r,c,k);
R(1:2:r, 1:2:c, :) = 4*I; % increase size 2 times; the padding is now 2 pixels wide

% interpolate, convolve with separable filter
R = imfilter(R,filter);     %horizontal
R = imfilter(R,filter');    %vertical

% remove the border
R = R(3:r - 2 - odd(1), 3:c - 2 - odd(2), :);

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function       blk_arr = Block_Matching( X, par)
% record the indexs of patches similar to the seed patch
blk_arr   =  zeros(par.nlsp, par.lenrc, 'single');
for  i  =  1 : par.lenrc
    seed = X(:, par.SelfIndex(i));
    neighbor = X(:, par.NeighborIndex(1:par.NumIndex(i), i));
    dis = sum(bsxfun(@minus, neighbor, seed).^2, 1);
    [~,ind]   =  sort(dis);
    indc        =  par.NeighborIndex( ind( 1:par.nlsp ), i );
%     indc(indc == par.SelfIndex(i)) = indc(1); % added on 08/01/2017
%     indc(1) = par.SelfIndex(i); % to make sure the first one of indc equals to off
    blk_arr(:, i) = indc;
end

author: csjunxu


text:
function [structOut, errors] = checkPrepDefaults(structIn, structOut, defaults)
% Check structIn input parameters against defaults fields of structOut
%
% Parameters:
%    structIn     Parameter structure to be checked
%    structOut    Base output structure for parameters.
%    defaults     Default information
%
% Note:  structOut contains fields that will always be there.
errors = cell(0);
fNames = fieldnames(defaults);
for k = 1:length(fNames)
    try
       nextValue = getStructureParameters(structIn, fNames{k}, ...
                                          defaults.(fNames{k}).value);
       validateattributes(nextValue, defaults.(fNames{k}).classes, ...
                         defaults.(fNames{k}).attributes);
       structOut.(fNames{k}) = nextValue;
    catch mex
        errors{end+1} = [fNames{k} ' parameter value invalid: ' mex.message]; %#ok<AGROW>
    end
end

author: VisLab

repo: EEG-Clean-Tools


text:
function [psiMapping,phiMapping,lies] = getMappings(Problem,indv,archive,option)
% Construct Psi-mapping or Phi-mapping or both

%------------------------------- 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".
%--------------------------------------------------------------------------

    psiApprox    = [];
    phiApprox    = [];
    constructPsi = 0;
    constructPhi = 0;

    if nargin < 4
        option = 3;
    end
    if option == 1
        constructPsi = 1;
    elseif option == 2
        constructPhi = 1;
    elseif option == 3
        constructPsi = 1;
        constructPhi = 1;
    end

    % Minimum number of Lower level points required for quadratic approximation of both mappings
    numberLowerLevelPointsTrain = (Problem.DU+1)*(Problem.DU+2)/2+2*(Problem.DU);
    numberLowerLevelPointsEval  = Problem.DU;
    minNumberLowerLevelPoints   = numberLowerLevelPointsTrain + numberLowerLevelPointsEval;

    PopDec = archive.decs;
    quadArchive.upper = PopDec(:,1:Problem.DU);
    quadArchive.lower = PopDec(:,Problem.DU+1:end);

    % Only the members close to the offspring are used for quadratic approximation
    distances = zeros(size(quadArchive.upper,1));
    for k = 1 : length(quadArchive)
        distances(k) = sum((indv - quadArchive.upper(k,:)).^2);
    end
    [~,I]  = sort(distances);
    I      = I(1:minNumberLowerLevelPoints);
    sizeI  = length(I);
    permut = randperm(sizeI);
    J      = permut(1:sizeI-numberLowerLevelPointsEval);
    quadApproxMembers = I(J);
    % I stores the members being considered from the archive
    % quadApproxMembers stores the members from I used for a quadratic
    setDiff = setdiff(I,quadApproxMembers);
    if constructPsi == 1
        psiApprox = cell(1,Problem.DL); sumMSE = 0;
        for j = 1 : Problem.DL
            psiApprox{j} = QuadApprox(quadArchive.lower(quadApproxMembers,j),quadArchive.upper(quadApproxMembers,:));
            %Calculating sum of MSE for all lower level variable approximations
            sumMSE = sumMSE + psiApprox{j}.mseNorm;
        end
        predictedLowerLevelVariables = zeros(length(setDiff),Problem.DL);
        for k = 1 : length(setDiff)
            for j = 1 : Problem.DL
                predictedLowerLevelVariables(k,j) = psiApprox{j}.constant + quadArchive.upper(setDiff(k),:)*psiApprox{j}.linear + quadArchive.upper(setDiff(k),:)*psiApprox{j}.sqmatrix*quadArchive.upper(setDiff(k),:)';
            end
        end
        psiMapping.function = psiApprox;
        psiMapping.sumMSE   = sumMSE;
        psiMapping.validMSE = mean(mean((predictedLowerLevelVariables - quadArchive.lower(setDiff,:)).^2));
    end

    if constructPhi == 1
        PopObj        = archive.objs;
        quadArchive.llFunctionValue = PopObj(:,2);
        llFunctionDim = size(quadArchive.llFunctionValue,2);
        sumMSE        = 0;
        phiApprox     = cell(1,llFunctionDim);
        for j = 1 : llFunctionDim
            phiApprox{j} = QuadApprox(quadArchive.llFunctionValue(quadApproxMembers,j), quadArchive.upper(quadApproxMembers,:));
            sumMSE       = sumMSE+phiApprox{j}.mseNorm;
        end
        predictedLowerLevelObjective = zeros(length(setDiff),llFunctionDim);
        for k = 1 : length(setDiff)
            for j = 1 : llFunctionDim
                predictedLowerLevelObjective(k,j) = phiApprox{j}.constant + quadArchive.upper(setDiff(k),:)*phiApprox{j}.linear + quadArchive.upper(setDiff(k),:)*phiApprox{j}.sqmatrix*quadArchive.upper(setDiff(k),:)';
            end
        end
        phiMapping.function = phiApprox;
        phiMapping.sumMSE   = sumMSE;
        phiMapping.validMSE = mean(mean((predictedLowerLevelObjective - quadArchive.llFunctionValue(setDiff,:)).^2));
    end

    %Checks if the offspring lies in between the training points and not outside
    lies = 1;
    for j = 1 : size(quadArchive.upper,2)
        if indv(j)>=max(quadArchive.upper(quadApproxMembers,j)) || indv(j)<=min(quadArchive.upper(quadApproxMembers,j))
            lies = 0;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% Pixelwise magnitude of a MxNx2 matrix
%
%    Copyright (C) 2013  Anestis Papazoglou
%
%    You can redistribute and/or modify this software for non-commercial use
%    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/>.
%
%    For commercial use, contact the author for licensing options.
%
%    Contact: a.papazoglou@sms.ed.ac.uk

function result = getMagnitude( input )

    if( ~isfloat( input ) )
        input = single( input );
    end
    result = sqrt( input( :, :, 1 ).^2 + input( :, :, 2 ).^2 );

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% theta = choose_orientations2(x,y,phi,wtheta,wphi)
%
% we will set the orientation to theta_t = phi_t + s_t * pi
% we want to choose s_t to minimize
% \sum_t cost(s_t|s_{t-1})
% cost(s_t|s_{t-1}) = [wtheta_t*d(\theta_t,\theta_{t-1}) +
%                      wphi(||v_t||^2)*d(\theta_t,angle(v_t))]
%
% we will find the most likely states s_t using the recursion
% cost_t(s_t) = min_{s_{t-1}} { cost_{t-1}(s_{t-1}) + cost(s_t|s_{t-1})
%
% Inputs:
% x: N x 1 vector where x(t) is the x-coordinate of the center of the fly
% at time t
% y: N x 1 vector where y(t) is the y-coordinate of the center of the fly
% at time t
% phi: N x 1 vector where phi(t) is the orientation of the fly at time t
% wtheta: N x 1 vector where wtheta(t) is the weight of the change in
% orientation term at time t
% wphi: N x 1 vector where wphi(t) is the weight of the velocity
% direction term at time t
%
function theta = choose_orientations2(x,y,theta,weight_theta,weight_phi)

inputsz = size(x);
x = x(:);
y = y(:);
theta = theta(:);
weight_phi = weight_phi(:);
weight_theta = weight_theta(:);

% number of frames
N = length(x);

% allocate space for storing the optimal path
stateprev = zeros(N-1,2);

% allocate space for computing costs
tmpcost = zeros(2,1);
costprevnew = zeros(2,1);

% initialize first frame
costprev = zeros(2,1);

% compute velocity
vx = [0;diff(x)];
vy = [0;diff(y)];

% compute angle of velocity
velocityangle = atan2(vy,vx);

% compute iteratively
for t = 2:N,

  % compute for both possible states
  for scurr = 1:2,

    % try both previous states
    thetacurr = theta(t) + (scurr-1)*pi;

    for sprev = 1:2,

      thetaprev = theta(t-1) + (sprev-1)*pi;
      costcurr = weight_theta(t)*angledist(thetaprev,thetacurr) + ...
        weight_phi(t)*angledist(thetacurr,velocityangle(t));
      tmpcost(sprev) = costprev(sprev) + costcurr;

    end

    % choose the minimum
    sprev = argmin(tmpcost);

    % set pointer for path
    stateprev(t-1,scurr) = sprev;

    % set cost
    costprevnew(scurr) = tmpcost(sprev);

  end

  % copy over
  costprev(:) = costprevnew(:);

end

% choose the best last state
scurr = argmin(costprev);

if scurr == 2,
  theta(end) = modrange(theta(end)+pi,-pi,pi);
end

% choose the best states
for t = N-1:-1:1,
  scurr = stateprev(t,scurr);
  if scurr == 2,
    theta(t) = modrange(theta(t)+pi,-pi,pi);
  end
end

theta = reshape(theta,inputsz);

author: kristinbranson

repo: JAABA


text:
function [out] = score_ridge(mdv, hilo, lambda, crossval)
% Calculates ridge score
%
% USAGE:
%
%    [out] = score_ridge(mdv, hilo, lambda, crossval)
%
% INPUTS:
%    mdv:         structure
%    hilo:        (0's and 1's), ideally there will be a similar # of each.
%
% OPTIONAL INPUTS:
%    lambda:      ridge parameter, default = .01
%    crossval:    whether to do cross validation. This severely slows
%                 down the computation, default is no.
%
% OUTPUT:
%    out:         score

if nargin < 4
    crossval = 0;  % don't cross validate by default
end
if nargin < 3
    lambda = .01;
end

hilo(hilo ==0) = -1; % make 0's into ones.
[nvars, npoints] = size(mdv);

if crossval == 0
    ybar = mean(hilo);
    b = inv(mdv*mdv' + lambda*eye(nvars))*mdv*hilo;

    yhat = mdv'*b+ybar;
    out = sum(sign(yhat) == sign(hilo))/npoints;
else
    results = zeros(npoints,1);
    for i = 1:npoints
        hilot = hilo([1:i-1,i+1:npoints]);
        mdvt = mdv(:,[1:i-1,i+1:npoints]);

        ybart = mean(hilot);
        b = inv(mdvt*mdvt' + lambda*eye(nvars))*mdvt*hilot;

        yhat = mdv(:,i)'*b+ybart;
        results(i)  = sign(yhat) == sign(hilo(i));
    end
    out = sum(results)/npoints;
end

author: opencobra

repo: cobratoolbox


text:
subset_map = containers.Map;
[a,b,c] = fileparts(image_list{1});
last_folder_all = a;
folder_names = strsplit(a,'\');
last_folder = [folder_names{end-1} '/' folder_names{end}];
subset_map(last_folder) = [1 0];
for i=2:length(image_list)
    [a,b,c] = fileparts(image_list{i});
    folder_names = strsplit(a,'\');
    if strcmp(last_folder_all, a) == 0
        start_end = subset_map(last_folder);
        subset_map(last_folder) = [start_end(1) i-1];
        last_folder = [folder_names{end-1} '/' folder_names{end}];
        subset_map(last_folder) = [i 0];
    end;
    last_folder_all = a;
end;
start_end = subset_map(last_folder);
subset_map(last_folder) = [start_end(1) length(image_list)];

author: happynear

repo: FaceVerification


text:
function [linGrid,gridSize,xx,yy,zz] = prtPlotUtilGenerateGrid(plotMins, plotMaxs, nSamplesPerDim, includePoints)
% [linGrid,gridSize,xx,yy,zz] = prtPlotUtilGenerateGrid(plotMins, plotMaxs, nSamplesPerDim, includePoints)
% Internal function,
% xxx Need Help xxx

yy = [];
zz = [];

nDims = length(plotMins);

if ~isnumeric(nSamplesPerDim)
    try
        nSamplesPerDim = nSamplesPerDim.nSamplesPerDim;
    catch  %#ok<CTCH>
        error('prt:prtPlotUtilGenerateGrid','Invalid nSamplesPerDim input');
    end
end

xx = linspace(plotMins(1),plotMaxs(1),nSamplesPerDim(nDims));
if nargin > 3 && ~isempty(includePoints)
    xx = cat(2,xx,includePoints(:,1)');
    xx = sort(xx,'ascend');
end

if nDims > 1
    yy = linspace(plotMins(2),plotMaxs(2),nSamplesPerDim(nDims));
    if nargin > 3 && ~isempty(includePoints)
        yy = cat(2,yy,includePoints(:,2)');
        yy = sort(yy,'ascend');
    end
end

if nDims > 2
    zz = linspace(plotMins(3),plotMaxs(3),nSamplesPerDim(nDims));
    if nargin > 3 && ~isempty(includePoints)
        zz = cat(2,zz,includePoints(:,3)');
        zz = sort(zz,'ascend');
    end
end

switch nDims
     case 1
        linGrid = xx(:);
        gridSize = size(xx);
        yy = ones(size(xx));
    case 2
        [xx,yy] = meshgrid(xx,yy);
        linGrid = [xx(:),yy(:)];
        gridSize = size(xx);
    case 3
        [xx,yy,zz] = meshgrid(xx,yy,zz);
        gridSize = size(xx);
        linGrid = [xx(:), yy(:), zz(:)];
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function fem1d_bvp_quadratic_test06 ( )

%*****************************************************************************80
%
%% FEM1D_BVP_QUADRATIC_TEST06 does an error analysis.
%
%  Discussion:
%
%    Use A6, C6, F6, EXACT6, EXACT_UX6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_BVP_QUADRATIC_TEST06\n' );
  fprintf ( 1, '  Solve -( A(x) U''(x) )'' + C(x) U(x) = F(x)\n' );
  fprintf ( 1, '  for 0 < x < 1, with U(0) = U(1) = 0.\n' );
  fprintf ( 1, '  A6(X)  = 1.0\n' );
  fprintf ( 1, '  C6(X)  = 0.0\n' );
  fprintf ( 1, '  F6(X)  = pi*pi*sin(pi*X)\n' );
  fprintf ( 1, '  U6(X)  = sin(pi*x)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compute L2 norm and seminorm of error for various N.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N        L1 error        L2 error      Seminorm error\n' );
  fprintf ( 1, '\n' );

  n = 11;
  for i = 0 : 4
%
%  Geometry definitions.
%
    x_first = 0.0;
    x_last = 1.0;
    x = linspace ( x_first, x_last, n );
    x = x(:);

    u = fem1d_bvp_quadratic ( n, @a6, @c6, @f6, x );

    e1 = l1_error ( n, x, u, @exact6 );
    e2 = l2_error_quadratic ( n, x, u, @exact6 );
    h1s = h1s_error_quadratic ( n, x, u, @exact_ux6 );

    fprintf ( 1, '  %4d  %14f  %14f  %14f\n', n, e1, e2, h1s );

    n = 2 * ( n - 1 ) + 1;

  end

  return
end
function value = a6 ( x )

%*****************************************************************************80
%
%% A6 evaluates A function #6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of A(X).
%
  value = 1.0;

  return
end
function value = c6 ( x )

%*****************************************************************************80
%
%% C6 evaluates C function #6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of C(X).
%
  value = 0.0;

  return
end
function value = exact6 ( x )

%*****************************************************************************80
%
%% EXACT6 returns exact solution #6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of U(X).
%
  value = sin ( pi * x );

  return
end
function value = exact_ux6 ( x )

%*****************************************************************************80
%
%% EXACT_UX6 evaluates the derivative of exact solution #6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of dUdX(X).
%
  value = pi * cos ( pi * x );

  return
end
function value = f6 ( x )

%*****************************************************************************80
%
%% F6 evaluates right hand side function #6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of F(X).
%
  value = pi * pi * sin ( pi * x );

  return
end


text:
function best_acc = LRSDL_top(dataset, N_train, k, k0, lambda1, lambda2, lambda3)
% function LRSDL_FDDL_top(dataset, N_train, k, k0, lambda1, lambda2, lambda3)

%% Dependencies
    addpath('utils');
    addpath('LRSDL_FDDL');
    addpath('ODL');
    %% test mode
    if nargin == 0
        dataset = 'myYaleB';
        N_train = 10;
        k = 8;
        k0 = 5;
        lambda1 = 0.001;
        lambda2 = 0.01;
        lambda3 = 0.02;
    end
    if ~exist('lambda3', 'var')
        lambda3 = 0;
    end
    %%
    t = getTimeStr();
    [dataset, Y_train, Y_test, label_train, label_test] = ...
        train_test_split(dataset, N_train);
    %% Parameter preparation
    fprintf('starting... %s\n', dataset)    ;
    [acc, rt] = LRSDL_wrapper(Y_train, label_train, Y_test , label_test, ...
                            k, k0, lambda1, lambda2, lambda3);
    fprintf('rt = %5.1f\n', rt);
    %% output filename
    if ~exist('results', 'dir')
        mkdir('results');
    end
    if k0 > 0
        if ~exist(fullfile('results', 'LRSDL'), 'dir')
            mkdir('results', 'LRSDL');
        end
        fn = fullfile('results', 'LRSDL', strcat(dataset, ...
            '_N_', num2str(N_train), '_k_', num2str(k), ...
            '_k0_',num2str(k0),'_l1_', num2str(lambda1), ...
            '_l2_', num2str(lambda2), '_l3_', num2str(lambda3), ...
            '_', t, '.mat'));
    else
        if ~exist(fullfile('results', 'FDDL'), 'dir')
            mkdir('results', 'FDDL');
        end
        fn = fullfile('results', 'FDDL', strcat(dataset, ...
            '_N_', num2str(N_train), '_k_', num2str(k),'_l1_', ...
            num2str(lambda1), '_l2_', num2str(lambda2),'_', t, '.mat'));
    end
    disp(fn);
    save(fn, 'acc', 'rt');
    best_acc = max(acc);
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function [output]=demosaicing(input)
% Bilinear Interpolation of the missing pixels
% Bayer CFA
%       R G R G
%       G B G B
%       R G R G
%       G B G B
%
% The input can be on one or three channels
%
% Output = a complete RGB image on 3 channels

im=double(input);

M = size(im, 1);
N = size(im, 2);
channel = size(size(im), 2);    % == 2 => one channel
                                % == 3 => three channel

red_mask = repmat([1 0; 0 0], M/2, N/2);
green_mask = repmat([0 1; 1 0], M/2, N/2);
blue_mask = repmat([0 0; 0 1], M/2, N/2);

if channel == 2
    R=im.*red_mask;
    G=im.*green_mask;
    B=im.*blue_mask;
elseif channel == 3
    R=im(:,:,1).*red_mask;
    G=im(:,:,2).*green_mask;
    B=im(:,:,3).*blue_mask;
end

% Interpolation for the green at the missing points
    G= G + imfilter(G, [0 1 0; 1 0 1; 0 1 0]/4);

% Interpolation for the blue at the missing points
% First, calculate the missing blue pixels at the red location
    B1 = imfilter(B,[1 0 1; 0 0 0; 1 0 1]/4);
% Second, calculate the missing blue pixels at the green locations
% by averaging the four neighouring blue pixels
    B2 = imfilter(B+B1,[0 1 0; 1 0 1; 0 1 0]/4);
    B = B + B1 + B2;

% Interpolation for the red at the missing points
% First, calculate the missing red pixels at the blue location
    R1 = imfilter(R,[1 0 1; 0 0 0; 1 0 1]/4);
% Second, calculate the missing red pixels at the green locations
    R2 = imfilter(R+R1,[0 1 0; 1 0 1; 0 1 0]/4);
    R = R + R1 + R2

    output(:,:,1)=R; output(:,:,2)=G; output(:,:,3)=B;

author: Sable

repo: mcbench-benchmarks


text:
function fk_body = autoGen_fk_body(q1,q2,q3)
%AUTOGEN_FK_BODY
%    FK_BODY = AUTOGEN_FK_BODY(Q1,Q2,Q3)

%    This function was generated by the Symbolic Math Toolbox version 8.4.
%    01-Jun-2020 11:59:03

t2 = cos(q1);
t3 = cos(q2);
t4 = sin(q1);
t5 = q2+q3;
t6 = cos(t5);
t7 = sin(t5);
t8 = t3.*3.2e+1;
t9 = t6.*5.7e+1;
t10 = t8+t9;
fk_body = reshape([t2.*t6,t4.*t6,-t7,0.0,-t4,t2,0.0,0.0,t2.*t7,t4.*t7,t6,0.0,(t2.*t10)./2.0e+2,(t4.*t10)./2.0e+2,t7.*(-5.7e+1./2.0e+2)-sin(q2).*(4.0./2.5e+1),1.0],[4,4]);

author: ShuoYangRobotics

repo: QuadrupedSim


text:
function L = fd_laplacian(side)
  % FD_LAPLACIAN  build a finite difference laplacian for a regular grid.
  %
  % L = fd_laplacian([h,w])
  % L = fd_laplacian([h,w,t])
  %
  % Inputs:
  %   dims  number of nodes along height (and width and depth)
  % Outputs:
  %   L   prod(side) by prod(side) Laplacian with negative diagonal.
  %
  % See also: fd_grad

  function B = vec(A)
    B = A(:)';
  end

  switch numel(side)
  case 2
    h = side(1);
    w = side(2);
    I = bsxfun(@plus,h*(0:w-1),(1:h)');
    L = sparse( ...
      [vec(I(2:end,:)) vec(I(:,2:end))], ...
      [vec(I(1:end-1,:)) vec(I(:,1:end-1))], ...
      1,w*h,w*h);
  case 3
    I = sub2ind(side([2 1 3]),1:prod(side))';
    I = reshape(I,side);
    L = sparse( ...
      [vec(I(2:end  ,:,:)) vec(I(:,2:end  ,:)) vec(I(:,:,2:end  ))], ...
      [vec(I(1:end-1,:,:)) vec(I(:,1:end-1,:)) vec(I(:,:,1:end-1))], ...
      1,prod(side),prod(side));
  end
  L = L + L';
  L = L-diag(sum(L,2));
end

author: alecjacobson

repo: gptoolbox


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function tfrbudt
%TFRBUDT Unit test for the function TFRBUD.

%       F. Auger, Dec. 1995 - O. Lemoine, March 1996.

% We test each property of the corresponding TFR :

N=128;

% Covariance by translation in time
t1=55; t2=70; f=0.3;
sig1=amgauss(N,t1).*fmconst(N,f,t1);
sig2=amgauss(N,t2).*fmconst(N,f,t2);
tfr1=tfrbud(sig1);
tfr2=tfrbud(sig2);
[tr,tc]=size(tfr1);
tfr=tfr1-tfr2(:,modulo((1:tc)-t1+t2,tc));
if any(any(abs(tfr)>sqrt(eps))),
 error('tfrbud test 1 failed');
end

% Reality of the TFR
sig=noisecg(N);
tfr=tfrbud(sig);
if sum(any(abs(imag(tfr))>sqrt(eps)))~=0,
 error('tfrbud test 2 failed');
end

% Energy conservation
sig=noisecg(N);
tfr=tfrbud(sig);
Es=norm(sig)^2;
Etfr=sum(mean(tfr));
if abs(Es-Etfr)>sqrt(eps),
 error('tfrbud test 3 failed');
end

% Time-marginal
sig=noisecg(N);
tfr=tfrbud(sig);
ip1=abs(sig).^2;
ip2=mean(tfr)';
if any(abs(ip1-ip2)>sqrt(eps)),
 error('tfrbud test 4 failed');
end

% Frequency-marginal
sig=noisecg(N);
tfr=tfrbud(sig,1:N,N,[1],ones(N-1,1));
FFT=fft(sig);
psd1=abs(FFT(2:N/2)).^2/(2*N);
psd2=mean(tfr(1:2:N,:)')';
if any(abs(psd1-psd2(2:N/2))>sqrt(eps)),
 error('tfrbud test 5 failed');
end

% Conservation of the time support (wide-sense)
sig=[zeros(N/4,1);fmlin(N/2);zeros(N/4,1)];
tfr=tfrbud(sig,1:N,N,[1],ones(N-1,1));
if sum(any(abs(tfr(:,1:N/4-1))>sqrt(eps))) | ...
   sum(any(abs(tfr(:,(3*N/4+1):N))>sqrt(eps))),
 error('tfrbud test 6 failed');
end

% time localization
t0=30; sig=((1:N)'==t0);
tfr=tfrbud(sig);
[ik,jk]=find(tfr~=0.0);
if any(jk~=t0)|any(ik'-(1:N)),
 error('tfrbud test 7 failed');
end;

% frequency localization
f0=10;
sig=fmconst(N+6,f0/N);
tfr=tfrbud(sig,N/2+2,N,tftb_window(11,'rect'),tftb_window(N+1,'rect'),1.2);
if (find(tfr>1/N)~=2*f0+1)|(abs(mean(tfr)-1.0)>2.0*eps),
 error('tfrbud test 8 failed');
end;

N=127;

% Covariance by translation in time
t1=55; t2=70; f=0.3;
sig1=amgauss(N,t1).*fmconst(N,f,t1);
sig2=amgauss(N,t2).*fmconst(N,f,t2);
tfr1=tfrbud(sig1);
tfr2=tfrbud(sig2);
[tr,tc]=size(tfr1);
tfr=tfr1-tfr2(:,modulo((1:tc)-t1+t2,tc));
if any(any(abs(tfr)>sqrt(eps))),
 error('tfrbud test 9 failed');
end

% Reality of the TFR
sig=noisecg(N);
tfr=tfrbud(sig);
if sum(any(abs(imag(tfr))>sqrt(eps)))~=0,
 error('tfrbud test 10 failed');
end

% Energy conservation
sig=noisecg(N);
tfr=tfrbud(sig);
Es=norm(sig)^2;
Etfr=sum(mean(tfr));
if abs(Es-Etfr)>sqrt(eps),
 error('tfrbud test 11 failed');
end

% Time-marginal
sig=noisecg(N);
tfr=tfrbud(sig);
ip1=abs(sig).^2;
ip2=mean(tfr)';
if any(abs(ip1-ip2)>sqrt(eps)),
 error('tfrbud test 12 failed');
end

% Frequency-marginal
sig=noisecg(N);
tfr=tfrbud(sig,1:N,N,[1],ones(N,1));
FFT=fft(sig);
psd1=abs(FFT(2:fix(N/2))).^2/(2*N);
psd2=mean(tfr(1:2:N,:)')';
if any(abs(psd1-psd2(2:fix(N/2)))>5e-2),
 error('tfrbud test 13 failed');
end

% Conservation of the time support (wide-sense)
sig=[zeros(round(N/4),1);fmlin(round(N/2));zeros(round(N/4),1)];
tfr=tfrbud(sig,1:N,N,[1],ones(N,1));
if sum(any(abs(tfr(:,1:round(N/4)-1))>sqrt(eps))) | ...
   sum(any(abs(tfr(:,(round(3*N/4)+2):N))>sqrt(eps))),
 error('tfrbud test 14 failed');
end

% time localization
t0=30; sig=((1:N)'==t0);
tfr=tfrbud(sig);
[ik,jk]=find(tfr~=0.0);
if any(jk~=t0)|any(ik'-(1:N)),
 error('tfrbud test 15 failed');
end;

% frequency localization
f0=10;
sig=fmconst(N+6,f0/N);
tfr=tfrbud(sig,ceil(N/2)+2,N,tftb_window(11,'rect'),tftb_window(N,'rect'),1.2);
if (find(tfr>1/N)~=2*f0+1)|(abs(mean(tfr)-1.0)>sqrt(eps)),
 error('tfrbud test 16 failed');
end;

author: HeLiangHIT


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
% use icp to align object together first
%

function [rvec, M] = match_icp(fvec,atlas,max_d);

dg = [0 max_d]; meshsize = -3;
[X, fs] = surf_spharm(atlas,dg,meshsize);
[P, fs] = surf_spharm(fvec,dg,meshsize);

% align P to X
[P,M] = align_icp(P,X);
rvec = (M(1:3,1:3)*fvec')';
rvec(1,:) = fvec(1,:) + M(1:3,4)'*2*sqrt(pi); % Y00 = 1/(2*sqrt(pi))

rmsd(1) = SPHARM_rmsd(fvec, atlas);
rmsd(2) = SPHARM_rmsd(rvec, atlas);

if rmsd(1)<rmsd(2)
    rvec = fvec;
end

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% -------------------------------------------------------------------------
function net = cnn_init_resnext(net,opts)
% -------------------------------------------------------------------------
% initialize classifier
net = dagnn.DagNN.loadobj(net) ;

% convs = find(arrayfun(@(a) isa(a.block, 'dagnn.Conv'), net.layers)==1);

fclayer = net.getLayer('classifier_0') ;
sizeW = size(net.params(fclayer.paramIndexes(1)).value);

% opts.nCls = 101;
nCls = opts.nCls ;
DropOutRate = opts.DropOutRate ;

net.params(fclayer.paramIndexes(1)).value = ...
  0.01 * randn([sizeW(1:3),nCls],'single') ;
net.params(fclayer.paramIndexes(2)).value = zeros(nCls,1,'single') ;

% change loss
softmax = find(arrayfun(@(a) isa(a.block, 'dagnn.SoftMax'), net.layers)==1);
if ~isempty(softmax)
  net.removeLayer(net.layers(softmax(1)).name) ;
end
% convs = find(arrayfun(@(a) isa(a.block, 'dagnn.Conv'), net.layers)==1);
fclayer = find(arrayfun(@(a) strcmp(a.name, 'classifier_0'), net.layers)==1);
net.renameVar(net.layers(fclayer(end)).name,'prediction') ;
net.renameVar('data','input') ;

%------------------------------------------------------------------------%
% configure appr-rank-pool
switch opts.pool1Type
  case 'arpool'
    if strcmp(opts.pool1Layer,'conv0')
      poolLyr1 = 1 ;
      net.addLayer('arpool',AppRankPooling('scale',0.1),{'input','VideoId1'},'DynImg');
      net.setLayerInputs(net.layers(1).name,{'DynImg'}) ;
    else
      poolLyr1 = find(arrayfun(@(a) strcmp(a.name, opts.pool1Layer), net.layers)==1);
      assert(~isempty(poolLyr1));
      net.addLayer('arpool',AppRankPooling('scale',0.1),{net.layers(poolLyr1).inputs{1},'VideoId1'},'DynImg');
      net.setLayerInputs(opts.pool1Layer,{'DynImg'}) ;
    end
  case 'ppool1'
    if strcmp(opts.pool1Layer,'conv0')
      poolLyr1 = 1 ;
    else
      poolLyr1 = find(arrayfun(@(a) strcmp(a.name, opts.pool1Layer), net.layers)==1);
    end
    net.addLayer('parampool',LinComb('pad',[1 1 10 1]),...
      {'features_4_0_merge','VideoId1'},'DynImg0',{'conv0f','conv0b'});

%     net.params(end-1).value = 0.1 * ones(1,1,10,1,'single');
    net.params(end-1).value = 0.1 * randn(1,1,10,1,'single');
    net.params(end).value = zeros(1,1,'single');

    net.addLayer('BnormDyn',dagnn.BatchNorm('numChannels',256),'DynImg0','DynImg',...
      {'dym','dyb','dybx'}) ;
    net.params(end-2).value =  ones(256,1,'single') ;
    net.params(end-1).value =  zeros(256,1,'single') ;
    net.params(end).value   =  zeros(256,2,'single') ;

%     net.addLayer('reluP',dagnn.ReLU(),...
%       {'DynImg1'},'DynImg');
    net.layers(16).inputs{1} = 'DynImg' ;
    for i=numel(net.params)-4:numel(net.params),
      net.params(i).learningRate = 0.1 * net.params(i).learningRate;
    end
  case 'none'
  otherwise
    error('Unknown pool type %s', opts.pool1Type) ;
end

net.rebuild() ;
%------------------------------------------------------------------------%
% second pool layer (max pooling)
% poolLyr2 = find(arrayfun(@(a) strcmp(a.name, 'pool5'), net.layers)==1);
poolLyr2 = find(arrayfun(@(a) strcmp(a.name, 'features_7_1_merge'), net.layers)==1);
net.addLayer('tempPoolMax',TemporalPooling('method','max'),...
  {net.layers(poolLyr2(1)).outputs{1},'VideoId2'},'tempPoolMax');

% change the input of fc last layer
% net.setLayerInputs(net.layers(convs(end)).name,'tempPoolMax') ;
% net.addLayer('bnar',dagnn.BatchNorm('numChannels',2048),{'tempPoolMax'},...
%   'tempPoolMaxbn',{'bnar_m','bnar_b','bnar_x'});
poolLyr2next = find(arrayfun(@(a) strcmp(a.name, 'features_7_1_id_relu'), net.layers)==1);
net.setLayerInputs(net.layers(poolLyr2next(1)).name,{'tempPoolMax'}) ;
net.rebuild() ;
%------------------------------------------------------------------------%
% add drop-out layers
if DropOutRate>0

  pool5 = find(arrayfun(@(a) strcmp(a.name, 'features_8'), net.layers)==1);
  oo = net.layers(pool5(1)).outputs{1};
  net.addLayer('drop_pool5',dagnn.DropOut('rate',DropOutRate),...
    oo,sprintf('drop_%s',oo),{});
  net.setLayerInputs('classifier_permute',{sprintf('drop_%s',oo)}) ;
end

%------------------------------------------------------------------------%
% add multi-class error
net.addLayer('errMC',ErrorMultiClass(),{'prediction','label'},'mcerr');

net.addLayer('loss', ...
             LossNormalized('loss', 'softmaxlog') ,...
             {'prediction', 'label'}, ...
             'objective') ;

%------------------------------------------------------------------------%
net.rebuild()

% replace standard matconvnet bnorm with my version
bns = find(arrayfun(@(a) strcmp(class(a.block), 'dagnn.BatchNorm'), net.layers)==1);
for i=1:numel(bns)
  bb = net.layers(bns(i)).block ;
  net.layers(bns(i)).block = BatchNormN('numChannels',bb.numChannels,...
  'epsilon',bb.epsilon,...
  'opts',bb.opts) ;
end

% dagMergeBatchNorm(net) ;
% dagRemoveLayersOfType(net, 'dagnn.BatchNorm') ;
net_ = net.saveobj ;
net = dagnn.DagNN.loadobj(net_) ;
net.meta.normalization.border = [32 32] ;

author: hbilen

repo: dynamic-image-nets


text:
function [x fx] = newtonArmijo(func,init,arguments,varargin)
%function [x fx] = newton(func,init,arguments,varargin)
% do a newton minimization
% func must return the likelihood and the value of the update -H^{-1}g

[max_its init_alpha thresh debug min_alpha max_inner_its] = process_options(varargin,'max-its',1000,'init-alpha',1,'thresh',1e-8,'debug',0,'min-alpha',1e-15,'max-inner-its',1000);
x = init;
beta = 1e-4; tau = 0.25;

iterator = newIterator(max_its,'thresh',thresh,'debug',debug);
alpha = init_alpha;
while ~iterator.done
    [old_score gradient step] = feval(func,x,arguments{:});

    %do armijo linear-searchto find step size
    new_score = feval(func,x+alpha*step,arguments{:});
    inner_its = 0;
    while inner_its < max_inner_its && alpha > min_alpha && (isnan(new_score) || isinf(new_score) || new_score > old_score + beta * gradient' * (alpha * step))
        %if debug, fprintf('alpha: %.2e -> %.2e\n',alpha,alpha*tau); end
        alpha = alpha * tau;
        new_score = feval(func,x+alpha*step,arguments{:});
        inner_its = inner_its + 1;
    end
    if alpha < min_alpha, alpha == 0; end
    try %try to take a step, if you can
        x = x + alpha * step;
        alpha = alpha ./ tau; %walk back one step
        fx(iterator.its+1) = new_score;
        iterator = updateIterator(iterator,-new_score);
    catch
        iterator.done = true;
    end
end

author: jacobeisenstein

repo: SAGE


text:
clear all;
%% Look at particular section
% Pick the gait example to observe
% load('TetrapodExample');
load('TripodExample');

% Plot the velocity of the leg tips
figure('pos',[153, 427, 560, 420]); hold on; axis tight; set(gcf,'color','w'); fontsize(16)
imagesc(example_vel);
xlabel('Time (seconds)')
xticklabels(xticks/Fs);
ylabel('Leg Tip')
yticks([1:6])
yticklabels({'RF','RM','RH','LF','LM','LH'})
ax1 = gca;
caxis([-10 10]);

% Plot the rasters
figure('pos',[850, 634, 848, 334]); hold on; axis tight; set(gcf,'color','w'); fontsize(16)
imagesc(example_stance); colormap('gray');
xlabel('Time (seconds)')
xticklabels(xticks/Fs);
ylabel('Leg Tip')
yticks([1:6])
yticklabels({'LF','LM','LH','RF','RM','RH'})
ax2 = gca;

% plot the forward velocity of the fly
figure('pos',[850, 359, 854, 186]); hold on; axis tight; set(gcf,'color','w'); fontsize(16)
plot(example_fv.*Fs)
xlabel('Time (seconds)')
xticklabels(xticks/Fs);
ylabel('Forward Velocity (mm/s)')
ax3 = gca;
ylim([0 40]);

% plot the raster of tripod, tetrapod, or non-canonical
figure('pos',[849, 114, 931, 150]); hold on; axis tight; set(gcf,'color','w'); fontsize(16)
example_gait = sum(example_stance,1);
example_gait(~(example_gait == 3 | example_gait == 4)) = 5;
imagesc(example_gait);colormap('jet');h = colorbar;
xlabel('Time (seconds)')
xticklabels(xticks/Fs);
yticks([])
ylabel(h,'Number of legs in stance')
ax4 = gca;
linkaxes([ax1,ax2,ax3,ax4],'x')

%% Plot the emission probabilities for each hidden states
load('GaitVectors3.mat');
emissions = hmm.ESTEMIT;
temp = emissions(2,:);
emissions(2,:) = emissions(3,:);
emissions(3,:) = temp;
figure; hold on;
imagesc(emissions);
for i = 1:size(emissions,1)
    for j = 1:size(emissions,2)
        caption = sprintf('%.2f',emissions(i,j));
        text(j,i,caption,'Fontsize',10,'FontWeight','bold','HorizontalAlignment','center','Color',[0 0 0]);
    end
end
axis ij;
axis tight
xlabel('Number of Legs in Stance');
yticks([1 2 3]);
yticklabels({'Tripod','Tetrapod','Non-canonical'})
xticklabels({'0','1','2','3','4','5','6'})
fontsize(16)

figure; hold on;
plot(emissions','LineWidth',3);
xlabel('Number of Legs in Stance');
xticklabels({'0','1','2','3','4','5','6'})
legend({'Tripod','Tetrapod','Non-canonical'})
ylabel('Emission Probability')

%% Plot the distribution of speeds
load('Gait_Speed_Distributions');
figure; hold on;
plot(edges1(1:end-1),N1,'LineWidth',3)
plot(edges2(1:end-1),N2,'LineWidth',3)
plot(edges3(1:end-1),N3,'LineWidth',3)
xlabel('Forward Velocity (ms)')
ylabel('Count')
xlim(speed_lim)
legend({'Tripod','Tetrapod','Non-canonical'})
fontsize(16)

%% Plot the velocity Distributions
load('Cluster_Velocity_Distributions');
figure; hold on;
cmap = spring(num_states);
for i = 1:num_states
    plot(edges{i}(1:end-1),N{i},'Color',cmap(num_states + 1 -i,:),'LineWidth',3);
end
grid on;
xlabel('Forward Velocity')
ylabel('Probability')
axis tight;
ylim([0 .2])

%% Plotting the mean and std with bounded lines
load('Swing_Velocity_Over_Time');
cmap = parula(numel(speed_levels));
p_lines = cell1(numel(speed_levels)-1);
figure('pos',[568, 186, 1036, 798]); figclosekey; set(gcf,'color','w'); hold on;
for i = 1:size(leg_vel_at_speed,1)
    yci = zeros(2,size(leg_vel_at_speed,2));
    yci(1,:) = leg_vel_std_at_speed(i,:);
    yci(2,:) = leg_vel_std_at_speed(i,:);
    [p_lines{i},~] = boundedline(win*10,leg_vel_at_speed(i,:),yci','alpha','cmap',cmap(i,:));
    p_lines{i}.LineWidth = 3;
end

% Legend
leg = cell([1 numel(speed_levels)-1]);
for i = 1:numel(speed_levels)-1
    leg{i} = sprintf('%d - %d mm/s',speed_levels(i),speed_levels(i+1));
end
l = legend([p_lines{:}],leg);
l.Position = [0.7503 0.6488 0.1573 0.3239];
fontsize(16)
xlabel('Time from swing onset (ms)')
ylabel('Swing velocity (mm/s)')
% export_fig('figs/Swing_Velocity_vs_Time_Confidences.png','-r300')

%% Plot the Swing_and_Stance_versus_Velocity
load('Swing_and_Stance_versus_Velocity')
figure; figclosekey, hold on;

yci = zeros(2,numel(stance_dur_std));
yci(1,:) = stance_dur_std;
yci(2,:) = stance_dur_std;
[bl1,~] = boundedline(stance_edges(2:end),stance_dur_mu',yci','alpha');

yci = zeros(2,numel(swing_dur_std));
yci(1,:) = swing_dur_std;
yci(2,:) = swing_dur_std;
[bl2,~] = boundedline(swing_edges(2:end),swing_dur_mu',yci','alpha','r');

yticklabels(round(yticks*10))
ylabel('Durations (ms)');
xlabel('Average Body Speed (mm/s)');
axis tight
legend([bl1,bl2],{'Stance','Swing'})
fontsize(16);

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
function p = convert_perspective_log(p)

p.kept = 1:length(p.c);
if isempty(p.evalMap)
    return
end

if ~any(p.variabletype == 4)
    return
end

variableIndex = [];
for i=1:length(p.evalMap)
    variableIndex = [variableIndex p.evalMap{i}.variableIndex];
end

removable = [];
for i = 1:length(p.evalMap)
    if strcmp(p.evalMap{i}.fcn,'log')
        argument = p.evalMap{i}.variableIndex;
        if length(argument) == 1
            if p.variabletype(argument)==4
                monoms = p.monomtable(argument,:);
                if nnz(monoms) == 2
                    k = find(monoms);
                    p1 = monoms(k(1));
                    p2 = monoms(k(2));
                    if isequal(sort([p1 p2]) , [-1 1])
                        if p2>p1
                            x = k(2);
                            y = k(1);
                        else
                            x = k(1);
                            y = k(2);
                        end
                        % Ok, so we have log(x/y)
                        % is this multiplied by x somewhere
                        logxy = p.evalMap{i}.computes;
                        enters_in = find(p.monomtable(:,logxy));
                        other = setdiff(enters_in,p.evalMap{i}.computes);
                        if length(nnz(other)) == 1
                            monomsxlog = p.monomtable(other,:);
                            if nnz(monomsxlog) == 2 & (monomsxlog(x) == 1)
                                % Hey, x*log(x/y)!
                                % we change this monomial variable to a
                                % callback variable
                                p.evalMap{i}.fcn = 'negated_perspective_log';
                                p.evalMap{i}.arg{1} = recover([x;y]);
                                p.evalMap{i}.arg{2} = [];
                                p.evalMap{i}.variableIndex = [x y];
                                p.evalMap{i}.computes = other;
                                p.evalMap{i}.properties.bounds = @nplog_bounds;
                                p.evalMap{i}.properties.convexhull = @nplog_convexhull;
                                p.evalMap{i}.properties.derivative = @nplog_derivative;
                                p.evalMap{i}.properties.inverse = [];
                                p.variabletype(other) = 0;
                                p.monomtable(other,:) = 0;
                                p.monomtable(other,other) = 1;
                                p.evalVariables(i) = other;
                                % Figure out if x/y can be removed
                                % This is possible if the x/y term never is
                                % used besides inside the log term
                                if nnz(p.F_struc(:,1+argument)) == 1 & p.c(argument) == 0 & nnz(argument == variableIndex) == 1
                                    removable = [removable argument];
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end
kept = 1:length(p.c);
kept = setdiff(kept,removable);
aux_used = zeros(1,length(p.c));
aux_used(p.aux_variables) = 1;
aux_used(removable)=[];
p.aux_variables = find(aux_used);
if length(removable) > 0
    kept = 1:length(p.c);
    kept = setdiff(kept,removable);
    [ii,jj,kk] = find(p.F_struc(:,1+removable));
    p.F_struc(:,1+removable) = [];
    p.F_struc(ii,:) = [];
    p.K.l = p.K.l - length(removable);
    p.c(removable) = [];
    p.Q(removable,:) = [];
    p.Q(:,removable) = [];
    p.variabletype(removable) = [];
    p.monomtable(:,removable) = [];
    p.monomtable(removable,:) = [];
    for i = 1:length(p.evalVariables)
        p.evalVariables(i) = find(p.evalVariables(i) == kept);
        for j = 1:length(p.evalMap{i}.variableIndex)
           p.evalMap{i}.variableIndex(j) = find(p.evalMap{i}.variableIndex(j) == kept);
        end
        for j = 1:length(p.evalMap{i}.computes)
           p.evalMap{i}.computes(j) = find(p.evalMap{i}.computes(j) == kept);
        end
    end
    p.lb(removable) = [];
    p.ub(removable) = [];
    p.used_variables(removable) = [];
end

function dp = nplog_derivative(x)
dp = [log(x(1)/x(2)) + 1;-x(1)/x(2)];

function [L,U] = nplog_bounds(xL,xU)
xU(isinf(xU)) = 1e12;
x1 = xL(1)*log(xL(1)/xL(2));
x2 = xU(1)*log(xU(1)/xL(2));
x3 = xL(1)*log(xL(1)/xU(2));
x4 = xU(1)*log(xU(1)/xU(2));
U = max([x1 x2 x3 x4]);
if (exp(-1)*xU(2) > xL(1)) & (exp(-1)*xU(2) < xU(1))
    L = -exp(-1)*xU(2);
else
    L = min([x1 x2 x3 x4]);
end

function [Ax,Ay,b,K] = nplog_convexhull(xL,xU);

x1 = [xL(1);xL(2)];
x2 = [xU(1);xL(2)];
x3 = [xL(1);xU(2)];
x4 = [xU(1);xU(2)];
x5 = (xL+xU)/2;

f1 = negated_perspective_log(x1);
f2 = negated_perspective_log(x2);
f3 = negated_perspective_log(x3);
f4 = negated_perspective_log(x4);
f5 = negated_perspective_log(x5);

df1 = nplog_derivative(x1);
df2 = nplog_derivative(x2);
df3 = nplog_derivative(x3);
df4 = nplog_derivative(x4);
df5 = nplog_derivative(x5);

[Ax,Ay,b,K] = convexhullConvex2D(x1,f1,df1,x2,f2,df2,x3,f3,df3,x4,f4,df4,x5,f5,df5);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function varargout=teager(varargin)
%
% For calling details please see v_teager.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_teager.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: teager.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_teager(varargin{:});
else
    v_teager(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function mylayer = riverplot_draw_layer(x_loc, n_rects, varargin)
%
% This is a subfunction of the riverplot toolbox.
%  - Draws a series of rects at an x location x_loc, number specified by n_rects
%  - rect1 and rect2 are structures created by riverplot_rect.m
%  - Ribbon covers a specified percentage of each rect
%  - Can start drawing from the top or bottom of each rect
%
% :Usage:
% ::
%
%     myribbon = riverplot_ribbon(rect1, rect2, [optional inputs])
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2016  Tor Wager
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
% ..
%
% :Inputs:
%
%   **rect1:**
%        Structure with handles and decription of imrect class object,
%        created with riverplot_rect.m.  This is always the LEFTMOST
%        rect of the pair.
%
%   **rect2:**
%        Same as rect1, but this is the RIGHTMOST rect of the pair
%
% :Optional Inputs:
%   **coverage:**
%        'coverage' followed by 2-element vector for [rect1 rect2], where
%        each value is between 0 and 1 and specifies the percentage of the
%        plot the ribbon will cover.
%
%   **color:**
%        followed by [r g b] color spec for ribbon
%
%   **ax:**
%        followed by axis handle to plot
%
%   **from_bottom:**
%        align ribbon with bottom edges of rects
%
% :Outputs:
%
%   **myribbon:**
%        structure with handles and info for ribbon
%
%
% :Examples:
% ::
%
% create_figure('riverplot'); set(gca, 'YDir', 'reverse');
% layer1 = riverplot_draw_layer(2, 5);
% layer2 = riverplot_draw_layer(5, 3, 'colors', bucknerlab_colors(3));
% ribbon1 = riverplot_ribbon(layer1{1}, layer2{2}, 'coverage', [.5 .2], 'color', [.2 .5 .2], 'from_bottom');
% ribbon2 = riverplot_ribbon(layer1{1}, layer2{3}, 'coverage', [1 .5], 'color', [.7 .5 .5]);
%
% :References:
%   CITATION(s) HERE
%
% :See also:
%   - list other functions related to this one, and alternatives*
%

% ..
%    Programmers' notes:
%    Created by Tor Wager, July 2, 2016
% ..

% ..
%    DEFAULTS AND INPUTS
% ..

ax = gca;
colors = custom_colors([.4 .2 0], [.8 .5 .2], n_rects); % series of cells with colors for each layer

y_loc = 0;
x_wid = .5;
y_wid = 1;

% optional inputs with default values
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case 'ax', ax = varargin{i+1}; varargin{i+1} = [];

            case {'color', 'colors'}, colors = varargin{i+1}; varargin{i+1} = [];

            case 'from_bottom', from_bottom = true;

            case 'coverage', coverage = varargin{i+1}; varargin{i+1} = [];

            case {'y_loc', 'y_start'}, y_loc = varargin{i+1}; varargin{i+1} = [];

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

axes(ax)

% set y starting point for each
totalwid = y_wid + .1*y_wid;     % y space occupied, total
y_loc = y_loc:totalwid:(y_loc + totalwid.* n_rects);

if length(colors) < n_rects
    error('Too few colors entered!');
end

for i = 1:n_rects
    mylayer{i} = riverplot_rect('color', colors{i}, 'position', [x_loc y_loc(i) x_wid y_wid]);
end

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function y2=K_q_escorTsallis(P1,P2,q)
  [M,N]=size(P1);
    y=ones(1,N);
    A1=ones(1,N);
    for n=1:N
    A1=sum(P1(:,n).^(1/q));
    A2=sum(P2(:,n).^(1/q));
       y(1,n)=sum( (P1(:,n)./((A1)^q)).*( ((P2(:,n).^(1/q))/A2).^(1-q) -((P1(:,n).^(1/q))/A1).^(1-q)) );
    end
    y2=y/(q-1);


text:
function VaR = garchkvar(data, model, ar, ma, p, q, max_forecast, alpha)
%{
-----------------------------------------------------------------------
 PURPOSE:
 Value-at-Risk Estimation for both long and short positions (Model Estimation)
-----------------------------------------------------------------------
 USAGE:
 VaR = garchvar(data, model, ar, ma, p, q, max_forecast, alpha)

 INPUTS:
 data:          (T x 1) vector of data
 model:         'GARCH', 'GJR', 'AGARCH','NAGARCH'
 ar:            positive scalar integer representing the order of AR
 am:            positive scalar integer representing the order of MA
 p:             positive scalar integer representing the order of ARCH
 q:             positive scalar integer representing the order of GARCH
 max_forecasts: maximum number of forecasts (i.e. 1-trading months 22 days)
 alpha:         confidence level

 OUTPUTS:
 VaR:           a vector of VaR forecasts
-----------------------------------------------------------------------
 Author:
 Alexandros Gabrielsen, a.gabrielsen@city.ac.uk
 Date:     11/2011
-----------------------------------------------------------------------
%}

if nargin == 0
    error('Data, GARCH Model, AR, MA, ARCH, GARCH, Maximum Number of Forecasts, Confidence Level')
end

if size(data,2) > 1
   error('Data vector should be a column vector')
end

if (length(ar) > 1) | (length(ma) > 1) | ar < 0 | ma < 0
    error('AR and MA should be positive scalars')
end

if (length(p) > 1) | (length(q) > 1) | p < 0 | q < 0 | alpha < 0 | alpha > 1
    error('P,Q and alpha should be positive scalars')
end

% Estimate the model
[parameters, stderrors, LLF, ht, kt, nu, resids] = garchk(data, strcat(model), ar, ma, 0, p, q, 0);

% Calling garchkfor2 function to pass the parameters
[MF, VF, KF, MC, VC] = garchkfor2(data, resids, ht, kt, parameters, strcat(model),ar, ma, p, q, max_forecast);

% Estimating Inverce-CDF using the degrees of freedom which are described
% as a function of the forecasted conditional kurtosis
cdfqtile1 = tinv(alpha,(4*KF-6)./(KF-3));
cdfqtile2 = tinv(1-alpha,(4*KF-6)./(KF-3));

% Estimating VaR
VaR = [MC + cdfqtile1.*VC, MC + cdfqtile2.*VC];

end


text:
function [ind,d] = find(v,w,epsilon,varargin)
% return index of all points in a epsilon neighborhood of a vector
%
% Syntax
%   ind = find(v,w,epsilon) % find all points out of v in a epsilon neighborhood of w
%   ind = find(v,w)         % find closest point out of v to w
%
% Input
%  v, w    - @vector3d
%  epsilon - double
%
% Options
%  antipodal - include <VectorsAxes.html antipodal symmetry>
%
% Output
%  ind     - int32

% compute distances
d = angle_outer(v,w,varargin{:});

% find neigbours
if nargin >= 3
  if epsilon == 1
    [d,ind] = min(d,[],1);
  else
    ind = d < epsilon;
  end
else
  [d,ind] = min(d,[],1);
end

author: mtex-toolbox

repo: mtex


text:
function K = sdlfmaXsdlfmvKernComputeBlock(lfmKern1, lfmKern2, t1, t2, ...
    kyy, kyv, kvy, kvv, i, j, generalConst)

% SDLFMAXSDLFMVKERNCOMPUTEBLOCK Computes SDLFM kernel matrix for block i,j
% FORMAT
% DESC computes the kernel matrix for the SDLFM kernel function in the
% block specified at indeces i,j. It assumes the computation for functions
% that describe acceleration (system 1) and velocity (system 2)
% ARG lfmKern1 : structure containing parameters for the system 1
% ARG lfmKern2 : structure containing parameters for the system 2
% ARG t1 : times at which the system 1 is evaluated
% ARG t2 : times at which the system 2 is evaluated
% ARG kyy : covariance for the initial conditions between position 1 and
% position 2 at block i,j
% ARG kyv : covariance for the initial conditions between position 1 and
% velocity 2 at block i,j
% ARG kvy : covariance for the initial conditions between velocity 1 and
% position 2 at block i,j
% ARG kvv : covariance for the initial conditions between velocity 1 and
% velocity 2 at block i,j
% ARG i : interval to be evaluated for system 1
% ARG j : interval to be evaluated for system 2
% ARG generalConstant : constants evaluated with sdlfmKernComputeConstant.m
% RETURN K : the kernel matrix portion of block i,j
%
% COPYRIGHT : Mauricio A. Alvarez, 2010.

% KERN

if nargin<11
    j = i;
    generalConst = [];
end

a1 = sdlfmaMeanCompute(lfmKern1(1), t1, 'Pos');
b1 = sdlfmaMeanCompute(lfmKern1(1), t1, 'Vel');
g2 = sdlfmvMeanCompute(lfmKern2(1), t2, 'Pos');
h2 = sdlfmvMeanCompute(lfmKern2(1), t2, 'Vel');

K = kyy*a1*g2.' + kyv*a1*h2.' + kvy*b1*g2.' + kvv*b1*h2.';

if i==j
    for k=1:length(lfmKern1)
        K  = K + lfmaXlfmvKernCompute(lfmKern1(k), lfmKern2(k), t1, t2);
    end
else
    if i>j
        PosVel = zeros(1, length(t2));
        VelVel = zeros(1, length(t2));
        for k=1:length(lfmKern1)
            PosVel = PosVel + lfmvXlfmKernCompute(lfmKern2(k), lfmKern1(k), t2, lfmKern2(k).limit).';
            VelVel = VelVel + lfmvXlfmvKernCompute(lfmKern1(k), lfmKern2(k), lfmKern2(k).limit, t2);
        end
        if isempty(generalConst{i,j})
            K = K + a1*PosVel + b1*VelVel;
        else
            K = K + (generalConst{i,j}(1,1)*a1 + generalConst{i,j}(2,1)*b1)*PosVel + ...
                (generalConst{i,j}(1,2)*a1 + generalConst{i,j}(2,2)*b1)*VelVel;
        end
    else
        AccelPos = zeros(length(t1),1);
        AccelVel = zeros(length(t1),1);
        for k =1:length(lfmKern1)
            AccelPos = AccelPos + lfmaXlfmKernCompute(lfmKern1(k), lfmKern2(k), t1, lfmKern1(k).limit);
            AccelVel = AccelVel + lfmaXlfmvKernCompute(lfmKern1(k), lfmKern2(k), t1, lfmKern1(k).limit);
        end
        if isempty(generalConst{i,j})
            K = K + AccelPos*g2.' + AccelVel*h2.';
        else
            K = K + AccelPos*(generalConst{i,j}(1,1)*g2.' + generalConst{i,j}(2,1)*h2.') + ...
                AccelVel*(generalConst{i,j}(1,2)*g2.' + generalConst{i,j}(2,2)*h2.');
        end
    end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text: %RDIVIDE Right dataset divide. Dataset overload

author: marianux

repo: ecg-kit

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


text:
% eeg_compare - compare EEG structures. Output are shown on the command
%               line regarding difference between the files.
% Usage:
%   eeg_compare(EEG1, EEG2);
%
% Input:
%  EEG1 - first EEGLAB structure
%  EEG2 - second EEGLAB structure

function eeg_compare(EEG, EEG2)

%% Assess difference between datasets
fields = fieldnames(EEG);
for iField = 1:length(fields)
    if ~isfield(EEG2, fields{iField})
        fprintf('Field %s missing in second dataset\n', fields{iField});
    else
        if ~isequal(EEG.(fields{iField}), EEG2.(fields{iField}))
            fprintf('Field %s differs\n', fields{iField});
        end
    end
end
if ~isequal(EEG.xmin, EEG2.xmin), fprintf('Difference between xmin is %1.6f sec\n', EEG.xmin-EEG2.xmin); end
if ~isequal(EEG.xmax, EEG2.xmax), fprintf('Difference between xmax is %1.6f sec\n', EEG.xmax-EEG2.xmax); end

% check chanlocs
[~,~,chanlocs1] = eeg_checkchanlocs( EEG.chanlocs, EEG.chaninfo);
[~,~,chanlocs2] = eeg_checkchanlocs( EEG2.chanlocs, EEG.chaninfo);
if length(chanlocs1) == length(chanlocs2)
    differ = 0;
    differLabel = 0;
    for iChan = 1:length(chanlocs1)
        if  sum(abs([ chanlocs1(iChan).X chanlocs1(iChan).Y chanlocs1(iChan).Z] - ...
                [ chanlocs2(iChan).X chanlocs2(iChan).Y chanlocs2(iChan).Z])) > 1e-12
            differ = differ+1;
        end
        if ~isequal(chanlocs1(iChan).labels, chanlocs2(iChan).labels)
            differLabel = differLabel+1;
        end
    end
    if differ
        fprintf('%d channel coordinates differ\n', differ);
    else
        disp('All channel coordinates are OK');
    end
    if differLabel
        fprintf('%d channel label(s) differ\n', differLabel);
    else
        disp('All channel labels are OK');
    end
else
    disp('Different numbers of channels');
end

% check events
if length(EEG.event) ~= length(EEG2.event)
    disp('Different numbers of events');
elseif isempty(EEG.event)
    disp('All events OK (empty)');
else
    fields1 = fieldnames(EEG.event);
    fields2 = fieldnames(EEG2.event);
    allFieldsOK = true;

    if ~isequal(sort(fields1), sort(fields2))
        disp('Not the same number of event fields');
        allFieldsOK = false;
    end

    for iField = 1:length(fields1)
        if isfield(EEG.event, fields1{iField}) && isfield(EEG2.event, fields1{iField})
            diffVal = zeros(1,length(EEG.event));
            if strcmpi(fields1{iField}, 'latency')
                for iEvent = 1:length(EEG.event)
                    diffVal(iEvent) = EEG.event(iEvent).(fields1{iField}) - EEG2.event(iEvent).(fields1{iField});
                end
            else
                for iEvent = 1:length(EEG.event)
                    diffVal(iEvent) = ~isequal(EEG.event(iEvent).(fields1{iField}), EEG2.event(iEvent).(fields1{iField}));
                end
            end
            if any(diffVal ~= 0)
                if strcmpi(fields1{iField}, 'latency')
                    fprintf('Event latency (%2.1f %%) are not OK (abs diff of these is %1.4f samples)\n', length(find(diffVal))/length(diffVal)*100, mean( abs(diffVal(diffVal ~=0 ))) );
                    fprintf(' ******** (see plot)\n');
                    figure; plot(diffVal);
                else
                    fprintf('Event fields "%s" are NOT OK (%2.1f %% of them)\n', fields1{iField}, length(find(diffVal))/length(diffVal)*100);
                end
                allFieldsOK = false;
            end
        end
    end
    disp('All other events OK (unless warning above)');
end

% check epochs
if length(EEG.epoch) == length(EEG2.epoch)
    if ~isempty(EEG.epoch)
        fields1 = fieldnames(EEG.epoch);
        fields2 = fieldnames(EEG2.epoch);
        allFieldsOK = true;
        if ~isequal(sort(fields1), sort(fields2))
            disp('Not the same number of event fields');
            allFieldsOK = false;
        else
            diffVal = [];
            for iField = 1:length(fields1)
                for iEpoch = 1:length(EEG.epoch)
                    diffVal(iEpoch) = ~isequal(EEG.epoch(iEpoch).(fields1{iField}), EEG2.epoch(iEpoch).(fields1{iField}));
                end
                if any(diffVal ~= 0)
                    fprintf('Epoch fields "%s" are NOT OK (%2.1f %% of them)\n', fields1{iField}, length(find(diffVal))/length(diffVal)*100);
                    allFieldsOK = false;
                end
            end
        end
        if allFieldsOK
            disp('All epoch and all epoch fields are OK');
        end
    end
elseif ~isempty(EEG.epoch)
    disp('Different numbers of epochs');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
close all;
clear all;
clc;

data_file_path = 'bin/mmv_phase_transition_snr_40db_s_4.mat';
options.export = true;
options.export_dir = 'bin';
options.export_name = 'mmv_snr_40_db_s_4';
options.chosen_ks = [2, 4, 8, 16, 32, 64];
options.subtitle = 'MMV, SNR=40dB s=4';
spx.pursuit.PhaseTransitionAnalysis.print_results(data_file_path, ...
    'CoSaMP', options);

author: indigits

repo: sparse-plex


text:
%%
clear all; close all; clc

%% set MATLAB path
% set the Explore2p code folder as current directory
if true % set the
    scriptName = mfilename('fullpath');
    [code_dir, filename, fileextension]= fileparts(scriptName);
    addpath(genpath(code_dir));
    cd(code_dir)

else
    % (hard-coded, customize here)
    folder = 'C:\Users\xiuye\Dropbox\!Research\2Pcode\Explore2p'; %#ok<UNRCH>
    addpath(genpath(folder));
    cd(folder)
end

%% Automatically load demo data
% isDemo = true;
%
% %% Start the GUI
% h = Explore2p(isDemo);

% to update the object handle in the workspace, go to the GUI Menu and press
% File\Export to workspace. See demo\demo_workspace_interactive.m

%% For Scripts: see demo\demo_script.m for using Explore2p in custom scripts

%% To load the GUI with non-demo data, simply call -
% Explore2p;

author: xiuyechen

repo: FishExplorer


text:
% Saving Movie Files Example
%
% This example demonstrates how to save the simulation animations as a
% movie. It builds on the Heterogeneous Propagation Medium Example.
%
% author: Bradley Treeby
% date: 30th June 2009
% last update: 17th October 2011
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave 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 Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all;

% create the computational grid
Nx = 128;           % number of grid points in the x (row) direction
Ny = 128;           % number of grid points in the y (column) direction
dx = 0.1e-3;        % grid point spacing in the x direction [m]
dy = 0.1e-3;        % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium
medium.sound_speed = 1500*ones(Nx, Ny); % [m/s]
medium.sound_speed(1:Nx/2, :) = 1800;   % [m/s]
medium.density = 1000*ones(Nx, Ny);     % [kg/m^3]
medium.density(:, Ny/4:end) = 1200;     % [kg/m^3]

% create initial pressure distribution using makeDisc
disc_magnitude = 5; % [Pa]
disc_x_pos = 50;    % [grid points]
disc_y_pos = 50;    % [grid points]
disc_radius = 8;    % [grid points]
disc_1 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

disc_magnitude = 3; % [Pa]
disc_x_pos = 80;    % [grid points]
disc_y_pos = 60;    % [grid points]
disc_radius = 5;    % [grid points]
disc_2 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

source.p0 = disc_1 + disc_2;

% define a centered circular sensor
sensor_radius = 4e-3;   % [m]
num_sensor_points = 50;
sensor.mask = makeCartCircle(sensor_radius, num_sensor_points);

% define the input arguments
input_args = {'PlotPML', false, 'RecordMovie', true, 'MovieName', 'example_movie_1', 'MovieType', 'image', 'PlotFreq', 5, 'MovieArgs', {'fps', 30}};

% run the simulation
kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

% define a second set of input arguments
input_args = {'MeshPlot', true, 'DisplayMask', 'off', 'PlotPML', false, 'RecordMovie', true, 'MovieName', 'example_movie_2', 'PlotFreq', 5, 'MovieArgs', {'fps', 30}};

% run the simulation without a sensor input
kspaceFirstOrder2D(kgrid, medium, source, [], input_args{:});

repo: Aero-matlab


text:
function out = MF_armax(y, orders, pTrain, numSteps)
% MF_armax  Statistics on a fitted ARMA model.
%
% Uses the functions iddata, armax, aic, and predict from Matlab's System
% Identification Toolbox
%
%---INPUTS:
%
% y, the input time series
%
% orders, a two-vector for p and q, the AR and MA components of the model,
%           respectively
%
% pTrain, the proportion of data to train the model on (the remainder is used
%           for testing)
%
% numSteps, number of steps to predict into the future for testing the model.
%
%
%---OUTPUTS: include the fitted AR and MA coefficients, the goodness of fit in
% the training data, and statistics on the residuals from using the fitted model
% to predict the testing data.

% ------------------------------------------------------------------------------
% 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/>.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%% Check that a System Identification Toolbox license is available:
% ------------------------------------------------------------------------------
BF_CheckToolbox('identification_toolbox')

% ------------------------------------------------------------------------------
%% Prepare Inputs
% ------------------------------------------------------------------------------
% (1) y, the time series as a column vector
if size(y,2) > size(y,1)
   y = y'; % ensure a column vector
end
N = length(y); % number of samples
% Convert y to time series object
y = iddata(y,[],1);

% orders; vector specifying the AR and MA components
if nargin < 2 || isempty(orders)
    orders = [3, 3]; % AR3, MA3
end
if nargin < 3 || isempty(pTrain)
    pTrain = 0.8; % train on 80% of the data
end
% if nargin < 4 || isempty(trainmode)
%     trainmode = 'first'; % trains on first pTrain proportion of the data.
% end
if nargin < 4 || isempty(numSteps)
    numSteps = 1; % one-step-ahead predictions
end

% ------------------------------------------------------------------------------
%% Fit the model
% ------------------------------------------------------------------------------

% Uses the System Identification Toolbox function armax
m = armax(y, orders);

% ------------------------------------------------------------------------------
%% Statistics on model
% ------------------------------------------------------------------------------

c_ar = m.a; % AR coefficients
c_ma = m.c; % MA coefficients
da = m.da; % must be uncertainties in AR coeffs
dc = m.dc; % must uncertainties in MA coeffs

% Make these outputs
if length(c_ar) > 1
    for i = 2:length(c_ar)
        out.(sprintf('AR_%u',i-1)) = c_ar(i);
    end
end
if length(c_ma) > 1
    for i = 2:length(c_ma)
        out.(sprintf('MA_%u',i-1)) = c_ma(i);
    end
end

if isempty(da)
    out.maxda = NaN;
else
    out.maxda = max(da);
end
if isempty(dc)
    out.maxdc = NaN;
else
    out.maxdc = max(dc);
end

% ------------------------------------------------------------------------------
% Fit statistics
% ------------------------------------------------------------------------------

% These three measures are basically equivalent -- default hctsa library
% only records fpe.
out.noisevar = m.NoiseVariance; % covariance matrix of noise source
% covmat = m.CovarianceMatrix; % covariance matrix for parameter vector
% parameters = m.ParameterVector; % parameter vector for model: initial values, I'd say...
out.lossfn = m.EstimationInfo.LossFcn;
out.fpe = m.EstimationInfo.FPE; % Final prediction error of model

% out.lastimprovement = m.EstimationInfo.LastImprovement; % Last improvement made in iteration
out.aic = aic(m); % ~ log(fpe)

% ------------------------------------------------------------------------------
%% Prediction
% ------------------------------------------------------------------------------

% Select first portion of data for estimation
% This could be any portion, actually... Maybe could look at robustness of
% model to different training sets...
ytrain = y(1:floor(pTrain*N));
% ytest = y;
ytest = y(floor(pTrain*N):end); % overlap

% Train the model on just this portion
mp = armax(ytrain, orders);

% Compute step-ahead predictions
% Maybe look at trends across different prediction horizons...
yp = predict(mp, ytest, numSteps, 'init', 'e'); % across whole dataset

mresiduals = ytest.y - yp.y;

% ------------------------------------------------------------------------------
% Get statistics on residuals
% ------------------------------------------------------------------------------
residout = MF_ResidualAnalysis(mresiduals);

% Convert these to local outputs in quick loop
% Note that default hctsa library does not include rmse field, which is highly
% correlated with the stde field
fields = fieldnames(residout);
for k = 1:length(fields);
    out.(fields{k}) = residout.(fields{k});
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
f{1}=dir('*.m')
f{2}=dir('visim/*.m')

fid=fopen('segymat-functions.xml','w');

fprintf(fid,'<sect1 id=\"%s\"><title>%s</title>\n','reference','M-file Reference');
for ff=1:length(f)

     for i=1:length(f{ff})
     [p,name,ext]=fileparts(f{ff}(i).name);
     disp(name)
     h=help(name);

     fprintf(fid,'<sect2 id=\"%s\"><title>%s</title>\n',['_',name],name);
     fprintf(fid,'<para><programlisting><![CDATA[%s]]></programlisting></para>\n',h);
     fprintf(fid,'</sect2>\n\n');

     end

end
fprintf(fid,'</sect1>');

fclose(fid);

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat