{"action":"created","issue":{"url":"https://api.github.com/repos/desihub/desisim/issues/357","repository_url":"https://api.github.com/repos/desihub/desisim","labels_url":"https://api.github.com/repos/desihub/desisim/issues/357/labels{/name}","comments_url":"https://api.github.com/repos/desihub/desisim/issues/357/comments","events_url":"https://api.github.com/repos/desihub/desisim/issues/357/events","html_url":"https://github.com/desihub/desisim/issues/357","id":319257814,"node_id":"MDU6SXNzdWUzMTkyNTc4MTQ=","number":357,"title":"chi2 and amplitude computation is slow and chi2 is wrong","user":{"login":"londumas","id":3166436,"node_id":"MDQ6VXNlcjMxNjY0MzY=","avatar_url":"https://avatars0.githubusercontent.com/u/3166436?v=4","gravatar_id":"","url":"https://api.github.com/users/londumas","html_url":"https://github.com/londumas","followers_url":"https://api.github.com/users/londumas/followers","following_url":"https://api.github.com/users/londumas/following{/other_user}","gists_url":"https://api.github.com/users/londumas/gists{/gist_id}","starred_url":"https://api.github.com/users/londumas/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/londumas/subscriptions","organizations_url":"https://api.github.com/users/londumas/orgs","repos_url":"https://api.github.com/users/londumas/repos","events_url":"https://api.github.com/users/londumas/events{/privacy}","received_events_url":"https://api.github.com/users/londumas/received_events","type":"User","site_admin":false},"labels":[],"state":"closed","locked":false,"assignee":null,"assignees":[],"milestone":null,"comments":3,"created_at":"2018-05-01T17:00:31Z","updated_at":"2018-06-29T14:39:01Z","closed_at":"2018-06-29T14:39:01Z","author_association":"CONTRIBUTOR","body":"https://github.com/desihub/desisim/blob/36b78ee41ad5d0f7e2bb4bc4b2058e761eeb3260/py/desisim/archetypes.py#L52\r\n\r\nThis line is very slow, a simple fix is using something like what follows.\r\nWe need to see how to add the error contribution.\r\nThis uses Bolton et al. 2012 (https://arxiv.org/pdf/1207.7326.pdf), see section 4.3.2.\r\nThis gives the same result for the amplitude but a factor 2 difference for the chi2.\r\nIt seems that my chi2 is good and that the one of `mathutils.quick_amplitude` is not.\r\n```\r\ndef my_compute_chi2(flux, ferr=None):  \r\n\r\n    nspec, npix = flux.shape\r\n\r\n    chi2 = np.zeros((nspec, nspec)).astype('f4')\r\n    amp = np.zeros((nspec, nspec)).astype('f4')\r\n\r\n    rescale = np.sqrt(npix/np.sum(flux**2,axis=1))\r\n    flux *= rescale[:,None]\r\n    if ferr is None:\r\n        ferr = np.ones(flux.shape)\r\n    else:\r\n        ferr *= rescale[:,None]\r\n\r\n    for ii in range(nspec):\r\n        if ii % 500 == 0 or ii == 0:\r\n            log.info('Computing chi2 matrix for spectra {}-{} out of {}.'.format(\r\n            (ii % 500) * 500, np.min(((ii+1) * (ii % 500), nspec-1)), nspec))\r\n        amp1 = np.sum(flux[ii]*flux,axis=1)/npix\r\n        chi21 = npix*(1.-amp1**2)\r\n        chi21[np.isclose(amp1,1.)] = 0.\r\n        chi2[ii, :] = chi21\r\n        amp[ii, :] = amp1\r\n\r\n    return chi2, amp\r\n```\r\n\r\nHere is a simple code to compare both results\r\n```\r\nimport scipy as sp\r\nimport matplotlib.pyplot as plt\r\n\r\nfrom desisim.archetypes import ArcheTypes, compute_chi2\r\n\r\ndef my_compute_chi2(flux):  \r\n\r\n    nspec, npix = flux.shape\r\n    ferr = sp.ones_like(flux)\r\n\r\n    chi2 = sp.zeros((nspec, nspec)).astype('f4')\r\n    amp = sp.zeros((nspec, nspec)).astype('f4')\r\n\r\n    for ii in range(nspec):\r\n        if ii % 500 == 0 or ii == 0:\r\n            print('Computing chi2 matrix for spectra {}-{} out of {}.'.format((ii % 500) * 500, min(((ii+1) * (ii % 500), nspec-1)), nspec))\r\n        amp1 = sp.sum(flux[ii]*flux,axis=1)/npix\r\n        chi21 = npix*(1.-amp1**2)\r\n        chi21[sp.isclose(amp1,1.)] = 0.\r\n        chi2[ii, :] = chi21\r\n        amp[ii, :] = amp1\r\n\r\n    return chi2, amp\r\n\r\ndef main():\r\n\r\n    ###\r\n    nspec = 2\r\n    npix = 100\r\n    flux = sp.zeros((nspec, npix))\r\n    flux[0,:] = 1.\r\n    flux[1,:] = 1.*sp.arange(npix)\r\n    for i in range(nspec):\r\n        flux[i,:] *= sp.sqrt(npix/sp.sum(flux[i,:]**2))\r\n\r\n    ###\r\n    print(\"\")\r\n    chi2, amp = compute_chi2(flux)\r\n    print(\"Chi2\")\r\n    print(chi2)\r\n    print(\"Amp\")\r\n    print(amp)\r\n\r\n    ###\r\n    print(\"\")\r\n    chi2, amp = my_compute_chi2(flux)\r\n    print(\"Chi2\")\r\n    print(chi2)\r\n    print(\"Amp\")\r\n    print(amp)\r\n\r\n    ### Chi2 by hand\r\n    chi2 = sp.sum( (flux[1,:]-amp[1,0]*flux[0,:])**2 )\r\n    print(chi2)\r\n\r\n    plt.plot(amp[1,0]*flux[0,:])\r\n    plt.plot(flux[1,:])\r\n    plt.grid()\r\n    plt.show()\r\n\r\n    return\r\n\r\nmain()\r\n```"},"comment":{"url":"https://api.github.com/repos/desihub/desisim/issues/comments/401374436","html_url":"https://github.com/desihub/desisim/issues/357#issuecomment-401374436","issue_url":"https://api.github.com/repos/desihub/desisim/issues/357","id":401374436,"node_id":"MDEyOklzc3VlQ29tbWVudDQwMTM3NDQzNg==","user":{"login":"londumas","id":3166436,"node_id":"MDQ6VXNlcjMxNjY0MzY=","avatar_url":"https://avatars0.githubusercontent.com/u/3166436?v=4","gravatar_id":"","url":"https://api.github.com/users/londumas","html_url":"https://github.com/londumas","followers_url":"https://api.github.com/users/londumas/followers","following_url":"https://api.github.com/users/londumas/following{/other_user}","gists_url":"https://api.github.com/users/londumas/gists{/gist_id}","starred_url":"https://api.github.com/users/londumas/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/londumas/subscriptions","organizations_url":"https://api.github.com/users/londumas/orgs","repos_url":"https://api.github.com/users/londumas/repos","events_url":"https://api.github.com/users/londumas/events{/privacy}","received_events_url":"https://api.github.com/users/londumas/received_events","type":"User","site_admin":false},"created_at":"2018-06-29T14:39:01Z","updated_at":"2018-06-29T14:39:01Z","author_association":"CONTRIBUTOR","body":"fixed"}}
{"actor":{"display_login":"londumas"}}