close issue #354; close issue #355 (#356)

* close issue #354; close issue  #355

Author: jverzani

Project.toml

@@ -1,6 +1,6 @@
 name = "SymPy"
 uuid = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
-version = "1.0.24"
+version = "1.0.25"
 
 
 [deps]

src/assumptions.jl

@@ -35,11 +35,13 @@ export ask
 ## We make a module Q to hold the assumptions
 ## this follows this page http://docs.sympy.org/0.7.5/_modules/sympy/assumptions/ask.html
 """
+    𝑄
+    SymPy.Q
 
-Documentation for the `Q` module.
+Documentation for the `SymPy.Q` module, exported as `𝑄`.
 
 SymPy allows for
-[assumptions](http://docs.sympy.org/latest/modules/sympy/assumptions/)
+[assumptions](https://docs.sympy.org/latest/modules/assumptions/index.html)
 on variables. These may be placed on free sympols at construction.
 
 For example, the following creates a real value variable `x` and a postive, real variable `y`:
@@ -53,11 +55,11 @@ The `Q` module exposes a means to *q*uery the assumptions on a
 variable. For example,
 
 ```
-ask(Q.positive(y))  # true
-ask(Q.negative(y))  # false
-ask(Q.positive(x))  # `nothing`
-ask(Q.positive(x^2)) # `nothing` -- might be 0
-ask(Q.positive(1 + x^2) # true  -- must be postive now.
+ask(𝑄.positive(y))  # true
+ask(𝑄.negative(y))  # false
+ask(SymPy.Q.positive(x))  # `nothing`
+ask(SymPy.Q.positive(x^2)) # `nothing` -- might be 0
+ask(SymPy.Q.positive(1 + x^2) # true  -- must be postive now.
 ```
 
 The ask function uses tri-state logic, returning one of 3 values:
@@ -67,14 +69,15 @@ The construction of predicates is done through `Q` methods. These can
 be combined logically. For example, this will be `true`:
 
 ```
-ask(Q.positive(y) & Q.negative(-x^2 - 1))
+ask(𝑄.positive(y) & 𝑄.negative(-x^2 - 1))
 ```
 
 The above use `&` as an infix operation for the binary operator
 `And`. Values can also be combined with `Or`, `Not`, `Xor`, `Nand`,
 `Nor`, `Implies`, `Equivalent`, and `satisfiable`.
 
-Note: As `SymPy.jl` converts symbolic matrices into Julia's `Array`
+!!! note: 
+    As `SymPy.jl` converts symbolic matrices into Julia's `Array`
 type and not as matrices within Python, the predicate functions from SymPy for
 matrices are not used, though a replacement is given.
 """
@@ -315,5 +318,22 @@ end
 
 end
 
+## Issue  #354; request to *not*  export  Q
 ## export
-export Q
+#export Q
+
+const 𝑄 = Q
+
+"""
+    𝑄
+
+Exported  symbol for  [`SymPy.Q`](@ref), a  Julia  module implementing `sympy.Q`. "Questions" can be asked through the patterns
+`𝑄.query(value)` (𝑄 is entered as  [slash]itQ[tab]) or `SymPy.Q.query(value)` *but  not* as `sympy.Q.query(value)`
+
+!!! note
+    At one time, the symbol `Q` was exported for this. To avoid namespace clutter, the unicode alternative is now used. Legacy code would need a definition like `const Q = SymPy.Q`  to work.
+
+"""
+𝑄
+export  𝑄
+

src/matrix.jl

@@ -81,7 +81,26 @@ function LinearAlgebra.eigvecs(a::Matrix{Sym})
     hcat((hcat(di[3]...) for di in ds)...)
 end
 
+# solve Ax=b for x, avoiding generic `lu`, which can be very slow for bigger n values
+# fix suggested by @olof3 in issue 355
+function LinearAlgebra.:\(A::AbstractArray{Sym,2}, b::AbstractArray{S,1}) where {S}
+
+    m,n  = size(A)
+    x =  Sym["x$i" for  i in 1:n]
+    out = solve(A*x-b, x)
+    isempty(out) && throw(SingularException(0)) # Could also return out here?
+    ret = Vector{Sym}(undef, n)
+    for (i,xᵢ)  in enumerate(x)
+        ret[i] =  get(out,  xᵢ, xᵢ)
+    end
+
+    return ret
 
+end
+
+function LinearAlgebra.:\(A::AbstractArray{T,2}, B::AbstractArray{S,2}) where {T <: Sym, S}
+    hcat([A \ bⱼ for bⱼ in eachcol(B)]...)
+end
 
 
 ##################################################

test/tests.jl

@@ -196,7 +196,28 @@ import PyCall
     @test_throws MethodError nsolve([z1^2-1, z1+z2z1-2], [z1,z2z1], [1,1])  # non symbolic first argument
     @test all(N.(sympy.nsolve([z1^2-1, z1+z2z1-2], [z1,z2z1], (1,1))) .≈ [1.0, 1.0])
 
+    ## issue 355: direct definition of LinearAlgebra.:\
+    @test_throws SingularException Sym[1 1; 1 1] \ [1, 2]
+    @vars  a b c d e f
+    A, b= [a b; c d], [e, f]
+    x = A \ b
+    @test  simplify.(A*x-b) == [0,0]
+    out =  Sym[1 1; 1 1] \ [1,1]
+    u =   free_symbols(out)[1]
+    @test out ==  [1-u,u]
+
+    
+    # Just a made-up example to test  if manageable
+    @vars a1 a2
+    n = 7
+    A = diagm(0 => ones(Int, n), 1 => fill(a1, n-1), 2 => fill(1, n-2), -1 => fill(a2, n-1))
+    b = Vector{Sym}(1:n)
+    x = A \ b
+    @test length(free_symbols(x)) ==  2
+    
+    
     ## limits
+    @vars x
     @test limit(x -> sin(x)/x, 0) == 1
     @test limit(sin(x)/x, x, 0) |> float == 1
     @test limit(sin(x)/x, x => 0) == 1
@@ -429,16 +450,16 @@ import PyCall
     end
 
     ## Assumptions
-    @test ask(Q.even(Sym(2))) == true
-    @test ask(Q.even(Sym(3))) == false
-    @test ask(Q.nonzero(Sym(3))) == true
+    @test ask(𝑄.even(Sym(2))) == true
+    @test ask(𝑄.even(Sym(3))) == false
+    @test ask(𝑄.nonzero(Sym(3))) == true
     @vars x_real real=true
     @vars x_real_positive real=true positive=true
-    @test ask(Q.positive(x_real)) == nothing
-    @test ask(Q.positive(x_real_positive)) == true
-    @test ask(Q.nonnegative(x_real^2)) == true
-    @test ask(Q.upper_triangular([x_real 1; 0 x_real])) == true
-    @test ask(Q.positive_definite([x_real 1; 1 x_real])) == nothing
+    @test ask(𝑄.positive(x_real)) == nothing
+    @test ask(𝑄.positive(x_real_positive)) == true
+    @test ask(𝑄.nonnegative(x_real^2)) == true
+    @test ask(𝑄.upper_triangular([x_real 1; 0 x_real])) == true
+    @test ask(𝑄.positive_definite([x_real 1; 1 x_real])) == nothing
 
 
               ## sets
@@ -549,8 +570,8 @@ end
     # issue 231 Q.complex
     @vars x_maybecomplex
     @vars x_imag imaginary=true
-    @test ask(Q.complex(x_maybecomplex)) == nothing
-    @test ask(Q.complex(x_imag)) == true
+    @test ask(𝑄.complex(x_maybecomplex)) == nothing
+    @test ask(𝑄.complex(x_imag)) == true
 
     # issue 242 lambdify and conjugate
     @vars x