Adds supports for gaps in the MC histories

README.md

@@ -27,6 +27,7 @@ robust and **exact error propagation even in iterative algorithms**.
   different parameters).
 - Support for **replicas** (i.e. several runs with the same simulation
   parameters). 
+- Irrelagular MC measurements are handled transparently.
 - Standalone **portable** implementation without any external
   dependencies.
 - **Fast** computation of autocorrelation functions with the [FFT

doc/aderrors.tex

@@ -813,6 +813,44 @@ they cannot be called on derived observables.
 
 \subsection{Externally accessible functions}
 
+\subsubsection{\texttt{aderr\_pobs(data, id\_msm, ivrep, nmsm)}}
+
+Returns a primary onservable with the measurements contained in
+\texttt{data(:)} corresponding to MC time \texttt{id\_msm(:)}. The
+total number of measurements in the ensemble and the replica vector
+are specified in \texttt{nmsm} and \texttt{id\_msm(:)} respectively. 
+
+Note that this call admits gaps in the MC history. For example if some
+ensemble has measurements at MC times $1,2,\dots,100$ for most
+observables, but some observable has only been measured at MC times
+$1,3,5,\dots,99$, this data can still be used in any expression by
+just calling
+\begin{verbatim}
+x = aderr_pobs(data(1:50),[(2*i-1,i=1,50)], [100], 100)
+\end{verbatim}
+will load in \texttt{x} the measurements stored in \texttt{data} that
+are \emph{only available for odd values of the MC time}. 
+
+Completely irregular MC histories are transparently dealt with in the code.
+
+\begin{description}
+\item[\texttt{data(:)}:] Type \texttt{real double precision} one
+  dimensional array. The MC data.
+
+\item[\texttt{id\_msm(:)}:] Type \texttt{integer} one dimensional
+  array. The measurements ID (i.e. ``trajectories'' in lattice jargon)
+  of each of the measurements.
+
+\item[\texttt{ivrep(:)}:] Type \texttt{integer} one simansional
+  array. The replica vector of the ensemble.
+
+\item[\texttt{nmsm}:] Type \texttt{integer}. The total number of
+  measurements in the ensemble. 
+
+\item[Output: ] Type \texttt{uwreal}. The primary observable. 
+\end{description}
+
+
 \subsubsection{\texttt{aderr\_new\_id(name, [id])}}
 
 Stores \texttt{name} in the database and associate it with ID

src/aderrors/aderrors.f90

@@ -254,6 +254,12 @@ module aderrors
        type (uwreal), intent  (in)    :: a(:)
        real (kind=DP), intent (in)    :: der(size(x),size(a)), mns(size(x))
      end subroutine
+
+     module subroutine addobs_ds(x, a, der, mns)
+       type (uwreal), intent  (inout) :: x
+       type (uwreal), intent  (in)    :: a(:)
+       real (kind=DP), intent (in)    :: der(:), mns
+     end subroutine
   end interface
 
   interface assignment (=)
@@ -296,14 +302,15 @@ module aderrors
      module procedure uwdobs_f, uwdobs_s
   end interface uwdobs
   interface addobs
-     module procedure addobs_f, addobs_s, addobs_d
+     module procedure addobs_f, addobs_s, addobs_d, addobs_ds
   end interface addobs
 
   public :: uwdobs, addobs
   
   interface
-     module subroutine uwerr(p)
+     module subroutine uwerr(p, ids_texp)
        class (uwreal), intent (inout) :: p
+       integer, intent (in), optional :: ids_texp(:)
      end subroutine
   end interface
 
@@ -604,7 +611,8 @@ module aderrors
      end function
   end interface
 
-  public :: aderr_set_mode, aderr_get_mode, adset_default_dsig, adset_default_stau
+  public :: aderr_set_mode, aderr_get_mode, adset_default_dsig, &
+       adset_default_stau, aderr_pobs
   
 contains
 
@@ -1247,7 +1255,33 @@ contains
 
     return
   end subroutine print_hist
-  
+
+! *****************************
+! *
+  function aderr_pobs(data, id_msm, ivrep, nmsm) result (p)
+! *
+! *****************************
+    type (uwreal) :: p
+    real (kind=DP), intent (in)   :: data(:)
+    integer, intent (in)          :: id_msm(:), ivrep(:), nmsm
+
+    integer :: i
+
+    if (size(data).ne.size(id_msm)) call module_error('[aderr_pobs]',&
+         ' Data and trajectories do not match')
+    
+    call p%init(1, [nmsm], [size(ivrep)], ivrep)
+
+    p%mean = sum(data)/real(size(data),kind=DP)
+    p%data = 0.0_DP
+    do i = 1, size(id_msm)
+       p%data(id_msm(i)) = data(i) - p%mean
+    end do
+    p%data = real(nmsm,kind=DP)/real(size(id_msm)) * p%data
+    
+    return
+  end function aderr_pobs
+
 ! ********************************
 ! *
     Subroutine module_error(routine, msg)