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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
This comment has been minimized.
E3V3A commentedMar 25, 2019
Very Nice! Thanks.