Showing posts with label Mathematica. Show all posts
Showing posts with label Mathematica. Show all posts

Dec 22, 2013

On the Mutual Languages of States

Languages lie beyond the boarder of nations. Normally, if two states speak the same language, cultural connection is implied. The notable example is the US. and UK, both are built by English speaking people. In the real world, one state may speak various languages. Each language is spoken by some portion of the population. To model the language similarity between two states, we introduce a function $\phi(state1, state2)$ such that: 1) $\phi = 1$ if state1 = state2; 2) $\phi(state1, state2) = \phi(state2, state1)$; 3) if no mutual language spoken between state1 and state2, then $\phi(state1, state2) = 0$.

A very intuitive definition of $\phi$ could be the fraction of people speaking any mutual language of two states. Let $p_1, p_2$ be the population of the two states, respectively. $f_i^{(1)}, f_i^{(2)}$ be the population fractions of the  $i$-th mutual language. Then
\[ \phi(state1, state2) := \frac{p_1 \sum_i f_i^{(1)} + p_2 \sum_i f_i^{(2)}}{p_1+p_2}.
\] For example, the mutual native languages spoken in the US and UK are English and Angloromani. In the US, the fractions of these two languages are: 89% and 0.043%. In the UK, the numbers are 98% and 0.16%. Taken into account the population of the two states, 308.476 millions and 61.899 millions, the similarity function $phi$ yields 0.9085.

Another interesting definition of $\phi$ (but 1). is not satified) would be the probability of two men from the two states respectively speak the same tongue.  It can be shown the probability is \[
\phi(state1, state2) := \sum_i f_i^{(1)} f_i^{(2)}.
\] In the case of US and UK, $\phi$ yields 0.8722.

Asia

With this definition, let's explore the international language connection. Fig. 1 shows the adjacent matrix of Asia states. Fig. 2 shows the graph representation of it wit the edge colored according to the language similarity.

Fig. 1, the language connection among Asian states: the adjacent matrix.

Fig. 2, the language connections among Asian states: the graph (We have removed the self-loops). The edges are colored according to their values.
Fig. 3 again shows the language connections among Asian states. But the edges are dyed by the type of the primary mutual language shared by the two vertex states. The red edges in Fig. 3 represent English. Compare with Fig. 2, English is not very popular in Asia but still widely used.
Fig. 3, the language connections among Asian states: the graph edges are colored according to the type of primary mutual languages shared by each pair of states. 
Notice that some states enjoy closer language relation with their neighbours than others. We can define a language capacity of a state within some group: \[
LC(s) := \frac{1}{|G|-1}\sum_{s' \in G, s'\ne s} \phi(s, s')
\] Fig. 4 shows the language capcities of Asian states.

Fig. 4, the Language Capacity of the states in Asia. The two countries that have the largest two language capacities are China and Singapore.
We can gather further information from the language relation graph. We can split the vertices into graph communities. Using the FindGraphCommunities function in Mathematica, we can identify 7 communities. The first community is Southeast Asia including China, Singapore, Brunei etc. The most popular language is Chinese, the most widely used language is English. The second and third communities are the crossroad countries including Russia, the central asia stans, Iran, Iraq, Syria, Isreal etc. The most popular and most widely used languages include Russia and Kurdish. The fourth zone is south asia, including India, Pakistan, Bangladesh Nepal Bhutan and Myanmar. The most popular language is Hindi and the most widely used language is Tibetan. The fifth region is the Arabic states. The most popular and most widely used language is Arabic. Next region is Japan and the Koreas. The most popular language is Japanese and the most widely used language is Korean. The last region is Turkey and Uzbekistan. Both countries speak Turkish.


Fig. 5, the language communities of the Asia states. The edges are colored by the communities.

Europe

Europe is another continent with flourish civilizations. The average language capacity in Europe is higher than in Asia.


The languages spoken in Asia may be very diverse, it is not true in other continent. Europe, for example,  is dominated by French, German, Italian, English and Romanian etc.
Finally, the world is dominated by four languages: English, French, Spanish and Portuguese - they may not be the most popular ones - they are widely used in communicating with other states. Apparently, a language beyond the boundary of a country is a proof of the cultural influence. As we know, the four dominant "foreign languages" are the results of the colonism.

Dec 7, 2013

Frenet Tube


Given a 3D curve $\vec{r}(t) = (x(t), y(t), z(t)), t \in I=[0,1]$,  Tubify[ {x, y, z},  a] turns it into a tube around the curve. The Mathematica Tube[] has similar function (works on lines). The essential thing is to find the unit circle in the normal plane for each point (Shown in Fig. 1). The tangential vector is simply $\vec{\alpha} = \vec{r}'/|\vec{r}'|$, where $\vec{r}' = \mathrm{d} \vec{r}/\mathrm{d}t$. The normal vector is $\vec \beta = \vec{\alpha}' / |\vec{\alpha}'| = \frac{\vec{r}'' |\vec{r}'|^2 - \vec{r}' \vec{r}'\cdot\vec{r}''}{|\vec{r}'| \sqrt{|\vec{r}''|^2 |\vec{r}'|^2 - (\vec{r}'\cdot\vec{r}'')^2} }$. The binormal vector is $\vec{\gamma} = \vec{\alpha} \times \vec{\beta} =  \frac{\vec{r}' \times \vec{r}''  }{\sqrt{|\vec{r}''|^2 |\vec{r}'|^2 - (\vec{r}'\cdot\vec{r}'')^2} } $.

Fig. 1: black, $\vec{\alpha}$, green $\vec{\beta}$, blue $\vec{\gamma}$. 

Tubify[{x_, y_, z_}, r_] :=
 ({x[#1], y[#1], z[#1]}
    + (r[#1] Cos[#2])/
      Norm[(Norm[{x'[#1], y'[#1], z'[#1]}]^2 {x''[#1], y''[#1], 
           z''[#1]} - {x'[#1], y'[#1], z'[#1]}.{x''[#1], y''[#1], 
            z''[#1]} {x'[#1], y'[#1], z'[#1]})] (Norm[{x'[#1], y'[#1],
            z'[#1]}]^2 {x''[#1], y''[#1], 
         z''[#1]} - {x'[#1], y'[#1], z'[#1]}.{x''[#1], y''[#1], 
          z''[#1]} {x'[#1], y'[#1], z'[#1]})
    + (r [#1] Sin[#2] )/
      Norm[{-z'[#1]*y''[#1] + z''[#1]*y'[#1], 
        z'[#1]*x''[#1] - z''[#1]*x'[#1], -y'[#1]*x''[#1] + 
         x'[#1]*y''[#1]}] {-z'[#1]*y''[#1] + z''[#1]*y'[#1], 
      z'[#1]*x''[#1] - z''[#1]*x'[#1], -y'[#1]*x''[#1] + 
       x'[#1]*y''[#1]}) &

All the arguments are functions.

Example 1:
Show[ParametricPlot3D[
  Tubify[{(2 + Cos[3 #]) Cos[ #] &, (2 + Cos[3 #]) Sin[#] &, 
     Sin[3 #] &}, 0.25 &][u, v],
  {u, 0, 2 Pi}, {v, 0, 2 Pi}, PlotStyle -> Opacity[0.2],
  Mesh -> None, Boxed -> True, Axes -> False, BoxRatios -> Automatic, 
  PerformanceGoal -> "Quality", PlotPoints -> 80, MaxRecursion -> 0, 
  PlotLabel -> 
   Style[TraditionalForm /@ {(2 + Cos[3 t]) Cos[ 
        t], (2 + Cos[3 t]) Sin[t], Sin[3 t]}, 20]
  ],
 ParametricPlot3D[{(2 + Cos[3 u]) Cos[ u], (2 + Cos[3 u]) Sin[u], 
   Sin[3 u]}, {u, 0, 2 Pi}, 
  PlotStyle -> Directive[{Red, Opacity[1], Thick}]], ImageSize -> 600
 ]

A similar code using Tube[]:

ParametricPlot3D[
  {(2 + Cos[3 u]) Cos[ u], (2 + Cos[3 u]) Sin[u], 
   Sin[3 u]}, {u, 0, 2 Pi}, PlotStyle -> Opacity[0.2],
  Mesh -> None, Boxed -> True, Axes -> False, BoxRatios -> Automatic, 
  PerformanceGoal -> "Quality", PlotPoints -> 80, MaxRecursion -> 0, 
  PlotLabel -> 
   Style[TraditionalForm /@ {(2 + Cos[3 t]) Cos[ 
        t], (2 + Cos[3 t]) Sin[t], Sin[3 t]}, 20]
  ]/.Line[pts_, rest___] :> Tube[pts, 0.1, rest]

Fig. 2

Example 2, trefoil knot:


Show[ParametricPlot3D[
  Tubify[{Sin[#] + 2 Sin[2 #] &, Cos[#] - 2 Cos[2 #] &, -Sin[3 #] &}, 
    0.25 &][u, v],
  {u, 0, 2 Pi}, {v, 0, 2 Pi}, PlotStyle -> Opacity[0.3],
  Mesh -> None, Boxed -> True, Axes -> False, BoxRatios -> Automatic, 
  PerformanceGoal -> "Quality", PlotPoints -> 100, MaxRecursion -> 0, 
  PlotLabel -> Style[TraditionalForm /@ trefoil[t], 20]
  ], ImageSize -> 600
 ]

Fig. 3 Trefoil Knot

Similarly, we can also define the Frenet's Ribbon:
Ribbonize[{x_, y_, z_}, r_] :=
 ({x[#1], y[#1], z[#1]}
    + r [#1] #2 /
      Norm[{-z'[#1]*y''[#1] + z''[#1]*y'[#1], 
        z'[#1]*x''[#1] - z''[#1]*x'[#1], -y'[#1]*x''[#1] + 
         x'[#1]*y''[#1]}] {-z'[#1]*y''[#1] + z''[#1]*y'[#1], 
      z'[#1]*x''[#1] - z''[#1]*x'[#1], -y'[#1]*x''[#1] + 
       x'[#1]*y''[#1]}) &

Fig. 4, Trefoil Ribbon


Fig. 5, Helix Ribbon

more examples:

Fig. 6, A Seashell Surface
Fig. 7, Double Helix (cf. DNA)
Fig. 8, A Figure-eight Knot

ParametricPlot3D[
 Tubify[{#/Pi Cos[#] &, #/Pi Sin[#] &, #/Pi &}, Sin[#/4] &][u, u*80],
 {u, 0, 4 Pi}, PlotStyle -> Directive[{Red, Thick, Opacity[0.8]}],
 Mesh -> None, Boxed -> True, Axes -> False, BoxRatios -> Automatic,
 PlotRange -> All, PerformanceGoal -> "Quality", MaxRecursion -> 0, 
 PlotPoints -> 2560, PlotLabel -> Style["", 20], 
 ImageSize -> 500
 ]
Fig. 9, a coilfied 3D curve

Nov 11, 2013

the graph of the international trade relations

Mathematica has provided excellent database of the countries in the function CountryData. I am particularly interested in the trade relations. CountryData[country, "ImportPartners"] and CountryData[country, "ExportPartners"] produce the leading import and export partners of the country. We can construct a graph based on the global trade relations.

trade = 
  Tally[Sort /@ 
    Flatten[Tuples[{{#}, 
         Select[Flatten[{CountryData[#, 
              "ImportPartners"] /. _Missing -> Null, 
            CountryData[#, "ExportPartners"] /. _Missing -> Null}], 
          MemberQ[CountryData[], #] &]}] & /@ CountryData[], 1]];

tradegraph = 
 Graph[Thread[
   trade[[All, 1, 1]] \[UndirectedEdge] 
    trade[[All, 1, 2]]], 
  EdgeWeight -> trade[[All, 2]]]

Fig. 1, the graph of the international trade relations.

We can make it more accessible by coloring the countries.

by continent:

Row[
 {Graph[Thread[
    trade[[All, 1, 1]] \[UndirectedEdge] 
     trade[[All, 1, 2]]], 
   EdgeWeight -> trade[[All, 2]],
   VertexShapeFunction -> ({ColorData[3, 
        First@First@
          Position[CountryData["Continents"], 
           CountryData[#2, "Continent"] ]], Opacity[.8], 
       EdgeForm[{Thin, LightGray}], Disk[#, .08], Darker@Red, 
       Opacity[1.0], 
       Text[Style[CountryData[#2, "CountryCode"], Bold], #1] } &), 
   ImageSize -> 800],
  Graphics[
   MapIndexed[{ColorData[3, First@#2], 
      Text[Style[
        StringJoin[
         Riffle[StringCases[#1, RegularExpression["[A-Z][a-z]+"]], 
          " "]], Bold, 22], {40, First@#2*60}]} &, 
    CountryData["Continents"]],
   ImageSize -> {Automatic, 300}]
  }]

Fig. 2,  (Colored by Continent) the graph of the international trade relations.

by trade values, 


Fig. 3,  (Colored by Continent) the graph of the international trade relations.

Fig. 4, colored by the import and export amount separated.

We can explore some further graph structures inside the international trade relations. What we focus here is the the graph communities. The communities characterize the clustering property of a graph. Apply the Mathematica function FindGraphCommunities[], we obtain several trade clusters. The leading 3 communities correspond to the three known trade zones: the Transatlantic Trade Zone including US. and the European Union; the Asia-Pacific Trade Zone including China, Japan, Korea Australia and the ASEANs; the South-America Trade Zone including Brazil, Argentina etc. Our conclusion is further justified if we color the countries inside the trade zone by their continent (Fig. 6).
Fig. 5, the graph communities of the international trade relations.
Fig. 6, trade communities and the geographic distribution. 

Sep 12, 2013

Tangram

Fig. 1 A standard tangram with rainbow colors: (red, orange, yellow, green, cyan, blue, purple)
The tangram (Chinese: 七巧板 - seven boards) is a set of boards (shown above) that can be put together to form specific shapes.  For example, different shapes of cat:
Fig. 2 several shapes of cat

Mathematica codes to draw the tangram shown in Fig. 1:
Graphics[{EdgeForm[{Darker@Gray, Thick, Opacity[0.4]}],
  Orange, Polygon[{{0, 0}, {0, 1}, {1/2, 1/2}}],
  Green, Polygon[{{0, 1}, {1, 1}, {1/2, 1/2}}],
  Red, Polygon[{{1, 1}, {1, 1/2}, {3/4, 3/4}}],
  Yellow, Polygon[{{3/4, 3/4}, {1, 1/2}, {3/4, 1/4}, {1/2, 1/2}}],
  Blue, Polygon[{{1/2, 1/2}, {3/4, 1/4}, {1/4, 1/4}}],
  Cyan, Polygon[{{1/4, 1/4}, {3/4, 1/4}, {1/2, 0}, {0, 0}}],
  Purple, Polygon[{{1/2, 0}, {1, 1/2}, {1, 0}}]
  }]

A Mathematica app to play with tangram:

Framed@DynamicModule[{O1 = {1/6, 1/2}, O2 = {1/2, 5/6}, 
   O3 = {11/12, 3/4}, O4 = {3/4, 1/2}, O5 = {1/2, 1/3}, 
   O6 = {3/8, 1/8}, O7 = {5/6, 1/6}, p1 = {-.25, -.25}, 
   p2 = {0, -.25}, p3 = {0.25, -.25}, p4 = {.5, -.25}, 
   p5 = {.75, -.25}, p6 = {1, -.25}, p7 = {1.25, -.25}, v1, v2, v3, 
   v4, v5, v6, v7, r = 0.08, v0 = {1, 0}, 
   g = Graphics[{Gray, PointSize[Medium], Point[{0, 0}]}], 
   gray = False},
  
  {v1, v2, v3, v4, v5, v6, 
    v7} = ({r, 0} + #) & /@ {p1, p2, p3, p4, p5, p6, p7};
  Column[{Row[{Text["check here to paint gray:"], 
      Checkbox[Dynamic[gray]]}], 
    LocatorPane[
     Dynamic[{v1, v2, v3, v4, v5, v6, v7, O1, O2, O3, O4, O5, O6, 
       O7}], Dynamic[
      Graphics[{EdgeForm[], If[gray, Gray, Orange], 
        Polygon[RotationTransform[{v0, Normalize[v1 - p1]}, 
           O1] /@ {O1 - {1/6, 1/2} + {0, 0}, O1 - {1/6, 1/2} + {0, 1},
            O1 - {1/6, 1/2} + {1/2, 1/2}}], If[gray, Gray, Green], 
        Polygon[RotationTransform[{v0, Normalize[v2 - p2]}, 
           O2] /@ {O2 - {1/2, 5/6} + {0, 1}, O2 - {1/2, 5/6} + {1, 1},
            O2 - {1/2, 5/6} + {1/2, 1/2}}], If[gray, Gray, Red], 
        Polygon[RotationTransform[{v0, Normalize[v3 - p3]}, 
           O3] /@ {O3 - {11/12, 3/4} + {1, 1}, 
           O3 - {11/12, 3/4} + {1, 1/2}, 
           O3 - {11/12, 3/4} + {3/4, 3/4}}], If[gray, Gray, Yellow], 
        Polygon[RotationTransform[{v0, Normalize[v4 - p4]}, 
           O4] /@ {O4 - {3/4, 1/2} + {3/4, 3/4}, 
           O4 - {3/4, 1/2} + {1, 1/2}, O4 - {3/4, 1/2} + {3/4, 1/4}, 
           O4 - {3/4, 1/2} + {1/2, 1/2}}], If[gray, Gray, Cyan], 
        Polygon[RotationTransform[{v0, Normalize[v5 - p5]}, 
           O5] /@ {O5 - {1/2, 1/3} + {1/2, 1/2}, 
           O5 - {1/2, 1/3} + {3/4, 1/4}, 
           O5 - {1/2, 1/3} + {1/4, 1/4}}], If[gray, Gray, Blue], 
        Polygon[RotationTransform[{v0, Normalize[v6 - p6]}, 
           O6] /@ {O6 - {3/8, 1/8} + {1/4, 1/4}, 
           O6 - {3/8, 1/8} + {3/4, 1/4}, O6 - {3/8, 1/8} + {1/2, 0}, 
           O6 - {3/8, 1/8} + {0, 0}}], If[gray, Gray, Purple], 
        Polygon[RotationTransform[{v0, Normalize[v7 - p7]}, 
           O7] /@ {O7 - {5/6, 1/6} + {1/2, 0}, 
           O7 - {5/6, 1/6} + {1, 1/2}, O7 - {5/6, 1/6} + {1, 0}}], 
        Opacity[0.6], Gray, PointSize[0.02], 
        Point[{O1, O2, O3, O4, O5, O6, O7}], Orange, Disk[p1, r], 
        Green, Disk[p2, r], Red, Disk[p3, r], Yellow, Disk[p4, r], 
        Cyan, Disk[p5, r], Blue, Disk[p6, r], Purple, Disk[p7, r], 
        Black, Arrowheads[Medium], 
        Arrow /@ {{p1, p1 + r Normalize[v1 - p1]}, {p2, 
           p2 + r Normalize[v2 - p2]}, {p3, 
           p3 + r Normalize[v3 - p3]}, {p4, 
           p4 + r Normalize[v4 - p4]}, {p5, 
           p5 + r Normalize[v5 - p5]}, {p6, 
           p6 + r Normalize[v6 - p6]}, {p7, 
           p7 + r Normalize[v7 - p7]}}}, 
       PlotRange -> {{-1/2, 3/2}, {-1/2, 3/2}}, ImageSize -> 600]], 
     Appearance -> None ]}]]

copy codes from gist: https://gist.github.com/anonymous/6564099
Fig. 3 use the app to form some shape

Jul 11, 2013

Connect to Remote MathKernels

Suppose you have a powerful workstation in your office and a laptop at home. Suppose you prefer to working in your comfortable apartment. But now you have some work with Mathematica to do. Your laptop is convenient to visualize data but incapable of doing the heavy computing tasks. You may think of ssh X11 forwarding:
ssh -X username@hostname.domain
The problem is, with limited network bandwith, X11 forwarding is simply not stable to do a typical Mathematica task. The good news is Mathematica allows you to do use remote MathKernel(s) but local  FrontEnd.

1. Connect to Remote MathKernel(s);




Unfortunately, the default (Basic Options or Advanced Options) configuration does not always work for typical family networks (router, firewall, NAT etc). This error message may pop up:


The problem is that your laptop's IP address is not known. For example, you may hide behind a wireless router and a firewall. The way to solve the problem is to tell MathKernels your port and IP address (your machine may not have one!). If you do have an IP,  replace `ipadress` with it.

You can get your public IP simply by googling "my ip". Of course, your public IP may not be the IP of your machine.

If you are using a dynamic IP address, remember your IP usually changes within a few days.

In the presence of firewall and router, your IP address may not exist. Instead, your router has a public IP (again, could be dynamic or static), say, 12.34.5.6 And your machine was assigned to a local IP address, say, 192.168.1.2 The remote MathKernel(s) will have problem talking to your local machine. What you need is a secure link.

VPN

The first choice is a VPN. With a VPN, you can happily work with the abstract address such as foo.vpn.thirdparty.com , regardless of the actual change of your IP address.


The mathssh is a Mathematica implementation of ssh secure shell built on Java. One can also use ssh directly: `java` -jar `mathssh` -> ssh


SSH Reverse Tunneling 

SSH reverse tunneling (aka. SSH tunnel) is secure method. Interested users can check out reference [3] and the man page of ssh -R option. A Free Package was provided by Sascha Kratky. Both Unix-like (Mac & Linux) and Windows versions are provided (See Ref [5]).

Port Forwarding

Port forwarding is another unsafe approach. Namely you can forward your router's port 22 (the very port that ssh uses), to your local IP say 192.168.1.2.

Warning: port forwarding is not recommended because it is not safe to expose your ssh port. Your computer will be under frequent attack (just like a public server). To improve security, you may want to disable the root login on the ssh server.

vim /etc/ssh/sshd_config
# Add or uncomment this line:
PermitRootLogin no

Your public IP can be checked out simply by googling "my ip". Your local IP can be obtained from (replace wlan0 by eth0 if you are using wired connection):



Ref
[0]: How to | Connect to a Remote Kernel
[1]: Remote kernel for Mathematica via home router
[2]: FAQ Remote Kernels
[3]: Reverse SSH tunnel or connecting to computer behind NAT router
[4]: Manually Creating a Remote Kernel Connection
[5]: Remote Kernel Strategies by Sascha Kratky


2. BTW: Connect to Remote Parallel Kernels

Setting up remote Parallel Kernels is much easier, if the master kernel is local. The master and slaves can simply communicate through ssh.


Jun 28, 2013

The OutputStream in Mathematica

See this link .

Format:

$stream = OpenWrite["file.ext"]

Export[ $stream, expr, "format"]

Close[$stream]



vs. Export[$file.ext", expr, "format"]



$stream can be used to export multiple expressions and to append to exist files.

$stream = OpenAppend["file.ext"]


May 4, 2013

Random Walk

I use Mathematica to play with 2D Brownian Motion a little bit. According to Smoluchowski and Einstein's theory, Brownian Motion is described by random walk (RW). There are several methods to generate random walk in Mathematica, such as the use of Accumulate + RandomChoice, which is an example in the documentation:
Graphics[
  Line[
    Accumulate[
      RandomChoice[{{-1, 0}, {1, 0}, {0, 1}, {0, -1}}, 1000]
    ]
  ]
]
Fig. 1,  Random walk on a 2D square lattice
Fig. 2, a random walk with 10^4 steps on the square lattice

Fig. 3, calculation of the fractal dimension using box counting method. $a$ is the size of the boxes, the vertical values counts the number of boxes that cover the random walk points. It gives a fractal dimension of ~1.602.
The mean distance after $N$ steps is proportional to $\sqrt{N}$. This is known as the Einstein's law. Fig. 4 shows the test of this law.

Fig. 4, the mean distance. The nonlinear fit shows that the mean distance agree with the Einstein's law.


Similarly, we can work out a random walk on a triangular lattice.
Graphics[
  Line[Accumulate[
    RandomChoice[{{-1, 0}, {1, 0}, {1/2, Sqrt[3]/2}, {-1/2, 
       Sqrt[3]/2}, {1/2, -Sqrt[3]/2}, {-1/2, -Sqrt[3]/2}}, 1000]]]]

Fig. 5, Random walk on a 2D triangular lattice


In order to incorporate more sophisticated random walks, I use Nest and NestList.

 
 v = 1;
 rw[r_] := r + v * Through[{Cos,Sin}[RandomReal[{0,2*Pi}]]];
 urw[r_] := r + RandomReal[{0,2*v}] * Through[{Cos,Sin}[RandomReal[{0,2*Pi}]]];
 mbrw[r_] := r + RandomVariate[MaxwellDistribution[v]]* Through[{Cos,Sin}[RandomReal[{0,2*Pi}]]];

 n = 10^5;
 Show[
  Graphics[{Thickness[Large],Circle[{0, 0}, Sqrt[n]], 
            Arrow[{{0, 0}, Sqrt[n] {Cos[Pi/4], Sin[Pi/4]}}], 
            Text[Style[Sqrt[N], 20], 0.5 Sqrt[n] {Cos[Pi/4], Sin[Pi/4]}],
            PointSize[0.01], Point[{0, 0}],
            Red, Opacity[0.5], Thin, Line[NestList[rw,{0,0},n]],
            Green, Opacity[0.5], Thin, Line[NestList[urw,{0,0},n]],
            Yellow, Opacity[0.5],Thin, Line[NestList[mbrw,{0,0},n]]}
  ],
 ImageSize -> 400]
 

Fig. 6, Colorlines: Red, random walk; Green, uniform distribution random walk; Yellow, Maxwell-Boltzmann random walk.

Fig. 7, RW in a circle.

Fig. 8, RW in a circle, accelerated in radial direction.
Now one can study the quantitative behavior of Brownian Motion, following the historical path of Langevin, Smoluchowski and Einstein. I'll put that part aside.

The next question is what will happen if self-interactions are added. One of these self-interactions is self-avoid (SA). In SARW, the new walk try to avoid its history. Depending on how to interpret the "history", one may avoid either the "path" or the "points". I am more interesting in the points here, because points may represent the "system configurations" (system states/phase space points/field configuration/functionals), whereas the path represents another "field" such as the gauge field/gauge links etc. Avoiding points is the simplest self-avoid self-interaction.

To implement SA, the whole list has to be stored.

 greedymbsarw[rn_] := Module[{new , eps = 1, v = 1},

  (* WriteString["stdout", Length[rn], " "];*)
  new = rn[[-1]] + RandomVariate[MaxwellDistribution[v]]* Through[{Cos, Sin}[RandomReal[{0, 2 Pi}]]];
  While[
   Or @@ Table[Norm[new - pt] < eps, {pt, rn}],
   new = rn[[-1]] + 
      RandomVariate[MaxwellDistribution[v]]* Through[{Cos, Sin}[RandomReal[{0, 2 Pi}]]];
   ];
  Append[rn, new]
  ];

  lazymbsarw[rn_] := 
    Module[{new , eps = 1, v = 1},
      new = rn[[-1]] + RandomVariate[MaxwellDistribution[v]]* Through[{Cos, Sin}[RandomReal[{0, 2 Pi}]]];
      If[ And @@ Table[Norm[new - pt] > eps, {pt, rn}], Append[rn, new], rn]  
    ];

 n = 10^5;
 Show[
  Graphics[{Thickness[Large],Circle[{0, 0}, Sqrt[n]], 
            Arrow[{{0, 0}, Sqrt[n] {Cos[Pi/4], Sin[Pi/4]}}], 
            Text[Style[Sqrt[N], 20], 0.5 Sqrt[n] {Cos[Pi/4], Sin[Pi/4]}],
            PointSize[0.01], Point[{0, 0}],
            Red, Opacity[0.5], Thin, Line[NestList[mbrw,{0,0},n]],
            Green, Opacity[0.5], Thin, Line[Nest[lazymbsarw,{{0,0}},n]],
            Blue, Opacity[0.5], Thin, Line[Nest[greedymbsarw, {{0, 0}}, n]]}
  ],
 ImageSize -> 400]

In these two implementations, greedy version will loop till a new legal walk is found; whereas the lazy version stations if new legal walk if found at the first proposal. It is important to use Maxwell-Boltzmann (or other global ) RW to avoid dead-lock. Note that for the greedy version may take a very long time. One thus may also add a upper bound for trials.
Colorlines: Red, Maxwell-Boltzmann RW; Green, lazy SARW; Blue, greedy SARW.
Again, I'll leave the discussion of quantitative feature later.

We can play with more interactions.





update:
OK, I spent sometime to workout the random self-avoid walk, with respect to avoiding the path instead of the points. Essentially, one only need to compute if two segments cross. The function to do that is,

IntersectQ[x1_, y1_, x2_, y2_, x1p_, y1p_, x2p_, y2p_] := 
 Module[{j1, j2, det}, 
  j1 = (x1 + x2 - x1p - x2p)*(y1p - y2p) - (y1 + y2 - y1p - 
       y2p)*(x1p - x2p);
  j2 = (x1 + x2 - x1p - x2p)*(y1 - y2) - (y1 + y2 - y1p - y2p)*(x1 - 
       x2);
  det = (x1 - x2)*(y1p - y2p) - (y1 - y2)*(x1p - x2p);
  Abs[j1] < Abs[det] && Abs[j2] < Abs[det]]

As one may check, $t=J_1/\det, s=J_2/\det$ is the solution (parameters) of the intersection of two lines. The explicit coordinate is $(x, y) = (\frac{x_1+x_2}{2}+t\frac{x_2-x_1}{2}, \frac{y_1+y_2}{2}+t\frac{y_2-y_1}{2})$. If $-1 < t < 1, -1 < s < 1$ then the segments intersect.

lazyrsaw[rn_] := Module[{new, eps = 1, l}, 
   l = Length[rn];
   new = rn[[-1]] + 
     RandomVariate[MaxwellDistribution[1]]*
      Through[{Cos, Sin}[RandomReal[{0, 2 Pi}]]];
   If[And @@ 
     Table[! (IntersectQ[rn[[-1, 1]], rn[[-1, 2]], new[[1]], new[[2]], 
         rn[[i, 1]], rn[[i, 2]], rn[[i + 1, 1]], rn[[i + 1, 2]]]), {i,
        1, l - 1}], Append[rn, new], rn
    ]
 ];

random SAW with lazy algorithm
It's clear that the RW easily gets stuck in the corner, though possibility exist for it to jump out of it. It's better to take a look at the typical behavior of random SAW:

$\bar{N}$ is the average steps of the RSAWs.




.

Feb 14, 2013

The Secret Formula for Love


Mathematica code:

(* a xkcd-style graph converter from http://mathematica.stackexchange.com/questions/11350/xkcd-style-graphs 
by Simon Woods *)

xkcdStyle = {FontFamily -> "Comic Sans MS", 16};

xkcdLabel[{str_, {x1_, y1_}, {xo_, yo_}}] := 
  Module[{x2, y2}, x2 = x1 + xo; y2 = y1 + yo;
   {Inset[
     Style[str, xkcdStyle], {x2, y2}, {1.2 Sign[x1 - x2], 
      Sign[y1 - y2] Boole[x1 == x2]}], Thick, 
    BezierCurve[{{0.9 x1 + 0.1 x2, 0.9 y1 + 0.1 y2}, {x1, y2}, {x2, 
       y2}}]}];

xkcdRules = {EdgeForm[ef : Except[None]] :> 
    EdgeForm[Flatten@{ef, Thick, Black}], 
   Style[x_, st_] :> Style[x, xkcdStyle], 
   Pane[s_String] :> 
    Pane[Style[s, xkcdStyle]], {h_Hue, l_Line} :> {Thickness[0.02], 
     White, l, Thick, h, l}, 
   Grid[{{g_Graphics, s_String}}] :> Grid[{{g, Style[s, xkcdStyle]}}],
    Rule[PlotLabel, lab_] :> Rule[PlotLabel, Style[lab, xkcdStyle]]};

xkcdShow[p_] := 
 Show[p, AxesStyle -> Thick, LabelStyle -> xkcdStyle] /. xkcdRules

xkcdShow[Labeled[p_, rest__]] := 
 Labeled[Show[p, AxesStyle -> Thick, LabelStyle -> xkcdStyle], 
   rest] /. xkcdRules

(* Mathematica Version 8.0+ *)
xkcdDistort[p_] := 
  Module[{r, ix, iy, rand}, 
   r = ImagePad[Rasterize@p, 10, Padding -> White];
   {ix, iy} = 
    Table[RandomImage[{-1, 1}, ImageDimensions@r]~ImageConvolve~
      GaussianMatrix[10], {2}];
   ImagePad[
    ImageTransformation[
     r, # + 15 {ImageValue[ix, #], ImageValue[iy, #]} &, 
     DataRange -> Full], -5]];

(* for Mathematica version < 7 *)
xkcdDistort[p_] := 
 Module[{r, id, ix, iy, samplepoints, funcs, channels},
  r = ImagePad[Rasterize@p, 10, Padding -> White]; 
  id = Reverse@ImageDimensions[r];
  {ix, iy} = Table[ListInterpolation[ImageData[
      Image@RandomReal[{-1, 1}, id]~ImageConvolve~GaussianMatrix[10]]], {2}]; 
  samplepoints = Table[{x + 15 ix[x, y], y + 15 iy[x, y]}, {x, id[[1]]}, {y, id[[2]]}]; 
  funcs = ListInterpolation[ImageData@#] & /@ ColorSeparate[r]; 
  channels = Apply[#, samplepoints, {2}] & /@ funcs; 
  ImagePad[ColorCombine[Image /@ channels], -10]]


xkcdConvert[x_] := xkcdDistort[xkcdShow[x]]


(* the heart plot *)

xkcdConvert[
 ContourPlot[x^2 + (y - (x^2)^(1/3))^2 == 1, {x, -1, 1}, {y, -1, 2}, 
  ContourStyle -> {Red, Thick}, ContourShading -> Red, 
  AspectRatio -> 1, ImageSize -> 400, Axes -> True, Frame -> False, 
  AxesLabel -> Map[Text@Style[#, 20, Italic] &, {"x", "y"}], 
  AxesStyle -> FontSize -> 15, 
  Epilog -> {Inset[
     Style[x^2 + (y - (x^2)^(1/3))^2 == 1, Medium, Italic], 
     Offset[{0, 0}], Scaled[{-0.1, -5.5}], 1], 
    Inset[Style["xkcd-style graph", xkcdStyle, FontSize -> 12 ], 
     Scaled[{0.15, 0.05}], Automatic]}, 
  PlotLabel -> 
   Style["The Secret Formula for Love", xkcdStyle, 
    FontSize -> 24]]]

Sep 7, 2012

Faster Mathematica II: List Manipulation


The most important datatype in math is list. Lists in Mathematica are like set in mathematics. List is also used to hold data. In the case of large amount of data in input/output or intermediate stage, fast list operations certainly would speedup the code. One surprising fact is, seemly the same list operations may have different performance.

1. List Construction

1-0 Table is fast when creating block list

The most simple way to construct a list is use Table. Table[] is the suitable choice in generic cases. But making use of knowledge of the list, you can get better performance.
Table[], Array[], and function apply

1-1 use build-in functions

Some build-in functions support generation of a list. They are usually faster that explicit call them in a list constructor like Table[].  
Build-in random number generator is way faster.

1-2 use Reap-Sow, instead of Insert, Append, AppendTo, Prepend or PrependTo

1-3 Use NestList for recursable list  

Recursable list means you can get the n+1-th from the n-th. 


2. List Traversal

2-1. Listable attribute

Many function as Listable attribute. This means, when hitting a list, the function automatically hit on each element of the list. And the result will be the list of all hits. Using Listable functions is the fastest  list traversal. 

One can also Map (/@), Apply (@@, @@@), Thread, MapThread the function (head) to a list if Listable attribute is not available. 

Note that some binary operations have been defined between lists and between list and scalar. This build-in operations between lists are fast and can be in parallel. 

2-2. iterate a list

If one has to explicitly traverse a list,  Mathematica allows iterating a list directly. It is faster than using iterators. The iteration concept is not limited to list. In Do and other iterating cases, list can also be iterated directly. 

 2-3. use build-in list manipulation

Mathematica has defined many build-in list operations. Using these build-in operations are generally faster.