<reponame>timtyree/bgmc<filename>python/lib/lib_care/routines/generate_tip_logs.py
#! /usr/bin/env python
# Generate Birth Death Rates for Given Initial Conditions
# Method = Find intersections of contours after computing entire contours for each observation (unavoidably slow).
# <NAME>
# 10.28.2020
from .. import * 
from ..measure.measure import *  
#^this import is not the best practice. however, minimalism reigns the pythonic.
import numpy as np, pandas as pd, matplotlib.pyplot as plt

#automate the boring stuff
# from IPython import utils
import time, os, sys, re
beep = lambda x: os.system("echo -n '\\a';sleep 0.2;" * x)
if not 'here_dir' in globals():
	here_dir = os.getcwd()
def generate_tip_logs_from_ic(initial_condition_dir, h, tmax,
	jump_threshold, size_threshold, V_threshold, pad, decimals,
	tmin_early_stopping, save_every_n_frames, round_output_decimals, 
	timing, printing, logging, asserting, beeping, saving, 
	data_dir_log, completed_ic_dir, print_log_dir, **kwargs):
	'''generates a log of tip locations on 2D grid with periodic boundary conditions.
	default key word arguments are returned by lib.routines.kwargs.get_kwargs(initial_condition_dir).'''
	
	# if logging, change the print statements to a .log file unique to ic
	if logging:
		log = open(print_log_dir, "a")
		sys.stdout = log

	if printing:
		print(f'loading initial conditions from: \n\t{initial_condition_dir}.')

	# os.chdir(here_dir)
	txt = load_buffer(initial_condition_dir)
	width, height = txt.shape[:2]
	zero_txt = txt.copy()*0.
	# width, height, channel_no = txt.shape

	kwargs.update({
		'width':width,
		'height':height
		})
	#reinitialize records
	time_start = 0.  #eval(buffer_fn[buffer_fn.find('time_')+len('time_'):-4])
	if asserting:
		assert (float(time_start) is not None)
	tip_state_lst = []
	t = time_start
	dict_out_lst = []  
	num_steps = int(np.around((tmax-t)/h))

	#precompute anything that needs precomputing
	if printing:
		#print(f"sigma is {sigma}, threshold is {threshold}.")
		#print(f"pad is {pad}, rejection_distance is edge_tolerance is {edge_tolerance}.")
		print(f"starting simulation.  integrating no further than time {tmax:.3f} milliseconds.")
	if timing:
		start = time.time()
	##########################################
	#run the simulation, measuring regularly
	##########################################
	step_count = 0
	n_tips=1 #to initialize loop invarient (n_tips > 0) to True
	while (t<tmax) & (n_tips > 0):
		if step_count%save_every_n_frames != 0:
			#forward Euler integration in time
			time_step(txt, h, zero_txt)
		else:
			#take measurements once every n frames
			#compute as discrete flow map dtexture_dt
			dtexture_dt = zero_txt.copy()
			get_time_step(txt, dtexture_dt)

			#compute the images to find isosurfaces of
			img    = txt[...,0]
			dimgdt = dtexture_dt[...,0]

			#compute both families of contours
			contours1 = find_contours(img,    level = 0.8)
			contours2 = find_contours(dimgdt, level = 0.0)

			#find_tips and measure tip topological/EP state
			s1_list, s2_list, x_lst, y_lst, v_lst, f_lst, s_lst = measure_system(contours1, contours2, width, height, txt, 
																				 jump_threshold = jump_threshold, 
																				 size_threshold = size_threshold, 
																				 pad=pad, decimals=decimals)
			n_tips = x_lst.size
			dict_out = {
				't': float(t),
				'n': int(n_tips),
				'x': tuple(x_lst),
				'y': tuple(y_lst),
				'n1': tuple(s1_list),
				'n2': tuple(s2_list),
				'v':v_lst,
				'f':f_lst,
				's':s_lst,
			}
			
			#record data for current time
			dict_out_lst.append(dict_out)

			#forward Euler integration in time
			txt += h*dtexture_dt

			#update progress bar after each measurement
			if not logging:
				if printing:
					printProgressBar(step_count, num_steps, prefix = 'Progress:', suffix = 'Complete', length = 50)
		#advance time by one step
		t   += h
		step_count += 1
	if not logging:
		if printing:
			printProgressBar(step_count, num_steps, prefix = 'Progress:', suffix = 'Complete', length = 50)

	if printing:
		#report the bottom line up front
		if timing:
			print(f"\ntime integration complete. run time was {time.time()-start:.2f} seconds in realtime")
		print(f"\ncurrent time is {t:.1f} ms in simulation time.")

		print(f"number of nan pixel voltages is {np.max(sum(np.isnan(txt[...,0])))}.")
		# print(f"current max voltage is {np.nanmax(txt[...,0]):.4f}.")
		# print(f"current max fast variable is {np.nanmax(txt[...,1]):.4f}.")
		# print(f"current max slow variable is {np.nanmax(txt[...,2]):.4f}.")
		print(f"number of tips is = {n_tips}.") 
		if n_tips==0:
			print(f"zero tips remaining at time t = {t:.1f} ms.")
	if beeping:
		beep(1)
	if printing:
		if t >= tmax:
			print( f"Caution! max_time was reached! Termination time not reached!  Consider rerunning with greater max_time!")
	if saving:
		df = pd.concat([pd.DataFrame(dict_out) for dict_out in dict_out_lst])
		df.reset_index(inplace=True, drop=True)
		#if the end of AF was indeed reachded, append a row recording this
		if n_tips==0:
			next_id = df.index.values[-1]+1
			df = pd.concat([df,pd.DataFrame({'t': float(save_every_n_frames*h+t),'n': int(n_tips)}, index = [next_id])])
		#save the recorded data
		df.round(round_output_decimals).to_csv(data_dir_log, index=False)
		if printing:
			print('saved to:')
			print(data_dir_log)

	#move the completed file to ic-out
	os.rename(initial_condition_dir,completed_ic_dir)
	#input ic moved to output
	if logging:
		if not log.closed:
			log.close()		

	return kwargs

if __name__=='__main__':
	for ic in sys.argv[1:]:
		kwargs = get_kwargs(ic)
		kwargs = generate_tip_logs_from_ic(ic, **kwargs)
		print(f"completed birth_death_rates_from_ic: {ic}")
		print(f"csv of spiral tip data stored in: {kwargs['completed_ic_dir']}")

#! /usr/bin/env python
# Generate Birth Death Rates for Given Initial Conditions
# Method = Find intersections of contours after computing entire contours for each observation (unavoidably slow).
# <NAME>
# 10.28.2020
from .. import * 
from ..measure.measure import *  
#^this import is not the best practice. however, minimalism reigns the pythonic.
import numpy as np, pandas as pd, matplotlib.pyplot as plt

#automate the boring stuff
# from IPython import utils
import time, os, sys, re
beep = lambda x: os.system("echo -n '\\a';sleep 0.2;" * x)
if not 'here_dir' in globals():
	here_dir = os.getcwd()
def generate_tip_logs_from_ic(initial_condition_dir, h, tmax,
	jump_threshold, size_threshold, V_threshold, pad, decimals,
	tmin_early_stopping, save_every_n_frames, round_output_decimals, 
	timing, printing, logging, asserting, beeping, saving, 
	data_dir_log, completed_ic_dir, print_log_dir, **kwargs):
	'''generates a log of tip locations on 2D grid with periodic boundary conditions.
	default key word arguments are returned by lib.routines.kwargs.get_kwargs(initial_condition_dir).'''
	
	# if logging, change the print statements to a .log file unique to ic
	if logging:
		log = open(print_log_dir, "a")
		sys.stdout = log

	if printing:
		print(f'loading initial conditions from: \n\t{initial_condition_dir}.')

	# os.chdir(here_dir)
	txt = load_buffer(initial_condition_dir)
	width, height = txt.shape[:2]
	zero_txt = txt.copy()*0.
	# width, height, channel_no = txt.shape

	kwargs.update({
		'width':width,
		'height':height
		})
	#reinitialize records
	time_start = 0.  #eval(buffer_fn[buffer_fn.find('time_')+len('time_'):-4])
	if asserting:
		assert (float(time_start) is not None)
	tip_state_lst = []
	t = time_start
	dict_out_lst = []  
	num_steps = int(np.around((tmax-t)/h))

	#precompute anything that needs precomputing
	if printing:
		#print(f"sigma is {sigma}, threshold is {threshold}.")
		#print(f"pad is {pad}, rejection_distance is edge_tolerance is {edge_tolerance}.")
		print(f"starting simulation.  integrating no further than time {tmax:.3f} milliseconds.")
	if timing:
		start = time.time()
	##########################################
	#run the simulation, measuring regularly
	##########################################
	step_count = 0
	n_tips=1 #to initialize loop invarient (n_tips > 0) to True
	while (t<tmax) & (n_tips > 0):
		if step_count%save_every_n_frames != 0:
			#forward Euler integration in time
			time_step(txt, h, zero_txt)
		else:
			#take measurements once every n frames
			#compute as discrete flow map dtexture_dt
			dtexture_dt = zero_txt.copy()
			get_time_step(txt, dtexture_dt)

			#compute the images to find isosurfaces of
			img    = txt[...,0]
			dimgdt = dtexture_dt[...,0]

			#compute both families of contours
			contours1 = find_contours(img,    level = 0.8)
			contours2 = find_contours(dimgdt, level = 0.0)

			#find_tips and measure tip topological/EP state
			s1_list, s2_list, x_lst, y_lst, v_lst, f_lst, s_lst = measure_system(contours1, contours2, width, height, txt, 
																				 jump_threshold = jump_threshold, 
																				 size_threshold = size_threshold, 
																				 pad=pad, decimals=decimals)
			n_tips = x_lst.size
			dict_out = {
				't': float(t),
				'n': int(n_tips),
				'x': tuple(x_lst),
				'y': tuple(y_lst),
				'n1': tuple(s1_list),
				'n2': tuple(s2_list),
				'v':v_lst,
				'f':f_lst,
				's':s_lst,
			}
			
			#record data for current time
			dict_out_lst.append(dict_out)

			#forward Euler integration in time
			txt += h*dtexture_dt

			#update progress bar after each measurement
			if not logging:
				if printing:
					printProgressBar(step_count, num_steps, prefix = 'Progress:', suffix = 'Complete', length = 50)
		#advance time by one step
		t   += h
		step_count += 1
	if not logging:
		if printing:
			printProgressBar(step_count, num_steps, prefix = 'Progress:', suffix = 'Complete', length = 50)

	if printing:
		#report the bottom line up front
		if timing:
			print(f"\ntime integration complete. run time was {time.time()-start:.2f} seconds in realtime")
		print(f"\ncurrent time is {t:.1f} ms in simulation time.")

		print(f"number of nan pixel voltages is {np.max(sum(np.isnan(txt[...,0])))}.")
		# print(f"current max voltage is {np.nanmax(txt[...,0]):.4f}.")
		# print(f"current max fast variable is {np.nanmax(txt[...,1]):.4f}.")
		# print(f"current max slow variable is {np.nanmax(txt[...,2]):.4f}.")
		print(f"number of tips is = {n_tips}.") 
		if n_tips==0:
			print(f"zero tips remaining at time t = {t:.1f} ms.")
	if beeping:
		beep(1)
	if printing:
		if t >= tmax:
			print( f"Caution! max_time was reached! Termination time not reached!  Consider rerunning with greater max_time!")
	if saving:
		df = pd.concat([pd.DataFrame(dict_out) for dict_out in dict_out_lst])
		df.reset_index(inplace=True, drop=True)
		#if the end of AF was indeed reachded, append a row recording this
		if n_tips==0:
			next_id = df.index.values[-1]+1
			df = pd.concat([df,pd.DataFrame({'t': float(save_every_n_frames*h+t),'n': int(n_tips)}, index = [next_id])])
		#save the recorded data
		df.round(round_output_decimals).to_csv(data_dir_log, index=False)
		if printing:
			print('saved to:')
			print(data_dir_log)

	#move the completed file to ic-out
	os.rename(initial_condition_dir,completed_ic_dir)
	#input ic moved to output
	if logging:
		if not log.closed:
			log.close()		

	return kwargs

if __name__=='__main__':
	for ic in sys.argv[1:]:
		kwargs = get_kwargs(ic)
		kwargs = generate_tip_logs_from_ic(ic, **kwargs)
		print(f"completed birth_death_rates_from_ic: {ic}")
		print(f"csv of spiral tip data stored in: {kwargs['completed_ic_dir']}")
#! /usr/bin/env python # Generate Birth Death Rates for Given Initial Conditions # Method = Find intersections of contours after computing entire contours for each observation (unavoidably slow). # <NAME> # 10.28.2020 #^this import is not the best practice. however, minimalism reigns the pythonic. #automate the boring stuff # from IPython import utils generates a log of tip locations on 2D grid with periodic boundary conditions. 	default key word arguments are returned by lib.routines.kwargs.get_kwargs(initial_condition_dir). # if logging, change the print statements to a .log file unique to ic # os.chdir(here_dir) # width, height, channel_no = txt.shape #reinitialize records #eval(buffer_fn[buffer_fn.find('time_')+len('time_'):-4]) #precompute anything that needs precomputing #print(f"sigma is {sigma}, threshold is {threshold}.") #print(f"pad is {pad}, rejection_distance is edge_tolerance is {edge_tolerance}.") ########################################## #run the simulation, measuring regularly ########################################## #to initialize loop invarient (n_tips > 0) to True #forward Euler integration in time #take measurements once every n frames #compute as discrete flow map dtexture_dt #compute the images to find isosurfaces of #compute both families of contours #find_tips and measure tip topological/EP state #record data for current time #forward Euler integration in time #update progress bar after each measurement #advance time by one step #report the bottom line up front # print(f"current max voltage is {np.nanmax(txt[...,0]):.4f}.") # print(f"current max fast variable is {np.nanmax(txt[...,1]):.4f}.") # print(f"current max slow variable is {np.nanmax(txt[...,2]):.4f}.") #if the end of AF was indeed reachded, append a row recording this #save the recorded data #move the completed file to ic-out #input ic moved to output