1 : #include <stdio.h>
2 : #include <stdlib.h>
3 : #include <string.h>
4 : #include <math.h>
5 : #include "tdlpack.h"
6 : #include "myerror.h"
7 : #include "meta.h"
8 : #include "memendian.h"
9 : #include "fileendian.h"
10 : #include "myassert.h"
11 : #include "myutil.h"
12 : #include "clock.h"
13 : #ifdef MEMWATCH
14 : #include "memwatch.h"
15 : #endif
16 :
17 : #define GRIB_UNSIGN_INT3(a,b,c) ((a<<16)+(b<<8)+c)
18 : #define GRIB_UNSIGN_INT2(a,b) ((a<<8)+b)
19 : #define GRIB_SIGN_INT3(a,b,c) ((1-(int) ((unsigned) (a & 0x80) >> 6)) * (int) (((a & 127) << 16)+(b<<8)+c))
20 : #define GRIB_SIGN_INT2(a,b) ((1-(int) ((unsigned) (a & 0x80) >> 6)) * (int) (((a & 0x7f) << 8) + b))
21 :
22 : /* *INDENT-OFF* */
23 : TDLP_TableType TDLP_B_Table[5] = {
24 : /* 0 */ {0, "Continuous field"},
25 : /* 1 */ {1, "Point Binary - cumulative from above"},
26 : /* 2 */ {2, "Point Binary - cumulative from below"},
27 : /* 3 */ {3, "Point Binary - discrete"},
28 : /* 4 */ {5, "Grid Binary - "},
29 : };
30 :
31 : TDLP_TableType TDLP_DD_Table[9] = {
32 : /* 0 */ {0, "Independent of model"},
33 : /* 1 */ {6, "NGM model"},
34 : /* 2 */ {7, "Eta model"},
35 : /* 3 */ {8, "AVN model"},
36 : /* 4 */ {9, "MRF model"},
37 : /* 5 */ {79, "NDFD forecast"},
38 : /* 6 */ {80, "MOS AEV forecast"},
39 : /* 7 */ {81, "Local AEV Firecasts"},
40 : /* 8 */ {82, "Obs matching AEV Forecasts"},
41 : };
42 :
43 : TDLP_TableType TDLP_V_Table[4] = {
44 : /* 0 */ {0, "No vertical processing"},
45 : /* 1 */ {1, "Difference Levels (UUUU - LLLL)"},
46 : /* 2 */ {2, "Sum Levels (UUUU + LLLL)"},
47 : /* 3 */ {3, "Mean Levels (UUUU + LLLL) / 2."},
48 : };
49 :
50 : TDLP_TableType TDLP_T_Table[3] = {
51 : /* 0 */ {0, "No nolinear tranform"},
52 : /* 1 */ {1, "Square transform"},
53 : /* 2 */ {2, "Square root tranform"},
54 : };
55 :
56 : TDLP_TableType TDLP_Oper_Table[9] = {
57 : /* 0 */ {0, "No time operator"},
58 : /* 1 */ {1, "Mean (Var 1, Var 2)"},
59 : /* 2 */ {2, "Difference (Var 1 - Var 2)"},
60 : /* 3 */ {3, "Maximum (Var 1, Var 2)"},
61 : /* 4 */ {4, "Minimum (Var 1, Var 2)"},
62 : /* 5 */ {5, "Mean (Var 1, Var 2)"},
63 : /* 6 */ {6, "Difference (Var 2 - Var 1)"},
64 : /* 7 */ {7, "Maximum (Var 1, Var 2)"},
65 : /* 8 */ {8, "Minimum (Var 1, Var 2)"},
66 : };
67 :
68 : TDLP_TableType TDLP_I_Table[4] = {
69 : /* 0 */ {0, "No interpolation"},
70 : /* 1 */ {1, "Bi-quadratic interpolation"},
71 : /* 2 */ {2, "Bi-linear interpolation"},
72 : /* 3 */ {3, "Special interpolation for QPF"},
73 : };
74 :
75 : TDLP_TableType TDLP_S_Table[6] = {
76 : /* 0 */ {0, "No smoothing"},
77 : /* 1 */ {1, "5-point smoothing"},
78 : /* 2 */ {2, "9-point smoothing"},
79 : /* 3 */ {3, "25-point smoothing"},
80 : /* 4 */ {4, "81-point smoothing"},
81 : /* 5 */ {5, "168-point smoothing"},
82 : };
83 : /* *INDENT-ON* */
84 :
85 : typedef struct {
86 : sInt4 min; /* Minimum value in a group. Can cause problems if
87 : * this is unsigned. */
88 : uChar bit; /* # of bits in a group. 2^31 is the largest # of
89 : * bits that can be used to hold # of bits in a
90 : * group. However, the # of bits for a number can't
91 : * be > 64, and is probably < 32, so bit is < 6 and
92 : * probably < 5. */
93 : uInt4 num; /* number of values in the group. May need this to be
94 : * signed. */
95 : sInt4 max; /* Max value for the group */
96 : uInt4 start; /* index in Data where group starts. */
97 : uChar f_trySplit; /* Flag, if we have tried to split this group. */
98 : uChar f_tryShift; /* Flag, if we have tried to shift this group. */
99 : } TDLGroupType;
100 :
101 : /*****************************************************************************
102 : * ReadTDLPSect1() --
103 : *
104 : * Arthur Taylor / MDL
105 : *
106 : * PURPOSE
107 : * Parses the TDLP "Product Definition Section" or section 1, filling out
108 : * the pdsMeta data structure.
109 : *
110 : * ARGUMENTS
111 : * pds = The compressed part of the message dealing with "PDS". (Input)
112 : * gribLen = The total length of the TDLP message. (Input)
113 : * curLoc = Current location in the TDLP message. (Output)
114 : * pdsMeta = The filled out pdsMeta data structure. (Output)
115 : * f_gds = boolean if there is a Grid Definition Section. (Output)
116 : * f_bms = boolean if there is a Bitmap Section. (Output)
117 : * DSF = Decimal Scale Factor for unpacking the data. (Output)
118 : * BSF = Binary Scale Factor for unpacking the data. (Output)
119 : *
120 : * FILES/DATABASES: None
121 : *
122 : * RETURNS: int (could use errSprintf())
123 : * 0 = OK
124 : * -1 = gribLen is too small.
125 : *
126 : * HISTORY
127 : * 10/2004 Arthur Taylor (MDL): Created.
128 : *
129 : * NOTES
130 : *****************************************************************************
131 : */
132 0 : static int ReadTDLPSect1 (uChar *pds, sInt4 tdlpLen, sInt4 *curLoc,
133 : pdsTDLPType * pdsMeta, char *f_gds, char *f_bms,
134 : short int *DSF, short int *BSF)
135 : {
136 : char sectLen; /* Length of section. */
137 : sInt4 li_temp; /* Temporary variable. */
138 : int W, XXXX, YY; /* Helps with Threshold calculation. */
139 : int year, t_year; /* The reference year, and a consistency test */
140 : uChar month, t_month; /* The reference month, and a consistency test */
141 : uChar day, t_day; /* The reference day, and a consistency test */
142 : uChar hour, t_hour; /* The reference hour, and a consistency test */
143 : uChar min; /* The reference minute */
144 : uShort2 project_hr; /* The projection in hours. */
145 : sInt4 tau; /* Used to cross check project_hr */
146 : int lenPL; /* Length of the Plain Language descriptor. */
147 :
148 0 : sectLen = *(pds++);
149 0 : *curLoc += sectLen;
150 0 : if (*curLoc > tdlpLen) {
151 0 : errSprintf ("Ran out of data in PDS (TDLP Section 1)\n");
152 0 : return -1;
153 : }
154 : myAssert (sectLen <= 71);
155 0 : if (sectLen < 39) {
156 0 : errSprintf ("TDLP Section 1 is too small.\n");
157 0 : return -1;
158 : }
159 0 : *f_bms = (GRIB2BIT_7 & *pds) ? 1 : 0;
160 0 : *f_gds = (GRIB2BIT_8 & *pds) ? 1 : 0;
161 0 : pds++;
162 0 : year = GRIB_UNSIGN_INT2 (*pds, pds[1]);
163 0 : pds += 2;
164 0 : month = *(pds++);
165 0 : day = *(pds++);
166 0 : hour = *(pds++);
167 0 : min = *(pds++);
168 0 : MEMCPY_BIG (&li_temp, pds, sizeof (sInt4));
169 0 : pds += 4;
170 0 : t_year = li_temp / 1000000L;
171 0 : li_temp -= t_year * 1000000L;
172 0 : t_month = li_temp / 10000L;
173 0 : li_temp -= t_month * 10000L;
174 0 : t_day = li_temp / 100;
175 0 : t_hour = li_temp - t_day * 100;
176 0 : if ((t_year != year) || (t_month != month) || (t_day != day) ||
177 : (t_hour != hour)) {
178 0 : errSprintf ("Error Inconsistant Times in ReadTDLPSect1.\n");
179 0 : return -1;
180 : }
181 0 : if (ParseTime (&(pdsMeta->refTime), year, month, day, hour, min, 0) != 0) {
182 0 : preErrSprintf ("Error In call to ParseTime in ReadTDLPSect1.\n");
183 0 : return -1;
184 : }
185 0 : MEMCPY_BIG (&(li_temp), pds, sizeof (sInt4));
186 0 : pds += 4;
187 0 : pdsMeta->ID1 = li_temp;
188 0 : pdsMeta->CCC = li_temp / 1000000L;
189 0 : li_temp -= pdsMeta->CCC * 1000000L;
190 0 : pdsMeta->FFF = li_temp / 1000;
191 0 : li_temp -= pdsMeta->FFF * 1000;
192 0 : pdsMeta->B = li_temp / 100;
193 0 : pdsMeta->DD = li_temp - pdsMeta->B * 100;
194 0 : MEMCPY_BIG (&(li_temp), pds, sizeof (sInt4));
195 0 : pds += 4;
196 0 : pdsMeta->ID2 = li_temp;
197 0 : pdsMeta->V = li_temp / 100000000L;
198 0 : li_temp -= pdsMeta->V * 100000000L;
199 0 : pdsMeta->LLLL = li_temp / 10000;
200 0 : pdsMeta->UUUU = li_temp - pdsMeta->LLLL * 10000;
201 0 : MEMCPY_BIG (&(li_temp), pds, sizeof (sInt4));
202 0 : pds += 4;
203 0 : pdsMeta->ID3 = li_temp;
204 0 : pdsMeta->T = li_temp / 100000000L;
205 0 : li_temp -= pdsMeta->T * 100000000L;
206 0 : pdsMeta->RR = li_temp / 1000000L;
207 0 : li_temp -= pdsMeta->RR * 1000000L;
208 0 : pdsMeta->Oper = li_temp / 100000L;
209 0 : li_temp -= pdsMeta->Oper * 100000L;
210 0 : pdsMeta->HH = li_temp / 1000;
211 0 : pdsMeta->ttt = li_temp - pdsMeta->HH * 1000;
212 0 : MEMCPY_BIG (&(li_temp), pds, sizeof (sInt4));
213 0 : pds += 4;
214 0 : pdsMeta->ID4 = li_temp;
215 0 : W = li_temp / 1000000000L;
216 0 : li_temp -= W * 1000000000L;
217 0 : XXXX = li_temp / 100000L;
218 0 : li_temp -= XXXX * 100000L;
219 0 : if (W) {
220 0 : XXXX = -1 * XXXX;
221 : }
222 0 : YY = li_temp / 1000L;
223 0 : li_temp -= YY * 1000L;
224 0 : if (YY >= 50) {
225 0 : YY = -1 * (YY - 50);
226 : }
227 0 : pdsMeta->thresh = (XXXX / 10000.) * pow (10.0, YY);
228 0 : pdsMeta->I = li_temp / 100;
229 0 : li_temp -= pdsMeta->I * 100L;
230 0 : pdsMeta->S = li_temp / 10;
231 0 : pdsMeta->G = li_temp - pdsMeta->S * 10;
232 0 : project_hr = GRIB_UNSIGN_INT2 (*pds, pds[1]);
233 0 : tau = pdsMeta->ID3 - ((pdsMeta->ID3 / 1000) * 1000);
234 0 : if (tau != project_hr) {
235 : printf ("Warning: Inconsistant Projections in hours in "
236 0 : "ReadTDLPSect1 (%d vs %d)\n", tau, project_hr);
237 : /*
238 : errSprintf ("Warning: Inconsistant Projections in hours in "
239 : "ReadTDLPSect1 (%ld vs %d)\n", tau, project_hr);
240 : */
241 0 : project_hr = tau;
242 : /* return -1; */
243 : }
244 0 : pds += 2;
245 0 : pdsMeta->project = project_hr * 3600 + (*(pds++)) * 60;
246 0 : pdsMeta->procNum = (*(pds++));
247 0 : pdsMeta->seqNum = (*(pds++));
248 0 : *DSF = (*pds > 128) ? 128 - (*(pds++)) : (*(pds++));
249 0 : *BSF = (*pds > 128) ? 128 - (*(pds++)) : (*(pds++));
250 0 : if ((*pds != 0) || (pds[1] != 0) || (pds[2] != 0)) {
251 0 : errSprintf ("Error Reserved was not set to 0 in ReadTDLPSect1.\n");
252 0 : return -1;
253 : }
254 0 : pds += 3;
255 0 : lenPL = (*(pds++));
256 0 : if (sectLen - lenPL != 39) {
257 : errSprintf ("Error sectLen(%d) - lenPL(%d) != 39 in ReadTDLPSect1.\n",
258 0 : sectLen, lenPL);
259 0 : return -1;
260 : }
261 0 : if (lenPL > 32) {
262 0 : lenPL = 32;
263 : }
264 0 : strncpy (pdsMeta->Descriptor, (char *) pds, lenPL);
265 0 : pdsMeta->Descriptor[lenPL] = '\0';
266 0 : strTrim (pdsMeta->Descriptor);
267 0 : return 0;
268 : }
269 :
270 : /*****************************************************************************
271 : * TDLP_TableLookUp() --
272 : *
273 : * Arthur Taylor / MDL
274 : *
275 : * PURPOSE
276 : * To look up TDL Ids information in a given table.
277 : *
278 : * ARGUMENTS
279 : * table = The Table to look up the Id in. (Input)
280 : * tableLen = The length of the Table. (Input)
281 : * index = The index into the Table. (Input)
282 : *
283 : * FILES/DATABASES: None
284 : *
285 : * RETURNS: char *
286 : * Returns the meta data that is associated with index in the table.
287 : *
288 : * HISTORY
289 : * 10/2004 Arthur Taylor (MDL): Created.
290 : *
291 : * NOTES
292 : *****************************************************************************
293 : */
294 0 : static const char *TDLP_TableLookUp(TDLP_TableType * table, int tableLen,
295 : int index)
296 : {
297 : int i; /* Loop counter. */
298 :
299 0 : for (i = 0; i < tableLen; i++) {
300 0 : if (table[i].index == index)
301 0 : return (table[i].data);
302 : }
303 0 : return "Unknown";
304 : }
305 :
306 : /*****************************************************************************
307 : * PrintPDS_TDLP() --
308 : *
309 : * Arthur Taylor / MDL
310 : *
311 : * PURPOSE
312 : * To generate the message for the Product Definition Sections of the TDLP
313 : * Message.
314 : *
315 : * ARGUMENTS
316 : * pds = The TDLP Product Definition Section to print. (Input)
317 : *
318 : * FILES/DATABASES: None
319 : *
320 : * RETURNS: void
321 : *
322 : * HISTORY
323 : * 10/2004 Arthur Taylor (MDL): Created.
324 : *
325 : * NOTES
326 : *****************************************************************************
327 : */
328 0 : void PrintPDS_TDLP (pdsTDLPType * pds)
329 : {
330 : char buffer[25]; /* Stores format of pds1->refTime. */
331 :
332 : /*
333 : strftime (buffer, 25, "%m/%d/%Y %H:%M:%S UTC", gmtime (&(pds->refTime)));
334 : */
335 0 : Clock_Print (buffer, 25, pds->refTime, "%m/%d/%Y %H:%M:%S UTC", 0);
336 :
337 0 : Print ("PDS-TDLP", "Reference Time", Prt_S, buffer);
338 0 : Print ("PDS-TDLP", "Plain Language", Prt_S, pds->Descriptor);
339 0 : sprintf (buffer, "%09d", pds->ID1);
340 0 : Print ("PDS-TDLP", "ID 1", Prt_S, buffer);
341 0 : sprintf (buffer, "%09d", pds->ID2);
342 0 : Print ("PDS-TDLP", "ID 2", Prt_S, buffer);
343 0 : sprintf (buffer, "%09d", pds->ID3);
344 0 : Print ("PDS-TDLP", "ID 3", Prt_S, buffer);
345 0 : Print ("PDS-TDLP", "ID 4", Prt_D, pds->ID4);
346 0 : Print ("PDS-TDLP", "Model or Process Number", Prt_D, pds->procNum);
347 0 : Print ("PDS-TDLP", "Sequence Number", Prt_D, pds->seqNum);
348 :
349 0 : sprintf (buffer, "%03d", pds->CCC);
350 0 : Print ("PDS-TDLP", "ID1-CCC", Prt_S, buffer);
351 0 : sprintf (buffer, "%03d", pds->FFF);
352 0 : Print ("PDS-TDLP", "ID1-FFF", Prt_S, buffer);
353 : Print ("PDS-TDLP", "ID1-B", Prt_DS, pds->B,
354 0 : TDLP_TableLookUp (TDLP_B_Table, sizeof (TDLP_B_Table), pds->B));
355 0 : sprintf (buffer, "%02d", pds->DD);
356 : Print ("PDS-TDLP", "ID1-DD", Prt_SS, buffer,
357 0 : TDLP_TableLookUp (TDLP_DD_Table, sizeof (TDLP_DD_Table), pds->DD));
358 :
359 : Print ("PDS-TDLP", "ID2-V", Prt_DS, pds->V,
360 0 : TDLP_TableLookUp (TDLP_V_Table, sizeof (TDLP_V_Table), pds->V));
361 0 : sprintf (buffer, "%04d", pds->LLLL);
362 0 : Print ("PDS-TDLP", "ID2-LLLL", Prt_S, buffer);
363 0 : sprintf (buffer, "%04d", pds->UUUU);
364 0 : Print ("PDS-TDLP", "ID2-UUUU", Prt_S, buffer);
365 :
366 0 : if (pds->Oper != 0) {
367 : Print ("PDS-TDLP", "ID3-T", Prt_DS, pds->T,
368 0 : TDLP_TableLookUp (TDLP_T_Table, sizeof (TDLP_T_Table), pds->T));
369 0 : sprintf (buffer, "%02d", pds->RR);
370 : Print ("PDS-TDLP", "ID3-RR", Prt_SS, buffer,
371 0 : "Run time offset in hours");
372 : Print ("PDS-TDLP", "ID3-Oper", Prt_DS, pds->Oper,
373 : TDLP_TableLookUp (TDLP_Oper_Table, sizeof (TDLP_Oper_Table),
374 0 : pds->Oper));
375 0 : sprintf (buffer, "%02d", pds->HH);
376 : Print ("PDS-TDLP", "ID3-HH", Prt_SS, buffer,
377 0 : "Number of hours between variables");
378 : } else {
379 : Print ("PDS-TDLP", "ID3-Oper", Prt_DS, pds->Oper,
380 : TDLP_TableLookUp (TDLP_Oper_Table, sizeof (TDLP_Oper_Table),
381 0 : pds->Oper));
382 : }
383 0 : sprintf (buffer, "%03d", pds->ttt);
384 0 : Print ("PDS-TDLP", "ID3-ttt", Prt_SS, buffer, "Forecast Projection");
385 :
386 0 : Print ("PDS-TDLP", "ID4-thresh", Prt_F, pds->thresh);
387 : Print ("PDS-TDLP", "ID4-I", Prt_DS, pds->I,
388 0 : TDLP_TableLookUp (TDLP_I_Table, sizeof (TDLP_I_Table), pds->I));
389 : Print ("PDS-TDLP", "ID4-S", Prt_DS, pds->S,
390 0 : TDLP_TableLookUp (TDLP_S_Table, sizeof (TDLP_S_Table), pds->S));
391 0 : Print ("PDS-TDLP", "ID4-G", Prt_D, pds->G);
392 0 : }
393 :
394 : /*****************************************************************************
395 : * TDLP_ElemSurfUnit() --
396 : *
397 : * Arthur Taylor / MDL
398 : *
399 : * PURPOSE
400 : * Deal with element name, unit, and comment for TDLP items. Also deals
401 : * with the surface information stored in ID2.
402 : *
403 : * ARGUMENTS
404 : * pds = The TDLP Product Definition Section to parse. (Input)
405 : * element = The resulting element name. (Output)
406 : * unitName = The resulting unit name. (Output)
407 : * comment = The resulting comment. (Output)
408 : * shortFstLevel = The resulting short forecast level. (Output)
409 : * longFstLevel = The resulting long forecast level. (Output)
410 : *
411 : * FILES/DATABASES: None
412 : *
413 : * RETURNS: void
414 : *
415 : * HISTORY
416 : * 10/2004 Arthur Taylor (MDL): Created.
417 : *
418 : * NOTES
419 : *****************************************************************************
420 : */
421 0 : static void TDLP_ElemSurfUnit (pdsTDLPType * pds, char **element,
422 : char **unitName, char **comment,
423 : char **shortFstLevel, char **longFstLevel)
424 : {
425 : char *ptr; /* Help guess unitName, and clean the elemName. */
426 : char *ptr2; /* Help guess unitName, and clean the elemName. */
427 :
428 : myAssert (*element == NULL);
429 : myAssert (*unitName == NULL);
430 : myAssert (*comment == NULL);
431 : myAssert (*shortFstLevel == NULL);
432 : myAssert (*longFstLevel == NULL);
433 :
434 0 : *element = (char *) malloc (1 + strlen (pds->Descriptor) * sizeof (char));
435 0 : strcpy (*element, pds->Descriptor);
436 0 : (*element)[strlen (pds->Descriptor)] = '\0';
437 :
438 0 : ptr = strchr (*element, '(');
439 0 : if (ptr != NULL) {
440 0 : ptr2 = strchr (ptr, ')');
441 0 : *ptr2 = '\0';
442 0 : if (strcmp (ptr + 1, "unofficial id") == 0) {
443 0 : *unitName = (char *) malloc ((1 + 3) * sizeof (char));
444 0 : strcpy (*unitName, "[-]");
445 : } else {
446 0 : reallocSprintf (unitName, "[%s]", ptr + 1);
447 : }
448 : /* Trim any parens from element. */
449 0 : *ptr = '\0';
450 0 : strTrimRight (*element, ' ');
451 : } else {
452 0 : *unitName = (char *) malloc ((1 + 3) * sizeof (char));
453 0 : strcpy (*unitName, "[-]");
454 : }
455 0 : ptr = *element;
456 0 : while (*ptr != '\0') {
457 0 : if (*ptr == ' ') {
458 0 : *ptr = '-';
459 : }
460 0 : ptr++;
461 : }
462 0 : strCompact (*element, '-');
463 :
464 : reallocSprintf (comment, "%09ld-%09ld-%09ld-%ld %s", pds->ID1,
465 0 : pds->ID2, pds->ID3, pds->ID4, *unitName);
466 0 : reallocSprintf (shortFstLevel, "%09ld", pds->ID2);
467 0 : reallocSprintf (longFstLevel, "%09ld", pds->ID2);
468 0 : }
469 :
470 : /*****************************************************************************
471 : * TDLP_Inventory() --
472 : *
473 : * Arthur Taylor / MDL
474 : *
475 : * PURPOSE
476 : * Parses the TDLP "Product Definition Section" for enough information to
477 : * fill out the inventory data structure so we can do a simple inventory on
478 : * the file in a similar way to how we did it for GRIB1 and GRIB2.
479 : *
480 : * ARGUMENTS
481 : * fp = An opened TDLP file already at the correct message. (Input)
482 : * tdlpLen = The total length of the TDLP message. (Input)
483 : * inv = The inventory data structure that we need to fill. (Output)
484 : *
485 : * FILES/DATABASES: None
486 : *
487 : * RETURNS: int (could use errSprintf())
488 : * 0 = OK
489 : * -1 = tdlpLen is too small.
490 : *
491 : * HISTORY
492 : * 10/2004 Arthur Taylor (MDL): Created.
493 : *
494 : * NOTES
495 : * Speed improvements...
496 : * 1) pds doen't need to be allocated each time.
497 : * 2) Not all data is needed, do something like TDLP_RefTime
498 : * 3) TDLP_ElemSurfUnit may be slow?
499 : *****************************************************************************
500 : */
501 0 : int TDLP_Inventory (DataSource &fp, sInt4 tdlpLen, inventoryType *inv)
502 : {
503 : sInt4 curLoc; /* Where we are in the current TDLP message. */
504 : int i_temp;
505 : uChar sectLen; /* Length of section. */
506 : uChar *pds; /* The part of the message dealing with the PDS. */
507 : pdsTDLPType pdsMeta; /* The pds parsed into a usable data structure. */
508 : char f_gds; /* flag if there is a GDS section. */
509 : char f_bms; /* flag if there is a BMS section. */
510 : short int DSF; /* Decimal Scale Factor for unpacking the data. */
511 : short int BSF; /* Binary Scale Factor for unpacking the data. */
512 :
513 0 : curLoc = 8;
514 0 : if ((i_temp = fp.DataSourceFgetc()) == EOF) {
515 0 : errSprintf ("Ran out of file in PDS (TDLP_Inventory).\n");
516 0 : return -1;
517 : }
518 0 : sectLen = (uChar) i_temp;
519 0 : curLoc += sectLen;
520 0 : if (curLoc > tdlpLen) {
521 0 : errSprintf ("Ran out of data in PDS (TDLP_Inventory)\n");
522 0 : return -1;
523 : }
524 0 : pds = (uChar *) malloc (sectLen * sizeof (uChar));
525 0 : *pds = sectLen;
526 0 : if (fp.DataSourceFread (pds + 1, sizeof (char), sectLen - 1) + 1 != sectLen) {
527 0 : errSprintf ("Ran out of file.\n");
528 0 : free (pds);
529 0 : return -1;
530 : }
531 :
532 0 : if (ReadTDLPSect1 (pds, tdlpLen, &curLoc, &pdsMeta, &f_gds, &f_bms,
533 : &DSF, &BSF) != 0) {
534 0 : preErrSprintf ("Inside TDLP_Inventory\n");
535 0 : free (pds);
536 0 : return -1;
537 : }
538 0 : free (pds);
539 :
540 0 : inv->element = NULL;
541 0 : inv->unitName = NULL;
542 0 : inv->comment = NULL;
543 0 : free (inv->shortFstLevel);
544 0 : inv->shortFstLevel = NULL;
545 0 : free (inv->longFstLevel);
546 0 : inv->longFstLevel = NULL;
547 : TDLP_ElemSurfUnit (&pdsMeta, &(inv->element), &(inv->unitName),
548 : &(inv->comment), &(inv->shortFstLevel),
549 0 : &(inv->longFstLevel));
550 :
551 0 : inv->refTime = pdsMeta.refTime;
552 0 : inv->validTime = pdsMeta.refTime + pdsMeta.project;
553 0 : inv->foreSec = pdsMeta.project;
554 0 : return 0;
555 : }
556 :
557 : /*****************************************************************************
558 : * TDLP_RefTime() --
559 : *
560 : * Arthur Taylor / MDL
561 : *
562 : * PURPOSE
563 : * Return just the reference time of a TDLP message.
564 : *
565 : * ARGUMENTS
566 : * fp = An opened TDLP file already at the correct message. (Input)
567 : * tdlpLen = The total length of the TDLP message. (Input)
568 : * refTime = The reference time of interest. (Output)
569 : *
570 : * FILES/DATABASES: None
571 : *
572 : * RETURNS: int (could use errSprintf())
573 : * 0 = OK
574 : * -1 = tdlpLen is too small.
575 : *
576 : * HISTORY
577 : * 10/2004 Arthur Taylor (MDL): Created.
578 : *
579 : * NOTES
580 : *****************************************************************************
581 : */
582 0 : int TDLP_RefTime (DataSource &fp, sInt4 tdlpLen, double *refTime)
583 : {
584 : int sectLen; /* Length of section. */
585 : sInt4 curLoc; /* Where we are in the current TDLP message. */
586 : int c_temp; /* Temporary variable for use with fgetc */
587 : short int si_temp; /* Temporary variable. */
588 : int year, t_year; /* The reference year, and a consistency test */
589 : uChar month, t_month; /* The reference month, and a consistency test */
590 : uChar day, t_day; /* The reference day, and a consistency test */
591 : uChar hour, t_hour; /* The reference hour, and a consistency test */
592 : uChar min; /* The reference minute */
593 : sInt4 li_temp; /* Temporary variable. */
594 :
595 0 : if ((sectLen = fp.DataSourceFgetc ()) == EOF)
596 0 : goto error;
597 0 : curLoc = 8 + sectLen;
598 0 : if (curLoc > tdlpLen) {
599 0 : errSprintf ("Ran out of data in PDS (TDLP_RefTime)\n");
600 0 : return -1;
601 : }
602 : myAssert (sectLen <= 71);
603 0 : if (sectLen < 39) {
604 0 : errSprintf ("TDLP Section 1 is too small.\n");
605 0 : return -1;
606 : }
607 :
608 0 : if ((c_temp = fp.DataSourceFgetc()) == EOF)
609 0 : goto error;
610 0 : if (FREAD_BIG (&si_temp, sizeof (short int), 1, fp) != 1)
611 0 : goto error;
612 0 : year = si_temp;
613 0 : if ((c_temp = fp.DataSourceFgetc()) == EOF)
614 0 : goto error;
615 0 : month = c_temp;
616 0 : if ((c_temp = fp.DataSourceFgetc()) == EOF)
617 0 : goto error;
618 0 : day = c_temp;
619 0 : if ((c_temp = fp.DataSourceFgetc()) == EOF)
620 0 : goto error;
621 0 : hour = c_temp;
622 0 : if ((c_temp = fp.DataSourceFgetc()) == EOF)
623 0 : goto error;
624 0 : min = c_temp;
625 :
626 0 : if (FREAD_BIG (&li_temp, sizeof (sInt4), 1, fp) != 1)
627 0 : goto error;
628 0 : t_year = li_temp / 1000000L;
629 0 : li_temp -= t_year * 1000000L;
630 0 : t_month = li_temp / 10000L;
631 0 : li_temp -= t_month * 10000L;
632 0 : t_day = li_temp / 100;
633 0 : t_hour = li_temp - t_day * 100;
634 :
635 0 : if ((t_year != year) || (t_month != month) || (t_day != day) ||
636 : (t_hour != hour)) {
637 0 : errSprintf ("Error Inconsistant Times in TDLP_RefTime.\n");
638 0 : return -1;
639 : }
640 0 : if (ParseTime (refTime, year, month, day, hour, min, 0) != 0) {
641 0 : preErrSprintf ("Error In call to ParseTime in TDLP_RefTime.\n");
642 0 : return -1;
643 : }
644 :
645 : /* Get to the end of the TDLP message. */
646 : /* (inventory.c : GRIB2Inventory), is responsible for this. */
647 : /* fseek (fp, gribLen - sectLen, SEEK_CUR); */
648 0 : return 0;
649 : error:
650 0 : errSprintf ("Ran out of file in PDS (TDLP_RefTime).\n");
651 0 : return -1;
652 : }
653 :
654 : /*****************************************************************************
655 : * ReadTDLPSect2() --
656 : *
657 : * Arthur Taylor / MDL
658 : *
659 : * PURPOSE
660 : * Parses the TDLP "Grid Definition Section" or section 2, filling out
661 : * the gdsMeta data structure.
662 : *
663 : * ARGUMENTS
664 : * gds = The compressed part of the message dealing with "GDS". (Input)
665 : * tdlpLen = The total length of the TDLP message. (Input)
666 : * curLoc = Current location in the TDLP message. (Output)
667 : * gdsMeta = The filled out gdsMeta data structure. (Output)
668 : *
669 : * FILES/DATABASES: None
670 : *
671 : * RETURNS: int (could use errSprintf())
672 : * 0 = OK
673 : * -1 = tdlpLen is too small.
674 : * -2 = unexpected values in gds.
675 : *
676 : * HISTORY
677 : * 10/2004 Arthur Taylor (MDL): Created.
678 : *
679 : * NOTES
680 : *****************************************************************************
681 : */
682 0 : static int ReadTDLPSect2 (uChar *gds, sInt4 tdlpLen, sInt4 *curLoc,
683 : gdsType *gdsMeta)
684 : {
685 : char sectLen; /* Length of section. */
686 : int gridType; /* Which type of grid. (see enumerated types). */
687 : sInt4 li_temp; /* Temporary variable. */
688 :
689 0 : sectLen = *(gds++);
690 0 : *curLoc += sectLen;
691 0 : if (*curLoc > tdlpLen) {
692 0 : errSprintf ("Ran out of data in GDS (TDLP Section 2)\n");
693 0 : return -1;
694 : }
695 : myAssert (sectLen == 28);
696 :
697 0 : gridType = *(gds++);
698 0 : gdsMeta->Nx = GRIB_UNSIGN_INT2 (*gds, gds[1]);
699 0 : gds += 2;
700 0 : gdsMeta->Ny = GRIB_UNSIGN_INT2 (*gds, gds[1]);
701 0 : gds += 2;
702 0 : gdsMeta->lat1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) / 10000.0;
703 0 : gds += 3;
704 0 : gdsMeta->lon1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) / 10000.0;
705 0 : gdsMeta->lon1 = 360 - gdsMeta->lon1;
706 0 : if (gdsMeta->lon1 < 0) {
707 0 : gdsMeta->lon1 += 360;
708 : }
709 0 : if (gdsMeta->lon1 > 360) {
710 0 : gdsMeta->lon1 -= 360;
711 : }
712 0 : gds += 3;
713 0 : gdsMeta->orientLon = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) / 10000.0;
714 0 : gdsMeta->orientLon = 360 - gdsMeta->orientLon;
715 0 : if (gdsMeta->orientLon < 0) {
716 0 : gdsMeta->orientLon += 360;
717 : }
718 0 : if (gdsMeta->orientLon > 360) {
719 0 : gdsMeta->orientLon -= 360;
720 : }
721 0 : gds += 3;
722 0 : MEMCPY_BIG (&li_temp, gds, sizeof (sInt4));
723 0 : gdsMeta->Dx = li_temp / 1000.0;
724 0 : gds += 4;
725 0 : gdsMeta->meshLat = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) / 10000.0;
726 0 : gds += 3;
727 0 : if ((*gds != 0) || (gds[1] != 0) || (gds[2] != 0) || (gds[3] != 0) ||
728 0 : (gds[4] != 0) || (gds[5] != 0)) {
729 0 : errSprintf ("Error Reserved was not set to 0 in ReadTDLPSect2.\n");
730 0 : return -1;
731 : }
732 :
733 0 : gdsMeta->numPts = gdsMeta->Nx * gdsMeta->Ny;
734 0 : gdsMeta->f_sphere = 1;
735 0 : gdsMeta->majEarth = 6371.2;
736 0 : gdsMeta->minEarth = 6371.2;
737 0 : gdsMeta->Dy = gdsMeta->Dx;
738 0 : gdsMeta->resFlag = 0;
739 0 : gdsMeta->center = 0;
740 0 : gdsMeta->scan = 64;
741 0 : gdsMeta->lat2 = 0;
742 0 : gdsMeta->lon2 = 0;
743 0 : gdsMeta->scaleLat1 = gdsMeta->meshLat;
744 0 : gdsMeta->scaleLat2 = gdsMeta->meshLat;
745 0 : gdsMeta->southLat = 0;
746 0 : gdsMeta->southLon = 0;
747 0 : switch (gridType) {
748 : case TDLP_POLAR:
749 0 : gdsMeta->projType = GS3_POLAR;
750 : /* 4/24/2006 Added the following. */
751 0 : gdsMeta->scaleLat1 = 90;
752 0 : gdsMeta->scaleLat2 = 90;
753 0 : break;
754 : case TDLP_LAMBERT:
755 0 : gdsMeta->projType = GS3_LAMBERT;
756 0 : break;
757 : case TDLP_MERCATOR:
758 0 : gdsMeta->projType = GS3_MERCATOR;
759 : /* 4/24/2006 Added the following. */
760 0 : gdsMeta->scaleLat1 = 0;
761 0 : gdsMeta->scaleLat2 = 0;
762 0 : break;
763 : default:
764 0 : errSprintf ("Grid projection number is %d\n", gridType);
765 0 : return -2;
766 : }
767 0 : return 0;
768 : }
769 :
770 : /*****************************************************************************
771 : * ReadTDLPSect3() --
772 : *
773 : * Arthur Taylor / MDL
774 : *
775 : * PURPOSE
776 : * Parses the TDLP "Bit Map Section" or section 3, filling out the bitmap
777 : * as needed.
778 : *
779 : * ARGUMENTS
780 : * bms = The compressed part of the message dealing with "BMS". (Input)
781 : * tdlpLen = The total length of the TDLP message. (Input)
782 : * curLoc = Current location in the TDLP message. (Output)
783 : * bitmap = The extracted bitmap. (Output)
784 : * NxNy = The total size of the grid. (Input)
785 : *
786 : * FILES/DATABASES: None
787 : *
788 : * RETURNS: int (could use errSprintf())
789 : * 0 = OK
790 : * -1 = tdlpLen is too small.
791 : * -2 = unexpected values in bms.
792 : *
793 : * HISTORY
794 : * 10/2004 Arthur Taylor (MDL): Created.
795 : *
796 : * NOTES
797 : *****************************************************************************
798 : */
799 0 : static int ReadTDLPSect3 (uChar *bms, sInt4 tdlpLen, sInt4 *curLoc,
800 : uChar *bitmap, sInt4 NxNy)
801 : {
802 0 : errSprintf ("Bitmap data is Not Supported\n");
803 0 : return -1;
804 : }
805 :
806 : /*****************************************************************************
807 : * ReadTDLPSect4() --
808 : *
809 : * Arthur Taylor / MDL
810 : *
811 : * PURPOSE
812 : * Unpacks the "Binary Data Section" or section 4.
813 : *
814 : * ARGUMENTS
815 : * bds = The compressed part of the message dealing with "BDS". (Input)
816 : * tdlpLen = The total length of the TDLP message. (Input)
817 : * curLoc = Current location in the TDLP message. (Output)
818 : * DSF = Decimal Scale Factor for unpacking the data. (Input)
819 : * BSF = Binary Scale Factor for unpacking the data. (Input)
820 : * data = The extracted grid. (Output)
821 : * meta = The meta data associated with the grid (Input/Output)
822 : * unitM = The M unit conversion value in equation y = Mx + B. (Input)
823 : * unitB = The B unit conversion value in equation y = Mx + B. (Input)
824 : *
825 : * FILES/DATABASES: None
826 : *
827 : * RETURNS: int (could use errSprintf())
828 : * 0 = OK
829 : * -1 = gribLen is too small.
830 : * -2 = unexpected values in bds.
831 : *
832 : * HISTORY
833 : * 4/2003 Arthur Taylor (MDL/RSIS): Created
834 : * 3/2004 AAT: Switched {# Pts * (# Bits in a Group) +
835 : * # of unused bits != # of available bits} to a warning from an
836 : * error.
837 : * 2/2005 AAT: Found bug: memBitRead grp[i].bit was sizeof sInt4 instead
838 : * of uChar.
839 : * 2/2005 AAT: Second order diff, no miss value bug (lastData - 1) should
840 : * be lastData.
841 : * 2/2005 AAT: Added test to see if the number of bits needed matches the
842 : * section length.
843 : *
844 : * NOTES
845 : * 1) See metaparse.c : ParseGrid()
846 : *****************************************************************************
847 : */
848 0 : static int ReadTDLPSect4 (uChar *bds, sInt4 tdlpLen, sInt4 *curLoc,
849 : short int DSF, short int BSF, double *data,
850 : grib_MetaData *meta, double unitM, double unitB)
851 : {
852 : uInt4 sectLen; /* Length in bytes of the current section. */
853 : uChar f_notGridPnt; /* Not Grid point data? */
854 : uChar f_complexPack; /* Complex packing? */
855 : uChar f_sndOrder; /* Second order differencing? */
856 : uChar f_primMiss; /* Primary missing value? */
857 : uChar f_secMiss; /* Secondary missing value? */
858 : uInt4 numPack; /* Number of points packed. */
859 : sInt4 li_temp; /* Temporary variable. */
860 : uInt4 uli_temp; /* Temporary variable. */
861 : uChar bufLoc; /* Keeps track of where to start getting more data
862 : * out of the packed data stream. */
863 : uChar f_negative; /* used to help with signs of numbers. */
864 : size_t numUsed; /* How many bytes were used in a given call to
865 : * memBitRead. */
866 0 : sInt4 origVal = 0; /* Original value. */
867 : uChar mbit; /* # of bits for abs (first first order difference) */
868 0 : sInt4 fstDiff = 0; /* First first order difference. */
869 0 : sInt4 diff = 0; /* general first order difference. */
870 : uChar nbit; /* # of bits for abs (overall min value) */
871 : sInt4 minVal; /* Minimum value. */
872 : size_t LX; /* Number of groups. */
873 : uChar ibit; /* # of bits for group min values. */
874 : uChar jbit; /* # of bits for # of bits for group. */
875 : uChar kbit; /* # of bits for # values in a group. */
876 : TDLGroupType *grp; /* Holds the info about each group. */
877 : size_t i, j; /* Loop counters. */
878 : uInt4 t_numPack; /* Used to total number of values in a group to check
879 : * the numPack value. */
880 : uInt4 t_numBits; /* Used to total number of bits used in the groups. */
881 : uInt4 t_numBytes; /* Used to total number of bytes used to compare to
882 : * sectLen. */
883 : sInt4 maxVal; /* The max value in a group. */
884 : uInt4 dataCnt; /* How many values (miss or othewise) we have read. */
885 : uInt4 lastData; /* Index to last actual data. */
886 : uInt4 numVal; /* # of actual (non-missing values) we have. */
887 : double scale; /* Amount to scale values by. */
888 : uInt4 dataInd; /* Index into data for this value (used to switch
889 : * from a11..a1n,a2n..a21 to normal grid of
890 : * a11..a1n,a21..a2n. */
891 : uChar f_missing; /* Used to help with primary missing values, and the
892 : * 0 bit possibility. */
893 : #ifdef DEBUG
894 : sInt4 t_UK1 = 0; /* Used to test theories about un defined values. */
895 : sInt4 t_UK2 = 0; /* Used to test theories about un defined values. */
896 : #endif
897 :
898 0 : sectLen = GRIB_UNSIGN_INT3 (*bds, bds[1], bds[2]);
899 0 : *curLoc += sectLen;
900 0 : if (*curLoc > tdlpLen) {
901 0 : errSprintf ("Ran out of data in BDS (TDLP Section 4)\n");
902 0 : return -1;
903 : }
904 0 : bds += 3;
905 0 : t_numBytes = 3;
906 0 : f_notGridPnt = (GRIB2BIT_4 & *bds) ? 1 : 0;
907 0 : f_complexPack = (GRIB2BIT_5 & *bds) ? 1 : 0;
908 0 : f_sndOrder = (GRIB2BIT_6 & *bds) ? 1 : 0;
909 0 : f_primMiss = (GRIB2BIT_7 & *bds) ? 1 : 0;
910 0 : f_secMiss = (GRIB2BIT_8 & *bds) ? 1 : 0;
911 :
912 0 : if (f_secMiss && (!f_primMiss)) {
913 0 : errSprintf ("Secondary missing value without a primary!\n");
914 0 : return -1;
915 : }
916 0 : if (f_complexPack) {
917 0 : if (!f_sndOrder) {
918 0 : meta->gridAttrib.packType = GS5_CMPLX;
919 : } else {
920 0 : meta->gridAttrib.packType = GS5_CMPLXSEC;
921 : }
922 : } else {
923 0 : errSprintf ("Simple pack is not supported at this time.\n");
924 0 : return -1;
925 : }
926 0 : bds++;
927 0 : MEMCPY_BIG (&numPack, bds, sizeof (sInt4));
928 0 : bds += 4;
929 0 : t_numBytes += 5;
930 0 : if (!f_notGridPnt) {
931 0 : if (numPack != meta->gds.numPts) {
932 : errSprintf ("Number packed %d != number of points %d\n", numPack,
933 0 : meta->gds.numPts);
934 0 : return -1;
935 : }
936 : }
937 0 : meta->gridAttrib.DSF = DSF;
938 0 : meta->gridAttrib.ESF = BSF;
939 0 : meta->gridAttrib.fieldType = 0;
940 0 : if (!f_primMiss) {
941 0 : meta->gridAttrib.f_miss = 0;
942 0 : meta->gridAttrib.missPri = 0;
943 0 : meta->gridAttrib.missSec = 0;
944 : } else {
945 0 : MEMCPY_BIG (&li_temp, bds, sizeof (sInt4));
946 0 : bds += 4;
947 0 : t_numBytes += 4;
948 0 : meta->gridAttrib.missPri = li_temp / 10000.0;
949 0 : if (!f_secMiss) {
950 0 : meta->gridAttrib.f_miss = 1;
951 0 : meta->gridAttrib.missSec = 0;
952 : } else {
953 0 : MEMCPY_BIG (&li_temp, bds, sizeof (sInt4));
954 0 : bds += 4;
955 0 : t_numBytes += 4;
956 0 : meta->gridAttrib.missSec = li_temp / 10000.0;
957 0 : meta->gridAttrib.f_miss = 2;
958 : }
959 : }
960 : /* Init the buffer location. */
961 0 : bufLoc = 8;
962 : /* The origValue and fstDiff are only present if sndOrder packed. */
963 0 : if (f_sndOrder) {
964 : memBitRead (&f_negative, sizeof (f_negative), bds, 1, &bufLoc,
965 0 : &numUsed);
966 0 : memBitRead (&uli_temp, sizeof (sInt4), bds, 31, &bufLoc, &numUsed);
967 : myAssert (numUsed == 4);
968 0 : bds += numUsed;
969 0 : t_numBytes += numUsed;
970 0 : origVal = (f_negative) ? -1 * uli_temp : uli_temp;
971 0 : memBitRead (&mbit, sizeof (mbit), bds, 5, &bufLoc, &numUsed);
972 : memBitRead (&f_negative, sizeof (f_negative), bds, 1, &bufLoc,
973 0 : &numUsed);
974 : myAssert (numUsed == 0);
975 : myAssert ((mbit > 0) && (mbit < 32));
976 0 : memBitRead (&uli_temp, sizeof (sInt4), bds, mbit, &bufLoc, &numUsed);
977 0 : bds += numUsed;
978 0 : t_numBytes += numUsed;
979 0 : fstDiff = (f_negative) ? -1 * uli_temp : uli_temp;
980 : }
981 0 : memBitRead (&nbit, sizeof (nbit), bds, 5, &bufLoc, &numUsed);
982 0 : bds += numUsed;
983 0 : t_numBytes += numUsed;
984 0 : memBitRead (&f_negative, sizeof (f_negative), bds, 1, &bufLoc, &numUsed);
985 0 : bds += numUsed;
986 0 : t_numBytes += numUsed;
987 : myAssert ((nbit > 0) && (nbit < 32));
988 0 : memBitRead (&uli_temp, sizeof (sInt4), bds, nbit, &bufLoc, &numUsed);
989 0 : bds += numUsed;
990 0 : t_numBytes += numUsed;
991 0 : minVal = (f_negative) ? -1 * uli_temp : uli_temp;
992 0 : memBitRead (&LX, sizeof (LX), bds, 16, &bufLoc, &numUsed);
993 0 : bds += numUsed;
994 0 : t_numBytes += numUsed;
995 0 : grp = (TDLGroupType *) malloc (LX * sizeof (TDLGroupType));
996 0 : memBitRead (&ibit, sizeof (ibit), bds, 5, &bufLoc, &numUsed);
997 0 : bds += numUsed;
998 0 : t_numBytes += numUsed;
999 0 : memBitRead (&jbit, sizeof (jbit), bds, 5, &bufLoc, &numUsed);
1000 : /* Following assert is because it is the # of bits of # of bits. Which
1001 : * means that # of bits of value that has a max of 64. */
1002 : myAssert (jbit < 6);
1003 0 : bds += numUsed;
1004 0 : t_numBytes += numUsed;
1005 0 : memBitRead (&kbit, sizeof (kbit), bds, 5, &bufLoc, &numUsed);
1006 0 : bds += numUsed;
1007 0 : t_numBytes += numUsed;
1008 : myAssert (ibit < 33);
1009 0 : for (i = 0; i < LX; i++) {
1010 0 : if (ibit == 0) {
1011 0 : grp[i].min = 0;
1012 : } else {
1013 0 : memBitRead (&(grp[i].min), sizeof (sInt4), bds, ibit, &bufLoc,
1014 0 : &numUsed);
1015 0 : bds += numUsed;
1016 0 : t_numBytes += numUsed;
1017 : }
1018 : }
1019 : myAssert (jbit < 8);
1020 0 : for (i = 0; i < LX; i++) {
1021 0 : if (jbit == 0) {
1022 0 : grp[i].bit = 0;
1023 : } else {
1024 : myAssert (jbit <= sizeof (uChar) * 8);
1025 0 : memBitRead (&(grp[i].bit), sizeof (uChar), bds, jbit, &bufLoc,
1026 0 : &numUsed);
1027 0 : bds += numUsed;
1028 0 : t_numBytes += numUsed;
1029 : }
1030 : myAssert (grp[i].bit < 32);
1031 : }
1032 : myAssert (kbit < 33);
1033 0 : t_numPack = 0;
1034 0 : t_numBits = 0;
1035 0 : for (i = 0; i < LX; i++) {
1036 0 : if (kbit == 0) {
1037 0 : grp[i].num = 0;
1038 : } else {
1039 0 : memBitRead (&(grp[i].num), sizeof (sInt4), bds, kbit, &bufLoc,
1040 0 : &numUsed);
1041 0 : bds += numUsed;
1042 0 : t_numBytes += numUsed;
1043 : }
1044 0 : t_numPack += grp[i].num;
1045 0 : t_numBits += grp[i].num * grp[i].bit;
1046 : }
1047 0 : if (t_numPack != numPack) {
1048 : errSprintf ("Number packed %d != number of values in groups %d\n",
1049 0 : numPack, t_numPack);
1050 0 : free (grp);
1051 0 : return -1;
1052 : }
1053 0 : if ((t_numBytes + ceil (t_numBits / 8.)) > sectLen) {
1054 : errSprintf ("# bytes in groups %ld (%ld + %ld / 8) > sectLen %ld\n",
1055 : (sInt4) (t_numBytes + ceil (t_numBits / 8.)),
1056 0 : t_numBytes, t_numBits, sectLen);
1057 0 : free (grp);
1058 0 : return -1;
1059 : }
1060 0 : dataCnt = 0;
1061 0 : dataInd = 0;
1062 :
1063 : #ifdef DEBUG
1064 : printf ("nbit %d, ibit %d, jbit %d, kbit %d\n", nbit, ibit, jbit, kbit);
1065 : if ((t_numBytes + ceil (t_numBits / 8.)) != sectLen) {
1066 : printf ("Caution: # bytes in groups %d (%d + %d / 8) != "
1067 : "sectLen %d\n", (sInt4) (t_numBytes + ceil (t_numBits / 8.)),
1068 : t_numBytes, t_numBits, sectLen);
1069 : }
1070 : #endif
1071 :
1072 : /* Binary scale factor in TDLP has reverse sign from GRIB definition. */
1073 0 : scale = pow (10.0, -1 * DSF) * pow (2.0, -1 * BSF);
1074 :
1075 0 : meta->gridAttrib.f_maxmin = 0;
1076 : /* Work with Second order complex packed data. */
1077 0 : if (f_sndOrder) {
1078 : /* *INDENT-OFF* */
1079 : /* The algorithm appears to be:
1080 : * Data: a1 a2 a3 a4 a5 ...
1081 : * 1st diff: 0 b2 b3 b4 b5 ...
1082 : * 2nd diff: UK1 UK2 c3 c4 c5 ...
1083 : * We already know a1 and b2, and unpack a stream of UK1 UK2 c3 c4
1084 : * The problem is that UK1, UK2 is undefined. Originally I thought
1085 : * this was 0, or c3, but it appears that if b2 != 0, then
1086 : * UK2 = c3 + 2 b2, and UK1 = c3 + 1 * b2, otherwise it appears that
1087 : * UK1 == UK2, and typically UK1 == c3 (but not always). */
1088 : /* *INDENT-ON* */
1089 : myAssert (numPack >= 2);
1090 0 : if (f_secMiss) {
1091 0 : numVal = 0;
1092 0 : lastData = 0;
1093 0 : for (i = 0; i < LX; i++) {
1094 0 : maxVal = (1 << grp[i].bit) - 1;
1095 0 : for (j = 0; j < grp[i].num; j++) {
1096 : /* signed int. */
1097 0 : memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1098 0 : &bufLoc, &numUsed);
1099 0 : bds += numUsed;
1100 0 : if (li_temp == maxVal) {
1101 0 : data[dataInd] = meta->gridAttrib.missPri;
1102 0 : } else if (li_temp == (maxVal - 1)) {
1103 0 : data[dataInd] = meta->gridAttrib.missSec;
1104 : } else {
1105 0 : if (numVal > 1) {
1106 : #ifdef DEBUG
1107 : if (numVal == 2) {
1108 : if (fstDiff != 0) {
1109 : /*
1110 : myAssert (t_UK1 == li_temp + fstDiff);
1111 : myAssert (t_UK2 == li_temp + 2 * fstDiff);
1112 : */
1113 : } else {
1114 : myAssert (t_UK1 == t_UK2);
1115 : }
1116 : }
1117 : #endif
1118 0 : diff += (li_temp + grp[i].min + minVal);
1119 0 : data[dataInd] = data[lastData] + diff * scale;
1120 0 : lastData = dataInd;
1121 0 : } else if (numVal == 1) {
1122 0 : data[dataInd] = (origVal + fstDiff) * scale;
1123 0 : lastData = dataInd;
1124 0 : diff = fstDiff;
1125 : #ifdef DEBUG
1126 : t_UK2 = li_temp;
1127 : #endif
1128 : } else {
1129 0 : data[dataInd] = origVal * scale;
1130 : #ifdef DEBUG
1131 : t_UK1 = li_temp;
1132 : #endif
1133 : }
1134 0 : numVal++;
1135 0 : if (!meta->gridAttrib.f_maxmin) {
1136 0 : meta->gridAttrib.min = data[dataInd];
1137 0 : meta->gridAttrib.max = data[dataInd];
1138 0 : meta->gridAttrib.f_maxmin = 1;
1139 : } else {
1140 0 : if (data[dataInd] < meta->gridAttrib.min) {
1141 0 : meta->gridAttrib.min = data[dataInd];
1142 : }
1143 0 : if (data[dataInd] > meta->gridAttrib.max) {
1144 0 : meta->gridAttrib.max = data[dataInd];
1145 : }
1146 : }
1147 : }
1148 0 : dataCnt++;
1149 : dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1150 : (2 * (dataCnt / meta->gds.Nx) + 1) *
1151 0 : meta->gds.Nx - dataCnt - 1 : dataCnt;
1152 : myAssert ((dataInd < numPack) || (dataCnt == numPack));
1153 : }
1154 : }
1155 0 : } else if (f_primMiss) {
1156 0 : numVal = 0;
1157 0 : lastData = 0;
1158 0 : for (i = 0; i < LX; i++) {
1159 0 : maxVal = (1 << grp[i].bit) - 1;
1160 0 : for (j = 0; j < grp[i].num; j++) {
1161 : /* signed int. */
1162 0 : memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1163 0 : &bufLoc, &numUsed);
1164 0 : bds += numUsed;
1165 0 : f_missing = 0;
1166 0 : if (li_temp == maxVal) {
1167 0 : data[dataInd] = meta->gridAttrib.missPri;
1168 : /* In the case of grp[i].bit == 0, if grp[i].min == 0, then
1169 : * it is the missing value, otherwise regular value. Only
1170 : * need to be concerned for primary missing values. */
1171 0 : f_missing = 1;
1172 0 : if ((grp[i].bit == 0) && (grp[i].min != 0)) {
1173 : #ifdef DEBUG
1174 : printf ("This doesn't happen often.\n");
1175 : printf ("%d %d %d\n", (int) i, grp[i].bit, grp[i].min);
1176 : #endif
1177 : myAssert (1 == 2);
1178 0 : f_missing = 0;
1179 : }
1180 : }
1181 0 : if (!f_missing) {
1182 0 : if (numVal > 1) {
1183 : #ifdef DEBUG
1184 : if (numVal == 2) {
1185 : if (fstDiff != 0) {
1186 : /*
1187 : myAssert (t_UK1 == li_temp + fstDiff);
1188 : myAssert (t_UK2 == li_temp + 2 * fstDiff);
1189 : */
1190 : } else {
1191 : myAssert (t_UK1 == t_UK2);
1192 : }
1193 : }
1194 : #endif
1195 0 : diff += (li_temp + grp[i].min + minVal);
1196 0 : data[dataInd] = data[lastData] + diff * scale;
1197 0 : lastData = dataInd;
1198 0 : } else if (numVal == 1) {
1199 0 : data[dataInd] = (origVal + fstDiff) * scale;
1200 0 : lastData = dataInd;
1201 0 : diff = fstDiff;
1202 : #ifdef DEBUG
1203 : t_UK2 = li_temp;
1204 : #endif
1205 : } else {
1206 0 : data[dataInd] = origVal * scale;
1207 : #ifdef DEBUG
1208 : t_UK1 = li_temp;
1209 : #endif
1210 : }
1211 0 : numVal++;
1212 0 : if (!meta->gridAttrib.f_maxmin) {
1213 0 : meta->gridAttrib.min = data[dataInd];
1214 0 : meta->gridAttrib.max = data[dataInd];
1215 0 : meta->gridAttrib.f_maxmin = 1;
1216 : } else {
1217 0 : if (data[dataInd] < meta->gridAttrib.min) {
1218 0 : meta->gridAttrib.min = data[dataInd];
1219 : }
1220 0 : if (data[dataInd] > meta->gridAttrib.max) {
1221 0 : meta->gridAttrib.max = data[dataInd];
1222 : }
1223 : }
1224 : }
1225 0 : dataCnt++;
1226 : dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1227 : (2 * (dataCnt / meta->gds.Nx) + 1) *
1228 0 : meta->gds.Nx - dataCnt - 1 : dataCnt;
1229 : myAssert ((dataInd < numPack) || (dataCnt == numPack));
1230 : }
1231 : }
1232 : } else {
1233 0 : lastData = 0;
1234 0 : for (i = 0; i < LX; i++) {
1235 0 : for (j = 0; j < grp[i].num; j++) {
1236 : /* signed int. */
1237 0 : memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1238 0 : &bufLoc, &numUsed);
1239 0 : bds += numUsed;
1240 0 : if (dataCnt > 1) {
1241 : #ifdef DEBUG
1242 : if (dataCnt == 2) {
1243 : if (fstDiff != 0) {
1244 : /*
1245 : myAssert (t_UK1 == li_temp + fstDiff);
1246 : myAssert (t_UK2 == li_temp + 2 * fstDiff);
1247 : */
1248 : } else {
1249 : myAssert (t_UK1 == t_UK2);
1250 : }
1251 : }
1252 : #endif
1253 0 : diff += (li_temp + grp[i].min + minVal);
1254 0 : data[dataInd] = data[lastData] + diff * scale;
1255 0 : lastData = dataInd;
1256 0 : } else if (dataCnt == 1) {
1257 0 : data[dataInd] = (origVal + fstDiff) * scale;
1258 0 : lastData = dataInd;
1259 0 : diff = fstDiff;
1260 : #ifdef DEBUG
1261 : t_UK2 = li_temp;
1262 : #endif
1263 : } else {
1264 0 : data[dataInd] = origVal * scale;
1265 : #ifdef DEBUG
1266 : t_UK1 = li_temp;
1267 : #endif
1268 : }
1269 : #ifdef DEBUG
1270 : /*
1271 : if (i >= 4153) {
1272 : */
1273 : /*
1274 : if ((data[dataInd] > 100) || (data[dataInd] < -100)) {
1275 : */
1276 : /*
1277 : if ((diff > 50) || (diff < -50)) {
1278 : printf ("li_temp :: %ld, diff = %ld\n", li_temp, diff);
1279 : printf ("data[dataInd] :: %f\n", data[dataInd]);
1280 : printf ("Group # %d element %d, grp[i].min %ld, "
1281 : "grp[i].bit %d, minVal %ld\n", i, j, grp[i].min,
1282 : grp[i].bit, minVal);
1283 : }
1284 : */
1285 : #endif
1286 0 : if (!meta->gridAttrib.f_maxmin) {
1287 0 : meta->gridAttrib.min = data[dataInd];
1288 0 : meta->gridAttrib.max = data[dataInd];
1289 0 : meta->gridAttrib.f_maxmin = 1;
1290 : } else {
1291 0 : if (data[dataInd] < meta->gridAttrib.min) {
1292 0 : meta->gridAttrib.min = data[dataInd];
1293 : }
1294 0 : if (data[dataInd] > meta->gridAttrib.max) {
1295 0 : meta->gridAttrib.max = data[dataInd];
1296 : }
1297 : }
1298 0 : dataCnt++;
1299 : dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1300 : (2 * (dataCnt / meta->gds.Nx) + 1) *
1301 0 : meta->gds.Nx - dataCnt - 1 : dataCnt;
1302 : myAssert ((dataInd < numPack) || (dataCnt == numPack));
1303 : }
1304 : }
1305 0 : numVal = dataCnt;
1306 : }
1307 :
1308 : /* Work with regular complex packed data. */
1309 : } else {
1310 : #ifdef DEBUG
1311 : /*
1312 : printf ("Work with regular complex packed data\n");
1313 : */
1314 : #endif
1315 0 : if (f_secMiss) {
1316 0 : numVal = 0;
1317 0 : for (i = 0; i < LX; i++) {
1318 0 : maxVal = (1 << grp[i].bit) - 1;
1319 0 : for (j = 0; j < grp[i].num; j++) {
1320 : /* signed int. */
1321 0 : memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1322 0 : &bufLoc, &numUsed);
1323 0 : bds += numUsed;
1324 0 : if (li_temp == maxVal) {
1325 0 : data[dataInd] = meta->gridAttrib.missPri;
1326 0 : } else if (li_temp == (maxVal - 1)) {
1327 0 : data[dataInd] = meta->gridAttrib.missSec;
1328 : } else {
1329 0 : data[dataInd] = (li_temp + grp[i].min + minVal) * scale;
1330 0 : numVal++;
1331 0 : if (!meta->gridAttrib.f_maxmin) {
1332 0 : meta->gridAttrib.min = data[dataInd];
1333 0 : meta->gridAttrib.max = data[dataInd];
1334 0 : meta->gridAttrib.f_maxmin = 1;
1335 : } else {
1336 0 : if (data[dataInd] < meta->gridAttrib.min) {
1337 0 : meta->gridAttrib.min = data[dataInd];
1338 : }
1339 0 : if (data[dataInd] > meta->gridAttrib.max) {
1340 0 : meta->gridAttrib.max = data[dataInd];
1341 : }
1342 : }
1343 : }
1344 0 : dataCnt++;
1345 : dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1346 : (2 * (dataCnt / meta->gds.Nx) + 1) *
1347 0 : meta->gds.Nx - dataCnt - 1 : dataCnt;
1348 : myAssert ((dataInd < numPack) || (dataCnt == numPack));
1349 : }
1350 : }
1351 0 : } else if (f_primMiss) {
1352 : #ifdef DEBUG
1353 : /*
1354 : printf ("Work with primary missing data\n");
1355 : */
1356 : #endif
1357 0 : numVal = 0;
1358 0 : for (i = 0; i < LX; i++) {
1359 0 : maxVal = (1 << grp[i].bit) - 1;
1360 0 : for (j = 0; j < grp[i].num; j++) {
1361 0 : memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1362 0 : &bufLoc, &numUsed);
1363 0 : bds += numUsed;
1364 0 : f_missing = 0;
1365 0 : if (li_temp == maxVal) {
1366 0 : data[dataInd] = meta->gridAttrib.missPri;
1367 : /* In the case of grp[i].bit == 0, if grp[i].min == 0, then
1368 : * it is the missing value, otherwise regular value. Only
1369 : * need to be concerned for primary missing values. */
1370 0 : f_missing = 1;
1371 0 : if ((grp[i].bit == 0) && (grp[i].min != 0)) {
1372 : #ifdef DEBUG
1373 : printf ("This doesn't happen often.\n");
1374 : printf ("%d %d %d\n", (int) i, grp[i].bit, grp[i].min);
1375 : myAssert (1 == 2);
1376 : #endif
1377 0 : f_missing = 0;
1378 : }
1379 : }
1380 0 : if (!f_missing) {
1381 0 : data[dataInd] = (li_temp + grp[i].min + minVal) * scale;
1382 0 : numVal++;
1383 0 : if (!meta->gridAttrib.f_maxmin) {
1384 0 : meta->gridAttrib.min = data[dataInd];
1385 0 : meta->gridAttrib.max = data[dataInd];
1386 0 : meta->gridAttrib.f_maxmin = 1;
1387 : } else {
1388 0 : if (data[dataInd] < meta->gridAttrib.min) {
1389 0 : meta->gridAttrib.min = data[dataInd];
1390 : }
1391 0 : if (data[dataInd] > meta->gridAttrib.max) {
1392 0 : meta->gridAttrib.max = data[dataInd];
1393 : }
1394 : }
1395 : }
1396 0 : dataCnt++;
1397 : dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1398 : (2 * (dataCnt / meta->gds.Nx) + 1) *
1399 0 : meta->gds.Nx - dataCnt - 1 : dataCnt;
1400 : myAssert ((dataInd < numPack) || (dataCnt == numPack));
1401 : }
1402 : }
1403 : } else {
1404 0 : for (i = 0; i < LX; i++) {
1405 0 : for (j = 0; j < grp[i].num; j++) {
1406 0 : memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1407 0 : &bufLoc, &numUsed);
1408 0 : bds += numUsed;
1409 0 : data[dataInd] = (li_temp + grp[i].min + minVal) * scale;
1410 0 : if (!meta->gridAttrib.f_maxmin) {
1411 0 : meta->gridAttrib.min = data[dataInd];
1412 0 : meta->gridAttrib.max = data[dataInd];
1413 0 : meta->gridAttrib.f_maxmin = 1;
1414 : } else {
1415 0 : if (data[dataInd] < meta->gridAttrib.min) {
1416 0 : meta->gridAttrib.min = data[dataInd];
1417 : }
1418 0 : if (data[dataInd] > meta->gridAttrib.max) {
1419 0 : meta->gridAttrib.max = data[dataInd];
1420 : }
1421 : }
1422 0 : dataCnt++;
1423 : dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1424 : (2 * (dataCnt / meta->gds.Nx) + 1) *
1425 0 : meta->gds.Nx - dataCnt - 1 : dataCnt;
1426 : myAssert ((dataInd < numPack) || (dataCnt == numPack));
1427 : }
1428 : }
1429 0 : numVal = dataCnt;
1430 : }
1431 : }
1432 0 : meta->gridAttrib.numMiss = dataCnt - numVal;
1433 0 : meta->gridAttrib.refVal = minVal * scale;
1434 :
1435 0 : free (grp);
1436 0 : return 0;
1437 : }
1438 :
1439 : /*****************************************************************************
1440 : * ReadTDLPRecord() --
1441 : *
1442 : * Arthur Taylor / MDL
1443 : *
1444 : * PURPOSE
1445 : * Reads in a TDLP message, and parses the data into various data
1446 : * structures, for use with other code.
1447 : *
1448 : * ARGUMENTS
1449 : * fp = An opened TDLP file already at the correct message. (Input)
1450 : * TDLP_Data = The read in TDLP data. (Output)
1451 : * tdlp_DataLen = Size of TDLP_Data. (Output)
1452 : * meta = A filled in meta structure (Output)
1453 : * IS = The structure containing all the arrays that the
1454 : * unpacker uses (Output)
1455 : * sect0 = Already read in section 0 data. (Input)
1456 : * tdlpLen = Length of the TDLP message. (Input)
1457 : * majEarth = Used to override the TDLP major axis of earth. (Input)
1458 : * minEarth = Used to override the TDLP minor axis of earth. (Input)
1459 : *
1460 : * FILES/DATABASES:
1461 : * An already opened file pointing to the desired TDLP message.
1462 : *
1463 : * RETURNS: int (could use errSprintf())
1464 : * 0 = OK
1465 : * -1 = Problems reading in the PDS.
1466 : * -2 = Problems reading in the GDS.
1467 : * -3 = Problems reading in the BMS.
1468 : * -4 = Problems reading in the BDS.
1469 : * -5 = Problems reading the closing section.
1470 : *
1471 : * HISTORY
1472 : * 10/2004 Arthur Taylor (MDL): Created
1473 : *
1474 : * NOTES
1475 : *****************************************************************************
1476 : */
1477 0 : int ReadTDLPRecord (DataSource &fp, double **TDLP_Data, uInt4 *tdlp_DataLen,
1478 : grib_MetaData *meta, IS_dataType *IS,
1479 : sInt4 sect0[SECT0LEN_WORD], uInt4 tdlpLen,
1480 : double majEarth, double minEarth)
1481 : {
1482 : sInt4 nd5; /* Size of TDLP message rounded up to the nearest *
1483 : * sInt4. */
1484 : uChar *c_ipack; /* A char ptr to the message stored in IS->ipack */
1485 : sInt4 curLoc; /* Current location in the GRIB message. */
1486 : char f_gds; /* flag if there is a GDS section. */
1487 : char f_bms; /* flag if there is a BMS section. */
1488 : short int DSF; /* Decimal Scale Factor for unpacking the data. */
1489 : short int BSF; /* Binary Scale Factor for unpacking the data. */
1490 : double *tdlp_Data; /* A pointer to TDLP_Data for ease of manipulation. */
1491 0 : double unitM = 1; /* M in y = Mx + B, for unit conversion. */
1492 0 : double unitB = 0; /* B in y = Mx + B, for unit conversion. */
1493 : sInt4 li_temp; /* Used to make sure section 5 is 7777. */
1494 : size_t pad; /* Number of bytes to pad the message to get to the
1495 : * correct byte boundary. */
1496 : char buffer[24]; /* Read the trailing bytes in the TDLPack record. */
1497 : uChar *bitmap; /* Would contain bitmap data if it was supported. */
1498 :
1499 : /* Make room for entire message, and read it in. */
1500 : /* nd5 needs to be tdlpLen in (sInt4) units rounded up. */
1501 0 : nd5 = (tdlpLen + 3) / 4;
1502 0 : if (nd5 > IS->ipackLen) {
1503 0 : IS->ipackLen = nd5;
1504 : IS->ipack = (sInt4 *) realloc ((void *) (IS->ipack),
1505 0 : (IS->ipackLen) * sizeof (sInt4));
1506 : }
1507 0 : c_ipack = (uChar *) IS->ipack;
1508 : /* Init last sInt4 to 0, to make sure that the padded bytes are 0. */
1509 0 : IS->ipack[nd5 - 1] = 0;
1510 : /* Init first 2 sInt4 to sect0. */
1511 0 : memcpy (c_ipack, sect0, SECT0LEN_WORD * 2);
1512 : /* Read in the rest of the message. */
1513 0 : if (fp.DataSourceFread(c_ipack + SECT0LEN_WORD * 2, sizeof (char),
1514 0 : (tdlpLen - SECT0LEN_WORD * 2)) + SECT0LEN_WORD * 2 != tdlpLen) {
1515 0 : errSprintf ("Ran out of file\n");
1516 0 : return -1;
1517 : }
1518 :
1519 : /* Preceeding was in degrib2, next part is specific to TDLP. */
1520 0 : curLoc = 8;
1521 0 : if (ReadTDLPSect1 (c_ipack + curLoc, tdlpLen, &curLoc, &(meta->pdsTdlp),
1522 : &f_gds, &f_bms, &DSF, &BSF) != 0) {
1523 0 : preErrSprintf ("Inside ReadGrib1Record\n");
1524 0 : return -1;
1525 : }
1526 :
1527 : /* Figure out some basic stuff about the grid. */
1528 0 : free (meta->element);
1529 0 : meta->element = NULL;
1530 0 : free (meta->unitName);
1531 0 : meta->unitName = NULL;
1532 0 : free (meta->comment);
1533 0 : meta->comment = NULL;
1534 0 : free (meta->shortFstLevel);
1535 0 : meta->shortFstLevel = NULL;
1536 0 : free (meta->longFstLevel);
1537 0 : meta->longFstLevel = NULL;
1538 : TDLP_ElemSurfUnit (&(meta->pdsTdlp), &(meta->element), &(meta->unitName),
1539 : &(meta->comment), &(meta->shortFstLevel),
1540 0 : &(meta->longFstLevel));
1541 0 : meta->center = 7; /* US NWS, NCEP */
1542 0 : meta->subcenter = 14; /* NWS Meteorological Development Laboratory */
1543 :
1544 : /* strftime (meta->refTime, 20, "%Y%m%d%H%M",
1545 : gmtime (&(meta->pdsTdlp.refTime)));
1546 : */
1547 0 : Clock_Print (meta->refTime, 20, meta->pdsTdlp.refTime, "%Y%m%d%H%M", 0);
1548 :
1549 : /*
1550 : validTime = meta->pdsTdlp.refTime + meta->pdsTdlp.project;
1551 : strftime (meta->validTime, 20, "%Y%m%d%H%M", gmtime (&(validTime)));
1552 : */
1553 : Clock_Print (meta->validTime, 20, meta->pdsTdlp.refTime +
1554 0 : meta->pdsTdlp.project, "%Y%m%d%H%M", 0);
1555 :
1556 0 : meta->deltTime = meta->pdsTdlp.project;
1557 :
1558 : /* Get the Grid Definition Section. */
1559 0 : if (f_gds) {
1560 0 : if (ReadTDLPSect2 (c_ipack + curLoc, tdlpLen, &curLoc,
1561 : &(meta->gds)) != 0) {
1562 0 : preErrSprintf ("Inside ReadGrib1Record\n");
1563 0 : return -2;
1564 : }
1565 : } else {
1566 0 : errSprintf ("Don't know how to handle vector data yet.\n");
1567 0 : return -2;
1568 : }
1569 :
1570 : /* Allow over ride of the earth radii. */
1571 0 : if ((majEarth > 6300) && (majEarth < 6400)) {
1572 0 : if ((minEarth > 6300) && (minEarth < 6400)) {
1573 0 : meta->gds.f_sphere = 0;
1574 0 : meta->gds.majEarth = majEarth;
1575 0 : meta->gds.minEarth = minEarth;
1576 0 : if (majEarth == minEarth) {
1577 0 : meta->gds.f_sphere = 1;
1578 : }
1579 : } else {
1580 0 : meta->gds.f_sphere = 1;
1581 0 : meta->gds.majEarth = majEarth;
1582 0 : meta->gds.minEarth = majEarth;
1583 : }
1584 : }
1585 :
1586 : /* Allocate memory for the grid. */
1587 0 : if (meta->gds.numPts > *tdlp_DataLen) {
1588 0 : *tdlp_DataLen = meta->gds.numPts;
1589 : *TDLP_Data = (double *) realloc ((void *) (*TDLP_Data),
1590 0 : (*tdlp_DataLen) * sizeof (double));
1591 : }
1592 0 : tdlp_Data = *TDLP_Data;
1593 :
1594 : /* Get the Bit Map Section. */
1595 0 : if (f_bms) {
1596 : /* errSprintf ("Bitmap data is Not Supported\n");*/
1597 : /* Need to allocate bitmap when this is implemented. */
1598 : ReadTDLPSect3 (c_ipack + curLoc, tdlpLen, &curLoc, bitmap,
1599 0 : meta->gds.numPts);
1600 0 : return -1;
1601 : }
1602 :
1603 : /* Read the GRID. */
1604 0 : if (ReadTDLPSect4 (c_ipack + curLoc, tdlpLen, &curLoc, DSF, BSF,
1605 : tdlp_Data, meta, unitM, unitB) != 0) {
1606 0 : preErrSprintf ("Inside ReadTDLPRecord\n");
1607 0 : return -4;
1608 : }
1609 :
1610 : /* Read section 5. If it is "7777" == 926365495 we are done. */
1611 0 : memcpy (&li_temp, c_ipack + curLoc, 4);
1612 0 : if (li_temp != 926365495L) {
1613 0 : errSprintf ("Did not find the end of the message.\n");
1614 0 : return -5;
1615 : }
1616 0 : curLoc += 4;
1617 : /* Read the trailing part of the message. */
1618 : /* TDLPack uses 4 bytes for FORTRAN record size, then another 8 bytes for
1619 : * the size of the record (so FORTRAN can see it), then the data rounded
1620 : * up to an 8 byte boundary, then a trailing 4 bytes for a final FORTRAN
1621 : * record size. However it only stores in the gribLen the non-rounded
1622 : * amount, so we need to take care of the rounding, and the trailing 4
1623 : * bytes here. */
1624 0 : pad = ((sInt4) (ceil (tdlpLen / 8.0))) * 8 - tdlpLen + 4;
1625 0 : if (fp.DataSourceFread(buffer, sizeof (char), pad) != pad) {
1626 0 : errSprintf ("Ran out of file\n");
1627 0 : return -1;
1628 : }
1629 0 : return 0;
1630 : }
1631 :
1632 : /*****************************************************************************
1633 : * TDL_ScaleData() --
1634 : *
1635 : * Arthur Taylor / MDL
1636 : *
1637 : * PURPOSE
1638 : * Deal with scaling while excluding the primary and secondary missing
1639 : * values. After this, dst should contain scaled data + primary or secondary
1640 : * missing values
1641 : * "tdlpack library"::pack2d.f line 257 or search for:
1642 : "the above statement"
1643 : *
1644 : * ARGUMENTS
1645 : * Src = The original data. (Input)
1646 : * Dst = The scaled data. (Output)
1647 : * numData = The number of elements in data. (Input)
1648 : * DSF = Decimal Scale Factor for scaling the data. (Input)
1649 : * BSF = Binary Scale Factor for scaling the data. (Input)
1650 : * f_primMiss = Flag saying if we have a primary missing value (In/Out)
1651 : * primMiss = primary missing value. (In/Out)
1652 : * f_secMiss = Flag saying if we have a secondary missing value (In/Out)
1653 : * secMiss = secondary missing value. (In/Out)
1654 : * f_min = Flag saying if we have the minimum value. (Output)
1655 : * min = minimum scaled value in the grid. (Output)
1656 : *
1657 : * FILES/DATABASES: None
1658 : *
1659 : * RETURNS: void
1660 : *
1661 : * HISTORY
1662 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
1663 : *
1664 : * NOTES
1665 : *****************************************************************************
1666 : */
1667 : #define SCALE_MISSING 10000L
1668 0 : static void TDL_ScaleData (double *Src, sInt4 *Dst, sInt4 numData,
1669 : int DSF, int BSF, char *f_primMiss,
1670 : double *primMiss, char *f_secMiss,
1671 : double *secMiss, char *f_min, sInt4 *min)
1672 : {
1673 : sInt4 cnt;
1674 0 : double *src = Src;
1675 0 : sInt4 *dst = Dst;
1676 0 : double scale = pow (10.0, -1 * DSF) * pow (2.0, -1 * BSF);
1677 0 : char f_actualPrim = 0;
1678 0 : char f_actualSec = 0;
1679 0 : sInt4 li_primMiss = (sInt4) (*primMiss * SCALE_MISSING + .5);
1680 0 : sInt4 li_secMiss = (sInt4) (*secMiss * SCALE_MISSING + .5);
1681 :
1682 0 : *f_min = 0;
1683 0 : for (cnt = 0; cnt < numData; cnt++) {
1684 0 : if (((*f_primMiss) || (*f_secMiss)) && (*src == *primMiss)) {
1685 0 : *(dst++) = li_primMiss;
1686 0 : src++;
1687 0 : f_actualPrim = 1;
1688 0 : } else if ((*f_secMiss) && (*src == *secMiss)) {
1689 0 : *(dst++) = li_secMiss;
1690 0 : src++;
1691 0 : f_actualSec = 1;
1692 : } else {
1693 0 : *(dst) = (long int) (floor ((*(src++) / scale) + .5));
1694 : /* Check if scaled value == primary missing value. */
1695 0 : if (((*f_primMiss) || (*f_secMiss)) && (*dst == li_primMiss)) {
1696 0 : *dst = *dst - 1;
1697 : }
1698 : /* Check if scaled value == secondary missing value. */
1699 0 : if ((*f_secMiss) && (*dst == li_secMiss)) {
1700 0 : *dst = *dst - 1;
1701 : /* Check if adjustment caused scaled value == primary missing. */
1702 0 : if (*dst == li_primMiss) {
1703 0 : *dst = *dst - 1;
1704 : }
1705 : }
1706 0 : if (!(*f_min)) {
1707 0 : *min = *dst;
1708 0 : *f_min = 1;
1709 0 : } else if (*min > *dst) {
1710 0 : *min = *dst;
1711 : }
1712 0 : dst++;
1713 : }
1714 : }
1715 0 : if ((*f_secMiss) && (!f_actualSec)) {
1716 0 : *f_secMiss = 0;
1717 : }
1718 0 : if (((*f_secMiss) || (*f_primMiss)) && (!f_actualPrim)) {
1719 0 : *f_primMiss = 0;
1720 : /* Check consistency. */
1721 0 : if (*f_secMiss) {
1722 0 : *f_secMiss = 0;
1723 0 : *f_primMiss = 1;
1724 0 : *primMiss = *secMiss;
1725 : }
1726 : }
1727 0 : }
1728 :
1729 : /*****************************************************************************
1730 : * TDL_ReorderGrid() --
1731 : *
1732 : * Arthur Taylor / MDL
1733 : *
1734 : * PURPOSE
1735 : * Loop through the data, so that
1736 : * data is: "a1,1 ... a1,n a2,n ... a2,1 ..."
1737 : * instead of: "a1,1 ... a1,n a2,1 ... a2,n ..."
1738 : *
1739 : * ARGUMENTS
1740 : * Src = The data. (Input/Output)
1741 : * NX = The number of X values. (Input)
1742 : * NY = The number of Y values. (Input)
1743 : *
1744 : * FILES/DATABASES: None
1745 : *
1746 : * RETURNS: void
1747 : *
1748 : * HISTORY
1749 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
1750 : *
1751 : * NOTES
1752 : *****************************************************************************
1753 : */
1754 0 : static void TDL_ReorderGrid (sInt4 *Src, short int NX, short int NY)
1755 : {
1756 : int i, j;
1757 : sInt4 *src1, *src2;
1758 : sInt4 li_temp;
1759 :
1760 0 : for (j = 1; j < NY; j += 2) {
1761 0 : src1 = Src + j * NX;
1762 0 : src2 = Src + (j + 1) * NX - 1;
1763 0 : for (i = 0; i < (NX / 2); i++) {
1764 0 : li_temp = *src1;
1765 0 : *(src1++) = *src2;
1766 0 : *(src2--) = li_temp;
1767 : }
1768 : }
1769 0 : }
1770 :
1771 : /*****************************************************************************
1772 : * TDL_GetSecDiff() --
1773 : *
1774 : * Arthur Taylor / MDL
1775 : *
1776 : * PURPOSE
1777 : * Get the second order difference where we have special values for missing,
1778 : * and for actual data we have the following scheme.
1779 : * Data: a1 a2 a3 a4 a5 ...
1780 : * 1st diff: 0 b2 b3 b4 b5 ...
1781 : * 2nd diff: UK1 UK2 c3 c4 c5 ...
1782 : * where UK1 = c3 + b2, and UK2 = c3 + 2 * b2. Note: The choice of UK1, and
1783 : * UK2 doesn't matter because of the following FORTRAN unpacking code:
1784 : * IWORK(1)=IFIRST
1785 : * IWORK(2)=IWORK(1)+IFOD
1786 : * ISUM=IFOD
1787 : * DO 385 K=3,IS4(3)
1788 : * ISUM=IWORK(K)+ISUM
1789 : * IWORK(K)=IWORK(K-1)+ISUM
1790 : * 385 CONTINUE
1791 : * So ISUM is a function of IWORK(3), not IWORK(1).
1792 : *
1793 : * ARGUMENTS
1794 : * Data = The data. (Input)
1795 : * numData = The number of elements in data. (Input)
1796 : * SecDiff = The secondary differences of the data. (Output)
1797 : * f_primMiss = Flag saying if we have a primary missing value (Input)
1798 : * li_primMiss = Scaled primary missing value. (Input)
1799 : * a1 = First non-missing value in the field. (Output)
1800 : * b2 = First non-missing value in the 1st order delta field (Out)
1801 : * min = The minimum value (Input).
1802 : *
1803 : * FILES/DATABASES: None
1804 : *
1805 : * RETURNS: int
1806 : * 0 = Success
1807 : * 1 = Couldn't find second differences (don't use).
1808 : *
1809 : * HISTORY
1810 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
1811 : *
1812 : * NOTES
1813 : *****************************************************************************
1814 : */
1815 0 : static int TDL_GetSecDiff (sInt4 *Data, int numData, sInt4 *SecDiff,
1816 : char f_primMiss, sInt4 li_primMiss,
1817 : sInt4 *a1, sInt4 *b2, sInt4 *min)
1818 : {
1819 : int i;
1820 0 : char f_min = 0;
1821 0 : sInt4 last = 0, before_last = 0;
1822 0 : int a1Index = -1;
1823 0 : int a2Index = -1;
1824 :
1825 0 : if (numData < 3) {
1826 0 : return 1;
1827 : }
1828 0 : if (f_primMiss) {
1829 0 : for (i = 0; i < numData; i++) {
1830 0 : if (Data[i] == li_primMiss) {
1831 0 : SecDiff[i] = li_primMiss;
1832 0 : } else if (a1Index == -1) {
1833 0 : a1Index = i;
1834 0 : *a1 = Data[a1Index];
1835 0 : } else if (a2Index == -1) {
1836 0 : a2Index = i;
1837 0 : *b2 = Data[a2Index] - Data[a1Index];
1838 0 : before_last = Data[a1Index];
1839 0 : last = Data[a2Index];
1840 : } else {
1841 0 : SecDiff[i] = Data[i] - 2 * last + before_last;
1842 0 : before_last = last;
1843 0 : last = Data[i];
1844 0 : if (!f_min) {
1845 : /* Set the UK1, UK2 values. */
1846 0 : *min = SecDiff[i];
1847 0 : f_min = 1;
1848 0 : SecDiff[a1Index] = SecDiff[i] + *b2;
1849 0 : SecDiff[a2Index] = SecDiff[i] + 2 * (*b2);
1850 0 : } else if (*min > SecDiff[i]) {
1851 0 : *min = SecDiff[i];
1852 : }
1853 : }
1854 : }
1855 0 : if (!f_min) {
1856 0 : return 1;
1857 : }
1858 : } else {
1859 0 : *a1 = Data[0];
1860 0 : *b2 = Data[1] - Data[0];
1861 0 : for (i = 3; i < numData; i++) {
1862 0 : SecDiff[i] = Data[i] - 2 * Data[i - 1] - Data[i - 2];
1863 0 : if (i == 3) {
1864 0 : *min = SecDiff[i];
1865 : /* Set the UK1, UK2 values. */
1866 0 : SecDiff[0] = SecDiff[i] + *b2;
1867 0 : SecDiff[1] = SecDiff[i] + 2 * (*b2);
1868 0 : } else if (*min > SecDiff[i]) {
1869 0 : *min = SecDiff[i];
1870 : }
1871 : }
1872 : }
1873 0 : return 0;
1874 : }
1875 :
1876 : /*****************************************************************************
1877 : * TDL_UseSecDiff_Prim() --
1878 : *
1879 : * Arthur Taylor / MDL
1880 : *
1881 : * PURPOSE
1882 : * Checks if the average range of 2nd order differences < average range of
1883 : * 0 order differnces, to determine if we should use second order differences
1884 : * This deals with the case when we have primary missing values.
1885 : *
1886 : * ARGUMENTS
1887 : * Data = The data. (Input)
1888 : * numData = The number of elements in data. (Input)
1889 : * SecDiff = The secondary differences of the data. (Input)
1890 : * li_primMiss = Scaled primary missing value. (Input)
1891 : * minGroup = The minimum group size. (Input)
1892 : *
1893 : * FILES/DATABASES: None
1894 : *
1895 : * RETURNS: int
1896 : * 0 = Don't use 2nd order differences.
1897 : * 1 = Use 2nd order differences.
1898 : *
1899 : * HISTORY
1900 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
1901 : *
1902 : * NOTES
1903 : *****************************************************************************
1904 : */
1905 0 : static int TDL_UseSecDiff_Prim (sInt4 *Data, sInt4 numData,
1906 : sInt4 *SecDiff, sInt4 li_primMiss,
1907 : int minGroup)
1908 : {
1909 : int i, locCnt;
1910 : int range0, range2;
1911 : int tot0, tot2;
1912 : char f_min;
1913 0 : sInt4 min = 0, max = 0;
1914 :
1915 0 : locCnt = 0;
1916 0 : range0 = 0;
1917 0 : tot0 = 0;
1918 0 : f_min = 0;
1919 : /* Compute scores for no differences */
1920 0 : for (i = 0; i < numData; i++) {
1921 0 : if (Data[i] != li_primMiss) {
1922 0 : if (!f_min) {
1923 0 : min = Data[i];
1924 0 : max = Data[i];
1925 0 : f_min = 1;
1926 : } else {
1927 0 : if (min > Data[i])
1928 0 : min = Data[i];
1929 0 : if (max < Data[i])
1930 0 : max = Data[i];
1931 : }
1932 : }
1933 0 : locCnt++;
1934 : /* Fake a "group" by using the minimum group size. */
1935 0 : if (locCnt == minGroup) {
1936 0 : if (f_min) {
1937 0 : range0 += (max - min);
1938 0 : tot0++;
1939 0 : f_min = 0;
1940 : }
1941 0 : locCnt = 0;
1942 : }
1943 : }
1944 0 : if (locCnt != 0) {
1945 0 : range0 += (max - min);
1946 0 : tot0++;
1947 : }
1948 :
1949 0 : locCnt = 0;
1950 0 : range2 = 0;
1951 0 : tot2 = 0;
1952 0 : f_min = 0;
1953 : /* Compute scores for second order differences */
1954 0 : for (i = 0; i < numData; i++) {
1955 0 : if (SecDiff[i] != li_primMiss) {
1956 0 : if (!f_min) {
1957 0 : min = SecDiff[i];
1958 0 : max = SecDiff[i];
1959 0 : f_min = 1;
1960 : } else {
1961 0 : if (min > SecDiff[i])
1962 0 : min = SecDiff[i];
1963 0 : if (max < SecDiff[i])
1964 0 : max = SecDiff[i];
1965 : }
1966 : }
1967 0 : locCnt++;
1968 : /* Fake a "group" by using the minimum group size. */
1969 0 : if (locCnt == minGroup) {
1970 0 : if (f_min) {
1971 0 : range2 += (max - min);
1972 0 : tot2++;
1973 0 : f_min = 0;
1974 : }
1975 0 : locCnt = 0;
1976 : }
1977 : }
1978 0 : if (locCnt != 0) {
1979 0 : range2 += (max - min);
1980 0 : tot2++;
1981 : }
1982 :
1983 : /* Compare average group size of no differencing to second order. */
1984 0 : if ((range0 / (tot0 + 0.0)) <= (range2 / (tot2 + 0.0))) {
1985 0 : return 0;
1986 : } else {
1987 0 : return 1;
1988 : }
1989 : }
1990 :
1991 : /*****************************************************************************
1992 : * TDL_UseSecDiff() --
1993 : *
1994 : * Arthur Taylor / MDL
1995 : *
1996 : * PURPOSE
1997 : * Checks if the average range of 2nd order differences < average range of
1998 : * 0 order differnces, to determine if we should use second order differences
1999 : *
2000 : * ARGUMENTS
2001 : * Data = The data. (Input)
2002 : * numData = The number of elements in data. (Input)
2003 : * SecDiff = The secondary differences of the data. (Input)
2004 : * minGroup = The minimum group size. (Input)
2005 : *
2006 : * FILES/DATABASES: None
2007 : *
2008 : * RETURNS: int
2009 : * 0 = Don't use 2nd order differences.
2010 : * 1 = Use 2nd order differences.
2011 : *
2012 : * HISTORY
2013 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
2014 : *
2015 : * NOTES
2016 : *****************************************************************************
2017 : */
2018 0 : static int TDL_UseSecDiff (sInt4 *Data, sInt4 numData,
2019 : sInt4 *SecDiff, int minGroup)
2020 : {
2021 : int i, locCnt;
2022 : int range0, range2;
2023 : int tot0, tot2;
2024 0 : sInt4 min = 0, max = 0;
2025 :
2026 0 : locCnt = 0;
2027 0 : range0 = 0;
2028 0 : tot0 = 0;
2029 : /* Compute scores for no differences */
2030 0 : for (i = 0; i < numData; i++) {
2031 0 : if (locCnt == 0) {
2032 0 : min = Data[i];
2033 0 : max = Data[i];
2034 : } else {
2035 0 : if (min > Data[i])
2036 0 : min = Data[i];
2037 0 : if (max < Data[i])
2038 0 : max = Data[i];
2039 : }
2040 0 : locCnt++;
2041 : /* Fake a "group" by using the minimum group size. */
2042 0 : if (locCnt == minGroup) {
2043 0 : range0 += (max - min);
2044 0 : tot0++;
2045 0 : locCnt = 0;
2046 : }
2047 : }
2048 0 : if (locCnt != 0) {
2049 0 : range0 += (max - min);
2050 0 : tot0++;
2051 : }
2052 :
2053 0 : locCnt = 0;
2054 0 : range2 = 0;
2055 0 : tot2 = 0;
2056 : /* Compute scores for second order differences */
2057 0 : for (i = 0; i < numData; i++) {
2058 0 : if (locCnt == 0) {
2059 0 : min = SecDiff[i];
2060 0 : max = SecDiff[i];
2061 : } else {
2062 0 : if (min > SecDiff[i])
2063 0 : min = SecDiff[i];
2064 0 : if (max < SecDiff[i])
2065 0 : max = SecDiff[i];
2066 : }
2067 0 : locCnt++;
2068 : /* Fake a "group" by using the minimum group size. */
2069 0 : if (locCnt == minGroup) {
2070 0 : range2 += (max - min);
2071 0 : tot2++;
2072 0 : locCnt = 0;
2073 : }
2074 : }
2075 0 : if (locCnt != 0) {
2076 0 : range2 += (max - min);
2077 0 : tot2++;
2078 : }
2079 :
2080 : /* Compare average group size of no differencing to second order. */
2081 0 : if ((range0 / (tot0 + 0.0)) <= (range2 / (tot2 + 0.0))) {
2082 0 : return 0;
2083 : } else {
2084 0 : return 1;
2085 : }
2086 : }
2087 :
2088 : /*****************************************************************************
2089 : * power() --
2090 : *
2091 : * Arthur Taylor / MDL
2092 : *
2093 : * PURPOSE
2094 : * Calculate the number of bits required to store a given positive number.
2095 : *
2096 : * ARGUMENTS
2097 : * val = The number to store (Input)
2098 : * extra = number of slots to allocate for prim/sec missing values (Input)
2099 : *
2100 : * FILES/DATABASES: None
2101 : *
2102 : * RETURNS: int
2103 : * The number of bits needed to store this number.
2104 : *
2105 : * HISTORY
2106 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
2107 : *
2108 : * NOTES
2109 : *****************************************************************************
2110 : */
2111 0 : static int power (uInt4 val, int extra)
2112 : {
2113 : int i;
2114 :
2115 0 : val += extra;
2116 0 : if (val == 0) {
2117 0 : return 1;
2118 : }
2119 0 : for (i = 0; val != 0; i++) {
2120 0 : val = val >> 1;
2121 : }
2122 0 : return i;
2123 : }
2124 :
2125 : /*****************************************************************************
2126 : * findMaxMin2() --
2127 : *
2128 : * Arthur Taylor / MDL
2129 : *
2130 : * PURPOSE
2131 : * Find the min/max value between start/stop index values in the data
2132 : * Assuming primary and secondary missing values.
2133 : *
2134 : * ARGUMENTS
2135 : * Data = The data. (Input)
2136 : * start = The starting index in data (Input)
2137 : * stop = The stoping index in data (Input)
2138 : * li_primMiss = scaled primary missing value (Input)
2139 : * li_secMiss = scaled secondary missing value (Input)
2140 : * min = The min value found (Output)
2141 : * max = The max value found (Output)
2142 : *
2143 : * FILES/DATABASES: None
2144 : *
2145 : * RETURNS: char
2146 : * Flag if min/max are valid.
2147 : *
2148 : * HISTORY
2149 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
2150 : *
2151 : * NOTES
2152 : *****************************************************************************
2153 : */
2154 0 : static char findMaxMin2 (sInt4 *Data, int start, int stop,
2155 : sInt4 li_primMiss, sInt4 li_secMiss,
2156 : sInt4 *min, sInt4 *max)
2157 : {
2158 0 : char f_min = 0; /* Flag if we found the max/min values */
2159 : int i; /* Loop counter. */
2160 :
2161 0 : *max = *min = Data[start];
2162 0 : for (i = start; i < stop; i++) {
2163 0 : if ((Data[i] != li_secMiss) && (Data[i] != li_primMiss)) {
2164 0 : if (!f_min) {
2165 0 : *max = Data[i];
2166 0 : *min = Data[i];
2167 0 : f_min = 1;
2168 : } else {
2169 0 : if (*max < Data[i]) {
2170 0 : *max = Data[i];
2171 0 : } else if (*min > Data[i]) {
2172 0 : *min = Data[i];
2173 : }
2174 : }
2175 : }
2176 : }
2177 0 : return f_min;
2178 : }
2179 :
2180 : /*****************************************************************************
2181 : * findMaxMin1() --
2182 : *
2183 : * Arthur Taylor / MDL
2184 : *
2185 : * PURPOSE
2186 : * Find the min/max value between start/stop index values in the data
2187 : * Assuming primary missing values.
2188 : *
2189 : * ARGUMENTS
2190 : * Data = The data. (Input)
2191 : * start = The starting index in data (Input)
2192 : * stop = The stoping index in data (Input)
2193 : * li_primMiss = scaled primary missing value (Input)
2194 : * min = The min value found (Output)
2195 : * max = The max value found (Output)
2196 : *
2197 : * FILES/DATABASES: None
2198 : *
2199 : * RETURNS: char
2200 : * Flag if min/max are valid.
2201 : *
2202 : * HISTORY
2203 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
2204 : *
2205 : * NOTES
2206 : *****************************************************************************
2207 : */
2208 0 : static char findMaxMin1 (sInt4 *Data, int start, int stop,
2209 : sInt4 li_primMiss, sInt4 *min, sInt4 *max)
2210 : {
2211 0 : char f_min = 0; /* Flag if we found the max/min values */
2212 : int i; /* Loop counter. */
2213 :
2214 0 : *max = *min = Data[start];
2215 0 : for (i = start; i < stop; i++) {
2216 0 : if (Data[i] != li_primMiss) {
2217 0 : if (!f_min) {
2218 0 : *max = Data[i];
2219 0 : *min = Data[i];
2220 0 : f_min = 1;
2221 : } else {
2222 0 : if (*max < Data[i]) {
2223 0 : *max = Data[i];
2224 0 : } else if (*min > Data[i]) {
2225 0 : *min = Data[i];
2226 : }
2227 : }
2228 : }
2229 : }
2230 0 : return f_min;
2231 : }
2232 :
2233 : /*****************************************************************************
2234 : * findMaxMin0() --
2235 : *
2236 : * Arthur Taylor / MDL
2237 : *
2238 : * PURPOSE
2239 : * Find the min/max value between start/stop index values in the data
2240 : * Assuming no missing values.
2241 : *
2242 : * ARGUMENTS
2243 : * Data = The data. (Input)
2244 : * start = The starting index in data (Input)
2245 : * stop = The stoping index in data (Input)
2246 : * min = The min value found (Output)
2247 : * max = The max value found (Output)
2248 : *
2249 : * FILES/DATABASES: None
2250 : *
2251 : * RETURNS: void
2252 : *
2253 : * HISTORY
2254 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
2255 : *
2256 : * NOTES
2257 : *****************************************************************************
2258 : */
2259 0 : static void findMaxMin0 (sInt4 *Data, int start, int stop, sInt4 *min,
2260 : sInt4 *max)
2261 : {
2262 : int i; /* Loop counter. */
2263 :
2264 0 : *max = *min = Data[start];
2265 0 : for (i = start + 1; i < stop; i++) {
2266 0 : if (*max < Data[i]) {
2267 0 : *max = Data[i];
2268 0 : } else if (*min > Data[i]) {
2269 0 : *min = Data[i];
2270 : }
2271 : }
2272 0 : }
2273 :
2274 : /*****************************************************************************
2275 : * findGroup2() --
2276 : *
2277 : * Arthur Taylor / MDL
2278 : *
2279 : * PURPOSE
2280 : * Find "split" so that the numbers between start and split are within
2281 : * "range" of each other... stops if it reaches "stop".
2282 : * Assumes primary and secondary missing values.
2283 : *
2284 : * ARGUMENTS
2285 : * Data = The data. (Input)
2286 : * start = The starting index in data (Input)
2287 : * stop = The stoping index in data (Input)
2288 : * li_primMiss = scaled primary missing value (Input)
2289 : * li_secMiss = scaled secondary missing value (Input)
2290 : * range = The range to use (Input)
2291 : * split = The first index that is out of the range (Output)
2292 : * min = The min value for the group. (Output)
2293 : * max = The max value for the group. (Output)
2294 : *
2295 : * FILES/DATABASES: None
2296 : *
2297 : * RETURNS: void
2298 : *
2299 : * HISTORY
2300 : * 1/2005 Arthur Taylor (MDL): Created
2301 : *
2302 : * NOTES
2303 : *****************************************************************************
2304 : */
2305 0 : static void findGroup2 (sInt4 *Data, int start, int stop,
2306 : sInt4 li_primMiss, sInt4 li_secMiss,
2307 : sInt4 range, int *split, sInt4 *min, sInt4 *max)
2308 : {
2309 0 : char f_min = 0; /* Flag if we found the max/min values */
2310 : int i; /* Loop counter. */
2311 :
2312 0 : *min = *max = 0;
2313 0 : for (i = start; i < stop; i++) {
2314 0 : if ((Data[i] != li_secMiss) && (Data[i] != li_primMiss)) {
2315 0 : if (!f_min) {
2316 0 : *max = *min = Data[i];
2317 0 : f_min = 1;
2318 : } else {
2319 0 : if (*max < Data[i]) {
2320 0 : if ((Data[i] - *min) > range) {
2321 0 : *split = i;
2322 0 : return;
2323 : }
2324 0 : *max = Data[i];
2325 0 : } else if (*min > Data[i]) {
2326 0 : if ((*max - Data[i]) > range) {
2327 0 : *split = i;
2328 0 : return;
2329 : }
2330 0 : *min = Data[i];
2331 : }
2332 : }
2333 : }
2334 : }
2335 0 : *split = stop;
2336 : }
2337 :
2338 : /*****************************************************************************
2339 : * findGroup1() --
2340 : *
2341 : * Arthur Taylor / MDL
2342 : *
2343 : * PURPOSE
2344 : * Find "split" so that the numbers between start and split are within
2345 : * "range" of each other... stops if it reaches "stop".
2346 : * Assumes primary missing values.
2347 : *
2348 : * ARGUMENTS
2349 : * Data = The data. (Input)
2350 : * start = The starting index in data (Input)
2351 : * stop = The stoping index in data (Input)
2352 : * li_primMiss = scaled primary missing value (Input)
2353 : * range = The range to use (Input)
2354 : * split = The first index that is out of the range (Output)
2355 : * min = The min value for the group. (Output)
2356 : * max = The max value for the group. (Output)
2357 : *
2358 : * FILES/DATABASES: None
2359 : *
2360 : * RETURNS: void
2361 : *
2362 : * HISTORY
2363 : * 1/2005 Arthur Taylor (MDL): Created
2364 : *
2365 : * NOTES
2366 : *****************************************************************************
2367 : */
2368 0 : static void findGroup1 (sInt4 *Data, int start, int stop,
2369 : sInt4 li_primMiss, sInt4 range, int *split,
2370 : sInt4 *min, sInt4 *max)
2371 : {
2372 0 : char f_min = 0; /* Flag if we found the max/min values */
2373 : int i; /* Loop counter. */
2374 :
2375 0 : *min = *max = 0;
2376 0 : for (i = start; i < stop; i++) {
2377 0 : if (Data[i] != li_primMiss) {
2378 0 : if (!f_min) {
2379 0 : *max = *min = Data[i];
2380 0 : f_min = 1;
2381 : } else {
2382 0 : if (*max < Data[i]) {
2383 0 : if ((Data[i] - *min) > range) {
2384 0 : *split = i;
2385 0 : return;
2386 : }
2387 0 : *max = Data[i];
2388 0 : } else if (*min > Data[i]) {
2389 0 : if ((*max - Data[i]) > range) {
2390 0 : *split = i;
2391 0 : return;
2392 : }
2393 0 : *min = Data[i];
2394 : }
2395 : }
2396 : }
2397 : }
2398 0 : *split = stop;
2399 : }
2400 :
2401 : /*****************************************************************************
2402 : * findGroup0() --
2403 : *
2404 : * Arthur Taylor / MDL
2405 : *
2406 : * PURPOSE
2407 : * Find "split" so that the numbers between start and split are within
2408 : * "range" of each other... stops if it reaches "stop".
2409 : * Assumes no missing values.
2410 : *
2411 : * ARGUMENTS
2412 : * Data = The data. (Input)
2413 : * start = The starting index in data (Input)
2414 : * stop = The stoping index in data (Input)
2415 : * range = The range to use (Input)
2416 : * split = The first index that is out of the range (Output)
2417 : * min = The min value for the group. (Output)
2418 : * max = The max value for the group. (Output)
2419 : *
2420 : * FILES/DATABASES: None
2421 : *
2422 : * RETURNS: void
2423 : *
2424 : * HISTORY
2425 : * 1/2005 Arthur Taylor (MDL): Created
2426 : *
2427 : * NOTES
2428 : *****************************************************************************
2429 : */
2430 0 : static void findGroup0 (sInt4 *Data, int start, int stop,
2431 : sInt4 range, int *split, sInt4 *min, sInt4 *max)
2432 : {
2433 : int i; /* Loop counter. */
2434 :
2435 0 : *max = *min = Data[0];
2436 0 : for (i = start + 1; i < stop; i++) {
2437 0 : if (*max < Data[i]) {
2438 0 : if ((Data[i] - *min) > range) {
2439 0 : *split = i;
2440 0 : return;
2441 : }
2442 0 : *max = Data[i];
2443 0 : } else if (*min > Data[i]) {
2444 0 : if ((*max - Data[i]) > range) {
2445 0 : *split = i;
2446 0 : return;
2447 : }
2448 0 : *min = Data[i];
2449 : }
2450 : }
2451 0 : *split = stop;
2452 : }
2453 :
2454 : /*****************************************************************************
2455 : * findGroupRev2() --
2456 : *
2457 : * Arthur Taylor / MDL
2458 : *
2459 : * PURPOSE
2460 : * Find "split" so that the numbers between split and stop are within
2461 : * "range" of each other... stops if it reaches "start".
2462 : * Assumes primary and secondary missing values.
2463 : *
2464 : * ARGUMENTS
2465 : * Data = The data. (Input)
2466 : * start = The starting index in data (Input)
2467 : * stop = The stoping index in data (Input)
2468 : * li_primMiss = scaled primary missing value (Input)
2469 : * li_secMiss = scaled secondary missing value (Input)
2470 : * range = The range to use (Input)
2471 : * split = The first index that is still in the range (Output)
2472 : * min = The min value for the group. (Output)
2473 : * max = The max value for the group. (Output)
2474 : *
2475 : * FILES/DATABASES: None
2476 : *
2477 : * RETURNS: void
2478 : *
2479 : * HISTORY
2480 : * 1/2005 Arthur Taylor (MDL): Created
2481 : *
2482 : * NOTES
2483 : *****************************************************************************
2484 : */
2485 0 : static void findGroupRev2 (sInt4 *Data, int start, int stop,
2486 : sInt4 li_primMiss, sInt4 li_secMiss,
2487 : sInt4 range, int *split, sInt4 *min, sInt4 *max)
2488 : {
2489 0 : char f_min = 0; /* Flag if we found the max/min values */
2490 : int i; /* Loop counter. */
2491 :
2492 0 : *min = *max = 0;
2493 0 : for (i = stop - 1; i >= start; i--) {
2494 0 : if ((Data[i] != li_secMiss) && (Data[i] != li_primMiss)) {
2495 0 : if (!f_min) {
2496 0 : *max = *min = Data[i];
2497 0 : f_min = 1;
2498 : } else {
2499 0 : if (*max < Data[i]) {
2500 0 : if ((Data[i] - *min) > range) {
2501 0 : *split = i + 1;
2502 0 : return;
2503 : }
2504 0 : *max = Data[i];
2505 0 : } else if (*min > Data[i]) {
2506 0 : if ((*max - Data[i]) > range) {
2507 0 : *split = i + 1;
2508 0 : return;
2509 : }
2510 0 : *min = Data[i];
2511 : }
2512 : }
2513 : }
2514 : }
2515 0 : *split = start;
2516 : }
2517 :
2518 : /*****************************************************************************
2519 : * findGroupRev1() --
2520 : *
2521 : * Arthur Taylor / MDL
2522 : *
2523 : * PURPOSE
2524 : * Find "split" so that the numbers between split and stop are within
2525 : * "range" of each other... stops if it reaches "start".
2526 : * Assumes primary missing values.
2527 : *
2528 : * ARGUMENTS
2529 : * Data = The data. (Input)
2530 : * start = The starting index in data (Input)
2531 : * stop = The stoping index in data (Input)
2532 : * li_primMiss = scaled primary missing value (Input)
2533 : * range = The range to use (Input)
2534 : * split = The first index that is still in the range (Output)
2535 : * min = The min value for the group. (Output)
2536 : * max = The max value for the group. (Output)
2537 : *
2538 : * FILES/DATABASES: None
2539 : *
2540 : * RETURNS: void
2541 : *
2542 : * HISTORY
2543 : * 1/2005 Arthur Taylor (MDL): Created
2544 : *
2545 : * NOTES
2546 : *****************************************************************************
2547 : */
2548 0 : static void findGroupRev1 (sInt4 *Data, int start, int stop,
2549 : sInt4 li_primMiss, sInt4 range, int *split,
2550 : sInt4 *min, sInt4 *max)
2551 : {
2552 0 : char f_min = 0; /* Flag if we found the max/min values */
2553 : int i; /* Loop counter. */
2554 :
2555 0 : *min = *max = 0;
2556 0 : for (i = stop - 1; i >= start; i--) {
2557 0 : if (Data[i] != li_primMiss) {
2558 0 : if (!f_min) {
2559 0 : *max = *min = Data[i];
2560 0 : f_min = 1;
2561 : } else {
2562 0 : if (*max < Data[i]) {
2563 0 : if ((Data[i] - *min) > range) {
2564 0 : *split = i + 1;
2565 0 : return;
2566 : }
2567 0 : *max = Data[i];
2568 0 : } else if (*min > Data[i]) {
2569 0 : if ((*max - Data[i]) > range) {
2570 0 : *split = i + 1;
2571 0 : return;
2572 : }
2573 0 : *min = Data[i];
2574 : }
2575 : }
2576 : }
2577 : }
2578 0 : *split = start;
2579 : }
2580 :
2581 : /*****************************************************************************
2582 : * findGroupRev0() --
2583 : *
2584 : * Arthur Taylor / MDL
2585 : *
2586 : * PURPOSE
2587 : * Find "split" so that the numbers between split and stop are within
2588 : * "range" of each other... stops if it reaches "start".
2589 : * Assumes no missing values.
2590 : *
2591 : * ARGUMENTS
2592 : * Data = The data. (Input)
2593 : * start = The starting index in data (Input)
2594 : * stop = The stoping index in data (Input)
2595 : * li_primMiss = scaled primary missing value (Input)
2596 : * range = The range to use (Input)
2597 : * split = The first index that is still in the range (Output)
2598 : * min = The min value for the group. (Output)
2599 : * max = The max value for the group. (Output)
2600 : *
2601 : * FILES/DATABASES: None
2602 : *
2603 : * RETURNS: void
2604 : *
2605 : * HISTORY
2606 : * 1/2005 Arthur Taylor (MDL): Created
2607 : *
2608 : * NOTES
2609 : *****************************************************************************
2610 : */
2611 0 : static void findGroupRev0 (sInt4 *Data, int start, int stop,
2612 : sInt4 range, int *split, sInt4 *min, sInt4 *max)
2613 : {
2614 : int i; /* Loop counter. */
2615 :
2616 0 : *max = *min = Data[stop - 1];
2617 0 : for (i = stop - 2; i >= start; i--) {
2618 0 : if (*max < Data[i]) {
2619 0 : if ((Data[i] - *min) > range) {
2620 0 : *split = i + 1;
2621 0 : return;
2622 : }
2623 0 : *max = Data[i];
2624 0 : } else if (*min > Data[i]) {
2625 0 : if ((*max - Data[i]) > range) {
2626 0 : *split = i + 1;
2627 0 : return;
2628 : }
2629 0 : *min = Data[i];
2630 : }
2631 : }
2632 0 : *split = start;
2633 : }
2634 :
2635 : /*****************************************************************************
2636 : * shiftGroup2() --
2637 : *
2638 : * Arthur Taylor / MDL
2639 : *
2640 : * PURPOSE
2641 : * Find "split" so that the numbers between split and start1 are all inside
2642 : * the range defined by max, min and bit. It allows max and min to change,
2643 : * as long as it doesn't exceed the range defined by bit.
2644 : * This is very similar to findGroupRev?() but here we already know
2645 : * information about the min/max values, and are just expanding the group a
2646 : * little, while the other one knew nothing about the group, and just wanted
2647 : * a group of a given range.
2648 : * Assumes primary and secondary missing values.
2649 : *
2650 : * ARGUMENTS
2651 : * Data = The data. (Input)
2652 : * start1 = The starting index in data (Input)
2653 : * start2 = The starting index of the earlier group (ie don't go to any
2654 : * earlier indicies than this. (Input)
2655 : * li_primMiss = scaled primary missing value (Input)
2656 : * li_secMiss = scaled secondary missing value (Input)
2657 : * bit = The range we are allowed to store this in. (Input)
2658 : * min = The min value for the group. (Input/Output)
2659 : * max = The max value for the group. (Input/Output)
2660 : * split = The first index that is still in the range (Output)
2661 : *
2662 : * FILES/DATABASES: None
2663 : *
2664 : * RETURNS: void
2665 : *
2666 : * HISTORY
2667 : * 1/2005 Arthur Taylor (MDL): Created
2668 : *
2669 : * NOTES
2670 : *****************************************************************************
2671 : */
2672 0 : static void shiftGroup2 (sInt4 *Data, int start1, int start2,
2673 : sInt4 li_primMiss, sInt4 li_secMiss, int bit,
2674 : sInt4 *min, sInt4 *max, size_t *split)
2675 : {
2676 : int i; /* Loop counter. */
2677 : int range; /* The range defined by bit. */
2678 :
2679 0 : range = (int) (pow (2.0, bit) - 1) - 1;
2680 : myAssert (start2 <= start1);
2681 0 : for (i = start1; i >= start2; i--) {
2682 0 : if ((Data[i] != li_primMiss) && (Data[i] != li_secMiss)) {
2683 0 : if (Data[i] > *max) {
2684 0 : if ((Data[i] - *min) <= range) {
2685 0 : *max = Data[i];
2686 : } else {
2687 0 : *split = i + 1;
2688 0 : return;
2689 : }
2690 0 : } else if (Data[i] < *min) {
2691 0 : if ((*max - Data[i]) <= range) {
2692 0 : *min = Data[i];
2693 : } else {
2694 0 : *split = i + 1;
2695 0 : return;
2696 : }
2697 : }
2698 : }
2699 : }
2700 0 : *split = start2;
2701 : }
2702 :
2703 : /*****************************************************************************
2704 : * shiftGroup1() --
2705 : *
2706 : * Arthur Taylor / MDL
2707 : *
2708 : * PURPOSE
2709 : * Find "split" so that the numbers between split and start1 are all inside
2710 : * the range defined by max, min and bit. It allows max and min to change,
2711 : * as long as it doesn't exceed the range defined by bit.
2712 : * This is very similar to findGroupRev?() but here we already know
2713 : * information about the min/max values, and are just expanding the group a
2714 : * little, while the other one knew nothing about the group, and just wanted
2715 : * a group of a given range.
2716 : * Assumes primary missing values.
2717 : *
2718 : * ARGUMENTS
2719 : * Data = The data. (Input)
2720 : * start1 = The starting index in data (Input)
2721 : * start2 = The starting index of the earlier group (ie don't go to any
2722 : * earlier indicies than this. (Input)
2723 : * li_primMiss = scaled primary missing value (Input)
2724 : * bit = The range we are allowed to store this in. (Input)
2725 : * min = The min value for the group. (Input/Output)
2726 : * max = The max value for the group. (Input/Output)
2727 : * split = The first index that is still in the range (Output)
2728 : *
2729 : * FILES/DATABASES: None
2730 : *
2731 : * RETURNS: void
2732 : *
2733 : * HISTORY
2734 : * 1/2005 Arthur Taylor (MDL): Created
2735 : *
2736 : * NOTES
2737 : *****************************************************************************
2738 : */
2739 0 : static void shiftGroup1 (sInt4 *Data, int start1, int start2,
2740 : sInt4 li_primMiss, int bit,
2741 : sInt4 *min, sInt4 *max, size_t *split)
2742 : {
2743 : int i; /* Loop counter. */
2744 : int range; /* The range defined by bit. */
2745 :
2746 0 : range = (int) (pow (2.0, bit) - 1) - 1;
2747 : myAssert (start2 <= start1);
2748 0 : for (i = start1; i >= start2; i--) {
2749 0 : if (Data[i] != li_primMiss) {
2750 0 : if (Data[i] > *max) {
2751 0 : if ((Data[i] - *min) <= range) {
2752 0 : *max = Data[i];
2753 : } else {
2754 0 : *split = i + 1;
2755 0 : return;
2756 : }
2757 0 : } else if (Data[i] < *min) {
2758 0 : if ((*max - Data[i]) <= range) {
2759 0 : *min = Data[i];
2760 : } else {
2761 0 : *split = i + 1;
2762 0 : return;
2763 : }
2764 : }
2765 : }
2766 : }
2767 0 : *split = start2;
2768 : }
2769 :
2770 : /*****************************************************************************
2771 : * shiftGroup0() --
2772 : *
2773 : * Arthur Taylor / MDL
2774 : *
2775 : * PURPOSE
2776 : * Find "split" so that the numbers between split and start1 are all inside
2777 : * the range defined by max, min and bit. It allows max and min to change,
2778 : * as long as it doesn't exceed the range defined by bit.
2779 : * This is very similar to findGroupRev?() but here we already know
2780 : * information about the min/max values, and are just expanding the group a
2781 : * little, while the other one knew nothing about the group, and just wanted
2782 : * a group of a given range.
2783 : * Assumes no missing values.
2784 : *
2785 : * ARGUMENTS
2786 : * Data = The data. (Input)
2787 : * start1 = The starting index in data (Input)
2788 : * start2 = The starting index of the earlier group (ie don't go to any
2789 : * earlier indicies than this. (Input)
2790 : * bit = The range we are allowed to store this in. (Input)
2791 : * min = The min value for the group. (Input/Output)
2792 : * max = The max value for the group. (Input/Output)
2793 : * split = The first index that is still in the range (Output)
2794 : *
2795 : * FILES/DATABASES: None
2796 : *
2797 : * RETURNS: void
2798 : *
2799 : * HISTORY
2800 : * 1/2005 Arthur Taylor (MDL): Created
2801 : *
2802 : * NOTES
2803 : *****************************************************************************
2804 : */
2805 0 : static void shiftGroup0 (sInt4 *Data, int start1, int start2, int bit,
2806 : sInt4 *min, sInt4 *max, size_t *split)
2807 : {
2808 : int i; /* Loop counter. */
2809 : int range; /* The range defined by bit. */
2810 :
2811 0 : range = (int) (pow (2.0, bit) - 1) - 0;
2812 : myAssert (start2 <= start1);
2813 0 : for (i = start1; i >= start2; i--) {
2814 0 : if (Data[i] > *max) {
2815 0 : if ((Data[i] - *min) <= range) {
2816 0 : *max = Data[i];
2817 : } else {
2818 0 : *split = i + 1;
2819 0 : return;
2820 : }
2821 0 : } else if (Data[i] < *min) {
2822 0 : if ((*max - Data[i]) <= range) {
2823 0 : *min = Data[i];
2824 : } else {
2825 0 : *split = i + 1;
2826 0 : return;
2827 : }
2828 : }
2829 : }
2830 0 : *split = start2;
2831 : }
2832 :
2833 : /*****************************************************************************
2834 : * doSplit() --
2835 : *
2836 : * Arthur Taylor / MDL
2837 : *
2838 : * PURPOSE
2839 : * Reduce the "bit range", and create groups that grab as much as they can
2840 : * to the right. Then reduce those groups if they improve the score.
2841 : *
2842 : * ARGUMENTS
2843 : * Data = The data. (Input)
2844 : * numData = The number of elements in data. (Input)
2845 : * G = The group to split. (Input)
2846 : * lclGroup = The resulting groups (Output)
2847 : * numLclGroup = The number of resulting groups. (Output)
2848 : * f_primMiss = Flag if we have a primary missing value (Input)
2849 : * li_primMiss = scaled primary missing value (Input)
2850 : * f_secMiss = Flag if we have a secondary missing value (Input)
2851 : * li_secMiss = scaled secondary missing value (Input)
2852 : * xFactor = Estimate of cost (in bits) of a group. (Input)
2853 : *
2854 : * FILES/DATABASES: None
2855 : *
2856 : * RETURNS: void
2857 : *
2858 : * HISTORY
2859 : * 1/2005 Arthur Taylor (MDL): Created.
2860 : *
2861 : * NOTES
2862 : *****************************************************************************
2863 : */
2864 0 : static void doSplit (sInt4 *Data, int numData, TDLGroupType * G,
2865 : TDLGroupType ** lclGroup, int *numLclGroup,
2866 : char f_primMiss, sInt4 li_primMiss,
2867 : char f_secMiss, sInt4 li_secMiss, int xFactor)
2868 : {
2869 : int start; /* Where to start the current group. */
2870 : int range; /* The range to make the groups. */
2871 : int final; /* One more than the last index in the group G. */
2872 : int split; /* Where to split the group. */
2873 : TDLGroupType G1; /* The current group to add. */
2874 : TDLGroupType G2; /* The group if we evaporated the previous group. */
2875 : int evaporate; /* How many groups we have "evaporated". */
2876 : int i; /* Loop counter to help "evaporate" groups. */
2877 : sInt4 scoreA; /* The original score for 2 groups */
2878 : sInt4 scoreB; /* The new score (having evaporated a group) */
2879 : int GroupLen; /* Actual alloc'ed group len. */
2880 :
2881 : /* *INDENT-OFF* */
2882 : /* The (pow (2, ..) -1) is because 2^n - 1 is max range of a group.
2883 : * Example n = 1, 2^1 -1 = 1 range of (1,0) is 1.
2884 : * Example n = 2, 2^2 -1 = 3 range of (3,2,1,0) is 3.
2885 : * The G.bit - 1 is because we are trying to reduce the range. */
2886 : /* *INDENT-ON* */
2887 0 : range = (int) (pow (2.0, G->bit - 1) - 1) - (f_secMiss + f_primMiss);
2888 0 : split = G->start;
2889 0 : start = G->start;
2890 0 : final = G->start + G->num;
2891 : myAssert (final <= numData);
2892 0 : *numLclGroup = 0;
2893 0 : GroupLen = 1;
2894 0 : *lclGroup = (TDLGroupType *) malloc (GroupLen * sizeof (TDLGroupType));
2895 0 : while (split < final) {
2896 0 : if (f_secMiss) {
2897 : findGroup2 (Data, start, final, li_primMiss, li_secMiss, range,
2898 0 : &split, &(G1.min), &(G1.max));
2899 0 : } else if (f_primMiss) {
2900 : findGroup1 (Data, start, final, li_primMiss, range, &split,
2901 0 : &(G1.min), &(G1.max));
2902 : } else {
2903 0 : findGroup0 (Data, start, final, range, &split, &(G1.min), &(G1.max));
2904 : }
2905 : G1.bit = (char) power ((uInt4) (G1.max - G1.min),
2906 0 : f_secMiss + f_primMiss);
2907 0 : G1.num = split - start;
2908 0 : G1.start = start;
2909 0 : G1.f_trySplit = 1;
2910 0 : G1.f_tryShift = 1;
2911 : /* Test if we should add to previous group, or create a new group. */
2912 0 : if (*numLclGroup == 0) {
2913 0 : *numLclGroup = 1;
2914 0 : (*lclGroup)[0] = G1;
2915 : } else {
2916 0 : G2.start = (*lclGroup)[*numLclGroup - 1].start;
2917 0 : G2.num = (*lclGroup)[*numLclGroup - 1].num + G1.num;
2918 0 : G2.min = ((*lclGroup)[*numLclGroup - 1].min < G1.min) ?
2919 0 : (*lclGroup)[*numLclGroup - 1].min : G1.min;
2920 0 : G2.max = ((*lclGroup)[*numLclGroup - 1].max > G1.max) ?
2921 0 : (*lclGroup)[*numLclGroup - 1].max : G1.max;
2922 : G2.bit = (char) power ((uInt4) (G2.max - G2.min),
2923 0 : f_secMiss + f_primMiss);
2924 0 : G2.f_trySplit = 1;
2925 0 : G2.f_tryShift = 1;
2926 0 : scoreA = ((*lclGroup)[*numLclGroup - 1].bit *
2927 0 : (*lclGroup)[*numLclGroup - 1].num) + xFactor;
2928 0 : scoreA += G1.bit * G1.num + xFactor;
2929 0 : scoreB = G2.bit * G2.num + xFactor;
2930 0 : if (scoreB < scoreA) {
2931 0 : (*lclGroup)[*numLclGroup - 1] = G2;
2932 : /* See if we can evaporate any of the old groups */
2933 0 : evaporate = 0;
2934 0 : for (i = *numLclGroup - 1; i > 0; i--) {
2935 0 : G1.start = (*lclGroup)[i - 1].start;
2936 0 : G1.num = (*lclGroup)[i].num + (*lclGroup)[i - 1].num;
2937 0 : G1.min = ((*lclGroup)[i].min < (*lclGroup)[i - 1].min) ?
2938 0 : (*lclGroup)[i].min : (*lclGroup)[i - 1].min;
2939 0 : G1.max = ((*lclGroup)[i].max > (*lclGroup)[i - 1].max) ?
2940 0 : (*lclGroup)[i].max : (*lclGroup)[i - 1].max;
2941 : G1.bit = (char) power ((uInt4) (G1.max - G1.min),
2942 0 : f_secMiss + f_primMiss);
2943 0 : G1.f_trySplit = 1;
2944 0 : G1.f_tryShift = 1;
2945 0 : scoreA = (*lclGroup)[i].bit * (*lclGroup)[i].num + xFactor;
2946 0 : scoreA += ((*lclGroup)[i - 1].bit * (*lclGroup)[i - 1].num +
2947 0 : xFactor);
2948 0 : scoreB = G1.bit * G1.num + xFactor;
2949 0 : if (scoreB < scoreA) {
2950 0 : evaporate++;
2951 0 : (*lclGroup)[i - 1] = G1;
2952 : } else {
2953 0 : break;
2954 : }
2955 : }
2956 0 : if (evaporate != 0) {
2957 0 : *numLclGroup = *numLclGroup - evaporate;
2958 : }
2959 : } else {
2960 0 : *numLclGroup = *numLclGroup + 1;
2961 0 : if (*numLclGroup > GroupLen) {
2962 0 : GroupLen = *numLclGroup;
2963 : *lclGroup = (TDLGroupType *) realloc ((void *) *lclGroup,
2964 : GroupLen *
2965 0 : sizeof (TDLGroupType));
2966 : }
2967 0 : (*lclGroup)[*numLclGroup - 1] = G1;
2968 : }
2969 : }
2970 0 : start = split;
2971 : }
2972 0 : }
2973 :
2974 : /*****************************************************************************
2975 : * doSplitRight() --
2976 : *
2977 : * Arthur Taylor / MDL
2978 : *
2979 : * PURPOSE
2980 : * Break into two groups right has range n - 1, left has range n.
2981 : *
2982 : * ARGUMENTS
2983 : * Data = The data. (Input)
2984 : * numData = The number of elements in data. (Input)
2985 : * G = The group to split. (Input)
2986 : * G1 = The right most group (Output)
2987 : * G2 = The remainder. (Output)
2988 : * f_primMiss = Flag if we have a primary missing value (Input)
2989 : * li_primMiss = scaled primary missing value (Input)
2990 : * f_secMiss = Flag if we have a secondary missing value (Input)
2991 : * li_secMiss = scaled secondary missing value (Input)
2992 : *
2993 : * FILES/DATABASES: None
2994 : *
2995 : * RETURNS: void
2996 : *
2997 : * HISTORY
2998 : * 1/2005 Arthur Taylor (MDL): Created.
2999 : *
3000 : * NOTES
3001 : *****************************************************************************
3002 : */
3003 0 : static void doSplitRight (sInt4 *Data, int numData, TDLGroupType * G,
3004 : TDLGroupType * G1, TDLGroupType * G2,
3005 : char f_primMiss, sInt4 li_primMiss,
3006 : char f_secMiss, sInt4 li_secMiss)
3007 : {
3008 : int range; /* The range to make the right most group. */
3009 : int final; /* One more than the last index in the group. */
3010 : int split; /* Where to split the group. */
3011 :
3012 : /* *INDENT-OFF* */
3013 : /* The (pow (2, ..) -1) is because 2^n - 1 is max range of a group.
3014 : * Example n = 1, 2^1 -1 = 1 range of (1,0) is 1.
3015 : * Example n = 2, 2^2 -1 = 3 range of (3,2,1,0) is 3.
3016 : * The G.bit - 1 is because we are trying to reduce the range. */
3017 : /* *INDENT-ON* */
3018 0 : range = (int) (pow (2.0, G->bit - 1) - 1) - (f_secMiss + f_primMiss);
3019 0 : final = G->start + G->num;
3020 0 : split = final;
3021 : myAssert (final <= numData);
3022 :
3023 0 : if (f_secMiss) {
3024 : findGroupRev2 (Data, G->start, final, li_primMiss, li_secMiss, range,
3025 0 : &split, &(G1->min), &(G1->max));
3026 : findMaxMin2 (Data, G->start, split, li_primMiss, li_secMiss,
3027 0 : &(G2->min), &(G2->max));
3028 0 : } else if (f_primMiss) {
3029 : findGroupRev1 (Data, G->start, final, li_primMiss, range, &split,
3030 0 : &(G1->min), &(G1->max));
3031 : findMaxMin1 (Data, G->start, split, li_primMiss, &(G2->min),
3032 0 : &(G2->max));
3033 : } else {
3034 : findGroupRev0 (Data, G->start, final, range, &split, &(G1->min),
3035 0 : &(G1->max));
3036 0 : findMaxMin0 (Data, G->start, split, &(G2->min), &(G2->max));
3037 : }
3038 :
3039 : G1->bit = (char) power ((uInt4) (G1->max - G1->min),
3040 0 : f_secMiss + f_primMiss);
3041 : G2->bit = (char) power ((uInt4) (G2->max - G2->min),
3042 0 : f_secMiss + f_primMiss);
3043 0 : G1->start = split;
3044 0 : G2->start = G->start;
3045 0 : G1->num = final - split;
3046 0 : G2->num = split - G->start;
3047 0 : G1->f_trySplit = 1;
3048 0 : G1->f_tryShift = 1;
3049 0 : G2->f_trySplit = 1;
3050 0 : G2->f_tryShift = 1;
3051 0 : }
3052 :
3053 : /*****************************************************************************
3054 : * ComputeGroupSize() --
3055 : *
3056 : * Arthur Taylor / MDL
3057 : *
3058 : * PURPOSE
3059 : * Compute the number of bits needed for the various elements of the groups
3060 : * as well as the total number of bits needed.
3061 : *
3062 : * ARGUMENTS
3063 : * group = Groups (Input)
3064 : * numGroup = Number of groups (Input)
3065 : * ibit = Number of bits needed for the minimum values of each group.
3066 : * Find max absolute value of group mins. Note: all group mins
3067 : * are positive (Output)
3068 : * jbit = Number of bits needed for the number of bits for each group.
3069 : * Find max absolute value of number of bits. (Output)
3070 : * kbit = Number of bits needed for the number of values for each group.
3071 : * Find max absolute value of number of values. (Output)
3072 : *
3073 : * FILES/DATABASES: None
3074 : *
3075 : * RETURNS: sInt4
3076 : * number of bits needed by the groups
3077 : *
3078 : * HISTORY
3079 : * 12/2004 Arthur Taylor (MDL): Updated from "tdlpack.c" in "C" tdlpack code
3080 : *
3081 : * NOTES
3082 : *****************************************************************************
3083 : */
3084 0 : static sInt4 ComputeGroupSize (TDLGroupType * group, int numGroup,
3085 : size_t *ibit, size_t *jbit, size_t *kbit)
3086 : {
3087 : int i; /* loop counter. */
3088 0 : sInt4 ans = 0; /* The number of bits needed. */
3089 0 : sInt4 maxMin = 0; /* The largest min value in the groups */
3090 0 : uChar maxBit = 0; /* The largest needed bits in the groups */
3091 0 : uInt4 maxNum = 0; /* The largest number of values in the groups. */
3092 :
3093 0 : for (i = 0; i < numGroup; i++) {
3094 0 : ans += group[i].bit * group[i].num;
3095 0 : if (group[i].min > maxMin) {
3096 0 : maxMin = group[i].min;
3097 : }
3098 0 : if (group[i].bit > maxBit) {
3099 0 : maxBit = group[i].bit;
3100 : }
3101 0 : if (group[i].num > maxNum) {
3102 0 : maxNum = group[i].num;
3103 : }
3104 : }
3105 : /* This only works for pos numbers... */
3106 0 : for (i = 0; (maxMin != 0); i++) {
3107 0 : maxMin = maxMin >> 1;
3108 : }
3109 : /* Allow 0 bits for min. Assumes that decoder allows 0 bits */
3110 0 : *ibit = i;
3111 : /* This only works for pos numbers... */
3112 0 : for (i = 0; (maxBit != 0); i++) {
3113 0 : maxBit = maxBit >> 1;
3114 : }
3115 : /* Allow 0 bits for min. Assumes that decoder allows 0 bits */
3116 0 : *jbit = i;
3117 : /* This only works for pos numbers... */
3118 0 : for (i = 0; (maxNum != 0); i++) {
3119 0 : maxNum = maxNum >> 1;
3120 : }
3121 : /* Allow 0 bits for min. Assumes that decoder allows 0 bits */
3122 0 : *kbit = i;
3123 0 : ans += ((*ibit) + (*jbit) + (*kbit)) * numGroup;
3124 0 : return ans;
3125 : }
3126 :
3127 : /*****************************************************************************
3128 : * splitGroup() --
3129 : *
3130 : * Arthur Taylor / MDL
3131 : *
3132 : * PURPOSE
3133 : * Tries to reduce (split) each group by 1 bit. It does this by:
3134 : * A) reduce the "bit range", and create groups that grab as much as they
3135 : * can to the right. Then reduce those groups if they improve the score.
3136 : * B) reduce the bit range and grab the left most group only, leaving the
3137 : * rest unchanged.
3138 : * C) reduce the bit range and grab the right most group only, leaving the
3139 : * rest unchanged.
3140 : *
3141 : * ARGUMENTS
3142 : * Data = The data. (Input)
3143 : * numData = The number of elements in data. (Input)
3144 : * group = The groups using the "best minGroup. (Input)
3145 : * numGroup = Number of groups (Input)
3146 : * lclGroup = The local copy of the groups (Output)
3147 : * numLclGroup = Number of local groups (Output)
3148 : * f_primMiss = Flag if we have a primary missing value (Input)
3149 : * li_primMiss = scaled primary missing value (Input)
3150 : * f_secMiss = Flag if we have a secondary missing value (Input)
3151 : * li_secMiss = scaled secondary missing value (Input)
3152 : * xFactor = Estimate of cost (in bits) of a group. (Input)
3153 : *
3154 : * FILES/DATABASES: None
3155 : *
3156 : * RETURNS: void
3157 : *
3158 : * HISTORY
3159 : * 1/2005 Arthur Taylor (MDL): Created.
3160 : *
3161 : * NOTES
3162 : *****************************************************************************
3163 : */
3164 0 : static int splitGroup (sInt4 *Data, int numData, TDLGroupType * group,
3165 : int numGroup, TDLGroupType ** lclGroup,
3166 : int *numLclGroup, char f_primMiss,
3167 : sInt4 li_primMiss, char f_secMiss,
3168 : sInt4 li_secMiss, size_t xFactor)
3169 : {
3170 : uInt4 minBit; /* The fewest # of bits, with no subdivision. */
3171 : TDLGroupType *subGroup; /* The subgroups that we tried splitting the
3172 : * primary group into. */
3173 : int numSubGroup; /* The number of groups in subGroup. */
3174 : sInt4 A_max; /* Max value of a given group. */
3175 : sInt4 A_min; /* Min value of a given group. */
3176 : sInt4 scoreA; /* The original score for a given group */
3177 : sInt4 scoreB; /* The new score */
3178 0 : int f_adjust = 0; /* Flag if group has changed. */
3179 : int f_keep; /* Flag to keep the subgroup instead of original */
3180 : int i; /* Loop counters */
3181 : int sub; /* loop counter over the sub group. */
3182 : int lclIndex; /* Used to help copy data from subGroup to answer. */
3183 : int GroupLen; /* Actual alloc'ed group len. */
3184 : int extra; /* Used to reduce number of allocs. */
3185 :
3186 : /* Figure out how few bits a group can have without being able to further
3187 : * divide it. */
3188 0 : if (f_secMiss) {
3189 : /* 11 = primMiss 10 = secMiss 01, 00 = data. */
3190 0 : minBit = 2;
3191 0 : } else if (f_primMiss) {
3192 : /* 1 = primMiss 0 = data. */
3193 : /* might try minBit = 1 here. */
3194 0 : minBit = 1;
3195 : } else {
3196 : /* 1, 0 = data. */
3197 0 : minBit = 1;
3198 : }
3199 :
3200 0 : *numLclGroup = 0;
3201 0 : *lclGroup = (TDLGroupType *) malloc (numGroup * sizeof (TDLGroupType));
3202 0 : GroupLen = numGroup;
3203 0 : extra = 0;
3204 0 : for (i = 0; i < numGroup; i++) {
3205 : /* Check if we have already tried to split this group, or it has too
3206 : * few members, or it doesn't have enough bits to split. If so, skip
3207 : * this group. */
3208 0 : if ((group[i].f_trySplit) && (group[i].num > xFactor) &&
3209 0 : (group[i].bit > minBit)) {
3210 0 : f_keep = 0;
3211 : doSplit (Data, numData, &(group[i]), &subGroup, &numSubGroup,
3212 0 : f_primMiss, li_primMiss, f_secMiss, li_secMiss, xFactor);
3213 0 : if (numSubGroup != 1) {
3214 0 : scoreA = group[i].bit * group[i].num + xFactor;
3215 0 : scoreB = 0;
3216 0 : for (sub = 0; sub < numSubGroup; sub++) {
3217 0 : scoreB += subGroup[sub].bit * subGroup[sub].num + xFactor;
3218 : }
3219 0 : if (scoreB < scoreA) {
3220 0 : f_keep = 1;
3221 0 : } else if (numSubGroup > 2) {
3222 : /* We can do "doSplitLeft" (which is breaking it into 2 groups
3223 : * the first having range n - 1, the second having range n,
3224 : * using what we know from doSplit. */
3225 0 : subGroup[1].num = group[i].num - subGroup[0].num;
3226 0 : if (f_secMiss) {
3227 0 : findMaxMin2 (Data, subGroup[1].start,
3228 0 : subGroup[1].start + subGroup[1].num,
3229 0 : li_primMiss, li_secMiss, &A_min, &A_max);
3230 0 : } else if (f_primMiss) {
3231 0 : findMaxMin1 (Data, subGroup[1].start,
3232 0 : subGroup[1].start + subGroup[1].num,
3233 0 : li_primMiss, &A_min, &A_max);
3234 : } else {
3235 0 : findMaxMin0 (Data, subGroup[1].start,
3236 0 : subGroup[1].start + subGroup[1].num,
3237 0 : &A_min, &A_max);
3238 : }
3239 0 : subGroup[1].min = A_min;
3240 0 : subGroup[1].max = A_max;
3241 0 : subGroup[1].bit = power ((uInt4) (A_max - A_min),
3242 0 : f_secMiss + f_primMiss);
3243 0 : subGroup[1].f_trySplit = 1;
3244 0 : subGroup[1].f_tryShift = 1;
3245 0 : numSubGroup = 2;
3246 0 : scoreB = subGroup[0].bit * subGroup[0].num + xFactor;
3247 0 : scoreB += subGroup[1].bit * subGroup[1].num + xFactor;
3248 0 : if (scoreB < scoreA) {
3249 0 : f_keep = 1;
3250 : }
3251 : }
3252 : }
3253 0 : if (!f_keep) {
3254 0 : if (numSubGroup == 1) {
3255 : subGroup = (TDLGroupType *) realloc (subGroup,
3256 : 2 *
3257 0 : sizeof (TDLGroupType));
3258 : }
3259 0 : numSubGroup = 2;
3260 : doSplitRight (Data, numData, &(group[i]), &(subGroup[1]),
3261 : &(subGroup[0]), f_primMiss, li_primMiss, f_secMiss,
3262 0 : li_secMiss);
3263 0 : scoreA = group[i].bit * group[i].num + xFactor;
3264 0 : scoreB = subGroup[0].bit * subGroup[0].num + xFactor;
3265 0 : scoreB += subGroup[1].bit * subGroup[1].num + xFactor;
3266 0 : if (scoreB < scoreA) {
3267 0 : f_keep = 1;
3268 : }
3269 : }
3270 0 : if (f_keep) {
3271 0 : lclIndex = *numLclGroup;
3272 0 : *numLclGroup = *numLclGroup + numSubGroup;
3273 0 : if (*numLclGroup > GroupLen) {
3274 0 : GroupLen += extra;
3275 0 : extra = 0;
3276 0 : if (*numLclGroup > GroupLen) {
3277 0 : GroupLen = *numLclGroup;
3278 : }
3279 : *lclGroup = (TDLGroupType *) realloc ((void *) *lclGroup,
3280 : GroupLen *
3281 0 : sizeof (TDLGroupType));
3282 : } else {
3283 0 : extra += numSubGroup - 1;
3284 : }
3285 : memcpy ((*lclGroup) + lclIndex, subGroup,
3286 0 : numSubGroup * sizeof (TDLGroupType));
3287 0 : f_adjust = 1;
3288 : } else {
3289 0 : *numLclGroup = *numLclGroup + 1;
3290 0 : if (*numLclGroup > GroupLen) {
3291 0 : GroupLen += extra;
3292 0 : extra = 0;
3293 0 : if (*numLclGroup > GroupLen) {
3294 0 : GroupLen = *numLclGroup;
3295 : }
3296 : *lclGroup = (TDLGroupType *) realloc ((void *) *lclGroup,
3297 : GroupLen *
3298 0 : sizeof (TDLGroupType));
3299 : }
3300 0 : (*lclGroup)[*numLclGroup - 1] = group[i];
3301 0 : (*lclGroup)[*numLclGroup - 1].f_trySplit = 0;
3302 : }
3303 0 : free (subGroup);
3304 0 : subGroup = NULL;
3305 : } else {
3306 0 : *numLclGroup = *numLclGroup + 1;
3307 0 : if (*numLclGroup > GroupLen) {
3308 0 : GroupLen += extra;
3309 0 : extra = 0;
3310 0 : if (*numLclGroup > GroupLen) {
3311 0 : GroupLen = *numLclGroup;
3312 : }
3313 : *lclGroup = (TDLGroupType *) realloc ((void *) *lclGroup,
3314 : GroupLen *
3315 0 : sizeof (TDLGroupType));
3316 : }
3317 0 : (*lclGroup)[*numLclGroup - 1] = group[i];
3318 0 : (*lclGroup)[*numLclGroup - 1].f_trySplit = 0;
3319 : }
3320 : }
3321 : myAssert (GroupLen == *numLclGroup);
3322 0 : return f_adjust;
3323 : }
3324 :
3325 : /*****************************************************************************
3326 : * shiftGroup() --
3327 : *
3328 : * Arthur Taylor / MDL
3329 : *
3330 : * PURPOSE
3331 : * Tries to shift / join the groups together. It does this by first
3332 : * calculating if a group should still exist. If it should, then it has
3333 : * each group grab as much as it can to the left without increasing its "bit
3334 : * range".
3335 : *
3336 : * ARGUMENTS
3337 : * Data = The data. (Input)
3338 : * numData = The number of elements in data. (Input)
3339 : * group = The resulting groups. (Output)
3340 : * numGroup = Number of groups (Output)
3341 : * f_primMiss = Flag if we have a primary missing value (Input)
3342 : * li_primMiss = scaled primary missing value (Input)
3343 : * f_secMiss = Flag if we have a secondary missing value (Input)
3344 : * li_secMiss = scaled secondary missing value (Input)
3345 : * xFactor = Estimate of cost (in bits) of a group. (Input)
3346 : *
3347 : * FILES/DATABASES: None
3348 : *
3349 : * RETURNS: void
3350 : *
3351 : * HISTORY
3352 : * 1/2005 Arthur Taylor (MDL): Created.
3353 : *
3354 : * NOTES
3355 : *****************************************************************************
3356 : */
3357 0 : static void shiftGroup (sInt4 *Data, int numData, TDLGroupType ** Group,
3358 : size_t *NumGroup, char f_primMiss, sInt4 li_primMiss,
3359 : char f_secMiss, sInt4 li_secMiss, int xFactor)
3360 : {
3361 0 : TDLGroupType *group = (*Group); /* Local pointer to Group. */
3362 0 : int numGroup = (*NumGroup); /* # elements in group. */
3363 : int i, j; /* loop counters. */
3364 : sInt4 A_max; /* Max value of a given group. */
3365 : sInt4 A_min; /* Min value of a given group. */
3366 : size_t begin; /* New start to the group. */
3367 : sInt4 scoreA; /* The original score for group i and i - 1 */
3368 : sInt4 scoreB; /* The new score for group i and i - 1 */
3369 : TDLGroupType G1; /* The "new" group[i - 1]. */
3370 : TDLGroupType G2; /* The "new" group[i]. */
3371 0 : int evaporate = 0; /* number of groups that "evaporated" */
3372 : int index; /* index while getting rid of "evaporated groups". */
3373 :
3374 0 : for (i = numGroup - 1; i > 0; i--) {
3375 : myAssert (group[i].num > 0);
3376 : /* See if we can evaporate the group n - 1 */
3377 0 : G1.start = group[i - 1].start;
3378 0 : G1.num = group[i].num + group[i - 1].num;
3379 0 : G1.min = (group[i].min < group[i - 1].min) ?
3380 0 : group[i].min : group[i - 1].min;
3381 0 : G1.max = (group[i].max > group[i - 1].max) ?
3382 0 : group[i].max : group[i - 1].max;
3383 : G1.bit = (char) power ((uInt4) (G1.max - G1.min),
3384 0 : f_secMiss + f_primMiss);
3385 0 : G1.f_trySplit = 1;
3386 0 : G1.f_tryShift = 1;
3387 0 : scoreA = group[i].bit * group[i].num + xFactor;
3388 0 : scoreA += group[i - 1].bit * group[i - 1].num + xFactor;
3389 0 : scoreB = G1.bit * G1.num + xFactor;
3390 0 : if (scoreB < scoreA) {
3391 : /* One of the groups evaporated. */
3392 0 : evaporate++;
3393 0 : group[i - 1] = G1;
3394 0 : group[i].num = 0;
3395 : /* See if that affects any of the previous groups. */
3396 0 : for (j = i + 1; j < numGroup; j++) {
3397 0 : if (group[j].num != 0) {
3398 0 : G1.start = G1.start;
3399 0 : G1.num = group[i - 1].num + group[j].num;
3400 0 : G1.min = (group[i - 1].min < group[j].min) ?
3401 0 : group[i - 1].min : group[j].min;
3402 0 : G1.max = (group[i - 1].max > group[j].max) ?
3403 0 : group[i - 1].max : group[j].max;
3404 : G1.bit = (char) power ((uInt4) (G1.max - G1.min),
3405 0 : f_secMiss + f_primMiss);
3406 0 : G1.f_trySplit = 1;
3407 0 : G1.f_tryShift = 1;
3408 0 : scoreA = group[i - 1].bit * group[i - 1].num + xFactor;
3409 0 : scoreA += group[j].bit * group[j].num + xFactor;
3410 0 : scoreB = G1.bit * G1.num + xFactor;
3411 0 : if (scoreB < scoreA) {
3412 0 : evaporate++;
3413 0 : group[i - 1] = G1;
3414 0 : group[j].num = 0;
3415 : } else {
3416 0 : break;
3417 : }
3418 : }
3419 : }
3420 0 : } else if (group[i].f_tryShift) {
3421 : /* Group did not evaporate, so do the "grabby" algorithm. */
3422 0 : if ((group[i].bit != 0) && (group[i - 1].bit >= group[i].bit)) {
3423 0 : if (f_secMiss) {
3424 0 : A_max = group[i].max;
3425 0 : A_min = group[i].min;
3426 0 : shiftGroup2 (Data, group[i].start - 1, group[i - 1].start,
3427 0 : li_primMiss, li_secMiss, group[i].bit, &A_min,
3428 0 : &A_max, &begin);
3429 0 : } else if (f_primMiss) {
3430 0 : A_max = group[i].max;
3431 0 : A_min = group[i].min;
3432 0 : shiftGroup1 (Data, group[i].start - 1, group[i - 1].start,
3433 0 : li_primMiss, group[i].bit, &A_min, &A_max,
3434 0 : &begin);
3435 : } else {
3436 0 : A_max = group[i].max;
3437 0 : A_min = group[i].min;
3438 0 : shiftGroup0 (Data, group[i].start - 1, group[i - 1].start,
3439 0 : group[i].bit, &A_min, &A_max, &begin);
3440 : }
3441 0 : if (begin != group[i].start) {
3442 : /* Re-Calculate min/max of group[i - 1], since it could have
3443 : * moved to group[i] */
3444 0 : G1 = group[i - 1];
3445 0 : G2 = group[i];
3446 0 : G2.min = A_min;
3447 0 : G2.max = A_max;
3448 0 : G1.num -= (group[i].start - begin);
3449 0 : if (f_secMiss) {
3450 : findMaxMin2 (Data, G1.start, G1.start + G1.num,
3451 0 : li_primMiss, li_secMiss, &A_min, &A_max);
3452 0 : } else if (f_primMiss) {
3453 : findMaxMin1 (Data, G1.start, G1.start + G1.num,
3454 0 : li_primMiss, &A_min, &A_max);
3455 : } else {
3456 : findMaxMin0 (Data, G1.start, G1.start + G1.num,
3457 0 : &A_min, &A_max);
3458 : }
3459 0 : if ((A_min != G1.min) || (A_max != G1.max)) {
3460 0 : G1.min = A_min;
3461 0 : G1.max = A_max;
3462 : G1.bit = (char) power ((uInt4) (A_max - A_min),
3463 0 : f_secMiss + f_primMiss);
3464 0 : G1.f_trySplit = 1;
3465 0 : G1.f_tryShift = 1;
3466 : }
3467 0 : G2.num += (group[i].start - begin);
3468 0 : G2.start = begin;
3469 0 : G2.f_trySplit = 1;
3470 0 : G2.f_tryShift = 1;
3471 0 : scoreA = group[i].bit * group[i].num + xFactor;
3472 0 : scoreA += group[i - 1].bit * group[i - 1].num + xFactor;
3473 0 : scoreB = G2.bit * G2.num + xFactor;
3474 0 : if (G1.num != 0) {
3475 0 : scoreB += G1.bit * G1.num + xFactor;
3476 : }
3477 : #ifdef DEBUG
3478 : if (scoreB > scoreA) {
3479 : printf ("Made score worse!\n");
3480 : }
3481 : #endif
3482 0 : if (begin == group[i - 1].start) {
3483 : /* Grabby algorithm evaporated a group. */
3484 0 : evaporate++;
3485 : myAssert (G1.num == 0);
3486 : /* Switch the evaporating group to other side so we have
3487 : * potential to continue evaporating the next group down. */
3488 0 : group[i - 1] = G2;
3489 0 : group[i] = G1;
3490 : } else {
3491 0 : group[i - 1] = G1;
3492 0 : group[i] = G2;
3493 : }
3494 : } else {
3495 0 : group[i].f_tryShift = 0;
3496 : }
3497 : }
3498 : }
3499 : }
3500 : /* Loop through the grid removing the evaporated groups. */
3501 0 : if (evaporate != 0) {
3502 0 : index = 0;
3503 0 : for (i = 0; i < numGroup; i++) {
3504 0 : if (group[i].num != 0) {
3505 0 : group[index] = group[i];
3506 0 : index++;
3507 : }
3508 : }
3509 0 : *NumGroup = numGroup - evaporate;
3510 : *Group = (TDLGroupType *) realloc ((void *) (*Group),
3511 0 : *NumGroup * sizeof (TDLGroupType));
3512 : }
3513 0 : }
3514 :
3515 : /*****************************************************************************
3516 : * GroupIt() --
3517 : *
3518 : * Arthur Taylor / MDL
3519 : *
3520 : * PURPOSE
3521 : * Attempts to find groups for packing the data. It starts by preparing
3522 : * the data, by removing the overall min value.
3523 : *
3524 : * Next it Creates any 0 bit groups (primary missing: missings are all 0
3525 : * bit groups, const values are 0 bit groups. No missing: const values are
3526 : * 0 bit groups.)
3527 : *
3528 : * Next it tries to reduce (split) each group by 1 bit. It does this by:
3529 : * A) reduce the "bit range", and create groups that grab as much as they
3530 : * can to the right. Then reduce those groups if they improve the score.
3531 : * B) reduce the bit range and grab the left most group only, leaving the
3532 : * rest unchanged.
3533 : * C) reduce the bit range and grab the right most group only, leaving the
3534 : * rest unchanged.
3535 : *
3536 : * Next it tries to shift / join those groups together. It does this by
3537 : * first calculating if a group should still exist. If it should, then it
3538 : * has each group grab as much as it can to the left without increasing its
3539 : * "bit range".
3540 : *
3541 : * ARGUMENTS
3542 : * OverallMin = The overall min value in the data. (Input)
3543 : * Data = The data. (Input)
3544 : * numData = The number of elements in data. (Input)
3545 : * group = The resulting groups. (Output)
3546 : * numGroup = Number of groups (Output)
3547 : * f_primMiss = Flag if we have a primary missing value (Input)
3548 : * li_primMiss = scaled primary missing value (Input)
3549 : * f_secMiss = Flag if we have a secondary missing value (Input)
3550 : * li_secMiss = scaled secondary missing value (Input)
3551 : * groupSize = How many bytes the groups and data will take. (Output)
3552 : * ibit = # of bits for largest minimum value in groups (Output)
3553 : * jbit = # of bits for largest # bits in groups (Output)
3554 : * kbit = # of bits for largest # values in groups (Output)
3555 : *
3556 : * FILES/DATABASES: None
3557 : *
3558 : * RETURNS: void
3559 : *
3560 : * HISTORY
3561 : * 1/2005 Arthur Taylor (MDL): Created.
3562 : *
3563 : * NOTES
3564 : * 1) Have not implemented const 0 bit groups for prim miss or no miss.
3565 : * 2) MAX_GROUP_LEN (found experimentally) is useful to keep cost of groups
3566 : * down.
3567 : *****************************************************************************
3568 : */
3569 : #define MAX_GROUP_LEN 255
3570 0 : static void GroupIt (sInt4 OverallMin, sInt4 *Data, size_t numData,
3571 : TDLGroupType ** group, size_t *numGroup, char f_primMiss,
3572 : sInt4 li_primMiss, char f_secMiss,
3573 : sInt4 li_secMiss, sInt4 *groupSize, size_t *ibit,
3574 : size_t *jbit, size_t *kbit)
3575 : {
3576 : sInt4 A_max; /* Max value of a given group. */
3577 : sInt4 A_min; /* Min value of a given group. */
3578 : TDLGroupType G; /* Used to init the groups. */
3579 : int f_adjust; /* Flag if we have changed the groups. */
3580 : TDLGroupType *lclGroup; /* A temporary copy of the groups. */
3581 : int numLclGroup; /* # of groups in lclGroup. */
3582 : size_t xFactor; /* Estimate of cost (in bits) of a group. */
3583 : size_t i; /* loop counter. */
3584 :
3585 : /* Subtract the Overall Min Value. */
3586 0 : if (OverallMin != 0) {
3587 0 : if (f_secMiss) {
3588 0 : for (i = 0; i < numData; i++) {
3589 0 : if ((Data[i] != li_secMiss) && (Data[i] != li_primMiss)) {
3590 0 : Data[i] -= OverallMin;
3591 : /* Check if we accidently adjusted to prim or sec, if so add
3592 : * 1. */
3593 0 : if ((Data[i] == li_secMiss) || (Data[i] == li_primMiss)) {
3594 : myAssert (1 == 2);
3595 0 : Data[i]++;
3596 0 : if ((Data[i] == li_secMiss) || (Data[i] == li_primMiss)) {
3597 : myAssert (1 == 2);
3598 0 : Data[i]++;
3599 : }
3600 : }
3601 : }
3602 : }
3603 0 : } else if (f_primMiss) {
3604 0 : for (i = 0; i < numData; i++) {
3605 0 : if (Data[i] != li_primMiss) {
3606 0 : Data[i] -= OverallMin;
3607 : /* Check if we accidently adjusted to prim or sec, if so add
3608 : * 1. */
3609 0 : if (Data[i] == li_primMiss) {
3610 : myAssert (1 == 2);
3611 0 : Data[i]++;
3612 : }
3613 : }
3614 : }
3615 : } else {
3616 0 : for (i = 0; i < numData; i++) {
3617 0 : Data[i] -= OverallMin;
3618 : }
3619 : }
3620 : }
3621 :
3622 : myAssert ((f_secMiss == 0) || (f_secMiss == 1));
3623 : myAssert ((f_primMiss == 0) || (f_primMiss == 1));
3624 :
3625 : /* Create zero groups. */
3626 0 : *numGroup = 0;
3627 0 : *group = NULL;
3628 0 : if (f_primMiss) {
3629 0 : G.min = Data[0];
3630 0 : G.max = Data[0];
3631 0 : G.num = 1;
3632 0 : G.start = 0;
3633 0 : for (i = 1; i < numData; i++) {
3634 0 : if (G.min == li_primMiss) {
3635 0 : if (Data[i] == li_primMiss) {
3636 0 : G.num++;
3637 0 : if (G.num == (MAX_GROUP_LEN + 1)) {
3638 : /* Close a missing group */
3639 0 : G.f_trySplit = 0;
3640 0 : G.f_tryShift = 1;
3641 0 : G.bit = 0;
3642 0 : G.min = 0;
3643 0 : G.max = 0;
3644 0 : G.num = MAX_GROUP_LEN;
3645 0 : (*numGroup)++;
3646 : *group = (TDLGroupType *) realloc ((void *) *group,
3647 : *numGroup *
3648 0 : sizeof (TDLGroupType));
3649 0 : (*group)[(*numGroup) - 1] = G;
3650 : /* Init a missing group. */
3651 0 : G.min = Data[i];
3652 0 : G.max = Data[i];
3653 0 : G.num = 1;
3654 0 : G.start = i;
3655 : }
3656 : } else {
3657 : /* Close a missing group */
3658 0 : G.f_trySplit = 0;
3659 0 : G.f_tryShift = 1;
3660 0 : G.bit = 0;
3661 0 : G.min = 0;
3662 0 : G.max = 0;
3663 0 : (*numGroup)++;
3664 : *group = (TDLGroupType *) realloc ((void *) *group,
3665 : *numGroup *
3666 0 : sizeof (TDLGroupType));
3667 0 : (*group)[(*numGroup) - 1] = G;
3668 : /* Init a non-missing group. */
3669 0 : G.min = Data[i];
3670 0 : G.max = Data[i];
3671 0 : G.num = 1;
3672 0 : G.start = i;
3673 : }
3674 : } else {
3675 0 : if (Data[i] == li_primMiss) {
3676 : /* Close a non-missing group */
3677 0 : G.f_trySplit = 1;
3678 0 : G.f_tryShift = 1;
3679 : G.bit = (char) power ((uInt4) (G.max - G.min),
3680 0 : f_secMiss + f_primMiss);
3681 : myAssert (G.bit != 0);
3682 0 : if ((G.min == 0) && (G.bit == 0) && (f_primMiss == 1)) {
3683 : printf ("Warning: potential confusion between const value "
3684 0 : "and prim-missing.\n");
3685 0 : G.bit = 1;
3686 : }
3687 0 : (*numGroup)++;
3688 : *group = (TDLGroupType *) realloc ((void *) *group,
3689 : *numGroup *
3690 0 : sizeof (TDLGroupType));
3691 0 : (*group)[(*numGroup) - 1] = G;
3692 : /* Init a missing group. */
3693 0 : G.min = Data[i];
3694 0 : G.max = Data[i];
3695 0 : G.num = 1;
3696 0 : G.start = i;
3697 : } else {
3698 0 : if (G.min > Data[i]) {
3699 0 : G.min = Data[i];
3700 0 : } else if (G.max < Data[i]) {
3701 0 : G.max = Data[i];
3702 : }
3703 0 : G.num++;
3704 : }
3705 : }
3706 : }
3707 0 : if (G.min == li_primMiss) {
3708 : /* Close a missing group */
3709 0 : G.f_trySplit = 0;
3710 0 : G.f_tryShift = 1;
3711 0 : G.bit = 0;
3712 0 : G.min = 0;
3713 0 : G.max = 0;
3714 0 : (*numGroup)++;
3715 : *group = (TDLGroupType *) realloc ((void *) *group,
3716 : *numGroup *
3717 0 : sizeof (TDLGroupType));
3718 0 : (*group)[(*numGroup) - 1] = G;
3719 : } else {
3720 : /* Close a non-missing group */
3721 0 : G.f_trySplit = 1;
3722 0 : G.f_tryShift = 1;
3723 : G.bit = (char) power ((uInt4) (G.max - G.min),
3724 0 : f_secMiss + f_primMiss);
3725 : myAssert (G.bit != 0);
3726 0 : if ((G.min == 0) && (G.bit == 0) && (f_primMiss == 1)) {
3727 : printf ("Warning: potential confusion between const value and "
3728 0 : "prim-missing.\n");
3729 0 : G.bit = 1;
3730 : }
3731 0 : (*numGroup)++;
3732 : *group = (TDLGroupType *) realloc ((void *) *group,
3733 : *numGroup *
3734 0 : sizeof (TDLGroupType));
3735 0 : (*group)[(*numGroup) - 1] = G;
3736 : }
3737 : } else {
3738 : /* Already handled the f_primMiss case */
3739 0 : if (f_secMiss) {
3740 : findMaxMin2 (Data, 0, numData, li_primMiss, li_secMiss, &A_min,
3741 0 : &A_max);
3742 : } else {
3743 0 : findMaxMin0 (Data, 0, numData, &A_min, &A_max);
3744 : }
3745 0 : G.start = 0;
3746 0 : G.num = numData;
3747 0 : G.min = A_min;
3748 0 : G.max = A_max;
3749 0 : G.bit = (char) power ((uInt4) (A_max - A_min), f_secMiss + f_primMiss);
3750 0 : G.f_trySplit = 1;
3751 0 : G.f_tryShift = 1;
3752 0 : *numGroup = 1;
3753 0 : *group = (TDLGroupType *) malloc (sizeof (TDLGroupType));
3754 0 : (*group)[0] = G;
3755 : }
3756 :
3757 0 : lclGroup = NULL;
3758 0 : numLclGroup = 0;
3759 0 : *groupSize = ComputeGroupSize (*group, *numGroup, ibit, jbit, kbit);
3760 0 : xFactor = *ibit + *jbit + *kbit;
3761 : #ifdef DEBUG
3762 : /*
3763 : printf ("NumGroup = %d: Bytes = %ld: XFactor %d\n", *numGroup,
3764 : (*groupSize / 8) + 1, xFactor);
3765 : */
3766 : #endif
3767 :
3768 0 : f_adjust = 1;
3769 0 : while (f_adjust) {
3770 : f_adjust = splitGroup (Data, numData, *group, *numGroup, &lclGroup,
3771 : &numLclGroup, f_primMiss, li_primMiss,
3772 0 : f_secMiss, li_secMiss, xFactor);
3773 0 : free (*group);
3774 0 : *group = lclGroup;
3775 0 : *numGroup = numLclGroup;
3776 :
3777 0 : if (f_adjust) {
3778 : shiftGroup (Data, numData, group, numGroup, f_primMiss,
3779 0 : li_primMiss, f_secMiss, li_secMiss, xFactor);
3780 0 : *groupSize = ComputeGroupSize (*group, *numGroup, ibit, jbit, kbit);
3781 0 : if (xFactor != *ibit + *jbit + *kbit) {
3782 0 : for (i = 0; i < *numGroup; i++) {
3783 0 : if (((*group)[i].num > *ibit + *jbit + *kbit) &&
3784 0 : ((*group)[i].num <= xFactor)) {
3785 0 : (*group)[i].f_trySplit = 1;
3786 : }
3787 : }
3788 : }
3789 0 : xFactor = *ibit + *jbit + *kbit;
3790 : #ifdef DEBUG
3791 : /*
3792 : printf ("NumGroup = %d: Bytes = %ld: XFactor %d\n", *numGroup,
3793 : (*groupSize / 8) + 1, xFactor);
3794 : fflush (stdout);
3795 : */
3796 : #endif
3797 : }
3798 : }
3799 0 : }
3800 :
3801 : /*****************************************************************************
3802 : * GroupPack() --
3803 : *
3804 : * Arthur Taylor / MDL
3805 : *
3806 : * PURPOSE
3807 : * To compute groups for packing the data using complex or second order
3808 : * complex packing.
3809 : *
3810 : * ARGUMENTS
3811 : * Src = The original data. (Input)
3812 : * Dst = The scaled data. (Output)
3813 : * numData = The number of elements in data. (Input)
3814 : * DSF = Decimal Scale Factor for scaling the data. (Input)
3815 : * BSF = Binary Scale Factor for scaling the data. (Input)
3816 : * f_primMiss = Flag saying if we have a primary missing value (In/Out)
3817 : * primMiss = primary missing value. (In/Out)
3818 : * f_secMiss = Flag saying if we have a secondary missing value (In/Out)
3819 : * secMiss = secondary missing value. (In/Out)
3820 : * f_grid = Flag if this is grid data (or vector) (Input)
3821 : * NX = The number of X values. (Input)
3822 : * NY = The number of Y values. (Input)
3823 : * f_sndOrder = Flag if we should do second order diffencing (Output)
3824 : * group = Resulting groups. (Output)
3825 : * numGroup = Number of groups. (Output)
3826 : * Min = Overall minimum. (Output)
3827 : * a1 = if f_sndOrder, the first first order difference (Output)
3828 : * b2 = if f_sndOrder, the first second order difference (Output)
3829 : * groupSize = How many bytes the groups and data will take. (Output)
3830 : * ibit = # of bits for largest minimum value in groups (Output)
3831 : * jbit = # of bits for largest # bits in groups (Output)
3832 : * kbit = # of bits for largest # values in groups (Output)
3833 : *
3834 : * FILES/DATABASES: None
3835 : *
3836 : * RETURNS: int (could use errSprintf())
3837 : * 0 = OK
3838 : * -1 = Primary or Secondary missing value == 0.
3839 : *
3840 : * HISTORY
3841 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
3842 : * 1/2005 AAT: Cleaned up.
3843 : *
3844 : * NOTES
3845 : *****************************************************************************
3846 : */
3847 0 : static int GroupPack (double *Src, sInt4 **Dst, sInt4 numData,
3848 : int DSF, int BSF, char *f_primMiss, double *primMiss,
3849 : char *f_secMiss, double *secMiss, char f_grid,
3850 : short int NX, short int NY, char *f_sndOrder,
3851 : TDLGroupType ** group, size_t *numGroup,
3852 : sInt4 *Min, sInt4 *a1, sInt4 *b2,
3853 : sInt4 *groupSize, size_t *ibit, size_t *jbit,
3854 : size_t *kbit)
3855 : {
3856 0 : sInt4 *SecDiff = NULL; /* Consists of the 2nd order differences if *
3857 : * requested. */
3858 : sInt4 *Data; /* The scaled data. */
3859 : char f_min; /* Flag saying overallMin is valid. */
3860 : sInt4 overallMin; /* The overall min of the scaled data. */
3861 : sInt4 secMin; /* The overall min of the 2nd order differences */
3862 0 : sInt4 li_primMiss = 0; /* The scaled primary missing value */
3863 0 : sInt4 li_secMiss = 0; /* The scaled secondary missing value */
3864 0 : int minGroup = 20; /* The minimum group size. Equivalent to xFactor?
3865 : * Chose 20 because that was a good estimate of
3866 : * XFactor. */
3867 :
3868 : /* Check consistency of f_primMiss and f_secMiss. */
3869 0 : if (*primMiss == *secMiss) {
3870 0 : *f_secMiss = 0;
3871 : }
3872 0 : if ((*f_secMiss) && (!(*f_primMiss))) {
3873 0 : *f_primMiss = *f_secMiss;
3874 0 : *primMiss = *secMiss;
3875 0 : *f_secMiss = 0;
3876 : }
3877 0 : if (*f_secMiss && (*secMiss == 0)) {
3878 0 : errSprintf ("Error: Secondary missing value not allowed to = 0.\n");
3879 0 : return -1;
3880 : }
3881 0 : if (*f_primMiss && (*primMiss == 0)) {
3882 0 : errSprintf ("Error: Primary missing value not allowed to = 0.\n");
3883 0 : return -1;
3884 : }
3885 :
3886 : /* Check minGroup size. */
3887 0 : if (minGroup > numData) {
3888 0 : minGroup = numData;
3889 : }
3890 :
3891 : /* Scale the data and check if we can change f_prim or f_sec. */
3892 : /* Note: if we use sec_diff, we have a different overall min. */
3893 0 : f_min = 0;
3894 0 : Data = (sInt4 *) malloc (numData * sizeof (sInt4));
3895 : TDL_ScaleData (Src, Data, numData, DSF, BSF, f_primMiss, primMiss,
3896 0 : f_secMiss, secMiss, &f_min, &overallMin);
3897 : /* Note: ScaleData also scales missing values. */
3898 0 : if (*f_primMiss) {
3899 0 : li_primMiss = (sInt4) (*primMiss * SCALE_MISSING + .5);
3900 : }
3901 0 : if (*f_secMiss) {
3902 0 : li_secMiss = (sInt4) (*secMiss * SCALE_MISSING + .5);
3903 : }
3904 :
3905 : /* Reason this is after TDL_ScaleData is we don't want to reoder the
3906 : * caller's copy of the data. */
3907 0 : if (f_grid) {
3908 0 : TDL_ReorderGrid (Data, NX, NY);
3909 : } else {
3910 : /* TDLPack requires the following (see pack2d.f and pack1d.f) */
3911 0 : *f_sndOrder = 0;
3912 : }
3913 : /* TDLPack requires the following (see pack2d.f) */
3914 0 : if (*f_secMiss) {
3915 0 : *f_sndOrder = 0;
3916 : }
3917 : /* If overallMin is "invalid" then they are all prim or sec values. */
3918 : /* See pack2d.f line 336 */
3919 : /* IF ALL VALUES ARE MISSING, JUST PACK; DON'T CONSIDER 2ND ORDER
3920 : * DIFFERENCES. */
3921 0 : if (!f_min) {
3922 0 : *f_sndOrder = 0;
3923 : }
3924 :
3925 : /* This has to be after TDL_ReorderGrid */
3926 0 : if (*f_sndOrder) {
3927 0 : SecDiff = (sInt4 *) malloc (numData * sizeof (sInt4));
3928 0 : if (TDL_GetSecDiff (Data, numData, SecDiff, *f_primMiss, li_primMiss,
3929 : a1, b2, &secMin)) {
3930 : /* Problem finding SecDiff, so we don't bother with it. */
3931 0 : *f_sndOrder = 0;
3932 : } else {
3933 : /* Check if it is worth doing second order differences. */
3934 0 : if (*f_primMiss) {
3935 : *f_sndOrder = TDL_UseSecDiff_Prim (Data, numData, SecDiff,
3936 0 : li_primMiss, minGroup);
3937 : } else {
3938 0 : *f_sndOrder = TDL_UseSecDiff (Data, numData, SecDiff, minGroup);
3939 : }
3940 : }
3941 : }
3942 :
3943 : /* Side affect of GroupIt2: it subtracts OverallMin from Data. */
3944 0 : if (!(*f_sndOrder)) {
3945 : GroupIt (overallMin, Data, numData, group, numGroup, *f_primMiss,
3946 : li_primMiss, *f_secMiss, li_secMiss, groupSize, ibit, jbit,
3947 0 : kbit);
3948 0 : *Min = overallMin;
3949 0 : *a1 = 0;
3950 0 : *b2 = 0;
3951 0 : *Dst = Data;
3952 0 : free (SecDiff);
3953 : } else {
3954 : GroupIt (secMin, SecDiff, numData, group, numGroup, *f_primMiss,
3955 : li_primMiss, *f_secMiss, li_secMiss, groupSize, ibit, jbit,
3956 0 : kbit);
3957 0 : *Min = secMin;
3958 0 : *Dst = SecDiff;
3959 0 : free (Data);
3960 : }
3961 0 : return 0;
3962 : }
3963 :
3964 : /*****************************************************************************
3965 : * WriteTDLPRecord() --
3966 : *
3967 : * Arthur Taylor / MDL
3968 : *
3969 : * PURPOSE
3970 : * Writes a TDLP message to file.
3971 : *
3972 : * ARGUMENTS
3973 : * fp = An opened TDLP file already at the correct location. (Input)
3974 : * Data = The data to write. (Input)
3975 : * DataLen = Length of Data. (Input)
3976 : * DSF = Decimal scale factor to apply to the data (Input)
3977 : * BSF = Binary scale factor to apply to the data (Input)
3978 : * f_primMiss = Flag saying if we have a primary missing value (Input)
3979 : * primMiss = primary missing value. (Input)
3980 : * f_secMiss = Flag saying if we have a secondary missing value (Input)
3981 : * secMiss = secondary missing value. (Input)
3982 : * gds = The grid definition section (Input)
3983 : * comment = Describes the kind of data (max 32 bytes). (Input)
3984 : * refTime = The reference (creation) time of this message. (Input)
3985 : * ID1 = TDLPack ID1 (Input)
3986 : * ID2 = TDLPack ID2 (Input)
3987 : * ID3 = TDLPack ID3 (Input)
3988 : * ID4 = TDLPack ID4 (Input)
3989 : * projSec = The projection in seconds (Input)
3990 : * processNum = The process number that created it (Input)
3991 : * seqNum = The sequence number that created it (Input)
3992 : *
3993 : * FILES/DATABASES:
3994 : * An already opened file pointing to the desired TDLP message.
3995 : *
3996 : * RETURNS: int (could use errSprintf())
3997 : * 0 = OK
3998 : * -1 = comment is too long.
3999 : * -2 = projHr is inconsistent with ID3.
4000 : * -3 = Type of map projection that TDLP can't handle.
4001 : * -4 = Primary or Secondary missing value == 0.
4002 : *
4003 : * HISTORY
4004 : * 12/2004 Arthur Taylor (MDL): Created
4005 : * 1/2005 AAT: Cleaned up.
4006 : *
4007 : * NOTES
4008 : *****************************************************************************
4009 : */
4010 0 : int WriteTDLPRecord (FILE * fp, double *Data, sInt4 DataLen, int DSF,
4011 : int BSF, char f_primMiss, double primMiss,
4012 : char f_secMiss, double secMiss, gdsType *gds,
4013 : char *comment, double refTime, sInt4 ID1,
4014 : sInt4 ID2, sInt4 ID3, sInt4 ID4,
4015 : sInt4 projSec, sInt4 processNum, sInt4 seqNum)
4016 : {
4017 : sInt4 *Scaled; /* The scaled data. */
4018 : TDLGroupType *group; /* The groups used to pack the data. */
4019 : size_t numGroup; /* Number of groups. */
4020 0 : char f_grid = 1; /* Flag if this is gridded data. In theory can handle
4021 : * vector data, but haven't tested it. */
4022 : char f_sndOrder; /* Flag if we should try second order packing. */
4023 : sInt4 overallMin; /* Overall min value of the scaled data. */
4024 : sInt4 a1; /* 2nd order difference: 1st value. */
4025 : sInt4 b2; /* 2nd order difference: 1st 1st order difference */
4026 : sInt4 li_primMiss; /* Scaled primary missing value. */
4027 : sInt4 li_secMiss; /* Scaled secondary missing value. */
4028 : int mbit; /* # of bits for b2. */
4029 : int nbit; /* # of bits for overallMin. */
4030 : size_t ibit; /* # of bits for largest minimum value in groups */
4031 : size_t jbit; /* # of bits for largest # bits in groups */
4032 : size_t kbit; /* # of bits for largest # values in groups */
4033 : int sec1Len; /* Length of section 1. */
4034 : size_t pad; /* Number of bytes to pad the message to get to the
4035 : * correct byte boundary. */
4036 : sInt4 groupSize; /* How many bytes the groups and data will take. */
4037 : sInt4 sec4Len; /* Length of section 4. */
4038 : sInt4 tdlRecSize; /* Size of the TDLP message. */
4039 : sInt4 recSize; /* Actual record size (including FORTRAN bytes). */
4040 : int commentLen; /* Length of comment */
4041 : short int projHr; /* The hours part of the forecast projection. */
4042 : char projMin; /* The minutes part of the forecast projection. */
4043 : sInt4 year; /* The reference year. */
4044 : int month, day; /* The reference month day. */
4045 : int hour, min; /* The reference hour minute. */
4046 : double sec; /* The reference second. */
4047 0 : char f_bitmap = 0; /* Bitmap flag: not implemented in specs. */
4048 0 : char f_simple = 0; /* Simple Pack flag: not implemented in specs. */
4049 : int gridType; /* Which type of grid. (Polar, Mercator, Lambert). */
4050 : int dataCnt; /* Keeps track of which element we are writing. */
4051 : sInt4 max0; /* The max value in a group. Represents primary or *
4052 : * secondary missing value depending on scheme. */
4053 : sInt4 max1; /* The next to max value in a group. Represents *
4054 : * secondary missing value. */
4055 : size_t i, j; /* loop counters */
4056 : sInt4 li_temp; /* Temporary variable (sInt4). */
4057 : short int si_temp; /* Temporary variable (short int). */
4058 : double d_temp; /* Temporary variable (double). */
4059 : char buffer[6]; /* Used to write reserved values */
4060 : uChar pbuf; /* A buffer of bits that weren't written to disk */
4061 : sChar pbufLoc; /* Where in pbuf to add more bits. */
4062 :
4063 0 : commentLen = strlen (comment);
4064 0 : if (commentLen > 32) {
4065 0 : errSprintf ("Error: '%s' is > 32 bytes long\n", comment);
4066 0 : return -1;
4067 : }
4068 0 : projHr = projSec / 3600;
4069 0 : projMin = (projSec % 3600) / 60;
4070 0 : if (projHr != (ID3 - ((ID3 / 1000) * 1000))) {
4071 : errSprintf ("Error: projHr = %d is inconsistent with ID3 = %ld\n",
4072 0 : projHr, ID3);
4073 0 : return -2;
4074 : }
4075 0 : if (f_grid) {
4076 0 : switch (gds->projType) {
4077 : case GS3_POLAR:
4078 0 : gridType = TDLP_POLAR;
4079 0 : break;
4080 : case GS3_LAMBERT:
4081 0 : gridType = TDLP_LAMBERT;
4082 0 : break;
4083 : case GS3_MERCATOR:
4084 0 : gridType = TDLP_MERCATOR;
4085 0 : break;
4086 : default:
4087 : errSprintf ("TDLPack can't handle GRIB projection type %d\n",
4088 0 : gds->projType);
4089 0 : return -3;
4090 : }
4091 : }
4092 :
4093 0 : if (GroupPack (Data, &Scaled, DataLen, DSF, BSF, &f_primMiss, &primMiss,
4094 : &f_secMiss, &secMiss, f_grid, gds->Nx, gds->Ny,
4095 : &f_sndOrder, &group, &numGroup, &overallMin, &a1, &b2,
4096 : &groupSize, &ibit, &jbit, &kbit) != 0) {
4097 0 : return -4;
4098 : }
4099 :
4100 : /* Make sure missing data is properly scaled. */
4101 0 : if (f_primMiss) {
4102 0 : li_primMiss = (sInt4) (primMiss * SCALE_MISSING + .5);
4103 : }
4104 0 : if (f_secMiss) {
4105 0 : li_secMiss = (sInt4) (secMiss * SCALE_MISSING + .5);
4106 : }
4107 :
4108 : /* Compute TDL record size. */
4109 : /* *INDENT-OFF* */
4110 : /* TDL Record size
4111 : * 8 (section 0),
4112 : * 39 + strlen(comment) (section 1),
4113 : * 0 or 28 (depending on if you have a section 2),
4114 : * 0 (section 3),
4115 : * 16 + 5 + 1 + nbit(min val) + 16 + 5 + 5 + 5 + GroupSize() (section 4)
4116 : * 4 (section 5)
4117 : * pad (to even 8 bytes...)
4118 : * Group size uses:
4119 : * ibit * num_groups + jbit * num_groups +
4120 : * kbit * num_groups + group.nbit * group.nvalues */
4121 : /* *INDENT-ON* */
4122 0 : sec1Len = 39 + commentLen;
4123 0 : if (overallMin < 0) {
4124 0 : nbit = power (-1 * overallMin, 0);
4125 : } else {
4126 0 : nbit = power (overallMin, 0);
4127 : }
4128 0 : if (!f_sndOrder) {
4129 0 : if (f_secMiss) {
4130 : sec4Len = 16 + (sInt4) (ceil ((5 + 1 + nbit + 16 + 5 + 5 + 5
4131 0 : + groupSize) / 8.));
4132 0 : } else if (f_primMiss) {
4133 : sec4Len = 12 + (sInt4) (ceil ((5 + 1 + nbit + 16 + 5 + 5 + 5
4134 0 : + groupSize) / 8.));
4135 : } else {
4136 : sec4Len = 8 + (sInt4) (ceil ((5 + 1 + nbit + 16 + 5 + 5 + 5
4137 0 : + groupSize) / 8.));
4138 : }
4139 : } else {
4140 0 : if (b2 < 0) {
4141 0 : mbit = power (-1 * b2, 0);
4142 : } else {
4143 0 : mbit = power (b2, 0);
4144 : }
4145 0 : if (f_secMiss) {
4146 : sec4Len =
4147 : 16 +
4148 : (sInt4) (ceil
4149 : ((32 + 5 + 1 + mbit + 5 + 1 + nbit + 16 + 5 + 5 + 5 +
4150 0 : groupSize) / 8.));
4151 0 : } else if (f_primMiss) {
4152 : sec4Len =
4153 : 12 +
4154 : (sInt4) (ceil
4155 : ((32 + 5 + 1 + mbit + 5 + 1 + nbit + 16 + 5 + 5 + 5 +
4156 0 : groupSize) / 8.));
4157 : } else {
4158 : sec4Len =
4159 : 8 +
4160 : (sInt4) (ceil
4161 : ((32 + 5 + 1 + mbit + 5 + 1 + nbit + 16 + 5 + 5 + 5 +
4162 0 : groupSize) / 8.));
4163 : }
4164 : }
4165 0 : if (f_grid) {
4166 0 : tdlRecSize = 8 + sec1Len + 28 + 0 + sec4Len + 4;
4167 : } else {
4168 0 : tdlRecSize = 8 + sec1Len + 0 + 0 + sec4Len + 4;
4169 : }
4170 : /* Actual recSize is 8 + round(tdlRecSize) to nearest 8 byte boundary. */
4171 0 : recSize = 8 + (sInt4) (ceil (tdlRecSize / 8.0)) * 8;
4172 0 : pad = (int) ((recSize - 8) - tdlRecSize);
4173 :
4174 : /* --- Start Writing record. --- */
4175 : /* First write FORTRAN record information */
4176 0 : fread(&(recSize), sizeof (sInt4), 1, fp);
4177 0 : li_temp = 0;
4178 0 : FWRITE_BIG (&(li_temp), sizeof (sInt4), 1, fp);
4179 0 : li_temp = recSize - 8; /* FORTRAN rec length. */
4180 0 : FWRITE_BIG (&(li_temp), sizeof (sInt4), 1, fp);
4181 :
4182 : /* Now write section 0. */
4183 0 : fwrite ("TDLP", sizeof (char), 4, fp);
4184 0 : FWRITE_ODDINT_BIG (&(tdlRecSize), 3, fp);
4185 : /* version number */
4186 0 : fputc (0, fp);
4187 :
4188 : /* Now write section 1. */
4189 0 : fputc (sec1Len, fp);
4190 : /* output type specification... */
4191 0 : i = 0;
4192 0 : if (f_grid) {
4193 0 : i |= 1;
4194 : }
4195 0 : if (f_bitmap) {
4196 0 : i |= 2;
4197 : }
4198 0 : fputc (i, fp);
4199 : /* tempTime = gmtime (&(refTime));*/
4200 0 : Clock_PrintDate (refTime, &year, &month, &day, &hour, &min, &sec);
4201 : /* year = tempTime->tm_year + 1900;
4202 : month = tempTime->tm_mon + 1;
4203 : day = tempTime->tm_mday;
4204 : hour = tempTime->tm_hour;
4205 : min = tempTime->tm_min;
4206 : */
4207 0 : si_temp = year;
4208 0 : FWRITE_BIG (&si_temp, sizeof (short int), 1, fp);
4209 0 : fputc (month, fp);
4210 0 : fputc (day, fp);
4211 0 : fputc (hour, fp);
4212 0 : fputc (min, fp);
4213 0 : li_temp = (year * 1000000L + month * 10000L + day * 100 + hour);
4214 0 : FWRITE_BIG (&li_temp, sizeof (sInt4), 1, fp);
4215 0 : FWRITE_BIG (&ID1, sizeof (sInt4), 1, fp);
4216 0 : FWRITE_BIG (&ID2, sizeof (sInt4), 1, fp);
4217 0 : FWRITE_BIG (&ID3, sizeof (sInt4), 1, fp);
4218 0 : FWRITE_BIG (&ID4, sizeof (sInt4), 1, fp);
4219 0 : FWRITE_BIG (&projHr, sizeof (short int), 1, fp);
4220 0 : fputc (projMin, fp);
4221 0 : fputc (processNum, fp);
4222 0 : fputc (seqNum, fp);
4223 0 : i = (DSF < 0) ? 128 - DSF : DSF;
4224 0 : fputc (i, fp);
4225 0 : i = (BSF < 0) ? 128 - BSF : BSF;
4226 0 : fputc (i, fp);
4227 : /* Reserved: 3 bytes of 0. */
4228 0 : li_temp = 0;
4229 0 : fwrite (&li_temp, sizeof (char), 3, fp);
4230 0 : fputc (commentLen, fp);
4231 0 : fwrite (comment, sizeof (char), commentLen, fp);
4232 :
4233 : /* Now write section 2. */
4234 0 : if (f_grid) {
4235 0 : fputc (28, fp);
4236 0 : fputc (gridType, fp);
4237 0 : si_temp = gds->Nx;
4238 0 : FWRITE_BIG (&si_temp, sizeof (short int), 1, fp);
4239 0 : si_temp = gds->Ny;
4240 0 : FWRITE_BIG (&si_temp, sizeof (short int), 1, fp);
4241 0 : li_temp = (sInt4) (gds->lat1 * 10000. + .5);
4242 0 : if (li_temp < 0) {
4243 0 : pbuf = 128;
4244 0 : li_temp = -1 * li_temp;
4245 : } else {
4246 0 : pbuf = 0;
4247 : }
4248 0 : pbufLoc = 7;
4249 0 : fileBitWrite (&(li_temp), sizeof (li_temp), 23, fp, &pbuf, &pbufLoc);
4250 : myAssert (pbufLoc == 8);
4251 : myAssert (pbuf == 0);
4252 0 : d_temp = 360 - gds->lon1;
4253 0 : if (d_temp < 0)
4254 0 : d_temp += 360;
4255 0 : if (d_temp > 360)
4256 0 : d_temp -= 360;
4257 0 : li_temp = (sInt4) (d_temp * 10000. + .5);
4258 0 : if (li_temp < 0) {
4259 0 : pbuf = 128;
4260 0 : li_temp = -1 * li_temp;
4261 : }
4262 0 : pbufLoc = 7;
4263 0 : fileBitWrite (&(li_temp), sizeof (li_temp), 23, fp, &pbuf, &pbufLoc);
4264 : myAssert (pbufLoc == 8);
4265 : myAssert (pbuf == 0);
4266 0 : d_temp = 360 - gds->orientLon;
4267 0 : if (d_temp < 0)
4268 0 : d_temp += 360;
4269 0 : if (d_temp > 360)
4270 0 : d_temp -= 360;
4271 0 : li_temp = (sInt4) (d_temp * 10000. + .5);
4272 0 : if (li_temp < 0) {
4273 0 : pbuf = 128;
4274 0 : li_temp = -1 * li_temp;
4275 : }
4276 0 : pbufLoc = 7;
4277 0 : fileBitWrite (&(li_temp), sizeof (li_temp), 23, fp, &pbuf, &pbufLoc);
4278 : myAssert (pbufLoc == 8);
4279 : myAssert (pbuf == 0);
4280 0 : li_temp = (sInt4) (gds->Dx * 1000. + .5);
4281 0 : FWRITE_BIG (&li_temp, sizeof (sInt4), 1, fp);
4282 0 : li_temp = (sInt4) (gds->meshLat * 10000. + .5);
4283 0 : if (li_temp < 0) {
4284 0 : pbuf = 128;
4285 0 : li_temp = -1 * li_temp;
4286 : }
4287 0 : pbufLoc = 7;
4288 0 : fileBitWrite (&(li_temp), sizeof (li_temp), 23, fp, &pbuf, &pbufLoc);
4289 : myAssert (pbufLoc == 8);
4290 : myAssert (pbuf == 0);
4291 0 : memset (buffer, 0, 6);
4292 0 : fwrite (buffer, sizeof (char), 6, fp);
4293 : }
4294 :
4295 : /* Now write section 3. */
4296 : /* Bitmap is not supported, skipping. */
4297 : myAssert (!f_bitmap);
4298 :
4299 : /* Now write section 4. */
4300 0 : FWRITE_ODDINT_BIG (&(sec4Len), 3, fp);
4301 0 : i = 0;
4302 0 : if (f_secMiss)
4303 0 : i |= 1;
4304 0 : if (f_primMiss)
4305 0 : i |= 2;
4306 0 : if (f_sndOrder)
4307 0 : i |= 4;
4308 0 : if (!f_simple)
4309 0 : i |= 8;
4310 0 : if (!f_grid)
4311 0 : i |= 16;
4312 0 : fputc (i, fp);
4313 0 : li_temp = DataLen;
4314 0 : FWRITE_BIG (&li_temp, sizeof (sInt4), 1, fp);
4315 0 : if (f_primMiss) {
4316 0 : FWRITE_BIG (&(li_primMiss), sizeof (sInt4), 1, fp);
4317 0 : if (f_secMiss) {
4318 0 : FWRITE_BIG (&(li_secMiss), sizeof (sInt4), 1, fp);
4319 : }
4320 : }
4321 0 : if (f_sndOrder) {
4322 0 : if (a1 < 0) {
4323 0 : pbuf = 128;
4324 0 : li_temp = -1 * a1;
4325 : } else {
4326 0 : pbuf = 0;
4327 0 : li_temp = a1;
4328 : }
4329 0 : pbufLoc = 7;
4330 0 : fileBitWrite (&(li_temp), sizeof (li_temp), 31, fp, &pbuf, &pbufLoc);
4331 : myAssert (pbufLoc == 8);
4332 : myAssert (pbuf == 0);
4333 0 : fileBitWrite (&mbit, sizeof (mbit), 5, fp, &pbuf, &pbufLoc);
4334 0 : if (b2 < 0) {
4335 0 : i = 1;
4336 0 : li_temp = -1 * b2;
4337 : } else {
4338 0 : i = 0;
4339 0 : li_temp = b2;
4340 : }
4341 0 : fileBitWrite (&i, sizeof (i), 1, fp, &pbuf, &pbufLoc);
4342 : myAssert (pbufLoc == 2);
4343 : fileBitWrite (&li_temp, sizeof (li_temp), (unsigned short int) mbit,
4344 0 : fp, &pbuf, &pbufLoc);
4345 : }
4346 0 : fileBitWrite (&nbit, sizeof (nbit), 5, fp, &pbuf, &pbufLoc);
4347 0 : if (overallMin < 0) {
4348 0 : i = 1;
4349 0 : li_temp = -1 * overallMin;
4350 : } else {
4351 0 : i = 0;
4352 0 : li_temp = overallMin;
4353 : }
4354 0 : fileBitWrite (&i, sizeof (i), 1, fp, &pbuf, &pbufLoc);
4355 : fileBitWrite (&li_temp, sizeof (li_temp), (unsigned short int) nbit, fp,
4356 0 : &pbuf, &pbufLoc);
4357 0 : fileBitWrite (&numGroup, sizeof (numGroup), 16, fp, &pbuf, &pbufLoc);
4358 0 : fileBitWrite (&ibit, sizeof (ibit), 5, fp, &pbuf, &pbufLoc);
4359 0 : fileBitWrite (&jbit, sizeof (jbit), 5, fp, &pbuf, &pbufLoc);
4360 0 : fileBitWrite (&kbit, sizeof (kbit), 5, fp, &pbuf, &pbufLoc);
4361 0 : for (i = 0; i < numGroup; i++) {
4362 0 : fileBitWrite (&(group[i].min), sizeof (sInt4),
4363 0 : (unsigned short int) ibit, fp, &pbuf, &pbufLoc);
4364 : }
4365 0 : for (i = 0; i < numGroup; i++) {
4366 0 : fileBitWrite (&(group[i].bit), sizeof (char),
4367 0 : (unsigned short int) jbit, fp, &pbuf, &pbufLoc);
4368 : }
4369 : #ifdef DEBUG
4370 : li_temp = 0;
4371 : #endif
4372 0 : for (i = 0; i < numGroup; i++) {
4373 0 : fileBitWrite (&(group[i].num), sizeof (sInt4),
4374 0 : (unsigned short int) kbit, fp, &pbuf, &pbufLoc);
4375 : #ifdef DEBUG
4376 : li_temp += group[i].num;
4377 : #endif
4378 : }
4379 : #ifdef DEBUG
4380 : /* Sanity check ! */
4381 : if (li_temp != DataLen) {
4382 : printf ("Total packed in groups %d != DataLen %d\n", li_temp,
4383 : DataLen);
4384 : }
4385 : myAssert (li_temp == DataLen);
4386 : #endif
4387 :
4388 0 : dataCnt = 0;
4389 : /* Start going through the data. Grid data has already been reordered.
4390 : * Already taken care of 2nd order differences. Data at this point should
4391 : * contain a stream of either 2nd order differences or 0 order
4392 : * differences. We have also removed overallMin from all non-missing
4393 : * elements in the stream */
4394 0 : if (f_secMiss) {
4395 0 : for (i = 0; i < numGroup; i++) {
4396 0 : max0 = (1 << group[i].bit) - 1;
4397 0 : max1 = (1 << group[i].bit) - 2;
4398 0 : for (j = 0; j < group[i].num; j++) {
4399 0 : if (Scaled[dataCnt] == li_primMiss) {
4400 0 : li_temp = max0;
4401 0 : } else if (Scaled[dataCnt] == li_secMiss) {
4402 0 : li_temp = max1;
4403 : } else {
4404 0 : li_temp = Scaled[dataCnt] - group[i].min;
4405 : }
4406 0 : fileBitWrite (&(li_temp), sizeof (sInt4), group[i].bit, fp,
4407 0 : &pbuf, &pbufLoc);
4408 0 : dataCnt++;
4409 : }
4410 : }
4411 0 : } else if (f_primMiss) {
4412 0 : for (i = 0; i < numGroup; i++) {
4413 : /* see what happens when bit == 0. */
4414 0 : max0 = (1 << group[i].bit) - 1;
4415 0 : for (j = 0; j < group[i].num; j++) {
4416 0 : if (group[i].bit != 0) {
4417 0 : if (Scaled[dataCnt] == li_primMiss) {
4418 0 : li_temp = max0;
4419 : } else {
4420 0 : li_temp = Scaled[dataCnt] - group[i].min;
4421 : }
4422 0 : fileBitWrite (&(li_temp), sizeof (sInt4), group[i].bit, fp,
4423 0 : &pbuf, &pbufLoc);
4424 : }
4425 0 : dataCnt++;
4426 : }
4427 : }
4428 : } else {
4429 0 : for (i = 0; i < numGroup; i++) {
4430 0 : for (j = 0; j < group[i].num; j++) {
4431 0 : li_temp = Scaled[dataCnt] - group[i].min;
4432 0 : if (group[i].bit != 0) {
4433 0 : fileBitWrite (&(li_temp), sizeof (sInt4), group[i].bit, fp,
4434 0 : &pbuf, &pbufLoc);
4435 : }
4436 0 : dataCnt++;
4437 : }
4438 : }
4439 : }
4440 : myAssert (dataCnt == DataLen);
4441 : /* flush the PutBit buffer... */
4442 0 : if (pbufLoc != 8) {
4443 0 : fputc ((int) pbuf, fp);
4444 : }
4445 :
4446 : /* Now write section 5. */
4447 0 : fwrite ("7777", sizeof (char), 4, fp);
4448 :
4449 : /* Deal with padding at end of record... */
4450 0 : for (i = 0; i < pad; i++) {
4451 0 : fputc (0, fp);
4452 : }
4453 : /* Now write FORTRAN record information */
4454 0 : FWRITE_BIG (&(recSize), sizeof (sInt4), 1, fp);
4455 :
4456 0 : free (Scaled);
4457 0 : free (group);
4458 0 : return 0;
4459 : }
|
