# This is a model of how wood density affects growth rate in tree species, and starts with the table
# of survival and growth rates BCIdata.rdata (described in WoodDensMortBayesWorkshop.r).

# True growth rate is defined as the lg[j], the log of diameter increment of species j. Assume lg,
# (which is log-transformed) is across species, with mean a+b*w and SD sigma. To account for sample size,
# let sigma[j]=tau/sqrt(N[j]), where N is the number of individuals of species j. That is, tau is a constant
# (to be fitted), while sigma varies depending on sample size.

# This is a short-cut for handling species abundance. A complete model would include the growth of every individual.


# AS OF 20 MARCH 2007: THESE FUNCTIONS ARE NOT FINALIZED. THEY MAY RUN, BUT ARE NOT CORRECT.
growth.density=function(data,steps,showstep=25)
{
 data=subset(data,!is.na(lg)&growth>0&!is.na(wsg))
 nospp=dim(data)[1]
 growthparam=matrix(nrow=steps,ncol=4)
 spmean=matrix(nrow=steps,ncol=nospp)

 growthparam[1,1]=mean(data$lg)
 growthparam[1,2]=0
 growthparam[1,3]=1
 growthparam[1,4]=1
 stepsize=rep(1,4)

 spmean[1,]=rep(0,nospp)
 spstep=rep(1,nospp)


 llikeSD=rep(1,steps)
# llike[1]=,log=TRUE))

 for(i in 2:steps)
  {
   for(j in 1:3)
    {
     param=growthparam[i-1,]
     testparam=param[j]
     if(j>1) param[1:(j-1)]=growthparam[i,1:(j-1)]

     metropResult=metrop1step(func=llike.growthparam,start.param=testparam,param=param,whichtest=j,scale.param=stepsize[j],spmean=spmean[i-1,],
                              w=data$wsg,adjust=1.01,target=0.25)
     growthparam[i,j]=metropResult[1]
     stepsize[j]=metropResult[2]
    }

   param=growthparam[i,]
   metropResult=metrop1step(func=llike.growthsd,start.param=growthparam[i-1,4],param=param,scale.param=stepsize[4],spmean=spmean[i-1,],w=data$wsg,
                            N=data$Ngrow,adjust=1.01,target=0.25)
   growthparam[i,4]=metropResult[1]
   stepsize[4]=metropResult[2]

   for(j in 1:nospp)
    {
     metropResult=metrop1step(func=llike.growthspmean,start.param=spmean[i-1,j],param=param,scale.param=spstep[j],w=data$wsg,N=data$Ngrow,
                              adjust=1.01,target=0.25)
     growthparam[i,4]=metropResult[1]
     stepsize[4]=metropResult[2]
    }

   if(i%%showstep==0)
     cat(i, ": a ", round(param[i,1],3), ": b ", round(param[i,2],3), ": sigmaBetween ", round(param[i,3],3), ": sigmaWithin ", round(param[i,4],3),"\n")
  }

}


llike.growthparam=function(testparam,param,whichtest,spmean,w)
{
 param[whichtest]=testparam
 a=param[1]
 b=param[2]
 sigma=param[3]

 if(sigma<=0) return(-Inf)

 meangrowth=a+b*w
 llike=dnorm(spmean,mean=meangrowth,sd=sigma,log=TRUE)

 return(sum(llike))
}



llike.growthspmean=function(testparam,param)
{
 spmean=testparam
 a=param[1]
 b=param[2]
 sigma=param[3]
 
 llike1=

}