Here's an interesting question for the coders among us.
I have two flavors of a simulation I'm running. The algorithm is identical other than running one using a For Each() Loop and the other running Parallel Assign. The For Each loop consistently takes much less time to do the same calculation set as the Parallel Assign, even with all 8 of my cores maxed out.
Anybody have thoughts as to why that might be? Code for each version of the simulation is below:
Parallel Assign Version
// JSL Code for Simulation of Harmonized Compendial Dissolution Criteria USP <711>
// Simulation function for USP dissolution - three stage test
// reps = number of simulated lots for simulation
// mean = the process mean
// sd = the process standard deviation
// Q = value as specified in the mongraph for the product
//
// Summary Matrix Key
// ==================
// 1 "PassS1"
// 2 "PassS2"
// 3 "PassS3"
// 4 "#<Q S1"
// 5 "<Q-25 S1"
// 6 "Avg12"
// 7 "#<Q-15 S2"
// 8 "#<Q-25 S2"
// 9 "Avg24"
// 10 "#<Q-15 S3"
// 11 "#<Q-25 S3"
// 12 "PassAll"
// 13 "GoS2"
// 14 "GoS3"
// helper functions for the simulation
Clear Log();
vecSize = Function( {start, end, step},
{default local},
m = start :: end :: step;
steps = Length( m );
steps;
);
stopWatch = Function( {stTick, endTick, units = "sec"},
{default local},
seconds = endTick - stTick;
result = Match( units,
"sec", seconds,
"min", seconds / 60,
"hr", seconds / 3600,
"day", seconds / 86400
);
result = Round( result, 3 );
);
// Parallelized Simulation
// define the size of the matrix ans the step size
rows = vecSize( 0.1, 12, 0.1 );
rowShape = [0.1, 12, 0.1];
cols = vecSize( -1, 10, 0.5 );
colShape = [-1, 10, 0.5];
reps = 1000;
ntests = rows * cols;
//create the Parallel test matrix
simMat = J( rows, cols, Empty() );
// run the simulation
start = Tick Seconds();
Parallel Assign(
{rs = rowShape, cs = colShape, repN = reps, simDissoTest = Function(
{reps = 100, mean = 100, sd = 5, Q = 80},
{default local},
// generate data following user chosen mean and sd
// generate 24 data values even if not needed (that is
// not going to all three stages)
data_matrix = J( nrow = reps, ncol = 24, Random Normal( mean, sd ) );
// Summary Matrix to hold simulation results
sum_matrix = J( reps, 14, Empty() );
// Populate the summary field
For Each( {e, i}, 1 :: reps,
v_6 = data_matrix[i, 1 :: 6];
v_12 = data_matrix[i, 1 :: 12];
v_24 = data_matrix[i, 1 :: 24];
sum_matrix[i, 4] = Sum( v_6[Loc( v_6 < (Q + 5) )] );
sum_matrix[i, 5] = Sum( v_6[Loc( v_6 < (Q - 25) )] );
sum_matrix[i, 6] = Mean( v_12 );
sum_matrix[i, 7] = Sum( v_12[Loc( v_12 < (Q - 15) )] );
sum_matrix[i, 8] = Sum( v_12[Loc( v_12 < (Q - 25) )] );
sum_matrix[i, 9] = Mean( V_24 );
sum_matrix[i, 10] = Sum( v_24[Loc( v_24 < (Q - 15) )] );
sum_matrix[i, 11] = Sum( v_24[Loc( v_24 < (Q - 25) )] );
// Stage 1
If( sum_matrix[i, 4] == 0,
sum_matrix[i, 1] = 1;
sum_matrix[i, 12] = 1;
sum_matrix[i, 13] = 0;
sum_matrix[i, 14] = 0;
,
sum_matrix[i, 1] = 0;
// Check to see if any results are below Q-25 Stage 1
// If there are any - do not proceed - fail test
If( sum_matrix[i, 5] == 0,
sum_matrix[i, 13] = 1,
sum_matrix[i, 12] = 0;
sum_matrix[i, 13] = 0;
sum_matrix[i, 14] = 0;
);
);
// Stage 2
If( sum_matrix[i, 13] == 1,
// Check to see if any results are below Q-25 in Stage 2
// and there are more than 2 units that are less than Q-15
// If there are any - do not proceed - fail test
If( sum_matrix[i, 6] >= Q & sum_matrix[i, 7] == 0,
sum_matrix[i, 2] = 1;
sum_matrix[i, 12] = 1;
sum_matrix[i, 14] = 0;
,
sum_matrix[i, 2] = 0;
If( sum_matrix[i, 7] <= 2 & sum_matrix[i, 8] == 0,
sum_matrix[i, 14] = 1,
sum_matrix[i, 14] = 0;
sum_matrix[i, 12] = 0;
);
)
);
// Stage 3
If( sum_matrix[i, 14] == 1,
If( sum_matrix[i, 9] >= Q & sum_matrix[i, 10] <= 2 & sum_matrix[i, 11] == 0,
sum_matrix[i, 3] = 1;
sum_matrix[i, 12] = 1;
,
sum_matrix[i, 3] = 0;
sum_matrix[i, 12] = 0;
)
);
);
// return the value
Mean( sum_matrix[0, 12] );
)},
simMat[a, b] = simDissoTest( repN, 80 + cs[1] + cs[3] * (b - 1), a * rs[3] )
);
stop = Tick Seconds();
time = stopWatch( start, stop, "min" );
Show( time, ntests );
For Each () Version
// JSL Code for Simulation of Harmonized Compendial Dissolution Criteria USP <711>
// Simulation function for USP dissolution - three stage test
// reps = number of simulated lots for simulation
// mean = the process mean
// sd = the process standard deviation
// Q = value as specified in the mongraph for the product
//
// Summary Matrix Key
// ==================
// 1 "PassS1"
// 2 "PassS2"
// 3 "PassS3"
// 4 "#<Q S1"
// 5 "<Q-25 S1"
// 6 "Avg12"
// 7 "#<Q-15 S2"
// 8 "#<Q-25 S2"
// 9 "Avg24"
// 10 "#<Q-15 S3"
// 11 "#<Q-25 S3"
// 12 "PassAll"
// 13 "GoS2"
// 14 "GoS3"
// helper functions for the simulation
Clear Log();
vecSize = Function( {start, end, step},
{default local},
m = start :: end :: step;
steps = Length( m );
steps;
);
stopWatch = Function( {stTick, endTick, units = "sec"},
{default local},
seconds = endTick - stTick;
result = Match( units,
"sec", seconds,
"min", seconds / 60,
"hr", seconds / 3600,
"day", seconds / 86400
);
result = Round( result, 3 );
);
disclaimer = Expr(
New Window( "License Agreement",
Modal,
Outline Box( "License Agreement for Corrective Code or Additional Functionality",
H List Box(
Spacer Box( Size( 50, 50 ) ),
V Scroll Box(
Size( 500 ),
Flexible( 1 ),
Text Box(
"JMP STATISTICAL DISCOVERY LLC (“JMP”) PERMITS THE USE OF THIS COMPUTER SOFTWARE CODE (“CODE”) ON AN AS-IS BASIS AND AUTHORIZES YOU TO USE THE CODE SUBJECT TO THE TERMS LISTED HEREIN. BY USING THE CODE, YOU AGREE TO THESE TERMS. YOUR USE OF THE CODE IS AT YOUR OWN RISK. JMP MAKES NO REPRESENTATION OR WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRIGEMENT, AND TITLE, WITH RESPECT TO THE CODE.
You may use the Code solely as part of a software product you currently have licensed from JMP, JMP’s parent company SAS Institute Inc. (“SAS”), or one of JMP’s or SAS’s subsidiaries or authorized agents (the “Software”), and not for any other purpose. The Code is designed to add functionality to the Software but has not necessarily been tested. Accordingly, JMP makes no representation or warranty that the Code will operate error-free. JMP is under no obligation to maintain, support, or continue to distribute the Code.
Neither JMP nor its licensors shall be liable to you or any third-party for any general, special, direct, indirect, consequential, incidental, or other damages whatsoever arising out of or related to your use or inability to use the Code, even if JMP has been advised of the possibility of such damages. Except as otherwise provided above, the Code is governed by the same agreement that governs the Software. If you do not have an existing agreement with JMP or SAS governing the Software, you may not use the Code.
JMP and all other JMP Statistical Discovery LLC product or service names are registered trademarks or trademarks of JMP Statistical Discovery LLC in the USA and other countries. ® indicates USA registration. Other brand and product names are registered trademarks or trademarks of their respective companies.
",
<<Set Wrap( 500 )
)
)
)
)
)
);
// main simulation function
simDissoTest = Function( {reps = 100, mean = 100, sd = 5, Q = 80},
{default local},
// generate data following user chosen mean and sd
// generate 24 data values even if not needed (that is
// not going to all three stages)
data_matrix = J( nrow = reps, ncol = 24, Random Normal( mean, sd ) );
// Summary Matrix to hold simulation results
sum_matrix = J( reps, 14, Empty() );
// Populate the summary field
For Each( {e, i}, 1 :: reps,
v_6 = data_matrix[i, 1 :: 6];
v_12 = data_matrix[i, 1 :: 12];
v_24 = data_matrix[i, 1 :: 24];
sum_matrix[i, 4] = Sum( v_6[Loc( v_6 < (Q + 5) )] );
sum_matrix[i, 5] = Sum( v_6[Loc( v_6 < (Q - 25) )] );
sum_matrix[i, 6] = Mean( v_12 );
sum_matrix[i, 7] = Sum( v_12[Loc( v_12 < (Q - 15) )] );
sum_matrix[i, 8] = Sum( v_12[Loc( v_12 < (Q - 25) )] );
sum_matrix[i, 9] = Mean( V_24 );
sum_matrix[i, 10] = Sum( v_24[Loc( v_24 < (Q - 15) )] );
sum_matrix[i, 11] = Sum( v_24[Loc( v_24 < (Q - 25) )] );
// Stage 1
If( sum_matrix[i, 4] == 0,
sum_matrix[i, 1] = 1;
sum_matrix[i, 12] = 1;
sum_matrix[i, 13] = 0;
sum_matrix[i, 14] = 0;
,
sum_matrix[i, 1] = 0;
// Check to see if any results are below Q-25 Stage 1
// If there are any - do not proceed - fail test
If( sum_matrix[i, 5] == 0,
sum_matrix[i, 13] = 1,
sum_matrix[i, 12] = 0;
sum_matrix[i, 13] = 0;
sum_matrix[i, 14] = 0;
);
);
// Stage 2
If( sum_matrix[i, 13] == 1,
// Check to see if any results are below Q-25 in Stage 2
// and there are more than 2 units that are less than Q-15
// If there are any - do not proceed - fail test
If( sum_matrix[i, 6] >= Q & sum_matrix[i, 7] == 0,
sum_matrix[i, 2] = 1;
sum_matrix[i, 12] = 1;
sum_matrix[i, 14] = 0;
,
sum_matrix[i, 2] = 0;
If( sum_matrix[i, 7] <= 2 & sum_matrix[i, 8] == 0,
sum_matrix[i, 14] = 1,
sum_matrix[i, 14] = 0;
sum_matrix[i, 12] = 0;
);
)
);
// Stage 3
If( sum_matrix[i, 14] == 1,
If( sum_matrix[i, 9] >= Q & sum_matrix[i, 10] <= 2 & sum_matrix[i, 11] == 0,
sum_matrix[i, 3] = 1;
sum_matrix[i, 12] = 1;
,
sum_matrix[i, 3] = 0;
sum_matrix[i, 12] = 0;
)
);
);
// return the value
Mean( sum_matrix[0, 12] );
);
//---------------------------------------------------------------------------------
// Parallelized Simulation
//---------------------------------------------------------------------------------
// define the size of the matrix ans the step size
rows = vecSize( 0.1, 12, 0.1 );
rowShape = [0.1, 12, 0.1];
cols = vecSize( -1, 10, 0.5 );
colShape = [-1, 10, 0.5];
reps = 100;
ntests = rows * cols;
//create the Parallel test matrix
simMat = J( rows, cols, Empty() );
// run the simulation (testing scalability)
//New Window( "test", lb = Lineup Box( N Col( 2 ) ) );
/*For Each( {repN, i}, {1, 50, 100, 500},
start = Tick Seconds();
For Each( {v, {i, j}}, simMat,
simMat[i, j] = simDissoTest( repN, 80 + colShape[1] + colShape[3] * (j - 1), i * rowShape[3] )
);
stop = Tick Seconds();
time = stopWatch( start, stop, "min" );
Show( repN, time, ntests );
//lb << Append( V List Box( Text Box( "Number of reps: " || Char( repN ) ), New Heat Image( simMat ) ) );
);*/
// run the simulation
start = Tick Seconds();
For Each( {v, {i, j}}, simMat,
simMat[i, j] = simDissoTest( reps, 80 + colShape[1] + colShape[3] * (j - 1), i * rowShape[3] )
);
stop = Tick Seconds();
time = stopWatch( start, stop, "min" );
Show( reps, time, ntests );
// make the data table from the matrix
tallMatrix = J( ntests, 3, Empty() );
k = 1; // an indexer
For Each( {v, {i, j}}, simMat,
tallMatrix[k, 1] = 80 + colShape[1] + colShape[3] * (j - 1);
tallMatrix[k, 2] = i * rowShape[3];
tallMatrix[k, 3] = v;
k++;
);
dt = As Table( tallMatrix );
Column( dt, 1 ) << Set Name( "Mean" );
Column( dt, 2 ) << Set Name( "Std Dev" );
Column( dt, 3 ) << Set Name( "Pass Rate" );
dt << New Script(
"Graph",
Graph Builder(
Size( 742, 817 ),
Show Control Panel( 0 ),
Fit to Window,
Variables( X( :Mean ), Y( :Std Dev ), Color( :Pass Rate ) ),
Elements( Contour( X, Y, Legend( 6 ), Number of Levels( 5 ), Smoothness( 0.146217717445672 ) ) ),
SendToReport(
Dispatch(
{},
"Mean",
ScaleBox,
{Min( 79.0036098075223 ), Max( 89.9788461538462 ), Inc( 2 ), Minor Ticks( 1 )}
),
Dispatch(
{},
"Std Dev",
ScaleBox,
{Min( 0.0888646043293861 ), Max( 11.971158234127 ), Inc( 2 ), Minor Ticks( 1 )}
),
Dispatch(
{},
"400",
ScaleBox,
{Legend Model(
6,
Properties(
0,
{gradient(
{Color Theme( "Green Yellow Red"(1) ), Scale Values( [0.8 .] ),
Reverse Gradient( 1 )}
)},
Item ID( "Pass Rate", 1 )
)
)}
)
)
)
);
dt << Run Script("Graph");
disclaimer;
