anchor.diff:  
       ! further compress data into 'sparse' multiplicities
       n_multi = count(multi0 > 0)
  -    allocate(hns(n_multi+1),stat=err)
  -    allocate(hzs(n_multi+1),stat=err)
  -    hns(1) = 0
  -    hzs(1) = n_empty_bins
  +    allocate(hn(n_multi+1),stat=err)
  +    allocate(hz(n_multi+1),stat=err)
  +    hn(1) = 0
  +    hz(1) = n_empty_bins
       k_ = 1
       do i_ = 1, nmax
          if (multi0(i_) > 0) then
             k_ = k_ + 1
  -          hns(k_) = i_
  -          hzs(k_) = multi0(i_)
  +          hn(k_) = i_
  +          hz(k_) = multi0(i_)
          end if
       end do
       deallocate(multi0)
   
       allocate(phi(n_multi+1), stat=err)
   
  -    n_data = sum(hzs * hns)
  +    n_data = sum(hz * hn)
   
     end subroutine initialize_from_counts
   
  +  subroutine initialize_from_multiplicities(hn1, hz1)
  +    ! set n_multi, hn, multi
  +    use constants
  +    real(real64), intent(in) :: hn1(:)
  +    real(real64), intent(in) :: hz1(:)
  +    integer(int32)              :: err
  +
  +
  +    allocate(hn, source=hn1, stat=err)
  +    allocate(hz, source=hz1, stat=err)
  +
  +    allocate(phi, mold=hn1, stat=err)
  +
  +    alphabet_size = sum(hz)
  +    n_data = sum(hz * hn)
  +
  +  end subroutine initialize_from_multiplicities
  +
     subroutine finalize()
   
  -    if (allocated(hzs)) then
  -       deallocate(hzs)
  +    if (allocated(hz)) then
  +       deallocate(hz)
       end if
   
  -    if (allocated(hns)) then
  -       deallocate(hns)
  +    if (allocated(hn)) then
  +       deallocate(hn)
       end if
   
       if (allocated(phi)) then
  
positive.code:     real(real64) :: c, nu, ni, xi, jsum
  
      nu = n_data + alpha * alphabet_size
      phi = digamma(hns + alpha + one) - &
           digamma(nu + two)
      c = trigamma(nu + two)
  
      h_var = 0.0
      do i_ = 0, size(hzs)-1
         ni = hns(i_) + alpha
         xi = phi(i_)
         jsum = sum(hzs * ni * (hns + alpha) * &
              (xi * phi - c))
         h_var = h_var + hzs(i_) * jsum
         h_var = h_var - hzs(i_) * ni**2 * (xi**2 - c)
         xi = xi + 1 / (ni + one)
         h_var = h_var + hzs(i_) * (ni + one) * ni * &
              (xi**2 + trigamma(ni + two) - c)
      end do
  
  
negative.code:          - alphabet_size * log_gamma(alpha) &
           - log_gamma(n_data + alpha * alphabet_size)
  
      wsum = sum(hz * (log_gamma(hn + alpha) - log_gamma(hn + one)))
  
      log_pna = log_pna + wsum
  
  