Skip to main content

Posts

Showing posts from May, 2015

programming - How to check if a 2D point is in a polygon?

Background: I use code from An Efficient Test For A Point To Be In A Convex Polygon Wolfram Demonstration to check if a point ( mouse pointer ) is in a ( convex ) polygon. Clearly this code fails for non-convex polygons. Question: I am looking for an efficient routine to check if a 2D-point is in a polygon. Answer Using the function winding from Heike's answer to a related question winding[poly_, pt_] := Round[(Total@ Mod[(# - RotateRight[#]) &@(ArcTan @@ (pt - #) & /@ poly), 2 Pi, -Pi]/2/Pi)] to modify the test function in this Wolfram Demonstration by R. Nowak to testpoint[poly_, pt_] := Round[(Total@ Mod[(# - RotateRight[#]) &@(ArcTan @@ (pt - #) & /@ poly), 2 Pi, -Pi]/2/Pi)] != 0 gives Update: Full code: Manipulate[With[{p = Rest@pts, pt = First@pts}, Graphics[{If[testpoint[p, pt], Pink, Orange], Polygon@p}, PlotRange -> 3 {{-1, 1}, {-1, 1}}, ImageSize -> {400, 475}, PlotLabel -> Text[Style[If[testpoint[p, pt], "True ...

performance tuning - Compiling Map over expression that yields a ragged array

I'm trying to speed up a function that looks in the neighborhood of each 3D point in a large dataset and finds all the points within 1 unit in each direction, x, y, z. I've started by using Select to find the points around a specified point and then Map this function over the dataset. First, lets make up some "data" and then define the findpoints function. data = RandomReal[10, {10^4, 3}]; findpoints[data_, pt_]:= Block[{x = pt[[1]], y = pt[[2]], z = pt[[3]]}, Select[data, x - 1 z - 1 ] Now lets map it over the random data. Map[findpoints[data, #] &, data]; // AbsoluteTiming This is quite slow, I gave up waiting, so I started by compiling the findpoint function. findpointsC = Compile[{{data, _Real, 2}, {pt, _Real, 1}}, Block[{x = pt[[1]], y = pt[[2]], z = pt[[3]]}, Select[data, x - 1 z - 1 ], CompilationTarget -> "C", RuntimeAttributes -> {Listable}, Parallelization -> True, RuntimeOptions -> "S...

plotting - How to save plots in grayscale

As sophomoric as this question seems, how should I save plots in grayscale in Mathematica? I generally like eps images for their scalability and I use ghostscript or other third party perl scripts to convert my images to grayscale. Answer One way would be to use ColorConvert to convert the RGB or Hue values to gray scale. Here's an example: Plot[{Sin[x], Cos[x], Exp[-x^2], Sinc[π x]}, {x, 0, π}] /. x : _RGBColor | _Hue | _CMYKColor :> ColorConvert[x, "Grayscale"] For 2D plots that accept a ColorFunction , you can simply use GrayLevel to get the plot in grayscale as: DensityPlot[ Sin[x ^2 + y^2], {x, 0, 3}, {y, 0, 3}, ColorFunction -> GrayLevel, PlotPoints -> 100 ] Typically, these grayscale plots are useful when submitting to journals that charge exorbitant prices just to print in colour. However, just a note of caution that discerning different shades of gray is not easy. For the most effect, it is recommended (at least in the journals I publish ...

compile - Compiling the VoigtDistribution PDF

According to List of compilable functions , Erf and Erfc are compilable functions. However, I want to make a compiled version of the PDF of a VoigtDistribution to use in a NonlinearModelFit , and it doesn't seem that the Erfc of a complex value will compile: funcReal = Compile[{{x, _Real}}, Erfc[x I], CompilationTarget -> "C", RuntimeOptions -> "Speed"] funcComplex = Compile[{{x, _Complex}}, Erfc[x I], CompilationTarget -> "C", RuntimeOptions -> "Speed"] Needs["CompiledFunctionTools`"]; CompilePrint[funcReal] (*same as funcComplex*) 1 argument 2 Real registers 4 Complex registers Underflow checking off Overflow checking off Integer overflow checking off RuntimeAttributes -> {} R0 = A1 C0 = 0. + 1. I R1 = 0. Result = C3 1 C1 = R0 + R1 I 2 C1 = C1 * C0 3 C2 = R0 + R1 I 4 C2 = C2 * C0 5 C3 = MainEvaluate[ Hold[Erfc][ C2]] 6 Return Note the call to MainEva...

scripting - Line by line processing with a WolframScript

I am having some difficulty understanding the correct syntax to enable a line by line processing of a file using a WolframScript . What I would like to do to mimic this simple python pipe in Mathematica. Take a python file myfunction.py : #!/usr/bin/env python import sys def myfunction(): for line in sys.stdin: data = line.strip().split("\t") print "{d[0]}\t{d[1]}".format(d=data) myfunction() What this script is doing is taking a text file with tab delimited columns, split at the tabs and then output the first two entries with a tab between them. Imagine a text file texfile.txt col11\tcol12\tcol13\tcol14 col21\tcol22\tcol23\tcol24 col31\tcol32\tcol33\tcol34 In a terminal run $> cat testfile.txt | ./myfunction.py which prints the output for each line in the file. So I would like to create a WolframScript file that does the same as myfunction.py and can be used in the same way. I've tried various combinations within the script such as: #!/u...

How can I implement object oriented programming in Mathematica?

Roman Maeder's object oriented programming package is nice, but I'm hoping someone can suggest a sleek and novel implementation that is easy to use. Answer Well, one obvious idea would be to build on the struct implementation by Bob Beretta. You would have to add information about methods and modify the implementation of --> to consider those as well, and for polymorphism, you'd also have to store the base class (or base classes, if multiple inheritance should be supported), and have --> look there if the field or method is not found.

performance tuning - File-backed lists/variables for handling large data

Background Currently I am working with some large data (most of it generated by Mathematica itself). I usually find it a hassle to do this. For example, I just exported a large amount of data to WDX on a machine with a lot of memory, just to find that I can't read it on my own machine (with little memory) because the file can only be read as a whole. It is also extremely slow to import (but using MX was not an option due to different architectures) Mathematica is excellent when working with in-memory data, as it's paradigms of operating on data as a whole ( Map , Outer , Table , etc.) are very convenient. But it is not great at working with data that is too large to fit into memory, and it is not good at sequential processing of on-disk data. There have been discussions about this (see the comment discussions on this and this question), and the following idea came up more than once: it would be great to be able to use Mathematica's native paradigms to work on large on-d...

plotting - ListPlot with lots of same couples of values

In my dataset of (x,y) values I have a lot of points with the exact same values on x and y: {{1,2},{3,4},{3,4},{3,4},{4,5}} So when I plot them with ListPlot it looks like I have only one point for {3,4}, which doesn't help me assess the real distribution of my data. What would be the best way to get a better sense of the distribution visually, to know that behind a single point are in fact hidden lots of other points? Answer try this also: t = {{1, 2}, {3, 4}, {3, 4}, {3, 4}, {4, 5}, {3, 4}, {1, 2}, {6, 7}}; ListPlot[Labeled[Style[#1, PointSize[#2/50], Hue[#2/10]], Style[#2, Bold, 12]] & @@@ (Tally[t])]

bugs - How to uncompress strings safely (without any evaluation)?

Compressed strings are often used to exchange Mathematica expressions on the Internet. It is, however, not easy to see what such expression will do after decompression. Resulting code can perform unwanted evaluations. Here is a (somewhat stupid) example: Warning: This will open a Youtube music video when decompressed! str = "1:eJx9jt0KgkAQhbUfiC6CegOvA70Pogu1CARBX6BNZ1Fad8Wd9e/p200ouvHmY87MnMNxniKhtmVZcqMRtoQpgkAXZnPQ8EVVC8XzsK8bkLIUnFrmtte4c4SGE/ZIyuyVCMYm20rj1pT5FGtUVEr8V+laD8bi/DyJYiB3egiAEsUwZFABR7qc6fIxbzVSwAAYGSCf62d/3weJUMU18PSoZYFYy5PndV3nDkKheoKbicrrCGbFpT0H0dg349UfYjkTn5r4/g2cBFpg"; To see the uncompressed content without evaluation, I have tried to use Defer (as suggested in this answer): Uncompress[str, Defer] but it does not prevent the video from being played. How can I uncompress strings from internet safely, i.e. without any evaluation taking place? Answer This is definitely a bug. Here's a slightly simpler workaround than those suggested in the commen...

calculus and analysis - How to make an organised investigation of branch cuts from a solution to a differential equation

I am attempting to solve two differential equations. The solution gives equations that have branch cuts. I need to choose appropriate branch cuts for my boundary conditions. How do I find the correct ones? The differential equations and the boundary conditions are ClearAll[a, b, x, ω, ν, U]; eqns = { ω b[x] - ν (a'')[x] == 0, U ω - ω a[x] - ν (b'')[x] == 0 }; bc1 = {a[0] == 0, b[0] == 0}; bc2 = {a[∞] == U, b[∞] == 0}; The boundary condition bc2 is ambitions and does not work if put into DSolve . However, with just boundary condition bc1 we can get a solution. sol = {a[x], b[x]} /. DSolve[Join[eqns, bc1], {a[x], b[x]}, {x, 0, ∞}] The solution is long and contains terms like (-1)^(3/4) which suggests four branch cuts. There are also constants of integration C[2] and C[4]. By playing around I find I can get a tidy solution by making substitutions and simplifying. I have replace C[2] with a normalised c[2] and similar for C[4]. I have replaced x with a norm...

Solving an integro-differential equation with Mathematica

I try to solve a nonlinear integro-differential equation with this code. Here i used a periodic condition. L=10; tmax = 2; NDSolve[{D[u[x, t], t] + u[x, t]*D[u[x, t], x] + D[u[x, t], {x, 2}] + D[u[x, t], {x, 4}] + 1/(2 L)*NIntegrate[D[u[xp, t],{xp, 3}]*Cot[\[Pi](x - xp)/(2*L)], {xp, -L, x, L}, Method -> {"PrincipalValue"}] == 0, u[-L, t] == u[L, t], u[x, 0] == 0.1*Cos[\[Pi]/L*x]}, u, {x, -L, L}, {t, 0, tmax}] which gives me NDSolve::delpde:Delay partial differential equations are not currently supported by NDSolve" The warning is understandable because the function u[xp, t] is still unknow when NIntegrate is evaluated. Note that we should use PrincipalValue here in NIntegrate because there is a singularity at $x=xp$ , which has been specified in the integration range.

algebraic manipulation - Factor a polynomial over the reals

I can't seem to find a function in Mathematica that factors a polynomial over the reals . Obviously, Factor doesn't work since it works over the integers, over $\mathbb{Z}_p$, or some algebraic fields extensions of the type $\mathbb{Q}(\alpha_1, \ldots, \alpha_n)$ where $\alpha_1, \ldots, \alpha_n$ are algebraic numbers. My question is general, although it arouse from my need to factor over $\mathbb{R}$ the polynomial $$ f(\lambda) = \lambda^4 + 17\lambda^2+12. $$ Of course, in this particular case I can find the roots and factor it myself. Answer One way is to find the roots, separate into real and complex, further separate the complex ones into conjugate pairs, then reform as a factorization (taking into account the leading coefficient). I'll illustrate with the given example. poly = x^4 + 17 x^2 + 12; roots = x /. NSolve[poly]; rroots = Select[rts, FreeQ[#, Complex] &] croots = Complement[roots, rroots]; topquad = Select[croots, Im[#] > 0 &] mult = Coeffici...

files and directories - Manipulate and FileNameSetter don't want to play nice together

In a Manipulate , one can specific the type of an individual control like so: Manipulate[v, {v, Red, ColorSetter}] I tried to do the same with FileNameSetter , i.e. Manipulate[v, {v, "no file chosen", FileNameSetter}] but that doesn't work, and the control is a slider. I can work around the issue like so: Manipulate[{v, FileNameSetter[Dynamic[v]]}, {v, "no file chosen"}, ControlType -> None] but I would rather understand why my second example above doesn’t behave as expected. Do you have any idea? Answer Manipulate does not directly support a control type called " FileNameSetter ", but fortunately it is possible to use custom controls (both in Manipulate and other functions), as described e.g. in the Control documentation: {u,func}        an arbitrary control from a function The trick to getting it working is that func must be a pure function. This is mentioned in the Advanced Manipulate Tutorial , in the section showing how to build a cus...

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...

plotting - Seeing coordinates in Graphics3D

I'm writing some structures built up from different strings that specify the coordinates of the lines. The code used to generate these is built up from a ton of user defined functions from my library, so posting them here doesn't work very well. However, I can provide an example here that'll suffice. What my question is about is how to show the coordinates in the structures that I draw, when looking at them with Graphics3D . When plotting functions, there's always this right mouse button option that has get coordinates, which will tell me exactly where I am on the plot. For Graphics3D this doesn't appear to be there. The underlying grid however definitely exists, as the lines are drawn at exact coordinates. So lets make an example. An example string based structure is x = {{Line[{{-273.5, 180., 0.}, {-273.5, 250., 0.}, {-237.5, 250., 0.}, {-237.5, 180., 0.}, {-242.5, 180., 0.}, {-242.5, 245., 0.}, {-268.5, 245., 0.}, {-268.5, 180., 0.}, {-273.5, 180., ...

files and directories - Path Names Longer Than 256 on Windows

Is there a way to work with pathnames that are longer than the typical 256 in Mathematica ? For example run the following in Cygwin echo 1 > 034567890345678903456789034567890345678903456789034567890345678903456789034567890345678903456789034567890345678903456789034567890345678903456789034567890345678903456789034567890345678903456789034567890345678903456789034567890345678903456789034567890 Then run FilePrint["034567890345678903456789034567890345678903456789034567890345678903456789034567890345678903456789034567890345678903456789034567890345678903456789034567890345678903456789034567890345678903456789034567890345678903456789034567890345678903456789034567890"] General::noopen: Cannot open ... Answer With a few rare exceptions, Mathematica is generally unable to work with long path names on Windows. This response presents two strategies to work around this difficulty: extended-length path syntax and short path names (documented in Naming Files, Paths, and Namespaces from t...

differential equations - how to simplify large expression with lots of special functions in it (BesselY, Hypergeometric, MeijerG etc...)

I saw this DE in Maple forum. When solving it using Mathematica 9.01, even though the result was correct (both solutions gave the same numerical answer for some random values), Mathematica's answer was much larger in size. Using Simplify or FullSimplify did not help much (expression was still order of magnitude larger than Maple's result). I think special simplification rules might be needed? I'd like to ask if there are other special simplifications to use, or commands I might be missing to reduce this result more. Below is the ODE and Maple and Mathematica results. Mathematica: ClearAll[y, x]; ode = y'''[x] + (1/x) y'[x] - (1/x^2) y[x] == 0; ic = {y[1] == 1, y'[1] == 0, y''[1] == 1}; sol = y[x] /. First@DSolve[{ode, ic}, y[x], x]; Simplify[sol] Result is too large fit in the margin of this answer, so I zoomed in Maple: ode:=diff(y(x),x$3)+(1/x)*diff(y(x),x)-y(x)/x^2=0; dsolve({ode, y(1)= 1, D(y)(1)=0, (D@@2)(y)(1)=1}); simplify(%,wronskian); ...

scripting - Running Mathematica Notebook files in command mode

The question is, how can I run a .nb file in the kernel mode of Mathematica? I am not an expert in Mathematica, but one of our users who use this program says that in the GUI mode, he selects all the cells (CTRL+A) and then evaluates the notebook (SHIFT+ENTER). However, he wants to run the program in background. When I test with math , the program quickly exits; however, in the GUI mode, the run time is very large actually. I read other documentation articles about that, but since I am not expert in Mathmatica, I have no idea! As an example, solve.nb file is an input to the command math -run . The output is also available here . I have no idea what the output means :| I simply tried to port the solution to Linux. So I wrote a solve.m file containing NotebookPauseForEvaluation[nb_] := Module[{}, While[NotebookEvaluatingQ[nb],Pause[1]]]; NotebookEvaluatingQ[nb_]:=Module[{}, SelectionMove[nb,All,Notebook]; Or@@Map["Evaluating"/.#&,Developer`CellInformation[nb]] ]; UsingFr...

Same name for functions in package and Global context

I am using a package that was written for Mathematica 3 while I am now working on Mathematica 8. I have issues with 2 functions in particular, Order and GraphComplement. They are a part of this package but are also present in Mathematica 8 by default but with distinctly different arguments and working. (I am assuming they weren't around in Mathematica 3) Now when I make use of some other functions in this package they make calls to both Order and GraphComplement and it keeps trying to use the in built wolfram version. I instead want it to use the one that comes with the package. Can anyone point me to the right direction? Answer If I understand the problem I think I can help, but what I propose is a bit weird. I think you have a package that uses the unqualified Symbol name Order internally, and you need the package not to see the System`Order Symbol while it is defined. To effect this you can temporarily change the Context of Order , then put it back after loading the package...

c++ - Link server WSTP API feature not receiving incoming connections

I'd like to apologize ahead of time that this may be a very long post but want to get all of this out there. Also, this post really only pertains to those with WSTP/MathLink experience. My problem is that I have written (note: many parts are straight from the documentation) a working program that creates a link server (i.e. WSLinkServer type) and starts listening for connecting links. Once the server receives one, it should then return control flow to my program and then activate, write to, and finally close the link and shutdown the server. Using the Command Prompt netstat command, I can clearly see my program listening, but when I try to connect to the server from a Mathematica notebook (or kernel for that matter) to test the server, such as link = LinkConnect["8000@127.0.0.1", LinkProtocol->"TCPIP"]; If[LinkReadyQ[link], LinkRead[link]] I then receive an error message and a $Failed return as seen here: LinkConnect::linkc: Unable to connect to LinkObject[...