MATLAB软件课程设计---MATLAB实现DNA序列的分类识别

装 订 线 装 订 线 数学软件实践 课程设计报告 题目：MATLAB 实现 DNA 序列的分类识别 学号： 姓名： 教师 日期： 论文题目 论文题目. 2 摘要. 3 1 问题重述. 4 2 问题分析. 4 3 模型假设. 4 4 符号说明. 5 5 模型的建立与求解. 5 5.1 问题一的模型. 5 5.1.1 模型建立. 5 模型 I. 6 5.1.2 模型求解. 6 5.2 问题二的模型. 7 5.2.1 模型建立. 7 模型 II . 8 5.2.2 模型求解. 9 5.3 问题三的模型. 10 5.3.1 模型建立. 10 模型 III. 10 5.3.2 模型求解. 11 6 模型的评价与改进. 12 附录. 13 A.1 附录描述 1. 13 A.2 附录描述 2. 20 A.3 附录描述 3. 23 摘要 1)DNA序列矩阵 a)首先调用函数I _importdata，将txt文档中的字符读入到workspace中，将 其中不是ATGC的字符用A来替换，存到矩阵DNA中。 b)给定DNA序列S s 1s2 L s N ，s i ACGT，A,T,G,C分别表示序列S的 前i(i 12 N)个元素中ACGT出现的次数.编写函数I _calculate计算 出DNA序列每一行的A,T,G,C。 c)根据定义的 欧几里得距离矩阵，为了求解ED方便 ,写编写一个函数 I _ ED1，求出四个矩阵Aij,Tij,Gij,Cij，然后编写函数I _ ED，便可以求出 AG，CT，AC，GT，AT，GC对应的ED。 d)根据定义的路径距离矩阵，编写函数I _ PD，根据不同的ED，求出不同的 PD。 e)根据定义的商矩阵， 编写函数I _EP， 由(b)和(c)求出的ED和PD计算出对 应的EP。 f)根据定义的商矩阵，编写函数I _ EG，由(b)求出的ED计算出对应的EG。 2)最大特征矩阵 对于给定的DNA序列， 令 ( 1 2 3 4 5 6 ),其中 1 2 3 4 5 6 分 别表示矩阵ED,EP,EG的标准化( i N)的最大特征根。将求得的结果放于 矩阵EGa1，EGa2，EGa3中。 3)调用给定程序作出DNA序列的聚类树图 可以根据给定的程序I _dendrogram，对于EGa1，EGa2，EGa3分别作图。 4)MATLAB的一些基础知识 关键字：DNA序列、最大特征根、聚类树图、MATLAB基础知识 1 问题重述 题目给的是关于DNA序列的问题 ,首先给出一个txt文档，是构成DNA序列的 ATGC字符，但其中有不是ATGC的字符，要求将txt文档中的所有字符序列读入到矩 阵DNA中， 并将其中不是ATGC的字符用ATGC中的任意一个换掉， 得到正确的DNA矩 阵。然后分别计算矩阵A,T,G,C (例如Ti表示前i个字符中含有的T的个数)。再根据 A,T,G,C四个矩阵不同的输入顺序求出 6 个不同的ED矩阵，对每个ED矩阵，可以求 出对应的PD,再由ED,PD求出对应的EP。对每个ED矩阵，也可以求出对应的EG矩 阵。 假设输入DNA序列的第一行，求出 6 个不同的ED,EP,EG后，每个矩阵可以求得一 个最大特征根，这样 6 个ED有 6 个最大特征根，放于矩阵EGa1 (EGa1为 24*6 的矩 阵)的第一行中，6 个EP有 6 个最大特征根，放于矩阵EGa2 (EGa2为 24*6 的矩阵) 的第一行中，6 个EG有 6 个最大特征根，放于矩阵EGa3 (EGa3为 24*6 的矩阵)的第 一行中。当把DNA序列的 24 行遍历完，EGa1,EGa2,EGa3三个矩阵便赋值完成。 提示中给了画DNA序列的聚类树图的程序，调用此程序，便可完成 3 个图的绘制。 2 问题分析 首先调用函数I _importdata，将txt文档中的字符读入到矩阵DNA中，将其中不是 ATGC的字符用A来替换。编写函数I _calculate计算出DNA序列每一行的A,T,G,C。 根据定义的欧几里得距离矩阵,为了求解ED方便,先由函数I _ ED1，求出四个矩阵 Aij,Tij,Gij,Cij，然后由函数I _ ED，便可以求出AG,CT,AC,GT,AT,GC对应的ED。根 据定义的路径距离矩阵，由函数I _ PD，根据不同的ED，求出不同的PD。根据定义 的商矩阵，由函数I _EP和求出的ED和PD计算出对应的EP。根据定义的商矩阵，由 函数I _ EG求出的ED计算出对应的EG。 编写函数I _ EGa1，假设输入 DNA 序列的第一行，求出 6 个不同的ED,EP,EG后， 每个矩阵对应一个最大特征根，6 个ED的 6 个最大特征根放于矩阵EGa1的第一行中， 6 个EP的 6 个最大特征根放于矩阵EGa2的第一行中，6 个EG的 6 个最大特征根放于 矩阵EGa3的第一行中。 遍历完DNA序列的 24 行， 得到EGa1,EGa2,EGa3三个最大特征 根矩阵。 由提示中给定的画DNA序列聚类树图的程序画出对应的聚类树图。 3 模型假设 1) 假设实验数据足够准确。 2) 假设实验过程中不存在仪器误差。 3) 假设提示中给的程序完整且正确。 4 符号说明 表 1 符号 DNA A T G C 描述 矩阵，存放从文档txt中取出的数据，并将非正常字符转换。 向量，其中Ai表示序列前i个中含有的A的个数 向量，其中Ti表示序列前i个中含有的T的个数 向量，其中Gi表示序列前i个中含有的G的个数 向量，其中Ci表示序列前i个中含有的C的个数 矩阵，A ij A j A i 矩阵，T ij T j T i 矩阵，Gij G j G i 矩阵，Cij C j C i 矩阵，根据ACGT输入顺序的不同，产生不同的ED 矩阵，PD ji PD ij ED ii1 ED i1i2 ED j1 j 矩阵，EPijE PijEDijPDiji j EPii 0 矩阵，EGijE GijEDij i j i j EGii 0 矩阵，存放ED的最大特征根 矩阵，存放EP的最大特征根 矩阵，存放EG的最大特征根 A ij T ij G ij C ij ED PD EP EG EGa1 EGa2 EGa3 5 模型的建立与求解 5.1 问题一的模型 读取数据到矩阵DNA中，并将其中不是ATGC的字符用A替换掉，得到处理过的 DNA矩阵。并根据输入的DNA序列的每一行，计算出 4 个矩阵A,T,G,C。 5.1.1 模型建立 首先，txt文档中的数据未经过处理，必须编写一个函数I _importdata，将txt中的 所有字符读入到工作空间中，其中共有 24 个序列，每一个序列都以开始，此时函数 中要有判断语句， 如果遇到就将它的下一行序列读入到矩阵DNA中， 读入完成以后， 要运用一段程序，判断其中的字符若不是ATGC，就用A替换。 这样就很清楚的得到了DNA序列矩阵。 提示中给定的A i C i G i ,T i 分别表示序列S的前i(i 12 N)个元素中ACGT出 现的次数。编写函数I _calculate，根据DNA中某一行的输入，得出需要的ATGC矩阵。 模型 I 1. 读取部分： P=importdata(TF24.txt); %读取 TF 中的数据到 a 中，P 是 cell 类型的列 向量 2. 判断部分： for i=1:m if double(Qi(1)=62 A(j)=i; %A 数组存放的的是首字符为的行数 j=j+1; end end 3. 取数据部分： for i=1:n DNAi=cell2mat(Q(A(i)+1:A(i+1)-1);%将介于两个之间的 段取出，然后转置，然后分开成单个字符，放于 cell 的 第一个中 end 4. 替换非正常数据部分： if DNAi(j)=A DNAi(j)=A; disp(转换完成！); end 5.1.2 模型求解 部分结果展示 图表 1 DNA 矩阵 图表 2输入DNA的第一行后得出的矩阵A的部分结果 5.2 问题二的模型 要求解ED,EP,EG的最大特征根矩阵，首先要把ED求解出来，而在求解ED时， 涉及到AG，CT，AC，GT，AT，GC这 6 个，而在求解它们是还有Aij,Tij,Gij,Cij这 4 个 矩阵，AG，CT，AC，GT，AT，GC的不同只是因为Aij,Tij,Gij,Cij输入顺序的不同而产 生的。所以先设计函数把Aij,Tij,Gij,Cij求解出来。再由函数求解不同的ED，根据求解 出来的ED，解相应的PD,EP,EG就变得简单多了。取出其中的ED及其对应的EP,EG， 分别求它们的最大特征根矩阵，放到三个数组EGa1,EGa2,EGa3中，再根据给定的作图 程序，做出DNA序列的聚类树图。 5.2.1 模型建立 1.为了求解ED方便 ,写编写一个函数I _ ED1，两层for 循环遍历求出 Aij,Tij,Gij,Cij。 2.由 函 数I _ ED， 其 中 输 入4个 形 参 ， 便 可 以 求 出 由函数I _ PD， 根据不同的ED，AG，CT，AC，GT，AT，GC对应的ED。 求出不同的PD，我们注意到，PD是对称矩阵，为了减少运算量，只求取 PD的上三角矩阵便可。 3.由函数I _EP，由于已经求出ED和PD，便可计算出对应的EP。此时求 解中注意到有除法运算， 而且最后对角线上全是 0。 为了避免除数为 0 的情 况，将对角阵的值全部变为 1，运算完成以后再将对角线赋值为 0 即可。 4.由函数I _ EG，根据求出的ED计算出对应的EG，只需要一个if else判 断即可。 5.由函数I _ EGa1， 求解ED,EP,EG的最大特征根矩阵， 将求得的结果放于矩 阵EGa1，EGa2，EGa3中。MATLAB有自带的求最大特征根的函数eigs， 只是它会求解出一列数据，用函数max取出最大的即可。EP,EG的最大特 征根求解中，只需要把I _ EGa1需要的 6 个参数换成相应的EP1EP6， EG1EG6即可。 6.调用给定程序作出DNA序列的聚类树图，调用过程中，变化的只是特征根 矩阵。 模型 II 1.I _ ED1 for i=1:n for j=1:n Aij(i,j)=A(j)-A(i); Tij(i,j)=T(j)-T(i); Gij(i,j)=G(j)-G(i); Cij(i,j)=C(j)-C(i); end end 2.I _ ED function ED=I_ED(x,y,z,w) %x y z w 是 Aij,Tij.Gij,Cij 中的一个 ED=sqrt(x-2*y).2+(z-2*w).2); 3.I _EP function EP=I_EP(ED,PD) n=length(PD); for i=1:n PD(i,i)=1; end EP=ED./PD; for i=1:n EP(i,i)=0; end 4.I _ EG function EG=I_EG(ED) n=length(ED); EG=zeros(n); for i=1:n for j=1:n if i=j EG(i,j)=ED(i,j)/abs(i-j); else EG(i,j)=0; end end end _ IEGa15. function y=I_EGa1(ED1,ED2,ED3,ED4,ED5,ED6) N=length(ED1); y(1)=max(eigs(ED1)/N); y(2)=max(eigs(ED2)/N); y(3)=max(eigs(ED3)/N); y(4)=max(eigs(ED4)/N); y(5)=max(eigs(ED5)/N); y(6)=max(eigs(ED6)/N); 6. _ Idendrogram EGe=zeros(24); for i=1:24 for j=i+1:24 EGe(i,j)=norm(EGa1(i,:)-EGa1(j,:); EGe(j,i)=EGe(i,j); end end EGe0=pdist(EGe,euclidean); % Pairwise distance EGtree=linkage(EGe0); % Create hierarchical cluster tree EPtreeShow=dendrogram(EGtree,orientation,left);% Dendrogram plot 5.2.2 模型求解 部分结果显示： A ij 图表 3图表 4 ED1 图表 5图表 6 EP1?EG1 图表 7 EGa1 图表 8 EGa2 5.3 问题三的模型 由于DNA序列的第 20 行数据比其他行多近一半的数据，一般的电脑运行时会出现 溢出现象，为了解决这个问题，第二十行的程序需要加以改进，使其能够适应一般的电 脑运行内存。 5.3.1 模型建立 运行时的问题从求解出ACGT以后就出现了。所以ED的函数不能用，现在我考虑，可 能是先求解出Aij,Tij,Gij,Cij导致占用了过多内存，我现在把A ij A j A i ,C ij C j C i , G ij G j G i ,T ij T j T i 的原始公式直接带入编写的函数中，发现运行速度变快。但是 这样6个ED需要依次求出， 然后依次运行PD1 I _PDED1，EP1 I _EPED1,PD1， u2 maxeigsEP1/ n，EG1 I _EGED1，u3 maxeigsEG1/ n，一遍运行完成 后，得到3个最大特征根矩阵第20行中的一个，运行6次，3个矩阵的第20行变求出来了。 模型 III 1. 求解EGa1第 20 行第一个值u1 n=length(DNA20); ED1=zeros(n,n); for i=1:n for j=1:n ED1(i,j)=sqrt(A(j)-A(i)-2*(G(j)-G(i)2+(T(j)-T(i) -2*(C(j)-C(i)2); end end u1=max(eigs(ED1)/n); 2. 求解EGa2第 20 行第一个值u2 PD1=I_PD(ED1); EP1=I_EP(ED1,PD1); u2=max(eigs(EP1)/n); 3. 求解EGa3第 20 行第一个值u3 EG1=I_EG(ED1) ; u3=max(eigs(EG1)/n); clear(ED1,PD1,EP1,EG1,u1,u2,u3); 5.3.2 模型求解 第 20 行结果展示： 图表 9 EGa1的 20 行 图表 10 EGa2的 20 行 图表 11 EGa3的 20 行 DNA序列聚类树图： 图表 12 EGa1的聚类树图 图表 13 EGa2的聚类树图 图表 14 EGa3的聚类树图 6 模型的评价与改进 此文章能较好的将DNA序列进行分类和鉴别，运行最终的主程序,便可清楚地给出 其聚类树图和相应的最大特征根矩阵。但是其在处理大量DNA序列的过程中，容易出 现溢出现象，表明程序仍然是不够完善的，毕竟编程过程中用到的两层for循环对于程 序的影响是很大的，可以考虑将其中的for循环尽量的简化，但是我还是没有更好的方 法。根据我对第 20 行的重新编码，可以处理的数据也可以很大，但是需要手动的运算 量增加，整体运行时间较长,消耗的计算机内存较大,这些都是需要改进的地方。 参考文献 1 张养安、刘林章，MATLAB软件在数据处理中的应用，水土保持研究，2005 年 8 月 2 姜长元，MATLAB和 Word的链接及应用，洛阳师范学院，2005 年第二期 3 王沫然，MATLAB与科学计算（第二版） ，电子出版社，2008 年 9 月。 附录 A.1 附录描述 1 问题一中的程序： P=importdata(TF24.txt); %读取 TF 中的数据到 a 中，P 是 cell 类型的列向量 Q=P;%将取出来的数据放到 Q 中 DNA=cell(24,1); m=length(Q); A=zeros(1,25); j=1;k=1; A(25)=m+1; for i=1:m if double(Qi(1)=62 A(j)=i; %A 数组存放的的是首字符为的行数 j=j+1; end end n=length(A)-1; for i=1:n DNAi=cell2mat(Q(A(i)+1:A(i+1)-1);%将介于两个之间的段取出， 然后转置，然后分开成单个字符，放于 cell 的第一个中 end % %找到不是 ATGC 的字符替换 % for i=1:24 for j=1:length(DNAi) if DNAi(j)=A i,j DNAi(j)=A; disp(转换完成！); end end end function A T G C=I_calculate(X) % X 是输入的数组,即 DNAi i=1,2,3.24 % A T G C 是输出数据,是行向量，A(i)表示前 i 个字符中含有的 A 的个数 % 刘海鹏 2012.08.25 A=zeros(1,length(X); T=zeros(1,length(X); G=zeros(1,length(X); C=zeros(1,length(X); A(1)=0; T(1)=0; G(1)=0; C(1)=0; if X(1)=A A(1)=A(1)+1; elseif X(1)=T T(1)=T(1)+1; elseif X(1)=G G(1)=G(1)+1; else C(1)=C(1)+1; end for i=2:length(X) if X(i)=A A(i)=A(i-1)+1; T(i)=T(i-1); G(i)=G(i-1); C(i)=C(i-1); elseif X(i)=T A(i)=A(i-1); T(i)=T(i-1)+1; G(i)=G(i-1); C(i)=C(i-1); elseif X(i)=G A(i)=A(i-1); T(i)=T(i-1); G(i)=G(i-1)+1; C(i)=C(i-1); else A(i)=A(i-1); T(i)=T(i-1); G(i)=G(i-1); C(i)=C(i-1)+1; end end % z1=0;z2=0;z3=0;z4=0; % for i=1:length(X) % if X(i)=A % z1=z1+1; % elseif X(i)=T % z2=z2+1; % elseif X(i)=G % z3=z3+1; % else z4=z4+1; % end % A(i)=z1; % T(i)=z2; % G(i)=z3; % C(i)=z4; % end 问题二中的程序 function Aij,Tij,Cij,Gij=I_ED1(A,T,G,C) n=length(A); Aij=zeros(n); Tij=zeros(n); Cij=zeros(n); Gij=zeros(n); for i=1:n for j=1:n Aij(i,j)=A(j)-A(i); Tij(i,j)=T(j)-T(i); Gij(i,j)=G(j)-G(i); Cij(i,j)=C(j)-C(i); end end function ED=I_ED(x,y,z,w) %x y z w 是 Aij,Tij.Gij,Cij 中的一个 ED=sqrt(x-2*y).2+(z-2*w).2); % ED1=I_ED(Aij,Gij,Tij,Cij); % ED2=I_ED(Cij,Tij,Aij,Gij); % ED3=I_ED(Aij,Cij,Tij,Gij); % ED4=I_ED(Gij,Tij,Aij,Cij); % ED5=I_ED(Aij,Tij,Cij,Gij); % ED6=I_ED(Gij,Cij,Aij,Tij); function PD=I_PD(ED) n=length(ED); for i=1:n for j=i+2:n ED(i,j)=ED(i,j-1)+ED(j-1,j); end ED(i,i)=0; end ED=triu(ED); PD=ED+ED; % PD1=I_PD(ED1); % PD2=I_PD(ED2); % PD3=I_PD(ED3); % PD4=I_PD(ED4); % PD5=I_PD(ED5); % PD6=I_PD(ED6); function EP=I_EP(ED,PD) n=length(PD); for i=1:n PD(i,i)=1; end EP=ED./PD; for i=1:n EP(i,i)=0; end % EP1=I_EP(ED1,PD1); % EP2=I_EP(ED2,PD2); % EP3=I_EP(ED3,PD3); % EP4=I_EP(ED4,PD4); % EP5=I_EP(ED5,PD5); % EP6=I_EP(ED6,PD6); function EG=I_EG(ED) n=length(ED); EG=zeros(n); for i=1:n for j=1:n if i=j EG(i,j)=ED(i,j)/abs(i-j); else EG(i,j)=0; end end end % EG1=I_EG(ED1); % EG2=I_EG(ED2); % EG3=I_EG(ED3); % EG4=I_EG(ED4); % EG5=I_EG(ED5); % EG6=I_EG(ED6); function y=I_EGa1(ED1,ED2,ED3,ED4,ED5,ED6) N=length(ED1); y(1)=max(eigs(ED1)/N); y(2)=max(eigs(ED2)/N); y(3)=max(eigs(ED3)/N); y(4)=max(eigs(ED4)/N); y(5)=max(eigs(ED5)/N); y(6)=max(eigs(ED6)/N); % EGa1(,1:6)=I_EGa1(ED1,ED2,ED3,ED4,ED5,ED6); % EGa2(,1:6)=I_EGa1(EP1,EP2,EP3,EP4,EP5,EP6); % EGa3(,1:6)=I_EGa1(EG1,EG2,EG3,EG4,EG5,EG6); EGe=zeros(24); for i=1:24 for j=i+1:24 EGe(i,j)=norm(EGa1(i,:)-EGa1(j,:); EGe(j,i)=EGe(i,j); end end EGe0=pdist(EGe,euclidean); % Pairwise distance EGtree=linkage(EGe0); % Create hierarchical cluster tree EPtreeShow=dendrogram(EGtree,orientation,left);% Dendrogram plot 问题 3 中的程序(第 20 行) n=length(DNA20); ED1=zeros(n,n); for i=1:n for j=1:n ED1(i,j)=sqrt(A(j)-A(i)-2*(G(j)-G(i)2+(T(j)-T(i)-2*(C(j )-C(i)2); end end u1=max(eigs(ED1)/n); PD1=I_PD(ED1); EP1=I_EP(ED1,PD1); u2=max(eigs(EP1)/n); EG1=I_EG(ED1) ; u3=max(eigs(EG1)/n); EGa1(20,1)=u1;EGa2(20,1)=u2;EGa3(20,1)=u3; clear(ED1,PD1,EP1,EG1,u1,u2,u3); % % ED2=CT % ED2=zeros(n); for i=1:n for j=1:n ED2(i,j)=sqrt(C(j)-C(i)-2*(T(j)-T(i)2+(A(j)-A(i)-2*(G(j )-G(i)2); end end u1=max(eigs(ED2)/n); PD2=I_PD(ED2); EP2=I_EP(ED2,PD2); u2=max(eigs(EP2)/n); EG2=I_EG(ED2) ; u3=max(eigs(EG2)/n); EGa1(20,2)=u1;EGa2(20,2)=u2;EGa3(20,2)=u3; clear(ED2,PD2,EP2,EG2,u1,u2,u3); % % ED3=AC % ED1=AG ED3=zeros(n); for i=1:n for j=1:n ED3(i,j)=sqrt(A(j)-A(i)-2*(C(j)-C(i)2+(T(j)-T(i)-2*(G( j)-G(i)2); end end u1=max(eigs(ED3)/n); PD3=I_PD(ED3); EP3=I_EP(ED3,PD3); u2=max(eigs(EP3)/n); EG3=I_EG(ED3) ; u3=max(eigs(EG3)/n); EGa1(20,3)=u1;EGa2(20,3)=u2;EGa3(20,3)=u3; clear(ED3,PD3,EP3,EG3,u1,u2,u3); % % ED4=GT % % ED1=AG ED4=zeros(n); for i=1:n for j=1:n ED4(i,j)=sqrt(G(j)-G(i)-2*(T(j)-T(i)2+(A(j)-A(i)-2*(C(j )-C(i)2); end end u1=max(eigs(ED4)/n); PD4=I_PD(ED4); EP4=I_EP(ED4,PD4); u2=max(eigs(EP4)/n); EG4=I_EG(ED4) ; u3=max(eigs(EG4)/n); EGa1(20,4)=u1;EGa2(20,4)=u2;EGa3(20,4)=u3; clear(ED4,PD4,EP4,EG4,u1,u2,u3); % % ED5=AT % ED1=AG ED5=zeros(n); for i=1:n for j=1:n ED5(i,j)=sqrt(A(j)-A(i)-2*(T(j)-T(i)2+(C(j)-C(i)-2*(G( j)-G(i)2); end end u1=max(eigs(ED5)/n); PD5=I_PD(ED5); EP5=I_EP(ED5,PD5); u2=max(eigs(EP5)/n); EG5=I_EG(ED5) ; u3=max(eigs(EG5)/n); EGa1(20,5)=u1;EGa2(20,5)=u2;EGa3(20,5)=u3; clear(ED5,PD5,EP5,EG5,u1,u2,u3); % % ED6=GC % ED6=zeros(n); for i=1:n for j=1:n ED6(i,j)=sqrt(G(j)-G(i)-2*(C(j)-C(i)2+(A(j)-A(i)-2*(T(j )-T(i)2); end end u1=max(eigs(ED6)/n); PD6=I_PD(ED6); EP6=I_EP(ED6,PD6); u2=max(eigs(EP6)/n); EG6=I_EG(ED6) ; u3=max(eigs(EG6)/n); EGa1(20,6)=u1;EGa2(20,6)=u2;EGa3(20,6)=u3; clear(ED6,PD6,EP6,EG6,u1,u2,u3,C,G,T,A); A.2 附录描述 2 主程序 clear; clc; profile on; profile clear; DNA=cell(24,1); I_importdata; clear(P,Q,m); EGa1=zeros(24,6); EGa2=zeros(24,6); EGa3=zeros(24,6); for i=1:19 A T G C=I_calculate(DNAi); Aij,Tij,Cij,Gij=I_ED1(A,T,G,C); clear(A,T,G,C); ED1=I_ED(Aij,Gij,Tij,Cij); ED2=I_ED(Cij,Tij,Aij,Gij); ED3=I_ED(Aij,Cij,Tij,Gij); ED4=I_ED(Gij,Tij,Aij,Cij); ED5=I_ED(Aij,Tij,Cij,Gij); ED6=I_ED(Gij,Cij,Aij,Tij); clear(Aij,Cij,Gij,Tij); EGa1(i,:)=I_EGa1(ED1,ED2,ED3,ED4,ED5,ED6); PD1=I_PD(ED1); PD2=I_PD(ED2); PD3=I_PD(ED3); PD4=I_PD(ED4); PD5=I_PD(ED5); PD6=I_PD(ED6); EP1=I_EP(ED1,PD1); EP2=I_EP(ED