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

funcpa.f90

Go to the documentation of this file.
00001 !------------------------------------------------------------------------------
00002 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00003 !                effects based on Ordinary Differential equations
00004 !------------------------------------------------------------------------------
00005 !
00006 ! VERSION : 1.0
00007 !
00008 ! MODULE: funcpa.f90 function funcpa
00042 
00043 function funcpa(b,npm2,id,thi,jd,thj)
00044     use WorkingSharedValues
00045     use mpimod
00046     implicit none
00047 
00048     double precision::funcpa
00049     double precision, dimension(npm2),intent(inout)::b
00050     integer,intent(in) ::npm2,id,jd
00051     double precision,intent(in)::thi,thj
00052 
00053     integer::ntest
00054 
00055     double precision :: tpscpu,det,ep,deter,rnrnor 
00056         ,rl,result,rltest,seuil,rlobs,rlobs1, 
00057         rlobs2,startobs,pena,rdmMarcq
00058     double precision,dimension(ndim*(ndim+3)/2) :: v,v2
00059     double precision,dimension(ndim)::starttest,rlabs2,startobs2,optimParam
00060     double precision,dimension(ndim,ndim)::scale,scaleinv,scaleinv2
00061     integer :: ier,i,j,k,ij,nql,idpos,istop,compteur,NDIM2,maxit,nimax, 
00062         debut,fin,errorflag
00063 
00065     integer :: typeOptyimization
00066 
00067 
00068     numparam=1
00069 
00070     ntest=funcpa_ntest
00071     nimax=funcpa_nimax
00072     seuil=funcpa_marcqThrehold
00073 
00074     nql=1
00075     b1=0.d0
00076     ier=0
00077     det=0
00078     ep=1.d-10
00079 
00080     b1=b
00081     if (id.ne.0) b1(id)=b1(id)+thi
00082     if (jd.ne.0) b1(jd)=b1(jd)+thj
00083 
00084     systeme=0
00085     ndim2=ndim
00086     funcpa=0.d0
00087     likelihoodERROR(:,:)=0
00088     likelihoodPRECISION(:,:)=0.0D0
00089 
00090     if(firstFuncpa) then
00091         patOK(:)=1
00092         listIdPatExcluded(:)=0
00093     end if
00094 
00095     call repartirSurCoeurs(debut,fin,nbpatient)
00096     do i=1,debut-1
00097         startsauv(:,i) = 0.D0
00098         extrema(:,i) = 0.D0
00099         likelihoodERROR(:,i)=0
00100         likelihoodPRECISION(:,i)=0.0D0
00101         startsauvind(:,i) = 0.D0
00102         scaleinvsauv(:,:,i) = 0.D0
00103         scalesauv(:,:,i) = 0.D0
00104         scaleinv2sauv(:,:,i) = 0.D0
00105         detersauv(i) = 0.D0
00106         vrais_obs(:)=0.0D0
00107     end do
00108 
00109     do i=fin+1,nbpatient
00110         startsauv(:,i) = 0.D0
00111         scalesauv(:,:,i) = 0.D0
00112         extrema(:,i) = 0.D0
00113         startsauvind(:,i) = 0.D0
00114         scaleinvsauv(:,:,i) = 0.D0
00115         scaleinv2sauv(:,:,i) = 0.D0
00116         detersauv(i) = 0.D0
00117         vrais_obs(:)=0.0D0
00118         likelihoodERROR(:,i)=0
00119         likelihoodPRECISION(:,i)=0.0D0
00120     end do
00121 
00122     if ((LOG_PRINT).and.(numproc.EQ.0).and.(writefuncpaFichier))write(0,*) &
00123         "  NUM  IDPAT       Log-Likelihood                   OPT   ODE   AGQ",&
00124         "               OPT               ODE               AGQ                    NBOBS "
00125     do i=debut,fin
00126 
00127         numpat1=i
00128         maxit=100
00129         adaptive=.false.
00130         adaptive2=.false.
00131         rl=0.d0
00132         compteur=0
00133         rltest=0.d0
00134         rlabs2=0.d0
00135         rlobs=0.d0
00136         rlobs1=0.d0
00137         j=0
00138 
00139         !
00140         ! *** First optimization of random effect
00141         !
00142         typeOptyimization=1
00143         rdmMarcq=100.0D0
00144         optimParam=0.0D0
00145         !optimParam=startsauvind(:,i)
00146         call marquardt(optimParam(:),ndim,maxit,v,rltest,istop,nimax,rdmMarcq)
00147         if((istop.NE.1).and.(seuil2 .le. seuil))rdmMarcq=-99
00148         if (rltest.ne.rltest) rltest=-1.d9
00149         ! *** Best value storage
00150         rlobs2=rltest
00151         startobs2(:)=optimParam(:)
00152         !
00153         ! *** Second optimization of random effect
00154         ! *** si log_likelihood_i/nbobs_i>5
00155         !
00156         if ((istop.NE.1) .and.(seuil2 .gt. seuil)) then
00157             !
00158             ! *** Unidimensional optimization
00159             !
00160             typeOptyimization=2
00161             if ((LOG_PRINT).and.(writefuncpaFichier)) write(0,*)i,"Unidimensional optimization"
00162             do k=1,ndim
00163                 rlobs=-1.d10
00164                 rlobs1=-1.d10
00165                 do j=1,ntest
00166                     rltest=0.d0
00167                     starttest(k)=0.d0
00168                     numeroalea=k
00169                     starttest(k)=2*rnrnor()
00170                     ndim2=1
00171                     if(j.eq.1) starttest(k)=startsauvind(k,i)
00172                     rdmMarcq=100.0D0
00173                     call marquardt(starttest(k),ndim2,maxit,v,rltest,istop,nimax,rdmMarcq)
00174                     if (LOG_PRINT.AND. debugFUNCPA) write(0,*)"j",k,j,seuil2
00175                     if (rltest.ne.rltest) rltest=-1.d10
00176                     if (rltest .ge. rlobs) then
00177                         rlobs1=rlobs
00178                         rlobs=rltest
00179                         startobs=starttest(k)
00180                     end if
00181                     if (seuil2 .le. seuil) exit
00182                 end do
00183                 rlabs2(k)=rlobs
00184                 extrema(k,i)=startobs
00185             end do
00186             rlabs2(1:ndim)=1/rlabs2(1:ndim)
00187             starttest(:)=extrema(:,i)
00188             numeroalea=0
00189             !
00190             ! *** Unidirectionnel summary : Random effects Optimization with unidimensionnal starting points
00191             !
00192             startsauvind(:,i)=starttest(:)
00193             rdmMarcq=100.0D0
00194             call marquardt(starttest(:),ndim,maxit,v,rltest,istop,nimax,rdmMarcq)
00195             if((istop.NE.1).and.(seuil2 .le. seuil))rdmMarcq=-99
00196             if (rltest.ne.rltest) rltest=-1.d-9
00197             if ((istop.EQ.1).or.(seuil2 .le. seuil))then
00198                 startsauv(:,i)=starttest(:)
00199                 startobs2=starttest(:)
00200                 rlobs2=rltest
00201             else
00202                 if (rltest .ge. rlobs2) then
00203                     rlobs2=rltest
00204                     startobs2=starttest(:)
00205                 end if
00206                 !
00207                 ! *** Random optimization
00208                 !
00209                 typeOptyimization=3
00210                 compteur=0
00211                 if ((LOG_PRINT).and.(writefuncpaFichier)) write(0,*)i,"Random optimization"
00212                 rltest=0.d0
00213                 starttest(:)=0.d0
00214                 do j=1,ntest
00215                     do k=1,ndim
00216                         starttest(k)=2*rnrnor()
00217                         if (compteur .ge. 200) starttest(k)=3*rnrnor()
00218                     end do
00219                     numeroalea=0
00220                     rdmMarcq=100.0D0
00221                     call marquardt(starttest(:),ndim,maxit,v,rltest,istop,nimax,rdmMarcq)
00222                     if((istop.NE.1).and.(seuil2 .le. seuil))rdmMarcq=-99
00223                     if (rltest.ne.rltest) rltest=-1.d-10
00224                     compteur=compteur+1
00225                     if (rltest.ge. rlobs2) then
00226                         rlobs1=rlobs2
00227                         rlobs2=rltest
00228                         startobs2(:)=starttest(:)
00229                     end if
00230                     if (seuil2 .le. seuil) then
00231                         startsauv(:,i)=startobs2(:)
00232                         exit
00233                     end if
00234                 end do
00235             end if
00236         end if
00237 
00238         !
00239         ! *** Derivative computation
00240         !
00241         startsauv(:,i)=startobs2(:)
00242         scale=0.d0
00243         scaleinv=0.d0
00244         scaleinv2=0.d0
00245         adaptive=.false.
00246         adaptive2=.false.
00247 
00248         call derivVRAISTOT(startsauv(:,i),ndim,v,rl)
00249 
00250         if (v(1) .eq. 0.d0) then
00251             if ((LOG_PRINT).and.(writefuncpaFichier)) write(0,*)"Null choleski => Identity taken"
00252             v=0.d0
00253             do k=1,ndim
00254                 j=k*(k+1)/2
00255                 v(j)= v(j)+1
00256             end do
00257         else
00258             do k=1,ndim
00259                 j=k*(k+1)/2
00260                 v(j)= v(j)+1
00261             end do
00262         end if
00263         v2=v
00264         nql=0
00265         call dchole(v,ndim,nql,idpos)
00266         if (idpos .ne. 0) then
00267             call inflateDiag(v2,ndim)
00268             v=v2
00269         end if
00270 
00271         scale=0.d0
00272         ij=1
00273         do j=1,ndim
00274             do k=1,j
00275                 scale(j,k)=v(ij)
00276                 if (j .ne. k) scale(k,j)=0.d0
00277                 ij=ij+1
00278             end do
00279         end do
00280         scale=transpose(scale)
00281         deter=1
00282         do j=1,ndim
00283             deter=scale(j,j)*deter
00284         end do
00285         scaleinv=0.d0
00286         errorflag=0
00287         call FINDInv(scale, scaleinv, ndim, errorflag)
00288         if (errorflag.EQ.-1) then
00289             write(999,*)'Random effects var-cov mat. inversion FAILURE patient ',i
00290             if ((LOG_PRINT)) write(0,*)'Random effects var-cov mat. inversion FAILURE patient ',i
00291         end if
00292 
00293         v=0.d0
00294         scaleinv2=matmul(scale,transpose(scale))
00295         ij=1
00296         do j=1,ndim
00297             do k=1,j
00298                 v(ij)=scaleinv2(j,k)
00299                 ij=ij+1
00300             end do
00301         end do
00302 
00303         ier=0
00304         call dsinv(v,ndim,ep,ier,det)
00305         if (ier .ne. 0) then
00306             if ((LOG_PRINT)) write(0,*)"echec inversion funcpa",i
00307         end if
00308         ij=1
00309         do j=1,ndim
00310             do k=1,j
00311                 scaleinv2(j,k)=-v(ij)
00312                 scaleinv2(k,j)=scaleinv2(j,k)
00313                 ij=ij+1
00314             end do
00315         end do
00316 
00317         do j=1,ndim
00318             scaleinv2(j,j)=1+scaleinv2(j,j)
00319         end do
00320 
00321         scalesauv(:,:,i)=scale(:,:)
00322         scaleinvsauv(:,:,i)=scaleinv(:,:)
00323         scaleinv2sauv(:,:,i)=scaleinv2(:,:)
00324         detersauv(i)=deter
00325 
00326         !
00327         ! *** Individual likelihood computation : integration over random effects
00328         !
00329         call vraisobs(b1,i,npm2)
00330 
00331         if (ISTOP.NE.1)likelihoodERROR(1,i)=1
00332         likelihoodPRECISION(1,i)=rdmMarcq
00333 
00334 
00335         if(writefuncpaFichier)write(0, &
00336             '(i6,"  Pat.",i6,"  LL",f20.6,"  ERR.",i6,i6,i6,"     PREC.", ES18.5, ES18.5, ES18.5,"  CONTR.",i6)')  &
00337             i,idpat(i),dlog(vrais_obs(i)),likelihoodERROR(:,i),likelihoodPRECISION(:,i),sum(nbobs(i,:))
00338         flush(0)
00339     end do
00340 
00341     if (MPIutilisation.EQ.1) then
00342         call synchroFUNCPA
00343     end if
00344 
00345     if(numproc.EQ.0) then
00346         if(writefuncpaFichier) then
00347             write(999,*)"  NUM  IDPAT       Log-Likelihood                   OPT   ODE   AGQ",&
00348                 "               OPT               ODE               AGQ                    NBOBS "
00349             do i=1, nbpatient
00350                 write(999,'(i6,"  Pat.",i6,"  LL",f20.6,"  ERR.",i6,i6,i6,"     PREC.", ES18.5, ES18.5, ES18.5,"  CONTR.",i6)')  &
00351                     i,idpat(i),dlog(vrais_obs(i)),likelihoodERROR(:,i),likelihoodPRECISION(:,i),sum(nbobs(i,:))
00352             end do
00353             if (estimationWanted)flush(999)
00354         end if
00355     end if
00356 
00357 
00358     call logLikelihood_penalization(b,pena)
00359 
00360     !
00361     ! *** LL computation over all patients
00362     !
00363     funcpa=0.0D0
00364     do i=1,nbpatient
00365         result=vrais_obs(i)
00366         if (result .le. 0.d0) then
00367             funcpa=-1.d9
00368             exit
00369         else
00370             funcpa=funcpa+log(result)
00371             recap(i)=log(result)
00372         end if
00373         if (funcpa.ne.funcpa) then
00374             funcpa=-1.d9
00375             exit
00376         end if
00377     end do
00378     loglike(1)=funcpa   ! NPLL ALL PATIENTS
00379     funcpa=funcpa+pena
00380     loglike(2)=funcpa   ! PLL ALL PATIENTS
00381     !
00382     ! *** LL computation over patients where all computations went good
00383     !
00384     LLnonPenalisee=loglike(1)
00385     LLPenalisee=loglike(2)
00386     if(withexclusion) then
00387         funcpa=0.0D0
00388         if(firstFuncpa) then
00389             k=1
00390             do i=1,nbpatient
00391                 result=vrais_obs(i)
00392                 ! Optimisation fails and RDM > absoluteRequestedToleranceOPT
00393                 if((likelihoodERROR(1,i).EQ.1).and. &
00394                     ((likelihoodPRECISION(1,i).GT.absoluteRequestedToleranceOPT) &
00395                     .or.(likelihoodPRECISION(1,i).NE.likelihoodPRECISION(1,i)))) then
00396                     patOK(i)=0
00397                     listIdPatExcluded(k)=idpat(i)
00398                     k=k+1
00399                     go to 943
00400                 end if
00401                 ! ODE fails and RDM > absoluteRequestedToleranceODE
00402                 if((likelihoodERROR(2,i).EQ.1).and. &
00403                     ((likelihoodPRECISION(2,i).GT.absoluteRequestedToleranceOPT) &
00404                     .or.(likelihoodPRECISION(2,i).NE.likelihoodPRECISION(2,i)))) then
00405                     patOK(i)=0
00406                     listIdPatExcluded(k)=idpat(i)
00407                     k=k+1
00408                     go to 943
00409                 end if
00410                 ! AGQ fails and RDM > absoluteRequestedToleranceAGD
00411                 if((likelihoodERROR(3,i).EQ.1).and. &
00412                     (likelihoodPRECISION(3,i).GT.absoluteRequestedToleranceAGQ) &
00413                     .or.((likelihoodPRECISION(3,i).NE.likelihoodPRECISION(3,i)))) then
00414                     patOK(i)=0
00415                     listIdPatExcluded(k)=idpat(i)
00416                     k=k+1
00417                     go to 943
00418                 end if
00419                 ! Individual LL undefined OR TOO SMALL
00420                 if((result.NE.result).or.(result.LE.1d-300)) then
00421                     patOK(i)=0
00422                     listIdPatExcluded(k)=idpat(i)
00423                     k=k+1
00424                     go to 943
00425                 end if
00426 943         continue
00427             end do
00428             nbpatOK=sum(patOK(1:nbpatient))
00429         end if
00430 
00431         do i=1,nbpatient
00432             result=vrais_obs(i)
00433             funcpa=funcpa+log(result)*dble(patOK(i))
00434             recap(i)=log(result)*dble(patOK(i))
00435             if (funcpa.ne.funcpa) then
00436                 funcpa=-1.d10
00437                 exit
00438             end if
00439         end do
00440         loglike(3)=funcpa   ! NPLL ALL PATIENTS
00441         funcpa=funcpa+pena
00442         loglike(4)=funcpa   ! PLL ALL PATIENTS
00443         LLnonPenalisee=loglike(3)
00444         LLPenalisee=loglike(4)
00445 
00446     end if
00447 
00448     if(withexclusion) then
00449         if(loglike(3).EQ.0.D0) funcpa=-1.D9
00450     end if
00451     if(funcpa.NE.funcpa) funcpa=-1.D9
00452 
00453 
00454     if((numproc.EQ.0)) then
00455         if(LOG_PRINT) then
00456             write(0,*)" *** OVER ALL PATIENTS *** "
00457             write(0,*)'Log-likelihood : Unpenalized=',loglike(1),'Penalized=',loglike(2), &
00458                 'Penalization=',pena
00459             write(0,*)'Standardized Log-likelihood : Unpenalized=',loglike(1)/nbpatient, &
00460                 'Penalized=',loglike(2)/nbpatient
00461             if(withexclusion) then
00462                 write(0,*)" *** FOR PATIENTS WITHOUT COMPUTATION ERRORS *** ",(1-(dble(nbpatOK)/dble(nbpatient)))*100, &
00463                     " % excluded","(",(nbpatient-nbpatOK),"pat.)"
00464                 if((1-(dble(nbpatOK)/dble(nbpatient))).NE.0.0D0) then
00465                     write(0,*)"Patients excluded :",(listIdPatExcluded(i),i=1,(nbpatient-nbpatOK))
00466                     write(0,*)'Log-likelihood : Unpenalized=',loglike(3),'Penalized=',loglike(4), &
00467                         'Penalization=',pena
00468                     write(0,*)'Standardized Log-likelihood : Unpenalized=',loglike(3)/nbpatOK, &
00469                         'Penalized=',loglike(3)/nbpatOK
00470                 end if
00471             else
00472                 write(0,*)'--> No exclusion requested'
00473             end if
00474             flush(0)
00475         end if
00476         write(999,*)" *** OVER ALL PATIENTS *** "
00477         write(999,*)'Log-likelihood : Unpenalized=',loglike(1),'Penalized=',loglike(2), &
00478             'Penalization=',pena
00479         write(999,*)'Standardized Log-likelihood : Unpenalized=',loglike(1)/nbpatient, &
00480             'Penalized=',loglike(2)/nbpatient
00481         if(withexclusion) then
00482             write(999,*)" *** FOR PATIENTS WITHOUT COMPUTATION ERRORS *** ",(1-(dble(nbpatOK)/dble(nbpatient)))*100, &
00483                 " % excluded","(",(nbpatient-nbpatOK),"pat.)"
00484             if((1-(dble(nbpatOK)/dble(nbpatient))).NE.0.0D0) then
00485                 write(999,*)"Patients excluded :",(listIdPatExcluded(i),i=1,(nbpatient-nbpatOK))
00486                 write(999,*)'Log-likelihood : Unpenalized=',loglike(3),'Penalized=',loglike(4), &
00487                     'Penalization=',pena
00488                 write(999,*)'Standardized Log-likelihood : Unpenalized=',loglike(3)/nbpatOK, &
00489                     'Penalized=',loglike(3)/nbpatOK
00490             end if
00491         else
00492             write(999,*)'--> No exclusion requested'
00493         end if
00494         flush(999)
00495     end if
00496 
00497 
00498     call cpu_time (tpscpu)
00499     if (numproc.EQ.0) write(0,*)"TIME : ",tpscpu
00500     flush(0)
00501 
00502 
00503     return
00504 end function funcpa
00505 
00506 
00507 
00508 !------------------------------------------------------------------------------
00509 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00510 !                effects based on Ordinary Differential equations
00511 !------------------------------------------------------------------------------
00512 !
00513 ! VERSION : 1.0
00514 !
00515 ! MODULE: funcpa.f90 function funcpaRandomEffect
00541 function funcpaRandomEffect(b,m,id,thi,jd,thj)
00542     use WorkingSharedValues
00543     implicit none
00544 
00545     double precision::funcpaRandomEffect
00546     double precision,dimension(m),intent(inout)::b
00547     integer,intent(in) :: jd,id
00548     double precision,intent(in) :: thi,thj
00549     integer,intent(in)::m
00550 
00551     double precision,dimension(ndim)::b2,b3
00552     double precision ::funvls
00553     integer :: k
00554 
00555 
00556     b2=0.d0
00557     funcpaRandomEffect=0.d0
00558     b3=0.d0
00559 
00560     do k=1,m
00561         b2(k)=b(k)
00562     end do
00563     if (id.ne.0) b2(id)=b2(id)+thi
00564     if (jd.ne.0) b2(jd)=b2(jd)+thj
00565     if (numeroalea .ne. 0) then
00566         if(debugFUNCPA2)print*,"b3 (dim), numeroalea",b3,numeroalea
00567         b3(numeroalea)=b2(1)
00568     else
00569         b3=b2
00570     end if
00571     if(debugFUNCPA2)print*,"avant vraistot",b3,numeroalea
00572     call vraistotEXP(m,b3,nf,FUNVLS)
00573     if(debugMARC)write(0,*)"apres vraistot",b3,FUNVLS
00574     if(debugFUNCPA2)print*,"apres vraistot",FUNVLS
00575 
00576     funcpaRandomEffect=FUNVLS+log(exp(-(5.d-1)*dot_product(b2,b2)))
00577     if (funvls.ne.funvls) then
00578         funvls=-1.d9
00579         funcpaRandomEffect=-1.d9
00580     end if
00581     seuil2=abs(FUNVLS)/sum(nbobs(numpat1,:))
00582     if(debugFUNCPA2)print*,"finFUNCPA2"
00583     return
00584 end function funcpaRandomEffect