1 | /* |
---|
2 | * ---------------------------------------------------------------------- |
---|
3 | * Rappture::DX |
---|
4 | * |
---|
5 | * Rappture DX object for file reading and interacting |
---|
6 | * with libDX and friends. |
---|
7 | * |
---|
8 | * ====================================================================== |
---|
9 | * AUTHOR: Derrick S. Kearney, Purdue University |
---|
10 | * Copyright (c) 2005-2008 Purdue Research Foundation |
---|
11 | * |
---|
12 | * See the file "license.terms" for information on usage and |
---|
13 | * redistribution of this file, and for a DISCLAIMER OF ALL WARRANTIES. |
---|
14 | * ====================================================================== |
---|
15 | */ |
---|
16 | #include "RpDX.h" |
---|
17 | #include <math.h> |
---|
18 | #include <stdio.h> |
---|
19 | #include <stdlib.h> |
---|
20 | #include <float.h> |
---|
21 | |
---|
22 | using namespace Rappture; |
---|
23 | |
---|
24 | DX::DX(Outcome &result, const char* filename) : |
---|
25 | _dataMin(FLT_MAX), |
---|
26 | _dataMax(-FLT_MAX), |
---|
27 | _nzero_min(FLT_MAX), |
---|
28 | _numAxis(0), |
---|
29 | _axisLen(NULL), |
---|
30 | _data(NULL), |
---|
31 | _n(0), |
---|
32 | _rank(0), |
---|
33 | _shape(0), |
---|
34 | _positions(NULL), |
---|
35 | _delta(NULL), |
---|
36 | _max(NULL), |
---|
37 | _origin(NULL) |
---|
38 | { |
---|
39 | Array dxpos; |
---|
40 | Array dxdata; |
---|
41 | Array dxgrid; |
---|
42 | // category and type are probably not needed |
---|
43 | // we keep them around in case they hold useful info for the future |
---|
44 | // we could replace them with NULL inthe DXGetArrayInfo fxn call |
---|
45 | Category category; |
---|
46 | Type type; |
---|
47 | |
---|
48 | if (filename == NULL) { |
---|
49 | result.addError("filename is NULL"); |
---|
50 | return; |
---|
51 | } |
---|
52 | // open the file with libdx |
---|
53 | fprintf(stdout, "Calling DXImportDX(%s)\n", filename); |
---|
54 | fflush(stdout); |
---|
55 | DXenable_locks(0); |
---|
56 | _dxobj = DXImportDX((char*)filename,NULL,NULL,NULL,NULL); |
---|
57 | if (_dxobj == NULL) { |
---|
58 | result.addError("can't import DX file \"%s\"", filename); |
---|
59 | return; |
---|
60 | } |
---|
61 | |
---|
62 | // parse out the positions array |
---|
63 | |
---|
64 | // FIXME: nanowire will need a different way to parse out the positions |
---|
65 | // array since it uses a productarray to store its positions. |
---|
66 | // Possibly use DXGetProductArray(). |
---|
67 | dxpos = (Array) DXGetComponentValue((Field) _dxobj, (char *)"positions"); |
---|
68 | if (dxpos == NULL) { |
---|
69 | result.addError("can't get component value of \"positions\""); |
---|
70 | return; |
---|
71 | } |
---|
72 | DXGetArrayInfo(dxpos, &_n, &type, &category, &_rank, &_shape); |
---|
73 | |
---|
74 | fprintf(stdout, "_n = %d\n",_n); |
---|
75 | if (type != TYPE_FLOAT) { |
---|
76 | result.addError("\"positions\" is not type float (type=%d)\n", type); |
---|
77 | return; |
---|
78 | } |
---|
79 | fprintf(stdout, "_rank = %d\n",_rank); |
---|
80 | fprintf(stdout, "_shape = %d\n",_shape); |
---|
81 | |
---|
82 | float* pos = NULL; |
---|
83 | pos = (float*) DXGetArrayData(dxpos); |
---|
84 | if (pos == NULL) { |
---|
85 | result.addError("DXGetArrayData failed to return positions array"); |
---|
86 | return; |
---|
87 | } |
---|
88 | |
---|
89 | // first call to get the number of axis needed |
---|
90 | dxgrid = (Array) DXGetComponentValue((Field) _dxobj, (char *)"connections"); |
---|
91 | if (DXQueryGridConnections(dxgrid, &_numAxis, NULL) == NULL) { |
---|
92 | // raise error, data is not a regular grid and we cannot handle it |
---|
93 | result.addError("DX says our grid is not regular, we cannot handle this data"); |
---|
94 | return; |
---|
95 | } |
---|
96 | |
---|
97 | _positions = new float[_n*_numAxis]; |
---|
98 | if (_positions == NULL) { |
---|
99 | // malloc failed, raise error |
---|
100 | result.addError("malloc of _positions array failed"); |
---|
101 | return; |
---|
102 | } |
---|
103 | memcpy(_positions,pos,sizeof(float)*_n*_numAxis); |
---|
104 | |
---|
105 | _axisLen = new int[_numAxis]; |
---|
106 | if (_axisLen == NULL) { |
---|
107 | // malloc failed, raise error |
---|
108 | result.addError("malloc of _axisLen array failed"); |
---|
109 | return; |
---|
110 | } |
---|
111 | memset(_axisLen, 0, _numAxis); |
---|
112 | |
---|
113 | _delta = new float[_numAxis*_numAxis]; |
---|
114 | if (_delta == NULL) { |
---|
115 | result.addError("malloc of _delta array failed"); |
---|
116 | return; |
---|
117 | } |
---|
118 | memset(_delta, 0, _numAxis*_numAxis); |
---|
119 | |
---|
120 | _origin = new float[_numAxis]; |
---|
121 | if (_origin == NULL) { |
---|
122 | result.addError("malloc of _origin array failed"); |
---|
123 | return; |
---|
124 | } |
---|
125 | memset(_origin, 0, _numAxis); |
---|
126 | |
---|
127 | _max = new float[_numAxis]; |
---|
128 | if (_max == NULL) { |
---|
129 | result.addError("malloc of _max array failed"); |
---|
130 | return; |
---|
131 | } |
---|
132 | memset(_max, 0, _numAxis); |
---|
133 | |
---|
134 | __findPosMax(); |
---|
135 | |
---|
136 | // parse out the gridconnections (length of each axis) array |
---|
137 | // dxgrid = (Array) DXQueryGridConnections(dxpos, &_numAxis, _axisLen); |
---|
138 | DXQueryGridPositions(dxpos, NULL, _axisLen, _origin, _delta); |
---|
139 | |
---|
140 | fprintf(stdout, "_max = [%g,%g,%g]\n",_max[0],_max[1],_max[2]); |
---|
141 | fprintf(stdout, "_delta = [%g,%g,%g]\n",_delta[0],_delta[1],_delta[2]); |
---|
142 | fprintf(stdout, " [%g,%g,%g]\n",_delta[3],_delta[4],_delta[5]); |
---|
143 | fprintf(stdout, " [%g,%g,%g]\n",_delta[6],_delta[7],_delta[8]); |
---|
144 | fprintf(stdout, "_origin = [%g,%g,%g]\n",_origin[0],_origin[1],_origin[2]); |
---|
145 | fprintf(stdout, "_axisLen = [%i,%i,%i]\n",_axisLen[0],_axisLen[1],_axisLen[2]); |
---|
146 | fflush(stdout); |
---|
147 | |
---|
148 | // grab the data array from the dx object and store it in _data |
---|
149 | dxdata = (Array) DXGetComponentValue((Field) _dxobj, (char *)"data"); |
---|
150 | DXGetArrayInfo(dxdata, NULL, &type, NULL, NULL, NULL); |
---|
151 | _data = new float[_n]; |
---|
152 | if (_data == NULL) { |
---|
153 | result.addError("malloc of _data array failed"); |
---|
154 | return; |
---|
155 | } |
---|
156 | |
---|
157 | switch (type) { |
---|
158 | case TYPE_FLOAT: |
---|
159 | float *float_data; |
---|
160 | |
---|
161 | float_data = (float *)DXGetArrayData(dxdata); |
---|
162 | memcpy(_data, float_data, sizeof(float)*_n); |
---|
163 | break; |
---|
164 | |
---|
165 | case TYPE_DOUBLE: |
---|
166 | double *double_data; |
---|
167 | double_data = (double *)DXGetArrayData(dxdata); |
---|
168 | for (int i = 0; i < _n; i++) { |
---|
169 | _data[i] = double_data[i]; |
---|
170 | } |
---|
171 | break; |
---|
172 | |
---|
173 | default: |
---|
174 | result.addError("don't know how to handle data of type %d\n", type); |
---|
175 | return; |
---|
176 | } |
---|
177 | |
---|
178 | // print debug info |
---|
179 | for (int lcv = 0, pt = 0; lcv < _n; lcv +=3, pt+=9) { |
---|
180 | fprintf(stdout, |
---|
181 | "(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n", |
---|
182 | _positions[pt],_positions[pt+1],_positions[pt+2], _data[lcv], |
---|
183 | _positions[pt+3],_positions[pt+4],_positions[pt+5],_data[lcv+1], |
---|
184 | _positions[pt+6],_positions[pt+7],_positions[pt+8],_data[lcv+2]); |
---|
185 | fflush(stdout); |
---|
186 | } |
---|
187 | __collectDataStats(); |
---|
188 | } |
---|
189 | |
---|
190 | DX::~DX() |
---|
191 | { |
---|
192 | delete[] _axisLen; |
---|
193 | delete[] _delta; |
---|
194 | delete[] _origin; |
---|
195 | delete[] _max; |
---|
196 | delete[] _data; |
---|
197 | delete[] _positions; |
---|
198 | } |
---|
199 | |
---|
200 | void |
---|
201 | DX::__findPosMax() |
---|
202 | { |
---|
203 | int lcv = 0; |
---|
204 | |
---|
205 | // initialize the max array and delta array |
---|
206 | // max holds the maximum value found for each index |
---|
207 | for (lcv = 0; lcv < _numAxis; lcv++) { |
---|
208 | _max[lcv] = _positions[lcv]; |
---|
209 | } |
---|
210 | |
---|
211 | for (lcv=lcv; lcv < _n*_numAxis; lcv++) { |
---|
212 | int maxIdx = lcv%_numAxis; |
---|
213 | if (_positions[lcv] > _max[maxIdx]) { |
---|
214 | _max[maxIdx] = _positions[lcv]; |
---|
215 | } |
---|
216 | } |
---|
217 | } |
---|
218 | |
---|
219 | void |
---|
220 | DX::__collectDataStats() |
---|
221 | { |
---|
222 | _dataMin = FLT_MAX; |
---|
223 | _dataMax = -FLT_MAX; |
---|
224 | _nzero_min = FLT_MAX; |
---|
225 | |
---|
226 | // populate _dataMin, _dataMax, _nzero_min |
---|
227 | for (int lcv = 0; lcv < _n; lcv++) { |
---|
228 | if (_data[lcv] < _dataMin) { |
---|
229 | _dataMin = _data[lcv]; |
---|
230 | } |
---|
231 | if (_data[lcv] > _dataMax) { |
---|
232 | _dataMax = _data[lcv]; |
---|
233 | } |
---|
234 | if ((_data[lcv] > 0.0f) && (_data[lcv] < _nzero_min)) { |
---|
235 | _nzero_min = _data[lcv]; |
---|
236 | } |
---|
237 | } |
---|
238 | if (_nzero_min == FLT_MAX) { |
---|
239 | fprintf(stderr, "could not find a positive minimum value\n"); |
---|
240 | fflush(stderr); |
---|
241 | } |
---|
242 | } |
---|
243 | |
---|
244 | /* |
---|
245 | * getInterpPos() |
---|
246 | * |
---|
247 | * creates a new grid of positions which can be used for interpolation. |
---|
248 | * we create the positions array using the function DXMakeGridPositionsV |
---|
249 | * which creates an n-dimensional grid of regularly spaced positions. |
---|
250 | * This function overwrites the original positions array |
---|
251 | */ |
---|
252 | void |
---|
253 | DX::__getInterpPos() |
---|
254 | { |
---|
255 | Array dxpos; |
---|
256 | float* pos = NULL; |
---|
257 | |
---|
258 | // gather the positions we want to interpolate over |
---|
259 | dxpos = DXMakeGridPositionsV(_numAxis, _axisLen, _origin, _delta); |
---|
260 | DXGetArrayInfo(dxpos, &_n, NULL, NULL, &_rank, &_shape); |
---|
261 | pos = (float*) DXGetArrayData(dxpos); |
---|
262 | if (pos == NULL) { |
---|
263 | fprintf(stdout, "DXGetArrayData failed to return positions array\n"); |
---|
264 | fflush(stdout); |
---|
265 | } |
---|
266 | |
---|
267 | if (_positions != NULL) { |
---|
268 | delete[] _positions; |
---|
269 | } |
---|
270 | _positions = new float[_n*_numAxis]; |
---|
271 | if (_positions == NULL) { |
---|
272 | // malloc failed, raise error |
---|
273 | fprintf(stdout, "malloc of _axisLen array failed"); |
---|
274 | fflush(stdout); |
---|
275 | } |
---|
276 | memcpy(_positions,pos,sizeof(float)*_n*_numAxis); |
---|
277 | |
---|
278 | pos = NULL; |
---|
279 | DXDelete((object*)dxpos); |
---|
280 | } |
---|
281 | |
---|
282 | /* |
---|
283 | * getInterpData() |
---|
284 | * |
---|
285 | * this function interpolates over a positions array to produce data for each |
---|
286 | * point in the positions array. we use the position data stored in _positions |
---|
287 | * array. |
---|
288 | */ |
---|
289 | void |
---|
290 | DX::__getInterpData() |
---|
291 | { |
---|
292 | int pts = _n; |
---|
293 | int interppts = pts; |
---|
294 | Interpolator interpolator; |
---|
295 | |
---|
296 | _data = new float[_n]; |
---|
297 | if (_data == NULL) { |
---|
298 | // malloc failed, raise error |
---|
299 | fprintf(stdout, "malloc of _data array failed"); |
---|
300 | fflush(stdout); |
---|
301 | } |
---|
302 | memset(_data,0,_n); |
---|
303 | |
---|
304 | // build the interpolator and interpolate |
---|
305 | fprintf(stdout, "creating DXNewInterpolator...\n"); |
---|
306 | fflush(stdout); |
---|
307 | interpolator = DXNewInterpolator(_dxobj,INTERP_INIT_IMMEDIATE,-1.0); |
---|
308 | fprintf(stdout,"_rank = %i\n",_rank); |
---|
309 | fprintf(stdout,"_shape = %i\n",_shape); |
---|
310 | fprintf(stdout,"_n = %i\n",_n); |
---|
311 | fprintf(stdout,"start interppts = %i\n",interppts); |
---|
312 | fflush(stdout); |
---|
313 | DXInterpolate(interpolator,&interppts,_positions,_data); |
---|
314 | fprintf(stdout,"interppts = %i\n",interppts); |
---|
315 | fflush(stdout); |
---|
316 | |
---|
317 | // print debug info |
---|
318 | for (int lcv = 0, pt = 0; lcv < pts; lcv+=3,pt+=9) { |
---|
319 | fprintf(stdout, |
---|
320 | "(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n", |
---|
321 | _positions[pt],_positions[pt+1],_positions[pt+2], _data[lcv], |
---|
322 | _positions[pt+3],_positions[pt+4],_positions[pt+5],_data[lcv+1], |
---|
323 | _positions[pt+6],_positions[pt+7],_positions[pt+8],_data[lcv+2]); |
---|
324 | fflush(stdout); |
---|
325 | } |
---|
326 | |
---|
327 | __collectDataStats(); |
---|
328 | } |
---|
329 | |
---|
330 | /* |
---|
331 | * interpolate() |
---|
332 | * |
---|
333 | * generate a positions array with optional new axis length and |
---|
334 | * interpolate to get the new data values at each point in |
---|
335 | * the positions array. this function currently only works if you |
---|
336 | * do not change the axis length (i.e. newAxisLen == NULL). |
---|
337 | */ |
---|
338 | DX& |
---|
339 | DX::interpolate(int* newAxisLen) |
---|
340 | { |
---|
341 | fprintf(stdout, "----begin interpolation----\n"); |
---|
342 | fflush(stdout); |
---|
343 | if (newAxisLen != NULL) { |
---|
344 | for (int i = 0; i < _numAxis; i++) { |
---|
345 | _axisLen[i] = newAxisLen[i]; |
---|
346 | } |
---|
347 | } |
---|
348 | __getInterpPos(); |
---|
349 | __getInterpData(); |
---|
350 | fprintf(stdout, "----end interpolation----\n"); |
---|
351 | fflush(stdout); |
---|
352 | return *this; |
---|
353 | } |
---|
354 | |
---|