Merge pull request ralna#84 from ralna/83-fix-params

83 fix params: Removes `params_box_type` and methods (pushed vars into new workspace)

libRALFit/doc/C/examples.rst

@@ -9,4 +9,4 @@ Examples
 
 This returns the following output:
 
-.. literalinclude:: ../../test/c_example_output.out
+.. literalinclude:: ../../test/nlls_c_example.output

libRALFit/doc/Fortran/examples.rst

@@ -9,4 +9,4 @@ Examples
 
 This returns the following output:
 
-.. literalinclude:: ../../test/fortran_example_output.out
+.. literalinclude:: ../../test/nlls_fortran_example.output

libRALFit/doc/common/method.rst

@@ -272,7 +272,7 @@ minimum of the model with a regularization term added:
 .. math::
    :label: regsub
    
-   \iter{\vs} = \arg \min_{\vs} \ \iter{m} (\vs)  + \frac{1}{\Delta_k}\cdot \frac{1}{p} \|\vs\|_B^p,
+   \iter{\vs} = \arg \min_{\vs} \ \iter{m} (\vs)  + \frac{1}{\Delta_k}\cdot \frac{1}{q} \|\vs\|_B^q,
 
 At present, only one method of solving
 this subproblem is supported:
@@ -281,7 +281,7 @@ this subproblem is supported:
   this solves :eq:`regsub` by first
   converting the problem into the form
 
-  .. math:: \min_\vp \vw^T \vp + \frac{1}{2} \vp^T \vD \vp + \frac{1}{p}\|\vp\|_2^p,
+  .. math:: \min_\vp \vw^T \vp + \frac{1}{2} \vp^T \vD \vp + \frac{1}{q}\|\vp\|_2^q,
 
   where :math:`\vD` is a diagonal matrix. We do this by performing an
   eigen-decomposition of the Hessian in the model. Then, we call the
@@ -410,10 +410,11 @@ The package supports two options:
 Incorporating the regularization term
 -------------------------------------
 
-If a non-zero regularization term is required in
-:eq:`lsq`, then this is handled by transforming the
-problem internally into a new non-linear least-squares problem. 
-The formulation used will depend on the value of ``regularization`` in |nlls_options|.
+If ``regularization = 0``, then no regularization is added.
+
+A non-zero regularization term can be used in :eq:`lsq` by setting ``regularization`` to be non-zero.
+This is done by transforming the problem internally into a new non-linear least-squares problem. 
+The formulation used will depend on the value of ``regularization`` in |nlls_options|, as described below.
 
 ``regularization = 1``
   **This is only supported if** :math:`\bf p = 2`.

libRALFit/doc/common/options.rst

@@ -96,7 +96,7 @@
 
 .. |hybrid_switch_its| replace:: if ``model=3``, sets how many iterates in a row must the condition in the definition of ``hybrid_tol`` hold before a switch.
 
-.. |reg_order| replace:: if ``type_of_method = 2``, the order of the regularization used (:math:`p` in  (:eq:`regsub`)).   If ``reg_order = 0.0``, then the algorithm chooses an appropriate value of :math:`p`. 
+.. |reg_order| replace:: if ``type_of_method = 2``, the order of the regularization used (:math:`q` in  (:eq:`regsub`)).   If ``reg_order = 0.0``, then the algorithm chooses an appropriate value of :math:`q`. 
 
 .. |inner_method| replace::  if ``nlls_method = 4``, specifies the method used to pass in the regularization parameter to the inner non-linear least squares solver.   Possible values are:
 

libRALFit/example/Fortran/LanczosBox.f90

@@ -5,7 +5,7 @@ module lanczos_box_module
   implicit none
   integer, parameter :: wp = kind(0d0)
 
-  type, extends(params_box_type) :: params_type
+  type, extends(params_base_type) :: params_type
      real(wp), dimension(:), allocatable :: t ! The m data points t_i
      real(wp), dimension(:), allocatable :: y ! The m data points y_i
   end type params_type
@@ -192,17 +192,17 @@ program lanczos_box
   options%box_tau_tr_step = 0.3_wp
   options%box_ls_step_maxit = 20
 
-  Call nlls_setup_bounds(params, n, blx, bux, options, inform)
-  if (inform%status/=0) then
-    Write(*,*) 'ERROR: nlls_setup_bounds failed, status=', inform%status
-    stop
-  End if
-
+!!$  Call nlls_setup_bounds(params, n, blx, bux, options, inform)
+!!$  if (inform%status/=0) then
+!!$    Write(*,*) 'ERROR: nlls_setup_bounds failed, status=', inform%status
+!!$    stop
+!!$  End if
 
 
   ! call fitting routine
   call cpu_time(tic)
-  call nlls_solve(n,m,x,eval_r, eval_J, eval_HF, params, options, inform)
+  call nlls_solve(n,m,x,eval_r, eval_J, eval_HF, params, options, inform, &
+       lower_bounds=blx, upper_bounds=bux)
   if(inform%status.ne.0) then
      print *, "ral_nlls() returned with error flag ", inform%status
      stop
@@ -212,15 +212,9 @@ program lanczos_box
   ! Print result
   Write(*,*) 'Solution: '
   Write(*,Fmt=99998) 'idx', 'low bnd', 'x', 'upp bnd'
-  if (params%iusrbox/=0) Then
-    Do i = 1, n
-      Write(*,Fmt=99999) i, blx(i), x(i), bux(i)
-    End Do
-  else
-    Do i = 1, n
-      Write(*,Fmt=99997) i, x(i)
-    End Do
-  End If
+  Do i = 1, n
+     Write(*,Fmt=99999) i, blx(i), x(i), bux(i)
+  End Do
   print *, ""
   print *, "Objective Value at solution    = ", inform%obj
   print *, "Objective Gradient at solution = ", inform%norm_g

libRALFit/src/ral_nlls_bounds.f90

@@ -1,103 +1,36 @@
 Module ral_nlls_bounds
 Implicit None
 Private
-! User routine
-Public :: nlls_setup_bounds
 ! Routines used by ral_nlls_internal
 Public :: box_proj, box_projdir
 
 Contains
 
-Subroutine nlls_setup_bounds(params, n, blx, bux, options, inform)
-    Use ral_nlls_workspaces
-    Implicit None
-    Integer, Intent(In)                       :: n
-    Type(NLLS_inform), Intent(InOut)          :: inform
-    Type(NLLS_options), Intent(In)            :: options
-    Class(params_box_type), Intent(InOut)     :: params
-    Real(Kind=wp), Intent(In)                 :: blx(n), bux(n)
-    Integer                                   :: i, iusrbox, ierr
-
-    Continue
-
-    params%prjchd = .False.
-    params%quad_i = 0
-    params%quad_c = 0.0_wp
-    params%quad_q = 0.0_wp
-
-    iusrbox = 0
-    Do i = 1, n
-      If ( blx(i) <= bux(i) .And. blx(i)==blx(i) .And. bux(i)==bux(i) ) Then
-        If ( blx(i) > -options%box_bigbnd .Or. options%box_bigbnd > bux(i)) Then
-          iusrbox = 1
-        End If
-      Else
-        inform%status = NLLS_ERROR_BAD_BOX_BOUNDS
-        Go To 100
-      End If
-    End Do
-
-    params%iusrbox = iusrbox
-    If (iusrbox==1) Then
-      params%iusrbox = 1
-      ! Clear all arrays...
-      If (allocated(params%blx)) deallocate(params%blx, Stat=ierr)
-      If (allocated(params%bux)) deallocate(params%bux, Stat=ierr)
-      If (allocated(params%pdir)) deallocate(params%pdir, Stat=ierr)
-      If (allocated(params%normFref)) deallocate(params%normFref, Stat=ierr)
-      If (allocated(params%sk)) deallocate(params%sk, Stat=ierr)
-      If (allocated(params%g)) deallocate(params%g, Stat=ierr)
-      Allocate(params%blx(n), params%bux(n), params%pdir(n), params%g(n),      &
-        params%normFref(options%box_nFref_max), params%sk(n), Stat=ierr)
-      if (ierr /= 0) Then
-        If (allocated(params%blx)) deallocate(params%blx, Stat=ierr)
-        If (allocated(params%bux)) deallocate(params%bux, Stat=ierr)
-        If (allocated(params%pdir)) deallocate(params%pdir, Stat=ierr)
-        If (allocated(params%normFref)) deallocate(params%normFref, Stat=ierr)
-        If (allocated(params%sk)) deallocate(params%sk, Stat=ierr)
-        If (allocated(params%g)) deallocate(params%g, Stat=ierr)
-        inform%status = NLLS_ERROR_ALLOCATION
-        inform%bad_alloc = 'ral_nlls_box'
-        Go To 100
-      end if
-      params%normfref(1:options%box_nFref_max) = -1.0e-20_wp
-      Do i = 1, n
-        params%blx(i) = max(blx(i), -options%box_bigbnd)
-        params%bux(i) = min(bux(i), options%box_bigbnd)
-      End Do
-    End If
-
-100 Continue
-
-End Subroutine nlls_setup_bounds
-
-
-
-  Subroutine box_proj(params, n, x, xnew, dir, alpha)
-    Use ral_nlls_workspaces, Only: params_box_type, wp
+Subroutine box_proj(w, n, x, xnew, dir, alpha)
+    Use ral_nlls_workspaces, Only: box_type, wp
 !   Two modes
 !   If xnew and d are present, then project x+alpha*dir, otherwise just 
 !   make x feasible (ignoring either dir or xnew) and return
 !   In either case flag in params%prjchd if the projection altered any entry
     Implicit None
     Integer, Intent(In)                    :: n
-    Class(params_box_type), Intent(InOut)  :: params
+    type(box_type), Intent(InOut)          :: w
     Real(Kind=wp), Intent(InOut)           :: x(n)
     Real(Kind=wp), Intent(InOut), Optional :: xnew(n)
     Real(Kind=wp), Intent(In), Optional    :: dir(n), alpha
     Integer                                :: i
     Real(Kind=wp)                          :: alp, xi
 
     Continue
-    params%prjchd = .False.
+    w%prjchd = .False.
     If (.not. (present(xnew) .And. present(dir))) Then
 !     Make feasible point x and return
-      If (params%iusrbox==1) Then
+      If (w%iusrbox==1) Then
         Do i = 1, n
           xi = x(i)
-          x(i) = max(min(params%bux(i), x(i)), params%blx(i))
+          x(i) = max(min(w%bux(i), x(i)), w%blx(i))
           If (xi/=x(i)) Then
-            params%prjchd = .True.
+            w%prjchd = .True.
           End If
         End Do
       End If
@@ -110,11 +43,11 @@ Subroutine box_proj(params, n, x, xnew, dir, alpha)
       alp = 1.0_wp
     End If
 
-    If (params%iusrbox==1) Then
+    If (w%iusrbox==1) Then
       Do i = 1, n
-        xnew(i) = max(min(params%bux(i), x(i)+alp*dir(i)), params%blx(i))
+        xnew(i) = max(min(w%bux(i), x(i)+alp*dir(i)), w%blx(i))
         If (xnew(i) /= x(i)+dir(i)) Then
-          params%prjchd = .True.
+          w%prjchd = .True.
         End If
       End Do
     Else
@@ -125,42 +58,44 @@ Subroutine box_proj(params, n, x, xnew, dir, alpha)
 
   End Subroutine box_proj
 
-  Subroutine box_projdir(params, n, x, dir, normg, sigma)
-    Use ral_nlls_workspaces, Only: params_box_type, wp
-!   Calculate the projected dir and it's two-norm
-!   Assumes dir = -fdx
-!   If there is no box, then normPD=normg and if pdir is allocated then 
-!   copy dir to it.
+  Subroutine box_projdir(w, n, x, dir, normg, sigma)
+    Use ral_nlls_workspaces, Only: wp, box_type
+    !   Calculate the projected dir and it's two-norm
+    !   Assumes dir = -fdx
+    !   If there is no box, then normPD=normg and if pdir is allocated then 
+    !   copy dir to it.
     Implicit None
+    type( box_type ), Intent(InOut)       :: w
     Integer, Intent(In)                   :: n
-    Class(params_box_type), Intent(InOut) :: params
     Real(Kind=wp), Intent(In)             :: x(n), dir(n), normg
     Real(Kind=wp), Intent(In), Optional   :: sigma
+
     Real(Kind=wp)                         :: alpb
     Integer                               :: i
-
-    If (params%iusrbox==1) Then
-      If (Present(sigma)) Then
-        alpb = sigma
-      Else
-        alpb = 1.0_wp
-      End If
-      params%normPD = 0.0_wp
-      do i = 1, n
-        If (params%bux(i)/=params%blx(i)) Then
-          params%pdir(i) = max(min(params%bux(i), x(i)+alpb*dir(i)), params%blx(i))-x(i)
-          params%normPD = params%normPD + (params%pdir(i))**2 
-        Else
-          params%pdir(i) = 0.0_wp
-        End If
-      end do
-      params%normPD = sqrt(params%normPD)
+    
+    If (w%iusrbox==1) Then
+       If (Present(sigma)) Then
+          alpb = sigma
+       Else
+          alpb = 1.0_wp
+       End If
+       w%normPD = 0.0_wp
+       do i = 1, n
+          If (w%bux(i)/=w%blx(i)) Then
+             w%pdir(i) = max(min(w%bux(i), x(i)+alpb*dir(i)), w%blx(i))-x(i)
+             w%normPD = w%normPD + (w%pdir(i))**2 
+          Else
+             w%pdir(i) = 0.0_wp
+          End If
+       end do
+       w%normPD = sqrt(w%normPD)
     Else
-      If (allocated(params%pdir)) Then
-        params%pdir(1:n) = dir(1:n)
-      End If
-      params%normPD = normg
+       If (allocated(w%pdir)) Then
+          w%pdir(1:n) = dir(1:n)
+       End If
+       w%normPD = normg
     End If
   End Subroutine box_projdir
-
+
+
 End Module ral_nlls_bounds

libRALFit/src/ral_nlls_double.f90

@@ -3,9 +3,8 @@
 module ral_nlls_double
 
   use ral_nlls_internal
-  use ral_nlls_workspaces, only : params_base_type, params_box_type,           &
+  use ral_nlls_workspaces, only : params_base_type,           &
     nlls_options, nlls_inform, nlls_workspace
-  use ral_nlls_bounds, Only: nlls_setup_bounds
 
   implicit none
 
@@ -50,7 +49,7 @@ end subroutine eval_hf_type
 
   public :: nlls_solve, nlls_iterate, nlls_finalize, nlls_strerror
   public :: nlls_options, nlls_inform, nlls_workspace
-  public :: params_base_type, params_box_type, nlls_setup_bounds
+  public :: params_base_type
 
 end module ral_nlls_double