InteriorPointBlock.hxx

Go to the documentation of this file.

#ifndef CONICBUNDLE_INTERIORPOINTBLOCK_HXX
#define CONICBUNDLE_INTERIORPOINTBLOCK_HXX

#include "symmat.hxx"
#include "CBout.hxx"

namespace ConicBundle {
  

class InteriorPointBlock: public virtual CBout
{
protected:
  int pol_le_zero_step(CH_Matrix_Classes::Real& stepsize,
                       CH_Matrix_Classes::Real q0,
                       CH_Matrix_Classes::Real q1,
                       CH_Matrix_Classes::Real q2,
                       CH_Matrix_Classes::Real q3,
                       CH_Matrix_Classes::Real q4,
                       CH_Matrix_Classes::Real abseps=1e-10) const;

  int minimize_pol_step(CH_Matrix_Classes::Real& stepsize,
                        CH_Matrix_Classes::Real q0,
                        CH_Matrix_Classes::Real q1,
                        CH_Matrix_Classes::Real q2,
                        CH_Matrix_Classes::Real q3,
                        CH_Matrix_Classes::Real q4,
                        CH_Matrix_Classes::Real abseps=1e-10) const;

  int control_nbh_step(CH_Matrix_Classes::Real& stepsize,
                       CH_Matrix_Classes::Real& max_nbh,
                       CH_Matrix_Classes::Real nbh_ubnd,
                       CH_Matrix_Classes::Real mu_xz,
                       CH_Matrix_Classes::Real mu_xdzpdxz,
                       CH_Matrix_Classes::Real mu_dxdz,
                       CH_Matrix_Classes::Real nrmsqr_xz,
                       CH_Matrix_Classes::Real nrmsqr_xdzpdxz,
                       CH_Matrix_Classes::Real nrmsqr_dxdz,
                       CH_Matrix_Classes::Real ip_xz_xdzpdxz,
                       CH_Matrix_Classes::Real ip_xz_dxdz,
                       CH_Matrix_Classes::Real ip_dxdz_xdzpdxz) const;

public:
  InteriorPointBlock(CBout* cb=0,int cbinc=-1):CBout(cb,cbinc){}

  virtual ~InteriorPointBlock(); 

  virtual CH_Matrix_Classes::Integer get_vecdim() const =0;
  
  virtual int center_x(CH_Matrix_Classes::Real val,bool add=false)=0;

  virtual int center_z(CH_Matrix_Classes::Real val,bool add=false)=0;

  virtual int set_x(const CH_Matrix_Classes::Matrix& vec,
                    CH_Matrix_Classes::Integer startindex,
                    CH_Matrix_Classes::Real& add_center_value)=0;

  virtual int set_z(const CH_Matrix_Classes::Matrix& vec,
                    CH_Matrix_Classes::Integer startindex,
                    CH_Matrix_Classes::Real& add_center_value)=0;

  virtual int vecgetsax(CH_Matrix_Classes::Matrix& vec,
                        CH_Matrix_Classes::Integer startindex,
                        CH_Matrix_Classes::Real a=1.,
                        bool add=false)=0;

  virtual int vecgetsaz(CH_Matrix_Classes::Matrix& vec,
                        CH_Matrix_Classes::Integer startindex,
                        CH_Matrix_Classes::Real a=1.,
                        bool add=false)=0;

  virtual int get_mu_info(CH_Matrix_Classes::Integer& mudim,
                          CH_Matrix_Classes::Real& tr_xz,
                          CH_Matrix_Classes::Real& tr_xdzpdxz,
                          CH_Matrix_Classes::Real& tr_dxdz,
                          CH_Matrix_Classes::Real& min_xz,
                          CH_Matrix_Classes::Real& max_xz) const=0;

  virtual int get_nbh_info(CH_Matrix_Classes::Integer mudim,
                           CH_Matrix_Classes::Real tr_xz,
                           CH_Matrix_Classes::Real tr_xdzpdxz,
                           CH_Matrix_Classes::Real tr_dxdz,
                           CH_Matrix_Classes::Real nbh_ubnd,
                           CH_Matrix_Classes::Real& alpha,
                           CH_Matrix_Classes::Real& max_nbh,
                           CH_Matrix_Classes::Real& nrmsqr_xz,
                           CH_Matrix_Classes::Real& nrmsqr_xdzpdxz,
                           CH_Matrix_Classes::Real& nrmsqr_dxdz,
                           CH_Matrix_Classes::Real& ip_xz_xdzpdxz,
                           CH_Matrix_Classes::Real& ip_xz_dxdz,
                           CH_Matrix_Classes::Real& ip_dxdz_xdzpdxz) const=0;


  virtual int linesearch(CH_Matrix_Classes::Real& alpha) const=0;

  virtual int add_muxinv(CH_Matrix_Classes::Matrix& rhs,
                         CH_Matrix_Classes::Integer startindex,
                         CH_Matrix_Classes::Real rhsmu,
                         CH_Matrix_Classes::Real rhscorr,
                         bool minus=false)=0;

  virtual int set_dx(const CH_Matrix_Classes::Matrix& rhs,CH_Matrix_Classes::Integer startindex)=0;

  virtual int set_dx_xizsolverhs(const CH_Matrix_Classes::Matrix& rhs,
                                CH_Matrix_Classes::Integer startindex)=0;

  virtual int apply_xizinv(CH_Matrix_Classes::Matrix& rhs,
                           CH_Matrix_Classes::Integer startindex,
                           bool minus=false)=0;

  virtual int apply_xiz(CH_Matrix_Classes::Matrix& rhs,
                        CH_Matrix_Classes::Integer startindex,
                        bool minus=false)=0;

  virtual int do_step(CH_Matrix_Classes::Real alpha)=0;


  virtual int add_AxizinvAt(const CH_Matrix_Classes::Matrix& A,
                            CH_Matrix_Classes::Symmatrix& globalsys,
                            bool minus=false,
                            bool Atrans=false)=0;

  virtual int add_xiz(CH_Matrix_Classes::Symmatrix& globalsys,
                      CH_Matrix_Classes::Integer startindex,
                      bool minus=false)=0;

  //---------------------- mainly for testing

  virtual int get_vecx(CH_Matrix_Classes::Matrix& vecx,CH_Matrix_Classes::Integer startindex)=0;
  
  virtual int get_vecz(CH_Matrix_Classes::Matrix& vecz,CH_Matrix_Classes::Integer startindex)=0;
  
  virtual int get_vecdx(CH_Matrix_Classes::Matrix& vecdx,CH_Matrix_Classes::Integer startindex)=0;
  
  virtual int get_vecdz(CH_Matrix_Classes::Matrix& vecdz,CH_Matrix_Classes::Integer startindex)=0;
  

};

  //--------------------------------------------------------
  inline void NNC_nbh_stepsize(CH_Matrix_Classes::Real x,
                               CH_Matrix_Classes::Real z,
                               CH_Matrix_Classes::Real dx,
                               CH_Matrix_Classes::Real dz,
                               CH_Matrix_Classes::Real mu_xz,
                               CH_Matrix_Classes::Real mu_xdzpdxz,
                               CH_Matrix_Classes::Real mu_dxdz,
                               CH_Matrix_Classes::Real mu_at_one,
                               CH_Matrix_Classes::Real nbh_ubnd,
                               CH_Matrix_Classes::Real& alpha,
                               CH_Matrix_Classes::Real& max_nbh,
                               CH_Matrix_Classes::Real& nrmsqr_xz,
                               CH_Matrix_Classes::Real& nrmsqr_xdzpdxz,
                               CH_Matrix_Classes::Real& nrmsqr_dxdz,
                               CH_Matrix_Classes::Real& ip_xz_xdzpdxz,
                               CH_Matrix_Classes::Real& ip_xz_dxdz,
                               CH_Matrix_Classes::Real& ip_dxdz_xdzpdxz)
  {
    CH_Matrix_Classes::Real xz=x*z-mu_xz;
    CH_Matrix_Classes::Real xdzpdxz=x*dz+dx*z-mu_xdzpdxz;
    CH_Matrix_Classes::Real dxdz=dx*dz-mu_dxdz;
    nrmsqr_xz+=CH_Matrix_Classes::sqr(xz);
    nrmsqr_xdzpdxz+=CH_Matrix_Classes::sqr(xdzpdxz);
    nrmsqr_dxdz+=CH_Matrix_Classes::sqr(dxdz);
    ip_xz_xdzpdxz+=xz*xdzpdxz;
    ip_xz_dxdz+=xz*dxdz;
    ip_dxdz_xdzpdxz+=dxdz*xdzpdxz;
    
    if (x+alpha*dx<CH_Matrix_Classes::eps_Real*mu_xz)
      alpha=CH_Matrix_Classes::min(alpha,CH_Matrix_Classes::max(0.,(CH_Matrix_Classes::eps_Real*mu_xz-x)/dx));
    if (z+alpha*dz<CH_Matrix_Classes::eps_Real*mu_xz)
      alpha=CH_Matrix_Classes::min(alpha,CH_Matrix_Classes::max(0.,(CH_Matrix_Classes::eps_Real*mu_xz-z)/dz));
    if (alpha>CH_Matrix_Classes::eps_Real){
      CH_Matrix_Classes::Real otheta=std::fabs(xz)/mu_xz;
      CH_Matrix_Classes::Real gtheta=nbh_ubnd;
      CH_Matrix_Classes::Real theta=CH_Matrix_Classes::max(otheta,nbh_ubnd);
      if (otheta>(1+1e-6)*nbh_ubnd){
        gtheta=std::sqrt(CH_Matrix_Classes::max(0.,otheta*otheta+2*(xz*xdzpdxz-otheta*otheta*mu_xz*mu_xdzpdxz)/CH_Matrix_Classes::max(mu_at_one,1e-6*mu_xz)));
        if (gtheta<.1*nbh_ubnd+0.9*otheta){
          gtheta=CH_Matrix_Classes::max(.1*gtheta+.9*otheta,nbh_ubnd);
        }
        else {
          gtheta=(1-1e-3)*otheta+1e-3*nbh_ubnd;
        }
      }

      //find the longest step for the positive branch of the absolute value
      CH_Matrix_Classes::Real q0 = CH_Matrix_Classes::min(xz - theta*mu_xz,0.);
      CH_Matrix_Classes::Real q1 = xdzpdxz - theta*mu_xdzpdxz-(gtheta-theta)*(mu_at_one);
      CH_Matrix_Classes::Real q2 = dxdz - theta*mu_dxdz;

      if (q0+1e-8*q1>CH_Matrix_Classes::eps_Real*mu_xz) {
        alpha=0.;
      }
      else if (q2>CH_Matrix_Classes::eps_Real*mu_xz){
        //strictly convex case with min >= zero, max is upper root
        CH_Matrix_Classes::Real p=q1/q2/2.;
        CH_Matrix_Classes::Real q=q0/q2;
        alpha=CH_Matrix_Classes::min(alpha,-p+std::sqrt(p*p-q));
      }
      else if (q2<-CH_Matrix_Classes::eps_Real*mu_xz){
        //strictly concave case, restrictions only for positive q1
        if (q1>0.) {
          CH_Matrix_Classes::Real p=q1/q2/2.;
          CH_Matrix_Classes::Real q=q0/q2;
          CH_Matrix_Classes::Real discr=p*p-q;
          if (discr>0){
            alpha=CH_Matrix_Classes::min(alpha,-p-std::sqrt(discr));
          }
        }
      }
      else {
        //q is affine
        if (q0+alpha*q1>0.) {
          alpha=CH_Matrix_Classes::min(alpha,CH_Matrix_Classes::max(0.,-q0/q1));
        }
      }

      //find the longest step for the negative branch of the absolute value
      q0 = -xz - theta*mu_xz;
      q1 = -xdzpdxz - theta*mu_xdzpdxz-(gtheta-theta)*(mu_at_one);
      q2 = -dxdz - theta*mu_dxdz;
      if (q0+1e-8*q1>CH_Matrix_Classes::eps_Real*mu_xz) {
        alpha=0;
      }
      else if (q2>CH_Matrix_Classes::eps_Real*mu_xz){
        //strictly convex case with min >= zero, max is upper root
        CH_Matrix_Classes::Real p=q1/q2/2.;
        CH_Matrix_Classes::Real q=q0/q2;
        alpha=CH_Matrix_Classes::min(alpha,-p+std::sqrt(p*p-q));
      }
      else if (q2<-CH_Matrix_Classes::eps_Real*mu_xz){
        //strictly concave case, restrictions only for positive q1
        if (q1>0.) {
          CH_Matrix_Classes::Real p=q1/q2/2.;
          CH_Matrix_Classes::Real q=q0/q2;
          CH_Matrix_Classes::Real discr=p*p-q;
          if (discr>0){
            alpha=CH_Matrix_Classes::min(alpha,-p-std::sqrt(discr));
          }
        }
      }
      else {
          //q is affine
        if (q0+alpha*q1>0.) {
          alpha=CH_Matrix_Classes::min(alpha,CH_Matrix_Classes::max(0.,-q0/q1));
        }
      }
    }

    CH_Matrix_Classes::Real step_nbh=std::fabs(xz+alpha*(xdzpdxz+alpha*dxdz))/CH_Matrix_Classes::max(mu_xz+alpha*(mu_xdzpdxz+alpha*(mu_dxdz)),1e-6*mu_xz);
    max_nbh=CH_Matrix_Classes::max(max_nbh,step_nbh);
  }



}

#endif