Java Code Examples for org.biojava.nbio.structure.align.model.AFPChain#getBlockNum()

public void setAFPChain(AFPChain afpChain) {

		this.afpChain = afpChain;
		coordManager.setAFPChain(afpChain);
		if ( afpChain != null) {
			selection = new BitSet (afpChain.getAlnLength());
			if ( afpChain.getBlockNum() > 1) {
				colorByAlignmentBlock = true;
			}
		}

	}
 

/**
 * Calculate the rotation angle for a structure
 * @param afpChain
 * @return The rotation angle, in radians
 * @throws StructureException If the alignment doesn't contain any blocks
 * @throws NullPointerException If the alignment doesn't have a rotation matrix set
 */
public static double getAngle(AFPChain afpChain) throws StructureException {
	if(afpChain.getBlockNum() < 1) {
		throw new StructureException("No aligned residues");
	}
	Matrix rotation = afpChain.getBlockRotationMatrix()[0];

	if(rotation == null) {
		throw new NullPointerException("AFPChain does not contain a rotation matrix");
	}
	return getAngle(rotation);
}
 

public double calNS(FatCatParameters params, AFPChain afpChain){
	int len1 = afpChain.getCa1Length();
	int len2 = afpChain.getCa2Length();
	double score =afpChain.getAlignScore();
	double rmsd = afpChain.getTotalRmsdOpt();
	int optLen  = afpChain.getOptLength();
	int r = afpChain.getBlockNum() -1;
	return calNS(len1, len2,score,rmsd,optLen,r);
}
 

/**
 * Fill the aligned Atom arrays with the equivalent residues in the afpChain.
 * @param afpChain
 * @param ca1
 * @param ca2
 * @param ca1aligned
 * @param ca2aligned
 */
public static void fillAlignedAtomArrays(AFPChain afpChain, Atom[] ca1,
		Atom[] ca2, Atom[] ca1aligned, Atom[] ca2aligned) {

	int pos=0;
	int[] blockLens = afpChain.getOptLen();
	int[][][] optAln = afpChain.getOptAln();
	assert(afpChain.getBlockNum() <= optAln.length);

	for (int block=0; block < afpChain.getBlockNum(); block++) {
		for(int i=0;i<blockLens[block];i++) {
			int pos1 = optAln[block][0][i];
			int pos2 = optAln[block][1][i];
			Atom a1 = ca1[pos1];
			Atom a2 = (Atom) ca2[pos2].clone();
			ca1aligned[pos] = a1;
			ca2aligned[pos] = a2;
			pos++;
		}
	}

	// this can happen when we load an old XML serialization which did not support modern ChemComp representation of modified residues.
	if (pos != afpChain.getOptLength()){
		logger.warn("AFPChainScorer getTMScore: Problems reconstructing alignment! nr of loaded atoms is " + pos + " but should be " + afpChain.getOptLength());
		// we need to resize the array, because we allocated too many atoms earlier on.
		ca1aligned = (Atom[]) resizeArray(ca1aligned, pos);
		ca2aligned = (Atom[]) resizeArray(ca2aligned, pos);
	}

}
 

/**
 * get the afp list and residue list for each block
 */

public static void blockInfo(AFPChain afpChain)
{
	int     i, j, k, a, n;

	int blockNum = afpChain.getBlockNum();

	int[] blockSize =afpChain.getBlockSize();
	int[] afpChainList = afpChain.getAfpChainList();
	int[] block2Afp = afpChain.getBlock2Afp();
	int[][][]blockResList = afpChain.getBlockResList();

	List<AFP>afpSet = afpChain.getAfpSet();
	int[] blockResSize = afpChain.getBlockResSize();

	for(i = 0; i < blockNum; i ++)  {
		n = 0;
		for(j = 0; j < blockSize[i]; j ++)      {
			//the index in afpChainList, not in the whole afp set
			a = afpChainList[block2Afp[i] + j];
			for(k = 0; k < afpSet.get(a).getFragLen(); k ++)     {
				blockResList[i][0][n] = afpSet.get(a).getP1() + k;
				blockResList[i][1][n] = afpSet.get(a).getP2() + k;
				n ++;
			}
		}
		blockResSize[i] = n;
	}

	afpChain.setBlockResSize(blockResSize);
	afpChain.setBlockSize(blockSize);
	afpChain.setAfpChainList(afpChainList);
	afpChain.setBlock2Afp(block2Afp);
	afpChain.setBlockResList(blockResList);
}
 

/**
 * Return a list of pdb Strings corresponding to the aligned positions of the molecule.
 * Only supports a pairwise alignment with the AFPChain DS.
 *
 * @param aligPos
 * @param afpChain
 * @param ca
 */
public static final List<String> getPDBresnum(int aligPos, AFPChain afpChain, Atom[] ca){
	List<String> lst = new ArrayList<String>();
	if ( aligPos > 1) {
		System.err.println("multiple alignments not supported yet!");
		return lst;
	}

	int blockNum = afpChain.getBlockNum();
	int[] optLen = afpChain.getOptLen();
	int[][][] optAln = afpChain.getOptAln();

	if ( optLen == null)
		return lst;

	for(int bk = 0; bk < blockNum; bk ++)       {

		for ( int i=0;i< optLen[bk];i++){

			int pos = optAln[bk][aligPos][i];
			if ( pos < ca.length) {
				String pdbInfo = JmolTools.getPdbInfo(ca[pos]);
				//lst.add(ca1[pos].getParent().getPDBCode());
				lst.add(pdbInfo);
			}
		}

	}
	return lst;
}
 

private static void reverseOptAln(AFPChain afpChain) {
	int ca2len = afpChain.getCa2Length();

	int[][][] optAln = afpChain.getOptAln();
	int[] optLen = afpChain.getOptLen();

	for (int block = 0; block < afpChain.getBlockNum(); block++) {
		for (int pos = 0; pos < optLen[block]; pos++) {
			optAln[block][1][pos] = ca2len - 1 - optAln[block][1][pos];
		}
	}

	afpChain.setOptAln(optAln);
}
 

/**
 *  Partitions an afpChain alignment into order blocks of aligned residues.
 *
 * @param afpChain
 * @param ca1
 * @param ca2
 * @param order
 * @return
 * @throws StructureException
 */
private static AFPChain partitionAFPchain(AFPChain afpChain,
		Atom[] ca1, Atom[] ca2, int order) throws StructureException{

	int[][][] newAlgn = new int[order][][];
	int repeatLen = afpChain.getOptLength()/order;

	//Extract all the residues considered in the first chain of the alignment
	List<Integer> alignedRes = new ArrayList<Integer>();
	for (int su=0; su<afpChain.getBlockNum(); su++){
		for (int i=0; i<afpChain.getOptLen()[su]; i++){
			alignedRes.add(afpChain.getOptAln()[su][0][i]);
		}
	}

	//Build the new alignment
	for (int su=0; su<order; su++){
		newAlgn[su] = new int[2][];
		newAlgn[su][0] = new  int[repeatLen];
		newAlgn[su][1] = new  int[repeatLen];
		for (int i=0; i<repeatLen; i++){
			newAlgn[su][0][i] = alignedRes.get(repeatLen*su+i);
			newAlgn[su][1][i] = alignedRes.get(
					(repeatLen*(su+1)+i)%alignedRes.size());
		}
	}

	return AlignmentTools.replaceOptAln(newAlgn, afpChain, ca1, ca2);
}
 

/**
 * Modifies the {@link AFPChain#setOptAln(int[][][]) optAln} of an AFPChain
 * by permuting the second protein.
 *
 * Sets residue numbers in the second protein to <i>(i-cp)%len</i>
 *
 * @param afpChain
 * @param cp Amount leftwards (or rightward, if negative) to shift the
 */
private static void permuteOptAln(AFPChain afpChain, int cp)
{
	int ca2len = afpChain.getCa2Length();

	if( ca2len <= 0) {
		throw new IllegalArgumentException("No Ca2Length specified in "+afpChain);
	}

	// Allow negative cp points for convenience.
	if(cp == 0) {
		return;
	}
	if(cp <= -ca2len || cp >= ca2len) {
		// could just take cp%ca2len, but probably its a bug if abs(cp)>=ca2len
		throw new ArrayIndexOutOfBoundsException( String.format(
				"Permutation point %d must be between %d and %d for %s",
				cp, 1-ca2len,ca2len-1, afpChain.getName2() ) );
	}
	if(cp < 0) {
		cp = cp + ca2len;
	}

	// the unprocessed alignment
	int[][][] optAln = afpChain.getOptAln();
	int[] optLen = afpChain.getOptLen();

	// the processed alignment
	List<List<List<Integer>>> blocks = new ArrayList<List<List<Integer>>>(afpChain.getBlockNum()*2);

	//Update residue indices
	// newi = (oldi-cp) % N
	for(int block = 0; block < afpChain.getBlockNum(); block++) {
		if(optLen[block]<1)
			continue;

		// set up storage for the current block
		List<List<Integer>> currBlock = new ArrayList<List<Integer>>(2);
		currBlock.add( new ArrayList<Integer>());
		currBlock.add( new ArrayList<Integer>());
		blocks.add(currBlock);

		// pos = 0 case
		currBlock.get(0).add( optAln[block][0][0] );
		currBlock.get(1).add( (optAln[block][1][0]+cp ) % ca2len);

		for(int pos = 1; pos < optLen[block]; pos++) {
			//check if we need to start a new block
			//this happens when the new alignment crosses the protein terminus
			if( optAln[block][1][pos-1]+cp<ca2len &&
					optAln[block][1][pos]+cp >= ca2len) {
				currBlock = new ArrayList<List<Integer>>(2);
				currBlock.add( new ArrayList<Integer>());
				currBlock.add( new ArrayList<Integer>());
				blocks.add(currBlock);
			}
			currBlock.get(0).add( optAln[block][0][pos] );
			currBlock.get(1).add( (optAln[block][1][pos]+cp ) % ca2len);
		}
	}

	// save permuted blocks to afpChain
	assignOptAln(afpChain,blocks);
}
 

/** Rotate the Atoms/Groups so they are aligned for the 3D visualisation
 *
 * @param afpChain
 * @param ca1
 * @param ca2
 * @return an array of Groups that are transformed for 3D display
 * @throws StructureException
 */
public static Group[] prepareGroupsForDisplay(AFPChain afpChain, Atom[] ca1, Atom[] ca2) throws StructureException{


	if ( afpChain.getBlockRotationMatrix().length == 0 ) {
		// probably the alignment is too short!
		System.err.println("No rotation matrix found to rotate 2nd structure!");
		afpChain.setBlockRotationMatrix(new Matrix[]{Matrix.identity(3, 3)});
		afpChain.setBlockShiftVector(new Atom[]{new AtomImpl()});
	}

	// List of groups to be rotated according to the alignment
	Group[] twistedGroups = new Group[ ca2.length];

	//int blockNum = afpChain.getBlockNum();

	int i = -1;

	// List of groups from the structure not included in ca2 (e.g. ligands)
	// Will be rotated according to first block
	List<Group> hetatms2 = StructureTools.getUnalignedGroups(ca2);

	if (  (afpChain.getAlgorithmName().equals(FatCatRigid.algorithmName) ) || (afpChain.getAlgorithmName().equals(FatCatFlexible.algorithmName) ) ){

		for (Atom a: ca2){
			i++;
			twistedGroups[i]=a.getGroup();

		}

		twistedGroups = AFPTwister.twistOptimized(afpChain, ca1, ca2);

		//} else  if  (( blockNum == 1 ) || (afpChain.getAlgorithmName().equals(CeCPMain.algorithmName))) {
	} else {

		Matrix m   =  afpChain.getBlockRotationMatrix()[ 0];
		Atom shift =  afpChain.getBlockShiftVector()   [ 0 ];

		shiftCA2(afpChain, ca2, m,shift, twistedGroups);

	}

	if ( afpChain.getBlockNum() > 0){

		// Superimpose ligands relative to the first block
		if( hetatms2.size() > 0 ) {

			if ( afpChain.getBlockRotationMatrix().length > 0 ) {

				Matrix m1      = afpChain.getBlockRotationMatrix()[0];
				//m1.print(3,3);
				Atom   vector1 = afpChain.getBlockShiftVector()[0];
				//System.out.println("shift vector:" + vector1);

				for ( Group g : hetatms2){
					Calc.rotate(g, m1);
					Calc.shift(g,vector1);
				}
			}
		}
	}

	return twistedGroups;
}
 

/**
 * to update the chaining score after block delete and merge processed
 * the blockScore value is important for significance evaluation
 */
public static void updateScore(FatCatParameters params, AFPChain afpChain)
{
	int     i, j, bknow, bkold, g1, g2;


	afpChain.setConn(0d);
	afpChain.setDVar(0d);

	int blockNum = afpChain.getBlockNum();
	int alignScoreUpdate = 0;
	double[] blockScore = afpChain.getBlockScore();
	int[] blockGap = afpChain.getBlockGap();
	int[] blockSize =afpChain.getBlockSize();
	int[] afpChainList = afpChain.getAfpChainList();
	List<AFP>afpSet = afpChain.getAfpSet();
	int[] block2Afp = afpChain.getBlock2Afp();

	double torsionPenalty = params.getTorsionPenalty();


	bkold = 0;
	for(i = 0; i < blockNum; i ++)  {
		blockScore[i] = 0;
		blockGap[i] = 0;
		for(j = 0; j < blockSize[i]; j ++)      {
			bknow = afpChainList[block2Afp[i] + j];
			if(j == 0)      {
				blockScore[i] = afpSet.get(bknow).getScore();
			}
			else    {
				AFPChainer.afpPairConn(bkold, bknow, params, afpChain); //note: j, i
				Double conn = afpChain.getConn();
				blockScore[i] += afpSet.get(bknow).getScore() + conn;
				g1 = afpSet.get(bknow).getP1() - afpSet.get(bkold).getP1() - afpSet.get(bkold).getFragLen();
				g2 = afpSet.get(bknow).getP2() - afpSet.get(bkold).getP2() - afpSet.get(bkold).getFragLen();
				blockGap[i] += (g1 > g2)?g1:g2;
			}
			bkold = bknow;
		}
		alignScoreUpdate += blockScore[i];
	}
	if(blockNum >= 2)       {
		alignScoreUpdate += (blockNum - 1) * torsionPenalty;
	}

	afpChain.setBlockGap(blockGap);
	afpChain.setAlignScoreUpdate(alignScoreUpdate);
	afpChain.setBlockScore(blockScore);
	afpChain.setBlockSize(blockSize);
	afpChain.setAfpChainList(afpChainList);
	afpChain.setBlock2Afp(block2Afp);
}
 

/**
 * in some special cases, there is no maginificent twists in the
final chaining result; however, their rmsd (original and after
optimizing) are very large. Therefore, a post-process is taken
to split the blocks further at the ralative bad connections (
with relative high distance variation)
to be tested:
  split or not according to optimized or initial chaining???
 */

private static void splitBlock(FatCatParameters params, AFPChain afpChain, Atom[] ca1, Atom[] ca2)
{
	if ( debug)
		System.err.println("AFPPostProcessor: splitBlock");
	int     i, a, bk, cut;
	double  maxs, maxt;
	int blockNum = afpChain.getBlockNum();
	int maxTra = params.getMaxTra();
	double badRmsd = params.getBadRmsd();

	int     blockNum0 = blockNum;

	double[] blockRmsd = afpChain.getBlockRmsd();
	int[] blockSize = afpChain.getBlockSize();
	int[] block2Afp = afpChain.getBlock2Afp();
	double[] afpChainTwiList = afpChain.getAfpChainTwiList();

	bk = 0;
	while(blockNum < maxTra + 1)    {
		maxs = 0;
		for(i = 0; i < blockNum; i ++)   {
			if(blockRmsd[i] > maxs && blockSize[i] > 2) { //according to the optimized alignment
				maxs = blockRmsd[i];
				bk = i;
			} //!(Note: optRmsd, not blockRmsd, according to the optimized alignment
		}
		if(maxs < badRmsd)      break;
		maxt = 0;
		cut = 0;
		for(i = 1; i < blockSize[bk]; i ++)     {
			a = i + block2Afp[bk];
			if(afpChainTwiList[a] > maxt)   {
				maxt = afpChainTwiList[a];
				cut = i;

			}
		}
		if(debug)
			System.out.println(String.format("block %d original size %d rmsd %.3f maxt %.2f cut at %d\n", bk, blockSize[bk], maxs, maxt, cut));
		for(i = blockNum - 1; i > bk; i --)     {
			block2Afp[i + 1] = block2Afp[i];
			blockSize[i + 1] = blockSize[i];
			blockRmsd[i + 1] = blockRmsd[i];
		} //update block information
		block2Afp[bk + 1] = cut + block2Afp[bk];
		blockSize[bk + 1] = blockSize[bk] - cut;
		blockSize[bk] = cut;

		if(debug)
			System.out.println(String.format("  split into %d and %d sizes\n", blockSize[bk], blockSize[bk + 1]));


		int[] afpChainList = afpChain.getAfpChainList();
		//int[] subrange1    = getSubrange(afpChainList, block2Afp[bk + 1] );
		blockRmsd[bk + 1]  = AFPChainer.calAfpRmsd(blockSize[bk + 1],  afpChainList, block2Afp[bk + 1] , afpChain, ca1, ca2);

		//int[] subrange2    = getSubrange(afpChainList, block2Afp[bk] );
		blockRmsd[bk]      = AFPChainer.calAfpRmsd(blockSize[bk],      afpChainList, block2Afp[bk], afpChain, ca1, ca2);

		//split a block at the biggest position
		blockNum ++;
		afpChain.setAfpChainList(afpChainList);
	}
	if(blockNum - blockNum0 > 0)    {
		if(debug)
			System.out.println(String.format("Split %d times:\n", blockNum - blockNum0));
		for(i = 0; i < blockNum; i ++)  {
			if(debug)
				System.out.println(String.format("  block %d size %d from %d rmsd %.3f\n", i, blockSize[i], block2Afp[i], blockRmsd[i]));
		}
	}


	afpChain.setBlockNum(blockNum);
	afpChain.setBlockSize(blockSize);
	afpChain.setBlockRmsd(blockRmsd);
	afpChain.setBlock2Afp(block2Afp);


}
 

/**
 * remove the artifical small rigid-body superimpose in the middle
 clust the similar superimpositions (caused by the small flexible
 region, which is detected as a seperate rigid superimposing region by adding
 two twists before and after it(artifically!)
 one possible solution: allowing long enough loops in the chaining process,
 which however increase the calculation complexity
 */
private static void deleteBlock(FatCatParameters params, AFPChain afpChain,Atom[] ca1, Atom[] ca2)
{
	int blockNum = afpChain.getBlockNum();
	List<AFP> afpSet = afpChain.getAfpSet();

	int[] afpChainList =afpChain.getAfpChainList();



	int[] block2Afp = afpChain.getBlock2Afp();
	int[] blockSize = afpChain.getBlockSize();

	double[] blockRmsd = afpChain.getBlockRmsd();

	int fragLen = params.getFragLen();

	//remove those blocks (both in terminals and in the middle) with only a AFP
	//but still keep those small blocks spaning large regions
	if(blockNum <= 1)       return;
	int     blockNumOld = blockNum;
	int     i, j, b1, b2, e1, e2, len;
	e1 = e2 = 0;
	for(i = 0; i < blockNum; i ++) {
		b1 = e1;
		b2 = e2;
		if(i < blockNum - 1)    {
			e1 = afpSet.get(afpChainList[block2Afp[i + 1]]).getP1();
			e2 = afpSet.get(afpChainList[block2Afp[i + 1]]).getP2();
		}
		else    {
			e1 = ca1.length;
			e2 = ca2.length;
		}
		if(blockSize[i] > 1)    continue;
		len = (e1 - b1) < (e2 - b2)?(e1 - b1):(e2 - b2);
		//if(i == blockNum - 1) blockNum --;
		if(len < 2 * fragLen)   {
			for(j = i; j < blockNum - 1; j ++)      {
				blockRmsd[j] = blockRmsd[j + 1];
				blockSize[j] = blockSize[j + 1];
				block2Afp[j] = block2Afp[j + 1];
			}
			blockNum --;
			i --;
		} //delete a block
	}
	if(blockNumOld > blockNum)
		if(debug)
			System.out.println(
					String.format("Delete %d small blocks\n", blockNumOld - blockNum)
			);


	if (debug)
		System.err.println("deleteBlock: end blockNum:"+ blockNum);
	afpChain.setBlock2Afp(block2Afp);
	afpChain.setBlockSize(blockSize);
	afpChain.setAfpChainList(afpChainList);
	afpChain.setBlockNum(blockNum);
	afpChain.setBlockRmsd(blockRmsd);
}
 

/**
 * Merge consecutive blocks with similar transformation
 */
private static  void mergeBlock(FatCatParameters params, AFPChain afpChain,Atom[] ca1,Atom[] ca2)
{

	int blockNum = afpChain.getBlockNum();
	double badRmsd = params.getBadRmsd();

	int[] block2Afp = afpChain.getBlock2Afp();
	int[] blockSize = afpChain.getBlockSize();

	double[] blockRmsd = afpChain.getBlockRmsd();

	int afpChainTwiNum = afpChain.getAfpChainTwiNum();

	//clustering the neighbor blocks if their transformations are similar
	int     i, j, b1, b2, minb1, minb2;
	double  minrmsd;
	int     merge = 0;
	int     blockNumOld = blockNum;
	double[][]  rmsdlist = null;
	if(blockNum > 1)        {
		rmsdlist = new double[blockNumOld][blockNumOld];
		for(b1 = 0; b1 < blockNum - 1; b1 ++)   {
			for(b2 = b1 + 1; b2 < blockNum; b2 ++)  {
				rmsdlist[b1][b2] = combineRmsd(b1, b2,afpChain,ca1,ca2);
			}
		}
	}
	minb1 = 0;
	while(blockNum > 1)     {
		minrmsd = 1000;
		for(i = 0; i < blockNum - 1; i ++)      {
			j = i + 1; //only consider neighbor blocks
			if(minrmsd > rmsdlist[i][j])    {
				minrmsd = rmsdlist[i][j];
				minb1 = i;
			}
		}
		minb2 = minb1 + 1; //merge those most similar blocks
		//maxrmsd = (blockRmsd[minb1] > blockRmsd[minb2])?blockRmsd[minb1]:blockRmsd[minb2];
		if(minrmsd < badRmsd)   {
			if(debug)
				System.out.println(String.format("merge block %d (rmsd %.3f) and %d (rmsd %.3f), total rmsd %.3f\n",
						minb1, blockRmsd[minb1], minb2, blockRmsd[minb2], minrmsd));
			blockSize[minb1] += blockSize[minb2];
			blockRmsd[minb1] = minrmsd;
			for(i = minb2; i < blockNum - 1; i ++)  {
				block2Afp[i] = block2Afp[i + 1];
				blockSize[i] = blockSize[i + 1];
				blockRmsd[i] = blockRmsd[i + 1];
			} //update block information
			afpChainTwiNum --;
			blockNum --;
			for(b1 = 0; b1 < blockNum - 1; b1 ++)   {
				for(b2 = b1 + 1; b2 < blockNum; b2 ++) {
					if(b1 == minb1 || b2 == minb1)  {
						rmsdlist[b1][b2] = combineRmsd(b1, b2, afpChain,ca1,ca2);
					}
					else if(b2 < minb1)     continue;
					else if(b1 < minb1)     {
						rmsdlist[b1][b2] = rmsdlist[b1][b2 + 1];
					}
					else    {
						rmsdlist[b1][b2] = rmsdlist[b1 + 1][b2 + 1];
					}
				}
			} //update the rmsdlist
			merge ++;
		} //merge two blocks
		else if(minrmsd >= 100) break;
		else    {
			rmsdlist[minb1][minb2] += 100;
		} //not merge, modify the rmsdlist so that this combination is not considered in next iteration
	}

	if(merge > 0)       {
		if(debug)
			System.out.println(String.format("Merge %d blocks, remaining %d blocks\n", merge, blockNum));
	}

	if (debug){
		System.err.println("AFPPostProcessor: mergeBlock end blocknum:" + blockNum);
	}
	afpChain.setBlock2Afp(block2Afp);
	afpChain.setBlockSize(blockSize);
	afpChain.setBlockNum(blockNum);
	afpChain.setBlockRmsd(blockRmsd);
	afpChain.setAfpChainTwiNum(afpChainTwiNum);
}
 

/**
 *
 * @param alignment The alignment to plot
 * @param background [Optional]A matrix of 'background colors' over which to draw the alignment.
 *
 *	Originally designed as a matrix of RMSD values between AFPs, so it is colorized
 *	accordingly from red (0) to black (>10).
 *
 *  If this set to null, the background is set to black.
 */
public DotPlotPanel(AFPChain alignment ){
	super();

	final double defaultBackground = 100.;

	// Convert the AFPChain alignment into the MatrixPanel format
	AlternativeAlignment[] aligns = new AlternativeAlignment[alignment.getBlockNum()];
	int alignNumber = 0;

	//One alternative alignment for each block
	int[][][] optAln = alignment.getOptAln(); // [block #][{0,1} chain index][pos]

	for(;alignNumber < optAln.length;alignNumber++) {
		List<int[]> alignPairs = new ArrayList<int[]>();
		for(int pos = 0; pos<optAln[alignNumber][0].length; pos++ ) {
			alignPairs.add( new int[] {
					optAln[alignNumber][0][pos],
					optAln[alignNumber][1][pos] }
			);
		}
		JointFragments frag = new JointFragments();
		frag.setIdxlist(alignPairs);
		aligns[alignNumber] = new AlternativeAlignment();
		aligns[alignNumber].apairs_from_idxlst(frag);

	}

	/* TODO After the AFPSet is fixed in CeMain#filterDuplicateAFPs, maybe include this again
	//add alignment for the AFPs
	List<AFP> afps = alignment.getAfpSet();
	List<int[]> alignPairs = new ArrayList<int[]>();
	for(AFP afp : afps) {
		int start1 = afp.getP1();
		int start2 = afp.getP2();
		for(int i=0;i<afp.getFragLen();i++) {
			alignPairs.add( new int[] { start1+i, start2+i } );
		}
	}
	JointFragments frag = new JointFragments();
	frag.setIdxlist(alignPairs);
	aligns[alignNumber] = new AlternativeAlignment();
	aligns[alignNumber].apairs_from_idxlst(frag);
	*/


	/* AFP boxes are unnecessary.
	// Calculate FragmentPairs based on alignment.
	// These are displayed as a small box around the start of each alignment.
	FragmentPair[] pairs = new FragmentPair[afps.size()];
	for(int i=0;i<pairs.length;i++) {
		AFP afp = afps.get(i);
		pairs[i] = new FragmentPair(afp.getFragLen(), afp.getP1(), afp.getP2());
		pairs[i].setRms(afp.getRmsd());
	}

	this.setFragmentPairs(pairs);
	*/


	// Now the alignments have been build; add it
	this.setAlternativeAligs(aligns);
	this.setSelectedAlignmentPos(0); //color white, not red

	Matrix background = alignment.getDistanceMatrix();
	//Fill with default black background if none given
	if(background == null) {
		background = new Matrix(alignment.getCa1Length(),alignment.getCa2Length());
		for(int i=0;i<background.getRowDimension();i++)
			for(int j=0;j<background.getColumnDimension(); j++) {
				background.set(i, j, defaultBackground);
			}
	}

	// Set parameters
	this.setMatrix(background);
}
 

private static final int getUngappedFatCatPos(AFPChain afpChain, int chainNr, int aligPos){
	char[] aseq;
	if ( chainNr == 0 )
		aseq = afpChain.getAlnseq1();
	else
		aseq = afpChain.getAlnseq2();

	if ( aligPos > aseq.length)
		return -1;
	if ( aligPos < 0)
		return -1;

	int blockNum = afpChain.getBlockNum();
	int[] optLen = afpChain.getOptLen();
	int[][][] optAln = afpChain.getOptAln();

	int p1, p2;
	int p1b = 0;
	int p2b = 0;
	int len = 0;


	for(int i = 0; i < blockNum; i ++)  {
		for(int j = 0; j < optLen[i]; j ++) {

			p1 = optAln[i][0][j];
			p2 = optAln[i][1][j];


			if(len > 0)     {

				int lmax = (p1 - p1b - 1)>(p2 - p2b - 1)?(p1 - p1b - 1):(p2 - p2b - 1);

				// lmax gives the length of an alignment gap

				//System.out.println("   p1-p2: " + p1 + " - " + p2 + " lmax: " + lmax + " p1b-p2b:"+p1b + " " + p2b);
				for(int k = 0; k < lmax; k ++)      {

					if(k >= (p1 - p1b - 1)) {
						// a gap position in chain 0
						if ( aligPos == len && chainNr == 0){
							return -1;
						}
					}
					else {
						if ( aligPos == len && chainNr == 0)
							return p1b+1+k;


					}
					if(k >= (p2 - p2b - 1)) {
						// a gap position in chain 1
						if ( aligPos == len && chainNr == 1){
							return -1;
						}
					}
					else  {
						if ( aligPos == len && chainNr == 1) {
							return p2b+1+k;
						}


					}
					len++;

				}
			}

			if ( aligPos == len && chainNr == 0)
				return p1;
			if ( aligPos == len && chainNr == 1)
				return p2;

			len++;
			p1b = p1;
			p2b = p2;


		}
	}


	// we did not find an aligned position
	return -1;
}
 

private static String getMultiBlockJmolScript(AFPChain afpChain, Atom[] ca1, Atom[] ca2) {

		int blockNum = afpChain.getBlockNum();
		int[] optLen = afpChain.getOptLen();
		int[][][] optAln = afpChain.getOptAln();

		if (optLen == null)
			return DEFAULT_SCRIPT;

		StringWriter jmol = new StringWriter();
		jmol.append(DEFAULT_SCRIPT);

		jmol.append("select */2; color lightgrey; model 2; ");

		for (int bk = 0; bk < blockNum; bk++) {

			printJmolScript4Block(ca1, ca2, blockNum, optLen, optAln, jmol, bk);
		}
		jmol.append("model 0;  ");

		jmol.append(LIGAND_DISPLAY_SCRIPT);
		// System.out.println(jmol);
		return jmol.toString();

	}
 

/**
 * Fundamentally, an alignment is just a list of aligned residues in each
 * protein. This method converts two lists of ResidueNumbers into an
 * AFPChain.
 *
 * <p>Parameters are filled with defaults (often null) or sometimes
 * calculated.
 *
 * <p>For a way to modify the alignment of an existing AFPChain, see
 * {@link AlignmentTools#replaceOptAln(AFPChain, Atom[], Atom[], Map)}
 * @param ca1 CA atoms of the first protein
 * @param ca2 CA atoms of the second protein
 * @param aligned1 A list of aligned residues from the first protein
 * @param aligned2 A list of aligned residues from the second protein.
 *  Must be the same length as aligned1.
 * @return An AFPChain representing the alignment. Many properties may be
 *  null or another default.
 * @throws StructureException if an error occured during superposition
 * @throws IllegalArgumentException if aligned1 and aligned2 have different
 *  lengths
 * @see AlignmentTools#replaceOptAln(AFPChain, Atom[], Atom[], Map)
 */
public static AFPChain createAFPChain(Atom[] ca1, Atom[] ca2,
									  ResidueNumber[] aligned1, ResidueNumber[] aligned2 ) throws StructureException {
	//input validation
	int alnLen = aligned1.length;
	if(alnLen != aligned2.length) {
		throw new IllegalArgumentException("Alignment lengths are not equal");
	}

	AFPChain a = new AFPChain(AFPChain.UNKNOWN_ALGORITHM);
	try {
		a.setName1(ca1[0].getGroup().getChain().getStructure().getName());
		if(ca2[0].getGroup().getChain().getStructure() != null) {
			// common case for cloned ca2
			a.setName2(ca2[0].getGroup().getChain().getStructure().getName());
		}
	} catch(Exception e) {
		// One of the structures wasn't fully created. Ignore
	}
	a.setBlockNum(1);
	a.setCa1Length(ca1.length);
	a.setCa2Length(ca2.length);

	a.setOptLength(alnLen);
	a.setOptLen(new int[] {alnLen});


	Matrix[] ms = new Matrix[a.getBlockNum()];
	a.setBlockRotationMatrix(ms);
	Atom[] blockShiftVector = new Atom[a.getBlockNum()];
	a.setBlockShiftVector(blockShiftVector);

	String[][][] pdbAln = new String[1][2][alnLen];
	for(int i=0;i<alnLen;i++) {
		pdbAln[0][0][i] = aligned1[i].getChainName()+":"+aligned1[i];
		pdbAln[0][1][i] = aligned2[i].getChainName()+":"+aligned2[i];
	}

	a.setPdbAln(pdbAln);

	// convert pdbAln to optAln, and fill in some other basic parameters
	AFPChainXMLParser.rebuildAFPChain(a, ca1, ca2);

	return a;

	// Currently a single block. Split into several blocks by sequence if needed
	//		return AlignmentTools.splitBlocksByTopology(a,ca1,ca2);
}
 

@Test
public void testCECP1() throws IOException, StructureException{

	String name1 = "PDP:3A2KAc";
	String name2 = "d1wy5a2";


	CeCPMain algorithm = new CeCPMain();

	AtomCache cache = new AtomCache();

	Atom[] ca1 = cache.getAtoms(name1);
	Atom[] ca2 = cache.getAtoms(name2);

	AFPChain afpChain = algorithm.align(ca1, ca2);
	CECalculator calculator = algorithm.getCECalculator();

	//               System.out.println(calculator.get);
	//StructureAlignmentJmol jmol =
	//StructureAlignmentDisplay.display(afpChain, ca1, ca2);
	if ( ! (afpChain.getBlockNum() == 1)){
		System.out.println(calculator.getLcmp());
		System.out.println(afpChain.toFatcat(ca1, ca2));
	}
	Assert.assertEquals(1, afpChain.getBlockNum());


}
 

private void testAlignment(String name1, String name2, Atom[] ca1, Atom[] ca2, boolean doRigid) throws StructureException,IOException{


		Atom[] ca3 = StructureTools.cloneAtomArray(ca2);


		AFPChain afpChain = doAlign(name1, name2, ca1,ca2,doRigid);



		String fatcat = afpChain.toFatcat(ca1, ca2);
		String xml 	  = AFPChainXMLConverter.toXML(afpChain,ca1,ca2);


		//System.out.println(xml);
		AFPChain newChain = AFPChainXMLParser.fromXML (xml, ca1, ca3);

		// test blockNum and optLen arrays
		int blockNum = afpChain.getBlockNum();
		int[] optLen = afpChain.getOptLen();

		assertTrue("The nr of aligned blocks is not the same! " + blockNum + " " + newChain.getBlockNum() , blockNum == newChain.getBlockNum());



		for ( int i =0 ; i < blockNum ; i++){
			int newLenI = newChain.getOptLen()[i];
			assertTrue("The values in the optLen field don't match! pos:" + i + " orig:" + optLen[i] + " new:" +  newLenI,optLen[i] == newLenI);
		}

		// test the internal optAlign data structure:

		int[][][] optAln1 = afpChain.getOptAln();
		int[][][] optAln2 = newChain.getOptAln();

		for(int i = 0; i < blockNum; i ++)  {
			for(int j = 0; j < optLen[i]; j ++) {
				int p1 = optAln1[i][0][j];
				int p2 = optAln1[i][1][j];

				int n1 = optAln2[i][0][j];
				int n2 = optAln2[i][1][j];

				assertTrue(p1 == n1);
				assertTrue(p2 == n2);

			}
		}

		String fatcat2 = newChain.toFatcat(ca1, ca3);
		//System.out.println("*** RESULT2 "+result2);

		assertEquals(fatcat,fatcat2);

		String xml2 = AFPChainXMLConverter.toXML(newChain, ca1, ca3);
		assertEquals(xml, xml2);

	}