diff --git a/.travis.yml b/.travis.yml index c103978..d7cbf4c 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,6 +1,12 @@ language: python +sudo: required + + python: - 2.7 + - 3.5 + - 3.6 + notifications: email: recipients: @@ -14,9 +20,15 @@ deploy: password: secure: lAM0ScaCdOTbixtJf908INAKIpVCIPLSu4ZpqSxdp0zzMsTJJoTLo00wlmzzjeQDsxqgpiUIvjnwlpdur7Q+IPmwFa9hkg4VX4Mqe0keoDaoUQaUzpypZqYF0w4egVgc62NnXCGg3B/JuqUFuyM692uvmjvokGmadXLevXvPaxM= +addons: + apt: + sources: + - ubuntu-toolchain-r-test + packages: + - gfortran-6 + # Setup anaconda before_install: - - sudo apt-get install gfortran - if [[ "$TRAVIS_PYTHON_VERSION" == "2.7" ]]; then wget https://repo.continuum.io/miniconda/Miniconda-latest-Linux-x86_64.sh -O miniconda.sh; else @@ -31,16 +43,17 @@ before_install: - conda info -a install: - - conda install python=$TRAVIS_PYTHON_VERSION numpy scipy mpmath nose astropy + - conda install python=$TRAVIS_PYTHON_VERSION numpy scipy mpmath nose astropy libgfortran - pip install python-coveralls + - pip install coverage - pip install emcee - - pip install mrpy - - CAMBURL=http://camb.info/CAMB_Mar13.tar.gz pip install git+git://github.com/steven-murray/pycamb.git - - pip install cached_property + #- CAMBURL=http://camb.info/CAMB_Mar13.tar.gz pip install git+git://github.com/steven-murray/pycamb.git + - sudo ln -s /usr/bin/gfortran-6 /usr/bin/gfortran + - pip install --egg camb - pip install coveralls script: - - nosetests --with-coverage --cover-package=hmf + - coverage run --source=hmf --omit=hmf/fitting/* $(which nosetests) # Calculate coverage diff --git a/CHANGELOG.rst b/CHANGELOG.rst index eb10ea3..bbdad83 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -1,9 +1,32 @@ Releases ======== +v3.0.0 [7th June 2017] +---------------------- +**Features** + +- Now provides compatibility with Python 3.x. Support for 2.x will be removed in hmf v3.1 (whenever that comes). +- Complete overhaul of the caching system. Should be invisible to the user, but streamlines writing of framework + code considerably. Removes required manual specification of dependencies between quantities, and adds ability + to specify parameter kinds (model, param, res, option, switch). + +**Bugfixes** + +- Fixed bug in Caroll1992 GrowthFactor class which affected high-redshift growth factors (thanks to Olmo Piana). +- Fixed astropy dependency to be >= 1.1 +- Fixed bug where Takahashi parameters were always passed through regardess of ``takahashi`` setting. +- Fixed small bug where the functional.get_label method returned differently ordered parameters because of dicts. +- Note that the fitting subpackage is temporarily unsupported and I discourage its use for the time being. + +**Enhancement** + +- Completely removes dependence on archaic pycamb package. Now supports natively supplied python interface to CAMB. + Install camb with ``pip install --egg camb``. This means that much more modern versions of CAMB can be used. +- Many new tests, to bring total coverage up to >80%, and continuous testing on Python 2.7, 3.5 and 3.6 + + v2.0.5 [12th January 2017] -------------------------- - **Bugfixes** - Fixed bug in GrowthFactor which gave ripples in functions of z when a coarse grid was used. Thanks to @mirochaj and @@ -15,6 +38,7 @@ v2.0.5 [12th January 2017] - Totally removed dependency on the Cache (super)class -- caching indexes now inherent to the called class. - More robust parameter information based on introspection. + v2.0.4 [11th November, 2016] ------ diff --git a/README.rst b/README.rst index b8e2e86..ae3993f 100644 --- a/README.rst +++ b/README.rst @@ -12,11 +12,15 @@ hmf :target: https://pypi.python.org/pypi/hmf .. image:: https://coveralls.io/repos/github/steven-murray/hmf/badge.svg?branch=master :target: https://coveralls.io/github/steven-murray/hmf?branch=master +.. image:: https://img.shields.io/pypi/pyversions/hmf.svg `hmf` is a python application that provides a flexible and simple way to calculate the Halo Mass Function for a range of varying parameters. It is also the backend to `HMFcalc `_, the online HMF calculator. +.. warning:: Due to the general trend of moving to Python 3 by important projects such as IPython and astropy, from + version 2.1, hmf is compatible with Python 3, and from version 3.0, it will drop support for Python 2. + Documentation ------------- `Read the docs. `_ diff --git a/development/.gitignore b/development/.gitignore new file mode 100644 index 0000000..f47cb20 --- /dev/null +++ b/development/.gitignore @@ -0,0 +1 @@ +*.out diff --git a/development/Plots/.gitignore b/development/Plots/.gitignore new file mode 100644 index 0000000..5b7e248 --- /dev/null +++ b/development/Plots/.gitignore @@ -0,0 +1,2 @@ +*.jpg +*.mp4 diff --git a/development/halofit_testing.ipynb b/development/halofit_testing.ipynb new file mode 100644 index 0000000..a72fab5 --- /dev/null +++ b/development/halofit_testing.ipynb @@ -0,0 +1,213 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "ExecuteTime": { + "end_time": "2017-05-22T05:12:46.410124", + "start_time": "2017-05-22T05:12:46.050906Z" + }, + "collapsed": false + }, + "outputs": [], + "source": [ + "from hmf.transfer import Transfer" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "ExecuteTime": { + "end_time": "2017-05-22T05:54:50.276527", + "start_time": "2017-05-22T05:54:50.269658Z" + }, + "collapsed": false + }, + "outputs": [], + "source": [ + "teh_nl_tk = Transfer(transfer_model=\"EH\", takahashi=True, z= 8.0)\n", + "teh_nl_ntk = Transfer(transfer_model=\"EH\", takahashi=False, z=8.0)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "ExecuteTime": { + "end_time": "2017-05-22T05:54:52.168860", + "start_time": "2017-05-22T05:54:52.164999Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "ExecuteTime": { + "end_time": "2017-05-22T05:56:18.289893", + "start_time": "2017-05-22T05:56:18.287974Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "teh_nl_tk.update(z=0)" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": { + "ExecuteTime": { + "end_time": "2017-05-22T05:57:46.801646", + "start_time": "2017-05-22T05:57:46.525151Z" + }, + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEACAYAAACtVTGuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8VfWZ7/HPNwkJhPs1cr9oEAPeCgV1BobWtqKnp7a2\nWuxobYtSptrazmumSuc14znTsWMvdg4z1VJandp6oVZry1havBG1FQHRgISLxnBJQggYSEIuEBKe\n88deoduY7OwkO+zslef9eu1X9vr91m+t52HHPK7Lby+ZGc4551x70pIdgHPOud7NC4VzzrmYvFA4\n55yLyQuFc865mLxQOOeci8kLhXPOuZi8UDjnnIvJC4VzzrmYvFA455yLyQuFc865mDKSHUAijBo1\nyqZMmZLsMNpVV1fHwIEDkx1Gt4UhjzDkAOHIIww5QGrnsWXLlnfNbHRH64WiUEyZMoXXXnst2WG0\nKz8/n4ULFyY7jG4LQx5hyAHCkUcYcoDUzkPSvnjW81NPzjnnYvJC4ZxzLqYOC4WkByUdkrS9nX5J\n+k9JRZK2SfpAVN8iSbuDvjuj2kdIelbS28HP4VF9y4P1d0u6orsJOuec6554jih+DiyK0X8lkBu8\nlgI/BpCUDtwX9OcB10vKC8bcCTxvZrnA88EyQf9iYGawz/uD7TjnnEuSDguFmb0EHImxytXALyzi\nVWCYpLHAXKDIzIrNrBFYHazbMuah4P1DwCej2leb2Qkz2wMUBdtxzjmXJIm462k8UBK1XBq0tdU+\nL3ifY2blwfuDQE7Utl5tY1vvI2kpkSMYcnJyyM/P73oGPay2trZXxxevMOQRhhwgHHmEIQcITx6x\nJP32WDMzSZ1+HquZrQJWAcyZM8d68+1pqXz7XLQw5BGGHCAceYQhB0hOHs2njLKjDew/Us/ArHQu\nnjS840HdkIhCUQZMjFqeELT1a6cdoELSWDMrD05THepgW8451+c0nzKKD9ey/UA1b1XUUny4luLD\ndeyrrKex+RQAV51/Fvf/7ewejSMRhWINcJuk1UROLVUHBeAwkCtpKpE/9ouBz0WNuQm4J/j5u6j2\nRyX9EBhH5AL5pgTE6Jxzvd6RukY27alk454jbCutZseBGhpONgPQL11MGpHNtNGD+PCMMUwdNZDJ\nIwcydVTPzwrvsFBIegxYCIySVArcReRoATNbCawFriJy4bke+GLQ1yTpNmAdkA48aGaFwWbvAR6X\ntATYB1wXjCmU9DiwA2gCbjWz5sSk6pzrLcyM401GVX0j2ZkZZGb0zSldzaeM1/Ye4dkdFfyp6F12\nHTwGQP9+aZw/fiiL505k1rihnD9hKNNGDSQjPTn/Th0WCjO7voN+A25tp28tkULSur0SuLydMXcD\nd3cUl3MuNVTUHOfV4kpe33eUXQePsbeyjkPHTmAGPPcsANmZ6eSOGcR5Y4ewYPpo/mb6aAZmJf0S\nao9oPmX8qehdnt56gOd3HeJIXSOZGWnMnTKCf7xiHJdMG8H544f1quIZzk/COZdUR+oaeWJLCb9/\n8yBbS6oAGJiZzoyxQ/jrc0Yzblh/Ksr2MSM3l/rGJirrGnmr4hh/2H6Q1ZtLGJSVwfVzJ3LLgmmM\nGdw/ydkkRsmReh5/rYQntpRSXn2cwf0zuHzGGD428ywWTB/NoF5cGHtvZM65lFNypJ4fvVDEUwVl\nNDad4vzxQ/nHK87lb6aPZsZZg99z6iQ/v5yFfz31PeObmk+xee9RHtu0nwf/vJfVm0v45hXncsMl\nk5F0ptNJiDdLq1n50jv84c3IjIAF00fzzx/P4/LzxpCVkRrzib1QOOe6rb6xiR+9UMTP/rSHNMG1\nsydw02VTmJ4zuFPbyUhP49KzR3Lp2SP5+kdy+ZffFfLPvytkQ3El3//MhSl1OmprSRXfW7eLPxdV\nMjgrg6ULzuamyyYzduiAZIfWaanzr+6c65UKD1Tz1cfeoPhwHddcPJ5/XHRuQv4YThs9iF8umcvP\nXt7Dv/9hJ/uP1PPIkksYmt0vAVH3nJIj9Xz3j7t4els5IwdmcueVM/jcvEkM6d+7447FC4VzrsvW\nFR7kq4+9wfDsfjx6yzwuO3tUQrcviVsWTOPsMQNZ9svXueGBjTx88zyGDuh9f3RPNp/iZy/v4f89\n9xZpEl/78DncsmAag1O4QLToPZfVnXMp5Yktpfzdw1vIGzuEtV+bn/AiEe3DM3L4yY2z2XWwhtse\nfZ2mYLJZb1F4oJpP/OjPfPePu/jQuWNY/w8L+fuPnRuKIgFeKJxzXbB+9yHueHIbl509ikdvmcfI\nQVk9vs8PzRjDv31yFi+//S7fWburx/cXDzPj2X0n+dR9r1BZe4Kf3DiblTfO5qyh4bhTq4WfenLO\ndcpbFce49ZHXmXHWYFbeOJvszDP3Z+SzH5zEzvJjPPjnPVx29kg+kpfT8aAeUnP8JH//q608t7OR\ny2eM4fvXXsiIgZlJi6cn+RGFcy5uJ5qauX11AQP6pfPfX/hgUu79X37VDPLGDuGOJ7dx+NiJM75/\ngH2VdVxz/yvk7z7E52Zk8rOb5oS2SIAXCudcJ6x47m12ltdwz6cvYMyQ5JxeycpIZ8Xii6g90cQ/\n/7bNB2/2qI3FlXzyvj/zbu0JfrlkHh+b0i9l53jEywuFcy4ue96tY9VLxVw7ewIfTeIpH4DcnMHc\n/pFc/lh4kOd3Vpyx/T67o4IbH9zE8IGZ/PYrf8WlZ488Y/tOJi8Uzrm4fO+Pu8jMSOObi2YkOxQA\nbpk/jek5g/iX3xVS39jU4/v7XUEZyx7ewnljh/DkssuYcga+tbW38ELhnOvQln1H+cP2g3x5wdmM\nHtzzdzjFo196Gnd/6nzKqhpY8dzbPbqvxzeX8PVfFTB3yggeuXkew0N8PaItXiiccx36cX4RIwdm\ncsuCqR2vfAZ9cMoIrpszgQf+tIe3K471yD7+Z+sB7vjNNhbkjua/v5icC/jJ5oXCORfTvso6nt91\niL+dN+mM3gobrzsWzSA7M5271hQSeepB4rywq4Jv/KqAD04ewcobZtO/X2p8iV+ixVUoJC2StFtS\nkaQ72+gfLukpSdskbZI0K6rvdknbJRVK+npU+68kFQSvvZIKgvYpkhqi+lYmIlHnXNf8/JW9ZKSJ\nGy6ZnOxQ2jRyUBb/cMW5vPJOJb8PvqE1ETa8U8nfPfw6540dwgNfmMOAzL5ZJCCOQiEpHbgPuBLI\nA66XlNdqtW8BBWZ2AfB5YEUwdhZwCzAXuBD4uKRzAMzss2Z2kZldBDwJ/CZqe++09JnZsm5l6Jzr\nsroTTfz6tVL+1/ljk3Y7bDz+dt5k8sYO4d+e3kndie5f2C4oqeLmhzYzaUQ2D31pbmi+iqOr4jmi\nmAsUmVmxmTUCq4GrW62TB7wAYGa7gCmScoDzgI1mVm9mTcCLwDXRAxW5Afk64LFuZeKcS7jndlZQ\ne6KJz83rnUcTLdLTxLc/OZODNcf5rxeKurWtXQdruOnBTYwclMXDN88L9US6eMVTKMYDJVHLpUFb\ntK0EBUDSXGAyMAHYDsyXNFJSNpFna09sNXY+UGFm0bctTA1OO70oaX7c2TjnEmpNwQHGDu3PnMnD\nkx1Kh2ZPHsGnPzCBB/5UzDuHa7u0jT3v1nHDzzbRv18aj9w8j5xefBR1Jqmjiz+SPgMsMrObg+Ub\ngXlmdlvUOkOInG66GHgTmAHcYmYFkpYAXwHqgELghJlFX6v4MZEjlnuD5SxgkJlVSpoN/BaYaWY1\nreJaCiwFyMnJmb169epu/DP0rNraWgYNGpTsMLotDHmEIQc4M3nUNhq3r6/no5P7sXhG4v+vuidy\nqD5h3PlyPdOGpvEPc/p3asb0uw2n+M7G4zQ2G8vnDWD8oPju9Unl36kPfehDW8xsTocrmlnMF3Ap\nsC5qeTmwPMb6AvYCQ9ro+w7wlajlDKACmBBje/nAnFgxzp4923qz9evXJzuEhAhDHmHIwezM5PHY\nxn02+Y6n7c3Sqh7Zfk/l8NAre2zyHU/bL17ZE/eYg9UNNv+7L9isu/7Y6XxT+XcKeM06qAFmFtep\np81ArqSpkjKBxcCa6BUkDQv6AG4GXrLgCEDSmODnJCKnpx6NGvoRYJeZlUZta3RwAR1J04BcoDiO\nOJ1zCfT7N8uZOmogM8cNSXYonXLDvMn8zfTRfPv3O9lZXtPh+gerj/O5n75KZe0JHvrSXGaNH3oG\nokwtHRYKi1yEvg1YB+wEHjezQknLJLXckXQesF3SbiJ3R90etYknJe0A/ge41cyqovoW8/6L2AuA\nbcHtsk8Ay8zsSBdyc851UUNjMxv3HOHyGWNS7gvv0tLEvdddyNAB/bjlF69RUXO83XWLDh3j0z9+\nhYPVx3nwCx/kA5N6/7WYZIhr9oyZrQXWtmpbGfV+AzC9nbHtXow2sy+00fYkkdtlnXNJ8uqeShqb\nTrFg+uhkh9IlowZl8cBNc7h+1at89icb+NlNH+ScMX+5jmBmPP5aCf/3f3aQnZnOr758qR9JxND7\nplk655LupbcOk5WRxtypI5IdSpddMGEYv1gyl6W/2MJVK17m2jkTuHjScI7UnWDN1gNsL6vhsrNH\ncu91FzJ26IBkh9ureaFwzr3PS28d5pJpI1P+KytmTx7B7782n3uf2c2Tr5fyyMb9AEzPGcT3Pn0B\nn5k9gbS01Dq1lgxeKJxz71F6tJ53Dtf1+kl28TpraH++f+2F/NunZnGg6jiDsjJ6zTfgpgovFM65\n93jlnUoA5ueOSnIkiZWVkc7UPvQMiUTyb491zr3H6/uOMiy7H+eMTs1JZC7xvFA4595jy76jfGDS\ncD93707zQuGcO626/iRvH6pldgp8t5M7c7xQOOdOe73kKIBPPHPv4YXCOXfa6/uOkp4mLpzok8/c\nX3ihcM6dtmXfUfLGDumVjzx1yeOFwjkHwKlTxtaSKi6eNCzZobhexguFcw6AfUfqqWts9u88cu/j\nhcI5B8COA5Gv5M4bm1pfK+56nhcK5xwAO8qryUgTuTk+0c69lxcK5xwQOaI4Z8wgsjJS+4sAXeJ5\noXDOAbCjvIa8FHuanTsz4ioUkhZJ2i2pSNKdbfQPl/SUpG2SNkmaFdV3u6TtkgolfT2q/f9IKpNU\nELyuiupbHuxrt6Qrupukcy62d2tPUFFzwq9PuDZ1WCiC51ffR+QRp3nA9ZLyWq32LaDAzC4APg+s\nCMbOAm4B5gIXAh+XdE7UuP8ws4uC19pgTB6RR6TOBBYB97c8Q9s51zNani3tRxSuLfEcUcwFisys\n2MwagdXA1a3WyQNeADCzXcAUSTlEnqW90czqg2dvvwhc08H+rgZWm9kJM9sDFAUxOOd6iN/x5GKJ\np1CMB0qilkuDtmhbCQqApLnAZGACsB2YL2mkpGzgKmBi1LivBqerHpTU8uUy8ezPOZdAuyuOcdaQ\n/gzLzkx2KK4XStQ8/XuAFZIKgDeBN4BmM9sp6bvAM0AdUAA0B2N+DHwbsODnvcCX4t2hpKXAUoCc\nnBzy8/MTk0kPqK2t7dXxxSsMeYQhB0h8Hm+808DIfpzRfxv/LFJHPIWijPceBUwI2k4zsxrgiwCS\nBOwBioO+B4AHgr7vEDlCwMwqWsZL+inwdLz7C8avAlYBzJkzxxYuXBhHKsmRn59Pb44vXmHIIww5\nQGLzMDMOr3+GT39gPAsXzup4QIL4Z5E64jn1tBnIlTRVUiaRC81roleQNCzoA7gZeCkoHkgaE/yc\nROT01KPB8tioTXyKyGkqgm0vlpQlaSqQC2zqSnLOuY5V1Jyg9kQT54zxiXaubR0eUZhZk6TbgHVA\nOvCgmRVKWhb0ryRy0fohSQYUAkuiNvGkpJHASeBWM6sK2r8n6SIip572Al8Otlco6XFgB9AUjGnG\nOdcjig7VAnC2P/rUtSOuaxTBratrW7WtjHq/AZjeztj57bTfGGN/dwN3xxObc6573jkcKRR+ROHa\n4zOznevjig7VMrh/BqMHZyU7FNdLeaFwro9753AtZ48eROQ+FOfezwuFc31c0aFaP+3kYvJC4Vwf\nVnP8JIeOnfAL2S4mLxTO9WF7DtcBMG30wCRH4nozLxTO9WH7jtQDMGWkFwrXPi8UzvVh+ysjRxST\nRmQnORLXm3mhcK4P21tZz5jBWQzI9G/yd+3zQuFcH7a/st5PO7kOeaFwrg/bd6SOSSP9tJOLzQuF\nc33U8ZPNVNScYLJfn3Ad8ELhXB+1P7jjafIoP/XkYvNC4VwftffdyB1PfkThOuKFwrk+6vQRhV+j\ncB3wQuFcH7Wvsp4h/TP8OdmuQ3EVCkmLJO2WVCTpzjb6h0t6StI2SZskzYrqu13SdkmFkr4e1f59\nSbuCMU9JGha0T5HUIKkgeK1svT/nXPftO1LPFL8+4eLQYaGQlA7cB1wJ5AHXS8prtdq3gAIzuwD4\nPLAiGDsLuAWYC1wIfFzSOcGYZ4FZwZi3gOVR23vHzC4KXsu6nJ1zrl2lR+uZMHxAssNwKSCeI4q5\nQJGZFZtZI7AauLrVOnnACwBmtguYIimHyCNSN5pZvZk1AS8SeW42ZvZM0AbwKjCh29k45+JiZhyo\namD8MC8UrmPxFIrxQEnUcmnQFm0rQQGQNBeYTOQP/3ZgvqSRkrKBq4CJbezjS8AfopanBqedXpTU\n5qNUnXNdV1nXyPGTp7xQuLjE9czsONwDrJBUALwJvAE0m9lOSd8FngHqgAKgOXqgpH8CmoBHgqZy\nYJKZVUqaDfxW0kwzq2k1bimwFCAnJ4f8/PwEpZJ4tbW1vTq+eIUhjzDkAN3Po7g68p/h0bJ3yM/f\nl6CoOsc/ixRiZjFfwKXAuqjl5cDyGOsL2AsMaaPvO8BXopa/AGwAsmNsLx+YEyvG2bNnW2+2fv36\nZIeQEGHIIww5mHU/j7XbDtjkO562wrLqxATUBf5ZJB/wmnVQA8wsrlNPm4FcSVMlZQKLgTXRK0ga\nFvQB3Ay8ZMERgKQxwc9JRE5PPRosLwK+CXzCzOqjtjU6uICOpGlALlAcR5zOuTiVVTUA+KknF5cO\nTz2ZWZOk24B1QDrwoJkVSloW9K8kctH6IUkGFAJLojbxpKSRwEngVjOrCtp/BGQBzwYPdX/VInc4\nLQD+VdJJ4BSwzMyOJCBX51yg9GgDg7IyGDIgUWefXZjF9VtiZmuBta3aVka93wBMb2dsmxejzeyc\ndtqfBJ6MJy7nXNe03PEU/E+aczH5zGzn+qCyqgbGDeuf7DBcivBC4VwfVFbVwHifbOfi5IXCuT6m\n7kQTVfUnGT/MvwzQxccLhXN9zIGWO578iMLFyQuFc31M6elbY/0ahYuPFwrn+piyoy2Fwk89ufh4\noXCujymraqBfuhgzOCvZobgU4YXCuT7mQFUDZw3tT1qaz6Fw8fFC4VwfU3bUv17cdY4XCuf6mLKq\nBr8+4TrFC4VzfcjJ5lNU1Bz3O55cp3ihcK4POVh9nFPmcyhc53ihcK4P+cvXi/upJxc/LxTO9SGn\n51D4EYXrBC8UzvUhLV/fMXaoX6Nw8fNC4VwfUl5znJEDM+nfLz3ZobgUElehkLRI0m5JRZLubKN/\nuKSnJG2TtEnSrKi+2yVtl1Qo6etR7SMkPSvp7eDn8Ki+5cG+dku6ortJOuciyoPJds51RoeFInh+\n9X3AlUAecL2kvFarfQsoMLMLgM8DK4Kxs4BbgLnAhcDHJbU82e5O4HkzywWeD5YJtr0YmAksAu5v\neYa2c657yquPM3aoX59wnRPPEcVcoMjMis2sEVgNXN1qnTzgBQAz2wVMkZRD5FnaG82s3syagBeB\na4IxVwMPBe8fAj4Z1b7azE6Y2R6gKIjBOddNkULhRxSuc+IpFOOBkqjl0qAt2laCAiBpLjAZmABs\nB+ZLGikpG7gKmBiMyTGz8uD9QSCnE/tzznVSfWMT1Q0nGeuT7VwnZSRoO/cAKyQVAG8CbwDNZrZT\n0neBZ4A6oABobj3YzEySdWaHkpYCSwFycnLIz8/vXgY9qLa2tlfHF68w5BGGHKBreZTXngKg6sBe\n8vNLeyCqzunLn0WqiadQlPGXowCIHCmURa9gZjXAFwEkCdgDFAd9DwAPBH3fIXKEAFAhaayZlUsa\nCxyKd3/BdlcBqwDmzJljCxcujCOV5MjPz6c3xxevMOQRhhyga3n8uehd+NNGPnzJxVwybWTPBNYJ\nffmzSDXxnHraDORKmiopk8iF5jXRK0gaFvQB3Ay8FBQPJI0Jfk4icnrq0WC9NcBNwfubgN9FtS+W\nlCVpKpALbOpKcs65vyivPg74HArXeR0eUZhZk6TbgHVAOvCgmRVKWhb0ryRy0fqh4PRRIbAkahNP\nShoJnARuNbOqoP0e4HFJS4B9wHXB9golPQ7sAJqCMe87XeWc65zyYLJdzhAvFK5z4rpGYWZrgbWt\n2lZGvd8ATG9n7Px22iuBy9vpuxu4O57YnHPx8cl2rqt8ZrZzfYRPtnNd5YXCuT7CJ9u5rvJC4Vwf\n4ZPtXFd5oXCuD/DJdq47vFA41wf4rbGuO7xQONcHHDxdKPwahes8LxTO9QH+wCLXHV4onOsDWo4o\nfLKd6wovFM71AQeqfbKd6zovFM71AQerG/yOJ9dlXiic6wPKq49z1hC/kO26xguFc31AefVxxvkR\nhesiLxTOhVzLZDv/nifXVV4onAu5lsl243wOhesiLxTOhVzLrbF+ROG6yguFcyHXMtnOjyhcV8VV\nKCQtkrRbUpGkO9voHy7pKUnbJG2SNCuq7xuSCiVtl/SYpP5B+68kFQSvvZIKgvYpkhqi+la23p9z\nLn6nJ9sNzUpyJC5VdfiEO0npwH3AR4FSYLOkNWa2I2q1bwEFZvYpSTOC9S+XNB74GpBnZg3BI04X\nAz83s89G7eNeoDpqe++Y2UXdTc45F5lsN2pQJlkZPtnOdU08RxRzgSIzKzazRmA1cHWrdfKAFwDM\nbBcwRVJO0JcBDJCUAWQDB6IHShKR52U/1uUsnHPtOljtT7Zz3RNPoRgPlEQtlwZt0bYC1wBImgtM\nBiaYWRnwA2A/UA5Um9kzrcbOByrM7O2otqnBaacXJbX5zG3nXHx8sp3rrg5PPcXpHmBFcJ3hTeAN\noFnScCJHH1OBKuDXkm4ws4ejxl7Pe48myoFJZlYpaTbwW0kzzawmeoeSlgJLAXJycsjPz09QKolX\nW1vbq+OLVxjyCEMO0Lk8SirrGJ/Z0Ovy7oufRcoys5gv4FJgXdTycmB5jPUF7AWGANcCD0T1fR64\nP2o5A6ggcvTR3vbygTmxYpw9e7b1ZuvXr092CAkRhjzCkINZ/HnUnThpk+942u5b/3bPBtQFfe2z\n6I2A16yDGmBmcZ162gzkSpoqKZPIxeg10StIGhb0AdwMvGSRI4D9wCWSsoNrEZcDO6OGfgTYZWal\nUdsaHVxAR9I0IBcojiNO51wrPtnOJUKHp57MrEnSbcA6IB140MwKJS0L+lcC5wEPSTKgEFgS9G2U\n9ATwOtBE5JTUqqjNL+b9F7EXAP8q6SRwClhmZke6kaNzfVZ5lU+2c90X1zUKM1sLrG3VtjLq/QZg\nejtj7wLuaqfvC220PQk8GU9czrnYyqt9sp3rPp+Z7VyI+WQ7lwheKJwLMZ9s5xLBC4VzIeaT7Vwi\neKFwLsR8sp1LBC8UzoXYwZrjjPUjCtdNXiicC6njJ5upqvcn27nu80LhXEidfmDREC8Urnu8UDgX\nUi2zsv3Uk+suLxTOhVRFTcscCi8Urnu8UDgXUuV+6skliBcK50KqouY4Q/pnMDArUU8TcH2VFwrn\nQqrcJ9u5BPFC4VxIHaw5wVn+ZYAuAbxQOBdS+yvrGD/MC4XrPi8UzoXQ0bpGjtaf5OzRA5MdiguB\nuAqFpEWSdksqknRnG/3DJT0laZukTZJmRfV9Q1KhpO2SHpPUP2j/P5LKJBUEr6uixiwP9rVb0hWJ\nSNS5vqT43VoApo7yQuG6r8NCETyW9D7gSiAPuF5SXqvVvgUUmNkFRJ6LvSIYOx74GpFnXs8i8oS8\nxVHj/sPMLgpea4MxecE6M4FFwP0tj0Z1zsWn+HAd4IXCJUY8RxRzgSIzKzazRmA1cHWrdfKAFwDM\nbBcwRVJO0JcBDJCUAWQDBzrY39XAajM7YWZ7gKIgBudcnEqO1CPBhOHZyQ7FhUA8hWI8UBK1XBq0\nRdsKXAMgaS4wGZhgZmXAD4D9QDlQbWbPRI37anC66kFJwzuxP+dcDCVHGxg7pD+ZGX4Z0nVfombi\n3AOskFQAvAm8ATQHf/yvBqYCVcCvJd1gZg8DPwa+DVjw817gS/HuUNJSYClATk4O+fn5CUol8Wpr\na3t1fPEKQx5hyAE6zmP7ngYGp9Grc+0rn0UYxFMoyoCJUcsTgrbTzKwG+CKAJAF7gGLgCmCPmR0O\n+n4DXAY8bGYVLeMl/RR4Ot79BftcBawCmDNnji1cuDCOVJIjPz+f3hxfvMKQRxhygI7zuOOV55if\nO5qFCy88c0F1Ul/5LMIgnuPSzUCupKmSMolcaF4TvYKkYUEfwM3AS0Hx2A9cIik7KCCXAzuDMWOj\nNvEpYHvwfg2wWFKWpKlALrCpa+k51/ccP9lMRc0JJvr1CZcgHR5RmFmTpNuAdUTuWnrQzAolLQv6\nVwLnAQ9JMqAQWBL0bZT0BPA60ETklNSqYNPfk3QRkVNPe4EvB2MKJT0O7AjG3GpmzQnK17nQK6tq\nAGDiCJ9s5xIjrmsUwa2ra1u1rYx6vwGY3s7Yu4C72mi/Mcb+7gbujic259x77T9SD8DEEX5E4RLD\nb4lwLmRKWwqFn3pyCeKFwrmQKa1qIDM9jTGDs5IdigsJLxTOhUzZ0QbGDutPWpqSHYoLCS8UzoVM\nWVWDf2usSygvFM6FzAEvFC7BvFA4FyKNTac4dOwE44d7oXCJ44XCuRApr27ADMb5EYVLIC8UzoVI\n2dHIZLsJXihcAnmhcC5EWmZl+6knl0heKJwLkdKjDUhw1tD+yQ7FhYgXCudCpOhQLZNGZJOV4Q+F\ndInjhcK5ENlZXsOMswYnOwwXMl4onAuJ+sYm9lTWMeOsIckOxYWMFwrnQmJbaTVmcNHEYckOxYWM\nFwrnQmJJXZ0UAAALZ0lEQVTLvqMAXDzJC4VLLC8UzoXE9rJqpozMZlh2ZscrO9cJcRUKSYsk7ZZU\nJOnONvqHS3pK0jZJmyTNiur7hqRCSdslPSapf9D+fUm7gjFPSRoWtE+R1CCpIHitbL0/59z7vVVx\njOk5fiHbJV6HhUJSOnAfcCWQB1wvKa/Vat8CCszsAuDzwIpg7Hjga8AcM5tF5FGqi4MxzwKzgjFv\nAcujtveOmV0UvJZ1OTvn+ojjJ5vZW1nPuX7Hk+sB8RxRzAWKzKzYzBqB1cDVrdbJA14AMLNdwBRJ\nOUFfBjBAUgaQDRwI1nvGzJqCdV4FJnQrE+f6sOLDdTSfMj+icD0inkIxHiiJWi4N2qJtBa4BkDQX\nmAxMMLMy4AfAfqAcqDazZ9rYx5eAP0QtTw1OO70oaX5cmTjXh7196BiAH1G4HpGRoO3cA6yQVAC8\nCbwBNEsaTuToYypQBfxa0g1m9nDLQEn/BDQBjwRN5cAkM6uUNBv4raSZZlYTvUNJS4GlADk5OeTn\n5ycolcSrra3t1fHFKwx5hCEHeH8ez77VSLpgf+FrHNiZGk+2C+tnEUpmFvMFXAqsi1peDiyPsb6A\nvcAQ4Frggai+zwP3Ry1/AdgAZMfYXj6Raxztxjh79mzrzdavX5/sEBIiDHmEIQez9+ex5Oeb7KM/\nzE9OMF0U1s8ilQCvWQc1wMziOvW0GciVNFVSJpGL0WuiV5A0LOgDuBl4ySJHAPuBSyRlSxJwObAz\nGLMI+CbwCTOrj9rW6OACOpKmAblAcRxxOtdnvVVR69cnXI/psFBY5ILzbcA6In/kHzezQknLJLXc\nkXQesF3SbiJ3R90ejN0IPAG8TuSUVBqwKhjzI2Aw8Gyr22AXANuC01hPAMvM7Ej3U3UunOobm9h/\npJ5zvVC4HhLXNQozWwusbdW2Mur9BmB6O2PvAu5qo/2cdtZ/Engynricc/B2RS0AuV4oXA/xmdnO\npbi3KvyOJ9ezvFA4l+LeqjhGVkYak0ZkJzsUF1JeKJxLcbsOHuOcMYNIT0uN22Jd6vFC4VwKMzMK\nD9Qwc5w/g8L1HC8UzqWw0qMNHKlr5PwJ/tXirud4oXAuhb2+P/IMiou8ULge5IXCuRT28tvvMnRA\nP/L81JPrQV4onEtRjU2nWL/rEPNzR/mFbNejvFA4l6Ke31lBZV0jn57t39DvepYXCudS1OrNJYwd\n2p8FuaOTHYoLOS8UzqWgI8dP8dLbh/nM7Al+2sn1OC8UzqWgP5c1YQbXzp6Y7FBcH+CFwrkUY2a8\nXNbEvKkjmDTSv7bD9TwvFM6lmE17jnCo3rhujh9NuDPDC4VzKebXW0rpnw5Xnn9WskNxfYQXCudS\nyLHjJ1n7Zjlzx2aQnZmoR947F1tchULSIkm7JRVJurON/uGSnpK0TdImSbOi+r4hqVDSdkmPSeof\ntI+Q9Kykt4Ofw6PGLA/2tVvSFYlI1LkweGJLKfWNzSyc6EXCnTkd/rYFz6++D/goUApslrTGzHZE\nrfYtoMDMPiVpRrD+5ZLGA18D8sysQdLjRJ65/XPgTuB5M7snKD53AndIygvWmQmMA56TNN3MmhOU\n82kVNcd5Yktpojf7PnuKGym0oh7fT08LQx6JyGF4dibXzplAv/Qze0B+6pTxiw37uHjSMKYNPXlG\n9+36tnj+t2QuUGRmxQCSVgNXA9GFIg+4B8DMdkmaIiknah8DJJ0EsoEDQfvVwMLg/UNAPnBH0L7a\nzE4AeyQVBTFs6EqCsRysPs731+1O9Gbb9vYZ2k9PC0MeCcjhFxv2MnFENiMHZjIi6nXW0P5MHJ7N\n2KH9yUhwIfltQRl73q3jv66/GI6+ldBtOxdLPIViPFAStVwKzGu1zlbgGuBlSXOBycAEM9si6QfA\nfqABeMbMngnG5JhZefD+INBSWMYDr7ba3/jWQUlaCiwFyMnJIT8/P45U3uuUGT/9WM/fXlhbW8eg\nQQN7fD89LQx5JCKH/JImXq+oY+f+Wo6dhNpGo9neu06aYER/MXqAGDcojYmD05gwOI2Jg9LIyuj8\nBLmSY6f4940NTBuaxsAju6mtq+vS73xvUltbm/I5QHjyiCVRJzrvAVZIKgDeBN4AmoPrDlcDU4Eq\n4NeSbjCzh6MHm5lJstYbjcXMVgGrAObMmWMLFy7sfhY9JD8/n94cX7zCkEcicvhoq2Uzo6ahicq6\nE5RXH6fkSD0lR+spPdrAvsp6Xj14jOf3NwIgwZSRA8kbN4SZ44Ywc9xQZo4bwqhBWW3uq6q+kafe\nKOOHW95iSHZ/fr70UiaNzPbPohcJSx6xxFMoyoDoG7YnBG2nmVkN8EUASQL2AMXAFcAeMzsc9P0G\nuAx4GKiQNNbMyiWNBQ7Fuz/nehNJDM3ux9DsfkwbPeh9/adOGWVVDewsr2Fn+TF2lFeztaSK328r\nP73OWUP6M2lkNsMG9CMzI426E03sq6xnb2UdpwwumTaCe6+7iPHDBpzJ1JwD4isUm4FcSVOJ/MFe\nDHwuegVJw4B6M2sEbgZeMrMaSfuBSyRlEzn1dDnwWjBsDXATkaORm4DfRbU/KumHRC5m5wKbup6i\nc8mVliYmjshm4ohsPjbzL3MfqutPUlhezY4DNRQeqOFAVeQIpLH5FAOz0pmeM5j/feE4PjxjDBdO\n9AcTueTpsFCYWZOk24B1QDrwoJkVSloW9K8EzgMeCk4fFQJLgr6Nkp4AXgeaiJySWhVs+h7gcUlL\ngH3AdcGYwuDuqB3BmFt74o4n55JtaHY/Ljt7FJedPSrZoTgXU1zXKMxsLbC2VdvKqPcbgOntjL0L\nuKuN9koiRxhtjbkbuDue2JxzzvUsn5ntnHMuJi8UzjnnYvJC4ZxzLiYvFM4552LyQuGccy4mLxTO\nOedi8kLhnHMuJpl16iuWeiVJh4lM2uutRgHvJjuIBAhDHmHIAcKRRxhygNTOY7KZje5opVAUit5O\n0mtmNifZcXRXGPIIQw4QjjzCkAOEJ49Y/NSTc865mLxQOOeci8kLxZmxquNVUkIY8ghDDhCOPMKQ\nA4Qnj3b5NQrnnHMx+RGFc865mLxQOOeci8kLhXPOuZi8UCSZpPmSVkr6maRXkh1PV0laKOnlIJeF\nyY6nKySdF8T/hKS/S3Y8XSVpmqQHgqdLpoxUjbu1sPweRfNC0Q2SHpR0SNL2Vu2LJO2WVCTpzljb\nMLOXzWwZ8DTwUE/G255E5AEYUAv0B0p7Ktb2JOiz2Bl8FtcBf9WT8bYnQXkUm9mSno00Pp3JpzfF\n3Von80j671HCmZm/uvgCFgAfALZHtaUD7wDTgExgK5AHnE+kGES/xkSNexwYnKp5AGnBuBzgkVTM\nIRjzCeAPwOdS9bOIGvdEMnLoaj69Ke7u5pHs36NEv+J6ZrZrm5m9JGlKq+a5QJGZFQNIWg1cbWb/\nDny8re1ImgRUm9mxHgy3XYnKI3AUyOqJOGNJVA5mtgZYI+n3wKM9F3HbEvxZJF1n8gF2nNno4tfZ\nPJL9e5Rofuop8cYDJVHLpUFbLEuA/+6xiLqmU3lIukbST4BfAj/q4dji1dkcFkr6zyCPtT0dXCd0\nNo+RklYCF0ta3tPBdUGb+aRA3K21l0dv/T3qMj+i6AXM7K5kx9BdZvYb4DfJjqM7zCwfyE9yGN1m\nZpXAsmTH0VmpGndrYfk9iuZHFIlXBkyMWp4QtKWaMOQRhhwgPHm0CEs+YcmjQ14oEm8zkCtpqqRM\nYDGwJskxdUUY8ghDDhCePFqEJZ+w5NEhLxTdIOkxYANwrqRSSUvMrAm4DVgH7AQeN7PCZMbZkTDk\nEYYcIDx5tAhLPmHJo6v8SwGdc87F5EcUzjnnYvJC4ZxzLiYvFM4552LyQuGccy4mLxTOOedi8kLh\nnHMuJi8UzjnnYvJC4ZxzLiYvFM4552L6//Io03NZj2KVAAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(teh_nl_tk.k, np.abs(teh_nl_ntk.nonlinear_power/teh_nl_tk.nonlinear_power -1))\n", + "#plt.plot(teh_nl_ntk.k, teh_nl_ntk.nonlinear_power)\n", + "#plt.plot(teh_nl_ntk.k, teh_nl_ntk.power)\n", + "\n", + "plt.xscale('log')\n", + "#plt.yscale('log')\n", + "plt.grid(True)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "ExecuteTime": { + "end_time": "2017-05-22T05:15:14.866949", + "start_time": "2017-05-22T05:15:14.857118Z" + }, + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "6.9077552789821368" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.log(1e3)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "ExecuteTime": { + "end_time": "2017-05-22T05:53:05.753761", + "start_time": "2017-05-22T05:53:05.751971Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "from astropy.cosmology import Planck15" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "ExecuteTime": { + "end_time": "2017-05-22T05:53:42.297282", + "start_time": "2017-05-22T05:53:42.294272Z" + }, + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "0.99009841609741556" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Planck15.Om(8.0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python [conda env:hmf3]", + "language": "python", + "name": "conda-env-hmf3-py" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.3" + }, + "latex_envs": { + "bibliofile": "biblio.bib", + "cite_by": "apalike", + "current_citInitial": 1, + "eqLabelWithNumbers": true, + "eqNumInitial": 0 + }, + "nav_menu": {}, + "toc": { + "navigate_menu": true, + "number_sections": true, + "sideBar": true, + "threshold": 6, + "toc_cell": false, + "toc_section_display": "block", + "toc_window_display": true + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/development/hmf_z.py b/development/hmf_z.py new file mode 100644 index 0000000..0a59ace --- /dev/null +++ b/development/hmf_z.py @@ -0,0 +1,113 @@ +""" +Info script for Steven. +Results show issues at higher redshift. +""" + +import hmf +from astropy.cosmology import LambdaCDM +from astropy.io import ascii +import matplotlib.pyplot as plt +from matplotlib.lines import Line2D +import numpy as np +import subprocess +import glob +import os + + +# --- Set parameter values. --- + +logMmin = 12.5 +logMmax = 16 +z = 0.01 + +hmf_model = "Tinker08" +delta_wrt = "crit" +delta_h = 500.0 + +H0 = 67.11 +Om0 = 0.3175 +Ob0 = 0.049 +Ode0 = 0.6825 + +sigma_8 = 0.8344 +n = 0.9624 +transfer_model = "EH" + +# ----------------------------- + + +# --- Create cosmology and set its parameter values. --- + +cosmo_model = LambdaCDM(H0=H0, Om0=Om0, Ob0=Ob0, Ode0=Ode0) +print "" +print "Cosmological model and parameter values:" +print cosmo_model +print "sigma_8 =", sigma_8, ", n =", n, ", transfer_model =", transfer_model + +# ----------------------------------------------------- + + +# --- Create mass function. --- + +mf = hmf.MassFunction(Mmin=logMmin, Mmax=logMmax, z=z, + hmf_model=hmf_model, delta_wrt=delta_wrt, delta_h=delta_h, + cosmo_model=cosmo_model, + sigma_8=sigma_8, n=n, transfer_model=transfer_model) + +print "" +print "Mass function and parameter values:" +print mf.parameter_values +print "" + +# ----------------------------- + +# Different zs: +#zmin = np.log10(0.01) +#zmax = np.log10(2) + +zmin = 0.01 +zmax = 2 +nz = 200 +outdir = 'Plots/' +if not os.path.exists(outdir): + print "The directory ", outdir, "does not exist. Creating it..." + os.makedirs(outdir) + +filnam = outdir + 'hmf_challenge' +exten_ima = '.jpg' +all_imas = filnam + "*" + exten_ima +all_imas_gl = glob.glob(all_imas) +#args = 'rm -f' + all_imas_gl +#outcmd = subprocess.Popen(args) +# Doesn't work for some reason (also don't want to set shell=True). + +for filerm in all_imas_gl: + os.remove(filerm) + +i = 0 +print "Calculating and plotting mass functions at", nz, "redshifts..." +#for z in np.logspace(zmin,zmax,nz): + +for z in np.linspace(zmin,zmax,nz): + mf.update(z=z) + plt.plot(mf.m,mf.dndm,color="darkblue",alpha=1) + outfil = filnam + str(i).zfill(3) + plt.xlim(1e12, 1e16) + plt.ylim(1e-50, 1e-15) + plt.xscale('log') + plt.yscale('log') + plt.xlabel(r"Mass, $[h^{-1}M_\odot]$") + plt.ylabel(r"$dn/dm$, $[h^{4}{\rm Mpc}^{-3}M_\odot^{-1}]$") + plt.text(5e12, 1e-40, '$z = $'+str(format(z, '6.4f')), fontsize=15) + plt.savefig(outfil + exten_ima) + plt.close() + i = i + 1 + +movie_nam = filnam + '.mp4' +ffmpeg_out = 'ffmpeg.out' +outcmd = subprocess.call(['rm', '-f', movie_nam]) +argus = ['ffmpeg', '-f', 'image2', '-r', '10', '-pattern_type', 'glob', '-i', all_imas, movie_nam] +f_std = open(ffmpeg_out, 'w') +outcmd = subprocess.Popen(argus, stdout=f_std, stderr=subprocess.STDOUT) +f_std.close() + diff --git a/docs/examples/deal_with_cosmology.ipynb b/docs/examples/deal_with_cosmology.ipynb index 4882220..82b3ccf 100644 --- a/docs/examples/deal_with_cosmology.ipynb +++ b/docs/examples/deal_with_cosmology.ipynb @@ -434,7 +434,7 @@ "metadata": { "hide_input": true, "kernelspec": { - "display_name": "Python 2", + "display_name": "Python [default]", "language": "python", "name": "python2" }, @@ -448,7 +448,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", - "version": "2.7.11" + "version": "2.7.13" }, "latex_envs": { "bibliofile": "biblio.bib", @@ -456,6 +456,16 @@ "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 0 + }, + "nav_menu": {}, + "toc": { + "navigate_menu": true, + "number_sections": true, + "sideBar": true, + "threshold": 6, + "toc_cell": false, + "toc_section_display": "block", + "toc_window_display": true } }, "nbformat": 4, diff --git a/docs/index.rst b/docs/index.rst index 465c967..11dc7e0 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -5,6 +5,7 @@ .. include:: ../README.rst + Features -------- * Calculate mass functions and related quantities extremely easily. @@ -26,6 +27,7 @@ Features * Nonlinear power spectra via HALOFIT * Functions for sampling the mass function. * CLI scripts both for producing any quantity included, or fitting any quantity. +* Python 2 and 3 compatible Installation ------------ @@ -35,14 +37,8 @@ installed when installing `hmf`, except for one. Explicitly, the dependencies ar numpy, scipy, scitools, cosmolopy and emcee. You will only need `emcee` if you are going to be using the fitting capabilities -of `hmf`. The final, optional, library is pycamb, which can not be installed -using pip currently. - -Please follow the guidelines on its `readme page `_. -installation instructions. - -.. note :: At present, versions of CAMB post March 2013 are not working with - `pycamb`. Please use earlier versions until further notice. +of `hmf`. The final, optional, library is `camb `_, +for which you should consult the repository for install instructions. Finally the `hmf` package needs to be installed: ``pip install hmf``. If you want to install the latest build (not necessarily stable), grab it `here diff --git a/hmf/__init__.py b/hmf/__init__.py index df56c41..00b4b5d 100644 --- a/hmf/__init__.py +++ b/hmf/__init__.py @@ -1,6 +1,6 @@ -__version__ = "2.0.5" -from hmf import MassFunction -import fitting_functions as fits -from cosmo import Cosmology -from transfer import Transfer -from sample import sample_mf +__version__ = "3.0.0" +from .hmf import MassFunction +from . import fitting_functions as fits +from .cosmo import Cosmology +from .transfer import Transfer +from .sample import sample_mf diff --git a/hmf/_cache.py b/hmf/_cache.py index ed96888..ac341e6 100644 --- a/hmf/_cache.py +++ b/hmf/_cache.py @@ -7,130 +7,126 @@ updated. """ from functools import update_wrapper +from copy import copy -def hidden_loc(obj,name): + +def hidden_loc(obj, name): """ Generate the location of a hidden attribute. Importantly deals with attributes beginning with an underscore. """ - return ("_" + obj.__class__.__name__ + "__"+ name).replace("___", "__") + return ("_" + obj.__class__.__name__ + "__" + name).replace("___", "__") + -def cached_property(*parents): +def cached_quantity(f): """ A robust property caching decorator. - This decorator only works when used with the entire system here.... + This decorator is intended for use with the complementary `parameter` decorator. + It caches the decorated quantity, contingent on the parameters it depends on. + When those parameters are modified, a further call to the quantity will result + in a recalculation. + + Examples + -------- + + >>> class CachedClass(object): + >>> @parameter + >>> def a_param(self,val): + >>> return val + >>> + >>> @cached_quantity + >>> def a_quantity(self): + >>> return 2*self.a_param + >>> + >>> @cached_quantity + >>> def a_child_quantity(self): + >>> return self.a_quantity**3 + + This code will calculate ``a_child_quantity`` on the first call, but return the cached + value on all subsequent calls. If `a_param` is modified, the + calculation of either `a_quantity` and `a_child_quantity` will be re-performed when requested. + """ - Usage:: - class CachedClass(Cache): + name = f.__name__ + def _get_property(self): + # Location of the property to be accessed + prop = hidden_loc(self, name) - @cached_property("parent_parameter") - def amethod(self): - ...calculations... - return final_value + # Locations of indexes [they are set up in the parameter decorator] + recalc = hidden_loc(self, "recalc") + recalc_prpa = hidden_loc(self, "recalc_prop_par") + recalc_prpa_static = hidden_loc(self, "recalc_prop_par_static") + recalc_papr = hidden_loc(self, "recalc_par_prop") - @cached_property("amethod") - def a_child_method(self): #method dependent on amethod - final_value = 3*self.amethod - return final_value + # First, if this property has already been indexed, + # we must ensure all dependent parameters are copied into active indexes, + # otherwise they will be lost to their parents. + if name in getattr(self, recalc): + for pr, v in getattr(self, recalc_prpa).items(): + v.update(getattr(self, recalc_prpa_static)[name]) - This code will calculate ``amethod`` on the first call, but return the cached - value on all subsequent calls. If any parent parameter is modified, the - calculation will be re-performed. - """ + # If this property already in recalc and doesn't need updating, just return. + if not getattr(self, recalc).get(name, True): + return getattr(self, prop) - def cache(f): - name = f.__name__ + # Otherwise, if its in recalc, and needs updating, just update it + elif name in getattr(self, recalc): + value = f(self) + setattr(self, prop, value) - def _get_property(self): - # Location of the property to be accessed - prop = hidden_loc(self,name) + # Ensure it doesn't need to be recalculated again + getattr(self, recalc)[name] = False - # Locations of indexes - recalc = hidden_loc(self,"recalc") - recalc_prpa = hidden_loc(self,"recalc_prop_par") - recalc_papr = hidden_loc(self,"recalc_par_prop") - - # If recalc is constructed, and it needs to be updated, - # or if recalc is NOT constructed, recalculate - if getattr(self, recalc).get(name, True): - value = f(self) - setattr(self, prop, value) - - # If recalc constructed and doesn't need updating, just return. - else: - return getattr(self, prop) - - # Figure out if it is a sub-class method - # To be sub-class, it must already be in the recalc dict, but some - # of the parents (if parameters) won't be in the papr dict, and also - # if some parents are properties, THEIR parents won't be in prpa. - # FIXME: this must be relatively slow, must be a better way. - is_sub = False - not_in_papr = any(p not in getattr(self, recalc_papr) for p in parents) - if not_in_papr and name in getattr(self, recalc): - all_pars = [] - if any(p in getattr(self, recalc_prpa) for p in parents): - for p in parents: - all_pars += list(getattr(self, recalc_prpa).get(p, [])) - if any(p not in getattr(self, recalc_prpa)[name] for p in all_pars): - is_sub = True - - - # At this point, the value has been calculated. - # If this quantity isn't in recalc, we need to construct an entry for it. - if name not in getattr(self, recalc) or is_sub: - final = set() - # For each of its parents, either get all its already known parameters, - # or just add the parent to a list (in this case it's a parameter itself). - for p in parents: - # Hit each parent to make sure it's evaluated - getattr(self, p) - if p in getattr(self, recalc_prpa): - final |= set(getattr(self, recalc_prpa)[p]) - else: - final.add(p) - - # final is a set of pure parameters that affect this quantity - if name in getattr(self, recalc_prpa): - # This happens in the case of inheritance - getattr(self, recalc_prpa)[name] |= final - else: - getattr(self, recalc_prpa)[name] = final + return value - # Go through each parameter and add the current quantity to its - # entry (inverse of prpa). - for e in final: - if e in getattr(self, recalc_papr): - getattr(self, recalc_papr)[e].add(name) - else: - getattr(self, recalc_papr)[e] = set([name]) + # Otherwise, we need to create its index for caching. + supered = name in getattr(self, recalc_prpa) # if name is already there, can only be because the method has been supered. + if not supered: + getattr(self, recalc_prpa)[name] = set() # Empty set to which parameter names will be added - getattr(self, recalc)[name] = False - return value - update_wrapper(_get_property, f) + # Go ahead and calculate the value -- each parameter accessed will add itself to the index. + value = f(self) + setattr(self, prop, value) - def _del_property(self): - # Locations of indexes - recalc = hidden_loc(self,"recalc") - recalc_prpa = hidden_loc(self,"recalc_prop_par") - recalc_papr = hidden_loc(self,"recalc_par_prop") + # Invert the index + for par in getattr(self, recalc_prpa)[name]: + getattr(self, recalc_papr)[par].add(name) - # Delete the property AND its recalc dicts - try: - prop = hidden_loc(self,name) - delattr(self, prop) - del getattr(self, recalc)[name] - del getattr(self, recalc_prpa)[name] - for e in getattr(self, recalc_papr): - if name in getattr(self, recalc_papr)[e]: - getattr(self, recalc_papr)[e].remove(name) - except AttributeError: - pass - return property(_get_property, None, _del_property) - return cache + # Copy index to static dict, and remove the index (so that parameters don't keep on trying to add themselves) + if not supered: # If super, don't want to remove the name just yet. + getattr(self, recalc_prpa_static)[name] = copy(getattr(self, recalc_prpa)[name]) + del getattr(self, recalc_prpa)[name] + + # Add entry to master recalc list + getattr(self, recalc)[name] = False + + return value + + update_wrapper(_get_property, f) + + def _del_property(self): + # Locations of indexes + recalc = hidden_loc(self, "recalc") + recalc_prpa = hidden_loc(self, "recalc_prop_par_static") + recalc_papr = hidden_loc(self, "recalc_par_prop") + + # Delete the property AND its recalc dicts + try: + prop = hidden_loc(self, name) + delattr(self, prop) + del getattr(self, recalc)[name] + del getattr(self, recalc_prpa)[name] + # for e in getattr(self, recalc_papr): + # if name in getattr(self, recalc_papr)[e]: + # getattr(self, recalc_papr)[e].remove(name) + except AttributeError: + pass + + return property(_get_property, None, _del_property) + def obj_eq(ob1, ob2): try: @@ -139,102 +135,120 @@ def obj_eq(ob1, ob2): else: return False except ValueError: - if (ob1 == ob2).all(): + if (ob1 == ob2).all(): return True else: return False -def parameter(f): +def parameter(kind): """ - A simple cached property which acts more like an input value. - - This cached property is intended to be used on values that are passed in - ``__init__``, and can possibly be reset later. It provides the opportunity - for complex setters, and also the ability to update dependent properties - whenever the value is modified. - - Usage:: - @set_property("amethod") - def parameter(self,val): - if isinstance(int,val): - return val - else: - raise ValueError("parameter must be an integer") - - @cached_property() - def amethod(self): - return 3*self.parameter - - Note that the definition of the setter merely returns the value to be set, - it doesn't set it to any particular instance attribute. The decorator - automatically sets ``self.__parameter = val`` and defines the get method - accordingly + A decorator which indicates a parameter of a calculation (i.e. something that must be input by user). + + This decorator is intended for use with the complementary `cached_quantity` decorator. + It provides the mechanisms by which the quantities are re-calculated intelligently. + Parameters should be set by the `__init__` call in any class, so that they are set before any + dependent quantity is accessed. + + Parameters + ---------- + kind : str + Either "param", "option", "model", "switch" or "res". Changes the behaviour of the parameter. + "param", "option", "model" and "res" all behave the same currently, while when a "switch" is modified, + all dependent quantities have their dependencies re-indexed. + + Examples + -------- + + >>> class CachedClass(object): + >>> @parameter + >>> def a_param(self,val): + >>> return val + >>> + >>> @cached_quantity + >>> def a_quantity(self): + >>> return 2*self.a_param + >>> + >>> @cached_quantity + >>> def a_child_quantity(self): + >>> return self.a_quantity**3 + + This code will calculate ``a_child_quantity`` on the first call, but return the cached + value on all subsequent calls. If `a_param` is modified, the + calculation of either `a_quantity` and `a_child_quantity` will be re-performed when requested. """ - name = f.__name__ - par_ext = "__parameters" - - def _set_property(self, val): - prop = hidden_loc(self, name) + def param(f): + name = f.__name__ - # The following does any complex setting that is written into the code - val = f(self, val) + def _set_property(self, val): + prop = hidden_loc(self, name) - # Locations of indexes - recalc = hidden_loc(self, "recalc") - recalc_prpa = hidden_loc(self, "recalc_prop_par") - recalc_papr = hidden_loc(self, "recalc_par_prop") - parameters = hidden_loc(self,"parameters") + # The following does any complex setting that is written into the code + val = f(self, val) - try: - # If the property has already been set, we can grab its old value -# old_val = getattr(self, prop) - old_val = self._get_property() - doset = False - except AttributeError: - # Otherwise, it has no old value. - old_val = None - doset = True + # Locations of indexes + recalc = hidden_loc(self, "recalc") + recalc_prpa = hidden_loc(self, "recalc_prop_par") + recalc_papr = hidden_loc(self, "recalc_par_prop") + recalc_prpa_static = hidden_loc(self, "recalc_prop_par_static") - # It's not been set before, so add it to our list of parameters try: - # Only works if something has been set before - parlist = getattr(self,parameters) - parlist += [name] - setattr(self,parameters,parlist) + # If the property has already been set, we can grab its old value + old_val = _get_property(self) + doset = False except AttributeError: - # If nothing has ever been set, start a list - setattr(self, parameters, [name]) - - # Given that *at least one* parameter must be set before properties are calculated, - # we can define the original empty indexes here. - setattr(self, recalc, {}) - setattr(self, recalc_prpa, {}) - setattr(self, recalc_papr, {}) - - # If either the new value is different from the old, or we never set it before - if not obj_eq(val, old_val) or doset: - # Then if its a dict, we update it - if isinstance(val, dict) and hasattr(self, prop) and val: - getattr(self, prop).update(val) - # Otherwise, just overwrite it. Note if dict is passed empty, it clears the whole dict. - else: - setattr(self, prop, val) - - # Make sure children are updated - for pr in getattr(self, recalc_papr).get(name, []): - getattr(self, recalc)[pr] = True - - update_wrapper(_set_property, f) + # Otherwise, it has no old value. + old_val = None + doset = True + + # It's not been set before, so add it to our list of parameters + try: + # Only works if something has been set before + getattr(self, recalc_papr)[name] = set() + + except AttributeError: + # Given that *at least one* parameter must be set before properties are calculated, + # we can define the original empty indexes here. + setattr(self, recalc, {}) + setattr(self, recalc_prpa, {}) + setattr(self, recalc_prpa_static, {}) + setattr(self, recalc_papr, {name: set()}) + + # If either the new value is different from the old, or we never set it before + if not obj_eq(val, old_val) or doset: + # Then if its a dict, we update it + if isinstance(val, dict) and hasattr(self, prop) and val: + getattr(self, prop).update(val) + # Otherwise, just overwrite it. Note if dict is passed empty, it clears the whole dict. + else: + setattr(self, prop, val) - def _get_property(self): - prop = hidden_loc(self,name) - return getattr(self, prop) + # Make sure children are updated + if kind != "switch" or doset: # Normal parameters just update dependencies + for pr in getattr(self, recalc_papr).get(name): + getattr(self, recalc)[pr] = True + else: # Switches mean that dependencies could depend on new parameters, so need to re-index + for pr in getattr(self, recalc_papr)[name]: + delattr(self, pr) + + update_wrapper(_set_property, f) + + def _get_property(self): + prop = hidden_loc(self, name) + recalc_prpa = hidden_loc(self, "recalc_prop_par") + + # Add parameter to any index that hasn't been finalised + for pr, v in getattr(self, recalc_prpa).items(): + v.add(name) + + return getattr(self, prop) + + # Here we set the documentation + doc = (f.__doc__ or "").strip() + if doc.startswith("\n"): + doc = doc[1:] - # Here we set the documentation - doc = (f.__doc__ or "").strip() - if doc.startswith("\n"): - doc = doc[1:] + return property(_get_property, _set_property, None, "**Parameter**: " + doc) - return property(_get_property, _set_property, None,"**Parameter**: "+doc)# \ No newline at end of file + return param diff --git a/hmf/_framework.py b/hmf/_framework.py index f46cdfb..ff78be7 100644 --- a/hmf/_framework.py +++ b/hmf/_framework.py @@ -84,7 +84,7 @@ def update(self, **kwargs): """ Update parameters of the framework with kwargs. """ - for k, v in kwargs.items(): + for k, v in list(kwargs.items()): if hasattr(self, k): setattr(self, k, v) del kwargs[k] @@ -96,7 +96,7 @@ def update(self, **kwargs): def get_all_parameter_names(cls): "Yield all parameter names in the class." K = cls() - return getattr(K,"_"+K.__class__.__name__+"__parameters") + return getattr(K,"_"+K.__class__.__name__+"__recalc_par_prop") @classmethod def get_all_parameter_defaults(cls,recursive=True): @@ -107,12 +107,12 @@ def get_all_parameter_defaults(cls,recursive=True): out[name] = getattr(K,name) if recursive: - for name,default in out.iteritems(): + for name,default in out.items(): if default == {} and name.endswith("_params"): try: out[name] = getattr(getattr(K,name.replace("_params","_model")),"_defaults") except Exception as e: - print e + print(e) pass return out @@ -121,7 +121,7 @@ def get_all_parameter_defaults(cls,recursive=True): def parameter_values(self): "Dictionary of all parameters and their current values" out = {} - for name in getattr(self,"_"+self.__class__.__name__+"__parameters"): + for name in getattr(self,"_"+self.__class__.__name__+"__recalc_par_prop"): out[name] = getattr(self,name) return out @@ -171,4 +171,4 @@ def parameter_info(cls,names=None): docs += "\n ".join(objdoc) +"\n\n" while "\n\n\n" in docs: docs.replace("\n\n\n","\n\n") - print(docs[:-1]) \ No newline at end of file + print((docs[:-1])) \ No newline at end of file diff --git a/hmf/cosmo.py b/hmf/cosmo.py index abe171b..96b32e3 100644 --- a/hmf/cosmo.py +++ b/hmf/cosmo.py @@ -11,9 +11,9 @@ may be used as inputs. """ -import _cache +from . import _cache from astropy.cosmology import Planck13, FLRW, WMAP5, WMAP7, WMAP9, Planck15 -import _framework +from . import _framework import sys import astropy.units as u @@ -50,7 +50,7 @@ def __init__(self, cosmo_model=Planck15, cosmo_params=None): #=========================================================================== # Parameters #=========================================================================== - @_cache.parameter + @_cache.parameter("model") def cosmo_model(self, val): """ The basis for the cosmology -- see astropy documentation. Can be a custom @@ -58,7 +58,7 @@ def cosmo_model(self, val): :type: instance of `astropy.cosmology.FLRW` subclass """ - if isinstance(val, basestring): + if isinstance(val, str): cosmo = get_cosmo(val) return cosmo @@ -67,7 +67,7 @@ def cosmo_model(self, val): else: return val - @_cache.parameter + @_cache.parameter("param") def cosmo_params(self, val): """ Parameters for the cosmology that deviate from the base cosmology passed. @@ -88,7 +88,7 @@ def cosmo_params(self, val): #=========================================================================== # DERIVED PROPERTIES AND FUNCTIONS #=========================================================================== - @_cache.cached_property("cosmo_params", "cosmo_model") + @_cache.cached_quantity def cosmo(self): """ Cosmographic object (:class:`astropy.cosmology.FLRW` object), with custom @@ -96,7 +96,7 @@ def cosmo(self): """ return self.cosmo_model.clone(**self.cosmo_params) - @_cache.cached_property("cosmo") + @_cache.cached_quantity def mean_density0(self): """ Mean density of universe at z=0, [Msun h^2 / Mpc**3] diff --git a/hmf/filters.py b/hmf/filters.py index 06692f7..7194e9b 100644 --- a/hmf/filters.py +++ b/hmf/filters.py @@ -7,8 +7,8 @@ from scipy.interpolate import InterpolatedUnivariateSpline as _spline import scipy.integrate as intg import collections -import _framework -import _utils +from . import _framework +from . import _utils import warnings class Filter(_framework.Component): diff --git a/hmf/fitting/cli_tools.py b/hmf/fitting/cli_tools.py index 36b252c..476206a 100644 --- a/hmf/fitting/cli_tools.py +++ b/hmf/fitting/cli_tools.py @@ -6,10 +6,12 @@ import sys import os -from ConfigParser import SafeConfigParser as cfg +from configparser import ConfigParser as cfg +from functools import reduce + cfg.optionxform = str import numpy as np -import fit +from . import fit import json import time import errno @@ -21,27 +23,34 @@ import copy import collections + def secondsToStr(t): return "%d:%02d:%02d.%03d" % \ - reduce(lambda ll,b : divmod(ll[0],b) + ll[1:], - [(t*1000,),1000,60,60]) + reduce(lambda ll, b: divmod(ll[0], b) + ll[1:], + [(t * 1000,), 1000, 60, 60]) + class CLIError(Exception): '''Generic exception to raise and log different fatal errors.''' + def __init__(self, msg): super(CLIError).__init__(type(self)) self.msg = "E: %s" % msg + def __str__(self): return self.msg + def __unicode__(self): return self.msg + def import_class(cl): d = cl.rfind(".") classname = cl[d + 1:len(cl)] m = __import__(cl[0:d], globals(), locals(), [classname]) return getattr(m, classname) + class CLIRunner(object): """ A class which imports and interprets a config file and runs a fit. @@ -64,7 +73,7 @@ def __init__(self, config, prefix="", restart=False, verbose=0): if self.outdir: try: os.makedirs(self.outdir) - except OSError, e: + except OSError as e: if e.errno != errno.EEXIST: raise @@ -86,7 +95,7 @@ def read_config(self, fname): config.read(fname) # Convert config to a dict - res = {s:dict(config.items(s)) for s in config.sections()} + res = {s: dict(config.items(s)) for s in config.sections()} if "outdir" not in res["IO"]: res["IO"]["outdir"] = "" if "covar_data" not in res["cosmo_paramsParams"]: @@ -115,14 +124,14 @@ def read_config(self, fname): self.burnin = json.loads(res["MCMC"].pop("burnin")) # Downhill-specific - self.downhill_kwargs = {k:json.loads(v) for k,v in res['Downhill'].iteritems()} + self.downhill_kwargs = {k: json.loads(v) for k, v in list(res['Downhill'].items())} del res["Downhill"] - #IO-specific + # IO-specific self.outdir = res["IO"].pop("outdir", None) self.chunks = int(res["IO"].pop("chunks")) - #Data-specific + # Data-specific self.data_file = res["Data"].pop("data_file") self.cov_file = res["Data"].pop("cov_file", None) @@ -135,9 +144,9 @@ def read_config(self, fname): except: pass - self.constraints = {k:json.loads(v) for k, v in res["Constraints"].iteritems()} + self.constraints = {k: json.loads(v) for k, v in list(res["Constraints"].items())} - param_dict = {k:res.pop(k) for k in res.keys() if k.endswith("Params")} + param_dict = {k: res.pop(k) for k in list(res.keys()) if k.endswith("Params")} return param_dict def get_data(self): @@ -218,13 +227,13 @@ def param_setup(self, params): # Deal specifically with cosmology priors, separating types original_cpars = copy.copy(params['cosmo_paramsParams']) - cosmo_priors = {k:json.loads(v) for k, v in params["cosmo_paramsParams"].iteritems()} + cosmo_priors = {k: json.loads(v) for k, v in list(params["cosmo_paramsParams"].items())} # the following rely on covdata - cov_vars = {k:v for k, v in cosmo_priors.iteritems() if v[0] == "cov"} - norm_vars = {k:v for k, v in cosmo_priors.iteritems() if (v[0] == "norm" and len(v) == 2)} + cov_vars = {k: v for k, v in list(cosmo_priors.items()) if v[0] == "cov"} + norm_vars = {k: v for k, v in list(cosmo_priors.items()) if (v[0] == "norm" and len(v) == 2)} # remove these to be left with normal stuff - for k in cov_vars.keys() + norm_vars.keys(): + for k in list(cov_vars.keys()) + list(norm_vars.keys()): del params["cosmo_paramsParams"][k] # sigma_8 and n are special cosmology parameters that don't nest @@ -241,22 +250,20 @@ def param_setup(self, params): if norm_vars: priors += cosmo_cov.get_normal_priors(*norm_vars) - # All non-cosmology-covariance-dependent stuff that is top-level otherparams = params["OtherParams"] - for param, val in otherparams.iteritems(): + for param, val in list(otherparams.items()): priors += self.set_prior(param, val) - # All non-cosmology-covariance-dependent stuff that is nested - for k, v in params.iteritems(): + for k, v in list(params.items()): if k != "OtherParams": - for kk, vv in v.iteritems(): + for kk, vv in list(v.items()): priors += self.set_prior(k[:-6] + ":" + kk, vv) # Create list of all the names of parameters (pure name without :) for prior in priors: - if isinstance(prior.name, basestring): + if isinstance(prior.name, str): keys += [prior.name] else: keys += prior.name @@ -266,24 +273,24 @@ def param_setup(self, params): params['cosmo_paramsParams'] = original_cpars guess = self.get_guess(params, keys, priors) - print "KEY NAMES: ", keys - print "INITIAL GUESSES: ", guess + print(("KEY NAMES: ", keys)) + print(("INITIAL GUESSES: ", guess)) return priors, keys, guess def get_guess(self, params, keys, priors): # Get all parameters to be set as a flat dictionary allparams = {} - for pset, vset in params.iteritems(): - for p, val in vset.iteritems(): + for pset, vset in list(params.items()): + for p, val in list(vset.items()): allparams[p] = val # Get the guesses guess = [] - for i,p in enumerate(priors): - if isinstance(p.name,list): + for i, p in enumerate(priors): + if isinstance(p.name, list): # get raw input - for j,k in enumerate(p.name): + for j, k in enumerate(p.name): val = json.loads(allparams[k.split(":")[-1]]) if val[-1] is None or len(val) == 1: guess.append(p.guess(k)) @@ -322,11 +329,12 @@ def _get_previous_sampler(self): Tries to find a pickled sampler in the current directory to use. """ try: - with open(self.prefix+"sampler.pickle") as f: - h = pickle.load(f) - #check that it lines up with current Parameters + with open(self.prefix + "sampler.pickle") as f: + h = pickle.load(f) + # check that it lines up with current Parameters if h.k != self.nwalkers: - warnings.warn("Imported previous chain had different number of walkers (%s) than specified (%s)"%(h.k,self.nwalkers)) + warnings.warn("Imported previous chain had different number of walkers (%s) than specified (%s)" % ( + h.k, self.nwalkers)) return None else: ## WE DO THE FOLLOWING IN FIT.PY, BUT HAVE TO DO IT HERE TO update @@ -378,7 +386,7 @@ def _setup_instance(self): return instance def run(self): - if self.fit_type=="opt": + if self.fit_type == "opt": self.run_downhill() else: self.run_mcmc() @@ -395,8 +403,8 @@ def run_downhill(self, instance=None): if instance is None: instance = self._setup_instance() - result = fitter.fit(instance,**self.downhill_kwargs) - print "Optimization Result: ", result + result = fitter.fit(instance, **self.downhill_kwargs) + print(("Optimization Result: ", result)) self._write_opt_log(result) return result @@ -410,7 +418,7 @@ def run_mcmc(self): prev_samples = self.sampler.iterations else: instance = self._setup_instance() - if self.fit_type=="both": + if self.fit_type == "both": optres = self.run_downhill(instance) if optres.success: self.guess = list(optres.x) @@ -425,31 +433,32 @@ def run_mcmc(self): start = time.time() if self.chunks == 0: - self.chunks = self.nsamples-prev_samples - nchunks = np.ceil((self.nsamples-prev_samples)/float(self.chunks)) - for i,s in enumerate(fitter.fit(self.sampler,instance, self.nwalkers, - self.nsamples-prev_samples,self.burnin,self.nthreads, - self.chunks)): + self.chunks = self.nsamples - prev_samples + nchunks = np.ceil((self.nsamples - prev_samples) / float(self.chunks)) + for i, s in enumerate(fitter.fit(self.sampler, instance, self.nwalkers, + self.nsamples - prev_samples, self.burnin, self.nthreads, + self.chunks)): # Write out files self.write_iter_pickle(s) - print "Done {0}%. Time per sample: {1}".format(100 * float(i + 1) / nchunks,(time.time() - start) / ((i + 1) * self.chunks*self.nwalkers)) + print(("Done {0}%. Time per sample: {1}".format(100 * float(i + 1) / nchunks, (time.time() - start) / ( + (i + 1) * self.chunks * self.nwalkers)))) total_time = time.time() - start - self._write_log_post(s,total_time) + self._write_log_post(s, total_time) self._write_data(s) - def write_iter_pickle(self,sampler): + def write_iter_pickle(self, sampler): """ Write out a pickle version of the sampler every chunk. """ - with open(self.full_prefix+"sampler.pickle",'w') as f: - pickle.dump(sampler,f) + with open(self.full_prefix + "sampler.pickle", 'w') as f: + pickle.dump(sampler, f) - def _write_opt_log(self,result): - with open(self.full_prefix+"opt.log",'w') as f: - for k,r in zip(self.keys,result.x): - f.write("%s: %s\n"%(k,r)) + def _write_opt_log(self, result): + with open(self.full_prefix + "opt.log", 'w') as f: + for k, r in zip(self.keys, result.x): + f.write("%s: %s\n" % (k, r)) f.write(str(result)) # if hasattr(result,"hess_inv"): # f.write("Covariance Matrix:\n") @@ -459,15 +468,15 @@ def _write_opt_log(self,result): # f.write("Func. Evaluations: %s\n"%result.nfev) # f.write("Message: %s\n"%result.message) - def _write_data(self,sampler): + def _write_data(self, sampler): """ Writes out chains and other data to longer-term readable files (ie ASCII) """ - with open(self.full_prefix+"chain",'w') as f: - np.savetxt(f,sampler.flatchain,header="\t".join(self.keys)) + with open(self.full_prefix + "chain", 'w') as f: + np.savetxt(f, sampler.flatchain, header="\t".join(self.keys)) - with open(self.full_prefix+"likelihoods",'w') as f: - np.savetxt(f,sampler.lnprobability.T) + with open(self.full_prefix + "likelihoods", 'w') as f: + np.savetxt(f, sampler.lnprobability.T) # We can write out any blobs that are parameters if self.blobs: @@ -479,10 +488,10 @@ def _write_data(self,sampler): sh = numblobs.shape numblobs = numblobs.reshape(sh[0] * sh[1], sh[2]) with open(self.full_prefix + "derived_parameters", "w") as f: - np.savetxt(f, numblobs,header="\t".join([self.blobs[ii] for ii in range(self.n_dparams)])) + np.savetxt(f, numblobs, header="\t".join([self.blobs[ii] for ii in range(self.n_dparams)])) def _write_log_pre(self): - with open(self.full_prefix + "log",'w') as f: + with open(self.full_prefix + "log", 'w') as f: f.write("Nsamples: %s\n" % self.nsamples) f.write("Nwalkers: %s\n" % self.nwalkers) f.write("Burnin: %s\n" % self.burnin) @@ -490,9 +499,10 @@ def _write_log_pre(self): def _write_log_post(self, sampler, total_time): with open(self.full_prefix + "log", 'a') as f: - f.write("Total Time: %s\n"%secondsToStr(total_time)) + f.write("Total Time: %s\n" % secondsToStr(total_time)) if isinstance(self.burnin, int): - f.write("Average time: %s\n" % (total_time / (self.nwalkers * self.nsamples + self.nwalkers * self.burnin))) + f.write( + "Average time: %s\n" % (total_time / (self.nwalkers * self.nsamples + self.nwalkers * self.burnin))) else: f.write("Average time (discounting burnin): %s\n" % (total_time / (self.nwalkers * self.nsamples))) f.write("Mean values = %s\n" % np.mean(sampler.chain, axis=0)) diff --git a/hmf/fitting/fit.py b/hmf/fitting/fit.py index a139777..5d94990 100644 --- a/hmf/fitting/fit.py +++ b/hmf/fitting/fit.py @@ -38,7 +38,7 @@ def should_pickle(k): class EnsembleSampler(es): def __getstate__(self): - return dict((k, v) for (k, v) in self.__dict__.iteritems() if should_pickle(k)) + return dict((k, v) for (k, v) in list(self.__dict__.items()) if should_pickle(k)) except ImportError: HAVE_EMCEE = False @@ -68,7 +68,7 @@ def model(parm, h, self): The log likelihood of the model at the given position. """ if self.verbose > 1: - print "Params: ", zip(self.attrs, parm) + print(("Params: ", list(zip(self.attrs, parm)))) ll = 0 p = copy.copy(parm) @@ -102,12 +102,12 @@ def model(parm, h, self): h.update(**param_dict) except ValueError as e: if self.relax: - print "WARNING: PARAMETERS FAILED ON UPDATE, RETURNING INF: ", zip(self.attrs, parm) - print e - print traceback.format_exc() + print(("WARNING: PARAMETERS FAILED ON UPDATE, RETURNING INF: ", list(zip(self.attrs, parm)))) + print(e) + print((traceback.format_exc())) return ret_arg(-np.inf, self.blobs) else: - print traceback.format_exc() + print((traceback.format_exc())) raise e # Get the quantity to compare (if exceptions are raised, treat properly) @@ -115,12 +115,12 @@ def model(parm, h, self): q = getattr(h, self.quantity) except Exception as e: if self.relax: - print "WARNING: PARAMETERS FAILED WHEN CALCULATING QUANTITY, RETURNING INF: ", zip(self.attrs, parm) - print e - print traceback.format_exc() + print(("WARNING: PARAMETERS FAILED WHEN CALCULATING QUANTITY, RETURNING INF: ", list(zip(self.attrs, parm)))) + print(e) + print((traceback.format_exc())) return ret_arg(-np.inf, self.blobs) else: - print traceback.format_exc() + print((traceback.format_exc())) raise e @@ -131,17 +131,17 @@ def model(parm, h, self): ll += np.sum(norm.logpdf(self.data, loc=q, scale=self.sigma)) # Add the likelihood of the contraints - for k, v in self.constraints.iteritems(): + for k, v in list(self.constraints.items()): ll += norm.logpdf(getattr(h, k), loc=v[0], scale=v[1]) if self.verbose > 2: - print "CONSTRAINT: ", k, getattr(h, k) + print(("CONSTRAINT: ", k, getattr(h, k))) if self.verbose: - print "Likelihood: ", ll + print(("Likelihood: ", ll)) if self.verbose > 1 : - print "Update Dictionary: ", param_dict + print(("Update Dictionary: ", param_dict)) if self.verbose > 2: - print "Final Quantity: ", q + print(("Final Quantity: ", q)) # Get blobs to return as well. if self.blobs is not None: @@ -207,7 +207,7 @@ def __init__(self, priors, data, quantity, constraints, sigma, guess=[], blobs=N # Save which attributes are updatable as a list self.attrs = [] for prior in self.priors: - if isinstance(prior.name, basestring): + if isinstance(prior.name, str): self.attrs += [prior.name] else: self.attrs += prior.name @@ -421,7 +421,7 @@ def _run_burnin(self,burnin,initial_pos): if self.verbose > 0: burnin = self.sampler.iterations - print "Used %s samples for burnin" % self.sampler.iterations + print(("Used %s samples for burnin" % self.sampler.iterations)) self.sampler.reset() return initial_pos, lnprob,rstate,blobs0 diff --git a/hmf/fitting_functions.py b/hmf/fitting_functions.py index d9d59a6..4c2df8e 100644 --- a/hmf/fitting_functions.py +++ b/hmf/fitting_functions.py @@ -8,10 +8,10 @@ import numpy as np from scipy.interpolate import InterpolatedUnivariateSpline as _spline import scipy.special as sp -import cosmo as csm -import _framework +from . import cosmo as csm +from . import _framework from copy import copy -import _utils +from . import _utils diff --git a/hmf/functional.py b/hmf/functional.py index 5de0372..7450914 100644 --- a/hmf/functional.py +++ b/hmf/functional.py @@ -11,7 +11,7 @@ """ import collections -import hmf +from . import hmf import itertools @@ -60,7 +60,7 @@ def get_best_param_order(kls, q="dndm", **kwargs): """ a = kls(**kwargs) - if isinstance(q, basestring): + if isinstance(q, str): getattr(a, q) else: for qq in q: @@ -68,7 +68,7 @@ def get_best_param_order(kls, q="dndm", **kwargs): final_list = [] final_num = [] - for k, v in getattr(a,"_"+a.__class__.__name__+"__recalc_par_prop").iteritems(): + for k, v in getattr(a,"_"+a.__class__.__name__+"__recalc_par_prop").items(): num = len(v) for i, l in enumerate(final_num): if l >= num: @@ -166,10 +166,10 @@ def get_hmf(req_qauntities, get_label=True, framework=hmf.MassFunction, >>> print [x[0][0]/1e10 for x in big_list] [8.531878308131338, 68.2550264650507, 230.36071431954613, 546.0402117204056, 1066.4847885164174, 1842.885714556369, 2926.434259689049, 4368.321693763245] """ - if isinstance(req_qauntities, basestring): + if isinstance(req_qauntities, str): req_qauntities = [req_qauntities] lists = {} - for k, v in kwargs.items(): + for k, v in list(kwargs.items()): if isinstance(v, (list, tuple)): if len(v) > 1: lists[k] = kwargs.pop(k) @@ -184,7 +184,7 @@ def get_hmf(req_qauntities, get_label=True, framework=hmf.MassFunction, yield [[getattr(x, a) for a in req_qauntities], x] if len(lists) == 1: - for k, v in lists.iteritems(): + for k, v in lists.items(): for vv in v: x.update(**{k: vv}) if get_label: @@ -205,12 +205,12 @@ def get_hmf(req_qauntities, get_label=True, framework=hmf.MassFunction, pass # # add the rest in any order (there shouldn't actually be any) - for k in lists.items(): + for k in list(lists.items()): if isinstance(lists[k], (list, tuple)): ordered_kwargs[k] = lists.pop(k) ordered_list = [ordered_kwargs[k] for k in ordered_kwargs] - final_list = [dict(zip(ordered_kwargs.keys(), v)) for v in itertools.product(*ordered_list)] + final_list = [collections.OrderedDict(list(zip(list(ordered_kwargs.keys()), v))) for v in itertools.product(*ordered_list)] for vals in final_list: x.update(**vals) @@ -223,11 +223,11 @@ def get_hmf(req_qauntities, get_label=True, framework=hmf.MassFunction, def _make_label(d): label = "" - for key, val in d.iteritems(): - if isinstance(val, basestring): + for key, val in d.items(): + if isinstance(val, str): label += val + ", " elif isinstance(val, dict): - for k, v in val.iteritems(): + for k, v in val.items(): label += "%s: %s, "%(k, v) else: label += "%s: %s, "%(key, val) diff --git a/hmf/growth_factor.py b/hmf/growth_factor.py index 6f8d0fe..f8486de 100644 --- a/hmf/growth_factor.py +++ b/hmf/growth_factor.py @@ -9,9 +9,9 @@ import numpy as np from scipy import integrate as intg -from _framework import Component as Cmpt +from ._framework import Component as Cmpt from scipy.interpolate import InterpolatedUnivariateSpline as _spline -from _utils import inherit_docstrings as _inherit +from ._utils import inherit_docstrings as _inherit class GrowthFactor(Cmpt): @@ -296,7 +296,7 @@ def _d_plus(self, z, getvec=False): Omega_m = om/denom Omega_L = self.cosmo.Ode0/denom coeff = 5.*Omega_m/(2./a) - term1 = Omega_m*(4./7.) + term1 = Omega_m**(4./7.) term3 = (1. + 0.5*Omega_m)*(1. + Omega_L/70.) return coeff/(term1 - Omega_L + term3) diff --git a/hmf/halofit.py b/hmf/halofit.py index abb9ecc..94bc73e 100644 --- a/hmf/halofit.py +++ b/hmf/halofit.py @@ -10,7 +10,7 @@ from scipy.integrate import simps as _simps import copy from scipy.interpolate import InterpolatedUnivariateSpline as _spline -from cosmo import Cosmology as csm +from .cosmo import Cosmology as csm def _get_spec(k, delta_k, sigma_8): """ @@ -60,7 +60,7 @@ def _get_spec(k, delta_k, sigma_8): try: lnsig1 = lnsig_old except: - print "WARNING: LOWEST R NOT LOW ENOUGH IN _GET_SPEC. ln(sig) starts below 0: ", lnsig1 + print("WARNING: LOWEST R NOT LOW ENOUGH IN _GET_SPEC. ln(sig) starts below 0: ", lnsig1) break lnsig_old = copy.copy(lnsig1) @@ -81,7 +81,7 @@ def _get_spec(k, delta_k, sigma_8): try: dev1, dev2 = sig_of_r.derivatives(np.log(1.0 / rknl))[1:3] except Exception as e: - print "HALOFIT WARNING: Requiring extra iterations to find derivatives of sigma at 1/rknl (this often happens at high redshift)." + print("HALOFIT WARNING: Requiring extra iterations to find derivatives of sigma at 1/rknl (this often happens at high redshift).") lnr = np.linspace(np.log(0.2 / rknl), np.log(5 / rknl), 100) lnsig = np.empty(100) diff --git a/hmf/hmf.py b/hmf/hmf.py index 72bc66e..72e5b02 100755 --- a/hmf/hmf.py +++ b/hmf/hmf.py @@ -11,14 +11,14 @@ import numpy as np import copy import logging -import fitting_functions as ff -import transfer -from _cache import parameter, cached_property -from integrate_hmf import hmf_integral_gtm as int_gtm +from . import fitting_functions as ff +from . import transfer +from ._cache import parameter, cached_quantity +from .integrate_hmf import hmf_integral_gtm as int_gtm from numpy import issubclass_ logger = logging.getLogger('hmf') -from filters import TopHat, Filter -from _framework import get_model +from .filters import TopHat, Filter +from ._framework import get_model from scipy.optimize import minimize from scipy.interpolate import InterpolatedUnivariateSpline as spline import warnings @@ -89,7 +89,7 @@ def __init__(self, Mmin=10, Mmax=15, dlog10m=0.01, hmf_model=ff.Tinker08, hmf_pa #=========================================================================== # PARAMETERS #=========================================================================== - @parameter + @parameter("res") def Mmin(self, val): """ Minimum mass at which to perform analysis [units :math:`\log_{10}M_\odot h^{-1}`]. @@ -98,7 +98,7 @@ def Mmin(self, val): """ return val - @parameter + @parameter("res") def Mmax(self, val): """ Maximum mass at which to perform analysis [units :math:`\log_{10}M_\odot h^{-1}`]. @@ -107,7 +107,7 @@ def Mmax(self, val): """ return val - @parameter + @parameter("res") def dlog10m(self, val): """ log10 interval between mass bins @@ -116,18 +116,18 @@ def dlog10m(self, val): """ return val - @parameter + @parameter("model") def filter_model(self, val): """ A model for the window/filter function. :type: str or :class:`hmf.filters.Filter` subclass """ - if not issubclass_(val, Filter) and not isinstance(val, basestring): + if not issubclass_(val, Filter) and not isinstance(val, str): raise ValueError("filter must be a Filter or string, got %s" % type(val)) return val - @parameter + @parameter("param") def filter_params(self, val): """ Model parameters for `filter_model`. @@ -136,7 +136,7 @@ def filter_params(self, val): """ return val - @parameter + @parameter("param") def delta_c(self, val): """ The critical overdensity for collapse, :math:`\delta_c`. @@ -155,18 +155,18 @@ def delta_c(self, val): return val - @parameter + @parameter("model") def hmf_model(self, val): """ A model to use as the fitting function :math:`f(\sigma)` :type: str or `hmf.fitting_functions.FittingFunction` subclass """ - if not issubclass_(val, ff.FittingFunction) and not isinstance(val, basestring): + if not issubclass_(val, ff.FittingFunction) and not isinstance(val, str): raise ValueError("hmf_model must be a ff.FittingFunction or string, got %s" % type(val)) return val - @parameter + @parameter("param") def hmf_params(self, val): """ Model parameters for `hmf_model`. @@ -177,7 +177,7 @@ def hmf_params(self, val): raise ValueError("hmf_params must be a dictionary") return val - @parameter + @parameter("param") def delta_h(self, val): """ The overdensity for the halo definition, with respect to :attr:`delta_wrt` @@ -198,7 +198,7 @@ def delta_h(self, val): - @parameter + @parameter("switch") def delta_wrt(self, val): """ Defines what the overdensity of a halo is with respect to, mean density @@ -251,15 +251,14 @@ def delta_wrt(self, val): #-------------------------------- PROPERTIES ------------------------------ - @cached_property("z", "mean_density0") + @cached_quantity def mean_density(self): """ Mean density of universe at redshift z """ return self.mean_density0 * (1 + self.z) ** 3 - @cached_property("hmf_model", "sigma", "z", "delta_halo", "nu", "m", "hmf_params", - "cosmo", "delta_c") + @cached_quantity def hmf(self): """ Instantiated model for the hmf fitting function. @@ -269,35 +268,35 @@ def hmf(self): delta_halo=self.delta_halo, omegam_z=self.cosmo.Om(self.z), delta_c=self.delta_c, n_eff=self.n_eff, ** self.hmf_params) - elif isinstance(self.hmf_model, basestring): + elif isinstance(self.hmf_model, str): return get_model(self.hmf_model, "hmf.fitting_functions", m=self.m, nu2=self.nu, z=self.z, delta_halo=self.delta_halo, omegam_z=self.cosmo.Om(self.z), delta_c=self.delta_c, n_eff=self.n_eff, ** self.hmf_params) - @cached_property("mean_density0", "filter_model", "delta_c", "k", "_unnormalised_power", "filter_params") + @cached_quantity def filter(self): """ Instantiated model for filter/window functions. """ if issubclass_(self.filter_model, Filter): return self.filter_model(self.k,self._unnormalised_power, **self.filter_params) - elif isinstance(self.filter_model, basestring): + elif isinstance(self.filter_model, str): return get_model(self.filter_model, "hmf.filters", k=self.k, power=self._unnormalised_power, **self.filter_params) - @cached_property("Mmin", "Mmax", "dlog10m") + @cached_quantity def m(self): """Masses""" return 10 ** np.arange(self.Mmin, self.Mmax, self.dlog10m) - @cached_property("m") + @cached_quantity def M(self): "Masses (alias of m, deprecated)" return self.m - @cached_property("delta_wrt", "delta_h", "z", "cosmo") + @cached_quantity def delta_halo(self): """ Overdensity of a halo w.r.t mean density""" if self.delta_wrt == 'mean': @@ -306,28 +305,28 @@ def delta_halo(self): elif self.delta_wrt == 'crit': return self.delta_h / self.cosmo.Om(self.z) - @cached_property("radii", "filter") + @cached_quantity def _unn_sigma0(self): """ Unnormalised mass variance at z=0 """ return self.filter.sigma(self.radii) - @cached_property("_normalisation", "_unn_sigma0") + @cached_quantity def _sigma_0(self): """ The normalised mass variance at z=0 :math:`\sigma` """ return self._normalisation * self._unn_sigma0 - @cached_property("m", "filter","mean_density0") + @cached_quantity def radii(self): """ The radii corresponding to the masses `m`. """ return self.filter.mass_to_radius(self.m,self.mean_density0) - @cached_property("radii", "filter") + @cached_quantity def _dlnsdlnm(self): """ The value of :math:`\left|\frac{\d \ln \sigma}{\d \ln m}\right|`, ``len=len(m)`` @@ -340,21 +339,21 @@ def _dlnsdlnm(self): """ return 0.5 * self.filter.dlnss_dlnm(self.radii) - @cached_property("_sigma_0", "growth_factor") + @cached_quantity def sigma(self): """ The mass variance at `z`, ``len=len(m)`` """ return self._sigma_0 * self.growth_factor - @cached_property("sigma", "delta_c") + @cached_quantity def nu(self): r""" The parameter :math:`\nu = \left(\frac{\delta_c}{\sigma}\right)^2`, ``len=len(m)`` """ return (self.delta_c / self.sigma) ** 2 - @cached_property("nu") + @cached_quantity def mass_nonlinear(self): """ The nonlinear mass, nu(Mstar) = 1. @@ -381,14 +380,14 @@ def mass_nonlinear(self): nu = spline(self.nu,self.m,k=5) return nu(1) - @cached_property("sigma") + @cached_quantity def lnsigma(self): """ Natural log of inverse mass variance, ``len=len(m)`` """ return np.log(1 / self.sigma) - @cached_property("_dlnsdlnm") + @cached_quantity def n_eff(self): """ Effective spectral index at scale of halo radius, ``len=len(m)`` @@ -404,14 +403,14 @@ def n_eff(self): """ return -3.0 * (2.0 * self._dlnsdlnm + 1.0) - @cached_property("hmf") + @cached_quantity def fsigma(self): """ The multiplicity function, :math:`f(\sigma)`, for `hmf_model`. ``len=len(m)`` """ return self.hmf.fsigma - @cached_property("fsigma", "mean_density0", "_dlnsdlnm", "m", "z" ) + @cached_quantity def dndm(self): """ The number density of haloes, ``len=len(m)`` [units :math:`h^4 M_\odot^{-1} Mpc^{-3}`] @@ -450,14 +449,14 @@ def dndm(self): return dndm - @cached_property("m", "dndm") + @cached_quantity def dndlnm(self): """ The differential mass function in terms of natural log of `m`, ``len=len(m)`` [units :math:`h^3 Mpc^{-3}`] """ return self.m * self.dndm - @cached_property("m", "dndm") + @cached_quantity def dndlog10m(self): """ The differential mass function in terms of log of `m`, ``len=len(m)`` [units :math:`h^3 Mpc^{-3}`] @@ -507,7 +506,7 @@ def _gtm(self, dndm, mass_density=False): # Since ngtm may have been extended, we cut it back return ngtm[:size] - @cached_property("m", "dndm") + @cached_quantity def ngtm(self): """ The cumulative mass function above `m`, ``len=len(m)`` [units :math:`h^3 Mpc^{-3}`] @@ -521,7 +520,7 @@ def ngtm(self): """ return self._gtm(self.dndm) - @cached_property("m", "dndm") + @cached_quantity def rho_gtm(self): """ Mass density in haloes `>m`, ``len=len(m)`` [units :math:`M_\odot h^2 Mpc^{-3}`] @@ -536,7 +535,7 @@ def rho_gtm(self): return self._gtm(self.dndm, mass_density=True) - @cached_property("mean_density0", 'rho_gtm') + @cached_quantity def rho_ltm(self): """ Mass density in haloes ` -15 or self.lnk_max < 9: @@ -240,26 +238,26 @@ def _unn_sig8(self): return filt.sigma(8.0)[0] - @cached_property("_unn_sig8", "sigma_8") + @cached_quantity def _normalisation(self): # Calculate the normalization factor return self.sigma_8 / self._unn_sig8 - @cached_property("_normalisation", "_unnormalised_power") + @cached_quantity def _power0(self): """ Normalised power spectrum at z=0 [units :math:`Mpc^3/h^3`] """ return self._normalisation ** 2 * self._unnormalised_power - @cached_property("sigma_8", "_unnormalised_lnT", "lnk", "mean_density0") - def _transfer(self): + @cached_quantity + def transfer_function(self): """ Normalised CDM log transfer function """ return self._normalisation * np.exp(self._unnormalised_lnT) - @cached_property("cosmo", "growth_model", "growth_params") + @cached_quantity def growth(self): "The instantiated growth model" if np.issubclass_(self.growth_model, gf.GrowthFactor): @@ -268,31 +266,28 @@ def growth(self): return get_model(self.growth_model, "hmf.growth_factor", cosmo=self.cosmo, **self.growth_params) - @cached_property("z", "growth") + @cached_quantity def growth_factor(self): r""" The growth factor """ - if self.z > 0: - return self.growth.growth_factor(self.z) - else: - return 1.0 + return self.growth.growth_factor(self.z) - @cached_property("growth_factor", "_power0") + @cached_quantity def power(self): """ Normalised log power spectrum [units :math:`Mpc^3/h^3`] """ return self.growth_factor ** 2 * self._power0 - @cached_property("k", "power") + @cached_quantity def delta_k(self): r""" Dimensionless power spectrum, :math:`\Delta_k = \frac{k^3 P(k)}{2\pi^2}` """ return self.k ** 3 * self.power / (2 * np.pi ** 2) - @cached_property("k", "nonlinear_delta_k") + @cached_quantity def nonlinear_power(self): """ Non-linear log power [units :math:`Mpc^3/h^3`] @@ -301,10 +296,10 @@ def nonlinear_power(self): """ return self.k ** -3 * self.nonlinear_delta_k * (2 * np.pi ** 2) - @cached_property("delta_k", "k", "z", "sigma_8", "cosmo", 'takahashi') + @cached_quantity def nonlinear_delta_k(self): r""" Dimensionless nonlinear power spectrum, :math:`\Delta_k = \frac{k^3 P_{\rm nl}(k)}{2\pi^2}` """ - return _hfit(self.k,self.delta_k,self.sigma_8,self.z,self.cosmo) + return _hfit(self.k,self.delta_k,self.sigma_8,self.z,self.cosmo,self.takahashi) diff --git a/hmf/transfer_models.py b/hmf/transfer_models.py index 7156529..394d4ea 100644 --- a/hmf/transfer_models.py +++ b/hmf/transfer_models.py @@ -6,10 +6,11 @@ ''' import numpy as np from scipy.interpolate import InterpolatedUnivariateSpline as spline -from _framework import Component +from ._framework import Component + try: - import pycamb -except ImportError: + import camb +except ImportError as e: pass @@ -52,9 +53,14 @@ def lnt(self, lnk): """ pass -class CAMB(TransferComponent): + + +class FromFile(TransferComponent): """ - Transfer function computed by CAMB. + Import a transfer function from file. + + .. note:: The file should be in the same format as output from CAMB, + or else in two-column ASCII format (k,T). Parameters ---------- @@ -66,30 +72,10 @@ class CAMB(TransferComponent): parameters are the following. To see their default values, check the :attr:`_defaults` class attribute. - :Scalar_initial_condition: Initial scalar perturbation mode (adiabatic=1, CDM iso=2, - Baryon iso=3,neutrino density iso =4, neutrino velocity iso = 5) - :lAccuracyBoost: Larger to keep more terms in the hierarchy evolution - :AccuracyBoost: Increase accuracy_boost to decrease time steps, use more k - values, etc.Decrease to speed up at cost of worse accuracy. - Suggest 0.8 to 3. - :w_perturb: Whether to perturb the dark energy equation of state. - :transfer__k_per_logint: Number of wavenumbers estimated per log interval by CAMB - Default of 11 gets best performance for requisite accuracy of mass function. - :transfer__kmax: Maximum value of the wavenumber. - Default of 0.25 is high enough for requisite accuracy of mass function. - :ThreadNum: Number of threads to use for calculation of transfer - function by CAMB. Default 0 automatically determines the number. - :scalar_amp: Amplitude of the power spectrum. It is not recommended to modify - this parameter, as the amplitude is typically set by sigma_8. + :fname: str + Location of the file to import. """ - _defaults = {"Scalar_initial_condition":1, - "lAccuracyBoost":1, - "AccuracyBoost":1, - "w_perturb":False, - "transfer__k_per_logint":11, - "transfer__kmax":5, - "ThreadNum":0, - "scalar_amp":1e-9} + _defaults = {"fname":""} def _check_low_k(self, lnk, lnT, lnkmin): """ @@ -132,23 +118,10 @@ def lnt(self, lnk): lnt : array_like The log of the transfer function at lnk. """ - - pycamb_dict = {"w_lam":self.cosmo.w(0), - "TCMB":self.cosmo.Tcmb0.value, - "Num_Nu_massless":self.cosmo.Neff, - "omegab":self.cosmo.Ob0, - "omegac":self.cosmo.Om0 - self.cosmo.Ob0, - "H0":self.cosmo.H0.value, - "omegav":self.cosmo.Ode0, - "omegak":self.cosmo.Ok0, - "omegan":self.cosmo.Onu0, -# "scalar_index":self.n, - } - - cdict = dict(pycamb_dict, - **self.params) - T = pycamb.transfers(**cdict)[1] - T = np.log(T[[0, 6], :, 0]) + try: + T = np.log(np.genfromtxt(self.params["fname"])[:, [0, 6]].T) + except IndexError: + T = np.log(np.genfromtxt(self.params["fname"])[:, [0, 1]].T) if lnk[0] < T[0, 0]: lnkout, lnT = self._check_low_k(T[0, :], T[1, :], lnk[0]) @@ -157,12 +130,10 @@ def lnt(self, lnk): lnT = T[1, :] return spline(lnkout, lnT, k=1)(lnk) -class FromFile(CAMB): - """ - Import a transfer function from file. - .. note:: The file should be in the same format as output from CAMB, - or else in two-column ASCII format (k,T). +class CAMB(FromFile): + """ + Transfer function computed by CAMB. Parameters ---------- @@ -170,14 +141,29 @@ class FromFile(CAMB): The cosmology used in the calculation \*\*model_parameters : unpack-dict - Parameters specific to this model. In this case, available - parameters are the following. To see their default values, - check the :attr:`_defaults` class attribute. + Parameters specific to this model. + + **camb_params:** An instantiated ``CAMBparams`` object, pre-set with desired accuracy options etc. - :fname: str - Location of the file to import. """ - _defaults = {"fname":""} + _defaults = {"camb_params": None} + + def __init__(self, *args, **kwargs): + super(CAMB, self).__init__(*args, **kwargs) + + # Save the CAMB object properly for use + # Set the cosmology + if self.params['camb_params'] is None: + self.params['camb_params'] = camb.CAMBparams() + + self.params['camb_params'].set_cosmology(H0=self.cosmo.H0.value, + ombh2=self.cosmo.Ob0 * self.cosmo.h ** 2, + omch2=(self.cosmo.Om0 - self.cosmo.Ob0) * self.cosmo.h ** 2, + omk=self.cosmo.Ok0, + nnu=self.cosmo.Neff, + standard_neutrino_neff=self.cosmo.Neff, + TCMB=self.cosmo.Tcmb0.value) + self.params['camb_params'].WantTransfer = True def lnt(self, lnk): """ @@ -193,18 +179,20 @@ def lnt(self, lnk): lnt : array_like The log of the transfer function at lnk. """ - try: - T = np.log(np.genfromtxt(self.params["fname"])[:, [0, 6]].T) - except IndexError: - T = np.log(np.genfromtxt(self.params["fname"])[:, [0, 1]].T) + + camb_transfers = camb.get_transfer_functions(self.params['camb_params']) + T = camb_transfers.get_matter_transfer_data().transfer_data + T = np.log(T[[0, 6], :, 0]) if lnk[0] < T[0, 0]: lnkout, lnT = self._check_low_k(T[0, :], T[1, :], lnk[0]) else: lnkout = T[0, :] lnT = T[1, :] + return spline(lnkout, lnT, k=1)(lnk) + class FromArray(FromFile): """ Use a spline over a given array to define the transfer function diff --git a/hmf/wdm.py b/hmf/wdm.py index 5a3e0d3..b1adc7f 100644 --- a/hmf/wdm.py +++ b/hmf/wdm.py @@ -8,10 +8,12 @@ ''' import numpy as np -from transfer import Transfer as _Tr -from hmf import MassFunction as _MF -from _cache import parameter, cached_property -from _framework import Component, get_model +from .transfer import Transfer as _Tr +from .hmf import MassFunction as _MF +from ._cache import parameter, cached_quantity +from ._framework import Component, get_model +from .cosmo import Planck15 + import astropy.units as u #=============================================================================== @@ -42,7 +44,7 @@ class WDM(Component): To see the default values, check the :attr:`_defaults` class attribute. """ - def __init__(self, mx, cosmo,z, **model_params): + def __init__(self, mx, cosmo = Planck15, z =0, **model_params): self.mx = mx self.cosmo = cosmo self.rho_mean = (1+z)**3 * (self.cosmo.Om0 * self.cosmo.critical_density0 / self.cosmo.h ** 2).to(u.MsolMass / u.Mpc ** 3).value * 1e6 @@ -65,7 +67,7 @@ def transfer(self, lnk): The WDM transfer function at `lnk`. """ - pass + raise NotImplementedError("You shouldn't call the WDM class, and any subclass should define the transfer method.") @@ -165,7 +167,7 @@ class WDMRecalibrateMF(Component): To see the default values, check the :attr:`_defaults` class attribute. """ - def __init__(self, m, dndm0, wdm, **model_parameters): + def __init__(self, m, dndm0, wdm=Viel05(mx=1.0), **model_parameters): self.m = m self.dndm0 = dndm0 self.wdm = wdm @@ -279,18 +281,18 @@ def __init__(self, wdm_mass=3.0, wdm_model=Viel05, wdm_params={}, #=========================================================================== # Parameters #=========================================================================== - @parameter + @parameter("model") def wdm_model(self, val): """ A model for the WDM effect on the transfer function. :type: str or :class:`WDM` subclass """ - if not np.issubclass_(val, WDM) and not isinstance(val, basestring): + if not np.issubclass_(val, WDM) and not isinstance(val, str): raise ValueError("wdm_model must be a WDM subclass or string, got %s" % type(val)) return val - @parameter + @parameter("param") def wdm_params(self, val): """ Parameters of the WDM model. @@ -299,7 +301,7 @@ def wdm_params(self, val): """ return val - @parameter + @parameter("param") def wdm_mass(self, val): """ Mass of the WDM particle. @@ -318,7 +320,7 @@ def wdm_mass(self, val): #=========================================================================== # Derived properties #=========================================================================== - @cached_property("cosmo", "wdm_mass", "wdm_model", "wdm_params") + @cached_quantity def wdm(self): """ The instantiated WDM model. @@ -328,11 +330,11 @@ def wdm(self): if np.issubclass_(self.wdm_model, WDM): return self.wdm_model(self.wdm_mass, self.cosmo,self.z, **self.wdm_params) - elif isinstance(self.wdm_transfer, basestring): + elif isinstance(self.wdm_transfer, str): return get_model(self.wdm_model, __name__, mx=self.wdm_mass, cosmo=self.cosmo, z=self.z,**self.wdm_params) - @cached_property("wdm") + @cached_quantity def _unnormalised_power(self): """ Normalised log power at :math:`z=0` @@ -358,26 +360,26 @@ def __init__(self, alter_dndm=None, alter_params=None, **kwargs): self.alter_dndm = alter_dndm self.alter_params = alter_params or {} - @parameter + @parameter("switch") def alter_dndm(self, val): """ A model for empirical recalibration of the HMF. :type: None, str, or :class`WDMRecalibrateMF` subclass. """ - if not np.issubclass_(val, WDMRecalibrateMF) and val is not None and not np.issubclass_(val, basestring): + if not np.issubclass_(val, WDMRecalibrateMF) and val is not None and not np.issubclass_(val, str): raise TypeError("alter_dndm must be a WDMRecalibrateMF subclass, string, or None") else: return val - @parameter + @parameter("param") def alter_params(self, val): """ Model parameters for `alter_dndm`. """ return val - @cached_property("alter_dndm", "M", "wdm", "alter_params") + @cached_quantity def dndm(self): """ The number density of haloes in WDM, ``len=len(m)`` [units :math:`h^4 M_\odot^{-1} Mpc^{-3}`] @@ -386,10 +388,10 @@ def dndm(self): if self.alter_dndm is not None: if np.issubclass_(self.alter_dndm, WDMRecalibrateMF): - alter = self.alter_dndm(M=self.M, dndm0=dndm, wdm=self.wdm, + alter = self.alter_dndm(m=self.m, dndm0=dndm, wdm=self.wdm, **self.alter_params) else: - alter = get_model(self.alter_dndm, __name__, M=self.M, + alter = get_model(self.alter_dndm, __name__, m=self.m, dndm0=dndm, wdm=self.wdm, **self.alter_params) dndm = alter.dndm_alter() diff --git a/setup.py b/setup.py index 54a5301..d869ec0 100644 --- a/setup.py +++ b/setup.py @@ -1,10 +1,12 @@ -from setuptools import setup +from setuptools import setup, find_packages import os import sys import re import io +test_suite = "nose.collector" + def read(*names, **kwargs): with io.open( os.path.join(os.path.dirname(__file__), *names), @@ -32,17 +34,20 @@ def find_version(*file_paths): setup( name="hmf", version=find_version("hmf", "__init__.py"), - packages=['hmf', 'hmf.fitting'], + packages=find_packages(), install_requires=["numpy>=1.6.2", "scipy>=0.12.0", - "astropy>=1.0"], + "astropy>=1.1"], scripts=["scripts/hmf", "scripts/hmf-fit"], author="Steven Murray", - author_email="steven.murray@uwa.edu.au", + author_email="steven.murray@curtin.edu.au", description="A halo mass function calculator", long_description=read('README.rst'), license="MIT", keywords="halo mass function", url="https://github.com/steven-murray/hmf", + classifiers=["Programming Language :: Python :: 2.7", + "Programming Language :: Python :: 3.6", + ] # could also include long_description, download_url, classifiers, etc. ) diff --git a/tests/_test_fit.py b/tests/_test_fit.py index ccde6cc..bbdf872 100644 --- a/tests/_test_fit.py +++ b/tests/_test_fit.py @@ -29,7 +29,7 @@ def test_circular_minimize(): guess=[0.9], blobs=None, verbose=0, store_class=False, relax=False) res = f.fit(h) - print "Diff: ", np.abs(res.x - 0.8) + print("Diff: ", np.abs(res.x - 0.8)) assert np.abs(res.x - 0.8) < 0.01 @@ -43,5 +43,5 @@ def test_circular_emcee(): verbose=0, store_class=False, relax=False) sampler = f.fit(h, nwalkers=16, nsamples=15, burnin=0, nthreads=0) - print "Diff: ", np.abs(np.mean(sampler.chain) - 0.8) + print("Diff: ", np.abs(np.mean(sampler.chain) - 0.8)) assert np.abs(np.mean(sampler.chain) - 0.8) < 0.01 diff --git a/tests/test_cosmo.py b/tests/test_cosmo.py index 796d95c..df3dee3 100644 --- a/tests/test_cosmo.py +++ b/tests/test_cosmo.py @@ -35,12 +35,12 @@ def test_cosmo_model(self): self.c.update(cosmo_model=WMAP7) assert self.c.cosmo.Om0 == 0.272 - print self.c.mean_density0 + print((self.c.mean_density0)) assert np.isclose(self.c.mean_density0, 75468972351.60081) def test_cosmo_params(self): self.c.update(cosmo_params={"H0":0.6}) - print self.c.cosmo.H0.value + print((self.c.cosmo.H0.value)) assert self.c.cosmo.H0.value == 0.6 self.c.update(cosmo_params={"Om0":0.2}) assert self.c.cosmo.Om0 == 0.2 diff --git a/tests/test_fcoll.py b/tests/test_fcoll.py index ca13801..a2d511c 100644 --- a/tests/test_fcoll.py +++ b/tests/test_fcoll.py @@ -14,6 +14,7 @@ from hmf import MassFunction from scipy.special import erfc + class TestFcoll(object): def check_fcoll(self, pert, fit): @@ -26,8 +27,8 @@ def check_fcoll(self, pert, fit): num = pert.rho_gtm / pert.mean_density0 err = np.abs((num - anl) / anl) - print np.max(err) - print num / anl - 1 + print(np.max(err)) + print(num / anl - 1) assert np.max(err) < 0.05 def test_fcolls(self): @@ -60,7 +61,7 @@ def check(self, hmf, minm, maxm): num = hmf.rho_gtm / hmf.mean_density0 err = np.abs((num - anl) / anl)[np.logical_and(hmf.m > 10 ** 10, hmf.m < 10 ** 15)] err = err[np.logical_not(np.isnan(err))] - print np.max(err) + print((np.max(err))) assert np.max(err) < TestCumulants.tol def test_ranges_not_cut(self): @@ -79,8 +80,8 @@ def test_ranges_cut(self): def check_mgtm(self, hmf, maxm): hmf.update(Mmin=0, Mmax=maxm, dlog10m=0.01) - print "rhogtm: ", hmf.rho_gtm - print "rhomean:", hmf.mean_density0 + print("rhogtm: ", hmf.rho_gtm) + print("rhomean:", hmf.mean_density0) assert np.abs(hmf.rho_gtm[0] / hmf.mean_density0 - 1) < 0.1 # THIS IS PRETTY BIG! @@ -91,7 +92,7 @@ def test_mgtm(self): def check_mltm(self, hmf, maxm): hmf.update(Mmin=3, Mmax=maxm, dlog10m=0.01) - print np.abs(hmf.rho_ltm[-1] / hmf.mean_density0 - 1) + print(np.abs(hmf.rho_ltm[-1] / hmf.mean_density0 - 1)) assert np.abs(hmf.rho_ltm[-1] / hmf.mean_density0 - 1) < 0.2 def test_mltm(self): diff --git a/tests/test_filters.py b/tests/test_filters.py index 258c671..1f9aa70 100644 --- a/tests/test_filters.py +++ b/tests/test_filters.py @@ -26,7 +26,7 @@ def test_sigma(self): R = 1.0 true =(9*R**2*sin(R)**2/2 + 9*R**2*cos(R)**2/2 + 9*R*sin(R)*cos(R)/2 - 9*sin(R)**2)/(2*pi**2*R**6) - print true,self.cls.sigma(R)**2 + print(true,self.cls.sigma(R)**2) assert np.isclose(self.cls.sigma(R)[0]**2,true) def test_sigma1(self): @@ -34,7 +34,7 @@ def test_sigma1(self): true = (3*R**2*sin(R)**2/2 + 3*R**2*cos(R)**2/2 + 9*R*sin(R)*cos(R)/2 - 9*sin(R)**2/4 + 45*cos(R)**2/4 - 45*sin(R)*cos(R)/(4*R))/(2*pi**2*R**6) - print true,self.cls.sigma(R,1)**2 + print(true,self.cls.sigma(R,1)**2) assert np.isclose(self.cls.sigma(R,1)[0]**2,true) def test_dwdlnkr(self): @@ -46,7 +46,7 @@ def test_dlnssdlnr(self): R = 1.0 true = 2*R**4*(-45*sin(R)**2/(4*R**2) - 27*cos(R)**2/(4*R**2) - 81*sin(R)*cos(R)/(4*R**3) + 27*sin(R)**2/R**4)/(9*R**2*sin(R)**2/2 + 9*R**2*cos(R)**2/2 + 9*R*sin(R)*cos(R)/2 - 9*sin(R)**2) - print true, self.cls.dlnss_dlnr(R) + print(true, self.cls.dlnss_dlnr(R)) assert np.isclose(self.cls.dlnss_dlnr(R),true) class TestSharpK(object): @@ -60,7 +60,7 @@ def test_sigma(self): t = 2 + 2 + 1 true = 1./(2*pi**2 * t * R**t) - print true,self.cls.sigma(R)**2 + print(true,self.cls.sigma(R)**2) assert np.isclose(self.cls.sigma(R)[0]**2,true) def test_sigma1(self): @@ -68,7 +68,7 @@ def test_sigma1(self): t = 4 + 2 + 1 true = 1./(2*pi**2 * t * R**t) - print true,self.cls.sigma(R,1)**2 + print(true,self.cls.sigma(R,1)**2) assert np.isclose(self.cls.sigma(R,1)[0]**2,true) @@ -78,7 +78,7 @@ def test_dlnssdlnr(self): sigma2 = 1./(2*pi**2 * t * R**t) true = -1./(2*pi**2*sigma2 * R**(3+2)) - print true, self.cls.dlnss_dlnr(R) + print(true, self.cls.dlnss_dlnr(R)) assert np.isclose(self.cls.dlnss_dlnr(R),true) @@ -87,7 +87,7 @@ def test_sigma_R3(self): t = 2 + 2 + 1 true = 1./(2*pi**2 * t * R**t) - print true,self.cls.sigma(R)**2 + print(true,self.cls.sigma(R)**2) assert np.isclose(self.cls.sigma(R)[0]**2,true) def test_sigma1_R3(self): @@ -95,7 +95,7 @@ def test_sigma1_R3(self): t = 4 + 2 + 1 true = 1./(2*pi**2 * t * R**t) - print true,self.cls.sigma(R,1)**2 + print(true,self.cls.sigma(R,1)**2) assert np.isclose(self.cls.sigma(R,1)[0]**2,true) @@ -105,7 +105,7 @@ def test_dlnssdlnr_R3(self): sigma2 = 1./(2*pi**2 * t * R**t) true = -1./(2*pi**2*sigma2 * R**(3+2)) - print true, self.cls.dlnss_dlnr(R) + print(true, self.cls.dlnss_dlnr(R)) assert np.isclose(self.cls.dlnss_dlnr(R),true) @@ -119,7 +119,7 @@ def test_sigma_Rhalf(self): with warnings.catch_warnings(record=True) as w: s2 = self.cls.sigma(R)[0]**2 assert w - print s2,true + print(s2,true) assert np.isclose(s2,true) @@ -135,7 +135,7 @@ def test_sigma1_Rhalf(self): s2 = self.cls.sigma(R,1)[0]**2 assert w - print s2,true + print(s2,true) assert np.isclose(s2,true) @@ -145,7 +145,7 @@ def test_dlnssdlnr_Rhalf(self): sigma2 = 1./(2*pi**2 * t * R**t) true = -1./(2*pi**2*sigma2 * R**(3+2)) - print true, self.cls.dlnss_dlnr(R) + print(true, self.cls.dlnss_dlnr(R)) assert np.isclose(self.cls.dlnss_dlnr(R),true) @@ -160,14 +160,14 @@ def test_sigma(self): R = 10.0 true = 3./(16*pi**(3./2.)*R**5) - print true,self.cls.sigma(R)**2 + print(true,self.cls.sigma(R)**2) assert np.isclose(self.cls.sigma(R)[0]**2,true) def test_sigma1(self): R = 10.0 true = 15/(32*pi**(3./2.)*R**7) - print true,self.cls.sigma(R,1)**2 + print(true,self.cls.sigma(R,1)**2) assert np.isclose(self.cls.sigma(R,1)[0]**2,true) @@ -175,5 +175,5 @@ def test_dlnssdlnr(self): R = 10.0 true = -5 - print true, self.cls.dlnss_dlnr(R) + print(true, self.cls.dlnss_dlnr(R)) assert np.isclose(self.cls.dlnss_dlnr(R),true) diff --git a/tests/test_fits.py b/tests/test_fits.py index 074b749..f3d269a 100644 --- a/tests/test_fits.py +++ b/tests/test_fits.py @@ -4,7 +4,7 @@ import inspect import os LOCATION = "/".join(os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))).split("/")[:-1]) -# from nose.tools import raises +from nose.tools import raises import sys sys.path.insert(0, LOCATION) from hmf import MassFunction @@ -34,6 +34,8 @@ def __init__(self): def test_max_lt_ps(self): for redshift in [0.0, 2.0]: for fit in allfits: + if fit is ff.PS: + continue # if fit is ff.AnguloBound: # continue yield self.check_form, fit, redshift @@ -44,3 +46,36 @@ def check_form(self,fit,redshift): assert self.ps_max >= self.hmf.fsigma[maxarg] assert np.all(np.diff(self.hmf.fsigma[:maxarg]) >= 0) assert np.all(np.diff(self.hmf.fsigma[maxarg:]) <= 0) + +def test_tinker08_dh(): + h = MassFunction(hmf_model="Tinker08", delta_h=200) + h1 = MassFunction(hmf_model="Tinker08", delta_h=200.1) + + assert np.allclose(h.fsigma,h1.fsigma,rtol=1e-2) + +def test_tinker10_dh(): + h = MassFunction(hmf_model="Tinker10", delta_h=200) + h1 = MassFunction(hmf_model="Tinker10", delta_h=200.1) + + assert np.allclose(h.fsigma,h1.fsigma,rtol=1e-2) + +@raises(ValueError) +def test_tinker10_neg_gam(): + h = MassFunction(hmf_model="Tinker10", hmf_params={"gamma_200":-1}) + h.fsigma + +@raises(ValueError) +def test_tinker10_neg_eta(): + h = MassFunction(hmf_model="Tinker10", hmf_params={"eta_200":-1}) + h.fsigma + +@raises(ValueError) +def test_tinker10_neg_etaphi(): + h = MassFunction(hmf_model="Tinker10", hmf_params={"eta_200":-1, "phi_200":0}) + h.fsigma + +@raises(ValueError) +def test_tinker10_neg_beta(): + h = MassFunction(hmf_model="Tinker10", hmf_params={"beta_200":-1}) + h.fsigma + diff --git a/tests/test_framework.py b/tests/test_framework.py new file mode 100644 index 0000000..4a80a2d --- /dev/null +++ b/tests/test_framework.py @@ -0,0 +1,46 @@ +import inspect +import os + +LOCATION = "/".join(os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))).split("/")[:-1]) +from nose.tools import raises +import sys +sys.path.insert(0, LOCATION) +from hmf import hmf + +@raises(TypeError) +def test_incorrect_argument(): + t = hmf.MassFunction(wrong_arg=3) + +@raises(ValueError) +def test_incorrect_update_arg(): + t = hmf.MassFunction() + t.update(wrong_arg=3) + + +class TestIntrospection(object): + def __init__(self): + self.cls = hmf.MassFunction + self.inst = self.cls(z=10) + + def test_parameter_names(self): + assert "cosmo_model" in self.cls.get_all_parameter_names() + + def test_parameter_defaults(self): + assert type(self.cls.get_all_parameter_defaults(recursive=False)) is dict + + assert self.cls.get_all_parameter_defaults()['z'] == 0 + + def test_parameter_default_rec(self): + pd = self.cls.get_all_parameter_defaults(recursive=True) + assert type(pd['cosmo_params']) is dict + + def test_param_values(self): + assert type(self.inst.parameter_values) is dict + assert self.inst.parameter_values['z'] == 10 + + def test_qnt_avail(self): + assert 'dndm' in self.cls.quantities_available() + + def test_parameter_info(self): + assert self.cls.parameter_info() is None + assert self.cls.parameter_info(names=['z']) is None \ No newline at end of file diff --git a/tests/test_functional.py b/tests/test_functional.py index 88195fc..6b786b2 100644 --- a/tests/test_functional.py +++ b/tests/test_functional.py @@ -23,8 +23,9 @@ def test_order(): "sigma.8: 0.7, PS, z: 2", "sigma.8: 0.8, PS, z: 2"] - for i,(quants, mf, label) in enumerate(tf.get_hmf(['dndm','ngtm'],z=range(3),hmf_model=["ST","PS"],sigma_8=[0.7,0.8])): - assert label==order[i] + for i,(quants, mf, label) in enumerate(tf.get_hmf(['dndm','ngtm'],z=list(range(3)),hmf_model=["ST","PS"],sigma_8=[0.7,0.8])): + assert len(label)==len(order[i]) + assert sorted(label)==sorted(order[i]) assert isinstance(mf,MassFunction) assert np.allclose(quants[0],mf.dndm) assert np.allclose(quants[1],mf.ngtm) diff --git a/tests/test_genmf.py b/tests/test_genmf.py index 6bbeae2..d728aa8 100644 --- a/tests/test_genmf.py +++ b/tests/test_genmf.py @@ -60,7 +60,7 @@ def rms_diff(vec1, vec2, tol): vec1 = vec1[mask] vec2 = vec2[mask] err = np.sqrt(np.mean(((vec1 - vec2) / vec2) ** 2)) - print "RMS Error: ", err, "(> ", tol, ")" + print("RMS Error: ", err, "(> ", tol, ")") return err < tol def max_diff_rel(vec1, vec2, tol): @@ -68,7 +68,7 @@ def max_diff_rel(vec1, vec2, tol): vec1 = vec1[mask] vec2 = vec2[mask] err = np.max(np.abs((vec1 - vec2) / vec2)) - print "Max Diff: ", err, "(> ", tol, ")" + print("Max Diff: ", err, "(> ", tol, ")") return err < tol def max_diff(vec1, vec2, tol): @@ -76,7 +76,7 @@ def max_diff(vec1, vec2, tol): vec1 = vec1[mask] vec2 = vec2[mask] err = np.max(np.abs((vec1 - vec2))) - print "Max Diff: ", err, "(> ", tol, ")" + print("Max Diff: ", err, "(> ", tol, ")") return err < tol #=============================================================================== diff --git a/tests/test_growth.py b/tests/test_growth.py index d0de88f..5d6d461 100644 --- a/tests/test_growth.py +++ b/tests/test_growth.py @@ -19,7 +19,7 @@ def __init__(self): def test_gf(self): for z in np.arange(0,8,0.5): - print self.gf.growth_factor(z),self.genf.growth_factor(z) + print(self.gf.growth_factor(z),self.genf.growth_factor(z)) assert np.isclose(self.gf.growth_factor(z),self.genf.growth_factor(z),rtol=1e-2 + z/500.0) def test_gr(self): @@ -31,14 +31,14 @@ def test_gfunc(self): gf_func = self.gf.growth_factor_fn(0.0) genf_func = self.genf.growth_factor_fn(0.0) - print gf_func(np.linspace(0,5,10)),genf_func(np.linspace(0,5,10)) + print(gf_func(np.linspace(0,5,10)),genf_func(np.linspace(0,5,10))) assert np.allclose(gf_func(np.linspace(0,5,10)),genf_func(np.linspace(0,5,10)),rtol=1e-2) def test_gr_func(self): gr_func = self.gf.growth_rate_fn(0.0) genf_func = self.genf.growth_rate_fn(0.0) - print gr_func(np.linspace(0,5,10)),genf_func(np.linspace(0,5,10)) + print(gr_func(np.linspace(0,5,10)),genf_func(np.linspace(0,5,10))) assert np.allclose(gr_func(np.linspace(0,5,10)),genf_func(np.linspace(0,5,10)),rtol=1e-2) def test_inverse(self): @@ -46,5 +46,5 @@ def test_inverse(self): genf_func = self.genf.growth_factor_fn(0.0,inverse=True) gf = np.linspace(0.15,0.99,10) - print gf_func(gf),genf_func(gf) + print(gf_func(gf),genf_func(gf)) assert np.allclose(gf_func(gf),genf_func(gf),rtol=1e-1) \ No newline at end of file diff --git a/tests/test_halofit.py b/tests/test_halofit.py new file mode 100644 index 0000000..035014e --- /dev/null +++ b/tests/test_halofit.py @@ -0,0 +1,44 @@ +import inspect +import os + +LOCATION = "/".join(os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))).split("/")[:-1]) +from nose.tools import raises +import sys +sys.path.insert(0, LOCATION) +from hmf import transfer +import numpy as np + + +def test_takahashi(): + t = transfer.Transfer(transfer_model="EH", takahashi=False, lnk_max=7) + tt = transfer.Transfer(transfer_model="EH", takahashi=True, lnk_max=7) + + assert np.isclose(t.nonlinear_power[0],tt.nonlinear_power[0],rtol=1e-4) + print(t.nonlinear_power[-1]/ tt.nonlinear_power[-1]) + assert np.logical_not(np.isclose(t.nonlinear_power[-1]/ tt.nonlinear_power[-1], 1, rtol=0.4)) + + +def test_takahashi_hiz(): + # This test should do the HALOFIT WARNING + t = transfer.Transfer(transfer_model="EH", takahashi=False, lnk_max=7,z=8.0) + tt = transfer.Transfer(transfer_model="EH", takahashi=True, lnk_max=7, z=8.0) + + assert np.isclose(t.nonlinear_power[0],tt.nonlinear_power[0],rtol=1e-4) + print(t.nonlinear_power[-1]/ tt.nonlinear_power[-1]) + assert np.logical_not(np.isclose(t.nonlinear_power[-1]/ tt.nonlinear_power[-1], 1, rtol=0.4)) + + t.update(z=0) + + assert np.logical_not(np.isclose(t.nonlinear_power[0]/ tt.nonlinear_power[0], 0.9, rtol=0.1)) + assert np.logical_not(np.isclose(t.nonlinear_power[-1]/ tt.nonlinear_power[-1], 0.99, rtol=0.1)) + + + +def test_halofit_high_s8(): + t = transfer.Transfer(transfer_model="EH", lnk_max=7,sigma_8=0.999) + thi = transfer.Transfer(transfer_model="EH", lnk_max=7, sigma_8=1.001) #just above threshold + + + print(t.nonlinear_power[0]/thi.nonlinear_power[0] -1, t.nonlinear_power[-1]/thi.nonlinear_power[-1] -1) + assert np.isclose(t.nonlinear_power[0],thi.nonlinear_power[0],rtol=2e-2) + assert np.isclose(t.nonlinear_power[-1], thi.nonlinear_power[-1], rtol=5e-2) diff --git a/tests/test_hmf.py b/tests/test_hmf.py index ed7033f..5bd30f6 100644 --- a/tests/test_hmf.py +++ b/tests/test_hmf.py @@ -5,6 +5,10 @@ import sys sys.path.insert(0, LOCATION) from hmf import MassFunction +from hmf.filters import TopHat +import numpy as np +import warnings + # @raises(ValueError) # def test_wrong_fit(): @@ -29,3 +33,73 @@ def test_delta_halo_mean(): def test_delta_halo_crit(): hmf = MassFunction(delta_h=180, delta_wrt="crit", cosmo_params={"Om0":0.3}) assert abs(hmf.delta_halo - 600.0) < 1e-3 + +@raises(ValueError) +def test_wrong_filter(): + h = MassFunction(filter_model=2) + + +@raises(ValueError) +def test_string_dc(): + h = MassFunction(delta_c="this") + + +@raises(ValueError) +def test_neg_dc(): + h = MassFunction(delta_c=-1) + + +@raises(ValueError) +def test_big_dc(): + h = MassFunction(delta_c=20.) + + +@raises(ValueError) +def test_wrong_fit(): + h = MassFunction(hmf_model=1) + + +@raises(ValueError) +def test_wrong_mf_par(): + h = MassFunction(hmf_params=2) + +@raises(ValueError) +def test_wrong_dh(): + h = MassFunction(delta_h="string") + + +@raises(ValueError) +def test_neg_dh(): + h = MassFunction(delta_h=0) + + +@raises(ValueError) +def test_big_dh(): + h = MassFunction(delta_h=1e5) + + +@raises(ValueError) +def test_delta_wrt(): + h = MassFunction(delta_wrt="this") + + +def test_str_filter(): + h = MassFunction(filter_model="TopHat") + h_ = MassFunction(filter_model="TopHat") + + assert np.allclose(h.sigma,h_.sigma) + + + + +def test_mass_nonlinear_outside_range(): + h = MassFunction(Mmin=8,Mmax=9) + with warnings.catch_warnings(record=True) as w: + assert h.mass_nonlinear>0 + assert len(w) + + + + + + diff --git a/tests/test_integrate_hmf.py b/tests/test_integrate_hmf.py index 56b7e59..3aba179 100644 --- a/tests/test_integrate_hmf.py +++ b/tests/test_integrate_hmf.py @@ -4,9 +4,19 @@ LOCATION = "/".join(os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))).split("/")[:-1]) import sys sys.path.insert(0, LOCATION) -from mrpy.special import gammainc +from mpmath import gammainc as _mp_ginc + from hmf.integrate_hmf import hmf_integral_gtm +def _flt(a): + try: + return a.astype('float') + except AttributeError: + return float(a) +_ginc_ufunc = np.frompyfunc(lambda z, x: _mp_ginc(z, x), 2, 1) +def gammainc(z,x): + return _flt(_ginc_ufunc(z,x)) + class TestAnalyticIntegral(object): def __init__(self): pass @@ -41,7 +51,7 @@ def test_high_z(self): dndm = self.tggd(m,9.0,-1.93,0.4) ngtm = self.anl_int(m,9.0,-1.93,0.4) - print ngtm/hmf_integral_gtm(m,dndm) + print(ngtm/hmf_integral_gtm(m,dndm)) assert np.allclose(ngtm,hmf_integral_gtm(m,dndm),rtol=0.03) # def test_low_mmax_z0(self): @@ -57,5 +67,5 @@ def test_low_mmax_high_z(self): dndm = self.tggd(m,9.0,-1.93,0.4) ngtm = self.anl_int(m,9.0,-1.93,0.4) - print ngtm/hmf_integral_gtm(m,dndm) + print(ngtm/hmf_integral_gtm(m,dndm)) assert np.allclose(ngtm,hmf_integral_gtm(m,dndm),rtol=0.03) diff --git a/tests/test_sample.py b/tests/test_sample.py index d35711d..555d136 100644 --- a/tests/test_sample.py +++ b/tests/test_sample.py @@ -15,7 +15,7 @@ def test_circular(): s = spline(np.log10(h.m),np.log10(h.dndm)) - print hist, 10**s(centres) + print(hist, 10**s(centres)) assert np.allclose(hist,10**s(centres),rtol=0.05) diff --git a/tests/test_transfer.py b/tests/test_transfer.py index 42f0c8f..f2d8e19 100644 --- a/tests/test_transfer.py +++ b/tests/test_transfer.py @@ -5,11 +5,11 @@ import sys sys.path.insert(0, LOCATION) from hmf.transfer import Transfer -from astropy.cosmology import LambdaCDM from hmf.transfer_models import EH_BAO + def rms(a): - print a - print "RMS: ", np.sqrt(np.mean(np.square(a))) + print(a) + print("RMS: ", np.sqrt(np.mean(np.square(a)))) return np.sqrt(np.mean(np.square(a))) def check_close(t, t2, fit): @@ -24,27 +24,39 @@ def test_updates(): t = Transfer() t2 = Transfer() for k, v in {"z":0.1, - "transfer_params":{"lAccuracyBoost":1.5, - "AccuracyBoost":1.5}, "sigma_8":0.82, "n":0.95, - "cosmo_params":{"H0":68.0}}.iteritems(): + "cosmo_params":{"H0":68.0}}.items(): yield check_update, t, t2, k, v def test_halofit(): t = Transfer(lnk_min=-20, lnk_max=20, dlnk=0.05, transfer_model="EH") - print EH_BAO._defaults - print "in test_transfer, params are: ", t.transfer_params + print(EH_BAO._defaults) + print("in test_transfer, params are: ", t.transfer_params) assert np.isclose(t.power[0],t.nonlinear_power[0]) assert 5 * t.power[-1] < t.nonlinear_power[-1] -def test_data(): - t = Transfer(cosmo_model=LambdaCDM(Om0=0.3, Ode0=0.7, H0=70.0, Ob0=0.05), sigma_8=0.8, - n=1, transfer_params={"transfer__k_per_logint":0, "transfer__kmax":100.0}, - lnk_min=np.log(1e-11), lnk_max=np.log(1e11)) - tdata = np.genfromtxt(LOCATION + "/data/transfer_for_hmf_tests.dat") - pdata = np.genfromtxt(LOCATION + "/data/power_for_hmf_tests.dat") - assert rms(t._unnormalised_lnT - np.log(tdata[:, 1])) < 0.05 # Does better than 0.001 on my system... - diff = t.power - pdata[:, 1] - print t._unnormalised_lnT[400], t._unnormalised_power[400], t._power0[400] - assert rms(t.power - pdata[:, 1]) < 0.001 +def test_ehnobao(): + t = Transfer(transfer_model="EH") + tnobao = Transfer(transfer_model="EH_NoBAO") + + assert np.isclose(t._unnormalised_lnT[0], tnobao._unnormalised_lnT[0],rtol=1e-5) + +def test_bondefs(): + t = Transfer(transfer_model="BondEfs") + print(np.exp(t._unnormalised_lnT)) + assert np.isclose(np.exp(t._unnormalised_lnT[0]),1,rtol=1e-5) +# Following test is too slow... and would need to be updated whenever CAMB is updated... + +# def test_data(): +# cp = camb.CAMBparams() +# cp.set_matter_power(kmax=100.) +# t = Transfer(cosmo_model=LambdaCDM(Om0=0.3, Ode0=0.7, H0=70.0, Ob0=0.05), sigma_8=0.8, +# n=1, transfer_params={"camb_params":cp}, +# lnk_min=np.log(1e-11), lnk_max=np.log(1e11)) +# tdata = np.genfromtxt(LOCATION + "/data/transfer_for_hmf_tests.dat") +# pdata = np.genfromtxt(LOCATION + "/data/power_for_hmf_tests.dat") +# #assert rms(t._unnormalised_lnT - np.log(tdata[:, 1])) < 0.05 # Does better than 0.001 on my system... +# diff = t.power - pdata[:, 1] +# #print(t._unnormalised_lnT[400], t._unnormalised_power[400], t._power0[400]) +# assert rms(t.power - pdata[:, 1]) < 0.001 diff --git a/tests/test_wdm.py b/tests/test_wdm.py new file mode 100644 index 0000000..dec7769 --- /dev/null +++ b/tests/test_wdm.py @@ -0,0 +1,97 @@ +import inspect +import os + +LOCATION = "/".join(os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))).split("/")[:-1]) +from nose.tools import raises +import sys +sys.path.insert(0, LOCATION) +from hmf import wdm, hmf +import numpy as np + + +@raises(NotImplementedError) +def test_null(): + w = wdm.WDM(mx=1.0) + w.transfer(1.) + + +class TestViel(object): + def __init__(self): + self.cls = wdm.Viel05(mx=1.0) + + def test_lowk_transfer(self): + assert np.isclose(self.cls.transfer(1e-5),1,rtol=1e-4) + + def test_lam_eff(self): + assert self.cls.lam_eff_fs > 0 + + def test_m_eff(self): + assert self.cls.m_fs > 0 + + def test_lam_hm(self): + assert self.cls.lam_hm > self.cls.lam_eff_fs + + def test_m_hm(self): + assert self.cls.m_hm > self.cls.m_fs + + +class TestBode(TestViel): + def __init__(self): + self.cls = wdm.Bode01(mx=1.0) + + +class TestSchneider12_vCDM(object): + def __init__(self): + self.cdm = hmf.MassFunction() + self.cls = wdm.Schneider12_vCDM(m = self.cdm.m, dndm0 = self.cdm.dndm) + + def test_high_m(self): + assert np.isclose(self.cls.dndm_alter()[-1], self.cdm.dndm[-1], rtol = 1e-3) + + +class TestSchneider12(TestSchneider12_vCDM): + def __init__(self): + self.cdm = hmf.MassFunction() + self.cls = wdm.Schneider12(m=self.cdm.m, dndm0=self.cdm.dndm) + +class TestLovell14(TestSchneider12_vCDM): + def __init__(self): + self.cdm = hmf.MassFunction() + self.cls = wdm.Lovell14(m=self.cdm.m, dndm0=self.cdm.dndm) + + + +class TestTransfer(object): + def __init__(self): + self.wdm = wdm.TransferWDM(wdm_mass=3.0, wdm_model=wdm.Viel05) + self.cdm = hmf.MassFunction() + + def test_wdm_model(self): + assert isinstance(self.wdm.wdm, wdm.Viel05) + + @raises(ValueError) + def test_wrong_model_type(self): + cls = wdm.TransferWDM(wdm_mass=3.0, wdm_model=3) + + def test_power(self): + print(self.wdm.power[0],self.cdm.power[0],self.wdm.power[0]/self.cdm.power[0]-1) + assert np.isclose(self.wdm.power[0], self.cdm.power[0],rtol=1e-4) + assert self.wdm.power[-1] < self.cdm.power[-1] + + + +class TestMassFunction(object): + + def __init__(self): + self.wdm = wdm.MassFunctionWDM(alter_dndm=None, wdm_mass=3.0, wdm_model=wdm.Viel05) + self.cdm = hmf.MassFunction() + + def test_dndm(self): + assert np.isclose(self.cdm.dndm[-1],self.wdm.dndm[-1], rtol=1e-3) + assert self.cdm.dndm[0] > self.wdm.dndm[0] + + +class TestMassFunctionAlter(TestMassFunction): + def __init__(self): + self.wdm = wdm.MassFunctionWDM(alter_dndm=wdm.Schneider12_vCDM, wdm_mass=3.0, wdm_model=wdm.Viel05) + self.cdm = hmf.MassFunction() \ No newline at end of file