File: | gromacs/fileio/libxdrf.c |
Location: | line 608, column 9 |
Description: | Value stored to 'lip' is never read |
1 | /* |
2 | * This file is part of the GROMACS molecular simulation package. |
3 | * |
4 | * Copyright (c) 1991-2000, University of Groningen, The Netherlands. |
5 | * Copyright (c) 2001-2004, The GROMACS development team. |
6 | * Copyright (c) 2013,2014, by the GROMACS development team, led by |
7 | * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, |
8 | * and including many others, as listed in the AUTHORS file in the |
9 | * top-level source directory and at http://www.gromacs.org. |
10 | * |
11 | * GROMACS is free software; you can redistribute it and/or |
12 | * modify it under the terms of the GNU Lesser General Public License |
13 | * as published by the Free Software Foundation; either version 2.1 |
14 | * of the License, or (at your option) any later version. |
15 | * |
16 | * GROMACS is distributed in the hope that it will be useful, |
17 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
18 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
19 | * Lesser General Public License for more details. |
20 | * |
21 | * You should have received a copy of the GNU Lesser General Public |
22 | * License along with GROMACS; if not, see |
23 | * http://www.gnu.org/licenses, or write to the Free Software Foundation, |
24 | * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. |
25 | * |
26 | * If you want to redistribute modifications to GROMACS, please |
27 | * consider that scientific software is very special. Version |
28 | * control is crucial - bugs must be traceable. We will be happy to |
29 | * consider code for inclusion in the official distribution, but |
30 | * derived work must not be called official GROMACS. Details are found |
31 | * in the README & COPYING files - if they are missing, get the |
32 | * official version at http://www.gromacs.org. |
33 | * |
34 | * To help us fund GROMACS development, we humbly ask that you cite |
35 | * the research papers on the package. Check out http://www.gromacs.org. |
36 | */ |
37 | #ifdef HAVE_CONFIG_H1 |
38 | #include <config.h> |
39 | #endif |
40 | |
41 | #include <limits.h> |
42 | #include <math.h> |
43 | #include <stdio.h> |
44 | #include <stdlib.h> |
45 | #include <string.h> |
46 | |
47 | #include "xdrf.h" |
48 | #include "xdr_datatype.h" |
49 | #include "gromacs/utility/futil.h" |
50 | |
51 | /* This is just for clarity - it can never be anything but 4! */ |
52 | #define XDR_INT_SIZE4 4 |
53 | |
54 | /* same order as the definition of xdr_datatype */ |
55 | const char *xdr_datatype_names[] = |
56 | { |
57 | "int", |
58 | "float", |
59 | "double", |
60 | "large int", |
61 | "char", |
62 | "string" |
63 | }; |
64 | |
65 | |
66 | /*___________________________________________________________________________ |
67 | | |
68 | | what follows are the C routine to read/write compressed coordinates together |
69 | | with some routines to assist in this task (those are marked |
70 | | static and cannot be called from user programs) |
71 | */ |
72 | #define MAXABS2147483647 -2 INT_MAX2147483647-2 |
73 | |
74 | #ifndef MIN |
75 | #define MIN(x, y)(((x)<(y))?(x):(y)) ((x) < (y) ? (x) : (y)) |
76 | #endif |
77 | #ifndef MAX |
78 | #define MAX(x, y)(((x)>(y))?(x):(y)) ((x) > (y) ? (x) : (y)) |
79 | #endif |
80 | #ifndef SQR |
81 | #define SQR(x)((x)*(x)) ((x)*(x)) |
82 | #endif |
83 | static const int magicints[] = { |
84 | 0, 0, 0, 0, 0, 0, 0, 0, 0, |
85 | 8, 10, 12, 16, 20, 25, 32, 40, 50, 64, |
86 | 80, 101, 128, 161, 203, 256, 322, 406, 512, 645, |
87 | 812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501, |
88 | 8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536, |
89 | 82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561, |
90 | 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304, 5284491, 6658042, |
91 | 8388607, 10568983, 13316085, 16777216 |
92 | }; |
93 | |
94 | #define FIRSTIDX9 9 |
95 | /* note that magicints[FIRSTIDX-1] == 0 */ |
96 | #define LASTIDX(sizeof(magicints) / sizeof(*magicints)) (sizeof(magicints) / sizeof(*magicints)) |
97 | |
98 | |
99 | /*____________________________________________________________________________ |
100 | | |
101 | | sendbits - encode num into buf using the specified number of bits |
102 | | |
103 | | This routines appends the value of num to the bits already present in |
104 | | the array buf. You need to give it the number of bits to use and you |
105 | | better make sure that this number of bits is enough to hold the value |
106 | | Also num must be positive. |
107 | | |
108 | */ |
109 | |
110 | static void sendbits(int buf[], int num_of_bits, int num) |
111 | { |
112 | |
113 | unsigned int cnt, lastbyte; |
114 | int lastbits; |
115 | unsigned char * cbuf; |
116 | |
117 | cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf); |
118 | cnt = (unsigned int) buf[0]; |
119 | lastbits = buf[1]; |
120 | lastbyte = (unsigned int) buf[2]; |
121 | while (num_of_bits >= 8) |
122 | { |
123 | lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/); |
124 | cbuf[cnt++] = lastbyte >> lastbits; |
125 | num_of_bits -= 8; |
126 | } |
127 | if (num_of_bits > 0) |
128 | { |
129 | lastbyte = (lastbyte << num_of_bits) | num; |
130 | lastbits += num_of_bits; |
131 | if (lastbits >= 8) |
132 | { |
133 | lastbits -= 8; |
134 | cbuf[cnt++] = lastbyte >> lastbits; |
135 | } |
136 | } |
137 | buf[0] = cnt; |
138 | buf[1] = lastbits; |
139 | buf[2] = lastbyte; |
140 | if (lastbits > 0) |
141 | { |
142 | cbuf[cnt] = lastbyte << (8 - lastbits); |
143 | } |
144 | } |
145 | |
146 | /*_________________________________________________________________________ |
147 | | |
148 | | sizeofint - calculate bitsize of an integer |
149 | | |
150 | | return the number of bits needed to store an integer with given max size |
151 | | |
152 | */ |
153 | |
154 | static int sizeofint(const int size) |
155 | { |
156 | int num = 1; |
157 | int num_of_bits = 0; |
158 | |
159 | while (size >= num && num_of_bits < 32) |
160 | { |
161 | num_of_bits++; |
162 | num <<= 1; |
163 | } |
164 | return num_of_bits; |
165 | } |
166 | |
167 | /*___________________________________________________________________________ |
168 | | |
169 | | sizeofints - calculate 'bitsize' of compressed ints |
170 | | |
171 | | given the number of small unsigned integers and the maximum value |
172 | | return the number of bits needed to read or write them with the |
173 | | routines receiveints and sendints. You need this parameter when |
174 | | calling these routines. Note that for many calls I can use |
175 | | the variable 'smallidx' which is exactly the number of bits, and |
176 | | So I don't need to call 'sizeofints for those calls. |
177 | */ |
178 | |
179 | static int sizeofints( const int num_of_ints, unsigned int sizes[]) |
180 | { |
181 | int i, num; |
182 | int bytes[32]; |
183 | unsigned int num_of_bytes, num_of_bits, bytecnt, tmp; |
184 | num_of_bytes = 1; |
185 | bytes[0] = 1; |
186 | num_of_bits = 0; |
187 | for (i = 0; i < num_of_ints; i++) |
188 | { |
189 | tmp = 0; |
190 | for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) |
191 | { |
192 | tmp = bytes[bytecnt] * sizes[i] + tmp; |
193 | bytes[bytecnt] = tmp & 0xff; |
194 | tmp >>= 8; |
195 | } |
196 | while (tmp != 0) |
197 | { |
198 | bytes[bytecnt++] = tmp & 0xff; |
199 | tmp >>= 8; |
200 | } |
201 | num_of_bytes = bytecnt; |
202 | } |
203 | num = 1; |
204 | num_of_bytes--; |
205 | while (bytes[num_of_bytes] >= num) |
206 | { |
207 | num_of_bits++; |
208 | num *= 2; |
209 | } |
210 | return num_of_bits + num_of_bytes * 8; |
211 | |
212 | } |
213 | |
214 | /*____________________________________________________________________________ |
215 | | |
216 | | sendints - send a small set of small integers in compressed format |
217 | | |
218 | | this routine is used internally by xdr3dfcoord, to send a set of |
219 | | small integers to the buffer. |
220 | | Multiplication with fixed (specified maximum ) sizes is used to get |
221 | | to one big, multibyte integer. Allthough the routine could be |
222 | | modified to handle sizes bigger than 16777216, or more than just |
223 | | a few integers, this is not done, because the gain in compression |
224 | | isn't worth the effort. Note that overflowing the multiplication |
225 | | or the byte buffer (32 bytes) is unchecked and causes bad results. |
226 | | |
227 | */ |
228 | |
229 | static void sendints(int buf[], const int num_of_ints, const int num_of_bits, |
230 | unsigned int sizes[], unsigned int nums[]) |
231 | { |
232 | |
233 | int i, num_of_bytes, bytecnt; |
234 | unsigned int bytes[32], tmp; |
235 | |
236 | tmp = nums[0]; |
237 | num_of_bytes = 0; |
238 | do |
239 | { |
240 | bytes[num_of_bytes++] = tmp & 0xff; |
241 | tmp >>= 8; |
242 | } |
243 | while (tmp != 0); |
244 | |
245 | for (i = 1; i < num_of_ints; i++) |
246 | { |
247 | if (nums[i] >= sizes[i]) |
248 | { |
249 | fprintf(stderrstderr, "major breakdown in sendints num %u doesn't " |
250 | "match size %u\n", nums[i], sizes[i]); |
251 | exit(1); |
252 | } |
253 | /* use one step multiply */ |
254 | tmp = nums[i]; |
255 | for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) |
256 | { |
257 | tmp = bytes[bytecnt] * sizes[i] + tmp; |
258 | bytes[bytecnt] = tmp & 0xff; |
259 | tmp >>= 8; |
260 | } |
261 | while (tmp != 0) |
262 | { |
263 | bytes[bytecnt++] = tmp & 0xff; |
264 | tmp >>= 8; |
265 | } |
266 | num_of_bytes = bytecnt; |
267 | } |
268 | if (num_of_bits >= num_of_bytes * 8) |
269 | { |
270 | for (i = 0; i < num_of_bytes; i++) |
271 | { |
272 | sendbits(buf, 8, bytes[i]); |
273 | } |
274 | sendbits(buf, num_of_bits - num_of_bytes * 8, 0); |
275 | } |
276 | else |
277 | { |
278 | for (i = 0; i < num_of_bytes-1; i++) |
279 | { |
280 | sendbits(buf, 8, bytes[i]); |
281 | } |
282 | sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]); |
283 | } |
284 | } |
285 | |
286 | |
287 | /*___________________________________________________________________________ |
288 | | |
289 | | receivebits - decode number from buf using specified number of bits |
290 | | |
291 | | extract the number of bits from the array buf and construct an integer |
292 | | from it. Return that value. |
293 | | |
294 | */ |
295 | |
296 | static int receivebits(int buf[], int num_of_bits) |
297 | { |
298 | |
299 | int cnt, num, lastbits; |
300 | unsigned int lastbyte; |
301 | unsigned char * cbuf; |
302 | int mask = (1 << num_of_bits) -1; |
303 | |
304 | cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf); |
305 | cnt = buf[0]; |
306 | lastbits = (unsigned int) buf[1]; |
307 | lastbyte = (unsigned int) buf[2]; |
308 | |
309 | num = 0; |
310 | while (num_of_bits >= 8) |
311 | { |
312 | lastbyte = ( lastbyte << 8 ) | cbuf[cnt++]; |
313 | num |= (lastbyte >> lastbits) << (num_of_bits - 8); |
314 | num_of_bits -= 8; |
315 | } |
316 | if (num_of_bits > 0) |
317 | { |
318 | if (lastbits < num_of_bits) |
319 | { |
320 | lastbits += 8; |
321 | lastbyte = (lastbyte << 8) | cbuf[cnt++]; |
322 | } |
323 | lastbits -= num_of_bits; |
324 | num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1); |
325 | } |
326 | num &= mask; |
327 | buf[0] = cnt; |
328 | buf[1] = lastbits; |
329 | buf[2] = lastbyte; |
330 | return num; |
331 | } |
332 | |
333 | /*____________________________________________________________________________ |
334 | | |
335 | | receiveints - decode 'small' integers from the buf array |
336 | | |
337 | | this routine is the inverse from sendints() and decodes the small integers |
338 | | written to buf by calculating the remainder and doing divisions with |
339 | | the given sizes[]. You need to specify the total number of bits to be |
340 | | used from buf in num_of_bits. |
341 | | |
342 | */ |
343 | |
344 | static void receiveints(int buf[], const int num_of_ints, int num_of_bits, |
345 | unsigned int sizes[], int nums[]) |
346 | { |
347 | int bytes[32]; |
348 | int i, j, num_of_bytes, p, num; |
349 | |
350 | bytes[0] = bytes[1] = bytes[2] = bytes[3] = 0; |
351 | num_of_bytes = 0; |
352 | while (num_of_bits > 8) |
353 | { |
354 | bytes[num_of_bytes++] = receivebits(buf, 8); |
355 | num_of_bits -= 8; |
356 | } |
357 | if (num_of_bits > 0) |
358 | { |
359 | bytes[num_of_bytes++] = receivebits(buf, num_of_bits); |
360 | } |
361 | for (i = num_of_ints-1; i > 0; i--) |
362 | { |
363 | num = 0; |
364 | for (j = num_of_bytes-1; j >= 0; j--) |
365 | { |
366 | num = (num << 8) | bytes[j]; |
367 | p = num / sizes[i]; |
368 | bytes[j] = p; |
369 | num = num - p * sizes[i]; |
370 | } |
371 | nums[i] = num; |
372 | } |
373 | nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24); |
374 | } |
375 | |
376 | /*____________________________________________________________________________ |
377 | | |
378 | | xdr3dfcoord - read or write compressed 3d coordinates to xdr file. |
379 | | |
380 | | this routine reads or writes (depending on how you opened the file with |
381 | | xdropen() ) a large number of 3d coordinates (stored in *fp). |
382 | | The number of coordinates triplets to write is given by *size. On |
383 | | read this number may be zero, in which case it reads as many as were written |
384 | | or it may specify the number if triplets to read (which should match the |
385 | | number written). |
386 | | Compression is achieved by first converting all floating numbers to integer |
387 | | using multiplication by *precision and rounding to the nearest integer. |
388 | | Then the minimum and maximum value are calculated to determine the range. |
389 | | The limited range of integers so found, is used to compress the coordinates. |
390 | | In addition the differences between succesive coordinates is calculated. |
391 | | If the difference happens to be 'small' then only the difference is saved, |
392 | | compressing the data even more. The notion of 'small' is changed dynamically |
393 | | and is enlarged or reduced whenever needed or possible. |
394 | | Extra compression is achieved in the case of GROMOS and coordinates of |
395 | | water molecules. GROMOS first writes out the Oxygen position, followed by |
396 | | the two hydrogens. In order to make the differences smaller (and thereby |
397 | | compression the data better) the order is changed into first one hydrogen |
398 | | then the oxygen, followed by the other hydrogen. This is rather special, but |
399 | | it shouldn't harm in the general case. |
400 | | |
401 | */ |
402 | |
403 | int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) |
404 | { |
405 | int *ip = NULL((void*)0); |
406 | int *buf = NULL((void*)0); |
407 | gmx_bool bRead; |
408 | |
409 | /* preallocate a small buffer and ip on the stack - if we need more |
410 | we can always malloc(). This is faster for small values of size: */ |
411 | unsigned prealloc_size = 3*16; |
412 | int prealloc_ip[3*16], prealloc_buf[3*20]; |
413 | int we_should_free = 0; |
414 | |
415 | int minint[3], maxint[3], mindiff, *lip, diff; |
416 | int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx; |
417 | int minidx, maxidx; |
418 | unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip; |
419 | int flag, k; |
420 | int smallnum, smaller, larger, i, is_small, is_smaller, run, prevrun; |
421 | float *lfp, lf; |
422 | int tmp, *thiscoord, prevcoord[3]; |
423 | unsigned int tmpcoord[30]; |
424 | |
425 | int bufsize, xdrid, lsize; |
426 | unsigned int bitsize; |
427 | float inv_precision; |
428 | int errval = 1; |
429 | int rc; |
430 | |
431 | bRead = (xdrs->x_op == XDR_DECODE); |
432 | bitsizeint[0] = bitsizeint[1] = bitsizeint[2] = 0; |
433 | prevcoord[0] = prevcoord[1] = prevcoord[2] = 0; |
434 | |
435 | if (!bRead) |
436 | { |
437 | /* xdrs is open for writing */ |
438 | |
439 | if (xdr_int(xdrs, size) == 0) |
440 | { |
441 | return 0; |
442 | } |
443 | size3 = *size * 3; |
444 | /* when the number of coordinates is small, don't try to compress; just |
445 | * write them as floats using xdr_vector |
446 | */ |
447 | if (*size <= 9) |
448 | { |
449 | return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3, |
450 | (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float)); |
451 | } |
452 | |
453 | if (xdr_float(xdrs, precision) == 0) |
454 | { |
455 | return 0; |
456 | } |
457 | |
458 | if (size3 <= prealloc_size) |
459 | { |
460 | ip = prealloc_ip; |
461 | buf = prealloc_buf; |
462 | } |
463 | else |
464 | { |
465 | we_should_free = 1; |
466 | bufsize = size3 * 1.2; |
467 | ip = (int *)malloc((size_t)(size3 * sizeof(*ip))); |
468 | buf = (int *)malloc((size_t)(bufsize * sizeof(*buf))); |
469 | if (ip == NULL((void*)0) || buf == NULL((void*)0)) |
470 | { |
471 | fprintf(stderrstderr, "malloc failed\n"); |
472 | exit(1); |
473 | } |
474 | } |
475 | /* buf[0-2] are special and do not contain actual data */ |
476 | buf[0] = buf[1] = buf[2] = 0; |
477 | minint[0] = minint[1] = minint[2] = INT_MAX2147483647; |
478 | maxint[0] = maxint[1] = maxint[2] = INT_MIN(-2147483647 -1); |
479 | prevrun = -1; |
480 | lfp = fp; |
481 | lip = ip; |
482 | mindiff = INT_MAX2147483647; |
483 | oldlint1 = oldlint2 = oldlint3 = 0; |
484 | while (lfp < fp + size3) |
485 | { |
486 | /* find nearest integer */ |
487 | if (*lfp >= 0.0) |
488 | { |
489 | lf = *lfp * *precision + 0.5; |
490 | } |
491 | else |
492 | { |
493 | lf = *lfp * *precision - 0.5; |
494 | } |
495 | if (fabs(lf) > MAXABS2147483647 -2) |
496 | { |
497 | /* scaling would cause overflow */ |
498 | errval = 0; |
499 | } |
500 | lint1 = lf; |
501 | if (lint1 < minint[0]) |
502 | { |
503 | minint[0] = lint1; |
504 | } |
505 | if (lint1 > maxint[0]) |
506 | { |
507 | maxint[0] = lint1; |
508 | } |
509 | *lip++ = lint1; |
510 | lfp++; |
511 | if (*lfp >= 0.0) |
512 | { |
513 | lf = *lfp * *precision + 0.5; |
514 | } |
515 | else |
516 | { |
517 | lf = *lfp * *precision - 0.5; |
518 | } |
519 | if (fabs(lf) > MAXABS2147483647 -2) |
520 | { |
521 | /* scaling would cause overflow */ |
522 | errval = 0; |
523 | } |
524 | lint2 = lf; |
525 | if (lint2 < minint[1]) |
526 | { |
527 | minint[1] = lint2; |
528 | } |
529 | if (lint2 > maxint[1]) |
530 | { |
531 | maxint[1] = lint2; |
532 | } |
533 | *lip++ = lint2; |
534 | lfp++; |
535 | if (*lfp >= 0.0) |
536 | { |
537 | lf = *lfp * *precision + 0.5; |
538 | } |
539 | else |
540 | { |
541 | lf = *lfp * *precision - 0.5; |
542 | } |
543 | if (fabs(lf) > MAXABS2147483647 -2) |
544 | { |
545 | /* scaling would cause overflow */ |
546 | errval = 0; |
547 | } |
548 | lint3 = lf; |
549 | if (lint3 < minint[2]) |
550 | { |
551 | minint[2] = lint3; |
552 | } |
553 | if (lint3 > maxint[2]) |
554 | { |
555 | maxint[2] = lint3; |
556 | } |
557 | *lip++ = lint3; |
558 | lfp++; |
559 | diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3); |
560 | if (diff < mindiff && lfp > fp + 3) |
561 | { |
562 | mindiff = diff; |
563 | } |
564 | oldlint1 = lint1; |
565 | oldlint2 = lint2; |
566 | oldlint3 = lint3; |
567 | } |
568 | if ( (xdr_int(xdrs, &(minint[0])) == 0) || |
569 | (xdr_int(xdrs, &(minint[1])) == 0) || |
570 | (xdr_int(xdrs, &(minint[2])) == 0) || |
571 | (xdr_int(xdrs, &(maxint[0])) == 0) || |
572 | (xdr_int(xdrs, &(maxint[1])) == 0) || |
573 | (xdr_int(xdrs, &(maxint[2])) == 0)) |
574 | { |
575 | if (we_should_free) |
576 | { |
577 | free(ip); |
578 | free(buf); |
579 | } |
580 | return 0; |
581 | } |
582 | |
583 | if ((float)maxint[0] - (float)minint[0] >= MAXABS2147483647 -2 || |
584 | (float)maxint[1] - (float)minint[1] >= MAXABS2147483647 -2 || |
585 | (float)maxint[2] - (float)minint[2] >= MAXABS2147483647 -2) |
586 | { |
587 | /* turning value in unsigned by subtracting minint |
588 | * would cause overflow |
589 | */ |
590 | errval = 0; |
591 | } |
592 | sizeint[0] = maxint[0] - minint[0]+1; |
593 | sizeint[1] = maxint[1] - minint[1]+1; |
594 | sizeint[2] = maxint[2] - minint[2]+1; |
595 | |
596 | /* check if one of the sizes is to big to be multiplied */ |
597 | if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) |
598 | { |
599 | bitsizeint[0] = sizeofint(sizeint[0]); |
600 | bitsizeint[1] = sizeofint(sizeint[1]); |
601 | bitsizeint[2] = sizeofint(sizeint[2]); |
602 | bitsize = 0; /* flag the use of large sizes */ |
603 | } |
604 | else |
605 | { |
606 | bitsize = sizeofints(3, sizeint); |
607 | } |
608 | lip = ip; |
Value stored to 'lip' is never read | |
609 | luip = (unsigned int *) ip; |
610 | smallidx = FIRSTIDX9; |
611 | while (smallidx < LASTIDX(sizeof(magicints) / sizeof(*magicints)) && magicints[smallidx] < mindiff) |
612 | { |
613 | smallidx++; |
614 | } |
615 | if (xdr_int(xdrs, &smallidx) == 0) |
616 | { |
617 | if (we_should_free) |
618 | { |
619 | free(ip); |
620 | free(buf); |
621 | } |
622 | return 0; |
623 | } |
624 | |
625 | maxidx = MIN(LASTIDX, smallidx + 8)((((sizeof(magicints) / sizeof(*magicints)))<(smallidx + 8 ))?((sizeof(magicints) / sizeof(*magicints))):(smallidx + 8)); |
626 | minidx = maxidx - 8; /* often this equal smallidx */ |
627 | smaller = magicints[MAX(FIRSTIDX, smallidx-1)(((9)>(smallidx-1))?(9):(smallidx-1))] / 2; |
628 | smallnum = magicints[smallidx] / 2; |
629 | sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx]; |
630 | larger = magicints[maxidx] / 2; |
631 | i = 0; |
632 | while (i < *size) |
633 | { |
634 | is_small = 0; |
635 | thiscoord = (int *)(luip) + i * 3; |
636 | if (smallidx < maxidx && i >= 1 && |
637 | abs(thiscoord[0] - prevcoord[0]) < larger && |
638 | abs(thiscoord[1] - prevcoord[1]) < larger && |
639 | abs(thiscoord[2] - prevcoord[2]) < larger) |
640 | { |
641 | is_smaller = 1; |
642 | } |
643 | else if (smallidx > minidx) |
644 | { |
645 | is_smaller = -1; |
646 | } |
647 | else |
648 | { |
649 | is_smaller = 0; |
650 | } |
651 | if (i + 1 < *size) |
652 | { |
653 | if (abs(thiscoord[0] - thiscoord[3]) < smallnum && |
654 | abs(thiscoord[1] - thiscoord[4]) < smallnum && |
655 | abs(thiscoord[2] - thiscoord[5]) < smallnum) |
656 | { |
657 | /* interchange first with second atom for better |
658 | * compression of water molecules |
659 | */ |
660 | tmp = thiscoord[0]; thiscoord[0] = thiscoord[3]; |
661 | thiscoord[3] = tmp; |
662 | tmp = thiscoord[1]; thiscoord[1] = thiscoord[4]; |
663 | thiscoord[4] = tmp; |
664 | tmp = thiscoord[2]; thiscoord[2] = thiscoord[5]; |
665 | thiscoord[5] = tmp; |
666 | is_small = 1; |
667 | } |
668 | |
669 | } |
670 | tmpcoord[0] = thiscoord[0] - minint[0]; |
671 | tmpcoord[1] = thiscoord[1] - minint[1]; |
672 | tmpcoord[2] = thiscoord[2] - minint[2]; |
673 | if (bitsize == 0) |
674 | { |
675 | sendbits(buf, bitsizeint[0], tmpcoord[0]); |
676 | sendbits(buf, bitsizeint[1], tmpcoord[1]); |
677 | sendbits(buf, bitsizeint[2], tmpcoord[2]); |
678 | } |
679 | else |
680 | { |
681 | sendints(buf, 3, bitsize, sizeint, tmpcoord); |
682 | } |
683 | prevcoord[0] = thiscoord[0]; |
684 | prevcoord[1] = thiscoord[1]; |
685 | prevcoord[2] = thiscoord[2]; |
686 | thiscoord = thiscoord + 3; |
687 | i++; |
688 | |
689 | run = 0; |
690 | if (is_small == 0 && is_smaller == -1) |
691 | { |
692 | is_smaller = 0; |
693 | } |
694 | while (is_small && run < 8*3) |
695 | { |
696 | if (is_smaller == -1 && ( |
697 | SQR(thiscoord[0] - prevcoord[0])((thiscoord[0] - prevcoord[0])*(thiscoord[0] - prevcoord[0])) + |
698 | SQR(thiscoord[1] - prevcoord[1])((thiscoord[1] - prevcoord[1])*(thiscoord[1] - prevcoord[1])) + |
699 | SQR(thiscoord[2] - prevcoord[2])((thiscoord[2] - prevcoord[2])*(thiscoord[2] - prevcoord[2])) >= smaller * smaller)) |
700 | { |
701 | is_smaller = 0; |
702 | } |
703 | |
704 | tmpcoord[run++] = thiscoord[0] - prevcoord[0] + smallnum; |
705 | tmpcoord[run++] = thiscoord[1] - prevcoord[1] + smallnum; |
706 | tmpcoord[run++] = thiscoord[2] - prevcoord[2] + smallnum; |
707 | |
708 | prevcoord[0] = thiscoord[0]; |
709 | prevcoord[1] = thiscoord[1]; |
710 | prevcoord[2] = thiscoord[2]; |
711 | |
712 | i++; |
713 | thiscoord = thiscoord + 3; |
714 | is_small = 0; |
715 | if (i < *size && |
716 | abs(thiscoord[0] - prevcoord[0]) < smallnum && |
717 | abs(thiscoord[1] - prevcoord[1]) < smallnum && |
718 | abs(thiscoord[2] - prevcoord[2]) < smallnum) |
719 | { |
720 | is_small = 1; |
721 | } |
722 | } |
723 | if (run != prevrun || is_smaller != 0) |
724 | { |
725 | prevrun = run; |
726 | sendbits(buf, 1, 1); /* flag the change in run-length */ |
727 | sendbits(buf, 5, run+is_smaller+1); |
728 | } |
729 | else |
730 | { |
731 | sendbits(buf, 1, 0); /* flag the fact that runlength did not change */ |
732 | } |
733 | for (k = 0; k < run; k += 3) |
734 | { |
735 | sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]); |
736 | } |
737 | if (is_smaller != 0) |
738 | { |
739 | smallidx += is_smaller; |
740 | if (is_smaller < 0) |
741 | { |
742 | smallnum = smaller; |
743 | smaller = magicints[smallidx-1] / 2; |
744 | } |
745 | else |
746 | { |
747 | smaller = smallnum; |
748 | smallnum = magicints[smallidx] / 2; |
749 | } |
750 | sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx]; |
751 | } |
752 | } |
753 | if (buf[1] != 0) |
754 | { |
755 | buf[0]++; |
756 | } |
757 | /* buf[0] holds the length in bytes */ |
758 | if (xdr_int(xdrs, &(buf[0])) == 0) |
759 | { |
760 | if (we_should_free) |
761 | { |
762 | free(ip); |
763 | free(buf); |
764 | } |
765 | return 0; |
766 | } |
767 | |
768 | |
769 | rc = errval * (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0])); |
770 | if (we_should_free) |
771 | { |
772 | free(ip); |
773 | free(buf); |
774 | } |
775 | return rc; |
776 | |
777 | } |
778 | else |
779 | { |
780 | |
781 | /* xdrs is open for reading */ |
782 | |
783 | if (xdr_int(xdrs, &lsize) == 0) |
784 | { |
785 | return 0; |
786 | } |
787 | if (*size != 0 && lsize != *size) |
788 | { |
789 | fprintf(stderrstderr, "wrong number of coordinates in xdr3dfcoord; " |
790 | "%d arg vs %d in file", *size, lsize); |
791 | } |
792 | *size = lsize; |
793 | size3 = *size * 3; |
794 | if (*size <= 9) |
795 | { |
796 | *precision = -1; |
797 | return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3, |
798 | (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float)); |
799 | } |
800 | if (xdr_float(xdrs, precision) == 0) |
801 | { |
802 | return 0; |
803 | } |
804 | |
805 | if (size3 <= prealloc_size) |
806 | { |
807 | ip = prealloc_ip; |
808 | buf = prealloc_buf; |
809 | } |
810 | else |
811 | { |
812 | we_should_free = 1; |
813 | bufsize = size3 * 1.2; |
814 | ip = (int *)malloc((size_t)(size3 * sizeof(*ip))); |
815 | buf = (int *)malloc((size_t)(bufsize * sizeof(*buf))); |
816 | if (ip == NULL((void*)0) || buf == NULL((void*)0)) |
817 | { |
818 | fprintf(stderrstderr, "malloc failed\n"); |
819 | exit(1); |
820 | } |
821 | } |
822 | |
823 | buf[0] = buf[1] = buf[2] = 0; |
824 | |
825 | if ( (xdr_int(xdrs, &(minint[0])) == 0) || |
826 | (xdr_int(xdrs, &(minint[1])) == 0) || |
827 | (xdr_int(xdrs, &(minint[2])) == 0) || |
828 | (xdr_int(xdrs, &(maxint[0])) == 0) || |
829 | (xdr_int(xdrs, &(maxint[1])) == 0) || |
830 | (xdr_int(xdrs, &(maxint[2])) == 0)) |
831 | { |
832 | if (we_should_free) |
833 | { |
834 | free(ip); |
835 | free(buf); |
836 | } |
837 | return 0; |
838 | } |
839 | |
840 | sizeint[0] = maxint[0] - minint[0]+1; |
841 | sizeint[1] = maxint[1] - minint[1]+1; |
842 | sizeint[2] = maxint[2] - minint[2]+1; |
843 | |
844 | /* check if one of the sizes is to big to be multiplied */ |
845 | if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) |
846 | { |
847 | bitsizeint[0] = sizeofint(sizeint[0]); |
848 | bitsizeint[1] = sizeofint(sizeint[1]); |
849 | bitsizeint[2] = sizeofint(sizeint[2]); |
850 | bitsize = 0; /* flag the use of large sizes */ |
851 | } |
852 | else |
853 | { |
854 | bitsize = sizeofints(3, sizeint); |
855 | } |
856 | |
857 | if (xdr_int(xdrs, &smallidx) == 0) |
858 | { |
859 | if (we_should_free) |
860 | { |
861 | free(ip); |
862 | free(buf); |
863 | } |
864 | return 0; |
865 | } |
866 | |
867 | maxidx = MIN(LASTIDX, smallidx + 8)((((sizeof(magicints) / sizeof(*magicints)))<(smallidx + 8 ))?((sizeof(magicints) / sizeof(*magicints))):(smallidx + 8)); |
868 | minidx = maxidx - 8; /* often this equal smallidx */ |
869 | smaller = magicints[MAX(FIRSTIDX, smallidx-1)(((9)>(smallidx-1))?(9):(smallidx-1))] / 2; |
870 | smallnum = magicints[smallidx] / 2; |
871 | sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx]; |
872 | larger = magicints[maxidx]; |
873 | |
874 | /* buf[0] holds the length in bytes */ |
875 | |
876 | if (xdr_int(xdrs, &(buf[0])) == 0) |
877 | { |
878 | if (we_should_free) |
879 | { |
880 | free(ip); |
881 | free(buf); |
882 | } |
883 | return 0; |
884 | } |
885 | |
886 | |
887 | if (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]) == 0) |
888 | { |
889 | if (we_should_free) |
890 | { |
891 | free(ip); |
892 | free(buf); |
893 | } |
894 | return 0; |
895 | } |
896 | |
897 | |
898 | |
899 | buf[0] = buf[1] = buf[2] = 0; |
900 | |
901 | lfp = fp; |
902 | inv_precision = 1.0 / *precision; |
903 | run = 0; |
904 | i = 0; |
905 | lip = ip; |
906 | while (i < lsize) |
907 | { |
908 | thiscoord = (int *)(lip) + i * 3; |
909 | |
910 | if (bitsize == 0) |
911 | { |
912 | thiscoord[0] = receivebits(buf, bitsizeint[0]); |
913 | thiscoord[1] = receivebits(buf, bitsizeint[1]); |
914 | thiscoord[2] = receivebits(buf, bitsizeint[2]); |
915 | } |
916 | else |
917 | { |
918 | receiveints(buf, 3, bitsize, sizeint, thiscoord); |
919 | } |
920 | |
921 | i++; |
922 | thiscoord[0] += minint[0]; |
923 | thiscoord[1] += minint[1]; |
924 | thiscoord[2] += minint[2]; |
925 | |
926 | prevcoord[0] = thiscoord[0]; |
927 | prevcoord[1] = thiscoord[1]; |
928 | prevcoord[2] = thiscoord[2]; |
929 | |
930 | |
931 | flag = receivebits(buf, 1); |
932 | is_smaller = 0; |
933 | if (flag == 1) |
934 | { |
935 | run = receivebits(buf, 5); |
936 | is_smaller = run % 3; |
937 | run -= is_smaller; |
938 | is_smaller--; |
939 | } |
940 | if (run > 0) |
941 | { |
942 | thiscoord += 3; |
943 | for (k = 0; k < run; k += 3) |
944 | { |
945 | receiveints(buf, 3, smallidx, sizesmall, thiscoord); |
946 | i++; |
947 | thiscoord[0] += prevcoord[0] - smallnum; |
948 | thiscoord[1] += prevcoord[1] - smallnum; |
949 | thiscoord[2] += prevcoord[2] - smallnum; |
950 | if (k == 0) |
951 | { |
952 | /* interchange first with second atom for better |
953 | * compression of water molecules |
954 | */ |
955 | tmp = thiscoord[0]; thiscoord[0] = prevcoord[0]; |
956 | prevcoord[0] = tmp; |
957 | tmp = thiscoord[1]; thiscoord[1] = prevcoord[1]; |
958 | prevcoord[1] = tmp; |
959 | tmp = thiscoord[2]; thiscoord[2] = prevcoord[2]; |
960 | prevcoord[2] = tmp; |
961 | *lfp++ = prevcoord[0] * inv_precision; |
962 | *lfp++ = prevcoord[1] * inv_precision; |
963 | *lfp++ = prevcoord[2] * inv_precision; |
964 | } |
965 | else |
966 | { |
967 | prevcoord[0] = thiscoord[0]; |
968 | prevcoord[1] = thiscoord[1]; |
969 | prevcoord[2] = thiscoord[2]; |
970 | } |
971 | *lfp++ = thiscoord[0] * inv_precision; |
972 | *lfp++ = thiscoord[1] * inv_precision; |
973 | *lfp++ = thiscoord[2] * inv_precision; |
974 | } |
975 | } |
976 | else |
977 | { |
978 | *lfp++ = thiscoord[0] * inv_precision; |
979 | *lfp++ = thiscoord[1] * inv_precision; |
980 | *lfp++ = thiscoord[2] * inv_precision; |
981 | } |
982 | smallidx += is_smaller; |
983 | if (is_smaller < 0) |
984 | { |
985 | smallnum = smaller; |
986 | if (smallidx > FIRSTIDX9) |
987 | { |
988 | smaller = magicints[smallidx - 1] /2; |
989 | } |
990 | else |
991 | { |
992 | smaller = 0; |
993 | } |
994 | } |
995 | else if (is_smaller > 0) |
996 | { |
997 | smaller = smallnum; |
998 | smallnum = magicints[smallidx] / 2; |
999 | } |
1000 | sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx]; |
1001 | } |
1002 | } |
1003 | if (we_should_free) |
1004 | { |
1005 | free(ip); |
1006 | free(buf); |
1007 | } |
1008 | return 1; |
1009 | } |
1010 | |
1011 | |
1012 | |
1013 | /****************************************************************** |
1014 | |
1015 | XTC files have a relatively simple structure. |
1016 | They have a header of 16 bytes and the rest are |
1017 | the compressed coordinates of the files. Due to the |
1018 | compression 00 is not present in the coordinates. |
1019 | The first 4 bytes of the header are the magic number |
1020 | 1995 (0x000007CB). If we find this number we are guaranteed |
1021 | to be in the header, due to the presence of so many zeros. |
1022 | The second 4 bytes are the number of atoms in the frame, and is |
1023 | assumed to be constant. The third 4 bytes are the frame number. |
1024 | The last 4 bytes are a floating point representation of the time. |
1025 | |
1026 | ********************************************************************/ |
1027 | |
1028 | /* Must match definition in xtcio.c */ |
1029 | #ifndef XTC_MAGIC1995 |
1030 | #define XTC_MAGIC1995 1995 |
1031 | #endif |
1032 | |
1033 | static const int header_size = 16; |
1034 | |
1035 | /* Check if we are at the header start. |
1036 | At the same time it will also read 1 int */ |
1037 | static int xtc_at_header_start(FILE *fp, XDR *xdrs, |
1038 | int natoms, int * timestep, float * time) |
1039 | { |
1040 | int i_inp[3]; |
1041 | float f_inp[10]; |
1042 | int i; |
1043 | gmx_off_t off; |
1044 | |
1045 | |
1046 | if ((off = gmx_ftell(fp)) < 0) |
1047 | { |
1048 | return -1; |
1049 | } |
1050 | /* read magic natoms and timestep */ |
1051 | for (i = 0; i < 3; i++) |
1052 | { |
1053 | if (!xdr_int(xdrs, &(i_inp[i]))) |
1054 | { |
1055 | gmx_fseek(fp, off+XDR_INT_SIZE4, SEEK_SET0); |
1056 | return -1; |
1057 | } |
1058 | } |
1059 | /* quick return */ |
1060 | if (i_inp[0] != XTC_MAGIC1995) |
1061 | { |
1062 | if (gmx_fseek(fp, off+XDR_INT_SIZE4, SEEK_SET0)) |
1063 | { |
1064 | return -1; |
1065 | } |
1066 | return 0; |
1067 | } |
1068 | /* read time and box */ |
1069 | for (i = 0; i < 10; i++) |
1070 | { |
1071 | if (!xdr_float(xdrs, &(f_inp[i]))) |
1072 | { |
1073 | gmx_fseek(fp, off+XDR_INT_SIZE4, SEEK_SET0); |
1074 | return -1; |
1075 | } |
1076 | } |
1077 | /* Make a rigourous check to see if we are in the beggining of a header |
1078 | Hopefully there are no ambiguous cases */ |
1079 | /* This check makes use of the fact that the box matrix has 3 zeroes on the upper |
1080 | right triangle and that the first element must be nonzero unless the entire matrix is zero |
1081 | */ |
1082 | if (i_inp[1] == natoms && |
1083 | ((f_inp[1] != 0 && f_inp[6] == 0) || |
1084 | (f_inp[1] == 0 && f_inp[5] == 0 && f_inp[9] == 0))) |
1085 | { |
1086 | if (gmx_fseek(fp, off+XDR_INT_SIZE4, SEEK_SET0)) |
1087 | { |
1088 | return -1; |
1089 | } |
1090 | *time = f_inp[0]; |
1091 | *timestep = i_inp[2]; |
1092 | return 1; |
1093 | } |
1094 | if (gmx_fseek(fp, off+XDR_INT_SIZE4, SEEK_SET0)) |
1095 | { |
1096 | return -1; |
1097 | } |
1098 | return 0; |
1099 | } |
1100 | |
1101 | static int |
1102 | xtc_get_next_frame_number(FILE *fp, XDR *xdrs, int natoms) |
1103 | { |
1104 | gmx_off_t off; |
1105 | int step; |
1106 | float time; |
1107 | int ret; |
1108 | |
1109 | if ((off = gmx_ftell(fp)) < 0) |
1110 | { |
1111 | return -1; |
1112 | } |
1113 | |
1114 | /* read one int just to make sure we dont read this frame but the next */ |
1115 | xdr_int(xdrs, &step); |
1116 | while (1) |
1117 | { |
1118 | ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time); |
1119 | if (ret == 1) |
1120 | { |
1121 | if (gmx_fseek(fp, off, SEEK_SET0)) |
1122 | { |
1123 | return -1; |
1124 | } |
1125 | return step; |
1126 | } |
1127 | else if (ret == -1) |
1128 | { |
1129 | if (gmx_fseek(fp, off, SEEK_SET0)) |
1130 | { |
1131 | return -1; |
1132 | } |
1133 | } |
1134 | } |
1135 | return -1; |
1136 | } |
1137 | |
1138 | |
1139 | static float xtc_get_next_frame_time(FILE *fp, XDR *xdrs, int natoms, |
1140 | gmx_bool * bOK) |
1141 | { |
1142 | gmx_off_t off; |
1143 | float time; |
1144 | int step; |
1145 | int ret; |
1146 | *bOK = 0; |
1147 | |
1148 | if ((off = gmx_ftell(fp)) < 0) |
1149 | { |
1150 | return -1; |
1151 | } |
1152 | /* read one int just to make sure we dont read this frame but the next */ |
1153 | xdr_int(xdrs, &step); |
1154 | while (1) |
1155 | { |
1156 | ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time); |
1157 | if (ret == 1) |
1158 | { |
1159 | *bOK = 1; |
1160 | if (gmx_fseek(fp, off, SEEK_SET0)) |
1161 | { |
1162 | *bOK = 0; |
1163 | return -1; |
1164 | } |
1165 | return time; |
1166 | } |
1167 | else if (ret == -1) |
1168 | { |
1169 | if (gmx_fseek(fp, off, SEEK_SET0)) |
1170 | { |
1171 | return -1; |
1172 | } |
1173 | return -1; |
1174 | } |
1175 | } |
1176 | return -1; |
1177 | } |
1178 | |
1179 | |
1180 | static float |
1181 | xtc_get_current_frame_time(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK) |
1182 | { |
1183 | gmx_off_t off; |
1184 | int step; |
1185 | float time; |
1186 | int ret; |
1187 | *bOK = 0; |
1188 | |
1189 | if ((off = gmx_ftell(fp)) < 0) |
1190 | { |
1191 | return -1; |
1192 | } |
1193 | |
1194 | while (1) |
1195 | { |
1196 | ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time); |
1197 | if (ret == 1) |
1198 | { |
1199 | *bOK = 1; |
1200 | if (gmx_fseek(fp, off, SEEK_SET0)) |
1201 | { |
1202 | *bOK = 0; |
1203 | return -1; |
1204 | } |
1205 | return time; |
1206 | } |
1207 | else if (ret == -1) |
1208 | { |
1209 | if (gmx_fseek(fp, off, SEEK_SET0)) |
1210 | { |
1211 | return -1; |
1212 | } |
1213 | return -1; |
1214 | } |
1215 | else if (ret == 0) |
1216 | { |
1217 | /*Go back.*/ |
1218 | if (gmx_fseek(fp, -2*XDR_INT_SIZE4, SEEK_CUR1)) |
1219 | { |
1220 | return -1; |
1221 | } |
1222 | } |
1223 | } |
1224 | return -1; |
1225 | } |
1226 | |
1227 | /* Currently not used, just for completeness */ |
1228 | static int |
1229 | xtc_get_current_frame_number(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK) |
1230 | { |
1231 | gmx_off_t off; |
1232 | int ret; |
1233 | int step; |
1234 | float time; |
1235 | *bOK = 0; |
1236 | |
1237 | if ((off = gmx_ftell(fp)) < 0) |
1238 | { |
1239 | return -1; |
1240 | } |
1241 | |
1242 | |
1243 | while (1) |
1244 | { |
1245 | ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time); |
1246 | if (ret == 1) |
1247 | { |
1248 | *bOK = 1; |
1249 | if (gmx_fseek(fp, off, SEEK_SET0)) |
1250 | { |
1251 | *bOK = 0; |
1252 | return -1; |
1253 | } |
1254 | return step; |
1255 | } |
1256 | else if (ret == -1) |
1257 | { |
1258 | if (gmx_fseek(fp, off, SEEK_SET0)) |
1259 | { |
1260 | return -1; |
1261 | } |
1262 | return -1; |
1263 | |
1264 | } |
1265 | else if (ret == 0) |
1266 | { |
1267 | /*Go back.*/ |
1268 | if (gmx_fseek(fp, -2*XDR_INT_SIZE4, SEEK_CUR1)) |
1269 | { |
1270 | return -1; |
1271 | } |
1272 | } |
1273 | } |
1274 | return -1; |
1275 | } |
1276 | |
1277 | |
1278 | |
1279 | static gmx_off_t xtc_get_next_frame_start(FILE *fp, XDR *xdrs, int natoms) |
1280 | { |
1281 | int inp; |
1282 | gmx_off_t res; |
1283 | int ret; |
1284 | int step; |
1285 | float time; |
1286 | /* read one int just to make sure we dont read this frame but the next */ |
1287 | xdr_int(xdrs, &step); |
1288 | while (1) |
1289 | { |
1290 | ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time); |
1291 | if (ret == 1) |
1292 | { |
1293 | if ((res = gmx_ftell(fp)) >= 0) |
1294 | { |
1295 | return res - XDR_INT_SIZE4; |
1296 | } |
1297 | else |
1298 | { |
1299 | return res; |
1300 | } |
1301 | } |
1302 | else if (ret == -1) |
1303 | { |
1304 | return -1; |
1305 | } |
1306 | } |
1307 | return -1; |
1308 | } |
1309 | |
1310 | |
1311 | static |
1312 | float |
1313 | xdr_xtc_estimate_dt(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK) |
1314 | { |
1315 | float res; |
1316 | float tinit; |
1317 | gmx_off_t off; |
1318 | |
1319 | *bOK = 0; |
1320 | if ((off = gmx_ftell(fp)) < 0) |
1321 | { |
1322 | return -1; |
1323 | } |
1324 | |
1325 | tinit = xtc_get_current_frame_time(fp, xdrs, natoms, bOK); |
1326 | |
1327 | if (!(*bOK)) |
1328 | { |
1329 | return -1; |
1330 | } |
1331 | |
1332 | res = xtc_get_next_frame_time(fp, xdrs, natoms, bOK); |
1333 | |
1334 | if (!(*bOK)) |
1335 | { |
1336 | return -1; |
1337 | } |
1338 | |
1339 | res -= tinit; |
1340 | if (0 != gmx_fseek(fp, off, SEEK_SET0)) |
1341 | { |
1342 | *bOK = 0; |
1343 | return -1; |
1344 | } |
1345 | return res; |
1346 | } |
1347 | |
1348 | |
1349 | int |
1350 | xdr_xtc_seek_frame(int frame, FILE *fp, XDR *xdrs, int natoms) |
1351 | { |
1352 | gmx_off_t low = 0; |
1353 | gmx_off_t high, pos; |
1354 | |
1355 | |
1356 | /* round to 4 bytes */ |
1357 | int fr; |
1358 | gmx_off_t offset; |
1359 | if (gmx_fseek(fp, 0, SEEK_END2)) |
1360 | { |
1361 | return -1; |
1362 | } |
1363 | |
1364 | if ((high = gmx_ftell(fp)) < 0) |
1365 | { |
1366 | return -1; |
1367 | } |
1368 | |
1369 | /* round to 4 bytes */ |
1370 | high /= XDR_INT_SIZE4; |
1371 | high *= XDR_INT_SIZE4; |
1372 | offset = ((high/2)/XDR_INT_SIZE4)*XDR_INT_SIZE4; |
1373 | |
1374 | if (gmx_fseek(fp, offset, SEEK_SET0)) |
1375 | { |
1376 | return -1; |
1377 | } |
1378 | |
1379 | while (1) |
1380 | { |
1381 | fr = xtc_get_next_frame_number(fp, xdrs, natoms); |
1382 | if (fr < 0) |
1383 | { |
1384 | return -1; |
1385 | } |
1386 | if (fr != frame && abs(low-high) > header_size) |
1387 | { |
1388 | if (fr < frame) |
1389 | { |
1390 | low = offset; |
1391 | } |
1392 | else |
1393 | { |
1394 | high = offset; |
1395 | } |
1396 | /* round to 4 bytes */ |
1397 | offset = (((high+low)/2)/4)*4; |
1398 | |
1399 | if (gmx_fseek(fp, offset, SEEK_SET0)) |
1400 | { |
1401 | return -1; |
1402 | } |
1403 | } |
1404 | else |
1405 | { |
1406 | break; |
1407 | } |
1408 | } |
1409 | if (offset <= header_size) |
1410 | { |
1411 | offset = low; |
1412 | } |
1413 | |
1414 | if (gmx_fseek(fp, offset, SEEK_SET0)) |
1415 | { |
1416 | return -1; |
1417 | } |
1418 | |
1419 | if ((pos = xtc_get_next_frame_start(fp, xdrs, natoms)) < 0) |
1420 | { |
1421 | /* we probably hit an end of file */ |
1422 | return -1; |
1423 | } |
1424 | |
1425 | if (gmx_fseek(fp, pos, SEEK_SET0)) |
1426 | { |
1427 | return -1; |
1428 | } |
1429 | |
1430 | return 0; |
1431 | } |
1432 | |
1433 | |
1434 | |
1435 | int xdr_xtc_seek_time(real time, FILE *fp, XDR *xdrs, int natoms, gmx_bool bSeekForwardOnly) |
1436 | { |
1437 | float t; |
1438 | float dt; |
1439 | gmx_bool bOK = FALSE0; |
1440 | gmx_off_t low = 0; |
1441 | gmx_off_t high, offset, pos; |
1442 | int res; |
1443 | int dt_sign = 0; |
1444 | |
1445 | if (bSeekForwardOnly) |
1446 | { |
1447 | low = gmx_ftell(fp); |
1448 | } |
1449 | if (gmx_fseek(fp, 0, SEEK_END2)) |
1450 | { |
1451 | return -1; |
1452 | } |
1453 | |
1454 | if ((high = gmx_ftell(fp)) < 0) |
1455 | { |
1456 | return -1; |
1457 | } |
1458 | /* round to int */ |
1459 | high /= XDR_INT_SIZE4; |
1460 | high *= XDR_INT_SIZE4; |
1461 | offset = (((high-low) / 2) / XDR_INT_SIZE4) * XDR_INT_SIZE4; |
1462 | |
1463 | if (gmx_fseek(fp, offset, SEEK_SET0)) |
1464 | { |
1465 | return -1; |
1466 | } |
1467 | |
1468 | |
1469 | /* |
1470 | * No need to call xdr_xtc_estimate_dt here - since xdr_xtc_estimate_dt is called first thing in the loop |
1471 | dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK); |
1472 | |
1473 | if (!bOK) |
1474 | { |
1475 | return -1; |
1476 | } |
1477 | */ |
1478 | |
1479 | while (1) |
1480 | { |
1481 | dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK); |
1482 | if (!bOK) |
1483 | { |
1484 | return -1; |
1485 | } |
1486 | else |
1487 | { |
1488 | if (dt > 0) |
1489 | { |
1490 | if (dt_sign == -1) |
1491 | { |
1492 | /* Found a place in the trajectory that has positive time step while |
1493 | other has negative time step */ |
1494 | return -2; |
1495 | } |
1496 | dt_sign = 1; |
1497 | } |
1498 | else if (dt < 0) |
1499 | { |
1500 | if (dt_sign == 1) |
1501 | { |
1502 | /* Found a place in the trajectory that has positive time step while |
1503 | other has negative time step */ |
1504 | return -2; |
1505 | } |
1506 | dt_sign = -1; |
1507 | } |
1508 | } |
1509 | t = xtc_get_next_frame_time(fp, xdrs, natoms, &bOK); |
1510 | if (!bOK) |
1511 | { |
1512 | return -1; |
1513 | } |
1514 | |
1515 | /* If we are before the target time and the time step is positive or 0, or we have |
1516 | after the target time and the time step is negative, or the difference between |
1517 | the current time and the target time is bigger than dt and above all the distance between high |
1518 | and low is bigger than 1 frame, then do another step of binary search. Otherwise stop and check |
1519 | if we reached the solution */ |
1520 | if ((((t < time && dt_sign >= 0) || (t > time && dt_sign == -1)) || ((t |
1521 | - time) >= dt && dt_sign >= 0) |
1522 | || ((time - t) >= -dt && dt_sign < 0)) && (abs(low - high) |
1523 | > header_size)) |
1524 | { |
1525 | if (dt >= 0 && dt_sign != -1) |
1526 | { |
1527 | if (t < time) |
1528 | { |
1529 | low = offset; |
1530 | } |
1531 | else |
1532 | { |
1533 | high = offset; |
1534 | } |
1535 | } |
1536 | else if (dt <= 0 && dt_sign == -1) |
1537 | { |
1538 | if (t >= time) |
1539 | { |
1540 | low = offset; |
1541 | } |
1542 | else |
1543 | { |
1544 | high = offset; |
1545 | } |
1546 | } |
1547 | else |
1548 | { |
1549 | /* We should never reach here */ |
1550 | return -1; |
1551 | } |
1552 | /* round to 4 bytes and subtract header*/ |
1553 | offset = (((high + low) / 2) / XDR_INT_SIZE4) * XDR_INT_SIZE4; |
1554 | if (gmx_fseek(fp, offset, SEEK_SET0)) |
1555 | { |
1556 | return -1; |
1557 | } |
1558 | } |
1559 | else |
1560 | { |
1561 | if (abs(low - high) <= header_size) |
1562 | { |
1563 | break; |
1564 | } |
1565 | /* re-estimate dt */ |
1566 | if (xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK) != dt) |
1567 | { |
1568 | if (bOK) |
1569 | { |
1570 | dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK); |
1571 | } |
1572 | } |
1573 | if (t >= time && t - time < dt) |
1574 | { |
1575 | break; |
1576 | } |
1577 | } |
1578 | } |
1579 | |
1580 | if (offset <= header_size) |
1581 | { |
1582 | offset = low; |
1583 | } |
1584 | |
1585 | gmx_fseek(fp, offset, SEEK_SET0); |
1586 | |
1587 | if ((pos = xtc_get_next_frame_start(fp, xdrs, natoms)) < 0) |
1588 | { |
1589 | return -1; |
1590 | } |
1591 | |
1592 | if (gmx_fseek(fp, pos, SEEK_SET0)) |
1593 | { |
1594 | return -1; |
1595 | } |
1596 | return 0; |
1597 | } |
1598 | |
1599 | float |
1600 | xdr_xtc_get_last_frame_time(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK) |
1601 | { |
1602 | float time; |
1603 | gmx_off_t off; |
1604 | int res; |
1605 | *bOK = 1; |
1606 | off = gmx_ftell(fp); |
1607 | if (off < 0) |
1608 | { |
1609 | *bOK = 0; |
1610 | return -1; |
1611 | } |
1612 | |
1613 | if ( (res = gmx_fseek(fp, -3*XDR_INT_SIZE4, SEEK_END2)) != 0) |
1614 | { |
1615 | *bOK = 0; |
1616 | return -1; |
1617 | } |
1618 | |
1619 | time = xtc_get_current_frame_time(fp, xdrs, natoms, bOK); |
1620 | if (!(*bOK)) |
1621 | { |
1622 | return -1; |
1623 | } |
1624 | |
1625 | if ( (res = gmx_fseek(fp, off, SEEK_SET0)) != 0) |
1626 | { |
1627 | *bOK = 0; |
1628 | return -1; |
1629 | } |
1630 | return time; |
1631 | } |
1632 | |
1633 | |
1634 | int |
1635 | xdr_xtc_get_last_frame_number(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK) |
1636 | { |
1637 | int frame; |
1638 | gmx_off_t off; |
1639 | int res; |
1640 | *bOK = 1; |
1641 | |
1642 | if ((off = gmx_ftell(fp)) < 0) |
1643 | { |
1644 | *bOK = 0; |
1645 | return -1; |
1646 | } |
1647 | |
1648 | |
1649 | if (gmx_fseek(fp, -3*XDR_INT_SIZE4, SEEK_END2)) |
1650 | { |
1651 | *bOK = 0; |
1652 | return -1; |
1653 | } |
1654 | |
1655 | frame = xtc_get_current_frame_number(fp, xdrs, natoms, bOK); |
1656 | if (!bOK) |
1657 | { |
1658 | return -1; |
1659 | } |
1660 | |
1661 | |
1662 | if (gmx_fseek(fp, off, SEEK_SET0)) |
1663 | { |
1664 | *bOK = 0; |
1665 | return -1; |
1666 | } |
1667 | |
1668 | return frame; |
1669 | } |