active_subspaces/test_base2.py at master · aerialhedgehog/active_subspaces · 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
89
90
91
92
93
import numpy as np
import active_subspaces as asm

def quad_fun(x):
    A = np.array([[ 0.2406659045776698, -0.3159904335007421, -0.1746908591702878],
       [-0.3159904335007421,  0.5532215729009683,  0.3777995408101305],
       [-0.1746908591702878,  0.3777995408101305,  0.3161125225213613]])
    f = 0.5*np.dot(x,np.dot(A,x.T))
    return f

def quad_grad(x):
    A = np.array([[ 0.2406659045776698, -0.3159904335007421, -0.1746908591702878],
       [-0.3159904335007421,  0.5532215729009683,  0.3777995408101305],
       [-0.1746908591702878,  0.3777995408101305,  0.3161125225213613]])
    df = np.dot(A,x.T)
    return df

if __name__ == '__main__':
    print 'Here we go'
    m = 3
    bflag = True

    # initialize the model
    model = asm.ActiveSubspaceModel(bflag=bflag)

    # build the model from an interface
    model.build_from_interface(m,quad_fun)
    model.build_from_interface(m,quad_fun,dfun=quad_grad)

    """
    X = np.random.uniform(-1.0,1.0,size=(100,3))
    f = np.zeros((100,1))
    df = np.zeros((100,3))
    for i in range(100):
        f[i] = quad_fun(X[i,:].T)
        v = quad_grad(X[i,:].T)
        df[i,:] = v.T
    np.savez('quad_data',X=X,f=f,df=df)
    """

    # build the model from data
    qd = np.load('active_subspaces/tests/data/quad_data.npz')
    X,f,df = qd['X'],qd['f'],qd['df']
    model.build_from_data(X,f)
    model.build_from_data(X,f,df=df)

    # make some plots
    #model.diagnostics()

    # trying 2d
    model.subspaces.partition(n=1)

    # just for tests -- these ran without error
    model.set_domain()
    model.set_response_surface()

    """
    XX = np.random.uniform(-1.0,1.0,size=(10,3))
    np.savez('quad_predict',XX=XX)

    """
    # check response surface predictions
    qd = np.load('active_subspaces/tests/data/quad_predict.npz')
    XX = qd['XX']
    ff,dff,vv = model.predict(XX)
    ff,dff,vv = model.predict(XX,compvar=True)
    ff,dff,vv = model.predict(XX,compgrad=True)
    ff,dff,vv = model.predict(XX,compgrad=True,compvar=True)
    print 'Predictions'
    print ff
    print 'Gradients'
    print dff
    print 'Variances'
    print vv

    # check statistics
    mu = model.mean(N=20)
    print 'Mean'
    print mu
    sig = model.variance(N=20)
    print 'Variance'
    print sig
    np.random.seed(seed=1000)
    p = model.probability(0.0,0.05)
    print 'Probability'
    print p

    # trying optimization
    fmin,xmin = model.minimum()
    print 'Minimum'
    print fmin
    print 'Argmin'
    print xmin