marquardt.f90

Go to the documentation of this file.

!------------------------------------------------------------------------------
! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
!                effects based on Ordinary Differential equations
!------------------------------------------------------------------------------
!
! VERSION : 1.0
!
! MODULE: marquardt.f90 function marquardt
subroutine marquardt(b,m,ni,v,rl,istop,maxiter,dd)

    use ToleranceThreshold
    use VectorLength
    USE m_vpropre
    use mpimod

    IMPLICIT NONE
    !   variables globales
    integer,intent(in) :: m,maxiter
    integer,intent(inout)::ni,istop
    double precision,dimension(m*(m+3)/2),intent(out)::v
    double precision,intent(out)::rl
    double precision,dimension(m),intent(inout)::b
    double precision,intent(out)::dd
    double precision::ca,cb

    !   variables locales
    integer::nql,ii,nfmax,idpos,ncount,id,jd,j,i,ij,integer,ier
    double precision,dimension(m*(m+3)/2)::fu
    double precision,dimension(m)::delta,b1,bh
    double precision::da,dm,ga,tr
    double precision::GHG,det,step,eps,vw,fi,maxt, 
        z,rl1,th,ep,epsa,epsb,epsd
    double precision,external::funcpaRandomEffect


    epsa=marquardt_ca
    epsb=marquardt_cb
    epsd=marquardt_rdm
    eps=gonfle_eps
    da=gonfle_da
    dm=gonfle_dm
    ep=dsinv_ep
    th=searpas_th

    id=0
    jd=0
    z=0.d0


    nfmax=m*(m+1)/2
    ca=epsa+1.d0
    cb=epsb+1.d0
    rl1=-1.d+10
    ni=0
    istop=0

    nql=1


    Main:Do
        if(debugMARC)write(0,*)
        if(debugMARC)write(0,*)'***     Iteration',ni,'      ***'
        if(debugMARC) write(0,*)'B=',(b(j),j=1,m)
        call derivRANDOMEFFECT(b,m,v,rl)
        if(debugMARC)write(0,*)'Log-likelihood=',rl


        rl1=rl

        dd = 0.d0
        fu=0.D0
        do i=1,m
            do j=i,m
                ij=(j-1)*j/2+i
                fu(ij)=v(ij)
            end do
        end do

        if(debugMARC)write(0,*)"mat. hessian",fu

        ier=0
        call dsinv(fu,m,ep,ier,det)
        if (ier.eq.-1) then
            dd=epsd+1.d0
        else
            GHG = 0.d0
            do i=1,m
                do j=1,m
                    if(j.ge.i) then
                        ij=(j-1)*j/2+i
                    else
                        ij=(i-1)*i/2+j
                    end if
                    GHG = GHG + v(nfmax+i)*fu(ij)*V(nfmax+j)
                end do
            end do
            dd=GHG/dble(m)
        end if

        if(ca.lt.epsa.and.cb.lt.epsb.and.dd.lt.epsd) exit main

        tr=0.d0
        do i=1,m
            ii=i*(i+1)/2
            tr=tr+dabs(v(ii))
        end do
        tr=tr/dble(m)

        ncount=0
        ga=gonfle_ga
400 continue
    do i=1,nfmax+m
        fu(i)=v(i)
    end do
    do i=1,m
        ii=i*(i+1)/2
        if (v(ii).ne.0) then
            fu(ii)=v(ii)+da*((1.d0-ga)*dabs(v(ii))+ga*tr)
        else
            fu(ii)=da*ga*tr
        endif
    end do
    call dchole(fu,m,nql,idpos)
    if (idpos.ne.0) then
        ncount=ncount+1
        if (ncount.le.3.or.ga.ge.1.d0) then
            da=da*dm
        else
            ga=ga*dm
            if (ga.gt.1.d0) ga=1.d0
        endif
        if (ncount .LE. 100000)goto 400
    else
        vw=1.0D0
        do i=1,m
            delta(i)=fu(nfmax+i)
            b1(i)=b(i)+delta(i)
        end do
        rl=funcpaRandomEffect(b1,m,id,z,jd,z)
        if (rl1.lt.rl) then
            if(da.lt.eps) then
                da=eps
            else
                da=da/(dm+2.d0)
            endif
            goto 800
        endif
    endif
    call dmaxt(maxt,delta,m)
    if(maxt.eq.0.D0) then
        vw=th
    else
        vw=th/maxt
    endif
    step=dlog(1.5d0)
    call searpas(vw,step,b,bh,m,delta,fi,funcpaRandomEffect)
    rl=-fi
    if(rl.eq.-1.D9) then
        istop=4
        goto 110
    end if

    do i=1,m
        delta(i)=vw*delta(i)
    end do
    da=(dm-3.d0)*da

800 cb=dabs(rl1-rl)
    ca=0.d0
    do i=1,m
        ca=ca+delta(i)*delta(i)
    end do
    if(debugMARC)write(0,*) 'ca =',ca,' cb =',cb,' dd =',dd
    if(debugMARC)write(0,*)"delta= ",delta,rl
    if(debugMARC)write(0,*)"vw= ",vw
    if(debugMARC.and.stepbystep)pause
    do i=1,m
        b(i)=b(i)+delta(i)
    end do

    ni=ni+1
    if (ni.ge.maxiter) then
        istop=2
        if(debugMARC)write(0,*) 'maximum number of iteration reached'
        goto 110
    end if
End do Main
v=0.D0
v(1:m*(m+1)/2)=fu(1:m*(m+1)/2)
istop=1

110 continue
    return
end subroutine marquardt