! Mark 28.6 Release. NAG Copyright 2022.
Module e04tcfe_mod
! .. Use Statements ..
Use iso_c_binding, Only: c_f_pointer, c_ptr
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Accessibility Statements ..
Private
Public :: lsqfun, lsqgrd
! .. Derived Type Definitions ..
Type, Public :: usr_data
Integer :: nres = 0
Real (Kind=nag_wp), Allocatable :: t(:), y(:)
End Type usr_data
Contains
Subroutine lsqfun(nvar,x,nres,rx,inform,iuser,ruser,cpuser)
! .. Implicit None Statement ..
Implicit None
! .. 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 ..
Type (usr_data), Pointer :: ud
Integer :: i
! .. Intrinsic Procedures ..
Intrinsic :: sin
! .. Executable Statements ..
! Unpack the user data from cpuser
Call c_f_pointer(cpuser,ud)
Do i = 1, nres
rx(i) = (x(1)*ud%t(i)**2+x(2)*ud%t(i)+x(3)+x(4)*sin(x(5)*ud%t(i))) - &
ud%y(i)
End Do
End Subroutine lsqfun
Subroutine lsqgrd(nvar,x,nres,nnzrd,rdx,inform,iuser,ruser,cpuser)
! .. Implicit None Statement ..
Implicit None
! .. 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 ..
Type (usr_data), Pointer :: ud
Integer :: i
! .. Intrinsic Procedures ..
Intrinsic :: cos, sin
! .. Executable Statements ..
! Unpack the user data from cpuser
Call c_f_pointer(cpuser,ud)
Do i = 1, nres
rdx((i-1)*nvar+1) = ud%t(i)**2
rdx((i-1)*nvar+2) = ud%t(i)
rdx((i-1)*nvar+3) = 1.0_nag_wp
rdx((i-1)*nvar+4) = sin(x(5)*ud%t(i))
rdx((i-1)*nvar+5) = x(4)*ud%t(i)*cos(x(5)*ud%t(i))
End Do
End Subroutine lsqgrd
End Module e04tcfe_mod
Program e04tcfe
! This example shows how to use the NAG Optimization Modeling Suite
! to edit a nonlinear least squares problem.
! In particular, the outlier residuals are disabled without having to
! redefine the full problem and comparing the solutions is easy.
!
! we try to fit the following model with 5 parameters:
! f(t) = at^2 + bt + c + d sin(omega t)
! to some noisy data (30 points) with 2 outliers (point 10 and 20)
! The data was generated using the values:
! a = 0.3
! b = 1.0
! c = 0.01
! d = 0.2
! omega = 5.0
! .. Use Statements ..
Use e04tcfe_mod, Only: lsqfun, lsqgrd, usr_data
Use iso_c_binding, Only: c_loc, c_null_ptr, c_ptr
Use nag_library, Only: e04ffu, e04ggf, e04ggu, e04ggv, e04raf, e04rmf, &
e04rzf, e04tbf, e04tcf, e04tdf, e04zmf, nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Integer, Parameter :: nin = 5, nout = 6
! .. Local Scalars ..
Type (c_ptr) :: cpuser, handle
Type (usr_data), Target :: ud
Integer :: ifail, isparse, nnzrd, nres, nvar
! .. Local Arrays ..
Real (Kind=nag_wp) :: rinfo(100), ruser(1), stats(100)
Real (Kind=nag_wp), Allocatable :: rx(:), x(:)
Integer :: icolrd(1), irowrd(1), iuser(1)
! .. Executable Statements ..
Write (nout,*) 'E04TCF Example Program Results'
Write (nout,*)
handle = c_null_ptr
cpuser = c_null_ptr
! Skip Header in data file
Read (nin,*)
! Read number of residuals
Read (nin,*) nres
ud%nres = nres
! Allocate memory
Allocate (ud%t(nres),ud%y(nres))
! Read observations
Read (nin,*) ud%t(1:nres)
Read (nin,*) ud%y(1:nres)
! try to fit the model
! f(t) = at^2 + bt + c + d sin(omega t)
! To the data {(t_i, y_i)}
nvar = 5
! Initialize the NAG optimization handle
ifail = 0
Call e04raf(handle,nvar,ifail)
! Define a dense nonlinear least-squares objective function
! (isparse = 0 => the sparsity pattern of the Jacobian doesn't need to be
! defined)
isparse = 0
nnzrd = 1
Call e04rmf(handle,nres,isparse,nnzrd,irowrd,icolrd,ifail)
! Set some optional parameters to control the output of the solver
Call e04zmf(handle,'Print Options = No',ifail)
Call e04zmf(handle,'Print Solution = X',ifail)
Call e04zmf(handle,'Print Level = 1',ifail)
Write (nout,*) 'First solve the problem with the outliers'
Write (nout,*)
Write (nout,*) &
'--------------------------------------------------------'
! Call the solver
Allocate (x(nvar),rx(nres))
x(1:nvar) = 1.0_nag_wp
cpuser = c_loc(ud)
ifail = -1
Call e04ggf(handle,lsqfun,lsqgrd,e04ggu,e04ggv,e04ffu,nvar,x,nres,rx, &
rinfo,stats,iuser,ruser,cpuser,ifail)
Write (nout,*) &
'--------------------------------------------------------'
Write (nout,*)
Write (nout,*) &
'Now remove the outlier residuals from the problem handle'
Write (nout,*)
Write (nout,*) &
'--------------------------------------------------------'
! Disable the two outlier residuals
Call e04tcf(handle,'NLS',2,(/10,20/),ifail)
! Solve the problem again
x(1:nvar) = 1.0_nag_wp
ifail = -1
Call e04ggf(handle,lsqfun,lsqgrd,e04ggu,e04ggv,e04ffu,nvar,x,nres,rx, &
rinfo,stats,iuser,ruser,cpuser,ifail)
Write (nout,*) &
'--------------------------------------------------------'
Write (nout,*)
Write (nout,*) 'Assuming the outliers points are measured again'
Write (nout,*) 'we can enable the residuals and adjust the values'
Write (nout,*)
Write (nout,*) &
'--------------------------------------------------------'
! Fix the first variable to its known value of 0.3
! enable the residuals and adjust the values in the data
Call e04tdf(handle,'variable',1,0.3_nag_wp,0.3_nag_wp,ifail)
Call e04tbf(handle,'NLS',2,(/10,20/),ifail)
ud%y(10) = -0.515629_nag_wp
ud%y(20) = 0.54920_nag_wp
! Solve the problem
x(1:nvar) = 1.0_nag_wp
ifail = -1
Call e04ggf(handle,lsqfun,lsqgrd,e04ggu,e04ggv,e04ffu,nvar,x,nres,rx, &
rinfo,stats,iuser,ruser,cpuser,ifail)
! Free the handle memory
ifail = 0
Call e04rzf(handle,ifail)
End Program e04tcfe