Skip to main content

Posts

Showing posts from February, 2017

import - Converting XMLElement objects to plain 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 extr

linear algebra - How to know the usage of undocumented function like LinearAlgebra`BLAS?

BLAS is not documented in mathematica. Using ?LinearAlgebra`BLAS`* gives But None of the function has a detailed usage information Click any of the function for example, GEMM , gives At first I thought, BLAS in mma is belong to MKL, so I look up the usage in MKL reference manual, it says call gemm ( a , b , c [ , transa ][ , transb ] [ , alpha ][ , beta ]) the last four parameters are all optional. But in fact, if I run LinearAlgebra`BLAS`GEMM[a, b, c] mma tells me, it needs 7 arguments LinearAlgebra BLAS GEMM::argrx: LinearAlgebra BLAS GEMM called with 3 arguments; 7 arguments are expected. if I run LinearAlgebra`BLAS`GEMM[a, b, c, "N", "N", 1., 0.] mma tells LinearAlgebra BLAS GEMM::blnsetst: The argument a at position 1 is not a string starting with one of the letters from the set NTCntc. so the order of the arguments is not the same as MKL reference!! How should I know the correct order of arguments without trying several times? Are there detailed usage infor

differential equations - How to specify PDE Boundary condition on a B-spline?

Context I would like to solve a PDE on a boundary which is parametrized as a BSpline . I am trying to solve the force-free Grad-Shafranov equation on a boundary whose shape I do not know in advance. Specifically I need to solve for the toroidal flux of the magnetic field above an accretion disc. The Grad-Shafranov equation reads (in cylindrical coordinates) R D[P[R, z], {R, 2}] + R D[P[R, z], {z, 2}] - D[P[R, z], R] == - R/2; and I am seeking solution satisfying P==0 on a spline, see below. This question is related to the physical context of that question , where we try in to explain astrophysical jets like this: Eventually I would like to optimize the problem while changing the shape of the spline. First attempt I define my region via a BSpline: ff0 = BSplineFunction[pts = {{1, 0}, {1.2, 2}, {0, 2}}] So the upper envelope of the jet looks like this: pl0 = ParametricPlot[ ff0[t] // Release, {t, 0, 1}, Frame -> False, Axes -> False, PlotPoints -> 15, ImageSize -> Small

text - Skip header lines on import

When importing a data file what are the comment symbols for Mathematica ? That is, given a file like this blabla bulbul 1 2 6 54 7 ... .. what symbol do I have to put in front of header lines so Mathematica skips them and starts reading at the line 1 2 6 54 7 ... . I tried # , which works in gnuplot, but that did not work. I know that I could just tell Mathematica to skip the lines, but as I can control the file output, it would be nicer to use some kind of a tag. Answer Here is an approach that handles interspersed comments in addition to "headers" FilePrint["test.txt"] #comment #comment #comment 1 2 3 #c2 4 5 6 7 8 9 ImportString[ StringReplace[Import["test.txt", "Text"], StartOfLine ~~ "#" ~~ Shortest[___] ~~ EndOfLine ~~ "\n" -> ""], "Table"] {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}} of course you can invent whatever convention you want or even a mix, eg.. `StartOfLine ~~ {"#",

differential equations - Only final result from NDSolve

Finally, I started to play with differential equations in Mathematica . And I have faced the problem, which seems to me so basic that I'm afraid this question is going to be closed soon. However, I've failed in finding solution here or in the documentation center. My question is: how to set that NDSolve will not save whole InterpolationFunction for the result? I'm only interested in final coordinates or for example each 100th. Is there simple way to achieve that? Anticipating questions: I know I can do something like r /@ Range[0.1, 1, .1] /. sol at the end, but still, there is whole interpolating function in memory. I want to avoid it because my final goal is to do N-Body simulation where N is huge and I will run out of the memory quite fast. What is important to me is only the set of coordinates as far in the future as it can be, not intermediate values. I can write something using Do or Nest but I want to avoid it since NDSolve allows us to implement different solv

calculus and analysis - Problem with evaluation of the integral

I am trying to evaluate the following integral: Integrate[(1 + k^2)/(2 k Sqrt[1 + k^2 + k^4]) - 1/(2 k ), k] which does not work. Mathematica is just unable to find the primitive function. However, if I separate this integral into two integrals Integrate[(1 + k^2)/(2 k Sqrt[1 + k^2 + k^4]), k] + Integrate[-1/(2 k ), k] Mathematica yields 1/4 (ArcSinh[(1 + 2 k^2)/Sqrt[3]] - Log[2 + k^2 + 2 Sqrt[1 + k^2 + k^4]]) I don't understand why the first formula does not work. Am I missing some conditions or something? Answer Here are a few ways that work. These give answers that differ by a constant. Substitutions ClearAll[trysub]; SetAttributes[trysub, HoldFirst]; trysub[Integrate[int_, x_], sub_Equal, u_] := Module[{sol, newint}, sol = Solve[sub, x]; (* should check Solve *) newint = int*Dt[x] /. Last[sol] /. Dt[u] -> 1 // Simplify; Integrate[newint, u] /. Last@Solve[sub, u] // Simplify (* should check Solve *) ]; Try the expressi

plotting - Plot is discontinuous (it shouldn't be)

Plot[{(E^x)^(-2 I π) Hypergeometric2F1[-2 I π - 2 I Sqrt[2] π, -2 I π + 2 I Sqrt[2] π, 1 - 4 I π, -E^x] + (E^x)^ (2 I π) Hypergeometric2F1[2 I π - 2 I Sqrt[2] π, 2 I π + 2 I Sqrt[2] π, 1 + 4 I π, -E^x]}, {x, -10, 10}] Why? How to fix? Answer The numerical evaluation of your argument function leads to very small imaginary parts in the result that are due to numerical inaccuracy. Remove them by wrapping the argument of Plot in Chop (see its documentation ): Plot[Chop[(E^x)^(-2 I π) Hypergeometric2F1[-2 I π - 2 I Sqrt[2] π, -2 I π + 2 I Sqrt[2] π, 1 - 4 I π, -E^x] + (E^x)^(2 I π) Hypergeometric2F1[ 2 I π - 2 I Sqrt[2] π, 2 I π + 2 I Sqrt[2] π, 1 + 4 I π, -E^x]], {x, -10, 10}]

printing - How to Print a Cell Landscape in a Portrait Orientation Notebook?

I have a few cells in a report (notebook or CDF) that need to be printed "Landscape" where the rest of the report is printed "Portrait" . These cells are in different sections/locations of the report. I found the PrintingOptions option and its "PaperOrientation" setting. There is no information on where to use it so I tried it on ExpressionCell with no effect. ExpressionCell[ Plot[x^2, {x, 0, 1}, AspectRatio -> 5/10.5, Frame -> True, ImageSize -> { QuantityMagnitude@UnitConvert[Quantity[10.5, "Inches"], "DesktopPublishingPoints"], Automatic}], "Output", PrintingOptions -> {"PaperOrientation" -> "Landscape"}] The above cell contains a plot $10.5\times5$ inches. It prints portrait and chops off a portion of the right-hand side of the plot. I also tried adding page breaks before and after the cell. However, Print Preview still shows all cells printed "Portrait" .

graphics - Plot Ellipse based on EigenSystem

Please consider the following : covMatrice = {{34.925, -10.21}, {-10.21, 22.462}}; COG = Mean@(Nfixations["d"][[1]])[[All, 1]]; fixations = {{19.4688, 17.4281}, {18.0563, 21.7156}, {13.0219, 24.7219}, {22.9594,25.5219}, {28.5406, 24.6719}, {27.0688, 17.1656}, {27.6781,16.325}, {28.9281, 10.7719}, {16.025, 13.6625}, {19.1313, 17.1094}}; With[{ eigVec = Eigenvectors[covMatrice], eigVal = Eigenvalues[covMatrice]}, Graphics[{ White, Rectangle @@ frmXY, Black, Disk[#, .5] & /@ fixations, Red, Line[(# + COG) & /@ {eigVec[[1]]*eigVal[[1]]/5, {0, 0}, eigVec[[2]]*eigVal[[2]]/5}]}]] How could I draw an ellipse representing on the EigenSystem given that Neither Disk or Circle enable to implement an "orientation" ? A rough example of my desired output as drawn using PPT : Answer You can use Rotate

dynamic - How to create animated snowfall?

Well, the title is self-explanatory. What sorts of snowfall can we generate using Mathematica ? There are two options I suggest to consider: 1) Continuous GIF animations with smallest possible number of frames. 2) Dynamic -based animations. Answer My simple version using Image : size = 300; r = ListConvolve[DiskMatrix[#], RandomInteger[BernoulliDistribution[0.001], {5 size, size}], {1, 1}] & /@ {1.5, 2, 3}; Dynamic[Image[(r[[#]] = RotateRight[r[[#]], #]) & /@ {1, 2, 3}; Total[r[[All, ;; size]]]]] Update A slightly prettier version, same basic idea but now with flakes. flake := Module[{arm}, arm = Accumulate[{{0, 0.1}}~Join~RandomReal[{-1, 1}, {5, 2}]]; arm = arm.Transpose@RotationMatrix[{arm[[-1]], {0, 1}}]; arm = arm~Join~Rest@Reverse[arm.{{-1, 0}, {0, 1}}]; Polygon[Flatten[arm.RotationMatrix[# \[Pi]/3] & /@ Range[6], 1]]]; snowfield[flakesize_, size_, num_] := Module[{x = 100/flakesize}, ImageData@ Image[Graphics[{White, Table[Translate[

image processing - Measure a DensityHistogram[] pair similarity

I study human vision and more specifically eye-movements. "If we display 2 symmetrical patterns (20 min one after the other), will our gaze distribution be symmetric is my research question." The 2 figures at the bottom row below is what is displayed to subjects for 3 seconds. One pattern, then, later in the experiment, its symmetrical transform. Above are their respective Gaze Density histograms . That is the distribution of where their eyes were while observing the stimuli. The Blue square is the Center Of Gravity of the stimuli. How can I measure their similarity ? If I have some ideas, I think Mathematica offers great means of image analysis that could be used here. You could find here the data : allSymFix : 93 sublist for the 93 stimuli pairs I present, along with a manipulate to see all the histograms allSymFix[[1,1]] are all the gaze observed on stimuli 1 original version. allSymFix[[1,2]] on its symmetrical transform How can I measure the similarity within each al

formatting - Best way to give presentations with Mathematica

I have typically used PowerPoint or plain PDFs of slides to give presentations, but with heavy mathematical content, it can be tedious to create these presentations and make them look good. How can I best make use of Mathematica to give presentations? (I would prefer a slide-by-slide type format to what I've seen a few people do—using a regular notebook with the font size pumped up and collapsing/expanding sections as they go along.) Answer You can create SlideShows using Mathematica and run it to demonstrate presentation. Main advantage of using such Slideshow over Powerpoint / PDF is that you can play dynamic content. With CDF format available with Mathematica now , Presentation can be saved in cdf format and can be presented using any browser in which CDF Player is installed Quick tips for Inpatients ! Create Slide Show File -> New -> Slide Show Open Slide Show Palette Palettes -> Slide Show Run Slide Show View Environment -> SlideShow Run in Full Screen Mode

webservices - Wolfram Cloud and reCAPTCHA

Is there a piece of Wolfram Language code out there that integrates reCAPTCHA (for human verification) into the Wolfram Cloud environment? This is what the reCaptcha service looks like (in Swedish): When you press the checkbox, metrics such as how long it took to mark the checkbox and how the mouse moved are sent to Google's server. This data is used to determine whether it was a human that marked the checkbox. It is not the same as the generic term "captcha".

differential equations - Defining a non-linear coordinate transformation involving DSolve

This question is related to an other question I've asked yesterday here . Havin not received a reply, I believe there isn't a ready-made function for non-linear coordinate transformation in order to obtain a normal form of an affine state system as in nonlinear control theory. I wish to address here some of the problems in implementing such a function. We would have the following non-linear control system: x' = f(x) + g(x) u y = h(x) First step As is well explained in this book Isidori , what is firstly necessary is to define a nonlinear coordinate transformation. In order to do this, we ought to firstly define a Lie derivative. I've implemented the following code to do this: QDim[a_, b_] := TrueQ[Dimensions[a] == Dimensions[b]] Lie[l_, f_, x_] := With[{}, If[QDim[f, x] == True, dl = Table[D[l, x[[i]]], {i, 1, Dimensions[x][[1]]}]; Sum[dl[[i]] f[[i]], {i, 1, Dimensions[x][[1]]}], Print["Error: dimensional mismatch."]]] DLie[l_, f_, x_, k_] := With[

plotting - Animate ParametricPlot3D for two different parametric equations

I want to write a computer program to generate 3-D plots (with x-y-t axes) to demonstrate the evolution of the two missiles in the following case: i. $x = 100t$, $y = 80t-16t^2$ for $0\leq t\leq 5$ for the incoming missile ii. $x= 500-200(t-2)$, $y = 80(t-2)-16(t-2)^2$ for $2\leq t\leq 7$ for the interceptor missile I coded it like this: ParametricPlot3D[{{100 t, 80 t - 16 t^2, 60 t}, {500 - 200 (s - 2), 80 (s - 2) - 16 (s - 2)^2, s}}, {t, 0, 5}, {s, 2, 7}] The problem is that How can I animate this two functions in 3D form (with x-y-t axes)? I want it to display two moving objects,so i can see if there's any collision.By the way,please help me to correct the code if I'm wrong.

evaluation - Does function call via @ ignore HoldFirst attribute?

I was trying to test whether using func[x,y] is the same as func[#,y]&@x : SetAttributes[test, HoldFirst] test[1 + 2 - 3, 5 - 5] test[#, 5 - 5] &@(1 + 2 - 3) test[1 + 2 - 3, 0] test[0, 0] So apparently, HoldFirst is ignored in the second case. Why is it so? Is the argument of @ evaluated before actually being passed? Answer andre has given a good answer. I will add a little further clarification. The function f[#]& is a different function than the function f and does not have its attributes. It is the short form of Function[x, f[x]] Assume Clear[f]; SetAttributes[f, HoldFirst] has been evaluated. Even though f is HoldFirst , both Function[x, f[x]][1 + 1] and Function[x, f[x]] @ (1 + 1) give f[2] When the pure function is also given the attribute HoldFirst , evaluation goes as you expect. Function[x, f[x], HoldFirst] @ (1 + 1) f[1 + 1]

finite element method - Solving Navier-Stokes equations for a steady-state compressible viscous flow in a 2D axisymmetric step

Note: you may apply or follow the edits on the code here in this GitHub Gist I'm trying to follow this post to solve Navier-Stokes equations for a compressible viscous flow in a 2D axisymmetric step. The geometry is : lc = 0.03; rc = 0.01; xp = 0.01; c = 0.005; rp = rc - c; lp = lc - xp; Subscript[T, 0] = 300; Subscript[\[Eta], 0] = 1.846*10^-5; Subscript[P, 1] = 6*10^5 ; Subscript[P, 0] = 10^5; Subscript[c, P] = 1004.9; Subscript[c, \[Nu]] = 717.8; Subscript[R, 0] = Subscript[c, P] - Subscript[c, \[Nu]]; \[CapitalOmega] = RegionDifference[ Rectangle[{0, 0}, {lc, rc}], Rectangle[{xp, 0}, {xp + lp, rp}]]; And meshing: Needs["NDSolve`FEM`"]; mesh = ToElementMesh[\[CapitalOmega], "MaxBoundaryCellMeasure" -> 0.00001, MaxCellMeasure -> {"Length" -> 0.0008}, "MeshElementConstraint" -> 20, MeshQualityGoal -> "Maximal"][ "Wireframe"] Where the model is axisymmetric around the x axis, the governing

plotting - how to plot the dispersion curve in first brillouin zone?

Those who know some solid state physics should know what the first brillouin zone is. How do I plot the dispersion relation in the 1st brillouin zone so that the curves can 'fold back'? For instance the code I have now is: OmegaP=8; epM[Omega_]:=1-OmegaP^2/Omega^^2; epD[Omega_]:= 1; e2alpha[Omega_]:=(epM[Omega]-epD[Omega])/(epM[Omega]+epD[Omega]) ContourPlot[{Exp[0.5*ak]==e2alpha[Omega],Exp[0.5*ak]==-e2alpha[Omega]},{ak,0,0.5},{Omega,0.05,8}]; I have rescaled the axis so that 1st Brillouin zone is 0-0.5, 2nd Brillouin zone is 0.5-1.0, and so on. The above above only generates the lower part of the whole dispersion relation. If I extend the range of ak to, say {ak,0,2.5} , then I get a lot more missing curves, but how to I 'fold' the part from 0.5 to 2.5 into the 1st Brillouin zone? PS: Please click here (last page, Fig.16) for an image showing what the '1st Brillouin zone' is. Answer The "folding in" of the reduced zone scheme is achieved by tra

table - Undocumented TableView function

Seeking information on undocumented function TableView. Exempli gratia , tab=Table[Table[If[Mod[i,Prime[j]]==1,1,0],{i,1,10}],{j,1,2}] tab//TableView $ $ {{1, 0, 1, 0, 1, 0, 1, 0, 1, 0}, {1, 0, 0, 1, 0, 0, 1, 0, 0, 1}}

image processing - Improve TextRecognize[] on numbers

I have an image contain only numbers, and TextRecognize fail to recognize some numbers: img= ; TextRecognize[img] (*826718*) The documentation says that "The quality of recognized text may improve by enlarging the image", but no luck on this example TextRecognize[ImageResize[img, Scaled[2]]] (*826718*) also tried different language, also no help TextRecognize[ImageResize[img, Scaled[2]], Language -> "French"] (*826718*) I also tried Walfram|Alpha, it also gave the same results as Mathematica: Are there some ways to solve the problem? Answer TextRecognize seems to be a work in progress, consider the following Rasterize[Graphics[Text[Style["3", 100]]]] // TextRecognize Rasterize[Graphics[Text[Style["a", 100]]]] // TextRecognize Rasterize[Graphics[Text[Style["123", 100]]]] // TextRecognize Rasterize[Graphics[Text[Style["1234", 100]]]] // TextRecognize Rasterize[Graphics[Text[Style["hello", 100]]]] // TextRecognize

Plotting two functions in one graph, with different value ranges

Plot[2x, {x,0,4}]; Plot[x^2, {x,4,8}]; How do I merge these two graphs into one? Answer Here's a way that may serve you also for other purposes: p[x_, left_, right_] := HeavisideTheta[x - left] HeavisideTheta[right - x] Plot[{2 x p[x, 0, 4], x^2 p[x, 4, 8]}, {x, 0, 8}] Another example: tab = Table[x^(1/n) p[x, n, n + 1], {n, 1, 10}]; Plot[tab, {x, 0, 8}, PlotStyle -> Thick]

bugs - Cannot get nearest neighbors from Nearest

Bug introduced in 10 and fixed in 10.3 There still persists the failure to return 30.72 (another bug according to Daniel Lichtblau but from the developer's point of view it may be a feature). Nearest[data,x,{n,r}] is supposed to give the n or fewer nearest neighbors to x that are within radius r of x . Example: Nearest[{1, 2, 3, 4, 5, 6}, 3.6, {10, 1.5}] (* gives {4, 3, 5} *) Nearest[{1, 2, 3, 4, 5, 6}, 3, {10, 2}] (* gives {3, 2, 4, 1, 5}, so it's inclusive *) But this is not working for my data: v = {10.38,17.77,21.25,20.38,14.34,15.7,19.98,20.83,21.82,24.04,23.24,17.89,24.8,23.95,22.61,27.54,20.13,20.68,22.15,14.36,15.71,12.44,14.26,23.04,21.38,16.4,21.53,20.25,25.27,15.05,25.11,18.7,23.98,26.47,17.88,21.59,21.72,18.42,25.2,20.82,21.58,21.35,24.81,20.28,21.81,17.6,16.84,18.66,14.63,22.3,21.6,16.34,18.24,18.7,22.02,18.75,18.57,21.59,19.31,11.79,14.88,20.98,22.15,13.86,23.84,23.94,21.01,19.04,17.33,16.49,21.31,14.64,24.52,15.79,16.52,19.65,10.94,16.15,23.97,18,20.97,15.86

plotting - automatic processing of numerical results in `Plot`

First I want to solve an equation $F(x,y)=0$ for $y$ by supplying a value of $x$. (suppose obtaining the analytic form of $y(x)$ is too difficult) Then I want to plot root $y$ (numerically calculated) as a function of $x$ by using the following: Plot[y /. FindRoot[.../.{x->x0},{y, 0.2}],{x0, 0, 1}] and I got something like the following I omit ... here since it is terribly complicated. Update : F(x,y) is of the form $\sum_{n,m}a_{m,n}x^ny^m$ and $m$ and $n$ can be as high as 19, which basically makes Solve impractical. The result is satisfying except a small number of points near 1.0. Setting another initial starting value of $y$ in FindRoot might be an option but it is very tedious and often I cannot find a value of $y$ that fits the whole range of $x$. My question is: suppose I stick with the initial value of $y$, is there a way to just eliminate that anomaly point after I Plot the numerical results? Or is there a better way to deal with this kind of numerical problem in gener

import - Why does ImageData need four times more memory?

Background We are dealing with very large TIFF image files that are imported, processed and exported using Mathematica 9.0.1. Many of our algorithms only work with an array representation of the images, where pixel values are stored in the form of nested lists. For this, we are currently using the function ImageData to directly convert images into an array representation. Problem I am struggeling with the increase in memory usage when converting images to the array representation using ImageData . What I am wondering about is that when I convert images with ImageData the allocated memory is four times higher compared to having only the images in memory. Example Consider the following example: I first create a stack of 50 images that is exported to a TIFF file. Then I import the file again with the option "Data" . I check the memory before and after Import . Compared to importing the file without any option, this takes up four times more memory. Export["testFile.tif&qu

equation solving - Why does the integral not converge?

I am trying to solve for $\Omega$ this nonlinear integral equation: $$1+\dfrac{z}{k^2}-\dfrac{z^2}{K_2(z)} \dfrac{\Omega}{k^3} \displaystyle\int_{1}^{\infty} \gamma^2\, \text{ArcTanh} \left(\sqrt{\frac{\gamma ^2-1}{\gamma ^2}} \dfrac{k}{\Omega}\right)\, e^{-\text{z$\gamma $}} \, d\gamma=0$$ where $K_2(z)$ is the modified Bessel function of the second kind, $\Omega$ and $k$ are reals, $z> 0$. The Mathematica input is: f[Ω_?NumericQ, k_?NumericQ, z_?NumericQ] := 1 + (z/k^2) - (z^2/BesselK[2, z]) (Ω/k^3) NIntegrate[γ^2 (ArcTanh[Sqrt[(γ^2 - 1)/γ^2] (k/Ω)]) Exp[-z γ], {γ, 1, Infinity}, MaxRecursion -> 500] The solutions $\Omega (k, z)$ are given by : W[k_,z_] := Re[Ω /. FindRoot[f[Ω, k, z], {Ω, .895}]] I need to calculate all solutions $\Omega(k, z)$ for $k = 3$ when $0.1 My problem is that I get accurate solutions for $1 For example: Block[{k=3},Table[{z, W[k,z]},{z,{100,10,5,1,0.5,0.1}}]] {{100, 1.139642672363642}, {10, 1.96313715768855}, {5, 2.393983432376982}, {1, 2.90563

searching - Find the position of a string element, knowing only part of the string

I have a list, e.g.: list = {1, 2, 3, "Element 1", 4, 5, "Element 2", "Something else 1", "etcetera"} Now, I want all elements starting with "El" . Using Position , I could find the position of e.g. "Element 1" : Position[list, "Element 1"] But, I would like to know the positions of "Element 1" and "Element 2" , as both start with "El" . So, I would like to have something like Position[list, "El"_] I just can't get something like this to work. Thanks for any help. Answer But I would like to know the positions of "Element 1" and "Element 2" ... You can still use Position[] ; things are a little more elaborate, though, due to the strings: Position[list, s_String /; StringMatchQ[s, "El*"]] {{4}, {7}} Extract[list, %] {"Element 1", "Element 2"}

graphs and networks - Efficient solution for a discrete assignment problem with pairwise costs

Let's consider a simple graph with N vertices, and a corresponding set of N items. The goal of the problem is to assign every item to a vertex on the graph so that sum of a per-edge (that is item-pairwise) cost function over vertices is minimized. This problem clearly has a discrete search space of N! candidates. Exhaustive search of the optimal solution becomes practically impossible with very low values of N, and I'd expect that there would be methods to put more efficient algorithms at work on this problem. A brute-force toy attempt on the problem is presented below. Here g is the graph, i are the item values, and w is an extremely simple pairwise cost function: With[{ g = GridGraph[{3, 3}], i = {1, 4, 4, 9, 9, 16, 16, 32, 64}, w = Apply[Abs@*Subtract]}, First@TakeSmallestBy[ SetProperty[g, VertexLabels -> MapThread[Rule, {Range@VertexCount@g, #}]] & /@ Permutations@i, Function[g, Total[w[PropertyValue[{g, #}, VertexLabels] & /@ #] &a

pattern matching - Variable scoping problem when mapping over delayed replacement

This is something that got me curious while I was playing around with Mathematica . Consider the following (contrived) example: Clear[x,y]; {{1, 2}, {3, 4}, {5, 6}} /. {x_, y_} :> {x, #} & /@ {y, y^2, y^3} (* {{{1, 2}, {3, 4}, {5, 6}}, {{1, 4}, {3, 16}, {5, 36}}, {{1, 8}, {3, 64}, {5, 216}}} *) No problems there. Now suppose y has been defined globally, so y = 10; {{1, 2}, {3, 4}, {5, 6}} /. {x_, y_} :> {x, #} & /@ {y, y^2, y^3} (* {{{1, 10}, {3, 10}, {5, 10}}, {{1, 100}, {3, 100}, {5, 100}}, {{1, 1000}, {3, 1000}, {5, 1000}}} *) The problem is that y gets evaluated to 10 before getting mapped to the delayed rule's rhs. I realised I could do {{1, 2}, {3, 4}, {5, 6}} /. {x_, y_} :> {x, ReleaseHold[#]} & /@ Thread[Hold[{y, y^2, y^3}]] (* {{{1, 2}, {3, 4}, {5, 6}}, {{1, 4}, {3, 16}, {5, 36}}, {{1, 8}, {3, 64}, {5, 216}}} *) which looks kinda complicated, but works. So, is there any better way to "shield" the y in the list being mapped over, from the g