kern_reg.f90

The subroutines glkern.f and lokern.f use an efficient and fast algorithm for automatically adapti

      IF(s(j-1) < wwl) THEN
        jl=j+1
      ELSE
        IF(s(j) > wwr) EXIT
        CALL dreg(sw,a1,a2,iord,x(j),s(j-1),s(j),tt(i),wid,1)
      END IF
    END DO
    jr=j-1
    wr=wwr
  END IF
!-
!------ IF BANDWIDTH IS TOO SMALL NO SMOOTHING
  6666 IF(wwl >= s(jr-1) .AND. wwr <= s(jr)) THEN
    y(i)=x(jr)
    IF(nue > 0) y(i)=0.d0
  ELSE
!-
!------ ADD FIRST AND LAST POINT OF THE SMOOTHING INTERVAL
    DO  k=1,iord
      xf(k)=sw(k)
    END DO
    IF(jl /= 1) CALL dreg(xf,a1,a2,iord,x(jl-1),wwl,s(jl-1),tt(i),wid,1)
    IF(jr /= n) CALL dreg(xf,a1,a2,iord,x(jr+1),s(jr),wwr,tt(i),wid,1)
!-
!------ NOW THE SUMS ARE BUILT THAT ARE NEEDED TO COMPUTE THE ESTIMATE
    CALL freg(xf,nue,kord,iboun,y(i),c,icall,a)
    IF(nue > 0) y(i)=y(i)/(wid**nue)
  END IF
!-
!------ NEW INITIALISATION ?
  IF(jl > jr .OR. wwl > wr .OR. init > 100) init=0
  wido=wid
!-
END DO
!-
RETURN
END SUBROUTINE kernfa



SUBROUTINE kernfp(t, x, n, b, nue, kord, ny, s, tt, m, y)
!-----------------------------------------------------------------------
!       SHORT-VERSION: JANUARY, 1997

!       PURPOSE:

!       COMPUTATION OF KERNEL ESTIMATE USING O(N) ALGORITHM BASED ON
!       LEGENDRE POLYNOMIALS, GENERAL SPACED DESIGN AND LOCAL
!       BANDWIDTH ALLOWED. (NEW INITIALISATIONS OF THE LEGENDRE SUMS
!       FOR NUMERICAL REASONS) WITHOUT BOUNDARY KERNELS, JUST NORMALIZING

!  PARAMETERS :

!  INPUT    T(N)         INPUT GRID
!  INPUT    X(N)         DATA, GIVEN ON T(N)
!  INPUT    N            LENGTH OF X
!  INPUT    B            ONE SIDED BANDWIDTH
!  INPUT    NUE          ORDER OF DERIVATIVE (0-4)
!  INPUT    KORD         ORDER OF KERNEL (<=6)
!  INPUT    NY           0, GLOBAL BANDWIDTH; 1, LOCAL BANDWIDTH IN Y
!  INPUT    S(0:N)       HALF POINT INTERPOLATION SEQUENCE
!  INPUT    TT(M)        OUTPUT GRID
!  INPUT    M            NUMBER OF POINTS TO ESTIMATE
!  INPUT    Y(M)         BANDWITH SEQUENCE FOR NY=1, DUMMY FOR NY=0
!  OUTPUT   Y(M)         ESTIMATED FUNCTION

!-----------------------------------------------------------------------

REAL (dp), INTENT(IN)      :: t(:)
REAL (dp), INTENT(IN)      :: x(:)
INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN)      :: b
INTEGER, INTENT(IN)        :: nue
INTEGER, INTENT(IN)        :: kord
INTEGER, INTENT(IN)        :: ny
REAL (dp), INTENT(IN)      :: s(0:n)
REAL (dp), INTENT(IN)      :: tt(:)
INTEGER, INTENT(IN)        :: m
REAL (dp), INTENT(IN OUT)  :: y(:)


INTEGER    :: j, k, iord, init, icall, i, iboun
INTEGER    :: jl, jr, jnr, jnl

REAL (dp)  :: c(7), sw(7), xf(7), dold, qq, q, xnor
REAL (dp)  :: a(7,7), a1(7), a2(7), a3(7,7), cm(7,6)
REAL (dp)  :: bmin, bmax, bb, s0, sn, wwl, wwr, wid, wr, wido
!-
!------ COMPUTE CONSTANTS FOR LATER USE
s0=1.5*t(1)-0.5*t(2)
sn=1.5*t(n)-0.5*t(n-1)
bmin=(sn-s0)*.6D0/DBLE(n)*DBLE(kord-1)
bmax=(s(n)-s(0))*.5
IF(kord == 2) bmin=bmin*0.1D0
iord=kord+1
CALL coffi(nue,kord,c)
DO  k=3,iord
  a1(k)=DBLE(2*k-1)/DBLE(k)
  a2(k)=DBLE(1-k)/DBLE(k)
END DO
!-
init=0
icall=0
dold=0.d0
!-
!------ SMOOTHING LOOP
DO  i=1,m
  bb=b
  IF (ny == 1) bb=y(i)
  IF(bb < bmin) bb=bmin
  IF(bb > bmax) bb=bmax
  iboun=0
!-
!------ COMPUTE LEFT BOUNDARY
  IF(tt(i) < s(0)+bb) THEN
    wwl=s(0)
    wwr=s(0)+bb+bb
    wid=wwr-tt(i)
    iboun=1
  END IF
!-
!------ COMPUTE RIGHT BOUNDARY
  IF(tt(i)+bb > s(n)) THEN
    wwl=s(n)-(bb+bb)
    wwr=s(n)
    wid=tt(i)-wwl
    iboun=-1
  END IF
!-
!------ NO BOUNDARY
  IF(iboun == 0) THEN
    wid=bb
    wwl=tt(i)-bb
    wwr=tt(i)+bb
    xnor=1.d0
  END IF
!-
!------ COMPUTE NORMALIZING CONSTANT
  IF(iboun /= 0) THEN
    IF(iboun == 1) q=(tt(i)-s(0))/wid
    IF(iboun == -1) q=(s(n)-tt(i))/wid
    qq=q*q
    xnor=c(1)*(1.d0+q)
    DO  k=3,iord,2
      q=q*qq
      xnor=xnor+c(k)*(1.d0+q)
    END DO
    iboun=0
  END IF
!-
!------ INITIALISATION FOR INIT=0
  IF(init == 0) THEN
    DO  k=1,iord
      sw(k)=0.
    END DO
    jl=1
    DO  j=1,n
      IF(s(j-1) < wwl) THEN
        jl=j+1
      ELSE
        IF(s(j) > wwr) EXIT
        CALL dreg(sw,a1,a2,iord,x(j),s(j-1),s(j),tt(i),wid,1)
      END IF
    END DO
    jr=j-1
    wr=wwr
    init=1
    GO TO 6666
  ELSE
    init=init+1
  END IF
!-
!------ COMPARE OLD SUM WITH NEW SMOOTHING INTERVALL TT(I)-B,TT(I)+B
  IF(s(jr-1) >= wwl) THEN
    jnr=jr
    jnl=jl
    IF(s(jr) > wwr) THEN
      DO  j=jr,jl,-1
        CALL dreg(sw, a1, a2, iord, x(j), s(j-1), s(j), tt(i-1), wido, -1)
        jnr=j-1
        IF(s(jnr) <= wwr) EXIT
      END DO
    END IF
    IF(s(jl-1) < wwl) THEN
      DO  j=jl,jr
        CALL dreg(sw, a1, a2, iord, x(j), s(j-1), s(j), tt(i-1), wido, -1)
        jnl=j+1
        IF(s(j) >= wwl) EXIT
      END DO
    END IF
!-
!------ UPDATING OF SW
    CALL lreg(sw, a3, iord, (tt(i)-tt(i-1))/wid, dold, wido/wid, cm)
    IF(jnr == jr) THEN
      DO  j=jr+1,n
        IF(s(j) > wwr) EXIT
        CALL dreg(sw, a1, a2, iord, x(j), s(j-1), s(j), tt(i), wid, 1)
        jnr=j
      END DO
    END IF
    jr=jnr
    IF(jl == jnl) THEN
      DO   j=jl-1,1,-1
        IF(s(j-1) < wwl) EXIT
        CALL dreg(sw,a1,a2,iord,x(j),s(j-1),s(j),tt(i),wid,1)
        jnl=j
      END DO
    END IF
    jl=jnl
  ELSE
!-
!------ NEW INITIALISATION OF SW
    sw(1:iord)=0.
    DO  j=jr,n
      IF(s(j-1) < wwl) THEN
        jl=j+1
      ELSE
        IF(s(j) > wwr) EXIT
        CALL dreg(sw,a1,a2,iord,x(j),s(j-1),s(j),tt(i),wid,1)
      END IF
    END DO
    jr=j-1
    wr=wwr
  END IF
!-
!------ IF BANDWIDTH IS TOO SMALL NO SMOOTHING
  6666 IF(wwl >= s(jr-1) .AND. wwr <= s(jr)) THEN
    y(i)=x(jr)
    IF(nue > 0) y(i)=0.d0
  ELSE
!-
!------ ADD FIRST AND LAST POINT OF THE SMOOTHING INTERVAL
    DO  k=1,iord
      xf(k)=sw(k)
    END DO
    IF(jl /= 1) CALL dreg(xf,a1,a2,iord,x(jl-1),wwl,s(jl-1),tt(i),wid,1)
    IF(jr /= n) CALL dreg(xf,a1,a2,iord,x(jr+1),s(jr),wwr,tt(i),wid,1)
!-
!------ NOW THE SUMS ARE BUILT THAT ARE NEEDED TO COMPUTE THE ESTIMATE
    CALL freg(xf,nue,kord,iboun,y(i),c,icall,a)
    IF(nue > 0) y(i)=y(i)/(wid**nue)
    IF(xnor /= 1.d0) y(i)=y(i)/xnor
  END IF
!-
!------ NEW INITIALISATION ?
  IF(jl > jr .OR. wwl > wr .OR. init > 100) init=0
  wido=wid
!-
END DO
!-
RETURN
END SUBROUTINE kernfp



SUBROUTINE dreg(sw, a1, a2, iord, x, sl, sr, t, b, iflop)
!-----------------------------------------------------------------------
!    VERSION: MAY, 1995

!    PURPOSE:

!    COMPUTES NEW LEGENDRE SUMS (FOR REGRESSION)

!    PARAMETERS:
!                 **************   INPUT   *******************

!     SW(IORD)  :   OLD SUM OF DATA WEIGHTS FOR LEGENDRE POLYNOM.
!     A1(7)     :   CONSTANTS OF RECURSIVE FORMULA FOR LEGENDRE POL.
!     A2(7)     :                             "
!     IORD      :   ORDER OF KERNEL POLYNOMIAL
!     X         :   DATA POINT
!     SL        :   LEFT S-VALUE
!     SR        :   RIGHT S-VALUE
!     T         :   POINT WHERE THE SMOOTHED VALUE IS TO BE ESTIMATED
!     B         :   BANDWIDTH
!     IFLOP     :   1: ADDITION, ELSE SUBTRACTION


!                 **************   OUTPUT   *******************

!     SW(IORD)  :   NEW SUM OF DATA WEIGHTS FOR LEGENDRE POLYNOM.

!-----------------------------------------------------------------------

REAL (dp), INTENT(OUT)  :: sw(7)
REAL (dp), INTENT(IN)   :: a1(7)
REAL (dp), INTENT(IN)   :: a2(7)
INTEGER, INTENT(IN)     :: iord
REAL (dp), INTENT(IN)   :: x
REAL (dp), INTENT(IN)   :: sl
REAL (dp), INTENT(IN)   :: sr
REAL (dp), INTENT(IN)   :: t
REAL (dp), INTENT(IN)   :: b
INTEGER, INTENT(IN)     :: iflop


INTEGER    :: k
REAL (dp)  :: p(7,2)
!-
!------  COMPUTE LEGENDRE POLYNOMIALS
p(1,1)=(t-sl)/b
p(1,2)=(t-sr)/b
p(2,1)=1.5D0*p(1,1)*p(1,1)-.5D0
p(2,2)=1.5D0*p(1,2)*p(1,2)-.5D0
DO  k=3,iord
  p(k,1)=a1(k)*p(k-1,1)*p(1,1)+a2(k)*p(k-2,1)
  p(k,2)=a1(k)*p(k-1,2)*p(1,2)+a2(k)*p(k-2,2)
END DO
!-
!------  COMPUTE NEW LEGENDRE SUMS
IF(iflop == 1) THEN
  DO  k=1,iord
    sw(k)=sw(k)+(p(k,1)-p(k,2))*x
  END DO
ELSE
  DO  k=1,iord
    sw(k)=sw(k)+(p(k,2)-p(k,1))*x
  END DO
END IF
RETURN
END SUBROUTINE dreg



SUBROUTINE lreg(sw, a3, iord, d, dold, q, c)
!------------------------------------------------------------------
!    VERSION: MAY, 1995

!    PURPOSE:

!    UPDATE OF SW-SEQUENCE ACCORDING TO NEW BANDWIDTH AND NEW DATA
!    (VERSION FOR REGRESSION)

!    PARAMETERS :
!                       **************   INPUT   *******************

!            SW(IORD)  :  SUM OF DATA WEIGHTS FOR LEGENDRE POLYNOM.
!            IORD      :  ORDER OF KERNEL POLYNOMIAL
!            D         :  DIST. TO THE NEXT POINT DIVIDED BY BANDW.
!            DOLD      :  D       PREVIOUS STEP
!            Q         :  NEW BANDWIDTH DIVIDED BY OLD BANDWIDTH
!                       **************    WORK   *******************

!            A3(7,7)   :  MATRIX (P*Q*P)**(-1)
!            C(7,6)    :  MATRIX OF COEFFICIENTS
!                       **************   OUTPUT   ******************

!            SW(IORD)  :  UPDATED VERSION OF SW

!---------------------------------------------------------------------

REAL (dp), INTENT(IN OUT)  :: sw(7)
REAL (dp), INTENT(OUT)     :: a3(7,7)
INTEGER, INTENT(IN)        :: iord
REAL (dp), INTENT(IN)      :: d
REAL (dp), INTENT(IN OUT)  :: dold
REAL (dp), INTENT(IN)      :: q
REAL (dp), INTENT(OUT)     :: c(7,6)


INTEGER    :: k, i, l
REAL (dp)  :: dd, ww, qq, xx

!-
!- BUILD UP MATRIX
IF(dold /= d .OR. dold == 0) THEN
  dold=d
  dd=d*d
!-
  IF(iord == 7) THEN
    c(7,6)=13.d0*d
    c(7,5)=71.5D0*dd
    c(7,4)=(214.5D0*dd+9.d0)*d
    c(7,3)=(375.375D0*dd+77.d0)*dd
    c(7,2)=((375.375D0*dd+247.5D0)*dd+5.d0)*d
    c(7,1)=((187.6875D0*dd+346.5D0)*dd+40.5D0)*dd
  END IF
!-
  IF(iord >= 6) THEN
    c(6,5)=11.d0*d
    c(6,4)=49.5D0*dd
    c(6,3)=(115.5D0*dd+7.d0)*d
    c(6,2)=(144.375D0*dd+45.d0)*dd
    c(6,1)=((86.625D0*dd+94.5D0)*dd+3.d0)*d
  END IF
!-
  IF(iord >= 5) THEN
    c(5,4)=9.d0*d
    c(5,3)=31.5D0*dd
    c(5,2)=(52.5D0*dd+5.d0)*d
    c(5,1)=(39.375D0*dd+21.d0)*dd
  END IF
!-
  IF(iord >= 4) THEN
    c(4,3)=7.d0*d
    c(4,2)=17.5D0*dd
    c(4,1)=(17.5D0*dd+3.d0)*d
  END IF
!-
  IF(iord >= 3) THEN
    c(3,2)=5.d0*d
    c(3,1)=7.5D0*dd
  END IF
!-
  c(2,1)=3.d0*d
END IF
IF(q < .9999 .OR. q > 1.0001) THEN
!-
!------- BUILT UP MATRIX A3=P*Q*P**-1
  a3(1,1)=q
  DO  k=2,iord
    a3(k,k)=a3(k-1,k-1)*q
  END DO
  ww=q*q-1.d0
  DO  k=1,iord-2
    ww=ww*q
    a3(k+2,k)=(k+.5D0)*ww
  END DO
!-
  IF(iord >= 5) THEN
    qq=a3(2,2)
    a3(5,1)=q*(1.875D0+qq*(-5.25D0+qq*3.375D0))
  END IF
  IF(iord >= 6) a3(6,2)=qq*(4.375D0+qq*(-11.25D0+qq*6.875D0))
  IF(iord == 7) THEN
    a3(7,1)=q*(-2.1875D0+qq*(11.8125D0+qq*(-18.5625D0+qq* 8.9375D0)))
    a3(7,3)=q*qq*(7.875D0+qq*(-19.25D0+qq*11.375D0))
  END IF
!-
!------- COMPUTE A*C AND NEW LEGENDRE SUMS
  DO  i=iord,2,-1
    xx=0.
    DO  k=1,i
      ww=0.
      DO  l=k,i-1,2
        ww=ww+a3(l,k)*c(i,l)
      END DO