! Mark 30.2 Release. NAG Copyright 2024.
Module e04gnfe_mod
! .. Use Statements ..
Use iso_c_binding, Only: c_ptr
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Accessibility Statements ..
Private
Public :: lsqfun, lsqgrd
Contains
Subroutine lsqfun(nvar,x,nres,rx,inform,iuser,ruser,cpuser)
! .. Scalar Arguments ..
Type (c_ptr), Intent (In) :: cpuser
Integer, Intent (Inout) :: inform
Integer, Intent (In) :: nres, nvar
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
Real (Kind=nag_wp), Intent (Out) :: rx(nres)
Real (Kind=nag_wp), Intent (In) :: x(nvar)
Integer, Intent (Inout) :: iuser(*)
! .. Local Scalars ..
Integer :: i
! .. Intrinsic Procedures ..
Intrinsic :: sin
! .. Executable Statements ..
rx(:) = 0.0_nag_wp
Do i = 1, nres
rx(i) = ruser(nres+i) - x(1)*ruser(i)*sin(-x(2)*ruser(i))
End Do
inform = 0
End Subroutine lsqfun
Subroutine lsqgrd(nvar,x,nres,nnzrd,rdx,inform,iuser,ruser,cpuser)
! .. Scalar Arguments ..
Type (c_ptr), Intent (In) :: cpuser
Integer, Intent (Inout) :: inform
Integer, Intent (In) :: nnzrd, nres, nvar
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: rdx(nnzrd), ruser(*)
Real (Kind=nag_wp), Intent (In) :: x(nvar)
Integer, Intent (Inout) :: iuser(*)
! .. Local Scalars ..
Integer :: i
! .. Intrinsic Procedures ..
Intrinsic :: cos, sin
! .. Executable Statements ..
rdx(:) = 0.0_nag_wp
Do i = 1, nres
rdx((i-1)*nvar+1) = -ruser(i)*sin(-x(2)*ruser(i))
rdx((i-1)*nvar+2) = ruser(i)**2*x(1)*cos(x(2)*ruser(i))
End Do
inform = 0
End Subroutine lsqgrd
End Module e04gnfe_mod
Program e04gnfe
! .. Use Statements ..
Use e04gnfe_mod, Only: lsqfun, lsqgrd
Use iso_c_binding, Only: c_null_ptr, c_ptr
Use nag_library, Only: e04gnf, e04gnu, e04gnx, e04gny, e04raf, e04rhf, &
e04rjf, e04rmf, e04rzf, e04zmf, nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Real (Kind=nag_wp), Parameter :: infbnd = 1.0E20_nag_wp
Integer, Parameter :: nout = 6
! .. Local Scalars ..
Type (c_ptr) :: cpuser, handle
Integer :: idlc, ifail, isparse, nclin, nnza, &
nnzrd, nres, nvar, x_idx
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: a(:), bla(:), blx(:), bua(:), &
bux(:), ruser(:), rx(:), x(:)
Real (Kind=nag_wp) :: rinfo(100), stats(100)
Integer, Allocatable :: icola(:), irowa(:)
Integer :: icolrd(0), irowrd(0), iuser(0)
! .. Executable Statements ..
Write (nout,*) 'e04gnf Example Program Results'
Write (nout,*)
Flush (nout)
cpuser = c_null_ptr
! Problem size
nvar = 2
! Residual quantity
nres = 24
Allocate (ruser(2*nres))
! Data for model y = f(t;A,w) := A*t*sin(-w*t)
! Data generated using A = 0.1 and w = 0.8 and added random noise
! Residual 4, 12, 16 and 20 are outliers
! t(:) =
ruser(1:nres) = (/1.0E+00_nag_wp,2.0E+00_nag_wp,3.0E+00_nag_wp, &
4.0E+00_nag_wp,5.0E+00_nag_wp,6.0E+00_nag_wp,7.0E+00_nag_wp, &
8.0E+00_nag_wp,9.0E+00_nag_wp,1.0E+01_nag_wp,1.1E+01_nag_wp, &
1.2E+01_nag_wp,1.3E+01_nag_wp,1.4E+01_nag_wp,1.5E+01_nag_wp, &
1.6E+01_nag_wp,1.7E+01_nag_wp,1.8E+01_nag_wp,1.9E+01_nag_wp, &
2.0E+01_nag_wp,2.1E+01_nag_wp,2.2E+01_nag_wp,2.3E+01_nag_wp, &
2.4E+01_nag_wp/)
! y(:) =
ruser(nres+1:nres*2) = (/0.0523_nag_wp,-0.1442_nag_wp,-0.0422_nag_wp, &
1.8106_nag_wp,0.3271_nag_wp,0.4684_nag_wp,0.4593_nag_wp, &
-0.0169_nag_wp,-0.7811_nag_wp,-1.1356_nag_wp,-0.5343_nag_wp, &
-3.0043_nag_wp,1.1832_nag_wp,1.5153_nag_wp,0.7120_nag_wp, &
-2.2923_nag_wp,-1.4871_nag_wp,-1.7083_nag_wp,-0.9936_nag_wp, &
-5.2873_nag_wp,1.7555_nag_wp,2.0642_nag_wp,0.9499_nag_wp, &
-0.6234_nag_wp/)
iuser(:) = 0
! Initialize handle
ifail = 0
Call e04raf(handle,nvar,ifail)
! Define residuals structure, isparse=0 means the residual structure is
! dense => irowrd and icolrd are not accessed
isparse = 0
nnzrd = 0
ifail = 0
Call e04rmf(handle,nres,isparse,nnzrd,irowrd,icolrd,ifail)
! Set options: use least-squares loss function
Call e04zmf(handle,'NLDF Loss Function Type = L2',ifail)
! Change printed output verbosity
Call e04zmf(handle,'Print Level = 1',ifail)
Call e04zmf(handle,'Print Options = No',ifail)
! Define starting point
Allocate (x(nvar),rx(nres))
x(1:nvar) = (/0.3_nag_wp,0.7_nag_wp/)
! Define bounds
Allocate (blx(nvar),bux(nvar))
blx(1) = -1.0_nag_wp
bux(1) = infbnd
blx(2) = 0.0_nag_wp
bux(2) = 1.0_nag_wp
ifail = 0
Call e04rhf(handle,nvar,blx,bux,ifail)
! Add linear constraints
nclin = 1
nnza = 2
Allocate (bla(nclin),bua(nclin),irowa(nnza),icola(nnza),a(nnza))
bla(1) = 0.2_nag_wp
bua(1) = 1.0_nag_wp
irowa(1:nnza) = (/1,1/)
icola(1:nnza) = (/1,2/)
a(1:nnza) = (/1.0_nag_wp,1.0_nag_wp/)
idlc = 0
ifail = 0
Call e04rjf(handle,nclin,bla,bua,nnza,irowa,icola,a,idlc,ifail)
Write (nout,*) &
'First solve the problem using least squares loss function'
Write (nout,*) &
'---------------------------------------------------------'
! Call the solver
ifail = -1
Call e04gnf(handle,lsqfun,lsqgrd,e04gnx,e04gny,e04gnu,nvar,x,nres,rx, &
rinfo,stats,iuser,ruser,cpuser,ifail)
! Print solution if optimal or suboptimal solution found
If (ifail==0 .Or. ifail==50) Then
Write (nout,*)
Write (nout,99999) 'Optimal parameters X (Least Squares):'
Write (nout,99997) 'x_idx', ' Value '
Do x_idx = 1, nvar
Write (nout,99998) x_idx, x(x_idx)
End Do
End If
Write (nout,*)
Write (nout,*) 'Now solve the problem using SmoothL1 loss function'
Write (nout,*) &
'---------------------------------------------------------'
! Use robust nonlinear regression and resolve the model
Call e04zmf(handle,'NLDF Loss Function Type = SmoothL1',ifail)
! Call the solver to resolve
x(1:nvar) = (/0.3_nag_wp,0.7_nag_wp/)
ifail = -1
Call e04gnf(handle,lsqfun,lsqgrd,e04gnx,e04gny,e04gnu,nvar,x,nres,rx, &
rinfo,stats,iuser,ruser,cpuser,ifail)
! Print solution if optimal or suboptimal solution found
If (ifail==0 .Or. ifail==50) Then
Write (nout,*)
Write (nout,99999) 'Optimal parameters X (SmoothL1):'
Write (nout,99997) 'x_idx', ' Value '
Do x_idx = 1, nvar
Write (nout,99998) x_idx, x(x_idx)
End Do
End If
! Free the handle memory
ifail = 0
Call e04rzf(handle,ifail)
99999 Format (1X,A)
99998 Format (2X,I5,3X,Es7.2e1)
99997 Format (2X,A5,3X,A7)
End Program e04gnfe