/* Minimum size mixed level orthogonal arrays */ /* To replicate examples of table 1 define your problem providing the levels of all the factors and the required strength*/ /* example n.1: 2^6 16 */ %let lev={2 2 2 2 2 2 16}; /* six 2-level factors and one 16-level factors */ %let t=2; /*strength*/ /* example n.2: 4^3 8 */ *%let lev={4 4 4 8}; /* three 4-level factor, one 8-level factor*/ *%let t=2; /*strength*/ /* example n.3: 3^4 12 */ *%let lev={3 3 3 3 12}; /* four 3-level factor and one 12-level factor*/ *%let t=2; /*strength*/ /* example n.4: 4^3 12 */ *%let lev={4 4 4 12}; /* three 4-level factor and one 12-level factor*/ *%let t=2; /*strength*/ /* example n.5: 3 6 9 */ *%let lev={3 6 9}; /* one 3-level factor, one 6-level factor and 9-level factor*/ *%let t=2; /*strength*/ /* example n.6: 4^3 16 */ *%let lev={4 4 4 16}; /* three 4-level factor and one 16-level factor*/ *%let t=2; /*strength*/ /* example n.7: 3^4 9 */ *%let lev={3 3 3 3 9}; /* four 3-level factor and one 9-level factor*/ *%let t=2; /*strength*/ /* CONSTRAINTS in LPSOLVE format */ /*Note: in some examples for the optimization phase is preferable to use lpsolve because it is better from the point of view of efficiency*/ /* If you want to build an external file that contains the constraints in lpsolve format go to the line labeled by YYY and write the folder in which your file should be written */ /* To replicate the examples of Minimum size orthogonal fractional factorial designs OFFD reported in section 6.1: */ /* First define your problem providing the levels of the factors*/ *%let lev={3 3 3 3}; /* 4 factors, 3 levels */ /* then replace the line labeled by XXX with the following */ *if norm<=2 | (alpha[4]=0 & alpha[1:3]^= 0) | (alpha[1]=0 & alpha[2:4]^=0) then do; /*main effects D1,D2,D3,D4 and the interaction D2*D3 estimable and orthogonal*/ /************************** CONSTRUCTION OF THE CONSTRAINTS ***********************/ proc iml; /********* Full factorial design with n_1... n_m levels *******/ start ffmix(lev); nfac=ncol(lev); /*number of factors*/ npt=1; do i=1 to nfac; npt=npt*lev[i]; /* cardinality of the design*/ end; ff=j(npt,nfac,0); do i=0 to npt-1; tmp=i; do j=1 to nfac; ff[i+1,nfac-j+1]= mod(tmp,lev[nfac-j+1]); /*conversion base algorithm*/ tmp=int(tmp/lev[nfac-j+1]); end; end; return(ff); finish; /*********************** c_alpha=0 ***********************/ start c_null(L,lev,t) global(cs,smax); do i=2 to nrow(L); alpha=L[i,]; c=loc(alpha); norm=ncol(c); /*norm of alpha, number of non-null elements of alpha*/ if norm<=t then do; /*line XXX: to replicate the generic OFFD example replace this line with the following one:*/ *if norm<=2 | (alpha[4]=0 & alpha[1:3]^= 0) | (alpha[1]=0 & alpha[2:4]^=0) then do; col=j(nrow(L),1,0); /* Simple term: norm of alpha=1 */ if norm=1 then do; s = lev[c]/gcd(lev[c],alpha[c]); col = mod(alpha[c]*L[,c],lev[c])*(s/lev[c]); end; /* Interaction term: norm of alpha>1 */ else do; s=1; do j=1 to norm; cnum=c[j]; sj = lev[cnum]/gcd(lev[cnum],alpha[cnum]); s =lcm(s,sj); end; do j=1 to norm; cnum = c[j]; X = mod(alpha[cnum]*L[,cnum],lev[cnum])*(s/lev[cnum]); col = mod(col+X,s); end; end; b=b||col; cs=cs||s; end; end; smax=max(cs); return(b); finish; /************ Complex conjugate **************/ start cc(b,cs); bc=j(nrow(b),ncol(b),.); do i=1 to ncol(b); do k=1 to nrow(b); bc[k,i]=mod(cs[i]-b[k,i],cs[i]); end; end; return(bc); finish; /****************** CYCLOTOMIC POLYNOMIALS AND REMAINDERS ***************/ /* Polynomial division: dividend and divisor have numerical coefficients */ start f(d1,d2); n=ncol(d1); i=1; k=1; do while (d1[i]=0); /* Leading term of the dividend */ i=i+1; end; do while (d2[k]=0); /* Leading term of the divisor */ k=k+1; end; ltd2=d2[,k]; zero=j(1,n,0); q=zero; do while(i<=k); e=k-i; ltd1=d1[,i]; q[,n-e]=ltd1/ltd2; /* Quotient */ c=loc(d2); d2s=zero; /* Remainder */ do m=1 to ncol(c); p=c[m]; d2s[,p-e]=d2[,p]*q[,n-e]; end; d1=d1-d2s; u=ncol(unique(d1)); /* Verify if the remainder is null */ if u=1 then return (q); else do; i=1; /*Leading term of the dividend*/ do while (d1[i]=0); i=i+1; end; end; end; finish; /* Find all divisors of a number */ start dvs(n); do i=1 to n; if mod(n,i)=0 then d=d||i; end; return(d); finish; /* Polynomial division: the dividend P(z) has non-numerical coefficients */ start dp(s,c); /* s=degree of P, c=divisor polynomial*/ R=I(s); zero=j(s,1,0); i=1; k=1; do while (c[k]=0); /* Leading term of the divisor */ k=k+1; end; do while(i<=k); null=j(s,s,0); /* Leading term of the dividend */ ltR=null; ltR[,i]=R[,i]; e=k-i; cp=loc(c); /* Quotient */ cs=zero; do m=1 to ncol(cp); p=cp[m]; cs[p-e,]=c[p]; end; cc=loc(cs); /* Remainder */ qG=null; do m=1 to ncol(cc); p=cc[m]; qG[,p]=R[,i]*cs[p]; end; R=R-qG; zero=j(s,1,0); /* Leading term of the dividend */ i=1; do while (R[,i]=zero); i=i+1; end; end; return (R); finish; /* Remainders and Pointers */ start rmd(smax) global(pointers); do s=2 to smax; /* Calculate the cyclotomic polynomial phi_s by dividing X^s-1 by the cyclotomic polynomials of the proper divisors of s */ d=dvs(s); /* Divisors of s*/ nd=ncol(d); if nd>2 then do; /* If s is not prime */ c=j(nd,s+1,0); c[1,s+1]=-1; /* phi_1*/ c[1,s]=1; do k=2 to nd; dv=d[,k]; g=dvs(dv); num=c[k,]; num[,s+1]=-1; /* Dividend */ num[,s+1-dv]=1; if ncol(g)=2 then do; /* If it is prime */ q=f(num,c[1,]); /* Quotient */ c[k,]=q; end; else do; /* If it is not prime */ do h=1 to ncol(g)-1; gg= g[h]; l=loc(d=gg); q=f(num,c[l,]); /* Quotient */ num=q; end; c[k,]=q; end; end; y=c[nd,2:s+1]; phi_s=t(y); /* Cyclotomic polynomial phi_s*/ /*The remainder matrix*/ r=dp(s,phi_s); /* We exclude null coefficients*/ free remainderstmp; do k=1 to ncol(r); tmp=ncol(unique(r[,k])); if tmp^=1 | (tmp=1 & unique(r[,k])^=0) then remainderstmp=remainderstmp//t(r[,k]); end; end; else do; /* if s is prime */ I=I(s-1); ones=j(s-1,1,-1); free remainderstmp; remainderstmp=ones||I; end; /* dots in the boxes that are empty */ if smax-ncol(remainderstmp)>0 then do; miss=j(nrow(remainderstmp),smax-ncol(remainderstmp),.); remainderstmp=remainderstmp||miss; end; pointers=pointers//(s||nrow(remainders)+1||nrow(remainders)+nrow(remainderstmp)); remainders=remainders//remainderstmp; end; pointers=pointers[2:smax,]; return(remainders); finish; /*********************** CONSTRAINTS: AY=0 ***********************/ start cns(cs,bc,pointers,remainders); free bs; do i=1 to ncol(bc); u=unique(bc[,i]); s=cs[i]; free co; free r; co=pointers[s-1,2:3]; do j=co[1] to co[2]; r=r||t(remainders[j,1:s]); /* Remainder */ end; do j=1 to ncol(r); tr=j(nrow(bc),1,0); do k=1 to s; tr[loc(bc[,i]=u[k]),]=r[s-k+1,j]; /* Constraints*/ end; bs=bs||tr; end; end; A=t(bs); return(A); finish; /*********************** MAIN ***********************/ lev=&lev.; t=&t.; L=ffmix(lev); /* Full factorial design */ print L; create L from L; append from L; close L; cs=0; /* c_alpha=0 and vector of s */ smax=0; b=c_null(L,lev,t); cs=cs[,2:ncol(cs)]; bc=cc(b,cs); /* Complex conjugate */ free remainders; /* Remainders and pointers */ pointers=j(1,3,0); remainders=rmd(smax); print remainders; A=cns(cs,bc,pointers,remainders); /*Constraints*/ print A; /**************** to Optmodel ****************/ nc=1:nrow(A); nc=t(nc); create nrow from nc[colname="R"]; append from nc; close nrow; /*Variables names*/ y=j(1,ncol(A)," "); do i=1 to ncol(y); y[i]=compress(concat("x",char(i))); end; /*Constraints dataset*/ create cs from A [colname=y]; append from A; close cs; /* Complete dataset */ data cs_compl; merge nrow cs; run; quit; proc contents data=cs out=tmp noprint; run; data _NULL_; set tmp; call symput('nvar',trim(left(_n_))); run; %put "NVAR" &nvar.; /* Number of variables */ /*********************** to lpsolve ***********************/ %let length=%eval(8*&nvar.+5); %put &length; data firstrow; keep str; length str $ &length.; str="min: S;"; run; data cnts; keep str; length str $ &length.; set cs; array x {&nvar.} x1-x&nvar.; str=""; do i=1 to &nvar.; if x[i]^=0 then do; if x[i]<0 then a=trim(left(put(x[i],$8.))); else a='+'||trim(left(put(x[i],$8.))); ii=trim(left(put(i,$8.))); str=compress(str||a||"x"||ii); end; end; str=compress(str||"=0;"); output; run; data objt; keep str; length str $ &length.; array x {&nvar.} x1-x&nvar.; str=""; do i=1 to &nvar.; ii=trim(left(put(i,$8.))); str=compress(str||"+x"||ii); end; str=compress(str||"=S;"); output; run; data v; keep str; length str $ &length.; str="S>=1;"; run; data int; keep str; length str $ &length.; array x {&nvar.} x1-x&nvar.; str="int_"; do i=1 to &nvar.; ii=trim(left(put(i,$8.))); str=compress(str||"x"||ii||","); end; str=compress(str||"S;"); output; run; data lpsolve; set firstrow cnts objt v int; str=translate(str," ","_"); run; /* line YYY: Specify the location of the output file (&your_path) and make the proc export working by removing the comment signs */ *%let yourpath= ; /* e.g. %let yourpath=C:\mydoc */ /* PROC EXPORT DATA= WORK.Lpsolve OUTFILE= "yourpath\lpsolve_input.txt" DBMS=TAB REPLACE; PUTNAMES=NO; RUN; */ /*********************** Optimization phase ***********************/ proc optmodel; number n=&nvar.; set variables=1..n; set RowName; number C{RowName,Variables}; read data cs_compl into RowName=[R] {i in variables} < C[R,i]=col("x"||i) >; var x{variables}>=0 integer; min t= sum{i in variables} x[i]; con S:sum{i in variables} x[i]>=1; con Constraint{k in Rowname}: sum {i in variables} C[k,i]*x[i]=0; solve with milp / presolver = automatic heuristics = automatic inttol=1e-10; print x; create data solution from [i] solution=x; /* Fraction */ quit; /* Cardinality of the fraction */ proc means data=solution sum; var solution; run; proc iml; use L; read all into L; close L; use solution; read all into solution; close solution; epsilon=.001; tmp=L||solution[,2]; fraction=tmp[loc(solution[,2]>=1-epsilon),]; create fraction from fraction; append from fraction; close fraction; print "Number of different points in the fraction" (nrow(fraction)); print "Size of the fraction" (fraction[+,ncol(fraction)]);