title:- A=['% design0.pl : experimental design using extended Galois field' ,'% By Kenryo Indo (Kanto Gakuen University)' ,'% file: design0.pl ( for SWI-Prolog 5.0.9)' ,'% created: 20-21 Aug 2003' ,'% modified: 15 Jul 2004 (corrected the positions of dynamic declarations)' ,'% Reference:' ,'% Lecture Note of Simulation (by Prof. Iwaro Takahashi, University of Tsukuba, 1987). ' ,'% import: matrix01.pl(20 Aug 2003) matrix operations.' ,'%-- the main has prolonged in the last part of the program ---------' ,'% othogonal_table/2.' ,'% start_design/0.' ] ,forall(member(B,A),(nl,write(B))). :- [matrix01]. % an expermental design problem %--------------------------------- :- dynamic number_of_categories/1. :- dynamic number_of_factors/1. :- dynamic number_of_trials/1. number_of_categories(2). number_of_factors(7). number_of_trials(8). % experimental design using GF(S). %--------------------------------- % careate theta(TH,K), K=1, ..., L, a minimal % L-dimensional independent vectors on GF(2) % which covers all factors above. design_pattern(S,S):-member(S,[2,3]). design_pattern(4,2^2). all_auxiliary_matrix(using_GF(S),dimension(L),aux(AQ)):- design_pattern(S,S1), findall(Q, auxiliary_vector(using_GF(S1),dimension(L),aux(Q)), AQ). minimal_auxiliary_matrix(design(S,M,N),dimension(L),aux(AQ),min(AP)):- number_of_categories(S), minimal_dimension(L,design(S,M,N)), !, all_auxiliary_matrix(using_GF(S),dimension(L),aux(AQ)), findall([(Max,Sum)|P], ( member(P,AQ), has_a_non_zero_member(P,Sum), max(Max,member(Max,P)) ), AP0), sort(AP0,AP1), findall_up_to(M,P,member([_|P],AP1),AP). findall_up_to(N,Target,Goal,Bag):- findall(Target,Goal,Bag0), findall(Target,(kth(K,Bag0,Target),K== Dim. covering_dimension_1(using_GF(S),factors(M),trials(N),dimension(L)):- length(F,M), nth1(L,F,L), findall(P, auxiliary_vector(using_GF(S),dimension(L),aux(P)), Ps1), length(Ps1,D), M =< D, N >= D. auxiliary_vector(using_GF(S),dimension(L),aux(P)):- member(S,[2,3]), number(L), anum_seq(A,S), bag0(P,A,L). auxiliary_vector(using_GF(2^2),dimension(L),aux(P)):- number(L), bag0(P,[0,1,b,1+b],L). aux(S,L,P):- auxiliary_vector(using_GF(S),dimension(L),aux(P)). %-------------------------------------- % create othogonal table on GF(S) %-------------------------------------- :- dynamic table1 / 2. othogonal_table(design(S,M,N),AQ*AP=W):- write('abolish tables ? >'), (read(y)->abolish(table1/2);true), minimal_auxiliary_matrix(design(S,M,N),dimension(_L),aux(AQ),min(AP)), matrix_transpose(AP -> APT), design_pattern(S,S1), matrix_multiply_GF(S1,AQ * APT = Z), matrix_mod(W is (Z mod S1)), assert(table1(design(S,M,N),AQ*AP=W)), write('display tables ? >'), (read(y)->pp_table(S1);true). pp_table(2^2):- table1(design(S,M,N),AQ*AP=W), nl,write(table1(design(S,M,N))), matrix_transpose(AQ -> AQT), matrix_transpose(AP -> APT), MM is M * 3, nl,write_n_times('-',MM+3), write_table_GF(AQT), nl,write_n_times('-',MM), write_table_GF(APT), nl,write_n_times('-',MM), %forall(member(X,W),(nl,write(X))), write_table_GF(W), matrix_transpose(W -> WT), sum_table_GF(WT,Sum), nl,write_n_times('-',MM), nl,write(Sum). pp_table(S):- member(S,[2,3]), table1(design(S,M,N),AQ*AP=W), nl,write(table1(design(S,M,N))), matrix_transpose(AQ -> AQT), matrix_transpose(AP -> APT), MM is M * 3, nl,write_n_times('-',MM+3), write_table(AQT), nl,write_n_times('-',MM), write_table(APT), nl,write_n_times('-',MM), write_table(W), matrix_transpose(W -> WT), sum_table(WT,Sum), nl,write_n_times('-',MM), nl,write(Sum). sum_table(W,Sum):- findall(S,(member(Q,W), sum(Q,S)),Sum). sum_table_GF(W,Sum):- findall(S,(member(Q,W), sum_GF(Q,S)),Sum). write_table(W):- forall(member(Q,W), (sum(Q,S),nl,write(Q),tab(2),write((S)))). write_table_GF(W):- forall(member(Q,W), (sum_GF(Q,S),nl,write(Q),tab(2),write((S)))). write_n_times(Q,N):- N1 is N, len(A,N1), forall(member(Q,A), (write(Q))). %----------------------- % utilities %----------------------- has_a_non_zero_member(P,Sum):- \+ ( member(AA,P), \+ number(AA) ), sum(P,Sum), Sum > 0, !. has_a_non_zero_member(P,Sum):- ( member(AA,P), \+ number(AA) ), sum_GF(P,Sum). sum_GF([],0). sum_GF([A|P],Sum):- sum_GF(P,S0), is_GF(Sum, A+S0). min(X,Goal):- max(-X,Goal). max(X,Goal):- % X: the objective variable, % Goal: the objective function and constraints, setof((X,Goal),Goal,Z), forall(member(A,Z),(nl,write(A))), max_0(X,Z). max_0(X,Z):- \+ ( member((AA,_),Z), \+ number(AA) ), member((X,_Goal),Z), %nl,write((case1,X,_Goal)), \+ ( member((Y,_),Z), Y > X ). % modified: 22 Aug 2003 a cut !. max_0(X,Z):- member((AA,Goal),Z), \+ number(AA), !, %nl,write((case2,AA,Goal)), sort(Z,Z1), rev(Z1,[(X,_)|_]). %--------------------------------- % operations on extended GF %--------------------------------- operation_GF(E):- member(X,[0,1,b,1+b]), member(Y,[0,1,b,1+b]), member(F,[+,-,*]), E=..[F,X,Y]. :- dynamic is_GF/2. operations(on_GF(2^2),RULE):- % primitive reduced polynominal: 1+x+x^2=0 RULE=[ % (is_GF(A,B):- A is B), is_GF(0, 0+0), is_GF(0, 0-0), is_GF(0, 0*0), is_GF(1, 0+1), is_GF(1, 0-1), is_GF(0, 0*1), is_GF(b, 0+b), is_GF(b, 0-b), is_GF(0, 0*b), is_GF(1+b, 0+ (1+b)), is_GF(1+b, 0- (1+b)), is_GF(0, 0* (1+b)), is_GF(1, 1+0), is_GF(1, 1-0), is_GF(0, 1*0), is_GF(0, 1+1), is_GF(0, 1-1), is_GF(1, 1*1), is_GF(1+b, 1+b), is_GF(1+b, 1-b), is_GF(b, 1*b), is_GF(b, 1+ (1+b)), is_GF(b, 1- (1+b)), is_GF(1+b, 1* (1+b)), is_GF(b, b+0), is_GF(b, b-0), is_GF(0, b*0), is_GF(1+b, b+1), is_GF(1+b, b-1), % X+1=b-1+1->X=(1+b) is_GF(b, b*1), is_GF(0, b+b), is_GF(0, b-b), is_GF(1+b, b*b), % since p.r.p. b^2 = -(1+b)=1+b is_GF(1, b+ (1+b)), is_GF(1, b- (1+b)), is_GF(1, b* (1+b)), % b*(1+b)=b+b^2=b+(1+ b)=1+(b+b)=1 is_GF(1+b, 1+b+0), is_GF(1+b, 1+b-0), is_GF(0, (1+b)*0), is_GF(b, 1+b+1), is_GF(b, 1+b-1), is_GF(1+b, (1+b)*1), is_GF(1, 1+b+b), is_GF(1, 1+b-b), is_GF(1, (1+b)*b), % (1+b)*b=b+b^2=b+(1+b)=1 is_GF(0, 1+b+ (1+b)), is_GF(0, 1+b- (1+b)), is_GF(b, (1+b)* (1+b)) %(1+b)* (1+b)=1+2b+b^2=1+(1+b)=b ]. assert_GF(2^2):- operations(on_GF(2^2),RULE), forall(member(X,RULE),assertz((X:-!))). :- assert_GF(2^2). % matrix multiplication for GF(2^2) % ----------------------------------------------------------- % matrix_multiply_GF(2^2,X * Y = Z):- !, \+ number(X), (matrix(X,A);X=A), (matrix(Y,B);Y=B), matrix_of_size(A,N,L0), matrix_of_size(B,L0,M), matrix_of_size(Z,N,M), matrix_transpose(Y -> C),!, findall(ZK, ( kth(K,Z,ZK), findall(ZKL, ( kth(L,ZK,ZKL), matrix_rows(A,[K],[AK]), matrix_rows(C,[L],[CL]), eq_sum_product_GF(AK,CL,_,ZKL,_) ), ZK) ), Z). matrix_multiply_GF(_S,X * Y = Z):- matrix_multiply(X * Y = Z). eq_sum_product_GF([],[],[],0,0). eq_sum_product_GF([P|Q],[A|B],[E|F],V,Vq):- len(Q,N), len(B,N), eq_sum_product_GF(Q,B,F,V1,Vq1), is_GF(E, P * A), is_GF(V, V1 + E), (Vq1=0 -> Vq = EP*EA; Vq = Vq1 + (EP * EA)). % problems %-------------------------------------- problem(a, number_of_categories(2), number_of_factors(7), number_of_trials(8) ). problem(b, number_of_categories(3), number_of_factors(12), number_of_trials(27) ). problem(c, number_of_categories(2), number_of_factors(10), number_of_trials(16) ). problem(d, number_of_categories(4), number_of_factors(5), number_of_trials(16) ). % setting the problem %--------------------------------- :- dynamic current_problem/1. select_problem(A):- problem(A,B,C,D), initialize_current_problem, update_current_problem(A,[B,C,D]). initialize_current_problem:- abolish(current_problem/1), abolish(number_of_categories/1), abolish(number_of_factors/1), abolish(number_of_trials/1). update_current_problem(A,Target):- assert(current_problem(A)), forall(member(Y,Target),assert(Y)). start_design:- CREATE_TABLE= othogonal_table(design(S,M,N),_), nl, write_n_times('*',15), write(' welcome '), write_n_times('*',15), nl, write(' this is prolog program of the experimental design problems. '), nl, write(' available samples as follows:'), P=problem(A, number_of_categories(S), number_of_factors(M), number_of_trials(N) ), forall(P,(nl,tab(2),write((A,categories:S,factors:M,trials:N)))), nl, nl, write(' please select a problem above (default is problem a):'), read(U), (select_problem(U)->true; (nl,write('no such problem.'),start_design0) ), nl, write(' ok. '), nl, write(' start to create table? :'), read(y), CREATE_TABLE, nl, write(' the table has created (refer table1/2).'), nl, write(' continue design ? :'), (read(y)->start_design;!,fail). start_design. % auto start. :- start_design. %------------------------------------------------