pintos-os.org Git - pspp/blob - src/math/chart-geometry.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2004, 2015 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #include <config.h>
  18 #include <math.h>
  19 #include <float.h>
  20 #include <assert.h>
  21
  22 #include "chart-geometry.h"
  23 #include "decimal.h"
  24 #include <stdlib.h>
  25
  26 #include "gl/xalloc.h"
  27 #include "gl/minmax.h"
  28 #include "gl/xvasprintf.h"
  29
  30 static const double standard_tick[] = {1, 2, 5, 10};
  31
  32 /* Adjust tick to be a sensible value
  33    ie:  ... 0.1,0.2,0.5,   1,2,5,  10,20,50 ... */
  34 void
  35 chart_rounded_tick (double tick, struct decimal *result)
  36 {
  37   int i;
  38
  39   struct decimal ddif = {1, 1000};
  40
  41   /* Avoid arithmetic problems with very small values */
  42   if (fabs (tick) < DBL_EPSILON)
  43     {
  44       result->ordinate = 0;
  45       result->mantissa = 0;
  46       return;
  47     }
  48
  49   struct decimal dt;
  50   decimal_from_double (&dt, tick);
  51
  52   double expd = dec_log10 (&dt) - 1;
  53
  54   for (i = 0  ; i < 4 ; ++i)
  55     {
  56       struct decimal candidate;
  57       struct decimal delta;
  58
  59       decimal_init (&candidate, standard_tick[i], expd);
  60
  61       delta = dt;
  62       decimal_subtract (&delta, &candidate);
  63       delta.ordinate = llabs (delta.ordinate);
  64
  65       if (decimal_cmp (&delta, &ddif) < 0)
  66         {
  67           ddif = delta;
  68           *result = candidate;
  69         }
  70     }
  71 }
  72
  73 /*
  74    Find a set {LOWER, INTERVAL, N_TICKS} such that:
  75
  76    LOWER <= LOWDBL,
  77    LOWER + INTERVAL > LOWDBL,
  78
  79    LOWER + N_TICKS * INTERVAL < HIGHDBL
  80    LOWER + (N_TICKS + 1) * INTERVAL >= HIGHDBL
  81
  82    INTERVAL = X * 10^N
  83     where: N is integer
  84     and    X is an element of {1, 2, 5}
  85
  86    In other words:
  87
  88          INTERVAL
  89          >      <
  90      |....+....+....+.      .+....|
  91    LOWER  1    2    3     N_TICKS
  92         ^LOWDBL                 ^HIGHDBL
  93 */
  94 void
  95 chart_get_scale (double highdbl, double lowdbl,
  96                  struct decimal *lower,
  97                  struct decimal *interval,
  98                  int *n_ticks)
  99 {
 100   int i;
 101   double fitness = DBL_MAX;
 102   *n_ticks = 0;
 103   struct decimal high;
 104   struct decimal low;
 105
 106   assert (highdbl >= lowdbl);
 107
 108   decimal_from_double (&high, highdbl);
 109   decimal_from_double (&low, lowdbl);
 110
 111   struct decimal diff = high;
 112   decimal_subtract (&diff, &low);
 113
 114   double expd = dec_log10 (&diff) - 2;
 115
 116   /* Find the most pleasing interval */
 117   for (i = 1; i < 4; ++i)
 118     {
 119       struct decimal clbound = low;
 120       struct decimal cubound = high;
 121       struct decimal candidate;
 122       decimal_init (&candidate, standard_tick[i], expd);
 123
 124       decimal_divide (&clbound, &candidate);
 125       int fl = decimal_floor (&clbound);
 126       decimal_int_multiply (&candidate, fl);
 127       clbound = candidate;
 128
 129
 130       decimal_init (&candidate, standard_tick[i], expd);
 131       decimal_divide (&cubound, &candidate);
 132       int fu = decimal_ceil (&cubound);
 133       decimal_int_multiply (&candidate, fu);
 134
 135       cubound = candidate;
 136
 137       decimal_init (&candidate, standard_tick[i], expd);
 138       decimal_subtract (&cubound, &clbound);
 139       decimal_divide (&cubound, &candidate);
 140
 141
 142       ord_t n_int = decimal_floor (&cubound);
 143
 144       /* We prefer to have between 5 and 10 tick marks on a scale */
 145       double f = 1 - exp (-0.2 *  fabs (n_int - 7.5) / 7.5);
 146
 147       if (f < fitness)
 148         {
 149           fitness = f;
 150           *lower = clbound;
 151           *interval = candidate;
 152           *n_ticks = n_int;
 153         }
 154     }
 155 }
 156
 157 /*
 158  * Compute the optimum format string and the scaling
 159  * for the tick drawing on a chart axis
 160  * Input:  max:     the maximum value of the range
 161  *         min:     the minimum value of the range
 162  *         nticks:  the number of tick intervals (bins) on the axis
 163  * Return: fs:      format string for printf to print the tick value
 164  *         scale:   scaling factor for the tick value
 165  * The format string has to be freed after usage.
 166  * An example format string and scalefactor:
 167  * Non Scientific: "%.3lf", scale=1.00
 168  * Scientific:     "%.2lfe3", scale = 0.001
 169  * Usage example:
 170  *   fs = chart_get_ticks_format(95359943.3,34434.9,8,&scale);
 171  *   printf(fs,value*scale);
 172  *   free(fs);
 173  */
 174 char *
 175 chart_get_ticks_format (const double max, const double min,
 176                         const unsigned int nticks, double *scale)
 177 {
 178   assert(max > min);
 179   double interval = (max - min)/nticks;
 180   double logmax = log10(fmax(fabs(max),fabs(min)));
 181   double logintv = log10(interval);
 182   int logshift = 0;
 183   char *format_string = NULL;
 184   int nrdecs = 0;
 185
 186   if (logmax > 0.0 && logintv < 0.0)
 187     {
 188       nrdecs = MIN(6,(int)(fabs(logintv))+1);
 189       logshift = 0;
 190       format_string = xasprintf("%%.%dlf",nrdecs);
 191     }
 192   else if (logmax > 0.0) /*logintv is > 0*/
 193     {
 194       if (logintv < 5.0)
 195         {
 196           logshift = 0; /* No scientific format */
 197           nrdecs = 0;
 198           format_string = xstrdup("%.0lf");
 199         }
 200       else
 201         {
 202           logshift = (int)logmax;
 203           nrdecs = MIN(6,(int)(logmax-logintv)+1);
 204           format_string = xasprintf("%%.%dlf&#8729;10<sup>%d</sup>",nrdecs,logshift);
 205         }
 206     }
 207   else /* logmax and logintv are < 0 */
 208     {
 209       if (logmax > -4.0)
 210         {
 211           logshift = 0; /* No scientific format */
 212           nrdecs = (int)(-logintv) + 1;
 213           format_string = xasprintf("%%.%dlf",nrdecs);
 214         }
 215       else
 216         {
 217           logshift = (int)logmax-1;
 218           nrdecs = MIN(6,(int)(logmax-logintv)+1);
 219           format_string = xasprintf("%%.%dlf&#8729;10<sup>%d</sup>",nrdecs,logshift);
 220         }
 221     }
 222   *scale = pow(10.0,-(double)logshift);
 223   return format_string;
 224 }