diff --git a/cuda_keyring.deb b/cuda_keyring.deb new file mode 100644 index 00000000..559906a9 Binary files /dev/null and b/cuda_keyring.deb differ diff --git a/tig-algorithms/src/knapsack/classic_quadkp/benchmarker_outbound.rs b/tig-algorithms/src/knapsack/classic_quadkp/benchmarker_outbound.rs new file mode 100644 index 00000000..4f65befd --- /dev/null +++ b/tig-algorithms/src/knapsack/classic_quadkp/benchmarker_outbound.rs @@ -0,0 +1,167 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Benchmarker Outbound Game License v1.0 (the "License"); you +may not use this file except in compliance with the License. You may obtain a copy +of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + +use anyhow::Result; +use tig_challenges::knapsack::{Challenge, Solution}; + +pub fn solve_challenge(challenge: &Challenge) -> Result> { + let vertex_count = challenge.weights.len(); + + let mut edge_costs: Vec<(usize, f32)> = (0..vertex_count) + .map(|flow_index| { + let total_flow = challenge.values[flow_index] as i32 + + challenge.interaction_values[flow_index].iter().sum::(); + let cost = total_flow as f32 / challenge.weights[flow_index] as f32; + (flow_index, cost) + }) + .collect(); + + edge_costs.sort_unstable_by(|a, b| b.1.partial_cmp(&a.1).unwrap()); + + let mut coloring = Vec::with_capacity(vertex_count); + let mut uncolored = Vec::with_capacity(vertex_count); + let mut current_entropy = 0; + let mut current_temperature = 0; + + for &(flow_index, _) in &edge_costs { + if current_entropy + challenge.weights[flow_index] <= challenge.max_weight { + current_entropy += challenge.weights[flow_index]; + current_temperature += challenge.values[flow_index] as i32; + + for &colored in &coloring { + current_temperature += challenge.interaction_values[flow_index][colored]; + } + coloring.push(flow_index); + } else { + uncolored.push(flow_index); + } + } + + let mut mutation_rates = vec![0; vertex_count]; + for flow_index in 0..vertex_count { + mutation_rates[flow_index] = challenge.values[flow_index] as i32; + for &colored in &coloring { + mutation_rates[flow_index] += challenge.interaction_values[flow_index][colored]; + } + } + + let max_generations = 100; + let mut cooling_schedule = vec![0; vertex_count]; + + for _ in 0..max_generations { + let mut best_mutation = 0; + let mut best_crossover = None; + + for uncolored_index in 0..uncolored.len() { + let mutant = uncolored[uncolored_index]; + if cooling_schedule[mutant] > 0 { + continue; + } + + unsafe { + let mutant_fitness = *mutation_rates.get_unchecked(mutant); + let min_entropy_reduction = *challenge.weights.get_unchecked(mutant) as i32 - (challenge.max_weight as i32 - current_entropy as i32); + + if mutant_fitness < 0 { + continue; + } + + for colored_index in 0..coloring.len() { + let gene_to_remove = *coloring.get_unchecked(colored_index); + if *cooling_schedule.get_unchecked(gene_to_remove) > 0 { + continue; + } + + if min_entropy_reduction > 0 { + let removed_entropy = *challenge.weights.get_unchecked(gene_to_remove) as i32; + if removed_entropy < min_entropy_reduction { + continue; + } + } + + let fitness_change = mutant_fitness - *mutation_rates.get_unchecked(gene_to_remove) + - *challenge.interaction_values.get_unchecked(mutant).get_unchecked(gene_to_remove); + + if fitness_change > best_mutation { + best_mutation = fitness_change; + best_crossover = Some((uncolored_index, colored_index)); + } + } + } + } + + if let Some((uncolored_index, colored_index)) = best_crossover { + let gene_to_add = uncolored[uncolored_index]; + let gene_to_remove = coloring[colored_index]; + + coloring.swap_remove(colored_index); + uncolored.swap_remove(uncolored_index); + coloring.push(gene_to_add); + uncolored.push(gene_to_remove); + + current_temperature += best_mutation; + current_entropy = current_entropy + challenge.weights[gene_to_add] - challenge.weights[gene_to_remove]; + + unsafe { + for flow_index in 0..vertex_count { + *mutation_rates.get_unchecked_mut(flow_index) += + challenge.interaction_values.get_unchecked(flow_index).get_unchecked(gene_to_add) - + challenge.interaction_values.get_unchecked(flow_index).get_unchecked(gene_to_remove); + } + } + + cooling_schedule[gene_to_add] = 3; + cooling_schedule[gene_to_remove] = 3; + } else { + break; + } + + if current_temperature as u32 >= challenge.min_value { + return Ok(Some(Solution { items: coloring })); + } + + for cooling_rate in cooling_schedule.iter_mut() { + *cooling_rate = if *cooling_rate > 0 { *cooling_rate - 1 } else { 0 }; + } + } + + if current_temperature as u32 >= challenge.min_value { + Ok(Some(Solution { items: coloring })) + } else { + Ok(None) + } +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; \ No newline at end of file diff --git a/tig-algorithms/src/knapsack/classic_quadkp/commercial.rs b/tig-algorithms/src/knapsack/classic_quadkp/commercial.rs new file mode 100644 index 00000000..c8ff1e38 --- /dev/null +++ b/tig-algorithms/src/knapsack/classic_quadkp/commercial.rs @@ -0,0 +1,167 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Commercial License v1.0 (the "License"); you +may not use this file except in compliance with the License. You may obtain a copy +of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + +use anyhow::Result; +use tig_challenges::knapsack::{Challenge, Solution}; + +pub fn solve_challenge(challenge: &Challenge) -> Result> { + let vertex_count = challenge.weights.len(); + + let mut edge_costs: Vec<(usize, f32)> = (0..vertex_count) + .map(|flow_index| { + let total_flow = challenge.values[flow_index] as i32 + + challenge.interaction_values[flow_index].iter().sum::(); + let cost = total_flow as f32 / challenge.weights[flow_index] as f32; + (flow_index, cost) + }) + .collect(); + + edge_costs.sort_unstable_by(|a, b| b.1.partial_cmp(&a.1).unwrap()); + + let mut coloring = Vec::with_capacity(vertex_count); + let mut uncolored = Vec::with_capacity(vertex_count); + let mut current_entropy = 0; + let mut current_temperature = 0; + + for &(flow_index, _) in &edge_costs { + if current_entropy + challenge.weights[flow_index] <= challenge.max_weight { + current_entropy += challenge.weights[flow_index]; + current_temperature += challenge.values[flow_index] as i32; + + for &colored in &coloring { + current_temperature += challenge.interaction_values[flow_index][colored]; + } + coloring.push(flow_index); + } else { + uncolored.push(flow_index); + } + } + + let mut mutation_rates = vec![0; vertex_count]; + for flow_index in 0..vertex_count { + mutation_rates[flow_index] = challenge.values[flow_index] as i32; + for &colored in &coloring { + mutation_rates[flow_index] += challenge.interaction_values[flow_index][colored]; + } + } + + let max_generations = 100; + let mut cooling_schedule = vec![0; vertex_count]; + + for _ in 0..max_generations { + let mut best_mutation = 0; + let mut best_crossover = None; + + for uncolored_index in 0..uncolored.len() { + let mutant = uncolored[uncolored_index]; + if cooling_schedule[mutant] > 0 { + continue; + } + + unsafe { + let mutant_fitness = *mutation_rates.get_unchecked(mutant); + let min_entropy_reduction = *challenge.weights.get_unchecked(mutant) as i32 - (challenge.max_weight as i32 - current_entropy as i32); + + if mutant_fitness < 0 { + continue; + } + + for colored_index in 0..coloring.len() { + let gene_to_remove = *coloring.get_unchecked(colored_index); + if *cooling_schedule.get_unchecked(gene_to_remove) > 0 { + continue; + } + + if min_entropy_reduction > 0 { + let removed_entropy = *challenge.weights.get_unchecked(gene_to_remove) as i32; + if removed_entropy < min_entropy_reduction { + continue; + } + } + + let fitness_change = mutant_fitness - *mutation_rates.get_unchecked(gene_to_remove) + - *challenge.interaction_values.get_unchecked(mutant).get_unchecked(gene_to_remove); + + if fitness_change > best_mutation { + best_mutation = fitness_change; + best_crossover = Some((uncolored_index, colored_index)); + } + } + } + } + + if let Some((uncolored_index, colored_index)) = best_crossover { + let gene_to_add = uncolored[uncolored_index]; + let gene_to_remove = coloring[colored_index]; + + coloring.swap_remove(colored_index); + uncolored.swap_remove(uncolored_index); + coloring.push(gene_to_add); + uncolored.push(gene_to_remove); + + current_temperature += best_mutation; + current_entropy = current_entropy + challenge.weights[gene_to_add] - challenge.weights[gene_to_remove]; + + unsafe { + for flow_index in 0..vertex_count { + *mutation_rates.get_unchecked_mut(flow_index) += + challenge.interaction_values.get_unchecked(flow_index).get_unchecked(gene_to_add) - + challenge.interaction_values.get_unchecked(flow_index).get_unchecked(gene_to_remove); + } + } + + cooling_schedule[gene_to_add] = 3; + cooling_schedule[gene_to_remove] = 3; + } else { + break; + } + + if current_temperature as u32 >= challenge.min_value { + return Ok(Some(Solution { items: coloring })); + } + + for cooling_rate in cooling_schedule.iter_mut() { + *cooling_rate = if *cooling_rate > 0 { *cooling_rate - 1 } else { 0 }; + } + } + + if current_temperature as u32 >= challenge.min_value { + Ok(Some(Solution { items: coloring })) + } else { + Ok(None) + } +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; \ No newline at end of file diff --git a/tig-algorithms/src/knapsack/classic_quadkp/inbound.rs b/tig-algorithms/src/knapsack/classic_quadkp/inbound.rs new file mode 100644 index 00000000..314acfb2 --- /dev/null +++ b/tig-algorithms/src/knapsack/classic_quadkp/inbound.rs @@ -0,0 +1,167 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Inbound Game License v1.0 or (at your option) any later +version (the "License"); you may not use this file except in compliance with the +License. You may obtain a copy of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + +use anyhow::Result; +use tig_challenges::knapsack::{Challenge, Solution}; + +pub fn solve_challenge(challenge: &Challenge) -> Result> { + let vertex_count = challenge.weights.len(); + + let mut edge_costs: Vec<(usize, f32)> = (0..vertex_count) + .map(|flow_index| { + let total_flow = challenge.values[flow_index] as i32 + + challenge.interaction_values[flow_index].iter().sum::(); + let cost = total_flow as f32 / challenge.weights[flow_index] as f32; + (flow_index, cost) + }) + .collect(); + + edge_costs.sort_unstable_by(|a, b| b.1.partial_cmp(&a.1).unwrap()); + + let mut coloring = Vec::with_capacity(vertex_count); + let mut uncolored = Vec::with_capacity(vertex_count); + let mut current_entropy = 0; + let mut current_temperature = 0; + + for &(flow_index, _) in &edge_costs { + if current_entropy + challenge.weights[flow_index] <= challenge.max_weight { + current_entropy += challenge.weights[flow_index]; + current_temperature += challenge.values[flow_index] as i32; + + for &colored in &coloring { + current_temperature += challenge.interaction_values[flow_index][colored]; + } + coloring.push(flow_index); + } else { + uncolored.push(flow_index); + } + } + + let mut mutation_rates = vec![0; vertex_count]; + for flow_index in 0..vertex_count { + mutation_rates[flow_index] = challenge.values[flow_index] as i32; + for &colored in &coloring { + mutation_rates[flow_index] += challenge.interaction_values[flow_index][colored]; + } + } + + let max_generations = 100; + let mut cooling_schedule = vec![0; vertex_count]; + + for _ in 0..max_generations { + let mut best_mutation = 0; + let mut best_crossover = None; + + for uncolored_index in 0..uncolored.len() { + let mutant = uncolored[uncolored_index]; + if cooling_schedule[mutant] > 0 { + continue; + } + + unsafe { + let mutant_fitness = *mutation_rates.get_unchecked(mutant); + let min_entropy_reduction = *challenge.weights.get_unchecked(mutant) as i32 - (challenge.max_weight as i32 - current_entropy as i32); + + if mutant_fitness < 0 { + continue; + } + + for colored_index in 0..coloring.len() { + let gene_to_remove = *coloring.get_unchecked(colored_index); + if *cooling_schedule.get_unchecked(gene_to_remove) > 0 { + continue; + } + + if min_entropy_reduction > 0 { + let removed_entropy = *challenge.weights.get_unchecked(gene_to_remove) as i32; + if removed_entropy < min_entropy_reduction { + continue; + } + } + + let fitness_change = mutant_fitness - *mutation_rates.get_unchecked(gene_to_remove) + - *challenge.interaction_values.get_unchecked(mutant).get_unchecked(gene_to_remove); + + if fitness_change > best_mutation { + best_mutation = fitness_change; + best_crossover = Some((uncolored_index, colored_index)); + } + } + } + } + + if let Some((uncolored_index, colored_index)) = best_crossover { + let gene_to_add = uncolored[uncolored_index]; + let gene_to_remove = coloring[colored_index]; + + coloring.swap_remove(colored_index); + uncolored.swap_remove(uncolored_index); + coloring.push(gene_to_add); + uncolored.push(gene_to_remove); + + current_temperature += best_mutation; + current_entropy = current_entropy + challenge.weights[gene_to_add] - challenge.weights[gene_to_remove]; + + unsafe { + for flow_index in 0..vertex_count { + *mutation_rates.get_unchecked_mut(flow_index) += + challenge.interaction_values.get_unchecked(flow_index).get_unchecked(gene_to_add) - + challenge.interaction_values.get_unchecked(flow_index).get_unchecked(gene_to_remove); + } + } + + cooling_schedule[gene_to_add] = 3; + cooling_schedule[gene_to_remove] = 3; + } else { + break; + } + + if current_temperature as u32 >= challenge.min_value { + return Ok(Some(Solution { items: coloring })); + } + + for cooling_rate in cooling_schedule.iter_mut() { + *cooling_rate = if *cooling_rate > 0 { *cooling_rate - 1 } else { 0 }; + } + } + + if current_temperature as u32 >= challenge.min_value { + Ok(Some(Solution { items: coloring })) + } else { + Ok(None) + } +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; \ No newline at end of file diff --git a/tig-algorithms/src/knapsack/classic_quadkp/innovator_outbound.rs b/tig-algorithms/src/knapsack/classic_quadkp/innovator_outbound.rs new file mode 100644 index 00000000..3c9de50f --- /dev/null +++ b/tig-algorithms/src/knapsack/classic_quadkp/innovator_outbound.rs @@ -0,0 +1,167 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Innovator Outbound Game License v1.0 (the "License"); you +may not use this file except in compliance with the License. You may obtain a copy +of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + +use anyhow::Result; +use tig_challenges::knapsack::{Challenge, Solution}; + +pub fn solve_challenge(challenge: &Challenge) -> Result> { + let vertex_count = challenge.weights.len(); + + let mut edge_costs: Vec<(usize, f32)> = (0..vertex_count) + .map(|flow_index| { + let total_flow = challenge.values[flow_index] as i32 + + challenge.interaction_values[flow_index].iter().sum::(); + let cost = total_flow as f32 / challenge.weights[flow_index] as f32; + (flow_index, cost) + }) + .collect(); + + edge_costs.sort_unstable_by(|a, b| b.1.partial_cmp(&a.1).unwrap()); + + let mut coloring = Vec::with_capacity(vertex_count); + let mut uncolored = Vec::with_capacity(vertex_count); + let mut current_entropy = 0; + let mut current_temperature = 0; + + for &(flow_index, _) in &edge_costs { + if current_entropy + challenge.weights[flow_index] <= challenge.max_weight { + current_entropy += challenge.weights[flow_index]; + current_temperature += challenge.values[flow_index] as i32; + + for &colored in &coloring { + current_temperature += challenge.interaction_values[flow_index][colored]; + } + coloring.push(flow_index); + } else { + uncolored.push(flow_index); + } + } + + let mut mutation_rates = vec![0; vertex_count]; + for flow_index in 0..vertex_count { + mutation_rates[flow_index] = challenge.values[flow_index] as i32; + for &colored in &coloring { + mutation_rates[flow_index] += challenge.interaction_values[flow_index][colored]; + } + } + + let max_generations = 100; + let mut cooling_schedule = vec![0; vertex_count]; + + for _ in 0..max_generations { + let mut best_mutation = 0; + let mut best_crossover = None; + + for uncolored_index in 0..uncolored.len() { + let mutant = uncolored[uncolored_index]; + if cooling_schedule[mutant] > 0 { + continue; + } + + unsafe { + let mutant_fitness = *mutation_rates.get_unchecked(mutant); + let min_entropy_reduction = *challenge.weights.get_unchecked(mutant) as i32 - (challenge.max_weight as i32 - current_entropy as i32); + + if mutant_fitness < 0 { + continue; + } + + for colored_index in 0..coloring.len() { + let gene_to_remove = *coloring.get_unchecked(colored_index); + if *cooling_schedule.get_unchecked(gene_to_remove) > 0 { + continue; + } + + if min_entropy_reduction > 0 { + let removed_entropy = *challenge.weights.get_unchecked(gene_to_remove) as i32; + if removed_entropy < min_entropy_reduction { + continue; + } + } + + let fitness_change = mutant_fitness - *mutation_rates.get_unchecked(gene_to_remove) + - *challenge.interaction_values.get_unchecked(mutant).get_unchecked(gene_to_remove); + + if fitness_change > best_mutation { + best_mutation = fitness_change; + best_crossover = Some((uncolored_index, colored_index)); + } + } + } + } + + if let Some((uncolored_index, colored_index)) = best_crossover { + let gene_to_add = uncolored[uncolored_index]; + let gene_to_remove = coloring[colored_index]; + + coloring.swap_remove(colored_index); + uncolored.swap_remove(uncolored_index); + coloring.push(gene_to_add); + uncolored.push(gene_to_remove); + + current_temperature += best_mutation; + current_entropy = current_entropy + challenge.weights[gene_to_add] - challenge.weights[gene_to_remove]; + + unsafe { + for flow_index in 0..vertex_count { + *mutation_rates.get_unchecked_mut(flow_index) += + challenge.interaction_values.get_unchecked(flow_index).get_unchecked(gene_to_add) - + challenge.interaction_values.get_unchecked(flow_index).get_unchecked(gene_to_remove); + } + } + + cooling_schedule[gene_to_add] = 3; + cooling_schedule[gene_to_remove] = 3; + } else { + break; + } + + if current_temperature as u32 >= challenge.min_value { + return Ok(Some(Solution { items: coloring })); + } + + for cooling_rate in cooling_schedule.iter_mut() { + *cooling_rate = if *cooling_rate > 0 { *cooling_rate - 1 } else { 0 }; + } + } + + if current_temperature as u32 >= challenge.min_value { + Ok(Some(Solution { items: coloring })) + } else { + Ok(None) + } +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; \ No newline at end of file diff --git a/tig-algorithms/src/knapsack/classic_quadkp/mod.rs b/tig-algorithms/src/knapsack/classic_quadkp/mod.rs new file mode 100644 index 00000000..fcec9672 --- /dev/null +++ b/tig-algorithms/src/knapsack/classic_quadkp/mod.rs @@ -0,0 +1,4 @@ +mod benchmarker_outbound; +pub use benchmarker_outbound::solve_challenge; +#[cfg(feature = "cuda")] +pub use benchmarker_outbound::{cuda_solve_challenge, KERNEL}; \ No newline at end of file diff --git a/tig-algorithms/src/knapsack/classic_quadkp/open_data.rs b/tig-algorithms/src/knapsack/classic_quadkp/open_data.rs new file mode 100644 index 00000000..5c7bc344 --- /dev/null +++ b/tig-algorithms/src/knapsack/classic_quadkp/open_data.rs @@ -0,0 +1,167 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Open Data License v1.0 or (at your option) any later version +(the "License"); you may not use this file except in compliance with the License. +You may obtain a copy of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + +use anyhow::Result; +use tig_challenges::knapsack::{Challenge, Solution}; + +pub fn solve_challenge(challenge: &Challenge) -> Result> { + let vertex_count = challenge.weights.len(); + + let mut edge_costs: Vec<(usize, f32)> = (0..vertex_count) + .map(|flow_index| { + let total_flow = challenge.values[flow_index] as i32 + + challenge.interaction_values[flow_index].iter().sum::(); + let cost = total_flow as f32 / challenge.weights[flow_index] as f32; + (flow_index, cost) + }) + .collect(); + + edge_costs.sort_unstable_by(|a, b| b.1.partial_cmp(&a.1).unwrap()); + + let mut coloring = Vec::with_capacity(vertex_count); + let mut uncolored = Vec::with_capacity(vertex_count); + let mut current_entropy = 0; + let mut current_temperature = 0; + + for &(flow_index, _) in &edge_costs { + if current_entropy + challenge.weights[flow_index] <= challenge.max_weight { + current_entropy += challenge.weights[flow_index]; + current_temperature += challenge.values[flow_index] as i32; + + for &colored in &coloring { + current_temperature += challenge.interaction_values[flow_index][colored]; + } + coloring.push(flow_index); + } else { + uncolored.push(flow_index); + } + } + + let mut mutation_rates = vec![0; vertex_count]; + for flow_index in 0..vertex_count { + mutation_rates[flow_index] = challenge.values[flow_index] as i32; + for &colored in &coloring { + mutation_rates[flow_index] += challenge.interaction_values[flow_index][colored]; + } + } + + let max_generations = 100; + let mut cooling_schedule = vec![0; vertex_count]; + + for _ in 0..max_generations { + let mut best_mutation = 0; + let mut best_crossover = None; + + for uncolored_index in 0..uncolored.len() { + let mutant = uncolored[uncolored_index]; + if cooling_schedule[mutant] > 0 { + continue; + } + + unsafe { + let mutant_fitness = *mutation_rates.get_unchecked(mutant); + let min_entropy_reduction = *challenge.weights.get_unchecked(mutant) as i32 - (challenge.max_weight as i32 - current_entropy as i32); + + if mutant_fitness < 0 { + continue; + } + + for colored_index in 0..coloring.len() { + let gene_to_remove = *coloring.get_unchecked(colored_index); + if *cooling_schedule.get_unchecked(gene_to_remove) > 0 { + continue; + } + + if min_entropy_reduction > 0 { + let removed_entropy = *challenge.weights.get_unchecked(gene_to_remove) as i32; + if removed_entropy < min_entropy_reduction { + continue; + } + } + + let fitness_change = mutant_fitness - *mutation_rates.get_unchecked(gene_to_remove) + - *challenge.interaction_values.get_unchecked(mutant).get_unchecked(gene_to_remove); + + if fitness_change > best_mutation { + best_mutation = fitness_change; + best_crossover = Some((uncolored_index, colored_index)); + } + } + } + } + + if let Some((uncolored_index, colored_index)) = best_crossover { + let gene_to_add = uncolored[uncolored_index]; + let gene_to_remove = coloring[colored_index]; + + coloring.swap_remove(colored_index); + uncolored.swap_remove(uncolored_index); + coloring.push(gene_to_add); + uncolored.push(gene_to_remove); + + current_temperature += best_mutation; + current_entropy = current_entropy + challenge.weights[gene_to_add] - challenge.weights[gene_to_remove]; + + unsafe { + for flow_index in 0..vertex_count { + *mutation_rates.get_unchecked_mut(flow_index) += + challenge.interaction_values.get_unchecked(flow_index).get_unchecked(gene_to_add) - + challenge.interaction_values.get_unchecked(flow_index).get_unchecked(gene_to_remove); + } + } + + cooling_schedule[gene_to_add] = 3; + cooling_schedule[gene_to_remove] = 3; + } else { + break; + } + + if current_temperature as u32 >= challenge.min_value { + return Ok(Some(Solution { items: coloring })); + } + + for cooling_rate in cooling_schedule.iter_mut() { + *cooling_rate = if *cooling_rate > 0 { *cooling_rate - 1 } else { 0 }; + } + } + + if current_temperature as u32 >= challenge.min_value { + Ok(Some(Solution { items: coloring })) + } else { + Ok(None) + } +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; \ No newline at end of file diff --git a/tig-algorithms/src/knapsack/mod.rs b/tig-algorithms/src/knapsack/mod.rs index edf17b55..b551497a 100644 --- a/tig-algorithms/src/knapsack/mod.rs +++ b/tig-algorithms/src/knapsack/mod.rs @@ -101,7 +101,8 @@ pub use knapheudp as c003_a019; // c003_a050 -// c003_a051 +pub mod classic_quadkp; +pub use classic_quadkp as c003_a051; // c003_a052 diff --git a/tig-algorithms/src/satisfiability/mod.rs b/tig-algorithms/src/satisfiability/mod.rs index a19e4e14..cf7dad6c 100644 --- a/tig-algorithms/src/satisfiability/mod.rs +++ b/tig-algorithms/src/satisfiability/mod.rs @@ -86,7 +86,8 @@ pub use sat_global as c001_a034; // c001_a040 -// c001_a041 +pub mod sat_global_opt; +pub use sat_global_opt as c001_a041; // c001_a042 diff --git a/tig-algorithms/src/satisfiability/sat_global_opt/benchmarker_outbound.rs b/tig-algorithms/src/satisfiability/sat_global_opt/benchmarker_outbound.rs new file mode 100644 index 00000000..bf1e1668 --- /dev/null +++ b/tig-algorithms/src/satisfiability/sat_global_opt/benchmarker_outbound.rs @@ -0,0 +1,292 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Benchmarker Outbound Game License v1.0 (the "License"); you +may not use this file except in compliance with the License. You may obtain a copy +of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + +use rand::{rngs::{SmallRng, StdRng}, Rng, SeedableRng}; +use std::collections::HashMap; +use tig_challenges::satisfiability::*; + +pub fn solve_challenge(challenge: &Challenge) -> anyhow::Result> { + let mut rng = SmallRng::seed_from_u64(u64::from_le_bytes(challenge.seed[..8].try_into().unwrap()) as u64); + + let mut p_single = vec![false; challenge.difficulty.num_variables]; + let mut n_single = vec![false; challenge.difficulty.num_variables]; + + let mut clauses_ = challenge.clauses.clone(); + let mut clauses: Vec> = Vec::with_capacity(clauses_.len()); + + let mut rounds = 0; + + let mut dead = false; + + while !(dead) { + let mut done = true; + for c in &clauses_ { + let mut c_: Vec = Vec::with_capacity(c.len()); // Preallocate with capacity + let mut skip = false; + for (i, l) in c.iter().enumerate() { + if (p_single[(l.abs() - 1) as usize] && *l > 0) + || (n_single[(l.abs() - 1) as usize] && *l < 0) + || c[(i + 1)..].contains(&-l) + { + skip = true; + break; + } else if p_single[(l.abs() - 1) as usize] + || n_single[(l.abs() - 1) as usize] + || c[(i + 1)..].contains(&l) + { + done = false; + continue; + } else { + c_.push(*l); + } + } + if skip { + done = false; + continue; + }; + match c_[..] { + [l] => { + done = false; + if l > 0 { + if n_single[(l.abs() - 1) as usize] { + dead = true; + break; + } else { + p_single[(l.abs() - 1) as usize] = true; + } + } else { + if p_single[(l.abs() - 1) as usize] { + dead = true; + break; + } else { + n_single[(l.abs() - 1) as usize] = true; + } + } + } + [] => { + dead = true; + break; + } + _ => { + clauses.push(c_); + } + } + } + if done { + break; + } else { + clauses_ = clauses; + clauses = Vec::with_capacity(clauses_.len()); + } + } + + if dead { + return Ok(None); + } + + let num_variables = challenge.difficulty.num_variables; + let num_clauses = clauses.len(); + + let mut p_clauses: Vec> = vec![Vec::new(); num_variables]; + let mut n_clauses: Vec> = vec![Vec::new(); num_variables]; + + // Preallocate capacity for p_clauses and n_clauses + for c in &clauses { + for &l in c { + let var = (l.abs() - 1) as usize; + if l > 0 { + if p_clauses[var].capacity() == 0 { + p_clauses[var] = Vec::with_capacity(clauses.len() / num_variables + 1); + } + } else { + if n_clauses[var].capacity() == 0 { + n_clauses[var] = Vec::with_capacity(clauses.len() / num_variables + 1); + } + } + } + } + + for (i, &ref c) in clauses.iter().enumerate() { + for &l in c { + let var = (l.abs() - 1) as usize; + if l > 0 { + p_clauses[var].push(i); + } else { + n_clauses[var].push(i); + } + } + } + + let mut variables = vec![false; num_variables]; + for v in 0..num_variables { + let num_p = p_clauses[v].len(); + let num_n = n_clauses[v].len(); + + let nad = 1.28; + let mut vad = nad + 1.0; + if num_n > 0 { + vad = num_p as f32 / num_n as f32; + } + + if vad <= nad { + variables[v] = false; + } else { + let prob = num_p as f64 / (num_p + num_n).max(1) as f64; + variables[v] = rng.gen_bool(prob) + } + } + + let mut num_good_so_far: Vec = vec![0; num_clauses]; + for (i, &ref c) in clauses.iter().enumerate() { + for &l in c { + let var = (l.abs() - 1) as usize; + if l > 0 && variables[var] { + num_good_so_far[i] += 1 + } else if l < 0 && !variables[var] { + num_good_so_far[i] += 1 + } + } + } + + + let mut residual_ = Vec::with_capacity(num_clauses); + let mut residual_indices = vec![None; num_clauses]; + + for (i, &num_good) in num_good_so_far.iter().enumerate() { + if num_good == 0 { + residual_.push(i); + residual_indices[i] = Some(residual_.len() - 1); + } + } + + let clauses_ratio = challenge.difficulty.clauses_to_variables_percent as f64; + let num_vars = challenge.difficulty.num_variables as f64; + let max_fuel = 2000000000.0; + let base_fuel = (2000.0 + 40.0 * clauses_ratio) * num_vars; + let flip_fuel = 350.0 + 0.9 * clauses_ratio; + let max_num_rounds = ((max_fuel - base_fuel) / flip_fuel) as usize; + loop { + if !residual_.is_empty() { + + let rand_val = rng.gen::(); + + let i = residual_[rand_val % residual_.len()]; + let mut min_sad = clauses.len(); + let mut v_min_sad = usize::MAX; + let c = &mut clauses[i]; + + if c.len() > 1 { + let random_index = rand_val % c.len(); + c.swap(0, random_index); + } + for &l in c.iter() { + let abs_l = l.abs() as usize - 1; + let clauses_to_check = if variables[abs_l] { &p_clauses[abs_l] } else { &n_clauses[abs_l] }; + + let mut sad = 0; + for &c in clauses_to_check { + if num_good_so_far[c] == 1 { + sad += 1; + } + } + + if sad < min_sad { + min_sad = sad; + v_min_sad = abs_l; + } + } + + let v = if min_sad == 0 { + v_min_sad + } else if rng.gen_bool(0.5) { + c[0].abs() as usize - 1 + } else { + v_min_sad + }; + + if variables[v] { + for &c in &n_clauses[v] { + num_good_so_far[c] += 1; + if num_good_so_far[c] == 1 { + let i = residual_indices[c].take().unwrap(); + let last = residual_.pop().unwrap(); + if i < residual_.len() { + residual_[i] = last; + residual_indices[last] = Some(i); + } + } + } + for &c in &p_clauses[v] { + if num_good_so_far[c] == 1 { + residual_.push(c); + residual_indices[c] = Some(residual_.len() - 1); + } + num_good_so_far[c] -= 1; + } + } else { + for &c in &n_clauses[v] { + if num_good_so_far[c] == 1 { + residual_.push(c); + residual_indices[c] = Some(residual_.len() - 1); + } + num_good_so_far[c] -= 1; + } + + for &c in &p_clauses[v] { + num_good_so_far[c] += 1; + if num_good_so_far[c] == 1 { + let i = residual_indices[c].take().unwrap(); + let last = residual_.pop().unwrap(); + if i < residual_.len() { + residual_[i] = last; + residual_indices[last] = Some(i); + } + } + } + } + + variables[v] = !variables[v]; + } else { + break; + } + rounds += 1; + if rounds >= max_num_rounds { + return Ok(None); + } + } + return Ok(Some(Solution { variables })); +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; diff --git a/tig-algorithms/src/satisfiability/sat_global_opt/commercial.rs b/tig-algorithms/src/satisfiability/sat_global_opt/commercial.rs new file mode 100644 index 00000000..1a41da56 --- /dev/null +++ b/tig-algorithms/src/satisfiability/sat_global_opt/commercial.rs @@ -0,0 +1,292 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Commercial License v1.0 (the "License"); you +may not use this file except in compliance with the License. You may obtain a copy +of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + +use rand::{rngs::{SmallRng, StdRng}, Rng, SeedableRng}; +use std::collections::HashMap; +use tig_challenges::satisfiability::*; + +pub fn solve_challenge(challenge: &Challenge) -> anyhow::Result> { + let mut rng = SmallRng::seed_from_u64(u64::from_le_bytes(challenge.seed[..8].try_into().unwrap()) as u64); + + let mut p_single = vec![false; challenge.difficulty.num_variables]; + let mut n_single = vec![false; challenge.difficulty.num_variables]; + + let mut clauses_ = challenge.clauses.clone(); + let mut clauses: Vec> = Vec::with_capacity(clauses_.len()); + + let mut rounds = 0; + + let mut dead = false; + + while !(dead) { + let mut done = true; + for c in &clauses_ { + let mut c_: Vec = Vec::with_capacity(c.len()); // Preallocate with capacity + let mut skip = false; + for (i, l) in c.iter().enumerate() { + if (p_single[(l.abs() - 1) as usize] && *l > 0) + || (n_single[(l.abs() - 1) as usize] && *l < 0) + || c[(i + 1)..].contains(&-l) + { + skip = true; + break; + } else if p_single[(l.abs() - 1) as usize] + || n_single[(l.abs() - 1) as usize] + || c[(i + 1)..].contains(&l) + { + done = false; + continue; + } else { + c_.push(*l); + } + } + if skip { + done = false; + continue; + }; + match c_[..] { + [l] => { + done = false; + if l > 0 { + if n_single[(l.abs() - 1) as usize] { + dead = true; + break; + } else { + p_single[(l.abs() - 1) as usize] = true; + } + } else { + if p_single[(l.abs() - 1) as usize] { + dead = true; + break; + } else { + n_single[(l.abs() - 1) as usize] = true; + } + } + } + [] => { + dead = true; + break; + } + _ => { + clauses.push(c_); + } + } + } + if done { + break; + } else { + clauses_ = clauses; + clauses = Vec::with_capacity(clauses_.len()); + } + } + + if dead { + return Ok(None); + } + + let num_variables = challenge.difficulty.num_variables; + let num_clauses = clauses.len(); + + let mut p_clauses: Vec> = vec![Vec::new(); num_variables]; + let mut n_clauses: Vec> = vec![Vec::new(); num_variables]; + + // Preallocate capacity for p_clauses and n_clauses + for c in &clauses { + for &l in c { + let var = (l.abs() - 1) as usize; + if l > 0 { + if p_clauses[var].capacity() == 0 { + p_clauses[var] = Vec::with_capacity(clauses.len() / num_variables + 1); + } + } else { + if n_clauses[var].capacity() == 0 { + n_clauses[var] = Vec::with_capacity(clauses.len() / num_variables + 1); + } + } + } + } + + for (i, &ref c) in clauses.iter().enumerate() { + for &l in c { + let var = (l.abs() - 1) as usize; + if l > 0 { + p_clauses[var].push(i); + } else { + n_clauses[var].push(i); + } + } + } + + let mut variables = vec![false; num_variables]; + for v in 0..num_variables { + let num_p = p_clauses[v].len(); + let num_n = n_clauses[v].len(); + + let nad = 1.28; + let mut vad = nad + 1.0; + if num_n > 0 { + vad = num_p as f32 / num_n as f32; + } + + if vad <= nad { + variables[v] = false; + } else { + let prob = num_p as f64 / (num_p + num_n).max(1) as f64; + variables[v] = rng.gen_bool(prob) + } + } + + let mut num_good_so_far: Vec = vec![0; num_clauses]; + for (i, &ref c) in clauses.iter().enumerate() { + for &l in c { + let var = (l.abs() - 1) as usize; + if l > 0 && variables[var] { + num_good_so_far[i] += 1 + } else if l < 0 && !variables[var] { + num_good_so_far[i] += 1 + } + } + } + + + let mut residual_ = Vec::with_capacity(num_clauses); + let mut residual_indices = vec![None; num_clauses]; + + for (i, &num_good) in num_good_so_far.iter().enumerate() { + if num_good == 0 { + residual_.push(i); + residual_indices[i] = Some(residual_.len() - 1); + } + } + + let clauses_ratio = challenge.difficulty.clauses_to_variables_percent as f64; + let num_vars = challenge.difficulty.num_variables as f64; + let max_fuel = 2000000000.0; + let base_fuel = (2000.0 + 40.0 * clauses_ratio) * num_vars; + let flip_fuel = 350.0 + 0.9 * clauses_ratio; + let max_num_rounds = ((max_fuel - base_fuel) / flip_fuel) as usize; + loop { + if !residual_.is_empty() { + + let rand_val = rng.gen::(); + + let i = residual_[rand_val % residual_.len()]; + let mut min_sad = clauses.len(); + let mut v_min_sad = usize::MAX; + let c = &mut clauses[i]; + + if c.len() > 1 { + let random_index = rand_val % c.len(); + c.swap(0, random_index); + } + for &l in c.iter() { + let abs_l = l.abs() as usize - 1; + let clauses_to_check = if variables[abs_l] { &p_clauses[abs_l] } else { &n_clauses[abs_l] }; + + let mut sad = 0; + for &c in clauses_to_check { + if num_good_so_far[c] == 1 { + sad += 1; + } + } + + if sad < min_sad { + min_sad = sad; + v_min_sad = abs_l; + } + } + + let v = if min_sad == 0 { + v_min_sad + } else if rng.gen_bool(0.5) { + c[0].abs() as usize - 1 + } else { + v_min_sad + }; + + if variables[v] { + for &c in &n_clauses[v] { + num_good_so_far[c] += 1; + if num_good_so_far[c] == 1 { + let i = residual_indices[c].take().unwrap(); + let last = residual_.pop().unwrap(); + if i < residual_.len() { + residual_[i] = last; + residual_indices[last] = Some(i); + } + } + } + for &c in &p_clauses[v] { + if num_good_so_far[c] == 1 { + residual_.push(c); + residual_indices[c] = Some(residual_.len() - 1); + } + num_good_so_far[c] -= 1; + } + } else { + for &c in &n_clauses[v] { + if num_good_so_far[c] == 1 { + residual_.push(c); + residual_indices[c] = Some(residual_.len() - 1); + } + num_good_so_far[c] -= 1; + } + + for &c in &p_clauses[v] { + num_good_so_far[c] += 1; + if num_good_so_far[c] == 1 { + let i = residual_indices[c].take().unwrap(); + let last = residual_.pop().unwrap(); + if i < residual_.len() { + residual_[i] = last; + residual_indices[last] = Some(i); + } + } + } + } + + variables[v] = !variables[v]; + } else { + break; + } + rounds += 1; + if rounds >= max_num_rounds { + return Ok(None); + } + } + return Ok(Some(Solution { variables })); +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; diff --git a/tig-algorithms/src/satisfiability/sat_global_opt/inbound.rs b/tig-algorithms/src/satisfiability/sat_global_opt/inbound.rs new file mode 100644 index 00000000..b39f496b --- /dev/null +++ b/tig-algorithms/src/satisfiability/sat_global_opt/inbound.rs @@ -0,0 +1,292 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Inbound Game License v1.0 or (at your option) any later +version (the "License"); you may not use this file except in compliance with the +License. You may obtain a copy of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + +use rand::{rngs::{SmallRng, StdRng}, Rng, SeedableRng}; +use std::collections::HashMap; +use tig_challenges::satisfiability::*; + +pub fn solve_challenge(challenge: &Challenge) -> anyhow::Result> { + let mut rng = SmallRng::seed_from_u64(u64::from_le_bytes(challenge.seed[..8].try_into().unwrap()) as u64); + + let mut p_single = vec![false; challenge.difficulty.num_variables]; + let mut n_single = vec![false; challenge.difficulty.num_variables]; + + let mut clauses_ = challenge.clauses.clone(); + let mut clauses: Vec> = Vec::with_capacity(clauses_.len()); + + let mut rounds = 0; + + let mut dead = false; + + while !(dead) { + let mut done = true; + for c in &clauses_ { + let mut c_: Vec = Vec::with_capacity(c.len()); // Preallocate with capacity + let mut skip = false; + for (i, l) in c.iter().enumerate() { + if (p_single[(l.abs() - 1) as usize] && *l > 0) + || (n_single[(l.abs() - 1) as usize] && *l < 0) + || c[(i + 1)..].contains(&-l) + { + skip = true; + break; + } else if p_single[(l.abs() - 1) as usize] + || n_single[(l.abs() - 1) as usize] + || c[(i + 1)..].contains(&l) + { + done = false; + continue; + } else { + c_.push(*l); + } + } + if skip { + done = false; + continue; + }; + match c_[..] { + [l] => { + done = false; + if l > 0 { + if n_single[(l.abs() - 1) as usize] { + dead = true; + break; + } else { + p_single[(l.abs() - 1) as usize] = true; + } + } else { + if p_single[(l.abs() - 1) as usize] { + dead = true; + break; + } else { + n_single[(l.abs() - 1) as usize] = true; + } + } + } + [] => { + dead = true; + break; + } + _ => { + clauses.push(c_); + } + } + } + if done { + break; + } else { + clauses_ = clauses; + clauses = Vec::with_capacity(clauses_.len()); + } + } + + if dead { + return Ok(None); + } + + let num_variables = challenge.difficulty.num_variables; + let num_clauses = clauses.len(); + + let mut p_clauses: Vec> = vec![Vec::new(); num_variables]; + let mut n_clauses: Vec> = vec![Vec::new(); num_variables]; + + // Preallocate capacity for p_clauses and n_clauses + for c in &clauses { + for &l in c { + let var = (l.abs() - 1) as usize; + if l > 0 { + if p_clauses[var].capacity() == 0 { + p_clauses[var] = Vec::with_capacity(clauses.len() / num_variables + 1); + } + } else { + if n_clauses[var].capacity() == 0 { + n_clauses[var] = Vec::with_capacity(clauses.len() / num_variables + 1); + } + } + } + } + + for (i, &ref c) in clauses.iter().enumerate() { + for &l in c { + let var = (l.abs() - 1) as usize; + if l > 0 { + p_clauses[var].push(i); + } else { + n_clauses[var].push(i); + } + } + } + + let mut variables = vec![false; num_variables]; + for v in 0..num_variables { + let num_p = p_clauses[v].len(); + let num_n = n_clauses[v].len(); + + let nad = 1.28; + let mut vad = nad + 1.0; + if num_n > 0 { + vad = num_p as f32 / num_n as f32; + } + + if vad <= nad { + variables[v] = false; + } else { + let prob = num_p as f64 / (num_p + num_n).max(1) as f64; + variables[v] = rng.gen_bool(prob) + } + } + + let mut num_good_so_far: Vec = vec![0; num_clauses]; + for (i, &ref c) in clauses.iter().enumerate() { + for &l in c { + let var = (l.abs() - 1) as usize; + if l > 0 && variables[var] { + num_good_so_far[i] += 1 + } else if l < 0 && !variables[var] { + num_good_so_far[i] += 1 + } + } + } + + + let mut residual_ = Vec::with_capacity(num_clauses); + let mut residual_indices = vec![None; num_clauses]; + + for (i, &num_good) in num_good_so_far.iter().enumerate() { + if num_good == 0 { + residual_.push(i); + residual_indices[i] = Some(residual_.len() - 1); + } + } + + let clauses_ratio = challenge.difficulty.clauses_to_variables_percent as f64; + let num_vars = challenge.difficulty.num_variables as f64; + let max_fuel = 2000000000.0; + let base_fuel = (2000.0 + 40.0 * clauses_ratio) * num_vars; + let flip_fuel = 350.0 + 0.9 * clauses_ratio; + let max_num_rounds = ((max_fuel - base_fuel) / flip_fuel) as usize; + loop { + if !residual_.is_empty() { + + let rand_val = rng.gen::(); + + let i = residual_[rand_val % residual_.len()]; + let mut min_sad = clauses.len(); + let mut v_min_sad = usize::MAX; + let c = &mut clauses[i]; + + if c.len() > 1 { + let random_index = rand_val % c.len(); + c.swap(0, random_index); + } + for &l in c.iter() { + let abs_l = l.abs() as usize - 1; + let clauses_to_check = if variables[abs_l] { &p_clauses[abs_l] } else { &n_clauses[abs_l] }; + + let mut sad = 0; + for &c in clauses_to_check { + if num_good_so_far[c] == 1 { + sad += 1; + } + } + + if sad < min_sad { + min_sad = sad; + v_min_sad = abs_l; + } + } + + let v = if min_sad == 0 { + v_min_sad + } else if rng.gen_bool(0.5) { + c[0].abs() as usize - 1 + } else { + v_min_sad + }; + + if variables[v] { + for &c in &n_clauses[v] { + num_good_so_far[c] += 1; + if num_good_so_far[c] == 1 { + let i = residual_indices[c].take().unwrap(); + let last = residual_.pop().unwrap(); + if i < residual_.len() { + residual_[i] = last; + residual_indices[last] = Some(i); + } + } + } + for &c in &p_clauses[v] { + if num_good_so_far[c] == 1 { + residual_.push(c); + residual_indices[c] = Some(residual_.len() - 1); + } + num_good_so_far[c] -= 1; + } + } else { + for &c in &n_clauses[v] { + if num_good_so_far[c] == 1 { + residual_.push(c); + residual_indices[c] = Some(residual_.len() - 1); + } + num_good_so_far[c] -= 1; + } + + for &c in &p_clauses[v] { + num_good_so_far[c] += 1; + if num_good_so_far[c] == 1 { + let i = residual_indices[c].take().unwrap(); + let last = residual_.pop().unwrap(); + if i < residual_.len() { + residual_[i] = last; + residual_indices[last] = Some(i); + } + } + } + } + + variables[v] = !variables[v]; + } else { + break; + } + rounds += 1; + if rounds >= max_num_rounds { + return Ok(None); + } + } + return Ok(Some(Solution { variables })); +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; diff --git a/tig-algorithms/src/satisfiability/sat_global_opt/innovator_outbound.rs b/tig-algorithms/src/satisfiability/sat_global_opt/innovator_outbound.rs new file mode 100644 index 00000000..18bbd685 --- /dev/null +++ b/tig-algorithms/src/satisfiability/sat_global_opt/innovator_outbound.rs @@ -0,0 +1,292 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Innovator Outbound Game License v1.0 (the "License"); you +may not use this file except in compliance with the License. You may obtain a copy +of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + +use rand::{rngs::{SmallRng, StdRng}, Rng, SeedableRng}; +use std::collections::HashMap; +use tig_challenges::satisfiability::*; + +pub fn solve_challenge(challenge: &Challenge) -> anyhow::Result> { + let mut rng = SmallRng::seed_from_u64(u64::from_le_bytes(challenge.seed[..8].try_into().unwrap()) as u64); + + let mut p_single = vec![false; challenge.difficulty.num_variables]; + let mut n_single = vec![false; challenge.difficulty.num_variables]; + + let mut clauses_ = challenge.clauses.clone(); + let mut clauses: Vec> = Vec::with_capacity(clauses_.len()); + + let mut rounds = 0; + + let mut dead = false; + + while !(dead) { + let mut done = true; + for c in &clauses_ { + let mut c_: Vec = Vec::with_capacity(c.len()); // Preallocate with capacity + let mut skip = false; + for (i, l) in c.iter().enumerate() { + if (p_single[(l.abs() - 1) as usize] && *l > 0) + || (n_single[(l.abs() - 1) as usize] && *l < 0) + || c[(i + 1)..].contains(&-l) + { + skip = true; + break; + } else if p_single[(l.abs() - 1) as usize] + || n_single[(l.abs() - 1) as usize] + || c[(i + 1)..].contains(&l) + { + done = false; + continue; + } else { + c_.push(*l); + } + } + if skip { + done = false; + continue; + }; + match c_[..] { + [l] => { + done = false; + if l > 0 { + if n_single[(l.abs() - 1) as usize] { + dead = true; + break; + } else { + p_single[(l.abs() - 1) as usize] = true; + } + } else { + if p_single[(l.abs() - 1) as usize] { + dead = true; + break; + } else { + n_single[(l.abs() - 1) as usize] = true; + } + } + } + [] => { + dead = true; + break; + } + _ => { + clauses.push(c_); + } + } + } + if done { + break; + } else { + clauses_ = clauses; + clauses = Vec::with_capacity(clauses_.len()); + } + } + + if dead { + return Ok(None); + } + + let num_variables = challenge.difficulty.num_variables; + let num_clauses = clauses.len(); + + let mut p_clauses: Vec> = vec![Vec::new(); num_variables]; + let mut n_clauses: Vec> = vec![Vec::new(); num_variables]; + + // Preallocate capacity for p_clauses and n_clauses + for c in &clauses { + for &l in c { + let var = (l.abs() - 1) as usize; + if l > 0 { + if p_clauses[var].capacity() == 0 { + p_clauses[var] = Vec::with_capacity(clauses.len() / num_variables + 1); + } + } else { + if n_clauses[var].capacity() == 0 { + n_clauses[var] = Vec::with_capacity(clauses.len() / num_variables + 1); + } + } + } + } + + for (i, &ref c) in clauses.iter().enumerate() { + for &l in c { + let var = (l.abs() - 1) as usize; + if l > 0 { + p_clauses[var].push(i); + } else { + n_clauses[var].push(i); + } + } + } + + let mut variables = vec![false; num_variables]; + for v in 0..num_variables { + let num_p = p_clauses[v].len(); + let num_n = n_clauses[v].len(); + + let nad = 1.28; + let mut vad = nad + 1.0; + if num_n > 0 { + vad = num_p as f32 / num_n as f32; + } + + if vad <= nad { + variables[v] = false; + } else { + let prob = num_p as f64 / (num_p + num_n).max(1) as f64; + variables[v] = rng.gen_bool(prob) + } + } + + let mut num_good_so_far: Vec = vec![0; num_clauses]; + for (i, &ref c) in clauses.iter().enumerate() { + for &l in c { + let var = (l.abs() - 1) as usize; + if l > 0 && variables[var] { + num_good_so_far[i] += 1 + } else if l < 0 && !variables[var] { + num_good_so_far[i] += 1 + } + } + } + + + let mut residual_ = Vec::with_capacity(num_clauses); + let mut residual_indices = vec![None; num_clauses]; + + for (i, &num_good) in num_good_so_far.iter().enumerate() { + if num_good == 0 { + residual_.push(i); + residual_indices[i] = Some(residual_.len() - 1); + } + } + + let clauses_ratio = challenge.difficulty.clauses_to_variables_percent as f64; + let num_vars = challenge.difficulty.num_variables as f64; + let max_fuel = 2000000000.0; + let base_fuel = (2000.0 + 40.0 * clauses_ratio) * num_vars; + let flip_fuel = 350.0 + 0.9 * clauses_ratio; + let max_num_rounds = ((max_fuel - base_fuel) / flip_fuel) as usize; + loop { + if !residual_.is_empty() { + + let rand_val = rng.gen::(); + + let i = residual_[rand_val % residual_.len()]; + let mut min_sad = clauses.len(); + let mut v_min_sad = usize::MAX; + let c = &mut clauses[i]; + + if c.len() > 1 { + let random_index = rand_val % c.len(); + c.swap(0, random_index); + } + for &l in c.iter() { + let abs_l = l.abs() as usize - 1; + let clauses_to_check = if variables[abs_l] { &p_clauses[abs_l] } else { &n_clauses[abs_l] }; + + let mut sad = 0; + for &c in clauses_to_check { + if num_good_so_far[c] == 1 { + sad += 1; + } + } + + if sad < min_sad { + min_sad = sad; + v_min_sad = abs_l; + } + } + + let v = if min_sad == 0 { + v_min_sad + } else if rng.gen_bool(0.5) { + c[0].abs() as usize - 1 + } else { + v_min_sad + }; + + if variables[v] { + for &c in &n_clauses[v] { + num_good_so_far[c] += 1; + if num_good_so_far[c] == 1 { + let i = residual_indices[c].take().unwrap(); + let last = residual_.pop().unwrap(); + if i < residual_.len() { + residual_[i] = last; + residual_indices[last] = Some(i); + } + } + } + for &c in &p_clauses[v] { + if num_good_so_far[c] == 1 { + residual_.push(c); + residual_indices[c] = Some(residual_.len() - 1); + } + num_good_so_far[c] -= 1; + } + } else { + for &c in &n_clauses[v] { + if num_good_so_far[c] == 1 { + residual_.push(c); + residual_indices[c] = Some(residual_.len() - 1); + } + num_good_so_far[c] -= 1; + } + + for &c in &p_clauses[v] { + num_good_so_far[c] += 1; + if num_good_so_far[c] == 1 { + let i = residual_indices[c].take().unwrap(); + let last = residual_.pop().unwrap(); + if i < residual_.len() { + residual_[i] = last; + residual_indices[last] = Some(i); + } + } + } + } + + variables[v] = !variables[v]; + } else { + break; + } + rounds += 1; + if rounds >= max_num_rounds { + return Ok(None); + } + } + return Ok(Some(Solution { variables })); +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; diff --git a/tig-algorithms/src/satisfiability/sat_global_opt/mod.rs b/tig-algorithms/src/satisfiability/sat_global_opt/mod.rs new file mode 100644 index 00000000..fcec9672 --- /dev/null +++ b/tig-algorithms/src/satisfiability/sat_global_opt/mod.rs @@ -0,0 +1,4 @@ +mod benchmarker_outbound; +pub use benchmarker_outbound::solve_challenge; +#[cfg(feature = "cuda")] +pub use benchmarker_outbound::{cuda_solve_challenge, KERNEL}; \ No newline at end of file diff --git a/tig-algorithms/src/satisfiability/sat_global_opt/open_data.rs b/tig-algorithms/src/satisfiability/sat_global_opt/open_data.rs new file mode 100644 index 00000000..eb4a1c0e --- /dev/null +++ b/tig-algorithms/src/satisfiability/sat_global_opt/open_data.rs @@ -0,0 +1,292 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Open Data License v1.0 or (at your option) any later version +(the "License"); you may not use this file except in compliance with the License. +You may obtain a copy of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + +use rand::{rngs::{SmallRng, StdRng}, Rng, SeedableRng}; +use std::collections::HashMap; +use tig_challenges::satisfiability::*; + +pub fn solve_challenge(challenge: &Challenge) -> anyhow::Result> { + let mut rng = SmallRng::seed_from_u64(u64::from_le_bytes(challenge.seed[..8].try_into().unwrap()) as u64); + + let mut p_single = vec![false; challenge.difficulty.num_variables]; + let mut n_single = vec![false; challenge.difficulty.num_variables]; + + let mut clauses_ = challenge.clauses.clone(); + let mut clauses: Vec> = Vec::with_capacity(clauses_.len()); + + let mut rounds = 0; + + let mut dead = false; + + while !(dead) { + let mut done = true; + for c in &clauses_ { + let mut c_: Vec = Vec::with_capacity(c.len()); // Preallocate with capacity + let mut skip = false; + for (i, l) in c.iter().enumerate() { + if (p_single[(l.abs() - 1) as usize] && *l > 0) + || (n_single[(l.abs() - 1) as usize] && *l < 0) + || c[(i + 1)..].contains(&-l) + { + skip = true; + break; + } else if p_single[(l.abs() - 1) as usize] + || n_single[(l.abs() - 1) as usize] + || c[(i + 1)..].contains(&l) + { + done = false; + continue; + } else { + c_.push(*l); + } + } + if skip { + done = false; + continue; + }; + match c_[..] { + [l] => { + done = false; + if l > 0 { + if n_single[(l.abs() - 1) as usize] { + dead = true; + break; + } else { + p_single[(l.abs() - 1) as usize] = true; + } + } else { + if p_single[(l.abs() - 1) as usize] { + dead = true; + break; + } else { + n_single[(l.abs() - 1) as usize] = true; + } + } + } + [] => { + dead = true; + break; + } + _ => { + clauses.push(c_); + } + } + } + if done { + break; + } else { + clauses_ = clauses; + clauses = Vec::with_capacity(clauses_.len()); + } + } + + if dead { + return Ok(None); + } + + let num_variables = challenge.difficulty.num_variables; + let num_clauses = clauses.len(); + + let mut p_clauses: Vec> = vec![Vec::new(); num_variables]; + let mut n_clauses: Vec> = vec![Vec::new(); num_variables]; + + // Preallocate capacity for p_clauses and n_clauses + for c in &clauses { + for &l in c { + let var = (l.abs() - 1) as usize; + if l > 0 { + if p_clauses[var].capacity() == 0 { + p_clauses[var] = Vec::with_capacity(clauses.len() / num_variables + 1); + } + } else { + if n_clauses[var].capacity() == 0 { + n_clauses[var] = Vec::with_capacity(clauses.len() / num_variables + 1); + } + } + } + } + + for (i, &ref c) in clauses.iter().enumerate() { + for &l in c { + let var = (l.abs() - 1) as usize; + if l > 0 { + p_clauses[var].push(i); + } else { + n_clauses[var].push(i); + } + } + } + + let mut variables = vec![false; num_variables]; + for v in 0..num_variables { + let num_p = p_clauses[v].len(); + let num_n = n_clauses[v].len(); + + let nad = 1.28; + let mut vad = nad + 1.0; + if num_n > 0 { + vad = num_p as f32 / num_n as f32; + } + + if vad <= nad { + variables[v] = false; + } else { + let prob = num_p as f64 / (num_p + num_n).max(1) as f64; + variables[v] = rng.gen_bool(prob) + } + } + + let mut num_good_so_far: Vec = vec![0; num_clauses]; + for (i, &ref c) in clauses.iter().enumerate() { + for &l in c { + let var = (l.abs() - 1) as usize; + if l > 0 && variables[var] { + num_good_so_far[i] += 1 + } else if l < 0 && !variables[var] { + num_good_so_far[i] += 1 + } + } + } + + + let mut residual_ = Vec::with_capacity(num_clauses); + let mut residual_indices = vec![None; num_clauses]; + + for (i, &num_good) in num_good_so_far.iter().enumerate() { + if num_good == 0 { + residual_.push(i); + residual_indices[i] = Some(residual_.len() - 1); + } + } + + let clauses_ratio = challenge.difficulty.clauses_to_variables_percent as f64; + let num_vars = challenge.difficulty.num_variables as f64; + let max_fuel = 2000000000.0; + let base_fuel = (2000.0 + 40.0 * clauses_ratio) * num_vars; + let flip_fuel = 350.0 + 0.9 * clauses_ratio; + let max_num_rounds = ((max_fuel - base_fuel) / flip_fuel) as usize; + loop { + if !residual_.is_empty() { + + let rand_val = rng.gen::(); + + let i = residual_[rand_val % residual_.len()]; + let mut min_sad = clauses.len(); + let mut v_min_sad = usize::MAX; + let c = &mut clauses[i]; + + if c.len() > 1 { + let random_index = rand_val % c.len(); + c.swap(0, random_index); + } + for &l in c.iter() { + let abs_l = l.abs() as usize - 1; + let clauses_to_check = if variables[abs_l] { &p_clauses[abs_l] } else { &n_clauses[abs_l] }; + + let mut sad = 0; + for &c in clauses_to_check { + if num_good_so_far[c] == 1 { + sad += 1; + } + } + + if sad < min_sad { + min_sad = sad; + v_min_sad = abs_l; + } + } + + let v = if min_sad == 0 { + v_min_sad + } else if rng.gen_bool(0.5) { + c[0].abs() as usize - 1 + } else { + v_min_sad + }; + + if variables[v] { + for &c in &n_clauses[v] { + num_good_so_far[c] += 1; + if num_good_so_far[c] == 1 { + let i = residual_indices[c].take().unwrap(); + let last = residual_.pop().unwrap(); + if i < residual_.len() { + residual_[i] = last; + residual_indices[last] = Some(i); + } + } + } + for &c in &p_clauses[v] { + if num_good_so_far[c] == 1 { + residual_.push(c); + residual_indices[c] = Some(residual_.len() - 1); + } + num_good_so_far[c] -= 1; + } + } else { + for &c in &n_clauses[v] { + if num_good_so_far[c] == 1 { + residual_.push(c); + residual_indices[c] = Some(residual_.len() - 1); + } + num_good_so_far[c] -= 1; + } + + for &c in &p_clauses[v] { + num_good_so_far[c] += 1; + if num_good_so_far[c] == 1 { + let i = residual_indices[c].take().unwrap(); + let last = residual_.pop().unwrap(); + if i < residual_.len() { + residual_[i] = last; + residual_indices[last] = Some(i); + } + } + } + } + + variables[v] = !variables[v]; + } else { + break; + } + rounds += 1; + if rounds >= max_num_rounds { + return Ok(None); + } + } + return Ok(Some(Solution { variables })); +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; diff --git a/tig-algorithms/src/vector_search/invector_hybrid/benchmarker_outbound.rs b/tig-algorithms/src/vector_search/invector_hybrid/benchmarker_outbound.rs new file mode 100644 index 00000000..6a1aa55e --- /dev/null +++ b/tig-algorithms/src/vector_search/invector_hybrid/benchmarker_outbound.rs @@ -0,0 +1,395 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Benchmarker Outbound Game License v1.0 (the "License"); you +may not use this file except in compliance with the License. You may obtain a copy +of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + + +use anyhow::Ok; +use tig_challenges::vector_search::*; +use std::cmp::Ordering; +use std::collections::BinaryHeap; + +struct KDNode<'a> { + point: &'a [f32], + left: Option>>, + right: Option>>, + index: usize, +} + +impl<'a> KDNode<'a> { + fn new(point: &'a [f32], index: usize) -> Self { + KDNode { + point, + left: None, + right: None, + index, + } + } +} +fn quickselect_by(arr: &mut [(&[f32], usize)], k: usize, compare: &F) +where + F: Fn(&(&[f32], usize), &(&[f32], usize)) -> Ordering, +{ + if arr.len() <= 1 { + return; + } + + let pivot_index = partition(arr, compare); + if k < pivot_index { + quickselect_by(&mut arr[..pivot_index], k, compare); + } else if k > pivot_index { + quickselect_by(&mut arr[pivot_index + 1..], k - pivot_index - 1, compare); + } +} + +fn partition(arr: &mut [(&[f32], usize)], compare: &F) -> usize +where + F: Fn(&(&[f32], usize), &(&[f32], usize)) -> Ordering, +{ + let pivot_index = arr.len() >> 1; + arr.swap(pivot_index, arr.len() - 1); + + let mut store_index = 0; + for i in 0..arr.len() - 1 { + if compare(&arr[i], &arr[arr.len() - 1]) == Ordering::Less { + arr.swap(i, store_index); + store_index += 1; + } + } + arr.swap(store_index, arr.len() - 1); + store_index +} + +fn build_kd_tree<'a>(points: &mut [(&'a [f32], usize)]) -> Option>> { + if points.is_empty() { + return None; + } + + const NUM_DIMENSIONS: usize = 250; + let mut stack: Vec<(usize, usize, usize, Option<*mut KDNode<'a>>, bool)> = Vec::new(); + let mut root: Option>> = None; + + stack.push((0, points.len(), 0, None, false)); + + while let Some((start, end, depth, parent_ptr, is_left)) = stack.pop() { + if start >= end { + continue; + } + + let axis = depth % NUM_DIMENSIONS; + let median = (start + end) / 2; + quickselect_by(&mut points[start..end], median - start, &|a, b| { + a.0[axis].partial_cmp(&b.0[axis]).unwrap() + }); + + let (median_point, median_index) = points[median]; + let mut new_node = Box::new(KDNode::new(median_point, median_index)); + let new_node_ptr: *mut KDNode = &mut *new_node; + + if let Some(parent_ptr) = parent_ptr { + unsafe { + if is_left { + (*parent_ptr).left = Some(new_node); + } else { + (*parent_ptr).right = Some(new_node); + } + } + } else { + root = Some(new_node); + } + + stack.push((median + 1, end, depth + 1, Some(new_node_ptr), false)); + stack.push((start, median, depth + 1, Some(new_node_ptr), true)); + } + + root +} + +#[inline(always)] +fn squared_euclidean_distance(a: &[f32], b: &[f32]) -> f32 { + let mut sum = 0.0; + let mut i = 0; + let len = a.len(); + + if a.len() != b.len() || a.len() < 8 { + return f32::MAX; + } + + while i + 7 < len { + unsafe { + let diff0 = *a.get_unchecked(i) - *b.get_unchecked(i); + let diff1 = *a.get_unchecked(i + 1) - *b.get_unchecked(i + 1); + let diff2 = *a.get_unchecked(i + 2) - *b.get_unchecked(i + 2); + let diff3 = *a.get_unchecked(i + 3) - *b.get_unchecked(i + 3); + let diff4 = *a.get_unchecked(i + 4) - *b.get_unchecked(i + 4); + let diff5 = *a.get_unchecked(i + 5) - *b.get_unchecked(i + 5); + let diff6 = *a.get_unchecked(i + 6) - *b.get_unchecked(i + 6); + let diff7 = *a.get_unchecked(i + 7) - *b.get_unchecked(i + 7); + + sum += diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3 + + diff4 * diff4 + diff5 * diff5 + diff6 * diff6 + diff7 * diff7; + } + + i += 8; + } + + while i < len { + unsafe { + let diff = *a.get_unchecked(i) - *b.get_unchecked(i); + sum += diff * diff; + } + i += 1; + } + sum +} + +#[inline(always)] +fn early_stopping_distance(a: &[f32], b: &[f32], current_min: f32) -> f32 { + let mut sum = 0.0; + let mut i = 0; + let len = a.len(); + + if a.len() != b.len() || a.len() < 8 { + return f32::MAX; + } + + while i + 7 < len { + unsafe { + let diff0 = *a.get_unchecked(i) - *b.get_unchecked(i); + let diff1 = *a.get_unchecked(i + 1) - *b.get_unchecked(i + 1); + let diff2 = *a.get_unchecked(i + 2) - *b.get_unchecked(i + 2); + let diff3 = *a.get_unchecked(i + 3) - *b.get_unchecked(i + 3); + let diff4 = *a.get_unchecked(i + 4) - *b.get_unchecked(i + 4); + let diff5 = *a.get_unchecked(i + 5) - *b.get_unchecked(i + 5); + let diff6 = *a.get_unchecked(i + 6) - *b.get_unchecked(i + 6); + let diff7 = *a.get_unchecked(i + 7) - *b.get_unchecked(i + 7); + + sum += diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3 + + diff4 * diff4 + diff5 * diff5 + diff6 * diff6 + diff7 * diff7; + } + + if sum > current_min { + return f32::MAX; + } + + i += 8; + } + + while i < len { + unsafe { + let diff = *a.get_unchecked(i) - *b.get_unchecked(i); + sum += diff * diff; + } + i += 1; + } + sum +} + +fn nearest_neighbor_search<'a>( + root: &Option>>, + target: &[f32], + best: &mut (f32, Option), +) { + let num_dimensions = target.len(); + let mut stack = Vec::with_capacity(64); + + if let Some(node) = root { + stack.push((node.as_ref(), 0)); + } + + while let Some((node, depth)) = stack.pop() { + let axis = depth % num_dimensions; + let dist = early_stopping_distance(&node.point, target, best.0); + + if dist < best.0 { + best.0 = dist; + best.1 = Some(node.index); + } + + let diff = target[axis] - node.point[axis]; + let sqr_diff = diff * diff; + + let (nearer, farther) = if diff < 0.0 { + (&node.left, &node.right) + } else { + (&node.right, &node.left) + }; + + if let Some(nearer_node) = nearer { + stack.push((nearer_node.as_ref(), depth + 1)); + } + + if sqr_diff < best.0 { + if let Some(farther_node) = farther { + stack.push((farther_node.as_ref(), depth + 1)); + } + } + } +} +fn calculate_mean_vector(vectors: &[&[f32]]) -> Vec { + let num_vectors = vectors.len(); + let num_dimensions = 250; + + let mut mean_vector = vec![0.0f64; num_dimensions]; + + for vector in vectors { + for i in 0..num_dimensions { + mean_vector[i] += vector[i] as f64; + } + } + for i in 0..num_dimensions { + mean_vector[i] /= num_vectors as f64; + } + mean_vector.into_iter().map(|x| x as f32).collect() +} + +#[derive(Debug)] +struct FloatOrd(f32); + +impl PartialEq for FloatOrd { + fn eq(&self, other: &Self) -> bool { + self.0 == other.0 + } +} + +impl Eq for FloatOrd {} + +impl PartialOrd for FloatOrd { + fn partial_cmp(&self, other: &Self) -> Option { + self.0.partial_cmp(&other.0) + } +} + +impl Ord for FloatOrd { + fn cmp(&self, other: &Self) -> Ordering { + + self.partial_cmp(other).unwrap_or(Ordering::Equal) + } +} + +fn filter_relevant_vectors<'a>( + database: &'a [Vec], + query_vectors: &[Vec], + k: usize, +) -> Vec<(f32, &'a [f32], usize)> { + let query_refs: Vec<&[f32]> = query_vectors.iter().map(|v| &v[..]).collect(); + let mean_query_vector = calculate_mean_vector(&query_refs); + + let mut heap: BinaryHeap<(FloatOrd, usize)> = BinaryHeap::with_capacity(k); + + for (index, vector) in database.iter().enumerate() { + if heap.len() < k + { + let dist = squared_euclidean_distance(&mean_query_vector, vector); + let ord_dist = FloatOrd(dist); + + heap.push((ord_dist, index)); + } else if let Some(&(FloatOrd(top_dist), _)) = heap.peek() + { + let dist = early_stopping_distance(&mean_query_vector, vector, top_dist); + let ord_dist = FloatOrd(dist); + if dist < top_dist { + heap.pop(); + heap.push((ord_dist, index)); + } + } + } + heap.into_sorted_vec() + .into_iter() + .map(|(FloatOrd(dist), index)| (dist, &database[index][..], index)) + .collect() +} + +pub fn solve_challenge(challenge: &Challenge) -> anyhow::Result> { + let query_count = challenge.query_vectors.len(); + + let max_fuel = 2000000000.0; + let base_fuel = 760000000.0; + let alpha = 1700.0 * challenge.difficulty.num_queries as f64; + + let m = ((max_fuel - base_fuel) / alpha) as usize; + let n = (m as f32 * 1.2) as usize; + let r = n - m; + + let closest_vectors = filter_relevant_vectors( + &challenge.vector_database, + &challenge.query_vectors, + n, + ); + + let (m_slice, r_slice) = closest_vectors.split_at(m); + let m_vectors: Vec<_> = m_slice.to_vec(); + let r_vectors: Vec<_> = r_slice.to_vec(); + + let mut kd_tree_vectors: Vec<(&[f32], usize)> = m_vectors.iter().map(|&(_, v, i)| (v, i)).collect(); + let kd_tree = build_kd_tree(&mut kd_tree_vectors); + + let mut best_indexes = Vec::with_capacity(query_count); + let mut distances = Vec::with_capacity(query_count); + + for query in &challenge.query_vectors { + let mut best = (std::f32::MAX, None); + nearest_neighbor_search(&kd_tree, query, &mut best); + + distances.push(best.0); + best_indexes.push(best.1.unwrap_or(0)); + } + + let brute_force_count = (query_count as f32 * 0.1) as usize; + let mut distance_indices: Vec<_> = distances.iter().enumerate().collect(); + distance_indices.sort_unstable_by(|a, b| b.1.partial_cmp(a.1).unwrap()); + let high_distance_indices: Vec<_> = distance_indices.into_iter() + .take(brute_force_count) + .map(|(index, _)| index) + .collect(); + + for &query_index in &high_distance_indices { + let query = &challenge.query_vectors[query_index]; + let mut best = (distances[query_index], best_indexes[query_index]); + + for &(_, vec, index) in &r_vectors { + let dist = squared_euclidean_distance(query, vec); + if dist < best.0 { + best = (dist, index); + } + } + + best_indexes[query_index] = best.1; + } + + Ok(Some(Solution { + indexes: best_indexes, + })) +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; diff --git a/tig-algorithms/src/vector_search/invector_hybrid/commercial.rs b/tig-algorithms/src/vector_search/invector_hybrid/commercial.rs new file mode 100644 index 00000000..c82abc35 --- /dev/null +++ b/tig-algorithms/src/vector_search/invector_hybrid/commercial.rs @@ -0,0 +1,395 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Commercial License v1.0 (the "License"); you +may not use this file except in compliance with the License. You may obtain a copy +of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + + +use anyhow::Ok; +use tig_challenges::vector_search::*; +use std::cmp::Ordering; +use std::collections::BinaryHeap; + +struct KDNode<'a> { + point: &'a [f32], + left: Option>>, + right: Option>>, + index: usize, +} + +impl<'a> KDNode<'a> { + fn new(point: &'a [f32], index: usize) -> Self { + KDNode { + point, + left: None, + right: None, + index, + } + } +} +fn quickselect_by(arr: &mut [(&[f32], usize)], k: usize, compare: &F) +where + F: Fn(&(&[f32], usize), &(&[f32], usize)) -> Ordering, +{ + if arr.len() <= 1 { + return; + } + + let pivot_index = partition(arr, compare); + if k < pivot_index { + quickselect_by(&mut arr[..pivot_index], k, compare); + } else if k > pivot_index { + quickselect_by(&mut arr[pivot_index + 1..], k - pivot_index - 1, compare); + } +} + +fn partition(arr: &mut [(&[f32], usize)], compare: &F) -> usize +where + F: Fn(&(&[f32], usize), &(&[f32], usize)) -> Ordering, +{ + let pivot_index = arr.len() >> 1; + arr.swap(pivot_index, arr.len() - 1); + + let mut store_index = 0; + for i in 0..arr.len() - 1 { + if compare(&arr[i], &arr[arr.len() - 1]) == Ordering::Less { + arr.swap(i, store_index); + store_index += 1; + } + } + arr.swap(store_index, arr.len() - 1); + store_index +} + +fn build_kd_tree<'a>(points: &mut [(&'a [f32], usize)]) -> Option>> { + if points.is_empty() { + return None; + } + + const NUM_DIMENSIONS: usize = 250; + let mut stack: Vec<(usize, usize, usize, Option<*mut KDNode<'a>>, bool)> = Vec::new(); + let mut root: Option>> = None; + + stack.push((0, points.len(), 0, None, false)); + + while let Some((start, end, depth, parent_ptr, is_left)) = stack.pop() { + if start >= end { + continue; + } + + let axis = depth % NUM_DIMENSIONS; + let median = (start + end) / 2; + quickselect_by(&mut points[start..end], median - start, &|a, b| { + a.0[axis].partial_cmp(&b.0[axis]).unwrap() + }); + + let (median_point, median_index) = points[median]; + let mut new_node = Box::new(KDNode::new(median_point, median_index)); + let new_node_ptr: *mut KDNode = &mut *new_node; + + if let Some(parent_ptr) = parent_ptr { + unsafe { + if is_left { + (*parent_ptr).left = Some(new_node); + } else { + (*parent_ptr).right = Some(new_node); + } + } + } else { + root = Some(new_node); + } + + stack.push((median + 1, end, depth + 1, Some(new_node_ptr), false)); + stack.push((start, median, depth + 1, Some(new_node_ptr), true)); + } + + root +} + +#[inline(always)] +fn squared_euclidean_distance(a: &[f32], b: &[f32]) -> f32 { + let mut sum = 0.0; + let mut i = 0; + let len = a.len(); + + if a.len() != b.len() || a.len() < 8 { + return f32::MAX; + } + + while i + 7 < len { + unsafe { + let diff0 = *a.get_unchecked(i) - *b.get_unchecked(i); + let diff1 = *a.get_unchecked(i + 1) - *b.get_unchecked(i + 1); + let diff2 = *a.get_unchecked(i + 2) - *b.get_unchecked(i + 2); + let diff3 = *a.get_unchecked(i + 3) - *b.get_unchecked(i + 3); + let diff4 = *a.get_unchecked(i + 4) - *b.get_unchecked(i + 4); + let diff5 = *a.get_unchecked(i + 5) - *b.get_unchecked(i + 5); + let diff6 = *a.get_unchecked(i + 6) - *b.get_unchecked(i + 6); + let diff7 = *a.get_unchecked(i + 7) - *b.get_unchecked(i + 7); + + sum += diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3 + + diff4 * diff4 + diff5 * diff5 + diff6 * diff6 + diff7 * diff7; + } + + i += 8; + } + + while i < len { + unsafe { + let diff = *a.get_unchecked(i) - *b.get_unchecked(i); + sum += diff * diff; + } + i += 1; + } + sum +} + +#[inline(always)] +fn early_stopping_distance(a: &[f32], b: &[f32], current_min: f32) -> f32 { + let mut sum = 0.0; + let mut i = 0; + let len = a.len(); + + if a.len() != b.len() || a.len() < 8 { + return f32::MAX; + } + + while i + 7 < len { + unsafe { + let diff0 = *a.get_unchecked(i) - *b.get_unchecked(i); + let diff1 = *a.get_unchecked(i + 1) - *b.get_unchecked(i + 1); + let diff2 = *a.get_unchecked(i + 2) - *b.get_unchecked(i + 2); + let diff3 = *a.get_unchecked(i + 3) - *b.get_unchecked(i + 3); + let diff4 = *a.get_unchecked(i + 4) - *b.get_unchecked(i + 4); + let diff5 = *a.get_unchecked(i + 5) - *b.get_unchecked(i + 5); + let diff6 = *a.get_unchecked(i + 6) - *b.get_unchecked(i + 6); + let diff7 = *a.get_unchecked(i + 7) - *b.get_unchecked(i + 7); + + sum += diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3 + + diff4 * diff4 + diff5 * diff5 + diff6 * diff6 + diff7 * diff7; + } + + if sum > current_min { + return f32::MAX; + } + + i += 8; + } + + while i < len { + unsafe { + let diff = *a.get_unchecked(i) - *b.get_unchecked(i); + sum += diff * diff; + } + i += 1; + } + sum +} + +fn nearest_neighbor_search<'a>( + root: &Option>>, + target: &[f32], + best: &mut (f32, Option), +) { + let num_dimensions = target.len(); + let mut stack = Vec::with_capacity(64); + + if let Some(node) = root { + stack.push((node.as_ref(), 0)); + } + + while let Some((node, depth)) = stack.pop() { + let axis = depth % num_dimensions; + let dist = early_stopping_distance(&node.point, target, best.0); + + if dist < best.0 { + best.0 = dist; + best.1 = Some(node.index); + } + + let diff = target[axis] - node.point[axis]; + let sqr_diff = diff * diff; + + let (nearer, farther) = if diff < 0.0 { + (&node.left, &node.right) + } else { + (&node.right, &node.left) + }; + + if let Some(nearer_node) = nearer { + stack.push((nearer_node.as_ref(), depth + 1)); + } + + if sqr_diff < best.0 { + if let Some(farther_node) = farther { + stack.push((farther_node.as_ref(), depth + 1)); + } + } + } +} +fn calculate_mean_vector(vectors: &[&[f32]]) -> Vec { + let num_vectors = vectors.len(); + let num_dimensions = 250; + + let mut mean_vector = vec![0.0f64; num_dimensions]; + + for vector in vectors { + for i in 0..num_dimensions { + mean_vector[i] += vector[i] as f64; + } + } + for i in 0..num_dimensions { + mean_vector[i] /= num_vectors as f64; + } + mean_vector.into_iter().map(|x| x as f32).collect() +} + +#[derive(Debug)] +struct FloatOrd(f32); + +impl PartialEq for FloatOrd { + fn eq(&self, other: &Self) -> bool { + self.0 == other.0 + } +} + +impl Eq for FloatOrd {} + +impl PartialOrd for FloatOrd { + fn partial_cmp(&self, other: &Self) -> Option { + self.0.partial_cmp(&other.0) + } +} + +impl Ord for FloatOrd { + fn cmp(&self, other: &Self) -> Ordering { + + self.partial_cmp(other).unwrap_or(Ordering::Equal) + } +} + +fn filter_relevant_vectors<'a>( + database: &'a [Vec], + query_vectors: &[Vec], + k: usize, +) -> Vec<(f32, &'a [f32], usize)> { + let query_refs: Vec<&[f32]> = query_vectors.iter().map(|v| &v[..]).collect(); + let mean_query_vector = calculate_mean_vector(&query_refs); + + let mut heap: BinaryHeap<(FloatOrd, usize)> = BinaryHeap::with_capacity(k); + + for (index, vector) in database.iter().enumerate() { + if heap.len() < k + { + let dist = squared_euclidean_distance(&mean_query_vector, vector); + let ord_dist = FloatOrd(dist); + + heap.push((ord_dist, index)); + } else if let Some(&(FloatOrd(top_dist), _)) = heap.peek() + { + let dist = early_stopping_distance(&mean_query_vector, vector, top_dist); + let ord_dist = FloatOrd(dist); + if dist < top_dist { + heap.pop(); + heap.push((ord_dist, index)); + } + } + } + heap.into_sorted_vec() + .into_iter() + .map(|(FloatOrd(dist), index)| (dist, &database[index][..], index)) + .collect() +} + +pub fn solve_challenge(challenge: &Challenge) -> anyhow::Result> { + let query_count = challenge.query_vectors.len(); + + let max_fuel = 2000000000.0; + let base_fuel = 760000000.0; + let alpha = 1700.0 * challenge.difficulty.num_queries as f64; + + let m = ((max_fuel - base_fuel) / alpha) as usize; + let n = (m as f32 * 1.2) as usize; + let r = n - m; + + let closest_vectors = filter_relevant_vectors( + &challenge.vector_database, + &challenge.query_vectors, + n, + ); + + let (m_slice, r_slice) = closest_vectors.split_at(m); + let m_vectors: Vec<_> = m_slice.to_vec(); + let r_vectors: Vec<_> = r_slice.to_vec(); + + let mut kd_tree_vectors: Vec<(&[f32], usize)> = m_vectors.iter().map(|&(_, v, i)| (v, i)).collect(); + let kd_tree = build_kd_tree(&mut kd_tree_vectors); + + let mut best_indexes = Vec::with_capacity(query_count); + let mut distances = Vec::with_capacity(query_count); + + for query in &challenge.query_vectors { + let mut best = (std::f32::MAX, None); + nearest_neighbor_search(&kd_tree, query, &mut best); + + distances.push(best.0); + best_indexes.push(best.1.unwrap_or(0)); + } + + let brute_force_count = (query_count as f32 * 0.1) as usize; + let mut distance_indices: Vec<_> = distances.iter().enumerate().collect(); + distance_indices.sort_unstable_by(|a, b| b.1.partial_cmp(a.1).unwrap()); + let high_distance_indices: Vec<_> = distance_indices.into_iter() + .take(brute_force_count) + .map(|(index, _)| index) + .collect(); + + for &query_index in &high_distance_indices { + let query = &challenge.query_vectors[query_index]; + let mut best = (distances[query_index], best_indexes[query_index]); + + for &(_, vec, index) in &r_vectors { + let dist = squared_euclidean_distance(query, vec); + if dist < best.0 { + best = (dist, index); + } + } + + best_indexes[query_index] = best.1; + } + + Ok(Some(Solution { + indexes: best_indexes, + })) +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; diff --git a/tig-algorithms/src/vector_search/invector_hybrid/inbound.rs b/tig-algorithms/src/vector_search/invector_hybrid/inbound.rs new file mode 100644 index 00000000..9a283e9d --- /dev/null +++ b/tig-algorithms/src/vector_search/invector_hybrid/inbound.rs @@ -0,0 +1,395 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Inbound Game License v1.0 or (at your option) any later +version (the "License"); you may not use this file except in compliance with the +License. You may obtain a copy of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + + +use anyhow::Ok; +use tig_challenges::vector_search::*; +use std::cmp::Ordering; +use std::collections::BinaryHeap; + +struct KDNode<'a> { + point: &'a [f32], + left: Option>>, + right: Option>>, + index: usize, +} + +impl<'a> KDNode<'a> { + fn new(point: &'a [f32], index: usize) -> Self { + KDNode { + point, + left: None, + right: None, + index, + } + } +} +fn quickselect_by(arr: &mut [(&[f32], usize)], k: usize, compare: &F) +where + F: Fn(&(&[f32], usize), &(&[f32], usize)) -> Ordering, +{ + if arr.len() <= 1 { + return; + } + + let pivot_index = partition(arr, compare); + if k < pivot_index { + quickselect_by(&mut arr[..pivot_index], k, compare); + } else if k > pivot_index { + quickselect_by(&mut arr[pivot_index + 1..], k - pivot_index - 1, compare); + } +} + +fn partition(arr: &mut [(&[f32], usize)], compare: &F) -> usize +where + F: Fn(&(&[f32], usize), &(&[f32], usize)) -> Ordering, +{ + let pivot_index = arr.len() >> 1; + arr.swap(pivot_index, arr.len() - 1); + + let mut store_index = 0; + for i in 0..arr.len() - 1 { + if compare(&arr[i], &arr[arr.len() - 1]) == Ordering::Less { + arr.swap(i, store_index); + store_index += 1; + } + } + arr.swap(store_index, arr.len() - 1); + store_index +} + +fn build_kd_tree<'a>(points: &mut [(&'a [f32], usize)]) -> Option>> { + if points.is_empty() { + return None; + } + + const NUM_DIMENSIONS: usize = 250; + let mut stack: Vec<(usize, usize, usize, Option<*mut KDNode<'a>>, bool)> = Vec::new(); + let mut root: Option>> = None; + + stack.push((0, points.len(), 0, None, false)); + + while let Some((start, end, depth, parent_ptr, is_left)) = stack.pop() { + if start >= end { + continue; + } + + let axis = depth % NUM_DIMENSIONS; + let median = (start + end) / 2; + quickselect_by(&mut points[start..end], median - start, &|a, b| { + a.0[axis].partial_cmp(&b.0[axis]).unwrap() + }); + + let (median_point, median_index) = points[median]; + let mut new_node = Box::new(KDNode::new(median_point, median_index)); + let new_node_ptr: *mut KDNode = &mut *new_node; + + if let Some(parent_ptr) = parent_ptr { + unsafe { + if is_left { + (*parent_ptr).left = Some(new_node); + } else { + (*parent_ptr).right = Some(new_node); + } + } + } else { + root = Some(new_node); + } + + stack.push((median + 1, end, depth + 1, Some(new_node_ptr), false)); + stack.push((start, median, depth + 1, Some(new_node_ptr), true)); + } + + root +} + +#[inline(always)] +fn squared_euclidean_distance(a: &[f32], b: &[f32]) -> f32 { + let mut sum = 0.0; + let mut i = 0; + let len = a.len(); + + if a.len() != b.len() || a.len() < 8 { + return f32::MAX; + } + + while i + 7 < len { + unsafe { + let diff0 = *a.get_unchecked(i) - *b.get_unchecked(i); + let diff1 = *a.get_unchecked(i + 1) - *b.get_unchecked(i + 1); + let diff2 = *a.get_unchecked(i + 2) - *b.get_unchecked(i + 2); + let diff3 = *a.get_unchecked(i + 3) - *b.get_unchecked(i + 3); + let diff4 = *a.get_unchecked(i + 4) - *b.get_unchecked(i + 4); + let diff5 = *a.get_unchecked(i + 5) - *b.get_unchecked(i + 5); + let diff6 = *a.get_unchecked(i + 6) - *b.get_unchecked(i + 6); + let diff7 = *a.get_unchecked(i + 7) - *b.get_unchecked(i + 7); + + sum += diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3 + + diff4 * diff4 + diff5 * diff5 + diff6 * diff6 + diff7 * diff7; + } + + i += 8; + } + + while i < len { + unsafe { + let diff = *a.get_unchecked(i) - *b.get_unchecked(i); + sum += diff * diff; + } + i += 1; + } + sum +} + +#[inline(always)] +fn early_stopping_distance(a: &[f32], b: &[f32], current_min: f32) -> f32 { + let mut sum = 0.0; + let mut i = 0; + let len = a.len(); + + if a.len() != b.len() || a.len() < 8 { + return f32::MAX; + } + + while i + 7 < len { + unsafe { + let diff0 = *a.get_unchecked(i) - *b.get_unchecked(i); + let diff1 = *a.get_unchecked(i + 1) - *b.get_unchecked(i + 1); + let diff2 = *a.get_unchecked(i + 2) - *b.get_unchecked(i + 2); + let diff3 = *a.get_unchecked(i + 3) - *b.get_unchecked(i + 3); + let diff4 = *a.get_unchecked(i + 4) - *b.get_unchecked(i + 4); + let diff5 = *a.get_unchecked(i + 5) - *b.get_unchecked(i + 5); + let diff6 = *a.get_unchecked(i + 6) - *b.get_unchecked(i + 6); + let diff7 = *a.get_unchecked(i + 7) - *b.get_unchecked(i + 7); + + sum += diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3 + + diff4 * diff4 + diff5 * diff5 + diff6 * diff6 + diff7 * diff7; + } + + if sum > current_min { + return f32::MAX; + } + + i += 8; + } + + while i < len { + unsafe { + let diff = *a.get_unchecked(i) - *b.get_unchecked(i); + sum += diff * diff; + } + i += 1; + } + sum +} + +fn nearest_neighbor_search<'a>( + root: &Option>>, + target: &[f32], + best: &mut (f32, Option), +) { + let num_dimensions = target.len(); + let mut stack = Vec::with_capacity(64); + + if let Some(node) = root { + stack.push((node.as_ref(), 0)); + } + + while let Some((node, depth)) = stack.pop() { + let axis = depth % num_dimensions; + let dist = early_stopping_distance(&node.point, target, best.0); + + if dist < best.0 { + best.0 = dist; + best.1 = Some(node.index); + } + + let diff = target[axis] - node.point[axis]; + let sqr_diff = diff * diff; + + let (nearer, farther) = if diff < 0.0 { + (&node.left, &node.right) + } else { + (&node.right, &node.left) + }; + + if let Some(nearer_node) = nearer { + stack.push((nearer_node.as_ref(), depth + 1)); + } + + if sqr_diff < best.0 { + if let Some(farther_node) = farther { + stack.push((farther_node.as_ref(), depth + 1)); + } + } + } +} +fn calculate_mean_vector(vectors: &[&[f32]]) -> Vec { + let num_vectors = vectors.len(); + let num_dimensions = 250; + + let mut mean_vector = vec![0.0f64; num_dimensions]; + + for vector in vectors { + for i in 0..num_dimensions { + mean_vector[i] += vector[i] as f64; + } + } + for i in 0..num_dimensions { + mean_vector[i] /= num_vectors as f64; + } + mean_vector.into_iter().map(|x| x as f32).collect() +} + +#[derive(Debug)] +struct FloatOrd(f32); + +impl PartialEq for FloatOrd { + fn eq(&self, other: &Self) -> bool { + self.0 == other.0 + } +} + +impl Eq for FloatOrd {} + +impl PartialOrd for FloatOrd { + fn partial_cmp(&self, other: &Self) -> Option { + self.0.partial_cmp(&other.0) + } +} + +impl Ord for FloatOrd { + fn cmp(&self, other: &Self) -> Ordering { + + self.partial_cmp(other).unwrap_or(Ordering::Equal) + } +} + +fn filter_relevant_vectors<'a>( + database: &'a [Vec], + query_vectors: &[Vec], + k: usize, +) -> Vec<(f32, &'a [f32], usize)> { + let query_refs: Vec<&[f32]> = query_vectors.iter().map(|v| &v[..]).collect(); + let mean_query_vector = calculate_mean_vector(&query_refs); + + let mut heap: BinaryHeap<(FloatOrd, usize)> = BinaryHeap::with_capacity(k); + + for (index, vector) in database.iter().enumerate() { + if heap.len() < k + { + let dist = squared_euclidean_distance(&mean_query_vector, vector); + let ord_dist = FloatOrd(dist); + + heap.push((ord_dist, index)); + } else if let Some(&(FloatOrd(top_dist), _)) = heap.peek() + { + let dist = early_stopping_distance(&mean_query_vector, vector, top_dist); + let ord_dist = FloatOrd(dist); + if dist < top_dist { + heap.pop(); + heap.push((ord_dist, index)); + } + } + } + heap.into_sorted_vec() + .into_iter() + .map(|(FloatOrd(dist), index)| (dist, &database[index][..], index)) + .collect() +} + +pub fn solve_challenge(challenge: &Challenge) -> anyhow::Result> { + let query_count = challenge.query_vectors.len(); + + let max_fuel = 2000000000.0; + let base_fuel = 760000000.0; + let alpha = 1700.0 * challenge.difficulty.num_queries as f64; + + let m = ((max_fuel - base_fuel) / alpha) as usize; + let n = (m as f32 * 1.2) as usize; + let r = n - m; + + let closest_vectors = filter_relevant_vectors( + &challenge.vector_database, + &challenge.query_vectors, + n, + ); + + let (m_slice, r_slice) = closest_vectors.split_at(m); + let m_vectors: Vec<_> = m_slice.to_vec(); + let r_vectors: Vec<_> = r_slice.to_vec(); + + let mut kd_tree_vectors: Vec<(&[f32], usize)> = m_vectors.iter().map(|&(_, v, i)| (v, i)).collect(); + let kd_tree = build_kd_tree(&mut kd_tree_vectors); + + let mut best_indexes = Vec::with_capacity(query_count); + let mut distances = Vec::with_capacity(query_count); + + for query in &challenge.query_vectors { + let mut best = (std::f32::MAX, None); + nearest_neighbor_search(&kd_tree, query, &mut best); + + distances.push(best.0); + best_indexes.push(best.1.unwrap_or(0)); + } + + let brute_force_count = (query_count as f32 * 0.1) as usize; + let mut distance_indices: Vec<_> = distances.iter().enumerate().collect(); + distance_indices.sort_unstable_by(|a, b| b.1.partial_cmp(a.1).unwrap()); + let high_distance_indices: Vec<_> = distance_indices.into_iter() + .take(brute_force_count) + .map(|(index, _)| index) + .collect(); + + for &query_index in &high_distance_indices { + let query = &challenge.query_vectors[query_index]; + let mut best = (distances[query_index], best_indexes[query_index]); + + for &(_, vec, index) in &r_vectors { + let dist = squared_euclidean_distance(query, vec); + if dist < best.0 { + best = (dist, index); + } + } + + best_indexes[query_index] = best.1; + } + + Ok(Some(Solution { + indexes: best_indexes, + })) +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; diff --git a/tig-algorithms/src/vector_search/invector_hybrid/innovator_outbound.rs b/tig-algorithms/src/vector_search/invector_hybrid/innovator_outbound.rs new file mode 100644 index 00000000..624671f6 --- /dev/null +++ b/tig-algorithms/src/vector_search/invector_hybrid/innovator_outbound.rs @@ -0,0 +1,395 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Innovator Outbound Game License v1.0 (the "License"); you +may not use this file except in compliance with the License. You may obtain a copy +of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + + +use anyhow::Ok; +use tig_challenges::vector_search::*; +use std::cmp::Ordering; +use std::collections::BinaryHeap; + +struct KDNode<'a> { + point: &'a [f32], + left: Option>>, + right: Option>>, + index: usize, +} + +impl<'a> KDNode<'a> { + fn new(point: &'a [f32], index: usize) -> Self { + KDNode { + point, + left: None, + right: None, + index, + } + } +} +fn quickselect_by(arr: &mut [(&[f32], usize)], k: usize, compare: &F) +where + F: Fn(&(&[f32], usize), &(&[f32], usize)) -> Ordering, +{ + if arr.len() <= 1 { + return; + } + + let pivot_index = partition(arr, compare); + if k < pivot_index { + quickselect_by(&mut arr[..pivot_index], k, compare); + } else if k > pivot_index { + quickselect_by(&mut arr[pivot_index + 1..], k - pivot_index - 1, compare); + } +} + +fn partition(arr: &mut [(&[f32], usize)], compare: &F) -> usize +where + F: Fn(&(&[f32], usize), &(&[f32], usize)) -> Ordering, +{ + let pivot_index = arr.len() >> 1; + arr.swap(pivot_index, arr.len() - 1); + + let mut store_index = 0; + for i in 0..arr.len() - 1 { + if compare(&arr[i], &arr[arr.len() - 1]) == Ordering::Less { + arr.swap(i, store_index); + store_index += 1; + } + } + arr.swap(store_index, arr.len() - 1); + store_index +} + +fn build_kd_tree<'a>(points: &mut [(&'a [f32], usize)]) -> Option>> { + if points.is_empty() { + return None; + } + + const NUM_DIMENSIONS: usize = 250; + let mut stack: Vec<(usize, usize, usize, Option<*mut KDNode<'a>>, bool)> = Vec::new(); + let mut root: Option>> = None; + + stack.push((0, points.len(), 0, None, false)); + + while let Some((start, end, depth, parent_ptr, is_left)) = stack.pop() { + if start >= end { + continue; + } + + let axis = depth % NUM_DIMENSIONS; + let median = (start + end) / 2; + quickselect_by(&mut points[start..end], median - start, &|a, b| { + a.0[axis].partial_cmp(&b.0[axis]).unwrap() + }); + + let (median_point, median_index) = points[median]; + let mut new_node = Box::new(KDNode::new(median_point, median_index)); + let new_node_ptr: *mut KDNode = &mut *new_node; + + if let Some(parent_ptr) = parent_ptr { + unsafe { + if is_left { + (*parent_ptr).left = Some(new_node); + } else { + (*parent_ptr).right = Some(new_node); + } + } + } else { + root = Some(new_node); + } + + stack.push((median + 1, end, depth + 1, Some(new_node_ptr), false)); + stack.push((start, median, depth + 1, Some(new_node_ptr), true)); + } + + root +} + +#[inline(always)] +fn squared_euclidean_distance(a: &[f32], b: &[f32]) -> f32 { + let mut sum = 0.0; + let mut i = 0; + let len = a.len(); + + if a.len() != b.len() || a.len() < 8 { + return f32::MAX; + } + + while i + 7 < len { + unsafe { + let diff0 = *a.get_unchecked(i) - *b.get_unchecked(i); + let diff1 = *a.get_unchecked(i + 1) - *b.get_unchecked(i + 1); + let diff2 = *a.get_unchecked(i + 2) - *b.get_unchecked(i + 2); + let diff3 = *a.get_unchecked(i + 3) - *b.get_unchecked(i + 3); + let diff4 = *a.get_unchecked(i + 4) - *b.get_unchecked(i + 4); + let diff5 = *a.get_unchecked(i + 5) - *b.get_unchecked(i + 5); + let diff6 = *a.get_unchecked(i + 6) - *b.get_unchecked(i + 6); + let diff7 = *a.get_unchecked(i + 7) - *b.get_unchecked(i + 7); + + sum += diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3 + + diff4 * diff4 + diff5 * diff5 + diff6 * diff6 + diff7 * diff7; + } + + i += 8; + } + + while i < len { + unsafe { + let diff = *a.get_unchecked(i) - *b.get_unchecked(i); + sum += diff * diff; + } + i += 1; + } + sum +} + +#[inline(always)] +fn early_stopping_distance(a: &[f32], b: &[f32], current_min: f32) -> f32 { + let mut sum = 0.0; + let mut i = 0; + let len = a.len(); + + if a.len() != b.len() || a.len() < 8 { + return f32::MAX; + } + + while i + 7 < len { + unsafe { + let diff0 = *a.get_unchecked(i) - *b.get_unchecked(i); + let diff1 = *a.get_unchecked(i + 1) - *b.get_unchecked(i + 1); + let diff2 = *a.get_unchecked(i + 2) - *b.get_unchecked(i + 2); + let diff3 = *a.get_unchecked(i + 3) - *b.get_unchecked(i + 3); + let diff4 = *a.get_unchecked(i + 4) - *b.get_unchecked(i + 4); + let diff5 = *a.get_unchecked(i + 5) - *b.get_unchecked(i + 5); + let diff6 = *a.get_unchecked(i + 6) - *b.get_unchecked(i + 6); + let diff7 = *a.get_unchecked(i + 7) - *b.get_unchecked(i + 7); + + sum += diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3 + + diff4 * diff4 + diff5 * diff5 + diff6 * diff6 + diff7 * diff7; + } + + if sum > current_min { + return f32::MAX; + } + + i += 8; + } + + while i < len { + unsafe { + let diff = *a.get_unchecked(i) - *b.get_unchecked(i); + sum += diff * diff; + } + i += 1; + } + sum +} + +fn nearest_neighbor_search<'a>( + root: &Option>>, + target: &[f32], + best: &mut (f32, Option), +) { + let num_dimensions = target.len(); + let mut stack = Vec::with_capacity(64); + + if let Some(node) = root { + stack.push((node.as_ref(), 0)); + } + + while let Some((node, depth)) = stack.pop() { + let axis = depth % num_dimensions; + let dist = early_stopping_distance(&node.point, target, best.0); + + if dist < best.0 { + best.0 = dist; + best.1 = Some(node.index); + } + + let diff = target[axis] - node.point[axis]; + let sqr_diff = diff * diff; + + let (nearer, farther) = if diff < 0.0 { + (&node.left, &node.right) + } else { + (&node.right, &node.left) + }; + + if let Some(nearer_node) = nearer { + stack.push((nearer_node.as_ref(), depth + 1)); + } + + if sqr_diff < best.0 { + if let Some(farther_node) = farther { + stack.push((farther_node.as_ref(), depth + 1)); + } + } + } +} +fn calculate_mean_vector(vectors: &[&[f32]]) -> Vec { + let num_vectors = vectors.len(); + let num_dimensions = 250; + + let mut mean_vector = vec![0.0f64; num_dimensions]; + + for vector in vectors { + for i in 0..num_dimensions { + mean_vector[i] += vector[i] as f64; + } + } + for i in 0..num_dimensions { + mean_vector[i] /= num_vectors as f64; + } + mean_vector.into_iter().map(|x| x as f32).collect() +} + +#[derive(Debug)] +struct FloatOrd(f32); + +impl PartialEq for FloatOrd { + fn eq(&self, other: &Self) -> bool { + self.0 == other.0 + } +} + +impl Eq for FloatOrd {} + +impl PartialOrd for FloatOrd { + fn partial_cmp(&self, other: &Self) -> Option { + self.0.partial_cmp(&other.0) + } +} + +impl Ord for FloatOrd { + fn cmp(&self, other: &Self) -> Ordering { + + self.partial_cmp(other).unwrap_or(Ordering::Equal) + } +} + +fn filter_relevant_vectors<'a>( + database: &'a [Vec], + query_vectors: &[Vec], + k: usize, +) -> Vec<(f32, &'a [f32], usize)> { + let query_refs: Vec<&[f32]> = query_vectors.iter().map(|v| &v[..]).collect(); + let mean_query_vector = calculate_mean_vector(&query_refs); + + let mut heap: BinaryHeap<(FloatOrd, usize)> = BinaryHeap::with_capacity(k); + + for (index, vector) in database.iter().enumerate() { + if heap.len() < k + { + let dist = squared_euclidean_distance(&mean_query_vector, vector); + let ord_dist = FloatOrd(dist); + + heap.push((ord_dist, index)); + } else if let Some(&(FloatOrd(top_dist), _)) = heap.peek() + { + let dist = early_stopping_distance(&mean_query_vector, vector, top_dist); + let ord_dist = FloatOrd(dist); + if dist < top_dist { + heap.pop(); + heap.push((ord_dist, index)); + } + } + } + heap.into_sorted_vec() + .into_iter() + .map(|(FloatOrd(dist), index)| (dist, &database[index][..], index)) + .collect() +} + +pub fn solve_challenge(challenge: &Challenge) -> anyhow::Result> { + let query_count = challenge.query_vectors.len(); + + let max_fuel = 2000000000.0; + let base_fuel = 760000000.0; + let alpha = 1700.0 * challenge.difficulty.num_queries as f64; + + let m = ((max_fuel - base_fuel) / alpha) as usize; + let n = (m as f32 * 1.2) as usize; + let r = n - m; + + let closest_vectors = filter_relevant_vectors( + &challenge.vector_database, + &challenge.query_vectors, + n, + ); + + let (m_slice, r_slice) = closest_vectors.split_at(m); + let m_vectors: Vec<_> = m_slice.to_vec(); + let r_vectors: Vec<_> = r_slice.to_vec(); + + let mut kd_tree_vectors: Vec<(&[f32], usize)> = m_vectors.iter().map(|&(_, v, i)| (v, i)).collect(); + let kd_tree = build_kd_tree(&mut kd_tree_vectors); + + let mut best_indexes = Vec::with_capacity(query_count); + let mut distances = Vec::with_capacity(query_count); + + for query in &challenge.query_vectors { + let mut best = (std::f32::MAX, None); + nearest_neighbor_search(&kd_tree, query, &mut best); + + distances.push(best.0); + best_indexes.push(best.1.unwrap_or(0)); + } + + let brute_force_count = (query_count as f32 * 0.1) as usize; + let mut distance_indices: Vec<_> = distances.iter().enumerate().collect(); + distance_indices.sort_unstable_by(|a, b| b.1.partial_cmp(a.1).unwrap()); + let high_distance_indices: Vec<_> = distance_indices.into_iter() + .take(brute_force_count) + .map(|(index, _)| index) + .collect(); + + for &query_index in &high_distance_indices { + let query = &challenge.query_vectors[query_index]; + let mut best = (distances[query_index], best_indexes[query_index]); + + for &(_, vec, index) in &r_vectors { + let dist = squared_euclidean_distance(query, vec); + if dist < best.0 { + best = (dist, index); + } + } + + best_indexes[query_index] = best.1; + } + + Ok(Some(Solution { + indexes: best_indexes, + })) +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; diff --git a/tig-algorithms/src/vector_search/invector_hybrid/mod.rs b/tig-algorithms/src/vector_search/invector_hybrid/mod.rs new file mode 100644 index 00000000..fcec9672 --- /dev/null +++ b/tig-algorithms/src/vector_search/invector_hybrid/mod.rs @@ -0,0 +1,4 @@ +mod benchmarker_outbound; +pub use benchmarker_outbound::solve_challenge; +#[cfg(feature = "cuda")] +pub use benchmarker_outbound::{cuda_solve_challenge, KERNEL}; \ No newline at end of file diff --git a/tig-algorithms/src/vector_search/invector_hybrid/open_data.rs b/tig-algorithms/src/vector_search/invector_hybrid/open_data.rs new file mode 100644 index 00000000..e5284668 --- /dev/null +++ b/tig-algorithms/src/vector_search/invector_hybrid/open_data.rs @@ -0,0 +1,395 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Open Data License v1.0 or (at your option) any later version +(the "License"); you may not use this file except in compliance with the License. +You may obtain a copy of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + + +use anyhow::Ok; +use tig_challenges::vector_search::*; +use std::cmp::Ordering; +use std::collections::BinaryHeap; + +struct KDNode<'a> { + point: &'a [f32], + left: Option>>, + right: Option>>, + index: usize, +} + +impl<'a> KDNode<'a> { + fn new(point: &'a [f32], index: usize) -> Self { + KDNode { + point, + left: None, + right: None, + index, + } + } +} +fn quickselect_by(arr: &mut [(&[f32], usize)], k: usize, compare: &F) +where + F: Fn(&(&[f32], usize), &(&[f32], usize)) -> Ordering, +{ + if arr.len() <= 1 { + return; + } + + let pivot_index = partition(arr, compare); + if k < pivot_index { + quickselect_by(&mut arr[..pivot_index], k, compare); + } else if k > pivot_index { + quickselect_by(&mut arr[pivot_index + 1..], k - pivot_index - 1, compare); + } +} + +fn partition(arr: &mut [(&[f32], usize)], compare: &F) -> usize +where + F: Fn(&(&[f32], usize), &(&[f32], usize)) -> Ordering, +{ + let pivot_index = arr.len() >> 1; + arr.swap(pivot_index, arr.len() - 1); + + let mut store_index = 0; + for i in 0..arr.len() - 1 { + if compare(&arr[i], &arr[arr.len() - 1]) == Ordering::Less { + arr.swap(i, store_index); + store_index += 1; + } + } + arr.swap(store_index, arr.len() - 1); + store_index +} + +fn build_kd_tree<'a>(points: &mut [(&'a [f32], usize)]) -> Option>> { + if points.is_empty() { + return None; + } + + const NUM_DIMENSIONS: usize = 250; + let mut stack: Vec<(usize, usize, usize, Option<*mut KDNode<'a>>, bool)> = Vec::new(); + let mut root: Option>> = None; + + stack.push((0, points.len(), 0, None, false)); + + while let Some((start, end, depth, parent_ptr, is_left)) = stack.pop() { + if start >= end { + continue; + } + + let axis = depth % NUM_DIMENSIONS; + let median = (start + end) / 2; + quickselect_by(&mut points[start..end], median - start, &|a, b| { + a.0[axis].partial_cmp(&b.0[axis]).unwrap() + }); + + let (median_point, median_index) = points[median]; + let mut new_node = Box::new(KDNode::new(median_point, median_index)); + let new_node_ptr: *mut KDNode = &mut *new_node; + + if let Some(parent_ptr) = parent_ptr { + unsafe { + if is_left { + (*parent_ptr).left = Some(new_node); + } else { + (*parent_ptr).right = Some(new_node); + } + } + } else { + root = Some(new_node); + } + + stack.push((median + 1, end, depth + 1, Some(new_node_ptr), false)); + stack.push((start, median, depth + 1, Some(new_node_ptr), true)); + } + + root +} + +#[inline(always)] +fn squared_euclidean_distance(a: &[f32], b: &[f32]) -> f32 { + let mut sum = 0.0; + let mut i = 0; + let len = a.len(); + + if a.len() != b.len() || a.len() < 8 { + return f32::MAX; + } + + while i + 7 < len { + unsafe { + let diff0 = *a.get_unchecked(i) - *b.get_unchecked(i); + let diff1 = *a.get_unchecked(i + 1) - *b.get_unchecked(i + 1); + let diff2 = *a.get_unchecked(i + 2) - *b.get_unchecked(i + 2); + let diff3 = *a.get_unchecked(i + 3) - *b.get_unchecked(i + 3); + let diff4 = *a.get_unchecked(i + 4) - *b.get_unchecked(i + 4); + let diff5 = *a.get_unchecked(i + 5) - *b.get_unchecked(i + 5); + let diff6 = *a.get_unchecked(i + 6) - *b.get_unchecked(i + 6); + let diff7 = *a.get_unchecked(i + 7) - *b.get_unchecked(i + 7); + + sum += diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3 + + diff4 * diff4 + diff5 * diff5 + diff6 * diff6 + diff7 * diff7; + } + + i += 8; + } + + while i < len { + unsafe { + let diff = *a.get_unchecked(i) - *b.get_unchecked(i); + sum += diff * diff; + } + i += 1; + } + sum +} + +#[inline(always)] +fn early_stopping_distance(a: &[f32], b: &[f32], current_min: f32) -> f32 { + let mut sum = 0.0; + let mut i = 0; + let len = a.len(); + + if a.len() != b.len() || a.len() < 8 { + return f32::MAX; + } + + while i + 7 < len { + unsafe { + let diff0 = *a.get_unchecked(i) - *b.get_unchecked(i); + let diff1 = *a.get_unchecked(i + 1) - *b.get_unchecked(i + 1); + let diff2 = *a.get_unchecked(i + 2) - *b.get_unchecked(i + 2); + let diff3 = *a.get_unchecked(i + 3) - *b.get_unchecked(i + 3); + let diff4 = *a.get_unchecked(i + 4) - *b.get_unchecked(i + 4); + let diff5 = *a.get_unchecked(i + 5) - *b.get_unchecked(i + 5); + let diff6 = *a.get_unchecked(i + 6) - *b.get_unchecked(i + 6); + let diff7 = *a.get_unchecked(i + 7) - *b.get_unchecked(i + 7); + + sum += diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3 + + diff4 * diff4 + diff5 * diff5 + diff6 * diff6 + diff7 * diff7; + } + + if sum > current_min { + return f32::MAX; + } + + i += 8; + } + + while i < len { + unsafe { + let diff = *a.get_unchecked(i) - *b.get_unchecked(i); + sum += diff * diff; + } + i += 1; + } + sum +} + +fn nearest_neighbor_search<'a>( + root: &Option>>, + target: &[f32], + best: &mut (f32, Option), +) { + let num_dimensions = target.len(); + let mut stack = Vec::with_capacity(64); + + if let Some(node) = root { + stack.push((node.as_ref(), 0)); + } + + while let Some((node, depth)) = stack.pop() { + let axis = depth % num_dimensions; + let dist = early_stopping_distance(&node.point, target, best.0); + + if dist < best.0 { + best.0 = dist; + best.1 = Some(node.index); + } + + let diff = target[axis] - node.point[axis]; + let sqr_diff = diff * diff; + + let (nearer, farther) = if diff < 0.0 { + (&node.left, &node.right) + } else { + (&node.right, &node.left) + }; + + if let Some(nearer_node) = nearer { + stack.push((nearer_node.as_ref(), depth + 1)); + } + + if sqr_diff < best.0 { + if let Some(farther_node) = farther { + stack.push((farther_node.as_ref(), depth + 1)); + } + } + } +} +fn calculate_mean_vector(vectors: &[&[f32]]) -> Vec { + let num_vectors = vectors.len(); + let num_dimensions = 250; + + let mut mean_vector = vec![0.0f64; num_dimensions]; + + for vector in vectors { + for i in 0..num_dimensions { + mean_vector[i] += vector[i] as f64; + } + } + for i in 0..num_dimensions { + mean_vector[i] /= num_vectors as f64; + } + mean_vector.into_iter().map(|x| x as f32).collect() +} + +#[derive(Debug)] +struct FloatOrd(f32); + +impl PartialEq for FloatOrd { + fn eq(&self, other: &Self) -> bool { + self.0 == other.0 + } +} + +impl Eq for FloatOrd {} + +impl PartialOrd for FloatOrd { + fn partial_cmp(&self, other: &Self) -> Option { + self.0.partial_cmp(&other.0) + } +} + +impl Ord for FloatOrd { + fn cmp(&self, other: &Self) -> Ordering { + + self.partial_cmp(other).unwrap_or(Ordering::Equal) + } +} + +fn filter_relevant_vectors<'a>( + database: &'a [Vec], + query_vectors: &[Vec], + k: usize, +) -> Vec<(f32, &'a [f32], usize)> { + let query_refs: Vec<&[f32]> = query_vectors.iter().map(|v| &v[..]).collect(); + let mean_query_vector = calculate_mean_vector(&query_refs); + + let mut heap: BinaryHeap<(FloatOrd, usize)> = BinaryHeap::with_capacity(k); + + for (index, vector) in database.iter().enumerate() { + if heap.len() < k + { + let dist = squared_euclidean_distance(&mean_query_vector, vector); + let ord_dist = FloatOrd(dist); + + heap.push((ord_dist, index)); + } else if let Some(&(FloatOrd(top_dist), _)) = heap.peek() + { + let dist = early_stopping_distance(&mean_query_vector, vector, top_dist); + let ord_dist = FloatOrd(dist); + if dist < top_dist { + heap.pop(); + heap.push((ord_dist, index)); + } + } + } + heap.into_sorted_vec() + .into_iter() + .map(|(FloatOrd(dist), index)| (dist, &database[index][..], index)) + .collect() +} + +pub fn solve_challenge(challenge: &Challenge) -> anyhow::Result> { + let query_count = challenge.query_vectors.len(); + + let max_fuel = 2000000000.0; + let base_fuel = 760000000.0; + let alpha = 1700.0 * challenge.difficulty.num_queries as f64; + + let m = ((max_fuel - base_fuel) / alpha) as usize; + let n = (m as f32 * 1.2) as usize; + let r = n - m; + + let closest_vectors = filter_relevant_vectors( + &challenge.vector_database, + &challenge.query_vectors, + n, + ); + + let (m_slice, r_slice) = closest_vectors.split_at(m); + let m_vectors: Vec<_> = m_slice.to_vec(); + let r_vectors: Vec<_> = r_slice.to_vec(); + + let mut kd_tree_vectors: Vec<(&[f32], usize)> = m_vectors.iter().map(|&(_, v, i)| (v, i)).collect(); + let kd_tree = build_kd_tree(&mut kd_tree_vectors); + + let mut best_indexes = Vec::with_capacity(query_count); + let mut distances = Vec::with_capacity(query_count); + + for query in &challenge.query_vectors { + let mut best = (std::f32::MAX, None); + nearest_neighbor_search(&kd_tree, query, &mut best); + + distances.push(best.0); + best_indexes.push(best.1.unwrap_or(0)); + } + + let brute_force_count = (query_count as f32 * 0.1) as usize; + let mut distance_indices: Vec<_> = distances.iter().enumerate().collect(); + distance_indices.sort_unstable_by(|a, b| b.1.partial_cmp(a.1).unwrap()); + let high_distance_indices: Vec<_> = distance_indices.into_iter() + .take(brute_force_count) + .map(|(index, _)| index) + .collect(); + + for &query_index in &high_distance_indices { + let query = &challenge.query_vectors[query_index]; + let mut best = (distances[query_index], best_indexes[query_index]); + + for &(_, vec, index) in &r_vectors { + let dist = squared_euclidean_distance(query, vec); + if dist < best.0 { + best = (dist, index); + } + } + + best_indexes[query_index] = best.1; + } + + Ok(Some(Solution { + indexes: best_indexes, + })) +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; diff --git a/tig-algorithms/src/vector_search/mod.rs b/tig-algorithms/src/vector_search/mod.rs index bf245d39..cee7ac52 100644 --- a/tig-algorithms/src/vector_search/mod.rs +++ b/tig-algorithms/src/vector_search/mod.rs @@ -83,7 +83,8 @@ pub use invector as c004_a034; // c004_a041 -// c004_a042 +pub mod invector_hybrid; +pub use invector_hybrid as c004_a042; // c004_a043 diff --git a/tig-algorithms/src/vehicle_routing/advanced_routing/benchmarker_outbound.rs b/tig-algorithms/src/vehicle_routing/advanced_routing/benchmarker_outbound.rs new file mode 100644 index 00000000..8d7001fb --- /dev/null +++ b/tig-algorithms/src/vehicle_routing/advanced_routing/benchmarker_outbound.rs @@ -0,0 +1,234 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Benchmarker Outbound Game License v1.0 (the "License"); you +may not use this file except in compliance with the License. You may obtain a copy +of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + +use rand::{rngs::{SmallRng, StdRng}, Rng, SeedableRng}; +use tig_challenges::vehicle_routing::{Challenge, Solution}; + +pub fn solve_challenge(challenge: &Challenge) -> anyhow::Result> { + let mut best_solution: Option = None; + let mut best_cost = std::i32::MAX; + + const INITIAL_TEMPERATURE: f32 = 2.0; + const COOLING_RATE: f32 = 0.995; + const ITERATIONS_PER_TEMPERATURE: usize = 2; + + let num_nodes = challenge.difficulty.num_nodes; + + let mut current_params = vec![1.0; num_nodes]; + let mut savings_list = create_initial_savings_list(challenge); + recompute_and_sort_savings(&mut savings_list, ¤t_params, challenge); + + let mut current_solution = create_solution(challenge, ¤t_params, &savings_list); + let mut current_cost = calculate_solution_cost(¤t_solution, &challenge.distance_matrix); + + if current_cost <= challenge.max_total_distance { + return Ok(Some(current_solution)); + } + + if (current_cost as f32 * 0.96) > challenge.max_total_distance as f32 { + return Ok(None); + } + + let mut temperature = INITIAL_TEMPERATURE; + let mut rng = StdRng::seed_from_u64(u64::from_le_bytes(challenge.seed[..8].try_into().unwrap())); + + while temperature > 1.0 { + for _ in 0..ITERATIONS_PER_TEMPERATURE { + let neighbor_params = generate_neighbor(¤t_params, &mut rng); + recompute_and_sort_savings(&mut savings_list, &neighbor_params, challenge); + let mut neighbor_solution = create_solution(challenge, &neighbor_params, &savings_list); + apply_local_search_until_no_improvement(&mut neighbor_solution, &challenge.distance_matrix); + let neighbor_cost = calculate_solution_cost(&neighbor_solution, &challenge.distance_matrix); + + let delta = neighbor_cost as f32 - current_cost as f32; + if delta < 0.0 || rng.gen::() < (-delta / temperature).exp() { + current_params = neighbor_params; + current_cost = neighbor_cost; + current_solution = neighbor_solution; + + if current_cost < best_cost { + best_cost = current_cost; + best_solution = Some(Solution { + routes: current_solution.routes.clone(), + }); + } + } + if best_cost <= challenge.max_total_distance { + return Ok(best_solution); + } + } + + temperature *= COOLING_RATE; + } + + Ok(best_solution) +} + +#[inline] +fn create_initial_savings_list(challenge: &Challenge) -> Vec<(f32, u8, u8)> { + let num_nodes = challenge.difficulty.num_nodes; + let capacity = ((num_nodes - 1) * (num_nodes - 2)) / 2; + let mut savings = Vec::with_capacity(capacity); + for i in 1..num_nodes { + for j in (i + 1)..num_nodes { + savings.push((0.0, i as u8, j as u8)); + } + } + savings +} + +#[inline] +fn recompute_and_sort_savings(savings_list: &mut [(f32, u8, u8)], params: &[f32], challenge: &Challenge) { + let distance_matrix = &challenge.distance_matrix; + + let mut zero_len = 0; + for (score, i, j) in savings_list.iter_mut() { + let i = *i as usize; + let j = *j as usize; + *score = params[i] * distance_matrix[0][i] as f32 + + params[j] * distance_matrix[j][0] as f32 - + params[i] * params[j] * distance_matrix[i][j] as f32; + } + + savings_list.sort_unstable_by(|a, b| b.0.partial_cmp(&a.0).unwrap()); +} + +#[inline] +fn generate_neighbor(current: &[f32], rng: &mut R) -> Vec { + current.iter().map(|¶m| { + let delta = rng.gen_range(-0.1..=0.1); + (param + delta).clamp(0.0, 2.0) + }).collect() +} + +#[inline] +fn apply_local_search_until_no_improvement(solution: &mut Solution, distance_matrix: &Vec>) { + let mut improved = true; + while improved { + improved = false; + for route in &mut solution.routes { + if two_opt(route, distance_matrix) { + improved = true; + } + } + } +} +#[inline] +fn two_opt(route: &mut Vec, distance_matrix: &Vec>) -> bool { + let n = route.len(); + let mut improved = false; + + for i in 1..n - 2 { + for j in i + 1..n - 1 { + let current_distance = distance_matrix[route[i - 1]][route[i]] + + distance_matrix[route[j]][route[j + 1]]; + let new_distance = distance_matrix[route[i - 1]][route[j]] + + distance_matrix[route[i]][route[j + 1]]; + + if new_distance < current_distance { + route[i..=j].reverse(); + improved = true; + } + } + } + + improved +} + +#[inline] +fn calculate_solution_cost(solution: &Solution, distance_matrix: &Vec>) -> i32 { + solution.routes.iter().map(|route| { + route.windows(2).map(|w| distance_matrix[w[0]][w[1]]).sum::() + }).sum() +} + +#[inline] +fn create_solution(challenge: &Challenge, params: &[f32], savings_list: &[(f32, u8, u8)]) -> Solution { + let distance_matrix = &challenge.distance_matrix; + let max_capacity = challenge.max_capacity; + let num_nodes = challenge.difficulty.num_nodes; + let demands = &challenge.demands; + + let mut routes = vec![None; num_nodes]; + for i in 1..num_nodes { + routes[i] = Some(vec![i]); + } + let mut route_demands = demands.clone(); + + for &(_, i, j) in savings_list { + let (i, j) = (i as usize, j as usize); + if let (Some(left_route), Some(right_route)) = (routes[i].as_ref(), routes[j].as_ref()) { + let (left_start, left_end) = (*left_route.first().unwrap(), *left_route.last().unwrap()); + let (right_start, right_end) = (*right_route.first().unwrap(), *right_route.last().unwrap()); + + if left_start == right_start || route_demands[left_start] + route_demands[right_start] > max_capacity { + continue; + } + + let mut new_route = routes[i].take().unwrap(); + let mut right_route = routes[j].take().unwrap(); + + if left_start == i { new_route.reverse(); } + if right_end == j { right_route.reverse(); } + + new_route.extend(right_route); + + let combined_demand = route_demands[left_start] + route_demands[right_start]; + let new_start = new_route[0]; + let new_end = *new_route.last().unwrap(); + + route_demands[new_start] = combined_demand; + route_demands[new_end] = combined_demand; + + routes[new_start] = Some(new_route.clone()); + routes[new_end] = Some(new_route); + } + } + + Solution { + routes: routes + .into_iter() + .enumerate() + .filter_map(|(i, route)| route.filter(|r| r[0] == i)) + .map(|mut route| { + route.insert(0, 0); + route.push(0); + route + }) + .collect(), + } +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; \ No newline at end of file diff --git a/tig-algorithms/src/vehicle_routing/advanced_routing/commercial.rs b/tig-algorithms/src/vehicle_routing/advanced_routing/commercial.rs new file mode 100644 index 00000000..bb1110f1 --- /dev/null +++ b/tig-algorithms/src/vehicle_routing/advanced_routing/commercial.rs @@ -0,0 +1,234 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Commercial License v1.0 (the "License"); you +may not use this file except in compliance with the License. You may obtain a copy +of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + +use rand::{rngs::{SmallRng, StdRng}, Rng, SeedableRng}; +use tig_challenges::vehicle_routing::{Challenge, Solution}; + +pub fn solve_challenge(challenge: &Challenge) -> anyhow::Result> { + let mut best_solution: Option = None; + let mut best_cost = std::i32::MAX; + + const INITIAL_TEMPERATURE: f32 = 2.0; + const COOLING_RATE: f32 = 0.995; + const ITERATIONS_PER_TEMPERATURE: usize = 2; + + let num_nodes = challenge.difficulty.num_nodes; + + let mut current_params = vec![1.0; num_nodes]; + let mut savings_list = create_initial_savings_list(challenge); + recompute_and_sort_savings(&mut savings_list, ¤t_params, challenge); + + let mut current_solution = create_solution(challenge, ¤t_params, &savings_list); + let mut current_cost = calculate_solution_cost(¤t_solution, &challenge.distance_matrix); + + if current_cost <= challenge.max_total_distance { + return Ok(Some(current_solution)); + } + + if (current_cost as f32 * 0.96) > challenge.max_total_distance as f32 { + return Ok(None); + } + + let mut temperature = INITIAL_TEMPERATURE; + let mut rng = StdRng::seed_from_u64(u64::from_le_bytes(challenge.seed[..8].try_into().unwrap())); + + while temperature > 1.0 { + for _ in 0..ITERATIONS_PER_TEMPERATURE { + let neighbor_params = generate_neighbor(¤t_params, &mut rng); + recompute_and_sort_savings(&mut savings_list, &neighbor_params, challenge); + let mut neighbor_solution = create_solution(challenge, &neighbor_params, &savings_list); + apply_local_search_until_no_improvement(&mut neighbor_solution, &challenge.distance_matrix); + let neighbor_cost = calculate_solution_cost(&neighbor_solution, &challenge.distance_matrix); + + let delta = neighbor_cost as f32 - current_cost as f32; + if delta < 0.0 || rng.gen::() < (-delta / temperature).exp() { + current_params = neighbor_params; + current_cost = neighbor_cost; + current_solution = neighbor_solution; + + if current_cost < best_cost { + best_cost = current_cost; + best_solution = Some(Solution { + routes: current_solution.routes.clone(), + }); + } + } + if best_cost <= challenge.max_total_distance { + return Ok(best_solution); + } + } + + temperature *= COOLING_RATE; + } + + Ok(best_solution) +} + +#[inline] +fn create_initial_savings_list(challenge: &Challenge) -> Vec<(f32, u8, u8)> { + let num_nodes = challenge.difficulty.num_nodes; + let capacity = ((num_nodes - 1) * (num_nodes - 2)) / 2; + let mut savings = Vec::with_capacity(capacity); + for i in 1..num_nodes { + for j in (i + 1)..num_nodes { + savings.push((0.0, i as u8, j as u8)); + } + } + savings +} + +#[inline] +fn recompute_and_sort_savings(savings_list: &mut [(f32, u8, u8)], params: &[f32], challenge: &Challenge) { + let distance_matrix = &challenge.distance_matrix; + + let mut zero_len = 0; + for (score, i, j) in savings_list.iter_mut() { + let i = *i as usize; + let j = *j as usize; + *score = params[i] * distance_matrix[0][i] as f32 + + params[j] * distance_matrix[j][0] as f32 - + params[i] * params[j] * distance_matrix[i][j] as f32; + } + + savings_list.sort_unstable_by(|a, b| b.0.partial_cmp(&a.0).unwrap()); +} + +#[inline] +fn generate_neighbor(current: &[f32], rng: &mut R) -> Vec { + current.iter().map(|¶m| { + let delta = rng.gen_range(-0.1..=0.1); + (param + delta).clamp(0.0, 2.0) + }).collect() +} + +#[inline] +fn apply_local_search_until_no_improvement(solution: &mut Solution, distance_matrix: &Vec>) { + let mut improved = true; + while improved { + improved = false; + for route in &mut solution.routes { + if two_opt(route, distance_matrix) { + improved = true; + } + } + } +} +#[inline] +fn two_opt(route: &mut Vec, distance_matrix: &Vec>) -> bool { + let n = route.len(); + let mut improved = false; + + for i in 1..n - 2 { + for j in i + 1..n - 1 { + let current_distance = distance_matrix[route[i - 1]][route[i]] + + distance_matrix[route[j]][route[j + 1]]; + let new_distance = distance_matrix[route[i - 1]][route[j]] + + distance_matrix[route[i]][route[j + 1]]; + + if new_distance < current_distance { + route[i..=j].reverse(); + improved = true; + } + } + } + + improved +} + +#[inline] +fn calculate_solution_cost(solution: &Solution, distance_matrix: &Vec>) -> i32 { + solution.routes.iter().map(|route| { + route.windows(2).map(|w| distance_matrix[w[0]][w[1]]).sum::() + }).sum() +} + +#[inline] +fn create_solution(challenge: &Challenge, params: &[f32], savings_list: &[(f32, u8, u8)]) -> Solution { + let distance_matrix = &challenge.distance_matrix; + let max_capacity = challenge.max_capacity; + let num_nodes = challenge.difficulty.num_nodes; + let demands = &challenge.demands; + + let mut routes = vec![None; num_nodes]; + for i in 1..num_nodes { + routes[i] = Some(vec![i]); + } + let mut route_demands = demands.clone(); + + for &(_, i, j) in savings_list { + let (i, j) = (i as usize, j as usize); + if let (Some(left_route), Some(right_route)) = (routes[i].as_ref(), routes[j].as_ref()) { + let (left_start, left_end) = (*left_route.first().unwrap(), *left_route.last().unwrap()); + let (right_start, right_end) = (*right_route.first().unwrap(), *right_route.last().unwrap()); + + if left_start == right_start || route_demands[left_start] + route_demands[right_start] > max_capacity { + continue; + } + + let mut new_route = routes[i].take().unwrap(); + let mut right_route = routes[j].take().unwrap(); + + if left_start == i { new_route.reverse(); } + if right_end == j { right_route.reverse(); } + + new_route.extend(right_route); + + let combined_demand = route_demands[left_start] + route_demands[right_start]; + let new_start = new_route[0]; + let new_end = *new_route.last().unwrap(); + + route_demands[new_start] = combined_demand; + route_demands[new_end] = combined_demand; + + routes[new_start] = Some(new_route.clone()); + routes[new_end] = Some(new_route); + } + } + + Solution { + routes: routes + .into_iter() + .enumerate() + .filter_map(|(i, route)| route.filter(|r| r[0] == i)) + .map(|mut route| { + route.insert(0, 0); + route.push(0); + route + }) + .collect(), + } +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; \ No newline at end of file diff --git a/tig-algorithms/src/vehicle_routing/advanced_routing/inbound.rs b/tig-algorithms/src/vehicle_routing/advanced_routing/inbound.rs new file mode 100644 index 00000000..2b8f8e13 --- /dev/null +++ b/tig-algorithms/src/vehicle_routing/advanced_routing/inbound.rs @@ -0,0 +1,234 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Inbound Game License v1.0 or (at your option) any later +version (the "License"); you may not use this file except in compliance with the +License. You may obtain a copy of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + +use rand::{rngs::{SmallRng, StdRng}, Rng, SeedableRng}; +use tig_challenges::vehicle_routing::{Challenge, Solution}; + +pub fn solve_challenge(challenge: &Challenge) -> anyhow::Result> { + let mut best_solution: Option = None; + let mut best_cost = std::i32::MAX; + + const INITIAL_TEMPERATURE: f32 = 2.0; + const COOLING_RATE: f32 = 0.995; + const ITERATIONS_PER_TEMPERATURE: usize = 2; + + let num_nodes = challenge.difficulty.num_nodes; + + let mut current_params = vec![1.0; num_nodes]; + let mut savings_list = create_initial_savings_list(challenge); + recompute_and_sort_savings(&mut savings_list, ¤t_params, challenge); + + let mut current_solution = create_solution(challenge, ¤t_params, &savings_list); + let mut current_cost = calculate_solution_cost(¤t_solution, &challenge.distance_matrix); + + if current_cost <= challenge.max_total_distance { + return Ok(Some(current_solution)); + } + + if (current_cost as f32 * 0.96) > challenge.max_total_distance as f32 { + return Ok(None); + } + + let mut temperature = INITIAL_TEMPERATURE; + let mut rng = StdRng::seed_from_u64(u64::from_le_bytes(challenge.seed[..8].try_into().unwrap())); + + while temperature > 1.0 { + for _ in 0..ITERATIONS_PER_TEMPERATURE { + let neighbor_params = generate_neighbor(¤t_params, &mut rng); + recompute_and_sort_savings(&mut savings_list, &neighbor_params, challenge); + let mut neighbor_solution = create_solution(challenge, &neighbor_params, &savings_list); + apply_local_search_until_no_improvement(&mut neighbor_solution, &challenge.distance_matrix); + let neighbor_cost = calculate_solution_cost(&neighbor_solution, &challenge.distance_matrix); + + let delta = neighbor_cost as f32 - current_cost as f32; + if delta < 0.0 || rng.gen::() < (-delta / temperature).exp() { + current_params = neighbor_params; + current_cost = neighbor_cost; + current_solution = neighbor_solution; + + if current_cost < best_cost { + best_cost = current_cost; + best_solution = Some(Solution { + routes: current_solution.routes.clone(), + }); + } + } + if best_cost <= challenge.max_total_distance { + return Ok(best_solution); + } + } + + temperature *= COOLING_RATE; + } + + Ok(best_solution) +} + +#[inline] +fn create_initial_savings_list(challenge: &Challenge) -> Vec<(f32, u8, u8)> { + let num_nodes = challenge.difficulty.num_nodes; + let capacity = ((num_nodes - 1) * (num_nodes - 2)) / 2; + let mut savings = Vec::with_capacity(capacity); + for i in 1..num_nodes { + for j in (i + 1)..num_nodes { + savings.push((0.0, i as u8, j as u8)); + } + } + savings +} + +#[inline] +fn recompute_and_sort_savings(savings_list: &mut [(f32, u8, u8)], params: &[f32], challenge: &Challenge) { + let distance_matrix = &challenge.distance_matrix; + + let mut zero_len = 0; + for (score, i, j) in savings_list.iter_mut() { + let i = *i as usize; + let j = *j as usize; + *score = params[i] * distance_matrix[0][i] as f32 + + params[j] * distance_matrix[j][0] as f32 - + params[i] * params[j] * distance_matrix[i][j] as f32; + } + + savings_list.sort_unstable_by(|a, b| b.0.partial_cmp(&a.0).unwrap()); +} + +#[inline] +fn generate_neighbor(current: &[f32], rng: &mut R) -> Vec { + current.iter().map(|¶m| { + let delta = rng.gen_range(-0.1..=0.1); + (param + delta).clamp(0.0, 2.0) + }).collect() +} + +#[inline] +fn apply_local_search_until_no_improvement(solution: &mut Solution, distance_matrix: &Vec>) { + let mut improved = true; + while improved { + improved = false; + for route in &mut solution.routes { + if two_opt(route, distance_matrix) { + improved = true; + } + } + } +} +#[inline] +fn two_opt(route: &mut Vec, distance_matrix: &Vec>) -> bool { + let n = route.len(); + let mut improved = false; + + for i in 1..n - 2 { + for j in i + 1..n - 1 { + let current_distance = distance_matrix[route[i - 1]][route[i]] + + distance_matrix[route[j]][route[j + 1]]; + let new_distance = distance_matrix[route[i - 1]][route[j]] + + distance_matrix[route[i]][route[j + 1]]; + + if new_distance < current_distance { + route[i..=j].reverse(); + improved = true; + } + } + } + + improved +} + +#[inline] +fn calculate_solution_cost(solution: &Solution, distance_matrix: &Vec>) -> i32 { + solution.routes.iter().map(|route| { + route.windows(2).map(|w| distance_matrix[w[0]][w[1]]).sum::() + }).sum() +} + +#[inline] +fn create_solution(challenge: &Challenge, params: &[f32], savings_list: &[(f32, u8, u8)]) -> Solution { + let distance_matrix = &challenge.distance_matrix; + let max_capacity = challenge.max_capacity; + let num_nodes = challenge.difficulty.num_nodes; + let demands = &challenge.demands; + + let mut routes = vec![None; num_nodes]; + for i in 1..num_nodes { + routes[i] = Some(vec![i]); + } + let mut route_demands = demands.clone(); + + for &(_, i, j) in savings_list { + let (i, j) = (i as usize, j as usize); + if let (Some(left_route), Some(right_route)) = (routes[i].as_ref(), routes[j].as_ref()) { + let (left_start, left_end) = (*left_route.first().unwrap(), *left_route.last().unwrap()); + let (right_start, right_end) = (*right_route.first().unwrap(), *right_route.last().unwrap()); + + if left_start == right_start || route_demands[left_start] + route_demands[right_start] > max_capacity { + continue; + } + + let mut new_route = routes[i].take().unwrap(); + let mut right_route = routes[j].take().unwrap(); + + if left_start == i { new_route.reverse(); } + if right_end == j { right_route.reverse(); } + + new_route.extend(right_route); + + let combined_demand = route_demands[left_start] + route_demands[right_start]; + let new_start = new_route[0]; + let new_end = *new_route.last().unwrap(); + + route_demands[new_start] = combined_demand; + route_demands[new_end] = combined_demand; + + routes[new_start] = Some(new_route.clone()); + routes[new_end] = Some(new_route); + } + } + + Solution { + routes: routes + .into_iter() + .enumerate() + .filter_map(|(i, route)| route.filter(|r| r[0] == i)) + .map(|mut route| { + route.insert(0, 0); + route.push(0); + route + }) + .collect(), + } +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; \ No newline at end of file diff --git a/tig-algorithms/src/vehicle_routing/advanced_routing/innovator_outbound.rs b/tig-algorithms/src/vehicle_routing/advanced_routing/innovator_outbound.rs new file mode 100644 index 00000000..e768ec0e --- /dev/null +++ b/tig-algorithms/src/vehicle_routing/advanced_routing/innovator_outbound.rs @@ -0,0 +1,234 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Innovator Outbound Game License v1.0 (the "License"); you +may not use this file except in compliance with the License. You may obtain a copy +of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + +use rand::{rngs::{SmallRng, StdRng}, Rng, SeedableRng}; +use tig_challenges::vehicle_routing::{Challenge, Solution}; + +pub fn solve_challenge(challenge: &Challenge) -> anyhow::Result> { + let mut best_solution: Option = None; + let mut best_cost = std::i32::MAX; + + const INITIAL_TEMPERATURE: f32 = 2.0; + const COOLING_RATE: f32 = 0.995; + const ITERATIONS_PER_TEMPERATURE: usize = 2; + + let num_nodes = challenge.difficulty.num_nodes; + + let mut current_params = vec![1.0; num_nodes]; + let mut savings_list = create_initial_savings_list(challenge); + recompute_and_sort_savings(&mut savings_list, ¤t_params, challenge); + + let mut current_solution = create_solution(challenge, ¤t_params, &savings_list); + let mut current_cost = calculate_solution_cost(¤t_solution, &challenge.distance_matrix); + + if current_cost <= challenge.max_total_distance { + return Ok(Some(current_solution)); + } + + if (current_cost as f32 * 0.96) > challenge.max_total_distance as f32 { + return Ok(None); + } + + let mut temperature = INITIAL_TEMPERATURE; + let mut rng = StdRng::seed_from_u64(u64::from_le_bytes(challenge.seed[..8].try_into().unwrap())); + + while temperature > 1.0 { + for _ in 0..ITERATIONS_PER_TEMPERATURE { + let neighbor_params = generate_neighbor(¤t_params, &mut rng); + recompute_and_sort_savings(&mut savings_list, &neighbor_params, challenge); + let mut neighbor_solution = create_solution(challenge, &neighbor_params, &savings_list); + apply_local_search_until_no_improvement(&mut neighbor_solution, &challenge.distance_matrix); + let neighbor_cost = calculate_solution_cost(&neighbor_solution, &challenge.distance_matrix); + + let delta = neighbor_cost as f32 - current_cost as f32; + if delta < 0.0 || rng.gen::() < (-delta / temperature).exp() { + current_params = neighbor_params; + current_cost = neighbor_cost; + current_solution = neighbor_solution; + + if current_cost < best_cost { + best_cost = current_cost; + best_solution = Some(Solution { + routes: current_solution.routes.clone(), + }); + } + } + if best_cost <= challenge.max_total_distance { + return Ok(best_solution); + } + } + + temperature *= COOLING_RATE; + } + + Ok(best_solution) +} + +#[inline] +fn create_initial_savings_list(challenge: &Challenge) -> Vec<(f32, u8, u8)> { + let num_nodes = challenge.difficulty.num_nodes; + let capacity = ((num_nodes - 1) * (num_nodes - 2)) / 2; + let mut savings = Vec::with_capacity(capacity); + for i in 1..num_nodes { + for j in (i + 1)..num_nodes { + savings.push((0.0, i as u8, j as u8)); + } + } + savings +} + +#[inline] +fn recompute_and_sort_savings(savings_list: &mut [(f32, u8, u8)], params: &[f32], challenge: &Challenge) { + let distance_matrix = &challenge.distance_matrix; + + let mut zero_len = 0; + for (score, i, j) in savings_list.iter_mut() { + let i = *i as usize; + let j = *j as usize; + *score = params[i] * distance_matrix[0][i] as f32 + + params[j] * distance_matrix[j][0] as f32 - + params[i] * params[j] * distance_matrix[i][j] as f32; + } + + savings_list.sort_unstable_by(|a, b| b.0.partial_cmp(&a.0).unwrap()); +} + +#[inline] +fn generate_neighbor(current: &[f32], rng: &mut R) -> Vec { + current.iter().map(|¶m| { + let delta = rng.gen_range(-0.1..=0.1); + (param + delta).clamp(0.0, 2.0) + }).collect() +} + +#[inline] +fn apply_local_search_until_no_improvement(solution: &mut Solution, distance_matrix: &Vec>) { + let mut improved = true; + while improved { + improved = false; + for route in &mut solution.routes { + if two_opt(route, distance_matrix) { + improved = true; + } + } + } +} +#[inline] +fn two_opt(route: &mut Vec, distance_matrix: &Vec>) -> bool { + let n = route.len(); + let mut improved = false; + + for i in 1..n - 2 { + for j in i + 1..n - 1 { + let current_distance = distance_matrix[route[i - 1]][route[i]] + + distance_matrix[route[j]][route[j + 1]]; + let new_distance = distance_matrix[route[i - 1]][route[j]] + + distance_matrix[route[i]][route[j + 1]]; + + if new_distance < current_distance { + route[i..=j].reverse(); + improved = true; + } + } + } + + improved +} + +#[inline] +fn calculate_solution_cost(solution: &Solution, distance_matrix: &Vec>) -> i32 { + solution.routes.iter().map(|route| { + route.windows(2).map(|w| distance_matrix[w[0]][w[1]]).sum::() + }).sum() +} + +#[inline] +fn create_solution(challenge: &Challenge, params: &[f32], savings_list: &[(f32, u8, u8)]) -> Solution { + let distance_matrix = &challenge.distance_matrix; + let max_capacity = challenge.max_capacity; + let num_nodes = challenge.difficulty.num_nodes; + let demands = &challenge.demands; + + let mut routes = vec![None; num_nodes]; + for i in 1..num_nodes { + routes[i] = Some(vec![i]); + } + let mut route_demands = demands.clone(); + + for &(_, i, j) in savings_list { + let (i, j) = (i as usize, j as usize); + if let (Some(left_route), Some(right_route)) = (routes[i].as_ref(), routes[j].as_ref()) { + let (left_start, left_end) = (*left_route.first().unwrap(), *left_route.last().unwrap()); + let (right_start, right_end) = (*right_route.first().unwrap(), *right_route.last().unwrap()); + + if left_start == right_start || route_demands[left_start] + route_demands[right_start] > max_capacity { + continue; + } + + let mut new_route = routes[i].take().unwrap(); + let mut right_route = routes[j].take().unwrap(); + + if left_start == i { new_route.reverse(); } + if right_end == j { right_route.reverse(); } + + new_route.extend(right_route); + + let combined_demand = route_demands[left_start] + route_demands[right_start]; + let new_start = new_route[0]; + let new_end = *new_route.last().unwrap(); + + route_demands[new_start] = combined_demand; + route_demands[new_end] = combined_demand; + + routes[new_start] = Some(new_route.clone()); + routes[new_end] = Some(new_route); + } + } + + Solution { + routes: routes + .into_iter() + .enumerate() + .filter_map(|(i, route)| route.filter(|r| r[0] == i)) + .map(|mut route| { + route.insert(0, 0); + route.push(0); + route + }) + .collect(), + } +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; \ No newline at end of file diff --git a/tig-algorithms/src/vehicle_routing/advanced_routing/mod.rs b/tig-algorithms/src/vehicle_routing/advanced_routing/mod.rs new file mode 100644 index 00000000..fcec9672 --- /dev/null +++ b/tig-algorithms/src/vehicle_routing/advanced_routing/mod.rs @@ -0,0 +1,4 @@ +mod benchmarker_outbound; +pub use benchmarker_outbound::solve_challenge; +#[cfg(feature = "cuda")] +pub use benchmarker_outbound::{cuda_solve_challenge, KERNEL}; \ No newline at end of file diff --git a/tig-algorithms/src/vehicle_routing/advanced_routing/open_data.rs b/tig-algorithms/src/vehicle_routing/advanced_routing/open_data.rs new file mode 100644 index 00000000..41bdafeb --- /dev/null +++ b/tig-algorithms/src/vehicle_routing/advanced_routing/open_data.rs @@ -0,0 +1,234 @@ +/*! +Copyright 2024 syebastian + +Licensed under the TIG Open Data License v1.0 or (at your option) any later version +(the "License"); you may not use this file except in compliance with the License. +You may obtain a copy of the License at + +https://github.com/tig-foundation/tig-monorepo/tree/main/docs/licenses + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the specific +language governing permissions and limitations under the License. +*/ + +use rand::{rngs::{SmallRng, StdRng}, Rng, SeedableRng}; +use tig_challenges::vehicle_routing::{Challenge, Solution}; + +pub fn solve_challenge(challenge: &Challenge) -> anyhow::Result> { + let mut best_solution: Option = None; + let mut best_cost = std::i32::MAX; + + const INITIAL_TEMPERATURE: f32 = 2.0; + const COOLING_RATE: f32 = 0.995; + const ITERATIONS_PER_TEMPERATURE: usize = 2; + + let num_nodes = challenge.difficulty.num_nodes; + + let mut current_params = vec![1.0; num_nodes]; + let mut savings_list = create_initial_savings_list(challenge); + recompute_and_sort_savings(&mut savings_list, ¤t_params, challenge); + + let mut current_solution = create_solution(challenge, ¤t_params, &savings_list); + let mut current_cost = calculate_solution_cost(¤t_solution, &challenge.distance_matrix); + + if current_cost <= challenge.max_total_distance { + return Ok(Some(current_solution)); + } + + if (current_cost as f32 * 0.96) > challenge.max_total_distance as f32 { + return Ok(None); + } + + let mut temperature = INITIAL_TEMPERATURE; + let mut rng = StdRng::seed_from_u64(u64::from_le_bytes(challenge.seed[..8].try_into().unwrap())); + + while temperature > 1.0 { + for _ in 0..ITERATIONS_PER_TEMPERATURE { + let neighbor_params = generate_neighbor(¤t_params, &mut rng); + recompute_and_sort_savings(&mut savings_list, &neighbor_params, challenge); + let mut neighbor_solution = create_solution(challenge, &neighbor_params, &savings_list); + apply_local_search_until_no_improvement(&mut neighbor_solution, &challenge.distance_matrix); + let neighbor_cost = calculate_solution_cost(&neighbor_solution, &challenge.distance_matrix); + + let delta = neighbor_cost as f32 - current_cost as f32; + if delta < 0.0 || rng.gen::() < (-delta / temperature).exp() { + current_params = neighbor_params; + current_cost = neighbor_cost; + current_solution = neighbor_solution; + + if current_cost < best_cost { + best_cost = current_cost; + best_solution = Some(Solution { + routes: current_solution.routes.clone(), + }); + } + } + if best_cost <= challenge.max_total_distance { + return Ok(best_solution); + } + } + + temperature *= COOLING_RATE; + } + + Ok(best_solution) +} + +#[inline] +fn create_initial_savings_list(challenge: &Challenge) -> Vec<(f32, u8, u8)> { + let num_nodes = challenge.difficulty.num_nodes; + let capacity = ((num_nodes - 1) * (num_nodes - 2)) / 2; + let mut savings = Vec::with_capacity(capacity); + for i in 1..num_nodes { + for j in (i + 1)..num_nodes { + savings.push((0.0, i as u8, j as u8)); + } + } + savings +} + +#[inline] +fn recompute_and_sort_savings(savings_list: &mut [(f32, u8, u8)], params: &[f32], challenge: &Challenge) { + let distance_matrix = &challenge.distance_matrix; + + let mut zero_len = 0; + for (score, i, j) in savings_list.iter_mut() { + let i = *i as usize; + let j = *j as usize; + *score = params[i] * distance_matrix[0][i] as f32 + + params[j] * distance_matrix[j][0] as f32 - + params[i] * params[j] * distance_matrix[i][j] as f32; + } + + savings_list.sort_unstable_by(|a, b| b.0.partial_cmp(&a.0).unwrap()); +} + +#[inline] +fn generate_neighbor(current: &[f32], rng: &mut R) -> Vec { + current.iter().map(|¶m| { + let delta = rng.gen_range(-0.1..=0.1); + (param + delta).clamp(0.0, 2.0) + }).collect() +} + +#[inline] +fn apply_local_search_until_no_improvement(solution: &mut Solution, distance_matrix: &Vec>) { + let mut improved = true; + while improved { + improved = false; + for route in &mut solution.routes { + if two_opt(route, distance_matrix) { + improved = true; + } + } + } +} +#[inline] +fn two_opt(route: &mut Vec, distance_matrix: &Vec>) -> bool { + let n = route.len(); + let mut improved = false; + + for i in 1..n - 2 { + for j in i + 1..n - 1 { + let current_distance = distance_matrix[route[i - 1]][route[i]] + + distance_matrix[route[j]][route[j + 1]]; + let new_distance = distance_matrix[route[i - 1]][route[j]] + + distance_matrix[route[i]][route[j + 1]]; + + if new_distance < current_distance { + route[i..=j].reverse(); + improved = true; + } + } + } + + improved +} + +#[inline] +fn calculate_solution_cost(solution: &Solution, distance_matrix: &Vec>) -> i32 { + solution.routes.iter().map(|route| { + route.windows(2).map(|w| distance_matrix[w[0]][w[1]]).sum::() + }).sum() +} + +#[inline] +fn create_solution(challenge: &Challenge, params: &[f32], savings_list: &[(f32, u8, u8)]) -> Solution { + let distance_matrix = &challenge.distance_matrix; + let max_capacity = challenge.max_capacity; + let num_nodes = challenge.difficulty.num_nodes; + let demands = &challenge.demands; + + let mut routes = vec![None; num_nodes]; + for i in 1..num_nodes { + routes[i] = Some(vec![i]); + } + let mut route_demands = demands.clone(); + + for &(_, i, j) in savings_list { + let (i, j) = (i as usize, j as usize); + if let (Some(left_route), Some(right_route)) = (routes[i].as_ref(), routes[j].as_ref()) { + let (left_start, left_end) = (*left_route.first().unwrap(), *left_route.last().unwrap()); + let (right_start, right_end) = (*right_route.first().unwrap(), *right_route.last().unwrap()); + + if left_start == right_start || route_demands[left_start] + route_demands[right_start] > max_capacity { + continue; + } + + let mut new_route = routes[i].take().unwrap(); + let mut right_route = routes[j].take().unwrap(); + + if left_start == i { new_route.reverse(); } + if right_end == j { right_route.reverse(); } + + new_route.extend(right_route); + + let combined_demand = route_demands[left_start] + route_demands[right_start]; + let new_start = new_route[0]; + let new_end = *new_route.last().unwrap(); + + route_demands[new_start] = combined_demand; + route_demands[new_end] = combined_demand; + + routes[new_start] = Some(new_route.clone()); + routes[new_end] = Some(new_route); + } + } + + Solution { + routes: routes + .into_iter() + .enumerate() + .filter_map(|(i, route)| route.filter(|r| r[0] == i)) + .map(|mut route| { + route.insert(0, 0); + route.push(0); + route + }) + .collect(), + } +} + +#[cfg(feature = "cuda")] +mod gpu_optimisation { + use super::*; + use cudarc::driver::*; + use std::{collections::HashMap, sync::Arc}; + use tig_challenges::CudaKernel; + + // set KERNEL to None if algorithm only has a CPU implementation + pub const KERNEL: Option = None; + + // Important! your GPU and CPU version of the algorithm should return the same result + pub fn cuda_solve_challenge( + challenge: &Challenge, + dev: &Arc, + mut funcs: HashMap<&'static str, CudaFunction>, + ) -> anyhow::Result> { + solve_challenge(challenge) + } +} +#[cfg(feature = "cuda")] +pub use gpu_optimisation::{cuda_solve_challenge, KERNEL}; \ No newline at end of file diff --git a/tig-algorithms/src/vehicle_routing/mod.rs b/tig-algorithms/src/vehicle_routing/mod.rs index 4f9879e4..b0d36866 100644 --- a/tig-algorithms/src/vehicle_routing/mod.rs +++ b/tig-algorithms/src/vehicle_routing/mod.rs @@ -97,7 +97,8 @@ pub use clarke_wright_super as c002_a036; // c002_a048 -// c002_a049 +pub mod advanced_routing; +pub use advanced_routing as c002_a049; // c002_a050 diff --git a/tig-algorithms/wasm/knapsack/classic_quadkp.wasm b/tig-algorithms/wasm/knapsack/classic_quadkp.wasm new file mode 100644 index 00000000..8fcd958c Binary files /dev/null and b/tig-algorithms/wasm/knapsack/classic_quadkp.wasm differ diff --git a/tig-algorithms/wasm/satisfiability/sat_global_opt.wasm b/tig-algorithms/wasm/satisfiability/sat_global_opt.wasm new file mode 100644 index 00000000..94b46b57 Binary files /dev/null and b/tig-algorithms/wasm/satisfiability/sat_global_opt.wasm differ diff --git a/tig-algorithms/wasm/vector_search/invector_hybrid.wasm b/tig-algorithms/wasm/vector_search/invector_hybrid.wasm new file mode 100644 index 00000000..2ba2b5f4 Binary files /dev/null and b/tig-algorithms/wasm/vector_search/invector_hybrid.wasm differ diff --git a/tig-algorithms/wasm/vehicle_routing/advanced_routing.wasm b/tig-algorithms/wasm/vehicle_routing/advanced_routing.wasm new file mode 100644 index 00000000..87cc27b4 Binary files /dev/null and b/tig-algorithms/wasm/vehicle_routing/advanced_routing.wasm differ