source: trunk/cppad/ipopt/solve_callback.hpp @ 3275

Last change on this file since 3275 was 3275, checked in by bradbell, 6 years ago

Better checking for if an file already included.

  • Property svn:keywords set to Id
File size: 31.2 KB
Line 
1/* $Id: solve_callback.hpp 3275 2014-05-20 17:35:11Z bradbell $ */
2# ifndef CPPAD_SOLVE_CALLBACK_INCLUDED
3# define CPPAD_SOLVE_CALLBACK_INCLUDED
4/* --------------------------------------------------------------------------
5CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell
6
7CppAD is distributed under multiple licenses. This distribution is under
8the terms of the
9                    Eclipse Public License Version 1.0.
10
11A copy of this license is included in the COPYING file of this distribution.
12Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
13-------------------------------------------------------------------------- */
14
15# include <cppad/cppad.hpp>
16# include <coin/IpIpoptApplication.hpp>
17# include <coin/IpTNLP.hpp>
18# include <cppad/ipopt/solve_result.hpp>
19
20namespace CppAD { // BEGIN_CPPAD_NAMESPACE
21namespace ipopt {
22/*!
23\file solve_callback.hpp
24\brief Class that connects ipopt::solve to Ipopt
25*/
26
27/*!
28Class that Ipopt uses for obtaining information about this problem.
29
30
31\section Evaluation_Methods Evaluation Methods
32The set of evaluation methods for this class is
33\verbatim
34        { eval_f, eval_grad_f, eval_g, eval_jac_g, eval_h }
35\endverbatim
36Note that the bool return flag for the evaluations methods
37does not appear in the Ipopt documentation.
38Looking at the code, it seems to be a flag telling Ipopt to abort
39when the flag is false.
40*/
41template <class Dvector, class ADvector, class FG_eval>
42class solve_callback : public Ipopt::TNLP
43{
44private:
45        // ------------------------------------------------------------------
46        // Types used by this class
47        // ------------------------------------------------------------------
48        /// A Scalar value used by Ipopt
49        typedef Ipopt::Number                         Number;
50        /// An index value used by Ipopt
51        typedef Ipopt::Index                          Index;
52        /// Indexing style used in Ipopt sparsity structure
53        typedef Ipopt::TNLP::IndexStyleEnum           IndexStyleEnum;
54        // ------------------------------------------------------------------
55        // Values directly passed in to constuctor
56        // ------------------------------------------------------------------
57        /// dimension of the range space for f(x).
58        /// The objective is sum_i f_i (x).
59        /// Note that, at this point, there is no advantage having nf_ > 1.
60        const size_t                    nf_;
61        /// dimension of the domain space for f(x) and g(x)
62        const size_t                    nx_;
63        /// dimension of the range space for g(x)
64        const size_t                    ng_;
65        /// initial value for x
66        const Dvector&                  xi_;
67        /// lower limit for x
68        const Dvector&                  xl_;
69        /// upper limit for x
70        const Dvector&                  xu_;
71        /// lower limit for g(x)
72        const Dvector&                  gl_;
73        /// upper limit for g(x)
74        const Dvector&                  gu_;
75        /// object that evaluates f(x) and g(x)
76        FG_eval&                        fg_eval_;
77        /// should operation sequence be retaped for each new x.
78        bool                            retape_;
79        /// Should sparse methods be used to compute Jacobians and Hessians
80        /// with forward mode used for Jacobian.
81        bool                            sparse_forward_;
82        /// Should sparse methods be used to compute Jacobians and Hessians
83        /// with reverse mode used for Jacobian.
84        bool                            sparse_reverse_;
85        /// final results are returned to this structure
86        solve_result<Dvector>&          solution_;
87        // ------------------------------------------------------------------
88        // Values that are initilaized by the constructor
89        // ------------------------------------------------------------------
90        /// AD function object that evaluates x -> [ f(x) , g(x) ]
91        /// If retape is false, this object is initialzed by constructor
92        /// otherwise it is set by cache_new_x each time it is called.
93        CppAD::ADFun<double>            adfun_;
94        /// value of x corresponding to previous new_x
95        Dvector                         x0_;
96        /// value of fg corresponding to previous new_x
97        Dvector                         fg0_;
98        // ----------------------------------------------------------------------
99        // Jacobian information
100        // ----------------------------------------------------------------------
101        /// Sparsity pattern for Jacobian of [f(x), g(x) ].
102        /// If sparse is true, this pattern set by constructor and does not change.
103        /// Otherwise this vector has size zero.
104        CppAD::vector< std::set<size_t> > pattern_jac_;
105        /// Row indices of [f(x), g(x)] for Jacobian of g(x) in row order.
106        /// (Set by constructor and not changed.)
107        CppAD::vector<size_t>           row_jac_;
108        /// Column indices for Jacobian of g(x), same order as row_jac_.
109        /// (Set by constructor and not changed.)
110        CppAD::vector<size_t>           col_jac_;
111        /// col_order_jac_ sorts row_jac_ and col_jac_ in column order.
112        /// (Set by constructor and not changed.)
113        CppAD::vector<size_t>           col_order_jac_;
114        /// Work vector used by SparseJacobian, stored here to avoid recalculation.
115        CppAD::sparse_jacobian_work     work_jac_;
116        // ----------------------------------------------------------------------
117        // Hessian information
118        // ----------------------------------------------------------------------
119        /// Sparsity pattern for Hessian of Lagragian
120        /// \f[ L(x) = \sigma \sum_i f_i (x) + \sum_i \lambda_i  g_i (x) \f]
121        /// If sparse is true, this pattern set by constructor and does not change.
122        /// Otherwise this vector has size zero.
123        CppAD::vector< std::set<size_t> > pattern_hes_;
124        /// Row indices of Hessian lower left triangle in row order.
125        /// (Set by constructor and not changed.)
126        CppAD::vector<size_t>           row_hes_;
127        /// Column indices of Hessian left triangle in same order as row_hes_.
128        /// (Set by constructor and not changed.)
129        CppAD::vector<size_t>           col_hes_;
130        /// Work vector used by SparseJacobian, stored here to avoid recalculation.
131        CppAD::sparse_hessian_work      work_hes_;
132        // ------------------------------------------------------------------
133        // Private member functions
134        // ------------------------------------------------------------------
135        /*!
136        Cache information for a new value of x.
137
138        \param x
139        is the new value for x.
140
141        \par x0_
142        the elements of this vector are set to the new value for x.
143
144        \par fg0_
145        the elements of this vector are set to the new value for [f(x), g(x)]
146
147        \par adfun_
148        If retape is true, the operation sequence for this function
149        is changes to correspond to the argument x.
150        If retape is false, the operation sequence is not changed.
151        The zero order Taylor coefficients for this function are set
152        so they correspond to the argument x. 
153        */
154        void cache_new_x(const Number* x)
155        {       size_t i;
156                if( retape_ )
157                {       // make adfun_, as well as x0_ and fg0_ correspond to this x
158                        ADvector a_x(nx_), a_fg(nf_ + ng_);
159                        for(i = 0; i < nx_; i++)
160                        {       x0_[i] = x[i];
161                                a_x[i] = x[i];
162                        }
163                        CppAD::Independent(a_x);
164                        fg_eval_(a_fg, a_x);
165                        adfun_.Dependent(a_x, a_fg);
166                }
167                else
168                {       // make x0_ and fg0_ correspond to this x
169                        for(i = 0; i < nx_; i++)
170                                x0_[i] = x[i];
171                }
172                fg0_ = adfun_.Forward(0, x0_);
173        }
174public:
175        // ----------------------------------------------------------------------
176        /*!
177        Constructor for the interface between ipopt::solve and Ipopt
178       
179        \param nf
180        dimension of the range space for f(x)
181       
182        \param nx
183        dimension of the domain space for f(x) and g(x).
184       
185        \param ng
186        dimension of the range space for g(x)
187       
188        \param xi
189        initial value of x during the optimization procedure (size nx).
190       
191        \param xl
192        lower limit for x (size nx).
193       
194        \param xu
195        upper limit for x (size nx).
196       
197        \param gl
198        lower limit for g(x) (size ng).
199       
200        \param gu
201        upper limit for g(x) (size ng).
202
203        \param fg_eval
204        function object that evaluations f(x) and g(x) using fg_eval(fg, x)
205
206        \param retape
207        should the operation sequence be retaped for each argument value.
208
209        \param sparse_forward
210        should sparse matrix computations be used for Jacobians and Hessians
211        with forward mode for Jacobian.
212
213        \param sparse_reverse
214        should sparse matrix computations be used for Jacobians and Hessians
215        with reverse mode for Jacobian.
216        (sparse_forward and sparse_reverse cannot both be true).
217       
218        \param solution
219        object where final results are stored.
220        */
221        solve_callback(
222                size_t                 nf              ,
223                size_t                 nx              , 
224                size_t                 ng              ,
225                const Dvector&         xi              ,
226                const Dvector&         xl              ,
227                const Dvector&         xu              ,
228                const Dvector&         gl              ,
229                const Dvector&         gu              ,
230                FG_eval&               fg_eval         ,
231                bool                   retape          ,
232                bool                   sparse_forward  ,
233                bool                   sparse_reverse  ,
234                solve_result<Dvector>& solution ) : 
235        nf_ ( nf ),
236        nx_ ( nx ),
237        ng_ ( ng ),
238        xi_ ( xi ),
239        xl_ ( xl ),
240        xu_ ( xu ),
241        gl_ ( gl ),
242        gu_ ( gu ),
243        fg_eval_ ( fg_eval ),
244        retape_ ( retape ),
245        sparse_forward_ ( sparse_forward ),
246        sparse_reverse_ ( sparse_reverse ),
247        solution_ ( solution )
248        {       CPPAD_ASSERT_UNKNOWN( ! ( sparse_forward_ & sparse_reverse_ ) );
249
250                size_t i, j;
251                size_t nfg = nf_ + ng_;
252
253                // initialize x0_ and fg0_ wih proper dimensions and value nan
254                x0_.resize(nx);
255                fg0_.resize(nfg);
256                for(i = 0; i < nx_; i++)
257                        x0_[i] = CppAD::nan(0.0);
258                for(i = 0; i < nfg; i++)
259                        fg0_[i] = CppAD::nan(0.0);
260
261                if( ! retape_ )
262                {       // make adfun_ correspond to x -> [ f(x), g(x) ]
263                        ADvector a_x(nx_), a_fg(nfg);
264                        for(i = 0; i < nx_; i++)
265                                a_x[i] = xi_[i];
266                        CppAD::Independent(a_x);
267                        fg_eval_(a_fg, a_x);
268                        adfun_.Dependent(a_x, a_fg);
269                        // optimize because we will make repeated use of this tape
270                        adfun_.optimize();
271                }
272                if( sparse_forward_ | sparse_reverse_ ) 
273                {       CPPAD_ASSERT_UNKNOWN( ! retape );       
274                        // -----------------------------------------------------------
275                        // Jacobian
276                        pattern_jac_.resize(nf_ + ng_);
277                        if( nx_ <= nf_ + ng_ )
278                        {       // use forward mode to compute sparsity
279                                CppAD::vector< std::set<size_t> > r(nx_);
280                                for(i = 0; i < nx_; i++)
281                                {       CPPAD_ASSERT_UNKNOWN( r[i].empty() );
282                                        r[i].insert(i);
283                                }
284                                pattern_jac_ = adfun_.ForSparseJac(nx_, r);
285                        }
286                        else
287                        {       // use reverse mode to compute sparsity
288                                size_t m = nf_ + ng_;
289                                CppAD::vector< std::set<size_t> > s(m);
290                                for(i = 0; i < m; i++)
291                                {       CPPAD_ASSERT_UNKNOWN( s[i].empty() );
292                                        s[i].insert(i);
293                                }
294                                pattern_jac_ = adfun_.RevSparseJac(m, s);
295                        }
296                        // Set row and column indices in Jacoian of [f(x), g(x)]
297                        // for Jacobian of g(x). These indices are in row major order.
298                        std::set<size_t>:: const_iterator itr, end;
299                        for(i = nf_; i < nfg; i++)
300                        {       itr = pattern_jac_[i].begin();
301                                end = pattern_jac_[i].end();
302                                while( itr != end )
303                                {       j = *itr++;
304                                        row_jac_.push_back(i);
305                                        col_jac_.push_back(j);
306                                }
307                        }
308                        // -----------------------------------------------------------
309                        // Hessian
310                        pattern_hes_.resize(nx_);
311                        CppAD::vector< std::set<size_t> > r(nx_);
312                        for(i = 0; i < nx_; i++)
313                        {       CPPAD_ASSERT_UNKNOWN( r[i].empty() );
314                                r[i].insert(i);
315                        }
316                        // forward Jacobian sparsity for fg
317                        adfun_.ForSparseJac(nx_, r);
318
319                        // The Lagragian can use any of the components of f(x), g(x)
320                        CppAD::vector< std::set<size_t> > s(1);
321                        CPPAD_ASSERT_UNKNOWN( s[0].empty() );
322                        for(i = 0; i < nf_ + ng_; i++)
323                                s[0].insert(i);
324                        pattern_hes_ = adfun_.RevSparseHes(nx_, s);
325
326                        // Set row and column indices for Lower triangle of Hessian
327                        // of Lagragian.  These indices are in row major order.
328                        for(i = 0; i < nx_; i++)
329                        {       itr = pattern_hes_[i].begin();
330                                end = pattern_hes_[i].end();
331                                while( itr != end )
332                                {       j = *itr++;
333                                        if( j <= i )
334                                        {       row_hes_.push_back(i);
335                                                col_hes_.push_back(j);
336                                        }
337                                }
338                        }
339                }
340                else
341                {       // Set row and column indices in Jacoian of [f(x), g(x)]
342                        // for Jacobian of g(x). These indices are in row major order.
343                        for(i = nf_; i < nfg; i++)
344                        {       for(j = 0; j < nx_; j++)
345                                {       row_jac_.push_back(i);
346                                        col_jac_.push_back(j); 
347                                }
348                        }
349                        // Set row and column indices for lower triangle of Hessian.
350                        // These indices are in row major order.
351                        for(i = 0; i < nx_; i++)
352                        {       for(j = 0; j <= i; j++)
353                                {       row_hes_.push_back(i);
354                                        col_hes_.push_back(j);
355                                }
356                        }
357                }
358
359                // Column order indirect sort of the Jacobian indices
360                col_order_jac_.resize( col_jac_.size() );
361                index_sort( col_jac_, col_order_jac_ );
362        }
363        // -----------------------------------------------------------------------
364        /*!
365        Return dimension information about optimization problem.
366       
367        \param[out] n
368        is set to the value nx_.
369       
370        \param[out] m
371        is set to the value ng_.
372       
373        \param[out] nnz_jac_g
374        is set to ng_ * nx_ (sparsity not yet implemented)
375       
376        \param[out] nnz_h_lag
377        is set to nx_*(nx_+1)/2 (sparsity not yet implemented)
378       
379        \param[out] index_style
380        is set to C_STYLE; i.e., zeoro based indexing is used in the
381        information passed to Ipopt.
382        */
383        virtual bool get_nlp_info(
384                Index&          n            , 
385                Index&          m            , 
386                Index&          nnz_jac_g    ,
387                Index&          nnz_h_lag    , 
388                IndexStyleEnum& index_style  )
389        {
390                n         = static_cast<Index>(nx_);
391                m         = static_cast<Index>(ng_);
392                nnz_jac_g = static_cast<Index>(row_jac_.size());
393                nnz_h_lag = static_cast<Index>(row_hes_.size());
394       
395# ifndef NDEBUG
396                if( ! (sparse_forward_ | sparse_reverse_) )
397                {       size_t nnz = static_cast<size_t>(nnz_jac_g);
398                        CPPAD_ASSERT_UNKNOWN( nnz == ng_ * nx_);
399                        //
400                        nnz = static_cast<size_t>(nnz_h_lag);
401                        CPPAD_ASSERT_UNKNOWN( nnz == (nx_ * (nx_ + 1)) / 2 );
402                }
403# endif
404
405                // use the fortran index style for row/col entries
406                index_style = C_STYLE;
407       
408                return true;
409        }
410        // -----------------------------------------------------------------------
411        /*!
412        Return bound information about optimization problem.
413       
414        \param[in] n
415        is the dimension of the domain space for f(x) and g(x); i.e.,
416        it must be equal to nx_.
417       
418        \param[out] x_l
419        is a vector of size nx_.
420        The input value of its elements does not matter.
421        On output, it is a copy of the lower bound for \f$ x \f$; i.e.,
422        xl_.
423       
424        \param[out] x_u
425        is a vector of size nx_.
426        The input value of its elements does not matter.
427        On output, it is a copy of the upper bound for \f$ x \f$; i.e.,
428        xu_.
429       
430        \param[in] m
431        is the dimension of the range space for g(x). i.e.,
432        it must be equal to ng_.
433       
434        \param[out] g_l
435        is a vector of size ng_.
436        The input value of its elements does not matter.
437        On output, it is a copy of the lower bound for \f$ g(x) \f$; i.e., gl_.
438       
439        \param[out] g_u
440        is a vector of size ng_.
441        The input value of its elements does not matter.
442        On output, it is a copy of the upper bound for \f$ g(x) \f$; i.e, gu_.
443        */
444        virtual bool get_bounds_info(
445                Index       n        , 
446                Number*     x_l      , 
447                Number*     x_u      ,
448                Index       m        , 
449                Number*     g_l      , 
450                Number*     g_u      )
451        {       size_t i;
452                // here, the n and m we gave IPOPT in get_nlp_info are passed back
453                CPPAD_ASSERT_UNKNOWN(static_cast<size_t>(n) == nx_);
454                CPPAD_ASSERT_UNKNOWN(static_cast<size_t>(m) == ng_);
455       
456                // pass back bounds
457                for(i = 0; i < nx_; i++)
458                {       x_l[i] = xl_[i];
459                        x_u[i] = xu_[i];
460                }
461                for(i = 0; i < ng_; i++)
462                {       g_l[i] = gl_[i];
463                        g_u[i] = gu_[i];
464                }
465               
466                return true;
467        }
468        // -----------------------------------------------------------------------
469        /*!
470        Return initial x value where optimiation is started.
471       
472        \param[in] n
473        must be equal to the domain dimension for f(x) and g(x); i.e.,
474        it must be equal to nx_.
475       
476        \param[in] init_x
477        must be equal to true.
478       
479        \param[out] x
480        is a vector of size nx_.
481        The input value of its elements does not matter.
482        On output, it is a copy of the initial value for \f$ x \f$; i.e. xi_.
483       
484        \param[in] init_z
485        must be equal to false.
486       
487        \param z_L
488        is not used.
489       
490        \param z_U
491        is not used.
492       
493        \param[in] m
494        must be equal to the domain dimension for f(x) and g(x); i.e.,
495        it must be equal to ng_.
496       
497        \param init_lambda
498        must be equal to false.
499       
500        \param lambda
501        is not used.
502        */
503        virtual bool get_starting_point(
504                Index           n            , 
505                bool            init_x       , 
506                Number*         x            ,
507                bool            init_z       , 
508                Number*         z_L          , 
509                Number*         z_U          ,
510                Index           m            , 
511                bool            init_lambda  ,
512                Number*         lambda       )
513        {       size_t j;
514       
515                CPPAD_ASSERT_UNKNOWN(static_cast<size_t>(n) == nx_ );
516                CPPAD_ASSERT_UNKNOWN(static_cast<size_t>(m) == ng_ );
517                CPPAD_ASSERT_UNKNOWN(init_x == true);
518                CPPAD_ASSERT_UNKNOWN(init_z == false);
519                CPPAD_ASSERT_UNKNOWN(init_lambda == false);
520       
521                for(j = 0; j < nx_; j++)
522                        x[j] = xi_[j];
523       
524                return true;
525        }
526        // -----------------------------------------------------------------------
527        /*!
528        Evaluate the objective fucntion f(x).
529       
530        \param[in] n
531        is the dimension of the argument space for f(x); i.e., must be equal nx_.
532       
533        \param[in] x
534        is a vector of size nx_ containing the point at which to evaluate
535        the function sum_i f_i (x).
536       
537        \param[in] new_x
538        is false if the previous call to any one of the
539        \ref Evaluation_Methods used the same value for x.
540       
541        \param[out] obj_value
542        is the value of the objective sum_i f_i (x) at this value of x.
543       
544        \return
545        The return value is always true; see \ref Evaluation_Methods.
546        */
547        virtual bool eval_f(
548                Index          n           , 
549                const Number*  x           , 
550                bool           new_x       , 
551                Number&        obj_value   )
552        {       size_t i;
553                if( new_x )
554                        cache_new_x(x);
555                //
556                double sum = 0.0;
557                for(i = 0; i < nf_; i++)
558                        sum += fg0_[i];
559                obj_value = static_cast<Number>(sum);
560                return true;
561        }
562        // -----------------------------------------------------------------------
563        /*!
564        Evaluate the gradient of f(x).
565       
566        \param[in] n
567        is the dimension of the argument space for f(x); i.e., must be equal nx_.
568       
569        \param[in] x
570        has a vector of size nx_ containing the point at which to evaluate
571        the gradient of f(x).
572       
573        \param[in] new_x
574        is false if the previous call to any one of the
575        \ref Evaluation_Methods used the same value for x.
576       
577        \param[out] grad_f
578        is a vector of size nx_.
579        The input value of its elements does not matter.
580        The output value of its elements is the gradient of f(x)
581        at this value of.
582       
583        \return
584        The return value is always true; see \ref Evaluation_Methods.
585        */
586        virtual bool eval_grad_f(
587                Index           n        , 
588                const Number*   x        , 
589                bool            new_x    , 
590                Number*         grad_f   )
591        {       size_t i;
592                if( new_x )
593                        cache_new_x(x);
594                //
595                Dvector w(nf_ + ng_), dw(nx_);
596                for(i = 0; i < nf_; i++)
597                        w[i] = 1.0;
598                for(i = 0; i < ng_; i++)
599                        w[nf_ + i] = 0.0;
600                dw = adfun_.Reverse(1, w);
601                for(i = 0; i < nx_; i++)
602                        grad_f[i] = dw[i];
603                return true;
604        }
605        // -----------------------------------------------------------------------
606        /*!
607        Evaluate the function g(x).
608       
609        \param[in] n
610        is the dimension of the argument space for g(x); i.e., must be equal nx_.
611       
612        \param[in] x
613        has a vector of size n containing the point at which to evaluate
614        the gradient of g(x).
615       
616        \param[in] new_x
617        is false if the previous call to any one of the
618        \ref Evaluation_Methods used the same value for x.
619       
620        \param[in] m
621        is the dimension of the range space for g(x); i.e., must be equal to ng_.
622       
623        \param[out] g
624        is a vector of size ng_.
625        The input value of its elements does not matter.
626        The output value of its elements is
627        the value of the function g(x) at this value of x.
628       
629        \return
630        The return value is always true; see \ref Evaluation_Methods.
631        */
632        virtual bool eval_g(
633                Index   n            , 
634                const   Number* x    , 
635                bool    new_x        , 
636                Index   m            , 
637                Number* g            )
638        {       size_t i;
639                if( new_x )
640                        cache_new_x(x);
641                //
642                for(i = 0; i < ng_; i++)
643                        g[i] = fg0_[nf_ + i];
644                return true;
645        }
646        // -----------------------------------------------------------------------
647        /*!
648        Evaluate the Jacobian of g(x).
649       
650        \param[in] n
651        is the dimension of the argument space for g(x);
652        i.e., must be equal nx_.
653       
654        \param x
655        If values is not NULL,
656        x is a vector of size nx_ containing the point at which to evaluate
657        the gradient of g(x).
658       
659        \param[in] new_x
660        is false if the previous call to any one of the
661        \ref Evaluation_Methods used the same value for  x.
662       
663        \param[in] m
664        is the dimension of the range space for g(x);
665        i.e., must be equal to ng_.
666       
667        \param[in] nele_jac
668        is the number of possibly non-zero elements in the Jacobian of g(x);
669        i.e., must be equal to ng_ * nx_.
670       
671        \param iRow
672        if values is not NULL, iRow is not defined.
673        if values is NULL, iRow
674        is a vector with size nele_jac.
675        The input value of its elements does not matter.
676        On output,
677        For <tt>k = 0 , ... , nele_jac-1, iRow[k]</tt> is the
678        base zero row index for the
679        k-th possibly non-zero entry in the Jacobian of g(x).
680       
681        \param jCol
682        if values is not NULL, jCol is not defined.
683        if values is NULL, jCol
684        is a vector with size nele_jac.
685        The input value of its elements does not matter.
686        On output,
687        For <tt>k = 0 , ... , nele_jac-1, jCol[k]</tt> is the
688        base zero column index for the
689        k-th possibly non-zero entry in the Jacobian of g(x).
690       
691        \param values
692        if \c values is not \c NULL, \c values
693        is a vector with size \c nele_jac.
694        The input value of its elements does not matter.
695        On output,
696        For <tt>k = 0 , ... , nele_jac-1, values[k]</tt> is the
697        value for the
698        k-th possibly non-zero entry in the Jacobian of g(x).
699       
700        \return
701        The return value is always true; see \ref Evaluation_Methods.
702        */
703        virtual bool eval_jac_g(
704                Index n, 
705                const Number* x, 
706                bool new_x,
707                       
708                Index m, 
709                Index nele_jac, 
710                Index* iRow, 
711                Index *jCol,
712                       
713                Number* values)
714        {       size_t i, j, k, ell;
715                CPPAD_ASSERT_UNKNOWN(static_cast<size_t>(m)         == ng_ );
716                CPPAD_ASSERT_UNKNOWN(static_cast<size_t>(n)         == nx_ );
717                //
718                size_t nk = row_jac_.size();
719                CPPAD_ASSERT_UNKNOWN( static_cast<size_t>(nele_jac) == nk );
720                //
721                if( new_x )
722                        cache_new_x(x);
723
724                if( values == NULL )
725                {       for(k = 0; k < nk; k++)
726                        {       i = row_jac_[k];
727                                j = col_jac_[k];
728                                CPPAD_ASSERT_UNKNOWN( i >= nf_ );
729                                iRow[k] = static_cast<Index>(i - nf_);
730                                jCol[k] = static_cast<Index>(j);
731                        }
732                        return true;
733                }
734                //
735                if( nk == 0 )
736                        return true;
737                //
738                if( sparse_forward_ )
739                {       Dvector jac(nk);
740                        adfun_.SparseJacobianForward(
741                                x0_ , pattern_jac_, row_jac_, col_jac_, jac, work_jac_
742                        );
743                        for(k = 0; k < nk; k++)
744                                values[k] = jac[k];
745                }
746                else if( sparse_reverse_ )
747                {       Dvector jac(nk);
748                        adfun_.SparseJacobianReverse(
749                                x0_ , pattern_jac_, row_jac_, col_jac_, jac, work_jac_
750                        );
751                        for(k = 0; k < nk; k++)
752                                values[k] = jac[k];
753                }
754                else if( nx_ < ng_ )
755                {       // use forward mode
756                        Dvector x1(nx_), fg1(nf_ + ng_);
757                        for(j = 0; j < nx_; j++)
758                                x1[j] = 0.0;
759                        // index in col_order_jac_ of next entry
760                        ell = 0;
761                        k   = col_order_jac_[ell];
762                        for(j = 0; j < nx_; j++)
763                        {       // compute j-th column of Jacobian of g(x)
764                                x1[j] = 1.0;
765                                fg1 = adfun_.Forward(1, x1);
766                                while( ell < nk && col_jac_[k] <= j )
767                                {       CPPAD_ASSERT_UNKNOWN( col_jac_[k] == j );
768                                        i = row_jac_[k];
769                                        CPPAD_ASSERT_UNKNOWN( i >= nf_ )
770                                        values[k] = fg1[i];
771                                        ell++;
772                                        if( ell < nk )
773                                                k = col_order_jac_[ell];
774                                }
775                                x1[j] = 0.0;
776                        }
777                }
778                else
779                {       // user reverse mode   
780                        size_t nfg = nf_ + ng_;
781                        // user reverse mode
782                        Dvector w(nfg), dw(nx_);
783                        for(i = 0; i < nfg; i++)
784                                w[i] = 0.0;
785                        // index in row_jac_ of next entry
786                        k = 0;
787                        for(i = nf_; i < nfg; i++)
788                        {       // compute i-th row of Jacobian of g(x)
789                                w[i] = 1.0;
790                                dw = adfun_.Reverse(1, w);
791                                while( k < nk && row_jac_[k] <= i )
792                                {       CPPAD_ASSERT_UNKNOWN( row_jac_[k] == i );
793                                        j = col_jac_[k];
794                                        values[k] = dw[j];
795                                        k++;
796                                }
797                                w[i] = 0.0;
798                        }
799                }
800                return true;
801        }
802        // -----------------------------------------------------------------------
803        /*!
804        Evaluate the Hessian of the Lagragian
805       
806        \section The_Hessian_of_the_Lagragian The Hessian of the Lagragian
807        The Hessian of the Lagragian is defined as
808        \f[
809        H(x, \sigma, \lambda )
810        =
811        \sigma \nabla^2 f(x) + \sum_{i=0}^{m-1} \lambda_i \nabla^2 g(x)_i
812        \f]
813       
814        \param[in] n
815        is the dimension of the argument space for g(x);
816        i.e., must be equal nx_.
817       
818        \param x
819        if values is not NULL, x
820        is a vector of size nx_ containing the point at which to evaluate
821        the Hessian of the Lagragian.
822       
823        \param[in] new_x
824        is false if the previous call to any one of the
825        \ref Evaluation_Methods used the same value for x.
826       
827        \param[in] obj_factor
828        the value \f$ \sigma \f$ multiplying the Hessian of
829        f(x) in the expression for \ref The_Hessian_of_the_Lagragian.
830       
831        \param[in] m
832        is the dimension of the range space for g(x);
833        i.e., must be equal to ng_.
834       
835        \param[in] lambda
836        if values is not NULL, lambda
837        is a vector of size ng_ specifing the value of \f$ \lambda \f$
838        in the expression for \ref The_Hessian_of_the_Lagragian.
839       
840        \param[in] new_lambda
841        is true if the previous call to eval_h had the same value for
842        lambda and false otherwise.
843        (Not currently used.)
844       
845        \param[in] nele_hess
846        is the number of possibly non-zero elements in the
847        Hessian of the Lagragian;
848        i.e., must be equal to nx_*(nx_+1)/2.
849       
850        \param iRow
851        if values is not NULL, iRow is not defined.
852        if values is NULL, iRow
853        is a vector with size nele_hess.
854        The input value of its elements does not matter.
855        On output,
856        For <tt>k = 0 , ... , nele_hess-1, iRow[k]</tt> is the
857        base zero row index for the
858        k-th possibly non-zero entry in the Hessian fo the Lagragian.
859       
860        \param jCol
861        if values is not NULL, jCol is not defined.
862        if values is NULL, jCol
863        is a vector with size nele_hess.
864        The input value of its elements does not matter.
865        On output,
866        For <tt>k = 0 , ... , nele_hess-1, jCol[k]</tt> is the
867        base zero column index for the
868        k-th possibly non-zero entry in the Hessian of the Lagragian.
869       
870        \param values
871        if values is not NULL, it
872        is a vector with size nele_hess.
873        The input value of its elements does not matter.
874        On output,
875        For <tt>k = 0 , ... , nele_hess-1, values[k]</tt> is the
876        value for the
877        k-th possibly non-zero entry in the Hessian of the Lagragian.
878       
879        \return
880        The return value is always true; see \ref Evaluation_Methods.
881        */
882        virtual bool eval_h(
883                Index         n              ,
884                const Number* x              , 
885                bool          new_x          ,
886                Number        obj_factor     , 
887                Index         m              ,
888                const Number* lambda         ,
889                bool          new_lambda     , 
890                Index         nele_hess      , 
891                Index*        iRow           ,
892                Index*        jCol           , 
893                Number*       values         )
894        {       size_t i, j, k;
895                CPPAD_ASSERT_UNKNOWN(static_cast<size_t>(m) == ng_ );
896                CPPAD_ASSERT_UNKNOWN(static_cast<size_t>(n) == nx_ );
897                //
898                size_t nk = row_hes_.size();
899                CPPAD_ASSERT_UNKNOWN( static_cast<size_t>(nele_hess) == nk ); 
900                //
901                if( new_x )
902                        cache_new_x(x);
903                //
904                if( values == NULL )
905                {       for(k = 0; k < nk; k++)
906                        {       i = row_hes_[k];
907                                j = col_hes_[k];
908                                iRow[k] = static_cast<Index>(i);
909                                jCol[k] = static_cast<Index>(j);
910                        }
911                        return true;
912                }
913                //
914                if( nk == 0 )
915                        return true;
916
917                // weigting vector for Lagragian
918                Dvector w(nf_ + ng_);
919                for(i = 0; i < nf_; i++)
920                        w[i] = obj_factor;
921                for(i = 0; i < ng_; i++)
922                        w[i + nf_] = lambda[i];
923                //
924                if( sparse_forward_ | sparse_reverse_ )
925                {       Dvector hes(nk);
926                        adfun_.SparseHessian(
927                                x0_, w, pattern_hes_, row_hes_, col_hes_, hes, work_hes_
928                        );
929                        for(k = 0; k < nk; k++)
930                                values[k] = hes[k];
931                }
932                else
933                {       Dvector hes(nx_ * nx_);
934                        hes = adfun_.Hessian(x0_, w);
935                        for(k = 0; k < nk; k++)
936                        {       i = row_hes_[k];
937                                j = col_hes_[k];
938                                values[k] = hes[i * nx_ + j];
939                        }
940                }
941                return true;
942        }
943        // ----------------------------------------------------------------------
944        /*!
945        Pass solution information from Ipopt to users solution structure.
946       
947        \param[in] status
948        is value that the Ipopt solution status
949        which gets mapped to a correponding value for
950        \n
951        <tt>solution_.status</tt>
952       
953        \param[in] n
954        is the dimension of the domain space for f(x) and g(x); i.e.,
955        it must be equal to nx_.
956       
957        \param[in] x
958        is a vector with size nx_ specifing the final solution.
959        This is the output value for
960        \n
961        <tt>solution_.x</tt>
962       
963        \param[in] z_L
964        is a vector with size nx_ specifing the Lagragian multipliers for the
965        constraint \f$ x^l \leq x \f$.
966        This is the output value for
967        \n
968        <tt>solution_.zl</tt>
969       
970        \param[in] z_U
971        is a vector with size nx_ specifing the Lagragian multipliers for the
972        constraint \f$ x \leq x^u \f$.
973        This is the output value for
974        \n
975        <tt>solution_.zu</tt>
976       
977        \param[in] m
978        is the dimension of the range space for g(x). i.e.,
979        it must be equal to ng_.
980       
981        \param[in] g
982        is a vector with size ng_ containing the value of the constraint function
983        g(x) at the final solution x.
984        This is the output value for
985        \n
986        <tt>solution_.g</tt>
987       
988        \param[in] lambda
989        is a vector with size ng_ specifing the Lagragian multipliers for the
990        constraints \f$ g^l \leq g(x) \leq g^u \f$.
991        This is the output value for
992        \n
993        <tt>solution_.lambda</tt>
994       
995        \param[in] obj_value
996        is the value of the objective function f(x) at the final solution x.
997        This is the output value for
998        \n
999        <tt>solution_.obj_value</tt>
1000       
1001        \param[in] ip_data
1002        is unspecified (by Ipopt) and hence not used.
1003       
1004        \param[in] ip_cq
1005        is unspecified (by Ipopt) and hence not used.
1006       
1007        \par solution_[out]
1008        this is a reference to the solution argument
1009        in the constructor for solve_callback.
1010        The results are stored here
1011        (see documentation above).
1012        */
1013        virtual void finalize_solution(
1014                Ipopt::SolverReturn               status    ,
1015                Index                             n         , 
1016                const Number*                     x         , 
1017                const Number*                     z_L       , 
1018                const Number*                     z_U       ,
1019                Index                             m         , 
1020                const Number*                     g         , 
1021                const Number*                     lambda    ,
1022                Number                            obj_value ,
1023                const Ipopt::IpoptData*           ip_data   ,
1024                Ipopt::IpoptCalculatedQuantities* ip_cq
1025        )
1026        {       size_t i, j;
1027       
1028                CPPAD_ASSERT_UNKNOWN(static_cast<size_t>(n) == nx_ );
1029                CPPAD_ASSERT_UNKNOWN(static_cast<size_t>(m) == ng_ );
1030       
1031                switch(status)
1032                {       // convert status from Ipopt enum to solve_result<Dvector> enum
1033                        case Ipopt::SUCCESS:
1034                        solution_.status = solve_result<Dvector>::success;
1035                        break;
1036       
1037                        case Ipopt::MAXITER_EXCEEDED:
1038                        solution_.status = 
1039                                solve_result<Dvector>::maxiter_exceeded;
1040                        break;
1041       
1042                        case Ipopt::STOP_AT_TINY_STEP:
1043                        solution_.status = 
1044                                solve_result<Dvector>::stop_at_tiny_step;
1045                        break;
1046       
1047                        case Ipopt::STOP_AT_ACCEPTABLE_POINT:
1048                        solution_.status = 
1049                                solve_result<Dvector>::stop_at_acceptable_point;
1050                        break;
1051       
1052                        case Ipopt::LOCAL_INFEASIBILITY:
1053                        solution_.status = 
1054                                solve_result<Dvector>::local_infeasibility;
1055                        break;
1056       
1057                        case Ipopt::USER_REQUESTED_STOP:
1058                        solution_.status = 
1059                                solve_result<Dvector>::user_requested_stop;
1060                        break;
1061       
1062                        case Ipopt::DIVERGING_ITERATES:
1063                        solution_.status = 
1064                                solve_result<Dvector>::diverging_iterates;
1065                        break;
1066       
1067                        case Ipopt::RESTORATION_FAILURE:
1068                        solution_.status = 
1069                                solve_result<Dvector>::restoration_failure;
1070                        break;
1071       
1072                        case Ipopt::ERROR_IN_STEP_COMPUTATION:
1073                        solution_.status = 
1074                                solve_result<Dvector>::error_in_step_computation;
1075                        break;
1076       
1077                        case Ipopt::INVALID_NUMBER_DETECTED:
1078                        solution_.status = 
1079                                solve_result<Dvector>::invalid_number_detected;
1080                        break;
1081       
1082                        case Ipopt::INTERNAL_ERROR:
1083                        solution_.status = 
1084                                solve_result<Dvector>::internal_error;
1085                        break;
1086       
1087                        default:
1088                        solution_.status = 
1089                                solve_result<Dvector>::unknown;
1090                }
1091       
1092                solution_.x.resize(nx_);
1093                solution_.zl.resize(nx_);
1094                solution_.zu.resize(nx_);
1095                for(j = 0; j < nx_; j++)
1096                {       solution_.x[j]   = x[j];
1097                        solution_.zl[j]  = z_L[j];
1098                        solution_.zu[j]  = z_U[j];
1099                }
1100                solution_.g.resize(ng_);
1101                solution_.lambda.resize(ng_);
1102                for(i = 0; i < ng_; i++)
1103                {       solution_.g[i]      = g[i];
1104                        solution_.lambda[i] = lambda[i];
1105                }
1106                solution_.obj_value = obj_value;
1107                return;
1108        }
1109};
1110
1111} // end namespace ipopt
1112} // END_CPPAD_NAMESPACE
1113
1114# endif
Note: See TracBrowser for help on using the repository browser.