main.rs

extern crate rand;
extern crate pbr;

use rand::distributions::{IndependentSample, Range};
use pbr::ProgressBar;

// N = size of the board
const N: usize = 4;
// highest temperature for parallel tempering
const T_HIGH: f32 = 20f32;
// lowest temperature for parallel tempering
const T_LOW: f32 = 0.1f32;
// Number of images
const NT: usize = 20;
// number of moves to give up finding a solution after
const MAXMOVES: usize = 100000000;
// Probability of swapping two boards
const PSWAP: f32 = 0.05;

struct Queen {
    x: usize,
    y: usize,
}

fn main() {
    let mut totsteps: usize = 0;
    // Create RNG stuff here and pass it to run, so we don't
    // re-initialize RNG for every solution. This might not be
    // necessary.
    let mut rng = rand::weak_rng();
    let nrange = Range::new(0, N);
    let ntrange = Range::new(0, NT);
    let probrange = Range::new(0f32, 1f32);

    // Number of minimizations to perform
    let trials = 200000;

    let mut pb = ProgressBar::new(trials);
    for i in 0..trials {
        totsteps += run(&nrange, &ntrange, &probrange, &mut rng);
        pb.message(format!("Average number of steps: {} | ", totsteps as f32/((i+1) as f32)).as_str());
        pb.inc();
    };
    pb.finish_println("Done!\n");
}

fn calc_e(rows: &Vec<isize>, columns: &Vec<isize>, pdiags: &Vec<isize>, ndiags: &Vec<isize>) -> isize {
    let mut e: isize = 0;
    for row in rows {
        e += (row - 1).pow(2);
    }
    for column in columns {
        e += (column - 1).pow(2);
    }
    for pdiag in pdiags {
        e += (pdiag - 1).pow(2) * ((*pdiag != 0) as isize);
    }
    for ndiag in ndiags {
        e += (ndiag - 1).pow(2) * ((*ndiag != 0) as isize);
    }
    return e;
}

fn calc_de(rows: &Vec<isize>, columns: &Vec<isize>, pdiags: &Vec<isize>, ndiags: &Vec<isize>, oldqueen: &Queen, newqueen: &Queen) -> isize {
    let mut de: isize = 0;
    assert!(rows[oldqueen.x] >= 1);
    assert!(columns[oldqueen.y] >= 1);
    assert!(pdiags[oldqueen.x + oldqueen.y] >= 1);
    assert!(pdiags[(oldqueen.x + N) - (oldqueen.y + 1)] >= 1);
    de -= 2*rows[oldqueen.x] - 3;
    de -= 2*columns[oldqueen.x] - 3;
    de -= (2*pdiags[oldqueen.x + oldqueen.y] - 3) * (pdiags[oldqueen.x + oldqueen.y] != 1) as isize;
    de -= (2*ndiags[(oldqueen.x + N) - (oldqueen.y + 1)] - 3) * (ndiags[(oldqueen.x + N) - (oldqueen.y + 1)] != 1) as isize;
    de += 2*rows[newqueen.x] - 1;
    de += 2*columns[newqueen.x] - 1;
    de += (2*pdiags[newqueen.x + newqueen.y] - 1) * (pdiags[newqueen.x + newqueen.y] != 0) as isize;
    de += (2*ndiags[(newqueen.x + N) - (newqueen.y + 1)] - 1) * (ndiags[(newqueen.x + N) - (newqueen.y + 1)] != 0) as isize;
    return de;
}


fn run(nrange: &Range<usize>, ntrange: &Range<usize>, probrange: &Range<f32>, rng: &mut rand::XorShiftRng) -> usize {

    // Temperature increase factor between each adjacent image
    let tfact: f32 = (T_HIGH / T_LOW).powf((NT as f32).powi(-1));
    // Vec of beta values
    let mut bs: Vec<f32> = Vec::new();

    // boards is a Vec of Vecs of Queens
    let mut boards: Vec<Vec<Queen>> = Vec::new();
    let mut queens: Vec<Queen>;

    // "Energy" of the system. A solution will have e==0.
    let mut e: Vec<isize> = Vec::new();
    
    let mut rows_all: Vec<Vec<isize>> = Vec::new();
    let mut columns_all: Vec<Vec<isize>> = Vec::new();
    let mut pdiags_all: Vec<Vec<isize>> = Vec::new();
    let mut ndiags_all: Vec<Vec<isize>> = Vec::new();

    let mut rows: Vec<isize>;
    let mut columns: Vec<isize>;
    let mut pdiags: Vec<isize>;
    let mut ndiags: Vec<isize>;

    // Create NT NxN boards, each with N randomly placed queens
    for i in 0..NT {
        bs.push((T_LOW * tfact.powi(i as i32)).powi(-1));
        queens = Vec::new();
        rows = Vec::new();
        columns = Vec::new();
        pdiags = Vec::new();
        ndiags = Vec::new();
        for _ in 0..N {
            rows.push(0);
            columns.push(0);
        }
        for _ in 0..2*N-1 {
            pdiags.push(0);
            ndiags.push(0);
        }
        for _ in 0..N {
            let x = nrange.ind_sample(rng);
            let y = nrange.ind_sample(rng);
            queens.push(Queen { x: x, y: y});
            rows[x] += 1;
            columns[y] += 1;
            pdiags[x + y] += 1;
            ndiags[(x + N) - (y + 1)] += 1;
        };
        e.push(calc_e(&rows, &columns, &pdiags, &ndiags));
        println!{"energy: {}, board number: {}", e[i], i};
        boards.push(queens);
        rows_all.push(rows);
        columns_all.push(columns);
        pdiags_all.push(pdiags);
        ndiags_all.push(ndiags);
    };

    // Change of energy associated with a trial move
    let mut de: isize;
    // Change in beta value, used for board swapping
    let mut db: f32;

    for j in 1..MAXMOVES {
        // Pick a board to work with
        let bn = ntrange.ind_sample(rng);

        // Check to see if we're switching boards
        if probrange.ind_sample(rng) < PSWAP {
            // If we picked the last board and are swapping boards, just do
            // nothing and continue with the next loop. This effectively lowers
            // our probability of swapping boards, but it makes the code a lot
            // simpler.
            if bn == NT - 1 { continue };
            de = e[bn + 1] - e[bn];
            db = bs[bn + 1] - bs[bn];
            if (de as f32 * db < 0f32) || (probrange.ind_sample(rng) < (-de as f32 * db).exp()) {
                e.swap(bn, bn+1);
                boards.swap(bn, bn+1);
                rows_all.swap(bn, bn+1);
                columns_all.swap(bn, bn+1);
                pdiags_all.swap(bn, bn+1);
                ndiags_all.swap(bn, bn+1);
            };
            continue;
        };

        // Pick a random queen to move and a new position
        let qn = nrange.ind_sample(rng);
        let x = nrange.ind_sample(rng);
        let y = nrange.ind_sample(rng);
        let newqueen = Queen { x: x, y: y };
        
        de = calc_de(&rows_all[bn], &columns_all[bn], &pdiags_all[bn], &ndiags_all[bn], &boards[bn][qn], &newqueen);
        
        //println!{"{}", de};
        // Metropolis-Hastings algorithm. If de <= 0, accept the move.
        // Otherwise, accept the move with probability given by e^(-dE/kT)
        // (so, bigger de => less likely to accept move). Here we replace
        // 1/kT with B.
        // If we accept the move, replace the old queen with the new one
        // and update the total system energy.

        if (de <= 0) || (probrange.ind_sample(rng) < (-de as f32 * bs[bn]).exp()) {
            e[bn] += de;

            // If the energy is 0, we're done.
            if e[bn] == 0 {
                return j;
            };

            rows_all[bn][x] += 1;
            columns_all[bn][y] += 1;
            pdiags_all[bn][x + y] += 1;
            ndiags_all[bn][(x + N) - (y + 1)] += 1;

            rows_all[bn][boards[bn][qn].x] -= 1;
            columns_all[bn][boards[bn][qn].y] -= 1;
            pdiags_all[bn][&boards[bn][qn].x + &boards[bn][qn].y] -= 1;
            ndiags_all[bn][(&boards[bn][qn].x + N) - (&boards[bn][qn].y + 1)] -= 1;

            boards[bn][qn] = newqueen;
        };
    };
    // Rather than actually handling the failure case, just return MAXMOVES.
    // This should not affect the results unless the board is very big.
    return MAXMOVES;
}