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

scores.f90

Go to the documentation of this file.
00001 !------------------------------------------------------------------------------
00002 !
00003 ! VERSION : 1.0
00004 !
00005 ! MODULE: scores.f90 function scoreRVS
00030 
00031 subroutine scoreRVS(b,num_coeff,uscore,typeScore)
00032     use WorkingSharedValues
00033     use mpimod
00034     implicit none
00035 
00036     double precision,dimension(npm),intent(in) :: b
00037     double precision,dimension(nbpatienta),intent(out)::uscore
00038     integer, intent(in) :: num_coeff
00039     integer,intent(in) :: typeScore
00040 
00041     integer ::i, NDIM2, NF2, MINPTS, MAXPTS, RESTAR,NEVAL, IFAIL,debut,fin
00042     DOUBLE PRECISION ::  EPSABS,EPSREL,RESULT,ABSERR2
00043     DOUBLE PRECISION,dimension(10000) :: WORK
00044     external :: FUNSUB_Biological, FUNSUB_Error, FUNSUB_Random
00045 
00046     uscore=0.0D0
00047     ndim2=ndim
00048     nf2=nf
00049     numcoeff1=num_coeff
00050     MINPTS=scores_MINPTS
00051     MAXPTS=scores_MAXPTS
00052     epsrel=scores_EPSREL
00053     call repartirSurCoeurs(debut,fin,nbpatient)
00054     scorePRECISION(:,:,:)=0.0D0
00055     scoreERROR(:,:,:)=0.0D0
00056     if((LOG_PRINT).and.(numproc.EQ.0).and.(INFO_SCORE))write(0,*)"  NUM  IDPAT       Log-Likelihood", &
00057          "                   ODE   AGQ               ODE               AGQ     "
00058     if((numproc.EQ.0).and.(INFO_SCORE))write(999,*)"  NUM  IDPAT       Log-Likelihood", &
00059          "                   ODE   AGQ               ODE               AGQ     "
00060     do i=debut,fin
00061 
00062         numpat1=i
00063 
00064         epsabs=vrais_obs(i)*abserr1!1d-3!  !! Abserr1 correspond a precision sur la variance/ (2*nbpatients)
00065 
00066         RESTAR=0
00067         adaptive=.true.
00068         uscore(i)=0.d0
00069         if (typeScore.EQ.1) then
00070             numparam=numcoeff1
00071             call  hrmsym( NDIM2, NF2, MINPTS, MAXPTS, FUNSUB_Biological, EPSABS, &
00072                 EPSREL, RESTAR, RESULT, ABSERR2, NEVAL, IFAIL, WORK)
00073         else
00074             if (typeScore.EQ.2) then
00075                 numparam=numcoeff1+npmbio+npmexpl
00076                 call  hrmsym( NDIM2, NF2, MINPTS, MAXPTS, FUNSUB_Random, EPSABS, &
00077                     EPSREL, RESTAR, RESULT, ABSERR2, NEVAL, IFAIL, WORK)
00078             else
00079                 numparam=numcoeff1+npmbio+npmexpl+ndim
00080                 call  hrmsym( NDIM2, NF2, MINPTS, MAXPTS, FUNSUB_Error, EPSABS, &
00081                     EPSREL, RESTAR, RESULT, ABSERR2, NEVAL, IFAIL, WORK)
00082             end if
00083         end if
00084 
00085         uscore(i)=result/vrais_obs(i)
00086         if(debugSCORE)write(0,*)"Decomposition uscore pat",i,epsabs,result,vrais_obs(i),IFAIL
00087 
00088         scorePRECISION(2,numparam,i)=ABSERR2
00089         IF(IFAIL.NE.0)scoreERROR(2,numparam,i)=1
00090         if((LOG_PRINT).and.(INFO_SCORE)) write(0,'(i6,"  Pat.",i6,"  LL",f20.6,"  ERR.",i6,i6,"     PREC.", ES18.5, ES18.5)')  &
00091                       i,idpat(i),uscore(i),scoreERROR(:,numparam,i),scorePRECISION(:,numparam,i)
00092     end do
00093     if (MPIutilisation.EQ.1) then
00094         call synchroCALCULSCORES(uscore)
00095     end if
00096 
00097     if((numproc.EQ.0).and.(INFO_SCORE)) then
00098        write(999,*)"  NUM  IDPAT       Log-Likelihood                   ODE   AGQ",&
00099          "               ODE               AGQ                  "
00100        do i=1, nbpatient
00101           write(999,'(i6,"  Pat.",i6,"  LL",f20.6,"  ERR.",i6,i6,"     PREC.", ES18.5, ES18.5)')  &
00102                       i,idpat(i),uscore(i),scoreERROR(:,numparam,i),scorePRECISION(:,numparam,i)
00103        end do
00104     flush(999)
00105     end if
00106 
00107 
00108     return
00109 end subroutine scoreRVS
00110 
00111 !------------------------------------------------------------------------------
00112 !
00113 ! VERSION : 1.0
00114 !
00115 ! MODULE: scores.f90 function FUNSUB_Biological
00139 
00140 subroutine FUNSUB_Biological(NDIM2,X,NF2,FUNVLS)
00141     use WorkingSharedValues
00142     use Constante
00143 
00144     implicit none
00145 
00146     integer,intent(in)::ndim2,nf2
00147     double precision,dimension(ndim2),intent(in)::X
00148     double precision,intent(out) :: funvls
00149 
00150     double precision ::vrais_tot,funvls2,funvls1,fn_val,phid
00151     double precision, dimension(tdef2,npmcomp)::simul,deriv
00152     double precision,dimension(ndim2)::Xprime,Xprime2
00153     integer :: i,j,k
00154 
00155     systeme(1)=2
00156     systeme(2)=numcoeff1
00157 
00158     !
00159     ! *** Call to the ODE solver to have sensitivity equation systems trajectories
00160     !
00161     if (adaptive) then
00162         Xprime(1:ndim2)=startsauv(1:ndim2,numpat1) + &
00163             matmul(scaleinvsauv(1:ndim2,1:ndim2,numpat1),X(1:ndim2))
00164         call solution(xprime,simul,deriv(:,1:npmcomp))
00165     else
00166         call solution(x,simul,deriv(:,1:npmcomp))
00167     end if
00168     fn_val=0.d0
00169     funvls1=0.d0
00170     vrais_tot=1.d0
00171 
00172     do k=1,npmcomp
00173         do j=1,nbobs(numpat1,k)
00174             !
00175             ! *** If there is no censorship
00176             !
00177              if (censor(numpat1,j,k).NE.1) then
00178 
00179                 !
00180                 ! *** likelihood
00181                 !
00182                 funvls1=donnees(numpat1,j,k)-&
00183                     simul(schedule(numpat1,j,k)+1,k)
00184                 funvls2= (funvls1 /b1(npm-npmcomp+k))**2
00185                 funvls2=-funvls2
00186                 funvls2=funvls2/2
00187                 funvls2=exp(funvls2)/abs(b1(npm-npmcomp+k))
00188                 funvls2=funvls2/sqrt(2.d0*pigrec)
00189                 !
00190                 ! *** For the differentiated value
00191                 !
00192                 funvls1 = funvls1 /(b1(npm-npmcomp+k)**2)
00193 
00194                 if(debugFUNSUB)write(0,*)schedule(numpat1,j,k),donnees(numpat1,j,k), &
00195                     simul(schedule(numpat1,j,k)+1,k),deriv(schedule(numpat1,j,k)+1,k)
00196             !
00197             ! *** If there is censorship
00198             !
00199             else
00200                 if(debugFUNSUB)write(0,*)"censored",schedule(numpat1,j,k),donnees(numpat1,j,k),simul(schedule(numpat1,j,k)+1,k)
00201                 funvls1=-((1.d0)/(DABS(b1(npm-npmcomp+k))*DSQRT(2*pigrec)))*DEXP(-(5.d-1)* &
00202                     (((donnees(numpat1,j,k)-simul(schedule(numpat1,j,k)+1,k))/b1(npm-npmcomp+k))**2))
00203                 funvls1=funvls1/phid((donnees(numpat1,j,k)-simul(schedule(numpat1,j,k)+1,k))/ DABS(b1(npm-npmcomp+k)))
00204                 funvls2=phid((donnees(numpat1,j,k)-simul(schedule(numpat1,j,k)+1,k))/DABS(b1(npm-npmcomp+k)))
00205                 if (funvls1 .ne. funvls1) funvls1=0.d0
00206                 if (funvls2 .ne. funvls2) funvls2=0.d0
00207             end if
00208             funvls1=funvls1*deriv(schedule(numpat1,j,k)+1,k)
00209             fn_val=fn_val + funvls1
00210             vrais_tot=vrais_tot*funvls2
00211         end do
00212     end do
00213     fn_val=fn_val*vrais_tot
00214     funvls=fn_val
00215 
00216     if (adaptive) then
00217         Xprime2(1:ndim)=matmul(scaleinvsauv(1:ndim,1:ndim,numpat1),&
00218             X(1:ndim))
00219         Xprime(1:ndim)=startsauv(1:ndim,numpat1) + Xprime2(1:ndim)
00220         funvls=funvls*exp(-(sum(startsauv(1:ndim,numpat1)**2) )/2)/&
00221             detersauv(numpat1)
00222         funvls=funvls*exp(-(dot_product(Xprime2(1:ndim),&
00223             startsauv(1:ndim,numpat1))))
00224         funvls=funvls*&
00225             exp((5.d-1)*dot_product(X(1:ndim),&
00226             matmul(scaleinv2sauv(1:ndim,1:ndim,numpat1),X(1:ndim))))
00227         !if (abs(funvls) <1.d-100) funvls=1.d-100
00228     end if
00229     funvls=funvls
00230     !if (abs(funvls) <1.d-100) funvls=1.d-100
00231     if (funvls.NE.funvls) funvls=1.d-100
00232     return
00233 end subroutine  FUNSUB_Biological
00234 
00235 
00236 !------------------------------------------------------------------------------
00237 !
00238 ! VERSION : 1.0
00239 !
00240 ! MODULE: scores.f90 function FUNSUB_Random
00264 subroutine FUNSUB_Random(NDIM2,X,NF2,FUNVLS)
00265     use WorkingSharedValues
00266     use Constante
00267 
00268     implicit none
00269 
00270     integer,intent(in)::ndim2,nf2
00271     double precision,dimension(ndim2),intent(in)::X
00272     double precision,intent(out) :: funvls
00273 
00274     double precision ::vrais_tot,funvls3,funvls2,funvls1,fn_val,phid
00275     double precision, dimension(tdef2,npmcomp)::simul,deriv
00276     double precision,dimension(ndim2)::Xprime,Xprime2
00277     integer :: i,j,k
00278     simul=0.d0
00279     systeme(1)=2
00280     systeme(2)=numcoeff1
00281 
00282     !
00283     ! *** Call to the ODE solver to have sensitivity equation systems trajectories
00284     !
00285     if (adaptive) then
00286         Xprime(1:ndim)=startsauv(1:ndim,numpat1) + &
00287             matmul(scaleinvsauv(1:ndim,1:ndim,numpat1),X(1:ndim))
00288         call solution(xprime,simul,deriv)
00289     else
00290         call solution(x,simul,deriv)
00291     end if
00292 
00293     fn_val=0.d0
00294     funvls1=0.d0
00295     vrais_tot=1.d0
00296     do k=1,npmcomp
00297         do j=1,nbobs(numpat1,k)
00298             !
00299             ! *** If there is no censorship
00300             !
00301              if (censor(numpat1,j,k).NE.1) then
00302                 funvls1=0.d0
00303                 funvls1=donnees(numpat1,j,k) - simul(schedule(numpat1,j,k)+1,k)
00304                 if(debugFUNSUB)write(0,*)schedule(numpat1,j,k),donnees(numpat1,j,k),simul(schedule(numpat1,j,k)+1,k)
00305                 funvls3= (funvls1 /b1(npm-npmcomp+k))**2
00306                 funvls1 = funvls1 /(b1(npm-npmcomp+k)**2)
00307                 funvls3=-funvls3
00308                 funvls3=funvls3/2
00309                 funvls3=exp(funvls3)/dabs(b1(npm-npmcomp+k))
00310                 funvls3=funvls3/dsqrt(2.d0*pigrec)
00311             else
00312                 if(debugFUNSUB)write(0,*)"censored",schedule(numpat1,j,k),donnees(numpat1,j,k),simul(schedule(numpat1,j,k)+1,k)
00313                 funvls1=0.d0
00314                 funvls2=donnees(numpat1,j,k)-simul(schedule(numpat1,j,k)+1,k)
00315                 funvls2 = funvls2 /DABS(b1(npm-npmcomp+k))
00316                 funvls1=-((1.d0)/DABS(b1(npm-npmcomp+k)))/phid(funvls2)
00317                 funvls2=funvls2**2
00318                 funvls2=-funvls2/2
00319                 funvls1=((1.d0)/(DSQRT(2*pigrec)))*funvls1*DEXP(funvls2)
00320                 funvls3=phid((donnees(numpat1,j,k)-simul(schedule(numpat1,j,k)+1,k))/DABS(b1(npm-npmcomp+k)))
00321                 if (funvls1 .ne. funvls1) funvls1=0.d0
00322                 if (funvls2 .ne. funvls2) funvls2=0.d0
00323                 if (funvls3 .ne. funvls3) funvls3=0.d0
00324             end if
00325 
00326             vrais_tot=vrais_tot*funvls3
00327             funvls1=funvls1*deriv(schedule(numpat1,j,k)+1,k)
00328             fn_val=fn_val + funvls1
00329         end do
00330     end do
00331 
00332     if (adaptive) then
00333         fn_val=fn_val*Xprime(numcoeff1)
00334     else
00335         fn_val=fn_val*X(numcoeff1)
00336     end if
00337     fn_val=fn_val*vrais_tot
00338     funvls=fn_val
00339 
00340     if (adaptive) then
00341         Xprime2(1:ndim)=matmul(scaleinvsauv(1:ndim,1:ndim,numpat1),&
00342             X(1:ndim))
00343         Xprime(1:ndim)=startsauv(1:ndim,numpat1) + Xprime2(1:ndim)
00344         funvls=funvls*exp(-(sum(startsauv(1:ndim,numpat1)**2) )/2)/&
00345             detersauv(numpat1)
00346         funvls=funvls*exp(-(dot_product(Xprime2(1:ndim),&
00347             startsauv(1:ndim,numpat1))))
00348         funvls=funvls*&
00349             exp((5.d-1)*dot_product(X(1:ndim),&
00350             matmul(scaleinv2sauv(1:ndim,1:ndim,numpat1),X(1:ndim))))
00351         !if (abs(funvls) <1.d-100) funvls=1.d-100
00352     end if
00353     funvls=funvls
00354     !if (abs(funvls) <1.d-100) funvls=1.d-100
00355     if (funvls.NE.funvls) funvls=1.d-100
00356     return
00357 end subroutine  FUNSUB_Random
00358 
00359 
00360 !------------------------------------------------------------------------------
00361 !
00362 ! VERSION : 1.0
00363 !
00364 ! MODULE: scores.f90 function FUNSUB_Error
00388 subroutine FUNSUB_Error(NDIM2,X,NF2,FUNVLS)
00389     use WorkingSharedValues
00390     use Constante
00391 
00392     implicit none
00393 
00394     integer,intent(in)::ndim2,nf2
00395     double precision,dimension(ndim2),intent(in)::X
00396     double precision,intent(out) :: funvls
00397 
00398     double precision ::vrais_tot,funvls4,funvls3,funvls2, 
00399         funvls1,fn_val,phid
00400     double precision, dimension(tdef2,npmcomp)::simul,deriv
00401     double precision,dimension(ndim2)::Xprime,Xprime2
00402     integer :: i,j,k
00403     integer,dimension(npmcomp)::compteur
00404 
00405     simul=0.d0
00406     systeme=1
00407 
00408     !
00409     ! *** Call to the ODE solver to have sensitivity equation systems trajectories
00410     !
00411     if (adaptive) then
00412         Xprime(1:ndim)=startsauv(1:ndim,numpat1) + &
00413             matmul(scaleinvsauv(1:ndim,1:ndim,numpat1),X(1:ndim))
00414         call solution(xprime,simul(:,1:npmcomp),deriv(:,1:npmcomp))
00415     else
00416         call solution(x,simul(:,1:npmcomp),deriv(:,1:npmcomp))
00417     end if
00418     deriv=0.d0
00419 
00420     fn_val=0.d0
00421     funvls=0.d0
00422     funvls1=0.d0
00423     vrais_tot=1.d0
00424     compteur=0
00425     do k=1,npmcomp
00426         compteur(k)=0
00427         do j=1,nbobs(numpat1,k)
00428              if (censor(numpat1,j,k).NE.1) then
00429 
00430                 if(debugFUNSUB)write(0,*)schedule(numpat1,j,k),donnees(numpat1,j,k),simul(schedule(numpat1,j,k)+1,k)
00431                 funvls1=donnees(numpat1,j,k)-simul(schedule(numpat1,j,k)+1,k)
00432                 funvls4=donnees(numpat1,j,k)-simul(schedule(numpat1,j,k)+1,k)
00433                 funvls4= (funvls4 /b1(npm-npmcomp+k))**2
00434                 funvls1 = (funvls1**(2.d0)) /(abs(b1(npm-npmcomp+k))**3)
00435                 compteur(k)=compteur(k)+1
00436                 funvls4=-funvls4
00437                 funvls4=funvls4/2
00438                 funvls4=exp(funvls4)/abs(b1(npm-npmcomp+k))
00439                 funvls4=funvls4/sqrt(2.d0*pigrec)
00440             else
00441                 if(debugFUNSUB)write(0,*)"censored",schedule(numpat1,j,k),donnees(numpat1,j,k),simul(schedule(numpat1,j,k)+1,k)
00442                 funvls3=donnees(numpat1,j,k)-simul(schedule(numpat1,j,k)+1,k)
00443                 funvls1=-funvls3/(b1(npm-npmcomp+k)**2)
00444                 funvls3= funvls3/DABS(b1(npm-npmcomp+k))
00445                 funvls1=funvls1/phid(funvls3)
00446                 funvls1=((1.d0)/(DSQRT(2*pigrec)))*funvls1
00447                 funvls3=-(5.d-1)*(funvls3**2)
00448                 funvls1=funvls1*DEXP(funvls3)
00449                 funvls4=phid((donnees(numpat1,j,k)-simul(schedule(numpat1,j,k)+1,k))/DABS(b1(npm-npmcomp+k)))
00450                 if (funvls1 .ne. funvls1) funvls1=0.d0
00451                 if (funvls3 .ne. funvls3) funvls3=0.d0
00452                 if (funvls4 .ne. funvls4) funvls4=0.d0
00453             end if
00454             if (k .ne. numcoeff1) funvls1=0.d0
00455             fn_val=fn_val + funvls1
00456             vrais_tot=vrais_tot*funvls4
00457         end do
00458     end do
00459 
00460     fn_val=fn_val - &
00461         (compteur(numcoeff1))/abs(b1(npm-npmcomp+numcoeff1))
00462     fn_val=fn_val*vrais_tot
00463     funvls=fn_val
00464 
00465     if (adaptive) then
00466         Xprime2(1:ndim)=matmul(scaleinvsauv(1:ndim,1:ndim,numpat1),&
00467             X(1:ndim))
00468         Xprime(1:ndim)=startsauv(1:ndim,numpat1) + Xprime2(1:ndim)
00469         funvls=funvls*exp(-(sum(startsauv(1:ndim,numpat1)**2) )/2)/&
00470             detersauv(numpat1)
00471         funvls=funvls*exp(-(dot_product(Xprime2(1:ndim),&
00472             startsauv(1:ndim,numpat1))))
00473         funvls=funvls*&
00474             exp((5.d-1)*dot_product(X(1:ndim),&
00475             matmul(scaleinv2sauv(1:ndim,1:ndim,numpat1),X(1:ndim))))
00476         !if (abs(funvls)<1.d-100) funvls=1.d-100
00477     end if
00478     funvls=funvls
00479     if (funvls.NE.funvls) funvls=1.d-100
00480    ! if (abs(funvls)<1.d-100) funvls=1.d-100
00481     return
00482 end subroutine FUNSUB_Error
00483