File: otherclasses.xml

package info (click to toggle)
coinor-cbc 2.9.9+repack1-1
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 7,848 kB
  • ctags: 5,787
  • sloc: cpp: 104,337; sh: 8,921; xml: 2,950; makefile: 520; ansic: 491; awk: 197
file content (816 lines) | stat: -rw-r--r-- 36,092 bytes parent folder | download | duplicates (6)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
<?xml version="1.0" encoding="ISO-8859-1"?>
  <chapter id="otherclasses">
  <title>
  Selecting the Next Node in the Search Tree
  </title>
  <section id="comparison">
  <title>CbcCompare - Comparison Methods</title>
  <para>
  The order in which the nodes of the search tree are explored can strongly influence the performance of branch-and-cut algorithms. CBC give users complete control over the search order. The search order is controlled via the <classname>CbcCompare...</classname> class. CBC provides an abstract base class, <classname>CbcCompareBase</classname>,  and several commonly used instances which are described in <xref linkend="compareTable"/>.
  </para>
    <table frame="none" id="compareTable">
  <title>Compare Classes Provided</title>
    <tgroup cols="2">
    <thead>
    <row>
    <entry>
    Class name 
    </entry>
    <entry>
    Description
    </entry>
    </row>
    </thead>
    <tbody>
    <row>
      <entry align="left" valign="top">
      <classname>CbcCompareDepth</classname>
      </entry>
      <entry align="left" valign="top">
      This will always choose the node deepest in tree.  It gives minimum
      tree size but may take a long time to find the best solution.
      </entry>
    </row>
    <row>
      <entry align="left" valign="top">
      <classname>CbcCompareObjective</classname>
      </entry>
      <entry align="left" valign="top">
      This will always choose the node with the best objective value.  This may
      give a very large tree.  It is likely that the first solution found
      will be the best and the search should finish soon after the first solution
      is found.
      </entry>
    </row>
    <row>
      <entry align="left" valign="top">
      <classname>CbcCompareDefault</classname>
      </entry>
      <entry align="left" valign="top">
      This is designed to do a mostly depth-first search until a solution has
      been found. It then use estimates that are designed to give a slightly better solution.
      If a reasonable number of nodes have been explored (or a reasonable number of
      solutions found), then this class will adopt a breadth-first search (i.e., making a comparison based strictly on objective function values) unless the tree is very large when it will revert to depth-first search. A better description of <classname>CbcCompareUser</classname> is given below.
      </entry>
    </row>
    <row>
      <entry align="left" valign="top">
      <classname>CbcCompareEstimate</classname>
      </entry>
      <entry align="left" valign="top">
      When pseudo costs are invoked, they can be used to guess a solution.  This class uses the guessed solution.
      </entry>
    </row>
    </tbody>
  </tgroup>
  </table>
  <para>
  It is relatively simple for an experienced user to create new compare class instances. The code in <xref linkend="test"/> describes how to build a new comparison class and the reasoning behind it. The complete source can be found in <filename>CbcCompareUser.hpp</filename> and <filename>CbcCompareUser.cpp</filename>, located in the CBC Samples directory. See <xref linkend="moreexamples"/>. The key method in <classname>CbcCompare</classname> is <function>bool test(CbcNode* x, CbcNode* y))</function> which returns <parameter>true</parameter> if node <parameter>y</parameter> is preferred over node <parameter>x</parameter>. In the <function>test()</function> method, information from <classname>CbcNode</classname> can easily be used. <xref linkend="nodeTable"/> list some commonly used methods to access information at a node. 
  <table frame="none" id="nodeTable">
  <title>Information Available from <classname>CbcNode</classname></title>
    <tgroup cols="2">
    <tbody>
    <row>
      <entry align="left" valign="top">
      <function>double objectiveValue() const</function>
      </entry>
      <entry align="left" valign="top">
      Value of objective at the node.
      </entry>
    </row>
    <row>
      <entry align="left" valign="top">
      <function>int numberUnsatisfied() const</function>
      </entry>
      <entry align="left" valign="top">
      Number of unsatisfied integers (assuming branching 
      object is an integer - otherwise it might be number of unsatisfied sets).
      </entry>
    </row>
    <row>
      <entry align="left" valign="top">
      <function>int depth() const</function>
      </entry>
      <entry align="left" valign="top">
       Depth of the node in the search tree.
      </entry>
    </row>
    <row>
      <entry align="left" valign="top">
      <function>double guessedObjectiveValue() const</function>
      </entry>
     <entry align="left" valign="top"> 
     If user was setting this (e.g., if using pseudo costs).
      </entry>
    </row>
    <row>
      <entry align="left" valign="top">
      <function>int way() const</function>
      </entry>
      <entry align="left" valign="top">
       The way which branching would next occur from this node
       (for more advanced use).
      </entry>
    </row>
    <row>
      <entry align="left" valign="top">
      <function>int variable() const</function>
      </entry>
      <entry align="left" valign="top">
       The branching "variable" (associated with the <classname>CbcBranchingObject</classname> -- for more advanced use).
      </entry>
    </row>
    </tbody>
  </tgroup>
  </table>
</para>

<para>
The node desired in the tree is often a function of the how the search is progressing. In the design of CBC, there is no information on the state of the tree. The CBC is designed so that the method 
  <function>newSolution()</function> is called whenever a solution is found and the method <function>every1000Nodes()</function> is called every 1000 nodes.  When these methods are called, the user has the opportunity to modify the 
  behavior of <function>test()</function> by adjusting their common variables (e.g., <varname>weight_</varname>). Because <classname>CbcNode</classname> has a pointer to the model, the user can also influence the search through actions such as changing the maximum time CBC is allowed, once a solution has been found (e.g., <function>CbcModel::setMaximumSeconds(double value)</function>). In <filename>CbcCompareUser.cpp</filename> of the <filename>COIN/Cbc/Samples</filename> directory,  four items of data are used.
</para>
<itemizedlist>
  <listitem>
  <para>
1) The number of solutions found so far 
  </para>
  </listitem>
  <listitem>
  <para>
2) The size of the tree (defined to be the number of active nodes)
  </para>
  </listitem>
  <listitem>
  <para>
3) A weight, <varname>weight_</varname>, which is initialized to -1.0
  </para>
  </listitem>
  <listitem>
  <para>
4) A saved value of weight, <varname>saveWeight_</varname> (for when weight is set back to -1.0 for special reason)
  </para>
  </listitem>
</itemizedlist>
<para>
The full code for the <function>CbcCompareUser::test()</function> method is given in <xref linkend="test"/>.
</para>
  <example id="test">
  <title><function>CbcCompareUser::test()</function></title>
  <programlisting>
  <![CDATA[  
// Returns true if y better than x
bool 
CbcCompareUser::test (CbcNode * x, CbcNode * y)
{
  if (weight_==-1.0) {
    // before solution
    if (x->numberUnsatisfied() > y->numberUnsatisfied())
      return true;
    else if (x->numberUnsatisfied() < y->numberUnsatisfied())
      return false;
    else
      return x->depth() < y->depth();
  } else {
    // after solution. 
    // note: if weight_=0, comparison is based 
    //       solely on objective value
    double weight = CoinMax(weight_,0.0);
    return x->objectiveValue()+ weight*x->numberUnsatisfied() > 
      y->objectiveValue() + weight*y->numberUnsatisfied();
  }
}
  ]]>   
  </programlisting>
  </example>
<para>
Initially, <varname>weight</varname>_ is -1.0 and the search is biased towards depth first.  In 
fact, <function>test()</function> prefers <parameter>y</parameter> if <parameter>y</parameter> has fewer unsatisfied variables. In the case of a tie, <function>test()</function> prefers the node with the greater depth in tree. Once a solution is found, <function>newSolution()</function> is called. The method <function>newSolution()</function> interacts with <function>test()</function> by means of the variable <varname>weight_</varname>. If the solution was achieved by branching, <!-- how do can you determine that? --> a calculation is made to determine the cost per unsatisfied integer variable to go from the continuous solution to an integer solution.  The variable <varname>weight_</varname> is then set to aim at a slightly better solution.  From then on, <function>test()</function> returns <parameter>true</parameter> if it seems that <parameter>y</parameter> will lead to a better solution than <parameter>x</parameter>. This source for <function>newSolution()</function> in given in <xref linkend="newSolution"/>.
</para>
  <example id="newSolution">
  <title><function>CbcCompareUser::newSolution()</function></title>
  <programlisting>
  <![CDATA[  
// This allows the test() method to change behavior by resetting weight_.
// It is called after each new solution is found.
void 
CbcCompareUser::newSolution(CbcModel * model,
			       double objectiveAtContinuous,
			       int numberInfeasibilitiesAtContinuous) 
{
  if (model->getSolutionCount()==model->getNumberHeuristicSolutions())
    return; // solution was found by rounding so ignore it.

  // set weight_ to get close to this solution
  double costPerInteger = 
    (model->getObjValue()-objectiveAtContinuous)/
    ((double) numberInfeasibilitiesAtContinuous);
  weight_ = 0.98*costPerInteger;
  saveWeight_=weight_;
  numberSolutions_++;
  if (numberSolutions_>5)
    weight_ =0.0; // comparison in test() will be 
                  // based strictly on objective value.
}
  ]]>   
  </programlisting>
  </example>
<para>

As the search progresses, the comparison can be <!-- what? is "this"? -->modified. If many nodes (or many solutions) have been genereated, then <varname>weight_</varname> is set to 0.0 leading to a breadth-first search.  Breadth-first search can lead to an enormous tree. If the tree size is exceeds 10000, it may be desirable to return to a search biased towards depth first. Changing the behavior in this manner is done by the method <function>every1000Nodes</function> shown in <xref linkend="everyK"/>.
  </para>
  <example id="everyK">
  <title><function>CbcCompareUser::every1000Nodes()</function></title>
  <programlisting>
  <![CDATA[  
// This allows the test() method to change behavior every so often
bool 
CbcCompareUser::every1000Nodes(CbcModel * model, int numberNodes)
{
  if (numberNodes>10000)
    weight_ =0.0; // compare nodes based on objective value
  else if (numberNodes==1000&&weight_==-2.0)
    weight_=-1.0; // Go to depth first
  // get size of tree
  treeSize_ = model->tree()->size();
  if (treeSize_>10000) {
    // set weight to reduce size most of time
    if (treeSize_>20000)
      weight_=-1.0;
    else if ((numberNodes%4000)!=0)
      weight_=-1.0;
    else
      weight_=saveWeight_;
  }
  return numberNodes==11000; // resort if first time
}
  ]]>   
  </programlisting>
  </example>
  </section>
</chapter>
<chapter id="hueristicChap">
  <title>
  Getting Good Bounds in CBC
  </title>
  <section id="heuristics">
  <title>CbcHeuristic - Heuristic Methods</title>
  <para>
  In practice, it is very useful to get a good solution reasonably fast. Any MIP-feasible solution produces an upper bound, and a good bound will greatly reduce the run time. Good solutions can satisfy the user
  on very large problems where a complete search is impossible.  Obviously, heuristics are
  problem dependent, although some do have more general use.
  At present there is only one heuristic in CBC itself, <classname>CbcRounding</classname>. Hopefully, the number will grow. Other heuristics are in the <filename>COIN/Cbc/Samples</filename>
  directory.  A heuristic tries to obtain a solution to the original
  problem so it only needs to consider the original rows and does not have to use the
  current bounds. CBC provides an abstract base class <classname>CbcHeuristic</classname> and a rounding heuristic in CBC. <!-- rlh: we need an example of using an existing heuristic -->
  </para>
  <para>
  This chapter describes how to build a greedy heuristic for a set covering problem, e.g., the miplib problem fast0507. A more general (and efficient) version of the heuristic is in <filename>CbcHeuristicGreedy.hpp</filename> and <filename>CbcHeuristicGreedy.cpp</filename> located in the <filename>COIN/Cbc/Samples</filename> directory, see <xref linkend="moreexamples"/>.
</para>
<para>
  The greedy heuristic will leave all variables taking value one at this node of the
  tree at value one, and will initially set all other variable to value zero.  
  All variables are then sorted in order of their cost
  divided by the number of entries in rows which are not yet covered. (We may randomize that
  value a bit so that ties will be broken in different ways on different runs of the heuristic.)
  The best one is choosen, and set to one. The process is repeated. Because this is
  a set covering problem (i.e., all constraints are &ge;), the heuristic is guaranteed to find a solution (but not necessarily an improved solution). The speed of the heuristic could be improved by just redoing those affected, but for illustrative purposes we will keep it simple.(The speed could also be improved if all elements are 1.0). 
</para>
<para>
  The key <classname>CbcHeuristic</classname> method is <function>int&nbsp;solution(double &amp; solutionValue,
                                              double&nbsp;*&nbsp;betterSolution)</function>.
  The <function>solution()</function> method returns 0 if no solution found, and returns 1 if a solution is found, in which case it fills in the objective value and primal solution.  The code in <filename>CbcHeuristicGreedy.cpp</filename> is a little more complicated than this following example. For instance, the code here assumes all variables are integer.  The important bit of data is a copy of the matrix (stored by column) before any cuts have been made.  The data used are bounds, objective and the matrix plus two work arrays. 
  </para>
  <example>
  <title>Data</title>
  <programlisting>
  <![CDATA[  
  OsiSolverInterface * solver = model_->solver(); // Get solver from CbcModel
  const double * columnLower = solver->getColLower(); // Column Bounds
  const double * columnUpper = solver->getColUpper();
  const double * rowLower = solver->getRowLower(); // We know we only need lower bounds
  const double * solution = solver->getColSolution();
  const double * objective = solver->getObjCoefficients(); // In code we also use min/max
  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
  double primalTolerance;
  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
  int numberRows = originalNumberRows_; // This is number of rows when matrix was passed in
  // Column copy of matrix (before cuts)
  const double * element = matrix_.getElements();
  const int * row = matrix_.getIndices();
  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
  const int * columnLength = matrix_.getVectorLengths();

  // Get solution array for heuristic solution
  int numberColumns = solver->getNumCols();
  double * newSolution = new double [numberColumns];
  // And to sum row activities
  double * rowActivity = new double[numberRows];
  ]]>   
  </programlisting>
  </example>
<para>
The <varname>newSolution</varname> is then initialized to the rounded down solution.
</para>
  <example>
  <title>Initialize <varname>newSolution</varname></title>
  <programlisting>
  <![CDATA[  
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    CoinBigIndex j;
    double value = solution[iColumn];
    // Round down integer
    if (fabs(floor(value+0.5)-value)<integerTolerance) 
      value=floor(CoinMax(value+1.0e-3,columnLower[iColumn]));
    // make sure clean
    value = CoinMin(value,columnUpper[iColumn]);
    value = CoinMax(value,columnLower[iColumn]);
    newSolution[iColumn]=value;
    if (value) {
      double cost = direction * objective[iColumn];
      newSolutionValue += value*cost;
      for (j=columnStart[iColumn];
           j<columnStart[iColumn]+columnLength[iColumn];j++) {
        int iRow=row[j];
        rowActivity[iRow] += value*element[j];
      }
    }
  }
  ]]>   
  </programlisting>
  </example>
<para>
<!-- rlh: where is "direction" from? -->

At this point some row activities may be below their lower bound. To correct this infeasibility, the variable which is cheapest in reducing the sum of infeasibilities is found and updated, and the process repeats.  This is a finite process. (Theimplementation could be faster, but is kept simple for illustrative purposes.)
  </para>
  <example>
  <title>Create Feasible <varname>newSolution</varname> from Initial <varname>newSolution</varname></title>
  <programlisting>
  <![CDATA[  
  while (true) {
    // Get column with best ratio
    int bestColumn=-1;
    double bestRatio=COIN_DBL_MAX;
    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
      CoinBigIndex j;
      double value = newSolution[iColumn];
      double cost = direction * objective[iColumn];
      // we could use original upper rather than current
      if (value+0.99<columnUpper[iColumn]) {
        double sum=0.0; // Compute how much we will reduce infeasibility by
        for (j=columnStart[iColumn];
             j<columnStart[iColumn]+columnLength[iColumn];j++) {
          int iRow=row[j];
          double gap = rowLower[iRow]-rowActivity[iRow];
          if (gap>1.0e-7) {
            sum += CoinMin(element[j],gap);
          if (element[j]+rowActivity[iRow]<rowLower[iRow]+1.0e-7) {
            sum += element[j];
	  }
	}
        if (sum>0.0) {
	  double ratio = (cost/sum)*(1.0+0.1*CoinDrand48());
	  if (ratio<bestRatio) {
            bestRatio=ratio;
            bestColumn=iColumn;
          }
	}
      }
    }
    if (bestColumn<0)
      break; // we have finished
    // Increase chosen column
    newSolution[bestColumn] += 1.0;
    double cost = direction * objective[bestColumn];
    newSolutionValue += cost;
    for (CoinBigIndex j=columnStart[bestColumn];
         j<columnStart[bestColumn]+columnLength[bestColumn];j++) {
      int iRow = row[j];
      rowActivity[iRow] += element[j];
    }
  }
  ]]>   
  </programlisting>
  </example>
<para>
A solution value of <varname>newSolution</varname> is compared to the best solution value. If <varname>newSolution</varname> is an improvement, its feasibility is validated.
  </para>
  <example>
  <title>Check Solution Quality of <varname>newSolution</varname></title>
  <programlisting>
  <![CDATA[  
  returnCode=0; // 0 means no good solution
  if (newSolutionValue<solutionValue) { // minimization
    // check feasible
    memset(rowActivity,0,numberRows*sizeof(double));
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      CoinBigIndex j;
      double value = newSolution[iColumn];
      if (value) {
	for (j=columnStart[iColumn];
	     j<columnStart[iColumn]+columnLength[iColumn];j++) {
	  int iRow=row[j];
	  rowActivity[iRow] += value*element[j];
	}
      }
    }
    // check was approximately feasible
    bool feasible=true;
    for (iRow=0;iRow<numberRows;iRow++) {
      if(rowActivity[iRow]<rowLower[iRow]) {
	if (rowActivity[iRow]<rowLower[iRow]-10.0*primalTolerance)
	  feasible = false;
      }
    }
    if (feasible) {
      // new solution
      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
      solutionValue = newSolutionValue;
      // We have good solution
      returnCode=1;
    }
  }
  ]]>   
  </programlisting>
  </example>
  </section>
</chapter>

<chapter id="branchChapter">
  <title>
  Branching
 </title>
  <section id="branching">
  <title>Pseudo Cost Branching</title>
  <para>
<!-- rlh: need a description of "objects" -->
If the user declares variables as integer but does no more, then CBC will treat them
as simple integer variables.  In many cases the user would like to do some more fine tuning. This section shows how to create integer variables with pseudo costs.  When pseudo costs are given then 
it is assumed that if a variable is at 1.3 then the cost of branching that variable down will be 0.3 times the down pseudo cost and the cost of branching up would be 0.7 times the up pseudo cost.  Pseudo costs can be used both for branching and for choosing a node.
   The full code is in <filename>longthin.cpp</filename> located in the CBC Samples directory, see
  <xref linkend="moreexamples"/>.
</para>
<para>  
  The idea is simple for set covering problems.
  Branching up gets us much closer to an integer solution so we will encourage that direction by branch up if variable value > 0.333333.
  The expected cost of going up obviously depends on the cost of the
  variable. The pseudo costs are choosen to reflect that fact.
<!-- rlh: what's the difference between simple var objects and pseudo cost objects? explain. -->
  </para>
  <example id="pseudo">
  <title><classname>CbcSimpleIntegerPseudoCosts</classname></title>
  <programlisting>
  <![CDATA[  
  int iColumn;
  int numberColumns = solver3->getNumCols();
  // do pseudo costs
  CbcObject ** objects = new CbcObject * [numberColumns];
  // Point to objective
  const double * objective = model.getObjCoefficients();
  int numberIntegers=0;
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    if (solver3->isInteger(iColumn)) {
      double cost = objective[iColumn];
      CbcSimpleIntegerPseudoCost * newObject =
        new CbcSimpleIntegerPseudoCost(&model,numberIntegers,iColumn,
                                       2.0*cost,cost);
      newObject->setMethod(3);
      objects[numberIntegers++]= newObject;
    }
  }
  // Now add in objects (they will replace simple integers)
  model.addObjects(numberIntegers,objects);
  for (iColumn=0;iColumn<numberIntegers;iColumn++)
    delete objects[iColumn];
  delete [] objects;
  ]]>   
  </programlisting>
  </example>
<para>
The code in <xref linkend="pseudo"/> also tries to give more importance to variables with more 
coefficients.  Whether this sort of thing is worthwhile should be the subject of experimentation.
<!-- needs more explaination, e.g., why do you have to delete objects? -->
<!-- what does a pseudo cost object do for you? and how does it do it? -->
</para>
</section>
  <section id="followOn">
  <title>Follow-On Branching</title>
  <para>
In crew scheduling, the problems are long and thin. A problem may have a few rows but many thousands of variables.  Branching a variable to 1 is very powerful as it fixes many other variables to zero, but branching to zero is very weak as thousands of variables can increase from zero.  In crew scheduling problems, each constraint is a flight leg, e.g., JFK airport to DFW airport.
From DFW there may be several flights the crew could take next - suppose one flight is 
the 9:30 flight from DFW to LAX airport.  A binary branch is that the crew arriving 
at DFW either take the 9:30 flight to LAX or they don't.  This "follow-on" branching does not
fix individual variables. Instead this branching divides all the variables with entries in the JFK-DFW 
constraint into two groups - those with entries in the DFW-LAX constraint and those without
entries.
</para>
<para>
   The full sample code for follow-on brancing is in <filename>crew.cpp</filename>
  located in the CBC Samples directory, see
  <xref linkend="moreexamples"/>).  In this case, the simple integer
variables are left which may be necessary if other sorts of constraints exist. Follow-on branching rules are to be considered first, so the priorities are set to indicated the follow-on rules take precedence. Priority 1 is the highest priority.
<!-- rlh:need to explain how *priorities* work, and where they are used -->
</para>
  <example>
  <title><classname>CbcFollowOn</classname></title>
  <programlisting>
  <![CDATA[  
  int iColumn;
  int numberColumns = solver3->getNumCols();
  /* We are going to add a single follow-on object but we
     want to give low priority to existing integers 
     As the default priority is 1000 we don't actually need to give
     integer priorities but it is here to show how.
  */
  // Normal integer priorities
  int * priority = new int [numberColumns];
  int numberIntegers=0;
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    if (solver3->isInteger(iColumn)) {
      priority[numberIntegers++]= 100; // low priority
    }
  }
  /* Second parameter is set to true for objects,
     and false for integers. This indicates integers */
  model.passInPriorities(priority,false);
  delete [] priority;
  /* Add in objects before we can give them a priority.
     In this case just one object
     - but it shows the general method
  */
  CbcObject ** objects = new CbcObject * [1];
  objects[0]=new CbcFollowOn(&model);
  model.addObjects(1,objects);
  delete objects[0];
  delete [] objects;
  // High priority
  int followPriority=1;
  model.passInPriorities(&followPriority,true);
  ]]>   
  </programlisting>
  </example>
<!-- rlh: need to explain the CbcFollowOn object -->
  </section>
</chapter>
<chapter id="SolverChap">
<title>
  Advance Solver Uses
</title>
  <section id="solver">
  <title>Creating a Solver via Inheritance</title>
  <para>
  CBC uses a generic <classname>OsiSolverInterface</classname> and its <function>resolve</function> capability.
  This does not give much flexibility so advanced users can inherit from their interface
  of choice.  This section illustrates how to implement such a solver for a long thin problem, e.g., fast0507 again.  As with the other examples in the Guide, the sample code is not guaranteed to be the fastest way to solve the problem. The main purpose of the example is to illustrate techniques. The full source is in <filename>CbcSolver2.hpp</filename> and <filename>CbcSolver2.cpp</filename> located in the CBC Samples directory, see
  <xref linkend="moreexamples"/>. 
</para>
<para>
The method <function>initialSolve</function> is called a few times in CBC, and provides a convenient starting point. The <varname>modelPtr_</varname> derives from <classname>OsiClpSolverInterface</classname>.
  </para>
  <example id="initialSolve">
  <title><function>initialSolve()</function></title>
  <programlisting>
  <![CDATA[  
  // modelPtr_ is of type ClpSimplex *
  modelPtr_->setLogLevel(1); // switch on a bit of printout
  modelPtr_->scaling(0); // We don't want scaling for fast0507
  setBasis(basis_,modelPtr_); // Put basis into ClpSimplex
  // Do long thin by sprint
  ClpSolve options;
  options.setSolveType(ClpSolve::usePrimalorSprint);
  options.setPresolveType(ClpSolve::presolveOff);
  options.setSpecialOption(1,3,15); // Do 15 sprint iterations
  modelPtr_->initialSolve(options); // solve problem
  basis_ = getBasis(modelPtr_); // save basis
  modelPtr_->setLogLevel(0); // switch off printout
  ]]>   
  </programlisting>
  </example>
<!-- rlh: need a verbal description to go with the intialSolve() code -->
<para>
The <function>resolve()</function> method is more complicated than <function>initialSolve()</function>.  The main pieces of data are a counter <varname>count_</varname> which is incremented each solve and an integer array <varname>node_</varname> which stores the last time
a variable was active in a solution.  For the first few times, the normal Dual Simplex is called and
<varname>node_</varname> array is updated.
</para>
  <example>
  <title>First Few Solves</title>
  <programlisting>
  <![CDATA[  
  if (count_<10) {
    OsiClpSolverInterface::resolve(); // Normal resolve
    if (modelPtr_->status()==0) {
      count_++; // feasible - save any nonzero or basic
      const double * solution = modelPtr_->primalColumnSolution();
      for (int i=0;i<numberColumns;i++) {
	if (solution[i]>1.0e-6||modelPtr_->getStatus(i)==ClpSimplex::basic) {
	  node_[i]=CoinMax(count_,node_[i]);
	  howMany_[i]++;
	}
      }
    } else {
      printf("infeasible early on\n");
    }
  }
  ]]>   
  </programlisting>
  </example>
<!-- rlh: need more explanation of what's being illustrated in the code -->
<para>
After the first few solves, only those variables which took part in a solution in the last so many
solves are used.  As fast0507 is a set covering problem, any rows which are already covered can be taken out.
  </para>
  <example>
  <title>Create Small Problem</title>
  <programlisting>
  <![CDATA[  
    int * whichRow = new int[numberRows]; // Array to say which rows used
    int * whichColumn = new int [numberColumns]; // Array to say which columns used
    int i;
    const double * lower = modelPtr_->columnLower();
    const double * upper = modelPtr_->columnUpper();
    setBasis(basis_,modelPtr_); // Set basis
    int nNewCol=0; // Number of columns in small model
    // Column copy of matrix
    const double * element = modelPtr_->matrix()->getElements();
    const int * row = modelPtr_->matrix()->getIndices();
    const CoinBigIndex * columnStart = modelPtr_->matrix()->getVectorStarts();
    const int * columnLength = modelPtr_->matrix()->getVectorLengths();
    
    int * rowActivity = new int[numberRows]; // Number of columns with entries in each row
    memset(rowActivity,0,numberRows*sizeof(int));
    int * rowActivity2 = new int[numberRows]; // Lower bound on row activity for each row
    memset(rowActivity2,0,numberRows*sizeof(int));
    char * mark = (char *) modelPtr_->dualColumnSolution(); // Get some space to mark columns
    memset(mark,0,numberColumns);
    for (i=0;i<numberColumns;i++) {
      bool choose = (node_[i]>count_-memory_&&node_[i]>0); // Choose if used recently
      // Take if used recently or active in some sense
      if ((choose&&upper[i])
	  ||(modelPtr_->getStatus(i)!=ClpSimplex::atLowerBound&&
             modelPtr_->getStatus(i)!=ClpSimplex::isFixed)
	  ||lower[i]>0.0) {
        mark[i]=1; // mark as used
	whichColumn[nNewCol++]=i; // add to list
        CoinBigIndex j;
        double value = upper[i];
        if (value) {
          for (j=columnStart[i];
               j<columnStart[i]+columnLength[i];j++) {
            int iRow=row[j];
            assert (element[j]==1.0);
            rowActivity[iRow] ++; // This variable can cover this row
          }
          if (lower[i]>0.0) {
            for (j=columnStart[i];
                 j<columnStart[i]+columnLength[i];j++) {
              int iRow=row[j];
              rowActivity2[iRow] ++; // This row redundant
            }
          }
        }
      }
    }
    int nOK=0; // Use to count rows which can be covered
    int nNewRow=0; // Use to make list of rows needed
    for (i=0;i<numberRows;i++) {
      if (rowActivity[i])
        nOK++;
      if (!rowActivity2[i])
        whichRow[nNewRow++]=i; // not satisfied
      else
        modelPtr_->setRowStatus(i,ClpSimplex::basic); // make slack basic
    }
    if (nOK<numberRows) {
      // The variables we have do not cover rows - see if we can find any that do
      for (i=0;i<numberColumns;i++) {
        if (!mark[i]&&upper[i]) {
          CoinBigIndex j;
          int good=0;
          for (j=columnStart[i];
               j<columnStart[i]+columnLength[i];j++) {
            int iRow=row[j];
            if (!rowActivity[iRow]) {
              rowActivity[iRow] ++;
              good++;
            }
          }
          if (good) {
            nOK+=good; // This covers - put in list
            whichColumn[nNewCol++]=i;
          }
        }
      }
    }
    delete [] rowActivity;
    delete [] rowActivity2;
    if (nOK<numberRows) {
      // By inspection the problem is infeasible - no need to solve
      modelPtr_->setProblemStatus(1);
      delete [] whichRow;
      delete [] whichColumn;
      printf("infeasible by inspection\n");
      return;
    }
    // Now make up a small model with the right rows and columns
    ClpSimplex *  temp = new ClpSimplex(modelPtr_,nNewRow,whichRow,nNewCol,whichColumn);
  ]]>   
  </programlisting>
  </example>
<para>
If the variables cover the rows, then the problem is feasible (no cuts are being used). If the rows
were equality constraints, then this might not be the case. More work would be needed.  After the solution, the reduct costs are checked. If any reduced costs are negative, the code goes back to the full problem and cleans up with Primal Simplex.
  </para>
  <example>
  <title>Check Optimal Solution</title>
  <programlisting>
  <![CDATA[  
    temp->setDualObjectiveLimit(1.0e50); // Switch off dual cutoff as problem is restricted
    temp->dual(); // solve
    double * solution = modelPtr_->primalColumnSolution(); // put back solution
    const double * solution2 = temp->primalColumnSolution();
    memset(solution,0,numberColumns*sizeof(double));
    for (i=0;i<nNewCol;i++) {
      int iColumn = whichColumn[i];
      solution[iColumn]=solution2[i];
      modelPtr_->setStatus(iColumn,temp->getStatus(i));
    }
    double * rowSolution = modelPtr_->primalRowSolution();
    const double * rowSolution2 = temp->primalRowSolution();
    double * dual = modelPtr_->dualRowSolution();
    const double * dual2 = temp->dualRowSolution();
    memset(dual,0,numberRows*sizeof(double));
    for (i=0;i<nNewRow;i++) {
      int iRow=whichRow[i];
      modelPtr_->setRowStatus(iRow,temp->getRowStatus(i));
      rowSolution[iRow]=rowSolution2[i];
      dual[iRow]=dual2[i];
    }
    // See if optimal
    double * dj = modelPtr_->dualColumnSolution();
    // get reduced cost for large problem
    // this assumes minimization
    memcpy(dj,modelPtr_->objective(),numberColumns*sizeof(double));
    modelPtr_->transposeTimes(-1.0,dual,dj);
    modelPtr_->setObjectiveValue(temp->objectiveValue());
    modelPtr_->setProblemStatus(0);
    int nBad=0;
      
    for (i=0;i<numberColumns;i++) {
      if (modelPtr_->getStatus(i)==ClpSimplex::atLowerBound
          &&upper[i]>lower[i]&&dj[i]<-1.0e-5)
        nBad++;
    }
    // If necessary clean up with primal (and save some statistics)
    if (nBad) {
      timesBad_++;
      modelPtr_->primal(1);
      iterationsBad_ += modelPtr_->numberIterations();
    }
  ]]>   
  </programlisting>
  </example>
<para>
The array <varname>node_</varname> is updated, as for the first few solves.  To give some idea of the effect of this tactic, the problem fast0507 has 63,009 variables but the small problem never has more than 4,000 variables. In only about ten percent of solves was it necessary to resolve, and then the average number of iterations
on full problem was less than 20.
</para>
</section>
 <section id="quadratic">
  <title>Quadratic MIP</title>
<para>
To give another example - again only for illustrative purposes -- it is possible to do quadratic
MIP with CBC.  In this case, we make <function>resolve</function> the same as 
<function>initialSolve</function>.
   The full code is in <filename>ClpQuadInterface.hpp</filename> and
   <filename>ClpQuadInterface.cpp</filename> located in the CBC Samples directory, see
  <xref linkend="moreexamples"/>).
  </para>
  <example>
  <title>Solving a Quadratic MIP</title>
  <programlisting>
  <![CDATA[  
  // save cutoff
  double cutoff = modelPtr_->dualObjectiveLimit();
  modelPtr_->setDualObjectiveLimit(1.0e50);
  modelPtr_->scaling(0);
  modelPtr_->setLogLevel(0);
  // solve with no objective to get feasible solution
  setBasis(basis_,modelPtr_);
  modelPtr_->dual();
  basis_ = getBasis(modelPtr_);
  modelPtr_->setDualObjectiveLimit(cutoff);
  if (modelPtr_->problemStatus()) 
    return; // problem was infeasible 
  // Now pass in quadratic objective
  ClpObjective * saveObjective  = modelPtr_->objectiveAsObject();
  modelPtr_->setObjectivePointer(quadraticObjective_);
  modelPtr_->primal();
  modelPtr_->setDualObjectiveLimit(cutoff);
  if (modelPtr_->objectiveValue()>cutoff)
    modelPtr_->setProblemStatus(1);
  modelPtr_->setObjectivePointer(saveObjective);
  ]]>   
  </programlisting>
  </example>
  </section>
  </chapter>