
Hi there

I thought this would be of particular interest to people using 'optim' 
functions and perhaps people involved with R development.

I've been beaten down by R trying to get it to perform an optimization on a 
mass-balance model.  I've written the same program in excel, and using the 
'solver' function, it comes up with an answer for my variables (p, ACT, 
which I've assigned to q in R) that gives a solution to the function "f" in 
about 3 seconds, with a value of the function around 0.0004. R, on the 
other hand, appears to get stuck in local minima, and spits back an 
approximation that is close the the p, ACT values excel does, but not 
nearly precise enough for my needs, and not nearly as precise as excel, and 
it takes about 3 minutes.  Also, the solution for the value it returns for 
the function is about 8 orders of magnitude greater than the excel version, 
so I can't really say the function is approximating zero.  I was able to 
determine this using a  "trace" command on function f, which is listed below.

This is very likely due to the fact that I've made some coding error along 
the way, or have done something else wrong, but I don't know.  Either way, 
I am shocked and surprised that a program like excel is outperforming 
R.  I've attached my command file and the dataset "temp.dat" at the bottom 
of this e-mail for anyone who would like to fiddle around with it, and if 
you come up with something, PLEASE let me know- In the meantime, I've got 
to start fiddling with excel and figuring out how to automate the solver 
calculation.....

Briefly, the point of the program is to approximate the model output from 
an iterative calculation, Wtmod and Hgtmod, to user-specified endpoints Wt 
and Hgt, by seeking the optimal values of p, ACT involved in the iterative 
process.

Also, if your interested in recent correspondence that explains the point 
of the program a bit, and how the function ties in to the iterative 
process, search the R help forum for e-mails entitled "[R] problem with 
coding for 'optim' in R".  Thanks also to Roger Peng and numerous others 
for helping me get this far.

The whole point of me doing this in R was because it's supposed to be 
spectacularly fast at automating complex loops, but seems to be falling 
short for this application.  Hopefully it's something wrong with my coding 
and not with R itself.

Mike

R COMMAND FILE:

####################################
#    perch.R                       #
# Hewett and Johnson bioenergetics #
# model combined with              #
# Trudel MMBM to estimate          #
# Consumption in perch in R code   #
# Execute with                     #
# R --vanilla < perch.R > perch.out#
####################################

#USER INPUT BELOW

#Weight at time 0
Wo<- 9.2

#Hg concentration at time 0 (ugHg/g wet weight)
Hgo<- 0.08

#Weight at time t
Wt<- 32.2

#Hg concentration at time t (ugHg/g wet weight)
Hgt<- 0.110

#Prey methylmercury concentration (as constant)
Hgp<- 0.033

#Prey caloric value (as constant)
Pc<- 800

#Energy density of fish (as constant, calories)
Ef <- 1000

#Maturity status, 0=immature, 1=mature
Mat<- 0

#Sex, 1=male, 2=female
Sex<- 1

#USER INPUT ABOVE

#Bioenergetics parameters for perch
CA <- 0.25
CB <- 0.73  #same as 1+(-0.27)- convert g/g/d to g/d * Pc to get cal/d
CQ <- 2.3
CTO <- 23
CTM <- 28
Zc<- (log(CQ))*(CTM-CTO)
Yc<- (log(CQ))*(CTM-CTO+2)
Xc<- ((Zc^2)*(1+(1+40/Yc)^0.5)^2)/400

RA <- 34.992  #0.0108*3240 cal/g 02, converting weight of 02 to cal
RB <- 0.8   #same as 1+(-0.2) see above...
RQ <- 2.1
RTO <- 28
RTM <- 33
Za <- (log(RQ))*(RTM-RTO)
Ya<- (log(RQ))*(RTM-RTO+2)
Xa<- ((Za^2)*(1+(1+40/Ya)^0.5)^2)/400

S <- 0.172

FA <- 0.158
FB <- -0.222
FG <- 0.631

UA<- 0.0253
UB<- 0.58
UG<- -0.299

#Mass balance model parameters
EA <- 0.002938
EB <- -0.2
EQ <- 0.066
a <- 0.8

#Specifying sex-specific parameters

GSI<- NULL

if (Sex==1) GSI<-0.05 else
if (Sex==2) GSI<-0.17

# Define margin of error functions
#merror <- function(phat,M,alpha) # (1-alpha)*100% merror for a proportion
#    {
#    z <- qnorm(1-alpha/2)
#    merror <- z * sqrt(phat*(1-phat)/M)  # M is (Monte Carlo) sample size
#    merror
#    }

#Bring in temp file

temper <- scan("temp.dat", na.strings = ".", list(Day=0, jday=0, Temp=0))

Day<-temper$Day ; jday<-temper$jday ; Temp<-temper$Temp ;

temp<- cbind (Day, jday, Temp)
#Day = number of days modelled, jday=julian day, Temp = daily avg. temp.
#temp [,2]

Vc<-(CTM-(temp[,3]))/(CTM-CTO)
Vr<-(RTM-(temp[,3]))/(RTM-RTO)

comp<- cbind (Day, jday, Temp, Vc, Vr)

#comp

bio<-matrix(NA, ncol=13, nrow=length(Day))
W<-NULL
C<-NULL
ASMR<-NULL
SMR<-NULL
A<-NULL
F<-NULL
U<-NULL
SDA<-NULL
Gr<-NULL
Hg<-NULL
Ed<-NULL
GHg<-NULL
K<-NULL
Expegk<-NULL
EGK<-NULL
p<-NULL
ACT<-NULL

#starting values for p, ACT
p <- 1 #  0.558626306252032 #solution set for p, ACT from excel 'solver' f'n
ACT <- 2 #  1.66764519286918

q<-c(p,ACT)

#specify sttarting values
#q0<-c(p = 1, ACT = 1)

#introduce function to solve
f <- function (q)
{


M<- length(Day) #number of days iterated

for (i in 1:M)
{

#Bioenergetics model

if (Day[i]==1) W[i] <- Wo else
if (jday[i]==121 && Mat==1) W[i] <- (W[i-1]-(W[i-1]*GSI*1.2)) else
W[i] <- (W[i-1]+(Gr[i-1]/Ef))
#W

#W<-Wo

C[i]<- q[1]*CA*(W[i]^CB)*((comp[i,4])^Xc)*(exp(Xc*(1-(comp[i,4]))))*Pc

ASMR[i]<- q[2]*RA*(W[i]^RB)*((comp[i,5])^Xa)*(exp(Xa*(1-(comp[i,5]))))

SMR[i]<- (ASMR[i]/q[2])

A[i]<- (ASMR[i]-SMR[i])

F[i]<- (FA*((comp[i,3])^FB)*(exp(FG*p))*C[i])

U[i]<- (UA*((comp[i,3])^UB)*(exp(UG*p))*(C[i]-F[i]))

SDA[i]<- (S*(C[i]-F[i]))

Gr[i]<- (C[i]-(ASMR[i]+F[i]+U[i]+SDA[i]))

#Trudel MMBM

if (Day[i]==1) Hg[i] <- Hgo else Hg[i] <- 
a*Hgp*(C[i-1]/Pc/W[i-1])/EGK[i-1]*(1-Expegk[i-1])+(Hg[i-1]*Expegk[i-1])

Ed[i]<- EA*(W[i]^EB)*(exp(EQ*(comp[i,3])))

GHg[i] <- Gr[i]/Ef/W[i]

if (Sex==1) 
K[i]<-(((0.1681*(10^(1.3324+(0.000453*Hg[i])))/1000)/Hg[i])*GSI)/M else
if (Sex==2) K[i]<-(((0.1500*(10^(0.8840+(0.000903*Hg[i])))/1000)/Hg[i])*GSI)/M
# = dw/ww conversion * gonad ~ body conc'n function(ng/g) / convert to ug/g
# then express as Q times GSI gives K / M gives daily K

EGK[i] <- (Ed[i] + GHg[i] + (K[i]*Mat))

Expegk[i] <- exp(-1*EGK[i])

bio<- cbind(W, C, ASMR, SMR, A, F, U, SDA, Gr, Ed, GHg, EGK, Hg)

}

#warnings()

dimnames (bio) <-list(NULL, c("W", "C", "ASMR", "SMR", "A", "F", "U", 
"SDA", "Gr", "Ed", "GHg", "EGK", "Hg"))

bioday<-cbind(jday, W, C, ASMR, SMR, A, F, U, SDA, Gr, Ed, GHg, EGK, Hg)

dimnames (bioday) <-list(NULL, c("jday", "W", "C", "ASMR", "SMR", "A", "F", 
"U", "SDA", "Gr", "Ed", "GHg", "EGK", "Hg"))

#bioday

Wtmod<- bioday [length(W),2]
Wtmod

Hgtmod<- bioday [length(Hg),14]
Hgtmod

q

f <- 1000000000*(((((Wt-Wtmod)^2)/Wt) + (((Hgt-Hgtmod)^2)/Hgt))^2) ; f

#warnings()

#write.table (bioday, file = "perch.csv", append = FALSE, sep=",", na = NA, 
col.names = TRUE)



#nlm(f,c(1,1))
}

optim(q, f, method = "L-BFGS-B",
         lower = c(0.2, 1), upper=c(2, 3),
         control = list(fnscale = 0.001))


TRACE FUNCTION USED TO DETERMINE WHERE R IS GETTING STUCK;

  trace("f", quote(print(q)), at = 1, print = FALSE)

o <- optim(c(1,2), f, method = "L-BFGS-B", lower = c(0.2,1), upper = c(2,3))

DATA FOR TEMP.DAT:

1       153     9.4
2       154     9.6
3       155     9.8
4       156     10
5       157     10.2
6       158     10.4
7       159     10.6
8       160     10.8
9       161     11
10      162     11.2
11      163     11.4
12      164     11.6
13      165     11.8
14      166     12
15      167     12.3
16      168     12.5
17      169     12.7
18      170     12.9
19      171     13.1
20      172     13.4
21      173     13.6
22      174     13.8
23      175     14
24      176     14.2
25      177     14.5
26      178     14.7
27      179     14.9
28      180     15.1
29      181     15.4
30      182     15.6
31      183     15.8
32      184     16
33      185     16.2
34      186     16.5
35      187     16.7
36      188     16.9
37      189     17.1
38      190     17.3
39      191     17.5
40      192     17.7
41      193     17.9
42      194     18.1
43      195     18.3
44      196     18.5
45      197     18.7
46      198     18.9
47      199     19
48      200     19.2
49      201     19.4
50      202     19.5
51      203     19.7
52      204     19.9
53      205     20
54      206     20.2
55      207     20.3
56      208     20.4
57      209     20.5
58      210     20.7
59      211     20.8
60      212     20.9
61      213     21
62      214     21.1
63      215     21.2
64      216     21.3
65      217     21.3
66      218     21.4
67      219     21.5
68      220     21.5
69      221     21.6
70      222     21.6
71      223     21.6
72      224     21.7
73      225     21.7
74      226     21.7
75      227     21.7
76      228     21.7
77      229     21.7
78      230     21.7
79      231     21.6
80      232     21.6
81      233     21.6
82      234     21.5
83      235     21.5
84      236     21.4
85      237     21.3
86      238     21.3
87      239     21.2
88      240     21.1
89      241     21
90      242     20.9
91      243     20.8
92      244     20.7
93      245     20.5
94      246     20.4
95      247     20.3
96      248     20.2
97      249     20
98      250     19.9
99      251     19.7
100     252     19.5
101     253     19.4
102     254     19.2
103     255     19
104     256     18.9
105     257     18.7
106     258     18.5
107     259     18.3
108     260     18.1
109     261     17.9
110     262     17.7
111     263     17.5
112     264     17.3
113     265     17.1
114     266     16.9
115     267     16.7
116     268     16.5
117     269     16.2
118     270     16
119     271     15.8
120     272     15.6
121     273     15.4
122     274     15.1
123     275     14.9
124     276     14.7
125     277     14.5
126     278     14.2
127     279     14
128     280     13.8
129     281     13.6
130     282     13.4
131     283     13.1
132     284     12.9
133     285     12.7
134     286     12.5
135     287     12.3
136     288     12
137     289     11.8
138     290     11.6
139     291     11.4
140     292     11.2
141     293     11
142     294     10.8
143     295     10.6
144     296     10.4
145     297     10.2
146     298     10
147     299     9.8
148     300     9.6
149     301     9.4
150     302     9.3
151     303     9.1
152     304     8.9
153     305     8.7
154     306     8.6
155     307     8.4
156     308     8.2
157     309     8.1
158     310     7.9
159     311     7.8
160     312     7.6
161     313     7.5
162     314     7.3
163     315     7.2
164     316     7
165     317     6.9
166     318     6.8
167     319     6.7
168     320     6.5
169     321     6.4
170     322     6.3
171     323     6.2
172     324     6.1
173     325     6
174     326     5.8
175     327     5.7
176     328     5.6
177     329     5.5
178     330     5.5
179     331     5.4
180     332     5.3
181     333     5.2
182     334     5.1
183     335     5
184     336     5
185     337     4.9
186     338     4.8
187     339     4.7
188     340     4.7
189     341     4.6
190     342     4.5
191     343     4.5
192     344     4.4
193     345     4.4
194     346     4.3
195     347     4.3
196     348     4.2
197     349     4.2
198     350     4.1
199     351     4.1
200     352     4
201     353     4
202     354     4
203     355     3.9
204     356     3.9
205     357     3.8
206     358     3.8
207     359     3.8
208     360     3.8
209     361     3.7
210     362     3.7
211     363     3.7
212     364     3.6
213     365     3.6
214     366     3.6
215     1       3.2
216     2       3.2
217     3       3.2
218     4       3.2
219     5       3.2
220     6       3.2
221     7       3.2
222     8       3.2
223     9       3.2
224     10      3.2
225     11      3.2
226     12      3.2
227     13      3.2
228     14      3.2
229     15      3.2
230     16      3.2
231     17      3.2
232     18      3.2
233     19      3.2
234     20      3.2
235     21      3.2
236     22      3.2
237     23      3.2
238     24      3.2
239     25      3.2
240     26      3.2
241     27      3.2
242     28      3.2
243     29      3.2
244     30      3.2
245     31      3.2
246     32      3.2
247     33      3.2
248     34      3.2
249     35      3.2
250     36      3.2
251     37      3.2
252     38      3.2
253     39      3.2
254     40      3.2
255     41      3.2
256     42      3.2
257     43      3.2
258     44      3.2
259     45      3.2
260     46      3.2
261     47      3.2
262     48      3.2
263     49      3.2
264     50      3.2
265     51      3.2
266     52      3.2
267     53      3.2
268     54      3.3
269     55      3.3
270     56      3.3
271     57      3.3
272     58      3.3
273     59      3.3
274     60      3.3
275     61      3.3
276     62      3.3
277     63      3.3
278     64      3.3
279     65      3.3
280     66      3.3
281     67      3.3
282     68      3.3
283     69      3.3
284     70      3.3
285     71      3.4
286     72      3.4
287     73      3.4
288     74      3.4
289     75      3.4
290     76      3.4
291     77      3.4
292     78      3.4
293     79      3.5
294     80      3.5
295     81      3.5
296     82      3.5
297     83      3.5
298     84      3.5
299     85      3.6
300     86      3.6
301     87      3.6
302     88      3.6
303     89      3.6
304     90      3.7
305     91      3.7
306     92      3.7
307     93      3.8
308     94      3.8
309     95      3.8
310     96      3.8
311     97      3.9
312     98      3.9
313     99      4
314     100     4
315     101     4
316     102     4.1
317     103     4.1
318     104     4.2
319     105     4.2
320     106     4.3
321     107     4.3
322     108     4.4
323     109     4.4
324     110     4.5
325     111     4.5
326     112     4.6
327     113     4.7
328     114     4.7
329     115     4.8
330     116     4.9
331     117     5
332     118     5
333     119     5.1
334     120     5.2
335     121     5.3
336     122     5.4
337     123     5.5
338     124     5.5
339     125     5.6
340     126     5.7
341     127     5.8
342     128     6
343     129     6.1
344     130     6.2
345     131     6.3
346     132     6.4
347     133     6.5
348     134     6.7
349     135     6.8
350     136     6.9
351     137     7
352     138     7.2
353     139     7.3
354     140     7.5
355     141     7.6
356     142     7.8
357     143     7.9
358     144     8.1
359     145     8.2
360     146     8.4
361     147     8.6
362     148     8.7
363     149     8.9
364     150     9.1
365     151     9.3
366     152     9.3



Michael Rennie
M.Sc. Candidate
University of Toronto at Mississauga
3359 Mississauga Rd. N.
Mississauga, ON  L5L 1C6
Ph: 905-828-5452  Fax: 905-828-3792
	[[alternative HTML version deleted]]


