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