function [minim,yhat]=quadmin(x,y,yprime,A)
% fits a quadratic response surface to the y using the columns of x.
% then determines the minimum of this surface. yprime is estimate of derivatives at same points.
% if x is n by k, then yprime is n by k and column j is estimated gradient in direction j.
% the input A is such that A'*A=sigma^-1  and sigma is 4n by 4n
% A=msqrt(inv(sigma));
design = quaddesign3(x)
%sigma=cov([y yprime]);   % the covariance matrix of each observation
alpha=.05;
% Fit response surface model design
if size(design,2) > size(design,1)
    s = ['Your data is insufficient to fit a  quadratic model.'];
    error(s);
end
[Q, R] = qr(A*design,0);
if rcond(R) < 1E-12
    s = ['Your data is insufficient to fit a quadratic model.'];
    error(s);
end
y=reshape([y yprime]',prod(4*size(x,1)),1);
%beta=inv(design'*A'*A*design)*design'*A'*A*y    % these two lines are the same. 
beta = R\(Q'*A*y)
%[y design*[0 0 0 0 .5 0 0 1 1 1]']    % output to debug the program...remove later
H=[beta([8 5 6]) beta([5 9 7]) beta([6 7 10])];
b=beta([2 3 4]);
%for i=1:size(x,1)
 %  x(i,:)*H*(x(i,:))'+x(i,:)*b+beta(1)   % this is the value of the estimated function at x(i,:)
  % b+2*H*(x(i,:)')            %  gradient
%end

minim=-.5*inv(H)*b;

yhat = design*beta;
residuals = y - yhat

% p=length(beta)
%df = max(size(y)) - p;

%rmse = sqrt(sum(residuals.*residuals)./df);
%crit = sqrt(p*finv(1 - alpha,p,df))*rmse;

