function [Pf,g] = HOSRSM(limit_state,n,g_idx,tol,h_ord,h_reg,N,Zmcs)
% [Pf,g] = HOSRSM(limit_state,n,g_idx,tol,h_ord,h_reg,N,Zmcs)
%
% High Order Stochastic Response Surface Method (HOSRSM)
% 
% This program approximates the probability of failure of a system using
% the high order stochastic response surface method (HOSRSM). The limit
% state function is returned by the matlab function called limit_state.m
% as follows:  
%     function G = limit_state(X)
% where X is the vector X of standardized normal random variables.   
% Transformations from the standard normal form to the required 
% distribution for each random variable is carried out in the 
% limit_state() matlab function.
%
% G need not be a scalar.   If G is not a scalar, then the desired element
% of G is specified by the scalar integer input variable "g_idx".
%
% The HOSRSM approximates the true limit state function with a 
% polynomial reponse surface of arbitrary order,
%    g(X) = a + \sum_{i=1}^n \sum_{j=1}^{k_i) b_ij X_i^j + 
%           \sum_{q=1}^m c_q \prod_{i=1}^n X_i^{p_iq}.
% The first stage of the algorithm determines the correct polynomial order,
% k_i in the response surace. Then the formulation of the mixed terms 
% \sum_{q=1}^m c_q \prod_{i=1}^n X_i^{p_iq} are derived by the second stage
% based on previous results. In the third stage, the response surface 
% is approximated. In the last stage, a Monte Carlo simulation is carried
% out to determine the probability of failure.
%
% INPUT
% ------------
% limit_state  the limit state function in the form
%                function g = limit_state(z)
%              where z are sample values of the standardized normal
%              random variables involved in the limit state function
% n            number of random variables in the limit state function
% g_idx        index of the desired element of the limit state function
% tol          accuracy tolerance for response surface model  default = 0.01
% h_ord        parameter controlling the sampling domain in the polynomial
%              order determination (1st stage)        default = 3
% h_reg        parameter controlling the sampling domain in the response
%              surface approximation (3rd stage)      default = 3
% N            sample size for full scale Monte Carlo simulation (4th stage)
%                                                     default = 1e6
% Zmcs         optional matrix of Monte-Carlo samples of standardized r.v's
%              i.e., Zmcs = randn(N,n);
%
% OUTPUT
% ------------
% Pf           probability of failure
%  g           values of the limit state function evaluated using
%              Monte Carlo simulatiaon of the response surface
%
% Reference:
%    Gavin, HP and Yau SC, ``High order limit state functions in the 
%    response surface method for structural reliability analysis,''
%    submitted to Structural Safety, December 2005.     
%
% Department of Civil and Environmental Engineering
% Duke University
% Siu Chung Yau, Henri P. Gavin, January 2006, Updated October 2009

if nargin < 2
   help HOSRSM
   return
end

if nargin < 3, g_idx = 1;     else g_idx = floor(g_idx(1,1)); end
if nargin < 4, tol   = 0.01   else tol   = abs(tol);          end
if nargin < 5, h_ord = 3;     else h_ord = abs(h_ord);        end
if nargin < 6, h_reg = 3;     else h_reg = abs(h_reg);        end
if nargin < 7, N     = 1e6;   else N     = abs(N);            end

%warning off all;

% 0th Stage: Re-Seed Random Number Generator -------------------------------

time_now = clock;
randn('state',1000*time_now(6));  % reseed normal random number generator
                                     % with the current second
rand('state',1000*time_now(6));   % reseed uniform random number generator
                                     % with the current second
clear time_now

% 1st Stage: Polynomial Order Determination -------------------------------
%            By fitting the limit state function one-by-one with the 
%            one-dimensional Chebyshev polynomials,
%            d_0*T_0(X_i) + d_1*T_1(X_i) + ... + d_ki*T_ki(X_i),
%            along dimension X_i, the polynomial orders k_i are determined 
%            by statistically and numerically testing the significane of 
%            polynomial coefficients.

fprintf('\nHigh Order Stochastic Response Surface Method (HOSRSM)\n\n');
fprintf('1st Stage: Polynomial Order Determination ...\n');

order       =  1*ones(1,n);	% initial guess of response surface orders
order_time  =  zeros(1,n);
quality     =  zeros(1,n);
R2_order    = zeros(1,n);

no_fevals = 0;

if 1              % 1: do order determination; 0: skip order determination

  for i = 1:n     % determine the orders for each variable one by one 
    
    tic             % start recording the time required to determine the order
                
    test_order = 1;    % test the significance of the next higher order 
                       % (3rd order)
                       %
                       % 'test_next_ord' serves as an indicator to tell
                       % whether the significance test of the next higher
                       % order term is necessary

    ki = order(i);     % initial guess of the order of the variable X_i in the
                       % response surface

    no_fevals = 0;
    
    for ki = order(i):5		% loop over possible polynomial orders

        no_pts = ki+5;          % number of sample points (must be > ki+1)

	% sample points along dimension X_i within the domain [-1,1]
	% (roots of no_pts-th order Chevbyshev polynomial)
        x = [cos(pi*([1:no_pts]-0.5)./no_pts)']; 
                                
	% allocate memory for matrix of sampling points along all variables
	% (zeros represent sampling  points taken at mean values,
	% since all variables are standardized normal random variables)
        xmap = zeros(no_pts,n);

	% the sample points along X_i are mapped onto the domain [-h_ord,h_ord]
	% only the i-th column is non-zero since all other variables 
	% are kept constants at their mean values 
        xmap(:,i) = x*h_ord;

        g = zeros(no_pts,1);    % allocate memory for limit state values

        % evaluate the limit state function at the Chybyshev sampling points
        for k = 1:no_pts
            gg   = feval(limit_state,xmap(k,:));
            g(k) = gg(g_idx); 
        end
        no_fevals = no_fevals+no_pts;% increment number of function evaluations

	% values of 0-th to ki-th degree Chevbyshev polynomials at
	% the sampling points.
	% this matrix is used for the determination of coefficients 
        Tx = cos(acos(x)*(0:ki));   
                                    
        d = Tx'*g./diag(Tx'*Tx);	% coefficients by least squares method
                
	residuals = g-Tx*d;		% residuals of the 1-D curve-fit
	fit_error = norm(residuals) / norm(g);	% error of the curve fit
	R2_order(i) = 1 - ( norm(residuals) / norm(g - mean(g)) )^2;
                                    
        figure(1)
         clf
         plot ( x*h_ord,g,'ob', x*h_ord,Tx*d,'*r')
         ttl = sprintf('k_{%d}=%d; R^2 = %f; fit-error = %f;', i, ki, R2_order(i),fit_error);
         title(ttl);
         pause(3)

	if ( R2_order(i) > 1-tol && fit_error < tol ) % the 1D fit is accurate
	    break;
        end

    end

    order(i) = ki;		% save the identified polynomial order

    order_time(i) = toc;	% record the time required 

  end

end	            % 1: do order determination; 0: skip order determination

% output the results
fprintf('  Variable     Determined Order    Time (sec)     R_sq \n');
for i = 1:n
    fprintf('%10.0f %20.0f %11.3f       %9.6f\n', ...
           i, order(i), order_time(i), R2_order(i) );
end
fprintf('  %41.3f seconds \n\n', sum(order_time) );


% 2nd Stage: Mixed Terms Determination ------------------------------------
%            All the terms in the response surface, especially the mixed 
%            terms, are formulated.

fprintf('2nd Stage: Mixed Terms Determination ...\n');

 % The matrix 'term' indicates which mixed term power-products are present
 % in the response surface function.  
 % Each element in this matrix represents the power of a variable in a term.
 % The columns of 'term' correspond to individual variables. 
 % Each row indicates the powers present in the prodcut of a term.

 Rows = prod(order+1);

 term = zeros(Rows,n);              % allocate memory for the 'term' matrix
 row = zeros(1,n); row(1) = -1;     % starting value for the first row

 for r = 1:Rows                     % loop over all rows
     row(1) = row(1)+1;             % increment first column
     for c = 1:n                    % check every column in the row
         if row(c) > order(c)
            row(c) = 0;
            row(c+1) = row(c+1)+1;  % increment columns as needed
         end
         term(r,:) = row;           % save the row in the 'term' matrix
      end
 end

 % The power of a variable in each term can not be greater than 
 % the order of that variable alone.  
 % Remove the terms in which the total order is larger than 
 % the highest order term.  

 term = term( find(sum(term,2) <= max(order)) , : );

 % The number of rows is in the matrix 'term' is the number of
 % required terms in the response surface. 

 no_term = size(term,1); 

fprintf('  Number of Mixed Terms: %2d\n', no_term-sum(order)-1);
fprintf('  Total Number of Terms: %2d\n\n', no_term);

% 3rd Stage: Response Surface Approximation -------------------------------
%            The coefficients of the response surface are estimated by
%            the ordinary least squares method or singular value decomposition

fprintf('3rd Stage: Response Surface Approximation ...\n');

 tic     % start recording time required for response surface approximation

 % random sampling points are taken uniformly in the domain [-h_reg, h_reg] 
 x = (rand(2.0*no_term,n)-0.5)*2*h_reg;

 m = size(x,1);  % number of samples for response surface fit

 g = zeros(m,1); % initialize a vector which holds the limit state values

 fprintf('  Number of Limit State Evaluations: %2d\n',2*no_term);
                            
 fprintf('  Limit State Evaluation in Progress ');
                            
 for k = 1:m
     gg = feval(limit_state,x(k,:)); 
     g(k,1) = gg(g_idx); % evaluate limit state values at the sampling points
     fprintf('...%10d',k);
     fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b');
     no_fevals = no_fevals+1; % number of function evaluations
 end
 fprintf('... done\n');

 % initialize a matrix used for the determination of the coefficients
 % ... the 'system matrix'  
 %
 % g = X * coef, where g is a vector of the limit state function 
 % evaluated at various sample points, and 
 % coef is the vector of coefficients in the response surface.
 X = ones(m,no_term);

 % in the matrix X, 
 % columns correspond to terms in the response surface polynomial and 
 % rows correspond to the sampling points
 for j = 1:no_term
    X(:,j) = [prod(x.^(ones(m,1)*term(j,:)),2)];
 end

 % the first column are all ones,
 % since the first term in the response surface is a constant
 X(:,1) = ones(m,1);

 % determine the coefficients of the response surface
%coeff = (X'*X)\(X'*g);    % ... by the ordinary least squares method
 coeff = X \ g;            % ... by singular value decomposition

 % statistical analysis of coefficients 
 condX = cond(X);			% condition number
 residuals = (g - X*coeff);		% residuals of the fit
 R2 = 1 - ( norm(residuals) / norm(g - mean(g)) )^2;	% R-squared error
 Adj_R2 = ((m-1) * R2 - no_term ) / (m - (no_term+1)); % Adjusted R2

 % Standard error of the fit coefficients
 Std_Err_Coeff = sqrt( (residuals'*residuals) * diag(inv(X'*X))/(m-no_term+1) );

 RS_time = toc;  % record the time require


 fprintf('  condition number  = %6.1f \n\n', condX );
 fprintf('  Response Surface Coefficients \n' );
 fprintf('   i     c(i)   ');
 for j=1:n 
	fprintf(' z(%d)', j);
 end
 fprintf('  StdErr \n');

 for i=1:no_term
	fprintf('  %2d  %8.4f ', i, coeff(i) );
 	for j=1:n 
		fprintf(' %2d  ', term(i,j) );
	end
 	fprintf(' %8.4f \n', Std_Err_Coeff(i) );
 end
 fprintf('\n');

 fprintf('  R-square          = %6.3f \n', R2 );
 fprintf('  Adjusted R-square = %6.3f \n',  Adj_R2 );
 fprintf('  Time ... %32.3f seconds\n\n',RS_time);

 clear x g X
                        

% 4th Stage: Monte Carlo Simulation ---------------------------------------
%            A full scale Monte Carlo simulation is run on the response surface
%            in order to estimate the probability of failure

fprintf('4th Stage: Monte Carlo Simulation ...\n');
fprintf('  Sample Size: %6d\n',N);

 tic                         % start recording the time required for MCS

 if nargin < 8
     Zmcs = randn(N,n);             % normal random sampling points
 else
     if ~all( size(Zmcs) == [ N n ] )
         fprintf(' *** The size of Zmcs must be N by n ***')
         help HOSRSM
         fprintf(' *** The size of Zmcs must be N by n ***')
         return
     end
 end

 g = zeros(1,N);                   % allocate memory for limit state evaluations

 BlockSize = 1e4;                  % number of points in a data block
 NumBlocks = floor(N/BlockSize);   % number of blocks

 if NumBlocks < 10		% run all samples in a single block

    X = ones(N,no_term);     % allocate memory for matrix used for limit state

    for j = 1:no_term
        X(:,j) = prod(Zmcs.^(ones(N,1)*term(j,:)),2);
    end

    g = X*coeff;             % estimate limit state values from the R.S.

    clear X

 else                        % very many samples, run them in blocks ...

    for b=1:NumBlocks        % in order to avoid swapping memory to/from disk

        X = ones(BlockSize,no_term);

        ZmcsB = Zmcs((b-1)*BlockSize+1:b*BlockSize,:); 
   
        for j = 1:no_term
            X(:,j) = prod(ZmcsB.^(ones(BlockSize,1)*term(j,:)),2);
        end

        g((b-1)*BlockSize+1:b*BlockSize) = X*coeff;

        clear X ZmcsB

    end

    fprintf('  %6d Blocks of Size %6d\n', NumBlocks, BlockSize );

 end

   
 Pf = length(find(g<=0))/N;  % probability of failure, using Monte-Carlo on R.S.
   
 MCS_time = toc; % record the time required

fprintf('  Time ... %37.3f seconds\n',MCS_time);
fprintf('  Number of Limit State Evaluations ... %4d (total) \n', no_fevals );
fprintf('  Probability of Failure ... %24.4e \n\n', Pf);

% updated 1.29.2006, 2.21.2007, 3.6.2007, 10.14.2009
% --------------------------------------------------------------------- HOSRSM
