# Tag Info

38

General comments First, if you plan to use multi-dimensional integrals it is better to test with multi-dimensional integrals not with one dimensional ones. One might think that the test in the question is an appropriate one if multi-dimensional integration is done by the integrator in a recursive manner. This seems to be case for scipy.integrate.nquad (see ...

33

I don't think it is possible for FullSimplify to assess how far it is from a meaningful reduction of complexity. At every stage parts of the expression may after some transformation cancel (or not). Sometimes complexity has to be increased before it can be lowered (as I showed in this answer). Simplification is the process of minimizing complexity, where ...

24

First off, your function is very simple without any hard number-crunching, so it will always be hard to get a large speedup for the compiled version. Secondly, your Parallelization option for Compile is useless because it doesn't do any parallelization this way. Let me give slightly changed versions of your examples and explain how you can achieve a large ...

23

Chunks of weak compositions Here is slightly modified version of algorithm used in CombinatoricaNextComposition converted to a LibraryFunction. Needs["CCompilerDriver"] " #include \"WolframLibrary.h\" DLLEXPORT mint WolframLibrary_getVersion() { return WolframLibraryVersion; } DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData) { ...

23

I've always wondered about the scalability of MathLink (now officially "Wolfram Symbolic Transfer Protocol"). This is the protocol used by Mathematica to communicate between the front end and the kernel, and the basis of the Parallel` package. It has quite low bandwidth and high latency relative to, for example, MPI libraries. I also wonder how many MathLink ...

21

dat = {0.71, 0.685, 0.16, 0.82, 0.73, 0.44, 0.89, 0.02, 0.47, 0.65}; Module[{t = 0}, Split[dat, (t += #) <= 1 || (t = 0) &] ] {{0.71, 0.685}, {0.16, 0.82, 0.73}, {0.44, 0.89}, {0.02, 0.47, 0.65}} Credit to Simon Woods for getting me to think about using Or in applications like this. Performance I decided to make an attempt at a higher ...

19

It definitely has something to do with the Interpolation function. Evaluating tempdata = Import["http://www.inrim.it/~magni/cm.dat.gz", "Table"]; cmfunc = Interpolation[tempdata] we get the warning Interpolation::udeg: Interpolation on unstructured grids is currently only supported for InterpolationOrder->1 or InterpolationOrder->All. Order will ...

19

SatisfiableQ has three methods: "BDD": converts the expression to a BDD (binary decision diagram), "SAT": uses the Minisat library, "TREE": a branch-and-bound method based on the expression tree. SatisfiabilityCount counts instances by converting the expression to a BDD, so its timing should be close to SatisfiableQ with the "BDD" method (counting ...

18

I guess I should not have been surprised that there are actually many ways to estimate the Gaussian and mean curvature of a triangular mesh. I shall present here a slightly compacted implementation of the MDSB method, with a few of my own wrinkles added in. I should say outright that the current implementation has two weaknesses: the determination of the ...

17

Here is a summary of comments (before @ciao's best answer above), with a change in notation. These functions calculate the number of partitions of n into exactly k distinct parts of size at most m. NumberOfWays000[n_, k_, m_] := Count[Map[Length,Map[DeleteDuplicates, IntegerPartitions[n,{k},Range[m]]]], k] NumberOfWays001[n_, k_, m_] := Total[Boole[...

16

This seems pretty quick, particularly on larger cases / larger k, e.g. 451, 29, 101 finishes in a few seconds on the loungebook. N.B. - I have not tested this exhaustively, just thrown together from ideas... If[Min[#3, #1 - Tr@Range@(#2 - 1)] < 0, 0, SeriesCoefficient[QPochhammer[-x y, x, Min[#3, #1 - Tr@Range@(#2 - 1)]], {x , ...

16

You are right, it can be done in a fraction of second. One can explicitly construct an array of indexes blockArray[mat_] := SparseArray[ Tuples[Range@# - {1, 0, 0}].{Rest@#, {1, 0}, {0, 1}} &@Dimensions@mat -> Flatten@mat] Timings: matrices = RandomReal[1, {48, 128, 128}]; s1 = SparseArray@ ArrayFlatten@ReleaseHold@DiagonalMatrix[...

16

Vectorization will help a lot: a[x_?NumericQ] := N[Exp[-Abs[x]]]; x = Table[-10 + 0.02 (j - 1), {j, 1, 1001}]; A = Outer[a[#1 - #2] &, x, x]; // AbsoluteTiming (* {2.11988, Null} *) B = Exp[-Abs[x - #]] & /@ x; // AbsoluteTiming (* {0.016182, Null} *) A == B (* True *) Notice that I am doing arithmetic on vectors the size of x instead of scalars....

15

Here is a totally different approach based on the fact that successive products forming the generating function are due to multiplication by a binomial $1+t*z^j$. Form a matrix $v$ of zeros with $n+1$ rows and $k+1$ columns. Initialize the top left corner to 1. Iterate $v=v+w$ where $w$ is the matrix $v$ shifted down by $j$ rows and to the right by 1. The ...

15

The following is a much faster, but not optimal, recursive solution: pts = RandomReal[1, {10000, 2}]; f = Nearest[pts]; k[{}, r_] := r k[ptsaux_, r_: {}] := Module[{x = RandomChoice[ptsaux]}, k[Complement[ptsaux, f[x, {Infinity, .05}]], Append[r, x]]] ListPlot@k[pts] Some timings show this is two orders of magnitude faster ...

15

You are needlessly computing exactly the same BesselJ and BesselJZero function values over and over again. As an example, in the 100x100 case, BesselJZero[1 + l, 1] (l is a constant) is computed 551 times! You just need to compute each one once. The easy way to do that is to memoize: bj[n_?NumericQ, z_?NumericQ] := bj[n, z] = N[BesselJ[n, z]] bjz[n_?...

15

As is demonstrated very well in this post you can use a criteria for your pattern, thereby only applying your function as long as you are searching and not to all elements. Also there is a specific FirstPosition function. f[x_] := Module[{}, Pause[0.5]; 2 x] AbsoluteTiming[ Position[f /@ Range[10], 10, 1, 1] ] AbsoluteTiming[ FirstPosition[f /@ Range[...

14

This is best I can do so far. The system is linear so LinearSolve is a natural thing to try. arrays = CoefficientArrays[eqs, vars] (* {SparseArray[< 2 >, {9}], SparseArray[< 35 >, {9, 9}]} *) solv1 = Thread[ vars -> LinearSolve[arrays[[2]], -arrays[[1]], Method -> "CofactorExpansion"]]; // AbsoluteTiming (* {0.28347, Null} ...

14

I'm posting a whole new answer because I don't want to inherit any of the votes I received for my previous wrong answer. In formulating my new answer, I was aiming for correctness, simplicity, and reasonable (but not stellar) performance. Simplicity was achieved by taking a recursive approach, the clarity of which gives me confidence in the correctness of ...

14

pts = Partition[RandomReal[1, 10000], 2]; ListPlot[pts] Use SameTest option with Union pts2 = Union[pts, SameTest -> (Norm[#1 - #2] < 0.05 &)]; Length[pts2] 326 ListPlot[pts2]

14

You can accomplish this using DeleteDuplicatesBy, by first taking your two input lists and making a matrix out of them, and then deleting the rows where the last element (the element that came from myNewList) is a duplicate. Then you transpose back and assign the sublists of the reduced matrix to the new lists you want. {myListDuplicatesDeleted, ...

14

I am sure you can easily install also Linux on it and then you could contact Vladyslav Shtabovenko, the current maintainer of FeynCalc (https://github.com/vsht) and ask him about hard problems in High Energy Physics he would like to benchmark on such a King-Kong machine. Either him or somebody else could also provide you with more complicated examples of ...

14

It took me a while, but the suggestion of @Michael E2 was quite helpful, and especially the post (Optimize inner loops). For those of you (like me) who are new to this style of programming in Mathematica there are a few things that helped in my particular example. In my slow version I was looping over all vertices in the mesh list. For example in ...

13

Outer is highly optimized for several built-in functions (Plus, Times, List). Therefore Exp@-Abs@Outer[Plus, #, -#] &@Range[-10, 10, 0.02]; // RepeatedTiming (* {0.025, Null} *) gives ~50x speedup over Outer[#1 - #2&, #, #] and ~15x speedup over Outer[Subtract, #, #]. Also is a bit faster then Kuba's Exp[-Abs[x - # & /@ x]].

13

Select[ IntegerPartitions[24, {8}, Range[5]], #.# == 86 & ] {{5, 5, 4, 2, 2, 2, 2, 2}, {5, 5, 3, 3, 3, 2, 2, 1}, {5, 4, 4, 4, 2, 2, 2, 1}, {5, 4, 4, 3, 3, 3, 1, 1}, {4, 4, 4, 4, 4, 2, 1, 1}} Slightly more general approach (in case where IntegerPartitions is not what we need): ClearAll[ar, a]; ar = Array[a, 8] ar /. Solve[Flatten@{ Tr[ar]...

13

That's not entirely true: a // Dimensions {2251, 101, 2} So it's a 2d array of paris, and each row is treated as separate set to plot, ListPlot is trying to style them differently, etc. That's what takes time probably. You can keep pairs and do a fast plot reducing it to simple list of pairs. With Flatten[#,1] or Catenate: ListPlot[ a ] // ...

13

Here's a million points processed in half a second: SeedRandom[0]; (* updated for reproducibility *) xyCoordinateCentreCircle = RandomReal[1, {1*^6, 3}]; Map[ Mean[#[[All, -1]]] &, BinLists[xyCoordinateCentreCircle, {0, 1, 1/10}, {0, 1, 1/10}, {0, 1, 1}], {3}] // AbsoluteTiming (* {0.513721, {{{0.506174}, {0.497757},..., {0.50284}}, {{...

13

SelectFirst would work: SelectFirst[Range[10], f[#] == 10 &] You can tell that it only executes f on indexes 1-5 if you add a Print statement to the definition: f[x_] := (Print[x]; 2 x)

12

At the end of this post you'll find the code for a small benchmark to compare LibraryLink with standard passing vs LibraryLink with MathLink based passing on two counts: function call overhead; use a function that does nothing, has no arguments, returns nothing (llNone and mlNone for LibraryLink and MathLink, respectively) passing real arrays; use a ...

12

This seems to be about twice faster on large generic lists, and I made it somewhat faster still on integer lists (about 4-5x faster, per my tests): ClearAll[partitionToVectorLS]; partitionToVectorLS[list : {__Integer}, partitions_, sparsenessThreshold_: 10] := Module[{max = Max[list], min = Min[list], copy = partitions, sparseness, inds, nonsparseQ,...

Only top voted, non community-wiki answers of a minimum length are eligible