Merge pull request #67 from SCIInstitute/uncertainty

Uncertainty

PythonLibrary/ICP/icp.py

@@ -4,7 +4,7 @@
 
 @author: andrewcorbato
 """
-def icp(model,data,*args):
+def icp(model,data,*arg):
 #   ICP Iterative Closest Point Algorithm. Takes use of
 #   Delaunay tesselation of points in model.
 #
@@ -110,66 +110,83 @@ def icp(model,data,*args):
 #   To clear old global tesselation variables run: "clear global Tes ir jc" (tes_flag=1)
 #   or run: "clear global Tes Tesind Tesver" (tes_flag=2) in Command Window.
 #
-#   m-file can be downloaded for free at
+#   based on an m-file can be downloaded for free at
 #   http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=12627&objectType=FILE
 #
-#   icp version 1.4
 #
-#   written by Per Bergstrm 2007-03-07
 
     import numpy as np
     import scipy as sp
     from numpy import matlib as ml
     from numpy import matrix as mx
+    import sys
 
-    if model == False and data == False:
-        import sys
-        sys.exit('ERROR: Only one input. There must be at least two inputs.')
-    elif model == False:
-        import sys
-        sys.exit('ERROR: Model cannot be an empty matrix.')
+    if model.any() == False and data.any() == False:
+        raise ValueError('ERROR: Only one input. There must be at least two inputs.')
+    elif model.any() == False:
+        raise ValueError('ERROR: Model cannot be an empty matrix.')
     global MODEL, DATA, TR, TT, TS
 
     dims_model = np.shape(model)
     if dims_model[1]<dims_model[0]:
-        global MODEL
         MODEL = np.transpose(model)
     else:
-        global MODEL
         MODEL = model
 
     dims_data = np.shape(data)
     if dims_data[1]<dims_data[0]:
         data = np.transpose(data)
         dims_data = np.shape(data)
-        global DATA        
         DATA = data
     else:
-        global DATA
         DATA = data
 
     if dims_data[0] is not dims_model[0]:
-        import sys
-        sys.exit('ERROR: DATA and MODEL cannot have different dimensions.')
-# default options       
-    if len(args) == 0:
-        options = {'max_iter': 104, 'min_iter': 4,'fitting': 2,'thres': 1e-05, \
+        raise ValueError('ERROR: DATA and MODEL cannot have different dimensions.')
+# default options
+    if not len(arg)==0:
+        opt = arg[0]
+    else:
+        opt = {'max_iter': 104, 'min_iter': 4,'fitting': [2],'thres': 1e-05, \
         'init_flag': 1,'tes_flag': 1, 'delta': 10,'mode': 'rigid', \
         'refpnt': np.array([]),}
-    elif len(args) != 9:
-        import sys
-        sys.exit('ERROR: Options dictionary must have 9 entries.')
+
+    print(opt)
+
+    if not isinstance(opt, dict):
+        raise ValueError("opt must be dictionary of options")
+    if not 'max_iter' in opt:
+        opt['max_iter'] = 104;
+    if not 'min_iter' in opt:
+        opt['min_iter'] = 4;
+    if not 'fitting' in opt:
+        opt['fitting'] = [2];
+    if not 'thres' in opt:
+        opt['thres'] = 1e-5;
+    if not 'delta' in opt:
+        opt['delta'] = 10;
+    if not 'mode' in opt:
+        opt['mode'] = 'rigid';
+    if not 'refpnt' in opt:
+        opt['refpnt'] = np.array([]);
+    if not 'init_flag' in opt:
+        opt['init_flag'] = 1;
+    if not 'tes_flag' in opt:
+        opt['tes_flag'] = 1;
+
+
+
 
 # move options out of dictionary for ease of use
-    max_iter = args['max_iter']
-    min_iter = args['min_iter']
-    fitting = args['fitting']
-    thres = args['thres']
-    init_flag = args['init_flag']
-    tes_flag = args['tes_flag']
-    delta = args['delta']
-    mode = args['mode']
-    refpnt = args['refpnt']
+    max_iter = opt['max_iter']
+    min_iter = opt['min_iter']
+    fitting = opt['fitting']
+    thres = opt['thres']
+    init_flag = opt['init_flag']
+    tes_flag = opt['tes_flag']
+    delta = opt['delta']
+    mode = opt['mode']
+    refpnt = opt['refpnt']
 
 # input error checks
     if (tes_flag != 0 and tes_flag != 1 and tes_flag != 2 and tes_flag != 3):
@@ -178,7 +195,7 @@ def icp(model,data,*args):
 must be 0-3.')    
     if dims_model[1] == dims_model[0] and tes_flag is not 0:
         import sys
-        sys.exit('ERROR: This icp method requires the number of model points \
+        raise ValueError('ERROR: This icp method requires the number of model points \
 to be greater than the dimension')
     if max_iter < min_iter:
         max_iter = min_iter
@@ -309,9 +326,10 @@ def icp_init(init_flag,fitting):
         DATA = DATA + ml.repmat(TT,1,N) # apply transformation
     else:
         import sys
-        sys.exit('ERROR: Wrong init_flag')
+        raise ValueError('ERROR: Wrong init_flag')
 ###############################################################################    
 def icp_struct(tes_flag):
+    global Tes, ir, jc, Tesind, Tesver
     if tes_flag != 3:
         if tes_flag == 0:
             global ir
@@ -321,15 +339,15 @@ def icp_struct(tes_flag):
                 dims_Tesind = np.shape(Tesind)
                 if dims_Tesind[0] == 0:
                     import sys
-                    sys.exit('ERROR: No tesselation system exists')
+                    raise ValueError('ERROR: No tesselation system exists')
                 else:
                     tes_flag = 2
             else:
                 tes_flag = 1
     elif tes_flag == 3:
         return tes_flag
     else:
-        global MODEL, Tes
+
         [m,n] = np.shape(MODEL)
         if m == 1:
             ind1 = np.argsort(MODEL)
@@ -364,7 +382,6 @@ def icp_struct(tes_flag):
     Tes = np.sort(np.sort(Tes,axis=1),axis=1)
     [mT,nT] = np.shape(Tes)
     if tes_flag == 1:
-        global ir, jc
         num = np.zeros((1,n))
         for jj in range(0,mT,1):
             for kk in range(0,nT,1):
@@ -381,7 +398,6 @@ def icp_struct(tes_flag):
                 ir[ind[Tes[i,j]]] = i
                 ind[Tes[i,j]] = ind[Tes[i,j]]+1
     else: # tes_flag ==2
-        global Tesind, Tesver
         Tesind = np.zeros(mT,nT)
         Tesver = np.zeros(mT,nT)
         couvec = np.zeros(mT,1)
@@ -458,8 +474,8 @@ def icp_closest_start(tes_flag,fitting):
 # ICP_CLOSEST_START finds indexes of closest MODEL points for each point in DATA.
 # The _start version allocates memory for iclosest and finds the closest MODEL points
 # to the DATA points
+    global MODEL, DATA, Tes, ir, jc, iclosest
     if tes_flag == 3:
-        global MODEL, DATA, iclosest
         dims_MODEL = np.shape(MODEL)
         dims_DATA = np.shape(DATA)
         mm = dims_MODEL[1]
@@ -477,25 +493,24 @@ def icp_closest_start(tes_flag,fitting):
         ERROR += err(dist,fitting,ID)
 
     elif tes_flag == 1:
-        global MODEL, DATA, Tes, ir, jc, iclosest
         dims_DATA = np.shape(DATA)
         md = dims_DATA[1]
         dims_MODEL = np.shape(MODEL)
         iclosest = np.zeros((1,md))
         mid = np.round(md/2)
-        iclosest(mid) = np.round((dims_MODEL[1]/2))
+        iclosest[mid] = np.round((dims_MODEL[1]/2))
         bol = 1
         while bol:
             bol = np.logical_not(bol)
-            distc = np.linalg.norm((DATA[:,mid]-MODEL[:,iclosest(mid)]))
+            distc = np.linalg.norm((DATA[:,mid]-MODEL[:,iclosest[mid]]))
             distcc = 2*distc
             for i in range(ir[jc[iclosest[mid]]],ir[jc[iclosest[mid]+1]-1],1):
                 for ind in Tes[i,:]:
                     distcc = np.linalg.norm((DATA[:,mid]-MODEL[:,ind]))
                     if distcc<distc:
                         distc = distcc
                         bol = np.logical_not(bol)
-                        iclosest(mid) = ind
+                        iclosest[mid] = ind
                         break
                 if bol:
                     break
@@ -514,7 +529,7 @@ def icp_closest_start(tes_flag,fitting):
                         if distcc<distc:
                             distc = distcc
                             bol = np.logical_not(bol)
-                            iclosest(mid) = ind
+                            iclosest[mid] = ind
                             break
                     if bol:
                         break
@@ -533,10 +548,9 @@ def icp_closest_start(tes_flag,fitting):
                         if distcc<distc:
                             distc = distcc
                             bol = np.logical_not(bol)
-                            iclosest(mid) = ind
+                            iclosest[mid] = ind
                             break
     else: # tes_flag == 2
-        global MODEL, DATA, Tes, Tesind, Tesver, icTesind, iclosest
         dims_DATA = np.shape(DATA)
         md = dims_DATA[1]
         iclosest = np.zeros((1,md))
@@ -691,8 +705,8 @@ def icp_transformation(fitting,refpnt,delta,optmode):
                     b = np.concatenate(0,DATA[0,n],0,0,DATA[1,n],0,0,DATA[2,n],0,0,1,0, axis=1)
                     c = np.concatenate(0,0,DATA[0,n],0,0,DATA[1,n],0,0,DATA[2,n],0,0,1, axis=1)
                     P[ind,:] = np.concatenate(a,b,c, axis=0)
-                    q[ind] = MODEL[0:2,iclosest[p_ind[n]]
-                    ind += 3
+                    q[ind] = MODEL[0:2,iclosest[p_ind[n]]]
+                    ind = 3
     if mode == 'affine':
         theta = q/P
         a_r = np.concatenate(theta[0],theta[3],theta[6],axis=0)
@@ -734,7 +748,7 @@ def err(dist,fitting,ind):
         ERR = dist**2
     else:
         ERR = 0
-        print(WARNING: Unknown fitting value.)
+        print('WARNING: Unknown fitting value.')
     return ERR
 ###############################################################################        
 def weightfcn(distances):
@@ -745,4 +759,4 @@ def weightfcn(distances):
     else:
         weights = max_distances+min_distances-distances
     weights = weights/np.sum(weights)
-    return weights
+    return weights

PythonLibrary/uncertainty/demo.py

@@ -0,0 +1,125 @@
+# Shows some PCE analysis from samples:
+# - generation of samples
+# - construction of PCE coefficients
+# - sensitivity analysis of PCE coefficients
+# - comparison of PCE predictor versus truth
+
+import itertools
+import numpy as np
+from matplotlib import pyplot as plt
+
+import indexing, wafp, sensitivity
+from test_functions import genz_oscillatory
+from opolynd import opolynd_eval
+from recurrence import jacobi_recurrence
+
+# construction of our "expensive model". Here just a collection of
+# not-so-complicated test functions.
+Nx = 100                    # Number of "spatial" grid points (here, space is 1D)
+x = np.linspace(0, 1, Nx)   # spatial grid
+print(x)
+
+def expensive_model(zz,xx):
+    """
+    Evaluates a toy function as if it were an expensive model.
+
+    x is a "spatial" grid (numpy array)
+    z is a single d-variate parameter/random variable instance.
+
+    Returns an array shaped like xx.
+    """
+
+    r = 2. # Fixed input into genz oscillatory function
+
+    d = zz.size
+
+    f = np.zeros(xx.shape)
+    for ind,xval in np.ndenumerate(xx):
+        f[ind] = genz_oscillatory(xval*np.ones(d), r)(zz)
+
+    return f
+
+######## Step 1: generate samples
+
+d = 2                       # dimension of random/parameter space
+k = 7                       # polynomial degree (parameter space)
+poly_space = indexing.total_degree_indices
+
+lambdas = poly_space(d, k)
+N = lambdas.shape[0]
+
+M = N + 10  # Number of samples (expensive model runs) to perform. Must
+            # be at least N, but the +10 is arbitrary.
+            # Of course, large M  ===> better.
+
+candidate_mesh_size = int(2e2)
+
+print("Generating parameter mesh...")
+z = wafp.legendre_wafp(lambdas, M=candidate_mesh_size)
+z = wafp.legendre_wafp_enrichment(z, lambdas, M-N)
+
+# The samples are the array z, each row is a d-dimensional sample on the
+# hypercube [-1,1]^d. The particular way these points are generated is
+# random, so you'll get a different set of z each time this is fun, but
+# the "quality" of the grid has relatively low variance from run to run.
+# This takes a while, but these points may be stored for future use.
+
+######## Step 2: run "expensive" model
+
+u = np.zeros([Nx, M])
+
+print("Evaluating model on mesh...")
+for ind,zval in enumerate(z):
+    u[:,ind] = expensive_model(zval, x)
+
+print(u)
+
+# Each column of u is a model run for a fixed parameter value
+
+######## Step 3: compute PCE coefficients
+print("Assembling PCE coefficients...")
+ab = jacobi_recurrence(lambdas.max()+1, alpha=0., beta=0., probability=True)
+V = opolynd_eval(z, lambdas, ab)
+weights = np.sqrt(float(N)/float(M)) / np.sqrt(np.sum(V**2,axis=1))
+
+# The PCE coefficients are computed as a weighted discrete least-squares
+# estimator with the weights above.
+coeffs = np.linalg.lstsq( (V.T*weights).T, (u*weights).T)[0].T
+
+# Each row of coeffs contains PCE coefficients for a single x gridpoint.
+# Each column of coeffs contains a particular PCE coefficient for all
+# values of x.
+
+######## Step 4: whatever postprocessing you want
+print("Processing PCE coefficients...")
+# Compute total sensitivities
+total_sensitivities = sensitivity.pce_total_sensitivity(coeffs.T, lambdas, list(range(d)))
+
+# Compute global sensitivities
+# Compute main-effect and main-interaction sensitivities
+Js = [[j] for j in range(d)]
+[Js.append(comb) for comb in itertools.combinations(list(range(d)), 2)]
+
+global_sensitivities = sensitivity.pce_global_sensitivity(coeffs.T, lambdas, Js)
+
+# Compare surrogate discrepancy at validation points
+M_validation = 100
+z_validation = np.random.uniform(0, 1, [100, d])
+
+u_truth = np.zeros([Nx, M_validation])
+for ind, zval in enumerate(z_validation):
+    u_truth[:,ind] = expensive_model(zval, x)
+
+u_pce = np.zeros([Nx, M_validation])
+V = opolynd_eval(z_validation, lambdas, ab)
+u_pce = np.dot(V, coeffs.T).T
+
+# Compute l2- and max-errors for each grid point:
+l2_error = np.sqrt(np.sum((u_truth - u_pce)**2, axis=1))/np.sqrt(N)
+linf_error = np.max(np.abs(u_truth - u_pce), axis=1)
+
+print("L2 error on validation mesh: {0:1.3e}\nMaximum error on validation mesh: {1:1.3e}".format(np.linalg.norm(l2_error)/np.sqrt(Nx), np.max(linf_error.flatten())))
+
+if d == 2:
+    plt.plot(z[:,0], z[:,1], 'r.')
+    plt.show()