function [cost grad] = fineTuneLRCost(alltheta, netstruct, label, data, lambda)

if nargin < 5; lambda = 0; end; % lambda is l2 regularizer for LR, we don't consider decay for networks
label(label==-1) = 0; % label=double(label);
[net lrtheta] = Rolling(alltheta, netstruct);
assert(length(lrtheta) == length(net(end).b)+1);


[lastact record] = passThruNet(net, data); % feedforward, no W2

data = data'; m = size(data,2);
lrtheta = lrtheta'; % row vector

%% -- COST and gradient for LR ----
[junk prob] = lrRegressionPred(lrtheta',lastact);

cost = -(mean(label.*log(prob)+(1-label).*log(1-prob))) + 0.5*lambda*sum(lrtheta.^2); % cross-entropy + weight decay

% a1 = label.*log(prob); a1(isnan(a1))=0; a1(isinf(a1))=0;
% a2 = (1-label).*log(1-prob); a2(isnan(a2))=0; a2(isinf(a2))=0;
% a3 = 0.5*lambda*sum(lrtheta.^2); %a3(isnan(a3))=0;
% cost = -mean(a1+a2)+a3;
% if isnan(cost); 
%     disp('!!!!!!!!!!!!!');cost = 0; 
%     a=label.*log(prob);
%     b=(1-label).*log(1-prob);
%     [a b prob]
% 
% %     c=0.5*lambda*sum(lrtheta.^2)
%     pause;
% end;
expand_act = [lastact ones(m,1)];
lrtheta_grad = -expand_act'*(label-prob)/m + lambda*lrtheta';

%% -- gradients for deepnets by backprop to each hidden units----
partial = lrtheta(1:end-1);
nLayer = length(netstruct) - 1;
Grads = struct;

% if n layer net then n+1 layer units (n_l = nLayer+1)
nHidden = nLayer+1;
D_out = -((label - prob)*partial)'.*sigmoid_grad(record(nHidden).A);
% Grads(nHidden-1).Wgrad = D_out*record(nHidden-1).A'/m;
% Grads(nHidden-1).bgrad = mean(D_out,2);

upperD = D_out;
for ll = nHidden-1 : -1 : 1        
    Grads(ll).Wgrad = upperD*record(ll).A'/m;
    Grads(ll).bgrad = mean(upperD,2); 
    D = net(ll).W'*upperD.*sigmoid_grad(record(ll).A);       
%     Grads(ll-1).Wgrad = upperD.*record(ll-1).A'/m;
%     Grads(ll-1).bgrad = mean(upperD,2);
    upperD = D;
end
% Ds = zeros(nLayer+1,size(D_last,1),size(D_last,2));
% Ds(nLayer+1,:,:) = D_last;
% grad = [];
% 
% for ll = (nLayer-1): -1 : 2
%     Ds(ll,:,:) = net(ll).W'*Ds(ll+1,:,:).*sigmoid_grad(record(ll).A);
% end
% 

grad=[];
for ll = 1:nLayer   
    grad = [grad; Grads(ll).Wgrad(:); Grads(ll).bgrad(:)];
end
grad = [grad; lrtheta_grad];

%% test -- only one hidden layer, topped by an LR layer
% D2 = -((label - prob)*partial)'.*sigmoid_grad(record(2).A); % not Z since f'(z) = a*(1-a)
% % D1 = net(1).W'*D2.*sigmoid_grad(record;
% Wgrad = D2*record(1).A'/m;
% bgrad = mean(D2,2);
% 
% grad = [Wgrad(:); bgrad(:)];
% grad = [grad; lrtheta_grad];

