Logo Search packages:      
Sourcecode: octave-octgpr version File versions  Download package


c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
c Author: Jaroslav Hajek <highegg@gmail.com>
c This file is part of OctGPR.
c OctGPR is free software; you can redistribute it and/or modify
c it under the terms of the GNU General Public License as published by
c the Free Software Foundation; either version 3 of the License, or
c (at your option) any later version.
c This program is distributed in the hope that it will be useful,
c but WITHOUT ANY WARRANTY; without even the implied warranty of
c GNU General Public License for more details.
c You should have received a copy of the GNU General Public License
c along with this software; see the file COPYING.  If not, see
c <http://www.gnu.org/licenses/>.
      subroutine infgpr(ndim,nx,X,theta,nu,var,nlin,mu,RP,corr,
c purpose:      compute the prediction value, spatial derivatives 
c               and prediction variance of the GPR regressor in
c               a single spatial point. Should be used after
c               nllgpr, but *NOT* after nldgpr (on the same data).
c               note: to facilitate better precision, the vector
c               R \ y~ is not fully formed for value prediction; 
c               L \ y~ and L \ r are used instead. We use it for
c               derivatives, however, as the precision is less 
c               important there (and the converse implies multiple
c               traversals through the matrix R and thus a speed
c               tradeoff).
c arguments:
c ndim (in)     number of dimensions of input space
c nx (in)       number of training points
c X (in)        array of training input vectors
c theta (in)    length scales
c nu (in)       relative white noise. nu = sqrt(var_white/var)
c var (in)      MLE estimated global variance from nllgpr
c nlin (in)     number of linear trend variables
c mu (in)       at least nlin+1. Linear trend from nllgpr
c RP (in)       size at least 2*nx+nx*(nx+1)/2. Contains the factorization 
c               details as computed by nllgpr, after packing the
c               triangular matrix L.
c corr          subroutine to calculate correlation value and its
c               derivative. Must be declared as follows:
c               subroutine corr(t,f,d)
c               double precision t,f,d
c               f = correlation
c               d = derivative
c               end subroutine
c               The correlation should satisfy f(0) = 1 and f(+inf) = 0.
c               t >= 0 is the scaled squared norm of input vectors
c               difference, i.e. sum(theta*(X(:,i)-X(:,j))**2)
c x0 (in)       the spatial point to predict in
c y0 (out)      the prediction values.
c sig0 (out)    the prediction sigmas (noise included).
c nder(in)      number of derivatives requested. 0 to omit derivatives.
c yd0 (out)     the prediction derivatives. if nder <= 0, yd0 is not
c               referenced.
c work (out)    workspace; size at least nx*(1+min(nder,1))
      integer ndim,nx,nlin,nder
      real*8 X(ndim,nx),theta(ndim),nu,var,x0(ndim)
      real*8 mu(0:nlin),RP(*),y0,sig0,yd0(*),work(nx,*)
      external corr
      external dwdis2,dsdacc,dtpsv,dcopy,xerbla
      real*8 tmp,dwdis2
      integer i,k,info
c argument checks
      info = 0
      if (ndim < 0) then
        info = 1
      else if (nx < 1) then
        info = 2
      else if (nlin < 0 .or. nlin > ndim .or. nlin >= nx) then
        info = 7
      end if
      if (info /= 0) then
        call xerbla('infgpr',info)
      end if
c calculate correlation vector r
      do i = 1,nx
        call corr(dwdis2(ndim,theta,X(1,i),x0),work(i,1),tmp)
c only use last part of workspace if derivatives requested
        if (nder > 0) work(i,2) = tmp
      end do
c form L \ r
      call dtpsv('L','N','N',nx,RP(2*nx+1),work(1,1),1)
c accumulate sum((L\r).^2) and (L\y)'*(L\r)
      call dsdacc(nx,work(1,1),RP,sig0,tmp)
c add linear trend
      tmp = tmp + mu(0)
      do k = 1,nlin
        tmp = tmp + x0(k)*mu(k)
      end do
      y0 = tmp
c get deviation
      sig0 = sqrt((1+nu**2-sig0) * var)
c calc derivatives only if necessary
      if (nder == 0) return
      do k = 1,nder
        yd0(k) = 0
      end do
      do i = 1,nx
c calculate the primary part
        tmp = RP(nx+i)*work(i,2)
c apply chain rule
        do k = 1,nder
          yd0(k) = yd0(k) + tmp*(x0(k) - X(k,i))
        end do
      end do
c correction
      do k = 1,nder
        yd0(k) = yd0(k) * (2 * theta(k)**2)
        if (k <= nlin) yd0(k) = yd0(k) + mu(k)
      end do
      end subroutine



Generated by  Doxygen 1.6.0   Back to index