% Test sequence
M_h=3; N_h=3; M_w=7; N_w=7;
x_1=load_texture(5,128,128);  x_1=x_1-mean(x_1(:)); m_x1=mean(x_1(:));
x_2=load_texture(92,128,128); x_2=x_2-mean(x_2(:)); m_x2=mean(x_2(:));
r_xx1=xcorrf2(x_1); r_xx1=r_xx1((-22:22)+128,(-22:22)+128);
r_xx2=xcorrf2(x_2); r_xx2=r_xx2((-22:22)+128,(-22:22)+128);
h=randn(M_h,N_h); h_=reshape(h.',M_h*N_h,1);
w=ones(7)/7^2;

% Object function gradient: OK
J =fisher_gopt_J(h_,w,m_x1,m_x2,r_xx1,r_xx2,M_h,N_h);
DJ_m=fisher_gopt_DJ(h_,w,m_x1,m_x2,r_xx1,r_xx2,M_h,N_h);
DJ_m=reshape(DJ_m,N_h,M_h).';
DJ=zeros(M_h,N_h);
waitbar(0,'Numerical gradient estimation');
for i=1:M_h
  for j=1:N_h
    dh=zeros(M_h,N_h); dh(i,j)=1e-10; dh_=reshape(dh',M_h*N_h,1);
    J2=fisher_gopt_J(h_+dh_,w,m_x1,m_x2,r_xx1,r_xx2,M_h,N_h);
    DJ(i,j)=(J2-J)./dh(i,j);
    waitbar(100*i/(M_h*N_h));
  end
end
close
disp('Criterion gradient: modelled and numerically estimated:');
disp([DJ_m DJ])

% Mean gradient: OK
[m_1,Dm_1,var_1,Dvar_1]=Dfstats(m_x1,r_xx1,h,w);
Dm=0;
for i=1:M_h
  for j=1:N_h
    dh=zeros(M_h,N_h); dh(i,j)=1e-10;
    m=fstats(m_x1,r_xx1,h+dh,w);
    Dm(i,j)=(m-m_1)./dh(i,j);
  end
end
disp('Mean gradient: modelled and numerically estimated:');
disp([Dm_1, Dm])

% Variance gradient: OK
Dvar=zeros(M_h,N_h);
for i=1:M_h
  for j=1:N_h
    dh=zeros(M_h,N_h); dh(i,j)=1e-10;
    [m,var]=fstats(m_x1,r_xx1,h+dh,w);
    Dvar(i,j)=(var-var_1)./dh(i,j);
  end
end
disp('Variance gradient: modelled and numerically estimated:');
disp([Dvar_1, Dvar])

 
 
 