Skip to main content

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