function [val,dval_dx] = post_processing(valin,dvalin_dx,model,net,pp)
if nargout==2
  calcder = true;
  dval_dx = dvalin_dx;
else
  calcder = false;
end
[~,nOutputs] = size(valin);
rhogauss = model.rgrid.rhogauss;

%% Clipping to minimum and maximum
% clip values to minimum/maximum (smoothly to avoid jumps in derivative)
if strcmp(net.name,'QLKNN4Dkin')
  % Threshold Matching was required only in 4Dkin network
  imposeThresholdMatching = net.imposeThresholdMatching;
else
  imposeThresholdMatching = false(1,size(valin,2));
end

[lammin,dlammin_dval] = apply_clipping(valin,pp.minTCoeff,'min',imposeThresholdMatching);
valmin = lammin.*valin + bsxfun(@times,(1-lammin),pp.minTCoeff);
[lammax,dlammax_dval] = apply_clipping(valmin,pp.maxTCoeff,'max',imposeThresholdMatching);
valmax = lammax.*valmin + bsxfun(@times,(1-lammax),pp.maxTCoeff);
val = valmax;

%% Core and edge patch
iedge = rhogauss>pp.rhoedge;
switch pp.edgeStyle
  case 'interpolate'
    lamedge = 1-(rhogauss(iedge)-pp.rhoedge).^2/(1-pp.rhoedge).^2;
    edgeValues = pp.edgeValues;
    
  case 'constant'
    lamedge = zeros(sum(iedge),1);
    edgeValues = pp.edgeValues;
    
  case 'hold'
    lamedge = zeros(sum(iedge),1);
    ilast = find(~iedge,1,'last'); % last value still in core
    edgeValues = repmat(val(ilast,:),sum(iedge),1); % last value;
  otherwise
    error('invalid method');
end
val(iedge,:)  = bsxfun(@times,lamedge,val(iedge,:))  + bsxfun(@times,(1-lamedge) , edgeValues);

% core
icore = rhogauss<pp.rhocore;
switch pp.coreStyle
  case 'interpolate'
    lamcore = 1-(rhogauss(icore)-pp.rhocore).^2/(pp.rhocore).^2;
    coreValues = pp.coreValues;
  case 'constant'
    lamcore = zeros(sum(icore),1);
    coreValues = pp.coreValues;
  case 'hold'
    lamcore = zeros(sum(icore),1);
    ifirst = find(~icore,1,'first'); % first value still in core
    coreValues = val(ifirst,:); % last value;
  case 'gaussian'
    c = 0; % Force centered at rho=0
    H = pp.coreValues;
    w = pp.rhocore / 2; % Derive width from boundary of core. 3 sigma ~ 99.7%
    blendw = 0.1;
    iblend = icore & rhogauss > pp.rhocore - blendw;
    lamcore = (rhogauss(iblend)-pp.rhocore + blendw).^2 ./ (blendw).^2;
    lamcore = [zeros([sum(icore) - sum(iblend) 1]); lamcore];
    coreValues = bsxfun(@times, H, exp(-((rhogauss(icore) - c)/w).^2)) + val(icore,:); % additive
  otherwise
    error('invalid method');
end
val(icore,:) = bsxfun(@times,lamcore,val(icore,:)) + bsxfun(@times,(1-lamcore), coreValues);

%% smoothing
if strcmp(model.hmode.modeltype, 'imposed')
  % Do not smooth over pedestal
  iped = model.rgrid.rhogauss > model.hmode.rhoped;
  out_sm_kernel = get_smkernel(rhogauss(~iped),pp.outSmoothing);
else
  out_sm_kernel = get_smkernel(rhogauss,pp.outSmoothing);
end
valpatch = val;
val = out_sm_kernel*valpatch;
%%
if calcder
  for iOut = 1:nOutputs
    % must apply chain rule in order of operations on original data
    
    % clip minimum
    dlammin_dx = bsxfun(@times,dlammin_dval(:,iOut),dvalin_dx(:,:,iOut));
    % d(lammin.*valin)/dx
    dval_dx_i = bsxfun(@times,lammin(:,iOut),dvalin_dx(:,:,iOut)) + ...
      bsxfun(@times,valin(:,iOut),dlammin_dx);
    % d((1-lammin) * minOutput)/dx
    dval_dx_i = dval_dx_i - pp.minTCoeff(iOut) * dlammin_dx;
    
    % clip maximum
    dlammax_dx = bsxfun(@times,dlammax_dval(:,iOut),dval_dx_i);
    dval_dx_i = bsxfun(@times,lammax(:,iOut),dval_dx_i) + ...
      bsxfun(@times,pp.maxTCoeff(iOut)-valmin(:,iOut),dlammax_dx);
    
    % core patch
    switch pp.coreStyle
      case 'hold'
        dvalcore_dx_i = dval_dx_i(ifirst,:);
      case 'gaussian'
        dvalcore_dx_i = dval_dx_i(icore, :);
      otherwise
        dvalcore_dx_i = zeros(1,model.dims.nx);
    end
    dval_dx_i(icore,:) = bsxfun(@times,lamcore,dval_dx_i(icore,:)) + ...
      bsxfun(@times,1-lamcore,dvalcore_dx_i);
    
    % edge patch
    switch pp.edgeStyle
      case 'hold'
        dvaledge_dx_i = dval_dx_i(ilast,:);
      otherwise
        dvaledge_dx_i = zeros(1,model.dims.nx);
    end
    dval_dx_i(iedge,:)  = bsxfun(@times,lamedge,dval_dx_i(iedge,:)) + ...
      bsxfun(@times,1-lamedge,dvaledge_dx_i);
    
    % smoothing
    dval_dx_i = out_sm_kernel*dval_dx_i;
    
    dval_dx(:,:,iOut) = dval_dx_i;
  end
end

if pp.doplot
  plot_postprocessing(val,valpatch,valmax,valin,model)
end
end

function plot_postprocessing(val,valpatch,valmax,valin,model)
nplot = size(val,2);
nrow = 3;
ncol = ceil(nplot/nrow);

figname = 'QLKNN postprocessing debug plot';
persistent haxpp
if isempty(haxpp) || ~all(all(ishandle(haxpp)))
  hf = findall(0,'Name',figname,'Type','figure');
  if isempty(hf) || ~ishandle(hf)
      hf = figure('Name',figname,'NumberTitle','off','handlevisibility','callback');
  end
  clf(hf);
  haxpp = reshape(multiaxes(hf,ncol,nrow),ncol,nrow)';
end

titstr = {'\chi_e','\chi_i','De','Ve','Di?','Vi?'};
for ii=1:nplot
  ax = subplot(nrow,ncol,ii,'parent',hf);
  hp=plot(ax,model.rgrid.rhogauss(1:numel(val(:,ii))),[val(:,ii),valpatch(:,ii),valmax(:,ii),valin(:,ii)]);
  if ii==1
    legend(hp,{'smoothed','patched','clipped','input'},'location','northwest','box','off');
  end
  title(ax,titstr{ii})
end
drawnow
end

function [lam,dlam_dval] = apply_clipping(values,clipvalues,minmax,iThresholdMatching_network)
minsm = 10; % 1/smoothing for clipping to minimum
maxsm = 10; % 1/smoothing for clipping to maximum
[nEval,nOutputs] = size(values);

% compute all interpolation coefficients for clipping
switch minmax
  case 'min'
    smsm = minsm;
    valexp = min(exp(smsm*(bsxfun(@minus,clipvalues,values))),1/eps.^2);
    %iThresholdMatching = iThresholdMatching_network;
    iThresholdMatching = false(1,nOutputs); % temp turn off
  case 'max'
    smsm = maxsm;
    valexp = min(exp(smsm*(bsxfun(@minus,values,clipvalues))),1/eps.^2); % set maximum for numerics
    iThresholdMatching = false(1,nOutputs); % temp turn off
end
% individual interpolation lambdas for each output
lam = 1./(1+valexp); % close to 0 when output lower than minimum or when larger than max
dlam_dval = (smsm.*valexp.*lam.^2);

if any(iThresholdMatching)
  %% TURNED OFF TEMPORARILY, REMOVE IF UNNECESSARY
  % multiply lambdas to ensure same clipping to minimum for all channels.
  % (to ensure matching the threshold)
  % take product over individual lam values for each output
  lamTM = lam(:,iThresholdMatching);
  lamcombined = prod(lamTM,2);
  nImpose = sum(int8(iThresholdMatching));
  
  dlamcombined_dval = zeros(nEval,1);
  dlamTM = dlam_dval(:,iThresholdMatching);
  for iImp=1:nImpose
    der_thisOut = dlamTM(:,iImp);
    prod_others = prod(lamTM(:,setdiff(1:nImpose,iImp)),2); % product of others
    dlamcombined_dval = dlamcombined_dval + (der_thisOut.*prod_others); % product rule
  end
  lam(:,iThresholdMatching) = repmat(lamcombined,1,nImpose);
  dlam_dval(:,iThresholdMatching) = repmat(dlamcombined_dval,1,nImpose); % substitute for those affected
end

end
