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

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);
}
 

/**
 * Delete an alignment position from the original alignment object.
 *
 * @param afpChain
 *            original alignment, will be modified
 * @param ca1
 *            atom array, will not be modified
 * @param ca2
 *            atom array, will not be modified
 * @param block
 *            block of the alignment position
 * @param pos
 *            position index in the block
 * @return the original alignment, with the alignment position removed
 * @throws StructureException
 */
public static AFPChain deleteColumn(AFPChain afpChain, Atom[] ca1,
		Atom[] ca2, int block, int pos) throws StructureException {

	// Check validity of the inputs
	if (afpChain.getBlockNum() <= block) {
		throw new IndexOutOfBoundsException(String.format(
				"Block index requested (%d) is higher than the total number of AFPChain blocks (%d).",
				block, afpChain.getBlockNum()));
	}
	if (afpChain.getOptAln()[block][0].length <= pos) {
		throw new IndexOutOfBoundsException(String.format(
				"Position index requested (%d) is higher than the total number of aligned position in the AFPChain block (%d).",
				block, afpChain.getBlockSize()[block]));
	}

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

	int[] newPos0 = new int[optAln[block][0].length - 1];
	int[] newPos1 = new int[optAln[block][1].length - 1];

	int position = 0;
	for (int p = 0; p < optAln[block][0].length; p++) {

		if (p == pos)
			continue;

		newPos0[position] = optAln[block][0][p];
		newPos1[position] = optAln[block][1][p];

		position++;
	}

	optAln[block][0] = newPos0;
	optAln[block][1] = newPos1;

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

/**
 * Find the alignment position with the highest atomic distance between the
 * equivalent atomic positions of the arrays and remove it from the
 * alignment.
 *
 * @param afpChain
 *            original alignment, will be modified
 * @param ca1
 *            atom array, will not be modified
 * @param ca2
 *            atom array, will not be modified
 * @return the original alignment, with the alignment position at the
 *         highest distance removed
 * @throws StructureException
 */
public static AFPChain deleteHighestDistanceColumn(AFPChain afpChain,
		Atom[] ca1, Atom[] ca2) throws StructureException {

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

	int maxBlock = 0;
	int maxPos = 0;
	double maxDistance = Double.MIN_VALUE;

	for (int b = 0; b < optAln.length; b++) {
		for (int p = 0; p < optAln[b][0].length; p++) {
			Atom ca2clone = ca2[optAln[b][1][p]];
			Calc.rotate(ca2clone, afpChain.getBlockRotationMatrix()[b]);
			Calc.shift(ca2clone, afpChain.getBlockShiftVector()[b]);

			double distance = Calc.getDistance(ca1[optAln[b][0][p]],
					ca2clone);
			if (distance > maxDistance) {
				maxBlock = b;
				maxPos = p;
				maxDistance = distance;
			}
		}
	}

	return deleteColumn(afpChain, ca1, ca2, maxBlock, maxPos);
}
 

/**
 * 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;
}
 

public static void printXMLEQRInferPositions(PrettyXMLWriter xml,
		AFPChain afpChain, int bk, Atom[] ca1, Atom[] ca2)  throws IOException{

	int[] optLen       = afpChain.getOptLen();

	if ( optLen == null)
		return;

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


	for ( int pos=0;pos< optLen[bk];pos++){
		int pos1 = optAln[bk][0][pos];
		int pos2 = optAln[bk][1][pos];
		xml.openTag("eqr");
		xml.attribute("eqrNr", String.valueOf(pos));
		xml.attribute("pdbres1",ca1[pos1].getGroup().getResidueNumber().toString());
		xml.attribute("chain1", ca1[pos1].getGroup().getChain().getName());
		xml.attribute("pdbres2",ca2[pos2].getGroup().getResidueNumber().toString());
		xml.attribute("chain2", ca2[pos2].getGroup().getChain().getName());

		xml.closeTag("eqr");
		//System.out.println("aligned position: " + pos1  + ":" + pos2 +
		//" pdbresnum " + ca1[pos1].getGroup().getResidueNumber().toString() + " " +
		//ca1[pos1].getParent().getPDBName()+":" +
		//ca2[pos2].getGroup().getResidueNumber().toString() + " " + ca2[pos2].getParent().getPDBName());

	}

}
 

/**
 * Print an AFPChain manually for debugging
 * @param cpAlignment
 */
@SuppressWarnings("unused")
private static void printOptAln(AFPChain cpAlignment) {
	int[] optLen = cpAlignment.getOptLen();
	int[][][] optAln = cpAlignment.getOptAln();

	for(int block=0;block<cpAlignment.getBlockNum();block++) {
		for(int pos=0;pos<optLen[block]; pos++) {
			System.out.format("%s\t%s\n", optAln[block][0][pos],
					optAln[block][1][pos]);
		}
		System.out.println();
	}
}
 

/**
 * 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);
}
 

public static Atom[] getAlignedAtoms1(AFPChain afpChain,Atom[] ca1){
	List<Atom> atoms = new ArrayList<Atom>();

	int blockNum = afpChain.getBlockNum();

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


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


		for ( int i=0;i< optLen[bk];i++){
			int pos = optAln[bk][0][i];
			atoms.add(ca1[pos]);
		}

	}
	return atoms.toArray(new Atom[atoms.size()]);
}
 

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();

	}
 

/** get the block number for an aligned position
 *
 * @param afpChain
 * @param aligPos
 * @return
 */
public static int getBlockNrForAlignPos(AFPChain afpChain, int aligPos){
	// moved here from DisplayAFP;

	int blockNum = afpChain.getBlockNum();

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

	int len = 0;
	int p1b=0;
	int p2b=0;

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

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

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

			if (len != 0) {
				// check for gapped region
				int lmax = (p1 - p1b - 1)>(p2 - p2b - 1)?(p1 - p1b - 1):(p2 - p2b - 1);
				for(int k = 0; k < lmax; k ++)      {
					len++;
				}
			}


			p1b = p1;
			p2b = p2;
			if ( len >= aligPos) {

				return i;
			}
			len++;
		}
	}

	return blockNum;

}
 

@Test
public void testIsSequential() throws StructureException, IOException {
	AtomCache cache = new AtomCache();

	String name1, name2;
	Atom[] ca1, ca2;
	AFPChain afpChain;
	StructureAlignment ce;


	// CP case
	name1="1QDM.A"; // swaposin
	name2="1NKL"; // saposin

	ca1=cache.getAtoms(name1);
	ca2=cache.getAtoms(name2);

	ce = StructureAlignmentFactory.getAlgorithm(CeCPMain.algorithmName);
	afpChain = ce.align(ca1,ca2);

	Assert.assertFalse("CeCPMain should give non-sequential alignments (between blocks).", AlignmentTools.isSequentialAlignment(afpChain, false));
	Assert.assertFalse("CeCPMain should give non-sequential alignments (within blocks).", AlignmentTools.isSequentialAlignment(afpChain, true));

	// linear case
	ce = StructureAlignmentFactory.getAlgorithm(CeMain.algorithmName);
	afpChain = ce.align(ca1,ca2);

	Assert.assertTrue("CeMain should give sequential alignments (within blocks).", AlignmentTools.isSequentialAlignment(afpChain, true));
	Assert.assertTrue("CeMain should give sequential alignments (between blocks).", AlignmentTools.isSequentialAlignment(afpChain, false));

	// now change the block interior a bit

	int[][][] optAln = afpChain.getOptAln();
	int tmp;
	tmp = optAln[0][0][2];
	optAln[0][0][2] = optAln[0][0][1];
	optAln[0][0][1] = tmp;
	tmp = optAln[0][1][2];
	optAln[0][1][2] = optAln[0][1][1];
	optAln[0][1][1] = tmp;

	Assert.assertTrue("Modifying block interior shouldn't effect block sequence.", AlignmentTools.isSequentialAlignment(afpChain, false));
	Assert.assertFalse("Modifying block interior should be not sequential.", AlignmentTools.isSequentialAlignment(afpChain, true));

}
 

/**
 * Swaps the order of structures in an AFPChain
 * @param a
 * @return
 */
public AFPChain invertAlignment(AFPChain a) {
	String name1 = a.getName1();
	String name2 = a.getName2();
	a.setName1(name2);
	a.setName2(name1);

	int len1 = a.getCa1Length();
	a.setCa1Length( a.getCa2Length() );
	a.setCa2Length( len1 );

	int beg1 = a.getAlnbeg1();
	a.setAlnbeg1(a.getAlnbeg2());
	a.setAlnbeg2(beg1);

	char[] alnseq1 = a.getAlnseq1();
	a.setAlnseq1(a.getAlnseq2());
	a.setAlnseq2(alnseq1);

	Matrix distab1 = a.getDisTable1();
	a.setDisTable1(a.getDisTable2());
	a.setDisTable2(distab1);

	int[] focusRes1 = a.getFocusRes1();
	a.setFocusRes1(a.getFocusRes2());
	a.setFocusRes2(focusRes1);

	//What are aftIndex and befIndex used for? How are they indexed?
	//a.getAfpAftIndex()


	String[][][] pdbAln = a.getPdbAln();
	if( pdbAln != null) {
		for(int block = 0; block < a.getBlockNum(); block++) {
			String[] paln1 = pdbAln[block][0];
			pdbAln[block][0] = pdbAln[block][1];
			pdbAln[block][1] = paln1;
		}
	}

	int[][][] optAln = a.getOptAln();
	if( optAln != null ) {
		for(int block = 0; block < a.getBlockNum(); block++) {
			int[] aln1 = optAln[block][0];
			optAln[block][0] = optAln[block][1];
			optAln[block][1] = aln1;
		}
	}
	a.setOptAln(optAln); // triggers invalidate()

	Matrix distmat = a.getDistanceMatrix();
	if(distmat != null)
		a.setDistanceMatrix(distmat.transpose());


	// invert the rotation matrices
	Matrix[] blockRotMat = a.getBlockRotationMatrix();
	Atom[] shiftVec = a.getBlockShiftVector();
	if( blockRotMat != null) {
		for(int block = 0; block < a.getBlockNum(); block++) {
			if(blockRotMat[block] != null) {
				// if y=x*A+b, then x=y*inv(A)-b*inv(A)
				blockRotMat[block] = blockRotMat[block].inverse();

				Calc.rotate(shiftVec[block],blockRotMat[block]);
				shiftVec[block] = Calc.invert(shiftVec[block]);
			}
		}
	}

	return a;
}
 

public  static double getTMScore(AFPChain align, Atom[] ca1, Atom[] ca2, boolean normalizeMin) throws StructureException
{
	if ( align.getNrEQR() == 0)
		return -1;


	// Create new arrays for the subset of atoms in the alignment.
	Atom[] ca1aligned = new Atom[align.getOptLength()];
	Atom[] ca2aligned = new Atom[align.getOptLength()];
	int pos=0;
	int[] blockLens = align.getOptLen();
	int[][][] optAln = align.getOptAln();
	assert(align.getBlockNum() <= optAln.length);

	for(int block=0;block< align.getBlockNum();block++) {

		if ( ! ( blockLens[block] <= optAln[block][0].length)) {
			logger.warn("AFPChainScorer getTMScore: errors reconstructing alignment block [" + block + "]. Length is " + blockLens[block] + " but should be <=" + optAln[block][0].length);
		}

		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 != align.getOptLength()){
		logger.warn("AFPChainScorer getTMScore: Problems reconstructing alignment! nr of loaded atoms is " + pos + " but should be " + align.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);
	}
	//Superimpose
	Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(ca1aligned),
			Calc.atomsToPoints(ca2aligned));

	Calc.transform(ca2aligned, trans);

	return Calc.getTMScore(ca1aligned, ca2aligned, ca1.length, ca2.length, normalizeMin);
}
 

/**
 * superimposing according to the optimized alignment
 *
 * @param afpChain
 * @param ca1
 * @param ca2
 * @return Group array twisted.
 * @throws StructureException
 */
public static Group[] twistOptimized(AFPChain afpChain, Atom[] ca1,
		Atom[] ca2) throws StructureException {

	Atom[] optTwistPdb = new Atom[ca2.length];

	int gPos = -1;
	for (Atom a : ca2) {
		gPos++;
		optTwistPdb[gPos] = a;
	}

	int blockNum = afpChain.getBlockNum();

	int b2 = 0;
	int e2 = 0;
	int focusResn = 0;
	int[] focusRes1 = afpChain.getFocusRes1();
	int[] focusRes2 = afpChain.getFocusRes2();

	if (focusRes1 == null) {
		focusRes1 = new int[afpChain.getCa1Length()];
		afpChain.setFocusRes1(focusRes1);
	}
	if (focusRes2 == null) {
		focusRes2 = new int[afpChain.getCa2Length()];
		afpChain.setFocusRes2(focusRes2);
	}

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

	for (int bk = 0; bk < blockNum; bk++) {
		// THIS IS TRANSFORMING THE ORIGINAL ca2 COORDINATES, NO CLONING...
		// copies the atoms over to iniTwistPdb later on in modifyCod
		transformOrigPDB(optLen[bk], optAln[bk][0], optAln[bk][1], ca1,
				ca2, afpChain, bk);

		// transform pro2 according to comparison of pro1 and pro2 at give
		// residues
		if (bk > 0) {
			b2 = e2;
		}
		if (bk < blockNum - 1) { // bend at the middle of two consecutive
									// blocks
			e2 = optAln[bk][1][optLen[bk] - 1];
			e2 = (optAln[bk + 1][1][0] - e2) / 2 + e2;
		} else {
			e2 = ca2.length;
		}
		cloneAtomRange(optTwistPdb, ca2, b2, e2);
		for (int i = 0; i < optLen[bk]; i++) {
			focusRes1[focusResn] = optAln[bk][0][i];
			focusRes2[focusResn] = optAln[bk][1][i];
			focusResn++;
		}
	}
	int totalLenOpt = focusResn;
	logger.debug("calrmsdopt for {} residues", focusResn);
	double totalRmsdOpt = calCaRmsd(ca1, optTwistPdb, focusResn, focusRes1,
			focusRes2);
	logger.debug("got opt RMSD: {}", totalRmsdOpt);
	int optLength = afpChain.getOptLength();

	if (totalLenOpt != optLength) {
		logger.warn("Final alignment length is different {} {}",
				totalLenOpt, optLength);
	}
	logger.debug("final alignment length {}, rmsd {}", focusResn,
			totalRmsdOpt);

	afpChain.setTotalLenOpt(totalLenOpt);
	afpChain.setTotalRmsdOpt(totalRmsdOpt);

	return StructureTools.cloneGroups(optTwistPdb);

}
 

@Test
public void testAFPconversion() throws Exception{

	//Fill an AFPChain with the general information
	AFPChain afp = new AFPChain("algorithm");
	afp.setName1("name1");
	afp.setName2("name2");
	afp.setVersion("1.0");
	afp.setCalculationTime(System.currentTimeMillis());
	//Generate a optimal alignment with 3 blocks and 5 residues per block
	int[][][] optAln = new int[3][][];
	for (int b=0; b<optAln.length; b++){
		int[][] block = new int[2][];
		for (int c=0; c<block.length; c++){
			int[] residues = {b+5,b+6,b+7,b+8,b+9};
			block[c] = residues;
		}
		optAln[b] = block;
	}
	afp.setOptAln(optAln);
	afp.setBlockNum(optAln.length);
	//Set the rotation matrix and shift to random numbers
	double[][] mat = {{0.13,1.5,0.84},{1.3,0.44,2.3},{1.0,1.2,2.03}};
	Matrix rot = new Matrix(mat);
	Atom shift = new AtomImpl();
	shift.setX(0.44);
	shift.setY(0.21);
	shift.setZ(0.89);
	Matrix[] blockRot = {rot,rot,rot};
	afp.setBlockRotationMatrix(blockRot);
	Atom[] blockShift = {shift,shift,shift};
	afp.setBlockShiftVector(blockShift);

	//Convert the AFPChain into a MultipleAlignment (without Atoms)
	MultipleAlignmentEnsemble ensemble =
			new MultipleAlignmentEnsembleImpl(afp,null,null,true);
	MultipleAlignment msa = ensemble.getMultipleAlignment(0);

	//Test for all the information
	assertEquals(afp.getName1(),ensemble.getStructureIdentifiers().get(0).getIdentifier());
	assertEquals(afp.getName2(), ensemble.getStructureIdentifiers().get(1).getIdentifier());
	assertEquals(afp.getAlgorithmName(), ensemble.getAlgorithmName());
	assertEquals(afp.getVersion(),ensemble.getVersion());
	assertTrue(ensemble.getCalculationTime().equals(
			afp.getCalculationTime()));
	assertEquals(afp.getBlockNum(), msa.getBlockSets().size());
	for (int b = 0; b<afp.getBlockNum(); b++){
		assertEquals(Calc.getTransformation(
				afp.getBlockRotationMatrix()[b],
				afp.getBlockShiftVector()[b]),
				msa.getBlockSet(b).getTransformations().get(1));
	}

	//Test for the scores
	assertEquals(msa.getScore(MultipleAlignmentScorer.CE_SCORE),
			(Double) afp.getAlignScore());
	assertEquals(msa.getScore(MultipleAlignmentScorer.AVGTM_SCORE),
			(Double) afp.getTMScore());
	assertEquals(msa.getScore(MultipleAlignmentScorer.RMSD),
			(Double) afp.getTotalRmsdOpt());


	//Test for the optimal alignment
	for (int b=0; b<3; b++){
		for (int c=0; c<2; c++){
			for (int res=0; res<5; res++){
				Integer afpRes = afp.getOptAln()[b][c][res];
				assertEquals(afpRes, msa.getBlock(b).
						getAlignRes().get(c).get(res));
			}
		}
	}
}
 

@Test
public void testDuplicateIndices() throws IOException, StructureException {
	String[] algorithmIDs = {CeMain.algorithmName, FatCatRigid.algorithmName};

	AtomCache cache = new AtomCache();

	// 3j47 is a bunch of a-helices, so there are many valid ways to align chains
	// structurally between each other.
	List<Chain> chains = cache.getStructure("3j47").getChains();

	for(String algorithmID:algorithmIDs) {
		StructureAlignment algorithm = StructureAlignmentFactory.getAlgorithm(algorithmID);
		System.out.println("Testing "+algorithmID);
		for (int c1 = 0; c1<chains.size()-1;c1++) {
			for (int c2=chains.size()-1;c2>c1;c2--) {
				Atom[] ca1 = StructureTools.getAtomCAArray(chains.get(c1));
				Atom[] ca2 = StructureTools.getAtomCAArray(chains.get(c2));

				AFPChain afpChain_fc = algorithm.align(ca1, ca2);
				assertEquals(1,afpChain_fc.getOptAln().length);

				int[][] optAln = afpChain_fc.getOptAln()[0];
				// two chains aligned
				assertEquals(2,optAln.length);

				//same number of aligned residues between the chains
				assertEquals(optAln[0].length,optAln[1].length);

				// no indices duplicated in the alignments
				for (int[] optAlnSeq : optAln) {
					long count_unique = Arrays.stream(optAlnSeq).distinct().count();
					long count_all = optAlnSeq.length;
					assertEquals(count_unique, count_all);
				}

			}
		}

	}

}
 

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;
}
 

/**
 *
 * @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);
}
 

public static String getJmolScript4Block(AFPChain afpChain, Atom[] ca1, Atom[] ca2, int blockNr) {
	int blockNum = afpChain.getBlockNum();

	if (blockNr >= blockNum)
		return DEFAULT_SCRIPT;

	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; ");

	printJmolScript4Block(ca1, ca2, blockNum, optLen, optAln, jmol, blockNr);

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

}
 

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);

	}