2025年11月30日日曜日

Cuopt評価

 COPTに対して、Cuoptを評価しました。インスタンスは、instance24ですが、LB値を求める途中のものを使用しています。

下記に結果を示しますが、評価環境が違うことに注意してください。COPTの評価環境は、CPU/GPU共Cuoptの評価環境のおよそ倍以上の性能が違います。Cuopt、BarrierGPUでは、インスタンス23でも数値エラーとなっていました。インスタンス22では、計測できましたが、PDLPに比べかなり遅い結果となっており、仮に数値演算エラーが無くなったとしても期待できない、と考えます。ちなみにBarrierSolverは、Cutoptに限らず、数値演算エラーが生じ易いようで、数値演算エラーが出たときの救済策は考えておく必要があります。




<COPT PDLPとCOPT PDLPの差>

恐らくは、COPTが未だPDLPXを実装していない、のが原因と思われます。ライセンス評価期間内に実装されることを望みます。(サポートと話をする機会があれば、聞いてみたいと思います。)

GPUの性能が約3.7倍であることを考慮すると、Cuopt PDLPは、COPT Barrierよりも高速という結論になります。

Cuoptは、DOCKERで計算サーバを立ち上げることもサポートしているので、Windowsの実装のままWSL2上にローカルサーバを構築し、REST APIでLP計算を依頼する実装も可能です。また、それとは別に、PDLPXは、C++実装が公開され、WarmStartもサポートしています。(Branch操作において、WarmStartは必須です。)これをC++の開発中ソルバに組み込むことも可能です。COPTのREST APIを実装する手間と、C++組み込みの手間とどちらがよいか悩みところですが、PDLPX化の課題は、それだけではなく、実務的には、crossoverを実装する必要があります。幸い同じ著者による実装が公開されていますが、この辺のところまでコントロールするなると結構大変です。全てが揃っているCOPTがやってくれれば、余計なことに悩む必要はありません。

https://schedule-nurse.blogspot.com/2025/12/copt-pdlppdlpx.html

以下は、評価時のログです。Cuoptでは、LPファイルをサポートしていないようでしたので、COPTで再びMPSにフォーマットし直して評価しています。



COPT> set lpmethod 2
Setting parameter 'LpMethod' to 2
COPT> set crossover 0
Setting parameter 'Crossover' to 0
COPT> read instance24.lp
Reading from 'C:\Users\.PC\highs_test\instance24.lp'
Reading finished (0.62s)
COPT> opt
Model fingerprint: 11cebb4f

Using Cardinal Optimizer v8.0.1 on Windows
Hardware has 8 cores and 16 threads. Using instruction set X86_NATIVE (1)
Minimizing an LP problem

The original problem has:
    14135 rows, 84693 columns and 11800327 non-zero elements
The presolved problem has:
    13108 rows, 81336 columns and 11541827 non-zero elements

Starting barrier solver using 8 CPU threads

Problem info:
    Dualized in presolve:            No
    Range of matrix coefficients:    [6e-02,5e+01]
    Range of rhs coefficients:       [4e-01,7e+01]
    Range of bound coefficients:     [1e+00,1e+00]
    Range of cost coefficients:      [1e+00,4e+04]

Factor info:
    Number of free columns:          0
    Number of dense columns:         0
    Number of matrix entries:        6.485e+07
    Number of factor entries:        7.216e+07
    Number of factor flops:          5.209e+11

Iter       Primal.Obj         Dual.Obj      Compl  Primal.Inf  Dual.Inf    Time
   0  +6.89041521e+07  -1.47285422e+08   5.11e+08    2.93e+03  3.86e+04  20.29s
   1  +2.06731774e+07  -4.05798355e+07   7.96e+07    2.95e+02  4.35e+03  23.00s
   2  +1.51129327e+06  -2.46845951e+06   4.38e+06    1.09e+01  1.98e+02  25.97s
   3  +2.34855347e+05  -3.81229124e+05   6.46e+05    1.21e+00  2.72e+01  28.98s
   4  +9.01673726e+04  -6.33129847e+04   1.57e+05    2.42e-01  5.99e+00  31.98s
   5  +6.76374779e+04  -8.77178644e+03   7.77e+04    1.24e-01  2.80e+00  34.96s
   6  +5.47063786e+04  +2.14338969e+04   3.37e+04    5.75e-02  1.08e+00  37.73s
   7  +5.11497402e+04  +2.99006339e+04   2.15e+04    4.20e-02  6.24e-01  40.31s
   8  +4.69504087e+04  +3.57171024e+04   1.13e+04    2.28e-02  3.05e-01  43.07s
   9  +4.48897023e+04  +3.80587025e+04   6.88e+03    1.38e-02  1.85e-01  45.68s
  10  +4.45573602e+04  +3.88240627e+04   5.77e+03    1.23e-02  1.48e-01  48.20s
  11  +4.31919111e+04  +4.03641692e+04   2.84e+03    5.75e-03  7.25e-02  50.90s
  12  +4.25249991e+04  +4.13820700e+04   1.15e+03    2.52e-03  2.42e-02  53.93s
  13  +4.21968343e+04  +4.16899379e+04   5.09e+02    1.06e-03  1.03e-02  56.74s
  14  +4.20201086e+04  +4.18738306e+04   1.47e+02    2.80e-04  2.79e-03  59.76s
  15  +4.19956492e+04  +4.19029428e+04   9.30e+01    1.81e-04  1.73e-03  62.30s
  16  +4.19603624e+04  +4.19413862e+04   1.90e+01    3.60e-05  3.04e-04  65.31s
  17  +4.19527894e+04  +4.19491700e+04   3.63e+00    7.10e-06  4.88e-05  68.31s
  18  +4.19519892e+04  +4.19499133e+04   2.08e+00    4.18e-06  2.72e-05  70.81s
  19  +4.19510585e+04  +4.19506456e+04   4.14e-01    7.27e-07  5.90e-06  73.72s
  20  +4.19508865e+04  +4.19508330e+04   5.36e-02    8.65e-08  7.75e-07  76.74s
  21  +4.19508736e+04  +4.19508441e+04   2.96e-02    4.07e-08  4.80e-07  79.34s
  22  +4.19508687e+04  +4.19508558e+04   1.29e-02    2.31e-08  1.68e-07  81.92s
  23  +4.19508649e+04  +4.19508586e+04   6.33e-03    1.10e-08  8.93e-08  84.48s
  24  +4.19508620e+04  +4.19508616e+04   4.04e-04    6.78e-09  6.45e-09  87.32s

Barrier status:                  OPTIMAL
Primal objective:                4.19508620e+04
Dual objective:                  4.19508616e+04
Duality gap (abs/rel):           4.03e-04 / 9.61e-09
Primal infeasibility (abs/rel):  6.78e-09 / 7.53e-10
Dual infeasibility (abs/rel):    6.45e-09 / 1.68e-13
Postsolving

Solving finished
Status: Optimal  Objective: 4.1950862012e+04  Iterations: 24(0)  Time: 87.50s
COPT> set gpumode 1
Setting parameter 'GPUMode' to 1
COPT> opt
Model fingerprint: 11cebb4f

Using Cardinal Optimizer v8.0.1 on Windows
Hardware has 8 cores and 16 threads. Using instruction set X86_NATIVE (1)
Minimizing an LP problem

The original problem has:
    14135 rows, 84693 columns and 11800327 non-zero elements
The presolved problem has:
    13108 rows, 81336 columns and 11541827 non-zero elements

Hardware has 1 supported GPU device with CUDA 12.8
  GPU 0: NVIDIA GeForce RTX 3080 (CUDA capability 8.6)

Starting barrier solver using 8 CPU threads and GPU 0

Problem info:
    Dualized in presolve:            No
    Range of matrix coefficients:    [6e-02,5e+01]
    Range of rhs coefficients:       [4e-01,7e+01]
    Range of bound coefficients:     [1e+00,1e+00]
    Range of cost coefficients:      [1e+00,4e+04]

Factor info:
    Number of free columns:          0
    Number of dense columns:         0
    Number of matrix entries:        6.485e+07
    Number of factor entries:        7.216e+07
    Number of factor flops:          5.209e+11

Iter       Primal.Obj         Dual.Obj      Compl  Primal.Inf  Dual.Inf    Time
   0  +6.89041521e+07  -1.47285422e+08   5.11e+08    2.93e+03  3.86e+04  20.43s
   1  +2.06731774e+07  -4.05798355e+07   7.96e+07    2.95e+02  4.35e+03  22.78s
   2  +1.51129327e+06  -2.46845951e+06   4.38e+06    1.09e+01  1.98e+02  25.28s
   3  +2.08981694e+05  -3.21353877e+05   5.55e+05    1.01e+00  2.23e+01  27.87s
   4  +1.50797642e+05  -1.91581467e+05   3.55e+05    6.55e-01  1.43e+01  30.14s
   5  +7.69542411e+04  -3.84123392e+04   1.18e+05    1.86e-01  4.59e+00  32.55s
   6  +5.72886861e+04  +1.30200724e+04   4.49e+04    7.17e-02  1.59e+00  34.95s
   7  +5.28199307e+04  +2.48832838e+04   2.83e+04    5.10e-02  9.35e-01  37.21s
   8  +4.93139363e+04  +3.14684082e+04   1.80e+04    3.57e-02  5.57e-01  39.46s
   9  +4.64869445e+04  +3.53299716e+04   1.13e+04    2.24e-02  3.45e-01  41.75s
  10  +4.49554866e+04  +3.72739519e+04   7.74e+03    1.51e-02  2.41e-01  44.03s
  11  +4.42270257e+04  +3.85967370e+04   5.67e+03    1.13e-02  1.74e-01  46.28s
  12  +4.34448741e+04  +3.97930041e+04   3.67e+03    7.37e-03  1.12e-01  48.55s
  13  +4.31180677e+04  +4.02015265e+04   2.93e+03    5.75e-03  9.12e-02  50.78s
  14  +4.25586573e+04  +4.10818743e+04   1.48e+03    2.85e-03  4.42e-02  53.20s
  15  +4.24566883e+04  +4.13332744e+04   1.13e+03    2.36e-03  3.17e-02  55.43s
  16  +4.21090974e+04  +4.17650782e+04   3.45e+02    7.22e-04  8.69e-03  58.03s
  17  +4.20776772e+04  +4.18154490e+04   2.63e+02    5.81e-04  6.30e-03  60.29s
  18  +4.20175668e+04  +4.18758590e+04   1.42e+02    3.05e-04  3.40e-03  62.64s
  19  +4.19982883e+04  +4.18983411e+04   1.00e+02    2.16e-04  2.39e-03  64.88s
  20  +4.19592054e+04  +4.19397701e+04   1.95e+01    3.48e-05  4.62e-04  67.39s
  21  +4.19527727e+04  +4.19480929e+04   4.69e+00    7.30e-06  1.13e-04  69.84s
  22  +4.19511738e+04  +4.19504306e+04   7.45e-01    1.06e-06  1.75e-05  72.33s
  23  +4.19510182e+04  +4.19506430e+04   3.76e-01    5.26e-07  8.95e-06  74.61s
  24  +4.19508776e+04  +4.19508442e+04   3.35e-02    4.36e-08  7.53e-07  77.14s
  25  +4.19508632e+04  +4.19508604e+04   2.87e-03    1.74e-08  6.06e-08  79.71s
  26  +4.19508622e+04  +4.19508614e+04   8.25e-04    3.16e-08  1.97e-08  82.09s
  27  +4.19508619e+04  +4.19508615e+04   4.11e-04    2.53e-08  1.29e-08  84.39s

Barrier status:                  OPTIMAL
Primal objective:                4.19508619e+04
Dual objective:                  4.19508615e+04
Duality gap (abs/rel):           4.15e-04 / 9.90e-09
Primal infeasibility (abs/rel):  2.53e-08 / 2.81e-09
Dual infeasibility (abs/rel):    1.29e-08 / 3.37e-13
Postsolving

Solving finished
Status: Optimal  Objective: 4.1950861950e+04  Iterations: 27(0)  Time: 84.64s
COPT> set lpmethod 6
Setting parameter 'LpMethod' to 6
COPT> opt
Model fingerprint: 11cebb4f

Using Cardinal Optimizer v8.0.1 on Windows
Hardware has 8 cores and 16 threads. Using instruction set X86_NATIVE (1)
Minimizing an LP problem

The original problem has:
    14135 rows, 84693 columns and 11800327 non-zero elements
The presolved problem has:
    13108 rows, 81336 columns and 11541827 non-zero elements

Hardware has 1 supported GPU device with CUDA 12.8
  GPU 0: NVIDIA GeForce RTX 3080 (CUDA capability 8.6)

Starting PDLP solver using GPU 0

Problem info:
    Range of matrix coefficients:    [3e-04,5e-01]
    Range of rhs coefficients:       [1e-01,4e+00]
    Range of bound coefficients:     [2e-01,1e+02]
    Range of cost coefficients:      [9e-03,2e+04]

Iterations       Primal.Obj         Dual.Obj        Gap  Primal.Inf  Dual.Inf    Time
         0  +0.00000000e+00  +0.00000000e+00  +0.00e+00    1.11e+04  0.00e+00   4.15s
      4000  +4.17407991e+04  +3.04897437e+04  +1.13e+04    2.98e+00  3.50e-01   6.63s
      8000  +4.19070214e+04  +4.08247925e+04  +1.08e+03    1.13e+00  3.70e-02   8.82s
     12000  +4.19414758e+04  +4.16358914e+04  +3.06e+02    4.70e-01  1.04e-02  11.00s
     16000  +4.19482878e+04  +4.18529461e+04  +9.53e+01    1.90e-01  3.07e-03  13.17s
     20000  +4.19509475e+04  +4.18808838e+04  +7.01e+01    1.31e-01  2.14e-03  15.36s
     24000  +4.19503849e+04  +4.18968156e+04  +5.36e+01    8.06e-02  1.69e-03  17.55s
     28000  +4.19513357e+04  +4.19022223e+04  +4.91e+01    6.54e-02  1.52e-03  19.72s
     32000  +4.19500141e+04  +4.19077636e+04  +4.23e+01    5.78e-02  1.34e-03  21.93s
     36000  +4.19503340e+04  +4.18930356e+04  +5.73e+01    4.79e-02  1.88e-03  24.13s
     40000  +4.19512214e+04  +4.18980661e+04  +5.32e+01    4.64e-02  1.95e-03  26.32s
     44000  +4.19517251e+04  +4.19044864e+04  +4.72e+01    4.24e-02  1.47e-03  28.51s
     48000  +4.19519873e+04  +4.19109102e+04  +4.11e+01    3.67e-02  1.31e-03  30.71s
     52000  +4.19507950e+04  +4.18687117e+04  +8.21e+01    4.05e-02  2.65e-03  32.93s
     56000  +4.19510083e+04  +4.18951665e+04  +5.58e+01    3.67e-02  3.33e-01  35.14s
     60000  +4.19508867e+04  +4.19031857e+04  +4.77e+01    2.64e-02  1.52e-03  37.36s
     64000  +4.19502261e+04  +4.19103260e+04  +3.99e+01    2.44e-02  1.25e-03  39.58s
     68000  +4.19505520e+04  +4.19165477e+04  +3.40e+01    2.03e-02  1.08e-03  41.80s
     72000  +4.19510967e+04  +4.19211335e+04  +3.00e+01    2.19e-02  9.25e-04  44.01s
     76000  +4.19509285e+04  +4.19249344e+04  +2.60e+01    1.64e-02  7.93e-04  46.22s
     80000  +4.19510760e+04  +4.19055577e+04  +4.55e+01    2.56e-02  1.48e-03  48.45s
     84000  +4.19506571e+04  +4.19217738e+04  +2.89e+01    1.92e-02  8.98e-04  50.66s
     88000  +4.19508379e+04  +4.19269419e+04  +2.39e+01    1.67e-02  7.64e-04  52.87s
     92000  +4.19508950e+04  +4.19301711e+04  +2.07e+01    1.52e-02  6.42e-04  55.08s
     96000  +4.19509516e+04  +4.19322505e+04  +1.87e+01    1.42e-02  5.79e-04  57.31s
    100000  +4.19510837e+04  +4.19340930e+04  +1.70e+01    1.26e-02  5.27e-04  59.51s
    104000  +4.19508114e+04  +4.19355736e+04  +1.52e+01    1.19e-02  4.77e-04  61.73s
    108000  +4.19508323e+04  +4.19366070e+04  +1.42e+01    1.02e-02  4.35e-04  63.95s
    112000  +4.19507764e+04  +4.19381225e+04  +1.27e+01    9.66e-03  4.19e-04  66.17s
    116000  +4.19510914e+04  +4.19389582e+04  +1.21e+01    1.01e-02  3.63e-04  68.39s
    120000  +4.19509229e+04  +4.19401513e+04  +1.08e+01    9.40e-03  3.40e-04  70.60s
    124000  +4.19508567e+04  +4.19402380e+04  +1.06e+01    7.96e-03  3.75e-04  72.82s
    128000  +4.19508746e+04  +4.19446876e+04  +6.19e+00    6.75e-03  1.96e-04  75.05s
    132000  +4.19507884e+04  +4.19456953e+04  +5.09e+00    6.51e-03  1.61e-04  77.28s
    136000  +4.19508935e+04  +4.19462747e+04  +4.62e+00    6.52e-03  1.41e-04  79.51s
    140000  +4.19508247e+04  +4.19466312e+04  +4.19e+00    6.01e-03  1.29e-04  81.72s
    144000  +4.19507159e+04  +4.19470167e+04  +3.70e+00    6.23e-03  1.21e-04  83.96s
    148000  +4.19508886e+04  +4.19471295e+04  +3.76e+00    4.95e-03  1.14e-04  86.19s
    152000  +4.19508882e+04  +4.19474299e+04  +3.46e+00    5.94e-03  1.06e-04  88.43s
    156000  +4.19507225e+04  +4.19474303e+04  +3.29e+00    4.48e-03  1.11e-04  90.66s
    160000  +4.19508571e+04  +4.19476056e+04  +3.25e+00    5.01e-03  9.95e-05  92.89s
    164000  +4.19509950e+04  +4.19476968e+04  +3.30e+00    4.79e-03  1.00e-04  95.13s
    168000  +4.19509529e+04  +4.19478663e+04  +3.09e+00    4.70e-03  9.31e-05  97.35s
    172000  +4.19508620e+04  +4.19479273e+04  +2.93e+00    3.99e-03  9.24e-05  99.58s
    176000  +4.19507600e+04  +4.19480613e+04  +2.70e+00    4.59e-03  7.03e-03    101s
    180000  +4.19509239e+04  +4.19481562e+04  +2.77e+00    4.17e-03  8.41e-05    104s
    184000  +4.19508979e+04  +4.19482315e+04  +2.67e+00    3.92e-03  8.42e-05    106s
    188000  +4.19508813e+04  +4.19483136e+04  +2.57e+00    4.21e-03  1.13e-04    108s
    192000  +4.19508572e+04  +4.19487905e+04  +2.07e+00    2.65e-03  7.23e-05    110s
    196000  +4.19508643e+04  +4.19493215e+04  +1.54e+00    2.82e-03  4.93e-05    112s
    200000  +4.19508545e+04  +4.19494812e+04  +1.37e+00    2.43e-03  4.27e-05    115s
    204000  +4.19508242e+04  +4.19496231e+04  +1.20e+00    2.41e-03  3.84e-05    117s
    208000  +4.19508278e+04  +4.19497299e+04  +1.10e+00    2.33e-03  3.53e-05    119s
    212000  +4.19508570e+04  +4.19497641e+04  +1.09e+00    2.14e-03  3.43e-05    121s
    216000  +4.19508813e+04  +4.19498278e+04  +1.05e+00    2.27e-03  3.17e-05    124s
    220000  +4.19508783e+04  +4.19498143e+04  +1.06e+00    2.03e-03  6.76e-05    126s
    224000  +4.19509001e+04  +4.19499121e+04  +9.88e-01    2.38e-03  2.90e-05    128s
    228000  +4.19508655e+04  +4.19499468e+04  +9.19e-01    2.06e-03  2.88e-05    130s
    232000  +4.19508844e+04  +4.19499751e+04  +9.09e-01    2.24e-03  2.73e-05    133s
    236000  +4.19508559e+04  +4.19499898e+04  +8.66e-01    2.06e-03  2.74e-05    135s
    240000  +4.19508181e+04  +4.19500076e+04  +8.10e-01    1.79e-03  2.68e-05    137s
    244000  +4.19507979e+04  +4.19500302e+04  +7.68e-01    2.04e-03  2.58e-05    139s
    248000  +4.19508312e+04  +4.19500505e+04  +7.81e-01    1.80e-03  2.54e-05    142s
    252000  +4.19508430e+04  +4.19500623e+04  +7.81e-01    1.68e-03  2.52e-05    144s
    256000  +4.19508648e+04  +4.19500766e+04  +7.88e-01    1.54e-03  2.44e-05    146s
    260000  +4.19508884e+04  +4.19500939e+04  +7.94e-01    1.87e-03  3.94e-05    148s
    264000  +4.19508957e+04  +4.19500869e+04  +8.09e-01    1.65e-03  2.38e-05    150s
    268000  +4.19508912e+04  +4.19501136e+04  +7.78e-01    1.64e-03  2.36e-05    153s
    272000  +4.19508724e+04  +4.19501426e+04  +7.30e-01    1.68e-03  2.28e-05    155s
    276000  +4.19508528e+04  +4.19499680e+04  +8.85e-01    1.60e-03  2.97e-05    157s
    280000  +4.19508641e+04  +4.19502685e+04  +5.96e-01    1.23e-03  1.88e-05    159s
    284000  +4.19508608e+04  +4.19503494e+04  +5.11e-01    1.08e-03  1.58e-05    162s
    288000  +4.19508616e+04  +4.19504107e+04  +4.51e-01    1.20e-03  1.39e-05    164s
    292000  +4.19508758e+04  +4.19504364e+04  +4.39e-01    1.08e-03  1.34e-05    166s
    293600  +4.19508627e+04  +4.19507789e+04  +8.38e-02    3.67e-04  2.53e-06    167s

PDLP status:                     OPTIMAL
PDLP iterations:                 293600
Primal objective:                4.19508627e+04
Dual objective:                  4.19507789e+04
Primal infeasibility (abs/rel):  3.67e-04 / 3.32e-08
Dual infeasibility (abs/rel):    2.53e-06 / 2.91e-11
Duality gap (abs/rel):           8.38e-02 / 9.99e-07

Postsolving

Solving finished
Status: Optimal  Objective: 4.1950862671e+04  Iterations: 293600(0)  Time: 167.63s


(.venv) root@dev1:/home/# cuopt_cli --method=1 --presolve=true --relative-primal-tolerance=1e-6 test/instance24.mps
Setting parameter method to 1
Setting parameter presolve to true
Setting parameter relative_primal_tolerance to 1.000000e-06
Reading file instance24.mps
cuOpt version: 25.10.1, git hash: 876fcfc, host arch: x86_64, device archs: 75-real,80-real,86-real,90a-real,100f-real,120a-real,120
CPU: AMD Ryzen 5 2600 Six-Core Processor, threads (physical/logical): 6/12, RAM: 5.71 GiB
CUDA 13.0, device: NVIDIA GeForce GTX 1660 (ID 0), VRAM: 6.00 GiB
CUDA device UUID: ffffff9bffffffbcffffffa225-ffffffa72

Solving a problem with 14135 constraints, 84693 variables (0 integers), and 11800327 nonzeros
Problem scaling:
Objective coefficents range:          [1e+00, 3e+02]
Constraint matrix coefficients range: [1e+00, 3e+03]
Constraint rhs / bounds range:        [1e+00, 1e+03]
Variable bounds range:                [0e+00, 1e+02]

Original problem: 14135 constraints, 84693 variables, 11800327 nonzeros
Calling Papilo presolver
Disabling the presolver methods that do not support dual postsolve
Presolve status: reduced the problem
Presolve removed: 1027 constraints, 13985 variables, 269128 nonzeros
Presolved problem: 13108 constraints, 70708 variables, 11531199 nonzeros
Papilo presolve time: 4.054019
Objective offset 42112200.000000 scaling_factor 1.000000
   Iter    Primal Obj.      Dual Obj.    Gap        Primal Res.  Dual Res.   Time
      0 +4.21122000e+07 +4.21122000e+07  0.00e+00   1.22e+01     7.36e+07   4.211s
   1000 +7.10231166e+04 +8.87823939e+04  1.78e+04   1.42e+01     2.36e+04   6.327s
   2000 +8.22431388e+04 +5.63061842e+04  2.59e+04   1.64e+00     5.21e+03   8.163s
   3000 +8.14735222e+04 +4.57217637e+04  3.58e+04   1.28e+00     1.37e+03   10.001s
   4000 +5.85322830e+04 +4.30805984e+04  1.55e+04   2.28e+00     5.00e+02   11.839s
   5000 +5.18175772e+04 +4.24257375e+04  9.39e+03   1.16e+00     2.30e+02   13.674s
   6000 +4.45638298e+04 +4.20909159e+04  2.47e+03   1.38e+00     1.10e+02   15.509s
   7000 +4.40521924e+04 +4.19842329e+04  2.07e+03   1.21e+00     4.65e+01   17.346s
   8000 +4.40021685e+04 +4.19670984e+04  2.04e+03   6.31e-01     2.77e+01   19.183s
   9000 +4.36845327e+04 +4.19595270e+04  1.73e+03   4.43e-01     2.02e+01   21.024s
  10000 +4.36066182e+04 +4.19466670e+04  1.66e+03   3.40e-01     1.34e+01   22.862s
  11000 +4.34340425e+04 +4.19473069e+04  1.49e+03   2.98e-01     1.10e+01   24.686s
  12000 +4.33650676e+04 +4.19487325e+04  1.42e+03   2.04e-01     1.08e+01   26.510s
  13000 +4.32280007e+04 +4.19498534e+04  1.28e+03   2.09e-01     1.01e+01   28.336s
  14000 +4.31635430e+04 +4.19502678e+04  1.21e+03   1.57e-01     9.38e+00   30.160s
  15000 +4.29279300e+04 +4.19528406e+04  9.75e+02   1.49e-01     1.43e+01   31.988s
  16000 +4.26113016e+04 +4.19548002e+04  6.57e+02   1.44e-01     2.02e+01   33.815s
  17000 +4.23855126e+04 +4.19542966e+04  4.31e+02   1.58e-01     2.21e+01   35.639s
  18000 +4.22919321e+04 +4.19555036e+04  3.36e+02   1.36e-01     2.10e+01   37.467s
  19000 +4.22261674e+04 +4.19541632e+04  2.72e+02   1.21e-01     1.87e+01   39.296s
  20000 +4.21221562e+04 +4.19535454e+04  1.69e+02   1.07e-01     1.61e+01   41.125s
  21000 +4.21004838e+04 +4.19532853e+04  1.47e+02   8.99e-02     1.41e+01   42.954s
  22000 +4.20546042e+04 +4.19533260e+04  1.01e+02   6.65e-02     1.20e+01   44.781s
  23000 +4.19625290e+04 +4.19515572e+04  1.10e+01   6.75e-02     9.74e+00   46.613s
  24000 +4.19716938e+04 +4.19517660e+04  1.99e+01   5.25e-02     9.30e+00   48.442s
  25000 +4.19650576e+04 +4.19516001e+04  1.35e+01   4.44e-02     8.43e+00   50.269s
  26000 +4.19744741e+04 +4.19512414e+04  2.32e+01   3.01e-02     7.38e+00   52.100s
  27000 +4.19498816e+04 +4.19510921e+04  1.21e+00   2.62e-02     6.18e+00   53.928s
  28000 +4.19202749e+04 +4.19510622e+04  3.08e+01   2.82e-02     5.27e+00   55.756s
  29000 +4.19200158e+04 +4.19510685e+04  3.11e+01   2.35e-02     4.64e+00   57.587s
  30000 +4.19339749e+04 +4.19510884e+04  1.71e+01   2.00e-02     4.13e+00   59.415s
  31000 +4.19372488e+04 +4.19510732e+04  1.38e+01   2.05e-02     3.65e+00   61.245s
  32000 +4.19450298e+04 +4.19509837e+04  5.95e+00   1.86e-02     3.32e+00   63.073s
  33000 +4.19377997e+04 +4.19509953e+04  1.32e+01   1.68e-02     3.08e+00   64.905s
  34000 +4.19116770e+04 +4.19509917e+04  3.93e+01   2.92e-02     2.66e+00   66.736s
  35000 +4.19171306e+04 +4.19510260e+04  3.39e+01   3.93e-02     2.08e+00   68.563s
  36000 +4.19018061e+04 +4.19509506e+04  4.91e+01   4.33e-02     1.68e+00   70.392s
  37000 +4.19288702e+04 +4.19509034e+04  2.20e+01   4.23e-02     1.34e+00   72.220s
  38000 +4.19244557e+04 +4.19508715e+04  2.64e+01   4.03e-02     1.10e+00   74.051s
  39000 +4.19356704e+04 +4.19508763e+04  1.52e+01   3.77e-02     9.55e-01   75.887s
  40000 +4.19165846e+04 +4.19508814e+04  3.43e+01   3.66e-02     8.39e-01   77.722s
  41000 +4.19203713e+04 +4.19508940e+04  3.05e+01   3.27e-02     7.23e-01   79.555s
  42000 +4.19118138e+04 +4.19508570e+04  3.90e+01   2.82e-02     6.40e-01   81.390s
  43000 +4.19373738e+04 +4.19508528e+04  1.35e+01   2.25e-02     5.75e-01   83.226s
  44000 +4.19374643e+04 +4.19508963e+04  1.34e+01   2.25e-02     5.16e-01   85.059s
  45000 +4.19347222e+04 +4.19509136e+04  1.62e+01   1.94e-02     4.71e-01   86.897s
  46000 +4.19361199e+04 +4.19508517e+04  1.47e+01   1.72e-02     4.46e-01   88.731s
  47000 +4.19325144e+04 +4.19508694e+04  1.84e+01   1.76e-02     4.14e-01   90.565s
  48000 +4.19313102e+04 +4.19508622e+04  1.96e+01   1.68e-02     3.63e-01   92.401s
  49000 +4.19345851e+04 +4.19508529e+04  1.63e+01   1.51e-02     3.24e-01   94.236s
  50000 +4.19391543e+04 +4.19508554e+04  1.17e+01   1.36e-02     3.02e-01   96.069s
  51000 +4.19414688e+04 +4.19508667e+04  9.40e+00   1.21e-02     2.89e-01   97.903s
  52000 +4.19430661e+04 +4.19508659e+04  7.80e+00   1.12e-02     2.79e-01   99.736s
  53000 +4.19416350e+04 +4.19508654e+04  9.23e+00   1.05e-02     2.68e-01   101.570s
  53800 +4.19428379e+04 +4.19508666e+04  8.03e+00   9.98e-03     2.58e-01   103.037s
LP Solver status:                Optimal
Primal objective:                +4.19428379e+04
Dual objective:                  +4.19508666e+04
Duality gap (abs/rel):           +8.03e+00 / +9.57e-05
Primal infeasibility (abs/rel):  +9.98e-03 / +9.02e-07
Dual infeasibility (abs/rel):    +2.58e-01 / +3.51e-09
PDLP finished
Status: Optimal   Objective: 4.19428379e+04  Iterations: 53800  Time: 98.968s, Total time 103.042s
Post-solve status: Post solved solution violates constraints. This is most likely due to different tolerances.
(.venv) root@dev1:/home/#

2025年11月29日土曜日

COPTとHighsをSubProblem用の求解機として使う

 上記を検討してみました。(COPTについては、前回のOSI実装を使っていますが、一応Simplexを指定しています。)

インスタンス24での解き始めの状態での比較になります。

下記がそのログになりますが、意外なことに、このインスタンスでは、Highsの方が5倍程度速いという結果になりました。どちらにしても、遅すぎて使えません、ということは変わりがありません。予想された結果ではありますが、これ速度では1年走らせても結果は得られないでしょう。(1週間程度のマシンRunningで厳密解が求まることを目指しています)

問題点としては、

1)1Rosterあたり1000秒以上かかっている。(目標10秒程度に対して遅すぎる)

2)マルチスレッドが効かない(CPU使用率でみると1スレッド分程度しか動いていない)

が挙げられます。マルチスレッドが効かないのは恐らくメモリ消費が3次キャッシュ容量を上回るからだと思われます。3次キャッシュがさらに大きなCPU(Ryzen 3D)を使えば、改善が期待できますが、現在でも32MBあり、最大でも4倍程度しか期待できません。COPTに僅かな期待をしてだけに残念な結果です。

いずれにせよ、既存のアプローチでは、現実的な時間内には求まらない、ということです。最難関の最終問題に相応しいインスタンスとなっています。


		COPT Barrier is used.
		2025-11-27 17:35:32: Roster::Highs MipSolver.
	Elapsed separation time 1657(msec) threads=1
		Loading time= 0(sec)
Load Problem
		subproblem[0]  person=0 obj=25.000000 10001(msec)
		subproblem[1]  person=1 obj=20.000000 7894(msec)
		subproblem[2]  person=2 obj=13.000000 5633(msec)
		subproblem[3]  person=3 obj=18.000000 6261(msec)
		subproblem[4]  person=4 obj=33.000000 12388(msec)
		subproblem[5]  person=5 obj=3.000000 5155(msec)
		subproblem[6]  person=6 obj=20.000000 9581(msec)
		subproblem[7]  person=7 obj=3.000000 6529(msec)
		subproblem[8]  person=8 obj=9.000000 4178(msec)
		subproblem[9]  person=9 obj=10.000000 5616(msec)
		subproblem[10]  person=10 obj=38.000000 5690(msec)
		subproblem[11]  person=11 obj=18.000000 7879(msec)
		subproblem[12]  person=12 obj=28.000000 11218(msec)
		subproblem[13]  person=13 obj=26.000000 10218(msec)
		subproblem[14]  person=14 obj=21.000000 6852(msec)
		subproblem[15]  person=15 obj=11.000000 5137(msec)
		subproblem[16]  person=16 obj=30.000000 27560(msec)
		subproblem[17]  person=17 obj=17.000000 9766(msec)
		subproblem[18]  person=18 obj=31.000000 4360(msec)
		subproblem[19]  person=19 obj=15.000000 5875(msec)
		subproblem[20]  person=20 obj=9.000000 8356(msec)
		subproblem[21]  person=21 obj=24.000000 8414(msec)
		subproblem[22]  person=22 obj=9.000000 10560(msec)
		subproblem[23]  person=23 obj=14.000000 6596(msec)
		subproblem[24]  person=24 obj=8.000000 6775(msec)
		subproblem[25]  person=25 obj=23.000000 6456(msec)
		subproblem[26]  person=26 obj=10.000000 11239(msec)
		subproblem[27]  person=27 obj=31.000000 7665(msec)
		subproblem[28]  person=28 obj=15.000000 11035(msec)
		subproblem[29]  person=29 obj=35.000000 21868(msec)
		subproblem[30]  person=30 obj=11.000000 4637(msec)
		subproblem[31]  person=31 obj=28.000000 6491(msec)
		subproblem[32]  person=32 obj=9.000000 15341(msec)
		subproblem[33]  person=33 obj=14.000000 5671(msec)
		subproblem[34]  person=34 obj=8.000000 4791(msec)
		subproblem[35]  person=35 obj=19.000000 10790(msec)
		subproblem[36]  person=36 obj=23.000000 6463(msec)
		subproblem[37]  person=37 obj=13.000000 7081(msec)
		subproblem[38]  person=38 obj=14.000000 10974(msec)
		subproblem[39]  person=39 obj=19.000000 4303(msec)
		subproblem[40]  person=40 obj=29.000000 5925(msec)
		subproblem[41]  person=41 obj=28.000000 4107(msec)
		subproblem[42]  person=42 obj=21.000000 6398(msec)
		subproblem[43]  person=43 obj=0.000000 7649(msec)
		subproblem[44]  person=44 obj=29.000000 8120(msec)
		subproblem[45]  person=45 obj=13.000000 11866(msec)
		subproblem[46]  person=46 obj=14.000000 6839(msec)
		subproblem[47]  person=47 obj=4.000000 7414(msec)
		subproblem[48]  person=48 obj=24.000000 9587(msec)
		subproblem[49]  person=49 obj=9.000000 6736(msec)
		subproblem[50]  person=50 obj=26.000000 6741(msec)
		subproblem[51]  person=51 obj=23.000000 7883(msec)
		subproblem[52]  person=52 obj=17.000000 10404(msec)
		subproblem[53]  person=53 obj=5.000000 7677(msec)
		subproblem[54]  person=54 obj=19.000000 3667(msec)
		subproblem[55]  person=55 obj=31.000000 12995(msec)
		subproblem[56]  person=56 obj=15.000000 6512(msec)
		subproblem[57]  person=57 obj=21.000000 6353(msec)
		subproblem[58]  person=58 obj=8.000000 5703(msec)
		subproblem[59]  person=59 obj=11.000000 4727(msec)
		subproblem[60]  person=60 obj=9.000000 5728(msec)
		subproblem[61]  person=61 obj=22.000000 6158(msec)
		subproblem[62]  person=62 obj=17.000000 7775(msec)
		subproblem[63]  person=63 obj=10.000000 6722(msec)
		subproblem[64]  person=64 obj=16.000000 6836(msec)
		subproblem[65]  person=65 obj=5.000000 11161(msec)
		subproblem[66]  person=66 obj=9.000000 8700(msec)
		subproblem[67]  person=67 obj=25.000000 4381(msec)
		subproblem[68]  person=68 obj=11.000000 9666(msec)
		subproblem[69]  person=69 obj=34.000000 12477(msec)
		subproblem[70]  person=70 obj=12.000000 10325(msec)
		subproblem[71]  person=71 obj=10.000000 10009(msec)
		subproblem[72]  person=72 obj=9.000000 7502(msec)
		subproblem[73]  person=73 obj=5.000000 4069(msec)
		subproblem[74]  person=74 obj=18.000000 6886(msec)
		subproblem[75]  person=75 obj=3.000000 7277(msec)
		subproblem[76]  person=76 obj=36.000000 6856(msec)
		subproblem[77]  person=77 obj=9.000000 4452(msec)
		subproblem[78]  person=78 obj=6.000000 5786(msec)
		subproblem[79]  person=79 obj=7.000000 6621(msec)
		subproblem[80]  person=80 obj=22.000000 6587(msec)
		subproblem[81]  person=81 obj=17.000000 10949(msec)
		subproblem[82]  person=82 obj=14.000000 12060(msec)
		subproblem[83]  person=83 obj=17.000000 4703(msec)
		subproblem[84]  person=84 obj=6.000000 4434(msec)
		subproblem[85]  person=85 obj=18.000000 4609(msec)
		subproblem[86]  person=86 obj=19.000000 5494(msec)
		subproblem[87]  person=87 obj=6.000000 7838(msec)
		subproblem[88]  person=88 obj=14.000000 9134(msec)
		subproblem[89]  person=89 obj=31.000000 6052(msec)
		subproblem[90]  person=90 obj=16.000000 8831(msec)
		subproblem[91]  person=91 obj=17.000000 9203(msec)
		subproblem[92]  person=92 obj=29.000000 12865(msec)
		subproblem[93]  person=93 obj=26.000000 4516(msec)
		subproblem[94]  person=94 obj=32.000000 11168(msec)
		subproblem[95]  person=95 obj=11.000000 9341(msec)
		subproblem[96]  person=96 obj=32.000000 15117(msec)
		subproblem[97]  person=97 obj=33.000000 5997(msec)
		subproblem[98]  person=98 obj=20.000000 26592(msec)
		subproblem[99]  person=99 obj=20.000000 5890(msec)
		subproblem[100]  person=100 obj=7.000000 10616(msec)
		subproblem[101]  person=101 obj=12.000000 4615(msec)
		subproblem[102]  person=102 obj=12.000000 12695(msec)
		subproblem[103]  person=103 obj=19.000000 6504(msec)
		subproblem[104]  person=104 obj=27.000000 10119(msec)
		subproblem[105]  person=105 obj=16.000000 7708(msec)
		subproblem[106]  person=106 obj=2.000000 5476(msec)
		subproblem[107]  person=107 obj=22.000000 12512(msec)
		subproblem[108]  person=108 obj=40.000000 9899(msec)
		subproblem[109]  person=109 obj=13.000000 9304(msec)
		subproblem[110]  person=110 obj=35.000000 3791(msec)
		subproblem[111]  person=111 obj=28.000000 6512(msec)
		subproblem[112]  person=112 obj=21.000000 6191(msec)
		subproblem[113]  person=113 obj=34.000000 4444(msec)
		subproblem[114]  person=114 obj=19.000000 7856(msec)
		subproblem[115]  person=115 obj=15.000000 4704(msec)
		subproblem[116]  person=116 obj=39.000000 7727(msec)
		subproblem[117]  person=117 obj=49.000000 4418(msec)
		subproblem[118]  person=118 obj=33.000000 6544(msec)
		subproblem[119]  person=119 obj=16.000000 5703(msec)
		subproblem[120]  person=120 obj=15.000000 2734(msec)
		subproblem[121]  person=121 obj=16.000000 4920(msec)
		subproblem[122]  person=122 obj=20.000000 4625(msec)
		subproblem[123]  person=123 obj=17.000000 5508(msec)
		subproblem[124]  person=124 obj=17.000000 7560(msec)
		subproblem[125]  person=125 obj=14.000000 3940(msec)
		subproblem[126]  person=126 obj=3.000000 3658(msec)
		subproblem[127]  person=127 obj=23.000000 4504(msec)
		subproblem[128]  person=128 obj=17.000000 7188(msec)
		subproblem[129]  person=129 obj=21.000000 4372(msec)
		subproblem[130]  person=130 obj=13.000000 4677(msec)
		subproblem[131]  person=131 obj=25.000000 4635(msec)
		subproblem[132]  person=132 obj=18.000000 5054(msec)
		subproblem[133]  person=133 obj=18.000000 5787(msec)
		subproblem[134]  person=134 obj=9.000000 5279(msec)
		subproblem[135]  person=135 obj=8.000000 11605(msec)
		subproblem[136]  person=136 obj=6.000000 7623(msec)
		subproblem[137]  person=137 obj=17.000000 7733(msec)
		subproblem[138]  person=138 obj=20.000000 3769(msec)
		subproblem[139]  person=139 obj=18.000000 6162(msec)
		subproblem[140]  person=140 obj=45.000000 1560(msec)
		subproblem[141]  person=141 obj=61.000000 1319(msec)
		subproblem[142]  person=142 obj=45.000000 1224(msec)
		subproblem[143]  person=143 obj=70.000000 1556(msec)
		subproblem[144]  person=144 obj=63.000000 980(msec)
		subproblem[145]  person=145 obj=52.000000 1006(msec)
		subproblem[146]  person=146 obj=61.000000 1212(msec)
		subproblem[147]  person=147 obj=19.000000 4651(msec)
		subproblem[148]  person=148 obj=34.000000 2927(msec)
		subproblem[149]  person=149 obj=19.000000 4289(msec)
	*********UB=839642(0)  1237.04(cpu sec)
		 scale=0.0001
		ov=3048.54 ext=150 depth=0 LB=1.79769e+308 UB=839642 2234.011000(cpu sec) R_calc_time=996.539204(cpu sec)
		ov=3040.83 ext=300 depth=0 LB=1.79769e+308 UB=839642 3236.103000(cpu sec) R_calc_time=1001.626574(cpu sec)
		ov=3035.14 ext=450 depth=0 LB=1.79769e+308 UB=839642 4493.486000(cpu sec) R_calc_time=1256.906835(cpu sec)
		ov=3031.17 ext=600 depth=0 LB=1.79769e+308 UB=839642 5552.075000(cpu sec) R_calc_time=1058.024451(cpu sec)
		ov=3027.87 ext=750 depth=0 LB=1.79769e+308 UB=839642 6713.695000(cpu sec) R_calc_time=1161.034803(cpu sec)
		ov=3025.16 ext=900 depth=0 LB=1.79769e+308 UB=839642 7842.875000(cpu sec) R_calc_time=1128.522348(cpu sec)
		ov=3022.98 ext=1050 depth=0 LB=1.79769e+308 UB=839642 9062.871000(cpu sec) R_calc_time=1219.282774(cpu sec)
		ov=3020.72 ext=1199 depth=0 LB=1.79769e+308 UB=839642 10465.942000(cpu sec) R_calc_time=1402.321249(cpu sec)
		ov=3018.66 ext=1347 depth=0 LB=1.79769e+308 UB=839642 11975.743000(cpu sec) R_calc_time=1508.981705(cpu sec)
		ov=3017.02 ext=1496 depth=0 LB=1.79769e+308 UB=839642 13784.329000(cpu sec) R_calc_time=1807.010644(cpu sec)
		ov=3015.07 ext=1645 depth=0 LB=1.79769e+308 UB=839642 15227.393000(cpu sec) R_calc_time=1442.084932(cpu sec)
		ov=3013.5 ext=1794 depth=0 LB=1.79769e+308 UB=839642 16704.044000(cpu sec) R_calc_time=1475.661454(cpu sec)
		ov=3011.94 ext=1944 depth=0 LB=1.79769e+308 UB=839642 18026.764000(cpu sec) R_calc_time=1321.644897(cpu sec)
		ov=3010.72 ext=2094 depth=0 LB=1.79769e+308 UB=839642 19267.014000(cpu sec) R_calc_time=1239.082048(cpu sec)
		ov=3009.51 ext=2244 depth=0 LB=1.79769e+308 UB=839642 20821.664000(cpu sec) R_calc_time=1553.365347(cpu sec)
		ov=3008.39 ext=2394 depth=0 LB=1.79769e+308 UB=839642 22056.175000(cpu sec) R_calc_time=1233.069065(cpu sec)
		ov=3007.44 ext=2541 depth=0 LB=1.79769e+308 UB=839642 24252.082000(cpu sec) R_calc_time=2194.327395(cpu sec)
		ov=3006.48 ext=2691 depth=0 LB=1.79769e+308 UB=839642 25847.171000(cpu sec) R_calc_time=1593.460037(cpu sec)
		ov=3005.57 ext=2840 depth=0 LB=1.79769e+308 UB=839642 27241.316000(cpu sec) R_calc_time=1392.388769(cpu sec)
		ov=3004.78 ext=2989 depth=0 LB=1.79769e+308 UB=839642 28661.593000(cpu sec) R_calc_time=1416.340474(cpu sec)
		ov=3004.11 ext=3139 depth=0 LB=1.79769e+308 UB=839642 29974.232000(cpu sec) R_calc_time=1310.504639(cpu sec)
		ov=3003.46 ext=3289 depth=0 LB=1.79769e+308 UB=839642 31712.428000(cpu sec) R_calc_time=1735.905981(cpu sec)
		ov=3002.87 ext=3439 depth=0 LB=1.79769e+308 UB=839642 33545.185000(cpu sec) R_calc_time=1830.363651(cpu sec)
	中止要求を受け付けました。
コンパイルの準備中
ソルバを呼び出し中です。
解探索が終了しました。 34452 (秒)
解を見つけることが出来ませんでした。
コンパイルの準備中
ソルバを呼び出し中です。
	制約をコンパイル中です。
	Python プロパティファイルの生成が終わりました。
		COPT Barrier is used.
		2025-11-28 03:09:59: Roster::COPT MIP Solver.
	Elapsed separation time 1654(msec) threads=1
		Loading time= 0(sec)
Load Problem
		subproblem[0]  person=0 obj=25.000000 29306(msec)
		subproblem[1]  person=1 obj=20.000000 36191(msec)
		subproblem[2]  person=2 obj=13.000000 38065(msec)
		subproblem[3]  person=3 obj=18.000000 31776(msec)
		subproblem[4]  person=4 obj=33.000000 35095(msec)
		subproblem[5]  person=5 obj=3.000000 37539(msec)
		subproblem[6]  person=6 obj=20.000000 37131(msec)
		subproblem[7]  person=7 obj=3.000000 36146(msec)
		subproblem[8]  person=8 obj=9.000000 36742(msec)
		subproblem[9]  person=9 obj=10.000000 35021(msec)
		subproblem[10]  person=10 obj=38.000000 21270(msec)
		subproblem[11]  person=11 obj=18.000000 45688(msec)
		subproblem[12]  person=12 obj=28.000000 38858(msec)
		subproblem[13]  person=13 obj=26.000000 36826(msec)
		subproblem[14]  person=14 obj=21.000000 37010(msec)
		subproblem[15]  person=15 obj=11.000000 34413(msec)
		subproblem[16]  person=16 obj=30.000000 29329(msec)
		subproblem[17]  person=17 obj=17.000000 34596(msec)
		subproblem[18]  person=18 obj=31.000000 40679(msec)
		subproblem[19]  person=19 obj=15.000000 38274(msec)
		subproblem[20]  person=20 obj=9.000000 39720(msec)
		subproblem[21]  person=21 obj=24.000000 37629(msec)
		subproblem[22]  person=22 obj=9.000000 40969(msec)
		subproblem[23]  person=23 obj=14.000000 37951(msec)
		subproblem[24]  person=24 obj=8.000000 33056(msec)
		subproblem[25]  person=25 obj=23.000000 36085(msec)
		subproblem[26]  person=26 obj=10.000000 34793(msec)
		subproblem[27]  person=27 obj=31.000000 37320(msec)
		subproblem[28]  person=28 obj=15.000000 24469(msec)
		subproblem[29]  person=29 obj=35.000000 60948(msec)
		subproblem[30]  person=30 obj=11.000000 37137(msec)
		subproblem[31]  person=31 obj=28.000000 38854(msec)
		subproblem[32]  person=32 obj=9.000000 39115(msec)
		subproblem[33]  person=33 obj=14.000000 40320(msec)
		subproblem[34]  person=34 obj=8.000000 40818(msec)
		subproblem[35]  person=35 obj=19.000000 38553(msec)
		subproblem[36]  person=36 obj=23.000000 35170(msec)
		subproblem[37]  person=37 obj=13.000000 12698(msec)
		subproblem[38]  person=38 obj=14.000000 40098(msec)
		subproblem[39]  person=39 obj=19.000000 39767(msec)
		subproblem[40]  person=40 obj=29.000000 69876(msec)
		subproblem[41]  person=41 obj=28.000000 42457(msec)
		subproblem[42]  person=42 obj=21.000000 69665(msec)
		subproblem[43]  person=43 obj=0.000000 39870(msec)
		subproblem[44]  person=44 obj=29.000000 58514(msec)
		subproblem[45]  person=45 obj=13.000000 39546(msec)
		subproblem[46]  person=46 obj=14.000000 47153(msec)
		subproblem[47]  person=47 obj=4.000000 47897(msec)
		subproblem[48]  person=48 obj=24.000000 23603(msec)
		subproblem[49]  person=49 obj=9.000000 42054(msec)
		subproblem[50]  person=50 obj=26.000000 45192(msec)
		subproblem[51]  person=51 obj=23.000000 55105(msec)
		subproblem[52]  person=52 obj=17.000000 56554(msec)
		subproblem[53]  person=53 obj=5.000000 14245(msec)
		subproblem[54]  person=54 obj=19.000000 46160(msec)
		subproblem[55]  person=55 obj=31.000000 46660(msec)
		subproblem[56]  person=56 obj=15.000000 44011(msec)
		subproblem[57]  person=57 obj=21.000000 29265(msec)
		subproblem[58]  person=58 obj=8.000000 38043(msec)
		subproblem[59]  person=59 obj=11.000000 39295(msec)
		subproblem[60]  person=60 obj=9.000000 40919(msec)
		subproblem[61]  person=61 obj=22.000000 61843(msec)
		subproblem[62]  person=62 obj=17.000000 47783(msec)
		subproblem[63]  person=63 obj=10.000000 22808(msec)
		subproblem[64]  person=64 obj=16.000000 49058(msec)
		subproblem[65]  person=65 obj=5.000000 21661(msec)
		subproblem[66]  person=66 obj=9.000000 63697(msec)
		subproblem[67]  person=67 obj=25.000000 48226(msec)
		subproblem[68]  person=68 obj=11.000000 53102(msec)
		subproblem[69]  person=69 obj=34.000000 88390(msec)
		subproblem[70]  person=70 obj=12.000000 49690(msec)
		subproblem[71]  person=71 obj=10.000000 43996(msec)
		subproblem[72]  person=72 obj=9.000000 44293(msec)
		subproblem[73]  person=73 obj=5.000000 37455(msec)
		subproblem[74]  person=74 obj=18.000000 51609(msec)
		subproblem[75]  person=75 obj=3.000000 12732(msec)
		subproblem[76]  person=76 obj=36.000000 69722(msec)
		subproblem[77]  person=77 obj=9.000000 8186(msec)
		subproblem[78]  person=78 obj=6.000000 47858(msec)
		subproblem[79]  person=79 obj=7.000000 40488(msec)
		subproblem[80]  person=80 obj=22.000000 64778(msec)
		subproblem[81]  person=81 obj=17.000000 66153(msec)
		subproblem[82]  person=82 obj=14.000000 55547(msec)
		subproblem[83]  person=83 obj=17.000000 46231(msec)
		subproblem[84]  person=84 obj=6.000000 26654(msec)
		subproblem[85]  person=85 obj=18.000000 44884(msec)
		subproblem[86]  person=86 obj=19.000000 52952(msec)
		subproblem[87]  person=87 obj=6.000000 16447(msec)
		subproblem[88]  person=88 obj=14.000000 73343(msec)
		subproblem[89]  person=89 obj=31.000000 74674(msec)
		subproblem[90]  person=90 obj=16.000000 34939(msec)
		subproblem[91]  person=91 obj=17.000000 42085(msec)
		subproblem[92]  person=92 obj=29.000000 64338(msec)
		subproblem[93]  person=93 obj=26.000000 64250(msec)
		subproblem[94]  person=94 obj=32.000000 47622(msec)
		subproblem[95]  person=95 obj=11.000000 27532(msec)
		subproblem[96]  person=96 obj=32.000000 49902(msec)
		subproblem[97]  person=97 obj=33.000000 55576(msec)
		subproblem[98]  person=98 obj=20.000000 44937(msec)
		subproblem[99]  person=99 obj=20.000000 41679(msec)
		subproblem[100]  person=100 obj=7.000000 41017(msec)
		subproblem[101]  person=101 obj=12.000000 67954(msec)
		subproblem[102]  person=102 obj=12.000000 60802(msec)
		subproblem[103]  person=103 obj=19.000000 16946(msec)
		subproblem[104]  person=104 obj=27.000000 62952(msec)
		subproblem[105]  person=105 obj=16.000000 66402(msec)
		subproblem[106]  person=106 obj=2.000000 12460(msec)
		subproblem[107]  person=107 obj=22.000000 57913(msec)
		subproblem[108]  person=108 obj=40.000000 68054(msec)
		subproblem[109]  person=109 obj=13.000000 34645(msec)
		subproblem[110]  person=110 obj=35.000000 56692(msec)
		subproblem[111]  person=111 obj=28.000000 76231(msec)
		subproblem[112]  person=112 obj=21.000000 50606(msec)
		subproblem[113]  person=113 obj=34.000000 70873(msec)
		subproblem[114]  person=114 obj=19.000000 60110(msec)
		subproblem[115]  person=115 obj=15.000000 50182(msec)
		subproblem[116]  person=116 obj=39.000000 69674(msec)
		subproblem[117]  person=117 obj=49.000000 74751(msec)
		subproblem[118]  person=118 obj=33.000000 58608(msec)
		subproblem[119]  person=119 obj=16.000000 105830(msec)
		subproblem[120]  person=120 obj=15.000000 55606(msec)
		subproblem[121]  person=121 obj=16.000000 55364(msec)
		subproblem[122]  person=122 obj=20.000000 64397(msec)
		subproblem[123]  person=123 obj=17.000000 51171(msec)
		subproblem[124]  person=124 obj=17.000000 47927(msec)
		subproblem[125]  person=125 obj=14.000000 17714(msec)
		subproblem[126]  person=126 obj=3.000000 24450(msec)
		subproblem[127]  person=127 obj=23.000000 50607(msec)
		subproblem[128]  person=128 obj=17.000000 38429(msec)
		subproblem[129]  person=129 obj=21.000000 60636(msec)
		subproblem[130]  person=130 obj=13.000000 46441(msec)
		subproblem[131]  person=131 obj=25.000000 58358(msec)
		subproblem[132]  person=132 obj=18.000000 62294(msec)
		subproblem[133]  person=133 obj=18.000000 61637(msec)
		subproblem[134]  person=134 obj=9.000000 23406(msec)
		subproblem[135]  person=135 obj=8.000000 62305(msec)
		subproblem[136]  person=136 obj=6.000000 41591(msec)
		subproblem[137]  person=137 obj=17.000000 21567(msec)
		subproblem[138]  person=138 obj=20.000000 54606(msec)
		subproblem[139]  person=139 obj=18.000000 60731(msec)
		subproblem[140]  person=140 obj=45.000000 3296(msec)
		subproblem[141]  person=141 obj=61.000000 64201(msec)
		subproblem[142]  person=142 obj=45.000000 44822(msec)
		subproblem[143]  person=143 obj=70.000000 89750(msec)
		subproblem[144]  person=144 obj=63.000000 1774(msec)
		subproblem[145]  person=145 obj=52.000000 1950(msec)
		subproblem[146]  person=146 obj=61.000000 2634(msec)
		subproblem[147]  person=147 obj=19.000000 46097(msec)
		subproblem[148]  person=148 obj=34.000000 57115(msec)
		subproblem[149]  person=149 obj=19.000000 66850(msec)
	*********UB=1.08457e+06(0)  6879.25(cpu sec)
		 scale=0.0001
		ov=3071.77 ext=119 depth=0 LB=1.79769e+308 UB=1.08457e+06 13553.670000(cpu sec) R_calc_time=6673.995054(cpu sec)
		ov=3060.75 ext=216 depth=0 LB=1.79769e+308 UB=1.08457e+06 23727.656000(cpu sec) R_calc_time=10173.523976(cpu sec)
		ov=3050.32 ext=332 depth=0 LB=1.79769e+308 UB=1.08457e+06 32283.480000(cpu sec) R_calc_time=8555.165989(cpu sec)
		ov=3043.18 ext=458 depth=0 LB=1.79769e+308 UB=1.08457e+06 39469.386000(cpu sec) R_calc_time=7185.224768(cpu sec)
		ov=3037.18 ext=597 depth=0 LB=1.79769e+308 UB=1.08457e+06 45862.829000(cpu sec) R_calc_time=6392.649065(cpu sec)
		ov=3033.24 ext=735 depth=0 LB=1.79769e+308 UB=1.08457e+06 53124.663000(cpu sec) R_calc_time=7260.942034(cpu sec)

2025年11月28日金曜日

COPT Revised OSI implementaion

 Next, we move on to the OSI implementation that enables COPT to be used not only as an LPSolver but also as a MIPSolver. To use COPT as a MIPSolver, it must be explicitly instructed, and for that purpose, we have added a mode "mip" since the previous update. One note is we restricted one thread during mip mode.


#ifdef COPT

#include "constraint_function.h"

#include 
static copt_env* env = NULL;

int OsiCOPTSolverInterface::error_handling(int errcode) const{
    char errmsg[COPT_BUFFSIZE];

    COPT_GetRetcodeMsg(errcode, errmsg, COPT_BUFFSIZE);
    string str = errmsg;
    str+="\n";
    //printf("ERROR %d: %s\n", errcode, errmsg);
    send_status_to_gui(str);
    return 0;
}
OsiCOPTSolverInterface::OsiCOPTSolverInterface()
{
    lpmethod = 2;//barrier
    int errcode = COPT_CreateEnv(&env);
    if (errcode) {
        error_handling(errcode);

    }
    else {
        errcode = COPT_CreateProb(env, &prob);
        

        if (errcode) {
            error_handling(errcode);
        }
        errcode = COPT_SetIntParam(prob, COPT_INTPARAM_LOGTOCONSOLE, 0);// COPT_INTPARAM_LOGGING, 0);
        if (errcode) {
            error_handling(errcode);
        }
        errcode=COPT_SetIntParam(prob, COPT_INTPARAM_CROSSOVER, 0);//crossover 禁止
        if (errcode) {
            error_handling(errcode);
        }
        errcode= COPT_SetIntParam(prob, COPT_INTPARAM_LPMETHOD, lpmethod);//barrier
        if (errcode) {
            error_handling(errcode);
        }

    }
    //*get_model() = env.CreateModel("lp_ex1");
    mip_obj = DBL_MAX;
    running = false;
    crossover = true;
    log_cnt = 0;
    mip = false;

}
OsiCOPTSolverInterface&
OsiCOPTSolverInterface::operator=(const OsiCOPTSolverInterface& rhs)
{

    int errcode=COPT_CreateCopy(rhs.get_prob(), &prob);
    objOffset = rhs.objOffset;

    running = rhs.running;
    crossover = rhs.crossover;

    col_lower_vec = rhs.col_lower_vec;
    col_upper_vec=rhs.col_upper_vec;

    row_lower_vec = rhs.row_lower_vec;
    row_upper_vec =rhs.row_upper_vec;
    col_cost_vec=rhs.col_cost_vec;
    lpmethod = rhs.lpmethod;
    mip = rhs.mip;

    errcode = COPT_SetIntParam(prob, COPT_INTPARAM_LOGTOCONSOLE, 0);// COPT_INTPARAM_LOGGING, 0);
    if (errcode) {
        error_handling(errcode);
    }
    

    return *this;


}

OsiCOPTSolverInterface::OsiCOPTSolverInterface(const OsiCOPTSolverInterface& original) {
    *this = original;
   
}
OsiCOPTSolverInterface::~OsiCOPTSolverInterface()
{
    // Delete problem and environment
    if (prob) {
        COPT_DeleteProb(&prob);
    }
    if (env) {
        COPT_DeleteEnv(&env);
    }
    mip_obj = DBL_MAX;
    running = false;
    crossover = true;
    log_cnt = 0;


}
void OsiCOPTSolverInterface::set_mip() {
    mip = true;
    int errcode = COPT_SetIntParam(prob, COPT_INTPARAM_THREADS,1);

    if (errcode) {
        error_handling(errcode);
    }
}

/// Get objective function sense (1 for min (default), -1 for max)
double OsiCOPTSolverInterface::getObjSense() const {
    return 1;//min
}
void OsiCOPTSolverInterface::set_time_limit(int sec) {
    //HighsStatus hs = this->highs->setOptionValue("time_limit", sec);
    int errcode = COPT_SetDblParam(prob, COPT_DBLPARAM_TIMELIMIT, sec);
    if (errcode) {
        error_handling(errcode);
    }
}
void OsiCOPTSolverInterface::setObjSense(double s) {

    ObjSense pass_sense = ObjSense::kMinimize;
    if (s == (double)ObjSense::kMaximize) pass_sense = ObjSense::kMaximize;
    if (pass_sense == ObjSense::kMinimize) {
        int errcode=COPT_SetObjSense(prob, COPT_MINIMIZE);
        if (errcode) {
            error_handling(errcode);
        }
    }
    else {
        int errcode = COPT_SetObjSense(prob, COPT_MAXIMIZE);
        if (errcode) {
            error_handling(errcode);
        }
    }
}
void OsiCOPTSolverInterface::writeMps(const char* filename, const char* extension, double objSense) const
{
    std::string fullname = std::string(filename) + "." + std::string(extension);
    int errcode=COPT_WriteMps( prob, fullname.c_str());
    if (errcode) {
        error_handling(errcode);
    }

}
void OsiCOPTSolverInterface::writeLp(const char* filename, const char* extension, double objSense) const
{
    std::string fullname = std::string(filename) + "." + std::string(extension);
    int errcode = COPT_WriteLp(prob, fullname.c_str());
    if (errcode) {
        error_handling(errcode);
    }

}
void OsiCOPTSolverInterface::setRowType(HighsInt index, char sense,
    double rightHandSide, double range) {
    // Assign arbitrary values so that compilation is clean
    double lo = 0;
    double hi = 1e200;
    this->convertSenseToBound(sense, rightHandSide, range, lo, hi);
    this->setRowBounds(index, lo, hi);
}
void OsiCOPTSolverInterface::setObjCoeff(HighsInt elementIndex,
    double elementValue) {

    int errcode=COPT_SetColObj(prob, 1,&elementIndex, &elementValue);
    if (errcode) {
        error_handling(errcode);
    }
    //this->highs->changeColCost(elementIndex, elementValue);
}
void OsiCOPTSolverInterface::setObjCoeffSet(const int* indexFirst, const int* indexLast, const double* coeffList) {

    OsiSolverInterface::setObjCoeffSet(indexFirst, indexLast - 1, coeffList);

}
void OsiCOPTSolverInterface::setColSetBounds(const HighsInt* indexFirst,
    const HighsInt* indexLast,
    const double* boundList) {

    OsiSolverInterface::setColSetBounds(indexFirst, indexLast - 1, boundList);
}
void OsiCOPTSolverInterface::setRowBounds(HighsInt elementIndex, double lower,
    double upper) {

    int errcode=COPT_SetRowLower(prob, 1, &elementIndex, &lower);
    if (errcode) {
        error_handling(errcode);
    }
    errcode = COPT_SetRowUpper(prob, 1, &elementIndex, &upper);
    if (errcode) {
        error_handling(errcode);
    }
    
}
void OsiCOPTSolverInterface::setColBounds(HighsInt elementIndex, double lower,
    double upper) {

    int errcode = COPT_SetColLower(prob, 1, &elementIndex, &lower);
    if (errcode) {
        error_handling(errcode);
    }
    errcode = COPT_SetColUpper(prob, 1, &elementIndex, &upper);
    if (errcode) {
        error_handling(errcode);
    }

}
void OsiCOPTSolverInterface::setRowSetBounds(const HighsInt* indexFirst,
    const HighsInt* indexLast,
    const double* boundList) {

    OsiSolverInterface::setRowSetBounds(indexFirst, indexLast - 1, boundList);
}


void OsiCOPTSolverInterface::branchAndBound() {
    
    // TODO
}
void OsiCOPTSolverInterface::set_lpmethod(int m,double d)
{
    lpmethod = m;
    if (lpmethod == 6) {
        send_status_to_gui("\t\tPDLP is selected.\n");
        int errcode = COPT_SetDblParam(prob, COPT_DBLPARAM_PDLPTOL, d);//Barrier
        if (errcode) {
            error_handling(errcode);
        }
        errcode = COPT_SetIntParam(prob, COPT_INTPARAM_GPUMODE, 1);//GPU mode=1
        if (errcode) {
            error_handling(errcode);
        }
    }
    else if (lpmethod == 2) {
        send_status_to_gui("\t\tBarrier is selected.\n");
        int errcode = COPT_SetIntParam(prob, COPT_INTPARAM_GPUMODE, 0);//GPU mode=0
        if (errcode) {
            error_handling(errcode);
        }
        //errcode = COPT_SetIntParam(prob, COPT_INTPARAM_GPUMODE, 1);//GPU mode=1
    }

}
void OsiCOPTSolverInterface::initialSolve() {
    if (mip) {
        int errcode = COPT_Solve(prob);
        if (errcode) {
            error_handling(errcode);
        }
        int lpstat = COPT_MIPSTATUS_UNSTARTED;
        errcode = COPT_GetIntAttr(prob, COPT_INTATTR_MIPSTATUS, &lpstat);
        if (errcode) {
            error_handling(errcode);
        }
        if (lpstat == COPT_MIPSTATUS_OPTIMAL) {
            sol.col_value.resize(getNumCols());
            int errcode = COPT_GetSolution(prob, &sol.col_value[0]);// 
            //int COPT_GetLpSolution(copt_prob * prob, double* value, double * slack(activity of constraints), double* rowDual(row price), double* redCost)
            if (errcode) {
                error_handling(errcode);
            }
        }
        return;
    }
    set_running(true);
    int errcode = COPT_SetIntParam(prob, COPT_INTPARAM_LPMETHOD, lpmethod);//Barrier
    if (errcode) {
        error_handling(errcode);
    }
Restart:
    errcode = COPT_SolveLp(prob);
    if (errcode) {
        error_handling(errcode);
    }
    int lpstat = COPT_LPSTATUS_UNSTARTED;
    double lpobjval;
    double* lpsol = NULL;
    int* colstat = NULL;
    errcode = COPT_GetIntAttr(prob, COPT_INTATTR_LPSTATUS, &lpstat);
    if (errcode) {
         error_handling(errcode);
    }
    if (lpstat == COPT_LPSTATUS_OPTIMAL) {
        sol.col_value.resize(getNumCols());
        sol.col_dual.resize(getNumCols());
        sol.row_dual.resize(getNumRows());
        sol.row_value.resize(getNumRows()); //getRowActivity
        int errcode = COPT_GetLpSolution(prob, &sol.col_value[0], &sol.row_value[0], &sol.row_dual[0], &sol.col_dual[0]);// 
        //int COPT_GetLpSolution(copt_prob * prob, double* value, double * slack(activity of constraints), double* rowDual(row price), double* redCost)
        if (errcode) {
            error_handling(errcode);
        }
    }
    else if (lpstat == COPT_LPSTATUS_NUMERICAL) {
        
        string str = "Numarical trouble encounterd.\n";
        send_status_to_gui(str);
        //writeMps("copt_trouble1");
        errcode = COPT_SetIntParam(prob, COPT_INTPARAM_LPMETHOD,0);//Simplex
        if (errcode) {
            error_handling(errcode);
        }
        goto Restart;

        
        errcode = COPT_SetIntParam(prob, COPT_INTPARAM_LPMETHOD, lpmethod);//Barrier
        if (errcode) {
            error_handling(errcode);
        }
    }
    else
        assert(0);
    errcode = COPT_SetIntParam(prob, COPT_INTPARAM_LPMETHOD, lpmethod);//Barrier
    if (errcode) {
        error_handling(errcode);
    }
    //this->status = this->highs->run();
    //const HighsSolution const_sol = this->highs->getSolution();
    //sol = const_sol;
    set_running(false);
}
void OsiCOPTSolverInterface::resolve() {

    initialSolve();
}
const double* OsiCOPTSolverInterface::getRowRange() const {

    if (this->rowRange != NULL) {
        delete[] this->rowRange;
    }

    HighsInt nrows = this->getNumRows();

    if (nrows == 0) {
        return this->rowRange;
    }

    this->rowRange = new double[nrows];

    for (HighsInt i = 0; i < nrows; i++) {
        // compute range for row i
        double lo = this->getRowLower()[i];
        double hi = this->getRowUpper()[i];
        double t1;
        char t2;
        this->convertBoundToSense(lo, hi, t2, t1, this->rowRange[i]);
    }

    return this->rowRange;
}
const double* OsiCOPTSolverInterface::getRightHandSide() const {
#
    if (this->rhs != NULL) {
        delete[] this->rhs;
    }

    HighsInt nrows = this->getNumRows();

    if (nrows == 0) {
        return this->rhs;
    }

    this->rhs = new double[nrows];

    for (HighsInt i = 0; i < nrows; i++) {
        // compute rhs for row i
        double lo = this->getRowLower()[i];
        double hi = this->getRowUpper()[i];
        double t1;
        char t2;
        this->convertBoundToSense(lo, hi, t2, this->rhs[i], t1);
    }

    return this->rhs;
}
const char* OsiCOPTSolverInterface::getRowSense() const {


    if (this->rowSense != NULL) {
        delete[] this->rowSense;
    }

    HighsInt nrows = this->getNumRows();

    if (nrows == 0) {
        return this->rowSense;
    }

    this->rowSense = new char[nrows];

    for (HighsInt i = 0; i < nrows; i++) {
        // compute sense for row i
        double lo = this->getRowLower()[i];
        double hi = this->getRowUpper()[i];
        double t1, t2;
        this->convertBoundToSense(lo, hi, this->rowSense[i], t1, t2);
    }

    return this->rowSense;
}
void OsiCOPTSolverInterface::make_crossover() {
    double d1 = getObjValue();
    {
        if (d1 == 0) {
            resolve();
        }

        int errcode = COPT_SetIntParam(prob, COPT_INTPARAM_CROSSOVER, 1);//crossover Enable
        if (errcode) {
            error_handling(errcode);
        }
        errcode = COPT_SetIntParam(prob, COPT_INTPARAM_LPMETHOD, 3);//crossover
        if (errcode) {
            error_handling(errcode);
        }
    }
    /*
    1: Dual simplex.
        2 : Barrier.
        3 : Crossover.
        4 : Concurrent(Use simplex and barrier simultaneously).
        5 : Choose between simplex and barrier automatically(Based on features such
            as sparsity and /or coefficients ranges).
        6 : First - order method(PDLP)
    */
    resolve();
    double d2 = getObjValue();
    if (d1) {
       // assert(fabs(d1 - d2) <= 1e-3);
    } {
        int errcode = COPT_SetIntParam(prob, COPT_INTPARAM_CROSSOVER, 0);//crossover 禁止
        if (errcode) {
            error_handling(errcode);
        }
        errcode = COPT_SetIntParam(prob, COPT_INTPARAM_LPMETHOD, lpmethod);//barrier
        if (errcode) {
            error_handling(errcode);
        }
    }
}

const CoinPackedMatrix* OsiCOPTSolverInterface::getMatrixByCol() const {

    if (this->matrixByCol != NULL) {
        delete this->matrixByCol;
    }

    HighsInt nrows = this->getNumRows();
    HighsInt ncols = this->getNumCols();
    HighsInt nelements = this->getNumElements();

    HighsInt* len = new int[ncols];
    HighsInt* start = new int[ncols + 1];
    HighsInt* index = new int[nelements];
    double* value = new double[nelements];

    int pReqSize = 0;
    /*
    int errcode = COPT_GetCols(prob, ncols,
        0,//const int* list,
        0,//int* colMatBeg,
        0,//int* colMatCnt,
        0,// int* colMatIdx,
        0,//double* colMatElem,
        nelements,
        &pReqSize);*/
    //vector list;
    //for (auto i = 0; i < ncols; i++) {
    //   list.push_back(i);
    //}
    int errcode = COPT_GetCols(prob, ncols,
        0,//const int* list,
        start,//int* colMatBeg,
        0,// lenは0で渡すFORMATにする,//int* colMatCnt,
        index,// int* colMatIdx,
        value,//double* colMatElem,
        nelements,
        &pReqSize);
    // copy data

    //assert(0);//TODO
    /*
    memcpy(start, &(this->highs->getLp().a_matrix_.start_[0]),
        (ncols + 1) * sizeof(HighsInt));
    memcpy(index, &(this->highs->getLp().a_matrix_.index_[0]),
        nelements * sizeof(HighsInt));
    memcpy(value, &(this->highs->getLp().a_matrix_.value_[0]),
        nelements * sizeof(double));
        */
    if (errcode) {
        error_handling(errcode);
    }
    for (HighsInt i = 0; i < ncols; i++) {//
        assert(start[i + 1] >= 0);
        len[i] = start[i + 1] - start[i];
    }
    

    this->matrixByCol = new CoinPackedMatrix();

    this->matrixByCol->assignMatrix(true, nrows, ncols, nelements, value, index,
        start, len);
    assert(this->matrixByCol->getNumCols() == ncols);
    assert(this->matrixByCol->getNumRows() == nrows);

    return this->matrixByCol;
}

const CoinPackedMatrix* OsiCOPTSolverInterface::getMatrixByRow() const {

    if (this->matrixByRow != NULL) {
        delete this->matrixByRow;
    }
    this->matrixByRow = new CoinPackedMatrix();
    this->matrixByRow->reverseOrderedCopyOf(*this->getMatrixByCol());

    return this->matrixByRow;
}
void OsiCOPTSolverInterface::assignProblem(CoinPackedMatrix*& matrix,
    double*& collb, double*& colub,
    double*& obj, double*& rowlb,
    double*& rowub) {

    loadProblem(*matrix, collb, colub, obj, rowlb, rowub);
    delete matrix;
    matrix = 0;
    delete[] collb;
    collb = 0;
    delete[] colub;
    colub = 0;
    delete[] obj;
    obj = 0;
    delete[] rowlb;
    rowlb = 0;
    delete[] rowub;
    rowub = 0;
}

void OsiCOPTSolverInterface::loadProblem(const CoinPackedMatrix& matrix,
    const double* collb,
    const double* colub,
    const double* obj, const char* rowsen,
    const double* rowrhs,
    const double* rowrng) {

    HighsInt numRow = matrix.getNumRows();

    double* rowlb = new double[numRow];
    double* rowub = new double[numRow];

    char* myrowsen = (char*)rowsen;
    bool rowsennull = false;
    double* myrowrhs = (double*)rowrhs;
    bool rowrhsnull = false;
    double* myrowrng = (double*)rowrng;
    bool rowrngnull = false;

    if (rowsen == NULL) {
        rowsennull = true;
        myrowsen = new char[numRow];
        for (HighsInt i = 0; i < numRow; i++) {
            myrowsen[i] = 'G';
        }
    }

    if (rowrhs == NULL) {
        rowsennull = true;
        myrowrhs = new double[numRow];
        for (HighsInt i = 0; i < numRow; i++) {
            myrowrhs[i] = 0.0;
        }
    }

    if (rowrng == NULL) {
        rowrngnull = true;
        myrowrng = new double[numRow];
        for (HighsInt i = 0; i < numRow; i++) {
            myrowrng[i] = 0.0;
        }
    }

    for (HighsInt i = 0; i < numRow; i++) {
        this->convertSenseToBound(myrowsen[i], myrowrhs[i], myrowrng[i], rowlb[i],
            rowub[i]);
    }

    this->loadProblem(matrix, collb, colub, obj, rowlb, rowub);

    delete[] rowlb;
    delete[] rowub;

    if (rowsennull) {
        delete[] myrowsen;
    }

    if (rowrhsnull) {
        delete[] myrowrhs;
    }

    if (rowrngnull) {
        delete[] myrowrng;
    }
}

void OsiCOPTSolverInterface::assignProblem(CoinPackedMatrix*& matrix,
    double*& collb, double*& colub,
    double*& obj, char*& rowsen,
    double*& rowrhs, double*& rowrng) {

    loadProblem(*matrix, collb, colub, obj, rowsen, rowrhs, rowrng);
    delete matrix;
    matrix = 0;
    delete[] collb;
    collb = 0;
    delete[] colub;
    colub = 0;
    delete[] obj;
    obj = 0;
    delete[] rowsen;
    rowsen = 0;
    delete[] rowrhs;
    rowrhs = 0;
    delete[] rowrng;
    rowrng = 0;
}




void OsiCOPTSolverInterface::loadProblem(
    const HighsInt numcols, const HighsInt numrows, const CoinBigIndex* start,
    const HighsInt* index, const double* value, const double* collb,
    const double* colub, const double* obj, const double* rowlb,
    const double* rowub) {

    double oldObjSense = this->getObjSense();
   
    //obj
    if (obj == 0) {
        col_cost_vec.clear();
    }
    else {
        col_cost_vec.resize(numcols);
        for (auto i = 0; i < numcols; i++) {
            col_cost_vec[i] = obj[i];
        }
    }
    //col
    col_lower_vec.resize(numcols);
    col_upper_vec.resize(numcols);
    for (auto i = 0; i < numcols; i++) {
        col_lower_vec[i] = collb[i];
        col_upper_vec[i] = colub[i];
    }
    //row
    row_lower_vec.resize(numrows);
    row_upper_vec.resize(numrows);
    for (auto i = 0; i < numrows; i++) {
        row_lower_vec[i] = rowlb[i];
        row_upper_vec[i] = rowub[i];
    }


    int errcode = COPT_LoadProb(prob,
        numcols,
        numrows,
        COPT_MINIMIZE,//The optimization sense, either COPT_MAXIMIZE or COPT_MINIMIZE.
        objOffset,//The constant part of the objective function.
        obj,//Objective coefficients of variables
        start,
        0,//const int* colMatCnt, //If colMatCnt is NULL, colMatBeg will need to have length of nCol+1,and the begin and end pointers to the i - th matrix column coefficients are defined using colMatBeg[i] and colMatBeg[i + 1]
        index,//const int* colMatIdx,
        value,//const double* colMatElem,
        0,//const char* colType, //If colType is NULL, all variables will be continuous.
        collb, // const double* colLower,
        colub,//const double* colUpper,
        0,//const char* rowSense,//If rowSense is NULL, then rowBound and rowUpper will be treated as lower  and upper bounds for constraints.T
        rowlb,//const double* rowBound,
        rowub,// double* rowUpper,
        0,//char const* const* colNames,//
        0);// char const* const* rowNames);Names of variables and constraints. Can be NULL
    
    if (errcode) {
        error_handling(errcode);
    }
    
}
void OsiCOPTSolverInterface::loadProblem(
    const HighsInt numcols, const HighsInt numrows, const CoinBigIndex* start,
    const HighsInt* index, const double* value, const double* collb,
    const double* colub, const double* obj, const char* rowsen,
    const double* rowrhs, const double* rowrng) {

    double* rowlb = new double[numrows];
    double* rowub = new double[numrows];

    for (HighsInt i = 0; i < numrows; i++) {
        this->convertSenseToBound(rowsen[i], rowrhs[i], rowrng[i], rowlb[i],
            rowub[i]);
    }

    this->loadProblem(numcols, numrows, start, index, value, collb, colub, obj,
        rowlb, rowub);

    delete[] rowlb;
    delete[] rowub;
}

void OsiCOPTSolverInterface::loadProblem(
    const CoinPackedMatrix& matrix, const double* collb, const double* colub,
    const double* obj, const double* rowlb, const double* rowub) {

    bool transpose = false;
    if (!matrix.isColOrdered()) {
        transpose = true;
        // ToDo: remove this hack
        //((CoinPackedMatrix *)&matrix)->transpose();
        ((CoinPackedMatrix*)&matrix)->reverseOrdering();
    }

    HighsInt numCol = matrix.getNumCols();
    HighsInt numRow = matrix.getNumRows();
    HighsInt num_nz = matrix.getNumElements();

    HighsInt* start = new int[numCol + 1];
    HighsInt* index = new int[num_nz];
    double* value = new double[num_nz];

    // get matrix data
    // const CoinBigIndex *vectorStarts = matrix.getVectorStarts();
    const HighsInt* vectorLengths = matrix.getVectorLengths();
    const double* elements = matrix.getElements();
    const HighsInt* indices = matrix.getIndices();

    // set matrix in HighsLp
    start[0] = 0;
    HighsInt nz = 0;
    for (HighsInt i = 0; i < numCol; i++) {
        start[i + 1] = start[i] + vectorLengths[i];
        CoinBigIndex first = matrix.getVectorFirst(i);
        for (HighsInt j = 0; j < vectorLengths[i]; j++) {
            index[nz] = indices[first + j];
            value[nz] = elements[first + j];
            nz++;
        }
    }
    assert(num_nz == nz);

    this->loadProblem(numCol, numRow, start, index, value, collb, colub, obj,
        rowlb, rowub);

    if (transpose) {
        //((CoinPackedMatrix)matrix).transpose();
        ((CoinPackedMatrix*)&matrix)->reverseOrdering();
    }

    delete[] start;
    delete[] index;
    delete[] value;
}
void OsiCOPTSolverInterface::addRow(const CoinPackedVectorBase& vec,
    const double rowlb, const double rowub) {


    int errcode = COPT_AddRow(prob,  vec.getNumElements(), vec.getIndices(), vec.getElements(), 0, rowlb, rowub, 0);
    if (errcode) {
        error_handling(errcode);
    }
    else {
        row_lower_vec.push_back(rowlb);
        row_upper_vec.push_back(rowub);
    }

 
}

void OsiCOPTSolverInterface::addRow(const CoinPackedVectorBase& vec,
    const char rowsen, const double rowrhs,
    const double rowrng) {

    // Assign arbitrary values so that compilation is clean
    double lb = 0;
    double ub = 1e200;
    this->convertSenseToBound(rowsen, rowrhs, rowrng, lb, ub);
    this->addRow(vec, lb, ub);
}

void OsiCOPTSolverInterface::addCol(const CoinPackedVectorBase& vec,
    const double collb, const double colub,
    const double obj) {

    int elements = vec.getNumElements();
    const int* index_ptr = vec.getIndices();
    const double* val_ptr = vec.getElements();
    
    
    int errcode = COPT_AddCol(prob, obj, elements, index_ptr,val_ptr, COPT_CONTINUOUS, collb, colub, 0);//Variable Type: 'C'
    if (errcode) {
        error_handling(errcode);
    }
    else {
        col_cost_vec.push_back(obj);
        col_lower_vec.push_back(collb);
        col_upper_vec.push_back(colub);
    }
    
}

void OsiCOPTSolverInterface::deleteCols(const HighsInt num,
    const HighsInt* colIndices) {
    int originalcols = getNumCols();
    assert(col_cost_vec.size() == getNumCols());

    int errcode = COPT_DelCols(prob, num, colIndices);
    if (errcode) {
        error_handling(errcode);
    }
    else {
        getColLowerCOPT();
        getColUpperCOPT();
        getObjCoefficientsCOPT();
        /*
        set del_set;
        for (int i = 0; i < num; i++) {
            int ind = colIndices[i];
            del_set.insert(ind);
        }
        vector new_col_cost_vec;
        vector new_col_upper, new_col_lower;
        for (int i = 0; i < originalcols; i++) {
            auto iter = del_set.find(i);
            if (iter == del_set.end()) {
                new_col_cost_vec.push_back(col_cost_vec[i]);
                new_col_lower.push_back(col_lower_vec[i]);
                new_col_upper.push_back(col_upper_vec[i]);
            }
        }
        assert(num == del_set.size());
        assert(col_cost_vec.size() == new_col_cost_vec.size() + del_set.size());
        
        col_cost_vec = new_col_cost_vec;
        col_lower_vec = new_col_lower;
        col_upper_vec = new_col_upper;
        assert(col_cost_vec.size() == col_lower_vec.size());
        assert(col_cost_vec.size() == col_upper_vec.size());
        assert(col_cost_vec.size() == getNumCols());
        */
    }

}

void OsiCOPTSolverInterface::deleteRows(const HighsInt num,
    const HighsInt* rowIndices) {
    int originalrows = getNumRows();
    int errcode = COPT_DelRows(prob, num, rowIndices);
    if (errcode) {
        error_handling(errcode);
    } else {
        getRowLowerCOPT();
        getRowUpperCOPT();
        /*
        set del_set;
        for (int i = 0; i < num; i++) {
            int ind = rowIndices[i];
            del_set.insert(ind);
        }
        
        vector new_row_upper, new_row_lower;
        for (int i = 0; i < originalrows; i++) {
            auto iter = del_set.find(i);
            if (iter == del_set.end()) {
                
                new_row_lower.push_back(row_lower_vec[i]);
                new_row_upper.push_back(row_upper_vec[i]);
            }
        }
        assert(num == del_set.size());
        assert(row_lower_vec.size() == new_row_lower.size() + del_set.size());

        
        row_lower_vec = new_row_lower;
        row_upper_vec = new_row_upper;
        assert(row_lower_vec.size() == getNumRows());
        */
    }
    
}

bool OsiCOPTSolverInterface::isContinuous(int colNumber) const
{
    char b;
    int errcode=COPT_GetColType(prob,1, &colNumber, &b);
    if (errcode) {
        error_handling(errcode);
    }


    return b== COPT_CONTINUOUS;
}
double OsiCOPTSolverInterface::getObjValue() const
{
    double lpobjval;
    if (mip) {
        int errcode = COPT_GetDblAttr(prob, COPT_DBLATTR_BESTOBJ, &lpobjval);
        if (errcode) {
            error_handling(errcode);
        }
        return lpobjval;
    }
    int errcode = COPT_GetDblAttr(prob, COPT_DBLATTR_LPOBJVAL, &lpobjval);
    if (errcode) {
        error_handling(errcode);
    }
    return lpobjval;
}
bool OsiCOPTSolverInterface::isAbandoned() const {

    return false;
}

bool OsiCOPTSolverInterface::isProvenOptimal() const {
    if (mip) {
        int lpstat = 0;
        int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_MIPSTATUS, &lpstat);
        if (errcode) {
            error_handling(errcode);
        }
        if (lpstat == COPT_MIPSTATUS_OPTIMAL) {
            return true;
        }
        assert(0);
        return false;
    }
    int* colstat = NULL;
    int lpstat = 0;
    int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_LPSTATUS, &lpstat);
    if (errcode) {
        error_handling(errcode);
    }
    return lpstat== COPT_LPSTATUS_OPTIMAL;
}

bool OsiCOPTSolverInterface::isProvenPrimalInfeasible() const {
    if (mip) {
        int lpstat = 0;
        int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_MIPSTATUS, &lpstat);
        if (errcode) {
            error_handling(errcode);
        }
        if (lpstat == COPT_MIPSTATUS_INFEASIBLE) {
            return true;
        }
        
        return false;
    }
    int* colstat = NULL;
    int lpstat = 0;
    int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_LPSTATUS, &lpstat);
    if (errcode) {
        error_handling(errcode);
    }
    return lpstat == COPT_LPSTATUS_INFEASIBLE;

    
}

bool OsiCOPTSolverInterface::isProvenDualInfeasible() const {
    assert(!mip);
    int* colstat = NULL;
    int lpstat = 0;
    int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_LPSTATUS, &lpstat);
    if (errcode) {
        error_handling(errcode);
    }
    return lpstat == COPT_LPSTATUS_INFEASIBLE;
    
}

bool OsiCOPTSolverInterface::isPrimalObjectiveLimitReached() const {

    return false;
}

bool OsiCOPTSolverInterface::isDualObjectiveLimitReached() const {
   

    return false;
}

bool OsiCOPTSolverInterface::isIterationLimitReached() const {
    int* colstat = NULL;
    int lpstat = 0;
    int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_LPSTATUS, &lpstat);
    if (errcode) {
        error_handling(errcode);
    }
    return lpstat == COPT_LPSTATUS_ITERLIMIT;
    
}
HighsInt OsiCOPTSolverInterface::getIterationCount() const {

    
    int iteration_count = 0;
    int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_BARRIERITER, &iteration_count);
    if (errcode) {
        error_handling(errcode);
    }

    
    return iteration_count;
}

HighsInt OsiCOPTSolverInterface::getNumCols() const {
    int elems = 0;
    int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_COLS, &elems);
    if (errcode) {
        error_handling(errcode);
    }
    return elems;
}

HighsInt OsiCOPTSolverInterface::getNumRows() const {
    int elems = 0;
    int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_ROWS, &elems);
    if (errcode) {
        error_handling(errcode);
    }
    return elems;
}

HighsInt OsiCOPTSolverInterface::getNumElements() const {
    int elems = 0;
    int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_ELEMS, &elems);
    if (errcode) {
        error_handling(errcode);
    }
    return elems;
}
void OsiCOPTSolverInterface::setRowLower(HighsInt elementIndex,
    double elementValue) {
    int errcode = COPT_SetRowLower(prob, 1, &elementIndex, &elementValue);
    if (errcode) {
        error_handling(errcode);
    }
    else {
        assert(elementIndex < row_lower_vec.size());
        row_lower_vec[elementIndex] = elementValue;
    }
   
}

void OsiCOPTSolverInterface::setRowUpper(HighsInt elementIndex,
    double elementValue) {

    int errcode = COPT_SetRowUpper(prob, 1, &elementIndex, &elementValue);
    if (errcode) {
        error_handling(errcode);
    }
    else {
        assert(elementIndex < row_upper_vec.size());
        row_upper_vec[elementIndex] = elementValue;
    }

}

void OsiCOPTSolverInterface::setColLower(HighsInt elementIndex,
    double elementValue) {
    int errcode = COPT_SetColLower(prob, 1, &elementIndex, &elementValue);
    if (errcode) {
        error_handling(errcode);
    }
    else {
        assert(elementIndex < col_lower_vec.size());
        col_lower_vec[elementIndex] = elementValue;
    }

}

void OsiCOPTSolverInterface::setColUpper(HighsInt elementIndex,
    double elementValue) {
    int errcode = COPT_SetColUpper(prob, 1, &elementIndex, &elementValue);
    if (errcode) {
        error_handling(errcode);
    }
    else {
        assert(elementIndex < col_upper_vec.size());
        col_upper_vec[elementIndex] = elementValue;
    }
}
void OsiCOPTSolverInterface::setContinuous(int index) {
    const char b = COPT_CONTINUOUS;
    int errcode=COPT_SetColType(prob, 1, &index, &b);
    if (errcode) {
        error_handling(errcode);
    }

}

/// Set the index-th variable to be an integer variable
void OsiCOPTSolverInterface::setInteger(int index) {
    const char b = COPT_INTEGER;
    int errcode = COPT_SetColType(prob, 1, &index, &b);
    if (errcode) {
        error_handling(errcode);
    }
}

const double* OsiCOPTSolverInterface::getColLower() const {
        
    assert(col_lower_vec.size());
    return &col_lower_vec[0];
}
const double* OsiCOPTSolverInterface::getColLowerCOPT() {

    col_lower_vec.resize(getNumCols(),0);

    int errcode = COPT_GetColInfo(prob, COPT_DBLINFO_LB, getNumCols(),
        0, &col_lower_vec[0]);
    if (errcode) {
        error_handling(errcode);
    }
    return &col_lower_vec[0];
}

const double* OsiCOPTSolverInterface::getColUpper() const {
    assert(col_upper_vec.size());
    return &col_upper_vec[0];
}

const double* OsiCOPTSolverInterface::getColUpperCOPT() {
    
    //assert(col_upper_vec.size() == getNumCols());
    col_upper_vec.resize(getNumCols(), 0);
    int errcode = COPT_GetColInfo(prob, COPT_DBLINFO_UB, getNumCols(),
        0, &col_upper_vec[0]);
    if (errcode) {
        error_handling(errcode);
    }
    return &col_upper_vec[0];
}
const double* OsiCOPTSolverInterface::getRowLowerCOPT() {
    //assert(row_lower_vec.size()==getNumRows());
    row_lower_vec.resize(getNumRows(), 0);
    int errcode = COPT_GetRowInfo(prob, COPT_DBLINFO_LB, getNumRows(),
        0, &row_lower_vec[0]);
    if (errcode) {
        error_handling(errcode);
    }


    return &row_lower_vec[0];
}
bool  OsiCOPTSolverInterface::compare( OsiCOPTSolverInterface& copt)
{
    if (getNumCols() != copt.getNumCols()) {
        return false;
    }
    if (getNumRows() != copt.getNumRows()) {
        return false;
    }

    getObjCoefficientsCOPT();
    copt.getObjCoefficientsCOPT();

    getColUpperCOPT();
    copt.getColUpperCOPT();
    
    getColLowerCOPT();
    copt.getColLowerCOPT();
    
    getRowUpperCOPT();
    copt.getRowUpperCOPT();
    
    getRowLowerCOPT();
    copt.getRowLowerCOPT();
    const CoinPackedMatrix* cp = getMatrixByCol();
    const CoinPackedMatrix* cp1 = copt.getMatrixByCol();
    bool identical2 = cp->isEquivalent(*cp1);
    if (!identical2) {
        assert(0);
        return false;
    }
    for (auto i = 0; i < getNumRows(); i++) {

        double obj1 = getRowLower()[i];
        double obj2 = copt.getRowLower()[i];
        if (fabs(obj1 - obj2) >= 1e-3) {
            assert(0);
            return false;
        }
        obj1 = getRowUpper()[i];
        obj2 = copt.getRowUpper()[i];
        if (fabs(obj1 - obj2) >= 1e-3) {
            assert(0);
            return false;
        }
    }


    for (auto i = 0; i < getNumCols(); i++) {
        double obj1 = getObjCoefficients()[i];
        double obj2 = copt.getObjCoefficients()[i];
        if (fabs(obj1 - obj2) >= 1e-3) {
            assert(0);
            return false;
        }
        obj1 = getColLower()[i];
        obj2 = copt.getColLower()[i];
        if (fabs(obj1 - obj2) >= 1e-3) {
            assert(0);
            return false;
        }
        obj1 = getColUpper()[i];
        obj2 = copt.getColUpper()[i];
        if (obj1 >= 1e30 && obj2 >= 1e30) {
            continue;
        }
        if (fabs(obj1 - obj2) >= 1e-3) {
            assert(0);
            return false;
        }
    }
    
    return true;

}
const double* OsiCOPTSolverInterface::getRowUpperCOPT()  {

    row_upper_vec.resize(getNumRows(), 0);
    //assert(row_upper_vec.size() == getNumRows());
    int errcode = COPT_GetRowInfo(prob, COPT_DBLINFO_UB, getNumRows(),
        0, &row_upper_vec[0]);
    if (errcode) {
        error_handling(errcode);
    }

    return &row_upper_vec[0];
}

const double* OsiCOPTSolverInterface::getObjCoefficientsCOPT() {
   // assert(col_cost_vec.size());
    //assert(col_cost_vec.size() == getNumCols());
    col_cost_vec.resize(getNumCols(), 0);
    int errcode = COPT_GetColInfo(prob, COPT_DBLINFO_OBJ, getNumCols(),
        0, &col_cost_vec[0]);
    if (errcode) {
        error_handling(errcode);
    }
    return &col_cost_vec[0];

    
}
const double* OsiCOPTSolverInterface::getRowLower() const {
    assert(row_lower_vec.size());
    return &row_lower_vec[0];
}

const double* OsiCOPTSolverInterface::getRowUpper() const {
    assert(row_upper_vec.size());
    return &row_upper_vec[0];
}

const double* OsiCOPTSolverInterface::getObjCoefficients() const {
    assert(col_cost_vec.size());
    return &col_cost_vec[0];
}

// TODO: review: 10^20?
double OsiCOPTSolverInterface::getInfinity() const {
    
    return kHighsInf;
}
const double* OsiCOPTSolverInterface::getColSolution() const {

    const double* ptr = &sol.col_value[0];
    
    return ptr;
}

const double* OsiCOPTSolverInterface::getRowPrice() const {

    return &sol.row_dual[0];
}

const double* OsiCOPTSolverInterface::getReducedCost() const {

    return &sol.col_dual[0];
}

const double* OsiCOPTSolverInterface::getRowActivity() const {


    return &sol.row_value[0];
}
int OsiCOPTSolverInterface::setBasisStatus(const int*colbasis, const int*rowbasis) {
    
    int errcode = COPT_SetBasis(prob, colbasis, rowbasis);
    if (errcode) {
        error_handling(errcode);
    }
    return 0;
}

void OsiCOPTSolverInterface::getBasisStatus(HighsInt* cstat,
    HighsInt* rstat) const {
    
    int errcode=COPT_GetBasis(prob,cstat,rstat);
    if (errcode) {
        error_handling(errcode);
    }


}

OsiSolverInterface* OsiCOPTSolverInterface::clone(bool copyData) const {

    if (!copyData) {
        OsiCOPTSolverInterface* cln = new OsiCOPTSolverInterface();
        return cln;

    }
    else {
        OsiCOPTSolverInterface* cln = new OsiCOPTSolverInterface(*this);
        

        return cln;
    }
}
void OsiCOPTSolverInterface::setRowPrice(const double* rowprice) {
    assert(getNumRows() == sol.row_dual.size());
    for (auto i = 0; i < getNumRows(); i++) {
        sol.row_value[i] = rowprice[i];
    }
    int errcode = COPT_SetLpSolution(prob, &sol.col_value[0], &sol.row_value[0], &sol.row_dual[0], &sol.col_dual[0]);// 
    //int COPT_GetLpSolution(copt_prob * prob, double* value, double * slack(activity of constraints), double* rowDual(row price), double* redCost)
    if (errcode) {
        error_handling(errcode);
    }
}

void OsiCOPTSolverInterface::setColSolution(const double* colsol) {
    assert(getNumCols() == sol.col_value.size());
    for (auto i = 0; i < getNumCols(); i++) {
        sol.col_value[i] = colsol[i];
    }
    int errcode = COPT_SetLpSolution(prob, &sol.col_value[0], &sol.row_value[0], &sol.row_dual[0], &sol.col_dual[0]);// 
    //int COPT_GetLpSolution(copt_prob * prob, double* value, double * slack(activity of constraints), double* rowDual(row price), double* redCost)
    if (errcode) {
        error_handling(errcode);
    }
   
}

std::vector OsiCOPTSolverInterface::getDualRays(HighsInt maxNumRays,
    bool fullRay) const {
    assert(0);
    return std::vector(0);
}

std::vector OsiCOPTSolverInterface::getPrimalRays(
    HighsInt maxNumRays) const {
    assert(0);
    return std::vector(0);
}

CoinWarmStart* OsiCOPTSolverInterface::getEmptyWarmStart() const {
    return 0;
}

CoinWarmStart* OsiCOPTSolverInterface::getWarmStart() const {
    return 0;
}

bool OsiCOPTSolverInterface::setWarmStart(const CoinWarmStart* warmstart) {
    //const HighsOptions& options = this->highs->getOptions();
    //highsLogDev(options.log_options, HighsLogType::kInfo,
    //    "Calling OsiHiGHSSolverInterface::setWarmStart()\n");
    return false;
}
#endif

2025年11月27日木曜日

Q.副師長の割り当てを記述してみましたが、正しいでしょうか?

 

Ans.

記述上は、問題ないと思います。多分、

今月平日=眼科外来月水金+眼科外来火木

という関係になっているので、現在の記述は、今月平日で一制約にまとめることが出来ますが、将来、優先度を各曜日で変えたい意図があって分けているのかな?と思いました。


2025年11月26日水曜日

Q.夜勤1年目の組合わせを禁止 を記述してみましたが正しいでしょうか?

 


Ans.

記述は問題ないと思います。ただし、ハード制約としていますが、ソフト制約にした方が良いと思います。

ハード制約は、絶対にKeepするという指示ですが、その使用は、最小限にしないとハードエラー(解がない事態)となってしまう確率が上がります。(夜勤者数は多分絶対ですが、新人はそれよりは優先度が低いと思われるので。その他記述も同様です。)

<ハード制約とするのは>
次の2種類に限ります。
1)絶対にキープしたい制約
2)入りの後の明け等、決まっている絶対的ルール

それ以外は、ソフト制約とします。ユーザ様のプロジェクトで、ハード制約だらけのプロジェクトを散見します。こういうプロジェクトでは、予定が入ったときに解が無い事態に陥り易くなります。絶対にキープしたい制約であるという気持ちは理解しますが、ハード制約は、同時に絶対に出来る制約である必要もあります。ハード列制約を見渡して、予定がフルに入っても全てのハード制約が全部出来るだろうか?と自問自動してみてください。少しでも出来ない可能性があるならば、ソフト制約とするべきです。

2025年11月25日火曜日

Q.ある人のペア禁止を複数の別な人にも適用したい

 現在は、以下で、シフト夜勤および長日勤のペア禁止を2組行っています。


その実装は、




つまり、ペア集合内で、そのシフトは最大1名、を記述すればOKです。
これを例えば、Staff6について、別なStaff8とのペアも禁止したいという場合は、項目を追加し、

のように記述します。

この実装は、青部追加とすればOKです。





2025年11月24日月曜日

Q.特定のスタッフだけ、夜勤回数を絶対にキープしたい

  Ans.

特定のスタッフだけは、夜勤回数を絶対にKeepし、その他スタッフはソフト制約としたい、そんな場合のシフト制約の仕方について述べます。

夜勤回数をハード制約として記述を追加します。ただし、ハード制約として作用するのは、指定スタッフだけとし、その他スタッフは、従来のソフト制約の適用とします。

下記のように夜勤回数設定プロパティシートに追加します。


グループ定義を追加します。


グループ集合を追加します
夜勤回数のグループ集合を変更、ハード制約を追加します。


集合を確認します。



以上で実装終了です。


以上で実装終了です。
以上です

2025年11月23日日曜日

COPT OSI implementation

Below is my COPT OSI implementation. Though I cannot include header files, I believe this approach will save time compared to building everything from scratch. The implementation is based on my Highs OSI work. For your reference only..
#ifdef COPT

#include "constraint_function.h"

#include 
static copt_env* env = NULL;

int OsiCOPTSolverInterface::error_handling(int errcode) const{
    char errmsg[COPT_BUFFSIZE];

    COPT_GetRetcodeMsg(errcode, errmsg, COPT_BUFFSIZE);
    string str = errmsg;
    str+="\n";
    //printf("ERROR %d: %s\n", errcode, errmsg);
    send_status_to_gui(str);
    return 0;
}
OsiCOPTSolverInterface::OsiCOPTSolverInterface()
{
    int errcode = COPT_CreateEnv(&env);
    if (errcode) {
        error_handling(errcode);

    }
    else {
        errcode = COPT_CreateProb(env, &prob);
        

        if (errcode) {
            error_handling(errcode);
        }
        errcode = COPT_SetIntParam(prob, COPT_INTPARAM_LOGTOCONSOLE, 0);// COPT_INTPARAM_LOGGING, 0);
        if (errcode) {
            error_handling(errcode);
        }
        errcode=COPT_SetIntParam(prob, COPT_INTPARAM_CROSSOVER, 0);//crossover 禁止
        if (errcode) {
            error_handling(errcode);
        }
        errcode= COPT_SetIntParam(prob, COPT_INTPARAM_LPMETHOD, 2);//barrier
        if (errcode) {
            error_handling(errcode);
        }

    }
    //*get_model() = env.CreateModel("lp_ex1");
    mip_obj = DBL_MAX;
    running = false;
    crossover = true;
    log_cnt = 0;


}
OsiCOPTSolverInterface&
OsiCOPTSolverInterface::operator=(const OsiCOPTSolverInterface& rhs)
{

    int errcode=COPT_CreateCopy(rhs.get_prob(), &prob);
    objOffset = rhs.objOffset;

    running = rhs.running;
    crossover = rhs.crossover;

    col_lower_vec = rhs.col_lower_vec;
    col_upper_vec=rhs.col_upper_vec;

    row_lower_vec = rhs.row_lower_vec;
    row_upper_vec =rhs.row_upper_vec;
    col_cost_vec=rhs.col_cost_vec;
    errcode = COPT_SetIntParam(prob, COPT_INTPARAM_LOGTOCONSOLE, 0);// COPT_INTPARAM_LOGGING, 0);
    if (errcode) {
        error_handling(errcode);
    }
    

    return *this;


}

OsiCOPTSolverInterface::OsiCOPTSolverInterface(const OsiCOPTSolverInterface& original) {
    *this = original;
   
}
OsiCOPTSolverInterface::~OsiCOPTSolverInterface()
{
    // Delete problem and environment
    if (prob) {
        COPT_DeleteProb(&prob);
    }
    if (env) {
        COPT_DeleteEnv(&env);
    }
    mip_obj = DBL_MAX;
    running = false;
    crossover = true;
    log_cnt = 0;


}








/// Get objective function sense (1 for min (default), -1 for max)
double OsiCOPTSolverInterface::getObjSense() const {
    return 1;//min
}

void OsiCOPTSolverInterface::setObjSense(double s) {

    ObjSense pass_sense = ObjSense::kMinimize;
    if (s == (double)ObjSense::kMaximize) pass_sense = ObjSense::kMaximize;
    if (pass_sense == ObjSense::kMinimize) {
        int errcode=COPT_SetObjSense(prob, COPT_MINIMIZE);
        if (errcode) {
            error_handling(errcode);
        }
    }
    else {
        int errcode = COPT_SetObjSense(prob, COPT_MAXIMIZE);
        if (errcode) {
            error_handling(errcode);
        }
    }
}
void OsiCOPTSolverInterface::writeMps(const char* filename, const char* extension, double objSense) const
{
    std::string fullname = std::string(filename) + "." + std::string(extension);
    int errcode=COPT_WriteMps( prob, fullname.c_str());
    if (errcode) {
        error_handling(errcode);
    }

}
void OsiCOPTSolverInterface::setRowType(HighsInt index, char sense,
    double rightHandSide, double range) {
    // Assign arbitrary values so that compilation is clean
    double lo = 0;
    double hi = 1e200;
    this->convertSenseToBound(sense, rightHandSide, range, lo, hi);
    this->setRowBounds(index, lo, hi);
}
void OsiCOPTSolverInterface::setObjCoeff(HighsInt elementIndex,
    double elementValue) {

    int errcode=COPT_SetColObj(prob, 1,&elementIndex, &elementValue);
    if (errcode) {
        error_handling(errcode);
    }
    //this->highs->changeColCost(elementIndex, elementValue);
}
void OsiCOPTSolverInterface::setObjCoeffSet(const int* indexFirst, const int* indexLast, const double* coeffList) {

    OsiSolverInterface::setObjCoeffSet(indexFirst, indexLast - 1, coeffList);

}
void OsiCOPTSolverInterface::setColSetBounds(const HighsInt* indexFirst,
    const HighsInt* indexLast,
    const double* boundList) {

    OsiSolverInterface::setColSetBounds(indexFirst, indexLast - 1, boundList);
}
void OsiCOPTSolverInterface::setRowBounds(HighsInt elementIndex, double lower,
    double upper) {

    int errcode=COPT_SetRowLower(prob, 1, &elementIndex, &lower);
    if (errcode) {
        error_handling(errcode);
    }
    errcode = COPT_SetRowUpper(prob, 1, &elementIndex, &upper);
    if (errcode) {
        error_handling(errcode);
    }
    
}
void OsiCOPTSolverInterface::setColBounds(HighsInt elementIndex, double lower,
    double upper) {

    int errcode = COPT_SetColLower(prob, 1, &elementIndex, &lower);
    if (errcode) {
        error_handling(errcode);
    }
    errcode = COPT_SetColUpper(prob, 1, &elementIndex, &upper);
    if (errcode) {
        error_handling(errcode);
    }

}
void OsiCOPTSolverInterface::setRowSetBounds(const HighsInt* indexFirst,
    const HighsInt* indexLast,
    const double* boundList) {

    OsiSolverInterface::setRowSetBounds(indexFirst, indexLast - 1, boundList);
}



void OsiCOPTSolverInterface::initialSolve() {
   
    set_running(true);
    Restart:
    int errcode = COPT_SolveLp(prob);
    if (errcode) {
        error_handling(errcode);
    }
    int lpstat = COPT_LPSTATUS_UNSTARTED;
    double lpobjval;
    double* lpsol = NULL;
    int* colstat = NULL;
    errcode = COPT_GetIntAttr(prob, COPT_INTATTR_LPSTATUS, &lpstat);
    if (errcode) {
         error_handling(errcode);
    }
    if (lpstat == COPT_LPSTATUS_OPTIMAL) {
        sol.col_value.resize(getNumCols());
        sol.col_dual.resize(getNumCols());
        sol.row_dual.resize(getNumRows());
        sol.row_value.resize(getNumRows()); //getRowActivity
        int errcode = COPT_GetLpSolution(prob, &sol.col_value[0], &sol.row_value[0], &sol.row_dual[0], &sol.col_dual[0]);// 
        //int COPT_GetLpSolution(copt_prob * prob, double* value, double * slack(activity of constraints), double* rowDual(row price), double* redCost)
        if (errcode) {
            error_handling(errcode);
        }
    }
    else if (lpstat == COPT_LPSTATUS_NUMERICAL) {//work around: if the barrier fails, we try with simplex.
        
        string str = "Numerical trouble encounterd.\n";
        send_status_to_gui(str);
        //writeMps("copt_trouble1");
        errcode = COPT_SetIntParam(prob, COPT_INTPARAM_LPMETHOD,0);//Simplex
        if (errcode) {
            error_handling(errcode);
        }
        goto Restart;

        
        errcode = COPT_SetIntParam(prob, COPT_INTPARAM_LPMETHOD, 2);//Barrier
        if (errcode) {
            error_handling(errcode);
        }
    }
    else
        assert(0);
    errcode = COPT_SetIntParam(prob, COPT_INTPARAM_LPMETHOD, 2);//Barrier
    if (errcode) {
        error_handling(errcode);
    }
    //this->status = this->highs->run();
    //const HighsSolution const_sol = this->highs->getSolution();
    //sol = const_sol;
    set_running(false);
}
void OsiCOPTSolverInterface::resolve() {

    initialSolve();
}
const double* OsiCOPTSolverInterface::getRowRange() const {

    if (this->rowRange != NULL) {
        delete[] this->rowRange;
    }

    HighsInt nrows = this->getNumRows();

    if (nrows == 0) {
        return this->rowRange;
    }

    this->rowRange = new double[nrows];

    for (HighsInt i = 0; i < nrows; i++) {
        // compute range for row i
        double lo = this->getRowLower()[i];
        double hi = this->getRowUpper()[i];
        double t1;
        char t2;
        this->convertBoundToSense(lo, hi, t2, t1, this->rowRange[i]);
    }

    return this->rowRange;
}
const double* OsiCOPTSolverInterface::getRightHandSide() const {
#
    if (this->rhs != NULL) {
        delete[] this->rhs;
    }

    HighsInt nrows = this->getNumRows();

    if (nrows == 0) {
        return this->rhs;
    }

    this->rhs = new double[nrows];

    for (HighsInt i = 0; i < nrows; i++) {
        // compute rhs for row i
        double lo = this->getRowLower()[i];
        double hi = this->getRowUpper()[i];
        double t1;
        char t2;
        this->convertBoundToSense(lo, hi, t2, this->rhs[i], t1);
    }

    return this->rhs;
}
const char* OsiCOPTSolverInterface::getRowSense() const {


    if (this->rowSense != NULL) {
        delete[] this->rowSense;
    }

    HighsInt nrows = this->getNumRows();

    if (nrows == 0) {
        return this->rowSense;
    }

    this->rowSense = new char[nrows];

    for (HighsInt i = 0; i < nrows; i++) {
        // compute sense for row i
        double lo = this->getRowLower()[i];
        double hi = this->getRowUpper()[i];
        double t1, t2;
        this->convertBoundToSense(lo, hi, this->rowSense[i], t1, t2);
    }

    return this->rowSense;
}
}

const CoinPackedMatrix* OsiCOPTSolverInterface::getMatrixByCol() const {

    if (this->matrixByCol != NULL) {
        delete this->matrixByCol;
    }

    HighsInt nrows = this->getNumRows();
    HighsInt ncols = this->getNumCols();
    HighsInt nelements = this->getNumElements();

    HighsInt* len = new int[ncols];
    HighsInt* start = new int[ncols + 1];
    HighsInt* index = new int[nelements];
    double* value = new double[nelements];

    int pReqSize = 0;
    /*
    int errcode = COPT_GetCols(prob, ncols,
        0,//const int* list,
        0,//int* colMatBeg,
        0,//int* colMatCnt,
        0,// int* colMatIdx,
        0,//double* colMatElem,
        nelements,
        &pReqSize);*/
    //vector list;
    //for (auto i = 0; i < ncols; i++) {
    //   list.push_back(i);
    //}
    int errcode = COPT_GetCols(prob, ncols,
        0,//const int* list,
        start,//int* colMatBeg,
        0,// lenは0で渡すFORMATにする,//int* colMatCnt,
        index,// int* colMatIdx,
        value,//double* colMatElem,
        nelements,
        &pReqSize);
    // copy data

    //assert(0);//TODO
    /*
    memcpy(start, &(this->highs->getLp().a_matrix_.start_[0]),
        (ncols + 1) * sizeof(HighsInt));
    memcpy(index, &(this->highs->getLp().a_matrix_.index_[0]),
        nelements * sizeof(HighsInt));
    memcpy(value, &(this->highs->getLp().a_matrix_.value_[0]),
        nelements * sizeof(double));
        */
    if (errcode) {
        error_handling(errcode);
    }
    for (HighsInt i = 0; i < ncols; i++) {//
        assert(start[i + 1] >= 0);
        len[i] = start[i + 1] - start[i];
    }
    

    this->matrixByCol = new CoinPackedMatrix();

    this->matrixByCol->assignMatrix(true, nrows, ncols, nelements, value, index,
        start, len);
    assert(this->matrixByCol->getNumCols() == ncols);
    assert(this->matrixByCol->getNumRows() == nrows);

    return this->matrixByCol;
}

const CoinPackedMatrix* OsiCOPTSolverInterface::getMatrixByRow() const {

    if (this->matrixByRow != NULL) {
        delete this->matrixByRow;
    }
    this->matrixByRow = new CoinPackedMatrix();
    this->matrixByRow->reverseOrderedCopyOf(*this->getMatrixByCol());

    return this->matrixByRow;
}
void OsiCOPTSolverInterface::assignProblem(CoinPackedMatrix*& matrix,
    double*& collb, double*& colub,
    double*& obj, double*& rowlb,
    double*& rowub) {

    loadProblem(*matrix, collb, colub, obj, rowlb, rowub);
    delete matrix;
    matrix = 0;
    delete[] collb;
    collb = 0;
    delete[] colub;
    colub = 0;
    delete[] obj;
    obj = 0;
    delete[] rowlb;
    rowlb = 0;
    delete[] rowub;
    rowub = 0;
}

void OsiCOPTSolverInterface::loadProblem(const CoinPackedMatrix& matrix,
    const double* collb,
    const double* colub,
    const double* obj, const char* rowsen,
    const double* rowrhs,
    const double* rowrng) {

    HighsInt numRow = matrix.getNumRows();

    double* rowlb = new double[numRow];
    double* rowub = new double[numRow];

    char* myrowsen = (char*)rowsen;
    bool rowsennull = false;
    double* myrowrhs = (double*)rowrhs;
    bool rowrhsnull = false;
    double* myrowrng = (double*)rowrng;
    bool rowrngnull = false;

    if (rowsen == NULL) {
        rowsennull = true;
        myrowsen = new char[numRow];
        for (HighsInt i = 0; i < numRow; i++) {
            myrowsen[i] = 'G';
        }
    }

    if (rowrhs == NULL) {
        rowsennull = true;
        myrowrhs = new double[numRow];
        for (HighsInt i = 0; i < numRow; i++) {
            myrowrhs[i] = 0.0;
        }
    }

    if (rowrng == NULL) {
        rowrngnull = true;
        myrowrng = new double[numRow];
        for (HighsInt i = 0; i < numRow; i++) {
            myrowrng[i] = 0.0;
        }
    }

    for (HighsInt i = 0; i < numRow; i++) {
        this->convertSenseToBound(myrowsen[i], myrowrhs[i], myrowrng[i], rowlb[i],
            rowub[i]);
    }

    this->loadProblem(matrix, collb, colub, obj, rowlb, rowub);

    delete[] rowlb;
    delete[] rowub;

    if (rowsennull) {
        delete[] myrowsen;
    }

    if (rowrhsnull) {
        delete[] myrowrhs;
    }

    if (rowrngnull) {
        delete[] myrowrng;
    }
}

void OsiCOPTSolverInterface::assignProblem(CoinPackedMatrix*& matrix,
    double*& collb, double*& colub,
    double*& obj, char*& rowsen,
    double*& rowrhs, double*& rowrng) {

    loadProblem(*matrix, collb, colub, obj, rowsen, rowrhs, rowrng);
    delete matrix;
    matrix = 0;
    delete[] collb;
    collb = 0;
    delete[] colub;
    colub = 0;
    delete[] obj;
    obj = 0;
    delete[] rowsen;
    rowsen = 0;
    delete[] rowrhs;
    rowrhs = 0;
    delete[] rowrng;
    rowrng = 0;
}




void OsiCOPTSolverInterface::loadProblem(
    const HighsInt numcols, const HighsInt numrows, const CoinBigIndex* start,
    const HighsInt* index, const double* value, const double* collb,
    const double* colub, const double* obj, const double* rowlb,
    const double* rowub) {

    double oldObjSense = this->getObjSense();
   
    //obj
    if (obj == 0) {
        col_cost_vec.clear();
    }
    else {
        col_cost_vec.resize(numcols);
        for (auto i = 0; i < numcols; i++) {
            col_cost_vec[i] = obj[i];
        }
    }
    //col
    col_lower_vec.resize(numcols);
    col_upper_vec.resize(numcols);
    for (auto i = 0; i < numcols; i++) {
        col_lower_vec[i] = collb[i];
        col_upper_vec[i] = colub[i];
    }
    //row
    row_lower_vec.resize(numrows);
    row_upper_vec.resize(numrows);
    for (auto i = 0; i < numrows; i++) {
        row_lower_vec[i] = rowlb[i];
        row_upper_vec[i] = rowub[i];
    }


    int errcode = COPT_LoadProb(prob,
        numcols,
        numrows,
        COPT_MINIMIZE,//The optimization sense, either COPT_MAXIMIZE or COPT_MINIMIZE.
        objOffset,//The constant part of the objective function.
        obj,//Objective coefficients of variables
        start,
        0,//const int* colMatCnt, //If colMatCnt is NULL, colMatBeg will need to have length of nCol+1,and the begin and end pointers to the i - th matrix column coefficients are defined using colMatBeg[i] and colMatBeg[i + 1]
        index,//const int* colMatIdx,
        value,//const double* colMatElem,
        0,//const char* colType, //If colType is NULL, all variables will be continuous.
        collb, // const double* colLower,
        colub,//const double* colUpper,
        0,//const char* rowSense,//If rowSense is NULL, then rowBound and rowUpper will be treated as lower  and upper bounds for constraints.T
        rowlb,//const double* rowBound,
        rowub,// double* rowUpper,
        0,//char const* const* colNames,//
        0);// char const* const* rowNames);Names of variables and constraints. Can be NULL
    
    if (errcode) {
        error_handling(errcode);
    }
    
}
void OsiCOPTSolverInterface::loadProblem(
    const HighsInt numcols, const HighsInt numrows, const CoinBigIndex* start,
    const HighsInt* index, const double* value, const double* collb,
    const double* colub, const double* obj, const char* rowsen,
    const double* rowrhs, const double* rowrng) {

    double* rowlb = new double[numrows];
    double* rowub = new double[numrows];

    for (HighsInt i = 0; i < numrows; i++) {
        this->convertSenseToBound(rowsen[i], rowrhs[i], rowrng[i], rowlb[i],
            rowub[i]);
    }

    this->loadProblem(numcols, numrows, start, index, value, collb, colub, obj,
        rowlb, rowub);

    delete[] rowlb;
    delete[] rowub;
}

void OsiCOPTSolverInterface::loadProblem(
    const CoinPackedMatrix& matrix, const double* collb, const double* colub,
    const double* obj, const double* rowlb, const double* rowub) {

    bool transpose = false;
    if (!matrix.isColOrdered()) {
        transpose = true;
        // ToDo: remove this hack
        //((CoinPackedMatrix *)&matrix)->transpose();
        ((CoinPackedMatrix*)&matrix)->reverseOrdering();
    }

    HighsInt numCol = matrix.getNumCols();
    HighsInt numRow = matrix.getNumRows();
    HighsInt num_nz = matrix.getNumElements();

    HighsInt* start = new int[numCol + 1];
    HighsInt* index = new int[num_nz];
    double* value = new double[num_nz];

    // get matrix data
    // const CoinBigIndex *vectorStarts = matrix.getVectorStarts();
    const HighsInt* vectorLengths = matrix.getVectorLengths();
    const double* elements = matrix.getElements();
    const HighsInt* indices = matrix.getIndices();

    // set matrix in HighsLp
    start[0] = 0;
    HighsInt nz = 0;
    for (HighsInt i = 0; i < numCol; i++) {
        start[i + 1] = start[i] + vectorLengths[i];
        CoinBigIndex first = matrix.getVectorFirst(i);
        for (HighsInt j = 0; j < vectorLengths[i]; j++) {
            index[nz] = indices[first + j];
            value[nz] = elements[first + j];
            nz++;
        }
    }
    assert(num_nz == nz);

    this->loadProblem(numCol, numRow, start, index, value, collb, colub, obj,
        rowlb, rowub);

    if (transpose) {
        //((CoinPackedMatrix)matrix).transpose();
        ((CoinPackedMatrix*)&matrix)->reverseOrdering();
    }

    delete[] start;
    delete[] index;
    delete[] value;
}
void OsiCOPTSolverInterface::addRow(const CoinPackedVectorBase& vec,
    const double rowlb, const double rowub) {


    int errcode = COPT_AddRow(prob,  vec.getNumElements(), vec.getIndices(), vec.getElements(), 0, rowlb, rowub, 0);
    if (errcode) {
        error_handling(errcode);
    }
    else {
        row_lower_vec.push_back(rowlb);
        row_upper_vec.push_back(rowub);
    }

 
}

void OsiCOPTSolverInterface::addRow(const CoinPackedVectorBase& vec,
    const char rowsen, const double rowrhs,
    const double rowrng) {

    // Assign arbitrary values so that compilation is clean
    double lb = 0;
    double ub = 1e200;
    this->convertSenseToBound(rowsen, rowrhs, rowrng, lb, ub);
    this->addRow(vec, lb, ub);
}

void OsiCOPTSolverInterface::addCol(const CoinPackedVectorBase& vec,
    const double collb, const double colub,
    const double obj) {

    int elements = vec.getNumElements();
    const int* index_ptr = vec.getIndices();
    const double* val_ptr = vec.getElements();
    
    
    int errcode = COPT_AddCol(prob, obj, elements, index_ptr,val_ptr, COPT_CONTINUOUS, collb, colub, 0);//Variable Type: 'C'
    if (errcode) {
        error_handling(errcode);
    }
    else {
        col_cost_vec.push_back(obj);
        col_lower_vec.push_back(collb);
        col_upper_vec.push_back(colub);
    }
    
}

void OsiCOPTSolverInterface::deleteCols(const HighsInt num,
    const HighsInt* colIndices) {
    int originalcols = getNumCols();
    assert(col_cost_vec.size() == getNumCols());

    int errcode = COPT_DelCols(prob, num, colIndices);
    if (errcode) {
        error_handling(errcode);
    }
    else {
        getColLowerCOPT();
        getColUpperCOPT();
        getObjCoefficientsCOPT();
        /*
        set del_set;
        for (int i = 0; i < num; i++) {
            int ind = colIndices[i];
            del_set.insert(ind);
        }
        vector new_col_cost_vec;
        vector new_col_upper, new_col_lower;
        for (int i = 0; i < originalcols; i++) {
            auto iter = del_set.find(i);
            if (iter == del_set.end()) {
                new_col_cost_vec.push_back(col_cost_vec[i]);
                new_col_lower.push_back(col_lower_vec[i]);
                new_col_upper.push_back(col_upper_vec[i]);
            }
        }
        assert(num == del_set.size());
        assert(col_cost_vec.size() == new_col_cost_vec.size() + del_set.size());
        
        col_cost_vec = new_col_cost_vec;
        col_lower_vec = new_col_lower;
        col_upper_vec = new_col_upper;
        assert(col_cost_vec.size() == col_lower_vec.size());
        assert(col_cost_vec.size() == col_upper_vec.size());
        assert(col_cost_vec.size() == getNumCols());
        */
    }

}

void OsiCOPTSolverInterface::deleteRows(const HighsInt num,
    const HighsInt* rowIndices) {
    int originalrows = getNumRows();
    int errcode = COPT_DelRows(prob, num, rowIndices);
    if (errcode) {
        error_handling(errcode);
    } else {
        getRowLowerCOPT();
        getRowUpperCOPT();
        /*
        set del_set;
        for (int i = 0; i < num; i++) {
            int ind = rowIndices[i];
            del_set.insert(ind);
        }
        
        vector new_row_upper, new_row_lower;
        for (int i = 0; i < originalrows; i++) {
            auto iter = del_set.find(i);
            if (iter == del_set.end()) {
                
                new_row_lower.push_back(row_lower_vec[i]);
                new_row_upper.push_back(row_upper_vec[i]);
            }
        }
        assert(num == del_set.size());
        assert(row_lower_vec.size() == new_row_lower.size() + del_set.size());

        
        row_lower_vec = new_row_lower;
        row_upper_vec = new_row_upper;
        assert(row_lower_vec.size() == getNumRows());
        */
    }
    
}

bool OsiCOPTSolverInterface::isContinuous(int colNumber) const
{
    char b;
    int errcode=COPT_GetColType(prob,1, &colNumber, &b);
    if (errcode) {
        error_handling(errcode);
    }


    return b== COPT_CONTINUOUS;
}
double OsiCOPTSolverInterface::getObjValue() const
{
    double lpobjval;
    int errcode = COPT_GetDblAttr(prob, COPT_DBLATTR_LPOBJVAL, &lpobjval);
    if (errcode) {
        error_handling(errcode);
    }
    return lpobjval;
}
bool OsiCOPTSolverInterface::isAbandoned() const {

    return false;
}

bool OsiCOPTSolverInterface::isProvenOptimal() const {
    int* colstat = NULL;
    int lpstat = 0;
    int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_LPSTATUS, &lpstat);
    if (errcode) {
        error_handling(errcode);
    }
    return lpstat== COPT_LPSTATUS_OPTIMAL;
}

bool OsiCOPTSolverInterface::isProvenPrimalInfeasible() const {
    int* colstat = NULL;
    int lpstat = 0;
    int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_LPSTATUS, &lpstat);
    if (errcode) {
        error_handling(errcode);
    }
    return lpstat == COPT_LPSTATUS_INFEASIBLE;

    
}

bool OsiCOPTSolverInterface::isProvenDualInfeasible() const {
    int* colstat = NULL;
    int lpstat = 0;
    int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_LPSTATUS, &lpstat);
    if (errcode) {
        error_handling(errcode);
    }
    return lpstat == COPT_LPSTATUS_INFEASIBLE;
    
}

bool OsiCOPTSolverInterface::isPrimalObjectiveLimitReached() const {

    return false;
}

bool OsiCOPTSolverInterface::isDualObjectiveLimitReached() const {
   

    return false;
}

bool OsiCOPTSolverInterface::isIterationLimitReached() const {
    int* colstat = NULL;
    int lpstat = 0;
    int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_LPSTATUS, &lpstat);
    if (errcode) {
        error_handling(errcode);
    }
    return lpstat == COPT_LPSTATUS_ITERLIMIT;
    
}
HighsInt OsiCOPTSolverInterface::getIterationCount() const {

    
    int iteration_count = 0;
    int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_BARRIERITER, &iteration_count);
    if (errcode) {
        error_handling(errcode);
    }

    
    return iteration_count;
}

HighsInt OsiCOPTSolverInterface::getNumCols() const {
    int elems = 0;
    int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_COLS, &elems);
    if (errcode) {
        error_handling(errcode);
    }
    return elems;
}

HighsInt OsiCOPTSolverInterface::getNumRows() const {
    int elems = 0;
    int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_ROWS, &elems);
    if (errcode) {
        error_handling(errcode);
    }
    return elems;
}

HighsInt OsiCOPTSolverInterface::getNumElements() const {
    int elems = 0;
    int errcode = COPT_GetIntAttr(prob, COPT_INTATTR_ELEMS, &elems);
    if (errcode) {
        error_handling(errcode);
    }
    return elems;
}
void OsiCOPTSolverInterface::setRowLower(HighsInt elementIndex,
    double elementValue) {
    int errcode = COPT_SetRowLower(prob, 1, &elementIndex, &elementValue);
    if (errcode) {
        error_handling(errcode);
    }
    else {
        assert(elementIndex < row_lower_vec.size());
        row_lower_vec[elementIndex] = elementValue;
    }
   
}

void OsiCOPTSolverInterface::setRowUpper(HighsInt elementIndex,
    double elementValue) {

    int errcode = COPT_SetRowUpper(prob, 1, &elementIndex, &elementValue);
    if (errcode) {
        error_handling(errcode);
    }
    else {
        assert(elementIndex < row_upper_vec.size());
        row_upper_vec[elementIndex] = elementValue;
    }

}

void OsiCOPTSolverInterface::setColLower(HighsInt elementIndex,
    double elementValue) {
    int errcode = COPT_SetColLower(prob, 1, &elementIndex, &elementValue);
    if (errcode) {
        error_handling(errcode);
    }
    else {
        assert(elementIndex < col_lower_vec.size());
        col_lower_vec[elementIndex] = elementValue;
    }

}

void OsiCOPTSolverInterface::setColUpper(HighsInt elementIndex,
    double elementValue) {
    int errcode = COPT_SetColUpper(prob, 1, &elementIndex, &elementValue);
    if (errcode) {
        error_handling(errcode);
    }
    else {
        assert(elementIndex < col_upper_vec.size());
        col_upper_vec[elementIndex] = elementValue;
    }
}
void OsiCOPTSolverInterface::setContinuous(int index) {
    const char b = COPT_CONTINUOUS;
    int errcode=COPT_SetColType(prob, 1, &index, &b);
    if (errcode) {
        error_handling(errcode);
    }

}

/// Set the index-th variable to be an integer variable
void OsiCOPTSolverInterface::setInteger(int index) {
    const char b = COPT_INTEGER;
    int errcode = COPT_SetColType(prob, 1, &index, &b);
    if (errcode) {
        error_handling(errcode);
    }
}

const double* OsiCOPTSolverInterface::getColLower() const {
        
    assert(col_lower_vec.size());
    return &col_lower_vec[0];
}
const double* OsiCOPTSolverInterface::getColLowerCOPT() {

    col_lower_vec.resize(getNumCols(),0);

    int errcode = COPT_GetColInfo(prob, COPT_DBLINFO_LB, getNumCols(),
        0, &col_lower_vec[0]);
    if (errcode) {
        error_handling(errcode);
    }
    return &col_lower_vec[0];
}

const double* OsiCOPTSolverInterface::getColUpper() const {
    assert(col_upper_vec.size());
    return &col_upper_vec[0];
}

const double* OsiCOPTSolverInterface::getColUpperCOPT() {
    
    //assert(col_upper_vec.size() == getNumCols());
    col_upper_vec.resize(getNumCols(), 0);
    int errcode = COPT_GetColInfo(prob, COPT_DBLINFO_UB, getNumCols(),
        0, &col_upper_vec[0]);
    if (errcode) {
        error_handling(errcode);
    }
    return &col_upper_vec[0];
}
const double* OsiCOPTSolverInterface::getRowLowerCOPT() {
    //assert(row_lower_vec.size()==getNumRows());
    row_lower_vec.resize(getNumRows(), 0);
    int errcode = COPT_GetRowInfo(prob, COPT_DBLINFO_LB, getNumRows(),
        0, &row_lower_vec[0]);
    if (errcode) {
        error_handling(errcode);
    }


    return &row_lower_vec[0];
}
bool  OsiCOPTSolverInterface::compare( OsiCOPTSolverInterface& copt)
{
    if (getNumCols() != copt.getNumCols()) {
        return false;
    }
    if (getNumRows() != copt.getNumRows()) {
        return false;
    }

    getObjCoefficientsCOPT();
    copt.getObjCoefficientsCOPT();

    getColUpperCOPT();
    copt.getColUpperCOPT();
    
    getColLowerCOPT();
    copt.getColLowerCOPT();
    
    getRowUpperCOPT();
    copt.getRowUpperCOPT();
    
    getRowLowerCOPT();
    copt.getRowLowerCOPT();
    const CoinPackedMatrix* cp = getMatrixByCol();
    const CoinPackedMatrix* cp1 = copt.getMatrixByCol();
    bool identical2 = cp->isEquivalent(*cp1);
    if (!identical2) {
        return false;
    }
    for (auto i = 0; i < getNumRows(); i++) {

        double obj1 = getRowLower()[i];
        double obj2 = copt.getRowLower()[i];
        if (fabs(obj1 - obj2) >= 1e-3) {
            return false;
        }
        obj1 = getRowUpper()[i];
        obj2 = copt.getRowUpper()[i];
        if (fabs(obj1 - obj2) >= 1e-3) {
            return false;
        }
    }


    for (auto i = 0; i < getNumCols(); i++) {
        double obj1 = getObjCoefficients()[i];
        double obj2 = copt.getObjCoefficients()[i];
        if (fabs(obj1 - obj2) >= 1e-3) {
            return false;
        }
        obj1 = getColLower()[i];
        obj2 = copt.getColLower()[i];
        if (fabs(obj1 - obj2) >= 1e-3) {
            return false;
        }
        obj1 = getColUpper()[i];
        obj2 = copt.getColUpper()[i];
        if (obj1 >= 1e30 && obj2 >= 1e30) {
            continue;
        }
        if (fabs(obj1 - obj2) >= 1e-3) {
            return false;
        }
    }
    
    return true;

}
const double* OsiCOPTSolverInterface::getRowUpperCOPT()  {

    row_upper_vec.resize(getNumRows(), 0);
    //assert(row_upper_vec.size() == getNumRows());
    int errcode = COPT_GetRowInfo(prob, COPT_DBLINFO_UB, getNumRows(),
        0, &row_upper_vec[0]);
    if (errcode) {
        error_handling(errcode);
    }

    return &row_upper_vec[0];
}

const double* OsiCOPTSolverInterface::getObjCoefficientsCOPT() {
   // assert(col_cost_vec.size());
    //assert(col_cost_vec.size() == getNumCols());
    col_cost_vec.resize(getNumCols(), 0);
    int errcode = COPT_GetColInfo(prob, COPT_DBLINFO_OBJ, getNumCols(),
        0, &col_cost_vec[0]);
    if (errcode) {
        error_handling(errcode);
    }
    return &col_cost_vec[0];

    
}
const double* OsiCOPTSolverInterface::getRowLower() const {
    assert(row_lower_vec.size());
    return &row_lower_vec[0];
}

const double* OsiCOPTSolverInterface::getRowUpper() const {
    assert(row_upper_vec.size());
    return &row_upper_vec[0];
}

const double* OsiCOPTSolverInterface::getObjCoefficients() const {
    assert(col_cost_vec.size());
    return &col_cost_vec[0];
}

// TODO: review: 10^20?
double OsiCOPTSolverInterface::getInfinity() const {
    
    return kHighsInf;
}
const double* OsiCOPTSolverInterface::getColSolution() const {

    const double* ptr = &sol.col_value[0];
    
    return ptr;
}

const double* OsiCOPTSolverInterface::getRowPrice() const {

    return &sol.row_dual[0];
}

const double* OsiCOPTSolverInterface::getReducedCost() const {

    return &sol.col_dual[0];
}

const double* OsiCOPTSolverInterface::getRowActivity() const {


    return &sol.row_value[0];
}
int OsiCOPTSolverInterface::setBasisStatus(const int*colbasis, const int*rowbasis) {
    
    int errcode = COPT_SetBasis(prob, colbasis, rowbasis);
    if (errcode) {
        error_handling(errcode);
    }
    return 0;
}

void OsiCOPTSolverInterface::getBasisStatus(HighsInt* cstat,
    HighsInt* rstat) const {
    
    int errcode=COPT_GetBasis(prob,cstat,rstat);
    if (errcode) {
        error_handling(errcode);
    }


}

OsiSolverInterface* OsiCOPTSolverInterface::clone(bool copyData) const {

    if (!copyData) {
        OsiCOPTSolverInterface* cln = new OsiCOPTSolverInterface();
        return cln;

    }
    else {
        OsiCOPTSolverInterface* cln = new OsiCOPTSolverInterface(*this);
        

        return cln;
    }
}
void OsiCOPTSolverInterface::setRowPrice(const double* rowprice) {
    assert(getNumRows() == sol.row_dual.size());
    for (auto i = 0; i < getNumRows(); i++) {
        sol.row_value[i] = rowprice[i];
    }
    int errcode = COPT_SetLpSolution(prob, &sol.col_value[0], &sol.row_value[0], &sol.row_dual[0], &sol.col_dual[0]);// 
    //int COPT_GetLpSolution(copt_prob * prob, double* value, double * slack(activity of constraints), double* rowDual(row price), double* redCost)
    if (errcode) {
        error_handling(errcode);
    }
}

void OsiCOPTSolverInterface::setColSolution(const double* colsol) {
    assert(getNumCols() == sol.col_value.size());
    for (auto i = 0; i < getNumCols(); i++) {
        sol.col_value[i] = colsol[i];
    }
    int errcode = COPT_SetLpSolution(prob, &sol.col_value[0], &sol.row_value[0], &sol.row_dual[0], &sol.col_dual[0]);// 
    //int COPT_GetLpSolution(copt_prob * prob, double* value, double * slack(activity of constraints), double* rowDual(row price), double* redCost)
    if (errcode) {
        error_handling(errcode);
    }
   
}

std::vector OsiCOPTSolverInterface::getDualRays(HighsInt maxNumRays,
    bool fullRay) const {
    assert(0);
    return std::vector(0);
}

std::vector OsiCOPTSolverInterface::getPrimalRays(
    HighsInt maxNumRays) const {
    assert(0);
    return std::vector(0);
}


#endif