是否可以对R中的optimize()函数进行矢量化处理，如果可以，该如何进行处理? [英] Is it possible to vectorize the optimize() function in R and if so how?

问题描述

我有一些我想加快速度的代码.在optimize()调用的函数中有一个瓶颈.我想传入参数的optimize()向量，在该向量中对向量的每个元素进行优化.这有意义吗?可以吗?

I have some code which I am trying to speed up. There is a bottle neck inside a function called by optimize(). I would like to pass into optimize() vectors of parameters where optim runs for each element of the vector. Does this make sense, can it be done and if so how?

我还应该指出，我在此问题中所指的优化瓶颈之一已在optimbreach中进行了说明，尽管我很想知道自己在做的任何事情效率低下.最初调用的主要函数是plotresults(). rundisplay()和displayProbEff()并不重要，在运行plotresults()时不会被调用.

I should also point out that one of the optimize bottlenecks I am referring to in this question is illustrated inside optimbreach, although I would love to know anything that I am doing that is inefficient. The main function to call initially is plotresults(). rundisplay() and displayProbEff() aren't important and won't be called when running plotresults().

library(scatterplot3d)
library(rgl)

plotresults = function()

{

r_start = 0.5
r_end = 1
r_step = 0.1
nsteps = (r_end-r_start)/r_step+1

riskstart = 0.01
riskend = 1
risksteps = 0.005

plotcols=NULL

for(i in 1:((r_end-r_start)/r_step+1))
  {
  plotcols = c(plotcols, rep(i, (riskend-riskstart)/risksteps+1))
  }

risk=NULL
rout=NULL
foregone=NULL
colour = 

for (r in seq(r_start,r_end,r_step))
  {
  output = species2(r, riskstart, riskend, risksteps)
  risk = c(risk, output$risk)
  rout = c(rout,rep(r, length(output$risk)))
  foregone = c(foregone, output$foregone)   
  #]mod = lm(results$foregone ~ poly(results$risk, 3))  
  }
r= rout  

plot3d(risk, r, foregone, col=plotcols, size=3)
}

species2 = function(r, riskstart,riskend,risksteps)

{
  #set how many years to run the model for when testing risk
  nYears = 3

  #setup the variables for risk values

  #init outputs
  #effort=NULL
  catch=numeric((riskend-riskstart)/risksteps)

  #read in the parameters for the 2 species
  SpeciesParams = as.matrix(read.csv("C:/Documents and Settings/localuser/My Documents/DefineIt/Code/R/speciesparams.csv", sep=",", header=T))

  # set the species parameter for species B to r (this enables us to get results for a range of r's
  SpeciesParams[2,4]=r

  #Calc optimum catch for species A
  maxcatch = optimcatch(SpeciesParams[1,])[[2]]

  #Create vector for risk
  risk = seq(riskstart, riskend, risksteps)

  i=1
    #browser() 
  #Determine the level of effort that subjects species B to a given level of risk
  for(risk in seq(riskstart,riskend,risksteps))
    {
    #calc the effort that would subject vulnerable species (B) to each level of risk
    effort=optimbreach(nYears,risk, SpeciesParams[2,])[[1]]

    #calc the catch of target species (A) for effort calculated in previous line
    catch[i]=maxcatch-meancatch(effort, SpeciesParams[1,])

    #if(catch[length(catch)]>200){browser()}
    i=i+1
    }

    #plot(x=seq(riskstart,riskend,risksteps), y=catch, xlab = "risk", ylab="Foregone Yield")
    output = list(risk=seq(riskstart,riskend,risksteps), foregone=catch)
    return (output)



}


RunDisplay = function(nTimeSteps, e) # e=effort, nTimeSteps=number of timesteps per run
#Runs the model and outputs the biomass and catches as two timeseries plots
  {

  output = RunModelnTimeSteps(nTimeSteps, stddevr_global, r_global, Init_B_global, k_global, q_global, e, 0)

  split.screen(c(1,2))
  screen(1)
  plot(output[[2]], ylim=c(0,k_global), type="l", ylab="Biomass", xlab="Timestep")
  screen(2)
  plot(output[[3]], ylim=c(0,k_global), type="l", ylab="Catch", xlab="TimeStep")

  }

displayProbEff = function()
#A simple subroutine that graphs the change in probability of breaching Bpa for different efforts
  {

  prob_lessBpa = NULL
  for(x in seq(0,2,0.02))
    {
    prob_lessBpa=c(prob_lessBpa,calcprob(x,3,Bpa_global))
    }
  plot(x=seq(0,2,0.02), y=prob_lessBpa, ylab="Probability", xlab="Effort")

  }

optimbreach = function(nTimeSteps, opt_prob, spec_params)
#Calculates the level of effort that meets management objective
#nTimeSteps is length of time to run simulation while checking not falling beneath Bpa
#opt_prob is the maximum probability of going beneath Bpa in a given time period (see previous line)
#that is acceptable by the management objective

  {

  diffprob = function(e, params, spec_params)
    {
    nTimeSteps=params[1]
    opt_prob=params[2]
    prob_breach = calcprob(e, nTimeSteps, spec_params)
    (prob_breach-opt_prob)^2
    }


  params=c(nTimeSteps, opt_prob)  
  a=optimise(diffprob, c(0,1.5), params, spec_params, maximum=FALSE)
  #browser()
  a

  }

optimcatch = function(species_params)
#Calculates the effort that gives the optimum catch
  {

  optimise(meancatch, c(0,2), species_params, maximum=TRUE)

  }

meancatch = function(ef, species_params)
  {
  output=RunModelnTimeSteps(nTimeSteps=1000, species_params, ef)  
  mean(output[[3]])
  }

calcmeans = function(effort)
#Calculates the mean catch and biomass over 1000 timesteps
  {

  output = RunModelnTimeSteps(1000, species_params, e)
  print(paste("Mean Biomass = ", mean(output[[2]])))
  print(paste("Mean Catch = ", mean(output[[3]])))

  }

RunModelnTimeSteps = function(nTimeSteps, sparams, e)
#Runs the surplus production model for nTimeSteps
  {
  #indexes in sparams for various parameters
  #Bpa = 2
  #sd_ = 3
  #r = 4
  #B = 5
  #k = 6
  #q_ = 7 

  P = 0                         #initial production

  Breached = FALSE              #This records whether the biomass fell beneath Bpa
  Biomass = numeric(nTimeSteps)                #biomass
  Catch = numeric(nTimeSteps)          #catch

  Biomass[1]=sparams[5]
  Catch[1]=2

  for (x in 2:nTimeSteps)
    {
    ModOutput = RunTimeStep(sparams[3], sparams[4], Biomass[x-1], sparams[6], sparams[7], e)  #Run the model for 1 timestep

    Biomass[x] = ModOutput[1]  #Get the biomass at end of timestep

    Catch[x] = ModOutput[2]  #Get the catch for this timestep
    #print(Biomass[x])
    #print(sparams[2])
    if (Biomass[x]<=sparams[2])  #Check if biomass has fallen beneath Bpa
      {                            
      Breached=TRUE
      if(Biomass[x]<0){Biomass[x]=0}
      }      
    }
  #if(Bpa==600){browser()   }
  list(Breached, Biomass, Catch)

  }


RunTimeStep = function(sd_, r, B, k, q_, e)
#Calculates one timestep of the model
#outputs the biomass[1] at the end of the timestep and catches[2] through timestep
  {

  change_r = rnorm(n=1, mean=0, sd=sd_) #random change in r
  P = (r + change_r) * B * (1-B/k)  #calc production
  Catch = q_ * e * B  #calc catch
  B = B + P - q_*e*B  #calc results biomass
  c(B,Catch)  #output biomass at end of timestep and catch

  }


calcprob = function(e, nTimeSteps, species_params) # e=effort, Bpa=precautionary biomass limit, nTimeSteps=number of timesteps per run
#Calculates the probability of going beneath Bpa within nTimeSteps
  {

  nIterations = 200            #number of times the simulation is run
  nFails = 0                    #Counts the number of times the biomass goes benath Bpa

  nFails = RunModelnTimeSteps2(nTimeSteps, species_params, e, nIterations)

  nFails/nIterations

  }


RunModelnTimeSteps2 = function(nTimeSteps, sparams, e, nIterations)
#Runs the surplus production model for nTimeSteps
  {
  #indexes in sparams for various parameters
  #Bpa = 2
  #sd_ = 3
  #r = 4
  #B = 5
  #k = 6
  #q_ = 7 

  P = 0                         #initial production

  Breached = rep(F,nIterations)              #This records whether the biomass fell beneath Bpa
  Biomass = matrix(nrow=nTimeSteps, ncol=nIterations)                #biomass
  Catch = matrix(nrow=nTimeSteps, ncol=nIterations)          #catch
  randomvalues = matrix(rnorm(nTimeSteps*nIterations,0,sparams[3]),ncol=nIterations)

  Biomass[1,]=sparams[5]
  Catch[1,]=2

  #ModOutput = RunTimeStep2(sparams[3], sparams[4], Biomass[x-1,], sparams[6], sparams[7], e)

  for (x in 2:nTimeSteps)
    {
    #ModOutput = RunTimeStep2(sparams[3], sparams[4], Biomass[x-1,], sparams[6], sparams[7], e) #Run the model for 1 timestep

    P = (sparams[4] + randomvalues[x,]) * Biomass[x-1] * (1-Biomass[x-1]/sparams[6])  #calc production
    Catch[x,] = sparams[7] * e * Biomass[x-1]
    Biomass[x,] = Biomass[x-1,] + P - Catch[x,]

    Breached = ifelse (Biomass[x,]<=sparams[2], T, Breached)
    ifelse (Biomass[x,]<0,0,Biomass[x,]) 

    }

  #if(Bpa==600){browser()   }

  sum(Breached)
  }

推荐答案

这是使用parallel获得速度提升的一种方法，但这不是向量化.如果您发布了一些代码，那么人们可能会注意到其他加快它的方法.只需重新编写代码即可将其提供给并行应用功能之一，如下所示:

Here's a way you can get speed gains using parallel, but it's not vectorization. If you posted some code, then folks might notice other ways to speed it up. Just rework your code to be able to feed it to one of the parallel apply functions, like this:

    # some data:
    set.seed(1)
    a <- 1:100
    b <- rnorm(100, mean = 3, sd = 2) + runif(100, 3, 4) * a

    # some starts:
    pars <- c(a = 1, b = 1)

    # a function to optimize
    FunctionToOptimize <- function(pars, a, b){
        sum(((pars[[1]] + pars[[2]] * a) - b) ^ 2)
    }

    optim(pars,FunctionToOptimize,a = a, b = b)$par

    # and in parallel?
    # a matrix holding your parameters:
    parsMat <- matrix(sample(2:5, 100, replace = TRUE), ncol = 2, nrow = 50)

    # a function that includes a call to optim()
    FunctionToOptimizePar <- function(Pars, a, b){

        FunctionToOptimize <- function(pars, a, b){
            sum(((pars[[1]] + pars[[2]] * a) - b) ^ 2)
        } 

        optim(Pars, FunctionToOptimize, a = a, b = b)$par
    }

    library(parallel)
    # (I just have 2 cores)
    cl <- makeCluster(2)
    # will spit out the optimized parameters in the same order as given in parsMat
    someresults <- matrix(parRapply(cl, parsMat, FunctionToOptimizePar, a = a, b = b), ncol = 2, byrow = TRUE)
    stopCluster(cl)

    head(someresults)

这篇关于是否可以对R中的optimize()函数进行矢量化处理，如果可以，该如何进行处理?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！