Download or view systemSolver2.frink in plain text format
// This program solves systems of equations.
// See "systemSolverTest.frink" or
// https://frinklang.org/fsp/solve2.fsp
// for examples of using it.
use allTransforms.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]
{
// If the unknowns in j are a proper subset of the unknowns in
// i, then it's probably a simpler form of the equation. The
// counterexample is with single solutions or degenerate cases
// like 0 where we want to see all solutions. For now, if the
// right-hand-side is zero, we don't consider it overconstrained
overconstrained = isProperSubset[unknowns@j, unknowns@i] && child[eqArray@j, 1] != 0 || (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["FunctionCall", ["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"
Download or view 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 19966 days, 21 hours, 59 minutes ago.