1use std::collections::{HashMap, HashSet};
2use std::fmt;
3use std::mem::{size_of, size_of_val};
4
5use omeco::{
6 CodeOptimizer, EinCode as OmecoEinCode, Initializer, NestedEinsum, ScoreFunction, TreeSA,
7};
8
9use crate::cache::{saturating_sum, vec_of_vec_retained_bytes, vec_retained_bytes};
10use crate::planning::plan::{compile_step_plans, DiagPlan, GemmPlan, ReducePlan, StepPlan};
11use crate::syntax::subscripts::Subscripts;
12use crate::util::{build_size_dict, intermediate_subs};
13use crate::{Error, Result};
14
15pub(crate) struct ContractionStep {
17 pub(crate) left: usize,
18 pub(crate) right: usize,
19}
20
21#[derive(Debug, Clone)]
27pub struct ContractionOptimizerOptions {
28 pub betas: Vec<f64>,
30 pub ntrials: usize,
32 pub niters: usize,
34 pub score: ScoreFunction,
36}
37
38impl Default for ContractionOptimizerOptions {
39 fn default() -> Self {
40 Self {
41 betas: Vec::new(),
42 ntrials: 1,
43 niters: 0,
44 score: ScoreFunction::default(),
45 }
46 }
47}
48
49impl ContractionOptimizerOptions {
50 fn to_treesa(&self) -> TreeSA {
51 TreeSA::new(
52 self.betas.clone(),
53 self.ntrials,
54 self.niters,
55 Initializer::Greedy,
56 self.score.clone(),
57 )
58 }
59
60 pub(crate) fn validate(&self) -> Result<()> {
61 if self.ntrials == 0 {
62 return Err(Error::planning(
63 "contraction optimizer ntrials must be at least 1",
64 ));
65 }
66 if self.betas.iter().any(|value| value.is_nan()) {
67 return Err(Error::planning(
68 "contraction optimizer betas must not contain NaN",
69 ));
70 }
71 if self.score.tc_weight.is_nan()
72 || self.score.sc_weight.is_nan()
73 || self.score.rw_weight.is_nan()
74 || self.score.sc_target.is_nan()
75 {
76 return Err(Error::planning(
77 "contraction optimizer score fields must not contain NaN",
78 ));
79 }
80 Ok(())
81 }
82}
83
84pub struct ContractionTree {
96 pub(crate) subscripts: Subscripts,
98 pub(crate) steps: Vec<ContractionStep>,
100 pub(crate) size_dict: HashMap<u32, usize>,
102 pub(crate) operand_subs: Vec<Vec<u32>>,
104 pub(crate) step_plans: Vec<StepPlan>,
106}
107
108impl fmt::Debug for ContractionTree {
109 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
110 f.debug_struct("ContractionTree")
111 .field("input_count", &self.subscripts.inputs.len())
112 .field("output_rank", &self.subscripts.output.len())
113 .field("steps_len", &self.steps.len())
114 .field("size_dict_len", &self.size_dict.len())
115 .field("operand_subs_len", &self.operand_subs.len())
116 .field("step_plans_len", &self.step_plans.len())
117 .finish_non_exhaustive()
118 }
119}
120
121impl ContractionTree {
122 pub fn optimize(subscripts: &Subscripts, shapes: &[&[usize]]) -> Result<Self> {
138 Self::optimize_with_options(subscripts, shapes, &ContractionOptimizerOptions::default())
139 }
140
141 pub fn optimize_with_options(
154 subscripts: &Subscripts,
155 shapes: &[&[usize]],
156 options: &ContractionOptimizerOptions,
157 ) -> Result<Self> {
158 options.validate()?;
159 let input_count = subscripts.inputs.len();
160 if input_count <= 1 {
161 return Self::from_pairs(subscripts, shapes, &[]);
162 }
163
164 let size_dict = build_size_dict(subscripts, shapes, None)?;
165 let pairs =
166 if let Some(omeco_pairs) = optimize_omeco_pairs(subscripts, &size_dict, options)? {
167 omeco_pairs
168 } else {
169 optimize_self_greedy_pairs(subscripts, &size_dict)?
170 };
171 Self::from_pairs(subscripts, shapes, &pairs)
172 }
173
174 pub fn from_pairs(
208 subscripts: &Subscripts,
209 shapes: &[&[usize]],
210 pairs: &[(usize, usize)],
211 ) -> Result<Self> {
212 let input_count = subscripts.inputs.len();
213 let required_steps = input_count.saturating_sub(1);
214 if pairs.len() != required_steps {
215 return Err(Error::planning(format!(
216 "explicit contraction path for {input_count} operands must have {required_steps} steps, got {}",
217 pairs.len()
218 )));
219 }
220 let size_dict = build_size_dict(subscripts, shapes, None)?;
221
222 let mut operand_subs: Vec<Vec<u32>> = subscripts.inputs.clone();
223 let mut live = vec![false; input_count + pairs.len()];
224 for slot in live.iter_mut().take(input_count) {
225 *slot = true;
226 }
227 let mut steps = Vec::new();
228
229 for (step_idx, &(left, right)) in pairs.iter().enumerate() {
230 let next_idx = input_count + step_idx;
231 if left == right {
232 return Err(Error::planning(format!(
233 "pair ({left}, {right}) must reference two distinct live operands"
234 )));
235 }
236 if left >= next_idx || right >= next_idx {
237 return Err(Error::planning(format!(
238 "pair ({left}, {right}) references non-existent operand"
239 )));
240 }
241 if !live[left] || !live[right] {
242 return Err(Error::planning(format!(
243 "pair ({left}, {right}) references an operand or intermediate that is no longer live"
244 )));
245 }
246
247 let mut needed: HashSet<u32> = subscripts.output.iter().copied().collect();
249 for (idx, subs) in operand_subs.iter().enumerate() {
250 if idx != left && idx != right && live[idx] {
251 needed.extend(subs.iter().copied());
252 }
253 }
254
255 let new_subs = intermediate_subs(&operand_subs[left], &operand_subs[right], &needed);
256 operand_subs.push(new_subs);
257 live[left] = false;
258 live[right] = false;
259 live[next_idx] = true;
260 steps.push(ContractionStep { left, right });
261 }
262
263 let live_count = live.iter().filter(|&&is_live| is_live).count();
264 if live_count != 1 {
265 return Err(Error::planning(format!(
266 "explicit contraction path must leave exactly one live result, got {live_count}"
267 )));
268 }
269
270 let mut tree = Self {
271 subscripts: subscripts.clone(),
272 steps,
273 size_dict,
274 operand_subs,
275 step_plans: Vec::new(),
276 };
277 tree.step_plans = compile_step_plans(&tree)?;
278 Ok(tree)
279 }
280
281 #[must_use]
298 pub fn step_count(&self) -> usize {
299 self.steps.len()
300 }
301
302 #[must_use]
322 pub fn step_pair(&self, step_idx: usize) -> Option<(usize, usize)> {
323 self.steps.get(step_idx).map(|step| (step.left, step.right))
324 }
325
326 #[must_use]
349 pub fn step_subscripts(&self, step_idx: usize) -> Option<(&[u32], &[u32], &[u32])> {
350 let input_count = self.subscripts.inputs.len();
351 let step = self.steps.get(step_idx)?;
352 let result_idx = input_count + step_idx;
353 let output_subs = if step_idx + 1 == self.steps.len() {
354 &self.subscripts.output
355 } else {
356 &self.operand_subs[result_idx]
357 };
358 Some((
359 &self.operand_subs[step.left],
360 &self.operand_subs[step.right],
361 output_subs,
362 ))
363 }
364
365 #[must_use]
378 pub fn step_plan(&self, step_idx: usize) -> Option<crate::lowering::PairwiseStepPlan<'_>> {
379 self.step_plans
380 .get(step_idx)
381 .map(crate::lowering::PairwiseStepPlan::new)
382 }
383
384 #[doc(hidden)]
385 #[must_use]
386 pub(crate) fn retained_bytes_for_cache_stats(&self) -> usize {
387 saturating_sum([
388 size_of::<Self>(),
389 subscripts_retained_bytes(&self.subscripts),
390 self.steps
391 .capacity()
392 .saturating_mul(size_of::<ContractionStep>()),
393 self.size_dict
394 .capacity()
395 .saturating_mul(size_of::<u32>().saturating_add(size_of::<usize>())),
396 vec_of_vec_retained_bytes(&self.operand_subs),
397 self.step_plans
398 .capacity()
399 .saturating_mul(size_of::<StepPlan>()),
400 saturating_sum(self.step_plans.iter().map(step_plan_retained_bytes)),
401 ])
402 }
403}
404
405fn subscripts_retained_bytes(subscripts: &Subscripts) -> usize {
406 saturating_sum([
407 vec_of_vec_retained_bytes(&subscripts.inputs),
408 vec_retained_bytes(&subscripts.output),
409 ])
410}
411
412fn reduce_plan_retained_bytes(plan: &ReducePlan) -> usize {
413 saturating_sum([
414 vec_retained_bytes(&plan.original_subs),
415 vec_retained_bytes(&plan.kept_subs),
416 vec_retained_bytes(&plan.out_shape),
417 ])
418}
419
420fn diag_plan_retained_bytes(plan: &DiagPlan) -> usize {
421 saturating_sum([
422 vec_retained_bytes(&plan.stages),
423 saturating_sum(plan.stages.iter().map(|stage| {
424 saturating_sum([
425 vec_retained_bytes(&stage.axis_pairs),
426 vec_retained_bytes(&stage.result_subs),
427 ])
428 })),
429 vec_retained_bytes(&plan.result_subs),
430 ])
431}
432
433fn gemm_plan_retained_bytes(plan: &GemmPlan) -> usize {
434 saturating_sum([
435 plan.reduce_a.as_ref().map_or(0, reduce_plan_retained_bytes),
436 plan.reduce_b.as_ref().map_or(0, reduce_plan_retained_bytes),
437 vec_retained_bytes(&plan.subs_a),
438 vec_retained_bytes(&plan.subs_b),
439 vec_retained_bytes(&plan.lo_modes),
440 vec_retained_bytes(&plan.ro_modes),
441 vec_retained_bytes(&plan.sum_modes),
442 vec_retained_bytes(&plan.lo_sizes),
443 vec_retained_bytes(&plan.ro_sizes),
444 vec_retained_bytes(&plan.sum_sizes),
445 vec_retained_bytes(&plan.batch_sizes),
446 vec_retained_bytes(&plan.target_a),
447 vec_retained_bytes(&plan.target_b),
448 vec_retained_bytes(&plan.c_gemm_shape),
449 vec_retained_bytes(&plan.expanded_shape),
450 vec_retained_bytes(&plan.canonical_modes),
451 vec_retained_bytes(&plan.a_gemm_shape),
452 vec_retained_bytes(&plan.b_gemm_shape),
453 ])
454}
455
456fn step_plan_retained_bytes(plan: &StepPlan) -> usize {
457 saturating_sum([
458 plan.diag_a.as_ref().map_or(0, diag_plan_retained_bytes),
459 plan.diag_b.as_ref().map_or(0, diag_plan_retained_bytes),
460 plan.strict_binary.as_ref().map_or(0, size_of_val),
461 gemm_plan_retained_bytes(&plan.gemm),
462 ])
463}
464
465fn optimize_omeco_pairs(
466 subscripts: &Subscripts,
467 size_dict: &HashMap<u32, usize>,
468 options: &ContractionOptimizerOptions,
469) -> Result<Option<Vec<(usize, usize)>>> {
470 let code = OmecoEinCode::new(subscripts.inputs.clone(), subscripts.output.clone());
471 let optimizer = options.to_treesa();
472 let Some(nested) = optimizer.optimize(&code, size_dict) else {
473 return Ok(None);
474 };
475
476 let mut next_operand = subscripts.inputs.len();
477 let mut pairs = Vec::with_capacity(subscripts.inputs.len().saturating_sub(1));
478 nested_to_pairs(&nested, &mut next_operand, &mut pairs)?;
479 Ok(Some(pairs))
480}
481
482fn nested_to_pairs(
483 nested: &NestedEinsum<u32>,
484 next_operand: &mut usize,
485 pairs: &mut Vec<(usize, usize)>,
486) -> Result<usize> {
487 match nested {
488 NestedEinsum::Leaf { tensor_index } => Ok(*tensor_index),
489 NestedEinsum::Node { args, .. } => {
490 if args.len() != 2 {
491 return Err(Error::planning(format!(
492 "omeco returned non-binary contraction node with {} children",
493 args.len()
494 )));
495 }
496 let left = nested_to_pairs(&args[0], next_operand, pairs)?;
497 let right = nested_to_pairs(&args[1], next_operand, pairs)?;
498 pairs.push((left, right));
499 let result_idx = *next_operand;
500 *next_operand += 1;
501 Ok(result_idx)
502 }
503 }
504}
505
506fn build_operand_label_sets(operand_subs: &[Vec<u32>]) -> Vec<HashSet<u32>> {
507 operand_subs
508 .iter()
509 .map(|subs| subs.iter().copied().collect())
510 .collect()
511}
512
513fn build_needed_label_counts(
514 output_subs: &[u32],
515 available: &[usize],
516 operand_label_sets: &[HashSet<u32>],
517) -> HashMap<u32, usize> {
518 let mut counts = HashMap::new();
519 for &label in output_subs {
520 counts.entry(label).or_insert(1);
521 }
522 for &idx in available {
523 add_labels_to_counts(&mut counts, &operand_label_sets[idx]);
524 }
525 counts
526}
527
528fn add_labels_to_counts(counts: &mut HashMap<u32, usize>, labels: &HashSet<u32>) {
529 for &label in labels {
530 *counts.entry(label).or_insert(0) += 1;
531 }
532}
533
534fn remove_labels_from_counts(counts: &mut HashMap<u32, usize>, labels: &HashSet<u32>) {
535 for &label in labels {
536 match counts.get(&label).copied() {
537 Some(1) => {
538 counts.remove(&label);
539 }
540 Some(count) => {
541 counts.insert(label, count - 1);
542 }
543 None => {}
544 }
545 }
546}
547
548fn candidate_label_is_needed(
549 label: u32,
550 left: usize,
551 right: usize,
552 operand_label_sets: &[HashSet<u32>],
553 needed_label_counts: &HashMap<u32, usize>,
554) -> bool {
555 let mut selected_count = 0;
556 if operand_label_sets[left].contains(&label) {
557 selected_count += 1;
558 }
559 if operand_label_sets[right].contains(&label) {
560 selected_count += 1;
561 }
562 needed_label_counts.get(&label).copied().unwrap_or(0) > selected_count
563}
564
565fn collect_candidate_intermediate_subs(
566 subs_left: &[u32],
567 subs_right: &[u32],
568 left: usize,
569 right: usize,
570 operand_label_sets: &[HashSet<u32>],
571 needed_label_counts: &HashMap<u32, usize>,
572 output: &mut Vec<u32>,
573) {
574 output.clear();
575 for &label in subs_left.iter().chain(subs_right.iter()) {
576 if candidate_label_is_needed(label, left, right, operand_label_sets, needed_label_counts)
577 && !output.contains(&label)
578 {
579 output.push(label);
580 }
581 }
582}
583
584#[derive(Clone, Copy)]
585struct CandidateCostContext<'a> {
586 operand_label_sets: &'a [HashSet<u32>],
587 needed_label_counts: &'a HashMap<u32, usize>,
588 size_dict: &'a HashMap<u32, usize>,
589}
590
591fn candidate_contraction_cost(
592 subs_left: &[u32],
593 subs_right: &[u32],
594 left: usize,
595 right: usize,
596 context: CandidateCostContext<'_>,
597 candidate_subs: &mut Vec<u32>,
598) -> Result<usize> {
599 collect_candidate_intermediate_subs(
600 subs_left,
601 subs_right,
602 left,
603 right,
604 context.operand_label_sets,
605 context.needed_label_counts,
606 candidate_subs,
607 );
608 let mut cost = 1usize;
609 for &label in candidate_subs.iter() {
610 let size = context.size_dict.get(&label).copied().ok_or_else(|| {
611 Error::planning(format!(
612 "unknown size for label {label} in contraction cost"
613 ))
614 })?;
615 cost = cost.saturating_mul(size);
616 }
617 Ok(cost.max(1))
618}
619
620fn optimize_self_greedy_pairs(
621 subscripts: &Subscripts,
622 size_dict: &HashMap<u32, usize>,
623) -> Result<Vec<(usize, usize)>> {
624 let input_count = subscripts.inputs.len();
625 let mut available: Vec<usize> = (0..input_count).collect();
626 let mut operand_subs: Vec<Vec<u32>> = subscripts.inputs.clone();
627 let mut operand_label_sets = build_operand_label_sets(&operand_subs);
628 let mut needed_label_counts =
629 build_needed_label_counts(&subscripts.output, &available, &operand_label_sets);
630 let mut candidate_subs = Vec::new();
631 let mut pairs: Vec<(usize, usize)> = Vec::new();
632
633 while available.len() > 1 {
634 let mut best_i = 0;
635 let mut best_j = 1;
636 let mut best_cost = usize::MAX;
637
638 for i in 0..available.len() {
639 for j in (i + 1)..available.len() {
640 let li = available[i];
641 let lj = available[j];
642 let cost = candidate_contraction_cost(
643 &operand_subs[li],
644 &operand_subs[lj],
645 li,
646 lj,
647 CandidateCostContext {
648 operand_label_sets: &operand_label_sets,
649 needed_label_counts: &needed_label_counts,
650 size_dict,
651 },
652 &mut candidate_subs,
653 )?;
654 if cost < best_cost {
655 best_cost = cost;
656 best_i = i;
657 best_j = j;
658 }
659 }
660 }
661
662 let left = available[best_i];
663 let right = available[best_j];
664 pairs.push((left, right));
665
666 let mut new_subs = Vec::new();
667 collect_candidate_intermediate_subs(
668 &operand_subs[left],
669 &operand_subs[right],
670 left,
671 right,
672 &operand_label_sets,
673 &needed_label_counts,
674 &mut new_subs,
675 );
676 let new_idx = operand_subs.len();
677 let new_label_set: HashSet<u32> = new_subs.iter().copied().collect();
678 remove_labels_from_counts(&mut needed_label_counts, &operand_label_sets[left]);
679 remove_labels_from_counts(&mut needed_label_counts, &operand_label_sets[right]);
680 add_labels_to_counts(&mut needed_label_counts, &new_label_set);
681 operand_subs.push(new_subs);
682 operand_label_sets.push(new_label_set);
683 available.remove(best_j);
684 available.remove(best_i);
685 available.push(new_idx);
686 }
687
688 Ok(pairs)
689}
690
691#[cfg(test)]
692mod tests;