Names Default to Here(1); // function to make a design matrix row. designRow = Function( {a,b}, {Default Local}, Matrix( {1, a, b, a*b, a^2, b^2} ) ); calcAverageVariance = Function( {n,x1,x2,Moments}, {Default Local}, X = J( n, 1, 1 ) || x1 || x2 || x1 :* x2 || x1^2 || x2^2; av = trace(inverse(X`*X)*Moments); ); calcDeterminant = Function( {n,x1,x2}, {Default Local}, X = J( n, 1, 1 ) || x1 || x2 || x1 :* x2 || x1^2 || x2^2; d = det(X`*X); ); calcInverseInfo = Function( {x1,x2}, {Default Local}, n = N Row(x1); X = J( n, 1, 1 ) || x1 || x2 || x1 :* x2 || x1^2 || x2^2; xpxi = inverse(X`*X); ); makeScatterPlot = Function( {x1,x2,oldx1,oldx2,detV,iVar}, {Default Local}, g1 = Graph Box( Frame Size( 300, 300 ), X Name( "X1" ), X Scale( -1.1, 1.1 ), Y Name( "X2" ), Y Scale( -1.1, 1.1 ), Double Buffer, Marker Size( 2 ); Pen Color( "Red" ); For( l=1, l<=n, l++, Arrow( {oldx1[l], oldx2[l]}, {x1[l], x2[l]} ) ); Marker( x1, x2 ); Text( Center Justified, Erased, {0, 0.75}, "Determinant = " || Char( detV, 6, 99 ) ); Text( Center Justified, Erased, {0, -0.75}, "Average Variance = " || Char( iVar, 6, 99 ) ); ); // graph ); makeContourPlot = Function({x1,x2}, g2 = Graph Box( Frame Size( 300, 300 ), X Name( "X1" ), X Scale( -1.1, 1.1 ), Y Name( "X2" ), Y Scale( -1.1, 1.1 ), Double Buffer, Contour Function( xh = designRow( tx, ty ); xh` * xpxi * xh , tx, ty, ((2::12)-1) / 10, Z Color( (1::11) + 2 ) ); Marker( x1, x2 ); ) ); dlg = Dialog( Line Up( 2, "Number of runs", n = Edit Number( 9 ), "Pause (sec)", w = Edit Number( 0.05 ) ), V List( "Optimality criterion", type = Radio Buttons( "D-optimal", "I-optimal" ) ), pause = Check Box( "Pause at each step", 1 ), H List( Button( "OK" ), Button( "Cancel" ) ) ); If( dlg["Button"] == -1, Throw( "User cancelled" ) ); Remove From( dlg ); Eval List( dlg ); n = Floor( n ); If( n < 6, Throw( Dialog( "Insufficient number of runs (6 minimum)", "", Button( "OK" ) ) ) ); //Start with uniformly random points. x1 = J( n, 1, round(2 * Random Uniform() - 1,2)); x2 = J( n, 1, round(2 * Random Uniform() - 1,2)); // display design in data table. dt = New Table( "DOE from Script" ); c1 = dt << New Column( "X1", Values( x1 ) ); c2 = dt << New Column( "X2", Values( x2 ) ); dt << Size Window( 500, 400 ); dt << Move Window( 50, 500 ); moments = [1 0 0 0 0.333 0.333, 0 0.333 0 0 0 0, 0 0 0.333 0 0 0, 0 0 0 0.111 0 0, 0.333 0 0 0 0.2 0.111, 0.333 0 0 0 0.111 0.2]; oldx1 = x1; oldx2 = x2; c = J( 6, 6, 1 ); iVar = calcAverageVariance(n,x1,x2,moments); avgVar = iVar; detV = calcDeterminant(n,x1,x2); xpxi = calcInverseInfo(x1,x2); nw = New Window( "Coordinate Exchange - Two Factor RSM", H List Box( Outline Box( "Coordinate Exchange Demo", g1 = makeScatterPlot(x1,x2,oldx1,oldx2,detV,iVar), Text Box( Choose( type, "D-optimal", "I-optimal" ) || " design for " || Char( n ) || " runs" ) ), // outline Outline Box( "Relative Prediction Variance", g2 = makeContourPlot(x1,x2), ) // outline ) // h list ); If( pause, Wait(w) ); nw << Move Window( 50, 25 ); dMax = detV; minVar = iVar; maxIndex = 0; done = 0; While( Not( done ), If( pause, Wait( w ) ); done = 1; For( i=1, i<=n, i++, For( j=1, j<3, j++, u = x1; v = x2; If ( type == 1, For( k=-1, k<2, k++, r = 1; If( j == 1, u[i] = k, v[i] = k ); detM = calcDeterminant(n,u,v); If( detM > dMax, dMax = detM; maxIndex = k; ); ), For( k=-10, k<11, k++, r = 1; If( j == 1, u[i] = k/10.0, v[i] = k/10.0 ); avgVar = calcAverageVariance(n,u,v,moments); If ( avgVar < minVar, minVar = avgVar; maxIndex = k / 10.0; ) ); ); If( ( (type == 1) & (dMax > detV) ) | ( (type == 2) & (minVar < iVar) ), done = 0; If( j == 1, u[i] = maxIndex, v[i] = maxIndex ); detV = calcDeterminant(n,u,v); iVar = calcAverageVariance(n,u,v,moments); oldx1[i] = x1[i]; oldx2[i] = x2[i]; dx = u[i] - oldx1[i]; g1 << Reshow; If( Abs( dx ) < 0.0001, dy = v[i] - oldx2[i]; For( r=1, r<6, r++, ychange = r * dy / 5.0; x2[i] = oldx2[i] + ychange; g1 << Reshow; If( pause, Wait( w ) ); ), For( r=1, r<6, r++, xchange = r * dx / 5.0; x1[i] = oldx1[i] + xchange; g1 << Reshow; If( pause, Wait( w ) ); ) ); x1[i] = u[i]; x2[i] = v[i]; oldx1[i] = x1[i]; oldx2[i] = x2[i]; c1 << Values( x1 ); c2 << Values( x2 ); If( j==1, Row State(i) = Color State(3), Row State(i) = Color State(1) ); If( pause, Wait( w ) ); xpxi = calcInverseInfo(x1,x2); {g1,g2} << Reshow; ); // end if (Dmax > detV) ); // end loop on j (columns) ); // end loop on i (rows) ); // end while