Index: issm/trunk/src/m/utils/DACE/Contents.m
===================================================================
--- issm/trunk/src/m/utils/DACE/Contents.m	(revision 1080)
+++ issm/trunk/src/m/utils/DACE/Contents.m	(revision 1080)
@@ -0,0 +1,36 @@
+%   DACE (Design and Analysis of Computer Experiments) Toolbox
+%   Version 2.5, September 4, 2002
+%   Copyright (c) 2002 by Hans Bruun Nielsen and IMM. 
+% 
+%   Model construction
+%     dacefit - Constrained non-linear least-squares fit of a given
+%               correlation model to the provided data set and
+%               regression model.
+% 
+%   Model prediction
+%     predictor - Model predictor with mean squared error estimate.
+% 
+%   Regression functions
+%     regpoly0  - Zero order polynomial.
+%     regpoly1  - First order polynomial.
+%     regpoly2  - Second order polynomial.
+% 
+%   Correlation functions
+%     corrcubic     - Local support, cubic polynomial
+%     correxp       - Exponential.
+%     correxpg      - General exponential.
+%     corrgauss     - Gaussian.
+%     corrlin       - Local support, linear.
+%     corrspherical - Local support, spherical.
+%     corrspline    - Local support, cubic spline.
+% 
+%   Experimental Design
+%     gridsamp - Points in a regular grid.
+%     lhsamp   - Latin hypercube distributed random numbers.
+%
+%   Auxiliary functions
+%     dsmerge  - Merge data for multiple design sites.
+%     
+%   Data files
+%     data1.mat - Example data S and Y 
+
Index: issm/trunk/src/m/utils/DACE/changelog
===================================================================
--- issm/trunk/src/m/utils/DACE/changelog	(revision 1080)
+++ issm/trunk/src/m/utils/DACE/changelog	(revision 1080)
@@ -0,0 +1,125 @@
+
+ c :  comment
+ * :  bug fixed
+ + :  new feature/function
+ x :  changed functionality
+
+--------------------------------------------------------------------
+Version 2.5, September 4, 2002
+
+ *  REGPOLY2  did not work properly for dimension greater than 2
+    Bug fixed.
+
+--------------------------------------------------------------------
+Version 2.4, September 3, 2002
+
+ x  DACEFIT  gives an error return in case of a rank deficient matrix F.  
+    May arise with a poor combination of choice of regression model 
+    and distribution of design sites.
+
+--------------------------------------------------------------------
+Version 2.3, August 30, 2002
+
+ x  DACEFIT  gives an error return in case of zero columns in matrix F.  
+    May arise with a poor combination of choice of regression model 
+    and distribution of design sites.
+
+--------------------------------------------------------------------
+Version 2.2, August 27, 2002
+
+ x  DACEFIT  allows design data with standard deviation zero 
+    in one or more directions.
+
+--------------------------------------------------------------------
+Version 2.1, August 26, 2002
+
+ *  PREDICTOR  returned erroneous mse-values in case of single 
+    response, multiple entry.  This bug is fixed.
+
+--------------------------------------------------------------------
+Version 2.0, August 5, 2002
+
+ x  The element  dmodel.regress  is renamed to  dmodel.regr
+    in  DACEFIT, lines 19, 106
+    in  PREDICTOR, lines 58, 90, 107
+    
+ c  Documentation is adapted to the new version. A separate report
+    documenting the algorithm and other special parts of the toolbox
+    is available.
+
+--------------------------------------------------------------------
+Version 1.3, July 3, 2002
+
+ *  A developer debug statement in PREDICTOR line 52:
+
+      save NX x
+
+    is removed.
+
+ +  DSMERGE is a new function for merging data sets with multiple
+    design sites (i.e., sites sampled more than one time).
+    See 'help dsmerge' for more information.
+
+ +  PREDICTOR now return Jacobian matrices for model and mean
+    squared error also for cases with multiple response functions.
+
+ c  We would like to acknowledge the comments and suggestions from
+    Thierry Dalon, from Siemens VDO Automotive AG, which is a great
+    help in improving this software.
+
+--------------------------------------------------------------------
+Version 1.2, June 25, 2002
+
+ +  DACEFIT now has a second output parameter
+
+        [dmodel, perf] = dacefit(S, Y, regr, corr, theta0)
+
+    or        
+
+        [dmodel, perf] = dacefit(S, Y, regr, corr, theta0, lob, upb)
+        
+    The second output parameter gives the convergence history of the
+    optimization - see 'help dacefit' for more information.
+
+ +  DACEFIT now allows the second input argument to be an m*q matrix,
+    i.e., allows q response functions to be modelled, the elements of
+    dmodel: beta, gamma, sigma2, Ysc are changed accordingly.
+    
+ c  The optimization method in DACEFIT is changed. A separate report
+    documenting the algorithm and other special parts of the toolbox
+    will be available in August 2002.
+    
+ +  PREDICTOR is adapted to the changes in DACEFIT, though in case
+    of multiple response functions it does not yet return the
+    Jacobian matrices for the predicted models and their respective
+    mean squared errors.
+    
+ +  CORRCUBIC is the new name of the function named CORRSPLINE
+    in version 1.1. See 'help corrcubic' for more information.
+    
+ x  CORRSPLINE gives a correlation function based on a cubic
+    spline. See 'help corrspline' for more information.
+    The above mentioned report will also discuss this new 
+    correlation model.
+    
+ x  The  parameter vector q in GRIDSAMP is changed so that now it 
+    is the number of equidistant points in each coordinate direction.
+    This is in accordance with the similar parameter in LHSAMP
+    and LINSPACE.
+
+--------------------------------------------------------------------
+Version 1.1, June 7, 2002
+	
+ *  A sign error in the gradient calculation in CORRGAUSS is
+    corrected: line 33 in the file is changed from
+   
+      dr = repmat(-2*theta(:).',m,1) .* d .* repmat(r,1,n);
+
+    to
+
+      dr = -2*repmat(-theta(:).',m,1) .* d .* repmat(r,1,n);
+
+--------------------------------------------------------------------
+Version 1.0, April 16, 2002
+
+ c  First version released
Index: issm/trunk/src/m/utils/DACE/corrcubic.m
===================================================================
--- issm/trunk/src/m/utils/DACE/corrcubic.m	(revision 1080)
+++ issm/trunk/src/m/utils/DACE/corrcubic.m	(revision 1080)
@@ -0,0 +1,42 @@
+function  [r, dr] = corrcubic(theta, d)
+%CORRCUBIC  Cubic correlation function,
+%
+%           n
+%   r_i = prod max(0, 1 - 3(theta_j*d_ij)^2 + 2(theta_j*d_ij)^3) ,  i = 1,...,m
+%          j=1
+%
+% If length(theta) = 1, then the model is isotropic:
+% all theta_j = theta.
+%
+% Call:    r = corrcubic(theta, d)
+%          [r, dr] = corrcubic(theta, d)
+%
+% theta :  parameters in the correlation function
+% d     :  m*n matrix with differences between given data points
+% r     :  correlation
+% dr    :  m*n  matrix with the Jacobian of r at x. It is
+%          assumed that x is given implicitly by d(i,:) = x - S(i,:), 
+%          where S(i,:) is the i'th design site. 
+
+% hbn@imm.dtu.dk  
+% Last update June 25, 2002
+
+[m n] = size(d);  % number of differences and dimension of data
+if  length(theta) == 1
+  theta = repmat(theta,1,n);
+elseif  length(theta) ~= n
+  error(sprintf('Length of theta must be 1 or %d',n))
+else
+  theta = theta(:).';
+end
+td = min(abs(d) .* repmat(theta,m,1), 1);
+ss = 1 - td.^2 .* (3 - 2*td);
+r = prod(ss, 2);
+
+if  nargout > 1
+  dr = zeros(m,n);
+  for  j = 1 : n
+    dd = 6*theta(j) * sign(d(:,j)) .* td(:,j) .* (td(:,j) - 1);
+    dr(:,j) = prod(ss(:,[1:j-1 j+1:n]),2) .* dd;
+  end
+end
Index: issm/trunk/src/m/utils/DACE/correxp.m
===================================================================
--- issm/trunk/src/m/utils/DACE/correxp.m	(revision 1080)
+++ issm/trunk/src/m/utils/DACE/correxp.m	(revision 1080)
@@ -0,0 +1,38 @@
+function  [r, dr] = correxp(theta, d)
+%CORREXP  Exponential correlation function
+%
+%           n
+%   r_i = prod exp(-theta_j * |d_ij|)
+%          j=1
+%
+% If length(theta) = 1, then the model is isotropic: 
+% theta_j = theta(1), j=1,...,n
+%
+% Call:    r = correxp(theta, d)
+%          [r, dr] = correxp(theta, d)
+%
+% theta :  parameters in the correlation function
+% d     :  m*n matrix with differences between given data points
+% r     :  correlation
+% dr    :  m*n matrix with the Jacobian of r at x. It is
+%          assumed that x is given implicitly by d(i,:) = x - S(i,:), 
+%          where S(i,:) is the i'th design site. 
+
+% hbn@imm.dtu.dk  
+% Last update April 12, 2002
+
+[m n] = size(d);  % number of differences and dimension of data
+lt = length(theta);
+if  lt == 1,  theta = repmat(theta,1,n);
+elseif  lt ~= n
+  error(sprintf('Length of theta must be 1 or %d',n))
+else
+  theta = theta(:).';
+end
+
+td = abs(d) .* repmat(-theta, m, 1);
+r = exp(sum(td,2));
+
+if  nargout > 1
+  dr = repmat(-theta,m,1) .* sign(d) .* repmat(r,1,n);
+end
Index: issm/trunk/src/m/utils/DACE/correxpg.m
===================================================================
--- issm/trunk/src/m/utils/DACE/correxpg.m	(revision 1080)
+++ issm/trunk/src/m/utils/DACE/correxpg.m	(revision 1080)
@@ -0,0 +1,40 @@
+function  [r, dr] = correxpg(theta, d)
+%CORREXPG  General exponential correlation function
+%
+%           n
+%   r_i = prod exp(-theta_j * d_ij^theta_n+1)
+%          j=1
+%
+% If n > 1 and length(theta) = 2, then the model is isotropic: 
+% theta_j = theta(1), j=1,...,n;  theta_(n+1) = theta(2) 
+%
+% Call:    r = correxpg(theta, d)
+%          [r, dr] = correxpg(theta, d)
+%
+% theta :  parameters in the correlation function
+% d     :  m*n matrix with differences between given data points
+% r     :  correlation
+% dr    :  m*n matrix with the Jacobian of r at x. It is
+%          assumed that x is given implicitly by d(i,:) = x - S(i,:), 
+%          where S(i,:) is the i'th design site. 
+
+% hbn@imm.dtu.dk  
+% Last update April 12, 2002
+
+[m n] = size(d);  % number of differences and dimension of data
+lt = length(theta);
+if  n > 1 & lt == 2
+  theta = [repmat(theta(1),1,n)  theta(2)];
+elseif  lt ~= n+1
+  error(sprintf('Length of theta must be 2 or %d',n+1))
+else
+  theta = theta(:).';
+end
+
+pow = theta(end);   tt = repmat(-theta(1:n), m, 1);
+td = abs(d).^pow .* tt;
+r = exp(sum(td,2));
+
+if  nargout > 1
+  dr = pow  * tt .* sign(d) .* (abs(d) .^ (pow-1)) .* repmat(r,1,n);
+end
Index: issm/trunk/src/m/utils/DACE/corrgauss.m
===================================================================
--- issm/trunk/src/m/utils/DACE/corrgauss.m	(revision 1080)
+++ issm/trunk/src/m/utils/DACE/corrgauss.m	(revision 1080)
@@ -0,0 +1,36 @@
+function  [r, dr] = corrgauss(theta, d)
+%CORRGAUSS  Gaussian correlation function,
+%
+%           n
+%   r_i = prod exp(-theta_j * d_ij^2) ,  i = 1,...,m
+%          j=1
+%
+% If length(theta) = 1, then the model is isotropic:
+% all  theta_j = theta .
+%
+% Call:    r = corrgauss(theta, d)
+%          [r, dr] = corrgauss(theta, d)
+%
+% theta :  parameters in the correlation function
+% d     :  m*n matrix with differences between given data points
+% r     :  correlation
+% dr    :  m*n matrix with the Jacobian of r at x. It is
+%          assumed that x is given implicitly by d(i,:) = x - S(i,:), 
+%          where S(i,:) is the i'th design site. 
+
+% hbn@imm.dtu.dk  
+% Last update June 2, 2002
+
+[m n] = size(d);  % number of differences and dimension of data
+if  length(theta) == 1
+  theta = repmat(theta,1,n);
+elseif  length(theta) ~= n
+  error(sprintf('Length of theta must be 1 or %d',n))
+end
+
+td = d.^2 .* repmat(-theta(:).',m,1);
+r = exp(sum(td, 2));
+
+if  nargout > 1
+  dr = repmat(-2*theta(:).',m,1) .* d .* repmat(r,1,n);
+end
Index: issm/trunk/src/m/utils/DACE/corrlin.m
===================================================================
--- issm/trunk/src/m/utils/DACE/corrlin.m	(revision 1080)
+++ issm/trunk/src/m/utils/DACE/corrlin.m	(revision 1080)
@@ -0,0 +1,39 @@
+function  [r, dr] = corrlin(theta, d)
+%CORRLIN  Linear correlation function,
+%
+%           n
+%   r_i = prod max(0, 1 - theta_j * d_ij) ,  i = 1,...,m
+%          j=1
+%
+% If length(theta) = 1, then the model is isotropic:
+% all theta_j = theta .
+%
+% Call:    r = corrlin(theta, d)
+%          [r, dr] = corrlin(theta, d)
+%
+% theta :  parameters in the correlation function
+% d     :  m*n matrix with differences between given data points
+% r     :  correlation
+% dr    :  m*n matrix with the Jacobian of r at x. It is
+%          assumed that x is given implicitly by d(i,:) = x - S(i,:), 
+%          where S(i,:) is the i'th design site. 
+
+% hbn@imm.dtu.dk  
+% Last update April 12, 2002
+
+[m n] = size(d);  % number of differences and dimension of data
+if  length(theta) == 1
+  theta = repmat(theta,1,n);
+elseif  length(theta) ~= n
+  error(sprintf('Length of theta must be 1 or %d',n))
+end
+
+td = max(1 - abs(d) .* repmat(theta(:).',m,1), 0);
+r = prod(td, 2);
+
+if  nargout > 1
+  dr = zeros(m,n);
+  for  j = 1 : n
+    dr(:,j) = prod(td(:,[1:j-1 j+1:n]),2) .* (-theta(j) * sign(d(:,j)));
+  end
+end
Index: issm/trunk/src/m/utils/DACE/corrspherical.m
===================================================================
--- issm/trunk/src/m/utils/DACE/corrspherical.m	(revision 1080)
+++ issm/trunk/src/m/utils/DACE/corrspherical.m	(revision 1080)
@@ -0,0 +1,42 @@
+function  [r, dr] = corrspherical(theta, d)
+%CORRSPHERICAL  Spherical correlation function,
+%
+%           n
+%   r_i = prod max(0, 1 - 1.5(theta_j*d_ij) + .5(theta_j*d_ij)^3) ,  i = 1,...,m
+%          j=1
+%
+% If length(theta) = 1, then the model is isotropic:
+% all  theta_j = theta .
+%
+% Call:    r = corrspherical(theta, d)
+%          [r, dr] = corrspherical(theta, d)
+%
+% theta :  parameters in the correlation function
+% d     :  m*n matrix with differences between given data points
+% r     :  correlation
+% dr    :  m*n matrix with the Jacobian of r at x. It is
+%          assumed that x is given implicitly by d(i,:) = x - S(i,:), 
+%          where S(i,:) is the i'th design site. 
+
+% hbn@imm.dtu.dk  
+% Last update April 12, 2002
+
+[m n] = size(d);  % number of differences and dimension of data
+if  length(theta) == 1
+  theta = repmat(theta,1,n);
+elseif  length(theta) ~= n
+  error(sprintf('Length of theta must be 1 or %d',n))
+else
+  theta = theta(:).';
+end
+td = min(abs(d) .* repmat(theta,m,1), 1);
+ss = 1 - td .* (1.5 - .5*td.^2);
+r = prod(ss, 2);
+
+if  nargout > 1
+  dr = zeros(m,n);
+  for  j = 1 : n
+    dd = 1.5*theta(j) * sign(d(:,j)).*(td(:,j).^2 - 1);
+    dr(:,j) = prod(ss(:,[1:j-1 j+1:n]),2) .* dd;
+  end
+end  
Index: issm/trunk/src/m/utils/DACE/corrspline.m
===================================================================
--- issm/trunk/src/m/utils/DACE/corrspline.m	(revision 1080)
+++ issm/trunk/src/m/utils/DACE/corrspline.m	(revision 1080)
@@ -0,0 +1,67 @@
+function  [r, dr] = corrspline(theta, d)
+%CORRSPLINE  Cubic spline correlation function,
+%
+%           n
+%   r_i = prod S(theta_j*d_ij) ,  i = 1,...,m
+%          j=1
+%
+% with
+%           1 - 15x^2 + 30x^3   for   0 <= x <= 0.5
+%   S(x) =  1.25(1 - x)^3       for  0.5 < x < 1
+%           0                   for    x >= 1
+% If length(theta) = 1, then the model is isotropic:
+% all  theta_j = theta.
+%
+% Call:    r = corrspline(theta, d)
+%          [r, dr] = corrspline(theta, d)
+%
+% theta :  parameters in the correlation function
+% d     :  m*n matrix with differences between given data points
+% r     :  correlation
+% dr    :  m*n matrix with the Jacobian of r at x. It is
+%          assumed that x is given implicitly by d(i,:) = x - S(i,:), 
+%          where S(i,:) is the i'th design site. 
+
+% hbn@imm.dtu.dk  
+% Last update May 30, 2002
+
+[m n] = size(d);  % number of differences and dimension of data
+if  length(theta) == 1
+  theta = repmat(theta,1,n);
+elseif  length(theta) ~= n
+  error(sprintf('Length of theta must be 1 or %d',n))
+else
+  theta = theta(:).';
+end
+mn = m*n;   ss = zeros(mn,1);
+xi = reshape(abs(d) .* repmat(theta,m,1), mn,1);
+% Contributions to first and second part of spline
+i1 = find(xi <= 0.2);
+i2 = find(0.2 < xi & xi < 1);
+if  ~isempty(i1)
+  ss(i1) = 1 - xi(i1).^2 .* (15  - 30*xi(i1));
+end
+if  ~isempty(i2)
+  ss(i2) = 1.25 * (1 - xi(i2)).^3;
+end
+% Values of correlation
+ss = reshape(ss,m,n);
+r = prod(ss, 2);
+
+if  nargout > 1  % get Jacobian
+  u = reshape(sign(d) .* repmat(theta,m,1), mn,1);
+  dr = zeros(mn,1);
+  if  ~isempty(i1)
+    dr(i1) = u(i1) .* ( (90*xi(i1) - 30) .* xi(i1) );
+  end
+  if  ~isempty(i2)
+    dr(i2) = -3.75 * u(i2) .* (1 - xi(i2)).^2;
+  end
+  ii = 1 : m;
+  for  j = 1 : n
+    sj = ss(:,j);  ss(:,j) = dr(ii);
+    dr(ii) = prod(ss,2);
+    ss(:,j) = sj;   ii = ii + m;
+  end
+  dr = reshape(dr,m,n);
+end  % Jacobian
Index: issm/trunk/src/m/utils/DACE/dacefit.m
===================================================================
--- issm/trunk/src/m/utils/DACE/dacefit.m	(revision 1080)
+++ issm/trunk/src/m/utils/DACE/dacefit.m	(revision 1080)
@@ -0,0 +1,298 @@
+function  [dmodel, perf] = dacefit(S, Y, regr, corr, theta0, lob, upb)
+%DACEFIT Constrained non-linear least-squares fit of a given correlation
+% model to the provided data set and regression model
+%
+% Call
+%   [dmodel, perf] = dacefit(S, Y, regr, corr, theta0)
+%   [dmodel, perf] = dacefit(S, Y, regr, corr, theta0, lob, upb)
+%
+% Input
+% S, Y    : Data points (S(i,:), Y(i,:)), i = 1,...,m
+% regr    : Function handle to a regression model
+% corr    : Function handle to a correlation function
+% theta0  : Initial guess on theta, the correlation function parameters
+% lob,upb : If present, then lower and upper bounds on theta
+%           Otherwise, theta0 is used for theta
+%
+% Output
+% dmodel  : DACE model: a struct with the elements
+%    regr   : function handle to the regression model
+%    corr   : function handle to the correlation function
+%    theta  : correlation function parameters
+%    beta   : generalized least squares estimate
+%    gamma  : correlation factors
+%    sigma2 : maximum likelihood estimate of the process variance
+%    S      : scaled design sites
+%    Ssc    : scaling factors for design arguments
+%    Ysc    : scaling factors for design ordinates
+%    C      : Cholesky factor of correlation matrix
+%    Ft     : Decorrelated regression matrix
+%    G      : From QR factorization: Ft = Q*G' .
+% perf    : struct with performance information. Elements
+%    nv     : Number of evaluations of objective function
+%    perf   : (q+2)*nv array, where q is the number of elements 
+%             in theta, and the columns hold current values of
+%                 [theta;  psi(theta);  type]
+%             |type| = 1, 2 or 3, indicate 'start', 'explore' or 'move'
+%             A negative value for type indicates an uphill step
+
+% hbn@imm.dtu.dk  
+% Last update September 3, 2002
+
+% Check design points
+[m n] = size(S);  % number of design sites and their dimension
+sY = size(Y);
+if  min(sY) == 1,  Y = Y(:);   lY = max(sY);  sY = size(Y);
+else,              lY = sY(1); end
+if m ~= lY
+  error('S and Y must have the same number of rows'), end
+
+% Check correlation parameters
+lth = length(theta0);
+if  nargin > 5  % optimization case
+  if  length(lob) ~= lth | length(upb) ~= lth
+    error('theta0, lob and upb must have the same length'), end
+  if  any(lob <= 0) | any(upb < lob)
+    error('The bounds must satisfy  0 < lob <= upb'), end
+else  % given theta
+  if  any(theta0 <= 0)
+    error('theta0 must be strictly positive'), end
+end
+
+% Normalize data
+mS = mean(S);   sS = std(S);
+mY = mean(Y);   sY = std(Y);
+% 02.08.27: Check for 'missing dimension'
+j = find(sS == 0);
+if  ~isempty(j),  sS(j) = 1; end
+j = find(sY == 0);
+if  ~isempty(j),  sY(j) = 1; end
+S = (S - repmat(mS,m,1)) ./ repmat(sS,m,1);
+Y = (Y - repmat(mY,m,1)) ./ repmat(sY,m,1);
+
+% Calculate distances D between points
+mzmax = m*(m-1) / 2;        % number of non-zero distances
+ij = zeros(mzmax, 2);       % initialize matrix with indices
+D = zeros(mzmax, n);        % initialize matrix with distances
+ll = 0;
+for k = 1 : m-1
+  ll = ll(end) + (1 : m-k);
+  ij(ll,:) = [repmat(k, m-k, 1) (k+1 : m)']; % indices for sparse matrix
+  D(ll,:) = repmat(S(k,:), m-k, 1) - S(k+1:m,:); % differences between points
+end
+if  min(sum(abs(D),2) ) == 0
+  error('Multiple design sites are not allowed'), end
+
+% Regression matrix
+F = feval(regr, S);  [mF p] = size(F);
+if  mF ~= m, error('number of rows in  F  and  S  do not match'), end
+if  p > mF,  error('least squares problem is underdetermined'), end
+
+% parameters for objective function
+par = struct('corr',corr, 'regr',regr, 'y',Y, 'F',F, ...
+  'D', D, 'ij',ij, 'scS',sS);
+
+% Determine theta
+if  nargin > 5
+  % Bound constrained non-linear optimization
+  [theta f fit perf] = boxmin(theta0, lob, upb, par);
+  if  isinf(f)
+    error('Bad parameter region.  Try increasing  upb'), end
+else
+  % Given theta
+  theta = theta0(:);   
+  [f  fit] = objfunc(theta, par);
+  perf = struct('perf',[theta; f; 1], 'nv',1);
+  if  isinf(f)
+    error('Bad point.  Try increasing theta0'), end
+end
+
+% Return values
+dmodel = struct('regr',regr, 'corr',corr, 'theta',theta.', ...
+  'beta',fit.beta, 'gamma',fit.gamma, 'sigma2',sY.^2.*fit.sigma2, ...
+  'S',S, 'Ssc',[mS; sS], 'Ysc',[mY; sY], ...
+  'C',fit.C, 'Ft',fit.Ft, 'G',fit.G);
+
+% >>>>>>>>>>>>>>>>   Auxiliary functions  ====================
+
+function  [obj, fit] = objfunc(theta, par)
+% Initialize
+obj = inf; 
+fit = struct('sigma2',NaN, 'beta',NaN, 'gamma',NaN, ...
+    'C',NaN, 'Ft',NaN, 'G',NaN);
+m = size(par.F,1);
+% Set up  R
+r = feval(par.corr, theta, par.D);
+idx = find(r > 0);   o = (1 : m)';   
+mu = (10+m)*eps;
+R = sparse([par.ij(idx,1); o], [par.ij(idx,2); o], ...
+  [r(idx); ones(m,1)+mu]);  
+% Cholesky factorization with check for pos. def.
+[C rd] = chol(R);
+if  rd,  return, end % not positive definite
+
+% Get least squares solution
+C = C';   Ft = C \ par.F;
+[Q G] = qr(Ft,0);
+if  rcond(G) < 1e-10
+  % Check   F  
+  if  cond(par.F) > 1e15 
+    T = sprintf('F is too ill conditioned\nPoor combination of regression model and design sites');
+    error(T)
+  else  % Matrix  Ft  is too ill conditioned
+    return 
+  end 
+end
+Yt = C \ par.y;   beta = G \ (Q'*Yt);
+rho = Yt - Ft*beta;  sigma2 = sum(rho.^2)/m;
+detR = prod( full(diag(C)) .^ (2/m) );
+obj = sum(sigma2) * detR;
+if  nargout > 1
+  fit = struct('sigma2',sigma2, 'beta',beta, 'gamma',rho' / C, ...
+    'C',C, 'Ft',Ft, 'G',G');
+end
+
+% --------------------------------------------------------
+
+function  [t, f, fit, perf] = boxmin(t0, lo, up, par)
+%BOXMIN  Minimize with positive box constraints
+
+% Initialize
+[t, f, fit, itpar] = start(t0, lo, up, par);
+if  ~isinf(f)
+  % Iterate
+  p = length(t);
+  if  p <= 2,  kmax = 2; else,  kmax = min(p,4); end
+  for  k = 1 : kmax
+    th = t;
+    [t, f, fit, itpar] = explore(t, f, fit, itpar, par);
+    [t, f, fit, itpar] = move(th, t, f, fit, itpar, par);
+  end
+end
+perf = struct('nv',itpar.nv, 'perf',itpar.perf(:,1:itpar.nv));
+
+% --------------------------------------------------------
+
+function  [t, f, fit, itpar] = start(t0, lo, up, par)
+% Get starting point and iteration parameters
+
+% Initialize
+t = t0(:);  lo = lo(:);   up = up(:);   p = length(t);
+D = 2 .^ ([1:p]'/(p+2));
+ee = find(up == lo);  % Equality constraints
+if  ~isempty(ee)
+  D(ee) = ones(length(ee),1);   t(ee) = up(ee); 
+end
+ng = find(t < lo | up < t);  % Free starting values
+if  ~isempty(ng)
+  t(ng) = (lo(ng) .* up(ng).^7).^(1/8);  % Starting point
+end
+ne = find(D ~= 1);
+
+% Check starting point and initialize performance info
+[f  fit] = objfunc(t,par);   nv = 1;
+itpar = struct('D',D, 'ne',ne, 'lo',lo, 'up',up, ...
+  'perf',zeros(p+2,200*p), 'nv',1);
+itpar.perf(:,1) = [t; f; 1];
+if  isinf(f)    % Bad parameter region
+  return
+end
+
+if  length(ng) > 1  % Try to improve starting guess
+  d0 = 16;  d1 = 2;   q = length(ng);
+  th = t;   fh = f;   jdom = ng(1);  
+  for  k = 1 : q
+    j = ng(k);    fk = fh;  tk = th;
+    DD = ones(p,1);  DD(ng) = repmat(1/d1,q,1);  DD(j) = 1/d0;
+    alpha = min(log(lo(ng) ./ th(ng)) ./ log(DD(ng))) / 5;
+    v = DD .^ alpha;   tk = th;
+    for  rept = 1 : 4
+      tt = tk .* v; 
+      [ff  fitt] = objfunc(tt,par);  nv = nv+1;
+      itpar.perf(:,nv) = [tt; ff; 1];
+      if  ff <= fk 
+        tk = tt;  fk = ff;
+        if  ff <= f
+          t = tt;  f = ff;  fit = fitt; jdom = j;
+        end
+      else
+        itpar.perf(end,nv) = -1;   break
+      end
+    end
+  end % improve
+  
+  % Update Delta  
+  if  jdom > 1
+    D([1 jdom]) = D([jdom 1]); 
+    itpar.D = D;
+  end
+end % free variables
+
+itpar.nv = nv;
+
+% --------------------------------------------------------
+
+function  [t, f, fit, itpar] = explore(t, f, fit, itpar, par)
+% Explore step
+
+nv = itpar.nv;   ne = itpar.ne;
+for  k = 1 : length(ne)
+  j = ne(k);   tt = t;   DD = itpar.D(j);
+  if  t(j) == itpar.up(j)
+    atbd = 1;   tt(j) = t(j) / sqrt(DD);
+  elseif  t(j) == itpar.lo(j)
+    atbd = 1;  tt(j) = t(j) * sqrt(DD);
+  else
+    atbd = 0;  tt(j) = min(itpar.up(j), t(j)*DD);
+  end
+  [ff  fitt] = objfunc(tt,par);  nv = nv+1;
+  itpar.perf(:,nv) = [tt; ff; 2];
+  if  ff < f
+    t = tt;  f = ff;  fit = fitt;
+  else
+    itpar.perf(end,nv) = -2;
+    if  ~atbd  % try decrease
+      tt(j) = max(itpar.lo(j), t(j)/DD);
+      [ff  fitt] = objfunc(tt,par);  nv = nv+1;
+      itpar.perf(:,nv) = [tt; ff; 2];
+      if  ff < f
+        t = tt;  f = ff;  fit = fitt;
+      else
+        itpar.perf(end,nv) = -2;
+      end
+    end
+  end
+end % k
+
+itpar.nv = nv;
+
+% --------------------------------------------------------
+
+function  [t, f, fit, itpar] = move(th, t, f, fit, itpar, par)
+% Pattern move
+
+nv = itpar.nv;   ne = itpar.ne;   p = length(t);
+v = t ./ th;
+if  all(v == 1)
+  itpar.D = itpar.D([2:p 1]).^.2;
+  return
+end
+
+% Proper move
+rept = 1;
+while  rept
+  tt = min(itpar.up, max(itpar.lo, t .* v));  
+  [ff  fitt] = objfunc(tt,par);  nv = nv+1;
+  itpar.perf(:,nv) = [tt; ff; 3];
+  if  ff < f
+    t = tt;  f = ff;  fit = fitt;
+    v = v .^ 2;
+  else
+    itpar.perf(end,nv) = -3;
+    rept = 0;
+  end
+  if  any(tt == itpar.lo | tt == itpar.up), rept = 0; end
+end
+
+itpar.nv = nv;
+itpar.D = itpar.D([2:p 1]).^.25;
Index: issm/trunk/src/m/utils/DACE/dsmerge.m
===================================================================
--- issm/trunk/src/m/utils/DACE/dsmerge.m	(revision 1080)
+++ issm/trunk/src/m/utils/DACE/dsmerge.m	(revision 1080)
@@ -0,0 +1,141 @@
+function  [mS, mY] = dsmerge(S, Y, ds, nms, wtds, wtdy)
+%DSMERGE  Merge data for multiple design sites.
+%
+% Call
+%   [mS, mY] = dsmerge(S, Y)
+%   [mS, mY] = dsmerge(S, Y, ds)
+%   [mS, mY] = dsmerge(S, Y, ds, nms)
+%   [mS, mY] = dsmerge(S, Y, ds, nms, wtds)
+%   [mS, mY] = dsmerge(S, Y, ds, nms, wtds, wtdy)
+%
+% Input
+% S, Y : Data points (S(i,:), Y(i,:)), i = 1,...,m
+% ds   : Threshold for equal, normalized sites. Default is 1e-14.
+% nms  : Norm, in which the distance is measured.
+%        nms =  1 : 1-norm (sum of absolute coordinate differences)
+%               2 : 2-norm (Euclidean distance) (default)
+%        otherwise: infinity norm (max coordinate difference)      
+% wtds : What to do with the S-values in case of multiple points.
+%        wtds = 1 : return the mean value (default)
+%               2 : return the median value
+%               3 : return the 'cluster center'
+% wtdy : What to do with the Y-values in case of multiple points.
+%        wtdy = 1 : return the mean value (default)
+%               2 : return the median value
+%               3 : return the 'cluster center' value
+%               4 : return the minimum value
+%               5 : return the maximum value    
+%
+% Output
+% mS : Compressed design sites, with multiple points merged
+%      according to wtds
+% mY : Responses, compressed according to wtdy
+
+% hbn@imm.dtu.dk  
+% Last update July 3, 2002
+
+% Check design points
+[m n] = size(S);  % number of design sites and their dimension
+sY = size(Y);
+if  min(sY) == 1,  Y = Y(:);   lY = max(sY);  sY = size(Y);
+else,              lY = sY(1); end
+if m ~= lY
+  error('S and Y must have the same number of rows'), end
+
+% Threshold
+if  nargin < 3
+  ds = 1e-14;
+elseif  (ds < 0) | (ds > .5)
+  error('ds must be in the range [0, 0.5]'), end
+
+% Which measure
+if  nargin < 4
+  nms = 2;
+elseif  (nms ~= 1) & (nms ~= 2)
+  nms = Inf;
+end
+
+% What to do
+if  nargin < 5
+  wtds = 1;
+else
+  wtds = round(wtds);
+  if  (wtds < 1) | (wtds > 3)
+    error('wtds must be in the range [1, 3]'), end
+end
+if  nargin < 6
+  wtdy = 1;
+else
+  wtdy = round(wtdy);
+  if  (wtdy < 1) | (wtdy > 5)
+    error('wtdy must be in the range [1, 5]'), end
+end
+
+% Process data
+more = 1;
+ladr = zeros(1,ceil(m/2));
+while more
+  m = size(S,1);
+  D = zeros(m,m);
+  
+  % Normalize sites
+  mS = mean(S);   sS = std(S);
+  scS = (S - repmat(mS,m,1)) ./ repmat(sS,m,1);
+  
+  % Calculate distances D (upper triangle of the symetric matrix)
+  for k = 1 : m-1
+    kk = k+1 : m;
+    dk = abs(repmat(scS(k,:), m-k, 1) - scS(kk,:));
+    if  nms == 1,      D(kk,k) = sum(dk,2);
+    elseif  nms == 2,  D(kk,k) = sqrt(sum(dk.^2,2));
+    else,              D(kk,k) = max(dk,[],2);      end
+  end
+  D = D + D'; % make D symetric
+  
+  % Check distances
+  mult = zeros(1,m);
+  for  j = 1 : m
+    % Find the number of multiple sites in each column of D
+    mult(j) = length(find(D(:,j) < ds));
+  end
+  % Find the first column with the maximum number of multiple sites
+  [mmult jj] = max(mult);
+  
+  if  mmult == 1,  more = 0;
+  else
+    nm = 0;
+    while  mmult > 1
+      nm = nm + 1;  % no. of points to merge
+      ladr(nm) = jj;
+      
+      % Merge point no jj and its neighbours, note that jj is the center
+      % of the cluster, as it has the most neighbors (among the multiple sites)
+      ngb = find(D(:,jj) < ds);
+      
+      switch  wtds
+      case 1,  S(jj,:) = mean(S(ngb,:));
+      case 2,  S(jj,:) = median(S(ngb,:));
+      case 3,  S(jj,:) = S(jj,:);
+      end
+      
+      switch  wtdy
+      case 1,  Y(jj,:) = mean(Y(ngb,:));
+      case 2,  Y(jj,:) = median(Y(ngb,:));
+      case 3,  Y(jj,:) = Y(jj,:);
+      case 4,  Y(jj,:) = min(Y(ngb,:));
+      case 5,  Y(jj,:) = max(Y(ngb,:));
+      end
+      
+      % Delete from list
+      mult(ngb) = 0;
+      [mmult jj] = max(mult);
+    end
+
+    % Reduced data set
+    act = [find(mult > 0)  ladr(1:nm)];
+    S = S(act,:);    Y = Y(act,:);
+  end % multiple
+end % loop
+
+% Return reduced set
+mS = S;   mY = Y;
Index: issm/trunk/src/m/utils/DACE/gridsamp.m
===================================================================
--- issm/trunk/src/m/utils/DACE/gridsamp.m	(revision 1080)
+++ issm/trunk/src/m/utils/DACE/gridsamp.m	(revision 1080)
@@ -0,0 +1,45 @@
+function  S = gridsamp(range, q)
+%GRIDSAMP  n-dimensional grid over given range
+%
+% Call:    S = gridsamp(range, q)
+%
+% range :  2*n matrix with lower and upper limits
+% q     :  n-vector, q(j) is the number of points
+%          in the j'th direction.
+%          If q is a scalar, then all q(j) = q
+% S     :  m*n array with points, m = prod(q)
+
+% hbn@imm.dtu.dk  
+% Last update June 25, 2002
+
+[mr n] = size(range);    dr = diff(range);
+if  mr ~= 2 | any(dr < 0)
+  error('range must be an array with two rows and range(1,:) <= range(2,:)')
+end 
+sq = size(q);
+if  min(sq) > 1 | any(q <= 0)
+  error('q must be a vector with non-negative elements')
+end
+p = length(q);   
+if  p == 1,  q = repmat(q,1,n); 
+elseif  p ~= n
+  error(sprintf('length of q must be either 1 or %d',n))
+end 
+
+% Check for degenerate intervals
+i = find(dr == 0);
+if  ~isempty(i),  q(i) = 0*q(i); end
+
+% Recursive computation
+if  n > 1
+  A = gridsamp(range(:,2:end), q(2:end));  % Recursive call
+  [m p] = size(A);   q = q(1);
+  S = [zeros(m*q,1) repmat(A,q,1)];
+  y = linspace(range(1,1),range(2,1), q);
+  k = 1:m;
+  for  i = 1 : q
+    S(k,1) = repmat(y(i),m,1);  k = k + m;
+  end
+else    
+  S = linspace(range(1,1),range(2,1), q).';
+end
Index: issm/trunk/src/m/utils/DACE/lhsamp.m
===================================================================
--- issm/trunk/src/m/utils/DACE/lhsamp.m	(revision 1080)
+++ issm/trunk/src/m/utils/DACE/lhsamp.m	(revision 1080)
@@ -0,0 +1,23 @@
+function S = lhsamp(m, n)
+%LHSAMP  Latin hypercube distributed random numbers
+%
+% Call:    S = lhsamp
+%          S = lhsamp(m)
+%          S = lhsamp(m, n)
+%
+% m : number of sample points to generate, if unspecified m = 1
+% n : number of dimensions, if unspecified n = m
+%
+% S : the generated n dimensional m sample points chosen from
+%     uniform distributions on m subdivions of the interval (0.0, 1.0)
+
+% hbn@imm.dtu.dk  
+% Last update April 12, 2002
+
+if nargin < 1, m = 1; end
+if nargin < 2, n = m; end
+
+S = zeros(m,n);
+for i = 1 : n
+  S(:, i) = (rand(1, m) + (randperm(m) - 1))' / m;
+end
Index: issm/trunk/src/m/utils/DACE/predictor.m
===================================================================
--- issm/trunk/src/m/utils/DACE/predictor.m	(revision 1080)
+++ issm/trunk/src/m/utils/DACE/predictor.m	(revision 1080)
@@ -0,0 +1,132 @@
+function  [y, or1, or2, dmse] = predictor(x, dmodel)
+%PREDICTOR  Predictor for y(x) using the given DACE model.
+%
+% Call:   y = predictor(x, dmodel)
+%         [y, or] = predictor(x, dmodel)
+%         [y, dy, mse] = predictor(x, dmodel) 
+%         [y, dy, mse, dmse] = predictor(x, dmodel) 
+%
+% Input
+% x      : trial design sites with n dimensions.  
+%          For mx trial sites x:
+%          If mx = 1, then both a row and a column vector is accepted,
+%          otherwise, x must be an mx*n matrix with the sites stored
+%          rowwise.
+% dmodel : Struct with DACE model; see DACEFIT
+%
+% Output
+% y    : predicted response at x.
+% or   : If mx = 1, then or = gradient vector/Jacobian matrix of predictor
+%        otherwise, or is an vector with mx rows containing the estimated
+%                   mean squared error of the predictor
+% Three or four results are allowed only when mx = 1,
+% dy   : Gradient of predictor; column vector with  n elements
+% mse  : Estimated mean squared error of the predictor;
+% dmse : Gradient vector/Jacobian matrix of mse
+
+% hbn@imm.dtu.dk
+% Last update August 26, 2002
+ 
+  or1 = NaN;   or2 = NaN;  dmse = NaN;  % Default return values
+  if  isnan(dmodel.beta)
+    y = NaN;   
+    error('DMODEL has not been found')
+  end
+
+  [m n] = size(dmodel.S);  % number of design sites and number of dimensions
+  sx = size(x);            % number of trial sites and their dimension
+  if  min(sx) == 1 & n > 1 % Single trial point 
+    nx = max(sx);
+    if  nx == n 
+      mx = 1;  x = x(:).';
+    end
+  else
+    mx = sx(1);  nx = sx(2);
+  end
+  if  nx ~= n
+    error(sprintf('Dimension of trial sites should be %d',n))
+  end
+  
+  % Normalize trial sites  
+  x = (x - repmat(dmodel.Ssc(1,:),mx,1)) ./ repmat(dmodel.Ssc(2,:),mx,1);
+  q = size(dmodel.Ysc,2);  % number of response functions
+  y = zeros(mx,q);         % initialize result
+  
+  if  mx == 1  % one site only
+    dx = repmat(x,m,1) - dmodel.S;  % distances to design sites
+    if  nargout > 1                 % gradient/Jacobian wanted
+      [f df] = feval(dmodel.regr, x);
+      [r dr] = feval(dmodel.corr, dmodel.theta, dx);
+      % Scaled Jacobian
+      dy = (df * dmodel.beta).' + dmodel.gamma * dr;
+      % Unscaled Jacobian
+      or1 = dy .* repmat(dmodel.Ysc(2, :)', 1, nx) ./ repmat(dmodel.Ssc(2,:), q, 1);
+      if q == 1
+        % Gradient as a column vector
+        or1 = or1';
+      end
+      if  nargout > 2  % MSE wanted
+        
+        rt = dmodel.C \ r;
+        u = dmodel.Ft.' * rt - f.';
+        v = dmodel.G \ u;
+        or2 = repmat(dmodel.sigma2,mx,1) .* repmat((1 + sum(v.^2) - sum(rt.^2))',1,q);
+        
+        if  nargout > 3  % gradient/Jacobian of MSE wanted
+          % Scaled gradient as a row vector
+          Gv = dmodel.G' \ v;
+          g = (dmodel.Ft * Gv - rt)' * (dmodel.C \ dr) - (df * Gv)';
+          % Unscaled Jacobian
+          dmse = repmat(2 * dmodel.sigma2',1,nx) .* repmat(g ./ dmodel.Ssc(2,:),q,1);
+          if q == 1
+            % Gradient as a column vector
+            dmse = dmse';
+          end
+        end
+        
+      end
+      
+    else  % predictor only
+      f = feval(dmodel.regr, x);
+      r = feval(dmodel.corr, dmodel.theta, dx);
+    end
+    
+    % Scaled predictor
+    sy = f * dmodel.beta + (dmodel.gamma*r).';
+    % Predictor
+    y = (dmodel.Ysc(1,:) + dmodel.Ysc(2,:) .* sy)';
+    
+  else  % several trial sites
+    % Get distances to design sites  
+    dx = zeros(mx*m,n);  kk = 1:m;
+    for  k = 1 : mx
+      dx(kk,:) = repmat(x(k,:),m,1) - dmodel.S;
+      kk = kk + m;
+    end
+    % Get regression function and correlation
+    f = feval(dmodel.regr, x);
+    r = feval(dmodel.corr, dmodel.theta, dx);
+    r = reshape(r, m, mx);
+    
+    % Scaled predictor 
+    sy = f * dmodel.beta + (dmodel.gamma * r).';
+    % Predictor
+    y = repmat(dmodel.Ysc(1,:),mx,1) + repmat(dmodel.Ysc(2,:),mx,1) .* sy;
+    
+    if  nargout > 1   % MSE wanted
+      rt = dmodel.C \ r;
+      u = dmodel.G \ (dmodel.Ft.' * rt - f.');
+      or1 = repmat(dmodel.sigma2,mx,1) .* repmat((1 + colsum(u.^2) - colsum(rt.^2))',1,q);
+      if  nargout > 2
+        disp('WARNING from PREDICTOR.  Only  y  and  or1=mse  are computed')
+      end
+    end
+    
+  end % of several sites
+  
+% >>>>>>>>>>>>>>>>   Auxiliary function  ====================
+
+function  s = colsum(x)
+% Columnwise sum of elements in  x
+if  size(x,1) == 1,  s = x; 
+else,                s = sum(x);  end
Index: issm/trunk/src/m/utils/DACE/regpoly0.m
===================================================================
--- issm/trunk/src/m/utils/DACE/regpoly0.m	(revision 1080)
+++ issm/trunk/src/m/utils/DACE/regpoly0.m	(revision 1080)
@@ -0,0 +1,18 @@
+function  [f, df] = regpoly0(S)
+%REGPOLY0  Zero order polynomial regression function
+%
+% Call:    f = regpoly0(S)
+%          [f, df] = regpoly0(S)
+%
+% S  : m*n matrix with design sites
+% f  : ones(m,1)
+% df : Jacobian at the first point (first row in S) 
+
+% hbn@imm.dtu.dk  
+% Last update  April 12, 2002
+
+[m n] = size(S);
+f = ones(m,1);
+if  nargout > 1
+  df = zeros(n,1);
+end
Index: issm/trunk/src/m/utils/DACE/regpoly1.m
===================================================================
--- issm/trunk/src/m/utils/DACE/regpoly1.m	(revision 1080)
+++ issm/trunk/src/m/utils/DACE/regpoly1.m	(revision 1080)
@@ -0,0 +1,18 @@
+function  [f, df] = regpoly1(S)
+%REGPOLY1  First order polynomial regression function
+%
+% Call:    f = regpoly1(S)
+%          [f, df] = regpoly1(S)
+%
+% S : m*n matrix with design sites
+% f = [1  s]
+% df : Jacobian at the first point (first row in S) 
+
+% hbn@imm.dtu.dk  
+% Last update April 12, 2002
+
+[m n] = size(S);
+f = [ones(m,1)  S];
+if  nargout > 1
+  df = [zeros(n,1) eye(n)];
+end
Index: issm/trunk/src/m/utils/DACE/regpoly2.m
===================================================================
--- issm/trunk/src/m/utils/DACE/regpoly2.m	(revision 1080)
+++ issm/trunk/src/m/utils/DACE/regpoly2.m	(revision 1080)
@@ -0,0 +1,31 @@
+function  [f, df] = regpoly2(S)
+%REGPOLY2  Second order polynomial regression function
+% Call:    f = regpoly2(S)
+%          [f, df] = regpoly2(S)
+%
+% S : m*n matrix with design sites
+% f =  [1 S S(:,1)*S S(:,2)S(:,2:n) ... S(:,n)^2]
+% df : Jacobian at the first point (first row in S) 
+
+% hbn@imm.dtu.dk  
+% Last update September 4, 2002
+
+[m n] = size(S);
+nn = (n+1)*(n+2)/2;  % Number of columns in f  
+% Compute  f
+f = [ones(m,1) S zeros(m,nn-n-1)];
+j = n+1;   q = n;
+for  k = 1 : n
+  f(:,j+(1:q)) = repmat(S(:,k),1,q) .* S(:,k:n);
+  j = j+q;   q = q-1;
+end
+
+if  nargout > 1
+  df = [zeros(n,1)  eye(n)  zeros(n,nn-n-1)];
+  j = n+1;   q = n; 
+  for  k = 1 : n
+    df(k,j+(1:q)) = [2*S(1,k) S(1,k+1:n)];
+    for i = 1 : n-k,  df(k+i,j+1+i) = S(1,k); end
+    j = j+q;   q = q-1;
+  end
+end 
