added documentation, tests, and rounding to closes 0.05 value

common/lib/capa/capa/chem/miller.py

+""" Calculation of Miller indices """
+
 import numpy as np
 import math
 import fractions as fr
@@ -7,150 +9,256 @@ import json
 
 
 def lcm(a, b):
-    """ return lcm of a,b """
+    """
+    Returns least common multiple of a, b
+
+    Args:
+        a, b: floats
+
+    Returns:
+        float
+    """
     return a * b / fr.gcd(a, b)
 
 
-def section_to_fraction(distance):
-    """ Convert float distance, that plane cut on axis
-    to fraction. Return inverted fraction
+def segment_to_fraction(distance):
+    """
+    Converts lengths of which the plane cuts the axes to fraction.
+
+    Tries convert distance to closest nicest fraction with denominator less or
+    equal than 10. It is
+    purely for simplicity and clearance of learning purposes. Jenny: 'In typical
+    courses students usually do not encounter indices any higher than 6'.
+
+    If distance is not a number (numpy nan), it means that plane is parallel to
+    axis or contains it. Inverted fraction to nan (nan is 1/0) = 0 / 1 is
+    returned
+
+    Generally (special cases):
+
+    a) if distance is smaller than some constant, i.g. 0.01011,
+    than fraction's denominator usually much greater than 10.
+
+    b) Also, if student will set point on 0.66 -> 1/3, so it is 333 plane,
+    But if he will slightly move the mouse and click on 0.65 -> it will be
+    (16,15,16) plane. That's why we are doing adjustments for points coordinates,
+    to the closest tick, tick + tick / 2 value. And now UI sends to server only
+    values multiple to 0.05  (half of tick). Same rounding is implemented for
+    unittests.
+
+    But if one will want to calculate miller indices with exact coordinates and
+    with nice fractions (which produce small Miller indices), he may want shift
+    to new origin if segments are like S = (0.015, > 0.05, >0.05) - close to zero
+    in one coordinate. He may update S to (0, >0.05, >0.05) and shift origin.
+    In this way he can recieve nice small fractions. Also there is can be
+    degenerated case when S = (0.015, 0.012, >0.05) - if update S to (0, 0, >0.05) -
+    it is a line. This case  should be considered separately. Small nice Miller
+    numbers and possibility to create very small segments can not be implemented
+    at same time).
+
+
+    Args:
+        distance: float distance that plane cuts on axis, it must not be 0.
+        Distance is multiple of 0.05.
+
+    Returns:
+         Inverted fraction.
+         0 / 1 if distance is nan
 
     """
-    # import ipdb; ipdb.set_trace()
-    if np.isnan(distance):  # plane || to axis (or contains axis)
-        print distance, 0
-        # return inverted fration to a == nan == 1/0  => 0 / 1
+    if np.isnan(distance):
         return fr.Fraction(0, 1)
-    elif math.fabs(distance) <= 0.05:  # plane goes through origin, 0.02 - UI delta
-        return fr.Fraction(1 if distance >= 0 else -1, 1)  # ERROR, need shift of coordinates
     else:
-        # limit_denominator to closest nicest fraction
-        # import ipdb; ipdb.set_trace()
-        fract = fr.Fraction(distance).limit_denominator(10)  # 5 / 2 : numerator / denominator
-        print 'Distance', distance, 'Inverted fraction', fract
-        # return inverted fraction
+        fract = fr.Fraction(distance).limit_denominator(10)
         return fr.Fraction(fract.denominator, fract.numerator)
 
 
-def sub_miller(sections):
-    ''' Calculate miller indices.
-        Plane does not intersect origin
+def sub_miller(segments):
     '''
-    fracts = [section_to_fraction(section) for section in sections]
+    Calculates Miller indices from segments.
 
-    print sections, fracts
+    Algorithm:
 
-    common_denominator = reduce(lcm, [fract.denominator for fract in fracts])
-    print 'common_denominator:', common_denominator
+    1. Obtain inverted fraction from segments
+
+    2. Find common denominator of inverted fractions
+
+    3. Lead fractions to common denominator and throws denominator away.
+
+    4. Return obtained values.
 
-    # 2) lead to a common denominator
-    # 3) throw away denominator
-    miller = [fract.numerator * math.fabs(common_denominator) / fract.denominator for fract in fracts]
+    Args:
+        List of 3 floats, meaning distances that plane cuts on x, y, z axes.
+        Any float not equals zero, it means that plane does not intersect origin,
+        i. e. shift of origin has already been done.
 
-    # import ipdb; ipdb.set_trace()
-    # nice output:
-    # output = '(' + ''.join(map(str, map(decimal.Decimal, miller))) + ')'
-    # import ipdb; ipdb.set_trace()
-    output = '(' + ','.join(map(str, map(decimal.Decimal, miller))) + ')'
-    # print 'Miller indices:', output
-    return output
+    Returns:
+        String that represents Miller indices, e.g: (-6,3,-6) or (2,2,2)
+    '''
+    fracts = [segment_to_fraction(segment) for segment in segments]
+    common_denominator = reduce(lcm, [fract.denominator for fract in fracts])
+    miller = ([fract.numerator * math.fabs(common_denominator) /
+                                fract.denominator for fract in fracts])
+    return'(' + ','.join(map(str, map(decimal.Decimal, miller))) + ')'
 
 
 def miller(points):
-    """Calculate miller indices of plane
     """
+    Calculates Miller indices from points.
 
-    # print "\nCalculating miller indices:"
-    # print 'Points:\n', points
-    # calculate normal to plane
-    N = np.cross(points[1] - points[0], points[2] - points[0])
-    # print "Normal:", N
+    Algorithm:
 
-    # origin
-    O = np.array([0, 0, 0])
+    1. Calculate normal vector to a plane that goes trough all points.
+
+    2. Set origin.
+
+    3. Create Cartesian coordinate system (Ccs).
 
-    # point of plane
-    P = points[0]
-
-    # equation of a line for axes: O + (B-O) * t
-    # t - parameters, B = [Bx, By, Bz]:
-    B = map(np.array, [[1.0, 0, 0], [0, 1.0, 0], [0, 0, 1.0]])
-
-    # coordinates of intersections with axis:
-    sections = [np.dot(P - O, N) / np.dot(B_axis, N) if np.dot(B_axis, N) != 0 else np.nan for B_axis in B]
-    # import ipdb; ipdb.set_trace()
-
-    if any(x == 0 for x in sections):  # Plane goes through origin.
-        # Need to shift plane out of origin.
-        # For this change origin position
-
-        # 1) find cube vertex, not crossed by plane
-        vertices = [
-            # top
-            np.array([1.0, 1.0, 1.0]),
-            np.array([0.0, 0.0, 1.0]),
-            np.array([1.0, 0.0, 1.0]),
-            np.array([0.0, 1.0, 1.0]),
-            # bottom, except 0,0,0
-            np.array([1.0, 0.0, 0.0]),
-            np.array([0.0, 1.0, 0.0]),
-            np.array([1.0, 1.0, 1.0]),
-                    ]
+    4. Find the lengths of segments of which the plane cuts the axes. Equation
+       of a line for axes: Origin + (Coordinate_vector - Origin) * parameter.
 
+    5. If plane goes trough Origin:
+
+       a) Find new random origin: find unit cube vertex, not crossed by a plane.
+
+       b) Repeat 2-4.
+
+       c) Fix signs of segments after Origin shift. This means to consider
+          original directions of axes. I.g.: Origin was 0,0,0 and became
+          new_origin. If new_origin has same Y coordinate as Origin, then segment
+          does not change its sign. But if new_origin has another Y coordinate than
+          origin (was 0, became 1), than segment has to change its sign (it now
+          lies on negative side of Y axis). New Origin 0 value of X or Y or Z
+          coordinate means that segment does not change sign, 1 value -> does
+          change. So new sign  is (1 - 2 * new_origin): 0 -> 1, 1 -> -1
+
+    6. Run function that calculates miller indices from segments.
+
+    Args:
+        List of points. Each point is list of float coordinates. Order of
+        coordinates in point's list: x, y, z. Points are different!
+
+    Returns:
+        String that represents Miller indices, e.g: (-6,3,-6) or (2,2,2)
+    """
+
+    N = np.cross(points[1] - points[0], points[2] - points[0])
+    O = np.array([0, 0, 0])
+    P = points[0]  # point of plane
+    Ccs = map(np.array, [[1.0, 0, 0], [0, 1.0, 0], [0, 0, 1.0]])
+    segments = ([np.dot(P - O, N) / np.dot(ort, N) if np.dot(ort, N) != 0 else
+                                            np.nan for ort in Ccs])
+    if any(x == 0 for x in segments):  # Plane goes through origin.
+        vertices = [  # top:
+                    np.array([1.0, 1.0, 1.0]),
+                    np.array([0.0, 0.0, 1.0]),
+                    np.array([1.0, 0.0, 1.0]),
+                    np.array([0.0, 1.0, 1.0]),
+                    # bottom, except 0,0,0:
+                    np.array([1.0, 0.0, 0.0]),
+                    np.array([0.0, 1.0, 0.0]),
+                    np.array([1.0, 1.0, 1.0]),
+                ]
         for vertex in vertices:
-            if np.dot(vertex - O, N) != 0:  # vertex in plane
+            if np.dot(vertex - O, N) != 0:  # vertex not in plane
                 new_origin = vertex
                 break
-
-        # get axis with center in new origin
+        # obtain new axes with center in new origin
         X = np.array([1 - new_origin[0], new_origin[1], new_origin[2]])
         Y = np.array([new_origin[0], 1 - new_origin[1], new_origin[2]])
         Z = np.array([new_origin[0], new_origin[1], 1 - new_origin[2]])
+        new_Ccs = [X - new_origin, Y - new_origin, Z - new_origin]
+        segments = ([np.dot(P - new_origin, N) / np.dot(ort, N) if
+                    np.dot(ort, N) != 0 else np.nan for ort in new_Ccs])
+        # fix signs of indices: 0 -> 1, 1 -> -1 (
+        segments = (1 - 2 * new_origin) * segments
 
-        # 2) calculate miller indexes by new origin
-        new_axis = [X - new_origin, Y - new_origin, Z - new_origin]
-        # coordinates of intersections with axis:
-        sections = [np.dot(P - new_origin, N) / np.dot(B_axis, N) if np.dot(B_axis, N) != 0 else np.nan for B_axis in new_axis]
-        # 3) fix signs of indexes
-        #  0 -> 1, 1 -> -1
-        sections = (1 - 2 * new_origin) * sections
-    return sub_miller(sections)
+    return sub_miller(segments)
 
 
 def grade(user_input, correct_answer):
     '''
-    Format:
-    user_input: {"lattice":"sc","points":[["0.77","0.00","1.00"],["0.78","1.00","0.00"],["0.00","1.00","0.72"]]}
-    "lattice" is one of: "", "sc", "bcc", "fcc"
-    correct_answer: {'miller': '(00-1)', 'lattice': 'bcc'}
+    Grade crystallography problem.
+
+    Returns true if lattices are the same and Miller indices are same or minus
+    same. E.g. (2,2,2) = (2, 2, 2) or (-2, -2, -2). Because sign depends only
+    on student's selection of origin.
+
+    Args:
+        user_input, correct_answer: json. Format:
+
+        user_input: {"lattice":"sc","points":[["0.77","0.00","1.00"],
+        ["0.78","1.00","0.00"],["0.00","1.00","0.72"]]}
+
+        correct_answer: {'miller': '(00-1)', 'lattice': 'bcc'}
+
+        "lattice" is one of: "", "sc", "bcc", "fcc"
+
+    Returns:
+        True or false.
     '''
-    def negative(s):
-        # import ipdb; ipdb.set_trace()
+    def negative(m):
+        """
+        Change sign of Miller indices.
+
+        Args:
+            m: string with meaning of Miller indices. E.g.:
+            (-6,3,-6) -> (6, -3, 6)
+
+        Returns:
+            String with changed signs.
+        """
         output = ''
         i = 1
-        while i in range(1, len(s) - 1):
-            if s[i] in (',', ' '):
-                output += s[i]
-            elif s[i] not in ('-', '0'):
-                output += '-' + s[i]
-            elif s[i] == '0':
-                output += s[i]
+        while i in range(1, len(m) - 1):
+            if m[i] in (',', ' '):
+                output += m[i]
+            elif m[i] not in ('-', '0'):
+                output += '-' + m[i]
+            elif m[i] == '0':
+                output += m[i]
             else:
                 i += 1
-                output += s[i]
+                output += m[i]
             i += 1
-            # import ipdb; ipdb.set_trace()
         return '(' + output + ')'
 
+    def round0_25(point):
+        """
+        Rounds point coordinates to closest 0.5 value.
+
+        Args:
+            point: list of float coordinates. Order of coordinates: x, y, z.
+
+        Returns:
+            list of coordinates rounded to closes 0.5 value
+        """
+        rounded_points = []
+        for coord in point:
+            base = math.floor(coord * 10)
+            fractional_part = (coord * 10 - base)
+            aliquot0_25 = math.floor(fractional_part / 0.25)
+            if aliquot0_25 == 0.0:
+                rounded_points.append(base / 10)
+            if aliquot0_25 in (1.0, 2.0):
+                rounded_points.append(base / 10 + 0.05)
+            if aliquot0_25 == 3.0:
+                rounded_points.append(base / 10 + 0.1)
+        return rounded_points
+
     user_answer = json.loads(user_input)
 
     if user_answer['lattice'] != correct_answer['lattice']:
         return False
 
     points = [map(float, p) for p in user_answer['points']]
+
+    # round point to closes 0.05 value
+    points = [round0_25(point) for point in points]
+
     points = [np.array(point) for point in points]
-    # import ipdb; ipdb.set_trace()
-    print miller(points), (correct_answer['miller'].replace(' ', ''), negative(correct_answer['miller']).replace(' ', ''))
+    print miller(points), (correct_answer['miller'].replace(' ', ''),
+        negative(correct_answer['miller']).replace(' ', ''))
 
     if miller(points) in (correct_answer['miller'].replace(' ', ''), negative(correct_answer['miller']).replace(' ', '')):
         return True
@@ -159,7 +267,7 @@ def grade(user_input, correct_answer):
 
 
 class Test_Crystallography_Miller(unittest.TestCase):
-    ''' test crystallography grade function '''
+    ''' Tests  for crystallography grade function.'''
 
     def test_1(self):
         user_input = '{"lattice": "bcc", "points": [["0.50", "0.00", "0.00"], ["0.00", "0.50", "0.00"], ["0.00", "0.00", "0.50"]]}'
@@ -178,8 +286,10 @@ class Test_Crystallography_Miller(unittest.TestCase):
         self.assertTrue(grade(user_input, {'miller': '(-3, 3, -3)', 'lattice': 'bcc'}))
 
     def test_5(self):
+        """ return true only in case points coordinates are exact.
+        But if they transform to closest 0.05 value it is not true"""
         user_input = '{"lattice": "bcc", "points": [["0.33", "1.00", "0.00"], ["0.00", "0.33", "0.00"], ["0.00", "1.00", "0.33"]]}'
-        self.assertTrue(grade(user_input, {'miller': '(-6,3,-6)', 'lattice': 'bcc'}))
+        self.assertFalse(grade(user_input, {'miller': '(-6,3,-6)', 'lattice': 'bcc'}))
 
     def test_6(self):
         user_input = '{"lattice": "bcc", "points": [["0.00", "0.25", "0.00"], ["0.25", "0.00", "0.00"], ["0.00", "0.00", "0.25"]]}'
@@ -261,6 +371,20 @@ class Test_Crystallography_Miller(unittest.TestCase):
         user_input = u'{"lattice":"","points":[["0.00","0.00","0.01"],["1.00","1.00","0.01"],["0.00","1.00","1.00"]]}'
         self.assertTrue(grade(user_input, {'miller': '(1,-1,1)', 'lattice': ''}))
 
+    def test_26(self):
+        user_input = u'{"lattice":"","points":[["0.00","0.01","0.00"],["1.00","0.00","0.00"],["0.00","0.00","1.00"]]}'
+        self.assertTrue(grade(user_input, {'miller': '(0,-1,0)', 'lattice': ''}))
+
+    def test_27(self):
+        """ rounding to 0.35"""
+        user_input = u'{"lattice":"","points":[["0.33","0.00","0.00"],["0.00","0.33","0.00"],["0.00","0.00","0.33"]]}'
+        self.assertTrue(grade(user_input, {'miller': '(3,3,3)', 'lattice': ''}))
+
+    def test_28(self):
+        """ rounding to 0.30"""
+        user_input = u'{"lattice":"","points":[["0.30","0.00","0.00"],["0.00","0.30","0.00"],["0.00","0.00","0.30"]]}'
+        self.assertTrue(grade(user_input, {'miller': '(10,10,10)', 'lattice': ''}))
+
     def test_wrong_lattice(self):
         user_input = '{"lattice": "bcc", "points": [["0.00", "0.00", "0.00"], ["1.00", "0.00", "0.00"], ["1.00", "1.00", "1.00"]]}'
         self.assertFalse(grade(user_input, {'miller': '(3,3,3)', 'lattice': 'fcc'}))