Updates to uq_pce_coefficients, obs, dobs that allow for more generic response...

Updates to uq_pce_coefficients, obs, dobs that allow for more generic response functions.  uq_pce_coefficients_dU needs updated and included in demo_uq to speed up Jacobian evaluation.

matlab_development/uq/examples/linear_harmonic_oscillator/demo_uq.m

@@ -29,7 +29,7 @@ end
 mu_k = 3;
 sig_k = 0.2;
 
-quadrature_order = 3;
+quadrature_order = 9;
 stochastic_dimension = 1;
 orthogonal_polynomial_type = 'Hermite';
 
@@ -67,6 +67,7 @@ if not(coco_exist(run_name, 'run'))
     end
 
     sol = coll_read_solution('orig', last_run_name, lab);
+    % Reset Time
     t0 = sol.tbp - sol.tbp(1);
     x0 = sol.xbp;
     p0 = sol.p;
@@ -116,16 +117,20 @@ prob = coco_prob();
 prob = coco_set(prob, 'ode', 'autonomous', false);
 prob = coco_set(prob, 'coll', 'NTST', NTST, 'NCOL', NCOL);%, 'var', true);
 prob = coco_set(prob, 'cont', 'PtMX', 1000, 'h', 0.5, 'almax', 30);
-    
-
-t0 = cell(size(k_vals, 1),1);
-x0 = cell(size(k_vals, 1),1);
-p0 = cell(size(k_vals, 1),1);
-data = cell(size(k_vals, 1),1);
-uidx = cell(size(k_vals, 1),1);
-maps = cell(size(k_vals, 1),1);
-
-for i=1:size(k_vals, 1)
+ 
+t0 = cell(size(nds, 1),1);
+x0 = cell(size(nds, 1),1);
+p0 = cell(size(nds, 1),1);
+data = cell(size(nds, 1),1);
+uidx = cell(size(nds, 1),1);
+maps = cell(size(nds, 1),1);
+xbp_idx = cell(size(nds, 1), 1);
+xbp_shp = cell(size(nds, 1), 1);
+p_idx = cell(size(nds, 1), 1);
+
+% This for loop only works for a single stochastic parameter.  Need to
+% generalize it to work for more.
+for i=1:size(nds, 1)
     k_lab = strcat('UQ',int2str(i));
     lab = coco_bd_labs(bd, k_lab);
 
@@ -143,6 +148,9 @@ for i=1:size(k_vals, 1)
                                             strcat(pname,'.coll'), ...
                                             'data', 'uidx');
     maps{i} = data{i}.coll_seg.maps;
+    xbp_idx{i} = uidx{i}(maps{i}.xbp_idx);
+    xbp_shp{i} = maps{i}.xbp_shp;
+    p_idx{i}  = uidx{i}(maps{i}.p_idx);
 
     % Add Boundary Conditions
     prob = coco_add_func(prob, strcat(pname,'.po'), @linode_het_bc_v0, ...
@@ -162,13 +170,13 @@ for i=1:size(k_vals, 1)
 end
 
 % Collect indices for unequal stiffness values
-p_idx = [uidx{1}(maps{1}.p_idx)];
+p_idx_2 = [uidx{1}(maps{1}.p_idx)];
 for j=2:size(uidx,1)
-    p_idx = [p_idx; uidx{j}(maps{j}.p_idx([1,2]))];
+    p_idx_2 = [p_idx_2; uidx{j}(maps{j}.p_idx([1,2]))];
 end
-p_names = cell(1,size(p_idx,1));
+p_names = cell(1,size(p_idx_2,1));
 p_names(1:3) = {'k1', 'phi1', 'om'};
-for j = 4:size(p_idx,1)
+for j = 4:size(p_idx_2,1)
     if mod(j,2)
         p_names{j} = strcat('phi',int2str((j-1)/2));
     else
@@ -177,7 +185,7 @@ for j = 4:size(p_idx,1)
 end
 % Need to Split phis into separate 'active' parameters so they don't get 
 % displayed
-prob = coco_add_pars(prob, 'pars', p_idx, p_names, 'inactive');
+prob = coco_add_pars(prob, 'pars', p_idx_2, p_names, 'inactive');
 
 %% Step 6
 %  Introduce a monitor function to construct the PCE from the desired
@@ -186,9 +194,9 @@ prob = coco_add_pars(prob, 'pars', p_idx, p_names, 'inactive');
 % Observable function is the initial position x0 which, by construction, is
 % the peak amplitude of the curve.
 % Collect the initial position indices from the various trajectories
-for i=1:size(uidx,1)
-    x_0{i} = uidx{i}(maps{i}.x0_idx(1));
-end
+% for i=1:size(uidx,1)
+%     x_0{i} = uidx{i}(maps{i}.x0_idx(1));
+% end
 
 % Set up a data structure to pass to the uncertainty quantification
 % problem.  Since the psi matrix is based on the standard random variables,
@@ -196,10 +204,12 @@ end
 data = struct();
 data.M = quadrature_order;
 data.s = stochastic_dimension;
+data.xbp_idx = xbp_idx;
+data.xbp_shp = xbp_shp;
+data.p_idx = p_idx;
 
 max_poly_order = 10;
-psi_mat = uq_make_psi_mat(nds, max_poly_order, ...
-                               orthogonal_polynomial_type);
+psi_mat = uq_make_psi_mat(nds, max_poly_order, orthogonal_polynomial_type);
 data.P = max_poly_order;
 
 % The quadrature weights are similarly identical for each continuation
@@ -226,14 +236,30 @@ end
 % linear in the alphas and should converge after one step of the Newton
 % Iteration regardless.
 alpha_ig = zeros(size(alphas));
+data.xidx = ones(quadrature_order,2);
+for i=1:size(data.xbp_idx,1)
+    data.xidx(i,2) = size(xbp_idx{i},1);
+end
+data.xidx(:,2) = cumsum(data.xidx(:,2));
+data.xidx(2:end,1) = data.xidx(2:end,1) + data.xidx(1:end-1,2);
+
+data.pidx = ones(quadrature_order, 2);
+for i=1:size(data.p_idx,1)
+    data.pidx(i,2) = size(p_idx{i},1);
+end
+data.pidx(:,2) = cumsum(data.pidx(:,2));
+data.pidx(2:end,1) = data.pidx(2:end,1) + data.pidx(1:end-1,2);
+data.pidx = data.pidx + data.xidx(end);
 
 prob = coco_add_func(prob, 'pce', @uq_pce_coefficients, ...
-                     data, 'zero', 'uidx', [x_0{:}]', ...
+                     data, 'zero', ...
+                     'uidx', vertcat(xbp_idx{:}, p_idx{:}), ...
                      'u0', alpha_ig);
+
 [alpha_idx, alpha_data] = coco_get_func_data(prob, 'pce', 'uidx', 'data');
-% Grab the last Nt parameters and associate continuation parameters with
-% them.
+% Grab the last Nt parameters
 alpha_idx = alpha_idx(end-alpha_data.Nt+1:end);
+% Associate continuation parameters with them.
 % prob = coco_add_pars(prob, 'alphas', alpha_idx, alphas, 'active');
 
 %% Step 7
@@ -246,20 +272,20 @@ prob = coco_add_func(prob, 'PCE_mean', @uq_pce_mean, ...
 prob = coco_add_func(prob, 'PCE_variance', @uq_pce_variance, ...
                      @uq_pce_variance_dU, {}, 'inactive', 'variance', ...
                      'uidx', alpha_idx);
-          
-% ps = cell(1,Nt+quadrature_order+1);
+
 ps = cell(1,quadrature_order+1);
 ps{1} = 'om';
-% ps{2} = 'alpha0';
-% ps{3} = 'alpha1';
-% for i=1:Nt
-%     ps{i+1} = strcat('alpha', int2str(i-1));
-% end
 for i=1:quadrature_order
     % Phase Angles
     ps{i+1} = strcat('phi',int2str(i));
 end
 
+% ps = cell(1,quadrature_order);
+% for i=1:quadrature_order
+%     % Phase Angles
+%     ps{i} = strcat('phi',int2str(i));
+% end
+
 bd = coco(prob, run_name, [], 1, [ps, {'mean', 'variance'}], [0.1, 10]);
 %% 
 

matlab_development/uq/examples/linear_harmonic_oscillator/dobs.m

-function u = obs(u)
-u = ones(size(u))
+function dr = dobs(x, p)
+    dr = eye(size(x,1));
 end
 

matlab_development/uq/examples/linear_harmonic_oscillator/obs.m

-function u = obs(u)
-
+function r = obs(x, p)
+    r = x(1,1);
 end
 

matlab_development/uq/examples/linear_harmonic_oscillator/plot_sweep_uq.m

+M = 9;
+P = 10;
+run_name = strcat('UQ_M=',int2str(M),'_P=',int2str(P));
+fprintf(run_name);
+fprintf('\n')
+bd = coco_bd_read(run_name);
+labs = coco_bd_labs(bd);
+data_mat = zeros(size(labs,2), M+1);
+for lab=labs
+    for i=1:M
+        curve_name = strcat('UQ', int2str(i));
+        sol = coll_read_solution(curve_name, run_name, lab);
+        data_mat(lab, i+1) = sol.xbp(1,1);
+    end
+    data_mat(lab,1) = sol.p(3);
+end
\ No newline at end of file

matlab_development/uq/toolbox/uq_pce_coefficients.m

@@ -5,16 +5,29 @@ function [data, y] = uq_pce_coefficients(prob, data, u)
 
 % Separate out components
 % Need to generalize this to accept entire trajectories and parameter sets
-last_state = data.M^data.s;
-x = u(1:last_state);
-alphas = u(last_state+1:end);
-
-r = data.fhan(x);
-y = alphas - data.wtd_psi_mat'*r;
-
-% mu = alphas(1);
-% var = sum(alphas(2:end).^2,2);
-% 
-% y = [mu; var];
+% Current thought:
+for i=1:data.M   % M is the number of quadrature points
+     x = u(data.xidx(i,1):data.xidx(i,2));
+     x = reshape(x, data.xbp_shp{i});
+     p = u(data.pidx(i,1):data.pidx(i,2));
+     if i == 1
+         % Set the size of the response array on the first trip through
+         % the loop and assign the first value
+         r1 = data.fhan(x, p);
+         r = zeros(data.M, size(r1,2));
+         r(i,:) = r1;
+     else
+         % Only Calculate the function value and assign it to correct spot
+         r(i,:) = data.fhan(x, p);
+     end
+end
+% last_state = data.M^data.s;
+% x = u(1:last_state);
+% alphas = u(last_state+1:end);
+% p = [];
+% r = data.fhan(x, p);
+alphas = u(data.pidx(end)+1:end);
+R = data.wtd_psi_mat'*r;
+y = alphas - R;
 
 end

matlab_development/uq/toolbox/uq_pce_coefficients_dU.m

-function [data, alphas] = uq_pce_coefficients_dU(prob, data, u)
+function [data, J] = uq_pce_coefficients_dU(prob, data, u)
 %UQ_TEMP Summary of this function goes here
 %   Detailed explanation goes here
 
-drdu = data.dfhan(u);
-alphas = data.wtd_psi_mat'*drdu;
 
-% mu = alphas(1);
-% var = sum(alphas(2:end).^2,2);
-% 
-% y = [mu;
-%      var];
+% Separate out components
+% Need to generalize this to accept entire trajectories and parameter sets
+% general idea:
+last_state = data.M^data.s;
+x = u(1:last_state);
+alphas = u(last_state+1:end);
+p = [];
+J = zeros(size(alphas,1), size(u,1));
+drdx = data.dfhan(x, p);
+J(:,1:last_state) = -data.wtd_psi_mat'*drdx;
+J(:,last_state+1:end) = eye(size(alphas,1));
 
 end