source: trunk/postgis/geography_estimate.c @ 5625

Last change on this file since 5625 was 5625, checked in by pramsey, 6 years ago

Remove crash for case when all geographies are on the outer edges of the histobox, causing all to be classified as "deviants" by the stdev code line (#474)

File size: 44.4 KB
Line 
1/**********************************************************************
2 * $Id: geography_estimate.c 4211 2009-06-25 03:32:43Z robe $
3 *
4 * PostGIS - Spatial Types for PostgreSQL
5 * http://postgis.refractions.net
6 * Copyright 2001-2006 Refractions Research Inc.
7 * Copyright 2009 Mark Cave-Ayland
8 *
9 * This is free software; you can redistribute and/or modify it under
10 * the terms of the GNU General Public Licence. See the COPYING file.
11 *
12 **********************************************************************/
13
14
15#include "postgres.h"
16#include "commands/vacuum.h"
17#include "nodes/relation.h"
18#include "parser/parsetree.h"
19#include "utils/lsyscache.h"
20#include "utils/syscache.h"
21
22#include "libgeom.h"
23#include "lwgeom_pg.h"
24
25/* Prototypes */
26Datum geography_gist_selectivity(PG_FUNCTION_ARGS);
27Datum geography_gist_join_selectivity(PG_FUNCTION_ARGS);
28Datum geography_analyze(PG_FUNCTION_ARGS);
29
30
31/**
32* Place holder selectivity calculations to make the index work in
33* the absence of real selectivity calculations.
34*/
35
36#define DEFAULT_GEOGRAPHY_SEL 0.000005
37
38/**
39 *      Assign a number to the postgis statistics kind
40 *
41 *      tgl suggested:
42 *
43 *      1-100:  reserved for assignment by the core Postgres project
44 *      100-199: reserved for assignment by PostGIS
45 *      200-9999: reserved for other globally-known stats kinds
46 *      10000-32767: reserved for private site-local use
47 *
48 */
49#define STATISTIC_KIND_GEOGRAPHY 101
50
51/*
52 * Define this if you want to use standard deviation based
53 * histogram extent computation. If you do, you can also
54 * tweak the deviation factor used in computation with
55 * SDFACTOR.
56 */
57#define USE_STANDARD_DEVIATION 1
58#define SDFACTOR 3.25
59
60
61/* Information about the dimensions stored in the sample */
62struct dimensions
63{
64        char axis;
65        double min;
66        double max;
67};
68
69/* Main geography statistics structure, stored in pg_statistics */
70typedef struct GEOG_STATS_T
71{
72        /* Dimensionality of this column */
73        float4 dims;
74
75        /* x . y . z = total boxes in grid */
76        float4 unitsx;
77        float4 unitsy;
78        float4 unitsz;
79
80        /* average feature coverage of not-null features:
81           this is either an area for the case of a 2D column
82           or a volume in the case of a 3D column */
83        float4 avgFeatureCoverage;
84
85        /*
86         * average number of histogram cells
87         * covered by the sample not-null features
88         */
89        float4 avgFeatureCells;
90
91        /* histogram extent */
92        float4 xmin, ymin, zmin, xmax, ymax, zmax;
93
94        /* No. of geometries within this column (approx) */
95        float4 totalrows;
96
97        /*
98         * variable length # of floats for histogram
99         */
100        float4 value[1];
101}
102GEOG_STATS;
103
104
105/**
106 * This function returns an estimate of the selectivity
107 * of a search_box looking at data in the GEOG_STATS
108 * structure.
109 * */
110static float8
111estimate_selectivity(GBOX *box, GEOG_STATS *geogstats)
112{
113        int x, y, z;
114        int x_idx_min, x_idx_max, y_idx_min, y_idx_max, z_idx_min, z_idx_max;
115        double intersect_x, intersect_y, intersect_z;
116        double AOI = 1.0;
117        double cell_coverage = 1.0;
118        double sizex, sizey, sizez;     /* dimensions of histogram */
119        int unitsx, unitsy, unitsz; /* histogram grid size */
120        double value;
121        float overlapping_cells;
122        float avg_feat_cells;
123        double gain;
124        float8 selectivity;
125
126
127        /*
128         * Search box completely miss histogram extent
129         */
130        if ( box->xmax < geogstats->xmin || box->xmin > geogstats->xmax ||
131                box->ymax < geogstats->ymin || box->ymin > geogstats->ymax ||
132                box->zmax < geogstats->zmin || box->zmin > geogstats->zmax)
133        {
134                POSTGIS_DEBUG(3, " search_box does not overlaps histogram, returning 0");
135
136                return 0.0;
137        }
138
139        /*
140         * Search box completely contains histogram extent
141         */
142        if ( box->xmax >= geogstats->xmax && box->xmin <= geogstats->xmin &&
143                box->ymax >= geogstats->ymax && box->ymin <= geogstats->ymin &&
144                box->zmax >= geogstats->zmax && box->zmin <= geogstats->zmin)
145        {
146                POSTGIS_DEBUG(3, " search_box contains histogram, returning 1");
147
148                return 1.0;
149        }
150
151        sizex = geogstats->xmax - geogstats->xmin;
152        sizey = geogstats->ymax - geogstats->ymin;
153        sizez = geogstats->zmax - geogstats->zmin;
154
155        unitsx = geogstats->unitsx;
156        unitsy = geogstats->unitsy;
157        unitsz = geogstats->unitsz;
158
159        POSTGIS_DEBUGF(3, " histogram has %d unitsx, %d unitsy, %d unitsz", unitsx, unitsy, unitsz);
160        POSTGIS_DEBUGF(3, " histogram geosize is %fx%fx%f", sizex, sizey, sizez);
161
162        /* Work out the coverage depending upon the number of dimensions */
163        switch ((int)geogstats->dims)
164        {
165        case 0:
166        case 1:
167                /* Non-existent or multiple single points */
168                cell_coverage = 1;
169                break;
170
171        case 2:
172                /* Work in areas for 2 dimensions: we have to work slightly
173                   harder here to work out which 2 dimensions are valid */
174
175                if (sizez == 0)
176                        /* X and Y */
177                        cell_coverage = (sizex * sizey) / (unitsx * unitsy);
178                else if (sizey == 0)
179                        /* X and Z */
180                        cell_coverage = (sizex * sizez) / (unitsx * unitsz);
181                else if (sizex == 0)
182                        /* Y and Z */
183                        cell_coverage = (sizey * sizez) / (unitsy * unitsz);
184                break;
185
186        case 3:
187                /* Work in volumes for 3 dimensions */
188                cell_coverage = (sizex * sizey * sizey) / (unitsx * unitsy * unitsz);
189                break;
190        }
191
192        value = 0;
193
194        /* Find first overlapping x unit */
195        x_idx_min = (box->xmin-geogstats->xmin) / sizex * unitsx;
196        if (x_idx_min < 0)
197        {
198                POSTGIS_DEBUGF(3, " search_box overlaps %d x units outside (low) of histogram grid", -x_idx_min);
199
200                /* should increment the value somehow */
201                x_idx_min = 0;
202        }
203        if (x_idx_min >= unitsx)
204        {
205                POSTGIS_DEBUGF(3, " search_box overlaps %d x units outside (high) of histogram grid", x_idx_min-unitsx+1);
206
207                /* should increment the value somehow */
208                x_idx_min = unitsx-1;
209        }
210
211        /* Find first overlapping y unit */
212        y_idx_min = (box->ymin-geogstats->ymin) / sizey * unitsy;
213        if (y_idx_min <0)
214        {
215                POSTGIS_DEBUGF(3, " search_box overlaps %d y units outside (low) of histogram grid", -y_idx_min);
216
217                /* should increment the value somehow */
218                y_idx_min = 0;
219        }
220        if (y_idx_min >= unitsy)
221        {
222                POSTGIS_DEBUGF(3, " search_box overlaps %d y units outside (high) of histogram grid", y_idx_min-unitsy+1);
223
224                /* should increment the value somehow */
225                y_idx_min = unitsy-1;
226        }
227
228        /* Find first overlapping z unit */
229        z_idx_min = (box->zmin-geogstats->zmin) / sizez * unitsz;
230        if (z_idx_min <0)
231        {
232                POSTGIS_DEBUGF(3, " search_box overlaps %d z units outside (low) of histogram grid", -z_idx_min);
233
234                /* should increment the value somehow */
235                z_idx_min = 0;
236        }
237        if (z_idx_min >= unitsz)
238        {
239                POSTGIS_DEBUGF(3, " search_box overlaps %d z units outside (high) of histogram grid", z_idx_min-unitsz+1);
240
241                /* should increment the value somehow */
242                z_idx_min = unitsz-1;
243        }
244
245        /* Find last overlapping x unit */
246        x_idx_max = (box->xmax-geogstats->xmin) / sizex * unitsx;
247        if (x_idx_max <0)
248        {
249                /* should increment the value somehow */
250                x_idx_max = 0;
251        }
252        if (x_idx_max >= unitsx )
253        {
254                /* should increment the value somehow */
255                x_idx_max = unitsx-1;
256        }
257
258        /* Find last overlapping y unit */
259        y_idx_max = (box->ymax-geogstats->ymin) / sizey * unitsy;
260        if (y_idx_max <0)
261        {
262                /* should increment the value somehow */
263                y_idx_max = 0;
264        }
265        if (y_idx_max >= unitsy)
266        {
267                /* should increment the value somehow */
268                y_idx_max = unitsy-1;
269        }
270
271        /* Find last overlapping z unit */
272        z_idx_max = (box->zmax-geogstats->zmin) / sizez * unitsz;
273        if (z_idx_max <0)
274        {
275                /* should increment the value somehow */
276                z_idx_max = 0;
277        }
278        if (z_idx_max >= unitsz)
279        {
280                /* should increment the value somehow */
281                z_idx_max = unitsz-1;
282        }
283
284        /*
285         * the {x,y,z}_idx_{min,max}
286         * define the grid squares that the box intersects
287         */
288        for (z=z_idx_min; z<=z_idx_max; z++)
289        {
290                for (y=y_idx_min; y<=y_idx_max; y++)
291                {
292                        for (x=x_idx_min; x<=x_idx_max; x++)
293                        {
294                                double val;
295                                double gain;
296
297                                val = geogstats->value[x + y * unitsx + z * unitsx * unitsy];
298
299                                /*
300                                * Of the cell value we get
301                                * only the overlap fraction.
302                                */
303
304                                intersect_x = LW_MIN(box->xmax, geogstats->xmin + (x+1) * sizex / unitsx) - LW_MAX(box->xmin, geogstats->xmin + x * sizex / unitsx);
305                                intersect_y = LW_MIN(box->ymax, geogstats->ymin + (y+1) * sizey / unitsy) - LW_MAX(box->ymin, geogstats->ymin + y * sizey / unitsy);
306                                intersect_z = LW_MIN(box->zmax, geogstats->zmin + (z+1) * sizez / unitsz) - LW_MAX(box->zmin, geogstats->zmin + z * sizez / unitsz);
307
308
309                                switch ((int)geogstats->dims)
310                                {
311                                case 0:
312                                        AOI = 1;
313                                case 1:
314                                        /* Working in 1 dimension: work out the value we need
315                                           for AOI */
316                                        if (sizex == 0 && sizey == 0)
317                                                AOI = intersect_z;
318                                        else if (sizex == 0 && sizez == 0)
319                                                AOI = intersect_y;
320                                        else if (sizey == 0 && sizez == 0)
321                                                AOI = intersect_x;
322                                        break;
323
324                                case 2:
325                                        /* Working in 2 dimensions: work out which 2 values we
326                                           need for AOI */
327                                        if (sizex == 0)
328                                                AOI = intersect_y * intersect_z;
329                                        else if (sizey == 0)
330                                                AOI = intersect_x * intersect_z;
331                                        else if (sizez == 0)
332                                                AOI = intersect_x * intersect_y;
333                                        break;
334
335                                case 3:
336                                        /* Working in 3 dimensions: use all 3 values */
337                                        AOI = intersect_x * intersect_y * intersect_z;
338                                        break;
339                                }
340
341                                gain = AOI/cell_coverage;
342
343                                POSTGIS_DEBUGF(4, " [%d,%d,%d] cell val %.15f",
344                                               x, y, z, val);
345                                POSTGIS_DEBUGF(4, " [%d,%d,%d] AOI %.15f",
346                                               x, y, z, AOI);
347                                POSTGIS_DEBUGF(4, " [%d,%d,%d] gain %.15f",
348                                               x, y, z, gain);
349
350                                val *= gain;
351
352                                POSTGIS_DEBUGF(4, " [%d,%d,%d] adding %.15f to value",
353                                               x, y, z, val);
354
355                                value += val;
356                        }
357                }
358        }
359
360
361        /*
362         * If the search_box is a point, it will
363         * overlap a single cell and thus get
364         * it's value, which is the fraction of
365         * samples (we can presume of row set also)
366         * which bumped to that cell.
367         *
368         * If the table features are points, each
369         * of them will overlap a single histogram cell.
370         * Our search_box value would then be correctly
371         * computed as the sum of the bumped cells values.
372         *
373         * If both our search_box AND the sample features
374         * overlap more then a single histogram cell we
375         * need to consider the fact that our sum computation
376         * will have many duplicated included. E.g. each
377         * single sample feature would have contributed to
378         * raise the search_box value by as many times as
379         * many cells in the histogram are commonly overlapped
380         * by both searc_box and feature. We should then
381         * divide our value by the number of cells in the virtual
382         * 'intersection' between average feature cell occupation
383         * and occupation of the search_box. This is as
384         * fuzzy as you understand it :)
385         *
386         * Consistency check: whenever the number of cells is
387         * one of whichever part (search_box_occupation,
388         * avg_feature_occupation) the 'intersection' must be 1.
389         * If sounds that our 'intersaction' is actually the
390         * minimun number between search_box_occupation and
391         * avg_feat_occupation.
392         *
393         */
394        overlapping_cells = (x_idx_max-x_idx_min+1) * (y_idx_max-y_idx_min+1) * (z_idx_max-z_idx_min+1);
395        avg_feat_cells = geogstats->avgFeatureCells;
396
397        POSTGIS_DEBUGF(3, " search_box overlaps %f cells", overlapping_cells);
398        POSTGIS_DEBUGF(3, " avg feat overlaps %f cells", avg_feat_cells);
399
400        if ( ! overlapping_cells )
401        {
402                POSTGIS_DEBUG(3, " no overlapping cells, returning 0.0");
403
404                return 0.0;
405        }
406
407        gain = 1 / LW_MIN(overlapping_cells, avg_feat_cells);
408        selectivity = value * gain;
409
410        POSTGIS_DEBUGF(3, " SUM(ov_histo_cells)=%f", value);
411        POSTGIS_DEBUGF(3, " gain=%f", gain);
412        POSTGIS_DEBUGF(3, " selectivity=%f", selectivity);
413
414        /* prevent rounding overflows */
415        if (selectivity > 1.0) selectivity = 1.0;
416        else if (selectivity < 0) selectivity = 0.0;
417
418        return selectivity;
419}
420
421
422/**
423 * This function should return an estimation of the number of
424 * rows returned by a query involving an overlap check
425 * ( it's the restrict function for the && operator )
426 *
427 * It can make use (if available) of the statistics collected
428 * by the geometry analyzer function.
429 *
430 * Note that the good work is done by estimate_selectivity() above.
431 * This function just tries to find the search_box, loads the statistics
432 * and invoke the work-horse.
433 *
434 */
435PG_FUNCTION_INFO_V1(geography_gist_selectivity);
436Datum geography_gist_selectivity(PG_FUNCTION_ARGS)
437{
438        PlannerInfo *root = (PlannerInfo *) PG_GETARG_POINTER(0);
439
440        /* Oid operator = PG_GETARG_OID(1); */
441        List *args = (List *) PG_GETARG_POINTER(2);
442        /* int varRelid = PG_GETARG_INT32(3); */
443        Oid relid;
444        HeapTuple stats_tuple;
445        GEOG_STATS *geogstats;
446        /*
447         * This is to avoid casting the corresponding
448         * "type-punned" pointer, which would break
449         * "strict-aliasing rules".
450         */
451        GEOG_STATS **gsptr=&geogstats;
452        int geogstats_nvalues = 0;
453        Node *other;
454        Var *self;
455        GSERIALIZED *serialized;
456        LWGEOM *geometry;
457        GBOX search_box;
458        float8 selectivity = 0;
459
460        POSTGIS_DEBUG(2, "geography_gist_selectivity called");
461
462        /* Fail if not a binary opclause (probably shouldn't happen) */
463        if (list_length(args) != 2)
464        {
465                POSTGIS_DEBUG(3, "geography_gist_selectivity: not a binary opclause");
466
467                PG_RETURN_FLOAT8(DEFAULT_GEOGRAPHY_SEL);
468        }
469
470        /*
471         * This selectivity function is invoked by a clause of the form <arg> && <arg>
472         *
473         * In typical usage, one argument will be a column reference, while the other will
474         * be a geography constant; set self to point to the column argument and other
475         * to point to the constant argument.
476         */
477        other = (Node *) linitial(args);
478        if ( ! IsA(other, Const) )
479        {
480                self = (Var *)other;
481                other = (Node *) lsecond(args);
482        }
483        else
484        {
485                self = (Var *) lsecond(args);
486        }
487
488        if ( ! IsA(other, Const) )
489        {
490                POSTGIS_DEBUG(3, " no constant arguments - returning default selectivity");
491
492                PG_RETURN_FLOAT8(DEFAULT_GEOGRAPHY_SEL);
493        }
494
495        /*
496         * We are working on two constants..
497         * TODO: check if expression is true,
498         *       returned set would be either
499         *       the whole or none.
500         */
501        if ( ! IsA(self, Var) )
502        {
503                POSTGIS_DEBUG(3, " no variable argument ? - returning default selectivity");
504
505                PG_RETURN_FLOAT8(DEFAULT_GEOGRAPHY_SEL);
506        }
507
508        /*
509         * Convert the constant to a GBOX
510         */
511        serialized = (GSERIALIZED *)PG_DETOAST_DATUM( ((Const*)other)->constvalue );
512        geometry = lwgeom_from_gserialized(serialized);
513
514        /* Convert coordinates to 3D geodesic */
515        FLAGS_SET_GEODETIC(search_box.flags, 1);
516        if (!lwgeom_calculate_gbox_geodetic(geometry, &search_box))
517        {
518                POSTGIS_DEBUG(3, " search box cannot be calculated");
519
520                PG_RETURN_FLOAT8(0.0);
521        }
522
523        POSTGIS_DEBUGF(4, " requested search box is : %.15g %.15g %.15g, %.15g %.15g %.15g",
524                       search_box.xmin, search_box.ymin, search_box.zmin,
525                       search_box.xmax, search_box.ymax, search_box.zmax);
526
527        /*
528         * Get pg_statistic row
529         */
530        relid = getrelid(self->varno, root->parse->rtable);
531
532        stats_tuple = SearchSysCache(STATRELATT, ObjectIdGetDatum(relid), Int16GetDatum(self->varattno), 0, 0);
533        if ( ! stats_tuple )
534        {
535                POSTGIS_DEBUG(3, " No statistics, returning default estimate");
536
537                PG_RETURN_FLOAT8(DEFAULT_GEOGRAPHY_SEL);
538        }
539
540
541        if ( ! get_attstatsslot(stats_tuple, 0, 0, STATISTIC_KIND_GEOGRAPHY, InvalidOid, NULL, NULL,
542#if POSTGIS_PGSQL_VERSION >= 85
543                                NULL,
544#endif
545                                (float4 **)gsptr, &geogstats_nvalues) )
546        {
547                POSTGIS_DEBUG(3, " STATISTIC_KIND_GEOGRAPHY stats not found - returning default geography selectivity");
548
549                ReleaseSysCache(stats_tuple);
550                PG_RETURN_FLOAT8(DEFAULT_GEOGRAPHY_SEL);
551        }
552
553        POSTGIS_DEBUGF(4, " %d read from stats", geogstats_nvalues);
554
555        POSTGIS_DEBUGF(4, " histo: xmin,ymin,zmin: %f,%f,%f", geogstats->xmin, geogstats->ymin, geogstats->zmin);
556        POSTGIS_DEBUGF(4, " histo: xmax,ymax: %f,%f,%f", geogstats->xmax, geogstats->ymax, geogstats->zmax);
557        POSTGIS_DEBUGF(4, " histo: unitsx: %f", geogstats->unitsx);
558        POSTGIS_DEBUGF(4, " histo: unitsy: %f", geogstats->unitsy);
559        POSTGIS_DEBUGF(4, " histo: unitsz: %f", geogstats->unitsz);
560        POSTGIS_DEBUGF(4, " histo: avgFeatureCoverage: %f", geogstats->avgFeatureCoverage);
561        POSTGIS_DEBUGF(4, " histo: avgFeatureCells: %f", geogstats->avgFeatureCells);
562
563        /*
564         * Do the estimation
565         */
566        selectivity = estimate_selectivity(&search_box, geogstats);
567
568        POSTGIS_DEBUGF(3, " returning computed value: %f", selectivity);
569
570        free_attstatsslot(0, NULL, 0, (float *)geogstats, geogstats_nvalues);
571        ReleaseSysCache(stats_tuple);
572        PG_RETURN_FLOAT8(selectivity);
573}
574
575
576/**
577* JOIN selectivity in the GiST && operator
578* for all PG versions
579*/
580PG_FUNCTION_INFO_V1(geography_gist_join_selectivity);
581Datum geography_gist_join_selectivity(PG_FUNCTION_ARGS)
582{
583        PlannerInfo *root = (PlannerInfo *) PG_GETARG_POINTER(0);
584
585        /* Oid operator = PG_GETARG_OID(1); */
586        List *args = (List *) PG_GETARG_POINTER(2);
587        JoinType jointype = (JoinType) PG_GETARG_INT16(3);
588
589        Node *arg1, *arg2;
590        Var *var1, *var2;
591        Oid relid1, relid2;
592
593        HeapTuple stats1_tuple, stats2_tuple;
594        GEOG_STATS *geogstats1, *geogstats2;
595        /*
596        * These are to avoid casting the corresponding
597        * "type-punned" pointers, which would break
598        * "strict-aliasing rules".
599        */
600        GEOG_STATS **gs1ptr=&geogstats1, **gs2ptr=&geogstats2;
601        int geogstats1_nvalues = 0, geogstats2_nvalues = 0;
602        float8 selectivity1 = 0.0, selectivity2 = 0.0;
603        float4 num1_tuples = 0.0, num2_tuples = 0.0;
604        float4 total_tuples = 0.0, rows_returned = 0.0;
605        GBOX search_box;
606
607
608        /**
609        * Join selectivity algorithm. To calculation the selectivity we
610        * calculate the intersection of the two column sample extents,
611        * sum the results, and then multiply by two since for each
612        * geometry in col 1 that intersects a geometry in col 2, the same
613        * will also be true.
614        */
615
616        POSTGIS_DEBUGF(3, "geography_gist_join_selectivity called with jointype %d", jointype);
617
618        /*
619        * We'll only respond to an inner join/unknown context join
620        */
621        if (jointype != JOIN_INNER)
622        {
623                elog(NOTICE, "geography_gist_join_selectivity called with incorrect join type");
624                PG_RETURN_FLOAT8(DEFAULT_GEOGRAPHY_SEL);
625        }
626
627        /*
628        * Determine the oids of the geometry columns we are working with
629        */
630        arg1 = (Node *) linitial(args);
631        arg2 = (Node *) lsecond(args);
632
633        if (!IsA(arg1, Var) || !IsA(arg2, Var))
634        {
635                elog(DEBUG1, "geography_gist_join_selectivity called with arguments that are not column references");
636                PG_RETURN_FLOAT8(DEFAULT_GEOGRAPHY_SEL);
637        }
638
639        var1 = (Var *)arg1;
640        var2 = (Var *)arg2;
641
642        relid1 = getrelid(var1->varno, root->parse->rtable);
643        relid2 = getrelid(var2->varno, root->parse->rtable);
644
645        POSTGIS_DEBUGF(3, "Working with relations oids: %d %d", relid1, relid2);
646
647        /* Read the stats tuple from the first column */
648        stats1_tuple = SearchSysCache(STATRELATT, ObjectIdGetDatum(relid1), Int16GetDatum(var1->varattno), 0, 0);
649        if ( ! stats1_tuple )
650        {
651                POSTGIS_DEBUG(3, " No statistics, returning default geometry join selectivity");
652
653                PG_RETURN_FLOAT8(DEFAULT_GEOGRAPHY_SEL);
654        }
655
656        if ( ! get_attstatsslot(stats1_tuple, 0, 0, STATISTIC_KIND_GEOGRAPHY, InvalidOid, NULL, NULL,
657#if POSTGIS_PGSQL_VERSION >= 85
658                                NULL,
659#endif
660                                (float4 **)gs1ptr, &geogstats1_nvalues) )
661        {
662                POSTGIS_DEBUG(3, " STATISTIC_KIND_GEOGRAPHY stats not found - returning default geometry join selectivity");
663
664                ReleaseSysCache(stats1_tuple);
665                PG_RETURN_FLOAT8(DEFAULT_GEOGRAPHY_SEL);
666        }
667
668
669        /* Read the stats tuple from the second column */
670        stats2_tuple = SearchSysCache(STATRELATT, ObjectIdGetDatum(relid2), Int16GetDatum(var2->varattno), 0, 0);
671        if ( ! stats2_tuple )
672        {
673                POSTGIS_DEBUG(3, " No statistics, returning default geometry join selectivity");
674
675                free_attstatsslot(0, NULL, 0, (float *)geogstats1, geogstats1_nvalues);
676                ReleaseSysCache(stats1_tuple);
677                PG_RETURN_FLOAT8(DEFAULT_GEOGRAPHY_SEL);
678        }
679
680        if ( ! get_attstatsslot(stats2_tuple, 0, 0, STATISTIC_KIND_GEOGRAPHY, InvalidOid, NULL, NULL,
681#if POSTGIS_PGSQL_VERSION >= 85
682                                NULL,
683#endif
684                                (float4 **)gs2ptr, &geogstats2_nvalues) )
685        {
686                POSTGIS_DEBUG(3, " STATISTIC_KIND_GEOGRAPHY stats not found - returning default geometry join selectivity");
687
688                free_attstatsslot(0, NULL, 0, (float *)geogstats1, geogstats1_nvalues);
689                ReleaseSysCache(stats2_tuple);
690                ReleaseSysCache(stats1_tuple);
691                PG_RETURN_FLOAT8(DEFAULT_GEOGRAPHY_SEL);
692        }
693
694
695        /**
696        * Setup the search box - this is the intersection of the two column
697        * extents.
698        */
699        search_box.xmin = LW_MAX(geogstats1->xmin, geogstats2->xmin);
700        search_box.ymin = LW_MAX(geogstats1->ymin, geogstats2->ymin);
701        search_box.zmin = LW_MAX(geogstats1->zmin, geogstats2->zmin);
702        search_box.xmax = LW_MIN(geogstats1->xmax, geogstats2->xmax);
703        search_box.ymax = LW_MIN(geogstats1->ymax, geogstats2->ymax);
704        search_box.zmax = LW_MIN(geogstats1->zmax, geogstats2->zmax);
705
706        /* If the extents of the two columns don't intersect, return zero */
707        if (search_box.xmin > search_box.xmax || search_box.ymin > search_box.ymax ||
708                search_box.zmin > search_box.zmax)
709                PG_RETURN_FLOAT8(0.0);
710
711        POSTGIS_DEBUGF(3, " -- geomstats1 box: %.15g %.15g %.15g, %.15g %.15g %.15g", geogstats1->xmin, geogstats1->ymin, geogstats1->zmin, geogstats1->xmax, geogstats1->ymax, geogstats1->zmax);
712        POSTGIS_DEBUGF(3, " -- geomstats2 box: %.15g %.15g %.15g, %.15g %.15g %.15g", geogstats2->xmin, geogstats2->ymin, geogstats2->zmin, geogstats2->xmax, geogstats2->ymax, geogstats2->zmax);
713        POSTGIS_DEBUGF(3, " -- calculated intersection box is : %.15g %.15g %.15g, %.15g %.15g %.15g", search_box.xmin, search_box.ymin, search_box.zmin, search_box.xmax, search_box.ymax, search_box.zmax);
714
715
716        /* Do the selectivity */
717        selectivity1 = estimate_selectivity(&search_box, geogstats1);
718        selectivity2 = estimate_selectivity(&search_box, geogstats2);
719
720        POSTGIS_DEBUGF(3, "selectivity1: %.15g   selectivity2: %.15g", selectivity1, selectivity2);
721
722        /*
723        * OK, so before we calculate the join selectivity we also need to
724        * know the number of tuples in each of the columns since
725        * estimate_selectivity returns the number of estimated tuples
726        * divided by the total number of tuples.
727        */
728        num1_tuples = geogstats1->totalrows;
729        num2_tuples = geogstats2->totalrows;
730
731        /* Free the statistic tuples */
732        free_attstatsslot(0, NULL, 0, (float *)geogstats1, geogstats1_nvalues);
733        ReleaseSysCache(stats1_tuple);
734
735        free_attstatsslot(0, NULL, 0, (float *)geogstats2, geogstats2_nvalues);
736        ReleaseSysCache(stats2_tuple);
737
738        /*
739        * Finally calculate the estimate of the number of rows returned
740        *
741        *    = 2 * (nrows from col1 + nrows from col2) /
742        *       total nrows in col1 x total nrows in col2
743        *
744        * The factor of 2 accounts for the fact that for each tuple in
745        * col 1 matching col 2,
746        * there will be another match in col 2 matching col 1
747        */
748        total_tuples = num1_tuples * num2_tuples;
749        rows_returned = 2 * ((num1_tuples * selectivity1) + (num2_tuples * selectivity2));
750
751        POSTGIS_DEBUGF(3, "Rows from rel1: %f", num1_tuples * selectivity1);
752        POSTGIS_DEBUGF(3, "Rows from rel2: %f", num2_tuples * selectivity2);
753        POSTGIS_DEBUGF(3, "Estimated rows returned: %f", rows_returned);
754
755        /*
756        * One (or both) tuple count is zero...
757        * We return default selectivity estimate.
758        * We could probably attempt at an estimate
759        * w/out looking at tables tuple count, with
760        * a function of selectivity1, selectivity2.
761        */
762        if ( ! total_tuples )
763        {
764                POSTGIS_DEBUG(3, "Total tuples == 0, returning default join selectivity");
765
766                PG_RETURN_FLOAT8(DEFAULT_GEOGRAPHY_SEL);
767        }
768
769        if ( rows_returned > total_tuples )
770                PG_RETURN_FLOAT8(1.0);
771
772        PG_RETURN_FLOAT8(rows_returned / total_tuples);
773}
774
775
776/**
777 * This function is called by the analyze function iff
778 * the geography_analyze() function give it its pointer
779 * (this is always the case so far).
780 * The geography_analyze() function is also responsible
781 * of deciding the number of "sample" rows we will receive
782 * here. It is able to give use other 'custom' data, but we
783 * won't use them so far.
784 *
785 * Our job is to build some statistics on the sample data
786 * for use by operator estimators.
787 *
788 * Currently we only need statistics to estimate the number of rows
789 * overlapping a given extent (estimation function bound
790 * to the && operator).
791 *
792 */
793
794static void
795compute_geography_stats(VacAttrStats *stats, AnalyzeAttrFetchFunc fetchfunc,
796                        int samplerows, double totalrows)
797{
798        MemoryContext old_context;
799        GEOG_STATS *geogstats;
800
801        GBOX gbox;
802
803        GBOX *sample_extent = NULL;
804        GBOX **sampleboxes;
805        GBOX histobox;
806        int histocells;
807        double sizex, sizey, sizez, edgelength;
808        int unitsx = 0, unitsy = 0, unitsz = 0;
809        int geog_stats_size;
810        struct dimensions histodims[3];
811        int ndims;
812
813        double total_width = 0;
814        int notnull_cnt = 0, examinedsamples = 0, total_count_cells=0, total_cells_coverage = 0;
815
816#if USE_STANDARD_DEVIATION
817        /* for standard deviation */
818        double avgLOWx, avgLOWy, avgLOWz, avgHIGx, avgHIGy, avgHIGz;
819        double sumLOWx = 0, sumLOWy = 0, sumLOWz = 0, sumHIGx = 0, sumHIGy = 0, sumHIGz = 0;
820        double sdLOWx = 0, sdLOWy = 0, sdLOWz = 0, sdHIGx = 0, sdHIGy = 0, sdHIGz = 0;
821        GBOX *newhistobox = NULL;
822#endif
823
824        bool isnull;
825        int i;
826
827        POSTGIS_DEBUG(2, "compute_geography_stats called");
828
829        /*
830         * We'll build an histogram having from 40 to 400 boxesPerSide
831         * Total number of cells is determined by attribute stat
832         * target. It can go from  1600 to 160000 (stat target: 10,1000)
833         */
834        histocells = 160 * stats->attr->attstattarget;
835
836        /*
837         * Memory to store the bounding boxes from all of the sampled rows
838         */
839        sampleboxes = palloc(sizeof(GBOX *) * samplerows);
840
841        /* Mark the GBOX as being geodetic */
842        FLAGS_SET_GEODETIC(gbox.flags, 1);
843
844        /*
845         * First scan:
846         *  o find extent of the sample rows
847         *  o count null-infinite/not-null values
848         *  o compute total_width
849         *  o compute total features's box area (for avgFeatureArea)
850         *  o sum features box coordinates (for standard deviation)
851         */
852        for (i = 0; i < samplerows; i++)
853        {
854                Datum datum;
855                GSERIALIZED *serialized;
856                LWGEOM *geometry;
857
858                /* Fetch the datum and cast it into a geography */
859                datum = fetchfunc(stats, i, &isnull);
860
861                /* Skip nulls */
862                if (isnull)
863                        continue;
864
865                serialized = (GSERIALIZED *)PG_DETOAST_DATUM(datum);
866                geometry = lwgeom_from_gserialized(serialized);
867
868                /* Convert coordinates to 3D geodesic */
869                if (!lwgeom_calculate_gbox_geodetic(geometry, &gbox))
870                {
871                        /* Unable to obtain or calculate a bounding box */
872                        POSTGIS_DEBUGF(3, "skipping geometry at position %d", i);
873
874                        continue;
875                }
876
877                /*
878                 * Skip infinite geoms
879                 */
880                if ( ! finite(gbox.xmin) || ! finite(gbox.xmax) ||
881                        ! finite(gbox.ymin) || ! finite(gbox.ymax) ||
882                        ! finite(gbox.zmin) || ! finite(gbox.zmax) )
883                {
884                        POSTGIS_DEBUGF(3, " skipped infinite geometry at position %d", i);
885
886                        continue;
887                }
888
889                /*
890                 * Store bounding box in array
891                 */
892                sampleboxes[notnull_cnt] = palloc(sizeof(GBOX));
893                memcpy(sampleboxes[notnull_cnt], &gbox, sizeof(GBOX));
894
895                /*
896                 * Add to sample extent union
897                 */
898                if ( ! sample_extent )
899                {
900                        sample_extent = palloc(sizeof(GBOX));
901                        memcpy(sample_extent, &gbox, sizeof(GBOX));
902                }
903                else
904                {
905                        sample_extent->xmax = LWGEOM_Maxf(sample_extent->xmax, gbox.xmax);
906                        sample_extent->ymax = LWGEOM_Maxf(sample_extent->ymax, gbox.ymax);
907                        sample_extent->zmax = LWGEOM_Maxf(sample_extent->zmax, gbox.zmax);
908                        sample_extent->xmin = LWGEOM_Minf(sample_extent->xmin, gbox.xmin);
909                        sample_extent->ymin = LWGEOM_Minf(sample_extent->ymin, gbox.ymin);
910                        sample_extent->zmin = LWGEOM_Minf(sample_extent->zmin, gbox.zmin);
911                }
912
913                /** TODO: ask if we need geom or bvol size for stawidth */
914                total_width += serialized->size;
915
916#if USE_STANDARD_DEVIATION
917                /*
918                 * Add bvol coordinates to sum for standard deviation
919                 * computation.
920                 */
921                sumLOWx += gbox.xmin;
922                sumLOWy += gbox.ymin;
923                sumLOWz += gbox.zmin;
924                sumHIGx += gbox.xmax;
925                sumHIGy += gbox.ymax;
926                sumHIGz += gbox.zmax;
927#endif
928
929                notnull_cnt++;
930
931                /* give backend a chance of interrupting us */
932                vacuum_delay_point();
933        }
934
935        POSTGIS_DEBUG(3, "End of 1st scan:");
936        POSTGIS_DEBUGF(3, " Sample extent (min, max): (%g %g %g), (%g %g %g)", sample_extent->xmin, sample_extent->ymin,
937                       sample_extent->zmin, sample_extent->xmax, sample_extent->ymax, sample_extent->zmax);
938        POSTGIS_DEBUGF(3, " No. of geometries sampled: %d", samplerows);
939        POSTGIS_DEBUGF(3, " No. of non-null geometries sampled: %d", notnull_cnt);
940
941        if ( ! notnull_cnt )
942        {
943                elog(NOTICE, " no notnull values, invalid stats");
944                stats->stats_valid = false;
945                return;
946        }
947
948#if USE_STANDARD_DEVIATION
949
950        POSTGIS_DEBUG(3, "Standard deviation filter enabled");
951
952        /*
953         * Second scan:
954         *  o compute standard deviation
955         */
956        avgLOWx = sumLOWx / notnull_cnt;
957        avgLOWy = sumLOWy / notnull_cnt;
958        avgLOWz = sumLOWz / notnull_cnt;
959        avgHIGx = sumHIGx / notnull_cnt;
960        avgHIGy = sumHIGy / notnull_cnt;
961        avgHIGz = sumHIGz / notnull_cnt;
962
963        for (i = 0; i < notnull_cnt; i++)
964        {
965                GBOX *box;
966                box = (GBOX *)sampleboxes[i];
967
968                sdLOWx += (box->xmin - avgLOWx) * (box->xmin - avgLOWx);
969                sdLOWy += (box->ymin - avgLOWy) * (box->ymin - avgLOWy);
970                sdLOWz += (box->zmin - avgLOWz) * (box->zmin - avgLOWz);
971                sdHIGx += (box->xmax - avgHIGx) * (box->xmax - avgHIGx);
972                sdHIGy += (box->ymax - avgHIGy) * (box->ymax - avgHIGy);
973                sdHIGz += (box->zmax - avgHIGz) * (box->zmax - avgHIGz);
974        }
975        sdLOWx = sqrt(sdLOWx / notnull_cnt);
976        sdLOWy = sqrt(sdLOWy / notnull_cnt);
977        sdLOWz = sqrt(sdLOWz / notnull_cnt);
978        sdHIGx = sqrt(sdHIGx / notnull_cnt);
979        sdHIGy = sqrt(sdHIGy / notnull_cnt);
980        sdHIGz = sqrt(sdHIGz / notnull_cnt);
981
982        POSTGIS_DEBUG(3, " standard deviations:");
983        POSTGIS_DEBUGF(3, "  LOWx - avg:%f sd:%f", avgLOWx, sdLOWx);
984        POSTGIS_DEBUGF(3, "  LOWy - avg:%f sd:%f", avgLOWy, sdLOWy);
985        POSTGIS_DEBUGF(3, "  LOWz - avg:%f sd:%f", avgLOWz, sdLOWz);
986        POSTGIS_DEBUGF(3, "  HIGx - avg:%f sd:%f", avgHIGx, sdHIGx);
987        POSTGIS_DEBUGF(3, "  HIGy - avg:%f sd:%f", avgHIGy, sdHIGy);
988        POSTGIS_DEBUGF(3, "  HIGz - avg:%f sd:%f", avgHIGz, sdHIGz);
989
990        histobox.xmin = LW_MAX((avgLOWx - SDFACTOR * sdLOWx), sample_extent->xmin);
991        histobox.ymin = LW_MAX((avgLOWy - SDFACTOR * sdLOWy), sample_extent->ymin);
992        histobox.zmin = LW_MAX((avgLOWz - SDFACTOR * sdLOWz), sample_extent->zmin);
993        histobox.xmax = LW_MIN((avgHIGx + SDFACTOR * sdHIGx), sample_extent->xmax);
994        histobox.ymax = LW_MIN((avgHIGy + SDFACTOR * sdHIGy), sample_extent->ymax);
995        histobox.zmax = LW_MIN((avgHIGz + SDFACTOR * sdHIGz), sample_extent->zmax);
996
997        POSTGIS_DEBUGF(3, " sd_extent: xmin, ymin, zmin: %f, %f, %f",
998                       histobox.xmin, histobox.ymin, histobox.zmin);
999        POSTGIS_DEBUGF(3, " sd_extent: xmax, ymax, zmax: %f, %f, %f",
1000                       histobox.xmax, histobox.ymax, histobox.zmax);
1001
1002        /*
1003         * Third scan:
1004         *   o skip hard deviants
1005         *   o compute new histogram box
1006         */
1007        for (i = 0; i < notnull_cnt; i++)
1008        {
1009                GBOX *box;
1010                box = (GBOX *)sampleboxes[i];
1011
1012                if ( box->xmin > histobox.xmax || box->xmax < histobox.xmin ||
1013                        box->ymin > histobox.ymax || box->ymax < histobox.ymin ||
1014                        box->zmin > histobox.zmax || box->zmax < histobox.zmin)
1015                {
1016                        POSTGIS_DEBUGF(4, " feat %d is an hard deviant, skipped", i);
1017
1018                        sampleboxes[i] = NULL;
1019                        continue;
1020                }
1021
1022                if ( ! newhistobox )
1023                {
1024                        newhistobox = palloc(sizeof(GBOX));
1025                        memcpy(newhistobox, box, sizeof(GBOX));
1026                }
1027                else
1028                {
1029                        if ( box->xmin < newhistobox->xmin )
1030                                newhistobox->xmin = box->xmin;
1031                        if ( box->ymin < newhistobox->ymin )
1032                                newhistobox->ymin = box->ymin;
1033                        if ( box->zmin < newhistobox->zmin )
1034                                newhistobox->zmin = box->zmin;
1035                        if ( box->xmax > newhistobox->xmax )
1036                                newhistobox->xmax = box->xmax;
1037                        if ( box->ymax > newhistobox->ymax )
1038                                newhistobox->ymax = box->ymax;
1039                        if ( box->zmax > newhistobox->zmax )
1040                                newhistobox->zmax = box->zmax;
1041                }
1042        }
1043
1044        /* If everything was a deviant, the new histobox is the same as the old histobox */
1045        if ( ! newhistobox )
1046        {
1047                newhistobox = palloc(sizeof(GBOX));
1048                memcpy(newhistobox, &histobox, sizeof(GBOX));
1049        }
1050
1051        /*
1052         * Set histogram extent as the intersection between
1053         * standard deviation based histogram extent
1054         * and computed sample extent after removal of
1055         * hard deviants (there might be no hard deviants).
1056         */
1057        if ( histobox.xmin < newhistobox->xmin )
1058                histobox.xmin = newhistobox->xmin;
1059        if ( histobox.ymin < newhistobox->ymin )
1060                histobox.ymin = newhistobox->ymin;
1061        if ( histobox.zmin < newhistobox->zmin )
1062                histobox.zmin = newhistobox->zmin;
1063        if ( histobox.xmax > newhistobox->xmax )
1064                histobox.xmax = newhistobox->xmax;
1065        if ( histobox.ymax > newhistobox->ymax )
1066                histobox.ymax = newhistobox->ymax;
1067        if ( histobox.zmax > newhistobox->zmax )
1068                histobox.zmax = newhistobox->zmax;
1069
1070#else /* ! USE_STANDARD_DEVIATION */
1071
1072        /*
1073        * Set histogram extent box
1074        */
1075        histobox.xmin = sample_extent->xmin;
1076        histobox.ymin = sample_extent->ymin;
1077        histobox.zmin = sample_extent->zmin;
1078        histobox.xmax = sample_extent->xmax;
1079        histobox.ymax = sample_extent->ymax;
1080        histobox.zmax = sample_extent->zmax;
1081
1082#endif /* USE_STANDARD_DEVIATION */
1083
1084
1085        POSTGIS_DEBUGF(3, " histogram_extent: xmin, ymin, zmin: %f, %f, %f",
1086                       histobox.xmin, histobox.ymin, histobox.zmin);
1087        POSTGIS_DEBUGF(3, " histogram_extent: xmax, ymax, zmax: %f, %f, %f",
1088                       histobox.xmax, histobox.ymax, histobox.zmax);
1089
1090        /* Calculate the size of each dimension */
1091        sizex = histobox.xmax - histobox.xmin;
1092        sizey = histobox.ymax - histobox.ymin;
1093        sizez = histobox.zmax - histobox.zmin;
1094
1095        /* In order to calculate a suitable aspect ratio for the histogram, we need
1096           to work out how many dimensions exist within our sample data (which we
1097           assume is representative of the whole data) */
1098        ndims = 0;
1099        if (sizex != 0)
1100        {
1101                histodims[ndims].axis = 'X';
1102                histodims[ndims].min = histobox.xmin;
1103                histodims[ndims].max = histobox.xmax;
1104                ndims++;
1105        }
1106
1107        if (sizey != 0)
1108        {
1109                histodims[ndims].axis = 'Y';
1110                histodims[ndims].min = histobox.ymin;
1111                histodims[ndims].max = histobox.ymax;
1112
1113                ndims++;
1114        }
1115
1116        if (sizez != 0)
1117        {
1118                histodims[ndims].axis = 'Z';
1119                histodims[ndims].min = histobox.zmin;
1120                histodims[ndims].max = histobox.zmax;
1121
1122                ndims++;
1123        }
1124
1125        /* Based upon the number of dimensions, we now work out the number of units in each dimension.
1126           The number of units is defined as the number of cell blocks in each dimension which make
1127           up the total number of histocells; i.e. unitsx * unitsy * unitsz = histocells */
1128
1129        /* Note: geodetic data is currently indexed in 3 dimensions; however code for remaining dimensions
1130           is also included to allow for indexing 3D cartesian data at a later date */
1131
1132        POSTGIS_DEBUGF(3, "Number of dimensions in sample set: %d", ndims);
1133
1134        switch (ndims)
1135        {
1136        case 0:
1137                /* An empty column, or multiple points in exactly the same
1138                   position in space */
1139                unitsx = 1;
1140                unitsy = 1;
1141                unitsz = 1;
1142                histocells = 1;
1143                break;
1144
1145        case 1:
1146                /* Sample data all lies on a single line, so set the correct
1147                   units variables depending upon which axis is in use */
1148                for (i = 0; i < ndims; i++)
1149                {
1150                        if ( (histodims[i].max - histodims[i].min) != 0)
1151                        {
1152                                /* We've found the non-zero dimension, so set the
1153                                   units variables accordingly */
1154                                switch (histodims[i].axis)
1155                                {
1156                                case 'X':
1157                                        unitsx = histocells;
1158                                        unitsy = 1;
1159                                        unitsz = 1;
1160                                        break;
1161
1162                                case 'Y':
1163                                        unitsx = 1;
1164                                        unitsy = histocells;
1165                                        unitsz = 1;
1166                                        break;
1167
1168                                case 'Z':
1169                                        unitsx = 1;
1170                                        unitsy = 1;
1171                                        unitsz = histocells;
1172                                        break;
1173                                }
1174                        }
1175                }
1176                break;
1177
1178        case 2:
1179                /* Sample data lies within 2D space: divide the total area by the total
1180                   number of cells, and thus work out the edge size of the unit block */
1181                edgelength = sqrt(
1182                                 LW_ABS(histodims[0].max - histodims[0].min) *
1183                                 LW_ABS(histodims[1].max - histodims[1].min) / (double)histocells
1184                             );
1185
1186                /* The calculation is easy; the harder part is to work out which dimensions
1187                   we actually have to set the units variables appropriately */
1188                if (histodims[0].axis == 'X' && histodims[1].axis == 'Y')
1189                {
1190                        /* X and Y */
1191                        unitsx = LW_ABS(histodims[0].max - histodims[0].min) / edgelength;
1192                        unitsy = LW_ABS(histodims[1].max - histodims[1].min) / edgelength;
1193                        unitsz = 1;
1194                }
1195                else if (histodims[0].axis == 'Y' && histodims[1].axis == 'X')
1196                {
1197                        /* Y and X */
1198                        unitsx = LW_ABS(histodims[1].max - histodims[1].min) / edgelength;
1199                        unitsy = LW_ABS(histodims[0].max - histodims[0].min) / edgelength;
1200                        unitsz = 1;
1201                }
1202                else if (histodims[0].axis == 'X' && histodims[1].axis == 'Z')
1203                {
1204                        /* X and Z */
1205                        unitsx = LW_ABS(histodims[0].max - histodims[0].min) / edgelength;
1206                        unitsy = 1;
1207                        unitsz = LW_ABS(histodims[1].max - histodims[1].min) / edgelength;
1208                }
1209                else if (histodims[0].axis == 'Z' && histodims[1].axis == 'X')
1210                {
1211                        /* Z and X */
1212                        unitsx = LW_ABS(histodims[0].max - histodims[0].min) / edgelength;
1213                        unitsy = 1;
1214                        unitsz = LW_ABS(histodims[1].max - histodims[1].min) / edgelength;
1215                }
1216                else if (histodims[0].axis == 'Y' && histodims[1].axis == 'Z')
1217                {
1218                        /* Y and Z */
1219                        unitsx = 1;
1220                        unitsy = LW_ABS(histodims[0].max - histodims[0].min) / edgelength;
1221                        unitsz = LW_ABS(histodims[1].max - histodims[1].min) / edgelength;
1222                }
1223                else if (histodims[0].axis == 'Z' && histodims[1].axis == 'Y')
1224                {
1225                        /* Z and X */
1226                        unitsx = 1;
1227                        unitsy = LW_ABS(histodims[1].max - histodims[1].min) / edgelength;
1228                        unitsz = LW_ABS(histodims[0].max - histodims[0].min) / edgelength;
1229                }
1230
1231                break;
1232
1233        case 3:
1234                /* Sample data lies within 3D space: divide the total volume by the total
1235                   number of cells, and thus work out the edge size of the unit block */
1236                edgelength = pow(
1237                                 LW_ABS(histodims[0].max - histodims[0].min) *
1238                                 LW_ABS(histodims[1].max - histodims[1].min) *
1239                                 LW_ABS(histodims[2].max - histodims[2].min) / (double)histocells,
1240                                 (double)1/3);
1241
1242                /* Units are simple in 3 dimensions */
1243                unitsx = LW_ABS(histodims[0].max - histodims[0].min) / edgelength;
1244                unitsy = LW_ABS(histodims[1].max - histodims[1].min) / edgelength;
1245                unitsz = LW_ABS(histodims[2].max - histodims[2].min) / edgelength;
1246
1247                break;
1248        }
1249
1250        POSTGIS_DEBUGF(3, " computed histogram grid size (X,Y,Z): %d x %d x %d (%d out of %d cells)", unitsx, unitsy, unitsz, unitsx * unitsy * unitsz, histocells);
1251
1252        /*
1253         * Create the histogram (GEOG_STATS)
1254         */
1255        old_context = MemoryContextSwitchTo(stats->anl_context);
1256        geog_stats_size = sizeof(GEOG_STATS) + (histocells - 1) * sizeof(float4);
1257        geogstats = palloc(geog_stats_size);
1258        MemoryContextSwitchTo(old_context);
1259
1260        geogstats->dims = ndims;
1261        geogstats->xmin = histobox.xmin;
1262        geogstats->ymin = histobox.ymin;
1263        geogstats->zmin = histobox.zmin;
1264        geogstats->xmax = histobox.xmax;
1265        geogstats->ymax = histobox.ymax;
1266        geogstats->zmax = histobox.zmax;
1267        geogstats->unitsx = unitsx;
1268        geogstats->unitsy = unitsy;
1269        geogstats->unitsz = unitsz;
1270        geogstats->totalrows = totalrows;
1271
1272        /* Initialize all values to 0 */
1273        for (i = 0; i < histocells; i++)
1274                geogstats->value[i] = 0;
1275
1276
1277        /*
1278         * Fourth scan:
1279         *  o fill histogram values with the number of
1280         *    features' bbox overlaps: a feature's bvol
1281         *    can fully overlap (1) or partially overlap
1282         *    (fraction of 1) an histogram cell.
1283         *
1284         *  o compute total cells occupation
1285         *
1286         */
1287
1288        POSTGIS_DEBUG(3, "Beginning histogram intersection calculations");
1289
1290        for (i = 0; i < notnull_cnt; i++)
1291        {
1292                GBOX *box;
1293
1294                /* Note these array index variables are zero-based */
1295                int x_idx_min, x_idx_max, x;
1296                int y_idx_min, y_idx_max, y;
1297                int z_idx_min, z_idx_max, z;
1298                int numcells = 0;
1299
1300                box = (GBOX *)sampleboxes[i];
1301                if ( ! box ) continue; /* hard deviant.. */
1302
1303                /* give backend a chance of interrupting us */
1304                vacuum_delay_point();
1305
1306                POSTGIS_DEBUGF(4, " feat %d box is %f %f %f, %f %f %f",
1307                               i, box->xmax, box->ymax, box->zmax, box->xmin, box->ymin, box->zmin);
1308
1309                /* Find first overlapping unitsx cell */
1310                x_idx_min = (box->xmin - geogstats->xmin) / sizex * unitsx;
1311                if (x_idx_min <0) x_idx_min = 0;
1312                if (x_idx_min >= unitsx) x_idx_min = unitsx - 1;
1313
1314                /* Find first overlapping unitsy cell */
1315                y_idx_min = (box->ymin - geogstats->ymin) / sizey * unitsy;
1316                if (y_idx_min <0) y_idx_min = 0;
1317                if (y_idx_min >= unitsy) y_idx_min = unitsy - 1;
1318
1319                /* Find first overlapping unitsz cell */
1320                z_idx_min = (box->zmin - geogstats->zmin) / sizez * unitsz;
1321                if (z_idx_min <0) z_idx_min = 0;
1322                if (z_idx_min >= unitsz) z_idx_min = unitsz - 1;
1323
1324                /* Find last overlapping unitsx cell */
1325                x_idx_max = (box->xmax - geogstats->xmin) / sizex * unitsx;
1326                if (x_idx_max <0) x_idx_max = 0;
1327                if (x_idx_max >= unitsx ) x_idx_max = unitsx - 1;
1328
1329                /* Find last overlapping unitsy cell */
1330                y_idx_max = (box->ymax - geogstats->ymin) / sizey * unitsy;
1331                if (y_idx_max <0) y_idx_max = 0;
1332                if (y_idx_max >= unitsy) y_idx_max = unitsy - 1;
1333
1334                /* Find last overlapping unitsz cell */
1335                z_idx_max = (box->zmax - geogstats->zmin) / sizez * unitsz;
1336                if (z_idx_max <0) z_idx_max = 0;
1337                if (z_idx_max >= unitsz) z_idx_max = unitsz - 1;
1338
1339                POSTGIS_DEBUGF(4, " feat %d overlaps unitsx %d-%d, unitsy %d-%d, unitsz %d-%d",
1340                               i, x_idx_min, x_idx_max, y_idx_min, y_idx_max, z_idx_min, z_idx_max);
1341
1342                /* Calculate the feature coverage - this of course depends upon the number of dims */
1343                switch (ndims)
1344                {
1345                case 1:
1346                        total_cells_coverage++;
1347                        break;
1348
1349                case 2:
1350                        total_cells_coverage += (box->xmax - box->xmin) * (box->ymax - box->ymin);
1351                        break;
1352
1353                case 3:
1354                        total_cells_coverage += (box->xmax - box->xmin) * (box->ymax - box->ymin) *
1355                                                (box->zmax - box->zmin);
1356                        break;
1357                }
1358
1359                /*
1360                 * the {x,y,z}_idx_{min,max}
1361                 * define the grid squares that the box intersects
1362                 */
1363
1364                for (z = z_idx_min; z <= z_idx_max; z++)
1365                {
1366                        for (y = y_idx_min; y <= y_idx_max; y++)
1367                        {
1368                                for (x = x_idx_min; x <= x_idx_max; x++)
1369                                {
1370                                        geogstats->value[x + y * unitsx + z * unitsx * unitsy] += 1;
1371                                        numcells++;
1372                                }
1373                        }
1374                }
1375
1376                /*
1377                 * before adding to the total cells
1378                 * we could decide if we really
1379                 * want this feature to count
1380                 */
1381                total_count_cells += numcells;
1382
1383                examinedsamples++;
1384        }
1385
1386        POSTGIS_DEBUGF(3, " examined_samples: %d/%d", examinedsamples, samplerows);
1387
1388        if ( ! examinedsamples )
1389        {
1390                elog(NOTICE, " no examined values, invalid stats");
1391                stats->stats_valid = false;
1392
1393                POSTGIS_DEBUG(3, " no stats have been gathered");
1394
1395                return;
1396        }
1397
1398        /** TODO: what about null features (TODO) */
1399        geogstats->avgFeatureCells = (float4)total_count_cells / examinedsamples;
1400        geogstats->avgFeatureCoverage = total_cells_coverage / examinedsamples;
1401
1402        POSTGIS_DEBUGF(3, " histo: total_boxes_cells: %d", total_count_cells);
1403        POSTGIS_DEBUGF(3, " histo: avgFeatureCells: %f", geogstats->avgFeatureCells);
1404        POSTGIS_DEBUGF(3, " histo: avgFeatureCoverage: %f", geogstats->avgFeatureCoverage);
1405
1406        /*
1407         * Normalize histogram
1408         *
1409         * We divide each histogram cell value
1410         * by the number of samples examined.
1411         *
1412         */
1413        for (i = 0; i < histocells; i++)
1414                geogstats->value[i] /= examinedsamples;
1415
1416#if POSTGIS_DEBUG_LEVEL >= 4
1417        /* Dump the resulting histogram for analysis */
1418        {
1419                int x, y, z;
1420                for (x = 0; x < unitsx; x++)
1421                {
1422                        for (y = 0; y < unitsy; y++)
1423                        {
1424                                for (z = 0; z < unitsz; z++)
1425                                {
1426                                        POSTGIS_DEBUGF(4, " histo[%d,%d,%d] = %.15f", x, y, z,
1427                                                       geogstats->value[x + y * unitsx + z * unitsx * unitsy]);
1428                                }
1429                        }
1430                }
1431        }
1432#endif
1433
1434        /*
1435         * Write the statistics data
1436         */
1437        stats->stakind[0] = STATISTIC_KIND_GEOGRAPHY;
1438        stats->staop[0] = InvalidOid;
1439        stats->stanumbers[0] = (float4 *)geogstats;
1440        stats->numnumbers[0] = geog_stats_size/sizeof(float4);
1441
1442        stats->stanullfrac = (float4)(samplerows - notnull_cnt)/samplerows;
1443        stats->stawidth = total_width/notnull_cnt;
1444        stats->stadistinct = -1.0;
1445
1446        POSTGIS_DEBUGF(3, " out: slot 0: kind %d (STATISTIC_KIND_GEOGRAPHY)",
1447                       stats->stakind[0]);
1448        POSTGIS_DEBUGF(3, " out: slot 0: op %d (InvalidOid)", stats->staop[0]);
1449        POSTGIS_DEBUGF(3, " out: slot 0: numnumbers %d", stats->numnumbers[0]);
1450        POSTGIS_DEBUGF(3, " out: null fraction: %d/%d=%g", (samplerows - notnull_cnt), samplerows, stats->stanullfrac);
1451        POSTGIS_DEBUGF(3, " out: average width: %d bytes", stats->stawidth);
1452        POSTGIS_DEBUG(3, " out: distinct values: all (no check done)");
1453
1454        stats->stats_valid = true;
1455
1456}
1457
1458
1459/**
1460 * This function will be called when the ANALYZE command is run
1461 * on a column of the "geography" type.
1462 *
1463 * It will need to return a stats builder function reference
1464 * and a "minimum" sample rows to feed it.
1465 * If we want analysis to be completely skipped we can return
1466 * FALSE and leave output vals untouched.
1467 *
1468 * What we know from this call is:
1469 *
1470 *      o The pg_attribute row referring to the specific column.
1471 *        Could be used to get reltuples from pg_class (which
1472 *        might quite inexact though...) and use them to set an
1473 *        appropriate minimum number of sample rows to feed to
1474 *        the stats builder. The stats builder will also receive
1475 *        a more accurate "estimation" of the number or rows.
1476 *
1477 *      o The pg_type row for the specific column.
1478 *        Could be used to set stat builder / sample rows
1479 *        based on domain type (when postgis will be implemented
1480 *        that way).
1481 *
1482 * Being this experimental we'll stick to a static stat_builder/sample_rows
1483 * value for now.
1484 *
1485 */
1486
1487PG_FUNCTION_INFO_V1(geography_analyze);
1488Datum geography_analyze(PG_FUNCTION_ARGS)
1489{
1490        VacAttrStats *stats = (VacAttrStats *)PG_GETARG_POINTER(0);
1491        Form_pg_attribute attr = stats->attr;
1492
1493        POSTGIS_DEBUG(2, "geography_analyze called");
1494
1495        /* If the attstattarget column is negative, use the default value */
1496        /* NB: it is okay to scribble on stats->attr since it's a copy */
1497        if (attr->attstattarget < 0)
1498                attr->attstattarget = default_statistics_target;
1499
1500        POSTGIS_DEBUGF(3, " attribute stat target: %d", attr->attstattarget);
1501
1502        /* Setup the minimum rows and the algorithm function */
1503        stats->minrows = 300 * stats->attr->attstattarget;
1504        stats->compute_stats = compute_geography_stats;
1505
1506        POSTGIS_DEBUGF(3, " minrows: %d", stats->minrows);
1507
1508        /* Indicate we are done successfully */
1509        PG_RETURN_BOOL(true);
1510}
Note: See TracBrowser for help on using the repository browser.