parameters: (self,dialyticeqs,ileftvar,returnmatrix=False,getsubs=None)
return_statement: return exportcoeffeqs,origmonoms
docstring: Return the coefficients to solve equations dialytically (Salmon 1885) leaving out variable index ileftvar.
  
          Extract the coefficients of 1, leftvar**1, leftvar**2, ... of every equation
          every len(dialyticeqs)*len(monoms) coefficients specify one degree of all the equations (order of monoms is specified in exportmonomorder
          there should be len(dialyticeqs)*len(monoms)*maxdegree coefficients
  
          Method also checks if the equations are linearly dependent
docstring_summary: Return the coefficients to solve equations dialytically (Salmon 1885) leaving out variable index ileftvar.
docstring_tokens: Return the coefficients to solve equations dialytically ( Salmon 1885 ) leaving out variable index ileftvar .
function: def solveDialytically(self,dialyticeqs,ileftvar,returnmatrix=False,getsubs=None):
          """ Return the coefficients to solve equations dialytically (Salmon 1885) leaving out variable index ileftvar.
  
          Extract the coefficients of 1, leftvar**1, leftvar**2, ... of every equation
          every len(dialyticeqs)*len(monoms) coefficients specify one degree of all the equations (order of monoms is specified in exportmonomorder
          there should be len(dialyticeqs)*len(monoms)*maxdegree coefficients
  
          Method also checks if the equations are linearly dependent
          """
          self._CheckPreemptFn(progress=0.12)
          if len(dialyticeqs) == 0:
              raise self.CannotSolveError('solveDialytically given zero equations')
          
          allmonoms = set()
          origmonoms = set()
          maxdegree = 0
          leftvar = dialyticeqs[0].gens[ileftvar]
          extradialyticeqs = []
          for peq in dialyticeqs:
              if sum(peq.degree_list()) == 0:
                  log.warn('solveDialytically: polynomial %s degree is 0',peq)
                  continue
              for m in peq.monoms():
                  mlist = list(m)
                  maxdegree=max(maxdegree,mlist.pop(ileftvar))
                  allmonoms.add(tuple(mlist))
                  origmonoms.add(tuple(mlist))
                  mlist[0] += 1
                  allmonoms.add(tuple(mlist))
              
              # check if any monoms are not expressed in this poly, and if so, add another poly with the monom multiplied, will this give bad solutions?
              for igen in range(len(peq.gens)):
                  if all([m[igen]==0 for m in peq.monoms()]):
                      log.debug('adding extra equation multiplied by %s', peq.gens[igen])
                      extradialyticeqs.append(peq*peq.gens[igen])
                      # multiply by peq.gens[igen]
                      for m in peq.monoms():
                          mlist = list(m)
                          mlist[igen] += 1
                          maxdegree=max(maxdegree,mlist.pop(ileftvar))
                          allmonoms.add(tuple(mlist))
                          origmonoms.add(tuple(mlist))
                          mlist[0] += 1
                          allmonoms.add(tuple(mlist))
          
          dialyticeqs = list(dialyticeqs) + extradialyticeqs # dialyticeqs could be a tuple
          allmonoms = list(allmonoms)
          allmonoms.sort()
          origmonoms = list(origmonoms)
          origmonoms.sort()
  
          if len(origmonoms) == 0 or len(allmonoms) == 0:
              raise self.CannotSolveError('solveDialytically has no equations')
          
          if len(allmonoms)<2*len(dialyticeqs):
              log.warn('solveDialytically equations %d > %d, should be equal...', 2*len(dialyticeqs),len(allmonoms))
              # TODO not sure how to select the equations
              N = len(allmonoms)//2
              dialyticeqs = dialyticeqs[:N]
          if len(allmonoms) == 0 or len(allmonoms)>2*len(dialyticeqs):
              raise self.CannotSolveError('solveDialytically: more unknowns than equations %d>%d'%(len(allmonoms), 2*len(dialyticeqs)))
          
          Mall = [zeros((2*len(dialyticeqs),len(allmonoms))) for i in range(maxdegree+1)]
          Mallindices = [-ones((2*len(dialyticeqs),len(allmonoms))) for i in range(maxdegree+1)]
          exportcoeffeqs = [S.Zero]*(len(dialyticeqs)*len(origmonoms)*(maxdegree+1))
          for ipeq,peq in enumerate(dialyticeqs):
              for m,c in peq.terms():
                  mlist = list(m)
                  degree=mlist.pop(ileftvar)
                  exportindex = degree*len(origmonoms)*len(dialyticeqs) + len(origmonoms)*ipeq+origmonoms.index(tuple(mlist))
                  assert(exportcoeffeqs[exportindex] == S.Zero)
                  exportcoeffeqs[exportindex] = c
                  Mall[degree][len(dialyticeqs)+ipeq,allmonoms.index(tuple(mlist))] = c
                  Mallindices[degree][len(dialyticeqs)+ipeq,allmonoms.index(tuple(mlist))] = exportindex
                  mlist[0] += 1
                  Mall[degree][ipeq,allmonoms.index(tuple(mlist))] = c
                  Mallindices[degree][ipeq,allmonoms.index(tuple(mlist))] = exportindex
  
              # check if any monoms are not expressed in this poly, and if so, add another poly with the monom multiplied, will this give bad solutions?
              for igen in range(len(peq.gens)):
                  if all([m[igen]==0 for m in peq.monoms()]):
                      for m,c in peq.terms():
                          mlist = list(m)
                          mlist[igen] += 1
                          degree=mlist.pop(ileftvar)
  
                          if tuple(mlist) not in origmonoms:
                              raise self.CannotSolveError('equations too simple, monom %r is not in %r'%(mlist, origmonoms))
                          
                          exportindex = degree*len(origmonoms)*len(dialyticeqs) + len(origmonoms)*ipeq+origmonoms.index(tuple(mlist))
                          assert(exportcoeffeqs[exportindex] == S.Zero)
                          exportcoeffeqs[exportindex] = c
                          Mall[degree][len(dialyticeqs)+ipeq,allmonoms.index(tuple(mlist))] = c
                          Mallindices[degree][len(dialyticeqs)+ipeq,allmonoms.index(tuple(mlist))] = exportindex
                          mlist[0] += 1
                          Mall[degree][ipeq,allmonoms.index(tuple(mlist))] = c
                          Mallindices[degree][ipeq,allmonoms.index(tuple(mlist))] = exportindex
  
          # have to check that the determinant is not zero for several values of ileftvar! It is very common that
          # some equations are linearly dependent and not solvable through this method.
          if self.testconsistentvalues is not None:
              linearlyindependent = False
              for itest,subs in enumerate(self.testconsistentvalues):
                  if getsubs is not None:
                      # have to explicitly evaluate since testsubs can be very complex
                      subsvals = [(s,v.evalf()) for s,v in subs]
                      try:
                          subs = subsvals+getsubs(subsvals)
                      except self.CannotSolveError as e:
                          # getsubs failed (sometimes it requires solving inverse matrix), so go to next set
                          continue
                  # have to sub at least twice with the global symbols
                  A = Mall[maxdegree].subs(subs)
                  for i in range(A.shape[0]):
                      for j in range(A.shape[1]):
                          A[i,j] = self._SubstituteGlobalSymbols(A[i,j]).subs(subs).evalf()
                  eps = 10**-(self.precision-3)
                  try:
                      Anumpy = numpy.array(numpy.array(A), numpy.float64)
                  except ValueError as e:
                      log.warn(u'could not convert to numpy array: %s', e)
                      continue
                  
                  if numpy.isnan(numpy.sum(Anumpy)):
                      log.info('A has NaNs')
                      break
                  eigenvals = numpy.linalg.eigvals(Anumpy)
                  if all([Abs(f) > eps for f in eigenvals]):
                      try:
                          Ainv = A.inv(method='LU')
                      except ValueError as e:
                          log.error('error when taking inverse: %s', e)
                          continue
                      B = Ainv*Mall[1].subs(subs)
                      for i in range(B.shape[0]):
                          for j in range(B.shape[1]):
                              B[i,j] = self._SubstituteGlobalSymbols(B[i,j]).subs(subs).evalf()
                      C = Ainv*Mall[0].subs(subs).evalf()
                      for i in range(C.shape[0]):
                          for j in range(C.shape[1]):
                              C[i,j] = self._SubstituteGlobalSymbols(C[i,j]).subs(subs).evalf()
                      A2 = zeros((B.shape[0],B.shape[0]*2))
                      for i in range(B.shape[0]):
                          A2[i,B.shape[0]+i] = S.One
                      A2=A2.col_join((-C).row_join(-B))
                      eigenvals2,eigenvecs2 = numpy.linalg.eig(numpy.array(numpy.array(A2),numpy.float64))
                      # check if solutions can actually be extracted
                      # find all the zero eigenvalues
                      roots = []
                      numrepeating = 0
                      for ieig,eigenvalue in enumerate(eigenvals2):
                          if abs(numpy.imag(eigenvalue)) < 1e-12:
                              if abs(numpy.real(eigenvalue)) > 1:
                                  ev = eigenvecs2[A.shape[0]:,ieig]
                              else:
                                  ev = eigenvecs2[:A.shape[0],ieig]
                              if abs(ev[0]) < 1e-14:
                                  continue
                              br = ev[1:] / ev[0]
                              dists = abs(numpy.array(roots) - numpy.real(eigenvalue))
                              if any(dists<1e-7):
                                  numrepeating += 1
                              roots.append(numpy.real(eigenvalue))
                      if numrepeating > 0:
                          log.info('found %d repeating roots in solveDialytically matrix: %s',numrepeating,roots)
                          # should go on even if there's repeating roots?
                          continue
                      Atotal = None
                      for idegree in range(maxdegree+1):
                          Adegree = Mall[idegree].subs(subs)
                          for i in range(Adegree.shape[0]):
                              for j in range(Adegree.shape[1]):
                                  Adegree[i,j] = self._SubstituteGlobalSymbols(Adegree[i,j]).subs(subs).evalf()
                          if Atotal is None:
                              Atotal = Adegree
                          else:
                              Atotal += Adegree*leftvar**idegree
                      # make sure the determinant of Atotal is not-zero for at least several solutions
                      leftvarvalue = leftvar.subs(subs).evalf()
                      hasnonzerodet = False
                      for testvalue in [-10*S.One, -S.One,-0.5*S.One, 0.5*S.One, S.One, 10*S.One]:
                          Atotal2 = Atotal.subs(leftvar,leftvarvalue+testvalue).evalf()
                          detvalue = Atotal2.det()
                          if abs(detvalue) > 1e-10:
                              hasnonzerodet = True
                      if not hasnonzerodet:
                          log.warn('has zero det, so failed')
                      else:
                          linearlyindependent = True
                      break
                  else:
                      log.info('not all abs(eigenvalues) > %e. min is %e', eps, min([Abs(f) for f in eigenvals if Abs(f) < eps]))
              if not linearlyindependent:
                  raise self.CannotSolveError('equations are not linearly independent')
  
          if returnmatrix:
              return Mall,allmonoms
  
          return exportcoeffeqs,origmonoms
function_tokens: def solveDialytically ( self , dialyticeqs , ileftvar , returnmatrix = False , getsubs = None ) : self . _CheckPreemptFn ( progress = 0.12 ) if len ( dialyticeqs ) == 0 : raise self . CannotSolveError ( 'solveDialytically given zero equations' ) allmonoms = set ( ) origmonoms = set ( ) maxdegree = 0 leftvar = dialyticeqs [ 0 ] . gens [ ileftvar ] extradialyticeqs = [ ] for peq in dialyticeqs : if sum ( peq . degree_list ( ) ) == 0 : log . warn ( 'solveDialytically: polynomial %s degree is 0' , peq ) continue for m in peq . monoms ( ) : mlist = list ( m ) maxdegree = max ( maxdegree , mlist . pop ( ileftvar ) ) allmonoms . add ( tuple ( mlist ) ) origmonoms . add ( tuple ( mlist ) ) mlist [ 0 ] += 1 allmonoms . add ( tuple ( mlist ) ) # check if any monoms are not expressed in this poly, and if so, add another poly with the monom multiplied, will this give bad solutions? for igen in range ( len ( peq . gens ) ) : if all ( [ m [ igen ] == 0 for m in peq . monoms ( ) ] ) : log . debug ( 'adding extra equation multiplied by %s' , peq . gens [ igen ] ) extradialyticeqs . append ( peq * peq . gens [ igen ] ) # multiply by peq.gens[igen] for m in peq . monoms ( ) : mlist = list ( m ) mlist [ igen ] += 1 maxdegree = max ( maxdegree , mlist . pop ( ileftvar ) ) allmonoms . add ( tuple ( mlist ) ) origmonoms . add ( tuple ( mlist ) ) mlist [ 0 ] += 1 allmonoms . add ( tuple ( mlist ) ) dialyticeqs = list ( dialyticeqs ) + extradialyticeqs # dialyticeqs could be a tuple allmonoms = list ( allmonoms ) allmonoms . sort ( ) origmonoms = list ( origmonoms ) origmonoms . sort ( ) if len ( origmonoms ) == 0 or len ( allmonoms ) == 0 : raise self . CannotSolveError ( 'solveDialytically has no equations' ) if len ( allmonoms ) < 2 * len ( dialyticeqs ) : log . warn ( 'solveDialytically equations %d > %d, should be equal...' , 2 * len ( dialyticeqs ) , len ( allmonoms ) ) # TODO not sure how to select the equations N = len ( allmonoms ) // 2 dialyticeqs = dialyticeqs [ : N ] if len ( allmonoms ) == 0 or len ( allmonoms ) > 2 * len ( dialyticeqs ) : raise self . CannotSolveError ( 'solveDialytically: more unknowns than equations %d>%d' % ( len ( allmonoms ) , 2 * len ( dialyticeqs ) ) ) Mall = [ zeros ( ( 2 * len ( dialyticeqs ) , len ( allmonoms ) ) ) for i in range ( maxdegree + 1 ) ] Mallindices = [ - ones ( ( 2 * len ( dialyticeqs ) , len ( allmonoms ) ) ) for i in range ( maxdegree + 1 ) ] exportcoeffeqs = [ S . Zero ] * ( len ( dialyticeqs ) * len ( origmonoms ) * ( maxdegree + 1 ) ) for ipeq , peq in enumerate ( dialyticeqs ) : for m , c in peq . terms ( ) : mlist = list ( m ) degree = mlist . pop ( ileftvar ) exportindex = degree * len ( origmonoms ) * len ( dialyticeqs ) + len ( origmonoms ) * ipeq + origmonoms . index ( tuple ( mlist ) ) assert ( exportcoeffeqs [ exportindex ] == S . Zero ) exportcoeffeqs [ exportindex ] = c Mall [ degree ] [ len ( dialyticeqs ) + ipeq , allmonoms . index ( tuple ( mlist ) ) ] = c Mallindices [ degree ] [ len ( dialyticeqs ) + ipeq , allmonoms . index ( tuple ( mlist ) ) ] = exportindex mlist [ 0 ] += 1 Mall [ degree ] [ ipeq , allmonoms . index ( tuple ( mlist ) ) ] = c Mallindices [ degree ] [ ipeq , allmonoms . index ( tuple ( mlist ) ) ] = exportindex # check if any monoms are not expressed in this poly, and if so, add another poly with the monom multiplied, will this give bad solutions? for igen in range ( len ( peq . gens ) ) : if all ( [ m [ igen ] == 0 for m in peq . monoms ( ) ] ) : for m , c in peq . terms ( ) : mlist = list ( m ) mlist [ igen ] += 1 degree = mlist . pop ( ileftvar ) if tuple ( mlist ) not in origmonoms : raise self . CannotSolveError ( 'equations too simple, monom %r is not in %r' % ( mlist , origmonoms ) ) exportindex = degree * len ( origmonoms ) * len ( dialyticeqs ) + len ( origmonoms ) * ipeq + origmonoms . index ( tuple ( mlist ) ) assert ( exportcoeffeqs [ exportindex ] == S . Zero ) exportcoeffeqs [ exportindex ] = c Mall [ degree ] [ len ( dialyticeqs ) + ipeq , allmonoms . index ( tuple ( mlist ) ) ] = c Mallindices [ degree ] [ len ( dialyticeqs ) + ipeq , allmonoms . index ( tuple ( mlist ) ) ] = exportindex mlist [ 0 ] += 1 Mall [ degree ] [ ipeq , allmonoms . index ( tuple ( mlist ) ) ] = c Mallindices [ degree ] [ ipeq , allmonoms . index ( tuple ( mlist ) ) ] = exportindex # have to check that the determinant is not zero for several values of ileftvar! It is very common that # some equations are linearly dependent and not solvable through this method. if self . testconsistentvalues is not None : linearlyindependent = False for itest , subs in enumerate ( self . testconsistentvalues ) : if getsubs is not None : # have to explicitly evaluate since testsubs can be very complex subsvals = [ ( s , v . evalf ( ) ) for s , v in subs ] try : subs = subsvals + getsubs ( subsvals ) except self . CannotSolveError as e : # getsubs failed (sometimes it requires solving inverse matrix), so go to next set continue # have to sub at least twice with the global symbols A = Mall [ maxdegree ] . subs ( subs ) for i in range ( A . shape [ 0 ] ) : for j in range ( A . shape [ 1 ] ) : A [ i , j ] = self . _SubstituteGlobalSymbols ( A [ i , j ] ) . subs ( subs ) . evalf ( ) eps = 10 ** - ( self . precision - 3 ) try : Anumpy = numpy . array ( numpy . array ( A ) , numpy . float64 ) except ValueError as e : log . warn ( u'could not convert to numpy array: %s' , e ) continue if numpy . isnan ( numpy . sum ( Anumpy ) ) : log . info ( 'A has NaNs' ) break eigenvals = numpy . linalg . eigvals ( Anumpy ) if all ( [ Abs ( f ) > eps for f in eigenvals ] ) : try : Ainv = A . inv ( method = 'LU' ) except ValueError as e : log . error ( 'error when taking inverse: %s' , e ) continue B = Ainv * Mall [ 1 ] . subs ( subs ) for i in range ( B . shape [ 0 ] ) : for j in range ( B . shape [ 1 ] ) : B [ i , j ] = self . _SubstituteGlobalSymbols ( B [ i , j ] ) . subs ( subs ) . evalf ( ) C = Ainv * Mall [ 0 ] . subs ( subs ) . evalf ( ) for i in range ( C . shape [ 0 ] ) : for j in range ( C . shape [ 1 ] ) : C [ i , j ] = self . _SubstituteGlobalSymbols ( C [ i , j ] ) . subs ( subs ) . evalf ( ) A2 = zeros ( ( B . shape [ 0 ] , B . shape [ 0 ] * 2 ) ) for i in range ( B . shape [ 0 ] ) : A2 [ i , B . shape [ 0 ] + i ] = S . One A2 = A2 . col_join ( ( - C ) . row_join ( - B ) ) eigenvals2 , eigenvecs2 = numpy . linalg . eig ( numpy . array ( numpy . array ( A2 ) , numpy . float64 ) ) # check if solutions can actually be extracted # find all the zero eigenvalues roots = [ ] numrepeating = 0 for ieig , eigenvalue in enumerate ( eigenvals2 ) : if abs ( numpy . imag ( eigenvalue ) ) < 1e-12 : if abs ( numpy . real ( eigenvalue ) ) > 1 : ev = eigenvecs2 [ A . shape [ 0 ] : , ieig ] else : ev = eigenvecs2 [ : A . shape [ 0 ] , ieig ] if abs ( ev [ 0 ] ) < 1e-14 : continue br = ev [ 1 : ] / ev [ 0 ] dists = abs ( numpy . array ( roots ) - numpy . real ( eigenvalue ) ) if any ( dists < 1e-7 ) : numrepeating += 1 roots . append ( numpy . real ( eigenvalue ) ) if numrepeating > 0 : log . info ( 'found %d repeating roots in solveDialytically matrix: %s' , numrepeating , roots ) # should go on even if there's repeating roots? continue Atotal = None for idegree in range ( maxdegree + 1 ) : Adegree = Mall [ idegree ] . subs ( subs ) for i in range ( Adegree . shape [ 0 ] ) : for j in range ( Adegree . shape [ 1 ] ) : Adegree [ i , j ] = self . _SubstituteGlobalSymbols ( Adegree [ i , j ] ) . subs ( subs ) . evalf ( ) if Atotal is None : Atotal = Adegree else : Atotal += Adegree * leftvar ** idegree # make sure the determinant of Atotal is not-zero for at least several solutions leftvarvalue = leftvar . subs ( subs ) . evalf ( ) hasnonzerodet = False for testvalue in [ - 10 * S . One , - S . One , - 0.5 * S . One , 0.5 * S . One , S . One , 10 * S . One ] : Atotal2 = Atotal . subs ( leftvar , leftvarvalue + testvalue ) . evalf ( ) detvalue = Atotal2 . det ( ) if abs ( detvalue ) > 1e-10 : hasnonzerodet = True if not hasnonzerodet : log . warn ( 'has zero det, so failed' ) else : linearlyindependent = True break else : log . info ( 'not all abs(eigenvalues) > %e. min is %e' , eps , min ( [ Abs ( f ) for f in eigenvals if Abs ( f ) < eps ] ) ) if not linearlyindependent : raise self . CannotSolveError ( 'equations are not linearly independent' ) if returnmatrix : return Mall , allmonoms return exportcoeffeqs , origmonoms