Octave

Octave er et programmeringsprog til numeriske beregninger. Det minder meget om MATLAB, der benyttes meget i akademia til computerdrevne beregniner, men hvor MATLAB er kommercielt, er Octave opensource og gratis. På denne side vil vi kigge lidt nærmere på Octave.

Har du programmeret før, er det formentlig piece of cake for dig at sætte dig ned og bruge Octave. Har du ikke programmeret før, kan det ikke skade at lære lidt om det. At kunne programmere er en af de mest anvendelige færdigheder nu om stunder. Det er næsten ligegyldigt hvilket programmeringsprog man lærer i første omgang; for når man kan det grundlæggende i et kan man læse og forstå det grundlæggende i de fleste andre.

Download

Octave kan downloades her:

https://www.gnu.org/software/octave/

Tutorials

Der eksisterer allerede en del guides og tutorials for Octave rundt omkring på nettet; og endnu mere for MATLAB, der ligner så meget, at det meste virker i begge sprog. En guide, jeg umiddelbart mener, virker fornuftig er denne:

http://www-mdp.eng.cam.ac.uk/web/CD/engapps/octave/octavetut.pdf

Tag dig tid og prøv kommandoerne, han beskriver, af. Den bedste måde at lære et programmeringssprog at kende på er 1) at lege med det, 2) at have et problem, man forsøger at løse med det. Søg på Google, hvis der er noget, du er i tvivl om, hvordan du skal gøre. Alle programmører søger hjælp online næsten hele tiden.

Jeg lavede et nogle scripts til Octave i 2016 samt optog et par videoer. Disse kan findes på det gamle website her.

AIR Tools

Når du er føler dig klar på at eksperimentere mere har Professor Per Christian Hansen fra DTU Compute en software pakke AIR Tools, du kan hente. AIR Tools, hvilket står for Algebraic Image Reconstruction Tools, er godt nok skrevet til MATLAB, men virker umiddelbart også i Octave. Oplever du problemer, som ikke er adresseret på siden her, så send mig en mail med problemet, og jeg skal forsøge at hjælpe med at løse det.

http://www2.compute.dtu.dk/~pcha/AIRtools/

Prøv først at køre fx ARTdemo og kig i koden for denne.

En anden brugbar funktion er paralleltomo. Dette funktion laver matricen $ A $ som blev omtalt i Lineær Algebra afsnittet, hvis du giver den grid-størrelsen $ N $, vinklerne $ \theta $ og hvor mange parallelle stråle du vil have $ p $. Pas på med at vælge $ N $ for stor, problemet vokser meget hurtigt i størelse. Når du har lavet matricen $ A $, så prøv eventuelt at køre spy(A) for at se strukturen af $ A $. Spy viser dig hvor matricen er forskellig fra 0. Her er et stykke eksempelkode:

N = 50;           % grid-størrelse
theta = 0:5:179;  % vinkler
p = 75;           % antal beams
A = paralleltomo(N, theta, p);
spy(A);

FAQ

Warning, division by zero
Octave giver en advarsel når man dividerer med 0, hvilket producerer resultatet inf (eller evt. -inf). Per Christian benytter et trick i sin kode, hvor der nogle gange divideres med 0. For at få Octave til at stoppe med at give advarslen kan: warning ("off", "Octave:divide-by-zero"); kommandoen benyttes.
error: '<...>' undefined near ...
Denne fejl skyldes formentlig, at du forsøger at køre et script (eller en funktion) uden 1) at være i folderen, hvor scriptet ligger, eller 2) at have tilføjes folderen, hvor scriptet ligger, til din path. Octave's GUI har en "Current Directory" bar øverst, hvor folderen du pt. opholder dig i er angivet. Ellers kan kommandoen pwd også fortælle dig det. Brug GUI'en eller cd kommandoen til at komme ind i den rette folder. Alternativt kan du tilføje stien til folderen, hvor scriptet du vil køre er i, til din path. Det kan gøres med addpath kommandoen. Se alle dine nuværende paths med path.

Scripts

Følgende script tegner et grid og en stråle, samt udskriver distancerne strålen bevæger sig gennem de forskellige kasser og hvilke kasser det er. Du kan lege med initialiseringen for forskellige grid og stråler.

clc        % Clear CLI
close all  % Close all windows
clear all  % Clear all variables

% Initialize :: Grid
N = 10;  x = linspace(-5,5,N+1)';
M = 10;  y = linspace(-5,5,M+1)';
% Initialize :: Beam
% Theta should be between 0 and pi.
theta = pi/4; rho = 0;

% Draw the grid and beam
figure; hold on;
for i = 1:length(x)
    plot(x(i)*[1,1],y([1,end]),'b-')
end
for i = 1:length(y)
    plot(x([1,end]),y(i)*[1,1],'b-')
end
t = linspace(-sqrt(y(1)^2 + x(1)^2), sqrt(y(end)^2 + x(end)^2), 2);
plot(cos(theta)*t - rho*sin(theta), sin(theta)*t + rho*cos(theta), 'r--')
axis equal

% Compute grid intersection times
tx = (x + rho*sin(theta))/cos(theta);
ty = (y - rho*cos(theta))/sin(theta);

% Gather times
t = [tx; ty];
% Gather coordinates
xxy = [x; cos(theta)*ty-rho*sin(theta)];
yyx = [sin(theta)*tx+rho*cos(theta); y];

% Filter out points outside grid
idx = x(1) <= xxy & xxy <= x(end) & y(1) <= yyx & yyx <= y(end);
t_ = t(idx);
xxy_ = xxy(idx);
yyx_ = yyx(idx);

% Sort points by intersection time
[t__,idx] = sort(t_);
xxy__ = xxy_(idx);
yyx__ = yyx_(idx);

% Find distances in each grid box
d = diff(t__);

% Find repeated points and remove them
idx = d < 1e-10;
d(idx) = [];
xxy___ = xxy__([~idx;true]);
yyx___ = yyx__([~idx;true]);

% List distances in boxes
disp('Dinstances: d = ')
disp(d') % <-- displays distances
%disp([ xxy___-x(1), yyx___-y(1) ]') <-- displays intersection points

% Compute grid index for boxes intersected by the beam
Lx = x(end)-x(1);
Ly = y(end)-y(1);
j = floor( (xxy___-x(1))*N/Lx )-1*(theta*2>pi);
k = floor( (yyx___-y(1))*M/Ly );
%disp( [j,k]' ) % <-- displays coordinates

% Makes grid-matrix and fills in distances
G = zeros(N,M);
for i = 1:length(d)
    G(1+j(i)+N*k(i)) = d(i);
end
G = rot90(G);

% Show G in CLI
disp('Grid matrix: G = ')
disp(G) <-- displays grid matrix