BoxIPBundleBlock.hxx

Go to the documentation of this file.

#ifndef CONICBUNDLE_BOXIPBUNDLEBLOCK_HXX
#define CONICBUNDLE_BOXIPBUNDLEBLOCK_HXX

#include "InteriorPointBundleBlock.hxx"

namespace ConicBundle {


class BoxIPBundleBlock: public virtual InteriorPointBundleBlock
{
protected:
  CH_Matrix_Classes::Integer vecdim;  
  CH_Matrix_Classes::Integer scalvecdim; 
  CH_Matrix_Classes::Matrix lb; 
  CH_Matrix_Classes::Matrix ub; 
  bool use_scaling;                
  CH_Matrix_Classes::Real scalub; 
  
  CH_Matrix_Classes::Matrix x;     
  CH_Matrix_Classes::Matrix lz;     
  CH_Matrix_Classes::Matrix uz;     
  CH_Matrix_Classes::Matrix dx;    
  CH_Matrix_Classes::Matrix dlz;    
  CH_Matrix_Classes::Matrix duz;    

  CH_Matrix_Classes::Real s;       
  CH_Matrix_Classes::Real lt;       
  CH_Matrix_Classes::Real ut;       
  CH_Matrix_Classes::Real ds;      
  CH_Matrix_Classes::Real dlt;      
  CH_Matrix_Classes::Real dut;      

  CH_Matrix_Classes::Matrix xiz;        
  CH_Matrix_Classes::Matrix lxinv;      
  CH_Matrix_Classes::Matrix uxinv;      
  CH_Matrix_Classes::Real hats;       
  CH_Matrix_Classes::Matrix hatb;       
  CH_Matrix_Classes::Matrix sqrt_xiz; 
  CH_Matrix_Classes::Real Chols;      
  CH_Matrix_Classes::Matrix Cholinvb; 
  
  CH_Matrix_Classes::Matrix compl_lrhs;  
  CH_Matrix_Classes::Matrix compl_urhs;  
  CH_Matrix_Classes::Real compl_ltrhs; 
  CH_Matrix_Classes::Real compl_utrhs; 
  CH_Matrix_Classes::Real sys_srhs; 
  
  CH_Matrix_Classes::Real last_rhs_mu;  
  CH_Matrix_Classes::Real mu;           
  CH_Matrix_Classes::Real old_mu;       

  CH_Matrix_Classes::Real last_alpha;   
  
  CH_Matrix_Classes::Matrix oldx;       
  CH_Matrix_Classes::Real olds;          
  CH_Matrix_Classes::Matrix oldlz;       
  CH_Matrix_Classes::Matrix olduz;       
  CH_Matrix_Classes::Real oldlt;         
  CH_Matrix_Classes::Real oldut;         
  
  CH_Matrix_Classes::Matrix Bt;       
  CH_Matrix_Classes::Matrix Boffset; 
  CH_Matrix_Classes::Matrix sqrBnorms; 

  mutable CH_Matrix_Classes::Matrix tmpvec;  
  mutable CH_Matrix_Classes::Matrix tmpmat;  
 
  CH_Matrix_Classes::Indexmatrix map_to_old; 

  bool test_myself_call; 

  int find_inactive_indices(CH_Matrix_Classes::Indexmatrix& inactive,
                            CH_Matrix_Classes::Real trace_bound,
                            bool cautious) const;
  
  void point_changed();

  int compute_NTscaling(void);

  int given_dx_compute_dz(void);
  
public:
  void clear(const CH_Matrix_Classes::Matrix& lb,
             const CH_Matrix_Classes::Matrix& ub,
             bool use_scaling,
             CH_Matrix_Classes::Real scalub);

  BoxIPBundleBlock(const CH_Matrix_Classes::Matrix& lb,
                 const CH_Matrix_Classes::Matrix& ub,
                 bool use_scaling=false,
                 CH_Matrix_Classes::Real scalub=-1.,
                 CBout* cb=0,int cbinc=-1);

  ~BoxIPBundleBlock();

  virtual InteriorPointBundleBlock* clone();

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

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

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

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

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

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

  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;

  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;

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

  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);
  
  
  virtual int set_dx(const CH_Matrix_Classes::Matrix& rhs,CH_Matrix_Classes::Integer startindex);

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

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

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

  virtual int do_step(CH_Matrix_Classes::Real alpha);

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

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



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

  virtual int get_vecx(CH_Matrix_Classes::Matrix& vecx,CH_Matrix_Classes::Integer startindex)
  {
    CH_Matrix_Classes::mat_xey(vecdim,vecx.get_store()+startindex,x.get_store());
    if (use_scaling)
      vecx(startindex+vecdim)=s;
    return 0; }
  
  virtual int get_vecz(CH_Matrix_Classes::Matrix& vecz,CH_Matrix_Classes::Integer startindex)
  {
    CH_Matrix_Classes::mat_xey(vecdim,vecz.get_store()+startindex,lz.get_store());
    CH_Matrix_Classes::mat_xmey(vecdim,vecz.get_store()+startindex,uz.get_store());
    if (use_scaling)
      vecz(startindex+vecdim)=lt-(scalub>0.?ut:0.);
    return 0;
  }
  
  virtual int get_vecdx(CH_Matrix_Classes::Matrix& vecdx,CH_Matrix_Classes::Integer startindex)
  {
    if (dx.dim()!=vecdim)
      return 1;
    CH_Matrix_Classes::mat_xey(vecdim,vecdx.get_store()+startindex,dx.get_store());
    if (use_scaling)
      vecdx(startindex+vecdim)=ds;
    return 0;
  }
  
  
  virtual int get_vecdz(CH_Matrix_Classes::Matrix& vecdz,CH_Matrix_Classes::Integer startindex)
  {
    if (dlz.dim()!=vecdim)
      return 1;
    CH_Matrix_Classes::mat_xey(vecdim,vecdz.get_store()+startindex,dlz.get_store());
    CH_Matrix_Classes::mat_xmey(vecdim,vecdz.get_store()+startindex,duz.get_store());
    if (use_scaling)
      vecdz(startindex+vecdim)=dlt-(scalub>0.?dut:0.);
    return 0;
  }

    virtual CH_Matrix_Classes::Integer dim_bundle() const
  {return bundle_dim;}
  
  virtual CH_Matrix_Classes::Real evaluate_trace_x();

  virtual CH_Matrix_Classes::Real evaluate_trace_z();

  virtual CH_Matrix_Classes::Real evaluate_trace_dx();

  virtual int add_trace(CH_Matrix_Classes::Matrix& vec,
                        CH_Matrix_Classes::Real alpha,
                        CH_Matrix_Classes::Integer startindex);
  
  virtual int set_xizinvsqrt_trace(CH_Matrix_Classes::Matrix& vec,
                                   CH_Matrix_Classes::Integer startindex);
  
  // after the bundle subproblem is solved, this retrieves the local linear solution vector; if linx_activity is set, the values between zero and one indicate the guess on the coefficients activity level 
  // virtual int get_linx(CH_Matrix_Classes::Matrix& linx,
  //                   CH_Matrix_Classes::Matrix* linx_activity=0) const;

  virtual CH_Matrix_Classes::Matrix&
  B_times(const CH_Matrix_Classes::Matrix& A,
          CH_Matrix_Classes::Matrix& C,
          CH_Matrix_Classes::Real alpha,
          CH_Matrix_Classes::Real beta,
          int Btrans,
          int Atrans,
          CH_Matrix_Classes::Integer startindex_model,
          MinorantBundle& globalbundle,
          CH_Matrix_Classes::Integer startindex_bundle);
                                                   
  virtual CH_Matrix_Classes::Matrix&
  times_B(const CH_Matrix_Classes::Matrix& A,
          CH_Matrix_Classes::Matrix& C,
          CH_Matrix_Classes::Real alpha,
          CH_Matrix_Classes::Real beta,
          int Atrans,
          int Btrans,
          CH_Matrix_Classes::Integer startindex_model,
          MinorantBundle& globalbundle,
          CH_Matrix_Classes::Integer startindex_bundle);

  virtual CH_Matrix_Classes::Symmatrix&
  add_BDBt(const CH_Matrix_Classes::Matrix& diagvec,
           CH_Matrix_Classes::Symmatrix& bigS,
           bool minus,
           CH_Matrix_Classes::Integer startindex,
           CH_Matrix_Classes::Matrix& Bt,
           CH_Matrix_Classes::Integer startindex_model,
           MinorantBundle& globalbundle,
           CH_Matrix_Classes::Integer startindex_bundle);

  virtual CH_Matrix_Classes::Matrix&
  get_Bt(CH_Matrix_Classes::Matrix& Bt,
         CH_Matrix_Classes::Integer startindex_model,
         MinorantBundle& global_bundle,
         CH_Matrix_Classes::Integer startindex_bundle);
  
  virtual int add_modelx_aggregate(CH_Matrix_Classes::Real& val,
                                   CH_Matrix_Classes::Matrix& vec,
                                   MinorantBundle& global_bundle,
                                   CH_Matrix_Classes::Integer startindex_bundle);

  virtual int get_sysviol_model(CH_Matrix_Classes::Matrix& sysviol_model,
                                CH_Matrix_Classes::Integer startindex_model,
                                const CH_Matrix_Classes::Matrix& dy,
                                const CH_Matrix_Classes::Real deltatrdual,
                                MinorantBundle& global_bundle,
                                CH_Matrix_Classes::Integer startindex_bundle);

  virtual int set_bundle_z(const CH_Matrix_Classes::Matrix& y,
                    MinorantBundle& global_bundle,
                    CH_Matrix_Classes::Integer startindex_bundle,
                    CH_Matrix_Classes::Real& add_center_value);

  virtual int add_trace_to_diff_model(CH_Matrix_Classes::Real trace_dual);

  virtual CH_Matrix_Classes::Real dualviol_2normsqr();

  virtual int do_bundle_step(CH_Matrix_Classes::Real alpha,
                      const CH_Matrix_Classes::Matrix& y,
                      MinorantBundle& globalbundle,
                      CH_Matrix_Classes::Integer startindex_bundle,
                      CH_Matrix_Classes::Real tracedual,
                      CH_Matrix_Classes::Real trace_rhs);

  virtual int set_modelrhs(CH_Matrix_Classes::Matrix& globalrhs, 
                           CH_Matrix_Classes::Real rhsmu,
                           CH_Matrix_Classes::Real rhscorr,
                           CH_Matrix_Classes::Integer startindex_model);

  virtual int add_BtinvsysB(CH_Matrix_Classes::Symmatrix& globalsys,
                            const MinorantBundle& bundle,
                            CH_Matrix_Classes::Integer startindex_bundle) ;

  int Schur_transform_bundle(CH_Matrix_Classes::Matrix& glob_lowrank,
                             MinorantBundle& globalbundle,
                             CH_Matrix_Classes::Integer startindex_bundle,
                             CH_Matrix_Classes::Matrix& trafotrace,
                             CH_Matrix_Classes::Integer startindex_trace);
  
  virtual int add_bundle_xizinv_diagonal(CH_Matrix_Classes::Matrix& diagonal,
                                         CH_Matrix_Classes::Matrix& ipBtrvec,
                                         MinorantBundle& globalbundle,
                                         CH_Matrix_Classes::Integer startindex_bundle,
                                         const CH_Matrix_Classes::Matrix& trafotrace,
                                         CH_Matrix_Classes::Integer startindex_trace);
                                         

  virtual int add_pcsubspace(CH_Matrix_Classes::Matrix& lowrank,
                             CH_Matrix_Classes::Matrix& sigma_guess,
                             const CH_Matrix_Classes::Matrix& Diag_inv,
                             CH_Matrix_Classes::Real minval,
                             CH_Matrix_Classes::Real diaginvval,
                             CH_Matrix_Classes::Matrix & minus_trmult,
                             CH_Matrix_Classes::Real schur_trace,
                             MinorantBundle& globalbundle,
                             CH_Matrix_Classes::Integer startindex_bundle);


  int add_bundle_xizinvsqrt_projection(CH_Matrix_Classes::Matrix& glob_lowrank,
                                       CH_Matrix_Classes::Matrix& subspace,
                                       CH_Matrix_Classes::Integer startindex_subsspace,
                                       MinorantBundle& globalbundle,
                                       CH_Matrix_Classes::Integer startindex_bundle);

  virtual int add_BtinvsysB_times(const CH_Matrix_Classes::Matrix& in_vec,
                                  CH_Matrix_Classes::Matrix& out_vec,
                                  CH_Matrix_Classes::Real zeta_inval,
                                  CH_Matrix_Classes::Real* zeta_outval,
                                  MinorantBundle& globalbundle,
                                  CH_Matrix_Classes::Integer startindex_bundle);

  virtual int set_dx_xizsolvestep(const CH_Matrix_Classes::Matrix& step_y,
                                 const CH_Matrix_Classes::Real step_trdual,
                                 MinorantBundle& globalbundle,
                                 CH_Matrix_Classes::Integer startindex_bundle);


  int get_boxx(CH_Matrix_Classes::Matrix& linx,
               CH_Matrix_Classes::Matrix* linx_activity,
               CH_Matrix_Classes::Real trace_rhs,
               bool cautious) const;
  
  CH_Matrix_Classes::Real constraints_cost()
  {
    if (use_scaling){
      if (scalub>0.) {
        return scalub*ut;
      }
      else {
        return 0.;
      }
    }
    CH_Matrix_Classes::Real val= CH_Matrix_Classes::ip(ub,uz)-CH_Matrix_Classes::ip(lb,lz);
    return val;
  }
  
};

  

}

#endif