MathUtils.java

Go to the documentation of this file.

/*
 *   DENOPTIM
 *   Copyright (C) 2019 Vishwesh Venkatraman <vishwesh.venkatraman@ntnu.no> and
 *   Marco Foscato <marco.foscato@uib.no>
 * 
 *   This program is free software: you can redistribute it and/or modify
 *   it under the terms of the GNU Affero General Public License as published
 *   by the Free Software Foundation, either version 3 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU Affero General Public License for more details.
 *
 *   You should have received a copy of the GNU Affero General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
package denoptim.utils;
 
import java.math.BigDecimal;
import java.math.RoundingMode;
import java.util.ArrayList;
import java.util.Collection;
import java.util.List;
 
import javax.vecmath.Point3d;
import javax.vecmath.Vector3d;
 
import denoptim.constants.DENOPTIMConstants;
 
public class MathUtils
{
 
//------------------------------------------------------------------------------
 
    public static double calculateAngle (double[] a, double[] b)
    {
        double alen = length(a);
        double blen = length(b);
        double abdot = computeDotProduct(a, b);
        double angle = Math.acos(abdot/(alen*blen));
        return Math.toDegrees(angle);
    }
    
//------------------------------------------------------------------------------    
    public static void norm(double[] d, double[] ret) 
    {
        double length;
        length = length(d);
        ret[0] = d[0] / length;
        ret[1] = d[1] / length;
        ret[2] = d[2] / length;
    }
 
//------------------------------------------------------------------------------
 
    public static Vector3d normDist(Point3d p1, Point3d p2)
    {
    Vector3d v = new Vector3d(p1.x - p2.x,
                  p1.y - p2.y,
                  p1.z - p2.z);
    v.normalize();
    return v;
    }
 
//------------------------------------------------------------------------------
    public static double[] scale(double[] d, double a)
    {
    double[] scaled = new double[3];
        scaled[0] = d[0] * a;
        scaled[1] = d[1] * a;
        scaled[2] = d[2] * a;
    return scaled;
    }
    
//------------------------------------------------------------------------------    
    public static void scalar(double[] d, double a, double[] ret) 
    {
        ret[0] = d[0] * a;
        ret[1] = d[1] * a;
        ret[2] = d[2] * a;
    }    
 
//------------------------------------------------------------------------------
 
    public static double length (double[] v)
    {
        return Math.sqrt (v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
    }
 
//------------------------------------------------------------------------------
 
    public static double computeDotProduct (double[] v0, double[] v1)
    {
        return v0[0] * v1[0] + v0[1] * v1[1] + v0[2] * v1[2];
    }
 
//------------------------------------------------------------------------------
 
    public static double[] computeCrossProduct (double[] v0, double[] v1)
    {
        double[] crossProduct = new double[3];
        crossProduct[0] = v0[1] * v1[2] - v0[2] * v1[1];
        crossProduct[1] = v0[2] * v1[0] - v0[0] * v1[2];
        crossProduct[2] = v0[0] * v1[1] - v0[1] * v1[0];
        return crossProduct;
    }
 
//------------------------------------------------------------------------------
 
    public static double[] subtract (double[] v0, double[] v1)
    {
        double[] res = new double[3];
        res[0] = v0[0] - v1[0];
        res[1] = v0[1] - v1[1];
        res[2] = v0[2] - v1[2];
        return res;
    }
 
//------------------------------------------------------------------------------
 
    public static double[] add (double[] v0, double[] v1)
    {
        double[] res = new double[3];
        res[0] = v0[0] + v1[0];
        res[1] = v0[1] + v1[1];
        res[2] = v0[2] + v1[2];
        return res;
    }
 
//------------------------------------------------------------------------------
 
    public static double computeDihedralAngle(Point3d p0, Point3d p1, 
                        Point3d p2, Point3d p3)
    {
        double[] v0 = new double[] {p0.x, p0.y, p0.z};
        double[] v1 = new double[] {p1.x, p1.y, p1.z};
        double[] v2 = new double[] {p2.x, p2.y, p2.z};
        double[] v3 = new double[] {p3.x, p3.y, p3.z};
 
    return computeDihedralAngle(v0, v1, v2, v3);
    }
 
//------------------------------------------------------------------------------
 
    public static double computeDihedralAngle (double[] v0, double[] v1,
                                                double[] v2, double[] v3)
    {
        double[] r1 = subtract(v1, v0);
        double[] r2 = subtract(v2, v1);
        double[] r3 = subtract(v3, v2);
 
        double[] n1 = computeCrossProduct(r1, r2);
        double[] n2 = computeCrossProduct(r2, r3);
 
        double psin = computeDotProduct(n1, r3) *
                            Math.sqrt(computeDotProduct(r2, r2));
        double pcos = computeDotProduct(n1, n2);
 
        return Math.toDegrees(Math.atan2(psin, pcos));
    }
 
//------------------------------------------------------------------------------
 
    public static double log2(double x)
    {
        return Math.log(x)/Math.log(2);
    }
 
 
//------------------------------------------------------------------------------
 
    public static double roundValue(double val, int decimalPlaces)
    {
        BigDecimal bd = new BigDecimal(val);
        bd = bd.setScale(decimalPlaces, RoundingMode.UP);
        return bd.doubleValue();
    }
 
//------------------------------------------------------------------------------
 
    public static double angle(Point3d a, Point3d b, Point3d c)
    {
        double xba = a.x - b.x;
        double xbc = c.x - b.x;
        double yba = a.y - b.y;
        double ybc = c.y - b.y;
        double zba = a.z - b.z;
        double zbc = c.z - b.z;
 
        double ba = Math.sqrt(xba*xba + yba*yba + zba*zba);
        double bc = Math.sqrt(xbc*xbc + ybc*ybc + zbc*zbc);
 
        double dot = xba*xbc + yba*ybc + zba*zbc;
        dot /= (ba * bc);
 
        double angle;
    if (dot < 0.0 && 
       Math.abs(Math.abs(dot)-1.0) < DENOPTIMConstants.FLOATCOMPARISONTOLERANCE)
    {
        angle = 180.0;
    }
    else
    {
        angle = Math.toDegrees(Math.acos(dot));
    }
 
        return angle;
    }
    
//------------------------------------------------------------------------------
    
    public static double getAngle(double[] A, double[] B, double[] C)
    {
        double[] vBA = {A[0]-B[0], A[1]-B[1], A[2]-B[2]}; // Vector from B to A
        double[] vBC = {C[0]-B[0], C[1]-B[1], C[2]-B[2]}; // Vector from B to C
 
    return calculateAngle(vBA, vBC);
    }
 
    
//------------------------------------------------------------------------------
 
    public static double distance(Point3d a, Point3d b)
    {
        double dx = a.x - b.x, dy = a.y - b.y, dz = a.z - b.z;
        return Math.sqrt(dx*dx + dy*dy + dz*dz);
    }
 
//------------------------------------------------------------------------------
 
    private static final double[] tp1 = new double[3];
    private static final double[] tp2 = new double[3];
    private static final double[] tp3 = new double[3];
    private static final double[] tp4 = new double[3];
 
    public static double torsion(Point3d p1, Point3d p2, Point3d p3, Point3d p4)
    {
        p1.get(tp1);
        p2.get(tp2);
        p3.get(tp3);
        p4.get(tp4);
 
    return computeDihedralAngle(tp1, tp2, tp3, tp4);
    }
    
//------------------------------------------------------------------------------    
    
    public static double sum(double[] a)
    {
        double sum = 0.0;
        for (double v : a) {
            sum += v;
        }
        return sum;
    }
 
//------------------------------------------------------------------------------        
    
    public static double mean(double[] a) 
    {
        if (a.length == 0) 
            return Double.NaN;
        double sum = sum(a);
        return sum / a.length;
    }
    
//------------------------------------------------------------------------------
    
    public static double mean(List<Double> vals) 
    {
        double[] arr = new double[vals.size()];
        for (int i=0; i<vals.size(); i++)
        {
            arr[i] = vals.get(i);
        }
        
        return mean(arr);
    }
    
//------------------------------------------------------------------------------    
 
    public static Vector3d getNormalDirection(Vector3d vecA)
    {
        Vector3d normalDir = new Vector3d();
        ArrayList<Vector3d> candidates = new ArrayList<>();
        candidates.add(new Vector3d(1.0, 0.0, 0.0));
        candidates.add(new Vector3d(0.0, 1.0, 0.0));
        candidates.add(new Vector3d(0.0, 0.0, 1.0));
        candidates.add(new Vector3d(1.0/Math.sqrt(2.0), 
                    1.0/Math.sqrt(2.0),
                    0.0));
        candidates.add(new Vector3d(1.0/Math.sqrt(2.0), 
                    0.0,
                    1.0/Math.sqrt(2.0)));
        candidates.add(new Vector3d(0.0, 
                    1.0/Math.sqrt(2.0), 
                    1.0/Math.sqrt(2.0)));
 
        for (Vector3d candidate : candidates) {
            double res = vecA.dot(candidate);
            if (Math.abs(res) > 0.1 && Math.abs(res) < 0.9) {
                normalDir.cross(vecA, candidate);
                break;
            }
        }
    normalDir.normalize();
 
        return normalDir;
    }
    
//------------------------------------------------------------------------------
 
    public static double[] centroidOf(Collection<double[]> points, int dimension) 
    {
        double[] centroid = new double[dimension];
        for (double[] point : points) 
        {
            for (int i = 0; i < centroid.length; i++) 
            {
                centroid[i] += point[i];
            }
        }
        for (int i = 0; i < centroid.length; i++) 
        {
            centroid[i] /= points.size();
        }
        return centroid;
    }
 
//------------------------------------------------------------------------------
}