The R Project SVN R

Rev

Rev 79263 | Rev 80633 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 79263 Rev 79296
Line 21... Line 21...
21 
 
 21 
 
22 
#ifdef HAVE_CONFIG_H
 22 
#ifdef HAVE_CONFIG_H
23 
# include <config.h>
 23 
# include <config.h>
24 
#endif
 24 
#endif
25 
 
 25 
 
26 
#include <Defn.h>
 26 
#include <Defn.h>   // Rexp10  (et al)
27 
#include <float.h>  /* for DBL_MAX */
 27 
#include <float.h>  /* for DBL_MAX */
28 
#include <Graphics.h>
 28 
#include <Graphics.h>
29 
#include <Print.h>
 29 
#include <Print.h>
30 
#include <Rmath.h> // for Rexp10, imax2
 30 
#include <Rmath.h> // for imax2
31 
 
 31 
 
32 
/* used in graphics and grid */
 32 
/* used in graphics and grid */
33 
SEXP CreateAtVector(double *axp, double *usr, int nint, Rboolean logflag)
 33 
SEXP CreateAtVector(double *axp, double *usr, int nint, Rboolean logflag)
34 
{
 34 
{
35 
/*	Create an  'at = ...' vector for  axis(.)
 35 
/*	Create an  'at = ...' vector for  axis(.)
Line 42... Line 42...
42 
 
 42 
 
43 
 *	The resulting REAL vector must have length >= 1, ideally >= 2
 43 
 *	The resulting REAL vector must have length >= 1, ideally >= 2
44 
 */
 44 
 */
45 
    SEXP at = R_NilValue;/* -Wall*/
 45 
    SEXP at = R_NilValue;/* -Wall*/
46 
    double umin, umax, dn, rng, small;
 46 
    double umin, umax, dn, rng, small;
47 
    int i, n, ne;
 47 
    int i, n;
48 
    if (!logflag || axp[2] < 0) { /* --- linear axis --- Only use axp[] arg. */
 48 
    if (!logflag || axp[2] < 0) { /* --- linear axis --- Only use axp[] arg. */
49 
	n = (int)(fabs(axp[2]) + 0.25);/* >= 0 */
 49 
	n = (int)(fabs(axp[2]) + 0.25);/* >= 0 */
50 
	dn = imax2(1, n);
 50 
	dn = imax2(1, n);
51 
	rng = axp[1] - axp[0];
 51 
	rng = axp[1] - axp[0];
52 
	small = fabs(rng)/(100.*dn);
 52 
	small = fabs(rng)/(100.*dn);
Line 64... Line 64...
64 
	/* {xy}axp[2] for 'log': GLpretty() [./graphics.c] sets
 64 
	/* {xy}axp[2] for 'log': GLpretty() [./graphics.c] sets
65 
	   n < 0: very small scale ==> linear axis, above, or
 65 
	   n < 0: very small scale ==> linear axis, above, or
66 
	   n = 1,2,3.  see switch() below */
 66 
	   n = 1,2,3.  see switch() below */
67 
	umin = usr[0];
 67 
	umin = usr[0];
68 
	umax = usr[1];
 68 
	umax = usr[1];
- 
 
 69 
#ifdef DEBUG_axis
- 
 
 70 
	REprintf("CreateAtVector(axp=(%g,%g,%g), usr=(%g,%g), _log_):",
- 
 
 71 
		 axp[0],axp[1],axp[2],  usr[0],usr[1]);
- 
 
 72 
#endif
69 
	if (umin > umax) {
 73 
	if (umin > umax) {
70 
	    reversed = (axp[0] > axp[1]);
 74 
	    reversed = (axp[0] > axp[1]);
71 
	    if (reversed) {
 75 
	    if (reversed) {
72 
		/* have *reversed* log axis -- whereas
 76 
		/* have *reversed* log axis -- whereas
73 
		 * the switch(n) { .. } below assumes *increasing* values
 77 
		 * the switch(n) { .. } below assumes *increasing* values
Line 97... Line 101...
97 
	/* You get the 3 cases below by
 101 
	/* You get the 3 cases below by
98 
	 *  for (y in 1e-5*c(1,2,8))  plot(y, log = "y")
 102 
	 *  for (y in 1e-5*c(1,2,8))  plot(y, log = "y")
99 
	 */
 103 
	 */
100 
	switch(n) {
 104 
	switch(n) {
101 
	case 1: /* large range:	1	 * 10^k */
 105 
	case 1: /* large range:	1	 * 10^k */
- 
 
 106 
	{
102 
	    i = (int)(floor(log10(axp[1])) - ceil(log10(axp[0])) + 0.25);
 107 
	    i = (int)(floor(log10(axp[1])) - ceil(log10(axp[0])) + 0.25);
- 
 
 108 
	    // want nint intervals, i.e. typically nint+1 breaks :
103 
	    ne = i / nint + 1;
 109 
	    int ne = i / nint;
- 
 
 110 
	    /* for nint breaks, i.e. typically nint-1 intervals, would be
- 
 
 111 
	     * ne = i / imax2(1, nint - 1);  *PLUS* replace s/nint/nint-1/ below !! */
104 
#ifdef DEBUG_axis
 112 
#ifdef DEBUG_axis
105 
	    REprintf("CreateAtVector [log-axis(), case 1]: (nint, ne) = (%d,%d)\n",
 113 
	    REprintf(" .. case 1: umin,umax= %g,%g;\n  (nint=%d, ne=%d); ",
106 
		     nint, ne);
 114 
		     umin, umax, nint, ne);
- 
 
 115 
	    if (ne < 1) {
- 
 
 116 
		REprintf("ne = %d <= 0 !!\n\t axp[0:1]=(%g,%g) ==> i = %d, nint = %d; ",
- 
 
 117 
			 ne, axp[0],axp[1], i, nint);
- 
 
 118 
	    }
107 
#endif
 119 
#endif
- 
 
 120 
	    double l10_max = log10(umax),
- 
 
 121 
		d0 = l10_max - log10(dn);
- 
 
 122 
#ifdef DEBUG_axis
- 
 
 123 
	    REprintf("exponent diff d0=%g\n", d0);
- 
 
 124 
#endif
108 
	    if (ne < 1)
 125 
	    if(ne < 1) ne = 1;
- 
 
 126 
	    else // if ne is too large, i.e, the "final tick" is beyond umax, reduce it :
- 
 
 127 
		while(ne > 1 && nint*ne > d0) {
- 
 
 128 
		    ne--;
- 
 
 129 
#ifdef DEBUG_axis
109 
		error("log - axis(), 'at' creation, _LARGE_ range: "
 130 
		    REprintf(" last > umax ==> ne--: ne=%d\n", ne);
- 
 
 131 
#endif
- 
 
 132 
		}
- 
 
 133 
	    int k = 1 + ne / 308; // >= 1, typically == 1.
- 
 
 134 
	    if(k > 1) {// i.e. ne > 308: 10^ne overflows; must split the multiplication
110 
		      "ne = %d <= 0 !!\n"
 135 
		ne = k*(ne/k); // <= ne_{previous}
- 
 
 136 
#ifdef DEBUG_axis
- 
 
 137 
		REprintf(" original ne > 308: split in k=%d parts; new ne=%d\n", k,ne);
- 
 
 138 
#endif
- 
 
 139 
	    }
- 
 
 140 
	    /* Now, still in exponent-10 range: nint*ne <= d0 = l10_max - log10(dn)
- 
 
 141 
	     * If difference (=: d1) is "large", say > 3, increase the first at[] =: d0
- 
 
 142 
	     */
- 
 
 143 
	    double d1 = d0 - nint*ne; // >= 0
- 
 
 144 
#ifdef DEBUG_axis
111 
		      "\t axp[0:1]=(%g,%g) ==> i = %d;	nint = %d",
 145 
	    REprintf("expo.diff d0 - nint*ne =: d1=%g\n", d1);
- 
 
 146 
#endif
- 
 
 147 
	    d0 = dn;
- 
 
 148 
	    if(d1 > 3) {
112 
		      ne, axp[0],axp[1], i, nint);
 149 
		d0 = dn * Rexp10(floor(d1/2));
- 
 
 150 
#ifdef DEBUG_axis
- 
 
 151 
		REprintf(" d0 := dn * 10 ^ fl(d1/2) = dn * 10^%d = %g\n",
- 
 
 152 
			 (int)floor(d1/2), d0);
- 
 
 153 
#endif
- 
 
 154 
	    }
113 
	    rng = Rexp10((double)ne); /* >= 10 */
 155 
	    rng = Rexp10((double)ne/k); // = 10^(ne/k) >= 10
114 
	    n = 0;
 156 
	    n=0;
- 
 
 157 
	    dn=d0;
115 
	    while (dn < umax) {
 158 
	    while(dn < umax) {
- 
 
 159 
		for(int j=0; j < k; j++)
- 
 
 160 
		    dn *= rng;
116 
		n++;
 161 
		n++;
117 
		dn *= rng;
 - 
 
118 
	    }
 162 
	    }
- 
 
 163 
#ifdef DEBUG_axis
- 
 
 164 
	    REprintf(" rng:=10^(ne/%d) = %g => n=%d, final dn=%g\n", k, rng, n, dn);
- 
 
 165 
#endif
119 
	    if (!n)
 166 
	    if (!n)
120 
		error("log - axis(), 'at' creation, _LARGE_ range: "
 167 
		error("log - axis(), 'at' creation, _LARGE_ range: "
121 
		      "invalid {xy}axp or par; nint=%d\n"
 168 
		      "invalid {xy}axp or par; nint=%d\n"
122 
		      "	 axp[0:1]=(%g,%g), usr[0:1]=(%g,%g); i=%d, ni=%d",
 169 
		      "	 axp[0:1]=(%g,%g), usr[0:1]=(%g,%g); i=%d, ni=%d",
123 
		      nint, axp[0],axp[1], umin,umax, i,ne);
 170 
		      nint, axp[0],axp[1], umin,umax, i,ne);
124 
	    at = allocVector(REALSXP, n);
 171 
	    at = allocVector(REALSXP, n);
125 
	    dn = axp[0];
 - 
 
126 
	    n = 0;
 172 
	    dn=d0;
127 
	    while (dn < umax) {
 173 
	    for(int i=0; i < n; i++) {
128 
		REAL(at)[n++] = dn;
 174 
		REAL(at)[i] = dn;
- 
 
 175 
		for(int j=0; j < k; j++)
129 
		dn *= rng;
 176 
		    dn *= rng;
130 
	    }
 177 
	    }
131 
	    break;
 178 
	    break;
132 
 
 179 
	}
133 
	case 2: /* medium range:  1, 5	  * 10^k */
 180 
	case 2: /* medium range:  1, 5	  * 10^k */
134 
	    n = 0;
 181 
	    n = 0;
135 
	    if (0.5 * dn >= umin) n++;
 182 
	    if (0.5 * dn >= umin) n++;
- 
 
 183 
#ifdef DEBUG_axis
- 
 
 184 
	    REprintf(" .. case 2: (dn, umin,umax, n) = (%g, %g,%g, %d)\n",
- 
 
 185 
		     dn, umin, umax, n);
- 
 
 186 
#endif
136 
	    for (;;) {
 187 
	    for (;;) {
137 
		if (dn > umax) break;
 188 
		if (dn > umax) break;
138 
		n++;
 189 
		n++;
139 
		if (5 * dn > umax) break;
 190 
		if (5 * dn > umax) break;
140 
		n++;
 191 
		n++;
Line 170... Line 221...
170 
		n++;
 221 
		n++;
171 
		if (5 * dn > umax) break;
 222 
		if (5 * dn > umax) break;
172 
		n++;
 223 
		n++;
173 
		dn *= 10;
 224 
		dn *= 10;
174 
	    }
 225 
	    }
- 
 
 226 
#ifdef DEBUG_axis
- 
 
 227 
	    REprintf(" .. case 3: (umin,umax) = (%g, %g); n=%d, dn=%g\n",
- 
 
 228 
		     umin, umax, n, dn);
- 
 
 229 
#endif
175 
	    if (!n)
 230 
	    if (!n)
176 
		error("log - axis(), 'at' creation, _SMALL_ range: "
 231 
		error("log - axis(), 'at' creation, _SMALL_ range: "
177 
		      "invalid {xy}axp or par;\n"
 232 
		      "invalid {xy}axp or par;\n"
178 
		      "	 axp[0]= %g, usr[0:1]=(%g,%g)",
 233 
		      "	 axp[0]= %g, usr[0:1]=(%g,%g)",
179 
		      axp[0], umin,umax);
 234 
		      axp[0], umin,umax);