Added some methods in the testing area

testing/gpc/unittest_gpc_moments.m

@@ -0,0 +1,21 @@
+function unittest_gpc_moments
+% UNITTEST_GPC_MOMENTS Test the GPC_MOMENTS function.
+%
+% Example (<a href="matlab:run_example unittest_gpc_moments">run</a>)
+%   unittest_gpc_moments
+%
+% See also GPC_MOMENTS, TESTSUITE 
+
+%   Elmar Zander
+%   Copyright 2013, Inst. of Scientific Computing, TU Braunschweig
+%
+%   This program is free software: you can redistribute it and/or modify it
+%   under the terms of the GNU General Public License as published by the
+%   Free Software Foundation, either version 3 of the License, or (at your
+%   option) any later version. 
+%   See the GNU General Public License for more details. You should have
+%   received a copy of the GNU General Public License along with this
+%   program.  If not, see <http://www.gnu.org/licenses/>.
+
+munit_set_function( 'gpc_moments' );
+

testing/gpc/unittest_gpcindex_combine.m

@@ -19,7 +19,64 @@
 
 munit_set_function( 'gpcindex_combine' );
 
+% the gpcindex_combine function is not yet written, but the tests it should
+% fulfill are already here.
 
+function foo
+p=4;
+I1 = multiindex(3,2)
+I2 = multiindex(2,4)
+Iz=multiindex_product({I1, I2})
+Iz0=multiindex_product_maxp({I1, I2},0)
+Iz1=multiindex_product_maxp({I1, I2},1)
+Iz2=multiindex_product_maxp({I1, I2},2)
+Iz3=multiindex_product_maxp({I1, I2},3)
+Iz4=multiindex_product_maxp({I1, I2},4)
+Iz5=multiindex_product_maxp({I1, I2},5)
+Iz6=multiindex_product_maxp({I1, I2},6)
+1;
+
+%different
+
+% product
+function I = multiindex_sum(Is)
+I1 = Is{1};
+I2 = Is{2};
+n1 = size(I1,1);
+n2 = size(I2,1);
+
+[ind1,ind2] = meshgrid(1:n1, 1:n2);
+I = [I1(ind1,:), I2(ind2,:)];
+
+function I = multiindex_product(Is)
+I1 = Is{1};
+I2 = Is{2};
+n1 = size(I1,1);
+n2 = size(I2,1);
+
+[ind1,ind2] = meshgrid(1:n1, 1:n2);
+I = [I1(ind1,:), I2(ind2,:)];
+
+% product_maxp
+function I=multiindex_product_maxp(Is, p)
+I1 = Is{1};
+I2 = Is{2};
+d1 = multiindex_order(I1);
+d2 = multiindex_order(I2);
+[D1,D2] = meshgrid(d1, d2);
+ind=(D1+D2<=p);
+[ind1,ind2] = meshgrid(1:size(I1,1), 1:size(I2,1));
+I = [I1(ind1(ind),:), I2(ind2(ind),:)];
+
+% sort 
+
+% sum
+
+% union, 
+
+
+return
+function foo
 I1a = multiindex(3,2);
 I1b = multiindex(3,4);
 I2 = multiindex(2,4);
@@ -30,12 +87,6 @@
 gpcindex_combine( {'H', I1a}, {'H', I1b}, 'product_mp');
 
 
-
-% gpcspace_create
-% gpcspace_product
-% gpcspace_direct_sum
-
-
 % 1. add two random variables on same V1a, V1b subset V1
 % 2. add two random variables on V1, V2 independent, minimum V3
 % 3. same as 2, "full" V3
@@ -55,3 +106,6 @@
 % gpc_integrate
 
 
+% gpcspace_create
+% gpcspace_product
+% gpcspace_direct_sum

testing/test_grid.m

@@ -0,0 +1,21 @@
+format compact
+format short g
+
+clf
+f = @clenshaw_curtis_legendre_rule;
+x = smolyak_grid(2, 5, f);
+size(x)
+plot(x(1,:), x(2,:), 'k.');
+axis square
+
+f = @gauss_hermite_rule;
+smolyak_grid(4,7,f);
+
+
+% for d=1:5
+%     fprintf('\n\n');
+%     underline(sprintf('d=%d', d));
+%     for s=1:7
+%         fprintf('s=%d  -  %4d  %4d  -  %4d %4d\n', [s, s^d, size(full_tensor_grid(d, s, f),2),  multiindex_size(d,s-1) size(smolyak_grid(d, s, f),2)]);
+%     end
+% end

testing/test_nd_exact.m

@@ -0,0 +1,39 @@
+function test_nd_exact
+
+pt = 13;
+for p=1:4
+    subplot(3,4,p)
+    test(13, p, @gauss_legendre_rule, 'smolyak')
+    subplot(3,4,p+4)
+    test(13, p, @gauss_legendre_rule, 'full_tensor')
+    subplot(3,4,p+8)
+    test(13, p, @clenshaw_curtis_legendre_rule, 'smolyak')
+end
+
+function test(pt, p, rule_func, grid)
+I = multiindex(2, pt, 'full', true);
+q = integrate_nd({@kernel, {I,[]}, {2,3}}, rule_func, 2, p, 'grid', grid);
+q_ex = prod(2./(I+1),2);
+ok = (abs(q-q_ex)<1e-10);
+plot(I(ok,1), I(ok,2), 'b.'); hold on;
+plot(I(~ok,1), I(~ok,2), 'r.'); hold off;
+xlim( [min(I(:,1))-0.5, max(I(:,1))+0.5] );
+ylim( [min(I(:,2))-0.5, max(I(:,2))+0.5] );
+axis square
+
+%[multiindex_order(I),q,q_ex,q-q_ex]
+%keyboard
+
+function y = kernel(x, I, weight_func)
+
+if ~isempty(weight_func)
+    w = funcall(weight_func, x);
+else
+    w = ones(size(x,2),1);
+end
+
+y = gpc_eval_basis({'M', I}, 0.5*(x+1)) * diag(w);
+
+
+
+

testing/test_sample_beta_gpc.m

@@ -0,0 +1,35 @@
+%%
+clear
+dist_func={@beta_stdnor, {4, 2}, {2, 3}};
+[a_i_alpha, I] = pce_expand_1d(dist_func, 4)
+V = {'H', I}
+
+N=100000;
+xi=gpc_sample(V, N);
+y=gpc_evaluate(a_i_alpha, V, xi);
+kde(y)
+
+
+%%
+clear
+
+dist = 'beta';
+params = {4, 2};
+params = {0.5, 0.7};
+[shift, scale]=gendist_fix_moments(dist, params, 3.2, 0.24);
+x=linspace(-1,5);
+plot(x,gendist_pdf(x, dist, params, shift, scale))
+hold all
+dist_func={@gendist_stdnor, {dist, params, shift, scale}, {2,3,4,5}}
+
+for p=1:6
+    [a_i_alpha, I] = pce_expand_1d(dist_func, p)
+    V = {'H', I}
+
+    N=100000;
+    xi=gpc_sample(V, N);
+    y=gpc_evaluate(a_i_alpha, V, xi);
+    %kde(y,100)
+    empirical_density(y);
+    legend('pdf', '1', '2', '3', '4', '5', '6' )
+end

testing/testing_mc1.m

@@ -0,0 +1,34 @@
+state = electrical_network_init();
+
+V = {'pp', multiindex(2, 0)};
+% gpcspace_create('uniform', 2, 0)
+
+N = 1000;
+xi = gpc_sample(V, N);
+plot(xi(1,:), xi(2,:), 'x')
+
+n = state.num_vars;
+u = zeros(n, N);
+for i=1:N
+   p = xi(:, i);
+    u(:,i) = electrical_network_solve(state, p);
+end
+
+fprintf('u_%d = %5.3g +- %5.3g\n', [1:n; mean(u, 2)'; std(u, [], 2)']);
+
+
+for i=1:n
+    for j=1:n
+        if i==j
+            continue
+        end
+        subplot(n, n, i+(j-1)*n);
+        plot(u(i,:), u(j,:), 'x')
+        xlim([.22,.39]);
+        ylim([.22,.39]);
+    end
+end
+
+
+
+

testing/testing_mc2.m

@@ -0,0 +1,17 @@
+state = electrical_network_init();
+
+V = {'pp', multiindex(2, 0)};
+
+N = 1000;
+n = state.num_vars;
+u_mean = [];
+u_var = [];
+for i=1:N
+    % p = gpc_sample(V, 'rand', @rand);
+    p = gpc_sample(V);
+    u = electrical_network_solve(state, p);
+    [u_mean, u_var] = mean_var_update(i, u, u_mean, u_var);
+end
+
+fprintf('u_%d = %5.3g +- %5.3g\n', [1:n; u_mean'; sqrt(u_var)']);
+

testing/testing_sparse2.m

@@ -0,0 +1,28 @@
+d = 3;
+l = 4;
+
+
+[x, w]  = full_tensor_grid(d, l, @clenshaw_curtis_legendre_rule);
+subplot(2,2,1)
+plot3(x(1,:), x(2,:), x(3,:), '.')
+axis square
+
+
+[x, w]  = smolyak_grid(d, l, @clenshaw_curtis_legendre_rule);
+subplot(2,2,2)
+plot3(x(1,:), x(2,:), x(3,:), '.')
+axis square
+
+l = 2^l+1;
+l = 9;
+
+[x, w]  = full_tensor_grid(d, l, @gauss_legendre_rule);
+subplot(2,2,3)
+plot3(x(1,:), x(2,:), x(3,:), '.')
+axis square
+
+
+[x, w]  = smolyak_grid(d, l, @gauss_legendre_rule);
+subplot(2,2,4)
+plot3(x(1,:), x(2,:), x(3,:), '.')
+axis square

testing/testing_sparse_int.m

@@ -0,0 +1,26 @@
+state = electrical_network_init();
+
+V = {'p', multiindex(2, 0)};
+
+
+xw = gpc_integrate([], V, 6, 'grid', 'full_tensor');
+%xw = gpc_integrate([], V, 6, 'grid', 'smolyak');
+x = xw{1};
+w = xw{2};
+plot(x(1,:), x(2,:), 'x')
+axis square
+
+n = state.num_vars;
+u_mean = zeros(n, 1);
+u_m2 = zeros(n, 1);
+
+for i=1:length(w)
+    p = x(:,i);
+    u = electrical_network_solve(state, p);
+
+    u_mean = u_mean + w(i) * u;
+    u_m2   = u_m2   + w(i) * u.^2;
+end
+u_var = u_m2 - u_mean.^2;
+
+fprintf('u_%d = %5.3g +- %5.3g\n', [1:n; u_mean'; sqrt(u_var)']);