[seqfan] Re: Chewing an old bone
David Wilson
davidwwilson at comcast.net
Fri Oct 14 05:18:26 CEST 2011
I downloaded R and took a crack at it.
There are three types of programming languages
- Ones in which it is easy to say what you want (Perl).
- Ones in which it is difficult to say what you want (C++).
- Ones that don't let you say what you want (R, apparently).
I downloaded an R IDE off the web and fired it up. I was able to tweak
the program
below to get it to generate a graph from the OEIS b-file for the Fibs.
At any rate, I tried the obvious and replaced log10(abs(v)+1) with my
asinh expression.
Amazingly, the asinh() function was defined, but it failed for large
negative numbers.
I suspect at bottom asinh is just ln(x+sqrt(x*x+1)), and for large
enough negative values,
x = sqrt(x*x+1) to double precision, resulting in an ln(0).
I also tried writing a more stable qlog function, but R has some
endearing features
like vectors as the only atomic type, and no concept of local variables.
I could not
get a simple function, or even a simple "if" statement, to work.
So I'm waving the white flag. Perhaps we could rewrite this program in
some slightly
less user-hostile language (e.g. Perl, C++, FORTRAN, TECO, Coleco Adam
BASIC, INTERCAL,
or perhaps MC68000 assembler).
On 10/13/2011 9:29 AM, Russ Cox wrote:
> This suggestion sounds interesting. The best way to explore it
> would be if someone could implement it in R, and then we can try
> both for a little while.
>
> I am including the R program that generates the current graphs below.
> The input file read by read.table is the standard b file format,
> without any comments or blank lines.
>
> If someone can adapt this program to generate these alternate
> plots, then I will be happy to set it up so that we can see both.
>
> Thanks.
> Russ
>
>
> # R CMD BATCH --no-restore --no-save --vanilla -- script.r
>
> infile = "b123.txt"
> outfile = "b123.png"
> seqname = "Name Of Sequence"
>
> x = read.table(infile)
>
> # Remove infinite values
> lgl = is.finite(x[,2])
> x = cbind(x[,1][lgl], x[,2][lgl])
>
> png(file=outfile, width=600, height=800)
> par(mar=c(5,5,2,2)+0.1, mfcol=c(2,1))
>
> k = min(200, nrow(x))
> xsub = x[1:k,]
> plot (xsub,
> ylim = c(min(0, min(xsub[,2])), max(xsub[,2])),
> xlab = "n",
> ylab=paste(seqname, "(n)", sep=""),
> bty="l",
> type="n",
> main = paste("Pin plot of", seqname))
> segments(xsub[,1], rep(0, k), xsub[,1], xsub[,2])
>
> # Log the data for the scatterplot?
>
> logplot = T
> tval = 100.
> if (nrow(x)> 10) {
> fv = fivenum(abs(x[,2]))
> if (all(is.finite(fv))) {
> if (((fv[5]-fv[4]) / (fv[2]-fv[1]+1))< tval) {
> logplot = F
> }
> }
> }
>
> if (logplot) {
> v = x[,2]
> if (all((v+1)>= 0)) {
> plot(x[,1], log10(v+1), pch="+", cex=.5,
> xlab = "n",
> ylab=paste("log(", seqname, "(n)+1)", sep=""),
> main=paste("Scatterplot of log(", seqname, "(n)+1)", sep=""),
> bty="l")
> } else {
> plot(x[,1], log10(abs(v)+1), pch="+", cex=.5,
> xlab = "n",
> ylab=paste("log(|", seqname, "(n)|+1)", sep=""),
> main=paste("Scatterplot of log(", seqname, "(n))", sep=""),
> bty="l")
> }
> } else {
> plot(x, pch="+", cex=.5,
> xlab = "n",
> ylab = paste(seqname, "(n)", sep=""),
> main = paste("Scatterplot of ", seqname, "(n)", sep=""),
> bty="l")
> }
>
> dev.off()
>
> _______________________________________________
>
> Seqfan Mailing list - http://list.seqfan.eu/
>
More information about the SeqFan
mailing list