Orthogonal table 实现理论

Orthogonal table

 

概述
正交试验中使用的正交表的生成过程其实就是一个加法Abel群的生成过程。对于m水平表而言,为了方便起见我们对水平表的每个元素均做的处理,任取其中的任意两列则一定具有完全对:

这说明a和b都是一个分量包含0~m-1的列向量,并且有:

也就是说在 上 或者 也都具有完全对。

至此我们可以定义加法:

对于标准表而言,若a和b是标准表的 2 个列,那么也一定是标准表的一个列。

有限Abel群的基底
考虑到标准表中任取 2 列的完全对出现的次数都是均等的 , 因此的列总是以

中的元素为基础关于运算生成的，其中 是矩阵的Kronecker乘积。

有限Abel群的子群

令G是以为基底，关于运算生成的Abel群，则G明显是拥有众多Abel子群的，例如在确定G的基底是

的情况下，显然：

是G的子群。现在我们扩展一下这个系列。令：

此时都是G的m阶循环子群。

构造

Step1.K的全体元素就是正交表的所有列。正交试验标准表可按此法构造。

Step2.若有 ,则说明第i列和第j列对应放置交互作用的列应该是第列。正交试验表头可按照此法构造。

后记

此法没有涉及任何的枚举过程，计算出的正交表和表头表的生成速度比正常枚举法要快非常多，正常的单线程运行此算法，一个行数和列数均在200以上的正交表只需要1秒不到就能算出。

值得注意的是Step2中所描述的正交表表头并非是一定能以此法完整的构造出的，因为Step2中的等号并非一定成立。为什么会出现这种情况呢，原因就是因为的构造方式并不能保证对于任意的i和j 都能找到使得Step2中等号都的成立。按照道理来说我们需要找一个只由构造出的循环群，如果这样的都能满足Step2中的等式，则正交表头可完整构造出来。

Python Code(import numpy)

class OrthogonalTableMap():
    
    def __init__(self, level, time):
        self.m = level
        self.time = time
        n = self.m**(time + 1)
        k = (n - 1)//(self.m - 1)
        
        # 构造self.base，self.base的元素与加法群的基底一一对应
        self.base = [numpy.zeros(n, numpy.int)]
        for j in range(0, time + 1):
            mat1 = numpy.ones(self.m**j, numpy.int)
            mat2 = numpy.linspace(1, self.m, self.m, dtype = numpy.int)
            mat3 = numpy.ones(self.m**(time - j), numpy.int)
            self.base.append(numpy.kron(numpy.kron(mat1, mat2), mat3) - 1)
        self.base = self.base[1:len(self.base)+1]
        
        # 构造正交表self.table和记录每列群属性的self.dir
        self.dir = [numpy.zeros((time + 1, self.m - 1), numpy.int)]
        self.dir[0][0, 0] = 1
        for j in range(1, self.m - 1):
            self.dir[0][:, j] = (j + 1)*self.dir[0][:, 0]%self.m
        
        self.table = [self.base[0]]
        for j in range(1, time + 1):
            self.table.append(self.base[j])
            
            self.dir.append(numpy.zeros((time + 1, self.m - 1), numpy.int))
            self.dir[-1][j, 0] = 1
            for k in range(1, self.m - 1):
                self.dir[-1][:, k] = (k + 1)*self.dir[-1][:, 0]%self.m
            
            for k in range(0, len(self.table) - 1):
                for s in range(1, self.m):
                    self.table.append((s*self.table[k] + self.base[j])%self.m)
                    
                    self.dir.append(numpy.zeros((time + 1, self.m - 1), numpy.int))
                    self.dir[-1][:, 0] = s*self.dir[k][:,0]
                    self.dir[-1][j, 0] += 1
                    self.dir[-1][:, 0] = self.dir[-1][:, 0]%self.m
                    for t in range(1, self.m - 1):
                        self.dir[-1][:, t] = (t + 1)*self.dir[-1][:, 0]%self.m
                        
        # 将self.base和self.table以矩阵形式表示，并且+1还原
        temp = numpy.zeros((n, len(self.base)), numpy.int)
        for j in range(0, len(self.base)):
            temp[:, j] = self.base[j]
        self.base = temp + 1
        temp = numpy.zeros((n, len(self.table)), numpy.int)
        for j in range(0, len(self.table)):
            temp[:, j] = self.table[j]
        self.table = temp + 1
        
        # self.hash_dir的一个元素都是self.table的一个列对应的群成员的hash值集合
        self.hash_dir = self.hash_function()
        
        self.error = 0
        
    def hash_function(self):
        # 制作hash函数将self.dir所包含的信息处理成数据类型为numpy.int的基本元素的集合簇输出
        mat = numpy.ones((1, self.time + 1), numpy.int)
        for j in range(1, self.time + 1):
            mat[0, j] *= self.m**j
        hash_dir = self.dir.copy()
        for j in range(0, len(self.dir)):
            hash_dir[j] = set(numpy.dot(mat, self.dir[j]).reshape(self.m - 1))
            
        return hash_dir
    
    def hash_add(self, var_set, format_check = True):
        # 制作与向量加法对应在hash映射下的加法hash_add
        if format_check:
            if len(var_set)<2:
                print("Error:var_set should at least has 2 elements.")
            if max(var_set)>len(self.hash_dir) or min(var_set)<0:
                print("Error:var_set's elements should between 0~len(hash_dir).")
        arg_set = numpy.array(var_set)
        
        var_out = 0
        for j in range(0, self.time + 1):
            temp = 0
            for k in range(0, len(var_set)):
                temp += arg_set[k]//self.m**(self.time - j)
            var_out += (temp%self.m)*self.m**(self.time - j)
            arg_set = arg_set%self.m**(self.time - j)
            
        return var_out
    
    def hash_set_add(self, sets_ind, end_set = None):
        # 定义hash后的self.dir的集合对应在hash映射下的加法hash_set_add
        if len(sets_ind)==0:
            set_num = list()
            for j in range(0, len(self.hash_dir)):
                if len(self.hash_dir[j]&end_set)>0:
                    set_num.append(j)
            return set_num
        
        if end_set is None:
            if len(sets_ind)<2:
                print("Error:set_ind should at least has 2 elements.")
            end_set = self.hash_dir[sets_ind[-1]]
            sets_ind.pop()
        
        temp_set = list()
        for j in self.hash_dir[sets_ind[-1]]:
            for k in end_set:
                temp_set.append(self.hash_add([j, k], format_check = False))
        end_set = set(temp_set)
        sets_ind.pop()
        
        return self.hash_set_add(sets_ind, end_set)
    
    def title(self):
        # 使用矩阵title表示带有多因素交互作用的表头，若没有成功找到所有子群则修改self.error = 1
        title = numpy.zeros((len(self.hash_dir)*(self.m - 1), 
                                  len(self.hash_dir)), numpy.int)
        try:
            for j in range(0, len(self.hash_dir)):
                title[(self.m - 1)*j:(self.m - 1)*(j + 1), j] = j + 1
                for k in range(j + 1, len(self.hash_dir)):
                    temp = numpy.array(self.hash_set_add([j, k])).reshape(self.m -1)
                    title[(self.m - 1)*j:(self.m - 1)*(j + 1), k] = temp +1
            self.error = 0
        except ValueError:
            self.error = 1
        return title

Demo

3水平13因素正交表(表头数据依赖于列排序，列排序不一定和文献资料一样，如需对照请自行重新排序对照)

ex = OrthogonalTableMap(3,2)

print("table = ")
print(ex.table)
print("title = ")
print(ex.title())

声明

本文由“Hamilton算符”原创，未经博主允许不得转载！