N.I.M.R.O.D.  

optimization.f90

Go to the documentation of this file.
00001 !------------------------------------------------------------------------------
00002 !
00003 ! VERSION : 1.0
00004 !
00005 ! MODULE: optimization.f90 function optimization
00081 
00082 subroutine optim(b,m,ni,v,rl,istop,conv,llikelihood,LCV_aOUT)
00083     use WorkingSharedValues
00084     use mpimod
00085     use M_vpropre
00086     IMPLICIT NONE
00087     double precision,dimension(m),intent(inout) :: b
00088     double precision,intent(out) :: rl
00089     double precision,dimension(m*(m+3)/2),intent(out) :: v
00090     integer, intent(in) :: m
00091     integer, intent(out) :: ni,istop
00092     double precision, dimension(3), intent(out) :: conv
00093     double precision, intent(out) ::llikelihood
00094     double precision,dimension(2), intent(out) ::LCV_aOUT
00095 
00096     double precision,dimension(m*(m+1)/2+m) :: fu
00097     double precision da,dm,ga,trca,cb,epsa,epsb,tpscpu,funcpa, 
00098         det,step,eps,epsd,vw,fi,maxt,z,rl1,th,ep, 
00099         dd,rdm,ca,tr,cbsigne,vwtemp,rlBestLineSearch,LCV_a,LCV_aRVS,LCV_aMARC
00100     double precision,dimension(m) ::delta,VP,D,WI,b2,bh, 
00101         delta_gradient,deltaBestLineSearch
00102     double precision,dimension(m,m)::v2,matcorr,vecp,Vtemp
00103     integer :: nql,ii,nfmax,idpos,ier,ncount,id,j,i,ij,k,compteurMARC,jd,m1, 
00104         countsearpas,ierMARC,ierRVS
00105     integer, dimension(m) :: intc
00106     double precision,dimension(m*(m+3)/2)::vRVS,vMARC
00107 
00108     double precision, dimension(npmcomp) :: errorevaluation
00109     logical :: rvs
00110     double precision :: choldet,seuilSEARPAS
00111     integer ::Check_Matrix
00112 
00113 
00114     ! VARIABLES INITIALISATION
00115     id=0
00116     jd=0
00117     z=0.d0
00118 
00119     th=searpas_th
00120     eps=gonfle_eps
00121     epsa=controleCA
00122     epsb=controleCB
00123     epsd=controleRDM
00124     ep=dsinv_ep
00125     da=gonfle_da
00126     dm=gonfle_dm
00127 
00128     ca=epsa+1.d0
00129     cb=epsb+1.d0
00130 
00131     rl1=-1.d+10
00132     ni=0
00133     istop=0
00134 
00135     nql=1
00136     nfmax=m*(m+1)/2
00137 
00138 
00139     !########################################
00140     abserrfuncpa=5.d-1 !1.d-3 !1.d-3 !1 !1.d-3 5.d-1
00141     abserr(1:m)=1.d-1 !3 !1 !1.d-3*dd/(2*nbpatient*epsd)
00142     compteurMARC=0
00143     !########################################
00144 
00145 
00146     Main:Do
00147 
00148         ! RESTRICT THE OUTPUT FILE LENGTH TO NumberKeptOutput INTERATIONS
00149         if(mod(ni,NumberKeptOutput).EQ.0) then
00150             close(999)
00151             open(999,file="./"//trim(adjustl(OutputFolder))//"/"//trim(adjustl(nomfileOut)))
00152         end if
00153 
00154         if(numproc.EQ.0) then
00155             if (LOG_PRINT) then
00156                 write(0,*)
00157                 write(0,*)'******************************************'
00158                 write(0,*)'Iteration :', ni
00159                 write(0,*)'Parameters:'
00160                 write(0,*)"b(1:npm)=(/"
00161                 write(0,'(20(ES30.15,","))')(b(i),i=1,m)
00162                 write(0,*)"/)"
00163                 flush(0)
00164             end if
00165             write(999,*)
00166             write(999,*)'******************************************'
00167             write(999,*)'Iteration :', ni
00168             write(999,*)'Parameters:'
00169             write(999,*)"b(1:npm)=(/"
00170             write(999,'(20(ES30.15,","))')(b(i),i=1,m)
00171             write(999,*)"/)"
00172             flush(999)
00173         end if
00174 
00175         !########################################
00176         !test_funcpa=.false.
00177         !########################################
00178 
00179         !
00180         ! *** Log Likelihood, gradients and hessian computation
00181         !
00182         if(marqONLY) then
00183             if ((numproc.EQ.0).and.(LOG_PRINT)) write(0,*)"OPTIMIZATION TYPE : MARQUARDT"
00184             if(numproc.EQ.0) write(999,*)"OPTIMIZATION TYPE : MARQUARDT"
00185             call derivMARC(b,m,v,rl)
00186             rvs=.false.
00187         else
00188             if ((rdm.gt.epsd).and. &
00189                 (cb.LT.switchMarquardtCb).and. &
00190                 (ca.LT.switchMarquardtCa) .and. &
00191                 SwitchMarquartd .and. &
00192                 compteurMARC.LE.5) then
00193                 if ((numproc.EQ.0).and.(LOG_PRINT)) write(0,*)"OPTIMIZATION TYPE : MARQUARDT"
00194                 if(numproc.EQ.0) write(999,*)"OPTIMIZATION TYPE : MARQUARDT"
00195                 call derivMARC(b,m,v,rl)
00196                 rvs=.false.
00197                 compteurMARC=compteurMARC+1
00198             else
00199                 if ((numproc.EQ.0).and.(LOG_PRINT)) write(0,*)"OPTIMIZATION TYPE : ROBUST VARIANCE SCORING (RVS)"
00200                 if(numproc.EQ.0) write(999,*)"OPTIMIZATION TYPE : ROBUST VARIANCE SCORING (RVS)"
00201                 call derivRVS(b,m,v,rl)
00202                 rvs=.true.
00203                 compteurMARC=0
00204             end if
00205         end if
00206         if(rl.EQ.0.D0) rl=-1.D9
00207         cb=rl-rl1
00208         rl1=rl
00209         call LCVa(b,m,LCV_a,llikelihood)
00210 
00211 
00212 
00213         !########################################
00214         abserr(1:m)=1.d0/(nbpatient)
00215         abserrfuncpa=1.d-3
00216         !########################################
00217 
00218         !cbsigne=rl1-rl
00219         !cb=dabs(rl1-rl)
00220 
00221 
00222 
00223         do i=1,m
00224             do j=i,m
00225                 ij=(j-1)*j/2+i
00226                 Vtemp(i,j)=v(ij)
00227                 Vtemp(j,i)=v(ij)
00228             end do
00229         END DO
00230         call cpu_time (tpscpu)
00231         if(numproc.EQ.0) then
00232             if (LOG_PRINT) then
00233                 write(0,*)"----> LL, GRADIENT, HESSIAN SUMMARY"
00234                 write(0,*)"TIME :",tpscpu
00235                 write(0,*)'Target-function to maximise value =',rl1,' LCV_a =',LCV_a
00236                 write(0,'("Penalised Gradients =",200(ES30.15,",",1x))')(v(i),i=m*(m+1)/2+1,m*(m+3)/2)
00237                 if(rvs) then
00238                     write(0,*)"APPROXIMATION OF THE HESSIAN : G"
00239                 else
00240                     write(0,*)"HESSIAN : H"
00241                 end if
00242                 do i=1,m
00243                     write(0,'(200(ES11.3,1x))') Vtemp(i,1:i)
00244                 end do
00245                 flush(0)
00246             end if
00247             write(999,*)"----> LL, GRADIENT, HESSIAN SUMMARY"
00248             write(999,*)"TIME :",tpscpu
00249             write(999,*)'Target-function to maximise value =',rl1,' LCV_a =',LCV_a
00250             write(999,'("Penalised Gradients =",200(ES30.15,",",1x))')(v(i),i=m*(m+1)/2+1,m*(m+3)/2)
00251             if(rvs) then
00252                 write(999,*)"APPROXIMATION OF THE HESSIAN : G"
00253             else
00254                 write(999,*)"HESSIAN : H"
00255             end if
00256             do i=1,m
00257                 write(999,'(200(ES11.3,1x))') Vtemp(i,1:i)
00258             end do
00259             flush(999)
00260         end if
00261 
00262 
00263 
00264         !
00265         ! *** RDM Computation : Main Convergence criteria
00266         !
00267         rdm=0.D0
00268         fu=0.0D0
00269         do i=1,m
00270             do j=i,m
00271                 ij=(j-1)*j/2+i
00272                 fu(ij)=v(ij)
00273             end do
00274         END DO
00275 
00276         ier=0
00277         call dsinv(fu,m,ep,ier,det)
00278 
00279         if (ier.eq.-1) then
00280             rdm=epsd+1.d0
00281         else
00282             dd = 0.d0
00283             do i=1,m
00284                 do j=1,m
00285                     if(j.ge.i) then
00286                         ij=(j-1)*j/2+i
00287                     else
00288                         ij=(i-1)*i/2+j
00289                     end if
00290                     dd = dd + v(nfmax+i)*fu(ij)*V(nfmax+j)
00291                 end do
00292             end do
00293             rdm=dd/dble(m)
00294         end if
00295         do i=1,m
00296             do j=i,m
00297                 ij=(j-1)*j/2+i
00298                 Vtemp(i,j)=fu(ij) ! Mat Var-Cov.
00299                 Vtemp(j,i)=fu(ij)
00300             end do
00301         END DO
00302 
00303         if(numproc.EQ.0) then
00304             if (LOG_PRINT) then
00305                 write(0,*)"----> CONVERGENCE SUMMARY"
00306                 if(ier.EQ.0) then
00307                     write(0,*)'Inversion OK RDM=',rdm
00308                     write(0,*)"Variance-Covariance Matrix"
00309                     do i=1,m
00310                         write(0,'(200(ES11.3,1x))') Vtemp(i,1:i)
00311                     end do
00312                 else
00313                     write(0,*)'Inversion Failure of H/G : RDM=', rdm
00314                 end if
00315                 flush(0)
00316             end if
00317             if(ier.EQ.0) then
00318                 write(999,*)"----> CONVERGENCE SUMMARY"
00319                 write(999,*)'Inversion OK RDM=',rdm
00320                 write(999,*)"Variance-Covariance Matrix"
00321                 do i=1,m
00322                     write(999,'(200(ES11.3,1x))') Vtemp(i,1:i)
00323                 end do
00324             else
00325                 write(999,*)'Inversion Failure of H/G : RDM=', rdm
00326             end if
00327             flush(999)
00328         end if
00329 
00330         !########################################
00331         abserr(1:m)= 1.d-3*rdm*m/nbpatient
00332         abserrfuncpa= 5.d-2*rdm*m/nbpatient
00333         !########################################
00334 
00335         conv(1)=ca
00336         conv(2)=cb
00337         conv(3)=rdm
00338         if(numproc.EQ.0) then
00339             if (LOG_PRINT) then
00340                 write(0,*) 'ca =',ca,' cb =',cb, 'RDM=',rdm
00341                 flush(0)
00342             end if
00343             write(999,*) 'ca =',ca,' cb =',cb, 'RDM=',rdm
00344             flush(999)
00345             write(1234,*)ni,loglike(1),loglike(2),loglike(3),loglike(4),nbpatOK, &
00346                 conv(1),conv(2),conv(3),LCV_a,ier, &
00347                 rvs,b(1:m),(Vtemp(i,i),i=1,m)
00348             flush(1234)
00349         end if
00350 
00351         !
00352         ! *** Eigen vector and value /// Specific convergence assessment
00353         !
00354 
00355         if(HessianDiagnostic) then
00356             do i=1,m
00357                 do j=i,m
00358                     ij=(j-1)*j/2+i
00359                     Vtemp(i,j)=v(ij) ! Mat Var-Cov.
00360                     Vtemp(j,i)=v(ij)
00361                 end do
00362             END DO
00363             call eigen(m,m,Vtemp,VP,wi,vecp,D,INTC)
00364             if(numproc.EQ.0) then
00365                 if (LOG_PRINT) then
00366                     write(0,*)"----> HESSIAN DIAGNOSTIC SUMMARY"
00367                     write(0,*)"Eigen Value Test : Negative and extremely small values are a problem"
00368                     write(0,*)"                   See associated eigen vector to identify the parameter"
00369                     write(0,*)"                   with low identifiability."
00370                     write(0,'("H/G Eigen values :",200(ES30.15,1x))')(vp(i),i=1,m)
00371                     write(0,*)"Eigen Vector Test : High values are associated with parameters of low identifiability."
00372                     do i=1,m
00373                         write(0,'("H/G Eigen Vector",i4," : ",200(ES30.15,1x))')i,vecp(1:m,i)
00374                     end do
00375                     flush(0)
00376                 end if
00377                 write(999,*)"----> HESSIAN DIAGNOSTIC SUMMARY"
00378                 write(999,*)"Eigen Value Test : Negative and extremely small values are a problem"
00379                 write(999,*)"                   See associated eigen vector to identify the parameter"
00380                 write(999,*)"                   with low identifiability."
00381                 write(999,'("H/G Eigen values :",200(ES30.15,1x))')(vp(i),i=1,m)
00382                 write(999,*)"Eigen Vector Test : High values are associated with parameters of low identifiability."
00383                 do i=1,m
00384                     write(999,'("H/G Eigen Vector",i4," : ",200(ES30.15,1x))')i,vecp(1:m,i)
00385                 end do
00386                 flush(999)
00387             end if
00388         end if
00389 
00390         !
00391         ! *** Algorithm convergence assessment and Restarts assessment
00392         !
00393         if (ca .lt. epsa .and. dabs(cb) .lt. epsb .and. rdm.lt.epsd ) exit main
00394         if (ca .lt. seuilStuckca .and. dabs(cb) .lt. seuilStuckcb .and. rdm.gt.epsd .and. ni.ge.nbiterationMIN) goto 445
00395         !flush(0)
00396 
00397 
00398         !
00399         ! *** Diagonal term of the hessian are inflated if needed to find the movement direction. After a sufficient number
00400         ! *** of step the inflation is neglectable.
00401         !
00402         if(.not.(estimationMeasurmentError))call computeMeasurmentError(b,errorevaluation)
00403         tr=0.d0
00404         do i=1,m
00405             ii=i*(i+1)/2
00406             tr=tr+dabs(v(ii))
00407         end do
00408         tr=tr/dble(m)
00409 
00410         ncount=0
00411         ga=gonfle_ga
00412 
00413 400 continue
00414     do i=1,(m*(m+1)/2)+m
00415         fu(i)=v(i)
00416     end do
00417     !if(ncount.NE.0) then
00418     do i=1,m
00419         ii=i*(i+1)/2
00420         if (v(ii).ne.0) then
00421             fu(ii)=v(ii)+da*((1.d0-ga)*dabs(v(ii))+ga*tr)
00422         else
00423             fu(ii)=da*ga*tr
00424         endif
00425     end do
00426     !end if
00427     call dchole(fu,m,nql,idpos)
00428     if (idpos.ne.0) then
00429         ncount=ncount+1
00430         !if(ncount.GT.1) then
00431         if (ncount.le.3.or.ga.ge.1.d0) then
00432             da=da*dm
00433         else
00434             ga=ga*dm
00435             if (ga.gt.1.d0) ga=1.d0
00436         endif
00437         !end if
00438         if (ncount .LE. 100000) goto 400
00439     else
00440         vw=1.D0
00441         do i=1,m
00442             delta(i)=fu(nfmax+i)
00443             if(.not.(estimationMeasurmentError))delta(npmbio+npmexpl+ndim+1:npm) &
00444                                =errorevaluation(1:npmcomp)-b(npmbio+npmexpl+ndim+1:npm)
00445             b2(i)=b(i)+delta(i)
00446         end do
00447 
00448         !########################################
00449         !test_funcpa=.true.
00450         !########################################
00451         if(nbpatOK.EQ.0)firstFuncpa=.true.
00452         rl=funcpa(b2,m,id,z,id,z)
00453         if(nbpatOK.EQ.0)firstFuncpa=.false.
00454 
00455         if(numproc.EQ.0) then
00456             if (LOG_PRINT) then
00457                 if (rl1.ge.rl) then
00458                     write(0,*)"----> DISPLACEMENT SUMMARY : No amelioration"
00459                 else
00460                     write(0,*)"----> DISPLACEMENT SUMMARY : WITH amelioration"
00461                 end if
00462                 write(0,*)"Inflation parameters : ncount=",ncount,"da=",da,"dm=",dm,"ga=",ga
00463                 write(0,*)'Proposed Parameters:'
00464                 write(0,*)"bPROPOSED(1:npm)=(/"
00465                 write(0,'(20(ES30.15,","))')(b2(i),i=1,m)
00466                 write(0,*)"/)"
00467                 write(0,*)'New Target-function to maximise value = ',rl
00468                 flush(0)
00469             end if
00470             if (rl1.ge.rl) then
00471                 write(999,*)"----> DISPLACEMENT SUMMARY : No amelioration"
00472             else
00473                 write(999,*)"----> DISPLACEMENT SUMMARY : WITH amelioration"
00474             end if
00475             write(999,*)"Inflation parameters : ncount=",ncount,"da=",da,"dm=",dm,"ga=",ga
00476             write(999,*)'Proposed Parameters:'
00477             write(999,*)"bPROPOSED(1:npm)=(/"
00478             write(999,'(20(ES30.15,","))')(b2(i),i=1,m)
00479             write(999,*)"/)"
00480             write(999,*)'New Target-function to maximise value = ',rl
00481             flush(999)
00482         end if
00483         if (rl1.lt.rl) then
00484             if(da.lt.eps) then
00485                 da=eps
00486             else
00487                 da=da/(dm+2.d0)
00488             endif
00489             deltaBestLineSearch=delta
00490             goto 800
00491         else
00492             rlBestLineSearch=rl
00493             deltaBestLineSearch=delta
00494         endif
00495     endif
00496 
00497 
00498     !
00499     ! *** First test for new combinaison of parameters
00500     !
00501     !if ((numproc.EQ.0).and.(LOG_PRINT)) write(0,*)"b2",b2(1:m)
00502 
00503     !deltabest=delta
00504     !if ((numproc.EQ.0).and.(LOG_PRINT)) write(0,*) 'rl1 =',rl1,' rl =',rl
00505 
00506 
00507 
00508     !if (rl1.gt.rl) da=(dm-3.d0)*da
00509 
00510     !
00511     ! *** Second test : Line search performing to find a better combination of parameters
00512     !
00513     seuilSEARPAS=0.0D0
00514     if(rdm.GT.seuilSEARPASRDM) seuilSEARPAS=seuilSEARPASOK
00515 
00516     call dmaxt(maxt,delta,m)
00517     if(maxt.eq.0.D0) then
00518         vw=th
00519     else
00520         vw=DMAX1(th/maxt,searpas_th)
00521     endif
00522     vwtemp=vw
00523     step=dlog(1.5d0)
00524     call searpas(vw,step,b,bh,m,delta,fi,funcpa)
00525     rl=-fi
00526     if(numproc.EQ.0) then
00527         if (LOG_PRINT) then
00528             if ((rl1.lt.rl).and.(abs(rl1-rl).GT.seuilSEARPAS)) then
00529                 write(0,*)"----> FIRST LINE SEARCH SUMMARY : WITH amelioration"
00530             else
00531                 write(0,*)"----> FIRST LINE SEARCH SUMMARY : No amelioration"
00532             end if
00533             write(0,*)"Line Search parameters : vw.in=",vwtemp,"vw.out=",vw,"step=",step
00534             write(0,*)'Proposed Parameters:'
00535             write(0,*)"bPROPOSED(1:npm)=(/"
00536             write(0,'(20(ES30.15,","))')(b(i)+vw*delta(i),i=1,m)
00537             write(0,*)"/)"
00538             write(0,*)'New Target-function to maximise value = ',rl
00539             flush(0)
00540         end if
00541         if ((rl1.lt.rl).and.(abs(rl1-rl).GT.seuilSEARPAS)) then
00542             write(999,*)"----> FIRST LINE SEARCH SUMMARY : WITH amelioration"
00543         else
00544             write(999,*)"----> FIRST LINE SEARCH SUMMARY : No amelioration"
00545         end if
00546         write(999,*)"Line Search parameters : vw.in=",vwtemp,"vw.out=",vw,"step=",step
00547         write(999,*)'Proposed Parameters:'
00548         write(999,*)"bPROPOSED(1:npm)=(/"
00549         write(999,'(20(ES30.15,","))')(b(i)+vw*delta(i),i=1,m)
00550         write(999,*)"/)"
00551         write(999,*)'New Target-function to maximise value = ',rl
00552         flush(999)
00553     end if
00554     if ((rl1.lt.rl).and.(abs(rl1-rl).GT.seuilSEARPAS)) then
00555         deltaBestLineSearch=vw*delta
00556         da=(dm-3.d0)*da
00557         go to 800
00558     else
00559         if(rlBestLineSearch.lt.rl) then
00560             rlBestLineSearch=rl
00561             deltaBestLineSearch=vw*delta
00562         end if
00563     end if
00564 
00565     if(subsequentSearpas) then
00566         vw=subsequentsearpas_th
00567         vwtemp=vw
00568         countsearpas=0
00569 298 continue
00570     countsearpas=countsearpas+1
00571     vw=vwtemp*10.D0
00572     step=dlog(1.5d0)
00573     vwtemp=vw
00574     call searpas(vw,step,b,bh,m,delta,fi,funcpa)
00575     rl=-fi
00576     if(numproc.EQ.0) then
00577         if (LOG_PRINT) then
00578             if ((rl1.lt.rl).and.(abs(rl1-rl).GT.seuilSEARPAS)) then
00579                 write(0,*)"----> SUBSEQUENT",countsearpas," LINE SEARCH SUMMARY : WITH amelioration"
00580             else
00581                 write(0,*)"----> SUBSEQUENT",countsearpas," LINE SEARCH SUMMARY : No amelioration"
00582             end if
00583             write(0,*)"Line Search parameters : vw.in=",vwtemp,"vw.out=",vw,"step=",step
00584             write(0,*)'Proposed Parameters:'
00585             write(0,*)"bPROPOSED(1:npm)=(/"
00586             write(0,'(20(ES30.15,","))')(b(i)+vw*delta(i),i=1,m)
00587             write(0,*)"/)"
00588             write(0,*)'New Target-function to maximise value = ',rl
00589             flush(0)
00590         end if
00591         if ((rl1.lt.rl).and.(abs(rl1-rl).GT.seuilSEARPAS)) then
00592             write(999,*)"----> SUBSEQUENT",countsearpas," LINE SEARCH SUMMARY : WITH amelioration"
00593         else
00594             write(999,*)"----> SUBSEQUENT",countsearpas," LINE SEARCH SUMMARY : No amelioration"
00595         end if
00596         write(999,*)"Line Search parameters : vw.in=",vwtemp,"vw.out=",vw,"step=",step
00597         write(999,*)'Proposed Parameters:'
00598         write(999,*)"bPROPOSED(1:npm)=(/"
00599         write(999,'(20(ES30.15,","))')(b(i)+vw*delta(i),i=1,m)
00600         write(999,*)"/)"
00601         write(999,*)'New Target-function to maximise value = ',rl
00602         flush(999)
00603     end if
00604     if ((rl1.lt.rl).and.(abs(rl1-rl).GT.seuilSEARPAS)) then
00605         deltaBestLineSearch=vw*delta
00606         da=(dm-3.d0)*da
00607         go to 800
00608     else
00609         if(rlBestLineSearch.lt.rl) then
00610             rlBestLineSearch=rl
00611             deltaBestLineSearch=vw*delta
00612         end if
00613         if(countsearpas.LE.maxsubsequentSEARPAS)go to 298
00614     end if
00615 end if
00616 
00617 da=(dm-3.d0)*da
00618 
00619 800 continue
00620     ! Parameters displacement
00621     ca=0.d0
00622     do i=1,m
00623         ca=ca+deltaBestLineSearch(i)*deltaBestLineSearch(i)
00624     end do
00625     ! New proposal for new iteration
00626     do i=1,m
00627         b(i)=b(i)+deltaBestLineSearch(i)
00628     end do
00629 
00630     ni=ni+1
00631     if (ni .gt. maxiteration) then
00632         exit main
00633     end if
00634      !if(.NOT.(rvs))call derivRVS(b,m,v,rl)
00635 ! ier=0
00636 ! call dsinv(v,m,ep,ier,det)
00637 ! do i=1,m
00638 !      if((numproc.EQ.0).and.(LOG_PRINT))write(0,*)"parametre ",i,b(i),sqrt(v(i*(i+1)/2))
00639 !      if(numproc.EQ.0)write(999,*)"parametre ",i,b(i),sqrt(v(i*(i+1)/2))
00640 !     if(numproc.EQ.0)write(999,*)"parametreLOG ",i,exp(b(i)),exp(b(i))*sqrt(v(i*(i+1)/2))
00641 ! end do
00642 !  call LCVa(b,m,LCV_a,llikelihood)
00643 !  if((numproc.EQ.0).and.(LOG_PRINT))write(0,*)"LCVa ",LCV_a
00644 ! if(numproc.EQ.0)write(999,*)"LCVa ",LCV_a
00645 
00646 ! goto 110
00647 
00648 End do Main
00649 
00650 if(ni .LE. maxiteration) then
00651     istop=1
00652 else
00653     istop=2
00654 end if
00655 withexclusion=.false.
00656 
00657 call derivRVS(b,m,v,rl)
00658 ier=0
00659 call dsinv(v,m,ep,ierRVS,det)
00660 vRVS=v
00661 call LCVa(b,m,LCV_aRVS,llikelihood)
00662 
00663 
00664 call derivMARC(b,m,v,rl)
00665 ier=0
00666 call dsinv(v,m,ep,ierMARC,det)
00667 vMARC=v
00668 
00669 if ((ierRVS.NE.0).and.(ierMARC.EQ.0)) then
00670     istop=3
00671 end if
00672 call LCVa(b,m,LCV_aMARC,llikelihood)
00673 
00674 LCV_aOUT(1)=LCV_aRVS
00675 LCV_aOUT(2)=LCV_aMARC
00676 v=vRVS
00677 445 continue
00678     if(istop.EQ.0)istop=4 ! Algorithm got stuck it must be restarted.
00679 110 continue
00680     if(numproc.EQ.0) then
00681         if (LOG_PRINT) then
00682             write(0,*)"----> ***FINAL*** SUMMARY : "
00683             write(0,*)"Convergence Status :",istop
00684             if (istop.EQ.4)  then
00685                 write(0,*)'Estimation completely failed : Target-function to maximise value =',rl
00686             else
00687                 if (istop.EQ.3) then
00688                     write(0,*)'Hessian Non-Inversible'
00689                 else
00690                     if (istop.EQ.2) write(0,*)'Maximum number of iterations reached BUT at this point: '
00691                     write(0,*) 'Inversibility Hessian (RVS,MARC)=',1-ierRVS,1-ierMARC
00692                     write(0,*) 'LCVa (RVS,MARC)=',LCV_aRVS,LCV_aMARC
00693                     do i=1,m
00694                         write(0,*)"Parameter estimation (mean/sd.) :",i,b(i),sqrt(vRVS(i*(i+1)/2)),sqrt(vMARC(i*(i+1)/2))
00695                     end do
00696                 end if
00697             end if
00698             flush(0)
00699         end if
00700         write(999,*)"----> ***FINAL*** SUMMARY : "
00701         write(999,*)"Convergence Status :",istop
00702         if (istop.EQ.4)  then
00703             write(999,*)'Estimation completely failed : Target-function to maximise value =',rl
00704         else
00705             if (istop.EQ.3) then
00706                 write(999,*)'Hessian Non-Inversible'
00707             else
00708                 if (istop.EQ.2) write(999,*)'Maximum number of iterations reached BUT at this point: '
00709                 write(999,*) 'Inversibility Hessian (RVS,MARC)=',1-ierRVS,1-ierMARC
00710                 write(999,*) 'LCVa (RVS,MARC)=',LCV_aRVS,LCV_aMARC
00711                 do i=1,m
00712                     write(999,*)"Parameter estimation (mean/sd.) :",i,b(i),sqrt(vRVS(i*(i+1)/2)),sqrt(vMARC(i*(i+1)/2))
00713                 end do
00714             end if
00715         end if
00716         flush(999)
00717     end if
00718 
00719     return
00720 end
00721 
00722 
00723 !------------------------------------------------------------------------------
00724 !
00725 ! VERSION : 1.0
00726 !
00727 ! MODULE: optimization.f90 function LCV_a
00757 subroutine LCVa(b,m,LCV_a,llikelihood)
00758     use WorkingSharedValues
00759     use mpimod
00760     implicit none
00761     double precision,dimension(m),intent(in) :: b
00762     integer,intent(in) :: m
00763     double precision,intent(inout) :: llikelihood
00764     double precision :: det,ep,trace
00765     double precision, intent(out)::LCV_a
00766     integer :: i,j,ij,ier
00767     double precision, dimension(m,m)::vpinv, vnonp,hh
00768     double precision,dimension(m*(m+1)/2+m) :: fu
00769     !
00770     ! *** Penalized and unpenalized hessian computation
00771         do i=1,m
00772             do j=i,m
00773                 ij=(j-1)*j/2+i
00774                 vnonp(i,j)=hessianNonPenalisee(ij)
00775                 vnonp(j,i)=hessianNonPenalisee(ij)
00776             end do
00777         END DO
00778     !
00779     ! *** Inversions
00780     !
00781     do i=1,npm
00782         do j=i,npm
00783             ij=(j-1)*j/2+i
00784             fu(ij)=hessianPenalisee(ij)
00785         end do
00786     end do
00787     ier=0
00788     call dsinv(fu,npm,ep,ier,det)
00789 
00790         do i=1,m
00791             do j=i,m
00792                 ij=(j-1)*j/2+i
00793                 vpinv(i,j)=fu(ij)
00794                 vpinv(j,i)=fu(ij)
00795             end do
00796         END DO
00797     !
00798     ! *** Computation of the product Hp-1H
00799     !
00800     hh=matmul(vpinv,vnonp)
00801 
00802     !
00803     ! *** Computation of the trace of Hp-1H
00804     !
00805     do i=1,npm
00806           trace=trace+(hh(i,i))
00807     end do
00808     !
00809     ! *** LCV_a formula
00810     !
00811     if(withexclusion) then
00812         LCV_a=-(LLnonPenalisee-trace)/nbpatOK
00813     else
00814         LCV_a=-(LLnonPenalisee-trace)/nbpatient
00815     end if
00816     return
00817 end subroutine LCVa
00818 
00819 
00820 subroutine computeMeasurmentError(b,err)
00821     use WorkingSharedValues
00822     implicit none
00823     integer :: i,k,j
00824     double precision,dimension(npm),intent(in)::b
00825 
00826     double precision,dimension(npmcomp),intent(out):: err
00827     double precision, dimension(ndim)::X
00828     double precision,dimension(tdef2,npmcomp):: deriv
00829 
00830     do i=1,nbpatient
00831         emp_bayes(:,i)=b(:)
00832     end do
00833     ! PEB calculation for parameters with random effects
00834     do i=1,ndim
00835         emp_bayes(i,1:nbpatient)=b(i)+startsauv(i,1:nbpatient)*b(npmbio+npmexpl+i)
00836     end do
00837 
00838     do i=1,nbpatient
00839         numpat1=i
00840         systeme=1
00841         X=0.d0
00842         b1=0.0D0
00843         b1(1:npmbio+npmexpl)=emp_bayes(1:npmbio+npmexpl,i)
00844         call solution(x,bayes_data(i,:,1:npmcomp),deriv)
00845     end do
00846 
00847 
00848     err=0.0D0
00849     do k=1, npmcomp
00850         do i=1, nbpatient
00851             do j=1,nbobs(i,k)
00852                 err(k)=err(k)+(donnees(i,j,k)-bayes_data(i,schedule(numpat1,j,k)+1,k))**2
00853             end do
00854         end do
00855         err(k)=sqrt(err(k)/sum(nbobs(:,k)))
00856     end do
00857     return
00858 end subroutine computeMeasurmentError
00859