conmin_class Derived Type

    type, public :: conmin_class

        !! main class.

        ! private  ! for now these are public so they can be set by user [should replace with an init method]

        real(wp) :: small = 1.0e-20_wp !! Small number
        real(wp) :: large = 1.0e+20_wp !! Large number

        ! cnmn1 variables
        real(wp) :: delfun = 0.001_wp   !! Minimum relative change in the objective
                                        !! function to indicate convergence.  If in `ITRM` consecutive
                                        !! iterations, `ABS(1.0-OBJ(J-1)/OBJ(J))<DELFUN` and the current
                                        !! design is feasible (all `G(J)<=ABS(CT)`), the minimization
                                        !! process is terminated.  If the current design is infeasible
                                        !! (some `G(J)>ABS(CT)`), five iterations are required to
                                        !! terminate and this situation indicates that a feasible design
                                        !! may not exist.
        real(wp) :: dabfun = 0.0_wp !! Default value = 0.001 times the initial function value.  Same
                                    !! as DELFUN except comparison is on absolute change in the
                                    !! objective function, `ABS(OBJ(J)-OBJ(J-1))`, instead of relative
                                    !! change.
        real(wp) :: fdch = 0.01_wp   !! Not used if `NFDG = 0`.  Relative change in
                                     !! decision variable `X(I)` in calculating finite difference
                                     !! gradients.  For example, `FDCH = 0.01` corresponds to a finite
                                     !! difference step of one percent of the value of the decision
                                     !! variable.
        real(wp) :: fdchm = 0.01_wp !! Not used if `NFDG = 0`.  Minimum absolute
                                    !! step in finite difference gradient calculations.  FDCHM applies
                                    !! to the unscaled variable values.
        real(wp) :: ct = -0.1_wp    !! Not used if `NCON = NSIDE = 0`.
                                    !! Constraint thickness parameter.  If `CT<=G(J)<=ABS(CT)`,
                                    !! `G(J)` is defined as active.  If `G(J)>ABS(CT)`, `G(J)` is said to
                                    !! be violated.  If `G(J)<CT`, `G(J)` is not active.  CT is
                                    !! sequentially reduced in magnitude during the optimization
                                    !! process.  If `ABS(CT)` is very small, one or more constraints
                                    !! may be active on one iteration and inactive on the next,
                                    !! only to become active again on a subsequent iteration.
                                    !! This is often referred to as "zigzagging" between constraints.
                                    !! A wide initial value of the constraint thickness is desirable
                                    !! for highly nonlinear problems so that when a constraint
                                    !! becomes active it tends to remain active, thus reducing the
                                    !! zigzagging problem.  The default value is usually adequate.
        real(wp) :: ctmin = 0.004_wp    !! Not used if `NCON = NSIDE = 0`.  Minimum
                                        !! absolute value of CT considered in the optimization process.
                                        !! CTMIN may be considered as "numerical zero" since it may not be
                                        !! meaningful to compare numbers smaller than CTMIN.  The value of
                                        !! CTMIN is chosen to indicate that satisfaction of a constraint
                                        !! within this tolerance is acceptable.  The default value is usually
                                        !! adequate.
        real(wp) :: ctl = -0.01_wp  !! Not used if `NCON = NSIDE = 0`.
                                    !! Constraint thickness parameter for linear and side constraints.
                                    !! CTL is smaller in magnitude than CT because the zigzagging
                                    !! problem is avoided with linear and side constraints.  The default
                                    !! value is usually adequate.
        real(wp) :: ctlmin = 0.001_wp   !! Not used if `NCON = NSIDE = 0`.  Minimum
                                        !! absolute value of CTL considered in the optimization process.
                                        !! The default value is usually adequate.
        real(wp) :: alphax = 0.1_wp  !! the maximum fractional change in any
                                     !! component of X as an initial estimate for ALPHA in the one-dimensional
                                     !! search. That is, the initial ALPHA will be such that no component of X is changed by
                                     !! more than this amount. This only applies to those X(i) of magnitude greater than 0.1.
                                     !! If an optimization run shows numerous ALPHA = 0 results for the one-dimensional search,
                                     !! it may help to try ALPHAX less than the default. ALPHAX is changed by CONMIN depending
                                     !! on the progress of the optimization.
        real(wp) :: abobj1 = 0.1_wp  !! the fractional change attempted as a first step in
                                     !! the one-dimensional search and is based on a linear approximation.
                                     !! ABOBJ1 is updated during the optimization, depending on progress.
                                     !! The initial step in the one-dimensional search is taken as the amount
                                     !! necessary to change `OBJ` by `ABOBJ1*ABS(OBJ)` or to change some `X(i)` by `ALPHAX*ABS( X(i) )`,
                                     !! whichever is less.
        real(wp) :: theta = 1.0_wp !! Not used if `NCON = NSIDE = 0`.  Mean value
                                    !! of the push-off factor in the method of feasible directions.
                                    !! A larger value of THETA is desirable if the constraints, `G(J)`,
                                    !! are known to be highly nonlinear, and a smaller value may be
                                    !! used if all `G(J)` are known to be nearly linear.  The actual
                                    !! value of the push-off factor used in the program is a quadratic
                                    !! function of each `G(J)`, varying from 0.0 for `G(J) = CT` to `4.0*THETA`
                                    !! for `G(J) = ABS(CT)`.  A value of `THETA = 0.0` is used in the
                                    !! program for constraints which are identified by the user to be
                                    !! strictly linear.  THETA is called a "push-off" factor because
                                    !! it pushes the design away from the active constraints into the
                                    !! feasible region.  The default value is usually adequate.
        real(wp) :: obj !! Value of objective function for the current decision variables,
                        !! `X(I), I = 1, NDV` contained in vector `X`.  Calculate OBJ if
                        !! INFO = 1 or INFO = 2.
        integer :: ndv !! Number of decision variables, `X(I)`, contained in vector `X`.
        integer :: ncon !! Number of constraint functions, `G(J)`.  NCON may be zero.
        integer :: nside    !! Side constraint parameter.  `NSIDE = 0` signifies that the
                            !! variables `X(I)` do not have lower or upper bounds.  `NSIDE>0`
                            !! signifies that all variables `X(I)` have lower and upper bounds
                            !! defined by `VLB(I)` and `VUB(I)` respectively.  If one or more
                            !! variables are not bounded while others are, the values of the
                            !! lower and upper bounds on the unbounded variables must be taken
                            !! as very large negative and positive values respectively
                            !! (i.e., `VLB(I) = -1.0E+10`, `VUB(I) = 1.0E+10`).
        integer :: iunit = output_unit !! Unit number for printing. Should be non-zero.
        integer :: iprint = 0   !! Print control.  All printing is done on unit number `iunit`.
                                !!
                                !!  * `iprint=0`:  Print nothing.
                                !!  * `iprint=1`:  Print initial and final function information.
                                !!  * `iprint=2`:  1st debug level.  Print all of above plus control
                                !!      parameters.  Print function value and X-vector at each
                                !!      iteration.
                                !!  * `iprint=3`:  2nd. debug level.  Print all of above plus all constraint
                                !!      values, numbers of active or violated constraints, direction
                                !!      vectors, move parameters and miscellaneous information.  The
                                !!      constraint parameter, BETA, printed under this option
                                !!      approaches zero as the optimum objective is achieved.
                                !!  * `iprint=4`:  Complete debug.  Print all of above plus gradients of
                                !!      objective function, active or violated constraint functions
                                !!      and miscellaneous information.
                                !!  * `iprint=5`:  all of above plus each proposed design vector, objective
                                !!      and constraints during the one-dimensional search.
        integer :: nfdg !! the finite difference gradient parameter:
                        !!
                        !!  * `NFDG = 0`:  all gradient information is calculated by finite difference
                        !!    within [[conmin]].
                        !!  * `NFDG = 1`:  all gradient information is supplied by the user.
                        !!  * `NFDG = 2`:  the gradient of `OBJ` is supplied by the user and the
                        !!    gradients of constraints are calculated by finite
                        !!    difference within [[conmin]].
        integer :: nscal    !! Scaling control parameter.  The decision variables will be
                            !! scaled linearly.
                            !!
                            !!  * `NSCAL<0`:  Scale variables `X(I) `by dividing by `SCAL(I)`, where
                            !!                vector SCAL is defined by the user.
                            !!  * `NSCAL==0`: Do not scale the variables.
                            !!  * `NSCAL>0`:  Scale the variables every NSCAL iterations.
                            !!                Variables are normalized so that scaled
                            !!                `X(I) = X(I)/ABS(X(I))`.  When using this option, it
                            !!                is desirable that `NSCAL = ICNDIR` if ICNDIR is input
                            !!                as nonzero, and `NSCAL = NDV + 1` in ICNDIR is input
                            !!                as zero.
        integer :: linobj = 0   !! Not used if NCON = NSIDE = 0.  Linear objective function
                                !! identifier.  If the objective, OBJ, is specifically known to
                                !! be a strictly linear function of the decision variables, `X(I)`,
                                !! set `LINOBJ = 1`.  If OBJ is a general nonlinear function, set
                                !! `LINOBJ = 0`.
        integer :: itmax = 10 !! Maximum number of iterations in the
                              !! minimization process.  If `NFDG==0` each iteration requires one
                              !! set of gradient computations (INFO = 3 or 4) and approximately
                              !! three function evaluations (INFO = 1 or 2).  If `NFDG>0`
                              !! each iteration requires approximately NDV + 3 function
                              !! evaluations (INFO = 1 or 2).
        integer :: itrm = 3 !! Number of consecutive iterations to indicate
                            !! convergence by relative or absolute changes, `DELFUN` or `DABFUN`.
        integer :: icndir !! Default value = `NDV + 1`.  Conjugate direction restart parameter.
                          !! If the function is currently unconstrained, (all `G(J)<CT` or
                          !! `NCON = NSIDE = 0`), Fletcher-Reeves conjugate direction method will
                          !! be restarted with a steepest descent direction every ICNDIR
                          !! iterations.  If `ICNDIR = 1` only steepest descent will be used.
        integer :: igoto = 0 !! Reverse communication flag.
                             !! CONMIN executes according to the parameter IGOTO which must be initialized to zero.
        integer :: nac  !! Number of active and violated constraints (`G(J)>=CT`).
                        !! Calculate `NAC` if `INFO = 4` and `NFDG = 0`.
        integer :: info !! * `INFO = 1`:  calculate OBJ and G(I), I = 1, NCON
                        !! * `INFO = 2`:  calculate NAC, IC(I), I = 1,NAC, the gradient of OBJ, and
                        !!   the gradient of G(J), where J = IC(I), I = 1,NAC.
                        !!   Store the gradients of G in columns of A.
        integer :: infog !! * `INFOG = 0`:   same as when INFOG was not used.
                         !! * `INFOG = 1`:   only those constraints identified as active or violated
                         !!                in array IC(I), I = 1, NAC need be evaluated.  This is
                         !!                only meaningful if finite difference gradients are
                         !!                calculated, and allows the user to avoid
                         !!                calculating non-essential information.  If it is
                         !!                convenient to evaluate all constraints each time,
                         !!                variable INFOG may be ignored.
        integer :: iter !! Iteration number.  The optimization process is iterative so
                        !! that the vector of decision variables at the Kth iteration
                        !! is defined by `X(K) = X(K - 1) + ALPHA*S(K)`, where in this case
                        !! `K` refers to the iteration number and the components `X(I)` are
                        !! all changed simultaneously.  `ALPHA` is defined as the move
                        !! parameter and is printed if the print control `IPRINT>=3`.
                        !! `S` is the move direction.

        ! consav variables
        ! [note: removed dm1, dm2, idm1, etc. as no longer necessary]
        real(wp) :: dct, dctl, abobj, cta, ctam, ctbm, obj1, &
                    slope, dx, dx1, fi, xi, dftdf1, alp, fff, a1, a2, a3, &
                    a4, f1, f2, f3, f4, cv1, cv2, cv3, cv4, app, alpca, &
                    alpfes, alpln, alpmin, alpnc, alpsav, alpsid, alptot, &
                    rspace
        real(wp) :: phi = 5.0_wp    !! Not used if `NCON = NSIDE = 0`.
                                    !! Participation coefficient, used if a design is infeasible
                                    !! (one or more `G(J)>ABS(CT)`).  PHI is a measure of how hard
                                    !! the design will be "pushed" towards the feasible region and
                                    !! is, in effect, a penalty parameter.  If in a given problem, a
                                    !! feasible solution cannot be obtained with the default value,
                                    !! PHI should be increased, and the problem run again.  If a
                                    !! feasible solution cannot be obtained with `PHI = 100`, it is
                                    !! probable that no feasible solution exists.  The default value
                                    !! is usually adequate.
        integer :: jdir, iobj, kobj, kcount, &
                   nfeas, mscal, ncobj, nvc, kount, icount, igood1, igood2, &
                   igood3, igood4, ibest, iii, nlnc, jgoto, ispace(2)
        integer :: ncal(2)  !! Bookkeeping information.  NCAL(1) gives the number of times
                            !! external routine SUB1 was called with INFO = 1.  NCAL(2) gives
                            !! the number of times INFO = 2.  NCAL(3) gives the number of times
                            !! INFO = 3 and NCAL(4) gives the number of times INFO = 4.

        procedure(report_f),pointer :: report => null() !! a function to call once per iteration to report the solution

    contains
        private
        procedure, public :: solve => conmin
        procedure :: cnmn01
        procedure :: cnmn02
        procedure :: cnmn03
        procedure :: cnmn04
        procedure :: cnmn05
        procedure :: cnmn06
        procedure :: cnmn07
    end type conmin_class