Zenn
😺

Rustで焼きなまし法のライブラリを作ってみた

2024/12/13に公開
Rust
tech

TL;DR;

この記事では焼きなまし法のPythonライブラリsimannealの遅さを克服するために作ったRust製ライブラリrusty-simannealを開発した経緯とその内容を紹介しています。
Rustで実装して早くしましただけでは無く、高速化の工夫（状態復元や効率的なエネルギー計算）を取り入れ、TSP問題で実装にRustのナイーブ実装と比べ最大3倍程度の速度向上を実現した話も書いてあります。

はじめに

焼きなまし法とは

焼きなまし法（Simulated Annealing, SA）は、最適化問題を解くための確率的手法の一つで、特に組み合わせ最適化問題で広く用いられます。焼きなましという名前は、金属の熱処理過程「焼きなまし（annealing）」に由来しています。

焼きなまし法の基本的な仕組み

  1. 初期状態の設定:

    • 問題の初期解をランダムまたは適切な方法で生成します。

  2. エネルギー関数の定義:

    • 最適化問題において、解の品質を評価するための目的関数をエネルギー関数として定義します。

  3. 温度（T）の設定:

    • 温度は探索の許容範囲を制御するパラメータで、高いほど広範囲を探索し、低いほど局所探索に絞られます。

  4. 遷移の試行:

    • 現在の解を基に、ランダムに近傍解を生成します。
    • 新しい解がエネルギー関数の値（目的関数）を改善する場合はその解を採用します。
    • 改善しない場合でも、以下の確率 ( P = \exp(-\Delta E / T) ) に基づいて受け入れることがあります。
      • ( \Delta E ) はエネルギー差、新解のエネルギー - 現解のエネルギー
      • ( T ) は温度

  5. 温度の低下:

    • 一定のルール（冷却スケジュール）に従い、温度を徐々に下げます。これにより、初期は広範囲を探索し、最終的には局所解に集中します。

  6. 停止条件の確認:

    • 一定回数の反復、温度が閾値以下になる、改善が見られなくなるなどの条件で探索を終了します。

  7. 最適解の出力:

    • 最終的に得られた解を出力します。

動機

仕事でsimannealというPyhtonのライブラリを利用していました。

このライブラリは以下のように書くと

from simanneal import Annealer

class TravellingSalesmanProblem(Annealer):
    """Test annealer with a travelling salesman problem."""

    def move(self):
        """Swaps two cities in the routes"""
        a = random.randint(0, len(self.state) - 1)
        b = random.randint(0, len(self.state) - 1)
        self.state[a], self.state[b] = self.state[b], self.state[a]
    
    def energy(self):
        """Calculates the length of the route."""
        e = 0
        for i in range(len(self.state)):
            e += self.distance(cities[self.state[i - 1]],
                          cities[self.state[i]])
        return e

比較的簡単に焼きなまし法を利用できる便利なライブラリなのですが以下のような問題がありました。

長い間メンテナンスがされていない

標準エラー出力がopt-out出来ない等の細かい問題がありforkして利用している状況でした。

単純に遅い

焼きなまし法は説明の通り、解を少しずつ変えながら反復的に解を改善していくアルゴリズムです。
実行速度が遅い場合は、時間を固定したときに解の質が悪くなることになります。

ここで、Pythonの実行速度の遅さが仇となっていました。

ということもあり、書き味のシンプルさを残したままRustに書き換えました。

設計

焼きなまし法の構成要素を分解すると以下のようになります。

解の状態

現在の解

Transtion

現在の解を近傍解に移動させるための要素

フローチャート

この解の状態とを表すtraitを用意して、ユーザ側はそれを実装するような方針にします。

具体的には

/// 解の評価が出来ることを表すtrait
pub trait EnergyMeasurable: Sized + Clone + Debug {
    type Energy: PartialOrd + Clone + Copy + Debug + Sub + Signed + Num + Neg + Into<f64>;
    type Context;

    fn energy(&self, ctx: &Self::Context) -> Self::Energy;
}

/// 乱数から初期解を作成するtrait
pub trait InitialState {
    type Context;

    fn initial_state<G: Rng>(&self, rng: &mut G, ctx: &Self::Context) -> Self;
}

/// Transitionを受け取って近傍解へ遷移する関数
pub trait AnnealingState: EnergyMeasurable {
    type Transition: Transition<Context = Self::Context, State = Self> + Debug;

    fn apply(&mut self, ctx: &Self::Context, op: &Self::Transition) -> Option<()>;
}

というtraitと

/// 乱数からTransitionを生成するtrait
pub trait Transition: Sized + Clone + Copy {
    type Context;
    type State;

    fn choose<G: Rng>(rng: &mut G, ctx: &Self::Context, state: &Self::State) -> Self;
}

というtraitを作る。

その上で

/// Simulated Annealing algorithm
/// minimize f(x) where x is a state
pub struct Annealer<S: EnergyMeasurable, C: Schedule> {
    pub state: S,
    pub ctx: S::Context,
    pub schedule: C,
    pub metrics: Vec<Metrics>,
}

impl<S: AnnealingState, C: Schedule> Annealer<S, C> {
    pub fn anneal<G: Rng>(&mut self, rng: &mut G) -> S {
        let mut best_state = self.state.clone();
        let mut best_energy = self.state.energy(&self.ctx);
        let mut current_energy = best_energy;
        let mut progress = Progress::zero();

        // 探索ループ
        while self.schedule.should_continue(&progress) {
            let prev_state = self.state.clone();
            let op = S::Transition::choose(rng, &self.ctx, &self.state);

            let (accept, improvement) = if let Some(_restore) = self.state.apply(&self.ctx, &op) {
                let temperature = self.schedule.temperature(&progress);
                let new_energy = self.state.energy(&self.ctx);

                let improvement = if new_energy < best_energy {
                    best_energy = new_energy;
                    best_state = self.state.clone();
                    true
                } else {
                    false
                };

                let delta = (new_energy - current_energy).into();
                let p = rng.gen_range(0.0..=1.0);
                // 乱数から遷移するかを判定する
                if delta.is_sign_positive() && (-delta / temperature).exp() < p {
                    // reject
                    debug!("reject {} -> {}", current_energy.into(), new_energy.into());
                    self.state = prev_state;
                    (false, improvement)
                } else {
                    // accept
                    debug!("accept {} -> {}", current_energy.into(), new_energy.into());
                    current_energy = new_energy;
                    (true, improvement)
                }
            } else {
                (false, false)
            };

            progress.update();
        }

        best_state
    }
}

というような、上のtraitの実装をまとめるような構造体と焼きなまし法を実行する関数を作成します。

使い方

簡単な2次関数 f(x) = a x^2 + b x + c の最小となるxを見つける問題を定義してみます。

// solve f(x) = a x^2 + b x + c
#[derive(Debug, Clone)]
struct QuadraticFunction {
    a: f64,
    b: f64,
    c: f64,
}

#[derive(Debug, Clone)]
struct QuadraticFunctionState {
    x: f64,
}
impl EnergyMeasurable for QuadraticFunctionState {
    type Energy = f64;
    type Context = QuadraticFunction;

    fn energy(&self, ctx: &Self::Context) -> Self::Energy {
        let f = ctx.a * self.x * self.x + ctx.b * self.x + ctx.c;
        f * ctx.a.signum()
    }
}

impl AnnealingState for QuadraticFunctionState {
    type Transition = QuadraticFunctionTransition;

    fn apply(&mut self, _ctx: &Self::Context, op: &Self::Transition) -> Option<()> {
        match op {
            QuadraticFunctionTransition::Add(x) => {
                self.x += x;
                Some(())
            }
            QuadraticFunctionTransition::Mul(x) => {
                self.x *= x;
                Some(())
            }
        }
    }
}

#[derive(Debug, Clone, Copy)]
enum QuadraticFunctionTransition {
    Add(f64),
    Mul(f64),
}

impl Transition for QuadraticFunctionTransition {
    type Context = QuadraticFunction;
    type State = QuadraticFunctionState;

    fn choose<G: Rng>(rng: &mut G, _ctx: &Self::Context, _state: &Self::State) -> Self {
        match rng.gen_range(0..=1) {
            0 => Self::Add(rng.gen_range(-10.0..=10.0)),
            1 => Self::Mul(rng.gen_range(0.3..=1.1)),
            _ => unreachable!(),
        }
    }
}

ちなみに、Pythonだと以下のような感じです。
少し長くなっちゃいましたね。

Python実装
from simanneal import Annealer

class QuadraticFunctionState(Annealer):
    def __init__(self, a, b, c):
        self.a = a
        self.b = b
        self.c = c
        self.x = rando.uniform(-100, 100)

    def move(self):
        add_mul = random.randint(0, len(self.state) - 1)
        if add_mul == 0:
            self.x += random.uniform(-10.0, 10.0)
        if add_mul == 1:
            self.x *= random.uniform(3.0, 1.1)

    
    def energy(self):
        return a * x ** 2 + b * x + c

探索を開始するには

fn main() {
     let mut annealer = Annealer::new(
        QuadraticFunctionState { x: 100.0 },
        QuadraticFunction {
            a: 1.0,
            b: 10.0,
            c: 30.0,
        },
        schedule::LinearStepSchedule::new(1000.0, 0.01, 10000),
    );

    let state = annealer.anneal::<_, false>(&mut rand::thread_rng());
}

というような感じで実行します。
このプログラムだと10000回ループを回して実行します。

高速化

序盤にも書いたように、焼きなまし法の実装において高速化は大事です。

rusty-simannealの動機にも高速な実行があったので高速化を行っていきます。

現在の実装の問題点

上記のサンプルの場合は状態はxという1つの変数になっているので問題になりにくいですが、実際はStateを表す構造体はある程度のサイズがあることが一般的です。

例として、よく例としてあげられるTSPのStateは以下のようになっていて頂点数と同じだけの配列を持ちます。

#[derive(Debug, Clone)]
struct TspState {
    route: Vec<usize>,
}

現在の実装だと

...
while self.schedule.should_continue(&progress) {
    let prev_state = self.state.clone();
    let op = S::Transition::choose(rng, &self.ctx, &self.state);
...

の部分で状態をcloneしているため重くなってしまうという問題があります。

現在のtraitの仕組みではこのclone回避できないためユーザ側が頑張って実装すればこのコストを回避できる方法を用意していきます。
このアイディアは高速なビームサーチが欲しい!!!を参考にしています。

そのtraitが

  • AnnealingStateBack
  • AnnealingStatePeeking

になります。
ユーザ側はこの二つのうちのどちらかを実装するとcloneを本質的に避けて焼きなまし法を実行できるようになります。

AnnealingStateBack

このtraitはTransitionで遷移するときに(比較的に小さい)Restoreと呼ばれる解を元に戻すための情報を出力する遷移関数とRestoreを受け取って解の状態を元に戻す関数を実装します。

/// AnnealingStateBack is implemented when the state can be back processed more efficiently than clone.
pub trait AnnealingStateBack: AnnealingState {
    type Restore;

    /// Apply the transition and restore information for returning to the previous state
    fn apply_with_restore(
        &mut self,
        ctx: &Self::Context,
        op: &Self::Transition,
    ) -> Option<Self::Restore>;

    /// Restore the state to the previous state
    fn back(&mut self, ctx: &Self::Context, restore: &Self::Restore);
}

これがState実装されている場合は以下のような実装が追加されanneal_backが実行できるようになります。
anneal_backはcloneを利用せずに前述のRestoreを用いて状態を元に戻すことでcloneを避けることが可能になっています。

impl<S: AnnealingStateBack, C: Schedule> Annealer<S, C> {
    /// Simulated Annealing algorithm
    /// minimize f(x) where x is a state
    /// Use BACK instead of CLONE when you want to abort and back to the state.
    pub fn anneal_back<G: Rng, const METRICS: bool>(&mut self, rng: &mut G) -> S {
        let mut best_state = self.state.clone();
        let mut best_energy = self.state.energy(&self.ctx);
        let mut current_energy = best_energy;
        let mut progress = Progress::zero();

        while self.schedule.should_continue(&progress) {
            let op = Transition::choose(rng, &self.ctx, &self.state);
            if let Some(restore) = self.state.apply_with_restore(&self.ctx, &op) {
                let temperature = self.schedule.temperature(&progress);
                let new_energy = self.state.energy(&self.ctx);
                let delta = (new_energy - current_energy).into();
                let p = rng.gen_range(0.0..=1.0);
                if delta.is_sign_positive() && (-delta / temperature).exp() < p {
                    self.state.back(&self.ctx, &restore);
                } else {
                    current_energy = new_energy;
                    if current_energy < best_energy {
                        best_energy = current_energy;
                        best_state = self.state.clone();
                    }
                }
            }
            progress.update();
        }

        best_state
    }
}

AnnealingStatePeeking

このtraitは状態を遷移させなくても次のenergyを計算することが出来る場合に実装します。

/// AnnealingStatePeeking is a trait to be implemented when the energy of the next state can be calculated efficiently without updating the state.
pub trait AnnealingStatePeeking: AnnealingState {
    /// Peek the energy of the state after applying the transition
    fn peek_energy(
        &self,
        ctx: &Self::Context,
        op: &Self::Transition,
        current_energy: Self::Energy,
    ) -> Option<Self::Energy>;
}

これがState実装されている場合は以下のような実装が追加されanneal_peekが実行できるようになります。
anneal_peekでは、次の状態に遷移させずにenergyを計算して判定関数に掛け必要な場合にのみ遷移させることによりcloneを削減させることが出来ています。

impl<S: AnnealingStatePeeking, C: Schedule> Annealer<S, C> {
    /// Simulated Annealing algorithm
    /// minimize f(x) where x is a state
    /// Use peek_energy instead of apply when the energy of the next state can be calculated efficiently without updating the state.
    pub fn anneal_peek<G: Rng, const METRICS: bool>(&mut self, rng: &mut G) -> S {
        let mut best_state = self.state.clone();
        let mut best_energy = self.state.energy(&self.ctx);
        let mut current_energy = best_energy;
        let mut progress = Progress::zero();

        while self.schedule.should_continue(&progress) {
            let op = Transition::choose(rng, &self.ctx, &self.state);
            if let Some(new_energy) = self.state.peek_energy(&self.ctx, &op, current_energy) {
                let temperature = self.schedule.temperature(&progress);
                let delta = (new_energy - current_energy).into();
                let p = rng.gen_range(0.0..=1.0);
                if !(delta.is_sign_positive() && (-delta / temperature).exp() < p) {
                    // accept
                    self.state.apply(&self.ctx, &op);
                    current_energy = new_energy;
                    if current_energy < best_energy {
                        best_energy = current_energy;
                        best_state = self.state.clone();
                    }
                }
            }
            progress.update();
        }

        best_state
    }
}

速度検証

試しによくある47都道府県の直線距離TSPの問題で(2-opt)速度を検証していきます。

実装の全体
use std::time::Instant;

use rand::{Rng, SeedableRng};
use rand::rngs::SmallRng;
use serde::Deserialize;

use rusty_simanneal::{
    AnnealingState, AnnealingStateBack, AnnealingStatePeeking, EnergyMeasurable, Transition,
};
use rusty_simanneal::schedule::LinearStepSchedule;

// from simanneal import Annealer
// class TravellingSalesmanProblem(Annealer):
//
//     """Test annealer with a travelling salesman problem.
//     """
//
//     # pass extra data (the distance matrix) into the constructor
//     def __init__(self, state, distance_matrix):
//         self.distance_matrix = distance_matrix
//         super(TravellingSalesmanProblem, self).__init__(state)  # important!
//
//     def move(self):
//         """Swaps two cities in the route."""
//         # no efficiency gain, just proof of concept
//         # demonstrates returning the delta energy (optional)
//         initial_energy = self.energy()
//
//         a = random.randint(0, len(self.state) - 1)
//         b = random.randint(0, len(self.state) - 1)
//         self.state[a], self.state[b] = self.state[b], self.state[a]
//
//         return self.energy() - initial_energy
//
//     def energy(self):
//         """Calculates the length of the route."""
//         e = 0
//         for i in range(len(self.state)):
//             e += self.distance_matrix[self.state[i-1]][self.state[i]]
//         return e

#[derive(Deserialize, Debug, Clone)]
struct Row {
    #[serde(rename = "Town")]
    town: String,
    #[serde(rename = "Longitude")]
    longitude: f64,
    #[serde(rename = "Latitude")]
    latitude: f64,
}

#[derive(Debug, Clone)]
struct TspContext {
    distance_matrix: Vec<Vec<f64>>,
}

#[derive(Debug, Clone)]
struct TspState {
    route: Vec<usize>,
}

#[derive(Debug, Clone, Copy)]
struct TspTransition {
    a: usize,
    b: usize,
}

impl EnergyMeasurable for TspState {
    type Energy = f64;
    type Context = TspContext;

    fn energy(&self, ctx: &Self::Context) -> f64 {
        let mut e = 0.0;
        for i in 1..self.route.len() {
            e += ctx.distance_matrix[self.route[i - 1]][self.route[i]];
        }
        e += ctx.distance_matrix[self.route[self.route.len() - 1]][self.route[0]];
        e
    }
}

impl AnnealingState for TspState {
    type Transition = TspTransition;

    fn apply(&mut self, _ctx: &Self::Context, op: &Self::Transition) -> Option<()> {
        self.route.swap(op.a, op.b);
        Some(())
    }
}

impl AnnealingStateBack for TspState {
    type Restore = Self::Transition;

    fn apply_with_restore(
        &mut self,
        _ctx: &Self::Context,
        op: &Self::Transition,
    ) -> Option<Self::Restore> {
        let restore = *op;
        self.route.swap(op.a, op.b);
        Some(restore)
    }

    fn back(&mut self, _ctx: &Self::Context, restore: &Self::Restore) {
        self.route.swap(restore.a, restore.b);
    }
}

impl AnnealingStatePeeking for TspState {
    fn peek_energy(
        &self,
        ctx: &Self::Context,
        op: &Self::Transition,
        current_energy: Self::Energy,
    ) -> Option<Self::Energy> {
        let mut e = current_energy;

        let a = i64::try_from(op.a).unwrap();
        let b = i64::try_from(op.b).unwrap();

        let (a, b) = if a < b { (a, b) } else { (b, a) };
        let (a, b) = if a == 0 && b == (self.route.len() as i64 - 1) {
            (b, a)
        } else {
            (a, b)
        };

        let prev_a = (a - 1).rem_euclid(self.route.len() as i64) as usize;
        let next_a = (a + 1).rem_euclid(self.route.len() as i64) as usize;
        let prev_b = (b - 1).rem_euclid(self.route.len() as i64) as usize;
        let next_b = (b + 1).rem_euclid(self.route.len() as i64) as usize;

        if (a - b).abs().min(b + self.route.len() as i64 - a).abs() > 1 {
            e -= ctx.distance_matrix[self.route[prev_a]][self.route[a as usize]];
            e -= ctx.distance_matrix[self.route[a as usize]][self.route[next_a]];
            e -= ctx.distance_matrix[self.route[prev_b]][self.route[b as usize]];
            e -= ctx.distance_matrix[self.route[b as usize]][self.route[next_b]];

            e += ctx.distance_matrix[self.route[prev_a]][self.route[b as usize]];
            e += ctx.distance_matrix[self.route[b as usize]][self.route[next_a]];
            e += ctx.distance_matrix[self.route[prev_b]][self.route[a as usize]];
            e += ctx.distance_matrix[self.route[a as usize]][self.route[next_b]];
        } else {
            e -= ctx.distance_matrix[self.route[prev_a]][self.route[a as usize]];
            e -= ctx.distance_matrix[self.route[b as usize]][self.route[next_b]];

            e += ctx.distance_matrix[self.route[prev_a]][self.route[b as usize]];
            e += ctx.distance_matrix[self.route[a as usize]][self.route[next_b]];
        }

        Some(e)
    }
}

impl Transition for TspTransition {
    type Context = TspContext;
    type State = TspState;

    fn choose<G: Rng>(rng: &mut G, ctx: &Self::Context, state: &Self::State) -> Self {
        let a = rng.gen_range(0..state.route.len());
        let b = rng.gen_range(0..state.route.len());
        Self { a, b }
    }
}

fn main() {
    let mut reader = include_bytes!("data/location.txt").to_vec();
    let rdr = csv::ReaderBuilder::new().from_reader(&reader[..]);

    let rows: Result<Vec<_>, _> = rdr.into_deserialize::<Row>().collect();
    let rows = rows.unwrap();

    let mut distance_matrix = vec![vec![0.0; rows.len()]; rows.len()];
    for i in 0..rows.len() {
        for j in 0..rows.len() {
            let x = rows[i].longitude - rows[j].longitude;
            let y = rows[i].latitude - rows[j].latitude;
            distance_matrix[i][j] = (x * x + y * y).sqrt();
        }
    }

    let ctx = TspContext { distance_matrix };
    let state = TspState {
        route: (0..rows.len()).collect::<Vec<_>>().try_into().unwrap(),
    };

    {
        println!("start normal(clone) annealing");
        let state = state.clone();
        let mut rng = SmallRng::seed_from_u64(0);
        let mut annealer = rusty_simanneal::Annealer::new(
            state,
            ctx.clone(),
            LinearStepSchedule::new(100.0, 0.01, 10_000_000),
        );

        let start = Instant::now();
        let best_state = annealer.anneal::<_, false>(&mut rng);

        println!("process time {}ms", start.elapsed().as_millis());
        println!("{:?}", best_state);
        println!("{:?}", best_state.energy(&ctx));
    }

    {
        println!("start state back annealing");
        let state = state.clone();
        let mut rng = SmallRng::seed_from_u64(0);
        let mut annealer = rusty_simanneal::Annealer::new(
            state,
            ctx.clone(),
            LinearStepSchedule::new(100.0, 0.01, 10_000_000),
        );

        let start = Instant::now();
        let best_state = annealer.anneal_back::<_, false>(&mut rng);

        println!("process time {}ms", start.elapsed().as_millis());
        println!("{:?}", best_state);
        println!("{:?}", best_state.energy(&ctx));
    }

    {
        println!("start state peek energy annealing");
        let state = state.clone();
        let mut rng = SmallRng::seed_from_u64(0);
        let mut annealer = rusty_simanneal::Annealer::new(
            state,
            ctx.clone(),
            LinearStepSchedule::new(100.0, 0.01, 10_000_000),
        );

        let start = Instant::now();
        let best_state = annealer.anneal_peek::<_, true>(&mut rng);

        println!("process time {}ms", start.elapsed().as_millis());
        println!("{:?}", best_state);
        println!("{:?}", best_state.energy(&ctx));
    }
}

結果は以下のようになりました。

Annealing Method Process Time (ms)
Normal (Clone) Annealing 1136
State Back Annealing 609
State Peek Energy Annealing 334

ナイーブな実装に比べると、anneal_backが倍程度、anneal_peekがさらにその倍早いことがわかります。
ただし、実装を見てみると以下のようになっており

AnnealingState < AnnealingStateBack < AnnealingStatePeeking

の順番で実装難易度が難しくなっていることもわかります。
AnnealingStatePeeking最強と言うよりは、実験フェーズ・高速化フェーズと相談しながら使い分けていく感じになるかと思います。

impl AnnealingState for TspState {
    type Transition = TspTransition;

    fn apply(&mut self, _ctx: &Self::Context, op: &Self::Transition) -> Option<()> {
        self.route.swap(op.a, op.b);
        Some(())
    }
}

impl AnnealingStateBack for TspState {
    type Restore = Self::Transition;

    fn apply_with_restore(
        &mut self,
        _ctx: &Self::Context,
        op: &Self::Transition,
    ) -> Option<Self::Restore> {
        let restore = *op;
        self.route.swap(op.a, op.b);
        Some(restore)
    }

    fn back(&mut self, _ctx: &Self::Context, restore: &Self::Restore) {
        self.route.swap(restore.a, restore.b);
    }
}

impl AnnealingStatePeeking for TspState {
    fn peek_energy(
        &self,
        ctx: &Self::Context,
        op: &Self::Transition,
        current_energy: Self::Energy,
    ) -> Option<Self::Energy> {
        let mut e = current_energy;

        let a = i64::try_from(op.a).unwrap();
        let b = i64::try_from(op.b).unwrap();

        let (a, b) = if a < b { (a, b) } else { (b, a) };
        let (a, b) = if a == 0 && b == (self.route.len() as i64 - 1) {
            (b, a)
        } else {
            (a, b)
        };

        let prev_a = (a - 1).rem_euclid(self.route.len() as i64) as usize;
        let next_a = (a + 1).rem_euclid(self.route.len() as i64) as usize;
        let prev_b = (b - 1).rem_euclid(self.route.len() as i64) as usize;
        let next_b = (b + 1).rem_euclid(self.route.len() as i64) as usize;

        if (a - b).abs().min(b + self.route.len() as i64 - a).abs() > 1 {
            e -= ctx.distance_matrix[self.route[prev_a]][self.route[a as usize]];
            e -= ctx.distance_matrix[self.route[a as usize]][self.route[next_a]];
            e -= ctx.distance_matrix[self.route[prev_b]][self.route[b as usize]];
            e -= ctx.distance_matrix[self.route[b as usize]][self.route[next_b]];

            e += ctx.distance_matrix[self.route[prev_a]][self.route[b as usize]];
            e += ctx.distance_matrix[self.route[b as usize]][self.route[next_a]];
            e += ctx.distance_matrix[self.route[prev_b]][self.route[a as usize]];
            e += ctx.distance_matrix[self.route[a as usize]][self.route[next_b]];
        } else {
            e -= ctx.distance_matrix[self.route[prev_a]][self.route[a as usize]];
            e -= ctx.distance_matrix[self.route[b as usize]][self.route[next_b]];

            e += ctx.distance_matrix[self.route[prev_a]][self.route[b as usize]];
            e += ctx.distance_matrix[self.route[a as usize]][self.route[next_b]];
        }

        Some(e)
    }
}

まとめ

今回は自作のRust製ライブラリのrusty-annealの紹介とともに、焼きなまし法とその実装の高速化などを話していきました。

紹介したライブラリのリポジトリはこれになります。
まだ、自分の実験用に利用しているだけなのと、割と困ってないのでドキュメントが充実してないですがスターいただけると励みになりますのでよろしくお願いします！

Discussion