R type conversion expression() function()

Question

I've been trying to write a program in R that implements Newton's method. I've been mostly successful, but there are two little snags that have been bothering me. Here's my code:

Newton<-function(f,f.,guess){
    #f <- readline(prompt="Function? ")
    #f. <- readline(prompt="Derivative? ")
    #guess <- as.numeric(readline(prompt="Guess? "))
    a <- rep(NA, length=1000)
    a[1] <- guess
    a[2] <- a[1] - f(a[1]) / f.(a[1])
    for(i in 2:length(a)){
        if(a[i] == a[i-1]){
           break
        } 
        else{
           a[i+1] <- a[i] - f(a[i]) / f.(a[i])
        }
    }   
    a <- a[complete.cases(a)]
    return(a)
}

1. I can't get R to recognize the functions f and f. if I try using readline() to prompt for user input. I get the error "Error in Newton() : could not find function "f."" However, if I comment out the readlines (as above), define f and f. beforehand, then everything works fine.

2. I've been trying to make R calculate the derivative of a function. The problem is that the class object with which R can take symbolic derivatives is expression(), but I want to take the derivative of a function() and have it give me a function(). In short, I'm having trouble with type conversion between expression() and function().

I have an ugly but effective solution for going from function() to expression(). Given a function f, D(body(f)[[2]],"x") will give the derivative of f. However, this output is an expression(), and I haven't been able to turn it back into a function(). Do I need to use eval() or something? I've tried subsetting, but to no avail. For instance:

g <- expression(sin(x))
g[[1]]
sin(x)
f <- function(x){g[[1]]}
f(0)
sin(x)

when what I want is f(0) = 0 since sin(0) = 0.

EDIT: Thanks all! Here's my new code:

Newton<-function(f,f.,guess){
    g<-readline(prompt="Function? ")
    g<-parse(text=g)
    g.<-D(g,"x")
    f<-function(x){eval(g[[1]])}
    f.<-function(x){eval(g.)}
    guess<-as.numeric(readline(prompt="Guess? "))
    a<-rep(NA, length=1000)
    a[1]<-guess
    a[2]<-a[1]-f(a[1])/f.(a[1])
    for(i in 2:length(a)){
        if(a[i]==a[i-1]){break
        }else{
        a[i+1]<-a[i]-f(a[i])/f.(a[i])
        }
    }   
a<-a[complete.cases(a)]
#a<-a[1:(length(a)-1)]
return(a)
}

Illation

Answer 1

1. This first problem arises because readline reads in a text string, whereas what you need is an expression. You can use parse() to convert the text string to an expression:

   f <-readline(prompt="Function? ")
   sin(x)
   f
   # [1] "sin(x)"
   
   f <- parse(text = f)
   f
   # expression(sin(x))
   
   g <- D(f, "x")
   g
   # cos(x)
   

2. To pass in values for the arguments in the function call in the expression (whew!), you can eval() it in an environment containing the supplied values. Nicely, R will allow you to supply those values in a list supplied to the envir= argument of eval():

   > eval(f, envir=list(x=0))
   # [1] 0
   

Panaggio

Thank you! Is there a better alternative to readline()? I didn't know about parse(), though I had tried as.expression unsuccessfully. Also, is there a better option for going from function() to expression() other than body()? – Illation

readline() is the right function for taking typed input from users, and there is no better way than parse() to convert text to expressions. Not sure what you mean by going from function() to expression(), though. Is function() a function call, or a function definition? – Feature

I mean things that are in the function() class, and which are specifically mathematical functions. For example, say I have the function f<-function(x){3x*sin(x/3) + 2*cos(4*x)} and I want to find its derivative. The only way I know to do this is to use body() to obtain an expression, and then use D(body(f)[[2]],"x"). – Illation

OK. There are a few problems here. (1) comparing class(f) and class(body(f)) shows that body(f) is a function call, not a function. (That's OK of course, since for D() you really want an object of class expression, not function.) (2) Unless I'm mistaken you want to be taking the derivative of something that looks more like body(f) than body(f)[[2]]. (3) You probably want to use as.expression to do this: D(as.expression(body(f)), "x"). Does that get you closer to where you want to be? – Feature

I think the subsetting (i.e. using the [[2]]) is necessary, otherwise the curly brackets cause problems. For instance: > f<-function(x){3*x*sin(x/3) + 2*cos(4*x)} > class(body(f)) [1] "{" > class(f) [1] "function" > class(body(f)) [1] "{" > class(body(f)[[2]]) [1] "call" And, trying out your example, I get an error: > D(as.expression(body(f)), "x") Error in D(as.expression(body(f)), "x") : Function '{' is not in the derivatives table I study math, so I apologize for not being too familiar with the terminology. – Illation

Got'cha. I had typed f <- function(x) 3*x*sin(x/3) + 2*cos(4*x) (i.e., without the brackets), for which D(as.expression(body(f)), "x") does work. – Feature

Haha, and today I learned that I don't need to surround the function's body with brackets. Thanks for all your help. – Illation

Answer 2

Josh has answered your question

For part 2 you could have used

g <- expression( sin(x) )

g[[1]]
# sin(x)

f <- function(x){ eval( g[[1]] ) }

f(0)
# [1] 0
f(pi/6)
# [1] 0.5

Lineate

Thank you! One more wrinkle: what if I want to get rid of the expression variable (g in your example)? Really, I'd like to do f <- expression(sin(x)) f <- function(x){eval(f[[1]])} but that creates a circular dependence since in 'eval(f[[1]])' the 'f[[1]]' isn't replaced by what's actuallly in f (namely 'sin(x)'). – Illation

I am no longer clear what you are asking for. Perhaps your are looking for something like this which does not have g: f <- function(x){ eval( expression( sin(x) )[[1]] ) }; f(pi/6) or perhaps something which is not a function x <- pi/4 ; eval( expression (sin(x))) – Lineate

Exactly--I want to eliminate g. I'd like to be able to extract the contents of g and then get rid of g. Using my example: '> g<-expression(sin(x)) > f<-function(x){eval(g[[1]])} > f function(x){eval(g[[1]])}' When what I really want is '> f function(x){sin(x)}' In other words, I want to remove the dependence of f on g. Does that make sense? – Illation

Answer 3

BTW, having recently written a toy which calculates fractal patterns based on root convergence of Newton's method in the complex plane, I can recommend you toss in some code like the following (where the main function's argument list includes "func" and "varname" ).

func<- gsub(varname, 'zvar', func)
    funcderiv<- try( D(parse(text=func), 'zvar') )
    if(class(funcderiv) == 'try-error') stop("Can't calculate derivative")

If you're more cautious, you could include a an argument "funcderiv" , and wrap my code in

if(missing(funcderiv)){blah blah}

Ahh, why not: here's my complete function for all to use and enjoy:-)

# build Newton-Raphson fractal
#define: f(z)  the convergence per Newton's method is 
# zn+1 = zn - f(zn)/f'(zn)
#record which root each starting z0 converges to, 
# and to get even nicer coloring, record the number of iterations to get there.
# Inputs:
#   func: character string, including the variable. E.g., 'x+ 2*x^2' or 'sin(x)'
#   varname: character string indicating the variable name
#   zreal: vector(preferably) of Re(z)
#   zim: vector of Im(z)
#   rootprec: convergence precision for the NewtonRaphson algorithm
#   maxiter: safety switch, maximum iterations, after which throw an error
#
nrfrac<-function(func='z^5 - 1 ', varname = 'z', zreal= seq(-5,5,by=.1), zim, rootprec=1.0e-5, maxiter=1e4, drawplot=T, drawiterplot=F, ...) {
    zreal<-as.vector(zreal)
    if(missing(zim)) zim <- as.vector(zreal)
# precalculate F/F' 
    # check for differentiability (in R's capability)
    # and make sure to get the correct variable name into the function
    func<- gsub(varname, 'zvar', func)
    funcderiv<- try( D(parse(text=func), 'zvar') )
    if(class(funcderiv) == 'try-error') stop("Can't calculate derivative")  
# Interesting "feature" of deparse : default is to limit each string to 60 or64
# chars.  Need to avoid that here.  Doubt I'd ever see a derivative w/ more
# than 500 chars, the max allowed by deparse. To do it right, 
# need sum(nchar(funcderiv)) as width, and even then need to do some sort of
# paste(deparse(...),collapse='') to get a single string
    nrfunc <- paste(text='(',func,')/(',deparse(funcderiv, width=500),')', collapse='')
# first arg to outer()  will give rows
# Stupid Bug: I need to REVERSE zim to get proper axis orientation
    zstart<- outer(rev(zim*1i), zreal, "+")
    zindex <- 1:(length(zreal)*length(zim))
    zvec <- data.frame(zdata=as.vector(zstart), zindex=zindex,     itermap=rep(0,length(zindex)), badroot=rep(0,length(zindex)), rooterr=rep(0,length(zindex)) )

#initialize data.frame for zout.  
    zout=data.frame(zdata=rep(NA,length(zstart)), zindex=rep(NA,length(zindex)),     itermap=rep(0,length(zindex)), badroot=rep(0,length(zindex)), rooterr=rep(0,length(zindex)))
    # a value for rounding purposes later on; yes it works for  rootprec >1 
    logprec <-  -floor(log10(rootprec))
    newtparam <- function(zvar) {}
    body(newtparam)[2]  <- parse(text=paste('newz<-', nrfunc, collapse=''))
    body(newtparam)[3] <- parse(text=paste('return(invisible(newz))'))
    iter <- 1
    zold <- zvec  # save zvec so I can return original values
    zoutind <- 1 #initialize location to write solved values
    while (iter <= maxiter & length(zold$zdata)>0 ) {
        zold$rooterr <- newtparam(zold$zdata)
        zold$zdata <- zold$zdata - zold$rooterr
        rooterr <- abs(zold$rooterr)
        zold$badroot[!is.finite(rooterr)] <- 1
        zold$zdata[!is.finite(rooterr)] <- NA
# what if solvind = FFFFFFF? -- can't write 'nothing' to zout
        solvind <- (zold$badroot >0 | rooterr<rootprec)
            if( sum(solvind)>0 ) zout[zoutind:(zoutind-1+sum(solvind)),] <- zold[solvind,]
    #update zout index to next 'empty' row
        zoutind<-zoutind + sum(solvind)
# update the iter count for remaining elements:
        zold$itermap <- iter
# and reduce the size of the matrix being fed back to loop
        zold<-zold[!solvind,]
        iter <- iter +1
    # just wonder if a gc() call here would make any difference
# wow -- it sure does
        gc()
    }  # end of while
# Now, there may be some nonconverged values, so:
#  badroot[]  is set to 2  to distinguish from Inf/NaN locations
        if( zoutind < length(zindex) ) { # there are nonconverged values
#  fill the remaining rows, i.e. zout.index:length(zindex)
            zout[(zoutind:length(zindex)),] <- zold # all of it
            zold$badroot[] <- 2 # yes this is safe for length(badroot)==0
            zold$zdata[]<-NA #keeps nonconverged values from messing up results
            }
#  be sure to properly re-order everything...
    zout<-zout[order(zout$zindex),]
    zout$zdata <- complex(re=round(Re(zout$zdata),logprec), im=round(Im(zout$zdata),logprec) )
    rootvec <- factor(as.vector(zout$zdata), labels=c(1:length(unique(na.omit(as.vector(zout$zdata))))))
    #convert from character, too!
    rootIDmap<-matrix(as.numeric(rootvec), nr=length(zim))
# to colorize very simply:  
    if(drawplot) {
             colorvec<-rainbow(length(unique(as.vector(rootIDmap))))
        imagemat<-rootIDmap
        imagemat[,]<-colorvec[imagemat]  #now has color strings
        dev.new()
# all '...' arguments used to set up plot
        plot(range((zreal)),range((zim)), t='n',xlab='real',ylab='imaginary',... ) 
        rasterImage(imagemat, range(zreal)[1], range(zim)[1], range(zreal)[2], range(zim)[2], interp=F)     
        }

    outs <- list(rootIDmap=rootIDmap, zvec=zvec, zout=zout, nrfunc=nrfunc)
    return(invisible(outs))
}

Kano

Thanks for posting your code. I haven't read all the way through it yet, but I'm sure I can improve my program by emulating your code. – Illation