text:
"""ATORtf_RFtri.py

# Author:
Richard Bruce Baxter - Copyright (c) 2021 Baxter AI (baxterai.com)

# License:
MIT License

# Installation:
See ATORtf.tf

# Usage:
See ATORtf.tf

# Description:
ATORtf RF Tri - generate tri (as represented by 3 feature points) receptive fields

"""

import tensorflow as tf
import numpy as np
import cv2
import copy

import ATORtf_RFproperties
import ATORtf_operations

#RF Tri can be applied to luminosity or contrast maps

debugSmallIterations = False

pointFeatureRFinsideRadius = 0.5
pointFeatureRFopponencyAreaFactor = 2	#floats unsupported by opencv ellipse draw

generatePointFeatureCorners = True
if(generatePointFeatureCorners):
	minimumCornerOpponencyPosition = -1	#CHECKTHIS
	maximumCornerOpponencyPosition = 1		#CHECKTHIS
	cornerOpponencyPositionResolution = 1	#CHECKTHIS
else:
	#only generate simple point features;
	minimumCornerOpponencyPosition = 0
	maximumCornerOpponencyPosition = 0
	cornerOpponencyPositionResolution = 1

#match ATORtf_RFellipse algorithm;
matchRFellipseAlgorithm = False

#match ATORtf_RFellipse algorithm;
RFnormaliseLocalEquilateralTriangle = True

#match ATORtf_RFellipse algorithm;
ellipseNormalisedAngle = 0.0
ellipseNormalisedCentreCoordinates = 0.0
ellipseNormalisedAxesLength = 1.0

#match ATORtf_RFellipse algorithm;
# ellipse axesLength definition (based on https://docs.opencv.org/4.5.3/d6/d6e/group__imgproc__draw.html#ga28b2267d35786f5f890ca167236cbc69)
#  ___
# / | \ minimumEllipseAxisLength1 | (axis1: axes width)
# | --| minimumEllipseAxisLength2 - (axis2: axis height)
# \___/
#
if(matchRFellipseAlgorithm):
	receptiveFieldOpponencyAreaFactorTri = ATORtf_RFproperties.receptiveFieldOpponencyAreaFactorEllipse	#2.0
	maximumAxisLengthMultiplierTri = 1
	maximumAxisLengthMultiplier = ATORtf_RFproperties.maximumAxisLengthMultiplierDefault
	minimumEllipseAxisLength1 = ATORtf_RFproperties.minimumEllipseAxisLength*2	#minimum elongation is required	#can be set to 1 as (axesLength1 > axesLength2) condition is enforced for RFellipse creation
	minimumEllipseAxisLength2 = ATORtf_RFproperties.minimumEllipseAxisLength
else:
	receptiveFieldOpponencyAreaFactorTri = 1.0
	maximumAxisLengthMultiplierTri = 2	#tri RF support 2x larger coordinates than ellipse RF
	maximumAxisLengthMultiplier = ATORtf_RFproperties.maximumAxisLengthMultiplierDefault*maximumAxisLengthMultiplierTri
	minimumEllipseAxisLength1 = ATORtf_RFproperties.minimumEllipseAxisLength*2*maximumAxisLengthMultiplierTri	#4	#minimum elongation is required	#can be set to 1 as (axesLength1 > axesLength2) condition is enforced for RFellipse creation
	minimumEllipseAxisLength2 = ATORtf_RFproperties.minimumEllipseAxisLength*maximumAxisLengthMultiplierTri	#2

if(ATORtf_RFproperties.lowResultFilterPosition):
	ellipseAxesLengthResolution = 1*maximumAxisLengthMultiplierTri	#2*maximumAxisLengthMultiplierTri	#2*	#pixels (at resolution r)	#use course grain resolution to decrease number of filters #OLD: 1
	if(debugSmallIterations):
		ellipseAxesLengthResolution = ellipseAxesLengthResolution*2	#increase resolution higher than is allowed to be drawn
else:
	ellipseAxesLengthResolution = 1*maximumAxisLengthMultiplierTri	#pixels (at resolution r)
ellipseAngleResolution = 30	#degrees
ellipseColourResolution = 64	#bits

def printTriProperties(triProperties):
	print("vertexCoordinatesRelative = ", triProperties.vertexCoordinatesRelative)

#~match ATORtf_RFellipse algorithm;
def normaliseGlobalTriProperties(ellipseProperties, resolutionFactor):
	resolutionFactor = ellipseProperties.resolutionFactor
	ellipsePropertiesNormalised = copy.deepcopy(ellipseProperties)
	ellipsePropertiesNormalised.centerCoordinates = (ellipsePropertiesNormalised.centerCoordinates[0]*resolutionFactor, ellipsePropertiesNormalised.centerCoordinates[1]*resolutionFactor)
	ellipsePropertiesNormalised.axesLength = (ellipsePropertiesNormalised.axesLength[0]*resolutionFactor, ellipsePropertiesNormalised.axesLength[1]*resolutionFactor)
	return ellipsePropertiesNormalised

#match ATORtf_RFellipse algorithm;
def normaliseLocalTriProperties(RFproperties):
	#normalise ellipse respect to major/minor ellipticity axis orientation (WRT self)
	RFpropertiesNormalised = copy.deepcopy(RFproperties)
	RFpropertiesNormalised.angle = ellipseNormalisedAngle	#CHECKTHIS
	RFpropertiesNormalised.centerCoordinates = (ellipseNormalisedCentreCoordinates, ellipseNormalisedCentreCoordinates, ellipseNormalisedCentreCoordinates)
	if(RFnormaliseLocalEquilateralTriangle):
		RFpropertiesNormalised.axesLength = ATORtf_operations.getEquilateralTriangleAxesLength(ellipseNormalisedAxesLength)
	else:
		RFpropertiesNormalised.axesLength = (ellipseNormalisedAxesLength, ellipseNormalisedAxesLength)

	return RFpropertiesNormalised

#match ATORtf_RFellipse algorithm;
def calculateFilterPixels(filterSize, numberOfDimensions):
	internalFilterSize = getInternalFilterSize(filterSize, numberOfDimensions)	#CHECKTHIS: only consider contribution of positive (additive) pixels
	#print("internalFilterSize = ", internalFilterSize)

	if(numberOfDimensions == 2):
		numberOfFilterPixels = internalFilterSize[0]*internalFilterSize[1] * 3	#3 features
	elif(numberOfDimensions == 3):
		numberOfFilterPixels = internalFilterSize[0]*internalFilterSize[1] * 3	#3 features	#CHECKTHIS
	#print("numberOfFilterPixels = ", numberOfFilterPixels)
	return numberOfFilterPixels

#match ATORtf_RFellipse algorithm;
def getInternalFilterSize(filterSize, numberOfDimensions):
	if(numberOfDimensions == 2):
		internalFilterSize = (int(pointFeatureRFinsideRadius*pointFeatureRFopponencyAreaFactor), int(pointFeatureRFinsideRadius*pointFeatureRFopponencyAreaFactor))
	elif(numberOfDimensions == 3):
		internalFilterSize = (int(pointFeatureRFinsideRadius*pointFeatureRFopponencyAreaFactor), int(pointFeatureRFinsideRadius*pointFeatureRFopponencyAreaFactor))	#CHECKTHIS
	return internalFilterSize

def generateRFfiltersTri(resolutionProperties, RFfiltersList, RFfiltersPropertiesList):
	#generate filter types

	#2D code;

	#filters are generated based on human magnocellular/parvocellular/koniocellular wavelength discrimination in LGN and VX (double/opponent receptive fields)
	filterTypeIndex = 0

	#magnocellular filters (monochromatic);
	colourH = (255, 255, 255)	#high
	colourL = (-255, -255, -255)	#low
	RFfiltersHL, RFpropertiesHL = generateRotationalInvariantRFfilters(resolutionProperties, False, colourH, colourL, filterTypeIndex)
	filterTypeIndex+=1

	RFfiltersList.append(RFfiltersHL)

	RFfiltersPropertiesList.append(RFpropertiesHL)

def generateRotationalInvariantRFfilters(resolutionProperties, isColourFilter, filterInsideColour, filterOutsideColour, filterTypeIndex):

	RFfiltersList2 = []
	RFfiltersPropertiesList2 = []

	#FUTURE: consider storing filters in n dimensional array and finding local minima of filter matches across all dimensions

	#reduce max size of ellipse at each res
	axesLengthMax, filterRadius, filterSize = ATORtf_RFproperties.getFilterDimensions(resolutionProperties, maximumAxisLengthMultiplier, receptiveFieldOpponencyAreaFactorTri)

	#print("axesLengthMax = ", axesLengthMax)

	for axesLength1 in range(minimumEllipseAxisLength1, axesLengthMax[0]+1, ellipseAxesLengthResolution):
		for axesLength2 in range(minimumEllipseAxisLength2, axesLengthMax[1]+1, ellipseAxesLengthResolution):
			if(axesLength1 > axesLength2):	#ensure that ellipse is always alongated towards axis1 (required for consistent normalisation)
				if((axesLength1 < filterRadius) and (axesLength2 < filterRadius)):
					for angle in range(0, 360, ellipseAngleResolution):	#degrees
						for corner1OpponencyPosition1 in range(minimumCornerOpponencyPosition, maximumCornerOpponencyPosition+1, cornerOpponencyPositionResolution):
							for corner1OpponencyPosition2 in range(minimumCornerOpponencyPosition, maximumCornerOpponencyPosition+1, cornerOpponencyPositionResolution):
								for corner2OpponencyPosition1 in range(minimumCornerOpponencyPosition, maximumCornerOpponencyPosition+1, cornerOpponencyPositionResolution):
									for corner2OpponencyPosition2 in range(minimumCornerOpponencyPosition, maximumCornerOpponencyPosition+1, cornerOpponencyPositionResolution):
										for corner3OpponencyPosition1 in range(minimumCornerOpponencyPosition, maximumCornerOpponencyPosition+1, cornerOpponencyPositionResolution):
											for corner3OpponencyPosition2 in range(minimumCornerOpponencyPosition, maximumCornerOpponencyPosition+1, cornerOpponencyPositionResolution):

												axesLength = (axesLength1, axesLength2)

												axesLengthInside = (pointFeatureRFinsideRadius, pointFeatureRFinsideRadius)	#estimation
												axesLengthOutside = (int(pointFeatureRFinsideRadius*pointFeatureRFopponencyAreaFactor), int(pointFeatureRFinsideRadius*pointFeatureRFopponencyAreaFactor))	#estimation
												angleInside = 0.0	#circular
												angleOutside = 0.0	#circular
												filterCenterCoordinates = (0, 0)

												RFpropertiesInside = ATORtf_RFproperties.RFpropertiesClass(resolutionProperties.resolutionIndex, resolutionProperties.resolutionFactor, filterSize, ATORtf_RFproperties.RFtypeEllipse, filterCenterCoordinates, axesLengthInside, angleInside, filterInsideColour)
												RFpropertiesOutside = ATORtf_RFproperties.RFpropertiesClass(resolutionProperties.resolutionIndex, resolutionProperties.resolutionFactor, filterSize, ATORtf_RFproperties.RFtypeEllipse, filterCenterCoordinates, axesLengthOutside, angleOutside, filterOutsideColour)
												RFpropertiesInside.isColourFilter = isColourFilter
												RFpropertiesOutside.isColourFilter = isColourFilter

												#number of corner filters;
												#3^6 = 729 filters (or 9^3)

												#example corner point filter;
												#++-0
												#++--
												#----
												#0--0

												vertexCoordinatesRelative = ATORtf_RFproperties.deriveTriVertexCoordinatesFromArtificialEllipseProperties(axesLength, angle)
												#print("vertexCoordinatesRelative = ", vertexCoordinatesRelative)
												vertexCoordinatesRelativeInside = copy.deepcopy(vertexCoordinatesRelative)
												vertexCoordinatesRelativeOutside = copy.deepcopy(vertexCoordinatesRelative)
												vertexCoordinatesRelativeOutside[0][0] = vertexCoordinatesRelativeOutside[0][0] + corner1OpponencyPosition1
												vertexCoordinatesRelativeOutside[0][1] = vertexCoordinatesRelativeOutside[0][1] + corner1OpponencyPosition2
												vertexCoordinatesRelativeOutside[1][0] = vertexCoordinatesRelativeOutside[1][0] + corner2OpponencyPosition1
												vertexCoordinatesRelativeOutside[1][1] = vertexCoordinatesRelativeOutside[1][1] + corner2OpponencyPosition2
												vertexCoordinatesRelativeOutside[2][0] = vertexCoordinatesRelativeOutside[2][0] + corner3OpponencyPosition1
												vertexCoordinatesRelativeOutside[2][1] = vertexCoordinatesRelativeOutside[2][1] + corner3OpponencyPosition2

												#vertexCoordinatesRelativeInside[0][0] = vertexCoordinatesRelativeInside[0][0] + corner1OpponencyPosition1
												#vertexCoordinatesRelativeInside[0][1] = vertexCoordinatesRelativeInside[0][1] + corner1OpponencyPosition2
												#vertexCoordinatesRelativeInside[1][0] = vertexCoordinatesRelativeInside[1][0] + corner2OpponencyPosition1
												#vertexCoordinatesRelativeInside[1][1] = vertexCoordinatesRelativeInside[1][1] + corner2OpponencyPosition2
												#vertexCoordinatesRelativeInside[2][0] = vertexCoordinatesRelativeInside[2][0] + corner3OpponencyPosition1
												#vertexCoordinatesRelativeInside[2][1] = vertexCoordinatesRelativeInside[2][1] + corner3OpponencyPosition2

												RFfilter = generateRFfilter(resolutionProperties.resolutionIndex, isColourFilter, RFpropertiesInside, RFpropertiesOutside, vertexCoordinatesRelativeInside, vertexCoordinatesRelativeOutside)
												RFfiltersList2.append(RFfilter)

												RFproperties = ATORtf_RFproperties.RFpropertiesClass(resolutionProperties.resolutionIndex, resolutionProperties.resolutionFactor, filterSize, ATORtf_RFproperties.RFtypeTri, filterCenterCoordinates, axesLength, angle, filterInsideColour)
												#RFproperties.centerCoordinates = centerCoordinates 	#centerCoordinates are set after filter is applied to imageSegment
												RFproperties.isColourFilter = isColourFilter
												RFfiltersPropertiesList2.append(RFproperties)	#CHECKTHIS: use RFpropertiesInside not RFpropertiesOutside

												#debug:
												#print(RFfilter.shape)
												if(resolutionProperties.debugVerbose):
													ATORtf_RFproperties.printRFproperties(RFproperties)
												#print("RFfilter = ", RFfilter)

												RFfilterImageFilename = "RFfilterResolutionIndex" + str(resolutionProperties.resolutionIndex) + "filterTypeIndex" + str(filterTypeIndex) + "axesLength1" + str(axesLength1) + "axesLength2" + str(axesLength2) + "angle" + str(angle) + "corner1OpponencyPosition1" + str(corner1OpponencyPosition1) + "corner1OpponencyPosition2" + str(corner1OpponencyPosition2) + "corner2OpponencyPosition1" + str(corner2OpponencyPosition1) + "corner2OpponencyPosition2" + str(corner2OpponencyPosition2) + "corner3OpponencyPosition1" + str(corner3OpponencyPosition1) + "corner3OpponencyPosition2" + str(corner3OpponencyPosition2) + ".png"
												ATORtf_RFproperties.saveRFFilterImage(RFfilter, RFfilterImageFilename)

	#create 3D tensor (for hardware accelerated test/application of filters)
	RFfiltersTensor = tf.stack(RFfiltersList2, axis=0)

	return RFfiltersTensor, RFfiltersPropertiesList2

def generateRFfilter(resolutionProperties, isColourFilter, RFpropertiesInside, RFpropertiesOutside, vertexCoordinatesRelativeInside, vertexCoordinatesRelativeOutside):

	# RF filter example (RFfilterTF) - generatePointFeatureCorners=False:
	#
	# 0 0 0 0 0 0 0 0 0 0 0 0
	# 0 0 - 0 0 0 0 - 0 0 0 0
	# 0 - + - 0 0 - + - 0 0 0
	# 0 0 - 0 0 0 0 - 0 0 0 0
	# 0 0 0 0 0 0 0 0 0 0 0 0
	# 0 0 0 0 0 0 0 0 0 0 0 0
	# 0 0 0 0 - 0 0 0 0 0 0 0
	# 0 0 0 - + - 0 0 0 0 0 0
	# 0 0 0 0 - 0 0 0 0 0 0 0
	# 0 0 0 0 0 0 0 0 0 0 0 0
	#
	# where "-" = -RFcolourOutside [R G B], "+" = +RFcolourInside [R G B], and "0" = [0, 0, 0]

	#generate ellipse on blank canvas
	blankArray = np.full((RFpropertiesInside.imageSize[1], RFpropertiesInside.imageSize[0], ATORtf_operations.rgbNumChannels), 0, np.uint8)	#rgb
	RFfilterTF = tf.convert_to_tensor(blankArray, dtype=tf.float32)

	RFpropertiesInside1 = copy.deepcopy(RFpropertiesInside)
	RFpropertiesOutside1 = copy.deepcopy(RFpropertiesOutside)
	RFpropertiesInside2 = copy.deepcopy(RFpropertiesInside)
	RFpropertiesOutside2 = copy.deepcopy(RFpropertiesOutside)
	RFpropertiesInside3 = copy.deepcopy(RFpropertiesInside)
	RFpropertiesOutside3 = copy.deepcopy(RFpropertiesOutside)

	RFpropertiesInside1.centerCoordinates = (RFpropertiesInside1.centerCoordinates[0]+vertexCoordinatesRelativeInside[0][0], RFpropertiesInside1.centerCoordinates[1]+vertexCoordinatesRelativeInside[0][1])
	RFpropertiesOutside1.centerCoordinates = (RFpropertiesOutside1.centerCoordinates[0]+vertexCoordinatesRelativeOutside[0][0], RFpropertiesOutside1.centerCoordinates[1]+vertexCoordinatesRelativeOutside[0][1])
	RFpropertiesInside2.centerCoordinates = (RFpropertiesInside2.centerCoordinates[0]+vertexCoordinatesRelativeInside[1][0], RFpropertiesInside2.centerCoordinates[1]+vertexCoordinatesRelativeInside[1][1])
	RFpropertiesOutside2.centerCoordinates = (RFpropertiesOutside2.centerCoordinates[0]+vertexCoordinatesRelativeOutside[1][0], RFpropertiesOutside2.centerCoordinates[1]+vertexCoordinatesRelativeOutside[1][1])
	RFpropertiesInside3.centerCoordinates = (RFpropertiesInside3.centerCoordinates[0]+vertexCoordinatesRelativeInside[2][0], RFpropertiesInside3.centerCoordinates[1]+vertexCoordinatesRelativeInside[2][1])
	RFpropertiesOutside3.centerCoordinates = (RFpropertiesOutside3.centerCoordinates[0]+vertexCoordinatesRelativeOutside[2][0], RFpropertiesOutside3.centerCoordinates[1]+vertexCoordinatesRelativeOutside[2][1])

	#print("RFpropertiesInside1.centerCoordinates = ", RFpropertiesInside1.centerCoordinates)
	#print("RFpropertiesOutside1.centerCoordinates = ", RFpropertiesOutside1.centerCoordinates)

	if(generatePointFeatureCorners):
		drawFeatureType = ATORtf_RFproperties.RFfeatureTypeCorner
	else:
		drawFeatureType = ATORtf_RFproperties.RFfeatureTypePoint

	RFfilterTF = ATORtf_RFproperties.drawRF(RFfilterTF, RFpropertiesInside1, RFpropertiesOutside1, drawFeatureType, True)
	RFfilterTF = ATORtf_RFproperties.drawRF(RFfilterTF, RFpropertiesInside2, RFpropertiesOutside2, drawFeatureType, True)
	RFfilterTF = ATORtf_RFproperties.drawRF(RFfilterTF, RFpropertiesInside3, RFpropertiesOutside3, drawFeatureType, True)

	#print("RFfilterTF = ", RFfilterTF)

	if(ATORtf_operations.storeRFfiltersValuesAsFractions):
		RFfilterTF = tf.divide(RFfilterTF, ATORtf_operations.rgbMaxValue)

	if(not isColourFilter):
		RFfilterTF = tf.image.rgb_to_grayscale(RFfilterTF)

	#print("RFfilterTF.shape = ", RFfilterTF.shape)
	#print("RFfilterTF = ", RFfilterTF)

	return RFfilterTF

max_stars_repo_name: baxterai/ATORtf

max_issues_repo_name: baxterai/ATORtf

max_forks_repo_name: baxterai/ATORtf


text:
import numpy as np
import matplotlib.pyplot as plt

def main():
    temp = [90, 95, 100, 105, 110, 115, 120, 125]
    x = [e*e for e in temp]
    y = [-59.86, -56.57, -52.55, -50.54, -46.88, -43.23, -39.75, -35.55]

    y = [np.sin(np.radians(e)) for e in y]

    # fig, axs = plt.subplots(nrows=1, ncols=1)
    # fig.suptitle('C3 Constant - Angle Vs Input Signal')
    plt.title('C3 Constant - Angle Vs Input Signal -> 3.696e-05 x - 1.17')
    plt.plot(x,y)
    plt.xlabel('Input Signal')
    plt.ylabel('Angle')

    z = np.poly1d(np.polyfit(x,y,1))
    y_p = z(x)

    # slope = ((x*y_p).mean(axis=1) - x.mean()*y_p.mean(axis=1)) / ((x**2).mean() - (x.mean())**2)

    print(z)
    plt.plot(x,y_p)

    plt.show()

if __name__ == '__main__':
    main()

ext: py

reason: import numpy


text:
import numpy as np
from numpy import pi
from collections import namedtuple

def return_muscle_settings(PreselectedMuscles=None):
	"""
	Notes:
	Coefficients from observation, Ramsay; 2009, FVC, Holtzbaur, Pigeon, Kuechle, or Banks. Optimal Muscle Length given in mm. Optimal tendon/muscle lengths and PCSA were taken from Garner and Pandy (2003)
	"""

	# Coefficients from observation, Ramsay; 2009, Pigeon, FVC, Holtzbaur, Garner & Pandy, or Banks.

	MA_Settings = namedtuple("MA_Settings",["Values","Source","Units","Equation_Number","Threshold","DOF"])
	Spindle_Settings = namedtuple("Spindle_Settings",["ActualNumber",'CorrectedNumber','RelativeAbundance',"Source"])
	Input_Source = namedtuple("Source_Settings",["Values","Source","Units"])

	def Pigeon_coeff_conversion(Coefficients):
		"""
		Takes in Coefficient values from Pigeon (1996) -- which take in angles in degrees -- and coverts them into the properly scaled coefficients for radians, additionally scaled by the magnitude listed in the paper.

		Note that the coefficients listed in Pigeon (1996) are given in decending order (i.e., c₅,c₄,c₃,c₂,c₁,c₀). However to maintain continuity with the equations given in Ramsay; 2009 (2009), we list coefficients in order of increasing power (i.e., c₀,c₁,c₂,c₃,c₄,c₅).
		"""
		import numpy as np
		assert len(Coefficients)==6, 'For Pigeon (1996) the list of Coefficients must be 6 elements long. Insert zeros (0) for any additional empty coefficients.'
		assert type(Coefficients)==list, 'Coefficients must be a 6 element list.'
		Rad_Conversion = np.multiply(Coefficients,\
				np.array([1,(180/np.pi),(180/np.pi)**2,(180/np.pi)**3,(180/np.pi)**4,(180/np.pi)**5],dtype = 'float64'))
		new_Coefficients =\
			np.multiply(Rad_Conversion,np.array([1,1e-1,1e-3,1e-5,1e-7,1e-9],dtype='float64'))
		return(new_Coefficients)

	PC_Settings = {\
		'Notes' : [\
						'This is the *clavicular* portion of the pectoralis major.',\
						'Banks and Garner & Pandy are parameter values for the entire muscle. Pigeon and Holzbaur have the values for the clavicular portion only.',\
						'Holzbaur parameters for shoulder are for frontal plane (ABD/ADD) only! This explains the relatively small values for some measures.'\
					],\
		'Shoulder MA' : {	"Primary Source" : "Pigeon; 1996",\
		 					"Sources" : \
								[\
									MA_Settings([50.80,0,0,0,0,0], 'Pigeon; 1996', "mm", None, None, "Shoulder"),\
									MA_Settings([2,0,0,0,0,0], 'Holzbaur; 2005', "mm", None, None, "Shoulder")\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Est", \
		 				"Sources" : \
							[\
								MA_Settings(0, "Est", "m", None, None, 'Elbow')\
							]}, \
		'Spindle' : Spindle_Settings(450,389.7,1.2,"Banks; 2006"),\
		'Mass' : {	"Primary Source" : "Banks; 2006", \
					"Sources" : \
						[\
							Input_Source(295.6, "Banks; 2006", 'g')\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(14.4, 'Holzbaur; 2005', 'cm'),\
											Input_Source(150, 'Est', 'mm')\
										]}, \
		'Optimal Tendon Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(0.3, 'Holzbaur; 2005', 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(17, 'Holzbaur; 2005', 'degrees')\
									]}, \
		'PCSA' : {	"Primary Source" : "Garner & Pandy; 2003", \
					'Sources' : \
						[\
							Input_Source(36.20,'Garner & Pandy; 2003','sq cm'),\
							Input_Source(2.6,'Holzbaur; 2005','sq cm')\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Garner & Pandy; 2003", \
										'Sources': \
											[\
												Input_Source(1175.01,'Garner & Pandy; 2003','N'),\
												Input_Source(364.4,'Holzbaur; 2005','N')\
											]}\
		}

	DELTa_Settings = {\
		'Notes' : [\
					"SFE MA is listed as 33.02 mm in Pigeon and estimated as 19 mm. Using Pigeon Coefficients convention, Kuechle (1997) has the DELTp MA for [-140,90] as Pigeon_coeff_conversion([ 13.4293189,  2.0316226, -0.2339031,  2.7807828,  0.,  0.]). This will yield a piecewise function that creates jumps with the new velocity formulation. Instead, we are going to try Pigeon_coeff_conversion([ 12.7928795,  2.0480346,  0.8917734,  3.2207214, -2.3928223,  0.]) so that the function is within range and continuous during the ROM. Threshold (pi/2) has been removed for this new MA function.",\
					"Garner & Pandy have much larger PCSA and Peak Force Values but only consider the entire Deltoid.",\
					"Holzbaur parameters for shoulder are for frontal plane (ABD/ADD) only! This explains the relatively small values for some measures.",\
					"Banks only had mass and spindle settings for the entire deltoid. As a result the parameters are divided by 3 as an estimate of the individual muscles. Will need to do sensitivity analysis as a result."\
					], \
		'Shoulder MA' : {	"Primary Source" : "Kuechle; 1997",\
		 					"Sources" : \
								[\
									MA_Settings(Pigeon_coeff_conversion([12.7928795,  2.0480346,  0.8917734,  3.2207214, -2.3928223,  0.]), 'Kuechle; 1997', "mm", None, None, 'Shoulder'),\
									MA_Settings([33.02,0,0,0,0,0], 'Pigeon; 1996', "mm", None, None, 'Shoulder'),\
									MA_Settings([1.9,0,0,0,0,0], 'Holzbaur; 2005', "cm", None, None, 'Shoulder'),\
									MA_Settings(19, "Est", "mm", None, None, 'Shoulder')\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Est", \
		 				"Sources" : \
							[\
								MA_Settings(0, "Est", "m", None, None, 'Elbow')\
							]}, \
		'Spindle' : Spindle_Settings(182/3,426.3/3,0.43,"Banks; 2006"),\
		'Mass' : {	"Primary Source" : "Banks; 2006", \
					"Sources" : \
						[\
							Input_Source(355.7/3, "Banks; 2006", 'g')\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(9.8, 'Holzbaur; 2005', 'cm')\
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(9.3, 'Holzbaur; 2005', 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(22, 'Holzbaur; 2005', 'degrees')\
									]}, \
		'PCSA' : {	"Primary Source" : "Holzbaur; 2005", \
					'Sources' : \
						[\
							Input_Source(82.98,'Garner & Pandy; 2003','sq cm'),\
							Input_Source(8.2,'Holzbaur; 2005','sq cm')\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Holzbaur; 2005", \
										'Sources': \
											[\
												Input_Source(2044.65,'Garner & Pandy; 2003','N'),\
												Input_Source(1142.6,'Holzbaur; 2005','N')\
											]}\
		}

	CB_Settings = {\
		'Notes' : [\
						"Holzbaur parameters for shoulder are for frontal plane (ABD/ADD) only! This explains the relatively small values for some measures. MA is negative as a result.",\
						"Garner & Pandy values for muscle length, PCSA, and peak force are very different from those reported in Wood (1989), Veeger (1991), Bassett (1990), Chen (1988), Keating (1993), Veeger (1997), An (1981), and Cutts (1991)." \
					],\
		'Shoulder MA' : {	"Primary Source" : "Est",\
		 					"Sources" : \
								[\
									MA_Settings(20, "Est", "mm", None, None, "Shoulder"),\
									MA_Settings([-20,0,0,0,0,0], "Holzbaur; 2005", "mm", None, None, "Shoulder")\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Est", \
		 				"Sources" : \
							[\
								MA_Settings(0, "Est", "m", None, None, 'Elbow')\
							]}, \
		'Spindle' : Spindle_Settings(123,147.3,0.83,"Banks; 2006"),\
		'Mass' : {	"Primary Source" : 'Banks; 2006',\
					'Sources' : \
						[\
							Input_Source(39.8, "Banks; 2006", 'g')\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(9.3, 'Holzbaur; 2005', 'cm'),\
											Input_Source(17.60, 'Garner & Pandy; 2003', 'cm')\
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : 'Holzbaur; 2005', \
									"Sources" : \
										[\
											Input_Source(9.7, 'Holzbaur; 2005', 'cm'),\
											Input_Source(4.23, 'Garner & Pandy; 2003', 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(27, "Holzbaur; 2005", 'deg')\
									]}, \
		'PCSA' : {	"Primary Source" : "Holzbaur; 2005", \
					'Sources' : \
						[\
							Input_Source(1.7,"Holzbaur; 2005","sq cm"),\
							Input_Source(4.55,"Garner & Pandy; 2003","sq cm")\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Holzbaur; 2005", \
										'Sources': \
											[\
												Input_Source(242.5, "Holzbaur; 2005", "N"),\
												Input_Source(150.02, "Garner & Pandy; 2003", "N")\
											]}\
		}

	DELTp_Settings = {\
		'Notes' : [\
						"DELTp SFE MA is listed as -78.74 mm in Pigeon. Using Pigeon Coefficients convention, Kuechle (1997) has the DELTp MA for [-140,90] as Pigeon_coeff_conversion([ 22.8547177,  3.9721238, -3.3900829, -3.6146546,  0.,  0.]). This will yield a piecewise function that creates jumps with the new velocity formulation. Instead, we are going to try Pigeon_coeff_conversion([-23.8165173, -4.486164 ,  5.8655808,  6.5003255, -8.2736695,2.0812998]) so that the function is within range and continuous during the ROM. Threshold (pi/2) has been removed for this new MA function.",\
						"Holzbaur parameters for shoulder are for frontal plane (ABD/ADD) only! This explains the relatively small values for some measures. MA is negative as a result.",\
						"Garner & Pandy values for muscle length, PCSA, and peak force are very different from those reported in Wood (1989), Veeger (1991), Bassett (1990), Chen (1988), Keating (1993), Veeger (1997), An (1981), and Cutts (1991). Also, they do no distinguish between ant, mid, post.",\
						"Holzbaur parameters for shoulder are for frontal plane (ABD/ADD) only! This explains the relatively small values for some measures.",\
						"Banks only had mass and spindle settings for the entire deltoid. As a result the parameters are divided by 3 as an estimate of the individual muscles. Will need to do sensitivity analysis as a result."\
					],\
		'Shoulder MA' : {	"Primary Source" : "Kuechle; 1997",\
		 					"Sources" : \
								[\
									MA_Settings(Pigeon_coeff_conversion([-23.8165173, -4.486164 ,  5.8655808,  6.5003255, -8.2736695,2.0812998]), 'Kuechle; 1997', "mm", None, None, "Shoulder"),\
									MA_Settings([-78.74,0,0,0,0,0], 'Pigeon; 1996', "mm", None, None, "Shoulder"),\
									MA_Settings([-8,0,0,0,0,0], 'Holzbaur; 2005', "mm", None, None, "Shoulder")
								]}, \
		'Elbow MA' : {	"Primary Source" : "Est", \
		 				"Sources" : \
							[\
								MA_Settings(0, "Est", "m", None, None, 'Elbow')\
							]}, \
		'Spindle' : Spindle_Settings(182/3,426.3/3,0.43,"Banks; 2006"),\
		'Mass' : {	"Primary Source" : "Banks; 2006", \
					"Sources" : \
						[\
							Input_Source(355.7/3, "Banks; 2006", 'g')\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(13.7, 'Holzbaur; 2005', 'cm'),\
											Input_Source(12.8, 'Garner & Pandy; 2003', 'cm')\
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(3.8, "Holzbaur; 2005", "cm"),\
											Input_Source(5.38, 'Garner & Pandy; 2003', 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(18, "Holzbaur; 2005", 'deg')\
									]}, \
		'PCSA' : {	"Primary Source" : "Holzbaur; 2005", \
					'Sources' : \
						[\
							Input_Source(1.9, "Holzbaur; 2005", 'sq cm'),\
							Input_Source(81.98,"Garner & Pandy; 2003","sq cm")\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Holzbaur; 2005", \
										'Sources': \
											[\
												Input_Source(259.9,"Holzbaur; 2005","N"),\
												Input_Source(2044.65,"Garner & Pandy; 2003","N")\
											]}\
		}

	BIC_Settings = {\
		'Notes' : [\
					"BIC EFE MA for Ramsay; 2009 has R² = 0.985 whereas Pigeon has R² = 0.9918. Pigeon, however, only takes elbow angle into account, whereas Ramsay; 2009 takes in variable PS angles. It appears that because Pigeon uses an average of fully pronated and fully supinated MAs, the BIC moment arm is similar but subject to variation as the level of PS is changed. (NOTE: BIC becomes slightly negative when q2 > 3.021. If trajectory has elbow angles exceding this value, enter a threshold of 3.021 into the model.)",\
					"Note: Only using the long head for optimal length, see Holzbaur (2005) for additional head parameters. Adding when logical."
					],\
		'Shoulder MA' : {	"Primary Source" : "Pigeon; 1996",\
		 					"Sources" : \
								[\
									MA_Settings([29.21,0,0,0,0,0], 'Pigeon; 1996', "mm", None, None, "Shoulder"),\
									MA_Settings(15, 'Est', "mm", None, None, "Shoulder")\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Ramsay; 2009", \
	 					"Sources" : \
							[\
								MA_Settings([8.4533,36.6147,2.4777,-19.432,2.0571,0,13.6502,0,0,-5.6172,0,-2.0854,0,0,0,0], 'Ramsay; 2009', "mm", 2, 3.021, 'Elbow'),\
								MA_Settings(Pigeon_coeff_conversion([14.660,4.5322,1.8047,-2.9883,0,0]), 'Pigeon; 1996', "mm", None, 2.9326, 'Elbow'),\
								MA_Settings([36,0,0,0,0,0], 'Holzbaur; 2005', "mm", None, None, "Elbow")
							]}, \
		'Spindle' : Spindle_Settings(320,292.6,1.1,"Banks; 2006"),\
		'Mass' : {	"Primary Source" : "Banks; 2006", \
					"Sources" : \
						[\
							Input_Source(163.8,"Banks; 2006", 'g')\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(11.6, 'Holzbaur; 2005', 'cm'),\
											Input_Source(14.22, "Garner & Pandy; 2003", "cm")
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(27.2, "Holzbaur; 2005", "cm"),\
											Input_Source(22.98, "Garner & Pandy; 2003", "cm")
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(0, "Holzbaur; 2005", 'deg')\
									]}, \
		'PCSA' : {	"Primary Source" : "Holzbaur; 2005", \
					'Sources' : \
						[\
							Input_Source((4.5+3.1), "Holzbaur; 2005", "sq cm"),\
							Input_Source(25.90, "Garner & Pandy; 2003", "sq cm")
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Holzbaur; 2005", \
										'Sources': \
											[\
												Input_Source((624.3+435.6), "Holzbaur; 2005", "N"),\
												Input_Source(849.29, "Garner & Pandy; 2003", "N")
											]}\
		}

	TRI_Settings = {\
		'Notes' : [\
					"TRI EFE MA for Ramsay; 2009 has R² = 0.997 whereas Pigeon has R² = 0.9904. Pigeon appears to really fail when the elbow angle is greater than 140°. For this reason, Ramsay; 2009 should be used. However the approach of fixing the MA for values greater than 140° can be adopted for completeness. Coefficients and equation number/type are listed below to test either implementation.",\
					"Note: Only using the long head for optimal length, see Holzbaur (2005) for additional head parameters.",\
					"Banks had the parameters for each head of the triceps, values were added.",\
					"Holzbaur settings only utilizes the long head of the TRI."\
					],
		'Shoulder MA' : {	"Primary Source" : "Pigeon; 1996",\
		 					"Sources" : \
								[\
									MA_Settings([-25.40,0,0,0,0,0], 'Pigeon; 1996', "mm", None, None, "Shoulder"),\
									MA_Settings(-15, 'Est', "mm", None, None, "Shoulder")\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Ramsay; 2009", \
	 					"Sources" : \
							[\
								MA_Settings([-24.5454,-8.8691,9.3509,-1.7518,0], 'Ramsay; 2009', 'mm', 1, None, 'Elbow'),\
								MA_Settings(Pigeon_coeff_conversion([-23.287,-3.0284,12.886,-19.092,13.277,-3.5171]), 'Pigeon; 1996', 'mm', None, None, 'Elbow'),\
								MA_Settings([-21,0,0,0,0,0], 'Holzbaur; 2005', 'mm', None, None, 'Elbow')
							]}, \
		'Spindle' : Spindle_Settings((200+222+98),(223.7+269.6+221.8),(0.89+0.82+0.44)/3,"Banks; 2006"),\
		'Mass' : {	"Primary Source" : "Banks; 2006",\
					"Sources" : \
						[\
							Input_Source((94.2+138.4+92.5), "Banks; 2006", 'g')\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005",\
		 							"Sources" : \
										[\
											Input_Source(13.4, 'Holzbaur; 2005', 'cm'),\
											Input_Source(8.77, 'Garner & Pandy; 2003', 'cm')\
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(14.3, "Holzbaur; 2005", 'cm'),\
											Input_Source(19.05, 'Garner & Pandy; 2003', 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(12, "Holzbaur; 2005", 'deg')\
									]}, \
		'PCSA' : {	"Primary Source" : "Holzbaur; 2005", \
					'Sources' : \
						[\
							Input_Source((5.7+4.5+4.5),"Holzbaur; 2005",'sq cm'),\
							Input_Source(76.30, 'Garner & Pandy; 2003', 'sq cm')\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Holzbaur; 2005", \
										'Sources': \
											[\
												Input_Source((798.5+624.3+624.3),"Holzbaur; 2005",'N'),\
												Input_Source(2332.92, 'Garner & Pandy; 2003', 'N')\
											]}\
		}

	BRA_Settings = {\
		"Notes" : [\
					"BRA (Brachialis) EFE MA for Ramsay; 2009 has R² = 0.990 whereas Pigeon has R² = 0.9988. Curve appears to be a better fit, as it experiences its smallest MA when Elbow angle = 0. Coefficients and equation number/type are listed below to test either implementation."\
					],
		'Shoulder MA' : {	"Primary Source" : "Est",\
		 					"Sources" : \
								[\
									MA_Settings(0, 'Est', 'm', None, None, 'Shoulder')\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Ramsay; 2009", \
	 					"Sources" : \
							[\
								MA_Settings([16.1991,-16.1463,24.5512,-6.3335,0], 'Ramsay; 2009', 'mm', 1, None, 'Elbow'),\
								MA_Settings(Pigeon_coeff_conversion([5.5492,2.3080,2.3425,-2.0530,0,0]), 'Pigeon; 1996', 'mm', None, None, 'Elbow'),\
								MA_Settings([18,0,0,0,0,0], 'Holzbaur; 2005', 'mm', None, None, 'Elbow')\
							]}, \
		'Spindle' : Spindle_Settings(256,272.1,0.94,"Banks; 2006"),\
		'Mass' : {	"Primary Source" : "Banks; 2006", \
					"Sources" : \
						[\
							Input_Source(141, "Banks; 2006", 'g')\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(8.6, 'Holzbaur; 2005', 'cm'),\
											Input_Source(10.28, 'Holzbaur; 2005', 'cm')\
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(5.4, "Holzbaur; 2005", 'cm'),\
											Input_Source(1.75, "Holzbaur; 2005", 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(0, "Holzbaur; 2005", 'deg')\
									]}, \
		'PCSA' : {	"Primary Source" : "Holzbaur; 2005", \
					'Sources' : \
						[\
							Input_Source(7.1,"Holzbaur; 2005",'sq cm'),\
							Input_Source(25.88,"Garner & Pandy; 2003",'sq cm')\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Holzbaur; 2005", \
										'Sources': \
											[\
												Input_Source(987.3,"Holzbaur; 2005","N"),\
												Input_Source(583.76,"Garner & Pandy; 2003","N")\
											]}\
		}

	BRD_Settings = {\
		"Notes" : [\
					"BRD (Brachioradialis) for Ramsay; 2009 has R² = 0.988 whereas Pigeon has R² = 0.9989. Pigeon, however, only takes elbow angle into account, whereas Ramsay; 2009 takes in variable PS angles. Coefficients and equation number/type are listed below to test either implementation."\
					],
		'Shoulder MA' : {	"Primary Source" : "Est",\
		 					"Sources" : \
								[\
									MA_Settings(0, 'Est', 'm', None, None, 'Shoulder')\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Ramsay; 2009", \
	 					"Sources" : \
							[\
								MA_Settings(	[15.2564,-11.8355,2.8129,-5.7781,44.8143,0,2.9032,0,0,-13.4956,0,-0.3940,0,0,0,0], 'Ramsay; 2009', 'mm', 2, None, 'Elbow'), \
								MA_Settings(	Pigeon_coeff_conversion([19.490,1.6681,10.084,-6.5171,0,0]), 'Pigeon; 1996', 'mm', None, None, 'Elbow')\
							]}, \
		'Spindle' : Spindle_Settings(70,190.2,0.37,"Banks; 2006"),\
		'Mass' : {	"Primary Source" : "Banks; 2006",\
					"Sources" : \
						[\
							Input_Source(64.7,"Banks; 2006", 'g')\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(17.3, 'Holzbaur; 2005', 'cm'),\
											Input_Source(27.03, 'Garner & Pandy; 2003', 'cm')\
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(13.3, "Holzbaur; 2005", 'cm'),\
											Input_Source(6.04, "Garner & Pandy; 2003", 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(0, "Holzbaur; 2005", 'deg')\
									]}, \
		'PCSA' : {	"Primary Source" : "Holzbaur; 2005", \
					'Sources' : \
						[\
							Input_Source(1.9,"Holzbaur; 2005",'sq cm'),\
							Input_Source(3.08,"Garner & Pandy; 2003",'sq cm')\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Holzbaur; 2005", \
										'Sources': \
											[\
												Input_Source(261.3,"Holzbaur; 2005",'N'),\
												Input_Source(101.56,"Garner & Pandy; 2003",'N')\
											]}\
		}

	PRO_Settings = {\
		'Notes' : [\
					""\
					],
		'Shoulder MA' : {	"Primary Source" : "Est",\
		 					"Sources" : \
								[\
									MA_Settings(0, 'Est', 'm', None, None, 'Shoulder')\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Ramsay; 2009", \
	 					"Sources" : \
							[\
								MA_Settings(	[11.0405,-1.0079,0.3933,-10.4824,-12.1639,-0.4369,36.9174,3.5232,-10.4223,21.2604,-37.2444,10.2666,-11.0060,14.5974,-3.9919,1.7526,-2.0089,0.5460], 'Ramsay; 2009', 'mm', 3, None, 'Elbow')\
							]}, \
		'Spindle' : Spindle_Settings(187.6,185.5,1.3,"Banks; 2006"),\
		'Mass' : {	"Primary Source" : "Banks; 2006",\
					"Sources" : \
						[\
							Input_Source(38.8, "Banks; 2006", 'g')\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(4.9, 'Holzbaur; 2005', 'cm')\
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(9.8, "Holzbaur; 2005", 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(10, "Holzbaur; 2005", 'deg')\
									]}, \
		'PCSA' : {	"Primary Source" : "Holzbaur; 2005", \
					'Sources' : \
						[\
							Input_Source(4.0,"Holzbaur; 2005",'sq cm')\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Holzbaur; 2005", \
										'Sources': \
											[\
												Input_Source(566.2,"Holzbaur; 2005",'N')\
											]}\
		}

	FCR_Settings = {\
		'Notes' : [\
					"FCR EFE MA is not listed in Ramsay; 2009 but Pigeon has a quadratic function with R² = 0.9975. Pigeon only takes elbow angle into account. Coefficients and equation number/type are listed below to test either implementation. EFE MA was estimated to be constant and 10 mm for this muscle. If you use Pigeon, make sure to only accept positive moment arm values, as this model fails outside the ROM. One option is to set the MA to the constant value (i.e., MA[theta>140°] = MA[140°])"\
					],\
		'Shoulder MA' : {	"Primary Source" : "Est",\
		 					"Sources" : \
								[\
									MA_Settings(0, 'Est', 'm', None, None, 'Shoulder')\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Pigeon; 1996", \
	 					"Sources" : \
							[\
								MA_Settings(Pigeon_coeff_conversion([0.9351,0.5375,-0.3627,0,0,0]), 'Pigeon; 1996', 'mm', None, 2.86, 'Elbow'),\
								MA_Settings(1, 'Est', 'cm', None, None, 'Elbow')\
							]}, \
		'Spindle' : Spindle_Settings(129,125.7,1.0,"Banks; 2006"),\
		'Mass' : (28.7, "Banks; 2006", 'g'),\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(6.3, 'Holzbaur; 2005', 'cm'),\
											Input_Source(5.10, 'Garner & Pandy; 2003', 'cm')\
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(24.4, "Holzbaur; 2005", 'cm'),\
											Input_Source(27.08, "Garner & Pandy; 2003", 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(3, "Holzbaur; 2005", 'deg')
									]}, \
		'PCSA' : {	"Primary Source" : "Holzbaur; 2005", \
					'Sources' : \
						[\
							Input_Source(1.6,"Holzbaur; 2005",'sq cm'),\
							Input_Source(11.16,"Garner & Pandy; 2003",'sq cm')\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Holzbaur; 2005", \
										'Sources': \
											[\
												Input_Source(74.0,"Holzbaur; 2005","N"),\
												Input_Source(368.63,"Garner & Pandy; 2003","N")\
											]}\
		}

	ECRB_Settings = {\
		'Notes' : [\
					"Garner and Pandy do not distinguish between ERCB and ECRL. Parameter values were included for completeness."\
					],\
		'Shoulder MA' : {	"Primary Source" : "Est",\
		 					"Sources" : \
								[\
									MA_Settings(0, 'Est', 'm', None, None, 'Shoulder')\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Ramsay; 2009", \
	 					"Sources" : \
							[\
								MA_Settings([-11.256,17.8548,1.6398,-0.5073,-2.8827,0,-0.0942,0,0,0,0,0,0,0,0,0], 'Ramsay; 2009', 'mm', 2, None, 'Elbow')\
							]}, \
		'Spindle' : Spindle_Settings(102,132.7,0.77,"Banks; 2006"),\
		'Mass' : {	"Primary Source" : "Banks; 2006",\
					"Sources" : \
						[\
							Input_Source(32.1, "Banks; 2006", 'g')\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(5.9, 'Holzbaur; 2005', 'cm'),\
											Input_Source(7.28, 'Garner & Pandy; 2003', 'cm')\
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(22.2, "Holzbaur; 2005", 'cm'),\
											Input_Source(26.80, "Garner & Pandy; 2003", 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(9, "Holzbaur; 2005", 'deg')\
									]}, \
		'PCSA' : {	"Primary Source" : "Holzbaur; 2005", \
					'Sources' : \
						[\
							Input_Source(2.2,"Holzbaur; 2005",'sq cm'),\
							Input_Source(24.89,"Garner & Pandy; 2003",'sq cm')\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Holzbaur; 2005", \
										'Sources': \
											[\
												Input_Source(100.5,"Holzbaur; 2005","N"),\
												Input_Source(755.76,"Garner & Pandy; 2003","N")\
											]}\
		}

	ECRL_Settings = {\
		'Notes' : [\
					"ECRL EFE MA for Ramsay; 2009 has R² = 0.978 whereas Pigeon has R² = 0.9986. Pigeon, however, only takes elbow angle into account, whereas Ramsay; 2009 takes in variable PS angles. Additionally, Pigeon only considers elbow angles between 0 and 140 degrees and exhibits a decrease in MA as elbow angle approaches the upper bound of the ROM. This should (intiutively speaking) make the extensors MA largest, but Pigeon exhibits a drop off that may make it less favorable for movements at the boundary of the ROM. Coefficients and equation number/type are listed below to test either implementation.",\
					"Garner and Pandy do not distinguish between ERCB and ECRL. Parameter values were included for completeness. Might need to divide by 2 to 'split' muscle."\
					],
		'Shoulder MA' : {	"Primary Source" : "Est",\
		 					"Sources" : \
								[\
									MA_Settings(0, 'Est', 'm', None, None, 'Shoulder')\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Pigeon; 1996", \
	 					"Sources" : \
							[\
								MA_Settings(Pigeon_coeff_conversion([4.7304,1.2590,4.4347,-3.0229,0,0]), 'Pigeon; 1996', 'mm', None, None, 'Elbow'),\
								MA_Settings([-7.7034,16.3913,7.4361,-1.7566,0,-1.3336,0,0.0742,0,0,0,0,0,0,0,0], 'Ramsay; 2009', 'mm', 2, None, 'Elbow')\
							]}, \
		'Spindle' : Spindle_Settings(74,155.2,0.48,"Banks; 2006"),\
		'Mass' : { "Primary Source" : "Banks; 2006", \
					"Sources" : \
						[\
							Input_Source(44.3, "Banks; 2006", 'g')\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(8.1, 'Holzbaur; 2005', 'cm'),\
											Input_Source(7.28, 'Garner & Pandy; 2003', 'cm')\
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(22.4, "Holzbaur; 2005", 'cm'),\
											Input_Source(26.80, "Garner & Pandy; 2003", 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(0, "Holzbaur; 2005", 'deg')\
									]}, \
		'PCSA' : {	"Primary Source" : "Holzbaur; 2005", \
					'Sources' : \
						[\
							Input_Source(2.2,"Holzbaur; 2005",'sq cm'),\
							Input_Source(24.89,"Garner & Pandy; 2003",'sq cm')\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Holzbaur; 2005", \
										'Sources': \
											[\
												Input_Source(304.9,"Holzbaur; 2005",'N'),\
												Input_Source(755.76,"Garner & Pandy; 2003",'N')\
											]}\
		}

	FCU_Settings = {\
		'Notes' : [\
					""\
					],
		'Shoulder MA' : {	"Primary Source" : "Est",\
		 					"Sources" : \
								[\
									MA_Settings(0, 'Est', 'm', None, None, 'Shoulder')\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Est", \
	 					"Sources" : \
							[\
								MA_Settings(1, 'Est', 'cm', None, None, 'Elbow')\
							]}, \
		'Spindle' : Spindle_Settings(175,141.2,1.2,"Banks; 2006"),\
		'Mass' : {	"Primary Source" : "Banks; 2006", \
					"Sources" : \
						[\
							Input_Source(36.5,"Banks; 2006", 'g')\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(5.1, 'Holzbaur; 2005', 'cm'),\
											Input_Source(3.98, 'Garner & Pandy; 2003', 'cm')\
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(26.5, "Holzbaur; 2005", 'cm'),\
											Input_Source(27.14, "Garner & Pandy; 2003", 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(12, "Holzbaur; 2005", 'deg')\
									]}, \
		'PCSA' : {	"Primary Source" : "Holzbaur; 2005", \
					'Sources' : \
						[\
							Input_Source(2.9,"Holzbaur; 2005",'sq cm'),\
							Input_Source(16.99,"Garner & Pandy; 2003",'sq cm')\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Holzbaur; 2005", \
										'Sources': \
											[\
												Input_Source(128.9,"Holzbaur; 2005",'N'),\
												Input_Source(561.00,"Garner & Pandy; 2003",'N')\
											]}\
		}

	FDS_Settings = {\
		'Notes' : [\
					"Note: only they muscle for the second digit was used for the FDS muscle. NEED TO DETERMINE IF THIS SHOULD BE A SUM OR AN AVERAGE FOR MEASURES LIKE PCSA, F_MAX, ETC.",\
					"As we are only considering one digit, do we need to sum all of the peak forces in order to get a better representation of its force producing capabilities?"\
					],
		'Shoulder MA' : {	"Primary Source" : "Est",\
		 					"Sources" : \
								[\
									MA_Settings(0, 'Est', 'm', None, None, 'Shoulder')\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Est", \
	 					"Sources" : \
							[\
								MA_Settings(1, 'Est', 'cm', None, None, 'Elbow')\
							]}, \
		'Spindle' : Spindle_Settings(356,224.9,1.6,"Banks; 2006"),\
		'Mass' : {	"Primary Source" : "Banks; 2006",\
					"Sources" : \
						[\
							Input_Source(95.2,"Banks; 2006", 'g')\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(8.4, 'Holzbaur; 2005', 'cm')\
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(27.5, "Holzbaur; 2005", 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(6, "Holzbaur; 2005", 'deg')\
									]}, \
		'PCSA' : {	"Primary Source" : "Holzbaur; 2005", \
					'Sources' : \
						[\
							Input_Source(1.4,"Holzbaur; 2005",'sq cm')\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Holzbaur; 2005", \
										'Sources': \
											[\
												Input_Source(61.2,"Holzbaur; 2005",'N')\
											]}\
		}

	PL_Settings = {\
		'Notes' : [\
					""\
					],
		'Shoulder MA' : {	"Primary Source" : "Est",\
		 					"Sources" : \
								[\
									MA_Settings(0, 'Est', 'm', None, None, 'Shoulder')\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Est", \
	 					"Sources" : \
							[\
								MA_Settings(1, 'Est', 'cm', None, None, 'Elbow')\
							]}, \
		'Spindle' : Spindle_Settings(None,None,None,None),\
		'Mass' : {	"Primary Source" : None,\
					"Sources" : \
						[\
							Input_Source(None, "N/A","N/A")\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(6.4, 'Holzbaur; 2005', 'cm')\
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(26.9, "Holzbaur; 2005", 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(4, "Holzbaur; 2005", 'deg')\
									]}, \
		'PCSA' : {	"Primary Source" : "Holzbaur; 2005", \
					'Sources' : \
						[\
							Input_Source(0.6,"Holzbaur; 2005",'sq cm')\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Holzbaur; 2005", \
										'Sources': \
											[\
												Input_Source(26.7,"Holzbaur; 2005",'N')\
											]}\
		}

	ECU_Settings = {\
		'Notes' : [\
					"ECU EFE MA is not listed in Ramsay; 2009 but Pigeon has a quadratic function with R² = 0.9966. Pigeon only takes elbow angle into account. Coefficients and equation number/type are listed below to test either implementation. EFE MA was estimated to be constant and -10 mm for this muscle. If you use Pigeon, make sure to only accept negative moment arm values, as this model fails outside the ROM. One option is to set the MA to the constant value (i.e., MA[theta>140°] = MA[140°])"\
					],
		'Shoulder MA' : {	"Primary Source" : "Est",\
		 					"Sources" : \
								[\
									MA_Settings(0, 'Est', 'm', None, None, 'Shoulder')\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Pigeon; 1996", \
	 					"Sources" : \
							[\
								MA_Settings(Pigeon_coeff_conversion([-2.1826,-1.7386,1.1491,0,0,0]), 'Pigeon; 1996', 'mm', None, None, 'Elbow'),\
								MA_Settings(-1, 'Est', 'cm', None, None, 'Elbow')\
							]}, \
		'Spindle' : Spindle_Settings(157,118,1.3,"Banks; 2006"),\
		'Mass' : {	"Primary Source" : "Banks; 2006",\
					"Sources" : \
						[\
							Input_Source(25.2,"Banks; 2006", 'g')\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(6.2, 'Holzbaur; 2005', 'cm'),\
											Input_Source(3.56, 'Garner & Pandy; 2003', 'cm')\
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(26.5, "Holzbaur; 2005", 'cm'),\
											Input_Source(28.18, "Garner & Pandy; 2003", 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(12, "Holzbaur; 2005", 'deg')\
									]}, \
		'PCSA' : {	"Primary Source" : "Holzbaur; 2005", \
					'Sources' : \
						[\
							Input_Source(2.9,"Holzbaur; 2005",'sq cm'),\
							Input_Source(8.04,"Garner & Pandy; 2003",'sq cm')\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Holzbaur; 2005", \
										'Sources': \
											[\
												Input_Source(128.9,"Holzbaur; 2005",'N'),\
												Input_Source(265.58,"Garner & Pandy; 2003",'N')\
											]}\
		}

	EDM_Settings = {\
		'Notes' : [\
					""\
					],
		'Shoulder MA' : {	"Primary Source" : "Est",\
		 					"Sources" : \
								[\
									MA_Settings(0, 'Est', 'm', None, None, 'Shoulder')\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Est", \
	 					"Sources" : \
							[\
								MA_Settings(-1, 'Est', 'cm', None, None, 'Elbow')\
							]}, \
		'Spindle' : Spindle_Settings(53,59.8,0.89,"Banks; 2006"),\
		'Mass' : {	"Primary Source" : "Banks; 2006", \
					"Sources" : \
						[\
							Input_Source(6.2, "Banks; 2006", 'g')\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(6.8, 'Holzbaur; 2005', 'cm')\
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : 'Holzbaur; 2005', \
									"Sources" : \
										[\
											Input_Source(32.2, 'Holzbaur; 2005', 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : 'Holzbaur; 2005', \
								"Sources" : \
									[\
										Input_Source(3, 'Holzbaur; 2005', 'deg')\
									]}, \
		'PCSA' : {	"Primary Source" : 'Holzbaur; 2005', \
					'Sources' : \
						[\
							Input_Source(0.6,'Holzbaur; 2005','sq cm')\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : 'Holzbaur; 2005', \
										'Sources': \
											[\
												Input_Source(25.3,'Holzbaur; 2005','N')\
											]}\
		}

	EDC_Settings = {\
		'Notes' : [\
					"Note: only they muscle for the second digit was used for the EDC muscle. NEED TO DETERMINE IF THIS SHOULD BE A SUM OR AN AVERAGE FOR MEASURES LIKE PCSA, F_MAX, ETC."\
					],
		'Shoulder MA' : {	"Primary Source" : "Est",\
		 					"Sources" : \
								[\
									MA_Settings(0, 'Est', 'm', None, None, 'Shoulder')\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Est", \
	 					"Sources" : \
							[\
								MA_Settings(-1, 'Est', 'cm', None, None, 'Elbow')\
							]}, \
		'Spindle' : Spindle_Settings(219,152.6,1.4,"Banks; 2006"),\
		'Mass' : {	"Primary Source" : "Banks; 2006", \
					"Sources" : \
						[\
							Input_Source(42.8, "Banks; 2006", 'g')\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(7.0, 'Holzbaur; 2005', 'cm')\
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(32.2, "Holzbaur; 2005", 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(3, "Holzbaur; 2005", 'deg')\
									]}, \
		'PCSA' : {	"Primary Source" : "Holzbaur; 2005", \
					'Sources' : \
						[\
							Input_Source(18.3,"Holzbaur; 2005",'sq cm')\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Holzbaur; 2005", \
										'Sources': \
											[\
												Input_Source(0.4,"Holzbaur; 2005",'N')\
											]}\
		}

	AN_Settings = {\
		'Notes' : [\
					""\
					],
		'Shoulder MA' : {	"Primary Source" : "Est",\
		 					"Sources" : \
								[\
									MA_Settings(0, 'Est', 'm', None, None, 'Shoulder')\
								]}, \
		'Elbow MA' : {	"Primary Source" : "Pigeon; 1996", \
	 					"Sources" : \
							[\
								MA_Settings(Pigeon_coeff_conversion([-5.3450,-2.2841,8.4297,-14.329,10.448,-2.736]), 'Pigeon; 1996', 'mm', None, None, 'Elbow')\
							]}, \
		'Spindle' : Spindle_Settings(None,None,None,None),\
		'Mass' : {	"Primary Source" : None, \
					"Sources" : \
						[\
							Input_Source(None, None, None)\
						]},\
		'Optimal Muscle Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(2.7,"Holzbaur; 2005",'cm')\
										]},\
		'Optimal Tendon Length' : {	"Primary Source" : "Holzbaur; 2005", \
									"Sources" : \
										[\
											Input_Source(1.8, "Holzbaur; 2005", 'cm')\
										]}, \
		'Pennation Angle' : {	"Primary Source" : "Holzbaur; 2005", \
								"Sources" : \
									[\
										Input_Source(0, "Holzbaur; 2005", 'deg')\
									]}, \
		'PCSA' : {	"Primary Source" : "Holzbaur; 2005", \
					'Sources' : \
						[\
							Input_Source(2.5,"Holzbaur; 2005",'sq cm')\
						]}, \
		'Maximum Isometric Force': {	"Primary Source" : "Holzbaur; 2005", \
										'Sources': \
											[\
												Input_Source(350.0,"Holzbaur; 2005",'N')\
											]}\
		}

	AllAvailableMuscles =[	"PC", "DELTa", "CB", "DELTp", "BIC", \
							"TRI", "BRA", "BRD", "PRO", "FCR",\
	 						"ECRB", "ECRL", "FCU", "FDS", "PL",\
	  						"ECU", "EDM", "EDC", "AN"]
	AllMuscleSettings = {	'PC': PC_Settings, 'DELTa' : DELTa_Settings, \
							'CB' : CB_Settings, 'DELTp' : DELTp_Settings,\
							'BIC' : BIC_Settings, 'TRI' : TRI_Settings, \
							'BRA' : BRA_Settings, 'BRD' : BRD_Settings, \
							'PRO' : PRO_Settings, 'FCR' : FCR_Settings, \
							'ECRB' : ECRB_Settings, 'ECRL' : ECRL_Settings, \
							'FCU' : FCU_Settings, 'FDS' : FDS_Settings, \
							'PL' : PL_Settings,'ECU' : ECU_Settings, \
							'EDM' : EDM_Settings, 'EDC' : EDC_Settings,\
							'AN' : AN_Settings}
	if PreselectedMuscles is None:
		ValidResponse_1 = False
		while ValidResponse_1 == False:
			MuscleSelectionType = input("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\nMuscle Selection:\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n (1) - Default\n (2) - Custom\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\nResponse: ")
			# import ipdb; ipdb.set_trace()
			print("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n")
			if MuscleSelectionType not in ['1','2','']:
				print('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\nInvalid Response! Please try again.')
				ValidResponse_1 = False
			elif MuscleSelectionType == '' or MuscleSelectionType == '1':
				for Muscle in ["PRO","AN"]:
					del(AllMuscleSettings[Muscle])
				ValidResponse_1 = True
			elif MuscleSelectionType == '2':
				ValidResponse_2 = False
				while ValidResponse_2 == False:
					MuscleListString = ""
					for i in range(len(AllAvailableMuscles)):
						MuscleListString += " (" + str(i+1) + ") - " + AllAvailableMuscles[i] + "\n"
					MuscleSelectionNumbers = input("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\nSelect Muscle Number(s)\n(separated by commas & groups with hyphens):\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n" + MuscleListString + "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\nMuscle Number(s): ")
					print("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n")
					MuscleSelectionNumbers = [el.strip() for el in MuscleSelectionNumbers.split(",")]
					for el in MuscleSelectionNumbers:
						if "-" in el:
							temp = el.split("-")
							MuscleSelectionNumbers.remove(el)
							[MuscleSelectionNumbers.append(str(i)) \
											for i in range(int(temp[0]),int(temp[1])+1)]
					if np.array([el in [str(i+1) for i in range(len(AllAvailableMuscles))] \
										for el in MuscleSelectionNumbers]).all() == False:
						print('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\nInvalid Response! Please check muscle numbers and try again.')
						ValidResponse_2 = False
					else:
						SelectedMuscles = [AllAvailableMuscles[int(el)-1] \
												for el in MuscleSelectionNumbers]
						MusclesToBeDeleted = [Muscle for Muscle in AllAvailableMuscles \
													if Muscle not in SelectedMuscles]
						for Muscle in MusclesToBeDeleted:
							del(AllMuscleSettings[Muscle])
						ValidResponse_2 = True
				ValidResponse_1 = True
	else:
		# assert type(PreselectedMuscles)==list and len(PreselectedMuscles)==8, "PreselectedMuscles, when used, must be a list of 8 numbers."
		assert np.array([type(MuscleNumber)==int for MuscleNumber in PreselectedMuscles]).all(),\
			"PreselectedMuscles must be a list of muscle numbers (ints)."
		assert np.array([MuscleNumber in range(1,len(AllAvailableMuscles)+1) \
			for MuscleNumber in PreselectedMuscles]).all(), \
				"PreselectedMuscles contains a muscle number outside the available muscles."
		SelectedMuscles = [AllAvailableMuscles[int(el)-1] \
								for el in PreselectedMuscles]
		MusclesToBeDeleted = [Muscle for Muscle in AllAvailableMuscles \
									if Muscle not in SelectedMuscles]
		for Muscle in MusclesToBeDeleted:
			del(AllMuscleSettings[Muscle])
	return(AllMuscleSettings)

def return_primary_source(Settings):
	assert Settings["Primary Source"]!=None, "No sources were found for this setting."
	TotalSources = Settings["Sources"]
	PrimarySource = Settings["Primary Source"]
	assert PrimarySource in [settings.Source for settings in TotalSources], "Error! Primary Source is not referenced."
	return(TotalSources[np.where([settings.Source == PrimarySource for settings in TotalSources])[0][0]])

def unit_conversion(Params):
    assert hasattr(Params,"Units"), "Params must have attr 'Units' in order for unit_conversion to work."
    assert hasattr(Params,"Values"), "Params must have attr 'Values' in order for unit_conversion to work."

    Units = Params.Units
    Value = Params.Values
    assert Units.capitalize() in \
            ["Degrees","Deg","Degree",
            "Radians","Radian","Rad",
            "Rads","In","Inches",
            "Cm","Centimeters","Centimeter",
            "Mm","Millimeters","Millimeter",
            "Meters","Meter","M",
            "Sq in","Squared inches","Inches squared",
            "In sq","Cm sq","Sq cm",
            "Centimeters squared","Squared centimeters","Mm sq",
            "Sq mm","Millimeters squared","Squared millimeters",
            "Meters squared","Squared meter","M sq",
            "Sq m","Lbs","Lb",
            "Pounds","G","Grams",
            "Gram","Kg","Kilograms",
            "Kilogram","N","Newton",
            "Newtons"], \
        "Improper Units Value. Please use appropriate Units."

    if Units.capitalize() in \
            ["Degrees","Deg","Degree",
            "Radians","Radian","Rad",
            "Rads"]:
    	if Units.capitalize() in ["Radians","Radian","Rad","Rads"]:
    		return(Value)
    	elif Units.capitalize() in ["Degrees","Deg","Degree"]:
    		if type(Value)==list:
    			return(list(np.array(Value)*np.pi/180))
    		else:
    			return(Value*np.pi/180)

    elif Units.capitalize() in \
            ["In","Inches","Cm",
            "Centimeters","Centimeter","Mm",
    		"Millimeters","Millimeter","Meters",
            "Meter","M"]:
    	if Units.capitalize() in ["Meter","Meters","M"]:
    		return(Value)
    	elif Units.capitalize() in ["In","Inches"]:
    		if type(Value)==list:
    			return(list(np.array(Value)*2.54/100))
    		else:
    			return(Value*2.54/100)
    	elif Units.capitalize() in  ["Cm","Centimeters","Centimeter"]:
    		if type(Value)==list:
    			return(list(np.array(Value)/100))
    		else:
    			return(Value/100)
    	elif Units.capitalize() in  ["Mm","Millimeters","Millimeter"]:
    		if type(Value)==list:
    			return(list(np.array(Value)/1000))
    		else:
    			return(Value/1000)

    elif Units.capitalize() in \
            ["Sq in","Squared inches","Inches squared",
            "In sq","Cm sq","Sq cm",
    		"Centimeters squared","Squared centimeters","Mm sq",
            "Sq mm","Millimeters squared","Squared millimeters",
            "Meters squared","Squared meter","M sq",
            "Sq m"]:
    	if Units.capitalize() in ["Meters squared","Squared meter","M sq","Sq m"]:
    		return(Value)
    	elif Units.capitalize() in ["Sq in","Squared inches","Inches squared","In sq"]:
    		if type(Value)==list:
    			return(list(np.array(Value)*((2.54/100)**2)))
    		else:
    			return(Value*((2.54/100)**2))
    	elif Units.capitalize() in ["Cm sq","Sq cm","Centimeters squared","Squared centimeters"]:
    		if type(Value)==list:
    			return(list(np.array(Value)/(100**2)))
    		else:
    			return(Value/(100**2))
    	elif Units.capitalize() in ["Mm sq","Sq mm","Millimeters squared","Squared millimeters"]:
    		if type(Value)==list:
    			return(list(np.array(Value)/(1000**2)))
    		else:
    			return(Value/(1000**2))

    elif Units.capitalize() in ["Lbs","Lb","Pounds","G","Grams","Gram","Kg","Kilograms","Kilogram"]:
    	if Units.capitalize() in ["Kg","Kilograms","Kilogram"]:
    		return(Value)
    	elif Units.capitalize() in ["G","Grams","Gram"]:
    		if type(Value)==list:
    			return(list(np.array(Value)/1000))
    		else:
    			return(Value/1000)
    	elif Units.capitalize() in  ["Lbs","Lb","Pounds"]:
    		if type(Value)==list:
    			return(list(np.array(Value)*0.45359237))
    		else:
    			return(Value*0.45359237)

    elif Units.capitalize() in ["N","Newton","Newtons"]:
    	if Units.capitalize() in ["N","Newton","Newtons"]:
    		return(Value)

max_stars_repo_name: danhagen/NonlinearControl

max_issues_repo_name: danhagen/NonlinearControl

max_forks_repo_name: danhagen/NonlinearControl

reason: import numpy,from numpy


text:
# Copyright 2019 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     https://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

from functools import partial

from absl.testing import absltest
import numpy as onp

from jax import lax
from jax import test_util as jtu
import jax.numpy as np

from examples import control

from jax.config import config
config.parse_flags_with_absl()
FLAGS = config.FLAGS

def one_step_lqr(dim, T):
  Q = np.stack(T * (np.eye(dim),))
  q = np.zeros((T, dim))
  R = np.zeros((T, dim, dim))
  r = np.zeros((T, dim))
  M = np.zeros((T, dim, dim))
  A = np.stack(T * (np.eye(dim),))
  B = np.stack(T * (np.eye(dim),))
  return control.LqrSpec(Q, q, R, r, M, A, B)

def control_from_lqr(lqr):
  T, dim, _ = lqr.Q.shape
  dot = np.dot

  def cost(t, x, u):
    return (
        dot(dot(lqr.Q[t], x), x) + dot(lqr.q[t], x) +
        dot(dot(lqr.R[t], u), u) + dot(lqr.r[t], u) +
        dot(dot(lqr.M[t], u), x))

  def dynamics(t, x, u):
    return dot(lqr.A[t], x) + dot(lqr.B[t], u)

  return control.ControlSpec(cost, dynamics, T, dim, dim)

def one_step_control(dim, T):

  def cost(t, x, u):
    return np.dot(x, x)

  def dynamics(t, x, u):
    return x + u

  return control.ControlSpec(cost, dynamics, T, dim, dim)

class ControlExampleTest(jtu.JaxTestCase):

  def testTrajectoryCyclicIntegerCounter(self):
    num_states = 3

    def dynamics(t, x, u):
      return (x + u) % num_states

    T = 10

    U = np.ones((T, 1))
    X = control.trajectory(dynamics, U, np.zeros(1))
    expected = np.arange(T + 1) % num_states
    expected = np.reshape(expected, (T + 1, 1))
    self.assertAllClose(X, expected, check_dtypes=False)

    U = 2 * np.ones((T, 1))
    X = control.trajectory(dynamics, U, np.zeros(1))
    expected = np.cumsum(2 * np.ones(T)) % num_states
    expected = np.concatenate((np.zeros(1), expected))
    expected = np.reshape(expected, (T + 1, 1))
    self.assertAllClose(X, expected, check_dtypes=False)

  def testTrajectoryTimeVarying(self):
    T = 6

    def clip(x, lo, hi):
      return np.minimum(hi, np.maximum(lo, x))

    def dynamics(t, x, u):
      # computes `(x + u) if t > T else 0`
      return (x + u) * clip(t - T, 0, 1)

    U = np.ones((2 * T, 1))
    X = control.trajectory(dynamics, U, np.zeros(1))
    expected = np.concatenate((np.zeros(T + 1), np.arange(T)))
    expected = np.reshape(expected, (2 * T + 1, 1))
    self.assertAllClose(X, expected, check_dtypes=True)

  def testTrajectoryCyclicIndicator(self):
    num_states = 3

    def position(x):
      '''finds the index of a standard basis vector, e.g. [0, 1, 0] -> 1'''
      x = np.cumsum(x)
      x = 1 - x
      return np.sum(x, dtype=np.int32)

    def dynamics(t, x, u):
      '''moves  the next standard basis vector'''
      idx = (position(x) + u[0]) % num_states
      return lax.dynamic_slice_in_dim(np.eye(num_states), idx, 1)[0]

    T = 8

    U = np.ones((T, 1), dtype=np.int32)
    X = control.trajectory(dynamics, U, np.eye(num_states, dtype=np.int32)[0])
    expected = np.vstack((np.eye(num_states),) * 3)
    self.assertAllClose(X, expected, check_dtypes=True)

  def testLqrSolve(self):
    dim, T = 2, 10
    p = one_step_lqr(dim, T)
    K, k = control.lqr_solve(p)
    K_ = -np.stack(T * (np.eye(dim),))
    self.assertAllClose(K, K_, check_dtypes=True, atol=1e-6, rtol=1e-6)
    self.assertAllClose(k, np.zeros((T, dim)), check_dtypes=True)

  def testLqrPredict(self):
    randn = onp.random.RandomState(0).randn
    dim, T = 2, 10
    p = one_step_lqr(dim, T)
    x0 = randn(dim)
    X, U = control.lqr_predict(p, x0)
    self.assertAllClose(X[0], x0, check_dtypes=True)
    self.assertAllClose(U[0], -x0, check_dtypes=True,
                        atol=1e-6, rtol=1e-6)
    self.assertAllClose(X[1:], np.zeros((T, 2)), check_dtypes=True,
                        atol=1e-6, rtol=1e-6)
    self.assertAllClose(U[1:], np.zeros((T - 1, 2)), check_dtypes=True,
                        atol=1e-6, rtol=1e-6)

  def testIlqrWithLqrProblem(self):
    randn = onp.random.RandomState(0).randn
    dim, T, num_iters = 2, 10, 3
    lqr = one_step_lqr(dim, T)
    p = control_from_lqr(lqr)
    x0 = randn(dim)
    X, U = control.ilqr(num_iters, p, x0, np.zeros((T, dim)))
    self.assertAllClose(X[0], x0, check_dtypes=True)
    self.assertAllClose(U[0], -x0, check_dtypes=True)
    self.assertAllClose(X[1:], np.zeros((T, 2)), check_dtypes=True)
    self.assertAllClose(U[1:], np.zeros((T - 1, 2)), check_dtypes=True)

  def testIlqrWithLqrProblemSpecifiedGenerally(self):
    randn = onp.random.RandomState(0).randn
    dim, T, num_iters = 2, 10, 3
    p = one_step_control(dim, T)
    x0 = randn(dim)
    X, U = control.ilqr(num_iters, p, x0, np.zeros((T, dim)))
    self.assertAllClose(X[0], x0, check_dtypes=True)
    self.assertAllClose(U[0], -x0, check_dtypes=True)
    self.assertAllClose(X[1:], np.zeros((T, 2)), check_dtypes=True)
    self.assertAllClose(U[1:], np.zeros((T - 1, 2)), check_dtypes=True)

  def testIlqrWithNonlinearProblem(self):
    def cost(t, x, u):
      return (x[0] ** 2. + 1e-3 * u[0] ** 2.) / (t + 1.)

    def dynamics(t, x, u):
      return (x ** 2. - u ** 2.) / (t + 1.)

    T, num_iters, d = 10, 7, 1
    p = control.ControlSpec(cost, dynamics, T, d, d)

    x0 = np.array([0.2])
    X, U = control.ilqr(num_iters, p, x0, 1e-5 * np.ones((T, d)))
    assert_close = partial(self.assertAllClose, atol=1e-2, check_dtypes=True)
    assert_close(X[0], x0)
    assert_close(U[0] ** 2., x0 ** 2.)
    assert_close(X[1:], np.zeros((T, d)))
    assert_close(U[1:], np.zeros((T - 1, d)))

  def testMpcWithLqrProblem(self):
    randn = onp.random.RandomState(0).randn
    dim, T, num_iters = 2, 10, 3
    lqr = one_step_lqr(dim, T)
    p = control_from_lqr(lqr)
    x0 = randn(dim)
    solver = partial(control.ilqr, num_iters)
    X, U = control.mpc_predict(solver, p, x0, np.zeros((T, dim)))
    self.assertAllClose(X[0], x0, check_dtypes=True)
    self.assertAllClose(U[0], -x0, check_dtypes=True)
    self.assertAllClose(X[1:], np.zeros((T, 2)), check_dtypes=True)
    self.assertAllClose(U[1:], np.zeros((T - 1, 2)), check_dtypes=True)

  def testMpcWithLqrProblemSpecifiedGenerally(self):
    randn = onp.random.RandomState(0).randn
    dim, T, num_iters = 2, 10, 3
    p = one_step_control(dim, T)
    x0 = randn(dim)
    solver = partial(control.ilqr, num_iters)
    X, U = control.mpc_predict(solver, p, x0, np.zeros((T, dim)))
    self.assertAllClose(X[0], x0, check_dtypes=True)
    self.assertAllClose(U[0], -x0, check_dtypes=True)
    self.assertAllClose(X[1:], np.zeros((T, 2)), check_dtypes=True)
    self.assertAllClose(U[1:], np.zeros((T - 1, 2)), check_dtypes=True)

  def testMpcWithNonlinearProblem(self):
    def cost(t, x, u):
      return (x[0] ** 2. + 1e-3 * u[0] ** 2.) / (t + 1.)

    def dynamics(t, x, u):
      return (x ** 2. - u ** 2.) / (t + 1.)

    T, num_iters, d = 10, 7, 1
    p = control.ControlSpec(cost, dynamics, T, d, d)

    x0 = np.array([0.2])
    solver = partial(control.ilqr, num_iters)
    X, U = control.mpc_predict(solver, p, x0, 1e-5 * np.ones((T, d)))
    assert_close = partial(self.assertAllClose, atol=1e-2, check_dtypes=True)
    assert_close(X[0], x0)
    assert_close(U[0] ** 2., x0 ** 2.)
    assert_close(X[1:], np.zeros((T, d)))
    assert_close(U[1:], np.zeros((T - 1, d)))

if __name__ == '__main__':
  absltest.main()

max_stars_repo_name: danielsuo/jax

max_issues_repo_name: danielsuo/jax

max_forks_repo_name: danielsuo/jax

reason: import numpy,import jax,from jax


text:
"""
This script 'biases' predictions based on the order of the examples. Naughty!

Usage: smoother.py <predictions_file> <smoothed_predictions_file>
"""

import numpy as np
import scipy.signal as sig
import sys

LABELS_PATH = "Y_train.npy"

MA_FACTOR = 110.0 #  100.0 # 200.0 # amount of moving average smoothing
SCALE = 0.3 # 1.0 # amount of biasing (additive)

if len(sys.argv) != 3:
    sys.exit("Usage: smoother.py <predictions_file> <smoothed_predictions_file>")

src_path = sys.argv[1]
tgt_path = sys.argv[2]

print "Loading predictions"
with open(src_path, 'r') as f:
    lines = f.readlines()

scores = np.array([float(l.strip()) for l in lines])

print "Loading training labels"
Y = np.load(LABELS_PATH)

print "Smoothing and rescaling training labels"
Ys = np.convolve(Y, np.ones(MA_FACTOR) / MA_FACTOR, mode='same')
Ysr = sig.resample(Ys, scores.shape[0])

print "Biasing predictions"
scores_biased = scores + SCALE * (Ysr - 0.5)

print "Store biased predictions"
with open(tgt_path, 'w') as f:
    for score in scores_biased:
        f.write(str(score) + '\n')

print "  stored in %s" % tgt_path

max_stars_repo_name: Keesiu/meta-kaggle

max_issues_repo_name: Keesiu/meta-kaggle

max_forks_repo_name: Keesiu/meta-kaggle

reason: import numpy,import scipy


text:
#Trabalhando com sorted no Pandas
import pandas as pd
import numpy as np
unsorted_df = pd.DataFrame(np.random.randn(10,2),index=[1,4,6,2,3,5,9,8,0,7],columns = ['col2','col1'])
sorted_df=unsorted_df.sort_index()
print (sorted_df)


text:
#!/usr/bin/env python
# coding: utf-8

# # DESIGN OF EXPERIMENTS

# This notebook describes the computations needed to do a 1-Way Design.

# # Overview

# This lecture illustrates a workflow for developing codes that run studies to do a one-factor-at-a-time (OFT) experimental design.
# The running example is the model of glycolytic oscillations developed in J. Wolf et al., Biochem., 2000.

# The function ``runStudy`` encapsulates the full workflow. This consists of:
# 1. running simulations for the factors and levels to be studied
# 1. calculating response values from the simulation results
# 1. constructing a table (dataframe) that reports the responses $y$ by factor level
# 1. calculating the $\mu$, $\alpha$ for design of experiments
#
# The notebook contains codes for (1)-(3) for the Wolf model.
# Responses are the oscillation frequency for glucose.

# # Preliminaries

# In[51]:

IS_COLAB = True

# In[52]:

# Code installs. Only needed once.
if IS_COLAB:
    pass

# In[53]:

import urllib.request
def getSharedCodes(moduleName):
  """
  Obtains common codes from the github repository.

  Parameters
  ----------
  moduleName: str
      name of the python module in the src directory
  """
  if IS_COLAB:
      url = "https://github.com/sys-bio/network-modeling-summer-school-2021/raw/main/src/%s.py" % moduleName
      local_python = "python.py"
      _, _ = urllib.request.urlretrieve(url=url, filename=local_python)
  else:
      local_python = "../../src/%s.py" % moduleName
  with open(local_python, "r") as fd:
    codeStr = "".join(fd.readlines())
  print(codeStr)
  exec(codeStr, globals())

# Acquire codes
getSharedCodes("util")

# TESTS
assert(isinstance(LINEAR_PATHWAY_DF, pd.DataFrame))

# In[54]:

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import numpy as np
import pandas as pd
from scipy import fftpack
import tellurium as te

# In[55]:

# Constants
AMPLIDX = 1  # Index of amplitude values
FONTSIZE = 16
FREQIDX = 0  # Index of frequency values
TIME = "time"
SMALLEST_PCT = -100  # Smallest percent change in a parameter value

# # Wolf Model

# We begin with an exploration of the Wolf model to understand its parameters and to determine what responses should be reported by an experiment.

# ## Constants for Model

# In[56]:

print(WOLF_MODEL)

# In[57]:

# Constants used in the analysis of the Wolf model
START = 0
END = 5
NUMPT = 300  # number of points to simulate
WOLF_FACTORS = [
  "J0_inputFlux", "J1_k1", "J1_Ki", "J1_n", "J2_k", "J3_k", "J4_kg", "J4_kp",
 "J4_ka", "J4_kk", "J5_k", "J6_k", "J7_k", "J8_k1", "J8_k2", "J9_k", "J10_k",
]

# ## Exploration of the Model

# In[58]:

def timePlots(data=WOLF_ARR, title=""):
    fig, ax = plt.subplots(1, figsize=(8, 4))
    colnames = WOLF_DF.columns
    p = ax.plot(data[:, 0], data[:,1:])
    _ = ax.legend(p, colnames[1:], bbox_to_anchor=(1.05, 1), loc='upper left')
    _ = ax.set_title(title)
timePlots()

# # Planning the Study

# 1. What are the research questions? In this study, we want to understand how the kinetic contants affect oscillations.
# 2. What are the factors and levels? The factors are the kinetic constants. The levels are percentage changes in the values of the contants compared with the baseline.
# 1. What are the responses? We calculate the Fast Fourier Transform (FFT) to obtain pairs of frequency and amplitude (the fourier coefficient).

# # Running Simulations With Changes in Parameter Values

# The first thing we need to do is to be able to run an experiment when a parameter is changed.
# We consider multiple levels for the parameter value.

# In[59]:

def runSimulations(parameter, percents, isPlot=True):
    """
    Runs experiments and produces plots for the fractional changes in parameter values.

    Parameters
    ----------
    parameter: str
    percents: list-float
        list of percent change in values
          (smallest value is -100)

    Returns
    -------
    dict: key=pct, value=data
    """
    roadrunner = te.loada(WOLF_MODEL)
    baseValue = roadrunner[parameter]
    if not parameter in roadrunner.keys():
        raise ValueError("Unknown parameter name: %s" % parameter)
    resultDct = {}
    for percent in percents:
        roadrunner.resetAll()
        percent = max(percent, SMALLEST_PCT)
        frac = 0.01*percent
        newValue = baseValue*(1 + frac)
        roadrunner[parameter] = newValue
        data = roadrunner.simulate(START, END, NUMPT)
        title = "%s: %f (%d%%)" % (parameter, newValue, percent)
        resultDct[percent] = data
        if isPlot:
            timePlots(data=data, title=title)
    return resultDct

# Tests
percents = [-10, 0, 10]
resultDct = runSimulations("J1_Ki", percents, isPlot=False)
assert(len(set(percents).symmetric_difference(resultDct.keys())) == 0)

# Here's an example of using these codes.

# In[60]:

# Experiments
changes = 2*np.array(range(5))
length = len(changes)
_ = runSimulations("J1_Ki", changes)

# **Questions**
# 1. What happens when J1_Ki is reduced?
# 1. Do other parameters behave in a similar way in that the onset of cyclic behavior is delayed with a larger change?
# 1. What metric(s) should we use to quantify the effect on cyclic behavior?

# # Calculating Responses

# One way to quantify the response to changes in parameters is to use the Fast Fourier Transform (FFT).
# This section develops codes that calculate responses.

# In[61]:

def calculateFft(molecule, data, offset=100):
    """
    Calculate the FFT for a molecule in the simulation output.
    The calculation does not include amplitudes at a frequency of 0.

    Parameters
    ----------
    molecule: str
    data: NamedArray/DataFrame
    offset: int
        Initial data that are not included in the FFT calculation

    Returns
    -------
    list-float, list-float
        freqs, fftValues
    """
    # Returns frequencies and abs(fft) for a chemical species (molecule)
    if isinstance(data, pd.DataFrame):
        data = data.copy()
        colnames = ["[" + c + "]" for c in data.columns]
        data.columns = colnames
        times = np.array(list(data.index))
    else:
        times = data["time"]
    col = "[%s]" % molecule
    values = np.array(list(data[col]))
    numPoint = len(data)
    count = numPoint - offset
    endTime = times[-1]
    freqs = fftpack.fftfreq(count, endTime/numPoint)
    fftValues = np.abs(fftpack.fft(values[offset:]))
    # Eliminate frequency of 0
    freqs = freqs[1:]
    fftValues = fftValues[1:]
    return freqs, fftValues

# Tests
freqs, fftValues = calculateFft("Glucose", WOLF_DF)
assert(max(fftValues) > 90)  # Top frequency should have a large magnitude

# In[62]:

def plotFft(molecule, data=WOLF_DF, ax=None, **kwargs):
    """
    Plot the FFT for a molecule.

    Parameters
    ----------
    data: NamedArray
    ax: Matplotlib.Axes
    kwargs: dict
        optional arguments passed to FFT calculation
    """
    if ax is None:
        _, ax = plt.subplots(1,1)
    freqs, fftValues = calculateFft(molecule, data=data, **kwargs)
    ax.plot(freqs, fftValues)
    ax.set_title(molecule)
    ax.set_ylim([0, 130])

# Tests
plotFft("Glucose")  # "Smoke" test

# In[63]:

def calculatePeakFrequency(molecule, data, numDelete=0, offset=100):
    """
    Calculate the FFT for a molecule in the simulation output.
    The calculation does not include amplitudes at a frequency of 0.

    Parameters
    ----------
    molecule: str
    data: NamedArray/DataFrame
    offset: int
        Initial data that are not included in the FFT calculation

    Returns
    -------
    float: peak frequency
    """
    def getIndexOfMax(arr):
        maxx = max(arr)
        lst = list(arr)
        return lst.index(maxx)
    def deleteArrayElement(arr, idx):
        lst = list(arr)
        lst1 = lst[:idx]
        lst2 = lst[idx+1:]
        newLst = lst1.extend(lst2)
        return np.array(lst1)
    #
    frequencies, amplitudes = calculateFft(molecule, data)
    # Delete the largest N values
    for _ in range(numDelete):
        maxIdx = getIndexOfMax(amplitudes)
        frequencies = deleteArrayElement(frequencies, maxIdx)
        amplitudes = deleteArrayElement(amplitudes, maxIdx)
    # Report the next highest value
    maxIdx = getIndexOfMax(amplitudes)
    return frequencies[maxIdx]

# TESTS
assert(calculatePeakFrequency("Glucose", WOLF_DF) > 0)

# # Calculate Responses

# In[64]:

def calculateFrequencyResponses(factors=WOLF_FACTORS, percents=[0], molecule="Glucose"):
    """
    Calculates the frequency responses for a 1WD for the factors and
    levels expressed as percents.

    Parameters
    ----------
    factors: list-str (kinetic constants)
    percents: float

    Returns
    -------
    pd.DataFrame
       rows: factor
       columns: percents
       values: frequency
    """
    percentDct = {p: [] for p in percents}
    for factor in factors:
        dataDct = runSimulations(factor, percents, isPlot=False)
        # Calculate the FFT for all percents for this factor
        for percent in dataDct.keys():
            frequency = calculatePeakFrequency(molecule, dataDct[percent])
            percentDct[percent].append(frequency)
    return pd.DataFrame(percentDct, index=factors)

df = calculateFrequencyResponses(percents=[-2, 0, 2])
assert(isinstance(df, pd.DataFrame))

# In[65]:

PERCENTS = [-20, -10, -5, -2, 0, 2, 5, 10, 20]
RESPONSE_DF = calculateFrequencyResponses(percents=PERCENTS)
RESPONSE_DF

# **Questions**
# 1. Which parameters show an effect on frequency for the levels considered?
# 1. Of these, which parameter values are inversely related to frequency?

# # Calculating $\mu$, $\alpha$

# $\mu$ is the baseline for the simulation.
# $\alpha_{i, k_i} = y_{i, k_i} - \mu$

# In[66]:

def calculateMuAlpha(responseDF):
    """
    Calculates the design of experiment parameters for a 1-Way Design.
    Assumes that there is a level of 0.

    Parameters
    ----------
    responseDF: pd.DataFrame
        column: levels
        index: factors
        values: y_i,k_i

    Returns
    -------
    float: Mu
    pd.DataFrame
        columns: levels (k_i)
        index: factors (i)
        values: \alpha_i,k_i
    """
    mu = responseDF.loc["J1_k1", 0]
    resultDF = responseDF - mu
    return mu, resultDF

# TESTS
mu, responseDF = calculateMuAlpha(RESPONSE_DF)
assert(np.abs(mu - 5.1) < 0.1)
assert(isinstance(responseDF, pd.DataFrame))

# In[67]:

mu, responseDF = calculateMuAlpha(RESPONSE_DF)
responseDF

# **Questions**
# 1. In many cases $\alpha_{i,k_i}$ is 0. What does this mean?

# # Running a Study

# A study is a collection of experiments in which one or more factors is set to level different from the baseline.
# For a 1WD, only one factor at a time is changed.

# In[68]:

def runStudy(factors=WOLF_FACTORS, percents=[0], molecule="Glucose"):
    """
    Runs a 1WD for the factors and levels expressed as percents.

    Parameters
    ----------
    factors: list-str (kinetic constants)
    percents: float

    Returns
    -------
    float: mu
    pd.DataFrame: alpha_i, k_i
       rows: factor
       columns: percents
       values: frequency
    """
    df = calculateFrequencyResponses(factors=factors, percents=percents,
                                     molecule=molecule)
    return calculateMuAlpha(df)

mu, df = runStudy(percents=[-2, 0, 2])
assert(isinstance(df, pd.DataFrame))
assert(isinstance(mu, float))

# Below we calculate the effect of parameter changes on NADH.

# In[69]:

mu, df = runStudy(percents=[-20, -10, -5, -2, 0, 2, 5, 10, 20], molecule="NADH")
df

# # Exercise

# Compare the effect of parameters on Glucose with their effects on NADH. Which parameters have more of an effect on Glucose than on NADH?

# ## Answer

# In[70]:

PERCENTS = [-20, -10, -5, -2, 0, 2, 5, 10, 20]
_, df1 = runStudy(percents=PERCENTS, molecule="Glucose")
_, df2 = runStudy(percents=PERCENTS, molecule="NADH")

# In[71]:

# Effects
df1 - df2

reason: import numpy,from scipy


text:
import os
from subprocess import getoutput

import numpy as np
import numbers
from matplotlib import pyplot as plt
from matplotlib.font_manager import FontProperties
from sklearn import metrics

DPI = 100

def mpl_fontprop():
    fc_list = [x.strip() for x in getoutput('fc-list :lang=zh').strip().split('\n')]
    fc = None
    for item in fc_list:
        if 'CN' in item:
            fc = item.strip().split()[0][:-1]
    if fc:
        fontp = FontProperties(fname=fc)
    else:
        fontp = FontProperties()
    return fontp

def smooth_line(x, y):
    newx = np.linspace(np.min(x), np.max(x), 100)
    newy = np.interp(newx, x, y)
    return newx, newy

def optimize_curve(a, b):
    optimal_idxs = np.where(np.r_[True,
                                  np.logical_or(np.diff(a, 2),
                                                np.diff(b, 2)),
                                  True])[0]
    a = a[optimal_idxs]
    b = b[optimal_idxs]
    return a, b, optimal_idxs

def roc_curve_ci(auc, pos_num, neg_num, ci=0.95):
    from scipy import stats as st
    q0 = auc * (1 - auc)
    q1 = auc / (2 - auc) - auc ** 2
    q2 = 2 * (auc ** 2) / (1 + auc) - auc ** 2
    se = np.sqrt((q0 + (pos_num - 1) * q1 + (neg_num - 1) * q2) / (pos_num * neg_num))
    alpha = (1 - ci) / 2
    z_crit = st.norm.ppf(1 - alpha)
    lower = auc - z_crit * se
    upper = auc + z_crit * se
    lower = max(lower, 0.)
    upper = min(upper, 1.)
    return lower, upper

def plot_curve(xs, ys, labels, xlabel='x', ylabel='y', title=None, save_path=None, diagonal=False, smooth=True, fscore=False,
               pos_nums=None, neg_nums=None):
    fig = plt.figure(dpi=DPI)
    fontprop = mpl_fontprop()
    ax = fig.add_subplot(111)
    ax.set_xlabel(xlabel, fontproperties=fontprop)
    ax.set_ylabel(ylabel, fontproperties=fontprop)
    if title is not None:
        ax.set_title(title, fontproperties=fontprop)
    auc_key = 'AUC' if diagonal else 'AP'
    for i, (x, y) in enumerate(zip(xs, ys)):
        y = np.asarray(y)
        x = np.asarray(x)
        # ROC
        auc_l = np.nan_to_num(metrics.auc(x, y))
        if diagonal == True:
            # x = np.concatenate([np.array([0]), x])
            # y = np.concatenate([np.array([0]), y])
            y[x == 0] = 0
        else:
            # sortx
            order = x.argsort()[::-1]
            x = x[order]
            y = y[order]
            y[x == 1] = 0
        x, y, thresh = optimize_curve(x, y)
        # sort end
        if smooth:
            x, y = smooth_line(x, y)
        if pos_nums is not None and neg_nums is not None:
            lower_l, upper_l = roc_curve_ci(auc_l, pos_nums[i], neg_nums[i])
            plt.plot(x, y, label='{}: {} {:0.3f} \n(95%CI:{:0.3f}-{:0.3f})'.format(labels[i], auc_key, auc_l, lower_l,
                                                                                upper_l),
                     linestyle='-',
                     # c='#03DAC5'
                     )
        else:
            plt.plot(x, y, label='{}={:.03f}'.format(labels[i], auc_l), linestyle='-')
    if diagonal:
        plt.plot([0, 1], [0, 1], color='gray', lw=2, linestyle='--', alpha=.8)
    if fscore:
        f_scores = np.linspace(0.2, 0.8, num=4)
        lines, labels = [], []
        for f_score in f_scores:
            x = np.linspace(0.01, 1)
            y = f_score * x / (2 * x - f_score)
            (l,) = plt.plot(x[y >= 0], y[y >= 0], color="gray", alpha=0.2)
            plt.annotate("f1={0:0.1f}".format(f_score), xy=(0.9, y[45] + 0.02))
            ax.set_xlim([0.0, 1.0])
            ax.set_ylim([0.0, 1.05])
    # plt.legend(prop=fontprop, loc=4)
    plt.legend(prop=fontprop, loc='best')
    if save_path is not None:
        plt.savefig(os.path.join(save_path, '{}.png'.format(title)))

def dca(preds, y_trues, labels, title=None):
    bottom = -0.05
    op_bottom = -0.1
    N = len(y_trues[0])
    pos_num = np.sum(y_trues[0])
    neg_num = N - pos_num
    for p, y, l in zip(preds, y_trues, labels):
        threshs = []
        nbs = []
        lower = 1.
        for th in np.unique(p):
            tp = np.sum((p >= th) * y)
            fp = np.sum((p >= th) * (y != 1))
            net_benefit = (tp - fp * th / (1 - th)) / N
            threshs.append(th)
            # nbs.append(net_benefit)
            lower = min(lower, net_benefit)
            if lower >= bottom:
                nbs.append(lower)
            else:
                break

        # remove null thresh curve
        if len(threshs[:len(nbs)]) < 1:
            print(f"{l} dca is invalid")
            continue

        threshs, nbs = smooth_line(threshs[:len(nbs)], nbs)
        plt.plot(threshs, nbs, label=l)

    # operation curve
    xs = np.linspace(0, 0.5, 20).tolist()
    operations = []
    for idx, NRI_th in enumerate(xs):
        NRI = pos_num / N - neg_num * (NRI_th / (1 - NRI_th)) / N
        if NRI >= op_bottom:
            operations.append(NRI)
        else:
            break
    plt.plot(xs[:len(operations)], operations, label='all with operation')
    plt.plot([0, 1], [0, 0], 'k:', label='none operation')
    if title is not None:
        plt.title(title)
    plt.legend()

def binary_tfpn(match, scores):
    desc_score_indices = np.argsort(scores, kind="mergsort")[::-1]
    match = match[desc_score_indices]
    scores = scores[desc_score_indices]
    # diff means gradient. remove duplicated score values as below
    distinct_value_indices = np.where(np.diff(scores))[0]
    threshold_idxs = np.r_[distinct_value_indices, match.size - 1]
    tps = np.cumsum(match)[threshold_idxs]
    fps = 1 + threshold_idxs - tps
    thresholds = scores[threshold_idxs]
    total_gt_pos = np.sum(match == 1)
    total_gt_neg = np.sum(match == 0)
    fns = total_gt_pos - tps
    tns = total_gt_neg - fps

    return tps, fps, tns, fns, thresholds

def simple_metric(gts, preds, cutoff=0.5):
    """

    :param gts: ground truth [0, 1] array
    :param preds: predictions 0-1 float array
    :param cutoff: could be a number as model threshold or statistics index such as f1-score/g-mean
    :return:
    """
    if not isinstance(preds, np.ndarray):
        preds = np.asarray(preds)
    if not isinstance(gts, np.ndarray):
        gts = np.asarray(gts)
    tps, fps, tns, fns, threshs = binary_tfpn(gts, preds)
    sim_metr = {"tps": tps,
                "fps": fps,
                "tns": tns,
                "fns": fns,
                "sensitivity": tps / (tps + fns),
                "precision": tps / (tps + fps),
                "specificity": tns / (tns + fps),
                "accuracy": (tps+tns)/(tps+tns+fps+fns)}
    sim_metr["f1-score"] = 2 * (sim_metr["precision"] * sim_metr["sensitivity"]) / (
            sim_metr["precision"] + sim_metr["sensitivity"])
    sim_metr["g-mean"] = np.sqrt(sim_metr["sensitivity"] * sim_metr["specificity"])
    sim_metr["youden-index"] = sim_metr["sensitivity"] + sim_metr["specificity"] - 1
    if isinstance(cutoff, numbers.Number):
        target_ind = np.max(np.where(threshs > cutoff)[0])
    elif isinstance(cutoff, str):
        if cutoff in sim_metr.keys():
            target_ind = np.argsort(sim_metr[cutoff], kind="mergsort")[::-1][0]
            # print(f"Set threshold to {threshs[target_ind]}")
            cutoff = threshs[target_ind]
        else:
            raise ValueError(f"{cutoff} is not available. keys are {sim_metr.keys()}")
    else:
        raise ValueError(f"{cutoff} is not available. number of str is required.")
    results = {}
    for k, v in sim_metr.items():
        results[k] = v[target_ind]
    results['youden-index'] = np.max(sim_metr['youden-index'])
    results['auroc'] = metrics.roc_auc_score(gts, preds)
    results['average_precision'] = metrics.average_precision_score(gts, preds)
    return results, cutoff