001/* This class is based on the original FATCAT implementation by
002 * <pre>
003 * Yuzhen Ye & Adam Godzik (2003)
004 * Flexible structure alignment by chaining aligned fragment pairs allowing twists.
005 * Bioinformatics vol.19 suppl. 2. ii246-ii255.
006 * https://www.ncbi.nlm.nih.gov/pubmed/14534198
007 * </pre>
008 *
009 * Thanks to Yuzhen Ye and A. Godzik for granting permission to freely use and redistribute this code.
010 *
011 * This code may be freely distributed and modified under the
012 * terms of the GNU Lesser General Public Licence.  This should
013 * be distributed with the code.  If you do not have a copy,
014 * see:
015 *
016 *      http://www.gnu.org/copyleft/lesser.html
017 *
018 * Copyright for this code is held jointly by the individual
019 * authors.  These should be listed in @author doc comments.
020 *
021 *
022 * Created on Jun 17, 2009
023 * Created by Andreas Prlic - RCSB PDB
024 *
025 */
026
027package org.biojava.nbio.structure.align.fatcat.calc;
028
029import org.biojava.nbio.structure.Atom;
030import org.biojava.nbio.structure.StructureException;
031import org.biojava.nbio.structure.align.model.AFP;
032import org.biojava.nbio.structure.align.model.AFPChain;
033
034import java.util.List;
035
036public class AFPOptimizer
037{
038
039        public static final boolean debug = FatCatAligner.debug;
040
041        /**
042         * optimize the alignment by dynamic programming
043         */
044        public static void optimizeAln(FatCatParameters params, AFPChain afpChain,Atom[] ca1, Atom[] ca2) throws StructureException
045        {
046
047                int minLen = afpChain.getMinLen();
048                int fragLen = params.getFragLen();
049
050
051                long optStart = System.currentTimeMillis();
052                int     i, a, k, p1, p2, bk, b1, b2, e1, e2, a1, a2;
053                int     iniLen;
054                int[][]     iniSet = new int[2][minLen];
055                int     maxi = 100;
056
057                int[][][] optAln = afpChain.getOptAln();
058                int[] optLen = afpChain.getOptLen();
059
060                int maxTra = params.getMaxTra();
061                double[] optRmsd = afpChain.getOptRmsd();
062
063                int blockNum = afpChain.getBlockNum();
064
065                if(optAln == null)      {
066                        optAln     = new int[maxTra+1][2][minLen];
067                        afpChain.setOptAln(optAln);
068                        optLen     = new int[maxTra+1];
069                        afpChain.setOptLen(optLen);
070                        optRmsd    = new double[maxTra+1];
071                        afpChain.setOptRmsd(optRmsd);
072                }
073
074                List<AFP> afpSet = afpChain.getAfpSet();
075
076                int optLength         = afpChain.getOptLength();
077                int[] afpChainList    = afpChain.getAfpChainList();
078                int[] block2Afp       = afpChain.getBlock2Afp();
079                int[] blockSize       = afpChain.getBlockSize();
080
081
082                if (debug)
083                        System.out.println("AFPOptimizer got blockNum: " +blockNum);
084                //optimize each alignment defined by a block
085                b1 = b2 = e1 = e2 = optLength = 0;
086                for(bk = 0; bk < blockNum; bk ++)       {
087                        //initial aligned position
088                        iniLen = 0;
089                        if(bk > 0)      {
090                                b1 = e1;
091                                b2 = e2;
092                        }
093                        if(bk < blockNum - 1)   {
094                                a1 = afpChainList[block2Afp[bk] + blockSize[bk] - 1]; //the last AFP in current block
095                                a2 = afpChainList[block2Afp[bk + 1]];  //the first AFP in next block
096                                e1 = (afpSet.get(a1).getP1() + fragLen +  afpSet.get(a2).getP1()) / 2;
097                                e2 = (afpSet.get(a1).getP2() + fragLen +  afpSet.get(a2).getP2()) / 2;
098                        } //use the middle point of the current and next AFPs. old (starting point of next AFP)
099                        else    {
100                                e1 = ca1.length;
101                                e2 = ca2.length;
102                        }
103
104
105                        for(i = block2Afp[bk]; i < block2Afp[bk] + blockSize[bk]; i ++) {
106                                a = afpChainList[i];
107                                p1 = afpSet.get(a).getP1();
108                                p2 = afpSet.get(a).getP2();
109                                for(k = 0; k < afpSet.get(a).getFragLen(); k ++)     {
110                                        iniSet[0][iniLen] = p1 + k - b1; //note -b1
111                                        iniSet[1][iniLen] = p2 + k - b2; //note -b2
112                                        iniLen ++;
113
114                                }
115                        }
116                        //optimize the align by dynamic programming & constraint the optimization region
117                        if(debug) {
118                                System.err.println(String.format("optimize block %d (%d afp), region %d-%d(len %d), %d-%d(len %d)\n",
119                                                bk, blockSize[bk], b1, e1, e1-b1, b2, e2, e2-b2));
120
121                                System.err.println(" initial alignment Length: " + iniLen );
122                        }
123
124                        StructureAlignmentOptimizer opt = new StructureAlignmentOptimizer(b1,e1, ca1, b2,e2, ca2, iniLen, iniSet);
125                        opt.runOptimization(maxi);
126                        optRmsd[bk] = opt.optimizeResult(optLen,bk,optAln[bk]);
127
128                        //System.out.println(optRmsd[bk]);
129                        // SALNOPT *opt = new SALNOPT(e1-b1, &pro1->caCod[3 * b1], e2-b2, &pro2->caCod[3 * b2], iniLen, iniSet, maxi);
130                        // optRmsd[bk] = opt->OptimizeResult(&optLen[bk], optAln[bk]);
131
132                        if(debug)
133                                System.out.println(String.format(" optimized len=%d, rmsd %f\n", optLen[bk], optRmsd[bk]));
134
135                        for(i = 0; i < optLen[bk]; i ++)        {
136                                optAln[bk][0][i] += b1; //restore the position
137                                optAln[bk][1][i] += b2; //restore the position
138                        }
139                        optLength += optLen[bk];
140
141
142                }
143
144                long optEnd = System.currentTimeMillis();
145                if(debug)       System.out.println("complete AlignOpt " + (optEnd-optStart) +"\n");
146
147                afpChain.setBlockNum(blockNum);
148                afpChain.setOptLength(optLength);
149                afpChain.setAfpChainList(afpChainList);
150                afpChain.setBlock2Afp(block2Afp);
151                afpChain.setBlockSize(blockSize);
152
153
154        }
155
156        /**
157         * get the afp list and residue list for each block
158         */
159
160        public static void blockInfo(AFPChain afpChain)
161        {
162                int     i, j, k, a, n;
163
164                int blockNum = afpChain.getBlockNum();
165
166                int[] blockSize =afpChain.getBlockSize();
167                int[] afpChainList = afpChain.getAfpChainList();
168                int[] block2Afp = afpChain.getBlock2Afp();
169                int[][][]blockResList = afpChain.getBlockResList();
170
171                List<AFP>afpSet = afpChain.getAfpSet();
172                int[] blockResSize = afpChain.getBlockResSize();
173
174                for(i = 0; i < blockNum; i ++)  {
175                        n = 0;
176                        for(j = 0; j < blockSize[i]; j ++)      {
177                                //the index in afpChainList, not in the whole afp set
178                                a = afpChainList[block2Afp[i] + j];
179                                for(k = 0; k < afpSet.get(a).getFragLen(); k ++)     {
180                                        blockResList[i][0][n] = afpSet.get(a).getP1() + k;
181                                        blockResList[i][1][n] = afpSet.get(a).getP2() + k;
182                                        n ++;
183                                }
184                        }
185                        blockResSize[i] = n;
186                }
187
188                afpChain.setBlockResSize(blockResSize);
189                afpChain.setBlockSize(blockSize);
190                afpChain.setAfpChainList(afpChainList);
191                afpChain.setBlock2Afp(block2Afp);
192                afpChain.setBlockResList(blockResList);
193        }
194
195        /**
196         * to update the chaining score after block delete and merge processed
197         * the blockScore value is important for significance evaluation
198         */
199        public static void updateScore(FatCatParameters params, AFPChain afpChain)
200        {
201                int     i, j, bknow, bkold, g1, g2;
202
203
204                afpChain.setConn(0d);
205                afpChain.setDVar(0d);
206
207                int blockNum = afpChain.getBlockNum();
208                int alignScoreUpdate = 0;
209                double[] blockScore = afpChain.getBlockScore();
210                int[] blockGap = afpChain.getBlockGap();
211                int[] blockSize =afpChain.getBlockSize();
212                int[] afpChainList = afpChain.getAfpChainList();
213                List<AFP>afpSet = afpChain.getAfpSet();
214                int[] block2Afp = afpChain.getBlock2Afp();
215
216                double torsionPenalty = params.getTorsionPenalty();
217
218
219                bkold = 0;
220                for(i = 0; i < blockNum; i ++)  {
221                        blockScore[i] = 0;
222                        blockGap[i] = 0;
223                        for(j = 0; j < blockSize[i]; j ++)      {
224                                bknow = afpChainList[block2Afp[i] + j];
225                                if(j == 0)      {
226                                        blockScore[i] = afpSet.get(bknow).getScore();
227                                }
228                                else    {
229                                        AFPChainer.afpPairConn(bkold, bknow, params, afpChain); //note: j, i
230                                        Double conn = afpChain.getConn();
231                                        blockScore[i] += afpSet.get(bknow).getScore() + conn;
232                                        g1 = afpSet.get(bknow).getP1() - afpSet.get(bkold).getP1() - afpSet.get(bkold).getFragLen();
233                                        g2 = afpSet.get(bknow).getP2() - afpSet.get(bkold).getP2() - afpSet.get(bkold).getFragLen();
234                                        blockGap[i] += (g1 > g2)?g1:g2;
235                                }
236                                bkold = bknow;
237                        }
238                        alignScoreUpdate += blockScore[i];
239                }
240                if(blockNum >= 2)       {
241                        alignScoreUpdate += (blockNum - 1) * torsionPenalty;
242                }
243
244                afpChain.setBlockGap(blockGap);
245                afpChain.setAlignScoreUpdate(alignScoreUpdate);
246                afpChain.setBlockScore(blockScore);
247                afpChain.setBlockSize(blockSize);
248                afpChain.setAfpChainList(afpChainList);
249                afpChain.setBlock2Afp(block2Afp);
250        }
251
252
253
254}