gsOptProblem.h

Go to the documentation of this file.

 
#include <gsCore/gsLinearAlgebra.h>
 
#pragma once
 
namespace gismo
{
template <typename T>
class gsOptProblem
{
 
public:
 
    // /** default constructor */
    // gsOptProblem();
 
    // /** default destructor */
    // virtual ~gsOptProblem();
 
public:
 
    virtual T evalObj( const gsAsConstVector<T> &) const
    {GISMO_NO_IMPLEMENTATION }
 
    virtual void gradObj_into ( const gsAsConstVector<T> & u, gsAsVector<T> & result) const
    {
        const index_t n = u.rows();
        //GISMO_ASSERT((index_t)m_numDesignVars == n*m, "Wrong design.");
 
        gsMatrix<T> uu = u;//copy
        gsAsVector<T> tmp(uu.data(), n);
        gsAsConstVector<T> ctmp(uu.data(), n);
        index_t c = 0;
 
        // for all partial derivatives (column-wise)
        for ( index_t i = 0; i!=n; i++ )
        {
            // to do: add m_desLowerBounds m_desUpperBounds check
            tmp[i]  += T(0.00001);
            const T e1 = this->evalObj(ctmp);
            tmp[i]   = u[i] + T(0.00002);
            const T e3 = this->evalObj(ctmp);
            tmp[i]   = u[i] - T(0.00001);
            const T e2 = this->evalObj(ctmp);
            tmp[i]   = u[i] - T(0.00002);
            const T e4 = this->evalObj(ctmp);
            tmp[i]   = u[i];
            result[c++]= ( 8 * (e1 - e2) + e4 - e3 ) / T(0.00012);
        }
    }
 
 
    virtual void evalCon_into ( const gsAsConstVector<T> &, gsAsVector<T> &) const
    {GISMO_NO_IMPLEMENTATION }
 
    virtual void jacobCon_into( const gsAsConstVector<T> &, gsAsVector<T> &) const
    {GISMO_NO_IMPLEMENTATION }
 
    virtual void hessLagr_into( const gsAsConstVector<T> &, gsAsVector<T> &) const
    {GISMO_NO_IMPLEMENTATION }
 
    virtual void computeJacStructure()
    {
        // full Jacobian
        m_conJacRows.resize(m_numConstraints*m_numDesignVars);
        m_conJacCols.resize(m_numConstraints*m_numDesignVars);
        index_t c = 0;
        for ( index_t j = 0; j < m_numDesignVars; ++j)
            for ( index_t i = 0; i < m_numConstraints; ++i)
            {
                m_conJacRows[c  ] = i;
                m_conJacCols[c++] = j;
            }
 
        m_numConJacNonZero = m_conJacRows.size();
    }
 
    //void evalSingleCon_into(int k, gsMatrix<T> & result);
    //void hassObj_into( gsMatrix<T> & result);
 
 
    // virtual bool intermediateCallback() { return true;}
 
public:
 
    int numDesignVars () const { return m_curDesign.size(); }
 
    int numConstraints() const { return m_numConstraints; }
 
    int numConJacNonZero() const { return m_numConJacNonZero; }
 
    int numConDerivs  () const { return m_conJacRows.size(); }
 
    T lowerCon(int k) const { return m_conLowerBounds[k]; }
 
    T upperCon(int k) const { return m_conUpperBounds[k]; }
 
    T lowerDesignVar(int k) const { return m_desLowerBounds[k]; }
 
    T upperDesignVar(int k) const { return m_desUpperBounds[k]; }
 
    const gsVector<T> & desLowerBounds() const {return m_desLowerBounds;}
    const gsVector<T> & desUpperBounds() const {return m_desUpperBounds;}
 
    const gsVector<T> & conLowerBounds() const {return m_conLowerBounds;}
    const gsVector<T> & conUpperBounds() const {return m_conUpperBounds;}
 
    const std::vector<index_t> & conJacRows() const { return m_conJacRows; }
 
    const std::vector<index_t> & conJacCols() const { return m_conJacCols; }
 
public:
 
    // void solve ();
 
    // std::ostream &print(std::ostream &os) const
    // {
    //     os << "Design variables: " << m_numDesignVars
    //        << "\nNumber of constraints: " << m_numConstraints
    //         //<< "\ndesign lower:" << m_desLowerBounds.transpose()
    //         //<< "\ndesign upper:" << m_desUpperBounds.transpose()
    //         //<< "\nconstr lower:" << m_conLowerBounds.transpose()
    //         //<< "\nconstr upper:" << m_conUpperBounds.transpose()
    //        << "\nNumber of active sensitivities: " << m_numConJacNonZero
    //        << "\nSparsity of sensitivities: " <<std::fixed<<std::setprecision(2)
    //        << (100.0*m_numConJacNonZero)/(m_numDesignVars*m_numConstraints) << " %"
    //         //<< "\nm_conJacRows:" << m_conJacRows.transpose()
    //         //<< "\nm_conJacCols:" << m_conJacCols.transpose()
    //         //<< "\nm_curDesign:" << m_curDesign.transpose()
    //         ;
    //     return os;
    // }
 
 
protected:
 
    int m_numDesignVars;
 
    int m_numConstraints;
 
    int m_numConJacNonZero;
 
    gsVector<T> m_desLowerBounds;
 
    gsVector<T> m_desUpperBounds;
 
    gsVector<T> m_conLowerBounds;
 
    gsVector<T> m_conUpperBounds;
 
    std::vector<index_t> m_conJacRows;
 
    std::vector<index_t> m_conJacCols;
 
    gsMatrix<T> m_curDesign;
 
    gsMatrix<T> m_lambda;
};
 
 
} // end namespace gismo
 
// note: must be statically compiled in header-only mode
#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsOptProblem.hpp)
#endif