UQ-Project/eval2Dmodel.m at master · TimWildey/UQ-Project · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
function varargout = eval2Dmodel(p)

N = size(p,1);

model = 'nlinv';

switch model
    case 'rosenbrock'
        a = 1;
        b = 100;
        qvals = (a-p(:,1)).^2 + b*(p(:,2)-p(:,1).^2).^2;
    case 'peaks'
        p1 = gaussian(2.0, p(:,1), p(:,2), 0.25, 0.75, 0.15,0.15);
        p2 = gaussian(3.0, p(:,1), p(:,2), 0.75, 0.75, 0.2,0.2);
        p3 = gaussian(2.5, p(:,1), p(:,2), 0.33, 0.33, 0.1,0.1);
        p4 = gaussian(-1.0, p(:,1), p(:,2), 0.8, 0.4, 0.1,0.2);

        qvals = p1+p2+p3+p4;
    case 'linear1'

        A = [ 1 2 ; -1 1];
        qvals = A*p';
        qvals = qvals(1,:)';

    case 'poly'

%         qvals = (1*p(:,1) - 2*p(:,1).^2 - 3*p(:,1).^3).*(2*p(:,2) + p(:,2).^2 - 2*p(:,2).^4);
        qvals = (1*p(:,1) + sin(3*pi*p(:,1))).*(2*p(:,2) + p(:,2).^2 - 2*p(:,2).^4);

    case 'linear'

        qvals = 0.0*p(:,1) - 2.0*p(:,2);

    case 'hump'

        qvals = tanh(40*(p(:,1)-1/2));

    case 'nlinv'
        qvals = zeros(N,2);
        for k = 1:N
            x = [1;1];

            resid = [1;1] - [p(k,1)*x(1)^2 + x(2)^2;x(1)^2 - p(k,2)*x(2)^2];
            J = [2*p(k,1)*x(1) 2*x(2);2*x(1) -2*p(k,2)*x(2)];

            iter = 0;
            while norm(resid) > 1e-10 && iter < 20

                iter = iter + 1;
                dx = J\resid;
                x = x+dx;
                resid = [1;1] - [p(k,1)*x(1)^2 + x(2)^2;x(1)^2 - p(k,2)*x(2)^2];
                J = [2*p(k,1)*x(1) 2*x(2);2*x(1) -2*p(k,2)*x(2)];

            end

            qvals(k,1) = x(1);
            qvals(k,2) = x(2);
        end

    case 'discont'
        x1 = 0+1*p(:,1);
        x2 = 0+1*p(:,2);

        f1 = exp(-x1.^2-x2.^2) - x1.^3-x2.^3;
        f2 = 1+f1+1/8*x2.^2;
        f = zeros(N,1);

        for j = 1:N
            if ((3*x1(j)+2*x2(j))>=0) && ((-x1(j)+0.3*x2(j))<0)
                f(j) = f1(j)-2;
            elseif  ((3*x1(j)+2*x2(j))>=0) && ((-x1(j)+0.3*x2(j))>=0)
                f(j) = 2*f2(j);
            elseif (x1(j)+1)^2+(x2(j)+1)^2 <= 0.95^2
                f(j) = 2*f1(j)+4;
            else
                f(j) = f1(j);
            end
        end

        qvals = f;
end

if nargout == 1
    varargout{1} = qvals;
    return;
end