nekolib/algo/
larsch.rs

1//! LARSCH algorithm。
2
3use std::cell::RefCell;
4use std::ops::Add;
5use std::rc::{Rc, Weak};
6
7enum Reduce<T, F, W> {
8    Reduce0E,
9    Reduce1E(Reduce1<T, F, W>),
10    Reduce2E(Reduce2<T, F, W>),
11}
12use Reduce::{Reduce0E, Reduce1E, Reduce2E};
13
14struct Fetcher<T, F, W> {
15    dp: RefCell<Vec<Option<T>>>,
16    map: F,
17    trans: W,
18}
19
20struct Reduce1<T, F, W> {
21    n: usize,
22    i: usize,
23    j0: usize,
24    j1: usize,
25    reduce: Box<Reduce<T, F, W>>,
26    row_shift: u32,
27    col_map: RefCell<Vec<usize>>,
28    c: RefCell<Vec<usize>>,
29    fetcher: Weak<RefCell<Fetcher<T, F, W>>>,
30}
31
32struct Reduce2<T, F, W> {
33    n: usize,
34    i: usize,
35    j: usize,
36    srow: Vec<usize>,
37    scol: Vec<usize>,
38    reduce: Box<Reduce<T, F, W>>,
39    row_shift: u32,
40    col_map: RefCell<Vec<usize>>,
41    c: RefCell<Vec<usize>>,
42    fetcher: Weak<RefCell<Fetcher<T, F, W>>>,
43}
44
45/// LARSCH algorithm。
46///
47/// $E\[i] = F\[j] + W(j, i)$ ($0\\le j<i<n$) と書ける DP を計算する。
48/// ここで、$F\[j]$ は $E\[j]$ から高速に計算できる値とし、$W$ は concave QI
49/// を満たすとする。
50///
51/// `todo!()` ちゃんと書く
52///
53/// $\\gdef\\DP{\operatorname{dp}}$
54/// 行列 $M\[i, j]$ は、$\\DP\[j]$ から $\\DP\[i]$ に遷移するときのコストを表し、
55/// 値がオンラインに定まる。この行列は concave totally monotone になっている。
56///
57/// この行列の row minima を求めるのに相当する。
58///
59/// `(argmin, min): (Vec<usize>, Vec<T>)` を返すので、`i = n - 1` から始めて
60/// `argmin[i] -> i` の遷移を辿ることで復元も可能。
61///
62/// # Notes
63///
64/// SMAWK のオンライン版。SMAWK もそのうち書く。
65///
66/// # Examples
67/// ```
68/// use nekolib::algo::Larsch;
69///
70/// let a = vec![1_i64, 2, 3, 4, 5];
71/// let c = 6;
72/// let n = a.len();
73///
74/// let map = |&x: &i64| x;
75/// let trans = |il: usize, ir: usize| (a[il] - a[ir]).pow(2) + c;
76/// let (argmin, min) = Larsch::new(n - 1, n - 1, 0, map, trans).solve();
77///
78/// assert_eq!(argmin, &[0, 0, 0, 0, 2]);
79/// assert_eq!(min, &[0, 7, 10, 15, 20]);
80/// ```
81///
82/// 愚直に書いたときの DP がどんな感じかを書いておくと役立つ場合がありそう。`todo!()`
83///
84/// # Complexity
85/// $E\[\\bullet]$ の計算を $O(1)$ time として $O(n)$ time。
86///
87/// # References
88/// - Larmore, Lawrence L., and Baruch Schieber. "On-line dynamic programming with applications to the prediction of RNA secondary structure." *Journal of Algorithms* 12, no. 3 (1991): 490--515.
89///
90/// ## See also
91/// - <https://noshi91.github.io/Library/algorithm/larsch.cpp>
92///
93/// # Applications
94/// [Examples](#examples) にある EDPC-Z の例の他、[AOJ 3086](https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=3086)
95/// などにも利用できる。
96///
97/// ```
98/// use std::cmp::Reverse;
99///
100/// use nekolib::{algo::Larsch, ds::N1Rmq};
101///
102/// let a = vec![1, 1, 5, 5, 1, 1];
103/// let l = 2;
104/// let n = a.len();
105///
106/// let a: Vec<_> = a.into_iter().map(|ai| Reverse(ai)).collect();
107/// let rmq: N1Rmq<_> = a.into();
108/// let oo = 10_i64.pow(18);
109///
110/// let map = |&x: &i64| x;
111/// let trans = |il: usize, ir: usize| {
112///     if ir < il + l { oo } else { -rmq.min(il, ir).0 }
113/// };
114/// let (argmin, min) = Larsch::new(n, n, 0, map, trans).solve();
115/// let max: Vec<_> = min.into_iter().map(|x| -x).collect();
116///
117/// assert_eq!(argmin, &[0, 0, 0, 0, 2, 3, 3]);
118/// assert_eq!(max, &[0, -oo, 1, 5, 6, 10, 10]);
119/// ```
120pub struct Larsch<T, F, W> {
121    n: usize,
122    reduce: Box<Reduce<T, F, W>>,
123    fetcher: Rc<RefCell<Fetcher<T, F, W>>>,
124}
125
126impl<T, F, W> Reduce<T, F, W>
127where
128    F: Fn(&T) -> T,
129    W: Fn(usize, usize) -> T,
130    T: Add<Output = T> + Ord,
131{
132    fn new(
133        n: usize,
134        m: usize,
135        row_shift: u32,
136        col_map: &[usize],
137        fetcher: Weak<RefCell<Fetcher<T, F, W>>>,
138    ) -> Self {
139        if n == 0 && m == 1 {
140            Reduce0E
141        } else if n >= m {
142            Reduce1E(Reduce1::new(n, m, row_shift, col_map, fetcher))
143        } else {
144            Reduce2E(Reduce2::new(n, m, row_shift, col_map, fetcher))
145        }
146    }
147
148    fn next(&mut self) -> Option<usize> {
149        match self {
150            Reduce0E => None,
151            Reduce1E(r) => r.next(),
152            Reduce2E(r) => r.next(),
153        }
154    }
155
156    fn set_col_map(&mut self, i_par: usize, i_child: usize) {
157        match self {
158            Reduce0E => {}
159            Reduce1E(r) => r.set_col_map(i_par, i_child),
160            Reduce2E(r) => r.set_col_map(i_par, i_child),
161        }
162    }
163
164    fn set_c(&mut self, i: usize, j: usize) {
165        match self {
166            Reduce0E => {}
167            Reduce1E(r) => r.set_c(i, j),
168            Reduce2E(r) => r.set_c(i, j),
169        }
170    }
171}
172
173impl<T, F, W> Fetcher<T, F, W>
174where
175    F: Fn(&T) -> T,
176    W: Fn(usize, usize) -> T,
177    T: Add<Output = T> + Ord,
178{
179    fn new(dp: Vec<Option<T>>, map: F, trans: W) -> Self {
180        let dp = RefCell::new(dp);
181        Self { dp, map, trans }
182    }
183
184    fn fetch(&self, i: usize, j: usize) -> T {
185        let f = (self.map)(self.dp.borrow()[j].as_ref().unwrap());
186        let w = (self.trans)(j, i);
187        let res = f + w;
188        res
189    }
190
191    fn set(&self, i: usize, val: T) { self.dp.borrow_mut()[i] = Some(val); }
192
193    fn finish(self) -> Vec<T> {
194        self.dp.into_inner().into_iter().map(Option::unwrap).collect()
195    }
196}
197
198impl<T, F, W> Larsch<T, F, W>
199where
200    F: Fn(&T) -> T,
201    W: Fn(usize, usize) -> T,
202    T: Add<Output = T> + Ord,
203{
204    pub fn new(n: usize, m: usize, init: T, map: F, trans: W) -> Self {
205        let dp: Vec<_> = (0..=n).map(|_| None).collect();
206        let fetcher = Fetcher::new(dp, map, trans);
207        fetcher.set(0, init);
208        let fetcher = Rc::new(RefCell::new(fetcher));
209        let row_shift = 0;
210        let col_map: Vec<_> = (0..m).collect();
211        let reduce =
212            Reduce::new(n, m, row_shift, &col_map, Rc::downgrade(&fetcher));
213        Self { n, reduce: Box::new(reduce), fetcher }
214    }
215
216    pub fn solve(mut self) -> (Vec<usize>, Vec<T>) {
217        let mut argmin = vec![0; self.n + 1];
218        for i in 1..=self.n {
219            self.reduce.set_c(i, i - 1);
220            argmin[i] = self.reduce.next().unwrap();
221            let x = self.fetcher.borrow().fetch(i, argmin[i]);
222            self.fetcher.borrow().set(i, x);
223        }
224
225        let dp: Vec<_> =
226            Rc::try_unwrap(self.fetcher).ok().unwrap().into_inner().finish();
227
228        (argmin, dp)
229    }
230}
231
232impl<T, F, W> Reduce1<T, F, W>
233where
234    F: Fn(&T) -> T,
235    W: Fn(usize, usize) -> T,
236    T: Add<Output = T> + Ord,
237{
238    fn new(
239        n: usize,
240        m: usize,
241        row_shift: u32,
242        col_map: &[usize],
243        fetcher: Weak<RefCell<Fetcher<T, F, W>>>,
244    ) -> Self {
245        let reduce = Reduce::new(
246            n / 2,
247            m,
248            row_shift + 1,
249            col_map,
250            Weak::clone(&fetcher),
251        );
252
253        Self {
254            n,
255            i: 1,
256            j0: 0,
257            j1: 0,
258            reduce: Box::new(reduce),
259            row_shift,
260            col_map: RefCell::new(col_map.to_vec()),
261            c: RefCell::new(vec![0; n + 1]),
262            fetcher,
263        }
264    }
265
266    fn c(&self, i: usize) -> usize { self.c.borrow()[i] }
267
268    fn set_col_map(&mut self, i_par: usize, i_child: usize) {
269        self.col_map.borrow_mut()[i_child] = i_par;
270        self.reduce.set_col_map(i_par, i_child);
271    }
272
273    fn set_c(&mut self, i: usize, j: usize) {
274        self.c.borrow_mut()[i] = j;
275        if i / 2 + 1 <= self.n / 2 {
276            self.reduce.set_c(i / 2 + 1, j);
277        }
278    }
279
280    fn fetch(&self, i: usize, j: usize) -> T {
281        let i = i << self.row_shift;
282        let j = self.col_map.borrow()[j];
283        self.fetcher.upgrade().unwrap().borrow().fetch(i, j)
284    }
285
286    fn next(&mut self) -> Option<usize> {
287        if self.i > self.n {
288            return None;
289        }
290
291        let j_range = if self.i % 2 == 1 {
292            if self.i % 2 == 1 && self.i / 2 + 1 <= self.n / 2 {
293                self.reduce.set_c(self.i / 2 + 1, self.c(self.i));
294            }
295            self.j1 = self.reduce.next().unwrap_or_else(|| self.c(self.i));
296            (self.j0..=self.j1).chain(None)
297        } else {
298            (self.c(self.i - 1) + 1..=self.c(self.i)).chain(Some(self.j1))
299        };
300
301        let j = j_range.min_by_key(|&j| (self.fetch(self.i, j), j)).unwrap();
302        if self.i % 2 == 0 {
303            self.j0 = j;
304        }
305
306        self.i += 1;
307        Some(j)
308    }
309}
310
311impl<T, F, W> Reduce2<T, F, W>
312where
313    F: Fn(&T) -> T,
314    W: Fn(usize, usize) -> T,
315    T: Add<Output = T> + Ord,
316{
317    fn new(
318        n: usize,
319        _m: usize,
320        row_shift: u32,
321        col_map: &[usize],
322        fetcher: Weak<RefCell<Fetcher<T, F, W>>>,
323    ) -> Self {
324        let reduce =
325            Reduce::new(n, n, row_shift, &vec![0; n], Weak::clone(&fetcher));
326
327        Self {
328            n,
329            i: 1,
330            j: 0,
331            srow: vec![0],
332            scol: vec![0],
333            reduce: Box::new(reduce),
334            row_shift,
335            col_map: RefCell::new(col_map.to_vec()),
336            c: RefCell::new(vec![0; n + 1]),
337            fetcher,
338        }
339    }
340
341    fn c(&self, i: usize) -> usize { self.c.borrow()[i] }
342
343    fn is_finite(&self, i: usize, j: usize) -> bool {
344        j <= self.c(i.min(self.i))
345    }
346
347    fn set_col_map(&mut self, i_par: usize, i_child: usize) {
348        self.col_map.borrow_mut()[i_child] = i_par;
349    }
350
351    fn set_c(&mut self, i: usize, j: usize) { self.c.borrow_mut()[i] = j; }
352
353    fn fetch(&self, i: usize, j: usize) -> Option<T> {
354        if !self.is_finite(i, j) {
355            return None;
356        }
357
358        let i = i << self.row_shift;
359        let j = self.col_map.borrow()[j];
360
361        Some(self.fetcher.upgrade().unwrap().borrow().fetch(i, j))
362    }
363
364    fn is_lt(
365        &self,
366        (il, jl): (usize, usize),
367        (ir, jr): (usize, usize),
368    ) -> bool {
369        if (il, jl) == (ir, jr) {
370            return false;
371        }
372        match (self.fetch(il, jl), self.fetch(ir, jr)) {
373            (Some(fl), Some(fr)) => (fl, jl) < (fr, jr),
374            (_, None) => true,  // finite < oo
375            (None, _) => false, // oo < finite
376        }
377    }
378
379    fn next(&mut self) -> Option<usize> {
380        if self.i > self.n {
381            return None;
382        }
383
384        // c[0] := -1
385        let jl = if self.i == 1 { 0 } else { self.c(self.i - 1) + 1 };
386        let jr = self.c(self.i);
387        for j in jl..=jr {
388            loop {
389                let r = *self.srow.last().unwrap();
390                let c = *self.scol.last().unwrap();
391                if !self.is_lt((r, j), (r, c)) {
392                    break;
393                }
394                self.srow.pop();
395                self.scol.pop();
396            }
397            if *self.srow.last().unwrap() < self.n {
398                let r = *self.srow.last().unwrap();
399                let i_ = (r + 1..).find(|&i_| self.is_finite(i_, j)).unwrap();
400                self.srow.push(i_);
401                self.scol.push(j);
402            }
403        }
404        if self.j + 1 < self.srow.len() && self.srow[self.j + 1] == self.i {
405            self.j += 1;
406            let tmp = self.col_map.borrow()[self.scol[self.j]];
407            self.reduce.set_col_map(tmp, self.j - 1);
408        }
409
410        self.reduce.set_c(self.i, self.j - 1);
411        let j = self.reduce.next().unwrap();
412        self.i += 1;
413        Some(self.scol[j + 1])
414    }
415}
nekolib/algo/larsch.rs

nekolib/algo/
larsch.rs
