Converting an equation with units to one without units

Question

I may have a solution but it is not slick! See below

Doing numerical calculations with functions that use units can be very slow. I would like to have a function or operator that can turn my functions into unitless functions. I have a partial solution but it doesn't work in all situations.

The blackbody radiation equation:

Clear[BB];
BB[l_Quantity, T_Quantity] := (
  2 Quantity["PlanckConstant"] Quantity[
    "SpeedOfLight"]^2)/(l^5 Quantity[1, "Steradians"]) 1/(
  Exp[Quantity[
      "PlanckConstant"] Quantity[
       "SpeedOfLight"]/(l Quantity["BoltzmannConstant"] T)] - 
   1)

Now if I want to do many calculations with this equation it is very slow as can be seen from just making a plot:

AbsoluteTiming[
  Plot[QuantityMagnitude[
    UnitConvert[
     BB[Quantity[l, "Micrometers"], Quantity[333, "Kelvins"]],
      "Microflicks"]], {l, .1, 15}
   , Frame -> True
   , FrameLabel -> {"\[Lambda] (\[Mu]m)", 
     "\[ScriptCapitalL] (\[Mu]flicks)"}
   , ImageSize -> Medium
   ]
  ]

Please keep in mind that I want to do much more than just make plots quickly. I need to do Monte Carlo simulations where I will make compiled versions of my equations. Also in my real work, the equations are much more complicated and automatically generated.

The function that I have come up with that works some of the time is the following:

Clear[NoUnits]
NoUnits[fun_, 
   outunits_] := ((UnitConvert[#, outunits] & /@ 
       fun) /. {UnitConvert[a_, outunits] :>  a}) /. {Quantity[q_, 
      u_] :> q};

I can apply this function to BB[l,T] to get:

Clear[BBnu];
BBnu[l_, T_] := 
 Evaluate[NoUnits[
   BB[\[Lambda], 
     T] /. {l -> Quantity[l, "Micrometers"], 
     T -> Quantity[T, "Kelvins"]}, "Microflicks"]]

Where I have to tell my NoUnits[] function what units to use for each variable and the units I want for the final output. The new function BBnu[l,T] is just a function of numbers and gives a numerical output. Now the calculations go much faster!

AbsoluteTiming[Plot[BBnu[l, 333], {l, .1, 15}
  , Frame -> True
  , FrameLabel -> {"\[Lambda] (\[Mu]m)", 
    "\[ScriptCapitalL] (\[Mu]flicks)"}
  , ImageSize -> Medium
  ]
 ]

But now if I try my function NoUnits[] on something more complex like a linear combination of blackbodies it fails:

NoUnits[(a BB[l, T1] + b BB[l, T2]) /. {l -> 
    Quantity[l, "Micrometers"], T1 -> Quantity[T1, "Kelvins"], 
   T2 -> Quantity[T2, "Kelvins"]}, "Microflicks"]

I have tried what I have learned from: How to create fast functions based on units? Converting the units of an equation or expression (how to become 3 times heavier…) Yet these do not solve my full problem.

Edit The following are some stressing testing functions:

The blackbody equation again B[l,T]:

BB[l_Quantity, T_Quantity] := (
   2 Quantity["PlanckConstant"] Quantity["SpeedOfLight"]^2)/(
   l^5 Quantity[1, "Steradians"]) 1/(
   Exp[Quantity[
       "PlanckConstant"] Quantity[
        "SpeedOfLight"]/(l Quantity["BoltzmannConstant"] T)] - 1);

A made up trasnsmsion equation trans[l]:

trans[l_Quantity] := 
  1/(1 + Exp[- 
       Quantity[1, ("Micrometers")^-1] (l - 
        Quantity[4, "Micrometers"])]) 1/(
   1 + Exp[Quantity[3, (
       "Micrometers")^-1] (l - Quantity[10, "Micrometers"])]);

A madeup source equation that includes an interpolation function with compound units Ladd[l,T]:

With[{
   fun = Interpolation[{{4/300, 0}, {5/300, 
       1 10^4}, {6/300, .5 10^6}, {8/300, 8 10^6}, {10/300, 
       4 10^6}, {12/300, 2 10^6}, {14/300, 1 10^6}, {15/300, 0}}, 
     InterpolationOrder -> 3]
   },
  Ladd[l_Quantity, T_Quantity] := Piecewise[{
     {Quantity[0, 
       "Kilograms" ("Meters")^-1 ("Seconds")^-3 ("Steradians")^-1], 
      l/T <= Quantity[4, "Micrometers"]/Quantity[300, "Kelvins"]}
     , {Quantity[
       fun[QuantityMagnitude[
         UnitConvert[l, "Micrometers"]/UnitConvert[T, "Kelvins"]]], 
       "Kilograms" ("Meters")^-1 ("Seconds")^-3 ("Steradians")^-1], 
      Quantity[4, "Micrometers"]/Quantity[300, "Kelvins"] < l/T <= 
       Quantity[15, "Micrometers"]/Quantity[300, "Kelvins"]}
     , {Quantity[0, 
       "Kilograms" ("Meters")^-1 ("Seconds")^-3 ("Steradians")^-1], 
      l/T > Quantity[15, "Micrometers"]/Quantity[300, "Kelvins"]}
     }, Quantity[0, 
     "Kilograms" ("Meters")^-1 ("Seconds")^-3 ("Steradians")^-1]]
  ];

Another made up function: LaddTwo[l]

With[{
   fun = Interpolation[{{4, 0}, {5, 1 10^4}, {6, .5 10^6}, {8, 
       8 10^6}, {10, 4 10^6}, {12, 2 10^6}, {14, 1 10^6}, {15, 0}}, 
     InterpolationOrder -> 3]
   },
  LaddTwo[l_Quantity] := Piecewise[{
     {Quantity[0, 
       "Kilograms" ("Meters")^-1 ("Seconds")^-3 ("Steradians")^-1], 
      l <= Quantity[4, "Micrometers"]}
     , {Quantity[
       fun[QuantityMagnitude[UnitConvert[l, "Micrometers"]]], 
       "Kilograms" ("Meters")^-1 ("Seconds")^-3 ("Steradians")^-1], 
      Quantity[4, "Micrometers"] < l <= Quantity[15, "Micrometers"]}
     , {Quantity[0, 
       "Kilograms" ("Meters")^-1 ("Seconds")^-3 ("Steradians")^-1], 
      l > Quantity[15, "Micrometers"]}
     }, Quantity[0, 
     "Kilograms" ("Meters")^-1 ("Seconds")^-3 ("Steradians")^-1]]
  ];

and finally, a full system equation that brings it all together Lfull[l,T]:

Lfull[l_Quantity, T_Quantity, a_, b_] := 
  a trans[l] BB[l, T] + b trans[l] Ladd[l, T]+ a b LaddTwo[l];

I'm not trying to make this overly completed but this is the smallest full-featured example.

The goal would be to get this "units to unitless" operator to work on the Lfull[l,T,a,b] function.

My poor solution I'm still very interested in a better solution! I can think of some cases where this may fail. Update: It has not failed yet but I still do not like the implementation. I have updated my function complete with some documentation:

    Clear[NoUnits];
NoUnits::usage = 
  "NoUnits[fun, rp, outunits] Will take a function fun that is \
defined with units and convert it into just a numerical function. Of \
course the user is now responsible for providing the numbers in the \
correct \"Units\".   rp is a list of variables and the corresponding \
units the user wants to use for the resulting function.  The units on 
 the output is the outunits.
  
  Options:
  \"MaxPower\" What maximum power to use when looking for match \
units. Default is 3. If you have something like l^5/a^2 in your \
equation you may wat to use 5 for your maximum power.
  
  Example:
  fun[v_Quantity,m_Quantity,c_]:= c m \!\(\*SuperscriptBox[\(v\), \(2\
\)]\);
  Kinetic energy equation were we have to put in c=0.5;
  Test it for some inputs:
  
 e=fun[Quantity[2,\"Meters\"/\"Seconds\"],Quantity[12,\"Kilograms\"],\
.5];
  UnitConvert[e,\"Ergs\"]
  \!\(\*TemplateBox[{\"2.4`*^8\",\"\\\"ergs\\\"\",\"ergs\",\"\\\"Ergs\
\\\"\"},\n\"Quantity\"]\)
  Lets say we want just a numbers equation where we input mph, pounds \
and want Ergs out Then do:
  eq=NoUnits[fun[v,m,c],{{v,\"Miles\"/\"Hours\"},{m,\"Pounds\"}},\"\
Ergs\"]
  \!\(\*FractionBox[\(88523495162333\\\ c\\\ m\\\ \
\*SuperscriptBox[\(v\), \(2\)]\), \(97656250\)]\)
  Test it:
  eq/.{c\[Rule].5,m\[Rule] 10, v\[Rule] 55}
  1.3710518930742134`*^10
  Use the original equation:
  e=fun[Quantity[55,\"Miles\"/\"Hours\"],Quantity[10,\"Pounds\"],.5];
  UnitConvert[e,\"Ergs\"]
  \!\(\*TemplateBox[{\"1.3710518930742136`*^10\",\"\\\"ergs\\\"\",\"\
ergs\",\"\\\"Ergs\\\"\"},\n\"Quantity\"]\)";
Options[NoUnits] = {
   "MaxPower" -> Automatic
   };
NoUnits[fun_, rp_List, outunits_, opts : OptionsPattern[NoUnits]] := 
  Module[{eq, list, u, pow, unit, a, b, mp},
   mp = OptionValue["MaxPower"];
   If[Not[NumericQ[mp] && IntegerQ[mp]], mp = 3];
   If[mp <= 0, mp = 3];
   eq = fun /. 
     Evaluate[Rule[#[[1]], Quantity[#[[1]], #[[2]]]] & /@ rp];
   eq = eq /. {Quantity[a_, b_] :> 
        UnitConvert[Quantity[a, b]]} /. {Quantity[a_, 
        "PureUnities"] :> a, Quantity[a_, "DimensionlessUnit"] :> a};
   eq = eq /. {Quantity[a_, b_] /; CompatibleUnitQ[b, outunits] :> 
       QuantityMagnitude[UnitConvert[Quantity[a, b], outunits]]};
   
   list = 
    DeleteDuplicates[Cases[eq, Quantity[a_, b_] :> b, \[Infinity]]];
   u = rp[[All, 2]];
   pow = Tuples[Range[-mp, mp], Length[u]];
   pow = SortBy[pow, Total[Abs[#]] &];
   pow = {u, #}\[Transpose] & /@ pow;
   Do[
    unit = 
     Select[Times[
          Sequence @@ #] & /@ (pow /. {{a_String, b_} :> 
            Power[a, b]}), CompatibleUnitQ[#, ele] &, 1][[1]];
    eq = eq /. {Quantity[a_, b_] /; CompatibleUnitQ[b, unit] :> 
        QuantityMagnitude[UnitConvert[Quantity[a, b], unit]]};
    , {ele, list}];
   
   Return[eq]
   ];

Here is an example of the usage:

Answer 1

Instead of using UnitConvert, you can just divide the output by the relevant output unit. Here's one implementation of the idea, based on the OP's version:

ClearAll @ NoUnits;
SetAttributes[NoUnits, HoldAll];

NoUnits[fun_, quants : {__Quantity}, output : Quantity[1, out_]] := 
   Expand[fun/output /. Thread[quants[[All, 1]] -> quants]] /. x_ y_Piecewise :> PiecewiseExpand[x y] /. Quantity[x_, _] :> x

Usage:

LFullNU[l_, T_, a_, b_] = NoUnits2[Lfull[l, T, a, b], {Quantity[l, "Micrometers"], Quantity[T, "Kelvins"]}, Quantity[1, "Microflicks"]];

Seems to work. I have as inputs the full Quantitys to make the replacements easier.

As a check, you can try

BBNU[l_, T_] = NoUnits2[BB[l, T], {Quantity[l, "Micrometers"], Quantity[T, "Kelvins"]}, Quantity[1, "Microflicks"]];

which yields the same result as in the OP.

Answer 2

Based on Mathematica documentation, and alternative solution could be implement a simple replacement rule as follows in the next function:

NoUnits[f_] := f /. q_Quantity :> QuantityMagnitude[q]

Then, implement the Plot as a pure function to increase drastically the execution time:

    AbsoluteTiming[Plot[
    #,
    {l, .1, 15}, Frame -> True, 
    FrameLabel -> {"\[Lambda] (\[Mu]m)", 
      "\[ScriptCapitalL] (\[Mu]flicks)"}, ImageSize -> Medium]] &@ NoUnits[BB[Quantity[l, "Micrometers"], Quantity[333, "Kelvins"]]]

The result will be:

Answer 3

The issue is that your definition of NoUnits uses Map which only maps at level one.

An alternative solution is:

Clear[NoUnits]
NoUnits[fun_, 
  outunits_] := ((fun /. 
      HoldPattern[Quantity[x__]] :> 
       UnitConvert[Quantity[x], outunits]) /. {UnitConvert[a_, 
       outunits] :> a}) /. {Quantity[q_, u_] :> q}

Answer 4

If already in the input of the actual notebook: KnownUnitQ["Microflicks"] gives True. So there is need for someaction to get this referentially connceted to other known unit entities.

Mathematica V 12.0.0 has FormularLookup and FormulaData for lot of the questions purposes.

FormulaLookup["planck's law"]

The unit system in Physics is in need to be closed both operatinal and pur informal. As shown on the documentation page for Quantity.

(* True *)

Normal[Quantity[10, "Percent"]] (* 1/10 *)

UnitDimensions["PartsPerMillion"] (* {} *)

For example

N[Quantity[1, "GravitationalConstant"], 20] (* 1 G *)

UnitConvert[N[Quantity[1, "GravitationalConstant"], 20]] (* Quantity[6.674*10^-11, ("Meters")^3/("Kilograms" ("Seconds")^2)] *)

This input does the required job:

equation = FormulaData[{"PlanckRadiationLaw", "Wavelength"},
    {"T" -> Quantity[5000, "Kelvins"], 
     "\[Lambda]" -> Quantity[wl, "Micrometers"]}
    ][[2, 2]];

Plot[equation, {wl, 0.1, 5}, 
 AxesLabel -> {"Wavelength [\[Mu]m]", 
   "Spectral radiance [W \!\(\*SuperscriptBox[\(sr\), \
\(-1\)]\)\!\(\*SuperscriptBox[\(m\), \(-3\)]\)]"}]

Like with the Unit it is not necessary that the forumla is a Physics one. If the running kernel knows it suffices.

To works with Units or Quantity is already optimized as a built-in. It is generalzied to match the requirement of a physical CAS completely. You have QuantityUnit and QuantityMagnitude for operatinal optimization. QuantityMagnitude works always for separation of magnitude from unit at most for visualization purposes.

You have only this built-in at hand. The solutions of the others answerer make use of that by the use of Quantity.

On the documentation page of QuantityMagnitude is an section about how Mathematica deals with unknown units. It is interpreted semantically.

QuantityMagnitude[Quantity[100, "Yards"], "Meterss"]

(* 2286/25 *)

There is the built-in IndependentUnit represents a unit string with no relationship to other units within a Quantity.

KnownUnitQ returns True for valid IndependentUnit specifications:

KnownUnitQ[IndependentUnit["Flicks"]]

(* True *)

Names of untis have to be unique. And can be checked with CompatibleUnitQ. UnitDimensions is the logical center of the system:

(* {{"ElectricCurrentUnit", -1}, {"LengthUnit", 2}, {"MassUnit", 1}, {"TimeUnit", -3}} *)

If this properly resolves Your Flicks You can do everything even get the input QuantityUnits out of the equation. And You can check this with it too.

But

UnitDimensions["Flicks"]

gives

{{"AngleUnit", -2}, {"LengthUnit", -1}, {"MassUnit", 
  1}, {"TimeUnit", -3}}

This proves my first solution in depths.

But take are:

UnitDimensions["1"]

(* UnitDimensions["1"] *)

So Your target is

UnitDimensions[1]

(* {} *)

**Numerical values are considered dimensionless!**

Simplifyication works this way:

UnitSimplify!

Interrogation of units works this way:

Composition[QuantityUnit, UnitConvert, Quantity] /@ {"PlanckConstant",
   "BoltzmannConstant", "MolarGasConstant", "AvogadroNumber"}

{("Kilograms" ("Meters")^2)/("Seconds"), ("Kilograms" ("Meters")^2)/(
 "Kelvins" ("Seconds")^2), ("Kilograms" ("Meters")^2)/(
 "Kelvins" "Moles" ("Seconds")^2), "DimensionlessUnit"}

The selection of fundamental constants can be replace by Your interests.

Good practice to work with units in graphical representation is

s = Quantity[100, "Meters"] - Quantity[9.8, "Meters/Seconds^2"]*t^2;
factor = QuantityMagnitude[Quantity[1, "Meters"], "Feet"];
Plot[s*factor, {t, Quantity[0, "Seconds"], Quantity[3, "Seconds"]}, 
 FrameLabel -> Automatic, PlotTheme -> "Web"]

from this answers: how do i properly use quantities units in plots?.

To get more free space for units make use of getting useful units for combinations of physical constants like on wolframalpha.

This should do the job: specify set of base units to use in unitconvert as I understand the question.

Implementation:

(*a set of standard units that are used when not specified*)
siUnits = {"Seconds", "Meters", "Kilograms", "Kelvins", 
   "KelvinsDifference", "Amperes", "Candelas", "Moles", "Radians"};
siUnitDimensions = UnitDimensions[#][[1, 1]] & /@ siUnits;

makeUnitSystem::dependent = 
  "The unit system `1` is overcomplete. Please remove some unit.";
makeUnitSystem[] = Thread[siUnitDimensions -> siUnits];
makeUnitSystem[L_List] := 
 Module[{M, n, u},(*convert the desired unit system to base units*)
  M = Lookup[#, siUnitDimensions, 0] & /@ 
    Apply[Rule, UnitDimensions /@ L, {2}];
  If[MatrixRank[M] < Length[L], 
   Message[makeUnitSystem::dependent, L];
   Return[$Failed]];
  (*check which base units cannot be expressed in this system*)
  n = Position[Diagonal[PseudoInverse[M].M], Except[1], {1}, 
    Heads -> False];
  (*extend the unit system if necessary*)
  If[Length[n] > 0, 
   Return[makeUnitSystem[Append[L, siUnits[[n[[1, 1]]]]]]]];
  (*find the compound units that represent the base units*)
  u = Times @@@ Transpose[L^Transpose[PseudoInverse[M]]];
  (*return replacement list*)Thread[siUnitDimensions -> u]]

unitConvert[x_Quantity, 
  unitSystem_ /; VectorQ[unitSystem, Head[#] === Rule &]] := 
 UnitConvert[x, Times @@ Power @@@ (UnitDimensions[x] /. unitSystem)]

With this You may check whether the dimensionless parameters suit into the systemes international de unites for a complete system. This is inspired by the answer of Roman Maeder, @roman.

This makes a consistent dimensionless system of dimensions provable. I does not make the choice. It does not make any calculations faster. This is done internally and at fastest if consistency is True.

This is for converting units, an example:

ClearAll[withUnits];
SetAttributes[withUnits, HoldAll];
withUnits[code_] :=
  Function[Null,
     Block[{Quantity},
       SetAttributes[Quantity, HoldRest];
       Quantity /: UnitConvert[arg_, Quantity[_, unit_]] :=
          UnitConvert[arg, unit];
       Quantity /: Times[0, Quantity[_, unit_]] :=
          Quantity[0, unit];
       With[{
          m = Quantity[1, "Meters"], 
          s = Quantity[1, "Seconds"],       
          min =  Quantity[1, "Minutes"],
          km = Quantity[1, "Kilometers"]
        },
       #]],
    HoldAll][code];

withUnits[UnitConvert[1 m/s^2*(1 min)^2,km]]

from simpler input for the new unit support by @leonid-shifrin.

k[uqty_] := Quantity[uqty, "Kelvins"]
um[uqty_] := Quantity[uqty, "Micrometers"]

FormulaData[{"PlanckRadiationLaw", "Wavelength"}, 
{"T" -> k@5000, "\[Lambda]" -> um@\[Lambda]}]

gives a formular in dimensionless [Lambda]: