{"action":"created","issue":{"url":"https://api.github.com/repos/markovmodel/PyEMMA/issues/759","repository_url":"https://api.github.com/repos/markovmodel/PyEMMA","labels_url":"https://api.github.com/repos/markovmodel/PyEMMA/issues/759/labels{/name}","comments_url":"https://api.github.com/repos/markovmodel/PyEMMA/issues/759/comments","events_url":"https://api.github.com/repos/markovmodel/PyEMMA/issues/759/events","html_url":"https://github.com/markovmodel/PyEMMA/issues/759","id":146090393,"number":759,"title":"Some difficulties with minRMSD clustering: chunk size shape mismatch","user":{"login":"jchodera","id":3656088,"avatar_url":"https://avatars.githubusercontent.com/u/3656088?v=3","gravatar_id":"","url":"https://api.github.com/users/jchodera","html_url":"https://github.com/jchodera","followers_url":"https://api.github.com/users/jchodera/followers","following_url":"https://api.github.com/users/jchodera/following{/other_user}","gists_url":"https://api.github.com/users/jchodera/gists{/gist_id}","starred_url":"https://api.github.com/users/jchodera/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/jchodera/subscriptions","organizations_url":"https://api.github.com/users/jchodera/orgs","repos_url":"https://api.github.com/users/jchodera/repos","events_url":"https://api.github.com/users/jchodera/events{/privacy}","received_events_url":"https://api.github.com/users/jchodera/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":15,"created_at":"2016-04-05T19:42:34Z","updated_at":"2016-04-11T14:34:38Z","closed_at":null,"body":"I'm not sure if this is a silly mistake on my part or an issue with the chunking, but I've run into an odd error in minRMSD clustering of a large dataset:\r\n```\r\n[chodera@mskcc-ln1 ~]$ cat cluster-CK2.o7065714\r\nWarning: no access to tty (Bad file descriptor).\r\nThus no job control in this shell.\r\nMon Apr  4 17:10:36 EDT 2016\r\nloading reference topology...\r\nInitializing featurizer...\r\nThere are 988951 frames total in 557 trajectories.\r\nClustering...\r\nTraceback (most recent call last):\r\n  File \"cluster.py\", line 51, in <module>\r\n    dtrajs = clustering.dtrajs\r\n  File \"/cbio/jclab/home/chodera/miniconda/lib/python2.7/site-packages/pyEMMA-2.0.4-py2.7-linux-x86_64.egg/pyemma/coordinates/clustering/interface.py\", line 73, in dtrajs\r\n    self._dtrajs = self.assign(stride=1)\r\n  File \"/cbio/jclab/home/chodera/miniconda/lib/python2.7/site-packages/pyEMMA-2.0.4-py2.7-linux-x86_64.egg/pyemma/coordinates/clustering/interface.py\", line 180, in assign\r\n    mapped = self.get_output(stride=stride)\r\n  File \"/cbio/jclab/home/chodera/miniconda/lib/python2.7/site-packages/pyEMMA-2.0.4-py2.7-linux-x86_64.egg/pyemma/coordinates/transform/transformer.py\", line 784, in get_output\r\n    trajs[itraj][t:t + L, :] = chunk[:, dimensions]\r\nValueError: could not broadcast input array from shape (100,1) into shape (60,1)\r\n```\r\nHere is the simple clustering script:\r\n```python\r\n#!/usr/bin/env python\r\n\r\nimport pyemma\r\nimport numpy as np\r\nimport mdtraj\r\nimport os\r\n\r\n# Source directory\r\nsource_directory = '/cbio/jclab/projects/fah/fah-data/munged-with-time/no-solvent/11406' # CK2\r\n\r\n################################################################################\r\n# Load reference topology\r\n################################################################################\r\n\r\nprint ('loading reference topology...')\r\nreference_pdb_filename = 'protein.pdb'\r\nreference_trajectory = os.path.join(source_directory, 'run0-clone0.h5')\r\ntraj = mdtraj.load(reference_trajectory)\r\ntraj[0].save_pdb(reference_pdb_filename)\r\n\r\n################################################################################\r\n# Initialize featurizer\r\n################################################################################\r\n\r\nprint('Initializing featurizer...')\r\nimport pyemma.coordinates\r\nfeaturizer = pyemma.coordinates.featurizer(reference_pdb_filename)\r\nfeaturizer.add_all()\r\n\r\n################################################################################\r\n# Define coordinates source\r\n################################################################################\r\n\r\nimport pyemma.coordinates\r\nfrom glob import glob\r\ntrajectory_filenames = glob(os.path.join(source_directory, 'run*-clone*.h5'))\r\ncoordinates_source = pyemma.coordinates.source(trajectory_filenames, features=featurizer)\r\nprint(\"There are %d frames total in %d trajectories.\" % (coordinates_source.n_frames_total(), coordinates_source.number_of_trajectories()))\r\n\r\n################################################################################\r\n# Cluster\r\n################################################################################\r\n\r\nprint('Clustering...')\r\ngenerator_ratio = 100\r\nnframes = coordinates_source.n_frames_total()\r\nnstates = int(nframes / generator_ratio)\r\nstride = 4\r\nmetric = 'minRMSD'\r\nclustering = pyemma.coordinates.cluster_uniform_time(data=coordinates_source, k=nstates, stride=stride, metric=metric)\r\ndtrajs = clustering.dtrajs\r\n\r\n# Save discrete trajectories.\r\ndtrajs_dir = 'dtrajs'\r\nclustering.save_dtrajs(output_dir=dtrajs_dir, output_format='npy', extension='.npy')\r\n```\r\nUnfortunately, the dataset is enormous, but I'm happy to get someone an account on our cluster if that would make things easier.\r\n\r\nI'm trying to find a minimum example that demonstrates the same issue.\r\n\r\nI'm using `pyemma` 2.0.4"},"comment":{"url":"https://api.github.com/repos/markovmodel/PyEMMA/issues/comments/208376903","html_url":"https://github.com/markovmodel/PyEMMA/issues/759#issuecomment-208376903","issue_url":"https://api.github.com/repos/markovmodel/PyEMMA/issues/759","id":208376903,"user":{"login":"clonker","id":1685266,"avatar_url":"https://avatars.githubusercontent.com/u/1685266?v=3","gravatar_id":"","url":"https://api.github.com/users/clonker","html_url":"https://github.com/clonker","followers_url":"https://api.github.com/users/clonker/followers","following_url":"https://api.github.com/users/clonker/following{/other_user}","gists_url":"https://api.github.com/users/clonker/gists{/gist_id}","starred_url":"https://api.github.com/users/clonker/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/clonker/subscriptions","organizations_url":"https://api.github.com/users/clonker/orgs","repos_url":"https://api.github.com/users/clonker/repos","events_url":"https://api.github.com/users/clonker/events{/privacy}","received_events_url":"https://api.github.com/users/clonker/received_events","type":"User","site_admin":false},"created_at":"2016-04-11T14:34:38Z","updated_at":"2016-04-11T14:34:38Z","body":"`_itraj` indeed can become `== self.number_of_trajectories()`, however the subsequent access to `self._data_source_filenames[self._itraj]` happens only if:\r\n\r\n`_itraj < self.number_of_trajectories() (== self._data_source.ntraj)`.\r\n\r\nThe mistake is rather that in the case of `_itraj >= self.number_of_trajectories()` no `StopIteration()` is raised, which could be fixed in a general way like @marscher suggested above."}}