Changeset 1017


Ignore:
Timestamp:
Jan 16, 2008 7:21:58 PM (12 years ago)
Author:
pbelotti
Message:

rewrote expr{Quad,Group} with vectors of exprVar*. Changed interface to procedures for independence of CouenneProblem?. Use variables instead of indices in branching objects. Changed integrality check. Fixed yet another standardization bug. Fixed reading of integer variables' indices. Fixed loose bounds in CouenneProblem::getIntegerCandidate()

Files:
4 added
117 edited
1 moved

Legend:

Unmodified
Added
Removed
  • branches/Couenne/Couenne/Makefile.in

    r976 r1017  
    2121# Author:  Andreas Waechter           IBM    2006-04-13
    2222
    23 # Copyright (C) 2006 International Business Machines and others.
     23# Copyright (C) 2006, 2007 International Business Machines and others.
    2424# All Rights Reserved.
    2525# This file is distributed under the Common Public License.
     
    6363        $(srcdir)/Makefile.am $(srcdir)/Makefile.in \
    6464        $(srcdir)/couenne_addlibs.txt.in $(top_srcdir)/configure \
    65         $(top_srcdir)/inc/config_couenne.h.in AUTHORS COPYING TODO
     65        $(top_srcdir)/inc/config_couenne.h.in AUTHORS COPYING
    6666@HAVE_EXTERNALS_TRUE@am__append_2 = Externals
    6767@HAVE_EXTERNALS_TRUE@am__append_3 = .Externals-stamp
     
    297297DISTCLEANFILES = $(am__append_3) $(VPATH_DISTCLEANFILES)
    298298DocFiles = README AUTHORS LICENSE
    299 DocInstallDir = $(prefix)/share/doc/$(PACKAGE)
     299DocInstallDir = $(prefix)/share/doc/coin/$(PACKAGE_NAME)
    300300all: all-recursive
    301301
     
    690690install-data-am:
    691691
     692install-exec-am: install-exec-local
     693
    692694install-info: install-info-recursive
    693695
     
    713715
    714716ps-am:
     717
     718uninstall-am: uninstall-info-am uninstall-local
    715719
    716720uninstall-info: uninstall-info-recursive
     
    757761        rm -f $(DESTDIR)$(libdir)/$(addlibsfile)
    758762
    759 install-exec-am:
     763install-doc: $(DocFiles)
    760764        test -z "$(DocInstallDir)" || $(mkdir_p) "$(DESTDIR)$(DocInstallDir)"
    761765        for file in $(DocFiles); do \
     
    764768        done
    765769
    766 uninstall-am:
     770uninstall-doc:
    767771        for file in $(DocFiles); do \
    768772          rm -f "$(DESTDIR)$(DocInstallDir)/$$file"; \
     
    805809@HAVE_EXTERNALS_TRUE@@MAINTAINER_MODE_TRUE@update-externals: .Externals-stamp
    806810@HAVE_EXTERNALS_TRUE@@MAINTAINER_MODE_TRUE@     cd $(srcdir); svn update
     811
     812.PHONY: install-doc uninstall-doc update-externals
    807813# Tell versions [3.59,3.63) of GNU make to not export all variables.
    808814# Otherwise a system limit (for SysV at least) may be exceeded.
  • branches/Couenne/Couenne/configure

    r915 r1017  
    10791079  --enable-doscompile     Under Cygwin, compile so that executables run under
    10801080                          DOS. Set to mingw to use gcc/g++/ld with
    1081                           -mno-cygwin. Set to msvc to use cl/link. Default
    1082                           when mentioned: mingw. Default when not mentioned:
    1083                           disabled.
     1081                          -mno-cygwin. Set to msvc to use cl/link (or
     1082                          icl/link). Default when mentioned: mingw. Default
     1083                          when not mentioned: disabled.
    10841084  --enable-static[=PKGS]
    10851085                          build static libraries [default=no]
     
    18561856if test x"$CXX" != x; then
    18571857  case "$CXX" in
    1858     cl* | */cl* | CL* | */CL*)
     1858    cl* | */cl* | CL* | */CL* | icl* | */icl* | ICL* | */ICL*)
    18591859      if test x"$CC" = x; then
    18601860        CC="$CXX"
     
    18831883  *-cygwin* | *-mingw*)
    18841884             if test "$enable_doscompile" = msvc ; then
    1885                comps="cl"
     1885               comps="icl cl"
    18861886             else
    18871887               comps="gcc cl"
     
    27692769            coin_dbg_cflags='-MTd'
    27702770            ;;
     2771          icl* | */icl* | ICL* | */ICL*)
     2772            coin_opt_cflags='-MT -Ox'
     2773            coin_add_cflags='-nologo -D_CRT_SECURE_NO_DEPRECATE'
     2774            coin_dbg_cflags='-MTd -debug'
     2775            ;;
    27712776        esac
    27722777        ;;
     
    30753080# Correct ADDLIBS initialization if we are using the MS compiler
    30763081case "$CC" in
    3077   cl* | */cl* | CL* | */CL*)
     3082  cl* | */cl* | CL* | */CL* | icl* | */icl* | ICL* | */ICL*)
    30783083    ADDLIBS=
    30793084    case $build in
     
    31243129  *-cygwin* | *-mingw*)
    31253130             if test "$enable_doscompile" = msvc ; then
    3126                comps="cl"
     3131               comps="icl cl"
    31273132             else
    31283133               comps="g++ cl"
     
    35013506# It seems that we need to cleanup something here for the Windows
    35023507case "$CXX" in
    3503   cl* | */cl* | CL* | */CL* )
     3508  cl* | */cl* | CL* | */CL* | icl* | */icl* | ICL* | */ICL*)
    35043509    sed -e 's/^void exit (int);//' confdefs.h >> confdefs.hh
    35053510    mv confdefs.hh confdefs.h
     
    36093614            coin_add_cxxflags='-nologo -EHsc -GR -wd4996 -D_CRT_SECURE_NO_DEPRECATE'
    36103615            coin_dbg_cxxflags='-MTd'
     3616            ;;
     3617          icl* | */icl* | ICL* | */ICL*)
     3618            # The MT and MTd options are mutually exclusive
     3619            coin_opt_cxxflags='-MT -Ox'
     3620            coin_add_cxxflags='-nologo -EHsc -GR -D_CRT_SECURE_NO_DEPRECATE'
     3621            coin_dbg_cxxflags='-MTd -debug'
    36113622            ;;
    36123623        esac
     
    39623973if test x"$CXX" != x; then
    39633974  case "$CXX" in
    3964     cl* | */cl* | CL* | */CL*)
     3975    cl* | */cl* | CL* | */CL* | icl* | */icl* | ICL* | */ICL*)
    39653976      if test x"$CC" = x; then
    39663977        CC="$CXX"
     
    39894000  *-cygwin* | *-mingw*)
    39904001             if test "$enable_doscompile" = msvc ; then
    3991                comps="cl"
     4002               comps="icl cl"
    39924003             else
    39934004               comps="gcc cl"
     
    46744685            coin_dbg_cflags='-MTd'
    46754686            ;;
     4687          icl* | */icl* | ICL* | */ICL*)
     4688            coin_opt_cflags='-MT -Ox'
     4689            coin_add_cflags='-nologo -D_CRT_SECURE_NO_DEPRECATE'
     4690            coin_dbg_cflags='-MTd -debug'
     4691            ;;
    46764692        esac
    46774693        ;;
     
    49804996# Correct ADDLIBS initialization if we are using the MS compiler
    49814997case "$CC" in
    4982   cl* | */cl* | CL* | */CL*)
     4998  cl* | */cl* | CL* | */CL* | icl* | */icl* | ICL* | */ICL*)
    49834999    ADDLIBS=
    49845000    case $build in
     
    67406756*-*-irix6*)
    67416757  # Find out which ABI we are using.
    6742   echo '#line 6742 "configure"' > conftest.$ac_ext
     6758  echo '#line 6758 "configure"' > conftest.$ac_ext
    67436759  if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5
    67446760  (eval $ac_compile) 2>&5
     
    78747890
    78757891# Provide some information about the compiler.
    7876 echo "$as_me:7876:" \
     7892echo "$as_me:7892:" \
    78777893     "checking for Fortran 77 compiler version" >&5
    78787894ac_compiler=`set X $ac_compile; echo $2`
     
    89418957   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    89428958   -e 's:$: $lt_compiler_flag:'`
    8943    (eval echo "\"\$as_me:8943: $lt_compile\"" >&5)
     8959   (eval echo "\"\$as_me:8959: $lt_compile\"" >&5)
    89448960   (eval "$lt_compile" 2>conftest.err)
    89458961   ac_status=$?
    89468962   cat conftest.err >&5
    8947    echo "$as_me:8947: \$? = $ac_status" >&5
     8963   echo "$as_me:8963: \$? = $ac_status" >&5
    89488964   if (exit $ac_status) && test -s "$ac_outfile"; then
    89498965     # The compiler can only warn and ignore the option if not recognized
     
    92099225   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    92109226   -e 's:$: $lt_compiler_flag:'`
    9211    (eval echo "\"\$as_me:9211: $lt_compile\"" >&5)
     9227   (eval echo "\"\$as_me:9227: $lt_compile\"" >&5)
    92129228   (eval "$lt_compile" 2>conftest.err)
    92139229   ac_status=$?
    92149230   cat conftest.err >&5
    9215    echo "$as_me:9215: \$? = $ac_status" >&5
     9231   echo "$as_me:9231: \$? = $ac_status" >&5
    92169232   if (exit $ac_status) && test -s "$ac_outfile"; then
    92179233     # The compiler can only warn and ignore the option if not recognized
     
    93139329   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    93149330   -e 's:$: $lt_compiler_flag:'`
    9315    (eval echo "\"\$as_me:9315: $lt_compile\"" >&5)
     9331   (eval echo "\"\$as_me:9331: $lt_compile\"" >&5)
    93169332   (eval "$lt_compile" 2>out/conftest.err)
    93179333   ac_status=$?
    93189334   cat out/conftest.err >&5
    9319    echo "$as_me:9319: \$? = $ac_status" >&5
     9335   echo "$as_me:9335: \$? = $ac_status" >&5
    93209336   if (exit $ac_status) && test -s out/conftest2.$ac_objext
    93219337   then
     
    1165811674  lt_status=$lt_dlunknown
    1165911675  cat > conftest.$ac_ext <<EOF
    11660 #line 11660 "configure"
     11676#line 11676 "configure"
    1166111677#include "confdefs.h"
    1166211678
     
    1175811774  lt_status=$lt_dlunknown
    1175911775  cat > conftest.$ac_ext <<EOF
    11760 #line 11760 "configure"
     11776#line 11776 "configure"
    1176111777#include "confdefs.h"
    1176211778
     
    1410214118   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    1410314119   -e 's:$: $lt_compiler_flag:'`
    14104    (eval echo "\"\$as_me:14104: $lt_compile\"" >&5)
     14120   (eval echo "\"\$as_me:14120: $lt_compile\"" >&5)
    1410514121   (eval "$lt_compile" 2>conftest.err)
    1410614122   ac_status=$?
    1410714123   cat conftest.err >&5
    14108    echo "$as_me:14108: \$? = $ac_status" >&5
     14124   echo "$as_me:14124: \$? = $ac_status" >&5
    1410914125   if (exit $ac_status) && test -s "$ac_outfile"; then
    1411014126     # The compiler can only warn and ignore the option if not recognized
     
    1420614222   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    1420714223   -e 's:$: $lt_compiler_flag:'`
    14208    (eval echo "\"\$as_me:14208: $lt_compile\"" >&5)
     14224   (eval echo "\"\$as_me:14224: $lt_compile\"" >&5)
    1420914225   (eval "$lt_compile" 2>out/conftest.err)
    1421014226   ac_status=$?
    1421114227   cat out/conftest.err >&5
    14212    echo "$as_me:14212: \$? = $ac_status" >&5
     14228   echo "$as_me:14228: \$? = $ac_status" >&5
    1421314229   if (exit $ac_status) && test -s out/conftest2.$ac_objext
    1421414230   then
     
    1577615792   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    1577715793   -e 's:$: $lt_compiler_flag:'`
    15778    (eval echo "\"\$as_me:15778: $lt_compile\"" >&5)
     15794   (eval echo "\"\$as_me:15794: $lt_compile\"" >&5)
    1577915795   (eval "$lt_compile" 2>conftest.err)
    1578015796   ac_status=$?
    1578115797   cat conftest.err >&5
    15782    echo "$as_me:15782: \$? = $ac_status" >&5
     15798   echo "$as_me:15798: \$? = $ac_status" >&5
    1578315799   if (exit $ac_status) && test -s "$ac_outfile"; then
    1578415800     # The compiler can only warn and ignore the option if not recognized
     
    1588015896   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    1588115897   -e 's:$: $lt_compiler_flag:'`
    15882    (eval echo "\"\$as_me:15882: $lt_compile\"" >&5)
     15898   (eval echo "\"\$as_me:15898: $lt_compile\"" >&5)
    1588315899   (eval "$lt_compile" 2>out/conftest.err)
    1588415900   ac_status=$?
    1588515901   cat out/conftest.err >&5
    15886    echo "$as_me:15886: \$? = $ac_status" >&5
     15902   echo "$as_me:15902: \$? = $ac_status" >&5
    1588715903   if (exit $ac_status) && test -s out/conftest2.$ac_objext
    1588815904   then
     
    1808718103   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    1808818104   -e 's:$: $lt_compiler_flag:'`
    18089    (eval echo "\"\$as_me:18089: $lt_compile\"" >&5)
     18105   (eval echo "\"\$as_me:18105: $lt_compile\"" >&5)
    1809018106   (eval "$lt_compile" 2>conftest.err)
    1809118107   ac_status=$?
    1809218108   cat conftest.err >&5
    18093    echo "$as_me:18093: \$? = $ac_status" >&5
     18109   echo "$as_me:18109: \$? = $ac_status" >&5
    1809418110   if (exit $ac_status) && test -s "$ac_outfile"; then
    1809518111     # The compiler can only warn and ignore the option if not recognized
     
    1835518371   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    1835618372   -e 's:$: $lt_compiler_flag:'`
    18357    (eval echo "\"\$as_me:18357: $lt_compile\"" >&5)
     18373   (eval echo "\"\$as_me:18373: $lt_compile\"" >&5)
    1835818374   (eval "$lt_compile" 2>conftest.err)
    1835918375   ac_status=$?
    1836018376   cat conftest.err >&5
    18361    echo "$as_me:18361: \$? = $ac_status" >&5
     18377   echo "$as_me:18377: \$? = $ac_status" >&5
    1836218378   if (exit $ac_status) && test -s "$ac_outfile"; then
    1836318379     # The compiler can only warn and ignore the option if not recognized
     
    1845918475   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    1846018476   -e 's:$: $lt_compiler_flag:'`
    18461    (eval echo "\"\$as_me:18461: $lt_compile\"" >&5)
     18477   (eval echo "\"\$as_me:18477: $lt_compile\"" >&5)
    1846218478   (eval "$lt_compile" 2>out/conftest.err)
    1846318479   ac_status=$?
    1846418480   cat out/conftest.err >&5
    18465    echo "$as_me:18465: \$? = $ac_status" >&5
     18481   echo "$as_me:18481: \$? = $ac_status" >&5
    1846618482   if (exit $ac_status) && test -s out/conftest2.$ac_objext
    1846718483   then
     
    2119921215    *-cygwin* | *-mingw*)
    2120021216       case "$CXX" in
    21201     cl* | */cl* | CL* | */CL*)
     21217    cl* | */cl* | CL* | */CL* | icl* | */icl* | ICL* | */ICL*)
    2120221218      { echo "$as_me:$LINENO: Applying patches to libtool for cl compiler" >&5
    2120321219echo "$as_me: Applying patches to libtool for cl compiler" >&6;}
     
    2186621882# Determine the name of the ASL library
    2186721883case "$CC" in
    21868   cl* | */cl* | CL* | */CL*)
     21884  cl* | */cl* | CL* | */CL* | icl* | */icl* | ICL* | */ICL*)
    2186921885    ampllib=amplsolv.lib ;;
    2187021886  *)
     
    2208022096fi
    2208122097case "$CC" in
    22082   cl* | */cl* | CL* | */CL*)
     22098  cl* | */cl* | CL* | */CL* | icl* | */icl* | ICL* | */ICL*)
    2208322099    coin_link_input_cmd=cp ;;
    2208422100esac
     
    2209822114
    2209922115  case "$CC" in
    22100     cl* | */cl* | CL* | */CL*)
     22116    cl* | */cl* | CL* | */CL* | icl* | */icl* | ICL* | */ICL*)
    2210122117         LIBEXT=lib ;;
    2210222118      *) LIBEXT=a ;;
  • branches/Couenne/Couenne/src/Makefile.am

    r999 r1017  
    118118        bound_tightening/tightenBounds.cpp \
    119119        bound_tightening/impliedBounds.cpp \
     120        bound_tightening/operators/impliedBounds-mul.cpp \
     121        bound_tightening/operators/impliedBounds-sum.cpp \
    120122        bound_tightening/operators/impliedBounds-exprSum.cpp \
    121123        bound_tightening/operators/impliedBounds-exprDiv.cpp \
    122124        bound_tightening/operators/impliedBounds-exprMul.cpp \
    123125        bound_tightening/operators/impliedBounds-exprQuad.cpp \
     126        bound_tightening/operators/impliedBounds-exprPow.cpp \
    124127        util/drawCuts.cpp \
    125128        util/rootQ.cpp
     
    141144                         readnl/nl2e.cpp \
    142145                         readnl/invmap.cpp \
    143                          problem/readASLfg.cpp
     146                         readnl/readASLfg.cpp
    144147endif
    145148
     
    156159        -I`$(CYGPATH_W) $(COUENNESRCDIR)/src/util` \
    157160        -I`$(CYGPATH_W) $(COUENNESRCDIR)/src/convex` \
    158         -I`$(CYGPATH_W) $(COUENNESRCDIR)/src/convex/sdpcuts` \
     161        -I`$(CYGPATH_W) $(COUENNESRCDIR)/src/standardize` \
    159162        -I`$(CYGPATH_W) $(COUENNESRCDIR)/src/problem` \
    160163        -I`$(CYGPATH_W) $(COUENNESRCDIR)/src/problem/depGraph` \
     
    205208        branch/projections.hpp \
    206209        branch/CouObjStats.hpp \
     210        standardize/lqelems.hpp \
    207211        expression/operators/exprAbs.hpp \
    208212        expression/operators/exprExp.hpp \
  • branches/Couenne/Couenne/src/Makefile.in

    r999 r1017  
    4949@COIN_HAS_ASL_TRUE@                      readnl/nl2e.cpp \
    5050@COIN_HAS_ASL_TRUE@                      readnl/invmap.cpp \
    51 @COIN_HAS_ASL_TRUE@                      problem/readASLfg.cpp
     51@COIN_HAS_ASL_TRUE@                      readnl/readASLfg.cpp
    5252
    5353
     
    150150        bound_tightening/tightenBounds.cpp \
    151151        bound_tightening/impliedBounds.cpp \
     152        bound_tightening/operators/impliedBounds-mul.cpp \
     153        bound_tightening/operators/impliedBounds-sum.cpp \
    152154        bound_tightening/operators/impliedBounds-exprSum.cpp \
    153155        bound_tightening/operators/impliedBounds-exprDiv.cpp \
    154156        bound_tightening/operators/impliedBounds-exprMul.cpp \
    155157        bound_tightening/operators/impliedBounds-exprQuad.cpp \
     158        bound_tightening/operators/impliedBounds-exprPow.cpp \
    156159        util/drawCuts.cpp util/rootQ.cpp readnl/readnl.cpp \
    157         readnl/nl2e.cpp readnl/invmap.cpp problem/readASLfg.cpp
     160        readnl/nl2e.cpp readnl/invmap.cpp readnl/readASLfg.cpp
    158161@COIN_HAS_ASL_TRUE@am__objects_1 = readnl.lo nl2e.lo invmap.lo \
    159162@COIN_HAS_ASL_TRUE@     readASLfg.lo
     
    188191        boundTightening.lo aggressiveBT.lo fake_tightening.lo obbt.lo \
    189192        obbt_iter.lo tightenBounds.lo impliedBounds.lo \
     193        impliedBounds-mul.lo impliedBounds-sum.lo \
    190194        impliedBounds-exprSum.lo impliedBounds-exprDiv.lo \
    191         impliedBounds-exprMul.lo impliedBounds-exprQuad.lo drawCuts.lo \
    192         rootQ.lo $(am__objects_1)
     195        impliedBounds-exprMul.lo impliedBounds-exprQuad.lo \
     196        impliedBounds-exprPow.lo drawCuts.lo rootQ.lo $(am__objects_1)
    193197libCouenne_la_OBJECTS = $(am_libCouenne_la_OBJECTS)
    194198depcomp = $(SHELL) $(top_srcdir)/../depcomp
     
    467471        bound_tightening/tightenBounds.cpp \
    468472        bound_tightening/impliedBounds.cpp \
     473        bound_tightening/operators/impliedBounds-mul.cpp \
     474        bound_tightening/operators/impliedBounds-sum.cpp \
    469475        bound_tightening/operators/impliedBounds-exprSum.cpp \
    470476        bound_tightening/operators/impliedBounds-exprDiv.cpp \
    471477        bound_tightening/operators/impliedBounds-exprMul.cpp \
    472478        bound_tightening/operators/impliedBounds-exprQuad.cpp \
     479        bound_tightening/operators/impliedBounds-exprPow.cpp \
    473480        util/drawCuts.cpp util/rootQ.cpp $(am__append_1)
    474481COINLIBS = \
     
    496503        $(COUENNESRCDIR)/src/util` -I`$(CYGPATH_W) \
    497504        $(COUENNESRCDIR)/src/convex` -I`$(CYGPATH_W) \
    498         $(COUENNESRCDIR)/src/convex/sdpcuts` -I`$(CYGPATH_W) \
     505        $(COUENNESRCDIR)/src/standardize` -I`$(CYGPATH_W) \
    499506        $(COUENNESRCDIR)/src/problem` -I`$(CYGPATH_W) \
    500507        $(COUENNESRCDIR)/src/problem/depGraph` -I`$(CYGPATH_W) \
     
    535542        branch/projections.hpp \
    536543        branch/CouObjStats.hpp \
     544        standardize/lqelems.hpp \
    537545        expression/operators/exprAbs.hpp \
    538546        expression/operators/exprExp.hpp \
     
    733741@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/impliedBounds-exprDiv.Plo@am__quote@
    734742@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/impliedBounds-exprMul.Plo@am__quote@
     743@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/impliedBounds-exprPow.Plo@am__quote@
    735744@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/impliedBounds-exprQuad.Plo@am__quote@
    736745@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/impliedBounds-exprSum.Plo@am__quote@
     746@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/impliedBounds-mul.Plo@am__quote@
     747@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/impliedBounds-sum.Plo@am__quote@
    737748@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/impliedBounds.Plo@am__quote@
    738749@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/infeasibility.Plo@am__quote@
     
    14771488@am__fastdepCXX_FALSE@  $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o impliedBounds.lo `test -f 'bound_tightening/impliedBounds.cpp' || echo '$(srcdir)/'`bound_tightening/impliedBounds.cpp
    14781489
     1490impliedBounds-mul.lo: bound_tightening/operators/impliedBounds-mul.cpp
     1491@am__fastdepCXX_TRUE@   if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT impliedBounds-mul.lo -MD -MP -MF "$(DEPDIR)/impliedBounds-mul.Tpo" -c -o impliedBounds-mul.lo `test -f 'bound_tightening/operators/impliedBounds-mul.cpp' || echo '$(srcdir)/'`bound_tightening/operators/impliedBounds-mul.cpp; \
     1492@am__fastdepCXX_TRUE@   then mv -f "$(DEPDIR)/impliedBounds-mul.Tpo" "$(DEPDIR)/impliedBounds-mul.Plo"; else rm -f "$(DEPDIR)/impliedBounds-mul.Tpo"; exit 1; fi
     1493@AMDEP_TRUE@@am__fastdepCXX_FALSE@      source='bound_tightening/operators/impliedBounds-mul.cpp' object='impliedBounds-mul.lo' libtool=yes @AMDEPBACKSLASH@
     1494@AMDEP_TRUE@@am__fastdepCXX_FALSE@      DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
     1495@am__fastdepCXX_FALSE@  $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o impliedBounds-mul.lo `test -f 'bound_tightening/operators/impliedBounds-mul.cpp' || echo '$(srcdir)/'`bound_tightening/operators/impliedBounds-mul.cpp
     1496
     1497impliedBounds-sum.lo: bound_tightening/operators/impliedBounds-sum.cpp
     1498@am__fastdepCXX_TRUE@   if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT impliedBounds-sum.lo -MD -MP -MF "$(DEPDIR)/impliedBounds-sum.Tpo" -c -o impliedBounds-sum.lo `test -f 'bound_tightening/operators/impliedBounds-sum.cpp' || echo '$(srcdir)/'`bound_tightening/operators/impliedBounds-sum.cpp; \
     1499@am__fastdepCXX_TRUE@   then mv -f "$(DEPDIR)/impliedBounds-sum.Tpo" "$(DEPDIR)/impliedBounds-sum.Plo"; else rm -f "$(DEPDIR)/impliedBounds-sum.Tpo"; exit 1; fi
     1500@AMDEP_TRUE@@am__fastdepCXX_FALSE@      source='bound_tightening/operators/impliedBounds-sum.cpp' object='impliedBounds-sum.lo' libtool=yes @AMDEPBACKSLASH@
     1501@AMDEP_TRUE@@am__fastdepCXX_FALSE@      DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
     1502@am__fastdepCXX_FALSE@  $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o impliedBounds-sum.lo `test -f 'bound_tightening/operators/impliedBounds-sum.cpp' || echo '$(srcdir)/'`bound_tightening/operators/impliedBounds-sum.cpp
     1503
    14791504impliedBounds-exprSum.lo: bound_tightening/operators/impliedBounds-exprSum.cpp
    14801505@am__fastdepCXX_TRUE@   if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT impliedBounds-exprSum.lo -MD -MP -MF "$(DEPDIR)/impliedBounds-exprSum.Tpo" -c -o impliedBounds-exprSum.lo `test -f 'bound_tightening/operators/impliedBounds-exprSum.cpp' || echo '$(srcdir)/'`bound_tightening/operators/impliedBounds-exprSum.cpp; \
     
    15051530@am__fastdepCXX_FALSE@  $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o impliedBounds-exprQuad.lo `test -f 'bound_tightening/operators/impliedBounds-exprQuad.cpp' || echo '$(srcdir)/'`bound_tightening/operators/impliedBounds-exprQuad.cpp
    15061531
     1532impliedBounds-exprPow.lo: bound_tightening/operators/impliedBounds-exprPow.cpp
     1533@am__fastdepCXX_TRUE@   if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT impliedBounds-exprPow.lo -MD -MP -MF "$(DEPDIR)/impliedBounds-exprPow.Tpo" -c -o impliedBounds-exprPow.lo `test -f 'bound_tightening/operators/impliedBounds-exprPow.cpp' || echo '$(srcdir)/'`bound_tightening/operators/impliedBounds-exprPow.cpp; \
     1534@am__fastdepCXX_TRUE@   then mv -f "$(DEPDIR)/impliedBounds-exprPow.Tpo" "$(DEPDIR)/impliedBounds-exprPow.Plo"; else rm -f "$(DEPDIR)/impliedBounds-exprPow.Tpo"; exit 1; fi
     1535@AMDEP_TRUE@@am__fastdepCXX_FALSE@      source='bound_tightening/operators/impliedBounds-exprPow.cpp' object='impliedBounds-exprPow.lo' libtool=yes @AMDEPBACKSLASH@
     1536@AMDEP_TRUE@@am__fastdepCXX_FALSE@      DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
     1537@am__fastdepCXX_FALSE@  $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o impliedBounds-exprPow.lo `test -f 'bound_tightening/operators/impliedBounds-exprPow.cpp' || echo '$(srcdir)/'`bound_tightening/operators/impliedBounds-exprPow.cpp
     1538
    15071539drawCuts.lo: util/drawCuts.cpp
    15081540@am__fastdepCXX_TRUE@   if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT drawCuts.lo -MD -MP -MF "$(DEPDIR)/drawCuts.Tpo" -c -o drawCuts.lo `test -f 'util/drawCuts.cpp' || echo '$(srcdir)/'`util/drawCuts.cpp; \
     
    15401572@am__fastdepCXX_FALSE@  $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o invmap.lo `test -f 'readnl/invmap.cpp' || echo '$(srcdir)/'`readnl/invmap.cpp
    15411573
    1542 readASLfg.lo: problem/readASLfg.cpp
    1543 @am__fastdepCXX_TRUE@   if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT readASLfg.lo -MD -MP -MF "$(DEPDIR)/readASLfg.Tpo" -c -o readASLfg.lo `test -f 'problem/readASLfg.cpp' || echo '$(srcdir)/'`problem/readASLfg.cpp; \
     1574readASLfg.lo: readnl/readASLfg.cpp
     1575@am__fastdepCXX_TRUE@   if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT readASLfg.lo -MD -MP -MF "$(DEPDIR)/readASLfg.Tpo" -c -o readASLfg.lo `test -f 'readnl/readASLfg.cpp' || echo '$(srcdir)/'`readnl/readASLfg.cpp; \
    15441576@am__fastdepCXX_TRUE@   then mv -f "$(DEPDIR)/readASLfg.Tpo" "$(DEPDIR)/readASLfg.Plo"; else rm -f "$(DEPDIR)/readASLfg.Tpo"; exit 1; fi
    1545 @AMDEP_TRUE@@am__fastdepCXX_FALSE@      source='problem/readASLfg.cpp' object='readASLfg.lo' libtool=yes @AMDEPBACKSLASH@
    1546 @AMDEP_TRUE@@am__fastdepCXX_FALSE@      DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
    1547 @am__fastdepCXX_FALSE@  $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o readASLfg.lo `test -f 'problem/readASLfg.cpp' || echo '$(srcdir)/'`problem/readASLfg.cpp
     1577@AMDEP_TRUE@@am__fastdepCXX_FALSE@      source='readnl/readASLfg.cpp' object='readASLfg.lo' libtool=yes @AMDEPBACKSLASH@
     1578@AMDEP_TRUE@@am__fastdepCXX_FALSE@      DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
     1579@am__fastdepCXX_FALSE@  $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o readASLfg.lo `test -f 'readnl/readASLfg.cpp' || echo '$(srcdir)/'`readnl/readASLfg.cpp
    15481580
    15491581mostlyclean-libtool:
     
    16231655
    16241656distdir: $(DISTFILES)
    1625         $(mkdir_p) $(distdir)/../inc $(distdir)/branch $(distdir)/convex $(distdir)/expression $(distdir)/expression/operators $(distdir)/expression/operators/bounds $(distdir)/problem $(distdir)/problem/depGraph $(distdir)/util
     1657        $(mkdir_p) $(distdir)/../inc $(distdir)/branch $(distdir)/convex $(distdir)/expression $(distdir)/expression/operators $(distdir)/expression/operators/bounds $(distdir)/problem $(distdir)/problem/depGraph $(distdir)/standardize $(distdir)/util
    16261658        @srcdirstrip=`echo "$(srcdir)" | sed 's|.|.|g'`; \
    16271659        topsrcdirstrip=`echo "$(top_srcdir)" | sed 's|.|.|g'`; \
  • branches/Couenne/Couenne/src/bound_tightening/boundTightening.cpp

    r1003 r1017  
    5050      // set infeasibility through a cut 1 <= x0 <= -1
    5151      Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING,
    52                           "#### infeasible node at BT\n");
     52                          "infeasible node at BT\n");
    5353      return false;
    5454    }
  • branches/Couenne/Couenne/src/bound_tightening/fake_tightening.cpp

    r1003 r1017  
    170170#endif
    171171
    172       // apply bound
    173       if (direction) {oub[index] = ub_[index] = fb; chg_bds [index].setUpper (t_chg_bounds::CHANGED);}
    174       else           {olb[index] = lb_[index] = fb; chg_bds [index].setLower (t_chg_bounds::CHANGED);}
     172      /*bool do_not_tighten = false;
     173
     174      // check if cut best known solution
     175      if (optimum_) {
     176        if (direction) {
     177          if ((oub [index] > optimum_ [index]) &&
     178              (fb          < optimum_ [index])) {
     179            printf ("aggressive bt cuts optimum ub %d: %g < %g < %g\n",
     180                    index, fb, optimum_ [index], oub [index]);
     181            do_not_tighten = true;
     182          }
     183        } else {
     184          if ((olb [index] < optimum_ [index]) &&
     185              (fb          > optimum_ [index])) {
     186            printf ("aggressive bt cuts optimum lb %d: %g < %g < %g\n",
     187                    index, olb [index], optimum_ [index], fb);
     188            do_not_tighten = true;
     189          }
     190        }
     191        }*/
     192
     193      //if (!do_not_tighten) {
     194
     195        // apply bound
     196      if (direction) {oub[index]=ub_[index] = fb; chg_bds [index].setUpper (t_chg_bounds::CHANGED);}
     197      else           {olb[index]=lb_[index] = fb; chg_bds [index].setLower (t_chg_bounds::CHANGED);}
    175198
    176199      outer = fb; // we have at least a tightened bound, save it
    177200
    178201      tightened = true;
     202        //}
    179203
    180204      // restore initial bound
  • branches/Couenne/Couenne/src/bound_tightening/impliedBounds.cpp

    r1003 r1017  
    1616  int nchg = 0; //< number of bounds changed for propagation
    1717
    18   if (Jnlst()->ProduceOutput(Ipopt::J_VECTOR, J_BOUNDTIGHTENING)) { 
     18  /*if (Jnlst()->ProduceOutput(Ipopt::J_VECTOR, J_BOUNDTIGHTENING)) { 
    1919    Jnlst()->Printf(Ipopt::J_VECTOR, J_BOUNDTIGHTENING,"=====================implied\n");
    2020    int j=0;
     
    2828      }
    2929    if (j % 6) Jnlst()->Printf(Ipopt::J_VECTOR, J_BOUNDTIGHTENING,"\n");
    30   }
     30    }*/
    3131
    3232  for (int ii = nVars (); ii--;) {
     
    3838      if (lb_ [i] > ub_ [i] + COUENNE_EPS) {
    3939        Jnlst()->Printf(Ipopt::J_DETAILED, J_BOUNDTIGHTENING,
    40                         "#### w_%d has infeasible bounds [%g,%g]\n",
     40                        "implied bounds: w_%d has infeasible bounds [%g,%g]\n",
    4141                        i, lb_ [i], ub_ [i]);
    4242        return -1;
     
    6666      // variables, have changed. If not, skip
    6767
    68       CouNumber
     68      /*CouNumber
    6969        l0 = lb_ [i],
    70         u0 = ub_ [i];
     70        u0 = ub_ [i];*/
    7171
    7272      if (variables_ [i] -> Image () -> impliedBound
    7373          (variables_ [i] -> Index (), lb_, ub_, chg_bds)) {
    7474
    75         if (Jnlst()->ProduceOutput(Ipopt::J_VECTOR, J_BOUNDTIGHTENING)) {
     75        /*if (Jnlst()->ProduceOutput(Ipopt::J_VECTOR, J_BOUNDTIGHTENING)) {
    7676          // todo: send all output through journalist
    7777          Jnlst()->Printf(Ipopt::J_VECTOR, J_BOUNDTIGHTENING,
     
    8383          variables_ [i] -> Image () -> print (std::cout);
    8484          Jnlst()->Printf(Ipopt::J_VECTOR, J_BOUNDTIGHTENING,"\n");
    85         }
     85          }*/
    8686
    8787        /*for (int i=0; i<nVars (); i++)
  • branches/Couenne/Couenne/src/bound_tightening/obbt.cpp

    r999 r1017  
    9191  int nimprov = 0, ni;
    9292 
    93   Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING, ":::::: OBBT on originals ----------------\n");
     93  //  Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING, ":::::: OBBT on originals ----------------\n");
    9494
    9595  if ((ni = call_iter (csi, chg_bds, warmstart, babInfo, objcoe, VAR,  1)) < 0) {
     
    109109  nimprov += ni;
    110110
    111   Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING, ":::::: OBBT on auxiliaries --------------\n");
     111  //  Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING, ":::::: OBBT on auxiliaries --------------\n");
    112112
    113113  if ((ni = call_iter (csi, chg_bds, warmstart, babInfo, objcoe, AUX,  1)) < 0) {
     
    125125  }
    126126
    127   Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING, ":::::: ---------------------------------\n");
     127  //  Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING, ":::::: ---------------------------------\n");
    128128
    129129  free (objcoe);
  • branches/Couenne/Couenne/src/bound_tightening/obbt_iter.cpp

    r1003 r1017  
    6767  // non-decreasing or non-increasing
    6868
    69   //#ifdef DEBUG
    70   static int iter = 0;
    71   //#endif
     69  //  static int iter = 0;
    7270
    7371  std::set <int> deplist;
     
    7573
    7674  bool issimple = false;
    77 
    78   //  CouenneProblem *p = Problem ();
    7975
    8076  exprVar *var = Var (index);
     
    8783
    8884  if ((var -> Type  () == AUX) &&
    89       ((deplistsize = var -> Image () -> DepList (deplist, STOP_AT_AUX, this)) <= 1)) {
     85      ((deplistsize = var -> Image () -> DepList (deplist, STOP_AT_AUX)) <= 1)) {
    9086
    9187    if (!deplistsize) { // funny, the expression is constant...
     
    173169    // m{in,ax}imize xi on csi
    174170
     171    /*
    175172    if (Jnlst()->ProduceOutput(J_MOREVECTOR, J_BOUNDTIGHTENING)) {
    176173      Jnlst()->Printf(J_MOREVECTOR, J_BOUNDTIGHTENING,
     
    181178        char fname [20];
    182179        sprintf (fname, "m%s_w%03d_%03d", (sense == 1) ? "in" : "ax", index, iter);
    183         Jnlst()->Printf(J_MOREVECTOR, J_BOUNDTIGHTENING,"writing %s\n", fname);
     180        //Jnlst()->Printf(J_MOREVECTOR, J_BOUNDTIGHTENING,"writing %s\n", fname);
    184181        csi -> writeLp (fname);
    185182      }
    186183    }
     184    */
    187185
    188186    csi -> setWarmStart (warmstart);
     
    274272      if (doFBBT_ && !(boundTightening (chg_bds, babInfo))) {
    275273        Jnlst()->Printf(J_VECTOR, J_BOUNDTIGHTENING,
    276                         "##### infeasible, bound tightening after OBBT\n");
     274                        "node is infeasible after post-OBBT tightening\n");
    277275        return -1; // tell caller this is infeasible
    278276      }
     
    280278      nImprov++;
    281279    }
    282     Jnlst()->Printf(J_VECTOR, J_BOUNDTIGHTENING,"\n");
     280
     281    //    Jnlst()->Printf(J_VECTOR, J_BOUNDTIGHTENING,"\n");
    283282
    284283    // if we solved the problem on the objective function's
  • branches/Couenne/Couenne/src/bound_tightening/operators/impliedBounds-exprDiv.cpp

    r1003 r1017  
    1818bool exprDiv::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg) {
    1919
     20  //return false; // !!!
     21
    2022  bool resx, resy = resx = false;
    2123
    22   // deal with the "y is a constant"
     24  // y is a constant
    2325  if (arglist_ [1] -> Type () == CONST) {
    2426
     
    2628
    2729    if (ind < 0) {
    28       printf ("exprDiv::impliedBound: Error, w=c/d constants\n");
     30      printf ("exprDiv::impliedBound: Error, w = c/d constants\n");
    2931      exit (-1);
    3032    }
  • branches/Couenne/Couenne/src/bound_tightening/operators/impliedBounds-exprMul.cpp

    r988 r1017  
    1616
    1717bool exprMul::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg) {
     18
     19  //return false; // !!!
    1820
    1921  bool resL, resU = resL = false;
  • branches/Couenne/Couenne/src/bound_tightening/operators/impliedBounds-exprQuad.cpp

    r895 r1017  
    1616//#define DEBUG
    1717
    18 // Compute bounds of exprQuad taking into account the unbounded variables
    19 void computeQuadFiniteBound (const exprQuad *,
    20                              CouNumber &, CouNumber &,
    21                              CouNumber *, CouNumber *,
    22                              int &, int &);
    23 
    2418
    2519/// implied bound processing for quadratic form upon change in lower-
     
    2822bool exprQuad::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg) {
    2923
    30   // return false;
     24  //return false; // !!!
     25
     26  // three implied bound computations, of increasing complexity. The
     27  // first reduces this sum to a sum of aux and tries, after
     28  // tightening bounds of each aux, to infer new bounds
    3129
    3230#ifdef DEBUG
     
    3836
    3937  // get all variables
    40 
     38  /*
    4139  std::set <int>
    4240    lin (index_,   index_   + nlterms_),
     
    4442    quj (qindexJ_, qindexJ_ + nqterms_),
    4543    indUnion;
     44  */
     45
     46  std::set <int> indexSet;
     47
     48  for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el)
     49    indexSet.insert (el -> first -> Index ());
     50
     51  for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
     52
     53    indexSet.insert (row -> first -> Index ());
     54
     55    for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col)
     56      indexSet.insert (col -> first -> Index ());
     57  }
     58
    4659
    4760  // CAUTION: this relies on the first version of bound expressions
    4861  // for quadratic form, i.e. the sum of bounds of independent terms
    49 
     62  /*
    5063  // find all indices appearing in the expression
    5164  std::set_union (qui .begin (), qui .end (),
     
    5669                  lin      .begin (), lin      .end (),
    5770                  std::inserter (indUnion, indUnion.begin ()));
     71  */
    5872
    5973  CouNumber qMin, qMax;
     
    6276
    6377  // get bounds for nonlinear part of the sum
    64   expression *lb, *ub;
    65   exprSum::getBounds (lb, ub);
    66 
    67   qMin = (*lb) ();
    68   qMax = (*ub) ();
    69 
    70   delete lb;
    71   delete ub;
     78  expression *qlb, *qub;
     79  exprSum::getBounds (qlb, qub);
     80  qMin = (*qlb) (); delete qlb;
     81  qMax = (*qub) (); delete qub;
     82  //exprSum::getBounds (qMin, qMax);
    7283
    7384  if (qMin < -COUENNE_INFINITY) indInfLo = -2;
     
    8091#ifdef DEBUG
    8192  printf ("1st phase... inf=(%d,%d) q=[%g,%g].\n", indInfLo, indInfUp, qMin, qMax);
    82   for (std::set <int>:: iterator iter = indUnion.begin ();
    83        iter != indUnion.end (); ++iter)
     93  for (std::set <int>:: iterator iter = indexSet.begin ();
     94       iter != indexSet.end (); ++iter)
    8495    printf ("%4d [%+6g %+6g]\n", *iter, l [*iter], u [*iter]);
    8596#endif
     
    8798  // compute bound on expression using only finite variable bounds
    8899
    89   computeQuadFiniteBound (this, qMin, qMax, l, u, indInfLo, indInfUp);
     100  computeQuadFiniteBound (qMin, qMax, l, u, indInfLo, indInfUp);
    90101
    91102  // similar to impliedBounds-exprGroup. indInf* are -1 if no
     
    111122  // prepare data structure for scanning all variables
    112123  int
    113     minindex = *(indUnion.begin  ()), // minimum index
    114     maxindex = *(indUnion.rbegin ()), // maximum index
     124    minindex = *(indexSet.begin  ()), // minimum index
     125    maxindex = *(indexSet.rbegin ()), // maximum index
    115126    nvars = maxindex - minindex + 1;  // upper bound on # variables involved
    116127
     
    130141
    131142  // assume all coefficients are finite
    132   for (int i=0; i<nlterms_; i++) {
    133 
    134     int ind = index_ [i];
     143  for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el) {
     144    //  for (int i=0; i<nlterms_; i++) {
     145
     146    int ind = el -> first -> Index ();//index_ [i];
    135147
    136148    CouNumber
    137       coe = coeff_ [i],
     149      coe = el -> second, //coeff_ [i],
    138150      li  = l [ind],
    139151      ui  = u [ind];
     
    155167#ifdef DEBUG
    156168  printf ("linear filling (%d,%d): -----------------------\n", minindex, maxindex);
    157   for (std::set <int>:: iterator iter = indUnion.begin ();
    158        iter != indUnion.end (); ++iter)
     169  for (std::set <int>:: iterator iter = indexSet.begin ();
     170       iter != indexSet.end (); ++iter)
    159171    printf ("%4d [%+6g %+6g] [%+6g %+6g]\n", *iter,
    160172            linCoeMin [*iter - minindex], linCoeMax [*iter - minindex],
     
    163175
    164176  // fill in remaining linear coefficients and quadratic ones
    165   for (int i=0; i<nqterms_; i++) {
    166 
    167     int qi = qindexI_ [i],
    168         qj = qindexJ_ [i];
    169 
    170     CouNumber coe = qcoeff_ [i],
    171       li = l [qi], lj = l [qj],
    172       ui = u [qi], uj = u [qj];
    173 
    174     if (qi == qj) { // quadratic term
    175 
    176       qi -= minindex;
    177 
    178       qii [qi] = coe; // quadratic term
    179 
    180       CouNumber
    181         maxbUb = CoinMax (fabs (li), fabs (ui)),
    182         maxbLb = (li >= 0) ? (li) : (ui <= 0) ? (ui) : 0;
    183 
    184       if (maxbUb > COUENNE_INFINITY) maxbUb = 0;
    185 
    186       maxbUb *= maxbUb * coe;
    187       maxbLb *= maxbLb * coe;
    188 
    189       if (coe > 0) {
    190         bCutUb [qi] += maxbUb;
    191         bCutLb [qi] += maxbLb;
    192       } else {
    193         bCutUb [qi] += maxbLb;
    194         bCutLb [qi] += maxbUb;
    195       }
    196     } else { // product term
    197 
    198       coe *= 2;
    199 
    200       CouNumber *b1, *b2;
    201 
    202       if (coe > 0) {b1 = l; b2 = u;}
    203       else         {b1 = u; b2 = l;}
    204 
    205       b1 += minindex;
    206       b2 += minindex;
    207 
    208       qi -= minindex;
    209       qj -= minindex;
    210 
    211       linCoeMin [qi] += coe * b1 [qj];
    212       linCoeMin [qj] += coe * b1 [qi];
    213 
    214       linCoeMax [qi] += coe * b2 [qj];
    215       linCoeMax [qj] += coe * b2 [qi];
    216 
    217       CouNumber
    218         addLo = CoinMin (CoinMin (li*lj, ui*uj),
    219                          CoinMin (ui*lj, li*uj)),
    220         addUp = CoinMax (CoinMax (li*lj, ui*uj),
    221                          CoinMax (ui*lj, li*uj));
    222 
    223       if (addLo < -COUENNE_INFINITY) addLo = 0;
    224       if (addUp >  COUENNE_INFINITY) addUp = 0;
    225 
    226       addLo *= coe;
    227       addUp *= coe;
    228 
    229       if (coe > 0) {
    230         bCutLb [qi] += addLo; bCutUb [qi] += addUp;
    231         bCutLb [qj] += addLo; bCutUb [qj] += addUp;
    232       } else {
    233         bCutLb [qi] += addUp; bCutUb [qi] += addLo;
    234         bCutLb [qj] += addUp; bCutUb [qj] += addLo;
     177  for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
     178
     179    for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
     180
     181      int
     182        qi = row -> first -> Index (),
     183        qj = col -> first -> Index (); //qindexJ_ [i];
     184
     185      CouNumber coe = col -> second, //qcoeff_ [i],
     186        li = l [qi], lj = l [qj],
     187        ui = u [qi], uj = u [qj];
     188
     189      if (qi == qj) { // quadratic term
     190
     191        qi -= minindex;
     192
     193        qii [qi] = coe; // quadratic term
     194
     195        CouNumber
     196          maxbUb = CoinMax (fabs (li), fabs (ui)),
     197          maxbLb = (li >= 0) ? (li) : (ui <= 0) ? (ui) : 0;
     198
     199        if (maxbUb > COUENNE_INFINITY) maxbUb = 0;
     200
     201        maxbUb *= maxbUb * coe;
     202        maxbLb *= maxbLb * coe;
     203
     204        if (coe > 0) {
     205          bCutUb [qi] += maxbUb;
     206          bCutLb [qi] += maxbLb;
     207        } else {
     208          bCutUb [qi] += maxbLb;
     209          bCutLb [qi] += maxbUb;
     210        }
     211      } else { // product term
     212
     213        coe *= 2;
     214
     215        CouNumber *b1, *b2;
     216
     217        if (coe > 0) {b1 = l; b2 = u;}
     218        else         {b1 = u; b2 = l;}
     219
     220        b1 += minindex;
     221        b2 += minindex;
     222
     223        qi -= minindex;
     224        qj -= minindex;
     225
     226        linCoeMin [qi] += coe * b1 [qj];
     227        linCoeMin [qj] += coe * b1 [qi];
     228
     229        linCoeMax [qi] += coe * b2 [qj];
     230        linCoeMax [qj] += coe * b2 [qi];
     231
     232        CouNumber
     233          addLo = CoinMin (CoinMin (li*lj, ui*uj),
     234                           CoinMin (ui*lj, li*uj)),
     235          addUp = CoinMax (CoinMax (li*lj, ui*uj),
     236                           CoinMax (ui*lj, li*uj));
     237
     238        if (addLo < -COUENNE_INFINITY) addLo = 0;
     239        if (addUp >  COUENNE_INFINITY) addUp = 0;
     240
     241        addLo *= coe;
     242        addUp *= coe;
     243
     244        if (coe > 0) {
     245          bCutLb [qi] += addLo; bCutUb [qi] += addUp;
     246          bCutLb [qj] += addLo; bCutUb [qj] += addUp;
     247        } else {
     248          bCutLb [qi] += addUp; bCutUb [qi] += addLo;
     249          bCutLb [qj] += addUp; bCutUb [qj] += addLo;
     250        }
    235251      }
    236252    }
     
    239255#ifdef DEBUG
    240256  printf ("quad filling: -----------------------\n");
    241   for (std::set <int>:: iterator iter = indUnion.begin ();
    242        iter != indUnion.end (); ++iter)
     257  for (std::set <int>:: iterator iter = indexSet.begin ();
     258       iter != indexSet.end (); ++iter)
    243259    printf ("%4d [%+6g %+6g] [%+6g %+6g]\n", *iter,
    244260            linCoeMin [*iter - minindex], linCoeMax [*iter - minindex],
     
    247263
    248264  // Done filling vectors /////////////////////////////////////////////////////////////
    249   // Now improve each independent variable in the set indUnion
     265  // Now improve each independent variable in the set indexSet
    250266
    251267  bool updated = false;
    252268
    253   for (std::set <int>:: iterator iter = indUnion.begin ();
    254        iter != indUnion.end (); ++iter) {
     269  for (std::set <int>:: iterator iter = indexSet.begin ();
     270       iter != indexSet.end (); ++iter) {
    255271
    256272    int ind  = *iter,
  • branches/Couenne/Couenne/src/bound_tightening/operators/impliedBounds-exprSum.cpp

    r983 r1017  
    88 */
    99
     10#include "CoinHelperFunctions.hpp"
    1011#include "exprSum.hpp"
    1112#include "exprGroup.hpp"
    1213
    13 
    1414/// vector operation to find bound to variable in a sum
    15 static long double scanBounds (int, int, int *, CouNumber *, CouNumber *, int *);
     15static CouNumber scanBounds (int, int, int *, CouNumber *, CouNumber *, int *);
    1616
    1717
     
    2020
    2121bool exprSum::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg) {
     22
     23  //return false; // !!!
    2224
    2325  /**
     
    4749  // compute number of pos/neg coefficients and sum up constants
    4850
    49   int nterms = nargs_; // # nonconstant terms
    50   int nlin   = 0;      // # linear terms
    51 
    52   long double a0 = 0.;   // constant term in the sum
    53 
    54   if (code () == COU_EXPRGROUP) { // if exprGroup, compute no. linear terms
    55 
    56     a0 +=      dynamic_cast <exprGroup *> (this) -> getc0      ();
    57     int *ind = dynamic_cast <exprGroup *> (this) -> getIndices ();
    58 
    59     while (*ind++ >= 0) // count linear terms
    60       nlin++;
     51  int
     52    nterms = nargs_, // # nonconstant terms
     53    nlin   = 0;      // # linear terms
     54
     55  CouNumber a0 = 0.;   // constant term in the sum
     56
     57  exprGroup *eg = NULL;
     58
     59  bool isExprGroup = (code () == COU_EXPRGROUP);
     60
     61  if (isExprGroup) { // if exprGroup, compute no. linear terms
     62
     63    eg = dynamic_cast <exprGroup *> (this);
     64
     65    a0   += eg -> getc0 ();
     66    nlin += eg -> lcoeff (). size ();
    6167  }
    6268
     
    7278            *I2 = (int       *) malloc (nlin   * sizeof (int));
    7379
     80  /*CoinFillN (C1, nterms, 0.);
     81  CoinFillN (C2, nlin,   0.);
     82
     83  CoinFillN (I1, nterms, 0);
     84  CoinFillN (I2, nlin,   0);*/
     85
    7486  int ipos, ineg = ipos = 0; // #pos and #neg terms
    7587
     
    7890  for (int i = nargs_; i--;) {
    7991
    80     register int index = arglist_ [i] -> Index ();
    81 
    82     if (index == -1)
     92    int index = arglist_ [i] -> Index ();
     93
     94    if (index < 0) // assume a non variable is a constant
    8395      a0 += arglist_ [i] -> Value ();
    8496    else {
     
    90102  // classify terms as pos/neg or constant for the remaining of exprGroup
    91103
    92   if (code () == COU_EXPRGROUP) {
    93 
    94     int       *ind = dynamic_cast <exprGroup *> (this) -> getIndices ();
    95     CouNumber *coe = dynamic_cast <exprGroup *> (this) -> getCoeffs  ();
    96 
    97     for (;*ind >= 0; ind++, coe++)
    98       if      (*coe >   COUENNE_EPS) {I1 [ipos] = *ind; C1 [ipos++] = *coe;}
    99       else if (*coe < - COUENNE_EPS) {I2 [ineg] = *ind; C2 [ineg++] = *coe;}
     104  if (isExprGroup) {
     105
     106    exprGroup::lincoeff &lcoe = eg -> lcoeff ();
     107
     108    for (exprGroup::lincoeff::iterator el = lcoe.begin ();
     109         el != lcoe.end (); ++el) {
     110
     111      CouNumber coe = el -> second;
     112      int       ind = el -> first -> Index ();
     113
     114      if      (coe >  0.) {I1 [ipos] = ind; C1 [ipos++] = coe;}
     115      else if (coe < -0.) {I2 [ineg] = ind; C2 [ineg++] = coe;}
     116    }
    100117  }
    101118
    102119  // realloc to save some memory
    103120
    104   I1 = (int *) realloc (I1, ipos * sizeof (int));
     121  /*I1 = (int *) realloc (I1, ipos * sizeof (int));
    105122  I2 = (int *) realloc (I2, ineg * sizeof (int));
    106123
    107124  C1 = (CouNumber *) realloc (C1, ipos * sizeof (CouNumber));
    108   C2 = (CouNumber *) realloc (C2, ineg * sizeof (CouNumber));
     125  C2 = (CouNumber *) realloc (C2, ineg * sizeof (CouNumber));*/
     126
     127  /*printf ("  a0 = %g\n", a0);
     128  printf ("  pos: "); for (int i=0; i<ipos; i++) printf ("%g x%d [%g,%g] ", C1 [i], I1 [i], l [I1 [i]], u [I1 [i]]); printf ("\n");
     129  printf ("  neg: "); for (int i=0; i<ineg; i++) printf ("%g x%d [%g,%g] ", C2 [i], I2 [i], l [I2 [i]], u [I2 [i]]); printf ("\n");*/
    109130
    110131  // now we have correct  I1, I2, C1, C2, ipos, ineg, and a0
     
    125146  //             + \sum_{i in I_2: a_i > -\infinity} a_i l_i$
    126147
    127   long double
     148  CouNumber
    128149
    129150    lower = a0 + scanBounds (ipos, -1, I1, C1, l, &infLo1)
     
    132153    upper = a0 + scanBounds (ipos, +1, I1, C1, u, &infUp1)
    133154               + scanBounds (ineg, -1, I2, C2, l, &infLo2);
    134  
     155
    135156  // Now compute lower bound for all or for some of the variables:
    136157  // There is a bound for all variables only if both infUp1 and infLo2
     
    173194        if (updateBound (+1, u + wind, upper)) {
    174195          chg [wind].setUpper(t_chg_bounds::CHANGED);
     196          tighter = true;
    175197        }
    176198
     
    286308/// sum bounds depending on coefficients
    287309
    288 static long double scanBounds (int        num,      /// cardinality of the set (I1 or I2)
    289                                int        sign,     /// +1: check for +inf, -1: check for -inf
    290                                int       *indices,  /// indices in the set, $i\in I_k$
    291                                CouNumber *coeff,    /// coefficients in the set
    292                                CouNumber *bounds,   /// variable bounds to check
    293                                int       *infnum) { /// return value: index of variable with infinite
     310static CouNumber scanBounds (int        num,      /// cardinality of the set (I1 or I2)
     311                             int        sign,     /// +1: check for +inf, -1: check for -inf
     312                             int       *indices,  /// indices in the set, $i\in I_k$
     313                             CouNumber *coeff,    /// coefficients in the set
     314                             CouNumber *bounds,   /// variable bounds to check
     315                             int       *infnum) { /// return value: index of variable with infinite
    294316                                                  /// bound, or -2 if more than one exist
    295317  CouNumber bound = 0.;
     
    301323    // be sensitive here, check for bounds a little within the finite realm
    302324
    303     if (((sign > 0) ? bd : -bd) > COUENNE_INFINITY / 1e10 - 1) {
    304 
    305       bounds [indices [i]] = (sign > 0) ? 1e300 : -1e300;
     325    if (((sign > 0) ? bd : -bd) > COUENNE_INFINITY / 1e10) {
     326
     327      bounds [indices [i]] = (sign > 0) ? COIN_DBL_MAX : -COIN_DBL_MAX;
     328      //bounds [indices [i]] = (sign > 0) ? 1e300 : -1e300;
    306329
    307330      // this variable has an infinite bound, mark it
  • branches/Couenne/Couenne/src/bound_tightening/tightenBounds.cpp

    r1003 r1017  
    5252
    5353    if (lb_ [i] > ub_ [i] + COUENNE_EPS) {
    54       //#ifdef DEBUG
     54
    5555      if (Jnlst()->ProduceOutput(J_VECTOR, J_BOUNDTIGHTENING)) {
    5656        Jnlst()->Printf(J_VECTOR, J_BOUNDTIGHTENING,
     
    6161        Jnlst()->Printf(J_VECTOR, J_BOUNDTIGHTENING,"\n");
    6262      }
    63       //#endif
     63
    6464      return -1; // declare this node infeasible
    6565    }
     
    8585
    8686      if (ll > uu + COUENNE_EPS) {
    87         //#ifdef DEBUG
     87
    8888        if (Jnlst()->ProduceOutput(J_VECTOR, J_BOUNDTIGHTENING)) {
    8989          Jnlst()->Printf(J_VECTOR, J_BOUNDTIGHTENING,
     
    9494          Jnlst()->Printf(J_VECTOR, J_BOUNDTIGHTENING,"\n");
    9595        }
    96         //#endif
     96
    9797        return -1; // declare this node infeasible
    9898      }
    9999
    100       //      bool chg = false;
    101100
    102101      // check if lower bound got higher
     
    107106           (fabs (ll / (lb_ [i]) - 1) > COUENNE_EPS)) ) {
    108107
    109         /*printf ("update lbound %d: %.10e >= %.10e + %.12e\n",
    110           i+j, ll, lb_ [i+j], COUENNE_EPS);*/
     108        /*if (Jnlst()->ProduceOutput(J_VECTOR, J_BOUNDTIGHTENING)) {
    111109
    112         /*printf ("update lbound %d: %g >= %g\n",
    113           i+j, ll, lb_ [i+j]);*/
    114 
    115         if (Jnlst()->ProduceOutput(J_VECTOR, J_BOUNDTIGHTENING)) {
    116           // ToDo: Pipe all output through journalist
    117110          Jnlst()->Printf(J_VECTOR, J_BOUNDTIGHTENING,
    118111                          "propa %2d [%g,(%g)] -> [%g,(%g)] (%g) ",
     
    135128            Jnlst()->Printf(J_VECTOR, J_BOUNDTIGHTENING,"\n");
    136129          }
    137         }
     130          }*/
    138131
    139132        lb_ [i] = ll;
     
    156149
    157150        if (Jnlst()->ProduceOutput(J_VECTOR, J_BOUNDTIGHTENING)) {
    158           Jnlst()->Printf(J_VECTOR, J_BOUNDTIGHTENING,
     151          /*Jnlst()->Printf(J_VECTOR, J_BOUNDTIGHTENING,
    159152                          "prop %2d [(%g),%g] -> [(%g),%g] (%g) ",
    160153                          i, lb_ [i], ub_ [i], ll, uu, ub_ [i] - uu);
     
    162155          Jnlst()->Printf(J_VECTOR, J_BOUNDTIGHTENING," := ");
    163156          Var (i) -> Image () -> print (std::cout);
    164           Jnlst()->Printf(J_VECTOR, J_BOUNDTIGHTENING,"\n");
     157          Jnlst()->Printf(J_VECTOR, J_BOUNDTIGHTENING,"\n");*/
    165158
    166159          if (optimum_ &&
  • branches/Couenne/Couenne/src/branch/CouenneBranchingObject.cpp

    r978 r1017  
    1919CouNumber CouenneBranchingObject::alpha_ = 0.25;
    2020
     21#define AGGR_MUL 2
     22
    2123
    2224/// computes a not-too-bad point where to branch, in the "middle" of an interval
     
    3234
    3335  CouNumber alpha = CouenneBranchingObject::Alpha ();
    34  
    35   if ((l > -COUENNE_INFINITY) && ((x-l) / (u-l) < COUENNE_NEAR_BOUND) ||
    36       (u <  COUENNE_INFINITY) && ((u-x) / (u-l) < COUENNE_NEAR_BOUND))
    37     if      (u <   COUENNE_INFINITY)
    38       if    (l > - COUENNE_INFINITY) x = alpha * x + (1. - alpha) * (l + u) / 2.;
    39       else                           x = -1 + (u<0) ? u*2 : u/2;
    40     else if (l > - COUENNE_INFINITY) x = +1 + (l>0) ? l*2 : l/2;
    41     else                             x = 0;
    42   //  else retval = x;//alpha * x + (1. - alpha) * (l + u) / 2.;
    4336
    44   return x;
     37  if   (l < -COUENNE_INFINITY / 10)
     38    if (u >  COUENNE_INFINITY / 10) return x; // 0.                                    // ]-inf,+inf[
     39    else                            return ((x < -COUENNE_EPS) ? (AGGR_MUL * (-1+x)) : // ]-inf,u]
     40                                            (x >  COUENNE_EPS) ? 0. : -AGGR_MUL);
     41  else
     42    if (u >  COUENNE_INFINITY / 10) return ((x >  COUENNE_EPS) ? (AGGR_MUL *  (1+x)) : // [l,+inf[
     43                                            (x < -COUENNE_EPS) ? 0. :  AGGR_MUL);
     44    else                            return alpha * x + (1. - alpha) * (l + u) / 2.;    // [l,u]
    4545}
    4646
     
    5252*/
    5353
    54 CouenneBranchingObject::CouenneBranchingObject (JnlstPtr jnlst, int index, int way,
    55                                                 CouNumber brpoint, bool isint):
    56 
    57   index_   (index),
    58   integer_ (isint),
    59   jnlst_   (jnlst) {
     54CouenneBranchingObject::CouenneBranchingObject (JnlstPtr jnlst, expression *var,
     55                                                int way, CouNumber brpoint):
     56  variable_ (var),
     57  //  integer_  (isint),
     58  jnlst_    (jnlst) {
    6059
    6160  firstBranch_ =  (way == TWO_LEFT)      ? 0 :
     
    6362                 ((CoinDrand48 () < 0.5) ? 0 : 1));
    6463
    65   assert (index_ >= 0);
     64  //  assert (index_ >= 0);
    6665
    67   CouNumber x = expression::Variable (index_); // current solution
     66  CouNumber x = (*variable_) ();//expression::Variable (index_); // current solution
    6867
    6968  if (fabs (brpoint) < COUENNE_INFINITY)
     
    8887  //  assert (fabs (u-l) > COUENNE_EPS);
    8988
    90   value_ = midInterval (x, expression::Lbound (index_), expression::Ubound (index_));
     89  CouNumber lb, ub;
     90
     91  var -> getBounds (lb, ub);
     92
     93  value_ = midInterval (x, lb, ub);//(*(var -> Lb ())) (), (*(var -> Ub ())) ());
    9194  //  value_ = midInterval (x, expression::Lbound (index_), expression::Ubound (index_));
    9295
     96  //  if (jnlst_ -> ProduceOutput (J_DETAILED, J_BRANCHING)) {
    9397  jnlst_ -> Printf (J_DETAILED, J_BRANCHING,
    94                     "=== x%d will branch on %g (at %g) [%g,%g]\n  with firstBranch_ = %d\n",
    95                     index_, value_,
    96                     expression::Variable (index_),
    97                     expression::Lbound   (index_),
    98                     expression::Ubound   (index_),
     98                    "Branch: x%-3d will branch on %g (at %g) [%g,%g]; firstBranch_ = %d\n",
     99                    variable_ -> Index (),
     100                    value_,
     101                    x, lb, ub,
     102                    //expression::Variable (index_),
     103                    //expression::Lbound   (index_),
     104                    //expression::Ubound   (index_),
    99105                    firstBranch_);
    100106}
     
    102108
    103109/** \brief Execute the actions required to branch, as specified by the
    104            current state of the branching object, and advance the
    105            object's state.
    106            Returns change in guessed objective on next branch
    107 */
     110 *         current state of the branching object, and advance the
     111 *         object's state.
     112 *
     113 *         Returns change in guessed objective on next branch
     114 */
    108115
    109116double CouenneBranchingObject::branch (OsiSolverInterface * solver) {
     
    112119  //       1 if ">=" node
    113120
    114   int way = (!branchIndex_) ? firstBranch_ : !firstBranch_;
     121  int
     122    way   = (!branchIndex_) ? firstBranch_ : !firstBranch_,
     123    index = variable_ -> Index ();
     124
     125  bool integer = variable_ -> isInteger ();
    115126
    116127  CouNumber
    117     l    = solver -> getColLower () [index_],
    118     u    = solver -> getColUpper () [index_],
     128    l    = solver -> getColLower () [index],
     129    u    = solver -> getColUpper () [index],
    119130    brpt = value_;
    120131
    121132  if (jnlst_->ProduceOutput(J_DETAILED, J_BRANCHING)) {
    122133    if (way) {
    123       if      (value_ < l)             jnlst_->Printf(J_DETAILED, J_BRANCHING, "Nonsense up-br: [ %.8f ,(%.8f)] -> %.8f\n", l,u,value_);
    124       else if (value_ < l+COUENNE_EPS) jnlst_->Printf(J_DETAILED, J_BRANCHING, "## WEAK  up-br: [ %.8f ,(%.8f)] -> %.8f\n", l,u,value_);
     134      if      (value_ < l)             
     135        jnlst_->Printf(J_DETAILED, J_BRANCHING,
     136                       "Nonsense up-br: [ %.8f ,(%.8f)] -> %.8f\n", l,u,value_);
     137      else if (value_ < l+COUENNE_EPS)
     138        jnlst_->Printf(J_DETAILED, J_BRANCHING,
     139                       "## WEAK  up-br: [ %.8f ,(%.8f)] -> %.8f\n", l,u,value_);
    125140    } else {
    126       if      (value_ > u)             jnlst_->Printf(J_DETAILED, J_BRANCHING, "Nonsense dn-br: [(%.8f), %.8f ] -> %.8f\n", l,u,value_);
    127       else if (value_ > u+COUENNE_EPS) jnlst_->Printf(J_DETAILED, J_BRANCHING, "## WEAK  dn-br: [(%.8f), %.8f ] -> %.8f\n", l,u,value_);
     141      if      (value_ > u)             
     142        jnlst_->Printf(J_DETAILED, J_BRANCHING,
     143                       "Nonsense dn-br: [(%.8f), %.8f ] -> %.8f\n", l,u,value_);
     144      else if (value_ > u+COUENNE_EPS)
     145        jnlst_->Printf(J_DETAILED, J_BRANCHING,
     146                       "## WEAK  dn-br: [(%.8f), %.8f ] -> %.8f\n", l,u,value_);
    128147    }
    129148  }
     
    132151  if (brpt > u) brpt = u;
    133152
    134   if (!way) solver -> setColUpper (index_, integer_ ? floor (brpt) : brpt); // down branch
    135   else      solver -> setColLower (index_, integer_ ? ceil  (brpt) : brpt); // up   branch
     153  if (!way) solver -> setColUpper (index, integer ? floor (brpt) : brpt); // down branch
     154  else      solver -> setColLower (index, integer ? ceil  (brpt) : brpt); // up   branch
    136155
    137156  // TODO: apply bound tightening to evaluate change in dual bound
    138157
    139   jnlst_->Printf(J_DETAILED, J_BRANCHING, "### Branch: x%d %c= %g\n",
    140           index_, way ? '>' : '<', integer_ ? (way ? ceil : floor) (brpt) : brpt);
     158  jnlst_ -> Printf (J_DETAILED, J_BRANCHING, "Branching: x%-3d %c= %g\n",
     159          index, way ? '>' : '<', integer ? (way ? ceil : floor) (brpt) : brpt);
    141160
    142161  branchIndex_++;
  • branches/Couenne/Couenne/src/branch/CouenneBranchingObject.hpp

    r969 r1017  
    4343
    4444  /// Constructor
    45   CouenneBranchingObject (JnlstPtr jnlst, int, int, CouNumber = - COIN_DBL_MAX, bool = false);
     45  CouenneBranchingObject (JnlstPtr jnlst, expression *, int, CouNumber = - COIN_DBL_MAX);
    4646
    4747  /// Copy constructor
    4848  CouenneBranchingObject (const CouenneBranchingObject &src):
    4949    OsiTwoWayBranchingObject (src),
    50     index_   (src.index_),
    51     integer_ (src.integer_),
    52     jnlst_   (src.jnlst_) {}
     50    variable_ (src.variable_),
     51    //    integer_ (src.integer_),
     52    jnlst_    (src.jnlst_) {}
    5353
    5454  /// Cloning method
     
    6969  /// either x or y, chosen previously with a call to getFixVar()
    7070  /// expression *reference_;
    71   int index_;
     71  expression *variable_;
    7272
    7373  /// True if the related variable is integer
    74   bool integer_;
     74  //  bool integer_;
    7575
    7676  /// Global value for convex combination between current point and
     
    8282};
    8383
    84 /// returns a point "inside enough" a given interval, or x if it is already
    85 CouNumber midInterval (CouNumber x, CouNumber l, CouNumber u);
    86 
    8784#endif
  • branches/Couenne/Couenne/src/branch/CouenneChooseVariable.cpp

    r895 r1017  
    4444
    4545  /// LEAVE THIS HERE. Need update of current point to evaluate infeasibility
    46   problem_ -> update (info -> solution_,
    47                       info -> lower_,
    48                       info -> upper_);
     46  problem_ -> initAuxs (info -> solution_,
     47                        info -> lower_,
     48                        info -> upper_);
     49
     50  problem_ -> update (info -> solution_, NULL, NULL);
     51
     52  /*for (int i = problem_ -> nVars (); i--;) {
     53    info -> lower_ [i] = problem_ -> Lb (i);
     54    info -> upper_ [i] = problem_ -> Ub (i);
     55    }*/
     56
     57  //CoinCopyN ((double *) (problem_ -> Lb ()), problem_ -> nVars (), info -> lower_);
     58  //CoinCopyN ((double *) (problem_ -> Ub ()), problem_ -> nVars (), info -> upper_);
    4959
    5060  // Make it stable, in OsiChooseVariable::setupList() numberObjects must be 0.
  • branches/Couenne/Couenne/src/branch/CouenneChooseVariable.hpp

    r895 r1017  
    1414#include "OsiChooseVariable.hpp"
    1515#include "CouenneProblem.hpp"
    16 //#include "BonBabSetupBase.hpp"
    1716
    1817// Forward declaration
  • branches/Couenne/Couenne/src/branch/CouenneObject.cpp

    r999 r1017  
    1616#include "exprGroup.hpp"
    1717#include "exprQuad.hpp"
    18 
    19 /// global index for CouenneObjects
    20 //int CouObjStats::Index_;
     18#include "lqelems.hpp"
    2119
    2220
     
    2624
    2725  reference_ (ref),
    28   brVarInd_  (-1),
     26  brVar_     (NULL),
    2927  brPts_     (NULL),
    3028  whichWay_  (BRANCH_NONE),
     
    5250CouenneObject::CouenneObject (const CouenneObject &src):
    5351  reference_ (src.reference_),
    54   brVarInd_  (src.brVarInd_),
     52  brVar_     (src.brVar_),
    5553  brPts_     (NULL),
    5654  whichWay_  (src.whichWay_),
     
    127125
    128126    exprGroup *e = dynamic_cast <exprGroup *> (expr);
    129     int *indices = e -> getIndices ();
    130 
    131     for (int n = e -> getnLTerms (); n--; indices++) {
    132       val = info -> solution_ [*indices];
    133       solver -> setColLower (*indices, val-TOL);
    134       solver -> setColUpper (*indices, val+TOL);
     127
     128    const exprGroup::lincoeff &lcoe = e -> lcoeff ();
     129
     130    for (int n = lcoe.size(), i=0; n--; i++) {
     131      //exprGroup::lincoeff::iterator el = lcoe.begin (); el != lcoe.end (); ++el) {
     132      int index = lcoe [i]. first -> Index ();
     133      val = info -> solution_ [index];
     134      solver -> setColLower (index, val-TOL);
     135      solver -> setColUpper (index, val+TOL);
    135136    }
    136137
     
    140141      exprQuad *e = dynamic_cast <exprQuad *> (expr);
    141142
    142       int *qi = e -> getQIndexI (),
    143           *qj = e -> getQIndexJ ();
    144 
    145       for (int n = e -> getnQTerms (); n--; qi++, qj++) {
    146 
    147         val = info -> solution_ [*qi];
    148         solver -> setColLower (*qi, val-TOL);
    149         solver -> setColUpper (*qi, val+TOL);
    150 
    151         val = info -> solution_ [*qj];
    152         solver -> setColLower (*qj, val-TOL);
    153         solver -> setColUpper (*qj, val+TOL);
     143      exprQuad::sparseQ q = e -> getQ ();
     144
     145      for (exprQuad::sparseQ::iterator row = q.begin ();
     146           row != q.end (); ++row) {
     147
     148        int xind = row -> first -> Index ();
     149
     150        val = info -> solution_ [xind];
     151        solver -> setColLower (xind, val-TOL);
     152        solver -> setColUpper (xind, val+TOL);
     153
     154        for (exprQuad::sparseQcol::iterator col = row -> second.begin ();
     155             col != row -> second.end (); ++col) {
     156
     157          int yind = col -> first -> Index ();
     158
     159          val = info -> solution_ [yind];
     160          solver -> setColLower (yind, val-TOL);
     161          solver -> setColUpper (yind, val+TOL);
     162        }
    154163      }
    155164    }
    156165  }
     166
     167  // TODO: better value through one run of btCore ()
    157168
    158169  return 0.;
     
    165176                                                 int way) const {
    166177
    167   bool isint = (brVarInd_ >= 0) && (si -> isInteger (brVarInd_));
     178  bool isint = (brVar_) && brVar_ -> isInteger ();
     179  // (brVarInd_ >= 0) && (si -> isInteger (brVarInd_));
    168180
    169181  // way has suggestion from CouenneObject::infeasibility(), but not
     
    208220
    209221
    210   if (brVarInd_ >= 0) // if applied latest selectBranching
     222  if (brVar_) // if applied latest selectBranching
    211223
    212224    switch (way) {
     
    215227    case TWO_RIGHT:
    216228    case TWO_RAND:
    217       jnlst_->Printf(J_DETAILED, J_BRANCHING, "2way Branch x%d at %g [%d] (%d)\n", brVarInd_, *brPts_, way, isint);
    218       return new CouenneBranchingObject (jnlst_, brVarInd_, way, *brPts_, isint);
     229      jnlst_->Printf(J_DETAILED, J_BRANCHING,
     230                     "2way Branch x%d at %g [%d] (%d)\n",
     231                     brVar_ -> Index (), *brPts_, way, isint);
     232      return new CouenneBranchingObject (jnlst_, brVar_, way, *brPts_);
    219233    case THREE_LEFT:
    220234    case THREE_CENTER:
    221235    case THREE_RIGHT:
    222236    case THREE_RAND:
    223       jnlst_->Printf(J_DETAILED, J_BRANCHING, "3Way Branch x%d @ %g ][ %g [%d] (%d)\n", brVarInd_, *brPts_, brPts_ [1], way, isint);
    224       return new CouenneThreeWayBranchObj (jnlst_, brVarInd_, brPts_ [0], brPts_ [1], way, isint);
     237      jnlst_->Printf(J_DETAILED, J_BRANCHING,
     238                     "3Way Branch x%d @ %g ][ %g [%d] (%d)\n",
     239                     brVar_ -> Index (), *brPts_, brPts_ [1], way, isint);
     240      return new CouenneThreeWayBranchObj (jnlst_, brVar_, brPts_ [0], brPts_ [1], way);
    225241    default:
    226242      printf ("CouenneObject::createBranch(): way=%d has no sense\n", way);
     
    234250    jnlst_->Printf(J_DETAILED, J_BRANCHING, "CO::createBranch: ");
    235251    reference_ -> print (std::cout);                              printf (" = ");
    236     reference_ -> Image () -> print (std::cout);                 printf (" --> branch on ");
     252    reference_ -> Image () -> print (std::cout); fflush (stdout); printf (" --> branch on ");
    237253    reference_ -> Image () -> getFixVar () -> print (std::cout);  printf ("\n");
    238254  }
     
    287303        || (fabs (ur-xr) < COUENNE_EPS)
    288304        || (fabs (ur-lr) < COUENNE_EPS))
    289       return new CouenneBranchingObject (jnlst_, index, way, x, depvar -> isInteger ()); 
    290   }
    291 
    292   return new CouenneBranchingObject (jnlst_, ref_ind, way, xr, reference_ -> isInteger ());
     305      return new CouenneBranchingObject (jnlst_, depvar, way, x); 
     306  }
     307
     308  return new CouenneBranchingObject (jnlst_, reference_, way, xr);
    293309}
  • branches/Couenne/Couenne/src/branch/CouenneObject.hpp

    r969 r1017  
    2020#include "CouenneJournalist.hpp"
    2121
    22 //#include "CouObjStats.hpp"
    23 
    24 //class CouObjStats;
    2522
    2623/// Define what kind of branching (two- or three-way) and where to
    27 /// start from: left, (center,) or right. The last is to help diversify
    28 /// branching through randomization, which may help when the same
    29 /// variable is branched upon in several points of the BB tree.
     24/// start from: left, (center,) or right. The last is to help
     25/// diversify branching through randomization, which may help when the
     26/// same variable is branched upon in several points of the BB tree.
    3027
    3128enum {TWO_LEFT,                 TWO_RIGHT,   TWO_RAND,
    3229      THREE_LEFT, THREE_CENTER, THREE_RIGHT, THREE_RAND, BRANCH_NONE};
    3330
     31//
     32class funtriplet;
     33CouNumber minMaxDelta (funtriplet *ft, CouNumber x, CouNumber y, CouNumber lb, CouNumber ub);
     34CouNumber maxHeight   (funtriplet *ft, CouNumber x, CouNumber y, CouNumber lb, CouNumber ub);
     35
     36/// returns a point "inside enough" a given interval, or x if it is already
     37CouNumber midInterval (CouNumber x, CouNumber l, CouNumber u);
    3438
    3539/// OsiObject for auxiliary variables $w=f(x)$.
     
    8286  {return strategy_;}
    8387
     88  /// pick branching point based on current strategy
     89  CouNumber getBrPoint (funtriplet *ft, CouNumber x0, CouNumber y0, CouNumber l, CouNumber u) const {
     90
     91    switch (strategy_) {
     92
     93    case CouenneObject::MIN_AREA:     return maxHeight   (ft, x0, y0, l, u); break;
     94    case CouenneObject::BALANCED:     return minMaxDelta (ft, x0, y0, l, u); break;
     95    case CouenneObject::MID_INTERVAL:
     96    default:                          return midInterval (    x0,     l, u); break;
     97    }
     98  }
     99
    84100protected:
    85101
    86   /// The (auxiliary) variable which this branching object refers
    87   /// to. If the expression is w=f(x,y), this is w, as opposed to
     102  /// The (auxiliary) variable this branching object refers to. If the
     103  /// expression is w=f(x,y), this is w, as opposed to
    88104  /// CouenneBranchingObject, where it would be either x or y.
    89105  exprVar *reference_;
    90106
    91107  /// index on the branching variable
    92   mutable int brVarInd_;
     108  mutable expression *brVar_;
    93109
    94110  /// where to branch. It is a vector in the event we want to use a
     
    103119  mutable int whichWay_;
    104120
    105   /// statistics on use of this objects
    106   //  mutable CouObjStats *stats_;
    107 
    108121  /// Branching point selection strategy
    109122  enum brSelStrat strategy_;
     
    113126};
    114127
    115 //
    116 class funtriplet;
    117 CouNumber minMaxDelta (funtriplet *ft, CouNumber x, CouNumber y, CouNumber lb, CouNumber ub);
    118 CouNumber maxHeight   (funtriplet *ft, CouNumber x, CouNumber y, CouNumber lb, CouNumber ub);
    119 
    120128#endif
  • branches/Couenne/Couenne/src/branch/CouenneThreeWayBranchObj.cpp

    r969 r1017  
    1717/// a call to operator () of that exprAux.
    1818CouenneThreeWayBranchObj::CouenneThreeWayBranchObj (JnlstPtr jnlst,
    19                                                     int index,
     19                                                    expression *brVar,
    2020                                                    CouNumber lcrop,
    2121                                                    CouNumber rcrop,
    22                                                     int way,
    23                                                     bool isint):
    24   index_   (index),
     22                                                    int way
     23                                                    //bool isint
     24                                                    ):
     25  brVar_   (brVar),
    2526  lcrop_   (lcrop),
    2627  rcrop_   (rcrop),
    27   integer_ (isint),
    2828  jnlst_   (jnlst) {
    2929
     
    6161  case 1: way = (firstBranch_ == 0) ? 1 : 0; break;
    6262  case 2: way = (firstBranch_ == 2) ? 1 : 2; break;
    63   default: jnlst_->Printf(J_WARNING, J_BRANCHING, "Warning, branchIndex_ has a strange value (%d)\n", branchIndex_);
     63  default: jnlst_->Printf(J_WARNING, J_BRANCHING,
     64                          "Warning, branchIndex_ has a strange value (%d)\n", branchIndex_);
    6465  }
    6566
    6667  // set lower or upper bound (round if this variable is integer)
    6768
     69  int  index   = brVar_ -> Index ();
     70  bool integer = brVar_ -> isInteger ();
     71
    6872  CouNumber
    69     l = solver -> getColLower () [index_],
    70     u = solver -> getColUpper () [index_];
     73    l = solver -> getColLower () [index],
     74    u = solver -> getColUpper () [index];
    7175
    7276  if (lcrop_ < l) lcrop_ = l;
     
    7579  switch (--way) { // from {0,1,2} to {-1,0,1}
    7680
    77   case -1: solver -> setColUpper (index_, integer_ ? floor (lcrop_) : lcrop_); break; // left
    78   case  0: solver -> setColLower (index_, integer_ ? ceil  (lcrop_) : lcrop_);
    79            solver -> setColUpper (index_, integer_ ? floor (rcrop_) : rcrop_); break; // central
    80   case  1: solver -> setColLower (index_, integer_ ? ceil  (rcrop_) : rcrop_); break; // right
     81  case -1: solver -> setColUpper (index, integer ? floor (lcrop_) : lcrop_); break; // left
     82  case  0: solver -> setColLower (index, integer ? ceil  (lcrop_) : lcrop_);
     83           solver -> setColUpper (index, integer ? floor (rcrop_) : rcrop_); break; // central
     84  case  1: solver -> setColLower (index, integer ? ceil  (rcrop_) : rcrop_); break; // right
    8185  default: jnlst_->Printf(J_WARNING, J_BRANCHING, "Couenne: branching on nonsense way %d\n", way);
    8286  }
     
    8690  if (jnlst_->ProduceOutput(J_DETAILED, J_BRANCHING)) { 
    8791    switch (way) {
    88     case -1: jnlst_->Printf(J_DETAILED, J_BRANCHING, "#3# Branch: x%d <= %g\n",               index_, lcrop_); break; // left
    89     case  0: jnlst_->Printf(J_DETAILED, J_BRANCHING, "#3# Branch: %g <= x%d <= %g\n", lcrop_, index_, rcrop_); break; // center
    90     case  1: jnlst_->Printf(J_DETAILED, J_BRANCHING, "#3# Branch: x%d >= %g\n",               index_, rcrop_); break; // right
     92    case -1: jnlst_->Printf(J_DETAILED, J_BRANCHING,
     93                            "#3# Branch: x%d <= %g\n",               index, lcrop_); break; // left
     94    case  0: jnlst_->Printf(J_DETAILED, J_BRANCHING,
     95                            "#3# Branch: %g <= x%d <= %g\n", lcrop_, index, rcrop_); break; // center
     96    case  1: jnlst_->Printf(J_DETAILED, J_BRANCHING,
     97                            "#3# Branch: x%d >= %g\n",               index, rcrop_); break; // right
    9198    default: jnlst_->Printf(J_DETAILED, J_BRANCHING, "Couenne: branching on nonsense way %d\n", way);
    9299    }
  • branches/Couenne/Couenne/src/branch/CouenneThreeWayBranchObj.hpp

    r969 r1017  
    3131  /// Constructor
    3232  CouenneThreeWayBranchObj (JnlstPtr jnlst,
    33                             int,
     33                            expression *,
    3434                            CouNumber,
    3535                            CouNumber,
    36                             int  = THREE_CENTER,
    37                             bool = false);
     36                            int  = THREE_CENTER
     37                            //bool = false
     38                            );
    3839
    3940  /// Copy constructor
    4041  CouenneThreeWayBranchObj (const CouenneThreeWayBranchObj &src):
    4142    OsiBranchingObject (src),
    42     index_ (src.index_),
     43    brVar_ (src.brVar_),
    4344    lcrop_ (src.lcrop_),
    4445    rcrop_ (src.rcrop_),
    4546    firstBranch_ (src.firstBranch_),
    46     integer_     (src.integer_),
    4747    jnlst_ (src.jnlst_){}
    4848
     
    6464  /// corresponding CouenneObject was created on w=f(x,y), it is
    6565  /// either x or y.
    66   int index_;
     66  expression *brVar_;
    6767
    6868  CouNumber lcrop_;  ///< left divider
     
    7373
    7474  /// True if the associated variable is integer
    75   bool integer_;
     75  //  bool integer_;
    7676
    7777  /// SmartPointer to the Journalist
  • branches/Couenne/Couenne/src/branch/infeasibility.cpp

    r1005 r1017  
    1414#include "CouenneThreeWayBranchObj.hpp"
    1515
    16 #include "exprGroup.hpp"
    17 #include "exprQuad.hpp"
    18 
    1916//#define DEBUG
    20 
    2117
    2218/// return difference between current value
     
    2723  // if selectBranch not called, choose one at random
    2824  whichWay_ = whichWay = TWO_LEFT;
    29   brVarInd_ = -1;
     25  brVar_ = NULL;
    3026
    3127  // infeasibility is always null for linear expressions
     
    4844
    4945  if ((delta                           < COUENNE_EPS) ||
    50       (delta / (1 + fabs (var + expr)) < COUENNE_EPS))
    51     {
     46      (delta / (1 + fabs (var + expr)) < COUENNE_EPS)) {
     47
    5248#if BR_TEST_LOG >= 0 && defined DEBUG
    53       if (reference_ -> Image () -> code () == COU_EXPRLOG) {
    54         printf ("---- found feasible point on curve: ");
    55         reference_ -> print (); printf (" := ");
    56         reference_ -> Image () -> print ();
    57         printf ("\n");
    58       }
     49    if (reference_ -> Image () -> code () == COU_EXPRLOG) {
     50      printf ("---- found feasible point on curve: ");
     51      reference_ -> print (); printf (" := ");
     52      reference_ -> Image () -> print ();
     53      printf ("\n");
     54    }
    5955#elif defined DEBUG
    60    printf ("----|%+g - %+g| = %+g  (delta=%+g) way %d, ind %d. ",  ////[%.2f,%.2f]
    61            var, expr,
    62            //       expression::Lbound (reference_ -> Index ()),
    63            //       expression::Ubound (reference_ -> Index ()),
    64            fabs (var - expr), delta, whichWay, brVarInd_);
    65    reference_             -> print (std::cout); std::cout << " = ";
    66    reference_ -> Image () -> print (std::cout); printf ("\n");
    67    return 0.;
     56    printf ("----|%+g - %+g| = %+g  (delta=%+g) way %d, ind %d. ",  ////[%.2f,%.2f]
     57            var, expr,
     58            //      expression::Lbound (reference_ -> Index ()),
     59            //      expression::Ubound (reference_ -> Index ()),
     60            fabs (var - expr), delta, whichWay, reference_ -> Index ());
     61    reference_             -> print (std::cout); std::cout << " = ";
     62    reference_ -> Image () -> print (std::cout); printf ("\n");
    6863#endif
    69       return 0.;
    70     }
     64    return 0.;
     65  }
    7166
    7267  // a nonlinear constraint w = f(x) is violated. The infeasibility
     
    10196
    10297  CouNumber improv = reference_ -> Image () ->
    103     selectBranch (this, info,                   // input parameters
    104                   brVarInd_, brPts_, whichWay); // result: who, where, and how to branch
     98    selectBranch (this, info,                // input parameters
     99                  brVar_, brPts_, whichWay); // result: who, where, and how to branch
    105100
    106   if (brVarInd_ >= 0) {
     101  if (brVar_) {
    107102
    108103#ifdef DEBUG
     
    114109#endif
    115110
    116     delta = improv;
     111    if (improv > COUENNE_EPS) delta = improv;
    117112    whichWay_ = whichWay;
    118113
    119114    // TODO: test this, if commented gives A LOT of problems in nvs24
    120115
    121     if (info -> lower_ [brVarInd_] >=
    122         info -> upper_ [brVarInd_] - COUENNE_EPS) {
     116    int index = brVar_ -> Index ();
     117
     118    if (info -> lower_ [index] >=
     119        info -> upper_ [index] - COUENNE_EPS) {
    123120
    124121#ifdef DEBUG
    125122      printf ("### warning, tiny bounding box [%g,%g] for x_%d\n",
    126               info -> lower_ [brVarInd_],
    127               info -> upper_ [brVarInd_], brVarInd_);
     123              info -> lower_ [index],
     124              info -> upper_ [index], index);
    128125#endif
    129126
    130       delta = 0;
     127      delta = 0.;
    131128    }
    132129  }
     
    142139#endif
    143140          var, expr,
    144           //        expression::Lbound (reference_ -> Index ()),
    145           //        expression::Ubound (reference_ -> Index ()),
    146           fabs (var - expr), delta, whichWay, brVarInd_
     141          fabs (var - expr), delta, whichWay, brVar_ -> Index ()
    147142#if BR_TEST_LOG >= 0
    148143          , info -> lower_ [BR_TEST_LOG]
  • branches/Couenne/Couenne/src/branch/operators/branchExprAbs.cpp

    r900 r1017  
    1616/// set up branching object by evaluating branching points for each
    1717/// expression's arguments. For an exprAbs, simply branch at zero.
    18 CouNumber exprAbs::selectBranch (const CouenneObject *obj, 
     18CouNumber exprAbs::selectBranch (const CouenneObject *obj,
    1919                                 const OsiBranchingInformation *info,
    20                                  int &ind,
    21                                  double * &brpts, 
     20                                 expression * &var,
     21                                 double * &brpts,
    2222                                 int &way) {
    23   ind = argument_ -> Index ();
     23
     24  var = argument_;
     25
     26  int ind = var -> Index ();
    2427
    2528  assert ((ind >= 0) && (obj -> Reference () -> Index () >= 0));
  • branches/Couenne/Couenne/src/branch/operators/branchExprDiv.cpp

    r895 r1017  
    1717CouNumber exprDiv::selectBranch (const CouenneObject *obj,
    1818                                 const OsiBranchingInformation *info,
    19                                  int &ind,
     19                                 expression *&var,
    2020                                 double * &brpts,
    2121                                 int &way) {
     
    4040  if ((yl < 0) && (yu > 0)) {
    4141
    42     ind = yi;
     42    var = arglist_ [1];
     43    //ind = yi;
    4344
    4445    way = TWO_RAND;
     
    6667      (yu >  COUENNE_INFINITY)) {
    6768
    68     ind = yi;
     69    var = arglist_ [1];//ind = yi;
    6970    brpts = (double *) realloc (brpts, sizeof (CouNumber));
    7071
     
    9192      (wu >  COUENNE_INFINITY)) {
    9293
    93     ind = wi;
     94    var = obj -> Reference ();//ind = wi;
    9495
    9596    if ((wl < -COUENNE_INFINITY) &&
     
    147148
    148149  if (dx > dy)
    149     if (dx > dw) {ind = xi; *brpts = (xl + xu) / 2.; return fabs (x0 - y0*w0);} // dx maximum
    150     else         {ind = wi; *brpts = (wl + wu) / 2.; return fabs (w0 - x0/y0);} // dw
     150    if (dx > dw) {var = arglist_[0];      *brpts = (xl+xu)/2.; return fabs (x0-y0*w0);}
     151    else         {var = obj->Reference(); *brpts = (wl+wu)/2.; return fabs (w0-x0/y0);}
    151152  else
    152     if (dy > dw) {ind = yi; *brpts = (yl + yu) / 2.; return fabs (y0 - x0/w0);} // dy
    153     else         {ind = wi; *brpts = (wl + wu) / 2.; return fabs (w0 - x0/y0);} // dw
     153    if (dy > dw) {var = arglist_[1];      *brpts = (yl+yu)/2.; return fabs (y0-x0/w0);}
     154    else         {var = obj->Reference(); *brpts = (wl+wu)/2.; return fabs (w0-x0/y0);}
    154155}
  • branches/Couenne/Couenne/src/branch/operators/branchExprExp.cpp

    r895 r1017  
    2121CouNumber exprExp::selectBranch (const CouenneObject *obj,
    2222                                 const OsiBranchingInformation *info,
    23                                  int &ind,
     23                                 expression *&var,
    2424                                 double * &brpts,
    2525                                 int &way) {
     
    4040  // propagated bounds will do the rest.
    4141
    42   ind    = argument_           -> Index ();
    43   int wi = obj -> Reference () -> Index ();
     42  var = argument_;
     43
     44  int
     45    ind = var -> Index (),
     46    wi  = obj -> Reference () -> Index ();
    4447
    4548  assert ((ind >= 0) && (wi >= 0));
     
    110113
    111114  simpletriplet ft (exp, exp, exp, log);
    112 
    113   switch (obj -> Strategy ()) {
    114 
    115   case CouenneObject::MIN_AREA:     *brpts = maxHeight   (&ft, x0, y0, l, u); break;
    116   case CouenneObject::BALANCED:     *brpts = minMaxDelta (&ft, x0, y0, l, u); break;
    117   case CouenneObject::MID_INTERVAL:
    118   default:                          *brpts = midInterval (     x0,     l, u); break;
    119   }
    120 
    121   /*  simpletriplet ft (exp, exp, exp, log);
    122 
    123   *brpts = (u < l + COUENNE_EPS) ?
    124     midInterval (powNewton (x0, y0, exp, exp, exp), l, u) : // find closest point on curve
    125     maxHeight (&ft, l, u); // apply minarea
    126   // minmaxdistance:  minMaxDelta (&ft, x0, l, u);
    127   */
     115  *brpts = obj -> getBrPoint (&ft, x0, y0, l, u);
    128116
    129117  way = TWO_RAND;
     
    132120  return CoinMin (projectSeg (x0, y0, l, exp (l), *brpts, exp (*brpts),             -1),
    133121                  projectSeg (x0, y0,             *brpts, exp (*brpts), u, exp (u), -1));
    134 
    135   /*  x0 -= *brpts;
    136   y0 -= exp (*brpts);
    137   return sqrt (x0*x0 + y0*y0);*/
    138122}
  • branches/Couenne/Couenne/src/branch/operators/branchExprInv.cpp

    r895 r1017  
    1818
    1919/// generic approach for negative powers (commom with exprInv::selectBranch
    20 CouNumber negPowSelectBranch (const CouenneObject *obj, int ind,
     20CouNumber negPowSelectBranch (const CouenneObject *obj,
    2121                              double * &brpts,
    2222                              int &way,
     
    147147
    148148  powertriplet ft (k);
    149 
    150   switch (obj -> Strategy ()) {
     149  *brpts = obj -> getBrPoint (&ft, x0, y0, l, u);
     150
     151  /*switch (obj -> Strategy ()) {
    151152
    152153  case CouenneObject::MIN_AREA:     *brpts = maxHeight   (&ft, x0, y0, l, u); break;
     
    154155  case CouenneObject::MID_INTERVAL:
    155156  default:                          *brpts = midInterval (     x0,     l, u); break;
    156   }
     157  }*/
    157158
    158159  /*
     
    185186CouNumber exprInv::selectBranch (const CouenneObject *obj,
    186187                                 const OsiBranchingInformation *info,
    187                                  int &ind,
     188                                 expression *&var,
    188189                                 double * &brpts,
    189190                                 int &way) {
    190191
    191   ind    = argument_           -> Index ();
    192   int wi = obj -> Reference () -> Index ();
     192  var = argument_;
     193
     194  int
     195    ind = argument_           -> Index (),
     196    wi  = obj -> Reference () -> Index ();
    193197
    194198  assert ((ind >= 0) && (wi >= 0));
     
    199203            u  = info -> upper_    [ind];
    200204
    201   return negPowSelectBranch (obj, ind, brpts, way, -1, x0, y0, l,  u);
     205  return negPowSelectBranch (obj, brpts, way, -1, x0, y0, l,  u);
    202206}
  • branches/Couenne/Couenne/src/branch/operators/branchExprLog.cpp

    r959 r1017  
    2323CouNumber exprLog::selectBranch (const CouenneObject *obj,
    2424                                 const OsiBranchingInformation *info,
    25                                  int &ind,
     25                                 expression *&var,
    2626                                 double * &brpts,
    2727                                 int &way) {
     
    4444  // propagated bounds will do the rest.
    4545
    46   ind    = argument_ -> Index ();
    47   int wi = obj -> Reference () -> Index ();
     46  var = argument_;
     47
     48  int
     49    ind = var -> Index (),
     50    wi  = obj -> Reference () -> Index ();
    4851
    4952  assert ((ind >= 0) && (wi >= 0));
     
    5659
    5760  if (u < COUENNE_EPS) { // strange case, return default branching rule
    58     ind = -1;
    59     return 0;
     61    var = NULL;
     62    return 0.;
    6063  }
    6164
     
    134137  simpletriplet ft (log, inv, oppInvSqr, inv);
    135138
    136   switch (obj -> Strategy ()) {
    137 
    138   case CouenneObject::MIN_AREA:     *brpts = maxHeight   (&ft, x0, y0, l, u); break;
    139   case CouenneObject::BALANCED:     *brpts = minMaxDelta (&ft, x0, y0, l, u); break;
    140   case CouenneObject::MID_INTERVAL:
    141   default:                          *brpts = midInterval (     x0,     l, u); break;
    142   }
     139  *brpts = obj -> getBrPoint (&ft, x0, y0, l, u);
    143140
    144141  //  *brpts = midInterval (powNewton (x0, y0, log, inv, oppInvSqr), l, u);
  • branches/Couenne/Couenne/src/branch/operators/branchExprMul.cpp

    r915 r1017  
    2222CouNumber exprMul::selectBranch (const CouenneObject *obj,
    2323                                 const OsiBranchingInformation *info,
    24                                  int &ind,
     24                                 expression *&var,
    2525                                 double * &brpts,
    2626                                 int &way) {
     
    4141  // make them pretty hard to deal with
    4242
    43   if (((ind = xi) >= 0) && (xl < - COUENNE_INFINITY) && (xu > COUENNE_INFINITY) ||
    44       ((ind = yi) >= 0) && (yl < - COUENNE_INFINITY) && (yu > COUENNE_INFINITY)) {
     43  if (((var = arglist_ [0]) -> Index() >= 0) && (xl < -COUENNE_INFINITY) && (xu > COUENNE_INFINITY) ||
     44      ((var = arglist_ [1]) -> Index() >= 0) && (yl < -COUENNE_INFINITY) && (yu > COUENNE_INFINITY)) {
    4545
    4646    // branch around current point. If it is also at a crazy value,
     
    4848
    4949    brpts = (double *) realloc (brpts, 2 * sizeof (double));
    50     CouNumber curr = expression::Variable (ind);
     50    CouNumber curr = (*var) ();//expression::Variable (ind);
    5151
    5252    if (fabs (curr) >= LARGE_BOUND) curr = 0;
     
    6969  // don't privilege xi over yi
    7070
     71  int ind;
     72
    7173  if (xl<-COUENNE_INFINITY) {ind=xi; *brpts=midInterval(((x0<0)?2:0.5)*x0,xl,xu); way=TWO_RIGHT;} else
    7274  if (xu> COUENNE_INFINITY) {ind=xi; *brpts=midInterval(((x0>0)?2:0.5)*x0,xl,xu); way=TWO_LEFT;}  else
     
    7678                            {ind=xi; *brpts=midInterval(x0, xl, xu);              way=TWO_RAND;}
    7779
     80  var = arglist_ [(ind == xi) ? 0 : 1];
     81
    7882  return fabs (w0 - x0*y0);
    7983}
  • branches/Couenne/Couenne/src/branch/operators/branchExprPow.cpp

    r959 r1017  
    1313#include "CouenneObject.hpp"
    1414#include "CouenneBranchingObject.hpp"
    15 
    1615#include "projections.hpp"
    1716#include "funtriplets.hpp"
     
    1918
    2019/// generic approach for negative powers (commom with exprInv::selectBranch
    21 CouNumber negPowSelectBranch (const CouenneObject *obj, int ind,
     20CouNumber negPowSelectBranch (const CouenneObject *obj,
    2221                              double * &brpts,
    2322                              int &way,
     
    3130CouNumber exprPow::selectBranch (const CouenneObject *obj,
    3231                                 const OsiBranchingInformation *info,
    33                                  int &ind,
     32                                 expression *&var,
    3433                                 double * &brpts,
    3534                                 int &way) {
     
    3837  // of the form x^k
    3938
    40   ind    = arglist_ [0]        -> Index ();
    41   int wi = obj -> Reference () -> Index ();
     39  var = arglist_ [0];
     40
     41  int
     42    ind = arglist_ [0]        -> Index (),
     43    wi  = obj -> Reference () -> Index ();
    4244
    4345  assert ((ind >= 0) && (wi >= 0) && (arglist_ [1] -> Type () == CONST));
     
    4951            l  = info -> lower_    [ind],
    5052            u  = info -> upper_    [ind];
     53
     54  /*printf ("selbra for "); print ();
     55  printf ("%g [%g,%g] -> %g [%g,%g]\n",
     56          x0, l, u, y0,
     57          info -> lower_    [wi],
     58          info -> upper_    [wi]);*/
    5159
    5260  if      (x0 < l) x0 = l;
     
    6169
    6270  // case 1: k negative, resort to method similar to exprInv:: ///////////////////////////////
    63   if (k<0) return negPowSelectBranch (obj, ind, brpts, way, k, x0, y0, l, u);
     71  if (k<0) return negPowSelectBranch (obj, brpts, way, k, x0, y0, l, u);
    6472
    6573  int intk = 0;
     
    8694        y0 -= pow (*brpts, k);
    8795
    88         return sqrt (x0*x0 + y0 * y0); // exact distance
     96        return sqrt (x0*x0 + y0*y0); // exact distance
    8997      }
    9098
     
    146154
    147155    powertriplet ft (k);
    148     *brpts = maxHeight (&ft, x0, y0, l, u);
     156    //*brpts = maxHeight (&ft, x0, y0, l, u);
     157    *brpts = obj -> getBrPoint (&ft, x0, y0, l, u);
    149158
    150159    way = (x0 < *brpts) ? TWO_LEFT : TWO_RIGHT;
    151160
    152     /*
    153     w -> print (); printf (" := "); dynamic_cast <exprAux *> (w) -> Image () -> print ();
     161    //    w -> print (); printf (" := "); dynamic_cast <exprAux *> (w) -> Image () -> print ();
     162    /*print ();
    154163    printf (" (%g,%g) [%g,%g] --> %g, distances = %g,%g\n",
    155164            x0, y0, l, u, *brpts,
    156165            projectSeg (x0, y0, l,      safe_pow (l,k), *brpts, safe_pow (*brpts,k), 0),
    157             projectSeg (x0, y0, *brpts, safe_pow (*brpts,k),      u, safe_pow (u,k), 0));
    158     */
     166            projectSeg (x0, y0, *brpts, safe_pow (*brpts,k),      u, safe_pow (u,k), 0));*/
    159167
    160168    // Min area  -- exact distance
    161     return CoinMin (projectSeg (x0, y0, l,      safe_pow (l,k), *brpts, safe_pow (*brpts,k), 0),
    162                     projectSeg (x0, y0, *brpts, safe_pow (*brpts,k),      u, safe_pow (u,k), 0));
     169    CouNumber retval =
     170      CoinMin (projectSeg (x0, y0, l,      safe_pow (l,k), *brpts, safe_pow (*brpts,k), 0),
     171               projectSeg (x0, y0, *brpts, safe_pow (*brpts,k),      u, safe_pow (u,k), 0));
     172
     173    //printf ("returning %.10e\n", retval);
     174
     175    return retval;
    163176
    164177    /*      brpts = (double *) realloc (brpts, sizeof (double));
     
    324337  // picks the default variable/branchingpoint for the expression
    325338
    326   ind = -1;
     339  var = NULL;
    327340  return 0.;
    328341}
  • branches/Couenne/Couenne/src/branch/operators/branchExprQuad.cpp

    r895 r1017  
    1414#include "CouenneBranchingObject.hpp"
    1515
    16 
    1716//#define DEBUG
    1817
     
    2120CouNumber exprQuad::selectBranch (const CouenneObject *obj,
    2221                                  const OsiBranchingInformation *info,
    23                                   int &ind,
     22                                  expression *&var,
    2423                                  double * &brpts,
    2524                                  int &way) {
    2625
    27   const CouNumber *x = info -> solution_,
    28                   *l = info -> lower_,
    29                   *u = info -> upper_;
     26  int ind = -1;
    3027
    31   CouNumber delta    = (*(obj -> Reference ())) () - (*this) (),
    32            *alpha    = (delta < 0) ? dCoeffLo_ : dCoeffUp_,
    33             maxcontr = -COUENNE_INFINITY;
     28  // use a combination of eigenvectors and bounds
    3429
    35   int bestVar =
    36     nqterms_ ? *qindexI_ :
    37     nlterms_ ? *index_   : -1; // initialize it to something very default
     30  CouNumber delta = (*(obj -> Reference ())) () - (*this) ();
    3831
    39   if (!alpha) { // no convexification available,
    40                 // find the largest interval
     32  /*printf ("infeasibility: ");
     33  obj -> Reference () -> print ();
     34  printf (" [%g=%g] := ",
     35          (*(obj -> Reference ())) (), info -> solution_ [obj -> Reference () -> Index ()]);
     36          print (); printf (" [%g]\n", (*this) ());*/
    4137
    42     CouNumber maxcontr = -COUENNE_INFINITY; // maximum contribution to current value
     38  brpts = (double *) realloc (brpts, sizeof (double));
     39  way = TWO_RAND;
    4340
    44     int bestVar = -1;
     41  // depending on where the current point is w.r.t. the curve,
     42  // branching is done on the eigenvector corresponding to the minimum
     43  // or maximum eigenvalue
    4544
    46     if (dIndex_) { // use indices in dIndex_ -- TODO: dIndex_ should
    47                    // be created by constructor
     45  std::vector <std::pair <CouNumber,
     46    std::vector <std::pair <exprVar *, CouNumber> > > >::iterator         fi = eigen_.begin ();
    4847
    49       for (int i=nDiag_; i--;) {
     48  std::vector <std::pair <CouNumber,
     49    std::vector <std::pair <exprVar *, CouNumber> > > >::reverse_iterator ri = eigen_.rbegin ();
    5050
    51         int ind = dIndex_ [i];
    52         CouNumber diff;
     51  CouNumber max_span = -COUENNE_INFINITY;
    5352
    54 #ifdef DEBUG
    55         printf ("[%10g %10g %10g] %4d %10g\n", x [ind], l [ind], u [ind], bestVar, maxcontr);
    56 #endif
     53  bool changed_sign = false;
    5754
    58         //if ((diff = CoinMin (xi - l [qi], u [qi] - xi)) > maxcontr) {bestVar = qi; maxcontr = diff;}
    59         if ((diff = u [ind] - l [ind])                  > maxcontr) {bestVar = ind; maxcontr = diff;}
     55  /////////////////////////////////////////////////////////
     56
     57  for (;
     58       ((delta < 0.) && (fi != eigen_. end  ()) || // && (fi -> first < 0.) ||
     59        (delta > 0.) && (ri != eigen_. rend ()));  // && (ri -> first > 0.));
     60
     61       ++fi, ++ri) {
     62
     63    std::vector <std::pair <exprVar *, CouNumber> > &ev =
     64      (delta < 0.) ? fi -> second : ri -> second;
     65
     66    for (std::vector <std::pair <exprVar *, CouNumber> >::iterator j = ev.begin ();
     67         j != ev.end (); ++j) {
     68
     69      int index = j -> first -> Index ();
     70
     71      CouNumber
     72        lb = info -> lower_ [index],
     73        ub = info -> upper_ [index];
     74
     75      if ((fabs (ub-lb) > COUENNE_EPS) ||
     76          // no variable was found but the eigenvalue is already
     77          // positive (negative)
     78          ((changed_sign = true) &&
     79           ((delta < 0.) && (fi -> first > 0.) ||
     80            (delta > 0.) && (ri -> first < 0.)) &&
     81           (max_span < 0.))) {
     82
     83        CouNumber span = -1;
     84
     85        if ((ub-lb > COUENNE_INFINITY) ||
     86            ((span = (ub-lb) * fabs (j -> second)) > max_span + COUENNE_EPS)) {
     87
     88          ind = index;
     89          var = j -> first;
     90
     91          *brpts = midInterval (info -> solution_ [index], lb, ub);
     92
     93          if (changed_sign)
     94            break;
     95
     96          if (span >= 0) max_span = span; // finite, keep searching
     97          else break;                     // span unbounded, stop searching
     98        }
    6099      }
    61     } else
    62       for (int i=nqterms_; i--;) {
     100    }
     101  }
    63102
    64         // find largest (in abs. value) coefficient with at least one
    65         // index within bounds
     103  if ((eigen_.size () == 0) // if no eigenvalue has been computed yet
     104      || (ind < 0)) {       // or no index was found, pick largest interval
    66105
    67         int qi = qindexI_ [i],
    68             qj = qindexJ_ [i];
     106    CouNumber max_span = -COUENNE_INFINITY;
    69107
    70         CouNumber
    71           xi = x [qi],
    72           xj = x [qj], diff;
     108    for (std::map <exprVar *, std::pair <CouNumber, CouNumber> >::iterator i = bounds_. begin ();
     109         i != bounds_. end (); ++i) {
    73110
    74         if ((diff = CoinMin (xi - l [qi], u [qi] - xi)) > maxcontr) {bestVar = qi; maxcontr = diff;}
    75         if ((diff = CoinMin (xj - l [qj], u [qj] - xj)) > maxcontr) {bestVar = qj; maxcontr = diff;}
     111      CouNumber
     112        lb = i -> second.first,
     113        ub = i -> second.second,
     114        span = ub - lb;
     115
     116      if ((span > COUENNE_EPS) &&
     117          (span > max_span)) {
     118
     119        var = i -> first;
     120        ind = var -> Index ();
     121        *brpts = midInterval (info -> solution_ [ind], lb, ub);   
    76122      }
     123    }
    77124
    78     ind = bestVar;
    79     brpts = (double *) realloc (brpts, sizeof (double));
    80     //    *brpts = x [bestVar];
    81     *brpts = (l [bestVar] + u [bestVar]) / 2;
    82     way = TWO_RAND;
     125    if (ind < 0) {
    83126
    84 #ifdef DEBUG
    85     printf ("brExprQuad (%s): |delta| = %g, brpt = %g (%g), var = x%d, way = %d -- NO alpha\n",
    86             (delta < 0) ? "lower" : "upper",
    87             fabs (delta), *brpts, x [bestVar], bestVar, way);
    88 #endif
     127      var = obj -> Reference ();
     128      ind = var -> Index ();
     129
     130      *brpts = midInterval (info -> solution_ [ind],
     131                            info -> lower_ [ind],
     132                            info -> upper_ [ind]);
     133    }
    89134
    90135    return fabs (delta);
    91136  }
    92137
    93   int *indices = dIndex_;
    94 
    95   // there is a convexification already, find i = argmin {alpha_i (x_i
    96   // - l_i) (u_i - x_i)}
    97 
    98   for (int i=0; i < nDiag_; i++, indices++) {
    99    
    100     CouNumber curx = x [*indices],
    101       contrib = *alpha++ *
    102       (curx         - l [*indices]) *
    103       (u [*indices] - curx);
    104 
    105     if (fabs (contrib) > maxcontr) {
    106       bestVar  = *indices;
    107       maxcontr = contrib;
    108     }
    109   }
    110 
    111   ind = bestVar;
    112   way = TWO_RAND;
    113 
    114   brpts = (double *) realloc (brpts, sizeof (double));
    115 
    116   *brpts = midInterval (x [bestVar], l [bestVar], u [bestVar]);
    117 
    118   /*  if ((*brpts > ub - COUENNE_NEAR_BOUND) ||
    119       (*brpts < lb + COUENNE_NEAR_BOUND))
    120 
    121       *brpts = 0.5 * (lb + ub);*/
    122 
    123 #ifdef DEBUG
    124   printf ("brExprQuad: |delta| = %g, brpt = %g (%g), var = x%d, way = %d\n",
    125           fabs (delta), *brpts, x [bestVar], bestVar, way);
    126 #endif
    127 
    128138  return fabs (delta);
    129139}
  • branches/Couenne/Couenne/src/branch/operators/branchExprSinCos.cpp

    r895 r1017  
    1818CouNumber trigSelBranch (const CouenneObject *obj,
    1919                         const OsiBranchingInformation *info,
    20                          int &ind,
     20                         expression *&var,
    2121                         double * &brpts,
    2222                         int &way,
     
    2525  // for now, apply default branching rule
    2626
    27   ind = -1;
     27  var = NULL;
    2828  return 0.;
    2929
    30   // minarea
     30  // TODO: minarea
    3131}
  • branches/Couenne/Couenne/src/branch/operators/minMaxDelta.cpp

    r915 r1017  
    1111
    1212//#define TEST_BRPTS
    13 #ifdef  TEST_BRPTS
     13/*#ifdef  TEST_BRPTS
    1414
    1515#include <stdio.h>
     
    135135
    136136
    137 #else
    138 #include "CouenneBranchingObject.hpp"
     137#else */
     138#include "CouenneObject.hpp"
    139139#include "funtriplets.hpp"
    140 #endif
     140//#endif
    141141
    142142
     
    259259CouNumber maxHeight (funtriplet *ft,
    260260                     CouNumber x, CouNumber y,
    261                      CouNumber lb, CouNumber ub)
    262 {
     261                     CouNumber lb, CouNumber ub) {
    263262  /* fprintf (stderr,"slope is (%g - %g) / (%g - %g) = %g / %g = %g ----> inverse is %g\n",
    264263          ft -> F (ub),
     
    269268          (ft -> F (ub) - ft -> F (lb)) / (ub - lb),
    270269          ft -> FpInv ((ft -> F (ub) - ft -> F (lb)) / (ub - lb)));*/
    271 return (ft -> FpInv ((ft -> F (ub) - ft -> F (lb)) / (ub - lb)));}
     270  return (ft -> FpInv ((ft -> F (ub) - ft -> F (lb)) / (ub - lb)));
     271}
    272272
    273273
     
    275275// testing stuff ////////////////////////////////////////////////////////////////////////
    276276
    277 #ifdef TEST_BRPTS
     277/*#ifdef TEST_BRPTS
    278278
    279279inline CouNumber safe_pow (register CouNumber base,
     
    367367  fprintf (stderr, "maxHeight = %g\n", b);
    368368}
    369 #endif
     369#endif*/
  • branches/Couenne/Couenne/src/branch/projections.cpp

    r895 r1017  
    1010#include "CouenneTypes.hpp"
    1111#include "CouennePrecisions.hpp"
    12 
    13 //#define DEBUG
    1412
    1513/*  compute projection of point (x0, y0) on the segment defined by
  • branches/Couenne/Couenne/src/convex/operators/alphaConvexify.cpp

    r959 r1017  
    88 */
    99
    10 #include "exprQuad.hpp"
    1110
    1211#include "CoinHelperFunctions.hpp"
    13 
    1412#include "OsiSolverInterface.hpp"
    1513#include "IpLapack.hpp"
    1614
     15#include "exprQuad.hpp"
     16#include "CouenneProblem.hpp"
     17
    1718//#define DEBUG
    18 
    19 // fill in one of the two dCoeff vectors
    20 void fill_dCoeff (CouNumber * &, CouNumber, CouNumber *, int);
    21 
    2219
    2320/** Computes alpha coefficients for an alpha under- and overestimator
     
    3734 * If Q is positive semidefinite, then dCoeffLo_ will be NULL.
    3835 * If Q is negative semidefinite, then dCoeffUp_ will be NULL.
    39 */
    40 
    41 void exprQuad::alphaConvexify (const OsiSolverInterface &si) {
    42 
    43   if (nqterms_ == 0) {
    44     nDiag_ = 0;
    45     return;
    46   }
    47 
    48   // inverse of dIndex_ mapping: for each variable tell me the
    49   // index that it will have in dIndex_, or -1 if not there
    50 
    51   int* indexmap = new int [si.getNumCols ()];
     36 *
     37 * Return true if a convexification is there, false otherwise.
     38 */
     39
     40bool exprQuad::alphaConvexify (const CouenneProblem *p,
     41                               const OsiSolverInterface &si) {
     42
     43  if (matrix_.size () == 0)
     44    return false;
     45
     46  // inverse of dIndex_ mapping: for each variable tell me the index
     47  // that it will have in dIndex_, or -1 if not there
     48
     49  int k=0,
     50     nDiag    = bounds_.size (),
     51    *indexmap = new int [si.getNumCols ()],
     52    *indices  = new int [nDiag];
     53
    5254  CoinFillN (indexmap, si.getNumCols (), -1);
    5355
    54   if (dIndex_ == NULL)
    55     make_dIndex (si.getNumCols (), indexmap);
    56   else
    57     // build indexmap as inverse of dIndex_
    58     for (int i=0; i<nDiag_; ++i)
    59       indexmap [dIndex_ [i]] = i;
    60 
    6156  // box diameter
    62   double *diam = new double [nDiag_];
    63 
    64   const double
    65     *lower = si.getColLower (),
    66     *upper = si.getColUpper ();
    67 
    68   for (int i=0; i<nDiag_; ++i) {
    69 
    70     int di = dIndex_ [i];
    71     diam [i] = upper [di] - lower [di];
     57  double *diam = new double [nDiag];
     58
     59  bool changed_bounds = false;
     60
     61  for (std::map <exprVar *, std::pair <CouNumber, CouNumber> >::iterator i = bounds_.begin ();
     62       i != bounds_.end (); ++i, k++) {
     63
     64#ifdef DEBUG
     65  printf ("b%04d. [%20g, %20g]\n", i->first->Index(), i->second.first, i->second.second);
     66#endif
     67
     68    int index = i -> first -> Index ();
     69    indexmap [index] = k;
     70    indices [k] = index;
     71
     72    CouNumber
     73      lb = (*(i -> first -> Lb ())) (),//lower [di],
     74      ub = (*(i -> first -> Ub ())) ();//upper [di];
     75
     76    // if one variable unbounded, bail out
     77    if ((lb < -COUENNE_INFINITY) ||
     78        (ub >  COUENNE_INFINITY)) {
     79
     80      delete [] diam;
     81      delete [] indexmap;
     82      delete [] indices;
     83
     84      return false;
     85    }
     86
     87    // if no variable has changed bounds, no need to convexify
     88    if (fabs (lb - i->second.first)  > COUENNE_EPS) {i -> second.first  = lb; changed_bounds = true;}
     89    if (fabs (ub - i->second.second) > COUENNE_EPS) {i -> second.second = ub; changed_bounds = true;}
     90
     91    diam [k] = ub - lb;
     92#ifdef DEBUG
     93    printf ("diam %4d - %4d = %g - %g = %g\n", index, k, ub, lb, diam [k]);
     94#endif
     95  }
     96
     97  if (!changed_bounds) {
     98
     99    delete [] diam;
     100    delete [] indexmap;
     101    delete [] indices;
     102
     103    return true;
    72104  }
    73105
    74106  // lower triangular of quadratic term matrix, scaled by box diameter
    75107
    76   double *matrix = new double [nDiag_ * nDiag_];
    77 
    78   CoinFillN (matrix, nDiag_ * nDiag_, 0.);
    79 
    80   for (int i=0; i < nqterms_; ++i) {
    81 
    82     int row = indexmap [qindexI_ [i]];
    83     int col = indexmap [qindexJ_ [i]];
    84 
    85     // compute value of matrix entry = q_ij * (u_i-l_i) * (u_j-l_j)
    86     // I (Stefan) do not understand the Lapack docu; it says it needs
    87     // only the lower triangular but it seem to need both parts to
    88     // work correct
    89     double cell = qcoeff_ [i] * diam [row] * diam [col];
    90 
    91     matrix          [col * nDiag_ + row] = cell;
    92     if (row != col)
    93       matrix        [row * nDiag_ + col] = cell;
    94     //    printf("row %d, col %d: %f\n", row, col, matrix[col*nDiag_+row]);
    95   }
     108  double *matrix = new double [nDiag * nDiag];
     109
     110  CoinFillN (matrix, nDiag * nDiag, 0.);
     111
     112  for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
     113
     114    int
     115      xind = row -> first -> Index (),
     116      irow = indexmap [xind];
     117
     118    for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
     119
     120      int
     121        yind = col -> first -> Index (),
     122        icol = indexmap [yind];
     123
     124      double cell = col -> second * diam [irow] * diam [icol];
     125
     126      matrix          [icol * nDiag + irow] = cell;
     127      if (irow != icol)
     128        matrix        [irow * nDiag + icol] = cell;
     129    }
     130  }
     131
     132  // compute value of matrix entry = q_ij * (u_i-l_i) * (u_j-l_j)
     133  // I (Stefan) do not understand the Lapack docu; it says it needs
     134  // only the lower triangular but it seem to need both parts to
     135  // work correct
    96136
    97137  delete [] indexmap;
    98138
    99139  // compute minimum and maximum eigenvalue of matrix
    100   double* eigval = new double [nDiag_];
     140  double* eigval = new double [nDiag];
    101141  int info;
    102142
    103   Ipopt::IpLapackDsyev (false,  // do not compute eigenvector
    104                         nDiag_, // dimension
     143  /*printf ("nDiag = %d\n", nDiag);
     144  for (int i=0; i<nDiag; i++) {
     145    for (int j=0; j<nDiag; j++)
     146      printf ("%6.2f ", matrix [i*nDiag + j]);
     147    printf ("\n");
     148    }*/
     149
     150  Ipopt::IpLapackDsyev (true,   // compute eigenvector
     151                        nDiag,  // dimension
    105152                        matrix, // matrix
    106                         nDiag_, // "leading dimension" (number of columns, I think)
     153                        nDiag, // "leading dimension" (number of columns, I think)
    107154                        eigval, // output vector to store eigenvalues
    108155                        info);  // output status variable
     156
     157  if (info != 0) {
     158    printf ("exprQuad::alphaConvexify, warning: problem computing eigenvalue, info=%d\n", info);
     159    return false;
     160    //TODO error handling
     161  }
     162
     163  // clean eigenvector structure
     164  eigen_.erase (eigen_.begin (), eigen_.end ());
     165
     166  for (int i=0; i<nDiag; i++) {
     167
     168    std::pair <CouNumber, std::vector <std::pair <exprVar *, CouNumber> > > eigenCoord;
     169
     170    eigenCoord. first = eigval [i];
     171
     172    for (int j=0; j<nDiag; j++) {
     173
     174      CouNumber elem = matrix [i * nDiag + j];
     175
     176      if (fabs (elem) > COUENNE_EPS)
     177        eigenCoord. second. push_back (std::pair <exprVar *, CouNumber>
     178                                       (p -> Var (indices [j]), elem));
     179    }
     180
     181    eigen_.push_back (eigenCoord);
     182  }
     183
     184#ifdef DEBUG
     185  for (std::vector <std::pair <CouNumber,
     186         std::vector <std::pair <exprVar *, CouNumber> > > >::iterator i = eigen_.begin ();
     187       i != eigen_.end (); ++i) {
     188    printf (" [%g] -- ", i -> first);
     189    for (std::vector <std::pair <exprVar *, CouNumber> >::iterator j = i -> second. begin();
     190         j != i -> second. end (); ++j)
     191      printf ("(%d,%g) ", j -> first -> Index (), j -> second);
     192    printf ("\n");
     193  }
     194#endif
     195
     196  delete [] indices;
    109197  delete [] matrix;
    110 
    111   if (info != 0) {
    112     printf ("exprQuad::alphaConvexify: problem computing eigenvalue, info=%d\n", info);
    113     exit (-1);
    114     //TODO error handling
    115   }
    116 
    117   // if min. eigenvalue negative, setup dCoeffLo_
    118   if (eigval [0] < 0)
    119     fill_dCoeff (dCoeffLo_, eigval [0], diam, nDiag_);
    120   else  // quadratic term is convex, no convexification needed
    121     if (dCoeffLo_) {
    122       delete dCoeffLo_;
    123       dCoeffLo_ = NULL;
    124     }
    125 
    126   // if max. eigenvalue is positive, setup dCoeffUp_
    127   if (eigval [nDiag_ - 1] > 0)
    128     fill_dCoeff (dCoeffUp_, eigval [nDiag_ - 1], diam, nDiag_);
    129   else // quadratic term is concave, no "concavification" needed
    130     if (dCoeffUp_) {
    131       delete dCoeffUp_;
    132       dCoeffUp_ = NULL;
    133     }
    134 
    135198  delete [] diam;
    136199  delete [] eigval;
     200
     201  return true;
    137202}
    138 
    139 
    140 // fill diagonal vector using eigenvalue passed as parameter
    141 void fill_dCoeff (CouNumber * &dCoeff, CouNumber eigval, CouNumber *diam, int n) {
    142 
    143   if (dCoeff == NULL)
    144     dCoeff = new CouNumber [n];
    145 
    146   for (int i=n; i--;) {
    147     CouNumber di = diam [i];
    148     dCoeff [i] = (fabs (di) < COUENNE_EPS) ?
    149       0. :
    150       eigval / (di * di);
    151   }
    152 }
  • branches/Couenne/Couenne/src/convex/operators/conv-exprAbs.cpp

    r975 r1017  
    1717// generate convexification cut for constraint w = |x|
    1818
    19 void exprAbs::generateCuts (exprAux *w, const OsiSolverInterface &si,
     19void exprAbs::generateCuts (expression *w, const OsiSolverInterface &si,
    2020                            OsiCuts &cs, const CouenneCutGenerator *cg,
    2121                            t_chg_bounds *chg, int wind,
  • branches/Couenne/Couenne/src/convex/operators/conv-exprDiv.cpp

    r861 r1017  
    99
    1010#include "CouenneTypes.hpp"
     11#include "expression.hpp"
     12#include "exprAux.hpp"
    1113#include "exprOp.hpp"
    1214#include "exprDiv.hpp"
     
    2628
    2729// generate convexification cut for constraint w = x/y
    28 void exprDiv::generateCuts (exprAux *w, const OsiSolverInterface &si,
     30void exprDiv::generateCuts (expression *w, const OsiSolverInterface &si,
    2931                            OsiCuts &cs, const CouenneCutGenerator *cg,
    3032                            t_chg_bounds *chg, int wind,
  • branches/Couenne/Couenne/src/convex/operators/conv-exprExp.cpp

    r861 r1017  
    2020// generate convexification cut for constraint w = this
    2121
    22 void exprExp::generateCuts (exprAux *aux, const OsiSolverInterface &si,
     22void exprExp::generateCuts (expression *aux, const OsiSolverInterface &si,
    2323                            OsiCuts &cs,  const CouenneCutGenerator *cg,
    2424                            t_chg_bounds *chg, int wind,
  • branches/Couenne/Couenne/src/convex/operators/conv-exprGroup.cpp

    r975 r1017  
    2323  expression *lbnl, *ubnl;
    2424
     25  // TODO: do not aggregate members of exprSum
     26
    2527  // compute lower/upper bound of nonlinear part
    2628  exprSum::getBounds (lbnl, ubnl);
    2729
    2830  // count linear and constant terms
    29   int nlin = 0;
     31  int nlin = lcoeff_.size();
    3032  if (fabs (c0_) > COUENNE_EPS) nlin++;
    31   for (register int *ind = index_; *ind++>=0; nlin++);
     33  //  for (register int *ind = index_; *ind++>=0; nlin++);
    3234
    3335  expression
     
    4143  }
    4244
    43   for (register int *ind = index_, i=0; *ind>=0;) {
     45  // derive linear part (obtain constant)
     46  for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el) {
     47    //    c0 += el -> second;
    4448
    45     CouNumber coeff = coeff_ [i++];
     49  /*  // derive quadratic part (obtain linear part)
     50  for (sparseQ::iterator row = q_.begin (); row != q_.end (); ++row) {
    4651
    47     expression *l = new exprLowerBound (*ind),
    48                *u = new exprUpperBound (*ind++);
     52    int xind = row -> first -> Index ();
     53
     54    for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
     55  */
     56
     57    //  for (register int *ind = index_, i=0; *ind>=0;) {
     58
     59    CouNumber coeff = el -> second;//coeff_ [i++];
     60    int         ind = el -> first -> Index ();
     61
     62    expression *l = new exprLowerBound (ind),
     63               *u = new exprUpperBound (ind);
    4964
    5065    if (fabs (coeff - 1.) < COUENNE_EPS) {
     
    7590
    7691// generate equality between *this and *w
    77 void exprGroup::generateCuts (exprAux *w, const OsiSolverInterface &si,
     92void exprGroup::generateCuts (expression *w, const OsiSolverInterface &si,
    7893                              OsiCuts &cs, const CouenneCutGenerator *cg,
    7994                              t_chg_bounds *chg,
     
    87102
    88103  // there is one linear term so far: -w
    89   int nterms = 0;
     104  int nterms = lcoeff_.size ();
    90105
    91106  OsiRowCut *cut = new OsiRowCut;
    92107
    93108  // count terms in linear part
    94   for (register int *ind = index_; *ind++ >= 0; nterms++);
     109  //  for (register int *ind = index_; *ind++ >= 0; nterms++);
    95110
    96111  int displacement = (wind < 0) ? 1: 0;
     
    110125
    111126  // now add linear terms
    112   for (register int i=0; i<nterms; i++) {
     127  lincoeff::iterator el = lcoeff_.begin ();
     128  for (int i=0; el != lcoeff_.end (); ++el) {
     129    //  for (register int i=0; i<nterms; i++) {
    113130
    114     coeff [i + displacement] = coeff_ [i];
    115     index [i + displacement] = index_ [i];
     131    coeff [i   + displacement] = el -> second;
     132    index [i++ + displacement] = el -> first -> Index ();
    116133  }
    117134
    118135  // scan arglist for (aux) variables and constants
    119   for (register int i=0; i<nargs_; i++) {
     136  for (int i=0; i<nargs_; i++) {
    120137
    121138    expression *curr = arglist_ [i];
  • branches/Couenne/Couenne/src/convex/operators/conv-exprInv.cpp

    r915 r1017  
    5151// generate convexification cut for constraint w = 1/x
    5252
    53 void exprInv::generateCuts (exprAux *aux, const OsiSolverInterface &si,
     53void exprInv::generateCuts (expression *aux, const OsiSolverInterface &si,
    5454                            OsiCuts &cs, const CouenneCutGenerator *cg,
    5555                            t_chg_bounds *chg, int wind,
  • branches/Couenne/Couenne/src/convex/operators/conv-exprLog.cpp

    r959 r1017  
    2222// generate convexification cut for constraint w = this
    2323
    24 void exprLog::generateCuts (exprAux *aux, const OsiSolverInterface &si,
     24void exprLog::generateCuts (expression *aux, const OsiSolverInterface &si,
    2525                            OsiCuts &cs, const CouenneCutGenerator *cg,
    2626                            t_chg_bounds *chg, int wind,
  • branches/Couenne/Couenne/src/convex/operators/conv-exprMul-genCuts.cpp

    r988 r1017  
    99
    1010#include "CouenneTypes.hpp"
     11//#include "expression.hpp"
     12//#include "exprAux.hpp"
    1113#include "exprMul.hpp"
    12 #include "exprPow.hpp"
    13 #include "exprDiv.hpp"
    14 #include "CouenneProblem.hpp"
     14//#include "exprPow.hpp"
     15//#include "exprDiv.hpp"
     16//#include "CouenneProblem.hpp"
    1517#include "CouenneCutGenerator.hpp"
    1618
     
    1820/// generate convexification cut for constraint w = x*y
    1921
    20 void exprMul::generateCuts (exprAux *w, const OsiSolverInterface &si,
     22void exprMul::generateCuts (expression *w, const OsiSolverInterface &si,
    2123                            OsiCuts &cs, const CouenneCutGenerator *cg,
    2224                            t_chg_bounds *chg, int wind,
  • branches/Couenne/Couenne/src/convex/operators/conv-exprOpp.cpp

    r988 r1017  
    1717
    1818// generate equality between *this and *w
    19 void exprOpp::generateCuts (exprAux *w, const OsiSolverInterface &si,
     19void exprOpp::generateCuts (expression *w, const OsiSolverInterface &si,
    2020                            OsiCuts &cs, const CouenneCutGenerator *cg,
    2121                            t_chg_bounds *chg,
  • branches/Couenne/Couenne/src/convex/operators/conv-exprPow.cpp

    r988 r1017  
    6666// generate convexification cut for constraint w = x^k
    6767
    68 void exprPow::generateCuts (exprAux *aux, const OsiSolverInterface &si,
     68void exprPow::generateCuts (expression *aux, const OsiSolverInterface &si,
    6969                            OsiCuts &cs, const CouenneCutGenerator *cg,
    7070                            t_chg_bounds *chg, int wind,
  • branches/Couenne/Couenne/src/convex/operators/conv-exprQuad.cpp

    r988 r1017  
    1616#include "exprQuad.hpp"
    1717#include "exprBQuad.hpp"
    18 
     18#include "CouenneCutGenerator.hpp"
    1919
    2020/// Get lower and upper bound of an expression (if any)
     
    2727
    2828// generate equality between *this and *w
    29 void exprQuad::generateCuts (exprAux *w, const OsiSolverInterface &si,
     29void exprQuad::generateCuts (expression *w, const OsiSolverInterface &si,
    3030                             OsiCuts &cs, const CouenneCutGenerator *cg,
    3131                             t_chg_bounds *chg,
    3232                             int wind, CouNumber lb, CouNumber ub) {
    3333
    34   // check if we really need a convexification cut
    35   if (fabs ((*this) () - (*w) ()) < COUENNE_EPS)
     34  if ((!(cg -> isFirst ())) &&                    // unless a convexification was never created,
     35      (fabs ((*this) () - (*w) ()) < COUENNE_EPS) // do we really need a convexification cut?
     36      || !alphaConvexify (cg -> Problem (), si))  // ... or a new alpha-convexification?
    3637    return;
    37 
    38   // see if it is necessary to create/renew the alpha-convexification
    39   alphaConvexify (si);
    4038
    4139  // generate linear cuts for convex quadratic [upper|lower]-envelope
  • branches/Couenne/Couenne/src/convex/operators/conv-exprSinCos.cpp

    r979 r1017  
    2727/// convex cuts for sine or cosine
    2828int trigEnvelope (const CouenneCutGenerator *, OsiCuts &,
    29                    exprAux *, expression *, enum cou_trig);
     29                   expression *, expression *, enum cou_trig);
    3030
    3131
    3232/// generate convexification cut for constraint w = sin (this)
    3333
    34 void exprSin::generateCuts (exprAux *w, const OsiSolverInterface &si,
     34void exprSin::generateCuts (expression *w, const OsiSolverInterface &si,
    3535                            OsiCuts &cs, const CouenneCutGenerator *cg,
    3636                            t_chg_bounds *chg, int wind,
     
    5757/// generate convexification cut for constraint w = cos (this)
    5858
    59 void exprCos::generateCuts (exprAux *w, const OsiSolverInterface &si,
     59void exprCos::generateCuts (expression *w, const OsiSolverInterface &si,
    6060                            OsiCuts &cs, const CouenneCutGenerator *cg,
    6161                            t_chg_bounds *chg, int wind,
     
    8989int trigEnvelope (const CouenneCutGenerator *cg, // cut generator that has called us
    9090                   OsiCuts &cs,                  // cut set to be enriched
    91                    exprAux *w,
     91                   expression *w,
    9292                   expression *arg,
    9393                   enum cou_trig which_trig) {
  • branches/Couenne/Couenne/src/convex/operators/conv-exprSub.cpp

    r861 r1017  
    1414
    1515// generate equality between *this and *w
    16 void exprSub::generateCuts (exprAux *w, const OsiSolverInterface &si,
     16void exprSub::generateCuts (expression *w, const OsiSolverInterface &si,
    1717                            OsiCuts &cs, const CouenneCutGenerator *cg,
    1818                            t_chg_bounds *chg,
  • branches/Couenne/Couenne/src/convex/operators/conv-exprSum.cpp

    r861 r1017  
    1818
    1919// generate equality between *this and *w
    20 void exprSum::generateCuts (exprAux *w, const OsiSolverInterface &si,
    21                               OsiCuts &cs, const CouenneCutGenerator *cg,
    22                               t_chg_bounds *chg,
    23                               int wind, CouNumber lb, CouNumber ub) {
     20void exprSum::generateCuts (expression *w, const OsiSolverInterface &si,
     21                            OsiCuts &cs, const CouenneCutGenerator *cg,
     22                            t_chg_bounds *chg,
     23                            int wind, CouNumber lb, CouNumber ub) {
    2424  if (!(cg -> isFirst ()))
    2525    return;
  • branches/Couenne/Couenne/src/convex/operators/quadCuts.cpp

    r988 r1017  
    1717//#define DEBUG
    1818
    19 void exprQuad::quadCuts (exprAux *w, OsiCuts &cs, const CouenneCutGenerator *cg){
    20 
    21   assert (dIndex_ != NULL);
    22 
    23 #ifdef DEBUG
    24   std::cout<<"Expression has "<<nlterms_<<" linear terms and "
    25            <<nqterms_<<" quadratic terms. " << std::endl;
    26 
    27   printf ("Q\n");
    28   for (int i=0; i<nqterms_; i++)
    29     printf ("<%d,%d,%g>\n",  qindexI_ [i], qindexJ_ [i], qcoeff_ [i]);
    30 
    31   printf ("b\n");
    32   for (int i=0; i < nlterms_; i++)
    33     printf ("<%d,%g>\n",  index_ [i], coeff_ [i]);
     19
     20void exprQuad::quadCuts (expression *w, OsiCuts &cs, const CouenneCutGenerator *cg){
     21
     22#ifdef DEBUG
     23  std::cout<<"Expression has "<< lcoeff_.size () <<" linear terms and "
     24           << nqterms_ <<" quadratic terms." << std::endl;
     25
     26  printf ("Q:");
     27  for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
     28    int xind = row -> first -> Index ();
     29    for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col)
     30      printf (" <%d,%d,%g>",  xind, col -> first -> Index (), col -> second);
     31  }
     32
     33  printf ("\nb:");
     34  for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el)
     35    //for (int i=0; i < nlterms_; i++)
     36    printf (" <%d,%g>",  el -> first -> Index (), el -> second);//index_ [i], coeff_ [i]);
    3437
    3538  if (c0_)
    36     printf ("<c0 = %g>\n", c0_);
    37 
    38   if (dCoeffLo_ && dCoeffUp_ && dIndex_) {
    39     printf ("alpha\n");
    40     for (int i=0; i<nDiag_; i++)
    41       printf ("[%d,%g,%g]\n", dIndex_ [i], dCoeffLo_ [i], dCoeffUp_ [i]);
     39    printf ("; <c0 = %g>", c0_);
     40
     41  printf ("\nBounds:\n");
     42
     43  int index = 0;
     44
     45  for (std::map <exprVar *, std::pair <CouNumber, CouNumber> >::iterator i = bounds_.begin ();
     46       i != bounds_.end (); ++i, index++) {
     47
     48    printf ("%3d:\t", index);
     49    i -> first -> print (); printf ("\t");
     50    printf (" %13g [%13g, %13g]",
     51            (*(i -> first)) (), i -> second.first, i -> second.second);
     52
     53    CouNumber
     54      lb = cg -> Problem () -> Lb (i -> first -> Index ()),
     55      ub = cg -> Problem () -> Ub (i -> first -> Index ());
     56
     57    if ((eigen_.size () > 0) &&
     58        (fabs (ub-lb) > COUENNE_EPS))
     59      printf (" --> %20g, scaled %20g",
     60              eigen_.begin () -> first,// / (ub-lb),
     61              eigen_.begin () -> first / (ub-lb));
     62
     63    printf ("\n");
    4264  }
    4365#endif
     
    4971  // Get on which side constraint is violated to get the good lambda
    5072
    51   double
    52      varVal  = (*w)    (),
    53      exprVal = (*this) (),
    54     *lambda  = (varVal < exprVal) ?
    55       dCoeffLo_ : // Use under-estimator
    56       dCoeffUp_,  // Use  over-estimator
    57      convVal = 0;
     73  CouNumber
     74    varVal  = (*w)    (),
     75    exprVal = (*this) (),
     76    lambda  =
     77    (eigen_.size () == 0) ? 0 :
     78    (varVal < exprVal) ?
     79      CoinMin (0., eigen_.begin  () -> first) : // Use under-estimator
     80      CoinMax (0., eigen_.rbegin () -> first),  // Use  over-estimator
     81    convVal = 0;
    5882
    5983  const CouenneProblem& problem = *(cg -> Problem ());
     
    6892  // the current point
    6993
    70   if (dIndex_) {
     94  if (fabs (lambda) > COUENNE_EPS) {
    7195
    7296    convVal = exprVal;
    7397
    74     if (lambda) { // there is a convexification, check if out of
    75                   // current point
    76 
    77       for (int i=0; i<nDiag_; i++) {
    78 
    79         int ind = dIndex_ [i];
    80         CouNumber xi = colsol [ind];
    81         convVal += lambda [i] * (xi - lower [ind]) * (upper [ind] - xi);
    82       }
    83 
    84       if (varVal < exprVal) {if (convVal < varVal) return;}
    85       else                  {if (convVal > varVal) return;}
    86     }
    87   }
    88 
    89 #ifdef DEBUG
    90   std::cout << "Point to cut: "; for (int i = 0 ; i < numcols ; i++) std::cout << colsol [i] << ", ";
    91   printf (" (w,f(x),c) = (%g,%g,%g) -- lambda = [", (*w) (), exprVal, convVal);
    92   for (int i=0; i < nDiag_ ; i++) printf ("%g ", lambda [i]); printf ("]\n");
     98    // there is a convexification, check if out of current point
     99
     100    for (std::map <exprVar *, std::pair <CouNumber, CouNumber> >::iterator i = bounds_.begin ();
     101         i != bounds_.end (); ++i) {
     102
     103      int ind = i -> first -> Index ();
     104
     105      CouNumber
     106        xi = colsol [ind],
     107        lb = lower  [ind],
     108        ub = upper  [ind],
     109        delta = ub-lb;
     110
     111      if (fabs (delta) > COUENNE_EPS)
     112        convVal += lambda * (xi-lb) * (ub-xi) / (delta * delta);
     113    }
     114
     115    if (varVal < exprVal) {if (convVal < varVal) return;}
     116    else                  {if (convVal > varVal) return;}
     117  }
     118
     119#ifdef DEBUG
     120  std::cout << "Point to cut: ";
     121  for (int i = 0 ; i < numcols ; i++) std::cout << colsol [i] << ", ";
     122  printf (" (w,f(x),c) = (%g,%g,%g) -- lambda = %g\n", (*w) (), exprVal, convVal, lambda);
    93123#endif
    94124
     
    98128
    99129  // Compute 2 * Q x^*.
    100   for (int k = 0 ; k < nqterms_ ; k++) {
    101 
    102     int qi = qindexI_ [k],
    103         qj = qindexJ_ [k];
    104 
    105     CouNumber qc = qcoeff_ [k];
    106 
    107     if (qi != qj) {
    108       Qxs [qi] += 2 * qc * colsol [qj]; // contribution of element $q_{ij}$ to (Qx)_i
    109       Qxs [qj] += 2 * qc * colsol [qi]; //                         $q_{ij}$    (Qx)_j
    110     }
    111     // elements on the diagonal are not halved upon reading
    112     else Qxs [qi] += 2 * qc * colsol [qi];
     130  for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
     131
     132    int qi = row -> first -> Index ();
     133
     134    for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
     135
     136      int qj = col -> first -> Index ();//qindexJ_ [k];
     137
     138      CouNumber qc = 2 * col -> second;//qcoeff_ [k];
     139
     140      if (qi != qj) {
     141        Qxs [qi] += qc * colsol [qj]; // contribution of element $q_{ij}$ to (Qx)_i
     142        Qxs [qj] += qc * colsol [qi]; //                         $q_{ij}$    (Qx)_j
     143      }
     144      else Qxs [qi] += qc * colsol [qi]; // elements on the diagonal are not halved upon reading
     145    }
    113146  }
    114147
     
    118151
    119152  // Add a to it.
    120   for (int i = 0 ; i < nlterms_ ; i++)
    121     Qxs [index_ [i]] += coeff_ [i];
     153  for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el)
     154    Qxs [el -> first -> Index ()] += el -> second; //coeff_ [i];
    122155
    123156  // multiply Qx^* by x^*^T again and store the result for the lower
    124157  // bound into constant term
    125158
    126   double a0 = 0;//- c0_; // constant term
     159  double a0 = -c0_; // constant term
    127160
    128161  for (int i = 0; i < numcols; i++){
     
    140173  a0 -= exprVal;
    141174
    142   if (lambda != NULL) // Now the part which depends on lambda, if there is one
    143 
    144     for (int k = 0 ; k < nDiag_ ; k++) {
    145 
    146       int ind = dIndex_ [k];
    147       CouNumber xi = colsol [ind],
    148         coeff = lambda [k] * (lower  [ind] + upper  [ind] - 2 * xi);
    149 
    150       a0 += coeff * xi - lambda [k] * (xi - lower [ind]) * (upper [ind] - xi);
     175  if (fabs (lambda) > COUENNE_EPS) // Now the part which depends on lambda, if there is one
     176
     177    for (std::map <exprVar *, std::pair <CouNumber, CouNumber> >::iterator i = bounds_.begin ();
     178         i != bounds_.end (); ++i) {
     179
     180      int ind = i -> first -> Index ();
     181
     182      CouNumber
     183        xi    = colsol [ind],
     184        lb    = lower [ind],
     185        ub    = upper [ind],
     186        delta = ub-lb,
     187        coeff = lambda / (delta*delta) * (lb + ub - 2 * xi);
     188
     189      if (fabs (delta) > COUENNE_EPS)
     190        a0 += coeff * xi - lambda / (delta*delta) * (xi - lb) * (ub - xi);
    151191      //      a0 += lambda [k] * lower  [ind] * upper  [ind];
    152192      //      a0 -= lambda [k] * colsol [ind] * colsol [ind];
     
    163203
    164204#ifdef DEBUG
    165   printf ("2Qx = (");
    166   for(int i = 0; i < numcols; i++)
    167     printf ("%g ", Qxs [i]);
    168   printf (")[%g], %d nz\n", a0, nnz);
     205  printf ("2Qx = (");for(int i=0;i<numcols;i++)printf("%g ",Qxs[i]);printf (")[%g], %d nz\n",a0,nnz);
    169206#endif
    170207
     
    198235
    199236  delete [] Qxs;
    200 
    201   if (lambda == dCoeffLo_) {
     237   
     238  if (varVal < exprVal) { //(lambda == dCoeffLo_) {
     239
    202240     cut.setUb (a0);
    203241
     
    215253  }
    216254  else {
     255
    217256    cut.setLb (a0);
    218257#ifdef DEBUG
  • branches/Couenne/Couenne/src/expression/CouenneTypes.hpp

    r915 r1017  
    5555enum dig_type {ORIG_ONLY, STOP_AT_AUX, TAG_AND_RECURSIVE, COUNT};
    5656
    57 /** integrality type of an auxiliary variable: unset, continuous, integer */
    58 enum integer_type {AUX_UNSET=-1, AUX_CONTINUOUS, AUX_INTEGER};
    59 
    6057/** status of lower/upper bound of a variable, to be checked/modified
    6158    in bound tightening */
  • branches/Couenne/Couenne/src/expression/exprAux.cpp

    r1001 r1017  
    88 */
    99
    10 #include "CouenneCutGenerator.hpp"
    11 #include "CouenneTypes.hpp"
    12 #include "expression.hpp"
    1310#include "exprAux.hpp"
    14 #include "exprVar.hpp"
    15 #include "CouenneProblem.hpp"
     11#include "exprBound.hpp"
     12#include "exprMax.hpp"
     13#include "exprMin.hpp"
     14#include "CglCutGenerator.hpp"
     15
     16class CouenneCutGenerator;
    1617
    1718//#define DEBUG
    1819
    1920// auxiliary expression Constructor
    20 exprAux::exprAux (expression *image, int index, int rank, bool integer):
     21exprAux::exprAux (expression *image, int index, int rank, enum intType isInteger):
    2122
    2223  exprVar       (index),
     
    2425  rank_         (rank),
    2526  multiplicity_ (1),
    26   integer_      (integer) {
     27  integer_      (isInteger) {
    2728
    2829  // do this later, in standardize()
     
    5051  rank_         (-1),
    5152  multiplicity_ (0),
    52   integer_      (false) {}
    53 
     53  integer_      (Unset) {}
     54//(image -> isInteger () ? Integer : Continuous)
    5455
    5556/// Copy constructor
     
    119120
    120121/// I/O
    121 void exprAux::print (std::ostream &out, bool descend, CouenneProblem *p) const {
     122void exprAux::print (std::ostream &out, bool descend) const {
    122123
    123124  if (descend)
    124     image_ -> print (out, descend, p);
     125    image_ -> print (out, descend);
    125126  else {
    126127    if (integer_) out << "u_"; // TODO: should be isInteger instead of
     
    136137/// expression
    137138int exprAux::DepList (std::set <int> &deplist,
    138                       enum dig_type type,
    139                       const CouenneProblem *p) {
     139                      enum dig_type type) {
    140140
    141141  if (type == ORIG_ONLY)   
    142     return image_ -> DepList (deplist, type, p);
     142    return image_ -> DepList (deplist, type);
    143143
    144144  if (deplist.find (varIndex_) == deplist.end ())
     
    149149    return 1;
    150150
    151   return 1 + image_ -> DepList (deplist, type, p);
     151  return 1 + image_ -> DepList (deplist, type);
    152152}
    153153
     
    207207      while (n--) {
    208208        if (fabs (el [n]) > COU_MAX_COEFF)  {
    209           printf ("### Warning: coefficient too large: %g x%d [", el [n], ind [n]);
     209          printf ("Couenne, warning: coefficient too large %g x%d: ", el [n], ind [n]);
    210210          cut -> print ();
    211211          warned_large_coeff = true;
     
    214214
    215215        if (fabs (rhs) > COU_MAX_COEFF) {
    216           printf ("rhs too large %g: ", rhs);
     216          printf ("Couenne, warning: rhs too large (%g): ", rhs);
    217217          cut -> print ();
    218218          warned_large_coeff = true;
     
    226226      (ncc < cs.sizeColCuts ()))
    227227    {
    228       printf ("----------------Generated cut for ");
     228      printf ("---------------- ConvCut: ");
    229229      print (std::cout);  printf (" := ");
    230230      image_ -> print (std::cout);
  • branches/Couenne/Couenne/src/expression/exprAux.hpp

    r1001 r1017  
    1313
    1414#include <iostream>
     15
     16#include "expression.hpp"
    1517#include "exprVar.hpp"
    16 #include "exprBound.hpp"
    17 #include "exprMax.hpp"
    18 #include "exprMin.hpp"
    19 #include "CouenneTypes.hpp"
    20 #include "CglCutGenerator.hpp"
    2118
     19class CouenneCutGenerator;
    2220
    2321/** Auxiliary variable
     
    3129
    3230class exprAux: public exprVar {
     31
     32 public:
     33
     34  /// integrality type of an auxiliary variable: unset, continuous, integer
     35  enum intType {Unset=-1, Continuous, Integer};
    3336
    3437 protected:
     
    5659
    5760  /// is this variable integer?
    58   bool integer_;
     61  enum intType integer_;
    5962
    6063 public:
     
    6568
    6669  /// Constructor
    67   exprAux (expression *, int, int, bool = false);
     70  exprAux (expression *, int, int, intType = Unset);
    6871
    6972  /// Constructor to be used with standardize ([...], false)
     
    8588  /// Print expression
    8689  virtual void print (std::ostream & = std::cout,
    87                       bool = false, CouenneProblem * = NULL) const;
     90                      bool = false) const;
    8891
    8992  /// The expression associated with this auxiliary variable
     
    102105  /// expression
    103106  int DepList (std::set <int> &deplist,
    104                enum dig_type type = ORIG_ONLY,
    105                const CouenneProblem *p = NULL);
     107               enum dig_type type = ORIG_ONLY);
    106108
    107109  /// simplify
     
    117119
    118120  /// Get lower and upper bound of an expression (if any) -- real values
    119   void getBounds (CouNumber &lb, CouNumber &ub)
    120     {expression::getBounds (lb, ub);}
     121  //void getBounds (CouNumber &lb, CouNumber &ub)
     122  //{expression::getBounds (lb, ub);}
    121123
    122124  /// set bounds depending on both branching rules and propagated
     
    132134
    133135  /// used in rank-based branching variable choice
    134   virtual inline int rank (CouenneProblem *p = NULL)
     136  virtual inline int rank ()
    135137    {return rank_;}
    136138
    137139  /// is this expression integer?
    138140  virtual inline bool isInteger () {
    139     return ((integer_ == AUX_INTEGER) ||
    140             (integer_ == AUX_UNSET) && ((integer_ = image_ -> isInteger ()) == AUX_INTEGER));
     141
     142    if ((integer_ == Integer) ||
     143        (integer_ == Unset) &&
     144        ((integer_ = (image_ -> isInteger ()) ?
     145          Integer : Continuous) == Integer))
     146      return true;
     147
     148    CouNumber lb = (*(Lb ())) ();
     149    return (::isInteger (lb) && (fabs (lb - (*(Ub ())) ()) < COUENNE_EPS));
    141150  }
    142151
  • branches/Couenne/Couenne/src/expression/exprBound.hpp

    r988 r1017  
    5454  /// Print to iostream
    5555  void print (std::ostream &out = std::cout,
    56               bool = false, CouenneProblem * = NULL) const
     56              bool = false) const
    5757    {out << "l_" << varIndex_;}
    5858
     
    6666
    6767  /// dependence on variable set
    68   inline int dependsOn (int *, int)
     68  inline int dependsOn (int *, int, enum dig_type type = STOP_AT_AUX)
    6969    {return 0;}
    7070
     
    103103  /// Print to iostream
    104104  void print (std::ostream &out = std::cout,
    105               bool = false, CouenneProblem * = NULL) const
     105              bool = false) const
    106106    {out << "u_" << varIndex_;}
    107107
     
    115115
    116116  /// dependence on variable set
    117   inline int dependsOn (int *, int)
     117  inline int dependsOn (int *, int, enum dig_type type = STOP_AT_AUX)
    118118    {return 0;}
    119119
  • branches/Couenne/Couenne/src/expression/exprClone.hpp

    r973 r1017  
    4242  /// I/O
    4343  void print (std::ostream &out = std::cout,
    44               bool descend      = false,
    45               CouenneProblem *p = NULL) const
    46     {copy_ -> Original () -> print (out, descend, p);}
     44              bool descend      = false) const
     45    {copy_ -> Original () -> print (out, descend);}
    4746
    4847  /// value
  • branches/Couenne/Couenne/src/expression/exprConst.hpp

    r960 r1017  
    5151  /// I/O
    5252  void print (std::ostream &out = std::cout,
    53               bool = false, CouenneProblem * = NULL) const
     53              bool = false) const
    5454    {out << value_;}
    5555
     
    6363
    6464  /// dependence on variable set
    65   int dependsOn (int *ind, int n,
    66                  CouenneProblem *p = NULL,
    67                  enum dig_type   type = STOP_AT_AUX)
     65  int dependsOn (int *ind, int n, enum dig_type type = STOP_AT_AUX)
    6866    {return 0;}
    6967
     
    7977
    8078  /// generate convexification cut for constraint w = this
    81   void generateCuts (exprAux *, const OsiSolverInterface &,
     79  void generateCuts (expression *, const OsiSolverInterface &,
    8280                     OsiCuts &, const CouenneCutGenerator *,
    8381                     t_chg_bounds * = NULL, int = -1,
     
    9189  /// is this expression integer?
    9290  virtual bool isInteger ()
    93   {return (fabs (value_ - COUENNE_round (value_)) < COUENNE_EPS);}
     91  {return ::isInteger (value_);}
    9492
    9593  /// used in rank-based branching variable choice
    96   virtual int rank (CouenneProblem *p)
     94  virtual int rank ()
    9795    {return 0;}
    9896};
  • branches/Couenne/Couenne/src/expression/exprCopy.hpp

    r1001 r1017  
    8585  /// I/O
    8686  virtual void print (std::ostream &out = std::cout,
    87                       bool descend      = false,
    88                       CouenneProblem *p = NULL) const
    89     {copy_ -> Original () -> print (out, descend, p);}
     87                      bool descend      = false) const
     88    {copy_ -> Original () -> print (out, descend);}
    9089
    9190  /// value
     
    111110  /// expression
    112111  inline int DepList (std::set <int> &deplist,
    113                       enum dig_type   type = ORIG_ONLY,
    114                       const CouenneProblem *p    = NULL)
    115     {return copy_ -> DepList (deplist, type, p);}
     112                      enum dig_type   type = ORIG_ONLY)
     113    {return copy_ -> DepList (deplist, type);}
    116114
    117115  /// simplify expression (useful for derivatives)
     
    135133
    136134  /// generate convexification cut for constraint w = this
    137   inline void generateCuts (exprAux *w, const OsiSolverInterface &si,
     135  inline void generateCuts (expression *w, const OsiSolverInterface &si,
    138136                            OsiCuts &cs, const CouenneCutGenerator *cg,
    139137                            t_chg_bounds *chg = NULL, int wind= -1,
     
    161159
    162160  /// used in rank-based branching variable choice
    163   int rank (CouenneProblem *p)
    164     {return copy_ -> rank (p);}
     161  int rank ()
     162    {return copy_ -> rank ();}
    165163
    166164  /// implied bound processing
     
    175173  /// Set up branching object by evaluating many branching points for each expression's arguments.
    176174  /// Return estimated improvement in objective function
    177   virtual CouNumber selectBranch (const CouenneObject *obj,
     175  virtual CouNumber selectBranch (const CouenneObject *obj, 
    178176                                  const OsiBranchingInformation *info,
    179                                   int     &ind,
    180                                   double *&brpts,
    181                                   int     &way)
    182   {return copy_ -> selectBranch (obj, info, ind, brpts, way);}
     177                                  expression * &var,
     178                                  double * &brpts,
     179                                  int &way)
     180
     181  {return copy_ -> selectBranch (obj, info, var, brpts, way);}
    183182
    184183  /// replace occurrence of a variable with another variable
    185   void replace (exprVar *, exprVar *);
     184  virtual void replace (exprVar *, exprVar *);
    186185
    187186  /// fill in dependence structure
  • branches/Couenne/Couenne/src/expression/exprIVar.hpp

    r988 r1017  
    1717#include "exprConst.hpp"
    1818#include "exprVar.hpp"
    19 
    20 class CouenneProblem;
    2119
    2220
     
    4139
    4240  /// print
    43   virtual void print (std::ostream &out = std::cout, bool = false,
    44                       CouenneProblem * = NULL) const
     41  virtual void print (std::ostream &out = std::cout, bool = false) const
    4542    {out << "y_" << varIndex_;}
    4643
  • branches/Couenne/Couenne/src/expression/exprOp.cpp

    r988 r1017  
    88 */
    99
    10 #include "CouenneProblem.hpp"
    11 
    1210#include "expression.hpp"
    1311#include "exprAux.hpp"
     
    1614#include "exprQuad.hpp"
    1715
     16class CouenneProblem;
    1817
    1918// General N-ary function destructor
     
    4241
    4342void exprOp::print (std::ostream &out,
    44                     bool descend,
    45                     CouenneProblem *p) const {
     43                    bool descend) const {
    4644 
    4745  if (printPos () == PRE)
     
    5149  for (int i=0; i<nargs_; i++) {
    5250    if (arglist_ [i])
    53       arglist_ [i] -> print (out, descend, p);
     51      arglist_ [i] -> print (out, descend);
    5452    fflush (stdout);
    5553    if (i < nargs_ - 1) {
     
    6664/// compare general n-ary expressions
    6765
    68 int exprOp::compare (exprOp  &e1) {
     66int exprOp::compare (exprOp &e1) {
    6967
    70   int c0 = code (),
     68  int c0 =     code (),
    7169      c1 = e1. code ();
    7270
     
    113111/// used in rank-based branching variable choice
    114112
    115 int exprOp::rank (CouenneProblem *p) {
     113int exprOp::rank () {
    116114
    117115  int maxrank = -1;
     
    119117  for (register expression **al = arglist_ + nargs_;
    120118       al-- > arglist_;) {
    121     register int r = (*al) -> rank (p);
     119    register int r = (*al) -> rank ();
    122120    if (r > maxrank) maxrank = r;
    123121  }
     
    137135// and the like) will do the part for its own object
    138136
    139 exprAux *exprOp::standardize (register CouenneProblem *p, bool addAux) {
     137exprAux *exprOp::standardize (CouenneProblem *p, bool addAux) {
    140138
    141139  register exprVar *subst;
     
    157155  for (register int i = nargs_; i--; al++)
    158156
    159     switch ((*al) -> Type () == VAR) {
     157    switch ((*al) -> Type ()) {
    160158
    161159    case VAR:
     
    182180  for (int i = nargs_; i--;)
    183181
    184     if (!(arglist_ [i] -> isInteger ())) { // this argument is not integer
     182    if (!(arglist_ [i] -> isInteger ())) {
    185183
    186       // last chance: check if constant and integer
     184      // this argument is not integer: check if constant and integer
    187185
    188       expression *lb, *ub;
     186      CouNumber lb, ub;
     187      arglist_ [i] -> getBounds (lb, ub);
    189188
    190       arglist_ [i] -> getBounds (lb, ub);
    191       CouNumber lv = (*lb) ();
    192 
    193       if ((fabs (lv - (*ub) ()) > COUENNE_EPS) ||
    194           (fabs (COUENNE_round (lv) - lv) < COUENNE_EPS))
    195       return false;
     189      if ((fabs (lb - ub) > COUENNE_EPS) ||
     190          !::isInteger (lb))
     191        return false;
    196192    }
    197193
  • branches/Couenne/Couenne/src/expression/exprOp.hpp

    r1001 r1017  
    1515#include "expression.hpp"
    1616#include "CouenneTypes.hpp"
    17 #include "exprUnary.hpp"
    1817
     18class CouenneProblem;
    1919
    2020/// general n-ary operator-type expression: requires argument
     
    2727 protected:
    2828
    29   expression **arglist_; //< argument list is an array of pointers to other expressions
    30   int          nargs_;   //< number of arguments (cardinality of arglist)
     29  expression **arglist_; ///< argument list is an array of pointers to other expressions
     30  int          nargs_;   ///< number of arguments (cardinality of arglist)
    3131
    3232 public:
     
    7171  /// I/O
    7272  virtual void print (std::ostream &out = std::cout,
    73                       bool = false, CouenneProblem * = NULL) const;
     73                      bool = false) const;
    7474
    7575  /// print position (PRE, INSIDE, POST)
     
    8181    {return "??";}
    8282
    83   /// function for the evaluation of the expression
    84   //  virtual inline CouNumber operator () ();
    85 
    86   /// dependence on variable set
    87   //  virtual int dependsOn (int * = NULL, int = 1);
    88 
    8983  /// fill in the set with all indices of variables appearing in the
    9084  /// expression
    9185  virtual inline int DepList (std::set <int> &deplist,
    92                               enum dig_type type = ORIG_ONLY,
    93                               const CouenneProblem *p = NULL) {
     86                              enum dig_type type = ORIG_ONLY) {
    9487    int tot = 0;
    9588    for (int i = nargs_; i--;)
    96       tot += arglist_ [i] -> DepList (deplist, type, p);
     89      tot += arglist_ [i] -> DepList (deplist, type);
    9790    return tot;
    9891  }
     
    137130
    138131  /// used in rank-based branching variable choice
    139   virtual int rank (CouenneProblem *);
     132  virtual int rank ();
    140133
    141134  /// fill in dependence structure
     
    150143};
    151144
    152 
    153 /// expression evaluation -- n-ary operator (non-variable, non-constant)
    154 
    155 /*
    156 inline CouNumber exprOp::operator () () {
    157 
    158   /// Fetch argument list and compute it "recursively" (the operator()
    159   /// of the elements in the list is called) to fill in the vector
    160   /// containing the numerical value of the argument list.
    161 
    162   register expression **al = arglist_;
    163 
    164   for (register int i = nargs_; i--;)
    165     *++sp = (**al++) ();
    166 
    167   return 0;
    168 }
    169 */
    170 
    171145#endif
  • branches/Couenne/Couenne/src/expression/exprStore.hpp

    r973 r1017  
    4545    {return new exprStore (*this);}
    4646
    47   /// value (the saved one)
    48   //  virtual inline CouNumber Value () const
    49   //    {return value_;}
    50 
    5147  /// null function for evaluating the expression
    5248  virtual inline CouNumber operator () ()
  • branches/Couenne/Couenne/src/expression/exprUnary.cpp

    r988 r1017  
    1616// print unary expression
    1717void exprUnary::print (std::ostream &out,
    18                        bool descend,
    19                        CouenneProblem *p) const {
     18                       bool descend) const {
    2019
    2120  if (printPos () == PRE)  out << printOp ();
    2221  out << "(";
    23   argument_ -> print (out, descend, p);
     22  argument_ -> print (out, descend);
    2423  out << ")";
    2524  if (printPos () == POST) out << printOp ();
     
    7574  // the corresponding evaluated expression is integer.
    7675
    77   register expression *al, *au;
     76  CouNumber al, au;
     77  argument_ -> getBounds (al, au);
    7878
    79   argument_ -> getBounds (al, au);
    80   CouNumber val = (*al) ();
     79  if (fabs (al - au) < COUENNE_EPS) { // argument is constant
    8180
    82   if (fabs (val - (*au) ()) < COUENNE_EPS) {
    83     // argument is constant
    84 
    85     register CouNumber fval = (F ()) (val);
     81    register CouNumber fval = (F ()) (al);
    8682
    8783    // check if f(lb=ub) is integer
  • branches/Couenne/Couenne/src/expression/exprUnary.hpp

    r1001 r1017  
    1515#include "expression.hpp"
    1616#include "CouenneTypes.hpp"
    17 #include "exprOp.hpp"
    1817
    1918
     
    7473  /// print this expression to iostream
    7574  virtual void print (std::ostream &out = std::cout,
    76                       bool = false, CouenneProblem * = NULL) const;
     75                      bool = false) const;
    7776
    7877  /// print position (PRE, INSIDE, POST)
     
    9190  /// expression
    9291  virtual inline int DepList (std::set <int> &deplist,
    93                               enum dig_type type = ORIG_ONLY,
    94                               const CouenneProblem *p = NULL)
    95     {return argument_ -> DepList (deplist, type, p);}
     92                              enum dig_type type = ORIG_ONLY)
     93    {return argument_ -> DepList (deplist, type);}
    9694
    9795  /// simplification
     
    123121
    124122  /// used in rank-based branching variable choice
    125   virtual int rank (CouenneProblem *p)
    126     {return (argument_ -> rank (p));}
     123  virtual int rank ()
     124    {return (argument_ -> rank ());}
    127125
    128126  /// fill in dependence structure
  • branches/Couenne/Couenne/src/expression/exprVar.cpp

    r988 r1017  
    99
    1010#include "CouenneCutGenerator.hpp"
    11 #include "CouenneTypes.hpp"
    12 #include "expression.hpp"
    1311#include "exprAux.hpp"
    14 #include "exprOp.hpp"
    15 #include "exprUnary.hpp"
    1612#include "exprVar.hpp"
    1713#include "exprBound.hpp"
    18 
    19 #include "CouenneProblem.hpp"
    2014#include "depGraph.hpp"
    2115
     
    3024
    3125// generate convexification cut for constraint w = this
    32 void exprVar::generateCuts (exprAux *w, const OsiSolverInterface &si,
     26void exprVar::generateCuts (expression *w, const OsiSolverInterface &si,
    3327                            OsiCuts &cs, const CouenneCutGenerator *cg,
    3428                            t_chg_bounds *chg, int,
     
    4438
    4539  bool res = false;
    46   if (updateBound (-1, l + varIndex_, l [wind])) {res = true; chg [varIndex_].setLower(t_chg_bounds::CHANGED);}
    47   if (updateBound (+1, u + varIndex_, u [wind])) {res = true; chg [varIndex_].setUpper(t_chg_bounds::CHANGED);}
     40
     41  if (updateBound (-1, l + varIndex_, l [wind]))
     42    {res = true; chg [varIndex_].setLower(t_chg_bounds::CHANGED);}
     43  if (updateBound (+1, u + varIndex_, u [wind]))
     44    {res = true; chg [varIndex_].setUpper(t_chg_bounds::CHANGED);}
     45
    4846  return res;
    4947}
     
    5553
    5654
    57 /// Bound get
    58 expression *exprVar::Lb () {return new exprLowerBound (varIndex_);}
    59 expression *exprVar::Ub () {return new exprUpperBound (varIndex_);}
     55expression *exprVar::Lb () {return new exprLowerBound (varIndex_);}///< lower bound of a variable
     56expression *exprVar::Ub () {return new exprUpperBound (varIndex_);}///< upper bound of a variable
  • branches/Couenne/Couenne/src/expression/exprVar.hpp

    r1001 r1017  
    1717#include "expression.hpp"
    1818#include "exprConst.hpp"
    19 
    20 class CouenneProblem;
    2119
    2220
     
    5856
    5957  /// for compatibility with exprAux
    60   virtual inline expression *Image () const
    61     {return NULL;}
     58  //virtual inline expression *Image () const
     59  //{return NULL;}
    6260
    63   // Bound get
     61  // Bounds
    6462  virtual expression *Lb (); ///< get lower bound expression
    6563  virtual expression *Ub (); ///< get upper bound expression
    6664
     65  // Bounds
     66  virtual CouNumber lb () {return expression::lbounds_ [varIndex_];} ///< lower bound
     67  virtual CouNumber ub () {return expression::ubounds_ [varIndex_];} ///< upper bound
     68
    6769  /// print
    6870  virtual void print (std::ostream &out = std::cout,
    69                       bool = false, CouenneProblem * = NULL) const
     71                      bool = false) const
    7072    {out << "x_" << varIndex_;}
    7173
     
    7375  virtual inline CouNumber operator () ()
    7476    {return expression::variables_ [varIndex_];}
    75 
    76   /// return the value of the variable
    77   //  inline CouNumber Value ()
    78   //    {return currValue_;}
    7977
    8078  /// differentiation
     
    8583  /// expression
    8684  virtual inline int DepList (std::set <int> &deplist,
    87                               enum dig_type type = ORIG_ONLY,
    88                               const CouenneProblem * = NULL) {
     85                              enum dig_type type = ORIG_ONLY) {
    8986
    9087    if (deplist.find (varIndex_) == deplist.end ()) {
     
    107104    {return LINEAR;}
    108105
    109   /// is this expression integer?
    110   virtual inline bool isInteger ()
    111     {return false;}
     106  /// is this variable integer?
     107  virtual inline bool isInteger () {
     108    CouNumber lb = (*(Lb ())) ();
     109    return ((fabs (lb - (*(Ub ())) ()) < COUENNE_EPS) && (::isInteger (lb)));
     110  }
    112111
    113112  /// Get lower and upper bound of an expression (if any)
     
    122121
    123122  /// generate convexification cut for constraint w = this
    124   virtual void generateCuts (exprAux *w, const OsiSolverInterface &si,
     123  virtual void generateCuts (expression *w,
     124                             const OsiSolverInterface &si,
    125125                             OsiCuts &cs, const CouenneCutGenerator *cg,
    126126                             t_chg_bounds * = NULL, int = -1,
     
    139139
    140140  /// rank of an original variable is always one
    141   virtual int rank (CouenneProblem *p)
     141  virtual int rank ()
    142142    {return 1;}
    143143
    144144  /// update dependence set with index of this variable
    145145  virtual void fillDepSet (std::set <DepNode *, compNode> *, DepGraph *);
     146
     147  /// is this variable fixed?
     148  virtual inline bool isFixed ()
     149  {return (fabs ((*(Lb ())) () - (*(Ub ())) ()) < COUENNE_EPS);}
    146150};
    147151
  • branches/Couenne/Couenne/src/expression/expression.cpp

    r973 r1017  
    99
    1010#include "CouenneCutGenerator.hpp"
    11 #include "CouenneProblem.hpp"
    1211
    1312#include "CouenneTypes.hpp"
     
    5049// generate one cut for a constant
    5150
    52 void exprConst::generateCuts (exprAux *w, const OsiSolverInterface &si,
     51void exprConst::generateCuts (expression *w, const OsiSolverInterface &si,
    5352                              OsiCuts &cs, const CouenneCutGenerator *cg,
    5453                              t_chg_bounds *chg, int,
     
    132131/// dependence on variable set: return cardinality of subset of the
    133132/// set of indices in first argument which occur in expression.
    134 int expression::dependsOn (int *ind, int n,
    135                            CouenneProblem *p,
    136                            enum dig_type type) {
     133int expression::dependsOn (int *ind, int n, enum dig_type type) {
    137134
    138135  std::set <int>
     
    141138    intersectn;
    142139
    143   DepList (deplist, type, p);
     140  DepList (deplist, type);
    144141
    145142  std::set_intersection (indlist .begin (), indlist .end (),
  • branches/Couenne/Couenne/src/expression/expression.hpp

    r1003 r1017  
    120120    {return EMPTY;}
    121121
     122  /// return pointer to self
     123  virtual inline expression *Image () const
     124    {return NULL;}
     125
    122126  /// value (empty)
    123127  virtual inline CouNumber Value () const
     
    132136  /// print expression to iostream
    133137  virtual void print (std::ostream &s  = std::cout,   /// output stream
    134                       bool             = false,       /// descend into auxiliaries' image?
    135                       CouenneProblem * = NULL) const  /// problem pointer (in exprGroup)
     138                      bool             = false) const /// descend into auxiliaries' image?
    136139    {s << '?';}
    137140
     
    147150  /// set of indices in first argument which occur in expression.
    148151  virtual int dependsOn (int *ind, int n,
    149                          CouenneProblem *p = NULL,
    150                          enum dig_type   type = STOP_AT_AUX);
     152                         enum dig_type type = STOP_AT_AUX);
     153
     154  /// version with one index only
     155  inline int dependsOn (int singleton,
     156                         enum dig_type type = STOP_AT_AUX)
     157  {return dependsOn (&singleton, 1, type);}
    151158
    152159  /// fill std::set with indices of variables on which this expression
     
    154161  /// pointers (exprGroup, exprQuad)
    155162  virtual inline int DepList (std::set <int> &deplist,
    156                               enum dig_type   type = ORIG_ONLY,
    157                               const CouenneProblem *p    = NULL)
     163                              enum dig_type   type = ORIG_ONLY)
    158164    {return 0;}
    159165
     
    192198
    193199  /// generate convexification cut for constraint w = this
    194   virtual void generateCuts (exprAux *w, const OsiSolverInterface &si,
     200  virtual void generateCuts (expression *w, const OsiSolverInterface &si,
    195201                             OsiCuts &cs, const CouenneCutGenerator *cg,
    196202                             t_chg_bounds *chg = NULL, int wind = -1,
     
    221227  /// have rank 1; auxiliary w=f(x) has rank r(w) = r(x)+1; finally,
    222228  /// auxiliary w=f(x1,x2...,xk) has rank r(w) = 1+max{r(xi):i=1..k}.
    223   virtual int rank (CouenneProblem *p = NULL)
    224     {return -1;} // return null rank
     229  virtual int rank ()
     230  {return -1;} // return null rank
    225231
    226232  /// does a backward implied bound processing on every expression,
     
    231237  /// bound on the variables on which the expression depends.
    232238  virtual bool impliedBound (int, CouNumber *, CouNumber *, t_chg_bounds *)
    233     {return false;}
     239  {return false;}
    234240
    235241  /// multiplicity of a variable
    236242  virtual int Multiplicity ()
    237     {return 1;}
     243  {return 1;}
    238244
    239245  /// set up branching object by evaluating many branching points for
     
    242248  virtual CouNumber selectBranch (const CouenneObject *obj,
    243249                                  const OsiBranchingInformation *info,
    244                                   int &ind,
     250                                  expression * &var,
    245251                                  double * &brpts,
    246252                                  int &way)
    247     {ind = -1; return 0.;}
     253  {var = NULL; return 0.;}
    248254
    249255  /// replace expression with another
     
    283289}
    284290
     291/// is this number integer?
     292inline bool isInteger (CouNumber x)
     293{return (fabs (COUENNE_round (x) - x) < COUENNE_EPS);}
     294
    285295#endif
  • branches/Couenne/Couenne/src/expression/operators/bounds/exprBQuad.cpp

    r854 r1017  
    1313//#define DEBUG
    1414
    15 CouNumber computeQBound (int sign, exprQuad *e) {
     15CouNumber exprQuad::computeQBound (int sign) {
    1616
    17   // compute lower (if sign == -1) or upper (sign == +1) bound of an
     17  // 1) loose: disaggregated bound
     18  // 2) tighter: aggregate coefficient per variable
     19  // 3) tightest: solve (convex) QP on available alpha-convexification
     20
     21  // 1) loose bound
     22  //
     23  // w = a0 + a'x + x'Qx means that its lower bound is
     24  //
     25  // w_l = a0 + a'_+ x_l + a'_- x_u + \sum_
     26
     27
     28  // 2) tighter bound -- TODO
     29  //
     30  // Compute lower (if sign == -1) or upper (sign == +1) bound of an
    1831  // exprQuad based on the information obtained through
    1932  // alpha-convexification, if any, or as follows:
    2033  //
    2134  // w = a0 + a'x + x'Qx =
    22   //   = a0 + sum{i} [(a_i + sum {j} q_ij * x_j) * x_i] =
    23   //   = a0 + sum{i} [                       z_i * x_i]
     35  //   = a0 + sum{i} [(a_i + sum {j>=i} q_ij * x_j) * x_i] =
     36  //   = a0 + sum{i} [                          z_i * x_i]
    2437  //
    25   // So some bound on z_i can be computed and a bound on the whole
     38  // Thus, some bound on z_i can be computed and a bound on the whole
    2639  // expression should be better than what can be obtained by summing
    2740  // all bounds separately.
     
    3043  // the convexification after some updates in the variable bounds
    3144  // without updating the convexification. Notice also that the
    32   // direction can also be vertical, not only horizontal
    33 
    34   int
    35     nlt = e -> getnLTerms (),
    36     *li = e -> getIndices (),
    37 
    38     nqt = e -> getnQTerms (),
    39     *qi = e -> getQIndexI (),
    40     *qj = e -> getQIndexJ ();
     45  // direction can also be vertical, not only horizontal.
    4146
    4247  CouNumber
    43     *lc = e -> getCoeffs  (),
    44     *qc = e -> getQCoeffs (),
     48    bound = c0_,
    4549    *lb = expression::Lbounds (),
    4650    *ub = expression::Ubounds (),
    47     bound = e -> getc0 (),
    4851    term;
    4952
    50 #ifdef DEBUG
    51   printf ("\n");
    52   for (int i=0; i<12; i++) printf ("%3d [%g,%g]\n",i, lb [i], ub [i]);
    53   e -> print ();
    54   printf ("\n (%g)\n ", bound);
    55 #endif
     53  // derive linear part (obtain constant)
     54  for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el) {
    5655
    57   if (sign < 0) { // compute lower bound ////////////////////////////////////////////////
     56    CouNumber coe = el -> second, term = 0.;
     57    int ind = el -> first -> Index ();
    5858
    59     while (nlt--) {
     59    if ((coe < 0.) && (sign < 0) ||
     60        (coe > 0.) && (sign > 0))
     61      {if    ((term=ub[ind]) > COUENNE_INFINITY) return (sign<0)? -COUENNE_INFINITY:COUENNE_INFINITY;}
     62    else {if ((term=lb[ind]) <-COUENNE_INFINITY) return (sign<0)? -COUENNE_INFINITY:COUENNE_INFINITY;}
    6063
    61 #ifdef DEBUG
    62       printf ("lin %d %g %g ", *li, *lc, (*lc < 0) ? ub [*li] : lb [*li]);
    63 #endif
     64    bound += coe * term;
     65  }
    6466
    65       if (*lc < 0) {if ((term = ub [*li++]) >  COUENNE_INFINITY) return -COUENNE_INFINITY;}
    66       else         {if ((term = lb [*li++]) < -COUENNE_INFINITY) return -COUENNE_INFINITY;}
     67  // derive quadratic part (obtain linear part)
     68  for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
    6769
    68       bound += *lc++ * term;
    69 
    70 #ifdef DEBUG
    71       printf (" --> %g\n", bound);
    72 #endif
    73     }
    74 
    75     while (nqt--) {
    76 
    77       int i = *qi++,
    78           j = *qj++;
     70    int xind = row -> first -> Index ();
    7971
    8072      CouNumber
    81         coe = *qc++,
    82         lbi = lb [i],
    83         ubi = ub [i];
     73        lbi = lb [xind],
     74        ubi = ub [xind];
    8475
    85       if (i==j) {
     76    for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
    8677
    87         if (coe > 0) term = (ubi <= 0) ? (ubi * ubi) : (lbi >= 0) ? (lbi * lbi) : 0;
    88         else if ((term = CoinMax (lbi*lbi, ubi*ubi)) > COUENNE_INFINITY)
    89           return -COUENNE_INFINITY;
     78      int yind = col -> first -> Index ();
     79
     80      CouNumber coe = col -> second;
     81
     82      if (xind == yind) { // term of the form q_ii x_i^2
     83
     84        if ((coe > 0.) && (sign < 0) ||
     85            (coe < 0.) && (sign > 0))
     86          term = (ubi < 0) ? (ubi * ubi) : (lbi > 0) ? (lbi * lbi) : 0.; //min{xi^2: xi in [lbi,ubi]
     87        else
     88          if ((term = CoinMax (lbi*lbi, ubi*ubi)) > COUENNE_INFINITY)
     89            return (sign < 0) ? -COUENNE_INFINITY : COUENNE_INFINITY;
    9090
    9191        term *= coe;
    9292
    9393#ifdef DEBUG
    94         printf ("Qii %d %g %g -> %g\n", i, coe, term, bound + term);
     94        printf ("Qii %d %g %g -> %g\n", xind, coe, term, bound + term);
    9595#endif
    9696      } else {
     
    9999
    100100        CouNumber
    101           lbj = lb [j], ubj = ub [j],
     101          lbj = lb [yind],
     102          ubj = ub [yind],
    102103          b1 = coe * lbi * lbj,
    103104          b2 = coe * lbi * ubj,
     
    110111        if (fabs (ubj) == 0) b2 = b4 = 0;
    111112
    112         if ((term = CoinMin (CoinMin (b1, b2), CoinMin (b3, b4))) < -COUENNE_INFINITY)
    113           return -COUENNE_INFINITY;
     113        if (sign < 0) {
     114          if ((term = CoinMin (CoinMin (b1, b2), CoinMin (b3, b4))) < -COUENNE_INFINITY)
     115            return -COUENNE_INFINITY;
     116        } else
     117          if ((term = CoinMax (CoinMax (b1, b2), CoinMax (b3, b4))) >  COUENNE_INFINITY)
     118            return  COUENNE_INFINITY;
    114119
    115120#ifdef DEBUG
    116         printf ("Qij %d %d %g %g -> %g\n", i, j, coe, term, bound + term);
     121        printf ("Qij %d %d %g %g -> %g\n", xind, yind, coe, term, bound + term);
    117122#endif
    118123      }
     
    121126      bound += term;
    122127    }
    123   } else { // compute upper bound /////////////////////////////////////////////////////////////
    124 
    125     while (nlt--) { // linear part
    126 
    127 #ifdef DEBUG
    128       printf ("lin %d %g %g %g\n", *li, *lc, (*lc < 0) ? ub [*li] : lb [*li], bound);
    129 #endif
    130 
    131       if (*lc > 0) {if ((term = ub [*li++]) >  COUENNE_INFINITY) return COUENNE_INFINITY;}
    132       else         {if ((term = lb [*li++]) < -COUENNE_INFINITY) return COUENNE_INFINITY;}
    133 
    134       bound += *lc++ * term;
    135     }
    136 
    137     while (nqt--) { // quadratic part
    138 
    139       int i = *qi++,
    140           j = *qj++;
    141 
    142       CouNumber
    143         coe = *qc++,
    144         lbi = lb [i],
    145         ubi = ub [i];
    146 
    147       if (i==j) {
    148 
    149         if (coe > 0) term = CoinMax (lbi * lbi, ubi * ubi);
    150         else         term = (ubi <= 0) ? (ubi * ubi) : (lbi >= 0) ? (lbi * lbi) : 0;
    151 
    152         if (term > COUENNE_INFINITY)
    153           return COUENNE_INFINITY;
    154 
    155         term *= coe;
    156 
    157 #ifdef DEBUG
    158         printf ("Qii %d %g %g --> %g ", i, coe, term, bound + term);
    159 #endif
    160       } else {
    161 
    162         coe *= 2;
    163 
    164         CouNumber
    165           lbj = lb [j],
    166           ubj = ub [j],
    167           b1 = coe * lbi * lbj,
    168           b2 = coe * lbi * ubj,
    169           b3 = coe * ubi * lbj,
    170           b4 = coe * ubi * ubj;
    171 
    172         // I hate this... but when you see
    173         //
    174         //   CoinMax (CoinMax (-0, nan), CoinMax (nan, -inf)) =
    175         // = CoinMax (nan, -inf) = -inf
    176         //
    177         // you feel changed.
    178 
    179         if (fabs (lbi) == 0) b1 = b2 = 0;
    180         if (fabs (lbj) == 0) b1 = b3 = 0;
    181         if (fabs (ubi) == 0) b3 = b4 = 0;
    182         if (fabs (ubj) == 0) b2 = b4 = 0;
    183 
    184         if ((term = CoinMax (CoinMax (b1, b2), CoinMax (b3, b4))) > COUENNE_INFINITY)
    185           return COUENNE_INFINITY;
    186 
    187 #ifdef DEBUG
    188         printf ("Qij %d %d %g %g -> %g [%g,%g,%g,%g] (%g,%g,%g,%g) {%g,%g,%g}\n",
    189                 i, j, coe, term, bound + term,
    190                 lbi, ubi, lbj, ubj,
    191                 b1, b2, b3, b4,
    192                 CoinMax (b1, b2), CoinMax (b3, b4),
    193                 CoinMax (CoinMax (b1, b2), CoinMax (b3, b4)));
    194 #endif
    195       }
    196 
    197       bound += term;
    198     }
    199128  }
    200129
  • branches/Couenne/Couenne/src/expression/operators/bounds/exprBQuad.hpp

    r960 r1017  
    1515
    1616
    17 /// method to actually compute the bound
    18 CouNumber computeQBound (int, exprQuad *);
    19 
    20 
    2117/// class to compute lower bound of a fraction based on the bounds of
    2218/// both numerator and denominator
     
    2420class exprLBQuad: public expression {
    2521
    26   exprQuad *ref_; //< quadratic form, reference expression
     22  exprQuad *ref_; ///< quadratic form, reference expression
    2723
    2824 public:
    2925
    3026  /// Constructor
    31   exprLBQuad (exprQuad *ref): ref_ (ref) {}
     27  exprLBQuad (exprQuad *ref):
     28    ref_ (ref) {}
    3229
    3330  /// copy constructor
     
    4441  /// function for the evaluation of the expression
    4542  inline CouNumber operator () ()
    46   {return computeQBound (-1, ref_);}
     43  {return ref_ -> computeQBound (-1);}
    4744
    4845  /// I/O
    4946  virtual void print (std::ostream &s = std::cout,     //< output stream
    50                       bool descend = false,            //< descend into auxiliaries' image?
    51                       CouenneProblem *p = NULL) const  //< problem pointer (in exprGroup)
     47                      bool descend = false) const      //< descend into auxiliaries' image?
    5248
    53   {s << "quadLower("; ref_ -> print (s, descend, p); s << ')';}
     49  {s << "quadLower("; ref_ -> print (s, descend); s << ')';}
    5450};
    5551
     
    6056class exprUBQuad: public expression {
    6157
    62   exprQuad *ref_; //< quadratic form, reference expression
     58  exprQuad *ref_; ///< quadratic form, reference expression
    6359
    6460 public:
    6561
    6662  /// Constructor
    67   exprUBQuad (exprQuad *ref): ref_ (ref) {}
     63  exprUBQuad (exprQuad *ref):
     64    ref_ (ref) {}
    6865
    6966  /// copy constructor
     
    8077  /// function for the evaluation of the expression
    8178  inline CouNumber operator () ()
    82   {return computeQBound (1, ref_);}
     79  {return ref_ -> computeQBound (1);}
    8380
    8481  /// I/O
    8582  virtual void print (std::ostream &s = std::cout,     //< output stream
    86                       bool descend = false,            //< descend into auxiliaries' image?
    87                       CouenneProblem *p = NULL) const  //< problem pointer (in exprGroup)
     83                      bool descend = false) const      //< descend into auxiliaries' image?
    8884
    89   {s << "quadUpper("; ref_ -> print (s, descend, p); s << ')';}
     85  {s << "quadUpper("; ref_ -> print (s, descend); s << ')';}
    9086};
    9187
  • branches/Couenne/Couenne/src/expression/operators/compQuadFinBounds.cpp

    r854 r1017  
    1212
    1313
    14 void updateInf (CouNumber, CouNumber, CouNumber, int &, int &, int );
     14void updateInf (CouNumber coe, CouNumber lb, CouNumber ub, int &indInf1, int &indInf2, int i) {
    1515
     16  if (coe > 0) {
     17    if (lb < 0) indInf1 = (indInf1 == -1) ? i : -2;
     18    if (ub > 0) indInf2 = (indInf2 == -1) ? i : -2;
     19  } else {
     20    if (lb < 0) indInf2 = (indInf2 == -1) ? i : -2;
     21    if (ub > 0) indInf1 = (indInf1 == -1) ? i : -2;
     22  }
     23}
    1624
    17 void computeQuadFiniteBound (const exprQuad *e,
    18                              CouNumber &qMin, CouNumber &qMax,
    19                              CouNumber *l, CouNumber *u,
    20                              int &indInfLo, int &indInfUp) {
     25// compute bound of a quadratic expression neglecting the infinite
     26// bounds of single variables
    2127
    22   exprQuad *eq = const_cast <exprQuad *> (e);
    23 
    24   int
    25     nlt = eq -> getnLTerms (),
    26     *li = eq -> getIndices (),
    27 
    28     nqt = eq -> getnQTerms (),
    29     *qi = eq -> getQIndexI (),
    30     *qj = eq -> getQIndexJ ();
    31 
    32   CouNumber
    33     *lc = eq -> getCoeffs  (),
    34     *qc = eq -> getQCoeffs ();
     28void exprQuad::computeQuadFiniteBound (CouNumber &qMin,
     29                                       CouNumber &qMax,
     30                                       CouNumber *l,
     31                                       CouNumber *u,
     32                                       int &indInfLo,
     33                                       int &indInfUp) {
    3534
    3635  // linear terms ///////////////////////////////////////////////////
    3736
    38   while (nlt--) {
     37  for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el) {
    3938
    40     int ind = *li++;
     39    int ind = el -> first -> Index ();
    4140
    4241    CouNumber
    43       coe = *lc++,
     42      coe = el -> second,
    4443      li  = l [ind],
    4544      ui  = u [ind];
     
    7069  }
    7170
    72   // quadratic terms ////////////////////////////////////////////////
     71  // quadratic terms ///////////////////////////////////////////////////
    7372
    74   while (nqt--) {
     73  for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
    7574
    76     int i = *qi++,
    77         j = *qj++;
     75    int i = row -> first -> Index ();
    7876
    79     CouNumber
    80       coe = *qc++,
    81       lbi = l [i],
    82       ubi = u [i];
     77    for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
    8378
    84     if (i==j) { // term of the form q_{ii} x_i^2
     79      int j = col -> first -> Index ();
    8580
    86       CouNumber
    87         tmin = (ubi <= 0) ? (ubi * ubi) : (lbi >= 0) ? (lbi * lbi) : 0,
    88         tmax = CoinMax (lbi*lbi, ubi*ubi);
     81      CouNumber
     82        coe = col -> second,
     83        lbi = l [i],
     84        ubi = u [i];
    8985
    90       if (coe < 0) { // negative coefficient, q_{ii} < 0
    91         qMax += coe * tmin;
    92         if (indInfLo > -2) {
    93           if (tmax > COUENNE_INFINITY) indInfLo = (indInfLo == -1) ? i : -2;
    94           else qMin += coe * tmax;
     86      if (i==j) { // term of the form q_{ii} x_i^2
     87
     88        CouNumber
     89          tmin = (ubi <= 0) ? (ubi * ubi) : (lbi >= 0) ? (lbi * lbi) : 0,
     90          tmax = CoinMax (lbi*lbi, ubi*ubi);
     91
     92        if (coe < 0) { // negative coefficient, q_{ii} < 0
     93          qMax += coe * tmin;
     94          if (indInfLo > -2) {
     95            if (tmax > COUENNE_INFINITY) indInfLo = (indInfLo == -1) ? i : -2;
     96            else qMin += coe * tmax;
     97          }
     98        } else { // positive coefficient
     99          qMin += coe * tmin;
     100          if (indInfUp > -2) {
     101            if (tmax > COUENNE_INFINITY) indInfUp = (indInfUp == -1) ? i : -2;
     102            else qMax += coe * tmax;
     103          }
    95104        }
    96       } else { // positive coefficient
    97         qMin += coe * tmin;
    98         if (indInfUp > -2) {
    99           if (tmax > COUENNE_INFINITY) indInfUp = (indInfUp == -1) ? i : -2;
    100           else qMax += coe * tmax;
    101         }
     105      } else { // term of the form q_{ij} x_i x_j, j\neq i
     106
     107        CouNumber
     108          lbj = l [j],
     109          ubj = u [j];
     110
     111        coe *= 2;
     112
     113        // check if index i wrings unbounded value in both directions
     114
     115        if (lbi < -COUENNE_INFINITY) updateInf (coe, lbj, ubj, indInfUp, indInfLo, i);
     116        if (lbj < -COUENNE_INFINITY) updateInf (coe, lbi, ubi, indInfUp, indInfLo, j);
     117
     118        if (ubi >  COUENNE_INFINITY) updateInf (coe, lbj, ubj, indInfLo, indInfUp, i);
     119        if (ubj >  COUENNE_INFINITY) updateInf (coe, lbi, ubi, indInfLo, indInfUp, j);
     120
     121        CouNumber term,
     122          b1 = coe * lbi * lbj,
     123          b2 = coe * lbi * ubj,
     124          b3 = coe * ubi * lbj,
     125          b4 = coe * ubi * ubj;
     126
     127        if ((term = CoinMin (CoinMin (b1, b2), CoinMin (b3, b4))) > -COUENNE_INFINITY) qMin += term;
     128        if ((term = CoinMax (CoinMax (b1, b2), CoinMax (b3, b4))) <  COUENNE_INFINITY) qMax += term;
    102129      }
    103     } else { // term of the form q_{ij} x_i x_j, j\neq i
    104 
    105       CouNumber
    106         lbj = l [j],
    107         ubj = u [j];
    108 
    109       coe *= 2;
    110 
    111       // check if index i wrings unbounded value in both directions
    112 
    113       if (lbi < -COUENNE_INFINITY) updateInf (coe, lbj, ubj, indInfUp, indInfLo, i);
    114       if (lbj < -COUENNE_INFINITY) updateInf (coe, lbi, ubi, indInfUp, indInfLo, j);
    115 
    116       if (ubi >  COUENNE_INFINITY) updateInf (coe, lbj, ubj, indInfLo, indInfUp, i);
    117       if (ubj >  COUENNE_INFINITY) updateInf (coe, lbi, ubi, indInfLo, indInfUp, j);
    118 
    119       CouNumber term,
    120         b1 = coe * lbi * lbj,
    121         b2 = coe * lbi * ubj,
    122         b3 = coe * ubi * lbj,
    123         b4 = coe * ubi * ubj;
    124 
    125       if ((term = CoinMin (CoinMin (b1, b2), CoinMin (b3, b4))) > -COUENNE_INFINITY) qMin += term;
    126       if ((term = CoinMax (CoinMax (b1, b2), CoinMax (b3, b4))) <  COUENNE_INFINITY) qMax += term;
    127130    }
    128131  }
    129132}
    130 
    131 
    132 void updateInf (CouNumber coe, CouNumber lb, CouNumber ub, int &indInf1, int &indInf2, int i) {
    133 
    134   if (coe > 0) {
    135     if (lb < 0) indInf1 = (indInf1 == -1) ? i : -2;
    136     if (ub > 0) indInf2 = (indInf2 == -1) ? i : -2;
    137   } else {
    138     if (lb < 0) indInf2 = (indInf2 == -1) ? i : -2;
    139     if (ub > 0) indInf1 = (indInf1 == -1) ? i : -2;
    140   }
    141 }
  • branches/Couenne/Couenne/src/expression/operators/exprAbs.hpp

    r960 r1017  
    4545
    4646  /// generate equality between *this and *w
    47   void generateCuts (exprAux *w, const OsiSolverInterface &si,
     47  void generateCuts (expression *w, const OsiSolverInterface &si,
    4848                     OsiCuts &cs, const CouenneCutGenerator *cg,
    4949                     t_chg_bounds * = NULL, int = -1,
     
    5555
    5656  /// is this expression integer?
    57   bool isInteger ()
     57  inline bool isInteger ()
    5858  {return argument_ -> isInteger ();}
    5959
     
    6363  /// set up branching object by evaluating many branching points for
    6464  /// each expression's arguments
    65   CouNumber selectBranch (const CouenneObject *obj, const OsiBranchingInformation *,
    66                           int &, double * &, int &);
     65  virtual CouNumber selectBranch (const CouenneObject *obj,
     66                                  const OsiBranchingInformation *info,
     67                                  expression * &var,
     68                                  double * &brpts,
     69                                  int &way);
    6770};
    6871
  • branches/Couenne/Couenne/src/expression/operators/exprCos.hpp

    r979 r1017  
    4242
    4343  /// generate equality between *this and *w
    44   void generateCuts (exprAux *w, const OsiSolverInterface &si,
     44  void generateCuts (expression *w, const OsiSolverInterface &si,
    4545                     OsiCuts &cs, const CouenneCutGenerator *cg,
    4646                     t_chg_bounds * = NULL, int = -1,
     
    5858  /// Set up branching object by evaluating many branching points for
    5959  /// each expression's arguments
    60   CouNumber selectBranch (const CouenneObject *obj, const OsiBranchingInformation *info,
    61                           int &ind, double * &brpts, int &way)
    62   {return trigSelBranch (obj, info, ind, brpts, way, COU_COSINE);}
     60  virtual CouNumber selectBranch (const CouenneObject *obj,
     61                                  const OsiBranchingInformation *info,
     62                                  expression * &var,
     63                                  double * &brpts,
     64                                  int &way)
     65  {return trigSelBranch (obj, info, var, brpts, way, COU_COSINE);}
    6366};
    6467
  • branches/Couenne/Couenne/src/expression/operators/exprDiv.cpp

    r960 r1017  
    7878expression *exprDiv::differentiate (int index) {
    7979
    80   if (!(arglist_ [0] -> dependsOn (&index, 1))  &&
    81       !(arglist_ [1] -> dependsOn (&index, 1)))
     80  if (!(arglist_ [0] -> dependsOn (index))  &&
     81      !(arglist_ [1] -> dependsOn (index)))
    8282    return new exprConst (0.);
    8383
  • branches/Couenne/Couenne/src/expression/operators/exprDiv.hpp

    r988 r1017  
    6565
    6666  /// Generate equality between *this and *w
    67   void generateCuts (exprAux *w, const OsiSolverInterface &si,
     67  void generateCuts (expression *w, const OsiSolverInterface &si,
    6868                     OsiCuts &cs, const CouenneCutGenerator *cg,
    6969                     t_chg_bounds * = NULL, int = -1,
     
    8686  /// Set up branching object by evaluating many branching points for
    8787  /// each expression's arguments
    88   CouNumber selectBranch (const CouenneObject *, const OsiBranchingInformation *,
    89                           int &, double * &, int &);
     88  virtual CouNumber selectBranch (const CouenneObject *obj,
     89                                  const OsiBranchingInformation *info,
     90                                  expression * &var,
     91                                  double * &brpts,
     92                                  int &way);
    9093};
    9194
  • branches/Couenne/Couenne/src/expression/operators/exprExp.hpp

    r960 r101