fn glauber_color(graph: &AdjacencyList, unif_colors: u32, mask: &[bool]) -> (u32, Vec<u32>) {
    let (ncolors, mut colors) = greedy_color(&graph);
    assert!(ncolors < unif_colors);

    // https://www.math.cmu.edu/~af1p/Texfiles/colorbdd.pdf
    // https://www.math.cmu.edu/~af1p/Teaching/MCC17/Papers/colorJ.pdf
    // run glauber markov chain on a coloring
    // chain sampling can be parallel with some simple conflict detection

    let max_degree = graph.max_degree(mask);
    let avg_degree = graph.avg_degree(mask);
    let subgraph_sz = mask.iter().copied().filter(|x| *x).count();

    // by expectation version of birthday problem, don't want more than this
    // much parallelism
    let max_parallel = ((subgraph_sz as f64).sqrt() / avg_degree as f64) as usize;
    let nthreads = (rayon::current_num_threads() as usize).min(max_parallel.max(16));

    println!(
        "chose parallelism {} (recommendation {}) for {}-coloring",
        nthreads, max_parallel, unif_colors
    );

    use std::sync::atomic::{AtomicU32, Ordering};

    let samples = max_degree * subgraph_sz;

    let colors = colors
        .into_iter()
        .map(|x| AtomicU32::new(x))
        .collect::<Vec<_>>();

    let conflicts = (0..nthreads)
        .into_par_iter()
        .map(|i| {
            let samples = (samples / nthreads).max(1);
            let mut rng = StdRng::seed_from_u64((1234 + i) as u64); // no idea how to do this better in rust
            let mut adjacent_color: Vec<bool> = vec![false; unif_colors as usize];
            let mut neighbor_colors: Vec<u32> = Vec::with_capacity(max_degree as usize);
            let mut conflicts = 0;

            for _ in 0..samples {
                // loop invariant is that none of adjacent_color elements are true

                let reservation = (unif_colors + 1) as u32;

                loop {
                    let v: u32 = rng.gen_range(0, graph.nvertices() as u32);
                    if !mask[v as usize] {
                        continue;
                    }

                    // "claim" this vertex
                    let prev = colors[v as usize].swap(reservation, Ordering::Relaxed);
                    if prev == reservation {
                        conflicts += 1;
                        continue;
                    }

                    let mut nadjacent_colors = 0;
                    let mut fail = v;
                    for w in graph.masked_neighbors(v, mask) {
                        let c = colors[w as usize].swap(reservation, Ordering::Relaxed);
                        if c == reservation {
                            fail = w;
                            break;
                        }
                        if !adjacent_color[c as usize] {
                            nadjacent_colors += 1;
                            adjacent_color[c as usize] = true;
                        }
                        neighbor_colors.push(c);
                    }
                    if fail != v {
                        // unlock
                        colors[v as usize].store(prev, Ordering::Relaxed);
                        graph
                            .masked_neighbors(v, mask)
                            .take_while(|&w| w != fail)
                            .zip(neighbor_colors.iter().copied())
                            .for_each(|(w, c)| colors[w as usize].store(c, Ordering::Relaxed));
                        conflicts += 1;
                        for i in adjacent_color.iter_mut() {
                            *i = false;
                        }
                        neighbor_colors.clear();
                        continue;
                    }

                    let chosen = adjacent_color
                        .iter()
                        .copied()
                        .enumerate()
                        .filter(|(_, x)| !x)
                        .map(|(i, _)| i)
                        .nth(rng.gen_range(0, (unif_colors as usize) - nadjacent_colors))
                        .expect("nth");

                    graph
                        .masked_neighbors(v, mask)
                        .zip(neighbor_colors.iter().copied())
                        .for_each(|(w, c)| colors[w as usize].store(c, Ordering::Relaxed));

                    colors[v as usize].store(chosen as u32, Ordering::Relaxed);

                    for i in adjacent_color.iter_mut() {
                        *i = false;
                    }
                    neighbor_colors.clear();
                    break;
                }
            }
            conflicts
        })
        .sum::<usize>();

    println!(
        "sampled {}-coloring with {} MCMC steps, {:.1}% overhead on {} threads",
        unif_colors,
        samples,
        100.0 * conflicts as f64 / samples as f64,
        nthreads
    );

    let mut colors: Vec<_> = colors.into_iter().map(|x| x.into_inner()).collect();

    let mut ctr = unif_colors;
    for (c, m) in colors.iter_mut().zip(mask.iter().copied()) {
        if !m {
            *c = ctr;
            ctr += 1;
        }
    }

    (ctr, colors)
}
/// Sample a coloring uniformly (and in parallel), with respect to the given graph.
///
/// Accepts a mask to operate on the induced subgraph of the given graph. Vertices
/// that are masked out are given unique colors. Uses `unif_colors` for the coloring,
/// which should be at least the number of colors used for a greedy coloring.
Sample a coloring uniformly (and in parallel), with respect to the given graph.
Accepts a mask to operate on the induced subgraph of the given graph. Vertices
that are masked out are given unique colors. Uses `unif_colors` for the coloring,
which should be at least the number of colors used for a greedy coloring.
Sample a coloring uniformly ( and in parallel ) with respect to the given graph . Accepts a mask to operate on the induced subgraph of the given graph . Vertices that are masked out are given unique colors . Uses ` unif_colors ` for the coloring which should be at least the number of colors used for a greedy coloring .
fn glauber_color(graph: &AdjacencyList, unif_colors: u32, mask: &[bool]) -> (u32, Vec<u32>) {
    let (ncolors, mut colors) = greedy_color(&graph);
    assert!(ncolors < unif_colors);
    // https://www.math.cmu.edu/~af1p/Texfiles/colorbdd.pdf
    // https://www.math.cmu.edu/~af1p/Teaching/MCC17/Papers/colorJ.pdf
    // run glauber markov chain on a coloring
    // chain sampling can be parallel with some simple conflict detection
    let max_degree = graph.max_degree(mask);
    let avg_degree = graph.avg_degree(mask);
    let subgraph_sz = mask.iter().copied().filter(|x| *x).count();
    // by expectation version of birthday problem, don't want more than this
    // much parallelism
    let max_parallel = ((subgraph_sz as f64).sqrt() / avg_degree as f64) as usize;
    let nthreads = (rayon::current_num_threads() as usize).min(max_parallel.max(16));
    println!(
        "chose parallelism {} (recommendation {}) for {}-coloring",
        nthreads, max_parallel, unif_colors
    );
    use std::sync::atomic::{AtomicU32, Ordering};
    let samples = max_degree * subgraph_sz;
    let colors = colors
        .into_iter()
        .map(|x| AtomicU32::new(x))
        .collect::<Vec<_>>();
    let conflicts = (0..nthreads)
        .into_par_iter()
        .map(|i| {
            let samples = (samples / nthreads).max(1);
            let mut rng = StdRng::seed_from_u64((1234 + i) as u64); // no idea how to do this better in rust
            let mut adjacent_color: Vec<bool> = vec![false; unif_colors as usize];
            let mut neighbor_colors: Vec<u32> = Vec::with_capacity(max_degree as usize);
            let mut conflicts = 0;
            for _ in 0..samples {
                // loop invariant is that none of adjacent_color elements are true
                let reservation = (unif_colors + 1) as u32;
                loop {
                    let v: u32 = rng.gen_range(0, graph.nvertices() as u32);
                    if !mask[v as usize] {
                        continue;
                    }
                    // "claim" this vertex
                    let prev = colors[v as usize].swap(reservation, Ordering::Relaxed);
                    if prev == reservation {
                        conflicts += 1;
                        continue;
                    }
                    let mut nadjacent_colors = 0;
                    let mut fail = v;
                    for w in graph.masked_neighbors(v, mask) {
                        let c = colors[w as usize].swap(reservation, Ordering::Relaxed);
                        if c == reservation {
                            fail = w;
                            break;
                        }
                        if !adjacent_color[c as usize] {
                            nadjacent_colors += 1;
                            adjacent_color[c as usize] = true;
                        }
                        neighbor_colors.push(c);
                    }
                    if fail != v {
                        // unlock
                        colors[v as usize].store(prev, Ordering::Relaxed);
                        graph
                            .masked_neighbors(v, mask)
                            .take_while(|&w| w != fail)
                            .zip(neighbor_colors.iter().copied())
                            .for_each(|(w, c)| colors[w as usize].store(c, Ordering::Relaxed));
                        conflicts += 1;
                        for i in adjacent_color.iter_mut() {
                            *i = false;
                        }
                        neighbor_colors.clear();
                        continue;
                    }
                    let chosen = adjacent_color
                        .iter()
                        .copied()
                        .enumerate()
                        .filter(|(_, x)| !x)
                        .map(|(i, _)| i)
                        .nth(rng.gen_range(0, (unif_colors as usize) - nadjacent_colors))
                        .expect("nth");
                    graph
                        .masked_neighbors(v, mask)
                        .zip(neighbor_colors.iter().copied())
                        .for_each(|(w, c)| colors[w as usize].store(c, Ordering::Relaxed));
                    colors[v as usize].store(chosen as u32, Ordering::Relaxed);
                    for i in adjacent_color.iter_mut() {
                        *i = false;
                    }
                    neighbor_colors.clear();
                    break;
                }
            }
            conflicts
        })
        .sum::<usize>();
    println!(
        "sampled {}-coloring with {} MCMC steps, {:.1}% overhead on {} threads",
        unif_colors,
        samples,
        100.0 * conflicts as f64 / samples as f64,
        nthreads
    );
    let mut colors: Vec<_> = colors.into_iter().map(|x| x.into_inner()).collect();
    let mut ctr = unif_colors;
    for (c, m) in colors.iter_mut().zip(mask.iter().copied()) {
        if !m {
            *c = ctr;
            ctr += 1;
        }
    }
    (ctr, colors)
}
fn glauber_color ( graph : & AdjacencyList , unif_colors : u32 , mask : & [ bool ] ) -> ( u32 , Vec < u32 > ) { let ( ncolors , mut colors ) = greedy_color ( & graph ) ; assert ! ( ncolors unif_colors ) ; // https://www.math.cmu.edu/~af1p/Texfiles/colorbdd.pdf // https://www.math.cmu.edu/~af1p/Teaching/MCC17/Papers/colorJ.pdf // run glauber markov chain on a coloring // chain sampling can be parallel with some simple conflict detection let max_degree = graph . max_degree ( mask ) ; let avg_degree = graph . avg_degree ( mask ) ; let subgraph_sz = mask . iter ( ) . copied ( ) . filter ( | x | * x ) . count ( ) ; // by expectation version of birthday problem, don't want more than this // much parallelism let max_parallel = ( ( subgraph_sz as f64 ) . sqrt ( ) / avg_degree as f64 ) as usize ; let nthreads = ( rayon :: current_num_threads ( ) as usize ) . min ( max_parallel . max ( 16 ) ) ; println ! ( " " nthreads max_parallel unif_colors ) ; use std :: sync :: atomic :: { AtomicU32 , Ordering } ; let samples = max_degree * subgraph_sz ; let colors = colors . into_iter ( ) . map ( | x | AtomicU32 :: new ( x ) ) . collect :: < Vec < _ > > ( ) ; let conflicts = ( 0 .. nthreads ) . into_par_iter ( ) . map ( | i | { let samples = ( samples / nthreads ) . max ( 1 ) ; let mut rng = StdRng :: seed_from_u64 ( ( 1234 + i ) as u64 ) ; // no idea how to do this better in rust let mut adjacent_color : Vec < bool > = vec ! [ false unif_colors as usize ] ; let mut neighbor_colors : Vec < u32 > = Vec :: with_capacity ( max_degree as usize ) ; let mut conflicts = 0 ; for _ in 0 .. samples { // loop invariant is that none of adjacent_color elements are true let reservation = ( unif_colors + 1 ) as u32 ; loop { let v : u32 = rng . gen_range ( 0 , graph . nvertices ( ) as u32 ) ; if ! mask [ v as usize ] { continue ; } // "claim" this vertex let prev = colors [ v as usize ] . swap ( reservation , Ordering :: Relaxed ) ; if prev == reservation { conflicts += 1 ; continue ; } let mut nadjacent_colors = 0 ; let mut fail = v ; for w in graph . masked_neighbors ( v , mask ) { let c = colors [ w as usize ] . swap ( reservation , Ordering :: Relaxed ) ; if c == reservation { fail = w ; break ; } if ! adjacent_color [ c as usize ] { nadjacent_colors += 1 ; adjacent_color [ c as usize ] = true ; } neighbor_colors . push ( c ) ; } if fail != v { // unlock colors [ v as usize ] . store ( prev , Ordering :: Relaxed ) ; graph . masked_neighbors ( v , mask ) . take_while ( | & w | w != fail ) . zip ( neighbor_colors . iter ( ) . copied ( ) ) . for_each ( | ( w , c ) | colors [ w as usize ] . store ( c , Ordering :: Relaxed ) ) ; conflicts += 1 ; for i in adjacent_color . iter_mut ( ) { * i = false ; } neighbor_colors . clear ( ) ; continue ; } let chosen = adjacent_color . iter ( ) . copied ( ) . enumerate ( ) . filter ( | ( _ , x ) | ! x ) . map ( | ( i , _ ) | i ) . nth ( rng . gen_range ( 0 , ( unif_colors as usize ) - nadjacent_colors ) ) . expect ( " " ) ; graph . masked_neighbors ( v , mask ) . zip ( neighbor_colors . iter ( ) . copied ( ) ) . for_each ( | ( w , c ) | colors [ w as usize ] . store ( c , Ordering :: Relaxed ) ) ; colors [ v as usize ] . store ( chosen as u32 , Ordering :: Relaxed ) ; for i in adjacent_color . iter_mut ( ) { * i = false ; } neighbor_colors . clear ( ) ; break ; } } conflicts } ) . sum :: < usize > ( ) ; println ! ( " " unif_colors samples 100.0 conflicts as f64 samples as f64 nthreads ) ; let mut colors : Vec < _ > = colors . into_iter ( ) . map ( | x | x . into_inner ( ) ) . collect ( ) ; let mut ctr = unif_colors ; for ( c , m ) in colors . iter_mut ( ) . zip ( mask . iter ( ) . copied ( ) ) { if ! m { * c = ctr ; ctr += 1 ; } } ( ctr , colors ) }
Sample a coloring uniformly (and in parallel), with respect to the given graph.
Sample a coloring uniformly ( and in parallel ) with respect to the given graph .
graph
&AdjacencyList
unif_colors
u32
mask
&[bool]
graph
&AdjacencyList
unif_colors
u32
mask
&[bool]