- GiNaC-list - ginac.de

lsolve() fails with "system is not linear"
by Vitaly Magerya 06 Aug '18

06 Aug '18

Here's what I get with the latest GiNaC release, as well as the latest commit: $ ginsh ginsh - GiNaC Interactive Shell (GiNaC V1.7.3) [...] > lsolve({(1+c)^(-1)*(1+x+x*c)==0},{x}); lsolve: system is not linear The problem is that the linearity of the system is tested in lsolve() by subtracting 'x' times it's coefficient from the equation, and seeing if there are any more 'x' remaining. In this case such subtraction results in '-(1+c)^(-1)-(1+c)^(-1)*c*x- (1+c)^(-1)*x+x', where 'x' is not explicitly canceled, thus the error. One workaround that solves this particular case is using 'numer()' on each equation, like so: diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp index db40dba6..b83b2548 100644 --- a/ginac/inifcns.cpp +++ b/ginac/inifcns.cpp @@ -1094,7 +1094,7 @@ ex lsolve(const ex &eqns, const ex &symbols, unsigned options) matrix vars(symbols.nops(),1); for (size_t r=0; r<eqns.nops(); r++) { - const ex eq = eqns.op(r).op(0)-eqns.op(r).op(1); // lhs-rhs==0 + const ex eq = numer(eqns.op(r).op(0)-eqns.op(r).op(1)); // lhs-rhs==0 const exset syms = symbolset(eq); ex linpart = eq; for (size_t c=0; c<symbols.nops(); c++) { ... this doesn't look performant to me though, and doesn't obviously solve the general case. Any alternative ideas?

1 0

Debugging with GiNaC
by Patrick Schulz 05 Aug '18

05 Aug '18

Hi all, I'm trying to debug with GiNaC objects. I am aware of the dbgprint() function, however, this function can't be found for a numeric or a symbol. GiNaC::basic has this function, but it can't be found for numeric or symbol. The strange thing also is that gdb can't find any function listed in the class reference of, say, numeric. So the basic question is: How can i print the contents of GiNaC::numeric and GiNaC::symbol while debugging? Thank you very much, kind regards, Patrick

3 2

Using polynomial power series for solving PDEs
by Foad Sojoodi Farimani 03 Aug '18

03 Aug '18

Hello guys again, I'm not sure if you have received my former email but I'm gonna reply to that so you can read it as well: While ago I was trying to use polynomial power series to solve a system of partial differential and algebraic equations, when realized there is no implementation of the idea. There is only Mathematica's AsymptoticDSolveValue which is just for ODEs. So I decided to implement it myself. Thanks to the Sympy community we now have some progress. I have one implementation over here <http://nbviewer.jupyter.org/gist/celliern/b38158d04d9dc3d8079dc44e3b747ac8> by Nicolas CELLIER <https://github.com/celliern> and some ideas over here <https://cs.stackexchange.com/questions/95886/algorithm-for-using-power-seri…> by me and some of SymPy developers. I was wondering if we could join forces to come of with a general algorithm, then implementionation on diffrent languages shouldn't be that difficult. I was wondering if you could take a look at this question <https://cs.stackexchange.com/questions/95886/algorithm-for-using-power-seri…> and help us out. FYI I have made a issue on SymPy's github over here <https://github.com/sympy/sympy/issues/15015#> which you can follow. Best, Foad On Thu, Aug 2, 2018 at 9:40 AM Foad Sojoodi Farimani <f.s.farimani(a)gmail.com> wrote: > Hello everyone, > > I want to generate a multivariate polynomial given an array of nonnegative > integers D=[d1,...,dm]. I want to replicate what I have been suggested > in SymPy, in GiNaC, to see if I can get a better performance. (see here > <https://groups.google.com/forum/#!topic/sympy/V4KRzA-3doI>, here > <https://gist.github.com/Foadsf/07d4649f94f4d71876612a0ec6d9939b>and here > <https://stackoverflow.com/questions/51628270/how-to-generate-a-symbolic-mul…>) > I would appreciate if you could help me with this. > > Best, > Foad >

1 0

[patch] Improve automatic algorithm selection in matrix::solve()
by Vitaly Magerya 18 Jun '18

18 Jun '18

Hi, folks. Since we've recently added solve_algo::markowitz, I thought it would be a good idea to update the automatic elimination algorithm selection logic we have for matrix::solve() (currently living in matrix::echelon_form()). To this end I've constructed a benchmark that generates random-ish matrices and measures how fast matrix::solve() can solve them. The matrix generation works by starting with a known simple solution and then permuting and mixing rows and columns randomly. You can find the code for this procedure at [1]. With that code I've generated a bunch of data points and run a quick analysis on it. Overall, the fastest (on average) algorithm seems to be: For numeric matrices: If ncells > 200 && density < 0.5: Markowitz Otherwise: Gauss For symbolic matrices: If density > 0.6 && ncells <= 12: Divfree If density > 0.6 && ncells < 120: Bareiss Otherwise: Markowitz You can find the details (and some pretty charts) at [2]. The attached patch implements this recommendation. Note that matrix::determinant() also has a similar automatic algorithm selection logic, but I didn't touch it for now. [1] https://github.com/magv/ginac-matrix-solve-bench/blob/master/mx-solve-bench… [2] https://github.com/magv/ginac-matrix-solve-bench/blob/master/analysis.ipynb

2 2

[patch] Add 'algo' parameter to matrix::rank()
by Vitaly Magerya 12 Jun '18

12 Jun '18

Hi, folks. Currently matrix::rank() is hardcoded to use Bareiss elimination, which for some matrices is way too slow (and I had to work around this limitation in my own code). I propose adding an argument that will allow users to select which elimination scheme to use. The argument is the same as matrix::solve() takes (so, solve_algo::*). I've tried to mimic how a similar argument was added to matrix::inverse() by keeping the current set of rank() functions and defining separate ones with the 'algo' argument, with the old ones calling the new ones with solve_algo::automatic. This is instead of using default argument values. I've also put the code that makes the automatic selection of the elimination method into a separate function so that both matrix::solve() and matrix::rank() could share it. There's a bit of similar code in matrix::determinant(), but since the self of determinant algorithms is different from the elimination algorithms, it's not clear if these two could be unified.

2 2

Re: [GiNaC-list] [patch] Add solve_algo::markowitz for large sparse systems
by Vitaly Magerya 08 Jun '18

08 Jun '18

On 05/31/2018 05:55 PM, Richard B. Kreckel wrote: > Thanks a lot for your contribution! After doing some (probably too > simplistic) checks with your implementation of Markowitz-ordered > Gaussian elimination and finding no issues I've committed it. > > Regarding style, I've taken the liberty to adjust indentation a little > bit to comply with <https://www.cebix.net/indentation.html>. Thanks; I tried to mimic the style of the rest of the code, but it turned out to be challenging. Is there a set of vim options that implement this style? > Also, > please let's keep the continue statements in conditionals on a line of > their own. This improves debuggability. Apart from that, it's perfect. Sure.

1 0

ldegree changes after expand
by Matthias Naaf 04 Jun '18

04 Jun '18

Dear All, I've just hit an unexpected behavior of the ex.ldegree() method in my application. The problem is that ldegree() gives 0 before and 1 (the expected value) after calling expand() on the expression. I did not expect this, as the tutorial says that degree() and ldegree() "work reliably on non-expanded input polynomials". A small running example and its output are attached. The expression is the one I encountered in my application. This looks like a bug to me, or is it necessary to call expand() before using ldegree()? A related question: What about coeff(), should I call expand() before using coeff()? GiNaC version: 1.7.4 g++ version: 7.3.1 openSUSE (tumbleweed) Thanks in advance for any help on these questions Matthias

2 2

[patch] Allow exprseq in the arguments to lsolve()
by Vitaly Magerya 31 May '18

31 May '18

Hi, folks. Currently lsolve(eqns, vars) takes two lists as arguments, but then only accesses them via op(), so essentially any other container should work equally well. Here's a small patch to make it support exprseq in addition to lst. The motivation is that std::vector<ex>, upon which exprseq is based, can be more convenient to work with than std::list<ex>. Note that there's one additional places where lst is required, but exprseq would work too: ex::subs(). It would make sense to support exprseq there too. Also note that testing for lst is done via two different patterns in the source code: sometimes it is 'info(info_flags::list)', other times it is 'is_a<lst>'. It's not clear to me why this difference exists, or what should be the preferred way.

2 1

[patch] Add solve_algo::markowitz for large sparse systems
by Vitaly Magerya 29 May '18

29 May '18

Hi, folks. I've been trying to apply GiNaC to solving systems of sparse linear equations: it is currently too slow, and I have an improvement. The setup of the issue is this: I have a system of 2279 equations in 749 variables, which you can find in 'sys749.txt' at [1]. When I give this system directly to lsolve(), like this: // compile with: // c++ -std=c++11 -o sys749 sys749.cpp -lginac -lcln -ldl #include <ginac/ginac.h> #include <ginac/parser.h> #include <fstream> using namespace GiNaC; using namespace std; int main(int argc, char *argv[]) { symtab table; parser reader(table); ifstream i("sys749.txt"); ex data = reader(i); lst eqns; for (auto &&e : data.op(0)) { if (!e.is_zero()) eqns.append(expand(e) == 0); } cout << lsolve(eqns, data.op(1), solve_algo::gauss); } ... after ~5 hours of time and 1.5GB of RAM, I conclude that this will never terminate. In contrast, it takes ~2.5 minutes for Maple to solve this system, so I know we can do better. The main slowdown for GiNaC comes from the classical 'fill-in' problem during Gaussian elimination. Consider a sparse matrix like this: [* * * * *] [* * 0 0 0] [* 0 * 0 0] [* 0 0 * 0] [* 0 0 0 *] When you apply naive Gaussian elimination to this matrix, on the first step you will subtract the first row from the rest of the matrix -- this will make all the zero entries non-zero (this is the "fill-in" of zeros), thus making the matrix dense, and killing all the performance because each subsequent row elimination step will require N normal() calls. On the other hand, if you where to permute the matrix like so: [* 0 0 0 *] [0 * 0 0 *] [0 0 * 0 *] [0 0 0 * *] [* * * * *] ... then naive Gaussian elimination would fully preserve the sparsity structure, and at each step most of the operations you'll need to perform would amount to subtracting zeros, and could be optimized away (so, nothing to normal()ize). This sort of permutation to reduce fill-in is an essential part of solving sparse linear systems, and is the main reason for Maple's performance, I think. In fact, if you'll just sort the list of equations by the number of terms before passing that list into lsolve(), this already allows GiNaC to solve sys749.txt, in under an hour (don't remember the exact timing). We can do better than this, though. There's an overview of the current state of the art in this sort of techniques at [2]. Out of all the tricks described there, the most straightforward to me is the Markowitz ordering, which is a heuristic that says: at each step of Gaussian elimination, choose the pivot (i,j) that minimizes the amount of fill-in you'll need to perform in this step, approximated as (r[i]-1)*(c[j]-1), where r[i] is the number of non-zero entries in row i, and c[j] is the number of non-zero entries in column j. This heuristic is expensive to calculate, but it is well worth it for sys749. The patch I'm attaching implements Markowitz-ordered Gaussian elimination as a new solve_algo::markowitz. An additional part of the slowness is the part of lsolve() that constructs the coefficient matrix: it does so by calling coeff() repeatedly for every variable -- this alone takes about 10 minutes for sys749. To optimize this part I first obtain the full set of all symbols in an equation, and then call coeff() only if the symbol in question belongs to that set. This changes 10 minutes into 10 seconds. Still, I think there's more room for improvement here: it should be possible to collect the whole coefficient row in one pass over the equation, without calling coeff() multiple times. Since it takes about 1.3 seconds to read the whole expression from disk, it should take no more than that to determine the full coefficient matrix, I think. If you know how to achieve that though, I'd be interested in learning. In summary, with the proposed optimizations GiNaC can solve sys749 in 60-90 seconds (differs every time; at one point I was getting a consistent 45 seconds, but I've probably changed something since then), which is actually faster than Maple, at least for this particular system. I would appreciate some correctness testing on this new algorithm, because Markowitz ordering requires full pivoting, and that is trickier to get right than the partial pivoting which GiNaC has had used until now, especially when 'vars' in matrix::solve() is wider than n*1. I've modified the check_lsolve test to also test this algorithm, but it doesn't look comprehensive. Style changes or algorithmic suggestions are also welcome. [1] http://tx97.net/~magv/sys749.txt [2] http://faculty.cse.tamu.edu/davis/publications_files/survey_tech_report.pdf

1 0

psi (digamma) solving
by Matej Usaj 06 Mar '18

06 Mar '18

Hi All, I'm trying to solve function that contains a call to psi (digamma) but my call to fsolve fails every time, ie: > f=psi(x/2)==0; psi(1/2*x)==0 > fsolve(f,x,2.7,3.1); fsolve(): function does not evaluate numerically but I'd expect the result to be x~=2.923... Am I missing something here? Thanks for your help! Cheers, Matej

2 1