Skip to content
Snippets Groups Projects
Commit 067bf8f0 authored by jcandrsn's avatar jcandrsn
Browse files

More work on adapting the pce coefficient calc to generic response functions. ...

More work on adapting the pce coefficient calc to generic response functions.  Still limited to scalar response functions.  Renamed obs.m and dobs.m and split up the state variable and parameter Jacobian between two .m files.
parent 601d687a
No related branches found
No related tags found
1 merge request!1Development
......@@ -217,8 +217,9 @@ data.P = max_poly_order;
% function for repeated use.
data.wtd_psi_mat = diag(wts)*psi_mat;
data.fhan = @obs;
data.dfhan = @dobs;
data.rhan = @x10;
data.drdxhan = @x10_drdx;
data.drdphan = @x10_drdp;
% Nt = Total Number of columns of Psi matrix and number of terms in the PCE
% expansion. Casting to int64 because the factorial function was returning
......@@ -252,7 +253,7 @@ 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', ...
@uq_pce_coefficients_dU, data, 'zero', ...
'uidx', vertcat(xbp_idx{:}, p_idx{:}), ...
'u0', alpha_ig);
......
function dr = dobs(x, p)
dr = eye(size(x,1));
end
function r = obs(x, p)
r = x(1,1);
end
function r = x10(x, p)
% A = zeros(size(x));
% A(1,1) = 1;
% r = A*x;
r(1) = x(1,1);
end
function Jp = x10_drdp(x, p)
% for example problem, x is 2x300, p is 3x1
% r is 1x1
Jp = sparse(zeros(size(p)));
end
\ No newline at end of file
function Jx = x10_drdx(x, p)
% for example problem, x is 2x300, p is 3x1
% r is 1x1
% A = zeros(size(x));
% A(1,1) = 1;
% r = A*x;
% J = A;
Jx = sparse(zeros(size(x)));
Jx(1) = 1;
end
\ No newline at end of file
......@@ -3,22 +3,21 @@ function [data, y] = uq_pce_coefficients(prob, data, u)
% Evaluate Response function and calculate coefficients of polynomial
% chaos expansion.
% Separate out components
% Need to generalize this to accept entire trajectories and parameter sets
% Current thought:
for i=1:data.M % M is the number of quadrature points
% Separate out components
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);
r1 = data.rhan(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);
r(i,:) = data.rhan(x, p);
end
end
% last_state = data.M^data.s;
......@@ -26,6 +25,7 @@ end
% 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;
......
function [data, J] = uq_pce_coefficients_dU(prob, data, u)
%UQ_TEMP Summary of this function goes here
% Detailed explanation goes here
% u contains all Basis point and parameter values for M trajectories as
% well as continuation variables for the alphas where M is the numerical
% 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);
% J_alpha = eye(size(alpha);
% J = [J_x J_p J_alpha];
%
pidx_shift = data.pidx - data.xidx(end);
dr.dx = sparse(1,data.xidx(end));
dr.dp = sparse(1,pidx_shift(end));
% 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));
for i=data.M:-1:1 % M is the number of quadrature points, going backwards
% to keep Matlab from squawking about memory allocation
% Cycling through like this to account for a potential case where the
% trajectories are allowed to have discretizations and thus will have
% varying numbers of basis point values.
% Separate out components
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));
dr(i).dx = sparse(1, data.xidx(i,1):data.xidx(i,2), ...
reshape(data.drdxhan(x,p),1,[]), 1, data.xidx(end));
dr(i).dp = sparse(1, pidx_shift(i,1):pidx_shift(i,2), ...
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
speye(data.P+1)]; % alphas are linear
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment