{"action":"created","comment":{"url":"https://api.github.com/repos/scipy/scipy/pulls/comments/28204431","id":28204431,"diff_hunk":"@@ -658,58 +660,101 @@ def lsim(system, U, T, X0=None, interp=1):\n         sys = system\n     else:\n         sys = lti(*system)\n-    U = atleast_1d(U)\n     T = atleast_1d(T)\n-    if len(U.shape) == 1:\n-        U = U.reshape((U.shape[0], 1))\n-    sU = U.shape\n     if len(T.shape) != 1:\n         raise ValueError(\"T must be a rank-1 array.\")\n-    if sU[0] != len(T):\n-        raise ValueError(\"U must have the same number of rows \"\n-                         \"as elements in T.\")\n-    if sU[1] != sys.inputs:\n-        raise ValueError(\"System does not define that many inputs.\")\n \n+    A, B, C, D = (np.asarray(M) for M in (sys.A, sys.B, sys.C, sys.D))\n+    n_states = A.shape[0]\n+    n_inputs = B.shape[1]\n+\n+    n_steps = len(T)\n     if X0 is None:\n-        X0 = zeros(sys.B.shape[0], sys.A.dtype)\n+        X0 = zeros(n_states, sys.A.dtype)\n+    xout = zeros((n_steps, n_states), sys.A.dtype)\n+\n+    if T[0] == 0:\n+        xout[0] = X0\n+    elif T[0] > 0:\n+        # step forward to initial time, with zero input\n+        xout[0] = dot(X0, linalg.expm(transpose(A) * T[0]))\n+    else:\n+        raise ValueError(\"Initial time must be nonnegative\")\n+\n+    no_input = (U is None\n+                or (isinstance(U, (int, float)) and U == 0.)\n+                or np.allclose(U, 0))\n+\n+    if n_steps == 1:\n+        yout = squeeze(dot(xout, transpose(C)))\n+        if not no_input:\n+            yout += squeeze(dot(U, transpose(D)))\n+        return T, squeeze(yout), squeeze(xout)\n \n-    xout = zeros((len(T), sys.B.shape[0]), sys.A.dtype)\n-    xout[0] = X0\n-    A = sys.A\n-    AT, BT = transpose(sys.A), transpose(sys.B)\n     dt = T[1] - T[0]\n-    lam, v = linalg.eig(A)\n-    vt = transpose(v)\n-    vti = linalg.inv(vt)\n-    GT = dot(dot(vti, diag(numpy.exp(dt * lam))), vt)\n-    GT = _cast_to_array_dtype(GT, xout)\n-\n-    ATm1 = linalg.inv(AT)\n-    ATm2 = dot(ATm1, ATm1)\n-    I = eye(A.shape[0], dtype=A.dtype)\n-    GTmI = GT - I\n-    F1T = dot(dot(BT, GTmI), ATm1)\n-    if interp:\n-        F2T = dot(BT, dot(GTmI, ATm2) / dt - ATm1)\n-\n-    for k in xrange(1, len(T)):\n-        dt1 = T[k] - T[k - 1]\n-        if dt1 != dt:\n-            dt = dt1\n-            GT = dot(dot(vti, diag(numpy.exp(dt * lam))), vt)\n-            GT = _cast_to_array_dtype(GT, xout)\n-            GTmI = GT - I\n-            F1T = dot(dot(BT, GTmI), ATm1)\n-            if interp:\n-                F2T = dot(BT, dot(GTmI, ATm2) / dt - ATm1)\n-\n-        xout[k] = dot(xout[k - 1], GT) + dot(U[k - 1], F1T)\n-        if interp:\n-            xout[k] = xout[k] + dot((U[k] - U[k - 1]), F2T)\n-\n-    yout = (squeeze(dot(U, transpose(sys.D))) +\n-            squeeze(dot(xout, transpose(sys.C))))\n+    if not np.allclose(T[1:] - T[:-1], dt):\n+        raise ValueError(\"Time values must be equally spaced.\")\n+\n+    if no_input:\n+        # Zero input: just use matrix exponential\n+        # take transpose because state is a row vector\n+        expAT_dt = linalg.expm(transpose(A) * dt)\n+        for i in xrange(1, n_steps):\n+            xout[i] = dot(xout[i-1], expAT_dt)\n+        yout = squeeze(dot(xout, transpose(C)))\n+        return T, squeeze(yout), squeeze(xout)\n+\n+    # Nonzero input\n+    U = atleast_1d(U)\n+    if len(U.shape) == 1:\n+        U = U.reshape((U.shape[0], 1))\n+    sU = U.shape\n+    if sU[0] != n_steps:\n+        raise ValueError(\"U must have the same number of rows \"\n+                         \"as elements in T.\")\n+    if sU[1] != n_inputs:\n+        raise ValueError(\"System does not define that many inputs.\")\n+    if not interp:\n+        # Zero-order hold\n+        # Algorithm: to integrate from time 0 to time dt, we solve\n+        #   xdot = A x + B u,  x(0) = x0\n+        #   udot = 0,          u(0) = u0.\n+        #\n+        # Solution is\n+        #   [ x(dt) ]       [ A*dt   B*dt ] [ x0 ]\n+        #   [ u(dt) ] = exp [  0     0    ] [ u0 ]\n+        M = np.bmat([[A * dt, B * dt],\n+                     [np.zeros((n_inputs, n_states + n_inputs))]])\n+        # transpose everything because the state and input are row vectors\n+        expMT = linalg.expm(transpose(M))\n+        Ad = expMT[:n_states, :n_states]\n+        Bd = expMT[n_states:, :n_states]\n+        for i in xrange(1, n_steps):\n+            xout[i] = dot(xout[i-1], Ad) + dot(U[i-1], Bd)\n+    else:\n+        # Linear interpolation between steps\n+        # Algorithm: to integrate from time 0 to time dt, with linear\n+        # interpolation between inputs u(0) = u0 and u(dt) = u1, we solve\n+        #   xdot = A x + B u,        x(0) = x0\n+        #   udot = (u1 - u0) / dt,   u(0) = u0.\n+        #\n+        # Solution is\n+        #   [ x(dt) ]       [ A*dt  B*dt  0 ] [  x0   ]\n+        #   [ u(dt) ] = exp [  0     0    I ] [  u0   ]\n+        #   [u1 - u0]       [  0     0    0 ] [u1 - u0]\n+        M = np.bmat([[A * dt, B * dt, np.zeros((n_states, n_inputs))],\n+                     [np.zeros((n_inputs, n_states + n_inputs)),\n+                      np.identity(n_inputs)],\n+                     [np.zeros((n_inputs, n_states + 2 * n_inputs))]])\n+        expMT = linalg.expm(transpose(M))\n+        Ad = expMT[:n_states, :n_states]\n+        Bd1 = expMT[n_states+n_inputs:, :n_states]\n+        Bd0 = expMT[n_states:n_states + n_inputs, :n_states] - Bd1\n+        for i in xrange(1, n_steps):\n+            xout[i] = (dot(xout[i-1], Ad) + dot(U[i-1], Bd0) +\n+                       dot(U[i], Bd1))","path":"scipy/signal/ltisys.py","position":163,"original_position":163,"commit_id":"1dbecec6b552512cc3689fefd0311520837f1dda","original_commit_id":"1dbecec6b552512cc3689fefd0311520837f1dda","user":{"login":"rgommers","id":98330,"avatar_url":"https://avatars.githubusercontent.com/u/98330?v=3","gravatar_id":"","url":"https://api.github.com/users/rgommers","html_url":"https://github.com/rgommers","followers_url":"https://api.github.com/users/rgommers/followers","following_url":"https://api.github.com/users/rgommers/following{/other_user}","gists_url":"https://api.github.com/users/rgommers/gists{/gist_id}","starred_url":"https://api.github.com/users/rgommers/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/rgommers/subscriptions","organizations_url":"https://api.github.com/users/rgommers/orgs","repos_url":"https://api.github.com/users/rgommers/repos","events_url":"https://api.github.com/users/rgommers/events{/privacy}","received_events_url":"https://api.github.com/users/rgommers/received_events","type":"User","site_admin":false},"body":"and style nitpick: please insert a blank line below and between consecutive if-statements in a few places above","created_at":"2015-04-12T15:37:27Z","updated_at":"2015-04-12T15:37:27Z","html_url":"https://github.com/scipy/scipy/pull/4675#discussion_r28204431","pull_request_url":"https://api.github.com/repos/scipy/scipy/pulls/4675","_links":{"self":{"href":"https://api.github.com/repos/scipy/scipy/pulls/comments/28204431"},"html":{"href":"https://github.com/scipy/scipy/pull/4675#discussion_r28204431"},"pull_request":{"href":"https://api.github.com/repos/scipy/scipy/pulls/4675"}}},"pull_request":{"url":"https://api.github.com/repos/scipy/scipy/pulls/4675","id":32188994,"html_url":"https://github.com/scipy/scipy/pull/4675","diff_url":"https://github.com/scipy/scipy/pull/4675.diff","patch_url":"https://github.com/scipy/scipy/pull/4675.patch","issue_url":"https://api.github.com/repos/scipy/scipy/issues/4675","number":4675,"state":"open","locked":false,"title":"Fix lsim","user":{"login":"cwrowley","id":2033951,"avatar_url":"https://avatars.githubusercontent.com/u/2033951?v=3","gravatar_id":"","url":"https://api.github.com/users/cwrowley","html_url":"https://github.com/cwrowley","followers_url":"https://api.github.com/users/cwrowley/followers","following_url":"https://api.github.com/users/cwrowley/following{/other_user}","gists_url":"https://api.github.com/users/cwrowley/gists{/gist_id}","starred_url":"https://api.github.com/users/cwrowley/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/cwrowley/subscriptions","organizations_url":"https://api.github.com/users/cwrowley/orgs","repos_url":"https://api.github.com/users/cwrowley/repos","events_url":"https://api.github.com/users/cwrowley/events{/privacy}","received_events_url":"https://api.github.com/users/cwrowley/received_events","type":"User","site_admin":false},"body":"This pull request provides a rewrite of `signal.lsim`, and fixes #2652 and #1104.\r\nThe new implementation fixes several known problems, in particular for systems with a zero eigenvalue, or with repeated eigenvalues.  It is also much faster than `lsim2`, which calls a generic ODE integrator.  See the commit messages for more details.","created_at":"2015-03-29T22:01:07Z","updated_at":"2015-04-12T15:37:27Z","closed_at":null,"merged_at":null,"merge_commit_sha":"c7b13aa6961dce2a8777012149752d679252ef5b","assignee":null,"milestone":null,"commits_url":"https://api.github.com/repos/scipy/scipy/pulls/4675/commits","review_comments_url":"https://api.github.com/repos/scipy/scipy/pulls/4675/comments","review_comment_url":"https://api.github.com/repos/scipy/scipy/pulls/comments{/number}","comments_url":"https://api.github.com/repos/scipy/scipy/issues/4675/comments","statuses_url":"https://api.github.com/repos/scipy/scipy/statuses/1dbecec6b552512cc3689fefd0311520837f1dda","head":{"label":"cwrowley:fix-lsim","ref":"fix-lsim","sha":"1dbecec6b552512cc3689fefd0311520837f1dda","user":{"login":"cwrowley","id":2033951,"avatar_url":"https://avatars.githubusercontent.com/u/2033951?v=3","gravatar_id":"","url":"https://api.github.com/users/cwrowley","html_url":"https://github.com/cwrowley","followers_url":"https://api.github.com/users/cwrowley/followers","following_url":"https://api.github.com/users/cwrowley/following{/other_user}","gists_url":"https://api.github.com/users/cwrowley/gists{/gist_id}","starred_url":"https://api.github.com/users/cwrowley/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/cwrowley/subscriptions","organizations_url":"https://api.github.com/users/cwrowley/orgs","repos_url":"https://api.github.com/users/cwrowley/repos","events_url":"https://api.github.com/users/cwrowley/events{/privacy}","received_events_url":"https://api.github.com/users/cwrowley/received_events","type":"User","site_admin":false},"repo":{"id":32708781,"name":"scipy","full_name":"cwrowley/scipy","owner":{"login":"cwrowley","id":2033951,"avatar_url":"https://avatars.githubusercontent.com/u/2033951?v=3","gravatar_id":"","url":"https://api.github.com/users/cwrowley","html_url":"https://github.com/cwrowley","followers_url":"https://api.github.com/users/cwrowley/followers","following_url":"https://api.github.com/users/cwrowley/following{/other_user}","gists_url":"https://api.github.com/users/cwrowley/gists{/gist_id}","starred_url":"https://api.github.com/users/cwrowley/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/cwrowley/subscriptions","organizations_url":"https://api.github.com/users/cwrowley/orgs","repos_url":"https://api.github.com/users/cwrowley/repos","events_url":"https://api.github.com/users/cwrowley/events{/privacy}","received_events_url":"https://api.github.com/users/cwrowley/received_events","type":"User","site_admin":false},"private":false,"html_url":"https://github.com/cwrowley/scipy","description":"Scipy library main repository","fork":true,"url":"https://api.github.com/repos/cwrowley/scipy","forks_url":"https://api.github.com/repos/cwrowley/scipy/forks","keys_url":"https://api.github.com/repos/cwrowley/scipy/keys{/key_id}","collaborators_url":"https://api.github.com/repos/cwrowley/scipy/collaborators{/collaborator}","teams_url":"https://api.github.com/repos/cwrowley/scipy/teams","hooks_url":"https://api.github.com/repos/cwrowley/scipy/hooks","issue_events_url":"https://api.github.com/repos/cwrowley/scipy/issues/events{/number}","events_url":"https://api.github.com/repos/cwrowley/scipy/events","assignees_url":"https://api.github.com/repos/cwrowley/scipy/assignees{/user}","branches_url":"https://api.github.com/repos/cwrowley/scipy/branches{/branch}","tags_url":"https://api.github.com/repos/cwrowley/scipy/tags","blobs_url":"https://api.github.com/repos/cwrowley/scipy/git/blobs{/sha}","git_tags_url":"https://api.github.com/repos/cwrowley/scipy/git/tags{/sha}","git_refs_url":"https://api.github.com/repos/cwrowley/scipy/git/refs{/sha}","trees_url":"https://api.github.com/repos/cwrowley/scipy/git/trees{/sha}","statuses_url":"https://api.github.com/repos/cwrowley/scipy/statuses/{sha}","languages_url":"https://api.github.com/repos/cwrowley/scipy/languages","stargazers_url":"https://api.github.com/repos/cwrowley/scipy/stargazers","contributors_url":"https://api.github.com/repos/cwrowley/scipy/contributors","subscribers_url":"https://api.github.com/repos/cwrowley/scipy/subscribers","subscription_url":"https://api.github.com/repos/cwrowley/scipy/subscription","commits_url":"https://api.github.com/repos/cwrowley/scipy/commits{/sha}","git_commits_url":"https://api.github.com/repos/cwrowley/scipy/git/commits{/sha}","comments_url":"https://api.github.com/repos/cwrowley/scipy/comments{/number}","issue_comment_url":"https://api.github.com/repos/cwrowley/scipy/issues/comments{/number}","contents_url":"https://api.github.com/repos/cwrowley/scipy/contents/{+path}","compare_url":"https://api.github.com/repos/cwrowley/scipy/compare/{base}...{head}","merges_url":"https://api.github.com/repos/cwrowley/scipy/merges","archive_url":"https://api.github.com/repos/cwrowley/scipy/{archive_format}{/ref}","downloads_url":"https://api.github.com/repos/cwrowley/scipy/downloads","issues_url":"https://api.github.com/repos/cwrowley/scipy/issues{/number}","pulls_url":"https://api.github.com/repos/cwrowley/scipy/pulls{/number}","milestones_url":"https://api.github.com/repos/cwrowley/scipy/milestones{/number}","notifications_url":"https://api.github.com/repos/cwrowley/scipy/notifications{?since,all,participating}","labels_url":"https://api.github.com/repos/cwrowley/scipy/labels{/name}","releases_url":"https://api.github.com/repos/cwrowley/scipy/releases{/id}","created_at":"2015-03-23T03:18:21Z","updated_at":"2015-03-23T03:18:29Z","pushed_at":"2015-03-29T21:53:32Z","git_url":"git://github.com/cwrowley/scipy.git","ssh_url":"git@github.com:cwrowley/scipy.git","clone_url":"https://github.com/cwrowley/scipy.git","svn_url":"https://github.com/cwrowley/scipy","homepage":"http://scipy.org/scipylib/","size":54105,"stargazers_count":0,"watchers_count":0,"language":"Python","has_issues":false,"has_downloads":true,"has_wiki":true,"has_pages":false,"forks_count":0,"mirror_url":null,"open_issues_count":0,"forks":0,"open_issues":0,"watchers":0,"default_branch":"master"}},"base":{"label":"scipy:master","ref":"master","sha":"52354d4b0e9049ad9df86de934ed340741dd40a8","user":{"login":"scipy","id":288277,"avatar_url":"https://avatars.githubusercontent.com/u/288277?v=3","gravatar_id":"","url":"https://api.github.com/users/scipy","html_url":"https://github.com/scipy","followers_url":"https://api.github.com/users/scipy/followers","following_url":"https://api.github.com/users/scipy/following{/other_user}","gists_url":"https://api.github.com/users/scipy/gists{/gist_id}","starred_url":"https://api.github.com/users/scipy/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/scipy/subscriptions","organizations_url":"https://api.github.com/users/scipy/orgs","repos_url":"https://api.github.com/users/scipy/repos","events_url":"https://api.github.com/users/scipy/events{/privacy}","received_events_url":"https://api.github.com/users/scipy/received_events","type":"Organization","site_admin":false},"repo":{"id":1460385,"name":"scipy","full_name":"scipy/scipy","owner":{"login":"scipy","id":288277,"avatar_url":"https://avatars.githubusercontent.com/u/288277?v=3","gravatar_id":"","url":"https://api.github.com/users/scipy","html_url":"https://github.com/scipy","followers_url":"https://api.github.com/users/scipy/followers","following_url":"https://api.github.com/users/scipy/following{/other_user}","gists_url":"https://api.github.com/users/scipy/gists{/gist_id}","starred_url":"https://api.github.com/users/scipy/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/scipy/subscriptions","organizations_url":"https://api.github.com/users/scipy/orgs","repos_url":"https://api.github.com/users/scipy/repos","events_url":"https://api.github.com/users/scipy/events{/privacy}","received_events_url":"https://api.github.com/users/scipy/received_events","type":"Organization","site_admin":false},"private":false,"html_url":"https://github.com/scipy/scipy","description":"Scipy library main repository","fork":false,"url":"https://api.github.com/repos/scipy/scipy","forks_url":"https://api.github.com/repos/scipy/scipy/forks","keys_url":"https://api.github.com/repos/scipy/scipy/keys{/key_id}","collaborators_url":"https://api.github.com/repos/scipy/scipy/collaborators{/collaborator}","teams_url":"https://api.github.com/repos/scipy/scipy/teams","hooks_url":"https://api.github.com/repos/scipy/scipy/hooks","issue_events_url":"https://api.github.com/repos/scipy/scipy/issues/events{/number}","events_url":"https://api.github.com/repos/scipy/scipy/events","assignees_url":"https://api.github.com/repos/scipy/scipy/assignees{/user}","branches_url":"https://api.github.com/repos/scipy/scipy/branches{/branch}","tags_url":"https://api.github.com/repos/scipy/scipy/tags","blobs_url":"https://api.github.com/repos/scipy/scipy/git/blobs{/sha}","git_tags_url":"https://api.github.com/repos/scipy/scipy/git/tags{/sha}","git_refs_url":"https://api.github.com/repos/scipy/scipy/git/refs{/sha}","trees_url":"https://api.github.com/repos/scipy/scipy/git/trees{/sha}","statuses_url":"https://api.github.com/repos/scipy/scipy/statuses/{sha}","languages_url":"https://api.github.com/repos/scipy/scipy/languages","stargazers_url":"https://api.github.com/repos/scipy/scipy/stargazers","contributors_url":"https://api.github.com/repos/scipy/scipy/contributors","subscribers_url":"https://api.github.com/repos/scipy/scipy/subscribers","subscription_url":"https://api.github.com/repos/scipy/scipy/subscription","commits_url":"https://api.github.com/repos/scipy/scipy/commits{/sha}","git_commits_url":"https://api.github.com/repos/scipy/scipy/git/commits{/sha}","comments_url":"https://api.github.com/repos/scipy/scipy/comments{/number}","issue_comment_url":"https://api.github.com/repos/scipy/scipy/issues/comments{/number}","contents_url":"https://api.github.com/repos/scipy/scipy/contents/{+path}","compare_url":"https://api.github.com/repos/scipy/scipy/compare/{base}...{head}","merges_url":"https://api.github.com/repos/scipy/scipy/merges","archive_url":"https://api.github.com/repos/scipy/scipy/{archive_format}{/ref}","downloads_url":"https://api.github.com/repos/scipy/scipy/downloads","issues_url":"https://api.github.com/repos/scipy/scipy/issues{/number}","pulls_url":"https://api.github.com/repos/scipy/scipy/pulls{/number}","milestones_url":"https://api.github.com/repos/scipy/scipy/milestones{/number}","notifications_url":"https://api.github.com/repos/scipy/scipy/notifications{?since,all,participating}","labels_url":"https://api.github.com/repos/scipy/scipy/labels{/name}","releases_url":"https://api.github.com/repos/scipy/scipy/releases{/id}","created_at":"2011-03-09T18:52:03Z","updated_at":"2015-04-12T09:29:41Z","pushed_at":"2015-04-12T13:57:19Z","git_url":"git://github.com/scipy/scipy.git","ssh_url":"git@github.com:scipy/scipy.git","clone_url":"https://github.com/scipy/scipy.git","svn_url":"https://github.com/scipy/scipy","homepage":"http://scipy.org/scipylib/","size":86236,"stargazers_count":1582,"watchers_count":1582,"language":"Python","has_issues":true,"has_downloads":true,"has_wiki":true,"has_pages":false,"forks_count":853,"mirror_url":null,"open_issues_count":800,"forks":853,"open_issues":800,"watchers":1582,"default_branch":"master"}},"_links":{"self":{"href":"https://api.github.com/repos/scipy/scipy/pulls/4675"},"html":{"href":"https://github.com/scipy/scipy/pull/4675"},"issue":{"href":"https://api.github.com/repos/scipy/scipy/issues/4675"},"comments":{"href":"https://api.github.com/repos/scipy/scipy/issues/4675/comments"},"review_comments":{"href":"https://api.github.com/repos/scipy/scipy/pulls/4675/comments"},"review_comment":{"href":"https://api.github.com/repos/scipy/scipy/pulls/comments{/number}"},"commits":{"href":"https://api.github.com/repos/scipy/scipy/pulls/4675/commits"},"statuses":{"href":"https://api.github.com/repos/scipy/scipy/statuses/1dbecec6b552512cc3689fefd0311520837f1dda"}}}}