Skip to content

Commit c5d21f7

Browse files
committed
Ver 0.39.7
2 parents 09ead45 + 0311bdb commit c5d21f7

12 files changed

Lines changed: 249 additions & 77 deletions

File tree

Cargo.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
[package]
22
name = "peroxide"
3-
version = "0.39.6"
3+
version = "0.39.7"
44
authors = ["axect <axect@outlook.kr>"]
55
edition = "2018"
66
description = "Rust comprehensive scientific computation library contains linear algebra, numerical analysis, statistics and machine learning tools with farmiliar syntax"

RELEASES.md

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,13 @@
1+
# Release 0.39.7 (2025-05-27)
2+
3+
- Add some methods for `DataFrame`
4+
- `filter_by<F: Fn(Scalar) -> bool>(&self, column: &str, f: F) -> anyhow::Result<DataFrame>`
5+
: Filter rows by a condition on a specific column
6+
- `mask(&self, mask: &Series) -> anyhow::Result<DataFrame>`
7+
: Filter rows by a boolean mask
8+
- `select_rows(&self, indices: &[usize]) -> DataFrame`
9+
: Select specific rows by indices
10+
111
# Release 0.39.6 (2025-05-16)
212

313
- New ODESolver: `RKF78`

examples/ode_test_orbit.rs

Lines changed: 19 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,8 @@
11
use peroxide::fuga::*;
22

3-
pub const MU: f64 = 398600.4418; // Standard gravitational parameter of Earth
4-
pub const R_EARTH: f64 = 6378.137; // Radius of Earth in km
5-
pub const J2: f64 = 1.08262668e-3; // J2 coefficient of Earth
3+
pub const MU: f64 = 398600.4418; // Standard gravitational parameter of Earth
4+
pub const R_EARTH: f64 = 6378.137; // Radius of Earth in km
5+
pub const J2: f64 = 1.08262668e-3; // J2 coefficient of Earth
66

77
fn main() -> Result<(), Box<dyn std::error::Error>> {
88
let selected_orbit = OrbitType::Molniya.create_orbit();
@@ -25,18 +25,8 @@ fn main() -> Result<(), Box<dyn std::error::Error>> {
2525

2626
let y0 = Vec::from(initial_state);
2727
y0.print();
28-
let (t, y_gl4) = gl4_solver.solve(
29-
&problem,
30-
(t0, tf),
31-
dt,
32-
&y0,
33-
)?;
34-
let (_, y_rk4) = rk4_solver.solve(
35-
&problem,
36-
(t0, tf),
37-
dt,
38-
&y0,
39-
)?;
28+
let (t, y_gl4) = gl4_solver.solve(&problem, (t0, tf), dt, &y0)?;
29+
let (_, y_rk4) = rk4_solver.solve(&problem, (t0, tf), dt, &y0)?;
4030

4131
let y_gl4 = py_matrix(y_gl4);
4232
let y_rk4 = py_matrix(y_rk4);
@@ -127,7 +117,8 @@ impl ToString for OrbitType {
127117
OrbitType::LEO => "LEO",
128118
OrbitType::GEO => "GEO",
129119
OrbitType::Molniya => "Molniya",
130-
}.to_string()
120+
}
121+
.to_string()
131122
}
132123
}
133124

@@ -157,18 +148,18 @@ impl OrbitType {
157148
raan: 0.0,
158149
w: 270.0f64.to_radians(),
159150
ta: 0.0,
160-
}
151+
},
161152
}
162153
}
163154
}
164155

165156
pub struct Orbit {
166-
pub a: f64, // Semi-major axis
167-
pub e: f64, // Eccentricity
168-
pub i: f64, // Inclination
169-
pub raan: f64, // Right ascension of ascending node
170-
pub w: f64, // Argument of perigee
171-
pub ta: f64, // True anomaly
157+
pub a: f64, // Semi-major axis
158+
pub e: f64, // Eccentricity
159+
pub i: f64, // Inclination
160+
pub raan: f64, // Right ascension of ascending node
161+
pub w: f64, // Argument of perigee
162+
pub ta: f64, // True anomaly
172163
}
173164

174165
impl Orbit {
@@ -178,17 +169,13 @@ impl Orbit {
178169

179170
#[allow(non_snake_case)]
180171
pub fn initial_state(&self) -> State {
181-
let r_pf = vec![
182-
self.r() * self.ta.cos(),
183-
self.r() * self.ta.sin(),
184-
0f64
185-
];
172+
let r_pf = vec![self.r() * self.ta.cos(), self.r() * self.ta.sin(), 0f64];
186173

187174
let p_orbit = self.a * (1.0 - self.e.powi(2));
188175
let v_pf = vec![
189-
- (MU / p_orbit).sqrt() * self.ta.sin(),
176+
-(MU / p_orbit).sqrt() * self.ta.sin(),
190177
(MU / p_orbit).sqrt() * (self.e + self.ta.cos()),
191-
0f64
178+
0f64,
192179
];
193180

194181
let Q = perifocal_to_eci_matrix(&self);
@@ -208,20 +195,12 @@ impl Orbit {
208195

209196
pub fn rot_x(theta: f64) -> Matrix {
210197
let (s, c) = theta.sin_cos();
211-
matrix(vec![
212-
1f64, 0f64, 0f64,
213-
0f64, c, -s,
214-
0f64, s, c
215-
], 3, 3, Row)
198+
matrix(vec![1f64, 0f64, 0f64, 0f64, c, -s, 0f64, s, c], 3, 3, Row)
216199
}
217200

218201
pub fn rot_z(theta: f64) -> Matrix {
219202
let (s, c) = theta.sin_cos();
220-
matrix(vec![
221-
c, -s, 0f64,
222-
s, c, 0f64,
223-
0f64, 0f64, 1f64
224-
], 3, 3, Row)
203+
matrix(vec![c, -s, 0f64, s, c, 0f64, 0f64, 0f64, 1f64], 3, 3, Row)
225204
}
226205

227206
#[allow(non_snake_case)]

src/numerical/ode.rs

Lines changed: 82 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -86,10 +86,10 @@
8686
//! }
8787
//! ```
8888
89-
use anyhow::{bail, Result};
9089
use crate::fuga::ConvToMat;
90+
use crate::traits::math::{InnerProduct, Norm, Normed, Vector};
9191
use crate::util::non_macro::eye;
92-
use crate::traits::math::{Norm, Normed, InnerProduct, Vector};
92+
use anyhow::{bail, Result};
9393

9494
/// Trait for defining an ODE problem.
9595
///
@@ -367,7 +367,10 @@ impl<BU: ButcherTableau> ODEIntegrator for BU {
367367
/// In MATLAB, it is called `ode3`.
368368
#[derive(Debug, Clone, Copy, Default)]
369369
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
370-
#[cfg_attr(feature = "rkyv", derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize))]
370+
#[cfg_attr(
371+
feature = "rkyv",
372+
derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
373+
)]
371374
pub struct RALS3;
372375

373376
impl ButcherTableau for RALS3 {
@@ -383,7 +386,10 @@ impl ButcherTableau for RALS3 {
383386
/// It calculates four intermediate values (k1, k2, k3, k4) to estimate the next step solution.
384387
#[derive(Debug, Clone, Copy, Default)]
385388
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
386-
#[cfg_attr(feature = "rkyv", derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize))]
389+
#[cfg_attr(
390+
feature = "rkyv",
391+
derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
392+
)]
387393
pub struct RK4;
388394

389395
impl ButcherTableau for RK4 {
@@ -398,7 +404,10 @@ impl ButcherTableau for RK4 {
398404
/// This fourth order method is known as minimum truncation error RK4.
399405
#[derive(Debug, Clone, Copy, Default)]
400406
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
401-
#[cfg_attr(feature = "rkyv", derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize))]
407+
#[cfg_attr(
408+
feature = "rkyv",
409+
derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
410+
)]
402411
pub struct RALS4;
403412

404413
impl ButcherTableau for RALS4 {
@@ -418,7 +427,10 @@ impl ButcherTableau for RALS4 {
418427
/// This integrator uses the 5th order Runge-Kutta method to numerically integrate the ODE system.
419428
#[derive(Debug, Clone, Copy, Default)]
420429
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
421-
#[cfg_attr(feature = "rkyv", derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize))]
430+
#[cfg_attr(
431+
feature = "rkyv",
432+
derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
433+
)]
422434
pub struct RK5;
423435

424436
impl ButcherTableau for RK5 {
@@ -478,7 +490,10 @@ impl ButcherTableau for RK5 {
478490
/// - `max_step_iter`: The maximum number of iterations per step.
479491
#[derive(Debug, Clone, Copy)]
480492
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
481-
#[cfg_attr(feature = "rkyv", derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize))]
493+
#[cfg_attr(
494+
feature = "rkyv",
495+
derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
496+
)]
482497
pub struct BS23 {
483498
pub tol: f64,
484499
pub safety_factor: f64,
@@ -560,7 +575,10 @@ impl ButcherTableau for BS23 {
560575
/// - `max_step_iter`: The maximum number of iterations per step.
561576
#[derive(Debug, Clone, Copy)]
562577
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
563-
#[cfg_attr(feature = "rkyv", derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize))]
578+
#[cfg_attr(
579+
feature = "rkyv",
580+
derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
581+
)]
564582
pub struct RKF45 {
565583
pub tol: f64,
566584
pub safety_factor: f64,
@@ -663,7 +681,10 @@ impl ButcherTableau for RKF45 {
663681
/// - `max_step_iter`: The maximum number of iterations per step.
664682
#[derive(Debug, Clone, Copy)]
665683
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
666-
#[cfg_attr(feature = "rkyv", derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize))]
684+
#[cfg_attr(
685+
feature = "rkyv",
686+
derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
687+
)]
667688
pub struct DP45 {
668689
pub tol: f64,
669690
pub safety_factor: f64,
@@ -785,7 +806,10 @@ impl ButcherTableau for DP45 {
785806
/// - Ch. Tsitouras, Comput. Math. Appl. 62 (2011) 770-780.
786807
#[derive(Debug, Clone, Copy)]
787808
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
788-
#[cfg_attr(feature = "rkyv", derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize))]
809+
#[cfg_attr(
810+
feature = "rkyv",
811+
derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
812+
)]
789813
pub struct TSIT45 {
790814
pub tol: f64,
791815
pub safety_factor: f64,
@@ -916,7 +940,10 @@ impl ButcherTableau for TSIT45 {
916940
/// - Meysam Mahooti (2025). [Runge-Kutta-Fehlberg (RKF78)](https://www.mathworks.com/matlabcentral/fileexchange/61130-runge-kutta-fehlberg-rkf78), MATLAB Central File Exchange.
917941
#[derive(Debug, Clone, Copy)]
918942
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
919-
#[cfg_attr(feature = "rkyv", derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize))]
943+
#[cfg_attr(
944+
feature = "rkyv",
945+
derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
946+
)]
920947
pub struct RKF78 {
921948
pub tol: f64,
922949
pub safety_factor: f64,
@@ -986,11 +1013,35 @@ impl ButcherTableau for RKF78 {
9861013
// k6
9871014
&[1.0 / 20.0, 0.0, 0.0, 5.0 / 20.0, 4.0 / 20.0],
9881015
// k7
989-
&[-25.0 / 108.0, 0.0, 0.0, 125.0 / 108.0, -260.0 / 108.0, 250.0 / 108.0],
1016+
&[
1017+
-25.0 / 108.0,
1018+
0.0,
1019+
0.0,
1020+
125.0 / 108.0,
1021+
-260.0 / 108.0,
1022+
250.0 / 108.0,
1023+
],
9901024
// k8
991-
&[31.0 / 300.0, 0.0, 0.0, 0.0, 61.0 / 225.0, -2.0 / 9.0, 13.0 / 900.0],
1025+
&[
1026+
31.0 / 300.0,
1027+
0.0,
1028+
0.0,
1029+
0.0,
1030+
61.0 / 225.0,
1031+
-2.0 / 9.0,
1032+
13.0 / 900.0,
1033+
],
9921034
// k9
993-
&[2.0, 0.0, 0.0, -53.0 / 6.0, 704.0 / 45.0, -107.0 / 9.0, 67.0 / 90.0, 3.0],
1035+
&[
1036+
2.0,
1037+
0.0,
1038+
0.0,
1039+
-53.0 / 6.0,
1040+
704.0 / 45.0,
1041+
-107.0 / 9.0,
1042+
67.0 / 90.0,
1043+
3.0,
1044+
],
9941045
// k10
9951046
&[
9961047
-91.0 / 108.0,
@@ -1051,7 +1102,7 @@ impl ButcherTableau for RKF78 {
10511102
// BU_i = BE_i (8th order) - ErrorCoeff_i
10521103
// ErrorCoeff_i = [-41/840, 0, ..., 0, -41/840 (for k11), 41/840 (for k12), 41/840 (for k13)]
10531104
const BU: &'static [f64] = &[
1054-
41.0 / 420.0, // 41/840 - (-41/840)
1105+
41.0 / 420.0, // 41/840 - (-41/840)
10551106
0.0,
10561107
0.0,
10571108
0.0,
@@ -1110,7 +1161,10 @@ impl ButcherTableau for RKF78 {
11101161
/// Currently, there are two options: fixed-point iteration and Broyden's method.
11111162
#[derive(Debug, Clone, Copy)]
11121163
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
1113-
#[cfg_attr(feature = "rkyv", derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize))]
1164+
#[cfg_attr(
1165+
feature = "rkyv",
1166+
derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
1167+
)]
11141168
pub enum ImplicitSolver {
11151169
FixedPoint,
11161170
Broyden,
@@ -1130,7 +1184,10 @@ pub enum ImplicitSolver {
11301184
/// - `max_step_iter`: The maximum number of iterations for the implicit solver per step.
11311185
#[derive(Debug, Clone, Copy)]
11321186
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
1133-
#[cfg_attr(feature = "rkyv", derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize))]
1187+
#[cfg_attr(
1188+
feature = "rkyv",
1189+
derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
1190+
)]
11341191
pub struct GL4 {
11351192
pub solver: ImplicitSolver,
11361193
pub tol: f64,
@@ -1203,7 +1260,7 @@ impl ODEIntegrator for GL4 {
12031260
break;
12041261
}
12051262
}
1206-
},
1263+
}
12071264
ImplicitSolver::Broyden => {
12081265
let m = 2 * n;
12091266
let mut U = vec![0.0; m];
@@ -1219,7 +1276,7 @@ impl ODEIntegrator for GL4 {
12191276
y2.copy_from_slice(y);
12201277
problem.rhs(t + c * dt, &y1, &mut k1)?;
12211278
problem.rhs(t + d * dt, &y2, &mut k2)?;
1222-
for i in 0 .. n {
1279+
for i in 0..n {
12231280
U[i] = k1[i];
12241281
U[n + i] = k2[i];
12251282
}
@@ -1233,12 +1290,14 @@ impl ODEIntegrator for GL4 {
12331290
let mut J_inv = eye(m);
12341291

12351292
// Repeat Broyden's method
1236-
for _ in 0 .. self.max_step_iter {
1293+
for _ in 0..self.max_step_iter {
12371294
// delta = - J_inv * F_vec
12381295
let mut delta = (&J_inv * &F_vec).mul_scalar(-1.0);
12391296

12401297
// U <- U + delta
1241-
U.iter_mut().zip(delta.iter_mut()).for_each(|(u, d)| *u += *d);
1298+
U.iter_mut()
1299+
.zip(delta.iter_mut())
1300+
.for_each(|(u, d)| *u += *d);
12421301

12431302
let mut F_new = vec![0.0; m];
12441303
compute_F(problem, t, y, dt, c, d, sqrt3, &U, &mut F_new)?;
@@ -1268,11 +1327,11 @@ impl ODEIntegrator for GL4 {
12681327
}
12691328

12701329
// Extract k1 and k2
1271-
for i in 0 .. n {
1330+
for i in 0..n {
12721331
k1[i] = U[i];
12731332
k2[i] = U[n + i];
12741333
}
1275-
},
1334+
}
12761335
}
12771336

12781337
for i in 0..n {

0 commit comments

Comments
 (0)