text:
function testRB1
%testRB1: test the RBio toolbox.
%
% Example:
%   testRB1
%
% See also UFget, RBread, RBreade, testRB2.

% Copyright 2006, Timothy A. Davis

files = {
'bcsstk01.rb'
'farm.rb'
'lap_25.pse'
'lap_25.rb'
'west0479.rb'
'west0479.rua' } ;

for k = 1:length(files)
    file = files {k} ;
    fprintf ('%s : ', file) ;
    if (file (end) == 'e')
	[A Z] = RBreade (file) ;
    else
	[A Z] = RBread (file) ;
    end
    fprintf ('%s\n', RBtype (A)) ;
    RBwrite ('temp.rb', A, Z) ;
    [A2 Z2] = RBread ('temp.rb') ;
    if (~isequal (A, A2))
	fprintf ('test failed: %s (A differs)\n', file) ;
	error ('!') ;
    end
    if (~isequal (Z, Z2))
	fprintf ('test failed: %s (Z differs)\n', file) ;
	error ('!') ;
    end
end

load west0479
C = west0479 ;
RBwrite ('mywest', C, 'WEST0479 chemical eng. problem', 'west0479') ;
A = RBread ('mywest') ;
if (~isequal (A, C))
    error ('test failed: west0479 (MATLAB version)') ;
end
fprintf ('west0479 (MATLAB matrix) : %s\n', RBtype (A)) ;

if (~strcmp (RBtype (A), 'rua'))
    error ('test failed: RBtype(A)\n') ;
end
if (~strcmp (RBtype (spones (A)), 'pua'))
    error ('test failed: RBtype(spones(A))\n') ;
end
if (~strcmp (RBtype (2*spones (A)), 'iua'))
    error ('test failed: RBtype(2*spones(A))\n') ;
end
C = A+A' ;
if (~strcmp (RBtype (C), 'rsa'))
    error ('test failed: RBtype(A+A'')\n') ;
end

delete ('mywest') ;
delete ('temp.rb') ;
fprintf ('testRB1: all tests passed\n') ;

author: fieldtrip

repo: fieldtrip


text:
function[varargout]=materncov(varargin)
%MATERNCOV  Autocovariance of the Matern random process and variations.
%
%   [TAU,R]=MATERNCOV(DT,N,SIGMA,ALPHA,LAMBDA) returns the autocovariance
%   function R of a length N complex-valued Matern random process having
%   variance SIGMA^2, slope parameter ALPHA, and damping parameter LAMBDA.
%
%   [TAU,R]=MATERNCOV(...,'real') instead forms the covariance of a real-
%   valued Matern process.
%
%   DT is the sample interval.  Note that LAMBDA is understood to have the
%   same units as the inverse sample interval 1/DT.
%
%   TAU is an array of time lags at which R is computed, and is given by
%   TAU=DT*[0,1,...,N-1].
%
%   By definition, R is one-sided theoretical autocovariance at
%   non-negative time lags.  See below for the relationship between this
%   and the full, length (2N-1) theoretical autocovariance function.
%
%   Note that for LAMBDA=0, the case of fractional Brownian motion, R will
%   contain only INFs because the autocovariance function is unbounded.
%
%   The input parameters SIGMA, ALPHA, and LAMBDA, may all either be
%   scalars or arrays of the same length M.  If the latter, then the output
%   autocovariance function R will be a matrix with N rows and M columns.
%
%   [TAU,R]=MATERNCOV(DT,N,SIGMA,ALPHA,LAMBDA,NU,MU) returns the
%   autocovariance function of various extensions of the Matern process.
%   MATERNOISE(...,'composite') also works.  See MATERNSPEC for details.
%
%   See MATERNSPEC for a more thorough discussion of the Matern process.
%
%   For details on the Matern process and its autocovariance function, see:
%
%     Lilly, Sykulski, Early, and Olhede, (2017).  Fractional Brownian
%        motion, the Matern process, and stochastic modeling of turbulent
%        dispersion.  Nonlinear Processes in Geophysics, 24: 481--514.
%   __________________________________________________________________
%
%   Relationship to full autocovariance
%
%   For a time series of length N, the full autocovariance function RF is
%   length 2N-1, defined at time lags -N+1,-N+2...,-1,0,1,...,N-2,N-1.
%
%   The one-sided autocovariance R contains the full autocovariance RF at
%   positive time lags. Negative lags are given by Hermitian symmetry.
%
%   [TAUF,RF]=MATERNCOV(...,'full') returns the full (two-sided)
%   autocovariance RF and the corresponding two-sided time array TAUF.
%
%   RF is constructed from R as RF=[FLIPUD(CONJ(R(2:end,:));R].
%   __________________________________________________________________
%
%   Real-valued processes
%
%   By default MATERNCOV returns the autocovariance of a complex-valued
%   process.
%
%   MATERNCOV(...,'real') instead returns the autocovariance of a real-
%   valued process. This also works with any of extended versions.
%   __________________________________________________________________
%
%   Computational notes
%
%   The autocovariances for the generalized and composite Matern spectra do
%   not have analytic forms. Rather, these are approximated to high
%   precision by inverse Fourier transforming the oversampled spectrum with
%   10 x oversampling over a 10 x longer time period, and then decimating.
%
%   MATERNCOV(...,'general',M,P) or  MATERNCOV(...,'composite',M,P)
%   specifies the numerical oversampling parameters in the calculation.
%
%   The spectrum is then computed over a time window of M times the
%   required  duration, and P times the required sampling density, for a
%   total of M*P time more points.  These values may set to optimize the
%   tradeoff between speed and accuracy.
%
%   This computation method is expected to minimize aliasing effects and
%   resolution errors.
%   __________________________________________________________________
%
%   See also MATERNSPEC, MATERNIMP, MATERNOISE, MATERNFIT.
%
%   'materncov --t' runs some tests.
%
%   Usage: [tau,R]=materncov(dt,N,sigma,alpha,lambda);
%          [tau,R]=materncov(dt,N,sigma,alpha,lambda,nu);
%          [tau,R]=materncov(dt,N,sigma,alpha,lambda,mu,'composite');
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2013--2021 J.M. Lilly --- type 'help jlab_license' for details

%   No longer supported, sorry
%   __________________________________________________________________
%
%   Reverse Matern
%
%   [TAU,R]=MATERNCOV(N,SIGMA,ALPHA,LAMBDA,'reverse') returns the autocovariance
%   of the reverse Matern process, in which the roles of the spectrum and
%   the autocovariance functions are swapped.  The autocovariance is
%
%      R(TAU) = SIGMA^2 / (1+LAMBDA^2*TAU^2)^ALPHA
%
%   Note that the parameter LAMBDA has been inverted so that it still has units
%   of frequency, as in the case of the forward Matern process.

%   __________________________________________________________________
%
%   Extensions
%
%   Several extensions of the basic Matern form are supported, all of which
%   are discussed in more detail in MATERNSPEC.
%
%   Oscillatory Matern
%       [F,SPP,SNN]=MATERNCOV(DT,N,SIGMA,ALPHA,LAMBDA,NU)
%   Generalized Matern (+ optional oscillations)
%       [F,SPP,SNN]=MATERNCOV(DT,N,SIGMA,ALPHA,LAMBDA,GAMMA,'general')
%       [F,SPP,SNN]=MATERNCOV(DT,N,SIGMA,ALPHA,LAMBDA,GAMMA,NU,'general')
%   Extended Matern  (+ optional oscillations)
%       [F,SPP,SNN]=MATERNCOV(DT,N,SIGMA,ALPHA,LAMBDA,MU,'extended')
%       [F,SPP,SNN]=MATERNCOV(DT,N,SIGMA,ALPHA,LAMBDA,MU,NU,'extended')
%   Composite Matern
%       [F,SPP,SNN]=MATERNCOV(DT,N,SIGMA,ALPHA,LAMBDA,MU,NU,'composite')

if strcmpi(varargin{1}, '--t')
    materncov_test,return
end

sid='one';
model='standard';
flag='complex';

M=10;  %Specifying oversampling rates for numerical computations
P=10;
for i=1:3
    if ischar(varargin{end})
        if strcmpi(varargin{end}(1:3),'rev')||strcmpi(varargin{end}(1:3),'com')||strcmpi(varargin{end}(1:3),'gen')||strcmpi(varargin{end}(1:3),'ext')
            model=varargin{end}(1:3);
        elseif strcmpi(varargin{end}(1:3),'rea')
            flag=varargin{end};
        else
            sid=varargin{end};
        end
        varargin=varargin(1:end-1);
    elseif ischar(varargin{end-2})  %For inputting M and P for testing
        if strcmpi(varargin{end-2}(1:3),'com')||strcmpi(varargin{end-2}(1:3),'gen')
            model=varargin{end-2};
            M=varargin{end-1};
            P=varargin{end};
            varargin=varargin(1:end-3);
        end
    end
end

if strcmpi(model(1:3),'com')&&strcmpi(model(1:3),'rev')
    error('Sorry, the COMPOSITE and REVERSE options cannot be used together.')
end

dt=double(varargin{1});
N=floor(varargin{2});
varargin=varargin(3:end);

if length(varargin)==1
    params=varargin{1};  %for input from materfit
else
    %this will be easier as a matrix
    Ncell=max(cellength(varargin));
    params=zeros(Ncell,length(varargin));
    for i=1:length(varargin)
        params(:,i)=varargin{i};
    end
end

tau=dt*[0:N-1]';
R=zeros(N,size(params,1));
for i=1:size(params,1)
    R(:,i)=materncov_one(tau,params(i,:),model,M,P,flag);
end

if strcmpi(sid(1:3),'ful')
    tau=[-flipud(tau(2:end,:));tau];
    R=[flipud(conj(R(2:end,:)));R];
end
varargout{1}=tau;
varargout{2}=R;

function[R]=materncov_one(tau,params,model,M,P,flag)

sigma=params(1);
alpha=params(2);
lambda=params(3);

nu=0;
if strcmpi(model(1:3),'com')
    mu=params(4);
    nu=params(5);

    N=length(tau);
    dt=tau(2)-tau(1);
    [f,Spp,Snn]=maternspec(dt,M*N*P,sigma,alpha,lambda/P,mu/P,nu/P,'composite',flag);
    S=[flipud(Snn(2:end));Spp];%figure,plot(S)
    Ri=ifft(ifftshift(S))./dt;  %Make sure it's ifftshift not fftshift
    Ri=Ri(1:P:end);
    R=Ri(1:N);
    R=R./R(1);
elseif strcmpi(model(1:3),'gen')
    %Generalized Matern form
    gamma=params(4);
    if length(params)==5
        nu=params(5);
    end
    %sigma,alpha,lambda/P,gamma,nu,flag
    N=length(tau);
    dt=tau(2)-tau(1);
    [f,Spp,Snn]=maternspec(dt,M*N*P,sigma,alpha,lambda/P,gamma,'generalized',flag);
    S=[flipud(Snn(2:end));Spp];%figure,plot(S)
    %figure,plot(f,[Spp Snn]),xlog,ylog,aresame(Spp,Snn,1e-13)
    Ri=ifft(ifftshift(S))./dt;  %Make sure it's ifftshift not fftshift
    %figure,plot(Ri)
    Ri=Ri(1:P:end);
    R=Ri(1:N);
    R=R./R(1);
elseif strcmpi(model(1:3),'ext')
    %Extended Matern, reduces to the standard Matern form for mu=0
    mu=params(4);
    if length(params)==5
        nu=params(5);
    end
    t=lambda.*sqrt(tau.^2+mu.^2);  %This is to support the extended Matern form
    R=frac(maternfun(alpha,t),maternfun(alpha,mu*lambda));
elseif strcmpi(model(1:3),'sta')
    %Standard Matern form
    if length(params)==4
        nu=params(4);
    end
    t=lambda.*tau;
    R=maternfun(alpha,t);
end
if nu~=0 && ~strcmpi(model(1:3),'com')
    R=R.*exp(sqrt(-1)*tau.*nu);
end
R=R.*sigma.^2;

if strcmpi(flag(1:3),'rea')
     R=real(R);
end

% Extended Matern form reduces to this Bessel function form, test is below
% if alpha=-1/2
%     tnorm=sqrt(tau.^2+mu.^2);
%     fact=besselk(1,lambda.*mu);
%     R1=frac(1,fact).*frac(1,tnorm./mu).*besselk(1,lambda.*tnorm);
% end

% function[]=maternfun_fig
%
% alpha=[-4:0.1:4];
% [tau,R]=materncov(dt,10000,1,alpha,0,0,1000);
% figure,plot(tau,R)
%
% hold on,h=plot(r1,xi(:,1)); linestyle -h h 3D
% hold on,h=plot(r1,xi(:,end)); linestyle -h h 3k
% title('The Matern function from \alpha = 1/2 (gray) to \alpha = 4 (black)')

function[]=materncov_test

tic
alpha=1.5;
h=0.1;
sigma=1;

N=1000;
[tau,R]=materncov(1,N,sigma,alpha,h);
[f,Spp,Snn]=maternspec(1,N,sigma,alpha,h);
[f,Spp2,Snn2]=blurspec(1,R,'aliased');

b1=aresame(Spp2./maxmax(Spp),Spp./maxmax(Spp),1e-4);
b2=aresame(Snn2./maxmax(Snn),Snn./maxmax(Snn),1e-4);

reporttest('MATERNCOV Fourier transforms to MATERNSPEC for case with negligible aliasing, even N',b1&&b2)

N=1000;
[tau,R]=materncov(1,N,sigma,alpha,h,h/2);
[f,Spp,Snn]=maternspec(1,N,sigma,alpha,h,h/2);
[f,Spp2,Snn2]=blurspec(1,R,'aliased');

b1=aresame(Spp2./maxmax(Spp),Spp./maxmax(Spp),1e-4);
b2=aresame(Snn2./maxmax(Snn),Snn./maxmax(Snn),1e-4);

reporttest('MATERNCOV Fourier transforms to MATERNSPEC for oscillatory Matern with negligible aliasing, even N',b1&&b2)

N=1000;
[tau,R]=materncov(1,N,sigma,alpha,h,h/2,'real');
[f,Spp,Snn]=maternspec(1,N,sigma,alpha,h,h/2,'real');
[f,Spp2,Snn2]=blurspec(1,R,'aliased');

b1=aresame(Spp2./maxmax(Spp),Spp./maxmax(Spp),1e-4);
b2=aresame(Snn2./maxmax(Snn),Snn./maxmax(Snn),1e-4);

reporttest('MATERNCOV Fourier transforms to MATERNSPEC for real-valued oscillatory Matern with negligible aliasing, even N',b1&&b2)

N=1000;
dt=7;
[tau,R]=materncov(dt,N,sigma,alpha,h/dt);
[f,Spp,Snn]=maternspec(dt,N,sigma,alpha,h/dt);
[f,Spp2,Snn2]=blurspec(dt,R,'aliased');

b1=aresame(Spp2./maxmax(Spp),Spp./maxmax(Spp),1e-4);
b2=aresame(Snn2./maxmax(Snn),Snn./maxmax(Snn),1e-4);

reporttest('MATERNCOV Fourier transforms to MATERNSPEC ..., even N, non-unit sample rate',b1&&b2)

N=1000;
[f,Spp,Snn]=maternspec(2,N,sigma,alpha,h/5,h,h*5,'composite');
[tau,R]=materncov(2,N,sigma,alpha,h/5,h,h*5,'composite');
[f,Spp2,Snn2]=blurspec(2,R,'aliased');

%[tau,R]=materncov(2,N,sigma,alpha,h/5,h*5,h,'composite');
%[f,Spp3,Snn3]=blurspec(2,R,'aliased');

%z=maternoise(N,sigma,alpha,h/5-1i*h*5,h,'composite');

%fact=2*sigma.^2.*(h/5).*((h*5).^2+h.^2).^alpha;
%2*sigma.^2.*H.*(nu.^2+G.^2).^alpha;

b1=aresame(Spp2./maxmax(Spp),Spp./maxmax(Spp),2e-2);
b2=aresame(Snn2./maxmax(Snn),Snn./maxmax(Snn),2e-2);

%this test is currently not working ... Spp2 is off by a factor of sqrt(3),
%related to the fact that R(1) is not unity in the above code after calling mspec
%reporttest('MATERNCOV Fourier transforms to MATERNSPEC for case composite Matern with negligible aliasing, even N',b1&&b2)

N=1000;
[tau,R]=materncov(1,N,sigma,-1/2,h,alpha*10,'extended');
[f,Spp,Snn]=maternspec(1,N,sigma,-1/2,h,alpha*10,'extended');
[f,Spp2,Snn2]=blurspec(1,R,'aliased');

b1=aresame(Spp2./maxmax(Spp),Spp./maxmax(Spp),1e-4);
b2=aresame(Snn2./maxmax(Snn),Snn./maxmax(Snn),1e-4);

reporttest('MATERNCOV Fourier transforms to MATERNSPEC for case with negligible aliasing, even N, exponential',b1&&b2)

N=1000;
dt=7;
[tau,R]=materncov(dt,N,sigma,-1/2,h/dt,(alpha*10*dt),'extended');
[f,Spp,Snn]=maternspec(dt,N,sigma,-1/2,h/dt,(alpha*10*dt),'extended');
[f,Spp2,Snn2]=blurspec(dt,R,'aliased');

b1=aresame(Spp2./maxmax(Spp),Spp./maxmax(Spp),1e-4);
b2=aresame(Snn2./maxmax(Snn),Snn./maxmax(Snn),1e-4);

reporttest('MATERNCOV Fourier transforms to MATERNSPEC ..., even N, exponential, non-unit sample rate',b1&&b2)

[tau,R]=materncov(1,1000,7,1.5,1/10,20,'extended');
[f,Spp,Snn]=maternspec(1,1000,7,1.5,1/10,20,'extended');
[f,Spp2,Snn2]=blurspec(1,R,'aliased');

b1=aresame(Spp2./maxmax(Spp),Spp./maxmax(Spp),1e-10);
b2=aresame(Snn2./maxmax(Snn),Snn./maxmax(Snn),1e-10);

reporttest('MATERNCOV Fourier transforms to MATERNSPEC for case with negligible aliasing, even N, extended Matern',b1&&b2)

dt=1.3;
[tau,R]=materncov(dt,1000,7,1.5,1/10/dt,20/dt,'extended');
[f,Spp,Snn]=maternspec(dt,1000,7,1.5,1/10/dt,20/dt,'extended');
%[f,Spp2,Snn2]=blurspec(dt,R);
[f,Spp2,Snn2]=blurspec(dt,R,'aliased');

b1=aresame(Spp2./maxmax(Spp),Spp./maxmax(Spp),1e-10);
b2=aresame(Snn2./maxmax(Snn),Snn./maxmax(Snn),1e-10);

reporttest('MATERNCOV Fourier transforms to MATERNSPEC ..., even N, extended Matern, non-unit sample rate',b1&&b2)

N=1000-1;
[tau,R]=materncov(1,N,sigma,alpha,h,2*h,'extended');
[f,Spp,Snn]=maternspec(1,N,sigma,alpha,h,2*h,'extended');
[f,Spp2,Snn2]=blurspec(1,R,'aliased');

b1=aresame(Spp2./maxmax(Spp),Spp./maxmax(Spp),1e-4);
b2=aresame(Snn2./maxmax(Spp),Snn./maxmax(Spp),1e-4);

reporttest('MATERNCOV Fourier transforms to MATERNSPEC for case with negligible aliasing, frequency shift case, odd N',b1&&b2)

N=1000;
dt=7;
lambda=h/dt;
mu=(alpha*10*dt);
[tau,R]=materncov(dt,N,sigma,-1/2,lambda,mu,'extended');
tnorm=sqrt(tau.^2+mu.^2);
fact=besselk(1,lambda.*mu);%See notes for this form
R1=frac(1,fact).*frac(1,tnorm./mu).*besselk(1,lambda.*tnorm);
reporttest('MATERNCOV ALPHA = -1/2 case matches expected form',aresame(R,R1,1e-10))

toc
% [tau,R]=materncov(1,N,sigma,alpha,h,'reverse');
% [f,Spp,Snn]=maternspec(1,N,sigma,alpha,h,'reverse');
% R(1)=R(1)/2;
% S=2*real(fft(R));
% Spp2=S(1:length(Spp));
% S=flipud(S);
% Snn2=[S(end);S(1:length(Snn)-1)];
%
% b1=aresame(Spp2./maxmax(Spp),Spp./maxmax(Spp),1e-4);
% b2=aresame(Snn2./maxmax(Snn),Snn./maxmax(Snn),1e-4);
%
% reporttest('MATERNCOV Fourier transforms to MATERNSPEC for case with negligible aliasing, even N, reverse Matern',b1&&b2)

%Playing around with asymptotic results... which don't work for alpha > 3/2
% sigma=7;
% alpha=[0.6 3/4 1.001 1.25 3/2-0.001 1.75];% 1.75 2 3 4 10];
% %alpha=[1.75 2 3 4 10];% 1.75 2 3 4 10];
% lambda=1/10000;
% N=10000;
% [tau,R]=materncov(N,sigma,alpha,lambda);
% Rapprox=zeros(size(R));
% for i=1:length(alpha)
%     %numer=(lambda.*tau).^(2*alpha(i)-1);
%     %denom=2*materncfun(alpha(i)).*cos(pi.*alpha(i)).*gamma(2*alpha(i));
%     numer1=-(lambda.*tau/2).^(2*alpha(i)-1)*gamma(3/2-alpha(i));
%     denom1=gamma(1/2+alpha(i));
%     numer2=(lambda.*tau/2).^2*gamma(3/2-alpha(i));
%     denom2=gamma(5/2-alpha(i));
%     Rapprox1(:,i)=sigma.^2.*(1+frac(numer,denom));
%     %Rapprox2(:,i)=sigma.^2.*(1+frac(numer1,denom2)+frac(numer2,denom2));
% end
% Rapprox(Rapprox<0)=nan;
% Rapprox(Rapprox>2*sigma.^2)=nan;

% alpha=1.5;
% h=0.1;
% sigma=1;
% nu=0.1;
%
% N=1000;
% figure
% [tau,R]=materncov(N,sigma,1,h-1i*nu);
% [f,Spp,Snn]=blurspec(R);
% plot(-f,Snn),hold on,plot(f,Spp)
% [tau,R]=materncov(N,sigma,1,h+1i*nu);
% [f,Spp,Snn]=blurspec(R);
% plot(-f,Snn),hold on,plot(f,Spp)
%
% [f,Spp,Snn]=maternspec(N,sigma,alpha,h);

% N=1000;
% r1=[0:.00001:0.00001]';
% alpha=[1/2:0.1:4];
% R=materncov(r1,1,alpha(:),1/100);
% dr=(R(2,:)-R(1,:))./(r1(2)-r1(1));
% dr2=-frac(gamma(abs(alpha-3/2)),gamma(alpha-1/2)).*frac(1/100,2);
% %dr2(alpha==1)=-1/100;
%
% plot(alpha,-dr,'r*'),hold on,plot(alpha,-dr2,'o')
%
% dxi=-frac(r,2).*frac(gamma(alpha-1),gamma(alpha)).*maternfun(alpha-1,r);
% reporttest('MATERNFUN iterative expression for derivative matches direct expression',aresame(xi1(2:end,:),dxi(2:end,:),1e-10))
%

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function qwarning(s)
% warning with option to stop

disp('To interrupt press Ctr-C');
warning(s);
pause;

author: mtex-toolbox

repo: mtex

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


text:
function [wc,w0,a0,ak,bk,c0,ck]=get_harmonics(y,pas)
%[wc,w0,a0,ak,bk,c0,ck]=get_harmonics(y,pas)
% given a signal x(t),y=x(t)
% the function get_harmoniques returns the pulse w0(rd/s) of the signal y(t)
% the coefficients ak (a1,a2,...)  and bk (b1,b2,...)  of the trigonomical
% representation  of the  signal y(t):
% y(t)= a0/2+ a1.cos(w0t)+b1.sin(w0t)+a2.cos(2w0t) +b2*sin(t*w0)+...
% it also returns the coefficients c0 and ck of the complex representation
% y(t)= ...+c(-1).exp(-jw0t)+c(0)+c(1).exp(jw0t)+c(2).exp(2jw0t)+...
% wc is the interval pulse;
% to plot the amplitude spectrum: stem(wc,abs(ck))
%[wc,w0,a0,ak,bk,c0,ck]=get_harmoniques(y,pas)
% be a signal x(t).
% For more theorical details look at www.azzimalik.com
N=length(y);
Y=fft(y);
if size(Y,2)==1;
    Y=Y';
end
w=2*(0:(N-1))*pi/(N*pas);
w0=w(2);
nc2=fix(N/2);
if nc2~=N/2;
    nc3=nc2+1;
    wc=[fliplr(-w(2:nc2+1)) w(1:nc2+1)];
    Yc=[Y(nc2+2:N) Y(1:nc2+1)];Yc=Yc/N;
else
    nc3=nc2;
    wc=[fliplr(-w(2:nc2+1)) w(1:nc2+1)];
    Yc=[Y(nc2+1:N) Y(1:nc2+1)];Yc=Yc/N;
end
%-------------------------------------------------------------
N1=length(Yc);
for l=2:nc2+1
    ak(l-1)=real(Y(l)+Y(N-l+2))/N;
    bk(l-1)=-imag(Y(l)-Y(N-l+2))/N;
    ckp(l-1)=(ak(l-1)-j*bk(l-1))/2;
    ckm(l-1)=(ak(l-1)+j*bk(l-1))/2;
end
a0=2*Y(1)/N;
c0=a0/2;ck=[fliplr(ckm) c0 ckp];

author: Sable

repo: mcbench-benchmarks


text:
function [ T1temp MP2RAGEcorrected] = T1B1correctpackageTFL( B1img,MP2RAGEimg,T1,MP2RAGE,brain,varargin)
% usage
%
% [  T1corr MP2RAGEcorr] = T1B1correctpackage(B1,MP2RAGEimg,T1,MP2RAGE,brain,varargin)
%
% B1 and MP2RAGEimg (and T1) are the nii structures resulting from loading
% the MP2RAGE (MP2RAGEimg, T1) and the result of some B1 mapping technique
%  with load_nii or load_untouch_nii
%
% the variable B1 is compulsory
%
% Only MP2RAGEimg or the T1 have to be loaded (I usually use the MP2RAGEimg)
%
% the variables that are not loaded can be simply left empty []
%
% MP2RAGE  variable contains all the relevant sequence
% information as delailed below
%
% B1 will be given in relative units => 1 if it was correct; values can varie
% from 0-2
%
%     MP2RAGE.B0=7;           % in Tesla
%     MP2RAGE.TR=6;           % MP2RAGE TR in seconds
%     MP2RAGE.TRFLASH=6.7e-3; % TR of the GRE readout
%     MP2RAGE.TIs=[800e-3 2700e-3];% inversion times - time between middle of refocusing pulse and excitatoin of the k-space center encoding
%     MP2RAGE.NZslices=[40 80];% Slices Per Slab * [PartialFourierInSlice-0.5  0.5]
%     MP2RAGE.FlipDegrees=[4 5];% Flip angle of the two readouts in degrees
%
% brain can be an image in the same space as the MP2RAGE and the
% Sa2RAGE that has zeros where there is no need to do any T1/B1 calculation
% (can be a binary mask or not). if left empty the calculation is done
% everywhere
%
% additionally the inversion efficiency of the adiabatic inversion can be
% set as a last optional variable. Ideally it should be 1.
% In the first implementation of the MP2RAGE the inversino efficiency was
% measured to be ~0.96
%
% outputs are:
%  T1corr - T1map corrected for B1 bias
%  MP2RAGEcorr - MP2RAGE image corrected for B1 bias
%
%
% please cite:
% Marques, J.P., Gruetter, R., 2013. New Developments and Applications of the MP2RAGE Sequence - Focusing the Contrast and High Spatial Resolution R1 Mapping. PLoS ONE 8. doi:10.1371/journal.pone.0069294
% Marques, J.P., Kober, T., Krueger, G., van der Zwaag, W., Van de Moortele, P.-F., Gruetter, R., 2010a. MP2RAGE, a self bias-field corrected sequence for improved segmentation and T1-mapping at high field. NeuroImage 49, 1271�1281. doi:10.1016/j.neuroimage.2009.10.002
%

% Check B1 range
b1med = median(prctile(B1img(:),90));

if b1med > 5

     warning(sprintf(['=============== mp2rage::b1correction ==========='...
        '\n B1 data is not in [0-2] range. B1map magnitude will be scaled down.'  ...
        '\n ===========================================================' ...
        ]));

    if b1med > 10 &&  b1med < 500

        B1img = double(B1img)./100;

    elseif b1med>500 && b1med<1500

        B1img = double(B1img)./1000;
    end
end

b1med = median(prctile(B1img(:),90));

if ~(b1med > 0.5) && ~(b1med < 1.5)

    warning(sprintf(['=============== mp2rage::b1correction ==========='...
        '\n B1 data may not be in the required [0-2] range'  ...
        '\n ===========================================================' ...
        ]));

end

if nargin==6

    invEFF=varargin{1};

else

    invEFF=0.96;

end

if isempty(brain)

    if isempty(MP2RAGEimg)

        brain=T1;

        brain.img=ones(size(brain.img));

    else

        brain=MP2RAGEimg;

        brain.img=ones(size(brain.img));

    end;

end;

B1_vector=0.005:0.05:1.9;

T1_vector=0.5:0.05:5.2;

[MP2RAGE.Intensity MP2RAGE.T1vector ]=MP2RAGE_lookuptable(2,MP2RAGE.TR,MP2RAGE.TIs,MP2RAGE.FlipDegrees,MP2RAGE.NZslices,MP2RAGE.TRFLASH,'normal',invEFF);

if isempty(MP2RAGEimg)

    T1.img=double(T1.img)/1000;

    MP2RAGEimg.img=reshape(interp1(MP2RAGE.T1vector,MP2RAGE.Intensity,T1.img(:)),size(B1img));

    MP2RAGEimg.img(isnan(MP2RAGEimg.img))=-0.5;

else

    MP2RAGEimg.img=double(MP2RAGEimg.img)/4095-0.5;

end;

%% now the fun starts

% creates a lookup table of MP2RAGE intensities as a function of B1 and T1

k=0;

for b1val=B1_vector

    k=k+1;

    [Intensity T1vector ]=MP2RAGE_lookuptable(2,MP2RAGE.TR,MP2RAGE.TIs,b1val*MP2RAGE.FlipDegrees,MP2RAGE.NZslices,MP2RAGE.TRFLASH,'normal',invEFF);

    MP2RAGEmatrix(k,:)=interp1(T1vector,Intensity,T1_vector);

end;

k=0;

%% make the matrix  MP2RAGEMatrix into T1_matrix(B1,ratio)

npoints=40;

MP2RAGE_vector=linspace(-0.5,0.5,npoints);

k=0;

for b1val=B1_vector

    k=k+1;

    try

        T1matrix(k,:)=interp1(MP2RAGEmatrix(k,:),T1_vector,MP2RAGE_vector,'pchirp');

    catch

        temp=MP2RAGEmatrix(k,:);

        temp(isnan(temp))=linspace(-0.5-eps,-1,length(find(isnan(temp)==1)));

        temp=interp1(temp,T1_vector,MP2RAGE_vector);

        T1matrix(k,:)=temp;

    end;

end;

%% correcting the estimates of T1 and B1 iteratively

T1temp=MP2RAGEimg;

brain.img(B1img==0)=0;

brain.img(MP2RAGEimg.img==0)=0;

T1temp.img(brain.img==0)=0;

T1temp.img(brain.img==1)=0;

B1img(brain.img==0)=0;

temp=squeeze(T1temp.img(:,end/2,:));

T1temp.img(brain.img~=0)=interp2(MP2RAGE_vector,B1_vector,T1matrix,MP2RAGEimg.img(brain.img~=0),B1img(brain.img~=0));

T1temp.img(isnan(T1temp.img))=4;

temp2=squeeze(T1temp.img(:,end/2,:));

%% creates an MP2RAGEcorrected image and puts both the B1 and T1 in the ms scale

[MP2RAGE.Intensity MP2RAGE.T1vector ]=MP2RAGE_lookuptable(2,MP2RAGE.TR,MP2RAGE.TIs,MP2RAGE.FlipDegrees,MP2RAGE.NZslices,MP2RAGE.TRFLASH,'normal',invEFF);

MP2RAGEcorrected=MP2RAGEimg;

MP2RAGEcorrected.img=reshape(interp1(MP2RAGE.T1vector,MP2RAGE.Intensity,T1temp.img(:)),size(T1temp.img));

MP2RAGEcorrected.img(isnan(MP2RAGEcorrected.img))=-0.5;

MP2RAGEcorrected.img=round(4095*(MP2RAGEcorrected.img+0.5));

%T1temp.img=(T1temp.img)*1000;

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function res = SumKL(p,K,L)
%SUMKL        Summation 'as if' computed in K-fold precision and stored in L results
%
%   res = SumKL(p,K,L)
%
%On return, sum(res) approximates sum(p) with accuracy as if computed
%  in K-fold precision, where res comprises of L elements.
%  Default for L is 1.
%
%Implements algorithm SumKL from
%  S.M. Rump: Inversion of extremely ill-conditioned matrices in floating-point,
%      Japan J. Indust. Appl. Math. (JJIAM), 26:249-277, 2009.
%
%Reference implementation! Slow due to interpretation!
%

% written  06/23/08     S.M. Rump
% modified 05/09/09     S.M. Rump  rounding to nearest
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if nargin==2
    L = 1;
  end

  n = length(p);
  for i=1:K-L
    p = VecSum(p);
  end
  res = zeros(1,L);
  for k=0:L-2
    p(1:n-k) = VecSum(p(1:n-k));
    res(k+1) = p(n-k);
  end
  res(L) = sum(p(1:n-L+1));

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function a = lietzke_inverse ( n )

%*****************************************************************************80
%
%% LIETZKE_INVERSE returns the inverse of the LIETZKE matrix.
%
%  Example:
%
%    N = 5
%
%   0.5833   -0.5000         0         0    0.0833
%  -0.5000    1.0000   -0.5000         0         0
%        0   -0.5000    1.0000   -0.5000         0
%        0         0   -0.5000    1.0000   -0.5000
%   0.0833         0         0   -0.5000    0.5833
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of rows and columns
%    of the matrix.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  a(1,1) = ( n + 2 ) / ( 2 * n + 2 );
  for i = 2 : n - 1
    a(i,i) = 1.0;
  end
  a(n,n) = ( n + 2 ) / ( 2 * n + 2 );

  if ( n == 2 )

    for i = 1 : n - 1
      a(i,i+1) = - 1.0 / 3.0;
    end

    for i = 2 : n
      a(i,i-1) = - 1.0 / 3.0;
    end

  else

    for i = 1 : n - 1
      a(i,i+1) = - 0.5;
    end

    for i = 2 : n
      a(i,i-1) = - 0.5;
    end

  end

  a(1,n) = 1.0 / ( 2 * n + 2 );
  a(n,1) = 1.0 / ( 2 * n + 2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ehat,v_e,etahat,v_eta] = VBA_getNoise(posterior,out)
% returns the Laplace approximation to the innovations posterior density
% function [ehat,v_e,etahat,v_eta] = VBA_getNoise(posterior)
% This may be useful if one is interested in recovering, e.g., the state
% noise that enters and perturbs the system. In particular, posterior
% covariances are attached to state noise, which means that statistical
% inference can be performed in the usual way...
% IN:
%   - posterior/out: output structures of VBA_NLStateSpaceModel.m
% OUT:
%   - ehat/etahat: the 1st-order moment of the posterior density on the
%   measurement (resp. state) noise
%   - v_e/v_eta: the second-order moment of the posterior density on the
%   measurement (resp. state) noise

dim = out.dim;
options = out.options;
u = out.u;
y = out.y;

ehat = zeros(dim.p,dim.n_t);
v_e = cell(1,dim.n_t);
etahat = zeros(dim.n,dim.n_t);
v_eta = cell(1,dim.n_t);

% initial condition
if dim.n > 0
    [fx,dfdx,dfdp] = VBA_evalFun('f',posterior.muX0,posterior.muTheta,u(:,1),options,dim,1);
    etahat(:,1) = posterior.muX(:,1) - fx;
    if isinf(posterior.a_alpha) && isequal(posterior.b_alpha,0)
        v_eta{1} = zeros(dim.n,dim.n);
    else
        v_eta{1} = posterior.SigmaX.current{1} ...
            + dfdx'*posterior.SigmaX0*dfdx;
        if dim.n_theta > 0
            v_eta{1} = v_eta{1} + dfdp'*posterior.SigmaTheta*dfdp;
        end
    end
end
[gx,dgdx,dgdp] = VBA_evalFun('g',posterior.muX(:,1),posterior.muPhi,u(:,1),options,dim,1);
ehat(:,1) = y(:,1) - gx;
v_e{1} = zeros(dim.p,dim.p);
if dim.n > 0
    v_e{1} = dgdx'*posterior.SigmaX.current{1}*dgdx;
end
if dim.n_phi > 0
    v_e{1} = v_e{1} + dgdp'*posterior.SigmaPhi*dgdp;
end

% loop over time samples
for t = 2:dim.n_t
    if dim.n > 0
        [fx,dfdx,dfdp] = VBA_evalFun('f',posterior.muX(:,t-1),posterior.muTheta,u(:,t),options,dim,t);
        etahat(:,t) = posterior.muX(:,t) - fx;
        if isinf(posterior.a_alpha) && isequal(posterior.b_alpha,0)
            v_eta{t} = zeros(dim.n,dim.n);
        else
            P = [-dfdx',eye(dim.n)];
            jointCov = ...
                [ posterior.SigmaX.current{t}  posterior.SigmaX.inter{t-1}'
                posterior.SigmaX.inter{t-1}  posterior.SigmaX.current{t-1} ];
            v_eta{t} = P*jointCov*P';
            if dim.n_theta > 0
                v_eta{t} = v_eta{t} + dfdp'*posterior.SigmaTheta*dfdp;
            end
        end
    end
    [gx,dgdx,dgdp] = VBA_evalFun('g',posterior.muX(:,t),posterior.muPhi,u(:,t),options,dim,t);
    ehat(:,t) = y(:,t) - gx;
    v_e{t} = zeros(dim.p,dim.p);
    if dim.n > 0
        v_e{t} = dgdx'*posterior.SigmaX.current{t}*dgdx;
    end
    if dim.n_phi > 0
        v_e{t} = v_e{t} + dgdp'*posterior.SigmaPhi*dgdp;
    end
end

author: MBB-team

repo: VBA-toolbox


text:
function B = tomatrix(A)
    % Flatten the symbolic expression

    % Convert inputs to SymExpression
    % A = SymExpression(A);

    % construct the operation string
    sstr = ['ToMatrixForm[' A.s ']'];

    % create a new object with the evaluated string
    B = SymExpression(sstr);

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function [z, mu, s] = zscore(x, varargin)

% [Z, MU, S] = ZSCORE(X, NORMALIZEFLAG, DIM, FLAG) computes the zscore, across all cells in x along
% the dimension dim, normalising by the total number of samples
%
% X should be an linear cell-array of matrices for which the size in at
% least one of the dimensions should be the same for all cells. If flag==1, the mean will
% be subtracted first (default behavior, but to save time on already demeaned data, it
% can be set to 0). MU and S are vectors containing the mean and standard deviations used
% for zscoring

if numel(varargin)==0, varargin{1} = []; end
if numel(varargin)==1, varargin{2} = []; end
if numel(varargin)==2, varargin{3} = 1;  end

if varargin{3}
  mu = nanmean(x, varargin{2});
  x  = cellvecadd(x, -mu);
end
s = nanstd(x, varargin{1:end});
z = cellvecmult(x, 1./s);

if ~varargin{3}
  mu = zeros(size(s))+nan;
end

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [] = aps_support_plot(technique)
% plotting the support information used for the APS correction methods
% technique flag
% when:
% 1 powerlaw
% 2 ERA-I
% 3 Meris
%
%     Copyright (C) 2015  Bekaert David - University of Leeds
%     Email: eedpsb@leeds.ac.uk or davidbekaert.com
%
%     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 Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
% By Bekaert David
%

fontsize = 15;
load tca_support

%% power- law
if technique==1
    if exist('aps_p','dir')~=7
        mkdir('aps_p')
    end
    if strcmpi(getparm_aps('powerlaw_ridge_constraint'),'y')

        deminfo = powerlaw_ridges.deminfo;
        InSAR_convexhull = powerlaw_ridges.InSAR_convexhull;
        mountain_ridge = powerlaw_ridges.mountain_ridge;
        hard_ridge_flag = powerlaw_ridges.hard_ridge_flag;

        % plotting hte powerlaw ridges
        fprintf('Plotting power-law ridges \n')

        h1 = figure('name','Power-law ridges','position', [ 200   243   610   603]);
        plot( mean([deminfo.xmin deminfo.xmax]),mean([deminfo.ymax deminfo.ymin]),'g-','linewidth',2)
        hold on
        plot( mean([deminfo.xmin deminfo.xmax]),mean([deminfo.ymax deminfo.ymin]),'r-','linewidth',2)
        hold on
        plot( mean([deminfo.xmin deminfo.xmax]),mean([deminfo.ymax deminfo.ymin]),'r--','linewidth',2)

        imagesc([deminfo.xmin deminfo.xmax],[deminfo.ymax deminfo.ymin],deminfo.dem_smooth)
        cc = colorbar;
        view(0,90)
        axis equal
        axis tight
        colormap(flipud(gray))
        axis xy
        box on

        hold on
        plot(InSAR_convexhull(:,1),InSAR_convexhull(:,2),'g-','linewidth',2)

        % plotting the line on the figure
        for counter=1:length(mountain_ridge)
            if hard_ridge_flag(counter)==1
                plot(mountain_ridge{counter}(:,1),mountain_ridge{counter}(:,2),'r-','linewidth',2)
            else
                plot(mountain_ridge{counter}(:,1),mountain_ridge{counter}(:,2),'r--','linewidth',2)
            end
        end
        box on
        xlabel(cc,'[m]','fontsize',fontsize)
        ylabel(cc,'Topography','fontsize',fontsize)
        legend('Location','NorthOutside','InSAR region','Hard ridge','Subregion boundary')
        title('Ridge definition','fontsize',fontsize)
        set(gca,'fontsize',fontsize)
        print(h1,'-dpng',['aps_p/mountain_ridges.png'])
        print(h1,'-depsc',['aps_p/mountain_ridges.eps'])

    else
       % plotting the power-law windows
        fprintf('Plotting power-law windows \n')

        % see if the real dem is available otherzie use the poitn heights
        dem_flag=1;
        if exist('powerlaw_ridges','var')==1
            deminfo = powerlaw_ridges.deminfo;
        elseif exist('era','var')==1
            deminfo = era.deminfo;
        else
            hgt_file = getparm_aps('hgt_matfile');
            ll_matfile = getparm_aps('ll_matfile');

            hgt =  load(hgt_file);
            ll = load(ll_matfile);
            if strcmpi(getparm_aps('stamps_processed'),'y')
                hgt = hgt.hgt;
                ll = ll.lonlat;
                dem_flag=0;
            end
        end
        % making sure the dem and the other data have the same origin in
        % longitude
        if dem_flag== 1
            if deminfo.xmax-max(InSAR_convexhull(:,1))>100
                deminfo.xmax = deminfo.xmax-360;
                deminfo.xmin = deminfo.xmin-360;
            end
        end

        % laoding of the data
        window_box_ll= powerlaw_windows.window_box_ll;
        window_box_center_ll = powerlaw_windows.window_box_center_ll;

        % plotting the results
        h1 = figure('name','Power-law windows','position', [ 200   243   610   603]);
        % for the legend
        plot(window_box_ll{round(0.5*(length(window_box_ll)))}(:,1),window_box_ll{round(0.5*(length(window_box_ll)))}(:,2),'r-','linewidth',2)
        hold on
        plot(window_box_center_ll(1,1),window_box_center_ll(1,2),'wo','markeredgecolor','k','markerfacecolor','w','markersize',15)
        hold on
        plot(InSAR_convexhull(:,1),InSAR_convexhull(:,2),'g-','linewidth',2)
        hold on

        if dem_flag==1
             % plot the topogrpahy
            imagesc([deminfo.xmin deminfo.xmax],[deminfo.ymax deminfo.ymin],deminfo.dem)
            view(0,90)
            axis xy
        else
            scatter3(ll(:,1),ll(:,2),hgt,3,hgt,'filled');
            view(0,90)
        end
        cc=colorbar;
        colormap(flipud(gray))
        xlabel(cc,'[m]','fontsize',fontsize)
        ylabel(cc,'Topography','fontsize',fontsize)
        axis equal
        axis tight

        hold on
        plot3(window_box_ll{round(0.5*(length(window_box_ll)))}(:,1),window_box_ll{round(0.5*(length(window_box_ll)))}(:,2),99999.*ones(size(window_box_ll{1}(:,2))),'r-','linewidth',2)
        hold on
        plot3(window_box_center_ll(:,1),window_box_center_ll(:,2),10000000.*ones(size(window_box_center_ll(:,2))),'wo','markeredgecolor','k','markerfacecolor','w','markersize',15)
        hold on
        plot3(InSAR_convexhull(:,1),InSAR_convexhull(:,2),99999.*ones(size(InSAR_convexhull(:,2))),'g-','linewidth',2)

        legend('Location','NorthOutside','Window example','Window centers','InSAR region')
        title('Window definition','fontsize',fontsize)
        set(gca,'fontsize',fontsize)

        grid off
        box on

        print(h1,'-dpng',['aps_p/window_locations.png'])
        print(h1,'-depsc',['aps_p/window_locations.eps'])

    end

%% ERA interim
elseif technique==2
    if exist('aps_e','dir')~=7
        mkdir('aps_e')
    end

    % loading the ERA data
    deminfo = era.deminfo;
    era_lonlat = era.era_lonlat;

    % check if the InSAR outline exist, if so plot it too
    if exist('InSAR_convexhull','var')==1;
        insar_flag =1;
    else
        insar_flag=0;
    end
    if insar_flag== 1
        if deminfo.xmax-max(InSAR_convexhull(:,1))>100
            InSAR_convexhull(:,1) = InSAR_convexhull(:,1)+360;
        end
    end

    % plotting the figure
    hfig = figure('name','ERA points and InSAR region','position', [ 200   243   610   603]);
    % for the legend generation
    plot(mean([deminfo.xmin deminfo.xmax]),mean([deminfo.ymax deminfo.ymin]),'wo','markeredgecolor','k','markerfacecolor','w','markersize',15)
    hold on
    % check if the insar convex hull can be plotted
    if insar_flag==1
        plot(mean([deminfo.xmin deminfo.xmax]),mean([deminfo.ymax deminfo.ymin]),'g-','linewidth',2)
    end
    % plot the topogrpahy
    imagesc([deminfo.xmin deminfo.xmax],[deminfo.ymax deminfo.ymin],deminfo.dem)
    cc=colorbar;
    view(0,90)
    colormap(flipud(gray))
    axis xy
    xlabel(cc,'[m]','fontsize',fontsize)
    ylabel(cc,'Topography','fontsize',fontsize)
    hold on
    plot(era_lonlat(:,1),era_lonlat(:,2),'wo','markeredgecolor','k','markerfacecolor','w','markersize',15)
    hold on
    if insar_flag==1
        plot(InSAR_convexhull(:,1),InSAR_convexhull(:,2),'g-','linewidth',2)
        legend('location','northoutside','Used ERA-I locations','InSAR region')
    else
        legend('location','northoutside','Used ERA-I locations')
    end
    title('ERA points distribution','fontsize',fontsize)
    set(gca,'fontsize',fontsize)
    axis equal
    axis tight

    print(hfig,'-dpng',['aps_e/era_datapoints.png'])
    print(hfig,'-depsc',['aps_e/era_datapoints.eps'])

elseif technique==3 % MERIS
    load psver
    ps = load(['ps' num2str(psver) '.mat']);
    stamps_processed = getparm_aps('stamps_processed');
    small_baseline_flag  = getparm('small_baseline_flag');
    if strcmpi(small_baseline_flag,'y')
        load('tca_sb2.mat','ph_tropo_meris');
        % for baselines you need ps of PS directory
        bperp = load(['..' filesep 'ps1.mat'],'bperp');
        bperp = bperp.bperp;
        bperp = bperp(ps.ifgday_ix);
        dates=ps.ifgday;
        ix_ifgs_keep = 1:ps.n_ifg;
    else
        load('tca2.mat','ph_tropo_meris');
        dates = [ps.day   repmat(ps.master_day,length(ps.day),1)];
        bperp = [ps.bperp zeros([length(ps.day) 1])];
        ix_ifgs_keep = 1:ps.n_ifg;
    end

    % plot the SB network and show for whcih connection there is meris data
    if strcmp(stamps_processed,'y')

        % remove that that have not been unwrapped before
        ix_dropped = getparm('drop_ifg');

        % interferograms with a meris correction
        temp = ix_ifgs_keep;
        temp(ix_dropped)=[];
        ifgs_good = find(sum(ph_tropo_meris(:,temp)~=0,1)~=0)

        if strcmp(small_baseline_flag,'n')
            ix_dropped = unique([ix_dropped ps.master_ix]);
        end
        ix_ifgs_keep(ix_dropped)=[];
        bperp(ix_dropped,:)=[];
        dates(ix_dropped,:)=[];

        % keep the original network for reference
        dates_all = dates;
        bperp_all = bperp;

        % plotting the network of the interferograms used in the RMSE computation
        h_baselineplot = figure('name','Processed network');
        % plotting all ifgs that were considered
        for ifgs_counter=1:size(bperp_all,1)
            hold on
            plot([dates_all(ifgs_counter,1)  dates_all(ifgs_counter,2)],  [bperp_all(ifgs_counter,1) bperp_all(ifgs_counter,2)] ,'k-','linewidth',1)
        end

        % plotting the network for which we have an APS correction
        for ifgs_counter=1:length(ifgs_good)
            ifgs_counter
            hold on
            plot([dates(ifgs_good(ifgs_counter),1)  dates(ifgs_good(ifgs_counter),2)],  [bperp(ifgs_good(ifgs_counter),1) bperp(ifgs_good(ifgs_counter),2)] ,'k-','linewidth',2)
            text(mean(dates(ifgs_good(ifgs_counter),:)),mean(bperp(ifgs_good(ifgs_counter),:))+20,num2str(ifgs_good(ifgs_counter)),'fontsize',fontsize-2,'backgroundColor',[1 1 1],'Margin',0.01)
        end
        hold on
        % [dates_unique,ix_unique] = unique(dates(ifgs_good,:));
        % bperp_unique = bperp(ifgs_good,:);
        % bperp_unique = bperp_unique(ix_unique);
        [dates_unique,ix_unique] = unique(dates_all);
        bperp_unique = bperp_all;
        bperp_unique = bperp_unique(ix_unique);
        plot(dates_unique,bperp_unique,'ko','markerfacecolor','r','markersize',7)
        clear dates_unique ix_unique bperp_unique
        hold on
        % set(gca,'XTick',dates_num)
        datetick('x','mmm yy')
        set(gca,'fontsize',fontsize)
        box on
        ylabel('Bperp [m]','fontsize',fontsize)
    end

end

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
function [EC,ec,degij] = edge_nei_overlap_bd(CIJ)
% EDGE_NEI_OVERLAP_BD     Overlap amongst neighbors of two adjacent nodes
%
%   [EC,ec,degij] = edge_nei_bd(CIJ);
%
%   This function determines the neighbors of two nodes that are linked by
%   an edge, and then computes their overlap.  Connection matrix must be
%   binary and directed.  Entries of 'EC' that are 'inf' indicate that no
%   edge is present.  Entries of 'EC' that are 0 denote "local bridges",
%   i.e. edges that link completely non-overlapping neighborhoods.  Low
%   values of EC indicate edges that are "weak ties".
%
%   If CIJ is weighted, the weights are ignored. Neighbors of a node can be
%   linked by incoming, outgoing, or reciprocal connections.
%
%   Inputs:     CIJ,      directed (binary/weighted) connection matrix
%
%   Outputs:    EC,     edge neighborhood overlap matrix
%               ec,     edge neighborhood overlap per edge, in vector format
%               degij,  degrees of node pairs connected by each edge
%
%   Reference:
%
%       Easley and Kleinberg (2010) Networks, Crowds, and Markets.
%           Cambridge University Press, Chapter 3
%
%   Olaf Sporns, Indiana University, 2012

[ik,jk,ck] = find(CIJ);
lel = length(ck);
N = size(CIJ,1);

[~,~,deg] = degrees_dir(CIJ);

ec = zeros(1,lel);
degij = zeros(2,lel);
for e=1:lel
    neiik = setdiff(union(find(CIJ(ik(e),:)),find(CIJ(:,ik(e))')),[ik(e) jk(e)]);
    neijk = setdiff(union(find(CIJ(jk(e),:)),find(CIJ(:,jk(e))')),[ik(e) jk(e)]);
    ec(e) = length(intersect(neiik,neijk))/length(union(neiik,neijk));
    degij(:,e) = [deg(ik(e)) deg(jk(e))];
end;

ff = find(CIJ);
EC = 1./zeros(N);
EC(ff) = ec;                            %#ok<FNDSB>

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function blas1_z_test06 ( )

%*****************************************************************************80
%
%% TEST06 tests ZAXPY.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2006
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  x = [ ...
      2.0 -1.0 * i, ...
     -4.0 -2.0 * i, ...
      3.0  1.0 * i, ...
      2.0  2.0 * i, ...
     -1.0 -1.0 * i ];
  y = [ ...
     -1.0,  0.0 * i, ...
      0.0, -3.0 * i, ...
      4.0,  0.0 * i, ...
     -3.0,  4.0 * i, ...
     -2.0,  0.0 * i ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST06\n' );
  fprintf ( 1, '  ZAXPY adds a multiple of one complex vector to another.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  X = \n' );
  fprintf ( 1, '\n' );
  for j = 1 : n
    fprintf ( 1, '  %6d  %10f  %10f\n', j, real ( x(j) ), imag ( x(j) ) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Y = \n' );
  fprintf ( 1, '\n' );
  for j = 1 : n
    fprintf ( 1, '  %6d  %10f  %10f\n', j, real ( y(j) ), imag ( y(j) ) );
  end

  s = 0.50 - 1.00 * i;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The scalar multiplier is: ( %f  %f )\n', ...
    real ( s ), imag ( s ) );

  y = zaxpy ( n, s, x, 1, y, 1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  A * X + Y = \n' );
  fprintf ( 1, '\n' );
  for j = 1 : n
    fprintf ( 1, '  %6d  %10f  %10f\n', j, real ( y(j) ), imag ( y(j) ) );
  end

  return
end


text:
function isargchar(varargin)
%ISARGCHAR tests if the given arg is a char and returns an error otherwise
%
%   Usage: isargstruct(arg1,arg2,...)
%
%   Input options:
%       args        - list of args
%
%   ISARGCHAR(args) tests if all given args are a char and returns
%   an error otherwise.
%
%   see also: isargstruct

%*****************************************************************************
% Copyright (c) 2010-2012 Quality & Usability Lab                            *
%                         Deutsche Telekom Laboratories, TU Berlin           *
%                         Ernst-Reuter-Platz 7, 10587 Berlin, Germany        *
%                                                                            *
% This file is part of the Sound Field Synthesis-Toolbox (SFS).              *
%                                                                            *
% The SFS 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.                                                         *
%                                                                            *
% The SFS 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/>.            *
%                                                                            *
% The SFS is a toolbox for Matlab/Octave to  simulate and  investigate sound *
% field  synthesis  methods  like  wave  field  synthesis  or  higher  order *
% ambisonics.                                                                *
%                                                                            *
% http://dev.qu.tu-berlin.de/projects/sfs-toolbox       sfstoolbox@gmail.com *
%*****************************************************************************

% AUTHOR: Hagen Wierstorf
% $LastChangedDate: 2012-04-26 10:12:51 +0200 (Thu, 26 Apr 2012) $
% $LastChangedRevision: 710 $
% $LastChangedBy: wierstorf.hagen $

%% ===== Checking for struct =============================================
for ii = 1:nargin
    if ~ischar(varargin{ii})
        error('%s need to be a string.',inputname(ii));
    end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MRiLab auto generated file: DO NOT EDIT!     %
% Generated by MRiLab "DoWriteXML2m" Generator %
% MRiLab Version 1.3                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [rfAmp,rfPhase,rfFreq,rfCoil,GzAmp,GyAmp,GxAmp,ADC,Ext,uts,ts,flags]=PSD_SPGR3DMT
global VCtl
global VVar
CV1=1e-3;
CV10=0;
CV11=0;
CV12=0;
CV13=0;
CV14=0;
CV2=2e-3;
CV3=abs(rem(VVar.TRCount,2)*2-1);
CV4=18e-3;
CV5=0;
CV6=250e3;
CV7=0;
CV8=0;
CV9=0;
rfAmpAll=[];
rfPhaseAll=[];
rfFreqAll=[];
rfCoilAll=[];
GzAmpAll=[];
GyAmpAll=[];
GxAmpAll=[];
ADCAll=[];
ExtAll=[];
rfTimeAll=[];
GzTimeAll=[];
GyTimeAll=[];
GxTimeAll=[];
ADCTimeAll=[];
ExtTimeAll=[];
SEtAll=[];
uts=[];
ts=[];
flags=[];
if VCtl.PlotSeq == 1
rfAmp=[];
rfPhase=[];
rfFreq=[];
rfCoil=[];
GzAmp=[];
GyAmp=[];
GxAmp=[];
ADC=[];
Ext=[];
Freq=1;
Notes='regular TR section';
AttributeOpt={'on','off'};
Switch=AttributeOpt{1};
TREnd=Inf;
TRStart=1;
tE=VCtl.TR;
tS=0;
if VVar.TRCount<TRStart | VVar.TRCount>TREnd | mod(VVar.TRCount-TRStart,Freq)~=0 | strcmp(Switch,'off')
% do nothing
else
ts = [ts tS tE];
end
ts = [0 max(ts)-min(ts)];
return;
end
%==============Pulses 1==============
rfAmp=[];
rfPhase=[];
rfFreq=[];
rfCoil=[];
GzAmp=[];
GyAmp=[];
GxAmp=[];
ADC=[];
Ext=[];
rfTime=[];
GzTime=[];
GyTime=[];
GxTime=[];
ADCTime=[];
ExtTime=[];
Freq=1;
Notes='regular TR section';
AttributeOpt={'on','off'};
Switch=AttributeOpt{1};
TREnd=Inf;
TRStart=1;
tE=VCtl.TR;
tS=0;
if isempty(tS) | isempty(tE) | (tS>=tE)
error('SE setting is incorrect for Pulses 1!');
end
if VVar.TRCount<TRStart | VVar.TRCount>TREnd | mod(VVar.TRCount-TRStart,Freq)~=0 | strcmp(Switch,'off')
% do nothing
else
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{1};
AttributeOpt={'Non','Hamming','Hanning'};
p.Apod=AttributeOpt{2};
p.CoilID=1;
p.DupSpacing=0;
p.Duplicates=1;
p.FA=VCtl.FlipAng;
p.Notes='sinc rf pulse';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.TBP=4;
p.dt=20e-6;
p.rfFreq=0;
p.rfPhase=0;
p.tEnd=CV4+0.6e-3;
p.tStart=CV4+0.1e-3;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp1,rfPhase1,rfFreq1,rfCoil1,rfTime1]=rfSinc(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil1(1)
rfAmp=[rfAmp rfAmp1];
rfPhase=[rfPhase rfPhase1];
rfFreq=[rfFreq rfFreq1];
rfCoil=[rfCoil rfCoil1];
rfTime=[rfTime rfTime1];
end
else
rfAmp=[rfAmp rfAmp1];
rfPhase=[rfPhase rfPhase1];
rfFreq=[rfFreq rfFreq1];
rfCoil=[rfCoil rfCoil1];
rfTime=[rfTime rfTime1];
end
end
p=[];
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{2};
AttributeOpt={'Non','Hamming','Hanning'};
p.Apod=AttributeOpt{2};
p.CoilID=1;
p.DupSpacing=0;
p.Duplicates=1;
p.FA=0;
p.Notes='sinc rf pulse';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.TBP=4;
p.dt=0.2e-3;
p.rfFreq=0;
p.rfPhase=0;
p.tEnd=VCtl.TR;
p.tStart=CV4+0.7e-3;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp2,rfPhase2,rfFreq2,rfCoil2,rfTime2]=rfSinc(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil2(1)
rfAmp=[rfAmp rfAmp2];
rfPhase=[rfPhase rfPhase2];
rfFreq=[rfFreq rfFreq2];
rfCoil=[rfCoil rfCoil2];
rfTime=[rfTime rfTime2];
end
else
rfAmp=[rfAmp rfAmp2];
rfPhase=[rfPhase rfPhase2];
rfFreq=[rfFreq rfFreq2];
rfCoil=[rfCoil rfCoil2];
rfTime=[rfTime rfTime2];
end
end
p=[];
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{2};
p.CoilID=1;
p.DupSpacing=0;
p.Duplicates=1;
p.FA=CV5;
p.Notes='Fermi rf pulse';
p.PW=5e-3;
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.dt=100e-6;
p.rfFreq=CV6;
p.rfPhase=0;
p.tEnd=CV4;
p.tStart=0;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp3,rfPhase3,rfFreq3,rfCoil3,rfTime3]=rfFermi(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil3(1)
rfAmp=[rfAmp rfAmp3];
rfPhase=[rfPhase rfPhase3];
rfFreq=[rfFreq rfFreq3];
rfCoil=[rfCoil rfCoil3];
rfTime=[rfTime rfTime3];
end
else
rfAmp=[rfAmp rfAmp3];
rfPhase=[rfPhase rfPhase3];
rfFreq=[rfFreq rfFreq3];
rfCoil=[rfCoil rfCoil3];
rfTime=[rfTime rfTime3];
end
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Gz1Sign=1;
p.Gz2Sign=0;
p.Notes='cartesian phase';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.t1End=CV4+CV2;
p.t1Start=CV4+CV1;
p.t2End=VCtl.TE+CV2;
p.t2Start=VCtl.TE+CV1;
p.tRamp=100e-6;
if strcmp(p.Switch,'on')
[GzAmp1,GzTime1]=GzCartesian(p);
GzAmp=[GzAmp GzAmp1];
GzTime=[GzTime GzTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Gy1Sign=1;
p.Gy2Sign=0;
p.Notes='cartesian phase';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.t1End=CV4+CV2;
p.t1Start=CV4+CV1;
p.t2End=VCtl.TE+CV2;
p.t2Start=VCtl.TE+CV1;
p.tRamp=100e-6;
if strcmp(p.Switch,'on')
[GyAmp1,GyTime1]=GyCartesian(p);
GyAmp=[GyAmp GyAmp1];
GyTime=[GyTime GyTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Gx1Sign=-1;
p.Gx2Sign=1;
p.Gx3Sign=0;
p.Notes='cartesian frequency';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.t1Start=CV4+CV1;
p.t2Middle=VCtl.TE;
p.t3Start=VCtl.TE+CV1;
p.tRamp=100e-6;
if strcmp(p.Switch,'on')
[GxAmp1,GxTime1]=GxCartesian(p);
GxAmp=[GxAmp GxAmp1];
GxTime=[GxTime GxTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Notes='cartesian readout';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tMiddle=VCtl.TE;
if strcmp(p.Switch,'on')
[ADC1,ADCTime1]=ADCCartesian(p);
ADC=[ADC ADC1];
ADCTime=[ADCTime ADCTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=5;
p.Notes='calculate remaining scan time';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=0;
if strcmp(p.Switch,'on')
[Ext1,ExtTime1]=ExtBit(p);
Ext=[Ext Ext1];
ExtTime=[ExtTime ExtTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=1;
p.Notes='reset K space location';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=CV4+0.6e-3;
if strcmp(p.Switch,'on')
[Ext2,ExtTime2]=ExtBit(p);
Ext=[Ext Ext2];
ExtTime=[ExtTime ExtTime2];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=6;
p.Notes='dephase Mxy';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=VCtl.TR*(99/100);
if strcmp(p.Switch,'on')
[Ext3,ExtTime3]=ExtBit(p);
Ext=[Ext Ext3];
ExtTime=[ExtTime ExtTime3];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=6;
p.Notes='dephase Mxy';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=CV4+0.05e-3;
if strcmp(p.Switch,'on')
[Ext4,ExtTime4]=ExtBit(p);
Ext=[Ext Ext4];
ExtTime=[ExtTime ExtTime4];
end
p=[];
%--------------------
SEt=[tS tE];
rfAmp(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfPhase(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfFreq(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfCoil(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
GzAmp(GzTime<0 | GzTime>SEt(2)-SEt(1)) = [];
GyAmp(GyTime<0 | GyTime>SEt(2)-SEt(1)) = [];
GxAmp(GxTime<0 | GxTime>SEt(2)-SEt(1)) = [];
ADC(ADCTime<0 | ADCTime>SEt(2)-SEt(1)) = [];
Ext(ExtTime<0 | ExtTime>SEt(2)-SEt(1)) = [];
rfTime(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
GzTime(GzTime<0 | GzTime>SEt(2)-SEt(1)) = [];
GyTime(GyTime<0 | GyTime>SEt(2)-SEt(1)) = [];
GxTime(GxTime<0 | GxTime>SEt(2)-SEt(1)) = [];
ADCTime(ADCTime<0 | ADCTime>SEt(2)-SEt(1)) = [];
ExtTime(ExtTime<0 | ExtTime>SEt(2)-SEt(1)) = [];
rfAmp(abs(rfAmp)<eps) = 0;
rfTime = rfTime + SEt(1);
GzTime = GzTime + SEt(1);
GyTime = GyTime + SEt(1);
GxTime = GxTime + SEt(1);
ADCTime = ADCTime + SEt(1);
ExtTime = ExtTime + SEt(1);
rfAmpAll=[rfAmpAll rfAmp];
rfPhaseAll=[rfPhaseAll rfPhase];
rfFreqAll=[rfFreqAll rfFreq];
rfCoilAll=[rfCoilAll rfCoil];
GzAmpAll=[GzAmpAll GzAmp];
GyAmpAll=[GyAmpAll GyAmp];
GxAmpAll=[GxAmpAll GxAmp];
ADCAll=[ADCAll ADC];
ExtAll=[ExtAll Ext];
rfTimeAll=[rfTimeAll rfTime];
GzTimeAll=[GzTimeAll GzTime];
GyTimeAll=[GyTimeAll GyTime];
GxTimeAll=[GxTimeAll GxTime];
ADCTimeAll=[ADCTimeAll ADCTime];
ExtTimeAll=[ExtTimeAll ExtTime];
SEtAll=[SEtAll SEt];
end
%====================================
if isempty(rfTimeAll)
error('rf sequence line can not be empty! Master Tx coil element must be used.');
end
if isempty(GzTimeAll)
error('GzSS sequence line can not be empty!');
end
if isempty(GyTimeAll)
error('GyPE sequence line can not be empty!');
end
if isempty(GxTimeAll)
error('GxR sequence line can not be empty!');
end
if isempty(ADCTimeAll)
error('ADC sequence line can not be empty!');
end
if isempty(ExtTimeAll)
error('Ext sequence line can not be empty!');
end
SEflag=repmat([0 0 0 0 0 0]',[1 2]);
rfflag=repmat([1 0 0 0 0 0]',[1 max(size(rfTimeAll))]);
Gzflag=repmat([0 1 0 0 0 0]',[1 max(size(GzTimeAll))]);
Gyflag=repmat([0 0 1 0 0 0]',[1 max(size(GyTimeAll))]);
Gxflag=repmat([0 0 0 1 0 0]',[1 max(size(GxTimeAll))]);
ADCflag=repmat([0 0 0 0 1 0]',[1 max(size(ADCTimeAll))]);
Extflag=repmat([0 0 0 0 0 1]',[1 max(size(ExtTimeAll))]);
ts=[[min(SEtAll) max(SEtAll)] rfTimeAll GzTimeAll GyTimeAll GxTimeAll ADCTimeAll ExtTimeAll]-min(SEtAll);
flags=[SEflag rfflag Gzflag Gyflag Gxflag ADCflag Extflag];
[ts,ind]=sort(ts);
uts=unique(ts);
flags=flags(:,ind);
[rfTime,ind]=sort(rfTimeAll-min(SEtAll));
rfAmp=rfAmpAll(:,ind);
rfPhase=rfPhaseAll(:,ind);
rfFreq=rfFreqAll(:,ind);
rfCoil=rfCoilAll(:,ind);
[GzTime,ind]=sort(GzTimeAll-min(SEtAll));
GzAmp=GzAmpAll(:,ind);
[GyTime,ind]=sort(GyTimeAll-min(SEtAll));
GyAmp=GyAmpAll(:,ind);
[GxTime,ind]=sort(GxTimeAll-min(SEtAll));
GxAmp=GxAmpAll(:,ind);
[ADCTime,ind]=sort(ADCTimeAll-min(SEtAll));
ADC=ADCAll(:,ind);
[ExtTime,ind]=sort(ExtTimeAll-min(SEtAll));
Ext=ExtAll(:,ind);
rfAmp(1) = 0;
rfPhase(1) = 0;
rfFreq(1) = 0;
GzAmp(1) = 0;
GyAmp(1) = 0;
GxAmp(1) = 0;
ADC(1) = 0;
Ext(1) = 0;
rfAmp(end) = 0;
rfPhase(end) = 0;
rfFreq(end) = 0;
GzAmp(end) = 0;
GyAmp(end) = 0;
GxAmp(end) = 0;
ADC(end) = 0;
Ext(end) = 0;
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [Hkt, cx, cy, Wxk] = loglikGaPreadparamWfix(x,qqq)
% [Hkt, cx, cy, Wxk] = loglikGaPreadparam(varargin);

sigpsf = qqq{2};  %std deviation of the PSF gaussian approx
peval =  qqq{5}; %parameters
% Hkt_linear=x(1:end-peval.ncomp*2); % intensities
% cx=x(end-peval.ncomp*2+1:end-peval.ncomp); %x-coordinates of the centers
% cy=x(end-peval.ncomp+1:end); % y-coordinates of the centers
% Hkt_linear=qqq{6}; % intensities
Hkt_linear=x; % intensities
cx= qqq{6};
cy= qqq{7};

sigpsf_vec=repmat(sigpsf,peval.ncomp,1); %all psfs same sigma

% this is anoying as it depends when the simulation was performed/... grr
cxy_vec=[cx'+1, cy'+1]; %different notation of the dip_image/new version
% cxy_vec=[cx'+0.5, cy'+0.5]; %different notation of the dip_image/old version
a_vec=1./(sigpsf_vec.^2*2*pi); % all normalised to 1

Hkt=reshape(Hkt_linear, peval.ncomp, peval.nt);
% generate PSFs from given parameters:
Wxkpix=gauss2dmultislice([peval.nx, peval.ny, peval.ncomp], cxy_vec, sigpsf_vec, a_vec);
Wxkpix=normalizePSF(Wxkpix); %normalize PSFs to 1
Wxk=reshape(Wxkpix,peval.nx*peval.ny, peval.ncomp);
end

function Wnorm=normalizePSF(W)
sw=size(W);
if numel(sw)<3 %if ncomp=1...
    sw(3) = 1;
end
Wr=reshape(W, sw(1)*sw(2),sw(3));
q=squeeze(sum(Wr,1));
Wrnorm=Wr./repmat(q,sw(1)*sw(2),1);
Wnorm=reshape(Wrnorm,sw(1), sw(2), sw(3));
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function annotate(aline)
    %%
%     s = regexp(aline, '/', 'split');
%     if strcmp(s{4}, 'ShadowImages')
%         im = imread(aline);
%         file_name = [s{1},'/',s{2},'/',s{3},'/annotation/',s{5}];
%         file_name = [file_name(1:length(file_name)-4), '.mat'];
%         disp 'Segmenting'
%         [dummy seg] = edison_wrapper(im, @RGB2Luv, ...
%             'SpatialBandWidth', 9, 'RangeBandWidth', 15, ...
%             'MinimumRegionArea', 200);
%          seg = seg + 1;
%          s{4} = 'ShadowMasks'
%          mask_name = [s{1},'/',s{2},'/',s{3},'/',s{4},'/',s{5}];
%          mask_name = [mask_name(1:length(mask_name)-4), '.png'];
%          mask = imread(mask_name);
%          shadow = [];
%          lit = [];
%          segnum = max(max(seg));
%          temp = zeros(2, segnum);
%          for i = 1:size(mask, 1)
%              for j = 1:size(mask, 2)
%                  if mask(i,j) == 0
%                      temp(1, seg(i,j)) = temp(1, seg(i,j)) + 1;
%                  else
%                      temp(2, seg(i,j)) = temp(2, seg(i,j)) + 1;
%                  end
%              end
%          end
%          for i = 1:size(temp, 2)
%              if temp(1, i) <= temp(2, i)
%                  shadow = [shadow, i];
%              else
%                  lit = [lit, i];
%              end
%          end
%          save(file_name, 'seg', 'shadow', 'lit', 'segnum');
%     end

    %%
    s = regexp(aline, '/', 'split');
    if strcmp(s{3}, 'Our_test') && strcmp(s{4}, 'original')
        aline
        im = imread(aline);
        file_name = [s{1},'/',s{2},'/',s{3},'/annotation/',s{5}];
        file_name = [file_name(1:length(file_name)-4), '.mat'];
        disp 'Segmenting'
        [dummy seg] = edison_wrapper(im, @RGB2Luv, ...
            'SpatialBandWidth', 9, 'RangeBandWidth', 15, ...
            'MinimumRegionArea', 200);
         seg = seg + 1;
         s{4} = 'gt';
         mask_name = [s{1},'/',s{2},'/',s{3},'/',s{4},'/',s{5}];
         mask_name = [mask_name(1:length(mask_name)-4), '.png'];
         mask = imread(mask_name);
         mask = imresize(mask, [size(im,1), size(im,2)]);
         shadow = [];
         lit = [];
         segnum = max(max(seg));
         temp = zeros(2, segnum);
         for i = 1:size(mask, 1)
             for j = 1:size(mask, 2)
                 if mask(i,j) == 0
                     temp(1, seg(i,j)) = temp(1, seg(i,j)) + 1;
                 else
                     temp(2, seg(i,j)) = temp(2, seg(i,j)) + 1;
                 end
             end
         end
         for i = 1:size(temp, 2)
             if temp(1, i) <= temp(2, i)
                 shadow = [shadow, i];
             else
                 lit = [lit, i];
             end
         end
         save(file_name, 'seg', 'shadow', 'lit', 'segnum');
    end
end

author: kittenish

repo: Image-Shadow-Detection-and-Removal


text:
function [val,idx] = of_MARE(obs,sim,idx)
% of_MARE Calculates the Mean Absolute Relative Error of simulated
% streamflow. Ignores time steps with negative flow values.  Adds a constant
% e of 1/100 of mean(obs) to avoid issues with zero flows (Pushpalatha et
% al., 2012).

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% In:
% obs       - time series of observations       [nx1]
% sim       - time series of simulations        [nx1]
% idx       - optional vector of indices to use for calculation, can be
%               logical vector [nx1] or numeric vector [mx1], with m <= n
%
% Out:
% val       - objective function value          [1x1]
% idx       - vector of indices used for calculation

% Pushpalatha, R.; Perrin, C.; le Moine, N. and Andréassian V. (2012). "A
% review of efficiency criteria suitable for evaluating low-flow
% simulations". Journal of Hydrology. 420-421, 171-182.
% doi:10.1016/j.jhydrol.2011.11.055

%% Check inputs and select timesteps
if nargin < 2
    error('Not enugh input arguments')
end

if nargin < 3; idx = []; end
[sim, obs, idx] = check_and_select(sim, obs, idx);

%% Find the constant e
e = mean(obs)/100;

%% Apply constant and transform flows
obs = obs+e;
sim = sim+e;

%% Calculate metric
val = mean(abs((sim-obs)./obs));
end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function varargout = charGUI2(varargin)
% CHARGUI2 M-file for charGUI2.fig
%      CHARGUI2, by itself, creates a new CHARGUI2 or raises the existing
%      singleton*.
%
%      H = CHARGUI2 returns the handle to a new CHARGUI2 or the handle to
%      the existing singleton*.
%
%      CHARGUI2('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in CHARGUI2.M with the given input arguments.
%
%      CHARGUI2('Property','Value',...) creates a new CHARGUI2 or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before charGUI_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to charGUI2_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help charGUI2

% Last Modified by GUIDE v2.5 03-Jul-2008 17:55:56

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @charGUI2_OpeningFcn, ...
                   'gui_OutputFcn',  @charGUI2_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before charGUI2 is made visible.
function charGUI2_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to charGUI2 (see VARARGIN)
load data;
assignin('base','net',net);
% Choose default command line output for charGUI2
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes charGUI2 wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = charGUI2_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in pbLoad.
function pbLoad_Callback(hObject, eventdata, handles)
% hObject    handle to pbLoad (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[filename, pathname] = uigetfile({'*.bmp';'*.jpg';'*.gif';'*.*'}, 'Pick an Image File');
S = imread([pathname,filename]);
axes(handles.axes1);
imshow(S);

handles.S = S;
guidata(hObject, handles);

% --- Executes on button press in pbSelect.
function pbSelect_Callback(hObject, eventdata, handles)
% hObject    handle to pbSelect (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
S = handles.S;
axes(handles.axes1);

% Selection of location
if isfield(handles,'api')
    handles.api.delete();
    rmfield(handles,'api');
    rmfield(handles,'hRect');
    axes(handles.axes1);
    imshow(S);
end

axes(handles.axes1);
sz = size(S);
handles.hRect = imrect(gca,[round(sz(2)/2) round(sz(1)/2) 20 20]); % Select object
handles.api = iptgetapi(handles.hRect);
guidata(hObject, handles);

% img_crop = imcrop(S);
% axes(handles.axes2);
% imshow(img_crop);
%
% handles.img_crop = img_crop;

guidata(hObject, handles);

% --- Executes on button press in pbCrop.
function pbCrop_Callback(hObject, eventdata, handles)
% hObject    handle to pbCrop (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

handles.loc = handles.api.getPosition();
axes(handles.axes1);
S = handles.S;

handles.img_crop = imcrop(S,handles.loc);
axes(handles.axes2);
imshow(handles.img_crop);

% img_crop = imcrop(S);
% axes(handles.axes2);
% imshow(img_crop);
%
% handles.img_crop = img_crop;

guidata(hObject, handles);

% --- Executes on button press in pbPreprocess.
function pbPreprocess_Callback(hObject, eventdata, handles)
% hObject    handle to pbPreprocess (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

img_crop = handles.img_crop;
imgGray = rgb2gray(img_crop);
bw = im2bw(img_crop,graythresh(imgGray));
axes(handles.axes3);
imshow(bw);
bw2 = edu_imgcrop(bw);
axes(handles.axes4);
imshow(bw2);
handles.bw2 = bw2;
guidata(hObject, handles);

% --- Executes on button press in pbExtract.
function pbExtract_Callback(hObject, eventdata, handles)
% hObject    handle to pbExtract (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
bw2  = handles.bw2;
charvec = edu_imgresize(bw2);
axes(handles.axes5);
plotchar(charvec);
handles.charvec = charvec;
guidata(hObject, handles);

% --- Executes on button press in pbRecognize.
function pbRecognize_Callback(hObject, eventdata, handles)
% hObject    handle to pbRecognize (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
charvec = handles.charvec;
selected_net = get(handles.editNN,'string');

selected_net = evalin('base',selected_net);

result = sim(selected_net,charvec);
[val, num] = max(result);
set(handles.editResult, 'string',num);

% --- Executes on button press in pbNN.
function pbNN_Callback(hObject, eventdata, handles)
% hObject    handle to pbNN (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

function editNN_Callback(hObject, eventdata, handles)
% hObject    handle to editNN (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editNN as text
%        str2double(get(hObject,'String')) returns contents of editNN as a double

% --- Executes during object creation, after setting all properties.
function editNN_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editNN (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function editResult_Callback(hObject, eventdata, handles)
% hObject    handle to editResult (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of editResult as text
%        str2double(get(hObject,'String')) returns contents of editResult as a double

% --- Executes during object creation, after setting all properties.
function editResult_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editResult (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


text:
function mini_batch = extract_mini_batch(data_set, start_i, n_cases)
    mini_batch.inputs = data_set.inputs(:, start_i : start_i + n_cases - 1);
    mini_batch.targets = data_set.targets(:, start_i : start_i + n_cases - 1);
end

repo: neural-nets


text:
function outkey = filter(inkey,modlist,seglist,keep)
% Removes some of the information in a key.  Useful for creating a
% gender specific key from a pooled gender key.  Depending on the
% value of 'keep', the two input lists indicate the strings to
% retain or the strings to discard.
% Inputs:
%   inkey: A Key object.
%   modlist: A cell array of strings which will be compared with
%     the modelset of 'inkey'.
%   seglist: A cell array of strings which will be compared with
%     the segset of 'inkey'.
%   keep: A boolean indicating whether modlist and seglist are the
%     models to keep or discard.
% Outputs:
%   outkey: A filtered version of 'inkey'.

if nargin == 0
    test_this();
    return
else
    assert(nargin==4)
end

assert(isa(inkey,'Key'))
assert(iscell(modlist))
assert(iscell(seglist))
assert(inkey.validate())

if keep
    keepmods = modlist;
    keepsegs = seglist;
else
    keepmods = setdiff(inkey.modelset,modlist);
    keepsegs = setdiff(inkey.segset,seglist);
end

keepmodidx = ismember(inkey.modelset,keepmods);
keepsegidx = ismember(inkey.segset,keepsegs);

outkey = Key();
outkey.modelset = inkey.modelset(keepmodidx);
outkey.segset = inkey.segset(keepsegidx);
outkey.tar = inkey.tar(keepmodidx,keepsegidx);
outkey.non = inkey.non(keepmodidx,keepsegidx);

assert(outkey.validate())

if length(inkey.modelset) > length(outkey.modelset)
    log_info('Number of models reduced from %d to %d.\n',length(inkey.modelset),length(outkey.modelset));
end
if length(inkey.segset) > length(outkey.segset)
    log_info('Number of test segments reduced from %d to %d.\n',length(inkey.segset),length(outkey.segset));
end

end

function test_this()

key = Key();
key.modelset = {'aaa','bbb','ccc','ddd'};
key.segset = {'11','22','33','44','55'};
key.tar = logical([1,0,0,1,0;0,1,0,1,1;0,0,0,1,0;1,1,0,0,0]);
key.non = logical([0,1,0,0,0;1,0,0,0,0;1,1,1,0,0;0,0,1,1,1]);

fprintf('key.modelset\n');
disp(key.modelset)
fprintf('key.segset\n');
disp(key.segset)
fprintf('key.tar\n');
disp(key.tar)
fprintf('key.non\n');
disp(key.non)

modlist = {'bbb','ddd'}
seglist = {'11','55'}

keep = true

out = Key.filter(key,modlist,seglist,keep);

fprintf('out.modelset\n');
disp(out.modelset)
fprintf('out.segset\n');
disp(out.segset)
fprintf('out.tar\n');
disp(out.tar)
fprintf('out.non\n');
disp(out.non)

keep = false

out = Key.filter(key,modlist,seglist,keep);

fprintf('out.modelset\n');
disp(out.modelset)
fprintf('out.segset\n');
disp(out.segset)
fprintf('out.tar\n');
disp(out.tar)
fprintf('out.non\n');
disp(out.non)

modlist = {'bbb','ddd','eee'}
seglist = {'11','66','77','55'}

keep = true

out = Key.filter(key,modlist,seglist,keep);

fprintf('out.modelset\n');
disp(out.modelset)
fprintf('out.segset\n');
disp(out.segset)
fprintf('out.tar\n');
disp(out.tar)
fprintf('out.non\n');
disp(out.non)

keep = false

out = Key.filter(key,modlist,seglist,keep);

fprintf('out.modelset\n');
disp(out.modelset)
fprintf('out.segset\n');
disp(out.segset)
fprintf('out.tar\n');
disp(out.tar)
fprintf('out.non\n');
disp(out.non)

end

author: nesl


text:
% This file is part of the project NILM-Eval (https://github.com/beckel/nilm-eval).
% Licence: GPL 2.0 (http://www.gnu.org/licenses/gpl-2.0.html)
% Copyright: ETH Zurich, 2014
% Author: Romano Cicchetti

function HMMs = generateHMMsFromSnippets(snippets)

    % generate HMMs using snippets

    HMMs = struct;
    numOfSnippets = length(snippets.mean);
    HMMs.mean = cell(numOfSnippets,1);
    HMMs.std = cell(numOfSnippets,1);
    HMMs.transition = cell(numOfSnippets,1);
    snippetsNumOfStates = cellfun('length', snippets.mean) + 1;
    for i = 1:length(snippets.mean)
        numOfStates = snippetsNumOfStates(i);
        HMMs.mean{i} = [cell2mat(snippets.mean(i)); 0];
        HMMs.std{i} = [cell2mat(snippets.std(i))'; 1];
        transition = zeros(numOfStates);
        for state = 1:numOfStates
            if state == numOfStates
               transition(state, state) = 1;
            else
               duration = cell2mat(snippets.duration(i));
               transition(state, state) = 1 - 1/duration(state);
               transition(state, mod(state, numOfStates) + 1) = 1/duration(state);
            end
        end
        HMMs.transition{i} = transition;
    end

end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function out=subsasgn(x,varargin)

if isempty(x)
 x=mp(x);
end
out=builtin('subsasgn',x,varargin{:});

% This has to be fixed to allow the incoming doubles (or other mp's)
% to be promoted to the right precision
%'ssssssss',kb

author: opencobra

repo: cobratoolbox


text:
function [ res, jac ] = opt14_rj ( x, flag )

%*****************************************************************************80
%
%% OPT14_RJ evaluates RES and JAC for test case #14.
%
%  Modified:
%
%    07 January 2008
%
%  Author:
%
%    Jeff Borggaard,
%    Gene Cliff,
%    Virginia Tech.
%
%  Reference:
%
%    John Dennis, Robert Schnabel,
%    Numerical Methods for Unconstrained Optimization
%    and Nonlinear Equations,
%    SIAM, 1996,
%    ISBN13: 978-0-898713-64-0,
%    LC: QA402.5.D44.
%
%  Parameters:
%
%    Input, real X(3), the evaluation point.
%
%    Input, string FLAG, indicates what must be computed.
%    'f' means only the value of RES is needed,
%    'g' means only the value of JAC is needed,
%    'all' means RES and JAC are needed.
%    It is acceptable to behave as though FLAG was 'all'
%    on every call.
%
%    Output, real RES(3,1), the function column vector.
%
%    Output, real JAC(3,3), the Jacobian matrix.
%
  n = length ( x );

  if ( n ~= 3 )
    fprintf ( '\n' );
    fprintf ( 'OPT14_RJ - Fatal error!\n' );
    fprintf ( '  The input vector X should have length 3.\n'),
    fprintf ( '  Instead, it has length = %d.\n', n );
    keyboard
  end

  res = zeros(n,1);

  res(1,1) = x(1)^2*x(2) + x(1)*x(2)^2;
  res(2,1) = 3 * x(1) * x(2)^2 * x(3) - x(1) * x(3) - 1;
  res(3,1) = x(1)*x(3) - 2;

  jac = zeros(n,n);

  jac(1,1) = 2 * x(1) * x(2) + x(2)^2;
  jac(1,2) = x(1)^2 + 2 * x(1) * x(2);
  jac(1,3) = 0;

  jac(2,1) = 3 * x(2)^2 * x(3) - x(3);
  jac(2,2) = 3 * 2 * x(1) * x(2) * x(3);
  jac(2,3) = 3 * x(1) * x(2)^2 - x(1);

  jac(3,1) = x(3);
  jac(3,2) = 0;
  jac(3,3) = x(1);


text:
% -------------------------------------------------------------------------------------------------------------------------
function sz = size_ndims(x, n)
% -------------------------------------------------------------------------------------------------------------------------
sz = cell(1, n);
[sz{1:n}] = size(x);
sz = cell2mat(sz);
end

author: bertinetto

repo: siamese-fc


text:
function M = fixedrankfactory_2factors_subspace_projection(m, n, k)
% Manifold of m-by-n matrices of rank k with two factor quotient geometry.
%
% function M = fixedrankfactory_2factors_subspace_projection(m, n, k)
%
% A point X on the manifold is represented as a structure with two
% fields: L and R. The matrix L (mxk) is orthonormal,
% while the matrix R (nxk) is a full column-rank
% matrix such that X = L*R'.
%
% Tangent vectors are represented as a structure with two fields: L, R.
%
% Note: L is orthonormal, i.e., columns are orthogonal to each other.
% Such a geometry might be of interest where the left factor has a
% subspace interpretation. A motivation is in Sections 3.3 and 6.4 of the
% paper below.
%
% Please cite the Manopt paper as well as the research paper:
%     @Article{mishra2014fixedrank,
%       Title   = {Fixed-rank matrix factorizations and {Riemannian} low-rank optimization},
%       Author  = {Mishra, B. and Meyer, G. and Bonnabel, S. and Sepulchre, R.},
%       Journal = {Computational Statistics},
%       Year    = {2014},
%       Number  = {3-4},
%       Pages   = {591--621},
%       Volume  = {29},
%       Doi     = {10.1007/s00180-013-0464-z}
%     }
%
% See also: fixedrankfactory_2factors fixedrankembeddedfactory fixedrankfactory_2factors_preconditioned

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, Dec. 30, 2012.
% Contributors:
% Change log:

    M.name = @() sprintf('LR'' quotient manifold of %dx%d matrices of rank %d', m, n, k);

    M.dim = @() (m+n-k)*k;

    % Some precomputations at the point X to be used in the inner product (and
    % pretty much everywhere else).
    function X = prepare(X)
        if ~all(isfield(X,{'RtR'}) == 1)
            X.RtR = X.R'*X.R;
        end
    end

    % The choice of the metric is motivated by symmetry and scale
    % invariance in the total space.
    M.inner = @iproduct;
    function ip = iproduct(X, eta, zeta)
        X = prepare(X);

        ip = eta.L(:).'*zeta.L(:)  + trace(X.RtR\(eta.R'*zeta.R));
    end

    M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));

    M.dist = @(x, y) error('fixedrankfactory_2factors_subspace_projection.dist not implemented yet.');

    M.typicaldist = @() 10*k;

    skew = @(X) .5*(X-X');
    symm = @(X) .5*(X+X');
    stiefel_proj = @(L, H) H - L*symm(L'*H);

    M.egrad2rgrad = @egrad2rgrad;
    function rgrad = egrad2rgrad(X, egrad)
        X = prepare(X);

        rgrad.L = stiefel_proj(X.L, egrad.L);
        rgrad.R = egrad.R*X.RtR;
    end

    M.ehess2rhess = @ehess2rhess;
    function Hess = ehess2rhess(X, egrad, ehess, eta)
        X = prepare(X);

        % Riemannian gradient.
        rgrad = egrad2rgrad(X, egrad);

        % Directional derivative of the Riemannian gradient.
        Hess.L = ehess.L - eta.L*symm(X.L'*egrad.L);
        Hess.L = stiefel_proj(X.L, Hess.L);

        Hess.R = ehess.R*X.RtR + 2*egrad.R*symm(eta.R'*X.R);

        % Correction factor for the non-constant metric on the factor R.
        Hess.R = Hess.R - rgrad.R*(X.RtR\(symm(X.R'*eta.R))) - eta.R*(X.RtR\(symm(X.R'*rgrad.R))) + X.R*(X.RtR\(symm(eta.R'*rgrad.R)));

        % Projection onto the horizontal space.
        Hess = M.proj(X, Hess);
    end

    M.proj = @projection;
    function etaproj = projection(X, eta)
        X = prepare(X);

        eta.L = stiefel_proj(X.L, eta.L); % On the tangent space.
        SS = X.RtR;
        AS1 = 2*X.RtR*skew(X.L'*eta.L)*X.RtR;
        AS2 = 2*skew(X.RtR*(X.R'*eta.R));
        AS  = skew(AS1 + AS2);

        Omega = nested_sylvester(SS,AS);
        etaproj.L = eta.L - X.L*Omega;
        etaproj.R = eta.R - X.R*Omega;
    end

    M.tangent = M.proj;
    M.tangent2ambient = @(X, eta) eta;

    M.retr = @retraction;
    function Y = retraction(X, eta, t)
        if nargin < 3
            t = 1.0;
        end
        Y.L = uf(X.L + t*eta.L);
        Y.R = X.R + t*eta.R;

        % These are reused in the computation of the gradient and Hessian.
        Y = prepare(Y);
    end

    M.exp = @exponential;
    function R = exponential(X, eta, t)
        if nargin < 3
            t = 1.0;
        end

        R = retraction(X, eta, t);
        warning('manopt:fixedrankfactory_2factors_subspace_projection:exp', ...
            ['Exponential for fixed rank ' ...
            'manifold not implemented yet. Lsed retraction instead.']);
    end

    M.hash = @(X) ['z' hashmd5([X.L(:) ; X.R(:)])];

    M.rand = @random;
    % Factors L lives on Stiefel manifold, hence we will reuse
    % its random generator.
    stiefelm = stiefelfactory(m, k);
    function X = random()
        X.L = stiefelm.rand();
        X.R = randn(n, k);
    end

    M.randvec = @randomvec;
    function eta = randomvec(X)
        eta.L = randn(m, k);
        eta.R = randn(n, k);
        eta = projection(X, eta);
        nrm = M.norm(X, eta);
        eta.L = eta.L / nrm;
        eta.R = eta.R / nrm;
    end

    M.lincomb = @lincomb;

    M.zerovec = @(X) struct('L', zeros(m, k),...
        'R', zeros(n, k));

    M.transp = @(x1, x2, d) projection(x2, d);

    % vec and mat are not isometries, because of the scaled inner metric.
    M.vec = @(X, U) [U.L(:) ; U.R(:)];
    M.mat = @(X, u) struct('L', reshape(u(1:(m*k)), m, k), ...
        'R', reshape(u((m*k+1):end), n, k));
    M.vecmatareisometries = @() false;

end

% Linear combination of tangent vectors.
function d = lincomb(x, a1, d1, a2, d2) %#ok<INLSL>

    if nargin == 3
        d.L = a1*d1.L;
        d.R = a1*d1.R;
    elseif nargin == 5
        d.L = a1*d1.L + a2*d2.L;
        d.R = a1*d1.R + a2*d2.R;
    else
        error('Bad use of fixedrankfactory_2factors_subspace_projection.lincomb.');
    end

end

function A = uf(A)
    [L, unused, R] = svd(A, 0); %#ok
    A = L*R';
end

function omega = nested_sylvester(sym_mat, asym_mat)
    % omega=nested_sylvester(sym_mat,asym_mat)
    % This function solves the system of nested Sylvester equations:
    %
    %     X*sym_mat + sym_mat*X = asym_mat
    %     Omega*sym_mat+sym_mat*Omega = X
    % Mishra, Meyer, Bonnabel and Sepulchre, 'Fixed-rank matrix factorizations and Riemannian low-rank optimization'

    % Uses built-in lyap function, but does not exploit the fact that it's
    % twice the same sym_mat matrix that comes into play.

    X = lyap(sym_mat, -asym_mat);
    omega = lyap(sym_mat, -X);

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function [output, firstIndexPoint] = read_bti_hs( filename, outfile)

%read_hs_file Reads in BTI-Headshape files
%   filename: file with the headshape informations
%   outfile: if present, a ctf ".shape" file is written
%   output: if present, a 3xN matrix containing the headshape-points
%
%   (C) 2007 by Thomas Hartmann

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

if nargin == 1
    outfile = [];
end %if

fid       = fopen_or_error(filename, 'r', 'b');
version   = fread(fid, 1, '*uint32');
timestamp = fread(fid, 1, '*int32');
checksum  = fread(fid, 1, '*int32');
nPoints   = fread(fid, 1, '*int32');

firstIndexPoint = fread(fid, [3, 5], 'double')';

points = fread(fid, [3, double(nPoints)], 'double');

fclose(fid);

if(nargout > 0)
    output = points';
end %if

if(nargin == 2)
    fid = fopen_or_error(outfile, 'wt');
    fprintf(fid, '%d\n', nPoints);
    for i = 1:size(points, 2)
        fprintf(fid, '%.3f\t%.3f\t%.3f\n', points(1, i), points(2, i), points(3, i));
    end %for
    fclose(fid);

end %if

author: fieldtrip

repo: fieldtrip


text:
% VL_VLAD   VLAD feature encoding
%   ENC = VL_VLAD(X, MEANS, ASSIGNMENTS) computes the VLAD
%   encoding of the vectors X relative to cluster centers MEANS and
%   vector-to-cluster soft assignments ASSIGNMENTS.
%
%   X has one column per data vector (e.g. a SIFT descriptor), and
%   MEANS has one column per component. Usually one has one component
%   per KMeans cluster and MEANS are the KMeans centers. X and MEANS
%   have the same number of rows and the data class, which can be
%   either SINGLE or DOUBLE.
%
%   ASSIGNMENTS has as many rows as clusters and as many columns as
%   X. Its columns are non-negative and should sum to one,
%   representing the soft assignment of the corresponding vector in X
%   to each of the clusters. It is of the same class as X.
%
%   ENC is a vector of the same class of X of size equal to the
%   product of the data dimension and the number of clusters.
%
%   By default, ENC is L2 normalized. VL_VLAD() accepts the following
%   options:
%
%   Unnormalized::
%     If specified, no overall normalization is applied to ENC.
%
%   NormalizeComponents::
%     If specified, the part of the encoding corresponding to each
%     cluster is individually normalized.
%
%   NormalizeMass::
%     If specified, each component is re-normalized by the mass
%     of data vectors assigned to it. If NormalizedComponents is
%     also selected, this has no effect.
%
%   SquareRoot::
%     If specified, the signed square root function is applied to
%     ENC before normalization.
%
%   Verbose::
%     Increase the verbosity level (may be specified multiple times).
%
%   See: VLAD, VL_HELP().

% Authors: David Novotny and Andrea Vedaldi

% Copyright (C) 2013 David Novotny and Andrea Vedaldi
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function [feat,idxs_bbox_pair,cb1,cb2] = get_spatial_features_neighbour_locref(locations,idxs_bbox_pair,next_joints, locref)

cb1 = locations(idxs_bbox_pair(:,1), :);
cb2 = locations(idxs_bbox_pair(:,2), :);

delta = cb2 - cb1;

deltaReal_forward = delta + squeeze(locref(idxs_bbox_pair(:,2), 2, :));
deltaReal_backward = -delta + squeeze(locref(idxs_bbox_pair(:,1), 1, :));

deltaPred_forward = squeeze(next_joints(idxs_bbox_pair(:,1), 1, :));
deltaPred_backward = squeeze(next_joints(idxs_bbox_pair(:,2), 2, :));

feat = cat(2, deltaReal_forward, deltaReal_backward, ...
              deltaPred_forward, deltaPred_backward);

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function w = load_from_datasource(ds, chans, startt, endt, combine_waves, usewkaround)
   %load_from_datasource
   %
   % w = load_from_datasource(ds, chans, startt, endt, combine_waves)

   % determine the load function based upon the datasource's type

   if isVoidInterpreter(ds)
      ds = setinterpreter(ds, get_load_routine(ds, usewkaround));
   end
   request = packDataRequest(ds, chans, startt, endt, combine_waves);
   getter = get(ds,'interpreter');
   %w = getter(request, combine_waves);
   w = getter(request);

   % TODO: Trim waveform to request
   % maybe we're already done?
   return
   % -------------------------------------------------------------------
   % if the datasource is file based, or if it requires a user-defined
   % intepreter function, then do what follows.  Otherwise, we're done

   if ~isVoidInterpreter(ds)
      myLoadFunc = get(ds,'interpreter');

      %user_defined datasource
      ALLOW_MULTIPLE = true;
      ERROR_ON_EMPTY = true;
      for j = 1:numel(startt)
         myStartTime = startt(j);
         myEndTime = endt(j);

         % grab all files for date range, discarding duplicates
         fn = getfilename(ds, chans,subdivide_files_by_date(ds,myStartTime, myEndTime));
         fn = unique(fn);

         %load all waveforms for these files
         clear somew

         for i=1:numel(fn)
            w = load_from_file(myLoadFunc,fn{i},ALLOW_MULTIPLE, ERROR_ON_EMPTY);
            %w = w(ismember([w.cha_tag],[chans])); %keep only appropriate station/chan
            ct = w.cha_tag;
            w = w(ct.matching(chans)); %Not 100% sure about this one
            w = filter_by_time(w, myStartTime, myEndTime);
            if numel(w) > 0
               somew(i) = {w};
            end
         end %for fn
         if ~exist('somew','var')
            w = waveform; w = w([]);
         else
            if combine_waves,
               w = combine([somew{:}]);
            else
               w = [somew{:}];
            end
         end
         allw(j) = {w};
      end %each start time
      w = [allw{:}];
   end %~isVoidInterpreter
end

function tf = isVoidInterpreter(ds)
   tf = strcmpi(func2str(get(ds,'interpreter')),'void_interpreter');
end

author: geoscience-community-codes

repo: GISMO


text:
function R = rand_rotation(n)
  % Uniformly sample SO(n)
  %
  % R = rand_rotation(n)
  %
  % Inputs:
  %   n  dimension
  % Output:
  %   R  n by n rotation matrix
  %

  % "How to generate a random unitary matrix" [Maris Ozols 2006]
  % http://home.lu.lv/~sd20008/papers/essays/Random%20unitary%20[paper].pdf
  [Q,~] = qr(randn(n));
  s = (2*(rand>0.5)-1);
  Q(:,1) = Q(:,1)*s;
  %b = det(Q);
  b = s*(2*mod(n,2)-1);
  Q(:,2) = b*Q(:,2);
  assert(abs(det(Q)-1)<1e-10)
  % rename as R
  R = Q;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function data = readfile(path, varargin)
opts.asString = false;
opts = vl_argparse(opts, varargin);

assert(exist(path, 'file') == 2, 'File %s does not exist.', path);

fd = fopen(path, 'r');
data = textscan(fd, '%s', 'delimiter', '\n');
fclose(fd);

data = data{1};
data = cellfun(@strtrim, data, 'UniformOutput', false);
isNonAscii = cellfun(@(s) all(s < 128), data);
assert(all(isNonAscii), 'File %s contains non ascii characters.', path);
if opts.asString
  data = strjoin(data, '\n');
end

author: hpatches

repo: hpatches-benchmark


text:
function sizes = getVarSizes(obj, inputSizes)
%GETVARSIZES  Get the size of the variables
%   SIZES = GETVARSIZES(OBJ, INPUTSIZES) computes the SIZES of the
%   DagNN variables given the size of the inputs. `inputSizes` is
%   a cell array of the type `{'inputName', inputSize, ...}`
%   Returns a cell array with sizes of all network variables.
%
%   Example, compute the storage needed for a batch size of 256 for an
%   imagenet-like network:
%   ```
%   batch_size = 256; single_num_bytes = 4;
%   input_size = [net.meta.normalization.imageSize, batch_size];
%   var_sizes = net.getVarSizes({'data', input_size});
%   fprintf('Network activations will take %.2fMiB in single.\n', ...
%     sum(prod(cell2mat(var_sizes, 1))) * single_num_bytes ./ 1024^3);
%   ```

% Copyright (C) 2015 Andrea Vedaldi, Karel Lenc.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

nv = numel(obj.vars) ;
sizes = num2cell(NaN(nv, 4),2)' ;

for i = 1:2:numel(inputSizes)
  v = obj.getVarIndex(inputSizes{i}) ;
  sizes{v} = inputSizes{i+1}(:)' ;
end

for layer = obj.layers
  in = layer.inputIndexes ;
  out = layer.outputIndexes ;
  sizes(out) = layer.block.getOutputSizes(sizes(in)) ;
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function tests = test_ft_preproc_online_filter

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_preproc_online_filter_init ft_preproc_online_filter_apply

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

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

nchan   = 2;
nsample = 100;

result = {};

[B,A] = butter(6, 0.05, 'high');
state = ft_preproc_online_filter_init(B, A, zeros(nchan,1));

[state, result{end+1}] = ft_preproc_online_filter_apply(state, randn(nchan, nsample) + 10);
[state, result{end+1}] = ft_preproc_online_filter_apply(state, randn(nchan, nsample) + 10);
[state, result{end+1}] = ft_preproc_online_filter_apply(state, randn(nchan, nsample) + 10);
[state, result{end+1}] = ft_preproc_online_filter_apply(state, randn(nchan, nsample) + 10);
[state, result{end+1}] = ft_preproc_online_filter_apply(state, randn(nchan, nsample) + 10);

[B,A] = butter(6, 0.20, 'low');
state = ft_preproc_online_filter_init(B, A, zeros(nchan,1));

[state, result{end+1}] = ft_preproc_online_filter_apply(state, randn(nchan, nsample) + 10);
[state, result{end+1}] = ft_preproc_online_filter_apply(state, randn(nchan, nsample) + 10);
[state, result{end+1}] = ft_preproc_online_filter_apply(state, randn(nchan, nsample) + 10);
[state, result{end+1}] = ft_preproc_online_filter_apply(state, randn(nchan, nsample) + 10);
[state, result{end+1}] = ft_preproc_online_filter_apply(state, randn(nchan, nsample) + 10);

[B,A] = butter(6, [0.05 0.20]);
state = ft_preproc_online_filter_init(B, A, zeros(nchan,1));

[state, result{end+1}] = ft_preproc_online_filter_apply(state, randn(nchan, nsample) + 10);
[state, result{end+1}] = ft_preproc_online_filter_apply(state, randn(nchan, nsample) + 10);
[state, result{end+1}] = ft_preproc_online_filter_apply(state, randn(nchan, nsample) + 10);
[state, result{end+1}] = ft_preproc_online_filter_apply(state, randn(nchan, nsample) + 10);
[state, result{end+1}] = ft_preproc_online_filter_apply(state, randn(nchan, nsample) + 10);

% the first  part is high-pass filtered
% the second part is low-pass  filtered
% the third  part is band-pass filtered
plot(cat(2, result{:})');

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