puanrs/
linalg.rs

1//! # Linalg functions
2//! 
3//! Handy tools from linear algebra
4//! 
5use std::collections::HashMap;
6
7use itertools::Itertools;
8/// Data structure for matrix
9#[derive(Debug)]
10#[derive(Default)]
11pub struct Matrix {
12    /// `Vec` holding the values of the matrix. Note that `val.len()` must be equal to the product of `ncols` and `nrows`.
13    pub val: Vec<f64>,
14    /// Number of columns of the matrix
15    pub ncols: usize,
16    /// Number of rows of the matrix
17    pub nrows: usize
18}
19
20impl Clone for Matrix {
21    fn clone(&self) -> Self {
22        return Matrix {
23            val : self.val.to_vec(),
24            ncols: self.ncols,
25            nrows: self.nrows
26        }
27    }
28}
29
30impl PartialEq for Matrix {
31    fn eq(&self, other: &Self) -> bool {
32        (self.val == other.val) & (self.ncols == other.ncols) & (self.nrows == other.nrows)
33    }
34}
35
36/// Creates an identity matrix based on the input size `n`
37pub fn identity_matrix(n: usize) -> Matrix {
38    let mut t: Vec<f64> = Vec::with_capacity(n*n);
39    for r in 0..n {
40        for i in 0..n {
41            if r == i {
42                t.push(1.);
43            } else {
44                t.push(0.);
45            }
46        }
47    }
48    Matrix { val: t, ncols: n, nrows: n}
49}
50
51/// Calculates the dot product between two matrices
52pub fn dot(mat1: &Matrix, mat2: &Matrix) -> Matrix{
53    if !(mat1.ncols == mat2.nrows){
54        panic!("Dimensions does not match, cannot calculate the dot product between matrices of shapes ({},{}) and ({},{})", mat1.nrows, mat1.ncols, mat2.nrows, mat2.ncols);
55    }
56    let mut result: Vec<f64> = Vec::with_capacity(mat1.nrows*mat2.ncols);
57    for i in 0..mat1.nrows {
58        for j in 0..mat2.ncols {
59            result.push(mat1.val[i*mat1.ncols..(i+1)*mat1.ncols].to_vec().iter().zip(mat2.val.iter().skip(j).step_by(mat2.ncols)).map(|(x, y)| x * y).sum());
60        }
61    }
62    return Matrix {
63        val: result,
64        ncols: mat2.ncols,
65        nrows: mat1.nrows
66    }
67}
68
69/// Divides one matrix by another.
70/// 
71/// If the matrices have the same shape each element $a_{ij}$ in the numerator matrix is divided by the corresponding element $b_{ij}$ in the denominator matrix.
72/// It is also possible to divide a matrix with another matrix with only one row, if they share the same number of columns. In this case, each row in the nominator matrix is
73/// divided by the denominator matrix. 
74/// 
75/// # Example:
76/// 
77/// $$ \[ a \ b \] / \[ c \ d \] = \[a/c \ b/d \]$$
78/// 
79/// ```
80/// use puanrs::linalg::Matrix;
81/// use puanrs::linalg::divide;
82/// // Divide 2x3 matrix with 2x3 matrix
83/// let res = divide(
84///     &Matrix{
85///         val: vec![1., 1., 1., 1., 1., 1.],
86///         ncols: 3,
87///         nrows: 2},
88///     &Matrix{
89///         val: vec![1., 2., 3., -1., 1., 0.],
90///         ncols: 3,
91///         nrows: 2});
92/// assert_eq!(format!("{res:#?}"),
93/// "Matrix {
94///     val: [
95///         1.0,
96///         0.5,
97///         0.3333333333333333,
98///         -1.0,
99///         1.0,
100///         1.7976931348623157e308,
101///     ],
102///     ncols: 3,
103///     nrows: 2,
104/// }")
105/// ```
106/// ```
107/// use puanrs::linalg::Matrix;
108/// use puanrs::linalg::divide;
109/// // Divide 2x3 matrix with 1x3 matrix
110/// let res = divide(
111///     &Matrix{
112///         val: vec![1., 1., 1., 6., 6., 6.],
113///         ncols: 3,
114///         nrows: 2},
115///     &Matrix{
116///         val: vec![1., 2., 3.],
117///         ncols: 3,
118///         nrows: 1});
119/// assert_eq!(format!("{res:#?}"),
120/// "Matrix {
121///     val: [
122///         1.0,
123///         0.5,
124///         0.3333333333333333,
125///         6.0,
126///         3.0,
127///         2.0,
128///     ],
129///     ncols: 3,
130///     nrows: 2,
131/// }")
132/// ```
133
134pub fn divide(mat1: &Matrix, mat2: &Matrix) -> Matrix {
135    if !(mat1.ncols == mat2.ncols && (mat1.nrows == mat2.nrows || mat2.nrows == 1)) {
136        panic!("Dimensions does not match, cannot divide a matrix of shape ({},{}) with a matrix of shape ({},{})", mat1.nrows, mat1.ncols, mat2.nrows, mat2.ncols);
137    }
138    let mut result: Vec<f64> = Vec::with_capacity(mat1.ncols*mat1.nrows);
139    for i in 0..mat1.nrows{
140        for j in 0..mat1.ncols {
141            if mat2.nrows > 1 && mat2.val[i*mat1.ncols+j] != 0.0 {
142                result.push(mat1.val[i*mat1.ncols+j]/mat2.val[i*mat1.ncols+j]);
143            } else if mat2.nrows == 1 && mat2.val[j] != 0.0 {
144                result.push(mat1.val[i*mat1.ncols+j]/mat2.val[j]);
145            } else {
146                result.push(f64::MAX);
147            }
148        }
149    }
150    return Matrix { val: result, ncols: mat1.ncols, nrows: mat1.nrows }
151}
152
153/// Divides one matrix by another if the divisor is greater than zero, otherwise the value is set to `f64::MAX`. (See [divide])
154/// 
155/// 
156/// # Example:
157/// 
158/// ```
159/// use puanrs::linalg::Matrix;
160/// use puanrs::linalg::ge_divide;
161/// // Divide 2x3 matrix with 2x3 matrix
162/// let res = ge_divide(
163///     &Matrix{
164///         val: vec![1., 1., 1., 1., 1., 1.],
165///         ncols: 3,
166///         nrows: 2},
167///     &Matrix{
168///         val: vec![1., 2., 3., -1., 1., 0.],
169///         ncols: 3,
170///         nrows: 2});
171/// assert_eq!(format!("{res:#?}"),
172/// "Matrix {
173///     val: [
174///         1.0,
175///         0.5,
176///         0.3333333333333333,
177///         1.7976931348623157e308,
178///         1.0,
179///         1.7976931348623157e308,
180///     ],
181///     ncols: 3,
182///     nrows: 2,
183/// }")
184/// ```
185pub fn ge_divide(mat1: &Matrix, mat2: &Matrix) -> Matrix {
186    if !(mat1.ncols == mat2.ncols && (mat1.nrows == mat2.nrows || mat2.nrows == 1)) {
187        panic!("Dimensions does not match, cannot divide a matrix of shape ({},{}) with a matrix of shape ({},{})", mat1.nrows, mat1.ncols, mat2.nrows, mat2.ncols);
188    }
189    let mut result: Vec<f64> = Vec::with_capacity(mat1.ncols*mat1.nrows);
190    for i in 0..mat1.nrows{
191        for j in 0..mat1.ncols {
192            if mat2.nrows > 1 && mat2.val[i*mat1.ncols+j] > 0.0 {
193                result.push(mat1.val[i*mat1.ncols+j]/mat2.val[i*mat1.ncols+j]);
194            } else if mat2.nrows == 1 && mat2.val[j] > 0.0 {
195                result.push(mat1.val[i*mat1.ncols+j]/mat2.val[j]);
196            } else {
197                result.push(f64::MAX);
198            }
199        }
200    }
201    return Matrix { val: result, ncols: mat1.ncols, nrows: mat1.nrows }
202}
203
204/// Divides one matrix by another if the divisor is less than zero, otherwise the value is set to `f64::MAX`. (See [divide])
205/// 
206/// # Example:
207/// 
208/// ```
209/// use puanrs::linalg::Matrix;
210/// use puanrs::linalg::le_divide;
211/// // Divide 2x3 matrix with 2x3 matrix
212/// let res = le_divide(
213///     &Matrix{
214///         val: vec![1., 1., 1., -1., -1., -1.],
215///         ncols: 3,
216///         nrows: 2},
217///     &Matrix{
218///         val: vec![-1., -2., -3., 1., -1., 0.],
219///         ncols: 3,
220///         nrows: 2});
221/// assert_eq!(format!("{res:#?}"),
222/// "Matrix {
223///     val: [
224///         -1.0,
225///         -0.5,
226///         -0.3333333333333333,
227///         1.7976931348623157e308,
228///         1.0,
229///         1.7976931348623157e308,
230///     ],
231///     ncols: 3,
232///     nrows: 2,
233/// }")
234/// ```
235pub fn le_divide(mat1: &Matrix, mat2: &Matrix) -> Matrix {
236    if !(mat1.ncols == mat2.ncols && (mat1.nrows == mat2.nrows || mat2.nrows == 1)) {
237        panic!("Dimensions does not match, cannot divide a matrix of shape ({},{}) with a matrix of shape ({},{})", mat1.nrows, mat1.ncols, mat2.nrows, mat2.ncols);
238    }
239    let mut result: Vec<f64> = Vec::with_capacity(mat1.ncols*mat1.nrows);
240    for i in 0..mat1.nrows{
241        for j in 0..mat1.ncols {
242            if mat2.nrows > 1 && mat2.val[i*mat1.ncols+j] < 0.0 {
243                result.push(mat1.val[i*mat1.ncols+j]/mat2.val[i*mat1.ncols+j]);
244            } else if mat2.nrows == 1 && mat2.val[j] < 0.0 {
245                result.push(mat1.val[i*mat1.ncols+j]/mat2.val[j]);
246            } else {
247                result.push(f64::MAX);
248            }
249        }
250    }
251    return Matrix { val: result, ncols: mat1.ncols, nrows: mat1.nrows }
252}
253
254/// Adds one matrix by another.
255/// 
256/// If the matrices have the same shape each element $a_{ij}$ in matrix 1 is added by the corresponding element $b_{ij}$ in matrix 2.
257/// It is also possible to add a matrix with a singled rowed matrix, if they share the same number of columns. In this case, each row in matrix 1 is
258/// added by matrix 2. 
259/// 
260/// # Example:
261/// 
262/// $$ \[ a \ b \] + \[ c \ d \] = \[a+c \ b+d \]$$
263/// 
264/// ```
265/// use puanrs::linalg::Matrix;
266/// use puanrs::linalg::add;
267/// // Add a 2x3 matrix with a 2x3 matrix
268/// let res = add(
269///     &Matrix{
270///         val: vec![1., 1., 1., 1., 1., 1.],
271///         ncols: 3,
272///         nrows: 2},
273///     &Matrix{
274///         val: vec![1., 2., 3., -1., 1., 0.],
275///         ncols: 3,
276///         nrows: 2});
277/// assert_eq!(format!("{res:#?}"),
278/// "Matrix {
279///     val: [
280///         2.0,
281///         3.0,
282///         4.0,
283///         0.0,
284///         2.0,
285///         1.0,
286///     ],
287///     ncols: 3,
288///     nrows: 2,
289/// }")
290/// ```
291/// ```
292/// use puanrs::linalg::Matrix;
293/// use puanrs::linalg::add;
294/// // Add a 2x3 matrix with a 1x3 matrix
295/// let res = add(
296///     &Matrix{
297///         val: vec![1., 1., 1., 6., 6., 6.],
298///         ncols: 3,
299///         nrows: 2},
300///     &Matrix{
301///         val: vec![1., 2., 3.],
302///         ncols: 3,
303///         nrows: 1});
304/// assert_eq!(format!("{res:#?}"),
305/// "Matrix {
306///     val: [
307///         2.0,
308///         3.0,
309///         4.0,
310///         7.0,
311///         8.0,
312///         9.0,
313///     ],
314///     ncols: 3,
315///     nrows: 2,
316/// }")
317/// ```
318pub fn add(mat1: &Matrix, mat2: &Matrix) -> Matrix {
319    if !(mat1.ncols == mat2.ncols && (mat1.nrows == mat2.nrows || mat2.nrows == 1)) {
320        panic!("Dimensions does not match, cannot add a matrix of shape ({},{}) with a matrix of shape ({},{})", mat1.nrows, mat1.ncols, mat2.nrows, mat2.ncols);
321    }
322    let mut result: Vec<f64> = Vec::with_capacity(mat1.ncols*mat1.nrows);
323    for i in 0..mat1.nrows{
324        for j in 0..mat1.ncols {
325            if mat2.nrows > 1 {
326                result.push(mat1.val[i*mat1.ncols+j]+mat2.val[i*mat1.ncols+j]);
327            } else {
328                result.push(mat1.val[i*mat1.ncols+j]+mat2.val[j]);
329            } 
330        }
331    }
332    return Matrix { val: result, ncols: mat1.ncols, nrows: mat1.nrows }
333}
334
335/// Subtracts one matrix by another.
336/// 
337/// If the matrices have the same shape each element $a_{ij}$ in matrix 1 is subtracted by the corresponding element $b_{ij}$ in matrix 2.
338/// It is also possible to subtract a matrix with another matrix with only one row, if they share the same number of columns. In this case, each row in matrix 1 is
339/// subtracted by matrix 2. 
340/// 
341/// # Example:
342/// 
343/// $$ \[ a \ b \] - \[ c \ d \] = \[a-c \ b-d \]$$
344/// 
345/// ```
346/// use puanrs::linalg::Matrix;
347/// use puanrs::linalg::subtract;
348/// // Subtract a 2x3 matrix with a 2x3 matrix
349/// let res = subtract(
350///     &Matrix{
351///         val: vec![1., 1., 1., 1., 1., 1.],
352///         ncols: 3,
353///         nrows: 2},
354///     &Matrix{
355///         val: vec![1., 2., 3., -1., 1., 0.],
356///         ncols: 3,
357///         nrows: 2});
358/// assert_eq!(format!("{res:#?}"),
359/// "Matrix {
360///     val: [
361///         0.0,
362///         -1.0,
363///         -2.0,
364///         2.0,
365///         0.0,
366///         1.0,
367///     ],
368///     ncols: 3,
369///     nrows: 2,
370/// }")
371/// ```
372/// ```
373/// use puanrs::linalg::Matrix;
374/// use puanrs::linalg::subtract;
375/// // Subtract a 2x3 matrix with a 1x3 matrix
376/// let res = subtract(
377///     &Matrix{
378///         val: vec![1., 1., 1., 6., 6., 6.],
379///         ncols: 3,
380///         nrows: 2},
381///     &Matrix{
382///         val: vec![1., 2., 3.],
383///         ncols: 3,
384///         nrows: 1});
385/// assert_eq!(format!("{res:#?}"),
386/// "Matrix {
387///     val: [
388///         0.0,
389///         -1.0,
390///         -2.0,
391///         5.0,
392///         4.0,
393///         3.0,
394///     ],
395///     ncols: 3,
396///     nrows: 2,
397/// }")
398/// ```
399pub fn subtract(mat1: &Matrix, mat2: &Matrix) -> Matrix {
400    if !(mat1.ncols == mat2.ncols && (mat1.nrows == mat2.nrows || mat2.nrows == 1)) {
401        panic!("Dimensions does not match, cannot subtract a matrix of shape ({},{}) with a matrix of shape ({},{})", mat1.nrows, mat1.ncols, mat2.nrows, mat2.ncols);
402    }
403    let mut result: Vec<f64> = Vec::with_capacity(mat1.ncols*mat1.nrows);
404    for i in 0..mat1.nrows{
405        for j in 0..mat1.ncols {
406            if mat2.nrows > 1 {
407                result.push(mat1.val[i*mat1.ncols+j]-mat2.val[i*mat1.ncols+j]);
408            } else {
409                result.push(mat1.val[i*mat1.ncols+j]-mat2.val[j]);
410            } 
411        }
412    }
413    return Matrix { val: result, ncols: mat1.ncols, nrows: mat1.nrows }
414}
415
416/// Transpose the input Matrix and returns the result as a new Matrix
417/// 
418/// # Example:
419/// 
420/// ```
421/// use puanrs::linalg::Matrix;
422/// use puanrs::linalg::transpose;
423/// let res = transpose(&Matrix{val: vec![1., 2., 3., 4., 5., 6.], ncols: 3, nrows: 2});
424/// assert_eq!(format!("{res:#?}"),
425/// "Matrix {
426///     val: [
427///         1.0,
428///         4.0,
429///         2.0,
430///         5.0,
431///         3.0,
432///         6.0,
433///     ],
434///     ncols: 2,
435///     nrows: 3,
436/// }")
437/// ```
438
439pub fn transpose(mat: &Matrix) -> Matrix{
440    let mut result = Vec::with_capacity(mat.val.len());
441    for i in 0..mat.ncols{
442        for j in 0..mat.nrows{
443            result.push(mat.val[j*mat.ncols + i])
444        }
445    }
446    return Matrix{val: result, nrows: mat.ncols, ncols: mat.nrows}
447}
448
449/// Eliminates column coefficient `col` of row `dst_row` by Gauss elemination with row `src_row`.
450pub fn gauss_elimination(mat: &Matrix, dst_row: usize, src_row: usize, col: usize) -> Matrix{
451    if mat.val[src_row*mat.ncols+col] == 0.0 {
452        panic!("Cannot perform Gauss elimination when 'col' coefficient in the 'src_row' is zero");
453    }
454    let corr = mat.val[dst_row*mat.ncols + col]/mat.val[src_row*mat.ncols+col];
455    row_subtraction(mat, dst_row, src_row, Some(corr))
456}
457
458/// Adds `src_row` to `dst_row`. If multiplier is given the `src_row` is multiplied by this number before the addition.
459pub fn row_addition(mat: &Matrix, dst_row: usize, src_row: usize, multiplier: Option<f64>) -> Matrix{
460    let mul: f64;
461    if multiplier.is_none() {
462        mul = 1.0;
463    } else {
464        mul = multiplier.unwrap();
465    }
466    let mut tmp = Vec::with_capacity(mat.val.len());
467    for i in 0..mat.val.len(){
468        if i >= dst_row*mat.ncols && i < dst_row*mat.ncols+mat.ncols{
469            tmp.push(mat.val[i]+mat.val[src_row*mat.ncols+i%mat.ncols]*mul);
470        } else {
471            tmp.push(mat.val[i]);
472        }
473    }
474    return Matrix{val: tmp, ncols: mat.ncols, nrows: mat.nrows}
475}
476
477/// Subtracts `dst_row` with `src_row`. If multiplier is given the `src_row` is multiplied by this number before the subtraction.
478pub fn row_subtraction(mat: &Matrix, dst_row: usize, src_row: usize, multiplier: Option<f64>) -> Matrix{
479    let mul: f64;
480    if multiplier.is_none() {
481        mul = 1.0;
482    } else {
483        mul = multiplier.unwrap();
484    }
485    let mut tmp = Vec::with_capacity(mat.val.len());
486    for i in 0..mat.val.len(){
487        if i >= dst_row*mat.ncols && i < dst_row*mat.ncols+mat.ncols{
488            tmp.push(mat.val[i]-mat.val[src_row*mat.ncols+i%mat.ncols]*mul);
489        } else {
490            tmp.push(mat.val[i]);
491        }
492    }
493    return Matrix{val: tmp, ncols: mat.ncols, nrows: mat.nrows}
494}
495
496/// Returns a new Matrix with the coluns specified by `ind`. 
497pub fn get_columns(mat: &Matrix, ind: &Vec<usize>) -> Matrix {
498    let mut result: Vec<f64> = Vec::with_capacity(ind.len());
499    for i in 0..mat.nrows {
500        result.extend(ind.iter().map(|j| mat.val[i*mat.ncols+j]).collect::<Vec<f64>>());
501    }
502    return Matrix { val: result, ncols: ind.len(), nrows: mat.nrows }   
503}
504
505/// Returns a new Matrix with column `ind` replaced with `v`. 
506pub fn update_column(mat: &Matrix, ind: usize, v: &Vec<f64>) -> Matrix{
507    if mat.nrows != v.len() {
508        panic!("Dimension does not match");
509    }
510    let mut result = mat.val.to_vec();
511    for i in 0..v.len() {
512        result[i*mat.ncols+ind] = v[i];
513    }
514    Matrix { val: result, ncols: mat.ncols, nrows: mat.nrows}
515}
516
517impl Matrix {
518    /// See [dot]
519    pub fn dot(&self, mat2: &Matrix) -> Matrix{
520        dot(self, mat2)
521    }
522    /// See [ge_divide]
523    pub fn ge_divide(&self, mat2: &Matrix) -> Matrix {
524        ge_divide(self, mat2)
525    }
526    /// See [le_divide]
527    pub fn le_divide(&self, mat2: &Matrix) -> Matrix {
528        le_divide(self, mat2)
529    }
530    /// See [divide]
531    pub fn divide(&self, mat2: &Matrix) -> Matrix {
532        divide(self, mat2)
533    }
534    /// See [subtract]
535    pub fn subtract(&self, mat2: &Matrix) -> Matrix {
536        subtract(self, mat2)
537    }
538    /// See [transpose]
539    pub fn transpose(&self) -> Matrix{
540        transpose(self)
541    }
542    pub fn insert_column(&self, column: usize, elem: Vec<f64>) -> Matrix {
543    if column > self.ncols {
544        panic!("Cannot insert column at {} to matrix with dimensions {}, {}", column, self.nrows, self.ncols);
545    }
546    let mut result = Vec::new();
547    let mut j = 0;
548    for (i, e) in self.val.iter().enumerate() {
549        if (i + self.ncols - column) > 0 && (i + self.ncols - column) % self.ncols == 0 {
550            result.push(elem[j].clone());
551            j = j + 1;
552        }
553        result.push(e.clone());
554    }
555    if (j + 1) == elem.len() {
556        result.push(elem[j])
557    }
558    return Matrix { val: result, ncols: self.ncols+1, nrows: self.nrows }
559    }
560    pub fn gauss_elimination(&self, row: usize, by: usize, col: usize) -> Matrix{
561        gauss_elimination(self, row, by, col)
562    }
563
564    pub fn row_addition(&self, dst_row: usize, src_row: usize, multiplier: Option<f64>) -> Matrix{
565        row_addition(self, dst_row, src_row, multiplier)
566    }
567
568    pub fn get_columns(&self, ind: &Vec<usize>) -> Matrix {
569        get_columns(self, ind)
570    }
571
572    pub fn update_column(&self, ind: usize, v: &Vec<f64>) -> Matrix{
573        update_column(self, ind, v)
574    }
575}
576
577
578#[derive(Default, Clone, Debug, PartialEq)]
579pub struct CsrMatrix {
580    /// The value for an index `i` in `row` is representing a row index in a virtual matrix
581    pub row: Vec<i64>,
582    /// The value for an index `j` in `col` is representing a column index in a virtual matrix
583    pub col: Vec<i64>,
584    /// The value for an element is the value for the cell (`i`,`j`) in a virtual matrix
585    pub val: Vec<f64>
586}
587
588impl CsrMatrix {
589
590    /// Converts into an ordinary Matrix
591    pub fn to_matrix(&self) -> Matrix {
592        let unique_rows: Vec<&i64> = self.row.iter().unique().sorted().collect_vec();
593        let unique_cols: Vec<&i64> = self.col.iter().unique().sorted().collect_vec();
594        let n_unique_rows = unique_rows.len();
595        let m_unique_cols = unique_cols.len();
596
597        let rows_map: HashMap<&i64, usize> = unique_rows.into_iter().zip(0..n_unique_rows).collect();
598        let cols_map: HashMap<&i64, usize> = unique_cols.into_iter().zip(0..m_unique_cols).collect();
599
600        let mut val: Vec<f64> = vec![0.0; n_unique_rows*m_unique_cols];
601        for (i,(j,v)) in self.row.iter().zip(self.col.iter().zip(self.val.iter())) {
602            let i_mapped = rows_map.get(i).expect("could not find row i in mapping");
603            let j_mapped = cols_map.get(j).expect("could not find col j in mapping");
604            val[((*i_mapped as usize) * m_unique_cols) + (*j_mapped)] = *v;
605        }
606        Matrix {
607            nrows: n_unique_rows,
608            ncols: m_unique_cols,
609            val: val,
610        }
611    }
612}
613
614/// Function converting a vector of groups into the smallest number for each group shadowing all previous numbers.
615/// 
616/// A group is defined by **consecutive** numbers sharing the same value.
617/// 
618/// # Example:
619/// [3, 3, 2, 3] converts to [1, 1, 3, 6]
620/// 
621/// The first group consisting of two 3s gets the smallest value 1. The second group, consisting of one 2,
622/// gets the smallest value shadowing the first group, i.e. $2\cdot1 + 1 = 3$. The last group, consisting of one 3, gets the smallest
623/// value shadowing all the previous groups, i.e. $2\cdot1 + 1\cdot3 + 1 = 6$.  
624pub fn optimized_bit_allocation_64(v: &Vec<i64>) -> Vec<i64>{
625    let mut res: Vec<i64> = Vec::with_capacity(v.len());
626    let mut grp_size= 0;
627    let mut val_adj = 1;
628    let mut current_val = v[0];
629    for i in 0..v.len(){
630        if v[i] == current_val {
631            grp_size += 1;
632        } else {
633            current_val = v[i];
634            val_adj = val_adj*grp_size + val_adj;
635            grp_size = 1;
636        }
637        res.push(val_adj);
638    }
639    res
640}
641
642#[cfg(test)]
643mod tests {
644    use super::*;
645    #[test]
646    fn test_dot_product() {
647        let m1 = Matrix {
648            val: vec![1.0,2.0,3.0,4.0,5.0,6.0],
649            ncols: 2,
650            nrows: 3,
651        };
652        let m2 = Matrix {
653            val: vec![7.0,8.0,9.0,10.0],
654            ncols: 2,
655            nrows: 2,
656        };
657        assert_eq!(m1.dot(&m2).val, Matrix{val: vec![25.0, 28.0, 57.0, 64.0, 89.0, 100.0], ncols: 2, nrows: 3}.val);
658    }
659    #[test]
660    fn test_optimized_bit_allocation(){
661        assert_eq!(optimized_bit_allocation_64(&vec![1,2,3,4]), vec![1,2,4,8]);
662        assert_eq!(optimized_bit_allocation_64(&vec![1,1,1,2]), vec![1,1,1,4]);
663        assert_eq!(optimized_bit_allocation_64(&vec![2,1,1,1,2]), vec![1,2,2,2,8]);
664        assert_eq!(optimized_bit_allocation_64(&vec![3, 3, 2, 3]), vec![1,1,3,6])
665    }
666
667}
puanrs/linalg.rs

puanrs/
linalg.rs
