diff --git a/common/lib/capa/capa/chem/miller.py b/common/lib/capa/capa/chem/miller.py index 8a98732874..f397e9d35e 100644 --- a/common/lib/capa/capa/chem/miller.py +++ b/common/lib/capa/capa/chem/miller.py @@ -1,3 +1,5 @@ +""" 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: + 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. + + 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. + + 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]) - print 'common_denominator:', common_denominator - - # 2) lead to a common denominator - # 3) throw away denominator - miller = [fract.numerator * math.fabs(common_denominator) / fract.denominator for fract in fracts] - - # 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 + 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. + + Algorithm: + + 1. Calculate normal vector to a plane that goes trough all points. + + 2. Set origin. + + 3. Create Cartesian coordinate system (Ccs). + + 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) """ - # 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 - - # origin O = np.array([0, 0, 0]) - - # 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]), - ] - + 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'}))