Simplifying loop manipulation

Initially, all the code was written with respect to system loops. Also, when manipulating loops either for removing stiffness or for any other reason, the loops themselves were manipulated. And this means the entire loop information, including the branch elements. All this only adds to overhead. Particularly since two other elements have been defined - branch_params which contains information of the branch - the elements, the resistance, inductance, voltages, and current. The other element is the system_loop_map which contains information of which branch exists in each loop.

So why not get rid of system_loops in further manipulations? All that is needed to be done is to extend the definition of system_loop_map. Before system_loop_map had only "stiff", "yes" and "no" as entries. Now it will have:

forward - branch is in the direction of loop
reverse - branch is against loop direction
stiff_forward - stiff branch but same direction
stiff_reverse - stiff branch but opposite direction

With this the entire information of the loop is contained in a far more compact manner. And all the functions were rewritten. For example, the remove_stiffness function is as follows (click on "view raw" below the code box to see the code in another window):

	def remove_stiffness(matrix_e, matrix_a, matrix_b, dt, sys_loops, branch_info, stiff_info, sys_loop_map, src_list):
	""" Reduces the stiffness of the ODE by limiting the
	appearance of stiff branches in the loops."""
	# First to find the stiff loops by locating the first branch
	# that is stiff. After that make all the loops below that loop
	# have no element for that stiff branch.
	# To first arrange the loops such that the first branch
	# that is stiff occurs in the first loop and so on.
	map_loop_count=0
	map_branch_count=0
	# Start with the first branch of the first loop.
	while map_loop_count<len(sys_loop_map):
	if sys_loop_map[map_loop_count][map_branch_count]=="stiff_forward" or sys_loop_map[map_loop_count][map_branch_count]=="stiff_reverse":
	loop_found="yes"
	else:
	loop_found="no"
	# If a loop is not found and the branches are not exhausted.
	while loop_found=="no" and map_branch_count<len(sys_loop_map[0]):
	# Check if the branch is stiff in the loop.
	if sys_loop_map[map_loop_count][map_branch_count]=="stiff_forward" or sys_loop_map[map_loop_count][map_branch_count]=="stiff_reverse":
	loop_found="yes"
	# If not, look at all the remaining loops corresponding to that branch.
	# The idea is that if a branch is stiff, it will appear in at least
	# one loop and so need to look through all loops to make sure it is
	# not missed.
	c1=map_loop_count+1
	while loop_found=="no" and c1<len(sys_loop_map):
	if sys_loop_map[c1][map_branch_count]=="stiff_forward" or sys_loop_map[c1][map_branch_count]=="stiff_reverse":
	for c2 in range(len(sys_loop_map[0])):
	# If a subsequenct loop is found to have the branch as stiff,
	# interchange the rows.
	sys_loop_map[c1][c2], sys_loop_map[map_loop_count][c2] = sys_loop_map[map_loop_count][c2], sys_loop_map[c1][c2]
	loop_found="yes"
	c1+=1
	# If all loops are exhausted and no element is found as stiff
	# for the branch, it means the branch is not stiff
	# or has been accounted in a previous loop.
	# So move on to the next branch.
	if loop_found=="no":
	map_branch_count+=1
	# If a loop has been found, the loops need to be made as upper triangular
	# So all stiff branches in the loops subsequent to the loop will have to eliminated.
	if loop_found=="yes":
	for c1 in range(map_loop_count+1, len(sys_loop_map)):
	if sys_loop_map[c1][map_branch_count]=="stiff_forward" or sys_loop_map[c1][map_branch_count]=="stiff_reverse":
	# This is simply row operation of loops.
	if sys_loop_map[map_loop_count][map_branch_count]==sys_loop_map[c1][map_branch_count]:
	for c2 in range(len(sys_loop_map[0])):
	if sys_loop_map[map_loop_count][c2]=="forward" and sys_loop_map[c1][c2]=="forward":
	sys_loop_map[c1][c2]="no"
	elif sys_loop_map[map_loop_count][c2]=="reverse" and sys_loop_map[c1][c2]=="reverse":
	sys_loop_map[c1][c2]="no"
	elif sys_loop_map[map_loop_count][c2]=="reverse" and sys_loop_map[c1][c2]=="no":
	sys_loop_map[c1][c2]="forward"
	elif sys_loop_map[map_loop_count][c2]=="forward" and sys_loop_map[c1][c2]=="no":
	sys_loop_map[c1][c2]="reverse"
	elif sys_loop_map[map_loop_count][c2]=="stiff_forward" and sys_loop_map[c1][c2]=="stiff_forward":
	sys_loop_map[c1][c2]="no"
	elif sys_loop_map[map_loop_count][c2]=="stiff_reverse" and sys_loop_map[c1][c2]=="stiff_reverse":
	sys_loop_map[c1][c2]="no"
	elif sys_loop_map[map_loop_count][c2]=="stiff_reverse" and sys_loop_map[c1][c2]=="no":
	sys_loop_map[c1][c2]="stiff_forward"
	elif sys_loop_map[map_loop_count][c2]=="stiff_forward" and sys_loop_map[c1][c2]=="no":
	sys_loop_map[c1][c2]="stiff_reverse"
	else:
	for c2 in range(len(sys_loop_map[0])):
	if sys_loop_map[map_loop_count][c2]=="forward" and sys_loop_map[c1][c2]=="reverse":
	sys_loop_map[c1][c2]="no"
	elif sys_loop_map[map_loop_count][c2]=="reverse" and sys_loop_map[c1][c2]=="forward":
	sys_loop_map[c1][c2]="no"
	elif sys_loop_map[map_loop_count][c2]=="reverse" and sys_loop_map[c1][c2]=="no":
	sys_loop_map[c1][c2]="reverse"
	elif sys_loop_map[map_loop_count][c2]=="forward" and sys_loop_map[c1][c2]=="no":
	sys_loop_map[c1][c2]="forward"
	elif sys_loop_map[map_loop_count][c2]=="stiff_forward" and sys_loop_map[c1][c2]=="stiff_reverse":
	sys_loop_map[c1][c2]="no"
	elif sys_loop_map[map_loop_count][c2]=="stiff_reverse" and sys_loop_map[c1][c2]=="stiff_forward":
	sys_loop_map[c1][c2]="no"
	elif sys_loop_map[map_loop_count][c2]=="stiff_reverse" and sys_loop_map[c1][c2]=="no":
	sys_loop_map[c1][c2]="stiff_reverse"
	elif sys_loop_map[map_loop_count][c2]=="stiff_forward" and sys_loop_map[c1][c2]=="no":
	sys_loop_map[c1][c2]="stiff_forward"
	# Increment the loop count and make the branch count equal that
	# So the default is a diagonal stiff element.
	map_loop_count+=1
	map_branch_count=map_loop_count
	# This is to make sure, that stiff loops
	# are connected to an input.
	for c1 in range(len(sys_loop_map)):
	is_loop_stiff="no"
	has_input="no"
	for c3 in range(len(branch_info)):
	#if branch_info[c3][:-1]==sys_loops[c1][c1][c2][:-1]:
	if sys_loop_map[c1][c3]=="stiff_forward" or sys_loop_map[c1][c3]=="stiff_reverse":
	is_loop_stiff="yes"
	# Check if any branch has a non-zero B matrix entry.
	for c4 in range(len(branch_info[c3][-1][1])):
	if branch_info[c3][-1][1][c4]:
	has_input="yes"
	# If loop is stiff and has no input
	if is_loop_stiff=="yes" and has_input=="no":
	c2=0
	flag_input="no"
	# Check all other loops whether they are stiff
	# and whether they have inputs.
	while flag_input=="no" and c2<len(sys_loop_map):
	if not c1==c2:
	is_loop_stiff="no"
	flag_input="no"
	for c4 in range(len(branch_info)):
	if sys_loop_map[c2][c4]=="stiff_forward" or sys_loop_map[c2][c4]=="stiff_reverse":
	is_loop_stiff="yes"
	if is_loop_stiff=="no":
	for c5 in range(len(branch_info[c4][-1][1])):
	if branch_info[c4][-1][1][c5]:
	flag_input="yes"
	c2=c2+1
	# Perform a row operation with a loop that is non-stiff
	# and has an input.
	if is_loop_stiff=="no" and flag_input=="yes":
	c2=c2-1
	#loop_manipulate(sys_loops, c1, c2, "add")
	# Update the sys_loop_map info
	for c3 in range(len(branch_info)):
	if sys_loop_map[c1][c3]=="forward" and sys_loop_map[c2][c3]=="reverse":
	sys_loop_map[c1][c3]="no"
	elif sys_loop_map[c1][c3]=="reverse" and sys_loop_map[c2][c3]=="forward":
	sys_loop_map[c1][c3]="no"
	elif sys_loop_map[c1][c3]=="forward" and sys_loop_map[c2][c3]=="no":
	sys_loop_map[c1][c3]="forward"
	elif sys_loop_map[c1][c3]=="reverse" and sys_loop_map[c2][c3]=="no":
	sys_loop_map[c1][c3]="reverse"
	elif sys_loop_map[c1][c3]=="stiff_forward" and sys_loop_map[c2][c3]=="stiff_reverse":
	sys_loop_map[c1][c3]="no"
	elif sys_loop_map[c1][c3]=="stiff_reverse" and sys_loop_map[c2][c3]=="stiff_forward":
	sys_loop_map[c1][c3]="no"
	elif sys_loop_map[c1][c3]=="stiff_forward" and sys_loop_map[c2][c3]=="no":
	sys_loop_map[c1][c3]="stiff_forward"
	elif sys_loop_map[c1][c3]=="stiff_reverse" and sys_loop_map[c2][c3]=="no":
	sys_loop_map[c1][c3]="stiff_reverse"
	# Get rid of repeating loops by setting all branches "no"
	for c1 in range(len(sys_loop_map)-1):
	for c2 in range(c1+1, len(sys_loop_map)):
	loop_repeats="yes"
	for c3 in range(len(sys_loop_map[c1])):
	if not sys_loop_map[c1][c3]==sys_loop_map[c2][c3]:
	loop_repeats="no"
	if loop_repeats=="yes":
	for c3 in range(len(sys_loop_map[c2])):
	sys_loop_map[c2][c3]="no"
	# Re-initialize the matrices A, B, and E
	matrix_a.zeros(len(sys_loops),len(sys_loops))
	matrix_e.zeros(len(sys_loops),len(sys_loops))
	matrix_b.zeros(len(sys_loops),matrix_b.columns)
	# Recalculate the matrices for the new system loops.
	recalculate_sys_matrices(sys_loop_map, branch_info, matrix_a, matrix_e, matrix_b)
	for c1 in range(matrix_a.rows):
	for c2 in range(matrix_a.columns):
	if abs(matrix_a.data[c1][c2])<1.0e-10:
	matrix_a.data[c1][c2]=0.0
	for c1 in range(matrix_e.rows):
	for c2 in range(matrix_e.columns):
	if abs(matrix_e.data[c1][c2])<1.0e-10:
	matrix_e.data[c1][c2]=0.0
	for c1 in range(matrix_b.rows):
	for c2 in range(matrix_b.columns):
	if abs(matrix_b.data[c1][c2])<1.0e-10:
	matrix_b.data[c1][c2]=0.0
	return

view raw remove_stiffness.py hosted with ❤ by GitHub

And the resultant code is much cleaner, the loop manipulations much easier. This will also change the code in the circuit elements (click on "view raw" below the code box to see the code in another window):

	def update_val(self, sys_loop_map, lbyr_ratio, mat_e, mat_a, mat_b, state_vec, mat_u, sys_branches, sys_events):
	""" This function calculates the actual current in the
	diode branch. With this, the branch voltage is found
	with respect to the existing diode resistance. The diode
	voltage is then used to decide the turn on condition. """
	for c1 in range(len(sys_branches)):
	if csv_tuple(self.pos) in sys_branches[c1]:
	branch_pos=c1
	self.current=0.0
	for c1 in range(len(sys_loop_map)):
	# If diode cathode is after the voltmeter
	# position, it means current is positive.
	if sys_branches[branch_pos].index(self.polrty)>sys_branches[branch_pos].index(csv_tuple(self.pos)):
	# Then check is the loop is aiding or opposing
	# the main loop.
	if sys_loop_map[c1][branch_pos]=="forward" or sys_loop_map[c1][branch_pos]=="stiff_forward":
	self.current+=state_vec.data[c1][0]
	elif sys_loop_map[c1][branch_pos]=="reverse" or sys_loop_map[c1][branch_pos]=="stiff_reverse":
	self.current-=state_vec.data[c1][0]
	else:
	if sys_loop_map[c1][branch_pos]=="forward" or sys_loop_map[c1][branch_pos]=="stiff_forward":
	self.current-=state_vec.data[c1][0]
	elif sys_loop_map[c1][branch_pos]=="reverse" or sys_loop_map[c1][branch_pos]=="stiff_reverse":
	self.current+=state_vec.data[c1][0]
	self.voltage=self.current*self.resistor
	# Diode will turn on when it is forward biased
	# and it was previously off.
	if self.status=="off" and self.voltage>1.0:
	sys_events[branch_pos]="yes"
	self.status="on"
	# Diode will turn off only when current becomes
	# negative.
	if self.status=="on" and self.current<-1.0e-5:
	sys_events[branch_pos]="yes"
	self.status="off"
	if self.status=="off":
	self.resistor=self.resistor_off
	else:
	self.resistor=self.resistor_on
	return

view raw update_value.py hosted with ❤ by GitHub