Add argument passing to nl regression routine

Author: jchristopherson

examples/nl_regression_example.f90

@@ -1,11 +1,12 @@
 module nl_example
     use iso_fortran_env
 contains
-    subroutine exfun(x, p, f, stop)
+    subroutine exfun(x, p, f, stop, args)
         ! Arguments
         real(real64), intent(in) :: x(:), p(:)
         real(real64), intent(out) :: f(:)
         logical, intent(out) :: stop
+        class(*), intent(inout), optional :: args
 
         ! Function
         f = p(4) * x**3 + p(3) * x**2 + p(2) * x + p(1)

src/fstats_regression.f90

@@ -137,7 +137,7 @@ module fstats_regression
     end type
 
     interface
-        subroutine regression_function(xdata, params, f, stop)
+        subroutine regression_function(xdata, params, f, stop, args)
             !! Defines the interface of a subroutine computing the function
             !! values at each of the N data points as part of a regression
             !! analysis.
@@ -154,6 +154,9 @@ subroutine regression_function(xdata, params, f, stop)
                 !! set to true, the iteration process will terminate.  If set
                 !! to false, the iteration process will continue along as 
                 !! normal.
+            class(*), intent(inout), optional :: args
+                !! An optional parameter that can be used to pass data 
+                !! in and out of the routine.
         end subroutine
 
         subroutine iteration_update(iter, funvals, resid, params, step)
@@ -686,7 +689,7 @@ function calculate_regression_statistics(resid, params, c, alpha, err) &
 
 ! ------------------------------------------------------------------------------
 subroutine jacobian(fun, xdata, params, &
-    jac, stop, f0, f1, step, err)
+    jac, stop, f0, f1, step, args, err)
     !! Computes the Jacobian matrix for a nonlinear regression problem.
     procedure(regression_function), intent(in), pointer :: fun
         !! A pointer to the regression_function to evaluate.
@@ -711,6 +714,9 @@ subroutine jacobian(fun, xdata, params, &
     real(real64), intent(in), optional :: step
         !! The differentiation step size.  The default is the square 
         !! root of machine precision.
+    class(*), intent(inout), optional :: args
+        !! An optional parameter that can be used to pass data 
+        !! in and out of the routine.
     class(errors), intent(inout), optional, target :: err
         !! A mechanism for communicating errors and warnings to the 
         !! caller.  Possible warning and error codes are as follows.
@@ -762,7 +768,7 @@ subroutine jacobian(fun, xdata, params, &
         allocate(f0a(m), stat = flag)
         if (flag /= 0) go to 20
         f0p(1:m) => f0a
-        call fun(xdata, params, f0p, stop)
+        call fun(xdata, params, f0p, stop, args)
         if (stop) return
     end if
     if (present(f1)) then
@@ -800,7 +806,7 @@ subroutine jacobian(fun, xdata, params, &
 ! ------------------------------------------------------------------------------
 subroutine nonlinear_least_squares(fun, x, y, params, ymod, &
     resid, weights, maxp, minp, stats, alpha, controls, settings, info, &
-    status, cov, err)
+    status, cov, args, err)
     !! Performs a nonlinear regression to fit a model using a version
     !! of the Levenberg-Marquardt algorithm.
     procedure(regression_function), intent(in), pointer :: fun
@@ -853,6 +859,9 @@ subroutine nonlinear_least_squares(fun, x, y, params, ymod, &
     real(real64), intent(out), optional, dimension(:,:) :: cov
         !! An optional N-by-N matrix that, if supplied, will be used to return
         !! the covariance matrix.
+    class(*), intent(inout), optional :: args
+        !! An optional parameter that can be used to pass data 
+        !! in and out of the [[fun]] routine.
     class(errors), intent(inout), optional, target :: err
         !! A mechanism for communicating errors and warnings to the 
         !! caller.  Possible warning and error codes are as follows.
@@ -1035,7 +1044,7 @@ subroutine nonlinear_least_squares(fun, x, y, params, ymod, &
 
     ! Process
     call lm_solve(fun, x, y, params, w, pmax, pmin, tol, opt, ymod, &
-        resid, JtWJ, inf, stop, errmgr, status)
+        resid, JtWJ, inf, stop, errmgr, status, args = args)
 
     ! Compute the covariance matrix
     if (present(stats) .or. present(cov)) then
@@ -1124,7 +1133,7 @@ subroutine lm_set_default_settings(x)
 ! - stop: A flag allowing the user to terminate model execution
 ! - work: A workspace array for the model parameters (N-by-1)
 subroutine jacobian_finite_diff(fun, xdata, params, f0, jac, f1, &
-    stop, step, work)
+    stop, step, work, args)
     ! Arguments
     procedure(regression_function), intent(in), pointer :: fun
     real(real64), intent(in) :: xdata(:), params(:)
@@ -1133,6 +1142,7 @@ subroutine jacobian_finite_diff(fun, xdata, params, f0, jac, f1, &
     real(real64), intent(out) :: f1(:), work(:)
     logical, intent(out) :: stop
     real(real64), intent(in) :: step
+    class(*), intent(inout), optional :: args
 
     ! Local Variables
     integer(int32) :: i, n
@@ -1147,7 +1157,7 @@ subroutine jacobian_finite_diff(fun, xdata, params, f0, jac, f1, &
     work = params
     do i = 1, n
         work(i) = work(i) + step
-        call fun(xdata, work, f1, stop)
+        call fun(xdata, work, f1, stop, args)
         if (stop) return
 
         jac(:,i) = (f1 - f0) / step
@@ -1218,7 +1228,7 @@ subroutine broyden_update(pOld, yOld, jac, p, y, dp, dy)
 ! - mwork: A workspace matrix (N-by-M)
 ! - update: Reset to false if a Jacobian evaluation was performed.
 subroutine lm_matrix(fun, xdata, ydata, pOld, yOld, dX2, jac, p, weights, &
-    neval, update, step, JtWJ, JtWdy, X2, yNew, stop, work, mwork)
+    neval, update, step, JtWJ, JtWdy, X2, yNew, stop, work, mwork, args)
     ! Arguments
     procedure(regression_function), pointer :: fun
     real(real64), intent(in) :: xdata(:), ydata(:), pOld(:), yOld(:), &
@@ -1231,6 +1241,7 @@ subroutine lm_matrix(fun, xdata, ydata, pOld, yOld, dX2, jac, p, weights, &
     real(real64), intent(out) :: X2, mwork(:,:), yNew(:)
     logical, intent(out) :: stop
     real(real64), intent(out), target :: work(:)
+    class(*), intent(inout), optional :: args
 
     ! Local Variables
     integer(int32) :: m, n
@@ -1243,15 +1254,15 @@ subroutine lm_matrix(fun, xdata, ydata, pOld, yOld, dX2, jac, p, weights, &
     w2(1:n) => work(m+1:n+m)
 
     ! Perform the next function evaluation
-    call fun(xdata, p, yNew, stop)
+    call fun(xdata, p, yNew, stop, args)
     neval = neval + 1
     if (stop) return
 
     ! Update or recompute the Jacobian matrix
     if (dX2 > 0 .or. update) then
         ! Recompute the Jacobian
         call jacobian_finite_diff(fun, xdata, p, yNew, jac, w1, &
-            stop, step, w2)
+            stop, step, w2, args)
         neval = neval + n
         if (stop) return
         update = .false.
@@ -1308,7 +1319,7 @@ subroutine lm_matrix(fun, xdata, ydata, pOld, yOld, dX2, jac, p, weights, &
 ! - err: An error handling mechanism
 subroutine lm_iter(fun, xdata, ydata, p, neval, niter, update, lambda, &
     maxP, minP, weights, JtWJ, JtWdy, h, pNew, deltaY, yNew, X2, X2Old, &
-    alpha, stop, iwork, err, status)
+    alpha, stop, iwork, err, status, args)
     ! Arguments
     procedure(regression_function), pointer :: fun
     real(real64), intent(in) :: xdata(:), ydata(:), p(:), maxP(:), &
@@ -1323,6 +1334,7 @@ subroutine lm_iter(fun, xdata, ydata, p, neval, niter, update, lambda, &
     integer(int32), intent(out) :: iwork(:)
     class(errors), intent(inout) :: err
     procedure(iteration_update), intent(in), pointer, optional :: status
+    class(*), intent(inout), optional :: args
 
     ! Local Variables
     integer(int32) :: i, n
@@ -1364,7 +1376,7 @@ subroutine lm_iter(fun, xdata, ydata, p, neval, niter, update, lambda, &
     end do
 
     ! Update the residual error
-    call fun(xdata, pNew, yNew, stop)
+    call fun(xdata, pNew, yNew, stop, args)
     neval = neval + 1
     deltaY = ydata - yNew
     if (stop) return
@@ -1382,7 +1394,7 @@ subroutine lm_iter(fun, xdata, ydata, p, neval, niter, update, lambda, &
             pNew(i) = min(max(minP(i), p(i) + h(i)), maxP(i))
         end do
 
-        call fun(xdata, pNew, yNew, stop)
+        call fun(xdata, pNew, yNew, stop, args)
         if (stop) return
         neval = neval + 1
         deltaY = ydata - yNew
@@ -1418,7 +1430,7 @@ subroutine lm_iter(fun, xdata, ydata, p, neval, niter, update, lambda, &
 ! - stop: A flag allowing the user to terminate model execution
 ! - err: An error handling object
 subroutine lm_solve(fun, xdata, ydata, p, weights, maxP, minP, controls, &
-    opt, y, resid, JtWJ, info, stop, err, status)
+    opt, y, resid, JtWJ, info, stop, err, status, args)
     ! Arguments
     procedure(regression_function), intent(in), pointer :: fun
     real(real64), intent(in) :: xdata(:), ydata(:), weights(:), maxP(:), &
@@ -1431,6 +1443,7 @@ subroutine lm_solve(fun, xdata, ydata, p, weights, maxP, minP, controls, &
     logical, intent(out) :: stop
     class(errors), intent(inout) :: err
     procedure(iteration_update), intent(in), pointer, optional :: status
+    class(*), intent(inout), optional :: args
 
     ! Local Variables
     logical :: update
@@ -1467,12 +1480,12 @@ subroutine lm_solve(fun, xdata, ydata, p, weights, maxP, minP, controls, &
     if (flag /= 0) go to 10
 
     ! Perform an initial function evaluation
-    call fun(xdata, p, y, stop)
+    call fun(xdata, p, y, stop, args)
     neval = 1
 
     ! Evaluate the problem matrices
     call lm_matrix(fun, xdata, ydata, pOld, yOld, 1.0d0, J, p, weights, &
-        neval, update, step, JtWJ, JtWdy, X2, y, stop, work, mwork)
+        neval, update, step, JtWJ, JtWdy, X2, y, stop, work, mwork, args)
     if (stop) go to 5
     X2Old = X2
     JtWJc = JtWJ
@@ -1492,15 +1505,15 @@ subroutine lm_solve(fun, xdata, ydata, p, weights, maxP, minP, controls, &
         ! update the new parameter estimates
         call lm_iter(fun, xdata, ydata, p, neval, niter, opt%method, &
             lambda, maxP, minP, weights, JtWJc, JtWdy, h, pTry, resid, &
-            yTemp, X2Try, X2Old, alpha, stop, iwork, err, status)
+            yTemp, X2Try, X2Old, alpha, stop, iwork, err, status, args)
         if (stop) go to 5
         if (err%has_error_occurred()) return
 
         ! Update the Chi-squared estimate, update the damping parameter
         ! lambda, and, if necessary, update the matrices
         call lm_update(fun, xdata, ydata, pOld, p, pTry, yOld, y, h, dX2, &
             X2Old, X2, X2Try, lambda, alpha, nu, JtWdy, JtWJ, J, weights, &
-            niter, neval, update, step, work, mwork, controls, opt, stop)
+            niter, neval, update, step, work, mwork, controls, opt, stop, args)
         if (stop) go to 5
         JtWJc = JtWJ
 
@@ -1545,7 +1558,7 @@ subroutine lm_solve(fun, xdata, ydata, p, weights, maxP, minP, controls, &
 !
 subroutine lm_update(fun, xdata, ydata, pOld, p, pTry, yOld, y, h, dX2, &
     X2old, X2, X2try, lambda, alpha, nu, JtWdy, JtWJ, J, weights, niter, &
-    neval, update, step, work, mwork, controls, opt, stop)
+    neval, update, step, work, mwork, controls, opt, stop, args)
     ! Arguments
     procedure(regression_function), intent(in), pointer :: fun
     real(real64), intent(in) :: xdata(:), ydata(:), X2try, h(:), step, &
@@ -1559,6 +1572,7 @@ subroutine lm_update(fun, xdata, ydata, pOld, p, pTry, yOld, y, h, dX2, &
     class(iteration_controls), intent(in) :: controls
     class(lm_solver_options), intent(in) :: opt
     logical, intent(out) :: stop
+    class(*), intent(inout), optional :: args
 
     ! Local Variables
     integer(int32) :: n
@@ -1585,7 +1599,7 @@ subroutine lm_update(fun, xdata, ydata, pOld, p, pTry, yOld, y, h, dX2, &
 
         ! Recompute the matrices
         call lm_matrix(fun, xdata, ydata, pOld, yOld, dX2, J, p, weights, &
-            neval, update, step, JtWJ, JtWdy, X2, y, stop, work, mwork)
+            neval, update, step, JtWJ, JtWdy, X2, y, stop, work, mwork, args)
         if (stop) return
 
         ! Decrease lambda
@@ -1605,7 +1619,7 @@ subroutine lm_update(fun, xdata, ydata, pOld, p, pTry, yOld, y, h, dX2, &
         if (mod(niter, 2 * n) /= 0) then
             call lm_matrix(fun, xdata, ydata, pOld, yOld, -1.0d0, J, p, &
                 weights, neval, update, step, JtWJ, JtWdy, dX2, y, stop, &
-                work, mwork)
+                work, mwork, args)
             if (stop) return
         end if