1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
use std::ops::{Index, IndexMut};
use super::tensors::{Matrix, Tensor, CovariantIndex, ContravariantIndex};
use typenum::Pow;
use typenum::uint::Unsigned;
use typenum::consts::U2;
use generic_array::{GenericArray, ArrayLength};
pub trait CoordinateSystem : Sized {
type Dimension: Unsigned + ArrayLength<f64> + ArrayLength<usize>;
fn small(_: &Point<Self>) -> f64 {
0.01
}
fn dimension() -> usize {
Self::Dimension::to_usize()
}
}
pub struct Point<T: CoordinateSystem> {
x: GenericArray<f64, T::Dimension>,
}
impl<T> Point<T> where T: CoordinateSystem
{
pub fn new(coords: GenericArray<f64, T::Dimension>) -> Point<T> {
Point { x: coords }
}
pub fn from_slice(coords: &[f64]) -> Point<T> {
Point { x: GenericArray::from_slice(coords) }
}
}
impl<T> Clone for Point<T> where T: CoordinateSystem
{
fn clone(&self) -> Point<T> {
Point::new(self.x.clone())
}
}
impl<T> Copy for Point<T> where T: CoordinateSystem, <T::Dimension as ArrayLength<f64>>::ArrayType: Copy {}
impl<T> Index<usize> for Point<T> where T: CoordinateSystem
{
type Output = f64;
fn index(&self, idx: usize) -> &f64 {
&self.x[idx]
}
}
impl<T> IndexMut<usize> for Point<T> where T: CoordinateSystem
{
fn index_mut(&mut self, idx: usize) -> &mut f64 {
&mut self.x[idx]
}
}
impl<T> PartialEq<Point<T>> for Point<T>
where T: CoordinateSystem
{
fn eq(&self, rhs: &Point<T>) -> bool {
(0..T::dimension()).all(|i| self[i] == rhs[i])
}
}
impl<T> Eq for Point<T> where T: CoordinateSystem {}
pub trait ConversionTo<T: CoordinateSystem + 'static> : CoordinateSystem
where T::Dimension: Pow<U2>,
<T::Dimension as Pow<U2>>::Output: ArrayLength<f64>
{
fn convert_point(p: &Point<Self>) -> Point<T>;
fn jacobian(p: &Point<Self>) -> Matrix<T> {
let d = Self::dimension();
let mut result = Matrix::new(Self::convert_point(p));
let h = Self::small(p);
for j in 0..d {
let mut x = p.clone();
x[j] = x[j] - h;
let y1 = Self::convert_point(&x);
x[j] = x[j] + h * 2.0;
let y2 = Self::convert_point(&x);
for i in 0..d {
let index = [i, j];
result[&index[..]] = (y2[i] - y1[i]) / (2.0 * h);
}
}
result
}
fn inv_jacobian(p: &Point<Self>) -> Tensor<T, (CovariantIndex, ContravariantIndex)> {
ConversionTo::<T>::jacobian(p).inverse().unwrap()
}
}