sursumCorda

1239 Reputation

15 Badges

2 years, 232 days

MaplePrimes Activity


These are replies submitted by sursumCorda

@Carl Love Many thanks. Surprisingly, I find that even if I don't compile it, the execution time of your implementation is still shorter than that of mine (well, tens of seconds faster).
Does loop-based code (instead of the so-called functional operation, which is used in my initial translation) typically provide better performance in similar cases?

To obtain some information about what the compiler is doing, you can set infolevel[Compiler] to a positive value (1 or 3). What can you find?

@dharr Many thanks. I believe that if GraphTheory:-AllPairsDistance (or networks:-allpairs) is fixed, this task will be more convenient to be finished.

@acer & @tomleslie Many thanks. I have just submitted a software change request.

@tomleslie According to updates,Maple2021,GraphTheoryGraphTheory['AllPairsDistance'] has had their efficiency improved. But somebody tells me that he cannot reproduce these in a "legacy" Maple. Can anyone compare the output from Maple 2020 and that from Maple2021 respectively?

@Axel Vogt & @Preben Alsholm Thanks. I believe that the built-in `sum` (as well as `int` mentioned in the previous question Integration problem) lacks something like automatic preprocessing and postprocessing, so the users have to complete such steps manually, which makes Maple‘s learning curve steeper.

@acer The much "ideal" result may be ; unfortunately, Maple cannot confirm this: 

:-simplify(diff(-x^(6-alpha)*Ei(alpha-5, beta*x), x)-x^(5-alpha)/exp(beta*x));
 = 
             (5 - alpha)                      
         -6 x            Ei(alpha - 5, beta x)

               (5 - alpha)                            
            + x            Ei(alpha - 5, beta x) alpha

               (6 - alpha)                            
            + x            beta Ei(-6 + alpha, beta x)

               (5 - alpha)             
            - x            exp(-beta x)

However, a less ideal result () can be obtained as follows: 

:-simplify(eval(convert(int(eval(exp(-beta*x)*x^(5-alpha), 5-alpha = a), x), GAMMA), a = 5-alpha));
 = 
          1   /        alpha  (-alpha)                   
      - ----- \(beta x)      x         (-GAMMA(6 - alpha)
            6                                            
        beta                                             

                                    \
         + GAMMA(6 - alpha, beta x))/
:-simplify(diff((GAMMA(-alpha + 6) - GAMMA(-alpha + 6, beta*x))*(beta*x)^alpha/(beta^6*x^alpha), x) - exp(-beta*x)*x^(5 - alpha));
 = 
                               0

@C_R 

What does the "*" between content and primpart do in a fucntional operator … 

It has been explained in evalapply (cf. `convert/algebraic`). But one may also use: 

(_ -> content(_) %* primpart(_))(sqrt(2*b - 4*a));
 = 
                    / (1/2)            (1/2)\
                  %*\2     , (-2 a + b)     /

 

@C_R As for "a command that factors sqrt(2) out of sqrt(2*a+2*b)", one may use: 

> (content*primpart)(sqrt(2*a+2*b));
 = 
                       (1/2)        (1/2)
                      2      (a + b)     

However, I don't think that these names are as common as factor (or the new command Physics:-Library:-FactorSum) is in general; Mathematica christens such symbols "FactorTerms", which is much easier to find for ordinary users (for example, if one searchs for factor, one will see "Factor", "FactorTerms", "FactorSquareFree" and so on).

@mmcdara Thanks for your advice. However, I find some strange things.
In the first approach, because

for matrices with non-numerical entries, linalg['definite'] returns a conjunction of Boolean expressions, all of which must be true if the matrix is to have the property specified by the parameter `kind`,

it is enough to search for a real instance satisfying the boolean expression. Strangely, the result is not correct.
As for the second approach, I attempt to resolve the conditions via brute-force methods: RealDomani[solve](And(RealDomain[solve](eig_1 >= 0), RealDomain[solve](eig_2 >= 0), RealDomain[solve](eig_3 >= 0), …)). Strangely, the result is still incorrect. 
 

restart``

mm:=<0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0|0,0,1,0,0,0,0,0,0,0,0,0,0,k__3-1/2,0,k__3-1/2,0,0,0,0,0,0,0,0,-1/2,0,1-2*k__3,0,-1/2,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0|0,0,0,0,0,1,0,k__3,0,-1/2,0,0,0,0,0,0,0,0,0,k__3,0,-2*k__3,0,0,0,0,0,0,0,0,-1/2,0,0,0,0,0|0,0,0,0,0,0,-2*k__3,0,k__3-1,0,-1/2,0,0,0,0,0,0,0,0,0,k__3,0,1-2*k__3,0,0,0,0,0,0,0,0,k__1,0,0,0,0|0,0,0,0,0,k__3,0,1-4*k__3,0,k__3-1,0,0,0,0,0,0,0,0,0,(1-k__2)/2,0,k__3-1/2,0,0,0,0,0,0,0,0,1-2*k__3,0,0,0,0,0|0,0,0,0,0,0,k__3-1,0,1-4*k__3,0,k__3,0,0,0,0,0,0,0,0,0,k__3-1/2,0,(1-k__2)/2,0,0,0,0,0,0,0,0,1-2*k__3,0,0,0,0|0,0,0,0,0,-1/2,0,k__3-1,0,-2*k__3,0,0,0,0,0,0,0,0,0,1-2*k__3,0,k__3,0,0,0,0,0,0,0,0,k__1,0,0,0,0,0|0,0,0,0,0,0,-1/2,0,k__3,0,1,0,0,0,0,0,0,0,0,0,-2*k__3,0,k__3,0,0,0,0,0,0,0,0,-1/2,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,-2*k__3,0,k__3,0,k__1,0,0,0,0,0,0,0,0,k__3-1,0,1-2*k__3,0,0,0,0,0,0,-1/2,0|0,0,k__3-1/2,0,0,0,0,0,0,0,0,0,0,k__2,0,-k__1+5*k__3-1/2,0,0,0,0,0,0,0,0,k__3-1/2,0,-k__1+5*k__3-1/2,0,1-2*k__3,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,k__3,0,1-2*k__3,0,k__3,0,0,0,0,0,0,0,0,k__3-1/2,0,k__3-1/2,0,0,0,0,0,0,-2*k__3,0|0,0,k__3-1/2,0,0,0,0,0,0,0,0,0,0,-k__1+5*k__3-1/2,0,k__2,0,0,0,0,0,0,0,0,1-2*k__3,0,-k__1+5*k__3-1/2,0,k__3-1/2,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,k__1,0,k__3,0,-2*k__3,0,0,0,0,0,0,0,0,1-2*k__3,0,k__3-1,0,0,0,0,0,0,-1/2,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0|0,0,0,0,0,k__3,0,(1-k__2)/2,0,1-2*k__3,0,0,0,0,0,0,0,0,0,1-4*k__3,0,k__3-1/2,0,0,0,0,0,0,0,0,k__3-1,0,0,0,0,0|0,0,0,0,0,0,k__3,0,k__3-1/2,0,-2*k__3,0,0,0,0,0,0,0,0,0,1-2*k__3,0,k__3-1/2,0,0,0,0,0,0,0,0,k__3,0,0,0,0|0,0,0,0,0,-2*k__3,0,k__3-1/2,0,k__3,0,0,0,0,0,0,0,0,0,k__3-1/2,0,1-2*k__3,0,0,0,0,0,0,0,0,k__3,0,0,0,0,0|0,0,0,0,0,0,1-2*k__3,0,(1-k__2)/2,0,k__3,0,0,0,0,0,0,0,0,0,k__3-1/2,0,1-4*k__3,0,0,0,0,0,0,0,0,k__3-1,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0|0,0,-1/2,0,0,0,0,0,0,0,0,0,0,k__3-1/2,0,1-2*k__3,0,0,0,0,0,0,0,0,1,0,k__3-1/2,0,-1/2,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,k__3-1,0,k__3-1/2,0,1-2*k__3,0,0,0,0,0,0,0,0,1-4*k__3,0,(1-k__2)/2,0,0,0,0,0,0,k__3,0|0,0,1-2*k__3,0,0,0,0,0,0,0,0,0,0,-k__1+5*k__3-1/2,0,-k__1+5*k__3-1/2,0,0,0,0,0,0,0,0,k__3-1/2,0,k__2,0,k__3-1/2,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,1-2*k__3,0,k__3-1/2,0,k__3-1,0,0,0,0,0,0,0,0,(1-k__2)/2,0,1-4*k__3,0,0,0,0,0,0,k__3,0|0,0,-1/2,0,0,0,0,0,0,0,0,0,0,1-2*k__3,0,k__3-1/2,0,0,0,0,0,0,0,0,-1/2,0,k__3-1/2,0,1,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0|0,0,0,0,0,-1/2,0,1-2*k__3,0,k__1,0,0,0,0,0,0,0,0,0,k__3-1,0,k__3,0,0,0,0,0,0,0,0,-2*k__3,0,0,0,0,0|0,0,0,0,0,0,k__1,0,1-2*k__3,0,-1/2,0,0,0,0,0,0,0,0,0,k__3,0,k__3-1,0,0,0,0,0,0,0,0,-2*k__3,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0|0,0,0,0,0,0,0,0,0,0,0,0,-1/2,0,-2*k__3,0,-1/2,0,0,0,0,0,0,0,0,k__3,0,k__3,0,0,0,0,0,0,1,0|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0>assuming'real':

(* First Approach *)

linalg:-definite(convert(mm, matrix), 'positive_semidef'): # unsolved constraints
SMTLIB:-Satisfy(`%`); # find an instance
LinearAlgebra:-IsDefinite(eval(mm, %), query = 'positive_semidefinite'); # False

{k__1 = -2, k__2 = 19, k__3 = -2}

 

false

(1)

(* Second Approach *)

RootOf~(factors(MTM:-poly(mm))[-1, .., 1], 'x'): # eigenvalues in implicit form
convert~(RealDomain:-solve~(% >=~ 0), `and`): # resolve: each of them ≥ 0
RealDomain:-solve(%, {k__ || (1 .. 3)}, 'maxsols' = 1);
if andmap(is, subs((temp := {k__1 = 5/3, k__2 = 35, k__3 = -7/3}), %)) then
    LinearAlgebra:-IsDefinite(eval(mm, temp), query = 'positive_semidefinite') # False again
fi;

{k__1 <= -2*(k__2*k__3+k__3^2-11*k__3+4)/(k__2-8*k__3+1), k__2 <= 6*k__3^2-3*k__3+1/2, 3*k__3^2+(9/2)*k__3+(1/2)*((4*k__3^2-4*k__3+9)*(3*k__3-2)^2)^(1/2) <= k__2, 6*k__3^2-k__2-k__3+1 <= k__1, k__3 < 1/2-(1/2)*17^(1/2)}

 

false

(2)

NULL


 

Download PSD_BUG.mw

 Is there something wrong?

 

@Carl Love Thanks for knowledge. The documentation of Alias gives me a feeling of being somewhat dangerous. Besides, some new commands in the ArrayTools package seem to compromise on efficiency (maybe for the sake of simplicity?).

@C_R Maple code is as follows: 

interface(printbytes = false): 
with(SignalProcessing, FFT): 
with(LinearAlgebra, HankelMatrix):
for n from 0 to 15 do
    m:=HankelMatrix(<$1+1..2**n+2**n>,datatype=complex[8],shape=[]);
    gc();print(n,time[real](FFT(m,normalization=none,inplace=true)))
od: # shape=[] is not optional

Now you can run Python directly (cf. Connectivity in Maple 2023); nevertheless, FFT is not built into Python.

from os import environ
environ['JAX_ENABLE_X64']="1"
import timeit,jax.numpy as jnp
for n in range(0x10):
	m=sum(jnp.indices((2**n,2**n),sparse=True,dtype='cdouble')+(1,1))
	print(n,timeit.timeit(lambda:jnp.fft.fft2(m),number=1),sep=",")

Here the better choice is using a standalone release.

@Carl Love "evalf[4](…)" is not enough to obtain the desired result. For example, the entry corresponding to 0.0 and 0.09 should be "0.5359" instead of "0.5358". So one has to use "evalf[4](evalf[4+1](…))" in the anonymous function.

@C_R Strangely, although the help page writes that

the code underlying these commands (in the Intel® IPP 2021.6.0 library) does not support Arrays of dimension greater than 1 directly,

this tutorial (updated on 10/26/2015) says that "there are 1D and 2D versions (to compute for real or complex data) in the library".
I believe that the documentation has some typos.

 

@C_R Thanks for your suggestion. But the Signal Processing package only works (and hence, only optimizes) on coerced Vectors with compatible datatype, because the documentation writes that 

when called with an Array of dimension greater than 1 (at most 5), SignalProcessing[FFT] calls the FourierTransform command in the DiscreteTransforms package instead.

Strangely, although it is said that "the code underlying these commands does not support Arrays of dimension two or higher directly", the Intel® IPP help page says that "there are 1D and 2D versions in the library".

First 11 12 13 14 15 16 17 Last Page 13 of 23