UQPConeModelBlock.hxx

Go to the documentation of this file.

#ifndef CONICBUNDLE_UQPCONEMODELBLOCK_HXX
#define CONICBUNDLE_UQPCONEMODELBLOCK_HXX

#include <vector>
#include "UQPModelBlock.hxx"


namespace ConicBundle {


  class UQPConeModelBlock: public UQPModelBlock, public QPConeModelDataObject
{
private:
  CH_Matrix_Classes::Matrix x;

  CH_Matrix_Classes::Matrix z;                 
    
  CH_Matrix_Classes::Matrix A;
  
  CH_Matrix_Classes::Real b;
  
  int less_or_equal;           
  CH_Matrix_Classes::Real s;      
  CH_Matrix_Classes::Real scost;  
  CH_Matrix_Classes::Real y;      
 
  CH_Matrix_Classes::Integer nnc_dim;    

  CH_Matrix_Classes::Indexmatrix soc_dim; 
  CH_Matrix_Classes::Indexmatrix soc_start; 

  std::vector<CH_Matrix_Classes::Symmatrix> Xp; 
  std::vector<CH_Matrix_Classes::Symmatrix> Zp; 
  std::vector<CH_Matrix_Classes::Symmatrix> Xinv; 
  //std::vector<CH_Matrix_Classes::Symmatrix> Xchol; //inverses of local primal psd variables X 
  //std::vector<CH_Matrix_Classes::Symmatrix> Zchol; //inverses of local primal psd variables Z

  CH_Matrix_Classes::Matrix lb;    
  CH_Matrix_Classes::Matrix ub;    
  bool box_scaling; 
  bool box_scaleub; 


  //CH_Matrix_Classes::Matrix lz;    // box dual variable for lower bounds is in z
  CH_Matrix_Classes::Matrix uz;    
  //CH_Matrix_Classes::Real box_lz;  //dual variable for lower bound on box_s is in z
  CH_Matrix_Classes::Real box_uz;  

  
  CH_Matrix_Classes::Integer box_dim;  

  CH_Matrix_Classes::Integer box_start;  


  CH_Matrix_Classes::Integer qp_xstart; 
  CH_Matrix_Classes::Integer qp_ystart; 

  CH_Matrix_Classes::Integer mu_dim;    
  CH_Matrix_Classes::Real restart_factor; 


  //---- copies of old variables for analyzing activity
  CH_Matrix_Classes::Real current_mu;   
  CH_Matrix_Classes::Real old_mu;  
  CH_Matrix_Classes::Matrix old_x;  
  CH_Matrix_Classes::Matrix old_z;  
  CH_Matrix_Classes::Matrix old_uz;   
  CH_Matrix_Classes::Real old_box_uz;  
  CH_Matrix_Classes::Real old_y;  
  CH_Matrix_Classes::Real old_s;  
  std::vector<CH_Matrix_Classes::Symmatrix> old_Xp; 
  std::vector<CH_Matrix_Classes::Symmatrix> old_Zp; 
  CH_Matrix_Classes::Real last_alpha; 

  //---- for forming or testing the corrector the old step is copied here
  CH_Matrix_Classes::Matrix zcorr;   
  CH_Matrix_Classes::Matrix uzcorr;  
  CH_Matrix_Classes::Real box_uzcorr; 
  CH_Matrix_Classes::Real scorr;    
  
  mutable CH_Matrix_Classes::Matrix dx;    
  mutable CH_Matrix_Classes::Real dy;      
  mutable CH_Matrix_Classes::Matrix dz;    
  mutable CH_Matrix_Classes::Matrix duz;   
  mutable CH_Matrix_Classes::Real box_duz; 
  mutable CH_Matrix_Classes::Real ds;      
  
  //---- temporary variables
  mutable CH_Matrix_Classes::Matrix tmpvec;  
  mutable CH_Matrix_Classes::Matrix tmpsvec; 
  mutable CH_Matrix_Classes::Symmatrix tmpsymmat; 
  mutable CH_Matrix_Classes::Symmatrix dX; 
 
  //the Block has to memorize and update its local x,y,z variables 
  //(and all additional ones that are not visible to the outside)
  //in particular it is assumed to keep a copy of the current variables
  //after the calls to the following functions:
  //  starting_x, starting_yz, restart_x, restart_yz, set_point

  CH_Matrix_Classes::Matrix mult_Arw(const CH_Matrix_Classes::Matrix& x,const CH_Matrix_Classes::Matrix& v) const;

  CH_Matrix_Classes::Matrix mult_Arwinv(const CH_Matrix_Classes::Matrix& x,const CH_Matrix_Classes::Matrix& v) const;

  CH_Matrix_Classes::Matrix mult_G(const CH_Matrix_Classes::Matrix& x,CH_Matrix_Classes::Real gamma,const CH_Matrix_Classes::Matrix& v) const;

  CH_Matrix_Classes::Matrix mult_Ginv(const CH_Matrix_Classes::Matrix& x,CH_Matrix_Classes::Real gamma,const CH_Matrix_Classes::Matrix& v) const;



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

  ~UQPConeModelBlock(){}  

  void clear();

  
  int init(const MinorantPointer& constant_minorant,
           const MinorantBundle& bundle,
           CH_Matrix_Classes::Integer nnc_dim,
           const CH_Matrix_Classes::Indexmatrix& soc_dim,
           const CH_Matrix_Classes::Indexmatrix& sdp_dim,
           const CH_Matrix_Classes::Matrix& box_lb,
           const CH_Matrix_Classes::Matrix& box_ub,
           CH_Matrix_Classes::Real b,
           FunctionTask ft,
           QPModelOracleDataObject* oracle_data=0,
           bool scale_box=true); 

  int adjust_trace(CH_Matrix_Classes::Real b); 

 CH_Matrix_Classes::Real evaluate_trace() const ; 

  CH_Matrix_Classes::Real get_trace(){return b;}
  
 int get_nncx(CH_Matrix_Classes::Matrix& nncx,
               CH_Matrix_Classes::Matrix* nncx_activity=0,
               bool cautious=false);

 int get_socx(CH_Matrix_Classes::Integer i,
               CH_Matrix_Classes::Matrix& socx,
               CH_Matrix_Classes::Real* socx_activity,
               bool cautious=false);
  
  int get_pscx(CH_Matrix_Classes::Integer i,
               CH_Matrix_Classes::Matrix& pscx_eigs,
               CH_Matrix_Classes::Matrix& pscx_vecs,
               CH_Matrix_Classes::Real& pscx_growthrate,
               CH_Matrix_Classes::Matrix& pscx_primalgrowth,
               CH_Matrix_Classes::Matrix& pscx_dualgrowth);
  
  int get_boxx(CH_Matrix_Classes::Matrix& /* linx */,
               CH_Matrix_Classes::Matrix* /* linx_activity */ =0 ,
               bool cautious=false);

  int add_modelx_aggregate(CH_Matrix_Classes::Real& offset,
                           CH_Matrix_Classes::Matrix& gradient);


  CH_Matrix_Classes::Real tracedual(CH_Matrix_Classes::Real* prec=0) const;

  int get_nncz(CH_Matrix_Classes::Matrix& vecz){ vecz.init(nnc_dim,1,z.get_store());return 0;}
  int get_socz(CH_Matrix_Classes::Integer i,CH_Matrix_Classes::Matrix& vecz)
  { vecz.init(soc_dim(i),1,z.get_store()+soc_start(i));return 0;}
  int get_X(CH_Matrix_Classes::Integer i,CH_Matrix_Classes::Symmatrix& X){X.init(Xp[(unsigned long)(i)]); return 0;}
  int get_Z(CH_Matrix_Classes::Integer i,CH_Matrix_Classes::Symmatrix& Z){Z.init(Zp[(unsigned long)(i)]); return 0;}
  CH_Matrix_Classes::Real get_y(void){return y;}
  CH_Matrix_Classes::Real get_s(void){if (less_or_equal) return s; return 0.;}
  CH_Matrix_Classes::Real get_mu(void){return current_mu;}
  

  int get_old_nncx(CH_Matrix_Classes::Matrix& vecx){ vecx.init(nnc_dim,1,old_x.get_store());return 0;}
  int get_old_nncz(CH_Matrix_Classes::Matrix& vecz){ vecz.init(nnc_dim,1,old_z.get_store());return 0;}
  int get_old_socx(CH_Matrix_Classes::Integer i,CH_Matrix_Classes::Matrix& vecx)
  { vecx.init(soc_dim(i),1,old_x.get_store()+soc_start(i));return 0;}
  int get_old_socz(CH_Matrix_Classes::Integer i,CH_Matrix_Classes::Matrix& vecz)
  { vecz.init(soc_dim(i),1,old_z.get_store()+soc_start(i));return 0;}
  int get_old_X(CH_Matrix_Classes::Integer i,CH_Matrix_Classes::Symmatrix& X){X.init(old_Xp[(unsigned long)(i)]); return 0;}
  int get_old_Z(CH_Matrix_Classes::Integer i,CH_Matrix_Classes::Symmatrix& Z){Z.init(old_Zp[(unsigned long)(i)]); return 0;}
  CH_Matrix_Classes::Real get_old_y(void){return old_y;}
  CH_Matrix_Classes::Real get_old_s(void){if (less_or_equal) return old_s; return 0.;}
  CH_Matrix_Classes::Real get_old_mu(void){return old_mu;}
  CH_Matrix_Classes::Real get_last_alpha(){return last_alpha;}

  int compute_local_directions(const CH_Matrix_Classes::Matrix& qp_dx,
                               const CH_Matrix_Classes::Matrix& qp_dy,
                               const CH_Matrix_Classes::Matrix& rhs_resid,
                               CH_Matrix_Classes::Matrix& dz,
                               CH_Matrix_Classes::Matrix& duz,
                               CH_Matrix_Classes::Real& box_ds,
                               CH_Matrix_Classes::Real& ds);
  
  int inner_line_search(CH_Matrix_Classes::Real& alpha,
                        const CH_Matrix_Classes::Matrix& qp_dx,
                        const CH_Matrix_Classes::Matrix& qp_dy,
                        const CH_Matrix_Classes::Matrix& dz,
                        const CH_Matrix_Classes::Matrix& duz,
                        CH_Matrix_Classes::Real box_ds,
                        CH_Matrix_Classes::Real ds);

  //-----------  QP_Block routines (see there)

  CH_Matrix_Classes::Integer xdim() const {return x.rowdim();}
  
  CH_Matrix_Classes::Integer ydim() const {return A.rowdim();}
  
  int set_qp_xstart(CH_Matrix_Classes::Integer x_start_index)
    {qp_xstart=x_start_index; return 0;}
  
  int set_qp_ystart(CH_Matrix_Classes::Integer y_start_index)
    {qp_ystart=y_start_index; return 0;}

  int starting_x(CH_Matrix_Classes::Matrix& qp_x);

  int starting_y(CH_Matrix_Classes::Matrix& qp_y,
                 const CH_Matrix_Classes::Matrix& qp_Qx,
                 const CH_Matrix_Classes::Matrix& qp_c);

  int get_Ab(CH_Matrix_Classes::Matrix& qp_A,CH_Matrix_Classes::Matrix &qp_b) const;

 int restart_x(CH_Matrix_Classes::Matrix& qp_x,
                const CH_Matrix_Classes::Matrix& qp_c,
                const CH_Matrix_Classes::Matrix& qp_dc);

  int restart_y(CH_Matrix_Classes::Matrix& qp_y,
                const CH_Matrix_Classes::Matrix& qp_Qx,
                const CH_Matrix_Classes::Matrix& qp_c,
                const CH_Matrix_Classes::Matrix& qp_dc);
    
 int add_xinv_kron_z(CH_Matrix_Classes::Symmatrix& barQ);

  int add_local_sys(CH_Matrix_Classes::Symmatrix& sysdy,
                    CH_Matrix_Classes::Matrix& rhs);

 CH_Matrix_Classes::Real get_local_primalcost() const 
  {return (less_or_equal? s*scost:0.);}

  CH_Matrix_Classes::Real get_local_dualcost() const; 

  int suggest_mu(CH_Matrix_Classes::Real& ip_xz,
                 CH_Matrix_Classes::Integer& mu_dim,
                 CH_Matrix_Classes::Real& sigma,
                 const CH_Matrix_Classes::Matrix& qp_dx,
                 const CH_Matrix_Classes::Matrix& qp_dy,
                 const CH_Matrix_Classes::Matrix& rhs_residual);

  int get_corr(CH_Matrix_Classes::Matrix& xcorr,
               CH_Matrix_Classes::Matrix& rhs,
               CH_Matrix_Classes::Real mu);

 int line_search(CH_Matrix_Classes::Real& alpha,
                  const CH_Matrix_Classes::Matrix& qp_dx,
                  const CH_Matrix_Classes::Matrix& qp_dy,
                  const CH_Matrix_Classes::Matrix& rhs_residual);

  int set_point(const CH_Matrix_Classes::Matrix& qp_x,
                const CH_Matrix_Classes::Matrix& qp_y,
                CH_Matrix_Classes::Real alpha);




  //---------------- for debugging purposes

  CH_Matrix_Classes::Matrix& add_Bs(CH_Matrix_Classes::Matrix& qp_vec) const;

 CH_Matrix_Classes::Matrix& subtract_z(CH_Matrix_Classes::Matrix& dual_residual,bool with_step=false) const;

   
};



}

#endif