diff --git a/examples/NONMEM/user/Example4/Cmaxppc.r b/examples/NONMEM/user/Example4/Cmaxppc.r index 05f0239..d62e7bc 100644 --- a/examples/NONMEM/user/Example4/Cmaxppc.r +++ b/examples/NONMEM/user/Example4/Cmaxppc.r @@ -2,7 +2,7 @@ library(dplyr) orgdata <- read.table("ORG.DAT",skip=1,header=TRUE) orgdata <- orgdata[orgdata$EVID==0&orgdata$TIME<=24,] # only first 24 hours simdata <- read.table("SIM.DAT",skip=1,header=TRUE) # only first 24 hours -simdata <- simdata[orgdata$EVID==0&orgdata$TIME<=24,] +simdata <- simdata[simdata$EVID==0&simdata$TIME<=24,] observed_maxes <- orgdata %>% group_by(ID) %>% summarise(max = max(DV)) sim_maxes <- simdata %>% group_by(ID) %>% diff --git a/src/darwin/ModelRun.py b/src/darwin/ModelRun.py index 4da95dc..6950cec 100644 --- a/src/darwin/ModelRun.py +++ b/src/darwin/ModelRun.py @@ -398,7 +398,7 @@ def run_model(self): if not failed: self.set_status('Finished model run') - if (options.isMOGA or self._post_run_r() and self._post_run_python()) and self._calc_fitness(): + if (self._post_run_r() and self._post_run_python()) and self._calc_fitness(): self.set_status('Done') def finish(self): diff --git a/src/darwin/algorithms/MOGA.py b/src/darwin/algorithms/MOGA.py index 224cefa..304b628 100644 --- a/src/darwin/algorithms/MOGA.py +++ b/src/darwin/algorithms/MOGA.py @@ -1,305 +1,188 @@ -import os -import shutil -import logging -import numpy as np -import warnings - -from darwin import Population -from darwin.Log import log -from darwin.options import options -import darwin.utils as utils -from darwin.ExecutionManager import keep_going -from darwin.ModelCode import ModelCode -from darwin.Template import Template -from darwin.ModelRun import ModelRun -from darwin.Population import Population -from darwin.ModelCache import get_model_cache -from darwin.ModelRunManager import rerun_models -from darwin.algorithms.run_downhill import do_moga_downhill_step -from darwin.ModelEngineAdapter import get_model_phenotype - -from .effect_limit import WeightedSampler - -from pymoo.algorithms.moo.nsga2 import NSGA2 -from pymoo.core.result import Result -from pymoo.core.population import pop_from_array_or_individual -from pymoo.operators.crossover.pntx import TwoPointCrossover -from pymoo.operators.mutation.bitflip import BitflipMutation -from pymoo.operators.sampling.rnd import BinaryRandomSampling -from pymoo.core.sampling import Sampling -from pymoo.core.problem import ElementwiseProblem - -warnings.filterwarnings('error', category=DeprecationWarning) -logger = logging.getLogger(__name__) - - -def _get_n_params(run: ModelRun) -> int: - model = run.model - - return model.estimated_omega_num + model.estimated_theta_num + model.estimated_sigma_num - - -class MogaProblem(ElementwiseProblem): - def __init__(self, n_var, run: ModelRun = None): - super(MogaProblem, self).__init__( - n_var=n_var, # number of bits - n_obj=2, - n_ieq_constr=0, - xl=np.zeros(n_var, dtype=int), - xu=np.ones(n_var, dtype=int), - # need this to send population and template to evaluate - requires_kwargs=True - ) - - self.run = run - - def _evaluate(self, x, out, *args, **kwargs): - f1 = self.run.result.ofv - f2 = 999999 if f1 >= options.crash_value else _get_n_params(self.run) - - out["F"] = [f1, f2] - - -def _get_front_runs(res: Result, template: Template, model_cache) -> list: - runs = [] - - maxes = template.gene_max - lengths = template.gene_length - - controls = {} - - for res_x in res.X: - cur_x = [int(x) for x in res_x.astype(int)] - mc = ModelCode.from_full_binary(cur_x, maxes, lengths) - run = model_cache.find_model_run(genotype=str(mc.IntCode)) \ - or model_cache.find_model_run(phenotype=get_model_phenotype(template, mc)) - - if run is None: - log.warn(f"Missing a front model: {cur_x}") - continue - - if run.control_file_name in controls: - controls[run.control_file_name] += 1 - continue - - controls[run.control_file_name] = 1 - - runs.append(run) - - runs = sorted(runs, key=lambda r: r.file_stem) - - for run in runs: - x = controls[run.control_file_name] - x = f"(x{x})" if x > 1 else '' - log.message(f"Model {run.control_file_name}, OFV = {run.result.ofv:.4f}, NEP = {_get_n_params(run)} {x}") - - reruns = [run for run in runs if not os.path.isdir(run.run_dir)] - - if reruns: - if options.rerun_front_models: - log.message('Re-running front models...') - - rerun_models(reruns) - else: - for run in reruns: - run.make_control_file() - run.output_results() - - return runs - - -def _front_str(runs: list) -> str: - lines = [f"OFV = {run.result.ofv}, NEP = {_get_n_params(run)}" for run in runs] - return '\n'.join(sorted(lines)) - - -class WeightedRandomSampling(Sampling): - def __init__(self, template: Template): - super(WeightedRandomSampling, self).__init__() - - self.sampler = WeightedSampler(template) - - def _do(self, problem, n_samples, **kwargs): - val = np.array([self.sampler.create_individual() for _ in range(n_samples)]) - return val.astype(bool) - - -class _MOGARunner: - def __init__(self, template: Template, pop_size: int): - self.n_var = sum(template.gene_length) - self.template = template - self.model_cache = get_model_cache() - self.pop = None - - problem = MogaProblem(n_var=self.n_var) - - self.algorithm = NSGA2( - pop_size=pop_size, - sampling=WeightedRandomSampling(template) if options.use_effect_limit else BinaryRandomSampling(), - crossover=TwoPointCrossover(prob=options.MOGA['crossover_rate']), - mutation=BitflipMutation(prob=options.MOGA['mutation_rate']), - eliminate_duplicates=True - ) - - self.algorithm.setup(problem, seed=options.random_seed, verbose=False) - - def has_next(self) -> bool: - return self.algorithm.has_next() - - def ask_population(self, n_gen: int) -> Population: - """ - Ask a population from the algorithm - It saves the moo population in self.pop for subsequent tell - - :param n_gen: current generation number - :return: Population - """ - self.pop = self.algorithm.ask() - - pop_full_bits = [[int(this_bit) for this_bit in this_ind.X] for this_ind in self.pop] - - return Population.from_codes(self.template, n_gen, pop_full_bits, ModelCode.from_full_binary) - - def tell_runs(self, runs: list) -> list: - """ - Tells the runs to the algorithm - It uses self.pop which is only saved by ask_population - The downhill step doesn't perform the _ask_, so pop is created from the runs - - :param runs: list of ModelRuns - :return: the front as a list of non-duplicated ModelRuns - """ - - pop = self.pop - - if pop is None: - infills = np.array([np.array(run.model.model_code.FullBinCode, dtype=bool) for run in runs]) - self.pop = pop = pop_from_array_or_individual(infills) - - for run, moo_ind in zip(runs, pop): - problem = MogaProblem(n_var=self.n_var, run=run) - - self.algorithm.evaluator.eval(problem, moo_ind) - - self.algorithm.tell(infills=pop) - - res = self.algorithm.result() - - self.pop = None - - log.message('Current Non Dominated models:') - - return _get_front_runs(res, self.template, self.model_cache) - - def run_moga_downhill(self, front: list, generation) -> tuple: - after = _front_str(front) - - changed = False - - for this_step in range(1, 100): # up to 99 steps - if not keep_going(): - break - - before = after - - downhill_runs = do_moga_downhill_step(self.template, front, generation, this_step) - - if not keep_going(): - break - - if not downhill_runs: - log.warn(f"Downhill step {generation}/{this_step} has nothing to add to the search, done with downhill") - break - - front = self.tell_runs(downhill_runs) - - after = _front_str(front) - - if before == after: - break - - changed = True - - if changed: - non_dominated_folder = os.path.join(options.non_dominated_models_dir, str(generation)) - - _copy_front_files(front, non_dominated_folder) - - return front, after - - def dump_res(self): - res = self.algorithm.result() - - log.message(f" MOGA best genome =\n{res.X.astype(int)},\n" - f" OFV and # of parameters =\n{res.F}") - - -def _copy_front_files(front: list, non_dominated_folder: str): - utils.remove_dir(non_dominated_folder) - os.mkdir(non_dominated_folder) - - for run in front: - if not os.path.isdir(run.run_dir): - log.warn('is not copied') - continue - - shutil.copytree(run.run_dir, os.path.join(non_dominated_folder, run.file_stem), dirs_exist_ok=True) - - -def run_moga(template: Template): - n_gens = options.num_generations - downhill_period = options.downhill_period - - n_gen = 0 - generations_no_change = 0 - after = '' - front = [] - - runner = _MOGARunner(template, options.population_size) - - while keep_going() and n_gen < n_gens: - if not runner.has_next(): - log.warn(f"MOGA finished before reaching generation {n_gens}") - break - - n_gen += 1 - - population = runner.ask_population(n_gen) - - population.run() - - if not keep_going(): - break - - before = after - - front = runner.tell_runs(population.runs) - - after = _front_str(front) - - non_dominated_folder = os.path.join(options.non_dominated_models_dir, str(n_gen)) - - _copy_front_files(front, non_dominated_folder) - - if downhill_period > 0 and n_gen % downhill_period == 0: - log.message(f"Starting downhill generation {n_gen}") - - front, after = runner.run_moga_downhill(front, n_gen) - - if before == after: - generations_no_change += 1 - log.message(f"No change in non dominated models for {generations_no_change} generations") - - if options.final_downhill_search and keep_going(): - log.message(f"Starting final downhill step") - front, after = runner.run_moga_downhill(front, 'FN') - - if not keep_going(): - return - - runner.dump_res() - - for run in front: - run.run_dir = options.output_dir - run.make_control_file(cleanup=False) - run.output_results() +from copy import copy +import time +import logging +import numpy as np +import random +import warnings + +import csv +import os + +from pymoo.core.problem import Problem +from pymoo.algorithms.moo.nsga3 import NSGA3 + + +from pymoo.optimize import minimize +from pymoo.util.ref_dirs import get_reference_directions +from pymoo.operators.crossover.pntx import SinglePointCrossover +from pymoo.operators.mutation.bitflip import BitflipMutation +from pymoo.operators.sampling.rnd import BinaryRandomSampling +from pymoo.indicators.hv import HV +from pymoo.util.nds.non_dominated_sorting import NonDominatedSorting + +import darwin.GlobalVars as GlobalVars + +from darwin.Log import log +from darwin.options import options +from darwin.ExecutionManager import keep_going +from darwin.ModelCode import ModelCode +from darwin.Population import Population +from darwin.Template import Template +from darwin.Model import Model +from darwin.ModelRun import ModelRun + + +warnings.filterwarnings('ignore', category=DeprecationWarning) + +logger = logging.getLogger(__name__) + +def _get_n_params(run: ModelRun) -> int: + model = run.model + + return model.estimated_omega_num + model.estimated_theta_num + model.estimated_sigma_num + + + +class MOGAProblem(Problem): + def __init__(self, n_var, genres, **kwargs): + super().__init__(n_var=n_var, + n_obj=3, + n_constr=3, + xl=np.zeros(n_var, dtype=int), + xu=np.ones(n_var, dtype=int), + requires_kwargs=True + ) + self.genres=genres + + + + def _evaluate(self, x, out, *args, **kwargs): + + f1 = self.genres[:,0] + f2 = self.genres[:,1] + f3 = self.genres[:,2] + + g1 = f1 - 999999 + g2 = 0.1-f2 + g3 = f3 - 999999 + + + + out["F"] = np.column_stack([f1, f2, f3]) + out["G"] = np.column_stack([g1, g2, g3]) + + + + +class _MOGARunner: + def __init__(self, template: Template, pop_size, num_generations): + self.generation = 0 + self.template = template + self.population = None + self.num_generations = num_generations + + + +def run_moga(model_template: Template) -> ModelRun: + # Determine the number of variables based on the model template's gene lengths. + n_var = int(np.sum(model_template.gene_length)) + # Retrieve genetic algorithm options. + ga_options = options.GA + crossover_probability = ga_options['crossover_rate'] + mutation_probability = ga_options['mutation_rate'] + attribute_mutation_probability = ga_options['attribute_mutation_probability'] + pop_size = options.population_size + n_gens = options.num_generations + # Create reference directions for the NSGA-III algorithm. + ref_dirs = get_reference_directions("das-dennis", 3, n_partitions=12) + # Initialize the NSGA-III runner. + runner = _MOGARunner(model_template, pop_size, options.num_generations) + runner.problem = MOGAProblem(n_var=n_var, genres=None) + # Configure the NSGA-III algorithm. + runner.algorithm = NSGA3(ref_dirs=ref_dirs, + pop_size=pop_size, + sampling=BinaryRandomSampling(), + crossover=SinglePointCrossover(prob=crossover_probability), + mutation=BitflipMutation(prob=mutation_probability,prob_var=attribute_mutation_probability), + eliminate_duplicates=True) + # Set up the algorithm with termination criteria and history tracking. + runner.algorithm.setup(runner.problem, termination=('n_gen', n_gens), seed=1, save_history=True, verbose=False) + # Initialize variables for tracking progress. + n_gen = 0 + pres = np.empty((0, 4)) + hvres = np.empty((0, 2)) + ref_point = np.array([100000000, 26, 100000000]) + # Main optimization loop. + while runner.algorithm.has_next(): + + pop = runner.algorithm.ask() + bitstrings = [ind.X.astype(int) for ind in pop] + n_gen += 1 + population = Population.from_codes(model_template,n_gen, bitstrings, + ModelCode.from_full_binary, + max_iteration=n_gens) + population.run() + run=population.runs + f1 = [r.result.ofv for r in population.runs] + f2 = [_get_n_params(r) for r in + population.runs] + f3 = [r.result.post_run_r_penalty for r in population.runs] + + genres = np.array([f1, f2, f3]).T + + + + n_gen = runner.algorithm.n_gen + n_eval = runner.algorithm.evaluator.n_eval + runner.problem = MOGAProblem(n_var=n_var, genres=genres) + runner.algorithm.evaluator.eval(runner.problem, pop) + runner.algorithm.tell(infills=pop) + parents_objectives = runner.algorithm.pop.get("F") + hv = HV(ref_point=ref_point).do(parents_objectives) + hv_values = np.array([hv, n_gen]).T + iteration_col = np.full((parents_objectives.shape[0], 1), n_gen) + parents_objectives2 = np.hstack((parents_objectives,iteration_col)) + pres = np.vstack((pres, parents_objectives2)) + hvres = np.vstack((hvres, hv_values)) + + # Extract the final results. + res = runner.algorithm.result() + all_X = [] + all_F = [] + for entry in res.history: + all_X.append(entry.pop.get("X")) + all_F.append(entry.pop.get("F")) + all_X = np.vstack(all_X) + all_X = all_X.astype(int) + all_F = np.vstack(all_F) + nds = NonDominatedSorting().do(all_F, only_non_dominated_front=True) + pareto_X = all_X[nds] + pareto_X = np.array([[int("".join(map(str, row)))] for row in pareto_X]) + pareto_F = all_F[nds] + allpareto = np.hstack([pareto_X, pareto_F]) + # Save Pareto fronts + pareto_dir = options.output_dir + '\\' + 'PARETO' + os.makedirs(pareto_dir, exist_ok=True) + output_path = os.path.join(pareto_dir, "allpareto.csv") + np.savetxt(output_path, allpareto, delimiter=",", header="Model,OFV,np,rp", comments="") + + # Save each generation's parents + genres_dir = options.output_dir + '\\' + 'GENRES' + os.makedirs(genres_dir, exist_ok=True) + output_path = os.path.join(genres_dir, "output.csv") + np.savetxt(output_path, pres, delimiter=",", header="f1,f2,f3,G", comments="") + # Save hyper volume + output_path = os.path.join(genres_dir, "hv.csv") + np.savetxt(output_path, hvres, delimiter=",", header="hv,G", comments="") + + log.message(f" NSGA3 ideal Pareto-front genome = {res.X.astype(int)},\n" + f" OFV, # of P, and R Penalty = {res.F}") + + for i in range(len(res.X)): + cur_x = [res.X[i].astype(int)] + cur_population = Population.from_codes(model_template, n_gen + 1, cur_x, + ModelCode.from_full_binary) + cur_run = cur_population.runs[0] + cur_run.run_dir = options.output_dir + '\\' + str(i) + cur_run.make_control_file() + cur_run.output_results() \ No newline at end of file