Mercurial > repos > bgruening > sklearn_train_test_eval

comparison search_model_validation.py @ 0:68aaa903052a 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 60f0fbc0eafd7c11bc60fb6c77f2937782efd8a9-dirty
author bgruening
date Fri, 09 Aug 2019 07:09:06 -0400
parents
children cc49634df38f
comparison
equal deleted inserted replaced
-1:000000000000 0:68aaa903052a
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 params_builder = params['search_schemes']['search_params_builder']
217
218 with open(infile_estimator, 'rb') as estimator_handler:
219 estimator = load_model(estimator_handler)
220 estimator_params = estimator.get_params()
221
222 # store read dataframe object
223 loaded_df = {}
224
225 input_type = params['input_options']['selected_input']
226 # tabular input
227 if input_type == 'tabular':
228 header = 'infer' if params['input_options']['header1'] else None
229 column_option = (params['input_options']['column_selector_options_1']
230 ['selected_column_selector_option'])
231 if column_option in ['by_index_number', 'all_but_by_index_number',
232 'by_header_name', 'all_but_by_header_name']:
233 c = params['input_options']['column_selector_options_1']['col1']
234 else:
235 c = None
236
237 df_key = infile1 + repr(header)
238 df = pd.read_csv(infile1, sep='\t', header=header,
239 parse_dates=True)
240 loaded_df[df_key] = df
241
242 X = read_columns(df, c=c, c_option=column_option).astype(float)
243 # sparse input
244 elif input_type == 'sparse':
245 X = mmread(open(infile1, 'r'))
246
247 # fasta_file input
248 elif input_type == 'seq_fasta':
249 pyfaidx = get_module('pyfaidx')
250 sequences = pyfaidx.Fasta(fasta_path)
251 n_seqs = len(sequences.keys())
252 X = np.arange(n_seqs)[:, np.newaxis]
253 for param in estimator_params.keys():
254 if param.endswith('fasta_path'):
255 estimator.set_params(
256 **{param: fasta_path})
257 break
258 else:
259 raise ValueError(
260 "The selected estimator doesn't support "
261 "fasta file input! Please consider using "
262 "KerasGBatchClassifier with "
263 "FastaDNABatchGenerator/FastaProteinBatchGenerator "
264 "or having GenomeOneHotEncoder/ProteinOneHotEncoder "
265 "in pipeline!")
266
267 elif input_type == 'refseq_and_interval':
268 path_params = {
269 'data_batch_generator__ref_genome_path': ref_seq,
270 'data_batch_generator__intervals_path': intervals,
271 'data_batch_generator__target_path': targets
272 }
273 estimator.set_params(**path_params)
274 n_intervals = sum(1 for line in open(intervals))
275 X = np.arange(n_intervals)[:, np.newaxis]
276
277 # Get target y
278 header = 'infer' if params['input_options']['header2'] else None
279 column_option = (params['input_options']['column_selector_options_2']
280 ['selected_column_selector_option2'])
281 if column_option in ['by_index_number', 'all_but_by_index_number',
282 'by_header_name', 'all_but_by_header_name']:
283 c = params['input_options']['column_selector_options_2']['col2']
284 else:
285 c = None
286
287 df_key = infile2 + repr(header)
288 if df_key in loaded_df:
289 infile2 = loaded_df[df_key]
290 else:
291 infile2 = pd.read_csv(infile2, sep='\t',
292 header=header, parse_dates=True)
293 loaded_df[df_key] = infile2
294
295 y = read_columns(
296 infile2,
297 c=c,
298 c_option=column_option,
299 sep='\t',
300 header=header,
301 parse_dates=True)
302 if len(y.shape) == 2 and y.shape[1] == 1:
303 y = y.ravel()
304 if input_type == 'refseq_and_interval':
305 estimator.set_params(
306 data_batch_generator__features=y.ravel().tolist())
307 y = None
308 # end y
309
310 optimizer = params['search_schemes']['selected_search_scheme']
311 optimizer = getattr(model_selection, optimizer)
312
313 # handle gridsearchcv options
314 options = params['search_schemes']['options']
315
316 if groups:
317 header = 'infer' if (options['cv_selector']['groups_selector']
318 ['header_g']) else None
319 column_option = (options['cv_selector']['groups_selector']
320 ['column_selector_options_g']
321 ['selected_column_selector_option_g'])
322 if column_option in ['by_index_number', 'all_but_by_index_number',
323 'by_header_name', 'all_but_by_header_name']:
324 c = (options['cv_selector']['groups_selector']
325 ['column_selector_options_g']['col_g'])
326 else:
327 c = None
328
329 df_key = groups + repr(header)
330 if df_key in loaded_df:
331 groups = loaded_df[df_key]
332
333 groups = read_columns(
334 groups,
335 c=c,
336 c_option=column_option,
337 sep='\t',
338 header=header,
339 parse_dates=True)
340 groups = groups.ravel()
341 options['cv_selector']['groups_selector'] = groups
342
343 splitter, groups = get_cv(options.pop('cv_selector'))
344 options['cv'] = splitter
345 options['n_jobs'] = N_JOBS
346 primary_scoring = options['scoring']['primary_scoring']
347 options['scoring'] = get_scoring(options['scoring'])
348 if options['error_score']:
349 options['error_score'] = 'raise'
350 else:
351 options['error_score'] = np.NaN
352 if options['refit'] and isinstance(options['scoring'], dict):
353 options['refit'] = primary_scoring
354 if 'pre_dispatch' in options and options['pre_dispatch'] == '':
355 options['pre_dispatch'] = None
356
357 # del loaded_df
358 del loaded_df
359
360 # handle memory
361 memory = joblib.Memory(location=CACHE_DIR, verbose=0)
362 # cache iraps_core fits could increase search speed significantly
363 if estimator.__class__.__name__ == 'IRAPSClassifier':
364 estimator.set_params(memory=memory)
365 else:
366 # For iraps buried in pipeline
367 for p, v in estimator_params.items():
368 if p.endswith('memory'):
369 # for case of `__irapsclassifier__memory`
370 if len(p) > 8 and p[:-8].endswith('irapsclassifier'):
371 # cache iraps_core fits could increase search
372 # speed significantly
373 new_params = {p: memory}
374 estimator.set_params(**new_params)
375 # security reason, we don't want memory being
376 # modified unexpectedly
377 elif v:
378 new_params = {p, None}
379 estimator.set_params(**new_params)
380 # For now, 1 CPU is suggested for iprasclassifier
381 elif p.endswith('n_jobs'):
382 new_params = {p: 1}
383 estimator.set_params(**new_params)
384 # for security reason, types of callbacks are limited
385 elif p.endswith('callbacks'):
386 for cb in v:
387 cb_type = cb['callback_selection']['callback_type']
388 if cb_type not in ALLOWED_CALLBACKS:
389 raise ValueError(
390 "Prohibited callback type: %s!" % cb_type)
391
392 param_grid = _eval_search_params(params_builder)
393 searcher = optimizer(estimator, param_grid, **options)
394
395 # do nested split
396 split_mode = params['outer_split'].pop('split_mode')
397 # nested CV, outer cv using cross_validate
398 if split_mode == 'nested_cv':
399 outer_cv, _ = get_cv(params['outer_split']['cv_selector'])
400
401 if options['error_score'] == 'raise':
402 rval = cross_validate(
403 searcher, X, y, scoring=options['scoring'],
404 cv=outer_cv, n_jobs=N_JOBS, verbose=0,
405 error_score=options['error_score'])
406 else:
407 warnings.simplefilter('always', FitFailedWarning)
408 with warnings.catch_warnings(record=True) as w:
409 try:
410 rval = cross_validate(
411 searcher, X, y,
412 scoring=options['scoring'],
413 cv=outer_cv, n_jobs=N_JOBS,
414 verbose=0,
415 error_score=options['error_score'])
416 except ValueError:
417 pass
418 for warning in w:
419 print(repr(warning.message))
420
421 keys = list(rval.keys())
422 for k in keys:
423 if k.startswith('test'):
424 rval['mean_' + k] = np.mean(rval[k])
425 rval['std_' + k] = np.std(rval[k])
426 if k.endswith('time'):
427 rval.pop(k)
428 rval = pd.DataFrame(rval)
429 rval = rval[sorted(rval.columns)]
430 rval.to_csv(path_or_buf=outfile_result, sep='\t',
431 header=True, index=False)
432 else:
433 if split_mode == 'train_test_split':
434 train_test_split = try_get_attr(
435 'galaxy_ml.model_validations', 'train_test_split')
436 # make sure refit is choosen
437 # this could be True for sklearn models, but not the case for
438 # deep learning models
439 if not options['refit'] and \
440 not all(hasattr(estimator, attr)
441 for attr in ('config', 'model_type')):
442 warnings.warn("Refit is change to `True` for nested "
443 "validation!")
444 setattr(searcher, 'refit', True)
445 split_options = params['outer_split']
446
447 # splits
448 if split_options['shuffle'] == 'stratified':
449 split_options['labels'] = y
450 X, X_test, y, y_test = train_test_split(X, y, **split_options)
451 elif split_options['shuffle'] == 'group':
452 if groups is None:
453 raise ValueError("No group based CV option was "
454 "choosen for group shuffle!")
455 split_options['labels'] = groups
456 if y is None:
457 X, X_test, groups, _ =\
458 train_test_split(X, groups, **split_options)
459 else:
460 X, X_test, y, y_test, groups, _ =\
461 train_test_split(X, y, groups, **split_options)
462 else:
463 if split_options['shuffle'] == 'None':
464 split_options['shuffle'] = None
465 X, X_test, y, y_test =\
466 train_test_split(X, y, **split_options)
467 # end train_test_split
468
469 # shared by both train_test_split and non-split
470 if options['error_score'] == 'raise':
471 searcher.fit(X, y, groups=groups)
472 else:
473 warnings.simplefilter('always', FitFailedWarning)
474 with warnings.catch_warnings(record=True) as w:
475 try:
476 searcher.fit(X, y, groups=groups)
477 except ValueError:
478 pass
479 for warning in w:
480 print(repr(warning.message))
481
482 # no outer split
483 if split_mode == 'no':
484 # save results
485 cv_results = pd.DataFrame(searcher.cv_results_)
486 cv_results = cv_results[sorted(cv_results.columns)]
487 cv_results.to_csv(path_or_buf=outfile_result, sep='\t',
488 header=True, index=False)
489
490 # train_test_split, output test result using best_estimator_
491 # or rebuild the trained estimator using weights if applicable.
492 else:
493 scorer_ = searcher.scorer_
494 if isinstance(scorer_, collections.Mapping):
495 is_multimetric = True
496 else:
497 is_multimetric = False
498
499 best_estimator_ = getattr(searcher, 'best_estimator_', None)
500 if not best_estimator_:
501 raise ValueError("GridSearchCV object has no "
502 "`best_estimator_` when `refit`=False!")
503
504 if best_estimator_.__class__.__name__ == 'KerasGBatchClassifier' \
505 and hasattr(estimator.data_batch_generator, 'target_path'):
506 test_score = best_estimator_.evaluate(
507 X_test, scorer=scorer_, is_multimetric=is_multimetric)
508 else:
509 test_score = _score(best_estimator_, X_test,
510 y_test, scorer_,
511 is_multimetric=is_multimetric)
512
513 if not is_multimetric:
514 test_score = {primary_scoring: test_score}
515 for key, value in test_score.items():
516 test_score[key] = [value]
517 result_df = pd.DataFrame(test_score)
518 result_df.to_csv(path_or_buf=outfile_result, sep='\t',
519 header=True, index=False)
520
521 memory.clear(warn=False)
522
523 if outfile_object:
524 best_estimator_ = getattr(searcher, 'best_estimator_', None)
525 if not best_estimator_:
526 warnings.warn("GridSearchCV object has no attribute "
527 "'best_estimator_', because either it's "
528 "nested gridsearch or `refit` is False!")
529 return
530
531 main_est = best_estimator_
532 if isinstance(best_estimator_, pipeline.Pipeline):
533 main_est = best_estimator_.steps[-1][-1]
534
535 if hasattr(main_est, 'model_') \
536 and hasattr(main_est, 'save_weights'):
537 if outfile_weights:
538 main_est.save_weights(outfile_weights)
539 del main_est.model_
540 del main_est.fit_params
541 del main_est.model_class_
542 del main_est.validation_data
543 if getattr(main_est, 'data_generator_', None):
544 del main_est.data_generator_
545 del main_est.data_batch_generator
546
547 with open(outfile_object, 'wb') as output_handler:
548 pickle.dump(best_estimator_, output_handler,
549 pickle.HIGHEST_PROTOCOL)
550
551
552 if __name__ == '__main__':
553 aparser = argparse.ArgumentParser()
554 aparser.add_argument("-i", "--inputs", dest="inputs", required=True)
555 aparser.add_argument("-e", "--estimator", dest="infile_estimator")
556 aparser.add_argument("-X", "--infile1", dest="infile1")
557 aparser.add_argument("-y", "--infile2", dest="infile2")
558 aparser.add_argument("-O", "--outfile_result", dest="outfile_result")
559 aparser.add_argument("-o", "--outfile_object", dest="outfile_object")
560 aparser.add_argument("-w", "--outfile_weights", dest="outfile_weights")
561 aparser.add_argument("-g", "--groups", dest="groups")
562 aparser.add_argument("-r", "--ref_seq", dest="ref_seq")
563 aparser.add_argument("-b", "--intervals", dest="intervals")
564 aparser.add_argument("-t", "--targets", dest="targets")
565 aparser.add_argument("-f", "--fasta_path", dest="fasta_path")
566 args = aparser.parse_args()
567
568 main(args.inputs, args.infile_estimator, args.infile1, args.infile2,
569 args.outfile_result, outfile_object=args.outfile_object,
570 outfile_weights=args.outfile_weights, groups=args.groups,
571 ref_seq=args.ref_seq, intervals=args.intervals,
572 targets=args.targets, fasta_path=args.fasta_path)