//*************************************************//
//                SciLab-Tutorial                  //
//   Numerical Algorithms for Visual Computing III //
//                 Kai Hagenburg                   //
//       based upon treatment by Peter Urban       //
//                copyright 2009                   //
//*************************************************//

// Some matrix and vector operations

// declaration of matrices and vectors
// matrix A (size 3x3)
A=[1 2 3; 2 3 4; 5 6 7];

// vector v1 of size (1x3)

v1=[1 2 3];

// vector v2 of size (3x1)

v2=[1;2;3];

// or alternatively

v2=v1';

// the '-operator gives the transpose of a matrix or a vector
// use semicolon if you don't want to output the 
// matrix again



// special functions:
// min(A)  -> minimum of a matrix
// max(A)  -> maximum of a matrix
// real(A) -> real part of a complex matrix
// imag(A) -> imaginary part of a complex matrix
// spec(A) -> gives eigenvalues of a matrix
// sum (A) -> gives sum of matrix
// max (A) -> gives maximum of matrix
// min (A) -> gives minimum of matrix

// sum(A,'c') -> computation of columnwise sum
// sum(A,'r') -> computation of rowwise sum

// special numbers
// %e   -> Euler's number
// %pi  -> Pi
// %i   -> something imaginary


// matrix-access:

//gives value a_{1,1} of the matrix
A(1,1);

//gives the entire first row of the matrix
A(1,1:3);

//gives the first column
A(1:3,1);

//give out the size of the dimensions of a matrix
//and store the results in values m,n
[m,n]=size(A);

//for vectors
length(v1);

//create linear entries in a matrix
B=[1:10];

//with some values in between
B=[1:0.5:10];

//with some values missing in between
B=[1:2:10];

//can also be done with the command linspace(a,b,n)
//a=left border
//b=right border
//n=sample points between borders
B=linspace(1,10,10);
B=linspace(1,10,19);

//Matrix Algebra Operations

//Addition:
A+A;

//Multiplication:
A*A;

//Powers
A^5;

//Transpose
A';

//Transpose for complex-valued matrices gives conjugate transpose
[1 2; 3+%i 4+%i]';

//More interesting matrix operations
//Solving linear systems of equations:
//A*x=b (if A is non-singular)
A=[1 2; 3 1];
b=[5 5]';

x=A\b;

//inverse of matrix
inv(A);

//diagonal entries of a matrix, stored in a vector
diag(A);


//Programming:

// functions:
// x is the output value
// if more than one value should be output: 
// function [x1,x2,x3,...] = yourfunction(...)

// gives out faculty (recursively)

function x = fac(n)
  if (n<>1)
    x=n*fac(n-1)
  else 
    x = 1
  end
endfunction


// if-Clause:
// if (argument)
//   statement
// else / elseif
//   statement
// else / elseif
//    ...
// end


// for argument:
// &  : and
// |  : or
// ~  : neg
// == : equal
// <  : less than
// >  : bigger than
// <= : less or equal than
// >= : bigger or equal than
// <> : inequality


// gives out faculty (iteratively)
// with for-loop
function x = faciter(n)
  x=1
  for i=2:n
    x=i*x
  end
endfunction

// gives out faculty (the fast Scilab way)
function x=facscilab(n)
  x=prod(1:n)
endfunction

//take function values of 1/x in the interval [1..10] and plot them
for i=1:10 vec1(i)=1/i; end
//plot(vec1);


//better refinement
//build more discretisation points
//blow up vector and rescale function

for i=1:100 vec2(i)=10/i; end
//plot(vec2);

//then build proper x_axis
//plot with x_axis
x_axis=[0.1:0.1:10];
//plot(vec2,x_axis);


//Now two different versions of SOR
//for more details on the numerical side,
//please refer to the NAVC2 script from
//Winter Semester 2008/2009, chapter 9



A=[2 -1; -1 2];
b=[4 1]'
x=[0 0]'
// solution=[3 2]'

// helper functions
function D = diago(A)
//first application of diag on matrix:
// returns a vector of diagonal entries
//second application of diag:
// returns diagonal matrix with entries given by the vector
D=diag(diag(A))
endfunction

function U=rightify(A)
  //triu gives upper triangle matrix with the diagonal entries
  U=triu(A)-diag(diag(A))
endfunction

function L=leftify(A)
  //triu gives lower triangle matrix with the diagonal entries
L = tril(A)-diag(diag(A))
endfunction


function x = SOR(A,b,iter,x)
  //linear system of equation given by Ax=b
  //x is here initialization, e.g. x=[0...0]'
  
  //starts a timer function
  timer()
  n = length(b)

  // compute triangular and diagonal matrices
  D = diago(A)
  R = rightify(A)
  L = leftify(A)

  // compute best value for omega parameter
  M1 = inv(D)*(D-A)
  rho=max(real(spec(M1)))  
  omega = 2/(1+sqrt(1-rho^2))
  
  
  printf("Best omega:%f\n",omega)
  M = inv(D+omega*L)*((1-omega)*D-omega*R)
  N = omega*inv(D+omega*L)
  printf("Maximal eigenvalue of the system matrix is %f\n",max(real(abs(spec(M)))))
  for i=1:iter
    x=M*x+N*b
end

//compute absolute error
sepsilon = abs(x-inv(A)*b)

//norm gives the L2-error
printf("L2-Error: %f\n",norm(x-inv(A)*b))
printf("time: %f\n",timer())
endfunction

function x = SOR_1(A,b,iter,x,omega)
n = length(b)
for l=1:iter
  tmp=x
  for i=1:n
    helper_1 = 0
    helper_2 = 0
   
    for j=1:i-1
      helper_1 = helper_1 + A(i,j)*tmp(j)
    end
    for j=i+1:n
      helper_2 = helper_2 +A(i,j)*x(j)
    end
    tmp(i)=(1-omega)*x(i)+omega/A(i,i)*(b(i)-helper_1 - helper_2)
end
x = tmp
end
endfunction



function x = SOR_2(A,b,x)
  //linear system of equation given by Ax=b
  //x is here initialization, e.g. x=[0...0]'
  
  //starts a timer function
  timer()
  n = length(b)

  // compute triangular and diagonal matrices
  D = diago(A)
  R = rightify(A)
  L = leftify(A)

  // compute best value for omega parameter
  M1 = inv(D)*(D-A)
  rho=max(real(spec(M1)))  
  omega = 2/(1+sqrt(1-rho^2))
  
  
  printf("Best omega:%f\n",omega)
  M = inv(D+omega*L)*((1-omega)*D-omega*R)
  N = omega*inv(D+omega*L)
  bla = A*b
 sol = inv(A)*b
  while (norm(x-sol) > 0.00000001)
    x=M*x+N*b
end

//compute absolute error
//sepsilon = abs(x-inv(A)*b)

//norm gives the L2-error
printf("L2-Error: %f\n",norm(x-sol))
printf("time: %f\n",timer())
endfunction




// create random vector : 
// uniform distribution of integer values on the interval [0,25]:
// vector of size (1,100)
// build a step-like signal with some random noise

a1 = 0
b1 = 25
noise = floor(a1 + (b1+1-a1)*rand(100,1))

signal(1:50)=100
signal(51:100)=200
signal = noise+signal


// Linear Diffusion
// Explicit scheme

function x = lindiff(x,iter,tau)
  n  = length(x)
  for t=1:iter
    temp(1)=(1-tau)*x(1)+tau*x(2)
    for i=2:n-1
      temp(i)= (1-2*tau)*x(i)+tau*x(i-1)+tau*x(i+1)
    end
    temp(n)=(1-tau)*x(n)+tau*x(n-1)
    x=temp
  end
endfunction


// Non-linear Diffusion
// Perona-Malik Diffusivity
// Semi-Implicit scheme
// solve linear system of equations
// please consult lectures on non-linear diffusion schemes in
// Differential Equations in Image Processing and Computer Vision
// and Numerical Algorithms for Visual Computing II for more information

function x = g(a,b,lambda)
s = (a-b)^2
x = 1/(1+s^2/lambda^2)
endfunction


function [B,x] = nonlindiff(x,iter,tau,lambda)
  n = length(x)
  A(1:n,1:n)=0
  // first plot of original signal
  plot2d(1:n,x,style=1)
  for t=1:iter

// build system matrix B
   for i=1:(n-1)
    A(i,i+1)=g(x(i),x(i+1),lambda)
    A(i+1,i)=g(x(i),x(i+1),lambda)
   end

   A(1,1)=-A(1,2)
   for i=2:n-1
    A(i,i)=-A(i,i+1)-A(i,i-1)
 end
   A(n,n)=-A(n,n-1)
   
   // build system matrix (I-tau*A)
   I = eye(n,n)  
   
   B = I-tau*A
   init(1:n) = 0
   
   x= SOR_2(B,x,init)
end
  plot2d(1:n,x,style=2)
endfunction