3
$\begingroup$

I want to make UnitConvert operate on custom units with special rules.

For instance, I would like to be able to convert between "Hartrees" (or any energy unit) and "Wavenumbers". This conversion is based off of assuming we're working with light, in which case there is a straightforward conversion between the two using the speed of light and Planck's constant.

Optimally, I could do something like:

UnitConvert[Quantity[..., "Hartrees"], "Wavenumbers"]

but of course Mathematica balks as the UnitDimensions aren't actually the same.

How can I do this? I know it won't be easy, but it at least ought to be possible.

$\endgroup$
3
$\begingroup$

I spent a while spelunking and came up with a solution. We can basically just dig through what UnitConvert is doing and although there are many pitfalls, some of which I probably have not avoided, I do have a solution.

Here goes. First, the most fundamental step is a function to register a new unit rule. There are many symbols involved here, so I'll post the code first then dissect it:

addUnitConvertRule[key_ -> Quantity[val_, unit_]] :=
 With[
  {
   customKey = Symbol["CalculateUnits`UnitCommonSymbols`" <> key]
   },
  If[Not@KeyExistsQ[QuantityUnits`Private`$UnitReplacementRules, 
     key],
   AppendTo[QuantityUnits`Private`$UnitReplacementRules, 
    key -> customKey]
   ];
  CalculateUnits`UnitCommonSymbols`KnownUnit0Q[customKey] = True;
  With[{uc = UnitConvert[unit]},
   With[{uval =
      val*QuantityMagnitude[uc]*QuantityUnit[uc] /. {
         s_String :> 
          Lookup[QuantityUnits`Private`$UnitReplacementRules, s]
         } //. u_Symbol?CalculateUnits`UnitCommonSymbols`KnownUnit0Q :>
        Replace[
         CalculateUnits`UnitCommonSymbols`UnitLookup[u, 
          "FundamentalUnitValue"],
         _CalculateUnits`UnitCommonSymbols`UnitLookup :> u
         ]
     },
    CalculateUnits`UnitCommonSymbols`UnitLookup[customKey, 
      "FundamentalUnitValue"] =
     uval;
    CalculateUnits`UnitTable`CUF0[customKey] =
     Evaluate[#*uval /. 
        s_Symbol?CalculateUnits`UnitCommonSymbols`KnownUnit0Q :> 
         1] &;
    CalculateUnits`UnitTable`CUF1[customKey] =
     Evaluate[#/uval /. 
         s_Symbol?CalculateUnits`UnitCommonSymbols`KnownUnit0Q :> 
          1]*# &
    ]
   ];
  CalculateUnits`UnitCommonSymbols`UnitLookup[customKey, 
    "UnitDimensions"] =
   With[{uc = UnitDimensions[unit]},
    Apply[Times, Power @@@ uc] /. 
     s_String :> Symbol["CalculateUnits`UnitCommonSymbols`" <> s]
    ];
  ]

Basically all unit conversions are registered on a number of different symbols as DownValues. The core is CalculateUnits`UnitCommonSymbols`UnitLookup which stores the most important data, but for finding this symbol we need to register it on QuantityUnits`Private`$UnitReplacementRules and for some internal registration stuff we need to set a rule on CalculateUnits`UnitCommonSymbols`KnownUnit0Q.

Next we register the core conversion on CalculateUnits`UnitCommonSymbols`UnitLookup via the "FundamentalUnitValue" subproperty. This must reduce down to other more basic unit symbols (in general it reduces down to SI units it seems).

Some conversions are "fast-tracked" it appears, and for that we need to register the fast tracked multiplicative conversion on CalculateUnits`UnitTable`CUF0 and CalculateUnits`UnitTable`CUF1. Why we need both I don't know, but experimentally it seems I need it.

Finally, we set the "UnitDimensions" on UnitLookup as well for compatibility determination.

Now, we probably don't always want our rules to overwrite built-in stuff, so we'll define a function that makes sure that rules are cleared afterwards:

copyFromLink[link_, 
   expr_] :=
  ((*For fast copies of mutable expressions*)

   LinkWrite[link, expr];
   LinkRead[link]
   );
copyFromLink[link_][expr_] :=
  copyFromLink[link, expr];

withCustomConvertRules[
   expr_,
   rules : {
     (_ -> _Quantity) ..
     },
   useCache : True | False : False
   ] :=
  With[{ll = LinkCreate[LinkMode -> Loopback]},
   Internal`WithLocalSettings[
    None,
    Block[
     {
      QuantityUnits`Private`$CompatibleUnitQCache =
       If[useCache,
        copyFromLink[ll]@ 
         QuantityUnits`Private`$CompatibleUnitQCache,
        System`Utilities`HashTable[]
        ],
      QuantityUnits`Private`$KnownUnitQCache =
       If[useCache,
        copyFromLink[ll]@QuantityUnits`Private`$KnownUnitQCache,
        System`Utilities`HashTable[]
        ],
      QuantityUnits`Private`$UnitDimensionsCache =
       If[useCache,
        copyFromLink[ll]@QuantityUnits`Private`$UnitDimensionsCache,
        System`Utilities`HashTable[]
        ],
       QuantityUnits`Private`$UnitReplacementRules =
       QuantityUnits`Private`$UnitReplacementRules
      },
     Internal`InheritedBlock[
      {
       CalculateUnits`UnitCommonSymbols`UnitLookup,
       CalculateUnits`UnitCommonSymbols`KnownUnit0Q,
       CalculateUnits`UnitTable`CUF0,
       CalculateUnits`UnitTable`CUF1
       },
      addUnitConvertRule /@ rules;
      expr
      ]
     ],
    LinkClose[ll];
    ]
   ];
withCustomConvertRules~SetAttributes~HoldFirst

Note that I have a parameter in there that specified whether we use the cached values or not. Using the caches is faster (and to make sure they don't get mutated I add that copyFromLink function), but that can lead to unexpected surprises if there is a cached value that we're trying to overwrite, so by default I reinitialize them as new empty hash tables.

Now we can check if this works:

withCustomConvertRules[
 {
  UnitConvert[QuantityArray[RandomReal[{-1, 1}, {3, 2}], "Dogs"], 
   "Hartrees"],
  UnitConvert[Quantity[-76.1, "Hartrees"/"Seconds"], 
   "EnergyWavenumber"*"Megahertz"]
  },
 {
  "Dogs" -> Quantity[\[Pi], "Hartrees"],
  "EnergyWavenumber" -> 
   Quantity[
    1.98644585714892864895281898824816039`7.6813769430313075*^-23, 
    "Joules"]
  }
 ]

{StructuredArray[QuantityArray, {3, 2}, 
StructuredArray`StructuredData[
  QuantityArray, {{-2.687447536766794, 
   2.032843871270052}, {-2.5325772236678157`, -1.141664370726879}, {
   2.8417903071482153`, -0.966862479353852}}, "Hartrees", {{1}, {
   2}}]], Quantity[-16.702, "EnergyWavenumber" "Megahertz"]}

And it certainly appears as if it did. I may have missed some subtleties here, but this will get you started.

If you want your rules to stick around permanently, just use addUnitConvertRule

Extension

Here's an extension useful for all the chemists out there. First I'll use part of a dataset of conversions that I had laying around for some of my stuff.

$chemUnitsMap =
  <|"JoulesToWavenumbersMagnitude" :> 
    QuantityMagnitude[
     UnitConvert[
      Quantity[1, ("Joules")/("PlanckConstant" "SpeedOfLight")], 
      "Wavenumbers"]], 
   "JoulesToWavenumbersUnits" :> 
    Quantity[1, ("Joules")/("PlanckConstant" "SpeedOfLight")], 
   "JoulesToWavenumbers" :> 
    UnitConvert[
     Quantity[1, ("Joules")/("PlanckConstant" "SpeedOfLight")], 
     "Wavenumbers"], 
   "JoulesToMegahertzMagnitude" :> 
    QuantityMagnitude[
     UnitConvert[Quantity[1, ("Joules")/("PlanckConstant")], 
      "Megahertz"]], 
   "JoulesToMegahertzUnit" :> 
    Quantity[1/1000, ("Joules")/("PlanckConstant")], 
   "JoulesToMegahertz" :> 
    UnitConvert[Quantity[1, ("Joules")/("PlanckConstant")], 
     "Megahertz"], 
   "InertialConstant" :> 
    UnitConvert[
     Quantity[1/(8 \[Pi]^2), ("PlanckConstant")/(
      "AtomicMassUnit" ("Angstroms")^2)], "Megahertz"]|>;

This could be cleaner if I just used some conversion mapping, but this is good enough. Next I'll define a function that will extract these, but put the prefix "Chem" before each custom unit:

getChemUnitsMap[] :=
 DeleteDuplicatesBy[First]@
  KeyValueMap[
   If[StringMatchQ[#, __ ~~ "To" ~~ __ ~~ "Magnitude"],
     With[
      {
       bt =
        Flatten@
         StringCases[#, 
          base__ ~~ "To" ~~ targ__ ~~ "Magnitude" :> {base, targ}]
       },
      "Chem" <> bt[[2]] -> Quantity[1/#2, bt[[1]]]
      ],
     Nothing
     ] &,
   $chemUnitsMap
   ]

Then a function that extracts these "Chem" units, replaces the target units with their custom names, applies the conversion, then reverts this:

withChemUnits[expr_] :=
 With[{cm = getChemUnitsMap[]},
  withCustomConvertRules[
    ReleaseHold[
     Hold[expr] /.
      Thread[StringTrim[Keys[cm], "Chem"] -> Keys[cm]]
     ],
    cm
    ] // ReplaceAll[
     # /. Thread[Keys[cm] -> StringTrim[Keys[cm], "Chem"]],
     HoldPattern@StructuredArray[QuantityArray, dim_, data_] :>
      StructuredArray[QuantityArray, dim, 
       ReplacePart[data,
        3 -> (data[[3]] /. 
           Thread[Keys[cm] -> StringTrim[Keys[cm], "Chem"]])
        ]
       ]
     ] &
  ]
withChemUnits~SetAttributes~HoldFirst

Finally we make a little utility function and check that it works:

chemUnitConvert[a_, targ_] :=
 With[{res =
    Quiet[
     Check[
      UnitConvert[a, targ], 
      withChemUnits@UnitConvert[a, targ], 
      Quantity::compat
      ],
     Quantity::compat  
     ]
   },
  res /; Head[res] =!= UnitConvert
  ]
chemUnitConvert[a_List, targ_] :=
  chemUnitConvert[#, targ] & /@ a;

chemUnitConvert[
 QuantityArray[RandomReal[{}, {3, 2}], "Hartrees"], "Wavenumbers"]

StructuredArray[QuantityArray, {3, 2}, 
StructuredArray`StructuredData[
 QuantityArray, {{189634.50413450395`, 5975.770140545627}, {
  156145.9752769196, 193341.39616016933`}, {206581.95882060917`, 
  100386.34314265262`}}, "Wavenumbers", {{1}, {2}}]]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.