Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Update II  Sample code for simulating boson-sampling experiments has been added (as an answer).

This code exploits new Mathematica capabilities relating to both empirical and smooth distributions; in particular KolmogorovSmirnovTest[__] finds use.

Update I  Multiple optimizations that were suggested by members "ssch" and Simon Woods have in aggregate yielded a ~5X code-speedup; and these optimizations now are incorporated in the example code.

Thank you both very much.

Further improvements are welcome, needless to say. In particular, for nxn matrix arguments, a further O(n) speedup can be achieved (in principle) by exploiting the Gray code structure of "[Delta]PermutationList". However, this would come at the cost of substantially increased code complexity and generally larger round-off error.


For research in BosonSampling (for example) it is desirable to compute matrix permanents by the fastest feasible algorithm. The appended Mathematica code uses Glynn's formula to compute the complex-valued matrix permanent. This code computes the permanent of a 20x20 matrix in ~250 ms (on a 2.93 GHz MacBook Pro laptop)

The Question Asked  Can further speed gains be achieved in numerical computation of the (complex-valued) matrix permanent?

The matrices of interest typically have dimension 10x10 to 25x25, and speed-of-execution for repeated permanent evaluations at fixed matrix-dimension is the sole figure-of-merit.

Suggestions for improvements will cheerfully be adopted!

--- code follows ---


BeginPackage["PermanentCode`"];

Permanent::usage = "\<\
Permanent[mArg_List/;MatrixQ[a]] is computed by Glynn's formula.

The algorithm requires O(m^2 2^m) operations, where m is the dimension 
of the matrix arg.

When the argument is numeric, compiled C-code is executed.

When the argument is non-numeric, a \"Permanent::symbolic\" message 
is issued, and the permanent is calculated symbolically.

Implementation Notes: 

(1) Glynn's formula is simplified with a view to speed-by-simplicity 
    (at negligible cost in formal efficiency); in brief the algorithm
    is implemented as a sequence of BLAS-compatible calls to built-in 
    Mathematica (BLAS) functions.

(2) At present the algorithm does not fully exploit the Gray-code 
    structure of the permutation list \[Delta]PermutationList.    

URL: http://en.wikipedia.org/wiki/Computing_the_permanent#Glynn_formula
URL: http://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms\>";

Permanent::symbolic = "\<\
ADVISORY: Permanent[_] argument is a non-numeric `1`\[Cross]`2` matrix\>";
On[Permanent::symbolic];

\[Delta]PermutationList::usage = "\<\
List of Gray-code permutations, saved in-memory 
for use by Permanent[_]'s Glynn formula.\>";

classicalPermanent::usage = "\<\
classicalPermanent[_] computes the matrix permanent
(slowly!) by expansion of the index permutation\
\>";

Begin["`Private`"];

classicalPermanent[mArg_] := Block[
    {rowList,colPerms},
    rowList = Table[i,{i,1,mArg//Length}];
    colPerms = rowList//Permutations;
    Map[
        (MapThread[mArg[[#1,#2]]&,{rowList,#}]//
          Times@@#&)&,
        colPerms
    ]//Plus@@#&
];

\[Delta]PermutationList[1] = {{1}};
\[Delta]PermutationList[m_Integer]/;(m>1) := (
    (* Conserve memory by purging irrelevant DownValues.
       These rules may exist for arbitrarily large arguments, 
       so a pattern-matched undefine "=." is applied *)
    (\[Delta]PermutationList//DownValues)[[All,1]]//
      ReplaceAll[#,HoldPattern[\[Delta]PermutationList[a_]]:>a]&//ReleaseHold//
        Select[#,(IntegerQ[#]&&(#!=1)&&(#!=m-1))&]&//
          Map[(\[Delta]PermutationList[#]=.;)&,#]&;
    (* now define-and-return the requested \[Delta]PermutationList *)
    \[Delta]PermutationList[m] = \[Delta]PermutationList[m-1]//
    (* idiom: the pipe holds \[Delta]PermutationList[m-1], so conserve 
       memory by deleting its DownValue immediately *)
    (If[m>2,\[Delta]PermutationList[m-1]=.;];#)&//(
          (* reflect DownValue[m-1] in Gray-code order *)
          Map[({1,1}~Join~(#//Rest))&,#] ~ Join ~
          Map[({1,-1}~Join~(#//Rest))&,#//Reverse]
      )&
);

Permanent[ (* numeric evaluation *)
    mArg_List/;(MatrixQ[mArg,NumericQ])
] := compiledGlynnAlgorithm[
        \[Delta]PermutationList[mArg//Length],
        mArg
    ]//Total[#[[1 ;; ;; 2]]] - Total[#[[2 ;; ;; 2]]]&//
      #/2^((mArg//Length)-1)&

compiledGlynnAlgorithm = Compile[{
        {d, _Integer, 1},
        {a, _Complex, 2}
    },
    Apply[Times,(d.a)],
    CompilationTarget -> "C", 
    RuntimeAttributes -> {Listable},
    Parallelization -> True
];

Permanent[ (* symbolic evaluation *)
    mArg_List/;
        (
            MatrixQ[mArg] &&
            (!MatrixQ[mArg,NumericQ]) &&
            (mArg//Length//Message[Permanent::symbolic,#,#]&;True)
        )
] := Map[
    Apply[Times,(#.mArg)]&,
    \[Delta]PermutationList[mArg//Length]
]//Total[#[[1 ;; ;; 2]]] - Total[#[[2 ;; ;; 2]]]&//
  #/2^((mArg//Length)-1)&;

End[];
EndPackage[];

Code to validate and benchmark

nPerm = 4;
Table[\[DoubleStruckCapitalC][i,j],{i,1,nPerm},{j,1,nPerm}]//
  Permanent[#]-classicalPermanent[#]&//
    Expand//
      If[
          #===0,
          Print["VALIDATED: ",nPerm,"\[Cross]",nPerm," symbolic permanent"];,
          Print["ERROR: ",nPerm,"\[Cross]",nPerm," symbolic permanent"];
      ]&;

nPerm = 5;
Table[\[DoubleStruckCapitalC][i,j],{i,1,nPerm},{j,1,nPerm}]//
  Permanent[#]-classicalPermanent[#]&//
    Expand//
      If[
          #===0,
          Print["VALIDATED: ",nPerm,"\[Cross]",nPerm," symbolic permanent"];,
          Print["ERROR: ",nPerm,"\[Cross]",nPerm," symbolic permanent"];
      ]&;

nPerm = 6;
nPerm//{#,#}&//(
            1*RandomVariate[NormalDistribution[0,1],#]+
            I*RandomVariate[NormalDistribution[0,1],#]
  )*1/Sqrt[2]&//
    {Permanent[#],classicalPermanent[#]}&//
      (#[[1]]-#[[2]])/Sqrt[#[[2]]\[Conjugate]*#[[2]]]&//
      If[
          Abs[#]<1000*10^(-$MachinePrecision),
          Print["VALIDATED: ",nPerm,"\[Cross]",nPerm," compiled numeric permanent"];,
          Print["ERROR: ",nPerm,"\[Cross]",nPerm," compiled numeric permanent"];
      ]&;

nPerm = 7;
nPerm//{#,#}&//(
            1*RandomVariate[NormalDistribution[0,1],#]+
            I*RandomVariate[NormalDistribution[0,1],#]
  )*1/Sqrt[2]&//
    {Permanent[#],classicalPermanent[#]}&//
      (#[[1]]-#[[2]])/Sqrt[#[[2]]\[Conjugate]*#[[2]]]&//
      If[
          Abs[#]<100*10^(-$MachinePrecision),
          Print["VALIDATED: ",nPerm,"\[Cross]",nPerm," compiled numeric permanent"];,
          Print["ERROR: ",nPerm,"\[Cross]",nPerm," compiled numeric permanent"];
      ]&;

Print["--------------"];
Print["*** first Permanent[_] evaluation ***"];
Do[
nPerm//{#,#}&//(
            1*RandomVariate[NormalDistribution[0,1],#]+
            I*RandomVariate[NormalDistribution[0,1],#]
  )*1/Sqrt[2]&//(
        (* first call stores Gray-code array *)
        (Permanent[#]//AbsoluteTiming)//First//1000*#&//Round//
          Print["Benchmark: Permanent[ ",nPerm,"\[Cross]",nPerm," ] took ",#," ms"]&;
    )&;,{nPerm,20,12,-1}];

Print["--------------"];
Print["*** second Permanent[_] evaluation ***"];
Do[
nPerm//{#,#}&//(
            1*RandomVariate[NormalDistribution[0,1],#]+
            I*RandomVariate[NormalDistribution[0,1],#]
  )*1/Sqrt[2]&//(
        (* second call runs fast *)
        Permanent[#]; 
        (Permanent[#]//AbsoluteTiming)//First//1000*#&//Round//
          Print["Benchmark: Permanent[ ",nPerm,"\[Cross]",nPerm," ] took ",#," ms"]&;
    )&;,{nPerm,20,12,-1}];

Print["--------------"];
Print["*** (large) Permanent[ 25\[Cross]25 ] evaluation ***"];

25//{#,#}&//(
            1*RandomVariate[NormalDistribution[0,1],#]+
            I*RandomVariate[NormalDistribution[0,1],#]
  )*1/Sqrt[2]&//(
        (Permanent[#]//AbsoluteTiming)//First//Round//
          Print["Benchmark: Permanent[ ",25,"\[Cross]",25," ] took ",#," s"]&;
        (Permanent[#]//AbsoluteTiming)//First//Round//
          Print["Benchmark: Permanent[ ",25,"\[Cross]",25," ] took ",#," s"]&;
    )&;

Results of validating and benchmarking

VALIDATED: 4\[Cross]4 symbolic permanent
VALIDATED: 5\[Cross]5 symbolic permanent
VALIDATED: 6\[Cross]6 compiled numeric permanent
VALIDATED: 7\[Cross]7 compiled numeric permanent
--------------
*** first Permanent[_] evaluation ***
Benchmark: Permanent[ 20\[Cross]20 ] took 1908 ms
Benchmark: Permanent[ 19\[Cross]19 ] took 926 ms
Benchmark: Permanent[ 18\[Cross]18 ] took 428 ms
Benchmark: Permanent[ 17\[Cross]17 ] took 235 ms
Benchmark: Permanent[ 16\[Cross]16 ] took 120 ms
Benchmark: Permanent[ 15\[Cross]15 ] took 52 ms
Benchmark: Permanent[ 14\[Cross]14 ] took 22 ms
Benchmark: Permanent[ 13\[Cross]13 ] took 13 ms
Benchmark: Permanent[ 12\[Cross]12 ] took 8 ms
--------------
*** second Permanent[_] evaluation ***
Benchmark: Permanent[ 20\[Cross]20 ] took 250 ms
Benchmark: Permanent[ 19\[Cross]19 ] took 118 ms
Benchmark: Permanent[ 18\[Cross]18 ] took 55 ms
Benchmark: Permanent[ 17\[Cross]17 ] took 27 ms
Benchmark: Permanent[ 16\[Cross]16 ] took 16 ms
Benchmark: Permanent[ 15\[Cross]15 ] took 22 ms
Benchmark: Permanent[ 14\[Cross]14 ] took 10 ms
Benchmark: Permanent[ 13\[Cross]13 ] took 4 ms
Benchmark: Permanent[ 12\[Cross]12 ] took 1 ms
--------------
*** (large) Permanent[ 25\[Cross]25 ] evaluation ***
Benchmark: Permanent[ 25\[Cross]25 ] took 75 s
Benchmark: Permanent[ 25\[Cross]25 ] took 11 s

Note that the initial evaluation is slower than subsequent evaluations, because initial evaluation creates Gray-code tables that are retained for subsequent use.

share|improve this question

4 Answers 4

up vote 12 down vote accepted

Looking at CompilePrint[compiledGlynnAlgorithm] there are some CopyTensor in it which aren't really needed. There's also a few CoerceTensor in there when it might be faster to just coerce the integer matrix once at the beginning.

By slightly adjusting the function all CopyTensor and CoerceTensor go away giving a small increase in speed:

compiledGlynnAlgorithmAlt = Compile[{
    {d, _Complex, 2}, {a, _Complex, 2}}, 
   Total@Map[Apply[Times, (#.a)*#] &, d],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True];


n = 20;
rc = RandomComplex[{-I - 1, I + 1}, {n, n}];
a = compiledGlynnAlgorithmAlt[δGrayCodeList[n], rc]; // AbsoluteTiming
b = compiledGlynnAlgorithm[δGrayCodeList[n], rc]; // AbsoluteTiming
a == b
(* {0.582192, Null} *)
(* {0.690600, Null} *)
(* True *)

Some more performance can be squeezed out by caching the resulting sign of each row in δGrayCodeList[n] the result is no longer exactly the same, but the relative difference is small:

δGrayCodeListSigns[n_] := δGrayCodeListSigns[n] = Times @@@ δGrayCodeList[n]

compiledGlynnAlgorithmKnownSign = 
  Compile[{{d, _Integer, 2}, {a, _Complex, 2}, {s, _Integer, 1}},
   s.Map[ Apply[Times, (#.a)] &, d]
   , CompilationTarget -> "C"
   , RuntimeAttributes -> {Listable}];

n = 20;
rc = RandomComplex[{-I - 1, I + 1}, {n, n}];

a = compiledGlynnAlgorithmAlt[δGrayCodeList[n], rc]; // AbsoluteTiming
b = compiledGlynnAlgorithm[δGrayCodeList[n], rc]; // AbsoluteTiming
c = compiledGlynnAlgorithmKnownSign[
      δGrayCodeList[n], rc, δGrayCodeListSigns[n]
    ]; // AbsoluteTiming

Abs[c - b]/Abs[b]

(* {0.565806, Null} *)
(* {0.614640, Null} *)
(* {0.430388, Null} *)
(* 2.49266*10^-13 *)
share|improve this answer
    
I affirm your speed-up of approximately 10% &amp; will award this the "answer" if nothing better appears. Thank you "ssh"! –  John Sidles Dec 2 '13 at 19:38
    
Your answer is "accepted". Thank you, "ssch". –  John Sidles Dec 3 '13 at 2:17
    
The code now cache's the sign, with a further speed-up that for Gray code ordering, the sign-list is a strictly alternating sequence of zeros and ones. –  John Sidles Dec 3 '13 at 10:58

You might get a speed up by restricting compiledGlynnAlgorithm to work on just one row of the Gray Code list, allowing the Listable and Parallelization to come into play. I say "might" because the speed up will depend on the details of your hardware.

Redefine compiledGlynnAlgorithm like so (note that it now takes a one dimensional list for d):

compiledGlynnAlgorithm = Compile[{{d, _Integer, 1}, {a, _Complex, 2}},
  Apply[Times, (d.a) d], 
   CompilationTarget -> "C", RuntimeAttributes -> {Listable}, Parallelization -> True]

And put the Total into Permanent

Permanent[mArg_List /; (MatrixQ[mArg, NumericQ])] := 
  Total@compiledGlynnAlgorithm[δGrayCodeList[mArg // Length], mArg] // 
    #/2^((mArg // Length) - 1) &;

a bit more speed

As ssch suggested, a little more performance can be squeezed out by exploiting the fact that the product of a given row of the Gray Code list is either 1 or -1. Furthermore, these occur alternately. So we can redefine compiledGlynnAlgorithm to remove the multiplication by d:

compiledGlynnAlgorithm = Compile[{{d, _Integer, 1}, {a, _Complex, 2}},
  Apply[Times, (d.a)],
  CompilationTarget -> "C", RuntimeAttributes -> {Listable}, Parallelization -> True]

and modify Permanent to Total the odd and even rows of the result separately:

Permanent[mArg_List /; (MatrixQ[mArg, NumericQ])] :=
 Module[{x},
  x = compiledGlynnAlgorithm[δGrayCodeList[mArg // Length], mArg];
  (Total[x[[;; ;; 2]]] - Total[x[[2 ;; ;; 2]]]) // #/2^((mArg // Length) - 1) &]

On my machine this gives about a factor of 3.5 speed increase over the original code for a 20x20 matrix.

share|improve this answer
    
Simon Woods, your optimizations are outstanding, and I have incorporated them (with thanks) into the code. Thank you very much. –  John Sidles Dec 3 '13 at 11:42

Here is a variant adapted from this MathGroup thread

permanentC = 
  Compile[{{m, _Real, 2}}, With[{len = Length[m]}, (-1)^len*Module[
      {s = {0.}, u = 0.},
      Do[
       s = N[IntegerDigits[n, 2, len]];
       u += (-1)^Round[Total[s]]*(Times @@ (m.s)),
       {n, 2^len - 1}];
      u]], CompilationTarget -> "C"];

I checked it on the test set below.

SeedRandom[11111];
testmats = Table[RandomInteger[1, {n, n}], {n, 8, 20, 2}];

It is slightly faster than Permament from the original post. It is also slightly wronger, so to speak. The issue is cancellation error, and for the larger dimensions it shows up in the last few places. The culprit is the Times @@ (m.s) part. Those get large and we lose digits on cancelling. Possibly there is a way to reorder things so as to avoid this numeric pitfall, but offhand I don't see it.

share|improve this answer
    
When I adapted this compiled code to accept _Complex input matrices, the resulting execution times were ~2.4X longer than with the "ssch" approach. Hence, this approach is not preferred (AFAICT). –  John Sidles Dec 2 '13 at 22:33

Having upgraded to Mathematica 10.0.0, it turns out that the Distribution library greatly facilitates the numerical simulation of boson-sampling experiments using the above PermanentCode` package.

The built-in function KolgomorovSmirnovTest is particularly valuable; the Wikipedia entry Kolmogorov–Smirnov test provides a good introduction.

The appended Mathematica code simulates the permanent-distribution of the output photons of the experiment describe by Lund et al "Boson sampling from a Gaussian state" (see PRL 2014 and arXiv:1305.4346).

Background

This example is preparatory to a planned TCS StackExchange question "Can Kolmogorov-Smirnov tests collapse the polynomial hierarchy?"

Example output

Sample code follows to produce the following results and graphics (for nPhoton = 20):

one-sample "A"    K-S test: p = 0.837777
one-sample "B"    K-S test: p = 0.362924
two-sample "AvsB" K-S test: p = 0.831093

Graphical output:

CDF of permanent-squared

Code

As usual, making the graphics look nice takes the most code:

(* ------------------------------------------------ *)
(* --- set boson-sampling simulation parameters --- *)
(* ------------------------------------------------ *)

nPhoton = 10;         (* number of photons detected: 
                         n = 10 runs fast; n = 20 runs slow *)
nSampleMax = 10^5;    (* upper bound to matrix samples; 
                         nSampleMax >= 10^5 is recommended *)
tSampleMax = 24*3600; (* time-used upper bound in seconds *)
kSample = 32;         (* number of Kolmogorov-Smirnov samples *)

(* --------------------------------------- *)
(* --- construct Haar-random unitaries --- *)
(* --------------------------------------- *)

mNode = nPhoton^2;
iSeed=2^nPhoton;

SeedRandom[iSeed];
Umatrix = RandomVariate[NormalDistribution[],{mNode,mNode}] + 
    I * RandomVariate[NormalDistribution[],{mNode,mNode}]//
  SingularValueDecomposition//
    #[[1]].ConjugateTranspose[#[[3]]]&;

(* ------------------------------------------------- *)
(* --- set the scale of the median |permanent|^2 --- *)
(* ------------------------------------------------- *)

PessoanPostulate::usage = "\<\
Per the boson-sampling experiments of Lund et al. \"Boson sampling 
from a gaussian state\" (PRL 2014, see Figure 1), let $n$ be the 
number of photons detected among $m=n^2$ output modes.  Then for 
a Haar-distributed unitary scattering the median value of the 
squared permanent is (empirically) $2n\\Gamma(n+1/2)/m^n$.\
\>";

PessoanPostulate = 2 nPhoton Gamma[nPhoton+1/2]/mNode^nPhoton;

(* ------------------------------------------------------ *)
(* --- sample combinatorically random output channels --- *)
(* ------------------------------------------------------ *)

permEstimatedRMS = Sqrt[PessoanPostulate//N];
amplitudeScale = permEstimatedRMS^(-1.0/nPhoton);

SeedRandom[iSeed+1];
sample$PermanentDistribution = For[
    iSample=0,
    iSample<nSampleMax,
    iSample++,
    Umatrix//RandomChoice[#,nPhoton]&//
      Transpose//RandomChoice[#,nPhoton]&//
        (* from a superabundance of caution, rescale 
           such that the computed permament is \[ScriptCapitalO](1) *)
        (#*amplitudeScale)&//Permanent//
          (#*permEstimatedRMS)&//Sow;
]//Hold//
  TimeConstrained[#//ReleaseHold,tSampleMax]&//
    Reap//Last//Last;

normedSortedData = sample$PermanentDistribution//
  Map[(1/PessoanPostulate)*(#*#\[Conjugate]//Re)&,#]&//Sort;

(* ----------------------------------------------------------- *)
(* --- construct distributions both empirical and smoothed --- *)
(* ----------------------------------------------------------- *)

combinatorialDistribution = normedSortedData//
  Log[10,#]&//
    EmpiricalDistribution;

empiricalDistribution = normedSortedData//
  Rule[#,(#//Log[10,#]&)]&//
    EmpiricalDistribution//Sow[#,"empiricalD"]&;

smoothDistribution = empiricalDistribution//
  RandomVariate[#,{10*(normedSortedData//Length)//Round}]&//
    SmoothKernelDistribution[#,0.1]&//Sow[#,"smoothD"]&;

(* ------------------------------------------------ *)
(* --- construct inverse distributions and CDFs --- *)
(* ------------------------------------------------ *)

inverseCDF = empiricalDistribution//InverseCDF;
forwardCDF = empiricalDistribution//CDF;
inverseSmoothCDF = smoothDistribution//InverseCDF;
inverseCombinatorialCDF = combinatorialDistribution//InverseCDF;

(* ------------------------------------------------------ *)
(* --- simulate k-sample experiments by Alice and Bob --- *)
(* ------------------------------------------------------ *)

SeedRandom[iSeed+2];
AliceBobSimulationData = smoothDistribution//
  RandomVariate[#,{2,kSample}]&;

{{"one-sample \"A\"   ","one-sample \"B\"   ","two-sample \"AvsB\""},{
    KolmogorovSmirnovTest[AliceBobSimulationData[[1]],smoothDistribution],
    KolmogorovSmirnovTest[AliceBobSimulationData[[2]],smoothDistribution],
    KolmogorovSmirnovTest[AliceBobSimulationData[[1]],AliceBobSimulationData[[2]]]
}}//Transpose//
  Map[Print[#[[1]]," K-S test: p = ",#[[2]]]&,#]&;

(* ----------------------------- *)
(* --- plot it all up nicely --- *)
(* ----------------------------- *)

nPlotPts = 500;

range = normedSortedData//Log[10,#]&//
  {#[[4;;5]]//Mean,#[[-5;;-4]]//Mean}&//
    Map[forwardCDF,#]&;

theSmoothSample = Range[1/2,nPlotPts]/nPlotPts//
  Select[#,(#>range[[1]])&]&//
    Select[#,(#<range[[2]])&]&//
      Map[{#,inverseSmoothCDF[#]}&,#//N]&;

smoothPlot = theSmoothSample//
  ListPlot[
    #,
    PlotJoined->True,
    PlotRange->{{0,1},{-3,3}},
    PlotStyle->{
      Directive[Black,AbsoluteThickness[1.8]]
    },
    AspectRatio->0.8,
    AxesOrigin->{0.5,0.0},
    Ticks->{{None,None},{None,None}},
    AxesStyle->Directive[Black,AbsoluteThickness[1.2]],
    Frame->True,
    FrameStyle->Directive[Black,AbsoluteThickness[1.2]],
    FrameTicks -> {
        {{
          Range[-3,3,1],
          {"0.001","0.01","0.1","1","10","100","1000"}
        }//Transpose,None},
        {{
          Range[0,1,0.2],
          {"0","0.2","0.4","0.6","0.8","1"}
        }//Transpose,None}
    },
    FrameTicksStyle->Directive[Black,AbsoluteThickness[0.6],FontSize->Medium],
    GridLines\[Rule]{
        Range[0,1,0.1],
        Outer[#1+#2&,Range[-3,2,1],Range[1,10,1]//Log[10,#]&]//
          Flatten
    },
    GridLinesStyle\[Rule]Directive[Black,AbsoluteThickness[0.6],Opacity[0.2]]
  ]&;

AliceBobInverseCDFPlot = AliceBobSimulationData//
  Map[((#//EmpiricalDistribution//
        InverseCDF[#]&)[x])&,#]&//
  Plot[#,{x,0,1},
      Exclusions -> None, Frame -> None, GridLines -> None, Axes -> None,
      PlotPoints->400,MaxRecursion->3,
      PlotRange->{{0,1},{-3,3}},
      PlotStyle->{
          Directive[RGBColor[0.8,0,0],AbsoluteThickness[1.2],Opacity[0.65]],
          Directive[RGBColor[0,0,0.8],AbsoluteThickness[1.2],Opacity[0.65]]
      }
  ]&;

AliceBobRegionPlot = AliceBobSimulationData//
  Map[((#//EmpiricalDistribution//
        InverseCDF[#]&)[x])&,#]&//
  RegionPlot[
        {
          y<#[[1]] && y>#[[2]],
          y<#[[2]] && y>#[[1]]
        },
        {x,0,1},{y,-3,3},
        Background -> None,
        Frame -> None, GridLines -> None, Axes -> None,
        PlotRange->{{0,1},{-3,3}},
        PlotPoints->200,MaxRecursion->2,
        BoundaryStyle -> None,
        PlotStyle -> {
          Directive[Red,Opacity[0.125]],
          Directive[Blue,Opacity[0.125]]
        }
  ]&;

Show[smoothPlot,AliceBobRegionPlot,AliceBobInverseCDFPlot,
    PlotLabel -> Style[
        "boson-sampling Kolmogorov-Smirnov analysis \n(n="<>
            (nPhoton//ToString)<>
            " photons, k=" <>
            (kSample//ToString)<>
            " detections, m=" <>
            (mNode//ToString)<>
            " modes)",
        FontSize->Medium,
        Black
    ],
    Background -> None,
    FrameLabel->{
      Style["cumulative probability",FontSize->Medium,Black],
      Style["(inverse CDF)/(2n\[CapitalGamma](n+1/2)m^-n)",FontSize->Medium,Black]
    }
]
share|improve this answer
    
Bugfix: in the above code RandomChoice[__] should be RandomSelection[__] ... the shape of the CDF doesn't alter much ... I'll post amended code in a day or two. –  John Sidles 17 hours ago

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.