clear all;

rng('default'); r = trnd(1,1000000,1);
a=-5; b=5; pts=100; stp=(b-a)/pts; c=a:stp:(b-stp); e=([a-stp c]+[c b+stp])/2;
H1 = histc(r,e)/length(r);
figure(1); stem([c c(end)+stp]',H1/max(H1));
hold on; plot(c,1./(1+c.^2),'r-');

a=0; b=10; pts=10000; stp=(b-a)/pts; x=a:stp:(b-stp);
y=(x/5+(1e-3)*((x-1).^2).*((x-5).^2).*((x-9).^2))+0.5*cos(2*pi*4*x); m=min(y);
y=(x/5+(1e-3)*((x-1).^2).*((x-5).^2).*((x-9).^2))+0.5*cos(2*pi*4*x)-m;
figure(2); plot(x,y,'k-');

% SA

rand('state',0); randn('state',0);
epsilon = 0.07; T0 = 0.5; K = 50; N = 1000; fim = 0; 

x = 5; xmin = x;
J = (x/5+(1e-3)*((x-1).^2).*((x-5).^2).*((x-9).^2))+0.5*cos(2*pi*4*x)-m; Jmin = J;
t = 1; T = T0/log2(t+1); n = 0; fim = 0; BKD = zeros(N*K,3);

while not(fim),
    n = n+1; 
    xhat = x + epsilon*randn(1); % Convergência mais lenta com T*epsilon*randn(1)
    Jhat = (xhat/5+(1e-3)*((xhat-1).^2).*((xhat-5).^2).*((xhat-9).^2))+0.5*cos(2*pi*4*xhat)-m;
    if rand(1) < exp((J-Jhat)/T), x = xhat; J = Jhat; end;
    if J<Jmin, Jmin = J; xmin = x; end;
    if rem(n,N) == 0, 
        t = t+1; T = T0/log2(t+1)
    end;
    if t > K, fim = 1; end;
    BKD(n,:) = [20*T J x];
end;

figure(3); plot(BKD); axis([0 5e4 -0.5 10]);

% FSA

rand('state',0); randn('state',0);
epsilon = 0.07; T0 = 0.5; K = 50; N = 1000; fim = 0; 

x = 9; xmin = x;
J = (x/5+(1e-3)*((x-1).^2).*((x-5).^2).*((x-9).^2))+0.5*cos(2*pi*4*x)-m; Jmin = J;
t = 1; T = T0; n = 0; fim = 0; BKD = zeros(N*K,3);

rng('default'); 
while not(fim),
    n = n+1; 
    xhat = x + epsilon*trnd(1); % Convergência mais lenta com T*epsilon*trnd(1)
    Jhat = (xhat/5+(1e-3)*((xhat-1).^2).*((xhat-5).^2).*((xhat-9).^2))+0.5*cos(2*pi*4*xhat)-m;
    if rand(1) < exp((J-Jhat)/T), x = xhat; J = Jhat; end;
    if J<Jmin, Jmin = J; xmin = x; end;
    if rem(n,N) == 0, 
        t = t+1; T = T0/t
    end;
    if t > K, fim = 1; end;
    BKD(n,:) = [20*T J x];
end;

figure(4); plot(BKD); axis([0 5e4 -0.5 10]);
