milk-0.1-alpha-0/0000755000175000017500000000000011302613311014313 5ustar luispedroluispedromilk-0.1-alpha-0/oldschoolclassify/0000755000175000017500000000000011302613311020037 5ustar luispedroluispedromilk-0.1-alpha-0/oldschoolclassify/nfoldcrossvalidation.py0000644000175000017500000000763211264757152024673 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008 Murphy Lab # Carnegie Mellon University # # Written by Luís Pedro Coelho # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published # by the Free Software Foundation; either version 2 of the License, # or (at your option) any later version. # # This program is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA # 02110-1301, USA. # # For additional information visit http://murphylab.web.cmu.edu or # send email to murphy@cmu.edu from __future__ import division import numpy import numpy as np from numpy import * from classify import defaultclassifier from classifier import normaliselabels __all__=['nfoldcrossvalidation'] def nfoldcrossvalidation(features, labels, nfolds=None, classifier=None, return_predictions=False, max_train=None): ''' Perform n-fold cross validation cmatrix,labelnames = nfoldcrossvalidation(features, labels, nfolds=10, classifier=None, return_predictions=False) cmatrix,labelnames,predictions = nfoldcrossvalidation(features, labels, nfolds=10, classifier=None, return_predictions=False) 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 labelnames[i] will correspond to the label name of class i @param features: a feature matrix or list of feature vectors @param labels: an array of labels, where label[i] is the label corresponding to features[i] @param nfolds: Nr of folds. @param return_predictions: whether to return predictions classifier should implement the train() and apply() methods ''' assert len(features) == len(labels), 'nfoldcrossvalidation: len(features) should match len(labels)' if classifier is None: classifier = defaultclassifier() labels,labelnames = normaliselabels(labels) predictions = labels*0-1 features = numpy.asanyarray(features) classcounts={} for L in labels: classcounts[L] = classcounts.get(L,0) + 1 min_class_count = min(classcounts.values()) if nfolds is None: nfolds = min(10,min_class_count) elif min_class_count < nfolds: from warnings import warn warn('pyslic.classify.nfoldcrossvalidation: Reducing the nr. of folds to %s (minimum class size).' % min_class_count) nfolds = min_class_count nclasses = len(classcounts) cmatrix = zeros((nclasses,nclasses)) for fold in xrange(nfolds): testingset=zeros(len(labels),bool) for L,C in classcounts.items(): idxs, = where(labels==L) N = len(idxs) perfold = N/nfolds start = floor(perfold*fold) end = floor(perfold*(fold+1)) idxs = idxs[start:end] testingset[idxs] = True trainingset= ~testingset if max_train is not None: for c in xrange(nclasses): if (labels[trainingset] == c).sum() > max_train: idxs, = np.where( (labels == L) & trainingset ) trainingset[idxs[max_train:]] = False classifier.train(features[trainingset],labels[trainingset]) prediction=classifier.apply(features[testingset]) predictions[testingset]=prediction for p, r in zip(prediction,labels[testingset]): cmatrix[r,p] += 1 if return_predictions: return cmatrix, labelnames, predictions return cmatrix, labelnames # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/oldschoolclassify/accuracy.py0000644000175000017500000000265511264757152022236 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008 Murphy Lab # Carnegie Mellon University # # Written by Luís Pedro Coelho # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published # by the Free Software Foundation; either version 2 of the License, # or (at your option) any later version. # # This program is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA # 02110-1301, USA. # # For additional information visit http://murphylab.web.cmu.edu or # send email to murphy@cmu.edu from __future__ import division from numpy import * __all__ = ['accuracy','waccuracy','rowsto1'] def rowsto1(cmatrix): return (cmatrix.T / cmatrix.sum(1)).T def accuracy(cmatrix): ''' acc = accuracy(cmatrix) Accuracy of cmatrix ''' return cmatrix.trace()/cmatrix.sum() def waccuracy(cmatrix): ''' wacc = waccuracy(cmatrix) Weighted accuracy of cmatrix ''' return (cmatrix.diagonal() / cmatrix.sum(1)).mean() # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/oldschoolclassify/classifywrap.py0000644000175000017500000000400211264757152023137 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008 Murphy Lab # Carnegie Mellon University # # Written by Luís Pedro Coelho # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published # by the Free Software Foundation; either version 2 of the License, # or (at your option) any later version. # # This program is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA # 02110-1301, USA. # # For additional information visit http://murphylab.web.cmu.edu or # send email to murphy@cmu.edu from __future__ import division __all__ = ['pretransformclassifier','concattransformers'] class concattransformers(object): def __init__(self,transformers): self.transformers = transformers def train(self,features,labels): for t in self.transformers: t.train(features,labels) features=t.apply(features) def apply(self,features): for t in self.transformers: features=t.apply(features) return features class pretransformclassifier(object): def __init__(self,transformer,classifier): if type(transformer) == list: transformer = concattransformers(transformer) self.transformer=transformer self.classifier=classifier def train(self,features,labels): self.transformer.train(features,labels) features=self.transformer.apply(features) self.classifier.train(features,labels) def apply(self,features): features=self.transformer.apply(features) return self.classifier.apply(features) # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/oldschoolclassify/gridsearch.py0000644000175000017500000000423511264757152022553 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008 Murphy Lab # Carnegie Mellon University # # Written by Luís Pedro Coelho # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published # by the Free Software Foundation; either version 2 of the License, # or (at your option) any later version. # # This program is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA # 02110-1301, USA. # # For additional information visit http://murphylab.web.cmu.edu or # send email to murphy@cmu.edu from __future__ import division from nfoldcrossvalidation import * def _allassignments(D): def allassignmentslist(L): if len(L) == 0: yield [] return key0,vals0=L[0] for v in vals0: for rest in allassignmentslist(L[1:]): yield [(key0,v)]+rest for A in allassignmentslist(list(D.items())): yield A def _set_assignment(obj,assignments): if not hasattr(obj,'set_option'): raise "Don't know how to set options" for k,v in assignments: obj.set_option(k,v) class gridsearch(object): def __init__(self,base,**kwargs): self.params = kwargs self.base = base self.best = None def is_multi_class(self): return self.base.is_multi_class() def train(self,features,labels): best_trace=-1 for assignement in _allassignments(self.params): _set_assignment(self.base,assignement) S,_=nfoldcrossvalidation(features,labels,classifier=self.base) if S.trace() > best_trace: self.best=assignement best_trace=S.trace() _set_assignment(self.base,self.best) def apply(self,features): return self.base.apply(features) milk-0.1-alpha-0/oldschoolclassify/libsvmclassifier.py0000644000175000017500000000577011264757152024006 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008 Murphy Lab # Carnegie Mellon University # # Written by Luís Pedro Coelho # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published # by the Free Software Foundation; either version 2 of the License, # or (at your option) any later version. # # This program is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA # 02110-1301, USA. # # For additional information visit http://murphylab.web.cmu.edu or # send email to murphy@cmu.edu from __future__ import division from numpy import array, asanyarray from classifier import classifier, classification_result from tempfile import NamedTemporaryFile _svm = None try: import svm _svm = svm except: pass try: from libsvm import svm _svm = svm except: pass __all__=['libsvmClassifier'] class libsvmClassifier(classifier): def __init__(self,probability = False, auto_weighting = True): classifier.__init__(self) if _svm is None: raise RuntimeError('SVM Library not found. Cannot use this classifier.') self.param = _svm.svm_parameter(kernel_type = svm.RBF, probability = probability) self.output_probability = probability self.auto_weighting = auto_weighting def set_option(self,optname,value): setattr(self.param,optname,value) def _dotrain(self,features,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 = _svm.svm_problem(labels.astype(float), features) self.model = _svm.svm_model(problem,self.param) def apply(self,feats): if len(feats.shape) == 2: return [self.apply(f) for f in feats] if self.output_probability: return self.model.predict_probability(feats) return self.model.predict(feats) 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._trained,self._labelnames def __setstate__(self,state): if _svm is None: raise RuntimeError('LibSVM Library not found. Cannot use this classifier.') S,self.output_probability,self._trained,self._labelnames = state N=NamedTemporaryFile() N.write(S) N.flush() self.model = _svm.svm_model(N.name) # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/oldschoolclassify/__init__.py0000644000175000017500000000231311264757152022172 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008 Murphy Lab # Carnegie Mellon University # # Written by Luís Pedro Coelho # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published # by the Free Software Foundation; either version 2 of the License, # or (at your option) any later version. # # This program is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA # 02110-1301, USA. # # For additional information visit http://murphylab.web.cmu.edu or # send email to murphy@cmu.edu from accuracy import accuracy, waccuracy, rowsto1 from classify import * from libsvmclassifier import libsvmClassifier from classifywrap import concattransformers, pretransformclassifier from nfoldcrossvalidation import nfoldcrossvalidation import normalise from normalise import * milk-0.1-alpha-0/oldschoolclassify/multi.py0000644000175000017500000000656211264757152021577 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008 Murphy Lab # Carnegie Mellon University # # Written by Luís Pedro Coelho # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published # by the Free Software Foundation; either version 2 of the License, # or (at your option) any later version. # # This program is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA # 02110-1301, USA. # # For additional information visit http://murphylab.web.cmu.edu or # send email to murphy@cmu.edu from __future__ import division from numpy import * from classifier import classifier import random class one_against_rest(classifier): __slots__ = ['_classifiers','_base','_base2'] def __init__(self,base): classifier.__init__(self) self._classifiers = None self._base = base self._base2 = None def is_multi_class(self): return True def _dotrain(self,features,labels): if labels.max() + 1 == 2: self._base2 = self._base() self._base2.train(features,labels) return #else self._classifiers=[] for i in xrange(labels.max()+1): s = self._base() s.train(features,labels==i) self._classifiers.append(s) def _doapply(self,feats): if self._classifiers is None: return self._base2.apply(feats) vals=array([c.apply(feats) for c in self._classifiers]) idxs, = where(vals == 1) if len(idxs) == 1: return idxs[0] elif len(idxs) == 0: return random.randint(0,self._nclasses) else: # choose at random return random.choice(idxs) class one_against_one(classifier): __slots__ = ['_classifiers','_base','_labelnames'] def __init__(self,factory): classifier.__init__(self) self._classifiers = None self._base = factory def _is_multiclass(self): return True def _dotrain(self,features,labels): nc = labels.max()+1 self._nclasses = nc if nc == 2: self._base = self._base() self._base.train(features,labels) return #else self._classifiers=[[None for i in xrange(nc)] for j in xrange(nc)] for i in xrange(nc): for j in xrange(i+1,nc): s=self._base() idxs=(labels == i) | (labels == j) s.train(features[idxs],labels[idxs]==i) self._classifiers[i][j]=s def _doapply(self,feats): if self._classifiers is None: return self._base.apply(feats) nc=self._nclasses votes=zeros(nc) for i in xrange(nc): for j in xrange(i+1,nc): c=self._classifiers[i][j].apply(feats) if c: votes[i] += 1 else: votes[j] += 1 return votes.argmax(0) # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/oldschoolclassify/pymlclassifier.py0000644000175000017500000000340711264757152023466 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008 Murphy Lab # Carnegie Mellon University # # Written by Luís Pedro Coelho # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published # by the Free Software Foundation; either version 2 of the License, # or (at your option) any later version. # # This program is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA # 02110-1301, USA. # # For additional information visit http://murphylab.web.cmu.edu or # send email to murphy@cmu.edu from numpy import array from PyML import * __ALL__=['PyMLSVM'] class PyMLSVM(object): def __init__(self): pass def train(self,feats,labels): allabels=dict([(L,None) for L in labels]) if len(allabels) > 2: self.s = multi.OneAgainstRest(svm.SVM()) else: self.s = svm.SVM() nlabels=[] for L in labels: nlabels.append(str(L)) data=datafunc.VectorDataSet(feats,L=nlabels) self.s.train(data) self.trained=True def apply(self,feats): assert self.trained if feats.ndim == 1: data=datafunc.VectorDataSet([feats]) return array(self.s.classify(data,0))[0] else: return array([self.apply(feats[i]) for i in xrange(feats.shape[0])]) # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/oldschoolclassify/classify.py0000644000175000017500000000605111264757152022253 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008 Murphy Lab # Carnegie Mellon University # # Written by Luís Pedro Coelho # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published # by the Free Software Foundation; either version 2 of the License, # or (at your option) any later version. # # This program is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA # 02110-1301, USA. # # For additional information visit http://murphylab.web.cmu.edu or # send email to murphy@cmu.edu from __future__ import division __ALL__ = ['learnclassifier','fastclassifier','applyclassifier','defaultclassifier'] def defaultclassifier(): ''' C = defaultclassifier() Returns the default classifier ''' from normalise import zscore_normalise,chkfinite,interval_normalise from gridsearch import gridsearch from classifywrap import pretransformclassifier from libsvmclassifier import libsvmClassifier from featureselection import sda_filter, remove_repeated_features, remove_linear_dependent_features from numpy import arange return pretransformclassifier( [chkfinite(),interval_normalise(),remove_linear_dependent_features(),sda_filter()], # pre-process gridsearch(libsvmClassifier(),C=2.**arange(-4,10), gamma=2.**arange(-7,2)) # classify ) def fastclassifier(): ''' C = fastclassifier() Returns the classifier, which is faster than the defaultclassifier() Currently, this means that this classifier does not perform a grid search for the best libsvm parameters ''' from normalise import zscore_normalise,chkfinite,interval_normalise from gridsearch import gridsearch from classifywrap import pretransformclassifier from libsvmclassifier import libsvmClassifier from pymlclassifier import PyMLSVM from featureselection import sda_filter, remove_repeated_features, remove_linear_dependent_features from numpy import arange return pretransformclassifier([chkfinite(),interval_normalise(),remove_linear_dependent_features(),sda_filter()],libsvmClassifier()) def learnclassifier(featmatrix,y): """ classifier=learnclassifier(featmatrix,y) Learn a classifier for the problem y[i]=f(featmatrix[i]) """ classifier=defaultclassifier() classifier.train(featmatrix,y) return classifier def applyclassifier(feats,model): """ y = applyclassifier(feats,models) output a label for feats according to the model feats can be either a feature array or a sequence of feature arrays """ return model.apply(feats) # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/oldschoolclassify/featureselection.py0000644000175000017500000002026511264757152024002 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008 Murphy Lab # Carnegie Mellon University # # Written by Luís Pedro Coelho # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published # by the Free Software Foundation; either version 2 of the License, # or (at your option) any later version. # # This program is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA # 02110-1301, USA. # # For additional information visit http://murphylab.web.cmu.edu or # send email to murphy@cmu.edu from __future__ import division from numpy import * from numpy.linalg import det from classifier import normaliselabels, classifier import scipy.stats import warnings _TOLERANCE = 0 _SIGNIFICANCE_IN = .15 _SIGNIFICANCE_OUT = .15 def _sweep(A,k,flag): N,_=A.shape Akk=A[k,k] B=zeros_like(A) try: from scipy import weave from scipy.weave import converters k=int(k) code = ''' #line 32 "featureselection.py" for (int i = 0; i != N; ++i) { for (int j = 0; j != N; ++j) { if (i == k) { if (j == k) { B(i,j) = - 1./A(k,k); } else { B(i,j)=flag*A(i,j)/A(k,k); } } else if (j == k) { B(i,j)=flag*A(i,j)/A(k,k); } else { B(i,j)=A(i,j) - A(i,k)*A(k,j)/A(k,k); } } } ''' weave.inline( code, ['A','B','k','Akk','flag','N'], type_converters=converters.blitz) except: for i in xrange(N): for j in xrange(N): if i == k: if j == k: B[i,j] = - 1./Akk else: B[i,j]=flag*A[i,j]/Akk elif j == k: B[i,j]=flag*A[i,j]/Akk else: B[i,j]=A[i,j] - A[i,k]*A[k,j]/Akk return B def sda(features,labels): ''' 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. ''' N, m = features.shape labels,labelsu = normaliselabels(labels) q=len(labelsu) mus=array([features[labels==i,:].mean(0) for i in xrange(q)]) mu=features.mean(0) W=zeros((m,m)) T=zeros((m,m)) try: from scipy import weave from scipy.weave import converters code=''' #line 68 "featureselection.py" for (int i = 0; i != m; ++i) { for (int j = 0; j != m; ++j) { for (int n = 0; n != N; ++n) { int g=labels(n); W(i,j) += (features(n,i)-mus(g,i))*(features(n,j)-mus(g,j)); T(i,j) += (features(n,i)-mu(i))*(features(n,j)-mu(j)); } } } ''' weave.inline( code, ['N','m','W','T','features','mu','mus','labels'], type_converters=converters.blitz) except: warnings.warn('scipy.weave failed. Resorting to (slow) Python code') for i in xrange(m): for j in xrange(m): for n in xrange(N): g=labels[n] W[i,j] += (features[n,i]-mus[g,i])*(features[n,j]-mus[g,j]) T[i,j] += (features[n,i]-mu[i])*(features[n,j]-mu[j]) ignoreidx = ( W.diagonal() == 0 ) if ignoreidx.any(): idxs, = where(~ignoreidx) F=sda(features[:,~ignoreidx],labels) return idxs[F] output=[] D=W.diagonal() df1 = q-1 while True: V = W.diagonal()/T.diagonal() W_d = W.diagonal() V_neg = (W_d < 0) p=V_neg.sum() df2 = N-p-q+1 if V_neg.any(): V_m = V[V_neg].min() k,=where(V == V_m) k=k[0] Fremove = (N-p-q+1)/(q-1)*(V_m-1) PrF = 1 - scipy.stats.f.cdf(Fremove,df1,df2) if PrF > _SIGNIFICANCE_OUT: #print 'removing ',k, 'V(k)', 1./V_m, 'Fremove', Fremove 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,=where(V==V_m) k=k[0] Fenter = (N-p-q)/(q-1) * (1-V_m)/V_m PrF = 1 - scipy.stats.f.cdf(Fenter,df1,df2) if PrF < _SIGNIFICANCE_IN: #print 'adding ',k, 'V(k)', 1./V_m, 'Fenter', Fenter W=_sweep(W,k,-1) T=_sweep(T,k,-1) if PrF < .0001: output.append((Fenter,k)) continue break output.sort(reverse=True) return array([x[1] for x in output]) def repeatedfeatures(featmatrix): ''' idxs = repeatedfeatures(feature_matrix) Returns a set of indices corresponding to repeated features in featmatrix. @see remove_repeated_features ''' featsig={} results=[] feat_sum=(featmatrix**2).sum(0) ordering=feat_sum.argsort() i = 0 #while i < len(ordering) and feat_sum[ordering[i]] == 0: # results.append(ordering[i]) # i += 1 while i < len(ordering) - 1: if feat_sum[ordering[i]] == feat_sum[ordering[i+1]]: if (featmatrix[:,ordering[i]] == featmatrix[:,ordering[i+1]]).all(): results.append(ordering[i+1]) i += 1 return array(results) def _rank(A,tol=1e-8): s = linalg.svd(A,compute_uv=0) return (s > tol).sum() def linearindependentfeatures(featmatrix): ''' Returns the indices of a set of linearly independent features (columns). indices = linearindependentfeatures(features) ''' independent=[] R=_rank(featmatrix) i=0 offset=0 while i < featmatrix.shape[1]: R_=_rank(delete(featmatrix,i,1)) if R_ == R: featmatrix=delete(featmatrix,i,1) offset+=1 else: independent.append(i+offset) i += 1 return array(independent) class sda_filter(object): __slots__ = ['idxs'] def is_multi_class(self): return True def __init__(self): pass def train(self,features,labels): self.idxs=sda(features,labels) if len(self.idxs) == 0: self.idxs = [0] def apply(self,features): return features[:,self.idxs] def __getstate__(self): return (self.idxs,) def __setstate__(self,state): self.idxs, = state class remove_repeated_features(object): __slots__ = ['repeats'] def is_multi_class(self): return True def __init__(self): pass def train(self,features,labels): self.repeats=repeatedfeatures(features) def apply(self,features): return delete(features,self.repeats,1) def __getstate__(self): return (self.repeats,) def __setstate__(self,state): self.repeats, = state class remove_linear_dependent_features(object): __slots__ = ['independent'] def is_multi_class(self): return True def __init__(self): pass def train(self,features,labels): self.independent=linearindependentfeatures(features) def apply(self,features): return features[:,self.independent] def __getstate__(self): return (self.independent,) def __setstate__(self,state): self.independent, = state milk-0.1-alpha-0/oldschoolclassify/classifier.py0000644000175000017500000000450011264757152022557 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008 Murphy Lab # Carnegie Mellon University # # Written by Luís Pedro Coelho # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published # by the Free Software Foundation; either version 2 of the License, # or (at your option) any later version. # # This program is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA # 02110-1301, USA. # # For additional information visit http://murphylab.web.cmu.edu or # send email to murphy@cmu.edu from __future__ import division from numpy import * __all__ = ['classification_result','normaliselabels','classifier'] class classification_result(object): __slots__ = ['label','p','funcvalue'] def __init__(self): pass def normaliselabels(labels): labelnames={} normalised=[] names=[] N=0 for L in labels: nr=labelnames.get(L,N) if nr == N: labelnames[L]=N names.append(L) N += 1 normalised.append(nr) return array(normalised),names class classifier(object): def __init__(self): self._trained = False def train(self,features,labels): nlabels,self._labelnames=normaliselabels(labels) self._nclasses = nlabels.max()+1 features=asanyarray(features) self._dotrain(features,nlabels) self._trained = True def apply(self,features): assert self._trained features=asanyarray(features) if features.ndim == 1: return self._labelnames[int(self._doapply(features))] else: return array([self._labelnames[int(self._doapply(features[i]))] for i in xrange(features.shape[0])]) class results_to_labels(classifier): def __init__(self): pass def _dotrain(self,features,labels): pass def _doapply(self,res): return res.label # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/oldschoolclassify/normalise.py0000644000175000017500000001012611264757152022425 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008 Murphy Lab # Carnegie Mellon University # # Written by Luís Pedro Coelho # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published # by the Free Software Foundation; either version 2 of the License, # or (at your option) any later version. # # This program is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA # 02110-1301, USA. # # For additional information visit http://murphylab.web.cmu.edu or # send email to murphy@cmu.edu from __future__ import division import numpy from numpy import * from bisect import bisect_right try: import ncreduce std=ncreduce.std except: std=numpy.std __all__ = ['zscore','zscore_normalise','interval_normalise','chkfinite','icdf_normalise'] def zscore(features): """ features = zscore(features) Returns a copy of features which has been normalised to zscores """ mu=features.mean(0) sigma=std(features,0) sigma[sigma == 0] = 1 return (features - mu) / sigma class zscore_normalise(object): ''' Normalise to z-scores A preprocessor that normalises features to z scores. ''' __slots__=['mu','sigma'] def __init__(self,features=None): if features: self.train(features) def train(self,features,labels): self.mu=features.mean(0) self.sigma=std(features,0) self.sigma[self.sigma == 0.]=1 # This makes the division a null op. def apply(self,features): return (features - self.mu)/self.sigma class interval_normalise(object): ''' Linearly scale to the interval [-1,1] (per libsvm recommendation) ''' __slots__=['mu','range'] def __init__(self,features=None): if features: self.train(features) def train(self,features,labels): self.mu=features.mean(0) self.range=features.max(0)-features.min(0) self.range[self.range == 0.]=1 # This makes the division a null op. def apply(self,features): return (features - self.mu)/self.range def __getstate__(self): return self.mu,self.range def __setstate__(self,state): self.mu,self.range = state class chkfinite(object): ''' Fill NaN & Inf values Replaces NaN & Inf values with zeros. ''' __slots__ = [] def __init__(self): pass def train(self,features,labels): pass def apply(self,features): nans=isnan(features)+isinf(features) if nans.any(): features=features.copy() features[nans]=0 return features def _do_one_icdf_normalise(values): values=values.copy() N=len(values) argsorted=values.argsort() corrected=numpy.empty(N) try: from scipy import weave from scipy.weave import converters code=''' for (int i = 0; i != N; ++i) { corrected(argsorted(i)) = i/double(N); } ''' weave.inline(code, ['corrected','argsorted','N'], type_converters=converters.blitz) except: corrected[argsorted]=arange(N)/N return corrected def icdf_normalise(fmatrix): ''' icdf_normalise(fmatrix) Normalise each feature so that each value is replaced by the inverse of the feature's cummulative distribution function. ''' if len(fmatrix.shape) == 1: return _do_one_icdf_normalise(fmatrix) elif len(fmatrix.shape) == 2: _,q=fmatrix.shape for i in xrange(q): fmatrix[:,i]=_do_one_icdf_normalise(fmatrix[:,i]) return fmatrix else: raise ValueError, "normalise_to_01 only works with 1 or 2 dimensional arrays!" # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/README0000644000175000017500000000075111302613267015210 0ustar luispedroluispedro============================== MILK: MACHINE LEARNING TOOLKIT ============================== Machine Learning in Python -------------------------- Milk is a machine learning toolkit in Python. Its focus is on supervised classification. Features -------- - SVMs. Using the libsvm solver with a pythonesque wrapper around it - Stepwise Discriminant Analysis for feature selection License: MIT Author: Luis Pedro Coelho (with code from LibSVM) Website: http://luispedro.org/software/milk milk-0.1-alpha-0/setup.py0000644000175000017500000000363411302613310016032 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2009, Luís 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 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 svm_ext = Extension('milk.supervised._svm', sources = ['milk/supervised/_svm.cpp']) ext_modules = [svm_ext] packages = setuptools.find_packages() if 'tests' in packages: packages.remove('tests') setup(name = 'milk', version = '0.1-alpha-0', description = 'Machine Learning Toolkit', author = u'Luis Pedro Coelho', author_email = 'lpc@mcu.edu', url = 'http://luispedro.org/software/milk', packages = packages, ext_modules = ext_modules, ) # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/PKG-INFO0000644000175000017500000000036211302613311015411 0ustar luispedroluispedroMetadata-Version: 1.0 Name: milk Version: 0.1-alpha-0 Summary: Machine Learning Toolkit Home-page: http://luispedro.org/software/milk Author: Luis Pedro Coelho Author-email: lpc@mcu.edu License: UNKNOWN Description: UNKNOWN Platform: UNKNOWN milk-0.1-alpha-0/milk.egg-info/0000755000175000017500000000000011302613311016741 5ustar luispedroluispedromilk-0.1-alpha-0/milk.egg-info/dependency_links.txt0000644000175000017500000000000111302613311023007 0ustar luispedroluispedro milk-0.1-alpha-0/milk.egg-info/SOURCES.txt0000644000175000017500000000274611302613311020636 0ustar luispedroluispedroREADME setup.py milk/__init__.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/roymccallum.py milk/active/uncertainty.py milk/measures/__init__.py milk/measures/measures.py milk/measures/nfoldcrossvalidation.py milk/supervised/__init__.py milk/supervised/_svm.cpp milk/supervised/classifier.py milk/supervised/defaultclassifier.py milk/supervised/featureselection.py milk/supervised/gridsearch.py milk/supervised/knn.py milk/supervised/libsvmclassifier.py milk/supervised/multi.py milk/supervised/normalise.py milk/supervised/parzen.py milk/supervised/setup.py milk/supervised/svm.py milk/supervised/tree.py milk/unsupervised/__init__.py milk/unsupervised/gaussianmixture.py milk/unsupervised/kmeans.py milk/unsupervised/normalise.py milk/unsupervised/pca.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 oldschoolclassify/__init__.py oldschoolclassify/accuracy.py oldschoolclassify/classifier.py oldschoolclassify/classify.py oldschoolclassify/classifywrap.py oldschoolclassify/featureselection.py oldschoolclassify/gridsearch.py oldschoolclassify/libsvmclassifier.py oldschoolclassify/multi.py oldschoolclassify/nfoldcrossvalidation.py oldschoolclassify/normalise.py oldschoolclassify/pymlclassifier.py tests/data/__init__.py tests/data/german/__init__.py tests/data/german/german.pymilk-0.1-alpha-0/milk.egg-info/top_level.txt0000644000175000017500000000003511302613311021471 0ustar luispedroluispedrooldschoolclassify tests milk milk-0.1-alpha-0/milk.egg-info/PKG-INFO0000644000175000017500000000036211302613311020037 0ustar luispedroluispedroMetadata-Version: 1.0 Name: milk Version: 0.1-alpha-0 Summary: Machine Learning Toolkit Home-page: http://luispedro.org/software/milk Author: Luis Pedro Coelho Author-email: lpc@mcu.edu License: UNKNOWN Description: UNKNOWN Platform: UNKNOWN milk-0.1-alpha-0/tests/0000755000175000017500000000000011302613311015455 5ustar luispedroluispedromilk-0.1-alpha-0/tests/data/0000755000175000017500000000000011302613311016366 5ustar luispedroluispedromilk-0.1-alpha-0/tests/data/__init__.py0000644000175000017500000000000011231643007020473 0ustar luispedroluispedromilk-0.1-alpha-0/tests/data/german/0000755000175000017500000000000011302613311017637 5ustar luispedroluispedromilk-0.1-alpha-0/tests/data/german/__init__.py0000644000175000017500000000000011231643007021744 0ustar luispedroluispedromilk-0.1-alpha-0/tests/data/german/german.py0000644000175000017500000001016611231643007021474 0ustar luispedroluispedro#! /usr/bin/env python # -*- coding: utf-8 -*- # Last Change: Tue Jul 17 05:00 PM 2007 J # The code and descriptive text is copyrighted and offered under the terms of # the BSD License from the authors; see below. However, the actual dataset may # have a different origin and intellectual property status. See the SOURCE and # COPYRIGHT variables for this information. # Copyright (c) 2007 David Cournapeau # # All rights reserved. # # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions are # met: # # * Redistributions of source code must retain the above copyright # notice, this list of conditions and the following disclaimer. # * Redistributions in binary form must reproduce the above copyright # notice, this list of conditions and the following disclaimer in # the documentation and/or other materials provided with the # distribution. # * Neither the author nor the names of any contributors may be used # to endorse or promote products derived from this software without # specific prior written permission. # # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS # "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED # TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR # PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR # CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, # EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, # PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; # OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, # WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR # OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF # ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. """Iris dataset.""" __docformat__ = 'restructuredtext' COPYRIGHT = """This is public domain. """ TITLE = "German Dataset" SOURCE = """ http://www.liacc.up.pt/ML/old/statlog/datasets.html Professor Dr. Hans Hofmann Institut für Statistik und Ökonometrie Universität Hamburg FB Wirtschaftswissenschaften Von-Melle-Park 5 2000 Hamburg 13 Two datasets are provided. the original dataset, in the form provided by Prof. Hofmann, contains categorical/symbolic attributes and is in the file "german.dat". For algorithms that need numerical attributes, Strathclyde University produced the file "german.numer". This file has been edited and several indicator variables added to make it suitable for algorithms which cannot cope with categorical variables. Several attributes that are ordered categorical (such as attribute 17) have been coded as integer. This was the form used by StatLog. Here (scikits), only the numeric datasets are provided. """ DESCRSHORT = """""" DESCRLONG = """""" NOTE = """ Number of Instances: 1000. 700 for class 0 (good credit) and 300 for class 1 (bad credit). Number of Attributes: 24. label: 0 for good credit, +1 for bad credit """ def load(): """load the german data and returns them. :returns: data: dict Contains the following values: 'data' : the actual data 'label' : label[i] is the label index of data[i] 'class' : class[label[i]] is the label name of data[i] """ import numpy import pickle import gzip from os.path import dirname, join feat,label = pickle.load(gzip.GzipFile(join(dirname(__file__), 'german.pp.gz'))) data = {} # Create the data data['data'] = numpy.array([[float(j) for j in f] for f in feat.values()]).T # Create the label data['label'] = numpy.array([numpy.int(i) for i in label], numpy.int) data['label'][data['label'] == -1] = 0 # Create the class classes = ['good credit', 'bad credit'] data['class'] = numpy.empty(2, 'S%d' % numpy.max([len(i) for i in classes])) data['class'][0] = 'good credit' data['class'][1] = 'bad credit' #data['data'] = numpy.core.records.fromrecords(data['data']) return data milk-0.1-alpha-0/milk/0000755000175000017500000000000011302613311015247 5ustar luispedroluispedromilk-0.1-alpha-0/milk/measures/0000755000175000017500000000000011302613311017073 5ustar luispedroluispedromilk-0.1-alpha-0/milk/measures/nfoldcrossvalidation.py0000644000175000017500000001166011277673337023731 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008, Luís 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 from collections import defaultdict from ..supervised.classifier import normaliselabels from ..supervised.defaultclassifier import defaultclassifier import numpy import numpy as np __all__ = ['foldgenerator', 'getfold', 'nfoldcrossvalidation'] def foldgenerator(labels,nfolds=None): ''' for test,train in foldgenerator(labels,nfolds=None) ... This generator breaks up the data into nfolds (default 10). Params ------ * labels: the labels * nfolds: nr of folds (default 10) ''' classcounts = defaultdict(int) for L in labels: classcounts[L] += 1 min_class_count = min(classcounts.values()) if nfolds is None: nfolds = min(10,min_class_count) elif min_class_count < nfolds: from warnings import warn warn('milk.measures.nfoldcrossvalidation: Reducing the nr. of folds to %s (minimum class size).' % min_class_count) nfolds = min_class_count testingset = np.empty(len(labels),bool) for fold in xrange(nfolds): testingset.fill(False) for L,C in classcounts.items(): idxs, = np.where(labels==L) N = len(idxs) perfold = N/nfolds start = np.floor(perfold*fold) end = np.floor(perfold*(fold+1)) idxs = idxs[start:end] testingset[idxs]=True trainingset = ~testingset yield trainingset,testingset def getfold(labels, fold, nfolds=None): ''' trainingset,testingset = getfold(labels, fold, nfolds=None) ''' _nfolds = (10 if nfolds is None else nfolds) assert fold < _nfolds, 'milk.getfold: Attempted to get fold %s out of %s' % (fold, nfolds) for i,(t,s) in enumerate(foldgenerator(labels,nfolds)): if i == fold: return t,s assert False, 'milk.getfold: Attempted to get fold %s but the number of actual folds was too small' % fold def nfoldcrossvalidation(features,labels,nfolds=None,classifier=None, return_predictions=False): ''' Perform n-fold cross validation cmatrix,labelnames = nfoldcrossvalidation(features, labels, nfolds=10, classifier=None, return_predictions=False) cmatrix,labelnames,predictions = nfoldcrossvalidation(features, labels, nfolds=10, classifier=None, return_predictions=False) 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 labelnames[i] will correspond to the label name of class i Arguments --------- * features: a feature matrix or list of feature vectors * labels: an array of labels, where label[i] is the label corresponding to features[i] * nfolds: Nr of folds. * return_predictions: whether to return predictions classifier should implement the train() method to return a model (something with an apply() method) ''' assert len(features) == len(labels), 'milk.measures.nfoldcrossvalidation: len(features) should match len(labels)' if classifier is None: classifier = defaultclassifier() labels,labelnames = normaliselabels(labels) predictions = np.zeros_like(labels)-1 features = numpy.asanyarray(features) nclasses = labels.max() + 1 cmatrix = np.zeros((nclasses,nclasses)) for trainingset,testingset in foldgenerator(labels, nfolds): model = classifier.train(features[trainingset], labels[trainingset]) prediction = np.array([model.apply(f) for f in features[testingset]]) predictions[testingset] = prediction.astype(predictions.dtype) for p, r in zip(prediction,labels[testingset]): cmatrix[r,p] += 1 if return_predictions: return cmatrix, labelnames, predictions return cmatrix, labelnames # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/measures/measures.py0000644000175000017500000000332711231643005021302 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008, Luís 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 __all__ = ['accuracy','waccuracy'] def accuracy(cmatrix): ''' acc = accuracy(cmatrix) Accuracy of confusion matrix ''' cmatrix = numpy.asanyarray(cmatrix) return cmatrix.trace()/cmatrix.sum() def waccuracy(cmatrix): ''' wacc = waccuracy(cmatrix) Weighted accuracy of cmatrix ''' cmatrix = numpy.asanyarray(cmatrix) return (cmatrix.diagonal() / cmatrix.sum(1)).mean() ## TODO: Implement http://en.wikipedia.org/wiki/Matthews_Correlation_Coefficient # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/measures/__init__.py0000644000175000017500000000002711256470614021222 0ustar luispedroluispedrofrom measures import * milk-0.1-alpha-0/milk/supervised/0000755000175000017500000000000011302613311017440 5ustar luispedroluispedromilk-0.1-alpha-0/milk/supervised/parzen.py0000644000175000017500000000415711231643006021325 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2009, Luís 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 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 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 # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/supervised/svm.py0000644000175000017500000003420511300556020020624 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008, Luís 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 from .classifier import normaliselabels from collections import deque import numpy import numpy as np import random from . import _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 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: def _k(q): return np.array([self.kernel(s, q) for s in self.svs]) self.kernelfunction = _k 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) MAJOR 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 MINOR 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): 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 Y,_ = normaliselabels(labels) assert numpy.all( (Y == 0) | (Y == 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. Major Arguments --------------- * F: Values of the function F * Y: Labels (in boolean format, ie, in (0,1)) Minor Arguments --------------- * max_iters: Maximum nr. of iterations * min_step: Minimum step * sigma: sigma * eps: A small number REFERENCE --------- 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 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): def train(self, features, labels): assert len(labels) >= 2, 'Cannot train from a single example' labelset = set(labels) assert len(labelset) == 2, 'milk.supervised.svm.svm_binary.train: Can only handle two class problems' c0,c1 = labelset if c0 != labels[0]: c0,c1 = c1, c0 return svm_binary_model( (c0,c1) ) 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): 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: the feature matrix. 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 N1 = int(N1) N2 = int(N2) try: from scipy import weave from scipy.weave import converters res = np.empty(3,np.double) code = ''' #line 400 "svm.py" double C1v = 0.; double C2v = 0.; double C12v = 0.; for (int i = 0; i != N1; ++i) for (int j = 0; j != N1; ++j) C1v += std::exp(C1(i,j)/sigma); for (int i = 0; i != N2; ++i) for (int j = 0; j != N2; ++j) C2v += std::exp(C2(i,j)/sigma); for (int i = 0; i != N1; ++i) for (int j = 0; j != N2; ++j) C12v += std::exp(C12(i,j)/sigma); res(0) = C1v/N1; res(1) = C2v/N2; res(2) = C12v/N1/N2; ''' weave.inline(code, ['C1','C2','C12','N1','N2','res','sigma'],type_converters = converters.blitz) C1v = res[0] C2v = res[1] C12v = res[2] except Exception, e: 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): 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) # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/supervised/knn.py0000644000175000017500000000442511277731540020625 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008, Luís 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 from collections import defaultdict import numpy as np class kNN(object): ''' k-Nearest Neighbour Classifier Naive implementation of a k-nearest neighbour classifier. C = kNN(k) Attributes: ----------- * k: the k to use ''' def __init__(self, k=1): self.k = k def train(self, features, labels, copy_features=False): features = np.asanyarray(features) labels = np.asanyarray(labels) if copy_features: features = features.copy() labels = labels.copy() return kNN_model(self.k, features, labels) class kNN_model(object): def __init__(self, k, features, labels): self.k = k self.features = features self.labels = labels def apply(self,features): diff2 = ( (self.features - features)**2 ).sum(1) 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 # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/supervised/setup.py0000644000175000017500000000272611231643006021166 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008, Luís 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 numpy.distutils.core import setup, Extension svm_ext= Extension('_svm', sources = ['_svm.cpp']) if __name__ == '__main__': setup (name = 'svm', version = '0.1', description = "Implement Platt's SMO Algorithm", ext_modules = [svm_ext] ) # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/supervised/defaultclassifier.py0000644000175000017500000000570711277666702023543 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2009, Luís 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 def defaultclassifier(mode='medium'): ''' C = defaultclassifier(mode='medium') Parameters ----------- * mode: One of ('fast','medium','slow'). This defines the speed accuracy trade-off. It essentially defines how large the SVM parameter range is. ''' 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.defaultclassifier: 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( lambda: ctransforms( \ svm.svm_raw(), \ svm.svm_binary() \ )), params={ (0,'C'): 2.**c_range, (0,'kernel'): [svm.rbf_kernel(2.**i) for i in sigma_range], } )) # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/supervised/gridsearch.py0000644000175000017500000000552611277673730022163 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2009, Luís 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 from ..measures.nfoldcrossvalidation import nfoldcrossvalidation def _allassignments(D): def allassignmentslist(L): if len(L) == 0: yield [] return key0,vals0 = L[0] for v in vals0: for rest in allassignmentslist(L[1:]): yield [(key0,v)]+rest for A in allassignmentslist(list(D.items())): yield A def _set_assignment(obj,assignments): for k,v in assignments: obj.set_option(k,v) class gridsearch(object): ''' G = gridsearch(base_classifier, measure=accuracy, param1=[1,2,3], param2=['a','b','c'],...) Perform a grid search for the best parameter values. Parameters ----------- * measure: a function which takes a confusion matrix and outputs a good-ness measure ''' def __init__(self, base, measure=None, params={}): self.params = params self.base = base self.best = None if measure is None: measure = np.trace self.measure = measure def is_multi_class(self): return self.base.is_multi_class() def train(self,features,labels): best_val = -1 for assignement in _allassignments(self.params): _set_assignment(self.base, assignement) #print 'gridsearch:', assignement S,_ = nfoldcrossvalidation(features, labels, classifier=self.base) cur = self.measure(S) if cur > best_val: self.best = assignement best_val = cur _set_assignment(self.base, self.best) return self.base.train(features, labels) # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/supervised/libsvmclassifier.py0000644000175000017500000000547311277734266023414 0ustar luispedroluispedro#s -*- coding: utf-8 -*- # Copyright (C) 2008 Murphy Lab # Carnegie Mellon University # # Written by Luís 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 from numpy import array, asanyarray _svm = None try: from libsvm import svm _svm = svm except: pass __all__=['libsvmClassifier'] class libsvmClassifier(object): def __init__(self,probability = False): classifier.__init__(self) if _svm is None: raise RuntimeError('SVM Library not found. Cannot use this classifier.') self.param = _svm.svm_parameter(kernel_type = svm.RBF, probability = probability) self.output_probability = False def set_option(self,optname,value): setattr(self.param,optname,value) def train(self,features,labels): problem = _svm.svm_problem(labels,features) return libsvmClassifier_model(_svm.svm_model(problem,self.param)) class libsvmClassifier_model(object): def __init__(self, model): self.model = model def apply(self,feats): if self.output_probability: return self.model.predict_probability(feats) return self.model.predict(feats) 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._trained,self._labelnames def __setstate__(self,state): if _svm is None: raise RuntimeError('LibSVM Library not found. Cannot use this classifier.') S,self.output_probability,self._trained,self._labelnames = state N=NamedTemporaryFile() N.write(S) N.flush() self.model = _svm.svm_model(N.name) # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/supervised/__init__.py0000644000175000017500000000234611256467602021600 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2009, Luís 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 .defaultclassifier import defaultclassifier # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/supervised/multi.py0000644000175000017500000001251711277734311021171 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008, Luís 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 from .classifier import normaliselabels import numpy as np import random 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 be a callable object that provides the base classifier to use. Example ------- multi = one_against_rest(milk.supervised.simple_svm) multi.train(training_features,labels) print multi.apply(testing_features) We are using a class as a base, but we can use any callable object: multi = one_against_rest(lambda : milk.supervised.tree(purity_measure='gini')) ... 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): labels, names = normaliselabels(labels) nclasses = labels.max() + 1 models = [] for i in xrange(nclasses): s = self.base() for k,v in self.options.iteritems(): s.set_option(k, v) model = s.train(features, labels == i) 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 == 1) if len(idxs) == 1: label = idxs[0] elif len(idxs) == 0: label = random.randint(0, self.nclasses - 1) else: label = random.choice(idxs) 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 be a callable object that provides the base classifier to use. Example ------- multi = one_against_one(milk.supervised.simple_svm) multi.train(training_features,labels) print multi.apply(testing_features) We are using a class as a base, but we can use any callable object: multi = one_against_one(lambda : milk.supervised.tree(purity_measure='gini')) ... 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): ''' one_against_one.train(objs,labels) ''' labels, names = normaliselabels(labels) 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): s = self.base() for k,v in self.options.iteritems(): s.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 = s.train(features[idxs],labels[idxs]==i) 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)] # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/supervised/tree.py0000644000175000017500000001247111277725552021004 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008, Luís 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. ''' ================ Tree Classifier ================ Decision tree based classifier ------------------------------ ''' from __future__ import division import numpy as np from collections import defaultdict import scipy.stats class _Leaf(object): def __init__(self, v, s): self.v = v/float(s) self.s = s def __repr__(self): return 'Leaf(%s,%s)' % (self.v, self.s) def _split(features, labels, criterion): N,f = features.shape best = None best_val = -1. for i in xrange(f): domain_i = np.unique(features[:,i]) for d in domain_i[1:]: cur_split = (features[:,i] < d) value = criterion(labels[cur_split],labels[~cur_split]) if value > best_val: best_val = value best = i,d return best def _entropy_set(labels): counts = defaultdict(float) for l in counts: counts[l] += 1. return scipy.stats.entropy(counts.values()) def information_gain(*args,**kwargs): ''' 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) ''' assert len(args) > 1, 'information_gain: only defined for two or more labels' H = 0. N = sum(len(A) for A in args) if kwargs.get('include_entropy',False): counts = defaultdict(float) for arg in args: for L in arg: counts[L] += 1. H = scipy.stats.entropy(counts.values()) for arg in args: H -= len(arg)/N * _entropy_set(arg) return H def build_tree(features, labels, criterion, min_split=4): ''' tree = build_tree(features, labels, criterion, min_split=4) Parameters ---------- * features: features to use * labels: labels * criterion: function to measure goodness of split * min_split: minimum size to split on ''' assert len(features) == len(labels), 'build_tree: Nr of labels does not match nr of features' features = np.asanyarray(features) labels = np.asanyarray(labels) if len(labels) < min_split: return _Leaf(labels.sum(), float(len(labels))) S = _split(features, labels, criterion) if S is None: return _Leaf(labels.sum(), float(len(labels))) i,thresh = S split = features[:,i] < thresh data0 = features[split] data1 = features[~split] labels0 = labels[split] labels1 = labels[~split] return (i, thresh, build_tree(data0, labels0, criterion, min_split), build_tree(data1, labels1, criterion, min_split)) 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 i,c,left,right = tree if features[i] < c: return apply_tree(left, features) return apply_tree(right, features) class tree_classifier(object): ''' tree = tree_classifier() tree.train(features,labels) predicted = tree.apply(testfeatures) A decision tree classifier (currently, implements the greedy ID3 algorithm without any pruning. Attributes ---------- * criterion: criterion to use for tree construction, this should be a function that receives a set of labels (default: information_gain). * min_split: minimum size to split on (default: 4). ''' def __init__(self, criterion=information_gain, min_split=4, return_label=True): self.criterion = criterion self.min_split = 4 self.return_label = return_label def train(self, features, labels): return tree_model(build_tree(features, labels, self.criterion, self.min_split), self.return_label) class tree_model(object): 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 # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/supervised/featureselection.py0000644000175000017500000001732211277667152023407 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2009, Luís 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 from numpy.linalg import det import scipy.stats from . classifier import normaliselabels _TOLERANCE = 0 _SIGNIFICANCE_IN = .15 _SIGNIFICANCE_OUT = .15 def _sweep(A, k, flag): N,_ = A.shape Akk = A[k,k] B = np.zeros_like(A) try: from scipy import weave from scipy.weave import converters k = int(k) code = ''' #line 42 "featureselection.py" for (int i = 0; i != N; ++i) { for (int j = 0; j != N; ++j) { if (i == k) { if (j == k) { B(i,j) = - 1./A(k,k); } else { B(i,j)=flag*A(i,j)/A(k,k); } } else if (j == k) { B(i,j)=flag*A(i,j)/A(k,k); } else { B(i,j)=A(i,j) - A(i,k)*A(k,j)/A(k,k); } } } ''' weave.inline( code, ['A','B','k','Akk','flag','N'], type_converters=converters.blitz) except: for i in xrange(N): for j in xrange(N): if i == k: if j == k: B[i,j] = -1./Akk else: B[i,j] = flag*A[i,j]/Akk elif j == k: B[i,j] = flag*A[i,j]/Akk else: B[i,j] = A[i,j] - A[i,k]*A[k,j]/Akk return B def sda(features,labels): ''' 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. ''' 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) mus = np.array([features[labels==i,:].mean(0) for i in xrange(q)]) mu = features.mean(0) W = np.zeros((m,m)) T = np.zeros((m,m)) try: from scipy import weave from scipy.weave import converters code=''' #line 106 "featureselection.py" for (int i = 0; i != m; ++i) { for (int j = 0; j != m; ++j) { for (int n = 0; n != N; ++n) { int g=labels(n); W(i,j) += (features(n,i)-mus(g,i))*(features(n,j)-mus(g,j)); T(i,j) += (features(n,i)-mu(i))*(features(n,j)-mu(j)); } } } ''' weave.inline( code, ['N','m','W','T','features','mu','mus','labels'], type_converters=converters.blitz) except ImportError: import warnings warnings.warn('scipy.weave failed. Resorting to (slow) Python code') for i in xrange(m): for j in xrange(m): for n in xrange(N): g = labels[n] W[i,j] += (features[n,i]-mus[g,i])*(features[n,j]-mus[g,j]) T[i,j] += (features[n,i]-mu[i])*(features[n,j]-mu[j]) ignoreidx = ( W.diagonal() == 0 ) if ignoreidx.any(): idxs, = np.where(~ignoreidx) F = sda(features[:,~ignoreidx],labels) return idxs[F] 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('pyslic.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 _rank(A,tol=1e-8): s = np.linalg.svd(A,compute_uv=0) return (s > tol).sum() def linear_independent_features(featmatrix, labels = None): ''' Returns the indices of a set of linearly independent features (columns). indices = linear_independent_features(features) ''' independent = [] R = _rank(featmatrix) i = 0 offset = 0 while i < featmatrix.shape[1]: R_ = _rank(np.delete(featmatrix,i,1)) if R_ == R: featmatrix = np.delete(featmatrix,i,1) offset += 1 else: independent.append(i+offset) i += 1 return np.array(independent) 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] 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): idxs = self.selector(features, labels) return filterfeatures(idxs) def sda_filter(): return featureselector(sda) milk-0.1-alpha-0/milk/supervised/classifier.py0000644000175000017500000000506211277674144022167 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008, Luís 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 __all__ = ['normaliselabels', 'ctransforms'] def normaliselabels(labels): ''' normalised, names = normaliselabels(labels) Normalises the labels to be integers from 0 through N-1 normalised is a numpy.array, while names is a dictionary mapping the indices to the old names. ''' labelnames={} normalised=[] names=[] N=0 for L in labels: nr=labelnames.get(L,N) if nr == N: labelnames[L]=N names.append(L) N += 1 normalised.append(nr) return numpy.array(normalised),names class ctransforms_model(object): 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): ''' ctransforms(c0, c1, c2, ...) Concatenate transforms. ''' def __init__(self,*args): self.transforms = args def train(self,features,labels): models = [] for T in self.transforms: model = T.train(features,labels) features = numpy.array([model.apply(f) for f in features]) models.append(model) return ctransforms_model(models) def set_option(self, opt, val): idx, opt = opt self.transforms[idx].set_option(opt,val) # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/supervised/_svm.cpp0000644000175000017500000007776011277121304021140 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 //using __gnu_debug::vector; using std::vector; 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_; std::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) { 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 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))); 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 (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) and (alpha2 < C) or (r2 > tol) and (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_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) { 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)); if (!(iter % 16)) check_for_interrupts(); // 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*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 = -G[t]; Gmax_idx = t; } } } else { if(!is_lower_bound(t)) { if(G[t] >= Gmax) { Gmax = G[t]; Gmax_idx = t; } } } } const int i = Gmax_idx; const double *Q_i = 0; const double *QDiag = cache_.get_diag(); if(i != -1) { // NULL Q_i not accessed: Gmax=-INF if i=-1 Q_i = cache_.get_kline(active_set[i],N); } for(int j=0;j= Gmax2) Gmax2 = G[j]; if (grad_diff > 0) { double obj_diff; double quad_coef=Q_i[active_set[i]]+QDiag[active_set[j]]-2*Y[i]*Q_i[active_set[j]]; if (quad_coef > 0) obj_diff = -(grad_diff*grad_diff)/quad_coef; else obj_diff = -(grad_diff*grad_diff)/tau; if (obj_diff <= obj_diff_min) { Gmin_idx=j; obj_diff_min = obj_diff; } } } } else { if (!is_upper_bound(j)) { const double grad_diff= Gmax-G[j]; if (-G[j] >= Gmax2) Gmax2 = -G[j]; if (grad_diff > 0) { double obj_diff; double quad_coef=Q_i[active_set[i]]+QDiag[active_set[j]]+2*Y[i]*Q_i[active_set[j]]; if (quad_coef > 0) obj_diff = -(grad_diff*grad_diff)/quad_coef; else obj_diff = -(grad_diff*grad_diff)/tau; if (obj_diff <= obj_diff_min) { Gmin_idx=j; obj_diff_min = obj_diff; } } } } } 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.1-alpha-0/milk/supervised/normalise.py0000644000175000017500000000735511277672654022047 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008, Luís 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 from collections import defaultdict import numpy as np try: import ncreduce _std = ncreduce.std except: _std = np.std __all__ = [ 'zscore', 'zscore_normalise', 'interval_normalise', 'chkfinite', ] def zscore(features): """ features = zscore(features) Returns a copy of features which has been normalised to zscores """ mu = features.mean(0) sigma = _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 class zscore_normalise(object): ''' Normalise to z-scores A preprocessor that normalises features to z scores. ''' def train(self,features,labels): shift = features.mean(0) factor = _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): D = features.ptp(0) shift = features.mean(0) + D/2. factor = D/2. return subtract_divide_model(shift, factor) 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 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): return self def apply(self, features): nans = np.isnan(features) + np.isinf(features) if nans.any(): features = features.copy() features[nans] = 0 return features # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/__init__.py0000644000175000017500000000242511256471214017377 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2009, Luís 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 measures.nfoldcrossvalidation as nfoldcrossvalidation # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/utils/0000755000175000017500000000000011302613311016407 5ustar luispedroluispedromilk-0.1-alpha-0/milk/utils/utils.py0000644000175000017500000000456111231643006020134 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2009, Luís 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 import random __all__ = [ 'get_nprandom', 'get_pyrandom', ] def get_nprandom(R): ''' R = get_nprandom(R) Returns a numpy.RandomState from R @param 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 ''' if R is None: return np.random.mtrand._rand if type(R) == int: return np.random.RandomState(R) return R def get_pyrandom(R): ''' R = get_pyrandom(R) Returns a random.Random object based on R @param R can be one of: * None : Returns the default python Random object * integer : Uses it as seed for constructing a new random generator * RandomState : Uses it as to generate a seed for a new random generator ''' if R is None: return random if type(R) is int: return random.Random(R) if type(R) is np.random.RandomState: return random.Random(R.randint(2**30)) raise TypeError,"get_pyrandom() does not know how to handle type %s." % type(R) # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/utils/__init__.py0000644000175000017500000000002411231643006020521 0ustar luispedroluispedrofrom utils import * milk-0.1-alpha-0/milk/unsupervised/0000755000175000017500000000000011302613311020003 5ustar luispedroluispedromilk-0.1-alpha-0/milk/unsupervised/kmeans.py0000644000175000017500000002356211255745030021656 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2009, Luís 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 import numpy from numpy import array, zeros, sqrt, inf, empty from ..utils import get_pyrandom from .normalise import zscore __all__ = ['kmeans','repeated_kmeans'] def _euclidean2(fmatrix, x, output=None): if output is None: output = np.zeros(N) try: from scipy import weave from scipy.weave import converters N,q = fmatrix.shape output[:] = 0 code = ''' for (int i = 0; i != N; ++i) { for (int j = 0; j != q; ++j) { output(i) += (fmatrix(i,j) - x(j))*(fmatrix(i,j)-x(j)); } } ''' weave.inline( code, ['fmatrix','N','q','x','output'], type_converters=converters.blitz) return output except: output[:] = ((fmatrix - x)**2).sum(1) return output def _mahalabonis2(fmatrix, x, icov, output=None): diff = (fmatrix-x) icov = (icov) # The expression below seems to be faster than looping over the elements and summing return dot(diff,dot(icov,diff.T)).diagonal() def centroid_errors(fmatrix,assignments,centroids): errors = [] for k in xrange(len(centroids)): errors.extend(fmatrix[assignments == k] - centroids[k]) return array(errors) def residual_sum_squares(fmatrix,assignments,centroids,distance='euclidean',**kwargs): if distance != 'euclidean': raise NotImplemented, "residual_sum_squares only implemented for 'euclidean' distance" return (centroid_errors(fmatrix,assignments,centroids)**2).sum() def kmeans(fmatrix,K,distance='euclidean',max_iter=1000,R=None,**kwargs): ''' assignments, centroids = kmean(fmatrix, K, distance='euclidean', R=None, icov=None, covmat=None) K-Means Clustering Parameters ========== * distance can be 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: Maximum number of iteration ''' fmatrix = np.asanyarray(fmatrix) if distance == 'euclidean': distfunction=_euclidean2 elif distance == 'seuclidean': fmatrix = zscore(fmatrix) distfunction=_euclidean2 elif distance == 'mahalanobis': icov = kwargs.get('icov',None) if icov is None: covmat = kwargs.get('covmat',None) if covmat is None: covmat = cov(fmatrix.T) icov = linalg.inv(covmat) distfunction = lambda f,x, out: _mahalabonis2(f, x, icov, output=out) else: raise 'Distance argument unknown (%s)' % distance R = get_pyrandom(R) N,q = fmatrix.shape centroids = np.array(R.sample(fmatrix,K), fmatrix.dtype) prev = np.zeros(N,np.int32) assignments = np.zeros(N,np.int32) dists = np.zeros(N, fmatrix.dtype) ndists = np.zeros(N, fmatrix.dtype) for i in xrange(max_iter): assignments[:] = 0 computed = False if distfunction is _euclidean2: try: from scipy import weave from scipy.weave import converters if fmatrix.flags['C_CONTIGUOUS'] and fmatrix.dtype in (np.float32,np.double): if fmatrix.dtype == np.float32: type = 'float' else: type = 'double' code = ''' #line 129 "kmeans.py" for (int i = 0; i != N; ++i) { %(type)s dist = std::numeric_limits<%(type)s>::infinity(); %(type)s* cd = centroids; for (int k = 0; k != K; ++k) { %(type)s ndist = 0.0; %(type)s* fd = fmatrix + i * q; for (int d = 0; d != q; ++d) { ndist += (*fd-*cd)*(*fd-*cd); ++fd; ++cd; } if (ndist < dist) { assignments[i] = k; dist = ndist; } } } ''' % { 'type' : type } weave.inline(code, ['fmatrix','N','q','K','assignments','centroids'], headers=['']) else: code = ''' for (int i = 0; i != N; ++i) { float dist = std::numeric_limits::infinity(); for (int k = 0; k != K; ++k) { float ndist = 0.0; for (int d = 0; d != q; ++d) { ndist += (fmatrix(i,d) - centroids(k,d))*(fmatrix(i,d) - centroids(k,d)); } if (ndist < dist) { assignments(i) = k; dist = ndist; } } } ''' weave.inline( code, ['fmatrix','N','q','K','assignments','centroids'], type_converters=converters.blitz) computed = True except: pass if not computed: dists[:] = np.inf for ci,C in enumerate(centroids): ndists = distfunction(fmatrix, C, output=ndists) try: from scipy import weave from scipy.weave import converters code = ''' for (int i = 0; i != N; ++i) { if (ndists(i) < dists(i)) { assignments(i) = ci; dists(i) = ndists(i); } } ''' weave.inline( code, ['dists', 'ndists', 'N', 'assignments', 'ci'], type_converters=converters.blitz) except Exception, e: print 'scipy.weave.inline failed. Resorting to Python code (Exception was "%s")' % e better = (ndists < dists) assignments[better] = ci dists[better] = ndists[better] if np.all(assignments == prev): break try: from scipy import weave from scipy.weave import converters centroids[:] = 0 counts = np.zeros(K,np.uint32) code = ''' for (int i = 0; i != N; ++i) { int c = assignments(i); ++counts(c); for (int j = 0; j != q; ++j) { centroids(c,j) += fmatrix(i,j); } } for (int i = 0; i != K; ++i) { for (int j = 0; j != q; ++j) { centroids(i,j) /= counts(i); } } ''' weave.inline( code, ['fmatrix','centroids','N','q','K','assignments','counts'], type_converters=converters.blitz) except Exception, e: print 'scipy.weave.inline failed. Resorting to Python code (Exception was "%s")' % e centroids = array([fmatrix[assignments == C].mean(0) for C in xrange(K)]) 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 kmeans ''' if distance == 'seuclidean': fmatrix = zscore(fmatrix) distance = 'euclidean' if distance != 'euclidean': raise NotImplementedError, "repeated_kmeans is only implemented for 'euclidean' or 'seuclidean' distance" best=+inf for i in xrange(iterations): A,C = kmeans(fmatrix, k, distance, max_iter=max_iter, R=R,**kwargs) rss = residual_sum_squares(fmatrix,A,C,distance,**kwargs) if rss < best: Ab,Cb = A,C best = rss return Ab,Cb # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/unsupervised/gaussianmixture.py0000644000175000017500000001124111231643006023611 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2009, Luís 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 from numpy import log, pi, array from numpy.linalg import det, inv from kmeans import residual_sum_squares, centroid_errors import scipy __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 = cov(diff.T) else: covm = covs[k] if covm.shape == (): covm = mat([[covm]]) icov = mat(inv(covm)) diff = mat(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'): 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, "nr_parameters: cannot handle model '%s'" % model def _compute(type, fmatrix, assignements, centroids, model='one_variance', covs=None): N,q = fmatrix.shape k = len(centroids) log_like = log_likelihood(fmatrix,assignements,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,assignements,centroids,model='one_variance',covs=None): ''' B = BIC(fmatrix,assignements,centroids,model) Compute Bayesian Information Criterion model can be one of: * 'one_variance': All features share the same variance parameter sigma^2 * 'full_covariance': Estimate a full covariance matrix or use covs[i] for centroid[i] @see AIC ''' return _compute('BIC', fmatrix, assignements, centroids, model, covs) def AIC(fmatrix,assignements,centroids,model='one_variance',covs=None): ''' A = AIC(fmatrix,assignements,centroids,model) Compute Akaike Information Criterion model can be one of: * 'one_variance': All features share the same variance parameter sigma^2 * 'full_covariance': Estimate a full covariance matrix or use covs[i] for centroid[i] @see BIC ''' return _compute('AIC', fmatrix, assignements, centroids, model, covs) # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/unsupervised/__init__.py0000644000175000017500000000236411231643006022126 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2009, Luís 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 kmeans import kmeans,repeated_kmeans from gaussianmixture import * from pca import pca from nnmf import * milk-0.1-alpha-0/milk/unsupervised/nnmf/0000755000175000017500000000000011302613311020741 5ustar luispedroluispedromilk-0.1-alpha-0/milk/unsupervised/nnmf/lee_seung.py0000644000175000017500000000536211231643006023274 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2009, Luís 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 from numpy import dot from ...utils import get_nprandom import scipy __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 @param cost can be one of: 'norm2' : minimise || X - AS ||_2 'i-div' : minimise D(X||AS), where D is I-divergence (generalisation of K-L divergence) @param max_iter Maximum number of iterations @param tol tolerance Threshold for early exit (when the update factor is with tol of 1., the function exits) @param R random seed @see get_nprandom 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 # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/unsupervised/nnmf/__init__.py0000644000175000017500000000246511231643006023066 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2009, Luís 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.1-alpha-0/milk/unsupervised/nnmf/hoyer.py0000644000175000017500000000714411231643006022454 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2009, Luís 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 import scipy from ...utils import get_nprandom __all__ = ['hoyer_sparse_nnmf'] def _norm2(x): return (x**2).sum() def _solve_alpha(s,m,L2): sm = s-m s2 = _norm2(s) sm2 = _norm2(sm) m2 = _norm2(m) dot = (m*sm).sum() alpha = (-dot + np.sqrt(dot**2 - sm2*(m2-L2)))/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)/(n-Z.sum()) s = s - c*(~Z) def _L1for(s,x,L2): ''' Solve for L1 in s = [ sqrt(n) - L1/sqrt(L2)] / [sqrt(n) - 1] ''' L2 = np.sqrt(L2) sn = np.sqrt(len(x)) return L2*s*((sn-1)-1) 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. 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 fixW(): for i in xrange(r): col = W[:,i] L2 = _norm2(col) W[:,i] = _project(col,_L1for(sparsenessW,col,L2),L2) def fixH(): for i in xrange(r): row = H[i,:] L2 = _norm2(row) H[i,:] = _project(row,_L1for(sparsenessH,col,L2),L2) 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 # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/unsupervised/pca.py0000644000175000017500000000335311231643006021131 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2009, Luís 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 from numpy import linalg from . import normalise 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: data matrix * zscore: whether to normalise to zscores (default: True) ''' 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 # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/unsupervised/normalise.py0000644000175000017500000000300311231643006022347 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2009, Luís 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.1-alpha-0/milk/active/0000755000175000017500000000000011302613311016522 5ustar luispedroluispedromilk-0.1-alpha-0/milk/active/roymccallum.py0000644000175000017500000000341011300347741021432 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008-2009, Luís 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 def expected_risk(labeled, pool, base): ''' goodnesses = expected_risk(classifier, pool) Roy & McCallum ''' features, labels = labeled labels, _ = normaliselabels(labeled) basic = baseclassifier() erisk = [] for i,val in enumerate(pool): classifier = loopclassifier() nfeatures[-1] = val for li in xrange(nclasses): nlabels[-1] = li classifier.train(nfeatures, nlabels) def _erisk(i): return np.array([_erisk(i) for i in xrange(len(pool))]) # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/active/__init__.py0000644000175000017500000000000011231643005020625 0ustar luispedroluispedromilk-0.1-alpha-0/milk/active/uncertainty.py0000644000175000017500000000424011231643005021445 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008, Luís 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. ''' 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(classifier, pool): ''' entropies = entropy(classifier,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(classifier.apply(u)) for u in pool] def one_minus_max(classifier,pool): ''' oneminus = one_minus_max(classifier,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(classifier.apply(u)) for u in pool] # vim: set ts=4 sts=4 sw=4 expandtab smartindent: milk-0.1-alpha-0/milk/active/eimpact.py0000644000175000017500000000457111231643005020531 0ustar luispedroluispedro# -*- coding: utf-8 -*- # Copyright (C) 2008, Luís 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.1-alpha-0/setup.cfg0000644000175000017500000000007311302613311016134 0ustar luispedroluispedro[egg_info] tag_build = tag_date = 0 tag_svn_revision = 0