% Przyklad sieci jednokierunkowej dwu-warstwowej trenowanej algorytmem
% propagacji wstecznej. 
% Siec ma: 
% 	R - wejsc,
%	S1 - neuronow ukrytych,
%	S2 - neuronow wyjsciowych.
% Neurony ukryte maja funkcje aktywacji tangensa hiperbolicznego.
% Neurony wyjsciowe sa liniowe. 
% Siec jest trenowana dla Q wektorow.
% Uwaga: Wektory interpretowane sa jako kolumny.


% DEFINICJA PROBLEMU 
%===================

% Macierz wektorow wejsciowych jest zapisana w P.
% Macierz wektorow wyjsciowych jest zapisana w T.

% INICJALIZACJA WAG POLACZEN SIECI 
%================================

% Ustalenie rozmiarow wektorow wejsciowych i wyjsciowych.

[R,Q] = size(P); [S2,Q] = size(T);	/ [liczba_wierszy, liczba_kolumn] = size(macierz);

% Ustalenie rozmiaru warstwy ukrytej (zmienna S1).

S1 = ?;		/ Mozna przyjac srednia arytmetyczna z R i S2;

% Inicjalizacja wag i przesuniec.

W10 - S1 x R - macierz wag polaczen z warstwy wejsciowej do ukrytej.
B10 - S1 x 1 - macierz przesuniec neuronow warstwy ukrytej.

W20 - S2 x S1 - macierz wag polaczen z warstwy ukrytej do wyjsciowej.
B20 - S2 x 1 - macierz przesuniec neuronow warstwy wyjsciowej.

[W10, B10] = rands(S1, R);	/ rands() - losowa inicjalizacja wag w zakresie (-1, 1);
[W20, B20] = rands(S2, S1);

% TRENING SIECI
%=============

% Parametry treningu

czestotliwosc-wyswietlania = 10;
max_liczba_epok = 500;
blad_docelowy = 0.02;
wsp_uczenia = 0.00008;

% Parametry sieci

W1 = W10; B1 = B10;	/ wartosci poczatkowe odpowiednich macierzy sa przepisywane do 
W2 = W20; B2 = B20;	/ W1, W2, B1, B2

% Faza prezentacji

A1 = tansig(W1*P,B1);		/ obliczanie aktywacji w warstwie ukrytej.
A2 = purelin(W2*A1,B2);		/ obliczanie aktywacji w warstwie wyjsciowej.

E = T-A2;			/ obliczanie bledu.
SSE = sumsqr(E);		/ obliczanie bledu srednio-kwadratowego (Sum Squared Error - SSE).

% WLASCIWY TRENING

TR = [SSE];	/ w macierzy TR - n x 1 - zapisywany jest przebieg treningu.

for (epok = 1; epok<=max_liczba_epok; epok++)
{
% Sprawdzenie wartosci bledu
if (SSE < blad_docelowy) {epok = epok-1; break; } / zakonczenie uczenia

% Faza uczenia

D2 = deltalin(A2,E);		/ dotyczy warstwy wyjsciowej
D1 = deltatan(A1,D2,W2);	/ dotyczy warstwy ukrytej

[dW1,dB1] = learnbp(P,D1,wsp_uczenia);	/ obliczenie wymaganych zmian wag
[dW2,dB2] = learnbp(A1,D2,wsp_uczenia);

W1 = W1 + dW1; B1 = B1 + dB1;
W2 = W2 + dW2; B2 = B2 + dB2;

A1 = tansig(W1*P,B1);		/ obliczanie aktywacji w warstwie ukrytej.
A2 = purelin(W2*A1,B2);		/ obliczanie aktywacji w warstwie wyjsciowej.

E = T-A2;			/ obliczanie bledu.
SSE = sumsqr(E);		/ obliczanie bledu SSE.

% Zapisanie przebiegu treningu w macierzy TR

TR = [TR SSE];

% Wyswietlanie przebiegu treningu

if (rem(epok,czestotliwosc_wyswietlania) == 0)	/ rem() - reszta z dzielenia
    printf('epok=%.0f  SSE=%g\n\n',epok,SSE);

}
________________________________________________


% WYSWIETLENIE WYNIKU TRENINGU
%===============================

% Wykres bledu
%==========

plot(TR);
pause;

% Inne wielkosci
% ==========

SSE = sumsqr(T-purelin(W2*tansig(W1*P,B1),B2));
printf('Liczba epok treningu %.0f \n',epok)
printf('Docelowy blad wynosi %g.\n',blad_docelowy);
printf('Koncowy blad SSE %g.\n',SSE);

______________________________________________________________

function d = deltalin(a,d,w)
%DELTALIN Funkcja obliczajaca delte dla neuronow liniowych
%         DELTALIN(A,E)
%           A - S1xQ macierz wektorow wyjsciowych
%           E - S1xQ macierz bledow
%         Zwraca macierz SxQ wartosci delty
%
%         DELTALIN(A,D,W)
%           D - S2xQ macierz wektorow delta w poprzedniej warstwy
%           W - S2xS1 macierz wag polaczen pomiedzy warstwami
%         Zawraca macierz SxQ wartosci delty dla warstwy ukrytej

if nargin == 2	/ liczba argumentow wejsciowych
  d = d;
else
  d = w'*d;	/ ' oznacza transpozycje
end
_____________________________________________________
function d = deltalog(a,d,w)
%DELTALOG Funkcja obliczajaca delte dla neuronow logistycznych
%         
%         DELTALOG(A,E)
%           A - S1xQ matrix of output vectors
%           E - S1xQ matrix of associated errors
%         Returns an SxQ matrix of output layer delta vectors.
%
%         DELTALOG(A,D,W)
%           D - S2xQ matrix of next layer delta vectors
%           W - S2xS1 weight matrix between layers.
%         Returns an SxQ matrix of hidden layer delta vectors.

if nargin == 2
  d = a.*(1-a).*d;
else
  d = a.*(1-a).*(w'*d);
end

__________________________________________________
function d = deltatan(a,d,w)
%DELTATAN Funkcja obliczajaca delte dla neuronow tangensa-hiperbolicznego
%         
%         DELTATAN(A,E)
%           A - S1xQ matrix of output vectors
%           E - S1xQ matrix of associated errors
%         Returns an SxQ matrix of output layer delta vectors.
%
%         DELTATAN(A,D,W)
%           D - S2xQ matrix of next layer delta vectors
%           W - S2xS1 weight matrix between layers.
%         Returns an SxQ matrix of hidden layer delta vectors.

if nargin == 2
  d = (ones-(a.*a)).*d;	/ ones oznacza 1; ones(N, M) oznacza macierz liczb 1 wymiaru NxM
else
  d = (ones-(a.*a)).*(w'*d);
end
__________________________________________________
function s = sumsqr(a)
%SUMSQR   Oblicza sume kwadratow elementow macierzy a.
%
%         SUMSQR(A)
%           A - a matrix.
%         Returns the sum of squared elements in A.

s = sum(sum(a.*a));	/ a.*a oznacza mnozenie tablic
			/ sum oznacza sumowanie - dla macierzy - wszystkich elemntow w wierszu
			/ dla kolumny - wszystkich elementow w kolumnie
_________________________________________________
function [dw,db] = learnbp(p,d,lr)
%LEARNBP Regula uczenia propagacji wstecznej:
%           dW(i,j) = wsp_uczenia*D(i)*P(j)
%
%         LEARNBP(P,D,LR)
%           P  - RxQ macierz wektorow wejsciowych
%           E  - SxQ macierz wektorow bledow
%           LR - wspolczynnik uczenia
%         Zwraca macierz zmian wag polaczen
%
%         [dW,dB] = LEARNBP(P,D,LR)
%         Zwraca:
%           dW - macierz zmian wag polaczen
%           dB - wektor zmian przesuniec

x = lr*d;
dw = x*p';
if nargout == 2
  [pr,pc] = size(p);
  db = x*ones(pc,1);
  end
_________________________________________________
function a = logsig(n,b)
%LOGSIG   Logistyczna funkcja aktywacji
% 
%         LOGSIG(N)
%           N - SxQ macierz wejsciowa
%         Zwraca obliczona funkcje logistyczna dla N bez przesuniecia
%
%         LOGSIG(N,B) 
%           B - Sx1 wektor przesuniec
%         Zwraca obliczona funkcje logistyczna dla N z przesunieciem (dla kazdego wiersza N 
%           odpowiedni wiersz B.

if nargin==1
  z = n;
else
  [nr,nc] = size(n);
  z = n + b*ones(1,nc);
  end

  a = 1 ./ (1+exp(-z));

_______________________________________________
function a = purelin(n,b)
%PURELIN  Liniowa funkcja aktywacji
% 
%         PURELIN(N)
%           N - SxQ matrix of net input vectors.
%         Returns N.
%
%         PURELIN(N,B) 
%           B - Sx1 vector of biases.
%         Returns the N with the bias in each
%           row of B added to the elements in
%           corresponding rows of N.

if nargin==1
  a = n;
else
  [nr,nc] = size(n);
  a = n+b*ones(1,nc);
  end
__________________________________________________________
function a = tansig(n,b)
%TANSIG   Funkcja aktywacji tangensa hiperbolicznego
% 
%         TANSIG(N)
%           N - SxQ matrix of net input vectors.
%         Returns the tan-sigmoid with bias of 0.
%
%         TANSIG(N,B)
%           B - Sx1 vector of biases.
%         Returns the tan-sigmoid using the biases
%           in each row of B for elements in
%           corresponding rows of N.
%

if nargin==1
   z = n;
else
  [nr,nc] = size(n);
  z = n + b*ones(1,nc);
  end

a = 2 ./ (1 + exp(-2*z)) - 1;