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

ode.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: ode.f90 function solution
00034 subroutine solution(x,simul,deriv)
00035     use WorkingSharedValues
00036     implicit none
00037 
00038     double precision, dimension(ndim),intent(in)::x
00039     double precision,dimension(tdef2,npmcomp),intent(out)::simul,deriv
00040 
00041     EXTERNAL  FEX, JEX,FEXcl,JEXcl,FEXV0,JEXV0
00042     EXTERNAL  FEXka,JEXka
00043 
00044     INTEGER  IOPT, IOUT, ISTATE, ITASK, ITOL, IWORK(50), LIW, LRW, MF, NEQ
00045     DOUBLE PRECISION,dimension(2*nbcomp)::  ATOL,Y,YDOT
00046     DOUBLE PRECISION :: RTOL, RWORK(10000), T, TOUT
00047     integer::i,ifail
00048 
00049 ! *** Parameter initialisation
00050     simul=0.d0
00051     deriv=0.d0
00052     t=0.d0
00053     Y=0.d0
00054 
00055     b0=b1
00056     do i=1,ndim
00057        b0(i)=b0(i)+X(i)*b1(npmbio+npmexpl+i)
00058    end do
00059 
00060     ITASK = 1
00061     ISTATE = 1
00062     IOPT = 0
00063     RWORK=10000
00064     LRW =10000
00065     LIW =100000
00066     MF = 21
00067     ITOL=1
00068     atol=ODESolveurATOL
00069     rtol=0.d0
00070     call transfoFixedInTime()
00071 !
00072 ! *** ODE Systems trajectories for b0 : time by time resolution
00073 !
00074     if (systeme(1) .eq. 1) then
00075         ! *** Initial values
00076         NEQ = nbcomp
00077         call initialPoints(t,Y,NEQ)
00078         Call SystemResolution(Y,1,simul,deriv,ifail)
00079         ! *** Time by time
00080         DO IOUT = 1,tailletime(numpat1)-1!tps3
00081           t=dble(listtime(numpat1,IOUT))
00082           tout=dble(listtime(numpat1,IOUT+1))
00083 
00084             CALL DLSODE (FEX, NEQ, Y,t, tout, ITOL, RTOL, ATOL, ITASK, &
00085                 ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JEX, MF)
00086             Call SystemResolution(Y,int(tout)+1,simul,deriv,ifail)
00087             IF ((ifail.EQ.1).OR.(ISTATE .LT. 0) ) exit
00088         END DO
00089 !
00090 ! *** ODE Systems Sensitivity equation trajectories for b0 : time by time resolution
00091 !
00092     else
00093         ! *** Initial values
00094         NEQ=2*nbcomp
00095         call initialPoints(t,Y,NEQ)
00096         call SystemResolution(Y,1,simul,deriv,ifail)
00097         ! *** Time by time, parameters must be given in the correct order.
00098         DO IOUT = 1,tailletime(numpat1)-1!tps3
00099           t=dble(listtime(numpat1,IOUT))
00100           tout=dble(listtime(numpat1,IOUT+1))
00101 
00102             ! *** cl is the first parameter
00103             if(systeme(2) .eq. 1 ) then
00104                !.or.(systeme(2) .eq. (npmbio+1) ) *** If the first explanatory variable is applied to cl
00105                 CALL DLSODE (FEXcl, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, &
00106                     ITASK,ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JEXcl, MF)
00107             else
00108                 ! *** V0 is the second parameter
00109                 if( systeme(2) .eq. 2 ) then
00110                     CALL DLSODE (FEXV0, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, &
00111                         ITASK,ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JEXV0, MF)
00112                 else
00113                     ! *** ka is the third parameter
00114                     if( systeme(2) .eq. 3 ) then
00115                         CALL DLSODE (FEXka, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, &
00116                             ITASK,ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JEXka, MF)
00117                     end if
00118                 end if
00119             end if
00120             call SystemResolution(Y,int(tout)+1,simul,deriv,ifail)
00121             IF ((ifail.EQ.1).OR.(ISTATE .LT. 0) ) exit
00122         end do
00123     end if
00124     return
00125 END subroutine solution
00126 
00127 !------------------------------------------------------------------------------
00128 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00129 !                effects based on Ordinary Differential equations
00130 !------------------------------------------------------------------------------
00131 !
00132 ! VERSION : 1.0
00133 !
00134 ! MODULE: ode.f90 function FEX
00158 SUBROUTINE  FEX (NEQ, T, Y, YDOT)
00159     use WorkingSharedValues
00160     implicit none
00161     INTEGER, intent(in)::  NEQ
00162     DOUBLE PRECISION,intent(in)::  T, Y(NEQ)
00163     DOUBLE PRECISION,intent(out):: YDOT(NEQ)
00164     call transfo(t,Y,NEQ)
00165     YDOT(1)=-ka*Y(1)
00166     YDOT(2)=(ka*Y(1)-cl*Y(2)*V0)/V0
00167     return
00168 END
00169 
00170 !------------------------------------------------------------------------------
00171 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00172 !                effects based on Ordinary Differential equations
00173 !------------------------------------------------------------------------------
00174 !
00175 ! VERSION : 1.0
00176 !
00177 ! MODULE: ode.f90 function JEX
00201 SUBROUTINE  JEX (NEQ, T, Y, ML, MU, PD, NRPD)
00202     use WorkingSharedValues
00203     implicit none
00204     INTEGER, intent(in)::  NEQ, ML, MU, NRPD
00205     DOUBLE PRECISION, intent(in)::  T, Y(NEQ)
00206     DOUBLE PRECISION, intent(out):: PD(nrpd,NEQ)
00207     call transfo(t,Y,NEQ)
00208     PD(1,1)=-ka
00209     PD(2,1)=ka/V0
00210     PD(2,2)=-cl
00211     return
00212 end
00213 
00214 !
00215 ! *** ODE System Sensitivity equations definition
00216 !
00217 
00218 !------------------------------------------------------------------------------
00219 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00220 !                effects based on Ordinary Differential equations
00221 !------------------------------------------------------------------------------
00222 !
00223 ! VERSION : 1.0
00224 !
00225 ! MODULE: ode.f90 function FEXcl
00249 SUBROUTINE  FEXcl(NEQ, T, Y, YDOT)
00250     use WorkingSharedValues
00251     implicit none
00252     INTEGER, intent(in)::  NEQ
00253     DOUBLE PRECISION,intent(in)::  T, Y(2*NEQ)
00254     DOUBLE PRECISION,intent(out):: YDOT(2*NEQ)
00255     call transfo(t,Y,NEQ)
00256     YDOT(1)=-ka*Y(1)
00257     YDOT(2)=(ka*Y(1)-cl*Y(2)*V0)/V0
00258     YDOT(3)=-ka*Y(3)
00259     YDOT(4)=(ka*Y(3)-dcl*Y(2)*V0-cl*Y(4)*V0)/V0
00260     return
00261 end
00262 
00263 !------------------------------------------------------------------------------
00264 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00265 !                effects based on Ordinary Differential equations
00266 !------------------------------------------------------------------------------
00267 !
00268 ! VERSION : 1.0
00269 !
00270 ! MODULE: ode.f90 function JEXcl
00294 SUBROUTINE  JEXcl (NEQ, T, Y, ML, MU, PD, NRPD)
00295     use WorkingSharedValues
00296     implicit none
00297     INTEGER, intent(in)::  NEQ, ML, MU, NRPD
00298     DOUBLE PRECISION, intent(in)::  T, Y(2*NEQ)
00299     DOUBLE PRECISION, intent(out):: PD(nrpd,2*NEQ)
00300     call transfo(t,Y,NEQ)
00301     PD(1,1)=-ka
00302     PD(2,1)=ka/V0
00303     PD(2,2)=-cl
00304     PD(3,3)=-ka
00305     PD(4,2)=-dcl
00306     PD(4,3)=ka/V0
00307     PD(4,4)=-cl
00308     return
00309 end
00310 
00311 !------------------------------------------------------------------------------
00312 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00313 !                effects based on Ordinary Differential equations
00314 !------------------------------------------------------------------------------
00315 !
00316 ! VERSION : 1.0
00317 !
00318 ! MODULE: ode.f90 function FEXV0
00342 SUBROUTINE  FEXV0 (NEQ, T, Y, YDOT)
00343     use WorkingSharedValues
00344     implicit none
00345     INTEGER, intent(in)::  NEQ
00346     DOUBLE PRECISION,intent(in)::  T, Y(2*NEQ)
00347     DOUBLE PRECISION,intent(out):: YDOT(2*NEQ)
00348     call transfo(t,Y,NEQ)
00349     YDOT(1)=-ka*Y(1)
00350     YDOT(2)=(ka*Y(1)-cl*Y(2)*V0)/V0
00351     YDOT(3)=-ka*Y(3)
00352     YDOT(4)=(ka*Y(3)-cl*Y(4)*V0-cl*Y(2)*dV0)/V0-(ka*Y(1)-cl*Y(2)*V0)*dV0/V0**2
00353 
00354     return
00355 end
00356 
00357 !------------------------------------------------------------------------------
00358 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00359 !                effects based on Ordinary Differential equations
00360 !------------------------------------------------------------------------------
00361 !
00362 ! VERSION : 1.0
00363 !
00364 ! MODULE: ode.f90 function JEXV0
00388 SUBROUTINE  JEXV0 (NEQ, T, Y, ML, MU, PD, NRPD)
00389     use WorkingSharedValues
00390     implicit none
00391     INTEGER, intent(in)::  NEQ, ML, MU, NRPD
00392     DOUBLE PRECISION, intent(in)::  T, Y(2*NEQ)
00393     DOUBLE PRECISION, intent(out):: PD(nrpd,2*NEQ)
00394     call transfo(t,Y,NEQ)
00395     PD(1,1)=-ka
00396     PD(2,1)=ka/V0
00397     PD(2,2)=-cl
00398     PD(3,3)=-ka
00399     PD(4,1)=-ka*dV0/V0**2
00400     PD(4,3)=ka/V0
00401     PD(4,4)=-cl
00402     return
00403 end
00404 
00405 !------------------------------------------------------------------------------
00406 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00407 !                effects based on Ordinary Differential equations
00408 !------------------------------------------------------------------------------
00409 !
00410 ! VERSION : 1.0
00411 !
00412 ! MODULE: ode.f90 function FEXka
00436 SUBROUTINE  FEXka (NEQ, T, Y, YDOT)
00437     use WorkingSharedValues
00438     implicit none
00439     INTEGER, intent(in)::  NEQ
00440     DOUBLE PRECISION,intent(in)::  T, Y(2*NEQ)
00441     DOUBLE PRECISION,intent(out):: YDOT(2*NEQ)
00442     call transfo(t,Y,NEQ)
00443     YDOT(1)=-ka*Y(1)
00444     YDOT(2)=(ka*Y(1)-cl*Y(2)*V0)/V0
00445     YDOT(3)=-dka*Y(1)-ka*Y(3)
00446     YDOT(4)=(dka*Y(1)+ka*Y(3)-cl*Y(4)*V0)/V0
00447     return
00448 end
00449 
00450 !------------------------------------------------------------------------------
00451 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00452 !                effects based on Ordinary Differential equations
00453 !------------------------------------------------------------------------------
00454 !
00455 ! VERSION : 1.0
00456 !
00457 ! MODULE: ode.f90 function JEXka
00481 SUBROUTINE  JEXka (NEQ, T, Y, ML, MU, PD, NRPD)
00482     use WorkingSharedValues
00483     implicit none
00484     INTEGER, intent(in)::  NEQ, ML, MU, NRPD
00485     DOUBLE PRECISION, intent(in)::  T, Y(2*NEQ)
00486     DOUBLE PRECISION, intent(out):: PD(nrpd,2*NEQ)
00487     call transfo(t,Y,NEQ)
00488     PD(1,1)=-ka
00489     PD(2,1)=ka/V0
00490     PD(2,2)=-cl
00491     PD(3,1)=-dka
00492     PD(3,3)=-ka
00493     PD(4,1)=dka/V0
00494     PD(4,3)=ka/V0
00495     PD(4,4)=-cl
00496     return
00497 end
00498 
00499 
00500 !------------------------------------------------------------------------------
00501 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00502 !                effects based on Ordinary Differential equations
00503 !------------------------------------------------------------------------------
00504 !
00505 ! VERSION : 1.0
00506 !
00507 ! MODULE: ode.f90 function initialPoints
00530 SUBROUTINE initialPoints(t,Y,NEQ)
00531    use WorkingSharedValues
00532    implicit none
00533    DOUBLE PRECISION, intent(in)::t
00534    INTEGER, intent(in)::  NEQ
00535    DOUBLE PRECISION,intent(out):: Y(NEQ)
00536    call transfo(t,Y,NEQ)
00537    Y(1)=dose0
00538    Y(2)=0.0D0
00539    if(systeme(1) .ne. 1) then
00540        ! *** derivative of Y with respect to the first parameter
00541        if (systeme(2) .eq. 1 ) then
00542                 Y(3)=0.0D0
00543                 Y(4)=0.0D0
00544             else
00545                 ! *** derivative of Y with respect to the second parameter
00546                 if (systeme(2) .eq. 2) then
00547                     Y(3)=0.0D0
00548                     Y(4)=0.0D0
00549                 else
00550                     ! *** derivative of Y with respect to the third parameter
00551                     if (systeme(2) .eq. 3 )then
00552                         Y(3)=0.0D0
00553                         Y(4)=0.0D0
00554                     end if
00555                 end if
00556             end if
00557         end if
00558    RETURN
00559 END
00560 
00561 
00562 
00563 
00564 
00565