Octave and Matlab Snippets

matlab-octave.m

%% Math %%

si = @(x) sin(x)./x; % cardinal sine without pi multiplied argument

hsin = @(x) 0.5*(1.0 - cos(x)); % haversed sine

hcos = @(x) 0.5*(1.0 + cos(x)); % haversed cosine

sigm = @(x,k) 0.5*tanh(0.5*k*x) + 0.5; % sigmoid function to (exp(-kx)+1)^-1

midr = @(x) 0.5*(max(x,[],2)+min(x,[],2)); % midrange

arcl = @(x,y) sum(sqrt(diff(x).^2+diff(y).^2)); % simple arc length approximation from coordinate vectors


%% Matrix %%

A(1:N+1:end) = 1.0; % set square N-dim matrix diagonal to one

a = A(A>0); % get all elements greater than zero

dinv = @(D) diag(1.0./(diag(D))); % invert diagonal matrix

[r,c] = find(A==max(A(:))); % get row and column index of maximum element in matrix

erand = @(k,l,m) -log(rand(k,l))./m; % exponential distribution

sperand = @(k,l,m) (-log(rand(k,l))./m).*(sprand(k,l,d)>0); % sparse exponential distribution

lrandn = @(k,l,m,v) exp(m + v*randn(k,l)); % log-normal distribution

splrandn = @(k,l,m,v,d) exp(m + v*randn(k,l)).*(sprand(k,l,d)>0); % sparse log-normal distribution

randab = @(k,l,a,b) a*( (a/b).^rand(k,l) ); % distribution between [a,b] uniformly over orders of magnitude

nans = @(r,c) NaN*zeros(r,c); % create matrix of NaNs

A(1000,2000) = 0; % pre-allocate matrix faster than with zeros

yey = @(n) fliplr(eye(n)); % anti-diagonal "identity"-matrix

adiag = @(A,k) diag(flipud(A),k); % get anti-diagonal entries

LINSPACE = @(a,b,n) a + (b-a).*linspace(0,1.0,n); % linspace variant for vector-valued arguments a,b

gramian = @(m) m*m'; % compute gramian matrix

isnormal = @(m) all(all(abs((m*m')-(m'*m))<1e-15)); % test if (small) matrix is normal

pad = @(m,r,c) [m,zeros(size(m,1),max(c-size(m,2),0));zeros(max(r-size(m,1),0),max(c,size(m,2)))]; % pad matrix with zeros

nanmax = @(a,b) max(a,b) + (a-a) + (b-b); % NaN propagating maximum

nanmin = @(a,b) min(a,b) + (a-a) + (b-b); % NaN propagating minimum


%% Linear Algebra %%

[U,D,V] = svd(A,'econ'); % faster singular value decomposition for non square matrices

[U,d,v] = svds(X,1); % principal POD mode (in U)

msvds = @(A,r) sqrt(eigs((A'*A),r)); % memory-economical sparse svd

trnorm = @(m) sum(svd(m)); % Trace norm (aka nuclear norm) of a matrix

L = chol(A+norm(A,1))-norm(A,1); % approximate cholesky decomposition for non positive definite matrices

logdet = @(m) 2.0*sum(log(diag(chol(m)))); % faster log(det(A))

prodtrace = @(a,b) sum(sum(a.*b')); % Faster than trace(a*b)

sympart = @(m) 0.5*(m+m'); % symmetric part of matrix

symnorm = @(m) norm(m-m',1); % check how un-symmetric a matrix is

unit = @(n,N) (1:N==n)'; % generate n-th N-dim unit vector

units = @(n,N) sparse(n,1,1,N,1); % generate sparse n-th N-dim unit vector 

specrad = @(A) abs(eigs(A,1,'lm')); % spectral radius

%% ODE Solver %%

deeval = @(t,x,k) interp1(t,x,k(:))'; % Interpolate ODE solution [t,x] at vector k time points 


%% Plotting %%

imshow(full(A)); % display dense or sparse matrix

imwrite(full(A),[filename,'.png']); % save dense or sparse matrix to a png

imshow(dicomread('filename'),[]); % display DICOM image with automatic range

gplot(sprand(N,N,d),[cos(2*pi/N:2*pi/N:2*pi)',sin(2*pi/N:2*pi/N:2*pi)']); % visualize networks adjacency matrix

n = 2.0^nextpow2(N); Z = fft(Y',n)/N; loglog(2*abs(Z(1:n/2+1))); % Bode-plot for a (outputs x steps) matrix timeseries

c = hsv(N); % get N colors from the hsv color map

set([gca; findall(gca,'Type','text')],'FontSize',16); % scale all fonts in a plot

figure('Name',mfilename,'NumberTitle','off'); % set figure title bar content

print('-dsvg',[mfilename,'.svg']); % save current figure as svg with running m-file as base name


%% System %%

exist('OCTAVE_VERSION','builtin') % test if gnu octave is used

warning off all; % supress all warnings

rand('seed',s); % set seed s of uniform random number generator

randn('seed',s); % set seed s of normal random number generator

wfk = @() eval('fprintf(''Any Key!\n''); pause;'); % wait for key


%% Functional %%

in = @(v,k) v(k); % functional vector component extraction 

vec = @(m) m(:); % useful function octave has but matlab has not


%% Utils %%

fprint('Hello World'); % disp without line-end

say = @(m) fprintf([m,'\n']); % emit text line

if(nargin<1 || isempty(a)), x = 0; end % default argument value for parameter a

system(['notify-send "',mfilename,':\ I am done!"']); % local notification

system('mutt -s "I am done!" me@host.tld -a result.svg < nohup.out'); % remote notification

varsize = @(m) whos(m).bytes/(1024*1024*1024); % memory footprint of variable (name as string) in BG

caller = @() (dbstack).name; % get current running function name

addpath(genpath('mypath')); % Add search path recursively for subfolders


%% Misc %%

minmax = @(A) [min(A(:)),max(A(:))]; % minimum and maximum of matrix

[r,n] = max(abs(diff(log(V)))); % find the largest change in consecutive elements of a vector

n = 10^round(log(a)); % round to next order of magnitude

ascii = @() [num2str((32:127)'),repmat(' ',[96,1]),char(32:127)']; % ASCII Table


%% System Theory %%

l1norm = @(y,h) h*norm(y(:),1); % L1-norm of time series y (states x time-steps (h))

l2norm = @(y,h) sqrt(h)*norm(y(:),2); % L2-norm of time series y (states x time-steps (h))

l8norm = @(y) norm(y(:),Inf); % Linfinity-norm of time series y (states x time-steps)

wc = @(A,B) sylvester(A,A',-B*B'); % algebraic controllability gramian 

wo = @(A,C) sylvester(A',A,-C'*C); % algebraic observability gramian

wx = @(A,B,C) sylvester(A,A,-B*C); % algebraic cross gramian


%% MATLAB %%

version('-blas') % display BLAS identifier and version

version('-lapack') % display LAPACK identifier and version

profile -memory on; % activate memory usage profiling

matlabpool(feature('numCores')); % allocate maximum local workers


## Octave ##

page_screen_output(0); # octave command to disable paged output

page_output_immediately(1); # octave command to enable immediate output

history_control("ignoredups"); # do not record duplicate entries in command history

graphics_toolkit("dumb"); # set render backend for in-terminal ascii plots

graphics_toolkit("gnuplot"); # set render backend to gnuplot

CURRENT_BLAS = flexiblas_current_backend # display current BLAS backend

#!/usr/bin/octave-cli # command-line octave shebang


## Remote ##

nohup matlab -nodisplay -r progname < /dev/null & # remote matlab execution

nohup octave-cli --eval progname & # remote octave execution


## Parallel ##

taskset -c 0,2,4,6 octave-cli # set CPU affinity to every second core (for SMT machines)

numactl -N 0 -m 0 octave-cli # pin to first node and its memory (for NUMA machines)

lstopo # print CPU topology (part of hwloc)

matlab2016b -nodisplay -singleCompThread # (use only single thread in terminal)


## Profiling ##

/usr/bin/time -f "%M KB" octave # record peak memory usage


## Custom Backends ##

LD_PRELOAD='/path/to/my/blas.so /path/to/lapack.so' # change BLAS and LAPACK backends for Octave

BLAS_VERSION='/path/to/my/blas.so' # change BLAS backend for Matlab

LAPACK_VERSION ='/path/to/my/lapack.so' # change LAPACK backend for Matlab


## Develop ##

http://wiki.octave.org/Nano ? syntax highlighting for the nano editor


% by: http://gramian.de