milk-0.3.10/0000755000175000017500000000000011562337213013330 5ustar luispedroluispedromilk-0.3.10/INSTALL.rst0000644000175000017500000000042211545146510015165 0ustar luispedroluispedro============= Building milk ============= The following should work:: python setup.py install A C++ compiler is required. On Windows, you might need to specify the compiler. For example:: python setup.py install --compiler=mingw32 If you have mingw installed. milk-0.3.10/MANIFEST.in0000644000175000017500000000005011545146510015060 0ustar luispedroluispedroinclude README.rst include INSTALL.rst milk-0.3.10/setup.py0000644000175000017500000000456411557532105015054 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2010, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. from __future__ import division try: import setuptools except: print ''' setuptools not found. On linux, the package is often called python-setuptools''' from sys import exit exit(1) from numpy.distutils.core import setup, Extension execfile('milk/milk_version.py') long_description = file('README.rst').read() kmeans_ext = Extension('milk.unsupervised._kmeans', sources = ['milk/unsupervised/_kmeans.cpp']) svm_ext = Extension('milk.supervised._svm', sources = ['milk/supervised/_svm.cpp']) som_ext = Extension('milk.unsupervised._som', sources = ['milk/unsupervised/_som.cpp']) tree_ext = Extension('milk.supervised._tree', sources = ['milk/supervised/_tree.cpp']) ext_modules = [kmeans_ext, svm_ext, som_ext, tree_ext] packages = filter(lambda p: p.startswith('milk'), setuptools.find_packages()) setup(name = 'milk', version = __version__, description = 'Machine Learning Toolkit', long_description = long_description, author = u'Luis Pedro Coelho', author_email = 'lpc@cmu.edu', url = 'http://luispedro.org/software/milk', license = 'MIT', packages = packages, ext_modules = ext_modules, test_suite = 'nose.collector', ) milk-0.3.10/PKG-INFO0000644000175000017500000001341111562337213014425 0ustar luispedroluispedroMetadata-Version: 1.0 Name: milk Version: 0.3.10 Summary: Machine Learning Toolkit Home-page: http://luispedro.org/software/milk Author: Luis Pedro Coelho Author-email: lpc@cmu.edu License: MIT Description: ============================== MILK: MACHINE LEARNING TOOLKIT ============================== Machine Learning in Python -------------------------- Milk is a machine learning toolkit in Python. Its focus is on supervised classification with several classifiers available: SVMs (based on libsvm), k-NN, random forests, decision trees. It also performs feature selection. These classifiers can be combined in many ways to form different classification systems. For unsupervised learning, milk supports k-means clustering and affinity propagation. Milk is flexible about its inputs. It optimised for numpy arrays, but can often handle anything (for example, for SVMs, you can use any dataype and any kernel and it does the right thing). There is a strong emphasis on speed and low memory usage. Therefore, most of the performance sensitive code is in C++. This is behind Python-based interfaces for convenience. To learn more, check the docs at `http://packages.python.org/milk/ `_ or the code demos included with the source at ``milk/demos/``. Examples -------- Here is how to test how well you can classify some ``features,labels`` data, measured by cross-validation:: import numpy as np import milk features = np.random.rand(100,10) # 2d array of features: 100 examples of 10 features each labels = np.zeros(100) features[50:] += .5 labels[50:] = 1 confusion_matrix, names = milk.nfoldcrossvalidation(features, labels) print 'Accuracy:', confusion_matrix.trace()/float(confusion_matrix.sum()) If want to use a classifier, you instanciate a *learner object* and call its ``train()`` method:: import numpy as np import milk features = np.random.rand(100,10) labels = np.zeros(100) features[50:] += .5 labels[50:] = 1 learner = milk.defaultclassifier() model = learner.train(features, labels) # Now you can use the model on new examples: example = np.random.rand(10) print model.apply(example) example2 = np.random.rand(10) example2 += .5 print model.apply(example2) There are several classification methods in the package, but they all use the same interface: ``train()`` returns a *model* object, which has an ``apply()`` method to execute on new instances. Details ------- License: MIT Author: Luis Pedro Coelho (with code from LibSVM and scikits.learn) API Documentation: `http://packages.python.org/milk/ `_ Mailing List: `http://groups.google.com/group/milk-users `__ Features -------- - SVMs. Using the libsvm solver with a pythonesque wrapper around it. - K-means using as little memory as possible. It can cluster millions of instances efficiently. - Random forests - Self organising maps - Stepwise Discriminant Analysis for feature selection. - Non-negative matrix factorisation - Affinity propagation Recent History -------------- The ChangeLog file contains a more complete history. New in 0.3.10 ~~~~~~~~~~~~~ - Add ext.jugparallel for integration with `jug `_ - parallel nfold crossvalidation using jug - parallel multiple kmeans runs using jug - cluster_agreement for non-ndarrays - Add histogram & normali(z|s)e options to ``milk.kmeans.assign_centroid`` - Fix bug in sda when features were constant for a class - Add select_best_kmeans - Added defaultlearner as a better name than defaultclassifier - Add ``measures.curves.precision_recall`` - Add ``unsupervised.parzen.parzen`` New in 0.3.9 ~~~~~~~~~~~~ - Add ``folds`` argument to ``nfoldcrossvalidation`` - Add ``assign_centroid`` function in milk.unsupervised.nfoldcrossvalidation - Improve speed of k-nearest neighbour (10x on scikits-learn benchmark) - Improve kmeans on newer numpy (works for larger datasets too) - Faster kmeans by coding centroid recalculation in C++ - Fix gridminize for low count labels - Fix bug with non-integer labels for tree learning New in 0.3.8 ~~~~~~~~~~~~ - Fix compilation on Windows New in 0.3.7 ~~~~~~~~~~~~ - Logistic regression - Source demos included (in source and documentation) - Add cluster agreement metrics - Fix nfoldcrossvalidation bug when using origins New in 0.3.6 ~~~~~~~~~~~~ - Unsupervised (1-class) kernel density modeling - Fix for when SDA returns empty - weights option to some learners - stump learner - Adaboost (result of above changes) New in 0.3.5 ~~~~~~~~~~~~ - Fixes for 64 bit machines - Functions in measures.py all have same interface now. Platform: UNKNOWN milk-0.3.10/README.rst0000644000175000017500000001067711562335661015037 0ustar luispedroluispedro============================== MILK: MACHINE LEARNING TOOLKIT ============================== Machine Learning in Python -------------------------- Milk is a machine learning toolkit in Python. Its focus is on supervised classification with several classifiers available: SVMs (based on libsvm), k-NN, random forests, decision trees. It also performs feature selection. These classifiers can be combined in many ways to form different classification systems. For unsupervised learning, milk supports k-means clustering and affinity propagation. Milk is flexible about its inputs. It optimised for numpy arrays, but can often handle anything (for example, for SVMs, you can use any dataype and any kernel and it does the right thing). There is a strong emphasis on speed and low memory usage. Therefore, most of the performance sensitive code is in C++. This is behind Python-based interfaces for convenience. To learn more, check the docs at `http://packages.python.org/milk/ `_ or the code demos included with the source at ``milk/demos/``. Examples -------- Here is how to test how well you can classify some ``features,labels`` data, measured by cross-validation:: import numpy as np import milk features = np.random.rand(100,10) # 2d array of features: 100 examples of 10 features each labels = np.zeros(100) features[50:] += .5 labels[50:] = 1 confusion_matrix, names = milk.nfoldcrossvalidation(features, labels) print 'Accuracy:', confusion_matrix.trace()/float(confusion_matrix.sum()) If want to use a classifier, you instanciate a *learner object* and call its ``train()`` method:: import numpy as np import milk features = np.random.rand(100,10) labels = np.zeros(100) features[50:] += .5 labels[50:] = 1 learner = milk.defaultclassifier() model = learner.train(features, labels) # Now you can use the model on new examples: example = np.random.rand(10) print model.apply(example) example2 = np.random.rand(10) example2 += .5 print model.apply(example2) There are several classification methods in the package, but they all use the same interface: ``train()`` returns a *model* object, which has an ``apply()`` method to execute on new instances. Details ------- License: MIT Author: Luis Pedro Coelho (with code from LibSVM and scikits.learn) API Documentation: `http://packages.python.org/milk/ `_ Mailing List: `http://groups.google.com/group/milk-users `__ Features -------- - SVMs. Using the libsvm solver with a pythonesque wrapper around it. - K-means using as little memory as possible. It can cluster millions of instances efficiently. - Random forests - Self organising maps - Stepwise Discriminant Analysis for feature selection. - Non-negative matrix factorisation - Affinity propagation Recent History -------------- The ChangeLog file contains a more complete history. New in 0.3.10 ~~~~~~~~~~~~~ - Add ext.jugparallel for integration with `jug `_ - parallel nfold crossvalidation using jug - parallel multiple kmeans runs using jug - cluster_agreement for non-ndarrays - Add histogram & normali(z|s)e options to ``milk.kmeans.assign_centroid`` - Fix bug in sda when features were constant for a class - Add select_best_kmeans - Added defaultlearner as a better name than defaultclassifier - Add ``measures.curves.precision_recall`` - Add ``unsupervised.parzen.parzen`` New in 0.3.9 ~~~~~~~~~~~~ - Add ``folds`` argument to ``nfoldcrossvalidation`` - Add ``assign_centroid`` function in milk.unsupervised.nfoldcrossvalidation - Improve speed of k-nearest neighbour (10x on scikits-learn benchmark) - Improve kmeans on newer numpy (works for larger datasets too) - Faster kmeans by coding centroid recalculation in C++ - Fix gridminize for low count labels - Fix bug with non-integer labels for tree learning New in 0.3.8 ~~~~~~~~~~~~ - Fix compilation on Windows New in 0.3.7 ~~~~~~~~~~~~ - Logistic regression - Source demos included (in source and documentation) - Add cluster agreement metrics - Fix nfoldcrossvalidation bug when using origins New in 0.3.6 ~~~~~~~~~~~~ - Unsupervised (1-class) kernel density modeling - Fix for when SDA returns empty - weights option to some learners - stump learner - Adaboost (result of above changes) New in 0.3.5 ~~~~~~~~~~~~ - Fixes for 64 bit machines - Functions in measures.py all have same interface now. milk-0.3.10/milk.egg-info/0000755000175000017500000000000011562337213015756 5ustar luispedroluispedromilk-0.3.10/milk.egg-info/dependency_links.txt0000644000175000017500000000000111562337213022024 0ustar luispedroluispedro milk-0.3.10/milk.egg-info/SOURCES.txt0000644000175000017500000000560211562337213017645 0ustar luispedroluispedroINSTALL.rst MANIFEST.in README.rst setup.py milk/__init__.py milk/milk_version.py milk/nfoldcrossvalidation.py milk.egg-info/PKG-INFO milk.egg-info/SOURCES.txt milk.egg-info/dependency_links.txt milk.egg-info/top_level.txt milk/active/__init__.py milk/active/eimpact.py milk/active/uncertainty.py milk/demos/__init__.py milk/demos/adaboost.py milk/demos/rf_wine_2d.py milk/demos/svm-decision-boundary.py milk/ext/__init__.py milk/ext/jugparallel.py milk/measures/__init__.py milk/measures/cluster_agreement.py milk/measures/curves.py milk/measures/measures.py milk/measures/nfoldcrossvalidation.py milk/supervised/__init__.py milk/supervised/_svm.cpp milk/supervised/_tree.cpp milk/supervised/adaboost.py milk/supervised/classifier.py milk/supervised/defaultclassifier.py milk/supervised/defaultlearner.py milk/supervised/featureselection.py milk/supervised/gridsearch.py milk/supervised/grouped.py milk/supervised/knn.py milk/supervised/logistic.py milk/supervised/multi.py milk/supervised/multi_view.py milk/supervised/normalise.py milk/supervised/parzen.py milk/supervised/precluster_learner.py milk/supervised/randomforest.py milk/supervised/set2binary_array.py milk/supervised/svm.py milk/supervised/tree.py milk/tests/__init__.py milk/tests/fast_classifier.py milk/tests/test_adaboost.py milk/tests/test_affinity.py milk/tests/test_basic.py milk/tests/test_curves.py milk/tests/test_defaultclassifier.py milk/tests/test_ext_jugparallel.py milk/tests/test_featureselection.py milk/tests/test_fisher.py milk/tests/test_gaussianmixture.py milk/tests/test_gridsearch.py milk/tests/test_grouped.py milk/tests/test_kmeans.py milk/tests/test_knn.py milk/tests/test_logistic.py milk/tests/test_measures.py milk/tests/test_measures_clusters.py milk/tests/test_multi.py milk/tests/test_multi_view.py milk/tests/test_nfoldcrossvalidation.py milk/tests/test_nfoldcrossvalidation_regression.py milk/tests/test_nnmf.py milk/tests/test_normalise.py milk/tests/test_normaliselabels.py milk/tests/test_parzen.py milk/tests/test_pca.py milk/tests/test_pdist.py milk/tests/test_precluster_learner.py milk/tests/test_regression.py milk/tests/test_regression_constant_features.py milk/tests/test_rf.py milk/tests/test_set2binary_array.py milk/tests/test_som.py milk/tests/test_svm.py milk/tests/test_tree.py milk/tests/test_utils.py milk/tests/data/__init__.py milk/tests/data/jugparallel_jugfile.py milk/tests/data/jugparallel_kmeans_jugfile.py milk/unsupervised/__init__.py milk/unsupervised/_kmeans.cpp milk/unsupervised/_som.cpp milk/unsupervised/affinity.py milk/unsupervised/gaussianmixture.py milk/unsupervised/kmeans.py milk/unsupervised/normalise.py milk/unsupervised/parzen.py milk/unsupervised/pca.py milk/unsupervised/pdist.py milk/unsupervised/som.py milk/unsupervised/nnmf/__init__.py milk/unsupervised/nnmf/hoyer.py milk/unsupervised/nnmf/lee_seung.py milk/utils/__init__.py milk/utils/utils.py milk/wrapper/__init__.py milk/wrapper/wraplibsvm.pymilk-0.3.10/milk.egg-info/top_level.txt0000644000175000017500000000000511562337213020503 0ustar luispedroluispedromilk milk-0.3.10/milk.egg-info/PKG-INFO0000644000175000017500000001341111562337213017053 0ustar luispedroluispedroMetadata-Version: 1.0 Name: milk Version: 0.3.10 Summary: Machine Learning Toolkit Home-page: http://luispedro.org/software/milk Author: Luis Pedro Coelho Author-email: lpc@cmu.edu License: MIT Description: ============================== MILK: MACHINE LEARNING TOOLKIT ============================== Machine Learning in Python -------------------------- Milk is a machine learning toolkit in Python. Its focus is on supervised classification with several classifiers available: SVMs (based on libsvm), k-NN, random forests, decision trees. It also performs feature selection. These classifiers can be combined in many ways to form different classification systems. For unsupervised learning, milk supports k-means clustering and affinity propagation. Milk is flexible about its inputs. It optimised for numpy arrays, but can often handle anything (for example, for SVMs, you can use any dataype and any kernel and it does the right thing). There is a strong emphasis on speed and low memory usage. Therefore, most of the performance sensitive code is in C++. This is behind Python-based interfaces for convenience. To learn more, check the docs at `http://packages.python.org/milk/ `_ or the code demos included with the source at ``milk/demos/``. Examples -------- Here is how to test how well you can classify some ``features,labels`` data, measured by cross-validation:: import numpy as np import milk features = np.random.rand(100,10) # 2d array of features: 100 examples of 10 features each labels = np.zeros(100) features[50:] += .5 labels[50:] = 1 confusion_matrix, names = milk.nfoldcrossvalidation(features, labels) print 'Accuracy:', confusion_matrix.trace()/float(confusion_matrix.sum()) If want to use a classifier, you instanciate a *learner object* and call its ``train()`` method:: import numpy as np import milk features = np.random.rand(100,10) labels = np.zeros(100) features[50:] += .5 labels[50:] = 1 learner = milk.defaultclassifier() model = learner.train(features, labels) # Now you can use the model on new examples: example = np.random.rand(10) print model.apply(example) example2 = np.random.rand(10) example2 += .5 print model.apply(example2) There are several classification methods in the package, but they all use the same interface: ``train()`` returns a *model* object, which has an ``apply()`` method to execute on new instances. Details ------- License: MIT Author: Luis Pedro Coelho (with code from LibSVM and scikits.learn) API Documentation: `http://packages.python.org/milk/ `_ Mailing List: `http://groups.google.com/group/milk-users `__ Features -------- - SVMs. Using the libsvm solver with a pythonesque wrapper around it. - K-means using as little memory as possible. It can cluster millions of instances efficiently. - Random forests - Self organising maps - Stepwise Discriminant Analysis for feature selection. - Non-negative matrix factorisation - Affinity propagation Recent History -------------- The ChangeLog file contains a more complete history. New in 0.3.10 ~~~~~~~~~~~~~ - Add ext.jugparallel for integration with `jug `_ - parallel nfold crossvalidation using jug - parallel multiple kmeans runs using jug - cluster_agreement for non-ndarrays - Add histogram & normali(z|s)e options to ``milk.kmeans.assign_centroid`` - Fix bug in sda when features were constant for a class - Add select_best_kmeans - Added defaultlearner as a better name than defaultclassifier - Add ``measures.curves.precision_recall`` - Add ``unsupervised.parzen.parzen`` New in 0.3.9 ~~~~~~~~~~~~ - Add ``folds`` argument to ``nfoldcrossvalidation`` - Add ``assign_centroid`` function in milk.unsupervised.nfoldcrossvalidation - Improve speed of k-nearest neighbour (10x on scikits-learn benchmark) - Improve kmeans on newer numpy (works for larger datasets too) - Faster kmeans by coding centroid recalculation in C++ - Fix gridminize for low count labels - Fix bug with non-integer labels for tree learning New in 0.3.8 ~~~~~~~~~~~~ - Fix compilation on Windows New in 0.3.7 ~~~~~~~~~~~~ - Logistic regression - Source demos included (in source and documentation) - Add cluster agreement metrics - Fix nfoldcrossvalidation bug when using origins New in 0.3.6 ~~~~~~~~~~~~ - Unsupervised (1-class) kernel density modeling - Fix for when SDA returns empty - weights option to some learners - stump learner - Adaboost (result of above changes) New in 0.3.5 ~~~~~~~~~~~~ - Fixes for 64 bit machines - Functions in measures.py all have same interface now. Platform: UNKNOWN milk-0.3.10/milk/0000755000175000017500000000000011562337213014264 5ustar luispedroluispedromilk-0.3.10/milk/ext/0000755000175000017500000000000011562337213015064 5ustar luispedroluispedromilk-0.3.10/milk/ext/__init__.py0000644000175000017500000000026411540474143017177 0ustar luispedroluispedro''' =============== Milk Extensions =============== These are modules whose functionality is not really part of the core functionality of milk, but which are useful with it. ''' milk-0.3.10/milk/ext/jugparallel.py0000644000175000017500000001211611561061601017733 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # License: MIT. See COPYING.MIT file in the milk distribution ''' ============ Jug Parallel ============ These are some functions that make it easier to take advantage of `jug `__ to perform tasks in parallel. All of the functions in this module return a jug.Task object and are usable with the ``CompoundTask`` interface in jug. ''' from __future__ import division import numpy as np import milk try: from jug import TaskGenerator, value from jug.utils import identity from jug.mapreduce import mapreduce except ImportError: raise ImportError('milk.ext.jugparallel requires jug (http://luispedro.org/software/jug') class _nfold_one(object): def __init__(self, features, labels, kwargs): self.features = features self.labels = labels self.kwargs = kwargs def __call__(self, i): return milk.nfoldcrossvalidation( value(self.features), value(self.labels), folds=[i], **value(self.kwargs)) def __jug_hash__(self): # jug.hash is only available in jug 0.9 # This is also the first version that would call __jug_hash__ # So, we import it here only. from jug import hash M = hash.new_hash_object() hash.hash_update(M,[ ('type', 'milk.nfoldcrossvalidation'), ('features', self.features), ('labels', self.labels), ('kwargs', self.kwargs), ]) hexdigest = M.hexdigest() self.__jug_hash__ = lambda: hexdigest return hexdigest def _nfold_reduce(a,b): cmat = a[0] + b[0] names = a[1] assert a[1] == b[1] if len(a) == 2: return cmat, names predictions = np.array([a[2],b[2]]) return cmat, names, predictions.max(0) def nfoldcrossvalidation(features, labels, **kwargs): ''' jug_task = nfoldcrossvalidation(features, labels, **kwargs) A jug Task that perform n-foldcrossvalidation N-fold cross validation is inherently parallel. This function returns a ``jug.Task`` which performs n-fold crossvalidation which jug can parallelise. Parameters ---------- features : sequence of features labels : sequence kwargs : any This will be passed down to ``milk.nfoldcrossvalidation`` Returns ------- jug_task : a jug.Task A Task object See Also -------- milk.nfoldcrossvalidation : The same functionality as a "normal" function jug.CompoundTask : This function can be used as argument to CompoundTask ''' mapper = _nfold_one(features, labels, kwargs) nfolds = kwargs.get('nfolds', 10) return mapreduce(_nfold_reduce, mapper, range(nfolds), map_step=1, reduce_step=(nfolds+1)) def _select_min(s0, s1): if s0[0] < s1[0]: return s0 else: return s1 def _evaluate_solution(args): features, results, method = args from milk.unsupervised.gaussianmixture import AIC, BIC if method == 'AIC': method = AIC elif method == 'BIC': method = BIC else: raise ValueError('milk.ext.jugparallel.kmeans_select_best: unknown method: %s' % method) assignments, centroids = results value = method(features, assignments, centroids) return value, results def _select_best(features, results, method): features = identity(features) return mapreduce(_select_min, _evaluate_solution, [(features,r,method) for r in results], reduce_step=32, map_step=8) def kmeans_select_best(features, ks, repeats=1, method='AIC', R=None, **kwargs): ''' assignments_centroids = kmeans_select_best(features, ks, repeats=1, method='AIC', R=None, **kwargs) Perform ``repeats`` calls to ``kmeans`` for each ``k`` in ``ks``, select the best one according to ``method.`` Note that, unlike a raw ``kmeans`` call, this is *always deterministic* even if ``R=None``. Parameters ---------- features : array-like 2D array ks : sequence of integers These will be the values of ``k`` to try repeats : integer, optional How many times to attempt each k (default: 1). method : str, optional Which method to use. Must be one of 'AIC' (default) or 'BIC'. R : random number source, optional If you do not pass a value, the result will be deterministic kwargs : other options These are passed transparently to ``kmeans`` Returns ------- assignments_centroids : jug.Task jug.Task which is the result of the best (as measured by ``method``) kmeans clustering. ''' from milk import kmeans from milk.utils import get_pyrandom kmeans = TaskGenerator(kmeans) if R is not None: start = get_pyrandom(R).randint(0,1024*1024) else: start = 7 results = [] for ki,k in enumerate(ks): for i in xrange(repeats): results.append(kmeans(features, k, R=(start+7*repeats*ki+i), **kwargs)) return _select_best(features, results, method)[1] milk-0.3.10/milk/nfoldcrossvalidation.py0000644000175000017500000000013011452411266021056 0ustar luispedroluispedrofrom .measures.nfoldcrossvalidation import foldgenerator, getfold, nfoldcrossvalidation milk-0.3.10/milk/measures/0000755000175000017500000000000011562337213016110 5ustar luispedroluispedromilk-0.3.10/milk/measures/nfoldcrossvalidation.py0000644000175000017500000001564211540476161022723 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # LICENSE: MIT from __future__ import division from ..supervised.classifier import normaliselabels import numpy as np __all__ = ['foldgenerator', 'getfold', 'nfoldcrossvalidation'] def foldgenerator(labels, nfolds=None, origins=None, folds=None): ''' for train,test in foldgenerator(labels, nfolds=None, origins=None) ... This generator breaks up the data into `n` folds (default 10). If `origins` is given, then all elements that share the same origin will either be in testing or in training (never in both). This is useful when you have several replicates that shouldn't be mixed together between training&testing but that can be otherwise be treated as independent for learning. Parameters ---------- labels : a sequence the labels nfolds : integer nr of folds (default 10 or minimum label size) origins : sequence, optional if present, must be an array of indices of the same size as labels. folds : sequence of int, optional which folds to generate Returns ------- iterator over `train, test`, two boolean arrays ''' labels,names = normaliselabels(labels) if origins is None: origins = np.arange(len(labels)) else: if len(origins) != len(labels): raise ValueError( 'milk.nfoldcrossvalidation.foldgenerator: origins must be of same size as labels') origins = np.asanyarray(origins) fmin = len(labels) for lab in set(labels): curmin = len(set(origins[labels == lab])) fmin = min(fmin, curmin) if fmin == 1: raise ValueError(''' milk.nfoldcrossvalidation.foldgenerator: nfolds was reduced to 1 because minimum class size was 1. If you passed in an origins parameter, it might be caused by having a class come from a single origin. ''') fold = np.zeros(len(labels)) fold -= 1 if nfolds is None: nfolds = min(fmin, 10) elif nfolds > fmin: from warnings import warn warn('milk.measures.nfoldcrossvalidation: Reducing the nr. of folds from %s to %s (minimum class size).' % (nfolds, fmin)) nfolds = fmin for lab in set(labels): locations = (labels == lab) usedorigins = np.unique(origins[locations]) weights = np.array([np.sum(origins == orig) for orig in usedorigins]) foldweight = np.zeros(nfolds, int) for w,orig in sorted(zip(weights, usedorigins)): f = np.argmin(foldweight) if np.any(fold[origins == orig] > -1): raise ValueError( 'milk.nfoldcrossvalidation.foldgenerator: something is wrong. Maybe origin %s is present in two labels.' % orig) fold[origins == orig] = f foldweight[f] += w for f in xrange(nfolds): if folds is not None and f not in folds: continue yield (fold != f), (fold == f) def getfold(labels, fold, nfolds=None, origins=None): ''' trainingset,testingset = getfold(labels, fold, nfolds=None, origins=None) Get the training and testing set for fold `fold` in `nfolds` Arguments are the same as for `foldgenerator` Parameters ---------- labels : ndarray of labels fold : integer nfolds : integer number of folds (default 10 or size of smallest class) origins : sequence, optional if given, then objects with same origin are *not* scattered across folds ''' if nfolds < fold: raise ValueError('milk.getfold: Attempted to get fold %s out of %s' % (fold, nfolds)) for i,(t,s) in enumerate(foldgenerator(labels, nfolds, origins)): if i == fold: return t,s raise ValueError('milk.getfold: Attempted to get fold %s but the number of actual folds was too small (%s)' % (fold,i)) def nfoldcrossvalidation(features, labels, nfolds=None, learner=None, origins=None, return_predictions=False, folds=None, classifier=None): ''' Perform n-fold cross validation cmatrix,names = nfoldcrossvalidation(features, labels, nfolds=10, learner={defaultclassifier()}, origins=None, return_predictions=False) cmatrix,names,predictions = nfoldcrossvalidation(features, labels, nfolds=10, learner={defaultclassifier()}, origins=None, return_predictions=True) cmatrix will be a N x N matrix, where N is the number of classes cmatrix[i,j] will be the number of times that an element of class i was classified as class j names[i] will correspond to the label name of class i Parameters ---------- features : a sequence labels : an array of labels, where label[i] is the label corresponding to features[i] nfolds : integer, optional Nr of folds. Default: 10 learner : learner object, optional learner should implement the train() method to return a model (something with an apply() method). defaultclassifier() by default This parameter used to be called `classifier` and that name is still supported origins : sequence, optional Origin ID (see foldgenerator) return_predictions : bool, optional whether to return predictions (default: False) folds : sequence of int, optional which folds to generate Returns ------- cmatrix : ndarray confusion matrix names : sequence sequence of labels so that cmatrix[i,j] corresponds to names[i], names[j] predictions : sequence predicted output for each element ''' import operator from .measures import confusion_matrix if len(features) != len(labels): raise ValueError('milk.measures.nfoldcrossvalidation: len(features) should match len(labels)') if classifier is not None: if learner is not None: raise ValueError('milk.nfoldcrossvalidation: Using both `learner` and `classifier` arguments. They are the same, but `learner` is preferred') learner = classifier if learner is None: from ..supervised.defaultclassifier import defaultclassifier learner = defaultclassifier() labels,names = normaliselabels(labels) if return_predictions: predictions = np.empty_like(labels) predictions.fill(-1) # This makes it clearer if there are bugs in the programme features = np.asanyarray(features) nclasses = labels.max() + 1 results = [] measure = confusion_matrix for trainingset,testingset in foldgenerator(labels, nfolds, origins=origins, folds=folds): model = learner.train(features[trainingset], labels[trainingset]) cur_preds = np.array([model.apply(f) for f in features[testingset]]) if return_predictions: predictions[testingset] = cur_preds results.append(measure(labels[testingset], cur_preds)) result = reduce(operator.add, results) if return_predictions: return result, names, predictions return result, names milk-0.3.10/milk/measures/cluster_agreement.py0000644000175000017500000000322211540761772022200 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division import numpy as np def rand_arand_jaccard(recovered, labels): ''' rand, a_rand, jaccard = rand_arand_jaccard(recovered, labels) Compute Rand, Adjusted Rand, and Jaccard indices These share most of the computation. Therefore, it is best to compute them together even if you are only going to use some. Parameters ---------- recovered : sequence of int The recovered clusters labels : sequence of int Underlying labels Returns ------- rand : float Rand index a_rand : float Adjusted Rand index jaccard : float Jaccard index References ---------- http://en.wikipedia.org/wiki/Rand_index http://en.wikipedia.org/wiki/Jaccard_index ''' from scipy.misc import comb recovered = np.asanyarray(recovered) labels = np.asanyarray(labels) contig,_,_ = np.histogram2d(recovered, labels,np.arange(max(recovered.max()+2,labels.max()+2))) A_0 = contig.sum(0) A_1 = contig.sum(1) Ai2 = np.sum(A_0*(A_0-1)/2.) Bi2 = np.sum(A_1*(A_1-1)/2.) n = A_0.sum() a = comb(contig.ravel(), 2).sum() b = comb(A_0, 2).sum()-a c = comb(A_1, 2).sum()-a d = comb(n, 2)-a-b-c rand = (a+d)/(a+b+c+d) jaccard = (a+d)/(b+c+d) index = np.sum(contig*(contig-1)/2) expected = Ai2*Bi2/n/(n-1)*2. maxindex = (Ai2+Bi2)/2. a_rand = (index-expected)/(maxindex-expected) return rand, a_rand, jaccard milk-0.3.10/milk/measures/measures.py0000644000175000017500000000722511467264577020334 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2010, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT from __future__ import division import numpy as np __all__ = [ 'accuracy', 'confusion_matrix', 'waccuracy', 'zero_one_loss', ] def accuracy(real, other=None, normalisedlabels=False, names=None): ''' acc = accuracy(real, predicted, normalisedlabels=False, names=None) Compute accuracy (fraction of correct predictions). Parameters ---------- real : sequence The real labels predicted : sequence The predicted sequence (must be same type as `real`) normalisedlabels : boolean, optional Whether labels have been normalised names : sequence The underlying names (unused) Returns ------- acc : float ''' if other is None: import warnings warnings.warn('milk.measures.accuracy: calling this with one argument is a deprecated interface.', DeprecationWarning) cmatrix = np.asanyarray(real) return cmatrix.trace()/cmatrix.sum() else: real = cmatrix return np.mean(np.asanyarray(real) != other) def zero_one_loss(real, predicted, normalisedlabels=False, names=None): ''' loss = zero_one_loss(real, predicted, normalisedlabels={unused}, names={unused}) Parameters ---------- real : sequence the underlying labels predicted : sequence the predicted labels normalisedlabels : unused names: unused Returns ------- loss : integer the number of instances where `real` differs from `predicted` ''' return np.sum(np.asanyarray(real) != np.asanyarray(predicted)) def waccuracy(real, predicted=None, normalisedlabels=False, names=None): ''' wacc = waccuracy(real, predicted, normalisedlabels={unused}, names={unused}) Weighted accuracy: average of accuracy for each (real) class. Can be very different from accuracy if the classes are unbalanced (in particular, if they are very unbalanced, you can get a high accuracy with a bad classifier). Parameters ---------- real : sequence the underlying labels predicted : sequence the predicted labels normalisedlabels : unused names: unused Returns ------- wacc : float the weighted accuracy ''' if predicted is None: import warnings warnings.warn('milk.measures.accuracy: calling this with one argument is a deprecated interface.', DeprecationWarning) cmatrix = np.asanyarray(real) else: cmatrix = confusion_matrix(real, predicted, normalisedlabels, names) return (cmatrix.diagonal() / cmatrix.sum(1)).mean() def confusion_matrix(real, predicted, normalisedlabels=False, names=None): ''' cmatrix = confusion_matrix(real, predicted, normalisedlabels=False, names=None) Computes the confusion matrix Parameters ---------- real : sequence The real labels predicted : sequence The predicted sequence (must be same type as `real`) normalisedlabels : boolean, optional Whether labels have been normalised names : sequence The underlying names (unused) Returns ------- cmatrix : 2 ndarray ''' if not normalisedlabels: from ..supervised.normalise import normaliselabels real, names = normaliselabels(real) predicted = map(names.index, predicted) n = np.max(real)+1 cmat = np.zeros((n,n), int) for r,p in zip(real, predicted): cmat[r,p] += 1 return cmat ## TODO: Implement http://en.wikipedia.org/wiki/Matthews_Correlation_Coefficient milk-0.3.10/milk/measures/curves.py0000644000175000017500000000373511561061601017773 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division import numpy as np def precision_recall(values, labels, mode='all', nr_steps=100): ''' precision, recall = precision_recall(values, labels, mode='all', nr_steps=100) Compute a precision-recall curve. For a given threshold ``T``, consider that the positions where ``values >= T`` are classified as True. Precision is defined as ``TP/(TP+FP)``, while recall is defined as ``TP/(TP+FN)``. Parameters ---------- values : sequence of numbers labels : boolean sequence mode : str, optional Which thresholds to consider. Either 'all' (i.e., use all values of `values` as possible thresholds), or 'step' (using `nr_steps` equidistant points from ``min(values)`` to ``max(values)``) nr_steps : integer, optional How many steps to use. Only meaningfule if ``mode == 'steps'`` Returns ------- precision : a sequence of floats recall : a sequence of floats Actually, ``2 x P`` array is returned. ''' values = np.asanyarray(values) labels = np.asanyarray(labels) if len(values) != len(labels): raise ValueError('milk.measures.precision_recall: `values` must be of same length as `labels`') if mode == 'all': points = list(set(values)) points.sort() elif mode == 'steps': points = np.linspace(values.min(), values.max(), nr_steps) else: raise ValueError('milk.measures.precision_recall: cannot handle mode: `%s`' % mode) true_pos = float(np.sum(labels)) precision_recall = np.empty((len(points),2), np.float) for i,p in enumerate(points): selected = (values >= p) selected = labels[selected] precision_recall[i] = (np.mean(selected), np.sum(selected)/true_pos) return precision_recall.T milk-0.3.10/milk/measures/__init__.py0000644000175000017500000000011211467264577020233 0ustar luispedroluispedrofrom measures import accuracy, waccuracy, zero_one_loss, confusion_matrix milk-0.3.10/milk/supervised/0000755000175000017500000000000011562337213016455 5ustar luispedroluispedromilk-0.3.10/milk/supervised/parzen.py0000644000175000017500000000136111511426617020330 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division import numpy as np def get_parzen_rbf_loocv(features,labels): xij = np.dot(features,features.T) f2 = np.sum(features**2,1) d = f2-2*xij d = d.T + f2 d_argsorted = d.argsort(1) d_sorted = d.copy() d_sorted.sort(1) e_d = np.exp(-d_sorted) labels_sorted = labels[d_argsorted].astype(np.double) labels_sorted *= 2 labels_sorted -= 1 def f(sigma): k = e_d ** (1./sigma) return (((k[:,1:] * labels_sorted[:,1:]).sum(1) > 0) == labels).mean() return f milk-0.3.10/milk/supervised/svm.py0000644000175000017500000003427511511426617017650 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division from .classifier import normaliselabels, ctransforms_model from collections import deque import numpy import numpy as np import random from . import _svm __all__ = [ 'rbf_kernel', 'polynomial_kernel', 'precomputed_kernel', 'dot_kernel', 'svm_raw', 'svm_binary', 'svm_to_binary', 'svm_sigmoidal_correction', 'sigma_value_fisher', 'fisher_tuned_rbf_svm', ] def _svm_apply(SVM, q): ''' f_i = _svm_apply(SVM, q) @internal: This is mostly used for testing ''' X,Y,Alphas,b,C,kernel=SVM N = len(X) s = 0.0 for i in xrange(N): s += Alphas[i] * Y[i] * kernel(q, X[i]) return s - b def svm_learn_smo(X,Y,kernel,C,eps=1e-4,tol=1e-2,cache_size=(1<<20)): ''' Learn a svm classifier X: data Y: labels in SVM format (ie Y[i] in (1,-1)) This is a very raw interface. In general, you should use a class like svm_classifier. Implements the Sequential Minimum Optimisation Algorithm from Platt's "Fast training of support vector machines using sequential minimal optimization" in Advances in kernel methods: support vector learning Pages: 185 - 208 Year of Publication: 1999 ISBN:0-262-19416-3 ''' assert numpy.all(numpy.abs(Y) == 1) assert len(X) == len(Y) N = len(Y) Y = Y.astype(numpy.int32) params = numpy.array([0,C,1e-3,1e-5],numpy.double) Alphas0 = numpy.zeros(N, numpy.double) _svm.eval_SMO(X,Y,Alphas0,params,kernel,cache_size) return Alphas0, params[0] def svm_learn_libsvm(X,Y,kernel,C,eps=1e-4,tol=1e-2,cache_size=(1<<20)): ''' Learn a svm classifier using LIBSVM optimiser X: data Y: labels in SVM format (ie Y[i] in (1,-1)) This is a very raw interface. In general, you should use a class like svm_classifier. This uses the LIBSVM optimisation algorithm ''' assert numpy.all(numpy.abs(Y) == 1) assert len(X) == len(Y) N = len(Y) Y = Y.astype(numpy.int32) p = -numpy.ones(N,numpy.double) params = numpy.array([0,C,eps,tol],numpy.double) Alphas0 = numpy.zeros(N, numpy.double) _svm.eval_LIBSVM(X,Y,Alphas0,p,params,kernel,cache_size) return Alphas0, params[0] class preprocessed_rbf_kernel(object): def __init__(self, X, sigma, beta): self.X = X self.Xsum = (X**2).sum(1) self.sigma = sigma self.beta = beta def __call__(self, q): minus_d2_sigma = np.dot(self.X,q) minus_d2_sigma *= 2. minus_d2_sigma -= self.Xsum minus_d2_sigma -= np.dot(q,q) minus_d2_sigma /= self.sigma return self.beta * np.exp(minus_d2_sigma) class rbf_kernel(object): ''' kernel = rbf_kernel(sigma,beta=1) Radial Basis Function kernel Returns a kernel (ie, a function that implements) beta * exp( - ||x1 - x2|| / sigma) ''' def __init__(self, sigma, beta=1): self.sigma = sigma self.beta = beta self.kernel_nr_ = 0 self.kernel_arg_ = float(sigma) def __call__(self, x1, x2): d2 = x1 - x2 d2 **= 2 d2 = d2.sum() res = self.beta*np.exp(-d2/self.sigma) return res def preprocess(self, X): return preprocessed_rbf_kernel(X, self.sigma, self.beta) class polynomial_kernel(object): ''' kernel = polynomial_kernel(d,c=1) returns a kernel (ie, a function) that implements: (+c)**d ''' def __init__(self, d, c=1): self.d = d self.c = c def __call__(self,x1,x2): return (np.dot(x1,x2)+self.c)**self.d class precomputed_kernel(object): ''' kernel = precomputed_kernel(kmatrix) A "fake" kernel which is precomputed. ''' def __init__(self, kmatrix, copy=False): kmatrix = np.ascontiguousarray(kmatrix, np.double, copy=copy) self.kernel_nr_ = 1 self.kernel_arg_ = 0. def __call__(self, x0, x1): return kmatrix[x0,x1] class _call_kernel(object): def __init__(self, k, svs): self.svs = svs self.kernel = k def __call__(self, q): return np.array([self.kernel(s, q) for s in self.svs]) class dot_kernel(object): def __init__(self): self.kernel_nr_ = 2 self.kernel_arg_ = 0. def __call__(self, x0, x1): return np.dot(x0, x1) class svm_raw_model(object): def __init__(self, svs, Yw, b, kernel): self.svs = svs self.Yw = Yw self.b = b self.kernel = kernel try: self.kernelfunction = self.kernel.preprocess(self.svs) except AttributeError: self.kernelfunction = _call_kernel(self.kernel, self.svs) def apply(self, q): Q = self.kernelfunction(q) return np.dot(Q, self.Yw) - self.b class svm_raw(object): ''' svm_raw: classifier classifier = svm_raw(kernel, C, eps=1e-3, tol=1e-8) Parameters ---------- kernel : the kernel to use. This should be a function that takes two data arguments see rbf_kernel and polynomial_kernel. C : the C parameter Other Parameters ---------------- eps : the precision to which to solve the problem (default 1e-3) tol : (|x| < tol) is considered zero ''' def __init__(self, kernel=None, C=1., eps=1e-3, tol=1e-8): self.C = C self.kernel = kernel self.eps = eps self.tol = tol self.algorithm = 'libsvm' def train(self, features, labels, normalisedlabels=False): assert self.kernel is not None, 'milk.supervised.svm_raw.train: kernel not set!' assert self.algorithm in ('libsvm','smo'), 'milk.supervised.svm_raw: unknown algorithm (%s)' % self.algorithm assert not (np.isinf(self.C) or np.isnan(self.C)), 'milk.supervised.svm_raw: setting C to NaN or Inf causes problems.' features = np.asanyarray(features) if normalisedlabels: Y = labels else: Y,_ = normaliselabels(labels) assert Y.max() == 1, 'milk.supervised.svm_raw can only handle binary problems' Y *= 2 Y -= 1 kernel = self.kernel try: kernel = (self.kernel.kernel_nr_, self.kernel.kernel_arg_) features = np.ascontiguousarray(features, np.double) except AttributeError: pass if self.algorithm == 'smo': alphas,b = svm_learn_smo(features,Y,kernel,self.C,self.eps,self.tol) else: alphas,b = svm_learn_libsvm(features,Y,kernel,self.C,self.eps,self.tol) svsi = (alphas != 0) svs = features[svsi] w = alphas[svsi] Y = Y[svsi] Yw = w * Y return svm_raw_model(svs, Yw, b, self.kernel) def get_params(self): return self.C, self.eps,self.tol def set_params(self,params): self.C,self.eps,self.tol = params def set_option(self, optname, value): setattr(self, optname, value) def learn_sigmoid_constants(F,Y, max_iters=None, min_step=1e-10, sigma=1e-12, eps=1e-5): ''' A,B = learn_sigmoid_constants(F,Y) This is a very low-level interface look into the svm_classifier class. Parameters ---------- F : Values of the function F Y : Labels (in boolean format, ie, in (0,1)) Other Parameters ---------------- max_iters : Maximum nr. of iterations min_step : Minimum step sigma : sigma eps : A small number Reference for Implementation ---------------------------- Implements the algorithm from "A Note on Platt's Probabilistic Outputs for Support Vector Machines" by Lin, Lin, and Weng. Machine Learning, Vol. 68, No. 3. (23 October 2007), pp. 267-276 ''' # the deci[i] array is called F[i] in this code F = np.asanyarray(F) Y = np.asanyarray(Y) assert len(F) == len(Y) assert numpy.all( (Y == 1) | (Y == 0) ) from numpy import log, exp N=len(F) if max_iters is None: max_iters = 1000 prior1 = Y.sum() prior0 = N-prior1 small_nr = 1e-4 hi_t = (prior1+1.)/(prior1+2.) lo_t = 1./(prior0+2.) T = Y*hi_t + (1-Y)*lo_t A = 0. B = log( (prior0+1.)/(prior1+1.) ) def target(A,B): fval = 0. for i in xrange(N): fApB = F[i]*A+B if fApB >= 0: fval += T[i]*fApB+log(1+exp(-fApB)) else: fval += (T[i]-1.)*fApB+log(1+exp(fApB)) return fval fval = target(A,B) for iter in xrange(max_iters): h11=sigma h22=sigma h21=0. g1=0. g2=0. for i in xrange(N): fApB = F[i]*A+B if (fApB >= 0): p = exp(-fApB)/(1.+exp(-fApB)) q = 1./(1.+exp(-fApB)) else: p = 1./(1.+exp(fApB)) q = exp(fApB)/(1.+exp(fApB)) d2 = p * q h11 += F[i]*F[i]*d2 h22 += d2 h21 += F[i]*d2 d1 = T[i] - p g1 += F[i]*d1 g2 += d1 if abs(g1) < eps and abs(g2) < eps: # Stopping criteria break det = h11*h22 - h21*h21 dA = - (h22*g1 - h21*g2)/det dB = - (h21*g1 + h11*g2)/det gd = g1*dA + g2*dB stepsize = 1. while stepsize >= min_step: newA = A + stepsize*dA newB = B + stepsize*dB newf = target(newA,newB) if newf < fval+eps*stepsize*gd: A = newA B = newB fval = newf break stepsize /= 2 if stepsize < min_step: print 'Line search fails' break return A,B class svm_binary_model(object): def __init__(self, classes): self.classes = classes def apply(self,f): return self.classes[f >= 0.] class svm_binary(object): ''' classifier = svm_binary() model = classifier.train(features, labels) assert model.apply(f) in labels ''' def train(self, features, labels, normalisedlabels=False): if normalisedlabels: return svm_binary_model( (0,1) ) assert len(labels) >= 2, 'Cannot train from a single example' names = sorted(set(labels)) assert len(names) == 2, 'milk.supervised.svm.svm_binary.train: Can only handle two class problems' return svm_binary_model(names) class svm_to_binary(object): ''' svm_to_binary(base_svm) A simple wrapper so that svm_to_binary(base_svm) is a model that takes the base_svm classifier and then binarises its model output. NOTE: This class does the same job as:: ctransforms(base_svm, svm_binary()) ''' def __init__(self, svm_base): ''' binclassifier = svm_to_binary(svm_base) a classifier that binarises the output of svm_base. ''' self.base = svm_base def train(self, features, labels, normalisedlabels=False): model = self.base.train(features, labels, normalisedlabels=normalisedlabels) binary = svm_binary() binary_model = binary.train(features, labels, normalisedlabels=normalisedlabels) return ctransforms_model([model, binary_model]) def set_option(self, opt, value): self.base.set_option(opt, value) class svm_sigmoidal_correction_model(object): def __init__(self, A, B): self.A = A self.B = B def apply(self,features): return 1./(1.+numpy.exp(features*self.A+self.B)) class svm_sigmoidal_correction(object): ''' svm_sigmoidal_correction : a classifier Sigmoidal approximation for obtaining a probability estimate out of the output of an SVM. ''' def __init__(self): self.max_iters = None def train(self, features, labels, normalisedlabels=False): A,B = learn_sigmoid_constants(features,labels,self.max_iters) return svm_sigmoidal_correction_model(A, B) def get_params(self): return self.max_iters def set_params(self,params): self.max_iters = params def sigma_value_fisher(features,labels): ''' f = sigma_value_fisher(features,labels) value_s = f(s) Computes a function which computes how good the value of sigma is for the features. This function should be *minimised* for a good value of sigma. Parameters ----------- features : features matrix as 2-ndarray. Returns ------- f : a function: float -> float this function should be minimised for a good `sigma` Reference ---------- Implements the measure in "Determination of the spread parameter in the Gaussian kernel for classification and regression" by Wenjian Wanga, Zongben Xua, Weizhen Luc, and Xiaoyun Zhanga ''' features = np.asanyarray(features) xij = np.dot(features,features.T) f2 = np.sum(features**2,1) d = f2-2*xij d = d.T + f2 N1 = (labels==0).sum() N2 = (labels==1).sum() C1 = -d[labels == 0][:,labels == 0] C2 = -d[labels == 1][:,labels == 1] C12 = -d[labels == 0][:,labels == 1] C1 = C1.copy() C2 = C2.copy() C12 = C12.copy() def f(sigma): sigma = float(sigma) N1 = C1.shape[0] N2 = C2.shape[0] if C12.shape != (N1,N2): raise ValueError C1v = np.sum(np.exp(C1/sigma))/N1 C2v = np.sum(np.exp(C2/sigma))/N2 C12v = np.sum(np.exp(C12/sigma))/N1/N2 return (N1 + N2 - C1v - C2v)/(C1v/N1+C2v/N2 - 2.*C12v) return f class fisher_tuned_rbf_svm(object): ''' F = fisher_tuned_rbf_svm(sigmas, base) Returns a wrapper classifier that uses RBF kernels automatically tuned using sigma_value_fisher. ''' def __init__(self, sigmas, base): self.sigmas = sigmas self.base = base def train(self, features, labels, normalisedlabels=False): f = sigma_value_fisher(features, labels) fs = [f(s) for s in self.sigmas] self.sigma = self.sigmas[np.argmin(fs)] self.base.set_option('kernel',rbf_kernel(self.sigma)) return self.base.train(features,labels) milk-0.3.10/milk/supervised/_tree.cpp0000644000175000017500000001162011525274201020253 0ustar luispedroluispedro// Copyright (C) 2010, Luis Pedro Coelho // License: MIT #include #include #include #include extern "C" { #include #include } namespace { template double set_entropy(const T* data, const int N, double* counts, long clen) { for (int i = 0; i != clen; ++i) counts[i] = 0.; for (int i = 0; i != N; ++i) { int value = data[i]; if (value >= clen) { PyErr_SetString(PyExc_RuntimeError, "milk.supervised.tree.set_entropy: label value too large. aborting"); return 0; } counts[value] += 1.; } // Here is the formula we use: // // H = - \sum px \log(px) // = - \sum (cx/N) \log( cx / N) // = - 1/N \sum { cx [ \log cx - \log N ] } // = - 1/N { (\sum cx \log cx ) - ( \sum cx \log N ) } // = - 1/N { (\sum cx \log cx ) - N \log N } // = ( - 1/N \sum cx \log cx ) + \log N double entropy = 0.; for (int i = 0; i != clen; ++i) { double cx = counts[i]; if (cx) entropy += cx * std::log(cx); } entropy /= -N; entropy += std::log(double(N)); return entropy; } double set_entropy(PyArrayObject* points, double* counts, long clen) { const int N = PyArray_SIZE(points); #define HANDLE(type_marker, type) \ if (PyArray_TYPE(points) == type_marker) return set_entropy(reinterpret_cast(PyArray_DATA(points)), N, counts, clen); HANDLE(NPY_INT, int) HANDLE(NPY_LONG, long) #undef HANDLE assert(false); return 0.; } PyObject* py_set_entropy(PyObject* self, PyObject* args) { const char* errmsg = "Arguments were not what was expected for set_entropy.\n" "This is an internal function: Do not call directly unless you know exactly what you're doing.\n"; PyArrayObject* labels; PyArrayObject* counts; if (!PyArg_ParseTuple(args, "OO", &labels, &counts)) { PyErr_SetString(PyExc_RuntimeError,errmsg); return 0; } if (!PyArray_Check(labels) || (PyArray_TYPE(labels) != NPY_INT && PyArray_TYPE(labels) != NPY_LONG) || !PyArray_ISCONTIGUOUS(labels) || !PyArray_Check(counts) || PyArray_TYPE(counts) != NPY_DOUBLE || !PyArray_ISCONTIGUOUS(counts)) { PyErr_SetString(PyExc_RuntimeError,errmsg); return 0; } double* cdata = reinterpret_cast(PyArray_DATA(counts)); const int nlabels = PyArray_DIM(counts, 0); double res = set_entropy(labels, cdata, nlabels); return PyFloat_FromDouble(res); } PyObject* py_information_gain(PyObject* self, PyObject* args) { const char* errmsg = "Arguments were not what was expected for information_gain.\n" "This is an internal function: Do not call directly unless you know exactly what you're doing.\n"; PyArrayObject* labels0; PyArrayObject* labels1; if (!PyArg_ParseTuple(args, "OO", &labels0, &labels1)) { PyErr_SetString(PyExc_RuntimeError,errmsg); return 0; } if (!PyArray_Check(labels0) || (PyArray_TYPE(labels0) != NPY_INT && PyArray_TYPE(labels0) != NPY_LONG) || !PyArray_ISCONTIGUOUS(labels0) || !PyArray_Check(labels1) || (PyArray_TYPE(labels1) != NPY_INT && PyArray_TYPE(labels1) != NPY_LONG) || !PyArray_ISCONTIGUOUS(labels1)) { PyErr_SetString(PyExc_RuntimeError,errmsg); std::cout << PyArray_TYPE(labels0) << " " << PyArray_TYPE(labels1) << '\n'; return 0; } double* counts; double counts_array[8]; long clen = 0; #define GET_MAX(index) \ const int N ## index = PyArray_DIM(labels ## index, 0); \ { \ for (int i = 0; i != N ## index; ++i) { \ long val = 0; \ if (PyArray_TYPE(labels ## index) == NPY_INT) val = *reinterpret_cast(PyArray_GETPTR1(labels ## index, i)); \ else if (PyArray_TYPE(labels ## index) == NPY_LONG) val = *reinterpret_cast(PyArray_GETPTR1(labels ## index, i)); \ if (val > clen) clen = val; \ } \ } GET_MAX(0); GET_MAX(1); ++clen; if (clen > 8) { counts = new(std::nothrow) double[clen]; if (!counts) { PyErr_NoMemory(); return 0; } } else { counts = counts_array; } const double N = N0 + N1; double H = - N0/N * set_entropy(labels0, counts, clen) - N1/N * set_entropy(labels1, counts, clen); if (clen > 8) delete [] counts; return PyFloat_FromDouble(H); } PyMethodDef methods[] = { {"set_entropy", py_set_entropy, METH_VARARGS , "Do NOT call directly.\n" }, {"information_gain", py_information_gain, METH_VARARGS , "Do NOT call directly.\n" }, {NULL, NULL,0,NULL}, }; const char * module_doc = "Internal Module.\n" "\n" "Do NOT use directly!\n"; } // namespace extern "C" void init_tree() { import_array(); (void)Py_InitModule3("_tree", methods, module_doc); } milk-0.3.10/milk/supervised/knn.py0000644000175000017500000000304711534275170017623 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division from collections import defaultdict import numpy as np __all__ = [ 'kNN', ] class kNN(object): ''' k-Nearest Neighbour Classifier Naive implementation of a k-nearest neighbour classifier. C = kNN(k) Attributes: ----------- k : integer number of neighbours to consider ''' def __init__(self, k=1): self.k = k def train(self, features, labels, normalisedlabels=False, copy_features=False): features = np.asanyarray(features) labels = np.asanyarray(labels) if copy_features: features = features.copy() labels = labels.copy() features2 = np.sum(features**2, axis=1) return kNN_model(self.k, features, features2, labels) class kNN_model(object): def __init__(self, k, features, features2, labels): self.k = k self.features = features self.f2 = features2 self.labels = labels def apply(self, features): features = np.asanyarray(features) diff2 = np.dot(self.features, (-2.)*features) diff2 += self.f2 neighbours = diff2.argsort()[:self.k] labels = self.labels[neighbours] votes = defaultdict(int) for L in labels: votes[L] += 1 v,L = max( (v,L) for L,v in votes.items() ) return L milk-0.3.10/milk/supervised/precluster_learner.py0000644000175000017500000000245011556634056022740 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division import numpy as np from milk.unsupervised.kmeans import select_best_kmeans, assign_centroids from milk import defaultlearner class precluster_model(object): def __init__(self, centroids, base): self.centroids = centroids self.base = base def apply(self, features): histogram = assign_centroids(features, self.centroids, histogram=True, normalise=1) return self.base.apply(histogram) class precluster_learner(object): ''' This learns a classifier by clustering the input features ''' def __init__(self, ks, base=None): if base is None: base = defaultlearner() self.ks = ks self.base = base def train(self, features, labels, **kwargs): allfeatures = np.concatenate(features) assignments, centroids = select_best_kmeans(allfeatures, self.ks, 1, "AIC") histograms = [assign_centroids(f, centroids, histogram=True, normalise=1) for f in features] base_model = self.base.train(histograms, labels, **kwargs) return precluster_model(centroids, base_model) milk-0.3.10/milk/supervised/randomforest.py0000644000175000017500000000450211514100570021522 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2010-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT. See COPYING.MIT file in the milk distribution ''' Random Forest ------------- Main elements ------------- rf_learner : A learner object ''' from __future__ import division import numpy as np import milk.supervised.tree from .normalise import normaliselabels __all__ = [ 'rf_learner', ] def _sample(features, labels, n, R): ''' features', labels' = _sample(features, labels, n, R) Sample n element from (features,labels) Parameters ---------- features : sequence labels : sequence Same size as labels n : integer R : random object Returns ------- features' : sequence labels' : sequence ''' N = len(features) sfeatures = [] slabels = [] for i in xrange(n): idx = R.randint(N) sfeatures.append(features[idx]) slabels.append(labels[idx]) return np.array(sfeatures), np.array(slabels) class rf_model(object): def __init__(self, forest, names): self.forest = forest self.names = names def apply(self, features): rf = len(self.forest) votes = sum(t.apply(features) for t in self.forest) return (votes > (rf//2)) class rf_learner(object): ''' Random Forest Learner learner = rf_learner(rf=101, frac=.7) Attributes ---------- rf : integer, optional Nr of trees to learn (default: 101) frac : float, optional Sample fraction ''' def __init__(self, rf=101, frac=.7): self.rf = rf self.frac = frac def train(self, features, labels, normalisedlabels=False, names=None, **kwargs): N,M = features.shape m = int(self.frac*M) n = int(self.frac*M) R = np.random tree = milk.supervised.tree.tree_learner() forest = [] if not normalisedlabels: labels,names = normaliselabels(labels) elif names is None: names = (0,1) for i in xrange(self.rf): forest.append( tree.train(*_sample(features, labels, n, R), **{'normalisedlabels' : True})) # This syntax is necessary for Python 2.5 return rf_model(forest, names) milk-0.3.10/milk/supervised/defaultlearner.py0000644000175000017500000001007511556634056022037 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division import numpy as np __all__ = [ 'defaultlearner', 'svm_simple', 'feature_selection_simple', ] def defaultlearner(mode='medium'): ''' learner = defaultlearner(mode='medium') Return the default classifier learner This is an SVM based classifier using the 1-vs-1 technique for multi-class problems. The features will be first cleaned up (normalised to [-1, +1]) and go through SDA feature selection. Parameters ----------- mode : string, optional One of ('fast','medium','slow', 'really-slow'). This defines the speed accuracy trade-off. It essentially defines how large the SVM parameter range is. Returns ------- learner : classifier learner object See Also -------- feature_selection_simple : Just perform the feature selection svm_simple : Perform classification ''' # These cannot be imported at module scope! # The reason is that they introduce a dependency loop: # gridsearch depends on nfoldcrossvalidation # nfoldcrossvalidation depends on defaultlearner # which cannot depend on gridsearch # # Importing at function level keeps all these issues at bay # from .classifier import ctransforms from .gridsearch import gridsearch from . import svm from .normalise import chkfinite, interval_normalise from .featureselection import sda_filter, featureselector, linear_independent_features from .multi import one_against_one assert mode in ('really-slow', 'slow', 'medium', 'fast'), \ "milk.supervised.defaultlearner: mode must be one of 'fast','slow','medium'." if mode == 'fast': c_range = np.arange(-2,4) sigma_range = np.arange(-2,3) elif mode == 'medium': c_range = np.arange(-2,4) sigma_range = np.arange(-4,4) elif mode == 'really-slow': c_range = np.arange(-4,10) sigma_range = np.arange(-7,7) else: # mode == 'slow' c_range = np.arange(-9,5) sigma_range = np.arange(-7,4) return ctransforms( chkfinite(), interval_normalise(), featureselector(linear_independent_features), sda_filter(), gridsearch(one_against_one(svm.svm_to_binary(svm.svm_raw())), params={ 'C': 2.**c_range, 'kernel': [svm.rbf_kernel(2.**i) for i in sigma_range], } )) def feature_selection_simple(): ''' selector = feature_selection_simple() Standard feature normalisation and selection This fills in NaNs and Infs (to 0 and large numbers), normalises features to [-1, +1] and uses SDA for feature selection. Returns ------- selector : supervised learner See Also -------- defaultlearner : perform feature selection *and* classification ''' from .classifier import ctransforms from .normalise import chkfinite, interval_normalise from .featureselection import sda_filter, featureselector, linear_independent_features return ctransforms( chkfinite(), interval_normalise(), featureselector(linear_independent_features), sda_filter(), ) def svm_simple(C, kernel): ''' learner = svm_simple(C, kernel) Returns a one-against-one SVM based classifier with `C` and `kernel` Parameters ---------- C : double C parameter kernel : kernel Kernel to use Returns ------- learner : supervised learner See Also -------- feature_selection_simple : Perform feature selection defaultlearner : feature selection and gridsearch for SVM parameters ''' from . import svm from .multi import one_against_one return one_against_one(svm.svm_to_binary(svm.svm_raw(C=C, kernel=kernel))) milk-0.3.10/milk/supervised/defaultclassifier.py0000644000175000017500000000012111556634056022522 0ustar luispedroluispedrofrom milk.supervised.defaultlearner import * defaultclassifier = defaultlearner milk-0.3.10/milk/supervised/gridsearch.py0000644000175000017500000001201011536473671021146 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division import numpy as np from .classifier import normaliselabels __all__ = [ 'gridminimise', 'gridsearch', ] def _allassignments(options): try: from itertools import product except ImportError: def product(*args, **kwds): # from http://docs.python.org/library/itertools.html#itertools.product pools = map(tuple, args) * kwds.get('repeat', 1) result = [[]] for pool in pools: result = [x+[y] for x in result for y in pool] for prod in result: yield tuple(prod) from itertools import repeat, izip for ks,vs in izip(repeat(options.keys()), product(*options.values())): yield zip(ks,vs) def _set_assignment(obj,assignments): for k,v in assignments: obj.set_option(k,v) def gridminimise(learner, features, labels, params, measure=None, nfolds=10): ''' best = gridminimise(learner, features, labels, params, measure={0/1 loss}) Grid search for the settings of parameters that maximises a given measure This function is equivalent to searching the grid, but does not actually search the whole grid. Parameters ---------- learner : a classifier object features : sequence of features labels : sequence of labels params : dictionary of sequences keys are the options to change, values are sequences of corresponding elements to try measure : function, optional a function that takes labels and outputs and returns the loss. Default: 0/1 loss. This must be an *additive* function. nfolds : integer, optional nr of folds to run, default: 10 Returns ------- best : a sequence of assignments ''' # The algorithm is as follows: # # for all assignments: error = 0, next_iteration = 0 # # at each iteration: # look for assignment with smallest error # if that is done: return it # else: perform one more iteration # # When the function returns, that assignment has the lowest error of all # assignments and all the iterations are done. Therefore, other assignments # could only be worse even if we never computed the whole error! from ..measures.nfoldcrossvalidation import foldgenerator if measure is None: def measure(real, preds): return np.sum(np.asarray(real) != np.asarray(preds)) labels,_ = normaliselabels(labels) allassignments = list(_allassignments(params)) N = len(allassignments) iteration = np.zeros(N, int) error = np.zeros(N, float) folds = [(Tr.copy(), Te.copy()) for Tr,Te in foldgenerator(labels, nfolds)] # foldgenerator might actually decide on a smaller number of folds, # depending on the distribution of class sizes: nfolds = len(folds) while True: next_pos = (error == error.min()) iter = iteration[next_pos].max() if iter == nfolds: (besti,) = np.where(next_pos & (iteration == iter)) besti = besti[0] return allassignments[besti] (ps,) = np.where(next_pos & (iteration == iter)) p = ps[0] _set_assignment(learner, allassignments[p]) train, test = folds[iter] model = learner.train(features[train], labels[train], normalisedlabels=True) preds = [model.apply(f) for f in features[test]] error[p] += measure(labels[test], preds) iteration[p] += 1 class gridsearch(object): ''' G = gridsearch(base, measure=accuracy, nfolds=10, params={ 'param1 : [...], param2 : [...]}) Perform a grid search for the best parameter values. When G.train() is called, then for each combination of p1 in param1, p2 in param2, ... it performs:: base_classifier.param1 = p1 base_classifier.param2 = p2 ... value[p1, p2,...] = measure(crossvaliation(base_classifier) it then picks the highest set of parameters and re-learns a model on the whole data. Parameters ----------- base_classifier : classifier to use measure : function, optional a function that takes labels and outputs and returns the loss. Default: 0/1 loss. This must be an *additive* function. nfolds : integer, optional Nr of folds params : dictionary ''' def __init__(self, base, measure=None, nfolds=10, params={}): self.params = params self.base = base self.nfolds = 10 self.measure = measure def is_multi_class(self): return self.base.is_multi_class() def train(self, features, labels, normalisedlabels=False): self.best = gridminimise(self.base, features, labels, self.params, self.measure, self.nfolds) _set_assignment(self.base, self.best) return self.base.train(features, labels, normalisedlabels=normalisedlabels) milk-0.3.10/milk/supervised/logistic.py0000644000175000017500000000704311524601354020646 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division import numpy as np from .normalise import normaliselabels __all__ = [ 'logistic_learner', ] @np.vectorize def _sigmoidal(z): if (z > 300): return 1. if z < -300: return 0. return 1./(1+np.exp(-z)) class logistic_model(object): def __init__(self, bs): self.bs = bs def apply(self, fs): return _sigmoidal(self.bs[0] + np.dot(fs, self.bs[1:])) class logistic_learner(object): ''' learner = logistic_learner(alpha=0.0) Logistic regression learner There are two implementations: 1. One which depends on ``scipy.optimize``. This is the default and is extremely fast. 2. If ``import scipy`` fails, then we fall back to a Python only gradient-descent. This gives good results, but is many times slower. Properties ---------- alpha : real, optional penalty for L2-normalisation. Default is zero, for no penalty. ''' def __init__(self, alpha=0.0): self.alpha = alpha def train(self, features, labels, normalisedlabels=False, names=None, **kwargs): def error(bs): response = bs[0] + np.dot(features, bs[1:]) response = _sigmoidal(response) diff = response - labels log_like = np.dot(diff, diff) L2_penalty = self.alpha * np.dot(bs, bs) return log_like + L2_penalty def error_prime(bs): fB = np.dot(features, bs[1:]) response = _sigmoidal(bs[0] + fB) sprime = response * (1-response) ds = (response - labels) * sprime b0p = np.sum(ds) b1p = np.dot(features.T, ds) bp = np.concatenate( ([b0p], b1p) ) return 2.*(bp + self.alpha*bs) features = np.asanyarray(features) if not normalisedlabels: labels, _ = normaliselabels(labels) N,f = features.shape bs = np.zeros(f+1) try: from scipy import optimize # Some testing revealed that this was a good combination # call fmin_cg twice first and then fmin # I do not understand why 100%, but there it is bs = optimize.fmin_cg(error, bs, error_prime, disp=False) bs = optimize.fmin_cg(error, bs, error_prime, disp=False) bs = optimize.fmin(error, bs, disp=False) except ImportError: import warnings warnings.warn('''\ milk.supervised.logistic.train: Could not import scipy.optimize. Fall back to very simple gradient descent (which is slow).''') bs = np.zeros(f+1) cur = 1.e-6 ebs = error(bs) for i in xrange(1000000): dir = error_prime(bs) step = (lambda e : bs - e *dir) enbs = ebs + 1 while enbs > ebs: cur /= 2. if cur == 0.: break nbs = step(cur) enbs = error(nbs) while cur < 10.: cur *= 2 nnbs = step(cur) ennbs = error(nnbs) if ennbs < enbs: nbs = nnbs enbs = ennbs else: break bs = nbs ebs = enbs return logistic_model(bs) milk-0.3.10/milk/supervised/__init__.py0000644000175000017500000000333411511426616020571 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT. See COPYING.MIT file in the milk distribution ''' milk.supervised This hold the supervised classification modules: Submodules ---------- - defaultclassifier: contains a default "good enough" classifier - svm: related to SVMs - grouped: contains objects to transform single object classifiers into group classifiers by voting - multi: transforms binary classifiers into multi-class classifiers (1-vs-1 or 1-vs-rest) - featureselection: feature selection - knn: k-nearest neighbours - tree: decision tree classifiers Classifiers ----------- All classifiers have a `train` function which takes 2 arguments: - features : sequence of features - labels : sequence of labels They return a `model` object, which has an `apply` function which takes a single input and returns its label. Note that there are always two objects: the learned and the model and they are independent. Every time you call learner.train() you get a new model. Both classifiers and models are pickle()able. Example ------- :: features = np.random.randn(100,20) features[:50] *= 2 labels = np.repeat((0,1), 50) classifier = milk.defaultclassifier() model = classifier.train(features, labels) new_label = model.apply(np.random.randn(100)) new_label2 = model.apply(np.random.randn(100)*2) ''' from .defaultclassifier import defaultclassifier, svm_simple from .classifier import normaliselabels from .gridsearch import gridsearch from .tree import tree_learner __all__ = [ 'normaliselabels', 'defaultclassifier', 'svm_simple', 'gridsearch', 'tree_learner', ] milk-0.3.10/milk/supervised/multi.py0000644000175000017500000001236711531303445020166 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division from .classifier import normaliselabels import numpy as np __all__ = [ 'one_against_rest', 'one_against_one', 'one_against_rest_multi', ] def _asanyarray(f): try: return np.asanyarray(f) except: return np.array(f, dtype=object) class one_against_rest(object): ''' Implements one vs. rest classification strategy to transform a binary classifier into a multi-class classifier. classifier = one_against_rest(base) base must obey the classifier interface Example ------- :: multi = one_against_rest(milk.supervised.simple_svm()) model = multi.train(training_features,labels) print model.apply(testing_features) See Also -------- one_against_one ''' def __init__(self,base): self.base = base self.is_multi_class = True self.options = {} def set_option(self, k, v): self.options[k] = v def train(self, features, labels, normalisedlabels=False): labels, names = normaliselabels(labels) nclasses = labels.max() + 1 models = [] for i in xrange(nclasses): for k,v in self.options.iteritems(): self.base.set_option(k, v) model = self.base.train(features, (labels == i).astype(int), normalisedlabels=True) models.append(model) return one_against_rest_model(models, names) class one_against_rest_model(object): def __init__(self, models, names): self.models = models self.nclasses = len(self.models) self.names = names def apply(self, feats): vals = np.array([c.apply(feats) for c in self.models]) (idxs,) = np.where(vals) if len(idxs) == 1: (label,) = idxs elif len(idxs) == 0: label = 0 else: label = idxs[0] return self.names[label] class one_against_one(object): ''' Implements one vs. one classification strategy to transform a binary classifier into a multi-class classifier. classifier = one_against_one(base) base must obey the classifier interface Example ------- :: multi = one_against_one(milk.supervised.simple_svm()) multi.train(training_features,labels) print multi.apply(testing_features) See Also -------- one_against_rest ''' def __init__(self, base): self.base = base self.is_multi_class = True self.options = {} def set_option(self, k, v): self.options[k] = v def train(self, features, labels, **kwargs): ''' one_against_one.train(objs,labels) ''' labels, names = normaliselabels(labels) features = _asanyarray(features) nclasses = labels.max() + 1 models = [ [None for i in xrange(nclasses)] for j in xrange(nclasses)] for i in xrange(nclasses): for j in xrange(i+1, nclasses): for k,v in self.options.iteritems(): self.base.set_option(k, v) idxs = (labels == i) | (labels == j) assert idxs.sum() > 0, 'milk.multi.one_against_one: Pair-wise classifier has no data' # Fixme: here I could add a Null model or something model = self.base.train(features[idxs], (labels[idxs]==i).astype(int), normalisedlabels=True) models[i][j] = model return one_against_one_model(models, names) class one_against_one_model(object): def __init__(self, models, names): self.models = models self.names = names self.nclasses = len(models) def apply(self,feats): ''' one_against_one.apply(objs) Classify one single object. ''' nc = self.nclasses votes = np.zeros(nc) for i in xrange(nc): for j in xrange(i+1,nc): c = self.models[i][j].apply(feats) if c: votes[i] += 1 else: votes[j] += 1 return self.names[votes.argmax(0)] class one_against_rest_multi_model(object): def __init__(self, models): self.models = models def apply(self, feats): return [lab for lab,model in self.models.iteritems() if model.apply(feats)] class one_against_rest_multi(object): ''' learner = one_against_rest_multi() model = learner.train(features, labels) classes = model.apply(f_test) This for multi-label problem (i.e., each instance can have more than one label). ''' def __init__(self, base): self.base = base if hasattr(base, 'set_option'): self.set_option = base.set_option def train(self, features, labels, normalisedlabels=False): ''' ''' import operator all_labels = set() for ls in labels: all_labels.update(ls) models = {} for label in all_labels: models[label] = self.base.train(features, [(label in ls) for ls in labels]) return one_against_rest_multi_model(models) milk-0.3.10/milk/supervised/set2binary_array.py0000644000175000017500000000157411511426617022317 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division import numpy as np __all__ = [ 'set2binary_array', ] class set2binary_array_model(object): def __init__(self, universe): self.universe = list(universe) def apply(self, elems): res = np.zeros(len(self.universe) + 1, bool) for e in elems: try: res[self.universe.index(e)] = True except : res[-1] = True return res class set2binary_array(object): def train(self, features, labels, normalisedlabels=False): allfeatures = set() for f in features: allfeatures.update(f) return set2binary_array_model(allfeatures) milk-0.3.10/milk/supervised/multi_view.py0000644000175000017500000000363611511426617021224 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT. See COPYING.MIT file in the milk distribution import numpy as np __all__ = [ 'multi_view_learner', ] class multi_view_model(object): def __init__(self, models): self.models = models def apply(self, features): if len(features) != len(self.models): raise ValueError('milk.supervised.two_view: Nr of features does not match training data (got %s, expected %s)' % (len(features) ,len(self.models))) Ps = np.array([model.apply(f) for model,f in zip(self.models, features)]) if np.any(Ps <= 0.): return False if np.any(Ps >= 1.): return True # This is binary only: # if \prod Pi > \prod (1-Pi) return 1 # is equivalent to # if \prod Pi/(1-Pi) > 1. return 1 # if \sum \log( Pi/(1-Pi) ) > 0. return 1 return np.sum( np.log(Ps/(1-Ps)) ) > 0 class multi_view_learner(object): ''' Multi View Learner This learner learns different classifiers on multiple sets of features and combines them for classification. ''' def __init__(self, bases): self.bases = bases def train(self, features, labels, normalisedlabels=False): features = zip(*features) if len(features) != len(self.bases): raise ValueError('milk.supervised.multi_view_learner: ' + 'Nr of features does not match classifiser construction (got %s, expected %s)' % (len(features) ,len(self.bases))) models = [] for basis,f in zip(self.bases, features): try: f = np.array(f) except: f = np.array(f, dtype=object) models.append(basis.train(f, labels)) return multi_view_model(models) multi_view_classifier = multi_view_learner milk-0.3.10/milk/supervised/grouped.py0000644000175000017500000000761411511426617020505 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2010-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT. See COPYING.MIT file in the milk distribution # -*- coding: utf-8 -*- from __future__ import division import numpy as np from collections import defaultdict from .classifier import normaliselabels __all__ = [ 'voting_learner', 'mean_learner', 'remove_outliers', 'filter_outliers', ] def _concatenate_features_labels(gfeatures, glabels): if type(gfeatures) == np.ndarray and gfeatures.dtype == object: gfeatures = list(gfeatures) features = np.concatenate(gfeatures) labels = [] for feats,label in zip(gfeatures, glabels): labels.extend( [label] * len(feats) ) return features, labels class voting_learner(object): ''' Implements a voting scheme for multiple sub-examples per example. classifier = voting_learner(base) base should be a binary classifier Example ------- :: voterlearn = voting_learner(milk.supervised.simple_svm()) voter = voterlearn.train(training_groups, labeled_groups) res = voter.apply([ [f0, f1, f3] ]) ''' def __init__(self, base): self.base = base def train(self, gfeatures, glabels, normalisedlabels=False): features, labels = _concatenate_features_labels(gfeatures, glabels) return voting_model(self.base.train(features, labels)) voting_classifier = voting_learner class voting_model(object): def __init__(self, base): self.base = base def apply(self, gfeatures): votes = defaultdict(int) for feats in gfeatures: votes[self.base.apply(feats)] += 1 best = None most_votes = 0 for k,v in votes.iteritems(): if v > most_votes: best = k most_votes = v return best class mean_learner(object): ''' Implements a mean scheme for multiple sub-examples per example. classifier = mean_learner(base) `base` should be a classifier that returns a numeric confidence value `classifier` will return the **mean** Example ------- :: meanlearner = mean_learner(milk.supervised.raw_svm()) model = meanlearner.train(training_groups, labeled_groups) res = model.apply([ [f0, f1, f3] ]) ''' def __init__(self, base): self.base = base if hasattr(base, 'set_option'): self.set_option = base.set_option def train(self, gfeatures, glabels, normalisedlabels=False): features, labels = _concatenate_features_labels(gfeatures, glabels) return mean_model(self.base.train(features, labels)) mean_classifier = mean_learner class mean_model(object): def __init__(self, base): self.base = base def apply(self, gfeatures): return np.mean([self.base.apply(feats) for feats in gfeatures]) def remove_outliers(features, limit, min_size): ''' features = remove_outliers(features, limit, min_size) ''' nsize = int(limit * len(features)) if nsize < min_size: return features normed = features - features.mean(0) std = normed.std(0) std[std == 0] = 1 normed /= std f2_sum1 = (normed**2).mean(1) values = f2_sum1.copy() values.sort() top = values[nsize] selected = f2_sum1 < top return features[selected] class filter_outliers_model(object): def __init__(self, limit, min_size): self.limit = limit self.min_size = min_size def apply(self, features): return remove_outliers(features, self.limit, self.min_size) class filter_outliers(object): def __init__(self, limit=.9, min_size=3): self.limit = limit self.min_size = min_size def train(self, features, labels, normalisedlabels=False): return filter_outliers_model(self.limit, self.min_size) milk-0.3.10/milk/supervised/tree.py0000644000175000017500000001624311537456345020006 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT. See COPYING.MIT file in the milk distribution ''' ================ Tree Classifier ================ Decision tree based classifier ------------------------------ ''' from __future__ import division import numpy as np from .classifier import normaliselabels __all__ = [ 'tree_learner', 'stump_learner', ] class Leaf(object): ''' v : value w : weight ''' def __init__(self, v, w): self.v = v self.w = w def __repr__(self): return 'Leaf(%s,%s)' % (self.v, self.w) class Node(object): # This could be replaced by a namedtuple def __init__(self, featid, featval, left, right): self.featid = featid self.featval = featval self.left = left self.right = right def _split(features, labels, weights, criterion, subsample, R): N,f = features.shape if subsample is not None: samples = np.array(R.sample(xrange(features.shape[1]), subsample)) features = features[:, samples] f = subsample best = None best_val = -1. for i in xrange(f): domain_i = sorted(set(features[:,i])) for d in domain_i[1:]: cur_split = (features[:,i] < d) if weights is not None: value = criterion(labels[cur_split], labels[~cur_split], weights[cur_split], weights[~cur_split]) else: value = criterion(labels[cur_split], labels[~cur_split]) if value > best_val: best_val = value if subsample is not None: ti = samples[i] else: ti = i best = ti,d return best from ._tree import set_entropy from ._tree import information_gain as _information_gain def information_gain(labels0, labels1, include_entropy=False): ''' ig = information_gain(labels0, labels1, include_entropy=False) Information Gain See http://en.wikipedia.org/wiki/Information_gain_in_decision_trees The function calculated here does not include the original entropy unless you explicitly ask for it (by passing include_entropy=True) ''' if include_entropy: return set_entropy(np.concatenate( (labels0, labels1) )) + \ _information_gain(labels0, labels1) return _information_gain(labels0, labels1) def z1_loss(labels0, labels1, weights0=None, weights1=None): ''' z = z1_loss(labels0, labels1) z = z1_loss(labels0, labels1, weights0, weights1) zero-one loss split for tree learning ''' def _acc(labels, weights): c = (labels.mean() > .5) if weights is not None: return -np.dot((labels != c), weights) return -np.sum(labels != c) return _acc(labels0, weights0) + _acc(labels1, weights1) def build_tree(features, labels, criterion, min_split=4, subsample=None, R=None, weights=None): ''' tree = build_tree(features, labels, criterion, min_split=4, subsample=None, R=None, weights={all 1s}) Parameters ---------- features : sequence features to use labels : sequence labels criterion : function {labels} x {labels} -> float function to measure goodness of split min_split : integer minimum size to split on subsample : integer, optional if given, then, at each step, choose R : source of randomness, optional See `milk.util.get_pyrandom` weights : sequence, optional weight of instance (default: all the same) Returns ------- tree : Tree ''' assert len(features) == len(labels), 'build_tree: Nr of labels does not match nr of features' features = np.asanyarray(features) labels = np.asanyarray(labels, dtype=np.int) if subsample is not None: if subsample <= 0: raise ValueError('milk.supervised.tree.build_tree: `subsample` must be > 0.\nDid you mean to use None to signal no subsample?') from ..utils import get_pyrandom R = get_pyrandom(R) def recursive(features, labels): N = float(len(labels)) if N < min_split: return Leaf(labels.sum()/N, N) S = _split(features, labels, weights, criterion, subsample, R) if S is None: return Leaf(labels.sum()/N, N) i,thresh = S split = features[:,i] < thresh return Node(featid=i, featval=thresh, left =recursive(features[ split], labels[ split]), right=recursive(features[~split], labels[~split])) return recursive(features, labels) def apply_tree(tree, features): ''' conf = apply_tree(tree, features) Applies the decision tree to a set of features. ''' if type(tree) is Leaf: return tree.v if features[tree.featid] < tree.featval: return apply_tree(tree.left, features) return apply_tree(tree.right, features) class tree_learner(object): ''' tree = tree_learner() model = tree.train(features, labels) model2 = tree.train(features, labels, weights=weights) predicted = model.apply(testfeatures) A decision tree classifier (currently, implements the greedy ID3 algorithm without any pruning). Attributes ---------- criterion : function, optional criterion to use for tree construction, this should be a function that receives a set of labels (default: information_gain). min_split : integer, optional minimum size to split on (default: 4). ''' def __init__(self, criterion=information_gain, min_split=4, return_label=True, subsample=None, R=None): self.criterion = criterion self.min_split = min_split self.return_label = return_label self.subsample = subsample self.R = R def train(self, features, labels, normalisedlabels=False, weights=None): if not normalisedlabels: labels,names = normaliselabels(labels) tree = build_tree(features, labels, self.criterion, self.min_split, self.subsample, self.R, weights) return tree_model(tree, self.return_label) tree_classifier = tree_learner class tree_model(object): ''' tree model ''' def __init__(self, tree, return_label): self.tree = tree self.return_label = return_label def apply(self,feats): value = apply_tree(self.tree, feats) if self.return_label: return value > .5 return value class stump_model(object): def __init__(self, idx, cutoff, names): self.names = names self.idx = idx self.cutoff = cutoff def apply(self, f): return self.names[f[self.idx] > self.cutoff] def __repr__(self): return '' % (self.idx, self.cutoff) class stump_learner(object): def __init__(self): pass def train(self, features, labels, normalisedlabels=False, weights=None, **kwargs): if not normalisedlabels: labels,names = normaliselabels(labels) idx,cutoff = _split(features, labels, weights, z1_loss, subsample=None, R=None) return stump_model(idx, cutoff, names) milk-0.3.10/milk/supervised/featureselection.py0000644000175000017500000001562411556634056022410 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division import numpy as np from numpy.linalg import det from . classifier import normaliselabels __all__ = [ 'sda', 'linearly_independent_subset', 'linear_independent_features', 'filterfeatures', 'featureselector', 'sda_filter', ] _TOLERANCE = 0 _SIGNIFICANCE_IN = .15 _SIGNIFICANCE_OUT = .15 def _sweep(A, k, flag): Akk = A[k,k] # cross[i,j] = A[i,k] * A[k,j] cross = (A[:,k][:, np.newaxis] * A[k]) B = A - cross/Akk # currently: B[i,j] = A[i,j] - A[i,k]*A[k,j]/Akk # Now fix row k and col k, followed by Bkk B[k] = flag * A[k]/A[k,k] B[:,k] = flag * A[:,k]/A[k,k] B[k,k] = -1./Akk return B def sda(features, labels, tolerance=None, significance_in=None, significance_out=None): ''' features_idx = sda(features,labels) Perform Stepwise Discriminant Analysis for feature selection Pre filter the feature matrix to remove linearly dependent features before calling this function. Behaviour is undefined otherwise. This implements the algorithm described in Jennrich, R.I. (1977), "Stepwise Regression" & "Stepwise Discriminant Analysis," both in Statistical Methods for Digital Computers, eds. K. Enslein, A. Ralston, and H. Wilf, New York; John Wiley & Sons, Inc. ''' import scipy.stats if tolerance is None: tolerance = _TOLERANCE if significance_in is None: significance_in = _SIGNIFICANCE_IN if significance_out is None: significance_out = _SIGNIFICANCE_OUT assert len(features) == len(labels), 'milk.supervised.featureselection.sda: length of features not the same as length of labels' N, m = features.shape labels,labelsu = normaliselabels(labels) q = len(labelsu) df = features - features.mean(0) T = np.dot(df.T, df) dfs = [(features[labels == i] - features[labels == i].mean(0)) for i in xrange(q)] W = np.sum(np.dot(d.T, d) for d in dfs) ignoreidx = ( W.diagonal() == 0 ) if ignoreidx.any(): idxs, = np.where(~ignoreidx) F = sda(features[:,~ignoreidx],labels) if len(F): return idxs[F] return ignoreidx output = [] D = W.diagonal() df1 = q-1 last_enter_k = -1 while True: V = W.diagonal()/T.diagonal() W_d = W.diagonal() V_neg = (W_d < 0) p = V_neg.sum() if V_neg.any(): V_m = V[V_neg].min() k, = np.where(V == V_m) k = k[0] Fremove = (N-p-q+1)/(q-1)*(V_m-1) df2 = N-p-q+1 PrF = 1 - scipy.stats.f.cdf(Fremove,df1,df2) if PrF > significance_out: #print 'removing ',k, 'V(k)', 1./V_m, 'Fremove', Fremove, 'df1', df1, 'df2', df2, 'PrF', PrF if k == last_enter_k: # We are going into an infinite loop. import warnings warnings.warn('milk.featureselection.sda: infinite loop detected (maybe bug?).') break W = _sweep(W,k,1) T = _sweep(T,k,1) continue ks = ( (W_d / D) > tolerance) if ks.any(): V_m = V[ks].min() k, = np.where(V==V_m) k = k[0] Fenter = (N-p-q)/(q-1) * (1-V_m)/V_m df2 = N-p-q PrF = 1 - scipy.stats.f.cdf(Fenter,df1,df2) if PrF < significance_in: #print 'adding ',k, 'V(k)', 1./V_m, 'Fenter', Fenter, 'df1', df1, 'df2', df2, 'PrF', PrF W = _sweep(W,k,-1) T = _sweep(T,k,-1) if PrF < .0001: output.append((Fenter,k)) last_enter_k = k continue break output.sort(reverse=True) return np.array([x[1] for x in output]) def linearly_independent_subset(V, threshold=1.e-5, return_orthogonal_basis=False): ''' subset = linearly_independent_subset(V, threshold=1.e-5) subset,U = linearly_independent_subset(V, threshold=1.e-5, return_orthogonal_basis=True) Discover a linearly independent subset of `V` Parameters ---------- V : sequence of input vectors threshold : float, optional vectors with 2-norm smaller or equal to this are considered zero (default: 1e.-5) return_orthogonal_basis : Boolean, optional whether to return orthogonal basis set Returns ------- subset : ndarray of integers indices used for basis U : 2-array orthogonal basis into span{V} Implementation Reference ------------------------ Use Gram-Schmidt with a check for when the v_k is close enough to zero to ignore See http://en.wikipedia.org/wiki/Gram-Schmidt_process ''' V = np.array(V, copy=True) orthogonal = [] used = [] for i,u in enumerate(V): for v in orthogonal: u -= np.dot(u,v)/np.dot(v,v) * v if np.dot(u,u) > threshold: orthogonal.append(u) used.append(i) if return_orthogonal_basis: return np.array(used),np.array(orthogonal) return np.array(used) def linear_independent_features(features, labels=None): ''' indices = linear_independent_features(features, labels=None) Returns the indices of a set of linearly independent features (columns). Parameters ---------- features : ndarray labels : ignored This argument is only here to conform to the learner interface. Returns ------- indices : ndarray of integers indices of features to keep See Also -------- `linearly_independent_subset` : this function is equivalent to `linearly_independent_subset(features.T)` ''' return linearly_independent_subset(features.T) class filterfeatures(object): ''' selector = filterfeatures(idxs) Returns a transformer which selects the features given by idxs ''' def __init__(self, idxs): self.idxs = idxs def apply(self, features): return features[self.idxs] def __repr__(self): return 'filterfeatures(%s)' % self.idxs class featureselector(object): ''' selector = featureselector(function) Returns a transformer which selects features according to selected_idxs = function(features,labels) ''' def __init__(self, selector): self.selector = selector def train(self, features, labels, normalisedlabels=False): idxs = self.selector(features, labels) if len(idxs) == 0: import warnings warnings.warn('milk.featureselection: No features selected! Using all features as fall-back.') idxs = np.arange(len(features[0])) return filterfeatures(idxs) def __repr__(self): return 'featureselector(%s)' % self.selector def sda_filter(): return featureselector(sda) milk-0.3.10/milk/supervised/adaboost.py0000644000175000017500000000411411511426616020623 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division import numpy as np from .normalise import normaliselabels ''' AdaBoost Simple implementation of Adaboost Learner ------- boost_learner ''' __all__ = [ 'boost_learner', ] def _adaboost(features, labels, base, max_iters): m = len(features) D = np.ones(m, dtype=float) D /= m Y = np.ones(len(labels), dtype=float) names = np.array([-1, +1]) Y = names[labels] H = [] A = [] for t in xrange(max_iters): Ht = base.train(features, labels, weights=D) train_out = np.array(map(Ht.apply, features)) train_out = names[train_out] Et = np.dot(D, (Y != train_out)) if Et > .5: # early return break At = .5 * np.log((1. + Et) / (1. - Et)) D *= np.exp((-At) * Y * train_out) D /= np.sum(D) A.append(At) H.append(Ht) return H, A class boost_model(object): def __init__(self, H, A, names): self.H = H self.A = A self.names = names def apply(self, f): v = sum((a*h.apply(f)) for h,a in zip(self.H, self.A)) v /= np.sum(self.A) return self.names[v > .5] class boost_learner(object): ''' learner = boost_learner(weak_learner_type(), max_iters=100) model = learner.train(features, labels) test = model.apply(f) AdaBoost learner Attributes ---------- base : learner Weak learner max_iters : integer Nr of iterations (default: 100) ''' def __init__(self, base, max_iters=100): self.base = base self.max_iters = max_iters def train(self, features, labels, normalisedlabels=False, names=(0,1), weights=None, **kwargs): if not normalisedlabels: labels,names = normaliselabels(labels) H,A = _adaboost(features, labels, self.base, self.max_iters) return boost_model(H, A, names) milk-0.3.10/milk/supervised/classifier.py0000644000175000017500000000650711521632340021155 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. from __future__ import division import numpy as np __all__ = ['normaliselabels', 'ctransforms'] def normaliselabels(labels): ''' normalised, names = normaliselabels(labels) Normalises the labels to be integers from 0 through N-1 `normalised` is a np.array, while `names` is a list mapping the indices to the old names. Parameters ---------- labels : any iterable of labels Returns ------ normalised : a numpy ndarray of integers 0 .. N-1 names : list of label names ''' names = sorted(set(labels)) normalised = map(names.index, labels) return np.array(normalised), names class threshold_model(object): ''' threshold_model Attributes ---------- threshold : float threshold value ''' def __init__(self, threshold=.5): self.threshold = .5 def apply(self, f): return f >= self.threshold class fixed_threshold_learner(object): def __init__(self, threshold=.5): self.threshold = threshold def train(self, features, labels, **kwargs): return threshold_model(self.threshold) class ctransforms_model(object): ''' model = ctransforms_model(models) A model that consists of a series of transformations. See Also -------- ctransforms ''' def __init__(self, models): self.models = models def apply(self,features): for T in self.models: features = T.apply(features) return features class ctransforms(object): ''' ctransf = ctransforms(c0, c1, c2, ...) Concatenate transforms. ''' def __init__(self,*args): self.transforms = args def train(self, features, labels, normalisedlabels=False): models = [] model = None for T in self.transforms: if model is not None: features = np.array([model.apply(f) for f in features]) model = T.train(features, labels, normalisedlabels=normalisedlabels) models.append(model) return ctransforms_model(models) def set_option(self, opt, val): idx, opt = opt self.transforms[idx].set_option(opt,val) milk-0.3.10/milk/supervised/_svm.cpp0000644000175000017500000010234111453701533020125 0ustar luispedroluispedro// Copyright (C) 2008, Luís Pedro Coelho // Copyright (c) 2000-2008 Chih-Chung Chang and Chih-Jen Lin (LIBSVM Code) // // Permission is hereby granted, free of charge, to any person obtaining a copy // of this software and associated documentation files (the "Software"), to deal // in the Software without restriction, including without limitation the rights // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell // copies of the Software, and to permit persons to whom the Software is // furnished to do so, subject to the following conditions: // // The above copyright notice and this permission notice shall be included in // all copies or substantial portions of the Software. // // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN // THE SOFTWARE. #include #include #include #include #include #include #include #include //#include //#include //using __gnu_debug::vector; //using __gnu_debug::list; using std::vector; using std::list; extern "C" { #include #include } namespace { const double INF = HUGE_VAL; // This is a boost function // Copied here for convenience. template inline Iter prior(Iter it) { --it; return it; } using std::max; using std::min; template T median(T bot, T mid, T hi) { if (mid < bot) return bot; if (mid > hi) return hi; return mid; } struct SMO_Exception { SMO_Exception(const char* msg): msg(msg) { } const char* msg; }; struct Python_Exception { }; void check_for_interrupts() { PyErr_CheckSignals(); if (PyErr_Occurred()) { throw Python_Exception(); } } class KernelComputation { public: virtual ~KernelComputation() { } virtual double do_kernel(int i1, int i2) const = 0; }; class KernelCache { public: KernelCache(std::auto_ptr computation_, int N, int cache_nr_doubles); virtual ~KernelCache(); double* get_kline(int idx, int size = -1); double* get_diag(); double kernel_apply(int i1, int i2) const; protected: virtual double do_kernel(int i1, int i2) const { return computation_->do_kernel(i1,i2); } const int N_; private: std::auto_ptr computation_; double ** cache_; double * dcache_; int cache_free_; list cache_lru; vector::iterator> cache_iter; }; KernelCache::KernelCache(std::auto_ptr computation, int N, int cache_nr_floats): N_(N), computation_(computation), dcache_(0) { cache_ = new double*[N_]; for (int i = 0; i != N_; ++i) cache_[i] = 0; cache_free_ = (cache_nr_floats/N_); cache_iter.resize(N_,cache_lru.end()); } KernelCache::~KernelCache() { for (int i = 0; i != N_; ++i) delete [] cache_[i]; delete [] cache_; delete [] dcache_; } double* KernelCache::get_kline(int idx, int s) { if (s == -1) s = N_; assert (s <= N_); if (!cache_[idx]) { if (!cache_free_) { int to_remove = cache_lru.front(); cache_lru.pop_front(); cache_[idx] = cache_[to_remove]; cache_[to_remove] = 0; } else { cache_[idx] = new double[N_]; --cache_free_; } for (int i = 0; i != N_; ++i) { if (i == idx && dcache_) cache_[i][i] = dcache_[i]; else if(i != idx && cache_[i]) cache_[idx][i] = cache_[i][idx]; else cache_[idx][i] = do_kernel(idx,i); } } else { cache_lru.erase(cache_iter[idx]); } cache_lru.push_back(idx); cache_iter[idx]=prior(cache_lru.end()); return cache_[idx]; } double* KernelCache::get_diag() { if (!dcache_) { dcache_ = new double[N_]; for (int i = 0; i != N_; ++i) { if (cache_[i]) dcache_[i] = cache_[i][i]; else dcache_[i] = do_kernel(i,i); } } return dcache_; } /*** * Returns the value of Kernel(X_i1, X_i2). * Uses the cache if possible, but does not update it. */ double KernelCache::kernel_apply(int i1, int i2) const { if (cache_[i1]) { assert(do_kernel(i1,i2) == cache_[i1][i2]); return cache_[i1][i2]; } if (cache_[i2]) { assert(do_kernel(i1,i2) == cache_[i2][i1]); return cache_[i2][i1]; } if (i1 == i2 && dcache_) { assert(do_kernel(i1,i2) == dcache_[i1]); return dcache_[i1]; } return do_kernel(i1,i2); } class PyKernel : public KernelComputation { PyKernel(const PyKernel&); public: PyKernel(PyObject* X, PyObject* kernel, int N); ~PyKernel(); protected: virtual double do_kernel(int i1, int i2) const; private: PyObject* const X_; PyObject* const pykernel_; const int N_; }; PyKernel::PyKernel(PyObject* X, PyObject* kernel, int N): X_(X), pykernel_(kernel), N_(N) { Py_INCREF(X); Py_INCREF(kernel); } PyKernel::~PyKernel() { Py_DECREF(X_); Py_DECREF(pykernel_); } double PyKernel::do_kernel(int i1, int i2) const { assert(i1 < N_); assert(i2 < N_); PyObject* obj1 = PySequence_GetItem(X_,i1); PyObject* obj2 = PySequence_GetItem(X_,i2); if (!obj1 || !obj2) { Py_XDECREF(obj1); Py_XDECREF(obj2); throw SMO_Exception("svm.eval_SMO: Unable to access element in X array"); } PyObject* arglist = Py_BuildValue("(OO)",obj1,obj2); PyObject* result = PyEval_CallObject(pykernel_,arglist); Py_XDECREF(obj1); Py_XDECREF(obj2); Py_DECREF(arglist); if (!result) { check_for_interrupts(); throw SMO_Exception("svm.eval_SMO: Unable to call kernel"); } double val = PyFloat_AsDouble(result); Py_DECREF(result); return val; } class RBFKernel : public KernelComputation { public: RBFKernel(PyArrayObject* X, double gamma); ~RBFKernel(); protected: virtual double do_kernel(int i1, int i2) const; private: PyArrayObject* X_; double ngamma_; const int N_; const int N1_; }; RBFKernel::RBFKernel(PyArrayObject* X, double gamma): X_(reinterpret_cast(X)), ngamma_(-1./gamma), N_(PyArray_DIM(X,0)), N1_(PyArray_DIM(X,1)) { if (!PyArray_Check(X)) { throw SMO_Exception("RBF Kernel used, but not with numpy array."); } if (!PyArray_ISCARRAY(X)) { throw SMO_Exception("RBF Kernel used but not with CARRAY."); } Py_INCREF(X); } RBFKernel::~RBFKernel() { Py_DECREF(X_); } double RBFKernel::do_kernel(int i1, int i2) const { assert(i1 < N_); assert(i2 < N_); const double* data1 = static_cast(PyArray_GETPTR1(X_,i1)); const double* data2 = static_cast(PyArray_GETPTR1(X_,i2)); double sumdiff = 0.; for (int i = 0; i != N1_; ++i) { double diff = data1[i]-data2[i]; sumdiff += diff * diff; } sumdiff *= ngamma_; double res = std::exp(sumdiff); return res; } class DotKernel : public KernelComputation { public: explicit DotKernel(PyArrayObject* X); ~DotKernel(); protected: virtual double do_kernel(int i1, int i2) const; private: PyArrayObject* X_; const int N_; const int N1_; }; DotKernel::DotKernel(PyArrayObject* X): X_(reinterpret_cast(X)), N_(PyArray_DIM(X,0)), N1_(PyArray_DIM(X,1)) { if (!PyArray_Check(X)) { throw SMO_Exception("Dot Kernel used, but not with numpy array."); } if (!PyArray_ISCARRAY(X)) { throw SMO_Exception("Dot Kernel used but not with CARRAY."); } Py_INCREF(X); } DotKernel::~DotKernel() { Py_DECREF(X_); } double DotKernel::do_kernel(int i1, int i2) const { assert(i1 < N_); assert(i2 < N_); const double* data1 = static_cast(PyArray_GETPTR1(X_,i1)); const double* data2 = static_cast(PyArray_GETPTR1(X_,i2)); double dotsum = 0.; for (int i = 0; i != N1_; ++i) { dotsum += data1[i] * data2[i]; } return dotsum; } class PrecomputedKernel : public KernelComputation { public: PrecomputedKernel(PyArrayObject* X); ~PrecomputedKernel(); protected: virtual double do_kernel(int i1, int i2) const; private: PyArrayObject* X_; }; PrecomputedKernel::PrecomputedKernel(PyArrayObject* X): X_(X) { if (!PyArray_Check(X)) { throw SMO_Exception("PrecomputedKernel used, but not with numpy array."); } if (!PyArray_ISCARRAY(X)) { throw SMO_Exception("Precomputed used but not with CARRAY."); } Py_INCREF(X); } PrecomputedKernel::~PrecomputedKernel() { Py_DECREF(X_); } double PrecomputedKernel::do_kernel(int i1, int i2) const { const double* data = static_cast(PyArray_GETPTR2(X_,i1,i2)); return *data; } std::auto_ptr get_kernel(PyObject* X, PyObject* kernel) { typedef std::auto_ptr res_type; if (PyCallable_Check(kernel)) return res_type(new PyKernel(X, kernel, PySequence_Length(X))); if (!PyTuple_Check(kernel) || PyTuple_Size(kernel) != 2) throw SMO_Exception("Cannot parse kernel."); PyObject* type = PyTuple_GET_ITEM(kernel,0); PyObject* arg = PyTuple_GET_ITEM(kernel,1); if (!PyInt_Check(type) || !PyFloat_Check(arg)) throw SMO_Exception("Cannot parse kernel (wrong types)"); long type_nr = PyInt_AsLong(type); double arg_value = PyFloat_AsDouble(arg); switch (type_nr) { case 0: return res_type(new RBFKernel(reinterpret_cast(X), arg_value)); case 1: return res_type(new PrecomputedKernel(reinterpret_cast(X))); case 2: return res_type(new DotKernel(reinterpret_cast(X))); default: throw SMO_Exception("Unknown kernel type!"); } } class LIBSVM_KernelCache : public KernelCache { public: LIBSVM_KernelCache(const int* Y, std::auto_ptr kernel, int N, int cache_nr_doubles) :KernelCache(kernel,N,cache_nr_doubles), Y_(Y) { } private: double do_kernel(int i, int j) const { double res = KernelCache::do_kernel(i,j); return res * Y_[i] * Y_[j]; } const int* Y_; }; class SMO { public: SMO(PyObject* X, int* Y, double* Alphas, double& b, double C, int N, PyObject* kernel, double eps, double tol, int cache_size): Alphas(Alphas), Y(Y), b(b), C(C), N(N), cache_(get_kernel(X,kernel),N,cache_size), eps(eps), tol(tol) { } void optimise(); double apply(int) const; bool take_step(int,int); bool examine_example(int); double get_error(int) const; private: double* Alphas; int* Y; double& b; double C; int N; mutable KernelCache cache_; const double eps; const double tol; }; double SMO::get_error(int j) const { return apply(j) - Y[j]; } double SMO::apply(int j) const { double sum = -b; double* Kernel_Line = cache_.get_kline(j); for (int i = 0; i != N; ++i) { if (Alphas[i] != C) { sum += Y[i] * Alphas[i] * Kernel_Line[i]; } } //std::cout << "SMO::apply( " << j << " ): " << sum << '\n'; return sum; } bool SMO::take_step(int i1, int i2) { //std::cout << "take_step( " << i1 << ", " << i2 << " );\n"; if (i1 == i2) return false; const double alpha1 = Alphas[i1]; const double alpha2 = Alphas[i2]; const int y1 = Y[i1]; const int y2 = Y[i2]; double L, H; if (y1 != y2) { L = max(0.,alpha2-alpha1); H = min(C,C+alpha2-alpha1); } else { L = max(0.,alpha1+alpha2-C); H = min(C,alpha1+alpha2); } if (L == H) return false; const int s = y1*y2; const double E1 = get_error(i1); const double E2 = get_error(i2); const double k11 = cache_.kernel_apply(i1,i1); const double k12 = cache_.kernel_apply(i1,i2); const double k22 = cache_.kernel_apply(i2,i2); const double eta = 2*k12-k11-k22; double a1,a2; if (eta < 0) { a2 = alpha2-y2*(E1-E2)/eta; a2 = median(L,a2,H); } else { double gamma = alpha1+s*alpha2; // Eq. (12.22) double v1=E1+y1+b-y1*alpha1*k11-y2*alpha2*k12; // Eq. (12.21) // Note that f(x1) = E1 + y1 double v2=E2+y2+b-y1*alpha1*k12-y2*alpha2*k22; // Eq. (12.21) double L_obj = gamma-s*L+L-.5*k11*(gamma-s*L)*(gamma-s*L)-.5*k22*L*L-s*k12*(gamma-s*L)*L-y1*(gamma-s*L)*v1-y2*L*v2; // # + W_const # Eq. (12.23) double H_obj = gamma-s*H+H-.5*k11*(gamma-s*H)*(gamma-s*L)-.5*k22*H*H-s*k12*(gamma-s*H)*H-y1*(gamma-s*H)*v1-y2*H*v2; // # + W_const # Eq. (12.23) if (L_obj > H_obj + eps) { a2 = L; } else if (L_obj < H_obj - eps) { a2 = H; } else { a2 = alpha2; } } if (a2 < tol) a2 = 0; else if (a2 > C-tol) a2 = C; if (std::abs(a2-alpha2) < eps*(a2+alpha2+eps)) return false; a1 = alpha1+s*(alpha2-a2); if (a1 < tol) a1 = 0; if (a1 > C-tol) a1 = C; // update everything Alphas[i1]=a1; Alphas[i2]=a2; double b1 = E1 + Y[i1]*(a1-alpha1)*k11+Y[i2]*(a2-alpha2)*k12+b; // Eq. (12.9) double b2 = E2 + Y[i1]*(a1-alpha1)*k12+Y[i2]*(a2-alpha2)*k22+b; // Eq. (12.10) const double new_b = (b1+b2)/2.; /* #for i in xrange(N): # if Alphas[i] in (0,C): # continue # elif i == i1 or i == i2: # E[i] = 0 # else: # E[i] += y1*(a1-alpha1)*kernel_apply(i1,i)+y2*(a2-alpha2)*kernel_apply(i2,i) + (b-new_b) # Eq. (12.11) #E[i1]=f_at(i1)-y1 #E[i2]=f_at(i2)-y2 */ b = new_b; return true; } bool SMO::examine_example(int i2) { //std::cout << "examine_example( " << i2 << " ) " << std::endl; const int y2 = Y[i2]; const double alpha2 = Alphas[i2]; const double E2 = get_error(i2); const double r2 = E2 * y2; //#print 'alpha2', alpha2, 'E2', E2, 'r2', r2 if ( ( (r2 < -tol) && (alpha2 < C) ) || ( (r2 > tol) && (alpha2 > 0) )){ int best_i1 = -1; double bestE = -1; for (int i = 0; i != N; ++i) { if (Alphas[i] != 0 && Alphas[i] != C) { double dE = E2-get_error(i); if (dE < 0.) dE = -dE; if (dE > bestE) { bestE=dE; best_i1 = i; } } } if (best_i1 != -1 && take_step(best_i1,i2)) return true; for (int i1 = 0; i1 != N; ++i1) { if (Alphas[i1] && Alphas[i1] != C && take_step(i1,i2)) return true; } for (int i1 = 0; i1 != N; ++i1) { if (take_step(i1,i2)) return true; } } return false; } void SMO::optimise() { b = 0; for (int i = 0; i != N; ++i) Alphas[i] = 0; int changed = 0; bool examineAll = true; //int iter = 0; while (changed || examineAll) { //std::cout << "SMO::optimize loop: " << iter++ << "\n"; check_for_interrupts(); changed = 0; for (int i = 0; i != N; ++i) { if (examineAll || (Alphas[i] != 0 && Alphas[i] != C)) { changed += examine_example(i); } } if (examineAll) examineAll = false; else if (!changed) examineAll = true; } } void assert_type_contiguous(PyArrayObject* array,int type) { if (!PyArray_Check(array) || PyArray_TYPE(array) != type || !PyArray_ISCONTIGUOUS(array)) { throw SMO_Exception("Arguments to eval_(SMO|LIBSVM) don't conform to expectation. Are you calling this directly? This is an internal function!"); } } PyObject* eval_SMO(PyObject* self, PyObject* args) { try { PyObject* X; PyArrayObject* Y; PyArrayObject* Alphas0; PyArrayObject* params; PyObject* kernel; int cache_size; if (!PyArg_ParseTuple(args, "OOOOOi", &X, &Y, &Alphas0, ¶ms,&kernel,&cache_size)) { const char* errmsg = "Arguments were not what was expected for eval_SMO.\n" "This is an internal function: Do not call directly unless you know what you're doing.\n"; PyErr_SetString(PyExc_RuntimeError,errmsg); return 0; } assert_type_contiguous(Y,NPY_INT32); assert_type_contiguous(Alphas0,NPY_DOUBLE); assert_type_contiguous(params,NPY_DOUBLE); if (PyArray_DIM(params,0) < 4) throw SMO_Exception("eval_SMO: Too few parameters"); int * Yv = static_cast(PyArray_DATA(Y)); double* Alphas = static_cast(PyArray_DATA(Alphas0)); unsigned N = PyArray_DIM(Y,0); double b = *static_cast(PyArray_GETPTR1(params,0)); double C = *static_cast(PyArray_GETPTR1(params,1)); double eps = *static_cast(PyArray_GETPTR1(params,2)); double tol = *static_cast(PyArray_GETPTR1(params,3)); SMO optimiser(X,Yv,Alphas,b,C,N,kernel,eps,tol,cache_size); optimiser.optimise(); *static_cast(PyArray_GETPTR1(params,0)) = b; // Write back b Py_RETURN_NONE; } catch (const Python_Exception&) { // if Python_Exception was thrown, then PyErr is already set. return 0; } catch (const SMO_Exception& exc) { PyErr_SetString(PyExc_RuntimeError,exc.msg); return 0; } catch (...) { PyErr_SetString(PyExc_RuntimeError,"Some sort of exception in eval_SMO."); return 0; } } // The code for LIBSVM_Solver is taken from LIBSVM and adapted to work well in milk // Copyright (c) 2000-2008 Chih-Chung Chang and Chih-Jen Lin // Changes were made to make it more similar to our formulation. // An SMO algorithm in Fan et al., JMLR 6(2005), p. 1889--1918 // Solves: // // min 0.5(\alpha^T Q \alpha) + p^T \alpha // // y^T \alpha = \delta // y_i = +1 or -1 // 0 <= alpha_i <= Cp for y_i = 1 // 0 <= alpha_i <= Cn for y_i = -1 // // Given: // // Q, p, y, Cp, Cn, and an initial feasible point \alpha // l is the size of vectors and matrices // eps is the stopping tolerance // // solution will be put in \alpha, objective value will be put in obj // #define info printf void info_flush() { } class LIBSVM_Solver { public: LIBSVM_Solver(PyObject* X, int* Y, double* Alphas, double* p, double& b, double C, int N, PyObject* kernel, double eps, double tol, int cache_size, bool shrinking): Alphas(Alphas), Y(Y), b(b), C(C), N(N), cache_(Y,get_kernel(X,kernel),N,cache_size), eps(eps), tol(tol), alpha_status(N), active_size(N), active_set(N), p(p), shrinking(shrinking), tau(eps) { } virtual ~LIBSVM_Solver() { } struct SolutionInfo { double obj; double rho; double upper_bound_p; double upper_bound_n; double r; // for LIBSVM_Solver_NU }; void optimise(); protected: double* Alphas; int* Y; double& b; double C; const int N; LIBSVM_KernelCache cache_; const double eps; const double tol; vector G; // gradient of objective function enum alpha_status_e { lower_bound, upper_bound, free}; vector alpha_status; int active_size; vector active_set; double *p; vector G_bar; // gradient, if we treat free variables as 0 bool unshrinked; // XXX bool shrinking; double tau; double get_C(int i) const { return C; //return (Y[i] > 0)? Cp : Cn; } void update_alpha_status(int i) { if(Alphas[i] >= get_C(i)) alpha_status[i] = upper_bound; else if(Alphas[i] <= 0) alpha_status[i] = lower_bound; else alpha_status[i] = free; } bool is_upper_bound(int i) const { return alpha_status[i] == upper_bound; } bool is_lower_bound(int i) const { return alpha_status[i] == lower_bound; } bool is_free(int i) const { return alpha_status[i] == free; } void swap_index(int i, int j); void reconstruct_gradient(); virtual bool select_working_set(int &i, int &j); virtual double calculate_rho(); virtual void do_shrinking(); private: bool be_shrunken(int i, double Gmax1, double Gmax2); void print_status() const; }; void LIBSVM_Solver::swap_index(int i, int j) { // We *do not* swap in the cache or kernel // Therefore *all acesses* to the cache or kernel // must be of the form cache_.get_kernel(active_set[i]) // instead of cache_.get_kernel(i)! std::swap(Y[i],Y[j]); std::swap(G[i],G[j]); std::swap(alpha_status[i],alpha_status[j]); std::swap(Alphas[i],Alphas[j]); std::swap(p[i],p[j]); std::swap(active_set[i],active_set[j]); std::swap(G_bar[i],G_bar[j]); } void LIBSVM_Solver::reconstruct_gradient() { // reconstruct inactive elements of G from G_bar and free variables if(active_size == N) return; for(int i=active_size;i= 0) && (i < active_size)); assert((j >= 0) && (j < active_size)); // update alpha[i] and alpha[j], handle bounds carefully const double *Q_i = cache_.get_kline(active_set[i], active_size); const double *Q_j = cache_.get_kline(active_set[j], active_size); const double C_i = get_C(i); const double C_j = get_C(j); const double old_alpha_i = Alphas[i]; const double old_alpha_j = Alphas[j]; if(Y[i]!=Y[j]) { double quad_coef = Q_i[active_set[i]]+Q_j[active_set[j]]+2.0*Q_i[active_set[j]]; if (quad_coef <= 0) quad_coef = tau; const double delta = (-G[i]-G[j])/quad_coef; const double diff = Alphas[i] - Alphas[j]; Alphas[i] += delta; Alphas[j] += delta; if(diff > 0) { if(Alphas[j] < 0) { Alphas[j] = 0; Alphas[i] = diff; } } else { if(Alphas[i] < 0) { Alphas[i] = 0; Alphas[j] = -diff; } } if(diff > C_i - C_j) { if(Alphas[i] > C_i) { Alphas[i] = C_i; Alphas[j] = C_i - diff; } } else { if(Alphas[j] > C_j) { Alphas[j] = C_j; Alphas[i] = C_j + diff; } } //print_status(); } else { double quad_coef = Q_i[active_set[i]]+Q_j[active_set[j]]-2*Q_i[active_set[j]]; if (quad_coef <= 0) quad_coef = tau; const double delta = (G[i]-G[j])/quad_coef; double sum = Alphas[i] + Alphas[j]; Alphas[i] -= delta; Alphas[j] += delta; if(sum > C_i) { if(Alphas[i] > C_i) { Alphas[i] = C_i; Alphas[j] = sum - C_i; } } else { if(Alphas[j] < 0) { Alphas[j] = 0; Alphas[i] = sum; } } if(sum > C_j) { if(Alphas[j] > C_j) { Alphas[j] = C_j; Alphas[i] = sum - C_j; } } else { if(Alphas[i] < 0) { Alphas[i] = 0; Alphas[j] = sum; } } } // update G const double delta_Alphas_i = Alphas[i] - old_alpha_i; const double delta_Alphas_j = Alphas[j] - old_alpha_j; // print_status(); for(int k=0;krho = calculate_rho(); b = calculate_rho(); // calculate objective value double v = 0; for(int i=0;iobj = v/2; // put back the solution for (int i = 0; i != N; ++i) { while (active_set[i] != i) { int j = active_set[i]; std::swap(Y[i],Y[j]); // It's not polite to clobber Y, so put it back std::swap(Alphas[i],Alphas[j]); std::swap(active_set[i],active_set[j]); } } //si->upper_bound_p = Cp; //si->upper_bound_n = Cn; } void LIBSVM_Solver::print_status() const { std::cout << " active_set_size: " << active_size << "\n"; for (int i = 0; i != N; ++i) { std::cout << " active_set[i]: " << active_set[i] << "\n"; std::cout << " p[i]: " << p[i] << "\n"; std::cout << " G[i]: " << G[i] << "\n"; std::cout << " G_bar[i]: " << G_bar[i] << "\n"; std::cout << " alpha_status[i]: " << alpha_status[i] << "\n"; std::cout << " Y[i]: " << Y[i] << "\n"; std::cout << " Alphas[i]: " << Alphas[i] << "\n"; std::cout << "\n"; } } // return true if already optimal, return false otherwise bool LIBSVM_Solver::select_working_set(int &out_i, int &out_j) { // return i,j such that // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha) // j: minimizes the decrease of obj value // (if quadratic coefficeint <= 0, replace it with tau) // -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha) //print_status(); double Gmax = -INF; double Gmax2 = -INF; int Gmax_idx = -1; int Gmin_idx = -1; double obj_diff_min = INF; for(int t=0;t= Gmax) { Gmax = -Y[t]*G[t]; Gmax_idx = t; } } } if (Gmax_idx == -1) return true; const int i = Gmax_idx; const double* Q_i = cache_.get_kline(active_set[i], N); const double* QDiag = cache_.get_diag(); for(int j=0;j= Gmax2) Gmax2 = YGj; if (grad_diff > 0) { const double quad_coef = Q_i[active_set[i]] + QDiag[active_set[j]] - 2*Y[j]*Y[i]*Q_i[active_set[j]]; const double obj_diff_factor = (quad_coef > 0) ? quad_coef : tau; const double obj_diff = -(grad_diff*grad_diff)/obj_diff_factor; if (obj_diff <= obj_diff_min) { Gmin_idx = j; obj_diff_min = obj_diff; } } } } if (Gmin_idx == -1) return true; if(Gmax+Gmax2 < eps) return true; out_i = Gmax_idx; out_j = Gmin_idx; return false; } bool LIBSVM_Solver::be_shrunken(int i, const double Gmax1, const double Gmax2) { if(is_upper_bound(i)) { if(Y[i]==+1) return -G[i] > Gmax1; return -G[i] > Gmax2; } else if(is_lower_bound(i)) { if (Y[i]==+1) return G[i] > Gmax2; return G[i] > Gmax1; } return false; } void LIBSVM_Solver::do_shrinking() { double Gmax1 = -INF; // max { -y_i * grad(f)_i | i in I_up(\alpha) } double Gmax2 = -INF; // max { y_i * grad(f)_i | i in I_low(\alpha) } // find maximal violating pair first for(int i=0;i i) { if (!be_shrunken(active_size, Gmax1, Gmax2)) { swap_index(i,active_size); break; } --active_size; } } } // unshrink, check all variables again before final iterations if(unshrinked || Gmax1 + Gmax2 > eps*10) return; unshrinked = true; reconstruct_gradient(); for(int i= N-1; i >= active_size; --i) { if (!be_shrunken(i, Gmax1, Gmax2)) { while (active_size < i) { if (be_shrunken(active_size, Gmax1, Gmax2)) { swap_index(i,active_size); break; } ++active_size; } ++active_size; } } } double LIBSVM_Solver::calculate_rho() { int nr_free = 0; double ub = INF, lb = -INF, sum_free = 0; for(int i=0;i0) return sum_free/nr_free; return (ub+lb)/2; } PyObject* eval_LIBSVM(PyObject* self, PyObject* args) { try { PyObject* X; PyArrayObject* Y; PyArrayObject* Alphas0; PyArrayObject* p; PyArrayObject* params; PyObject* kernel; int cache_size; if (!PyArg_ParseTuple(args, "OOOOOOi", &X, &Y, &Alphas0, &p, ¶ms,&kernel,&cache_size)) { const char* errmsg = "Arguments were not what was expected for eval_LIBSVM.\n" "This is an internal function: Do not call directly unless you know exactly what you're doing.\n"; PyErr_SetString(PyExc_RuntimeError,errmsg); return 0; } assert_type_contiguous(Y,NPY_INT32); assert_type_contiguous(Alphas0,NPY_DOUBLE); assert_type_contiguous(p,NPY_DOUBLE); assert_type_contiguous(params,NPY_DOUBLE); if (PyArray_DIM(params,0) < 4) throw SMO_Exception("eval_LIBSVM: Too few parameters"); int * Yv = static_cast(PyArray_DATA(Y)); double* Alphas = static_cast(PyArray_DATA(Alphas0)); double* pv = static_cast(PyArray_DATA(p)); unsigned N = PyArray_DIM(Y,0); double& b = *static_cast(PyArray_GETPTR1(params,0)); double C = *static_cast(PyArray_GETPTR1(params,1)); double eps = *static_cast(PyArray_GETPTR1(params,2)); double tol = *static_cast(PyArray_GETPTR1(params,3)); LIBSVM_Solver optimiser(X,Yv,Alphas,pv,b,C,N,kernel,eps,tol,cache_size, true); optimiser.optimise(); Py_RETURN_NONE; } catch (const Python_Exception&) { // if Python_Exception was thrown, then PyErr is already set. return 0; } catch (const SMO_Exception& exc) { PyErr_SetString(PyExc_RuntimeError,exc.msg); return 0; } catch (...) { PyErr_SetString(PyExc_RuntimeError,"Some sort of exception in eval_LIBSVM."); return 0; } } PyMethodDef methods[] = { {"eval_SMO",eval_SMO, METH_VARARGS , "Do NOT call directly.\n" }, {"eval_LIBSVM",eval_LIBSVM, METH_VARARGS , "Do NOT call directly.\n" }, {NULL, NULL,0,NULL}, }; const char * module_doc = "Internal SVM Module.\n" "\n" "Do NOT use directly!\n"; } // namespace extern "C" void init_svm() { import_array(); (void)Py_InitModule3("_svm", methods, module_doc); } milk-0.3.10/milk/supervised/normalise.py0000644000175000017500000000614311511426617021025 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division from collections import defaultdict import numpy as np from .classifier import normaliselabels __all__ = [ 'zscore', 'zscore_normalise', 'interval_normalise', 'chkfinite', 'sample_to_2min', 'normaliselabels' ] def zscore(features): """ features = zscore(features) Returns a copy of features which has been normalised to zscores """ mu = features.mean(0) sigma = np.std(features,0) sigma[sigma == 0] = 1 return (features - mu) / sigma class subtract_divide_model(object): def __init__(self, shift, factor): factor[factor == 0] = 1 # This makes the division a null op. self.shift = shift self.factor = factor def apply(self, features): return (features - self.shift)/self.factor def __repr__(self): return 'subtract_divide_model(%s, %s)' % (self.shift, self.factor) class zscore_normalise(object): ''' Normalise to z-scores A preprocessor that normalises features to z scores. ''' def train(self, features, labels, normalisedlabels=False): shift = features.mean(0) factor = np.std(features,0) return subtract_divide_model(shift, factor) class interval_normalise(object): ''' Linearly scale to the interval [-1,1] (per libsvm recommendation) ''' def train(self, features, labels, normalisedlabels=False): ptp = features.ptp(0) shift = features.min(0) + ptp/2. factor = ptp/2. return subtract_divide_model(shift, factor) def __repr__(self): return 'interval_normalise()' def sample_to_2min(labels): ''' selected = sample_to_2min(labels) Select examples so that the ratio of size of the largest class to the smallest class is at most two (i.e., min_label_count = min { (labels == L).sum() | for L in set(labels) } for L' in set(labels): assert (labels == L').sum() <= 2 * min_label_count ) Parameters ---------- * labels: sequence of labels Output ------ * selected: a Boolean numpy.ndarray ''' counts = defaultdict(int) for n in labels: counts[n] += 1 labels = np.asanyarray(labels) max_entries = np.min(counts.values())*2 selected = np.zeros(len(labels), bool) for c in counts.iterkeys(): p, = np.where(labels == c) p = p[:max_entries] selected[p] = 1 return selected class chkfinite(object): ''' Fill NaN & Inf values Replaces NaN & Inf values with zeros. ''' def __init__(self): pass def train(self, features, labels, normalisedlabels=False): return self def apply(self, features): nans = np.isnan(features) + np.isinf(features) if nans.any(): features = features.copy() features[nans] = 0 return features def __repr__(self): return 'chkfinite()' milk-0.3.10/milk/demos/0000755000175000017500000000000011562337213015373 5ustar luispedroluispedromilk-0.3.10/milk/demos/__init__.py0000644000175000017500000000000011513410326017463 0ustar luispedroluispedromilk-0.3.10/milk/demos/rf_wine_2d.py0000644000175000017500000000167711513656324020001 0ustar luispedroluispedrofrom milk.supervised import randomforest from milk.supervised.multi import one_against_one import milk.nfoldcrossvalidation import milk.unsupervised import pylab from milksets import wine # Load 'wine' dataset features, labels = wine.load() # random forest learner rf_learner = randomforest.rf_learner() # rf is a binary learner, so we transform it into a multi-class classifier learner = one_against_one(rf_learner) # cross validate with this learner and return predictions on left-out elements cmat,names, preds = milk.nfoldcrossvalidation(features, labels, classifier=learner, return_predictions=1) print 'cross-validation accuracy:', cmat.trace()/float(cmat.sum()) # dimensionality reduction for display x,v = milk.unsupervised.pca(features) colors = "rgb" # predicted colour marks = "xo" # whether the prediction was correct for (y,x),p,r in zip(x[:,:2], preds, labels): c = colors[p] m = marks[p == r] pylab.plot(y,x,c+m) pylab.show() milk-0.3.10/milk/demos/svm-decision-boundary.py0000644000175000017500000000222411457420270022165 0ustar luispedroluispedrofrom pylab import * import numpy as np from milksets.wine import load import milk.supervised import milk.unsupervised.pca import milk.supervised.svm features, labels = load() features = features[labels < 2] labels = labels[labels < 2] features,_ = milk.unsupervised.pca(features) features = features[:,:2] learner = milk.supervised.svm.svm_raw(kernel=np.dot, C=12) model = learner.train(features, labels) w = np.dot(model.svs.T, model.Yw) b = model.b x = np.linspace(-.5, .1, 100) y = -w[0]/w[1]*x + b/w[1] plot(features[labels == 1][:,0], features[labels == 1][:,1], 'bx') plot(features[labels == 0][:,0], features[labels == 0][:,1], 'ro') plot(x,y) clf() learner = milk.supervised.svm.svm_raw(kernel=milk.supervised.svm.rbf_kernel(1.), C=12) model = learner.train(features, labels) Y, X = (np.mgrid[:101,:101]-50)/12.5 values = [model.apply((y,x)) for y,x in zip(Y.ravel(),X.ravel())] values = np.array(values).reshape(Y.shape) sfeatures = features*12.5 sfeatures += 50 plot(sfeatures[labels == 0][:,0], sfeatures[labels == 0][:,1], 'bo') plot(sfeatures[labels == 1][:,0], sfeatures[labels == 1][:,1], 'ro') imshow(values.T) savefig('decision-boundary.png') milk-0.3.10/milk/demos/adaboost.py0000644000175000017500000000115011524607075017542 0ustar luispedroluispedroimport pylab as plt import milk.supervised.tree import milk.supervised.adaboost from milksets import wine import milk.supervised.multi weak = milk.supervised.tree.stump_learner() learner = milk.supervised.adaboost.boost_learner(weak) learner = milk.supervised.multi.one_against_one(learner) features, labels = wine.load() cmat,names,predictions = milk.nfoldcrossvalidation(features,labels, classifier=learner, return_predictions=True) colors = "rgb" codes = "xo" for y,x,r,p in zip(features.T[0], features.T[1], labels, predictions): code = codes[int(r == p)] plt.plot([y],[x], colors[p]+code) plt.show() milk-0.3.10/milk/milk_version.py0000644000175000017500000000002711562335225017337 0ustar luispedroluispedro__version__ = '0.3.10' milk-0.3.10/milk/__init__.py0000644000175000017500000000415511556634056016412 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2010, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. ''' Milk Machine learning in Python Toplevel functions ------------------ - nfoldcrossvalidation: n-fold crossvalidation - defaultclassifier: get a general purpose classifier - kmeans: kmeans clustering Modules ------- - supervised - unsupervised - measures Example ------- :: features = np.random.randn(100,20) features[:50] *= 2 labels = np.repeat((0,1), 50) classifier = milk.defaultclassifier() model = classifier.train(features, labels) new_label = model.apply(np.random.randn(100)) new_label2 = model.apply(np.random.randn(100)*2) ''' from .nfoldcrossvalidation import nfoldcrossvalidation from .supervised.defaultclassifier import defaultclassifier from .supervised.defaultlearner import defaultlearner from .unsupervised.kmeans import kmeans from milk_version import __version__ __all__ = [ '__version__', 'kmeans', 'defaultclassifier', 'defaultlearner', 'nfoldcrossvalidation', ] milk-0.3.10/milk/tests/0000755000175000017500000000000011562337213015426 5ustar luispedroluispedromilk-0.3.10/milk/tests/test_regression_constant_features.py0000644000175000017500000000037311556634056025041 0ustar luispedroluispedroimport milk import numpy as np def test_constant_features(): learner = milk.defaultclassifier() features = np.ones(20).reshape((-1,1)) labels = np.zeros(20) labels[10:] += 1 features[10:] *= -1 learner.train(features, labels) milk-0.3.10/milk/tests/test_svm.py0000644000175000017500000001477211452411266017656 0ustar luispedroluispedroimport milk.supervised.svm from milk.supervised.svm import svm_learn_smo, svm_learn_libsvm, _svm_apply, learn_sigmoid_constants, svm_sigmoidal_correction import numpy import random import numpy as np import milk.supervised.svm import pickle eps=1e-3 def approximate(a,b): a=numpy.asanyarray(a) b=numpy.asanyarray(b) return numpy.abs(a-b).max() < eps def assert_kkt(SVM): X,Y,Alphas,b,C,kernel=SVM N=len(Alphas) for i in xrange(N): if Alphas[i] == 0.: assert Y[i]*_svm_apply(SVM,X[i])+eps >= 1 elif Alphas[i] == C: assert Y[i]*_svm_apply(SVM,X[i])-eps <= 1 else: assert abs(Y[i]*_svm_apply(SVM,X[i])-1) <= eps def assert_all_correctly_classified(SVM,X,Y): N = len(X) for i in xrange(N): assert _svm_apply(SVM,X[i]) * Y[i] > 0 def assert_more_than_50(SVM,X,Y): N = len(X) correct = 0 for i in xrange(N): correct += (_svm_apply(SVM,X[i]) * Y[i] > 0) assert correct > N/2 def test_simplest(): X=numpy.array([ [1], [2], ]) Y=numpy.array([1,-1]) C=4. kernel=numpy.dot Alphas,b=svm_learn_smo(X,Y,kernel,C) SVM=(X,Y,Alphas,b,C,kernel) Alphas_,b_=svm_learn_libsvm(X,Y,kernel,C) assert approximate(Alphas,Alphas_) assert approximate(b,b_) assert approximate(Alphas,[2.,2.]) assert approximate(b,-3) assert_kkt(SVM) assert_all_correctly_classified(SVM,X,Y) def test_more_complex_kkt(): X=numpy.array([ [1,0], [2,1], [2,0], [4,2], [0,1], [0,2], [1,2], [2,8]]) Y=numpy.array([1,1,1,1,-1,-1,-1,-1]) C=4. kernel=numpy.dot Alphas,b=svm_learn_smo(X,Y,kernel,C) def verify(): SVM=(X,Y,Alphas,b,C,kernel) sv=numpy.array([1,1,0,0,1,0,1,0]) nsv=~sv computed_sv = (Alphas > 0) & (Alphas < C) computed_nsv = ~computed_sv assert_kkt(SVM) assert numpy.all((sv-computed_sv) >= 0) # computed_sv in sv assert numpy.all((computed_nsv-nsv) >= 0) # nsv in computed_nsv assert_all_correctly_classified(SVM,X,Y) verify() Alphas,b=svm_learn_libsvm(X,Y,kernel,C) verify() def rbf(xi,xj): return numpy.exp(-((xi-xj)**2).sum()) def test_rbf(): X=numpy.array([ [0,0,0], [1,1,1], ]) Y=numpy.array([ 1, -1 ]) C=10 Alphas,b=svm_learn_smo(X,Y,rbf,C) SVM=(X,Y,Alphas,b,C,rbf) assert_all_correctly_classified(SVM,X,Y) Alphas,b=svm_learn_libsvm(X,Y,rbf,C) SVM=(X,Y,Alphas,b,C,rbf) assert_all_correctly_classified(SVM,X,Y) def test_random(): R=numpy.random.RandomState(123) X=R.rand(10,3) X[:5]+=1. C=2 Y=numpy.ones(10) Y[5:] *= -1 Alphas,b=svm_learn_smo(X,Y,rbf,C) SVM=(X,Y,Alphas,b,C,rbf) assert_more_than_50(SVM,X,Y) Alphas,b=svm_learn_libsvm(X,Y,rbf,C) SVM=(X,Y,Alphas,b,C,rbf) assert_more_than_50(SVM,X,Y) def test_platt_correction(): F=numpy.ones(10) L=numpy.ones(10) F[:5] *= -1 L[:5] *= 0 A,B=learn_sigmoid_constants(F,L) assert 1./(1.+numpy.exp(+10*A+B)) > .99 assert 1./(1.+numpy.exp(-10*A+B)) <.01 def test_platt_correction_class(): F=numpy.ones(10) L=numpy.ones(10) F[:5] *= -1 L[:5] *= 0 corrector = svm_sigmoidal_correction() A,B = learn_sigmoid_constants(F,L) model = corrector.train(F,L) assert model.A == A assert model.B == B assert model.apply(10) > .99 assert model.apply(-10) < .01 def test_perfect(): data = numpy.zeros((10,2)) data[5:] = 1 labels = numpy.zeros(10) labels[5:] = 1 classifier=milk.supervised.svm.svm_raw(kernel=milk.supervised.svm.rbf_kernel(1),C=4) model = classifier.train(data,labels) assert numpy.all( (numpy.array([model.apply(data[i]) for i in xrange(10)]) > 0) == labels ) def test_smart_rbf(): import milksets.wine features,labels = milksets.wine.load() labels = (labels == 1) kernel = milk.supervised.svm.rbf_kernel(2.**-4) C = milk.supervised.svm.svm_raw(C=2.,kernel=kernel) model = C.train(features,labels) smartkernel = [model.apply(f) for f in features] del kernel.kernel_nr_ del kernel.kernel_arg_ C = milk.supervised.svm.svm_raw(C=2.,kernel=kernel) model = C.train(features,labels) dumbkernel = [model.apply(f) for f in features] np.allclose(smartkernel, dumbkernel) def test_preprockernel(): def test_preproc(kernel): X = np.random.rand(200,20) Qs = np.random.rand(100, 20) prekernel = kernel.preprocess(X) Qs = np.random.rand(10, 20) for q in Qs: preprocval = prekernel(q) procval = np.array([kernel(q,x) for x in X]) assert np.allclose(preprocval, procval) yield test_preproc, milk.supervised.svm.rbf_kernel(2.) def test_svm_to_binary(): base = milk.supervised.svm.svm_raw(kernel=np.dot, C=4.) bin = milk.supervised.svm.svm_to_binary(base) X = np.array([ [1,0], [2,1], [2,0], [4,2], [0,1], [0,2], [1,2], [2,8]]) Y = np.array([1,1,1,1,-1,-1,-1,-1]) model = bin.train(X,Y) modelbase = base.train(X,Y) assert np.all(modelbase.Yw == model.models[0].Yw) for x in X: assert model.apply(x) in set(Y) def test_fast_dotkernel(): base = milk.supervised.svm.svm_raw(kernel=np.dot, C=4.) X = np.array([ [1,0], [2,1], [2,0], [4,2], [0,1], [0,2], [1,2], [2,8]]) Y = np.array([1,1,1,1,-1,-1,-1,-1]) model = base.train(X,Y) base = milk.supervised.svm.svm_raw(kernel = milk.supervised.svm.dot_kernel(), C = 4.) model2 = base.train(X,Y) assert np.all( model2.Yw == model.Yw ) def count_common_elements(s0,s1): return len(s0.intersection(s1)) def test_not_ndarray_input(): elems = [ set([1,2,3]), set([2,3]), set([2,]), set([0,1]), set([0,2,3]), set([1,2,3]), set([2,3,45]), set([1,2,100])] has_one = [(1 in s) for s in elems] has_one[0] = not has_one[0] c = milk.supervised.svm.svm_raw(kernel = count_common_elements, C = 2.) model = c.train(elems, has_one) trainpreds = [model.apply(e) for e in elems] assert np.mean((np.array(trainpreds) > 0 ) == has_one ) > .5 model = pickle.loads(pickle.dumps(model)) trainpreds = [model.apply(e) for e in elems] assert np.mean((np.array(trainpreds) > 0 ) == has_one ) > .5 milk-0.3.10/milk/tests/test_pca.py0000644000175000017500000000073011231643010017566 0ustar luispedroluispedroimport numpy.random import milk.unsupervised.pca import numpy as np def test_pca(): numpy.random.seed(123) X = numpy.random.rand(10,4) X[:,1] += numpy.random.rand(10)**2*X[:,0] X[:,1] += numpy.random.rand(10)**2*X[:,0] X[:,2] += numpy.random.rand(10)**2*X[:,0] Y,V = milk.unsupervised.pca(X) Xn = milk.unsupervised.normalise.zscore(X) assert X.shape == Y.shape assert ((np.dot(V[:4].T,Y[:,:4].T).T-Xn)**2).sum()/(Xn**2).sum() < .3 milk-0.3.10/milk/tests/test_affinity.py0000644000175000017500000000055611452411266020655 0ustar luispedroluispedroimport milk.unsupervised.affinity import numpy as np def test_affinity(): np.random.seed(22) X = np.random.randn(100,10) X[:40] += .4 S = milk.unsupervised.pdist(X) clusters, labels = milk.unsupervised.affinity.affinity_propagation(S) assert labels.max()+1 == len(clusters) assert len(labels) == len(X) assert clusters.max() < len(X) milk-0.3.10/milk/tests/test_precluster_learner.py0000644000175000017500000000143611556634056022753 0ustar luispedroluispedroimport numpy as np from milk.supervised.precluster_learner import precluster_learner from milk.tests.fast_classifier import fast_classifier def test_precluster(): np.random.seed(22) learner = precluster_learner([2], base=fast_classifier()) def c0(): return np.random.rand(8) def c1(): return c0()+2.*np.ones(8) features = [] labels =[] for i in xrange(200): f = [] for j in xrange(40): use_0 = (i < 100 and j < 30) or (i >= 100 and j >= 30) if use_0: f.append(c0()) else: f.append(c1()) labels.append((i < 100)) features.append(f) model = learner.train(features,labels) assert model.apply([c0() for i in xrange(35)]) assert not model.apply([c1() for i in xrange(35)]) milk-0.3.10/milk/tests/test_knn.py0000644000175000017500000000155111277730267017640 0ustar luispedroluispedroimport numpy as np import milk.supervised.knn import numpy def test_simple(): X=np.array([ [0,0,0], [1,1,1], ]) Y=np.array([ 1, -1 ]) kNN = milk.supervised.knn.kNN(1) kNN = kNN.train(X,Y) assert kNN.apply(X[0]) == Y[0] assert kNN.apply(X[1]) == Y[1] assert kNN.apply([0,0,1]) == Y[0] assert kNN.apply([0,1,1]) == Y[1] def test_nnclassifier(): labels=[0,1] data=[[0.,0.],[1.,1.]] C = milk.supervised.knn.kNN(1) model = C.train(data,labels) assert model.apply(data[0]) == 0 assert model.apply(data[1]) == 1 assert model.apply([.01,.01]) == 0 assert model.apply([.99,.99]) == 1 assert model.apply([100,100]) == 1 assert model.apply([-100,-100]) == 0 assert model.apply([.9,.9]) == 1 middle = model.apply([.5,.5]) assert (middle == 0) or (middle == 1) milk-0.3.10/milk/tests/test_featureselection.py0000644000175000017500000000254411453673726022420 0ustar luispedroluispedroimport milk.supervised.featureselection import numpy as np def test_sda(): from milksets import wine features, labels = wine.load() selected = milk.supervised.featureselection.sda(features,labels) for sel in selected: assert sel <= features.shape[1] def test_linear_independent_features(): np.random.seed(122) X3 = np.random.rand(20,3) X = np.c_[X3,X3*2+np.random.rand(20,3)/20.,-X3*2+np.random.rand(20,3)/10.] X2 = np.c_[X3,X3*2,-X3*3e-3] assert len(milk.supervised.featureselection.linear_independent_features(X)) == 9 assert len(milk.supervised.featureselection.linear_independent_features(X2)) == 3 assert np.all (np.sort(milk.supervised.featureselection.linear_independent_features(X2) % 3) == np.arange(3)) def _rank(A,tol=1e-8): s = np.linalg.svd(A,compute_uv=0) return (s > tol).sum() def _slow_linear_independent_features(featmatrix): ''' Returns the indices of a set of linearly independent features (columns). indices = linear_independent_features(features) ''' independent = [0,] rank = 1 feat = [featmatrix[:,0]] for i,col in enumerate(featmatrix.T): feat.append(col) nrank = _rank(np.array(feat)) if nrank == rank: del feat[-1] else: rank = nrank independent.append(i) return np.array(independent) milk-0.3.10/milk/tests/test_parzen.py0000644000175000017500000000155611231643010020331 0ustar luispedroluispedrofrom __future__ import division import milk.supervised.normalise from milk.supervised.parzen import get_parzen_rbf_loocv import numpy as np import milksets def _slow_parzen(features, labels, sigma): correct = 0 N = len(features) labels = 2*labels - 1 def kernel(fi, fj): return np.exp(-((fi-fj)**2).sum()/sigma) for i in xrange(N): C = 0. for j in xrange(N): if i == j: continue C += labels[j] * kernel(features[i],features[j]) if (C*labels[i] > 0): correct += 1 return correct/N def test_parzen(): features,labels = milksets.wine.load() labels = (labels == 1) features = milk.supervised.normalise.zscore(features) f = get_parzen_rbf_loocv(features, labels) sigmas = 2.**np.arange(-4,4) for s in sigmas: assert abs(_slow_parzen(features, labels, s) - f(s)) < 1e-6 milk-0.3.10/milk/tests/test_measures_clusters.py0000644000175000017500000000102511522025760022602 0ustar luispedroluispedroimport milk.measures.cluster_agreement import numpy as np def test_rand_arand_jaccard(): np.random.seed(33) labels = np.repeat(np.arange(4),10) clusters = np.repeat(np.arange(4),10) a0,b0,c0= milk.measures.cluster_agreement.rand_arand_jaccard(clusters, labels) assert a0 == 1. assert b0 == 1. np.random.shuffle(clusters) a1,b1,c1= milk.measures.cluster_agreement.rand_arand_jaccard(clusters, labels) assert a1 >= 0. assert a1 < 1. assert b1 < 1. assert b1 >= 0. assert c1 < c0 milk-0.3.10/milk/tests/test_pdist.py0000644000175000017500000000074011452411266020162 0ustar luispedroluispedroimport numpy as np from milk.unsupervised import pdist, plike def test_pdist(): np.random.seed(222) X = np.random.randn(100,23) Y = np.random.randn(80,23) Dxx = pdist(X) for i in xrange(X.shape[0]): for j in xrange(X.shape[1]): assert np.allclose(Dxx[i,j], np.sum((X[i]-X[j])**2)) Dxy = pdist(X,Y) for i in xrange(X.shape[0]): for j in xrange(Y.shape[1]): assert np.allclose(Dxy[i,j], np.sum((X[i]-Y[j])**2)) milk-0.3.10/milk/tests/test_kmeans.py0000644000175000017500000000314611537446661020333 0ustar luispedroluispedroimport numpy as np import milk.unsupervised from milk.unsupervised.kmeans import assign_centroids def test_kmeans(): np.random.seed(132) features = np.r_[np.random.rand(20,3)-.5,.5+np.random.rand(20,3)] def test_distance(dist, kwargs={}): centroids, _ = milk.unsupervised.kmeans(features, 2, distance=dist, **kwargs) positions = [0]*20 + [1]*20 correct = (centroids == positions).sum() assert correct >= 38 or correct <= 2 yield test_distance, 'euclidean' yield test_distance, 'seuclidean' yield test_distance, 'mahalanobis', { 'icov' : np.eye(3) } def test_kmeans_centroids(): np.random.seed(132) features = np.random.rand(201,30) for k in [2,3,5,10]: indices,centroids = milk.unsupervised.kmeans(features, k) for i in xrange(k): if np.any(indices == i): assert np.allclose(centroids[i], features[indices == i].mean(0)) def test_assign_cids(): from milksets.wine import load features,_ = load() assigns, centroids = milk.unsupervised.kmeans(features, 3, R=2, max_iters=10) assert np.all(assign_centroids(features, centroids) == assigns) def test_non_contiguous_fmatrix(): from milksets.wine import load features,_ = load() features = features[:,::2] assigns, centroids = milk.unsupervised.kmeans(features, 3, R=2, max_iters=10) assert np.all(assign_centroids(features, centroids) == assigns) features = features.astype(np.int32) assigns, centroids = milk.unsupervised.kmeans(features, 3, R=2, max_iters=10) assert np.all(assign_centroids(features, centroids) == assigns) milk-0.3.10/milk/tests/test_rf.py0000644000175000017500000000060611462317362017452 0ustar luispedroluispedrofrom milk.supervised import randomforest import numpy as np def test_rf(): from milksets import wine features, labels = wine.load() features = features[labels < 2] labels = labels[labels < 2] learner = randomforest.rf_learner() model = learner.train(features[::5], labels[::5]) test = [model.apply(f) for f in features] assert np.mean(labels == test) > .7 milk-0.3.10/milk/tests/data/0000755000175000017500000000000011562337213016337 5ustar luispedroluispedromilk-0.3.10/milk/tests/data/jugparallel_jugfile.py0000644000175000017500000000056311540477424022731 0ustar luispedroluispedroimport milk.ext.jugparallel from milksets.wine import load from milk.tests.fast_classifier import fast_classifier features,labels = load() classified = milk.ext.jugparallel.nfoldcrossvalidation(features, labels, learner=fast_classifier()) classified_wpred = milk.ext.jugparallel.nfoldcrossvalidation(features, labels, learner=fast_classifier(), return_predictions=True) milk-0.3.10/milk/tests/data/__init__.py0000644000175000017500000000000011231643007020430 0ustar luispedroluispedromilk-0.3.10/milk/tests/data/jugparallel_kmeans_jugfile.py0000644000175000017500000000035511556634056024271 0ustar luispedroluispedroimport milk.ext.jugparallel from milksets.wine import load from milk.tests.fast_classifier import fast_classifier features,labels = load() clustered = milk.ext.jugparallel.kmeans_select_best(features, ks=(2,8), repeats=2, max_iters=6) milk-0.3.10/milk/tests/test_logistic.py0000644000175000017500000000060411503576164020661 0ustar luispedroluispedroimport milk.supervised.logistic import milksets.wine import numpy as np def test_better_than_random(): learner = milk.supervised.logistic.logistic_learner() features, labels = milksets.wine.load() model = learner.train(features, labels == 0) error = np.array([np.abs(model.apply(f)-(l == 0)) for f,l in zip(features, labels)]) assert error.mean() < .1 milk-0.3.10/milk/tests/test_nfoldcrossvalidation.py0000644000175000017500000001034511521632267023273 0ustar luispedroluispedroimport milk.measures.nfoldcrossvalidation from milk.measures.nfoldcrossvalidation import nfoldcrossvalidation, foldgenerator import milk.supervised.tree import numpy as np from fast_classifier import fast_classifier def test_foldgenerator(): labels = np.array([1]*20+[2]*30+[3]*20) for nf in [None,2,3,5,10,15,20]: assert np.array([test.copy() for _,test in foldgenerator(labels,nf)]).sum(0).max() == 1 assert np.array([test.copy() for _,test in foldgenerator(labels,nf)]).sum() == len(labels) assert np.array([(test&train).sum() for train,test in foldgenerator(labels,nf)]).sum() == 0 def test_foldgenerator_not_empty(): for nf in (None, 2, 3, 5, 10, 15, 20): for Tr,Te in foldgenerator([0] * 10 + [1] *10, nf, None): assert not np.all(Tr) assert not np.all(Te) def test_nfoldcrossvalidation_simple(): from milksets import wine features, labels = wine.load() features = features[::2] labels = labels[::2] cmat,clabels = nfoldcrossvalidation(features, labels, classifier=fast_classifier()) assert cmat.shape == (3,3) assert len(clabels) == 3 class test_classifier(object): def __init__(self,N): self.tested = np.zeros(N,bool) self.N = N def apply(self, features): cur = np.zeros(self.N,bool) cur[features] = True assert not np.any(cur & self.tested) self.tested |= cur return np.zeros_like(features) def train(self,f,l): return self def test_nfoldcrossvalidation_testall(): N = 121 C = test_classifier(N) features = np.arange(N) labels = np.zeros(N) cmat,clabels = nfoldcrossvalidation(features, labels, classifier=C) assert np.all(C.tested) def test_getfold(): A = np.zeros(20) A[:10] = 1 t,s = milk.measures.nfoldcrossvalidation.getfold(A,0,10) tt,ss = milk.measures.nfoldcrossvalidation.getfold(A,1,10) assert not np.any((~t)&(~tt)) def test_nfoldcrossvalidation_defaultclassifier(): np.random.seed(2233) X = np.random.rand(60,5) X[:20] += 4. X[-20:] -= 4. Y = np.ones(60) Y[:20] = 0 Y[-20:] = 2 Y += 100 cmat,clabels = milk.measures.nfoldcrossvalidation.nfoldcrossvalidation(X,Y) assert cmat.shape == (3,3) clabels.sort() assert np.all(clabels == [100,101,102]) def test_foldgenerator_origins(): def test_origins(labels, origins): for nf in (2,3,5,7): assert np.array([test.copy() for _,test in foldgenerator(labels, nf, origins)]).sum(0).max() == 1 assert np.array([test.copy() for _,test in foldgenerator(labels, nf, origins)]).sum() == len(labels) assert np.min([test.sum() for _,test in foldgenerator(labels, nf, origins)]) > 0 assert np.min([train.sum() for train,_ in foldgenerator(labels, nf, origins)]) > 0 for Tr,Te in foldgenerator(labels, nf, origins): assert not np.any(Tr&Te) in_test = set(origins[Te]) in_train = set(origins[Tr]) assert len(in_train.intersection(in_test)) == 0 tested = np.zeros(len(labels)) for Tr,Te in foldgenerator(labels, nf, origins): tested[Te] += 1 assert np.all(tested == 1) labels = np.zeros(120, np.uint8) labels[39:] += 1 labels[66:] += 1 origins = np.repeat(np.arange(40), 3) yield test_origins, labels, origins reorder = np.argsort(np.random.rand(len(labels))) labels = labels[reorder] origins = origins[reorder] yield test_origins, labels, origins def test_stringlabels(): np.random.seed(222) D = np.random.rand(100,10) D[:40] += np.random.rand(40,10)**2 labelnames = ['one'] * 40 + ['two'] * 60 cmat,Lo = nfoldcrossvalidation(D, labelnames, classifier=fast_classifier()) assert Lo[0] in labelnames assert Lo[1] in labelnames assert Lo[0] != Lo[1] in labelnames def test_predictions(): np.random.seed(222) D = np.random.rand(100,10) D[:40] += np.random.rand(40,10)**2 labels = [0] * 40 + [1] * 60 cmat,_,predictions = nfoldcrossvalidation(D, labels, classifier=fast_classifier(), return_predictions=1) assert np.all((predictions == 0)|(predictions == 1)) assert cmat.trace() == np.sum(predictions == labels) milk-0.3.10/milk/tests/test_adaboost.py0000644000175000017500000000076611470616152020644 0ustar luispedroluispedroimport numpy as np import milk.supervised.tree import milk.supervised.adaboost def test_learner(): from milksets import wine learner = milk.supervised.adaboost.boost_learner(milk.supervised.tree.stump_learner()) features, labels = wine.load() features = features[labels < 2] labels = labels[labels < 2] == 0 labels = labels.astype(int) model = learner.train(features, labels) train_out = np.array(map(model.apply, features)) assert (train_out == labels).mean() > .9 milk-0.3.10/milk/tests/test_gridsearch.py0000644000175000017500000000460711461573003021156 0ustar luispedroluispedroimport milk.supervised.gridsearch import milk.supervised.svm from milk.supervised.gridsearch import gridminimise, _allassignments import numpy as np def slow_gridminimise(learner, features, labels, params, measure=None): from ..measures.nfoldcrossvalidation import nfoldcrossvalidation if measure is None: measure = np.trace best_val = initial_value best = None for assignement in _allassignments(params): _set_assignment(learner, assignement) S,_ = nfoldcrossvalidation(features, labels, classifier=learner) cur = measure(S) if cur > best_val: best = assignement best_val = cur return best def test_gridsearch(): from milksets import wine features, labels = wine.load() selected = (labels < 2) features = features[selected] labels = labels[selected] G = milk.supervised.gridsearch( milk.supervised.svm.svm_raw(), params={'C':[.01,.1,1.,10.], 'kernel':[milk.supervised.svm.rbf_kernel(0.1),milk.supervised.svm.rbf_kernel(1.)] }) model = G.train(features,labels) reslabels = [model.apply(f) for f in features] assert len(reslabels) == len(features) test_gridsearch.slow = True def test_all_assignements(): assert len(list(_allassignments({'C': [0,1], 'kernel' : ['a','b','c']}))) == 2 * 3 class simple_model: def __init__(self, c): self.c = c def apply(self, f): return self.c def f(a,b,c): return a**2 + b**3 + c class simple_learner: def set_option(self, k, v): setattr(self, k, v) def train(self, fs, ls, normalisedlabels=False): return simple_model(f(self.a, self.b, self.c)) def test_gridminimise(): features = np.arange(100) labels = np.tile((0,1), 50) best = gridminimise(simple_learner(), features, labels, {'a': np.arange(4), 'b' : np.arange(-3,3), 'c' : np.linspace(2., 10) }, measure=(lambda _, p: p[0])) best = dict(best) val = f(best['a'], best['b'], best['c']) for a in np.arange(4): for b in np.arange(-3,3): for c in np.linspace(2., 10): assert val <= f(a,b,c) def test_gridminimise(): from milksets.wine import load features, labels = load() x = gridminimise(milk.supervised.svm_simple(kernel=np.dot, C=2.), features[::2], labels[::2] == 0, {'C' : (0.5,) }) cval, = x assert cval == ('C', .5) milk-0.3.10/milk/tests/test_normalise.py0000644000175000017500000000565211453127167021044 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2010, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. from __future__ import division import numpy import numpy as np from milk.supervised.normalise import sample_to_2min import milk.supervised.normalise def test_zscore_normalise(): I=milk.supervised.normalise.zscore_normalise() numpy.random.seed(1234) features = numpy.random.rand(20,100) L = numpy.zeros(100) model = I.train(features, L) transformed = np.array([model.apply(f) for f in features]) assert np.all( transformed.mean(0)**2 < 1e-7 ) assert np.all( np.abs(transformed.std(0) - 1) < 1e-3 ) def test_sample_to_2min(): A = np.zeros(256, np.int32) def test_one(A): selected = sample_to_2min(A) ratios = [] for l0 in set(A): for l1 in set(A): ratios.append( (A[selected] == l0).sum() / (A[selected] == l1).sum() ) assert np.max(ratios) <= 2.001 A[20:] = 1 yield test_one, A A[21:] = 1 yield test_one, A A[129:] = 2 yield test_one, A def test_sample_to_2min_list(): from collections import defaultdict def count(xs): counts = defaultdict(int) for x in xs: counts[x] += 1 return counts labels = ["A"]*8 + ["B"]*12 + ["C"]*16 + ["D"] * 24 + ["E"] * 1000 selected = sample_to_2min(labels) before = count(labels) after = count(np.array(labels)[selected]) assert max(after.values()) == min(before.values())*2 def test_interval_normalise(): interval = milk.supervised.normalise.interval_normalise() np.random.seed(105) features = np.random.randn(100, 5) model = interval.train(features, features[0] > 0) transformed = np.array([model.apply(f) for f in features]) assert np.allclose(transformed.min(0), -1) assert np.allclose(transformed.max(0), +1) milk-0.3.10/milk/tests/test_normaliselabels.py0000644000175000017500000000066711452411266022223 0ustar luispedroluispedrofrom milk.supervised.classifier import normaliselabels import numpy as np def test_normaliselabels(): np.random.seed(22) labels = np.zeros(120, np.uint8) labels[40:] += 1 labels[65:] += 1 reorder = np.argsort(np.random.rand(len(labels))) labels = labels[reorder] labels2,names = normaliselabels(labels) for new_n,old_n in enumerate(names): assert np.all( (labels == old_n) == (labels2 == new_n) ) milk-0.3.10/milk/tests/test_regression.py0000644000175000017500000000305011536473671021227 0ustar luispedroluispedroimport pickle import numpy as np import milk.supervised._svm from gzip import GzipFile from os import path import numpy as np from milksets.wine import load from milk.supervised import defaultclassifier import milk def test_svm_crash(): X,Y,kernel, C, eps ,tol, = pickle.load(GzipFile(path.dirname(__file__) + '/data/regression-2-Dec-2009.pp.gz')) X = X[2:-2,:].copy() Y = Y[2:-2].copy() N = len(Y) Y = Y.astype(np.int32) p = -np.ones(N,np.double) params = np.array([0,C,eps,tol],np.double) Alphas0 = np.zeros(N, np.double) cache_size = (1<<20) # The line below crashed milk: milk.supervised._svm.eval_LIBSVM(X,Y,Alphas0,p,params,kernel,cache_size) # HASN'T CRASHED! def test_nov2010(): # Bug submitted by Mao Ziyang # This was failing in 0.3.5 because SDA selected no features np.random.seed(222) features = np.random.randn(100,20) features[:50] *= 2 labels = np.repeat((0,1), 50) classifier = milk.defaultclassifier() model = classifier.train(features, labels) new_label = model.apply(np.random.randn(20)*2) new_label2 = model.apply(np.random.randn(20)) assert new_label == 0 assert new_label2 == 1 def test_default_small(): features, labels = load() selected = np.concatenate( [np.where(labels < 2)[0], np.where(labels == 2)[0][:6]] ) features = features[selected] labels = labels[selected] learner = defaultclassifier('fast') # For version 0.3.8, the line below led to an error milk.nfoldcrossvalidation(features, labels, classifier=learner) milk-0.3.10/milk/tests/test_grouped.py0000644000175000017500000000462411452411266020511 0ustar luispedroluispedroimport numpy as np import milk.supervised.svm from milk.supervised.svm import rbf_kernel import milk.supervised.multi import milk.supervised.grouped from milk.supervised.classifier import ctransforms import milksets.wine def group(features, labels, step): N = len(labels) i = 0 gfeatures = [] glabels = [] while i < N: next = i + step while next > N or labels[next-1] != labels[i]: next -= 1 gfeatures.append(features[i:next]) glabels.append(labels[i]) i = next return gfeatures, glabels def test_voting(): base = ctransforms(milk.supervised.svm.svm_raw(C=2.,kernel=milk.supervised.svm.rbf_kernel(2.**-3)),milk.supervised.svm.svm_binary()) base = milk.supervised.multi.one_against_rest(base) features,labels = milksets.wine.load() gfeatures, glabels = group(features, labels, 3) learner = milk.supervised.grouped.voting_classifier(base) learner.train(gfeatures, glabels) model = learner.train(gfeatures, glabels) assert ([model.apply(f) for f in gfeatures] == np.array(glabels)).mean() > .8 def test_filter_outliers(): np.random.seed(22) features = [np.random.randn(10,10) for i in xrange(20)] for f in features: f[0] *= 10 trainer = milk.supervised.grouped.filter_outliers(.9) model = trainer.train(features, [0] * len(features)) for f in features: ff = model.apply(f) assert np.all(ff == f[1:]) def test_nfoldcrossvalidation(): np.random.seed(22) features = np.array([np.random.rand(8+(i%3), 12)*(i//20) for i in xrange(40)], dtype=object) labels = np.zeros(40, int) labels[20:] = 1 classifier = milk.supervised.grouped.voting_classifier(milk.supervised.svm_simple(C=1., kernel=rbf_kernel(1./12))) cmat, names = milk.nfoldcrossvalidation(features, labels, classifier=classifier) assert cmat.shape == (2,2) assert sorted(names) == range(2) class identity_classifier(object): def train(self, features, labels): return identity_model() class identity_model(object): def apply(self, f): return f def test_meanclassif(): gfeatures = [np.arange(10), np.arange(10)%2] glabels = [0,1] meanclassif = milk.supervised.grouped.mean_classifier(identity_classifier()) model = meanclassif.train(gfeatures, glabels) assert model.apply(gfeatures[0]) == np.arange(10).mean() assert model.apply(gfeatures[1]) == .5 milk-0.3.10/milk/tests/fast_classifier.py0000644000175000017500000000126411464064572021152 0ustar luispedroluispedroimport numpy as np class fast_classifier(object): def __init__(self): pass def train(self, features, labels): examples = {} for f,lab in zip(features, labels): if lab not in examples: examples[lab] = f return fast_model(examples) class fast_model(object): def __init__(self, examples): self.examples = examples assert len(self.examples) def apply(self, f): best = None best_val = +np.inf for k,v in self.examples.iteritems(): dist = np.dot(v-f, v-f) if dist < best_val: best = k best_val = dist return best milk-0.3.10/milk/tests/__init__.py0000644000175000017500000000066411537727335017557 0ustar luispedroluispedrotry: import milksets del milksets except ImportError: import sys sys.stderr.write('''\ Could not import milksets. This companion package does not provide any functionality, but is necessary for some of the testing.''') def run(): import nose from os import path currentdir = path.dirname(__file__) updir = path.join(currentdir, '..') nose.run('milk', argv=['', '--exe', '-w', updir]) milk-0.3.10/milk/tests/test_ext_jugparallel.py0000644000175000017500000000153611557532105022227 0ustar luispedroluispedrotry: import jug from jug import value import jug.options from jug.tests.task_reset import task_reset except ImportError: from nose import SkipTest def task_reset(f): def g(): raise SkipTest() return g @task_reset def test_nfoldcrossvalidation(): store, space = jug.jug.init('milk/tests/data/jugparallel_jugfile.py', 'dict_store') options = jug.options.default_options jug.jug.execute(store, options) assert len(jug.value(space['classified'])) == 2 assert len(jug.value(space['classified_wpred'])) ==3 @task_reset def test_kmeans(): from jug.task import alltasks store, space = jug.jug.init('milk/tests/data/jugparallel_kmeans_jugfile.py', 'dict_store') options = jug.options.default_options jug.jug.execute(store, options) assert len(value(space['clustered'])) == 2 milk-0.3.10/milk/tests/test_multi.py0000644000175000017500000000404111453112412020157 0ustar luispedroluispedroimport numpy as np import random import milk.supervised.svm import milk.supervised.multi from milk.supervised.classifier import ctransforms from fast_classifier import fast_classifier import milksets.wine features,labels = milksets.wine.load() A = np.arange(len(features)) random.seed(9876543210) random.shuffle(A) features = features[A] labels = labels[A] labelset = set(labels) base = ctransforms(milk.supervised.svm.svm_raw(C=2.,kernel=milk.supervised.svm.rbf_kernel(2.**-3)),milk.supervised.svm.svm_binary()) def test_one_against_rest(): M = milk.supervised.multi.one_against_rest(base) M = M.train(features[:100,:],labels[:100]) tlabels = [M.apply(f) for f in features[100:]] for tl in tlabels: assert tl in labelset def test_one_against_one(): M = milk.supervised.multi.one_against_one(base) M = M.train(features[:100,:],labels[:100]) tlabels = [M.apply(f) for f in features[100:]] for tl in tlabels: assert tl in labelset def test_two_thirds(): np.random.seed(2345) C = milk.supervised.defaultclassifier('fast') X = np.random.rand(120,4) X[:40] += np.random.rand(40,4) X[:40] += np.random.rand(40,4) X[40:80] -= np.random.rand(40,4) X[40:80] -= np.random.rand(40,4) Y = np.repeat(np.arange(3), 40) model = C.train(X,Y) Y_ = np.array([model.apply(x) for x in X]) assert (Y_ == Y).mean() * 3 > 2 def test_multi_labels(): clabels = [[lab, lab+7] for lab in labels] multi_label = milk.supervised.multi.one_against_rest_multi(base) model = multi_label.train(features[::2], clabels[::2]) test_vals = [model.apply(f) for f in features[1::2]] for ts in test_vals: if 0.0 in ts: assert 7.0 in ts if 1.0 in ts: assert 8.0 in ts if 2.0 in ts: assert 9.0 in ts def test_classifier_no_set_options(): # Basically these should not raise an exception milk.supervised.multi.one_against_rest_multi(fast_classifier()) milk.supervised.multi.one_against_rest(fast_classifier()) milk.supervised.multi.one_against_one(fast_classifier()) milk-0.3.10/milk/tests/test_utils.py0000644000175000017500000000176011452411266020202 0ustar luispedroluispedroimport numpy as np from milk.utils.utils import get_nprandom, get_pyrandom def test_nprandom(): assert get_nprandom(None).rand() != get_nprandom(None).rand() assert get_nprandom(1).rand() != get_nprandom(2).rand() assert get_nprandom(1).rand() == get_nprandom(1).rand() r = get_nprandom(1) assert get_nprandom(r).rand() != r.rand() def test_pyrandom(): assert get_pyrandom(None).random() != get_pyrandom(None).random() assert get_pyrandom(1).random() != get_pyrandom(2).random() assert get_pyrandom(1).random() == get_pyrandom(1).random() r = get_pyrandom(1) assert get_pyrandom(r).random() != r.random() def test_cross_random(): assert get_pyrandom(get_nprandom(1)).random() == get_pyrandom(get_nprandom(1)).random() assert get_nprandom(get_pyrandom(1)).rand() == get_nprandom(get_pyrandom(1)).rand() def test_recursive(): def recurse(f): R = f(None) assert f(R) is R yield recurse, get_pyrandom yield recurse, get_nprandom milk-0.3.10/milk/tests/test_defaultclassifier.py0000644000175000017500000000141011452411266022523 0ustar luispedroluispedroimport numpy as np import milk.supervised.defaultclassifier import pickle def test_defaultclassifier(): from milksets import wine features, labels = wine.load() C = milk.supervised.defaultclassifier() model = C.train(features,labels) labelset = set(labels) for f in features: assert model.apply(f) in labelset test_defaultclassifier.slow = True def test_pickle(): np.random.seed(23232432) X = np.random.rand(100,10) labels = np.zeros(100) X[50:] += .5 labels[50:] = 1 classifier = milk.supervised.defaultclassifier() model = classifier.train(X, labels) s = pickle.dumps(model) model = pickle.loads(s) test = [model.apply(x) for x in X] test = np.array(test) assert (test == labels).mean() > .6 milk-0.3.10/milk/tests/test_basic.py0000644000175000017500000000004311256467633020127 0ustar luispedroluispedrodef test_import(): import milk milk-0.3.10/milk/tests/test_tree.py0000644000175000017500000000447711537456345020024 0ustar luispedroluispedroimport milk.supervised.tree import milk.supervised._tree from milk.supervised._tree import set_entropy from milk.supervised.tree import information_gain import numpy as np def test_tree(): from milksets import wine features, labels = wine.load() selected = (labels < 2) features = features[selected] labels = labels[selected] C = milk.supervised.tree.tree_classifier() model = C.train(features,labels) assert (np.array([model.apply(f) for f in features]) == labels).mean() > .5 def test_split_subsample(): import random from milksets import wine features, labels = wine.load() labels = labels.astype(np.int) seen = set() for i in xrange(20): random.seed(2) i,s = milk.supervised.tree._split(features[::10], labels[::10], None, milk.supervised.tree.information_gain, 2, random) seen.add(i) assert len(seen) <= 2 def test_set_entropy(): labels = np.arange(101)%3 counts = np.zeros(3) entropy = milk.supervised._tree.set_entropy(labels, counts) slow_counts = np.array([(labels == i).sum() for i in xrange(3)]) assert np.all(counts == slow_counts) px = slow_counts.astype(float)/ slow_counts.sum() slow_entropy = - np.sum(px * np.log(px)) assert np.abs(slow_entropy - entropy) < 1.e-8 def slow_information_gain(labels0, labels1): H = 0. N = len(labels0) + len(labels1) nlabels = 1+max(labels0.max(), labels1.max()) counts = np.empty(nlabels, np.double) for arg in (labels0, labels1): H -= len(arg)/float(N) * set_entropy(arg, counts) return H def test_information_gain(): np.random.seed(22) for i in xrange(8): labels0 = (np.random.randn(20) > .2).astype(int) labels1 = (np.random.randn(33) > .8).astype(int) fast = information_gain(labels0, labels1) slow = slow_information_gain(labels0, labels1) assert np.abs(fast - slow) < 1.e-8 def test_information_gain_small(): labels1 = np.array([0]) labels0 = np.array([0, 1]) assert information_gain(labels0, labels1) < 0. def test_z1_loss(): from milk.supervised.tree import z1_loss L0 = np.zeros(10) L1 = np.ones(10) L1[3] = 0 W0 = np.ones(10) W1 = np.ones(10) assert z1_loss(L0, L1) == z1_loss(L0, L1, W0, W1) assert z1_loss(L0, L1) != z1_loss(L0, L1, W0, .8*W1) milk-0.3.10/milk/tests/test_som.py0000644000175000017500000000374311454607320017643 0ustar luispedroluispedroimport numpy as np import random from milk.unsupervised import som from milk.unsupervised.som import putpoints, closest def _slow_putpoints(grid, points, L=.2): for point in points: dpoint = grid-point y,x = np.unravel_index(np.abs(dpoint).argmin(), dpoint.shape) for dy in xrange(-4, +4): for dx in xrange(-4, +4): ny = y + dy nx = x + dx if ny < 0 or ny >= grid.shape[0]: continue if nx < 0 or nx >= grid.shape[1]: continue L2 = L/(1+np.abs(dy)+np.abs(dx)) grid[ny,nx] *= 1. - L2 grid[ny,nx] += point*L2 def data_grid(): np.random.seed(22) data = np.arange(100000, dtype=np.float32) grid = np.array([data.flat[np.random.randint(0, data.size)] for i in xrange(64*64)]).reshape((64,64,1)) data = data.reshape((-1,1)) return grid, data def test_putpoints(): grid, points = data_grid() points = points[:100] grid2 = grid.copy() putpoints(grid, points, L=0., R=1) assert np.all(grid == grid2) putpoints(grid, points, L=.5, R=1) assert not np.all(grid == grid2) def test_against_slow(): grid, points = data_grid() grid2 = grid.copy() putpoints(grid, points[:10], shuffle=False) _slow_putpoints(grid2.reshape((64,64)), points[:10]) assert np.allclose(grid, grid2) def test_som(): N = 10000 np.random.seed(2) data = np.array([np.arange(N), N/4.*np.random.randn(N)]) data = data.transpose().copy() grid = som(data, (8,8), iterations=3, R=4) assert grid.shape == (8,8,2) y,x = closest(grid, data[0]) assert 0 <= y < grid.shape[0] assert 0 <= x < grid.shape[1] grid2 = grid.copy() np.random.shuffle(grid2) full = np.abs(np.diff(grid2[:,:,0], axis=0)).mean() obs = np.abs(np.diff(grid[:,:,0], axis=0)).mean() obs2 = np.abs(np.diff(grid[:,:,0], axis=1)).mean() assert obs + 4*np.abs(obs-obs2) < full milk-0.3.10/milk/tests/test_multi_view.py0000644000175000017500000000216611452411266021227 0ustar luispedroluispedroimport milk.supervised.multi_view import numpy as np import milk.supervised.svm from milk.supervised.defaultclassifier import feature_selection_simple def test_multi_view(): from milksets.wine import load features, labels = load() features0 = features[::10] features1 = features[1::10] features2 = features[2::10] labels0 = labels[::10] labels1 = labels[1::10] labels2 = labels[2::10] assert np.all(labels0 == labels1) assert np.all(labels1 == labels2) labels = labels0 train_features = zip(features0,features1,features2) test_features = zip(features[3::10], features[4::10], features[5::10]) base = milk.supervised.classifier.ctransforms( feature_selection_simple(), milk.supervised.svm.svm_raw(C=128, kernel=milk.supervised.svm.rbf_kernel(4.)), milk.supervised.svm.svm_sigmoidal_correction() ) classifier = milk.supervised.multi_view.multi_view_classifier([base,base,base]) model = classifier.train(train_features, labels == 0) assert ([model.apply(f) for f in test_features] == (labels == 0)).mean() > .9 milk-0.3.10/milk/tests/test_nfoldcrossvalidation_regression.py0000644000175000017500000002350711524577615025547 0ustar luispedroluispedroimport numpy as np from milk.measures.nfoldcrossvalidation import nfoldcrossvalidation, foldgenerator # Regression test in 2011-01-31 def test_getfoldgenerator(): labels = np.array([ 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ]) origins = np.array([ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 10, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 32, 32, 32, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 49, 49, 49, 49, 49, 49, 49, 49, 49, 49, 49, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 54, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 55, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56 ]) for Tr,Te in foldgenerator(labels, 3, origins): assert (np.array(labels)[Te] == 2).any() for Tr,Te in foldgenerator(labels[::2], 3, origins[::2]): assert (np.array(labels[::2])[Te] == 2).any() for Tr,Te in foldgenerator(labels[::3], 3, origins[::3]): assert (np.array(labels[::3])[Te] == 2).any() def test_getfoldgenerator_simplified(): # This is a cut-down version of the above labels = np.zeros(45, bool) labels[35:] = 1 origins = np.array([0, 0, 1, 2, 2, 2, 3, 4, 4, 4, 5, 6, 6, 6, 6, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 14, 14, 14, 14, 15, 15, 15, 16]) for Tr,Te in foldgenerator(labels, 3, origins): assert np.any(labels[Te]) for Tr,Te in foldgenerator(labels, 4, origins): assert np.any(labels[Te]) def test_getfoldgenerator_simplified_2(): # This is a cut-down version of the above labels = np.zeros(44, bool) labels[35:] = 1 origins = np.array([0, 0, 1, 2, 2, 2, 3, 4, 4, 4, 5, 6, 6, 6, 6, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 14, 14, 14, 14, 15, 15, 15]) for Tr,Te in foldgenerator(labels, 3, origins): assert np.any(labels[Te]) milk-0.3.10/milk/tests/test_gaussianmixture.py0000644000175000017500000000221711452411266022270 0ustar luispedroluispedroimport numpy as np from milk.unsupervised import gaussianmixture def _sq(x): return x*x def test_gm(): np.random.seed(22) centroids = np.repeat(np.arange(4), 4).reshape((4,4)) fmatrix = np.concatenate([(np.random.randn(12,4)+c) for c in centroids]) assignments = np.repeat(np.arange(4), 12) rss = sum(np.sum(_sq(fmatrix[i*12:(i+1)*12]-i)) for i in xrange(4)) assert np.abs(gaussianmixture.residual_sum_squares(fmatrix, assignments, centroids) - rss) < 1.e-12 assert gaussianmixture.BIC(fmatrix, assignments, centroids) > 0 assert gaussianmixture.AIC(fmatrix, assignments, centroids) > 0 assert gaussianmixture.BIC(fmatrix, assignments, centroids, model='full_covariance') > \ gaussianmixture.BIC(fmatrix, assignments, centroids, model='diagonal_covariance') > \ gaussianmixture.BIC(fmatrix, assignments, centroids, model='one_variance') assert gaussianmixture.AIC(fmatrix, assignments, centroids, model='full_covariance') > \ gaussianmixture.AIC(fmatrix, assignments, centroids, model='diagonal_covariance') > \ gaussianmixture.AIC(fmatrix, assignments, centroids, model='one_variance') milk-0.3.10/milk/tests/test_curves.py0000644000175000017500000000110611561061601020336 0ustar luispedroluispedrofrom milk.measures.curves import precision_recall import numpy as np def test_precision_recall(): labels = [0,1]*10 values = np.linspace(0,1,len(labels)) precision, recall = precision_recall(values, labels) assert np.min(recall) >= 0. assert np.max(recall) <= 1. assert np.max(precision) <= 1. assert np.min(precision) >= 0. labels = [0]*10 + [1] * 10 values = np.linspace(0,1.,20) precision,recall = precision_recall(values, labels, 'steps', 10) assert min(precision) >= .5 assert max(precision) == 1. assert max(recall) == 1. milk-0.3.10/milk/tests/test_nnmf.py0000644000175000017500000000302311452411266017772 0ustar luispedroluispedroimport milk.unsupervised import numpy as np def test_nnmf(): def test3(method): np.random.seed(8) X3 = np.random.rand(20,3) X = np.c_[ X3, X3*2+np.random.rand(20,3)/20., -X3*2+np.random.rand(20,3)/10.] W,V = method(X, 3, R=7) assert np.sum((np.dot(W,V)-X)**2)/np.sum(X**2) < .5 yield test3, milk.unsupervised.lee_seung yield test3, milk.unsupervised.sparse_nnmf def test_sparse_nnmf(): # This is really just a smoke test because the test case is not sparse!! from milk.unsupervised import sparse_nnmf np.random.seed(8) X3 = np.random.rand(20,3) X = np.c_[ X3, X3*2+np.random.rand(20,3)/20., -X3*2+np.random.rand(20,3)/10.] W,V = sparse_nnmf(X, 3, sparsenessW=.7, sparsenessH=.7, R=7) assert not np.any(np.isnan(W)) assert not np.any(np.isnan(V)) error = np.dot(W,V)-X assert error.var() < X.var() def test_hoyer_project(): from milk.unsupervised.nnmf.hoyer import _L1for, _project def sp(n, L1, L2): return (np.sqrt(n) - L1/L2)/(np.sqrt(n) - 1) sparseness = .6 n = 9. row = np.arange(int(n))/n L2 = np.sqrt(np.dot(row, row)) L1 = _L1for(sparseness, row, L2) assert np.abs(sp(n, L1, L2) - sparseness) < 1.e-4 row_ = _project(row, L1, L2) assert not np.any(np.isnan(row_)) assert np.all(row_ >= 0) L2 = np.sqrt(np.dot(row, row)) L1 = np.sum(np.abs(row_)) res = sp(n, L1, L2) assert np.abs(res - sparseness) < 1.e-4 milk-0.3.10/milk/tests/test_set2binary_array.py0000644000175000017500000000074111310124730022306 0ustar luispedroluispedroimport numpy as np from milk.supervised import set2binary_array def test_set2binary_array_len(): s2f = set2binary_array.set2binary_array() inputs = [ np.arange(1,3)*2, np.arange(4)**2, np.arange(6)+2 ] labels = [0,0,1] model = s2f.train(inputs,labels) assert len(model.apply(inputs[0])) == len(model.apply(inputs[1])) assert len(model.apply(inputs[0])) == len(model.apply(inputs[2])) assert len(model.apply(inputs[0])) == len(model.apply(range(128))) milk-0.3.10/milk/tests/test_measures.py0000644000175000017500000000200011464064572020661 0ustar luispedroluispedroimport milk.measures.measures import numpy as np import numpy from milk.measures import accuracy, waccuracy def test_100(): C=numpy.zeros((2,2)) C[0,0]=100 C[1,1]=50 assert accuracy(C) == 1. assert waccuracy(C) == 1. def test_0(): C = numpy.array([ [0, 10], [10, 0] ]) assert waccuracy(C) == 0. assert accuracy(C) == 0. def test_50(): C = numpy.array([ [10, 10], [10, 10] ]) assert accuracy(C) == .5 assert waccuracy(C) == .5 def test_unbalanced(): C = numpy.array([ [20, 10], [10, 0] ]) assert accuracy(C) == .5 assert waccuracy(C) == 1./3 def test_confusion_matrix(): np.random.seed(323) labels0 = np.arange(101)%3 labels1 = (labels0 + np.random.rand(101)*2).astype(np.int) % 3 cmat = milk.measures.measures.confusion_matrix(labels0, labels1) for i in xrange(3): for j in xrange(3): assert cmat[i,j] == np.sum( (labels0 == i) & (labels1 == j) ) milk-0.3.10/milk/tests/test_fisher.py0000644000175000017500000000233011452411266020314 0ustar luispedroluispedroimport milk.supervised.svm import milk.supervised.normalise import numpy as np import milk.supervised.svm def _slow_f(features,labels,kernel_or_sigma): try: kernel = kernel_or_sigma kernel(features[0],features[1]) except: kernel = milk.supervised.svm.rbf_kernel(kernel_or_sigma) N1 = (labels == 0).sum() N2 = (labels == 1).sum() x1 = features[labels == 0] x2 = features[labels == 1] dm = 0 for i in xrange(N1): for j in xrange(N1): dm += kernel(x1[i],x1[j])/N1/N1 for i in xrange(N2): for j in xrange(N2): dm += kernel(x2[i],x2[j])/N2/N2 for i in xrange(N1): for j in xrange(N2): dm -= 2*kernel(x1[i],x2[j])/N1/N2 s1 = N1 for i in xrange(N1): for j in xrange(N1): s1 -= kernel(x1[i],x1[j])/N1 s2 = N2 for i in xrange(N2): for j in xrange(N2): s2 -= kernel(x2[i],x2[j])/N2 return (s1 + s2)/dm def test_fisher_approx(): from milksets import wine features,labels = wine.load() f = milk.supervised.svm.sigma_value_fisher(features,labels) for sigma in (2.**-4,2.,16.,32.): assert abs(f(sigma) - _slow_f(features,labels,sigma)) < 1e-6 milk-0.3.10/milk/utils/0000755000175000017500000000000011562337213015424 5ustar luispedroluispedromilk-0.3.10/milk/utils/utils.py0000644000175000017500000000534511452411266017144 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2010, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. from __future__ import division import numpy as np import random __all__ = [ 'get_nprandom', 'get_pyrandom', ] def get_nprandom(R): ''' R' = get_nprandom(R) Returns a numpy.RandomState from R Parameters ---------- R : can be one of: None : Returns the default numpy global state integer : Uses it as a seed for constructing a new random generator RandomState : returns R Returns ------- R' : np.RandomState ''' if R is None: return np.random.mtrand._rand if type(R) == int: return np.random.RandomState(R) if type(R) is random.Random: return np.random.RandomState(R.randint(0, 2**30)) if type(R) is np.random.RandomState: return R raise TypeError,"get_nprandom() does not know how to handle type %s." % type(R) def get_pyrandom(R): ''' R = get_pyrandom(R) Returns a random.Random object based on R Parameters ---------- R : can be one of: None : Returns the default numpy global state integer : Uses it as a seed for constructing a new random generator RandomState : returns R Returns ------- R' : random.Random ''' if R is None: return random.seed.im_self if type(R) is int: return random.Random(R) if type(R) is np.random.RandomState: return random.Random(R.randint(2**30)) if type(R) is random.Random: return R raise TypeError,"get_pyrandom() does not know how to handle type %s." % type(R) milk-0.3.10/milk/utils/__init__.py0000644000175000017500000000002411231643006017522 0ustar luispedroluispedrofrom utils import * milk-0.3.10/milk/unsupervised/0000755000175000017500000000000011562337213017020 5ustar luispedroluispedromilk-0.3.10/milk/unsupervised/parzen.py0000644000175000017500000000551311561061716020676 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2010, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. from __future__ import division import numpy as np def get_parzen_1class_rbf_loocv(features): ''' f,fprime = get_parzen_1class_rbf_loocv(features) Leave-one-out crossvalidation value for 1-class Parzen window evaluator for features. Parameters ---------- features : ndarray feature matrix Returns ------- f : function: double -> double function which evaluates the value of a window value. Minize to get the best window value. fprime : function: double -> double function: df/dh ''' from milk.unsupervised.pdist import pdist D2 = -pdist(features) n = len(features) sumD2 = D2.sum() D2.flat[::(n+1)] = -np.inf def f(h): D2h = D2 / (2.*h) np.exp(D2h, D2h) val = D2h.sum() return val/np.sqrt(2*h*np.pi) def fprime(h): D2h = D2 / (2.*h) D2h.flat[::(n+1)] = 1. D2h *= np.exp(D2h) val = D2h.sum() - D2h.trace() val /= np.sqrt(2*h*np.pi) return -1./(4*np.pi*h)*f(h) + val return f,fprime def parzen(features, h): ''' f = parzen(features, h) Parzen window smoothing Parameters ---------- features : ndarray feature matrix h : double bandwidth Returns ------- f : callable (double^N -> double) density function ''' sum2 = np.array([np.dot(f,f) for f in features]) N = len(features) beta = np.sqrt(2*h*np.pi)/N def f(x): dist = np.dot(features, -2*x) dist += sum2 dist += np.dot(c,c) dist /= 2.*h np.exp(dist, dist) val = dist.sum() return val*beta return f milk-0.3.10/milk/unsupervised/kmeans.py0000644000175000017500000002531211557532105020654 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. from __future__ import division import numpy as np from numpy import linalg from . import _kmeans from ..utils import get_pyrandom from .normalise import zscore __all__ = [ 'kmeans', 'repeated_kmeans', ] try: # This tests for the presence of 3-argument np.dot # with the 3rd argument being the output argument _x = np.array([ [1., 1.], [0.,1.]] ) _y = np.array([2., 4.]) _r = np.array([0.,0.]) np.dot(_x, _y, _r) if _r[0] != 6 or _r[1] != 4: raise NotImplementedError _dot3 = np.dot except: def _dot3(x, y, _): return np.dot(x,y) finally: del _x del _y del _r def _mahalanobis2(fmatrix, x, icov): diff = fmatrix-x # The expression below seems to be faster than looping over the elements and summing return np.dot(diff, np.dot(icov, diff.T)).diagonal() def centroid_errors(fmatrix, assignments, centroids): ''' errors = centroid_errors(fmatrix, assignments, centroids) Computes the following:: for all i, j: ci = assignments[i] errors[i,j] = fmatrix[ci, j] - centroids[ci, j] Parameters ---------- fmatrix : 2D ndarray feature matrix assignments : 1D ndarray Assignments array centroids : 2D ndarray centroids Returns ------- errors : float Difference between fmatrix and corresponding centroid ''' errors = [] for k,c in enumerate(centroids): errors.append(fmatrix[assignments == k] - c) return np.concatenate(errors) def residual_sum_squares(fmatrix,assignments,centroids,distance='euclidean',**kwargs): ''' rss = residual_sum_squares(fmatrix, assignments, centroids, distance='euclidean', **kwargs) Computes residual sum squares Parameters ---------- fmatrix : 2D ndarray feature matrix assignments : 1D ndarray Assignments array centroids : 2D ndarray centroids Returns ------- rss : float residual sum squares ''' if distance != 'euclidean': raise NotImplemented, "residual_sum_squares only implemented for 'euclidean' distance" rss = 0.0 for k, c in enumerate(centroids): diff = fmatrix[assignments == k] - c diff = diff.ravel() rss += np.dot(diff, diff) return rss def assign_centroids(fmatrix, centroids, histogram=False, normalize=False, normalise=None): ''' cids = assign_centroids(fmatrix, centroids, histogram=False, normalize=False) Assigns a centroid to each element of fmatrix Parameters ---------- fmatrix : 2D ndarray feature matrix centroids : 2D ndarray centroids matrix histogram : boolean, optional If True, then the result is actually a histogram normalize : boolean, optional If True and ``histogram``, then the histogram is normalized to sum to one. Returns ------- cids : sequence ``cids[i]`` is the index of the centroid closes to ``fmatrix[i];`` or, if ``histogram``, then ``cids[i]`` is the number of points that were assigned to centroid ``i.`` ''' dists = np.dot(fmatrix, (-2)*centroids.T) dists += np.array([np.dot(c,c) for c in centroids]) cids = dists.argmin(1) if histogram: hist = np.array( [np.sum(cids == ci) for ci in xrange(len(centroids))], np.float) if (normalize or normalise) and len(fmatrix): hist /= hist.sum() return hist return cids def _pycomputecentroids(fmatrix, centroids, assignments, counts): k, Nf = centroids.shape bins = np.arange(k+1) ncounts,_ = np.histogram(assignments, bins) counts[:] = ncounts any_empty = False mean = None for ci,count in enumerate(counts): if count: where = (assignments.T == ci) mean = _dot3(where, fmatrix, mean) # mean = dot(fmatrix.T, where.T), but it is better to not cause copies mean /= count centroids[ci] = mean else: any_empty = True return any_empty def kmeans(fmatrix, k, distance='euclidean', max_iter=1000, R=None, **kwargs): ''' assignments, centroids = kmean(fmatrix, k, distance='euclidean', max_iter=1000, R=None, icov=None, covmat=None) k-Means Clustering Parameters ---------- fmatrix : ndarray 2-ndarray (Nelements x Nfeatures) distance: string, optional one of: - 'euclidean' : euclidean distance (default) - 'seuclidean' : standartised euclidean distance. This is equivalent to first normalising the features. - 'mahalanobis' : mahalanobis distance. This can make use of the following keyword arguments: + 'icov' (the inverse of the covariance matrix), + 'covmat' (the covariance matrix) If neither is passed, then the function computes the covariance from the feature matrix max_iter : integer, optional Maximum number of iteration (default: 1000) R : source of randomness, optional Returns ------- assignments : ndarray An 1-D array of size `len(fmatrix)` centroids : ndarray An array of `k'` centroids ''' fmatrix = np.asanyarray(fmatrix) if distance == 'seuclidean': fmatrix = zscore(fmatrix) distance = 'euclidean' if distance == 'euclidean': def distfunction(fmatrix, cs, dists): dists = _dot3(fmatrix, (-2)*cs.T, dists) dists += np.array([np.dot(c,c) for c in cs]) # For a distance, we'd need to add the fmatrix**2 components, but # it doesn't matter because we are going to perform an argmin() on # the result. return dists elif distance == 'mahalanobis': icov = kwargs.get('icov', None) if icov is None: covmat = kwargs.get('covmat', None) if covmat is None: covmat = np.cov(fmatrix.T) icov = linalg.inv(covmat) def distfunction(fmatrix, cs, _): return np.array([_mahalanobis2(fmatrix, c, icov) for c in cs]).T else: raise ValueError('milk.unsupervised.kmeans: `distance` argument unknown (%s)' % distance) if k < 2: raise ValueError('milk.unsupervised.kmeans `k` should be >= 2.') if fmatrix.dtype in (np.float32, np.float64) and fmatrix.flags['C_CONTIGUOUS']: computecentroids = _kmeans.computecentroids else: computecentroids = _pycomputecentroids R = get_pyrandom(R) centroids = np.array(R.sample(fmatrix,k), fmatrix.dtype) prev = np.zeros(len(fmatrix), np.int32) counts = np.empty(k, np.int32) dists = None for i in xrange(max_iter): dists = distfunction(fmatrix, centroids, dists) assignments = dists.argmin(1) if np.all(assignments == prev): break if computecentroids(fmatrix, centroids, assignments.astype(np.int32), counts): (empty,) = np.where(counts == 0) centroids = np.delete(centroids, empty, axis=0) k = len(centroids) counts = np.empty(k, np.int32) # This will cause new matrices to be allocated in the next iteration dists = None prev[:] = assignments return assignments, centroids def repeated_kmeans(fmatrix,k,iterations,distance='euclidean',max_iter=1000,R=None,**kwargs): ''' assignments,centroids = repeated_kmeans(fmatrix, k, repeats, distance='euclidean',max_iter=1000,**kwargs) Runs kmeans repeats times and returns the best result as evaluated according to distance See Also -------- kmeans : runs kmeans once Parameters ---------- fmatrix : feature matrix k : nr of centroids iterations : Nr of repetitions distance : 'euclidean' (default) or 'seuclidean' max_iter : Max nr of iterations per kmeans run R : random source Returns ------- assignments : 1-D array of assignments centroids : centroids These are the same returns as the kmeans function ''' kwargs['max_iter'] = max_iter return select_best_kmeans(fmatrix, [k], repeats=iterations, method='loglike', distance=distance, **kwargs) def select_best_kmeans(fmatrix, ks, repeats=1, method='AIC', R=None, **kwargs): ''' assignments,centroids = select_best_kmeans(fmatrix, ks, repeats=1, method='AIC', R=None, **kwargs) Runs kmeans repeats times and returns the best result as evaluated according to distance See Also -------- kmeans : runs kmeans once Parameters ---------- fmatrix : feature matrix ks : sequence of integers nr of centroids to try iterations : integer, optional Nr of repetitions for each value of k R : random source, optional Returns ------- assignments : 1-D array of assignments centroids : 2-D ndarray centroids These are the same returns as the kmeans function ''' best = None best_val = np.inf R = get_pyrandom(R) from milk.unsupervised.gaussianmixture import AIC, BIC if method == 'AIC': method = AIC elif method == 'BIC': method = BIC else: raise ValueError('milk.kmeans.select_best_kmeans: unknown method: %s' % method) if 'distance' in kwargs and kwargs['distance'] == 'seuclidean': fmatrix = zscore(fmatrix) for k in ks: for i in xrange(repeats): As,Cs = kmeans(fmatrix, k, R=R, **kwargs) value = method(fmatrix, As, Cs) if value < best_val: best_val = value best = As,Cs return best milk-0.3.10/milk/unsupervised/_som.cpp0000644000175000017500000000750211453677127020477 0ustar luispedroluispedro#include #include #include extern "C" { #include #include } namespace { struct SOM_Exception { SOM_Exception(const char* msg): msg(msg) { } const char* msg; }; void assert_type_contiguous(PyArrayObject* array,int type) { if (!PyArray_Check(array) || PyArray_TYPE(array) != type || !PyArray_ISCONTIGUOUS(array)) { throw SOM_Exception("Arguments to putpoints don't conform to expectation. Are you calling this directly? This is an internal function!"); } } void putpoints(PyArrayObject* grid, PyArrayObject* points, float L, int radius) { if (PyArray_NDIM(grid) != 3) throw SOM_Exception("grid should be three dimensional"); if (PyArray_NDIM(points) != 2) throw SOM_Exception("points should be two dimensional"); const int rows = PyArray_DIM(grid, 0); const int cols = PyArray_DIM(grid, 1); const int d = PyArray_DIM(grid, 2); const int n = PyArray_DIM(points, 0); if (PyArray_DIM(points, 1) != d) throw SOM_Exception("second dimension of points is not third dimension of grid"); Py_BEGIN_ALLOW_THREADS for (int i = 0; i != n; i++){ const float* p = static_cast(PyArray_GETPTR1(points, i)); int min_y = 0; int min_x = 0; float best = std::numeric_limits::max(); for (int y = 0; y != rows; ++y) { for (int x = 0; x != cols; ++x) { float dist = 0.; const float* gpoint = static_cast(PyArray_GETPTR2(grid, y, x)); for (int j = 0; j != d; ++j) { dist += (p[j] - gpoint[j])*(p[j] - gpoint[j]); } if (dist < best) { best = dist; min_y = y; min_x = x; } } } const int start_y = std::max(0, min_y - radius); const int start_x = std::max(0, min_x - radius); const int end_y = std::min(rows, min_y + radius); const int end_x = std::min(rows, min_x + radius); for (int y = start_y; y != end_y; ++y) { for (int x = start_x; x != end_x; ++x) { const float L2 = L /(1 + std::abs(min_y - y) + std::abs(min_x - x)); float* gpoint = static_cast(PyArray_GETPTR2(grid, y, x)); for (int j = 0; j != d; ++j) { gpoint[j] *= (1.-L2); gpoint[j] += L2 * p[j]; } } } } Py_END_ALLOW_THREADS } PyObject* py_putpoints(PyObject* self, PyObject* args) { try { PyArrayObject* grid; PyArrayObject* points; float L; int radius; if (!PyArg_ParseTuple(args, "OOfi", &grid, &points, &L, &radius)) { const char* errmsg = "Arguments were not what was expected for putpoints.\n" "This is an internal function: Do not call directly unless you know exactly what you're doing.\n"; PyErr_SetString(PyExc_RuntimeError,errmsg); return 0; } assert_type_contiguous(grid, NPY_FLOAT); assert_type_contiguous(points, NPY_FLOAT); putpoints(grid, points, L, radius); Py_RETURN_NONE; } catch (const SOM_Exception& exc) { PyErr_SetString(PyExc_RuntimeError,exc.msg); return 0; } catch (...) { PyErr_SetString(PyExc_RuntimeError,"Some sort of exception in putpoints."); return 0; } } PyMethodDef methods[] = { {"putpoints", py_putpoints, METH_VARARGS , "Do NOT call directly.\n" }, {NULL, NULL,0,NULL}, }; const char * module_doc = "Internal SOM Module.\n" "\n" "Do NOT use directly!\n"; } // namespace extern "C" void init_som() { import_array(); (void)Py_InitModule3("_som", methods, module_doc); } milk-0.3.10/milk/unsupervised/gaussianmixture.py0000644000175000017500000001236711525045463022635 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division import numpy as np from numpy import log, pi, array from numpy.linalg import det, inv from kmeans import residual_sum_squares, centroid_errors __all__ = [ 'BIC', 'AIC', 'log_likelihood', 'nr_parameters', ] def log_likelihood(fmatrix,assignments,centroids,model='one_variance',covs=None): ''' log_like = log_likelihood(feature_matrix, assignments, centroids, model='one_variance', covs=None) Compute the log likelihood of feature_matrix[i] being generated from centroid[i] ''' N,q = fmatrix.shape k = len(centroids) if model == 'one_variance': Rss = residual_sum_squares(fmatrix,assignments,centroids) #sigma2=Rss/N return -N/2.*log(2*pi*Rss/N)-N/2 elif model == 'diagonal_covariance': errors = centroid_errors(fmatrix,assignments,centroids) errors *= errors errors = errors.sum(1) Rss = np.zeros(k) counts = np.zeros(k) for i in xrange(fmatrix.shape[0]): c = assignments[i] Rss[c] += errors[i] counts[c] += 1 sigma2s = Rss/(counts+(counts==0)) return -N/2.*log(2*pi) -N/2. -1/2.*np.sum(counts*np.log(sigma2s+(counts==0))) elif model == 'full_covariance': res = -N*q/2.*log(2*pi) for k in xrange(len(centroids)): diff = (fmatrix[assignments == k] - centroids[k]) if covs is None: covm = np.cov(diff.T) else: covm = covs[k] if covm.shape == (): covm = np.matrix([[covm]]) icov = np.matrix(inv(covm)) diff = np.matrix(diff) Nk = diff.shape[0] res += -Nk/2.*log(det(covm)) + \ -.5 * (diff * icov * diff.T).diagonal().sum() return res raise ValueError, "log_likelihood: cannot handle model '%s'" % model def nr_parameters(fmatrix,k,model='one_variance'): ''' nr_p = nr_parameters(fmatrix, k, model='one_variance') Compute the number of parameters for a model of k clusters on Parameters ---------- fmatrix : 2d-array feature matrix k : integer nr of clusters model : str one of 'one_variance' (default), 'diagonal_covariance', or 'full_covariance' Returns ------- nr_p : integer Number of parameters ''' N,q = fmatrix.shape if model == 'one_variance': return k*q+1 elif model == 'diagonal_covariance': return k*(q+1) elif model == 'full_covariance': return k*+q*q raise ValueError, "milk.unsupervised.gaussianmixture.nr_parameters: cannot handle model '%s'" % model def _compute(type, fmatrix, assignments, centroids, model='one_variance', covs=None): N,q = fmatrix.shape k = len(centroids) log_like = log_likelihood(fmatrix, assignments, centroids, model, covs) n_param = nr_parameters(fmatrix,k,model) if type == 'BIC': return -2*log_like + n_param * log(N) elif type == 'AIC': return -2*log_like + 2 * n_param else: assert False def BIC(fmatrix, assignments, centroids, model='one_variance', covs=None): ''' B = BIC(fmatrix, assignments, centroids, model='one_variance', covs={From Data}) Compute Bayesian Information Criterion Parameters ---------- fmatrix : 2d-array feature matrix assignments : 2d-array Centroid assignments centroids : sequence Centroids model : str, optional one of 'one_variance' All features share the same variance parameter sigma^2. Default 'full_covariance' Estimate a full covariance matrix or use covs[i] for centroid[i] covs : sequence or matrix, optional Covariance matrices. If None, then estimate from the data. If scalars instead of matrices are given, then s stands for sI (i.e., the diagonal matrix with s along the diagonal). Returns ------- B : float BIC value See Also -------- AIC ''' return _compute('BIC', fmatrix, assignments, centroids, model, covs) def AIC(fmatrix,assignments,centroids,model='one_variance',covs=None): ''' A = AIC(fmatrix,assignments,centroids,model) Compute Akaike Information Criterion Parameters ---------- fmatrix : 2d-array feature matrix assignments : 2d-array Centroid assignments centroids : sequence Centroids model : str, optional one of 'one_variance' All features share the same variance parameter sigma^2. Default 'full_covariance' Estimate a full covariance matrix or use covs[i] for centroid[i] covs : sequence, optional Covariance matrices. If None, then estimate from the data. If scalars instead of matrices are given, then s stands for sI (i.e., the diagonal matrix with s along the diagonal). Returns ------- B : float AIC value See Also -------- BIC ''' return _compute('AIC', fmatrix, assignments, centroids, model, covs) milk-0.3.10/milk/unsupervised/som.py0000644000175000017500000000627511454604522020202 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2010, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division import numpy as np from ..utils import get_pyrandom from . import _som def putpoints(grid, points, L=.2, radius=4, iterations=1, shuffle=True, R=None): ''' putpoints(grid, points, L=.2, radius=4, iterations=1, shuffle=True, R=None) Feeds elements of `points` into the SOM `grid` Parameters ---------- grid : ndarray Self organising map points : ndarray data to feed to array L : float, optional How much to influence neighbouring points (default: .2) radius : integer, optional Maximum radius of influence (in L_1 distance, default: 4) iterations : integer, optional Number of iterations shuffle : boolean, optional Whether to shuffle the points before each iterations R : source of randomness ''' if radius is None: radius = 4 if type(L) != float: raise TypeError("milk.unsupervised.som: L should be floating point") if type(radius) != int: raise TypeError("milk.unsupervised.som: radius should be an integer") if grid.dtype != np.float32: raise TypeError('milk.unsupervised.som: only float32 arrays are accepted') if points.dtype != np.float32: raise TypeError('milk.unsupervised.som: only float32 arrays are accepted') if len(grid.shape) == 2: grid = grid.reshape(grid.shape+(1,)) if shuffle: random = get_pyrandom(R) for i in xrange(iterations): if shuffle: random.shuffle(points) _som.putpoints(grid, points, L, radius) def closest(grid, f): ''' y,x = closest(grid, f) Finds the coordinates of the closest point in the `grid` to `f` :: y,x = \\argmin_{y,x} { || grid[y,x] - f ||^2 } Parameters ---------- grid : ndarray of shape Y,X,J self-organised map f : ndarray of shape J point Returns ------- y,x : integers coordinates into `grid` ''' delta = grid - f delta **= 2 delta = delta.sum(2) return np.unravel_index(delta.argmin(), delta.shape) def som(data, shape, iterations=1000, L=.2, radius=4, R=None): ''' grid = som(data, shape, iterations=1000, L=.2, radius=4, R=None): Self-organising maps Parameters ---------- points : ndarray data to feed to array shape : tuple Desired shape of output. Must be 2-dimensional. L : float, optional How much to influence neighbouring points (default: .2) radius : integer, optional Maximum radius of influence (in L_1 distance, default: 4) iterations : integer, optional Number of iterations R : source of randomness Returns ------- grid : ndarray Map ''' R = get_pyrandom(R) d = data.shape[1] if data.dtype != np.float32: data = data.astype(np.float32) grid = np.array(R.sample(data, np.product(shape))).reshape(shape + (d,)) putpoints(grid, data, L=L, radius=radius, iterations=iterations, shuffle=True, R=R) return grid milk-0.3.10/milk/unsupervised/__init__.py0000644000175000017500000000135411513401272021125 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT. See COPYING.MIT file in the milk distribution ''' milk.unsupervised Unsupervised Learning --------------------- - kmeans: This is a highly optimised implementation of kmeans - PCA: Simple implementation - Non-negative matrix factorisation: both direct and with sparsity constraints ''' from kmeans import kmeans,repeated_kmeans from gaussianmixture import * from pca import pca import nnmf from nnmf import * from pdist import pdist, plike from .som import som __all__ = [ 'kmeans', 'repeated_kmeans', 'pca', 'pdist', 'plike', 'som', ] + \ nnmf.__all__ milk-0.3.10/milk/unsupervised/pdist.py0000644000175000017500000000311211513401272020503 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division import numpy as np __all__ = [ 'pdist', 'plike', ] def pdist(X, Y=None, distance='euclidean2'): ''' D = pdist(X, Y={X}, distance='euclidean2') Compute distance matrix:: D[i,j] == np.sum( (X[i] - Y[j])**2 ) Parameters ---------- X : feature matrix Y : feature matrix (default: use `X`) distance : one of 'euclidean' or 'euclidean2' (default) Returns ------- D : matrix of doubles ''' # Use Dij = np.dot(Xi, Xi) + np.dot(Xj,Xj) - 2.*np.dot(Xi,Xj) if Y is None: D = np.dot(X, X.T) x2 = D.diagonal() y2 = x2 else: D = np.dot(X, Y.T) x2 = np.array([np.dot(x,x) for x in X]) y2 = np.array([np.dot(y,y) for y in Y]) D *= -2. D += x2[:,np.newaxis] D += y2 # Because of numerical imprecision, we might get negative numbers # (which cause problems down the road, .e.g., when doing the sqrt): np.maximum(D, 0, D) if distance == 'euclidean': np.sqrt(D, D) return D def plike(X, sigma2=None): ''' L = plike(X, sigma2={guess based on X}) Parameters ---------- X : feature matrix sigma2 : bandwidth Returns ------- L : likelihood matrix ''' L = pdist(X) if sigma2 is None: sigma2 = np.median(L) L /= -sigma2 np.exp(L, L) return L milk-0.3.10/milk/unsupervised/_kmeans.cpp0000644000175000017500000001151011537446747021156 0ustar luispedroluispedro#include extern "C" { #include #include } namespace { struct Kmeans_Exception { Kmeans_Exception(const char* msg): msg(msg) { } const char* msg; }; void assert_type_contiguous(PyArrayObject* array,int type) { if (!PyArray_Check(array) || PyArray_TYPE(array) != type || !PyArray_ISCONTIGUOUS(array)) { throw Kmeans_Exception("Arguments to kmeans don't conform to expectation. Are you calling this directly? This is an internal function!"); } } template int computecentroids(ftype* f, ftype* centroids, PyArrayObject* a_assignments, PyArrayObject* a_counts, const int N, const int Nf, const int k) { int zero_counts = 0; Py_BEGIN_ALLOW_THREADS const npy_int32* assignments = static_cast(PyArray_DATA(a_assignments)); npy_int32* counts = static_cast(PyArray_DATA(a_counts)); std::fill(counts, counts + k, 0); std::fill(centroids, centroids + k*Nf, 0.0); for (int i = 0; i != N; i++){ const int c = assignments[i]; if (c >= k) continue; // throw Kmeans_Exception("wrong value in assignments"); ftype* ck = centroids + Nf*c; for (int j = 0; j != Nf; ++j) { *ck++ += *f++; } ++counts[c]; } for (int i = 0; i != k; ++i) { ftype* ck = centroids + Nf*i; const ftype c = counts[i]; if (c == 0) { ++zero_counts; } else { for (int j = 0; j != Nf; ++j) { *ck++ /= c; } } } Py_END_ALLOW_THREADS return zero_counts; } PyObject* py_computecentroids(PyObject* self, PyObject* args) { try { PyArrayObject* fmatrix; PyArrayObject* centroids; PyArrayObject* assignments; PyArrayObject* counts; if (!PyArg_ParseTuple(args, "OOOO", &fmatrix, ¢roids, &assignments, &counts)) { throw Kmeans_Exception("Wrong number of arguments for computecentroids."); } if (!PyArray_Check(fmatrix) || !PyArray_ISCONTIGUOUS(fmatrix)) throw Kmeans_Exception("fmatrix not what was expected."); if (!PyArray_Check(centroids) || !PyArray_ISCONTIGUOUS(centroids)) throw Kmeans_Exception("centroids not what was expected."); if (!PyArray_Check(counts) || !PyArray_ISCONTIGUOUS(counts)) throw Kmeans_Exception("counts not what was expected."); if (!PyArray_Check(assignments) || !PyArray_ISCONTIGUOUS(assignments)) throw Kmeans_Exception("assignments not what was expected."); if (PyArray_TYPE(counts) != NPY_INT32) throw Kmeans_Exception("counts should be int32."); //if (PyArray_TYPE(assignments) != NPY_INT32) throw Kmeans_Exception("assignments should be int32."); if (PyArray_TYPE(fmatrix) != PyArray_TYPE(centroids)) throw Kmeans_Exception("centroids and fmatrix should have same type."); if (PyArray_NDIM(fmatrix) != 2) throw Kmeans_Exception("fmatrix should be two dimensional"); if (PyArray_NDIM(centroids) != 2) throw Kmeans_Exception("centroids should be two dimensional"); if (PyArray_NDIM(assignments) != 1) throw Kmeans_Exception("assignments should be two dimensional"); const int N = PyArray_DIM(fmatrix, 0); const int Nf = PyArray_DIM(fmatrix, 1); const int k = PyArray_DIM(centroids, 0); if (PyArray_DIM(centroids, 1) != Nf) throw Kmeans_Exception("centroids has wrong number of features."); if (PyArray_DIM(assignments, 0) != N) throw Kmeans_Exception("assignments has wrong size."); if (PyArray_DIM(counts, 0) != k) throw Kmeans_Exception("counts has wrong size."); switch (PyArray_TYPE(fmatrix)) { #define TRY_TYPE(code, ftype) \ case code: \ if (computecentroids( \ static_cast(PyArray_DATA(fmatrix)), \ static_cast(PyArray_DATA(centroids)), \ assignments, \ counts, \ N, Nf, k)) { \ Py_RETURN_TRUE; \ } \ Py_RETURN_FALSE; TRY_TYPE(NPY_FLOAT, float); TRY_TYPE(NPY_DOUBLE, double); } throw Kmeans_Exception("Cannot handle this type."); } catch (const Kmeans_Exception& exc) { PyErr_SetString(PyExc_RuntimeError,exc.msg); return 0; } catch (...) { PyErr_SetString(PyExc_RuntimeError,"Some sort of exception in computecentroids."); return 0; } } PyMethodDef methods[] = { {"computecentroids", py_computecentroids, METH_VARARGS , "Do NOT call directly.\n" }, {NULL, NULL,0,NULL}, }; const char * module_doc = "Internal _kmeans Module.\n" "\n" "Do NOT use directly!\n"; } // namespace extern "C" void init_kmeans() { import_array(); (void)Py_InitModule3("_kmeans", methods, module_doc); } milk-0.3.10/milk/unsupervised/nnmf/0000755000175000017500000000000011562337213017756 5ustar luispedroluispedromilk-0.3.10/milk/unsupervised/nnmf/lee_seung.py0000644000175000017500000000573111452411266022303 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2010, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. from __future__ import division import numpy as np from numpy import dot from ...utils import get_nprandom __all__ = ['nnmf'] def nnmf(V, r, cost='norm2', max_iter=int(1e4), tol=1e-8, R=None): ''' A,S = nnmf(X, r, cost='norm2', tol=1e-8, R=None) Implement Lee & Seung's algorithm Parameters ---------- V : 2-ndarray input matrix r : integer nr of latent features cost : one of: 'norm2' : minimise || X - AS ||_2 (default) 'i-div' : minimise D(X||AS), where D is I-divergence (generalisation of K-L divergence) max_iter : integer, optional maximum number of iterations (default: 10000) tol : double tolerance threshold for early exit (when the update factor is with tol of 1., the function exits) R : integer, optional random seed Returns ------- A : 2-ndarray S : 2-ndarray Reference --------- "Algorithms for Non-negative Matrix Factorization" by Daniel D Lee, Sebastian H Seung (available at http://citeseer.ist.psu.edu/lee01algorithms.html) ''' # Nomenclature in the function follows lee & seung, while outside nomenclature follows eps = 1e-8 n,m = V.shape R = get_nprandom(R) W = R.standard_normal((n,r))**2 H = R.standard_normal((r,m))**2 for i in xrange(max_iter): if cost == 'norm2': updateH = dot(W.T,V)/(dot(dot(W.T,W),H)+eps) H *= updateH updateW = dot(V,H.T)/(dot(W,dot(H,H.T))+eps) W *= updateW elif cost == 'i-div': raise NotImplementedError,'I-Div not implemented in lee_seung.nnmf' if True or (i % 10) == 0: max_update = max(updateW.max(),updateH.max()) if abs(1.-max_update) < tol: break return W,H milk-0.3.10/milk/unsupervised/nnmf/__init__.py0000644000175000017500000000246411452411266022074 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2010, Luis Pedro Coelho # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. import lee_seung as ls from hoyer import sparse_nnmf import hoyer lee_seung = ls.nnmf __all__ = ['lee_seung','sparse_nnmf'] # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.3.10/milk/unsupervised/nnmf/hoyer.py0000644000175000017500000001023111452411266021452 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2010, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. from __future__ import division import numpy as np from ...utils import get_nprandom __all__ = ['hoyer_sparse_nnmf'] def sp(s): L2 = np.sqrt(np.dot(s,s)) L1 = np.abs(s).sum() sn = np.sqrt(len(s)) return (sn-L1/L2)/(sn-1) def _solve_alpha(s,m,L2): sm = s-m s2 = np.dot(s,s) sm2 = np.dot(sm, sm) m2 = np.dot(m, m) dot = np.dot(m, sm) alpha = (-dot + np.sqrt(dot**2 - sm2*(m2-L2**2)))/sm2 return alpha def _project(x,L1,L2): ''' Implement projection onto sparse space ''' x = np.asanyarray(x) n = len(x) s = x + (L1 - x.sum())/n Z = np.zeros(n,bool) while True: m = (~Z) * L1/(n-Z.sum()) alpha = _solve_alpha(s,m,L2) s = m + alpha * (s - m) negs = (s < 0) if not negs.any(): return s Z |= negs s[Z] = 0 c = (s.sum() - L1)/(~Z).sum() s -= c*(~Z) def _L1for(s,x,L2): ''' Solve for L1 in s = [ sqrt(n) - L1/L2] / [sqrt(n) - 1] ''' sn = np.sqrt(len(x)) return L2*(s+sn-s*sn) def sparse_nnmf(V, r, sparsenessW=None, sparsenessH=None, max_iter=10000, R=None): ''' W,H = hoyer.sparse_nnmf(V, r, sparsenessW = None, sparsenessH = None, max_iter=10000, R=None) Implement sparse nonnegative matrix factorisation. Parameters ---------- V : 2-D matrix input feature matrix r : integer number of latent features sparsenessW : double, optional sparseness contraint on W (default: no sparsity contraint) sparsenessH : double, optional sparseness contraint on H (default: no sparsity contraint) max_iter : integer, optional maximum nr of iterations (default: 10000) R : integer, optional source of randomness Returns ------- W : 2-ndarray H : 2-ndarray Reference --------- "Non-negative Matrix Factorisation with Sparseness Constraints" by Patrik Hoyer in Journal of Machine Learning Research 5 (2004) 1457--1469 ''' n,m = V.shape R = get_nprandom(R) mu_W = .15 mu_H = .15 eps = 1e-8 W = R.standard_normal((n,r))**2 H = R.standard_normal((r,m))**2 def fix(X, sparseness): for i in xrange(r): row = X[i] L2 = np.sqrt(np.dot(row, row)) row /= L2 X[i] = _project(row, _L1for(sparseness, row, 1.), 1.) def fixW(): fix(W.T, sparsenessW) def fixH(): fix(H, sparsenessH) if sparsenessW is not None: fixW() if sparsenessH is not None: fixH() for i in xrange(max_iter): if sparsenessW is not None: W -= mu_W * np.dot(np.dot(W,H)-V,H.T) fixW() else: updateW = np.dot(V,H.T)/(np.dot(W,np.dot(H,H.T))+eps) W *= updateW if sparsenessH is not None: H -= mu_H * np.dot(W.T,np.dot(W,H)-V) fixH() else: updateH = np.dot(W.T,V)/(np.dot(np.dot(W.T,W),H)+eps) H *= updateH return W,H milk-0.3.10/milk/unsupervised/pca.py0000644000175000017500000000172611540467121020141 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2011, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division import numpy as np from numpy import linalg from . import normalise __all__ = [ 'pca', ] def pca(X, zscore=True): ''' Y,V = pca(X, zscore=True) Principal Component Analysis Performs principal component analysis. Returns transformed matrix and principal components Parameters ---------- X : 2-dimensional ndarray data matrix zscore : boolean, optional whether to normalise to zscores (default: True) Returns ------- Y : ndarray Transformed matrix (of same dimension as X) V : ndarray principal components ''' if zscore: X = normalise.zscore(X) C = np.cov(X.T) w,v = linalg.eig(C) Y = np.dot(v,X.T).T return Y,v milk-0.3.10/milk/unsupervised/affinity.py0000644000175000017500000001054711513401272021203 0ustar luispedroluispedro# -*- coding: utf-8 -*- # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # Copyright (C) 2010-2011, # Luis Pedro Coelho , # Alexandre Gramfort , # Gael Varoquaux # # License: MIT. See COPYING.MIT file in the milk distribution """Algorithms for clustering : Meanshift and Affinity propagation Original Authors (for scikits.learn): Alexandre Gramfort alexandre.gramfort@inria.fr Gael Varoquaux gael.varoquaux@normalesup.org I have made this implementation more careful about allocating intermediate arrays. """ import numpy as np __all__ = [ 'affinity_propagation', ] def affinity_propagation(S, p=None, convit=30, maxit=200, damping=0.5, copy=True, R=0): """Perform Affinity Propagation Clustering of data Parameters ---------- S : array [n_points, n_points] Matrix of similarities between points p : array [n_points,] or float, optional Preferences for each point damping : float, optional Damping factor copy : boolean, optional If copy is False, the affinity matrix is modified inplace by the algorithm, for memory efficiency R : source of randomness Returns ------- cluster_centers_indices : array [n_clusters] index of clusters centers labels : array [n_points] cluster labels for each point Notes ----- See examples/plot_affinity_propagation.py for an example. Reference: Brendan J. Frey and Delbert Dueck, "Clustering by Passing Messages Between Data Points", Science Feb. 2007 """ if copy: # Copy the affinity matrix to avoid modifying it inplace S = S.copy() n_points = S.shape[0] assert S.shape[0] == S.shape[1] if p is None: p = np.median(S) if damping < 0.5 or damping >= 1: raise ValueError('damping must be >= 0.5 and < 1') random_state = np.random.RandomState(R) # Place preferences on the diagonal of S S.flat[::(n_points+1)] = p A = np.zeros((n_points, n_points)) R = np.zeros((n_points, n_points)) # Initialize messages # Remove degeneracies noise = random_state.randn(n_points, n_points) typeinfo = np.finfo(S.dtype) noise *= typeinfo.tiny*100 S += noise del noise # Execute parallel affinity propagation updates e = np.zeros((n_points, convit)) ind = np.arange(n_points) for it in range(maxit): Aold = A.copy() Rold = R.copy() A += S I = np.argmax(A, axis=1) Y = A[ind, I]#np.max(A, axis=1) A[ind, I] = typeinfo.min Y2 = np.max(A, axis=1) R = S - Y[:, np.newaxis] R[ind, I[ind]] = S[ind, I] - Y2 Rold *= damping R *= (1-damping) R += Rold # Compute availabilities Rd = R.diagonal().copy() np.maximum(R, 0, R) R.flat[::n_points+1] = Rd A = np.sum(R, axis=0)[np.newaxis, :] - R dA = np.diag(A) A = np.minimum(A, 0) A.flat[::n_points+1] = dA Aold *= damping A *= (1-damping) A += Aold # Check for convergence E = (np.diag(A) + np.diag(R)) > 0 e[:, it % convit] = E K = np.sum(E, axis=0) if it >= convit: se = np.sum(e, axis=1); unconverged = np.sum((se == convit) + (se == 0)) != n_points if (not unconverged and (K>0)) or (it==maxit): print "Converged after %d iterations." % it break else: print "Did not converge" I = np.where(np.diag(A+R) > 0)[0] K = I.size # Identify exemplars if K > 0: c = np.argmax(S[:, I], axis=1) c[I] = np.arange(K) # Identify clusters # Refine the final set of exemplars and clusters and return results for k in range(K): ii = np.where(c==k)[0] j = np.argmax(np.sum(S[ii, ii], axis=0)) I[k] = ii[j] c = np.argmax(S[:, I], axis=1) c[I] = np.arange(K) labels = I[c] # Reduce labels to a sorted, gapless, list cluster_centers_indices = np.unique(labels) labels = np.searchsorted(cluster_centers_indices, labels) else: labels = np.empty((n_points, 1)) cluster_centers_indices = None labels.fill(np.nan) return cluster_centers_indices, labels milk-0.3.10/milk/unsupervised/normalise.py0000644000175000017500000000300211452411266021355 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2010, Luis Pedro Coelho # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. from __future__ import division import numpy as np __all__ = ['zscore'] def zscore(features): """ features = zscore(features) Returns a copy of features which has been normalised to zscores """ mu = features.mean(0) sigma = np.std(features,0) sigma[sigma == 0] = 1 return (features - mu) / sigma # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.3.10/milk/wrapper/0000755000175000017500000000000011562337213015744 5ustar luispedroluispedromilk-0.3.10/milk/wrapper/__init__.py0000644000175000017500000000000011452411266020042 0ustar luispedroluispedromilk-0.3.10/milk/wrapper/wraplibsvm.py0000644000175000017500000000456411452411266020514 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2010, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # License: MIT. See COPYING.MIT file in the milk distribution from __future__ import division from milk.supervised.classifier import normaliselabels try: from libsvm import svm as libsvm except ImportError: try: import svm as libsvm except ImportError: libsvm = None from tempfile import NamedTemporaryFile class libsvmModel(object): def __init__(self, model, names, output_probability): self.model = model self.names = names self.output_probability = output_probability def apply(self,feats): if self.output_probability: return self.model.predict_probability(feats) res = self.model.predict(feats) return self.names[int(res)] def __getstate__(self): # This is really really really hacky, but it works N = NamedTemporaryFile() self.model.save(N.name) S = N.read() return S,self.output_probability,self.names def __setstate__(self,state): if libsvm is None: raise RuntimeError('LibSVM Library not found. Cannot use this classifier.') S,self.output_probability,self.names = state N = NamedTemporaryFile() N.write(S) N.flush() self.model = libsvm.svm_model(N.name) class libsvmClassifier(object): def __init__(self,probability = False, auto_weighting = True): if libsvm is None: raise RuntimeError('LibSVM Library not found. Cannot use this classifier.') self.param = libsvm.svm_parameter(kernel_type = libsvm.RBF, probability = probability) self.output_probability = probability self.auto_weighting = auto_weighting def set_option(self,optname,value): setattr(self.param, optname, value) def train(self, features, labels): labels,names = normaliselabels(labels) if self.auto_weighting: nlabels = labels.max() + 1 self.param.nr_weight = int(nlabels) self.param.weight_label = range(nlabels) self.param.weight = [(labels != i).mean() for i in xrange(nlabels)] problem = libsvm.svm_problem(labels.astype(float), features) model = libsvm.svm_model(problem, self.param) return libsvmModel(model, names, self.output_probability) milk-0.3.10/milk/active/0000755000175000017500000000000011562337213015537 5ustar luispedroluispedromilk-0.3.10/milk/active/__init__.py0000644000175000017500000000000011231643005017626 0ustar luispedroluispedromilk-0.3.10/milk/active/uncertainty.py0000644000175000017500000000420711452411266020460 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2010, Luis Pedro Coelho # vim: set ts=4 sts=4 sw=4 expandtab smartindent: # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. ''' Uncertainty ============ Implements uncertainty-based active learning strategies. These are strategies that are based on querying those elements in the pool which we are most uncertain about Functions ---------- * entropy * one_minus_max ''' from __future__ import division import numpy def entropy(model, pool): ''' entropies = entropy(model, pool) Returns the entropy of each classification output for members in the pool. ''' def _entropy(labels, ps): H = 0. for p in ps: if p > 1e-9: H += p * np.log(p) return H return [_entropy(model.apply(u)) for u in pool] def one_minus_max(model,pool): ''' oneminus = one_minus_max(model,pool) oneminus[i] = 1 - max_L { P(pool_i == L) } Returns one minus the probability for the best label guess. ''' def _minus1(labels, ps): return 1. - np.max(ps) return [_minus1(model.apply(u)) for u in pool] milk-0.3.10/milk/active/eimpact.py0000644000175000017500000000457511452411266017545 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2010, Luis Pedro Coelho # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. from __future__ import division import numpy from ..supervised import svm from ..supervised.classifier import ctransforms def expected_impacts(D,labels,U): ''' EIs = expected_impacts(D,labels,U) Compute Expected impact for each element of U Eis[i]: P(label[i] == 1) * IMPACT(label[i] == 1) + P(label[i] == 0) * IMPACT(label[i] == 0) ''' assert len(D) == len(labels), 'Nr of labeled examples should match lenght of labels vector' K = svm.rbf_kernel(20000) prob_classifier = ctransforms(svm.svm_raw(kernel=K,C=4),svm.svm_sigmoidal_correction()) label_classifier = ctransforms(svm.svm_raw(kernel=K,C=4),svm.svm_binary()) prob_classifier.train(D,labels) u_probs = prob_classifier(U) u_labels = (u_probs > .5) impacts = [] for u,p in zip(U,u_probs): print len(impacts) label_classifier.train(numpy.vstack((D,u)),numpy.hstack((labels,[0]))) u_labels_0 = label_classifier(U) label_classifier.train(numpy.vstack((D,u)),numpy.hstack((labels,[1]))) u_labels_1 = label_classifier(U) e_impact = (1.-p)*(u_labels != u_labels_0).sum() + p*(u_labels != u_labels_1).sum() impacts.append(e_impact) return impacts # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.3.10/setup.cfg0000644000175000017500000000007311562337213015151 0ustar luispedroluispedro[egg_info] tag_build = tag_date = 0 tag_svn_revision = 0