Empty .png figures generated by a function "PlotWithGrid" whereas no problem with Export function

Question

I have an issue about a .wls script that I can't manage to solve for the moment. I execute it with wolframscript interpeter on command line :

A function PlotWithGrid has been coded this way :

PlotWithGrid[func_,initValueList_?ListQ, intervalList_?ListQ, axesLabels_ , PlotType_, minVec_, exportDir_?StringQ] := Module[
{ tb, f,cons, m, n, pl},
Switch[ PlotType,
1,
f[m_,n_, i_, j_]:=Apply[func,initValueList/.{initValueList[[i]]-> m, initValueList[[j]]-> n}];
tb = ParallelTable[
Piecewise[{{ContourPlot[f[x, y, i, j],{x,intervalList[[j,1]],intervalList[[j,2]]},{y,intervalList[[i,1]],intervalList[[i,2]]},AspectRatio->1,Frame->True,Axes->False, FrameLabel->{axesLabels[[j]], axesLabels[[i]]}],j<i}},Null],
{i, 2, Length@initValueList },
{j, 1, Length@initValueList - 1}
],
2,
tb = ParallelTable[
Piecewise[{{
ListPlot[func[[All, {j, i}]],Frame->True,FrameLabel->{axesLabels[[j]], axesLabels[[i]]},PlotStyle->{Red,PointSize[0.003]},PlotRange->{{intervalList[[j,1]],intervalList[[j,2]]}, {intervalList[[i,1]],intervalList[[i,2]]}},
Axes->False,ImageSize->Medium,
FrameTicks->{{Automatic, None}, {fticks, None}},
Epilog->{{PointSize[Medium],Blue,Point[func[[1,{j, i}]]],{PointSize[Medium],Black,Point[minVec[[{j, i}]]]}}},
AspectRatio->1
]
,j<i}},Null],
{i, 2, Length@initValueList },
{j, 1, Length@initValueList - 1}
],
3,
tb = Table[
Piecewise[{{
Show[
Confid2D[func[[3,All,{j, i}]],Blue],
Confid2D[func[[2,All,{j, i}]],Green],
Confid2D[func[[1,All,{j, i}]],Red],
PlotRange->{{intervalList[[j,1]],intervalList[[j,2]]}, {intervalList[[i,1]],intervalList[[i,2]]}},
FrameTicks->{{Automatic, None}, {fticks, None}},
Epilog->{PointSize[Medium],Black,Point[minVec[[{j, i}]]]},
FrameLabel->{axesLabels[[j]], axesLabels[[i]]},
AspectRatio->1,
ImageSize-> Medium
]
,j<i}},Null],
{i, 2, Length@initValueList },
{j, 1, Length@initValueList - 1}
],
4,
tb = Table[
Piecewise[{{
Show[
ListPlot[func[[3,All,{j, i}]], PlotStyle->{Blue,PointSize[Medium]}, Axes->False, Frame->True],
ListPlot[func[[2,All,{j, i}]], PlotStyle->{Green,PointSize[Medium]}, Axes->False, Frame->True],
ListPlot[func[[1,All,{j, i}]], PlotStyle->{Red,PointSize[Medium]}, Axes->False, Frame->True],
PlotRange->{{intervalList[[j,1]],intervalList[[j,2]]}, {intervalList[[i,1]],intervalList[[i,2]]}},
FrameTicks->{{Automatic, None}, {fticks, None}},
Epilog->{PointSize[Medium],Black,Point[minVec[[{j, i}]]]},
FrameLabel->{axesLabels[[j]], axesLabels[[i]]},
AspectRatio->1,
ImageSize-> Medium
]
,j<i}},Null],
{i, 2, Length@initValueList },
{j, 1, Length@initValueList - 1}
],
5,
cons[i_, j_]:=Quantile[PDF[SmoothKernelDistribution[func[[All,{ i, j}]]],func[[All, {i, j}]]],1-Erf[{1,2}/Sqrt[2]]//N];
f[ m_, n_, i_, j_]:=PDF[SmoothKernelDistribution[func[[All, {i, j}]]],{m, n}];

tb = Table[
Piecewise[{{
f[ m_, n_] = f[ m, n, j, i];
ContourPlot[f[m, n],{m, intervalList[[j, 1]],intervalList[[j,2]]},{n, intervalList[[i, 1]],intervalList[[i,2]]},Contours->cons[i, j],ContourShading->{White,Hue[0.5,0.9,0.9],Hue[0.55,0.8,0.8]},FrameLabel->{axesLabels[[j]], axesLabels[[i]]},PlotRange->{{intervalList[[j,1]],intervalList[[j,2]]}, {intervalList[[i,1]],intervalList[[i,2]]}},ClippingStyle->None, 
FrameTicks->{{Automatic, None}, {fticks, None}},Epilog->{PointSize[Medium],Black,Point[minVec[[{j, i}]]]}]
,j<i}},Null],
{i, 2, Length@initValueList },
{j, 1, Length@initValueList - 1}
]
];

pl=ResourceFunction["PlotGrid"][
tb,
PlotRange-> Max,
ImageSize-> Large,
Method->{"AllCustomTicks"->Automatic, "FixFrameTicks"->True}
];
Print@pl;
If[exportDir != "",
Export[exportDir,pl,"AllowRasterization"->False];
];

Select[tb//Flatten, #=!=Null&]

];

1. If I use the function PlotWithGrid under this form :

   pl = PlotWithGrid[{data1s[[All, variableIdx]], data2s[[All, variableIdx]], data3s[[All, variableIdx]]}, {[CapitalOmega]m,[CapitalOmega]k,H0,[Phi]0,d[Phi]0,[Omega]BD}, varIntervals, axesLabels, 3,chainShort[[minimumChiPos]], "contour_1.png"];

The figure "contour_1.png" is empty and I don't understand why.

2. I have into my code another type of plot like this :

   pl=ListPlot[chain[[All,8]],Joined->True,PlotRange->All,Frame->True,FrameLabel->{"N","!(*SuperscriptBox[([Chi]), (2)])"}]; Print@pl; pics = Join[pics, {pl}]; Export["chi2.png",pl,"AllowRasterization"->False];

In this case, the figure "chi2.png" is well generated (not empty).

3. So I tried to use this Export function to my initial issue by doing :

   pl = PlotWithGrid[{data1s[[All, variableIdx]], data2s[[All, variableIdx]], data3s[[All, variableIdx]]}, {[CapitalOmega]m,[CapitalOmega]k,H0,[Phi]0,d[Phi]0,[Omega]BD}, varIntervals, axesLabels, 3,chainShort[[minimumChiPos]]]; pics = Join[pics, pl]; Export["contour_1.png",pl,"AllowRasterization"->False];

But with this snippet code, none figure "contour_1.png" is produced, nothing during the run.

If someone could indicate me a workaround to create correctly a none empty figure "contour_1.png" figure, this would be fine.

EDIT 1: Here a short example of .wls that don't work in direct command line : $ wolframswript short.wls

#!/usr/bin/env wolframscript
(* ::Package:: *)

(* SetOptions[EvaluationNotebook[],{AutoGeneratedPackage->Automatic,InitializationCellEvaluation->False,InitializationCellWarning->False,StyleDefinitions->Notebook[{Cell[StyleData[StyleDefinitions->"Default.nb"]],Cell[StyleData["Input"],InitializationCell->True]},Visible->False,StyleDefinitions->"PrivateStylesheetFormatting.nb"]}] *)


(* SetDirectory[NotebookDirectory[]]; *)
SetDirectory["."];
Needs["ComputationalGeometry`"];

fticks[min_,max_]:={#,Rotate[#,Pi/2],{ 0.01, 0}}&/@(N@FindDivisions[{min,max},6]);
PlotWithGrid[func_,initValueList_?ListQ, intervalList_?ListQ, axesLabels_ , PlotType_, minVec_, exportDir_?StringQ] := Module[
{ tb, f,cons, m, n, pl},
Switch[ PlotType,
2,
tb = ParallelTable[
Piecewise[{{
ListPlot[func[[All, {j, i}]],Frame->True,FrameLabel->{axesLabels[[j]], axesLabels[[i]]},PlotStyle->{Red,PointSize[0.003]},PlotRange->{{intervalList[[j,1]],intervalList[[j,2]]}, {intervalList[[i,1]],intervalList[[i,2]]}},
Axes->False,ImageSize->Medium,
FrameTicks->{{Automatic, None}, {fticks, None}},
Epilog->{{PointSize[Medium],Blue,Point[func[[1,{j, i}]]],{PointSize[Medium],Black,Point[minVec[[{j, i}]]]}}},
AspectRatio->1
]
,j<i}},Null],
{i, 2, Length@initValueList },
{j, 1, Length@initValueList - 1}
],
3,
tb = Table[
Piecewise[{{
Show[
Confid2D[func[[3,All,{j, i}]],Blue],
Confid2D[func[[2,All,{j, i}]],Green],
Confid2D[func[[1,All,{j, i}]],Red],
PlotRange->{{intervalList[[j,1]],intervalList[[j,2]]}, {intervalList[[i,1]],intervalList[[i,2]]}},
FrameTicks->{{Automatic, None}, {fticks, None}},
Epilog->{PointSize[Medium],Black,Point[minVec[[{j, i}]]]},
FrameLabel->{axesLabels[[j]], axesLabels[[i]]},
AspectRatio->1,
ImageSize-> Medium
]
,j<i}},Null],
{i, 2, Length@initValueList },
{j, 1, Length@initValueList - 1}
]
];

pl=ResourceFunction["PlotGrid"][
tb,
PlotRange-> Max,
ImageSize-> Large,
Method->{"AllCustomTicks"->Automatic, "FixFrameTicks"->True}
];
If[exportDir != "",
Export[exportDir,pl,"AllowRasterization"->False];
];

Select[tb//Flatten, #=!=Null&]

];

Confid2D[data_?ListQ,color_]:=Module[{cv,cv1,cv2},
cv=ConvexHull[data,AllPoints->True];
cv1=data[[#]]&/@cv;
If[Length@cv1 > 1, 
cv2=Join[cv1,{cv1[[1]],cv1[[2]]}],
cv2=Join[cv1,{cv1[[1]]}]];
Show[Graphics[{color,FilledCurve[BSplineCurve[cv2,SplineDegree->2,SplineKnots->"Unclamped"]],Antialiasing->True}],Frame->True]
];

dH[\[Rho]m_,\[Phi]_,u_,\[Omega]BD_, \[CapitalOmega]k_,du_, z_]:= -(((1+z) (16 N[\[Pi]] \[Rho]m-6 (1+z)^2 \[CapitalOmega]k \[Phi]) ((1+z) \[Omega]BD u^3-2 \[Omega]BD u ((1+z) du[\[Rho]m,\[Phi],u,\[Omega]BD, \[CapitalOmega]k,z]+u) \[Phi]-6 du[\[Rho]m,\[Phi],u,\[Omega]BD, \[CapitalOmega]k,z] \[Phi]^2)+(1/(1+z))6 \[Phi] (-8 N[\[Pi]] \[Rho]m+(1+z)^2 \[CapitalOmega]k ((1+z) u+2 \[Phi])) (6 \[Phi]^2-(1+z) u ((1+z) \[Omega]BD u+6 \[Phi])))/(2 Sqrt[(-16 N[\[Pi]] \[Rho]m+6 (1+z)^2 \[CapitalOmega]k \[Phi])/(6 (1+z) u+((1+z)^2 \[Omega]BD u^2)/\[Phi]-6 \[Phi])] ((1+z)^2 \[Omega]BD u^2+6 (1+z) u \[Phi]-6 \[Phi]^2)^2));
d\[Rho]m[\[Rho]m_,\[Phi]_,u_,z_]:=3/(1+z) \[Rho]m;
d\[Phi][\[Rho]m_,\[Phi]_,u_,z_]:= u;
du[\[Rho]m_,\[Phi]_,u_,\[Omega]BD_, \[CapitalOmega]k_, z_]:=-((24 N[\[Pi]] \[Rho]m \[Phi]^3+(1+z) u \[Phi]^2 (8 \[Pi] (-3+\[Omega]BD) \[Rho]m-3 (1+z)^2 (3+2 \[Omega]BD) \[CapitalOmega]k \[Phi])+3 (1+z)^2 u^2 \[Phi] (-4 N[\[Pi]] \[Omega]BD \[Rho]m+(1+z)^2 (3+2 \[Omega]BD) \[CapitalOmega]k \[Phi])+(1+z)^3 \[Omega]BD u^3 (-4 \[Pi] (1+\[Omega]BD) \[Rho]m+(1+z)^2 (3+2 \[Omega]BD) \[CapitalOmega]k \[Phi]))/((1+z)^2 (3+2 \[Omega]BD) \[Phi]^2 (-8 N[\[Pi]] \[Rho]m+3 (1+z)^2 \[CapitalOmega]k \[Phi])));

RK4Method[dH_,d\[Phi]_,d\[Rho]m_,du_,\[CapitalOmega]m_,\[CapitalOmega]k_,H0_,\[Phi]0_,d\[Phi]0_,\[Omega]BD_,zLine0_]:=Module[
{ h, Htable, \[Rho]mtable, \[Phi]table, utable, Hk1, Hk2, Hk3,Hk4, \[Rho]mk1, \[Rho]mk2, \[Rho]mk3,\[Rho]mk4, \[Phi]k1, \[Phi]k2,\[Phi]k3,\[Phi]k4, uk1,uk2, uk3,uk4, containsIndeterminate, containsComplex, Hval,zLine =Join[{0},zLine0]},
containsIndeterminate = False;
containsComplex = False;

Htable=Table[Infinity,{i,1,Length@zLine}];
Htable[[1]]=H0;
\[Rho]mtable=3 H0 H0 \[Phi]0 \[CapitalOmega]m/(8N[\[Pi]]);
utable=d\[Phi]0;
\[Phi]table=\[Phi]0;
Do[
h=(zLine[[i]]-zLine[[i - 1]]);
Hk1=dH[\[Rho]mtable,\[Phi]table,utable,\[Omega]BD, \[CapitalOmega]k,du,zLine[[i - 1]]];
\[Phi]k1=d\[Phi][\[Rho]mtable,\[Phi]table,utable,zLine[[i - 1]]];
\[Rho]mk1=d\[Rho]m[\[Rho]mtable,\[Phi]table,utable,zLine[[i - 1]]];
uk1=du[\[Rho]mtable,\[Phi]table,utable,\[Omega]BD, \[CapitalOmega]k,zLine[[i - 1]]];

Hk2=dH[\[Rho]mtable+h/2*\[Rho]mk1, \[Phi]table+h/2* \[Phi]k1,utable + h/2* uk1, \[Omega]BD, \[CapitalOmega]k,du, zLine[[i - 1]]+h/2];
\[Phi]k2=d\[Phi][\[Rho]mtable+h/2*\[Rho]mk1, \[Phi]table+h/2* \[Phi]k1,utable + h/2* uk1,zLine[[i - 1]]+h/2];
\[Rho]mk2=d\[Rho]m[\[Rho]mtable+h/2*\[Rho]mk1, \[Phi]table+h/2* \[Phi]k1,utable + h/2* uk1,zLine[[i - 1]]+h/2];
uk2=du[\[Rho]mtable+h/2*\[Rho]mk1, \[Phi]table+h/2* \[Phi]k1,utable+ h/2* uk1,\[Omega]BD, \[CapitalOmega]k,zLine[[i - 1]]+h/2];

Hk3=dH[\[Rho]mtable+h/2*\[Rho]mk2,\[Phi]table+h/2* \[Phi]k2,utable+h/2* uk2,\[Omega]BD, \[CapitalOmega]k,du,zLine[[i - 1]]+h/2];
\[Phi]k3=d\[Phi][\[Rho]mtable+h/2*\[Rho]mk2,\[Phi]table+h/2* \[Phi]k2,utable+h/2* uk2,zLine[[i - 1]]+h/2];
\[Rho]mk3=d\[Rho]m[\[Rho]mtable+h/2*\[Rho]mk2,\[Phi]table+h/2* \[Phi]k2,utable+h/2* uk2,zLine[[i - 1]]+h/2];
uk3=du[\[Rho]mtable+h/2*\[Rho]mk2,\[Phi]table+h/2* \[Phi]k2,utable+h/2* uk2,\[Omega]BD, \[CapitalOmega]k,zLine[[i - 1]]+h/2];

Hk4=dH[\[Rho]mtable+h*\[Rho]mk3,\[Phi]table+h* \[Phi]k3,utable+h* uk3,\[Omega]BD, \[CapitalOmega]k,du,zLine[[i ]]];
\[Phi]k4=d\[Phi][\[Rho]mtable+h*\[Rho]mk3,\[Phi]table+h* \[Phi]k3,utable+h* uk3,zLine[[i ]]];
\[Rho]mk4=d\[Rho]m[\[Rho]mtable+h*\[Rho]mk3,\[Phi]table+h* \[Phi]k3,utable+h* uk3,zLine[[i ]]];
uk4=du[\[Rho]mtable+h*\[Rho]mk3,\[Phi]table+h* \[Phi]k3,utable+h* uk3,\[Omega]BD, \[CapitalOmega]k,zLine[[i ]]];

Hval=Htable[[i-1]]+h (Hk1+2 Hk2+2 Hk3+Hk4)/6;

If[
Hval===Indeterminate,
containsIndeterminate = True;
Break[]
];

If[
!Element[Hval, Reals],
containsComplex = True;
Break[]
];
Htable[[i]] = Hval;
\[Rho]mtable=\[Rho]mtable+h (\[Rho]mk1+2 \[Rho]mk2+2 \[Rho]mk3+\[Rho]mk4)/6;
\[Phi]table=\[Phi]table+h (\[Phi]k1+2 \[Phi]k2+2 \[Phi]k3+\[Phi]k4)/6;
utable=utable+h (uk1+2 uk2+2 uk3+uk4)/6,

{i,2,Length@zLine }
];
{containsIndeterminate, containsComplex, Htable[[2;;]]}];

applyMCMC[variables_?ListQ, chi2old_, priors_?ListQ]:=Module[{\[CapitalOmega]m,\[CapitalOmega]k,H0,\[Phi]0,d\[Phi]0,\[Omega]BD, chi2prop, \[CapitalOmega]mNew,\[CapitalOmega]kNew,\[CapitalOmega]deNew,H0New,\[Phi]0New,d\[Phi]0New,\[Omega]BDNew, containsIndeterminate, containsComplex, vals, jumpprob, newVariables},
(* Propose a new point based on the jump function with different standard dev*)

{\[CapitalOmega]mNew,\[CapitalOmega]kNew,H0New,\[Phi]0New,d\[Phi]0New,\[Omega]BDNew} = variables +( RandomVariate[NormalDistribution[0,#]]&/@{0.01,0.001, 0.1,0.1, 0.001, 10 });

(*{\[CapitalOmega]mNew,\[CapitalOmega]kNew,H0New,\[Phi]0New,d\[Phi]0New,\[Omega]BDNew} = RandomVariate[MultinormalDistribution[variables,DiagonalMatrix[{0.01,0.01, 0.1,0.01, 0.001, 1.5 }]^2]]*);

\[CapitalOmega]deNew=1 - \[CapitalOmega]mNew - \[CapitalOmega]kNew;

{containsIndeterminate, containsComplex, vals} = Quiet@RK4Method[dH,d\[Phi],d\[Rho]m,du,\[CapitalOmega]mNew,\[CapitalOmega]kNew,H0New,\[Phi]0New,d\[Phi]0New,\[Omega]BDNew,zLine];

newVariables = {};

If[!containsIndeterminate ,

If[!containsComplex ,
(* This is the \[Chi]^2 at the proposed point *)
chi2prop=chi2[vals];
(* This is the probability to jump to the new point=Exp[-\[Chi]^2(prop)]/Exp[-\[Chi]^2(current)]*)
jumpprob=Exp[-(chi2prop-chi2old)/2];
(*Print@jumpprob;*)
(* Decide whether to keep the new point taking into account also priors. If so, store it in the chain, otherwise start over. *)
If[RandomReal[]<=Min[{1,jumpprob}]&& And@@ IntervalMemberQ[Interval@#&/@priors,{\[CapitalOmega]mNew,\[CapitalOmega]kNew,H0New,\[Phi]0New,d\[Phi]0New,\[Omega]BDNew}],
newVariables = {\[CapitalOmega]mNew,\[CapitalOmega]kNew,\[CapitalOmega]deNew,H0New,\[Phi]0New,d\[Phi]0New,\[Omega]BDNew,chi2prop};
];
];
];
newVariables
];

dataH=Import["H_All.txt","Table"];
dataH = DeleteCases[dataH, x_?(Length[#]==0&), 1];
ndata=Length[dataH];
zLine = dataH[[All, 1]];

(* Choose the total number of steps *)
totalsteps=30000;
(* Choose the initial point for BD osmology as {\[CapitalOmega]m,\[CapitalOmega]de,\[CapitalOmega]k=1.0-\[CapitalOmega]m-\[CapitalOmega]de,H0,\[Phi]0,d\[Phi]0,\[Omega]BD}*)
\[CapitalOmega]m=0.275;
\[CapitalOmega]k=-0.005;
\[CapitalOmega]de=1-\[CapitalOmega]m - \[CapitalOmega]k;
H0=71.0;
\[Phi]0=2.5;
d\[Phi]0=0.005;
\[Omega]BD=77.75;
burnin=500;
priors={{0.24,0.4},{-0.01,0.01}, {64,76}, {0,4}, {0., 0.01}, {-350,350}};
variables = {\[CapitalOmega]m,\[CapitalOmega]k,H0,\[Phi]0,d\[Phi]0,\[Omega]BD};
variableIdx = {1,2,4,5, 6, 7};
axesLabels = {"\[CapitalOmega]m","\[CapitalOmega]k","H0","\[Phi]0","d\[Phi]0","\[Omega]BD"};
chi2[h_?ListQ]:=Total[((dataH[[All,2]]-h)/dataH[[All,3]])^2];

(* This is the current point solution*)
{containsIndeterminate, containsComplex, vals}  = RK4Method[dH,d\[Phi],d\[Rho]m,du,\[CapitalOmega]m,\[CapitalOmega]k,H0,\[Phi]0,d\[Phi]0,\[Omega]BD,zLine];
(* This is the \[Chi]^2 at the current point*)
chi2old= chi2[vals];

(* Initialize the chain. This is were the points are stored *)
chain={};
chlen = 0;

(* Start the MCMC *)
(* Initialize the MCMC *)

step=0;
While[
step< totalsteps,
newVariables=Quiet@applyMCMC[variables, chi2old, priors];
If[Length@newVariables>1,
AppendTo[chain,newVariables];
variables = newVariables[[variableIdx]];
chi2old = newVariables[[Length@newVariables]];
chlen=Dimensions[chain][[1]];
];
step++
]

Export["chain.m",chain];
(* The values of \[CapitalDelta]\[Chi]^2 for 1,2,3\[Sigma] *)
dx2=N@Table[2 InverseGammaRegularized[1/2,1-Erf[x\[Sigma]/(\[Sqrt]2)]],{x\[Sigma],1,3}];
chi2List = chain[[All, Dimensions[chain][[2]]]];
Print[StringForm["Estimated Initial \!\(\*SuperscriptBox[\(\[Chi]\), \(2\)]\) is ``, \nEstimated Minimum \!\(\*SuperscriptBox[\(\[Chi]\), \(2\)]\) is ``",chi2List[[1]],Min@chi2List]];
(*sortedIdx = Reverse@Ordering[chi2List];
chain = chain[[sortedIdx, All]];*)
chi2min = Min@chi2List;
minimumChiPos = Flatten[ Position[chi2List, chi2min], 1]//First;
{\[CapitalOmega]m,\[CapitalOmega]k,\[CapitalOmega]de,H0,\[Phi]0,d\[Phi]0,\[Omega]BD,chi2current} = chain[[minimumChiPos, All]];
Print[StringForm["Parameters converged to the following values \n\[CapitalOmega]m -> ``,\n \[CapitalOmega]k -> ``,\n \[CapitalOmega]de -> ``,\n H0 -> ``,\n \[Phi]0 -> ``,\n d\[Phi]0 -> ``,\n \[Omega]BD -> ``",\[CapitalOmega]m,\[CapitalOmega]k,\[CapitalOmega]de,H0,\[Phi]0,d\[Phi]0,\[Omega]BD]];

chainShort = chain[[All,variableIdx ]];
varIntervals = MinMax/@Transpose@chainShort;

pl = PlotWithGrid[chainShort, {\[CapitalOmega]m,\[CapitalOmega]k,H0,\[Phi]0,d\[Phi]0,\[Omega]BD}, varIntervals, axesLabels, 2,chainShort[[minimumChiPos]], "scatter.png"];
(* Select the points in the 1,2,3 sigma region *)
data1s=Select[chain,#[[8]]<=(chi2min+dx2[[1]])&];
data2s=Select[chain,#[[8]]<=(chi2min+dx2[[2]])&];
data3s=Select[chain,#[[8]]<=(chi2min+dx2[[3]])&];
varIntervals = MinMax/@Transpose@data3s[[All, variableIdx]];

pl = PlotWithGrid[{data1s[[All, variableIdx]], data2s[[All, variableIdx]], data3s[[All, variableIdx]]}, {\[CapitalOmega]m,\[CapitalOmega]k,H0,\[Phi]0,d\[Phi]0,\[Omega]BD}, varIntervals, axesLabels, 3,chainShort[[minimumChiPos]], "contour_1.png"];

**Results : scatter.png  and contour_1.png are empty.**

The file "H_All.txt" is below :

0.0700 69.00 19.6
0.0900 69.00 12.0
0.1200 68.60 26.2
0.1700 83.00 8.00
0.1790 75.00 4.00
0.1990 75.00 5.00
0.2000 72.90 29.6
0.2400 79.69 2.65
0.2700 77.00 14.0
0.2800 88.80 36.6
0.3500 84.40 7.00
0.3520 83.00 14.0
0.3800 83.00 13.5
0.4000 95.00 17.0
0.4004 77.00 10.2
0.4247 87.10 11.2
0.4300 86.45 3.68
0.4400 82.60 7.80
0.4449 92.80 12.9
0.4783 80.90 9.00
0.4800 97.00 62.0
0.5700 92.40 4.50
0.5930 104.0 13.0
0.6000 87.90 6.10
0.6800 92.00 8.00
0.7300 97.30 7.00
0.7810 105.0 12.0
0.8750 125.0 17.0
0.8800 90.00 40.0
0.9000 117.0 23.0
1.0370 154.0 20.0
1.3000 168.0 17.0
1.3630 160.0 33.6
1.4300 177.0 18.0
1.5300 140.0 14.0
1.7500 202.0 40.0
1.9650 186.5 50.4
2.3000 224.0 8.00
2.3400 222.0 7.00
2.3600 226.0 8.00
    

and the error log I get during a direct launching by command line :

$ wolframscript BD_MCMC_2_short.wls

StringForm[Estimated Initial DisplayForm[SuperscriptBox[χ, 2]] is ``,
Estimated Minimum DisplayForm[SuperscriptBox[χ, 2]] is ``, 220.3044990025602, 21.58177375066642]
StringForm[Parameters converged to the following values
Ωm -> ``,
 Ωk -> ``,
 Ωde -> ``,
 H0 -> ``,
 ϕ0 -> ``,
 dϕ0 -> ``,
 ωBD -> ``, 0.24041269047499136, 0.0016398283879296882, 0.757947481137079, 64.06857800452705, 1.579419511447762, 0.00011446222576435193, 334.91426614586453]

FrontEndObject::notavail: A front end is not available; certain operations require a front end.

FrontEndObject::notavail: A front end is not available; certain operations require a front end.

FrontEndObject::notavail: A front end is not available; certain operations require a front end.

General::stop: Further output of FrontEndObject::notavail will be suppressed during this calculation.

ReplaceAll::reps: {$Failed} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

Symbol::argx: Symbol called with 0 arguments; 1 argument is expected.

Symbol::argx: Symbol called with 0 arguments; 1 argument is expected.

Function::slotn: Slot number 2 in {Reverse[MinMax[#1]], MinMax[#2]} &  cannot be filled from ({Reverse[MinMax[#1]], MinMax[#2]} & )[Missing[String[], Symbol[]]].

Function::slotn: Slot number 2 in {Reverse[MinMax[#1]], MinMax[#2]} &  cannot be filled from ({Reverse[MinMax[#1]], MinMax[#2]} & )[String[]].

Function::slotn: Slot number 2 in {Reverse[MinMax[#1]], MinMax[#2]} &  cannot be filled from ({Reverse[MinMax[#1]], MinMax[#2]} & )[Symbol[]].

General::stop: Further output of Function::slotn will be suppressed during this calculation.

Outer::heads: Heads List and Missing at positions 3 and 2 are expected to be the same.

...                                                                                                                                                                                                                Total[Transpose[Mean[MapThread[Switch[(Abs @* Subtract)[1], Inherited, (Abs @* Subtract)[2], Min, (Abs @* Subtract)[#3, 1], Max, (Abs @* Subtract)[#3, 2], (Abs @* Subtract)[], (Abs @* Subtract)[1]] & , {{{{Abs[Subtract[Max, Max, Max, Max, Max]], Abs[Subtract[Max, Max, Max, Max, Max]]}, {Abs[Subtract[Max, Max, Max, Max, Max]], Abs[Subtract[Max, Max, Max, Max, Max]]}}, {{Abs[Subtract[Max, Max, Max, Max, Max]], Abs[Subtract[Max, Max, Max, Max, Max]]}, {Abs[Subtract[Max, Max, Max, Max, Max]], Abs[Subtract[Max, Max, Max, Max, Max]]}}, {{Abs[Subtract[Max, Max, Max, Max, Max]], Abs[Subtract[Max, Max, Max, Max, Max]]}, {Abs[Subtract[Max, Max, Max, Max, Max]], Abs[Subtract[Max, Max, Max, Max, Max]]}}, {{Abs[Subtract[Max, Max, Max, Max, Max]], Abs[Subtract[Max, Max, Max, Max, Max]]}, {Abs[Subtract[Max, Max, Max, Max, Max]], Abs[Subtract[Max, Max, Max, Max, Max]]}}, {{Abs[Subtract[Max, Max, Max, Max, Max]], Abs[Subtract[Max, Max, Max, Max, Max]]}, {Abs[Subtract[Max, Max, Max, Max, Max]], Abs[Subtract[Max, Max, Max, Max, Max]]}}}, Missing[NotAvailable, PlotRange], {{{Abs[Subtract[#2]], Abs[Subtract[#1]]}[{Abs[String[] - Symbol[]], Abs[String[] - Symbol[]]}, 0], {Abs[Subtract[#2]], Abs[Subtract[#1]]}[{Abs[String[] - Symbol[]], Abs[String[] - Symbol[]]}, 0]}, {{Abs[Subtract[#2]], Abs[Subtract[#1]]}[{0, 0}, 0], {Abs[Subtract[#2]], Abs[Subtract[#1]]}[{0, 0}, 0]}}}, 4]]]]   Mean[1] + Mean[2] + Mean[3]  Mean[1] + Mean[2] + Mean[3]  Mean[1] + Mean[2] + Mean[3]              Total[Transpose[Mean[Missing[NotAvailable, ImageSize]]]]   Mean[1] + Mean[2] + Mean[3]  Mean[1] + Mean[2] + Mean[3]  Mean[1] + Mean[2] + Mean[3]
    
    Transpose::argt: Transpose called with 0 arguments; 1 or 2 arguments are expected.
    
    General::stop: Further output of Transpose::argt will be suppressed during this calculation.
    
    MapThread::list: List expected at position 2 in MapThread[If[#2 === Null, ##1] & , Missing[PartAbsent, 2], 2].
    
    Divide::argr: Divide called with 1 argument; 2 arguments are expected.
    
    ^C
    Interrupt> Your options are:
        abort (or a) to abort current calculation
        continue (or c) to continue
        exit (or quit) to exit Mathematica
        inspect (or i) to enter an interactive dialog

EDIT 2: Here a suggestion from Wolfram support :

" The appropriate syntax would be

UsingFrontEnd[pl=ResourceFunction["PlotGrid"][...]]

However, please note that ResourceFunction items are not part of the core Wolfram Language, and are instead created by community members. As such, we cannot provide assistance or guidance with any specific ResourceFunction, and any issues regarding PlotGrid not working with wolframscript or when run in certain ways should be directed to the function's authors via the feedback form "

So, i tried to do :

UsingFrontEnd[pl=ResourceFunction["PlotGrid"][
tb,
PlotRange-> Max,
ImageSize-> Large,
Method->{"AllCustomTicks"->Automatic, "FixFrameTicks"->True}
]];

But get still errors during execution of Wolframscript .wls. Any help is welcome.

Answer 1

Ok, finally managed to get it to run from within wolframscript...

I have created an updated version of PlotGrid that should work properly in non-standard environments such as wolframscript. It is not yet published, but you can get it as ResourceFunction@CloudObject["https://www.wolframcloud.com/obj/langl/DeployedResources/Function/PlotGrid/ResourceObject.wl"] in the meantime.

With this, your script produces plots as expected:

(* rest of script ... *)

pl=ResourceFunction[CloudObject["https://www.wolframcloud.com/obj/langl/DeployedResources/Function/PlotGrid/ResourceObject.wl"]][
   tb,
   PlotRange-> Max,
   ImageSize-> Large,
   Method->{"AllCustomTicks"->Automatic, "FixFrameTicks"->True}
];

(* rest of script ... *)

The resulting contour_1.png and scatter.png look like this for me:

What went wrong

For those who are curious: There are a few things that went wrong in the original version of PlotGrid

• UsingFrontEnd was required since PlotGrid uses CallFrontEnd internally to obtain the dimensions of the plots (see also this answer by @CarlWoll). I have now included UsingFrontEnd, so that it automatically starts a frond-end if required.
• The main issues all come from the fact that the front-end launched by UsingFrontEnd in turn launches another kernel that is used to evaluate GridLines, Ticks and FrameTicks. This has two consequences:

  • Using custom ticks functions doesn't work, since the other kernel doesn't know about them. This is fixed rather easily by transporting the relevant definitions to the other kernel using Language`ExtendedFullDefinition inside a FrameTicks specification.

  • Any tracking injected into these options can't communicate the results back to the main kernel. To complicate things further, the other kernel is launched in a sandboxed mode, heavily restricting the communication channels. The best solution I have found is a trick based on ExportPacket:

    FrontEndExecute@FrontEnd`ExportPacket[
      Cell@BoxData@DynamicBox@ToBoxes[expr],
      "PlainText"
    ]
    

    This evaluates expr on the other kernel and returns the result as a string, which we can then interpret using ToExpression. Starting in version 13.0, this communication is no longer required, since we can get all the necessary information from the newly-fixed AbsoluteOptions. The definitions however still need to be transported over.

In total, the strategy is thus as follows:

1. "Rasterize" a Graphics expression with FrameTicks->{{None,None},{None,(Language`ExtendedFullDefinition[]=def)&}. This transports the relevant definitions to the other kernel. Here, "rasterize" is effectively equivalent to Rasterize[...,"Regions"].
2. "Rasterize"/Measure the actual Graphics expressions. Here, we save some data about the plot range in the DownValues of SowData on the helper kernel. Originally, this was using Sow, but since this stuff is now evaluated on the other kernel, it didn't make its way back to the main kernel.
3. Use the ExportPacket trick to get the DownValues of SowData to the main kernel, and evaluate them, effectively replaying the Sow calls.