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

inAndOutInternal.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 : 2.0
00007 !
00008 ! MODULE: inAndOutInternal.f90 function OutputFileCreation
00036 subroutine OutputFileCreation(b,v,ni,istop,conv,Vrais,VraisNP,LCV_a,tpscpu )
00037     use WorkingSharedValues
00038     use mpimod
00039     implicit none
00040     double precision, dimension(npm),intent(in) ::b
00041     integer,intent(in) ::ni
00042     integer,intent(in) ::istop
00043     double precision, dimension(3),intent(in) ::conv
00044     double precision,intent(in) ::Vrais,VraisNP,LCV_a,tpscpu
00045     double precision,dimension(npm*(npm+3)/2),intent(in)::V
00046 
00047 
00048     integer :: i
00049     logical :: logic
00050 
00051     open(543,file="./"//trim(adjustl(OutputFolder))//"/Convergence_LATEX.txt")
00052 
00053     if (numproc.EQ.0) write(543,*)"\documentclass[12pt,a4paper]{article}"
00054     if (numproc.EQ.0) write(543,*)"\usepackage[english]{babel}"
00055     if (numproc.EQ.0) write(543,*)"\usepackage[dvips,dvipsnames,usenames]{color}"
00056     if (numproc.EQ.0) write(543,*)"\usepackage[margin=2.5cm, head=1.27cm]{geometry}"
00057     if (numproc.EQ.0) write(543,*)"\","begin{document}
00058 
00059 "    if (numproc.EQ.0) write(543,*)","section*{Model Specification :}"
00060     if (numproc.EQ.0) write(543,*)"Data used : $",biomarkersFile,"$ \","
00061     if (numproc.EQ.0) write(543,*)"Biological parameters : ",npmbio," \","
00062     if (numproc.EQ.0) write(543,*)"Explanatory variables : ",npmexpl," \","
00063     if (numproc.EQ.0) write(543,*)"Random Effect number  : ",ndim," \","
00064     if (numproc.EQ.0) write(543,*)"Observed compartments : ",npmcomp," \","
00065 
00066     if (numproc.EQ.0) write(543,*)"\","section*{Algorithm Specification :}
00067 
00068 
00069 "    logic=.true.    if(penalisationBiologique.EQ.0)logic=.false.    if (numproc.EQ.0) write(543,*)"Biological Penalisation : ",logic," ,
00070 
00071 
00072 "\"    logic=.true.    if(penalisationAll.EQ.0)logic=.false.    if (numproc.EQ.0) write(543,*)"All parameters Penalisation : ",logic," ,
00073 "\"    if (numproc.EQ.0) write(543,*)"Treshold for Parameter variation  : ",controleCA," ,
00074 "\"    if (numproc.EQ.0) write(543,*)"Treshold for Likelihood variation : ",controleCB," ,
00075 "\"    if (numproc.EQ.0) write(543,*)"Treshold for RDM : ",controleRDM," ,
00076 
00077 
00078 "\"    if (numproc.EQ.0) write(543,*)","section*{Results :}"
00079     if (numproc.EQ.0) write(543,*)"\","begin{tabular}{rrrrrrrrrr}
00080 "    if (numproc.EQ.0) write(543,*)"& &  \multicolumn{ 2}{c}{PRIOR} & \multicolumn{ 2}{c}{POSTERIOR} 
00081 ", &        " & \multicolumn{ 2}{c}{POSTERIOR} & wald"," ,
00082 "\"    if (numproc.EQ.0) write(543,*)"&Start &  \multicolumn{ 2}{c}{(log-scale)} & \multicolumn{ 2}{c}{(log-scale)} &  &
00083                                                 \multicolumn{ 2}{c}{(natural-scale)} &"    ," ,
00084 "\"    if (numproc.EQ.0) write(543,*)"&Point &  mean &  sd  &  mean &     sd &   mean &   sd &"," ,
00085 "\"    if (numproc.EQ.0) write(543,*)"\hline
00086 
00087 
00088 
00089 
00090 "    if(penalisationBiologique.EQ.1) then        do i=1,npmbio            if (numproc.EQ.0) write(543,'(3(A),7(F20.6,A),A)')"$,nomparams(i),"$&",binit(i),"&",esp_prior(i), &
00091                 "&",std_prior(i),"&",b(i),"&", &
00092                 sqrt(v(i*(i+1)/2)),"&",exp(b(i)),"&",exp(b(i))*sqrt(v(i*(i+1)/2)), " & \","
00093         end do
00094     else
00095         do i=1,npmbio
00096             if (numproc.EQ.0) write(543,'(3(A),5(F20.6,A),A)')"$\",nomparams(i),"$&",binit(i),"& - & - &",b(i),"&
00097 ", &                sqrt(v(i*(i+1)/2)),"&",exp(b(i)),"&",exp(b(i))*sqrt(v(i*(i+1)/2))," & ,
00098 
00099 
00100 
00101 
00102 
00103 "\"        end do    end if    if(penalisationAll.EQ.1) then        do i=npmbio+1,npmbio+npmexpl            if (numproc.EQ.0) write(543,'(3(A),6(F20.6,A),A)')"$,nomparams(i),"$&",binit(i),"&",esp_prior(i),"&", &
00104                 std_prior(i),"&",b(i),"&", &
00105                 sqrt(v(i*(i+1)/2)),"& - & - &",b(i)/sqrt(v(i*(i+1)/2))," \","
00106         end do
00107     else
00108         do i=npmbio+1,npmbio+npmexpl
00109             if (numproc.EQ.0) write(543,'(3(A),4(F20.6,A),A)')"$\",nomparams(i),"$&",binit(i),"& - & - &",b(i),"&
00110 ", &                sqrt(v(i*(i+1)/2)),"& - & - &",b(i)/sqrt(v(i*(i+1)/2))," ,
00111 
00112 
00113 
00114 
00115 
00116 "\"        end do    end if    if(penalisationAll.EQ.1) then        do i=npmbio+npmexpl+1,npmbio+npmexpl+ndim            if (numproc.EQ.0) write(543,'(3(A),5(F20.6,A),A)')"$,nomparams(i),"$ &",binit(i),"&",esp_prior(i),"& - &", &
00117                 b(i),"&", &
00118                 sqrt(v(i*(i+1)/2)),"& - & - &",dabs(b(i))/sqrt(v(i*(i+1)/2))," \","
00119         end do
00120     else
00121         do i=npmbio+npmexpl+1,npmbio+npmexpl+ndim
00122             if (numproc.EQ.0) write(543,'(3(A),4(F20.6,A),A)')"$\",nomparams(i),"$ &",binit(i),"& - & - &",b(i),"&
00123 ", &                sqrt(v(i*(i+1)/2)),"& - & - &",dabs(b(i))/sqrt(v(i*(i+1)/2))," ,
00124 
00125 
00126 
00127 
00128 "\"        end do    end if    do i=npmbio+npmexpl+ndim+1,npm        if (numproc.EQ.0) write(543,'(3(A),4(F20.6,A),A)')"$,nomparams(i),"$ &",binit(i),"& - & - &", &
00129             b(i),"&", &
00130             sqrt(v(i*(i+1)/2)),"& - & - &",dabs(b(i))/sqrt(v(i*(i+1)/2))," \","
00131     end do
00132 
00133 
00134     if (numproc.EQ.0) write(543,*)"\end{tabular}"," \","
00135     if (numproc.EQ.0) write(543,*)" \","
00136     if (numproc.EQ.0) write(543,*)" \","
00137     if (numproc.EQ.0) write(543,*)"Convergence status (1: OK; 2:Max. iter. reached; ", &
00138         " 3:non-inv. Hessian ; 4 got stuck) : ", istop," \","
00139     if (numproc.EQ.0) write(543,*)"Time : ", tpscpu ," sec. on",nbprocs+1," cores \","
00140     if (numproc.EQ.0) write(543,*)"Number of iterations : ",ni," \","
00141     if (numproc.EQ.0) write(543,*)"Parameter variation : ",conv(1)," \","
00142     if (numproc.EQ.0) write(543,*)"Likelihood variation : ",conv(2)," \","
00143     if (numproc.EQ.0) write(543,*)"RDM : ",conv(3)," \","
00144     if (numproc.EQ.0) write(543,'(A,F20.6,A,A)')"Penalized Log-Likelihood : ",LLPenalisee," \","
00145     if (numproc.EQ.0) write(543,'(A,F20.6,A,A)')"NON-Penalized Log-Likelihood : ",LLnonPenalisee," \","
00146     if (numproc.EQ.0) write(543,'(A,F20.6,A,A)')"LCVa : ",LCV_a," \","
00147 
00148     if (numproc.EQ.0) write(543,*)"\end{document}"
00149     flush(543)
00150     close(543)
00151     return
00152 end subroutine OutputFileCreation
00153 
00154 
00155 !------------------------------------------------------------------------------
00156 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00157 !                effects based on Ordinary Differential equations
00158 !------------------------------------------------------------------------------
00159 !
00160 ! VERSION : 2.0
00161 !
00162 ! MODULE: inAndOutInternal.f90 function InputReading
00181 subroutine InputReading()
00182     use WorkingSharedValues
00183     use mpimod
00184     implicit none
00185     integer :: i,j
00186     integer :: inputSelection
00187     open(221,file="Input.txt")
00188     read(221,*)
00189     read(221,*)
00190     ! Name and initial value reading (-99 : departure value is the prior mean)
00191     do i=1,npmbio
00192         read(221,*)nomparams(i),binit(i)
00193     end do
00194     do i=npmbio+1,npm
00195         read(221,*)nomparams(i),binit(i)
00196     end do
00197     read(221,*)
00198     ! I/O files Reading
00199     read(221,*)
00200     read(221,*)OutputFolder
00201     read(221,*)
00202     read(221,*)biomarkersFile
00203     read(221,*)
00204     read(221,*)treatmentFile
00205     read(221,*)
00206     ! Algorithm parameter Reading
00207     read(221,*)
00208     read(221,*)controleCA
00209     read(221,*)
00210     read(221,*)controleCB
00211     read(221,*)
00212     read(221,*)controleRDM
00213     read(221,*)
00214     read(221,*)maxiteration
00215     read(221,*)
00216     read(221,*)marqONLY
00217     read(221,*)
00218     read(221,*)SwitchMarquartd
00219     read(221,*)
00220     read(221,*)estimationWanted
00221     read(221,*)
00222     ! Priors Reading
00223     read(221,*)
00224     read(221,*)penalisationBiologique
00225     read(221,*)
00226     read(221,*)penalisationAll
00227     read(221,*)
00228     if (penalisationBiologique.EQ.1) then
00229         do i=1,npmbio
00230             read(221,*)esp_prior(i),std_prior(i)
00231         end do
00232     end if
00233     if (penalisationAll.EQ.1) then
00234         do i=npmbio+1,npmbio+npmexpl
00235             read(221,*)esp_prior(i),std_prior(i)
00236         end do
00237         do i=npmbio+npmexpl+1,npmbio+npmexpl+ndim
00238             read(221,*)esp_prior(i)
00239         end do
00240     end if
00241     read(221,*)
00242 
00243 
00244     do i=1,npm
00245         if(binit(i).EQ.-99)binit(i)=esp_prior(i)
00246     end do
00247 
00248     close(221)
00249 
00250     if((numproc.EQ.0).and.(LOG_PRINT))write(0,*)"! *** System size informations"
00251     if((numproc.EQ.0).and.(LOG_PRINT))write(0,*)"Data input file name :", biomarkersFile
00252     if((numproc.EQ.0).and.(LOG_PRINT))write(0,*)"Number of biological parameters, ","npmbio=",npmbio
00253     if((numproc.EQ.0).and.(LOG_PRINT))write(0,*)"Number of explanatory variables, ","npmexpl=",npmexpl
00254     if((numproc.EQ.0).and.(LOG_PRINT))write(0,*)"Number of random effects, ","ndim=",ndim
00255     if((numproc.EQ.0).and.(LOG_PRINT))write(0,*)"Number of observed compartments, ","npmcomp=",npmcomp
00256     if((numproc.EQ.0).and.(LOG_PRINT))write(0,*)"Total Number of compartments, ","nbcomp=",nbcomp
00257 
00258     return
00259 end subroutine InputReading
00260 
00261 !------------------------------------------------------------------------------
00262 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00263 !                effects based on Ordinary Differential equations
00264 !------------------------------------------------------------------------------
00265 !
00266 ! VERSION : 2.0
00267 !
00268 ! MODULE: inAndOutInternal.f90 function TraceR
00287 subroutine TraceR(b)
00288     use WorkingSharedValues
00289     use mpimod
00290    implicit none
00291    integer :: i
00292    double precision,dimension(npm),intent(in)::b
00293 
00294    open(555, file="./"//trim(adjustl(OutputFolder))//"/graphics/programR.txt")
00295    write(555,*)"DATA<-read.table('",trim(adjustl(biomarkersFile)),"',header=TRUE) "
00296    write(555,*)"id<-unique(DATA[,1])"
00297    write(555,*)"FIT<-read.table('./",trim(adjustl(OutputFolder)),"/PEBTrajectories.txt',header=FALSE)"
00298    write(555,*)"tmax<-",int(tdef2-15)
00299     write(555,*)"npmcomp<-",npmcomp
00300    write(555,*) "for (i in 1:",nbpatient,"){"
00301    write(555,*)"ylimax<-min(tmax,max(DATA[which(DATA[,1]==i),2])*1.10)"
00302    write(555,*)"xlim<-c(0,ylimax)"
00303 
00304    write(555,*)"idpat<-id[i]"
00305    write(555,*)"pdf(paste('./"//trim(adjustl(OutputFolder))//
00306 "/graphics/', &                      paste(paste('pat_',idpat,sep=''),'.pdf',sep=''),sep=''))"
00307    write(555,*) "par(mfrow=c(1,npmcomp))"
00308    write(555,*) "for (j in 1:npmcomp){"
00309    write(555,*)"infborne<-min(c(DATA[which((DATA[,1]==i)&(DATA[,2+j]!=-99)),2+j],FIT[which(FIT[,1]==i),2+npmcomp+j]))"
00310    write(555,*)"maxborne<-max(c(DATA[which((DATA[,1]==i)&(DATA[,2+j]!=-99)),2+j],FIT[which(FIT[,1]==i),2+2*npmcomp+j]))"
00311    write(555,*)"ylim<-c(infborne-0.25*abs(infborne),maxborne+0.25*abs(maxborne))"
00312    write(555,*)"plot(0,0,xlim=xlim,ylim=ylim,xlab='Time',ylab=paste('Biomarker ',j,sep=),",  &
00313                 "font.lab=2, font.axis=2, font=2,col='white')"
00314 
00315    write(555,*)
00316 "lines(FIT[which((FIT[,1]==i)&(FIT[,2]<tmax)),2], &                     FIT[which((FIT[,1]==i)&(FIT[,2]<tmax)),2+npmcomp+j],lwd=1,lty=2)"    
00317 
00318    write(555,*)
00319 "lines(FIT[which((FIT[,1]==i)&(FIT[,2]<tmax)),2], &                     FIT[which((FIT[,1]==i)&(FIT[,2]<tmax)),2+2*npmcomp+j],lwd=1,lty=2)"  
00320 
00321    write(555,*)"polygon(c(FIT[which((FIT[,1]==i)&(FIT[,2]<tmax)),2],", &
00322                   "rev(FIT[which((FIT[,1]==i)&(FIT[,2]<tmax)),2])),", &
00323                  "c(FIT[which((FIT[,1]==i)&(FIT[,2]<tmax)),2+2*npmcomp+j],", &
00324                   "rev(FIT[which((FIT[,1]==i)&(FIT[,2]<tmax)),2+npmcomp+j])),",&
00325                  "col = 'grey', border = NA)"
00326 
00327    write(555,*)
00328 "lines(FIT[which((FIT[,1]==i)&(FIT[,2]<tmax)),2], &                     FIT[which((FIT[,1]==i)&(FIT[,2]<tmax)),2+j],lwd=2)"
00329 
00330 
00331 
00332    write(555,*)"points(DATA[which((DATA[,1]==i)&(DATA[,2+j]!=-99)),2], ", &
00333                       "(DATA[which((DATA[,1]==i)&(DATA[,2+j]!=-99)),2+j]),pch=17,cex=2)"
00334 
00335    write(555,*) "}"
00336    write(555,*)"dev.off()"
00337    write(555,*) "}"
00338    close(555)
00339  end subroutine TraceR
00340 
00341 
00342 
00343  !------------------------------------------------------------------------------
00344 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00345 !                effects based on Ordinary Differential equations
00346 !------------------------------------------------------------------------------
00347 !
00348 ! VERSION : 2.0
00349 !
00350 ! MODULE: inAndOutInternal.f90 function CreateRelanceFile
00369 subroutine CreateRelanceFile(b)
00370     use WorkingSharedValues
00371     use mpimod
00372     implicit none
00373     double precision, dimension(npm),intent(in) ::b
00374     integer :: i,j
00375 
00376     open(221,file="./"//trim(adjustl(OutputFolder))//"/relance.txt")
00377     write(221,*)
00378     write(221,*)
00379     ! Name and initial value reading (-99 : departure value is the prior mean)
00380     do i=1,npmbio
00381         write(221,*)nomparams(i),b(i)
00382     end do
00383     do i=npmbio+1,npm
00384         write(221,*)nomparams(i),b(i)
00385     end do
00386     write(221,*)
00387     ! I/O files Reading
00388     write(221,*)
00389     write(221,*)OutputFolder
00390     write(221,*)
00391     write(221,*)biomarkersFile
00392     write(221,*)
00393     write(221,*)treatmentFile
00394     write(221,*)
00395     ! Algorithm parameter Reading
00396     write(221,*)
00397     write(221,*)controleCA
00398     write(221,*)
00399     write(221,*)controleCB
00400     write(221,*)
00401     write(221,*)controleRDM
00402     write(221,*)
00403     write(221,*)maxiteration
00404     write(221,*)
00405     write(221,*)marqONLY
00406     write(221,*)
00407     write(221,*)SwitchMarquartd
00408     write(221,*)
00409     write(221,*)estimationWanted
00410     write(221,*)
00411     ! Priors Reading
00412     write(221,*)
00413     write(221,*)penalisationBiologique
00414     write(221,*)
00415     write(221,*)penalisationAll
00416     write(221,*)
00417     if (penalisationBiologique.EQ.1) then
00418         do i=1,npmbio
00419             write(221,*)esp_prior(i),std_prior(i)
00420         end do
00421     end if
00422     if (penalisationAll.EQ.1) then
00423         do i=npmbio+1,npmbio+npmexpl
00424             write(221,*)esp_prior(i),std_prior(i)
00425         end do
00426         do i=npmbio+npmexpl+1,npmbio+npmexpl+ndim
00427             write(221,*)esp_prior(i)
00428         end do
00429     end if
00430     write(221,*)
00431     close(221)
00432     return
00433 end subroutine CreateRelanceFile