title: correct way to map a function on entries down the pages in a 3D matrix?

text: **updated:** added another method to the methods I tried below for completion. This is new way of doing it, but most likely not very efficient, did not time it I can only visualize 3D matrices in my mind as a book made up of pages. Each page is a 2D matrix. The third dimension of the matrix is the number of pages in it. This can be easily carried over to higher dimensions. 4D matrix, is a volume of books, and so on. Suppose I have a 3D matrix in Mathematica like this:               a = {{1, 2, 3, 4}, {5, 6, 7, 8}};     b = {{9, 10, 11, 12}, {13, 14, 15, 16}};     c = {{17, 18, 19, 20}, {21, 22, 23, 24}};     m = {a, b, c}          (* {  {{1, 2, 3, 4}, {5, 6, 7, 8}},            {{9, 10, 11, 12}, {13, 14, 15,16}},            {{17, 18, 19, 20}, {21, 22, 23, 24}}        }     *)      and want to map `f` on all elements in each row/col but going down across the pages. From first page to last page of the book, then come back and do the same for next entry and so on. So that the result will be               {          {f[{1, 9, 17}],  f[{2, 10, 18}], f[{3, 11, 19}], f[{4, 12, 20}]},        {f[{5, 13, 21}], f[{6, 14, 22}], f[{7, 15, 23}], f[{8, 16, 24}]}     }      ![Mathematica graphics](http://i.stack.imgur.com/e28gc.png) This is explained in this diagram, which I very skillfully just made using Microsoft Visio ;) ![Mathematica graphics](http://i.stack.imgur.com/acR3H.png) There are the two methods I came up with after a bit of struggle. I'd like to find if there is a better or more elegant method. For example, I could not figure how to use `Map[f,mat, ??]` of the above to produce this result in one step. **Method 1**               Transpose[m];     f /@ Transpose[#] & /@ %      **Method 2**               Map[f[#] &, Table[m[[All, i, j]], {i, 2}, {j, 4}], {2}]       **Method 3 (added 8/20/13)**               Map[f, Partition[Flatten[m, {2, 3}], 4], {2}]      ## Summary and conclusions This problem can be broken up to 2 parts. The first part is how transform the 3D matrix with the correct list in each entry. Once this is done, the rest is easy, since now one can just do `Map[f,mat,{2}]` on the resulting matrix since this maps `f` to each entry. So the trick was is how to do the first part, as the second part I knew. All the answers are very good and learned from them, I selected Leonid's answer since it showed me more general way do this. I did not realized that `Transpose` can be used as a general transformation using its second arguments. Here is the possible transformations showing the one I wanted. This is useful, since now if I want different transformation, I can look this up, and just use `Transpose` again with different permutation. I am adding this to my Mathematica sheet cheat. ![enter image description here](http://i.stack.imgur.com/gmeAJ.png) **Appendix** On a side note. The main difference for me when working with Mathematica 3D matrices, is to remember that the first index is the page number and not the row number as the normal case with 2D matrices. In Matlab, in 3D, the page number is the third index (which is a little bit more natural to me). So in Mathematica, `m[[1,2,3]]` is the second row, third column in the first page. While in Matab to reference the same entry, one writes `m(2,3,1)`. Notice the 3rd dimension index location difference. So in Matlab, the first index is always the row and the second index is the column number, regardless if the matrix is 2D or 3D. Here is the same matrix in Matlab:               a = [1  2  3  4 ;  5 6  7  8];     b = [9  10 11 12; 13 14 15 16];     c = [17 18 19 20; 21 22 23 24];          m = zeros(2,4,3);          m(:,:,1)=a;m(:,:,2)=b;m(:,:,3)=c;          m      gives               m(:,:,1) =              1     2     3     4          5     6     7     8          m(:,:,2) =              9    10    11    12         13    14    15    16          m(:,:,3) =             17    18    19    20         21    22    23    24      It is not a big deal, just something to remember that is all. Ofcourse Matlab does not have a nice to way to use Map's and such so easily as in Mathematica, so it all evens out at the end :)


title: Program file deleted while running

text: I have a program running for some time at one computer. Today, unintentionally, I edited the program file from another computer (I use Dropbox), and the program kept running. Trying to correct my mistake, I deleted the program file, and, from the computer it was running, saved it again. My question is: have I ruined the evaluation? Should I start over or Mathematica can deal with such situation? My concern is that I do not know how long the evaluation will take, and it can be some months. Thanks in advance


title: Behavior of 1.2 ∈ Rationals changes in v9?

text: A friend of mine using version 9 recently showed me this: ![enter image description here](http://i.stack.imgur.com/WJa9Y.png) It's strange because in version 8.0.4 `1.2 ∈ Rationals` does return `False`: ![enter image description here](http://i.stack.imgur.com/PFubW.jpg) I checked the online document of `Element` and `Rationals`, no modification is mentioned. So it's a silent changing, or a bug?


title: Numbering element in descending order

text: How can I numbering the elements in descending order of second elements of sublists? For example, if I have the list which is               m={{1, 5}, {2, 8}, {3, 9}, {4, 2}, {5, 9}, {6, 7}, {7, 9}, {8, 10}, {9, 5}, {10, 2}};      The smallest second element is 2 of `{10,2}` so its number is 1.   Second smallest second element is also 2 of `{4,2}` so its number will be 2.   The output that I want is then: `{4, 6, 9, 2, 8, 5, 7, 10, 3, 1}`.   How can I do this?


title: Mathematica Compute Time Is Extremely Long (Integrate simple function with parameter scan)

text: I was wondering if you can help. I have a set of equations I need to solve for a list of parameters; entering each parameter value and solving anew each time is simply not economical. I found a way to vary the parameter and have Mathematica solve the equation, but the calculation now takes extremely long and does not even return a solution. Here is my code below (very simple): ![enter image description here](http://i.stack.imgur.com/oB4pu.png) What am I doing wrong? What can be done to make the calculation faster? Thank you in advance!


title: List comparison in constraints for Minimize

text: I am trying to solve an optimization problem, a simplified version of which is as follows:               x = Table[Symbol["x" <> ToString[i]], {i, 7}];     A = {3, 2, 5, 1, 7, 9, 6};     Minimize[{A.x,       And @@ Thread[0 <= x <= 1] &&        Plus @@ x == 3 &&        x \[Element] Integers &&        x != {1, 1, 0, 1, 0, 0, 0}       }, x]      Here, the sought-after solution is for `x`. The issue is that the solution when all but the last constraint (`x != {1, 1, 0, 1, 0, 0, 0}`) are used is exactly `{1, 1, 0, 1, 0, 0, 0}`. When I introduce the last constraint, there is no effect - I still obtain the same solution, although that constraint should remove that particular solution and output the next best one. If I replace the last constraint with `x == {0, 1, 1, 1, 0, 0, 0}`, then the output solution is exactly `{0, 1, 1, 1, 0, 0, 0}`, so the comparison seems to be evaluated, but for some reason, the inequality given by `x != {1, 1, 0, 1, 0, 0, 0}` does not evaluate to `False`, although it should. Any thoughts on this would be appreciated.sov


title: Plot counting function semi primes

text: I am trying to generate a plottable semi-prime counting function. Have tried:               DiscretePlot[Gather[{a = PrimeOmega[Range[100]];     b = PrimeNu[Range[100]];     Count[Transpose[{a, b}], {2, 2}]}     +     {a = PrimeOmega[Range[100]];     b = PrimeNu[Range[100]];     Count[Transpose[{a, b}], {2, 1}]}],{x, 0, 15}, Filling -> Bottom]]      but really have no clue as to where to go from here!


title: Map several functions in one routine

text: I have several functions, let's assume they are:               func1[x_]=x;     func2[x_]=3*x-5;     func3[x_]=0.1*x^2;      and a lot more like these. For each and every one of these I want to do the following               xvalues = Range[0, 500, 2.5];     points1 = Map[func1, xvalues];     Do[If points1[[[i]] < 0, points1[[i]] = 0, points1[[i]] = points1[[i]],       {i, 1,  Length[points1]}]     table1 = Transpose[{xvalues, points1}];      Now seeing as I have a lot of these functiones, is there any way to automate this in some kind of routine? While answering, please be aware that I don't really have any extensive Mathematica knowledge.


title: How to visualize large graphs with directed edge bundling

text: I'm wondering if anyone has implemented edge bundled graph plotting? http://www.win.tue.nl/~dholten/papers/forcebundles_eurovis.pdf


title: Time series (stochastic process) estimating parameters using characteristic function

text: I have a time series of assets ${A_1, A_2, ..., A_n}$, which is described by a sophisticated distribution having the following characteristic function: $\phi(u; t;\theta)$, where $\theta$ is a vector of unknown parameters. I need to esimate vector of unknown parameters $\theta$ of characteristic function. I tried to find a PDF using the inverse Fourier transform to use the maximum likelihood method, but the characteristic function is too complicated for that. I also thought about building the empirical characteristic function using the time series of assets and to estimate parameters using the least square method, but I do not know how to build the empirical characteristic function, because time series is not just a sample of random variables, it is a random process that depends on time.


title: Approximate minimum degree permutation algorithm in Mathematica

text: In MATLAB there is a nice implementation of the so called AMD (approximate minimum degree permutation) algorithm named amd (see Online MATLAB Documentation). There is an alternative algorithm called COLAMD (see Online MATLAB Documentation). This algorithm is especially interesting when one wants to model large resistor networks e.g. calculate path resistances or reduce a network of resistors (for the theory see, e.g., Efficient Methods for Large Resistor Networks). There the AMD step is a preconditioning step for a large matrix before a Cholesky decomposition (available in _Mathematica_ through _CholeskyDecomposition_ ) is done. The original fortran coded algorithms can be found on the webpage from the autors. Is anybody aware of an implementation of the AMD or the COLAMD algorithm in _Mathematica_?


title: Calling Mathematica from Java

text: I'm trying to write a Java program that uses _Mathematica_ code, but when I compile it, a window appears asking for MathLink. I enter               c:\\program files\\wolfram research\\mathematica\\9.0\\mathkernel.exe      but this causes an exception to occur. I don't understand what the problem is. Here is the code I am using:               import com.wolfram.jlink.*;     public class SampleProgram {         public static void main(String[] argv) {             KernelLink ml = null;             try {                 ml = MathLinkFactory.createKernelLink(argv);             } catch (MathLinkException e) {                 System.out.println("Fatal error opening link: " + e.getMessage());                 return;             }                  try {                 // Get rid of the initial InputNamePacket the kernel will send                 // when it is launched.                 ml.discardAnswer();                      ml.evaluate("<<MyPackage.m");                 ml.discardAnswer();                      ml.evaluate("2+2");                 ml.waitForAnswer();                      int result = ml.getInteger();                 System.out.println("2 + 2 = " + result);                      // Here's how to send the same input, but not as a string:                 ml.putFunction("EvaluatePacket", 1);                 ml.putFunction("Plus", 2);                 ml.put(3);                 ml.put(3);                 ml.endPacket();                 ml.waitForAnswer();                 result = ml.getInteger();                 System.out.println("3 + 3 = " + result);                      // If you want the result back as a string, use evaluateToInputForm                 // or evaluateToOutputForm. The second arg for either is the                 // requested page width for formatting the string. Pass 0 for                 // PageWidth->Infinity. These methods get the result in one                 // step--no need to call waitForAnswer.                 String strResult = ml.evaluateToOutputForm("4+4", 0);                 System.out.println("4 + 4 = " + strResult);                  } catch (MathLinkException e) {                 System.out.println("MathLinkException occurred: " + e.getMessage());             } finally {                 ml.close();             }         }     }      My objective is to use _Mathematica_ in java to calculate the maximum of a function. I don't know how to do that, and the documentation on the net is few.


title: Detecting a collision in n-body simulation with NDSolve

text: Consider this three-body system:               g = 10; (* gravitational constant *)     m = {100, 1, 1}; (* masses *)     μ0 = g*m; (* gravitational parameters *)     r0 = {{0, 0, 0}, {1, 0, 0}, {-1.1, 0, 0}}; (* initial location vectors *)     v0 = {{0, 0, 0}, {0, 30, 0}, {0, -30.1, 0}}; (* initial velocity vectors *)      Differential equations for velocities and positions:               jmax = Length[m];     v = Evaluate[#] & /@      Table[ToExpression["v" <> ToString[j]], {j, jmax}];     r = Evaluate[#] & /@      Table[ToExpression["r" <> ToString[j]], {j, jmax}];     eqns[μ_List] := Flatten[Table[{     v[[j]]'[t] ==       Sum[(-Normalize[r[[j]][t] - r[[i]][t]]*μ[[i]])/       EuclideanDistance[r[[j]][t], r[[i]][t]]^2, {i,         Delete[Range[jmax], j]}],     r[[j]]'[t] == v[[j]][t],     v[[j]][0] == v0[[j]],     r[[j]][0] == r0[[j]]     }, {j, jmax}]]     fns = Flatten[Table[{v[[j]][t], r[[j]][t]}, {j, jmax}]]      The numerical integration stops when two bodies collide:               μ = μ0;     sol = NDSolve[eqns[μ], fns, {t, 0, 10}][[1]]     vel = (fns /. sol)[[1 ;; -2 ;; 2]];     pos = vel = (fns /. sol)[[2 ;; -1 ;; 2]];      > NDSolve::ndsz: At t == 2.393483348087306`, step size is effectively zero; > singularity or stiff system suspected. >> The solution:               tmax = 2.393483348087306;     ParametricPlot3D[pos, {t, tmax - .1, tmax},       PlotRange -> {All, All, {-.01, .01}}, ImageSize -> 800,       PlotPoints -> 500, MaxRecursion -> 15]      ![collision plot](http://i.stack.imgur.com/yAqA5.png) I would like to detect collisions and modify the computation so that the masses (and in effect, gravitational parameters) and velocities of the collided bodies change: one body's new mass becomes the sum of both masses and the new velocity becomes the sum of both velocities; the other body's new mass and velocity become zero. I have tried to do it with `WhenEvent`:               μ = μ0;     sol = NDSolve[     Append[eqns[μ],      WhenEvent[Length[$MessageList] > 0, Print[tmax = t];      dist =        Table[EuclideanDistance[r[[j + 1]][t], r[[j]][t]], {j,          Length[r] - 1}];      distDif = Table[Abs[dist[[n + 1]] - dist[[n]]], {n, Length[dist] - 1}];      p = Position[distDif, Min[distDif]][[1, 1]];      μ[[p]] += μ[[p + 1]]; μ[[p + 1]] = 0;      v[[p]][t] += v[[p + 1]][t]; v[[p + 1]][t] = 0;      "RestartIntegration"]], fns, {t, 0, 10}][[1]]     vel = (fns /. sol)[[1 ;; -2 ;; 2]];     pos = vel = (fns /. sol)[[2 ;; -1 ;; 2]];      However, this does not work. It seems that I do not catch the error properly, the program aborts before the message can be caught. How to catch this properly?


title: How to draw a directed graph with arrows showing vertically from bottom to top

text: If I have a directed graph, how can I draw it with arrows pointing vertically from bottom to top, like showing a class inheritance pattern in OOP(object oriented programming)?


title: Plotting solutions of a 4th order polynomial equation

text: I have a polynomial equation of the fourth order, which has $4$ roots depending on a variable parameter `s1`. For each `s1` I have $4$ solutions. I need a `LinePlot` of all roots to see how they move in the complex plane when the parameter `s1` changes. I can solve my equation but I don't know how to plot all the roots together on one diagram (points connected with line for each solution).               poly = -6.110000000000001`*^6 k^4 + 1000.` s1^2 + 60.335263000000005`(-5.` k + s1)^2;          Table[ NSolve[ poly == 0, k], {s1, {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9}}]


title: Suppressing negative roots in Mathematica

text: **Problem** Using Mathematica's `Solve` operator can sometimes lead to an output involving a positive and negative root (say when solving for a variable such as `N^2`. **What Is Happening** Here I simply have a defined function `Radial` followed by implementation (`Solve`) of that function.               Radial[n_,l_,r_]:= LaguerreL[n-l-1, 2*l+1, ((2*Z)*r)/(n*a)] * ((Z*r)/(n*a))^l Exp[-((Z*r)/(n*a))];     Solve[Integrate[             N^2 Radial[1,0,r] r^2,              {r,0,\[Infinity]},             Assumptions -> {Z>0,a>0,{Z,a},Reals}           ]== 1,N]          Out[1]= {{N -> -(Z^(3/2)/(Sqrt[2] a^(3/2)))}, {N -> Z^(3/2)/(Sqrt[2] a^(3/2))}}      _Mathematica_ is kind enough to give me both positive and negative roots. **What I Want** I would like to suppress the negative root, leaving only the positive root solution, i.e. `{N -> (2 Z^(3/2))/a^(3/2)}` **What I've Tried** Here I include a couple things I have tried (some suggested by the kind people in chat).   1. Inserting `N>=0` into `Assumptions` : Resulted in two roots.   2. Inserting `N>=0` and `N` into the `Reals` section. : Resulted in two roots.   3. Inserting `&& N>0` following the `==1` combined with all possible permutations of Tries (1) and (2). : Resulted in no answer given (nor any errors). Hopefully there is a simple way of fixing this (and not some convoluted input syntax). However, for the benefit of the community, any suggestion is welcomed. For my particular case, the simpler the better.


title: Manipulating a continuous stream of sounds

text: I need some audio-feedback when a locator is moved on a surface. I have the following code that produces a beep every time the dynamic controller is interacted with such that the pitch of the sound correlates with the _z_ -value of the underlying surface of the `LocatorPane`. The problem is that when the locator is moved, it does not provide a continuous soundstream, rather jerky, disconnected notes. If I increase the length of each sound primitive then they overlap, and cause severe lags.               pt = {0, 0}; (* initialize locator *)          func = Sin[\[Pi]*2*#1*#2] &; (* arbitrary surface *)          sound[{x_, y_}] :=       Sound[SampledSoundFunction[        Sin[0.15/(Rescale[func[x, y], {-1, 1}, {1, 0}] + 1) #] &, 4000,         20000], .1];          LocatorPane[      Dynamic[pt, (pt = Clip@#; EmitSound@sound@pt) &],      ContourPlot[func[x, y], {x, -1, 1}, {y, -1, 1}, ImageSize -> 200]      ]      ![Mathematica graphics](http://i.stack.imgur.com/clLMW.png)


title: Subsets of a list

text: Consider the following demand of products for the next four months.               data={1,20,3,40};      I could produce the whole demand in the first months, leading to the following "production vector": `{63,0,0,0}`. As you can imagine, producing everything in the first month leads to high costs of inventory. Another possibility would be {21,0,43,0}. Now I would like to get all possible combinations of the subsets of `data`, leading to:               results=     {{{1},{20,3,40}},     {{1,20},{3,40}},     {{1,20,3},{40}},     {{1},{20,3},{40}},     {{1},{20},{3},{40}},     {1,20,3,4}     }      It is important that `Table[Map[Flatten,results][[i]]==data,{i,Length@results}]` is `True` for all `i`. This means that each sublist of `results` has the same order of figures as in data.


title: Using Evaluate and Evaluated -> True in Plot

text: A few use cases: (1) This, as expected, returns (almost) immediately:               f[x_?NumericQ] := Evaluate[(Pause[1]; x)]     Plot[f[x], {x, 0, 1}]        (2) This takes forever:               f[x_?NumericQ] := (Pause[1]; x);     Plot[f[x], {x, 0, 1}, Evaluated -> True]      (3) This is much faster than (1):               f[x_?NumericQ] := {(Pause[1]; x), x x};     Plot[f[x][[1]], {x, 0, 1}, Evaluated -> True]      (4) This returns immediately ... with a "wrong" plot:               f[x_?NumericQ] := {(Pause[a = x + 1]; a), x};     Plot[f[x][[1]], {x, 0, 1}, Evaluated -> True]      (5) While this, as expected, returns immediately, with the "right" plot:               f[x_?NumericQ] := Evaluate[(Pause[a = x + 1]; a)]     Plot[f[x], {x, 0, 1}]      I would like an "authoritative" explanation of cases (3) and (4) ... the rest I think I'm able to understand, but posted for completeness, and to get a complete answer that could be useful for others too.


title: Finding a Graph Shortest Path with Custom Graph Distance Function

text: What is the best way to find a shortest path by using a custom distance function? I'm trying to play a little bit with the concept of using distributions instead of a deterministic weight value for each given edge. The results can be different, if we want, for example, to check what is the shortest path, with 95% "confidence". (The Flaw of Averages.) Our knowledge of the distances might not be a certainty and we may want to work with an interval or confidence range for the weights.               loc = Flatten[Solve[Mean[WeibullDistribution[2, 1, x]] == #, x] & /@ {1, 2, 3.01}];     {n1, n2, n3} = WeibullDistribution[2, 1, #] & /@ loc[[All, 2]];     network = {1 \[UndirectedEdge] 2, 1 \[UndirectedEdge] 4, 1 \[UndirectedEdge] 5,                2 \[UndirectedEdge] 3, 2 \[UndirectedEdge] 4, 4 \[UndirectedEdge] 3,                4 \[UndirectedEdge] 7, 4 \[UndirectedEdge] 6, 3 \[UndirectedEdge] 7,                5 \[UndirectedEdge] 6, 6 \[UndirectedEdge] 7};     networkWeight = {n2, n3, n3, n2, n1, n1, n3, n2, n2, n2, n1};     edgeWeights = Map[Mean, networkWeight] // N;     edgeLabels = Thread[Rule[network, edgeWeights]];     g = Graph[network, EdgeWeight -> edgeWeights, EdgeLabels -> edgeLabels,      VertexLabels -> "Name", VertexLabelStyle -> Directive[Red, Italic, 20], ImagePadding -> 20,      GraphLayout -> "SpringEmbedding"]     FindShortestPath[g, 1, 7]      This will tell me that the shortest path from 1 to 7 is `{1,2,3,7}`. But if I calculate the shortest path by hand using Monte Carlo,                  path147 = RandomVariate[n3, 10000] + RandomVariate[n3, 10000];        path1237 =RandomVariate[n2, 10000] + RandomVariate[n2, 10000] + RandomVariate[n2, 10000];        dist147 = EstimatedDistribution[path147, WeibullDistribution[α, β, μ]];     dist1237 = EstimatedDistribution[path1237, WeibullDistribution[α, β, μ]];     NSolve[CDF[dist1237, x] == 0.95 && x > 0, x, Reals] // Quiet     NSolve[CDF[dist147, x] == 0.95 && x > 0, x, Reals] // Quiet      we can see that the path `{1,4,7}` would be our choice in this case. Is there a way to, after assigning the distribution to the `EdgeWeight`, to define a function that will calculate the total weight using the same concept shown above? Many thanks.


title: Difference in integral character heights

text: ![enter image description here](http://i.stack.imgur.com/YkhKN.png) Above is an example of what happens. It seems that the first integral is longer than the second because the y-character is longer than the x-character as this doesn't happen with other similar-height characters: ![enter image description here](http://i.stack.imgur.com/Bdgvk.png) Is there a way to change this so that the integral defaults to the larger size regardless of whether x or y is used? If it's any help: the stylesheet is **Outline** and the screen environment is **bulleted**. **Edit** \- In response to the comments, this was in a text cell (part of a set of math notes) **Edit** \- By default it is in TraditionalForm. Converting to StandardForm, then back to TraditionalForm makes them look the ideal way: ![enter image description here](http://i.stack.imgur.com/y4cdX.png)


title: Step-by-step definite integration

text: I am interested in obtaining the intermediate steps for rather straightforward symbolic definite integrations (e.g., from an elementary calculus text). As an example, for $\int_0^1 x\,dx$ I'd like output along the lines of $$ \int_0^1 x\,dx={x^2\over 2}\Bigg|_0^1={1\over 2} $$ I have reviewed the responses to this question dealing with WolframAlpha-like step-by-step output as well as this one and this one. The answer in the first link above is in the spirit of what I'm after (especially the last answer given by FDSg), but it focuses on _differentiation_. However, when I try this on even elementary definite integrals like $\int_0^1 x\,dx$ it simply returns the answer without the steps. Perhaps I am overlooking something. I do not intend to duplicate that question and this distinction is key to that. In short, I'd be happy if `ShowSteps` (applied to `Integrate[x,{x,0,1}]`) would simply show the output in the displayed equation above (the emphasis being that I don't need to see _how_ the indefinite integral is computed, just what it is, and then completing the evaluation between limits).


title: Isomorphisms for graphs with loops and multiple edges

text: I am working with graphs with multiple edges and loops, and I want to eliminate all isomorphic graphs from a long list I've generated. The `FindGraphIsomorphism` function is very nice, but only works for simple graphs. I'm looking to find such a function or its equivalent for multigraphs. For example, ![two multigraphs](http://i.stack.imgur.com/iYDi3.png) Is there possibly a third-party package, or maybe even an external program that would find such graph isomorphisms? Also, is it possible to compare graphics objects in Mathematica, to see whether they are similar? I have tried `SameQ`, but to no avail. ![graphs with self-loops](http://i.stack.imgur.com/36yfX.png) I am assuming this is the result because Mathematica plots the graphs with floating-point values, so it is impossible to have two identical graphs like this. So is there a function that looks at the similarity between two graphics objects? I think something like that might do the trick for me in most cases.


title: How to add a label to a CopyEvaluate button?

text: I found a definition of a button               Button[Defer[ Information[ \[SelectionPlaceholder], LongForm -> False]],        None, BaseStyle -> "CopyEvaluate"]      that is doing a good job. How can I add a label to this button definition similar to `Button[label, action]`? I have trid to change the definition given above, but nothing I've done worked.


title: Export in Manipulate, infinite loop

text: I want to make export of a plot inside Manipulate, but only once. However, this code runs infinitely and does export infinitely:               Manipulate[Module[{pl},       pl = Plot[Sin[a x + b], {x, 0, 6}];       Export["pl.pdf", pl];       pl]      , {a, 1, 4}, {b, 0, 10}]      If I remove Export[] from Manipulate, problem disappears, but I need it inside Manipulate for other reason;) Can we control number of evaluations of certain code inside Manipulate?


title: Positive integrand giving negative answer

text: I'm integrating a positive function f(t) times sin(t) from 0 to pi/5 and get -38. Actually f is slightly negative for a short time (smallest value ~ -0.0005), but far from enough to explain this. Here is the code for the function:               NormFact[k_, S_] :=      Sqrt[(4 S - 2 k - 1) (4 S - k) (4 S - k + 1) (k +        1) (k + 2)/(4*Pi*S)];          xi[k_, j_, S_] := (-1)^j*      Binomial[k, j]*(4*S - k - 1)!/((j + 2)!*(4*S - k - j - 1)!);          Ofunc[t_, S_, coeffs_] :=      1 - Cos[t/2]^(4*S) +      coeffs.Table[      NormFact[k, S]*       Sum[xi[k, j, S]*Sin[t/2]^(2*j + 2)*Cos[t/2]^(4*S - 2*j - 2), {j,          0, k}], {k, 0, Length[coeffs] - 1}];          cfs = {-2.6155133, 0.9614036, 0.100279, -0.432464,      0.39887624, -0.2508507, 0.10555472, -0.0007824, -0.0526054,      0.0703274, -0.0688152, 0.05138949, -0.03833719,      0.02062684, -0.0018939, -0.0027808};      And here's the result when I integrate:               In:= Integrate[Ofunc[t, 58.5, cfs]*Sin[t], {t, 0, Pi/5}]     Out:= -38.132 + 0. I      Can anyone see my mistake?


title: can mathematica plot graph/visualize rotation matrices

text: just wondering is there any example online   that show plotting of 3d rotation matrices like in the picture below (source) ![enter image description here](http://i.stack.imgur.com/quZ1h.png)


title: Custom operators; custom pattern matching with pure functions

text: What I'm trying to achieve in _Mathematica_ is the creation of a binary operator whose operands are both pure functions over the natural numbers. The result of the operator should be another pure function over the natural numbers. To demonstrate concretely what I want, suppose I have the following functions defined:               f[n_Natural]:=2*n;     g[n_Natural]:=n-1;      (There is no Head called "Natural" so the above pattern matching won't work. But I want _f_ and _g_ to accept only natural numbers. This is **problem No. 1** [SOLVED]) I then want a binary operator defined like so:               Needs["Notation`"];     CombinedFunction[f_NaturalFunction,g_NaturalFunction]:={#}/.{{x_Natural}:>f[#]+g[#]}}&;     InfixNotation[ParsedBoxWrapper["\[CirclePlus]"], CombinedFunction];      Operating $f$ $\oplus$ $g$ yields a pure function $h$ that only takes a natural number as an argument. I have found a way of enforcing the domain of $h$ thanks to this thread, but I want to extend this to ensure that $\oplus$ itself is only defined for unary functions over the natural numbers. Seeing as there's no Head like 'NaturalFunction', I don't know how to do this. This is **problem No. 2**. As an additional issue, the operator (which currently yields a function defined over the integers) currently gives an unsimplified output:               Needs["Notation`"];     CombinedFunction[f_, g_] := {#} /. {{x_Integer} :> f[x] + g[x]} &;     AddInputAlias["4" -> ParsedBoxWrapper["\[CirclePlus]"]];     InfixNotation[ParsedBoxWrapper["\[CirclePlus]"], CombinedFunction];          f=1&;     g=#&;     h=f\[CirclePlus]g      > >     {#1} /. {{x$_Integer} :> (1 &)[x$] \[LeftRightArrow] (#1 &)[x$]} & >   I would have expected the output to be: > >     (1+#)& >   I'm unsure of the inner workings of what I've written so I don't know how to obtain a simplified result. I can now apply $h$ to an integer and it operates as expected. However:               h[3.5]      > >     {3.5} >   I want instead _Mathematica_ to behave as if the function was simply undefined for anything but an integer, just as it would do if I defined $h$ like so:               Clear[h]; h[x_Integer]:=x+1;     h[3.5]      > >     h[3.5] >


title: Identifying subsets with minimal standard deviation

text: Consider the following:               data = {10, 11, 15, 12, 9, 8, 7, 6, 8};          DataAll=StandardDeviation[data];     DataMost=StandardDeviation[Most@data];       DataRest=StandardDeviation[Rest@data];              DataOtherCombinations=           Table[                 Sqrt[                      (StandardDeviation@data[[1 ;; -i]])^2                     +(StandardDeviation@data[[-i + 1 ;; -1]])^2                 ],         {i, 3, 8}]          ListPlot[Join[{DataAll,DataMost,DataRest},DataOtherCombinations], Joined -> True]      The plot shows, that the combination `{{10, 11, 15, 12},{9, 8, 7, 6, 8}}` (`i=6` in `DataOtherCombinations`) has the lowest total standard deviation. (I define the total standard deviation as `Sqrt[StandardDeviation[a]^2+StandardDeviation[b]^2]`). I was hoping to find the same result with `FindClusters` but this is not the case:               In[1455]:= FindClusters[data,DistanceFunction->EuclideanDistance]     Out[1455]= {{10, 11, 15, 12, 9, 8, 7, 6, 8}}      Maybe there is a way of customising `FindClusters` or even better: an approach based on the minimisation of the total standard deviation exists already. By the way: `DataOtherCombinations` does not represent all possible combinations. I'm looking for all subsets of `data` with two restrictions: the order of numbers must not be changed (i.e. the subset must be a sequence of `data`) and the minimum length of one sequence must be greater than or equal to `3`. Hence, another possible combination would be `{{10, 11, 15},{12, 9, 8},{7, 6, 8}}`. (I tried it with `Subsets` but this function seems more to return all permutations of all subsets.)


title: Import image using InputField and Button

text: I've got a pretty basic problem, which is probably easily solved, but for all my efforts I haven't found a solution yet: I want to make an Image importer of the following form: An `InputField` in which the name (either on computer or on web) can be put in, and next to it a button saying `"Import"`, and clicking of which imports, and displays the picture in the output. I've tried the following, but I keep getting an error message                DynamicModule[{x = "image name"},                     {InputField[Dynamic[x]],                      Button["import", Print[Import[Dynamic[x]]]]                    }                   ]      > >     Import::chtype: >   What am I doing wrong? Might be something very basic, I'm still pretty new to _Mathematica_ and `Dynamic` functions in particular.


title: Sophistication of Series[...]

text: I'll give a concrete example and I hope that my general question will be clear. Say I have three variables, $f$, $g$, $h$, and I know that $f=\mathcal O(x)$, $g=\mathcal O(x^2)$, $h=\mathcal O(x^3)$ ($\mathcal O$-notation) or something of that sort. Now, I want to do a series expansion of expressions involving $f$, $g$, $h$ to a consistent order in $x$. For example, if I write `Series[..., {x,0,4}]` I want to keep terms like `f^3` and `f g` but to throw away terms like`f g^2` and `g h`. Is there a simple way to do so?


title: ParametricPlot Question

text: I want to graph more than two, lets say 5 different `ParametricPlot`-s all in terms of different time variables. Is there any way to do this, such that they all show up on the same plot? For example, I want to plot:   * `{3Cos[a],4Sin[a]}`.. on some interval,   * `{4Cos[b],5Sin[b]}`.. on some interval,   * `{5Cos[c],6Sin[c]}`.. on some interval,   * `{6Cos[d],7Sin[d]}`.. on some interval,   * `{7Cos[e],8Sin[e]}`.. on some interval, all on one picture. Thank You Very Much.


title: How can I efficiently "remap" an image?

text: Many image processing libraries like OpenCV, Intel Performance Primitives or Octave have a useful function called "remap", that takes an image, an array with X coordinates and an array with Y coordinates, and returns an image that "transforms" the image by that geometric mapping. Or, as the IPP documentation puts it: _"Pixel remapping is performed using pxMap and pyMap buffers to look- up the coordinates of the source image pixel that is written to the target destination image pixel:_ > dst_pixel[i, j] = src_pixel[pxMap[i, j], pyMap[i, j]] The closest thing in Mathematica that I'm aware of is `ImageTransformation`, but that takes a function, not an array. Using arrays can often save a lot of time (for example, you can apply the same mapping-array to multiple image, and arithmetic operations on arrays are very fast). The best I've come up with so far is to combine `ImageTransformation` and `ListInterpolation` to convert the arrays to functions:               img = ExampleData[{"TestImage", "Lena"}];     mapX = Table[i + j, {i, 500}, {j, 500}];     mapY = Table[i - j, {i, 500}, {j, 500}];      (Imagine some time-consuming operation here. I know that passing `i+j,i-j` as a function to `ImageTransformation` is probably much faster in this simple case, but that's not the point.)               {xFn, yFn} =        ListInterpolation[#, {{0, 1}, {0, 1}}, InterpolationOrder -> 1] & /@ {mapX, mapY};          Timing[ImageTransformation[img, {xFn @@ #, yFn @@ #} &, {500, 500},        PlotRange -> {{0, 1}, {0, 1}}, DataRange -> Full]]      This works, but is extremely slow (4.7s on my PC). The second idea I had was to use `ListInterpolation` on the image directly (in this case, only on the red channel):               redFn = ListInterpolation[ImageData[img][[All, All, 1]],         InterpolationOrder -> 1];         Timing[Image[redFn[mapY, mapX]]]      This takes 1.7s for one color channel, so it's even slower for 3 channels. For comparison: The IPP's remap function usually takes a few milliseconds.


title: Memory not freed after running Clear when using a Table

text: When I'm creating large lists with `Table` I see the memory footprint increasing.               MemoryInUse[]/1024.     (* Out-> 16992.4 *)          data = Table[x, {x, 0, 10000}, {y, 10000}];         MemoryInUse[]/1024.     (* Out-> 408372. *)      I would expect that when I generate the same command `data = Table[...]` just as above the memory would not increase but in fact it just allocates a new block of memory.               data = Table[x, {x, 0, 10000}, {y, 10000}];        MemoryInUse[]/1024.     (* Out -> 799041. *)      Why is that? Even when I do:               Clear[data];     MemoryInUse[]/1024.     (* Out-> 799044. *)      Is this expected behaviour or am I missing something? I can't seem to find a function that clears the memory. I'm using Mathematica 8.0.4.0 on Windows 7


title: Automatic formatting of isotope symbols in a text cell

text: `InputAutoReplacements` is an obvious choice for formatting isotope names in text cells. I'd like to be able to type `"1H"`, for example, and when I hit the space bar have it replaced with `Row[{Superscript["",1],"H"}]` inline. But my initial attempts at this have not been productive. Is it possible to use string patterns in `InputAutoReplacements`, such as the pattern in the following `StringReplace`?               StringReplace["14C",        n:NumberString~~atom:(CharacterRange["A","Z"]~~Repeated[CharacterRange["a","z"],{0,1}]        ~~WordBoundary) :>         ToString[Row[{Style[Superscript["",Style[n,Smaller]],                            ScriptBaselineShifts->{Automatic,.4}],atom}],                 TraditionalForm]]      I expected that entering that pattern into the Option Inspector (Global preferences for `InputAutoReplacements`, using `RuleDelayed`) would cause occurrences of `"1H "` typed into a text cell to raise the "1" to a superscript. Nothing happens. You might imaging that formatting the superscripts according to the above rule by hand is tedious.


title: Is the CompiledFunctionCall WVM opcode efficient?

text: There is example in Mathematica documentation about `Compile`               cf = Compile[{{x}}, Sin[x]];     cf5 = Compile[ {{x}}, cf[x^2],         CompilationOptions -> {"InlineExternalDefinitions" -> True,           "InlineCompiledFunctions" -> False}];      The above code makes `"InlineCompiledFunctions"` `False`, thus prevents inlining. If we looks at the `CompilePrint` output,               <<CompiledFunctionTools`     CompilePrint[cf5]      will be                       1 argument             3 Real registers             Underflow checking off             Overflow checking off             Integer overflow checking on             RuntimeAttributes -> {}                  R0 = A1             Result = R2          1   R1 = Square[ R0]     2   R2 = CompiledFunctionCall[ Hold[CompiledFunction[{x}, Sin[x], -CompiledCode-]][ R1]]     3   Return      There is a `CompiledFunctionCall` opcode. **I don't know whether this call is as efficient as the "inline" version or not. The documentation didn't say.** I quote a paragraph here which I don't quite understand. > Here the external definition is used, but the compiled function is not > inlined. Instead it uses an efficient instruction to allow one compiled > function to call another. This type of call is important since it could > allow a compiled function to call itself, and when parallel execution is > carried out in the compiler the call can be done without any synchronization > locking. What does it mean? It says "it could allow a compiled function to call itself". But there is a example in the same documentation where the compiled function is a recursive call. For example,               cFact = Compile[{{x, _Integer}}, If[x == 1, 1, x*cFact[x - 1]],        {{_cFact, _Integer}},         CompilationOptions -> {"InlineExternalDefinitions" -> True}];     CompilePrint[cFact]      And what is **"synchronization locking"**?


title: Maximizing a function with assumptions

text: Using               f[s_] := Log[(s/r)^α ((α - 2) n0 r^α  + 2 π Pmax ρ r^2)                 /((α - 2) n0 s^α + 2 π Pmax ρ s^2)]/s      When I run the following line:               Assuming[ s > r && r > 1 && Pmax > n0 && n0 > 0 && ρ > 1 && α > 2,      Maximize[f[s],s]]      I get the following output:               Maximize[Log[((s/r)^α (n0 r^α (-2 + α) + 2 π Pmax r^2 ρ))                 /(n0 s^α (-2 + α) + 2 π Pmax s^2 ρ)]/s, s]      Actually, Mathematica ignores the assumptions and outputs the function only. That is, nothing is done by Mathematica. I know that `Assuming[.]` works with `Refine`, `Simplify`, and `Integrate`. But, is there any way to use `Maximize` with assumptions? What do you think about the given function? Is it analytically solvable at all? **EDIT:** Let `g[x]=n0 β + 2 π ρ Pmax x^(-β)`, where `β = α - 2`. We can rewrite `f[s]` as:               f[s_] := Log[g[r]/g[s]]/s      for `β>0`. Therefore, the optimal point is found by maximizing `Log[(g[r]/g[s])^(1/s)]` or equivalently maximizing `(g[r]/g[s])^(1/s)`. Mathematica cannot solve the latter optimization problem as well. Even with some more simplification, `r=1` and `n0=1`, which results in the following problem, it cannot be solved:               Maximize[{((β + 2 π Pmax ρ)/(β + 2 π Pmax s^-β ρ))^(1/s), β > 0, Pmax > 1, ρ > 1}, s]


title: How to protect pattern or subexpression when distributing / expanding expression?

text: I've got an expression like               expr = (1-x)(a+b)      that I would like to distribute / expand while keeping factors of `(1-x)` intact, i.e. the result should in the above example should look like               (1-x)a + (1-x)b      I know that for the explicit example given here, `Expand[expr,(a+b)]` would yield the desired result. However, I would need a solution where `(a+b)` can be any arbitrary algebraic expressions that is distributed with `(1-x)` being left untouched. Is there maybe a way to define a pattern that matches my `(1-x)` terms that I can than hold while distributing?


title: Delete a style type from stylesheet

text: I have created a new stylesheet but now I'd like to delete a style type, for example "program" from it. How can I do it ?


title: Efficient circular buffer?

text: I wish to create an efficient circular buffer. That is, I wish to keep a fixed length list while appending a new value and dropping the oldest, repeatedly. As we know lists in _Mathematica_ are implemented as arrays and `Append` etc. are slow on long lists. For example:               big = Range@1*^7;          Do[big = Append[Rest@big, RandomInteger@99], {100}] // AbsoluteTiming      > {2.2100031, Null} `Internal`Bag` and related functions are appropriate for a continually accumulating list but do not appear to be applicable to this situation. Is there an efficient means to have a large circular buffer in _Mathematica_?


title: Using a list of tuples in a pure function

text: I want to use a list of tuples within a function to make assignments. Say I want to make assignments of the form               value[i,j] = val      What I have is a list of tuples that give the relevant indices (i,j):               tuples = {{i[1],j[1]},{i[1],j[2]}, ... , {i[imax],j[jmax]}}      What I have found is that this statement will work in Mathematica 8 :               Function[t, value[t[[1]], t[[2]]] = 0.] /@ tuples      But using slots does somehow not work:               value[#[[1]],#[[2]]] = 0. & /@ tuples      ... gives an error. Where is the mistake I am making or is there no way to do this using slots? (using `Part[#,1]` and `Part[#,2]` will not help either]? Kind regards, Guido


title: Make Entire Area of New Graphic Editable

text: Create a new graphic via the menu bar: `Graphics -> New Graphic`. The `FullForm` expression for this graphic is:               Graphics[{}, ContentSelectable -> True, ImageSize -> {480, 360}, PlotRange -> {{0, 4/3}, {0, 1}}]      Using the Drawing Tools palette, the entire area of the newly created graphic is editable: ![fullyeditable](http://i.stack.imgur.com/Nv1iC.png) Suppose now I want to make a new graphic but make it taller; that is, I want a new graphic but with dimension 480 x 800. I can do so with `Graphics[{}, ImageSize -> {480, 800}]` but only a small part of it is editable. ![noteditable](http://i.stack.imgur.com/GWpdf.png) If I add a `PlotRange`, for example, `Graphics[{}, ImageSize -> {480, 800}, PlotRange -> {{0, 1}, {0, 4/3}}]`, the editable area increases. ![noteditable2](http://i.stack.imgur.com/21Eip.png) To make the entire area editable, I can change `PlotRange` to `{{0, 1}, {0, 5/3}}`. ![noteidtable3](http://i.stack.imgur.com/sQmfY.png) If I change `ImageSize`, is there a systematic way to choose `PlotRange` to make the entire area editable? Other approaches are also welcome.


title: Converting XMLElement objects to plain text

text: **Short version:** How can I convert an `XMLElement` that represents part of a HTML document to plain text? **Long version:** The more general problem is extracting information from webpages. We can import the webpage as an `XMLObject` and extract the relevant part. But this may still be a complex expression of many nested `XMLElement`s (several paragraphs, links, emphasis, etc.), while I'm typically only interested in the text. Let's take a random example: extracting the text from this article. Using the developer tools of any modern browser it's easy to find out that the relevant part is in a `div` with `id="article-body-blocks"`. So we do               page = Import[        "http://www.guardian.co.uk/science/blog/2012/nov/13/science-enforced-humility",        "XMLObject"];          body = Cases[page,         XMLElement["div", {"id" -> "article-body-blocks"}, ___], Infinity];      The body is still a compound expression. Is there a built-in, direct way to extract the text? My workaround is               ImportString[ExportString[First@body, "XML"], "HTML"]      but this is a hack (i.e. fragile, likely to break in future versions or with an input I didn't anticipate). Is there anything specifically meant for this purpose?


title: Using free-form linguistic input in a dynamic input field

text: I've been playing around with dynamic input fields; here's a basic example, which just consists of an input field where the user can enter a function and then a button that can be clicked to graph the function. What I'm wondering is: is there a way to allow free-form linguistic input to be entered in the input field? Right now, I must enter regular Mathematica syntax like Cos[x] to get a graph, but I'd really like to be able to enter something like "cos x". (Ultimately, I'd like my dynamic module to do more than just graph the function, which is why I don't want to just use Wolfram Alpha directly.)               DynamicModule[{graph, f, button},      graph = Plot[1, {x, 0, 1}, PlotRange -> {0, 0.1}, Axes -> False,         ImageSize -> 480, Background -> White];      f = Cos[x];           button = Button["Graph!",        graph = Plot[f, {x, -3, 3}, ImageSize -> 480, Background -> White];        ];           Panel[       Grid[{         {InputField[Dynamic[f]], button},         {Dynamic[graph], SpanFromLeft}},        Spacings -> {2, 1}, Alignment -> {Left, Baseline}        ]       ]      ]      Thanks!


title: An elegant way to plot a numeric function that returns a list, and have each element in a different color

text: I have a function that takes a numeric argument and returns a list of numbers. I want to plot each element of the list in a different color. If I use this command,               Plot[f[x],{x,-1,1}]      all the elements are plotted in the same color. The function takes only a numeric argument (its definition is `f[x_?NumericQ]:=...`) so I can't use `Evaluate` like in this question. So far I've been using this command:               Plot[{f[x][[1]],f[x][[2]],f[x][[3]],f[x][[4]]},{x,-1,1}]      Which works fine (since the function evaluates very fast, I don't mind there are redundant evaluations here, see this question). However this is not very elegant, and considering I have 16 elements to plot, it gets downright ugly. **Is there a more elegant way to plot each element in a different color?**


title: How can I create a colour image using raw black and white data? (FITS data cube)

text: I have a data cube containing pictures of a galaxy across a thousand wavelengths. I want to take the data for the image slides in the red and blue parts of the spectrum, and have them represented in the appropriate colours. If I could then superimpose the images somehow to show the red and blue regions of the galaxy together that would be fantastic. Currently I have been working with the raw data for each slide and the pixel magnitudes. To plot the galaxy at any given wavelength I've just being using a list density plot. I don't know how to map the range of images with a gradual colour change or how to combine them?


title: Exporting lists with two extensions to files

text: Continuation of this question: How to export if I need to have two extensions instead of one. E.g., in the previous case one can generate lists labeled L[1], L[2],.....and so on using a loop. But if I have 2 loops and hence two variables and the lists are labeled L[1,1], L[1,2],… then L[2,1], L[2,2],… How do I export them in wdx format?


title: Plotting two vectors on the unit circle

text: I would like to plot a vector on the unit circle, depending on an angle `Pi/k`, and the same vector rotated by an angle `d`. I try the following:               k = 0.9; d = 0.3;     pt = {Sin[Pi/k], Cos[Pi/k]};     pt2 = {pt[[1]]*Cos[d] + pt[[2]]*Sin[d],      pt[[2]]*Cos[d] - pt[[1]]*Sin[d]};          Graphics[Arrow[{{0, 0}, pt, {0, 0}, pt2}], Axes -> True,      AxesLabel -> {x, z}, PlotRange -> {{-1.1, 1.1}, {-1.1, 1.1}}]      In some cases, like when `k=0.9` and `d=0.3` as above, this almost works, except that the first vector has no arrowhead, only a line (see figure). In other cases, as when `k=1`, Only the second vector is drawn. Can anyone see what the problem is? ![The plotted vectors](http://i.stack.imgur.com/qLPxE.png)


title: ParallelTable 70 times slower on 16 cores than Table on single core

text: Consider these code(sorry it's messy), why the `ParallelTable` version is 70 times slower than the `Table` version?               Quit[]     Clear["`*"]       $Version       (*     ==> "8.0 for Linux x86 (64-bit) (February 23, 2011)"     *)        Kernels[]       (*     ==> {}     *)        a = RandomReal[{0., 1.}, {401, 300000}];         b = RandomReal[{0., 1.}, {401, 300000}];       Developer`PackedArrayQ /@ {a, b}      (*     ==> {True, True}     *)          ht = 2*0.375*^-9;         DFT[A_, ht_] :=        RotateRight[        ht/Sqrt[2 \[Pi]]*         Fourier[RotateLeft[A, Length[A]/2 - 1],           FourierParameters -> {1, 1}], Length[A]/2 - 1];          SmoothDFT[A_, ht_, n_] :=        DFT[Table[0., {(n - 1)*Length[A]/2}]~Join~A~Join~         Table[0., {(n - 1)*Length[A]/2}], ht];          SelectbyWRange[A_, {WMin_, WMax_}, {TakeWMin_, TakeWMax_}] :=       Module[{lthA, nMax, nMin}, lthA = Length[A];       nMin = Round[-((-WMax + lthA WMin)/(           WMax - WMin)) - ((1 - lthA) TakeWMin)/(WMax - WMin)];       nMax = Round[-((-WMax + lthA WMin)/(           WMax - WMin)) - ((1 - lthA) TakeWMax)/(WMax - WMin)];       Transpose[{Table[          TakeWMin + n *(TakeWMax - TakeWMin)/(nMax - nMin), {n, 0,            nMax - nMin}], Take[A, {nMin, nMax}]}]       ]          Smtx1 =         Table[SelectbyWRange[-Im[            SmoothDFT[a[[n]], ht, 2]*             Conjugate[SmoothDFT[b[[n]], ht, 2]]], {-834., 834.}, {19.5,            20.5}], {n, 1, 2}]; // AbsoluteTiming          (*     ==> {0.404922, Null}     *)       Kernels[]       (*     ==> {}     *)     LaunchKernels[]         (*     ==> {KernelObject[1, "local"], KernelObject[2, "local"],       KernelObject[3, "local"], KernelObject[4, "local"],       KernelObject[5, "local"], KernelObject[6, "local"],       KernelObject[7, "local"], KernelObject[8, "local"],       KernelObject[9, "local"], KernelObject[10, "local"],       KernelObject[11, "local"], KernelObject[12, "local"],       KernelObject[13, "local"], KernelObject[14, "local"],       KernelObject[15, "local"], KernelObject[16, "local"]}     *)        Smtx2 =         ParallelTable[         SelectbyWRange[-Im[            SmoothDFT[a[[n]], ht, 2]*             Conjugate[SmoothDFT[b[[n]], ht, 2]]], {-834., 834.}, {19.5,            20.5}], {n, 1, 2}]; // AbsoluteTiming          (*     ==> {25.459674, Null}     *)      Note that I have 16 cores on the node and the table have only 2 elements, i.e. `{n,1,2}`, even if I change to `{n,1,16}`, the `ParallelTable` version is still 10 times slower than the Table version. If it is the overhead, why it has a such huge overhead? Thanks. These are the screenshots: ![enter image description here](http://i.stack.imgur.com/uuv3A.png) ![enter image description here](http://i.stack.imgur.com/AUNYk.png) **Update** 1.As OleksandrR point out that there is no `DistributeDefinitions`, but in the documentation of `ParallelTable` it promises to automatically distribute the calculations (version 9): > ParallelTable is a parallel version of Table which automatically distributes > different evaluations of expr among different kernels and processors. > > The default value is `DistributedContexts:>$DistributedContexts` with > `$DistributedContexts:=$Context`, which distributes definitions of all > symbols in the current context, but does not distribute definitions of > symbols from packages. 2.OleksandrR also gave an excellent analysis in this post, pointing out that the performance problem is the `MemberQ` function. Indeed if we turn on the unpack warning, we can see it unpacks the array in call to MemberQ               On["Packing"];     Smtx2 = AbsoluteTiming[ParallelTable[         SelectbyWRange[-Im[SmoothDFT[a[[n]], ht, 2]*             Conjugate[SmoothDFT[b[[n]], ht, 2]]], {-834., 834.}, {19.5, 20.5}], {n, 1, 2}];]      > Developer`FromPackedArray::unpack: Unpacking array in call to MemberQ. >>               (*{25.557433, Null}*)      However a second evaluation is much faster even the array still unpacks               ClearSystemCache[]     Smtx2 = AbsoluteTiming[ParallelTable[         SelectbyWRange[-Im[SmoothDFT[a[[n]], ht, 2]*             Conjugate[SmoothDFT[b[[n]], ht, 2]]], {-834., 834.}, {19.5, 20.5}], {n, 1, 2}];]      > Developer`FromPackedArray::unpack: Unpacking array in call to MemberQ. >>               (*{0.156995, Null}*)      If we use the temporary fix of `MemberQ` proposed by Szabolcs in the same post, then the unpack warning is gone, but the evaluation is still slow.               Quit[]     On["Packing"];          (*need to reevaluate all the definition code above*)          memberQ[list_, form_] := Or @@ (MatchQ[#, form] & /@ list)     ClearAll[fix]     SetAttributes[fix, HoldAll]     fix[expr_] := Block[{MemberQ = memberQ}, expr]          Smtx2 = fix@       AbsoluteTiming[ParallelTable[SelectbyWRange[-Im[             SmoothDFT[a[[n]], ht, 2]*              Conjugate[SmoothDFT[b[[n]], ht, 2]]], {-834., 834.}, {19.5,             20.5}], {n, 1, 2}];]     (*{25.192359,Null}*)      So how to fix this?


title: Dealing with errors and resuming the code after one

text: How can I write my _Mathematica_ code so that it resumes at a specific place if any `Message` is displayed during the execution of an expression? For example, in Visual Basic for Applications there's `on error resume next`, or `on error goto`. I'd like also to be able to log the issue in a list so I can reproduce and correct the problem later.


title: Basic Row Traversal Question

text: Right now I'm using this `For` loop to do an operation to each row of a $n \times n$ matrix, I was trying to figure out how to do this with `Map` or `Take` or something like that.               evecs = RandomReal[1, {5, 5}]     For[{i, 1, Length@evecs}, i++, evecs[[i, All]] = evecs[[i, All]]/Norm[evecs[[i, All]]]];


title: Evaluating derivatives for use in solve

text: I need to evaluate the derivative of a function. My function looks like               approx[vector_, x_] :=  Sum[vector[[n]]*Basis[n - 1, x], {n, 1, Length[vector]}]      and I'm defining my derivative like:               ddappx[vector_, x_] := D[approx[vector, x], {x, 2}]       However if I evaluate it at a specific x I get an error:               General::ivar: 1 is not a valid variable. >>      I searched enough to know that a solution when plotting this is to evaluate it first. Something about HoldAll and so on. What I need to do is use it in a set of simultaneous equations. So I should be able to have:               Solve[{ddappx[{a,b}, 1]==0,ddappx[{a,b}, -1]==0},{a,b}]      To solve for {a,b} by setting the 2nd derivative of approx ==0 at x=1 and x=-1. However I get the               General::ivar: 1 is not a valid variable. >>      Errors again. Evaluating like                Solve[{Evaluate[ddappx[{a,b}, 1]]==0,Evaluate[ddappx[{a,b}, -1]]==0},{a,b}]      isn't working. so, any ideas? cheers EDIT: Having made the suggested change, that error is solved. However although evaluating approx works perfectly well, evaluating ddappx gives 0, no matter what I do. there is no error, just 0. I've added the whole code which goes in front just in case you think its relevant               diffeq[y_, x_] := y''[x] + y[x]     f[x_] := Cos[2*x]     Basis[n_, x_] := ChebyshevT[n, x]     approx[vector_, x_] := Sum[vector[[n]]*Basis[n - 1, x], {n, 1, Length[vector]}]          ddappx[vector_, x_] := Derivative[0, 2][approx][vector, x]      ....               approx[{0, 0, 1, 1}, x]      returns correctly:               -1 - 3 x + 2 x^2 + 4 x^3                     ddappx[{0, 0, 1, 1, 1}, x]      returns incorrectly:               0


title: Differential Equation power spectra

text: ε = .1; A = 0; ω0 = 1; ωf = 0;      data =        NDSolve[         {x''[t] == ε (1 - x[t]^2) x'[t] - ω0^2*x[t] + A*Cos[ωf*t], x[0] == 1, x'[0] == 0},          {t, 0, 100}];     Periodogram[data, ScalingFunctions -> "Absolute", PlotRange -> All]      How would I create a suitable Power spectrum of the frequencies that this system has. For these set of values (at the top), there should ideally be only one unique frequency (ω0) since the other one is not in play. Mathematica's example of power spectra utilizes the periodogram funcition, although I don't really understand what it does mathematically, so there is room for using another type of graph.


title: Using transformation rules to manipulate expressions

text: I'm trying out Mathematica as a replacement for LaTeX for typesetting math. It works very well, but I would like to take advantage of Mathematica's computer algebra capabilities as well. Specifically, instead of manually figuring out and typing the result of a transformation to an expression, I would like Mathematica to do it. For example, suppose I am trying to prove that associativity of multiplication holds in the complex numbers. The first step might be the following:               (a+bi)((c+di)(e+fi)) = (a+bi)(ce + cfi + edi - df)      How can I automate this using transformation rules? I can do the following:               In:  (a+bi)((c+di)(e+fi)) /. x_*(y_*z_)->x+y+z     Out: a + bi + c + di + e + fi      which works as expected, so the pattern is matching properly (I've tried various other simple transformations as well). However, the following does not work as expected:               In:  (a+bi)((c+di)(e+fi)) /. x_*(y_*z_)->x*(Expand[y*z])     Out: (a+bi)(c+di)(e+fi)      What am I missing? Am I even on the right track, or is there a better way to accomplish this?


title: How to read data from IBM 32-bit floating point binary file ? (not the common IEEE !!!)

text: I am trying to read a binary file. This file records seismic data. Basically its first 3600 bytes are for the general file information, following are for seismic traces, each begins with 240 bytes for information about this trace. For the file I am reading, it has 3617 traces, each trace has 6240 bytes, with the first 240 bytes for trace headers. Following the trace headers are some velocity values for underground formation, it has 1500 data points. I can view this data using some commercial software, the values should be around 4999.99 feet/s(velocity of sea water) to 14800 feet/s. But what I read range from 590.125 to 743.25, which is obviously not correct. Later I found out that the binary file I am trying to read is IBM binary file. If I read it as IEEE binary file, for integer data, it is all right, but for floating point data, what I read is totally wrong. The trace head I read is about the geographic location of trace, and it is correctly read by mathematica, which is {0, 20, 40, 60, ..., 72320} feet. I don't need to read all the trace headers. What I care most is to read the correct velocity values. Could someone help me? Thanks a lot! (* The following code do not work right away, you need to download and unzip the file first and replace the file string path, the file name is timodel_vp.segy.gz, it is about 22 Megabytes after unzip and absolutely free to download. The address is: http://software.seg.org/datasets/2D/Hess_VTI/ Sorry for the inconvenience! *)               nByteFile =        FileByteCount[        "e:/SeismicDataset/HessVTI/VelocityModelFiles/timodel_c11.segy"];          nByteTrace = (nByteFile - 3600)/3617;          nTraceSample = (nByteTrace - 240)/4;               str = OpenRead[        "e:/SeismicDataset/HessVTI/VelocityModelFiles/timodel_vp.segy",         BinaryFormat -> True];     pos = Range[3600, nByteFile - 1, nByteTrace];          traces = {SetStreamPosition[str, # + 76];        BinaryRead[str, "Integer32", ByteOrdering -> +1],        SetStreamPosition[str, # + 240];        Table[BinaryRead[str, "Real32",          ByteOrdering -> +1], {nTraceSample}]} & /@ pos;


title: Exploiting a relation among constants

text: Say I have two arbitrary real constants $p,q$ that satisfy $p^2+q^2=1$. Is there any way to easily declare this relation in _Mathematica_? **Edit** More specifically I have this (among others) expression: $$\frac{\sqrt{\frac{p^4 q^4(p^2-q^2)^2}{p^4 + q^4}}}{\sqrt{p^{10}+p^2q^8+p^8(-2+q^2)+2p^4q^4(-2+q^2)+q^6(-1+q^2)^2+p^6(1+2q^4)}}$$               Divide[Sqrt[Divide[p^4 q^4 (p^2 - q^2)^2, p^4 + q^4]],       Sqrt[p^{10} + p^2 q^8 + p^8 (-2 + q^2) + 2 p^4 q^4 (-2 + q^2) +         q^6 (-1 + q^2)^2 + p^6 (1 + 2 q^4)]]      and am asking whether there is a way to have _Mathematica_ reduce it using the relation $\;q^2+p^2=1$.


title: Integrate function without specifying limits of integration

text: I've been having trouble to make mathematica evaluate an integral without specifying the values of the limits of integration. I mean, if I write               Integrate[x,{x,x1,x2}]      Mathematica returns the correct answer: $\frac{(x2)^2-(x1)^2}{2}$. But If I write something like               Integrate[x^2/(1-x),{x,x1,x2}]      Mathematica just keeps the evaluation forever... Is there something that I can do to avoid this behaviour?


title: Strange behavior of BinCounts

text: Just try this example:               lst = {0, 1, 1, Sqrt[2], Sqrt[5], 2 Sqrt[2], 2, Sqrt[13], Sqrt[5], Sqrt[6], Sqrt[11],             Sqrt[14], Sqrt[11], Sqrt[10], 3 Sqrt[2], Sqrt[17], Sqrt[13], 2 Sqrt[5], Sqrt[29],             2 Sqrt[2]};          ex1 = BinCounts[lst, {Range[0, 18]}]     ex2 = BinCounts[lst, {Range[0, 19]}]     ex3 = BinCounts[N@lst, {Range[0, 18]}]      > >     {1, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} >     {1, 3, 6, 6, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} >     {1, 3, 6, 6, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} >   Seems that `BinCounts` failed to notice the irrational numbers in `ex1`:               Select[lst, IntegerQ]      > >     (* This result is consistent with ex1. *) >     {0, 1, 1, 2} >   Is it a bug? I'm using _Mathematica 8.0.4, Windows Vista 32bit_.


title: What pattern will match accented letters in a regular expression?

text: How should I modify the following such that the accented letter `É` is matched and the expression returns `True`?               StringMatchQ["Éta", RegularExpression["[a-zA-Z]+"]]


title: The shortcut for Subsuperscript doesn't work?

text: The document for Subsuperscript says that: > To enter a subsuperscript in a notebook, use either `Ctrl`+`_` to begin a > regular subscript or `Ctrl`+`^` to begin a regular superscript. After typing > the first script, use `Ctrl`+`%` to move to the opposite script position. > `Ctrl`+`Space` moves out of the subscript or superscript position. This just create a expresion that **looks like** a symbol with subsuperscript! Specifically, I typed `Ctrl`+`-` / `Ctrl`+`6` and a regular subscript / superscript appeared, `Ctrl`+`5` and the opposite script appeared, and after I moved the cursor out of the symbol by pressing `→` and checked its `InputForm`, what I saw is **not** something like               Subsuperscript[a, x, y]      but               Subscript[a, x]^y      _Snapshot:_ ![Mathematica graphics](http://i.stack.imgur.com/O8Nyb.png) What's wrong with it, I misunderstood the document? I use _Mathematica 8.0.4_ and _Windows Vista Home Basic 32bit_.


title: Optimizing the fitting of Ordinary Differential Equation(NDsolve) involving interpolation functions

text: I am a newbie to this forum. I am trying to use NDsolve to solve some ODEs involving the interpolation function. Then, I will compare the fitted curve and data to determine the ODE coefficients. Basically, I measure the Voltage(V) and Current(I) for a circuit. Then I want to use these V and I to determine the value of electric components in the circuit. First, I make a interpolation function out of the waveform V data:               dataIndex = Range[2, 10001];     myVoltData = Transpose[{rawData[[dataIndex, 1]], rawData[[dataIndex, 4]]}];     myCurrData = Transpose[{rawData[[dataIndex, 1]], rawData[[dataIndex, 2]]}];     voltageFunc = Interpolation[myVoltData, InterpolationOrder -> 2];      Second, use `NDSolve` to solve the ODE, then use the solution to calculate the chi-sq. In the following code, the ODE parameters to be fit are `r0Fit, r1Fit, c1Fit, currOffset`, and the function `fitFunc` will eventually return the chi- sq:               fitFunc[r0Fit_?NumberQ, r1Fit_?NumberQ, c1Fit_?NumberQ, currOffsetFit_?NumberQ] :=       Block[ {sol, curr00, curr11},       sol = NDSolve[ {voltageFunc[t] == curr00[t]*r0Fit + curr11[t]*r1Fit ,               curr00[t] == (curr11[t] + r1Fit*c1Fit/10^6*curr11'[t]),               curr00[t0fit] == curr11[t0fit] == (voltageFunc[t0fit])/(r0Fit + r1Fit)},              {curr00, curr11}, {t, t0fit, t1fit}, AccuracyGoal -> 7][[1]];         Apply[Plus,      ((myCurrData[[fitDataRange, 2]] + 50*currOffsetFit/10^6      - 50*Flatten[curr00[t] /. sol /.t -> myCurrData[[fitDataRange, 1]]])     /myCurrDataErr[[fitDataRange, 2]])^2/dataPnts]          ]      Finally, use `NMinimize` find the fitting result with chi-sq minimum:               Timing[fitResult = NMinimize[ {fitFunc[r0f, r1f, c1f,  iOffset],      130 < r0f < 170, 3000 < r1f < 3200, 650 < c1f < 750, 60 < iOffset < 100},       {r0f, r1f, c1f,  iOffset},AccuracyGoal -> 0, Method -> "DifferentialEvolution"]]       My first problem is, since `voltageFunc[t]` is the interpolation function from data, it severely procrastinates the `NDSolve` command, as a result, the `NMinimize` is consuming vast amount of time(>20 mins). I wonder, why does interpolation function of ODE delay `NDSolve` so much? Is there any optimization method to make the `NMinimize` and `NDSolve` faster? Secondly, `NMinimize` occasionally fails to find the global minimum. Any suggestions to improve the fitting result's reliability. Any suggestion is appreciated. Thank you all.


title: How to color lines in Epilog

text: I often use Epilog to label lines within graphs, and I like to coordinate the color of the labels with the color of the lines they're labeling. This is no problem. However, I'm currently working on a project in which I need vertical connecting lines at selected places, and these should color match as well. For the life of me, I can't figure out how to do this.               Plot[{Sin[x], Sin[x + 0.5]}, {x, 0, 2*Pi}, PlotStyle -> {Blue, Red},       Epilog -> {        Line[{{2, 0}, {2, Sin[2]}}],        Line[{{3, 0}, {3, Sin[3.5]}}]}]      ![enter image description here](http://i.stack.imgur.com/U0XFe.png) Any help would be appreciated.


title: Exponential function Integral

text: I have tried to solve this integral:               Integrate[E^(-((a^2 b c^2)/(a^2 + b)))/(a^2 + b)^2,a]      Mathemathica is not able to solve it, I have tried the integration by parts and it did not work, as well as some substitutions. Any idea how to tackle this problem?


title: Getting message DSolve::conarg:

text: I want to solve a simple linear PDE about $p_0$ which is a function of ξ, ζ and τ. However, the boundary condition is defined at `ζ = h`, which is a function of ξ and τ. The other quantities are all constant. In fact, this PDE can be solved by hand. I want to find p0 as a function of `h` and ζ. I just want to verify it, then I will try another more complex one. ![enter image description here](http://i.stack.imgur.com/GYHRo.jpg) The solution I want is as follows which is obtained by hand. ![enter image description here](http://i.stack.imgur.com/Aw0ft.jpg) I tried the following code, but it does not work.               DSolve[      {D[Subscript[p, 0][          h[ξ, τ], ζ], ζ] == (R/            Pr)*(1 - ζ*(h[ξ, τ] + k)^-1) - G, (2/3)*         Subscript[p, 0][h[ξ, τ], h[ξ, τ]] ==         e^2* d^(-1)*((h[ξ, τ] + k)^-1)^2 -          2*s*D[h[ξ, τ], {ξ, 2}]       },      Subscript[p, 0][       h[ξ, τ], ζ], {h[ξ, τ], ζ}      ]      I get the message > DSolve::conarg: The arguments should be ordered consistently. >>


title: Minimize amount of evaluations with fold/table/thread or something else

text: I can't figure out how to use fold/table/map/distribute or whatever to get this:               {f[f[{},a],a],f[f[{},a],b],f[f[{},a],c],f[f[{},a],d],f[f[{},a],e],     f[f[{},b],a],f[f[{},b],b],f[f[{},b],c],f[f[{},b],d],f[f[{},b],e],     f[f[{},c],a],f[f[{},c],b],f[f[{},c],c],f[f[{},c],d],f[f[{},c],e],     f[f[{},d],a],f[f[{},d],b],f[f[{},d],c],f[f[{},d],d],f[f[{},d],e],     f[f[{},e],a],f[f[{},e],b],f[f[{},e],c],f[f[{},e],d],f[f[{},e],e]}      That is I need to evalute f[{},n] n=a,...,e. For each result I need do the same again and again. So total evaluations for i iterations should be:               Sum[Power[5,n],{n,i}]


title: Butterworthfilter data filter

text: I am currently doing a signal analysis and filtering using Matlab. The filter I am using is Butterworthfilter model to filter scattered data (vector data). I am using Matlab as follows:               function [y]=fButterworthfilter(input_data,f_cut)     rng default;     [b,a] = butter(8,f_cut/(200/2));     y=filtfilt(b,a,input_data); % zero-phase filtering      I want to see how can I do the Butterworthfilter filtering using Mathematica. I tried ButterworthFilterModel in the help but the input is a function not a vector data. Let us say I have the following noisy data:               data1 = Table[Exp[-i^2], {i, -5, 5, .01}];     noise = RandomReal /@ RandomReal[{-0.3, 1}, Length[data]];     data2 = data1+noise;      How can I do Butterworthfilter filtering for this data? Thanks


title: CUDALink ran out of available memory

text: I have a problem with automatic GPU memory management when using CUDA. The situation is the following: I create a simple CUDA function that operates on a 2D array (image), which does not take much memory space. This 2D array is passed as input then modified and output again. The function is loaded into mathematica using `CUDAFuncitonLoad` (pseudocode):               CUDATestFunction = CUDAFunctionLoad[...,{"Float", 2, "InputOutput"},{16,16}];      The function compiles fine and works. Then I create some sort of loop that involves calling the CUDA-function mentioned before, assigning the output to a new variable and then doing the whole thing again till some number of iterations is reached (again only pseudocode):               Array2Dold = imageData; (*assign raw image data*)     For[i=1,i<maxiter,i++,         Array2Dnew = CUDATestFunction[Array2Dold];         Array2Dold = Array2Dnew;        ]      The problem I have now is that when the number of iterations gets high I am running out of GPU memory and I get the error-message:"CUDAFunction::outmem: CUDALink ran out of available memory, possibly due to not freeing memory using the memory manager". **Does anybody know why that happens?** My guess is that mathematica for some reason keeps the `Array2Dold` or some modified version of it after calculation in CUDA in GPU memory and does not free the memory when the CUDA-funtion (CUDATestFunction) is finished. I know that there is a CUDALink Memory Manager, but from what I read in the mathematica documentation mathematica should at least be able to handle the above situation automatically in a way that does not lead to a memory overflow. Real Code for the actual problem with some modification to simplify the input (which does not change the occurring problem):               (*data initialization*)     imageData = Table[RandomReal[], {i, 200}, {j, 300}]; (*this would usually be some GrayScale ImageData*)     f = Table[RandomReal[], {i, 200}, {j, 300}]; (*this would be some function evaluated on the imageData*)     maxiter = 10000;     testImageData = Table[{RandomReal[], RandomReal[]}, {i, 200}, {j, 300}]; (*would be the gradient of the imageData*)          (*the for loop that leads to the memory overflow*)     For[i = 1, i < maxiter, i++,      test = First@CUDATestFunction[testImageData, 0.1, f, imageGradientNormalized, 1, Sequence @@ Dimensions[imageData]];      testImageData = test;      ]          (*CUDACode*)     CUDATestFunction = CUDAFunctionLoad["       __device__ float length(const float2& a)       {         return sqrtf(a.x*a.x+a.y*a.y);       }       __device__ float2 operator+(const float2& a, const float2& b)       {         return make_float2(a.x + b.x, a.y + b.y);       }       __device__ float2 operator-(const float2& a, const float2& b)       {         return make_float2(a.x - b.x, a.y - b.y);       }       __device__ float2 operator*(const float& a, const float2& b)       {         return make_float2(a * b.x, a * b.y);       }       __device__ float2 operator/(const float2& a, const float& b)       {         return make_float2(a.x / b, a.y / b);       }       __global__ void resolvFs(float* p, float sigma, float* f, float* imageGradientNormalized, float lambda1, mint width, mint height)       {         int xIndex = threadIdx.x + blockIdx.x * blockDim.x;           int yIndex = threadIdx.y + blockIdx.y * blockDim.y;           int index = 2*(xIndex + yIndex * width);         if(xIndex < width && yIndex < height)         {             float2 vecP = make_float2(p[index], p[index+1]);             float2 vecN = make_float2(imageGradientNormalized[index], imageGradientNormalized[index+1]);                  vecP = vecP + 2 * sqrtf(lambda1) * vecN;             vecP = vecP/fmaxf(1, length(vecP)/(2*sqrtf(f[index/2]+lambda1)));             vecP = vecP - 2 * sqrtf(lambda1) * vecN;                  p[index] = vecP.x;             p[index+1] = vecP.y;         }       }", "resolvFs", {{"Float", 3, "InputOutput"},         "Float", {"Float", 2, "Input"}, {"Float", 3, "Input"},         "Float", _Integer, _Integer}, {16, 16}]


title: How to visualise an iterative approximation to pi

text: I'm experimenting with different algorithms that approximate pi via iteration and comparing the result to pi. I want to both visualise and perhaps know the function (if any) that describes the increasing trend in accuracy as the number of iterations rises. For example, 1 iteration might give me 3.0, 10 iterations might give me 3.12, 50 iterations might give me 3.1409 etc. The more iterations, the better. Sometimes it might be a case of diminishing returns, running 1000 iterations to get an additional decimal point of correctness - is it worth it? Knowing the sweet spot at which to quit would be nice to know, for each pi approximation algorithm type. I'm a mature programmer who is rekindling an interest in math, and this project idea interests me. Trouble is, I don't know what the correct math and statistical terminology is to describe what I want to do. Could someone please locate my problem within the history of math and possibly illustrate with mathematica? **UPDATE** : To be more specific. My pi generation algorithm produces the following pairs of numbers, where the first number is the number of iterations it looped and the second number is the abs(result - pi) i.e. how much the result differs from pi. The more iterations, the smaller the second number.               data = {{10, 0.19809}, {50, 0.039984}, {100, 0.019998}, {500,      0.00399998}, {1000, 0.002}, {20000, 0.0001}, {100000,      0.00002}, {500000, 4.*10^-6}}      I want to visualize what is going on with the data and find a function that describes what is going on. I'd like to also know what is possible in this area with this sort of data. What I've tried:               ListPlot[data, Joined -> True]      gives me an an almost unreadable graph due to the small numbers involved.               GeneralizedLinearModelFit[data,{x},{x}]      gives me FittedModel[0.0405405 -<<22>> x] which when plotted Show[ListPlot[data,PlotStyle->Red],Plot[%36[x],{x,10,500000}]] gives me a straight line. Not helpful as it doesn't seem to tell me anything about the trend re gradually converging on pi. I think there would be a curve of some sort that would better fit the data, but functions like FindFit seem to require you supply an expression e.g. Log and various parameters a, b etc - but how do I know what expression and what parameters to supply? Should I be looking at Interpolation instead? **UPDATE 2** : Thanks for the answers so far. What's holding me back from accepting an answer is the fact that I'm expecting a curved graph - not a straight line. Because the whole point of more and more iterations is that it gets closer to pi and so the answer must be a curve. And that curve must have a formulae, which it would be nice to have defined. So I did some experimenting in mathematica and came up with a curve after all. Here's how I did it. My data came from a crude Wells / Gregory / Leibniz algorithm (see Gregory Series) where I get pi by alternately adding and subtracting 1/n where n is 3, 5, 7... for example               N[(1 - 1/3 + 1/5 - 1/7 + 1/9 - 1/11 + 1/13 - 1/15 + 1/17) * 4]      which gives me 3.25237 - not bad for 8 numbers in the series, or 8 iterations if its a program that's generating the result. So I went on to define a Mathematica function that takes the number of iterations you want and generates pi. The following function also happens to subtract from real pi and return the discrepancy difference. The higher the iterations _it__ passed to this function, the smaller the difference from pi and the smaller the result coming out of this function:               f[it_]:=Abs[N[Sum[(-1)^(n-1)/(2n-1),{n,it}]]*4-\[Pi]]      So my original data points were for 10, 50, 100, 1000, 20000 and 100000 iterations. I thought those would be representative. Thus I can now generate my original data in Mathematica with               data = {{10,f[10]},{50,f[50]},{100,f[100]},{500,f[500]},{1000,f[1000]},{20000,f[20000]},{100000,f[100000]}}     output: {{10,0.099753},{50,0.019998},{100,0.00999975},{500,0.002},{1000,0.001},{20000,0.00005},{100000,0.00001}}      then plot it using the recommended techniques offered as comments and answers on this page               ListLogLogPlot[data,Joined->True,GridLines->Automatic,PlotStyle->Thick]      and I of course get a straight line like everybody else. Hmmm - remember I am after a curve. My instinct tells me it should be a curve. So then I discover that if I generate a more gradual list of points like this:               data = Table[f[x],{x,20}]     ListPlot[data]      then I finally get my curve. ![curved graph](http://i.imgur.com/XkKxiHF.png) Yey! It looks exactly like I imagined it would - curving gradually towards perfection. And if I change the iterations from 20 to 1000 or any higher number, I get the same curve, only smoother. So this is the function I want a precise definition of. Interestingly, if I use the ListLogLogPlot everybody is recommending on my more gradually spaced data               ListLogLogPlot[data, Joined -> True, GridLines -> Automatic, PlotStyle -> Thick]      I get a straight line again! Grrr. ![straight line again](http://i.imgur.com/5tPbiig.png) What I suspect is going on here is that there _was_ probably a curve in the normal ListPlot of my original data but we couldn't see it due to the scaling, which is why we went for the ListLogLogPlot. But that probably turns curves into straight lines (hence the log in the name?) - which disappointed me because I was expecting a visual curve. And I don't think my original choice of x values (representing the number of iterations being fed into my pi calculation function) and how those x values are spaced makes a difference. My x values of 10, 50, 100, 1000, 20000 and 100000 certainly made it hard to see the curve in the graph, but the curve was still probably there. My second, gentler set of x data points (from the same pi generation function) of 1..20 makes the curve ridiculously easy to see. So what counts as an answer to this question?   * An exact definition of the curve being seen here.    * And perhaps some commentary on what functions describing curves typically arise when attempting to visualise pi using iterative techniques. I hope I haven't changed the goal posts in this question. I think I have always asked for an exact definition of the curve in my data plus some commentary of what might be typical / best practice in this sort of visualisation/graphing territory, where there is a known value we are iterating towards.


title: FindMaximum inconsistency

text: The code below seems to work for n<11\. But for n=11, and above, it outputs `newa` then just outputs "beep" sound. WhyTheBeep says "The kernel Local has quit (exited) during the course of an evaluation.". On two separate occasions, (after making trivial changes) it gave the correct output for n=11. I then increasing n and it failed again. After setting n=11 again, it it failed to output and gave the "beep" sound again. Nothing was altered other than n between success and failure in this instance. Why the inconsistency/problem with n>10? (I'm using Mathematica 8.0 on Windows7) Any help would be much appreciated.               n = 10     c = 1.     e = 5.          newa = Table[(2 - i) Subscript[p, m] + (i - 1)/n, {m, 2 n + 2}, {i, 2}]          FindMaximum[     {((Subscript[p, 2 n + 2] - c (1 - Subscript[p, 2 n + 1]))/(Subscript[p, 2 n + 1])),     Subscript[p, 2 n + 2] == Sum[((n + 1) - m) Subscript[p, m] Subscript[p, n + m], {m, 1, n}],     Subscript[p, 2 n + 1] == Sum[Subscript[p, m] Subscript[p, m + n], {m, 1, n}],     Table[Subscript[p,m] == (Subscript[p,2 n + 1] (1 - Subscript[p, m + n]))/(((e*n)/(1 +   e)) Subscript[p, n + m] - (1 - Subscript[p, n + m])^2), {m, n}],     Table[Subscript[p, m] >= 0, {m, 2 n}],     Table[Subscript[p, m] <= 1, {m, 2 n}],     Sum[Subscript[p, m], {m, 1, n}] == 1     }, newa, MaxIterations -> 10000]