function [xcm,ycm,xucm,yucm]=optrsm(b0,b,B,opt,lohi)
%OPTRSM.M
%[xcm,ycm,xucm,yucm]=optrsm(b0,b,B,'opt',lohi)
% OR  optrsm(b0,b,B,'opt',lohi)
%     OPTIMISE	RSM MODELS == USING CONSTR.M ALGORTHM
% Uses MAXRSM1.M and CONSTR.M
% Calculates CONSTRAINED/UNCONSTRAINED OPTIMUM for factors in RSM models
%
% IN: lohi  = n*2 matrix of lo and hi contraints for each factor x
%     b0,b,B =Parameter matrices from RSM Model out of LEASRSM(1).M
%    'opt ' = 'max' or 'min' depends   on problem
%
%
% OUT:	    Constrained/unconstrained optimum values of Factors --> xcm/xucm
%	    Optimum value of y(ycm,yucm) --> y= b0 +b.x +x'.B.x  (RSM Model)
%	    Canonical Model Info
[nlohi,junk]=size(lohi);
if length(b) ~= nlohi , error('Matrix lohi must have  #rows=#variables'),end
% B matrix from Leasrsm is lower diag: convert to full diag= BF
c_=B+B';
d_=diag(diag(c_));
BF=c_-d_/2;   % Full B Matrix
% set up parameters for maxrsm1
p1=b0;p2=b;p3=BF; p4=opt;
% set initial guess at mean( lo,hi)
x0_=mean(lohi')'; %gives mean values of each row of lohi
vlb=lohi(:,1);vub=lohi(:,2);	  %set lower/upper bounds
% for no equation constraints g=[], dont bother cal g in maxrsm1,set to []
options(1)=1; % display on
% parameter NAMES in CONSTR must be P1,P2 .. VALUES p1,p2 ..can be anything
xcm=constr('f=maxrsm1(x,P1,P2,P3,P4);g=[];',x0_,options,vlb,vub,[],p1,p2,p3,p4); % calc contrained opt
 ycm=maxrsm1(xcm,p1,p2,p3,p4); % calc corresponding constrained y
if opt=='max',ycm=-ycm;end % change sign back (see maxrsm)
if opt=='max',disp('Constrained Max'), else,disp('Constrained  Min'),end;
xcm
ycm
disp (' Unconstrained optimum --max or min')
xucm= -0.5*(BF\b') % watch out for \ precedence need ()
yucm=b0+0.5*xucm'*b'	% Ref Box&Draper p.337
%
% CANONICAL MODEL INFO
%
disp('Canonical Info')
disp(' Distance of unconstr/constr optimum from Design Center')
DUC=sqrt(xucm'*xucm) ,  DC=sqrt(xcm'*xcm)
disp('Eigen Vectors values')
[M,lam]=eig(BF)
disp(' Slopes A Canonical Form')
theta=M'*b'
disp('Unit Change in Response along each canonical direction')
r=sqrt(3*theta.^2 +diag(lam).^2)