Error in stats::integrate(newf, vi.from, vi.to[1], rel.tol = .tol) : maximum number of subdivisions reached

[chrisvwn]

I am getting this error when trying to perform an integration. Is it possible to find out the reason for this? Can I modify the tolerance that is passed on to integrate?

The formula I am trying to integrate is:

t1 <- antiD(as.formula(7.48 * (1.41172190749093 * 0.000571499999999999 * z/0.0028/(0.000571499999999999/(0.000571499999999999 + 
    0.00673091999999992) * d_) - (0.000571499999999999 * z/0.0028/(0.000571499999999999/(0.000571499999999999 + 
    0.00673091999999992) * d_))^2)/(1 + (1.41172190749093 - 2) * 
    0.000571499999999999 * z/0.0028/(0.000571499999999999/(0.000571499999999999 + 
    0.00673091999999992) * d_)) ~ d_))

t1(d_ = 76, z = 5.947891)

Error in stats::integrate(newf, vi.from, vi.to[1], rel.tol = .tol) :    maximum number of subdivisions reached

Any help would be appreciated.

[dtkaplan]

I think there is an error in the way the code was processed by the email. I translated your code into this form

t0 <- makeFun(7.48 (a b z/c/(b/(b + cc) x) - (b z/c/(b/(b + cc) x))^2) / (1 + (a - 2) b z/c/(b/(b + cc) * x) ) ~ x, a = 1.41172190749093, b = 0.000571499999999999, c = 0.0028, cc = 0.00673091999999992, z = 5.947891 )

t1 <- antiD(t0(x) ~ x)

First, I defined the function as t0() that's inside your antiD() call, then I called antiD() on the function t0().

There was a missing parenthesis in the formula. I needed to add one in to construct a proper expression, but of course I don't know if I added it in the place you intended. The one I added is right before the ~ .

With my added parentheses (which may make the rest of what I say irrelevant, if it's in the wrong place), the function has a singularity at zero and another near x = 9.125453. The numerical integration is choking at these singularities. A way to get around this is to integrate in pieces. Perhaps use the built-in integrator, which provides many convergence options, for integrating pieces around the singularity. I don't know if it will work, but first you should determine if I placed the parenthesis in the right place, Incidentally, there's no reason to call as.formula() with a formula as an argument. So I took it out.

Regards, Danny

On Mon, Oct 26, 2020 at 10:52 AM Chris Njuguna notifications@github.com wrote:

I am getting this error when trying to perform an integration. Is it possible to find out the reason for this? Can I modify the tolerance that is passed on to integrate?

The formula I am trying to integrate is:

t1 <- antiD(as.formula(7.48 (1.41172190749093 0.000571499999999999 z/0.0028/(0.000571499999999999/(0.000571499999999999 + 0.00673091999999992) d) - (0.000571499999999999 z/0.0028/(0.000571499999999999/(0.000571499999999999 + 0.00673091999999992) d))^2)/(1 + (1.41172190749093 - 2) 0.000571499999999999 z/0.0028/(0.000571499999999999/(0.000571499999999999 + 0.00673091999999992) * d)) ~ d))

t1(d_ = 76, z = 5.947891)

Error in stats::integrate(newf, vi.from, vi.to[1], rel.tol = .tol) : maximum number of subdivisions reached

Any help would be appreciated.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/ProjectMOSAIC/mosaicCalc/issues/5, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAECBLOFOBBODVD7IBCFVETSMWSLNANCNFSM4S7UP5OA .

[chrisvwn]

Thanks @dtkaplan

This syntax is a bit new to me. I am currently building the formula from variables using quote and substituting the values in which is why I have the as.formula. I think I may have mismatched the brackets in the process. The bracket you have added is in the right position.

So maybe it will be easier to explain the whole task. I am porting code from maple to R. In a particular function where the above code was extracted, we first integrate and then differentiate but with scenarios where we may have unresolved symbols. Different variables can be the target hence the crazy code (below) trying to avoid symbolic calculations and delay them for later when values are available. It works like a charm when different variables are the target except for two variables which appear in the integration itself. So the D(antiD()) chokes when using stats::integrate or cubature::cubintegrate. My only hope so far is mosaicCalc.

      fn0 <- function(g,q,d,s1,f1,a1,th1,th2,X,Y){
        g_ <- g
        q_ <- q
        d_ <- d
        s1_ <- s1
        f1_ <- f1
        a1_ <- a1
        th1_ <- th1
        th2_ <- th2

        fc_ <- if(!inherits(try(x<-a2*f1_,TRUE), "try-error")) x else quote(a2*f1_)
        e1_ <- if(!inherits(try(x<-a1_*(f1_*10)^(1/3)/10,TRUE), "try-error")) x else quote(a1_*(f1_*10)^(1/3)/10)

        if(0.7*(x_f1_k*10)^0.31/1000<0.0028){
          e2_ <- quote(0.7*f1_^0.31/1000)
        }else{
          e2_ <- 0.0028
        }

        k_ <- if(!inherits(try(x<-1.05*e1_*e2_/f1_,TRUE), "try-error")) x else quote(1.05*e1_*e2_/f1_)
        x_ <- if(!inherits(try(x<-x1/(x1+x2)*d_,TRUE), "try-error")) x else quote(x1/(x1+x2)*d_)

        #s1_ <- function(z)fc_*(k_*x1*z/e2_/x_-(x1*z/e2_/x_)^2)/(1+(k_-2)*x1*z/e2_/x_)
        s1_ <- quote(fc_*(k_*x1*z/e2_/x_-(x1*z/e2_/x_)^2)/(1+(k_-2)*x1*z/e2_/x_) ~ X)
        s1_ <- do.call(substitute, list(s1_, list(fc_=fc_,f1_=f1_,e1_=e1_,k_=k_,x1=x1,x2=x2,e2_=e2_,x_=x_,X=ensym(X))))
        s1_ <- do.call(substitute, list(s1_, list(fc_=fc_,f1_=f1_,e1_=e1_,k_=k_,x1=x1,x2=x2,e2_=e2_,x_=x_)))
        s1_ <- do.call(substitute, list(s1_, list(fc_=fc_,f1_=f1_,e1_=e1_,k_=k_,x1=x1,x2=x2,e2_=e2_,x_=x_)))

        s1A_ <- quote(z* (fc_*(k_*x1*z/e2_/x_-(x1*z/e2_/x_)^2)/(1+(k_-2)*x1*z/e2_/x_)) ~ X)
        s1A_ <- do.call(substitute, list(s1A_, list(fc_=fc_,f1_=f1_,e1_=e1_,k_=k_,x1=x1,x2=x2,e2_=e2_,x_=x_,X=ensym(X))))
        s1A_ <- do.call(substitute, list(s1A_, list(fc_=fc_,f1_=f1_,e1_=e1_,k_=k_,x1=x1,x2=x2,e2_=e2_,x_=x_)))
        s1A_ <- do.call(substitute, list(s1A_, list(fc_=fc_,f1_=f1_,e1_=e1_,k_=k_,x1=x1,x2=x2,e2_=e2_,x_=x_)))

        #commented out but kept to show how I was using cubintegrate before running into the symbol problem
        #Fc_ <- if(!inherits(try(x<-b*cubintegrate(s1_,0,x_)$integral,TRUE), "try-error")) x else quote(antiD(as.formula(s1_)))
        fns1_ <- antiD(as.formula(s1_), force.numeric = TRUE)
        Fc_ <- str2lang(paste0("fns1_(", X, " = ", X, ", z = z)"))

        #commented out but kept to show how I was using cubintegrate before running into the symbol problem
        #a_ <- if(!inherits(try(x<-x_-1/Fc_*b*cubintegrate(function(z)s1_(z)*z,0,x_)$integral,TRUE), "try-error")) x else quote(x_-1/Fc_*b*cubintegrate(function(z)z*function(z)fc_*(k_*x1*z/e2_/x_-(x1*z/e2_/x_)^2)/(1+(k_-2)*x1*z/e2_/x_),0,x_)$integral)
        fnA_ <- antiD(as.formula(s1A_))
        a_ <- str2lang(paste0("fnA_(",X, " = ", X, ", z = z)"))

        f2_ <- if(!inherits(try(x<-rhol*b*lp_d*s1_,TRUE),'try-error')) x else quote(rhol*b*lp_d*s1_)

        gv_ <- expression(th2_*(f2_*d_-Fc_*a_)-th1_*(g1*g_*l^2/8+g2*q_*l^2/8) ~ X)
        gv_ <- do.call(substitute, list(gv_[[1]], list(a_=a_,b=b,d_=d_,e2_=e2_,f2_=f2_,Fc_=Fc_,f2_=f2_,g_=g_,g1=g1,g2=g2,k_=k_,l=l,lp_d=lp_d,q_=q_,th2_=th2_,th1_=th1_,x_=x_,x1=x1,x2=x2,X=ensym(X))))
        gv_ <- do.call(substitute, list(gv_, list(a_=a_,b=b,d_=d_,e2_=e2_,f2_=f2_,Fc_=Fc_,f2_=f2_,g_=g_,g1=g1,g2=g2,k_=k_,l=l,lp_d=lp_d,q_=q_,th2_=th2_,th1_=th1_,x_=x_,x1=x1,x2=x2)))
        gv_ <- do.call(substitute, list(gv_, list(a_=a_,b=b,d_=d_,e2_=e2_,f2_=f2_,Fc_=Fc_,f2_=f2_,g_=g_,g1=g1,g2=g2,k_=k_,l=l,lp_d=lp_d,q_=q_,th2_=th2_,th1_=th1_,x_=x_,x1=x1,x2=x2)))
        ab <- D(as.formula(gv_))
        wab <- eval(str2lang(paste0("ab(", X, " = ", Y, ", z = ", x_, ")")))
      }

A couple of calls look like this:

This one where g_ is the target works easy:

        An_g <- fn0(g=quote(g_), q=0.02148409, d=76, s1=165.6479, f1=8.8, a1=9500, th1=1, th2=1, X="g_", Y=0.03069155);

This one is more problematic since d_ is the target of the integral:

        An_d <- fn0(g=0.03069155, q=0.02148409, d=quote(d_), s1=165.6479, f1=8.8, a1=9500, th1=1, th2=1, X="d_", Y=76);

I hope this is not too crazy. I can post the maple code as well if maybe there is an easier way to translate this kind of thing to R.