Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Octave and Matlab Snippets
%% 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,d) 0.5 * (max(x,[],d) + min(x,[],d)); % midrange along dimension d
arcl = @(x,y) sum(sqrt(diff(x).^2 + diff(y).^2)); % simple arc length approximation from coordinate vectors
eoc = @(y,y1,y2,h1,h2) log(norm(y1 - y,'fro')/norm(y2 - y,'fro')) / log(h1 / h2); % experimental order of convergence
coh = @(x,y) sum(x(:))*sum(y(:)) / sum(x(:).*y(:)); % coherence
tau = 2.0 * pi; % define tau constant
%% 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
ainv = @(m) ((m - diag(diag(m))) .* (-1.0./diag(m))) .* (1.0./diag(m)') + diag(1.0./diag(m)); % rough approximate inverse
[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* ( (b/a).^rand(k,l) ); % distribution between [a,b] uniformly over orders of magnitude
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
saxpy = @(a,x,y) a * x + y; % scalar a times x plus y
symtridiag = @(n,a,b) toeplitz([a,b,zeros(1,n-2)]); % symmetric tri-diagonal matrix
valmat = @(r,c,v) repmat(v,r,c); % create matrix with r rows and c columns filled with value v (also NaN)
%% 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
nornorm = @(m) norm(m*m' - m'*m,'fro'); % check how non-normal 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
gmean = @(A,d) prod(A,d).^(1/size(A,d)); % geometric mean
%% 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
ylim([10^floor(log10(min([y(:)]))-1),1]); % set lower log y limit based on data y
ytickformat = @(f) set(gca,'yticklabel',arrayfun(@(x) sprintf(f,x),get(gca,'ytick'),'UniformOutput',false)); % tick formatting
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
clim = @(a,b) set(gca,'CLim',[a,b]); % set colorbar limits
%% System %%
exist('OCTAVE_VERSION','builtin') % test if gnu octave is used
verLessThan('matlab','9.1') % test if mathworks matlab version is less than 2016b
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
fuse = onCleanup(@() clear myvar); % global destructor triggered by destruction of fuse
chars = @(n,c) arrayfun(@() c,blanks(n)); % create n-dim vector of characters c
%% Functional %%
vec = @(m) m(:); % useful function octave has but matlab has not
getelem = @(m,r,c) m(r,c); % get matrix element
setelem = @(m,r,c,v) m + sparse(r,c,-m(r,c)+v,size(m,1),size(m,2)); % set matrix element
cmov = @(c,varargin) varargin{2 - c}(); % conditional set: if c then varargin{1} else varargin{2}
foreach = @arrayfun; % note that Matlab's version of for-each is arrayfun (or cellfun)
%% 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 argument 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 GB
caller = @() (dbstack).name; % get current running function name
addpath(genpath('mypath')); % Add search path recursively for subfolders
picked = @(list,name) any(strcmp(list,name)); % Check if string exists inside cell array
%% Misc %%
[Amin,Amax] = bounds(A,d); % minimum and maximum of matrix A along dimension d
maxabs = @(m) max(abs(m(:))); % Maximum absolute value in matrix
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
expand = @(a) [zeros(isempty(a)),a]; % expand argument to zero if it empty else return argument
L = idiv(a,b); % Compute floor(a/b), useful ie for number of time steps from time horizon and step-width
[r,n] = max(abs(diff(log(V)))); % find the largest change in consecutive elements of a vector
n = 10.0.^round(log10(a)); % round to next order of magnitude
A(A == 0) = NaN; % change zero entries to NaN; useful for log plots
ascii = @() [num2str((32:127)'),repmat(' ',[96,1]),char(32:127)']; % ASCII Table
%% System Theory %%
l0norm = @(y) sum(prod(abs(y),1).^(1/size(y,1)),2); % L0-norm of time series y (states x time-steps)
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
%% Notes %%
% use not() instead of ~() inside if conditions
% use end%if, end%for, etc for Matlab compatible specific ends in Octave
% use "if" and "for" without parenthesis for better readability
% Numerical arrays are stored in continous memory, cell arrays are not.
% Every .m file should have a header with: project, version, authors, licence, summary
%% 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
svd_driver("gesvd"); # set safe LAPACK SVD 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(args)" > my.log & # 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
mlint('mycode.m','-cyc') % static code analysis and cyclomatic complexity in MATLAB
flamegraph(profdata) # use flame graph visualization for Octave profiler data ( https://git.io/flamegraph )
## 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 MATLAB BLAS backend
LAPACK_VERSION ='/path/to/my/lapack.so' # change MATLAB LAPACK backend
# use FlexiBLAS ( https://www.mpi-magdeburg.mpg.de/projects/flexiblas )
## Postprocessing ##
pdfcrop figure.pdf cropped.pdf # remove white-space framing plots
## Develop ##
http://wiki.octave.org/Nano # syntax highlighting for the nano editor
% by: https://gramian.de
@E3V3A

This comment has been minimized.

Copy link

E3V3A commented Mar 25, 2019

Very Nice! Thanks.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.