Modified output of uq_gauss_nodes to be similarly shaped to other coco...

Modified output of uq_gauss_nodes to be similarly shaped to other coco functions (rows are different parameters).  uq_make_psi_mat and other functions modified to accomodate this change.  Kept original functions as *_orig.m files and confirmed correct operation with transposing_nds_test.m

matlab_development/uq/toolbox/transform_uniform_sample.m

@@ -23,7 +23,7 @@ function sample = transform_uniform_sample(sample, mu, CoV)
 %  array([[ 1.8660254 ,  2.73205081,  3.59807621,  4.46410162],
 %         [ 5.33012702,  6.19615242,  7.06217783,  7.92820323]])
 %  Matlab test:
-%  >>transform_uniform_sample(a,1,0.5)
+%  >>transform_uniform_sample([[1,2,3,4];[5,6,7,8]]),1,0.5)
 % 
 %  ans =
 % 

matlab_development/uq/toolbox/transposing_nds_test.m

+quadrature_order=2;
+stochastic_dimension = 2;
+pce_order = 1;
+orthogonal_polynomial_type = 'Legendre';
+
+[nds, wts] = uq_gauss_nodes_orig(quadrature_order, ...
+                                stochastic_dimension, ...
+                                orthogonal_polynomial_type);
+psi_mat = uq_make_psi_mat_orig(nds, pce_order, orthogonal_polynomial_type);
+diag(wts)*psi_mat
+
+[nds, wts] = uq_gauss_nodes(quadrature_order, ...
+                                stochastic_dimension, ...
+                                orthogonal_polynomial_type);
+psi_mat = uq_make_psi_mat(nds, pce_order, ...
+                                  orthogonal_polynomial_type);
+
+psi_mat*diag(wts)
\ No newline at end of file

matlab_development/uq/toolbox/uq_gauss_nodes.m

@@ -27,8 +27,8 @@ for i=1:n
     nds(i,:) = reshape(nd_cells{i},1,[]);
     wts(i,:) = reshape(wt_cells{i},1,[]);
 end
-nds = nds';
-wts = prod(wts,1)';
+
+wts = prod(wts,1);
 
 end
 

matlab_development/uq/toolbox/uq_gauss_nodes_orig.m

+function [nds, wts] = uq_gauss_nodes(m, n, type)
+% Given a quadrature order, m (integer, is the same in each dimension), 
+% a number of , n (integer), and a quadrature rule type, return the tensor 
+% product of the quadrature node locations, nds, and appropriately
+% multiplied quadrature weight, wts.
+
+switch type
+    case {'Legendre','Le','LeN'}
+        [nds, wts] = gauss_legendre_nodes(m);
+    case {'Hermite','He','HeN'}
+        [nds, wts] = gauss_hermite_nodes(m);
+    otherwise
+        warning(strcat('Specified Type not supported, Providing',...
+                       ' weights and nodes for Gauss-Legendre',...
+                       ' quadrature'))
+        [nds, wts] = gauss_legendre_nodes(m);
+end
+
+nd_cells = repmat({nds}, 1,n);
+wt_cells = repmat({wts}, 1,n);
+[nd_cells{:}] = ndgrid(nd_cells{:});
+[wt_cells{:}] = ndgrid(wt_cells{:});
+
+nds = zeros(n, numel(nd_cells{1}));
+wts = zeros(n, numel(wt_cells{1}));
+for i=1:n
+    nds(i,:) = reshape(nd_cells{i},1,[]);
+    wts(i,:) = reshape(wt_cells{i},1,[]);
+end
+nds = nds';
+wts = prod(wts,1)';
+
+end
+
+function [nds, wts] = gauss_legendre_nodes(m)
+n = (1:m-1)';
+g = n.*sqrt(1./(4*n.^2-1)); 
+J = -diag(g,1)-diag(g,-1);
+[w, x] = eig(J);
+nds = diag(x)';
+% Normalizing the weights so they add up to one.
+% The weight function is 1/2 for a uniform distribution between -1 and 1
+% instead of the typical uniform distribution between 0 and 1.
+wts = (2*w(1,:).^2)/2;
+end
+
+function [nds, wts] = gauss_hermite_nodes(m)
+n = (1:m-1)';
+g = sqrt(n); 
+J = diag(g,1)+diag(g,-1);
+[w, x] = eig(J);
+[nds, idx] = sort(diag(x));
+nds = nds';
+wts = sqrt(2*pi)*w(1,:).^2;
+% Normalizing the weights so they add up to one
+% The weight function is the pdf of the standard normal distribution
+% (1/sqrt(2*pi))*exp((-x.^2)/2).
+wts = wts(idx)/(sqrt(2*pi));
+end
\ No newline at end of file

matlab_development/uq/toolbox/uq_linsolve_pce_coefficients.m

-function [data, alphas] = uq_pce_coefficients(prob, data, u)
-%UQ_TEMP Summary of this function goes here
-%   Detailed explanation goes here
-
-M = data.M;
-s = data.s;
-Nt = data.Nt;
-
-
-% Split out the alphas from the u's
-alphas = 
-alphas = data.wtd_psi_mat'*r;
-
-% mu = alphas(1);
-% var = sum(alphas(2:end).^2,2);
-% 
-% y = [mu;
-%      var];
-
-end
-

matlab_development/uq/toolbox/uq_make_psi_mat.m

@@ -8,7 +8,7 @@ function psi = uq_make_psi_mat(sample, max_order, poly_type)
 %  Number of Quadrature points, M = size(sample, 1) = Quadrature_Order^s
 %  psi is an M X Nt) matrix where Nt = (s + P)! / (s! P!).
 %  Quadrature_Order is set independently of this function.
-
+sample = sample';
 if nargin < 3
     poly_type = 'Legendre';
 end
@@ -43,3 +43,4 @@ psi = psi(:,idx)';
 psi = reshape(psi,d3,d4,d1);
 psi = prod(psi,2);
 psi = permute(psi,[3,1,2]);
+psi = psi';

matlab_development/uq/toolbox/uq_make_psi_mat_orig.m

+function psi = uq_make_psi_mat(sample, max_order, poly_type)
+%  Given a random sample, sample, and a max polynomial order, max_order,
+%  return the combination of polynomials whose product is lower or equal
+%  order than max_order.  poly_type defines the basis functions to use.
+%  ------------------------------------------------------------------------
+%  max_order = P = Maximum Total Polynomial Order (as described above)
+%  Number of Random Variables Describing the system, s = size(sample, 2)
+%  Number of Quadrature points, M = size(sample, 1) = Quadrature_Order^s
+%  psi is an M X Nt) matrix where Nt = (s + P)! / (s! P!).
+%  Quadrature_Order is set independently of this function.
+
+if nargin < 3
+    poly_type = 'Legendre';
+end
+
+n = size(sample,2);
+for k = 1:n
+    ps{k} = 0:max_order;
+end
+grid_cells = cell(size(ps));
+
+
+[grid_cells{:}] = ndgrid(ps{:});
+
+grids = zeros(n, numel(grid_cells{1}));
+for i=1:n
+    grids(i,:) = reshape(grid_cells{i},1,[]);
+end
+idx = sum(grids,1) <= max_order;
+grids = grids(:, idx)';
+[~,idx] = sort(sum(grids,2));
+grids = grids(idx,:);
+idx = max_order*(0:(n-1)) + (1:n);
+idx = grids + repmat(idx,size(grids,1),1);
+
+[d1, d2] = size(sample);
+[d3, d4] = size(idx);
+
+psi = uq_orthogonal_poly_vals_orig(sample, poly_type, max_order, 1);
+psi = permute(psi, [1,3,2]);
+psi = reshape(psi, [], d2*(max_order+1));
+psi = psi(:,idx)';
+psi = reshape(psi,d3,d4,d1);
+psi = prod(psi,2);
+psi = permute(psi,[3,1,2]);

matlab_development/uq/toolbox/uq_orthogonal_poly_vals_orig.m

+function psi = uq_orthogonal_poly_vals(x, type, max_order, norm)
+% For given evaluation locations, x, return a matrix of orthogonal
+% polynomial values of type, type, up to order, max_order.  The value of
+% norm determines whether the coefficients are normalized such that the
+% norm of a particular polynomial order equals one or not.  
+[d1, d2] = size(x);
+psi = zeros(d1, d2, max_order+1);
+psi(:,:,1) = 1;
+psi(:,:,2) = x;
+
+if nargin < 4
+    norm = 1;
+end
+
+switch type
+  case {'Legendre', 'Le'}
+    if norm
+      wts = permute(repmat(sqrt(2*(0:max_order)+1), d1, 1, d2), [1,3,2]);
+    else
+      wts = ones(d1, d2, max_order+1);
+    end
+      if max_order==0
+          % Special code for trivial case to get matrix dimensions to agree
+          psi = psi(:,:,1);
+      else    
+          for i=3:(max_order+1)
+              n = i-2;
+              % Three term recurrence for Legendre Orthogonal Polynomials
+              psi(:,:,i) = ((2*n+1)*x.*psi(:,:,i-1) - n*psi(:,:,i-2))/(n+1);
+          end
+      end
+
+  case {'Hermite', 'He'}
+    if norm
+      wts = permute(repmat(sqrt(factorial(0:max_order)).^-1, d1, 1, d2),...
+                    [1,3,2]);
+    else
+      wts = ones(d1, d2, max_order+1);
+    end
+      
+      if max_order==0
+          % Special code for trivial case to get matrix dimensions to agree
+          psi = psi(:,:,1);
+      else
+          % Three term recurrence for Hermite Orthogonal polynomials
+          for i=3:(max_order+1)
+              psi(:,:,i) = x.*psi(:,:,i-1) - (i-2)*psi(:,:,i-2);
+          end
+      end
+end
+
+psi = wts.*psi;
+
+end
\ No newline at end of file

matlab_development/uq/toolbox/uq_pce_coefficients.m

@@ -28,7 +28,7 @@ for i=1:data.M   % M is the number of quadrature points
 end
 
 alphas = u(data.pidx(end)+1:end);
-R = data.wtd_psi_mat'*r;
+R = data.wtd_psi_mat*r;
 y = alphas - R;
 
 end

matlab_development/uq/toolbox/uq_pce_coefficients_dU.m

@@ -5,9 +5,9 @@ function [data, J] = uq_pce_coefficients_dU(prob, data, u)
 %   integration order of the PCE.
 %   u = [ xbp1, xbp2, ... xbpM, p1, p2, ... pM ,alpha];
 %
-%   y = alpha - Psi'*r(x,p);
-%   J_x = dy/dx = (dy/dr)*(dr/dx) = -Psi'*dr/dx(x,p);
-%   J_p = dy/dp = (dy/dr)*(dr/dp) = -Psi'*dr/dp(x,p);
+%   y = alpha - w*Psi*r(x,p);
+%   J_x = dy/dx = (dy/dr)*(dr/dx) = -w*Psi*dr/dx(x,p);
+%   J_p = dy/dp = (dy/dr)*(dr/dp) = -w*Psi*dr/dp(x,p);
 %   J_alpha = eye(size(alpha);
 %   J = [J_x J_p J_alpha];
 %
@@ -35,8 +35,8 @@ for i=data.M:-1:1   % M is the number of quadrature points, going backwards
                       reshape(data.drdphan(x,p),1,[]), 1, pidx_shift(end));
 end
 
-J = [sparse(-1*data.wtd_psi_mat'*vertcat(dr.dx)) ... % Psi' * dr/dx
-     sparse(-1*data.wtd_psi_mat'*vertcat(dr.dp)) ... % Psi' * dr/dp
+J = [sparse(-1*data.wtd_psi_mat*vertcat(dr.dx)) ...  % w*Psi * dr/dx
+     sparse(-1*data.wtd_psi_mat*vertcat(dr.dp)) ...  % w*Psi * dr/dp
      speye(data.P+1)];                               % alphas are linear
 
 end