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

Revision 2c47b73, 5.8 KB checked in by Hal Finkel <hfinkel@…>, 6 years ago (diff)

more work on adding SZ (latest version)

  • Property mode set to 100644
RevLine 
[2c47b73]1/**
2 *  @file sz_double_ts.c
3 *  @author Sheng Di and Dingwen Tao
4 *  @date Aug, 2016
5 *  @brief
6 *  (C) 2016 by Mathematics and Computer Science (MCS), Argonne National Laboratory.
7 *      See COPYRIGHT in top-level directory.
8 */
9
10
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 "zlib.h"
22#include "rw.h"
23#include "sz_double_ts.h"
24
25unsigned int optimize_intervals_double_1D_ts(double *oriData, size_t dataLength, double* preData, double realPrecision)
26{       
27        size_t i = 0, radiusIndex;
28        double pred_value = 0, pred_err;
29        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
30        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
31        size_t totalSampleSize = dataLength/confparams_cpr->sampleDistance;
32        for(i=2;i<dataLength;i++)
33        {
34                if(i%confparams_cpr->sampleDistance==0)
35                {
36                        pred_value = preData[i];
37                        pred_err = fabs(pred_value - oriData[i]);
38                        radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
39                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
40                                radiusIndex = confparams_cpr->maxRangeRadius - 1;                       
41                        intervals[radiusIndex]++;
42                }
43        }
44        //compute the appropriate number
45        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
46        size_t sum = 0;
47        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
48        {
49                sum += intervals[i];
50                if(sum>targetCount)
51                        break;
52        }
53        if(i>=confparams_cpr->maxRangeRadius)
54                i = confparams_cpr->maxRangeRadius-1;
55               
56        unsigned int accIntervals = 2*(i+1);
57        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
58       
59        if(powerOf2<32)
60                powerOf2 = 32;
61       
62        free(intervals);
63        return powerOf2;
64}
65
66TightDataPointStorageD* SZ_compress_double_1D_MDQ_ts(double *oriData, size_t dataLength, sz_multisteps* multisteps,
67double realPrecision, double valueRangeSize, double medianValue_d)
68{
69double* preStepData = (double*)(multisteps->hist_data);
70        //store the decompressed data
71        double* decData = (double*)malloc(sizeof(double)*dataLength);
72        memset(decData, 0, sizeof(double)*dataLength);
73       
74        unsigned int quantization_intervals;
75        if(exe_params->optQuantMode==1)
76                quantization_intervals = optimize_intervals_double_1D_ts(oriData, dataLength, preStepData, realPrecision);
77        else
78                quantization_intervals = exe_params->intvCapacity;
79        updateQuantizationInfo(quantization_intervals); 
80
81        size_t i;
82        int reqLength;
83        double medianValue = medianValue_d;
84        short radExpo = getExponent_double(valueRangeSize/2);
85
86        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);     
87
88        int* type = (int*) malloc(dataLength*sizeof(int));
89               
90        double* spaceFillingValue = oriData; //
91       
92        DynamicIntArray *exactLeadNumArray;
93        new_DIA(&exactLeadNumArray, DynArrayInitLen);
94       
95        DynamicByteArray *exactMidByteArray;
96        new_DBA(&exactMidByteArray, DynArrayInitLen);
97       
98        DynamicIntArray *resiBitArray;
99        new_DIA(&resiBitArray, DynArrayInitLen);
100
101        unsigned char preDataBytes[8];
102        longToBytes_bigEndian(preDataBytes, 0);
103       
104        int reqBytesLength = reqLength/8;
105        int resiBitsLength = reqLength%8;
106
107        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
108        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));                       
109                               
110        //add the first data   
111        type[0] = 0;
112        compressSingleDoubleValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
113        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
114        memcpy(preDataBytes,vce->curBytes,8);
115        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
116        decData[0] = vce->data;
117               
118        //add the second data
119        type[1] = 0;
120        compressSingleDoubleValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
121        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
122        memcpy(preDataBytes,vce->curBytes,8);
123        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
124        decData[1] = vce->data; 
125       
126        int state = 0;
127        double checkRadius = 0;
128        double curData = 0;
129        double pred = 0;
130        double predAbsErr = 0;
131        checkRadius = (exe_params->intvCapacity-1)*realPrecision;
132        double interval = 2*realPrecision;
133
134        for(i=2;i<dataLength;i++)
135        {                               
136                curData = spaceFillingValue[i];
137                pred = preStepData[i];
138                predAbsErr = fabs(curData - pred);     
139                if(predAbsErr<=checkRadius)
140                {
141                        state = (predAbsErr/realPrecision+1)/2;
142                        if(curData>=pred)
143                        {
144                                type[i] = exe_params->intvRadius+state;
145                                pred = pred + state*interval;
146                        }
147                        else //curData<pred
148                        {
149                                type[i] = exe_params->intvRadius-state;
150                                pred = pred - state*interval;
151                        }
152                               
153                        continue;
154                }
155               
156                //unpredictable data processing
157                type[i] = 0;           
158                compressSingleDoubleValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
159                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
160                memcpy(preDataBytes,vce->curBytes,8);
161                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
162                decData[i] = vce->data;
163        }//end of for
164               
165        size_t exactDataNum = exactLeadNumArray->size;
166       
167        TightDataPointStorageD* tdps;
168                       
169        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum, 
170                        type, exactMidByteArray->array, exactMidByteArray->size, 
171                        exactLeadNumArray->array, 
172                        resiBitArray->array, resiBitArray->size, 
173                        resiBitsLength, 
174                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
175       
176        //free memory
177        free_DIA(exactLeadNumArray);
178        free_DIA(resiBitArray);
179        free(type);     
180        free(vce);
181        free(lce);     
182        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
183               
184        memcpy(preStepData, decData, dataLength*sizeof(float)); //update the data
185        free(decData);
186       
187        return tdps;
188}
189
190
Note: See TracBrowser for help on using the repository browser.