tau_exp can be changed when calling uwerr()

doc/aderrors.tex

@@ -58,7 +58,7 @@
 data. It uses the
 $\Gamma$-method together with techniques of automatic differentiation to provide
 a robust, efficient, portable, transparent and open source module for
-the analysis of Monte Carlo data~\cite{aderrors-mod}.
+the analysis of Monte Carlo data~\cite{Ramos:2018vgu}.
 \end{abstract}
 
 \tableofcontents
@@ -691,9 +691,12 @@ This example shows the calls to the different accesible components
 
 
 \begin{description}
-\item[\texttt{uwerr()}: ] Subroutine. Performs the error analysis. Only
+\item[\texttt{uwerr([ids\_texp])}: ] Subroutine. Performs the error analysis. Only
   after calling this method we can expect that \texttt{error(), derror(),
-    taui(n), dtaui(n), window(n)} will return the proper values.
+    taui(n), dtaui(n), window(n)} will return the proper values. The
+  optional argument \texttt{ids\_texp} select the ensemble ID's to
+  attach a tail to the autocorrelation function (by default all
+  ensembles with $\tau_{\rm exp}>0$ are analyzed in this way).
 \item[\texttt{error\_src(n)}: ] Function. Returns the contribution to the sum of
   errors in quadrature of the n$^{\rm th}$ MC chain contributing to
   the error.

src/aderrors/aderrors_cf.f90

@@ -79,16 +79,40 @@ contains
 
 ! ****************************************
 ! *
-  module subroutine uwerr(p)
+  module subroutine uwerr(p, ids_texp)
 ! *
 ! ****************************************
     class (uwreal), intent (inout) :: p
-
-    integer :: n, irst, irnd
-
+    integer, intent (in), optional :: ids_texp(:)
+    
+    integer :: n, k, irst, irnd
+    logical :: do_texp(p%nid)
+
+    if (present(ids_texp)) then
+       do_texp = .false.
+       do n = 1, p%nid
+          do k = 1, size(ids_texp)
+             if (p%id(n).eq.ids_texp(k)) then
+                if (p%texp(n).gt.0.0_DP) then
+                   do_texp(n) = .true.
+                else
+                   do_texp(n) = .false.
+                end if
+                exit
+             end if
+          end do
+       end do
+    else
+       where (p%texp.gt.0.0_DP)
+          do_texp = .true.
+       elsewhere
+          do_texp = .false.
+       end where
+    end if
+    
     call init_cf(p)
     call autocorrfunc_obs(p)
-    call wopt_texp(p)
+    call wopt_texp(p,do_texp)
     call tauint(p)
    
     p%var  = 0.0_DP
@@ -101,12 +125,19 @@ contains
           irst = sum(p%tmax(1:n-1))
           irnd = irst + p%iw(n)
           
-          p%var(n) = p%gamm(irst)/real(p%nd(n),kind=dp) + &
-               2.0_dp * sum(p%gamm(irst+1:irnd))/real(p%nd(n),kind=dp)
-          if (p%texp(n) > 0.0_dp) p%var(n) = p%var(n) + &
-               2.0_dp*p%texp(n)*p%gamm(irnd)/real(p%nd(n),kind=dp)
-
-          p%derr = p%derr + &
+          if ( (do_texp(n))   .and. &
+               (p%gamm(irnd) > 0.0_DP).and. &
+               (p%iw(n) > 0)         ) then
+             
+             p%var(n) = p%gamm(irst)/real(p%nd(n),kind=dp) + &
+                  2.0_dp * sum(p%gamm(irst+1:irnd-1))/real(p%nd(n),kind=dp) + &
+                  2.0_dp*p%texp(n)*p%gamm(irnd)/real(p%nd(n),kind=dp)
+          else
+             p%var(n) = p%gamm(irst)/real(p%nd(n),kind=dp) + &
+                  2.0_dp * sum(p%gamm(irst+1:irnd))/real(p%nd(n),kind=dp)
+          end if
+             
+          if (p%iw(n).gt.0) p%derr = p%derr + &
                p%var(n)*(real(p%iw(n),kind=DP)-0.5_DP) / real(p%nd(n),kind=DP)
        end if
     end do
@@ -263,30 +294,36 @@ contains
              rw  = real(i-ist,kind=DP)
              tiw = tiw + p%gamm(i)/p%gamm(ist)
              if (tiw < 0.5_DP) then
-                tau = 0.01_DP
+                tau =  0.01_DP
+                gw  = -1.0_DP
              else
-                tau = p%s(n)/log((2.0_DP*tiw+1.0_DP)/(2.0_DP*tiw-1.0_DP))
+                if (tiw.gt.0.5_DP) then
+                   tau = p%s(n)/log((2.0_DP*tiw+1.0_DP)/(2.0_DP*tiw-1.0_DP))
+                   gw  = exp(-rw/tau) - tau/sqrt(rw*real(p%nd(n),kind=DP))
+                else
+                   tau = 1.0_DP
+                   gw  = 1.0_DP
+                end if
              end if
-             gw  = exp(-rw/tau) - tau/sqrt(rw*real(p%nd(n),kind=DP))
 
              if (gw < 0.0_DP) exit
           end do
-          p%iw(n) = i-ist
+          p%iw(n) = min(i-ist,ind-ist)
        end if
     end do
-
+    
     return
   end subroutine wopt
   
 ! ********************************  
 ! *
-  subroutine wopt_texp(p)
+  subroutine wopt_texp(p, do_texp)
 ! *
 ! ********************************
     type (uwreal), intent (inout) :: p
+    logical, intent (in)          :: do_texp(p%nid)
     
     integer :: i=0, n, ist, ind
-    Real (kind=DP) :: tiw
 
     do n = 1, p%nid
        if (p%nd(n).eq.1) then
@@ -294,13 +331,15 @@ contains
           cycle
        end if
        call p%get_cf_offset(ist,ind,n)
-       
-       tiw = 0.5_DP
-       if (p%texp(n) > 0.0_DP) then
+
+       if (do_texp(n)) then
+          p%iw(n) = 0
+          if (p%gamm(ist).eq.0.0_DP) cycle
           do i = ist+1, ind
+             if (p%gamm(i).lt.0.0_DP)   exit
              if (p%dsig(n)*p%drho(i) > p%gamm(i)/p%gamm(ist)) exit
           end do
-          p%iw(n) = i-ist
+          p%iw(n) = i-ist-1
        end if
     end do
 

src/aderrors/aderrors_cov.f90

@@ -64,12 +64,20 @@ contains
              else
                 irst = sum(wsid%tmax(1:n-1))
                 irnd = irst + wsid%iw(n)
-                
-                mat(io1,io2) = mat(io1,io2) + &
-                     wscf%gamm(irst,io1,io2)/real(wsid%nd(n),kind=dp) + &
-                     2.0_DP * sum(wscf%gamm(irst+1:irnd,io1,io2))/real(wsid%nd(n),kind=dp)
-                if (wsid%texp(n) > 0.0_dp) mat(io1,io2) = mat(io1,io2) + &
-                     2.0_dp*wsid%texp(n)*wscf%gamm(irnd,io1,io2)/real(wsid%nd(n),kind=dp)
+
+                if ( (wsid%texp(n)            > 0.0_DP).and. &
+                     (wscf%gamm(irnd,io1,io2) > 0.0_DP).and. &
+                     (wsid%iw(n)              > 0)         ) then
+
+                   mat(io1,io2) = mat(io1,io2) + &
+                        wscf%gamm(irst,io1,io2)/real(wsid%nd(n),kind=dp) + &
+                        2.0_DP * sum(wscf%gamm(irst+1:irnd-1,io1,io2))/real(wsid%nd(n),kind=dp) + &
+                        2.0_dp*wsid%texp(n)*wscf%gamm(irnd,io1,io2)/real(wsid%nd(n),kind=dp)
+                else
+                   mat(io1,io2) = mat(io1,io2) + &
+                        wscf%gamm(irst,io1,io2)/real(wsid%nd(n),kind=dp) + &
+                        2.0_DP * sum(wscf%gamm(irst+1:irnd,io1,io2))/real(wsid%nd(n),kind=dp)
+                end if
              end if
           end do
        end do
@@ -91,7 +99,6 @@ contains
 ! ********************************
     integer, intent (in) :: is, ir, io
     
-    if ((wsid%ns.ge.is).and.(wsid%no.ge.io).and.(wsid%nr.gt.ir)) return
     if (allocated(wsid%map))   deallocate(wsid%map)
     if (allocated(wsid%ids))   deallocate(wsid%ids)
     if (allocated(wsid%nrep))  deallocate(wsid%nrep)
@@ -144,8 +151,9 @@ contains
     if ((maxrl.gt.wscf%maxrl).or.(nobs.gt.wscf%no)) then
        if (allocated(wscf%gamm)) deallocate(wscf%gamm)
        if (allocated(wscf%data)) deallocate(wscf%data)
-       allocate(wscf%gamm(0:maxrl-1,nobs,nobs), &
-            wscf%data(2*maxrl,nobs), wscf%drho(0:maxrl-1,nobs) )
+       if (allocated(wscf%drho)) deallocate(wscf%drho)
+       allocate(wscf%gamm(0:maxrl/2-1,nobs,nobs), &
+            wscf%data(maxrl,nobs), wscf%drho(0:maxrl/2-1,nobs) )
        wscf%maxrl = maxrl
     end if
     
@@ -239,9 +247,12 @@ contains
        do io = 1, wsid%no          
           tiw = 0.5_DP
           if (wsid%texp(n) > 0.0_DP) then
-             do i = ist+1, ind
-                if (SFT_FAC*wscf%drho(i,io) > wscf%gamm(i,io,io)/wscf%gamm(ist,io,io)) exit
-             end do
+             wsid%iw(n) = 0
+             if (wscf%gamm(ist,io,io).gt.0.0_DP) then 
+                do i = ist+1, ind
+                   if (SFT_FAC*wscf%drho(i,io) > wscf%gamm(i,io,io)/wscf%gamm(ist,io,io)) exit
+                end do
+             end if
              wsid%iw(n) = max(wsid%iw(n),i-ist)
           else
              if (wscf%gamm(ist,io,io)==0.0_DP) then
@@ -252,10 +263,16 @@ contains
                    tiw = tiw + wscf%gamm(i,io,io)/wscf%gamm(ist,io,io)
                    if (tiw < 0.5_DP) then
                       tau = 0.01_DP
+                      gw  = -1.0_DP
                    else
-                      tau = wsid%s(n)/log((2.0_DP*tiw+1.0_DP)/(2.0_DP*tiw-1.0_DP))
+                      if (tiw.gt.0.5_DP) then
+                         tau = wsid%s(n)/log((2.0_DP*tiw+1.0_DP)/(2.0_DP*tiw-1.0_DP))
+                         gw  = exp(-rw/tau) - tau/sqrt(rw*real(wsid%nd(n),kind=DP))
+                      else
+                         tau = 1.0_DP
+                         gw  = 1.0_DP
+                      end if
                    end if
-                   gw  = exp(-rw/tau) - tau/sqrt(rw*real(wsid%nd(n),kind=DP))
                    
                    if (gw < 0.0_DP) exit
                 end do
@@ -291,7 +308,7 @@ contains
     call combine_ids_multi(nid,p)
     wsid%nid = nid
 
-    maxrl = sum(wsid%tmax(1:nid))
+    maxrl = sum(wsid%nd(1:nid))
     call init_wscf(0, no, maxrl, n1, n2)
     if (maxrl.eq.0) return