/*
* BioJava development code
*
* This code may be freely distributed and modified under the
* terms of the GNU Lesser General Public Licence. This should
* be distributed with the code. If you do not have a copy,
* see:
*
* http://www.gnu.org/copyleft/lesser.html
*
* Copyright for this code is held jointly by the individual
* authors. These should be listed in @author doc comments.
*
* For more information on the BioJava project and its aims,
* or to join the biojava-l mailing list, visit the home page
* at:
*
* http://www.biojava.org/
*
* Created on Jan 28, 2006
*
*/
package org.biojava.nbio.structure.align.pairwise;
import org.biojava.nbio.structure.Atom;
import org.biojava.nbio.structure.Calc;
import org.biojava.nbio.structure.StructureException;
import org.biojava.nbio.structure.align.StrucAligParameters;
import org.biojava.nbio.structure.align.helper.AlignTools;
import org.biojava.nbio.structure.align.helper.JointFragments;
import org.biojava.nbio.structure.jama.Matrix;
import java.io.Serializable;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;
import java.util.List;
import java.util.logging.Logger;
/** Joins the initial Fragments together to larger Fragments
*
* @author Andreas Prlic
* @author Peter Lackner
* @since 1.5
* @version %I% %G%
*/
public class FragmentJoiner {
public static Logger logger = Logger.getLogger("org.biojava.nbio.structure.align");
public FragmentJoiner() {
super();
}
/**
* Reallocates an array with a new size, and copies the contents
* of the old array to the new array.
* @param oldArray the old array, to be reallocated.
* @param newSize the new array size.
* @return A new array with the same contents.
*/
@SuppressWarnings("rawtypes")
public static Object resizeArray (Object oldArray, int newSize) {
int oldSize = java.lang.reflect.Array.getLength(oldArray);
Class elementType = oldArray.getClass().getComponentType();
Object newArray = java.lang.reflect.Array.newInstance(
elementType,newSize);
int preserveLength = Math.min(oldSize,newSize);
if (preserveLength > 0)
System.arraycopy (oldArray,0,newArray,0,preserveLength);
return newArray;
}
/**
* In helices often many similar fragments can be found. To reduce these to a few
* representative ones this check can be used. It does a distance check between
* all known Fragments and a new one. If this one is on a similar diagonal and it
* has a lower rms, this one is a better representation. Note: shifts of one are
* not allowed.
*
* @param fragments
* @param f
* @param rmsmat
* @return true - if this is a better representant for a group of locala fragments.
*/
public static boolean reduceFragments(List<FragmentPair> fragments, FragmentPair f, Matrix rmsmat){
boolean doNotAdd =false;
int i = f.getPos1();
int j = f.getPos2();
for ( int p =0; p < fragments.size(); p++){
FragmentPair tmp = fragments.get(p);
int di1 = Math.abs(f.getPos1() - tmp.getPos1());
int di2 = Math.abs(f.getPos2() - tmp.getPos2());
if (( Math.abs(di1-di2) == 2)) {
double rms1 = rmsmat.get(tmp.getPos1(),tmp.getPos2());
double rms2 = rmsmat.get(i,j);
doNotAdd = true;
if ( rms2 < rms1){
fragments.remove(p);
fragments.add(f);
break;
}
p = fragments.size();
}
}
return doNotAdd;
}
public JointFragments[] approach_ap3(Atom[] ca1, Atom[] ca2, FragmentPair[] fraglst,
StrucAligParameters params) throws StructureException {
//the final list of joined fragments stores as apairs
List<JointFragments> fll = new ArrayList<JointFragments>();
FragmentPair[] tmpfidx = new FragmentPair[fraglst.length];
for ( int i=0 ; i < fraglst.length; i++){
tmpfidx[i] = (FragmentPair)fraglst[i].clone();
}
int n = tmpfidx.length;
// if two fragments after having been joint have rms < X
// keep the new fragment.
for (int i=0;i< fraglst.length;i++){
boolean[] used = new boolean[n];
int p1i = tmpfidx[i].getPos1();
int p1j = tmpfidx[i].getPos2();
int l1 = tmpfidx[i].getLength();
JointFragments f = new JointFragments();
int maxi=p1i+l1-1;
f.add(p1i,p1j,0,l1);
used[i] = true;
//n = tmpfidx.length;
for (int j=(i+1);j<n;j++){
if ( used[j])
continue;
int p2i = tmpfidx[j].getPos1();
int p2j = tmpfidx[j].getPos2();
int l2 = tmpfidx[j].getLength();
if ( p2i > maxi) {
// TODO: replace this with plo angle calculation
if ( params.isDoAngleCheck()){
// was 0.174
if (! angleCheckOk(tmpfidx[i], tmpfidx[j], 0.4f))
continue;
}
if ( params.isDoDistanceCheck()) {
if (! distanceCheckOk(tmpfidx[i],tmpfidx[j], params.getFragmentMiniDistance()))
continue;
}
if ( params.isDoDensityCheck()) {
if ( ! densityCheckOk(ca1,ca2, f.getIdxlist(), p2i,p2j, l2, params.getDensityCutoff()))
continue;
}
if ( params.isDoRMSCheck()) {
double rms = rmsCheck(ca1,ca2, f.getIdxlist(), p2i, p2j, l2);
if ( rms > params.getJoinRMSCutoff())
continue;
f.setRms(rms);
}
f.add(p2i,p2j,0,l2);
used[j] = true;
maxi = p2i+l2-1;
}
}
fll.add(f);
}
Comparator<JointFragments> comp = new JointFragmentsComparator();
Collections.sort(fll,comp);
Collections.reverse(fll);
int m = Math.min(params.getMaxrefine(),fll.size());
List<JointFragments> retlst = new ArrayList<JointFragments>();
for ( int i = 0 ; i < m ; i++){
JointFragments jf = fll.get(i);
retlst.add(jf);
}
return retlst.toArray(new JointFragments[retlst.size()]);
}
private boolean densityCheckOk(Atom[] aa1, Atom[] aa2, List<int[]> idxlist,
int p2i, int p2j, int l2,
double densityCutoff) throws StructureException {
JointFragments ftmp = new JointFragments();
ftmp.setIdxlist(idxlist);
ftmp.add(p2i,p2j,0,l2);
AlternativeAlignment ali = new AlternativeAlignment();
ali.apairs_from_idxlst(ftmp);
Atom[] aa3 = aa2.clone();
int[] idx1 = ali.getIdx1();
int[] idx2 = ali.getIdx2();
Atom[] ca1subset = AlignTools.getFragmentFromIdxList(aa1, idx1);
Atom[] ca2subset = AlignTools.getFragmentFromIdxList(aa3,idx2);
double density = getDensity(ca1subset, ca2subset);
return density <= densityCutoff;
}
/** this is probably useless
*
* @param ca1subset
* @param ca2subset
* @return a double
* @throws StructureException
*/
private double getDensity(Atom[] ca1subset, Atom[] ca2subset ) throws StructureException{
Atom centroid1 = Calc.getCentroid(ca1subset);
Atom centroid2 = Calc.getCentroid(ca2subset);
// get Average distance to centroid ...
double d1 = 0;
double d2 = 0;
for ( int i = 0 ; i < ca1subset.length;i++){
double dd1 = Calc.getDistance(centroid1, ca1subset[i]);
double dd2 = Calc.getDistance(centroid2, ca2subset[i]);
d1 += dd1;
d2 += dd2;
}
double avd1 = d1 / ca1subset.length;
double avd2 = d2 / ca2subset.length;
return Math.min(avd1,avd2);
}
private double rmsCheck(Atom[] a1, Atom[] a2,List<int[]> idxlist, int p2i, int p2j, int l2) throws StructureException {
//System.out.println(dd);
// check if a joint fragment has low rms ...
JointFragments ftmp = new JointFragments();
ftmp.setIdxlist(idxlist);
ftmp.add(p2i,p2j,0,l2);
Atom[] a3 = new Atom[a2.length];
for (int i=0;i < a2.length;i++){
a3[i] = (Atom)a2[i].clone();
}
return getRMS(a1,a3,ftmp);
}
/** get the RMS of the JointFragments pair frag
*
* @param ca1 the array of all atoms of structure1
* @param ca2 the array of all atoms of structure1
* @param frag the JointFragments object that contains the list of identical positions
* @return the rms
*/
public static double getRMS(Atom[] ca1, Atom[]ca2,JointFragments frag) throws StructureException {
// now svd ftmp and check if the rms is < X ...
AlternativeAlignment ali = new AlternativeAlignment();
ali.apairs_from_idxlst(frag);
double rms = 999;
int[] idx1 = ali.getIdx1();
int[] idx2 = ali.getIdx2();
Atom[] ca1subset = AlignTools.getFragmentFromIdxList(ca1, idx1);
Atom[] ca2subset = AlignTools.getFragmentFromIdxList(ca2,idx2);
ali.calculateSuperpositionByIdx(ca1,ca2);
Matrix rot = ali.getRotationMatrix();
Atom atom = ali.getShift();
for (Atom a : ca2subset) {
Calc.rotate(a, rot);
Calc.shift(a, atom);
}
rms = Calc.rmsd(ca1subset,ca2subset);
return rms;
}
public boolean angleCheckOk(FragmentPair a, FragmentPair b, float distcutoff){
double dist = -999;
Atom v1 = a.getUnitv();
Atom v2 = b.getUnitv();
dist = Calc.getDistance(v1,v2);
return dist <= distcutoff;
}
private boolean distanceCheckOk(FragmentPair a, FragmentPair b, float fragCompatDist){
double dd ;
Atom c1i = a.getCenter1();
Atom c1j = b.getCenter1();
Atom c2i = a.getCenter2();
Atom c2j = b.getCenter2();
dd = Calc.getDistance(c1i,c1j) - Calc.getDistance(c2i,c2j);
if ( dd < 0) dd = -dd;
return dd <= fragCompatDist;
}
/**
* Calculate the pairwise compatibility of fpairs.
Iterates through a list of fpairs and joins them if
they have compatible rotation and translation parameters.
* @param fraglst FragmentPair[] array
* @param angleDiff angle cutoff
* @param fragCompatDist distance cutoff
* @param maxRefine max number of solutions to keep
* @return JointFragments[]
*/
public JointFragments[] frag_pairwise_compat(FragmentPair[] fraglst, int angleDiff, float fragCompatDist, int maxRefine ){
FragmentPair[] tmpfidx = new FragmentPair[fraglst.length];
for ( int i=0 ; i < fraglst.length; i++){
tmpfidx[i] = (FragmentPair)fraglst[i].clone();
}
int n = tmpfidx.length;
//indicator if a fragment was already used
int[] used = new int[n];
//the final list of joined fragments stores as apairs
List<JointFragments> fll = new ArrayList<JointFragments>();
double adiff = angleDiff * Math.PI / 180d;
logger.finer("addiff" + adiff);
//distance between two unit vectors with angle adiff
double ddiff = Math.sqrt(2.0-2.0*Math.cos(adiff));
logger.finer("ddiff" + ddiff);
// the fpairs in the flist have to be sorted with respect to their positions
while (tmpfidx.length > 0){
int i = 0;
int j = 1;
used[i]=1;
JointFragments f = new JointFragments();
int p1i = tmpfidx[i].getPos1();
int p1j = tmpfidx[i].getPos2();
int maxi = p1i+tmpfidx[i].getLength();
f.add(p1i,p1j,0,tmpfidx[i].getLength());
n = tmpfidx.length;
while ( j < n) {
int p2i = tmpfidx[j].getPos1();
int p2j = tmpfidx[j].getPos2();
int l2 = tmpfidx[j].getLength();
if ( p2i > maxi) {
double dist = Calc.getDistance(tmpfidx[i].getUnitv(), tmpfidx[j].getUnitv());
if ( dist < ddiff) {
// centroids have to be closer than fragCompatDist
double dd = Calc.getDistance(tmpfidx[i].getCenter1(),tmpfidx[j].getCenter1()) -
Calc.getDistance(tmpfidx[i].getCenter2(),tmpfidx[j].getCenter2());
if ( dd < 0)
dd = -dd;
if ( dd < fragCompatDist){
maxi = p2i+l2;
used[j]=1;
f.add(p2i,p2j,0,tmpfidx[j].getLength());
}
}
}
j+=1;
}
int red = 0;
for (int k = 0 ; k < n ; k ++) {
if (used[k] == 0) {
tmpfidx[red] = tmpfidx[k];
red+=1;
}
}
used = new int[n];
tmpfidx = (FragmentPair[])resizeArray(tmpfidx,red);
fll.add(f);
}
Comparator<JointFragments> comp = new JointFragmentsComparator();
Collections.sort(fll,comp);
Collections.reverse(fll);
int m = Math.min(maxRefine,fll.size());
List<JointFragments> retlst = new ArrayList<JointFragments>();
for ( int i = 0 ; i < m ; i++){
JointFragments jf = fll.get(i);
retlst.add(jf);
}
return retlst.toArray(new JointFragments[retlst.size()]);
}
public void extendFragments(Atom[] ca1, Atom[] ca2 ,JointFragments[] fragments, StrucAligParameters params) throws StructureException {
for(JointFragments p : fragments){
extendFragments(ca1, ca2, p, params);
}
}
public void extendFragments(Atom[] ca1, Atom[] ca2 , JointFragments fragments, StrucAligParameters params) throws StructureException {
List<int[]> pos = fragments.getIdxlist();
int[] firstP = pos.get(0);
int pstart1 = firstP[0];
int pstart2 = firstP[1];
int[] lastP = pos.get(pos.size()-1);
int pend1 = lastP[0];
int pend2 = lastP[1];
boolean startReached = false;
boolean endReached = false;
while (! (startReached && endReached)){
if ( ! (startReached) && ((pstart1 <= 0) || (pstart2 <= 0))) {
startReached = true;
} else {
pstart1--;
pstart2--;
}
if ( ! (endReached) && (( pend1 >= (ca1.length -1) ) || ( pend2 >= ca2.length -1 )) ){
endReached = true;
} else {
pend1++;
pend2++;
}
if ( ! startReached){
double newRms1 = testAdd(ca1, ca2, fragments,pstart1,pstart2);
if (newRms1 < 3.7 ) {
fragments.add(pstart1,pstart2,0,1);
} else {
startReached = true;
}
}
if( ! endReached){
double newRms2 = testAdd(ca1, ca2, fragments, pend1, pend2);
if ( newRms2 < 3.7) {
fragments.add(pend1,pend2,0,1);
} else {
endReached = true;
}
}
}
}
private double testAdd(Atom[] ca1, Atom[] ca2, JointFragments fragments, int pstart1, int pstart2) throws StructureException {
JointFragments frag = new JointFragments();
frag.setIdxlist(fragments.getIdxlist());
frag.add(pstart1, pstart2, 0,1);
return FragmentJoiner.getRMS(ca1, ca2, frag);
}
}
class JointFragmentsComparator
implements Comparator<JointFragments>, Serializable {
private static final long serialVersionUID = 1;
@Override
public int compare(JointFragments one, JointFragments two) {
int s1 = one.getIdxlist().size();
int s2 = two.getIdxlist().size();
double rms1 = one.getRms();
double rms2 = two.getRms();
if ( s1 > s2 ) {
return 1 ;
}
else if ( s1 < s2){
return -1;
}
else{
if ( rms1 < rms2)
return 1;
if ( rms1 > rms2)
return -1;
return 0;
}
}
}