// Attempt to solve a system of equations.

use solvingTransformations.frink
showApproximations[false]
symbolicMode[true]

solveSystem[equations, solveForString] :=
{
   eqArray = new array
   size = length[equations]
   allUnknowns = new set
   for i=0 to size-1
   {
      syms = getSymbols[parseToExpression[equations@i]]
      equation = equations@i
//      println["Equation is $equation"]
//      println["Symbols are $syms"]
      eqArray.push[ [equation, syms] ]
//      println["eqArray is $eqArray"]
      allUnknowns = union[allUnknowns, syms]
//      println["allUnknowns is $allUnknowns"]
   }

//   println["All unknowns: $allUnknowns"]
   // Sort equations by number of unknowns (fewest first)
   sort[eqArray, {|a,b| length[a@1] <=> length[b@1]}]

   res = solveParts[eqArray, solveForString]
   res = eliminateDuplicates[res]
   res = eliminateOverconstrained[res]
//   res = eliminateSelfReferential[res]
   return res
}

// This is a recursive internal method to solve for other unknowns.
solveParts[eqArray, solveForString] :=
{
   results = new array
   size = length[eqArray]
   for i=0 to size-1
   {
      [eq, unknowns] = eqArray@i
      if unknowns.contains[solveForString]
      {
         oe = parseToExpression[eq]
         equation = parseToExpression["solve[$eq, $solveForString]"]
         solvedEq = transformExpression["solving", equation]
         results.push[solvedEq]
//         println["Result is $solvedEq"]
         otherEqs = eqArray.shallowCopy[]
         otherEqs.remove[i]

         otherSize = length[otherEqs]
         for unknown = unknowns
         {
            pUnknown = parseToExpression[unknown]
            if (unknown != solveForString)
               for j=0 to otherSize-1
               {
                  res2 = solveParts[otherEqs, unknown]
//                  println["Unknown is $unknown"]
//                  println["Res2 is $res2"]
                  for respart = res2
                  {
//                     println["Replacing in $solvedEq, $pUnknown becomes " + child[respart,1]]
                    // TODO:  Replace solving with simplification rules.
                     r = transformExpression["solving",substituteExpression[solvedEq, pUnknown, child[respart,1]]]
//                     println["  Result: $r"]
                     results.push[r]
                  }
               }
         }
      }
   }
   return flatten[results]
}


eliminateDuplicates[eqArray] :=
{
   i=0
   j=1
   while (i<length[eqArray])
   {
      ie = eqArray@i
      while (j<length[eqArray])
      {
         if structureEquals[ie, eqArray@j]
         {
//            println["Removing duplicate " + eqArray@j]
            eqArray.remove[j]   // Don't advance index in this case
         } else
         j=j+1
      }
      i=i+1
      j=i+1
   }
   return eqArray
}

// 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.
eliminateOverconstrained[eqArray] :=
{
   size = length[eqArray]
   unknowns = new array
   for i = 0 to size-1
      unknowns@i = getSymbols[child[eqArray@i,1]]

   res = new array

   isProper = true
   ILOOP:
   for i=0 to size-1
   {
      overconstrained = false
      j = 0
      do
      {
         overconstrained = isProperSubset[unknowns@j, unknowns@i]
         if (overconstrained)
            println[eqArray@j + " is a proper subset of " + eqArray@i]
         j=j+1
      } while (j < size) && ! overconstrained

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

   return res
}


symbolicMode[true]
//println[join["\n",solveSystem[["d === 2 r", "c === pi d", "a === pi r^2", "e===f", "g===h"], "r"]]]

// See https://web.archive.org/web/20210421140341/http://answers.yahoo.com/question/index?qid=20091120001614AAInec3
//
// Q. Given that P(A)=0.75, P(B|A)=0.8 and P(B|A')=0.6, calculate P(B) and P(A|B)
// Ans. 0.75, 0.8
// TODO:
//   eliminate pAprime and pAB, as those are what we want to solve for.
//println[join["\n",solveSystem[["pAB === pBA pA / pB", "pA === 1 - pAprime", "pAB === ( pBA pA ) / ( pBA pA + pBAprime pAprime ) "], "pB"]]]