1! Subroutine for computation of the approximating Gaussian model for non-Gaussian models
2
3subroutine approx(yt, ymiss, timevar, zt, tt, rtv, ht, qt, a1, p1,p1inf, p, n, m, r,&
4theta, u, ytilde, dist, maxiter, tol, rankp, convtol, diff, lik, info, expected, htol)
5
6    implicit none
7
8    integer, intent(in) ::  p,m, r, n,rankp, expected
9    integer, intent(in), dimension(n,p) :: ymiss
10    integer, intent(in), dimension(5) :: timevar
11    integer, intent(in), dimension(p) :: dist
12    integer, intent(inout) :: maxiter,info
13    integer ::  i, k,tvrqr,kk,jt,dt
14    double precision, intent(in) :: tol,convtol
15    double precision, intent(inout) :: htol
16    double precision, intent(in), dimension(n,p) :: u
17    double precision, intent(in), dimension(n,p) :: yt
18    double precision, intent(in), dimension(p,m,(n-1)*timevar(1)+1) :: zt
19    double precision, intent(in), dimension(m,m,(n-1)*timevar(3)+1) :: tt
20    double precision, intent(in), dimension(m,r,(n-1)*timevar(4)+1) :: rtv
21    double precision, intent(in), dimension(r,r,(n-1)*timevar(5)+1) :: qt
22    double precision, intent(in), dimension(m) :: a1
23    double precision, intent(in), dimension(m,m) ::  p1,p1inf
24    double precision, intent(inout), dimension(n,p) :: theta
25    double precision, intent(inout), dimension(n,p) :: ytilde
26    double precision, intent(inout), dimension(p,p,n) :: ht
27    double precision, intent(inout) :: diff
28    double precision, dimension(m,r) :: mr
29    double precision, dimension(m,m,(n-1)*max(timevar(4),timevar(5))+1) :: rqr
30    double precision, intent(inout) :: lik
31    double precision, external :: ddot
32    integer, external :: finitex
33    double precision dev, devold
34    double precision, dimension(n,p) :: thetanew, thetaold
35    double precision, dimension(p,n) :: ft,finf
36    double precision, dimension(m,p,n) :: kt,kinf
37
38    external dgemm, pytheta, pthetafirst, approxloop, pthetarest
39
40    !compute rqr
41    tvrqr = max(timevar(4),timevar(5))
42    do i=1, (n-1)*tvrqr+1
43        call dgemm('n','n',m,r,r,1.0d0,rtv(:,:,(i-1)*timevar(4)+1),m,qt(:,:,(i-1)*timevar(5)+1),r,0.0d0,mr,m)
44        call dgemm('n','t',m,m,r,1.0d0,mr,m,rtv(:,:,(i-1)*timevar(4)+1),m,0.0d0,rqr(:,:,i),m)
45    end do
46
47    ! compute logp(theta) for the first time, no need to compute kt/kinf and ft/finf in successive calls
48    ! in case of totally diffuse initialization term p(theta)=0 for all theta
49    if(rankp .NE. m) then
50        call pthetafirst(theta, timevar, zt, tt, rqr, a1, p1, p1inf, p, m, n, devold, tol,rankp,kt,kinf,ft,finf,dt,jt)
51    end if
52    thetaold = theta
53    devold = -huge(devold)
54    k=0
55    do while(k .LT. maxiter)
56
57        k=k+1
58        ! compute new guess thetanew
59        call approxloop(yt, ymiss, timevar, zt, tt, rtv, ht, qt, rqr, tvrqr, a1, p1,p1inf, p,n,m,r, &
60        theta, thetanew, u, ytilde, dist,tol,rankp,lik, expected)
61        ! and log(p(theta|y))
62        call pytheta(thetanew, dist, u, yt, ymiss, dev, p, n)
63        if(rankp .NE. m) then
64            call pthetarest(thetanew, timevar, zt, tt, a1, p, m, n, dev, kt,kinf,ft,finf,dt,jt)
65        end if
66        !non-finite value in linear predictor or muhat
67        if(finitex(sum(thetanew, MASK = ymiss.EQ.0)).EQ.0 .OR. &
68           finitex(maxval(exp(thetanew), MASK = ymiss.EQ.0)).EQ.0 ) then
69            if(k .GT. 1) then
70                kk = 0
71                do while(finitex(sum(thetanew, MASK = ymiss.EQ.0)).EQ.0 .OR. &
72                         finitex(maxval(exp(thetanew), MASK = ymiss.EQ.0)).EQ.0)
73                    kk = kk + 1
74                    if(kk .GT. maxiter) then
75                        info = 1
76                        return
77                    end if
78                    !backtrack
79                    theta = 0.5d0*(thetaold+theta)
80                    call approxloop(yt, ymiss, timevar, zt, tt, rtv, ht, qt, rqr, tvrqr, a1, p1,p1inf, p,n,m,r, &
81                    theta, thetanew, u, ytilde, dist,tol,rankp,lik, expected)
82
83                    call pytheta(thetanew, dist, u, yt, ymiss, dev, p, n)
84                    if(rankp .NE. m) then
85                        call pthetarest(thetanew, timevar, zt, tt, a1, p, m, n, dev, kt,kinf,ft,finf,dt,jt)
86                    end if
87                end do
88            else !cannot correct step size as we have just began
89                info = 1
90                return
91            end if
92        end if
93
94        if(finitex(dev).EQ.0) then !non-finite value of objective function
95            if(k .GT. 1) then
96                kk = 0
97                do while(finitex(dev).EQ.0)
98                    kk = kk + 1
99                    if(kk .GT. maxiter) then !did not find valid likelihood
100                        info = 2
101                        return
102                    end if
103
104                    theta = 0.5d0*(thetaold+theta)
105                    call approxloop(yt, ymiss, timevar, zt, tt, rtv, ht, qt, rqr, tvrqr, a1, p1,p1inf, p,n,m,r, &
106                    theta, thetanew, u, ytilde, dist,tol,rankp,lik, expected)
107
108                    call pytheta(thetanew, dist, u, yt, ymiss, dev, p, n)
109                    if(rankp .NE. m) then
110                        call pthetarest(thetanew, timevar, zt, tt, a1, p, m, n, dev, kt,kinf,ft,finf,dt,jt)
111                    end if
112
113                end do
114            else !cannot correct step size as we have just began
115                info = 2
116                return
117            end if
118        end if
119
120
121        ! decreasing deviance
122        if((dev - devold)/(0.1d0 + abs(dev)) .LT. -convtol .AND. k  .GT.  1) then
123            kk = 0
124            do while((dev - devold)/(0.1d0 + abs(dev)) .LT. convtol .AND. kk .LT. maxiter)
125                kk = kk + 1
126                ! previous theta produced too 'big' thetanew
127                ! new guess by halving the last try
128                theta = 0.5d0*(thetaold+theta)
129                call approxloop(yt, ymiss, timevar, zt, tt, rtv, ht, qt, rqr, tvrqr, a1, p1,p1inf, p,n,m,r, &
130                theta, thetanew, u, ytilde, dist,tol,rankp,lik, expected)
131
132                call pytheta(thetanew, dist, u, yt, ymiss, dev, p, n)
133                if(rankp .NE. m) then
134                    call pthetarest(thetanew, timevar, zt, tt, a1, p, m, n, dev, kt,kinf,ft,finf,dt,jt)
135                end if
136
137            end do
138        end if
139
140        diff = abs(dev - devold)/(0.1d0 + abs(dev))
141        if(diff .LT. convtol) then !convergence
142            theta=thetanew
143            info=0
144            exit
145        else
146            thetaold=theta
147            theta=thetanew
148            devold = dev
149        end if
150    end do
151    if(maxiter.EQ.k) then
152        info=3
153    end if
154    maxiter=k
155
156    ! check if the the approximation lead to potentially zero signal-to-noise ratio
157    if(maxval(ht) > htol) then
158      info = -5
159      htol = maxval(ht)
160      return
161    end if
162
163end subroutine approx
164