/* * R : A Computer Language for Statistical Data Analysis * Copyright (C) 1997--2021 The R Core Team * Copyright (C) 2002--2009 The R Foundation * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, a copy is available at * https://www.R-project.org/Licenses/ */ #ifdef HAVE_CONFIG_H # include #endif #include // Rexp10 (et al) #include /* for DBL_MAX */ #include #include #include // for imax2 /* used in graphics and grid */ SEXP CreateAtVector(double axp[], const double usr[], int nint, Rboolean logflag) { /* Create an 'at = ...' vector for axis(.) * i.e., the vector of tick mark locations, * when none has been specified (= default). * * axp[0:2] = (x1, x2, nInt), where x1..x2 are the extreme tick marks * {unless in log case, where nInt \in {1,2,3 ; -1,-2,....} * and the `nint' argument is used *instead*.} * * only if(logflag && axp[2] >= 0) * usr[0:1] is used, additionally * * The resulting REAL vector must have length >= 1, ideally >= 2 */ SEXP at = R_NilValue;/* -Wall*/ double dn, rng, small; int i, n; // "arbitrary" threshold: |delta_tick| / SMALL is "barely visible" in plot #define SMALL_F 100. if (!logflag || axp[2] < 0) { /* --- linear axis --- Only use axp[] arg. */ n = (int)(fabs(axp[2]) + 0.25);/* >= 0 */ dn = imax2(1, n); rng = axp[1] - axp[0]; at = allocVector(REALSXP, n + 1); double a_i; if(!R_FINITE(rng)) { // need to carefully work around overflow double at_ = axp[0]/dn; // 2021-07: "/dn" avoids overflow rng = axp[1]/dn - at_; small = fabs(rng)/SMALL_F; #ifdef DEBUG_axis REprintf("CreateAtVector(axp=(%g,%g, %g), log=F, diff(*)=Inf: at_=%g, rng=%g, small=%g\n", axp[0],axp[1], axp[2], at_, rng, small); #endif int n2 = n/2; // integer division for (i = 0; i <= n2; i++) { // from the left a_i = axp[0] + i * rng; // REprintf(" at[i=%2d]=%g\n", i+1, a_i); REAL(at)[i] = (fabs(a_i) < small) ? 0. : a_i; } for (int i2 = 0; i2 < n-n2; i2++) { // from the right i = n-i2; // for(i in n:k) where k = n-(n-n2-1) = n2+1 a_i = axp[1] - i2 * rng; // REprintf(" at[i=%2d]=%g\n", i+1, a_i); REAL(at)[i] = (fabs(a_i) < small) ? 0. : a_i; } } else { // rng is finite (normal case): small = fabs(rng)/SMALL_F/dn; for (i = 0; i <= n; i++) { a_i = axp[0] + (i / dn) * rng; REAL(at)[i] = (fabs(a_i) < small) ? 0. : a_i; } } } else { /* ------ log axis ----- */ Rboolean reversed = FALSE; double umin = usr[0], umax = usr[1]; n = (int)(axp[2] + 0.5); /* {xy}axp[2] for 'log': GLpretty() [./graphics.c] sets n < 0: very small scale ==> linear axis, above, or n = 1,2,3. see switch() below */ #ifdef DEBUG_axis REprintf("CreateAtVector(axp=(%g,%g,%g), usr=(%g,%g), _log_):", axp[0],axp[1],axp[2], usr[0],usr[1]); #endif if (umin > umax) { reversed = (axp[0] > axp[1]); if (reversed) { /* have *reversed* log axis -- whereas * the switch(n) { .. } below assumes *increasing* values * --> reverse axis direction here, and reverse back at end */ umin = usr[1]; umax = usr[0]; dn = axp[0]; axp[0] = axp[1]; axp[1] = dn; } else { /* can the following still happen... ? */ warning("CreateAtVector \"log\"(from axis()): " "usr[0] = %g > %g = usr[1] !", umin, umax); } } /* allow a fuzz (iff we don't under-/over-flow) since we will do things like 0.2*dn >= umin */ dn = 1 - 1e-12; if(fabs(umin*dn) > 0. ) umin *= dn; dn = 1 + 1e-12; if(fabs(umax*dn) <= DBL_MAX) umax *= dn; dn = axp[0]; if (dn < DBL_MIN) {/* was 1e-300; now seems too cautious */ if (dn <= 0) /* real trouble (once for Solaris) later on */ error("CreateAtVector [log-axis()]: axp[0] = %g < 0!", dn); else warning("CreateAtVector [log-axis()]: small axp[0] = %g", dn); } /* You get the 3 cases below by * for (y in 1e-5*c(1,2,8)) plot(y, log = "y") */ switch(n) { case 1: /* large range: 1 * 10^k */ { i = (int)(floor(log10(axp[1])) - ceil(log10(axp[0])) + 0.25); // want nint intervals, i.e. typically nint+1 breaks : int ne = i / nint; /* for nint breaks, i.e. typically nint-1 intervals, would be * ne = i / imax2(1, nint - 1); *PLUS* replace s/nint/nint-1/ below !! */ #ifdef DEBUG_axis REprintf(" .. case 1: umin,umax= %g,%g;\n (nint=%d, ne=%d); ", umin, umax, nint, ne); if (ne < 1) { REprintf("ne = %d <= 0 !!\n\t axp[0:1]=(%g,%g) ==> i = %d, nint = %d; ", ne, axp[0],axp[1], i, nint); } #endif double l10_max = log10(umax), d0 = l10_max - log10(dn); #ifdef DEBUG_axis REprintf("exponent diff d0=%g\n", d0); #endif if(ne < 1) ne = 1; else // if ne is too large, i.e, the "final tick" is beyond umax, reduce it : while(ne > 1 && nint*ne > d0) { ne--; #ifdef DEBUG_axis REprintf(" last > umax ==> ne--: ne=%d\n", ne); #endif } int k = 1 + ne / 308; // >= 1, typically == 1. if(k > 1) {// i.e. ne > 308: 10^ne overflows; must split the multiplication ne = k*(ne/k); // <= ne_{previous} #ifdef DEBUG_axis REprintf(" original ne > 308: split in k=%d parts; new ne=%d\n", k,ne); #endif } /* Now, still in exponent-10 range: nint*ne <= d0 = l10_max - log10(dn) * If difference (=: d1) is "large", say > 3, increase the first at[] =: d0 */ double d1 = d0 - nint*ne; // >= 0 #ifdef DEBUG_axis REprintf("expo.diff d0 - nint*ne =: d1=%g\n", d1); #endif d0 = dn; #define Large_D1 5 // === was '3' all up into R 4.1.0 if(d1 > Large_D1) { d0 = dn * Rexp10(floor(d1/2)); #ifdef DEBUG_axis REprintf("large d1 => d0 := dn * 10 ^ fl(d1/2) = dn * 10^%d = %g\n", (int)floor(d1/2), d0); #endif } rng = Rexp10((double)ne/k); // = 10^(ne/k) >= 10 n=0; dn=d0; while(dn < umax) { for(int j=0; j < k; j++) dn *= rng; n++; } #ifdef DEBUG_axis REprintf(" rng:=10^(ne/(k=%d)) = %g => n=%d, final dn=%g\n", k, rng, n, dn); #endif if (!n) error("log - axis(), 'at' creation, _LARGE_ range: " "invalid {xy}axp or par; nint=%d\n" " axp[0:1]=(%g,%g), usr[0:1]=(%g,%g); i=%d, ni=%d", nint, axp[0],axp[1], umin,umax, i,ne); at = allocVector(REALSXP, n); dn=d0; for(int i=0; i < n; i++) { REAL(at)[i] = dn; for(int j=0; j < k; j++) dn *= rng; } break; } case 2: /* medium range: 1, 5 * 10^k */ n = 0; if (0.5 * dn >= umin) n++; #ifdef DEBUG_axis REprintf(" .. case 2: (dn, umin,umax, n) = (%g, %g,%g, %d)\n", dn, umin, umax, n); #endif for (;;) { if (dn > umax) break; n++; if (5 * dn > umax) break; n++; dn *= 10; } if (!n) error("log - axis(), 'at' creation, _MEDIUM_ range: " "invalid {xy}axp or par;\n" " axp[0]= %g, usr[0:1]=(%g,%g)", axp[0], umin,umax); at = allocVector(REALSXP, n); dn = axp[0]; n = 0; if (0.5 * dn >= umin) REAL(at)[n++] = 0.5 * dn; for (;;) { if (dn > umax) break; REAL(at)[n++] = dn; if (5 * dn > umax) break; REAL(at)[n++] = 5 * dn; dn *= 10; } break; case 3: /* small range: 1,2,5,10 * 10^k */ n = 0; if (0.2 * dn >= umin) n++; if (0.5 * dn >= umin) n++; for (;;) { if (dn > umax) break; n++; if (2 * dn > umax) break; n++; if (5 * dn > umax) break; n++; dn *= 10; } #ifdef DEBUG_axis REprintf(" .. case 3: (umin,umax)-usr[*] = (%g, %g); n=%d, dn=%g\n", umin-usr[reversed? 1: 0], umax-usr[reversed? 0: 1], n, dn); #endif if (!n) error("log - axis(), 'at' creation, _SMALL_ range: " "invalid {xy}axp or par;\n" " axp[0]= %g, usr[0:1]=(%g,%g)", axp[0], umin,umax); at = allocVector(REALSXP, n); dn = axp[0]; n = 0; if (0.2 * dn >= umin) REAL(at)[n++] = 0.2 * dn; if (0.5 * dn >= umin) REAL(at)[n++] = 0.5 * dn; for (;;) { if (dn > umax) break; REAL(at)[n++] = dn; if (2 * dn > umax) break; REAL(at)[n++] = 2 * dn; if (5 * dn > umax) break; REAL(at)[n++] = 5 * dn; dn *= 10; } break; default: error("log - axis(), 'at' creation: INVALID {xy}axp[3] = %g", axp[2]); } if (reversed) {/* reverse back again - last assignment was at[n++]= . */ for (i = 0; i < n/2; i++) { /* swap( at[i], at[n-i-1] ) : */ dn = REAL(at)[i]; REAL(at)[i] = REAL(at)[n-i-1]; REAL(at)[n-i-1] = dn; } } } /* linear / log */ return at; }