

function [biasgp,ge] = capacity_deep(m,seed)
mall = [3:20,30:10:300];
N_0= 30;   
Ns = 2000;
M0 = mall(m);


P = 6300;
%Nall = fliplr(P./alphaall);

%N_1 = 200;
%sigma = sigmaall(seed1);
sigma = 1;
beta = 1e5;


%%
    load('mnist.mat')
    
    noise = 0.;

    rng(10,'twister');

    %% generating samples from noisy linear teacher neurons
    % random projection weights
trainX = double(trainX);
testX = double(testX);

% trainX = (trainX'-mean(trainX'))./std(trainX');
% testX = (testX'-mean(testX'))./std(testX');

 w0 = normrnd(0,1,N_0,784);


trainX = (trainX'-mean(trainX'))./std(trainX');
testX = (testX'-mean(testX'))./std(testX');
% % 
trainX = (w0*trainX>0).*(1/sqrt(N_0)*w0*trainX);
testX = (w0*testX>0).*(1/sqrt(N_0)*w0*testX);

trainX = (trainX-mean(trainX))./std(trainX);
testX = (testX-mean(testX))./std(testX);


     
     rng(2,'twister');


    indtrain = find(mod(trainY,2)==0);
    indtest = find(mod(testY,2)==0);
    indtrain1 = find(mod(trainY,2)==1);
    indtest1 = find(mod(testY,2)==1);
    index = [indtrain(1:P/2),indtrain1(1:P/2)];
    index_t = [indtest(1:Ns/2),indtest1(1:Ns/2)];

    x_0 = trainX(:,index);
    x_t = testX(:,index_t);

    y = 2*(mod(trainY(index),2)==0)-1;
    y_t = 2*(mod(testY(index_t),2)==0)-1;

    th = 0;


%     rng(0,'twister');
%     iteration = 100;
%     beta0 = 0.1;
%     [g,xc] = softkmeans(M, x_0, beta0, iteration, 3);
%     g = g(1:M,:);
% %g = (g>1/N/2);
% 
%     for i0 = 1:M
%     d(i0,:) = sum((x_t - xc(:, i0)).^2)';
%     expd(i0,:) = exp(-beta0*d(i0,:));
%     end
%     gp = expd./sum(expd);
%     gp = gp(1:M,:);
%gp = (gp>1/N/2);
     rng(seed,'twister');
     th = 0.;
    v = normrnd(0,1,M0,N_0);
    g = (v*x_0>th);
    gp = (v*x_t>th);
    
ind = find(sum(g)==0);
for i = 1:length(find(sum(g)==0))
    field = v*x_0(:,ind(i));
    maxind = find(field==max(field));
    g(maxind,ind(i)) = 1;
end
ind = find(sum(gp)==0);
for i = 1:length(find(sum(gp)==0))
    field = v*x_t(:,ind(i));
    maxind = find(field==max(field));
    gp(maxind,ind(i)) = 1;
end
    
 


% noise = 0;
% seed = 1;
%    iteration = 100;
%     beta0 = 0.1;
%     [g,xc] = softkmeans(M0, x_0, beta0, iteration, 3);
%     g = g(1:M0,:)+ noise*normrnd(0,1,M0,P);
%     %g(g>=1/M) = 1;
% 
% 
%     for i0 = 1:M0
%        d(i0,:) = sum((x_t - xc(:, i0)).^2)';
%        expd(i0,:) = exp(-beta0*d(i0,:));
%     end
%     gp = expd./sum(expd);
%     gp = gp(1:M0,:)+ noise*normrnd(0,1,M0,Ns);
    %gp(gp>=1/M) = 1;
  %v = normrnd(0,1,M0,N_0);
  %g =(v*x_train>0);
%  gp =(v*x>0);
%g = normrnd(0,1,M,P);
%gp = normrnd(0,1,M,N_s);
   % g = normrnd(0,1,M0,P);
   %g = (normrnd(0,1,M0,P)>0);

   % compute K0
   
    K0 = sigma^2/N_0*x_t'*x_t;
    K0g = sigma^2/N_0*x_0'*x_0;
    K00g = sigma^2/N_0*x_t'*x_0;

    for l = 1:3

     Kgp{l} = sigma^l*((gp'*gp)/M0).^l.*K0;
     Kg = sigma^l*((g'*g)/M0).^l.*K0g + 1/beta*eye(P,P);
     Kggp = sigma^l*((gp'*g)/M0).^l.*K00g;
     fmean = Kggp*pinv(Kg)*y';
     fvar = diag(Kgp{l} - Kggp*pinv(Kg)*Kggp');
     
     biasgp(l) = mean((fmean-y_t').^2);
     ge(l) = biasgp(l) + mean(fvar);
     gb(l) = mean((y_t'+1)/2+(-y_t').*erfc(-fmean./sqrt(fvar*2))/2);

     

    end
    
    save(['../capacity_deep/capacity_m',num2str(M0),'seed',num2str(seed),'.mat'],'biasgp','ge','gb');
end