systemSolver2.frink

View or download systemSolver2.frink in plain text format


// This program solves systems of equations.
// See "systemSolverTest.frink" for examples of using it.

use derivatives.frink
use integrals.frink
use solvingTransformations.frink

showApproximations[false]

class System
{
   // An array of equations, parsed to expressions.  These should use the
   // === operator to separate the sides of the equations.
   var equations

   // A dictionary mapping symbol names (as strings)
   // to an array of the equation numbers
   // in which they occur.  (Or should it be the equation objects?)
   var variableOccurrences

   // A list of symbols in each equation.
   var symbols

   // A set of symbols that we're not going to solve for, because
   // they're intended to be units.  This is a set of strings.
   var ignoreSymbolsStrings = new set

   // A set of symbols that we're not going to solve for, because
   // they're intended to be units.  This is a set of symbols
   var ignoreSymbols = new set
   
   // A dictionary containing solved equations.
   // The key is the variable name (as a string)
   // The value is an array of equations that comprise solutions for that
   // variable.
   var solutions

   // Create a new System with the given equations, passed in as a list.
   new[equationList, ignoreSet=[], debug=false] :=
   {
      if debug
         println["orig ignoreSet is " + ignoreSet]
      for sym = ignoreSet
      {
         if debug
            println["Adding $sym"]
         addIgnoreSymbol[sym]
      }

      if debug
         println["ignoreSymbolsStrings is " + ignoreSymbolsStrings]
      
      equationTemp = deepCopy[equationList]
      equationTemp = solveSimultaneous[equationTemp, debug]
      //      equationTemp = solveSimultaneous[equationTemp, debug]

      equations = new array
      variableOccurrences = new dict
      symbols = new array

      for e = equationTemp
         addEquation[e]

      solutions = new dict
   }

   // Adds a symbol to the ignore list.  This adds it both to the
   // ignoreSymbols and ignoreSymbolsStrings sets.
   addIgnoreSymbol[sym is string] :=
   {
      ignoreSymbolsStrings.put[sym]
      ignoreSymbols.put[constructExpression["Symbol", sym]]
   }

   // Adds a new equation to the system.  Normally we would want all of the
   // equations to be known at construction time, but this may allow us to
   // add more equations to the system later, to provide a "refinement" of
   // the system we're using.  We generally don't want people to add equations
   // after we've partially solved the system, or if we do, we'll have to
   // clear any cache (like the "solutions" dictionary) and re-solve.
   addEquation[e] :=
   {
      i = length[equations]
      equations.push[e]
      
      syms = getSymbols[e]
      symbols@i = syms
      for sym = syms
      {
         if variableOccurrences@sym == undef
            variableOccurrences@sym = new array

         (variableOccurrences@sym).push[i]
      }
   }

   // Find the solutions for the symbol x (specified as a string.)
   // This function sets up for a recursive call.
   solveFor[x, debug=false] :=
   {
      prev = solutions@x
      if prev
         return prev
      
      results = new array       // Array of expressions containing results.
      usedEquations = new set   // Set of numbers of used equations.
      usedSymbols = new set     // Set of names of used symbols.
      results = solveFor[x, results, usedEquations, usedSymbols, debug]
      
      solutions@x = flatten[results]     // Put into dictionary
      return results
   }

   // The private recursive call.
   solveFor[x, results, usedEquations, usedSymbols, debug=false] :=
   {
      usedSymbols.put[x]
      // println["Used symbol $x"]
      for ei = variableOccurrences@x  // For each equation that x is used in...
      {
         if usedEquations.contains[ei]
            next                // Skip equation if it's already been used.
            
         usedEquations.put[ei]
         // println["Used equation $ei"]
         eq = equations@ei

         results.push[solveSingle[eq, x]]
         // println["The equation is $solveEq"]
         // println[result]
      }

      // println["Results before recursion:\n" + results]

      // We count the number of results currently on the list and only
      // substitute into the existing results, as we'll be adding to the
      // results list.
      s = length[results]
      for n = 0 to s-1
      {
         res = results@n
         symbols = array[getSymbols[res]]

         for sym = symbols
            if ! usedSymbols.contains[sym]
            {
               results2 = new array
//               println["Solving for $sym"]
               results2 = solveFor[sym, results2, deepCopy[usedEquations], deepCopy[usedSymbols], debug]
//               println["Solved for $sym, got back " + results2]
               symSymbol = constructExpression["Symbol", sym]

               // TODO:  Only substitute below if the result is a fully solved
               // === expression with the appropriate form.
               for eq2 = flatten[results2]
               {
                  if structureEquals[_a === _b, eq2]
                  {
                     results.push[transformExpression[substituteExpression[res, symSymbol, child[eq2, 1]]]]
//                   println["Used $eq2"]
                  } else
                       if debug
                           println["Discarded unsolved $eq2"]
               }
            }
      }

      return eliminateOverconstrained[results, debug]
   }

   // Find the solutions for the symbol x (specified as a string.)
   // values is an array of [varname, value] pairs.
   solveForValues[x, values, allValues=false, debug=false] :=
   {
      origSolutions = solveFor[x, debug]
      newSolutions = new array
//      if debug
//         println["solveForValues: origSolutions are:" + origSolutions]
      
      for sol = origSolutions
      {
         newSol = sol
         for [symName, value] = values
         {
            if debug
            {
               println["\nReplacing $symName with $value"]
               println["old is:" + newSol]
            }
            
            newSol = replaceVar[newSol, symName, value]
            
            if debug
               println["new is:" + newSol]
         }

         newSolutions.push[transformExpression[newSol]]
      }

      if (allValues == false)
         newSolutions = eliminateOverconstrained[newSolutions, debug]

      solvedVals = new array
      allVals = new array
      for sol = newSolutions
      {
         right = getChild[sol,1]
         allVals.push[right]
         if length[setDifference[getSymbols[right], ignoreSymbolsStrings]] == 0 // Solved?
            solvedVals.push[right]
      }

      if (allValues==false and length[solvedVals] > 0)  // Have fully-solved?
         return solvedVals
      else
         return allVals
   }

   // Find the solutions for the symbol x (specified as a string.)
   // values is an array of [varname, value] pairs.
   // This just sets up arguments for the recursive call.
   solveForValuesOld[x, values, allValues=false, debug=false] :=
   {
      // Create a new set of values we've already substituted for.
      substituted = new set
      // Dictionary<set, equations>
      cachedSolutions = new dict
      return solveForValuesRecursive[x, values, substituted, cachedSolutions, allValues, debug]
   }

   // Puts the specified values into the cache.
   class putCache[cachedSolutions, varName, substituted, vals] :=
   {
      if ! cachedSolutions.containsKey[varName]
         cachedSolutions@varName = new dict

      cachedSolutions@varName@substituted = vals
   }

   // This is the "private" method that does recursive solution of
   // solutions for the symbol x.
   solveForValuesRecursive[x, values, substituted, cachedSolutions, allValues, debug=false] :=
   {
      if cachedSolutions.containsKey[x]
         if (cachedSolutions@x).containsKey[substituted]
            return cachedSolutions@x@substituted

      if debug
         println["Solving for $x, substituted is $substituted"]

      sols = solveFor[x, debug]
      valset = new set
      for [varname, val] = values
      {
         if varname == x        // User passed in solution?
         {
            putCache[cachedSolutions, x, substituted, [val]]
            return [val]          // TODO:  Warn as overspecified?
         }
//         sols = replaceVar[sols, varname, val]  // THINK ABOUT:  Do this later?
         valset.put[varname]
      }

      if debug
         println["Initial solutions for $x are $sols"]
      if (allValues == false)
         sols = eliminateOverconstrained[flatten[sols]]
      sort[sols, {|a,b| length[getSymbols[a]] <=> length[getSymbols[b]]}]

      // Ignore the variable we're solving for and things in the ignore list.
      ignoreVals = new set
      ignoreVals.put[x]
      ignoreVals = union[ignoreVals, ignoreSymbolsStrings]
      res = new array
      for sol = sols
         if length[setDifference[getSymbols[sol], ignoreVals]] == 0  // Solved?
            for s = array[sol]
               res.push[child[s,1]]  // Eval here?

      if length[res] > 0        // Got fully-solved solutions?
      {
         putCache[cachedSolutions, x, substituted, res]
         if debug
            println["Solutions for $x (fully solved) are $res"]
         return res
      } else    // Not fully solved, return symbolic solutions 
      {
          xSym = constructExpression["Symbol", x]
          allIgnores = deepCopy[ignoreSymbolsStrings]
          allIgnores = union[allIgnores, substituted]
          allIgnores = union[allIgnores, valset]
          allIgnores.put[x]

          // Get all symbols in all solutions, minus the things we're
          // ignoring.
          allSyms = setDifference[getSymbols[sols], allIgnores]

          // solsDict is a dictionary <solveFor, [equations]>
          solsDict = new dict
          for u = allSyms
          {
             newsub = deepCopy[substituted]
             newsub.put[x]
             solsDict@u = solveForValuesRecursive[u, values, newsub, cachedSolutions, allValues, debug]
          }

//          println["solsDict is $solsDict"]

          // Now iterate through each solution and substitute in each 
          for s = sols
             if structureEquals[_a === _b, s]
             {
                right = child[s,1]
                unknowns = setDifference[getSymbols[right], allIgnores]
                newsub = deepCopy[substituted]
                newsub.put[x]

                for u = unknowns
                {
                   res2 = solveForValuesRecursive[u, values, newsub, cachedSolutions, allValues, debug]
                   uSym = constructExpression["Symbol", u]
                   for r2 = flatten[res2]
                   {
                      if debug
                         println["Substituting in $s, from $u, to $r2"]
                      s = substituteExpression[s, uSym, r2]
                   }
                }
                
                s = transformExpression[s]
                
                g = getSymbols[child[s, 1]]
                if ! g.contains[x]  // Solution contains x?
                {
                   res.push[s]
                   if debug
                      println["Adding solution $s"]
                } else
                   if debug
                      println["Solution $s rejected because it contains $x"]
             } else
                if debug
                   println["Solution $s rejected because unsolved."]
             

          for [varname, val] = values
             res = replaceVar[res, varname, val]  // THINK ABOUT:  Do this later?

          if (allValues == false)
             res = eliminateOverconstrained[res, debug]
          size = length[res]
          for i=0 to size-1
          {
             child = child[res@i, 1]  // Return only right-hand-sides
             if isSubset[getSymbols[child], ignoreSymbolsStrings]
                res@i = child  // Eval here?
             else
                res@i = child
          }
          if debug
             println["Solutions for $x are $res"]

          putCache[cachedSolutions, x, substituted, res]
          return res  //, false, debug]
      }
   }
   
   // Solve for all unknowns in the system.
   // Setting fullReplace=true makes this work harder at performing full
   // substitutions of all variables in all solutions, but it may be
   // prohibitively slow for a system with a large number of solutions and
   // variables.
   solveAll[fullReplace=false, debug=false] :=
   {
      results = new array
      for [sym, appear] = variableOccurrences
         if ! ignoreSymbolsStrings.contains[sym]
            results.push[solveFor[sym, debug]]

      if debug
         println["Before fullReplace:\n" + join["\n", flatten[results]]]

         // TODO:  Comment back in fullReplace.
         //println["Results are" + results]

      if fullReplace
         return fullReplace[flatten[results], true, debug]
      else
         return flatten[results]
   }

   
   // Once all equations have been solved, this function will redistribute all
   // solutions through all equations.
   fullReplace[equations, allValues=false, debug=false] :=
   {
      // This is a dictionary with variable names as keys and a list of 
      // equations that are the solutions with those variables.  It will
      // replace the solutions member when we're complete.
      newSolutions = new dict
      solArray = new array

      for eq = equations
      {
         if debug
            println["Equation is: $eq"]
         left = getChild[eq, 0]
         right = getChild[eq, 1]
         leftStr = toString[left]
//         println["leftStr is $leftStr, type of leftStr is " + type[leftStr]]
         allSyms = getSymbols[right]
//         println["Symbols on right are " + allSyms]
//         println["ignoreSymbols are " + ignoreSymbols]
//         println["ignoreSymbolsStrings are " + ignoreSymbolsStrings]
         replaceSyms = setDifference[allSyms, ignoreSymbolsStrings]

         symList = toArray[replaceSyms]
         if debug
            println["symList is " + symList]
         stateList = new array
         solList = new array

         for sym = symList
         {
            stateList.push[new range[false, true]]
            
            // Make an array of right-hand-sides of solutions for that symbol
            solList2 = new array
            for sol = solutions@sym
            {
               solRight = child[sol, 1]
//               println["$sym, " + getSymbols[solRight]]
               if ! (getSymbols[solRight].contains[leftStr])
                  solList2.push[solRight]
//               else
//                  println["Rejected $sol because it contains $leftStr"]
            }
            solList.push[solList2]
         }

         numStates = length[symList]

         newEq = eq
         if debug
            println["Before replace: $newEq"]

         multifor states = stateList
         {
            tempSolArray = [eq]
            if debug
            {
               println["states is $states"]
               println["solList is $solList"]
            }

            newSolArray = new array
            for i = 0 to numStates-1
            {
               replaced = false
               if states@i == true
               {
                  for sol = solList@i
                  {
                     replaced = true
                     existingSols = length[tempSolArray]
                     
                     for k = 0 to existingSols-1
                     {
                        curreq = tempSolArray@k
                        curreq = substituteExpression[curreq, constructExpression["Symbol", symList@i], sol]
                        newSolArray.push[transformExpression[curreq]]
                     }
                  }
               }

               if replaced
                  tempSolArray = newSolArray
            }

            tempSolArray = eliminateOverconstrained[tempSolArray]
            if debug
               println["after replace:  $tempSolArray\n"]
            // TODO:  Change to set to avoid dups?
            solArray.pushAll[tempSolArray]
         }
      }
      
      println["Pre-Reduced is $solArray"]
      solArray = eliminateOverconstrained[solArray]
      
      if debug
         println["Reduced is $solArray"]

      for sol = solArray
      {
         key = toString[child[sol,0]]
         if newSolutions.containsKey[key]
            (newSolutions@key).push[sol]
         else
            newSolutions@key = [sol]
      }

      solutions = newSolutions
      
      return solArray
   }

   // Once all equations have been solved, this function will redistribute all
   // solutions through all equations.
   fullReplaceOld[equations, debug=false] :=
   {
      newSolutions = new dict
      res = deepCopy[equations]
      for eq = equations
      {
         if debug
            println["Working on equation $eq"]

         eqLeft = child[eq,0]
         eqRight = child[eq,1]
         lsymbols = getSymbols[eqLeft]
         lsymbolStr = array[lsymbols]@0  // Kludge to get symbol on left
         rsymbols = getSymbols[eqRight]
         syms = getSymbols[eqRight]  // Symbols on right-hand-side of equation
         syms = setDifference[syms, ignoreSymbolsStrings]
         symArray = array[syms]
         len = length[symArray]
         if debug
            println[" symArray is $symArray"]

         if len > 0
         {
            // Now go through and replace all permutations of variables.
            // This essentially counts through all the 2^n binary permutations
            states = new array 
            for i=0 to len-2
               states@i = false

            states@(len-1) = true

            i = len-1
            r = eq
            while i>=0
            {
               if (debug)
                  println["states is $states, $symArray, j is $j, symArray@j is " + symArray@j ]
               for repval = flatten[solveFor[symArray@j]]
               {
                  replaceVal = child[repval,1]
                  for j=0 to len-1
                     if states@j
                        if ! (getSymbols[replaceVal].contains[lsymbolStr])
                        {
                           if debug
                              println["    Replacing " + symArray@j + " with " + replaceVal + " from $repval in " + r]
                           r = transformExpression[replaceVar[r, symArray@j, replaceVal]]
                           if debug
                              println["     Solution is $r"]
                        }
               }

               // Eliminate x === x
               if ! structureEquals[eqLeft, child[r,1]]
                  res.push[r]

               if structureEquals[r, eq]
                  println["Expression $r unchanged."]
               else
                  println["Expression $eq changed to $r"]
  
               // Advance to next binary state
               flipped = false
               i = len-1
               while i>=0 and !flipped
               {
                  // Enter next state
                  if states@i == false
                  {
                     states@i = true
                     flipped = true
                  } else
                  {  // Carry
                     states@i = false
                     i = i - 1
                  }
               }
            }
         }
      }

      res = eliminateOverconstrained[res]
      println["\nAfter fullReplace:\n" + join["\n", res] + "\n"]

      return res
   }

   // Attempts to solve simultaneous equations in the system.
   solveSimultaneous[equationList, debug=true] :=
   {
      // Sort so that shortest equations are first.  This simplifies
      // substitutions and makes resulting equations shorter.
      // We can't do the latter because the sort function doesn't have
      // access to ignoreSymbols!  We have to create a closure to make this
      // happen?
      //      sort[equationList, {|a,b| length[setDifference[getSymbols[a],
      //      ignoreSymbols]] <=> length[setDifference[getSymbols[b],
      //      ignoreSymbols]]}]
      sort[equationList, {|a,b| length[getSymbols[a]] <=> length[getSymbols[b]]}]
      
      s = length[equationList]
      for i = 0 to s-2
      {
         ei = equationList@i
         ui = getSymbols[ei]
         aui = array[ui]
         lui = length[aui]
         
         JLOOP:
         for j = i+1 to s-1
         {
            ej = equationList@j
            uj = getSymbols[ej]

            commonSet = intersection[ui, uj]
            if length[commonSet] >= 1 or lui == 1
            {
               commonMinusIgnore = setDifference[commonSet, ignoreSymbols]
               if length[commonMinusIgnore] > 0
                  solveFor = sort[array[commonMinusIgnore]]@0
               else
                  solveFor = array[ui]@0
               
               eqSolved = solveSingle[ei, solveFor]

               if (debug)
                  println["Solving $ei for $solveFor, result is $eqSolved"]

               for part = flatten[eqSolved]
               {
                  // See if part is fully solved.
                  if structureEquals[_a === _b, part]
                  {
                     // Right-hand-side of solved equation
                     eqSolution = child[part,1]  

                     xSymbol = constructExpression["Symbol", solveFor]

                     // Substitute result into other equations.
                     for k=0 to s-1
                        if k != i     // Don't substitute into this equation.
                        {
                           orig = equationList@k
                           subst = substituteExpression[equationList@k, xSymbol, eqSolution]
                           if (orig != subst)
                           {
                              subst = transformExpression[subst]
                              // Only replace the equation if we actually simplified it.
                              if length[getSymbols[subst]] <= length[getSymbols[orig]]
                              {
                                 if debug
                                    println["Useful substitution: $part TO $subst"]
                                 equationList@k = subst
                              } else
                              if debug
                                 println["Unuseful substitution: $part TO $subst"]
                              }
                           }
                           break JLOOP
                  } else
                      println["Unsolved: $part"]
               }
            }
         }
      }
      return equationList
   }

   // This function eliminates overconstrained equations.  For example, a
   // system containing the solutions a===1/2 c r  and  a===c d^-1 r^2  is
   // overconstrained because a value can always be obtained with the first
   // equation.  The second is not necessary, and could lead to
   // inconsistent results.   This method ignores any symbols listed in the
   // ignoreSymbols list, (these are probably units,) eliminating them from
   // the equations.
   eliminateOverconstrained[eqArray, debug=false] :=
   {
      size = length[eqArray]
      unknowns = new array
      lefts = new array
      for i = 0 to size-1
      {
         lefts@i = child[eqArray@i, 0]
         unknowns@i = setDifference[getSymbols[child[eqArray@i,1]], ignoreSymbols]
      }

      res = new array

      // Check for duplicates.
      for i=0 to size-1
      {
         overconstrained = false
         j = 0
         do
         {
            if i != j and structureEquals[lefts@i, lefts@j]
            {
               overconstrained = isProperSubset[unknowns@j, unknowns@i] || (i<j and structureEquals[eqArray@i, eqArray@j])
               if (overconstrained)
                  if debug
                     println[eqArray@j + " is a proper subset or match of " + eqArray@i]
            }
            j=j+1
         } while (j < size) and ! overconstrained

         if (! overconstrained)
            res.push[eqArray@i]  // If we got here, no j is a proper subset of i.
      }

      return res
   }
}

// Solve a single equation for the specified symbol.
solveSingle[eq, symbol] :=
{
   xSymbol = constructExpression["Symbol", symbol]
   // We could use this in symbolic mode, otherwise it warns.
   // solveEq = solve[eq, xSymbol]
   solveEq = constructExpression["Function", ["solve", eq, xSymbol]]
   
   return transformExpression[solveEq]
}

// Replace a variable in the specified equation with the specified value.
// You should probably wrap the "eq" and "value" in a noEval[] block if
// you're not passing them in from a variable.
replaceVar[eq, varName, value] :=
{
   sym = constructExpression["Symbol", varName]
   res = substituteExpression[eq, sym, value]
   return res
}

/** This filters a list, returning only positive results. */
positive[x] :=
{
   result = new array
   for sol = x
   {
      rs = realSignum[sol]
      if rs != undef && rs > 0
         result.push[sol]
   }
   return result      
}

"systemSolver2.frink included OK"


View or download systemSolver2.frink in plain text format


This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 17592 days, 15 hours, 22 minutes ago.