Parameter space search routine is too fast?

I wrote a small script to comb through the parameter space of a function looking for where the function satisfies a certain relation everywhere on its (co-ordinate space) domain. However, since I am quantizing the parameter space, I end up with a lattice with something like 41 trillion lattice points. Moreover, at each point of this lattice, I am checking a condition on roughly 500 points of the co-ordinate space. So, in actuality, Mathematica should be looking through something like 20 quadrillion individual points.

So I would expect this to take a long time! However, when I execute the code, it finishes in something like 12 seconds! I've tested the script on much simpler functions and a much smaller parameter space and it behaves exactly as I expect it to. So unless Mathematica is extremely efficient or I'm a coding genius, I cant see how it completes the script this fast.

I incorporated parallel processing into the main search-routine of the script, but each kernel should still have to comb through basically quadrillion lattice points. Now, this number of lattice points is on the larger side than most reasonable functions I mess with, but its still a parameter space size I need to comb through.

Setup

I'll lay out the relevant parts of the code below: Heres the function I'm feeding through the script. I wont paste the function itself since its massive, but I will paste the code that generates the function:

\[Phi] = z (Sech[2*(x^2 + y^2 + (z)^2 + 1)] -
Sech[2*(x^2 + y^2 + (z)^2)]);
expr = e*(z - \[Alpha])*(
Tanh[s*(x^2 + y^2 + (z - \[Alpha])^2 + l)] -
Tanh[s*(x^2 + y^2 + (z - \[Alpha])^2 - l)])/(
2*Tanh[s*l]) {-y, x, 0} +
f*(y - \[Beta]) (
Tanh[r*(x^2 + (y - \[Beta])^2 + (z)^2 + m)] -
Tanh[r*(x^2 + (y - \[Beta])^2 + (z)^2 - m)])/(
2*Tanh[r*m]) {0, -z, y} +
g*(x - \[Chi])*(
Tanh[t*((x - \[Chi])^2 + (y)^2 + (z)^2 + n)] -
Tanh[t*((x - \[Chi])^2 + (y)^2 + (z)^2 - n)])/(
2*Tanh[t*n]) {z, 0, -x};
Curlexpr = Curl[expr, {x, y, z}];

func = (-Derivative[0, 0, 0, 1][Subscript[N, 3]][t, x, y, z]^2 - Derivative[0, 0, 1, 0][Subscript[N, 2]][t, x, y, z]^2 -
(Derivative[0, 0, 0, 1][Subscript[N, 2]][t, x, y, z] + Derivative[0, 0, 1, 0][Subscript[N, 3]][t, x, y, z])^2/2 -
Derivative[0, 1, 0, 0][Subscript[N, 1]][t, x, y, z]^2 + (Derivative[0, 0, 0, 1][Subscript[N, 3]][t, x, y, z] +
Derivative[0, 0, 1, 0][Subscript[N, 2]][t, x, y, z] + Derivative[0, 1, 0, 0][Subscript[N, 1]][t, x, y, z])^2 -
(Derivative[0, 0, 1, 0][Subscript[N, 1]][t, x, y, z] + Derivative[0, 1, 0, 0][Subscript[N, 2]][t, x, y, z])^2/2 -
(Derivative[0, 0, 0, 1][Subscript[N, 1]][t, x, y, z] + Derivative[0, 1, 0, 0][Subscript[N, 3]][t, x, y, z])^2/2)/.Table[Subscript[N, i] ->
Evaluate@Function[{t, x, y, z},
Evaluate@(D[\[Phi], {{x, y, z}, 1}] + Curlexpr)[[i]]], {i, 1, 3}]



so func is the function I am evaluating on the quantized spaces.

Heres the setup of the script.

function = Func; (*test function:   x*Exp[-x^2/\[Sigma]]-s;*)

Quantifier = function >= 0; (*what condition are we trying to satisfy over the entire domain*)

variables = {x, y, z};

Complement[Variables@Level[function, -1], variables];

Parameters = {e, f, g, l, m, n, r, s, t, \[Alpha], \[Beta], \[Chi]};

(*Complement[Variables@Level[function,-1], variables];
(*ORDERED*)*)(*extract the parameters from the test function, not including the co-ordinate variables*)

ParameterSpaceSizes = {{-3, 3}, {-3, 3}, {-3, 3}, {-3, 3}, {-3, 3}, {-3, 3}, {0.1, 4}, {0.1,4}, {0.1, 4}, {-1, 1}, {-1, 1}, {-1,1}}; (*ORDERED*)

CoordinateSpaceBounds = {{-2, 2}, {-2, 2}, {-2, 2}};(*ORDERED*)

ParameterSpaceResolution = 0.3; (* How accurate should the parameter search be*)

CoordinateSpaceResolution = 0.5;


The quantifier is the relation I want the function to satisfy on its entire domain. The parameterspacesizes are the ranges of the parameters, same for the co-ords. the resolutions are the lattice spacing for the parameter/co-ordinate space.

Functions

Here are some relevant functions. Meshify generates a lattice grid, given the input intervals as a list, e.g. {{a,b}, {c,d}} would represent the intervals [a,b] and [c,d], and the lattice spacing, also as a list.

Meshify[ParameterSpaceIntervals_, ParameterSpaceResolutions_]:=
(*Discretize the parameter space, specified by bounds of the individual parameters, with a given resolution*)

Module[{
ListofDiscretizedIntervals = Array[(Array[(N@#&), Round@((ParameterSpaceIntervals[[All,2]][[#1]]-ParameterSpaceIntervals[[All,1]][[#1]])/ParameterSpaceResolutions[[#1]]+1), {ParameterSpaceIntervals[[All,1]][[#1]], ParameterSpaceIntervals[[All,2]][[#1]]}]&), Length@ParameterSpaceResolutions] (*discretize the intervals and store the output in the temporary variable [ListofDiscretizedIntervals]*)
},

Return[Tuples@ListofDiscretizedIntervals, Module];
(*Take the cartesian product of the discretized intervals to get a discretized parameter space and return the output as the return value of the Module function*)
]


nPartitions just partitions a set into n partitions in a maximally even way, i.e. all partitions are roughly the same size.

nPartitions[set_,
n_] := (*Partition the input set in n maximally-even partitions*)

Module[{
residual = Mod[Length@set, n], LenSet = Length@set
},
If[residual != 0,
ReducedSet =
set[[1 ;;
LenSet -
residual]]; (*If the number of partitions doesnt divide the \
norm of the set,
then we partition a reduced set into n partitions and append the \
remaining elements to the nth partition in an even way,
thus creating a maximally even partition. Note,
the set wont be maximally evenly partitioned if the number of \
partitions is larger than half the norm of the set.*)

ReducedPartition = Partition[ReducedSet, Length@ReducedSet/n]; (*
define the reduced partition*)

Do[AppendTo[ReducedPartition[[-residual + ii]],
ReducedPartition[[-residual + ii + 1]][[1]]], {ii, 0,
residual - 2}];(*
Move over the first elements of one partition to the end of there \
left neighbor*)

Do[ReducedPartition =
Delete[ReducedPartition, {-residual + ii, 1}], {ii, 1,
residual -
1}];(*remove the shifted over elements from their previous \
position*)

Do[AppendTo[ReducedPartition[[-1]],
set[[ii]]], {ii, -residual, -1}]; (*append the residual elements \
to the last partition*)

Return[ReducedPartition,
Module](*return the maximally even partitioned set*)
,
Return[Partition[set, Length@set/n],
Module] (*if the number of partitions divides the norm of the \
set, then we can partition evenly*)
]
]


PartitionIntervals takes in a set of intervals and slices each interval in n pieces, then combines the i'th pieces together as a single partition.

PartitionIntervals[Intervals_,
nPartitions_] :=(* Partitions a set of intervals of the form \
{{a,b}, {c,d},{e,f},...} into nPartitions partitions*)
Array[
(Table[{Intervals[[ii]][[
1]] + (# - 1)*(Intervals[[ii]][[2]] - Intervals[[ii]][[1]])/
nPartitions,
Intervals[[ii]][[
1]] + #*(Intervals[[ii]][[2]] - Intervals[[ii]][[1]])/
nPartitions}, {ii, 1, Length@Intervals}] &), nPartitions
]


the scanObjectUntilCondition function is really the real meat and potatoes of this search routine. It takes in a quantifier as a set of relations (think 'for all x, f(x)>3' would translate to a lattice as, 'at each lattice point, f(lattice point)>3', so each position of the set corresponds to each lattice point), the parameters of the function, the intervals of values of the parameters, and the lattice spacing. It then loops through the number of values a single parameter can take (since the lattice spacing is constant, its just 1/lattice spacing, i.e. the number of lattice points on any edge of the lattice), for each pass of this loop, it loops through the number of lattice points on an edge (the number of possible values a single parameter can take), for each pass of THIS inner loop, the parameter values are plugged into the quantifier set. If any of the relations in the quantifier set evaluate to false, (i.e. the quantifier is not satisfied), it increments the next parameter by 1 lattice spacing in the next pass of the inner loop. (i.e. the inner loop proceeds as (1,1,1)-> (2,1,1) -> (2,2,1) -> (2,2,2), then the next iteration of the outer loop would go (3,2,2)-> (3,3,2)-> (3,3,3), etc. until we've scanned through the entire parameter space)

Thats a lot of words to basically say the function scans through the quantized parameter space looking for points where the quantifier is satisfied.

scanObjectUntilCondition[Quantifier_, params_, intervals_,
ress_] := (*Scan the parameters params through the intervals \
intervals with resolution ress until quantifier is satisfied at every \
element*)Module[{paramVals = intervals[[All, 1]],
EndPoints = intervals[[All, 2]],
StartPoints = intervals[[All, 1]], con = Quantifier},
Result = Catch[
For[\[Lambda]\[Lambda]\[Lambda] =
1, \[Lambda]\[Lambda]\[Lambda] <=
Round@(1/
ress), \[Lambda]\[Lambda]\[Lambda]++,(*Loop over the bins,
given by 1/ress*)

Array[(*loop over the parameters*)
(
If[
AllTrue[(*If the Quantifier is satisfied at each of its \
elements, then immediately return the point where this is achieved*)

con /. Thread[params -> paramVals], TrueQ],
Throw[{"Condition met at point: ",
paramVals =
ReplacePart[
paramVals, # ->
paramVals[[#]] + (EndPoints[[#]] - StartPoints[[#]])*

ress](*If the quantifier contains a point where it \
evaluates to False,
then increment the parameter values by one bin width and \
continue searching*)
]
) &, Length@intervals]
], o (*Catch Throw Tag*)
];
If[TrueQ[Result[[0]] == List], Return[Result, Module],
Return[$Failed, Module]] ]  Parallelization Now I set up the parameters of the parallelization routine. the QuantifierOverCoordSpace is the variable that evaluates the function over the quantized co-ordinate space, so then the parameters are the only free variables. This set is the quantized version of "For all x, f(x)>0" captured in a single variable. (*construct the discretized co-ordinate space and extract the number \ of CPU cores to run on *) NumParams = Length@Parameters; (*number of parameters*) NumParamSpacePartitions \ =$ProcessorCount; (*how many partitions should we make*)

DiscretizedCoordinateSpace =
Meshify[CoordinateSpaceBounds,
ConstantArray[CoordinateSpaceResolution,
Length@CoordinateSpaceBounds]];
PartitionParameterIntervals =
PartitionIntervals[ParameterSpaceSizes, NumParamSpacePartitions];

(*Setup parallelization*)

Off[LaunchKernels::nodef]
LaunchKernels[]; (*make sure multiple kernels are running *)
On[
LaunchKernels::nodef]
QuantifierOverCoordSpace =
ParallelMap[Quantifier /. Thread[variables -> #] &,
DiscretizedCoordinateSpace];
DistributeDefinitions[Off[General::munfl], Off[GreaterEqual::nord],
Parameters, PartitionParameterIntervals, ParameterSpaceResolution,
QuantifierOverCoordSpace, scanObjectUntilCondition];


And heres the part of the script that executes the parallelized search routine

Print["Executing over ", Length@Kernels[], " Kernels"]
ClearAll[result];

ParallelTry[
(result =
scanObjectUntilCondition[QuantifierOverCoordSpace, Parameters,
PartitionParameterIntervals[[#]], ParameterSpaceResolution];

Update 2 (slightly less problems)

So, per my comment below, the temperamental problem was an issue of critical sections. The If statement in the ParallelTry argument is actually redundant since my function scanObjectUntilCondition already returns $Failed if it doesnt find anything. A dumb mistake, I should probably know what my own code does! Nonetheless, that fixed the issue. So now it returns the expected result of the simple function above everytime. But it still does it WAY faster than I expect it too. A single subkernel is still scanning through roughly 500 trillion points in about 0.1 seconds. My CPU runs at 3.7 GHZ so a single core should realistically only be able to search 370 Million points in that time, several orders of magnitude lower than what it claims to be doing. I still cant figure out why, but I need to know so I can trust the future results on much larger and more important functions. Heres the updated code that runs the above sample functions. I'll leave the original code in this post in case it helps someone else find their problem.  function =x*Exp[-x^2/(\[Sigma] + \[Tau] + \[Upsilon] + \[Zeta])] + s; Quantifier = function >= 0; (*what condition are we trying to satisfy over the entire domain*) variables = {x}; Complement[Variables@Level[function, -1], variables]; Parameters = {s, \[Sigma], \[Tau], \[Upsilon], \[Zeta]}; (*Complement[Variables@Level[function,-1], variables]; \ (*ORDERED*)*)(*extract the parameters from the test function, not \ including the co-ordinate variables*) ParameterSpaceSizes = {{-10,2}, {0.1, 5}, {0.1, 5}, {0.1, 5}, {0.1,5}};(*ORDERED*)(* s can range from -5 to 5 and \[Sigma] and run \ from 0.001 to 5*) CoordinateSpaceBounds = {{-2, 2}}; {{-2, 2}, {-2,2}, {-2, 2}};(*ORDERED*) ParameterSpaceResolution = 0.01; (* How accurate should the parameter search be*) CoordinateSpaceResolution = 0.1; (*Some functions to use in setting up and discretizing the parameter space and coordinate space*) Meshify[ParameterSpaceIntervals_, ParameterSpaceResolutions_] := (*Discretize the parameter space, specified by bounds of the individual parameters, with a given resolution*) Module[{ListofDiscretizedIntervals = Array[(Array[(N@# &), Round@((ParameterSpaceIntervals[[All, 2]][[#1]] - ParameterSpaceIntervals[[All, 1]][[#1]]) ParameterSpaceResolutions[[#1]] + 1), {ParameterSpaceIntervals[[All, 1]][[#1]], ParameterSpaceIntervals[[All, 2]][[#1]]}] &), Length@ParameterSpaceResolutions] (*discretize the intervals and store the output in the temporary variable [ ListofDiscretizedIntervals]*) }, Return[Tuples@ListofDiscretizedIntervals, Module]; (*Take the cartesian product of the discretized intervals to get a discretized parameter space and return the output as the return value of the Module function*) ] nPartitions[set_, n_] := (*Partition the input set in n maximally-even partitions*) Module[{residual = Mod[Length@set, n], LenSet = Length@set}, If[residual != 0,ReducedSet = set[[1 ;; LenSet - residual]]; (*If the number of partitions doesnt divide the norm of the set, then we partition a reduced set into n partitions and append the remaining elements to the last few partitions in an even way, thus creating a maximally even partition. Note, the set wont be maximally evenly partitioned if the number of partitions is larger than half the norm of the set.*) ReducedPartition = Partition[ReducedSet, Length@ReducedSet/n]; (* define the reduced partition*) Do[AppendTo[ReducedPartition[[-residual + ii]], ReducedPartition[[-residual + ii + 1]][[1]]], {ii, 0, residual - 2}];(* Move over the first elements of one partition to the end of there left neighbor*) Do[ReducedPartition = Delete[ReducedPartition, {-residual + ii, 1}], {ii, 1, residual - 1}];(*remove the shifted over elements from their previous position*) Do[AppendTo[ReducedPartition[[-1]], set[[ii]]], {ii, -residual, -1}]; (*append the residual elements to the last partition*) Return[ReducedPartition, Module](*return the maximally even partitioned set*), Return[Partition[set, Length@set/n], Module] (*if the number of partitions divides the norm of the set, then we can partition evenly*) ] ] PartitionIntervals[Intervals_, nPartitions_] :=(* Partitions a set of intervals of the form {{a,b}, {c,d},{e,f},...} into nPartitions partitions*) Array[ (Table[{Intervals[[ii]][[ 1]] + (# - 1)*(Intervals[[ii]][[2]] - Intervals[[ii]][[1]]) nPartitions, Intervals[[ii]][[1]] + #*(Intervals[[ii]][[2]] - Intervals[[ii]][[1]])/nPartitions}, {ii, 1, Length@Intervals}] &), nPartitions ] scanObjectUntilCondition[Quantifier_, params_, intervals_, ress_] := (*Scan the parameters params through the intervals intervals with resolution ress until quantifier is satisfied at every element*) Module[{paramVals = intervals[[All, 1]], EndPoints = intervals[[All, 2]], StartPoints = intervals[[All, 1]], con = Quantifier}, Result = Check[ Catch[ For[\[Lambda]\[Lambda]\[Lambda] = 1, \[Lambda]\[Lambda]\[Lambda] <= Round@(1/ress), \[Lambda]\[Lambda]\[Lambda]++,(*Loop over the bins, given by 1/ress*) Array[(*loop over the parameters*) ( If[ AllTrue[(*If the Quantifier is satisfied at each of its elements, then immediately return the point where this is achieved*) con /. Thread[params -> paramVals], TrueQ], Throw[{"Condition met at point: ", Thread[params -> paramVals]}, o], paramVals = ReplacePart[paramVals, # -> paramVals[[#]] + (EndPoints[[#]] - StartPoints[[#]])*ress](*If the quantifier contains a point where it evaluates to False, then increment the parameter values by one bin width and continue searching*) ] ) &, Length@intervals] ], o (*Catch Throw Tag*) ], err ]; If[TrueQ[Result[[0]] == List], Return[Result, Module], Return[$Failed, Module]]
]

(*construct the discretized co-ordinate space and extract the number of CPU cores to run on *)

NumParams = Length@Parameters; (*number of parameters*)

NumParamSpacePartitions = $ProcessorCount; (*how many partitions should we make*) DiscretizedCoordinateSpace = Meshify[CoordinateSpaceBounds, ConstantArray[CoordinateSpaceResolution, Length@CoordinateSpaceBounds]]; PartitionParameterIntervals = PartitionIntervals[ParameterSpaceSizes, NumParamSpacePartitions]; (*Setup parallelization*) Off[LaunchKernels::nodef] LaunchKernels[]; (*make sure multiple kernels are running *) On[ LaunchKernels::nodef] QuantifierOverCoordSpace = ParallelMap[Quantifier /. Thread[variables -> #] &, DiscretizedCoordinateSpace]; DistributeDefinitions[Parameters, PartitionParameterIntervals, ParameterSpaceResolution, QuantifierOverCoordSpace, scanObjectUntilCondition]; Print["Executing over ", NumParamSpacePartitions, " Kernels"] (*Run the parallelized search routine*) ParallelizedResult = ParallelTry[ (scanObjectUntilCondition[QuantifierOverCoordSpace, Parameters, PartitionParameterIntervals[[#]], ParameterSpaceResolution]) &, Range[NumParamSpacePartitions] ] // AbsoluteTiming Print["Times in seconds for the parameter space to be searched with \ given resolution and bounds: \n\t\t", ParallelizedResult[[1]], "\nResult of search: \n\t\t", ParallelizedResult[[2]]] $$$$  • I see that you turned various messages off in your code. Try turning them back on: any error messages might hint at the reason for early termination. Other things that stand out are the extensive use of Subscript which may lead to subtle and sneaky malfunctions, and the use of AppendTo, which is often awful for performance (standard advice would be to consider Sow / Reap instead, but those have their own problems in parallel evaluation). Finally, what is your end game? Do you want to fix this code and then run it for a month or longer on your machine or find an alternative? – MarcoB Nov 30 '20 at 17:14 • Ahh,, Sorry I forgot to remove those from the pasted code. The code on my machine currently has those messages turned on. Those messages are because the search routine hits a certain value of the parameter that make the function undefined, like Coth[0]. The end game is to run the search routine on more reasonablly sized spaces, definitely not a 40 quadrillion point lattice, but more like tens of trillions, or trillions. I just need help finding why the routine is cutting itself short without an abort message like TimeConstrain or something of the sorts, that way I have more confidence in it. – shanedrum Nov 30 '20 at 17:21 • Oh interesting, I didnt know that about Subscript. That's worth checking into, thanks! Also, ya with AppendTo, I definitely am not using it on humongous sets, just on sets the size of the number of parameters. Moreover, the subkernels never call that function, only the master kernel outside of the search routine. Its only used when setting up everything. – shanedrum Nov 30 '20 at 17:27 • So I checked my code and Subscript is only used when calculating the function that goes into the routine. There's no instance of Subscript that the search routine sees at any level. – shanedrum Nov 30 '20 at 17:40 • Perhaps I should be using CriticalSections? It could be that a subkernel is overwriting the result variable in the ParallelTry argument before it can test if is in fact a result. – shanedrum Dec 1 '20 at 12:02 1 Answer Ahh, I figured it out. It is because the function scanObjectUntilCondition doesnt actually increment the parameters by one lattice spacing. Instead, It divides the length of the interval for each parameter by the "lattice spacing", i.e. the effective resolution size, and increments the parameter by this value. So this value is the actual bin width for each individual interval for the parameter. In other words, the interation goes like: if a is defined on the interval $$[-1,1]$$ and I specify a resolution size of 0.1, then a will be incremented like $$a = -1$$$$\downarrow$$$$-1 + (1-(-1))*0.1 =-.8$$$$\downarrow$$$$-0.8 + (1 - (-1))*0.1 = -0.6$$$$\downarrow$$$$etc.$$ This means the search routine isn't actually searching through quadrillions of lattice points, but a much smaller subset. Its funny, I knew the function would increment in this way when I coded it but just forgot about it when calculating the number of lattice points it was supposedly scanning over. Combined with the CriticalSection issue in my last update, my issues are fixed. Heres the updated scanObjectUntilCondition function if anyone would find a use for it :) scanObjectUntilCondition[Quantifier_, params_, intervals_, ress_] := (*Scan the parameters params through the intervals \ intervals with resolution ress until quantifier is satisfied at every \ element*)Module[{paramVals = intervals[[All, 1]], EndPoints = intervals[[All, 2]], StartPoints = intervals[[All, 1]], con = Quantifier, nPointsSearched = 0}, Result = Check[ Catch[ For[\[Lambda]\[Lambda]\[Lambda] = 1, \[Lambda]\[Lambda]\[Lambda] <= Round@(Max[EndPoints - StartPoints]/ ress), \[Lambda]\[Lambda]\[Lambda]++,(*Loop over the bins, given by 1/ress*) Array[(*loop over the parameters*) (If[ AllTrue[con /. Thread[params -> paramVals], TrueQ],(*If the Quantifier is satisfied at each of its \ elements, then immediately return the point where this is achieved*) Throw["Condition met at point: \n\t\t" <> ToString@Thread[params -> paramVals] <> "\n Number of Points Searched:\n\t\t" <> ToString@ nPointsSearched, o], If[\[Not] (paramVals[[#]] + ress > EndPoints[[#]]), paramVals = ReplacePart[paramVals, # -> paramVals[[#]] + ress]; nPointsSearched += 1 ](*If the quantifier contains a point where it evaluates \ to False, then increment the parameter values by one bin width and \ continue searching*) ] ) &, Length@intervals] ], o (*Catch Throw Tag*) ], err ]; If[TrueQ[Result[[0]] == String], Return[Result, Module], Return[$Failed, Module]]
]
`

Note, this now increments each parameter by the resolution size, so the lattice is no longer square but the true lattice spacing is now given by the user-defined resolution variable. Thats something to keep in mind when scanning over a large parameter space with a small resolution.

I basically fixed this issue myself, but I'll leave this post up in case it helps someone else out.