398 lines
13 KiB
FortranFixed
398 lines
13 KiB
FortranFixed
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
|
|
subroutine fluorescencejmax()
|
|
implicit none
|
|
include '../testarea/LeafGasParams.h'
|
|
include '../testarea/LeafGasHybridFit.h'
|
|
integer i,ndim,k,j,iderivative,iwrong,n,icompete,i2,isitnaninf,
|
|
&nave
|
|
double precision beta(4),sumsquare0,beta0(4),sumsquarecp,
|
|
&betacp(4),ftol,xtol,shortx(maxobs,2),shorty(maxobs),
|
|
&xvar(maxobs,2),weitx(maxobs,2),weity(maxobs),
|
|
&templflights0(maxobs),aparlights0(maxobs),termmin,termmax,
|
|
&ftol_relax,term1,term2,ran2,discount,history(2000,10),upper,lower,
|
|
&f1dim_flujmax,flujmax_pikaia
|
|
parameter(ftol=1.0d-7,xtol=1.0d-7)
|
|
external funkmin_flujmax,f1dim_flujmax,FCN_flujmax,flujmax_pikaia
|
|
!beta(1)=fjmax25
|
|
beta(1)=univparams(8)
|
|
betamin(1)=univparamsmin(8)
|
|
betamax(1)=univparamsmax(8)
|
|
!beta(2)=phifactor
|
|
beta(2)=univparams(11)
|
|
betamin(2)=univparamsmin(11)
|
|
betamax(2)=univparamsmax(11)
|
|
!beta(3)=thetafactor
|
|
beta(3)=univparams(12)
|
|
betamin(3)=univparamsmin(12)
|
|
betamax(3)=univparamsmax(12)
|
|
ndim=3
|
|
ntotlights=0
|
|
termmax=-1.0d+9
|
|
termmin=1.0d+9
|
|
do i=1,numALightcurves
|
|
do j=1,nALightPoints(i)
|
|
if(ALightchlflphips2(j,i).gt.0.0d0.and.
|
|
&j.le.nstartalight(i))then
|
|
!Only points before nstartalight are used because these points are apparently limited by RuBP regeneration and therefore
|
|
!the electron transport equation applies.
|
|
ntotlights=ntotlights+1
|
|
templflights(ntotlights)=ALighttempleaf(j,i)
|
|
if(templflights(ntotlights).lt.termmin)
|
|
&termmin=templflights(ntotlights)
|
|
if(templflights(ntotlights).gt.termmax)
|
|
&termmax=templflights(ntotlights)
|
|
aparlights(ntotlights)=ALightaPPFDlf(j,i)
|
|
flphips2lights(ntotlights)=ALightchlflphips2(j,i)
|
|
xvar(ntotlights,1)=aparlights(ntotlights)
|
|
xvar(ntotlights,2)=templflights(ntotlights)
|
|
weitx(ntotlights,1)=1.0d0
|
|
weitx(ntotlights,2)=1.0d0
|
|
weity(ntotlights)=1.0d0
|
|
templflights0(ntotlights)=templflights(ntotlights)
|
|
aparlights0(ntotlights)=aparlights(ntotlights)
|
|
endif
|
|
enddo
|
|
enddo
|
|
if((termmax-termmin).gt.2.0d0)then
|
|
ndim=4
|
|
!beta(4)=ha_jmax
|
|
beta(4)=univparams(17)
|
|
betamin(4)=univparamsmin(17)
|
|
betamax(4)=univparamsmax(17)
|
|
endif
|
|
if(ntotlights.lt.ndim)then
|
|
ntotlights=0
|
|
return
|
|
endif
|
|
isitbounded=1
|
|
call funkmin_flujmax(ndim,beta,flujmaxfval)
|
|
do i=1,ndim
|
|
beta0(i)=beta(i)
|
|
history(1,i)=beta(i)
|
|
enddo
|
|
sumsquare0=flujmaxfval
|
|
history(1,ndim+1)=flujmaxfval
|
|
!entrance counter
|
|
history(1,ndim+2)=1.0d0
|
|
!failure counter
|
|
history(1,ndim+3)=0.0d0
|
|
!Is it a competition among different initial guesses?
|
|
icompete=0
|
|
!j the total number of calls to nongradopt; k is the number of returns to the current best and reset
|
|
!to zero if a better minumum is found; n is the number of scouting points over the landscape of the cost function.
|
|
!The first initial guess provided by the user is always part of the set of scouting points.the rest consist of outcomes
|
|
!from calls to nongradopt if they are significantly different from the current best.
|
|
j=0
|
|
k=0
|
|
n=1
|
|
nave=1
|
|
ftol_relax=ftol*1000.0d0
|
|
discount=2.0d0
|
|
30 do i=1,ndim
|
|
betacp(i)=beta(i)
|
|
enddo
|
|
sumsquarecp=flujmaxfval
|
|
iderivative=0
|
|
iwrong=0
|
|
call odr_leastsquare(ndim,FCN_flujmax,beta,ntotlights,
|
|
&xvar(1:ntotlights,1:2),2,flphips2lights,1,weitx(1:ntotlights,1:2),
|
|
&weity,iderivative,shortx(1:ntotlights,1:2),shorty(1:ntotlights),
|
|
&flujmaxfval,iwrong)
|
|
!after odr_leastsquare, forcing variables are destroyed. restore to the origninals
|
|
do i=1,ntotlights
|
|
templflights(i)=templflights0(i)
|
|
aparlights(i)=aparlights0(i)
|
|
enddo
|
|
call funkmin_flujmax(ndim,beta,flujmaxfval)
|
|
if(isitnaninf(flujmaxfval).eq.1.or.flujmaxfval.gt.sumsquarecp)then
|
|
do i=1,ndim
|
|
beta(i)=betacp(i)
|
|
enddo
|
|
flujmaxfval=sumsquarecp
|
|
else
|
|
do i=1,ndim
|
|
betacp(i)=beta(i)
|
|
enddo
|
|
sumsquarecp=flujmaxfval
|
|
endif
|
|
call nongradopt(ndim,funkmin_flujmax,f1dim_flujmax,
|
|
&beta,betamin,betamax,ftol_relax,flujmaxfval)
|
|
if(isitnaninf(flujmaxfval).eq.1.or.flujmaxfval.gt.sumsquarecp)then
|
|
do i=1,ndim
|
|
beta(i)=betacp(i)
|
|
enddo
|
|
flujmaxfval=sumsquarecp
|
|
endif
|
|
if(flujmaxfval.gt.1.0d0)then
|
|
term1=flujmaxfval*ftol_relax
|
|
else
|
|
term1=ftol_relax*10.0d0
|
|
endif
|
|
if(flujmaxfval.gt.sumsquare0)then
|
|
!failure
|
|
if((flujmaxfval-sumsquare0).gt.term1)then
|
|
if(icompete.eq.1)history(1,ndim+3)=history(1,ndim+3)+1.5d0
|
|
!even though flujmaxfval is much worse than sumsquare0, it is an output of optimization after all so
|
|
!include it in the set if it has not already been included in the set.
|
|
i=1
|
|
i2=1
|
|
40 if(dabs(history(i2,i)-beta(i)).gt.ftol_relax)then
|
|
if(dabs(history(i2,ndim+1)-flujmaxfval).lt.term1)then
|
|
history(i2,ndim+3)=history(i2,ndim+3)+1.0d0
|
|
goto 60
|
|
endif
|
|
if(i2.ge.n)goto 50
|
|
i2=i2+1
|
|
i=1
|
|
goto 40
|
|
else
|
|
if(i.ge.ndim)goto 60
|
|
i=i+1
|
|
goto 40
|
|
endif
|
|
50 n=n+1
|
|
do i=1,ndim
|
|
history(n,i)=beta(i)
|
|
enddo
|
|
history(n,ndim+1)=flujmaxfval
|
|
history(n,ndim+2)=0.0d0
|
|
history(n,ndim+3)=0.0d0
|
|
!use average only when there is improvement
|
|
nave=n
|
|
else
|
|
!the difference is minimal even though flujmaxfval is larger than sumsquare0.
|
|
!Increment the counter for arriving at the same minimum.
|
|
if(icompete.eq.1)history(1,ndim+3)=history(1,ndim+3)+1.0d0
|
|
k=k+1
|
|
endif
|
|
60 do i=1,ndim
|
|
beta(i)=beta0(i)
|
|
enddo
|
|
flujmaxfval=sumsquare0
|
|
else
|
|
!success
|
|
if((sumsquare0-flujmaxfval).lt.term1)then
|
|
!negligible improvement. Increment the counter for arriving at the same minimum.
|
|
!no increment for the set of central initial guesses
|
|
if(icompete.eq.1)history(1,ndim+3)=history(1,ndim+3)+0.5d0
|
|
k=k+1
|
|
else
|
|
!reset the counter for arriving at a better minimum.
|
|
!Increment the set of central initial guesses
|
|
if(dabs(discount-2.0d0).lt.ftol)then
|
|
discount=dmax1(0.001d0,(sumsquare0-flujmaxfval)/1000.0d0)
|
|
endif
|
|
k=0
|
|
n=n+1
|
|
do i=1,ndim
|
|
history(n,i)=beta(i)
|
|
enddo
|
|
history(n,ndim+1)=flujmaxfval
|
|
history(n,ndim+2)=0.0d0
|
|
history(n,ndim+3)=0.0d0
|
|
endif
|
|
do i=1,ndim
|
|
beta0(i)=beta(i)
|
|
enddo
|
|
sumsquare0=flujmaxfval
|
|
endif
|
|
j=j+1
|
|
if(j.lt.200.and.k.lt.3)then
|
|
!first explore around the very first initial guess
|
|
if(j.lt.10)then
|
|
term1=0.05d0+dmin1(history(1,ndim+3)*0.1d0,0.9d0)
|
|
history(1,ndim+2)=history(1,ndim+2)+1.0d0
|
|
do i=1,ndim
|
|
lower=history(1,i)-term1*(history(1,i)-betamin(i))
|
|
upper=history(1,i)+term1*(betamax(i)-history(1,i))
|
|
beta(i)=lower+ran2()*(upper-lower)
|
|
enddo
|
|
icompete=1
|
|
goto 70
|
|
endif
|
|
!try average
|
|
if(n.gt.nave)then
|
|
term1=1.0d0/(history(1,ndim+1)+1.0d-5)
|
|
do i=2,n
|
|
term1=term1+1.0d0/(history(i,ndim+1)+1.0d-5)
|
|
enddo
|
|
do i=1,ndim
|
|
beta(i)=history(1,i)/(term1*(history(1,ndim+1)+1.0d-5))
|
|
do icompete=2,n
|
|
beta(i)=beta(i)+history(icompete,i)/
|
|
&(term1*(history(icompete,ndim+1)+1.0d-5))
|
|
enddo
|
|
enddo
|
|
nave=n
|
|
icompete=0
|
|
goto 70
|
|
endif
|
|
!try different initial guesses
|
|
if(ran2().gt.0.2d0)then
|
|
!guess around the best
|
|
icompete=1
|
|
term1=history(1,ndim+1)+
|
|
&discount*history(1,ndim+2)*history(1,ndim+3)
|
|
do i=2,n
|
|
term2=history(i,ndim+1)+
|
|
&discount*history(i,ndim+2)*history(i,ndim+3)
|
|
if(term2.le.term1)then
|
|
term1=term2
|
|
do i2=1,ndim+3
|
|
history(n+1,i2)=history(i,i2)
|
|
history(i,i2)=history(1,i2)
|
|
history(1,i2)=history(n+1,i2)
|
|
enddo
|
|
endif
|
|
enddo
|
|
term1=0.05d0+dmin1(history(1,ndim+2)*history(1,ndim+3)*
|
|
&0.015d0,0.9d0)
|
|
history(1,ndim+2)=history(1,ndim+2)+1.0d0
|
|
do i=1,ndim
|
|
lower=history(1,i)-term1*(history(1,i)-betamin(i))
|
|
upper=history(1,i)+term1*(betamax(i)-history(1,i))
|
|
beta(i)=lower+ran2()*(upper-lower)
|
|
enddo
|
|
else
|
|
!completely random guess
|
|
do i=1,ndim
|
|
beta(i)=betamin(i)+ran2()*(betamax(i)-betamin(i))
|
|
enddo
|
|
icompete=0
|
|
endif
|
|
70 call funkmin_flujmax(ndim,beta,flujmaxfval)
|
|
goto 30
|
|
else
|
|
if((ftol_relax-ftol).gt.ftol)then
|
|
if(k.le.1)then
|
|
n=n+1
|
|
do i=1,ndim+3
|
|
history(n,i)=history(1,i)
|
|
enddo
|
|
do i=1,ndim
|
|
history(1,i)=beta(i)
|
|
enddo
|
|
history(1,ndim+1)=flujmaxfval
|
|
history(1,ndim+2)=0.0d0
|
|
history(1,ndim+3)=0.0d0
|
|
do i=1,n
|
|
do icompete=1,ndim
|
|
betacp(icompete)=history(i,icompete)
|
|
enddo
|
|
sumsquarecp=history(i,ndim+1)
|
|
call RepeatCompassSearch(ndim,betacp,sumsquarecp,
|
|
&betamin,betamax,funkmin_flujmax,f1dim_flujmax,ftol_relax)
|
|
call funkmin_flujmax(ndim,betacp,sumsquarecp)
|
|
if(isitnaninf(sumsquarecp).eq.0.and.sumsquarecp.lt.
|
|
&flujmaxfval)then
|
|
do icompete=1,ndim
|
|
beta(icompete)=betacp(icompete)
|
|
enddo
|
|
flujmaxfval=sumsquarecp
|
|
endif
|
|
enddo
|
|
do i=1,ndim
|
|
beta0(i)=beta(i)
|
|
enddo
|
|
sumsquare0=flujmaxfval
|
|
j=0
|
|
icompete=1
|
|
else
|
|
icompete=0
|
|
endif
|
|
ftol_relax=ftol
|
|
goto 30
|
|
endif
|
|
endif
|
|
|
|
goto 110
|
|
|
|
call RepeatCompassSearch(ndim,beta,flujmaxfval,betamin,
|
|
&betamax,funkmin_flujmax,f1dim_flujmax,xtol)
|
|
call funkmin_flujmax(ndim,beta,flujmaxfval)
|
|
if(isitnaninf(flujmaxfval).eq.1.or.flujmaxfval.gt.sumsquare0)then
|
|
do i=1,ndim
|
|
beta(i)=beta0(i)
|
|
enddo
|
|
flujmaxfval=sumsquare0
|
|
else
|
|
do i=1,ndim
|
|
beta0(i)=beta(i)
|
|
enddo
|
|
sumsquare0=flujmaxfval
|
|
endif
|
|
do i=1,ndim
|
|
betacp(i)=(beta(i)-betamin(i))/(betamax(i)-betamin(i))
|
|
enddo
|
|
call pikaia(flujmax_pikaia,ndim,gacontrol,betacp,flujmaxfval,i)
|
|
if(i.eq.0)then
|
|
do i=1,ndim
|
|
betacp(i)=betamin(i)+betacp(i)*(betamax(i)-betamin(i))
|
|
enddo
|
|
call funkmin_flujmax(ndim,betacp,flujmaxfval)
|
|
if(isitnaninf(flujmaxfval).eq.0.and.flujmaxfval.lt.sumsquare0)
|
|
&then
|
|
do i=1,ndim
|
|
beta(i)=betacp(i)
|
|
beta0(i)=betacp(i)
|
|
enddo
|
|
sumsquare0=flujmaxfval
|
|
endif
|
|
endif
|
|
flujmaxfval=sumsquare0
|
|
iderivative=0
|
|
iwrong=0
|
|
call odr_leastsquare(ndim,FCN_flujmax,beta,ntotlights,
|
|
&xvar(1:ntotlights,1:2),2,flphips2lights,1,weitx(1:ntotlights,1:2),
|
|
&weity,iderivative,shortx(1:ntotlights,1:2),shorty(1:ntotlights),
|
|
&flujmaxfval,iwrong)
|
|
!after odr_leastsquare, forcing variables are destroyed. restore to the origninals
|
|
do i=1,ntotlights
|
|
templflights(i)=templflights0(i)
|
|
aparlights(i)=aparlights0(i)
|
|
enddo
|
|
call funkmin_flujmax(ndim,beta,flujmaxfval)
|
|
if(isitnaninf(flujmaxfval).eq.1.or.flujmaxfval.gt.sumsquare0)then
|
|
do i=1,ndim
|
|
beta(i)=beta0(i)
|
|
enddo
|
|
flujmaxfval=sumsquare0
|
|
endif
|
|
j=0
|
|
100 sumsquare0=flujmaxfval
|
|
do i=1,ndim
|
|
beta0(i)=beta(i)
|
|
enddo
|
|
call nongradopt(ndim,funkmin_flujmax,f1dim_flujmax,
|
|
&beta,betamin,betamax,ftol,flujmaxfval)
|
|
call funkmin_flujmax(ndim,beta,flujmaxfval)
|
|
if(flujmaxfval.eq.sumsquare0)return
|
|
if(isitnaninf(flujmaxfval).eq.1.or.flujmaxfval.gt.sumsquare0)then
|
|
do i=1,ndim
|
|
beta(i)=beta0(i)
|
|
enddo
|
|
flujmaxfval=sumsquare0
|
|
endif
|
|
sumsquarecp=flujmaxfval
|
|
do i=1,ndim
|
|
betacp(i)=beta(i)
|
|
enddo
|
|
call RepeatCompassSearch(ndim,betacp,sumsquarecp,betamin,
|
|
&betamax,funkmin_flujmax,f1dim_flujmax,xtol)
|
|
call funkmin_flujmax(ndim,betacp,sumsquarecp)
|
|
if(flujmaxfval.eq.sumsquarecp)return
|
|
if(isitnaninf(sumsquarecp).eq.0.and.flujmaxfval.gt.sumsquarecp)
|
|
&then
|
|
do i=1,ndim
|
|
beta(i)=betacp(i)
|
|
enddo
|
|
flujmaxfval=sumsquarecp
|
|
endif
|
|
j=j+1
|
|
if(j.le.2.and.(sumsquare0-flujmaxfval).gt.ftol)goto 100
|
|
!
|
|
!------------------------------------------------------
|
|
110 call funkmin_flujmax(ndim,beta,flujmaxfval)
|
|
return
|
|
END subroutine fluorescencejmax
|