wangjie212 / TSSOS

A sparse polynomial optimization tool based on the moment-SOS hierarchy.
MIT License
39 stars 14 forks source link

BoundsError when solution=true #5

Closed dibondar closed 2 years ago

dibondar commented 2 years ago

Dear Jie Wang, this is a great library! I am currently using in a variety of physics projects.

If run on Julia 1.7.2 and Mosek 10.0

using DynamicPolynomials
using TSSOS

@polyvar x
@polyvar y

obj = (x * y)^2

opt,sol,data = tssos_first(obj, variables(obj), solution = true)

I get the BoundsError when extracting solutions, as show below

************************TSSOS************************
TSSOS is launching...
Starting to compute the block structure...
------------------------------------------------------
The sizes of PSD blocks:
[1]
[2]
------------------------------------------------------
Obtained the block structure. The maximal size of blocks is 1.
Assembling the SDP...
There are 7 affine constraints.
Solving the SDP...
Problem
  Name                   :                 
  Objective sense        : max             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 7               
  Cones                  : 0               
  Scalar variables       : 3               
  Matrix variables       : 1               
  Integer variables      : 0               

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 2                 time                   : 0.00            
Lin. dep.  - tries                  : 1                 time                   : 0.00            
Lin. dep.  - number                 : 0               
Presolve terminated. Time: 0.00    
Problem
  Name                   :                 
  Objective sense        : max             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 7               
  Cones                  : 0               
  Scalar variables       : 3               
  Matrix variables       : 1               
  Integer variables      : 0               

Optimizer  - threads                : 56              
Optimizer  - solved problem         : the primal      
Optimizer  - Constraints            : 3
Optimizer  - Cones                  : 1
Optimizer  - Scalar variables       : 3                 conic                  : 2               
Optimizer  - Semi-definite variables: 1                 scalarized             : 6               
Factor     - setup time             : 0.00              dense det. time        : 0.00            
Factor     - ML order time          : 0.00              GP order time          : 0.00            
Factor     - nonzeros before factor : 6                 after factor           : 6               
Factor     - dense dim.             : 0                 flops                  : 1.40e+02        
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  
0   2.0e+00  1.0e+00  1.0e+00  0.00e+00   0.000000000e+00   0.000000000e+00   1.0e+00  0.00  
1   4.9e-01  2.5e-01  1.7e-01  1.09e+00   4.469925023e-01   2.740939884e-01   2.5e-01  0.00  
2   6.2e-02  3.1e-02  5.6e-03  1.16e+00   2.339419280e-02   2.159644755e-02   3.1e-02  0.00  
3   2.7e-05  1.4e-05  3.5e-08  1.03e+00   -1.023482533e-06  6.032912939e-06   1.4e-05  0.00  
4   1.4e-10  6.8e-11  1.8e-13  1.00e+00   1.358038542e-10   2.040584821e-10   6.8e-11  0.00  
Optimizer terminated. Time: 0.01    

SDP solving time: 2.040073584 seconds.
optimum = 1.3580385424324652e-10

BoundsError: attempt to access 2×2 Matrix{UInt8} at index [1:2, 3]

Stacktrace:
 [1] throw_boundserror(A::Matrix{UInt8}, I::Tuple{Base.Slice{Base.OneTo{Int64}}, Int64})
   @ Base ./abstractarray.jl:691
 [2] checkbounds
   @ ./abstractarray.jl:656 [inlined]
 [3] _getindex
   @ ./multidimensional.jl:838 [inlined]
 [4] getindex
   @ ./abstractarray.jl:1218 [inlined]
 [5] blockupop(n::Int64, supp::Matrix{UInt8}, coe::Vector{Int64}, basis::Matrix{UInt8}, blocks::Vector{Vector{Int64}}, cl::Int64, blocksize::Vector{Int64}; nb::Int64, solver::String, feasible::Bool, QUIET::Bool, solve::Bool, solution::Bool, MomentOne::Bool, Gram::Bool)
   @ TSSOS ~/.julia/packages/TSSOS/drkzz/src/blockpop_uncons.jl:478
 [6] tssos_first(f::Monomial{true}, x::Vector{PolyVar{true}}; nb::Int64, newton::Bool, reducebasis::Bool, TS::String, merge::Bool, feasible::Bool, md::Int64, solver::String, QUIET::Bool, solve::Bool, MomentOne::Bool, Gram::Bool, solution::Bool, tol::Float64)
   @ TSSOS ~/.julia/packages/TSSOS/drkzz/src/blockpop_uncons.jl:81
 [7] top-level scope
   @ In[3]:9
 [8] eval
   @ ./boot.jl:373 [inlined]
 [9] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1196
wangjie212 commented 2 years ago

Dear Dibondar, this issue is caused because the Newton polytope method removes some monomials that are necessary for solution extraction. You can try

opt,sol,data = tssos_first(obj, variables(obj), newton=false, solution = true)

It should work.

dibondar commented 2 years ago

Indeed, it does work! Thank you!