aurora-opensource / au

A C++14-compatible physical units library with no dependencies and a single-file delivery option. Emphasis on safety, accessibility, performance, and developer experience.
Apache License 2.0
329 stars 21 forks source link

Question: is there a way to cast a scalar to a `Quantity`? #248

Closed pdaoust closed 5 months ago

pdaoust commented 5 months ago

I'm writing a template function that takes a value (could be a Quantity, could be a bare naked scalar, trying to be as generic as possible) and does some mathy operations on it before casting it back to the type it was before. I don't know the value's type, but it's a template param so it's known at compile time.

Question: Does Quantity have any cast operator that lets me do the mathy stuff that potentially yields a bare scalar (e.g., Quantity<U, R> / Quantity<U, R>Rthen cast it back into the original Quantity<U, R> type? I know I'm verging into unsafe territory, but I'm writing a library for third-party consumption and don't want to assume they want to use au. (Though, I mean, who wouldn't want to use a units library? c'mon)

(edit: and yes, I'm currently quesitoning my assumption that a unitless value isn't, in fact appropriate in the context, but in the meantime I'll let this question stay up)

chiphogg commented 5 months ago

Hmm... I'm not quite sure I understand your use case. Could you flesh it out a bit more, maybe with an example?

It seems to me that if the operation results in a scalar, and you want to cast it back to its original Quantity type, then either the original unit U was equivalent to "Unos" (i.e., a dimensionless unit of magnitude 1), or the computation contains an error and it's good that we prevent it.

There is a way to get a (possibly mutable) reference to the underlying data, if that's what you want (but I don't think it is).

Also, there is a way to turn a raw number into a Quantity type that is "unitless" (i.e., dimensionless and of magnitude 1), if that's what you want (but I don't think it is).

Let me know! :slightly_smiling_face:

pdaoust commented 5 months ago

Thanks! I found that C++ and AU are both smart enough to do the right sort of inference if I use a lot of auto in my function. But for purpose of discussion, to help me understand this better, here's the thing I'm trying to do: a function that returns a function that can be used to calculate a moving average (single-pole IIR filter) over an iterable of measured values Iter<TMeas> and an iterable of timestamps Iter<TTime>. (Note that Iter is my own invention and is not a std::iterator or anything like that -- and that's just a stand-in for something more complicated)

Note that TMeas, TTime, and TInterval could be bare scalars, or AU units, or (in the case of TTime and TInterval) std::chrono::time_point and std::chrono::duration.

I had wanted to replace the auto in the alpha term with an actual explicit type that's inferred from the template parameters or something, because I'm allergic to auto whenever I can be more specific, buuuut... maybe I'm just being pedantic, because the below code works with auto.

template <typename TMeas, typename TTime, typename TInterval>
Iter<TMeas> createExponentialMovingAverageIter(
  Iter<TMeas> values,
  Iter<TTime> timestamps,
  TInterval timeConstant
) {
  // Combine the two iterables into an iterable of tuples.
  Iter<std::tuple<TMeas, TTime>> timestampedValues = Iter::zip(values, timestamps);

  // Calculate the moving average.
  Iter<std::tuple<TMeas, TTime>> movingAverage = Iter::reduce(
    values,
    [timeConstant](
      std::tuple<TMeas, TTime> prev,
      std::tuple<TMeas, TTime> next
    ) {
      TMeas prevValue = std::get<0>(prev);
      TTime prevTime = std::get<1>(prev);
      TMeas nextValue = std::get<0>(next);
      TTime nextTime = std::get<1>(next);

      TInterval timeDelta = nextTime - prevTime;

      // If TTime is unsigned long (common in Arduino projects)
      // TInterval and alpha will be unsigned long too.
      // If TTime is au::QuantityPoint<au::Seconds, float>
      // then TIntervals will be au::Quantity<au::seconds, float>
      // and the type of alpha will be (ideally) a bare float after one interval is divided by the other.
      auto alpha = 1 - pow(M_E, -timeDelta / timeConstant);

      // The type of integrated should also be TMeas, even though alpha is scalar.
      TMeas integrated = prevValue + alpha * (nextValue - prevValue);
      return std::tuple<TMeas, TTime>(integrated, nextTime);
    }
  );

  return Iter::map(movingAverage, [](std::tuple<TMeas, TTime> value) { return std::get<0>(value); });
}
chiphogg commented 5 months ago

OK, so I'm assuming that the auto alpha = ... line is the only line we're focusing on here --- the rest is for context --- and that replacing auto with a more descriptive type is the main goal. I'll proceed under those assumptions; let me know if I've missed the mark!

First: the comment says that alpha will (in some situations) be unsigned long. How can that be? It's basically $1 - e^{-x}$, where $x > 0$. If we store this result in a type unsigned long, won't it always be either exactly 0u or 1u? I assume the type for alpha will be basically dictated by the return type for pow, which, if it's the C function, is presumably double.

Maybe the comments about the type of alpha were meant to apply to the type of the second argument of pow. If that's right, then we're good already, because Au currently produces a raw numeric value as the result of dividing two Quantity types with the same unit.

NOTE: Actually, we plan to change this in a future version --- see #185! But even in this case, I think it would still work, because a Quantity whose unit is "unitless" (unos) will implicitly convert to its underlying type. There is a chance this would not trigger, but, as per the latest comment in that issue, that would be easy to fix because we could provide an as_raw_number(...) API that would work for both Quantity and non-Quantity types. So the tl;dr is that this should continue to work both now and indefinitely, with the worst possibility being potentially needing to add a wrapping function call when upgrading to a future version of Au.

Finally, there is one more wrinkle that your code snippet suggests, but which you evidently haven't encountered. If your argument types are integer-backed Quantity types, then the division will not work, because integer division is dangerous and error-prone when the denominator has units (see links in the OP for #249 for examples). It will be a compile-time error that tells you to use the integer_quotient function instead. Obviously, that solution wouldn't work in a generic context like this! In this case, it occurs to me that we can and perhaps should provide a carve-out and permit integer division when the units are identical, because it's hard for me to see how this could be problematic. I'm tracking that in #249.

Is there more help I can provide on this present issue?

pdaoust commented 5 months ago

OK, so I'm assuming that the auto alpha = ... line is the only line we're focusing on here --- the rest is for context --- and that replacing auto with a more descriptive type is the main goal.

yeah, sorry, I've been known to give a little too much context...

First: the comment says that alpha will (in some situations) be unsigned long. How can that be? It's basically $1 - e^{-x}$, where $x > 0$. If we store this result in a type unsigned long, won't it always be either exactly 0u or 1u?

Well crap, you're right. Before I refactored this function to use template params, I was indeed casting those values to floats. Maybe I need to go back to explicit casting to floating point again, just for that one line, cuz we can't leave it up to the API consumer's discretion. I think that resolves the bigger question that my casting question shadows.

I definitely encountered the integer division wrinkle too, and in my test code I'm using converting unsigned long timestamps to float just to get around it. Thanks so much for using this code snippet as motivation for supporting integer division for cases where it's safe!

Is there more help I can provide on this present issue?

Nope, you've been a super big help and prodded me to improve my codebase!

pdaoust commented 5 months ago

by the way, I don't know if au is a hobby project or if you're able to spend your company 20% time on it or something, but either way this has been an amazing amount of support I've experienced. Not just helpful but swift too. I want to get to a place where we offer the same level of dev support in our (commercially funded) open source project. Thank you!

pdaoust commented 5 months ago

Uh-oh, I've got a gnarly one, not quite sure how to sort it out and stay within safe unit math! I'm hoping you can counsel me on the 'right' way to do this while sticking to my goal of using template params. This is for a PID loop; I'll just show the guts of the calculation.

      // values.setpoint and values.processVariable may be scalars,
      // or they may be au::QuantityPoint<au::Kelvins, float> or something like that.
      // That means TError is au::Quantity<au::Kelvins, float>
      TError error = values.setpoint - values.processVariable;

      // TInterval is likely to be std::chrono::steady_duration, or it could be au::Quantity<au::Seconds, long long>
      TInterval timeDelta = values.timestamp - prevState.timestamp;
      // Similar spelling, but a different type. If using AU units it might end up being au::Quantity<au::UnitProduct<au::Kelvins, au::Seconds>, float>
      TIntegral integral = prevState.integral + error * timeDelta;
      // Similar thing going on here. So far all is good; we've got some composite units, but we can handle that.
      TDeriv derivative = (TCalc)(error - prevState.error) / (TCalc)timeDelta;
      // WHOA! Kp, Ki, and Kd are scalars, so that part is okay, but how do you add a temperature quantity, a temperature*duration quantity, and a temperature/duration quantity?
      auto control = values.weights.Kp * error + values.weights.Ki * integral + values.weights.Kd * derivative;
chiphogg commented 5 months ago

by the way, I don't know if au is a hobby project or if you're able to spend your company 20% time on it or something, but either way this has been an amazing amount of support I've experienced. Not just helpful but swift too. I want to get to a place where we offer the same level of dev support in our (commercially funded) open source project. Thank you!

We don't really have 20% time --- we're laser focused on turning the Aurora driver into an actual product that is worthy of running on public roads. I do spend some work time on this project, but it's in the low single digits of percent. :slightly_smiling_face: Sometimes you'll get lucky because I use GitHub notifications for my main work (we use GitHub enterprise), and if I happen to see one from Au, and it happens to be something I understand pretty well, then I can dash off a quick response. But most of my tinkering on the project is in personal hobby project time.

Uh-oh, I've got a gnarly one, not quite sure how to sort it out and stay within safe unit math! I'm hoping you can counsel me on the 'right' way to do this while sticking to my goal of using template params. This is for a PID loop; I'll just show the guts of the calculation.

      // values.setpoint and values.processVariable may be scalars,
      // or they may be au::QuantityPoint<au::Kelvins, float> or something like that.
      // That means TError is au::Quantity<au::Kelvins, float>
      TError error = values.setpoint - values.processVariable;

      // TInterval is likely to be std::chrono::steady_duration, or it could be au::Quantity<au::Seconds, long long>
      TInterval timeDelta = values.timestamp - prevState.timestamp;
      // Similar spelling, but a different type. If using AU units it might end up being au::Quantity<au::UnitProduct<au::Kelvins, au::Seconds>, float>
      TIntegral integral = prevState.integral + error * timeDelta;
      // Similar thing going on here. So far all is good; we've got some composite units, but we can handle that.
      TDeriv derivative = (TCalc)(error - prevState.error) / (TCalc)timeDelta;
      // WHOA! Kp, Ki, and Kd are scalars, so that part is okay, but how do you add a temperature quantity, a temperature*duration quantity, and a temperature/duration quantity?
      auto control = values.weights.Kp * error + values.weights.Ki * integral + values.weights.Kd * derivative;

Now this is a really interesting one! The inescapable fact is that Kp, Ki, and Kd all have units --- and that they necessarily all have different units. Those units are important for the semantic understanding of those coefficients, so in a sense, it's good that this is failing to compile. The only two ways to reconcile this are to change their types to reflect those units, or to treat the quantities as raw scalars and discard all unit information.

The former approach is more appealing, and it's probably easier than it sounds. Let's say the desired type of control is some typename Control. Then the declarations could use types that are something like:

using KpType = decltype(std::declval<Control>() / std::declval<TError>());
using KiType = decltype(std::declval<Control>() / std::declval<TIntegral>());
using KdType = decltype(std::declval<Control>() / std::declval<TDeriv>());

I think if you do something like this, it should give you all the control you need. (I'm assuming you have control of the definition of the type of values.) If I've missed any context, let me know!

pdaoust commented 5 months ago

huh, that's a really interesting idea to chew on. Kp, Ki, and Kd all have units. (Curious if you've been using PID loops in your work on self-driving semis!) I just sorta pictured the coefficients as unitless scalars -- arbitrary numbers. I can certainly see how a type of Control / TError etc would cause the calculations to be 'safe' when you then multiply it with a value of TError, just wasn't sure if that was uncovering some inherent semantics in the coefficients that I hadn't noticed before, or if it was just getting around the type system somehow.

The type of Control is indeed under my control (or rather, the control of the consumer of my API) as a template param. In fact I used to have that as a template param and removed it. But AFAICT it should usually be unitless -- the only meaning I can think of is that it'll often be a duty cycle of a PWM, so I guess that might mean a consumer would want an au percentage unit.

Anyhow, I appreciate the time you've taken to help me understand this. Part of the issue is that I'm also learning C++ at the same time that I'm creating this API, so I'm learning a huge pile and I'm not always storing my learnings in a tidy fashion in my head. I love type systems, but I'm more accustomed to Rust and C# than C++. How to do generic/template params well is a bit of a mind-bender right now.

[EDIT] Well heck, thinking about the calculus, I'm seeing it now... makes total sense that

the PID loop is in its happy point when all three terms are near zero. But I'm still trying to wrap my head around the idea of what the units of the coefficients mean. I'm on my way though!

chiphogg commented 5 months ago

Really glad I could help! And surprised that you're new to C++ --- it doesn't show.

I'm closing this before our team meeting because I think it's resolved, but please feel free to reopen if you need anything else.

pdaoust commented 5 months ago

well, I may be new to C++ but I'm finding a lot of the concepts from Rust and C# carry over. In fact, I'm surprised just how Rusty C++'s lifetime semantics are. And I'm a big fan of type systems as programming languages that 'make impossible situations impossible' in the words of my colleauges :) Thanks for the compliment! (And no, this comment doesn't mean I want to reopen the issue 😄 )