| 1 | """Least squares smooting and interpolation. |
|---|
| 2 | |
|---|
| 3 | measure the speed of least squares. |
|---|
| 4 | |
|---|
| 5 | Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou |
|---|
| 6 | Geoscience Australia, 2004. |
|---|
| 7 | """ |
|---|
| 8 | |
|---|
| 9 | |
|---|
| 10 | import unittest |
|---|
| 11 | from math import sqrt |
|---|
| 12 | import time |
|---|
| 13 | |
|---|
| 14 | |
|---|
| 15 | from pyvolution.least_squares import Interpolation |
|---|
| 16 | from Numeric import allclose, array, transpose |
|---|
| 17 | |
|---|
| 18 | from coordinate_transforms.geo_reference import Geo_reference |
|---|
| 19 | |
|---|
| 20 | def distance(x, y): |
|---|
| 21 | return sqrt( sum( (array(x)-array(y))**2 )) |
|---|
| 22 | |
|---|
| 23 | def linear_function(point): |
|---|
| 24 | point = array(point) |
|---|
| 25 | return point[:,0]+point[:,1] |
|---|
| 26 | |
|---|
| 27 | |
|---|
| 28 | class Test_Least_Squares(unittest.TestCase): |
|---|
| 29 | |
|---|
| 30 | def setUp(self): |
|---|
| 31 | pass |
|---|
| 32 | |
|---|
| 33 | def tearDown(self): |
|---|
| 34 | pass |
|---|
| 35 | |
|---|
| 36 | def ztest_expand_search2(self): |
|---|
| 37 | from random import seed, random |
|---|
| 38 | from pmesh.mesh import Mesh |
|---|
| 39 | |
|---|
| 40 | seed(2) |
|---|
| 41 | m = Mesh() |
|---|
| 42 | verts_per_sector = 100 #100 |
|---|
| 43 | num_of_points = 300 |
|---|
| 44 | sec_side = 9.0 |
|---|
| 45 | |
|---|
| 46 | vert_atts = [] |
|---|
| 47 | for vert in range(verts_per_sector): |
|---|
| 48 | m.addUserVertex(random()*sec_side, random()*sec_side) |
|---|
| 49 | vert_atts.append(random()) |
|---|
| 50 | for vert in range(verts_per_sector): |
|---|
| 51 | m.addUserVertex(90.0+random()*sec_side, 90.0+random()*sec_side) |
|---|
| 52 | vert_atts.append(random()) |
|---|
| 53 | for vert in range(verts_per_sector): |
|---|
| 54 | m.addUserVertex(random()*sec_side, 90.0+random()*sec_side) |
|---|
| 55 | vert_atts.append(random()) |
|---|
| 56 | for vert in range(verts_per_sector): |
|---|
| 57 | m.addUserVertex(90.0+random()*sec_side, random()*sec_side) |
|---|
| 58 | vert_atts.append(random()) |
|---|
| 59 | #for vert in range(verts_per_sector): |
|---|
| 60 | # m.addUserVertex(40.0+random()*sec_side, 40.0+random()*sec_side) |
|---|
| 61 | # vert_atts.append(random()) |
|---|
| 62 | |
|---|
| 63 | m.autoSegment(alpha = 100 ) |
|---|
| 64 | m.generateMesh(mode = "Q") |
|---|
| 65 | m.export_mesh_file("aaaa.tsh") |
|---|
| 66 | mesh_dict = m.Mesh2IOTriangulationDict() |
|---|
| 67 | |
|---|
| 68 | points = [] |
|---|
| 69 | point_atts = [] |
|---|
| 70 | for point in range(num_of_points): |
|---|
| 71 | points.append([random()*100, random()*100]) |
|---|
| 72 | point_atts.append(10.0) |
|---|
| 73 | m.autoSegment(alpha = 100 ) |
|---|
| 74 | |
|---|
| 75 | #print 'points', points |
|---|
| 76 | t0 = time.time() |
|---|
| 77 | |
|---|
| 78 | interp = Interpolation(mesh_dict['vertices'], |
|---|
| 79 | mesh_dict['triangles'], points, |
|---|
| 80 | alpha=0.2, expand_search=True, |
|---|
| 81 | verbose = False, |
|---|
| 82 | max_points_per_cell = 4) |
|---|
| 83 | calc = interp.fit_points(point_atts ) |
|---|
| 84 | #print "interp.expanded_quad_searches", interp.expanded_quad_searches |
|---|
| 85 | #print "calc", calc |
|---|
| 86 | print 'That took %.2f seconds' %(time.time()-t0) |
|---|
| 87 | |
|---|
| 88 | |
|---|
| 89 | def test_expand_search3(self): |
|---|
| 90 | from random import seed, random |
|---|
| 91 | from pmesh.mesh import Mesh |
|---|
| 92 | |
|---|
| 93 | num_of_points = 20000 |
|---|
| 94 | |
|---|
| 95 | seed(2) |
|---|
| 96 | m = Mesh() |
|---|
| 97 | m.addUserVertex(0,0) |
|---|
| 98 | m.addUserVertex(100,0) |
|---|
| 99 | m.addUserVertex(0,100) |
|---|
| 100 | m.addUserVertex(100,100) |
|---|
| 101 | |
|---|
| 102 | m.autoSegment(alpha = 100 ) |
|---|
| 103 | |
|---|
| 104 | dict = {} |
|---|
| 105 | dict['points'] = [[10,10],[90,20]] |
|---|
| 106 | dict['segments'] = [[0,1]] |
|---|
| 107 | dict['segment_tags'] = ['wall1'] |
|---|
| 108 | m.addVertsSegs(dict) |
|---|
| 109 | |
|---|
| 110 | dict = {} |
|---|
| 111 | dict['points'] = [[10,90],[40,20]] |
|---|
| 112 | dict['segments'] = [[0,1]] |
|---|
| 113 | dict['segment_tags'] = ['wall2'] |
|---|
| 114 | m.addVertsSegs(dict) |
|---|
| 115 | |
|---|
| 116 | dict = {} |
|---|
| 117 | dict['points'] = [[20,90],[60,60]] |
|---|
| 118 | dict['segments'] = [[0,1]] |
|---|
| 119 | dict['segment_tags'] = ['wall3'] |
|---|
| 120 | m.addVertsSegs(dict) |
|---|
| 121 | |
|---|
| 122 | dict = {} |
|---|
| 123 | dict['points'] = [[60,20],[90,90]] |
|---|
| 124 | dict['segments'] = [[0,1]] |
|---|
| 125 | dict['segment_tags'] = ['wall4'] |
|---|
| 126 | m.addVertsSegs(dict) |
|---|
| 127 | |
|---|
| 128 | m.addVertsSegs(dict) |
|---|
| 129 | m.generateMesh(mode = "Q", maxArea = "1000") |
|---|
| 130 | m.export_mesh_file("aaaa.tsh") |
|---|
| 131 | mesh_dict = m.Mesh2IOTriangulationDict() |
|---|
| 132 | |
|---|
| 133 | points = [] |
|---|
| 134 | point_atts = [] |
|---|
| 135 | for point in range(num_of_points): |
|---|
| 136 | points.append([random()*100, random()*100]) |
|---|
| 137 | point_atts.append(10.0) |
|---|
| 138 | #m.autoSegment(alpha = 100 ) |
|---|
| 139 | |
|---|
| 140 | #print 'points', points |
|---|
| 141 | t0 = time.time() |
|---|
| 142 | interp = Interpolation(mesh_dict['vertices'], |
|---|
| 143 | mesh_dict['triangles'], points, |
|---|
| 144 | alpha=0.2, expand_search=True, |
|---|
| 145 | verbose = True, |
|---|
| 146 | max_points_per_cell = 4) |
|---|
| 147 | calc = interp.fit_points(point_atts ) |
|---|
| 148 | print "interp.expanded_quad_searches", interp.expanded_quad_searches |
|---|
| 149 | #print "calc", calc |
|---|
| 150 | print 'That took %.2f seconds' %(time.time()-t0) |
|---|
| 151 | |
|---|
| 152 | #------------------------------------------------------------- |
|---|
| 153 | if __name__ == "__main__": |
|---|
| 154 | suite = unittest.makeSuite(Test_Least_Squares,'test') |
|---|
| 155 | runner = unittest.TextTestRunner(verbosity=1) |
|---|
| 156 | runner.run(suite) |
|---|
