vegas.mm

(* ::Package:: *)

BeginPackage["Vegas`",{"ErrorBarPlots`","ErrorBarLogPlots`"}];


MergePlots::usage="MergePlots[list] statistically combines a list of plots. MergePlots[list, True] also returns the bin-wise chisq";
IntegrateHistogram::usage="IntegrateHistogram[plot] integrates a histogram.";
CombinePlots::usage="CombinePlots[A,B] adds the histograms A and B, CombinePlots[A,B,op->f] combines them using the function f";
ApplyPlot::usage="ApplyPlot[a,op] applies the operator op to the histogram a";
ScalePlot::usage="Rescales a plot by a factor such that the integral remains unchanged";


MergeNumbers::usage="MergeNumbers[{ {y1, e1}, {y2, e2}, ..}] statistically combines a list of numbers";
PlusNumbers::usage="PlusNumbers[{y1,e1},{y2,e2},...] computes y1+y2+.. and the correct error";
DivideNumbers::usage="DivideNumbers[{y1,e1},{y2,e2}] computes y1/y2 and the correct error";
TimesNumbers::usage="TimesNumbers[{y1,e1},{y2,e2}] computes y1 y2 and the correct error";


CombinePlots::tdlen = "Objects of unequal length";
CombinePlots::wrongx = "X values do not match";
Options[CombinePlots] = {op -> Plus};


AddTarLists::usage="AddTarLists can add two lists from LoadTarFile";
MergeTarLists::usage="MergeTarLists can merge multiple lists from LoadTarFile";
ScaleTarList::usage="ScaleTarList rescales a list from LoadTarFile";


MergeBins::usage="MergeBins[plot, n] merges n adjacent bins into one, reducing the error";
ErrorBandPlot::usage="ErrorBandPlot produces an error band plot from a list (see ErrorListPlot)";

$LoadVerbosity={"\[Chi]n1"};


GuessGood\[Xi]::usage = "GuessGood\[Xi][filename_, flavour_, pieceR_, pieceV_, min_:-4] finds \
the intersection of the real and virtual \!\(\*SubscriptBox[\(\[Xi]\), \(cut\)]\) dependences";

LoadTarFile::usage="LoadTarFile[filename] loads all .vegas files from a tar file and returns a list of rules. ";
LoadTarFile::shamis="SHA checksums of `1` files do not match (`2`).";
LoadTarFile::nenough="`1` is not enough datasets";
(*Prefactor::usage="Prefactor is the prefactor restored during LoadTarFile";
Options[LoadTarFile]={Prefactor\[Rule]1(*conv \[Alpha]em^2*)};*)


BuildPattern::usage="BuildPattern[whichpiece, flavour] creates a regular expression to match the correct files.";
WithPattern::usage="WithPattern[whichpiece, ... , expr] creates an enviroment where \[Sigma] can be used to build the pattern. Avaible options include the folder to search (i.e. Folder->\"nnlo/out\") and the obserable (Observable->\"30\").";
Observable::usage="The observable as specified by the O variable of the vegas file";
Folder::usage="The folder to search. This can be a pattern.";

Begin["`Private`"];

ErrorListPlot[{{DirectedInfinity[-1], u__}, stuff___}] := ErrorListPlot[{stuff}]
ErrorListPlot[{stuff___, {DirectedInfinity[+1], o__}}] := ErrorListPlot[{stuff}]
ErrorListLogPlot[{{DirectedInfinity[-1], u__}, stuff___}] := ErrorListLogPlot[{stuff}]
ErrorListLogPlot[{stuff___, {DirectedInfinity[+1], o__}}] := ErrorListLogPlot[{stuff}]
ErrorListLogLinearPlot[{{DirectedInfinity[-1], u__}, stuff___}] := ErrorListLogLinearPlot[{stuff}]
ErrorListLogLinearPlot[{stuff___, {DirectedInfinity[+1], o__}}] := ErrorListLogLinearPlot[{stuff}]
ErrorListLogLogPlot[{{DirectedInfinity[-1], u__}, stuff___}] := ErrorListLogLogPlot[{stuff}]
ErrorListLogLogPlot[{stuff___, {DirectedInfinity[+1], o__}}] := ErrorListLogLogPlot[{stuff}]


BuildPattern[whichpiece_,flavour_,folder_String]:=RegularExpression[folder<>"/"<>whichpiece<>"_"<>flavour<>"_S\\d*X([\\d\\.]*.)D([\\d\\.]*).*"]->"{$1,$2}"
BuildPattern[whichpiece_,flavour_,opts___Rule]:=Module[{folder, obs},
    folder = Folder/.{opts}; If[folder === Folder, folder=".*"];
    obs = Observable/.{opts}; If[obs === Observable, obs=".*", obs="O"<>obs<>"\\.vegas"];
    RegularExpression[folder<>"/"<>whichpiece<>"_"<>flavour<>"_S\\d*X([\\d\\.]*.)D([\\d\\.]*)_.*"<>obs]->"{$1,$2}"
]
sigmasign = ToExpression@If[MemberQ[Contexts[], "X`"], "X`DiracS", "Global`\[Sigma]"]

Attributes[WithPattern]={HoldAll};
WithPattern[conf__,body_]:=body/.{sigmasign[x_]:>BuildPattern[x,conf]}


AddTarLists[{"value"->v1_,"chi2a"->c1_,"plots"->p1_},{"value"->v2_,"chi2a"->c2_,"plots"->p2_}]/; p1[[;;,1]]===p2[[;;,1]] :={
  "value"->PlusNumbers[v1,v2],"chi2a"->{c1,c2},
  "plots"->Table[
    p->CombinePlots[p/.p1,p/.p2],
    {p,p1[[;;,1]]}
  ]
}
AddTarLists[a:{(_List->_) ..},b:{(_List->_) ..}] /; a[[;;,1]]===b[[;;,1]] := Table[
  k->AddTarLists[k/.a,k/.b],
  {k,a[[;;,1]]}
]


ScaleTarList[{"value"->v1_,"chi2a"->c1_,"plots"->p1_},f_] :={
  "value"->f v1,"chi2a"->c1,
  "plots"->Table[
    p->ApplyPlot[p/.p1,# f&],
    {p,p1[[;;,1]]}
  ]
}
ScaleTarList[a:{(_List->_) ..},f_] := Table[
  k->ScaleTarList[k/.a,f],
  {k,a[[;;,1]]}
]


MergeTarLists[a:{(_List->_) ..}] /; Equal@@("plots"/.a[[;;,2]])[[;;,;;,1]] :=Module[{c,y,e,plotnames,plotdb},
  plotdb=("plots"/.a[[;;,2]]);
  plotnames=plotdb[[1,;;,1]];
  {c,y,e}=MergeNumbers["value"/.a[[;;,2]],True];
  If[MemberQ[$LoadVerbosity,"plot\[Xi]"],
    Print[ErrorListLogLinearPlot[Thread[{
      ToExpression[a[[;;,1]]][[;;,1,1]],
      ("value"/.a[[;;,2]])[[;;,1]],("value"/.a[[;;,2]])[[;;,2]]

    }],PlotRange->All]];
  ];
  If[MemberQ[$LoadVerbosity,"plot\[Delta]"],
    Print[ErrorListLogLinearPlot[Thread[{
      ToExpression[a[[;;,1]]][[;;,1,2]],
      ("value"/.a[[;;,2]])[[;;,1]],("value"/.a[[;;,2]])[[;;,2]]

    }],PlotRange->All]];
  ];
  If[MemberQ[$LoadVerbosity,"\[Chi]n1"],
    Print["Combination of ",a[[;;,1]],"has \!\(\*SuperscriptBox[\(\[Chi]\), \(2\)]\)=",c];
  ];
  {
    "value"->{y,e},
    "chi2a"->{c,"chi2a"/.a[[;;,2]]},
    "plots"->Table[
      p->MergePlots[p/.plotdb],
      {p,plotnames}
    ]
  }
]

If[10<$VersionNumber<11.4,
    addCompletion = FE`Evaluate[FEPrivate`AddSpecialArgCompletion[#]] &,
(* else *)
    addCompletion[x_] := Null;
]

addCompletion["LoadTarFile" -> {2}]

LoadTarFile[filename_String]:=Module[{filenames},
  filenames=Select[
    Import[filename],
    FileExtension[#]=="vegas"&
  ];
  Thread[
    filenames->Import[
      filename,{
      "TAR",
        filenames
    }]
  ]
]

LoadTarFile[matched_List]:=Module[{plotnames,plots,value},
  If[Length[matched]<2,Message[LoadTarFile::nenough,Length[matched]];Abort[]];
  If[Not[Equal@@("SHA"/.matched)],
    Message[LoadTarFile::shamis,Length[matched],StringRiffle[("SHA"/.matched),","]];
  ];
  plotnames=Complement[
    matched[[1,;;,1]],
    {
      "SHA","iteration","chi2a","value","runtime"
    }
  ];
  value=MergeNumbers["value"/.matched,True];
  plots=Table[
     name->MergePlots[name/.matched],
     {name,plotnames}
  ];

  {
    "value"->value[[2;;]],"chi2a"->{value[[1]],"chi2a"/.matched},
    "plots"->plots
  }
]


LoadTarFile[filename_,pat_RegularExpression]:=Module[{tar,matched},
  If[Head[filename]===String,
    tar=LoadTarFile[filename],
  (*else*)
    tar=filename
  ];
  LoadTarFile[Select[tar,StringMatchQ[#[[1]],pat]&][[;;,2]]]

]

LoadTarFile[filename_,pat_RegularExpression->ans_String]:=Module[{tar,matches,out},
  If[Head[filename]===String,
    tar=LoadTarFile[filename],
  (*else*)
    tar=filename
  ];
  out=SortBy[Normal@GroupBy[
    StringCases[#1,pat->ans]->#2&@@@Select[tar,StringMatchQ[#[[1]],pat]&],
    First,
    LoadTarFile[#[[;;,2]]]&
  ],First];
  If[Length[out]==1,
    out[[1,2]],
  (*else*)
    out
  ]
]

LoadTarFile[filename_,pats_CirclePlus|pats_CircleDot]:=Module[{tar,out,next,scale,npat},
  If[Head[filename]===String,
    tar=LoadTarFile[filename],
  (*else*)
    tar=filename
  ];
  out=None;
  Table[
    next=LoadTarFile[tar,pat];
    (*If[MatchQ[pat,f_ (np_CirclePlus|np_CircleDot)],
      scale=pat/.f_ (np_CirclePlus|np_CircleDot)\[RuleDelayed]f;
      npat=pat/.f_ (np_CirclePlus|np_CircleDot)\[RuleDelayed]np;
      next=ScaleTarList[next,scale];
      Print["Scale",Evaluate[next]];
      If[Or[Head[npat]===CirclePlus,Head[npat]===CircleDot],
        next=MergeTarLists[next];
      ]
    ];*)
    If[Head[pat]===CircleDot,
      next=MergeTarLists[next];
    ];
    If[out===None,
      out=next,
    (*else*)
      out=AddTarLists[out,next]
    ],
    {pat,List@@pats}
  ];
  If[Head[pats]===CircleDot,
    MergeTarLists[out],
  (*else*)
    out
  ]
]

LoadTarFile[filename_,f_ (pats_CirclePlus|pats_CircleDot|pats_RegularExpression|pats_Rule)]:=ScaleTarList[LoadTarFile[filename,pats],f]


addCompletion["GuessGood\[Xi]" -> {2}]
GuessGood\[Xi][filename_, flavour_, pieceR_, pieceV_, min_:-4]:=Module[{
  tar, rdata, vdata, fits, ans
},
  If[Head[filename]===String,
    tar=LoadTarFile[filename],
  (*else*)
    tar=filename
  ];
  rdata = Prepend["value" /. #2, ToExpression[#1][[1, 1]]] & @@@ LoadTarFile[tar, BuildPattern[pieceR, flavour]];
  vdata = Prepend["value" /. #2, ToExpression[#1][[1, 1]]] & @@@ LoadTarFile[tar, BuildPattern[pieceV, flavour]];

  fits = Fit[#[[;; , ;; 2]], {1, Log[Global`\[Xi]]}, Global`\[Xi]]&/@{
    rdata,vdata
  };

  ans = Global`\[Xi] /. First@Solve[Equal @@ fits];
  Print@LogLinearPlot[
      fits, {Global`\[Xi], 10^min, 1},
      GridLines -> {{ans}, {Plus @@ fits /. Global`\[Xi] -> ans}}
  ];
  ans
]



MergeNumbers[list_, quiet_: False] := Block[{chi = 0, weight = 0, value},
  weight = Plus @@ (1/(list[[All, 2]]^2));
  value = Plus @@ (list[[All, 1]]/(list[[All, 2]]^2))/weight;
  chi = 1/(Length[list] - 1) (Plus @@ (((list[[All, 1]] - value)/list[[All, 2]])^2));
  If[Not[quiet], Print[chi]];
  Quiet[If[quiet, {chi, value, Sqrt[1/weight]}, {value, Sqrt[1/weight]}]]
] /. Indeterminate -> 0


PlusNumbers[list__] := {Plus @@ {list}[[All, 1]], Sqrt[Plus @@ ({list}[[All, 2]]^2)]}
DivideNumbers[a_, b_] := {a[[1]]/b[[1]], Sqrt[(b[[2]]^2 a[[1]]^2)/b[[1]]^4 + a[[2]]^2/b[[1]]^2]}
TimesNumbers[a_, b_] := {a[[1]] b[[1]], Sqrt[b[[1]]^2 a[[2]]^2 + a[[1]]^2 b[[2]]^2]}


IntegrateHistogram[{{DirectedInfinity[-1], u_,e_}, stuff___}] := u+IntegrateHistogram[{stuff}]
IntegrateHistogram[{stuff___, {DirectedInfinity[+1], o_,e_}}] := o+IntegrateHistogram[{stuff}]
IntegrateHistogram[x_]/;FreeQ[x,DirectedInfinity|Infinity] := Mean[Differences[x[[All, 1]]]] (  Plus @@ x[[All, 2]] )


MergePlots[list_, chi_:False] := Quiet[Module[{out, chilist},
  out = list[[1]];
  (* When the bin is empty, the zero causes trouble *)
  out[[All, 2]] = Plus @@ (list[[All, All, 2]]/(list[[All, All, 3]]^2))/  Plus @@ (1/(list[[All, All, 3]]^2)) /. Indeterminate -> 0;

  (* When the bin is empty, the zero causes trouble *)
  out[[All, 3]] = Sqrt[1/Plus @@ (1/(list[[All, All, 3]]^2))] /. Indeterminate -> 0;

  chilist = Thread[{
    out[[;;,1]],
    1/(Length[list] - 1) * Plus @@@ ((
        (Transpose[list[[;; , ;; , 2]]] - out[[;; , 2]])^2
      )/(
       Transpose[list[[;; , ;; , 3]]]^2
    ))}];
  If[chi, {out, chilist}, out]
]]


CombinePlots[A_, B_, OptionsPattern[]] := Block[{opE, x1, x2},
   opE[e1_, y1_, e2_, y2_] := Sqrt[D[OptionValue[op][x1, x2], x1]^2 e1^2 + D[OptionValue[op][x1, x2], x2]^2 e2^2] /. {x1 -> y1, x2 -> y2};
   Select[
    Table[Table[
       If[A[[i]][[1]] == B[[j]][[1]], {A[[i]][[1]],
         OptionValue[op][A[[i]][[2]], B[[j]][[2]]],
         opE[A[[i]][[3]], A[[i]][[2]], B[[j]][[3]],
          B[[j]][[2]]]}, ## &[]], {j, Length[B]}], {i,
       Length[A]}] /. {x_} -> x, UnsameQ[#, {}] &]
];


ApplyPlot[A_, fY_] := Block[{x, dummy, opE},
   opE[e_, y_] := Abs[D[fY[x], x]] e /. x -> y;
   Table[{A[[i]][[1]], fY[A[[i]][[2]]],
     opE[A[[i]][[3]], A[[i]][[2]]]}, {i, Length[A]}]
];



ScalePlot[a_,s_]:={1/s,s,s}#&/@a
ScalePlot[a_,sx_,sy_]:={1/sx,sy,sy}#&/@a


MergeBins[{{DirectedInfinity[-1], u_,e_}, stuff___}, n_] := Prepend[
    MergeBins[{stuff}, n],
    {DirectedInfinity[-1], u,e}
]
MergeBins[{stuff___, {DirectedInfinity[+1], o_,e_}}, n_] := Append[
    MergeBins[{stuff}, n],
    {DirectedInfinity[+1], o,e}
]
MergeBins[plot_, n_]/;FreeQ[plot,DirectedInfinity|Infinity] := {
    Plus @@ #[[;; , 1]]/n, Plus @@ #[[;; , 2]]/n, Sqrt[
    Plus @@ (#[[;; , 3]]^2)]/n
} & /@ Partition[plot, n]


ClearAll[ErrorBandPlot]
Options[ErrorBandPlot] = Join[Complement[
    Options[ListLinePlot],
    {
        MaxPlotPoints -> Infinity, InterpolationOrder -> None,
        Joined -> True, Filling -> None, FillingStyle -> Automatic,
        LabelingFunction -> Automatic, LabelingSize -> Automatic,
        LabelStyle -> {}, Mesh -> None, MeshFunctions -> {#1 &},
        MeshShading -> None, MeshStyle -> Automatic,
        PlotLayout -> "Overlaid", PlotLabels -> None, PlotMarkers -> None
    }
],{
  Underflow->False, Overflow->False
}];

(* one plot with x *)
ErrorBandPlot[dat_, opts : OptionsPattern[]] /; Length[dat[[1]]] == 3 := Block[{newdat, ps, diff, s},
    ps = OptionValue[PlotStyle];
    If[ps === Automatic, ps = ColorData[97][1]];

    newdat = Select[dat,FreeQ[DirectedInfinity]];
    diff = First@Differences[newdat[[;; , 1]]];
    newdat = {
        {#1, #2}, ErrorBar[diff/2, #3]
    } & @@@ newdat;

    If[And[OptionValue[Underflow]=!=False, dat[[1,1]]==DirectedInfinity[-1]],
      s = 2;
      If[NumberQ[OptionValue[Underflow]], s = OptionValue[Underflow]];
      PrependTo[newdat, {
        {dat[[+2,1]]-(s+1/2)*diff, dat[[+1,2]]}, ErrorBar[s*diff, dat[[+1,3]]]
      }]
    ];
    If[And[OptionValue[Overflow ]=!=False, dat[[-1,1]]==DirectedInfinity[+1]],
      s = 2;
      If[NumberQ[OptionValue[Overflow]], s = OptionValue[Overflow]];
      AppendTo[newdat, {
        {dat[[-2,1]]+(s+1/2)*diff, dat[[-1,2]]}, ErrorBar[s*diff, dat[[-1,3]]]
      }]
    ];

    ErrorListPlot[
        newdat,
        ErrorBarFunction -> Function[{coords, errs}, {
            Opacity[1], EdgeForm[ps], ps,
            Rectangle[
                coords + {errs[[1, 1]], errs[[2, 1]]},
                coords + {errs[[1, 2]], errs[[2, 2]]}
            ]
        }],
        Epilog->{
            ps,
            Line@Flatten[{
                {#1[[1,1]]-#1[[2,1]],#1[[1,2]]},
                {#1[[1,1]]+#1[[2,1]],#1[[1,2]]}
            }&/@newdat,1]
        },
        PlotStyle -> None, FilterRules[{opts},Options[ErrorListPlot]]
    ]
]

(* one plot, no x *)
ErrorBandPlot[dat_, opts : OptionsPattern[]] /; Length[dat[[1]]] == 2 :=
    ErrorBandPlot[MapThread[Prepend, {dat, Range[Length[dat]]}], opts]


(* multiple plots *)
ErrorBandPlot[dat_, opts : OptionsPattern[]] /; Length[dat[[1]]] > 3 := Block[{newdat, ps,plots},
    ps = OptionValue[PlotStyle];
    If[ps === Automatic, ps = ColorData[97] /@ Range[Length[dat]]];
    If[Head[ps] =!= List, ps = ConstantArray[ps, Length[dat]]];
    plots = MapThread[
        ErrorBandPlot[#1, PlotStyle -> #2, opts] &,
        {dat, ps}
    ];
    Show[plots, Epilog->Join[Epilog/.AbsoluteOptions[plots,Epilog]],opts]
]

End[];
EndPackage[];


Begin["System`Convert`VEGAS`"]


ImportExport`RegisterImport["vegas",
  importer,
  "FunctionChannels"->{"Streams"},"BinaryFormat"->True
]
ImportExport`RegisterImport["vegasL",
  importer[#,8]&,
  "FunctionChannels"->{"Streams"},"BinaryFormat"->True
]

importer[filename_String,opts___]:=Module[{stream,ans},
  stream=OpenRead[filename,BinaryFormat->True];
  ans=importer[stream,opts];
  Close[stream];
  ans
]


ReadFortran[stream_,type_]:=Module[{len,ans,len2},
  len=BinaryRead[stream,"UnsignedInteger32"];
  If[len===EndOfFile, Return[$Failed]];
  ans=BinaryReadList[stream,type,len/Switch[type,
    "Byte"|"Character8"|"Integer8",1,
    "Integer16",2,
    "Integer32"|"Real32",4,
    "Integer64"|"Real64",8
    ]
  ];
  len2=BinaryRead[stream,"UnsignedInteger32"];
  If[len!=len2,Print["Something's fishy"]];
  If[Length[ans]==1,First[ans],ans]
]
WriteFortran[stream_, dat_, type_] := Module[{len},
  len = Switch[type,
    "Byte" | "Character8" | "Integer8", 1,
    "Integer16", 2,
    "Integer32" | "Real32", 4,
    "Integer64" | "Real64", 8
  ];
  If[Head[dat] == List,
    len = Length[dat]*len
  ];
  BinaryWrite[stream, len, "UnsignedInteger32"];
  BinaryWrite[stream, dat, type];
  BinaryWrite[stream, len, "UnsignedInteger32"];
]


GuessVegasVersion[stream_, intsize_:4]:=Module[{nrq, nrbins, inttype, time, magic,version},
  SetStreamPosition[stream,0];
  magic = StringJoin @@ Select[BinaryReadList[stream, "Character8", 17], # != "\000" &];
  If[magic === "\t McMule  \t",
    version=StringJoin @@ ReadFortran[stream, "Character8"];
    version = First@StringCases[
      version,
      RegularExpression["v(\\d+)([NL])"] :> {
        ToExpression["$1"],
        "Integer" <> ToString[8 If["$2" == "L", 8, 4]]
      }
    ];
    Return[version, Module];
  ];

  inttype = "Integer"<>ToString[8 intsize];
  (* Test v1 vs. v2 *)
  SetStreamPosition[stream,6893 + 3 intsize];
  {nrq, nrbins} = ReadFortran[stream, inttype];
  SetStreamPosition[stream,6925 + 5 intsize + 22 nrq + 16 nrbins nrq];
  time=ReadFortran[stream, "Real64"];
  SetStreamPosition[stream,0];
  If[FailureQ[time],
    Return[{1,inttype}, Module],
  (*else*)
    Return[{2,inttype}, Module]
  ];
]


importeverything=False;

importer[stream_InputStream,intsize_:4, opts___]:=Module[{
  len,len2,ans,
  sha,it,ndo,si,swgt,schi,xi,randy,nrq,nrbins,quant,minv,maxv,names, X, delta, msg, version,
  inttype, namelen, nrbinsuo,AddOverUnderflowBins
},
  {version, inttype} = GuessVegasVersion[stream, intsize];
  sha=StringJoin@@ReadFortran[stream,"Character8"];
  it=ReadFortran[stream,inttype];ndo=ReadFortran[stream,inttype];
  si=ReadFortran[stream,"Real64"];swgt=ReadFortran[stream,"Real64"];schi=ReadFortran[stream,"Real64"];

  xi=ReadFortran[stream,"Real64"];

  randy=ReadFortran[stream,inttype];
  If[version>2,
    {nrq,nrbins,namelen}=ReadFortran[stream,inttype],
  (*else*)
    {nrq,nrbins}=ReadFortran[stream,inttype];
    namelen = 6;
  ];
  {minv,maxv}=Partition[ReadFortran[stream,"Real64"],nrq];
  names=StringTrim[StringJoin@@@Partition[ReadFortran[stream,"Character8"],namelen]];
  quant=ReadFortran[stream,"Real64"];

  If[version>1,
    time=ReadFortran[stream,"Real64"];
    msg=StringJoin@@ReadFortran[stream,"Character8"];
    If[Or[StringContainsQ[msg,"Warning"],MemberQ[$LoadVerbosity,"msg"]],
      Print[StringTrim[msg]]
    ];
  ];


  If[version>2,
    AddOverUnderflowBins[u_,o_,l_List] := Join[{u},l,{o}];
    AddOverUnderflowBins[u_,o_,l_?NumberQ] := Join[{u},ConstantArray[l,nrbins],{o}],
  (*else*)
    AddOverUnderflowBins[u_,o_,l_] := l
  ];

  Flatten[{
    "SHA"->sha,"iteration"->it,
    "value"->{si/swgt,Sqrt[1/swgt]},
    "chi2a"->If[it>1,(schi-si^2/swgt)/(it-1),-1],
    If[version>1,"runtime"->time,{}],
    If[it>1,
      delta = (maxv-minv)/nrbins;
      nrbinsuo = If[version>2,nrbins + 2,nrbins];
      Table[
        names[[i]]->Thread[{
          (* X *)
          AddOverUnderflowBins[
            DirectedInfinity[-1],DirectedInfinity[+1],
            minv[[i]]+delta[[i]]*(Range[nrbins]-1/2)
          ],
          (* Y *)
          quant[[i;;nrq nrbinsuo;;nrq]]/it / AddOverUnderflowBins[1,1,delta[[i]]],
          (* E *)
          Sqrt[
            (
              quant[[nrq nrbinsuo+i;;;;nrq]]-quant[[i;;nrq nrbinsuo;;nrq]]^2/it
            )/(it (it-1))] / AddOverUnderflowBins[1,1,delta[[i]]]
        }
      ],{i,nrq}],
    (* else *)
      {}
    ],

    If[importeverything,
        {
            "si"->si,"swgt"->swgt,"schi"->schi,
            "xi"->xi,
            "randy"->randy
        },
    (* else *)
        {}
    ]
  }]
]

tzero = FromCharacterCode[0] <> FromCharacterCode[0] <> FromCharacterCode[0];



ImportExport`RegisterExport["vegas",
  exporter,
  "FunctionChannels"->{"Streams"},"BinaryFormat"->True
]
ImportExport`RegisterExport["vegasL",
  exporter[#,8]&,
  "FunctionChannels"->{"Streams"},"BinaryFormat"->True
]

exporter[filename_String,opts___]:=Module[{stream},
  stream=OpenWrite[filename,BinaryFormat->True];
  exporter[stream,opts];
  Close[stream];
]


exporter[stream_OutputStream, {content__Rule}, intsize_Integer: 4] := Module[{
    keys, inttype,
    sha, iteration, chi2a, y, e, runtime, ndo, xi, randy, msg,
    hkeys, histogram, delta, minv, maxv, nrbins, namelen, nrq, quant
},
  keys = {
    "SHA", "iteration", "chi2a", "value", "runtime", "ndo", "xi",
    "randy", "msg"
  };


  (* histograms *)
  hkeys = Complement[{content}[[;; , 1]], keys];

  (* Lengths of things *)
  nrq = Length[hkeys];
  namelen = Max[Max[StringLength /@ hkeys], 6];

  (* Calculate number of bins from first histogram *)
  histogram = hkeys[[1]] /. {content};
  If[Head[histogram[[1, 1]]] =!= DirectedInfinity,
    nrbins = Length[histogram],
  (*else*)
    nrbins = Length[histogram]-2;
  ];

  (* Prepare variables to be written later *)
  quant = ConstantArray[0, 2 (nrbins + 2) nrq];
  minv = {}; maxv = {};

  Table[
    histogram = hkeys[[i]] /. {content};
    If[Head[histogram[[1, 1]]] =!= DirectedInfinity,
      (* Add OU bins *)
      histogram = Join[
        {{DirectedInfinity[-1], 0, 0}}, histogram, {{DirectedInfinity[-1], 0, 0}}
      ];
    ];

    delta = histogram[[3, 1]] - histogram[[2, 1]];
    AppendTo[maxv, Chop[histogram[[2, 1]] - delta 1/2 + nrbins delta]];
    AppendTo[minv, Chop[histogram[[2, 1]] - delta 1/2]];

    (* delta to array with OU bins *)
    delta = Join[{1}, ConstantArray[delta, nrbins], {1}];

    quant[[i ;; nrq nrbins + 2 ;; nrq]] = delta iteration histogram[[;; , 2]];

    If[And @@ EqualTo[0] /@ histogram[[;; , 3]],
      quant[[nrq (nrbins + 2) + i ;; ;; nrq]] = quant[[i ;; nrq nrbins + 2 ;; nrq]]^2/iteration
    (*else*),
      quant[[nrq (nrbins + 2) + i ;; ;; nrq]] = delta^2 iteration (
        histogram[[;; , 3]]^2 (-1 + iteration) + histogram[[;; , 2]]^2
      )
    ],
    {i, Length@hkeys}
  ];


  (* other data *)
  {sha, iteration, chi2a, {y, e}, runtime, ndo, xi, randy, msg} = keys
    /. {content}
    /. {
      "SHA" ->
      StringTake[Hash[DownValues[In], "SHA", "HexString"], 5],
      "iteration" -> 2,
      "chi2a" -> 0,
      "runtime" :> TimeUsed[],
      "ndo" -> -1(* TODO *),
      "xi" -> ConstantArray[-1, 17 50](* TODO *),
      "randy" -> -1,
      "msg" -> "Warning: Generated with Mathematica"
  };
  e = e + $MachineEpsilon;
  inttype = "Integer" <> ToString[8 intsize];

  (* Magic *)
  BinaryWrite[stream, "\t" <> tzero <> " McMule  \t" <> tzero];
  (* Header *)
  BinaryWrite[stream, "\n" <> tzero <> "v3" <> If[intsize == 8, "L", "N"] <> "       \n" <> tzero];
  WriteFortran[stream, Characters@sha, "Character8"];

  (* Integrator state *)
  WriteFortran[stream, iteration, inttype];
  WriteFortran[stream, ndo, inttype];
  WriteFortran[stream, y/e^2, "Real64"];
  WriteFortran[stream, 1/e^2, "Real64"];
  WriteFortran[stream, chi2a (-1 + iteration) + y^2/e^2, "Real64"];
  WriteFortran[stream, xi, "Real64"];
  WriteFortran[stream, randy, inttype];

  (* Histograms *)

  (* Metadata *)
  WriteFortran[stream, {nrq, nrbins, namelen}, inttype];
  WriteFortran[stream, Join[minv, maxv], "Real64"];
  WriteFortran[stream, Flatten[Characters@StringPadRight[#, namelen] & /@ hkeys], "Character8"];

  (* Data *)
  WriteFortran[stream, quant, "Real64"];

  (* Misc *)
  WriteFortran[stream, runtime, "Real64"];
  WriteFortran[stream, Characters@msg, "Character8"];
]



ImportExport`RegisterImport["badpoints",
  bpimporter,
  "FunctionChannels"->{"Streams"},"BinaryFormat"->True
]

bpimporter[filename_String,opts___]:=Module[{stream,ans},
  stream=OpenRead[filename,BinaryFormat->True];
  ans=bpimporter[stream,opts];
  Close[stream];
  ans
]


bpimporter[stream_InputStream, opts___]:=Module[{len,msg,whichpiece,flavour,xicut,delcut,ndim,x,len2},
  Last@Reap@While[True,
    len=BinaryRead[stream,"UnsignedInteger32"];
    If[len===EndOfFile, Break[]];
    If[len!=30+10+8+1 8+1 8+1 4+17 8,Print["Failed!"]];
    msg=StringTrim[StringJoin@@BinaryReadList[stream,"Character8",30]];
    whichpiece=StringTrim[StringJoin@@BinaryReadList[stream,"Character8",10]];
    flavour=StringTrim[StringJoin@@BinaryReadList[stream,"Character8",8]];
    xicut=BinaryRead[stream,"Real64"];
    delcut=BinaryRead[stream,"Real64"];
    ndim=BinaryRead[stream,"Integer32"];
    x=BinaryReadList[stream,"Real64",17][[;;ndim]];
    len2=BinaryRead[stream,"UnsignedInteger32"];If[len!=len2,Print["Something's fishy"]];
    Sow[{
      "msg"->msg,"which_piece"->whichpiece,"flavour"->flavour,"xicut"->xicut,"delcut"->delcut,"x"->x
    }];
  ]
]


End[];


Unprotect[Import];
Import[name_String,opts___?OptionQ]/;(
  ToLowerCase[FileExtension[name]]
)==="vegas":=Import[name,"vegas",opts];
Import[name_String,opts___?OptionQ]/;(
  ToLowerCase[FileExtension[name]]
)==="bp":=Import[name,"badpoints",opts];
Protect[Import];


Unprotect[Export];
Export[name_String,data_List,opts___?OptionQ]/;(
  ToLowerCase[FileExtension[name]]
)==="vegas":=Export[name,data,"vegas",opts];
Protect[Export];

(* ::Input:: *)
(**)