% Poszukiwanie minimum funkcji bledu metoda symulowanego wyzarzania.
% Wybor nowego polozenia z rozkladem rownomiernym, akceptacja nastepuje z
% pewnym rozkladem prawdopodobienstwa zaleznym od wielkosci zmniejszenia
% bledu(dE) i temperatury(T):
%   h(dE,T) = 1/(1 + exp(dE/cT))
% oraz pewnej stalej c. Jesli c mala dodatnia to algorytm redukuje sie
% do szukania przypadkowego (akceptacja tylko modyfikacji zmniejszajacych blad),
% jesli c duza to do bladzenia przypadkowego (poprawka nie zalezy od bledu)

liczba_param = 2;
N = liczba_param;

liczba_cykli = 550000;

Rozw = rand(1,N)*20-10;                 % rozw. poczatkowe

T = 500;                                % pocz. krok iteracji
Tmin = 0.01;                            % minimalny krok iteracji
wT = 0.99999;                           % wsp zminy kroku
c = 0.1;                                % wsp. regulacji wplywu temperatury
e = 0.001;
znak = [-1,1];

Emin = 10e40;
Epop = 0;
format long;                            % 15 miejsc znaczacych
Droga = [Rozw];

tabE = zeros(2,liczba_cykli);

for cyk=1 : liczba_cykli
   %Nowa wartosc x
   Rozw2(1) = Rozw(1) + T*(-1+2*rand)*e;           % nowe rozwiazanie (powinno byc uzaleznine od temperatury)
   %Sprawdzenie czy wychodzi poza zakres
   while abs(Rozw2(1)) > 10
    if Rozw2(1) > 10
       Rozw2(1) = Rozw2(1) + (10 - Rozw2(1));
    end
    if Rozw2(1) < -10
       Rozw2(1) = Rozw2(1) + (-10 - Rozw2(1));
    end
   end

   %Nowa wartosc y
   Rozw2(2) = Rozw(2) + T*(-1+2*rand)*e;
   %Sprawdzenie czy wychodzi poza zakres
   while abs(Rozw2(2)) > 10
    if Rozw2(2) > 10
       Rozw2(2) = Rozw2(2) + (10 - Rozw2(2));
    end
    if Rozw2(2) < -10
       Rozw2(2) = Rozw2(2) + (10 - Rozw2(2));
    end
   end

   E = fun3(Rozw2(1),Rozw2(2));         % wartosc funkcji

   dE = E - Epop;                       % roznica wartosi funkcji
   if rand < 1/(1+exp(dE/(c*T)))        % warunek akceptacji
      Rozw = Rozw2;
      Epop = E;
   end

   if E < Emin                          % sprawdzenie czy jest rekord
      Emin = E;
      Rozwmin = Rozw2 ;
      Droga = [Droga ; Rozwmin];
   end
   tabE(1,cyk) = Emin;
   tabE(2,cyk) = cyk;

   T = T*wT;                            % nowa wartosc temperatury
   if T<Tmin
      T=Tmin;
   end
end

wykres_fun
hold on
for cyk = 1 : length(Droga)
    plot(Droga(1:cyk,1),Droga(1:cyk,2));
    hold on;
    plot(Droga(1:cyk,1),Droga(1:cyk,2),'ro');
    text(Droga(cyk,1)'+0.01, Droga(cyk,2)'+0.01, int2str(cyk), 'fontsize', 8);
    pause(0.1);
    drawnow;
end
hold on

plot(Droga(end,1)',Droga(end,2)','mo','LineWidth',3);
x = Droga(end,1)
y = Droga(end,2)
f = tabE(1,end)

% %plot(tabE(2,:),tabE(1,:),'LineWidth',1.5)