function  [muY,sigY,table] = squashGauss(muX,sigX,table)

%SQUASHGAUSS -- Mean and std.dev. of a truncated Gaussian variable.
%
%  [muY,sigY] = squashGauss(muX,sigX)
%  [muY,sigY] = squashGauss(muX,sigX,table)
%
%  Let X is a (matrix of) Gaussian variable(s) (cf. NORMPDF, NORMCDF,
%  and ERF) with mean(s) MUX and standard deviation(s) SIGX.
%  Construct a new variable Y = max(0,sigmoid(X)), where
%    sigmoid(X) = 1./(1+exp(-X))-0.5
%  is the sigmoidal function with gain 1.0 and asymptotic levels +/-0.5.
%  The new variable has mean 0<=MUY<=0.5 and std. deviation SIGY<=SIGX.
%  SQUASHGAUSS tabulates MUY and SIGY as a function of MUX and SIGX.
%  Typically, MUX is a row vector and SIGX is a scalar, in which case
%  both MUY and SIGY are row vectors of the same length as MUX.
%
%  The values are interpolated from a table built via Monte Carlo methods.
%  To generate the table, use the call
%    [muY,sigY,table] = squashGauss(muX,sigX,N) ,
%  where MUX defines the mu-grid (e.g. [0:0.05:10]), SIGX defines the
%  sigma-grid (e.g. [.1:.1:10]'), and N is an integer specifying the
%  number of samples. (MU- and SIG- grids must satisfy INTERP2 argument
%  conventions. Typically MU_GRID is a row vector and SIG_GRID a col-v.)
%  The TABLE is a structure with five fields: 'mu_grid', 'sig_grid',
%  'muY', 'sigY', and 'N'. It is stored in a persistent variable within
%  the local workspace so subsequent calls to SQUASHGAUSS can re-use it.
%  Alternatively, TABLE can be explicitly passed to the function or
%  loaded from a file. In this latter case, the second argument must
%  be a string indicating a filename. The default is 'squashGauss.mat'.
%
%  See also TRUNCGAUSS, INTERP2, RANDN, NORMINV.

%  Original coding by Alexander Petrov, UC Irvine.
%  $Revision: 1.0 $  $Date: 2002/08/05 12:32 $

persistent  SQUASHGAUSS_TABLE

%-- Make sure there is a valid table, one way or another
rebuild_p = 0 ;
if (nargin==2)    % nothing is supplied explicitly; the typical case
   if (isempty(SQUASHGAUSS_TABLE))
      disp('Loading the SQUASHGAUSS table from squashGauss.mat ...') ;
      load('squashGauss.mat') ;  % load SQUASHGAUSS_TABLE from disk
   %else the persistent variable is already set and everything is OK
   end
elseif (ischar(table))           % file name?
   disp(['Loading the SQUASHGAUSS table from ',table,' ...']) ;
   load(table) ;
elseif (isstruct(table))         % explicit structure?
   SQUASHGAUSS_TABLE = table ;   % update the persistent variable
elseif  (isnumeric(table))       % integer N ==> request to (re)build
   rebuild_p = 1 ;
   SQUASHGAUSS_TABLE = build_table(muX,sigX,table) ;
else
   error('TABLE must be a struct, a filename, or an integer.') ;
end
table = SQUASHGAUSS_TABLE ;

%-- Interpolate for the requested values of MUX and SIGX
if (rebuild_p)  % MUX and SIGX are really MU_GRID and SIG_GRID
   muY = table.muY ;
   sigY = table.sigY ;
else  % normal operation -- interpolate
   mu_grid = table.mu_grid ;
   sig_grid = table.sig_grid ;
   %- Check ranges
   if (min(min(muX)) < mu_grid(1) | max(max(muX)) > mu_grid(end) | ...
         min(min(sigX))<sig_grid(1) | max(max(sigX))>sig_grid(end))
      warning('Out of bounds. Please expand the SQUASHGAUSS table.') ;
   end
   %- Canonicalize sizes (as INTERP2 does not quite work as documented):
   [mu_rows,mu_cols] = size(muX) ;
   [sig_rows,sig_cols] = size(sigX) ;
   if (mu_rows==1 & sig_rows>1)  muX=repmat(muX,sig_rows,1) ; end
   if (sig_cols==1 & mu_cols>1)  sigX=repmat(sigX,1,mu_cols) ; end
   if (any(size(muX)~=size(sigX)))
      error('Incompatible sizes of MUX and SIGX.') ;
   end
   %- Now we are finally ready to interpolate
   muY  = interp2(mu_grid,sig_grid,table.muY ,muX,sigX,'*linear') ;
   sigY = interp2(mu_grid,sig_grid,table.sigY,muX,sigX,'*linear') ;
end

%-- Return MUY and SIGY
%%%%% End of the main function


%%%%%%%%%  SUBFUNCTIONS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  table = build_table(mu_grid,sig_grid,N)

disp('Building SQUASHGAUSS table...') ;

X_grid = norminv([1:N-1]./N)' ;  % a large set of Gaussian(0,1) deviates

%-- Canonicalize grids
[mu_rows,mu_cols] = size(mu_grid) ;
[sig_rows,sig_cols] = size(sig_grid) ;
if (mu_rows==1)
   replicate_rows_p = 1 ;
elseif (mu_rows~=sig_rows)
   error('MU_GRID and SIG_GRID should have the same number of rows.') ;
end
if (sig_cols==1)
   replicate_cols_p = 1 ;
elseif (sig_cols~=mu_cols)
   error('MU_GRID and SIG_GRID should have the same number of columns.') ;
end

%-- Main loop
muX = zeros(sig_rows,mu_cols) ;
sigX = zeros(sig_rows,mu_cols) ;
for r = 1:sig_rows
   for c = 1:mu_cols
      if replicate_rows_p  mu =mu_grid(c);  else mu =mu_grid(r,c);  end
      if replicate_cols_p  sig=sig_grid(r); else sig=sig_grid(r,c); end
      X = X_grid .* sig + mu ;      % Gaussian
      X = 1 ./ (1+exp(-X)) - 0.5 ;  % sigmoid transform: gain=1, M=+/-0.5
      Y = max(0,X) ;
      muY(r,c) = mean(Y) ;
      sigY(r,c) = std(Y,1) ;
   end
end

%-- Pack and return
table.mu_grid = mu_grid ;
table.sig_grid = sig_grid ;
table.muY = muY ;
table.sigY = sigY ;
table.N = N ;

%%%%% End of subfunction BUILD_TABLE

%%%%% End of file SQUASHGAUSS.M
