Context Navigation

source: thirdparty/SZ/sz/src/sz_double.c @ 2c47b73

Revision 2c47b73, 104.6 KB checked in by Hal Finkel <hfinkel@…>, 6 years ago (diff)
more work on adding SZ (latest version)
Property mode set to `100644`

Line
1	/**
2	* @file sz_double.c
3	* @author Sheng Di and Dingwen Tao
4	* @date Aug, 2016
5	* @brief SZ_Init, Compression and Decompression functions
6	* (C) 2016 by Mathematics and Computer Science (MCS), Argonne National Laboratory.
7	* See COPYRIGHT in top-level directory.
8	*/
9
10
11	#include <stdio.h>
12	#include <stdlib.h>
13	#include <string.h>
14	#include <unistd.h>
15	#include <math.h>
16	#include "sz.h"
17	#include "CompressElement.h"
18	#include "DynamicByteArray.h"
19	#include "DynamicIntArray.h"
20	#include "TightDataPointStorageD.h"
21	#include "sz_double.h"
22	#include "sz_double_pwr.h"
23	#include "szd_double.h"
24	#include "szd_double_pwr.h"
25	#include "zlib.h"
26	#include "rw.h"
27	#include "sz_double_ts.h"
28
29	unsigned char* SZ_skip_compress_double(double* data, size_t dataLength, size_t* outSize)
30	{
31	outSize = dataLengthsizeof(double);
32	unsigned char* out = (unsigned char)malloc(dataLengthsizeof(double));
33	memcpy(out, data, dataLength*sizeof(double));
34	return out;
35	}
36
37	void computeReqLength_double(double realPrecision, short radExpo, int* reqLength, double* medianValue)
38	{
39	short reqExpo = getPrecisionReqLength_double(realPrecision);
40	*reqLength = 12+radExpo - reqExpo; //radExpo-reqExpo == reqMantiLength
41	if(*reqLength<12)
42	*reqLength = 12;
43	if(*reqLength>64)
44	{
45	*reqLength = 64;
46	*medianValue = 0;
47	}
48	}
49
50	unsigned int optimize_intervals_double_1D(double *oriData, size_t dataLength, double realPrecision)
51	{
52	size_t i = 0, radiusIndex;
53	double pred_value = 0, pred_err;
54	size_t intervals = (size_t)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
55	memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
56	size_t totalSampleSize = dataLength/confparams_cpr->sampleDistance;
57	for(i=2;i<dataLength;i++)
58	{
59	if(i%confparams_cpr->sampleDistance==0)
60	{
61	//pred_value = 2*oriData[i-1] - oriData[i-2];
62	pred_value = oriData[i-1];
63	pred_err = fabs(pred_value - oriData[i]);
64	radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
65	if(radiusIndex>=confparams_cpr->maxRangeRadius)
66	radiusIndex = confparams_cpr->maxRangeRadius - 1;
67	intervals[radiusIndex]++;
68	}
69	}
70	//compute the appropriate number
71	size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
72	size_t sum = 0;
73	for(i=0;i<confparams_cpr->maxRangeRadius;i++)
74	{
75	sum += intervals[i];
76	if(sum>targetCount)
77	break;
78	}
79
80	if(i>=confparams_cpr->maxRangeRadius)
81	i = confparams_cpr->maxRangeRadius-1;
82	unsigned int accIntervals = 2*(i+1);
83	unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
84
85	if(powerOf2<32)
86	powerOf2 = 32;
87
88	free(intervals);
89	//printf("accIntervals=%d, powerOf2=%d\n", accIntervals, powerOf2);
90	return powerOf2;
91	}
92
93	unsigned int optimize_intervals_double_2D(double *oriData, size_t r1, size_t r2, double realPrecision)
94	{
95	size_t i,j, index;
96	size_t radiusIndex;
97	double pred_value = 0, pred_err;
98	size_t intervals = (size_t)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
99	memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
100	size_t totalSampleSize = (r1-1)*(r2-1)/confparams_cpr->sampleDistance;
101	for(i=1;i<r1;i++)
102	{
103	for(j=1;j<r2;j++)
104	{
105	if((i+j)%confparams_cpr->sampleDistance==0)
106	{
107	index = i*r2+j;
108	pred_value = oriData[index-1] + oriData[index-r2] - oriData[index-r2-1];
109	pred_err = fabs(pred_value - oriData[index]);
110	radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
111	if(radiusIndex>=confparams_cpr->maxRangeRadius)
112	radiusIndex = confparams_cpr->maxRangeRadius - 1;
113	intervals[radiusIndex]++;
114	}
115	}
116	}
117	//compute the appropriate number
118	size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
119	size_t sum = 0;
120	for(i=0;i<confparams_cpr->maxRangeRadius;i++)
121	{
122	sum += intervals[i];
123	if(sum>targetCount)
124	break;
125	}
126	if(i>=confparams_cpr->maxRangeRadius)
127	i = confparams_cpr->maxRangeRadius-1;
128	unsigned int accIntervals = 2*(i+1);
129	unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
130	//printf("confparams_cpr->maxRangeRadius = %d, accIntervals=%d, powerOf2=%d\n", confparams_cpr->maxRangeRadius, accIntervals, powerOf2);
131
132	if(powerOf2<32)
133	powerOf2 = 32;
134
135	free(intervals);
136	return powerOf2;
137	}
138
139	unsigned int optimize_intervals_double_3D(double *oriData, size_t r1, size_t r2, size_t r3, double realPrecision)
140	{
141	size_t i,j,k, index;
142	size_t radiusIndex;
143	size_t r23=r2*r3;
144	double pred_value = 0, pred_err;
145	size_t intervals = (size_t)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
146	memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
147	size_t totalSampleSize = (r1-1)(r2-1)(r3-1)/confparams_cpr->sampleDistance;
148	for(i=1;i<r1;i++)
149	{
150	for(j=1;j<r2;j++)
151	{
152	for(k=1;k<r3;k++)
153	{
154	if((i+j+k)%confparams_cpr->sampleDistance==0)
155	{
156	index = ir23+jr3+k;
157	pred_value = oriData[index-1] + oriData[index-r3] + oriData[index-r23]
158	- oriData[index-1-r23] - oriData[index-r3-1] - oriData[index-r3-r23] + oriData[index-r3-r23-1];
159	pred_err = fabs(pred_value - oriData[index]);
160	radiusIndex = (pred_err/realPrecision+1)/2;
161	if(radiusIndex>=confparams_cpr->maxRangeRadius)
162	radiusIndex = confparams_cpr->maxRangeRadius - 1;
163	intervals[radiusIndex]++;
164	}
165	}
166
167	}
168	}
169	//compute the appropriate number
170	size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
171	size_t sum = 0;
172	for(i=0;i<confparams_cpr->maxRangeRadius;i++)
173	{
174	sum += intervals[i];
175	if(sum>targetCount)
176	break;
177	}
178	if(i>=confparams_cpr->maxRangeRadius)
179	i = confparams_cpr->maxRangeRadius-1;
180
181	unsigned int accIntervals = 2*(i+1);
182	unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
183
184	if(powerOf2<32)
185	powerOf2 = 32;
186
187	free(intervals);
188	//printf("confparams_cpr->maxRangeRadius = %d, accIntervals=%d, powerOf2=%d\n", confparams_cpr->maxRangeRadius, accIntervals, powerOf2);
189	return powerOf2;
190	}
191
192	unsigned int optimize_intervals_double_4D(double *oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision)
193	{
194	size_t i,j,k,l, index;
195	size_t radiusIndex;
196	size_t r234=r2r3r4;
197	size_t r34=r3*r4;
198	double pred_value = 0, pred_err;
199	size_t intervals = (size_t)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
200	memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
201	size_t totalSampleSize = (r1-1)(r2-1)(r3-1)*(r4-1)/confparams_cpr->sampleDistance;
202	for(i=1;i<r1;i++)
203	{
204	for(j=1;j<r2;j++)
205	{
206	for(k=1;k<r3;k++)
207	{
208	for (l=1;l<r4;l++)
209	{
210	if((i+j+k+l)%confparams_cpr->sampleDistance==0)
211	{
212	index = ir234+jr34+k*r4+l;
213	pred_value = oriData[index-1] + oriData[index-r3] + oriData[index-r34]
214	- oriData[index-1-r34] - oriData[index-r4-1] - oriData[index-r4-r34] + oriData[index-r4-r34-1];
215	pred_err = fabs(pred_value - oriData[index]);
216	radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
217	if(radiusIndex>=confparams_cpr->maxRangeRadius)
218	radiusIndex = confparams_cpr->maxRangeRadius - 1;
219	intervals[radiusIndex]++;
220	}
221	}
222	}
223	}
224	}
225	//compute the appropriate number
226	size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
227	size_t sum = 0;
228	for(i=0;i<confparams_cpr->maxRangeRadius;i++)
229	{
230	sum += intervals[i];
231	if(sum>targetCount)
232	break;
233	}
234	if(i>=confparams_cpr->maxRangeRadius)
235	i = confparams_cpr->maxRangeRadius-1;
236
237	unsigned int accIntervals = 2*(i+1);
238	unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
239
240	if(powerOf2<32)
241	powerOf2 = 32;
242
243	free(intervals);
244	return powerOf2;
245	}
246
247	TightDataPointStorageD* SZ_compress_double_1D_MDQ(double *oriData,
248	size_t dataLength, double realPrecision, double valueRangeSize, double medianValue_d)
249	{
250	#ifdef HAVE_TIMECMPR
251	double* decData = NULL;
252	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
253	decData = (double*)(multisteps->hist_data);
254	#endif
255
256	unsigned int quantization_intervals;
257	if(exe_params->optQuantMode==1)
258	quantization_intervals = optimize_intervals_double_1D_opt(oriData, dataLength, realPrecision);
259	else
260	quantization_intervals = exe_params->intvCapacity;
261	updateQuantizationInfo(quantization_intervals);
262
263	size_t i;
264	int reqLength;
265	double medianValue = medianValue_d;
266	short radExpo = getExponent_double(valueRangeSize/2);
267
268	computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
269
270	int* type = (int) malloc(dataLengthsizeof(int));
271
272	double* spaceFillingValue = oriData; //
273
274	DynamicIntArray *exactLeadNumArray;
275	new_DIA(&exactLeadNumArray, DynArrayInitLen);
276
277	DynamicByteArray *exactMidByteArray;
278	new_DBA(&exactMidByteArray, DynArrayInitLen);
279
280	DynamicIntArray *resiBitArray;
281	new_DIA(&resiBitArray, DynArrayInitLen);
282
283	unsigned char preDataBytes[8];
284	longToBytes_bigEndian(preDataBytes, 0);
285
286	int reqBytesLength = reqLength/8;
287	int resiBitsLength = reqLength%8;
288	double last3CmprsData[3] = {0};
289
290	DoubleValueCompressElement vce = (DoubleValueCompressElement)malloc(sizeof(DoubleValueCompressElement));
291	LossyCompressionElement lce = (LossyCompressionElement)malloc(sizeof(LossyCompressionElement));
292
293	//add the first data
294	type[0] = 0;
295	compressSingleDoubleValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
296	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
297	memcpy(preDataBytes,vce->curBytes,8);
298	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
299	listAdd_double(last3CmprsData, vce->data);
300	#ifdef HAVE_TIMECMPR
301	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
302	decData[0] = vce->data;
303	#endif
304
305	//add the second data
306	type[1] = 0;
307	compressSingleDoubleValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
308	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
309	memcpy(preDataBytes,vce->curBytes,8);
310	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
311	listAdd_double(last3CmprsData, vce->data);
312	#ifdef HAVE_TIMECMPR
313	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
314	decData[1] = vce->data;
315	#endif
316	int state;
317	double checkRadius;
318	double curData;
319	double pred;
320	double predAbsErr;
321	checkRadius = (exe_params->intvCapacity-1)*realPrecision;
322	double interval = 2*realPrecision;
323
324	for(i=2;i<dataLength;i++)
325	{
326	//printf("%.30G\n",last3CmprsData[0]);
327	curData = spaceFillingValue[i];
328	//pred = 2*last3CmprsData[0] - last3CmprsData[1];
329	pred = last3CmprsData[0];
330	predAbsErr = fabs(curData - pred);
331	if(predAbsErr<=checkRadius)
332	{
333	state = (predAbsErr/realPrecision+1)/2;
334	if(curData>=pred)
335	{
336	type[i] = exe_params->intvRadius+state;
337	pred = pred + state*interval;
338	}
339	else //curData<pred
340	{
341	type[i] = exe_params->intvRadius-state;
342	pred = pred - state*interval;
343	}
344	listAdd_double(last3CmprsData, pred);
345	#ifdef HAVE_TIMECMPR
346	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
347	decData[i] = pred;
348	#endif
349	continue;
350	}
351
352	//unpredictable data processing
353	type[i] = 0;
354	compressSingleDoubleValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
355	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
356	memcpy(preDataBytes,vce->curBytes,8);
357	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
358
359	listAdd_double(last3CmprsData, vce->data);
360	#ifdef HAVE_TIMECMPR
361	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
362	decData[i] = vce->data;
363	#endif
364
365	}//end of for
366
367	int exactDataNum = exactLeadNumArray->size;
368
369	TightDataPointStorageD* tdps;
370
371	new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
372	type, exactMidByteArray->array, exactMidByteArray->size,
373	exactLeadNumArray->array,
374	resiBitArray->array, resiBitArray->size,
375	resiBitsLength,
376	realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
377
378	// printf("exactDataNum=%d, expSegmentsInBytes_size=%d, exactMidByteArray->size=%d\n",
379	// exactDataNum, expSegmentsInBytes_size, exactMidByteArray->size);
380
381	//free memory
382	free_DIA(exactLeadNumArray);
383	free_DIA(resiBitArray);
384	free(type);
385	free(vce);
386	free(lce);
387	free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
388
389	return tdps;
390	}
391
392	void SZ_compress_args_double_StoreOriData(double* oriData, size_t dataLength, TightDataPointStorageD* tdps,
393	unsigned char** newByteData, size_t *outSize)
394	{
395	int doubleSize = sizeof(double);
396	size_t k = 0, i;
397	tdps->isLossless = 1;
398	size_t totalByteLength = 3 + MetaDataByteLength + exe_params->SZ_SIZE_TYPE + 1 + doubleSize*dataLength;
399	newByteData = (unsigned char)malloc(totalByteLength);
400
401	unsigned char dsLengthBytes[8];
402	for (i = 0; i < 3; i++)//3
403	(*newByteData)[k++] = versionNumber[i];
404
405	if(exe_params->SZ_SIZE_TYPE==4)//1
406	(*newByteData)[k++] = 16; //00010000
407	else
408	(*newByteData)[k++] = 80; //01010000: 01000000 indicates the SZ_SIZE_TYPE=8
409
410	convertSZParamsToBytes(confparams_cpr, &((*newByteData)[k]));
411	k = k + MetaDataByteLength;
412
413	sizeToBytes(dsLengthBytes,dataLength);
414	for (i = 0; i < exe_params->SZ_SIZE_TYPE; i++)//ST: 4 or 8
415	(*newByteData)[k++] = dsLengthBytes[i];
416
417	if(sysEndianType==BIG_ENDIAN_SYSTEM)
418	memcpy((newByteData)+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE, oriData, dataLengthdoubleSize);
419	else
420	{
421	unsigned char* p = (*newByteData)+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE;
422	for(i=0;i<dataLength;i++,p+=doubleSize)
423	doubleToBytes(p, oriData[i]);
424	}
425	*outSize = totalByteLength;
426	}
427
428
429	char SZ_compress_args_double_NoCkRngeNoGzip_1D(unsigned char** newByteData, double *oriData,
430	size_t dataLength, double realPrecision, size_t *outSize, double valueRangeSize, double medianValue_d)
431	{
432	char compressionType = 0;
433	TightDataPointStorageD* tdps = NULL;
434	#ifdef HAVE_TIMECMPR
435	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
436	{
437	int timestep = sz_tsc->currentStep;
438	if(timestep % confparams_cpr->snapshotCmprStep != 0)
439	{
440	tdps = SZ_compress_double_1D_MDQ_ts(oriData, dataLength, multisteps, realPrecision, valueRangeSize, medianValue_d);
441	compressionType = 1; //time-series based compression
442	}
443	else
444	{
445	tdps = SZ_compress_double_1D_MDQ(oriData, dataLength, realPrecision, valueRangeSize, medianValue_d);
446	compressionType = 0; //snapshot-based compression
447	multisteps->lastSnapshotStep = timestep;
448	}
449	}
450	else
451	#endif
452	tdps = SZ_compress_double_1D_MDQ(oriData, dataLength, realPrecision, valueRangeSize, medianValue_d);
453
454	convertTDPStoFlatBytes_double(tdps, newByteData, outSize);
455
456	if(outSize>dataLengthsizeof(double))
457	SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
458
459	free_TightDataPointStorageD(tdps);
460	return compressionType;
461	}
462
463	TightDataPointStorageD* SZ_compress_double_2D_MDQ(double *oriData, size_t r1, size_t r2, double realPrecision, double valueRangeSize, double medianValue_d)
464	{
465	#ifdef HAVE_TIMECMPR
466	double* decData = NULL;
467	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
468	decData = (double*)(multisteps->hist_data);
469	#endif
470
471	unsigned int quantization_intervals;
472	if(exe_params->optQuantMode==1)
473	{
474	quantization_intervals = optimize_intervals_double_2D_opt(oriData, r1, r2, realPrecision);
475	updateQuantizationInfo(quantization_intervals);
476	}
477	else
478	quantization_intervals = exe_params->intvCapacity;
479	size_t i,j;
480	int reqLength;
481	double pred1D, pred2D;
482	double diff = 0.0;
483	double itvNum = 0;
484	double P0, P1;
485
486	size_t dataLength = r1*r2;
487
488	P0 = (double)malloc(r2sizeof(double));
489	memset(P0, 0, r2*sizeof(double));
490	P1 = (double)malloc(r2sizeof(double));
491	memset(P1, 0, r2*sizeof(double));
492
493	double medianValue = medianValue_d;
494	short radExpo = getExponent_double(valueRangeSize/2);
495	computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
496
497	int* type = (int) malloc(dataLengthsizeof(int));
498	//type[dataLength]=0;
499
500	double* spaceFillingValue = oriData; //
501
502	DynamicIntArray *exactLeadNumArray;
503	new_DIA(&exactLeadNumArray, DynArrayInitLen);
504
505	DynamicByteArray *exactMidByteArray;
506	new_DBA(&exactMidByteArray, DynArrayInitLen);
507
508	DynamicIntArray *resiBitArray;
509	new_DIA(&resiBitArray, DynArrayInitLen);
510
511	type[0] = 0;
512
513	unsigned char preDataBytes[8];
514	longToBytes_bigEndian(preDataBytes, 0);
515
516	int reqBytesLength = reqLength/8;
517	int resiBitsLength = reqLength%8;
518
519	DoubleValueCompressElement vce = (DoubleValueCompressElement)malloc(sizeof(DoubleValueCompressElement));
520	LossyCompressionElement lce = (LossyCompressionElement)malloc(sizeof(LossyCompressionElement));
521
522	/* Process Row-0 data 0*/
523	type[0] = 0;
524	compressSingleDoubleValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
525	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
526	memcpy(preDataBytes,vce->curBytes,8);
527	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
528	P1[0] = vce->data;
529	#ifdef HAVE_TIMECMPR
530	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
531	decData[0] = vce->data;
532	#endif
533
534	/* Process Row-0 data 1*/
535	pred1D = P1[0];
536	diff = spaceFillingValue[1] - pred1D;
537
538	itvNum = fabs(diff)/realPrecision + 1;
539
540	if (itvNum < exe_params->intvCapacity)
541	{
542	if (diff < 0) itvNum = -itvNum;
543	type[1] = (int) (itvNum/2) + exe_params->intvRadius;
544	P1[1] = pred1D + 2 * (type[1] - exe_params->intvRadius) * realPrecision;
545	}
546	else
547	{
548	type[1] = 0;
549	compressSingleDoubleValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
550	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
551	memcpy(preDataBytes,vce->curBytes,8);
552	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
553	P1[1] = vce->data;
554	}
555	#ifdef HAVE_TIMECMPR
556	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
557	decData[1] = P1[1];
558	#endif
559
560	/* Process Row-0 data 2 --> data r2-1 */
561	for (j = 2; j < r2; j++)
562	{
563	pred1D = 2*P1[j-1] - P1[j-2];
564	diff = spaceFillingValue[j] - pred1D;
565
566	itvNum = fabs(diff)/realPrecision + 1;
567
568	if (itvNum < exe_params->intvCapacity)
569	{
570	if (diff < 0) itvNum = -itvNum;
571	type[j] = (int) (itvNum/2) + exe_params->intvRadius;
572	P1[j] = pred1D + 2 * (type[j] - exe_params->intvRadius) * realPrecision;
573	}
574	else
575	{
576	type[j] = 0;
577	compressSingleDoubleValue(vce, spaceFillingValue[j], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
578	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
579	memcpy(preDataBytes,vce->curBytes,8);
580	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
581	P1[j] = vce->data;
582	}
583	#ifdef HAVE_TIMECMPR
584	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
585	decData[j] = P1[j];
586	#endif
587	}
588
589	/* Process Row-1 --> Row-r1-1 */
590	size_t index;
591	for (i = 1; i < r1; i++)
592	{
593	/* Process row-i data 0 */
594	index = i*r2;
595	pred1D = P1[0];
596	diff = spaceFillingValue[index] - pred1D;
597
598	itvNum = fabs(diff)/realPrecision + 1;
599
600	if (itvNum < exe_params->intvCapacity)
601	{
602	if (diff < 0) itvNum = -itvNum;
603	type[index] = (int) (itvNum/2) + exe_params->intvRadius;
604	P0[0] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
605	}
606	else
607	{
608	type[index] = 0;
609	compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
610	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
611	memcpy(preDataBytes,vce->curBytes,8);
612	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
613	P0[0] = vce->data;
614	}
615	#ifdef HAVE_TIMECMPR
616	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
617	decData[index] = P0[0];
618	#endif
619
620	/* Process row-i data 1 --> r2-1*/
621	for (j = 1; j < r2; j++)
622	{
623	index = i*r2+j;
624	pred2D = P0[j-1] + P1[j] - P1[j-1];
625
626	diff = spaceFillingValue[index] - pred2D;
627
628	itvNum = fabs(diff)/realPrecision + 1;
629
630	if (itvNum < exe_params->intvCapacity)
631	{
632	if (diff < 0) itvNum = -itvNum;
633	type[index] = (int) (itvNum/2) + exe_params->intvRadius;
634	P0[j] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
635	}
636	else
637	{
638	type[index] = 0;
639	compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
640	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
641	memcpy(preDataBytes,vce->curBytes,8);
642	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
643	P0[j] = vce->data;
644	}
645	#ifdef HAVE_TIMECMPR
646	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
647	decData[index] = P0[j];
648	#endif
649	}
650
651	double *Pt;
652	Pt = P1;
653	P1 = P0;
654	P0 = Pt;
655	}
656
657	if(r2!=1)
658	free(P0);
659	free(P1);
660	size_t exactDataNum = exactLeadNumArray->size;
661
662	TightDataPointStorageD* tdps;
663
664	new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
665	type, exactMidByteArray->array, exactMidByteArray->size,
666	exactLeadNumArray->array,
667	resiBitArray->array, resiBitArray->size,
668	resiBitsLength,
669	realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
670
671	/* int sum =0;
672	for(i=0;i<dataLength;i++)
673	if(type[i]==0) sum++;
674	printf("opt_quantizations=%d, exactDataNum=%d, sum=%d\n",quantization_intervals, exactDataNum, sum);
675
676	for(i=0;i<dataLength;i++)
677	printf("%d ", type[i]);
678	printf("\n");*/
679
680	// printf("exactDataNum=%d, expSegmentsInBytes_size=%d, exactMidByteArray->size=%d\n",
681	// exactDataNum, expSegmentsInBytes_size, exactMidByteArray->size);
682
683	// for(i = 3800;i<3844;i++)
684	// printf("exactLeadNumArray->array[%d]=%d\n",i,exactLeadNumArray->array[i]);
685
686	//free memory
687	free_DIA(exactLeadNumArray);
688	free_DIA(resiBitArray);
689	free(type);
690	free(vce);
691	free(lce);
692	free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
693
694	return tdps;
695	}
696
697	/**
698	*
699	* Note: @r1 is high dimension
700	* @r2 is low dimension
701	* */
702	char SZ_compress_args_double_NoCkRngeNoGzip_2D(unsigned char** newByteData, double oriData, size_t r1, size_t r2, double realPrecision, size_t outSize, double valueRangeSize, double medianValue_d)
703	{
704	size_t dataLength = r1*r2;
705	char compressionType = 0;
706	TightDataPointStorageD* tdps = NULL;
707	#ifdef HAVE_TIMECMPR
708	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
709	{
710	int timestep = sz_tsc->currentStep;
711	if(timestep % confparams_cpr->snapshotCmprStep != 0)
712	{
713	tdps = SZ_compress_double_1D_MDQ_ts(oriData, dataLength, multisteps, realPrecision, valueRangeSize, medianValue_d);
714	compressionType = 1; //time-series based compression
715	}
716	else
717	{
718	tdps = SZ_compress_double_2D_MDQ(oriData, r1, r2, realPrecision, valueRangeSize, medianValue_d);
719	compressionType = 0; //snapshot-based compression
720	multisteps->lastSnapshotStep = timestep;
721	}
722	}
723	else
724	#endif
725	tdps = SZ_compress_double_2D_MDQ(oriData, r1, r2, realPrecision, valueRangeSize, medianValue_d);
726
727	convertTDPStoFlatBytes_double(tdps, newByteData, outSize);
728
729	if(outSize>dataLengthsizeof(double))
730	SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
731
732	free_TightDataPointStorageD(tdps);
733	return compressionType;
734	}
735
736	TightDataPointStorageD* SZ_compress_double_3D_MDQ(double *oriData, size_t r1, size_t r2, size_t r3, double realPrecision, double valueRangeSize, double medianValue_d)
737	{
738	#ifdef HAVE_TIMECMPR
739	double* decData = NULL;
740	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
741	decData = (double*)(multisteps->hist_data);
742	#endif
743
744	unsigned int quantization_intervals;
745	if(exe_params->optQuantMode==1)
746	{
747	quantization_intervals = optimize_intervals_double_3D_opt(oriData, r1, r2, r3, realPrecision);
748	updateQuantizationInfo(quantization_intervals);
749	}
750	else
751	quantization_intervals = exe_params->intvCapacity;
752	size_t i,j,k;
753	int reqLength;
754	double pred1D, pred2D, pred3D;
755	double diff = 0.0;
756	double itvNum = 0;
757	double P0, P1;
758
759	size_t dataLength = r1r2r3;
760
761	size_t r23 = r2*r3;
762
763	P0 = (double)malloc(r23sizeof(double));
764	P1 = (double)malloc(r23sizeof(double));
765
766	double medianValue = medianValue_d;
767	short radExpo = getExponent_double(valueRangeSize/2);
768	computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
769
770	int* type = (int) malloc(dataLengthsizeof(int));
771	//type[dataLength]=0;
772
773	double* spaceFillingValue = oriData; //
774
775	DynamicIntArray *exactLeadNumArray;
776	new_DIA(&exactLeadNumArray, DynArrayInitLen);
777
778	DynamicByteArray *exactMidByteArray;
779	new_DBA(&exactMidByteArray, DynArrayInitLen);
780
781	DynamicIntArray *resiBitArray;
782	new_DIA(&resiBitArray, DynArrayInitLen);
783
784	type[0] = 0;
785
786	unsigned char preDataBytes[8];
787	longToBytes_bigEndian(preDataBytes, 0);
788
789	int reqBytesLength = reqLength/8;
790	int resiBitsLength = reqLength%8;
791
792	DoubleValueCompressElement vce = (DoubleValueCompressElement)malloc(sizeof(DoubleValueCompressElement));
793	LossyCompressionElement lce = (LossyCompressionElement)malloc(sizeof(LossyCompressionElement));
794
795
796	/////////////////////////// Process layer-0 ///////////////////////////
797	/* Process Row-0 data 0*/
798	type[0] = 0;
799	compressSingleDoubleValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
800	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
801	memcpy(preDataBytes,vce->curBytes,8);
802	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
803	P1[0] = vce->data;
804	#ifdef HAVE_TIMECMPR
805	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
806	decData[0] = P1[0];
807	#endif
808
809	/* Process Row-0 data 1*/
810	pred1D = P1[0];
811	diff = spaceFillingValue[1] - pred1D;
812
813	itvNum = fabs(diff)/realPrecision + 1;
814
815	if (itvNum < exe_params->intvCapacity)
816	{
817	if (diff < 0) itvNum = -itvNum;
818	type[1] = (int) (itvNum/2) + exe_params->intvRadius;
819	P1[1] = pred1D + 2 * (type[1] - exe_params->intvRadius) * realPrecision;
820	}
821	else
822	{
823	type[1] = 0;
824	compressSingleDoubleValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
825	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
826	memcpy(preDataBytes,vce->curBytes,8);
827	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
828	P1[1] = vce->data;
829	}
830	#ifdef HAVE_TIMECMPR
831	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
832	decData[1] = P1[1];
833	#endif
834
835	/* Process Row-0 data 2 --> data r3-1 */
836	for (j = 2; j < r3; j++)
837	{
838	pred1D = 2*P1[j-1] - P1[j-2];
839	diff = spaceFillingValue[j] - pred1D;
840
841	itvNum = fabs(diff)/realPrecision + 1;
842
843	if (itvNum < exe_params->intvCapacity)
844	{
845	if (diff < 0) itvNum = -itvNum;
846	type[j] = (int) (itvNum/2) + exe_params->intvRadius;
847	P1[j] = pred1D + 2 * (type[j] - exe_params->intvRadius) * realPrecision;
848	}
849	else
850	{
851	type[j] = 0;
852	compressSingleDoubleValue(vce, spaceFillingValue[j], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
853	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
854	memcpy(preDataBytes,vce->curBytes,8);
855	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
856	P1[j] = vce->data;
857	}
858	#ifdef HAVE_TIMECMPR
859	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
860	decData[j] = P1[j];
861	#endif
862	}
863
864	/* Process Row-1 --> Row-r2-1 */
865	size_t index;
866	for (i = 1; i < r2; i++)
867	{
868	/* Process row-i data 0 */
869	index = i*r3;
870	pred1D = P1[index-r3];
871	diff = spaceFillingValue[index] - pred1D;
872
873	itvNum = fabs(diff)/realPrecision + 1;
874
875	if (itvNum < exe_params->intvCapacity)
876	{
877	if (diff < 0) itvNum = -itvNum;
878	type[index] = (int) (itvNum/2) + exe_params->intvRadius;
879	P1[index] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
880	}
881	else
882	{
883	type[index] = 0;
884	compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
885	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
886	memcpy(preDataBytes,vce->curBytes,8);
887	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
888	P1[index] = vce->data;
889	}
890	#ifdef HAVE_TIMECMPR
891	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
892	decData[index] = P1[index];
893	#endif
894
895	/* Process row-i data 1 --> data r3-1*/
896	for (j = 1; j < r3; j++)
897	{
898	index = i*r3+j;
899	pred2D = P1[index-1] + P1[index-r3] - P1[index-r3-1];
900
901	diff = spaceFillingValue[index] - pred2D;
902
903	itvNum = fabs(diff)/realPrecision + 1;
904
905	if (itvNum < exe_params->intvCapacity)
906	{
907	if (diff < 0) itvNum = -itvNum;
908	type[index] = (int) (itvNum/2) + exe_params->intvRadius;
909	P1[index] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
910	}
911	else
912	{
913	type[index] = 0;
914	compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
915	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
916	memcpy(preDataBytes,vce->curBytes,8);
917	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
918	P1[index] = vce->data;
919	}
920	#ifdef HAVE_TIMECMPR
921	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
922	decData[index] = P1[index];
923	#endif
924	}
925	}
926
927
928	/////////////////////////// Process layer-1 --> layer-r1-1 ///////////////////////////
929
930	for (k = 1; k < r1; k++)
931	{
932	/* Process Row-0 data 0*/
933	index = k*r23;
934	pred1D = P1[0];
935	diff = spaceFillingValue[index] - pred1D;
936
937	itvNum = fabs(diff)/realPrecision + 1;
938
939	if (itvNum < exe_params->intvCapacity)
940	{
941	if (diff < 0) itvNum = -itvNum;
942	type[index] = (int) (itvNum/2) + exe_params->intvRadius;
943	P0[0] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
944	}
945	else
946	{
947	type[index] = 0;
948	compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
949	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
950	memcpy(preDataBytes,vce->curBytes,8);
951	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
952	P0[0] = vce->data;
953	}
954	#ifdef HAVE_TIMECMPR
955	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
956	decData[index] = P0[0];
957	#endif
958
959	/* Process Row-0 data 1 --> data r3-1 */
960	for (j = 1; j < r3; j++)
961	{
962	//index = kr2r3+j;
963	index ++;
964	pred2D = P0[j-1] + P1[j] - P1[j-1];
965	diff = spaceFillingValue[index] - pred2D;
966
967	itvNum = fabs(diff)/realPrecision + 1;
968
969	if (itvNum < exe_params->intvCapacity)
970	{
971	if (diff < 0) itvNum = -itvNum;
972	type[index] = (int) (itvNum/2) + exe_params->intvRadius;
973	P0[j] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
974	}
975	else
976	{
977	type[index] = 0;
978	compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
979	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
980	memcpy(preDataBytes,vce->curBytes,8);
981	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
982	P0[j] = vce->data;
983	}
984	#ifdef HAVE_TIMECMPR
985	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
986	decData[index] = P0[j];
987	#endif
988	}
989
990	/* Process Row-1 --> Row-r2-1 */
991	size_t index2D;
992	for (i = 1; i < r2; i++)
993	{
994	/* Process Row-i data 0 */
995	index = kr23 + ir3;
996	index2D = i*r3;
997	pred2D = P0[index2D-r3] + P1[index2D] - P1[index2D-r3];
998	diff = spaceFillingValue[index] - pred2D;
999
1000	itvNum = fabs(diff)/realPrecision + 1;
1001
1002	if (itvNum < exe_params->intvCapacity)
1003	{
1004	if (diff < 0) itvNum = -itvNum;
1005	type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1006	P0[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1007	}
1008	else
1009	{
1010	type[index] = 0;
1011	compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1012	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1013	memcpy(preDataBytes,vce->curBytes,8);
1014	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1015	P0[index2D] = vce->data;
1016	}
1017	#ifdef HAVE_TIMECMPR
1018	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1019	decData[index] = P0[index2D];
1020	#endif
1021
1022	/* Process Row-i data 1 --> data r3-1 */
1023	for (j = 1; j < r3; j++)
1024	{
1025	//index = kr2r3 + i*r3 + j;
1026	index ++;
1027	index2D = i*r3 + j;
1028	pred3D = P0[index2D-1] + P0[index2D-r3]+ P1[index2D] - P0[index2D-r3-1] - P1[index2D-r3] - P1[index2D-1] + P1[index2D-r3-1];
1029	diff = spaceFillingValue[index] - pred3D;
1030
1031	itvNum = fabs(diff)/realPrecision + 1;
1032
1033	if (itvNum < exe_params->intvCapacity)
1034	{
1035	if (diff < 0) itvNum = -itvNum;
1036	type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1037	P0[index2D] = pred3D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1038	}
1039	else
1040	{
1041	type[index] = 0;
1042	compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1043	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1044	memcpy(preDataBytes,vce->curBytes,8);
1045	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1046	P0[index2D] = vce->data;
1047	}
1048	#ifdef HAVE_TIMECMPR
1049	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1050	decData[index] = P0[index2D];
1051	#endif
1052	}
1053	}
1054
1055	double *Pt;
1056	Pt = P1;
1057	P1 = P0;
1058	P0 = Pt;
1059	}
1060	if(r23!=1)
1061	free(P0);
1062	free(P1);
1063	size_t exactDataNum = exactLeadNumArray->size;
1064
1065	TightDataPointStorageD* tdps;
1066
1067	new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
1068	type, exactMidByteArray->array, exactMidByteArray->size,
1069	exactLeadNumArray->array,
1070	resiBitArray->array, resiBitArray->size,
1071	resiBitsLength,
1072	realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
1073
1074	// printf("exactDataNum=%d, expSegmentsInBytes_size=%d, exactMidByteArray->size=%d\n",
1075	// exactDataNum, expSegmentsInBytes_size, exactMidByteArray->size);
1076
1077	// for(i = 3800;i<3844;i++)
1078	// printf("exactLeadNumArray->array[%d]=%d\n",i,exactLeadNumArray->array[i]);
1079
1080	//free memory
1081	free_DIA(exactLeadNumArray);
1082	free_DIA(resiBitArray);
1083	free(type);
1084	free(vce);
1085	free(lce);
1086	free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
1087
1088	return tdps;
1089	}
1090
1091
1092	char SZ_compress_args_double_NoCkRngeNoGzip_3D(unsigned char** newByteData, double oriData, size_t r1, size_t r2, size_t r3, double realPrecision, size_t outSize, double valueRangeSize, double medianValue_d)
1093	{
1094	size_t dataLength = r1r2r3;
1095	char compressionType = 0;
1096	TightDataPointStorageD* tdps = NULL;
1097	#ifdef HAVE_TIMECMPR
1098	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1099	{
1100	int timestep = sz_tsc->currentStep;
1101	if(timestep % confparams_cpr->snapshotCmprStep != 0)
1102	{
1103	tdps = SZ_compress_double_1D_MDQ_ts(oriData, dataLength, multisteps, realPrecision, valueRangeSize, medianValue_d);
1104	compressionType = 1; //time-series based compression
1105	}
1106	else
1107	{
1108	tdps = SZ_compress_double_3D_MDQ(oriData, r1, r2, r3, realPrecision, valueRangeSize, medianValue_d);
1109	compressionType = 0; //snapshot-based compression
1110	multisteps->lastSnapshotStep = timestep;
1111	}
1112	}
1113	else
1114	#endif
1115	tdps = SZ_compress_double_3D_MDQ(oriData, r1, r2, r3, realPrecision, valueRangeSize, medianValue_d);
1116
1117	convertTDPStoFlatBytes_double(tdps, newByteData, outSize);
1118
1119	if(outSize>dataLengthsizeof(double))
1120	SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
1121
1122	free_TightDataPointStorageD(tdps);
1123	return compressionType;
1124	}
1125
1126	TightDataPointStorageD* SZ_compress_double_4D_MDQ(double *oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision, double valueRangeSize, double medianValue_d)
1127	{
1128	unsigned int quantization_intervals;
1129	if(exe_params->optQuantMode==1)
1130	{
1131	quantization_intervals = optimize_intervals_double_4D(oriData, r1, r2, r3, r4, realPrecision);
1132	updateQuantizationInfo(quantization_intervals);
1133	}
1134	else
1135	quantization_intervals = exe_params->intvCapacity;
1136
1137	size_t i,j,k;
1138	int reqLength;
1139	double pred1D, pred2D, pred3D;
1140	double diff = 0.0;
1141	double itvNum = 0;
1142	double P0, P1;
1143
1144	size_t dataLength = r1r2r3*r4;
1145
1146	size_t r234 = r2r3r4;
1147	size_t r34 = r3*r4;
1148
1149	P0 = (double)malloc(r34sizeof(double));
1150	P1 = (double)malloc(r34sizeof(double));
1151
1152	double medianValue = medianValue_d;
1153	short radExpo = getExponent_double(valueRangeSize/2);
1154	computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
1155
1156	int* type = (int) malloc(dataLengthsizeof(int));
1157
1158	double* spaceFillingValue = oriData; //
1159
1160	DynamicIntArray *exactLeadNumArray;
1161	new_DIA(&exactLeadNumArray, DynArrayInitLen);
1162
1163	DynamicByteArray *exactMidByteArray;
1164	new_DBA(&exactMidByteArray, DynArrayInitLen);
1165
1166	DynamicIntArray *resiBitArray;
1167	new_DIA(&resiBitArray, DynArrayInitLen);
1168
1169	unsigned char preDataBytes[8];
1170	longToBytes_bigEndian(preDataBytes, 0);
1171
1172	int reqBytesLength = reqLength/8;
1173	int resiBitsLength = reqLength%8;
1174
1175	DoubleValueCompressElement vce = (DoubleValueCompressElement)malloc(sizeof(DoubleValueCompressElement));
1176	LossyCompressionElement lce = (LossyCompressionElement)malloc(sizeof(LossyCompressionElement));
1177
1178
1179	size_t l;
1180	for (l = 0; l < r1; l++)
1181	{
1182
1183	/////////////////////////// Process layer-0 ///////////////////////////
1184	/* Process Row-0 data 0*/
1185	size_t index = l*r234;
1186	size_t index2D = 0;
1187
1188	type[index] = 0;
1189	compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1190	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1191	memcpy(preDataBytes,vce->curBytes,8);
1192	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1193	P1[index2D] = vce->data;
1194
1195	/* Process Row-0 data 1*/
1196	index = l*r234+1;
1197	index2D = 1;
1198
1199	pred1D = P1[index2D-1];
1200	diff = spaceFillingValue[index] - pred1D;
1201
1202	itvNum = fabs(diff)/realPrecision + 1;
1203
1204	if (itvNum < exe_params->intvCapacity)
1205	{
1206	if (diff < 0) itvNum = -itvNum;
1207	type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1208	P1[index2D] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1209	}
1210	else
1211	{
1212	type[index] = 0;
1213	compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1214	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1215	memcpy(preDataBytes,vce->curBytes,8);
1216	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1217	P1[index2D] = vce->data;
1218	}
1219
1220	/* Process Row-0 data 2 --> data r4-1 */
1221	for (j = 2; j < r4; j++)
1222	{
1223	index = l*r234+j;
1224	index2D = j;
1225
1226	pred1D = 2*P1[index2D-1] - P1[index2D-2];
1227	diff = spaceFillingValue[index] - pred1D;
1228
1229	itvNum = fabs(diff)/realPrecision + 1;
1230
1231	if (itvNum < exe_params->intvCapacity)
1232	{
1233	if (diff < 0) itvNum = -itvNum;
1234	type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1235	P1[index2D] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1236	}
1237	else
1238	{
1239	type[index] = 0;
1240	compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1241	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1242	memcpy(preDataBytes,vce->curBytes,8);
1243	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1244	P1[index2D] = vce->data;
1245	}
1246	}
1247
1248	/* Process Row-1 --> Row-r3-1 */
1249	for (i = 1; i < r3; i++)
1250	{
1251	/* Process row-i data 0 */
1252	index = lr234+ir4;
1253	index2D = i*r4;
1254
1255	pred1D = P1[index2D-r4];
1256	diff = spaceFillingValue[index] - pred1D;
1257
1258	itvNum = fabs(diff)/realPrecision + 1;
1259
1260	if (itvNum < exe_params->intvCapacity)
1261	{
1262	if (diff < 0) itvNum = -itvNum;
1263	type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1264	P1[index2D] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1265	}
1266	else
1267	{
1268	type[index] = 0;
1269	compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1270	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1271	memcpy(preDataBytes,vce->curBytes,8);
1272	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1273	P1[index2D] = vce->data;
1274	}
1275
1276	/* Process row-i data 1 --> data r4-1*/
1277	for (j = 1; j < r4; j++)
1278	{
1279	index = lr234+ir4+j;
1280	index2D = i*r4+j;
1281
1282	pred2D = P1[index2D-1] + P1[index2D-r4] - P1[index2D-r4-1];
1283
1284	diff = spaceFillingValue[index] - pred2D;
1285
1286	itvNum = fabs(diff)/realPrecision + 1;
1287
1288	if (itvNum < exe_params->intvCapacity)
1289	{
1290	if (diff < 0) itvNum = -itvNum;
1291	type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1292	P1[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1293	}
1294	else
1295	{
1296	type[index] = 0;
1297	compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1298	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1299	memcpy(preDataBytes,vce->curBytes,8);
1300	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1301	P1[index2D] = vce->data;
1302	}
1303	}
1304	}
1305
1306
1307	/////////////////////////// Process layer-1 --> layer-r2-1 ///////////////////////////
1308
1309	for (k = 1; k < r2; k++)
1310	{
1311	/* Process Row-0 data 0*/
1312	index = lr234+kr34;
1313	index2D = 0;
1314
1315	pred1D = P1[index2D];
1316	diff = spaceFillingValue[index] - pred1D;
1317
1318	itvNum = fabs(diff)/realPrecision + 1;
1319
1320	if (itvNum < exe_params->intvCapacity)
1321	{
1322	if (diff < 0) itvNum = -itvNum;
1323	type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1324	P0[index2D] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1325	}
1326	else
1327	{
1328	type[index] = 0;
1329	compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1330	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1331	memcpy(preDataBytes,vce->curBytes,8);
1332	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1333	P0[index2D] = vce->data;
1334	}
1335
1336
1337	/* Process Row-0 data 1 --> data r4-1 */
1338	for (j = 1; j < r4; j++)
1339	{
1340	index = lr234+kr34+j;
1341	index2D = j;
1342
1343	pred2D = P0[index2D-1] + P1[index2D] - P1[index2D-1];
1344	diff = spaceFillingValue[index] - pred2D;
1345
1346	itvNum = fabs(diff)/realPrecision + 1;
1347
1348	if (itvNum < exe_params->intvCapacity)
1349	{
1350	if (diff < 0) itvNum = -itvNum;
1351	type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1352	P0[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1353	}
1354	else
1355	{
1356	type[index] = 0;
1357	compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1358	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1359	memcpy(preDataBytes,vce->curBytes,8);
1360	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1361	P0[index2D] = vce->data;
1362	}
1363	}
1364
1365	/* Process Row-1 --> Row-r3-1 */
1366	for (i = 1; i < r3; i++)
1367	{
1368	/* Process Row-i data 0 */
1369	index = lr234+kr34+i*r4;
1370	index2D = i*r4;
1371
1372	pred2D = P0[index2D-r4] + P1[index2D] - P1[index2D-r4];
1373	diff = spaceFillingValue[index] - pred2D;
1374
1375	itvNum = fabs(diff)/realPrecision + 1;
1376
1377	if (itvNum < exe_params->intvCapacity)
1378	{
1379	if (diff < 0) itvNum = -itvNum;
1380	type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1381	P0[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1382	}
1383	else
1384	{
1385	type[index] = 0;
1386	compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1387	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1388	memcpy(preDataBytes,vce->curBytes,8);
1389	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1390	P0[index2D] = vce->data;
1391	}
1392
1393	/* Process Row-i data 1 --> data r4-1 */
1394	for (j = 1; j < r4; j++)
1395	{
1396	index = lr234+kr34+i*r4+j;
1397	index2D = i*r4+j;
1398
1399	pred3D = P0[index2D-1] + P0[index2D-r4]+ P1[index2D] - P0[index2D-r4-1] - P1[index2D-r4] - P1[index2D-1] + P1[index2D-r4-1];
1400	diff = spaceFillingValue[index] - pred3D;
1401
1402
1403	itvNum = fabs(diff)/realPrecision + 1;
1404
1405	if (itvNum < exe_params->intvCapacity)
1406	{
1407	if (diff < 0) itvNum = -itvNum;
1408	type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1409	P0[index2D] = pred3D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1410	}
1411	else
1412	{
1413	type[index] = 0;
1414	compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1415	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1416	memcpy(preDataBytes,vce->curBytes,8);
1417	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1418	P0[index2D] = vce->data;
1419	}
1420	}
1421	}
1422
1423	double *Pt;
1424	Pt = P1;
1425	P1 = P0;
1426	P0 = Pt;
1427	}
1428	}
1429
1430	free(P0);
1431	free(P1);
1432	size_t exactDataNum = exactLeadNumArray->size;
1433
1434	TightDataPointStorageD* tdps;
1435
1436	new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
1437	type, exactMidByteArray->array, exactMidByteArray->size,
1438	exactLeadNumArray->array,
1439	resiBitArray->array, resiBitArray->size,
1440	resiBitsLength,
1441	realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
1442
1443	//free memory
1444	free_DIA(exactLeadNumArray);
1445	free_DIA(resiBitArray);
1446	free(type);
1447	free(vce);
1448	free(lce);
1449	free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
1450
1451	return tdps;
1452	}
1453
1454
1455	char SZ_compress_args_double_NoCkRngeNoGzip_4D(unsigned char** newByteData, double oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision, size_t outSize, double valueRangeSize, double medianValue_d)
1456	{
1457	TightDataPointStorageD* tdps = SZ_compress_double_4D_MDQ(oriData, r1, r2, r3, r4, realPrecision, valueRangeSize, medianValue_d);
1458
1459	convertTDPStoFlatBytes_double(tdps, newByteData, outSize);
1460
1461	size_t dataLength = r1r2r3*r4;
1462	if(outSize>dataLengthsizeof(double))
1463	SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
1464
1465	free_TightDataPointStorageD(tdps);
1466	return 0;
1467	}
1468
1469	void SZ_compress_args_double_withinRange(unsigned char** newByteData, double oriData, size_t dataLength, size_t outSize)
1470	{
1471	TightDataPointStorageD* tdps = (TightDataPointStorageD*) malloc(sizeof(TightDataPointStorageD));
1472	tdps->rtypeArray = NULL;
1473	tdps->typeArray = NULL;
1474	tdps->leadNumArray = NULL;
1475	tdps->residualMidBits = NULL;
1476
1477	tdps->allSameData = 1;
1478	tdps->dataSeriesLength = dataLength;
1479	tdps->exactMidBytes = (unsigned char)malloc(sizeof(unsigned char)8);
1480	tdps->pwrErrBoundBytes = NULL;
1481	tdps->isLossless = 0;
1482	double value = oriData[0];
1483	doubleToBytes(tdps->exactMidBytes, value);
1484	tdps->exactMidBytes_size = 8;
1485
1486	size_t tmpOutSize;
1487	//unsigned char *tmpByteData;
1488	convertTDPStoFlatBytes_double(tdps, newByteData, &tmpOutSize);
1489	//convertTDPStoFlatBytes_double(tdps, &tmpByteData, &tmpOutSize);
1490
1491	//newByteData = (unsigned char)malloc(sizeof(unsigned char)*16); //for floating-point data (1+3+4+4)
1492	//memcpy(*newByteData, tmpByteData, 16);
1493	*outSize = tmpOutSize;//12==3+1+8(double_size)+MetaDataByteLength
1494	free_TightDataPointStorageD(tdps);
1495	}
1496
1497	int SZ_compress_args_double_wRngeNoGzip(unsigned char** newByteData, double *oriData,
1498	size_t r5, size_t r4, size_t r3, size_t r2, size_t r1, size_t *outSize,
1499	int errBoundMode, double absErr_Bound, double relBoundRatio, double pwrErrRatio)
1500	{
1501	int status = SZ_SCES;
1502	size_t dataLength = computeDataLength(r5,r4,r3,r2,r1);
1503	double valueRangeSize = 0, medianValue = 0;
1504
1505	double min = computeRangeSize_double(oriData, dataLength, &valueRangeSize, &medianValue);
1506	double max = min+valueRangeSize;
1507	double realPrecision = getRealPrecision_double(valueRangeSize, errBoundMode, absErr_Bound, relBoundRatio, &status);
1508
1509	if(valueRangeSize <= realPrecision)
1510	{
1511	SZ_compress_args_double_withinRange(newByteData, oriData, dataLength, outSize);
1512	}
1513	else
1514	{
1515	if(r5==0&&r4==0&&r3==0&&r2==0)
1516	{
1517	if(errBoundMode>=PW_REL)
1518	{
1519	//SZ_compress_args_double_NoCkRngeNoGzip_1D_pwr(newByteData, oriData, realPrecision, r1, outSize, min, max);
1520	SZ_compress_args_double_NoCkRngeNoGzip_1D_pwrgroup(newByteData, oriData, r1, absErr_Bound, relBoundRatio, pwrErrRatio, valueRangeSize, medianValue, outSize);
1521	}
1522	else
1523	SZ_compress_args_double_NoCkRngeNoGzip_1D(newByteData, oriData, r1, realPrecision, outSize, valueRangeSize, medianValue);
1524	}
1525	else if(r5==0&&r4==0&&r3==0)
1526	{
1527	if(errBoundMode>=PW_REL)
1528	SZ_compress_args_double_NoCkRngeNoGzip_2D_pwr(newByteData, oriData, realPrecision, r2, r1, outSize, min, max);
1529	else
1530	SZ_compress_args_double_NoCkRngeNoGzip_2D(newByteData, oriData, r2, r1, realPrecision, outSize, valueRangeSize, medianValue);
1531	}
1532	else if(r5==0&&r4==0)
1533	{
1534	if(errBoundMode>=PW_REL)
1535	SZ_compress_args_double_NoCkRngeNoGzip_3D_pwr(newByteData, oriData, realPrecision, r3, r2, r1, outSize, min, max);
1536	else
1537	SZ_compress_args_double_NoCkRngeNoGzip_3D(newByteData, oriData, r3, r2, r1, realPrecision, outSize, valueRangeSize, medianValue);
1538	}
1539	else if(r5==0)
1540	{
1541	if(errBoundMode>=PW_REL)
1542	SZ_compress_args_double_NoCkRngeNoGzip_3D_pwr(newByteData, oriData, realPrecision, r4*r3, r2, r1, outSize, min, max);
1543	else
1544	SZ_compress_args_double_NoCkRngeNoGzip_3D(newByteData, oriData, r4*r3, r2, r1, realPrecision, outSize, valueRangeSize, medianValue);
1545	}
1546	}
1547	return status;
1548	}
1549
1550	int SZ_compress_args_double(unsigned char** newByteData, double *oriData,
1551	size_t r5, size_t r4, size_t r3, size_t r2, size_t r1, size_t *outSize,
1552	int errBoundMode, double absErr_Bound, double relBoundRatio, double pwRelBoundRatio)
1553	{
1554	confparams_cpr->errorBoundMode = errBoundMode;
1555	if(errBoundMode==PW_REL)
1556	{
1557	confparams_cpr->pw_relBoundRatio = pwRelBoundRatio;
1558	//confparams_cpr->pwr_type = SZ_PWR_MIN_TYPE;
1559	if(confparams_cpr->pwr_type==SZ_PWR_AVG_TYPE && r3 != 0 )
1560	{
1561	printf("Error: Current version doesn't support 3D data compression with point-wise relative error bound being based on pwrType=AVG\n");
1562	exit(0);
1563	return SZ_NSCS;
1564	}
1565	}
1566
1567	int status = SZ_SCES;
1568	size_t dataLength = computeDataLength(r5,r4,r3,r2,r1);
1569
1570	if(dataLength <= MIN_NUM_OF_ELEMENTS)
1571	{
1572	*newByteData = SZ_skip_compress_double(oriData, dataLength, outSize);
1573	return status;
1574	}
1575
1576	double valueRangeSize = 0, medianValue = 0;
1577
1578	double min = computeRangeSize_double(oriData, dataLength, &valueRangeSize, &medianValue);
1579	double max = min+valueRangeSize;
1580
1581	double realPrecision = 0;
1582
1583	if(confparams_cpr->errorBoundMode==PSNR)
1584	{
1585	confparams_cpr->errorBoundMode = ABS;
1586	realPrecision = confparams_cpr->absErrBound = computeABSErrBoundFromPSNR(confparams_cpr->psnr, (double)confparams_cpr->predThreshold, valueRangeSize);
1587	}
1588	else
1589	realPrecision = getRealPrecision_double(valueRangeSize, errBoundMode, absErr_Bound, relBoundRatio, &status);
1590
1591	if(valueRangeSize <= realPrecision)
1592	{
1593	SZ_compress_args_double_withinRange(newByteData, oriData, dataLength, outSize);
1594	}
1595	else
1596	{
1597	size_t tmpOutSize = 0;
1598	unsigned char* tmpByteData;
1599	if (r2==0)
1600	{
1601	if(confparams_cpr->errorBoundMode>=PW_REL)
1602	{
1603	//SZ_compress_args_double_NoCkRngeNoGzip_1D_pwr(&tmpByteData, oriData, realPrecision, r1, &tmpOutSize, min, max);
1604	SZ_compress_args_double_NoCkRngeNoGzip_1D_pwrgroup(&tmpByteData, oriData, r1, absErr_Bound, relBoundRatio, pwRelBoundRatio,
1605	valueRangeSize, medianValue, &tmpOutSize);
1606	}
1607	else
1608	#ifdef HAVE_TIMECMPR
1609	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1610	multisteps->compressionType = SZ_compress_args_double_NoCkRngeNoGzip_1D(&tmpByteData, oriData, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1611	else
1612	#endif
1613	SZ_compress_args_double_NoCkRngeNoGzip_1D(&tmpByteData, oriData, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1614	}
1615	else
1616	if (r3==0)
1617	{
1618	if(confparams_cpr->errorBoundMode>=PW_REL)
1619	SZ_compress_args_double_NoCkRngeNoGzip_2D_pwr(&tmpByteData, oriData, realPrecision, r2, r1, &tmpOutSize, min, max);
1620	else
1621	#ifdef HAVE_TIMECMPR
1622	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1623	multisteps->compressionType = SZ_compress_args_double_NoCkRngeNoGzip_2D(&tmpByteData, oriData, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1624	else
1625	#endif
1626	SZ_compress_args_double_NoCkRngeNoGzip_2D(&tmpByteData, oriData, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1627	}
1628	else
1629	if (r4==0)
1630	{
1631	if(confparams_cpr->errorBoundMode>=PW_REL)
1632	SZ_compress_args_double_NoCkRngeNoGzip_3D_pwr(&tmpByteData, oriData, realPrecision, r3, r2, r1, &tmpOutSize, min, max);
1633	else
1634	#ifdef HAVE_TIMECMPR
1635	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1636	multisteps->compressionType = SZ_compress_args_double_NoCkRngeNoGzip_3D(&tmpByteData, oriData, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1637	else
1638	#endif
1639	SZ_compress_args_double_NoCkRngeNoGzip_3D(&tmpByteData, oriData, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1640	}
1641	else
1642	if (r5==0)
1643	{
1644	if(confparams_cpr->errorBoundMode>=PW_REL)
1645	SZ_compress_args_double_NoCkRngeNoGzip_3D_pwr(&tmpByteData, oriData, realPrecision, r4*r3, r2, r1, &tmpOutSize, min, max);
1646	else
1647	#ifdef HAVE_TIMECMPR
1648	if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1649	multisteps->compressionType = SZ_compress_args_double_NoCkRngeNoGzip_4D(&tmpByteData, oriData, r4, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1650	else
1651	#endif
1652	SZ_compress_args_double_NoCkRngeNoGzip_4D(&tmpByteData, oriData, r4, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1653	}
1654	else
1655	{
1656	printf("Error: doesn't support 5 dimensions for now.\n");
1657	status = SZ_DERR;
1658	}
1659
1660	//Call Gzip to do the further compression.
1661	if(confparams_cpr->szMode==SZ_BEST_SPEED)
1662	{
1663	*outSize = tmpOutSize;
1664	*newByteData = tmpByteData;
1665	}
1666	else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION \|\| confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
1667	{
1668	*outSize = zlib_compress5(tmpByteData, tmpOutSize, newByteData, confparams_cpr->gzipMode);
1669	free(tmpByteData);
1670	}
1671	else
1672	{
1673	printf("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
1674	status = SZ_MERR;
1675	}
1676	}
1677
1678	return status;
1679	}
1680
1681	//TODO
1682	int SZ_compress_args_double_subblock(unsigned char* compressedBytes, double *oriData,
1683	size_t r5, size_t r4, size_t r3, size_t r2, size_t r1,
1684	size_t s5, size_t s4, size_t s3, size_t s2, size_t s1,
1685	size_t e5, size_t e4, size_t e3, size_t e2, size_t e1,
1686	size_t *outSize, int errBoundMode, double absErr_Bound, double relBoundRatio)
1687	{
1688	int status = SZ_SCES;
1689	double valueRangeSize = 0, medianValue = 0;
1690	computeRangeSize_double_subblock(oriData, &valueRangeSize, &medianValue, r5, r4, r3, r2, r1, s5, s4, s3, s2, s1, e5, e4, e3, e2, e1);
1691
1692	double realPrecision = getRealPrecision_double(valueRangeSize, errBoundMode, absErr_Bound, relBoundRatio, &status);
1693
1694	if(valueRangeSize <= realPrecision)
1695	{
1696	//TODO
1697	//SZ_compress_args_double_withinRange_subblock();
1698	}
1699	else
1700	{
1701	if (r2==0)
1702	{
1703	//TODO
1704	if(errBoundMode==PW_REL)
1705	{
1706	//TODO
1707	//SZ_compress_args_double_NoCkRngeNoGzip_1D_pwr_subblock();
1708	printf ("Current subblock version does not support point-wise relative error bound.\n");
1709	}
1710	else
1711	SZ_compress_args_double_NoCkRnge_1D_subblock(compressedBytes, oriData, realPrecision, outSize, valueRangeSize, medianValue, r1, s1, e1);
1712	}
1713	else
1714	if (r3==0)
1715	{
1716	if(errBoundMode==PW_REL)
1717	{
1718	//TODO
1719	//SZ_compress_args_double_NoCkRngeNoGzip_2D_pwr_subblock();
1720	printf ("Current subblock version does not support point-wise relative error bound.\n");
1721	}
1722	else
1723	SZ_compress_args_double_NoCkRnge_2D_subblock(compressedBytes, oriData, realPrecision, outSize, valueRangeSize, medianValue, r2, r1, s2, s1, e2, e1);
1724	}
1725	else
1726	if (r4==0)
1727	{
1728	if(errBoundMode==PW_REL)
1729	{
1730	//TODO
1731	//SZ_compress_args_double_NoCkRngeNoGzip_3D_pwr_subblock();
1732	printf ("Current subblock version does not support point-wise relative error bound.\n");
1733	}
1734	else
1735	SZ_compress_args_double_NoCkRnge_3D_subblock(compressedBytes, oriData, realPrecision, outSize, valueRangeSize, medianValue, r3, r2, r1, s3, s2, s1, e3, e2, e1);
1736	}
1737	else
1738	if (r5==0)
1739	{
1740	if(errBoundMode==PW_REL)
1741	{
1742	//TODO
1743	//SZ_compress_args_double_NoCkRngeNoGzip_4D_pwr_subblock();
1744	printf ("Current subblock version does not support point-wise relative error bound.\n");
1745	}
1746	else
1747	SZ_compress_args_double_NoCkRnge_4D_subblock(compressedBytes, oriData, realPrecision, outSize, valueRangeSize, medianValue, r4, r3, r2, r1, s4, s3, s2, s1, e4, e3, e2, e1);
1748	}
1749	else
1750	{
1751	printf("Error: doesn't support 5 dimensions for now.\n");
1752	status = SZ_DERR; //dimension error
1753	}
1754	}
1755	return status;
1756	}
1757
1758	void SZ_compress_args_double_NoCkRnge_1D_subblock(unsigned char* compressedBytes, double oriData, double realPrecision, size_t outSize, double valueRangeSize, double medianValue_d,
1759	size_t r1, size_t s1, size_t e1)
1760	{
1761	TightDataPointStorageD* tdps = SZ_compress_double_1D_MDQ_subblock(oriData, realPrecision, valueRangeSize, medianValue_d, r1, s1, e1);
1762
1763	if (confparams_cpr->szMode==SZ_BEST_SPEED)
1764	convertTDPStoFlatBytes_double_args(tdps, compressedBytes, outSize);
1765	else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION \|\| confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
1766	{
1767	unsigned char *tmpCompBytes;
1768	size_t tmpOutSize;
1769	convertTDPStoFlatBytes_double(tdps, &tmpCompBytes, &tmpOutSize);
1770	*outSize = zlib_compress3(tmpCompBytes, tmpOutSize, compressedBytes, confparams_cpr->gzipMode);
1771	free(tmpCompBytes);
1772	}
1773	else
1774	{
1775	printf ("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
1776	}
1777
1778	//TODO
1779	// if(outSize>dataLengthsizeof(double))
1780	// SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
1781
1782	free_TightDataPointStorageD(tdps);
1783	}
1784
1785	void SZ_compress_args_double_NoCkRnge_2D_subblock(unsigned char* compressedBytes, double oriData, double realPrecision, size_t outSize, double valueRangeSize, double medianValue_d,
1786	size_t r2, size_t r1, size_t s2, size_t s1, size_t e2, size_t e1)
1787	{
1788	TightDataPointStorageD* tdps = SZ_compress_double_2D_MDQ_subblock(oriData, realPrecision, valueRangeSize, medianValue_d, r2, r1, s2, s1, e2, e1);
1789
1790	if (confparams_cpr->szMode==SZ_BEST_SPEED)
1791	convertTDPStoFlatBytes_double_args(tdps, compressedBytes, outSize);
1792	else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION \|\| confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
1793	{
1794	unsigned char *tmpCompBytes;
1795	size_t tmpOutSize;
1796	convertTDPStoFlatBytes_double(tdps, &tmpCompBytes, &tmpOutSize);
1797	*outSize = zlib_compress3(tmpCompBytes, tmpOutSize, compressedBytes, confparams_cpr->gzipMode);
1798	free(tmpCompBytes);
1799	}
1800	else
1801	{
1802	printf ("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
1803	}
1804
1805	//TODO
1806	// if(outSize>dataLengthsizeof(double))
1807	// SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
1808
1809	free_TightDataPointStorageD(tdps);
1810	}
1811
1812	void SZ_compress_args_double_NoCkRnge_3D_subblock(unsigned char* compressedBytes, double oriData, double realPrecision, size_t outSize, double valueRangeSize, double medianValue_d,
1813	size_t r3, size_t r2, size_t r1, size_t s3, size_t s2, size_t s1, size_t e3, size_t e2, size_t e1)
1814	{
1815	TightDataPointStorageD* tdps = SZ_compress_double_3D_MDQ_subblock(oriData, realPrecision, valueRangeSize, medianValue_d, r3, r2, r1, s3, s2, s1, e3, e2, e1);
1816
1817	if (confparams_cpr->szMode==SZ_BEST_SPEED)
1818	convertTDPStoFlatBytes_double_args(tdps, compressedBytes, outSize);
1819	else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION \|\| confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
1820	{
1821	unsigned char *tmpCompBytes;
1822	size_t tmpOutSize;
1823	convertTDPStoFlatBytes_double(tdps, &tmpCompBytes, &tmpOutSize);
1824	*outSize = zlib_compress3(tmpCompBytes, tmpOutSize, compressedBytes, confparams_cpr->gzipMode);
1825	free(tmpCompBytes);
1826	}
1827	else
1828	{
1829	printf ("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
1830	}
1831
1832	//TODO
1833	// if(outSize>dataLengthsizeof(double))
1834	// SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
1835
1836	free_TightDataPointStorageD(tdps);
1837	}
1838
1839	void SZ_compress_args_double_NoCkRnge_4D_subblock(unsigned char* compressedBytes, double oriData, double realPrecision, size_t outSize, double valueRangeSize, double medianValue_d,
1840	size_t r4, size_t r3, size_t r2, size_t r1, size_t s4, size_t s3, size_t s2, size_t s1, size_t e4, size_t e3, size_t e2, size_t e1)
1841	{
1842	TightDataPointStorageD* tdps = SZ_compress_double_4D_MDQ_subblock(oriData, realPrecision, valueRangeSize, medianValue_d, r4, r3, r2, r1, s4, s3, s2, s1, e4, e3, e2, e1);
1843
1844	if (confparams_cpr->szMode==SZ_BEST_SPEED)
1845	convertTDPStoFlatBytes_double_args(tdps, compressedBytes, outSize);
1846	else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION \|\| confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
1847	{
1848	unsigned char *tmpCompBytes;
1849	size_t tmpOutSize;
1850	convertTDPStoFlatBytes_double(tdps, &tmpCompBytes, &tmpOutSize);
1851	*outSize = zlib_compress3(tmpCompBytes, tmpOutSize, compressedBytes, confparams_cpr->gzipMode);
1852	free(tmpCompBytes);
1853	}
1854	else
1855	{
1856	printf ("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
1857	}
1858
1859	//TODO
1860	// if(outSize>dataLengthsizeof(double))
1861	// SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
1862
1863	free_TightDataPointStorageD(tdps);
1864	}
1865
1866
1867	unsigned int optimize_intervals_double_1D_subblock(double *oriData, double realPrecision, size_t r1, size_t s1, size_t e1)
1868	{
1869	size_t dataLength = e1 - s1 + 1;
1870	oriData = oriData + s1;
1871
1872	size_t i = 0;
1873	unsigned long radiusIndex;
1874	double pred_value = 0, pred_err;
1875	int intervals = (int)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
1876	memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
1877	size_t totalSampleSize = dataLength/confparams_cpr->sampleDistance;
1878	for(i=2;i<dataLength;i++)
1879	{
1880	if(i%confparams_cpr->sampleDistance==0)
1881	{
1882	pred_value = 2*oriData[i-1] - oriData[i-2];
1883	//pred_value = oriData[i-1];
1884	pred_err = fabs(pred_value - oriData[i]);
1885	radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
1886	if(radiusIndex>=confparams_cpr->maxRangeRadius)
1887	radiusIndex = confparams_cpr->maxRangeRadius - 1;
1888	intervals[radiusIndex]++;
1889	}
1890	}
1891	//compute the appropriate number
1892	size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
1893	size_t sum = 0;
1894	for(i=0;i<confparams_cpr->maxRangeRadius;i++)
1895	{
1896	sum += intervals[i];
1897	if(sum>targetCount)
1898	break;
1899	}
1900
1901	if(i>=confparams_cpr->maxRangeRadius)
1902	i = confparams_cpr->maxRangeRadius-1;
1903	unsigned int accIntervals = 2*(i+1);
1904	unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
1905
1906	if(powerOf2<32)
1907	powerOf2 = 32;
1908
1909	free(intervals);
1910	return powerOf2;
1911	}
1912
1913	unsigned int optimize_intervals_double_2D_subblock(double *oriData, double realPrecision, size_t r1, size_t r2, size_t s1, size_t s2, size_t e1, size_t e2)
1914	{
1915	size_t R1 = e1 - s1 + 1;
1916	size_t R2 = e2 - s2 + 1;
1917
1918	size_t i,j, index;
1919	unsigned long radiusIndex;
1920	double pred_value = 0, pred_err;
1921	int intervals = (int)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
1922	memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
1923	size_t totalSampleSize = R1*R2/confparams_cpr->sampleDistance;
1924	for(i=s1+1;i<=e1;i++)
1925	{
1926	for(j=s2+1;j<=e2;j++)
1927	{
1928	if((i+j)%confparams_cpr->sampleDistance==0)
1929	{
1930	index = i*r2+j;
1931	pred_value = oriData[index-1] + oriData[index-r2] - oriData[index-r2-1];
1932	pred_err = fabs(pred_value - oriData[index]);
1933	radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
1934	if(radiusIndex>=confparams_cpr->maxRangeRadius)
1935	radiusIndex = confparams_cpr->maxRangeRadius - 1;
1936	intervals[radiusIndex]++;
1937	}
1938	}
1939	}
1940	//compute the appropriate number
1941	size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
1942	size_t sum = 0;
1943	for(i=0;i<confparams_cpr->maxRangeRadius;i++)
1944	{
1945	sum += intervals[i];
1946	if(sum>targetCount)
1947	break;
1948	}
1949	if(i>=confparams_cpr->maxRangeRadius)
1950	i = confparams_cpr->maxRangeRadius-1;
1951	unsigned int accIntervals = 2*(i+1);
1952	unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
1953
1954	if(powerOf2<32)
1955	powerOf2 = 32;
1956
1957	free(intervals);
1958	return powerOf2;
1959	}
1960
1961	unsigned int optimize_intervals_double_3D_subblock(double *oriData, double realPrecision, size_t r1, size_t r2, size_t r3, size_t s1, size_t s2, size_t s3, size_t e1, size_t e2, size_t e3)
1962	{
1963	size_t R1 = e1 - s1 + 1;
1964	size_t R2 = e2 - s2 + 1;
1965	size_t R3 = e3 - s3 + 1;
1966
1967	size_t r23 = r2*r3;
1968
1969	size_t i,j,k, index;
1970	unsigned long radiusIndex;
1971	double pred_value = 0, pred_err;
1972	int intervals = (int)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
1973	memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
1974	size_t totalSampleSize = R1R2R3/confparams_cpr->sampleDistance;
1975	for(i=s1+1;i<=e1;i++)
1976	{
1977	for(j=s2+1;j<=e2;j++)
1978	{
1979	for(k=s3+1;k<=e3;k++)
1980	{
1981	if((i+j+k)%confparams_cpr->sampleDistance==0)
1982	{
1983	index = ir23+jr3+k;
1984	pred_value = oriData[index-1] + oriData[index-r3] + oriData[index-r23]
1985	- oriData[index-1-r23] - oriData[index-r3-1] - oriData[index-r3-r23] + oriData[index-r3-r23-1];
1986	pred_err = fabs(pred_value - oriData[index]);
1987	radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
1988	if(radiusIndex>=confparams_cpr->maxRangeRadius)
1989	radiusIndex = confparams_cpr->maxRangeRadius - 1;
1990	intervals[radiusIndex]++;
1991	}
1992	}
1993
1994	}
1995	}
1996	//compute the appropriate number
1997	size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
1998	size_t sum = 0;
1999	for(i=0;i<confparams_cpr->maxRangeRadius;i++)
2000	{
2001	sum += intervals[i];
2002	if(sum>targetCount)
2003	break;
2004	}
2005	if(i>=confparams_cpr->maxRangeRadius)
2006	i = confparams_cpr->maxRangeRadius-1;
2007
2008	unsigned int accIntervals = 2*(i+1);
2009	unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
2010
2011	if(powerOf2<32)
2012	powerOf2 = 32;
2013
2014	free(intervals);
2015	return powerOf2;
2016	}
2017
2018	unsigned int optimize_intervals_double_4D_subblock(double *oriData, double realPrecision,
2019	size_t r1, size_t r2, size_t r3, size_t r4, size_t s1, size_t s2, size_t s3, size_t s4, size_t e1, size_t e2, size_t e3, size_t e4)
2020	{
2021	size_t R1 = e1 - s1 + 1;
2022	size_t R2 = e2 - s2 + 1;
2023	size_t R3 = e3 - s3 + 1;
2024	size_t R4 = e4 - s4 + 1;
2025
2026	size_t r34 = r3*r4;
2027	size_t r234 = r2r3r4;
2028
2029	size_t i,j,k,l, index;
2030	unsigned long radiusIndex;
2031	double pred_value = 0, pred_err;
2032	int intervals = (int)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
2033	memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
2034	size_t totalSampleSize = R1R2R3*R4/confparams_cpr->sampleDistance;
2035	for(i=s1+1;i<=e1;i++)
2036	{
2037	for(j=s2+1;j<=e2;j++)
2038	{
2039	for(k=s3+1;k<=e3;k++)
2040	{
2041	for(l=s4+1;l<=e4;l++)
2042	{
2043	if((i+j+k+l)%confparams_cpr->sampleDistance==0)
2044	{
2045	index = ir234+jr34+k*r4+l;
2046	pred_value = oriData[index-1] + oriData[index-r4] + oriData[index-r34]
2047	- oriData[index-1-r34] - oriData[index-r4-1] - oriData[index-r4-r34] + oriData[index-r4-r34-1];
2048	pred_err = fabs(pred_value - oriData[index]);
2049	radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
2050	if(radiusIndex>=confparams_cpr->maxRangeRadius)
2051	radiusIndex = confparams_cpr->maxRangeRadius - 1;
2052	intervals[radiusIndex]++;
2053	}
2054	}
2055	}
2056
2057	}
2058	}
2059	//compute the appropriate number
2060	size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
2061	size_t sum = 0;
2062	for(i=0;i<confparams_cpr->maxRangeRadius;i++)
2063	{
2064	sum += intervals[i];
2065	if(sum>targetCount)
2066	break;
2067	}
2068	if(i>=confparams_cpr->maxRangeRadius)
2069	i = confparams_cpr->maxRangeRadius-1;
2070
2071	unsigned int accIntervals = 2*(i+1);
2072	unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
2073
2074	if(powerOf2<32)
2075	powerOf2 = 32;
2076
2077	free(intervals);
2078	return powerOf2;
2079	}
2080
2081	TightDataPointStorageD* SZ_compress_double_1D_MDQ_subblock(double *oriData, double realPrecision, double valueRangeSize, double medianValue_d,
2082	size_t r1, size_t s1, size_t e1)
2083	{
2084	size_t dataLength = e1 - s1 + 1;
2085
2086	unsigned int quantization_intervals;
2087	if(exe_params->optQuantMode==1)
2088	quantization_intervals = optimize_intervals_double_1D_subblock(oriData, realPrecision, r1, s1, e1);
2089	else
2090	quantization_intervals = exe_params->intvCapacity;
2091	updateQuantizationInfo(quantization_intervals);
2092
2093	size_t i;
2094	int reqLength;
2095	double medianValue = medianValue_d;
2096	short radExpo = getExponent_double(valueRangeSize/2);
2097
2098	computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
2099
2100	int* type = (int) malloc(dataLengthsizeof(int));
2101
2102	double* spaceFillingValue = oriData + s1; //
2103
2104	DynamicIntArray *exactLeadNumArray;
2105	new_DIA(&exactLeadNumArray, DynArrayInitLen);
2106
2107	DynamicByteArray *exactMidByteArray;
2108	new_DBA(&exactMidByteArray, DynArrayInitLen);
2109
2110	DynamicIntArray *resiBitArray;
2111	new_DIA(&resiBitArray, DynArrayInitLen);
2112
2113	type[0] = 0;
2114
2115	unsigned char preDataBytes[8];
2116	longToBytes_bigEndian(preDataBytes, 0);
2117
2118	int reqBytesLength = reqLength/8;
2119	int resiBitsLength = reqLength%8;
2120	double last3CmprsData[3] = {0};
2121
2122	DoubleValueCompressElement vce = (DoubleValueCompressElement)malloc(sizeof(DoubleValueCompressElement));
2123	LossyCompressionElement lce = (LossyCompressionElement)malloc(sizeof(LossyCompressionElement));
2124
2125	//add the first data
2126	compressSingleDoubleValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2127	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2128	memcpy(preDataBytes,vce->curBytes,8);
2129	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2130	listAdd_double(last3CmprsData, vce->data);
2131
2132	//add the second data
2133	type[1] = 0;
2134	compressSingleDoubleValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2135	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2136	memcpy(preDataBytes,vce->curBytes,8);
2137	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2138	listAdd_double(last3CmprsData, vce->data);
2139
2140	int state;
2141	double checkRadius;
2142	double curData;
2143	double pred;
2144	double predAbsErr;
2145	checkRadius = (exe_params->intvCapacity-1)*realPrecision;
2146	double interval = 2*realPrecision;
2147
2148	for(i=2;i<dataLength;i++)
2149	{
2150	//printf("%.30G\n",last3CmprsData[0]);
2151	curData = spaceFillingValue[i];
2152	pred = 2*last3CmprsData[0] - last3CmprsData[1];
2153	//pred = last3CmprsData[0];
2154	predAbsErr = fabs(curData - pred);
2155	if(predAbsErr<=checkRadius)
2156	{
2157	state = (predAbsErr/realPrecision+1)/2;
2158	if(curData>=pred)
2159	{
2160	type[i] = exe_params->intvRadius+state;
2161	pred = pred + state*interval;
2162	}
2163	else //curData<pred
2164	{
2165	type[i] = exe_params->intvRadius-state;
2166	pred = pred - state*interval;
2167	}
2168	listAdd_double(last3CmprsData, pred);
2169	continue;
2170	}
2171
2172	//unpredictable data processing
2173	type[i] = 0;
2174	compressSingleDoubleValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2175	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2176	memcpy(preDataBytes,vce->curBytes,8);
2177	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2178
2179	listAdd_double(last3CmprsData, vce->data);
2180	}//end of for
2181
2182	size_t exactDataNum = exactLeadNumArray->size;
2183
2184	TightDataPointStorageD* tdps;
2185
2186	new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
2187	type, exactMidByteArray->array, exactMidByteArray->size,
2188	exactLeadNumArray->array,
2189	resiBitArray->array, resiBitArray->size,
2190	resiBitsLength,
2191	realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
2192
2193	//free memory
2194	free_DIA(exactLeadNumArray);
2195	free_DIA(resiBitArray);
2196	free(type);
2197	free(vce);
2198	free(lce);
2199	free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
2200
2201	return tdps;
2202	}
2203
2204
2205	TightDataPointStorageD* SZ_compress_double_2D_MDQ_subblock(double *oriData, double realPrecision, double valueRangeSize, double medianValue_d,
2206	size_t r1, size_t r2, size_t s1, size_t s2, size_t e1, size_t e2)
2207	{
2208	unsigned int quantization_intervals;
2209	if(exe_params->optQuantMode==1)
2210	{
2211	quantization_intervals = optimize_intervals_double_2D_subblock(oriData, realPrecision, r1, r2, s1, s2, e1, e2);
2212	updateQuantizationInfo(quantization_intervals);
2213	}
2214	else
2215	quantization_intervals = exe_params->intvCapacity;
2216
2217	size_t i,j;
2218	int reqLength;
2219	double pred1D, pred2D;
2220	double diff = 0.0;
2221	double itvNum = 0;
2222	double P0, P1;
2223
2224	size_t R1 = e1 - s1 + 1;
2225	size_t R2 = e2 - s2 + 1;
2226	size_t dataLength = R1*R2;
2227
2228	P0 = (double)malloc(R2sizeof(double));
2229	memset(P0, 0, R2*sizeof(double));
2230	P1 = (double)malloc(R2sizeof(double));
2231	memset(P1, 0, R2*sizeof(double));
2232
2233	double medianValue = medianValue_d;
2234	short radExpo = getExponent_double(valueRangeSize/2);
2235	computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
2236
2237	int* type = (int) malloc(dataLengthsizeof(int));
2238
2239	double* spaceFillingValue = oriData; //
2240
2241	DynamicIntArray *exactLeadNumArray;
2242	new_DIA(&exactLeadNumArray, DynArrayInitLen);
2243
2244	DynamicByteArray *exactMidByteArray;
2245	new_DBA(&exactMidByteArray, DynArrayInitLen);
2246
2247	DynamicIntArray *resiBitArray;
2248	new_DIA(&resiBitArray, DynArrayInitLen);
2249
2250	unsigned char preDataBytes[8];
2251	longToBytes_bigEndian(preDataBytes, 0);
2252
2253	int reqBytesLength = reqLength/8;
2254	int resiBitsLength = reqLength%8;
2255
2256	DoubleValueCompressElement vce = (DoubleValueCompressElement)malloc(sizeof(DoubleValueCompressElement));
2257	LossyCompressionElement lce = (LossyCompressionElement)malloc(sizeof(LossyCompressionElement));
2258
2259	/* Process Row-s1 data s2*/
2260	size_t gIndex;
2261	size_t lIndex;
2262
2263	gIndex = s1*r2+s2;
2264	lIndex = 0;
2265
2266	type[lIndex] = 0;
2267	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2268	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2269	memcpy(preDataBytes,vce->curBytes,8);
2270	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2271	P1[0] = vce->data;
2272
2273	/* Process Row-s1 data s2+1*/
2274	gIndex = s1*r2+(s2+1);
2275	lIndex = 1;
2276
2277	pred1D = P1[0];
2278	diff = spaceFillingValue[gIndex] - pred1D;
2279
2280	itvNum = fabs(diff)/realPrecision + 1;
2281
2282	if (itvNum < exe_params->intvCapacity)
2283	{
2284	if (diff < 0) itvNum = -itvNum;
2285	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2286	P1[1] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2287	}
2288	else
2289	{
2290	type[lIndex] = 0;
2291	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2292	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2293	memcpy(preDataBytes,vce->curBytes,8);
2294	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2295	P1[1] = vce->data;
2296	}
2297
2298	/* Process Row-s1 data s2+2 --> data e2 */
2299	for (j = 2; j < R2; j++)
2300	{
2301	gIndex = s1*r2+(s2+j);
2302	lIndex = j;
2303
2304	pred1D = 2*P1[j-1] - P1[j-2];
2305	diff = spaceFillingValue[gIndex] - pred1D;
2306
2307	itvNum = fabs(diff)/realPrecision + 1;
2308
2309	if (itvNum < exe_params->intvCapacity)
2310	{
2311	if (diff < 0) itvNum = -itvNum;
2312	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2313	P1[j] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2314	}
2315	else
2316	{
2317	type[lIndex] = 0;
2318	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2319	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2320	memcpy(preDataBytes,vce->curBytes,8);
2321	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2322	P1[j] = vce->data;
2323	}
2324	}
2325
2326	/* Process Row-s1+1 --> Row-e1 */
2327	for (i = 1; i < R1; i++)
2328	{
2329	/* Process row-s1+i data s2 */
2330	gIndex = (s1+i)*r2+s2;
2331	lIndex = i*R2;
2332
2333	pred1D = P1[0];
2334	diff = spaceFillingValue[gIndex] - pred1D;
2335
2336	itvNum = fabs(diff)/realPrecision + 1;
2337
2338	if (itvNum < exe_params->intvCapacity)
2339	{
2340	if (diff < 0) itvNum = -itvNum;
2341	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2342	P0[0] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2343	}
2344	else
2345	{
2346	type[lIndex] = 0;
2347	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2348	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2349	memcpy(preDataBytes,vce->curBytes,8);
2350	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2351	P0[0] = vce->data;
2352	}
2353
2354	/* Process row-s1+i data s2+1 --> e2 */
2355	for (j = 1; j < R2; j++)
2356	{
2357	gIndex = (s1+i)*r2+(s2+j);
2358	lIndex = i*R2+j;
2359
2360	pred2D = P0[j-1] + P1[j] - P1[j-1];
2361	diff = spaceFillingValue[gIndex] - pred2D;
2362
2363	itvNum = fabs(diff)/realPrecision + 1;
2364
2365	if (itvNum < exe_params->intvCapacity)
2366	{
2367	if (diff < 0) itvNum = -itvNum;
2368	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2369	P0[j] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2370	}
2371	else
2372	{
2373	type[lIndex] = 0;
2374	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2375	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2376	memcpy(preDataBytes,vce->curBytes,8);
2377	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2378	P0[j] = vce->data;
2379	}
2380	}
2381
2382	double *Pt;
2383	Pt = P1;
2384	P1 = P0;
2385	P0 = Pt;
2386	}
2387
2388	free(P0);
2389	free(P1);
2390	size_t exactDataNum = exactLeadNumArray->size;
2391
2392	TightDataPointStorageD* tdps;
2393
2394	new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
2395	type, exactMidByteArray->array, exactMidByteArray->size,
2396	exactLeadNumArray->array,
2397	resiBitArray->array, resiBitArray->size,
2398	resiBitsLength,
2399	realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
2400
2401	//free memory
2402	free_DIA(exactLeadNumArray);
2403	free_DIA(resiBitArray);
2404	free(type);
2405	free(vce);
2406	free(lce);
2407	free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
2408
2409	return tdps;
2410	}
2411
2412	TightDataPointStorageD* SZ_compress_double_3D_MDQ_subblock(double *oriData, double realPrecision, double valueRangeSize, double medianValue_d,
2413	size_t r1, size_t r2, size_t r3, size_t s1, size_t s2, size_t s3, size_t e1, size_t e2, size_t e3)
2414	{
2415	unsigned int quantization_intervals;
2416	if(exe_params->optQuantMode==1)
2417	{
2418	quantization_intervals = optimize_intervals_double_3D_subblock(oriData, realPrecision, r1, r2, r3, s1, s2, s3, e1, e2, e3);
2419	updateQuantizationInfo(quantization_intervals);
2420	}
2421	else
2422	quantization_intervals = exe_params->intvCapacity;
2423
2424	size_t i,j,k;
2425	int reqLength;
2426	double pred1D, pred2D, pred3D;
2427	double diff = 0.0;
2428	double itvNum = 0;
2429	double P0, P1;
2430
2431	size_t R1 = e1 - s1 + 1;
2432	size_t R2 = e2 - s2 + 1;
2433	size_t R3 = e3 - s3 + 1;
2434	size_t dataLength = R1R2R3;
2435
2436	size_t r23 = r2*r3;
2437	size_t R23 = R2*R3;
2438
2439	P0 = (double)malloc(R23sizeof(double));
2440	P1 = (double)malloc(R23sizeof(double));
2441
2442	double medianValue = medianValue_d;
2443	short radExpo = getExponent_double(valueRangeSize/2);
2444	computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
2445
2446	int* type = (int) malloc(dataLengthsizeof(int));
2447
2448	double* spaceFillingValue = oriData; //
2449
2450	DynamicIntArray *exactLeadNumArray;
2451	new_DIA(&exactLeadNumArray, DynArrayInitLen);
2452
2453	DynamicByteArray *exactMidByteArray;
2454	new_DBA(&exactMidByteArray, DynArrayInitLen);
2455
2456	DynamicIntArray *resiBitArray;
2457	new_DIA(&resiBitArray, DynArrayInitLen);
2458
2459	unsigned char preDataBytes[8];
2460	longToBytes_bigEndian(preDataBytes, 0);
2461
2462	int reqBytesLength = reqLength/8;
2463	int resiBitsLength = reqLength%8;
2464
2465	DoubleValueCompressElement vce = (DoubleValueCompressElement)malloc(sizeof(DoubleValueCompressElement));
2466	LossyCompressionElement lce = (LossyCompressionElement)malloc(sizeof(LossyCompressionElement));
2467
2468
2469	/////////////////////////// Process layer-s1 ///////////////////////////
2470	/* Process Row-s2 data s3*/
2471	size_t gIndex; //global index
2472	size_t lIndex; //local index
2473	size_t index2D; //local 2D index
2474
2475	gIndex = s1r23+s2r3+s3;
2476	lIndex = 0;
2477	index2D = 0;
2478
2479	type[lIndex] = 0;
2480	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2481	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2482	memcpy(preDataBytes,vce->curBytes,8);
2483	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2484	P1[index2D] = vce->data;
2485
2486	/* Process Row-s2 data s3+1*/
2487	gIndex = s1r23+s2r3+s3+1;
2488	lIndex = 1;
2489	index2D = 1;
2490
2491	pred1D = P1[index2D-1];
2492	diff = spaceFillingValue[gIndex] - pred1D;
2493
2494	itvNum = fabs(diff)/realPrecision + 1;
2495
2496	if (itvNum < exe_params->intvCapacity)
2497	{
2498	if (diff < 0) itvNum = -itvNum;
2499	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2500	P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2501	}
2502	else
2503	{
2504	type[lIndex] = 0;
2505	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2506	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2507	memcpy(preDataBytes,vce->curBytes,8);
2508	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2509	P1[index2D] = vce->data;
2510	}
2511
2512	/* Process Row-s2 data s3+2 --> data e3 */
2513	for (j = 2; j < R3; j++)
2514	{
2515	gIndex = s1r23+s2r3+s3+j;
2516	lIndex = j;
2517	index2D = j;
2518
2519	pred1D = 2*P1[index2D-1] - P1[index2D-2];
2520	diff = spaceFillingValue[gIndex] - pred1D;
2521
2522	itvNum = fabs(diff)/realPrecision + 1;
2523
2524	if (itvNum < exe_params->intvCapacity)
2525	{
2526	if (diff < 0) itvNum = -itvNum;
2527	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2528	P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2529	}
2530	else
2531	{
2532	type[lIndex] = 0;
2533	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2534	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2535	memcpy(preDataBytes,vce->curBytes,8);
2536	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2537	P1[index2D] = vce->data;
2538	}
2539	}
2540
2541	/* Process Row-s2+1 --> Row-e2 */
2542	for (i = 1; i < R2; i++)
2543	{
2544	/* Process row-s2+i data s3 */
2545	gIndex = s1r23+(s2+i)r3+s3;
2546	lIndex = i*R3;
2547	index2D = i*R3;
2548
2549	pred1D = P1[index2D-R3];
2550	diff = spaceFillingValue[gIndex] - pred1D;
2551
2552	itvNum = fabs(diff)/realPrecision + 1;
2553
2554	if (itvNum < exe_params->intvCapacity)
2555	{
2556	if (diff < 0) itvNum = -itvNum;
2557	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2558	P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2559	}
2560	else
2561	{
2562	type[lIndex] = 0;
2563	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2564	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2565	memcpy(preDataBytes,vce->curBytes,8);
2566	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2567	P1[index2D] = vce->data;
2568	}
2569
2570	/* Process row-s2+i data s3+1 --> data e3*/
2571	for (j = 1; j < R3; j++)
2572	{
2573	gIndex = s1r23+(s2+i)r3+s3+j;
2574	lIndex = i*R3+j;
2575	index2D = i*R3+j;
2576
2577	pred2D = P1[index2D-1] + P1[index2D-R3] - P1[index2D-R3-1];
2578	diff = spaceFillingValue[gIndex] - pred2D;
2579
2580	itvNum = fabs(diff)/realPrecision + 1;
2581
2582	if (itvNum < exe_params->intvCapacity)
2583	{
2584	if (diff < 0) itvNum = -itvNum;
2585	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2586	P1[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2587	}
2588	else
2589	{
2590	type[lIndex] = 0;
2591	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2592	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2593	memcpy(preDataBytes,vce->curBytes,8);
2594	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2595	P1[index2D] = vce->data;
2596	}
2597	}
2598	}
2599
2600
2601	/////////////////////////// Process layer-s1+1 --> layer-e1 ///////////////////////////
2602
2603	for (k = 1; k < R1; k++)
2604	{
2605	/* Process Row-s2 data s3*/
2606	gIndex = (s1+k)r23+s2r3+s3;
2607	lIndex = k*R23;
2608	index2D = 0;
2609
2610	pred1D = P1[index2D];
2611	diff = spaceFillingValue[gIndex] - pred1D;
2612
2613	itvNum = fabs(diff)/realPrecision + 1;
2614
2615	if (itvNum < exe_params->intvCapacity)
2616	{
2617	if (diff < 0) itvNum = -itvNum;
2618	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2619	P0[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2620	}
2621	else
2622	{
2623	type[lIndex] = 0;
2624	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2625	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2626	memcpy(preDataBytes,vce->curBytes,8);
2627	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2628	P0[index2D] = vce->data;
2629	}
2630
2631
2632	/* Process Row-s2 data s3+1 --> data e3 */
2633	for (j = 1; j < R3; j++)
2634	{
2635	gIndex = (s1+k)r23+s2r3+s3+j;
2636	lIndex = k*R23+j;
2637	index2D = j;
2638
2639	pred2D = P0[index2D-1] + P1[index2D] - P1[index2D-1];
2640	diff = spaceFillingValue[gIndex] - pred2D;
2641
2642	itvNum = fabs(diff)/realPrecision + 1;
2643
2644	if (itvNum < exe_params->intvCapacity)
2645	{
2646	if (diff < 0) itvNum = -itvNum;
2647	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2648	P0[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2649	}
2650	else
2651	{
2652	type[lIndex] = 0;
2653	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2654	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2655	memcpy(preDataBytes,vce->curBytes,8);
2656	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2657	P0[index2D] = vce->data;
2658	}
2659	}
2660
2661	/* Process Row-s2+1 --> Row-e2 */
2662	for (i = 1; i < R2; i++)
2663	{
2664	/* Process Row-s2+i data s3 */
2665	gIndex = (s1+k)r23+(s2+i)r3+s3;
2666	lIndex = kR23+iR3;
2667	index2D = i*R3;
2668
2669	pred2D = P0[index2D-R3] + P1[index2D] - P1[index2D-R3];
2670	diff = spaceFillingValue[gIndex] - pred2D;
2671
2672	itvNum = fabs(diff)/realPrecision + 1;
2673
2674	if (itvNum < exe_params->intvCapacity)
2675	{
2676	if (diff < 0) itvNum = -itvNum;
2677	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2678	P0[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2679	}
2680	else
2681	{
2682	type[lIndex] = 0;
2683	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2684	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2685	memcpy(preDataBytes,vce->curBytes,8);
2686	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2687	P0[index2D] = vce->data;
2688	}
2689
2690	/* Process Row-s2+i data s3+1 --> data e3 */
2691	for (j = 1; j < R3; j++)
2692	{
2693	gIndex = (s1+k)r23+(s2+i)r3+s3+j;
2694	lIndex = kR23+iR3+j;
2695	index2D = i*R3+j;
2696
2697	pred3D = P0[index2D-1] + P0[index2D-R3]+ P1[index2D] - P0[index2D-R3-1] - P1[index2D-R3] - P1[index2D-1] + P1[index2D-R3-1];
2698	diff = spaceFillingValue[gIndex] - pred3D;
2699
2700	itvNum = fabs(diff)/realPrecision + 1;
2701
2702	if (itvNum < exe_params->intvCapacity)
2703	{
2704	if (diff < 0) itvNum = -itvNum;
2705	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2706	P0[index2D] = pred3D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2707	}
2708	else
2709	{
2710	type[lIndex] = 0;
2711	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2712	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2713	memcpy(preDataBytes,vce->curBytes,8);
2714	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2715	P0[index2D] = vce->data;
2716	}
2717	}
2718	}
2719
2720	double *Pt;
2721	Pt = P1;
2722	P1 = P0;
2723	P0 = Pt;
2724	}
2725
2726	free(P0);
2727	free(P1);
2728	size_t exactDataNum = exactLeadNumArray->size;
2729
2730	TightDataPointStorageD* tdps;
2731
2732	new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
2733	type, exactMidByteArray->array, exactMidByteArray->size,
2734	exactLeadNumArray->array,
2735	resiBitArray->array, resiBitArray->size,
2736	resiBitsLength,
2737	realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
2738
2739	//free memory
2740	free_DIA(exactLeadNumArray);
2741	free_DIA(resiBitArray);
2742	free(type);
2743	free(vce);
2744	free(lce);
2745	free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
2746
2747	return tdps;
2748	}
2749
2750	TightDataPointStorageD* SZ_compress_double_4D_MDQ_subblock(double *oriData, double realPrecision, double valueRangeSize, double medianValue_d,
2751	size_t r1, size_t r2, size_t r3, size_t r4, size_t s1, size_t s2, size_t s3, size_t s4, size_t e1, size_t e2, size_t e3, size_t e4)
2752	{
2753	unsigned int quantization_intervals;
2754	if(exe_params->optQuantMode==1)
2755	{
2756	quantization_intervals = optimize_intervals_double_4D_subblock(oriData, realPrecision, r1, r2, r3, r4, s1, s2, s3, s4, e1, e2, e3, e4);
2757	updateQuantizationInfo(quantization_intervals);
2758	}
2759	else
2760	quantization_intervals = exe_params->intvCapacity;
2761
2762	size_t i,j,k;
2763	int reqLength;
2764	double pred1D, pred2D, pred3D;
2765	double diff = 0.0;
2766	double itvNum = 0;
2767	double P0, P1;
2768
2769	size_t R1 = e1 - s1 + 1;
2770	size_t R2 = e2 - s2 + 1;
2771	size_t R3 = e3 - s3 + 1;
2772	size_t R4 = e4 - s4 + 1;
2773
2774	size_t dataLength = R1R2R3*R4;
2775
2776	size_t r34 = r3*r4;
2777	size_t r234 = r2r3r4;
2778	size_t R34 = R3*R4;
2779	size_t R234 = R2R3R4;
2780
2781	P0 = (double)malloc(R34sizeof(double));
2782	P1 = (double)malloc(R34sizeof(double));
2783
2784	double medianValue = medianValue_d;
2785	short radExpo = getExponent_double(valueRangeSize/2);
2786	computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
2787
2788	int* type = (int) malloc(dataLengthsizeof(int));
2789
2790	double* spaceFillingValue = oriData; //
2791
2792	DynamicIntArray *exactLeadNumArray;
2793	new_DIA(&exactLeadNumArray, DynArrayInitLen);
2794
2795	DynamicByteArray *exactMidByteArray;
2796	new_DBA(&exactMidByteArray, DynArrayInitLen);
2797
2798	DynamicIntArray *resiBitArray;
2799	new_DIA(&resiBitArray, DynArrayInitLen);
2800
2801	unsigned char preDataBytes[8];
2802	longToBytes_bigEndian(preDataBytes, 0);
2803
2804	int reqBytesLength = reqLength/8;
2805	int resiBitsLength = reqLength%8;
2806
2807	DoubleValueCompressElement vce = (DoubleValueCompressElement)malloc(sizeof(DoubleValueCompressElement));
2808	LossyCompressionElement lce = (LossyCompressionElement)malloc(sizeof(LossyCompressionElement));
2809
2810	size_t l;
2811	for (l = 0; l < R1; l++)
2812	{
2813
2814	/////////////////////////// Process layer-s2 ///////////////////////////
2815	/* Process Row-s3 data s4*/
2816	size_t gIndex; //global index
2817	size_t lIndex; //local index
2818	size_t index2D; //local 2D index
2819
2820	gIndex = (s1+l)r234+s2r34+s3*r4+s4;
2821	lIndex = l*R234;
2822	index2D = 0;
2823
2824	type[lIndex] = 0;
2825	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2826	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2827	memcpy(preDataBytes,vce->curBytes,8);
2828	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2829	P1[index2D] = vce->data;
2830
2831	/* Process Row-s3 data s4+1*/
2832	gIndex = (s1+l)r234+s2r34+s3*r4+s4+1;
2833	lIndex = l*R234+1;
2834	index2D = 1;
2835
2836	pred1D = P1[index2D-1];
2837	diff = spaceFillingValue[gIndex] - pred1D;
2838
2839	itvNum = fabs(diff)/realPrecision + 1;
2840
2841	if (itvNum < exe_params->intvCapacity)
2842	{
2843	if (diff < 0) itvNum = -itvNum;
2844	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2845	P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2846	}
2847	else
2848	{
2849	type[lIndex] = 0;
2850	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2851	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2852	memcpy(preDataBytes,vce->curBytes,8);
2853	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2854	P1[index2D] = vce->data;
2855	}
2856
2857	/* Process Row-s3 data s4+2 --> data e4 */
2858	for (j = 2; j < R4; j++)
2859	{
2860	gIndex = (s1+l)r234+s2r34+s3*r4+s4+j;
2861	lIndex = l*R234+j;
2862	index2D = j;
2863
2864	pred1D = 2*P1[index2D-1] - P1[index2D-2];
2865	diff = spaceFillingValue[gIndex] - pred1D;
2866
2867	itvNum = fabs(diff)/realPrecision + 1;
2868
2869	if (itvNum < exe_params->intvCapacity)
2870	{
2871	if (diff < 0) itvNum = -itvNum;
2872	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2873	P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2874	}
2875	else
2876	{
2877	type[lIndex] = 0;
2878	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2879	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2880	memcpy(preDataBytes,vce->curBytes,8);
2881	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2882	P1[index2D] = vce->data;
2883	}
2884	}
2885
2886	/* Process Row-s3+1 --> Row-e3 */
2887	for (i = 1; i < R3; i++)
2888	{
2889	/* Process row-s2+i data s3 */
2890	gIndex = (s1+l)r234+s2r34+(s3+i)*r4+s4;
2891	lIndex = lR234+iR4;
2892	index2D = i*R4;
2893
2894	pred1D = P1[index2D-R4];
2895	diff = spaceFillingValue[gIndex] - pred1D;
2896
2897	itvNum = fabs(diff)/realPrecision + 1;
2898
2899	if (itvNum < exe_params->intvCapacity)
2900	{
2901	if (diff < 0) itvNum = -itvNum;
2902	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2903	P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2904	}
2905	else
2906	{
2907	type[lIndex] = 0;
2908	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2909	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2910	memcpy(preDataBytes,vce->curBytes,8);
2911	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2912	P1[index2D] = vce->data;
2913	}
2914
2915	/* Process row-s3+i data s4+1 --> data e4*/
2916	for (j = 1; j < R4; j++)
2917	{
2918	gIndex = (s1+l)r234+s2r34+(s3+i)*r4+s4+j;
2919	lIndex = lR234+iR4+j;
2920	index2D = i*R4+j;
2921
2922	pred2D = P1[index2D-1] + P1[index2D-R4] - P1[index2D-R4-1];
2923	diff = spaceFillingValue[gIndex] - pred2D;
2924
2925	itvNum = fabs(diff)/realPrecision + 1;
2926
2927	if (itvNum < exe_params->intvCapacity)
2928	{
2929	if (diff < 0) itvNum = -itvNum;
2930	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2931	P1[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2932	}
2933	else
2934	{
2935	type[lIndex] = 0;
2936	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2937	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2938	memcpy(preDataBytes,vce->curBytes,8);
2939	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2940	P1[index2D] = vce->data;
2941	}
2942	}
2943	}
2944
2945
2946	/////////////////////////// Process layer-s2+1 --> layer-e2 ///////////////////////////
2947
2948	for (k = 1; k < R2; k++)
2949	{
2950	/* Process Row-s3 data s4*/
2951	gIndex = (s1+l)r234+(s2+k)r34+s3*r4+s4;
2952	lIndex = lR234+kR34;
2953	index2D = 0;
2954
2955	pred1D = P1[index2D];
2956	diff = spaceFillingValue[gIndex] - pred1D;
2957
2958	itvNum = fabs(diff)/realPrecision + 1;
2959
2960	if (itvNum < exe_params->intvCapacity)
2961	{
2962	if (diff < 0) itvNum = -itvNum;
2963	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2964	P0[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2965	}
2966	else
2967	{
2968	type[lIndex] = 0;
2969	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2970	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2971	memcpy(preDataBytes,vce->curBytes,8);
2972	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2973	P0[index2D] = vce->data;
2974	}
2975
2976
2977	/* Process Row-s3 data s4+1 --> data e4 */
2978	for (j = 1; j < R4; j++)
2979	{
2980	gIndex = (s1+l)r234+(s2+k)r34+s3*r4+s4+j;
2981	lIndex = lR234+kR34+j;
2982	index2D = j;
2983
2984	pred2D = P0[index2D-1] + P1[index2D] - P1[index2D-1];
2985	diff = spaceFillingValue[gIndex] - pred2D;
2986
2987	itvNum = fabs(diff)/realPrecision + 1;
2988
2989	if (itvNum < exe_params->intvCapacity)
2990	{
2991	if (diff < 0) itvNum = -itvNum;
2992	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2993	P0[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2994	}
2995	else
2996	{
2997	type[lIndex] = 0;
2998	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2999	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
3000	memcpy(preDataBytes,vce->curBytes,8);
3001	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
3002	P0[index2D] = vce->data;
3003	}
3004	}
3005
3006	/* Process Row-s3+1 --> Row-e3 */
3007	for (i = 1; i < R3; i++)
3008	{
3009	/* Process Row-s3+i data s4 */
3010	gIndex = (s1+l)r234+(s2+k)r34+(s3+i)*r4+s4;
3011	lIndex = lR234+kR34+i*R4;
3012	index2D = i*R4;
3013
3014	pred2D = P0[index2D-R4] + P1[index2D] - P1[index2D-R4];
3015	diff = spaceFillingValue[gIndex] - pred2D;
3016
3017	itvNum = fabs(diff)/realPrecision + 1;
3018
3019	if (itvNum < exe_params->intvCapacity)
3020	{
3021	if (diff < 0) itvNum = -itvNum;
3022	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
3023	P0[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
3024	}
3025	else
3026	{
3027	type[lIndex] = 0;
3028	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
3029	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
3030	memcpy(preDataBytes,vce->curBytes,8);
3031	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
3032	P0[index2D] = vce->data;
3033	}
3034
3035	/* Process Row-s3+i data s4+1 --> data e4 */
3036	for (j = 1; j < R4; j++)
3037	{
3038	gIndex = (s1+l)r234+(s2+k)r34+(s3+i)*r4+s4+j;
3039	lIndex = lR234+kR34+i*R4+j;
3040	index2D = i*R4+j;
3041
3042	// printf ("global index = %d, local index = %d\n", gIndex, lIndex);
3043
3044	pred3D = P0[index2D-1] + P0[index2D-R4]+ P1[index2D] - P0[index2D-R4-1] - P1[index2D-R4] - P1[index2D-1] + P1[index2D-R4-1];
3045	diff = spaceFillingValue[gIndex] - pred3D;
3046
3047	itvNum = fabs(diff)/realPrecision + 1;
3048
3049	if (itvNum < exe_params->intvCapacity)
3050	{
3051	if (diff < 0) itvNum = -itvNum;
3052	type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
3053	P0[index2D] = pred3D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
3054	}
3055	else
3056	{
3057	type[lIndex] = 0;
3058	compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
3059	updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
3060	memcpy(preDataBytes,vce->curBytes,8);
3061	addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
3062	P0[index2D] = vce->data;
3063	}
3064	}
3065	}
3066
3067	double *Pt;
3068	Pt = P1;
3069	P1 = P0;
3070	P0 = Pt;
3071	}
3072	}
3073
3074	free(P0);
3075	free(P1);
3076	size_t exactDataNum = exactLeadNumArray->size;
3077
3078	TightDataPointStorageD* tdps;
3079
3080	new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
3081	type, exactMidByteArray->array, exactMidByteArray->size,
3082	exactLeadNumArray->array,
3083	resiBitArray->array, resiBitArray->size,
3084	resiBitsLength,
3085	realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
3086
3087	//free memory
3088	free_DIA(exactLeadNumArray);
3089	free_DIA(resiBitArray);
3090	free(type);
3091	free(vce);
3092	free(lce);
3093	free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
3094
3095	return tdps;
3096	}
3097
3098	/**
3099	*
3100	* This is a fast implementation for optimize_intervals_double_3D()
3101	* */
3102	unsigned int optimize_intervals_double_3D_opt(double *oriData, size_t r1, size_t r2, size_t r3, double realPrecision){
3103	size_t i;
3104	size_t radiusIndex;
3105	size_t r23=r2*r3;
3106	double pred_value = 0, pred_err;
3107	size_t intervals = (size_t)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
3108	memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
3109	size_t totalSampleSize = 0;
3110
3111	size_t offset_count = confparams_cpr->sampleDistance - 2; // count r3 offset
3112	size_t offset_count_2;
3113	double * data_pos = oriData + r23 + r3 + offset_count;
3114	size_t n1_count = 1, n2_count = 1; // count i,j sum
3115	size_t len = r1 * r2 * r3;
3116	while(data_pos - oriData < len){
3117	totalSampleSize++;
3118	pred_value = data_pos[-1] + data_pos[-r3] + data_pos[-r23] - data_pos[-1-r23] - data_pos[-r3-1] - data_pos[-r3-r23] + data_pos[-r3-r23-1];
3119	pred_err = fabs(pred_value - *data_pos);
3120	radiusIndex = (pred_err/realPrecision+1)/2;
3121	if(radiusIndex>=confparams_cpr->maxRangeRadius)
3122	{
3123	radiusIndex = confparams_cpr->maxRangeRadius - 1;
3124	//printf("radiusIndex=%d\n", radiusIndex);
3125	}
3126	intervals[radiusIndex]++;
3127	// printf("TEST: %ld, i: %ld\tj: %ld\tk: %ld\n", data_pos - oriData);
3128	// fflush(stdout);
3129	offset_count += confparams_cpr->sampleDistance;
3130	if(offset_count >= r3){
3131	n2_count ++;
3132	if(n2_count == r2){
3133	n1_count ++;
3134	n2_count = 1;
3135	data_pos += r3;
3136	}
3137	offset_count_2 = (n1_count + n2_count) % confparams_cpr->sampleDistance;
3138	data_pos += (r3 + confparams_cpr->sampleDistance - offset_count) + (confparams_cpr->sampleDistance - offset_count_2);
3139	offset_count = (confparams_cpr->sampleDistance - offset_count_2);
3140	if(offset_count == 0) offset_count ++;
3141	}
3142	else data_pos += confparams_cpr->sampleDistance;
3143	}
3144	// printf("sample_count: %ld\n", sample_count);
3145	// fflush(stdout);
3146	// if(max_freq < 0.15) max_freq *= 2;
3147	//compute the appropriate number
3148	size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
3149	size_t sum = 0;
3150	for(i=0;i<confparams_cpr->maxRangeRadius;i++)
3151	{
3152	sum += intervals[i];
3153	if(sum>targetCount)
3154	break;
3155	}
3156	if(i>=confparams_cpr->maxRangeRadius)
3157	i = confparams_cpr->maxRangeRadius-1;
3158	unsigned int accIntervals = 2*(i+1);
3159	unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
3160
3161	if(powerOf2<32)
3162	powerOf2 = 32;
3163	free(intervals);
3164	//printf("targetCount=%d, sum=%d, totalSampleSize=%d, ratio=%f, accIntervals=%d, powerOf2=%d\n", targetCount, sum, totalSampleSize, (double)sum/(double)totalSampleSize, accIntervals, powerOf2);
3165	return powerOf2;
3166	}
3167
3168	unsigned int optimize_intervals_double_2D_opt(double *oriData, size_t r1, size_t r2, double realPrecision)
3169	{
3170	size_t i;
3171	size_t radiusIndex;
3172	double pred_value = 0, pred_err;
3173	size_t intervals = (size_t)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
3174	memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
3175	size_t totalSampleSize = 0;//(r1-1)*(r2-1)/confparams_cpr->sampleDistance;
3176
3177	size_t offset_count = confparams_cpr->sampleDistance - 1; // count r2 offset
3178	size_t offset_count_2;
3179	double * data_pos = oriData + r2 + offset_count;
3180	size_t n1_count = 1; // count i sum
3181	size_t len = r1 * r2;
3182	while(data_pos - oriData < len){
3183	totalSampleSize++;
3184	pred_value = data_pos[-1] + data_pos[-r2] - data_pos[-r2-1];
3185	pred_err = fabs(pred_value - *data_pos);
3186	radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
3187	if(radiusIndex>=confparams_cpr->maxRangeRadius)
3188	radiusIndex = confparams_cpr->maxRangeRadius - 1;
3189	intervals[radiusIndex]++;
3190
3191	offset_count += confparams_cpr->sampleDistance;
3192	if(offset_count >= r2){
3193	n1_count ++;
3194	offset_count_2 = n1_count % confparams_cpr->sampleDistance;
3195	data_pos += (r2 + confparams_cpr->sampleDistance - offset_count) + (confparams_cpr->sampleDistance - offset_count_2);
3196	offset_count = (confparams_cpr->sampleDistance - offset_count_2);
3197	if(offset_count == 0) offset_count ++;
3198	}
3199	else data_pos += confparams_cpr->sampleDistance;
3200	}
3201
3202	//compute the appropriate number
3203	size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
3204	size_t sum = 0;
3205	for(i=0;i<confparams_cpr->maxRangeRadius;i++)
3206	{
3207	sum += intervals[i];
3208	if(sum>targetCount)
3209	break;
3210	}
3211	if(i>=confparams_cpr->maxRangeRadius)
3212	i = confparams_cpr->maxRangeRadius-1;
3213	unsigned int accIntervals = 2*(i+1);
3214	unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
3215
3216	if(powerOf2<32)
3217	powerOf2 = 32;
3218
3219	free(intervals);
3220	return powerOf2;
3221	}
3222
3223	unsigned int optimize_intervals_double_1D_opt(double *oriData, size_t dataLength, double realPrecision)
3224	{
3225	size_t i = 0, radiusIndex;
3226	double pred_value = 0, pred_err;
3227	size_t intervals = (size_t)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
3228	memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
3229	size_t totalSampleSize = 0;//dataLength/confparams_cpr->sampleDistance;
3230
3231	double * data_pos = oriData + 2;
3232	while(data_pos - oriData < dataLength){
3233	totalSampleSize++;
3234	//pred_value = 2*data_pos[-1] - data_pos[-2];
3235	pred_value = data_pos[-1];
3236	pred_err = fabs(pred_value - *data_pos);
3237	radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
3238	if(radiusIndex>=confparams_cpr->maxRangeRadius)
3239	radiusIndex = confparams_cpr->maxRangeRadius - 1;
3240	intervals[radiusIndex]++;
3241
3242	data_pos += confparams_cpr->sampleDistance;
3243	}
3244	//compute the appropriate number
3245	size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
3246	size_t sum = 0;
3247	for(i=0;i<confparams_cpr->maxRangeRadius;i++)
3248	{
3249	sum += intervals[i];
3250	if(sum>targetCount)
3251	break;
3252	}
3253	if(i>=confparams_cpr->maxRangeRadius)
3254	i = confparams_cpr->maxRangeRadius-1;
3255
3256	unsigned int accIntervals = 2*(i+1);
3257	unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
3258
3259	if(powerOf2<32)
3260	powerOf2 = 32;
3261
3262	free(intervals);
3263	//printf("accIntervals=%d, powerOf2=%d\n", accIntervals, powerOf2);
3264	return powerOf2;
3265	}

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

Download in other formats: