fastsim_core/
drive_cycle.rs

1pub mod maneuvers;
2pub mod manipulation_utils;
3
4use crate::drive_cycle::manipulation_utils::{
5    speed_for_constant_jerk, ConstantJerkTrajectory, CycleCache,
6};
7use crate::imports::*;
8use crate::prelude::*;
9use fastsim_2::cycle::RustCycle as Cycle2;
10use std::cmp;
11
12#[serde_api]
13#[derive(Clone, Debug, Deserialize, Serialize, PartialEq, Default)]
14#[non_exhaustive]
15#[serde(deny_unknown_fields)]
16#[cfg_attr(feature = "pyo3", pyclass(module = "fastsim", subclass, eq))]
17/// Container
18pub struct Cycle {
19    /// Name of cycle (can be left empty)
20    #[serde(default, skip_serializing_if = "String::is_empty")]
21    pub name: String,
22    /// inital elevation
23    pub init_elev: Option<si::Length>,
24    /// simulation time
25    pub time: Vec<si::Time>,
26    /// prescribed speed
27    #[serde(alias = "speed_mps")]
28    pub speed: Vec<si::Velocity>,
29    // TODO: consider trapezoidal integration scheme
30    /// calculated prescribed distance based on RHS integral of time and speed
31    #[serde(default, skip_serializing_if = "Vec::is_empty")]
32    pub dist: Vec<si::Length>,
33    /// road grade (expressed as a decimal, not percent)
34    #[serde(default, skip_serializing_if = "Vec::is_empty")]
35    pub grade: Vec<si::Ratio>,
36    // TODO: consider trapezoidal integration scheme
37    // TODO: @mokeefe, please check out how elevation is handled
38    /// calculated prescribed elevation based on RHS integral distance and grade
39    #[serde(default, skip_serializing_if = "Vec::is_empty")]
40    pub elev: Vec<si::Length>,
41    /// road charging/discharing capacity
42    #[serde(default, skip_serializing_if = "Vec::is_empty")]
43    pub pwr_max_chrg: Vec<si::Power>,
44    /// ambient air temperature w.r.t. to time (rather than spatial position)
45    #[serde(default, skip_serializing_if = "Vec::is_empty")]
46    pub temp_amb_air: Vec<si::Temperature>,
47    /// solar heat load w.r.t. to time (rather than spatial position)
48    #[serde(default, skip_serializing_if = "Vec::is_empty")]
49    pub pwr_solar_load: Vec<si::Power>,
50    // TODO: add provision for optional time-varying aux load
51    /// grade interpolator
52    #[serde(default, skip_serializing_if = "Option::is_none")]
53    pub grade_interp: Option<InterpolatorEnumOwned<f64>>,
54    /// elevation interpolator
55    #[serde(default, skip_serializing_if = "Option::is_none")]
56    pub elev_interp: Option<InterpolatorEnumOwned<f64>>,
57}
58
59#[pyo3_api]
60impl Cycle {
61    #[pyo3(name = "len")]
62    /// return the length of the cycle
63    fn len_py(&self) -> PyResult<usize> {
64        Ok(self.len_checked()?)
65    }
66
67    #[pyo3(name = "to_microtrips", signature=(stop_speed_m_per_s=None))]
68    /// convert cycle to a list of microtrips.
69    /// If stop speed is specified, it signifies the speed at or below which
70    /// a vehicle should be considered as stopped. This can be useful when
71    /// processing real-world data.
72    fn to_microtrips_py(&self, stop_speed_m_per_s: Option<f64>) -> PyResult<Vec<Cycle>> {
73        let stop_speed = stop_speed_m_per_s.map(|v| v * uc::MPS);
74        Ok(self.to_microtrips(stop_speed))
75    }
76
77    #[pyo3(name = "extend_time", signature=(absolute_time_s=None, time_fraction=None))]
78    /// extend cycle with idle time.
79    /// This is useful when a cycle's duration needs to be extended.
80    /// - absolute_time_s: optional time to extend the cycle
81    /// - time_fraction: optional fraction of cycle's duration to add to cycle.
82    ///
83    /// NOTE: if both absolute time and time fraction are specified, they
84    /// both add to extend the cycle. For example, if we have a 100 s cycle
85    /// and specify an absolute_time_s of 10 and time_fraction of 0.5, the
86    /// resulting cycle will have a duration of 160 s = 100.0 + (10 + 100.0 * 0.5)
87    fn extend_time_py(
88        &mut self,
89        absolute_time_s: Option<f64>,
90        time_fraction: Option<f64>,
91    ) -> PyResult<Cycle> {
92        let absolute_time = absolute_time_s.map(|t| t * uc::S);
93        let time_fraction = time_fraction.map(|f| f * uc::R);
94        Ok(self.extend_time(absolute_time, time_fraction))
95    }
96
97    #[pyo3(name = "dt_at_i")]
98    /// time step duration at step i.
99    pub fn dt_at_i_py(&self, i: usize) -> PyResult<f64> {
100        let i = std::cmp::max(1, i);
101        let dt = if i < self.time.len() {
102            self.time[i].get::<si::second>() - self.time[i - 1].get::<si::second>()
103        } else {
104            0.0
105        };
106        Ok(dt)
107    }
108
109    #[pyo3(name = "ending_idle_time_s")]
110    /// calculate and return the ending "idle" time of a cycle.
111    /// "Idle" time is defined as the amount of contiguous time
112    /// at the end of a cycle where the vehicle is not moving.
113    pub fn ending_idle_time_py(&self) -> PyResult<f64> {
114        let dt_end_idle = self.ending_idle_time();
115        Ok(dt_end_idle.get::<si::second>())
116    }
117
118    #[pyo3(name = "trim_ending_idle", signature=(idle_to_keep_s=None))]
119    /// trim ending "idle" time from a cycle.
120    /// The "idle" time is the time the vehicle is not moving.
121    /// - idle_to_keep_s: the amount of time to keep
122    ///
123    /// NOTE: if idle_to_keep_s is specified, the ending idle duration
124    /// will be UP TO this idle_to_keep_s amount but could be less if
125    /// there is insufficient idle time.
126    pub fn trim_ending_idle_py(&self, idle_to_keep_s: Option<f64>) -> PyResult<Cycle> {
127        let idle_to_keep = idle_to_keep_s.map(|idle| idle * uc::S);
128        Ok(self.trim_ending_idle(idle_to_keep))
129    }
130
131    #[pyo3(name = "average_speed_m_per_s", signature=(while_moving=None))]
132    /// calculate and return the average speed of the cycle in (m/s).
133    /// - while_moving: if specified and true, calculate the speed only
134    ///   while the vehicle is moving. Otherwise, calculate the average speed
135    ///   including stopped time.
136    pub fn average_speed_py(&self, while_moving: Option<bool>) -> PyResult<f64> {
137        let while_moving = while_moving.unwrap_or(false);
138        let vavg = self.average_speed(while_moving);
139        Ok(vavg.get::<si::meter_per_second>())
140    }
141
142    #[pyo3(name = "average_step_speeds_m_per_s")]
143    /// calculate and return the average speeds per time-step in (m/s).
144    pub fn average_step_speeds_py(&self) -> PyResult<Vec<f64>> {
145        Ok(self
146            .average_step_speeds()
147            .iter()
148            .map(|v| v.get::<si::meter_per_second>())
149            .collect())
150    }
151
152    #[pyo3(name = "average_step_speed_in_m_per_s_at")]
153    /// calculate the average step speed at the given step in (m/s).
154    pub fn average_step_speed_at_py(&self, i: usize) -> PyResult<f64> {
155        Ok(self.average_step_speed_at(i).get::<si::meter_per_second>())
156    }
157
158    #[pyo3(name = "resample")]
159    /// create a new cycle with the values resampled to the given time-step
160    /// duration.
161    pub fn resample_py(&self, time_step_s: f64) -> PyResult<Cycle> {
162        let time_step = time_step_s.max(0.01) * uc::S;
163        Ok(self.resample(time_step))
164    }
165}
166
167lazy_static! {
168    pub static ref ELEV_DEFAULT: si::Length = 400. * uc::FT;
169}
170
171impl Init for Cycle {
172    /// Sets `self.dist` and `self.elev`
173    /// # Assumptions
174    /// - if `init_elev.is_none()`, then defaults to [static@ELEV_DEFAULT]
175    fn init(&mut self) -> Result<(), Error> {
176        let _ = self
177            .len_checked()
178            .map_err(|err| Error::InitError(format_dbg!(err)))?;
179
180        if !self.temp_amb_air.is_empty() {
181            if self.temp_amb_air.len() != self.time.len() {
182                return Err(Error::InitError(format_dbg!()));
183            }
184        } else {
185            self.temp_amb_air = vec![*TE_STD_AIR; self.time.len()];
186        }
187
188        // calculate distance from RHS integral of speed and time
189        self.dist = {
190            self.time
191                .diff()
192                .iter()
193                .zip(&self.speed)
194                .scan(0. * uc::M, |dist, (dt, speed)| {
195                    *dist += *dt * *speed;
196                    Some(*dist)
197                })
198                .collect()
199        };
200
201        // populate grade if not provided
202        if self.grade.is_empty() {
203            self.grade = vec![
204                si::Ratio::ZERO;
205                self.len_checked()
206                    .map_err(|err| Error::InitError(format_dbg!(err)))?
207            ]
208        };
209        // calculate elevation from RHS integral of grade and distance
210        self.init_elev = self.init_elev.or_else(|| Some(*ELEV_DEFAULT));
211        self.elev = self
212            .grade
213            .iter()
214            .zip(&self.dist.diff())
215            .scan(
216                // already guaranteed to be `Some`
217                self.init_elev.unwrap(),
218                |elev, (grade, dist)| {
219                    *elev += *dist * grade.atan().sin();
220                    Some(*elev)
221                },
222            )
223            .collect();
224        let g0 = if !self.grade.is_empty() {
225            self.grade[0]
226        } else {
227            0.0 * uc::R
228        };
229        if self.grade.iter().all(|&g| g != g0) {
230            self.grade_interp = Some(
231                InterpolatorEnum::new_1d(
232                    self.dist.iter().map(|x| x.get::<si::meter>()).collect(),
233                    self.grade.iter().map(|y| y.get::<si::ratio>()).collect(),
234                    strategy::Linear,
235                    Extrapolate::Error,
236                )
237                .map_err(|e| Error::NinterpError(e.to_string()))?,
238            );
239
240            self.elev_interp = Some(
241                InterpolatorEnum::new_1d(
242                    self.dist.iter().map(|x| x.get::<si::meter>()).collect(),
243                    self.elev.iter().map(|y| y.get::<si::meter>()).collect(),
244                    strategy::Linear,
245                    Extrapolate::Error,
246                )
247                .map_err(|e| Error::NinterpError(e.to_string()))?,
248            );
249        } else {
250            self.grade_interp = Some(InterpolatorEnum::new_0d(g0.get::<si::ratio>()));
251            self.elev_interp = Some(InterpolatorEnum::new_0d(
252                self.init_elev.unwrap().get::<si::meter>(),
253            ));
254        }
255
256        Ok(())
257    }
258}
259
260impl SerdeAPI for Cycle {
261    const ACCEPTED_BYTE_FORMATS: &'static [&'static str] = &[
262        #[cfg(feature = "csv")]
263        "csv",
264        #[cfg(feature = "json")]
265        "json",
266        #[cfg(feature = "msgpack")]
267        "msgpack",
268        #[cfg(feature = "toml")]
269        "toml",
270        #[cfg(feature = "yaml")]
271        "yaml",
272    ];
273    const ACCEPTED_STR_FORMATS: &'static [&'static str] = &[
274        #[cfg(feature = "csv")]
275        "csv",
276        #[cfg(feature = "json")]
277        "json",
278        #[cfg(feature = "toml")]
279        "toml",
280        #[cfg(feature = "yaml")]
281        "yaml",
282    ];
283    #[cfg(feature = "resources")]
284    const RESOURCES_SUBDIR: &'static str = "cycles";
285
286    /// Write (serialize) an object into anything that implements [`std::io::Write`]
287    ///
288    /// # Arguments:
289    ///
290    /// * `wtr` - The writer into which to write object data
291    /// * `format` - The target format, any of those listed in [`ACCEPTED_BYTE_FORMATS`](`SerdeAPI::ACCEPTED_BYTE_FORMATS`)
292    ///
293    fn to_writer<W: std::io::Write>(&self, mut wtr: W, format: &str) -> Result<(), Error> {
294        match format.trim_start_matches('.').to_lowercase().as_str() {
295            #[cfg(feature = "csv")]
296            "csv" => {
297                let mut wtr = csv::Writer::from_writer(wtr);
298                for i in 0..self
299                    .len_checked()
300                    .map_err(|err| Error::SerdeError(format_dbg!(err)))?
301                {
302                    wtr.serialize(CycleElement {
303                        // unchecked indexing should be ok because of `self.len()`
304                        time: self.time[i],
305                        speed: self.speed[i],
306                        grade: if !self.grade.is_empty() {
307                            Some(self.grade[i])
308                        } else {
309                            None
310                        },
311                        pwr_max_charge: if !self.pwr_max_chrg.is_empty() {
312                            Some(self.pwr_max_chrg[i])
313                        } else {
314                            None
315                        },
316                        temp_amb_air: if !self.temp_amb_air.is_empty() {
317                            Some(self.temp_amb_air[i])
318                        } else {
319                            None
320                        },
321                        pwr_solar_load: if !self.pwr_solar_load.is_empty() {
322                            Some(self.pwr_solar_load[i])
323                        } else {
324                            None
325                        },
326                    })
327                    .map_err(|err| Error::SerdeError(format_dbg!(err)))?;
328                }
329                wtr.flush()
330                    .map_err(|err| Error::SerdeError(format_dbg!(err)))?
331            }
332            #[cfg(feature = "json")]
333            "json" => serde_json::to_writer(wtr, self)
334                .map_err(|err| Error::SerdeError(format_dbg!(err)))?,
335            #[cfg(feature = "toml")]
336            "toml" => {
337                let toml_string = self
338                    .to_toml()
339                    .map_err(|err| Error::SerdeError(format_dbg!(err)))?;
340                wtr.write_all(toml_string.as_bytes())
341                    .map_err(|err| Error::SerdeError(format_dbg!(err)))?;
342            }
343            #[cfg(feature = "yaml")]
344            "yaml" | "yml" => serde_yaml::to_writer(wtr, self)
345                .map_err(|err| Error::SerdeError(format_dbg!(err)))?,
346            _ => Err(Error::SerdeError(format!(
347                "Unsupported format {format:?}, must be one of {:?}",
348                Self::ACCEPTED_BYTE_FORMATS,
349            )))?,
350        }
351        Ok(())
352    }
353
354    /// Deserialize an object from anything that implements [`std::io::Read`]
355    ///
356    /// # Arguments:
357    ///
358    /// * `rdr` - The reader from which to read object data
359    /// * `format` - The source format, any of those listed in [`ACCEPTED_BYTE_FORMATS`](`SerdeAPI::ACCEPTED_BYTE_FORMATS`)
360    ///
361    fn from_reader<R: std::io::Read>(
362        rdr: &mut R,
363        format: &str,
364        skip_init: bool,
365    ) -> Result<Self, Error> {
366        let mut deserialized: Self =
367            match format.trim_start_matches('.').to_lowercase().as_str() {
368                #[cfg(feature = "csv")]
369                "csv" => {
370                    // Create empty cycle to be populated
371                    let mut cyc = Self::default();
372                    let mut rdr = csv::Reader::from_reader(rdr);
373                    for result in rdr.deserialize() {
374                        cyc.push(result.map_err(|err| Error::SerdeError(format_dbg!(err)))?)
375                            .map_err(|err| Error::SerdeError(format!("{err}")))?;
376                    }
377                    cyc
378                }
379                #[cfg(feature = "json")]
380                "json" => serde_json::from_reader(rdr)
381                    .map_err(|err| Error::SerdeError(format!("{err}")))?,
382                #[cfg(feature = "toml")]
383                "toml" => {
384                    let mut buf = String::new();
385                    rdr.read_to_string(&mut buf)
386                        .map_err(|err| Error::SerdeError(format_dbg!(err)))?;
387                    Self::from_toml(buf, skip_init)
388                        .map_err(|err| Error::SerdeError(format_dbg!(err)))?
389                }
390                #[cfg(feature = "yaml")]
391                "yaml" | "yml" => serde_yaml::from_reader(rdr)
392                    .map_err(|err| Error::SerdeError(format_dbg!(err)))?,
393                _ => {
394                    return Err(Error::SerdeError(format!(
395                        "Unsupported format {format:?}, must be one of {:?}",
396                        Self::ACCEPTED_BYTE_FORMATS
397                    )))
398                }
399            };
400        if !skip_init {
401            deserialized.init()?;
402        }
403        Ok(deserialized)
404    }
405
406    /// Write (serialize) an object into a string
407    ///
408    /// # Arguments:
409    ///
410    /// * `format` - The target format, any of those listed in [`ACCEPTED_STR_FORMATS`](`SerdeAPI::ACCEPTED_STR_FORMATS`)
411    ///
412    fn to_str(&self, format: &str) -> anyhow::Result<String> {
413        match format.trim_start_matches('.').to_lowercase().as_str() {
414            #[cfg(feature = "csv")]
415            "csv" => self.to_csv(),
416            #[cfg(feature = "json")]
417            "json" => self.to_json(),
418            #[cfg(feature = "toml")]
419            "toml" => self.to_toml(),
420            #[cfg(feature = "yaml")]
421            "yaml" | "yml" => self.to_yaml(),
422            _ => bail!(
423                "Unsupported format {format:?}, must be one of {:?}",
424                Self::ACCEPTED_STR_FORMATS
425            ),
426        }
427    }
428
429    /// Read (deserialize) an object from a string
430    ///
431    /// # Arguments:
432    ///
433    /// * `contents` - The string containing the object data
434    /// * `format` - The source format, any of those listed in [`ACCEPTED_STR_FORMATS`](`SerdeAPI::ACCEPTED_STR_FORMATS`)
435    ///
436    fn from_str<S: AsRef<str>>(contents: S, format: &str, skip_init: bool) -> anyhow::Result<Self> {
437        Ok(
438            match format.trim_start_matches('.').to_lowercase().as_str() {
439                #[cfg(feature = "csv")]
440                "csv" => Self::from_csv(contents, skip_init)?,
441                #[cfg(feature = "json")]
442                "json" => Self::from_json(contents, skip_init)?,
443                #[cfg(feature = "toml")]
444                "toml" => Self::from_toml(contents, skip_init)?,
445                #[cfg(feature = "yaml")]
446                "yaml" | "yml" => Self::from_yaml(contents, skip_init)?,
447                _ => bail!(
448                    "Unsupported format {format:?}, must be one of {:?}",
449                    Self::ACCEPTED_STR_FORMATS
450                ),
451            },
452        )
453    }
454}
455
456impl Cycle {
457    /// rust-internal time steps at i
458    pub fn dt_at_i(&self, i: usize) -> anyhow::Result<si::Time> {
459        Ok(*self.time.get(i).with_context(|| format_dbg!())?
460            - *self.time.get(i - 1).with_context(|| format_dbg!())?)
461    }
462
463    /// return the length of the cycle
464    pub fn len_checked(&self) -> anyhow::Result<usize> {
465        ensure!(
466            self.time.len() == self.speed.len(),
467            format!(
468                "{}\n`time` and `speed` fields do not have same `len()`",
469                format_dbg!()
470            )
471        );
472        ensure!(
473            self.dist.is_empty() || self.time.len() == self.dist.len(),
474            format!(
475                "{}\n`time` and `dist` fields do not have same `len()`",
476                format_dbg!()
477            )
478        );
479        ensure!(
480            self.grade.is_empty() || self.time.len() == self.grade.len(),
481            format!(
482                "{}\n`time` and `grade` fields do not have same `len()`",
483                format_dbg!()
484            )
485        );
486        ensure!(
487            self.elev.is_empty() || self.grade.len() == self.elev.len(),
488            format!(
489                "{}\n`grade` and `elev` fields do not have same `len()`",
490                format_dbg!()
491            )
492        );
493        ensure!(
494            self.pwr_max_chrg.is_empty() || self.time.len() == self.pwr_max_chrg.len(),
495            format!(
496                "{}\n`time` and `pwr_max_chrg` fields do not have same `len()`",
497                format_dbg!()
498            )
499        );
500        ensure!(
501            self.temp_amb_air.is_empty() || self.time.len() == self.temp_amb_air.len(),
502            format!(
503                "{}\n`time` and `temp_amb_air` fields do not have same `len()`",
504                format_dbg!()
505            )
506        );
507        Ok(self.time.len())
508    }
509
510    /// return true if the cycle is empty, else false
511    pub fn is_empty(&self) -> anyhow::Result<bool> {
512        Ok(self.len_checked().with_context(|| format_dbg!())? == 0)
513    }
514
515    /// append the given cycle element
516    pub fn push(&mut self, element: CycleElement) -> anyhow::Result<()> {
517        // TODO: maybe automate generation of this function as derive macro
518        // TODO: maybe automate `ensure!` that all vec fields are same length before returning result
519        // TODO: make sure all fields are being updated as appropriate
520        self.time.push(element.time);
521        self.speed.push(element.speed);
522        match element.grade {
523            Some(grade) => self.grade.push(grade),
524            None => self.grade.push(si::Ratio::ZERO),
525        }
526        match element.pwr_max_charge {
527            Some(pwr_max_chrg) => self.pwr_max_chrg.push(pwr_max_chrg),
528            None => self.pwr_max_chrg.push(si::Power::ZERO),
529        }
530        match element.temp_amb_air {
531            Some(temp_amb_air) => self.temp_amb_air.push(temp_amb_air),
532            None => self.temp_amb_air.push(*TE_STD_AIR),
533        }
534        match element.pwr_solar_load {
535            Some(pwr_solar_load) => self.pwr_solar_load.push(pwr_solar_load),
536            None => self.pwr_solar_load.push(si::Power::ZERO),
537        }
538        Ok(())
539    }
540
541    /// extend the cycle by a vector of elements
542    pub fn extend(&mut self, vec: Vec<CycleElement>) -> anyhow::Result<()> {
543        self.time.extend(vec.iter().map(|x| x.time).clone());
544        todo!();
545        // self.time.extend(vec.iter().map(|x| x.time).clone());
546        // match (&mut self.grade, vec.grade) {
547        //     (Some(grade_mut), Some(grade)) => grade_mut.push(grade),
548        //     (None, Some(_)) => {
549        //         bail!("Element and Cycle `grade` fields must both be `Some` or `None`")
550        //     }
551        //     (Some(_), None) => {
552        //         bail!("Element and Cycle `grade` fields must both be `Some` or `None`")
553        //     }
554        //     _ => {}
555        // }
556        // match (&mut self.pwr_max_chrg, vec.pwr_max_charge) {
557        //     (Some(pwr_max_chrg_mut), Some(pwr_max_chrg)) => pwr_max_chrg_mut.push(pwr_max_chrg),
558        //     (None, Some(_)) => {
559        //         bail!("Element and Cycle `pwr_max_chrg` fields must both be `Some` or `None`")
560        //     }
561        //     (Some(_), None) => {
562        //         bail!("Element and Cycle `pwr_max_chrg` fields must both be `Some` or `None`")
563        //     }
564        //     _ => {}
565        // }
566        // self.speed.push(vec.speed);
567        // Ok(())
568    }
569
570    /// trim the cycle to the given start_idx and end_idx.
571    ///
572    /// NOTE: ending cycle will include start_idx but NOT end_idx
573    pub fn trim(&mut self, start_idx: Option<usize>, end_idx: Option<usize>) -> anyhow::Result<()> {
574        let start_idx = start_idx.unwrap_or_default();
575        let len = self.len_checked().with_context(|| format_dbg!())?;
576        let end_idx = end_idx.unwrap_or(len);
577        ensure!(end_idx <= len, format_dbg!(end_idx <= len));
578
579        self.time = self.time[start_idx..end_idx].to_vec();
580        self.speed = self.speed[start_idx..end_idx].to_vec();
581        Ok(())
582    }
583
584    /// Write (serialize) cycle to a CSV string
585    #[cfg(feature = "csv")]
586    pub fn to_csv(&self) -> anyhow::Result<String> {
587        let mut buf = Vec::with_capacity(self.len_checked().with_context(|| format_dbg!())?);
588        self.to_writer(&mut buf, "csv")?;
589        Ok(String::from_utf8(buf)?)
590    }
591
592    /// Read (deserialize) an object from a CSV string
593    ///
594    /// # Arguments
595    ///
596    /// * `json_str` - JSON-formatted string to deserialize from
597    ///
598    #[cfg(feature = "csv")]
599    fn from_csv<S: AsRef<str>>(csv_str: S, skip_init: bool) -> anyhow::Result<Self> {
600        let mut csv_de = Self::from_reader(&mut csv_str.as_ref().as_bytes(), "csv", skip_init)?;
601        if !skip_init {
602            csv_de.init()?;
603        }
604        Ok(csv_de)
605    }
606
607    pub fn to_fastsim2(&self) -> anyhow::Result<Cycle2> {
608        let cyc2 = Cycle2 {
609            name: self.name.clone(),
610            time_s: self.time.iter().map(|t| t.get::<si::second>()).collect(),
611            mps: self
612                .speed
613                .iter()
614                .map(|s| s.get::<si::meter_per_second>())
615                .collect(),
616            grade: self.grade.iter().map(|g| g.get::<si::ratio>()).collect(),
617            orphaned: false,
618            road_type: vec![0.; self.len_checked().with_context(|| format_dbg!())?].into(),
619        };
620
621        Ok(cyc2)
622    }
623
624    /// convert cycle to a vector of CycleElement
625    pub fn to_elements(&self) -> Vec<CycleElement> {
626        let mut result = Vec::with_capacity(self.time.len());
627        for idx in 0..self.time.len() {
628            let element = CycleElement {
629                time: self.time[idx],
630                speed: self.speed[idx],
631                grade: if self.grade.is_empty() {
632                    None
633                } else {
634                    Some(self.grade[idx])
635                },
636                pwr_max_charge: if self.pwr_max_chrg.is_empty() {
637                    None
638                } else {
639                    Some(self.pwr_max_chrg[idx])
640                },
641                temp_amb_air: if self.temp_amb_air.is_empty() {
642                    None
643                } else {
644                    Some(self.temp_amb_air[idx])
645                },
646                pwr_solar_load: if self.pwr_solar_load.is_empty() {
647                    None
648                } else {
649                    Some(self.pwr_solar_load[idx])
650                },
651            };
652            result.push(element);
653        }
654        result
655    }
656
657    /// Convert cycle into a vector of "microtrips".
658    /// A microtrip is a start to a subsequent stop plus any idle time.
659    /// - stop_speed: the speed at or below which vehicle is considered "stopped"
660    ///
661    /// RETURN: vector of cycles with each cycle being a "microtrip".
662    pub fn to_microtrips(&self, stop_speed: Option<si::Velocity>) -> Vec<Cycle> {
663        let stop_speed = stop_speed.unwrap_or(1e-6 * uc::MPS);
664        let mut microtrips = Vec::new();
665        let mut current = Cycle {
666            name: self.name.clone(),
667            init_elev: self.init_elev,
668            time: vec![],
669            speed: vec![],
670            dist: vec![],
671            grade: vec![],
672            elev: vec![],
673            pwr_max_chrg: vec![],
674            temp_amb_air: vec![],
675            pwr_solar_load: vec![],
676            grade_interp: self.grade_interp.clone(),
677            elev_interp: self.elev_interp.clone(),
678        };
679        let elements = self.to_elements();
680        let mut moving: bool = false;
681        for element in &elements {
682            if element.speed > stop_speed && !moving && current.time.len() > 1 {
683                current.init().unwrap();
684                let last_idx = current.time.len() - 1;
685                let last_time = current.time[last_idx];
686                let last_speed = current.speed[last_idx];
687                let last_grade = if last_idx >= current.grade.len() {
688                    None
689                } else {
690                    Some(current.grade[last_idx])
691                };
692                let last_elevation = if last_idx >= current.elev.len() {
693                    None
694                } else {
695                    Some(current.elev[last_idx])
696                };
697                let last_temperature = if last_idx >= current.temp_amb_air.len() {
698                    None
699                } else {
700                    Some(current.temp_amb_air[last_idx])
701                };
702                let last_solar_load = if last_idx >= current.pwr_solar_load.len() {
703                    None
704                } else {
705                    Some(current.pwr_solar_load[last_idx])
706                };
707                let last_charge_power = if last_idx >= current.pwr_max_chrg.len() {
708                    None
709                } else {
710                    Some(current.pwr_max_chrg[last_idx])
711                };
712                current.time = current.time.iter().map(|t| *t - current.time[0]).collect();
713                microtrips.push(current.clone());
714                current = Cycle {
715                    name: self.name.clone(),
716                    init_elev: last_elevation,
717                    time: vec![last_time],
718                    speed: vec![last_speed],
719                    dist: vec![],
720                    grade: if let Some(g) = last_grade {
721                        vec![g]
722                    } else {
723                        vec![]
724                    },
725                    elev: vec![],
726                    pwr_max_chrg: if let Some(p) = last_charge_power {
727                        vec![p]
728                    } else {
729                        vec![]
730                    },
731                    temp_amb_air: if let Some(temp) = last_temperature {
732                        vec![temp]
733                    } else {
734                        vec![]
735                    },
736                    pwr_solar_load: if let Some(p) = last_solar_load {
737                        vec![p]
738                    } else {
739                        vec![]
740                    },
741                    grade_interp: self.grade_interp.clone(),
742                    elev_interp: self.elev_interp.clone(),
743                };
744            }
745            current
746                .push(element.clone())
747                .expect("Push shouldn't have an error path");
748            moving = element.speed > stop_speed;
749        }
750        if current.time.len() > 1 {
751            current.time = current.time.iter().map(|t| *t - current.time[0]).collect();
752            current.init().unwrap();
753            microtrips.push(current.clone());
754        }
755        microtrips
756    }
757
758    /// Determine average speed of cycle.
759    /// -- while_moving: if true, only takes average while moving
760    ///
761    /// RETURN: average speed
762    pub fn average_speed(&self, while_moving: bool) -> si::Velocity {
763        let mut d = si::Length::ZERO;
764        let mut t = si::Time::ZERO;
765        for idx in 1..self.speed.len() {
766            let dt = self.time[idx] - self.time[idx - 1];
767            let vavg = 0.5 * (self.speed[idx] + self.speed[idx - 1]);
768            let dd = vavg * dt;
769            let no_move = (dd.get::<si::meter>().ceil() as i32) == 0;
770            d += dd;
771            t += if while_moving && no_move {
772                si::Time::ZERO
773            } else {
774                dt
775            };
776        }
777        if t > si::Time::ZERO {
778            d / t
779        } else {
780            si::Velocity::ZERO
781        }
782    }
783
784    /// Return the average step speeds of the cycle as vector of velicities.
785    /// NOTE: the average speed from sample i-1 to i will appear as entry i.
786    /// RETURN: vector of velocities representing average step speeds.
787    pub fn average_step_speeds(&self) -> Vec<si::Velocity> {
788        let mut result = Vec::with_capacity(self.time.len());
789        result.push(0.0 * uc::MPS);
790        for i in 1..self.time.len() {
791            result.push(0.5 * (self.speed[i] + self.speed[i - 1]));
792        }
793        result
794    }
795
796    /// Calculate the average step speed at step i
797    /// (i.e., from sample point i-1 to i)
798    pub fn average_step_speed_at(&self, i: usize) -> si::Velocity {
799        if i >= self.speed.len() {
800            return 0.0 * uc::MPS;
801        }
802        0.5 * (self.speed[i] + self.speed[i - 1])
803    }
804
805    /// The distances traveled over each step using trapezoidal
806    /// integration.
807    pub fn trapz_step_distances(&self) -> Vec<si::Length> {
808        let mut result = Vec::with_capacity(self.time.len());
809        result.push(0.0 * uc::M);
810        for i in 1..self.time.len() {
811            let step_time = self.time[i] - self.time[i - 1];
812            let average_speed = 0.5 * (self.speed[i] + self.speed[i - 1]);
813            result.push(step_time * average_speed);
814        }
815        result
816    }
817
818    /// The elevation climb each step using trapezoidal integration.
819    // TODO: verify the height calculation is correct, see cycle init changes
820    pub fn trapz_step_elevations(&self) -> Vec<si::Length> {
821        let mut result = Vec::with_capacity(self.time.len());
822        result.push(0.0 * uc::M);
823        for i in 1..self.time.len() {
824            let step_time = self.time[i].get::<si::second>() - self.time[i - 1].get::<si::second>();
825            let average_speed = 0.5
826                * (self.speed[i].get::<si::meter_per_second>()
827                    + self.speed[i - 1].get::<si::meter_per_second>());
828            let step_dist = step_time * average_speed;
829            let gr = self.grade[i].get::<si::ratio>();
830            let dh = gr.atan().cos() * step_dist * gr;
831            result.push(dh * uc::M);
832        }
833        result
834    }
835
836    /// The distance traveled from start to the beginning of step i
837    /// (i.e., distance traveled up to sample point i-1)
838    pub fn trapz_step_start_distance(&self, step: usize) -> si::Length {
839        let mut distance = 0.0 * uc::M;
840        let step_max = cmp::min(step, self.time.len());
841        for i in 1..step_max {
842            let step_time = self.time[i] - self.time[i - 1];
843            let average_speed = 0.5 * (self.speed[i] + self.speed[i - 1]);
844            distance += step_time * average_speed;
845        }
846        distance
847    }
848
849    /// The distance traveled during the given step
850    /// (i.e., distance from sample point i-1 to i for step i)
851    pub fn trapz_distance_for_step(&self, step: usize) -> si::Length {
852        let average_speed = self.average_step_speed_at(step);
853        let elapsed_time = self.time[step] - self.time[step - 1];
854        average_speed * elapsed_time
855    }
856
857    /// Calculate the distance from step i_start to the start of step i_end
858    /// (i.e., distance from sample point i_start - 1 to i_end - 1)
859    pub fn trapz_distance_over_range(&self, step0: usize, step1: usize) -> si::Length {
860        let distances = self.trapz_step_distances();
861        let last_i = cmp::max(distances.len() - 1, 0);
862        let i_start = cmp::min(step0, last_i);
863        let i_end = cmp::min(step1, last_i);
864        let mut distance = 0.0 * uc::M;
865        for d in &distances[cmp::min(i_start, i_end)..cmp::max(i_start, i_end)] {
866            distance += *d;
867        }
868        distance
869    }
870
871    /// Calculate the time in a cycle spent moving
872    /// - stopped_speed_m_per_s: the speed above which we are considered to be moving
873    ///
874    /// RETURN: the time spent moving in seconds
875    pub fn time_spent_moving(&self, stopped_speed: Option<si::Velocity>) -> si::Time {
876        let stop_speed = stopped_speed.unwrap_or(0.0 * uc::MPS);
877        let mut result = 0.0 * uc::S;
878        for i in 1..self.time.len() {
879            let step_time = self.time[i] - self.time[i - 1];
880            if self.speed[i] > stop_speed || self.speed[i - 1] > stop_speed {
881                result += step_time;
882            }
883        }
884        result
885    }
886
887    /// Create distance and target speeds by microtrip.
888    /// Splits cycle into microtrips and returns a list of
889    /// 2-tuples of:
890    /// (distance from start in meters, target speed in m/s)
891    /// The distance is measured to the start of the microtrip.
892    ///
893    /// # Parameters
894    ///
895    /// * `blend_factor`: from 0.0 to 1.0
896    ///    - if 0.0, use the average speed of the microtrip
897    ///    - if 1.0, use the average speed while moving (i.e., no stopped time)
898    ///    - otherwise, something in between
899    /// * `min_target_speed`: the minimum target speed allowed
900    ///
901    /// # Result
902    ///
903    /// List of 2-tuple of (distance from start, target speed).
904    /// A tuple represents the distance from start of the start
905    /// of the given microtrip and its target speed.
906    ///
907    /// # Notes
908    ///
909    /// * target speed per microtrip is not allowed to be
910    ///   below the `min_target_speed`
911    pub fn distance_and_target_speeds_by_microtrip(
912        &self,
913        stop_speed: Option<si::Velocity>,
914        blend_factor: f64,
915        min_target_speed: si::Velocity,
916    ) -> Vec<(si::Length, si::Velocity)> {
917        let blend_factor = blend_factor.clamp(0.0, 1.0);
918        let mut result = Vec::new();
919        let microtrips = self.to_microtrips(stop_speed);
920        let mut distance_at_start = 0.0 * uc::M;
921        let t0 = 0.0 * uc::S;
922        let v0 = 0.0 * uc::MPS;
923        let d0 = 0.0 * uc::M;
924        for mt in microtrips {
925            let distance = mt
926                .trapz_step_distances()
927                .iter()
928                .fold(0.0 * uc::M, |total, dist| total + *dist);
929            let last_index = cmp::max(mt.time.len() - 1, 0);
930            let end_time = mt.time[last_index];
931            let start_time = mt.time[0];
932            let total_time = end_time - start_time;
933            let moving_time = mt.time_spent_moving(stop_speed);
934            let average_speed = if total_time > t0 {
935                distance / total_time
936            } else {
937                v0
938            };
939            let moving_average_speed = if moving_time > t0 {
940                distance / moving_time
941            } else {
942                v0
943            };
944            let target_speed =
945                blend_factor * (moving_average_speed - average_speed) + average_speed;
946            let target_speed = if target_speed > min_target_speed {
947                target_speed
948            } else {
949                min_target_speed
950            };
951            if distance > d0 {
952                result.push((distance_at_start, target_speed));
953                distance_at_start += distance;
954            }
955        }
956        result
957    }
958
959    /// Add idle time to Cycle.
960    /// By "idle" time, we mean "stopped" time (i.e., vehicle not moving).
961    pub fn extend_time(
962        &self,
963        absolute_time: Option<si::Time>,
964        time_fraction: Option<si::Ratio>,
965    ) -> Cycle {
966        let absolute_time = absolute_time.unwrap_or(0.0 * uc::S);
967        let time_fraction = time_fraction.unwrap_or(0.0 * uc::R);
968        let mut ts = self.time.clone();
969        let mut vs = self.speed.clone();
970        let mut gs = self.grade.clone();
971        let mut ps = self.pwr_max_chrg.clone();
972        let mut temps = self.temp_amb_air.clone();
973        let mut ss = self.pwr_solar_load.clone();
974        let t_end = *ts.last().unwrap();
975        let extra_time_s = (absolute_time.get::<si::second>()
976            + time_fraction.get::<si::ratio>() * t_end.get::<si::second>())
977        .round() as i32;
978        if extra_time_s == 0 {
979            return self.clone();
980        }
981        let dt = 1.0 * uc::S;
982        let dt_s = dt.get::<si::second>();
983        let mut idx = 1;
984        loop {
985            let dt_extra_s = dt_s * idx as f64;
986            if dt_extra_s > extra_time_s as f64 {
987                break;
988            }
989            ts.push(t_end + dt_extra_s * uc::S);
990            vs.push(0.0 * uc::MPS);
991            if !gs.is_empty() {
992                gs.push(0.0 * uc::R);
993            }
994            if !ps.is_empty() {
995                ps.push(*ps.last().unwrap());
996            }
997            if !temps.is_empty() {
998                temps.push(*temps.last().unwrap());
999            }
1000            if !ss.is_empty() {
1001                ss.push(*ss.last().unwrap());
1002            }
1003            idx += 1;
1004        }
1005        let mut cyc = Cycle {
1006            name: self.name.clone(),
1007            init_elev: self.init_elev,
1008            time: ts,
1009            speed: vs,
1010            dist: vec![],
1011            grade: gs,
1012            elev: vec![],
1013            pwr_max_chrg: vec![],
1014            grade_interp: self.grade_interp.clone(),
1015            elev_interp: self.elev_interp.clone(),
1016            temp_amb_air: temps,
1017            pwr_solar_load: ss,
1018        };
1019        cyc.init().unwrap();
1020        cyc
1021    }
1022
1023    /// Create a cache object for faster computations on Cycle.
1024    pub fn build_cache(&self) -> CycleCache {
1025        CycleCache::new(self)
1026    }
1027
1028    /// Returns the average grade over the given range of distances.
1029    /// - distance_start: the distance at start of evaluation area
1030    /// - delta_distance: distance traveled from distance_start
1031    /// - cache: optional CycleCache which can save computation time
1032    ///
1033    /// RETURN: average grade (rise over run) for the given range.
1034    ///
1035    /// NOTE: grade is assumed to be constant from just after the
1036    /// previous sample point until the current sample point (inclusive).
1037    /// That is, grade[i] applies from distance, d, of (d[i - 1], d[i]]
1038    pub fn average_grade_over_range(
1039        &self,
1040        distance_start: si::Length,
1041        delta_distance: si::Length,
1042        cache: Option<&CycleCache>,
1043    ) -> si::Ratio {
1044        let tol = 1e-6;
1045        match &cache {
1046            Some(cc) => {
1047                let dd_m = delta_distance.get::<si::meter>();
1048                if cc.grade_all_zero {
1049                    0.0 * uc::R
1050                } else if dd_m <= tol {
1051                    let dist_m = distance_start.get::<si::meter>();
1052                    cc.interp_grade(dist_m) * uc::R
1053                } else {
1054                    let dist0_m = distance_start.get::<si::meter>();
1055                    let dist1_m = dist0_m + dd_m;
1056                    let e0 = cc.interp_elevation(dist0_m);
1057                    let e1 = cc.interp_elevation(dist1_m);
1058                    ((e1 - e0) / dd_m).asin().tan() * uc::R
1059                }
1060            }
1061            None => {
1062                let zero_grade = 0.0 * uc::R;
1063                let grade_all_zero = {
1064                    let mut all0 = true;
1065                    for idx in 0..self.grade.len() {
1066                        if self.grade[idx] != zero_grade {
1067                            all0 = false;
1068                            break;
1069                        }
1070                    }
1071                    all0
1072                };
1073                if grade_all_zero {
1074                    0.0 * uc::R
1075                } else {
1076                    let delta_dists_m: Vec<f64> = self
1077                        .trapz_step_distances()
1078                        .iter()
1079                        .map(|dd| dd.get::<si::meter>())
1080                        .collect();
1081                    let trapz_distances_m = {
1082                        let mut d = 0.0;
1083                        let mut result = Vec::with_capacity(delta_dists_m.len());
1084                        for dd in &delta_dists_m {
1085                            d += *dd;
1086                            result.push(d);
1087                        }
1088                        result
1089                    };
1090                    let dist0_m = distance_start.get::<si::meter>();
1091                    let dd_m = delta_distance.get::<si::meter>();
1092                    let dist1_m = dist0_m + dd_m;
1093                    if dd_m < tol {
1094                        if dist0_m < trapz_distances_m[0] {
1095                            return self.grade[0];
1096                        }
1097                        let max_idx = self.grade.len() - 1;
1098                        if dist0_m > trapz_distances_m[max_idx] {
1099                            return self.grade[max_idx];
1100                        }
1101                        for idx in 1..self.time.len() {
1102                            if dist0_m > trapz_distances_m[idx - 1]
1103                                && dist0_m <= trapz_distances_m[idx]
1104                            {
1105                                return self.grade[idx];
1106                            }
1107                        }
1108                        self.grade[max_idx]
1109                    } else {
1110                        // NOTE: we use the following instead of delta_elev_m
1111                        // as it uses more precise trapezoidal diatance and
1112                        // elevation at sample points. This also uses the
1113                        // fully accurate trig functions in case we have large
1114                        // slope angles. This level of rigor may be overkill.
1115                        let trapz_elevations_m = {
1116                            let delta_elevs_m: Vec<f64> = self
1117                                .grade
1118                                .iter()
1119                                .zip(delta_dists_m)
1120                                .map(|(g, dd)| {
1121                                    let gr = g.get::<si::ratio>();
1122                                    gr.atan().cos() * dd * gr
1123                                })
1124                                .collect();
1125                            let mut result = Vec::with_capacity(delta_elevs_m.len());
1126                            let mut elev_m = 0.0;
1127                            for de in &delta_elevs_m {
1128                                elev_m += *de;
1129                                result.push(elev_m);
1130                            }
1131                            result
1132                        };
1133                        let interp: InterpolatorEnum<ndarray::OwnedRepr<f64>> =
1134                            InterpolatorEnum::new_1d(
1135                                trapz_distances_m.clone().into(),
1136                                trapz_elevations_m.clone().into(),
1137                                strategy::Linear,
1138                                Extrapolate::Clamp,
1139                            )
1140                            .unwrap();
1141                        let e0_m = interp.interpolate(&[dist0_m]).unwrap();
1142                        let e1_m = interp.interpolate(&[dist1_m]).unwrap();
1143                        ((e1_m - e0_m) / dd_m).asin().tan() * uc::R
1144                    }
1145                }
1146            }
1147        }
1148    }
1149
1150    /// Calculate the distance to next stop from `distance`.
1151    /// - distance: the distance to calculate distance-to-stop from
1152    ///
1153    /// RETURN: returns the distance to the next stop from `distance`
1154    ///
1155    /// NOTE: distance may be negative if we're beyond the last stop
1156    pub fn calc_distance_to_next_stop_from(
1157        &self,
1158        distance: si::Length,
1159        cache: Option<&CycleCache>,
1160    ) -> si::Length {
1161        let tol = 1e-6;
1162        let distance_m = distance.get::<si::meter>();
1163        match cache {
1164            Some(cc) => {
1165                for (&d_m, &v) in cc.trapz_distances_m.iter().zip(self.speed.iter()) {
1166                    let v_mps = v.get::<si::meter_per_second>();
1167                    if (v_mps < tol) && (d_m > (distance_m + tol)) {
1168                        return (d_m - distance_m) * uc::M;
1169                    }
1170                }
1171                (*cc.trapz_distances_m.last().unwrap_or(&0.0) * uc::M) - distance
1172            }
1173            None => {
1174                let ds_m = {
1175                    let mut result = Vec::with_capacity(self.time.len());
1176                    let mut d_m = 0.0;
1177                    for dd in self.trapz_step_distances() {
1178                        let dd_m = dd.get::<si::meter>();
1179                        d_m += dd_m;
1180                        result.push(d_m);
1181                    }
1182                    result
1183                };
1184                for (&d_m, &v) in ds_m.iter().zip(self.speed.iter()) {
1185                    let v_mps = v.get::<si::meter_per_second>();
1186                    if (v_mps < tol) && (d_m > (distance_m + tol)) {
1187                        return (d_m - distance_m) * uc::M;
1188                    }
1189                }
1190                *ds_m.last().unwrap_or(&0.0) * uc::M
1191            }
1192        }
1193    }
1194
1195    /// Modify the cycle using the given constant-jerk trajectory.
1196    /// - i: the index into the cycle to initiate modification
1197    ///   NOTE: THIS point is modified as trajectory is calculated as
1198    ///   starting at i-1
1199    /// - n: the number of steps ahead
1200    /// - jerk: the jerk (deriviative of acceleration with time)
1201    /// - accel0: the starting accelartion
1202    ///
1203    /// NOTE:
1204    /// - modifies the cycle in-place. Purpose is to allow hitting
1205    ///   a rendezvous point in time/speed in the future.
1206    /// - CAUTION: not robust against variable duration time-steps
1207    ///
1208    /// RETURN: the final modified speed
1209    pub fn modify_by_const_jerk_trajectory(
1210        &mut self,
1211        i: usize,
1212        n: usize,
1213        jerk: si::Jerk,
1214        accel0: si::Acceleration,
1215    ) -> si::Velocity {
1216        if n == 0 {
1217            return si::Velocity::ZERO;
1218        }
1219        let jerk_m_per_s3 = jerk.get::<si::meter_per_second_cubed>();
1220        let accel0_m_per_s2 = accel0.get::<si::meter_per_second_squared>();
1221        let num_samples = self.speed.len();
1222        if i >= num_samples {
1223            if num_samples > 0 {
1224                return self.speed[num_samples - 1];
1225            }
1226            return si::Velocity::ZERO;
1227        }
1228        let v0 = self.speed[i - 1].get::<si::meter_per_second>();
1229        let dt = (self.time[i] - self.time[i - 1]).get::<si::second>();
1230        let mut v = v0;
1231        for ni in 1..(n + 1) {
1232            let idx_to_set = (i - 1) + ni;
1233            if idx_to_set >= num_samples {
1234                break;
1235            }
1236            v = speed_for_constant_jerk(ni, v0, accel0_m_per_s2, jerk_m_per_s3, dt);
1237            self.speed[idx_to_set] = v.max(0.0) * uc::MPS;
1238        }
1239        self.init().unwrap();
1240        v * uc::MPS
1241    }
1242
1243    /// Modify cycle to add a braking trajectory that would cover the same
1244    /// distance as the given constant brake deceleration.
1245    /// - brake_accel: the brake acceleration (m/s2); must be negative
1246    /// - i: index where to initiate the stop trectory; start of the step
1247    /// - desired_distance_to_stop: the desired distance to stop within. If
1248    ///   not provided, it is calculated based on the braking deceleration.
1249    ///
1250    /// RETURN: (final speed of modified trajectory, number of steps to complete)
1251    /// - the final speed should be zero ideally
1252    /// - the number of time-steps required to complete the braking maneuver
1253    ///
1254    /// NOTE:
1255    /// - modifies the cycle in-place.
1256    pub fn modify_with_braking_trajectory(
1257        &mut self,
1258        brake_accel: si::Acceleration,
1259        i: usize,
1260        desired_distance_to_stop: Option<si::Length>,
1261    ) -> (si::Velocity, usize) {
1262        let brake_accel = if brake_accel > si::Acceleration::ZERO {
1263            -brake_accel
1264        } else {
1265            brake_accel
1266        };
1267        assert!(brake_accel < si::Acceleration::ZERO);
1268        if i >= self.time.len() {
1269            return (*self.speed.last().unwrap(), 0);
1270        }
1271        let i = if i < 1 { 1 } else { i };
1272        let v0 = self.speed[i - 1].get::<si::meter_per_second>();
1273        let dt = (self.time[i] - self.time[i - 1]).get::<si::second>();
1274        let brake_accel_m_per_s2 = brake_accel.get::<si::meter_per_second_squared>();
1275        // distance-to-stop (m)
1276        let dts_m = match desired_distance_to_stop {
1277            Some(value) => {
1278                let result = value.get::<si::meter>();
1279                if result > 0.0 {
1280                    result
1281                } else {
1282                    -0.5 * v0 * v0 / brake_accel_m_per_s2
1283                }
1284            }
1285            None => -0.5 * v0 * v0 / brake_accel_m_per_s2,
1286        };
1287        if dts_m <= 0.0 {
1288            return (v0 * uc::MPS, 0);
1289        }
1290        // time-to-stop (s)
1291        let tts_s = -v0 / brake_accel_m_per_s2;
1292        // number of steps to stop
1293        let n = (tts_s / dt).round() as usize;
1294        let n = if n < 2 { 2 } else { n }; // need at least 2 steps
1295        let traj =
1296            ConstantJerkTrajectory::from_speed_and_distance_targets(n, 0.0, v0, dts_m, 0.0, dt);
1297        let v_final = self.modify_by_const_jerk_trajectory(
1298            i,
1299            n,
1300            traj.jerk_m_per_s3 * uc::MPS3,
1301            traj.acceleration_m_per_s2 * uc::MPS2,
1302        );
1303        (v_final, n)
1304    }
1305
1306    /// Report the stopped time (i.e., idle) at the end of a cycle.
1307    ///
1308    /// RESULT: time vehicle is at zero speed at cycle end
1309    pub fn ending_idle_time(&self) -> si::Time {
1310        let mut result = si::Time::ZERO;
1311        let vzero = si::Velocity::ZERO;
1312        for idx in (1..self.time.len()).rev() {
1313            let v0 = self.speed[idx - 1];
1314            let v1 = self.speed[idx];
1315            if v0 != vzero || v1 != vzero {
1316                break;
1317            } else {
1318                let dt = self.time[idx] - self.time[idx - 1];
1319                result += dt;
1320            }
1321        }
1322        result
1323    }
1324
1325    /// Remove idel time at end of cycle except for the optionally
1326    /// specified duration.
1327    /// - idle_to_keep: optional duration of idle time to keep. Default is 0 s
1328    ///
1329    /// RESULT: a new cycle with idle time trimmed.
1330    pub fn trim_ending_idle(&self, idle_to_keep: Option<si::Time>) -> Cycle {
1331        let idle_to_keep = idle_to_keep.unwrap_or(si::Time::ZERO).max(si::Time::ZERO);
1332        let vzero = si::Velocity::ZERO;
1333        let mut idle_start_idx = 0;
1334        for idx in (1..self.time.len()).rev() {
1335            let v0 = self.speed[idx - 1];
1336            let v1 = self.speed[idx];
1337            if v0 != vzero || v1 != vzero {
1338                idle_start_idx = idx + 1;
1339                break;
1340            }
1341        }
1342        if idle_start_idx >= self.time.len() {
1343            return self.clone();
1344        }
1345        let end_idx = if idle_to_keep == si::Time::ZERO {
1346            idle_start_idx
1347        } else {
1348            let mut dt_idle = si::Time::ZERO;
1349            let mut idx_drop = idle_start_idx;
1350            for idx in idle_start_idx..self.time.len() {
1351                let dt = self.time[idx] - self.time[idx - 1];
1352                dt_idle += dt;
1353                if dt_idle > idle_to_keep {
1354                    idx_drop = idx;
1355                    break;
1356                }
1357            }
1358            idx_drop
1359        };
1360        let mut cyc = Cycle {
1361            name: self.name.clone(),
1362            time: self.time[0..end_idx].to_vec(),
1363            speed: self.speed[0..end_idx].to_vec(),
1364            init_elev: self.init_elev,
1365            grade: if self.grade.is_empty() {
1366                vec![]
1367            } else {
1368                self.grade[0..end_idx].to_vec()
1369            },
1370            dist: vec![],
1371            elev: vec![],
1372            pwr_max_chrg: if self.pwr_max_chrg.is_empty() {
1373                vec![]
1374            } else {
1375                self.pwr_max_chrg[0..end_idx].to_vec()
1376            },
1377            temp_amb_air: if self.temp_amb_air.is_empty() {
1378                vec![]
1379            } else {
1380                self.temp_amb_air[0..end_idx].to_vec()
1381            },
1382            pwr_solar_load: if self.pwr_solar_load.is_empty() {
1383                vec![]
1384            } else {
1385                self.pwr_solar_load[0..end_idx].to_vec()
1386            },
1387            grade_interp: None,
1388            elev_interp: None,
1389        };
1390        cyc.init().unwrap();
1391        cyc
1392    }
1393
1394    /// Resample cycle to a lower or higher frequency.
1395    /// - dt: the new step duration.
1396    ///
1397    /// RETURN: cycle
1398    /// NOTE: a value of dt <= 0 s implies to just clone the current cycle
1399    /// "as is"
1400    pub fn resample(&self, dt: si::Time) -> Cycle {
1401        if dt <= si::Time::ZERO {
1402            return self.clone();
1403        }
1404        let mut t = si::Time::ZERO;
1405        let speed_interp: InterpolatorEnum<OwnedRepr<f64>> = InterpolatorEnum::new_1d(
1406            self.time.iter().map(|x| x.get::<si::second>()).collect(),
1407            self.speed
1408                .iter()
1409                .map(|y| y.get::<si::meter_per_second>())
1410                .collect(),
1411            strategy::Linear,
1412            Extrapolate::Clamp,
1413        )
1414        .unwrap();
1415        let grade_interp: InterpolatorEnum<OwnedRepr<f64>> = InterpolatorEnum::new_1d(
1416            self.time.iter().map(|x| x.get::<si::second>()).collect(),
1417            self.grade.iter().map(|y| y.get::<si::ratio>()).collect(),
1418            strategy::RightNearest,
1419            Extrapolate::Clamp,
1420        )
1421        .unwrap();
1422        let temp_interp: Option<InterpolatorEnum<OwnedRepr<f64>>> =
1423            if self.temp_amb_air.len() == self.time.len() {
1424                Some(
1425                    InterpolatorEnum::new_1d(
1426                        self.time.iter().map(|t| t.get::<si::second>()).collect(),
1427                        self.temp_amb_air
1428                            .iter()
1429                            .map(|temp| temp.get::<si::kelvin_abs>())
1430                            .collect(),
1431                        strategy::Linear,
1432                        Extrapolate::Clamp,
1433                    )
1434                    .unwrap(),
1435                )
1436            } else {
1437                None
1438            };
1439        let solar_interp: Option<InterpolatorEnum<OwnedRepr<f64>>> =
1440            if self.pwr_solar_load.len() == self.time.len() {
1441                Some(
1442                    InterpolatorEnum::new_1d(
1443                        self.time.iter().map(|t| t.get::<si::second>()).collect(),
1444                        self.pwr_solar_load
1445                            .iter()
1446                            .map(|p| p.get::<si::kilowatt>())
1447                            .collect(),
1448                        strategy::Linear,
1449                        Extrapolate::Clamp,
1450                    )
1451                    .unwrap(),
1452                )
1453            } else {
1454                None
1455            };
1456        let chg_pwr_interp: Option<InterpolatorEnum<OwnedRepr<f64>>> =
1457            if self.pwr_max_chrg.len() == self.time.len() {
1458                Some(
1459                    InterpolatorEnum::new_1d(
1460                        self.time.iter().map(|t| t.get::<si::second>()).collect(),
1461                        self.pwr_max_chrg
1462                            .iter()
1463                            .map(|p| p.get::<si::kilowatt>())
1464                            .collect(),
1465                        strategy::Linear,
1466                        Extrapolate::Clamp,
1467                    )
1468                    .unwrap(),
1469                )
1470            } else {
1471                None
1472            };
1473        let mut ts = vec![];
1474        let mut vs = vec![];
1475        let mut gs = vec![];
1476        let mut pwr_chg = vec![];
1477        let mut temps = vec![];
1478        let mut solars = vec![];
1479        while t <= self.time[self.time.len() - 1] {
1480            ts.push(t);
1481            let t0 = t.get::<si::second>();
1482            let v = speed_interp.interpolate(&[t0]).unwrap();
1483            vs.push(v * uc::MPS);
1484            let g = grade_interp.interpolate(&[t0]).unwrap();
1485            gs.push(g * uc::R);
1486            if let Some(ref interp) = chg_pwr_interp {
1487                let pchg = interp.interpolate(&[t0]).unwrap();
1488                pwr_chg.push(pchg * uc::KW);
1489            }
1490            if let Some(ref interp) = temp_interp {
1491                let temp = interp.interpolate(&[t0]).unwrap();
1492                temps.push(temp * uc::KELVIN);
1493            }
1494            if let Some(ref interp) = solar_interp {
1495                let solar = interp.interpolate(&[t0]).unwrap();
1496                solars.push(solar * uc::KW);
1497            }
1498            t += dt;
1499        }
1500
1501        let mut cyc = Cycle {
1502            name: self.name.clone(),
1503            init_elev: self.init_elev,
1504            time: ts,
1505            speed: vs,
1506            dist: vec![],
1507            grade: gs,
1508            elev: vec![],
1509            pwr_max_chrg: pwr_chg,
1510            temp_amb_air: temps,
1511            pwr_solar_load: solars,
1512            grade_interp: None,
1513            elev_interp: None,
1514        };
1515        cyc.init().unwrap();
1516        cyc
1517    }
1518}
1519
1520impl TryFrom<CycleBuilder> for Cycle {
1521    type Error = anyhow::Error;
1522    fn try_from(value: CycleBuilder) -> anyhow::Result<Self, Self::Error> {
1523        let mut cyc = Self {
1524            name: value.name,
1525            init_elev: None,
1526            time: value.time,
1527            speed: value.speed,
1528            dist: Default::default(),
1529            grade: Default::default(),
1530            elev: Default::default(),
1531            pwr_max_chrg: Default::default(),
1532            temp_amb_air: Default::default(),
1533            pwr_solar_load: Default::default(),
1534            grade_interp: None,
1535            elev_interp: Default::default(),
1536        };
1537        cyc.init()?;
1538        Ok(cyc)
1539    }
1540}
1541
1542/// Trait for CycleBuilder and Cycle to support builder pattern
1543pub trait CBTrait {
1544    /// Return cycle with `grade`
1545    fn with_grade(&mut self, grade: Vec<si::Ratio>) -> anyhow::Result<Cycle>;
1546
1547    /// Return cycle with `temp_amb_air`
1548    fn with_temp_amb_air(&mut self, temp_amb_air: Vec<si::Temperature>) -> anyhow::Result<Cycle>;
1549
1550    // TODO: add more of these builder helpers
1551}
1552
1553impl CBTrait for Cycle {
1554    fn with_grade(&mut self, grade: Vec<si::Ratio>) -> anyhow::Result<Cycle> {
1555        ensure!(
1556            self.len_checked().with_context(|| format_dbg!())? == grade.len(),
1557            format!(
1558                "{}\n`self.len()`: `{}\n`grade.len()`",
1559                self.len_checked().with_context(|| format_dbg!())?,
1560                grade.len()
1561            )
1562        );
1563        self.grade = grade;
1564        Ok(self.clone())
1565    }
1566
1567    fn with_temp_amb_air(&mut self, temp_amb_air: Vec<si::Temperature>) -> anyhow::Result<Cycle> {
1568        ensure!(
1569            self.len_checked().with_context(|| format_dbg!())? == temp_amb_air.len(),
1570            format!(
1571                "{}\n`self.len()`: `{}\n`temp_amb_air.len()`",
1572                self.len_checked().with_context(|| format_dbg!())?,
1573                temp_amb_air.len()
1574            )
1575        );
1576        self.temp_amb_air = temp_amb_air;
1577        Ok(self.clone())
1578    }
1579}
1580
1581#[serde_api]
1582#[derive(Default, Debug, Serialize, Deserialize, PartialEq, Clone)]
1583#[non_exhaustive]
1584/// Simple cycle to be converted into [Cycle] with appropriate defaults
1585pub struct CycleBuilder {
1586    /// Name of cycle (can be left empty)
1587    #[serde(default, skip_serializing_if = "String::is_empty")]
1588    pub name: String,
1589    /// simulation time
1590    pub time: Vec<si::Time>,
1591    /// prescribed speed
1592    pub speed: Vec<si::Velocity>,
1593}
1594
1595impl CBTrait for CycleBuilder {
1596    fn with_grade(&mut self, grade: Vec<si::Ratio>) -> anyhow::Result<Cycle> {
1597        let mut cyc: Cycle = self.clone().try_into().with_context(|| format_dbg!())?;
1598        cyc.grade = grade;
1599        Ok(cyc)
1600    }
1601
1602    fn with_temp_amb_air(&mut self, temp_amb_air: Vec<si::Temperature>) -> anyhow::Result<Cycle> {
1603        let mut cyc: Cycle = self.clone().try_into().with_context(|| format_dbg!())?;
1604        cyc.temp_amb_air = temp_amb_air;
1605        Ok(cyc)
1606    }
1607}
1608
1609#[serde_api]
1610#[derive(Default, Debug, Serialize, Deserialize, PartialEq, Clone)]
1611#[non_exhaustive]
1612#[serde(deny_unknown_fields)]
1613#[cfg_attr(feature = "pyo3", pyclass(module = "fastsim", subclass, eq))]
1614/// Element of `Cycle`.  Used for vec-like operations.
1615pub struct CycleElement {
1616    /// simulation time \[s\]
1617    #[serde(alias = "cycSecs")]
1618    pub time: si::Time,
1619    /// simulation power \[W\]
1620    #[serde(alias = "speed_mps", alias = "cycMps")]
1621    pub speed: si::Velocity,
1622    // `dist` is not included here because it is derived in `Init::init`
1623    /// road grade
1624    #[serde(alias = "cycGrade")]
1625    pub grade: Option<si::Ratio>,
1626    // `elev` is not included here because it is derived in `Init::init`
1627    /// road charging/discharing capacity
1628    pub pwr_max_charge: Option<si::Power>,
1629    // TODO: make sure all fields in cycle are represented here, as appropriate
1630    /// ambient air temperature w.r.t. to time (rather than spatial position)
1631    pub temp_amb_air: Option<si::Temperature>,
1632    /// solar heat load w.r.t. to time (rather than spatial position)
1633    pub pwr_solar_load: Option<si::Power>,
1634}
1635
1636impl SerdeAPI for CycleElement {}
1637impl Init for CycleElement {}
1638
1639#[pyo3_api]
1640impl CycleElement {}
1641
1642#[cfg(test)]
1643mod tests {
1644    use super::{manipulation_utils::ConstantJerkTrajectory, *};
1645    /// Build, initialize, and return 2-element cycle
1646    fn mock_cyc_len_2() -> Cycle {
1647        let mut cyc = Cycle {
1648            name: String::new(),
1649            init_elev: None,
1650            time: (0..=2).map(|x| (x as f64) * uc::S).collect(),
1651            speed: (0..=2).map(|x| (x as f64) * uc::MPS).collect(),
1652            dist: vec![],
1653            grade: (0..=2).map(|x| (x as f64 * uc::R) / 100.).collect(),
1654            elev: vec![],
1655            pwr_max_chrg: vec![],
1656            grade_interp: Default::default(),
1657            elev_interp: Default::default(),
1658            temp_amb_air: Default::default(),
1659            pwr_solar_load: Default::default(),
1660        };
1661        cyc.init().unwrap();
1662        cyc
1663    }
1664
1665    fn make_two_triangles_cycle() -> Cycle {
1666        let mut cyc = Cycle {
1667            name: String::from("Two Triangles"),
1668            init_elev: Some(0.0 * uc::M),
1669            time: vec![
1670                0.0 * uc::S,
1671                10.0 * uc::S,
1672                20.0 * uc::S,
1673                30.0 * uc::S,
1674                40.0 * uc::S,
1675                50.0 * uc::S,
1676            ],
1677            speed: vec![
1678                0.0 * uc::MPS,
1679                4.0 * uc::MPS,
1680                0.0 * uc::MPS,
1681                0.0 * uc::MPS,
1682                5.0 * uc::MPS,
1683                0.0 * uc::MPS,
1684            ],
1685            dist: vec![],
1686            grade: vec![
1687                0.0 * uc::R,
1688                0.0 * uc::R,
1689                0.0 * uc::R,
1690                0.0 * uc::R,
1691                0.01 * uc::R,
1692                0.01 * uc::R,
1693            ],
1694            elev: vec![],
1695            pwr_max_chrg: vec![],
1696            grade_interp: Default::default(),
1697            elev_interp: Default::default(),
1698            temp_amb_air: Default::default(),
1699            pwr_solar_load: Default::default(),
1700        };
1701        cyc.init().unwrap();
1702        cyc
1703    }
1704
1705    #[test]
1706    fn test_init() {
1707        let cyc = mock_cyc_len_2();
1708        assert_eq!(
1709            cyc.dist,
1710            [0., 1., 3.] // meters
1711                .iter()
1712                .map(|x| *x * uc::M)
1713                .collect::<Vec<si::Length>>()
1714        );
1715        assert_eq!(
1716            cyc.elev,
1717            [121.92, 121.9299995000375, 121.9699915024367] // meters
1718                .iter()
1719                .map(|x| *x * uc::M)
1720                .collect::<Vec<si::Length>>()
1721        );
1722    }
1723
1724    #[test]
1725    fn test_to_elements() {
1726        let cyc = mock_cyc_len_2();
1727        let elements = cyc.to_elements();
1728        assert_eq!(elements.len(), 3);
1729        assert_eq!(elements[0].time, 0.0 * uc::S);
1730        assert_eq!(elements[2].time, cyc.time[2]);
1731        assert_eq!(elements[2].speed, cyc.speed[2]);
1732        assert_eq!(elements[2].grade.unwrap(), 0.02 * uc::R);
1733        assert!(elements[2].pwr_max_charge.is_none());
1734        assert_eq!(elements[2].temp_amb_air.unwrap(), *TE_STD_AIR);
1735        assert!(elements[2].pwr_solar_load.is_none());
1736    }
1737
1738    #[test]
1739    fn test_to_microtrips() {
1740        let cyc = make_two_triangles_cycle();
1741        let actual = cyc.to_microtrips(Some(0.01 * uc::MPH));
1742        assert_eq!(actual.len(), 2);
1743        let cyc0 = &actual[0];
1744        assert_eq!(
1745            cyc0.time,
1746            vec![0.0 * uc::S, 10.0 * uc::S, 20.0 * uc::S, 30.0 * uc::S]
1747        );
1748        assert_eq!(
1749            cyc0.speed,
1750            vec![0.0 * uc::MPS, 4.0 * uc::MPS, 0.0 * uc::MPS, 0.0 * uc::MPS]
1751        );
1752        assert_eq!(
1753            cyc0.grade,
1754            vec![0.0 * uc::R, 0.0 * uc::R, 0.0 * uc::R, 0.0 * uc::R]
1755        );
1756        let cyc1 = &actual[1];
1757        assert_eq!(cyc1.time, vec![0.0 * uc::S, 10.0 * uc::S, 20.0 * uc::S]);
1758        assert_eq!(
1759            cyc1.speed,
1760            vec![0.0 * uc::MPS, 5.0 * uc::MPS, 0.0 * uc::MPS]
1761        );
1762        assert_eq!(cyc1.grade, vec![0.0 * uc::R, 0.01 * uc::R, 0.01 * uc::R]);
1763    }
1764
1765    #[test]
1766    fn test_distance_and_target_speeds_by_microtrip() {
1767        let cyc = make_two_triangles_cycle();
1768        let expected = [
1769            (0.0 * uc::M, (40.0 / 20.0) * uc::MPS),
1770            (40.0 * uc::M, (50.0 / 20.0) * uc::MPS),
1771        ];
1772        let actual = cyc.distance_and_target_speeds_by_microtrip(None, 1.0, 0.0 * uc::MPS);
1773        assert_eq!(actual.len(), expected.len());
1774        for i in 0..expected.len() {
1775            assert_eq!(actual[i].0, expected[i].0);
1776            assert_eq!(actual[i].1, expected[i].1);
1777        }
1778        let expected = [
1779            (0.0 * uc::M, (40.0 / 30.0) * uc::MPS),
1780            (40.0 * uc::M, (50.0 / 20.0) * uc::MPS),
1781        ];
1782        let actual = cyc.distance_and_target_speeds_by_microtrip(None, 0.0, 0.0 * uc::MPS);
1783        assert_eq!(actual.len(), expected.len());
1784        for i in 0..expected.len() {
1785            assert_eq!(actual[i].0, expected[i].0);
1786            assert_eq!(actual[i].1, expected[i].1);
1787        }
1788    }
1789
1790    #[test]
1791    fn test_extending_cycle_time() {
1792        let cyc = make_two_triangles_cycle();
1793        let expected = {
1794            let mut c = Cycle {
1795                name: String::from("Two Triangles"),
1796                init_elev: Some(0.0 * uc::M),
1797                time: vec![
1798                    0.0 * uc::S,
1799                    10.0 * uc::S,
1800                    20.0 * uc::S,
1801                    30.0 * uc::S,
1802                    40.0 * uc::S,
1803                    50.0 * uc::S,
1804                    51.0 * uc::S,
1805                    52.0 * uc::S,
1806                    53.0 * uc::S,
1807                    54.0 * uc::S,
1808                    55.0 * uc::S,
1809                    56.0 * uc::S,
1810                    57.0 * uc::S,
1811                    58.0 * uc::S,
1812                ],
1813                speed: vec![
1814                    0.0 * uc::MPS,
1815                    4.0 * uc::MPS,
1816                    0.0 * uc::MPS,
1817                    0.0 * uc::MPS,
1818                    5.0 * uc::MPS,
1819                    0.0 * uc::MPS,
1820                    0.0 * uc::MPS,
1821                    0.0 * uc::MPS,
1822                    0.0 * uc::MPS,
1823                    0.0 * uc::MPS,
1824                    0.0 * uc::MPS,
1825                    0.0 * uc::MPS,
1826                    0.0 * uc::MPS,
1827                    0.0 * uc::MPS,
1828                ],
1829                dist: vec![],
1830                grade: vec![
1831                    0.0 * uc::R,
1832                    0.0 * uc::R,
1833                    0.0 * uc::R,
1834                    0.0 * uc::R,
1835                    0.01 * uc::R,
1836                    0.01 * uc::R,
1837                    0.0 * uc::R,
1838                    0.0 * uc::R,
1839                    0.0 * uc::R,
1840                    0.0 * uc::R,
1841                    0.0 * uc::R,
1842                    0.0 * uc::R,
1843                    0.0 * uc::R,
1844                    0.0 * uc::R,
1845                ],
1846                elev: vec![],
1847                pwr_max_chrg: vec![],
1848                grade_interp: Default::default(),
1849                elev_interp: Default::default(),
1850                temp_amb_air: Default::default(),
1851                pwr_solar_load: Default::default(),
1852            };
1853            c.init().unwrap();
1854            c
1855        };
1856        let absolute_time = Some(3.0 * uc::S);
1857        let time_fraction = Some(0.10 * uc::R);
1858        // extend by 3 s and 10% of existing time (i.e., 5 s)
1859        // = extend by 8 s
1860        let actual = cyc.extend_time(absolute_time, time_fraction);
1861        assert_eq!(actual, expected);
1862    }
1863
1864    /// Round the given number n to the given number of digits
1865    /// - n: the number to round
1866    /// - digits: the digits to round or defaults to 2; if not positive,
1867    fn round(n: f64, digits: Option<i32>) -> f64 {
1868        let digits = digits.unwrap_or(2);
1869        let digits = if digits < 0 { 0 } else { digits };
1870        let multiplier = 10.0_f64.powi(digits);
1871        (n * multiplier).round() / multiplier
1872    }
1873
1874    #[test]
1875    fn cycle_step_distances_are_as_expected() {
1876        let c = make_two_triangles_cycle();
1877        let expected = [
1878            0.0 * uc::M,
1879            20.0 * uc::M,
1880            20.0 * uc::M,
1881            0.0 * uc::M,
1882            25.0 * uc::M,
1883            25.0 * uc::M,
1884        ];
1885        let actual = c.trapz_step_distances();
1886        assert_eq!(actual.len(), expected.len());
1887        for i in 0..expected.len() {
1888            assert_eq!(actual[i], expected[i], "differ at step {i}");
1889        }
1890    }
1891
1892    #[test]
1893    fn cycle_elevations_are_as_expected() {
1894        let c = make_two_triangles_cycle();
1895        let dh = 0.01_f64.atan().cos() * 25.0_f64 * 0.01_f64;
1896        let expected = [
1897            0.0 * uc::M,
1898            0.0 * uc::M,
1899            0.0 * uc::M,
1900            0.0 * uc::M,
1901            dh * uc::M,
1902            dh * uc::M,
1903        ];
1904        let actual = c.trapz_step_elevations();
1905        assert_eq!(actual.len(), expected.len());
1906        for i in 0..expected.len() {
1907            assert_eq!(actual[i], expected[i], "differ at step {i}");
1908        }
1909    }
1910
1911    #[test]
1912    fn test_elevation_accumulation() {
1913        let mut cyc = Cycle {
1914            name: String::from("elevation test"),
1915            init_elev: Some(0.0 * uc::M),
1916            time: Vec::linspace(0., 1000., 1001)
1917                .iter()
1918                .map(|x| (*x as f64) * uc::S)
1919                .collect(),
1920            speed: vec![20.0 * uc::MPS; 1001],
1921            dist: vec![],
1922            grade: vec![0.05 * uc::R; 1001],
1923            elev: vec![],
1924            pwr_max_chrg: vec![],
1925            grade_interp: Default::default(),
1926            elev_interp: Default::default(),
1927            temp_amb_air: Default::default(),
1928            pwr_solar_load: Default::default(),
1929        };
1930        cyc.init().unwrap();
1931
1932        let delta_elev = cyc.elev.last().unwrap().get::<si::meter>();
1933        // Expected elevation change: 20 m/s * 1000 s * sin(atan(0.05)) = 998.7523388778305 m
1934        assert!(almost_eq(delta_elev, 998.7523388778305, None));
1935    }
1936
1937    #[test]
1938    fn cycle_cache_yields_same_results() {
1939        let c = make_two_triangles_cycle();
1940        let cache = c.build_cache();
1941        let dist_m = 0.0;
1942        let e0_expected = 0.0;
1943        let e0_actual = cache.interp_elevation(dist_m);
1944        assert_eq!(e0_actual, e0_expected);
1945        let dist_m = 65.0;
1946        let e1_expected = 0.01_f64.atan().cos() * 25.0_f64 * 0.01_f64;
1947        let e1_actual = cache.interp_elevation(dist_m);
1948        assert_eq!(e1_actual, e1_expected);
1949    }
1950
1951    #[test]
1952    fn average_grade_over_range_is_correct() {
1953        let c = make_two_triangles_cycle();
1954        let cache = c.build_cache();
1955        let d0 = 40.0 * uc::M;
1956        let dd = 50.0 * uc::M;
1957        let expected0 = 0.01 * uc::R;
1958        let actual00 = c.average_grade_over_range(d0, dd, None);
1959        let actual00 = round(actual00.get::<si::ratio>(), Some(6)) * uc::R;
1960        assert_eq!(actual00, expected0);
1961        let actual01 = c.average_grade_over_range(d0, dd, Some(&cache));
1962        let actual01 = round(actual01.get::<si::ratio>(), Some(6)) * uc::R;
1963        assert_eq!(actual01, expected0);
1964    }
1965
1966    #[test]
1967    fn distance_to_next_stop_is_correct() {
1968        let c = make_two_triangles_cycle();
1969        let cache = c.build_cache();
1970        let d = 20.0 * uc::M;
1971        let expected = 20.0 * uc::M;
1972        let actual = c.calc_distance_to_next_stop_from(d, None);
1973        assert_eq!(actual, expected);
1974        let actual = c.calc_distance_to_next_stop_from(d, Some(&cache));
1975        assert_eq!(actual, expected);
1976        let d = 65.0 * uc::M;
1977        let expected = 25.0 * uc::M;
1978        let actual = c.calc_distance_to_next_stop_from(d, None);
1979        assert_eq!(actual, expected);
1980        let actual = c.calc_distance_to_next_stop_from(d, Some(&cache));
1981        assert_eq!(actual, expected);
1982        let d = 0.0 * uc::M;
1983        let expected = 40.0 * uc::M;
1984        let actual = c.calc_distance_to_next_stop_from(d, None);
1985        assert_eq!(actual, expected);
1986        let actual = c.calc_distance_to_next_stop_from(d, Some(&cache));
1987        assert_eq!(actual, expected);
1988    }
1989
1990    #[test]
1991    fn modifying_a_cycle_with_trajectory() {
1992        let c0 = make_two_triangles_cycle();
1993        let mut c = c0.clone();
1994        let n = 3;
1995        let d0 = 20.0; // units: m
1996        let v0 = 4.0; // units: m/s
1997        let dr = 65.0; // units: m
1998        let vr = 5.0; // units: m/s
1999        let dt = 10.0; // units: s
2000        let traj = ConstantJerkTrajectory::from_speed_and_distance_targets(n, d0, v0, dr, vr, dt);
2001        c.modify_by_const_jerk_trajectory(
2002            2,
2003            n,
2004            traj.jerk_m_per_s3 * uc::MPS3,
2005            traj.acceleration_m_per_s2 * uc::MPS2,
2006        );
2007        let expected = {
2008            let mut cyc = Cycle {
2009                name: String::from("Two Triangles"),
2010                init_elev: Some(0.0 * uc::M),
2011                time: vec![
2012                    0.0 * uc::S,
2013                    10.0 * uc::S,
2014                    20.0 * uc::S,
2015                    30.0 * uc::S,
2016                    40.0 * uc::S,
2017                    50.0 * uc::S,
2018                ],
2019                speed: vec![
2020                    0.0 * uc::MPS,
2021                    4.0 * uc::MPS,
2022                    traj.speed_at_step(1) * uc::MPS,
2023                    traj.speed_at_step(2) * uc::MPS,
2024                    5.0 * uc::MPS,
2025                    0.0 * uc::MPS,
2026                ],
2027                dist: vec![],
2028                grade: vec![
2029                    0.0 * uc::R,
2030                    0.0 * uc::R,
2031                    0.0 * uc::R,
2032                    0.0 * uc::R,
2033                    0.01 * uc::R,
2034                    0.01 * uc::R,
2035                ],
2036                elev: vec![],
2037                pwr_max_chrg: vec![],
2038                grade_interp: Default::default(),
2039                elev_interp: Default::default(),
2040                temp_amb_air: Default::default(),
2041                pwr_solar_load: Default::default(),
2042            };
2043            cyc.init().expect("initializaiton should not throw");
2044            cyc
2045        };
2046        assert_eq!(c.time.len(), expected.time.len());
2047        assert_eq!(c.speed.len(), expected.speed.len());
2048        assert_eq!(c.dist.len(), expected.dist.len());
2049        assert_eq!(c.grade.len(), expected.grade.len());
2050        for idx in 0..c.speed.len() {
2051            assert_eq!(c.time[idx], expected.time[idx]);
2052            assert_eq!(c.speed[idx], expected.speed[idx]);
2053            assert_eq!(c.dist[idx], expected.dist[idx]);
2054            assert_eq!(c.grade[idx], expected.grade[idx]);
2055        }
2056    }
2057
2058    #[test]
2059    pub fn modify_with_braking_trajectory() {
2060        let mut actual = {
2061            let mut cyc = Cycle {
2062                name: String::from("Test"),
2063                init_elev: Some(0.0 * uc::M),
2064                time: vec![
2065                    0.0 * uc::S,
2066                    1.0 * uc::S,
2067                    2.0 * uc::S,
2068                    3.0 * uc::S,
2069                    4.0 * uc::S,
2070                    5.0 * uc::S,
2071                ],
2072                speed: vec![
2073                    0.0 * uc::MPS,
2074                    4.0 * uc::MPS,
2075                    4.0 * uc::MPS,
2076                    1.0 * uc::MPS,
2077                    1.0 * uc::MPS,
2078                    0.0 * uc::MPS,
2079                ],
2080                dist: vec![],
2081                grade: vec![],
2082                elev: vec![],
2083                pwr_max_chrg: vec![],
2084                grade_interp: Default::default(),
2085                elev_interp: Default::default(),
2086                temp_amb_air: Default::default(),
2087                pwr_solar_load: Default::default(),
2088            };
2089            cyc.init().expect("initializaiton should not throw");
2090            cyc
2091        };
2092        let precision = Some(6);
2093        let (v_end, n_steps) =
2094            actual.modify_with_braking_trajectory((-4.0 / 3.0) * uc::MPS2, 3, Some(4.0 * uc::M));
2095        let v_end = round(v_end.get::<si::meter_per_second>(), precision);
2096        assert_eq!(v_end, 0.0);
2097        assert_eq!(n_steps, 3);
2098        let expected = {
2099            let n = 3;
2100            let d0 = 0.0;
2101            let v0 = 4.0;
2102            let dr = 4.0;
2103            let vr = 0.0;
2104            let dt = 1.0;
2105            let traj =
2106                ConstantJerkTrajectory::from_speed_and_distance_targets(n, d0, v0, dr, vr, dt);
2107            let mut cyc = Cycle {
2108                name: String::from("Test"),
2109                init_elev: Some(0.0 * uc::M),
2110                time: vec![
2111                    0.0 * uc::S,
2112                    1.0 * uc::S,
2113                    2.0 * uc::S,
2114                    3.0 * uc::S,
2115                    4.0 * uc::S,
2116                    5.0 * uc::S,
2117                ],
2118                speed: vec![
2119                    0.0 * uc::MPS,
2120                    4.0 * uc::MPS,
2121                    4.0 * uc::MPS,
2122                    traj.speed_at_step(1) * uc::MPS,
2123                    traj.speed_at_step(2) * uc::MPS,
2124                    traj.speed_at_step(3) * uc::MPS,
2125                ],
2126                dist: vec![],
2127                grade: vec![],
2128                elev: vec![],
2129                pwr_max_chrg: vec![],
2130                grade_interp: Default::default(),
2131                elev_interp: Default::default(),
2132                temp_amb_air: Default::default(),
2133                pwr_solar_load: Default::default(),
2134            };
2135            cyc.init().expect("initializaiton should not throw");
2136            cyc
2137        };
2138        assert_eq!(actual.time.len(), expected.time.len());
2139        for i in 0..actual.time.len() {
2140            let at = round(actual.time[i].get::<si::second>(), precision);
2141            let et = round(expected.time[i].get::<si::second>(), precision);
2142            let av = round(actual.speed[i].get::<si::meter_per_second>(), precision);
2143            let ev = round(expected.speed[i].get::<si::meter_per_second>(), precision);
2144            let ad = round(actual.dist[i].get::<si::meter>(), precision);
2145            let ed = round(expected.dist[i].get::<si::meter>(), precision);
2146            assert_eq!(at, et, "time@t={et}&i={i}");
2147            assert_eq!(av, ev, "speed@t={et}&i={i}");
2148            assert_eq!(ad, ed, "dist@t={et}&i={i}");
2149        }
2150    }
2151
2152    #[test]
2153    pub fn test_trim() {
2154        let c = make_two_triangles_cycle();
2155        let cyc = c.extend_time(Some(10.0 * uc::S), None);
2156        let dt_idle = cyc.ending_idle_time();
2157        assert_eq!(dt_idle, 10.0 * uc::S);
2158        // NOTE: extend_time adds time by 1.0 s increments so 10 points
2159        assert_eq!(cyc.time.len(), c.time.len() + 10);
2160        assert_eq!(*cyc.time.iter().last().unwrap(), 60.0 * uc::S);
2161        let cyc_trimmed = cyc.trim_ending_idle(None);
2162        assert_eq!(cyc_trimmed.time.len(), c.time.len());
2163    }
2164    type StructWithResources = Cycle;
2165
2166    #[test]
2167    fn test_resources() {
2168        let resource_list = StructWithResources::list_resources().unwrap();
2169        assert!(!resource_list.is_empty());
2170
2171        // verify that resources can all load
2172        for resource in resource_list {
2173            StructWithResources::from_resource(resource.clone(), false)
2174                .with_context(|| format_dbg!(resource))
2175                .unwrap();
2176        }
2177    }
2178
2179    #[test]
2180    fn test_resample() {
2181        let cyc0 = {
2182            let mut c = Cycle {
2183                name: String::from("a test"),
2184                time: vec![0.0 * uc::S, 10.0 * uc::S, 20.0 * uc::S],
2185                speed: vec![0.0 * uc::MPS, 10.0 * uc::MPS, 0.0 * uc::MPS],
2186                grade: vec![0.01 * uc::R, 0.01 * uc::R, -0.01 * uc::R],
2187                init_elev: None,
2188                dist: vec![],
2189                elev: vec![],
2190                pwr_max_chrg: vec![],
2191                temp_amb_air: vec![],
2192                pwr_solar_load: vec![],
2193                grade_interp: None,
2194                elev_interp: None,
2195            };
2196            c.init().unwrap();
2197            c
2198        };
2199        let cyc1 = cyc0.resample(1.0 * uc::S);
2200        assert_eq!(21, cyc1.time.len());
2201        assert_eq!(
2202            cyc1.time[cyc1.time.len() - 1],
2203            cyc0.time[cyc0.time.len() - 1]
2204        );
2205        assert_eq!(cyc1.time[0], cyc0.time[0]);
2206        assert_eq!(cyc1.time[0], 0.0 * uc::S);
2207        assert_eq!(cyc1.time[5], 5.0 * uc::S);
2208        assert_eq!(cyc1.speed[5], 5.0 * uc::MPS);
2209        assert_eq!(cyc1.grade[5], 0.01 * uc::R);
2210        assert_eq!(cyc1.time[10], 10.0 * uc::S);
2211        assert_eq!(cyc1.speed[10], 10.0 * uc::MPS);
2212        assert_eq!(cyc1.grade[10], 0.01 * uc::R);
2213        assert_eq!(cyc1.time[11], 11.0 * uc::S);
2214        assert_eq!(cyc1.speed[11], 9.0 * uc::MPS);
2215        assert_eq!(cyc1.grade[11], -0.01 * uc::R);
2216        assert_eq!(cyc1.time[20], 20.0 * uc::S);
2217        assert_eq!(cyc1.speed[20], 0.0 * uc::MPS);
2218        assert_eq!(cyc1.grade[20], -0.01 * uc::R);
2219    }
2220}
2221
2222lazy_static! {
2223    pub static ref CYC_ACCEL: Cycle = Cycle::try_from(CycleBuilder {
2224        name: String::from("accel test"),
2225        time: (0..300)
2226            .map(|t| (t as f64) * uc::S)
2227            .collect::<Vec<si::Time>>(),
2228        speed: vec![90.0 * uc::MPH; 300],
2229    })
2230    .unwrap();
2231}