Homec4science
Diffusion G(roup)-RENAISSANCE (master) Browse

grenaissance/docker/work/1-renaissance.py
85021e462c8amaster

Last Change
Raw File

1-renaissance.py
View Options

# Native Python libraries
import sys
import os
import stat
import math
import time
from datetime import date
from argparse import ArgumentParser
from configparser import ConfigParser
# Public Python libraries
import numpy as np
import mpire as mp
# Custom Python libraries
import helper as hp
from evostrat.init_mlp import MLP
from evostrat.evolution_strategy import EvolutionStrategy
from kinetics.jacobian_solver import CheckJacobian
### Global functions for multiprocessing purpose.
def mp_singular_renaissance(shared_objects, state):
"""Parallelize Renaissance training loop over singular steady-states for singular state optimization
:param shared_objects (Ren aissance class): Renaissance object
:param state (string): state on which to parallelize the training loop for singular optimization of generators
"""
renaissance_object = shared_objects
renaissance_object.train(state)
return
def mp_jacobian(shared_objects, state):
"""Parallelize Renaissance reward function over steady-states Jacobian objects
:param shared_objects: shared parameters over multiprocessing workers (number of workers for this thread,
Jacobian dict, generated sets of kinetic parameters)
:param state: steady-state of this thread
:return: list of steady-state's Jacobian maximum (real-part) eigenvalues
"""
n_thread, jacobian_dict, generated_params = shared_objects
jacobian = jacobian_dict[state]
with mp.WorkerPool(n_jobs=n_thread, shared_objects=jacobian) as pool:
results = pool.map(mp_max_eigenvalues, generated_params)
return results
def mp_max_eigenvalues(jacobian, generated_params):
"""Parallelize max (real-part) eigenvalues computation for given model's Jacobian and set of generated params.
:param jacobian: model's Jacobian object
:param generated_params: list of sets of generated kinetic parameters to be parallelized over workers
:return: list of maximum (real-part) eigenvalues
"""
results = jacobian.calc_eigenvalues_recal_vmax(generated_params)
return results
def mp_fitness(renaissance_object, maximum_eigenvalues):
"""Parallelize fitness of steady-state given max (real-part) eigenvalues of model's Jacobian.
:param renaissance_object (Renaissance): Renaissance class object
:param maximum_eigenvalues (list): maximum (real-part) eigenvalues of model's Jacobian
:return: fitness score
"""
results = renaissance_object._fitness(maximum_eigenvalues)
return results
class RenParser(ArgumentParser):
def __init__(self, description):
ArgumentParser.__init__(self, description)
self.add_argument('-V', '--verbose', action='store_true', default=False,
help="Verbose mode, if used, then some additional prints will be displayed about"
"time of processes.",
dest='verbose')
self.add_argument('-S', '--states', nargs='+', default=None, type=int,
help="Steady-states list on which to run Renaissance. If not specified, read it from"
"configfile.ini.", dest='states')
self.add_argument('-M', '--mode', nargs='?', default='group', type=str,
choices=['group', 'singular', 'both'],
help="Optimization mode of Renaissance;\ngroup: consider the given steady-states list as "
"one group.\nsingular: consider the given steady-states list as singular steady-states."
"\nboth: do both optimizations.", dest='mode')
# renaissance class: allows to run an instance of ES for generating kinetic parameters
class Renaissance:
"""Class encapsulating all necessary functions to optimize generator weights using Natural Evolution Strategies
in order to generate dynamically relevant Michaelis constants (Kms) relatively to given (group of) steady-states.
"""
def __init__(self, config_file, optimization, steady_states=None, verbose=False):
configs = ConfigParser()
configs.read(config_file)
self.optimization = optimization
self.verbose = verbose
self.n_samples = int(configs['MLP']['n_samples'])
self.latent_dim = int(configs['MLP']['latent_dim'])
self.n_kms = int(configs['MLP']['no_kms'])
self.layers = [int(configs['MLP']['layer_1']), int(configs['MLP']['layer_2']),
int(configs['MLP']['layer_3'])]
self.lnminkm = float(configs['CONSTRAINTS']['min_km'])
self.lnmaxkm = float(configs['CONSTRAINTS']['max_km'])
self.repeats = int(configs['EVOSTRAT']['repeats'])
self.save_step = int(configs['EVOSTRAT']['save_step'])
self.generations = int(configs['EVOSTRAT']['generations'])
self.pop_size = int(configs['EVOSTRAT']['pop_size'])
self.noise = float(configs['EVOSTRAT']['noise'])
self.lr = float(configs['EVOSTRAT']['lr'])
self.decay = float(configs['EVOSTRAT']['decay'])
if steady_states:
self.steady_states = steady_states
else:
self.steady_states = [int(idx) for idx in configs['EVOSTRAT']['ss_idx'].split(',')]
print(f"Steady-states specified in configfile.ini are {self.steady_states}.\n")
self.n_threads = int(configs['EVOSTRAT']['n_threads'])
self.print_step = int(configs['EVOSTRAT']['print_step'])
self.output_path = configs['PATHS']['output_path']
self.met_model = configs['PATHS']['met_model']
self.names_km = hp.load_pkl(f'models/{self.met_model}/parameter_names_km_fdp1.pkl')
self.reward_flag = int(configs['REWARDS']['reward_flag'])
self.eig_partition = float(configs['REWARDS']['eig_partition'])
self.n_consider = int(configs['REWARDS']['n_consider'])
self.pf_flag = int(configs['PARAMETER_FIXING']['pf_flag'])
self.fixed_parameter_names = configs['PARAMETER_FIXING']['fixed_parameter_names'] if self.pf_flag else None
self.fixed_parameter_ranges = configs['PARAMETER_FIXING']['fixed_parameter_ranges'] if self.pf_flag else None
self.path_to_weights = configs['LOAD_AND_GENERATE']['path_to_weights']
self.save_path = self._save_parameters()
# Initialize a CheckJacobian objects per steady-state and store it in a dictionnary for multiprocessing purpose
self.jacobian_dict = {}
for state in self.steady_states:
jacobian = CheckJacobian(n_cpu=self.n_threads,
met_model=self.met_model,
exp_id='fdp1',
steady_states=state,
names_km=self.names_km)
self.jacobian_dict.update({state: jacobian})
print(f"Dict of {len(self.steady_states)} {CheckJacobian.__name__} objects"
f" created for multiprocessing purpose\n"
f"Total size used: {round(sys.getsizeof(self.jacobian_dict) / 1e3, 3)} kb")
def _fitness(self, max_eigenvalues):
"""Depending on the value of self.reward_flag (defined in configfile.ini), return the fitness of the given
set of maximum eigenvalues.
Keyword arguments:
:param max_eigenvalues: numpy array NxM of maximum eigenvalues;
N steady-states and M generated sets of kinetic parameters
:return fitness: fitness score
"""
# Make sure array is NxM, N sets of generated params and M steady-states
if max_eigenvalues.ndim == 1:
max_eigenvalues = max_eigenvalues.reshape((1, -1))
# Reward for generation of models with fast dynamics
if self.reward_flag == 0:
max_neg_eig = max_eigenvalues.min(axis=1)
if np.all(max_neg_eig > self.eig_partition):
fitness = 0.01 / (1 + np.exp(np.mean(max_neg_eig) - self.eig_partition))
else:
fitness = len(np.where(max_eigenvalues <= self.eig_partition)[0]) / max_eigenvalues.size
# Reward for maximum incidence of valid models
elif self.reward_flag == 1:
max_eigenvalues.sort(axis=1)
ss_considered_avg = [max_eig[:self.n_consider].mean() for max_eig in max_eigenvalues]
fitness = np.exp(-0.1 * np.mean(ss_considered_avg)) / 2
else:
fitness = 0
return fitness
def _rewards(self, generated_params, state=None, separate=False):
"""Compute the rewards in a parallel fashion (if group of steady-states):
parameterize the generated params to each Jacobian objects and compute the maximum (real-part) eigenvalues,
then call _fitness(...) to compute the fitness score of the generator.
:param generated_params: generated sets of kinetic parameters
:param state: if the optimization mode is currently 'singular', then indicates which steady-state to use
:param separate: boolean flag indicating if the method should return a separate list for each steady-state
:return: dict of rewards for the given weights (and number of steady-states if separate is True)
"""
# Make the generated parameters iterables for multiprocessing purpose
generated_params_iters = [[params] for params in generated_params]
# If current optimization mode is 'group' and more than 1 steady-state was given
if len(self.steady_states) > 1 and \
(self.optimization == 'group' or (self.optimization == 'both' and state is None)):
n_thread_pool = math.floor((self.n_threads - len(self.steady_states)) / len(self.steady_states))
with mp.WorkerPool(n_jobs=len(self.steady_states), daemon=False,
shared_objects=(n_thread_pool, self.jacobian_dict, generated_params_iters)) as pool:
max_eig = pool.map(mp_jacobian, self.steady_states)
# 1 steady-state given or current optimization mode is 'singular'
else:
if state is None:
state = self.steady_states[0]
n_thread_pool = math.floor((self.n_threads - len(self.steady_states)) / len(self.steady_states)) \
if self.optimization == 'singular' else self.n_threads
with mp.WorkerPool(n_jobs=n_thread_pool, shared_objects=self.jacobian_dict[state]) as pool:
max_eig = pool.map(mp_max_eigenvalues, generated_params_iters)
max_eig = np.array(max_eig).reshape((-1, self.n_samples))
# Compute the reward of each individual steady-state separately
if separate and len(self.steady_states) > 1 and \
(self.optimization == 'group' or (self.optimization == 'both' and state is None)):
with mp.WorkerPool(n_jobs=len(self.steady_states), shared_objects=self) as pool:
this_reward = pool.map(
mp_fitness,
[eig for eig in max_eig])
this_reward = dict(map(lambda i: (self.steady_states[i], this_reward[i]), range(len(this_reward))))
this_reward.update({'group': self._fitness(max_eig)})
else:
this_reward = {'group': self._fitness(max_eig)}
return this_reward
def _initialize_generator(self, mlp):
""" Initialize generator with custom weights if there exists some, else use default initializer.
:param mlp: Multilayer Perceptron generator to be initialized
"""
if self.path_to_weights:
weights = hp.load_pkl(self.path_to_weights)
mlp.generator.set_weights(weights)
print("MLP has been initialized with custom weights")
else:
print("Path to custom weights was empty.\nMLP weights are initialized with Keras Glorot-uniform ",
"initializer. See https://www.tensorflow.org/api_docs/python/tf/keras/initializers/GlorotUniform "
"for more information.")
mlp.reset_weights()
def _save_parameters(self, state=None):
"""Save to disk the run parameters.
:param state: given if current optimization mode is 'singular'
:return:
"""
# Create outputs directory
today = date.today().strftime("%Y-%m-%d")
# Six random digits identification number for the run
if state:
return str(state)
else:
seed = int(str(time.time())[:-10])
np.random.seed(seed)
this_run_name = str(np.random.randint(0, (10 ** 6) - 1))
if len(this_run_name) < 6:
this_run_name = (6 - len(this_run_name)) * '0' + this_run_name
mode = self.optimization.upper()[0]
save_path = f'{self.output_path}{today}_{mode}-SS{len(self.steady_states)}RF{self.reward_flag}' \
f'PF{self.pf_flag}_{this_run_name}/'
os.makedirs(save_path, exist_ok=True)
# Save config parameters
config = self.__dict__.copy()
try:
del config['jacobian_dict']
except KeyError:
pass
hp.save_pkl(save_path + 'config', config)
with open(save_path + 'steady_states.txt', 'wt') as f:
if (self.optimization == 'singular' or self.optimization == 'both') and state is not None:
f.write(str(state) + '\n')
else:
for state in self.steady_states:
f.write(str(state) + '\n')
return save_path
def run(self):
"""Start the optimization run, parallelize runs if current optimization mode is 'singular'.
"""
if self.optimization == 'group':
print(f"\ngroup optimization starting...\n")
self.train()
elif self.optimization == 'singular':
print(f"\nSingular optimization starting...\n")
with mp.WorkerPool(n_jobs=len(self.steady_states), daemon=False, shared_objects=self) as pool:
pool.map(mp_singular_renaissance, self.steady_states)
elif self.optimization == 'both':
# group optimization first, then parallel optimization of singular states inside the group
print(f"\ngroup optimization starting...\n")
self.train()
print(f"\nSingular optimization starting...\n")
with mp.WorkerPool(n_jobs=len(self.steady_states), daemon=False, shared_objects=self) as pool:
pool.map(mp_singular_renaissance, self.steady_states)
else:
raise ValueError(f"{self.optimization} is not a correct value, use -h/--help flag to see correct values.")
def train(self, state=None):
"""Training loop of generator optimization.
:param state: given if current optimization run is 'singular'
"""
if state:
save_path = os.path.join(self.save_path, self._save_parameters(state))
os.makedirs(save_path, exist_ok=True)
else:
save_path = self.save_path
# Initializing MLP object
mlp = MLP(n_samples=self.n_samples,
latent_dim=self.latent_dim,
layers=self.layers,
n_kms=self.n_kms,
cond_class=1,
min_x=self.lnminkm,
max_x=self.lnmaxkm,
names_km=self.names_km,
param_fixing=self.pf_flag,
path_fixed_names=self.fixed_parameter_names,
path_fixed_ranges=self.fixed_parameter_ranges)
first = True
try:
for rep in range(self.repeats):
rep_save_path = os.path.join(save_path, f'repeat_{rep}/')
os.makedirs(rep_save_path, exist_ok=True)
print(f"\n//////////////////////// CROSS-VALIDATION FOLD {rep + 1} ////////////////////////")
# Re-initialize generator weights if start of a new cross-validation fold
if first:
first = False
else:
self._initialize_generator(mlp)
es = EvolutionStrategy(mlp,
self._rewards,
rep_save_path,
self.path_to_weights,
self.save_step,
population_size=self.pop_size,
sigma=self.noise, # noise std deviation
learning_rate=self.lr,
decay=self.decay,
print_step=self.print_step,
verbose=self.verbose)
rewards = es.run(self.generations, state)
hp.save_pkl(f'{rep_save_path}/rewards', rewards)
except KeyboardInterrupt as ke:
print(ke)
finally:
# Give permissions to output directory
os.system("chmod -R 777 " + self.save_path)
print("Permissions changed for everyone to write, read and execute!")
if __name__ == "__main__":
parser = RenParser(description="Renaissance model: train MLP using Natural Selection Strategies to generate "
"efficiently valid kinetic parameters of metabolism.")
args = parser.parse_args(sys.argv[1:])
verbose = args.verbose
if not args.states:
states_info = f"Steady-states will be read from configfile.ini.\n"
else:
states_info = f"Steady-states are {args.states}.\n"
print(f"Optimization mode is {args.mode}.\n"
+ states_info
+ f"Verbose is set to {verbose}")
renaissance = Renaissance(config_file='configfile.ini',
optimization=args.mode,
steady_states=args.states,
verbose=verbose)
renaissance.run()
c4science · Help