/** This program tries to find a combination of basic constants that can
    produce the desired dimensions. */
use Matrix.frink
use formatEquation.frink

/** Finds a combination of base units that can be combined to produce a unit
    with the specified dimensions. */
planckify[unit, bases = ["c", "G", "hbar", "k", "electroncharge"]] :=
{
   [expMatrix, solMatrix, dims] = makeExponentMatrix[unit, bases]
   if exp == undef
      return undef
   
   println[formatTable[[[expMatrix.formatMatrix[], solMatrix.formatMatrix[]]]]]
   augment = expMatrix.augment[solMatrix] // Make augmented matrix
   augment.reduceRows[]   // Does the actual solve
   println[augment.formatMatrix[]]
   col = augment.getColumnAsArray[augment.getCols[]]  // Extract solution col
   str = ""
   for i = 0 to length[col]-1
   {
      ee = col@i
      if ee != 0
      {
         if isInteger[ee]
            if ee == 1
               expStr = ""
            else
               expStr = "^$ee"
         else
            expStr = "^(" + inputForm[ee] + ")"

         if str != ""
            str = str + " "
         
         str = str + bases@i + expStr
      }
   }

   println[str]
   println[]
   println[formatExpression[parseToExpression[str]]]
   println[]
   println[eval[str, false, true]]
}

makeExponentMatrix[unit, bases] :=
{
   dims = new set    // Set of base dimension names
   combos = new set  // Set of exponent combos

   for uu = bases
   {
      combo = dimensionsToArray[resolveUnit[uu]]
      if combos.contains[combo]
      {
         println["makeExponentMatrix:  overconstrained.  Exponent combo is repeated for $uu: $combo"]
         return undef
      } else
         combos.put[combo]
         
      for [dimname, exp] = combo
         dims.put[dimname]
   }

//   println["combos is $combos"]

   numBases = length[bases]
   numDims = length[dims]

   // Sort dimArray by the order that the exponents occur in standard order.
   dimArray = toArray[dims]
   sort[dimArray, {|a,b,data| data.indexOf[a] <=> data.indexOf[b]}, baseUnitNames[]]
//   println["dims is $dims"]
//   println["dimArray is $dimArray"]
//   println["Base units  " + baseUnitNames[]]

   // TODO: If unit contains a dimension not in dims, solution is impossible
   udims = dimensionsToArray[resolveUnit[unit]]
//   println["udims is $udims"]
   for [uname, exp] = udims
      if ! dims.contains[uname]
      {
         println["Unit $unit has base unit \"$uname\" which does not occur in the specified bases: $dimArray"]
         return undef
      }

   // Note that we transpose this matrix before returning it.   
   grid = new array[numBases]
   for uu = bases
   {
      resolved = resolveUnit[uu]
      row = new array[numDims]
      for dim = dimArray
         row.push[getExponent[resolved, dim]]

      grid.push[row]
   }

   // Note that we will transpose into a column this before returning it
   umatrix = new array
   ures = resolveUnit[unit]
   for dim = dimArray
      umatrix.push[getExponent[ures, dim]]

   return [new Matrix[grid.transpose[]],
           new Matrix[umatrix.transpose[]],
           dimArray]
}

/** Resolves a unit name to its value. */
resolveUnit[uu] :=
{
   if isUnit[uu]
      return uu
   
   if isString[uu]
      return unit[uu]

   println["resolveUnit:  Unhandled unit type " + type[uu]]
   return undef
}

println[planckify[psi]]