~ubuntu-branches/ubuntu/lucid/igraph/lucid

« back to all changes in this revision

Viewing changes to src/bliss_graph.cc

Committer: Bazaar Package Importer
Author(s): Mathieu Malaterre
Date: 2009-11-16 18:12:42 UTC
Revision ID: james.westby@ubuntu.com-20091116181242-mzv9p5fz9uj57xd1

Tags: upstream-0.5.3

Import upstream version 0.5.3

files added:

AUTHORS

COPYING

ChangeLog

INSTALL

Makefile.am

Makefile.in

NEWS

README

TODO

aclocal.m4

config.guess

config.h.in

config.sub

configure

configure.in

depcomp

doc/Makefile.in

doc/book

doc/book/Makefile.in

doc/igraph-docs.info

doc/igraph.3

examples

examples/simple

examples/simple/LINKS.NET

examples/simple/adjlist.c

examples/simple/ak-4102.max

examples/simple/bellman_ford.c

examples/simple/bellman_ford.out

examples/simple/cattributes.c

examples/simple/cattributes.out

examples/simple/cattributes2.c

examples/simple/cattributes2.out

examples/simple/d_indheap.c

examples/simple/d_indheap.out

examples/simple/dijkstra.c

examples/simple/dijkstra.out

examples/simple/dot.c

examples/simple/dot.out

examples/simple/dqueue.c

examples/simple/eigenvector_centrality.c

examples/simple/eigenvector_centrality.out

examples/simple/flow.c

examples/simple/foreign.c

examples/simple/foreign.out

examples/simple/gml.c

examples/simple/gml.out

examples/simple/graphml-hsa05010.xml

examples/simple/graphml.c

examples/simple/graphml.out

examples/simple/heap.c

examples/simple/igraph_add_edges.c

examples/simple/igraph_add_edges.out

examples/simple/igraph_add_vertices.c

examples/simple/igraph_adjacency.c

examples/simple/igraph_array.c

examples/simple/igraph_array.out

examples/simple/igraph_atlas.c

examples/simple/igraph_atlas.out

examples/simple/igraph_average_path_length.c

examples/simple/igraph_barabasi_game.c

examples/simple/igraph_bfs.c

examples/simple/igraph_bfs.out

examples/simple/igraph_bipartite_create.c

examples/simple/igraph_bipartite_create.out

examples/simple/igraph_bipartite_projection.c

examples/simple/igraph_cliques.c

examples/simple/igraph_cliques.out

examples/simple/igraph_cocitation.c

examples/simple/igraph_cocitation.out

examples/simple/igraph_community_edge_betweenness.c

examples/simple/igraph_community_edge_betweenness.out

examples/simple/igraph_community_fastgreedy.c

examples/simple/igraph_community_fastgreedy.out

examples/simple/igraph_community_label_propagation.c

examples/simple/igraph_community_label_propagation.out

examples/simple/igraph_community_leading_eigenvector.c

examples/simple/igraph_community_leading_eigenvector.out

examples/simple/igraph_complementer.c

examples/simple/igraph_complementer.out

examples/simple/igraph_compose.c

examples/simple/igraph_compose.out

examples/simple/igraph_convergence_degree.c

examples/simple/igraph_convergence_degree.out

examples/simple/igraph_convex_hull.c

examples/simple/igraph_convex_hull.out

examples/simple/igraph_copy.c

examples/simple/igraph_create.c

examples/simple/igraph_decompose.c

examples/simple/igraph_decompose.out

examples/simple/igraph_degree.c

examples/simple/igraph_degree.out

examples/simple/igraph_degree_sequence_game.c

examples/simple/igraph_delete_edges.c

examples/simple/igraph_delete_vertices.c

examples/simple/igraph_diameter.c

examples/simple/igraph_diameter.out

examples/simple/igraph_difference.c

examples/simple/igraph_difference.out

examples/simple/igraph_disjoint_union.c

examples/simple/igraph_disjoint_union.out

examples/simple/igraph_edge_betweenness.c

examples/simple/igraph_edge_betweenness.out

examples/simple/igraph_empty.c

examples/simple/igraph_erdos_renyi_game.c

examples/simple/igraph_es_adj.c

examples/simple/igraph_es_adj.out

examples/simple/igraph_es_fromto.c

examples/simple/igraph_es_fromto.out

examples/simple/igraph_es_pairs.c

examples/simple/igraph_es_path.c

examples/simple/igraph_full.c

examples/simple/igraph_get_eid.c

examples/simple/igraph_get_eid.out

examples/simple/igraph_get_eids.c

examples/simple/igraph_get_eids.out

examples/simple/igraph_get_shortest_paths.c

examples/simple/igraph_get_shortest_paths.out

examples/simple/igraph_get_shortest_paths_dijkstra.c

examples/simple/igraph_get_shortest_paths_dijkstra.out

examples/simple/igraph_girth.c

examples/simple/igraph_grg_game.c

examples/simple/igraph_growing_random_game.c

examples/simple/igraph_hashtable.c

examples/simple/igraph_hashtable.out

examples/simple/igraph_i_cutheap.c

examples/simple/igraph_i_cutheap.out

examples/simple/igraph_i_layout_sphere.c

examples/simple/igraph_independent_sets.c

examples/simple/igraph_independent_sets.out

examples/simple/igraph_intersection.c

examples/simple/igraph_intersection.out

examples/simple/igraph_is_directed.c

examples/simple/igraph_is_loop.c

examples/simple/igraph_is_loop.out

examples/simple/igraph_is_multiple.c

examples/simple/igraph_is_multiple.out

examples/simple/igraph_lattice.c

examples/simple/igraph_layout_lgl.c

examples/simple/igraph_layout_merge.c

examples/simple/igraph_layout_reingold_tilford.c

examples/simple/igraph_layout_reingold_tilford.in

examples/simple/igraph_lcf.c

examples/simple/igraph_local_transitivity.c

examples/simple/igraph_mincut.c

examples/simple/igraph_mincut.out

examples/simple/igraph_minimum_spanning_tree.c

examples/simple/igraph_minimum_spanning_tree.out

examples/simple/igraph_motifs_randesu.c

examples/simple/igraph_motifs_randesu.out

examples/simple/igraph_neighbors.c

examples/simple/igraph_neighbors.out

examples/simple/igraph_pagerank.c

examples/simple/igraph_pagerank.out

examples/simple/igraph_preference_game.c

examples/simple/igraph_psumtree.c

examples/simple/igraph_read_graph_graphdb.c

examples/simple/igraph_read_graph_graphdb.out

examples/simple/igraph_read_graph_lgl-1.lgl

examples/simple/igraph_read_graph_lgl-2.lgl

examples/simple/igraph_read_graph_lgl-3.lgl

examples/simple/igraph_read_graph_lgl.c

examples/simple/igraph_read_graph_lgl.out

examples/simple/igraph_reciprocity.c

examples/simple/igraph_rewire.c

examples/simple/igraph_rewire.out

examples/simple/igraph_ring.c

examples/simple/igraph_set.c

examples/simple/igraph_set.out

examples/simple/igraph_similarity.c

examples/simple/igraph_similarity.out

examples/simple/igraph_simplify.c

examples/simple/igraph_simplify.out

examples/simple/igraph_small.c

examples/simple/igraph_small.out

examples/simple/igraph_star.c

examples/simple/igraph_strvector.c

examples/simple/igraph_strvector.out

examples/simple/igraph_to_undirected.c

examples/simple/igraph_to_undirected.out

examples/simple/igraph_topological_sorting.c

examples/simple/igraph_topological_sorting.out

examples/simple/igraph_transitivity.c

examples/simple/igraph_tree.c

examples/simple/igraph_trie.c

examples/simple/igraph_trie.out

examples/simple/igraph_union.c

examples/simple/igraph_union.out

examples/simple/igraph_vs_nonadj.c

examples/simple/igraph_vs_nonadj.out

examples/simple/igraph_vs_seq.c

examples/simple/igraph_vs_seq.out

examples/simple/igraph_vs_vector.c

examples/simple/igraph_vs_vector.out

examples/simple/igraph_weighted_adjacency.c

examples/simple/igraph_weighted_adjacency.out

examples/simple/igraph_write_graph_lgl.c

examples/simple/igraph_write_graph_pajek.c

examples/simple/igraph_write_graph_pajek.out

examples/simple/indheap.c

examples/simple/iso_b03_m1000.A00

examples/simple/karate.gml

examples/simple/lineendings.c

examples/simple/lineendings.out

examples/simple/matrix.c

examples/simple/matrix.out

examples/simple/matrix2.c

examples/simple/matrix2.out

examples/simple/pajek1.net

examples/simple/pajek2.net

examples/simple/pajek3.net

examples/simple/pajek4.net

examples/simple/random_seed.c

examples/simple/spinglass.c

examples/simple/spmatrix.c

examples/simple/spmatrix.out

examples/simple/stack.c

examples/simple/test.gxl

examples/simple/topology.c

examples/simple/topology.out

examples/simple/vector.c

examples/simple/vector.out

examples/simple/vector2.c

examples/simple/vector2.out

examples/simple/vector_ptr.c

examples/simple/walktrap.c

examples/simple/walktrap.out

igraph.pc

igraph.pc.in

igraph_Info.plist.in

include

include/arpack.h

include/array.h

include/attributes.h

include/dqueue.h

include/error.h

include/heap.h

include/igraph.h

include/igraph_pmt.h

include/igraph_pmt_off.h

include/interrupt.h

include/matrix.h

include/stack.h

include/types.h

include/vector.h

install-sh

interfaces

interfaces/Makefile.am

interfaces/Makefile.in

interfaces/R

interfaces/R/DESCRIPTION.in

interfaces/R/Makefile.in

interfaces/R/configure.in.in

interfaces/R/src

interfaces/R/src/Makevars.win.in

interfaces/shell

interfaces/shell/Makefile.am

interfaces/shell/Makefile.in

interfaces/shell/interface.c

ltmain.sh

missing

src/DensityGrid.cpp

src/DensityGrid.h

src/DensityGrid_3d.cpp

src/DensityGrid_3d.h

src/Makefile.am

src/Makefile.in

src/NetDataTypes.cpp

src/NetDataTypes.h

src/NetRoutines.cpp

src/NetRoutines.h

src/adjlist.c

src/arpack

src/arpack.c

src/arpack.h

src/arpack/dgetv0.c

src/arpack/dlaqrb.c

src/arpack/dmout.c

src/arpack/dnaitr.c

src/arpack/dnapps.c

src/arpack/dnaup2.c

src/arpack/dnaupd.c

src/arpack/dnconv.c

src/arpack/dneigh.c

src/arpack/dneupd.c

src/arpack/dngets.c

src/arpack/dsaitr.c

src/arpack/dsapps.c

src/arpack/dsaup2.c

src/arpack/dsaupd.c

src/arpack/dsconv.c

src/arpack/dseigt.c

src/arpack/dsesrt.c

src/arpack/dseupd.c

src/arpack/dsgets.c

src/arpack/dsortc.c

src/arpack/dsortr.c

src/arpack/dstatn.c

src/arpack/dstats.c

src/arpack/dstqrb.c

src/arpack/dvout.c

src/arpack/ivout.c

src/arpack/second.c

src/arpack_internal.h

src/array.c

src/array.h

src/array.pmt

src/atlas-edges.h

src/atlas.c

src/attributes.c

src/attributes.h

src/basic_query.c

src/bfgs.c

src/bipartite.c

src/blas

src/blas/dasum.c

src/blas/daxpy.c

src/blas/dcopy.c

src/blas/ddot.c

src/blas/dgemv.c

src/blas/dger.c

src/blas/dnrm2.c

src/blas/drot.c

src/blas/dscal.c

src/blas/dswap.c

src/blas/dtrmm.c

src/blas/idamax.c

src/blas/lsame.c

src/bliss.cc

src/bliss_bignum.hh

src/bliss_defs.hh

src/bliss_eqrefhash.cc

src/bliss_eqrefhash.hh

src/bliss_graph.cc

src/bliss_graph.hh

src/bliss_heap.cc

src/bliss_heap.hh

src/bliss_kqueue.hh

src/bliss_kstack.hh

src/bliss_orbit.cc

src/bliss_orbit.hh

src/bliss_partition.cc

src/bliss_partition.hh

src/bliss_timer.cc

src/bliss_timer.hh

src/bliss_utils.cc

src/bliss_utils.hh

src/cattributes.c

src/centrality.c

src/cliques.c

src/clustertool.cpp

src/cocitation.c

src/community.c

src/components.c

src/conversion.c

src/cores.c

src/dqueue.c

src/dqueue.h

src/dqueue.pmt

src/drl_Node.h

src/drl_Node_3d.h

src/drl_graph.cpp

src/drl_graph.h

src/drl_graph_3d.cpp

src/drl_graph_3d.h

src/drl_layout.cpp

src/drl_layout.h

src/drl_layout_3d.cpp

src/drl_layout_3d.h

src/drl_parse.cpp

src/drl_parse.h

src/error.c

src/error.h

src/evolver_cit.c

src/f2c

src/f2c.h

src/f2c/d_lg10.c

src/f2c/d_sign.c

src/f2c/etime_.c

src/f2c/i_dnnt.c

src/f2c/pow_dd.c

src/f2c/pow_di.c

src/f2c/s_cmp.c

src/f2c/s_copy.c

src/fast_community.c

src/flow.c

src/foreign-gml-lexer.c

src/foreign-gml-lexer.l

src/foreign-gml-parser.c

src/foreign-gml-parser.h

src/foreign-gml-parser.y

src/foreign-graphml.c

src/foreign-lgl-lexer.c

src/foreign-lgl-lexer.l

src/foreign-lgl-parser.c

src/foreign-lgl-parser.h

src/foreign-lgl-parser.y

src/foreign-ncol-lexer.c

src/foreign-ncol-lexer.l

src/foreign-ncol-parser.c

src/foreign-ncol-parser.h

src/foreign-ncol-parser.y

src/foreign-pajek-lexer.c

src/foreign-pajek-lexer.l

src/foreign-pajek-parser.c

src/foreign-pajek-parser.h

src/foreign-pajek-parser.y

src/foreign.c

src/forestfire.c

src/games.c

src/gengraph_box_list.cpp

src/gengraph_box_list.h

src/gengraph_definitions.h

src/gengraph_degree_sequence.cpp

src/gengraph_degree_sequence.h

src/gengraph_graph_molloy_hash.cpp

src/gengraph_graph_molloy_hash.h

src/gengraph_graph_molloy_optimized.cpp

src/gengraph_graph_molloy_optimized.h

src/gengraph_hash.h

src/gengraph_header.h

src/gengraph_mr-connected.cpp

src/gengraph_powerlaw.cpp

src/gengraph_powerlaw.h

src/gengraph_qsort.h

src/gengraph_random.cpp

src/gengraph_random.h

src/gengraph_vertex_cover.h

src/gml_tree.c

src/gml_tree.h

src/heap.c

src/heap.h

src/heap.pmt

src/igraph.h

src/igraph_buckets.c

src/igraph_grid.c

src/igraph_hashtable.c

src/igraph_heap.c

src/igraph_math.h

src/igraph_pmt.h

src/igraph_pmt_off.h

src/igraph_psumtree.c

src/igraph_set.c

src/igraph_stack.c

src/igraph_strvector.c

src/igraph_trie.c

src/interrupt.c

src/interrupt.h

src/iterators.c

src/lapack

src/lapack/dgeqr2.c

src/lapack/dlabad.c

src/lapack/dlacon.c

src/lapack/dlacpy.c

src/lapack/dladiv.c

src/lapack/dlae2.c

src/lapack/dlaev2.c

src/lapack/dlaexc.c

src/lapack/dlahqr.c

src/lapack/dlaln2.c

src/lapack/dlamch.c

src/lapack/dlange.c

src/lapack/dlanhs.c

src/lapack/dlanst.c

src/lapack/dlanv2.c

src/lapack/dlapy2.c

src/lapack/dlarf.c

src/lapack/dlarfg.c

src/lapack/dlarfx.c

src/lapack/dlarnv.c

src/lapack/dlartg.c

src/lapack/dlaruv.c

src/lapack/dlascl.c

src/lapack/dlaset.c

src/lapack/dlasr.c

src/lapack/dlasrt.c

src/lapack/dlassq.c

src/lapack/dlasy2.c

src/lapack/dorm2r.c

src/lapack/dsteqr.c

src/lapack/dtrevc.c

src/lapack/dtrexc.c

src/lapack/dtrsen.c

src/lapack/dtrsyl.c

src/lapack/xerbla.c

src/layout.c

src/math.c

src/matrix.c

src/matrix.h

src/matrix.pmt

src/measure_dynamics.c

src/memory.c

src/memory.h

src/motifs.c

src/operators.c

src/other.c

src/pottsmodel_2.cpp

src/pottsmodel_2.h

src/progress.c

src/random.c

src/random.h

src/revolver_cit.c

src/revolver_grow.c

src/revolver_ml_cit.c

src/spectral_properties.c

src/spmatrix.c

src/stack.h

src/stack.pmt

src/structural_properties.c

src/structure_generators.c

src/topology.c

src/type_indexededgelist.c

src/types.h

src/vector.c

src/vector.h

src/vector.pmt

src/vector_ptr.c

src/visitors.c

src/walktrap.cpp

src/walktrap_communities.cpp

src/walktrap_communities.h

src/walktrap_graph.cpp

src/walktrap_graph.h

src/walktrap_heap.cpp

src/walktrap_heap.h

src/zeroin.c

tests

tests/Makefile.am

tests/Makefile.in

tests/arpack.at

tests/atlocal.in

tests/attributes.at

tests/basic.at

tests/bipartite.at

tests/cliques.at

tests/community.at

tests/components.at

tests/conversion.at

tests/eigen.at

tests/flow.at

tests/foreign.at

tests/iterators.at

tests/layout.at

tests/motifs.at

tests/operators.at

tests/other.at

tests/package.m4

tests/structural_properties.at

tests/structure_generators.at

tests/testsuite

tests/testsuite.at

tests/topology.at

tests/types.at

tests/visitors.at

ylwrap

Show diffs side-by-side

added added

removed removed

src/bliss_graph.cc

This program is free software; you can redistribute it and/or modify

it under the terms of the GNU General Public License version 2

as published by the Free Software Foundation.

This program is distributed in the hope that it will be useful,

but WITHOUT ANY WARRANTY; without even the implied warranty of

MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

GNU General Public License for more details.

You should have received a copy of the GNU General Public License

along with this program; if not, write to the Free Software

Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.

/* FSF address fixed in the above notice on 1 Oct 2009 by Tamas Nepusz */

#include <stdio.h>

#include <assert.h>

#include <ctype.h>

#include <set>

#include <list>

#include <algorithm>

#include "bliss_defs.hh"

#include "bliss_timer.hh"

#include "bliss_graph.hh"

#include "bliss_partition.hh"

#include <limits.h> // INT_MAX, etc

extern bool bliss_verbose;

extern FILE *bliss_verbstr;

namespace igraph {

static const bool should_not_happen = false;

/*-------------------------------------------------------------------------

* Constructor and destructor routines for the abstract graph class

*-------------------------------------------------------------------------*/

AbstractGraph::AbstractGraph()

{

/* Initialize stuff */

first_path_labeling = 0;

first_path_labeling_inv = 0;

best_path_labeling = 0;

best_path_labeling_inv = 0;

first_path_automorphism = 0;

best_path_automorphism = 0;

//certificate = 0;

in_search = false;

}

AbstractGraph::~AbstractGraph()

{

if(first_path_labeling) {

free(first_path_labeling); first_path_labeling = 0; }

if(first_path_labeling_inv) {

free(first_path_labeling_inv); first_path_labeling_inv = 0; }

if(best_path_labeling) {

free(best_path_labeling); best_path_labeling = 0; }

if(best_path_labeling_inv) {

free(best_path_labeling_inv); best_path_labeling_inv = 0; }

if(first_path_automorphism) {

free(first_path_automorphism); first_path_automorphism = 0; }

if(best_path_automorphism) {

free(best_path_automorphism); best_path_automorphism = 0; }

//if(certificate) {

// free(certificate); certificate = 0; }

while(!long_prune_fixed.empty())

{

delete long_prune_fixed.back();

long_prune_fixed.pop_back();

}

while(!long_prune_mcrs.empty())

{

delete long_prune_mcrs.back();

long_prune_mcrs.pop_back();

}

/*-------------------------------------------------------------------------

* Routines for refinement to equitable partition

*-------------------------------------------------------------------------*/

void AbstractGraph::refine_to_equitable()

100

{

101

assert(p.splitting_queue.is_empty());

102

103

for(Cell *cell = p.first_cell; cell; cell = cell->next)

104

{

105

p.add_in_splitting_queue(cell);

106

}

107

108

return do_refine_to_equitable();

109

}

110

111

112

void AbstractGraph::refine_to_equitable(Cell *cell1)

113

{

114

DEBUG_ASSERT(cell1->length == 1);

115

116

#ifdef EXPENSIVE_CONSISTENCY_CHECKS

117

for(Cell *cell = p.first_cell; cell; cell = cell->next) {

118

assert(cell->in_splitting_queue == false);

119

assert(cell->in_neighbour_heap == false);

120

}

121

#endif

122

123

assert(p.splitting_queue.is_empty());

124

125

p.add_in_splitting_queue(cell1);

126

127

return do_refine_to_equitable();

128

}

129

130

131

void AbstractGraph::refine_to_equitable(Cell *cell1, Cell *cell2)

132

{

133

DEBUG_ASSERT(cell1->length == 1);

134

DEBUG_ASSERT(cell2->length == 1);

135

136

#ifdef EXPENSIVE_CONSISTENCY_CHECKS

137

for(Cell *cell = p.first_cell; cell; cell = cell->next) {

138

assert(cell->in_splitting_queue == false);

139

assert(cell->in_neighbour_heap == false);

140

}

141

#endif

142

143

assert(p.splitting_queue.is_empty());

144

145

p.add_in_splitting_queue(cell1);

146

p.add_in_splitting_queue(cell2);

147

148

return do_refine_to_equitable();

149

}

150

151

152

void AbstractGraph::do_refine_to_equitable()

153

{

154

assert(!p.splitting_queue.is_empty());

155

assert(neighbour_heap.is_empty());

156

157

eqref_hash.reset();

158

159

while(!p.splitting_queue.is_empty())

160

{

161

Cell *cell = p.splitting_queue.pop_front();

162

DEBUG_ASSERT(cell->in_splitting_queue);

163

cell->in_splitting_queue = false;

164

165

if(cell->length == 1)

166

{

167

if(in_search) {

168

if(first_path_automorphism) {

169

/* Build the (potential) automorphism on-the-fly */

170

assert(first_path_labeling_inv);

171

first_path_automorphism[first_path_labeling_inv[cell->first]] =

172

p.elements[cell->first];

173

}

174

if(best_path_automorphism)

175

{

176

/* Build the (potential) automorphism on-the-fly */

177

assert(best_path_labeling_inv);

178

best_path_automorphism[best_path_labeling_inv[cell->first]] =

179

p.elements[cell->first];

180

}

181

}

182

183

bool worse = split_neighbourhood_of_unit_cell(cell);

184

if(in_search && worse)

185

goto worse_exit;

186

}

187

else

188

{

189

split_neighbourhood_of_cell(cell);

190

}

191

}

192

193

eqref_worse_than_certificate = false;

194

return;

195

196

worse_exit:

197

/* Clear splitting_queue */

198

p.clear_splitting_queue();

199

eqref_worse_than_certificate = true;

200

return;

201

}

202

203

204

205

206

207

/*-------------------------------------------------------------------------

208

209

* Routines for handling the canonical labeling

210

211

*-------------------------------------------------------------------------*/

212

213

214

void AbstractGraph::update_labeling(unsigned int * const labeling)

215

{

216

const unsigned int N = get_nof_vertices();

217

unsigned int *ep = p.elements;

218

for(unsigned int i = 0; i < N; i++, ep++)

219

labeling[*ep] = i;

220

}

221

222

223

void AbstractGraph::update_labeling_and_its_inverse(unsigned int * const labeling,

224

unsigned int * const labeling_inv)

225

{

226

const unsigned int N = get_nof_vertices();

227

unsigned int *ep = p.elements;

228

unsigned int *clip = labeling_inv;

229

230

for(unsigned int i = 0; i < N; ) {

231

labeling[*ep] = i;

232

i++;

233

*clip = *ep;

234

ep++;

235

clip++;

236

}

237

}

238

239

240

241

/*-------------------------------------------------------------------------

242

243

* Routines for handling automorphisms

244

245

*-------------------------------------------------------------------------*/

246

247

248

void AbstractGraph::reset_permutation(unsigned int *perm)

249

{

250

const unsigned int N = get_nof_vertices();

251

for(unsigned int i = 0; i < N; i++, perm++)

252

*perm = i;

253

}

254

255

256

bool AbstractGraph::is_automorphism(unsigned int * const perm)

257

{

258

assert(should_not_happen);

259

return false;

260

}

261

262

263

264

265

266

/*-------------------------------------------------------------------------

267

268

* Long prune code

269

270

*-------------------------------------------------------------------------*/

271

272

void AbstractGraph::long_prune_init()

273

{

274

const unsigned int N = get_nof_vertices();

275

long_prune_temp.clear();

276

long_prune_temp.resize(N);

277

#ifdef DEBUG

278

for(unsigned int i = 0; i < N; i++)

279

assert(long_prune_temp[i] == false);

280

#endif

281

const unsigned int nof_fitting_in_max_mem =

282

(long_prune_options_max_mem * 1024 * 1024) / (((N * 2) / 8)+1);

283

long_prune_max_stored_autss = long_prune_options_max_stored_auts;

284

/* Had some problems with g++ in using (a<b)?a:b when constants involved,

285

so had to make this in a stupid way...*/

286

if(nof_fitting_in_max_mem < long_prune_options_max_stored_auts)

287

long_prune_max_stored_autss = nof_fitting_in_max_mem;

288

289

while(!long_prune_fixed.empty())

290

{

291

delete long_prune_fixed.back();

292

long_prune_fixed.pop_back();

293

}

294

while(!long_prune_mcrs.empty())

295

{

296

delete long_prune_mcrs.back();

297

long_prune_mcrs.pop_back();

298

}

299

for(unsigned int i = 0; i < long_prune_max_stored_autss; i++)

300

{

301

long_prune_fixed.push_back(new std::vector<bool>(N));

302

long_prune_mcrs.push_back(new std::vector<bool>(N));

303

}

304

long_prune_begin = 0;

305

long_prune_end = 0;

306

}

307

308

void AbstractGraph::long_prune_swap(const unsigned int i, const unsigned int j)

309

{

310

assert(long_prune_begin <= long_prune_end);

311

assert(long_prune_end - long_prune_end <= long_prune_max_stored_autss);

312

assert(i >= long_prune_begin);

313

assert(i < long_prune_end);

314

assert(j >= long_prune_begin);

315

assert(j < long_prune_end);

316

const unsigned int real_i = i % long_prune_max_stored_autss;

317

const unsigned int real_j = j % long_prune_max_stored_autss;

318

std::vector<bool> * tmp = long_prune_fixed[real_i];

319

long_prune_fixed[real_i] = long_prune_fixed[real_j];

320

long_prune_fixed[real_j] = tmp;

321

tmp = long_prune_mcrs[real_i];

322

long_prune_mcrs[real_i] = long_prune_mcrs[real_j];

323

long_prune_mcrs[real_j] = tmp;

324

}

325

326

std::vector<bool> &AbstractGraph::long_prune_get_fixed(const unsigned int index)

327

{

328

assert(long_prune_begin <= long_prune_end);

329

assert(long_prune_end - long_prune_end <= long_prune_max_stored_autss);

330

assert(index >= long_prune_begin);

331

assert(index < long_prune_end);

332

return *long_prune_fixed[index % long_prune_max_stored_autss];

333

}

334

335

std::vector<bool> &AbstractGraph::long_prune_get_mcrs(const unsigned int index)

336

{

337

assert(long_prune_begin <= long_prune_end);

338

assert(long_prune_end - long_prune_end <= long_prune_max_stored_autss);

339

assert(index >= long_prune_begin);

340

assert(index < long_prune_end);

341

return *long_prune_mcrs[index % long_prune_max_stored_autss];

342

}

343

344

345

void AbstractGraph::long_prune_add_automorphism(const unsigned int *aut)

346

{

347

if(long_prune_max_stored_autss == 0)

348

return;

349

350

const unsigned int N = get_nof_vertices();

351

352

#ifdef DEBUG

353

assert(long_prune_temp.size() == N);

354

for(unsigned int i = 0; i < N; i++)

355

assert(long_prune_temp[i] == false);

356

#endif

357

358

DEBUG_ASSERT(long_prune_fixed.size() == long_prune_mcrs.size());

359

assert(long_prune_begin <= long_prune_end);

360

assert(long_prune_end - long_prune_end <= long_prune_max_stored_autss);

361

if(long_prune_end - long_prune_begin == long_prune_max_stored_autss)

362

{

363

long_prune_begin++;

364

}

365

long_prune_end++;

366

std::vector<bool> &fixed = long_prune_get_fixed(long_prune_end-1);

367

std::vector<bool> &mcrs = long_prune_get_mcrs(long_prune_end-1);

368

369

for(unsigned int i = 0; i < N; i++)

370

{

371

fixed[i] = (aut[i] == i);

372

if(!long_prune_temp[i])

373

{

374

mcrs[i] = true;

375

unsigned int j = aut[i];

376

while(j != i)

377

{

378

assert(i <= j);

379

long_prune_temp[j] = true;

380

j = aut[j];

381

}

382

}

383

else

384

{

385

mcrs[i] = false;

386

}

387

long_prune_temp[i] = false;

388

}

389

390

391

#ifdef DEBUG

392

for(unsigned int i = 0; i < N; i++)

393

assert(long_prune_temp[i] == false);

394

#endif

395

}

396

397

398

/*-------------------------------------------------------------------------

399

400

* Routines for handling orbit information

401

402

*-------------------------------------------------------------------------*/

403

404

405

void AbstractGraph::update_orbit_information(Orbit &o, const unsigned int *p)

406

{

407

const unsigned int N = get_nof_vertices();

408

for(unsigned int i = 0; i < N; i++)

409

if(p[i] != i)

410

o.merge_orbits(i, p[i]);

411

}

412

413

414

415

416

417

/*-------------------------------------------------------------------------

418

419

* Print a permutation in cycle notation

420

421

*-------------------------------------------------------------------------*/

422

423

424

void AbstractGraph::print_permutation(FILE *fp, const unsigned int *perm)

425

{

426

const unsigned int N = get_nof_vertices();

427

for(unsigned int i = 0; i < N; i++) {

428

unsigned int j = perm[i];

429

if(j == i)

430

continue;

431

bool is_first = true;

432

while(j != i) {

433

if(j < i) {

434

is_first = false;

435

break;

436

}

437

j = perm[j];

438

}

439

if(!is_first)

440

continue;

441

fprintf(fp, "(%u,", i);

442

j = perm[i];

443

while(j != i) {

444

fprintf(fp, "%u", j);

445

j = perm[j];

446

if(j != i)

447

fprintf(fp, ",");

448

}

449

fprintf(fp, ")");

450

}

451

}

452

453

454

455

456

457

/*-------------------------------------------------------------------------

458

459

* The actual backtracking search

460

461

*-------------------------------------------------------------------------*/

462

463

464

typedef struct {

465

int split_element;

466

unsigned int split_cell_first;

467

unsigned int refinement_stack_size;

468

unsigned int certificate_index;

469

470

bool in_first_path;

471

bool in_best_path;

472

bool equal_to_first_path;

473

int cmp_to_best_path;

474

475

bool needs_long_prune;

476

unsigned int long_prune_begin;

477

std::set<unsigned int, std::less<unsigned int> > long_prune_redundant;

478

479

EqrefHash eqref_hash;

480

unsigned int subcertificate_length;

481

} LevelInfo;

482

483

484

485

typedef struct t_path_info {

486

unsigned int splitting_element;

487

unsigned int certificate_index;

488

unsigned int subcertificate_length;

489

EqrefHash eqref_hash;

490

} PathInfo;

491

492

493

void AbstractGraph::search(const bool canonical, Stats &stats)

494

{

495

const unsigned int N = get_nof_vertices();

496

497

const bool write_automorphisms = 0;

498

499

unsigned int all_same_level = UINT_MAX;

500

501

p.graph = this;

502

503

504

* Must be done!

505

506

remove_duplicate_edges();

507

508

509

* Reset search statistics

510

511

stats.group_size.assign(1);

512

stats.nof_nodes = 1;

513

stats.nof_leaf_nodes = 1;

514

stats.nof_bad_nodes = 0;

515

stats.nof_canupdates = 0;

516

stats.nof_generators = 0;

517

stats.max_level = 0;

518

519

if(first_path_labeling)

520

{

521

free(first_path_labeling);

522

first_path_labeling = 0;

523

}

524

if(first_path_labeling_inv)

525

{

526

free(first_path_labeling_inv);

527

first_path_labeling_inv = 0;

528

}

529

if(first_path_automorphism)

530

{

531

free(first_path_automorphism);

532

first_path_automorphism = 0;

533

}

534

535

if(best_path_labeling)

536

{

537

free(best_path_labeling);

538

best_path_labeling = 0;

539

}

540

if(best_path_labeling_inv)

541

{

542

free(best_path_labeling_inv);

543

best_path_labeling_inv = 0;

544

}

545

if(best_path_automorphism)

546

{

547

free(best_path_automorphism);

548

best_path_automorphism = 0;

549

}

550

551

if(N == 0)

552

return;

553

554

p.init(N);

555

neighbour_heap.init(N);

556

557

in_search = false;

558

559

p.level = 0;

560

561

Timer t1;

562

t1.start();

563

564

make_initial_equitable_partition();

565

566

#if defined(VERIFY_EQUITABLEDNESS)

567

assert(is_equitable());

568

#endif

569

570

t1.stop();

571

if(bliss_verbose) {

572

fprintf(bliss_verbstr, "Initial partition computed in %.2fs\n",

573

t1.get_duration());

574

fflush(bliss_verbstr);

575

}

576

577

578

* Allocate space for the labelings

579

580

if(first_path_labeling)

581

free(first_path_labeling);

582

first_path_labeling = (unsigned int*)calloc(N, sizeof(unsigned int));

583

if(best_path_labeling)

584

free(best_path_labeling);

585

best_path_labeling = (unsigned int*)calloc(N, sizeof(unsigned int));

586

587

588

* Are there any non-singleton cells?

589

590

if(p.is_discrete())

591

{

592

update_labeling(best_path_labeling);

593

return;

594

}

595

596

//p.print_signature(stderr); fprintf(stderr, "\n");

597

598

599

* Allocate space for the inverses of the labelings

600

601

if(first_path_labeling_inv)

602

free(first_path_labeling_inv);

603

first_path_labeling_inv = (unsigned int*)calloc(N, sizeof(unsigned int));

604

if(best_path_labeling_inv)

605

free(best_path_labeling_inv);

606

best_path_labeling_inv = (unsigned int*)calloc(N, sizeof(unsigned int));

607

608

609

610

* Allocate space for the automorphisms

611

612

if(first_path_automorphism) free(first_path_automorphism);

613

first_path_automorphism = (unsigned int*)malloc(N * sizeof(unsigned int));

614

if(best_path_automorphism) free(best_path_automorphism);

615

best_path_automorphism = (unsigned int*)malloc(N * sizeof(unsigned int));

616

617

618

619

* Initialize orbit information

620

621

first_path_orbits.init(N);

622

best_path_orbits.init(N);

623

624

625

* Initialize certificate memory

626

627

initialize_certificate();

628

//assert(certificate);

629

assert(certificate_index == 0);

630

631

LevelInfo info;

632

std::vector<LevelInfo> search_stack;

633

std::vector<PathInfo> first_path_info;

634

std::vector<PathInfo> best_path_info;

635

636

search_stack.clear();

637

p.refinement_stack.clean();

638

assert(neighbour_heap.is_empty());

639

640

641

* Initialize long prune

642

643

long_prune_init();

644

645

646

* Build the first level info

647

648

info.split_cell_first = find_next_cell_to_be_splitted(p.first_cell)->first;

649

info.split_element = -1;

650

info.refinement_stack_size = p.refinement_stack.size();

651

info.certificate_index = 0;

652

info.in_first_path = false;

653

info.in_best_path = false;

654

info.long_prune_begin = 0;

655

search_stack.push_back(info);

656

657

658

* Set status and global flags for search related procedures

659

660

in_search = true;

661

refine_compare_certificate = false;

662

stats.nof_leaf_nodes = 0;

663

664

665

#ifdef PRINT_SEARCH_TREE_DOT

666

dotty_output = fopen("debug_stree.dot", "w");

667

fprintf(dotty_output, "digraph stree {\n");

668

fprintf(dotty_output, "\"n\" [label=\"");

669

fprintf(dotty_output, "M"); //p.print(dotty_output);

670

fprintf(dotty_output, "\"];\n");

671

#endif

672

673

p.consistency_check();

674

675

676

* The actual backtracking search

677

678

while(!search_stack.empty())

679

{

680

info = search_stack.back();

681

search_stack.pop_back();

682

683

p.consistency_check();

684

685

686

* Restore partition, certificate index, and split cell

687

688

p.unrefine(p.level, info.refinement_stack_size);

689

assert(info.certificate_index <= certificate_size);

690

certificate_index = info.certificate_index;

691

certificate_current_path.resize(certificate_index);

692

Cell * const cell = p.element_to_cell_map[p.elements[info.split_cell_first]];

693

assert(cell->length > 1);

694

695

p.consistency_check();

696

697

if(p.level > 0 && !info.in_first_path)

698

{

699

if(info.split_element == -1)

700

{

701

info.needs_long_prune = true;

702

}

703

else if(info.needs_long_prune)

704

{

705

info.needs_long_prune = false;

706

/* THIS IS A QUITE HORRIBLE HACK! */

707

unsigned int begin = (info.long_prune_begin>long_prune_begin)?info.long_prune_begin:long_prune_begin;

708

for(unsigned int i = begin; i < long_prune_end; i++)

709

{

710

const std::vector<bool> &fixed = long_prune_get_fixed(i);

711

bool fixes_all = true;

712

for(unsigned int l = 0; l < p.level; l++)

713

{

714

if(fixed[search_stack[l].split_element] == false)

715

{

716

fixes_all = false;

717

break;

718

}

719

}

720

if(!fixes_all)

721

{

722

long_prune_swap(begin, i);

723

begin++;

724

info.long_prune_begin = begin;

725

continue;

726

}

727

const std::vector<bool> &mcrs = long_prune_get_mcrs(i);

728

unsigned int *ep = p.elements + cell->first;

729

for(unsigned int j = cell->length; j > 0; j--, ep++) {

730

if(mcrs[*ep] == false)

731

{

732

info.long_prune_redundant.insert(*ep);

733

}

734

}

735

}

736

}

737

}

738

739

740

* Find the next smallest element in cell

741

742

unsigned int next_split_element = UINT_MAX;

743

unsigned int *next_split_element_pos = 0;

744

unsigned int *ep = p.elements + cell->first;

745

if(info.in_first_path)

746

{

747

/* Find the next larger splitting element that is a mor */

748

for(unsigned int i = cell->length; i > 0; i--, ep++) {

749

if((int)(*ep) > info.split_element &&

750

*ep < next_split_element &&

751

first_path_orbits.is_minimal_representative(*ep)) {

752

next_split_element = *ep;

753

next_split_element_pos = ep;

754

}

755

}

756

}

757

else if(info.in_best_path)

758

{

759

/* Find the next larger splitting element that is a mor */

760

for(unsigned int i = cell->length; i > 0; i--, ep++) {

761

if((int)(*ep) > info.split_element &&

762

*ep < next_split_element &&

763

best_path_orbits.is_minimal_representative(*ep) &&

764

(info.long_prune_redundant.find(*ep) ==

765

info.long_prune_redundant.end())) {

766

next_split_element = *ep;

767

next_split_element_pos = ep;

768

}

769

}

770

}

771

else

772

{

773

/* Find the next larger splitting element */

774

for(unsigned int i = cell->length; i > 0; i--, ep++) {

775

if((int)(*ep) > info.split_element &&

776

*ep < next_split_element &&

777

(info.long_prune_redundant.find(*ep) ==

778

info.long_prune_redundant.end())) {

779

next_split_element = *ep;

780

next_split_element_pos = ep;

781

}

782

}

783

}

784

if(next_split_element == UINT_MAX)

785

{

786

787

* No more splitting elements (unexplored children) in the cell

788

789

/* Update group size if required */

790

if(info.in_first_path == true) {

791

const unsigned int index =

792

first_path_orbits.orbit_size(first_path_info[p.level].splitting_element);

793

stats.group_size.multiply(index);

794

795

* Update all_same_level

796

797

if(index == cell->length && all_same_level == p.level+1)

798

all_same_level = p.level;

799

if (bliss_verbose) {

800

fprintf(stdout,

801

"Level %u: orbits=%u, index=%u/%u, all_same_level=%u\n",

802

p.level,

803

first_path_orbits.nof_orbits(),

804

index, cell->length,

805

all_same_level);

806

fflush(stdout);

807

}

808

}

809

/* Backtrack to the previous level */

810

p.level--;

811

continue;

812

}

813

814

/* Split on smallest */

815

info.split_element = next_split_element;

816

817

818

* Save the current search situation

819

820

search_stack.push_back(info);

821

822

823

* No more in the first path

824

825

info.in_first_path = false;

826

827

* No more in the best path

828

829

info.in_best_path = false;

830

831

p.level++;

832

stats.nof_nodes++;

833

if(p.level > stats.max_level)

834

stats.max_level = p.level;

835

836

p.consistency_check();

837

838

839

* Move the split element to be the last in the cell

840

841

*next_split_element_pos = p.elements[cell->first + cell->length - 1];

842

p.in_pos[*next_split_element_pos] = next_split_element_pos;

843

p.elements[cell->first + cell->length - 1] = next_split_element;

844

p.in_pos[next_split_element] = p.elements+ cell->first + cell->length -1;

845

846

* Split the cell in two:

847

* the last element in the cell (split element) forms a singleton cell

848

849

Cell * const new_cell = p.aux_split_in_two(cell, cell->length - 1);

850

p.element_to_cell_map[p.elements[new_cell->first]] = new_cell;

851

p.consistency_check();

852

853

854

const bool prev_equal_to_first_path = info.equal_to_first_path;

855

const int prev_cmp_to_best_path = info.cmp_to_best_path;

856

857

//assert(!(!info.equal_to_first_path && info.cmp_to_best_path < 0));

858

859

if(!first_path_info.empty())

860

{

861

refine_equal_to_first = info.equal_to_first_path;

862

if(refine_equal_to_first)

863

refine_first_path_subcertificate_end =

864

first_path_info[p.level-1].certificate_index +

865

first_path_info[p.level-1].subcertificate_length;

866

if(canonical)

867

{

868

refine_cmp_to_best = info.cmp_to_best_path;

869

if(refine_cmp_to_best == 0)

870

refine_best_path_subcertificate_end =

871

best_path_info[p.level-1].certificate_index +

872

best_path_info[p.level-1].subcertificate_length;

873

}

874

else

875

refine_cmp_to_best = -1;

876

}

877

878

* Refine the new partition to equitable

879

880

if(cell->length == 1)

881

refine_to_equitable(cell, new_cell);

882

else

883

refine_to_equitable(new_cell);

884

885

p.consistency_check();

886

887

888

#ifdef PRINT_SEARCH_TREE_DOT

889

fprintf(dotty_output, "\"n");

890

for(unsigned int i = 0; i < search_stack.size(); i++) {

891

fprintf(dotty_output, "%u", search_stack[i].split_element);

892

if(i < search_stack.size() - 1) fprintf(dotty_output, ".");

893

}

894

fprintf(dotty_output, "\"");

895

fprintf(dotty_output, " [label=\"");

896

fprintf(dotty_output, "%u",cell->first); /*p.print(dotty_output);*/

897

fprintf(dotty_output, "\"]");

898

if(!first_path_info.empty() && canonical && refine_cmp_to_best > 0) {

899

fprintf(dotty_output, "[color=green]");

900

}

901

fprintf(dotty_output, ";\n");

902

903

fprintf(dotty_output, "\"n");

904

for(unsigned int i = 0; i < search_stack.size() - 1; i++) {

905

fprintf(dotty_output, "%u", search_stack[i].split_element);

906

if(i < search_stack.size() - 2) fprintf(dotty_output, ".");

907

}

908

fprintf(dotty_output, "\" -> \"n");

909

for(unsigned int i = 0; i < search_stack.size(); i++) {

910

fprintf(dotty_output, "%u", search_stack[i].split_element);

911

if(i < search_stack.size() - 1) fprintf(dotty_output, ".");

912

}

913

fprintf(dotty_output, "\" [label=\"%d\"];\n", next_split_element);

914

#endif

915

916

917

if(prev_cmp_to_best_path == 0 && refine_cmp_to_best < 0)

918

fprintf(stderr, "BP- ");

919

if(prev_cmp_to_best_path == 0 && refine_cmp_to_best > 0)

920

fprintf(stderr, "BP+ ");

921

922

923

if(p.is_discrete())

924

{

925

/* Update statistics */

926

stats.nof_leaf_nodes++;

927

928

if(stats.nof_leaf_nodes % 100 == 0) {

929

fprintf(stdout, "Nodes: %lu, Leafs: %lu, Bad: %lu\n",

930

stats.nof_nodes, stats.nof_leaf_nodes,

931

stats.nof_bad_nodes);

932

fflush(stdout);

933

}

934

935

}

936

937

if(!first_path_info.empty())

938

{

939

/* We are no longer on the first path */

940

assert(best_path_info.size() > 0);

941

assert(certificate_current_path.size() >= certificate_index);

942

const unsigned int subcertificate_length =

943

certificate_current_path.size() - certificate_index;

944

if(refine_equal_to_first)

945

{

946

/* Was equal to the first path so far */

947

assert(first_path_info.size() >= p.level);

948

PathInfo &first_pinfo = first_path_info[p.level-1];

949

assert(first_pinfo.certificate_index == certificate_index);

950

if(subcertificate_length != first_pinfo.subcertificate_length)

951

{

952

refine_equal_to_first = false;

953

}

954

else if(first_pinfo.eqref_hash.cmp(eqref_hash) != 0)

955

{

956

refine_equal_to_first = false;

957

}

958

}

959

if(canonical && (refine_cmp_to_best == 0))

960

{

961

/* Was equal to the best path so far */

962

assert(best_path_info.size() >= p.level);

963

PathInfo &best_pinfo = best_path_info[p.level-1];

964

assert(best_pinfo.certificate_index == certificate_index);

965

if(subcertificate_length < best_pinfo.subcertificate_length)

966

{

967

refine_cmp_to_best = -1;

968

//fprintf(stderr, "BSCL- ");

969

}

970

else if(subcertificate_length > best_pinfo.subcertificate_length)

971

{

972

refine_cmp_to_best = 1;

973

//fprintf(stderr, "BSCL+ ");

974

}

975

else if(best_pinfo.eqref_hash.cmp(eqref_hash) > 0)

976

{

977

refine_cmp_to_best = -1;

978

//fprintf(stderr, "BHL- ");

979

}

980

else if(best_pinfo.eqref_hash.cmp(eqref_hash) < 0)

981

{

982

refine_cmp_to_best = 1;

983

//fprintf(stderr, "BHL+ ");

984

}

985

}

986

if(refine_equal_to_first == false &&

987

(!canonical || (refine_cmp_to_best < 0)))

988

{

989

/* Backtrack */

990

#ifdef PRINT_SEARCH_TREE_DOT

991

fprintf(dotty_output, "\"n");

992

for(unsigned int i = 0; i < search_stack.size(); i++) {

993

fprintf(dotty_output, "%u", search_stack[i].split_element);

994

if(i < search_stack.size() - 1) fprintf(dotty_output, ".");

995

}

996

fprintf(dotty_output, "\" [color=red];\n");

997

#endif

998

stats.nof_bad_nodes++;

999

if(search_stack.back().equal_to_first_path == true &&

1000

p.level > all_same_level)

1001

{

1002

assert(all_same_level >= 1);

1003

for(unsigned int i = all_same_level;

1004

i < search_stack.size();

1005

i++)

1006

{

1007

search_stack[i].equal_to_first_path = false;

1008

}

1009

}

1010

while(!search_stack.empty())

1011

{

1012

p.level--;

1013

LevelInfo &info2 = search_stack.back();

1014

if(!(info2.equal_to_first_path == false &&

1015

(!canonical || (info2.cmp_to_best_path < 0))))

1016

break;

1017

search_stack.pop_back();

1018

}

1019

continue;

1020

}

1021

}

1022

1023

#if defined(VERIFY_EQUITABLEDNESS)

1024

/* The new partition should be equitable */

1025

assert(is_equitable());

1026

#endif

1027

1028

info.equal_to_first_path = refine_equal_to_first;

1029

info.cmp_to_best_path = refine_cmp_to_best;

1030

1031

certificate_index = certificate_current_path.size();

1032

1033

search_stack.back().eqref_hash = eqref_hash;

1034

search_stack.back().subcertificate_length =

1035

certificate_index - info.certificate_index;

1036

1037

1038

if(!p.is_discrete())

1039

{

1040

1041

* An internal, non-leaf node

1042

1043

/* Build the next node info */

1044

/* Find the next cell to be splitted */

1045

assert(cell == p.element_to_cell_map[p.elements[info.split_cell_first]]);

1046

Cell * const next_split_cell = find_next_cell_to_be_splitted(cell);

1047

assert(next_split_cell);

1048

/* Copy current info to the search stack */

1049

search_stack.push_back(info);

1050

LevelInfo &new_info = search_stack.back();

1051

new_info.split_cell_first = next_split_cell->first;

1052

new_info.split_element = -1;

1053

new_info.certificate_index = certificate_index;

1054

new_info.refinement_stack_size = p.refinement_stack.size();

1055

new_info.long_prune_redundant.clear();

1056

new_info.long_prune_begin = info.long_prune_begin;

1057

continue;

1058

}

1059

1060

1061

* A leaf node

1062

1063

assert(certificate_index == certificate_size);

1064

1065

if(first_path_info.empty())

1066

{

1067

/* The first path, update first_path and best_path */

1068

//fprintf(stdout, "Level %u: FIRST\n", p.level); fflush(stdout);

1069

stats.nof_canupdates++;

1070

1071

* Update labelings and their inverses

1072

1073

update_labeling_and_its_inverse(first_path_labeling,

1074

first_path_labeling_inv);

1075

update_labeling_and_its_inverse(best_path_labeling,

1076

best_path_labeling_inv);

1077

1078

* Reset automorphism array

1079

1080

reset_permutation(first_path_automorphism);

1081

reset_permutation(best_path_automorphism);

1082

1083

* Reset orbit information

1084

1085

first_path_orbits.reset();

1086

best_path_orbits.reset();

1087

1088

* Reset group size

1089

1090

stats.group_size.assign(1);

1091

1092

* Reset all_same_level

1093

1094

all_same_level = p.level;

1095

1096

* Mark the current path to be the first and best one and save it

1097

1098

const unsigned int base_size = search_stack.size();

1099

assert(p.level == base_size);

1100

best_path_info.clear();

1101

//fprintf(stdout, " New base is: ");

1102

for(unsigned int i = 0; i < base_size; i++) {

1103

search_stack[i].in_first_path = true;

1104

search_stack[i].in_best_path = true;

1105

search_stack[i].equal_to_first_path = true;

1106

search_stack[i].cmp_to_best_path = 0;

1107

PathInfo path_info;

1108

path_info.splitting_element = search_stack[i].split_element;

1109

path_info.certificate_index = search_stack[i].certificate_index;

1110

path_info.eqref_hash = search_stack[i].eqref_hash;

1111

path_info.subcertificate_length = search_stack[i].subcertificate_length;

1112

first_path_info.push_back(path_info);

1113

best_path_info.push_back(path_info);

1114

//fprintf(stdout, "%u ", search_stack[i].split_element);

1115

}

1116

//fprintf(stdout, "\n"); fflush(stdout);

1117

certificate_first_path = certificate_current_path;

1118

certificate_best_path = certificate_current_path;

1119

1120

refine_compare_certificate = true;

1121

1122

* Backtrack to the previous level

1123

1124

p.level--;

1125

continue;

1126

}

1127

1128

DEBUG_ASSERT(first_path_info.size() > 0);

1129

1130

//fprintf(stdout, "Level %u: LEAF %d %d\n", p.level, info.equal_to_first_path, info.cmp_to_best_path); fflush(stdout);

1131

1132

if(info.equal_to_first_path)

1133

{

1134

1135

* An automorphism found: aut[i] = elements[first_path_labeling[i]]

1136

1137

assert(!info.in_first_path);

1138

//fprintf(stdout, "A"); fflush(stdout);

1139

1140

#ifdef PRINT_SEARCH_TREE_DOT

1141

fprintf(dotty_output, "\"n");

1142

for(unsigned int i = 0; i < search_stack.size(); i++) {

1143

fprintf(dotty_output, "%u", search_stack[i].split_element);

1144

if(i < search_stack.size() - 1) fprintf(dotty_output, ".");

1145

}

1146

fprintf(dotty_output, "\" [color=blue];\n");

1147

#endif

1148

1149

#if defined(DEBUG)

1150

/* Verify that the automorphism is correctly built */

1151

for(unsigned int i = 0; i < N; i++)

1152

assert(first_path_automorphism[i] ==

1153

p.elements[first_path_labeling[i]]);

1154

#endif

1155

1156

#if defined(VERIFY_AUTOMORPHISMS)

1157

/* Verify that it really is an automorphism */

1158

assert(is_automorphism(first_path_automorphism));

1159

#endif

1160

1161

long_prune_add_automorphism(first_path_automorphism);

1162

1163

1164

* Update orbit information

1165

1166

update_orbit_information(first_path_orbits, first_path_automorphism);

1167

1168

1169

* Compute backjumping level

1170

1171

unsigned int backjumping_level = 0;

1172

for(unsigned int i = search_stack.size(); i > 0; i--) {

1173

const unsigned int split_element =

1174

search_stack[backjumping_level].split_element;

1175

if(first_path_automorphism[split_element] != split_element)

1176

break;

1177

backjumping_level++;

1178

}

1179

assert(backjumping_level < p.level);

1180

1181

* Go back to backjumping_level

1182

1183

p.level = backjumping_level;

1184

search_stack.resize(p.level + 1);

1185

1186

if(write_automorphisms)

1187

{

1188

print_permutation(stdout, first_path_automorphism);

1189

fprintf(stdout, "\n");

1190

}

1191

stats.nof_generators++;

1192

continue;

1193

}

1194

1195

assert(canonical);

1196

assert(info.cmp_to_best_path >= 0);

1197

if(info.cmp_to_best_path > 0)

1198

{

1199

1200

* A new, better representative found

1201

1202

//fprintf(stdout, "Level %u: NEW BEST\n", p.level); fflush(stdout);

1203

stats.nof_canupdates++;

1204

1205

* Update canonical labeling and its inverse

1206

1207

update_labeling_and_its_inverse(best_path_labeling,

1208

best_path_labeling_inv);

1209

/* Reset best path automorphism */

1210

reset_permutation(best_path_automorphism);

1211

/* Reset best path orbit structure */

1212

best_path_orbits.reset();

1213

1214

* Mark the current path to be the best one and save it

1215

1216

const unsigned int base_size = search_stack.size();

1217

assert(p.level == base_size);

1218

best_path_info.clear();

1219

//fprintf(stdout, " New base is: ");

1220

for(unsigned int i = 0; i < base_size; i++) {

1221

search_stack[i].cmp_to_best_path = 0;

1222

search_stack[i].in_best_path = true;

1223

PathInfo path_info;

1224

path_info.splitting_element = search_stack[i].split_element;

1225

path_info.certificate_index = search_stack[i].certificate_index;

1226

path_info.eqref_hash = search_stack[i].eqref_hash;

1227

path_info.subcertificate_length = search_stack[i].subcertificate_length;

1228

best_path_info.push_back(path_info);

1229

//fprintf(stdout, "%u ", search_stack[i].split_element);

1230

}

1231

certificate_best_path = certificate_current_path;

1232

//fprintf(stdout, "\n"); fflush(stdout);

1233

1234

* Backtrack to the previous level

1235

1236

p.level--;

1237

continue;

1238

}

1239

1240

{

1241

//fprintf(stderr, "BAUT ");

1242

1243

* Equal to the previous best path

1244

1245

#if defined(DEBUG)

1246

/* Verify that the automorphism is correctly built */

1247

for(unsigned int i = 0; i < N; i++)

1248

assert(best_path_automorphism[i] ==

1249

p.elements[best_path_labeling[i]]);

1250

#endif

1251

1252

#if defined(VERIFY_AUTOMORPHISMS)

1253

/* Verify that it really is an automorphism */

1254

assert(is_automorphism(best_path_automorphism));

1255

#endif

1256

1257

unsigned int gca_level_with_first = 0;

1258

for(unsigned int i = search_stack.size(); i > 0; i--) {

1259

if((int)first_path_info[gca_level_with_first].splitting_element !=

1260

search_stack[gca_level_with_first].split_element)

1261

break;

1262

gca_level_with_first++;

1263

}

1264

assert(gca_level_with_first < p.level);

1265

1266

unsigned int gca_level_with_best = 0;

1267

for(unsigned int i = search_stack.size(); i > 0; i--) {

1268

if((int)best_path_info[gca_level_with_best].splitting_element !=

1269

search_stack[gca_level_with_best].split_element)

1270

break;

1271

gca_level_with_best++;

1272

}

1273

assert(gca_level_with_best < p.level);

1274

1275

long_prune_add_automorphism(best_path_automorphism);

1276

1277

1278

* Update orbit information

1279

1280

update_orbit_information(best_path_orbits, best_path_automorphism);

1281

1282

1283

* Update orbit information

1284

1285

const unsigned int nof_old_orbits = first_path_orbits.nof_orbits();

1286

update_orbit_information(first_path_orbits, best_path_automorphism);

1287

if(nof_old_orbits != first_path_orbits.nof_orbits())

1288

{

1289

if(write_automorphisms)

1290

{

1291

print_permutation(stdout, best_path_automorphism);

1292

fprintf(stdout, "\n");

1293

}

1294

stats.nof_generators++;

1295

}

1296

1297

1298

* Compute backjumping level

1299

1300

unsigned int backjumping_level = p.level - 1;

1301

if(!first_path_orbits.is_minimal_representative(search_stack[gca_level_with_first].split_element))

1302

{

1303

backjumping_level = gca_level_with_first;

1304

/*fprintf(stderr, "bj1: %u %u\n", p.level, backjumping_level);*/

1305

}

1306

else

1307

{

1308

assert(!best_path_orbits.is_minimal_representative(search_stack[gca_level_with_best].split_element));

1309

backjumping_level = gca_level_with_best;

1310

/*fprintf(stderr, "bj2: %u %u\n", p.level, backjumping_level);*/

1311

}

1312

/* Backtrack */

1313

search_stack.resize(backjumping_level + 1);

1314

p.level = backjumping_level;

1315

continue;

1316

}

1317

}

1318

1319

#ifdef PRINT_SEARCH_TREE_DOT

1320

fprintf(dotty_output, "}\n");

1321

fclose(dotty_output);

1322

#endif

1323

}

1324

1325

1326

1327

1328

void AbstractGraph::find_automorphisms(Stats &stats)

1329

{

1330

search(false, stats);

1331

1332

if(first_path_labeling)

1333

{

1334

free(first_path_labeling);

1335

first_path_labeling = 0;

1336

}

1337

if(best_path_labeling)

1338

{

1339

free(best_path_labeling);

1340

best_path_labeling = 0;

1341

}

1342

}

1343

1344

1345

const unsigned int *AbstractGraph::canonical_form(Stats &stats)

1346

{

1347

search(true, stats);

1348

1349

return best_path_labeling;

1350

}

1351

1352

1353

1354

1355

/*-------------------------------------------------------------------------

1356

1357

* Routines for undirected graphs

1358

1359

*-------------------------------------------------------------------------*/

1360

1361

Graph::Vertex::Vertex()

1362

{

1363

label = 1;

1364

nof_edges = 0;

1365

}

1366

1367

1368

Graph::Vertex::~Vertex()

1369

{

1370

;

1371

}

1372

1373

1374

void Graph::Vertex::add_edge(const unsigned int other_vertex)

1375

{

1376

edges.push_back(other_vertex);

1377

nof_edges++;

1378

DEBUG_ASSERT(nof_edges == edges.size());

1379

}

1380

1381

1382

void Graph::Vertex::remove_duplicate_edges(bool * const duplicate_array)

1383

{

1384

for(std::vector<unsigned int>::iterator iter = edges.begin();

1385

iter != edges.end(); )

1386

{

1387

const unsigned int dest_vertex = *iter;

1388

if(duplicate_array[dest_vertex] == true)

1389

{

1390

/* A duplicate edge found! */

1391

iter = edges.erase(iter);

1392

nof_edges--;

1393

DEBUG_ASSERT(nof_edges == edges.size());

1394

}

1395

else

1396

{

1397

/* Not seen earlier, mark as seen */

1398

duplicate_array[dest_vertex] = true;

1399

iter++;

1400

}

1401

}

1402

1403

/* Clear duplicate_array */

1404

for(std::vector<unsigned int>::iterator iter = edges.begin();

1405

iter != edges.end();

1406

iter++)

1407

{

1408

duplicate_array[*iter] = false;

1409

}

1410

}

1411

1412

1413

1414

1415

1416

/*-------------------------------------------------------------------------

1417

1418

* Constructor and destructor for undirected graphs

1419

1420

*-------------------------------------------------------------------------*/

1421

1422

1423

Graph::Graph(const unsigned int nof_vertices)

1424

{

1425

vertices.resize(nof_vertices);

1426

sh = sh_flm;

1427

}

1428

1429

1430

Graph::~Graph()

1431

{

1432

;

1433

}

1434

1435

1436

unsigned int Graph::add_vertex(const unsigned int new_label)

1437

{

1438

const unsigned int new_vertex_num = vertices.size();

1439

vertices.resize(new_vertex_num + 1);

1440

vertices.back().label = new_label;

1441

return new_vertex_num;

1442

}

1443

1444

1445

void Graph::add_edge(const unsigned int vertex1, const unsigned int vertex2)

1446

{

1447

//fprintf(stderr, "(%u,%u) ", vertex1, vertex2);

1448

assert(vertex1 < vertices.size());

1449

assert(vertex2 < vertices.size());

1450

vertices[vertex1].add_edge(vertex2);

1451

vertices[vertex2].add_edge(vertex1);

1452

}

1453

1454

1455

void Graph::change_label(const unsigned int vertex,

1456

const unsigned int new_label)

1457

{

1458

assert(vertex < vertices.size());

1459

vertices[vertex].label = new_label;

1460

}

1461

1462

1463

1464

1465

1466

/*-------------------------------------------------------------------------

1467

1468

* Read graph in the DIMACS format

1469

1470

*-------------------------------------------------------------------------*/

1471

1472

Graph *Graph::read_dimacs(FILE *fp)

1473

{

1474

Graph *g = 0;

1475

unsigned int nof_vertices, nof_edges;

1476

unsigned int line_num = 1;

1477

int c;

1478

1479

/* read comments and problem line*/

1480

while(1) {

1481

c = getc(fp);

1482

if(c == 'c') {

1483

while((c = getc(fp)) != '\n') {

1484

if(c == EOF) {

1485

fprintf(stderr, "error in line %u: not in DIMACS format\n",

1486

line_num);

1487

goto error_exit;

1488

}

1489

}

1490

line_num++;

1491

continue;

1492

}

1493

if(c == 'p') {

1494

if(fscanf(fp, " edge %u %u\n", &nof_vertices, &nof_edges) != 2) {

1495

fprintf(stderr, "error in line %u: not in DIMACS format\n",

1496

line_num);

1497

goto error_exit; }

1498

line_num++;

1499

break;

1500

}

1501

fprintf(stderr, "error in line %u: not in DIMACS format\n", line_num);

1502

goto error_exit;

1503

}

1504

1505

if(nof_vertices <= 0) {

1506

fprintf(stderr, "error: no vertices\n");

1507

goto error_exit;

1508

}

1509

#if 0

1510

if(nof_edges <= 0) {

1511

fprintf(stderr, "error: no edges\n");

1512

goto error_exit;

1513

}

1514

#endif

1515

if(bliss_verbose) {

1516

fprintf(bliss_verbstr, "Instance has %d vertices and %d edges\n",

1517

nof_vertices, nof_edges);

1518

fflush(bliss_verbstr);

1519

}

1520

1521

g = new Graph(nof_vertices);

1522

1523

1524

// Read vertex labels

1525

1526

if(bliss_verbose) {

1527

fprintf(bliss_verbstr, "Reading vertex labels...\n");

1528

fflush(bliss_verbstr); }

1529

while(1) {

1530

c = getc(fp);

1531

if(c != 'n') {

1532

ungetc(c, fp);

1533

break;

1534

}

1535

ungetc(c, fp);

1536

unsigned int vertex, label;

1537

if(fscanf(fp, "n %u %u\n", &vertex, &label) != 2) {

1538

fprintf(stderr, "error in line %u: not in DIMACS format\n",

1539

line_num);

1540

goto error_exit;

1541

}

1542

if(vertex > nof_vertices) {

1543

fprintf(stderr, "error in line %u: not in DIMACS format\n",

1544

line_num);

1545

goto error_exit;

1546

}

1547

line_num++;

1548

g->change_label(vertex - 1, label);

1549

}

1550

if(bliss_verbose) {

1551

fprintf(bliss_verbstr, "Done\n");

1552

fflush(bliss_verbstr); }

1553

1554

1555

// Read edges

1556

1557

if(bliss_verbose) {

1558

fprintf(bliss_verbstr, "Reading edges...\n");

1559

fflush(bliss_verbstr); }

1560

for(unsigned i = 0; i < nof_edges; i++) {

1561

unsigned int from, to;

1562

if(fscanf(fp, "e %u %u\n", &from, &to) != 2) {

1563

fprintf(stderr, "error in line %u: not in DIMACS format\n",

1564

line_num);

1565

goto error_exit;

1566

}

1567

if(from > nof_vertices || to > nof_vertices) {

1568

fprintf(stderr, "error in line %u: not in DIMACS format\n",

1569

line_num);

1570

goto error_exit;

1571

}

1572

line_num++;

1573

g->add_edge(from - 1, to - 1);

1574

}

1575

if(bliss_verbose) {

1576

fprintf(bliss_verbstr, "Done\n");

1577

fflush(bliss_verbstr);

1578

}

1579

1580

return g;

1581

1582

error_exit:

1583

if(g)

1584

delete g;

1585

return 0;

1586

1587

}

1588

1589

Graph *Graph::from_igraph(const igraph_t *graph) {

1590

1591

unsigned int nof_vertices= (unsigned int)igraph_vcount(graph);

1592

unsigned int nof_edges= (unsigned int)igraph_ecount(graph);

1593

Graph *g=new Graph(nof_vertices);

1594

// for (unsigned int i=0; i<nof_vertices; i++) {

1595

// g->change_label(i, i);

1596

// }

1597

for (unsigned int i=0; i<nof_edges; i++) {

1598

g->add_edge((unsigned int)IGRAPH_FROM(graph, i),

1599

(unsigned int)IGRAPH_TO(graph, i));

1600

}

1601

return g;

1602

}

1603

1604

void Graph::print_dimacs(FILE *fp)

1605

{

1606

unsigned int nof_edges = 0;

1607

for(unsigned int i = 0; i < get_nof_vertices(); i++)

1608

{

1609

Vertex &v = vertices[i];

1610

for(std::vector<unsigned int>::const_iterator ei = v.edges.begin();

1611

ei != v.edges.end();

1612

ei++)

1613

{

1614

const unsigned int dest_i = *ei;

1615

if(dest_i < i)

1616

continue;

1617

nof_edges++;

1618

}

1619

}

1620

1621

fprintf(fp, "p edge %u %u\n", get_nof_vertices(), nof_edges);

1622

for(unsigned int i = 0; i < get_nof_vertices(); i++)

1623

{

1624

Vertex &v = vertices[i];

1625

if(v.label != 1)

1626

{

1627

fprintf(fp, "n %u %u\n", i+1, v.label);

1628

}

1629

}

1630

for(unsigned int i = 0; i < get_nof_vertices(); i++)

1631

{

1632

Vertex &v = vertices[i];

1633

for(std::vector<unsigned int>::const_iterator ei = v.edges.begin();

1634

ei != v.edges.end();

1635

ei++)

1636

{

1637

const unsigned int dest_i = *ei;

1638

if(dest_i < i)

1639

continue;

1640

fprintf(fp, "e %u %u\n", i+1, dest_i+1);

1641

}

1642

}

1643

}

1644

1645

1646

1647

1648

Graph *Graph::permute(const unsigned int *perm)

1649

{

1650

Graph *g = new Graph(get_nof_vertices());

1651

for(unsigned int i = 0; i < get_nof_vertices(); i++)

1652

{

1653

Vertex &v = vertices[i];

1654

Vertex &permuted_v = g->vertices[perm[i]];

1655

permuted_v.label = v.label;

1656

for(std::vector<unsigned int>::const_iterator ei = v.edges.begin();

1657

ei != v.edges.end();

1658

ei++)

1659

{

1660

const unsigned int dest_v = *ei;

1661

permuted_v.add_edge(perm[dest_v]);

1662

}

1663

std::sort(permuted_v.edges.begin(), permuted_v.edges.end());

1664

}

1665

return g;

1666

}

1667

1668

1669

1670

1671

1672

/*-------------------------------------------------------------------------

1673

1674

* Print graph in graphviz format

1675

1676

*-------------------------------------------------------------------------*/

1677

1678

1679

void Graph::to_dot(const char *file_name)

1680

{

1681

FILE *fp = fopen(file_name, "w");

1682

if(fp)

1683

to_dot(fp);

1684

fclose(fp);

1685

}

1686

1687

1688

void Graph::to_dot(FILE *fp)

1689

{

1690

remove_duplicate_edges();

1691

1692

fprintf(fp, "graph g {\n");

1693

1694

unsigned int vnum = 0;

1695

for(std::vector<Vertex>::iterator vi = vertices.begin();

1696

vi != vertices.end();

1697

vi++, vnum++)

1698

{

1699

Vertex &v = *vi;

1700

fprintf(fp, "v%u [label=\"%u:%u\"];\n", vnum, vnum, v.label);

1701

for(std::vector<unsigned int>::const_iterator ei = v.edges.begin();

1702

ei != v.edges.end();

1703

ei++)

1704

{

1705

const unsigned int vnum2 = *ei;

1706

if(vnum2 > vnum)

1707

fprintf(fp, "v%u -- v%u\n", vnum, vnum2);

1708

}

1709

}

1710

1711

fprintf(fp, "}\n");

1712

}

1713

1714

1715

1716

1717

1718

void Graph::remove_duplicate_edges()

1719

{

1720

bool *duplicate_array = (bool*)calloc(vertices.size(), sizeof(bool));

1721

1722

for(std::vector<Vertex>::iterator vi = vertices.begin();

1723

vi != vertices.end();

1724

vi++)

1725

{

1726

#ifdef EXPENSIVE_CONSISTENCY_CHECKS

1727

for(unsigned int i = 0; i < vertices.size(); i++)

1728

assert(duplicate_array[i] == false);

1729

#endif

1730

Vertex &v = *vi;

1731

v.remove_duplicate_edges(duplicate_array);

1732

}

1733

1734

free(duplicate_array);

1735

}

1736

1737

1738

1739

1740

1741

/*-------------------------------------------------------------------------

1742

1743

* Partition independent invariants

1744

1745

*-------------------------------------------------------------------------*/

1746

1747

1748

unsigned int Graph::label_invariant(Graph *g, unsigned int v)

1749

{

1750

DEBUG_ASSERT(v < g->vertices.size());

1751

return g->vertices[v].label;

1752

}

1753

1754

1755

unsigned int Graph::degree_invariant(Graph *g, unsigned int v)

1756

{

1757

DEBUG_ASSERT(v < g->vertices.size());

1758

DEBUG_ASSERT(g->vertices[v].edges.size() ==

1759

g->vertices[v].nof_edges);

1760

return g->vertices[v].nof_edges;

1761

}

1762

1763

1764

1765

1766

1767

1768

1769

/*-------------------------------------------------------------------------

1770

1771

* Refine the partition p according to a partition independent invariant

1772

1773

*-------------------------------------------------------------------------*/

1774

1775

bool Graph::refine_according_to_invariant(unsigned int (*inv)(Graph * const g, unsigned int v))

1776

{

1777

bool refined = false;

1778

1779

for(Cell *cell = p.first_cell; cell; )

1780

{

1781

assert(cell->max_ival == 0);

1782

assert(cell->max_ival_count == 0);

1783

1784

Cell * const next_cell = cell->next;

1785

1786

if(cell->length == 1)

1787

{

1788

cell = next_cell;

1789

continue;

1790

}

1791

1792

const unsigned int *ep = p.elements + cell->first;

1793

for(unsigned int i = cell->length; i > 0; i--, ep++)

1794

{

1795

unsigned int ival = inv(this, *ep);

1796

p.invariant_values[*ep] = ival;

1797

if(ival > cell->max_ival) {

1798

cell->max_ival = ival;

1799

cell->max_ival_count = 1;

1800

}

1801

else if(ival == cell->max_ival) {

1802

cell->max_ival_count++;

1803

}

1804

}

1805

Cell * const last_new_cell = p.zplit_cell(cell, true);

1806

refined = (last_new_cell != cell);

1807

cell = next_cell;

1808

}

1809

1810

return refined;

1811

}

1812

1813

1814

1815

1816

1817

/*-------------------------------------------------------------------------

1818

1819

* Split the neighbourhood of a cell according to the equitable invariant

1820

1821

*-------------------------------------------------------------------------*/

1822

1823

1824

void Graph::split_neighbourhood_of_cell(Cell * const cell)

1825

{

1826

DEBUG_ASSERT(neighbour_heap.is_empty());

1827

DEBUG_ASSERT(cell->length > 1);

1828

1829

eqref_hash.update(cell->first);

1830

eqref_hash.update(cell->length);

1831

1832

unsigned int *ep = p.elements + cell->first;

1833

for(unsigned int i = cell->length; i > 0; i--)

1834

{

1835

const Vertex &v = vertices[*ep];

1836

ep++;

1837

1838

std::vector<unsigned int>::const_iterator ei = v.edges.begin();

1839

for(unsigned int j = v.nof_edges; j > 0; j--)

1840

{

1841

const unsigned int dest_vertex = *ei++;

1842

Cell * const neighbour_cell = p.element_to_cell_map[dest_vertex];

1843

if(neighbour_cell->length == 1)

1844

continue;

1845

const unsigned int ival = p.invariant_values[dest_vertex] + 1;

1846

p.invariant_values[dest_vertex] = ival;

1847

if(ival > neighbour_cell->max_ival) {

1848

neighbour_cell->max_ival = ival;

1849

neighbour_cell->max_ival_count = 1;

1850

}

1851

else if(ival == neighbour_cell->max_ival) {

1852

neighbour_cell->max_ival_count++;

1853

}

1854

if(!neighbour_cell->in_neighbour_heap) {

1855

neighbour_cell->in_neighbour_heap = true;

1856

neighbour_heap.insert(neighbour_cell->first);

1857

}

1858

}

1859

}

1860

1861

while(!neighbour_heap.is_empty())

1862

{

1863

const unsigned int start = neighbour_heap.remove();

1864

Cell * const neighbour_cell = p.element_to_cell_map[p.elements[start]];

1865

DEBUG_ASSERT(neighbour_cell->first == start);

1866

DEBUG_ASSERT(neighbour_cell->in_neighbour_heap);

1867

neighbour_cell->in_neighbour_heap = false;

1868

1869

DEBUG_ASSERT(neighbour_cell->length > 1);

1870

DEBUG_ASSERT(neighbour_cell->max_ival >= 1);

1871

DEBUG_ASSERT(neighbour_cell->max_ival_count >= 1);

1872

1873

eqref_hash.update(neighbour_cell->first);

1874

eqref_hash.update(neighbour_cell->length);

1875

eqref_hash.update(neighbour_cell->max_ival);

1876

eqref_hash.update(neighbour_cell->max_ival_count);

1877

1878

Cell * const last_new_cell = p.zplit_cell(neighbour_cell, true);

1879

/* Update hash */

1880

const Cell *c = neighbour_cell;

1881

while(1)

1882

{

1883

eqref_hash.update(c->first);

1884

eqref_hash.update(c->length);

1885

if(c == last_new_cell)

1886

break;

1887

c = c->next;

1888

}

1889

}

1890

}

1891

1892

1893

bool Graph::split_neighbourhood_of_unit_cell(Cell * const unit_cell)

1894

{

1895

DEBUG_ASSERT(neighbour_heap.is_empty());

1896

1897

DEBUG_ASSERT(unit_cell->length == 1);

1898

DEBUG_ASSERT(p.element_to_cell_map[p.elements[unit_cell->first]] == unit_cell);

1899

DEBUG_ASSERT(p.in_pos[p.elements[unit_cell->first]] ==

1900

p.elements + unit_cell->first);

1901

1902

eqref_hash.update(0x87654321);

1903

eqref_hash.update(unit_cell->first);

1904

eqref_hash.update(1);

1905

1906

const Vertex &v = vertices[p.elements[unit_cell->first]];

1907

1908

std::vector<unsigned int>::const_iterator ei = v.edges.begin();

1909

for(unsigned int j = v.nof_edges; j > 0; j--)

1910

{

1911

const unsigned int dest_vertex = *ei++;

1912

Cell * const neighbour_cell = p.element_to_cell_map[dest_vertex];

1913

DEBUG_ASSERT(*p.in_pos[dest_vertex] == dest_vertex);

1914

1915

if(neighbour_cell->length == 1) {

1916

DEBUG_ASSERT(!neighbour_cell->in_neighbour_heap);

1917

if(in_search) {

1918

neighbour_cell->in_neighbour_heap = true;

1919

neighbour_heap.insert(neighbour_cell->first);

1920

}

1921

continue;

1922

}

1923

if(!neighbour_cell->in_neighbour_heap) {

1924

neighbour_cell->in_neighbour_heap = true;

1925

neighbour_heap.insert(neighbour_cell->first);

1926

}

1927

neighbour_cell->max_ival_count++;

1928

DEBUG_ASSERT(neighbour_cell->max_ival_count <= neighbour_cell->length);

1929

1930

unsigned int * const swap_position =

1931

p.elements + neighbour_cell->first + neighbour_cell->length -

1932

neighbour_cell->max_ival_count;

1933

DEBUG_ASSERT(p.in_pos[dest_vertex] <= swap_position);

1934

*p.in_pos[dest_vertex] = *swap_position;

1935

p.in_pos[*swap_position] = p.in_pos[dest_vertex];

1936

*swap_position = dest_vertex;

1937

p.in_pos[dest_vertex] = swap_position;

1938

}

1939

1940

while(!neighbour_heap.is_empty())

1941

{

1942

const unsigned int start = neighbour_heap.remove();

1943

Cell *neighbour_cell = p.element_to_cell_map[p.elements[start]];

1944

DEBUG_ASSERT(neighbour_cell->in_neighbour_heap);

1945

neighbour_cell->in_neighbour_heap = false;

1946

1947

#ifdef DEBUG

1948

assert(neighbour_cell->first == start);

1949

if(neighbour_cell->length == 1) {

1950

assert(neighbour_cell->max_ival_count == 0);

1951

} else {

1952

assert(neighbour_cell->max_ival_count > 0);

1953

assert(neighbour_cell->max_ival_count <= neighbour_cell->length);

1954

}

1955

#endif

1956

1957

eqref_hash.update(neighbour_cell->first);

1958

eqref_hash.update(neighbour_cell->length);

1959

eqref_hash.update(neighbour_cell->max_ival_count);

1960

1961

if(neighbour_cell->length > 1 &&

1962

neighbour_cell->max_ival_count != neighbour_cell->length) {

1963

1964

p.consistency_check();

1965

1966

Cell * const new_cell = p.aux_split_in_two(neighbour_cell, neighbour_cell->length - neighbour_cell->max_ival_count);

1967

unsigned int *ep = p.elements + new_cell->first;

1968

unsigned int * const lp = p.elements+new_cell->first+new_cell->length;

1969

while(ep < lp) {

1970

DEBUG_ASSERT(p.in_pos[*ep] == ep);

1971

p.element_to_cell_map[*ep] = new_cell;

1972

ep++;

1973

}

1974

neighbour_cell->max_ival_count = 0;

1975

1976

p.consistency_check();

1977

1978

/* update hash */

1979

eqref_hash.update(neighbour_cell->first);

1980

eqref_hash.update(neighbour_cell->length);

1981

eqref_hash.update(0);

1982

eqref_hash.update(new_cell->first);

1983

eqref_hash.update(new_cell->length);

1984

eqref_hash.update(1);

1985

1986

/* Add cells in splitting_queue */

1987

DEBUG_ASSERT(!new_cell->in_splitting_queue);

1988

if(neighbour_cell->in_splitting_queue) {

1989

/* Both cells must be included in splitting_queue in order

1990

to have refinement to equitable partition */

1991

p.add_in_splitting_queue(new_cell);

1992

} else {

1993

Cell *min_cell, *max_cell;

1994

if(neighbour_cell->length <= new_cell->length) {

1995

min_cell = neighbour_cell;

1996

max_cell = new_cell;

1997

} else {

1998

min_cell = new_cell;

1999

max_cell = neighbour_cell;

2000

}

2001

/* Put the smaller cell in splitting_queue */

2002

p.add_in_splitting_queue(min_cell);

2003

if(max_cell->length == 1) {

2004

/* Put the "larger" cell also in splitting_queue */

2005

p.add_in_splitting_queue(max_cell);

2006

}

2007

}

2008

/* Update pointer for certificate generation */

2009

neighbour_cell = new_cell;

2010

}

2011

else

2012

neighbour_cell->max_ival_count = 0;

2013

2014

2015

* Build certificate if required

2016

2017

if(in_search)

2018

{

2019

for(unsigned int i = neighbour_cell->first,

2020

j = neighbour_cell->length,

2021

c_index = certificate_current_path.size();

2022

j > 0;

2023

j--, i++, c_index += 2)

2024

{

2025

if(refine_compare_certificate)

2026

{

2027

if(refine_equal_to_first)

2028

{

2029

if(c_index >= refine_first_path_subcertificate_end)

2030

refine_equal_to_first = false;

2031

else if(certificate_first_path[c_index] !=

2032

unit_cell->first)

2033

refine_equal_to_first = false;

2034

else if(certificate_first_path[c_index+1] != i)

2035

refine_equal_to_first = false;

2036

}

2037

if(refine_cmp_to_best == 0)

2038

{

2039

if(c_index >= refine_best_path_subcertificate_end)

2040

{

2041

refine_cmp_to_best = 1;

2042

}

2043

else if(unit_cell->first>certificate_best_path[c_index])

2044

{

2045

refine_cmp_to_best = 1;

2046

}

2047

else if(unit_cell->first<certificate_best_path[c_index])

2048

{

2049

refine_cmp_to_best = -1;

2050

}

2051

else if(i > certificate_best_path[c_index+1])

2052

{

2053

refine_cmp_to_best = 1;

2054

}

2055

else if(i < certificate_best_path[c_index+1])

2056

{

2057

refine_cmp_to_best = -1;

2058

}

2059

}

2060

if((refine_equal_to_first == false) &&

2061

(refine_cmp_to_best < 0))

2062

goto worse_exit;

2063

}

2064

certificate_current_path.push_back(unit_cell->first);

2065

certificate_current_path.push_back(i);

2066

}

2067

} /* if(in_search) */

2068

} /* while(!neighbour_heap.is_empty()) */

2069

2070

return false;

2071

2072

worse_exit:

2073

while(!neighbour_heap.is_empty())

2074

{

2075

const unsigned int start = neighbour_heap.remove();

2076

Cell * const neighbour_cell = p.element_to_cell_map[p.elements[start]];

2077

DEBUG_ASSERT(neighbour_cell->in_neighbour_heap);

2078

neighbour_cell->in_neighbour_heap = false;

2079

neighbour_cell->max_ival_count = 0;

2080

}

2081

return true;

2082

}

2083

2084

2085

2086

2087

2088

/*-------------------------------------------------------------------------

2089

2090

* Check whether the current partition p is equitable

2091

* Slow: use only for debugging purposes

2092

* Side effect: resets max_ival and max_ival_count fields in cells

2093

2094

*-------------------------------------------------------------------------*/

2095

2096

bool Graph::is_equitable()

2097

{

2098

bool result = true;

2099

2100

2101

* Max ival and max_ival_count are used for counting purposes,

2102

* they should be reset...

2103

2104

for(Cell *cell = p.first_cell; cell; cell = cell->next)

2105

{

2106

assert(cell->prev_next_ptr && *(cell->prev_next_ptr) == cell);

2107

assert(cell->max_ival == 0);

2108

assert(cell->max_ival_count == 0);

2109

}

2110

2111

2112

for(Cell *cell = p.first_cell; cell; cell = cell->next)

2113

{

2114

if(cell->length == 1)

2115

continue;

2116

2117

unsigned int *ep = p.elements + cell->first;

2118

Vertex &first_vertex = vertices[*ep++];

2119

2120

/* Count edges of the first vertex for cells in max_ival */

2121

std::vector<unsigned int>::const_iterator ei = first_vertex.edges.begin();

2122

for(unsigned int j = first_vertex.nof_edges; j > 0; j--)

2123

{

2124

p.element_to_cell_map[*ei++]->max_ival++;

2125

}

2126

2127

/* Count and compare edges of the other vertices */

2128

for(unsigned int i = cell->length; i > 1; i--)

2129

{

2130

Vertex &vertex = vertices[*ep++];

2131

std::vector<unsigned int>::const_iterator ei = vertex.edges.begin();

2132

for(unsigned int j = vertex.nof_edges; j > 0; j--)

2133

{

2134

p.element_to_cell_map[*ei++]->max_ival_count++;

2135

}

2136

for(Cell *cell2 = p.first_cell; cell2; cell2 = cell2->next)

2137

{

2138

if(cell2->max_ival != cell2->max_ival_count)

2139

{

2140

result = false;

2141

goto done;

2142

}

2143

cell2->max_ival_count = 0;

2144

}

2145

}

2146

/* Reset max_ival */

2147

for(Cell *cell2 = p.first_cell; cell2; cell2 = cell2->next)

2148

{

2149

cell2->max_ival = 0;

2150

assert(cell2->max_ival_count == 0);

2151

}

2152

}

2153

2154

done:

2155

2156

for(Cell *cell = p.first_cell; cell; cell = cell->next)

2157

{

2158

cell->max_ival = 0;

2159

cell->max_ival_count = 0;

2160

}

2161

2162

return result;

2163

}

2164

2165

2166

2167

2168

2169

/*-------------------------------------------------------------------------

2170

2171

* Build the initial equitable partition

2172

2173

*-------------------------------------------------------------------------*/

2174

2175

void Graph::make_initial_equitable_partition()

2176

{

2177

refine_according_to_invariant(&label_invariant);

2178

p.clear_splitting_queue();

2179

//p.print_signature(stderr); fprintf(stderr, "\n");

2180

2181

refine_according_to_invariant(&degree_invariant);

2182

p.clear_splitting_queue();

2183

//p.print_signature(stderr); fprintf(stderr, "\n");

2184

2185

/* To do: add loop invariant */

2186

2187

refine_to_equitable();

2188

p.refinement_stack.clean();

2189

//p.print_signature(stderr); fprintf(stderr, "\n");

2190

}

2191

2192

2193

2194

2195

2196

/*-------------------------------------------------------------------------

2197

2198

* Find the next cell to be splitted

2199

2200

*-------------------------------------------------------------------------*/

2201

2202

Cell *Graph::find_next_cell_to_be_splitted(Cell *cell)

2203

{

2204

assert(!p.is_discrete());

2205

switch(sh) {

2206

case sh_f:

2207

return sh_first(cell);

2208

case sh_fs:

2209

return sh_first_smallest(cell);

2210

case sh_fl:

2211

return sh_first_largest(cell);

2212

case sh_fm:

2213

return sh_first_max_neighbours(cell);

2214

case sh_fsm:

2215

return sh_first_smallest_max_neighbours(cell);

2216

case sh_flm:

2217

return sh_first_largest_max_neighbours(cell);

2218

default:

2219

assert(false && "Unknown splitting heuristics");

2220

}

2221

}

2222

2223

/* First nonsingleton cell */

2224

Cell *Graph::sh_first(Cell *cell)

2225

{

2226

return p.first_nonsingleton_cell;

2227

}

2228

2229

/* First smallest nonsingleton cell. */

2230

Cell *Graph::sh_first_smallest(Cell *cell)

2231

{

2232

Cell *best_cell = 0;

2233

unsigned int best_size = UINT_MAX;

2234

for(cell = p.first_nonsingleton_cell; cell; cell = cell->next_nonsingleton)

2235

{

2236

assert(cell->length > 1);

2237

if(cell->length < best_size)

2238

{

2239

best_size = cell->length;

2240

best_cell = cell;

2241

}

2242

}

2243

assert(best_cell);

2244

return best_cell;

2245

}

2246

2247

/* First largest nonsingleton cell. */

2248

Cell *Graph::sh_first_largest(Cell *cell)

2249

{

2250

Cell *best_cell = 0;

2251

unsigned int best_size = 0;

2252

for(cell = p.first_nonsingleton_cell; cell; cell = cell->next_nonsingleton)

2253

{

2254

assert(cell->length > 1);

2255

if(cell->length > best_size)

2256

{

2257

best_size = cell->length;

2258

best_cell = cell;

2259

}

2260

}

2261

assert(best_cell);

2262

return best_cell;

2263

}

2264

2265

/* First nonsingleton cell with max number of neighbouring

2266

* nonsingleton cells.

2267

* Assumes that the partition p is equitable.

2268

* Messes up in_neighbour_heap and max_ival fields of cells

2269

* (assumes they are false/0).

2270

2271

Cell *Graph::sh_first_max_neighbours(Cell *cell)

2272

{

2273

Cell *best_cell = 0;

2274

int best_value = -1;

2275

for(cell = p.first_nonsingleton_cell; cell; cell = cell->next_nonsingleton)

2276

{

2277

assert(cell->length > 1);

2278

2279

const Vertex &v = vertices[p.elements[cell->first]];

2280

std::vector<unsigned int>::const_iterator ei = v.edges.begin();

2281

std::list<Cell*> neighbour_cells_visited;

2282

for(unsigned int j = v.nof_edges; j > 0; j--)

2283

{

2284

const unsigned int dest_vertex = *ei++;

2285

Cell * const neighbour_cell = p.element_to_cell_map[dest_vertex];

2286

if(neighbour_cell->length == 1)

2287

continue;

2288

neighbour_cell->max_ival++;

2289

if(neighbour_cell->in_neighbour_heap)

2290

continue;

2291

neighbour_cell->in_neighbour_heap = true;

2292

neighbour_cells_visited.push_back(neighbour_cell);

2293

}

2294

int value = 0;

2295

while(!neighbour_cells_visited.empty())

2296

{

2297

Cell * const neighbour_cell = neighbour_cells_visited.front();

2298

neighbour_cells_visited.pop_front();

2299

assert(neighbour_cell->in_neighbour_heap);

2300

neighbour_cell->in_neighbour_heap = false;

2301

if(neighbour_cell->max_ival != neighbour_cell->length)

2302

value++;

2303

neighbour_cell->max_ival = 0;

2304

}

2305

if(value > best_value)

2306

{

2307

best_value = value;

2308

best_cell = cell;

2309

}

2310

}

2311

assert(best_cell);

2312

return best_cell;

2313

}

2314

/* First smallest nonsingleton cell with max number of neighbouring

2315

* nonsingleton cells.

2316

* Assumes that the partition p is equitable.

2317

* Messes up in_neighbour_heap and max_ival fields of cells

2318

* (assumes they are false).

2319

2320

Cell *Graph::sh_first_smallest_max_neighbours(Cell *cell)

2321

{

2322

Cell *best_cell = 0;

2323

int best_value = -1;

2324

int best_size = INT_MAX;

2325

for(cell = p.first_nonsingleton_cell; cell; cell = cell->next_nonsingleton)

2326

{

2327

assert(cell->length > 1);

2328

2329

const Vertex &v = vertices[p.elements[cell->first]];

2330

std::vector<unsigned int>::const_iterator ei = v.edges.begin();

2331

std::list<Cell*> neighbour_cells_visited;

2332

for(unsigned int j = v.nof_edges; j > 0; j--)

2333

{

2334

const unsigned int dest_vertex = *ei++;

2335

Cell * const neighbour_cell = p.element_to_cell_map[dest_vertex];

2336

if(neighbour_cell->length == 1)

2337

continue;

2338

neighbour_cell->max_ival++;

2339

if(neighbour_cell->in_neighbour_heap)

2340

continue;

2341

neighbour_cell->in_neighbour_heap = true;

2342

neighbour_cells_visited.push_back(neighbour_cell);

2343

}

2344

int value = 0;

2345

while(!neighbour_cells_visited.empty())

2346

{

2347

Cell * const neighbour_cell = neighbour_cells_visited.front();

2348

neighbour_cells_visited.pop_front();

2349

assert(neighbour_cell->in_neighbour_heap);

2350

neighbour_cell->in_neighbour_heap = false;

2351

if(neighbour_cell->max_ival != neighbour_cell->length)

2352

value++;

2353

neighbour_cell->max_ival = 0;

2354

}

2355

if((value > best_value) ||

2356

(value == best_value && (int)cell->length < best_size))

2357

{

2358

best_value = value;

2359

best_size = cell->length;

2360

best_cell = cell;

2361

}

2362

}

2363

assert(best_cell);

2364

return best_cell;

2365

}

2366

/* First largest nonsingleton cell with max number of neighbouring

2367

* nonsingleton cells.

2368

* Assumes that the partition p is equitable.

2369

* Messes up in_neighbour_heap and max_ival fields of cells

2370

* (assumes they are false/0).

2371

2372

Cell *Graph::sh_first_largest_max_neighbours(Cell *cell)

2373

{

2374

Cell *best_cell = 0;

2375

int best_value = -1;

2376

int best_size = -1;

2377

for(cell = p.first_nonsingleton_cell; cell; cell = cell->next_nonsingleton)

2378

{

2379

assert(cell->length > 1);

2380

2381

const Vertex &v = vertices[p.elements[cell->first]];

2382

std::vector<unsigned int>::const_iterator ei = v.edges.begin();

2383

std::list<Cell*> neighbour_cells_visited;

2384

for(unsigned int j = v.nof_edges; j > 0; j--)

2385

{

2386

const unsigned int dest_vertex = *ei++;

2387

Cell * const neighbour_cell = p.element_to_cell_map[dest_vertex];

2388

if(neighbour_cell->length == 1)

2389

continue;

2390

neighbour_cell->max_ival++;

2391

if(neighbour_cell->in_neighbour_heap)

2392

continue;

2393

neighbour_cell->in_neighbour_heap = true;

2394

neighbour_cells_visited.push_back(neighbour_cell);

2395

}

2396

int value = 0;

2397

while(!neighbour_cells_visited.empty())

2398

{

2399

Cell * const neighbour_cell = neighbour_cells_visited.front();

2400

neighbour_cells_visited.pop_front();

2401

assert(neighbour_cell->in_neighbour_heap);

2402

neighbour_cell->in_neighbour_heap = false;

2403

if(neighbour_cell->max_ival != neighbour_cell->length)

2404

value++;

2405

neighbour_cell->max_ival = 0;

2406

}

2407

if((value > best_value) ||

2408

(value == best_value && (int)cell->length > best_size))

2409

{

2410

best_value = value;

2411

best_size = cell->length;

2412

best_cell = cell;

2413

}

2414

}

2415

assert(best_cell);

2416

return best_cell;

2417

}

2418

2419

2420

2421

2422

2423

/*-------------------------------------------------------------------------

2424

2425

* Initialize the certificate size and memory

2426

2427

*-------------------------------------------------------------------------*/

2428

2429

void Graph::initialize_certificate()

2430

{

2431

certificate_size = 0;

2432

for(Cell *cell = p.first_cell; cell; cell = cell->next)

2433

{

2434

if(cell->length > 1) {

2435

certificate_size +=

2436

vertices[p.elements[cell->first]].nof_edges * 2 * cell->length;

2437

}

2438

}

2439

//if(certificate)

2440

// free(certificate);

2441

//certificate = (unsigned int*)malloc(certificate_size * sizeof(unsigned int));

2442

certificate_index = 0;

2443

2444

certificate_current_path.clear();

2445

certificate_first_path.clear();

2446

certificate_best_path.clear();

2447

}

2448

2449

2450

2451

2452

2453

/*-------------------------------------------------------------------------

2454

2455

* Check whether perm is an automorphism

2456

2457

*-------------------------------------------------------------------------*/

2458

2459

bool Graph::is_automorphism(unsigned int * const perm)

2460

{

2461

std::set<unsigned int, std::less<unsigned int> > edges1;

2462

std::set<unsigned int, std::less<unsigned int> > edges2;

2463

2464

bool result = true;

2465

2466

for(unsigned int i = 0; i < vertices.size(); i++)

2467

{

2468

Vertex &v1 = vertices[i];

2469

edges1.clear();

2470

for(std::vector<unsigned int>::iterator ei = v1.edges.begin();

2471

ei != v1.edges.end();

2472

ei++)

2473

edges1.insert(perm[*ei]);

2474

2475

Vertex &v2 = vertices[perm[i]];

2476

edges2.clear();

2477

for(std::vector<unsigned int>::iterator ei = v2.edges.begin();

2478

ei != v2.edges.end();

2479

ei++)

2480

edges2.insert(*ei);

2481

2482

if(!(edges1 == edges2))

2483

{

2484

result = false;

2485

goto done;

2486

}

2487

}

2488

2489

done:

2490

2491

return result;

2492

}

2493

2494

}

Older »