function [out,dout_din] = evaluate_4Dkin(net,pp,inputs_constrained)
% inputs_constrained: [nIn x nEval]
% out: matrix of outputs [nEval x nOut];
% dout_din: 3D matrix of Jacobians: [nEval x nIn x nOut]

%%% THIS SHOULD BE GENERALIZED for ND and generic networks
% net. structure should be modified accordingly

if nargout==2
  calcder=true;
else
  calcder=false;
end

%normalize input vector to match input normalization in NN
unorm = bsxfun(@times,bsxfun(@minus,inputs_constrained',net.dispfac),1./net.scalefac);
nrho = size(unorm,2);

%%shift R/LTi thgresholds for threshold matching
shift_b1_eef = pp.shiftOutputs(1)*net.IW_eef(:,4)./net.scalefac(4);
shift_b1_ief = pp.shiftOutputs(2)*net.IW_ief(:,4)./net.scalefac(4);

shift_b1_dfe = pp.shiftOutputs(3)*net.IW_dfe(:,4)./net.scalefac(4);
shift_b1_vce = pp.shiftOutputs(4)*net.IW_vce(:,4)./net.scalefac(4);
shift_b1_vte = pp.shiftOutputs(4)*net.IW_vte(:,4)./net.scalefac(4);

%calculate chie regression neural network piece by piece, for all radii
g_e=bsxfun(@plus,net.IW_eef *unorm,net.b1_eef+shift_b1_eef);
sg_e = sigmo(g_e);
f_e=bsxfun(@plus,net.L1W_eef *sg_e,net.b2_eef);
sf_e = sigmo(f_e);

chieGB = (net.L2W_eef*sf_e + net.b3_eef)';

% if the NN output normalization was to between 0 and 1
if net.prepros_eef==0
  chieGB= chieGB.*(net.max_qlk_eef-net.min_qlk_eef)+ones(nrho,1)*net.min_qlk_eef;
end

% if the NN output normalization was to between -1 a +1
if net.prepros_eef==-1
  chieGB= (chieGB+ones(size(chieGB))).*(ones(nrho,1)*net.max_qlk_eef-ones(nrho,1)*net.min_qlk_eef)/2 + ones(nrho,1)*net.min_qlk_eef';
end

%%%%%%%%
%calculate chii regression neural network piece by piece, for all radii
g_i=bsxfun(@plus,net.IW_ief *unorm,net.b1_ief+shift_b1_ief);
sg_i = sigmo(g_i);
f_i=bsxfun(@plus,net.L1W_ief *sg_i,net.b2_ief);
sf_i = sigmo(f_i);

chiiGB = (net.L2W_ief*sf_i + net.b3_ief)';

% if the NN output normalization was to between 0 and 1
if net.prepros_ief==0
  chiiGB= chiiGB.*(net.max_qlk_ief-net.min_qlk_ief)+ones(nrho,1)*net.min_qlk_ief;
end

% if the NN output normalization was to between -1 a +1
if net.prepros_ief==-1
  chiiGB= (chiiGB+ones(size(chiiGB))).*(ones(nrho,1)*net.max_qlk_ief-ones(nrho,1)*net.min_qlk_ief)/2 + ones(nrho,1)*net.min_qlk_ief';
end

if calcder
  % derivatives w.r.t. NN input parameters
  % evalation by chain rule
  
  %First for chie
  dsgdg = (sg_e+1).^2.*exp(-2*g_e); % d(sig(g))/dg
  dsfdf = (sf_e+1).^2.*exp(-2*f_e); % d(sig(f))/df
  dgdu = bsxfun(@times,net.IW_eef ,1./net.scalefac'); % dg/duphys
  
  dchieGB = zeros(nrho,4);
  for iu=1:4 % loop over NN inputs
    % chain rule
    if net.prepros_eef==0
      dchieGB(:,iu) = (net.max_qlk_eef-net.min_qlk_eef).*net.L2W_eef*...
        bsxfun(@times,dsfdf,net.L1W_eef *bsxfun(@times,dsgdg,dgdu(:,iu)));
    end
    if net.prepros_eef==-1
      dchieGB(:,iu) = (net.max_qlk_eef-net.min_qlk_eef)./2.*net.L2W_eef*...
        bsxfun(@times,dsfdf,net.L1W_eef*bsxfun(@times,dsgdg,dgdu(:,iu)));
    end
    
  end
  
  %%Now for chii
  dsgdg = (sg_i+1).^2.*exp(-2*g_i); % d(sig(g))/dg
  dsfdf = (sf_i+1).^2.*exp(-2*f_i); % d(sig(f))/df
  dgdu = bsxfun(@times,net.IW_ief ,1./net.scalefac'); % dg/duphys
  
  dchiiGB = zeros(nrho,4);
  for iu=1:4 % loop over NN inputs
    % chain rule
    if net.prepros_ief==0
      dchiiGB(:,iu) = (net.max_qlk_ief-net.min_qlk_ief).*net.L2W_ief*...
        bsxfun(@times,dsfdf,net.L1W_ief *bsxfun(@times,dsgdg,dgdu(:,iu)));
    end
    if net.prepros_ief==-1
      dchiiGB(:,iu) = (net.max_qlk_ief-net.min_qlk_ief)./2.*net.L2W_ief*...
        bsxfun(@times,dsfdf,net.L1W_ief*bsxfun(@times,dsgdg,dgdu(:,iu)));
    end  
  end
end

% Particle transport coefficients were (accidentally) trained in SI units
% We must revert back to GB units using the dimensional units employed in the training, and then revert to SI
% using the dimensional units in our current simulation
% NN training set values: Te = 8 keV , a = 1 m, B0 = 3 T, Amain = 2, n = 5e19 [m^-3], R/Ln=2
qel = 1.6e-19;
chiGBNN = (8e3.*qel).^1.5.*sqrt(2.*1.67e-27)./(qel.^2.*3.^2.*1);

%calculate De regression neural network piece by piece, for all radii
g_de = bsxfun(@plus,net.IW_dfe *unorm,net.b1_dfe+shift_b1_dfe);
sg_de = sigmo(g_de);
f_de = bsxfun(@plus,net.L1W_dfe *sg_de,net.b2_dfe);
sf_de = sigmo(f_de);

deSI = (net.L2W_dfe*sf_de + net.b3_dfe)';

% if the NN output normalization was to between 0 and 1
if net.prepros_dfe==0
  deSI = deSI.*(net.max_qlk_dfe-net.min_qlk_dfe)+ones(nrho,1)*net.min_qlk_dfe;
end

% if the NN output normalization was to between -1 a +1
if net.prepros_dfe==-1
  deSI = (deSI+ones(size(deSI))).*(ones(nrho,1)*net.max_qlk_dfe-ones(nrho,1)*net.min_qlk_dfe)/2 + ones(nrho,1)*net.min_qlk_dfe';
end

%set output vectors
deGB = deSI./chiGBNN; %convert from SI to GB with dimensional units used in training

%%%%%%%%
%calculate Ve regression neural network piece by piece, for all radii
g_vte  = bsxfun(@plus,net.IW_vte *unorm,net.b1_vte+shift_b1_vte);
sg_vte = sigmo(g_vte);
f_vte  = bsxfun(@plus,net.L1W_vte *sg_vte,net.b2_vte);
sf_vte = sigmo(f_vte);

vteSI  = (net.L2W_vte*sf_vte + net.b3_vte)';

% if the NN output normalization was to between 0 and 1
if net.prepros_vte==0
  vteSI = vteSI.*(net.max_qlk_vte-net.min_qlk_vte)+ones(nrho,1)*net.min_qlk_vte;
end

% if the NN output normalization was to between -1 a +1
if net.prepros_vte==-1
  vteSI = (vteSI+ones(size(vteSI))).*(ones(nrho,1)*net.max_qlk_vte-ones(nrho,1)*net.min_qlk_vte)/2 + ones(nrho,1)*net.min_qlk_vte';
end

%set output vectors
vteGB = vteSI./chiGBNN.*3;  %VGB=VSI*R/chiGB;

g_vce = bsxfun(@plus,net.IW_vce *unorm,net.b1_vce+shift_b1_vce);
sg_vce = sigmo(g_vce);
f_vce = bsxfun(@plus,net.L1W_vce *sg_vce,net.b2_vce);
sf_vce = sigmo(f_vce);

vceSI = (net.L2W_vce*sf_vce + net.b3_vce)';

% if the NN output normalization was to between 0 and 1
if net.prepros_vce==0
  vceSI = vceSI.*(net.max_qlk_vce-net.min_qlk_vce)+ones(nrho,1)*net.min_qlk_vce;
end

% if the NN output normalization was to between -1 a +1
if net.prepros_vce==-1
  vceSI = (vceSI+ones(size(vceSI))).*(ones(nrho,1)*net.max_qlk_vce-ones(nrho,1)*net.min_qlk_vce)/2 + ones(nrho,1)*net.min_qlk_vce';
end

%set output vectors
vceGB = vceSI./chiGBNN.*3;  %VGB=VSI*R/chiGB;
VeGB = vceGB + vteGB;

if calcder
  % derivatives w.r.t. NN input parameters
  % evalation by chain rule
  
  %First for chie
  dsgdg = (sg_de+1).^2.*exp(-2*g_de); % d(sig(g))/dg
  dsfdf = (sf_de+1).^2.*exp(-2*f_de); % d(sig(f))/df
  dgdu = bsxfun(@times,net.IW_dfe ,1./net.scalefac'); % dg/duphys
  
  ddeGB = zeros(nrho,4);
  for iu=1:4 % loop over NN inputs
    % chain rule
    if net.prepros_dfe==0
      ddeGB(:,iu) = ((net.max_qlk_dfe-net.min_qlk_dfe).*net.L2W_dfe*...
        bsxfun(@times,dsfdf,net.L1W_dfe *bsxfun(@times,dsgdg,dgdu(:,iu))) )./chiGBNN;
    end
    if net.prepros_dfe==-1
      ddeGB(:,iu) = ((net.max_qlk_dfe-net.min_qlk_dfe)./2.*net.L2W_dfe*...
        bsxfun(@times,dsfdf,net.L1W_dfe*bsxfun(@times,dsgdg,dgdu(:,iu))))./chiGBNN;
    end
    
  end
  
  %Now for vte
  dsgdg = (sg_vte+1).^2.*exp(-2*g_vte); % d(sig(g))/dg
  dsfdf = (sf_vte+1).^2.*exp(-2*f_vte); % d(sig(f))/df
  dgdu = bsxfun(@times,net.IW_vte ,1./net.scalefac'); % dg/duphys
  
  dvteGB = zeros(nrho,4);
  for iu=1:4 % loop over NN inputs
    % chain rule
    if net.prepros_vte==0
      dvteGB(:,iu) = ((net.max_qlk_vte-net.min_qlk_vte).*net.L2W_vte*...
        bsxfun(@times,dsfdf,net.L1W_vte *bsxfun(@times,dsgdg,dgdu(:,iu))) )./chiGBNN.*3;
    end
    if net.prepros_vte==-1
      dvteGB(:,iu) = ((net.max_qlk_vte-net.min_qlk_vte)./2.*net.L2W_vte*...
        bsxfun(@times,dsfdf,net.L1W_vte*bsxfun(@times,dsgdg,dgdu(:,iu))))./chiGBNN.*3;
    end
    
  end
  
  %Now for vce
  dsgdg = (sg_vce+1).^2.*exp(-2*g_vce); % d(sig(g))/dg
  dsfdf = (sf_vce+1).^2.*exp(-2*f_vce); % d(sig(f))/df
  dgdu = bsxfun(@times,net.IW_vce ,1./net.scalefac'); % dg/duphys
  
  dvceGB = zeros(nrho,4);
  for iu=1:4 % loop over NN inputs
    % chain rule
    if net.prepros_vce==0
      dvceGB(:,iu) = ((net.max_qlk_vce-net.min_qlk_vce).*net.L2W_vce*...
        bsxfun(@times,dsfdf,net.L1W_vce *bsxfun(@times,dsgdg,dgdu(:,iu))) )./chiGBNN.*3;
    end
    if net.prepros_vce==-1
      dvceGB(:,iu) = ((net.max_qlk_vce-net.min_qlk_vce)./2.*net.L2W_vce*...
        bsxfun(@times,dsfdf,net.L1W_vce*bsxfun(@times,dsgdg,dgdu(:,iu))))./chiGBNN.*3;
    end
  end
  dveGB = dvteGB + dvceGB;
end

out = [chieGB,chiiGB,deGB,VeGB];
%fprintf('[chie,chii,de,ve] = %3.3e,%3.3e,%3.3e,%3.3e \n', out)

if calcder
  %fprintf('dchieGB = %3.3e,%3.3e,%3.3e,%3.3e \n', dchieGB)
  %fprintf('dchiiGB = %3.3e,%3.3e,%3.3e,%3.3e \n', dchiiGB)
  %fprintf('dchiiGB = %3.3e,%3.3e,%3.3e,%3.3e \n', ddeGB)
  %fprintf('dchiiGB = %3.3e,%3.3e,%3.3e,%3.3e \n', dveGB)
  
  % Jacobians should be [nEval x nOutputs x nInputs]
  dout_din = assemble_Jac(dchieGB,dchiiGB,ddeGB,dveGB);
end

return


%nonlinear neuron function
function s = sigmo(x)
s= 2./(1+exp(-2.*x))-1; %define sigmoid function (nonlinear transfer functions in NN)
return

function Jac = assemble_Jac(dchieGB,dchiiGB,ddeGB,dveGB)
% assemble Jacobian following convention [nEval x nOut x nIn]
Jac = [resh(dchieGB),resh(dchiiGB),resh(ddeGB),resh(dveGB)];
return

function D=resh(A)
  % reshape [n x m] into [n x 1 x m]
D = reshape(A,size(A,1),1,size(A,2));
return