compute_symmetry_nauty.c

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the program and library             */
/*         SCIP --- Solving Constraint Integer Programs                      */
/*                                                                           */
/*  Copyright (c) 2002-2025 Zuse Institute Berlin (ZIB)                      */
/*                                                                           */
/*  Licensed under the Apache License, Version 2.0 (the "License");          */
/*  you may not use this file except in compliance with the License.         */
/*  You may obtain a copy of the License at                                  */
/*                                                                           */
/*      http://www.apache.org/licenses/LICENSE-2.0                           */
/*                                                                           */
/*  Unless required by applicable law or agreed to in writing, software      */
/*  distributed under the License is distributed on an "AS IS" BASIS,        */
/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
/*  See the License for the specific language governing permissions and      */
/*  limitations under the License.                                           */
/*                                                                           */
/*  You should have received a copy of the Apache-2.0 license                */
/*  along with SCIP; see the file LICENSE. If not visit scipopt.org.         */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
/**@file   compute_symmetry_nauty.c
 * @brief  interface for symmetry computations to nauty/traces
 * @author Marc Pfetsch
 */
 
/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
 
#include "compute_symmetry.h"
 
/* the following determines whether nauty or traces is used: */
#define NAUTY
 
/* include nauty/traces */
/* turn off warning (just for including nauty/traces) */
#pragma GCC diagnostic ignored "-Wunused-variable"
#pragma GCC diagnostic ignored "-Wredundant-decls"
#pragma GCC diagnostic ignored "-Wpedantic"
 
#ifdef NAUTY
#include "nauty/nauty.h"
#include "nauty/nausparse.h"
#else
#include "nauty/traces.h"
#endif
 
#pragma GCC diagnostic warning "-Wunused-variable"
#pragma GCC diagnostic warning "-Wredundant-decls"
#pragma GCC diagnostic warning "-Wpedantic"
 
#include "scip/symmetry_graph.h"
#include "scip/expr_var.h"
#include "scip/expr_sum.h"
#include "scip/expr_pow.h"
#include "scip/expr.h"
#include "scip/cons_nonlinear.h"
#include "scip/cons_linear.h"
#include "scip/scip_mem.h"
#include "tinycthread/tinycthread.h"
 
/** struct for nauty callback */
struct NAUTY_Data
{
   SCIP*                 scip;               /**< SCIP pointer */
   SYM_SYMTYPE           symtype;            /**< type of symmetries that need to be computed */
   int                   npermvars;          /**< number of variables for permutations */
   int                   nperms;             /**< number of permutations */
   int**                 perms;              /**< permutation generators as (nperms x npermvars) matrix */
   int                   nmaxperms;          /**< maximal number of permutations */
   int                   maxgenerators;      /**< maximal number of generators to be constructed (= 0 if unlimited) */
   SCIP_Bool             restricttovars;     /**< whether permutations shall be restricted to variables */
   int                   maxlevel;           /**< maximum depth level of nauty's search tree (-1: unlimited) */
};
typedef struct NAUTY_Data NAUTY_DATA;
 
/** static data for nauty callback */
#if defined(_Thread_local)
static _Thread_local NAUTY_DATA data_;
#else
static NAUTY_DATA data_;
#endif
 
/* ------------------- hook functions ------------------- */
 
#ifdef NAUTY
 
/** callback function for nauty */  /*lint -e{715}*/
static
void nautyhook(
   int                   count,              /**< ID of this generator */
   int*                  p,                  /**< generator (permutation) that nauty found */
   int*                  orbits,             /**< orbits generated by the group found so far */
   int                   numorbits,          /**< number of orbits */
   int                   stabvertex,         /**< stabilizing node */
   int                   n                   /**< number of nodes in the graph */
   )
{  /* lint --e{715} */
   SCIP_Bool isidentity = TRUE;
   int permlen;
   int* pp;
 
   assert( p != NULL );
 
   /* make sure we do not generate more than maxgenerators many permutations */
   if ( data_.maxgenerators != 0 && data_.nperms >= data_.maxgenerators )
   {
      /* request a kill from nauty */
      nauty_kill_request = 1;
      return;
   }
 
   /* copy first part of automorphism */
   if ( data_.restricttovars )
   {
      if ( data_.symtype == SYM_SYMTYPE_PERM )
         permlen = data_.npermvars;
      else
         permlen = 2 * data_.npermvars;
   }
   else
      permlen = n;
 
   /* check whether permutation is identity */
   for (int j = 0; j < permlen; ++j)
   {
      if ( (int) p[j] != j )
         isidentity = FALSE;
   }
 
   /* don't store identity permutations */
   if ( isidentity )
      return;
 
   /* check whether we should allocate space for perms */
   if ( data_.nmaxperms <= 0 )
   {
      if ( data_.maxgenerators == 0 )
         data_.nmaxperms = 100;   /* seems to cover many cases */
      else
         data_.nmaxperms = data_.maxgenerators;
 
      if ( SCIPallocBlockMemoryArray(data_.scip, &data_.perms, data_.nmaxperms) != SCIP_OKAY )
         return;
   }
   else if ( data_.nperms >= data_.nmaxperms )    /* check whether we need to resize */
   {
      int newsize;
 
      newsize = SCIPcalcMemGrowSize(data_.scip, data_.nperms + 1);
      assert( newsize >= data_.nperms );
      assert( data_.maxgenerators == 0 );
 
      if ( SCIPreallocBlockMemoryArray(data_.scip, &data_.perms, data_.nmaxperms, newsize) != SCIP_OKAY )
         return;
 
      data_.nmaxperms = newsize;
   }
 
   if ( SCIPduplicateBlockMemoryArray(data_.scip, &pp, p, permlen) != SCIP_OKAY )
      return;
   data_.perms[data_.nperms++] = pp;
}
 
/** callback function for nauty to terminate early */  /*lint -e{715}*/
static
void nautyterminationhook(
   graph*                g,                  /**< sparse graph for nauty */
   int*                  lab,                /**< labels of node */
   int*                  ptn,                /**< array indicating change of set in node parition of graph */
   int                   level,              /**< level of current node in nauty's tree */
   int                   numcells,           /**< number of cells in current partition */
   int                   tc,                 /**< index of target cells in lab if children need to be explored */
   int                   code,               /**< code produced by refinement and vertex-invariant procedures */
   int                   m,                  /**< number of edges in the graph */
   int                   n                   /**< number of nodes in the graph */
   )
{  /* lint --e{715} */
   /* add level limit to work around call stack overflow in nauty */
   if ( level > data_.maxlevel && data_.maxlevel >= 0 )
   {
      SCIPverbMessage(data_.scip, SCIP_VERBLEVEL_MINIMAL, NULL,
         "symmetry computation terminated early because Nauty level limit %d is exceeded\n",
         data_.maxlevel);
      SCIPverbMessage(data_.scip, SCIP_VERBLEVEL_MINIMAL, NULL,
         "for running full symmetry detection, increase value of parameter propagating/symmetry/nautymaxlevel\n");
      nauty_kill_request = 1;
   }
}
 
#else
 
/** callback function for traces */
static
void traceshook(
   int                   count,              /**< number of generator */
   int*                  p,                  /**< generator that traces found */
   int                   n                   /**< number of nodes in the graph */
   )
{
   SCIP_Bool isidentity = TRUE;
   int permlen;
   int* pp;
   int j;
 
   assert( p != NULL );
 
   /* make sure we do not generate more than maxgenerators many permutations */
   if ( data_.maxgenerators != 0 && data_.nperms >= data_.maxgenerators )
   {
      /* request a kill from traces */
      nauty_kill_request = 1;
      return;
   }
 
   /* copy first part of automorphism */
   if ( data_.restricttovars )
   {
      if ( data_.symtype == SYM_SYMTYPE_PERM )
         permlen = data_.npermvars;
      else
         permlen = 2 * data_.npermvars;
   }
   else
      permlen = n;
 
   /* check whether permutation is identity */
   for (j = 0; j < permlen; ++j)
   {
      if ( (int) p[j] != j )
         isidentity = FALSE;
   }
 
   /* ignore trivial generators, i.e. generators that only permute the constraints */
   if ( isidentity )
      return;
 
   /* check whether we should allocate space for perms */
   if ( data_.nmaxperms <= 0 )
   {
      if ( data_.maxgenerators == 0 )
         data_.nmaxperms = 100;   /* seems to cover many cases */
      else
         data_.nmaxperms = data_.maxgenerators;
 
      if ( SCIPallocBlockMemoryArray(data_.scip, &data_.perms, data_.nmaxperms) != SCIP_OKAY )
         return;
   }
   else if ( data_.nperms >= data_.nmaxperms )    /* check whether we need to resize */
   {
      int newsize;
 
      newsize = SCIPcalcMemGrowSize(data_.scip, data_.nperms + 1);
      assert( newsize >= data_.nperms );
      assert( data_.maxgenerators == 0 );
 
      if ( SCIPreallocBlockMemoryArray(data_.scip, &data_.perms, data_.nmaxperms, newsize) != SCIP_OKAY )
         return;
 
      data_.nmaxperms = newsize;
   }
 
   if ( SCIPduplicateBlockMemoryArray(data_.scip, &pp, p, permlen) != SCIP_OKAY )
      return;
   data_.perms[data_.nperms++] = pp;
}
 
#endif
 
/** returns whether an edge is considered in grouping process */
static
SCIP_Bool isEdgeGroupable(
   SYM_GRAPH*            symgraph,           /**< symmetry detection graph */
   int                   edgeidx,            /**< index of edge to be checked */
   SCIP_Bool             groupbycons         /**< whether edges are grouped by constraints */
   )
{
   int first;
   int second;
 
   assert(symgraph != NULL);
 
   first = SCIPgetSymgraphEdgeFirst(symgraph, edgeidx);
   second = SCIPgetSymgraphEdgeSecond(symgraph, edgeidx);
 
   /* uncolored edges are not grouped */
   if ( ! SCIPisSymgraphEdgeColored(symgraph, edgeidx) )
      return FALSE;
 
   /* two variable nodes are connected */
   if ( first < 0 && second < 0 )
      return FALSE;
 
   if ( ! groupbycons )
   {
      /* grouping by variables requires one variable node */
      if ( first < 0 || second < 0 )
         return TRUE;
   }
   else
   {
      /* check whether there is exactly one constraint node in edge */
      if ( first >= 0 && second >= 0 )
      {
         if ( (SCIPgetSymgraphNodeType(symgraph, first) == SYM_NODETYPE_CONS
               && SCIPgetSymgraphNodeType(symgraph, second) != SYM_NODETYPE_CONS)
            || (SCIPgetSymgraphNodeType(symgraph, first) != SYM_NODETYPE_CONS
               && SCIPgetSymgraphNodeType(symgraph, second) == SYM_NODETYPE_CONS) )
            return TRUE;
      }
      else if ( first >= 0 )
      {
         if ( SCIPgetSymgraphNodeType(symgraph, first) == SYM_NODETYPE_CONS )
            return TRUE;
      }
      else
      {
         if ( SCIPgetSymgraphNodeType(symgraph, second) == SYM_NODETYPE_CONS )
            return TRUE;
      }
   }
 
   return FALSE;
}
 
/** adds grouped edges all of which have one common endpoint to a graph
 *
 * The grouping mechanism works by sorting the edges according to their color. If two
 * edges have the same color, they share the same intermediate node, which is connected
 * to the common node and the other endpoints of equivalent edges.
 */
static
SCIP_RETCODE addOrDetermineEffectOfGroupedEdges(
   SCIP*                 scip,               /**< SCIP pointer */
   sparsegraph*          SG,                 /**< graph to be constructed */
   int*                  edgestartpos,       /**< initialized array of starting positions of new edges for each node */
   SCIP_Bool             determinesize,      /**< whether only the effect of grouping on the graph shall be checked */
   int*                  internodeid,        /**< (initialized) pointer to store the ID of the next intermediate node */
   int**                 degrees,            /**< pointer to array of degrees for nodes in G */
   int*                  maxdegrees,         /**< (initialized) pointer to maximum number of entries degrees can hold */
   int**                 colorstostore,      /**< pointer to array of colors of graph to be constructed */
   int*                  ncolorstostore,     /**< (initialized) pointer to maximum number of entries in colorstostore */
   int*                  nnodes,             /**< (initialized) pointer to store the number of */
   int*                  nedges,             /**< (initialized) pointer to store the number of */
   int                   commonnodeidx,      /**< index of common node in G */
   int*                  neighbors,          /**< neighbors of common node */
   int*                  colors,             /**< colors of edges to neighbors */
   int                   nneighbors,         /**< number of neighbors */
   int*                  naddednodes,        /**< pointer to store number of nodes added to G */
   int*                  naddededges         /**< pointer to store number of edges added to G */
   )
{
   int curcolor;
   int curstart;
   int e;
   int f;
 
   assert( SG != NULL || determinesize );
   assert( internodeid != NULL );
   assert( *internodeid >= 0 );
   assert( degrees != NULL );
   assert( maxdegrees != NULL );
   assert( *maxdegrees > 0 );
   assert( nnodes != NULL );
   assert( nedges != NULL );
   assert( neighbors != NULL );
   assert( colors != NULL );
   assert( naddednodes != NULL );
   assert( naddededges != NULL );
   assert( commonnodeidx <= *internodeid );
 
   *naddednodes = 0;
   *naddededges = 0;
 
   /* sort edges according to color */
   SCIPsortIntInt(colors, neighbors, nneighbors);
 
   /* iterate over groups of identical edges and group them, ignoring the last group */
   curcolor = colors[0];
   curstart = 0;
   for (e = 1; e < nneighbors; ++e)
   {
      /* if a new group started, add edges for previous group */
      if ( colors[e] != curcolor )
      {
         if ( determinesize )
         {
            SCIP_CALL( SCIPensureBlockMemoryArray(scip, degrees, maxdegrees, *internodeid + 1) );
            SCIP_CALL( SCIPensureBlockMemoryArray(scip, colorstostore, ncolorstostore, *internodeid + 1) );
            (*degrees)[*internodeid] = 1;
            ++(*degrees)[commonnodeidx];
            (*colorstostore)[*internodeid] = curcolor;
 
            for (f = curstart; f < e; ++f)
            {
               assert( neighbors[f] <= *internodeid );
               ++(*degrees)[*internodeid];
               ++(*degrees)[neighbors[f]];
            }
         }
         else
         {
            SG->e[edgestartpos[commonnodeidx]++] = *internodeid;
            SG->e[edgestartpos[*internodeid]++] = commonnodeidx;
 
            for (f = curstart; f < e; ++f)
            {
               SG->e[edgestartpos[neighbors[f]]++] = *internodeid;
               SG->e[edgestartpos[*internodeid]++] = neighbors[f];
            }
         }
         *naddednodes += 1;
         *naddededges += e - curstart + 1;
         ++(*internodeid);
 
         curcolor = colors[e];
         curstart = e;
      }
   }
 
   /* add edges of last group */
   if ( determinesize )
   {
      SCIP_CALL( SCIPensureBlockMemoryArray(scip, degrees, maxdegrees, *internodeid + 1) );
      SCIP_CALL( SCIPensureBlockMemoryArray(scip, colorstostore, ncolorstostore, *internodeid + 1) );
 
      (*degrees)[*internodeid] = 1;
      ++(*degrees)[commonnodeidx];
      (*colorstostore)[*internodeid] = curcolor;
 
      for (f = curstart; f < nneighbors; ++f)
      {
         assert( neighbors[f] <= *internodeid );
         ++(*degrees)[*internodeid];
         ++(*degrees)[neighbors[f]];
      }
   }
   else
   {
      SG->e[edgestartpos[commonnodeidx]++] = *internodeid;
      SG->e[edgestartpos[*internodeid]++] = commonnodeidx;
 
      for (f = curstart; f < nneighbors; ++f)
      {
         SG->e[edgestartpos[*internodeid]++] = neighbors[f];
         SG->e[edgestartpos[neighbors[f]]++] = *internodeid;
      }
   }
   *naddednodes += 1;
   *naddededges += nneighbors - curstart + 1;
   ++(*internodeid);
 
   return SCIP_OKAY;
}
 
/** either creates a graph or determines its size */
static
SCIP_RETCODE createOrDetermineSizeGraph(
   SCIP*                 scip,               /**< SCIP instance */
   SYM_GRAPH*            symgraph,           /**< symmetry detection graph */
   SCIP_Bool             determinesize,      /**< whether only the size of the graph shall be determined */
   sparsegraph*          SG,                 /**< graph to be constructed */
   int*                  nnodes,             /**< pointer to store the total number of nodes in graph */
   int*                  nedges,             /**< pointer to store the total number of edges in graph */
   int**                 degrees,            /**< pointer to store the degrees of the nodes */
   int*                  maxdegrees,         /**< pointer to store the maximal size of the degree array */
   int**                 colors,             /**< pointer to store the colors of the nodes */
   int*                  ncolors,            /**< pointer to store number of different colors in graph */
   SCIP_Bool*            success             /**< pointer to store whether the construction was successful */
   )
{ /*lint !e438*/
   SYM_SYMTYPE symtype;
   SYM_NODETYPE comparetype;
   SCIP_Bool groupByConstraints;
   int* groupfirsts = NULL;
   int* groupseconds = NULL;
   int* groupcolors = NULL;
   int* pos = NULL;
   int edgebegincnt = 0;
   int ngroupedges = 0;
   int nvarnodestoadd;
   int internodeid;
   int nsymvars;
   int nsymedges;
   int first;
   int second;
   int e;
   int j;
 
   assert( scip != NULL );
   assert( symgraph != NULL );
   assert( SG != NULL || determinesize );
   assert( nnodes != NULL );
   assert( nedges != NULL );
   assert( degrees != NULL );
   assert( maxdegrees != NULL );
   assert( colors != NULL );
   assert( ncolors != NULL );
   assert( success != NULL );
 
   *success = TRUE;
 
   /* collect basic information from symmetry detection graph */
   nsymvars = SCIPgetSymgraphNVars(symgraph);
   symtype = SCIPgetSymgraphSymtype(symgraph);
   switch ( symtype )
   {
   case SYM_SYMTYPE_PERM:
      nvarnodestoadd = nsymvars;
      break;
   default:
      assert( symtype == SYM_SYMTYPE_SIGNPERM );
      nvarnodestoadd = 2 * nsymvars;
   } /*lint !e788*/
 
   /* possibly find number of nodes in sassy graph */
   if ( determinesize )
   {
      int cnt = 0;
 
      /* initialize counters */
      *nnodes = SCIPgetSymgraphNNodes(symgraph) + nvarnodestoadd;
      *nedges = 0;
 
      /* possibly allocate memory for degrees (will grow dynamically) */
      *degrees = NULL;
      *maxdegrees = 0;
      SCIP_CALL( SCIPensureBlockMemoryArray(scip, degrees, maxdegrees, *nnodes + 100) );
      for (j = 0; j < *nnodes; ++j)
         (*degrees)[j] = 0;
 
      /* possibly allocate memory for colors (will grow dynamically) */
      *colors = NULL;
      *ncolors = 0;
      SCIP_CALL( SCIPensureBlockMemoryArray(scip, colors, ncolors, *nnodes + 100) );
      for (j = 0; j < nvarnodestoadd; ++j)
         (*colors)[cnt++] = SCIPgetSymgraphVarnodeColor(symgraph, j);
      for (j = 0; j < SCIPgetSymgraphNNodes(symgraph); ++j)
         (*colors)[cnt++] = SCIPgetSymgraphNodeColor(symgraph, j);
   }
   else
   {
      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &pos, *nnodes) );
 
      /* add nodes for variables and remaining (axuiliary) nodes in graph */
      for (j = 0; j < *nnodes; ++j)
      {
         SG->d[j] = (*degrees)[j];
         SG->v[j] = (size_t) (unsigned) edgebegincnt;
         pos[j] = edgebegincnt;
         edgebegincnt += (*degrees)[j];
      }
   }
 
   /* determine grouping depending on the number of rhs vs. variables */
   groupByConstraints = SCIPgetSymgraphNConsnodes(symgraph) < SCIPgetSymgraphNVars(symgraph);
 
   /* allocate arrays to collect edges to be grouped */
   nsymedges = SCIPgetSymgraphNEdges(symgraph);
   SCIP_CALL( SCIPallocBufferArray(scip, &groupfirsts, nsymedges) );
   SCIP_CALL( SCIPallocBufferArray(scip, &groupseconds, nsymedges) );
   SCIP_CALL( SCIPallocBufferArray(scip, &groupcolors, nsymedges) );
 
   /* loop through all edges of the symmetry detection graph and either get degrees of nodes or add edges */
   internodeid = SCIPgetSymgraphNNodes(symgraph) + nvarnodestoadd;
   for (e = 0; e < SCIPgetSymgraphNEdges(symgraph); ++e)
   {
      first = SCIPgetSymgraphEdgeFirst(symgraph, e);
      second = SCIPgetSymgraphEdgeSecond(symgraph, e);
 
      /* get the first and second node in edge (corrected by variable shift) */
      if ( first < 0 )
         first = -first - 1;
      else
         first += nvarnodestoadd;
      if ( second < 0 )
         second = -second - 1;
      else
         second += nvarnodestoadd;
      assert(first >= 0);
      assert(second >= 0);
 
      /* check whether edge is used for grouping */
      if ( ! SCIPhasGraphUniqueEdgetype(symgraph) && isEdgeGroupable(symgraph, e, groupByConstraints) )
      {
         /* store edge, first becomes the cons or var node */
         comparetype = groupByConstraints ? SYM_NODETYPE_CONS : SYM_NODETYPE_VAR;
 
         if ( SCIPgetSymgraphNodeType(symgraph, SCIPgetSymgraphEdgeFirst(symgraph, e)) == comparetype )
         {
            groupfirsts[ngroupedges] = first;
            groupseconds[ngroupedges] = second;
         }
         else
         {
            groupfirsts[ngroupedges] = second;
            groupseconds[ngroupedges] = first;
         }
         groupcolors[ngroupedges++] = SCIPgetSymgraphEdgeColor(symgraph, e);
      }
      else
      {
         /* immediately add edge or increase degrees */
         assert(0 <= first && first < *nnodes);
         assert(0 <= second && second < *nnodes);
 
         /* possibly split edge if it is colored */
         if ( ! SCIPhasGraphUniqueEdgetype(symgraph) && SCIPisSymgraphEdgeColored(symgraph, e) )
         {
            if ( determinesize )
            {
               SCIP_CALL( SCIPensureBlockMemoryArray(scip, degrees, maxdegrees, internodeid + 1) );
               SCIP_CALL( SCIPensureBlockMemoryArray(scip, colors, ncolors, internodeid + 1) );
               ++(*degrees)[first];
               ++(*degrees)[second];
               (*degrees)[internodeid] = 2;
               ++(*nnodes);
               *nedges += 2;
               (*colors)[internodeid] = SCIPgetSymgraphEdgeColor(symgraph, e);
               ++internodeid;
            }
            else
            {
               assert( internodeid < *nnodes );
 
               /* add (bidirected) edges */
               SG->e[pos[first]++] = internodeid;
               SG->e[pos[internodeid]++] = first;
               SG->e[pos[second]++] = internodeid;
               SG->e[pos[internodeid]++] = second;
 
               assert( first == *nnodes - 1 || pos[first] <= (int) SG->v[first+1] );
               assert( second == *nnodes - 1 || pos[second] <= (int) SG->v[second+1] );
               assert( internodeid == *nnodes - 1 || pos[internodeid] <= (int) SG->v[internodeid+1] );
 
               ++internodeid;
            }
         }
         else
         {
            if ( determinesize )
            {
               ++(*degrees)[first];
               ++(*degrees)[second];
               ++(*nedges);
            }
            else
            {
               /* add (bidirected) edge */
               SG->e[pos[first]++] = second;
               SG->e[pos[second]++] = first;
 
               assert( first == *nnodes - 1 || pos[first] <= (int) SG->v[first+1] );
               assert( second == *nnodes - 1 || pos[second] <= (int) SG->v[second+1] );
            }
         }
      }
   }
 
   /* possibly add groupable edges */
   if ( ngroupedges > 0 )
   {
      int firstidx = 0;
      int firstnodeidx;
      int naddednodes;
      int naddededges;
 
      /* sort edges according to their first nodes */
      SCIPsortIntIntInt(groupfirsts, groupseconds, groupcolors, ngroupedges);
      firstnodeidx = groupfirsts[0];
 
      for (j = 1; j < ngroupedges; ++j)
      {
         /* if a new first node has been found, group the edges of the previous first node; ignoring the last group */
         if ( groupfirsts[j] != firstnodeidx )
         {
            SCIP_CALL( addOrDetermineEffectOfGroupedEdges(scip, SG, pos, determinesize, &internodeid,
                  degrees, maxdegrees, colors, ncolors, nnodes, nedges, firstnodeidx, &groupseconds[firstidx],
                  &groupcolors[firstidx], j - firstidx, &naddednodes, &naddededges) );
 
            firstidx = j;
            firstnodeidx = groupfirsts[j];
 
            if ( determinesize )
            {
               *nnodes += naddednodes;
               *nedges += naddededges;
            }
         }
      }
 
      /* process the last group */
      SCIP_CALL( addOrDetermineEffectOfGroupedEdges(scip, SG, pos, determinesize, &internodeid,
            degrees, maxdegrees, colors, ncolors, nnodes, nedges, firstnodeidx, &groupseconds[firstidx],
            &groupcolors[firstidx], ngroupedges - firstidx, &naddednodes, &naddededges) );
 
      if ( determinesize )
      {
         *nnodes += naddednodes;
         *nedges += naddededges;
      }
   }
 
   SCIPfreeBufferArray(scip, &groupcolors);
   SCIPfreeBufferArray(scip, &groupseconds);
   SCIPfreeBufferArray(scip, &groupfirsts);
 
   /* for signed permutation, also add edges connecting a variable and its negation */
   switch ( SCIPgetSymgraphSymtype(symgraph) )
   {
   case SYM_SYMTYPE_SIGNPERM:
      if ( determinesize )
      {
         for (j = 0; j < nvarnodestoadd; ++j)
            ++(*degrees)[j];
         *nedges += nsymvars;
      }
      else
      {
         for (j = 0; j < nsymvars; ++j)
         {
            SG->e[pos[j]++] = j + nsymvars;
            SG->e[pos[j + nsymvars]++] = j;
 
            assert( pos[j] <= (int) SG->v[j+1] );
            assert( j + nsymvars == *nnodes - 1 || pos[j+nsymvars] <= (int) SG->v[j+nsymvars+1] );
         }
      }
      break;
   default:
      assert( SCIPgetSymgraphSymtype(symgraph) == SYM_SYMTYPE_PERM );
   } /*lint !e788*/
 
   if ( determinesize )
   {
      SCIPdebugMsg(scip, "#nodes: %d\n", *nnodes);
      SCIPdebugMsg(scip, "#nodes for variables: %d\n", nvarnodestoadd);
      SCIPdebugMsg(scip, "#nodes for rhs: %d\n", SCIPgetSymgraphNConsnodes(symgraph));
      SCIPdebugMsg(scip, "#edges: %d\n", *nedges);
   }
   else
   {
      SCIPfreeBlockMemoryArray(scip, &pos, *nnodes);
   }
 
   return SCIP_OKAY;
}
 
/** either creates a graph for checking symmetries or determines its size
 *
 *  The input are two graphs and the graph to be constructed consists of copies
 *  of the two input graphs, in which non-variable nodes are colored according
 *  to the colors used in symmetry detection. Each variable gets a unique color.
 *  Finally, a dummy node is introduced that is adjacent with all remaining nodes.
 */
static
SCIP_RETCODE createOrDetermineSizeGraphCheck(
   SCIP*                 scip,               /**< SCIP instance */
   SYM_GRAPH*            graph1,             /**< first symmetry detection graph */
   SYM_GRAPH*            graph2,             /**< second symmetry detection graph */
   SCIP_Bool             determinesize,      /**< whether only the size of the graph shall be determined */
   sparsegraph*          SG,                 /**< graph to be constructed */
   int*                  nnodes,             /**< pointer to store the total number of nodes in graph */
   int*                  nedges,             /**< pointer to store the total number of edges in graph */
   int**                 degrees,            /**< pointer to store the degrees of the nodes */
   int*                  maxdegrees,         /**< pointer to store the maximal size of the degree array */
   int**                 colors,             /**< pointer to store the colors of the nodes */
   int*                  ncolors,            /**< pointer to store number of different colors in graph */
   int*                  nusedvars,          /**< pointer to store number of variables used in one graph */
   int*                  nnodesfromgraph1,   /**< pointer to store number of nodes arising from graph1 (or NULL) */
   SCIP_Bool*            success             /**< pointer to store whether the graph could be built */
   )
{
   SYM_SYMTYPE symtype;
   SYM_NODETYPE comparetype;
   SCIP_Bool groupByConstraints;
   SYM_GRAPH* symgraph;
   int* nvarused1 = NULL;
   int* nvarused2 = NULL;
   int* varlabel = NULL;
   int* groupfirsts = NULL;
   int* groupseconds = NULL;
   int* groupcolors = NULL;
   int* pos = NULL;
   int nusdvars = 0;
   int edgebegincnt = 0;
   int ngroupedges = 0;
   int nodeshift;
   int curnnodes;
   int nvarnodestoadd;
   int internodeid;
   int nsymvars;
   int nsymedges;
   int first;
   int second;
   int color;
   int e;
   int i;
   int j;
 
   assert( scip != NULL );
   assert( graph1 != NULL );
   assert( graph2 != NULL );
   assert( SG != NULL || determinesize );
   assert( nnodes != NULL );
   assert( nedges != NULL );
   assert( degrees != NULL );
   assert( maxdegrees != NULL );
   assert( colors != NULL );
   assert( ncolors != NULL );
   assert( nusedvars != NULL );
   assert( ! determinesize || nnodesfromgraph1 != NULL );
   assert( success != NULL );
 
   *success = FALSE;
   if ( determinesize )
   {
      *degrees = NULL;
      *colors = NULL;
      *maxdegrees = 0;
      *ncolors = 0;
   }
 
   /* graphs cannot be symmetric */
   if ( SCIPgetSymgraphNEdges(graph1) != SCIPgetSymgraphNEdges(graph2)
      || SCIPgetSymgraphNVars(graph1) != SCIPgetSymgraphNVars(graph2) )
      return SCIP_OKAY;
 
   /* collect basic information from symmetry detection graph */
   nsymvars = SCIPgetSymgraphNVars(graph1);
   nsymedges = SCIPgetSymgraphNEdges(graph1);
   symtype = SCIPgetSymgraphSymtype(graph1);
   switch ( symtype )
   {
   case SYM_SYMTYPE_PERM:
      nvarnodestoadd = nsymvars;
      break;
   default:
      assert( symtype == SYM_SYMTYPE_SIGNPERM );
      nvarnodestoadd = 2 * nsymvars;
   } /*lint !e788*/
 
   /* find the variables that are contained in an edge */
   SCIP_CALL( SCIPallocClearBufferArray(scip, &nvarused1, nvarnodestoadd) );
   SCIP_CALL( SCIPallocClearBufferArray(scip, &nvarused2, nvarnodestoadd) );
   SCIP_CALL( SCIPallocBufferArray(scip, &varlabel, nvarnodestoadd) );
 
   for (e = 0; e < nsymedges; ++e)
   {
      first = SCIPgetSymgraphEdgeFirst(graph1, e);
      second = SCIPgetSymgraphEdgeSecond(graph1, e);
      if ( first < 0 )
         nvarused1[-first - 1] += 1;
      if ( second < 0 )
         nvarused1[-second - 1] += 1;
 
      first = SCIPgetSymgraphEdgeFirst(graph2, e);
      second = SCIPgetSymgraphEdgeSecond(graph2, e);
      if ( first < 0 )
         nvarused2[-first - 1] += 1;
      if ( second < 0 )
         nvarused2[-second - 1] += 1;
   }
 
   for (j = 0; j < nvarnodestoadd; ++j)
   {
      /* graphs cannot be identical */
      if ( nvarused1[j] != nvarused2[j] )
      {
         SCIPfreeBufferArray(scip, &varlabel);
         SCIPfreeBufferArray(scip, &nvarused2);
         SCIPfreeBufferArray(scip, &nvarused1);
 
         return SCIP_OKAY;
      }
 
      /* relabel variables by restricting to variables used in constraint (or their negation) */
      if ( nvarused1[j] > 0 || nvarused1[j % SCIPgetSymgraphNVars(graph1)] > 0 )
         varlabel[j] = nusdvars++;
      else
         varlabel[j] = -1;
   }
 
   /* possibly find number of nodes in sassy graph and allocate memory for dynamic array */
   if ( determinesize )
   {
      SCIP_CALL( SCIPensureBlockMemoryArray(scip, degrees, maxdegrees,
            SCIPgetSymgraphNNodes(graph1) + SCIPgetSymgraphNNodes(graph2) + 2 * nusdvars + 100) );
      SCIP_CALL( SCIPensureBlockMemoryArray(scip, colors, ncolors,
            SCIPgetSymgraphNNodes(graph1) + SCIPgetSymgraphNNodes(graph2) + 2 * nusdvars + 100) );
 
      *nnodes = 0;
      *nedges = 0;
   }
   else
   {
      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &pos, *nnodes) );
 
      /* add nodes for variables and remaining (axuiliary) nodes in graph */
      for (j = 0; j < *nnodes; ++j)
      {
         SG->d[j] = (*degrees)[j];
         SG->v[j] = (size_t) (unsigned) edgebegincnt;
         pos[j] = edgebegincnt;
         edgebegincnt += (*degrees)[j];
      }
   }
 
   /* determine grouping depending on the number of rhs vs. variables */
   groupByConstraints = SCIPgetSymgraphNConsnodes(graph1) < SCIPgetSymgraphNVars(graph1);
 
   /* allocate arrays to collect edges to be grouped */
   SCIP_CALL( SCIPallocBufferArray(scip, &groupfirsts, nsymedges) );
   SCIP_CALL( SCIPallocBufferArray(scip, &groupseconds, nsymedges) );
   SCIP_CALL( SCIPallocBufferArray(scip, &groupcolors, nsymedges) );
 
   /* collect information or generate graphs, we shift the node indices of the second graph when adding them to G */
   nodeshift = 0;
   for (i = 0; i < 2; ++i)
   {
      curnnodes = 0;
      symgraph = i == 0 ? graph1 : graph2;
      ngroupedges = 0;
 
      /* possibly add nodes for variables and remaining nodes, each variable gets a unique color */
      if ( determinesize )
      {
         /* add nodes for variables */
         for (j = 0; j < nvarnodestoadd; ++j)
         {
            if ( varlabel[j] >= 0 )
            {
               SCIP_CALL( SCIPensureBlockMemoryArray(scip, degrees, maxdegrees, *nnodes + 1) );
               (*degrees)[nodeshift + varlabel[j]] = 0;
               (*colors)[nodeshift + varlabel[j]] = SCIPgetSymgraphVarnodeColor(symgraph, j);
               ++(*nnodes);
               ++curnnodes;
            }
         }
 
         /* add nodes for remaining nodes of graph */
         for (j = 0; j < SCIPgetSymgraphNNodes(symgraph); ++j)
         {
            SCIP_CALL( SCIPensureBlockMemoryArray(scip, degrees, maxdegrees, *nnodes + 1) );
            (*degrees)[nodeshift + nusdvars + j] = 0;
            (*colors)[nodeshift + nusdvars + j] = SCIPgetSymgraphNodeColor(symgraph, j);
            ++(*nnodes);
            ++curnnodes;
         }
      }
      else
      {
         /* increase counter of nodes */
         for (j = 0; j < nvarnodestoadd; ++j)
         {
            if ( varlabel[j] >= 0 )
               ++curnnodes;
         }
         curnnodes += SCIPgetSymgraphNNodes(symgraph);
      }
 
      /* loop through all edges of the symmetry detection graph and either get degrees of nodes or add edges */
      internodeid = nodeshift + curnnodes;
      for (e = 0; e < nsymedges; ++e)
      {
         first = SCIPgetSymgraphEdgeFirst(symgraph, e);
         second = SCIPgetSymgraphEdgeSecond(symgraph, e);
 
         /* get the first and second node in edge (corrected by variable shift) */
         if ( first < 0 )
            first = varlabel[-first - 1];
         else
            first = nusdvars + first;
         if ( second < 0 )
            second = varlabel[-second - 1];
         else
            second = nusdvars + second;
 
         /* check whether edge is used for grouping */
         if ( ! SCIPhasGraphUniqueEdgetype(symgraph) && isEdgeGroupable(symgraph, e, groupByConstraints) )
         {
            /* store edge, first becomes the cons or var node */
            comparetype = groupByConstraints ? SYM_NODETYPE_CONS : SYM_NODETYPE_VAR;
 
            if ( SCIPgetSymgraphNodeType(symgraph, SCIPgetSymgraphEdgeFirst(symgraph, e)) == comparetype )
            {
               groupfirsts[ngroupedges] = nodeshift + first;
               groupseconds[ngroupedges] = nodeshift + second;
            }
            else
            {
               groupfirsts[ngroupedges] = nodeshift + second;
               groupseconds[ngroupedges] = nodeshift + first;
            }
            groupcolors[ngroupedges++] = nusdvars + SCIPgetSymgraphEdgeColor(symgraph, e);
         }
         else
         {
            /* immediately add edge or increase degrees */
            assert(0 <= first && first < *nnodes);
            assert(0 <= second && second < *nnodes);
 
            /* possibly split edge if it is colored */
            if ( ! SCIPhasGraphUniqueEdgetype(symgraph) && SCIPisSymgraphEdgeColored(symgraph, e) )
            {
               if ( determinesize )
               {
                  SCIP_CALL( SCIPensureBlockMemoryArray(scip, degrees, maxdegrees, nodeshift + internodeid + 1) );
                  SCIP_CALL( SCIPensureBlockMemoryArray(scip, colors, ncolors, nodeshift + internodeid + 1) );
 
                  ++(*degrees)[nodeshift + first];
                  ++(*degrees)[nodeshift + second];
                  (*degrees)[internodeid] = 2;
 
                  color = SCIPgetSymgraphEdgeColor(symgraph, e);
                  (*colors)[internodeid] = nusdvars + color;
 
                  ++(*nnodes);
                  *nedges += 2;
               }
               else
               {
                  assert( internodeid < *nnodes );
 
                  SG->e[pos[internodeid]++] = nodeshift + first;
                  SG->e[pos[internodeid]++] = nodeshift + second;
                  SG->e[pos[nodeshift + first]++] = internodeid;
                  SG->e[pos[nodeshift + second]++] = internodeid;
 
                  assert( internodeid == *nnodes - 1
                     || pos[internodeid] <= (int) SG->v[internodeid+1] );
                  assert( nodeshift + first == *nnodes - 1
                     || pos[nodeshift + first] <= (int) SG->v[nodeshift+first+1] );
                  assert( nodeshift + second == *nnodes - 1 ||
                     pos[nodeshift + second] <= (int) SG->v[nodeshift+second+1] );
               }
               ++internodeid;
               ++curnnodes;
            }
            else
            {
               if ( determinesize )
               {
                  ++(*degrees)[nodeshift + first];
                  ++(*degrees)[nodeshift + second];
                  ++(*nedges);
               }
               else
               {
                  SG->e[pos[nodeshift + first]++] = nodeshift + second;
                  SG->e[pos[nodeshift + second]++] = nodeshift + first;
 
                  assert( nodeshift+first == *nnodes - 1 || pos[nodeshift+first] <= (int) SG->v[nodeshift+first+1] );
                  assert( nodeshift+second == *nnodes - 1 || pos[nodeshift+second] <= (int) SG->v[nodeshift+second+1] );
               }
            }
         }
      }
 
      /* possibly add groupable edges */
      if ( ngroupedges > 0 )
      {
         int firstidx = 0;
         int firstnodeidx;
         int naddednodes;
         int naddededges;
 
         /* sort edges according to their first nodes */
         SCIPsortIntIntInt(groupfirsts, groupseconds, groupcolors, ngroupedges);
         firstnodeidx = groupfirsts[0];
 
         for (j = 1; j < ngroupedges; ++j)
         {
            /* if a new first node has been found, group the edges of the previous first node; ignoring the last group */
            if ( groupfirsts[j] != firstnodeidx )
            {
               SCIP_CALL( addOrDetermineEffectOfGroupedEdges(scip, SG, pos, determinesize, &internodeid,
                     degrees, maxdegrees, colors, ncolors, nnodes, nedges, firstnodeidx,
                     &groupseconds[firstidx], &groupcolors[firstidx], j - firstidx, &naddednodes, &naddededges) );
 
               firstidx = j;
               firstnodeidx = groupfirsts[j];
 
               if ( determinesize )
               {
                  *nnodes += naddednodes;
                  *nedges += naddededges;
               }
               curnnodes += naddednodes;
            }
         }
 
         /* process the last group */
         SCIP_CALL( addOrDetermineEffectOfGroupedEdges(scip, SG, pos, determinesize, &internodeid,
               degrees, maxdegrees, colors, ncolors, nnodes, nedges, firstnodeidx,
               &groupseconds[firstidx], &groupcolors[firstidx], ngroupedges - firstidx, &naddednodes, &naddededges) );
 
         if ( determinesize )
         {
            *nnodes += naddednodes;
            *nedges += naddededges;
         }
         curnnodes += naddednodes;
      }
 
      /* for signed permutation, also add edges connecting a variable and its negation */
      if ( SCIPgetSymgraphSymtype(graph1) == SYM_SYMTYPE_SIGNPERM )
      {
         if ( determinesize )
         {
            for (j = 0; j < nusdvars; ++j)
               ++(*degrees)[nodeshift + j];
            (*nedges) += nusdvars / 2;
         }
         else
         {
            for (j = 0; j < nusdvars/2; ++j)
            {
               SG->e[pos[nodeshift+j]++] = nodeshift + j + nusdvars/2;
               SG->e[pos[nodeshift + j + nusdvars/2]++] = nodeshift + j;
 
               assert( pos[nodeshift+j] <= (int) SG->v[nodeshift+j+1] );
               assert( nodeshift+j+nusdvars/2 == *nnodes - 1
                  || pos[nodeshift+j+nusdvars/2] <= (int) SG->v[nodeshift+j+nusdvars/2+1] );
            }
         }
      }
      nodeshift = curnnodes;
 
      /* possibly store number of nodes arising from first graph */
      if ( determinesize && i == 0 )
         *nnodesfromgraph1 = *nnodes;
   }
 
   /* add dummy node */
   if ( determinesize )
   {
      SCIP_CALL( SCIPensureBlockMemoryArray(scip, degrees, maxdegrees, *nnodes + 1) );
      SCIP_CALL( SCIPensureBlockMemoryArray(scip, colors, ncolors, *nnodes + 1) );
 
      ++(*nnodes);
      for (j = 0; j < *nnodes - 1; ++j)
         ++(*degrees)[j];
      (*degrees)[*nnodes - 1] = *nnodes - 1;
      (*nedges) += *nnodes - 1;
      (*colors)[*nnodes - 1] = 8;
   }
   else
   {
      for (j = 0; j < *nnodes - 1; ++j)
      {
         SG->e[pos[j]++] = *nnodes - 1;
         SG->e[pos[*nnodes-1]++] = j;
      }
      SCIPfreeBlockMemoryArray(scip, &pos, *nnodes);
   }
 
   SCIPfreeBufferArray(scip, &groupcolors);
   SCIPfreeBufferArray(scip, &groupseconds);
   SCIPfreeBufferArray(scip, &groupfirsts);
 
   SCIPfreeBufferArray(scip, &varlabel);
   SCIPfreeBufferArray(scip, &nvarused2);
   SCIPfreeBufferArray(scip, &nvarused1);
 
   *success = TRUE;
   if ( determinesize )
      *nusedvars = nusdvars;
 
   return SCIP_OKAY;
}
 
/** return whether symmetry can be computed */
SCIP_Bool SYMcanComputeSymmetry(void)
{
   return TRUE;
}
 
/** nauty/traces version string */
#ifdef NAUTY
static const char nautyname[] = {'N', 'a', 'u', 't', 'y', ' ', NAUTYVERSIONID/10000 + '0', '.', (NAUTYVERSIONID%10000)/1000 + '0', '.', (NAUTYVERSIONID%1000)/10 + '0', '\0'};
#else
static const char nautyname[] = {'T', 'r', 'a', 'c', 'e', 's', ' ', NAUTYVERSIONID/10000 + '0', '.', (NAUTYVERSIONID%10000)/1000 + '0', '.', (NAUTYVERSIONID%1000)/10 + '0', '\0'};
#endif
 
/** return name of external program used to compute generators */
const char* SYMsymmetryGetName(void)
{
   return nautyname;
}
 
/** return description of external program used to compute generators */
const char* SYMsymmetryGetDesc(void)
{
#ifdef NAUTY
   return "Computing Graph Automorphism Groups by Brendan D. McKay (https://users.cecs.anu.edu.au/~bdm/nauty/)";
#else
   return "Computing Graph Automorphism Groups by Adolfo Piperno (https://pallini.di.uniroma1.it/)";
#endif
}
 
/** return name of additional external program used for computing symmetries */
const char* SYMsymmetryGetAddName(void)
{
   return NULL;
}
 
/** return description of additional external program used to compute symmetries */
const char* SYMsymmetryGetAddDesc(void)
{
   return NULL;
}
 
/** compute generators of symmetry group */
SCIP_RETCODE SYMcomputeSymmetryGenerators(
   SCIP*                 scip,               /**< SCIP pointer */
   int                   maxgenerators,      /**< maximal number of generators constructed (= 0 if unlimited) */
   SYM_GRAPH*            symgraph,           /**< symmetry detection graph */
   int*                  nperms,             /**< pointer to store number of permutations */
   int*                  nmaxperms,          /**< pointer to store maximal number of permutations (needed for freeing storage) */
   int***                perms,              /**< pointer to store permutation generators as (nperms x npermvars) matrix */
   SCIP_Real*            log10groupsize,     /**< pointer to store size of group */
   SCIP_Real*            symcodetime         /**< pointer to store the time for symmetry code */
   )
{
   SCIP_Real oldtime;
   int nnodes;
   int nedges;
   int* degrees;
   int maxdegrees;
   int* colors;
   int ncolors;
   SCIP_Bool success;
   int v;
 
   /* nauty data structures */
   sparsegraph SG;
   int* lab;
   int* ptn;
   int* orbits;
 
#ifdef NAUTY
   DEFAULTOPTIONS_SPARSEGRAPH(options);
   statsblk stats;
#else
   static DEFAULTOPTIONS_TRACES(options);
   TracesStats stats;
#endif
 
   /* init */
   *nperms = 0;
   *nmaxperms = 0;
   *perms = NULL;
   *log10groupsize = 0;
   *symcodetime = 0.0;
 
   /* init options */
#ifdef NAUTY
   /* init callback functions for nauty (accumulate the group generators found by nauty) */
   options.writeautoms = FALSE;
   options.userautomproc = nautyhook;
   options.defaultptn = FALSE; /* use color classes */
   options.usernodeproc = nautyterminationhook;
#else
   /* init callback functions for traces (accumulate the group generators found by traces) */
   options.writeautoms = FALSE;
   options.userautomproc = traceshook;
   options.defaultptn = FALSE; /* use color classes */
#endif
 
   oldtime = SCIPgetSolvingTime(scip);
 
   /* determine the number of nodes and edges */
   SCIP_CALL( createOrDetermineSizeGraph(scip, symgraph, TRUE, NULL, &nnodes, &nedges,
         &degrees, &maxdegrees, &colors, &ncolors, &success) );
 
   if ( ! success )
   {
      SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, 0,
         "Stopped symmetry computation: Symmetry graph would become too large.\n");
      return SCIP_OKAY;
   }
 
   /* init graph */
   SG_INIT(SG);
 
   SG_ALLOC(SG, (size_t) nnodes, 2 * (size_t)(unsigned) nedges, "malloc"); /*lint !e647*//*lint !e774*//*lint !e571*/
 
   SG.nv = nnodes;                   /* number of nodes */
   SG.nde = (size_t) (unsigned) (2 * nedges);   /* number of directed edges */
 
   /* add the nodes for linear and nonlinear constraints to the graph */
   SCIP_CALL( createOrDetermineSizeGraph(scip, symgraph, FALSE, &SG, &nnodes, &nedges,
         &degrees, &maxdegrees, &colors, &ncolors, &success) );
 
   SCIPfreeBlockMemoryArray(scip, &degrees, maxdegrees);
 
   /* memory allocation for nauty/traces */
   SCIP_CALL( SCIPallocBufferArray(scip, &lab, nnodes) );
   SCIP_CALL( SCIPallocBufferArray(scip, &ptn, nnodes) );
   SCIP_CALL( SCIPallocBufferArray(scip, &orbits, nnodes) );
 
   /* fill in array with colors for variables */
   for (v = 0; v < nnodes; ++v)
      lab[v] = v;
 
   /* sort nodes according to colors */
   SCIPsortIntInt(colors, lab, nnodes);
 
   /* set up ptn marking new colors */
   for (v = 0; v < nnodes; ++v)
   {
      if ( v < nnodes-1 && colors[v] == colors[v+1] )
         ptn[v] = 1;  /* color class does not end */
      else
         ptn[v] = 0;  /* color class ends */
   }
 
   SCIPdebugMsg(scip, "Symmetry detection graph has %d nodes.\n", nnodes);
 
   data_.scip = scip;
   data_.npermvars = SCIPgetSymgraphNVars(symgraph);
   data_.nperms = 0;
   data_.nmaxperms = 0;
   data_.maxgenerators = maxgenerators;
   data_.perms = NULL;
   data_.symtype = SCIPgetSymgraphSymtype(symgraph);
   data_.restricttovars = TRUE;
   SCIP_CALL( SCIPgetIntParam(scip, "propagating/symmetry/nautymaxlevel", &data_.maxlevel) );
 
   /* call nauty/traces */
#ifdef NAUTY
   sparsenauty(&SG, lab, ptn, orbits, &options, &stats, NULL);
#else
   Traces(&SG, lab, ptn, orbits, &options, &stats, NULL);
#endif
   *symcodetime = SCIPgetSolvingTime(scip) - oldtime;
 
   SCIPfreeBufferArray(scip, &orbits);
   SCIPfreeBufferArray(scip, &ptn);
   SCIPfreeBufferArray(scip, &lab);
 
   SCIPfreeBlockMemoryArray(scip, &colors, ncolors);
 
   SG_FREE(SG); /*lint !e774*/
 
   /* prepare return values */
   if ( data_.nperms > 0 )
   {
      *perms = data_.perms;
      *nperms = data_.nperms;
      *nmaxperms = data_.nmaxperms;
   }
   else
   {
      assert( data_.perms == NULL );
      assert( data_.nmaxperms == 0 );
 
      *perms = NULL;
      *nperms = 0;
      *nmaxperms = 0;
   }
 
   /* determine log10 of symmetry group size */
   *log10groupsize = log10(stats.grpsize1 * pow(10.0, (SCIP_Real) stats.grpsize2));
 
   return SCIP_OKAY;
}
 
/** returns whether two given graphs are identical */
SCIP_Bool SYMcheckGraphsAreIdentical(
   SCIP*                 scip,               /**< SCIP pointer */
   SYM_SYMTYPE           symtype,            /**< type of symmetries to be checked */
   SYM_GRAPH*            G1,                 /**< first graph */
   SYM_GRAPH*            G2                  /**< second graph */
   )
{
   int nnodes;
   int nedges;
   int* degrees;
   int maxdegrees;
   int* colors;
   int ncolors;
   int nusedvars;
   SCIP_Bool success;
   int v;
   int nnodesfromG1;
 
   assert( scip != NULL );
   assert( G1 != NULL );
   assert( G2 != NULL );
 
   /* some simple checks */
   if ( G1->nnodes != G2->nnodes ||  G1->nopnodes != G2->nopnodes || G1->nvalnodes != G2->nvalnodes
      || G1->nconsnodes != G2->nconsnodes || G1->nedges != G2->nedges )
      return FALSE;
 
   SCIP_CALL_ABORT( createOrDetermineSizeGraphCheck(scip, G1, G2, TRUE, NULL, &nnodes, &nedges, &degrees, &maxdegrees,
         &colors, &ncolors, &nusedvars, &nnodesfromG1, &success) );
 
   if ( ! success )
   {
      SCIPfreeBlockMemoryArrayNull(scip, &degrees, maxdegrees);
      SCIPfreeBlockMemoryArrayNull(scip, &colors, ncolors);
 
      return FALSE;
   }
 
   /* nauty data structures */
   sparsegraph SG;
   int* lab;
   int* ptn;
   int* orbits;
 
#ifdef NAUTY
   DEFAULTOPTIONS_SPARSEGRAPH(options);
   statsblk stats;
#else
   static DEFAULTOPTIONS_TRACES(options);
   TracesStats stats;
#endif
 
   /* init options */
#ifdef NAUTY
   /* init callback functions for nauty (accumulate the group generators found by nauty) */
   options.writeautoms = FALSE;
   options.userautomproc = nautyhook;
   options.defaultptn = FALSE; /* use color classes */
#else
   /* init callback functions for traces (accumulate the group generators found by traces) */
   options.writeautoms = FALSE;
   options.userautomproc = traceshook;
   options.defaultptn = FALSE; /* use color classes */
#endif
 
   /* init graph */
   SG_INIT(SG);
 
   SG_ALLOC(SG, (size_t) nnodes, 2 * (size_t)(unsigned) nedges, "malloc"); /*lint !e647*//*lint !e774*//*lint !e571*/
 
   SG.nv = nnodes;                   /* number of nodes */
   SG.nde = (size_t) (unsigned) (2 * nedges);   /* number of directed edges */
 
   /* add the nodes for linear and nonlinear constraints to the graph */
   SCIP_CALL_ABORT( createOrDetermineSizeGraphCheck(scip, G1, G2, FALSE, &SG, &nnodes, &nedges, &degrees, &maxdegrees,
         &colors, &ncolors, &nusedvars, NULL, &success) );
   assert( success );
 
   SCIPfreeBlockMemoryArray(scip, &degrees, maxdegrees);
 
#ifdef SCIP_DISABLED_CODE
   /* print information about sparsegraph */
   SCIPinfoMessage(scip, NULL, "number of nodes: %d\n", SG.nv);
   SCIPinfoMessage(scip, NULL, "number of (directed) edges: %lu\n", SG.nde);
   SCIPinfoMessage(scip, NULL, "degrees\n");
   for (v = 0; v < SG.nv; ++v)
   {
      SCIPinfoMessage(scip, NULL, "node %d: %d\n", v, SG.d[v]);
   }
   SCIPinfoMessage(scip, NULL, "colors\n");
   for (v = 0; v < SG.nv; ++v)
   {
      SCIPinfoMessage(scip, NULL, "node %d: %d\n", v, colors[v]);
   }
   SCIPinfoMessage(scip, NULL, "edges\n");
   for (v = 0; v < SG.nv; ++v)
   {
      for (int w = 0; w < SG.d[v]; ++w)
      {
         SCIPinfoMessage(scip, NULL, "(%d,%d)\n", v, SG.e[SG.v[v] + w]);
      }
   }
#endif
 
   /* memory allocation for nauty/traces */
   SCIP_CALL_ABORT( SCIPallocBufferArray(scip, &lab, nnodes) );
   SCIP_CALL_ABORT( SCIPallocBufferArray(scip, &ptn, nnodes) );
   SCIP_CALL_ABORT( SCIPallocBufferArray(scip, &orbits, nnodes) );
 
   /* fill in array with colors for variables */
   for (v = 0; v < nnodes; ++v)
      lab[v] = v;
 
   /* sort nodes according to colors */
   SCIPsortIntInt(colors, lab, nnodes);
 
   /* set up ptn marking new colors */
   for (v = 0; v < nnodes; ++v)
   {
      if ( v < nnodes-1 && colors[v] == colors[v+1] )
         ptn[v] = 1;  /* color class does not end */
      else
         ptn[v] = 0;  /* color class ends */
   }
 
#ifdef SCIP_DISABLED_CODE
   /* print further information about sparsegraph */
   SCIPinfoMessage(scip, NULL, "lab (and ptn):\n");
   for (v = 0; v < SG.nv; ++v)
   {
      SCIPinfoMessage(scip, NULL, "%d (%d)\n", lab[v], ptn[v]);
   }
#endif
 
   /* compute automorphisms */
   data_.scip = scip;
   data_.npermvars = SCIPgetSymgraphNVars(G1);
   data_.nperms = 0;
   data_.nmaxperms = 0;
   data_.maxgenerators = 0;
   data_.perms = NULL;
   data_.symtype = symtype;
   data_.restricttovars = FALSE;
   SCIP_CALL( SCIPgetIntParam(scip, "propagating/symmetry/nautymaxlevel", &data_.maxlevel) ); /*lint !e641*//*lint !e732*/
 
   /* call nauty/traces */
#ifdef NAUTY
   sparsenauty(&SG, lab, ptn, orbits, &options, &stats, NULL);
#else
   Traces(&SG, lab, ptn, orbits, &options, &stats, NULL);
#endif
 
   SCIPfreeBufferArray(scip, &orbits);
   SCIPfreeBufferArray(scip, &ptn);
   SCIPfreeBufferArray(scip, &lab);
 
   SCIPfreeBlockMemoryArray(scip, &colors, ncolors);
 
   SG_FREE(SG); /*lint !e774*/
 
   /* G1 and G2 cannot be isomorphic */
   if ( data_.nperms == 0 )
      return FALSE;
 
   success = FALSE;
   for (int p = 0; p < data_.nperms; ++p)
   {
      for (int i = 0; i < nnodesfromG1; ++i)
      {
         if ( data_.perms[p][i] >= nnodesfromG1 )
         {
            success = TRUE;
            break;
         }
      }
   }
 
   /* free memory */
   for (int p = 0; p < data_.nperms; ++p)
   {
      SCIPfreeBlockMemoryArray(scip, &data_.perms[p], nnodes);
   }
   SCIPfreeBlockMemoryArrayNull(scip, &data_.perms, data_.nmaxperms);
 
   SG_FREE(SG); /*lint !e774*/
 
   return success;
}