Removed Python (deprecated) and renamed Python_Simplified -> Python a…

…s per #68; updated README to reflect this change

CaseStudies/glass/Implementations/Python/Implementation/calculations.py

@@ -8,58 +8,16 @@
 
 from . import interp
 import math
-import numpy as np
 
-
-def calc_q(w_array, data_sd, data_q, params):
-    """
-    Interpolates q (demand from IM3) from Figure 2 using wtnt and SD
-    """
-    idx, jdx, kdx, num_interp1, num_interp2 = interp.find_bounds(w_array, data_sd, params.wtnt, params.sd)
-    q = interp.interp(idx, jdx, kdx, num_interp1, num_interp2, w_array, data_sd, data_q, params.wtnt, params.sd)
-    return q
-
-
-def calc_j(j_array, data_asprat, data_qstar, q, params):
-    """
-     Interpolates J (stress distribution factor from DD3) and q_hat_tol (tolerable pressure, DD7) from Figure 3 using
-     the aspect ratio, Jtol (DD8) and q_hat (DD6).
-     q* is interpolated for every J and stored in q_vec.  The corresponding J values are stored in j_vec.  J and
-     q_hat_tol are interpolated from q_vec and j_vec using q_hat and Jtol.
-    """
+
+def calc_q_hat(q, params):
     q_hat = q * pow((params.a * params.b), 2) / (params.E * pow(params.h, 4)) * (1 / params.gtf)
-    j_tol = math.log((math.log(1/(1-params.pbtol)))*((pow((params.a/1000)*(params.b/1000), params.m-1))/(params.k *
-                                                (pow(params.E*1000*(pow((params.h/1000), 2)), params.m))*params.ldf)))
-    q_vec = []
-    j_vec = []
-    for i in range(len(j_array)):
-        idx_2 = i
-        jdx_2 = (np.abs(data_asprat[:, idx_2] - params.asprat)).argmin()
-        if data_asprat[jdx_2, idx_2] > params.asprat:
-            jdx_2 -= 1
-        if data_asprat[0, idx_2] > params.asprat:
-            continue
-        else:
-            q_star = interp.interp(idx_2, jdx_2, -1, 0, 1, j_array, data_asprat, data_qstar, -1, params.asprat)
-            q_vec.append(q_star)
-            j_vec.append(j_array[idx_2])
-    q_vec = np.array(q_vec)
-    j_vec = np.array(j_vec)
-    idx_3 = (np.abs(q_vec - q_hat)).argmin()
-    if q_hat < np.amin(q_vec):
-        raise SystemExit("Error: q_hat(DD6 in the SRS) is out of bounds. q_hat is less than the smallest value in q_vec. The input a, b might be too small while t might be too large. Please refer to the data definitions section and the data constraints section in the SRS.")
-    if q_hat > np.amax(q_vec):
-        raise SystemExit("Error: q_hat(DD6 in SRS) is out of bounds. q_hat is greater than the biggest value in q_vec. The input a, b might be too large while t might be too small. Please refer to the data definitions section and the data constraints section in the SRS.")
-    else:
-        if q_vec[idx_3] > q_hat:
-            idx_3 -= 1
-        j = interp.lin_interp(j_vec[idx_3], j_vec[idx_3+1], q_vec[idx_3], q_vec[idx_3+1], q_hat)
-        idx_4 = (np.abs(j_vec - j_tol)).argmin()
-        if j_vec[idx_4] > j_tol:
-            idx_4 -= 1
-        q_hat_tol = interp.lin_interp(q_vec[idx_4], q_vec[idx_4+1], j_vec[idx_4], j_vec[idx_4+1], j_tol)
-        return j, q_hat_tol
-
+    return q_hat
+
+def calc_j_tol(params):
+    j_tol = math.log((math.log(1/(1-params.pbtol)))*((pow((params.a/1000)*(params.b/1000), params.m-1))/(params.k *(pow(params.E*1000*(pow((params.h/1000), 2)), params.m))*params.ldf)))
+    return j_tol                                         
+
 
 def calc_pb(j, params):
     """
@@ -71,30 +29,24 @@ def calc_pb(j, params):
     return pb
 
 
-def calc_lr(q_hat_tol, params):
-    """
-    Calculates Capacity (LR, IM2).
-    """
+def calc_nfl(q_hat_tol, params):
     nfl = (q_hat_tol * params.E * pow(params.h, 4)) / (pow(params.a * params.b, 2))
+    return nfl
+
+def calc_lr(nfl, params):
     lr = nfl * params.gtf * params.lsf
-    return lr, nfl
-
+    return lr
 
-def is_safe(pb, lr, q, params):
-    """
-    Checks if Pb < Pbtol and if LR > q (T1 and T2).
-    If both are true, the glass is considered safe.
-    """
+def is_safe1(pb, params):
     if pb < params.pbtol:
         is_safe1 = True
     else:
         is_safe1 = False
+    return is_safe1
+
+def is_safe2(lr, q):
     if lr > q:
         is_safe2 = True
     else:
         is_safe2 = False
-    if is_safe1 == True and is_safe2 == True:
-        safe = 'For the given input parameters, the glass is considered safe'
-    else:
-        safe = 'For the given input parameters, the glass is NOT considered safe'
-    return is_safe1, is_safe2, safe
+    return is_safe2

CaseStudies/glass/Implementations/Python/Implementation/checkConstraints.py

@@ -10,11 +10,11 @@ def check_constraints(params):
         raise SystemExit("InputError: a and b must be greater than 0")
     if params.asprat < 1 or params.asprat > 5:
         raise SystemExit("InputError: a/b must be between 1 and 5")
+    if not (params.t in ["2.5","2.7","3","3.0","4","4.0","5","5.0","6","6.0","8","8.0","10","10.0","12","12.0","16","16.0","19","19.0","22","22.0"]):
+        raise SystemExit("InputError: t must be in [2.5,2.7,3.0,4.0,5.0,6.0,8.0,10.0,12.0,16.0,19.0,22.0]")
     if params.tnt <= 0:
         raise SystemExit("InputError: TNT must be greater than 0")
     if params.wtnt < 4.5 or params.wtnt > 910:
         raise SystemExit("InputError: wtnt must be between 4.5 and 910")
     if params.sd<6 or params.sd>130:
         raise SystemExit("InputError: SD must be between 6 and 130")
-    if not (params.t in ["2.5","2.7","3.0","4.0","5.0","6.0","8.0","10.0","12.0","16.0","19.0","22.0"]):
-        raise SystemExit("InputError: t must be in [2.5,2.7,3.0,4.0,5.0,6.0,8.0,10.0,12.0,16.0,19.0,22.0]")

CaseStudies/glass/Implementations/Python/Implementation/inputFormat.py

@@ -5,8 +5,6 @@
 input parameters module.
 """
 
-from numpy import float64
-
 
 def get_input(filename, params):
     """
@@ -17,23 +15,23 @@ def get_input(filename, params):
     infile = open(filename, "r")
     text = infile.readline()
     # dimensions of glass slab
-    params.a = float64(infile.readline())
-    params.b = float64(infile.readline())
+    params.a = float(infile.readline())
+    params.b = float(infile.readline())
     params.t = infile.readline().rstrip()
     text = infile.readline()
     # glass type
     params.gt = infile.readline().rstrip()
     # Weight of Charge
     text = infile.readline()
-    params.w = float64(infile.readline())
+    params.w = float(infile.readline())
     text = infile.readline()
-    params.tnt = float64(infile.readline())
+    params.tnt = float(infile.readline())
     text = infile.readline()
     # stand off distance coordinates
-    sdx = float64(infile.readline())
-    sdy = float64(infile.readline())
-    sdz = float64(infile.readline())
+    sdx = float(infile.readline())
+    sdy = float(infile.readline())
+    sdz = float(infile.readline())
     params.sdvect = (sdx, sdy, sdz)
     text = infile.readline()
-    params.pbtol = float64(infile.readline())
+    params.pbtol = float(infile.readline())
     infile.close()

CaseStudies/glass/Implementations/Python/Implementation/interp.py

@@ -5,121 +5,75 @@
 data and return an interpolated value.
 """
 
-import numpy as np
 
+class BoundError(Exception):
+    pass
 
-def lin_interp(y1, y2, x1, x2, input_param):
-    """
-    Defines the equation for the linear interpolation.
-    """
-    y0 = y1 + (y2 - y1)/(x2 - x1)*(input_param - x1)
-    return y0
+def lin_interp(x1, y1, x2, y2, x):
+    y = (y2 - y1)/(x2 - x1)*(x - x1) + y1
+    return y
 
+# arr => Sj (i.e. sequence of numbers representing possible standOff distances)
+# v   => SD calculated (i.e. number)
+def indInSeq(arr, v):
+    for i in range(len(arr) - 1):
+        if arr[i] <= v and v <= arr[i+1]:
+            return i
+    raise BoundError
 
-def proper_index(index1,index2,data,value):
-    if index1 == 0:
-        if data[index1,index2] == data[index1+1,index2]:
-            index1 += 1
-    elif index1 == (data[:,index2]).argmax():
-        index1 -= 1
-    else:
-        if data[index1,index2] > value:
-            index1 -= 1
-        elif data[index1,index2] < value:
-            if (index1+1 < (data[:,index2]).argmax()) and (data[index1,index2] == data[index1+1,index2]):
-                index1 += 1
-            elif (index1+1 == (data[:,index2]).argmax()) and (data[index1,index2] == data[index1+1,index2]):
-                index1 -= 1 
-    return index1 
+def matrixCol(mat, c):
+    col = []
+    for i in range(len(mat)):
+        col.append(mat[i][c])
+    return col
+
+def interpY(x_array, y_array, z_array, x, z):
+    # find index that bounds z
+    i = indInSeq(z_array, z)
+
+    # x and y values for the z curves at i and i+1 (z_1 and z_2)
+    x_z_1 = matrixCol(x_array, i)
+    y_z_1 = matrixCol(y_array, i)
+    x_z_2 = matrixCol(x_array, i + 1)
+    y_z_2 = matrixCol(y_array, i + 1)
+
+    # find indices that bound x for each z curve
+    try:
+        j = indInSeq(x_z_1, x)
+        k = indInSeq(x_z_2, x)
+    except BoundError:
+        raise SystemExit("Interpolation of y failed")
+
+    # interpolate y values on z curves z_1 and z_2 at x
+    y_1 = lin_interp(x_z_1[j], y_z_1[j], x_z_1[j+1], y_z_1[j+1], x)
+    y_2 = lin_interp(x_z_2[k], y_z_2[k], x_z_2[k+1], y_z_2[k+1], x)
 
-
-# find_bounds: numpy.ndarray numpy.ndarray Num Num -> Int Int Int Int Int
-def find_bounds(data1, data2, value1, value2):
-    """
-    Finds the closest values above and below the input parameter to be used in
-    the interpolation.  Also returns how many interpolations are to be performed.
-    """
-    # idx represents the index of the closest value of value1 in the array data1(e.g. for w_array,idx=0 if input w = 4.5)
-    idx = (np.abs(data1 - value1)).argmin()
-    # if the closest value is greater than the input parameter value1, we let idx=idx-1 to ensure the input parameter is 
-    #   within the range [closest value below value1(data1[idx]), closest value above value1(data1[idx+1])]
-    if data1[idx] > value1:
-        idx -= 1
-    # jdx represents the index of the closest value of value2 in the ith column of the array data2(data2[:,idx])
-    jdx = (np.abs(data2[:, idx] - value2)).argmin()
-    # kdx represents the index of the closest value of value2 in the (i+1)th column of the array data2(data2[:,idx+1])
-    kdx = (np.abs(data2[:, idx+1] - value2)).argmin() 
-    # Case1: value1 in data1; num_interp1 = 0; look for value2 in data2[:,idx]
-    if value1 in data1:
-        num_interp1 = 0
-        # Case1_1: value2 in data2[:,idx]; num_interp2 = 0
-        if value2 in data2[:,idx]:
-            num_interp2 = 0
-        # Case1_2: value2 NOT in data2[:,idx]; num_interp2 = 1
-        else:
-            num_interp2 = 1
-    # Case2: value1 NOT in data1; num_interp1 = 1; look for value2 in both data2[:,idx] and data2[:,idx+1]   
-    else:
-        num_interp1 = 1
-        # Case2_1: value2 in both data2[:,idx] and data2[:,idx+1]; num_interp2 = 0
-        if value2 in data2[:,idx] and value2 in data2[:,idx+1]:
-            num_interp2 = 0
-        # Case2_2: value2 in data2[:,idx]; num_interp2 = 1
-        elif value2 in data2[:,idx]:
-            num_interp2 = 1
-        # Case2_3: value2 in data2[:,idx+1]; num_interp2 = 2
-        elif value2 in data2[:, idx+1]:
-            num_interp2 = 2
-        # Case2_4: value2 NOT in data2[:,idx] or data2[:,idx+1]; num_interp2 = 3
-        else:
-            num_interp2 = 3
-    jdx = proper_index(jdx,idx,data2,value2)
-    kdx = proper_index(kdx,idx+1,data2,value2)
-    return idx, jdx, kdx, num_interp1, num_interp2
-
-
-# interp: Int Int Int Int Int Int numpy.ndarray numpy.ndarray numpy.ndarray Num Num -> Num
-def interp(idx, jdx, kdx, num_interp1, num_interp2, data1, data2, data3, value1, value2):
-    """
-    Performs the appropriate number of interpolations based on the output of find_bounds.
-    """
-    # Case1_1: value1 in data1 and value2 in data2[:,idx]; no need for interpolation, interp_value = data3[jdx,idx]
-    if num_interp1 == 0 and num_interp2 == 0:
-        interp_value = data3[jdx, idx]
-    # Case1_2: value1 in data1 but value2 NOT in data2[:,idx]; use the closest values above and below value2 in data2[:,idx] to do linear interpolation, 
-    #   where (x1,y1)=(data2[jdx,idx],data3[jdx,idx]), (x2,y2)=(data2[jdx+1,idx],data3[jdx+1,idx])
-    elif num_interp1 == 0 and num_interp2 == 1:
-        interp_value = lin_interp(data3[jdx, idx], data3[jdx+1, idx], data2[jdx, idx], data2[jdx+1, idx], value2)
-    # Case2_1: value1 NOT in data1 but value2 in both data2[:,idx] and data2[:,idx+1]; use the closest values above and below value1 to do linear
-    #   interpolation, where (x1,y1)=(data1[idx],data3[jdx,idx]), (x2,y2)=(data1[idx+1],data3[kdx,idx+1]) 
-    elif num_interp1 == 1 and num_interp2 == 0:
-        y0_1 = data3[jdx, idx]
-        y0_2 = data3[kdx, idx+1]
-        interp_value = lin_interp(y0_1, y0_2, data1[idx], data1[idx+1], value1)
-    # Case2_2: value1 NOT in data1 but value2 in data2[:,idx]; use the closest values below and above value2 in data2[:,idx+1] to do linear interpolation for 
-    #   y2(x2=data1[idx+1] but value2 not in data2[:,idx+1]); then use the closest values below and above value1 to do linear interpolation for interp_value, 
-    #   where (x1,y1)=(data1[idx],data3[jdx,idx]), (x2,y2)=(data1[idx+1],lin_interp(data3[kdx,idx+1],data3[kdx+1,idx+1],data2[kdx,idx+1],data2[kdx+1,idx+1],value2))
-    elif num_interp1 == 1 and num_interp2 == 1:
-        y0_1 = data3[jdx, idx]
-        y0_2 = lin_interp(data3[kdx, idx+1], data3[kdx+1, idx+1], data2[kdx, idx+1], data2[kdx+1, idx+1], value2)
-        interp_value = lin_interp(y0_1, y0_2, data1[idx], data1[idx+1], value1)
-    # Case2_3: value1 NOT in data1 but value2 in data2[:,idx+1]; use the closest values below and above value2 in data2[:,idx] to do linear interpolation for
-    #   y1(x1=data1[idx] but value2 not in data2[:,idx]); then use the closest values below and above value1 to do linear interpolation for interp_value, where
-    #   (x1,y1)=(data1[idx],lin_interp(data3[jdx,idx],data3[jdx+1,idx],data2[jdx,idx],data2[jdx+1.idx],value2)), (x2,y2)=(data1[idx+1], data3[kdx,idx+1])
-    elif num_interp1 == 1 and num_interp2 ==2:
-        y0_2 = data3[kdx, idx+1]
-        y0_1 = lin_interp(data3[jdx, idx], data3[jdx+1, idx], data2[jdx, idx], data2[jdx+1, idx], value2)
-        x1_1 = data1[idx]
-        x1_2 = data1[idx+1]
-        interp_value = lin_interp(y0_1, y0_2, x1_1, x1_2, value1)
-    # Case2_4: value1 NOT in data1 and value2 NOT in data2[:,idx] or data2[:,idx+1]; use the closest values below and above value2 in data2[:,idx] to do 
-    #   linear interpolation for y1 and the closest values below and above value2 in data2[:,idx+1] to do linear interpolation for y2; then use the closest
-    #   values of value1 in data1 to do linear interpolation for interp_value, where (x1,y1)=(data1[idx],lin_interp(data3[jdx,idx],data3[jdx+1,idx],
-    #   data2[jdx,idx],data2[jdx+1,idx],value2)), (x2,y2)=(data1[idx+1],lin_interp(data3[kdx,idx+1],data3[kdx+1,idx+1],data2[kdx,idx+1],data2[kdx+1,idx+1],value2))
-    elif num_interp1 == 1 and num_interp2 == 3:
-        y0_1 = lin_interp(data3[jdx, idx], data3[jdx+1, idx], data2[jdx, idx], data2[jdx+1, idx], value2)
-        y0_2 = lin_interp(data3[kdx, idx+1], data3[kdx+1, idx+1], data2[kdx, idx+1], data2[kdx+1, idx+1], value2)
-        x0_1 = data1[idx]
-        x0_2 = data1[idx+1]
-        interp_value = lin_interp(y0_1, y0_2, x0_1, x0_2, value1)
-    return interp_value
+    # interpolate y
+    return lin_interp(z_array[i], y_1, z_array[i+1], y_2, z)
+
+
+def interpZ(x_array, y_array, z_array, x, y):
+    for i in range(len(z_array) - 1):
+        # x and y values for two consecutive z curves 
+        x_z_1 = matrixCol(x_array, i)
+        y_z_1 = matrixCol(y_array, i)
+        x_z_2 = matrixCol(x_array, i + 1)
+        y_z_2 = matrixCol(y_array, i + 1)
+
+        # find indices that bound x for each z curve
+        try:
+            j = indInSeq(x_z_1, x)
+            k = indInSeq(x_z_2, x)
+        except BoundError:
+            continue
+
+        # interpolate y values at x on z curves
+        y_lower = lin_interp(x_z_1[j], y_z_1[j], x_z_1[j+1], y_z_1[j+1], x)
+        y_upper = lin_interp(x_z_2[k], y_z_2[k], x_z_2[k+1], y_z_2[k+1], x)
+
+        # check if y is bound
+        if y_lower <= y and y <= y_upper:
+            # interpolate z
+            return lin_interp(y_lower, z_array[i], y_upper, z_array[i+1], y)
+
+    raise SystemExit("Interpolation of z failed")

CaseStudies/glass/Implementations/Python/Implementation/mainfun.py

@@ -12,24 +12,34 @@
 from . import calculations
 from . import derivedValues
 from . import checkConstraints
+from . import interp
 
 
 def main(filename):
     params = param.Param()
     inputFormat.get_input(filename, params)
     derivedValues.derived_params(params)
-
     checkConstraints.check_constraints(params)
-    w_array, data_sd, data_q = readTable.read_table('TSD.txt')
-    j_array, data_asprat, data_qstar = readTable.read_table('SDF.txt')
+
+    w_array = readTable.read_z_array('TSD.txt')
+    data_sd = readTable.read_x_array('TSD.txt', len(w_array))
+    data_q = readTable.read_y_array('TSD.txt', len(w_array))
+
+    j_array = readTable.read_z_array('SDF.txt')
+    data_asprat = readTable.read_x_array('SDF.txt', len(j_array))
+    data_qstar = readTable.read_y_array('SDF.txt', len(j_array))
 
-    q = calculations.calc_q(w_array, data_sd, data_q, params)
-    j, q_hat_tol = calculations.calc_j(j_array, data_asprat, data_qstar, q, params)
+    q = interp.interpY(data_sd, data_q, w_array, params.sd, params.wtnt)
+    q_hat = calculations.calc_q_hat(q, params)
+    j_tol = calculations.calc_j_tol(params)
+    j = interp.interpZ(data_asprat, data_qstar, j_array, params.asprat, q_hat)
+    q_hat_tol = interp.interpY(data_asprat, data_qstar, j_array, params.asprat, j_tol)
     pb = calculations.calc_pb(j, params)
-    lr, nfl = calculations.calc_lr(q_hat_tol, params)
-    is_safe1, is_safe2, safe = calculations.is_safe(pb, lr, q, params)
-
-    outputFormat.display_output('outputfile.txt', q, j, q_hat_tol, pb, lr, nfl, is_safe1, is_safe2, safe, params)
+    nfl = calculations.calc_nfl(q_hat_tol, params)
+    lr = calculations.calc_lr(nfl, params)
+    is_safe1 = calculations.is_safe1(pb, params)
+    is_safe2 = calculations.is_safe2(lr, q)
+    outputFormat.display_output('outputfile.txt', q, j, q_hat_tol, pb, lr, nfl, is_safe1, is_safe2, params)
     print("Main has been executed and the results have been written to 'outputfile.txt'.")
 
 if __name__ == "__main__":