diff --git a/GPy/util/__init__.py b/GPy/util/__init__.py index db9b7362..7d71c83c 100644 --- a/GPy/util/__init__.py +++ b/GPy/util/__init__.py @@ -14,5 +14,16 @@ import visualize import decorators import classification import latent_space_visualizations +import maps + +try: + import sympy + _sympy_available = True + del sympy +except ImportError as e: + _sympy_available = False + +if _sympy_available: + import symbolic import netpbmfile diff --git a/GPy/util/block_matrices.py b/GPy/util/block_matrices.py new file mode 100644 index 00000000..8fd5f89d --- /dev/null +++ b/GPy/util/block_matrices.py @@ -0,0 +1,24 @@ +import numpy as np + +def get_blocks(A, blocksizes): + assert (A.shape[0]==A.shape[1]) and len(A.shape)==2, "can;t blockify this non-square matrix" + N = np.sum(blocksizes) + assert A.shape[0] == N, "bad blocksizes" + num_blocks = len(blocksizes) + B = np.empty(shape=(num_blocks, num_blocks), dtype=np.object) + count_i = 0 + for Bi, i in enumerate(blocksizes): + count_j = 0 + for Bj, j in enumerate(blocksizes): + B[Bi, Bj] = A[count_i:count_i + i, count_j : count_j + j] + count_j += j + count_i += i + return B + + + +if __name__=='__main__': + A = np.zeros((5,5)) + B = get_blocks(A,[2,3]) + B[0,0] += 7 + print B diff --git a/GPy/util/config.py b/GPy/util/config.py index 02796e0b..b0789fe0 100644 --- a/GPy/util/config.py +++ b/GPy/util/config.py @@ -8,8 +8,8 @@ config = ConfigParser.ConfigParser() home = os.getenv('HOME') or os.getenv('USERPROFILE') user_file = os.path.join(home,'.gpy_config.cfg') default_file = os.path.abspath(os.path.join(os.path.dirname( __file__ ), '..', 'gpy_config.cfg')) -print user_file, os.path.isfile(user_file) -print default_file, os.path.isfile(default_file) +# print user_file, os.path.isfile(user_file) +# print default_file, os.path.isfile(default_file) # 1. check if the user has a ~/.gpy_config.cfg if os.path.isfile(user_file): diff --git a/GPy/util/data_resources.json b/GPy/util/data_resources.json new file mode 100644 index 00000000..ca15bf2d --- /dev/null +++ b/GPy/util/data_resources.json @@ -0,0 +1,382 @@ +{ + "rogers_girolami_data":{ + "files":[ + [ + "firstcoursemldata.tar.gz" + ] + ], + "license":null, + "citation":"A First Course in Machine Learning. Simon Rogers and Mark Girolami: Chapman & Hall/CRC, ISBN-13: 978-1439824146", + "details":"Data from the textbook 'A First Course in Machine Learning'. Available from http://www.dcs.gla.ac.uk/~srogers/firstcourseml/.", + "urls":[ + "https://www.dropbox.com/sh/7p6tu1t29idgliq/_XqlH_3nt9/" + ], + "suffices":[ + [ + "?dl=1" + ] + ], + "size":21949154 + }, + "ankur_pose_data":{ + "files":[ + [ + "ankurDataPoseSilhouette.mat" + ] + ], + "citation":"3D Human Pose from Silhouettes by Relevance Vector Regression (In CVPR'04). A. Agarwal and B. Triggs.", + "license":null, + "urls":[ + "http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/ankur_pose_data/" + ], + "details":"Artificially generated data of silhouettes given poses. Note that the data does not display a left/right ambiguity because across the entire data set one of the arms sticks out more the the other, disambiguating the pose as to which way the individual is facing.", + "size":1 + }, + "osu_accad":{ + "files":[ + [ + "swagger1TXT.ZIP", + "handspring1TXT.ZIP", + "quickwalkTXT.ZIP", + "run1TXT.ZIP", + "sprintTXT.ZIP", + "dogwalkTXT.ZIP", + "camper_04TXT.ZIP", + "dance_KB3_TXT.ZIP", + "per20_TXT.ZIP", + "perTWO07_TXT.ZIP", + "perTWO13_TXT.ZIP", + "perTWO14_TXT.ZIP", + "perTWO15_TXT.ZIP", + "perTWO16_TXT.ZIP" + ], + [ + "connections.txt" + ] + ], + "license":"Data is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License (http://creativecommons.org/licenses/by-nc-sa/3.0/).", + "citation":"The Open Motion Data Project by The Ohio State University Advanced Computing Center for the Arts and Design, http://accad.osu.edu/research/mocap/mocap_data.htm.", + "details":"Motion capture data of different motions from the Open Motion Data Project at Ohio State University.", + "urls":[ + "http://accad.osu.edu/research/mocap/data/", + "http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/stick/" + ], + "size":15922790 + }, + "isomap_face_data":{ + "files":[ + [ + "face_data.mat" + ] + ], + "license":null, + "citation":"A Global Geometric Framework for Nonlinear Dimensionality Reduction, J. B. Tenenbaum, V. de Silva and J. C. Langford, Science 290 (5500): 2319-2323, 22 December 2000", + "details":"Face data made available by Tenenbaum, de Silva and Langford to demonstrate isomap, available from http://isomap.stanford.edu/datasets.html.", + "urls":[ + "http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/isomap_face_data/" + ], + "size":24229368 + }, + "boston_housing":{ + "files":[ + [ + "Index", + "housing.data", + "housing.names" + ] + ], + "license":null, + "citation":"Harrison, D. and Rubinfeld, D.L. 'Hedonic prices and the demand for clean air', J. Environ. Economics & Management, vol.5, 81-102, 1978.", + "details":"The Boston Housing data relates house values in Boston to a range of input variables.", + "urls":[ + "http://archive.ics.uci.edu/ml/machine-learning-databases/housing/" + ], + "size":51276 + }, + "cmu_mocap_full":{ + "files":[ + [ + "allasfamc.zip" + ] + ], + "license":"From http://mocap.cs.cmu.edu. This data is free for use in research projects. You may include this data in commercially-sold products, but you may not resell this data directly, even in converted form. If you publish results obtained using this data, we would appreciate it if you would send the citation to your published paper to jkh+mocap@cs.cmu.edu, and also would add this text to your acknowledgments section: The data used in this project was obtained from mocap.cs.cmu.edu. The database was created with funding from NSF EIA-0196217.", + "citation":"Please include this in your acknowledgements: The data used in this project was obtained from mocap.cs.cmu.edu.\nThe database was created with funding from NSF EIA-0196217.", + "details":"CMU Motion Capture data base. Captured by a Vicon motion capture system consisting of 12 infrared MX-40 cameras, each of which is capable of recording at 120 Hz with images of 4 megapixel resolution. Motions are captured in a working volume of approximately 3m x 8m. The capture subject wears 41 markers and a stylish black garment.", + "urls":[ + "http://mocap.cs.cmu.edu/subjects" + ], + "size":null + }, + "brendan_faces":{ + "files":[ + [ + "frey_rawface.mat" + ] + ], + "license":null, + "citation":"Frey, B. J., Colmenarez, A and Huang, T. S. Mixtures of Local Linear Subspaces for Face Recognition. Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition 1998, 32-37, June 1998. Computer Society Press, Los Alamitos, CA.", + "details":"A video of Brendan Frey's face popularized as a benchmark for visualization by the Locally Linear Embedding.", + "urls":[ + "http://www.cs.nyu.edu/~roweis/data/" + ], + "size":1100584 + }, + "olympic_marathon_men":{ + "files":[ + [ + "olympicMarathonTimes.csv" + ] + ], + "license":null, + "citation":null, + "details":"Olympic mens' marathon gold medal winning times from 1896 to 2012. Time given in pace (minutes per kilometer). Data is originally downloaded and collated from Wikipedia, we are not responsible for errors in the data", + "urls":[ + "http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/olympic_marathon_men/" + ], + "size":584 + }, + "pumadyn-32nm":{ + "files":[ + [ + "pumadyn-32nm.tar.gz" + ] + ], + "license":"Data is made available by the Delve system at the University of Toronto", + "citation":"Created by Zoubin Ghahramani using the Matlab Robotics Toolbox of Peter Corke. Corke, P. I. (1996). A Robotics Toolbox for MATLAB. IEEE Robotics and Automation Magazine, 3 (1): 24-32.", + "details":"Pumadyn non linear 32 input data set with moderate noise. See http://www.cs.utoronto.ca/~delve/data/pumadyn/desc.html for details.", + "urls":[ + "ftp://ftp.cs.toronto.edu/pub/neuron/delve/data/tarfiles/pumadyn-family/" + ], + "size":5861646 + }, + "ripley_prnn_data":{ + "files":[ + [ + "Cushings.dat", + "README", + "crabs.dat", + "fglass.dat", + "fglass.grp", + "pima.te", + "pima.tr", + "pima.tr2", + "synth.te", + "synth.tr", + "viruses.dat", + "virus3.dat" + ] + ], + "license":null, + "citation":"Pattern Recognition and Neural Networks by B.D. Ripley (1996) Cambridge University Press ISBN 0 521 46986 7", + "details":"Data sets from Brian Ripley's Pattern Recognition and Neural Networks", + "urls":[ + "http://www.stats.ox.ac.uk/pub/PRNN/" + ], + "size":93565 + }, + "three_phase_oil_flow":{ + "files":[ + [ + "DataTrnLbls.txt", + "DataTrn.txt", + "DataTst.txt", + "DataTstLbls.txt", + "DataVdn.txt", + "DataVdnLbls.txt" + ] + ], + "license":null, + "citation":"Bishop, C. M. and G. D. James (1993). Analysis of multiphase flows using dual-energy gamma densitometry and neural networks. Nuclear Instruments and Methods in Physics Research A327, 580-593", + "details":"The three phase oil data used initially for demonstrating the Generative Topographic mapping.", + "urls":[ + "http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/three_phase_oil_flow/" + ], + "size":712796 + }, + "robot_wireless":{ + "files":[ + [ + "uw-floor.txt" + ] + ], + "license":null, + "citation":"WiFi-SLAM using Gaussian Process Latent Variable Models by Brian Ferris, Dieter Fox and Neil Lawrence in IJCAI'07 Proceedings pages 2480-2485. Data used in A Unifying Probabilistic Perspective for Spectral Dimensionality Reduction: Insights and New Models by Neil D. Lawrence, JMLR 13 pg 1609--1638, 2012.", + "details":"Data created by Brian Ferris and Dieter Fox. Consists of WiFi access point strengths taken during a circuit of the Paul Allen building at the University of Washington.", + "urls":[ + "http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/robot_wireless/" + ], + "size":284390 + }, + "xw_pen":{ + "files":[ + [ + "xw_pen_15.csv" + ] + ], + "license":null, + "citation":"Michael E. Tipping and Neil D. Lawrence. Variational inference for Student-t models: Robust Bayesian interpolation and generalised component analysis. Neurocomputing, 69:123--141, 2005", + "details":"Accelerometer pen data used for robust regression by Tipping and Lawrence.", + "urls":[ + "http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/xw_pen/" + ], + "size":3410 + }, + "swiss_roll":{ + "files":[ + [ + "swiss_roll_data.mat" + ] + ], + "license":null, + "citation":"A Global Geometric Framework for Nonlinear Dimensionality Reduction, J. B. Tenenbaum, V. de Silva and J. C. Langford, Science 290 (5500): 2319-2323, 22 December 2000", + "details":"Swiss roll data made available by Tenenbaum, de Silva and Langford to demonstrate isomap, available from http://isomap.stanford.edu/datasets.html.", + "urls":[ + "http://isomap.stanford.edu/" + ], + "size":800256 + }, + "osu_run1":{ + "files":[ + [ + "run1TXT.ZIP" + ], + [ + "connections.txt" + ] + ], + "license":"Data is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License (http://creativecommons.org/licenses/by-nc-sa/3.0/).", + "citation":"The Open Motion Data Project by The Ohio State University Advanced Computing Center for the Arts and Design, http://accad.osu.edu/research/mocap/mocap_data.htm.", + "details":"Motion capture data of a stick man running from the Open Motion Data Project at Ohio State University.", + "urls":[ + "http://accad.osu.edu/research/mocap/data/", + "http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/stick/" + ], + "size":338103 + }, + "creep_rupture":{ + "files":[ + [ + "creeprupt.tar" + ] + ], + "license":null, + "citation":"Materials Algorithms Project Data Library: MAP_DATA_CREEP_RUPTURE. F. Brun and T. Yoshida.", + "details":"Provides 2066 creep rupture test results of steels (mainly of two kinds of steels: 2.25Cr and 9-12 wt% Cr ferritic steels). See http://www.msm.cam.ac.uk/map/data/materials/creeprupt-b.html.", + "urls":[ + "http://www.msm.cam.ac.uk/map/data/tar/" + ], + "size":602797 + }, + "olivetti_faces":{ + "files":[ + [ + "att_faces.zip" + ], + [ + "olivettifaces.mat" + ] + ], + "license":null, + "citation":"Ferdinando Samaria and Andy Harter, Parameterisation of a Stochastic Model for Human Face Identification. Proceedings of 2nd IEEE Workshop on Applications of Computer Vision, Sarasota FL, December 1994", + "details":"Olivetti Research Labs Face data base, acquired between December 1992 and December 1994 in the Olivetti Research Lab, Cambridge (which later became AT&T Laboratories, Cambridge). When using these images please give credit to AT&T Laboratories, Cambridge. ", + "urls":[ + "http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/olivetti_faces/", + "http://www.cs.nyu.edu/~roweis/data/" + ], + "size":8561331 + }, + "olivetti_glasses":{ + "files":[ + [ + "has_glasses.np" + ], + [ + "olivettifaces.mat" + ] + ], + "license":null, + "citation":"Information recorded in olivetti_faces entry. Should be used from there.", + "details":"Information recorded in olivetti_faces entry. Should be used from there.", + "urls":[ + "http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/olivetti_faces/", + "http://www.cs.nyu.edu/~roweis/data/" + ], + "size":4261047 + }, + "della_gatta":{ + "files":[ + [ + "DellaGattadata.mat" + ] + ], + "license":null, + "citation":"Direct targets of the TRP63 transcription factor revealed by a combination of gene expression profiling and reverse engineering. Giusy Della Gatta, Mukesh Bansal, Alberto Ambesi-Impiombato, Dario Antonini, Caterina Missero, and Diego di Bernardo, Genome Research 2008", + "details":"The full gene expression data set from della Gatta et al (http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2413161/) processed by RMA.", + "urls":[ + "http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/della_gatta/" + ], + "size":3729650 + }, + "epomeo_gpx":{ + "files":[ + [ + "endomondo_1.gpx", + "endomondo_2.gpx", + "garmin_watch_via_endomondo.gpx", + "viewranger_phone.gpx", + "viewranger_tablet.gpx" + ] + ], + "license":null, + "citation":"", + "details":"Five different GPS traces of the same run up Mount Epomeo in Ischia. The traces are from different sources. endomondo_1 and endomondo_2 are traces from the mobile phone app Endomondo, with a split in the middle. garmin_watch_via_endomondo is the trace from a Garmin watch, with a segment missing about 4 kilometers in. viewranger_phone and viewranger_tablet are traces from a phone and a tablet through the viewranger app. The viewranger_phone data comes from the same mobile phone as the Endomondo data (i.e. there are 3 GPS devices, but one device recorded two traces).", + "urls":[ + "http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/epomeo_gpx/" + ], + "size":2031872 + }, + "mauna_loa":{ + "files":[ + [ + "co2_mm_mlo.txt" + ] + ], + "license":"-------------------------------------------------------------------- USE OF NOAA ESRL DATA\n\n These data are made freely available to the public and the scientific community in the belief that their wide dissemination will lead to greater understanding and new scientific insights. The availability of these data does not constitute publication of the data. NOAA relies on the ethics and integrity of the user to insure that ESRL receives fair credit for their work. If the data are obtained for potential use in a publication or presentation, ESRL should be informed at the outset of the nature of this work. If the ESRL data are essential to the work, or if an important result or conclusion depends on the ESRL data, co-authorship may be appropriate. This should be discussed at an early stage in the work. Manuscripts using the ESRL data should be sent to ESRL for review before they are submitted for publication so we can insure that the quality and limitations of the data are accurately represented.\n\n Contact: Pieter Tans (303 497 6678; pieter.tans@noaa.gov)\n\n RECIPROCITY Use of these data implies an agreement to reciprocate. Laboratories making similar measurements agree to make their own data available to the general public and to the scientific community in an equally complete and easily accessible form. Modelers are encouraged to make available to the community, upon request, their own tools used in the interpretation of the ESRL data, namely well documented model code, transport fields, and additional information necessary for other scientists to repeat the work and to run modified versions. Model availability includes collaborative support for new users of the models.\n --------------------------------------------------------------------\n\n See www.esrl.noaa.gov/gmd/ccgg/trends/ for additional details.", + "citation":"Mauna Loa Data. Dr. Pieter Tans, NOAA/ESRL (www.esrl.noaa.gov/gmd/ccgg/trends/) and Dr. Ralph Keeling, Scripps Institution of Oceanography (scrippsco2.ucsd.edu/).", + "details":"The 'average' column contains the monthly mean CO2 mole fraction determined from daily averages. The mole fraction of CO2, expressed as parts per million (ppm) is the number of molecules of CO2 in every one million molecules of dried air (water vapor removed). If there are missing days concentrated either early or late in the month, the monthly mean is corrected to the middle of the month using the average seasonal cycle. Missing months are denoted by -99.99. The 'interpolated' column includes average values from the preceding column and interpolated values where data are missing. Interpolated values are computed in two steps. First, we compute for each month the average seasonal cycle in a 7-year window around each monthly value. In this way the seasonal cycle is allowed to change slowly over time. We then determine the 'trend' value for each month by removing the seasonal cycle; this result is shown in the 'trend' column. Trend values are linearly interpolated for missing months. The interpolated monthly mean is then the sum of the average seasonal cycle value and the trend value for the missing month.\n\nNOTE: In general, the data presented for the last year are subject to change, depending on recalibration of the reference gas mixtures used, and other quality control procedures. Occasionally, earlier years may also be changed for the same reasons. Usually these changes are minor.\n\nCO2 expressed as a mole fraction in dry air, micromol/mol, abbreviated as ppm \n\n (-99.99 missing data; -1 no data for daily means in month)", + "urls":[ + "ftp://aftp.cmdl.noaa.gov/products/trends/co2/" + ], + "size":46779 + }, + "boxjenkins_airline":{ + "files":[ + [ + "boxjenkins_airline.csv" + ] + ], + "license":"You may copy and redistribute the data. You may make derivative works from the data. You may use the data for commercial purposes. You may not sublicence the data when redistributing it. You may not redistribute the data under a different license. Source attribution on any use of this data: Must refer source.", + "citation":"Box & Jenkins (1976), in file: data/airpass, Description: International airline passengers: monthly totals in thousands. Jan 49 – Dec 60", + "details":"International airline passengers, monthly totals from January 1949 to December 1960.", + "urls":[ + "http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/boxjenkins_airline/" + ], + "size":46779 + }, + + "decampos_characters":{ + "files":[ + [ + "characters.npy", + "digits.npy" + ] + ], + "license":null, + "citation":"T. de Campos, B. R. Babu, and M. Varma. Character recognition in natural images. VISAPP 2009.", + "details":"Examples of hand written digits taken from the de Campos et al paper on Character Recognition in Natural Images.", + "urls":[ + "http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/decampos_digits/" + ], + "size":2031872 + } +} diff --git a/GPy/util/datasets.py b/GPy/util/datasets.py index 2d10645c..059a39c3 100644 --- a/GPy/util/datasets.py +++ b/GPy/util/datasets.py @@ -3,11 +3,10 @@ import numpy as np import GPy import scipy.io import cPickle as pickle -import urllib as url import zipfile import tarfile import datetime - +import json ipython_available=True try: import IPython @@ -15,137 +14,28 @@ except ImportError: ipython_available=False -import sys, urllib +import sys, urllib2 -def reporthook(a,b,c): +def reporthook(a,b,c): # ',' at the end of the line is important! #print "% 3.1f%% of %d bytes\r" % (min(100, float(a * b) / c * 100), c), #you can also use sys.stdout.write sys.stdout.write("\r% 3.1f%% of %d bytes" % (min(100, float(a * b) / c * 100), c)) sys.stdout.flush() - + # Global variables data_path = os.path.join(os.path.dirname(__file__), 'datasets') default_seed = 10000 overide_manual_authorize=False neil_url = 'http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/' -sam_url = 'http://www.cs.nyu.edu/~roweis/data/' -cmu_url = 'http://mocap.cs.cmu.edu/subjects/' -# Note: there may be a better way of storing data resources, for the -# moment we are storing them in a dictionary. -data_resources = {'ankur_pose_data' : {'urls' : [neil_url + 'ankur_pose_data/'], - 'files' : [['ankurDataPoseSilhouette.mat']], - 'license' : None, - 'citation' : """3D Human Pose from Silhouettes by Relevance Vector Regression (In CVPR'04). A. Agarwal and B. Triggs.""", - 'details' : """Artificially generated data of silhouettes given poses. Note that the data does not display a left/right ambiguity because across the entire data set one of the arms sticks out more the the other, disambiguating the pose as to which way the individual is facing."""}, - - 'boston_housing' : {'urls' : ['http://archive.ics.uci.edu/ml/machine-learning-databases/housing/'], - 'files' : [['Index', 'housing.data', 'housing.names']], - 'citation' : """Harrison, D. and Rubinfeld, D.L. 'Hedonic prices and the demand for clean air', J. Environ. Economics & Management, vol.5, 81-102, 1978.""", - 'details' : """The Boston Housing data relates house values in Boston to a range of input variables.""", - 'license' : None, - 'size' : 51276 - }, - 'brendan_faces' : {'urls' : [sam_url], - 'files': [['frey_rawface.mat']], - 'citation' : 'Frey, B. J., Colmenarez, A and Huang, T. S. Mixtures of Local Linear Subspaces for Face Recognition. Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition 1998, 32-37, June 1998. Computer Society Press, Los Alamitos, CA.', - 'details' : """A video of Brendan Frey's face popularized as a benchmark for visualization by the Locally Linear Embedding.""", - 'license': None, - 'size' : 1100584}, - 'cmu_mocap_full' : {'urls' : ['http://mocap.cs.cmu.edu'], - 'files' : [['allasfamc.zip']], - 'citation' : """Please include this in your acknowledgements: The data used in this project was obtained from mocap.cs.cmu.edu. -The database was created with funding from NSF EIA-0196217.""", - 'details' : """CMU Motion Capture data base. Captured by a Vicon motion capture system consisting of 12 infrared MX-40 cameras, each of which is capable of recording at 120 Hz with images of 4 megapixel resolution. Motions are captured in a working volume of approximately 3m x 8m. The capture subject wears 41 markers and a stylish black garment.""", - 'license' : """From http://mocap.cs.cmu.edu. This data is free for use in research projects. You may include this data in commercially-sold products, but you may not resell this data directly, even in converted form. If you publish results obtained using this data, we would appreciate it if you would send the citation to your published paper to jkh+mocap@cs.cmu.edu, and also would add this text to your acknowledgments section: The data used in this project was obtained from mocap.cs.cmu.edu. The database was created with funding from NSF EIA-0196217.""", - 'size' : None}, - 'creep_rupture' : {'urls' : ['http://www.msm.cam.ac.uk/map/data/tar/'], - 'files' : [['creeprupt.tar']], - 'citation' : 'Materials Algorithms Project Data Library: MAP_DATA_CREEP_RUPTURE. F. Brun and T. Yoshida.', - 'details' : """Provides 2066 creep rupture test results of steels (mainly of two kinds of steels: 2.25Cr and 9-12 wt% Cr ferritic steels). See http://www.msm.cam.ac.uk/map/data/materials/creeprupt-b.html.""", - 'license' : None, - 'size' : 602797}, - 'della_gatta' : {'urls' : [neil_url + 'della_gatta/'], - 'files': [['DellaGattadata.mat']], - 'citation' : 'Direct targets of the TRP63 transcription factor revealed by a combination of gene expression profiling and reverse engineering. Giusy Della Gatta, Mukesh Bansal, Alberto Ambesi-Impiombato, Dario Antonini, Caterina Missero, and Diego di Bernardo, Genome Research 2008', - 'details': "The full gene expression data set from della Gatta et al (http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2413161/) processed by RMA.", - 'license':None, - 'size':3729650}, - 'epomeo_gpx' : {'urls' : [neil_url + 'epomeo_gpx/'], - 'files': [['endomondo_1.gpx', 'endomondo_2.gpx', 'garmin_watch_via_endomondo.gpx','viewranger_phone.gpx','viewranger_tablet.gpx']], - 'citation' : '', - 'details': "Five different GPS traces of the same run up Mount Epomeo in Ischia. The traces are from different sources. endomondo_1 and endomondo_2 are traces from the mobile phone app Endomondo, with a split in the middle. garmin_watch_via_endomondo is the trace from a Garmin watch, with a segment missing about 4 kilometers in. viewranger_phone and viewranger_tablet are traces from a phone and a tablet through the viewranger app. The viewranger_phone data comes from the same mobile phone as the Endomondo data (i.e. there are 3 GPS devices, but one device recorded two traces).", - 'license':None, - 'size': 2031872}, - 'three_phase_oil_flow': {'urls' : [neil_url + 'three_phase_oil_flow/'], - 'files' : [['DataTrnLbls.txt', 'DataTrn.txt', 'DataTst.txt', 'DataTstLbls.txt', 'DataVdn.txt', 'DataVdnLbls.txt']], - 'citation' : 'Bishop, C. M. and G. D. James (1993). Analysis of multiphase flows using dual-energy gamma densitometry and neural networks. Nuclear Instruments and Methods in Physics Research A327, 580-593', - 'details' : """The three phase oil data used initially for demonstrating the Generative Topographic mapping.""", - 'license' : None, - 'size' : 712796}, - 'rogers_girolami_data' : {'urls' : ['https://www.dropbox.com/sh/7p6tu1t29idgliq/_XqlH_3nt9/'], - 'files' : [['firstcoursemldata.tar.gz']], - 'suffices' : [['?dl=1']], - 'citation' : 'A First Course in Machine Learning. Simon Rogers and Mark Girolami: Chapman & Hall/CRC, ISBN-13: 978-1439824146', - 'details' : """Data from the textbook 'A First Course in Machine Learning'. Available from http://www.dcs.gla.ac.uk/~srogers/firstcourseml/.""", - 'license' : None, - 'size' : 21949154}, - 'olivetti_faces' : {'urls' : [neil_url + 'olivetti_faces/', sam_url], - 'files' : [['att_faces.zip'], ['olivettifaces.mat']], - 'citation' : 'Ferdinando Samaria and Andy Harter, Parameterisation of a Stochastic Model for Human Face Identification. Proceedings of 2nd IEEE Workshop on Applications of Computer Vision, Sarasota FL, December 1994', - 'details' : """Olivetti Research Labs Face data base, acquired between December 1992 and December 1994 in the Olivetti Research Lab, Cambridge (which later became AT&T Laboratories, Cambridge). When using these images please give credit to AT&T Laboratories, Cambridge. """, - 'license': None, - 'size' : 8561331}, - 'olympic_marathon_men' : {'urls' : [neil_url + 'olympic_marathon_men/'], - 'files' : [['olympicMarathonTimes.csv']], - 'citation' : None, - 'details' : """Olympic mens' marathon gold medal winning times from 1896 to 2012. Time given in pace (minutes per kilometer). Data is originally downloaded and collated from Wikipedia, we are not responsible for errors in the data""", - 'license': None, - 'size' : 584}, - 'osu_run1' : {'urls': ['http://accad.osu.edu/research/mocap/data/', neil_url + 'stick/'], - 'files': [['run1TXT.ZIP'],['connections.txt']], - 'details' : "Motion capture data of a stick man running from the Open Motion Data Project at Ohio State University.", - 'citation' : 'The Open Motion Data Project by The Ohio State University Advanced Computing Center for the Arts and Design, http://accad.osu.edu/research/mocap/mocap_data.htm.', - 'license' : 'Data is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License (http://creativecommons.org/licenses/by-nc-sa/3.0/).', - 'size': 338103}, - 'osu_accad' : {'urls': ['http://accad.osu.edu/research/mocap/data/', neil_url + 'stick/'], - 'files': [['swagger1TXT.ZIP','handspring1TXT.ZIP','quickwalkTXT.ZIP','run1TXT.ZIP','sprintTXT.ZIP','dogwalkTXT.ZIP','camper_04TXT.ZIP','dance_KB3_TXT.ZIP','per20_TXT.ZIP','perTWO07_TXT.ZIP','perTWO13_TXT.ZIP','perTWO14_TXT.ZIP','perTWO15_TXT.ZIP','perTWO16_TXT.ZIP'],['connections.txt']], - 'details' : "Motion capture data of different motions from the Open Motion Data Project at Ohio State University.", - 'citation' : 'The Open Motion Data Project by The Ohio State University Advanced Computing Center for the Arts and Design, http://accad.osu.edu/research/mocap/mocap_data.htm.', - 'license' : 'Data is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License (http://creativecommons.org/licenses/by-nc-sa/3.0/).', - 'size': 15922790}, - 'pumadyn-32nm' : {'urls' : ['ftp://ftp.cs.toronto.edu/pub/neuron/delve/data/tarfiles/pumadyn-family/'], - 'files' : [['pumadyn-32nm.tar.gz']], - 'details' : """Pumadyn non linear 32 input data set with moderate noise. See http://www.cs.utoronto.ca/~delve/data/pumadyn/desc.html for details.""", - 'citation' : """Created by Zoubin Ghahramani using the Matlab Robotics Toolbox of Peter Corke. Corke, P. I. (1996). A Robotics Toolbox for MATLAB. IEEE Robotics and Automation Magazine, 3 (1): 24-32.""", - 'license' : """Data is made available by the Delve system at the University of Toronto""", - 'size' : 5861646}, - 'robot_wireless' : {'urls' : [neil_url + 'robot_wireless/'], - 'files' : [['uw-floor.txt']], - 'citation' : """WiFi-SLAM using Gaussian Process Latent Variable Models by Brian Ferris, Dieter Fox and Neil Lawrence in IJCAI'07 Proceedings pages 2480-2485. Data used in A Unifying Probabilistic Perspective for Spectral Dimensionality Reduction: Insights and New Models by Neil D. Lawrence, JMLR 13 pg 1609--1638, 2012.""", - 'details' : """Data created by Brian Ferris and Dieter Fox. Consists of WiFi access point strengths taken during a circuit of the Paul Allen building at the University of Washington.""", - 'license' : None, - 'size' : 284390}, - 'swiss_roll' : {'urls' : ['http://isomap.stanford.edu/'], - 'files' : [['swiss_roll_data.mat']], - 'details' : """Swiss roll data made available by Tenenbaum, de Silva and Langford to demonstrate isomap, available from http://isomap.stanford.edu/datasets.html.""", - 'citation' : 'A Global Geometric Framework for Nonlinear Dimensionality Reduction, J. B. Tenenbaum, V. de Silva and J. C. Langford, Science 290 (5500): 2319-2323, 22 December 2000', - 'license' : None, - 'size' : 800256}, - 'ripley_prnn_data' : {'urls' : ['http://www.stats.ox.ac.uk/pub/PRNN/'], - 'files' : [['Cushings.dat', 'README', 'crabs.dat', 'fglass.dat', 'fglass.grp', 'pima.te', 'pima.tr', 'pima.tr2', 'synth.te', 'synth.tr', 'viruses.dat', 'virus3.dat']], - 'details' : """Data sets from Brian Ripley's Pattern Recognition and Neural Networks""", - 'citation': """Pattern Recognition and Neural Networks by B.D. Ripley (1996) Cambridge University Press ISBN 0 521 46986 7""", - 'license' : None, - 'size' : 93565}, - 'isomap_face_data' : {'urls' : [neil_url + 'isomap_face_data/'], - 'files' : [['face_data.mat']], - 'details' : """Face data made available by Tenenbaum, de Silva and Langford to demonstrate isomap, available from http://isomap.stanford.edu/datasets.html.""", - 'citation' : 'A Global Geometric Framework for Nonlinear Dimensionality Reduction, J. B. Tenenbaum, V. de Silva and J. C. Langford, Science 290 (5500): 2319-2323, 22 December 2000', - 'license' : None, - 'size' : 24229368}, - } +# Read data resources from json file. +# Don't do this when ReadTheDocs is scanning as it breaks things +on_rtd = os.environ.get('READTHEDOCS', None) == 'True' #Checks if RTD is scanning +if not (on_rtd): + path = os.path.join(os.path.dirname(__file__), 'data_resources.json') + json_data=open(path).read() + data_resources = json.loads(json_data) def prompt_user(prompt): @@ -158,14 +48,14 @@ def prompt_user(prompt): print(prompt) choice = raw_input().lower() # would like to test for exception here, but not sure if we can do that without importing IPython - except: + except: print('Stdin is not implemented.') print('You need to set') print('overide_manual_authorize=True') print('to proceed with the download. Please set that variable and continue.') raise - + if choice in yes: return True elif choice in no: @@ -183,7 +73,7 @@ def data_available(dataset_name=None): if not os.path.exists(os.path.join(data_path, dataset_name, file)): return False return True - + def download_url(url, store_directory, save_name = None, messages = True, suffix=''): """Download a file from a url and save it to disk.""" i = url.rfind('/') @@ -194,7 +84,21 @@ def download_url(url, store_directory, save_name = None, messages = True, suffix print "Downloading ", url, "->", os.path.join(store_directory, file) if not os.path.exists(dir_name): os.makedirs(dir_name) - urllib.urlretrieve(url+suffix, save_name, reporthook) + try: + response = urllib2.urlopen(url+suffix) + except urllib2.URLError, e: + if not hasattr(e, "code"): + raise + response = e + if response.code > 399 and response.code<500: + raise ValueError('Tried url ' + url + suffix + ' and received client error ' + str(response.code)) + elif response.code > 499: + raise ValueError('Tried url ' + url + suffix + ' and received server error ' + str(response.code)) + # if we wanted to get more sophisticated maybe we should check the response code here again even for successes. + with open(save_name, 'wb') as f: + f.write(response.read()) + + #urllib.urlretrieve(url+suffix, save_name, reporthook) def authorize_download(dataset_name=None): """Check with the user that the are happy with terms and conditions for the data set.""" @@ -243,18 +147,20 @@ def download_data(dataset_name=None): for file in files: download_url(os.path.join(url,file), dataset_name, dataset_name) return True - + def data_details_return(data, data_set): """Update the data component of the data dictionary with details drawn from the data_resources.""" data.update(data_resources[data_set]) return data - + def cmu_urls_files(subj_motions, messages = True): ''' - Find which resources are missing on the local disk for the requested CMU motion capture motions. + Find which resources are missing on the local disk for the requested CMU motion capture motions. ''' - + dr = data_resources['cmu_mocap_full'] + cmu_url = dr['urls'][0] + subjects_num = subj_motions[0] motions_num = subj_motions[1] @@ -274,15 +180,15 @@ def cmu_urls_files(subj_motions, messages = True): motions[i].append(curMot) all_skels = [] - + assert len(subjects) == len(motions) - + all_motions = [] - + for i in range(len(subjects)): skel_dir = os.path.join(data_path, 'cmu_mocap') cur_skel_file = os.path.join(skel_dir, subjects[i] + '.asf') - + url_required = False file_download = [] if not os.path.exists(cur_skel_file): @@ -299,7 +205,7 @@ def cmu_urls_files(subj_motions, messages = True): url_required = True file_download.append(subjects[i] + '_' + motions[i][j] + '.amc') if url_required: - resource['urls'].append(cmu_url + subjects[i] + '/') + resource['urls'].append(cmu_url + '/' + subjects[i] + '/') resource['files'].append(file_download) return resource @@ -326,10 +232,10 @@ if gpxpy_available: points = [point for track in gpx.tracks for segment in track.segments for point in segment.points] data = [[(point.time-datetime.datetime(2013,8,21)).total_seconds(), point.latitude, point.longitude, point.elevation] for point in points] X.append(np.asarray(data)[::sample_every, :]) - gpx_file.close() + gpx_file.close() return data_details_return({'X' : X, 'info' : 'Data is an array containing time in seconds, latitude, longitude and elevation in that order.'}, data_set) -del gpxpy_available +#del gpxpy_available @@ -402,7 +308,7 @@ def oil(data_set='three_phase_oil_flow'): return data_details_return({'X': X, 'Y': Y, 'Xtest': Xtest, 'Ytest': Ytest, 'Xtest' : Xtest, 'Xvalid': Xvalid, 'Yvalid': Yvalid}, data_set) #else: # throw an error - + def oil_100(seed=default_seed, data_set = 'three_phase_oil_flow'): np.random.seed(seed=seed) data = oil() @@ -489,6 +395,18 @@ def silhouette(data_set='ankur_pose_data'): Ytest = mat_data['Z_test'] return data_details_return({'X': X, 'Y': Y, 'Xtest': Xtest, 'Ytest': Ytest}, data_set) +def decampos_digits(data_set='decampos_characters', which_digits=[0,1,2,3,4,5,6,7,8,9]): + if not data_available(data_set): + download_data(data_set) + path = os.path.join(data_path, data_set) + digits = np.load(os.path.join(path, 'digits.npy')) + digits = digits[which_digits,:,:,:] + num_classes, num_samples, height, width = digits.shape + Y = digits.reshape((digits.shape[0]*digits.shape[1],digits.shape[2]*digits.shape[3])) + lbls = np.array([[l]*num_samples for l in which_digits]).reshape(Y.shape[0], 1) + str_lbls = np.array([[str(l)]*num_samples for l in which_digits]) + return data_details_return({'Y': Y, 'lbls': lbls, 'str_lbls' : str_lbls, 'info': 'Digits data set from the de Campos characters data'}, data_set) + def ripley_synth(data_set='ripley_prnn_data'): if not data_available(data_set): download_data(data_set) @@ -498,7 +416,36 @@ def ripley_synth(data_set='ripley_prnn_data'): test = np.genfromtxt(os.path.join(data_path, data_set, 'synth.te'), skip_header=1) Xtest = test[:, 0:2] ytest = test[:, 2:3] - return data_details_return({'X': X, 'y': y, 'Xtest': Xtest, 'ytest': ytest, 'info': 'Synthetic data generated by Ripley for a two class classification problem.'}, data_set) + return data_details_return({'X': X, 'Y': y, 'Xtest': Xtest, 'Ytest': ytest, 'info': 'Synthetic data generated by Ripley for a two class classification problem.'}, data_set) + +def mauna_loa(data_set='mauna_loa', num_train=543, refresh_data=False): + path = os.path.join(data_path, data_set) + if data_available(data_set) and not refresh_data: + print 'Using cached version of the data set, to use latest version set refresh_data to True' + else: + download_data(data_set) + data = np.loadtxt(os.path.join(data_path, data_set, 'co2_mm_mlo.txt')) + print 'Most recent data observation from month ', data[-1, 1], ' in year ', data[-1, 0] + allX = data[data[:, 3]!=-99.99, 2:3] + allY = data[data[:, 3]!=-99.99, 3:4] + X = allX[:num_train, 0:1] + Xtest = allX[num_train:, 0:1] + Y = allY[:num_train, 0:1] + Ytest = allY[num_train:, 0:1] + return data_details_return({'X': X, 'Y': Y, 'Xtest': Xtest, 'Ytest': Ytest, 'info': "Mauna Loa data with " + str(num_train) + " values used as training points."}, data_set) + + +def boxjenkins_airline(data_set='boxjenkins_airline', num_train=96): + path = os.path.join(data_path, data_set) + if not data_available(data_set): + download_data(data_set) + data = np.loadtxt(os.path.join(data_path, data_set, 'boxjenkins_airline.csv'), delimiter=',') + Y = data[:num_train, 1:2] + X = data[:num_train, 0:1] + Xtest = data[num_train:, 0:1] + Ytest = data[num_train:, 1:2] + return data_details_return({'X': X, 'Y': Y, 'Xtest': Xtest, 'Ytest': Ytest, 'info': "Montly airline passenger data from Box & Jenkins 1976."}, data_set) + def osu_run1(data_set='osu_run1', sample_every=4): path = os.path.join(data_path, data_set) @@ -547,7 +494,7 @@ def simulation_BGPLVM(): Y = np.array(mat_data['Y'], dtype=float) S = np.array(mat_data['initS'], dtype=float) mu = np.array(mat_data['initMu'], dtype=float) - return data_details_return({'S': S, 'Y': Y, 'mu': mu}, mat_data) + #return data_details_return({'S': S, 'Y': Y, 'mu': mu}, data_set) return {'Y': Y, 'S': S, 'mu' : mu, 'info': "Simulated test dataset generated in MATLAB to compare BGPLVM between python and MATLAB"} @@ -591,6 +538,21 @@ def toy_linear_1d_classification(seed=default_seed): X = (np.r_[x1, x2])[:, None] return {'X': X, 'Y': sample_class(2.*X), 'F': 2.*X, 'seed' : seed} +def olivetti_glasses(data_set='olivetti_glasses', num_training=200, seed=default_seed): + path = os.path.join(data_path, data_set) + if not data_available(data_set): + download_data(data_set) + y = np.load(os.path.join(path, 'has_glasses.np')) + y = np.where(y=='y',1,0).reshape(-1,1) + faces = scipy.io.loadmat(os.path.join(path, 'olivettifaces.mat'))['faces'].T + np.random.seed(seed=seed) + index = np.random.permutation(faces.shape[0]) + X = faces[index[:num_training],:] + Xtest = faces[index[num_training:],:] + Y = y[index[:num_training],:] + Ytest = y[index[num_training:]] + return data_details_return({'X': X, 'Y': Y, 'Xtest': Xtest, 'Ytest': Ytest, 'seed' : seed, 'info': "ORL Faces with labels identifiying who is wearing glasses and who isn't. Data is randomly partitioned according to given seed. Presence or absence of glasses was labelled by James Hensman."}, 'olivetti_faces') + def olivetti_faces(data_set='olivetti_faces'): path = os.path.join(data_path, data_set) if not data_available(data_set): @@ -608,8 +570,16 @@ def olivetti_faces(data_set='olivetti_faces'): Y = np.asarray(Y) lbls = np.asarray(lbls)[:, None] return data_details_return({'Y': Y, 'lbls' : lbls, 'info': "ORL Faces processed to 64x64 images."}, data_set) - -def download_rogers_girolami_data(): + +def xw_pen(data_set='xw_pen'): + if not data_available(data_set): + download_data(data_set) + Y = np.loadtxt(os.path.join(data_path, data_set, 'xw_pen_15.csv'), delimiter=',') + X = np.arange(485)[:, None] + return data_details_return({'Y': Y, 'X': X, 'info': "Tilt data from a personalized digital assistant pen. Plot in original paper showed regression between time steps 175 and 275."}, data_set) + + +def download_rogers_girolami_data(data_set='rogers_girolami_data'): if not data_available('rogers_girolami_data'): download_data(data_set) path = os.path.join(data_path, data_set) @@ -675,7 +645,7 @@ def olympic_marathon_men(data_set='olympic_marathon_men'): Y = olympics[:, 1:2] return data_details_return({'X': X, 'Y': Y}, data_set) -def olympics(): +def olympic_sprints(data_set='rogers_girolami_data'): """All olympics sprint winning times for multiple output prediction.""" X = np.zeros((0, 2)) Y = np.zeros((0, 1)) @@ -693,7 +663,18 @@ def olympics(): data['X'] = X data['Y'] = Y data['info'] = "Olympics sprint event winning for men and women to 2008. Data is from Rogers and Girolami's First Course in Machine Learning." - return data + return data_details_return({ + 'X': X, + 'Y': Y, + 'info': "Olympics sprint event winning for men and women to 2008. Data is from Rogers and Girolami's First Course in Machine Learning.", + 'output_info': { + 0:'100m Men', + 1:'100m Women', + 2:'200m Men', + 3:'200m Women', + 4:'400m Men', + 5:'400m Women'} + }, data_set) # def movielens_small(partNo=1,seed=default_seed): # np.random.seed(seed=seed) @@ -786,15 +767,15 @@ def creep_data(data_set='creep_rupture'): X = all_data[:, features].copy() return data_details_return({'X': X, 'y': y}, data_set) -def cmu_mocap_49_balance(): +def cmu_mocap_49_balance(data_set='cmu_mocap'): """Load CMU subject 49's one legged balancing motion that was used by Alvarez, Luengo and Lawrence at AISTATS 2009.""" train_motions = ['18', '19'] test_motions = ['20'] - data = cmu_mocap('49', train_motions, test_motions, sample_every=4) + data = cmu_mocap('49', train_motions, test_motions, sample_every=4, data_set=data_set) data['info'] = "One legged balancing motions from CMU data base subject 49. As used in Alvarez, Luengo and Lawrence at AISTATS 2009. It consists of " + data['info'] return data -def cmu_mocap_35_walk_jog(): +def cmu_mocap_35_walk_jog(data_set='cmu_mocap'): """Load CMU subject 35's walking and jogging motions, the same data that was used by Taylor, Roweis and Hinton at NIPS 2007. but without their preprocessing. Also used by Lawrence at AISTATS 2007.""" train_motions = ['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12', @@ -802,7 +783,7 @@ def cmu_mocap_35_walk_jog(): '20', '21', '22', '23', '24', '25', '26', '28', '30', '31', '32', '33', '34'] test_motions = ['18', '29'] - data = cmu_mocap('35', train_motions, test_motions, sample_every=4) + data = cmu_mocap('35', train_motions, test_motions, sample_every=4, data_set=data_set) data['info'] = "Walk and jog data from CMU data base subject 35. As used in Tayor, Roweis and Hinton at NIPS 2007, but without their pre-processing (i.e. as used by Lawrence at AISTATS 2007). It consists of " + data['info'] return data @@ -814,7 +795,7 @@ def cmu_mocap(subject, train_motions, test_motions=[], sample_every=4, data_set= # Make sure the data is downloaded. all_motions = train_motions + test_motions resource = cmu_urls_files(([subject], [all_motions])) - data_resources[data_set] = data_resources['cmu_mocap_full'] + data_resources[data_set] = data_resources['cmu_mocap_full'].copy() data_resources[data_set]['files'] = resource['files'] data_resources[data_set]['urls'] = resource['urls'] if resource['urls']: @@ -884,3 +865,5 @@ def cmu_mocap(subject, train_motions, test_motions=[], sample_every=4, data_set= if sample_every != 1: info += ' Data is sub-sampled to every ' + str(sample_every) + ' frames.' return data_details_return({'Y': Y, 'lbls' : lbls, 'Ytest': Ytest, 'lblstest' : lblstest, 'info': info, 'skel': skel}, data_set) + + diff --git a/GPy/util/datasets/data_resources_create.py b/GPy/util/datasets/data_resources_create.py new file mode 100644 index 00000000..8ae62a85 --- /dev/null +++ b/GPy/util/datasets/data_resources_create.py @@ -0,0 +1,127 @@ +import json + +neil_url = 'http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/' +sam_url = 'http://www.cs.nyu.edu/~roweis/data/' +cmu_url = 'http://mocap.cs.cmu.edu/subjects/' + +data_resources = {'ankur_pose_data' : {'urls' : [neil_url + 'ankur_pose_data/'], + 'files' : [['ankurDataPoseSilhouette.mat']], + 'license' : None, + 'citation' : """3D Human Pose from Silhouettes by Relevance Vector Regression (In CVPR'04). A. Agarwal and B. Triggs.""", + 'details' : """Artificially generated data of silhouettes given poses. Note that the data does not display a left/right ambiguity because across the entire data set one of the arms sticks out more the the other, disambiguating the pose as to which way the individual is facing."""}, + + 'boston_housing' : {'urls' : ['http://archive.ics.uci.edu/ml/machine-learning-databases/housing/'], + 'files' : [['Index', 'housing.data', 'housing.names']], + 'citation' : """Harrison, D. and Rubinfeld, D.L. 'Hedonic prices and the demand for clean air', J. Environ. Economics & Management, vol.5, 81-102, 1978.""", + 'details' : """The Boston Housing data relates house values in Boston to a range of input variables.""", + 'license' : None, + 'size' : 51276 + }, + 'brendan_faces' : {'urls' : [sam_url], + 'files': [['frey_rawface.mat']], + 'citation' : 'Frey, B. J., Colmenarez, A and Huang, T. S. Mixtures of Local Linear Subspaces for Face Recognition. Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition 1998, 32-37, June 1998. Computer Society Press, Los Alamitos, CA.', + 'details' : """A video of Brendan Frey's face popularized as a benchmark for visualization by the Locally Linear Embedding.""", + 'license': None, + 'size' : 1100584}, + 'cmu_mocap_full' : {'urls' : ['http://mocap.cs.cmu.edu'], + 'files' : [['allasfamc.zip']], + 'citation' : """Please include this in your acknowledgements: The data used in this project was obtained from mocap.cs.cmu.edu. +The database was created with funding from NSF EIA-0196217.""", + 'details' : """CMU Motion Capture data base. Captured by a Vicon motion capture system consisting of 12 infrared MX-40 cameras, each of which is capable of recording at 120 Hz with images of 4 megapixel resolution. Motions are captured in a working volume of approximately 3m x 8m. The capture subject wears 41 markers and a stylish black garment.""", + 'license' : """From http://mocap.cs.cmu.edu. This data is free for use in research projects. You may include this data in commercially-sold products, but you may not resell this data directly, even in converted form. If you publish results obtained using this data, we would appreciate it if you would send the citation to your published paper to jkh+mocap@cs.cmu.edu, and also would add this text to your acknowledgments section: The data used in this project was obtained from mocap.cs.cmu.edu. The database was created with funding from NSF EIA-0196217.""", + 'size' : None}, + 'creep_rupture' : {'urls' : ['http://www.msm.cam.ac.uk/map/data/tar/'], + 'files' : [['creeprupt.tar']], + 'citation' : 'Materials Algorithms Project Data Library: MAP_DATA_CREEP_RUPTURE. F. Brun and T. Yoshida.', + 'details' : """Provides 2066 creep rupture test results of steels (mainly of two kinds of steels: 2.25Cr and 9-12 wt% Cr ferritic steels). See http://www.msm.cam.ac.uk/map/data/materials/creeprupt-b.html.""", + 'license' : None, + 'size' : 602797}, + 'della_gatta' : {'urls' : [neil_url + 'della_gatta/'], + 'files': [['DellaGattadata.mat']], + 'citation' : 'Direct targets of the TRP63 transcription factor revealed by a combination of gene expression profiling and reverse engineering. Giusy Della Gatta, Mukesh Bansal, Alberto Ambesi-Impiombato, Dario Antonini, Caterina Missero, and Diego di Bernardo, Genome Research 2008', + 'details': "The full gene expression data set from della Gatta et al (http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2413161/) processed by RMA.", + 'license':None, + 'size':3729650}, + 'epomeo_gpx' : {'urls' : [neil_url + 'epomeo_gpx/'], + 'files': [['endomondo_1.gpx', 'endomondo_2.gpx', 'garmin_watch_via_endomondo.gpx','viewranger_phone.gpx','viewranger_tablet.gpx']], + 'citation' : '', + 'details': "Five different GPS traces of the same run up Mount Epomeo in Ischia. The traces are from different sources. endomondo_1 and endomondo_2 are traces from the mobile phone app Endomondo, with a split in the middle. garmin_watch_via_endomondo is the trace from a Garmin watch, with a segment missing about 4 kilometers in. viewranger_phone and viewranger_tablet are traces from a phone and a tablet through the viewranger app. The viewranger_phone data comes from the same mobile phone as the Endomondo data (i.e. there are 3 GPS devices, but one device recorded two traces).", + 'license':None, + 'size': 2031872}, + 'three_phase_oil_flow': {'urls' : [neil_url + 'three_phase_oil_flow/'], + 'files' : [['DataTrnLbls.txt', 'DataTrn.txt', 'DataTst.txt', 'DataTstLbls.txt', 'DataVdn.txt', 'DataVdnLbls.txt']], + 'citation' : 'Bishop, C. M. and G. D. James (1993). Analysis of multiphase flows using dual-energy gamma densitometry and neural networks. Nuclear Instruments and Methods in Physics Research A327, 580-593', + 'details' : """The three phase oil data used initially for demonstrating the Generative Topographic mapping.""", + 'license' : None, + 'size' : 712796}, + 'rogers_girolami_data' : {'urls' : ['https://www.dropbox.com/sh/7p6tu1t29idgliq/_XqlH_3nt9/'], + 'files' : [['firstcoursemldata.tar.gz']], + 'suffices' : [['?dl=1']], + 'citation' : 'A First Course in Machine Learning. Simon Rogers and Mark Girolami: Chapman & Hall/CRC, ISBN-13: 978-1439824146', + 'details' : """Data from the textbook 'A First Course in Machine Learning'. Available from http://www.dcs.gla.ac.uk/~srogers/firstcourseml/.""", + 'license' : None, + 'size' : 21949154}, + 'olivetti_faces' : {'urls' : [neil_url + 'olivetti_faces/', sam_url], + 'files' : [['att_faces.zip'], ['olivettifaces.mat']], + 'citation' : 'Ferdinando Samaria and Andy Harter, Parameterisation of a Stochastic Model for Human Face Identification. Proceedings of 2nd IEEE Workshop on Applications of Computer Vision, Sarasota FL, December 1994', + 'details' : """Olivetti Research Labs Face data base, acquired between December 1992 and December 1994 in the Olivetti Research Lab, Cambridge (which later became AT&T Laboratories, Cambridge). When using these images please give credit to AT&T Laboratories, Cambridge. """, + 'license': None, + 'size' : 8561331}, + 'olympic_marathon_men' : {'urls' : [neil_url + 'olympic_marathon_men/'], + 'files' : [['olympicMarathonTimes.csv']], + 'citation' : None, + 'details' : """Olympic mens' marathon gold medal winning times from 1896 to 2012. Time given in pace (minutes per kilometer). Data is originally downloaded and collated from Wikipedia, we are not responsible for errors in the data""", + 'license': None, + 'size' : 584}, + 'osu_run1' : {'urls': ['http://accad.osu.edu/research/mocap/data/', neil_url + 'stick/'], + 'files': [['run1TXT.ZIP'],['connections.txt']], + 'details' : "Motion capture data of a stick man running from the Open Motion Data Project at Ohio State University.", + 'citation' : 'The Open Motion Data Project by The Ohio State University Advanced Computing Center for the Arts and Design, http://accad.osu.edu/research/mocap/mocap_data.htm.', + 'license' : 'Data is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License (http://creativecommons.org/licenses/by-nc-sa/3.0/).', + 'size': 338103}, + 'osu_accad' : {'urls': ['http://accad.osu.edu/research/mocap/data/', neil_url + 'stick/'], + 'files': [['swagger1TXT.ZIP','handspring1TXT.ZIP','quickwalkTXT.ZIP','run1TXT.ZIP','sprintTXT.ZIP','dogwalkTXT.ZIP','camper_04TXT.ZIP','dance_KB3_TXT.ZIP','per20_TXT.ZIP','perTWO07_TXT.ZIP','perTWO13_TXT.ZIP','perTWO14_TXT.ZIP','perTWO15_TXT.ZIP','perTWO16_TXT.ZIP'],['connections.txt']], + 'details' : "Motion capture data of different motions from the Open Motion Data Project at Ohio State University.", + 'citation' : 'The Open Motion Data Project by The Ohio State University Advanced Computing Center for the Arts and Design, http://accad.osu.edu/research/mocap/mocap_data.htm.', + 'license' : 'Data is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License (http://creativecommons.org/licenses/by-nc-sa/3.0/).', + 'size': 15922790}, + 'pumadyn-32nm' : {'urls' : ['ftp://ftp.cs.toronto.edu/pub/neuron/delve/data/tarfiles/pumadyn-family/'], + 'files' : [['pumadyn-32nm.tar.gz']], + 'details' : """Pumadyn non linear 32 input data set with moderate noise. See http://www.cs.utoronto.ca/~delve/data/pumadyn/desc.html for details.""", + 'citation' : """Created by Zoubin Ghahramani using the Matlab Robotics Toolbox of Peter Corke. Corke, P. I. (1996). A Robotics Toolbox for MATLAB. IEEE Robotics and Automation Magazine, 3 (1): 24-32.""", + 'license' : """Data is made available by the Delve system at the University of Toronto""", + 'size' : 5861646}, + 'robot_wireless' : {'urls' : [neil_url + 'robot_wireless/'], + 'files' : [['uw-floor.txt']], + 'citation' : """WiFi-SLAM using Gaussian Process Latent Variable Models by Brian Ferris, Dieter Fox and Neil Lawrence in IJCAI'07 Proceedings pages 2480-2485. Data used in A Unifying Probabilistic Perspective for Spectral Dimensionality Reduction: Insights and New Models by Neil D. Lawrence, JMLR 13 pg 1609--1638, 2012.""", + 'details' : """Data created by Brian Ferris and Dieter Fox. Consists of WiFi access point strengths taken during a circuit of the Paul Allen building at the University of Washington.""", + 'license' : None, + 'size' : 284390}, + 'swiss_roll' : {'urls' : ['http://isomap.stanford.edu/'], + 'files' : [['swiss_roll_data.mat']], + 'details' : """Swiss roll data made available by Tenenbaum, de Silva and Langford to demonstrate isomap, available from http://isomap.stanford.edu/datasets.html.""", + 'citation' : 'A Global Geometric Framework for Nonlinear Dimensionality Reduction, J. B. Tenenbaum, V. de Silva and J. C. Langford, Science 290 (5500): 2319-2323, 22 December 2000', + 'license' : None, + 'size' : 800256}, + 'ripley_prnn_data' : {'urls' : ['http://www.stats.ox.ac.uk/pub/PRNN/'], + 'files' : [['Cushings.dat', 'README', 'crabs.dat', 'fglass.dat', 'fglass.grp', 'pima.te', 'pima.tr', 'pima.tr2', 'synth.te', 'synth.tr', 'viruses.dat', 'virus3.dat']], + 'details' : """Data sets from Brian Ripley's Pattern Recognition and Neural Networks""", + 'citation': """Pattern Recognition and Neural Networks by B.D. Ripley (1996) Cambridge University Press ISBN 0 521 46986 7""", + 'license' : None, + 'size' : 93565}, + 'isomap_face_data' : {'urls' : [neil_url + 'isomap_face_data/'], + 'files' : [['face_data.mat']], + 'details' : """Face data made available by Tenenbaum, de Silva and Langford to demonstrate isomap, available from http://isomap.stanford.edu/datasets.html.""", + 'citation' : 'A Global Geometric Framework for Nonlinear Dimensionality Reduction, J. B. Tenenbaum, V. de Silva and J. C. Langford, Science 290 (5500): 2319-2323, 22 December 2000', + 'license' : None, + 'size' : 24229368}, + 'xw_pen' : {'urls' : [neil_url + 'xw_pen/'], + 'files' : [['xw_pen_15.csv']], + 'details' : """Accelerometer pen data used for robust regression by Tipping and Lawrence.""", + 'citation' : 'Michael E. Tipping and Neil D. Lawrence. Variational inference for Student-t models: Robust Bayesian interpolation and generalised component analysis. Neurocomputing, 69:123--141, 2005', + 'license' : None, + 'size' : 3410} + } + +with open('data_resources.json', 'w') as file: + json.dump(data_resources, file) diff --git a/GPy/util/diag.py b/GPy/util/diag.py new file mode 100644 index 00000000..3d6b4dc9 --- /dev/null +++ b/GPy/util/diag.py @@ -0,0 +1,114 @@ +''' +.. module:: GPy.util.diag + +.. moduleauthor:: Max Zwiessele + +''' +__updated__ = '2013-12-03' + +import numpy as np + +def view(A, offset=0): + """ + Get a view on the diagonal elements of a 2D array. + + This is actually a view (!) on the diagonal of the array, so you can + in-place adjust the view. + + :param :class:`ndarray` A: 2 dimensional numpy array + :param int offset: view offset to give back (negative entries allowed) + :rtype: :class:`ndarray` view of diag(A) + + >>> import numpy as np + >>> X = np.arange(9).reshape(3,3) + >>> view(X) + array([0, 4, 8]) + >>> d = view(X) + >>> d += 2 + >>> view(X) + array([ 2, 6, 10]) + >>> view(X, offset=-1) + array([3, 7]) + >>> subtract(X, 3, offset=-1) + array([[ 2, 1, 2], + [ 0, 6, 5], + [ 6, 4, 10]]) + """ + from numpy.lib.stride_tricks import as_strided + assert A.ndim == 2, "only implemented for 2 dimensions" + assert A.shape[0] == A.shape[1], "attempting to get the view of non-square matrix?!" + if offset > 0: + return as_strided(A[0, offset:], shape=(A.shape[0] - offset, ), strides=((A.shape[0]+1)*A.itemsize, )) + elif offset < 0: + return as_strided(A[-offset:, 0], shape=(A.shape[0] + offset, ), strides=((A.shape[0]+1)*A.itemsize, )) + else: + return as_strided(A, shape=(A.shape[0], ), strides=((A.shape[0]+1)*A.itemsize, )) + +def _diag_ufunc(A,b,offset,func): + dA = view(A, offset); func(dA,b,dA) + return A + +def times(A, b, offset=0): + """ + Times the view of A with b in place (!). + Returns modified A + Broadcasting is allowed, thus b can be scalar. + + if offset is not zero, make sure b is of right shape! + + :param ndarray A: 2 dimensional array + :param ndarray-like b: either one dimensional or scalar + :param int offset: same as in view. + :rtype: view of A, which is adjusted inplace + """ + return _diag_ufunc(A, b, offset, np.multiply) +multiply = times + +def divide(A, b, offset=0): + """ + Divide the view of A by b in place (!). + Returns modified A + Broadcasting is allowed, thus b can be scalar. + + if offset is not zero, make sure b is of right shape! + + :param ndarray A: 2 dimensional array + :param ndarray-like b: either one dimensional or scalar + :param int offset: same as in view. + :rtype: view of A, which is adjusted inplace + """ + return _diag_ufunc(A, b, offset, np.divide) + +def add(A, b, offset=0): + """ + Add b to the view of A in place (!). + Returns modified A. + Broadcasting is allowed, thus b can be scalar. + + if offset is not zero, make sure b is of right shape! + + :param ndarray A: 2 dimensional array + :param ndarray-like b: either one dimensional or scalar + :param int offset: same as in view. + :rtype: view of A, which is adjusted inplace + """ + return _diag_ufunc(A, b, offset, np.add) + +def subtract(A, b, offset=0): + """ + Subtract b from the view of A in place (!). + Returns modified A. + Broadcasting is allowed, thus b can be scalar. + + if offset is not zero, make sure b is of right shape! + + :param ndarray A: 2 dimensional array + :param ndarray-like b: either one dimensional or scalar + :param int offset: same as in view. + :rtype: view of A, which is adjusted inplace + """ + return _diag_ufunc(A, b, offset, np.subtract) + +if __name__ == '__main__': + import doctest + doctest.testmod() \ No newline at end of file diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index 213cd047..4cafe044 100644 --- a/GPy/util/linalg.py +++ b/GPy/util/linalg.py @@ -12,19 +12,34 @@ import ctypes from ctypes import byref, c_char, c_int, c_double # TODO # import scipy.lib.lapack import scipy +import warnings +import os +from config import * if np.all(np.float64((scipy.__version__).split('.')[:2]) >= np.array([0, 12])): import scipy.linalg.lapack as lapack else: from scipy.linalg.lapack import flapack as lapack -try: - _blaslib = ctypes.cdll.LoadLibrary(np.core._dotblas.__file__) # @UndefinedVariable - _blas_available = True - assert hasattr('dsyrk_',_blaslib) - assert hasattr('dsyr_',_blaslib) -except: - _blas_available = False + +if config.getboolean('anaconda', 'installed') and config.getboolean('anaconda', 'MKL'): + try: + anaconda_path = str(config.get('anaconda', 'location')) + mkl_rt = ctypes.cdll.LoadLibrary(os.path.join(anaconda_path, 'DLLs', 'mkl_rt.dll')) + dsyrk = mkl_rt.dsyrk + dsyr = mkl_rt.dsyr + _blas_available = True + except: + _blas_available = False +else: + try: + _blaslib = ctypes.cdll.LoadLibrary(np.core._dotblas.__file__) # @UndefinedVariable + dsyrk = _blaslib.dsyrk_ + dsyr = _blaslib.dsyr_ + _blas_available = True + except AttributeError as e: + _blas_available = False + warnings.warn("warning: caught this exception:" + str(e)) def dtrtrs(A, B, lower=0, trans=0, unitdiag=0): """ @@ -61,6 +76,14 @@ def dpotri(A, lower=0): """ return lapack.dpotri(A, lower=lower) +def pddet(A): + """ + Determinant of a positive definite matrix, only symmetric matricies though + """ + L = jitchol(A) + logdetA = 2*sum(np.log(np.diag(L))) + return logdetA + def trace_dot(a, b): """ Efficiently compute the trace of the matrix product of a and b @@ -205,7 +228,7 @@ def multiple_pdinv(A): return np.dstack(invs), np.array(halflogdets) -def PCA(Y, input_dim): +def pca(Y, input_dim): """ Principal component analysis: maximum likelihood solution by SVD @@ -218,7 +241,7 @@ def PCA(Y, input_dim): """ if not np.allclose(Y.mean(axis=0), 0.0): - print "Y is not zero mean, centering it locally (GPy.util.linalg.PCA)" + print "Y is not zero mean, centering it locally (GPy.util.linalg.pca)" # Y -= Y.mean(axis=0) @@ -229,6 +252,124 @@ def PCA(Y, input_dim): W *= v; return X, W.T +def ppca(Y, Q, iterations=100): + """ + EM implementation for probabilistic pca. + + :param array-like Y: Observed Data + :param int Q: Dimensionality for reduced array + :param int iterations: number of iterations for EM + """ + from numpy.ma import dot as madot + N, D = Y.shape + # Initialise W randomly + W = np.random.randn(D, Q) * 1e-3 + Y = np.ma.masked_invalid(Y, copy=0) + mu = Y.mean(0) + Ycentered = Y - mu + try: + for _ in range(iterations): + exp_x = np.asarray_chkfinite(np.linalg.solve(W.T.dot(W), madot(W.T, Ycentered.T))).T + W = np.asarray_chkfinite(np.linalg.solve(exp_x.T.dot(exp_x), madot(exp_x.T, Ycentered))).T + except np.linalg.linalg.LinAlgError: + #"converged" + pass + return np.asarray_chkfinite(exp_x), np.asarray_chkfinite(W) + +def ppca_missing_data_at_random(Y, Q, iters=100): + """ + EM implementation of Probabilistic pca for when there is missing data. + + Taken from + + .. math: + \\mathbf{Y} = \mathbf{XW} + \\epsilon \\text{, where} + \\epsilon = \\mathcal{N}(0, \\sigma^2 \mathbf{I}) + + :returns: X, W, sigma^2 + """ + from numpy.ma import dot as madot + import diag + from GPy.util.subarray_and_sorting import common_subarrays + import time + debug = 1 + # Initialise W randomly + N, D = Y.shape + W = np.random.randn(Q, D) * 1e-3 + Y = np.ma.masked_invalid(Y, copy=1) + nu = 1. + #num_obs_i = 1./Y.count() + Ycentered = Y - Y.mean(0) + + X = np.zeros((N,Q)) + cs = common_subarrays(Y.mask) + cr = common_subarrays(Y.mask, 1) + Sigma = np.zeros((N, Q, Q)) + Sigma2 = np.zeros((N, Q, Q)) + mu = np.zeros(D) + if debug: + import matplotlib.pyplot as pylab + fig = pylab.figure("FIT MISSING DATA"); + ax = fig.gca() + ax.cla() + lines = pylab.plot(np.zeros((N,Q)).dot(W)) + W2 = np.zeros((Q,D)) + + for i in range(iters): +# Sigma = np.linalg.solve(diag.add(madot(W,W.T), nu), diag.times(np.eye(Q),nu)) +# exp_x = madot(madot(Ycentered, W.T),Sigma)/nu +# Ycentered = (Y - exp_x.dot(W).mean(0)) +# #import ipdb;ipdb.set_trace() +# #Ycentered = mu +# W = np.linalg.solve(madot(exp_x.T,exp_x) + Sigma, madot(exp_x.T, Ycentered)) +# nu = (((Ycentered - madot(exp_x, W))**2).sum(0) + madot(W.T,madot(Sigma,W)).sum(0)).sum()/N + for csi, (mask, index) in enumerate(cs.iteritems()): + mask = ~np.array(mask) + Sigma2[index, :, :] = nu * np.linalg.inv(diag.add(W2[:,mask].dot(W2[:,mask].T), nu)) + #X[index,:] = madot((Sigma[csi]/nu),madot(W,Ycentered[index].T))[:,0] + X2 = ((Sigma2/nu) * (madot(Ycentered,W2.T).base)[:,:,None]).sum(-1) + mu2 = (Y - X.dot(W)).mean(0) + for n in range(N): + Sigma[n] = nu * np.linalg.inv(diag.add(W[:,~Y.mask[n]].dot(W[:,~Y.mask[n]].T), nu)) + X[n, :] = (Sigma[n]/nu).dot(W[:,~Y.mask[n]].dot(Ycentered[n,~Y.mask[n]].T)) + for d in range(D): + mu[d] = (Y[~Y.mask[:,d], d] - X[~Y.mask[:,d]].dot(W[:, d])).mean() + Ycentered = (Y - mu) + nu3 = 0. + for cri, (mask, index) in enumerate(cr.iteritems()): + mask = ~np.array(mask) + W2[:,index] = np.linalg.solve(X[mask].T.dot(X[mask]) + Sigma[mask].sum(0), madot(X[mask].T, Ycentered[mask,index]))[:,None] + W2[:,index] = np.linalg.solve(X.T.dot(X) + Sigma.sum(0), madot(X.T, Ycentered[:,index])) + #nu += (((Ycentered[mask,index] - X[mask].dot(W[:,index]))**2).sum(0) + W[:,index].T.dot(Sigma[mask].sum(0).dot(W[:,index])).sum(0)).sum() + nu3 += (((Ycentered[index] - X.dot(W[:,index]))**2).sum(0) + W[:,index].T.dot(Sigma.sum(0).dot(W[:,index])).sum(0)).sum() + nu3 /= N + nu = 0. + nu2 = 0. + W = np.zeros((Q,D)) + for j in range(D): + W[:,j] = np.linalg.solve(X[~Y.mask[:,j]].T.dot(X[~Y.mask[:,j]]) + Sigma[~Y.mask[:,j]].sum(0), madot(X[~Y.mask[:,j]].T, Ycentered[~Y.mask[:,j],j])) + nu2f = np.tensordot(W[:,j].T, Sigma[~Y.mask[:,j],:,:], [0,1]).dot(W[:,j]) + nu2s = W[:,j].T.dot(Sigma[~Y.mask[:,j],:,:].sum(0).dot(W[:,j])) + nu2 += (((Ycentered[~Y.mask[:,j],j] - X[~Y.mask[:,j],:].dot(W[:,j]))**2) + nu2f).sum() + for i in range(N): + if not Y.mask[i,j]: + nu += ((Ycentered[i,j] - X[i,:].dot(W[:,j]))**2) + W[:,j].T.dot(Sigma[i,:,:].dot(W[:,j])) + nu /= N + nu2 /= N + nu4 = (((Ycentered - X.dot(W))**2).sum(0) + W.T.dot(Sigma.sum(0).dot(W)).sum(0)).sum()/N + import ipdb;ipdb.set_trace() + if debug: + #print Sigma[0] + print "nu:", nu, "sum(X):", X.sum() + pred_y = X.dot(W) + for x, l in zip(pred_y.T, lines): + l.set_ydata(x) + ax.autoscale_view() + ax.set_ylim(pred_y.min(), pred_y.max()) + fig.canvas.draw() + time.sleep(.3) + return np.asarray_chkfinite(X), np.asarray_chkfinite(W), nu + def tdot_numpy(mat, out=None): return np.dot(mat, mat.T, out) @@ -264,7 +405,7 @@ def tdot_blas(mat, out=None): BETA = c_double(0.0) C = out.ctypes.data_as(ctypes.c_void_p) LDC = c_int(np.max(out.strides) / 8) - _blaslib.dsyrk_(byref(UPLO), byref(TRANS), byref(N), byref(K), + dsyrk(byref(UPLO), byref(TRANS), byref(N), byref(K), byref(ALPHA), A, byref(LDA), byref(BETA), C, byref(LDC)) symmetrify(out, upper=True) @@ -294,7 +435,7 @@ def DSYR_blas(A, x, alpha=1.): A_ = A.ctypes.data_as(ctypes.c_void_p) x_ = x.ctypes.data_as(ctypes.c_void_p) INCX = c_int(1) - _blaslib.dsyr_(byref(UPLO), byref(N), byref(ALPHA), + dsyr(byref(UPLO), byref(N), byref(ALPHA), x_, byref(INCX), A_, byref(LDA)) symmetrify(A, upper=True) @@ -325,7 +466,7 @@ def symmetrify(A, upper=False): """ N, M = A.shape assert N == M - + c_contig_code = """ int iN; for (int i=1; i= c and D >= d: + index.append(rec[0]) + shape_records.append(rec[1]) + else: + cond1 = A <= a and B <= b and C >= a and D >= b + cond2 = A <= c and B <= d and C >= c and D >= d + cond3 = A <= a and D >= d and C >= a and B <= d + cond4 = A <= c and D >= b and C >= c and B <= b + cond5 = a <= C and b <= B and d >= D + cond6 = c <= A and b <= B and d >= D + cond7 = d <= B and a <= A and c >= C + cond8 = b <= D and a <= A and c >= C + if cond1 or cond2 or cond3 or cond4 or cond5 or cond6 or cond7 or cond8: + index.append(rec[0]) + shape_records.append(rec[1]) + return index,shape_records + + +def plot_bbox(sf,bbox,inside_only=True): + """ + Plot the geometry of a shapefile within a bbox + + :param sf: shapefile + :type sf: shapefile object + :param bbox: bounding box + :type bbox: list of floats [x_min,y_min,x_max,y_max] + :inside_only: True if the objects returned are those that lie within the bbox and False if the objects returned are any that intersect the bbox + :type inside_only: Boolean + """ + index,shape_records = bbox_match(sf,bbox,inside_only) + A,B,C,D = bbox + plot(shape_records,xlims=[bbox[0],bbox[2]],ylims=[bbox[1],bbox[3]]) + +def plot_string_match(sf,regex,field): + """ + Plot the geometry of a shapefile whose fields match a regular expression given + + :param sf: shapefile + :type sf: shapefile object + :regex: regular expression to match + :type regex: string + :field: field number to be matched with the regex + :type field: integer + """ + index,shape_records = string_match(sf,regex,field) + plot(shape_records) + + +def new_shape_string(sf,name,regex,field=2,type=shapefile.POINT): + + newshp = shapefile.Writer(shapeType = sf.shapeType) + newshp.autoBalance = 1 + + index,shape_records = string_match(sf,regex,field) + + _fi = [sf.fields[j] for j in index] + for f in _fi: + newshp.field(name=f[0],fieldType=f[1],size=f[2],decimal=f[3]) + + _shre = shape_records + for sr in _shre: + _points = [] + _parts = [] + for point in sr.shape.points: + _points.append(point) + _parts.append(_points) + + newshp.line(parts=_parts) + newshp.records.append(sr.record) + print len(sr.record) + + newshp.save(name) + print index diff --git a/GPy/util/misc.py b/GPy/util/misc.py index 8cf5806f..1cb4c182 100644 --- a/GPy/util/misc.py +++ b/GPy/util/misc.py @@ -32,18 +32,6 @@ def chain_3(d3f_dg3, dg_dx, d2f_dg2, d2g_dx2, df_dg, d3g_dx3): """ return d3f_dg3*(dg_dx**3) + 3*d2f_dg2*dg_dx*d2g_dx2 + df_dg*d3g_dx3 -### make a parameter to its corresponding array: -def param_to_array(*param): - """ - Convert an arbitrary number of parameters to :class:ndarray class objects. This is for - converting parameter objects to numpy arrays, when using scipy.weave.inline routine. - In scipy.weave.blitz there is no automatic array detection (even when the array inherits - from :class:ndarray)""" - assert len(param) > 0, "At least one parameter needed" - if len(param) == 1: - return param[0].view(np.ndarray) - return map(lambda x: x.view(np.ndarray), param) - def opt_wrapper(m, **kwargs): """ This function just wraps the optimization procedure of a GPy @@ -171,7 +159,6 @@ def fast_array_equal(A, B): elif ((A == None) and (B != None)) or ((A != None) and (B == None)): return False elif A.shape == B.shape: - A, B = param_to_array(A, B) if A.ndim == 2: N, D = [int(i) for i in A.shape] value = weave.inline(code2, support_code=support_code, @@ -187,6 +174,7 @@ def fast_array_equal(A, B): return value + if __name__ == '__main__': import pylab as plt X = np.linspace(1,10, 100)[:, None] diff --git a/GPy/util/mocap.py b/GPy/util/mocap.py index 78f00955..58662cf9 100644 --- a/GPy/util/mocap.py +++ b/GPy/util/mocap.py @@ -67,14 +67,14 @@ class tree: for i in range(len(self.vertices)): if self.vertices[i].id == id: return i - raise Error, 'Reverse look up of id failed.' + raise ValueError('Reverse look up of id failed.') def get_index_by_name(self, name): """Give the index associated with a given vertex name.""" for i in range(len(self.vertices)): if self.vertices[i].name == name: return i - raise Error, 'Reverse look up of name failed.' + raise ValueError('Reverse look up of name failed.') def order_vertices(self): """Order vertices in the graph such that parents always have a lower index than children.""" @@ -433,6 +433,8 @@ class acclaim_skeleton(skeleton): lin = self.read_line(fid) while lin != ':DEGREES': lin = self.read_line(fid) + if lin == '': + raise ValueError('Could not find :DEGREES in ' + fid.name) counter = 0 lin = self.read_line(fid) @@ -443,9 +445,9 @@ class acclaim_skeleton(skeleton): if frame_no: counter += 1 if counter != frame_no: - raise Error, 'Unexpected frame number.' + raise ValueError('Unexpected frame number.') else: - raise Error, 'Single bone name ...' + raise ValueError('Single bone name ...') else: ind = self.get_index_by_name(parts[0]) bones[ind].append(np.array([float(channel) for channel in parts[1:]])) @@ -573,7 +575,7 @@ class acclaim_skeleton(skeleton): return lin = self.read_line(fid) else: - raise Error, 'Unrecognised file format' + raise ValueError('Unrecognised file format') self.finalize() def read_units(self, fid): diff --git a/GPy/util/plot_latent.py b/GPy/util/plot_latent.py index cecb811c..207a7974 100644 --- a/GPy/util/plot_latent.py +++ b/GPy/util/plot_latent.py @@ -20,8 +20,8 @@ def most_significant_input_dimensions(model, which_indices): input_1, input_2 = which_indices return input_1, input_2 -def plot_latent(model, labels=None, which_indices=None, - resolution=50, ax=None, marker='o', s=40, +def plot_latent(model, labels=None, which_indices=None, + resolution=50, ax=None, marker='o', s=40, fignum=None, plot_inducing=False, legend=True, aspect='auto', updates=False): """ @@ -38,11 +38,9 @@ def plot_latent(model, labels=None, which_indices=None, input_1, input_2 = most_significant_input_dimensions(model, which_indices) - X = np.asarray(model.X) - # first, plot the output variance as a function of the latent space - Xtest, xx, yy, xmin, xmax = util.plot.x_frame2D(X[:, [input_1, input_2]], resolution=resolution) - Xtest_full = np.zeros((Xtest.shape[0], X.shape[1])) + Xtest, xx, yy, xmin, xmax = util.plot.x_frame2D(model.X[:, [input_1, input_2]], resolution=resolution) + Xtest_full = np.zeros((Xtest.shape[0], model.X.shape[1])) def plot_function(x): Xtest_full[:, [input_1, input_2]] = x @@ -50,10 +48,10 @@ def plot_latent(model, labels=None, which_indices=None, var = var[:, :1] return np.log(var) view = ImshowController(ax, plot_function, - tuple(X[:, [input_1, input_2]].min(0)) + tuple(X[:, [input_1, input_2]].max(0)), + tuple(model.X[:, [input_1, input_2]].min(0)) + tuple(model.X[:, [input_1, input_2]].max(0)), resolution, aspect=aspect, interpolation='bilinear', cmap=pb.cm.binary) - + # ax.imshow(var.reshape(resolution, resolution).T, # extent=[xmin[0], xmax[0], xmin[1], xmax[1]], cmap=pb.cm.binary, interpolation='bilinear', origin='lower') @@ -76,11 +74,11 @@ def plot_latent(model, labels=None, which_indices=None, index = np.nonzero(labels == ul)[0] if model.input_dim == 1: - x = X[index, input_1] + x = model.X[index, input_1] y = np.zeros(index.size) else: - x = X[index, input_1] - y = X[index, input_2] + x = model.X[index, input_1] + y = model.X[index, input_2] ax.scatter(x, y, marker=m, s=s, color=util.plot.Tango.nextMedium(), label=this_label) ax.set_xlabel('latent dimension %i' % input_1) @@ -102,8 +100,8 @@ def plot_latent(model, labels=None, which_indices=None, raw_input('Enter to continue') return ax -def plot_magnification(model, labels=None, which_indices=None, - resolution=60, ax=None, marker='o', s=40, +def plot_magnification(model, labels=None, which_indices=None, + resolution=60, ax=None, marker='o', s=40, fignum=None, plot_inducing=False, legend=True, aspect='auto', updates=False): """ @@ -119,17 +117,16 @@ def plot_magnification(model, labels=None, which_indices=None, labels = np.ones(model.num_data) input_1, input_2 = most_significant_input_dimensions(model, which_indices) - X = np.asarray(model.X) - + # first, plot the output variance as a function of the latent space - Xtest, xx, yy, xmin, xmax = util.plot.x_frame2D(X[:, [input_1, input_2]], resolution=resolution) - Xtest_full = np.zeros((Xtest.shape[0], X.shape[1])) + Xtest, xx, yy, xmin, xmax = util.plot.x_frame2D(model.X[:, [input_1, input_2]], resolution=resolution) + Xtest_full = np.zeros((Xtest.shape[0], model.X.shape[1])) def plot_function(x): Xtest_full[:, [input_1, input_2]] = x mf=model.magnification(Xtest_full) return mf view = ImshowController(ax, plot_function, - tuple(X.min(0)[:, [input_1, input_2]]) + tuple(X.max(0)[:, [input_1, input_2]]), + tuple(model.X.min(0)[:, [input_1, input_2]]) + tuple(model.X.max(0)[:, [input_1, input_2]]), resolution, aspect=aspect, interpolation='bilinear', cmap=pb.cm.gray) @@ -152,11 +149,11 @@ def plot_magnification(model, labels=None, which_indices=None, index = np.nonzero(labels == ul)[0] if model.input_dim == 1: - x = X[index, input_1] + x = model.X[index, input_1] y = np.zeros(index.size) else: - x = X[index, input_1] - y = X[index, input_2] + x = model.X[index, input_1] + y = model.X[index, input_2] ax.scatter(x, y, marker=m, s=s, color=util.plot.Tango.nextMedium(), label=this_label) ax.set_xlabel('latent dimension %i' % input_1) diff --git a/GPy/util/subarray_and_sorting.py b/GPy/util/subarray_and_sorting.py new file mode 100644 index 00000000..49385771 --- /dev/null +++ b/GPy/util/subarray_and_sorting.py @@ -0,0 +1,56 @@ +''' +.. module:: GPy.util.subarray_and_sorting + +.. moduleauthor:: Max Zwiessele + +''' +__updated__ = '2013-12-02' + +import numpy as np + +def common_subarrays(X, axis=0): + """ + Find common subarrays of 2 dimensional X, where axis is the axis to apply the search over. + Common subarrays are returned as a dictionary of pairs, where + the subarray is a tuple representing the subarray and the index is the index + for the subarray in X, where index is the index to the remaining axis. + + :param :class:`np.ndarray` X: 2d array to check for common subarrays in + :param int axis: axis to apply subarray detection over. + When the index is 0, compare rows, columns, otherwise. + + Examples: + ========= + + In a 2d array: + >>> import numpy as np + >>> X = np.zeros((3,6), dtype=bool) + >>> X[[1,1,1],[0,4,5]] = 1; X[1:,[2,3]] = 1 + >>> X + array([[False, False, False, False, False, False], + [ True, False, True, True, True, True], + [False, False, True, True, False, False]], dtype=bool) + >>> d = common_subarrays(X,axis=1) + >>> len(d) + 3 + >>> X[:, d[tuple(X[:,0])]] + array([[False, False, False], + [ True, True, True], + [False, False, False]], dtype=bool) + >>> d[tuple(X[:,4])] == d[tuple(X[:,0])] == [0, 4, 5] + True + >>> d[tuple(X[:,1])] + [1] + """ + from collections import defaultdict + from itertools import count + from operator import iadd + assert X.ndim == 2 and axis in (0,1), "Only implemented for 2D arrays" + subarrays = defaultdict(list) + cnt = count() + np.apply_along_axis(lambda x: iadd(subarrays[tuple(x)], [cnt.next()]), 1-axis, X) + return subarrays + +if __name__ == '__main__': + import doctest + doctest.testmod() \ No newline at end of file diff --git a/GPy/util/symbolic.py b/GPy/util/symbolic.py index a49a507d..49c8c33a 100644 --- a/GPy/util/symbolic.py +++ b/GPy/util/symbolic.py @@ -1,4 +1,4 @@ -from sympy import Function, S, oo, I, cos, sin, asin, log, erf,pi,exp, sqrt +from sympy import Function, S, oo, I, cos, sin, asin, log, erf, pi, exp, sqrt, sign class ln_diff_erf(Function): @@ -10,7 +10,7 @@ class ln_diff_erf(Function): return -2*exp(-x1**2)/(sqrt(pi)*(erf(x0)-erf(x1))) elif argindex == 1: x0, x1 = self.args - return 2*exp(-x0**2)/(sqrt(pi)*(erf(x0)-erf(x1))) + return 2.*exp(-x0**2)/(sqrt(pi)*(erf(x0)-erf(x1))) else: raise ArgumentIndexError(self, argindex) @@ -19,25 +19,209 @@ class ln_diff_erf(Function): if x0.is_Number and x1.is_Number: return log(erf(x0)-erf(x1)) -class sim_h(Function): +class dh_dd_i(Function): nargs = 5 - - def fdiff(self, argindex=5): - t, tprime, d_i, d_j, l = self.args - if argindex == 5: - return -2*sin(t) - elif argindex == 4: - return 2*exp(d_i) - elif argindex == 3: - return 4*exp(d_j) - elif argindex == 2: - return 3*exp(l) - elif argindex == 1: - return 2*d_j - @classmethod def eval(cls, t, tprime, d_i, d_j, l): - return exp((d_j/2*l)**2)/(d_i+d_j)*(exp(-d_j*(tprime - t))*(erf((tprime-t)/l - d_j/2*l) + erf(t/l + d_j/2*l)) - exp(-(d_j*tprime + d_i))*(erf(tprime/l - d_j/2*l) + erf(d_j/2*l))) + if (t.is_Number + and tprime.is_Number + and d_i.is_Number + and d_j.is_Number + and l.is_Number): + + diff_t = (t-tprime) + l2 = l*l + h = h(t, tprime, d_i, d_j, l) + half_l_di = 0.5*l*d_i + arg_1 = half_l_di + tprime/l + arg_2 = half_l_di - (t-tprime)/l + ln_part_1 = ln_diff_erf(arg_1, arg_2) + arg_1 = half_l_di + arg_2 = half_l_di - t/l + sign_val = sign(t/l) + ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t/l) + + base = ((0.5*d_i*l2*(d_i+d_j)-1)*h + + (-diff_t*sign_val*exp(half_l_di*half_l_di + -d_i*diff_t + +ln_part_1) + +t*sign_val*exp(half_l_di*half_l_di + -d_i*t-d_j*tprime + +ln_part_2)) + + l/sqrt(pi)*(-exp(-diff_t*diff_t/l2) + +exp(-tprime*tprime/l2-d_i*t) + +exp(-t*t/l2-d_j*tprime) + -exp(-(d_i*t + d_j*tprime)))) + return base/(d_i+d_j) + +class dh_dd_j(Function): + nargs = 5 + @classmethod + def eval(cls, t, tprime, d_i, d_j, l): + if (t.is_Number + and tprime.is_Number + and d_i.is_Number + and d_j.is_Number + and l.is_Number): + diff_t = (t-tprime) + l2 = l*l + half_l_di = 0.5*l*d_i + h = h(t, tprime, d_i, d_j, l) + arg_1 = half_l_di + tprime/l + arg_2 = half_l_di - (t-tprime)/l + ln_part_1 = ln_diff_erf(arg_1, arg_2) + arg_1 = half_l_di + arg_2 = half_l_di - t/l + sign_val = sign(t/l) + ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t/l) + sign_val = sign(t/l) + base = tprime*sign_val*exp(half_l_di*half_l_di-(d_i*t+d_j*tprime)+ln_part_2)-h + return base/(d_i+d_j) + +class dh_dl(Function): + nargs = 5 + @classmethod + def eval(cls, t, tprime, d_i, d_j, l): + if (t.is_Number + and tprime.is_Number + and d_i.is_Number + and d_j.is_Number + and l.is_Number): + + diff_t = (t-tprime) + l2 = l*l + h = h(t, tprime, d_i, d_j, l) + return 0.5*d_i*d_i*l*h + 2./(sqrt(pi)*(d_i+d_j))*((-diff_t/l2-d_i/2.)*exp(-diff_t*diff_t/l2)+(-tprime/l2+d_i/2.)*exp(-tprime*tprime/l2-d_i*t)-(-t/l2-d_i/2.)*exp(-t*t/l2-d_j*tprime)-d_i/2.*exp(-(d_i*t+d_j*tprime))) + +class dh_dt(Function): + nargs = 5 + @classmethod + def eval(cls, t, tprime, d_i, d_j, l): + if (t.is_Number + and tprime.is_Number + and d_i.is_Number + and d_j.is_Number + and l.is_Number): + if (t is S.NaN + or tprime is S.NaN + or d_i is S.NaN + or d_j is S.NaN + or l is S.NaN): + return S.NaN + else: + half_l_di = 0.5*l*d_i + arg_1 = half_l_di + tprime/l + arg_2 = half_l_di - (t-tprime)/l + ln_part_1 = ln_diff_erf(arg_1, arg_2) + arg_1 = half_l_di + arg_2 = half_l_di - t/l + sign_val = sign(t/l) + ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t/l) + + + return (sign_val*exp(half_l_di*half_l_di + - d_i*(t-tprime) + + ln_part_1 + - log(d_i + d_j)) + - sign_val*exp(half_l_di*half_l_di + - d_i*t - d_j*tprime + + ln_part_2 + - log(d_i + d_j))).diff(t) + +class dh_dtprime(Function): + nargs = 5 + @classmethod + def eval(cls, t, tprime, d_i, d_j, l): + if (t.is_Number + and tprime.is_Number + and d_i.is_Number + and d_j.is_Number + and l.is_Number): + if (t is S.NaN + or tprime is S.NaN + or d_i is S.NaN + or d_j is S.NaN + or l is S.NaN): + return S.NaN + else: + half_l_di = 0.5*l*d_i + arg_1 = half_l_di + tprime/l + arg_2 = half_l_di - (t-tprime)/l + ln_part_1 = ln_diff_erf(arg_1, arg_2) + arg_1 = half_l_di + arg_2 = half_l_di - t/l + sign_val = sign(t/l) + ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t/l) + + + return (sign_val*exp(half_l_di*half_l_di + - d_i*(t-tprime) + + ln_part_1 + - log(d_i + d_j)) + - sign_val*exp(half_l_di*half_l_di + - d_i*t - d_j*tprime + + ln_part_2 + - log(d_i + d_j))).diff(tprime) + + +class h(Function): + nargs = 5 + def fdiff(self, argindex=5): + t, tprime, d_i, d_j, l = self.args + if argindex == 1: + return dh_dt(t, tprime, d_i, d_j, l) + elif argindex == 2: + return dh_dtprime(t, tprime, d_i, d_j, l) + elif argindex == 3: + return dh_dd_i(t, tprime, d_i, d_j, l) + elif argindex == 4: + return dh_dd_j(t, tprime, d_i, d_j, l) + elif argindex == 5: + return dh_dl(t, tprime, d_i, d_j, l) + + + @classmethod + def eval(cls, t, tprime, d_i, d_j, l): + # putting in the is_Number stuff forces it to look for a fdiff method for derivative. If it's left out, then when asking for self.diff, it just does the diff on the eval symbolic terms directly. We want to avoid that because we are looking to ensure everything is numerically stable. Maybe it's because of the if statement that this happens? + if (t.is_Number + and tprime.is_Number + and d_i.is_Number + and d_j.is_Number + and l.is_Number): + if (t is S.NaN + or tprime is S.NaN + or d_i is S.NaN + or d_j is S.NaN + or l is S.NaN): + return S.NaN + else: + half_l_di = 0.5*l*d_i + arg_1 = half_l_di + tprime/l + arg_2 = half_l_di - (t-tprime)/l + ln_part_1 = ln_diff_erf(arg_1, arg_2) + arg_1 = half_l_di + arg_2 = half_l_di - t/l + sign_val = sign(t/l) + ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t/l) + + + return (sign_val*exp(half_l_di*half_l_di + - d_i*(t-tprime) + + ln_part_1 + - log(d_i + d_j)) + - sign_val*exp(half_l_di*half_l_di + - d_i*t - d_j*tprime + + ln_part_2 + - log(d_i + d_j))) + + + # return (exp((d_j/2.*l)**2)/(d_i+d_j) + # *(exp(-d_j*(tprime - t)) + # *(erf((tprime-t)/l - d_j/2.*l) + # + erf(t/l + d_j/2.*l)) + # - exp(-(d_j*tprime + d_i)) + # *(erf(tprime/l - d_j/2.*l) + # + erf(d_j/2.*l)))) class erfc(Function): nargs = 1 @@ -53,52 +237,3 @@ class erfcx(Function): def eval(cls, arg): return erfc(arg)*exp(arg*arg) -class sinc_grad(Function): - nargs = 1 - - def fdiff(self, argindex=1): - if argindex==1: - # Strictly speaking this should be computed separately, as it won't work when x=0. See http://calculus.subwiki.org/wiki/Sinc_function - return ((2-x*x)*sin(self.args[0]) - 2*x*cos(x))/(x*x*x) - else: - raise ArgumentIndexError(self, argindex) - - - @classmethod - def eval(cls, x): - if x.is_Number: - if x is S.NaN: - return S.NaN - elif x is S.Zero: - return S.Zero - else: - return (x*cos(x) - sin(x))/(x*x) - -class sinc(Function): - - nargs = 1 - - def fdiff(self, argindex=1): - if argindex==1: - return sinc_grad(self.args[0]) - else: - raise ArgumentIndexError(self, argindex) - - - @classmethod - def eval(cls, arg): - if arg.is_Number: - if arg is S.NaN: - return S.NaN - elif arg is S.Zero: - return S.One - else: - return sin(arg)/arg - - if arg.func is asin: - x = arg.args[0] - return x / arg - - def _eval_is_real(self): - return self.args[0].is_real - diff --git a/GPy/util/univariate_Gaussian.py b/GPy/util/univariate_Gaussian.py index 5a5880d5..702ab25c 100644 --- a/GPy/util/univariate_Gaussian.py +++ b/GPy/util/univariate_Gaussian.py @@ -13,24 +13,32 @@ def std_norm_cdf(x): Cumulative standard Gaussian distribution Based on Abramowitz, M. and Stegun, I. (1970) """ + #Generalize for many x + x = np.asarray(x).copy() + cdf_x = np.zeros_like(x) + N = x.size support_code = "#include " code = """ - double sign = 1.0; - if (x < 0.0){ - sign = -1.0; - x = -x; + double sign, t, erf; + for (int i=0; i