The R Project SVN R

Rev

Rev 80633 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 r 1 
/*
564 maechler 2 
 *  R : A Computer Language for Statistical Data Analysis
80633 maechler 3 
 *  Copyright (C) 1997--2021  The R Core Team
75182 maechler 4 
 *  Copyright (C) 2002--2009  The R Foundation
2 r 5 
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
6 
 *
7 
 *  This program is free software; you can redistribute it and/or modify
8 
 *  it under the terms of the GNU General Public License as published by
80633 maechler 9 
 *  the Free Software Foundation; either version 3 of the License, or
2 r 10 
 *  (at your option) any later version.
11 
 *
12 
 *  This program is distributed in the hope that it will be useful,
13 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
14 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 
 *  GNU General Public License for more details.
16 
 *
42307 ripley 17 
 *  You should have received a copy of the GNU General Public License
18 
 *  along with this program; if not, a copy is available at
68947 ripley 19 
 *  https://www.R-project.org/Licenses/
2 r 20 
 */
21 
 
5187 hornik 22 
#ifdef HAVE_CONFIG_H
44031 ripley 23 
# include <config.h>
5187 hornik 24 
#endif
25 
 
79296 maechler 26 
#include <Defn.h>   // Rexp10  (et al)
51838 ripley 27 
#include <float.h>  /* for DBL_MAX */
11499 ripley 28 
#include <Graphics.h>
29 
#include <Print.h>
79296 maechler 30 
#include <Rmath.h> // for imax2
2 r 31 
 
59266 ripley 32 
/* used in graphics and grid */
80633 maechler 33 
SEXP CreateAtVector(double axp[], const double usr[], int nint, Rboolean logflag)
658 ihaka 34 
{
59266 ripley 35 
/*	Create an  'at = ...' vector for  axis(.)
2666 maechler 36 
 *	i.e., the vector of tick mark locations,
37 
 *	when none has been specified (= default).
38 
 *
21044 maechler 39 
 *	axp[0:2] = (x1, x2, nInt), where x1..x2 are the extreme tick marks
50323 maechler 40 
 *		   {unless in log case, where nInt \in {1,2,3 ; -1,-2,....}
41 
 *		    and the `nint' argument is used *instead*.}
80633 maechler 42 
 *
43 
 * only if(logflag && axp[2] >= 0)
44 
 *			usr[0:1] is used, additionally
45 
 *
2666 maechler 46 
 *	The resulting REAL vector must have length >= 1, ideally >= 2
47 
 */
987 maechler 48 
    SEXP at = R_NilValue;/* -Wall*/
80633 maechler 49 
    double dn, rng, small;
79296 maechler 50 
    int i, n;
80633 maechler 51 
    // "arbitrary" threshold: |delta_tick| / SMALL  is "barely visible" in plot
52 
#define SMALL_F 100.
21044 maechler 53 
    if (!logflag || axp[2] < 0) { /* --- linear axis --- Only use axp[] arg. */
59173 ripley 54 
	n = (int)(fabs(axp[2]) + 0.25);/* >= 0 */
8267 ripley 55 
	dn = imax2(1, n);
658 ihaka 56 
	rng = axp[1] - axp[0];
57 
	at = allocVector(REALSXP, n + 1);
80633 maechler 58 
	double a_i;
59 
	if(!R_FINITE(rng)) { // need to carefully work around overflow
60 
	    double at_ = axp[0]/dn; // 2021-07: "/dn" avoids overflow
61 
	    rng = axp[1]/dn - at_;
62 
	    small = fabs(rng)/SMALL_F;
63 
#ifdef DEBUG_axis
64 
	    REprintf("CreateAtVector(axp=(%g,%g, %g), log=F, diff(*)=Inf: at_=%g, rng=%g, small=%g\n",
65 
		     axp[0],axp[1], axp[2], at_, rng, small);
66 
#endif
67 
	    int n2 = n/2; // integer division
68 
	    for (i = 0; i <= n2; i++) { // from the left
69 
		a_i = axp[0] + i * rng;
70 
		// REprintf(" at[i=%2d]=%g\n", i+1, a_i);
71 
		REAL(at)[i] = (fabs(a_i) < small) ? 0. : a_i;
72 
	    }
73 
	    for (int i2 = 0; i2 < n-n2; i2++) { // from the right
74 
		i = n-i2; // for(i in n:k) where k = n-(n-n2-1) = n2+1
75 
		a_i = axp[1] - i2 * rng;
76 
		// REprintf(" at[i=%2d]=%g\n", i+1, a_i);
77 
		REAL(at)[i] = (fabs(a_i) < small) ? 0. : a_i;
78 
	    }
8552 maechler 79 
	}
80633 maechler 80 
	else { // rng is finite (normal case):
81 
	    small = fabs(rng)/SMALL_F/dn;
82 
	    for (i = 0; i <= n; i++) {
83 
		a_i = axp[0] + (i / dn) * rng;
84 
		REAL(at)[i] = (fabs(a_i) < small) ? 0. : a_i;
85 
	    }
86 
	}
658 ihaka 87 
    }
2669 maechler 88 
    else { /* ------ log axis ----- */
87891 ripley 89 
	bool reversed = false;
80633 maechler 90 
	double
91 
	    umin = usr[0],
92 
	    umax = usr[1];
59174 ripley 93 
	n = (int)(axp[2] + 0.5);
33335 maechler 94 
	/* {xy}axp[2] for 'log': GLpretty() [./graphics.c] sets
2666 maechler 95 
	   n < 0: very small scale ==> linear axis, above, or
96 
	   n = 1,2,3.  see switch() below */
79296 maechler 97 
#ifdef DEBUG_axis
98 
	REprintf("CreateAtVector(axp=(%g,%g,%g), usr=(%g,%g), _log_):",
99 
		 axp[0],axp[1],axp[2],  usr[0],usr[1]);
100 
#endif
34502 maechler 101 
	if (umin > umax) {
102 
	    reversed = (axp[0] > axp[1]);
103 
	    if (reversed) {
104 
		/* have *reversed* log axis -- whereas
105 
		 * the switch(n) { .. } below assumes *increasing* values
106 
		 * --> reverse axis direction here, and reverse back at end */
107 
		umin = usr[1];
108 
		umax = usr[0];
109 
		dn = axp[0]; axp[0] = axp[1]; axp[1] = dn;
110 
	    }
111 
	    else {
112 
		/* can the following still happen... ? */
113 
		warning("CreateAtVector \"log\"(from axis()): "
114 
			"usr[0] = %g > %g = usr[1] !", umin, umax);
115 
	    }
116 
	}
80633 maechler 117 
	/* allow a fuzz (iff we don't under-/over-flow) since we will do things like 0.2*dn >= umin */
118 
	dn = 1 - 1e-12; if(fabs(umin*dn) >     0.  ) umin *= dn;
119 
	dn = 1 + 1e-12; if(fabs(umax*dn) <= DBL_MAX) umax *= dn;
34502 maechler 120 
 
2690 maechler 121 
	dn = axp[0];
33405 maechler 122 
	if (dn < DBL_MIN) {/* was 1e-300; now seems too cautious */
123 
	    if (dn <= 0) /* real trouble (once for Solaris) later on */
124 
		error("CreateAtVector [log-axis()]: axp[0] = %g < 0!", dn);
75182 maechler 125 
	    else
126 
		warning("CreateAtVector [log-axis()]: small axp[0] = %g", dn);
33405 maechler 127 
	}
2669 maechler 128 
 
2666 maechler 129 
	/* You get the 3 cases below by
5507 ihaka 130 
	 *  for (y in 1e-5*c(1,2,8))  plot(y, log = "y")
2666 maechler 131 
	 */
658 ihaka 132 
	switch(n) {
2666 maechler 133 
	case 1: /* large range:	1	 * 10^k */
79296 maechler 134 
	{
59173 ripley 135 
	    i = (int)(floor(log10(axp[1])) - ceil(log10(axp[0])) + 0.25);
79296 maechler 136 
	    // want nint intervals, i.e. typically nint+1 breaks :
137 
	    int ne = i / nint;
138 
	    /* for nint breaks, i.e. typically nint-1 intervals, would be
139 
	     * ne = i / imax2(1, nint - 1);  *PLUS* replace s/nint/nint-1/ below !! */
50323 maechler 140 
#ifdef DEBUG_axis
79296 maechler 141 
	    REprintf(" .. case 1: umin,umax= %g,%g;\n  (nint=%d, ne=%d); ",
142 
		     umin, umax, nint, ne);
143 
	    if (ne < 1) {
144 
		REprintf("ne = %d <= 0 !!\n\t axp[0:1]=(%g,%g) ==> i = %d, nint = %d; ",
145 
			 ne, axp[0],axp[1], i, nint);
146 
	    }
50323 maechler 147 
#endif
79296 maechler 148 
	    double l10_max = log10(umax),
149 
		d0 = l10_max - log10(dn);
150 
#ifdef DEBUG_axis
151 
	    REprintf("exponent diff d0=%g\n", d0);
152 
#endif
153 
	    if(ne < 1) ne = 1;
154 
	    else // if ne is too large, i.e, the "final tick" is beyond umax, reduce it :
155 
		while(ne > 1 && nint*ne > d0) {
156 
		    ne--;
157 
#ifdef DEBUG_axis
158 
		    REprintf(" last > umax ==> ne--: ne=%d\n", ne);
159 
#endif
160 
		}
161 
	    int k = 1 + ne / 308; // >= 1, typically == 1.
162 
	    if(k > 1) {// i.e. ne > 308: 10^ne overflows; must split the multiplication
163 
		ne = k*(ne/k); // <= ne_{previous}
164 
#ifdef DEBUG_axis
165 
		REprintf(" original ne > 308: split in k=%d parts; new ne=%d\n", k,ne);
166 
#endif
167 
	    }
168 
	    /* Now, still in exponent-10 range: nint*ne <= d0 = l10_max - log10(dn)
169 
	     * If difference (=: d1) is "large", say > 3, increase the first at[] =: d0
170 
	     */
171 
	    double d1 = d0 - nint*ne; // >= 0
172 
#ifdef DEBUG_axis
173 
	    REprintf("expo.diff d0 - nint*ne =: d1=%g\n", d1);
174 
#endif
175 
	    d0 = dn;
80633 maechler 176 
 
177 
#define Large_D1 5
178 
//              === was '3' all up into R 4.1.0
179 
	    if(d1 > Large_D1) {
79296 maechler 180 
		d0 = dn * Rexp10(floor(d1/2));
181 
#ifdef DEBUG_axis
80633 maechler 182 
		REprintf("large d1 => d0 := dn * 10 ^ fl(d1/2) = dn * 10^%d = %g\n",
79296 maechler 183 
			 (int)floor(d1/2), d0);
184 
#endif
185 
	    }
186 
	    rng = Rexp10((double)ne/k); // = 10^(ne/k) >= 10
187 
	    n=0;
188 
	    dn=d0;
189 
	    while(dn < umax) {
190 
		for(int j=0; j < k; j++)
191 
		    dn *= rng;
658 ihaka 192 
		n++;
193 
	    }
79296 maechler 194 
#ifdef DEBUG_axis
80633 maechler 195 
	    REprintf(" rng:=10^(ne/(k=%d)) = %g => n=%d, final dn=%g\n", k, rng, n, dn);
79296 maechler 196 
#endif
5507 ihaka 197 
	    if (!n)
2666 maechler 198 
		error("log - axis(), 'at' creation, _LARGE_ range: "
33125 ripley 199 
		      "invalid {xy}axp or par; nint=%d\n"
5731 ripley 200 
		      "	 axp[0:1]=(%g,%g), usr[0:1]=(%g,%g); i=%d, ni=%d",
2666 maechler 201 
		      nint, axp[0],axp[1], umin,umax, i,ne);
658 ihaka 202 
	    at = allocVector(REALSXP, n);
79296 maechler 203 
	    dn=d0;
204 
	    for(int i=0; i < n; i++) {
205 
		REAL(at)[i] = dn;
206 
		for(int j=0; j < k; j++)
207 
		    dn *= rng;
658 ihaka 208 
	    }
209 
	    break;
79296 maechler 210 
	}
2666 maechler 211 
	case 2: /* medium range:  1, 5	  * 10^k */
658 ihaka 212 
	    n = 0;
5507 ihaka 213 
	    if (0.5 * dn >= umin) n++;
79296 maechler 214 
#ifdef DEBUG_axis
215 
	    REprintf(" .. case 2: (dn, umin,umax, n) = (%g, %g,%g, %d)\n",
216 
		     dn, umin, umax, n);
217 
#endif
5507 ihaka 218 
	    for (;;) {
69855 ripley 219 
		if (dn > umax) break;
220 
		n++;
221 
		if (5 * dn > umax) break;
222 
		n++;
2666 maechler 223 
		dn *= 10;
658 ihaka 224 
	    }
5507 ihaka 225 
	    if (!n)
2666 maechler 226 
		error("log - axis(), 'at' creation, _MEDIUM_ range: "
33125 ripley 227 
		      "invalid {xy}axp or par;\n"
5731 ripley 228 
		      "	 axp[0]= %g, usr[0:1]=(%g,%g)",
2666 maechler 229 
		      axp[0], umin,umax);
230 
 
658 ihaka 231 
	    at = allocVector(REALSXP, n);
232 
	    dn = axp[0];
233 
	    n = 0;
5507 ihaka 234 
	    if (0.5 * dn >= umin) REAL(at)[n++] = 0.5 * dn;
235 
	    for (;;) {
79260 ripley 236 
		if (dn > umax) break;
237 
		REAL(at)[n++] = dn;
238 
		if (5 * dn > umax) break;
239 
		REAL(at)[n++] = 5 * dn;
2666 maechler 240 
		dn *= 10;
658 ihaka 241 
	    }
242 
	    break;
2666 maechler 243 
 
244 
	case 3: /* small range:	 1,2,5,10 * 10^k */
658 ihaka 245 
	    n = 0;
5507 ihaka 246 
	    if (0.2 * dn >= umin) n++;
247 
	    if (0.5 * dn >= umin) n++;
248 
	    for (;;) {
69855 ripley 249 
		if (dn > umax) break;
250 
		n++;
251 
		if (2 * dn > umax) break;
252 
		n++;
253 
		if (5 * dn > umax) break;
254 
		n++;
2666 maechler 255 
		dn *= 10;
658 ihaka 256 
	    }
79296 maechler 257 
#ifdef DEBUG_axis
80633 maechler 258 
	    REprintf(" .. case 3: (umin,umax)-usr[*] = (%g, %g); n=%d, dn=%g\n",
259 
		     umin-usr[reversed? 1: 0],
260 
		     umax-usr[reversed? 0: 1], n, dn);
79296 maechler 261 
#endif
5507 ihaka 262 
	    if (!n)
2666 maechler 263 
		error("log - axis(), 'at' creation, _SMALL_ range: "
33125 ripley 264 
		      "invalid {xy}axp or par;\n"
5731 ripley 265 
		      "	 axp[0]= %g, usr[0:1]=(%g,%g)",
2666 maechler 266 
		      axp[0], umin,umax);
658 ihaka 267 
	    at = allocVector(REALSXP, n);
268 
	    dn = axp[0];
269 
	    n = 0;
5507 ihaka 270 
	    if (0.2 * dn >= umin) REAL(at)[n++] = 0.2 * dn;
271 
	    if (0.5 * dn >= umin) REAL(at)[n++] = 0.5 * dn;
272 
	    for (;;) {
79262 ripley 273 
		if (dn > umax) break;
274 
		REAL(at)[n++] = dn;
275 
		if (2 * dn > umax) break;
276 
		REAL(at)[n++] = 2 * dn;
277 
		if (5 * dn > umax) break;
278 
		REAL(at)[n++] = 5 * dn;
2666 maechler 279 
		dn *= 10;
658 ihaka 280 
	    }
281 
	    break;
2666 maechler 282 
	default:
33125 ripley 283 
	    error("log - axis(), 'at' creation: INVALID {xy}axp[3] = %g",
2666 maechler 284 
		  axp[2]);
658 ihaka 285 
	}
34502 maechler 286 
 
287 
	if (reversed) {/* reverse back again - last assignment was at[n++]= . */
288 
	    for (i = 0; i < n/2; i++) { /* swap( at[i], at[n-i-1] ) : */
289 
		dn = REAL(at)[i];
80633 maechler 290 
		     REAL(at)[i] = REAL(at)[n-i-1];
291 
		                   REAL(at)[n-i-1] = dn;
34502 maechler 292 
	    }
293 
	}
294 
    } /* linear / log */
658 ihaka 295 
    return at;
296 
}