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

derivation.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: derivation.f90 function derivRVS
00050 subroutine derivRVS(b,m,v,rl)
00051     use WorkingSharedValues
00052     use mpimod
00053     USE m_vpropre
00054 
00055     IMPLICIT NONE
00056     double precision, dimension(m),intent(inout)::b
00057     integer,intent(in)::m
00058     double precision,dimension(m*(m+3)/2),intent(out) :: v
00059     double precision,intent(out)::rl
00060 
00061     double precision, dimension(m) :: pena_grad,pena_hess
00062     double precision :: funcpa,z
00063     double precision,dimension(nbpatienta,m) :: uscore
00064     double precision,dimension(m):: scorederivb
00065     double precision,dimension(m,m):: variance, variance2
00066     integer i,k,k0,j,iii,ij
00067 
00068 
00069     double precision,dimension(m,m)::v2,vecp
00070     double precision,dimension(m) ::VP,D,WI
00071     integer, dimension(m) :: intc
00072 
00073     z=0.d0
00074     k0=0
00075     if ((numproc.EQ.0).and.(LOG_PRINT))write(0,*)" -> Likelihood Computation"
00076     writefuncpaFichier=.true.
00077     firstFuncpa=.true.
00078     rl=funcpa(b,m, k0,z,k0,z)
00079     writefuncpaFichier=.false.
00080     firstFuncpa=.false.
00081     v=0.d0
00082 
00083 
00084     if(debugSCOREHESSIANfichier) then
00085      if ((numproc.EQ.0))open(546,file="debugSCORE.txt")
00086     end if
00087 
00088     !
00089     ! *** Scores computation
00090     !
00091     if ((numproc.EQ.0).and.(LOG_PRINT))write(0,*)" -> Scores Computation"
00092     do i=1,npmbio
00093         abserr1=abserr(i)
00094         call scoreRVS(b,i,uscore(:,i),1)
00095         where ( (uscore .NE.uscore)  .or.(abs(uscore) .ge. 1.d4))
00096             uscore=0.d0
00097         end where
00098 
00099         if ((numproc.EQ.0).and.(LOG_PRINT)) write(0,'(A,i3,A,f20.8,A,f10.3,A,ES14.7)')"ScoreBIOL( ",i," ) ", &
00100             sum(uscore(1:nbpatient,i)),"       prec. not reached=",sum(scoreERROR(2,i,1:nbpatient))/dble(nbpatient)*100.D0, &
00101             "%   mean accuracy when fails=",MAX(0.0D0,DOT_PRODUCT(scorePRECISION(2,i,1:nbpatient),scoreERROR(2,i,1:nbpatient)) &
00102             /sum(scoreERROR(2,i,1:nbpatient)))
00103         flush(0)
00104 
00105 
00106         if ((numproc.EQ.0)) write(999,'(A,i3,A,f20.8,A,f10.3,A,ES14.7)')"ScoreBIOL( ",i," ) ", &
00107             sum(uscore(1:nbpatient,i)),"       prec. not reached=",sum(scoreERROR(2,i,1:nbpatient))/dble(nbpatient)*100.D0, &
00108             "%   mean accuracy when fails=",MAX(0.0D0,DOT_PRODUCT(scorePRECISION(2,i,1:nbpatient),scoreERROR(2,i,1:nbpatient)) &
00109             /sum(scoreERROR(2,i,1:nbpatient)))
00110         flush(999)
00111     end do
00112 
00113 
00114     do i=npmbio+1,npmbio+npmexpl
00115         abserr1=abserr(i)
00116         call scoreRVS(b,i,uscore(:,i),1)
00117         where ( (uscore .NE.uscore)  .or.(abs(uscore) .ge. 1.d4))
00118             uscore=0.d0
00119         end where
00120 
00121         if ((numproc.EQ.0).and.(LOG_PRINT)) write(0,'(A,i3,A,f20.8,A,f10.3,A,ES14.7)')"ScoreBIOL( ",i," ) ", &
00122             sum(uscore(1:nbpatient,i)),"       prec. not reached=",sum(scoreERROR(2,i,1:nbpatient))/dble(nbpatient)*100.D0, &
00123             "%   mean accuracy when fails=",MAX(0.0D0,DOT_PRODUCT(scorePRECISION(2,i,1:nbpatient),scoreERROR(2,i,1:nbpatient)) &
00124             /sum(scoreERROR(2,i,1:nbpatient)))
00125 
00126         if ((numproc.EQ.0)) write(999,'(A,i3,A,f20.8,A,f10.3,A,ES14.7)')"ScoreBIOL( ",i," ) ", &
00127             sum(uscore(1:nbpatient,i)),"       prec. not reached=",sum(scoreERROR(2,i,1:nbpatient))/dble(nbpatient)*100.D0, &
00128             "%   mean accuracy when fails=",MAX(0.0D0,DOT_PRODUCT(scorePRECISION(2,i,1:nbpatient),scoreERROR(2,i,1:nbpatient)) &
00129             /sum(scoreERROR(2,i,1:nbpatient)))
00130         flush(999)
00131     end do
00132     do i=1,ndim
00133         abserr1=abserr(i+npmbio+npmexpl)
00134         call scoreRVS(b,i,uscore(:,i+npmbio+npmexpl),2)
00135         where ( (uscore .NE.uscore)  .or.(abs(uscore) .ge. 1.d4))
00136             uscore=0.d0
00137         end where
00138 
00139         if ((numproc.EQ.0).and.(LOG_PRINT)) write(0,'(A,i3,A,f20.8,A,f10.3,A,ES14.7)')"ScoreBIOL( ",i+npmbio+npmexpl," ) ", &
00140             sum(uscore(1:nbpatient,i+npmbio+npmexpl)),"       prec. not reached=", &
00141             sum(scoreERROR(2,i+npmbio+npmexpl,1:nbpatient))/dble(nbpatient)*100.D0, &
00142             "%   mean accuracy when fails=",MAX(0.0D0,DOT_PRODUCT(scorePRECISION(2,i+npmbio+npmexpl,1:nbpatient), &
00143             scoreERROR(2,i+npmbio+npmexpl,1:nbpatient)) &
00144             /sum(scoreERROR(2,i+npmbio+npmexpl,1:nbpatient)))
00145 
00146         if ((numproc.EQ.0)) write(999,'(A,i3,A,f20.8,A,f10.3,A,ES14.7)')"ScoreBIOL( ",i+npmbio+npmexpl," ) ", &
00147             sum(uscore(1:nbpatient,i+npmbio+npmexpl)),"       prec. not reached=", &
00148             sum(scoreERROR(2,i+npmbio+npmexpl,1:nbpatient))/dble(nbpatient)*100.D0, &
00149             "%   mean accuracy when fails=",MAX(0.0D0,DOT_PRODUCT(scorePRECISION(2,i+npmbio+npmexpl,1:nbpatient), &
00150             scoreERROR(2,i+npmbio+npmexpl,1:nbpatient)) &
00151             /sum(scoreERROR(2,i+npmbio+npmexpl,1:nbpatient)))
00152         flush(999)
00153     end do
00154     do i=1,npmcomp
00155         uscore(:,m-npmcomp+i)=0.d0
00156         abserr1=abserr(m-npmcomp+i)
00157         call scoreRVS(b,i,uscore(:,m-npmcomp+i),3)
00158         where ( (uscore .NE.uscore)  .or.(abs(uscore) .ge. 1.d4).or.(abs(uscore) .le. -1.d10))
00159             uscore=0.d0
00160         end where
00161 
00162         if ((numproc.EQ.0).and.(LOG_PRINT)) write(0,'(A,i3,A,f20.8,A,f10.3,A,ES14.7)')"ScoreBIOL( ",m-npmcomp+i," ) ", &
00163             sum(uscore(1:nbpatient,m-npmcomp+i)),"       prec. not reached=", &
00164             sum(scoreERROR(2,m-npmcomp+i,1:nbpatient))/dble(nbpatient)*100.D0, &
00165             "%   mean accuracy when fails=",MAX(0.0D0,DOT_PRODUCT(scorePRECISION(2,m-npmcomp+i,1:nbpatient), &
00166             scoreERROR(2,m-npmcomp+i,1:nbpatient)) &
00167             /sum(scoreERROR(2,m-npmcomp+i,1:nbpatient)))
00168 
00169         if ((numproc.EQ.0)) write(999,'(A,i3,A,f20.8,A,f10.3,A,ES14.7)')"ScoreBIOL( ",m-npmcomp+i," ) ", &
00170             sum(uscore(1:nbpatient,m-npmcomp+i)),"       prec. not reached=", &
00171             sum(scoreERROR(2,m-npmcomp+i,1:nbpatient))/dble(nbpatient)*100.D0, &
00172             "%   mean accuracy when fails=",MAX(0.0D0,DOT_PRODUCT(scorePRECISION(2,m-npmcomp+i,1:nbpatient), &
00173             scoreERROR(2,m-npmcomp+i,1:nbpatient)) &
00174             /sum(scoreERROR(2,m-npmcomp+i,1:nbpatient)))
00175         flush(999)
00176 
00177     end do
00178 
00179 
00180     if(withexclusion) then
00181         do i=1,nbpatient
00182             uscore(i,:)=uscore(i,:)*dble(patOK(i))
00183         end do
00184     end if
00185 
00186     do i=1,m
00187         scorederivb(i)=sum(uscore(:,i))
00188     end do
00189 
00190     if(debugSCOREHESSIANfichier) then
00191     write(546,*) 'LES SCORES'
00192     do i=1,m
00193         if ((numproc.EQ.0).and.(LOG_PRINT)) write(546,'(i5,1x,200(D20.12,",",1x))'),i,uscore(1:nbpatient,i)
00194     end do
00195     end if
00196     !
00197     ! *** Hessian computation
00198     !
00199     if(withexclusion) then
00200         if(nbpatOK.GT.0) then
00201             k=1
00202             do i=1,m
00203                 do j=1,i
00204                     variance(i,j)=sum(uscore(:,i)*uscore(:,j))
00205                     variance2(i,j)=- (1.d0/nbpatOK)*sum(uscore(:,i)) &
00206                         *sum(uscore(:,j))
00207                     v(k)=variance(i,j)+variance2(i,j)
00208                     k=k+1
00209                 end do
00210             end do
00211         else
00212             v(1:(m*(m+1)/2))=0.0D0
00213             do i=1,m
00214                 v((i*(i+1)/2))=1.0D0
00215             end do
00216         end if
00217     else
00218         k=1
00219         do i=1,m
00220             do j=1,i
00221                 variance(i,j)=sum(uscore(:,i)*uscore(:,j))
00222                 variance2(i,j)=- (1.d0/nbpatient)*sum(uscore(:,i)) &
00223                     *sum(uscore(:,j))
00224                 v(k)=variance(i,j)+variance2(i,j)
00225                 k=k+1
00226             end do
00227         end do
00228     end if
00229 
00230     !
00231     ! *** Valeurs propres
00232     !
00233     do i=1,m
00234         do j=i,m
00235             ij=(j-1)*j/2+i
00236             V2(i,j)=v(ij)
00237             V2(j,i)=v(ij)
00238         end do
00239     END DO
00240 
00241     if(debugSCOREHESSIANfichier) then
00242     write(546,*) 'je suis sur V2'
00243     do i=1,m
00244         if ((numproc.EQ.0)) write(546,'(200(ES30.15,",",1x))') V2(i,1:m)
00245     end do
00246     end if
00247 
00248     do i=1,m
00249         v(i+(m*(m+1)/2))=scorederivb(i)
00250     end do
00251 
00252     !
00253     ! *** MAP if necessary
00254     !
00255     hessianNonPenalisee=v
00256     call gradients_penalization(b,pena_grad)
00257     call hessian_penalization(b,pena_hess)
00258     do i=1,npmbio+npmexpl+ndim+npmcomp
00259         v(i+(m*(m+1)/2))=v(i+(m*(m+1)/2))+pena_grad(i)
00260         k=i*(i+1)/2
00261         v(k)=v(k)-pena_hess(i)
00262     end do
00263     hessianPenalisee=v
00264 
00265     if(debugSCOREHESSIANfichier) then
00266     write(546,*) 'LES SCORES PENALISEES'
00267     do i=1,m
00268         if ((numproc.EQ.0).and.(LOG_PRINT)) write(546,'(i5,1x,200(D20.12,",",1x))'),i,uscore(1:nbpatient,i)+pena_grad(i)/35.d0
00269     end do
00270     end if
00271 
00272     do i=1,m
00273         do j=i,m
00274             ij=(j-1)*j/2+i
00275             V2(i,j)=v(ij)
00276             V2(j,i)=v(ij)
00277         end do
00278     END DO
00279 
00280     if(debugSCOREHESSIANfichier) then
00281     write(546,*) 'Penalization'
00282     do i=1,m
00283         if ((numproc.EQ.0)) write(546,'(200(ES30.15,",",1x))') V2(i,1:m)
00284     end do
00285     end if
00286     if(debugSCORE) then
00287      if ((numproc.EQ.0))close(546)
00288      pause
00289     end if
00290     return
00291 end subroutine derivRVS
00292 
00293 !------------------------------------------------------------------------------
00294 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00295 !                effects based on Ordinary Differential equations
00296 !------------------------------------------------------------------------------
00297 !
00298 ! VERSION : 1.0
00299 !
00300 ! MODULE: derivation.f90 function derivMARC
00325 subroutine derivMARC(b,m,v,rl)
00326     use WorkingSharedValues
00327     use mpimod
00328     IMPLICIT NONE
00329 
00330 
00331     integer,intent(in)::m
00332     double precision,intent(inout)::rl
00333     double precision,dimension(m),intent(in)::b
00334     double precision,dimension((m*(m+3)/2)),intent(out)::v
00335     double precision,dimension(m)::b3,fcith
00336     integer ::i0,m1,ll,i,k,j
00337     double precision::thn,th,z,vl,temp,thi,thj
00338     double precision,external::funcpa
00339     double precision, dimension(m) :: pena_grad,pena_hess
00340 
00341     th=1.d-5
00342     thn=-th
00343     z=0.d0
00344     i0=0
00345     writefuncpaFichier=.true.
00346     firstFuncpa=.true.
00347     rl=funcpa(b,m,i0,z,i0,z)
00348     firstFuncpa=.false.
00349     writefuncpaFichier=.false.
00350 
00351 
00352     do i=1,m
00353         th=DMAX1(1.d-7, 1.d-4 * DABS(b(i)))
00354         b3=b
00355         fcith(i)=funcpa(b3,m,i,th,i0,z)
00356         fth(:,i)=recap(:)
00357         thn=-DMAX1(1.d-7, 1.d-4 * DABS(b(i)))
00358         b3=b
00359         fcith(i)=funcpa(b3,m,i,thn,i0,z)
00360         fthn(:,i)=recap(:)
00361         vl=sum(fth(:,i)-fthn(:,i))/(2.d0*(th))
00362         score(:,i)=(fth(:,i)-fthn(:,i))/(2.d0*(th))
00363         if (numproc.EQ.0) write(0,*)"derivee du prm",i,vl
00364         v((m*(m+1)/2)+i)=vl
00365     end do
00366 
00367     k=0
00368     ll=m*(m+1)/2
00369     do i=1,m
00370         do j=1,i
00371             k=k+1
00372             v(k)=DOT_PRODUCT(score(:,i),score(:,j))
00373             v(k)=v(k)-(1.d0/nbpatient)*v((m*(m+1)/2)+i)*v((m*(m+1)/2)+j)
00374         end do
00375     end do
00376 
00377     !
00378     ! *** Penalization
00379     !
00380     hessianNonPenalisee=v
00381     call gradients_penalization(b,pena_grad)
00382     call hessian_penalization(b,pena_hess)
00383 
00384     do i=1,npmbio+npmexpl+ndim+npmcomp
00385         v(i+(m*(m+1)/2))=v(i+(m*(m+1)/2))+pena_grad(i)
00386         k=i*(i+1)/2
00387         v(k)=v(k)-pena_hess(i)
00388     end do
00389     hessianPenalisee=v
00390     return
00391 
00392 end subroutine derivMARC
00393 
00394 
00395 !------------------------------------------------------------------------------
00396 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00397 !                effects based on Ordinary Differential equations
00398 !------------------------------------------------------------------------------
00399 !
00400 ! VERSION : 1.0
00401 !
00402 ! MODULE: derivation.f90 function derivVRAISTOT
00425 subroutine derivVRAISTOT(b,m,v,rl)
00426     use VectorLength
00427     implicit none
00428     double precision, dimension(m),intent(inout)::b
00429     double precision, dimension(m*(m+3)/2),intent(out)::v
00430     integer,intent(in)::m
00431     double precision,intent(out)::rl
00432 
00433     double precision funvls,thn,th,vl,thi,thj
00434     double precision,dimension(m) :: b1,fcith
00435     integer i0,ll,i,k,j,ii
00436 
00437     call vraistotEXP(NDIM,b,NF,FUNVLS)
00438     rl=funvls
00439     do i=1,m
00440         th=DMAX1(1.d-7, 1.d-4 * DABS(b(i)))
00441         b1=b
00442         b1(i)=b(i)+th
00443         call vraistotEXP(NDIM,b1,NF,FUNVLS)
00444         fcith(i)=funvls
00445     end do
00446     k=0
00447 
00448     ll=m*(m+1)/2
00449     do i=1,m
00450         ll=ll+1
00451         thn=-DMAX1(1.d-7, 1.d-4 * DABS(b(i)))
00452         b1=b
00453         b1(i)=b(i)+thn
00454         call vraistotEXP(NDIM,b1,NF,FUNVLS)
00455         vl=(fcith(i)-FUNVLS)/(2.d0*(-thn))
00456         v(ll)=vl
00457         do j=1,i
00458             k=k+1
00459             thi=DMAX1(1.d-7, 1.d-4 * DABS(b(i)))
00460             thj=DMAX1(1.d-7, 1.d-4 * DABS(b(j)))
00461             b1=b
00462             b1(i)=b(i)+thi
00463             b1(j)=b1(j)+thj
00464             call vraistotEXP(NDIM,b1,NF,FUNVLS)
00465             v(k)=-(FUNVLS-fcith(j)-fcith(i)+rl)/(thi*thj)
00466         end do
00467     end do
00468 
00469     return
00470 end subroutine derivVRAISTOT
00471 
00472 
00473 
00474 !------------------------------------------------------------------------------
00475 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00476 !                effects based on Ordinary Differential equations
00477 !------------------------------------------------------------------------------
00478 !
00479 ! VERSION : 1.0
00480 !
00481 ! MODULE: derivation.f90 function derivRANDOMEFFECT
00503 subroutine derivRANDOMEFFECT(b,m,v,rl)
00504     use VectorLength
00505     implicit none
00506 
00507     integer,intent(in)::m
00508     double precision,intent(inout)::rl
00509     double precision,dimension(m),intent(in)::b
00510     double precision,dimension((m*(m+3)/2)),intent(out)::v
00511     double precision,dimension(m)::fcith
00512     integer ::i0,m1,ll,i,k,j
00513     double precision::thn,th,z,vl,temp,thi,thj
00514     double precision,external::funcpaRandomEffect
00515 
00516     z=0.d0
00517     i0=0
00518 
00519     rl=funcpaRandomEffect(b,m,i0,z,i0,z)
00520 
00521     if(rl.eq.-1.d9) then
00522         goto 123
00523     end if
00524 
00525     do i=1,m
00526         th=DMAX1(1.d-7, 1.d-4 * DABS(b(i)))
00527         fcith(i)=funcpaRandomEffect(b,m,i,th,i0,z)
00528         if(fcith(i).eq.-1.d9) then
00529             rl=-1.d9
00530             goto 123
00531         end if
00532     end do
00533 
00534     k=0
00535     m1=m*(m+1)/2
00536     ll=m1
00537     Main:do i=1,m
00538         ll=ll+1
00539         thn=-DMAX1(1.d-7, 1.d-4 * DABS(b(i)))
00540         temp=funcpaRandomEffect( b,m, i,thn,i0,z)
00541         if(temp.eq.-1.d9) then
00542             rl=-1.d9
00543             exit Main
00544         end if
00545         vl=(fcith(i)-temp)/(2.d0*(-thn))
00546         v(ll)=vl
00547         do j=1,i
00548             k=k+1
00549 
00550             thi=DMAX1(1.d-7, 1.d-4 * DABS(b(i)))
00551             thj=DMAX1(1.d-7, 1.d-4 * DABS(b(j)))
00552 
00553             temp=funcpaRandomEffect(b,m,i,thi,j,thj)
00554             if(temp.eq.-1.d9) then
00555                 rl=-1.d9
00556                 exit Main
00557             end if
00558             v(k)=-(temp-fcith(j)-fcith(i)+rl)/(thi*thj)
00559         end do
00560     end do Main
00561 123 continue
00562     return
00563 
00564 end subroutine derivRANDOMEFFECT
00565 
00566 
00567