Parameter space search routine is too fast?

Question

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: ", 
           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*)
    ];
  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];
    If[TrueQ[Result[[0]] == List], result, $Failed]
    ) &, Range[NumParamSpacePartitions]] // Timing

Problem

Now, when I execute all the above code, the Timing function says the code completes in 7 seconds! But the thing that puzzles me greatly is that when I try a much simpler function with 1 or two parameters and 1 co-ordinate, it executes just as expected, I get back the first set of parameters that that satisfy the quantifier relation.

For example, if you change the function variable to something like function =x*Exp[-x^2/\[Sigma]] - s,

the variables variable to variables={x},

the Parameters variable to Parameters={s,\[sigma]},

the ParameterSpaceSizes variable to ParameterSpaceSizes={{-10,10}, {0.1, 5}},

the CoordinateSpaceBounds variable to CoordinateSpaceBounds={-5,5},

the code executes perfectly and gives me the parameters that satisfy the relation function\[greaterequal] 0 on the range {x,-5,5}.

So unless the subkernels are somehow aborting the evaluation when executing over the huge parameter space without notifying me, Mathematica is blistering fast, or my code is supremely efficient (not likely), I cant see how Mathematica finishes scanning over 12 quadrillion points in ~7 seconds over my 6 core CPU.

I dont think its a bug in my code since the script works perfectly well for much smaller parameter spaces and function, but its entirely possible. Perhaps Mathematica just kills the evaluation when it sees the numbers it has to comb through are ginormous?

Heres a small back of the envelope calculation I did just to get a sense of the order of magnitudes these numbers are:

My 6 core CPU typically will reach about 3.7 gigahertz when its working on the script. The script takes something like 60 seconds to complete when working over the largest parameter space I have used. If each core checks 1 lattice point every cycle, then after 60 seconds at 3 gigahertz, the CPU will have checked about 1-2 trillion lattice points. This is 40 thousand times smaller than the number of lattice points of the entire space! In order to check the entire space at this clock speed, it will have to run for 1 month! But its finishing in only 1 minute. what is going on?

Update (still problem)

So this is really interesting. If i use the following simple function

$$F(x) = x*e^{-\frac{x^2}{\sigma + \tau + \zeta + \Upsilon}} + s$$

and let the gaussian weights vary on some unimportant interval, and let s vary on $[-10, 2]$, then we expect the quantifier relation $\forall x$, $ F(x)\geq 0 $ to be satisfied at the parameter point $s=1.7$ for example. Now, I find the search routine is tempermental. Sometimes it will spit out a result, but other times it will produce $Failed which shouldnt happen. The result flips every time I run the code. Moreover, there are ~200 Billion lattice points to individually search through but the routine finishes in about 0.06 seconds, regardless if it found a point or not.

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]]]
```

Answer 1

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.