[seqfan] Re: Chewing an old bone

[David Wilson]

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/
>