Context Navigation

source: thirdparty/SZ/sz/src/szd_double.c @ 9ee2ce3

Revision 9ee2ce3, 153.6 KB checked in by Hal Finkel <hfinkel@…>, 6 years ago (diff)
importing new SZ files
Property mode set to `100644`

Line
1	/**
2	* @file szd_double.c
3	* @author Sheng Di and Dingwen Tao
4	* @date Aug, 2016
5	* @brief
6	* (C) 2016 by Mathematics and Computer Science (MCS), Argonne National Laboratory.
7	* See COPYRIGHT in top-level directory.
8	*/
9
10	#include <stdlib.h>
11	#include <stdio.h>
12	#include <string.h>
13	#include "szd_double.h"
14	#include "TightDataPointStorageD.h"
15	#include "sz.h"
16	#include "Huffman.h"
17	#include "szd_double_pwr.h"
18	#include "szd_double_ts.h"
19	#include "utility.h"
20
21	int SZ_decompress_args_double(double** newData, size_t r5, size_t r4, size_t r3, size_t r2, size_t r1, unsigned char* cmpBytes, size_t cmpSize)
22	{
23	int status = SZ_SCES;
24	size_t dataLength = computeDataLength(r5,r4,r3,r2,r1);
25
26	//unsigned char* tmpBytes;
27	size_t targetUncompressSize = dataLength <<3; //i.e., *8
28	//tmpSize must be "much" smaller than dataLength
29	size_t i, tmpSize = 12+MetaDataByteLength+exe_params->SZ_SIZE_TYPE;
30	unsigned char* szTmpBytes;
31	if(cmpSize!=12+4+MetaDataByteLength && cmpSize!=12+8+MetaDataByteLength)
32	{
33	confparams_dec->losslessCompressor = is_lossless_compressed_data(cmpBytes, cmpSize);
34	if(confparams_dec->szMode!=SZ_TEMPORAL_COMPRESSION)
35	{
36	if(confparams_dec->losslessCompressor!=-1)
37	confparams_dec->szMode = SZ_BEST_COMPRESSION;
38	else
39	confparams_dec->szMode = SZ_BEST_SPEED;
40	}
41	if(confparams_dec->szMode==SZ_BEST_SPEED)
42	{
43	tmpSize = cmpSize;
44	szTmpBytes = cmpBytes;
45	}
46	else if(confparams_dec->szMode==SZ_BEST_COMPRESSION \|\| confparams_dec->szMode==SZ_DEFAULT_COMPRESSION \|\| confparams_dec->szMode==SZ_TEMPORAL_COMPRESSION)
47	{
48	if(targetUncompressSize<MIN_ZLIB_DEC_ALLOMEM_BYTES) //Considering the minimum size
49	targetUncompressSize = MIN_ZLIB_DEC_ALLOMEM_BYTES;
50	tmpSize = sz_lossless_decompress(confparams_dec->losslessCompressor, cmpBytes, (unsigned long)cmpSize, &szTmpBytes, (unsigned long)targetUncompressSize+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE);
51	//szTmpBytes = (unsigned char)malloc(sizeof(unsigned char)tmpSize);
52	//memcpy(szTmpBytes, tmpBytes, tmpSize);
53	//free(tmpBytes); //release useless memory
54	}
55	else
56	{
57	printf("Wrong value of confparams_dec->szMode in the double compressed bytes.\n");
58	status = SZ_MERR;
59	return status;
60	}
61	}
62	else
63	szTmpBytes = cmpBytes;
64	//TODO: convert szTmpBytes to double array.
65	TightDataPointStorageD* tdps;
66	int errBoundMode = new_TightDataPointStorageD_fromFlatBytes(&tdps, szTmpBytes, tmpSize);
67
68	int dim = computeDimension(r5,r4,r3,r2,r1);
69	int doubleSize = sizeof(double);
70	if(tdps->isLossless)
71	{
72	newData = (double)malloc(doubleSize*dataLength);
73	if(sysEndianType==BIG_ENDIAN_SYSTEM)
74	{
75	memcpy(newData, szTmpBytes+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE, dataLengthdoubleSize);
76	}
77	else
78	{
79	unsigned char* p = szTmpBytes+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE;
80	for(i=0;i<dataLength;i++,p+=doubleSize)
81	(*newData)[i] = bytesToDouble(p);
82	}
83	}
84	else
85	{
86	if(tdps->raBytes_size > 0) //v2.0
87	{
88	if (dim == 1)
89	getSnapshotData_double_1D(newData,r1,tdps, errBoundMode);
90	else if(dim == 2)
91	decompressDataSeries_double_2D_nonblocked_with_blocked_regression(newData, r2, r1, tdps->raBytes);
92	else if(dim == 3)
93	decompressDataSeries_double_3D_nonblocked_with_blocked_regression(newData, r3, r2, r1, tdps->raBytes);
94	else if(dim == 4)
95	decompressDataSeries_double_3D_nonblocked_with_blocked_regression(newData, r4*r3, r2, r1, tdps->raBytes);
96	else
97	{
98	printf("Error: currently support only at most 4 dimensions!\n");
99	status = SZ_DERR;
100	}
101	}
102	else //1.4.13
103	{
104	if (dim == 1)
105	getSnapshotData_double_1D(newData,r1,tdps, errBoundMode);
106	else
107	if (dim == 2)
108	getSnapshotData_double_2D(newData,r2,r1,tdps, errBoundMode);
109	else
110	if (dim == 3)
111	getSnapshotData_double_3D(newData,r3,r2,r1,tdps, errBoundMode);
112	else
113	if (dim == 4)
114	getSnapshotData_double_4D(newData,r4,r3,r2,r1,tdps, errBoundMode);
115	else
116	{
117	printf("Error: currently support only at most 4 dimensions!\n");
118	status = SZ_DERR;
119	}
120	}
121	}
122
123	free_TightDataPointStorageD2(tdps);
124	if(confparams_dec->szMode!=SZ_BEST_SPEED && cmpSize!=12+MetaDataByteLength+exe_params->SZ_SIZE_TYPE)
125	free(szTmpBytes);
126	return status;
127	}
128
129	void decompressDataSeries_double_1D(double** data, size_t dataSeriesLength, TightDataPointStorageD* tdps)
130	{
131	updateQuantizationInfo(tdps->intervals);
132	size_t i, j, k = 0, p = 0, l = 0; // k is to track the location of residual_bit
133	// in resiMidBits, p is to track the
134	// byte_index of resiMidBits, l is for
135	// leadNum
136	unsigned char* leadNum;
137	double interval = tdps->realPrecision*2;
138
139	convertByteArray2IntArray_fast_2b(tdps->exactDataNum, tdps->leadNumArray, tdps->leadNumArray_size, &leadNum);
140	data = (double)malloc(sizeof(double)*dataSeriesLength);
141
142	int* type = (int)malloc(dataSeriesLengthsizeof(int));
143
144	HuffmanTree* huffmanTree = createHuffmanTree(tdps->stateNum);
145	decode_withTree(huffmanTree, tdps->typeArray, dataSeriesLength, type);
146	SZ_ReleaseHuffman(huffmanTree);
147
148	unsigned char preBytes[8];
149	unsigned char curBytes[8];
150
151	memset(preBytes, 0, 8);
152
153	size_t curByteIndex = 0;
154	int reqBytesLength, resiBitsLength, resiBits;
155	unsigned char leadingNum;
156	double medianValue, exactData, predValue;
157
158	reqBytesLength = tdps->reqLength/8;
159	resiBitsLength = tdps->reqLength%8;
160	medianValue = tdps->medianValue;
161
162	int type_;
163	for (i = 0; i < dataSeriesLength; i++) {
164	type_ = type[i];
165	switch (type_) {
166	case 0:
167	// compute resiBits
168	resiBits = 0;
169	if (resiBitsLength != 0) {
170	int kMod8 = k % 8;
171	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
172	if (rightMovSteps > 0) {
173	int code = getRightMovingCode(kMod8, resiBitsLength);
174	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
175	} else if (rightMovSteps < 0) {
176	int code1 = getLeftMovingCode(kMod8);
177	int code2 = getRightMovingCode(kMod8, resiBitsLength);
178	int leftMovSteps = -rightMovSteps;
179	rightMovSteps = 8 - leftMovSteps;
180	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
181	p++;
182	resiBits = resiBits
183	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
184	} else // rightMovSteps == 0
185	{
186	int code = getRightMovingCode(kMod8, resiBitsLength);
187	resiBits = (tdps->residualMidBits[p] & code);
188	p++;
189	}
190	k += resiBitsLength;
191	}
192
193	// recover the exact data
194	memset(curBytes, 0, 8);
195	leadingNum = leadNum[l++];
196	memcpy(curBytes, preBytes, leadingNum);
197	for (j = leadingNum; j < reqBytesLength; j++)
198	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
199	if (resiBitsLength != 0) {
200	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
201	curBytes[reqBytesLength] = resiByte;
202	}
203
204	exactData = bytesToDouble(curBytes);
205	(*data)[i] = exactData + medianValue;
206	memcpy(preBytes,curBytes,8);
207	break;
208	default:
209	//predValue = 2 * (data)[i-1] - (data)[i-2];
210	predValue = (*data)[i-1];
211	(data)[i] = predValue + (type_-exe_params->intvRadius)interval;
212	break;
213	}
214	//printf("%.30G\n",(*data)[i]);
215	}
216
217	#ifdef HAVE_TIMECMPR
218	if(confparams_dec->szMode == SZ_TEMPORAL_COMPRESSION)
219	memcpy(multisteps->hist_data, (data), dataSeriesLengthsizeof(double));
220	#endif
221
222	free(leadNum);
223	free(type);
224	return;
225	}
226
227	void decompressDataSeries_double_2D(double** data, size_t r1, size_t r2, TightDataPointStorageD* tdps)
228	{
229	updateQuantizationInfo(tdps->intervals);
230	//printf("tdps->intervals=%d, exe_params->intvRadius=%d\n", tdps->intervals, exe_params->intvRadius);
231
232	size_t j, k = 0, p = 0, l = 0; // k is to track the location of residual_bit
233	// in resiMidBits, p is to track the
234	// byte_index of resiMidBits, l is for
235	// leadNum
236	size_t dataSeriesLength = r1*r2;
237	// printf ("%d %d\n", r1, r2);
238
239	unsigned char* leadNum;
240	double realPrecision = tdps->realPrecision;
241
242	convertByteArray2IntArray_fast_2b(tdps->exactDataNum, tdps->leadNumArray, tdps->leadNumArray_size, &leadNum);
243
244	data = (double)malloc(sizeof(double)*dataSeriesLength);
245
246	int* type = (int)malloc(dataSeriesLengthsizeof(int));
247
248	HuffmanTree* huffmanTree = createHuffmanTree(tdps->stateNum);
249	decode_withTree(huffmanTree, tdps->typeArray, dataSeriesLength, type);
250	SZ_ReleaseHuffman(huffmanTree);
251
252	unsigned char preBytes[8];
253	unsigned char curBytes[8];
254
255	memset(preBytes, 0, 8);
256
257	size_t curByteIndex = 0;
258	int reqBytesLength, resiBitsLength, resiBits;
259	unsigned char leadingNum;
260	double medianValue, exactData;
261	int type_;
262
263	reqBytesLength = tdps->reqLength/8;
264	resiBitsLength = tdps->reqLength%8;
265	medianValue = tdps->medianValue;
266
267	double pred1D, pred2D;
268	size_t ii, jj;
269
270	/* Process Row-0, data 0 */
271
272	// compute resiBits
273	resiBits = 0;
274	if (resiBitsLength != 0) {
275	int kMod8 = k % 8;
276	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
277	if (rightMovSteps > 0) {
278	int code = getRightMovingCode(kMod8, resiBitsLength);
279	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
280	} else if (rightMovSteps < 0) {
281	int code1 = getLeftMovingCode(kMod8);
282	int code2 = getRightMovingCode(kMod8, resiBitsLength);
283	int leftMovSteps = -rightMovSteps;
284	rightMovSteps = 8 - leftMovSteps;
285	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
286	p++;
287	resiBits = resiBits
288	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
289	} else // rightMovSteps == 0
290	{
291	int code = getRightMovingCode(kMod8, resiBitsLength);
292	resiBits = (tdps->residualMidBits[p] & code);
293	p++;
294	}
295	k += resiBitsLength;
296	}
297
298	// recover the exact data
299	memset(curBytes, 0, 8);
300	leadingNum = leadNum[l++];
301	memcpy(curBytes, preBytes, leadingNum);
302	for (j = leadingNum; j < reqBytesLength; j++)
303	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
304	if (resiBitsLength != 0) {
305	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
306	curBytes[reqBytesLength] = resiByte;
307	}
308
309	exactData = bytesToDouble(curBytes);
310	(*data)[0] = exactData + medianValue;
311	memcpy(preBytes,curBytes,8);
312
313	/* Process Row-0, data 1 */
314	type_ = type[1];
315	if (type_ != 0)
316	{
317	pred1D = (*data)[0];
318	(data)[1] = pred1D + 2 (type_ - exe_params->intvRadius) * realPrecision;
319	}
320	else
321	{
322	// compute resiBits
323	resiBits = 0;
324	if (resiBitsLength != 0) {
325	int kMod8 = k % 8;
326	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
327	if (rightMovSteps > 0) {
328	int code = getRightMovingCode(kMod8, resiBitsLength);
329	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
330	} else if (rightMovSteps < 0) {
331	int code1 = getLeftMovingCode(kMod8);
332	int code2 = getRightMovingCode(kMod8, resiBitsLength);
333	int leftMovSteps = -rightMovSteps;
334	rightMovSteps = 8 - leftMovSteps;
335	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
336	p++;
337	resiBits = resiBits
338	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
339	} else // rightMovSteps == 0
340	{
341	int code = getRightMovingCode(kMod8, resiBitsLength);
342	resiBits = (tdps->residualMidBits[p] & code);
343	p++;
344	}
345	k += resiBitsLength;
346	}
347
348	// recover the exact data
349	memset(curBytes, 0, 8);
350	leadingNum = leadNum[l++];
351	memcpy(curBytes, preBytes, leadingNum);
352	for (j = leadingNum; j < reqBytesLength; j++)
353	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
354	if (resiBitsLength != 0) {
355	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
356	curBytes[reqBytesLength] = resiByte;
357	}
358
359	exactData = bytesToDouble(curBytes);
360	(*data)[1] = exactData + medianValue;
361	memcpy(preBytes,curBytes,8);
362	}
363
364	/* Process Row-0, data 2 --> data r2-1 */
365	for (jj = 2; jj < r2; jj++)
366	{
367	type_ = type[jj];
368	if (type_ != 0)
369	{
370	pred1D = 2(data)[jj-1] - (*data)[jj-2];
371	(data)[jj] = pred1D + 2 (type_ - exe_params->intvRadius) * realPrecision;
372	}
373	else
374	{
375	// compute resiBits
376	resiBits = 0;
377	if (resiBitsLength != 0) {
378	int kMod8 = k % 8;
379	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
380	if (rightMovSteps > 0) {
381	int code = getRightMovingCode(kMod8, resiBitsLength);
382	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
383	} else if (rightMovSteps < 0) {
384	int code1 = getLeftMovingCode(kMod8);
385	int code2 = getRightMovingCode(kMod8, resiBitsLength);
386	int leftMovSteps = -rightMovSteps;
387	rightMovSteps = 8 - leftMovSteps;
388	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
389	p++;
390	resiBits = resiBits
391	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
392	} else // rightMovSteps == 0
393	{
394	int code = getRightMovingCode(kMod8, resiBitsLength);
395	resiBits = (tdps->residualMidBits[p] & code);
396	p++;
397	}
398	k += resiBitsLength;
399	}
400
401	// recover the exact data
402	memset(curBytes, 0, 8);
403	leadingNum = leadNum[l++];
404	memcpy(curBytes, preBytes, leadingNum);
405	for (j = leadingNum; j < reqBytesLength; j++)
406	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
407	if (resiBitsLength != 0) {
408	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
409	curBytes[reqBytesLength] = resiByte;
410	}
411
412	exactData = bytesToDouble(curBytes);
413	(*data)[jj] = exactData + medianValue;
414	memcpy(preBytes,curBytes,8);
415	}
416	}
417
418	size_t index;
419	/* Process Row-1 --> Row-r1-1 */
420	for (ii = 1; ii < r1; ii++)
421	{
422	/* Process row-ii data 0 */
423	index = ii*r2;
424
425	type_ = type[index];
426	if (type_ != 0)
427	{
428	pred1D = (*data)[index-r2];
429	(data)[index] = pred1D + 2 (type_ - exe_params->intvRadius) * realPrecision;
430	}
431	else
432	{
433	// compute resiBits
434	resiBits = 0;
435	if (resiBitsLength != 0) {
436	int kMod8 = k % 8;
437	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
438	if (rightMovSteps > 0) {
439	int code = getRightMovingCode(kMod8, resiBitsLength);
440	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
441	} else if (rightMovSteps < 0) {
442	int code1 = getLeftMovingCode(kMod8);
443	int code2 = getRightMovingCode(kMod8, resiBitsLength);
444	int leftMovSteps = -rightMovSteps;
445	rightMovSteps = 8 - leftMovSteps;
446	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
447	p++;
448	resiBits = resiBits
449	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
450	} else // rightMovSteps == 0
451	{
452	int code = getRightMovingCode(kMod8, resiBitsLength);
453	resiBits = (tdps->residualMidBits[p] & code);
454	p++;
455	}
456	k += resiBitsLength;
457	}
458
459	// recover the exact data
460	memset(curBytes, 0, 8);
461	leadingNum = leadNum[l++];
462	memcpy(curBytes, preBytes, leadingNum);
463	for (j = leadingNum; j < reqBytesLength; j++)
464	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
465	if (resiBitsLength != 0) {
466	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
467	curBytes[reqBytesLength] = resiByte;
468	}
469
470	exactData = bytesToDouble(curBytes);
471	(*data)[index] = exactData + medianValue;
472	memcpy(preBytes,curBytes,8);
473	}
474
475	/* Process row-ii data 1 --> r2-1*/
476	for (jj = 1; jj < r2; jj++)
477	{
478	index = ii*r2+jj;
479	pred2D = (data)[index-1] + (data)[index-r2] - (*data)[index-r2-1];
480
481	type_ = type[index];
482	if (type_ != 0)
483	{
484	(data)[index] = pred2D + 2 (type_ - exe_params->intvRadius) * realPrecision;
485	}
486	else
487	{
488	// compute resiBits
489	resiBits = 0;
490	if (resiBitsLength != 0) {
491	int kMod8 = k % 8;
492	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
493	if (rightMovSteps > 0) {
494	int code = getRightMovingCode(kMod8, resiBitsLength);
495	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
496	} else if (rightMovSteps < 0) {
497	int code1 = getLeftMovingCode(kMod8);
498	int code2 = getRightMovingCode(kMod8, resiBitsLength);
499	int leftMovSteps = -rightMovSteps;
500	rightMovSteps = 8 - leftMovSteps;
501	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
502	p++;
503	resiBits = resiBits
504	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
505	} else // rightMovSteps == 0
506	{
507	int code = getRightMovingCode(kMod8, resiBitsLength);
508	resiBits = (tdps->residualMidBits[p] & code);
509	p++;
510	}
511	k += resiBitsLength;
512	}
513
514	// recover the exact data
515	memset(curBytes, 0, 8);
516	leadingNum = leadNum[l++];
517	memcpy(curBytes, preBytes, leadingNum);
518	for (j = leadingNum; j < reqBytesLength; j++)
519	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
520	if (resiBitsLength != 0) {
521	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
522	curBytes[reqBytesLength] = resiByte;
523	}
524
525	exactData = bytesToDouble(curBytes);
526	(*data)[index] = exactData + medianValue;
527	memcpy(preBytes,curBytes,8);
528	}
529	}
530	}
531
532	#ifdef HAVE_TIMECMPR
533	if(confparams_dec->szMode == SZ_TEMPORAL_COMPRESSION)
534	memcpy(multisteps->hist_data, (data), dataSeriesLengthsizeof(double));
535	#endif
536
537	free(leadNum);
538	free(type);
539	return;
540	}
541
542	void decompressDataSeries_double_3D(double** data, size_t r1, size_t r2, size_t r3, TightDataPointStorageD* tdps)
543	{
544	updateQuantizationInfo(tdps->intervals);
545	size_t j, k = 0, p = 0, l = 0; // k is to track the location of residual_bit
546	// in resiMidBits, p is to track the
547	// byte_index of resiMidBits, l is for
548	// leadNum
549	size_t dataSeriesLength = r1r2r3;
550	size_t r23 = r2*r3;
551	// printf ("%d %d %d\n", r1, r2, r3);
552
553	unsigned char* leadNum;
554	double realPrecision = tdps->realPrecision;
555
556	convertByteArray2IntArray_fast_2b(tdps->exactDataNum, tdps->leadNumArray, tdps->leadNumArray_size, &leadNum);
557
558	data = (double)malloc(sizeof(double)*dataSeriesLength);
559
560	int* type = (int)malloc(dataSeriesLengthsizeof(int));
561
562	HuffmanTree* huffmanTree = createHuffmanTree(tdps->stateNum);
563	decode_withTree(huffmanTree, tdps->typeArray, dataSeriesLength, type);
564	SZ_ReleaseHuffman(huffmanTree);
565
566	unsigned char preBytes[8];
567	unsigned char curBytes[8];
568
569	memset(preBytes, 0, 8);
570
571	size_t curByteIndex = 0;
572	int reqBytesLength, resiBitsLength, resiBits;
573	unsigned char leadingNum;
574	double medianValue, exactData;
575	int type_;
576
577	reqBytesLength = tdps->reqLength/8;
578	resiBitsLength = tdps->reqLength%8;
579	medianValue = tdps->medianValue;
580
581	double pred1D, pred2D, pred3D;
582	size_t ii, jj, kk;
583
584	/////////////////////////// Process layer-0 ///////////////////////////
585	/* Process Row-0 data 0*/
586	// compute resiBits
587	resiBits = 0;
588	if (resiBitsLength != 0) {
589	int kMod8 = k % 8;
590	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
591	if (rightMovSteps > 0) {
592	int code = getRightMovingCode(kMod8, resiBitsLength);
593	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
594	} else if (rightMovSteps < 0) {
595	int code1 = getLeftMovingCode(kMod8);
596	int code2 = getRightMovingCode(kMod8, resiBitsLength);
597	int leftMovSteps = -rightMovSteps;
598	rightMovSteps = 8 - leftMovSteps;
599	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
600	p++;
601	resiBits = resiBits
602	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
603	} else // rightMovSteps == 0
604	{
605	int code = getRightMovingCode(kMod8, resiBitsLength);
606	resiBits = (tdps->residualMidBits[p] & code);
607	p++;
608	}
609	k += resiBitsLength;
610	}
611
612	// recover the exact data
613	memset(curBytes, 0, 8);
614	leadingNum = leadNum[l++];
615	memcpy(curBytes, preBytes, leadingNum);
616	for (j = leadingNum; j < reqBytesLength; j++)
617	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
618	if (resiBitsLength != 0) {
619	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
620	curBytes[reqBytesLength] = resiByte;
621	}
622
623	exactData = bytesToDouble(curBytes);
624	(*data)[0] = exactData + medianValue;
625	memcpy(preBytes,curBytes,8);
626
627	/* Process Row-0, data 1 */
628	pred1D = (*data)[0];
629
630	type_ = type[1];
631	if (type_ != 0)
632	{
633	(data)[1] = pred1D + 2 (type_ - exe_params->intvRadius) * realPrecision;
634	}
635	else
636	{
637	// compute resiBits
638	resiBits = 0;
639	if (resiBitsLength != 0) {
640	int kMod8 = k % 8;
641	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
642	if (rightMovSteps > 0) {
643	int code = getRightMovingCode(kMod8, resiBitsLength);
644	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
645	} else if (rightMovSteps < 0) {
646	int code1 = getLeftMovingCode(kMod8);
647	int code2 = getRightMovingCode(kMod8, resiBitsLength);
648	int leftMovSteps = -rightMovSteps;
649	rightMovSteps = 8 - leftMovSteps;
650	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
651	p++;
652	resiBits = resiBits
653	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
654	} else // rightMovSteps == 0
655	{
656	int code = getRightMovingCode(kMod8, resiBitsLength);
657	resiBits = (tdps->residualMidBits[p] & code);
658	p++;
659	}
660	k += resiBitsLength;
661	}
662
663	// recover the exact data
664	memset(curBytes, 0, 8);
665	leadingNum = leadNum[l++];
666	memcpy(curBytes, preBytes, leadingNum);
667	for (j = leadingNum; j < reqBytesLength; j++)
668	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
669	if (resiBitsLength != 0) {
670	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
671	curBytes[reqBytesLength] = resiByte;
672	}
673
674	exactData = bytesToDouble(curBytes);
675	(*data)[1] = exactData + medianValue;
676	memcpy(preBytes,curBytes,8);
677	}
678
679	/* Process Row-0, data 2 --> data r3-1 */
680	for (jj = 2; jj < r3; jj++)
681	{
682	pred1D = 2(data)[jj-1] - (*data)[jj-2];
683
684	type_ = type[jj];
685	if (type_ != 0)
686	{
687	(data)[jj] = pred1D + 2 (type_ - exe_params->intvRadius) * realPrecision;
688	}
689	else
690	{
691	// compute resiBits
692	resiBits = 0;
693	if (resiBitsLength != 0) {
694	int kMod8 = k % 8;
695	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
696	if (rightMovSteps > 0) {
697	int code = getRightMovingCode(kMod8, resiBitsLength);
698	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
699	} else if (rightMovSteps < 0) {
700	int code1 = getLeftMovingCode(kMod8);
701	int code2 = getRightMovingCode(kMod8, resiBitsLength);
702	int leftMovSteps = -rightMovSteps;
703	rightMovSteps = 8 - leftMovSteps;
704	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
705	p++;
706	resiBits = resiBits
707	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
708	} else // rightMovSteps == 0
709	{
710	int code = getRightMovingCode(kMod8, resiBitsLength);
711	resiBits = (tdps->residualMidBits[p] & code);
712	p++;
713	}
714	k += resiBitsLength;
715	}
716
717	// recover the exact data
718	memset(curBytes, 0, 8);
719	leadingNum = leadNum[l++];
720	memcpy(curBytes, preBytes, leadingNum);
721	for (j = leadingNum; j < reqBytesLength; j++)
722	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
723	if (resiBitsLength != 0) {
724	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
725	curBytes[reqBytesLength] = resiByte;
726	}
727
728	exactData = bytesToDouble(curBytes);
729	(*data)[jj] = exactData + medianValue;
730	memcpy(preBytes,curBytes,8);
731	}
732	}
733
734	size_t index;
735	/* Process Row-1 --> Row-r2-1 */
736	for (ii = 1; ii < r2; ii++)
737	{
738	/* Process row-ii data 0 */
739	index = ii*r3;
740	pred1D = (*data)[index-r3];
741
742	type_ = type[index];
743	if (type_ != 0)
744	{
745	(data)[index] = pred1D + 2 (type_ - exe_params->intvRadius) * realPrecision;
746	}
747	else
748	{
749	// compute resiBits
750	resiBits = 0;
751	if (resiBitsLength != 0) {
752	int kMod8 = k % 8;
753	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
754	if (rightMovSteps > 0) {
755	int code = getRightMovingCode(kMod8, resiBitsLength);
756	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
757	} else if (rightMovSteps < 0) {
758	int code1 = getLeftMovingCode(kMod8);
759	int code2 = getRightMovingCode(kMod8, resiBitsLength);
760	int leftMovSteps = -rightMovSteps;
761	rightMovSteps = 8 - leftMovSteps;
762	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
763	p++;
764	resiBits = resiBits
765	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
766	} else // rightMovSteps == 0
767	{
768	int code = getRightMovingCode(kMod8, resiBitsLength);
769	resiBits = (tdps->residualMidBits[p] & code);
770	p++;
771	}
772	k += resiBitsLength;
773	}
774
775	// recover the exact data
776	memset(curBytes, 0, 8);
777	leadingNum = leadNum[l++];
778	memcpy(curBytes, preBytes, leadingNum);
779	for (j = leadingNum; j < reqBytesLength; j++)
780	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
781	if (resiBitsLength != 0) {
782	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
783	curBytes[reqBytesLength] = resiByte;
784	}
785
786	exactData = bytesToDouble(curBytes);
787	(*data)[index] = exactData + medianValue;
788	memcpy(preBytes,curBytes,8);
789	}
790
791	/* Process row-ii data 1 --> r3-1*/
792	for (jj = 1; jj < r3; jj++)
793	{
794	index = ii*r3+jj;
795	pred2D = (data)[index-1] + (data)[index-r3] - (*data)[index-r3-1];
796
797	type_ = type[index];
798	if (type_ != 0)
799	{
800	(data)[index] = pred2D + 2 (type_ - exe_params->intvRadius) * realPrecision;
801	}
802	else
803	{
804	// compute resiBits
805	resiBits = 0;
806	if (resiBitsLength != 0) {
807	int kMod8 = k % 8;
808	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
809	if (rightMovSteps > 0) {
810	int code = getRightMovingCode(kMod8, resiBitsLength);
811	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
812	} else if (rightMovSteps < 0) {
813	int code1 = getLeftMovingCode(kMod8);
814	int code2 = getRightMovingCode(kMod8, resiBitsLength);
815	int leftMovSteps = -rightMovSteps;
816	rightMovSteps = 8 - leftMovSteps;
817	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
818	p++;
819	resiBits = resiBits
820	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
821	} else // rightMovSteps == 0
822	{
823	int code = getRightMovingCode(kMod8, resiBitsLength);
824	resiBits = (tdps->residualMidBits[p] & code);
825	p++;
826	}
827	k += resiBitsLength;
828	}
829
830	// recover the exact data
831	memset(curBytes, 0, 8);
832	leadingNum = leadNum[l++];
833	memcpy(curBytes, preBytes, leadingNum);
834	for (j = leadingNum; j < reqBytesLength; j++)
835	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
836	if (resiBitsLength != 0) {
837	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
838	curBytes[reqBytesLength] = resiByte;
839	}
840
841	exactData = bytesToDouble(curBytes);
842	(*data)[index] = exactData + medianValue;
843	memcpy(preBytes,curBytes,8);
844	}
845	}
846	}
847
848	/////////////////////////// Process layer-1 --> layer-r1-1 ///////////////////////////
849
850	for (kk = 1; kk < r1; kk++)
851	{
852	/* Process Row-0 data 0*/
853	index = kk*r23;
854	pred1D = (*data)[index-r23];
855
856	type_ = type[index];
857	if (type_ != 0)
858	{
859	(data)[index] = pred1D + 2 (type_ - exe_params->intvRadius) * realPrecision;
860	}
861	else
862	{
863	// compute resiBits
864	resiBits = 0;
865	if (resiBitsLength != 0) {
866	int kMod8 = k % 8;
867	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
868	if (rightMovSteps > 0) {
869	int code = getRightMovingCode(kMod8, resiBitsLength);
870	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
871	} else if (rightMovSteps < 0) {
872	int code1 = getLeftMovingCode(kMod8);
873	int code2 = getRightMovingCode(kMod8, resiBitsLength);
874	int leftMovSteps = -rightMovSteps;
875	rightMovSteps = 8 - leftMovSteps;
876	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
877	p++;
878	resiBits = resiBits
879	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
880	} else // rightMovSteps == 0
881	{
882	int code = getRightMovingCode(kMod8, resiBitsLength);
883	resiBits = (tdps->residualMidBits[p] & code);
884	p++;
885	}
886	k += resiBitsLength;
887	}
888
889	// recover the exact data
890	memset(curBytes, 0, 8);
891	leadingNum = leadNum[l++];
892	memcpy(curBytes, preBytes, leadingNum);
893	for (j = leadingNum; j < reqBytesLength; j++)
894	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
895	if (resiBitsLength != 0) {
896	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
897	curBytes[reqBytesLength] = resiByte;
898	}
899
900	exactData = bytesToDouble(curBytes);
901	(*data)[index] = exactData + medianValue;
902	memcpy(preBytes,curBytes,8);
903	}
904
905	/* Process Row-0 data 1 --> data r3-1 */
906	for (jj = 1; jj < r3; jj++)
907	{
908	index = kk*r23+jj;
909	pred2D = (data)[index-1] + (data)[index-r23] - (*data)[index-r23-1];
910
911	type_ = type[index];
912	if (type_ != 0)
913	{
914	(data)[index] = pred2D + 2 (type_ - exe_params->intvRadius) * realPrecision;
915	}
916	else
917	{
918	// compute resiBits
919	resiBits = 0;
920	if (resiBitsLength != 0) {
921	int kMod8 = k % 8;
922	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
923	if (rightMovSteps > 0) {
924	int code = getRightMovingCode(kMod8, resiBitsLength);
925	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
926	} else if (rightMovSteps < 0) {
927	int code1 = getLeftMovingCode(kMod8);
928	int code2 = getRightMovingCode(kMod8, resiBitsLength);
929	int leftMovSteps = -rightMovSteps;
930	rightMovSteps = 8 - leftMovSteps;
931	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
932	p++;
933	resiBits = resiBits
934	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
935	} else // rightMovSteps == 0
936	{
937	int code = getRightMovingCode(kMod8, resiBitsLength);
938	resiBits = (tdps->residualMidBits[p] & code);
939	p++;
940	}
941	k += resiBitsLength;
942	}
943
944	// recover the exact data
945	memset(curBytes, 0, 8);
946	leadingNum = leadNum[l++];
947	memcpy(curBytes, preBytes, leadingNum);
948	for (j = leadingNum; j < reqBytesLength; j++)
949	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
950	if (resiBitsLength != 0) {
951	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
952	curBytes[reqBytesLength] = resiByte;
953	}
954
955	exactData = bytesToDouble(curBytes);
956	(*data)[index] = exactData + medianValue;
957	memcpy(preBytes,curBytes,8);
958	}
959	}
960
961	/* Process Row-1 --> Row-r2-1 */
962	for (ii = 1; ii < r2; ii++)
963	{
964	/* Process Row-i data 0 */
965	index = kkr23 + iir3;
966	pred2D = (data)[index-r3] + (data)[index-r23] - (*data)[index-r23-r3];
967
968	type_ = type[index];
969	if (type_ != 0)
970	{
971	(data)[index] = pred2D + 2 (type_ - exe_params->intvRadius) * realPrecision;
972	}
973	else
974	{
975	// compute resiBits
976	resiBits = 0;
977	if (resiBitsLength != 0) {
978	int kMod8 = k % 8;
979	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
980	if (rightMovSteps > 0) {
981	int code = getRightMovingCode(kMod8, resiBitsLength);
982	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
983	} else if (rightMovSteps < 0) {
984	int code1 = getLeftMovingCode(kMod8);
985	int code2 = getRightMovingCode(kMod8, resiBitsLength);
986	int leftMovSteps = -rightMovSteps;
987	rightMovSteps = 8 - leftMovSteps;
988	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
989	p++;
990	resiBits = resiBits
991	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
992	} else // rightMovSteps == 0
993	{
994	int code = getRightMovingCode(kMod8, resiBitsLength);
995	resiBits = (tdps->residualMidBits[p] & code);
996	p++;
997	}
998	k += resiBitsLength;
999	}
1000
1001	// recover the exact data
1002	memset(curBytes, 0, 8);
1003	leadingNum = leadNum[l++];
1004	memcpy(curBytes, preBytes, leadingNum);
1005	for (j = leadingNum; j < reqBytesLength; j++)
1006	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1007	if (resiBitsLength != 0) {
1008	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1009	curBytes[reqBytesLength] = resiByte;
1010	}
1011
1012	exactData = bytesToDouble(curBytes);
1013	(*data)[index] = exactData + medianValue;
1014	memcpy(preBytes,curBytes,8);
1015	}
1016
1017	/* Process Row-i data 1 --> data r3-1 */
1018	for (jj = 1; jj < r3; jj++)
1019	{
1020	index = kkr23 + iir3 + jj;
1021	pred3D = (data)[index-1] + (data)[index-r3] + (*data)[index-r23]
1022	- (data)[index-r3-1] - (data)[index-r23-r3] - (data)[index-r23-1] + (data)[index-r23-r3-1];
1023
1024	type_ = type[index];
1025	if (type_ != 0)
1026	{
1027	(data)[index] = pred3D + 2 (type_ - exe_params->intvRadius) * realPrecision;
1028	}
1029	else
1030	{
1031	// compute resiBits
1032	resiBits = 0;
1033	if (resiBitsLength != 0) {
1034	int kMod8 = k % 8;
1035	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1036	if (rightMovSteps > 0) {
1037	int code = getRightMovingCode(kMod8, resiBitsLength);
1038	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1039	} else if (rightMovSteps < 0) {
1040	int code1 = getLeftMovingCode(kMod8);
1041	int code2 = getRightMovingCode(kMod8, resiBitsLength);
1042	int leftMovSteps = -rightMovSteps;
1043	rightMovSteps = 8 - leftMovSteps;
1044	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1045	p++;
1046	resiBits = resiBits
1047	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1048	} else // rightMovSteps == 0
1049	{
1050	int code = getRightMovingCode(kMod8, resiBitsLength);
1051	resiBits = (tdps->residualMidBits[p] & code);
1052	p++;
1053	}
1054	k += resiBitsLength;
1055	}
1056
1057	// recover the exact data
1058	memset(curBytes, 0, 8);
1059	leadingNum = leadNum[l++];
1060	memcpy(curBytes, preBytes, leadingNum);
1061	for (j = leadingNum; j < reqBytesLength; j++)
1062	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1063	if (resiBitsLength != 0) {
1064	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1065	curBytes[reqBytesLength] = resiByte;
1066	}
1067
1068	exactData = bytesToDouble(curBytes);
1069	(*data)[index] = exactData + medianValue;
1070	memcpy(preBytes,curBytes,8);
1071	}
1072	}
1073	}
1074	}
1075
1076	#ifdef HAVE_TIMECMPR
1077	if(confparams_dec->szMode == SZ_TEMPORAL_COMPRESSION)
1078	memcpy(multisteps->hist_data, (data), dataSeriesLengthsizeof(double));
1079	#endif
1080
1081	free(leadNum);
1082	free(type);
1083	return;
1084	}
1085
1086	void decompressDataSeries_double_4D(double** data, size_t r1, size_t r2, size_t r3, size_t r4, TightDataPointStorageD* tdps)
1087	{
1088	updateQuantizationInfo(tdps->intervals);
1089	size_t j, k = 0, p = 0, l = 0; // k is to track the location of residual_bit
1090	// in resiMidBits, p is to track the
1091	// byte_index of resiMidBits, l is for
1092	// leadNum
1093	size_t dataSeriesLength = r1r2r3*r4;
1094	size_t r234 = r2r3r4;
1095	size_t r34 = r3*r4;
1096	// printf ("%d %d %d\n", r1, r2, r3, r4);
1097
1098	unsigned char* leadNum;
1099	double realPrecision = tdps->realPrecision;
1100
1101	convertByteArray2IntArray_fast_2b(tdps->exactDataNum, tdps->leadNumArray, tdps->leadNumArray_size, &leadNum);
1102
1103	data = (double)malloc(sizeof(double)*dataSeriesLength);
1104
1105	int* type = (int)malloc(dataSeriesLengthsizeof(int));
1106
1107	HuffmanTree* huffmanTree = createHuffmanTree(tdps->stateNum);
1108	decode_withTree(huffmanTree, tdps->typeArray, dataSeriesLength, type);
1109	SZ_ReleaseHuffman(huffmanTree);
1110
1111	unsigned char preBytes[8];
1112	unsigned char curBytes[8];
1113
1114	memset(preBytes, 0, 8);
1115
1116	size_t curByteIndex = 0;
1117	int reqBytesLength, resiBitsLength, resiBits;
1118	unsigned char leadingNum;
1119	double medianValue, exactData;
1120	int type_;
1121
1122	reqBytesLength = tdps->reqLength/8;
1123	resiBitsLength = tdps->reqLength%8;
1124	medianValue = tdps->medianValue;
1125
1126	double pred1D, pred2D, pred3D;
1127	size_t ii, jj, kk, ll;
1128	size_t index;
1129
1130	for (ll = 0; ll < r1; ll++)
1131	{
1132
1133	/////////////////////////// Process layer-0 ///////////////////////////
1134	/* Process Row-0 data 0*/
1135	index = ll*r234;
1136
1137	// compute resiBits
1138	resiBits = 0;
1139	if (resiBitsLength != 0) {
1140	int kMod8 = k % 8;
1141	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1142	if (rightMovSteps > 0) {
1143	int code = getRightMovingCode(kMod8, resiBitsLength);
1144	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1145	} else if (rightMovSteps < 0) {
1146	int code1 = getLeftMovingCode(kMod8);
1147	int code2 = getRightMovingCode(kMod8, resiBitsLength);
1148	int leftMovSteps = -rightMovSteps;
1149	rightMovSteps = 8 - leftMovSteps;
1150	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1151	p++;
1152	resiBits = resiBits
1153	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1154	} else // rightMovSteps == 0
1155	{
1156	int code = getRightMovingCode(kMod8, resiBitsLength);
1157	resiBits = (tdps->residualMidBits[p] & code);
1158	p++;
1159	}
1160	k += resiBitsLength;
1161	}
1162
1163	// recover the exact data
1164	memset(curBytes, 0, 8);
1165	leadingNum = leadNum[l++];
1166	memcpy(curBytes, preBytes, leadingNum);
1167	for (j = leadingNum; j < reqBytesLength; j++)
1168	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1169	if (resiBitsLength != 0) {
1170	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1171	curBytes[reqBytesLength] = resiByte;
1172	}
1173
1174	exactData = bytesToDouble(curBytes);
1175	(*data)[index] = exactData + medianValue;
1176	memcpy(preBytes,curBytes,8);
1177
1178	/* Process Row-0, data 1 */
1179	index = ll*r234+1;
1180
1181	pred1D = (*data)[index-1];
1182
1183	type_ = type[index];
1184	if (type_ != 0)
1185	{
1186	(data)[index] = pred1D + 2 (type_ - exe_params->intvRadius) * realPrecision;
1187	}
1188	else
1189	{
1190	// compute resiBits
1191	resiBits = 0;
1192	if (resiBitsLength != 0) {
1193	int kMod8 = k % 8;
1194	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1195	if (rightMovSteps > 0) {
1196	int code = getRightMovingCode(kMod8, resiBitsLength);
1197	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1198	} else if (rightMovSteps < 0) {
1199	int code1 = getLeftMovingCode(kMod8);
1200	int code2 = getRightMovingCode(kMod8, resiBitsLength);
1201	int leftMovSteps = -rightMovSteps;
1202	rightMovSteps = 8 - leftMovSteps;
1203	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1204	p++;
1205	resiBits = resiBits
1206	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1207	} else // rightMovSteps == 0
1208	{
1209	int code = getRightMovingCode(kMod8, resiBitsLength);
1210	resiBits = (tdps->residualMidBits[p] & code);
1211	p++;
1212	}
1213	k += resiBitsLength;
1214	}
1215
1216	// recover the exact data
1217	memset(curBytes, 0, 8);
1218	leadingNum = leadNum[l++];
1219	memcpy(curBytes, preBytes, leadingNum);
1220	for (j = leadingNum; j < reqBytesLength; j++)
1221	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1222	if (resiBitsLength != 0) {
1223	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1224	curBytes[reqBytesLength] = resiByte;
1225	}
1226
1227	exactData = bytesToDouble(curBytes);
1228	(*data)[index] = exactData + medianValue;
1229	memcpy(preBytes,curBytes,8);
1230	}
1231
1232	/* Process Row-0, data 2 --> data r4-1 */
1233	for (jj = 2; jj < r4; jj++)
1234	{
1235	index = ll*r234+jj;
1236
1237	pred1D = 2(data)[index-1] - (*data)[index-2];
1238
1239	type_ = type[index];
1240	if (type_ != 0)
1241	{
1242	(data)[index] = pred1D + 2 (type_ - exe_params->intvRadius) * realPrecision;
1243	}
1244	else
1245	{
1246	// compute resiBits
1247	resiBits = 0;
1248	if (resiBitsLength != 0) {
1249	int kMod8 = k % 8;
1250	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1251	if (rightMovSteps > 0) {
1252	int code = getRightMovingCode(kMod8, resiBitsLength);
1253	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1254	} else if (rightMovSteps < 0) {
1255	int code1 = getLeftMovingCode(kMod8);
1256	int code2 = getRightMovingCode(kMod8, resiBitsLength);
1257	int leftMovSteps = -rightMovSteps;
1258	rightMovSteps = 8 - leftMovSteps;
1259	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1260	p++;
1261	resiBits = resiBits
1262	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1263	} else // rightMovSteps == 0
1264	{
1265	int code = getRightMovingCode(kMod8, resiBitsLength);
1266	resiBits = (tdps->residualMidBits[p] & code);
1267	p++;
1268	}
1269	k += resiBitsLength;
1270	}
1271
1272	// recover the exact data
1273	memset(curBytes, 0, 8);
1274	leadingNum = leadNum[l++];
1275	memcpy(curBytes, preBytes, leadingNum);
1276	for (j = leadingNum; j < reqBytesLength; j++)
1277	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1278	if (resiBitsLength != 0) {
1279	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1280	curBytes[reqBytesLength] = resiByte;
1281	}
1282
1283	exactData = bytesToDouble(curBytes);
1284	(*data)[index] = exactData + medianValue;
1285	memcpy(preBytes,curBytes,8);
1286	}
1287	}
1288
1289	/* Process Row-1 --> Row-r3-1 */
1290	for (ii = 1; ii < r3; ii++)
1291	{
1292	/* Process row-ii data 0 */
1293	index = llr234+iir4;
1294
1295	pred1D = (*data)[index-r4];
1296
1297	type_ = type[index];
1298	if (type_ != 0)
1299	{
1300	(data)[index] = pred1D + 2 (type_ - exe_params->intvRadius) * realPrecision;
1301	}
1302	else
1303	{
1304	// compute resiBits
1305	resiBits = 0;
1306	if (resiBitsLength != 0) {
1307	int kMod8 = k % 8;
1308	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1309	if (rightMovSteps > 0) {
1310	int code = getRightMovingCode(kMod8, resiBitsLength);
1311	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1312	} else if (rightMovSteps < 0) {
1313	int code1 = getLeftMovingCode(kMod8);
1314	int code2 = getRightMovingCode(kMod8, resiBitsLength);
1315	int leftMovSteps = -rightMovSteps;
1316	rightMovSteps = 8 - leftMovSteps;
1317	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1318	p++;
1319	resiBits = resiBits
1320	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1321	} else // rightMovSteps == 0
1322	{
1323	int code = getRightMovingCode(kMod8, resiBitsLength);
1324	resiBits = (tdps->residualMidBits[p] & code);
1325	p++;
1326	}
1327	k += resiBitsLength;
1328	}
1329
1330	// recover the exact data
1331	memset(curBytes, 0, 8);
1332	leadingNum = leadNum[l++];
1333	memcpy(curBytes, preBytes, leadingNum);
1334	for (j = leadingNum; j < reqBytesLength; j++)
1335	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1336	if (resiBitsLength != 0) {
1337	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1338	curBytes[reqBytesLength] = resiByte;
1339	}
1340
1341	exactData = bytesToDouble(curBytes);
1342	(*data)[index] = exactData + medianValue;
1343	memcpy(preBytes,curBytes,8);
1344	}
1345
1346	/* Process row-ii data 1 --> r4-1*/
1347	for (jj = 1; jj < r4; jj++)
1348	{
1349	index = llr234+iir4+jj;
1350
1351	pred2D = (data)[index-1] + (data)[index-r4] - (*data)[index-r4-1];
1352
1353	type_ = type[index];
1354	if (type_ != 0)
1355	{
1356	(data)[index] = pred2D + 2 (type_ - exe_params->intvRadius) * realPrecision;
1357	}
1358	else
1359	{
1360	// compute resiBits
1361	resiBits = 0;
1362	if (resiBitsLength != 0) {
1363	int kMod8 = k % 8;
1364	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1365	if (rightMovSteps > 0) {
1366	int code = getRightMovingCode(kMod8, resiBitsLength);
1367	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1368	} else if (rightMovSteps < 0) {
1369	int code1 = getLeftMovingCode(kMod8);
1370	int code2 = getRightMovingCode(kMod8, resiBitsLength);
1371	int leftMovSteps = -rightMovSteps;
1372	rightMovSteps = 8 - leftMovSteps;
1373	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1374	p++;
1375	resiBits = resiBits
1376	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1377	} else // rightMovSteps == 0
1378	{
1379	int code = getRightMovingCode(kMod8, resiBitsLength);
1380	resiBits = (tdps->residualMidBits[p] & code);
1381	p++;
1382	}
1383	k += resiBitsLength;
1384	}
1385
1386	// recover the exact data
1387	memset(curBytes, 0, 8);
1388	leadingNum = leadNum[l++];
1389	memcpy(curBytes, preBytes, leadingNum);
1390	for (j = leadingNum; j < reqBytesLength; j++)
1391	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1392	if (resiBitsLength != 0) {
1393	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1394	curBytes[reqBytesLength] = resiByte;
1395	}
1396
1397	exactData = bytesToDouble(curBytes);
1398	(*data)[index] = exactData + medianValue;
1399	memcpy(preBytes,curBytes,8);
1400	}
1401	}
1402	}
1403
1404	/////////////////////////// Process layer-1 --> layer-r2-1 ///////////////////////////
1405
1406	for (kk = 1; kk < r2; kk++)
1407	{
1408	/* Process Row-0 data 0*/
1409	index = llr234+kkr34;
1410
1411	pred1D = (*data)[index-r34];
1412
1413	type_ = type[index];
1414	if (type_ != 0)
1415	{
1416	(data)[index] = pred1D + 2 (type_ - exe_params->intvRadius) * realPrecision;
1417	}
1418	else
1419	{
1420	// compute resiBits
1421	resiBits = 0;
1422	if (resiBitsLength != 0) {
1423	int kMod8 = k % 8;
1424	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1425	if (rightMovSteps > 0) {
1426	int code = getRightMovingCode(kMod8, resiBitsLength);
1427	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1428	} else if (rightMovSteps < 0) {
1429	int code1 = getLeftMovingCode(kMod8);
1430	int code2 = getRightMovingCode(kMod8, resiBitsLength);
1431	int leftMovSteps = -rightMovSteps;
1432	rightMovSteps = 8 - leftMovSteps;
1433	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1434	p++;
1435	resiBits = resiBits
1436	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1437	} else // rightMovSteps == 0
1438	{
1439	int code = getRightMovingCode(kMod8, resiBitsLength);
1440	resiBits = (tdps->residualMidBits[p] & code);
1441	p++;
1442	}
1443	k += resiBitsLength;
1444	}
1445
1446	// recover the exact data
1447	memset(curBytes, 0, 8);
1448	leadingNum = leadNum[l++];
1449	memcpy(curBytes, preBytes, leadingNum);
1450	for (j = leadingNum; j < reqBytesLength; j++)
1451	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1452	if (resiBitsLength != 0) {
1453	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1454	curBytes[reqBytesLength] = resiByte;
1455	}
1456
1457	exactData = bytesToDouble(curBytes);
1458	(*data)[index] = exactData + medianValue;
1459	memcpy(preBytes,curBytes,8);
1460	}
1461
1462	/* Process Row-0 data 1 --> data r4-1 */
1463	for (jj = 1; jj < r4; jj++)
1464	{
1465	index = llr234+kkr34+jj;
1466
1467	pred2D = (data)[index-1] + (data)[index-r34] - (*data)[index-r34-1];
1468
1469	type_ = type[index];
1470	if (type_ != 0)
1471	{
1472	(data)[index] = pred2D + 2 (type_ - exe_params->intvRadius) * realPrecision;
1473	}
1474	else
1475	{
1476	// compute resiBits
1477	resiBits = 0;
1478	if (resiBitsLength != 0) {
1479	int kMod8 = k % 8;
1480	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1481	if (rightMovSteps > 0) {
1482	int code = getRightMovingCode(kMod8, resiBitsLength);
1483	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1484	} else if (rightMovSteps < 0) {
1485	int code1 = getLeftMovingCode(kMod8);
1486	int code2 = getRightMovingCode(kMod8, resiBitsLength);
1487	int leftMovSteps = -rightMovSteps;
1488	rightMovSteps = 8 - leftMovSteps;
1489	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1490	p++;
1491	resiBits = resiBits
1492	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1493	} else // rightMovSteps == 0
1494	{
1495	int code = getRightMovingCode(kMod8, resiBitsLength);
1496	resiBits = (tdps->residualMidBits[p] & code);
1497	p++;
1498	}
1499	k += resiBitsLength;
1500	}
1501
1502	// recover the exact data
1503	memset(curBytes, 0, 8);
1504	leadingNum = leadNum[l++];
1505	memcpy(curBytes, preBytes, leadingNum);
1506	for (j = leadingNum; j < reqBytesLength; j++)
1507	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1508	if (resiBitsLength != 0) {
1509	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1510	curBytes[reqBytesLength] = resiByte;
1511	}
1512
1513	exactData = bytesToDouble(curBytes);
1514	(*data)[index] = exactData + medianValue;
1515	memcpy(preBytes,curBytes,8);
1516	}
1517	}
1518
1519	/* Process Row-1 --> Row-r3-1 */
1520	for (ii = 1; ii < r3; ii++)
1521	{
1522	/* Process Row-i data 0 */
1523	index = llr234+kkr34+ii*r4;
1524
1525	pred2D = (data)[index-r4] + (data)[index-r34] - (*data)[index-r34-r4];
1526
1527	type_ = type[index];
1528	if (type_ != 0)
1529	{
1530	(data)[index] = pred2D + 2 (type_ - exe_params->intvRadius) * realPrecision;
1531	}
1532	else
1533	{
1534	// compute resiBits
1535	resiBits = 0;
1536	if (resiBitsLength != 0) {
1537	int kMod8 = k % 8;
1538	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1539	if (rightMovSteps > 0) {
1540	int code = getRightMovingCode(kMod8, resiBitsLength);
1541	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1542	} else if (rightMovSteps < 0) {
1543	int code1 = getLeftMovingCode(kMod8);
1544	int code2 = getRightMovingCode(kMod8, resiBitsLength);
1545	int leftMovSteps = -rightMovSteps;
1546	rightMovSteps = 8 - leftMovSteps;
1547	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1548	p++;
1549	resiBits = resiBits
1550	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1551	} else // rightMovSteps == 0
1552	{
1553	int code = getRightMovingCode(kMod8, resiBitsLength);
1554	resiBits = (tdps->residualMidBits[p] & code);
1555	p++;
1556	}
1557	k += resiBitsLength;
1558	}
1559
1560	// recover the exact data
1561	memset(curBytes, 0, 8);
1562	leadingNum = leadNum[l++];
1563	memcpy(curBytes, preBytes, leadingNum);
1564	for (j = leadingNum; j < reqBytesLength; j++)
1565	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1566	if (resiBitsLength != 0) {
1567	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1568	curBytes[reqBytesLength] = resiByte;
1569	}
1570
1571	exactData = bytesToDouble(curBytes);
1572	(*data)[index] = exactData + medianValue;
1573	memcpy(preBytes,curBytes,8);
1574	}
1575
1576	/* Process Row-i data 1 --> data r4-1 */
1577	for (jj = 1; jj < r4; jj++)
1578	{
1579	index = llr234+kkr34+ii*r4+jj;
1580
1581	pred3D = (data)[index-1] + (data)[index-r4] + (*data)[index-r34]
1582	- (data)[index-r4-1] - (data)[index-r34-r4] - (data)[index-r34-1] + (data)[index-r34-r4-1];
1583
1584	type_ = type[index];
1585	if (type_ != 0)
1586	{
1587	(data)[index] = pred3D + 2 (type_ - exe_params->intvRadius) * realPrecision;
1588	}
1589	else
1590	{
1591	// compute resiBits
1592	resiBits = 0;
1593	if (resiBitsLength != 0) {
1594	int kMod8 = k % 8;
1595	int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1596	if (rightMovSteps > 0) {
1597	int code = getRightMovingCode(kMod8, resiBitsLength);
1598	resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1599	} else if (rightMovSteps < 0) {
1600	int code1 = getLeftMovingCode(kMod8);
1601	int code2 = getRightMovingCode(kMod8, resiBitsLength);
1602	int leftMovSteps = -rightMovSteps;
1603	rightMovSteps = 8 - leftMovSteps;
1604	resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1605	p++;
1606	resiBits = resiBits
1607	\| ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1608	} else // rightMovSteps == 0
1609	{
1610	int code = getRightMovingCode(kMod8, resiBitsLength);
1611	resiBits = (tdps->residualMidBits[p] & code);
1612	p++;
1613	}
1614	k += resiBitsLength;
1615	}
1616
1617	// recover the exact data
1618	memset(curBytes, 0, 8);
1619	leadingNum = leadNum[l++];
1620	memcpy(curBytes, preBytes, leadingNum);
1621	for (j = leadingNum; j < reqBytesLength; j++)
1622	curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1623	if (resiBitsLength != 0) {
1624	unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1625	curBytes[reqBytesLength] = resiByte;
1626	}
1627
1628	exactData = bytesToDouble(curBytes);
1629	(*data)[index] = exactData + medianValue;
1630	memcpy(preBytes,curBytes,8);
1631	}
1632	}
1633	}
1634	}
1635	}
1636
1637	//I didn't implement time-based compression for 4D actually.
1638	//#ifdef HAVE_TIMECMPR
1639	// if(confparams_dec->szMode == SZ_TEMPORAL_COMPRESSION)
1640	// memcpy(multisteps->hist_data, (data), dataSeriesLengthsizeof(double));
1641	//#endif
1642
1643	free(leadNum);
1644	free(type);
1645	return;
1646	}
1647
1648	void getSnapshotData_double_1D(double** data, size_t dataSeriesLength, TightDataPointStorageD* tdps, int errBoundMode)
1649	{
1650	size_t i;
1651	if (tdps->allSameData) {
1652	double value = bytesToDouble(tdps->exactMidBytes);
1653	data = (double)malloc(sizeof(double)*dataSeriesLength);
1654	for (i = 0; i < dataSeriesLength; i++)
1655	(*data)[i] = value;
1656	} else {
1657	if (tdps->rtypeArray == NULL) {
1658	if(errBoundMode < PW_REL)
1659	{
1660	#ifdef HAVE_TIMECMPR
1661	if(confparams_dec->szMode == SZ_TEMPORAL_COMPRESSION)
1662	{
1663	if(multisteps->compressionType == 0) //snapshot
1664	decompressDataSeries_double_1D(data, dataSeriesLength, tdps);
1665	else
1666	decompressDataSeries_double_1D_ts(data, dataSeriesLength, multisteps, tdps);
1667	}
1668	else
1669	#endif
1670	decompressDataSeries_double_1D(data, dataSeriesLength, tdps);
1671	}
1672	else
1673	{
1674	decompressDataSeries_double_1D_pwr_pre_log(data, dataSeriesLength, tdps);
1675	//decompressDataSeries_double_1D_pwrgroup(data, dataSeriesLength, tdps);
1676	}
1677	return;
1678	} else {
1679	data = (double)malloc(sizeof(double)*dataSeriesLength);
1680	// insert the reserved values
1681	//int[] rtypes = TypeManager.convertByteArray2IntArray_fast_1b(
1682	// dataSeriesLength, rtypeArray);
1683	int* rtypes;
1684	int validLength = computeBitNumRequired(dataSeriesLength);
1685	decompressBitArraybySimpleLZ77(&rtypes, tdps->rtypeArray, tdps->rtypeArray_size, dataSeriesLength, validLength);
1686	size_t count = 0;
1687	for (i = 0; i < dataSeriesLength; i++) {
1688	if (rtypes[i] == 1)
1689	(*data)[i] = tdps->reservedValue;
1690	else
1691	count++;
1692	}
1693	// get the decompressed data
1694	double* decmpData;
1695	if(errBoundMode < PW_REL)
1696	decompressDataSeries_double_1D(&decmpData, dataSeriesLength, tdps);
1697	else
1698	//decompressDataSeries_double_1D_pwr(&decmpData, dataSeriesLength, tdps);
1699	decompressDataSeries_double_1D_pwr_pre_log(&decmpData, dataSeriesLength, tdps);
1700	// insert the decompressed data
1701	size_t k = 0;
1702	for (i = 0; i < dataSeriesLength; i++) {
1703	if (rtypes[i] == 0) {
1704	(*data)[i] = decmpData[k++];
1705	}
1706	}
1707	free(decmpData);
1708	free(rtypes);
1709	}
1710	}
1711	}
1712
1713	void getSnapshotData_double_2D(double** data, size_t r1, size_t r2, TightDataPointStorageD* tdps, int errBoundMode)
1714	{
1715	size_t i;
1716	size_t dataSeriesLength = r1*r2;
1717	if (tdps->allSameData) {
1718	double value = bytesToDouble(tdps->exactMidBytes);
1719	data = (double)malloc(sizeof(double)*dataSeriesLength);
1720	for (i = 0; i < dataSeriesLength; i++)
1721	(*data)[i] = value;
1722	} else {
1723	if (tdps->rtypeArray == NULL) {
1724	if(errBoundMode < PW_REL)
1725	{
1726	#ifdef HAVE_TIMECMPR
1727	if(confparams_dec->szMode == SZ_TEMPORAL_COMPRESSION)
1728	{
1729	if(multisteps->compressionType == 0) //snapshot
1730	decompressDataSeries_double_2D(data, r1, r2, tdps);
1731	else
1732	decompressDataSeries_double_1D_ts(data, dataSeriesLength, multisteps, tdps);
1733	}
1734	else
1735	#endif
1736	decompressDataSeries_double_2D(data, r1, r2, tdps);
1737	}
1738	else
1739	//decompressDataSeries_double_2D_pwr(data, r1, r2, tdps);
1740	decompressDataSeries_double_2D_pwr_pre_log(data, r1, r2, tdps);
1741	return;
1742	} else {
1743	data = (double)malloc(sizeof(double)*dataSeriesLength);
1744	// insert the reserved values
1745	//int[] rtypes = TypeManager.convertByteArray2IntArray_fast_1b(
1746	// dataSeriesLength, rtypeArray);
1747	int* rtypes;
1748	int validLength = computeBitNumRequired(dataSeriesLength);
1749	decompressBitArraybySimpleLZ77(&rtypes, tdps->rtypeArray, tdps->rtypeArray_size, dataSeriesLength, validLength);
1750	size_t count = 0;
1751	for (i = 0; i < dataSeriesLength; i++) {
1752	if (rtypes[i] == 1)
1753	(*data)[i] = tdps->reservedValue;
1754	else
1755	count++;
1756	}
1757	// get the decompressed data
1758	double* decmpData;
1759	if(errBoundMode < PW_REL)
1760	decompressDataSeries_double_2D(&decmpData, r1, r2, tdps);
1761	else
1762	//decompressDataSeries_double_2D_pwr(&decmpData, r1, r2, tdps);
1763	decompressDataSeries_double_2D_pwr_pre_log(&decmpData, r1, r2, tdps);
1764	// insert the decompressed data
1765	size_t k = 0;
1766	for (i = 0; i < dataSeriesLength; i++) {
1767	if (rtypes[i] == 0) {
1768	(*data)[i] = decmpData[k++];
1769	}
1770	}
1771	free(decmpData);
1772	free(rtypes);
1773	}
1774	}
1775	}
1776
1777	void getSnapshotData_double_3D(double** data, size_t r1, size_t r2, size_t r3, TightDataPointStorageD* tdps, int errBoundMode)
1778	{
1779	size_t i;
1780	size_t dataSeriesLength = r1r2r3;
1781	if (tdps->allSameData) {
1782	double value = bytesToDouble(tdps->exactMidBytes);
1783	data = (double)malloc(sizeof(double)*dataSeriesLength);
1784	for (i = 0; i < dataSeriesLength; i++)
1785	(*data)[i] = value;
1786	} else {
1787	if (tdps->rtypeArray == NULL) {
1788	if(errBoundMode < PW_REL)
1789	{
1790	#ifdef HAVE_TIMECMPR
1791	if(confparams_dec->szMode == SZ_TEMPORAL_COMPRESSION)
1792	{
1793	if(multisteps->compressionType == 0) //snapshot
1794	decompressDataSeries_double_3D(data, r1, r2, r3, tdps);
1795	else
1796	decompressDataSeries_double_1D_ts(data, dataSeriesLength, multisteps, tdps);
1797	}
1798	else
1799	#endif
1800	decompressDataSeries_double_3D(data, r1, r2, r3, tdps);
1801	}
1802	else
1803	//decompressDataSeries_double_3D_pwr(data, r1, r2, r3, tdps);
1804	decompressDataSeries_double_3D_pwr_pre_log(data, r1, r2, r3, tdps);
1805	return;
1806	} else {
1807	data = (double)malloc(sizeof(double)*dataSeriesLength);
1808	// insert the reserved values
1809	//int[] rtypes = TypeManager.convertByteArray2IntArray_fast_1b(
1810	// dataSeriesLength, rtypeArray);
1811	int* rtypes;
1812	int validLength = computeBitNumRequired(dataSeriesLength);
1813	decompressBitArraybySimpleLZ77(&rtypes, tdps->rtypeArray, tdps->rtypeArray_size, dataSeriesLength, validLength);
1814	size_t count = 0;
1815	for (i = 0; i < dataSeriesLength; i++) {
1816	if (rtypes[i] == 1)
1817	(*data)[i] = tdps->reservedValue;
1818	else
1819	count++;
1820	}
1821	// get the decompressed data
1822	double* decmpData;
1823	if(errBoundMode < PW_REL)
1824	decompressDataSeries_double_3D(&decmpData, r1, r2, r3, tdps);
1825	else
1826	//decompressDataSeries_double_3D_pwr(&decmpData, r1, r2, r3, tdps);
1827	decompressDataSeries_double_3D_pwr_pre_log(&decmpData, r1, r2, r3, tdps);
1828	// insert the decompressed data
1829	size_t k = 0;
1830	for (i = 0; i < dataSeriesLength; i++) {
1831	if (rtypes[i] == 0) {
1832	(*data)[i] = decmpData[k++];
1833	}
1834	}
1835	free(decmpData);
1836	free(rtypes);
1837	}
1838	}
1839	}
1840
1841	void getSnapshotData_double_4D(double** data, size_t r1, size_t r2, size_t r3, size_t r4, TightDataPointStorageD* tdps, int errBoundMode)
1842	{
1843	size_t i;
1844	size_t dataSeriesLength = r1r2r3*r4;
1845	if (tdps->allSameData) {
1846	double value = bytesToDouble(tdps->exactMidBytes);
1847	data = (double)malloc(sizeof(double)*dataSeriesLength);
1848	for (i = 0; i < dataSeriesLength; i++)
1849	(*data)[i] = value;
1850	} else {
1851	if (tdps->rtypeArray == NULL) {
1852	if(errBoundMode < PW_REL)
1853	{
1854	#ifdef HAVE_TIMECMPR
1855	if(confparams_dec->szMode == SZ_TEMPORAL_COMPRESSION)
1856	{
1857	if(multisteps->compressionType == 0)
1858	decompressDataSeries_double_4D(data, r1, r2, r3, r4, tdps);
1859	else
1860	decompressDataSeries_double_1D_ts(data, r1r2r3*r4, multisteps, tdps);
1861	}
1862	else
1863	#endif
1864	decompressDataSeries_double_4D(data, r1, r2, r3, r4, tdps);
1865	}
1866	else
1867	{
1868	//decompressDataSeries_double_3D_pwr(data, r1*r2, r3, r4, tdps);
1869	decompressDataSeries_double_3D_pwr_pre_log(data, r1*r2, r3, r4, tdps);
1870	//ToDO
1871	//decompressDataSeries_double_4D_pwr(data, r1, r2, r3, r4, tdps);
1872	}
1873	return;
1874	} else {
1875	data = (double)malloc(sizeof(double)*dataSeriesLength);
1876	int* rtypes;
1877	int validLength = computeBitNumRequired(dataSeriesLength);
1878	decompressBitArraybySimpleLZ77(&rtypes, tdps->rtypeArray, tdps->rtypeArray_size, dataSeriesLength, validLength);
1879	size_t count = 0;
1880	for (i = 0; i < dataSeriesLength; i++) {
1881	if (rtypes[i] == 1)
1882	(*data)[i] = tdps->reservedValue;
1883	else
1884	count++;
1885	}
1886	// get the decompressed data
1887	double* decmpData;
1888	if(errBoundMode < PW_REL)
1889	decompressDataSeries_double_4D(&decmpData, r1, r2, r3, r4, tdps);
1890	else
1891	//decompressDataSeries_double_3D_pwr(&decmpData, r1*r2, r3, r4, tdps);
1892	decompressDataSeries_double_3D_pwr_pre_log(&decmpData, r1*r2, r3, r4, tdps);
1893	//ToDo
1894	//decompressDataSeries_double_4D_pwr(&decmpData, r1, r2, r3, r4, tdps);
1895	// insert the decompressed data
1896	size_t k = 0;
1897	for (i = 0; i < dataSeriesLength; i++) {
1898	if (rtypes[i] == 0) {
1899	(*data)[i] = decmpData[k++];
1900	}
1901	}
1902	free(decmpData);
1903	free(rtypes);
1904	}
1905	}
1906	}
1907
1908	void decompressDataSeries_double_2D_nonblocked_with_blocked_regression(double** data, size_t r1, size_t r2, unsigned char* comp_data){
1909
1910	size_t dim0_offset = r2;
1911	size_t num_elements = r1 * r2;
1912
1913	data = (double)malloc(sizeof(double)*num_elements);
1914
1915	unsigned char * comp_data_pos = comp_data;
1916
1917	size_t block_size = bytesToInt_bigEndian(comp_data_pos);
1918	comp_data_pos += sizeof(int);
1919	// calculate block dims
1920	size_t num_x, num_y;
1921	SZ_COMPUTE_3D_NUMBER_OF_BLOCKS(r1, num_x, block_size);
1922	SZ_COMPUTE_3D_NUMBER_OF_BLOCKS(r2, num_y, block_size);
1923
1924	size_t split_index_x, split_index_y;
1925	size_t early_blockcount_x, early_blockcount_y;
1926	size_t late_blockcount_x, late_blockcount_y;
1927	SZ_COMPUTE_BLOCKCOUNT(r1, num_x, split_index_x, early_blockcount_x, late_blockcount_x);
1928	SZ_COMPUTE_BLOCKCOUNT(r2, num_y, split_index_y, early_blockcount_y, late_blockcount_y);
1929
1930	size_t num_blocks = num_x * num_y;
1931
1932	double realPrecision = bytesToDouble(comp_data_pos);
1933	comp_data_pos += sizeof(double);
1934	unsigned int intervals = bytesToInt_bigEndian(comp_data_pos);
1935	comp_data_pos += sizeof(int);
1936
1937	updateQuantizationInfo(intervals);
1938
1939	unsigned int tree_size = bytesToInt_bigEndian(comp_data_pos);
1940	comp_data_pos += sizeof(int);
1941
1942	int stateNum = 2*intervals;
1943	HuffmanTree* huffmanTree = createHuffmanTree(stateNum);
1944
1945	int nodeCount = bytesToInt_bigEndian(comp_data_pos);
1946
1947	node root = reconstruct_HuffTree_from_bytes_anyStates(huffmanTree,comp_data_pos+sizeof(int), nodeCount);
1948	comp_data_pos += sizeof(int) + tree_size;
1949
1950	double mean;
1951	unsigned char use_mean;
1952	memcpy(&use_mean, comp_data_pos, sizeof(unsigned char));
1953	comp_data_pos += sizeof(unsigned char);
1954	memcpy(&mean, comp_data_pos, sizeof(double));
1955	comp_data_pos += sizeof(double);
1956	size_t reg_count = 0;
1957
1958	unsigned char * indicator;
1959	size_t indicator_bitlength = (num_blocks - 1)/8 + 1;
1960	convertByteArray2IntArray_fast_1b(num_blocks, comp_data_pos, indicator_bitlength, &indicator);
1961	comp_data_pos += indicator_bitlength;
1962	for(size_t i=0; i<num_blocks; i++){
1963	if(!indicator[i]) reg_count ++;
1964	}
1965	//printf("reg_count: %ld\n", reg_count);
1966
1967	int coeff_intvRadius[3];
1968	int * coeff_result_type = (int ) malloc(num_blocks3*sizeof(int));
1969	int * coeff_type[3];
1970	double precision[3];
1971	double * coeff_unpred_data[3];
1972	if(reg_count > 0){
1973	for(int i=0; i<3; i++){
1974	precision[i] = bytesToDouble(comp_data_pos);
1975	comp_data_pos += sizeof(double);
1976	coeff_intvRadius[i] = bytesToInt_bigEndian(comp_data_pos);
1977	comp_data_pos += sizeof(int);
1978	unsigned int tree_size = bytesToInt_bigEndian(comp_data_pos);
1979	comp_data_pos += sizeof(int);
1980	int stateNum = 2coeff_intvRadius[i]2;
1981	HuffmanTree* huffmanTree = createHuffmanTree(stateNum);
1982	int nodeCount = bytesToInt_bigEndian(comp_data_pos);
1983	node root = reconstruct_HuffTree_from_bytes_anyStates(huffmanTree, comp_data_pos+sizeof(int), nodeCount);
1984	comp_data_pos += sizeof(int) + tree_size;
1985
1986	coeff_type[i] = coeff_result_type + i * num_blocks;
1987	size_t typeArray_size = bytesToSize(comp_data_pos);
1988	decode(comp_data_pos + sizeof(size_t), reg_count, root, coeff_type[i]);
1989	comp_data_pos += sizeof(size_t) + typeArray_size;
1990	int coeff_unpred_count = bytesToInt_bigEndian(comp_data_pos);
1991	comp_data_pos += sizeof(int);
1992	coeff_unpred_data[i] = (double *) comp_data_pos;
1993	comp_data_pos += coeff_unpred_count * sizeof(double);
1994	SZ_ReleaseHuffman(huffmanTree);
1995	}
1996	}
1997	double last_coefficients[3] = {0.0};
1998	int coeff_unpred_data_count[3] = {0};
1999	int coeff_index = 0;
2000	updateQuantizationInfo(intervals);
2001
2002	size_t total_unpred;
2003	memcpy(&total_unpred, comp_data_pos, sizeof(size_t));
2004	comp_data_pos += sizeof(size_t);
2005	double * unpred_data = (double *) comp_data_pos;
2006	comp_data_pos += total_unpred * sizeof(double);
2007
2008	int * result_type = (int ) malloc(num_elements sizeof(int));
2009	decode(comp_data_pos, num_elements, root, result_type);
2010	SZ_ReleaseHuffman(huffmanTree);
2011
2012	int intvRadius = exe_params->intvRadius;
2013
2014	int * type;
2015
2016	double * data_pos = *data;
2017	size_t offset_x, offset_y;
2018	size_t current_blockcount_x, current_blockcount_y;
2019	size_t cur_unpred_count;
2020
2021	unsigned char * indicator_pos = indicator;
2022	if(use_mean){
2023	type = result_type;
2024	for(size_t i=0; i<num_x; i++){
2025	for(size_t j=0; j<num_y; j++){
2026	offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
2027	offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
2028	data_pos = data + offset_x dim0_offset + offset_y;
2029
2030	current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
2031	current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
2032
2033	size_t current_block_elements = current_blockcount_x * current_blockcount_y;
2034	if(*indicator_pos){
2035	// decompress by SZ
2036
2037	double * block_data_pos = data_pos;
2038	double pred;
2039	size_t index = 0;
2040	int type_;
2041	// d11 is current data
2042	size_t unpredictable_count = 0;
2043	double d00, d01, d10;
2044	for(size_t ii=0; ii<current_blockcount_x; ii++){
2045	for(size_t jj=0; jj<current_blockcount_y; jj++){
2046	type_ = type[index];
2047	if(type_ == intvRadius){
2048	*block_data_pos = mean;
2049	}
2050	else if(type_ == 0){
2051	*block_data_pos = unpred_data[unpredictable_count ++];
2052	}
2053	else{
2054	d00 = d01 = d10 = 1;
2055	if(i == 0 && ii == 0){
2056	d00 = d01 = 0;
2057	}
2058	if(j == 0 && jj == 0){
2059	d00 = d10 = 0;
2060	}
2061	if(d00){
2062	d00 = block_data_pos[- dim0_offset - 1];
2063	}
2064	if(d01){
2065	d01 = block_data_pos[- dim0_offset];
2066	}
2067	if(d10){
2068	d10 = block_data_pos[- 1];
2069	}
2070	if(type_ < intvRadius) type_ += 1;
2071	pred = d10 + d01 - d00;
2072	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2073	}
2074	index ++;
2075	block_data_pos ++;
2076	}
2077	block_data_pos += dim0_offset - current_blockcount_y;
2078	}
2079	cur_unpred_count = unpredictable_count;
2080	}
2081	else{
2082	// decompress by regression
2083	{
2084	//restore regression coefficients
2085	double pred;
2086	int type_;
2087	for(int e=0; e<3; e++){
2088	type_ = coeff_type[e][coeff_index];
2089	if (type_ != 0){
2090	pred = last_coefficients[e];
2091	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
2092	}
2093	else{
2094	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
2095	coeff_unpred_data_count[e] ++;
2096	}
2097	}
2098	coeff_index ++;
2099	}
2100	{
2101	double * block_data_pos = data_pos;
2102	double pred;
2103	int type_;
2104	size_t index = 0;
2105	size_t unpredictable_count = 0;
2106	for(size_t ii=0; ii<current_blockcount_x; ii++){
2107	for(size_t jj=0; jj<current_blockcount_y; jj++){
2108	type_ = type[index];
2109	if (type_ != 0){
2110	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2];
2111	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2112	}
2113	else{
2114	*block_data_pos = unpred_data[unpredictable_count ++];
2115	}
2116
2117	index ++;
2118	block_data_pos ++;
2119	}
2120	block_data_pos += dim0_offset - current_blockcount_y;
2121	}
2122	cur_unpred_count = unpredictable_count;
2123	}
2124	}
2125
2126	type += current_block_elements;
2127	indicator_pos ++;
2128	unpred_data += cur_unpred_count;
2129	}
2130	}
2131	}
2132	else{
2133	type = result_type;
2134	for(size_t i=0; i<num_x; i++){
2135	for(size_t j=0; j<num_y; j++){
2136	offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
2137	offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
2138	data_pos = data + offset_x dim0_offset + offset_y;
2139
2140	current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
2141	current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
2142
2143	size_t current_block_elements = current_blockcount_x * current_blockcount_y;
2144	if(*indicator_pos){
2145	// decompress by SZ
2146
2147	double * block_data_pos = data_pos;
2148	double pred;
2149	size_t index = 0;
2150	int type_;
2151	// d11 is current data
2152	size_t unpredictable_count = 0;
2153	double d00, d01, d10;
2154	for(size_t ii=0; ii<current_blockcount_x; ii++){
2155	for(size_t jj=0; jj<current_blockcount_y; jj++){
2156	type_ = type[index];
2157	if(type_ == 0){
2158	*block_data_pos = unpred_data[unpredictable_count ++];
2159	}
2160	else{
2161	d00 = d01 = d10 = 1;
2162	if(i == 0 && ii == 0){
2163	d00 = d01 = 0;
2164	}
2165	if(j == 0 && jj == 0){
2166	d00 = d10 = 0;
2167	}
2168	if(d00){
2169	d00 = block_data_pos[- dim0_offset - 1];
2170	}
2171	if(d01){
2172	d01 = block_data_pos[- dim0_offset];
2173	}
2174	if(d10){
2175	d10 = block_data_pos[- 1];
2176	}
2177	pred = d10 + d01 - d00;
2178	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2179	}
2180	index ++;
2181	block_data_pos ++;
2182	}
2183	block_data_pos += dim0_offset - current_blockcount_y;
2184	}
2185	cur_unpred_count = unpredictable_count;
2186	}
2187	else{
2188	// decompress by regression
2189	{
2190	//restore regression coefficients
2191	double pred;
2192	int type_;
2193	for(int e=0; e<3; e++){
2194	type_ = coeff_type[e][coeff_index];
2195	if (type_ != 0){
2196	pred = last_coefficients[e];
2197	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
2198	}
2199	else{
2200	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
2201	coeff_unpred_data_count[e] ++;
2202	}
2203	}
2204	coeff_index ++;
2205	}
2206	{
2207	double * block_data_pos = data_pos;
2208	double pred;
2209	int type_;
2210	size_t index = 0;
2211	size_t unpredictable_count = 0;
2212	for(size_t ii=0; ii<current_blockcount_x; ii++){
2213	for(size_t jj=0; jj<current_blockcount_y; jj++){
2214	type_ = type[index];
2215	if (type_ != 0){
2216	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2];
2217	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2218	}
2219	else{
2220	*block_data_pos = unpred_data[unpredictable_count ++];
2221	}
2222	index ++;
2223	block_data_pos ++;
2224	}
2225	block_data_pos += dim0_offset - current_blockcount_y;
2226	}
2227	cur_unpred_count = unpredictable_count;
2228	}
2229	}
2230
2231	type += current_block_elements;
2232	indicator_pos ++;
2233	unpred_data += cur_unpred_count;
2234	}
2235	}
2236	}
2237	free(coeff_result_type);
2238
2239	free(indicator);
2240	free(result_type);
2241	}
2242
2243
2244	void decompressDataSeries_double_3D_nonblocked_with_blocked_regression(double** data, size_t r1, size_t r2, size_t r3, unsigned char* comp_data){
2245
2246	size_t dim0_offset = r2 * r3;
2247	size_t dim1_offset = r3;
2248	size_t num_elements = r1 * r2 * r3;
2249
2250	data = (double)malloc(sizeof(double)*num_elements);
2251
2252	unsigned char * comp_data_pos = comp_data;
2253
2254	size_t block_size = bytesToInt_bigEndian(comp_data_pos);
2255	comp_data_pos += sizeof(int);
2256	// calculate block dims
2257	size_t num_x, num_y, num_z;
2258	SZ_COMPUTE_3D_NUMBER_OF_BLOCKS(r1, num_x, block_size);
2259	SZ_COMPUTE_3D_NUMBER_OF_BLOCKS(r2, num_y, block_size);
2260	SZ_COMPUTE_3D_NUMBER_OF_BLOCKS(r3, num_z, block_size);
2261
2262	size_t split_index_x, split_index_y, split_index_z;
2263	size_t early_blockcount_x, early_blockcount_y, early_blockcount_z;
2264	size_t late_blockcount_x, late_blockcount_y, late_blockcount_z;
2265	SZ_COMPUTE_BLOCKCOUNT(r1, num_x, split_index_x, early_blockcount_x, late_blockcount_x);
2266	SZ_COMPUTE_BLOCKCOUNT(r2, num_y, split_index_y, early_blockcount_y, late_blockcount_y);
2267	SZ_COMPUTE_BLOCKCOUNT(r3, num_z, split_index_z, early_blockcount_z, late_blockcount_z);
2268
2269	size_t num_blocks = num_x * num_y * num_z;
2270
2271	double realPrecision = bytesToDouble(comp_data_pos);
2272	comp_data_pos += sizeof(double);
2273	unsigned int intervals = bytesToInt_bigEndian(comp_data_pos);
2274	comp_data_pos += sizeof(int);
2275
2276	updateQuantizationInfo(intervals);
2277
2278	unsigned int tree_size = bytesToInt_bigEndian(comp_data_pos);
2279	comp_data_pos += sizeof(int);
2280
2281	int stateNum = 2*intervals;
2282	HuffmanTree* huffmanTree = createHuffmanTree(stateNum);
2283
2284	int nodeCount = bytesToInt_bigEndian(comp_data_pos);
2285	node root = reconstruct_HuffTree_from_bytes_anyStates(huffmanTree,comp_data_pos+4, nodeCount);
2286	comp_data_pos += sizeof(int) + tree_size;
2287
2288	double mean;
2289	unsigned char use_mean;
2290	memcpy(&use_mean, comp_data_pos, sizeof(unsigned char));
2291	comp_data_pos += sizeof(unsigned char);
2292	memcpy(&mean, comp_data_pos, sizeof(double));
2293	comp_data_pos += sizeof(double);
2294	size_t reg_count = 0;
2295
2296	unsigned char * indicator;
2297	size_t indicator_bitlength = (num_blocks - 1)/8 + 1;
2298	convertByteArray2IntArray_fast_1b(num_blocks, comp_data_pos, indicator_bitlength, &indicator);
2299	comp_data_pos += indicator_bitlength;
2300	for(size_t i=0; i<num_blocks; i++){
2301	if(!indicator[i]) reg_count ++;
2302	}
2303
2304	int coeff_intvRadius[4];
2305	int * coeff_result_type = (int ) malloc(num_blocks4*sizeof(int));
2306	int * coeff_type[4];
2307	double precision[4];
2308	double * coeff_unpred_data[4];
2309	if(reg_count > 0){
2310	for(int i=0; i<4; i++){
2311	precision[i] = bytesToDouble(comp_data_pos);
2312	comp_data_pos += sizeof(double);
2313	coeff_intvRadius[i] = bytesToInt_bigEndian(comp_data_pos);
2314	comp_data_pos += sizeof(int);
2315	unsigned int tree_size = bytesToInt_bigEndian(comp_data_pos);
2316	comp_data_pos += sizeof(int);
2317	int stateNum = 2coeff_intvRadius[i]2;
2318	HuffmanTree* huffmanTree = createHuffmanTree(stateNum);
2319	int nodeCount = bytesToInt_bigEndian(comp_data_pos);
2320	node root = reconstruct_HuffTree_from_bytes_anyStates(huffmanTree, comp_data_pos+4, nodeCount);
2321	comp_data_pos += sizeof(int) + tree_size;
2322
2323	coeff_type[i] = coeff_result_type + i * num_blocks;
2324	size_t typeArray_size = bytesToSize(comp_data_pos);
2325	decode(comp_data_pos + sizeof(size_t), reg_count, root, coeff_type[i]);
2326	comp_data_pos += sizeof(size_t) + typeArray_size;
2327	int coeff_unpred_count = bytesToInt_bigEndian(comp_data_pos);
2328	comp_data_pos += sizeof(int);
2329	coeff_unpred_data[i] = (double *) comp_data_pos;
2330	comp_data_pos += coeff_unpred_count * sizeof(double);
2331	SZ_ReleaseHuffman(huffmanTree);
2332	}
2333	}
2334	double last_coefficients[4] = {0.0};
2335	int coeff_unpred_data_count[4] = {0};
2336	int coeff_index = 0;
2337	updateQuantizationInfo(intervals);
2338
2339	size_t total_unpred;
2340	memcpy(&total_unpred, comp_data_pos, sizeof(size_t));
2341	comp_data_pos += sizeof(size_t);
2342	double * unpred_data = (double *) comp_data_pos;
2343	comp_data_pos += total_unpred * sizeof(double);
2344
2345	int * result_type = (int ) malloc(num_elements sizeof(int));
2346	decode(comp_data_pos, num_elements, root, result_type);
2347	SZ_ReleaseHuffman(huffmanTree);
2348
2349	int intvRadius = exe_params->intvRadius;
2350
2351	int * type;
2352	double * data_pos = *data;
2353	size_t offset_x, offset_y, offset_z;
2354	size_t current_blockcount_x, current_blockcount_y, current_blockcount_z;
2355	size_t cur_unpred_count;
2356	unsigned char * indicator_pos = indicator;
2357	if(use_mean){
2358	// type = result_type;
2359
2360	// for(size_t i=0; i<num_x; i++){
2361	// for(size_t j=0; j<num_y; j++){
2362	// for(size_t k=0; k<num_z; k++){
2363	// offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
2364	// offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
2365	// offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
2366	// data_pos = data + offset_x dim0_offset + offset_y * dim1_offset + offset_z;
2367
2368	// current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
2369	// current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
2370	// current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
2371
2372	// // type_offset = offset_x * dim0_offset + offset_y * current_blockcount_x * dim1_offset + offset_z * current_blockcount_x * current_blockcount_y;
2373	// // type = result_type + type_offset;
2374	// size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
2375	// // index = i * num_y * num_z + j * num_z + k;
2376
2377	// // printf("i j k: %ld %ld %ld\toffset: %ld %ld %ld\tindicator: %ld\n", i, j, k, offset_x, offset_y, offset_z, indicator[index]);
2378	// if(*indicator_pos){
2379	// // decompress by SZ
2380	// // cur_unpred_count = decompressDataSeries_double_3D_blocked_nonblock_pred(data_pos, r1, r2, r3, current_blockcount_x, current_blockcount_y, current_blockcount_z, i, j, k, realPrecision, type, unpred_data);
2381	// double * block_data_pos = data_pos;
2382	// double pred;
2383	// size_t index = 0;
2384	// int type_;
2385	// // d111 is current data
2386	// size_t unpredictable_count = 0;
2387	// double d000, d001, d010, d011, d100, d101, d110;
2388	// for(size_t ii=0; ii<current_blockcount_x; ii++){
2389	// for(size_t jj=0; jj<current_blockcount_y; jj++){
2390	// for(size_t kk=0; kk<current_blockcount_z; kk++){
2391	// type_ = type[index];
2392	// if(type_ == intvRadius){
2393	// *block_data_pos = mean;
2394	// }
2395	// else if(type_ == 0){
2396	// *block_data_pos = unpred_data[unpredictable_count ++];
2397	// }
2398	// else{
2399	// d000 = d001 = d010 = d011 = d100 = d101 = d110 = 1;
2400	// if(i == 0 && ii == 0){
2401	// d000 = d001 = d010 = d011 = 0;
2402	// }
2403	// if(j == 0 && jj == 0){
2404	// d000 = d001 = d100 = d101 = 0;
2405	// }
2406	// if(k == 0 && kk == 0){
2407	// d000 = d010 = d100 = d110 = 0;
2408	// }
2409	// if(d000){
2410	// d000 = block_data_pos[- dim0_offset - dim1_offset - 1];
2411	// }
2412	// if(d001){
2413	// d001 = block_data_pos[- dim0_offset - dim1_offset];
2414	// }
2415	// if(d010){
2416	// d010 = block_data_pos[- dim0_offset - 1];
2417	// }
2418	// if(d011){
2419	// d011 = block_data_pos[- dim0_offset];
2420	// }
2421	// if(d100){
2422	// d100 = block_data_pos[- dim1_offset - 1];
2423	// }
2424	// if(d101){
2425	// d101 = block_data_pos[- dim1_offset];
2426	// }
2427	// if(d110){
2428	// d110 = block_data_pos[- 1];
2429	// }
2430	// if(type_ < intvRadius) type_ += 1;
2431	// pred = d110 + d101 + d011 - d100 - d010 - d001 + d000;
2432	// block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2433	// }
2434	// index ++;
2435	// block_data_pos ++;
2436	// }
2437	// block_data_pos += dim1_offset - current_blockcount_z;
2438	// }
2439	// block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2440	// }
2441	// cur_unpred_count = unpredictable_count;
2442	// }
2443	// else{
2444	// // decompress by regression
2445	// {
2446	// //restore regression coefficients
2447	// double pred;
2448	// int type_;
2449	// for(int e=0; e<4; e++){
2450	// // if(i == 0 && j == 0 && k == 19){
2451	// // printf("~\n");
2452	// // }
2453	// type_ = coeff_type[e][coeff_index];
2454	// if (type_ != 0){
2455	// pred = last_coefficients[e];
2456	// last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
2457	// }
2458	// else{
2459	// last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
2460	// coeff_unpred_data_count[e] ++;
2461	// }
2462	// if(fabs(last_coefficients[e]) > 10000){
2463	// printf("%d %d %d-%d: pred %.4f type %d precision %.4g last_coefficients %.4g\n", i, j, k, e, pred, type_, precision[e], last_coefficients[e]);
2464	// exit(0);
2465	// }
2466	// }
2467	// coeff_index ++;
2468	// }
2469	// {
2470	// double * block_data_pos = data_pos;
2471	// double pred;
2472	// int type_;
2473	// size_t index = 0;
2474	// size_t unpredictable_count = 0;
2475	// for(size_t ii=0; ii<current_blockcount_x; ii++){
2476	// for(size_t jj=0; jj<current_blockcount_y; jj++){
2477	// for(size_t kk=0; kk<current_blockcount_z; kk++){
2478	// if(block_data_pos - (*data) == 19470788){
2479	// printf("dec stop\n");
2480	// }
2481
2482	// type_ = type[index];
2483	// if (type_ != 0){
2484	// pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
2485	// block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2486	// }
2487	// else{
2488	// *block_data_pos = unpred_data[unpredictable_count ++];
2489	// }
2490	// index ++;
2491	// block_data_pos ++;
2492	// }
2493	// block_data_pos += dim1_offset - current_blockcount_z;
2494	// }
2495	// block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2496	// }
2497	// cur_unpred_count = unpredictable_count;
2498	// }
2499	// }
2500
2501	// type += current_block_elements;
2502	// indicator_pos ++;
2503	// unpred_data += cur_unpred_count;
2504	// // decomp_unpred += cur_unpred_count;
2505	// // printf("block comp done, data_offset from %ld to %ld: diff %ld\n", data, data_pos, data_pos - data);
2506	// // fflush(stdout);
2507	// }
2508	// }
2509	// }
2510
2511	type = result_type;
2512	// i == 0
2513	{
2514	// j == 0
2515	{
2516	// k == 0
2517	{
2518	data_pos = *data;
2519
2520	current_blockcount_x = early_blockcount_x;
2521	current_blockcount_y = early_blockcount_y;
2522	current_blockcount_z = early_blockcount_z;
2523	size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
2524	if(*indicator_pos){
2525	// decompress by SZ
2526	double * block_data_pos = data_pos;
2527	double pred;
2528	size_t index = 0;
2529	int type_;
2530	size_t unpredictable_count = 0;
2531	// ii == 0
2532	{
2533	// jj == 0
2534	{
2535	{
2536	// kk == 0
2537	type_ = type[index];
2538	if(type_ == intvRadius){
2539	*block_data_pos = mean;
2540	}
2541	else if(type_ == 0){
2542	*block_data_pos = unpred_data[unpredictable_count ++];
2543	}
2544	else{
2545	if(type_ < intvRadius) type_ += 1;
2546	pred = 0;
2547	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2548	}
2549	index ++;
2550	block_data_pos ++;
2551	}
2552	for(size_t kk=1; kk<current_blockcount_z; kk++){
2553	type_ = type[index];
2554	if(type_ == intvRadius){
2555	*block_data_pos = mean;
2556	}
2557	else if(type_ == 0){
2558	*block_data_pos = unpred_data[unpredictable_count ++];
2559	}
2560	else{
2561	if(type_ < intvRadius) type_ += 1;
2562	pred = block_data_pos[- 1];
2563	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2564	}
2565	index ++;
2566	block_data_pos ++;
2567	}
2568	block_data_pos += dim1_offset - current_blockcount_z;
2569	}
2570	for(size_t jj=1; jj<current_blockcount_y; jj++){
2571	{
2572	// kk == 0
2573	type_ = type[index];
2574	if(type_ == intvRadius){
2575	*block_data_pos = mean;
2576	}
2577	else if(type_ == 0){
2578	*block_data_pos = unpred_data[unpredictable_count ++];
2579	}
2580	else{
2581	if(type_ < intvRadius) type_ += 1;
2582	pred = block_data_pos[- dim1_offset];
2583	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2584	}
2585	index ++;
2586	block_data_pos ++;
2587	}
2588	for(size_t kk=1; kk<current_blockcount_z; kk++){
2589	type_ = type[index];
2590	if(type_ == intvRadius){
2591	*block_data_pos = mean;
2592	}
2593	else if(type_ == 0){
2594	*block_data_pos = unpred_data[unpredictable_count ++];
2595	}
2596	else{
2597	if(type_ < intvRadius) type_ += 1;
2598	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] - block_data_pos[- dim1_offset - 1];
2599	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2600	}
2601	index ++;
2602	block_data_pos ++;
2603	}
2604	block_data_pos += dim1_offset - current_blockcount_z;
2605	}
2606	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2607	}
2608	for(size_t ii=1; ii<current_blockcount_x; ii++){
2609	// jj == 0
2610	{
2611	{
2612	// kk == 0
2613	type_ = type[index];
2614	if(type_ == intvRadius){
2615	*block_data_pos = mean;
2616	}
2617	else if(type_ == 0){
2618	*block_data_pos = unpred_data[unpredictable_count ++];
2619	}
2620	else{
2621	if(type_ < intvRadius) type_ += 1;
2622	pred = block_data_pos[- dim0_offset];
2623	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2624	}
2625	index ++;
2626	block_data_pos ++;
2627	}
2628	for(size_t kk=1; kk<current_blockcount_z; kk++){
2629	type_ = type[index];
2630	if(type_ == intvRadius){
2631	*block_data_pos = mean;
2632	}
2633	else if(type_ == 0){
2634	*block_data_pos = unpred_data[unpredictable_count ++];
2635	}
2636	else{
2637	if(type_ < intvRadius) type_ += 1;
2638	pred = block_data_pos[- 1] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - 1];
2639	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2640	}
2641	index ++;
2642	block_data_pos ++;
2643	}
2644	block_data_pos += dim1_offset - current_blockcount_z;
2645	}
2646	for(size_t jj=1; jj<current_blockcount_y; jj++){
2647	{
2648	// kk == 0
2649	type_ = type[index];
2650	if(type_ == intvRadius){
2651	*block_data_pos = mean;
2652	}
2653	else if(type_ == 0){
2654	*block_data_pos = unpred_data[unpredictable_count ++];
2655	}
2656	else{
2657	if(type_ < intvRadius) type_ += 1;
2658	pred = block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - dim1_offset];
2659	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2660	}
2661	index ++;
2662	block_data_pos ++;
2663	}
2664	for(size_t kk=1; kk<current_blockcount_z; kk++){
2665	type_ = type[index];
2666	if(type_ == intvRadius){
2667	*block_data_pos = mean;
2668	}
2669	else if(type_ == 0){
2670	*block_data_pos = unpred_data[unpredictable_count ++];
2671	}
2672	else{
2673	if(type_ < intvRadius) type_ += 1;
2674	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
2675	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2676	}
2677	index ++;
2678	block_data_pos ++;
2679	}
2680	block_data_pos += dim1_offset - current_blockcount_z;
2681	}
2682	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2683	}
2684	cur_unpred_count = unpredictable_count;
2685	}
2686	else{
2687	// decompress by regression
2688	{
2689	//restore regression coefficients
2690	double pred;
2691	int type_;
2692	for(int e=0; e<4; e++){
2693	type_ = coeff_type[e][coeff_index];
2694	if (type_ != 0){
2695	pred = last_coefficients[e];
2696	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
2697	}
2698	else{
2699	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
2700	coeff_unpred_data_count[e] ++;
2701	}
2702	}
2703	coeff_index ++;
2704	}
2705	{
2706	double * block_data_pos = data_pos;
2707	double pred;
2708	int type_;
2709	size_t index = 0;
2710	size_t unpredictable_count = 0;
2711	for(size_t ii=0; ii<current_blockcount_x; ii++){
2712	for(size_t jj=0; jj<current_blockcount_y; jj++){
2713	for(size_t kk=0; kk<current_blockcount_z; kk++){
2714	type_ = type[index];
2715	if (type_ != 0){
2716	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
2717	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2718	}
2719	else{
2720	*block_data_pos = unpred_data[unpredictable_count ++];
2721	}
2722	index ++;
2723	block_data_pos ++;
2724	}
2725	block_data_pos += dim1_offset - current_blockcount_z;
2726	}
2727	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2728	}
2729	cur_unpred_count = unpredictable_count;
2730	}
2731	}
2732	indicator_pos ++;
2733	type += current_block_elements;
2734	unpred_data += cur_unpred_count;
2735	} // end k == 0
2736	// i == 0 j == 0 k != 0
2737	for(size_t k=1; k<num_z; k++){
2738	offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
2739	data_pos = *data + offset_z;
2740
2741	current_blockcount_x = early_blockcount_x;
2742	current_blockcount_y = early_blockcount_y;
2743	current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
2744
2745	size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
2746	if(*indicator_pos){
2747	// decompress by SZ
2748	double * block_data_pos = data_pos;
2749	double pred;
2750	size_t index = 0;
2751	int type_;
2752	size_t unpredictable_count = 0;
2753	// ii == 0
2754	{
2755	// jj == 0
2756	{
2757	for(size_t kk=0; kk<current_blockcount_z; kk++){
2758	type_ = type[index];
2759	if(type_ == intvRadius){
2760	*block_data_pos = mean;
2761	}
2762	else if(type_ == 0){
2763	*block_data_pos = unpred_data[unpredictable_count ++];
2764	}
2765	else{
2766	if(type_ < intvRadius) type_ += 1;
2767	pred = block_data_pos[- 1];
2768	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2769	}
2770	index ++;
2771	block_data_pos ++;
2772	}
2773	block_data_pos += dim1_offset - current_blockcount_z;
2774	}
2775	for(size_t jj=1; jj<current_blockcount_y; jj++){
2776	for(size_t kk=0; kk<current_blockcount_z; kk++){
2777	type_ = type[index];
2778	if(type_ == intvRadius){
2779	*block_data_pos = mean;
2780	}
2781	else if(type_ == 0){
2782	*block_data_pos = unpred_data[unpredictable_count ++];
2783	}
2784	else{
2785	if(type_ < intvRadius) type_ += 1;
2786	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] - block_data_pos[- dim1_offset - 1];
2787	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2788	}
2789	index ++;
2790	block_data_pos ++;
2791	}
2792	block_data_pos += dim1_offset - current_blockcount_z;
2793	}
2794	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2795	}
2796	for(size_t ii=1; ii<current_blockcount_x; ii++){
2797	// jj == 0
2798	{
2799	for(size_t kk=0; kk<current_blockcount_z; kk++){
2800	type_ = type[index];
2801	if(type_ == intvRadius){
2802	*block_data_pos = mean;
2803	}
2804	else if(type_ == 0){
2805	*block_data_pos = unpred_data[unpredictable_count ++];
2806	}
2807	else{
2808	if(type_ < intvRadius) type_ += 1;
2809	pred = block_data_pos[- 1] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - 1];
2810	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2811	}
2812	index ++;
2813	block_data_pos ++;
2814	}
2815	block_data_pos += dim1_offset - current_blockcount_z;
2816	}
2817	for(size_t jj=1; jj<current_blockcount_y; jj++){
2818	for(size_t kk=0; kk<current_blockcount_z; kk++){
2819	type_ = type[index];
2820	if(type_ == intvRadius){
2821	*block_data_pos = mean;
2822	}
2823	else if(type_ == 0){
2824	*block_data_pos = unpred_data[unpredictable_count ++];
2825	}
2826	else{
2827	if(type_ < intvRadius) type_ += 1;
2828	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
2829	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2830	}
2831	index ++;
2832	block_data_pos ++;
2833	}
2834	block_data_pos += dim1_offset - current_blockcount_z;
2835	}
2836	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2837	}
2838	cur_unpred_count = unpredictable_count;
2839	}
2840	else{
2841	// decompress by regression
2842	{
2843	//restore regression coefficients
2844	double pred;
2845	int type_;
2846	for(int e=0; e<4; e++){
2847	type_ = coeff_type[e][coeff_index];
2848	if (type_ != 0){
2849	pred = last_coefficients[e];
2850	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
2851	}
2852	else{
2853	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
2854	coeff_unpred_data_count[e] ++;
2855	}
2856	}
2857	coeff_index ++;
2858	}
2859	{
2860	double * block_data_pos = data_pos;
2861	double pred;
2862	int type_;
2863	size_t index = 0;
2864	size_t unpredictable_count = 0;
2865	for(size_t ii=0; ii<current_blockcount_x; ii++){
2866	for(size_t jj=0; jj<current_blockcount_y; jj++){
2867	for(size_t kk=0; kk<current_blockcount_z; kk++){
2868	type_ = type[index];
2869	if (type_ != 0){
2870	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
2871	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2872	}
2873	else{
2874	*block_data_pos = unpred_data[unpredictable_count ++];
2875	}
2876	index ++;
2877	block_data_pos ++;
2878	}
2879	block_data_pos += dim1_offset - current_blockcount_z;
2880	}
2881	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2882	}
2883	cur_unpred_count = unpredictable_count;
2884	}
2885	}
2886	indicator_pos ++;
2887	type += current_block_elements;
2888	unpred_data += cur_unpred_count;
2889	}
2890	}// end j==0
2891	for(size_t j=1; j<num_y; j++){
2892	// k == 0
2893	{
2894	offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
2895	data_pos = data + offset_y dim1_offset;
2896
2897	current_blockcount_x = early_blockcount_x;
2898	current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
2899	current_blockcount_z = early_blockcount_z;
2900	size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
2901	if(*indicator_pos){
2902	// decompress by SZ
2903	double * block_data_pos = data_pos;
2904	double pred;
2905	size_t index = 0;
2906	int type_;
2907	size_t unpredictable_count = 0;
2908	// ii == 0
2909	{
2910	for(size_t jj=0; jj<current_blockcount_y; jj++){
2911	{
2912	// kk == 0
2913	type_ = type[index];
2914	if(type_ == intvRadius){
2915	*block_data_pos = mean;
2916	}
2917	else if(type_ == 0){
2918	*block_data_pos = unpred_data[unpredictable_count ++];
2919	}
2920	else{
2921	if(type_ < intvRadius) type_ += 1;
2922	pred = block_data_pos[- dim1_offset];
2923	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2924	}
2925	index ++;
2926	block_data_pos ++;
2927	}
2928	for(size_t kk=1; kk<current_blockcount_z; kk++){
2929	type_ = type[index];
2930	if(type_ == intvRadius){
2931	*block_data_pos = mean;
2932	}
2933	else if(type_ == 0){
2934	*block_data_pos = unpred_data[unpredictable_count ++];
2935	}
2936	else{
2937	if(type_ < intvRadius) type_ += 1;
2938	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] - block_data_pos[- dim1_offset - 1];
2939	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2940	}
2941	index ++;
2942	block_data_pos ++;
2943	}
2944	block_data_pos += dim1_offset - current_blockcount_z;
2945	}
2946	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2947	}
2948	for(size_t ii=1; ii<current_blockcount_x; ii++){
2949	for(size_t jj=0; jj<current_blockcount_y; jj++){
2950	{
2951	// kk == 0
2952	type_ = type[index];
2953	if(type_ == intvRadius){
2954	*block_data_pos = mean;
2955	}
2956	else if(type_ == 0){
2957	*block_data_pos = unpred_data[unpredictable_count ++];
2958	}
2959	else{
2960	if(type_ < intvRadius) type_ += 1;
2961	pred = block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - dim1_offset];
2962	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2963	}
2964	index ++;
2965	block_data_pos ++;
2966	}
2967	for(size_t kk=1; kk<current_blockcount_z; kk++){
2968	type_ = type[index];
2969	if(type_ == intvRadius){
2970	*block_data_pos = mean;
2971	}
2972	else if(type_ == 0){
2973	*block_data_pos = unpred_data[unpredictable_count ++];
2974	}
2975	else{
2976	if(type_ < intvRadius) type_ += 1;
2977	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
2978	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
2979	}
2980	index ++;
2981	block_data_pos ++;
2982	}
2983	block_data_pos += dim1_offset - current_blockcount_z;
2984	}
2985	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2986	}
2987	cur_unpred_count = unpredictable_count;
2988	}
2989	else{
2990	// decompress by regression
2991	{
2992	//restore regression coefficients
2993	double pred;
2994	int type_;
2995	for(int e=0; e<4; e++){
2996	type_ = coeff_type[e][coeff_index];
2997	if (type_ != 0){
2998	pred = last_coefficients[e];
2999	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
3000	}
3001	else{
3002	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
3003	coeff_unpred_data_count[e] ++;
3004	}
3005	}
3006	coeff_index ++;
3007	}
3008	{
3009	double * block_data_pos = data_pos;
3010	double pred;
3011	int type_;
3012	size_t index = 0;
3013	size_t unpredictable_count = 0;
3014	for(size_t ii=0; ii<current_blockcount_x; ii++){
3015	for(size_t jj=0; jj<current_blockcount_y; jj++){
3016	for(size_t kk=0; kk<current_blockcount_z; kk++){
3017	type_ = type[index];
3018	if (type_ != 0){
3019	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
3020	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3021	}
3022	else{
3023	*block_data_pos = unpred_data[unpredictable_count ++];
3024	}
3025	index ++;
3026	block_data_pos ++;
3027	}
3028	block_data_pos += dim1_offset - current_blockcount_z;
3029	}
3030	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3031	}
3032	cur_unpred_count = unpredictable_count;
3033	}
3034	}
3035	indicator_pos ++;
3036	type += current_block_elements;
3037	unpred_data += cur_unpred_count;
3038	} // end k == 0
3039	for(size_t k=1; k<num_z; k++){
3040	offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
3041	offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
3042	data_pos = data + offset_y dim1_offset + offset_z;
3043
3044	current_blockcount_x = early_blockcount_x;
3045	current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
3046	current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
3047
3048	size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
3049	if(*indicator_pos){
3050	// decompress by SZ
3051	double * block_data_pos = data_pos;
3052	double pred;
3053	size_t index = 0;
3054	int type_;
3055	size_t unpredictable_count = 0;
3056	// ii == 0
3057	{
3058	for(size_t jj=0; jj<current_blockcount_y; jj++){
3059	for(size_t kk=0; kk<current_blockcount_z; kk++){
3060	type_ = type[index];
3061	if(type_ == intvRadius){
3062	*block_data_pos = mean;
3063	}
3064	else if(type_ == 0){
3065	*block_data_pos = unpred_data[unpredictable_count ++];
3066	}
3067	else{
3068	if(type_ < intvRadius) type_ += 1;
3069	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] - block_data_pos[- dim1_offset - 1];
3070	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3071	}
3072	index ++;
3073	block_data_pos ++;
3074	}
3075	block_data_pos += dim1_offset - current_blockcount_z;
3076	}
3077	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3078	}
3079	for(size_t ii=1; ii<current_blockcount_x; ii++){
3080	for(size_t jj=0; jj<current_blockcount_y; jj++){
3081	for(size_t kk=0; kk<current_blockcount_z; kk++){
3082	type_ = type[index];
3083	if(type_ == intvRadius){
3084	*block_data_pos = mean;
3085	}
3086	else if(type_ == 0){
3087	*block_data_pos = unpred_data[unpredictable_count ++];
3088	}
3089	else{
3090	if(type_ < intvRadius) type_ += 1;
3091	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
3092	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3093	}
3094	index ++;
3095	block_data_pos ++;
3096	}
3097	block_data_pos += dim1_offset - current_blockcount_z;
3098	}
3099	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3100	}
3101	cur_unpred_count = unpredictable_count;
3102	}
3103	else{
3104	// decompress by regression
3105	{
3106	//restore regression coefficients
3107	double pred;
3108	int type_;
3109	for(int e=0; e<4; e++){
3110	type_ = coeff_type[e][coeff_index];
3111	if (type_ != 0){
3112	pred = last_coefficients[e];
3113	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
3114	}
3115	else{
3116	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
3117	coeff_unpred_data_count[e] ++;
3118	}
3119	}
3120	coeff_index ++;
3121	}
3122	{
3123	double * block_data_pos = data_pos;
3124	double pred;
3125	int type_;
3126	size_t index = 0;
3127	size_t unpredictable_count = 0;
3128	for(size_t ii=0; ii<current_blockcount_x; ii++){
3129	for(size_t jj=0; jj<current_blockcount_y; jj++){
3130	for(size_t kk=0; kk<current_blockcount_z; kk++){
3131	type_ = type[index];
3132	if (type_ != 0){
3133	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
3134	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3135	}
3136	else{
3137	*block_data_pos = unpred_data[unpredictable_count ++];
3138	}
3139	index ++;
3140	block_data_pos ++;
3141	}
3142	block_data_pos += dim1_offset - current_blockcount_z;
3143	}
3144	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3145	}
3146	cur_unpred_count = unpredictable_count;
3147	}
3148	}
3149	indicator_pos ++;
3150	type += current_block_elements;
3151	unpred_data += cur_unpred_count;
3152	}
3153	}
3154	} // end i==0
3155	for(size_t i=1; i<num_x; i++){
3156	// j == 0
3157	{
3158	// k == 0
3159	{
3160	offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
3161	data_pos = data + offset_x dim0_offset;
3162
3163	current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
3164	current_blockcount_y = early_blockcount_y;
3165	current_blockcount_z = early_blockcount_z;
3166	size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
3167	if(*indicator_pos){
3168	// decompress by SZ
3169	double * block_data_pos = data_pos;
3170	double pred;
3171	size_t index = 0;
3172	int type_;
3173	size_t unpredictable_count = 0;
3174	for(size_t ii=0; ii<current_blockcount_x; ii++){
3175	// jj == 0
3176	{
3177	{
3178	// kk == 0
3179	type_ = type[index];
3180	if(type_ == intvRadius){
3181	*block_data_pos = mean;
3182	}
3183	else if(type_ == 0){
3184	*block_data_pos = unpred_data[unpredictable_count ++];
3185	}
3186	else{
3187	if(type_ < intvRadius) type_ += 1;
3188	pred = block_data_pos[- dim0_offset];
3189	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3190	}
3191	index ++;
3192	block_data_pos ++;
3193	}
3194	for(size_t kk=1; kk<current_blockcount_z; kk++){
3195	type_ = type[index];
3196	if(type_ == intvRadius){
3197	*block_data_pos = mean;
3198	}
3199	else if(type_ == 0){
3200	*block_data_pos = unpred_data[unpredictable_count ++];
3201	}
3202	else{
3203	if(type_ < intvRadius) type_ += 1;
3204	pred = block_data_pos[- 1] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - 1];
3205	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3206	}
3207	index ++;
3208	block_data_pos ++;
3209	}
3210	block_data_pos += dim1_offset - current_blockcount_z;
3211	}
3212	for(size_t jj=1; jj<current_blockcount_y; jj++){
3213	{
3214	// kk == 0
3215	type_ = type[index];
3216	if(type_ == intvRadius){
3217	*block_data_pos = mean;
3218	}
3219	else if(type_ == 0){
3220	*block_data_pos = unpred_data[unpredictable_count ++];
3221	}
3222	else{
3223	if(type_ < intvRadius) type_ += 1;
3224	pred = block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - dim1_offset];
3225	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3226	}
3227	index ++;
3228	block_data_pos ++;
3229	}
3230	for(size_t kk=1; kk<current_blockcount_z; kk++){
3231	type_ = type[index];
3232	if(type_ == intvRadius){
3233	*block_data_pos = mean;
3234	}
3235	else if(type_ == 0){
3236	*block_data_pos = unpred_data[unpredictable_count ++];
3237	}
3238	else{
3239	if(type_ < intvRadius) type_ += 1;
3240	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
3241	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3242	}
3243	index ++;
3244	block_data_pos ++;
3245	}
3246	block_data_pos += dim1_offset - current_blockcount_z;
3247	}
3248	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3249	}
3250	cur_unpred_count = unpredictable_count;
3251	}
3252	else{
3253	// decompress by regression
3254	{
3255	//restore regression coefficients
3256	double pred;
3257	int type_;
3258	for(int e=0; e<4; e++){
3259	type_ = coeff_type[e][coeff_index];
3260	if (type_ != 0){
3261	pred = last_coefficients[e];
3262	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
3263	}
3264	else{
3265	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
3266	coeff_unpred_data_count[e] ++;
3267	}
3268	}
3269	coeff_index ++;
3270	}
3271	{
3272	double * block_data_pos = data_pos;
3273	double pred;
3274	int type_;
3275	size_t index = 0;
3276	size_t unpredictable_count = 0;
3277	for(size_t ii=0; ii<current_blockcount_x; ii++){
3278	for(size_t jj=0; jj<current_blockcount_y; jj++){
3279	for(size_t kk=0; kk<current_blockcount_z; kk++){
3280	type_ = type[index];
3281	if (type_ != 0){
3282	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
3283	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3284	}
3285	else{
3286	*block_data_pos = unpred_data[unpredictable_count ++];
3287	}
3288	index ++;
3289	block_data_pos ++;
3290	}
3291	block_data_pos += dim1_offset - current_blockcount_z;
3292	}
3293	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3294	}
3295	cur_unpred_count = unpredictable_count;
3296	}
3297	}
3298	indicator_pos ++;
3299	type += current_block_elements;
3300	unpred_data += cur_unpred_count;
3301	} // end k == 0
3302	for(size_t k=1; k<num_z; k++){
3303	offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
3304	offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
3305	data_pos = data + offset_x dim0_offset + offset_z;
3306
3307	current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
3308	current_blockcount_y = early_blockcount_y;
3309	current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
3310	size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
3311	if(*indicator_pos){
3312	// decompress by SZ
3313	double * block_data_pos = data_pos;
3314	double pred;
3315	size_t index = 0;
3316	int type_;
3317	size_t unpredictable_count = 0;
3318	for(size_t ii=0; ii<current_blockcount_x; ii++){
3319	// jj == 0
3320	{
3321	for(size_t kk=0; kk<current_blockcount_z; kk++){
3322	type_ = type[index];
3323	if(type_ == intvRadius){
3324	*block_data_pos = mean;
3325	}
3326	else if(type_ == 0){
3327	*block_data_pos = unpred_data[unpredictable_count ++];
3328	}
3329	else{
3330	if(type_ < intvRadius) type_ += 1;
3331	pred = block_data_pos[- 1] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - 1];
3332	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3333	}
3334	index ++;
3335	block_data_pos ++;
3336	}
3337	block_data_pos += dim1_offset - current_blockcount_z;
3338	}
3339	for(size_t jj=1; jj<current_blockcount_y; jj++){
3340	for(size_t kk=0; kk<current_blockcount_z; kk++){
3341	type_ = type[index];
3342	if(type_ == intvRadius){
3343	*block_data_pos = mean;
3344	}
3345	else if(type_ == 0){
3346	*block_data_pos = unpred_data[unpredictable_count ++];
3347	}
3348	else{
3349	if(type_ < intvRadius) type_ += 1;
3350	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
3351	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3352	}
3353	index ++;
3354	block_data_pos ++;
3355	}
3356	block_data_pos += dim1_offset - current_blockcount_z;
3357	}
3358	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3359	}
3360	cur_unpred_count = unpredictable_count;
3361	}
3362	else{
3363	// decompress by regression
3364	{
3365	//restore regression coefficients
3366	double pred;
3367	int type_;
3368	for(int e=0; e<4; e++){
3369	type_ = coeff_type[e][coeff_index];
3370	if (type_ != 0){
3371	pred = last_coefficients[e];
3372	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
3373	}
3374	else{
3375	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
3376	coeff_unpred_data_count[e] ++;
3377	}
3378	}
3379	coeff_index ++;
3380	}
3381	{
3382	double * block_data_pos = data_pos;
3383	double pred;
3384	int type_;
3385	size_t index = 0;
3386	size_t unpredictable_count = 0;
3387	for(size_t ii=0; ii<current_blockcount_x; ii++){
3388	for(size_t jj=0; jj<current_blockcount_y; jj++){
3389	for(size_t kk=0; kk<current_blockcount_z; kk++){
3390	type_ = type[index];
3391	if (type_ != 0){
3392	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
3393	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3394	}
3395	else{
3396	*block_data_pos = unpred_data[unpredictable_count ++];
3397	}
3398	index ++;
3399	block_data_pos ++;
3400	}
3401	block_data_pos += dim1_offset - current_blockcount_z;
3402	}
3403	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3404	}
3405	cur_unpred_count = unpredictable_count;
3406	}
3407	}
3408	indicator_pos ++;
3409	type += current_block_elements;
3410	unpred_data += cur_unpred_count;
3411	}
3412	}// end j = 0
3413	for(size_t j=1; j<num_y; j++){
3414	// k == 0
3415	{
3416	offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
3417	offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
3418	data_pos = data + offset_x dim0_offset + offset_y * dim1_offset;
3419
3420	current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
3421	current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
3422	current_blockcount_z = early_blockcount_z;
3423	size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
3424	if(*indicator_pos){
3425	// decompress by SZ
3426	double * block_data_pos = data_pos;
3427	double pred;
3428	size_t index = 0;
3429	int type_;
3430	size_t unpredictable_count = 0;
3431	for(size_t ii=0; ii<current_blockcount_x; ii++){
3432	for(size_t jj=0; jj<current_blockcount_y; jj++){
3433	{
3434	// kk == 0
3435	type_ = type[index];
3436	if(type_ == intvRadius){
3437	*block_data_pos = mean;
3438	}
3439	else if(type_ == 0){
3440	*block_data_pos = unpred_data[unpredictable_count ++];
3441	}
3442	else{
3443	if(type_ < intvRadius) type_ += 1;
3444	pred = block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - dim1_offset];
3445	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3446	}
3447	index ++;
3448	block_data_pos ++;
3449	}
3450	for(size_t kk=1; kk<current_blockcount_z; kk++){
3451	type_ = type[index];
3452	if(type_ == intvRadius){
3453	*block_data_pos = mean;
3454	}
3455	else if(type_ == 0){
3456	*block_data_pos = unpred_data[unpredictable_count ++];
3457	}
3458	else{
3459	if(type_ < intvRadius) type_ += 1;
3460	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
3461	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3462	}
3463	index ++;
3464	block_data_pos ++;
3465	}
3466	block_data_pos += dim1_offset - current_blockcount_z;
3467	}
3468	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3469	}
3470	cur_unpred_count = unpredictable_count;
3471	}
3472	else{
3473	// decompress by regression
3474	{
3475	//restore regression coefficients
3476	double pred;
3477	int type_;
3478	for(int e=0; e<4; e++){
3479	type_ = coeff_type[e][coeff_index];
3480	if (type_ != 0){
3481	pred = last_coefficients[e];
3482	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
3483	}
3484	else{
3485	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
3486	coeff_unpred_data_count[e] ++;
3487	}
3488	}
3489	coeff_index ++;
3490	}
3491	{
3492	double * block_data_pos = data_pos;
3493	double pred;
3494	int type_;
3495	size_t index = 0;
3496	size_t unpredictable_count = 0;
3497	for(size_t ii=0; ii<current_blockcount_x; ii++){
3498	for(size_t jj=0; jj<current_blockcount_y; jj++){
3499	for(size_t kk=0; kk<current_blockcount_z; kk++){
3500	type_ = type[index];
3501	if (type_ != 0){
3502	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
3503	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3504	}
3505	else{
3506	*block_data_pos = unpred_data[unpredictable_count ++];
3507	}
3508	index ++;
3509	block_data_pos ++;
3510	}
3511	block_data_pos += dim1_offset - current_blockcount_z;
3512	}
3513	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3514	}
3515	cur_unpred_count = unpredictable_count;
3516	}
3517	}
3518	indicator_pos ++;
3519	type += current_block_elements;
3520	unpred_data += cur_unpred_count;
3521	} // end k == 0
3522	for(size_t k=1; k<num_z; k++){
3523	offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
3524	offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
3525	offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
3526	data_pos = data + offset_x dim0_offset + offset_y * dim1_offset + offset_z;
3527
3528	current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
3529	current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
3530	current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
3531
3532	size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
3533	if(*indicator_pos){
3534	// decompress by SZ
3535	double * block_data_pos = data_pos;
3536	double pred;
3537	size_t index = 0;
3538	int type_;
3539	size_t unpredictable_count = 0;
3540	for(size_t ii=0; ii<current_blockcount_x; ii++){
3541	for(size_t jj=0; jj<current_blockcount_y; jj++){
3542	for(size_t kk=0; kk<current_blockcount_z; kk++){
3543	type_ = type[index];
3544	if(type_ == intvRadius){
3545	*block_data_pos = mean;
3546	}
3547	else if(type_ == 0){
3548	*block_data_pos = unpred_data[unpredictable_count ++];
3549	}
3550	else{
3551	if(type_ < intvRadius) type_ += 1;
3552	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
3553	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3554	}
3555	index ++;
3556	block_data_pos ++;
3557	}
3558	block_data_pos += dim1_offset - current_blockcount_z;
3559	}
3560	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3561	}
3562	cur_unpred_count = unpredictable_count;
3563	}
3564	else{
3565	// decompress by regression
3566	{
3567	//restore regression coefficients
3568	double pred;
3569	int type_;
3570	for(int e=0; e<4; e++){
3571	type_ = coeff_type[e][coeff_index];
3572	if (type_ != 0){
3573	pred = last_coefficients[e];
3574	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
3575	}
3576	else{
3577	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
3578	coeff_unpred_data_count[e] ++;
3579	}
3580	}
3581	coeff_index ++;
3582	}
3583	{
3584	double * block_data_pos = data_pos;
3585	double pred;
3586	int type_;
3587	size_t index = 0;
3588	size_t unpredictable_count = 0;
3589	for(size_t ii=0; ii<current_blockcount_x; ii++){
3590	for(size_t jj=0; jj<current_blockcount_y; jj++){
3591	for(size_t kk=0; kk<current_blockcount_z; kk++){
3592	type_ = type[index];
3593	if (type_ != 0){
3594	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
3595	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3596	}
3597	else{
3598	*block_data_pos = unpred_data[unpredictable_count ++];
3599	}
3600	index ++;
3601	block_data_pos ++;
3602	}
3603	block_data_pos += dim1_offset - current_blockcount_z;
3604	}
3605	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3606	}
3607	cur_unpred_count = unpredictable_count;
3608	}
3609	}
3610	indicator_pos ++;
3611	type += current_block_elements;
3612	unpred_data += cur_unpred_count;
3613	}
3614	}
3615	}
3616	}
3617	else{
3618	type = result_type;
3619	// i == 0
3620	{
3621	// j == 0
3622	{
3623	// k == 0
3624	{
3625	data_pos = *data;
3626
3627	current_blockcount_x = early_blockcount_x;
3628	current_blockcount_y = early_blockcount_y;
3629	current_blockcount_z = early_blockcount_z;
3630	size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
3631	if(*indicator_pos){
3632	// decompress by SZ
3633	double * block_data_pos = data_pos;
3634	double pred;
3635	size_t index = 0;
3636	int type_;
3637	size_t unpredictable_count = 0;
3638	// ii == 0
3639	{
3640	// jj == 0
3641	{
3642	{
3643	// kk == 0
3644	type_ = type[index];
3645	if(type_ == 0){
3646	*block_data_pos = unpred_data[unpredictable_count ++];
3647	}
3648	else{
3649	pred = 0;
3650	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3651	}
3652	index ++;
3653	block_data_pos ++;
3654	}
3655	for(size_t kk=1; kk<current_blockcount_z; kk++){
3656	type_ = type[index];
3657	if(type_ == 0){
3658	*block_data_pos = unpred_data[unpredictable_count ++];
3659	}
3660	else{
3661	pred = block_data_pos[- 1];
3662	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3663	}
3664	index ++;
3665	block_data_pos ++;
3666	}
3667	block_data_pos += dim1_offset - current_blockcount_z;
3668	}
3669	for(size_t jj=1; jj<current_blockcount_y; jj++){
3670	{
3671	// kk == 0
3672	type_ = type[index];
3673	if(type_ == 0){
3674	*block_data_pos = unpred_data[unpredictable_count ++];
3675	}
3676	else{
3677	pred = block_data_pos[- dim1_offset];
3678	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3679	}
3680	index ++;
3681	block_data_pos ++;
3682	}
3683	for(size_t kk=1; kk<current_blockcount_z; kk++){
3684	type_ = type[index];
3685	if(type_ == 0){
3686	*block_data_pos = unpred_data[unpredictable_count ++];
3687	}
3688	else{
3689	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] - block_data_pos[- dim1_offset - 1];
3690	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3691	}
3692	index ++;
3693	block_data_pos ++;
3694	}
3695	block_data_pos += dim1_offset - current_blockcount_z;
3696	}
3697	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3698	}
3699	for(size_t ii=1; ii<current_blockcount_x; ii++){
3700	// jj == 0
3701	{
3702	{
3703	// kk == 0
3704	type_ = type[index];
3705	if(type_ == 0){
3706	*block_data_pos = unpred_data[unpredictable_count ++];
3707	}
3708	else{
3709	pred = block_data_pos[- dim0_offset];
3710	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3711	}
3712	index ++;
3713	block_data_pos ++;
3714	}
3715	for(size_t kk=1; kk<current_blockcount_z; kk++){
3716	type_ = type[index];
3717	if(type_ == 0){
3718	*block_data_pos = unpred_data[unpredictable_count ++];
3719	}
3720	else{
3721	pred = block_data_pos[- 1] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - 1];
3722	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3723	}
3724	index ++;
3725	block_data_pos ++;
3726	}
3727	block_data_pos += dim1_offset - current_blockcount_z;
3728	}
3729	for(size_t jj=1; jj<current_blockcount_y; jj++){
3730	{
3731	// kk == 0
3732	type_ = type[index];
3733	if(type_ == 0){
3734	*block_data_pos = unpred_data[unpredictable_count ++];
3735	}
3736	else{
3737	pred = block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - dim1_offset];
3738	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3739	}
3740	index ++;
3741	block_data_pos ++;
3742	}
3743	for(size_t kk=1; kk<current_blockcount_z; kk++){
3744	type_ = type[index];
3745	if(type_ == 0){
3746	*block_data_pos = unpred_data[unpredictable_count ++];
3747	}
3748	else{
3749	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
3750	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3751	}
3752	index ++;
3753	block_data_pos ++;
3754	}
3755	block_data_pos += dim1_offset - current_blockcount_z;
3756	}
3757	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3758	}
3759	cur_unpred_count = unpredictable_count;
3760	}
3761	else{
3762	// decompress by regression
3763	{
3764	//restore regression coefficients
3765	double pred;
3766	int type_;
3767	for(int e=0; e<4; e++){
3768	type_ = coeff_type[e][coeff_index];
3769	if (type_ != 0){
3770	pred = last_coefficients[e];
3771	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
3772	}
3773	else{
3774	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
3775	coeff_unpred_data_count[e] ++;
3776	}
3777	}
3778	coeff_index ++;
3779	}
3780	{
3781	double * block_data_pos = data_pos;
3782	double pred;
3783	int type_;
3784	size_t index = 0;
3785	size_t unpredictable_count = 0;
3786	for(size_t ii=0; ii<current_blockcount_x; ii++){
3787	for(size_t jj=0; jj<current_blockcount_y; jj++){
3788	for(size_t kk=0; kk<current_blockcount_z; kk++){
3789	type_ = type[index];
3790	if (type_ != 0){
3791	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
3792	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3793	}
3794	else{
3795	*block_data_pos = unpred_data[unpredictable_count ++];
3796	}
3797	index ++;
3798	block_data_pos ++;
3799	}
3800	block_data_pos += dim1_offset - current_blockcount_z;
3801	}
3802	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3803	}
3804	cur_unpred_count = unpredictable_count;
3805	}
3806	}
3807	indicator_pos ++;
3808	type += current_block_elements;
3809	unpred_data += cur_unpred_count;
3810	} // end k == 0
3811	// i == 0 j == 0 k != 0
3812	for(size_t k=1; k<num_z; k++){
3813	offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
3814	data_pos = *data + offset_z;
3815
3816	current_blockcount_x = early_blockcount_x;
3817	current_blockcount_y = early_blockcount_y;
3818	current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
3819
3820	size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
3821	if(*indicator_pos){
3822	// decompress by SZ
3823	double * block_data_pos = data_pos;
3824	double pred;
3825	size_t index = 0;
3826	int type_;
3827	size_t unpredictable_count = 0;
3828	// ii == 0
3829	{
3830	// jj == 0
3831	{
3832	for(size_t kk=0; kk<current_blockcount_z; kk++){
3833	type_ = type[index];
3834	if(type_ == 0){
3835	*block_data_pos = unpred_data[unpredictable_count ++];
3836	}
3837	else{
3838	pred = block_data_pos[- 1];
3839	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3840	}
3841	index ++;
3842	block_data_pos ++;
3843	}
3844	block_data_pos += dim1_offset - current_blockcount_z;
3845	}
3846	for(size_t jj=1; jj<current_blockcount_y; jj++){
3847	for(size_t kk=0; kk<current_blockcount_z; kk++){
3848	type_ = type[index];
3849	if(type_ == 0){
3850	*block_data_pos = unpred_data[unpredictable_count ++];
3851	}
3852	else{
3853	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] - block_data_pos[- dim1_offset - 1];
3854	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3855	}
3856	index ++;
3857	block_data_pos ++;
3858	}
3859	block_data_pos += dim1_offset - current_blockcount_z;
3860	}
3861	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3862	}
3863	for(size_t ii=1; ii<current_blockcount_x; ii++){
3864	// jj == 0
3865	{
3866	for(size_t kk=0; kk<current_blockcount_z; kk++){
3867	type_ = type[index];
3868	if(type_ == 0){
3869	*block_data_pos = unpred_data[unpredictable_count ++];
3870	}
3871	else{
3872	pred = block_data_pos[- 1] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - 1];
3873	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3874	}
3875	index ++;
3876	block_data_pos ++;
3877	}
3878	block_data_pos += dim1_offset - current_blockcount_z;
3879	}
3880	for(size_t jj=1; jj<current_blockcount_y; jj++){
3881	for(size_t kk=0; kk<current_blockcount_z; kk++){
3882	type_ = type[index];
3883	if(type_ == 0){
3884	*block_data_pos = unpred_data[unpredictable_count ++];
3885	}
3886	else{
3887	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
3888	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3889	}
3890	index ++;
3891	block_data_pos ++;
3892	}
3893	block_data_pos += dim1_offset - current_blockcount_z;
3894	}
3895	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3896	}
3897	cur_unpred_count = unpredictable_count;
3898	}
3899	else{
3900	// decompress by regression
3901	{
3902	//restore regression coefficients
3903	double pred;
3904	int type_;
3905	for(int e=0; e<4; e++){
3906	type_ = coeff_type[e][coeff_index];
3907	if (type_ != 0){
3908	pred = last_coefficients[e];
3909	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
3910	}
3911	else{
3912	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
3913	coeff_unpred_data_count[e] ++;
3914	}
3915	}
3916	coeff_index ++;
3917	}
3918	{
3919	double * block_data_pos = data_pos;
3920	double pred;
3921	int type_;
3922	size_t index = 0;
3923	size_t unpredictable_count = 0;
3924	for(size_t ii=0; ii<current_blockcount_x; ii++){
3925	for(size_t jj=0; jj<current_blockcount_y; jj++){
3926	for(size_t kk=0; kk<current_blockcount_z; kk++){
3927	type_ = type[index];
3928	if (type_ != 0){
3929	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
3930	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3931	}
3932	else{
3933	*block_data_pos = unpred_data[unpredictable_count ++];
3934	}
3935	index ++;
3936	block_data_pos ++;
3937	}
3938	block_data_pos += dim1_offset - current_blockcount_z;
3939	}
3940	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3941	}
3942	cur_unpred_count = unpredictable_count;
3943	}
3944	}
3945	indicator_pos ++;
3946	type += current_block_elements;
3947	unpred_data += cur_unpred_count;
3948	}
3949	}// end j==0
3950	for(size_t j=1; j<num_y; j++){
3951	// k == 0
3952	{
3953	offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
3954	data_pos = data + offset_y dim1_offset;
3955
3956	current_blockcount_x = early_blockcount_x;
3957	current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
3958	current_blockcount_z = early_blockcount_z;
3959	size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
3960	if(*indicator_pos){
3961	// decompress by SZ
3962	double * block_data_pos = data_pos;
3963	double pred;
3964	size_t index = 0;
3965	int type_;
3966	size_t unpredictable_count = 0;
3967	// ii == 0
3968	{
3969	for(size_t jj=0; jj<current_blockcount_y; jj++){
3970	{
3971	// kk == 0
3972	type_ = type[index];
3973	if(type_ == 0){
3974	*block_data_pos = unpred_data[unpredictable_count ++];
3975	}
3976	else{
3977	pred = block_data_pos[- dim1_offset];
3978	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3979	}
3980	index ++;
3981	block_data_pos ++;
3982	}
3983	for(size_t kk=1; kk<current_blockcount_z; kk++){
3984	type_ = type[index];
3985	if(type_ == 0){
3986	*block_data_pos = unpred_data[unpredictable_count ++];
3987	}
3988	else{
3989	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] - block_data_pos[- dim1_offset - 1];
3990	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
3991	}
3992	index ++;
3993	block_data_pos ++;
3994	}
3995	block_data_pos += dim1_offset - current_blockcount_z;
3996	}
3997	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3998	}
3999	for(size_t ii=1; ii<current_blockcount_x; ii++){
4000	for(size_t jj=0; jj<current_blockcount_y; jj++){
4001	{
4002	// kk == 0
4003	type_ = type[index];
4004	if(type_ == 0){
4005	*block_data_pos = unpred_data[unpredictable_count ++];
4006	}
4007	else{
4008	pred = block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - dim1_offset];
4009	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4010	}
4011	index ++;
4012	block_data_pos ++;
4013	}
4014	for(size_t kk=1; kk<current_blockcount_z; kk++){
4015	type_ = type[index];
4016	if(type_ == 0){
4017	*block_data_pos = unpred_data[unpredictable_count ++];
4018	}
4019	else{
4020	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
4021	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4022	}
4023	index ++;
4024	block_data_pos ++;
4025	}
4026	block_data_pos += dim1_offset - current_blockcount_z;
4027	}
4028	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4029	}
4030	cur_unpred_count = unpredictable_count;
4031	}
4032	else{
4033	// decompress by regression
4034	{
4035	//restore regression coefficients
4036	double pred;
4037	int type_;
4038	for(int e=0; e<4; e++){
4039	type_ = coeff_type[e][coeff_index];
4040	if (type_ != 0){
4041	pred = last_coefficients[e];
4042	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
4043	}
4044	else{
4045	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
4046	coeff_unpred_data_count[e] ++;
4047	}
4048	}
4049	coeff_index ++;
4050	}
4051	{
4052	double * block_data_pos = data_pos;
4053	double pred;
4054	int type_;
4055	size_t index = 0;
4056	size_t unpredictable_count = 0;
4057	for(size_t ii=0; ii<current_blockcount_x; ii++){
4058	for(size_t jj=0; jj<current_blockcount_y; jj++){
4059	for(size_t kk=0; kk<current_blockcount_z; kk++){
4060	type_ = type[index];
4061	if (type_ != 0){
4062	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
4063	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4064	}
4065	else{
4066	*block_data_pos = unpred_data[unpredictable_count ++];
4067	}
4068	index ++;
4069	block_data_pos ++;
4070	}
4071	block_data_pos += dim1_offset - current_blockcount_z;
4072	}
4073	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4074	}
4075	cur_unpred_count = unpredictable_count;
4076	}
4077	}
4078	indicator_pos ++;
4079	type += current_block_elements;
4080	unpred_data += cur_unpred_count;
4081	} // end k == 0
4082	for(size_t k=1; k<num_z; k++){
4083	offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
4084	offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
4085	data_pos = data + offset_y dim1_offset + offset_z;
4086
4087	current_blockcount_x = early_blockcount_x;
4088	current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
4089	current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
4090
4091	size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
4092	if(*indicator_pos){
4093	// decompress by SZ
4094	double * block_data_pos = data_pos;
4095	double pred;
4096	size_t index = 0;
4097	int type_;
4098	size_t unpredictable_count = 0;
4099	// ii == 0
4100	{
4101	for(size_t jj=0; jj<current_blockcount_y; jj++){
4102	for(size_t kk=0; kk<current_blockcount_z; kk++){
4103	type_ = type[index];
4104	if(type_ == 0){
4105	*block_data_pos = unpred_data[unpredictable_count ++];
4106	}
4107	else{
4108	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] - block_data_pos[- dim1_offset - 1];
4109	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4110	}
4111	index ++;
4112	block_data_pos ++;
4113	}
4114	block_data_pos += dim1_offset - current_blockcount_z;
4115	}
4116	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4117	}
4118	for(size_t ii=1; ii<current_blockcount_x; ii++){
4119	for(size_t jj=0; jj<current_blockcount_y; jj++){
4120	for(size_t kk=0; kk<current_blockcount_z; kk++){
4121	type_ = type[index];
4122	if(type_ == 0){
4123	*block_data_pos = unpred_data[unpredictable_count ++];
4124	}
4125	else{
4126	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
4127	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4128	}
4129	index ++;
4130	block_data_pos ++;
4131	}
4132	block_data_pos += dim1_offset - current_blockcount_z;
4133	}
4134	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4135	}
4136	cur_unpred_count = unpredictable_count;
4137	}
4138	else{
4139	// decompress by regression
4140	{
4141	//restore regression coefficients
4142	double pred;
4143	int type_;
4144	for(int e=0; e<4; e++){
4145	type_ = coeff_type[e][coeff_index];
4146	if (type_ != 0){
4147	pred = last_coefficients[e];
4148	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
4149	}
4150	else{
4151	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
4152	coeff_unpred_data_count[e] ++;
4153	}
4154	}
4155	coeff_index ++;
4156	}
4157	{
4158	double * block_data_pos = data_pos;
4159	double pred;
4160	int type_;
4161	size_t index = 0;
4162	size_t unpredictable_count = 0;
4163	for(size_t ii=0; ii<current_blockcount_x; ii++){
4164	for(size_t jj=0; jj<current_blockcount_y; jj++){
4165	for(size_t kk=0; kk<current_blockcount_z; kk++){
4166	type_ = type[index];
4167	if (type_ != 0){
4168	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
4169	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4170	}
4171	else{
4172	*block_data_pos = unpred_data[unpredictable_count ++];
4173	}
4174	index ++;
4175	block_data_pos ++;
4176	}
4177	block_data_pos += dim1_offset - current_blockcount_z;
4178	}
4179	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4180	}
4181	cur_unpred_count = unpredictable_count;
4182	}
4183	}
4184	indicator_pos ++;
4185	type += current_block_elements;
4186	unpred_data += cur_unpred_count;
4187	}
4188	}
4189	} // end i==0
4190	for(size_t i=1; i<num_x; i++){
4191	// j == 0
4192	{
4193	// k == 0
4194	{
4195	offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
4196	data_pos = data + offset_x dim0_offset;
4197
4198	current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
4199	current_blockcount_y = early_blockcount_y;
4200	current_blockcount_z = early_blockcount_z;
4201	size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
4202	if(*indicator_pos){
4203	// decompress by SZ
4204	double * block_data_pos = data_pos;
4205	double pred;
4206	size_t index = 0;
4207	int type_;
4208	size_t unpredictable_count = 0;
4209	for(size_t ii=0; ii<current_blockcount_x; ii++){
4210	// jj == 0
4211	{
4212	{
4213	// kk == 0
4214	type_ = type[index];
4215	if(type_ == 0){
4216	*block_data_pos = unpred_data[unpredictable_count ++];
4217	}
4218	else{
4219	pred = block_data_pos[- dim0_offset];
4220	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4221	}
4222	index ++;
4223	block_data_pos ++;
4224	}
4225	for(size_t kk=1; kk<current_blockcount_z; kk++){
4226	type_ = type[index];
4227	if(type_ == 0){
4228	*block_data_pos = unpred_data[unpredictable_count ++];
4229	}
4230	else{
4231	pred = block_data_pos[- 1] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - 1];
4232	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4233	}
4234	index ++;
4235	block_data_pos ++;
4236	}
4237	block_data_pos += dim1_offset - current_blockcount_z;
4238	}
4239	for(size_t jj=1; jj<current_blockcount_y; jj++){
4240	{
4241	// kk == 0
4242	type_ = type[index];
4243	if(type_ == 0){
4244	*block_data_pos = unpred_data[unpredictable_count ++];
4245	}
4246	else{
4247	pred = block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - dim1_offset];
4248	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4249	}
4250	index ++;
4251	block_data_pos ++;
4252	}
4253	for(size_t kk=1; kk<current_blockcount_z; kk++){
4254	type_ = type[index];
4255	if(type_ == 0){
4256	*block_data_pos = unpred_data[unpredictable_count ++];
4257	}
4258	else{
4259	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
4260	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4261	}
4262	index ++;
4263	block_data_pos ++;
4264	}
4265	block_data_pos += dim1_offset - current_blockcount_z;
4266	}
4267	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4268	}
4269	cur_unpred_count = unpredictable_count;
4270	}
4271	else{
4272	// decompress by regression
4273	{
4274	//restore regression coefficients
4275	double pred;
4276	int type_;
4277	for(int e=0; e<4; e++){
4278	type_ = coeff_type[e][coeff_index];
4279	if (type_ != 0){
4280	pred = last_coefficients[e];
4281	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
4282	}
4283	else{
4284	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
4285	coeff_unpred_data_count[e] ++;
4286	}
4287	}
4288	coeff_index ++;
4289	}
4290	{
4291	double * block_data_pos = data_pos;
4292	double pred;
4293	int type_;
4294	size_t index = 0;
4295	size_t unpredictable_count = 0;
4296	for(size_t ii=0; ii<current_blockcount_x; ii++){
4297	for(size_t jj=0; jj<current_blockcount_y; jj++){
4298	for(size_t kk=0; kk<current_blockcount_z; kk++){
4299	type_ = type[index];
4300	if (type_ != 0){
4301	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
4302	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4303	}
4304	else{
4305	*block_data_pos = unpred_data[unpredictable_count ++];
4306	}
4307	index ++;
4308	block_data_pos ++;
4309	}
4310	block_data_pos += dim1_offset - current_blockcount_z;
4311	}
4312	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4313	}
4314	cur_unpred_count = unpredictable_count;
4315	}
4316	}
4317	indicator_pos ++;
4318	type += current_block_elements;
4319	unpred_data += cur_unpred_count;
4320	} // end k == 0
4321	for(size_t k=1; k<num_z; k++){
4322	offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
4323	offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
4324	data_pos = data + offset_x dim0_offset + offset_z;
4325
4326	current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
4327	current_blockcount_y = early_blockcount_y;
4328	current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
4329
4330	size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
4331	if(*indicator_pos){
4332	// decompress by SZ
4333	double * block_data_pos = data_pos;
4334	double pred;
4335	size_t index = 0;
4336	int type_;
4337	size_t unpredictable_count = 0;
4338	for(size_t ii=0; ii<current_blockcount_x; ii++){
4339	// jj == 0
4340	{
4341	for(size_t kk=0; kk<current_blockcount_z; kk++){
4342	type_ = type[index];
4343	if(type_ == 0){
4344	*block_data_pos = unpred_data[unpredictable_count ++];
4345	}
4346	else{
4347	pred = block_data_pos[- 1] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - 1];
4348	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4349	}
4350	index ++;
4351	block_data_pos ++;
4352	}
4353	block_data_pos += dim1_offset - current_blockcount_z;
4354	}
4355	for(size_t jj=1; jj<current_blockcount_y; jj++){
4356	for(size_t kk=0; kk<current_blockcount_z; kk++){
4357	type_ = type[index];
4358	if(type_ == 0){
4359	*block_data_pos = unpred_data[unpredictable_count ++];
4360	}
4361	else{
4362	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
4363	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4364	}
4365	index ++;
4366	block_data_pos ++;
4367	}
4368	block_data_pos += dim1_offset - current_blockcount_z;
4369	}
4370	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4371	}
4372	cur_unpred_count = unpredictable_count;
4373	}
4374	else{
4375	// decompress by regression
4376	{
4377	//restore regression coefficients
4378	double pred;
4379	int type_;
4380	for(int e=0; e<4; e++){
4381	type_ = coeff_type[e][coeff_index];
4382	if (type_ != 0){
4383	pred = last_coefficients[e];
4384	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
4385	}
4386	else{
4387	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
4388	coeff_unpred_data_count[e] ++;
4389	}
4390	}
4391	coeff_index ++;
4392	}
4393	{
4394	double * block_data_pos = data_pos;
4395	double pred;
4396	int type_;
4397	size_t index = 0;
4398	size_t unpredictable_count = 0;
4399	for(size_t ii=0; ii<current_blockcount_x; ii++){
4400	for(size_t jj=0; jj<current_blockcount_y; jj++){
4401	for(size_t kk=0; kk<current_blockcount_z; kk++){
4402	type_ = type[index];
4403	if (type_ != 0){
4404	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
4405	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4406	}
4407	else{
4408	*block_data_pos = unpred_data[unpredictable_count ++];
4409	}
4410	index ++;
4411	block_data_pos ++;
4412	}
4413	block_data_pos += dim1_offset - current_blockcount_z;
4414	}
4415	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4416	}
4417	cur_unpred_count = unpredictable_count;
4418	}
4419	}
4420	indicator_pos ++;
4421	type += current_block_elements;
4422	unpred_data += cur_unpred_count;
4423	}
4424	}// end j = 0
4425	for(size_t j=1; j<num_y; j++){
4426	// k == 0
4427	{
4428	offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
4429	offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
4430	data_pos = data + offset_x dim0_offset + offset_y * dim1_offset;
4431
4432	current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
4433	current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
4434	current_blockcount_z = early_blockcount_z;
4435	size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
4436	if(*indicator_pos){
4437	// decompress by SZ
4438	double * block_data_pos = data_pos;
4439	double pred;
4440	size_t index = 0;
4441	int type_;
4442	size_t unpredictable_count = 0;
4443	for(size_t ii=0; ii<current_blockcount_x; ii++){
4444	for(size_t jj=0; jj<current_blockcount_y; jj++){
4445	{
4446	// kk == 0
4447	type_ = type[index];
4448	if(type_ == 0){
4449	*block_data_pos = unpred_data[unpredictable_count ++];
4450	}
4451	else{
4452	pred = block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - dim1_offset];
4453	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4454	}
4455	index ++;
4456	block_data_pos ++;
4457	}
4458	for(size_t kk=1; kk<current_blockcount_z; kk++){
4459	type_ = type[index];
4460	if(type_ == 0){
4461	*block_data_pos = unpred_data[unpredictable_count ++];
4462	}
4463	else{
4464	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
4465	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4466	}
4467	index ++;
4468	block_data_pos ++;
4469	}
4470	block_data_pos += dim1_offset - current_blockcount_z;
4471	}
4472	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4473	}
4474	cur_unpred_count = unpredictable_count;
4475	}
4476	else{
4477	// decompress by regression
4478	{
4479	//restore regression coefficients
4480	double pred;
4481	int type_;
4482	for(int e=0; e<4; e++){
4483	type_ = coeff_type[e][coeff_index];
4484	if (type_ != 0){
4485	pred = last_coefficients[e];
4486	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
4487	}
4488	else{
4489	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
4490	coeff_unpred_data_count[e] ++;
4491	}
4492	}
4493	coeff_index ++;
4494	}
4495	{
4496	double * block_data_pos = data_pos;
4497	double pred;
4498	int type_;
4499	size_t index = 0;
4500	size_t unpredictable_count = 0;
4501	for(size_t ii=0; ii<current_blockcount_x; ii++){
4502	for(size_t jj=0; jj<current_blockcount_y; jj++){
4503	for(size_t kk=0; kk<current_blockcount_z; kk++){
4504	type_ = type[index];
4505	if (type_ != 0){
4506	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
4507	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4508	}
4509	else{
4510	*block_data_pos = unpred_data[unpredictable_count ++];
4511	}
4512	index ++;
4513	block_data_pos ++;
4514	}
4515	block_data_pos += dim1_offset - current_blockcount_z;
4516	}
4517	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4518	}
4519	cur_unpred_count = unpredictable_count;
4520	}
4521	}
4522	indicator_pos ++;
4523	type += current_block_elements;
4524	unpred_data += cur_unpred_count;
4525	} // end k == 0
4526	for(size_t k=1; k<num_z; k++){
4527	offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
4528	offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
4529	offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
4530	data_pos = data + offset_x dim0_offset + offset_y * dim1_offset + offset_z;
4531
4532	current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
4533	current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
4534	current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
4535
4536	size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
4537	if(*indicator_pos){
4538	// decompress by SZ
4539	double * block_data_pos = data_pos;
4540	double pred;
4541	size_t index = 0;
4542	int type_;
4543	size_t unpredictable_count = 0;
4544	for(size_t ii=0; ii<current_blockcount_x; ii++){
4545	for(size_t jj=0; jj<current_blockcount_y; jj++){
4546	for(size_t kk=0; kk<current_blockcount_z; kk++){
4547	type_ = type[index];
4548	if(type_ == 0){
4549	*block_data_pos = unpred_data[unpredictable_count ++];
4550	}
4551	else{
4552	pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
4553	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4554	}
4555	index ++;
4556	block_data_pos ++;
4557	}
4558	block_data_pos += dim1_offset - current_blockcount_z;
4559	}
4560	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4561	}
4562	cur_unpred_count = unpredictable_count;
4563	}
4564	else{
4565	// decompress by regression
4566	{
4567	//restore regression coefficients
4568	double pred;
4569	int type_;
4570	for(int e=0; e<4; e++){
4571	type_ = coeff_type[e][coeff_index];
4572	if (type_ != 0){
4573	pred = last_coefficients[e];
4574	last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
4575	}
4576	else{
4577	last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
4578	coeff_unpred_data_count[e] ++;
4579	}
4580	}
4581	coeff_index ++;
4582	}
4583	{
4584	double * block_data_pos = data_pos;
4585	double pred;
4586	int type_;
4587	size_t index = 0;
4588	size_t unpredictable_count = 0;
4589	for(size_t ii=0; ii<current_blockcount_x; ii++){
4590	for(size_t jj=0; jj<current_blockcount_y; jj++){
4591	for(size_t kk=0; kk<current_blockcount_z; kk++){
4592	type_ = type[index];
4593	if (type_ != 0){
4594	pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
4595	block_data_pos = pred + 2 (type_ - intvRadius) * realPrecision;
4596	}
4597	else{
4598	*block_data_pos = unpred_data[unpredictable_count ++];
4599	}
4600	index ++;
4601	block_data_pos ++;
4602	}
4603	block_data_pos += dim1_offset - current_blockcount_z;
4604	}
4605	block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4606	}
4607	cur_unpred_count = unpredictable_count;
4608	}
4609	}
4610	indicator_pos ++;
4611	type += current_block_elements;
4612	unpred_data += cur_unpred_count;
4613	}
4614	}
4615	}
4616	}
4617
4618	free(coeff_result_type);
4619
4620	free(indicator);
4621	free(result_type);
4622	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: