package gdsc.smlm.model;
/*-----------------------------------------------------------------------------
* GDSC SMLM Software
*
* Copyright (C) 2013 Alex Herbert
* Genome Damage and Stability Centre
* University of Sussex, UK
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or
* (at your option) any later version.
*---------------------------------------------------------------------------*/
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import org.apache.commons.math3.random.RandomGenerator;
/**
* Contains a model for a compound moving molecule (contains multiple molecules held in a fixed configuration).
* <p>
* The coordinates of the object represent the current centre of mass. The coordinates of each molecule can be retrieved
* using the {@link #getCoordinates(int)} method;
*/
public class CompoundMoleculeModel extends MoleculeModel
{
static final double DEGREES_TO_RADIANS = 180.0 / Math.PI;
public static final double[] X_AXIS = new double[] { 1, 0, 0 };
public static final double[] Y_AXIS = new double[] { 0, 1, 0 };
public static final double[] Z_AXIS = new double[] { 0, 0, 1 };
/**
* Identity matrix for no rotation
*/
private static final double[] I = new double[] { 1, 0, 0, 0, 1, 0, 0, 0, 1 };
private List<? extends MoleculeModel> molecules;
/**
* The fraction of a population that this molecule represents.
* <p>
* This is used when creating a particle distribution of different types of molecules. This is included so that the
* population description can be serialised.
*/
private double fraction = 1;
/**
* The diffusion rate for the molecule
*/
private double diffusionRate = 0;
/**
* The diffusion type for the molecule
*/
private DiffusionType diffusionType = DiffusionType.RANDOM_WALK;
/**
* Create a new molecule
* <p>
* Note: molecules mass may be updated, see {@link #getMass()}
*
* @param id
* @param xyz
* [x,y,z]
* @param molecules
* @throws IllegalArgumentException
* If the input list contains nulls
*/
public CompoundMoleculeModel(int id, double[] xyz, List<? extends MoleculeModel> molecules)
throws IllegalArgumentException
{
this(id, xyz, molecules, true);
}
/**
* Create a new molecule
* <p>
* Note: molecules mass may be updated, see {@link #getMass()}
*
* @param id
* @param x
* @param y
* @param z
* @param molecules
* @throws IllegalArgumentException
* If the input list contains nulls
*/
public CompoundMoleculeModel(int id, double x, double y, double z, List<? extends MoleculeModel> molecules)
throws IllegalArgumentException
{
this(id, x, y, z, molecules, true);
}
/**
* Create a new molecule
* <p>
* Note: molecules mass may be updated, see {@link #getMass()}
*
* @param id
* @param xyz
* [x,y,z]
* @param molecules
* @param centre
* Shift molecules to have a centre-of-mass at 0,0,0.
* @throws IllegalArgumentException
* If the input list contains nulls
*/
public CompoundMoleculeModel(int id, double[] xyz, List<? extends MoleculeModel> molecules, boolean centre)
throws IllegalArgumentException
{
super(id, xyz);
setMolecules(molecules, centre);
}
/**
* Create a new molecule
* <p>
* Note: molecules mass may be updated, see {@link #getMass()}
*
* @param id
* @param x
* @param y
* @param z
* @param molecules
* @param centre
* Shift molecules to have a centre-of-mass at 0,0,0.
* @throws IllegalArgumentException
* If the input list contains nulls
*/
public CompoundMoleculeModel(int id, double x, double y, double z, List<? extends MoleculeModel> molecules,
boolean centre) throws IllegalArgumentException
{
super(id, x, y, z);
setMolecules(molecules, centre);
}
private void setMolecules(List<? extends MoleculeModel> molecules, boolean centre) throws IllegalArgumentException
{
if (molecules == null)
molecules = new ArrayList<MoleculeModel>(0);
this.molecules = molecules;
checkMass();
if (centre)
centre();
}
/**
* Check that all the molecules have a valid mass. If any are below zero then set to zero. If all are below zero
* then set to 1 so that the centre-of-mass is valid.
*/
public void checkMass()
{
int invalidMass = 0;
for (MoleculeModel m : molecules)
{
if (m == null)
throw new IllegalArgumentException("Input list contains null molecules");
if (m.mass <= 0 || Double.isNaN(m.mass))
{
invalidMass++;
//throw new IllegalArgumentException("Input list contains molecules with invalid mass: " + m.mass);
}
}
this.mass = 0;
if (invalidMass > 0)
{
final double resetMass = (invalidMass == molecules.size()) ? 1 : 0;
for (MoleculeModel m : molecules)
{
if (m.mass <= 0 || Double.isNaN(m.mass))
m.mass = resetMass;
this.mass += m.mass;
}
}
else
{
for (MoleculeModel m : molecules)
{
this.mass += m.mass;
}
}
}
/**
* Shift molecules to have a centre-of-mass at 0,0,0.
* Coordinates are therefore relative to the origin.
*/
public void centre()
{
if (molecules.isEmpty())
return;
if (this.mass == 0)
checkMass();
double[] com = new double[3];
for (MoleculeModel m : molecules)
{
for (int i = 0; i < 3; i++)
com[i] += m.xyz[i] * m.mass;
}
// Checked during initialisation
//if (w <= 0 || Double.isNaN(w))
// throw new IllegalArgumentException("Input list contains molecules with invalid mass");
for (int i = 0; i < 3; i++)
{
com[i] /= this.mass;
}
for (MoleculeModel m : molecules)
{
for (int i = 0; i < 3; i++)
m.xyz[i] -= com[i];
}
}
/**
* Rotate the molecule using a random axis and a random rotation angle. The rotation is around the centre-of-mass.
*
* @param maxAngle
* The maximum angle to rotate (in either direction) in degrees
* @param random
*/
public void rotateRandom(double maxAngle, RandomGenerator random)
{
if (maxAngle == 0 || getSize() < 2)
return;
if (this.mass == 0)
checkMass();
double[] axis = new double[3];
for (int i = 0; i < 3; i++)
{
axis[i] = random.nextGaussian();
}
final double angle = (-maxAngle + random.nextDouble() * 2.0 * maxAngle);
if (angle == 0)
return;
rotateMolecules(axis, angle);
}
/**
* Rotate the molecule using a specified axis and a random rotation angle. The rotation is around the
* centre-of-mass.
*
* @param axis
* The axis to rotate around
* @param maxAngle
* The maximum angle to rotate (in either direction) in degrees
* @param random
*/
public void rotateRandomAngle(double[] axis, double maxAngle, RandomGenerator random)
{
if (maxAngle == 0 || getSize() < 2)
return;
if (this.mass == 0)
checkMass();
final double angle = (-maxAngle + random.nextDouble() * 2.0 * maxAngle);
if (angle == 0)
return;
rotateMolecules(axis, angle);
}
/**
* Rotate the molecule using a random axis and a specified rotation angle. The rotation is around the
* centre-of-mass.
*
* @param axis
* The axis to rotate around
* @param maxAngle
* The maximum angle to rotate (in either direction) in degrees
* @param random
*/
public void rotateRandomAxis(double angle, RandomGenerator random)
{
if (angle == 0 || getSize() < 2)
return;
if (this.mass == 0)
checkMass();
double[] axis = new double[3];
for (int i = 0; i < 3; i++)
{
axis[i] = random.nextGaussian();
}
rotateMolecules(axis, angle);
}
/**
* Rotate the molecule using a specified axis and rotation angle. The rotation is around the centre-of-mass.
*
* @param angle
* The angle to rotate (in degrees)
* @param axis
* The axis to rotate around
*/
public void rotate(double[] axis, double angle)
{
if (angle == 0 || getSize() < 2 || axis == null || axis.length < 3)
return;
if (this.mass == 0)
checkMass();
rotateMolecules(axis, angle);
}
/**
* Rotate the molecule using a specified axis and rotation angle. The rotation is around the centre-of-mass.
*
* @param angle
* The angle to rotate (in degrees)
* @param axis
* The axis to rotate around
*/
private void rotateMolecules(double[] axis, double angle)
{
double[] r = getRotationMatrix(axis, angle);
if (r == I)
return;
// Use the actual centre of mass of the molecules to avoid drift over time
// (i.e. do not assume the centre is 0,0,0)
double[] com = new double[3];
for (MoleculeModel m : molecules)
{
for (int i = 0; i < 3; i++)
com[i] += m.xyz[i] * m.mass;
}
for (int i = 0; i < 3; i++)
{
com[i] /= this.mass;
}
for (MoleculeModel m : molecules)
{
final double xtmp = m.xyz[0] - com[0];
final double ytmp = m.xyz[1] - com[1];
final double ztmp = m.xyz[2] - com[2];
// Do not add back COM since the molecules should be centred around the origin
m.xyz[0] = xtmp * r[0] + ytmp * r[1] + ztmp * r[2];
m.xyz[1] = xtmp * r[3] + ytmp * r[4] + ztmp * r[5];
m.xyz[2] = xtmp * r[6] + ytmp * r[7] + ztmp * r[8];
}
}
/**
* Get the rotation matrix for a rotation around an axis
*
* @param axis
* axis
* @param angle
* angle (in degrees)
* @return The rotation matrix
*/
private double[] getRotationMatrix(double[] axis, double angle)
{
/* Set to unit length */
double length = 0;
for (int i = 0; i < 3; i++)
{
length += axis[i] * axis[i];
}
if (length == 0 || Double.isNaN(length))
return I;
length = Math.sqrt(length);
for (int i = 0; i < 3; i++)
{
axis[i] /= length;
}
//System.out.printf("Angle = %.1f, Axis = %.3f, %.3f, %.3f\n", angle, axis[0], axis[1], axis[2]);
/* Store the components and their squares */
final double u = axis[0];
final double u2 = u * u;
final double v = axis[1];
final double v2 = v * v;
final double w = axis[2];
final double w2 = w * w;
final double uv = u * v;
final double uw = u * w;
final double vw = v * w;
angle /= DEGREES_TO_RADIANS;
final double cost = Math.cos(angle);
final double sint = Math.sin(angle);
final double[] r = new double[9];
/* Set the rotation matrix */
r[0] = u2 + ((v2 + w2) * cost);
r[1] = uv - uv * cost - w * sint;
r[2] = uw - uw * cost + v * sint;
r[3] = uv - uv * cost + w * sint;
r[4] = v2 + ((u2 + w2) * cost);
r[5] = vw - vw * cost - u * sint;
r[6] = -1.0 * (uw * (-1.0 + cost)) - v * sint;
r[7] = -1.0 * (vw * (-1.0 + cost)) + u * sint;
r[8] = w2 + ((u2 + v2) * cost);
/*
* r[0] = u2 + ((v2 + w2)*cost);
* r[3] = uv - uv*cost - w*sint;
* r[6] = uw - uw*cost + v*sint;
* r[1] = uv - uv*cost + w*sint;
* r[4] = v2 + ((u2 + w2)*cost);
* r[7] = vw - vw*cost - u*sint;
* r[2] = -1.0*(uw*(-1.0 + cost)) - v*sint;
* r[5] = -1.0*(vw*(-1.0 + cost)) + u*sint;
* r[8] = w2 + ((u2 + v2)*cost);
*/
return r;
}
/**
* @return The number of molecules in the compound molecule
*/
public int getSize()
{
return molecules.size();
}
/**
* Get the current coordinates of the nth molecule in the compound molecule
*
* @param n
* The requested molecule (0 <= n < {@link #getSize()})
* @return The xyz coordinates
* @see #getSize()
* @throws IndexOutOfBoundsException
* If the requested molecule does not exist
*/
public double[] getCoordinates(int n) throws IndexOutOfBoundsException
{
double[] xyz = Arrays.copyOf(this.xyz, 3);
MoleculeModel m = molecules.get(n);
for (int i = 0; i < 3; i++)
xyz[i] += m.xyz[i];
return xyz;
}
/**
* Get the current coordinates of the nth molecule in the compound molecule relative to the centre-of-mass
*
* @param n
* The requested molecule (0 <= n < {@link #getSize()})
* @return The xyz coordinates
* @see #getSize()
* @throws IndexOutOfBoundsException
* If the requested molecule does not exist
*/
public double[] getRelativeCoordinates(int n) throws IndexOutOfBoundsException
{
MoleculeModel m = molecules.get(n);
return Arrays.copyOf(m.xyz, 3);
}
/**
* Get the nth molecule in the compound molecule
* <p>
* Note that the molecule coordinates are relative the centre-of-mass of the compound
*
* @param n
* The requested molecule (0 <= n < {@link #getSize()})
* @return The molecule
* @see #getSize()
* @throws IndexOutOfBoundsException
* If the requested molecule does not exist
*/
public MoleculeModel getMolecule(int n) throws IndexOutOfBoundsException
{
return molecules.get(n);
}
/**
* @return The fraction of a population that this molecule represents
*/
public double getFraction()
{
return fraction;
}
/**
* Set the fraction of a population that this molecule represents
*
* @param fraction
* the fraction to set
*/
public void setFraction(double fraction)
{
this.fraction = fraction;
}
/**
* Return the of all the molecules.
* <p>
* The mass is calculated once during initialisation. If any molecule has a mass less than zero it will be set to
* zero. If all molecules have a mass of zero then the mass for each will be reset to one. This allows the
* centre-of-mass calculation to function correctly. If a molecule is part of the compound and has no mass it will
* be rotated and moved but will not contribute to the COM.
*
* @return The mass of all the molecules
*/
public double getMass()
{
return mass;
}
/**
* Scale the molecules relative coordinates by the given factor
*
* @param factor
*/
public void scale(double factor)
{
for (MoleculeModel m : molecules)
{
for (int i = 0; i < 3; i++)
m.xyz[i] *= factor;
}
}
/**
* @return the diffusionRate
*/
public double getDiffusionRate()
{
return diffusionRate;
}
/**
* @param diffusionRate
* the diffusionRate to set
*/
public void setDiffusionRate(double diffusionRate)
{
this.diffusionRate = diffusionRate;
}
/**
* Get the diffusion type
*
* @return The diffusion type
*/
public DiffusionType getDiffusionType()
{
return diffusionType;
}
/**
* Set the diffusion type
*
* @param diffusionType
* The diffusion type
*/
public void setDiffusionType(DiffusionType diffusionType)
{
this.diffusionType = diffusionType;
}
}