Mercurial > repos > bgruening > sklearn_fitted_model_eval

comparison search_model_validation.py @ 0:eaddff553324 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/bgruening/galaxytools/tree/master/tools/sklearn commit eb703290e2589561ea215c84aa9f71bcfe1712c6"
author bgruening
date Fri, 01 Nov 2019 17:15:22 -0400
parents
children cf54bae8ad42
comparison
equal deleted inserted replaced
-1:000000000000 0:eaddff553324
1 import argparse
2 import collections
3 import imblearn
4 import joblib
5 import json
6 import numpy as np
7 import pandas as pd
8 import pickle
9 import skrebate
10 import sklearn
11 import sys
12 import xgboost
13 import warnings
14 from imblearn import under_sampling, over_sampling, combine
15 from scipy.io import mmread
16 from mlxtend import classifier, regressor
17 from sklearn.base import clone
18 from sklearn import (cluster, compose, decomposition, ensemble,
19 feature_extraction, feature_selection,
20 gaussian_process, kernel_approximation, metrics,
21 model_selection, naive_bayes, neighbors,
22 pipeline, preprocessing, svm, linear_model,
23 tree, discriminant_analysis)
24 from sklearn.exceptions import FitFailedWarning
25 from sklearn.model_selection._validation import _score, cross_validate
26 from sklearn.model_selection import _search, _validation
27
28 from galaxy_ml.utils import (SafeEval, get_cv, get_scoring, load_model,
29 read_columns, try_get_attr, get_module)
30
31
32 _fit_and_score = try_get_attr('galaxy_ml.model_validations', '_fit_and_score')
33 setattr(_search, '_fit_and_score', _fit_and_score)
34 setattr(_validation, '_fit_and_score', _fit_and_score)
35
36 N_JOBS = int(__import__('os').environ.get('GALAXY_SLOTS', 1))
37 CACHE_DIR = './cached'
38 NON_SEARCHABLE = ('n_jobs', 'pre_dispatch', 'memory', '_path',
39 'nthread', 'callbacks')
40 ALLOWED_CALLBACKS = ('EarlyStopping', 'TerminateOnNaN', 'ReduceLROnPlateau',
41 'CSVLogger', 'None')
42
43
44 def _eval_search_params(params_builder):
45 search_params = {}
46
47 for p in params_builder['param_set']:
48 search_list = p['sp_list'].strip()
49 if search_list == '':
50 continue
51
52 param_name = p['sp_name']
53 if param_name.lower().endswith(NON_SEARCHABLE):
54 print("Warning: `%s` is not eligible for search and was "
55 "omitted!" % param_name)
56 continue
57
58 if not search_list.startswith(':'):
59 safe_eval = SafeEval(load_scipy=True, load_numpy=True)
60 ev = safe_eval(search_list)
61 search_params[param_name] = ev
62 else:
63 # Have `:` before search list, asks for estimator evaluatio
64 safe_eval_es = SafeEval(load_estimators=True)
65 search_list = search_list[1:].strip()
66 # TODO maybe add regular express check
67 ev = safe_eval_es(search_list)
68 preprocessings = (
69 preprocessing.StandardScaler(), preprocessing.Binarizer(),
70 preprocessing.MaxAbsScaler(),
71 preprocessing.Normalizer(), preprocessing.MinMaxScaler(),
72 preprocessing.PolynomialFeatures(),
73 preprocessing.RobustScaler(), feature_selection.SelectKBest(),
74 feature_selection.GenericUnivariateSelect(),
75 feature_selection.SelectPercentile(),
76 feature_selection.SelectFpr(), feature_selection.SelectFdr(),
77 feature_selection.SelectFwe(),
78 feature_selection.VarianceThreshold(),
79 decomposition.FactorAnalysis(random_state=0),
80 decomposition.FastICA(random_state=0),
81 decomposition.IncrementalPCA(),
82 decomposition.KernelPCA(random_state=0, n_jobs=N_JOBS),
83 decomposition.LatentDirichletAllocation(
84 random_state=0, n_jobs=N_JOBS),
85 decomposition.MiniBatchDictionaryLearning(
86 random_state=0, n_jobs=N_JOBS),
87 decomposition.MiniBatchSparsePCA(
88 random_state=0, n_jobs=N_JOBS),
89 decomposition.NMF(random_state=0),
90 decomposition.PCA(random_state=0),
91 decomposition.SparsePCA(random_state=0, n_jobs=N_JOBS),
92 decomposition.TruncatedSVD(random_state=0),
93 kernel_approximation.Nystroem(random_state=0),
94 kernel_approximation.RBFSampler(random_state=0),
95 kernel_approximation.AdditiveChi2Sampler(),
96 kernel_approximation.SkewedChi2Sampler(random_state=0),
97 cluster.FeatureAgglomeration(),
98 skrebate.ReliefF(n_jobs=N_JOBS),
99 skrebate.SURF(n_jobs=N_JOBS),
100 skrebate.SURFstar(n_jobs=N_JOBS),
101 skrebate.MultiSURF(n_jobs=N_JOBS),
102 skrebate.MultiSURFstar(n_jobs=N_JOBS),
103 imblearn.under_sampling.ClusterCentroids(
104 random_state=0, n_jobs=N_JOBS),
105 imblearn.under_sampling.CondensedNearestNeighbour(
106 random_state=0, n_jobs=N_JOBS),
107 imblearn.under_sampling.EditedNearestNeighbours(
108 random_state=0, n_jobs=N_JOBS),
109 imblearn.under_sampling.RepeatedEditedNearestNeighbours(
110 random_state=0, n_jobs=N_JOBS),
111 imblearn.under_sampling.AllKNN(random_state=0, n_jobs=N_JOBS),
112 imblearn.under_sampling.InstanceHardnessThreshold(
113 random_state=0, n_jobs=N_JOBS),
114 imblearn.under_sampling.NearMiss(
115 random_state=0, n_jobs=N_JOBS),
116 imblearn.under_sampling.NeighbourhoodCleaningRule(
117 random_state=0, n_jobs=N_JOBS),
118 imblearn.under_sampling.OneSidedSelection(
119 random_state=0, n_jobs=N_JOBS),
120 imblearn.under_sampling.RandomUnderSampler(
121 random_state=0),
122 imblearn.under_sampling.TomekLinks(
123 random_state=0, n_jobs=N_JOBS),
124 imblearn.over_sampling.ADASYN(random_state=0, n_jobs=N_JOBS),
125 imblearn.over_sampling.RandomOverSampler(random_state=0),
126 imblearn.over_sampling.SMOTE(random_state=0, n_jobs=N_JOBS),
127 imblearn.over_sampling.SVMSMOTE(random_state=0, n_jobs=N_JOBS),
128 imblearn.over_sampling.BorderlineSMOTE(
129 random_state=0, n_jobs=N_JOBS),
130 imblearn.over_sampling.SMOTENC(
131 categorical_features=[], random_state=0, n_jobs=N_JOBS),
132 imblearn.combine.SMOTEENN(random_state=0),
133 imblearn.combine.SMOTETomek(random_state=0))
134 newlist = []
135 for obj in ev:
136 if obj is None:
137 newlist.append(None)
138 elif obj == 'all_0':
139 newlist.extend(preprocessings[0:35])
140 elif obj == 'sk_prep_all': # no KernalCenter()
141 newlist.extend(preprocessings[0:7])
142 elif obj == 'fs_all':
143 newlist.extend(preprocessings[7:14])
144 elif obj == 'decomp_all':
145 newlist.extend(preprocessings[14:25])
146 elif obj == 'k_appr_all':
147 newlist.extend(preprocessings[25:29])
148 elif obj == 'reb_all':
149 newlist.extend(preprocessings[30:35])
150 elif obj == 'imb_all':
151 newlist.extend(preprocessings[35:54])
152 elif type(obj) is int and -1 < obj < len(preprocessings):
153 newlist.append(preprocessings[obj])
154 elif hasattr(obj, 'get_params'): # user uploaded object
155 if 'n_jobs' in obj.get_params():
156 newlist.append(obj.set_params(n_jobs=N_JOBS))
157 else:
158 newlist.append(obj)
159 else:
160 sys.exit("Unsupported estimator type: %r" % (obj))
161
162 search_params[param_name] = newlist
163
164 return search_params
165
166
167 def main(inputs, infile_estimator, infile1, infile2,
168 outfile_result, outfile_object=None,
169 outfile_weights=None, groups=None,
170 ref_seq=None, intervals=None, targets=None,
171 fasta_path=None):
172 """
173 Parameter
174 ---------
175 inputs : str
176 File path to galaxy tool parameter
177
178 infile_estimator : str
179 File path to estimator
180
181 infile1 : str
182 File path to dataset containing features
183
184 infile2 : str
185 File path to dataset containing target values
186
187 outfile_result : str
188 File path to save the results, either cv_results or test result
189
190 outfile_object : str, optional
191 File path to save searchCV object
192
193 outfile_weights : str, optional
194 File path to save model weights
195
196 groups : str
197 File path to dataset containing groups labels
198
199 ref_seq : str
200 File path to dataset containing genome sequence file
201
202 intervals : str
203 File path to dataset containing interval file
204
205 targets : str
206 File path to dataset compressed target bed file
207
208 fasta_path : str
209 File path to dataset containing fasta file
210 """
211 warnings.simplefilter('ignore')
212
213 with open(inputs, 'r') as param_handler:
214 params = json.load(param_handler)
215
216 # conflict param checker
217 if params['outer_split']['split_mode'] == 'nested_cv' \
218 and params['save'] != 'nope':
219 raise ValueError("Save best estimator is not possible for nested CV!")
220
221 if not (params['search_schemes']['options']['refit']) \
222 and params['save'] != 'nope':
223 raise ValueError("Save best estimator is not possible when refit "
224 "is False!")
225
226 params_builder = params['search_schemes']['search_params_builder']
227
228 with open(infile_estimator, 'rb') as estimator_handler:
229 estimator = load_model(estimator_handler)
230 estimator_params = estimator.get_params()
231
232 # store read dataframe object
233 loaded_df = {}
234
235 input_type = params['input_options']['selected_input']
236 # tabular input
237 if input_type == 'tabular':
238 header = 'infer' if params['input_options']['header1'] else None
239 column_option = (params['input_options']['column_selector_options_1']
240 ['selected_column_selector_option'])
241 if column_option in ['by_index_number', 'all_but_by_index_number',
242 'by_header_name', 'all_but_by_header_name']:
243 c = params['input_options']['column_selector_options_1']['col1']
244 else:
245 c = None
246
247 df_key = infile1 + repr(header)
248 df = pd.read_csv(infile1, sep='\t', header=header,
249 parse_dates=True)
250 loaded_df[df_key] = df
251
252 X = read_columns(df, c=c, c_option=column_option).astype(float)
253 # sparse input
254 elif input_type == 'sparse':
255 X = mmread(open(infile1, 'r'))
256
257 # fasta_file input
258 elif input_type == 'seq_fasta':
259 pyfaidx = get_module('pyfaidx')
260 sequences = pyfaidx.Fasta(fasta_path)
261 n_seqs = len(sequences.keys())
262 X = np.arange(n_seqs)[:, np.newaxis]
263 for param in estimator_params.keys():
264 if param.endswith('fasta_path'):
265 estimator.set_params(
266 **{param: fasta_path})
267 break
268 else:
269 raise ValueError(
270 "The selected estimator doesn't support "
271 "fasta file input! Please consider using "
272 "KerasGBatchClassifier with "
273 "FastaDNABatchGenerator/FastaProteinBatchGenerator "
274 "or having GenomeOneHotEncoder/ProteinOneHotEncoder "
275 "in pipeline!")
276
277 elif input_type == 'refseq_and_interval':
278 path_params = {
279 'data_batch_generator__ref_genome_path': ref_seq,
280 'data_batch_generator__intervals_path': intervals,
281 'data_batch_generator__target_path': targets
282 }
283 estimator.set_params(**path_params)
284 n_intervals = sum(1 for line in open(intervals))
285 X = np.arange(n_intervals)[:, np.newaxis]
286
287 # Get target y
288 header = 'infer' if params['input_options']['header2'] else None
289 column_option = (params['input_options']['column_selector_options_2']
290 ['selected_column_selector_option2'])
291 if column_option in ['by_index_number', 'all_but_by_index_number',
292 'by_header_name', 'all_but_by_header_name']:
293 c = params['input_options']['column_selector_options_2']['col2']
294 else:
295 c = None
296
297 df_key = infile2 + repr(header)
298 if df_key in loaded_df:
299 infile2 = loaded_df[df_key]
300 else:
301 infile2 = pd.read_csv(infile2, sep='\t',
302 header=header, parse_dates=True)
303 loaded_df[df_key] = infile2
304
305 y = read_columns(
306 infile2,
307 c=c,
308 c_option=column_option,
309 sep='\t',
310 header=header,
311 parse_dates=True)
312 if len(y.shape) == 2 and y.shape[1] == 1:
313 y = y.ravel()
314 if input_type == 'refseq_and_interval':
315 estimator.set_params(
316 data_batch_generator__features=y.ravel().tolist())
317 y = None
318 # end y
319
320 optimizer = params['search_schemes']['selected_search_scheme']
321 optimizer = getattr(model_selection, optimizer)
322
323 # handle gridsearchcv options
324 options = params['search_schemes']['options']
325
326 if groups:
327 header = 'infer' if (options['cv_selector']['groups_selector']
328 ['header_g']) else None
329 column_option = (options['cv_selector']['groups_selector']
330 ['column_selector_options_g']
331 ['selected_column_selector_option_g'])
332 if column_option in ['by_index_number', 'all_but_by_index_number',
333 'by_header_name', 'all_but_by_header_name']:
334 c = (options['cv_selector']['groups_selector']
335 ['column_selector_options_g']['col_g'])
336 else:
337 c = None
338
339 df_key = groups + repr(header)
340 if df_key in loaded_df:
341 groups = loaded_df[df_key]
342
343 groups = read_columns(
344 groups,
345 c=c,
346 c_option=column_option,
347 sep='\t',
348 header=header,
349 parse_dates=True)
350 groups = groups.ravel()
351 options['cv_selector']['groups_selector'] = groups
352
353 splitter, groups = get_cv(options.pop('cv_selector'))
354 options['cv'] = splitter
355 options['n_jobs'] = N_JOBS
356 primary_scoring = options['scoring']['primary_scoring']
357 options['scoring'] = get_scoring(options['scoring'])
358 if options['error_score']:
359 options['error_score'] = 'raise'
360 else:
361 options['error_score'] = np.NaN
362 if options['refit'] and isinstance(options['scoring'], dict):
363 options['refit'] = primary_scoring
364 if 'pre_dispatch' in options and options['pre_dispatch'] == '':
365 options['pre_dispatch'] = None
366
367 # del loaded_df
368 del loaded_df
369
370 # handle memory
371 memory = joblib.Memory(location=CACHE_DIR, verbose=0)
372 # cache iraps_core fits could increase search speed significantly
373 if estimator.__class__.__name__ == 'IRAPSClassifier':
374 estimator.set_params(memory=memory)
375 else:
376 # For iraps buried in pipeline
377 for p, v in estimator_params.items():
378 if p.endswith('memory'):
379 # for case of `__irapsclassifier__memory`
380 if len(p) > 8 and p[:-8].endswith('irapsclassifier'):
381 # cache iraps_core fits could increase search
382 # speed significantly
383 new_params = {p: memory}
384 estimator.set_params(**new_params)
385 # security reason, we don't want memory being
386 # modified unexpectedly
387 elif v:
388 new_params = {p, None}
389 estimator.set_params(**new_params)
390 # For now, 1 CPU is suggested for iprasclassifier
391 elif p.endswith('n_jobs'):
392 new_params = {p: 1}
393 estimator.set_params(**new_params)
394 # for security reason, types of callbacks are limited
395 elif p.endswith('callbacks'):
396 for cb in v:
397 cb_type = cb['callback_selection']['callback_type']
398 if cb_type not in ALLOWED_CALLBACKS:
399 raise ValueError(
400 "Prohibited callback type: %s!" % cb_type)
401
402 param_grid = _eval_search_params(params_builder)
403 searcher = optimizer(estimator, param_grid, **options)
404
405 # do nested split
406 split_mode = params['outer_split'].pop('split_mode')
407 # nested CV, outer cv using cross_validate
408 if split_mode == 'nested_cv':
409 outer_cv, _ = get_cv(params['outer_split']['cv_selector'])
410
411 if options['error_score'] == 'raise':
412 rval = cross_validate(
413 searcher, X, y, scoring=options['scoring'],
414 cv=outer_cv, n_jobs=N_JOBS, verbose=0,
415 error_score=options['error_score'])
416 else:
417 warnings.simplefilter('always', FitFailedWarning)
418 with warnings.catch_warnings(record=True) as w:
419 try:
420 rval = cross_validate(
421 searcher, X, y,
422 scoring=options['scoring'],
423 cv=outer_cv, n_jobs=N_JOBS,
424 verbose=0,
425 error_score=options['error_score'])
426 except ValueError:
427 pass
428 for warning in w:
429 print(repr(warning.message))
430
431 keys = list(rval.keys())
432 for k in keys:
433 if k.startswith('test'):
434 rval['mean_' + k] = np.mean(rval[k])
435 rval['std_' + k] = np.std(rval[k])
436 if k.endswith('time'):
437 rval.pop(k)
438 rval = pd.DataFrame(rval)
439 rval = rval[sorted(rval.columns)]
440 rval.to_csv(path_or_buf=outfile_result, sep='\t',
441 header=True, index=False)
442 else:
443 if split_mode == 'train_test_split':
444 train_test_split = try_get_attr(
445 'galaxy_ml.model_validations', 'train_test_split')
446 # make sure refit is choosen
447 # this could be True for sklearn models, but not the case for
448 # deep learning models
449 if not options['refit'] and \
450 not all(hasattr(estimator, attr)
451 for attr in ('config', 'model_type')):
452 warnings.warn("Refit is change to `True` for nested "
453 "validation!")
454 setattr(searcher, 'refit', True)
455 split_options = params['outer_split']
456
457 # splits
458 if split_options['shuffle'] == 'stratified':
459 split_options['labels'] = y
460 X, X_test, y, y_test = train_test_split(X, y, **split_options)
461 elif split_options['shuffle'] == 'group':
462 if groups is None:
463 raise ValueError("No group based CV option was "
464 "choosen for group shuffle!")
465 split_options['labels'] = groups
466 if y is None:
467 X, X_test, groups, _ =\
468 train_test_split(X, groups, **split_options)
469 else:
470 X, X_test, y, y_test, groups, _ =\
471 train_test_split(X, y, groups, **split_options)
472 else:
473 if split_options['shuffle'] == 'None':
474 split_options['shuffle'] = None
475 X, X_test, y, y_test =\
476 train_test_split(X, y, **split_options)
477 # end train_test_split
478
479 # shared by both train_test_split and non-split
480 if options['error_score'] == 'raise':
481 searcher.fit(X, y, groups=groups)
482 else:
483 warnings.simplefilter('always', FitFailedWarning)
484 with warnings.catch_warnings(record=True) as w:
485 try:
486 searcher.fit(X, y, groups=groups)
487 except ValueError:
488 pass
489 for warning in w:
490 print(repr(warning.message))
491
492 # no outer split
493 if split_mode == 'no':
494 # save results
495 cv_results = pd.DataFrame(searcher.cv_results_)
496 cv_results = cv_results[sorted(cv_results.columns)]
497 cv_results.to_csv(path_or_buf=outfile_result, sep='\t',
498 header=True, index=False)
499
500 # train_test_split, output test result using best_estimator_
501 # or rebuild the trained estimator using weights if applicable.
502 else:
503 scorer_ = searcher.scorer_
504 if isinstance(scorer_, collections.Mapping):
505 is_multimetric = True
506 else:
507 is_multimetric = False
508
509 best_estimator_ = getattr(searcher, 'best_estimator_', None)
510 if not best_estimator_:
511 raise ValueError("GridSearchCV object has no "
512 "`best_estimator_` when `refit`=False!")
513
514 if best_estimator_.__class__.__name__ == 'KerasGBatchClassifier' \
515 and hasattr(estimator.data_batch_generator, 'target_path'):
516 test_score = best_estimator_.evaluate(
517 X_test, scorer=scorer_, is_multimetric=is_multimetric)
518 else:
519 test_score = _score(best_estimator_, X_test,
520 y_test, scorer_,
521 is_multimetric=is_multimetric)
522
523 if not is_multimetric:
524 test_score = {primary_scoring: test_score}
525 for key, value in test_score.items():
526 test_score[key] = [value]
527 result_df = pd.DataFrame(test_score)
528 result_df.to_csv(path_or_buf=outfile_result, sep='\t',
529 header=True, index=False)
530
531 memory.clear(warn=False)
532
533 if outfile_object:
534 best_estimator_ = getattr(searcher, 'best_estimator_', None)
535 if not best_estimator_:
536 warnings.warn("GridSearchCV object has no attribute "
537 "'best_estimator_', because either it's "
538 "nested gridsearch or `refit` is False!")
539 return
540
541 main_est = best_estimator_
542 if isinstance(best_estimator_, pipeline.Pipeline):
543 main_est = best_estimator_.steps[-1][-1]
544
545 if hasattr(main_est, 'model_') \
546 and hasattr(main_est, 'save_weights'):
547 if outfile_weights:
548 main_est.save_weights(outfile_weights)
549 del main_est.model_
550 del main_est.fit_params
551 del main_est.model_class_
552 del main_est.validation_data
553 if getattr(main_est, 'data_generator_', None):
554 del main_est.data_generator_
555
556 with open(outfile_object, 'wb') as output_handler:
557 pickle.dump(best_estimator_, output_handler,
558 pickle.HIGHEST_PROTOCOL)
559
560
561 if __name__ == '__main__':
562 aparser = argparse.ArgumentParser()
563 aparser.add_argument("-i", "--inputs", dest="inputs", required=True)
564 aparser.add_argument("-e", "--estimator", dest="infile_estimator")
565 aparser.add_argument("-X", "--infile1", dest="infile1")
566 aparser.add_argument("-y", "--infile2", dest="infile2")
567 aparser.add_argument("-O", "--outfile_result", dest="outfile_result")
568 aparser.add_argument("-o", "--outfile_object", dest="outfile_object")
569 aparser.add_argument("-w", "--outfile_weights", dest="outfile_weights")
570 aparser.add_argument("-g", "--groups", dest="groups")
571 aparser.add_argument("-r", "--ref_seq", dest="ref_seq")
572 aparser.add_argument("-b", "--intervals", dest="intervals")
573 aparser.add_argument("-t", "--targets", dest="targets")
574 aparser.add_argument("-f", "--fasta_path", dest="fasta_path")
575 args = aparser.parse_args()
576
577 main(args.inputs, args.infile_estimator, args.infile1, args.infile2,
578 args.outfile_result, outfile_object=args.outfile_object,
579 outfile_weights=args.outfile_weights, groups=args.groups,
580 ref_seq=args.ref_seq, intervals=args.intervals,
581 targets=args.targets, fasta_path=args.fasta_path)