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`

Rev	Line
[2c47b73]	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"
[9ee2ce3]	19	#include "utility.h"
[2c47b73]	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	{
[9ee2ce3]	33	confparams_dec->losslessCompressor = is_lossless_compressed_data(cmpBytes, cmpSize);
[2c47b73]	34	if(confparams_dec->szMode!=SZ_TEMPORAL_COMPRESSION)
	35	{
[9ee2ce3]	36	if(confparams_dec->losslessCompressor!=-1)
[2c47b73]	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;
[9ee2ce3]	50	tmpSize = sz_lossless_decompress(confparams_dec->losslessCompressor, cmpBytes, (unsigned long)cmpSize, &szTmpBytes, (unsigned long)targetUncompressSize+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE);
[2c47b73]	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	}
[9ee2ce3]	84	else
[2c47b73]	85	{
[9ee2ce3]	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
[2c47b73]	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	{
[9ee2ce3]	1674	decompressDataSeries_double_1D_pwr_pre_log(data, dataSeriesLength, tdps);
	1675	//decompressDataSeries_double_1D_pwrgroup(data, dataSeriesLength, tdps);
[2c47b73]	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
[9ee2ce3]	1698	//decompressDataSeries_double_1D_pwr(&decmpData, dataSeriesLength, tdps);
	1699	decompressDataSeries_double_1D_pwr_pre_log(&decmpData, dataSeriesLength, tdps);
[2c47b73]	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
[9ee2ce3]	1739	//decompressDataSeries_double_2D_pwr(data, r1, r2, tdps);
	1740	decompressDataSeries_double_2D_pwr_pre_log(data, r1, r2, tdps);
[2c47b73]	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
[9ee2ce3]	1762	//decompressDataSeries_double_2D_pwr(&decmpData, r1, r2, tdps);
	1763	decompressDataSeries_double_2D_pwr_pre_log(&decmpData, r1, r2, tdps);
[2c47b73]	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
[9ee2ce3]	1803	//decompressDataSeries_double_3D_pwr(data, r1, r2, r3, tdps);
	1804	decompressDataSeries_double_3D_pwr_pre_log(data, r1, r2, r3, tdps);
[2c47b73]	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
[9ee2ce3]	1826	//decompressDataSeries_double_3D_pwr(&decmpData, r1, r2, r3, tdps);
	1827	decompressDataSeries_double_3D_pwr_pre_log(&decmpData, r1, r2, r3, tdps);
[2c47b73]	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	{
[9ee2ce3]	1868	//decompressDataSeries_double_3D_pwr(data, r1*r2, r3, r4, tdps);
	1869	decompressDataSeries_double_3D_pwr_pre_log(data, r1*r2, r3, r4, tdps);
[2c47b73]	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
[9ee2ce3]	1891	//decompressDataSeries_double_3D_pwr(&decmpData, r1*r2, r3, r4, tdps);
	1892	decompressDataSeries_double_3D_pwr_pre_log(&decmpData, r1*r2, r3, r4, tdps);
[2c47b73]	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	}
[9ee2ce3]	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: