diff --git a/.travis.yml b/.travis.yml index 60f8eaa5..8fc78a32 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,35 +1,49 @@ +os: + - linux + - osx + language: c sudo: false -addons: - apt: - sources: - - ubuntu-toolchain-r-test + +branches: + only: + - master install: - - bash -x devtools/travis-ci/install.sh + - source devtools/travis-ci/install.sh - export PYTHONUNBUFFERED=true - - export PATH=$HOME/miniconda3/bin:$PATH -script: +before_script: + # Mimic X display + - "export DISPLAY=:99.0" + - #"sh -e /etc/init.d/xvfb start" + - sleep 3 # give xvfb some time to start - # Create an environment for testing +script: + # Add org channel and force update + - conda config --add channels omnia --add channels conda-forge + - conda update --yes --all + # Create a test environment - conda create --yes -n test python=$python + # Activate the test environment - source activate test - # Add omnia channel for packages - - conda config --add channels omnia - - # Install package dependencies - - conda install --yes --file devtools/travis-ci/requirements.txt - - python setup.py install - + # Build the recipe + - conda build devtools/conda-recipe + # Install the package + - conda install --yes --use-local ${PACKAGENAME}-dev + # Install testing dependencies + - conda install --yes --quiet nose nose-timer # Test the package - - nosetests -v + - cd devtools && nosetests $PACKAGENAME --nocapture --verbosity=2 --with-timer -a '!slow' && cd .. env: matrix: - python=2.7 CONDA_PY=27 - python=3.5 CONDA_PY=35 - - python=3.4 CONDA_PY=34 + - python=3.6 CONDA_PY=36 + global: + - MPLBACKEND="Agg" + - PACKAGENAME="bitc" after_success: - bash -x devtools/travis-ci/after_success.sh diff --git a/basesetup.py b/basesetup.py new file mode 100644 index 00000000..1fbd3847 --- /dev/null +++ b/basesetup.py @@ -0,0 +1,329 @@ +from __future__ import print_function, absolute_import +import os +import sys +import json +import string +import shutil +import subprocess +import tempfile +from distutils.dep_util import newer_group +from distutils.core import Extension +from distutils.errors import DistutilsExecError +from distutils.ccompiler import new_compiler +from distutils.sysconfig import customize_compiler, get_config_vars +from distutils.command.build_ext import build_ext as _build_ext + + +def find_packages(): + """Find all of bitc's python packages. + Adapted from IPython's setupbase.py. Copyright IPython + contributors, licensed under the BSD license. + """ + packages = ['bitc.scripts'] + for dir,subdirs,files in os.walk('bitc'): + package = dir.replace(os.path.sep, '.') + if '__init__.py' not in files: + # not a package + continue + packages.append(package) + return packages + + + +################################################################################ +# Detection of compiler capabilities +################################################################################ + +class CompilerDetection(object): + # Necessary for OSX. See https://github.com/assaytools/assaytools/issues/576 + # The problem is that distutils.sysconfig.customize_compiler() + # is necessary to properly invoke the correct compiler for this class + # (otherwise the CC env variable isn't respected). Unfortunately, + # distutils.sysconfig.customize_compiler() DIES on OSX unless some + # appropriate initialization routines have been called. This line + # has a side effect of calling those initialzation routes, and is therefor + # necessary for OSX, even though we don't use the result. + _DONT_REMOVE_ME = get_config_vars() + + def __init__(self, disable_openmp): + cc = new_compiler() + customize_compiler(cc) + + self.msvc = cc.compiler_type == 'msvc' + self._print_compiler_version(cc) + + if disable_openmp: + self.openmp_enabled = False + else: + self.openmp_enabled, openmp_needs_gomp = self._detect_openmp() + self.sse3_enabled = self._detect_sse3() if not self.msvc else True + self.sse41_enabled = self._detect_sse41() if not self.msvc else True + + self.compiler_args_sse2 = ['-msse2'] if not self.msvc else ['/arch:SSE2'] + self.compiler_args_sse3 = ['-mssse3'] if (self.sse3_enabled and not self.msvc) else [] + + self.compiler_args_sse41, self.define_macros_sse41 = [], [] + if self.sse41_enabled: + self.define_macros_sse41 = [('__SSE4__', 1), ('__SSE4_1__', 1)] + if not self.msvc: + self.compiler_args_sse41 = ['-msse4'] + + if self.openmp_enabled: + self.compiler_libraries_openmp = [] + + if self.msvc: + self.compiler_args_openmp = ['/openmp'] + else: + self.compiler_args_openmp = ['-fopenmp'] + if openmp_needs_gomp: + self.compiler_libraries_openmp = ['gomp'] + else: + self.compiler_libraries_openmp = [] + self.compiler_args_openmp = [] + + if self.msvc: + self.compiler_args_opt = ['/O2'] + else: + self.compiler_args_opt = ['-O3', '-funroll-loops'] + print() + + def _print_compiler_version(self, cc): + print("C compiler:") + try: + if self.msvc: + if not cc.initialized: + cc.initialize() + cc.spawn([cc.cc]) + else: + cc.spawn([cc.compiler[0]] + ['-v']) + except DistutilsExecError: + pass + + def hasfunction(self, funcname, include=None, libraries=None, extra_postargs=None): + # running in a separate subshell lets us prevent unwanted stdout/stderr + part1 = ''' +from __future__ import print_function +import os +import json +from distutils.ccompiler import new_compiler +from distutils.sysconfig import customize_compiler, get_config_vars + +FUNCNAME = json.loads('%(funcname)s') +INCLUDE = json.loads('%(include)s') +LIBRARIES = json.loads('%(libraries)s') +EXTRA_POSTARGS = json.loads('%(extra_postargs)s') + ''' % { + 'funcname': json.dumps(funcname), + 'include': json.dumps(include), + 'libraries': json.dumps(libraries or []), + 'extra_postargs': json.dumps(extra_postargs)} + + part2 = ''' +get_config_vars() # DON'T REMOVE ME +cc = new_compiler() +customize_compiler(cc) +for library in LIBRARIES: + cc.add_library(library) + +status = 0 +try: + with open('func.c', 'w') as f: + if INCLUDE is not None: + f.write('#include %s\\n' % INCLUDE) + f.write('int main(void) {\\n') + f.write(' %s;\\n' % FUNCNAME) + f.write('}\\n') + objects = cc.compile(['func.c'], output_dir='.', + extra_postargs=EXTRA_POSTARGS) + cc.link_executable(objects, 'a.out') +except Exception as e: + status = 1 +exit(status) + ''' + tmpdir = tempfile.mkdtemp(prefix='hasfunction-') + try: + curdir = os.path.abspath(os.curdir) + os.chdir(tmpdir) + with open('script.py', 'w') as f: + f.write(part1 + part2) + proc = subprocess.Popen( + [sys.executable, 'script.py'], + stderr=subprocess.PIPE, stdout=subprocess.PIPE) + proc.communicate() + status = proc.wait() + finally: + os.chdir(curdir) + shutil.rmtree(tmpdir) + + return status == 0 + + def _print_support_start(self, feature): + print('Attempting to autodetect {0:6} support...'.format(feature), end=' ') + + def _print_support_end(self, feature, status): + if status is True: + print('Compiler supports {0}'.format(feature)) + else: + print('Did not detect {0} support'.format(feature)) + + def _detect_openmp(self): + self._print_support_start('OpenMP') + hasopenmp = self.hasfunction('omp_get_num_threads()', extra_postargs=['-fopenmp', '/openmp']) + needs_gomp = hasopenmp + if not hasopenmp: + hasopenmp = self.hasfunction('omp_get_num_threads()', libraries=['gomp']) + needs_gomp = hasopenmp + self._print_support_end('OpenMP', hasopenmp) + return hasopenmp, needs_gomp + + def _detect_sse3(self): + "Does this compiler support SSE3 intrinsics?" + self._print_support_start('SSE3') + result = self.hasfunction('__m128 v; _mm_hadd_ps(v,v)', + include='', + extra_postargs=['-msse3']) + self._print_support_end('SSE3', result) + return result + + def _detect_sse41(self): + "Does this compiler support SSE4.1 intrinsics?" + self._print_support_start('SSE4.1') + result = self.hasfunction( '__m128 v; _mm_round_ps(v,0x00)', + include='', + extra_postargs=['-msse4']) + self._print_support_end('SSE4.1', result) + return result + +################################################################################ +# Writing version control information to the module +################################################################################ + +def git_version(): + # Return the git revision as a string + # copied from numpy setup.py + def _minimal_ext_cmd(cmd): + # construct minimal environment + env = {} + for k in ['SYSTEMROOT', 'PATH']: + v = os.environ.get(k) + if v is not None: + env[k] = v + # LANGUAGE is used on win32 + env['LANGUAGE'] = 'C' + env['LANG'] = 'C' + env['LC_ALL'] = 'C' + out = subprocess.Popen( + cmd, stdout=subprocess.PIPE, env=env).communicate()[0] + return out + + try: + out = _minimal_ext_cmd(['git', 'rev-parse', 'HEAD']) + GIT_REVISION = out.strip().decode('ascii') + except OSError: + GIT_REVISION = 'Unknown' + + return GIT_REVISION + + +def write_version_py(VERSION, ISRELEASED, filename='bitc/version.py'): + cnt = """ +# THIS FILE IS GENERATED FROM BAYESIAN ITC SETUP.PY +short_version = '%(version)s' +version = '%(version)s' +full_version = '%(full_version)s' +git_revision = '%(git_revision)s' +release = %(isrelease)s + +if not release: + version = full_version +""" + # Adding the git rev number needs to be done inside write_version_py(), + # otherwise the import of numpy.version messes up the build under Python 3. + FULLVERSION = VERSION + if os.path.exists('.git'): + GIT_REVISION = git_version() + else: + GIT_REVISION = 'Unknown' + + if not ISRELEASED: + FULLVERSION += '.dev-' + GIT_REVISION[:7] + + a = open(filename, 'w') + try: + a.write(cnt % {'version': VERSION, + 'full_version': FULLVERSION, + 'git_revision': GIT_REVISION, + 'isrelease': str(ISRELEASED)}) + finally: + a.close() + + +class StaticLibrary(Extension): + def __init__(self, *args, **kwargs): + self.export_include = kwargs.pop('export_include', []) + Extension.__init__(self, *args, **kwargs) + + +class build_ext(_build_ext): + + def build_extension(self, ext): + if isinstance(ext, StaticLibrary): + self.build_static_extension(ext) + else: + _build_ext.build_extension(self, ext) + + def build_static_extension(self, ext): + from distutils import log + + sources = ext.sources + if sources is None or not isinstance(sources, (list, tuple)): + raise DistutilsSetupError( + ("in 'ext_modules' option (extension '%s'), " + + "'sources' must be present and must be " + + "a list of source filenames") % ext.name) + sources = list(sources) + + ext_path = self.get_ext_fullpath(ext.name) + depends = sources + ext.depends + if not (self.force or newer_group(depends, ext_path, 'newer')): + log.debug("skipping '%s' extension (up-to-date)", ext.name) + return + else: + log.info("building '%s' extension", ext.name) + + extra_args = ext.extra_compile_args or [] + macros = ext.define_macros[:] + for undef in ext.undef_macros: + macros.append((undef,)) + objects = self.compiler.compile(sources, + output_dir=self.build_temp, + macros=macros, + include_dirs=ext.include_dirs, + debug=self.debug, + extra_postargs=extra_args, + depends=ext.depends) + self._built_objects = objects[:] + if ext.extra_objects: + objects.extend(ext.extra_objects) + extra_args = ext.extra_link_args or [] + + language = ext.language or self.compiler.detect_language(sources) + + libname = os.path.basename(ext_path).split(os.extsep)[0] + output_dir = os.path.dirname(ext_path) + + if (self.compiler.static_lib_format.startswith('lib') and + libname.startswith('lib')): + libname = libname[3:] + + if not os.path.exists(output_dir): + # necessary for windows + os.makedirs(output_dir) + + self.compiler.create_static_lib(objects, + output_libname=libname, + output_dir=output_dir, + target_lang=language) + + for item in ext.export_include: + shutil.copy(item, output_dir) diff --git a/devtools/conda-recipe/build.sh b/devtools/conda-recipe/build.sh new file mode 100644 index 00000000..c8a24fc8 --- /dev/null +++ b/devtools/conda-recipe/build.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +pip install . diff --git a/devtools/conda-recipe/meta.yaml b/devtools/conda-recipe/meta.yaml new file mode 100644 index 00000000..6633812a --- /dev/null +++ b/devtools/conda-recipe/meta.yaml @@ -0,0 +1,39 @@ +package: + name: bitc-dev + version: 0.0.0 + +source: + path: ../.. + +requirements: + build: + - setuptools + - python + - numpy + - nose + - pandas + - pymc + - pint + - docopt + - schema + - scikit-learn + + run: + - python + - numpy + - nose + - pandas + - pymc + - pint + - docopt + - schema + - scikit-learn + +test: + imports: + - bitc + +about: + home: https://github.com/choderalab/bayesian-itc/ + license: GNU General Public License v3 or later (GPLv3+) + summary: Python tools for the analysis and modeling of isothermal titration calorimetry (ITC) experiments (development snapshot) diff --git a/devtools/travis-ci/install.sh b/devtools/travis-ci/install.sh index a46f5d92..aec69a4a 100644 --- a/devtools/travis-ci/install.sh +++ b/devtools/travis-ci/install.sh @@ -2,20 +2,18 @@ pushd . cd $HOME -# Install Miniconda -MINICONDA=Miniconda3-latest-Linux-x86_64.sh -MINICONDA_HOME=$HOME/miniconda3 +if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then MINICONDA=Miniconda3-latest-MacOSX-x86_64.sh; fi +if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then MINICONDA=Miniconda3-latest-Linux-x86_64.sh; fi + MINICONDA_MD5=$(curl -s https://repo.continuum.io/miniconda/ | grep -A3 $MINICONDA | sed -n '4p' | sed -n 's/ *\(.*\)<\/td> */\1/p') -wget -q http://repo.continuum.io/miniconda/$MINICONDA -if [[ $MINICONDA_MD5 != $(md5sum $MINICONDA | cut -d ' ' -f 1) ]]; then - echo "Miniconda MD5 mismatch" - exit 1 -fi -bash $MINICONDA -b -p $MINICONDA_HOME +wget https://repo.continuum.io/miniconda/$MINICONDA +bash $MINICONDA -b +rm -f $MINICONDA export PATH=$HOME/miniconda3/bin:$PATH -# Configure miniconda -export PIP_ARGS="-U" -export PATH=$MINICONDA_HOME/bin:$PATH -conda update --yes conda +conda update -yq conda +conda install -yq conda-build jinja2 anaconda-client pip + +# Restore original directory +popd diff --git a/devtools/travis-ci/requirements.txt b/devtools/travis-ci/requirements.txt deleted file mode 100644 index 02dcfad2..00000000 --- a/devtools/travis-ci/requirements.txt +++ /dev/null @@ -1,8 +0,0 @@ -numpy -nose -pandas -pymc -pint -docopt -schema -scikit-learn diff --git a/examples/sampl4_notebook/models.py b/examples/sampl4_notebook/models.py index 523e1c2e..fd9da509 100644 --- a/examples/sampl4_notebook/models.py +++ b/examples/sampl4_notebook/models.py @@ -41,7 +41,7 @@ def __init__(self, dictionary, beta, max_scale=1.03, verbose=0, interval=100): # Store stochastics. self.dictionary = dictionary - + # Initialize superclass. pymc.StepMethod.__init__(self, dictionary.values(), verbose) @@ -60,30 +60,30 @@ def __init__(self, dictionary, beta, max_scale=1.03, verbose=0, interval=100): # Report if self.verbose: - print "Initialization..." - print "max_scale: ", self.max_scale + print("Initialization...") + print("max_scale: ", self.max_scale) def propose(self): # Choose trial scaling factor or its inverse with equal probability, so that proposal move is symmetric. factor = (self.max_scale - 1) * numpy.random.rand() + 1; - if (numpy.random.rand() < 0.5): + if (numpy.random.rand() < 0.5): factor = 1./factor; - # Scale thermodynamic parameters and variables with this factor. + # Scale thermodynamic parameters and variables with this factor. self.dictionary['Ls'].value = self.dictionary['Ls'].value * factor self.dictionary['P0'].value = self.dictionary['P0'].value * factor self.dictionary['DeltaH'].value = self.dictionary['DeltaH'].value / factor self.dictionary['DeltaG'].value = self.dictionary['DeltaG'].value + (1./self.beta) * numpy.log(factor); - return + return - def step(self): + def step(self): # Probability and likelihood for stochastic's current value: logp = sum([stochastic.logp for stochastic in self.stochastics]) loglike = self.loglike if self.verbose > 1: - print 'Current likelihood: ', logp+loglike - + print('Current likelihood: ', logp+loglike) + # Sample a candidate value self.propose() @@ -94,41 +94,41 @@ def step(self): logp_p = sum([stochastic.logp for stochastic in self.stochastics]) loglike_p = self.loglike if self.verbose > 2: - print 'Current likelihood: ', logp+loglike - + print('Current likelihood: ', logp+loglike) + if numpy.log(numpy.random.rand()) < logp_p + loglike_p - logp - loglike: accept = True self.accepted += 1 if self.verbose > 2: - print 'Accepted' + print('Accepted') else: self.rejected += 1 if self.verbose > 2: - print 'Rejected' + print('Rejected') except pymc.ZeroProbability: self.rejected += 1 logp_p = None loglike_p = None if self.verbose > 2: - print 'Rejected with ZeroProbability error.' + print('Rejected with ZeroProbability error.') if (not self._current_iter % self.interval) and self.verbose > 1: - print "Step ", self._current_iter - print "\tLogprobability (current, proposed): ", logp, logp_p - print "\tloglike (current, proposed): : ", loglike, loglike_p + print("Step ", self._current_iter) + print("\tLogprobability (current, proposed): ", logp, logp_p) + print("\tloglike (current, proposed): : ", loglike, loglike_p) for stochastic in self.stochastics: - print "\t", stochastic.__name__, stochastic.last_value, stochastic.value + print("\t", stochastic.__name__, stochastic.last_value, stochastic.value) if accept: - print "\tAccepted\t*******\n" + print("\tAccepted\t*******\n") else: - print "\tRejected\n" - print "\tAcceptance ratio: ", self.accepted/(self.accepted+self.rejected) - + print("\tRejected\n") + print("\tAcceptance ratio: ", self.accepted/(self.accepted+self.rejected)) + if not accept: self.reject() self._current_iter += 1 - + return @classmethod @@ -136,14 +136,14 @@ def competence(self, stochastic): if str(stochastic) in ['DeltaG', 'DeltaH', 'DeltaH_0', 'Ls', 'P0']: return 1 return 0 - + def reject(self): for stochastic in self.stochastics: # stochastic.value = stochastic.last_value stochastic.revert() def tune(self, verbose): - return False + return False #============================================================================================= # Binding models #============================================================================================= @@ -156,14 +156,14 @@ class BindingModel(object): def __init__(self): pass - + #============================================================================================= # Two-component binding model #============================================================================================= class TwoComponentBindingModel(BindingModel): def __init__(self, Ls_stated, P0_stated, q_n_observed, DeltaVn, temperature, V0): - + # Determine number of observations. self.N = q_n_observed.size @@ -178,7 +178,7 @@ def __init__(self, Ls_stated, P0_stated, q_n_observed, DeltaVn, temperature, V0) # Store temperature. self.temperature = temperature # temperature (kelvin) - self.beta = 1.0 / (kB * temperature) # inverse temperature 1/(kcal/mol) + self.beta = 1.0 / (kB * temperature) # inverse temperature 1/(kcal/mol) # Compute uncertainties in stated concentrations. dP0 = 0.10 * P0_stated # uncertainty in protein stated concentration (M) - 10% error @@ -189,7 +189,7 @@ def __init__(self, Ls_stated, P0_stated, q_n_observed, DeltaVn, temperature, V0) DeltaG_guess = -8.3 # kcal/mol DeltaH_guess = -12.0 # kcal/mol DeltaH_0_guess = q_n_observed[-1] # cal/injection - + # Determine min and max range for log_sigma log_sigma_min = log_sigma_guess - 10 log_sigma_max = log_sigma_guess + 5 @@ -202,7 +202,7 @@ def __init__(self, Ls_stated, P0_stated, q_n_observed, DeltaVn, temperature, V0) heat_interval = q_n_observed.max() - q_n_observed.min() DeltaH_0_min = q_n_observed.min() - heat_interval # (cal/mol) DeltaH_0_max = q_n_observed.max() + heat_interval # (cal/mol) - + # Define priors for concentrations. #self.P0 = pymc.Normal('P0', mu=P0_stated, tau=1.0/dP0**2, value=P0_stated) #self.Ls = pymc.Normal('Ls', mu=Ls_stated, tau=1.0/dLs**2, value=Ls_stated) @@ -233,24 +233,24 @@ def __init__(self, Ls_stated, P0_stated, q_n_observed, DeltaVn, temperature, V0) mcmc.use_step_method(RescalingStep, { 'Ls' : self.Ls, 'P0' : self.P0, 'DeltaH' : self.DeltaH, 'DeltaG' : self.DeltaG }, self.beta) self.mcmc = mcmc return - + def expected_injection_heats(self, P0, Ls, DeltaG, DeltaH, DeltaH_0, q_n_obs): """ Expected heats of injection for two-component binding model. - + ARGUMENTS - + DeltaG - free energy of binding (kcal/mol) DeltaH - enthalpy of binding (kcal/mol) DeltaH_0 - heat of injection (cal/mol) - + """ - + debug = False - + Kd = exp(self.beta * DeltaG) * C0 # dissociation constant (M) N = self.N - + # Compute complex concentrations. Pn = numpy.zeros([N], numpy.float64) # Pn[n] is the protein concentration in sample cell after n injections (M) Ln = numpy.zeros([N], numpy.float64) # Ln[n] is the ligand concentration in sample cell after n injections (M) @@ -266,35 +266,38 @@ def expected_injection_heats(self, P0, Ls, DeltaG, DeltaH, DeltaH_0, q_n_obs): PLn[n] = 0.5/self.V0 * ((P + L + Kd*self.V0) - sqrt((P + L + Kd*self.V0)**2 - 4*P*L)); # complex concentration (M) Pn[n] = P/self.V0 - PLn[n]; # free protein concentration in sample cell after n injections (M) Ln[n] = L/self.V0 - PLn[n]; # free ligand concentration in sample cell after n injections (M) - + # Compute expected injection heats. q_n = numpy.zeros([N], numpy.float64) # q_n_model[n] is the expected heat from injection n # Instantaneous injection model (perfusion) q_n[0] = -(1000.0*DeltaH) * self.V0 * PLn[0] - DeltaH_0 # first injection for n in range(1,N): - d = 1.0 - (self.DeltaVn[n] / self.V0) # dilution factor (dimensionless) + d = 1.0 - (self.DeltaVn[n] / self.V0) # dilution factor (dimensionless) q_n[n] = -(1000.0*DeltaH) * self.V0 * (PLn[n] - d*PLn[n-1]) - DeltaH_0 # subsequent injections # Debug output if debug: - print "DeltaG = %6.1f kcal/mol ; DeltaH = %6.1f kcal/mol ; DeltaH_0 = %6.1f ucal/injection" % (DeltaG, DeltaH, DeltaH_0*1e6) + print("DeltaG = %6.1f kcal/mol ; DeltaH = %6.1f kcal/mol ; DeltaH_0 = %6.1f ucal/injection" % (DeltaG, DeltaH, DeltaH_0*1e6)) + line = "" for n in range(N): - print "%6.1f" % (PLn[n]*1e6), - print "" + line += "%6.1f" % (PLn[n]*1e6) + print(line) + line = "" for n in range(N): - print "%6.1f" % (q_n[n]*1e6), - print "" + line += "%6.1f" % (q_n[n]*1e6) + print(line) + line = "" for n in range(N): - print "%6.1f" % (q_n_obs[n]*1e6), - print "" - print "" - + line += "%6.1f" % (q_n_obs[n]*1e6) + print(line) + print('') + return q_n def tau(self, log_sigma): """ Injection heat measurement precision. - + """ return exp(-2.0*log_sigma) @@ -307,13 +310,13 @@ class Experiment(object): A calorimetry experiment. """ - + def __init__(self, sample_cell_concentrations, syringe_concentrations, injection_volumes, injection_heats, temperature): """ Initialize a calorimetry experiment. - + ARGUMENTS - + sample_cell_concentrations (dict) - a dictionary of initial concentrations of each species in sample cell (M) syringe_concentrations (dict) - a dictionary of initial concentrations of each species in the syringe injection_volumes (list or numpy array of N floats) - injection volumes in L @@ -336,7 +339,7 @@ def __init__(self, sample_cell_concentrations, syringe_concentrations, injection """ # TODO: Do sanity checking to make sure number of injections matches up, etc. - + self.sample_cell_concentrations = sample_cell_concentrations self.syringe_concentrations = syringe_concentrations self.injection_volumes = numpy.array(injection_volumes) @@ -352,9 +355,9 @@ def __init__(self, sample_cell_concentrations, syringe_concentrations, injection class CompetitiveBindingModel(BindingModel): """ Competitive binding model. - + """ - + def __init__(self, experiments, receptor, V0, concentration_uncertainty=0.10, verbose=False): """ ARGUMENTS @@ -364,7 +367,7 @@ def __init__(self, experiments, receptor, V0, concentration_uncertainty=0.10, ve V0 (float) - calorimeter sample cell volume OPTIONAL ARGUMENTS - + concentration_uncertainty (float) - relative uncertainty in concentrations """ @@ -374,30 +377,30 @@ def __init__(self, experiments, receptor, V0, concentration_uncertainty=0.10, ve # Store temperature. # NOTE: Right now, there can only be one. self.temperature = experiments[0].temperature # temperature (kelvin) - self.beta = 1.0 / (kB * self.temperature) # inverse temperature 1/(kcal/mol) - + self.beta = 1.0 / (kB * self.temperature) # inverse temperature 1/(kcal/mol) + # Store copy of experiments. self.experiments = copy.deepcopy(experiments) - if verbose: print "%d experiments" % len(self.experiments) + if verbose: print("%d experiments" % len(self.experiments)) # Store sample cell volume. self.V0 = V0 - + # Store the name of the receptor. self.receptor = receptor - if verbose: print "species '%s' will be treated as receptor" % self.receptor + if verbose: print("species '%s' will be treated as receptor" % self.receptor) # Make a list of names of all molecular species. self.species = set() # all molecular species for experiment in experiments: self.species.update( experiment.sample_cell_concentrations.keys() ) self.species.update( experiment.syringe_concentrations.keys() ) - if verbose: print "species: ", self.species - + if verbose: print("species: ", self.species) + # Make a list of all ligands. self.ligands = copy.deepcopy(self.species) self.ligands.remove(receptor) - if verbose: print "ligands: ", self.ligands + if verbose: print("ligands: ", self.ligands) # Create a list of all stochastics. self.stochastics = list() @@ -416,10 +419,10 @@ def __init__(self, experiments, receptor, V0, concentration_uncertainty=0.10, ve name = "DeltaH of %s * %s" % (self.receptor, ligand) x = pymc.Uniform(name, lower=DeltaH_min, upper=DeltaH_max, value=0.0) self.thermodynamic_parameters[name] = x - self.stochastics.append(x) + self.stochastics.append(x) if verbose: - print "thermodynamic parameters:" - print self.thermodynamic_parameters + print("thermodynamic parameters:") + print(self.thermodynamic_parameters) # DEBUG: Set initial thermodynamic parameters to literature values. self.thermodynamic_parameters["DeltaG of HIV protease * acetyl pepstatin"].value = -9.0 @@ -428,7 +431,7 @@ def __init__(self, experiments, receptor, V0, concentration_uncertainty=0.10, ve self.thermodynamic_parameters["DeltaH of HIV protease * KNI-10033"].value = -8.200 self.thermodynamic_parameters["DeltaG of HIV protease * KNI-10075"].value = -14.620 self.thermodynamic_parameters["DeltaH of HIV protease * KNI-10075"].value = -12.120 - + # Determine min and max range for log_sigma (log of instrument heat measurement error) # TODO: This should depend on a number of factors, like integration time, heat signal, etc.? sigma_guess = 0.0 @@ -447,7 +450,7 @@ def __init__(self, experiments, receptor, V0, concentration_uncertainty=0.10, ve for (index, experiment) in enumerate(self.experiments): # Number of observations experiment.ninjections = experiment.observed_injection_heats.size - if verbose: print "Experiment %d has %d injections" % (index, experiment.ninjections) + if verbose: print("Experiment %d has %d injections" % (index, experiment.ninjections)) # Heat of dilution / mixing # We allow the heat of dilution/mixing to range in observed range of heats, plus a larger margin of the range of oberved heats. @@ -466,7 +469,7 @@ def __init__(self, experiments, receptor, V0, concentration_uncertainty=0.10, ve value=concentration) experiment.true_sample_cell_concentrations[species] = x self.stochastics.append(x) - + experiment.true_syringe_concentrations = dict() for species, concentration in experiment.syringe_concentrations.iteritems(): x = pymc.Lognormal("initial syringe concentration of %s in experiment %d" % (species, index), @@ -474,14 +477,14 @@ def __init__(self, experiments, receptor, V0, concentration_uncertainty=0.10, ve value=concentration) experiment.true_syringe_concentrations[species] = x self.stochastics.append(x) - + # Add species not explicitly listed with zero concentration. for species in self.species: if species not in experiment.true_sample_cell_concentrations: experiment.true_sample_cell_concentrations[species] = 0.0 if species not in experiment.true_syringe_concentrations: experiment.true_syringe_concentrations[species] = 0.0 - + # True injection heats experiment.true_injection_heats = pymc.Lambda("true injection heats for experiment %d" % index, lambda experiment=experiment, @@ -499,12 +502,12 @@ def __init__(self, experiments, receptor, V0, concentration_uncertainty=0.10, ve self.stochastics.append(experiment.observation) # Create sampler. - print "Creating sampler..." + print("Creating sampler...") mcmc = pymc.MCMC(self.stochastics, db='ram') #db = pymc.database.pickle.load('MCMC.pickle') # DEBUG #mcmc = pymc.MCMC(self.stochastics, db=db) for stochastic in self.stochastics: - print stochastic + print(stochastic) try: mcmc.use_step_method(pymc.Metropolis, stochastic) except: @@ -517,12 +520,12 @@ def __init__(self, experiments, receptor, V0, concentration_uncertainty=0.10, ve mcmc.use_step_method(RescalingStep, { 'Ls' : self.experiments[1].true_syringe_concentrations['KNI-10033'], 'P0' : self.experiments[1].true_sample_cell_concentrations['HIV protease'], 'DeltaH' : self.thermodynamic_parameters['DeltaH of HIV protease * KNI-10033'], - 'DeltaG' : self.thermodynamic_parameters['DeltaG of HIV protease * KNI-10033'] }, self.beta) + 'DeltaG' : self.thermodynamic_parameters['DeltaG of HIV protease * KNI-10033'] }, self.beta) mcmc.use_step_method(RescalingStep, { 'Ls' : self.experiments[2].true_syringe_concentrations['KNI-10075'], 'P0' : self.experiments[2].true_sample_cell_concentrations['HIV protease'], 'DeltaH' : self.thermodynamic_parameters['DeltaH of HIV protease * KNI-10075'], - 'DeltaG' : self.thermodynamic_parameters['DeltaG of HIV protease * KNI-10075'] }, self.beta) + 'DeltaG' : self.thermodynamic_parameters['DeltaG of HIV protease * KNI-10075'] }, self.beta) self.mcmc = mcmc @@ -536,7 +539,7 @@ def equilibrium_concentrations(self, Ka_n, C0_R, C0_Ln, V, c0=None): x_R (float) - the total number of moles of receptor in the sample volume x_n (numpy N-array of float) - x_n[n] is the total number of moles of ligand species n in the sample volume V (float) - the total sample volume (L) - + RETURNS C_n (numpy N-array of float) - C_n[n] is the concentration of complex of receptor with ligand species n @@ -550,37 +553,37 @@ def equilibrium_concentrations(self, Ka_n, C0_R, C0_Ln, V, c0=None): >>> C_PLn = equilibrium_concentrations(Ka_n, x_R, x_Ln, V) NOTES - + Each complex concentration C_n must obey the relation - + Ka_n[n] = C_RLn[n] / (C_R * C_Ln[n]) for n = 1..N - + with conservation of mass constraints - + V * (C_Ln[n] + C_RLn[n]) = x_Ln[n] for n = 1..N - + and - + V * (C_R + C_RLn[:].sum()) = x_R - + along with the constraints - + 0 <= V * C_RLn[n] <= min(x_Ln[n], x_R) for n = 1..N V * C_RLn[:].sum() <= x_R - + We can rearrange these expressions to give - + V * C_R * C_Ln[n] * Ka_n[n] - V * C_RLn[n] = 0 - + and eliminate C_Ln[n] and C_R to give - + V * (x_R/V - C_RLn[:].sum()) * (x_Ln[n]/V - C_RLn[n]) * Ka_n[n] - V * C_RLn[n] = 0 for n = 1..N - + """ x_R = C0_R * V x_Ln = C0_Ln * V - + nspecies = Ka_n.size #print "x_R = ", x_R #print "x_Ln = ", x_Ln @@ -592,7 +595,7 @@ def func(C_RLn): f_n = V * (x_R/V - C_RLn[:].sum()) * (x_Ln[:]/V - C_RLn[:]) * Ka_n[:] - V * C_RLn[:] #print "f_n = ", f_n return f_n - + def fprime(C_RLn): nspecies = C_RLn.size G_nm = numpy.zeros([nspecies,nspecies], numpy.float64) # G_nm[n,m] is the derivative of func[n] with respect to C_RLn[m] @@ -606,12 +609,12 @@ def sfunc(s): f_n = V * (x_R/V - (s[:]**2).sum()) * (x_Ln[:]/V - s[:]**2) * Ka_n[:] - V * s[:]**2 #print "f_n = ", f_n return f_n - + def sfprime(s): nspecies = s.size G_nm = numpy.zeros([nspecies,nspecies], numpy.float64) # G_nm[n,m] is the derivative of func[n] with respect to C_RLn[m] for n in range(nspecies): - G_nm[n,:] = - V * (x_Ln[:]/V - s[:]**2) * Ka_n[:] + G_nm[n,:] = - V * (x_Ln[:]/V - s[:]**2) * Ka_n[:] G_nm[n,n] -= V * (Ka_n[n] * (x_R/V - (s[:]**2).sum()) + 1.0) G_nm[n,:] *= 2. * s[n] return G_nm @@ -622,9 +625,9 @@ def sfprime(s): #x0 = (x_Ln / V).copy() #x = scipy.optimize.fsolve(func, x0, fprime=fprime) #C_RLn = x - + #x0 = numpy.sqrt(x_Ln / V).copy() - #x = scipy.optimize.fsolve(sfunc, x0, fprime=sfprime) + #x = scipy.optimize.fsolve(sfunc, x0, fprime=sfprime) #C_RLn = x**2 def objective(x): @@ -637,8 +640,8 @@ def objective(x): grad += 2 * f_n[n] * G_nm[n,:] return (obj, grad) - - #x0 = numpy.zeros([nspecies], numpy.float64) + + #x0 = numpy.zeros([nspecies], numpy.float64) #bounds = list() #for n in range(nspecies): # m = min(C0_R, C0_Ln[n]) @@ -649,7 +652,7 @@ def objective(x): def ode(c_n, t, Ka_n, x_Ln, x_R): dc_n = - c_n[:] + Ka_n[:] * (x_Ln[:]/V - c_n[:]) * (x_R/V - c_n[:].sum()) return dc_n - + def odegrad(c_n, t, Ka_n, x_Ln, x_R): N = c_n.size d2c = numpy.zeros([N,N], numpy.float64) @@ -684,9 +687,9 @@ def odegrad(c_n, t, Ka_n, x_Ln, x_R): #c[indices] = scipy.optimize.fsolve(ode, c[indices], fprime=odegrad, args=(0.0, Ka_n[indices], x_Ln[indices], x_R), xtol=1.0e-6, warning=False) c[indices] = scipy.optimize.fsolve(ode, c[indices], fprime=odegrad, args=(0.0, Ka_n[indices], x_Ln[indices], x_R), xtol=1.0e-6) C_RLn = c - + return C_RLn - + def expected_injection_heats(self, experiment, true_sample_cell_concentrations, true_syringe_concentrations, DeltaH_0, thermodynamic_parameters): """ Expected heats of injection for two-component binding model. @@ -694,24 +697,24 @@ def expected_injection_heats(self, experiment, true_sample_cell_concentrations, TODO - Make experiment a dict, or somehow tell it how to replace members of 'experiment'? - + ARGUMENTS - + sample_cell_concentrations (dict of floats) - concentrations[species] is the initial concentration of species in sample cell, or zero if absent (M) syringe_concentrations (dict of floats) - concentrations[species] is the initial concentration of species in sample cell, or zero if absent (M) thermodynamic_parameters (dict of floats) - thermodynamic_parameters[parameter] is the value of thermodynamic parameter (kcal/mol) e.g. for parameter 'DeltaG of receptor * species' V_n (numpy array of floats) - V_n[n] is injection volume of injection n (L) - + """ - + debug = False # Number of ligand species nspecies = len(self.ligands) # Compute association constants for receptor and each ligand species. - DeltaG_n = numpy.zeros([nspecies], numpy.float64) # + DeltaG_n = numpy.zeros([nspecies], numpy.float64) # for (n, ligand) in enumerate(self.ligands): name = "DeltaG of %s * %s" % (self.receptor, ligand) # determine name of free energy of binding for this ligand DeltaG_n[n] = thermodynamic_parameters[name] # retrieve free energy of binding @@ -735,18 +738,18 @@ def expected_injection_heats(self, experiment, true_sample_cell_concentrations, #print "x_R in mol:" #print x_Ri #print "x_Lin in mol: " - #print x_Lin + #print x_Lin # Solve for initial concentration. - x_R0 = true_sample_cell_concentrations[self.receptor] + x_R0 = true_sample_cell_concentrations[self.receptor] x_L0n = numpy.zeros([nspecies], numpy.float64) - C_RL0n = numpy.zeros([nspecies], numpy.float64) + C_RL0n = numpy.zeros([nspecies], numpy.float64) for (n, ligand) in enumerate(self.ligands): x_L0n[n] = true_sample_cell_concentrations[ligand] C_RL0n[:] = self.equilibrium_concentrations(Ka_n, x_R0, x_L0n[:], self.V0) #print "C_RL0n in uM:" #print C_RL0n * 1.e6 - + # Compute complex concentrations after each injection. # NOTE: The total cell volume would be modified for a cumulative model. C_RLin = numpy.zeros([experiment.ninjections,nspecies], numpy.float64) # C_RLin[i,n] is the concentration of complex RLn[n] after injection i @@ -765,29 +768,29 @@ def expected_injection_heats(self, experiment, true_sample_cell_concentrations, # NOTE: This is for an instantaneous injection / perfusion model. q_n = DeltaH_0 * numpy.ones([experiment.ninjections], numpy.float64) # q_n_model[n] is the expected heat from injection n d = 1.0 - (experiment.injection_volumes[0] / self.V0) # dilution factor (dimensionless) - for n in range(nspecies): + for n in range(nspecies): q_n[0] += (1000.0*DeltaH_n[n]) * V0 * (C_RLin[0,n] - d*C_RL0n[n]) # first injection for i in range(1,experiment.ninjections): - d = 1.0 - (experiment.injection_volumes[i] / self.V0) # dilution factor (dimensionless) + d = 1.0 - (experiment.injection_volumes[i] / self.V0) # dilution factor (dimensionless) for n in range(nspecies): q_n[i] += (1000.0*DeltaH_n[n]) * V0 * (C_RLin[i,n] - d*C_RLin[i-1,n]) # subsequent injections # Debug output debug = False if debug: - print experiment.name - print "DeltaG = ", DeltaG_n - print "DeltaH = ", DeltaH_n - print "DeltaH_0 = ", DeltaH_0 - print "model: ", + print(experiment.name) + print("DeltaG = ", DeltaG_n) + print("DeltaH = ", DeltaH_n) + print("DeltaH_0 = ", DeltaH_0) + print("model: ") for heat in q_n: - print "%6.1f" % (heat*1e6), - print "" - print "obs : ", + print("%6.1f" % (heat*1e6)) + print("") + print("obs : ") for heat in experiment.observed_injection_heats: - print "%6.1f" % (heat*1e6), - print "" - print "" + print("%6.1f" % (heat*1e6)) + print("") + print("") return q_n @@ -811,22 +814,22 @@ def expected_injection_heats(self, experiment, true_sample_cell_concentrations, Ls_stated = 384.e-6 # ligand syringe stated concentration (M) temperature = 298.15 # temperature (K) q_n = numpy.array([ - -13.343, -13.417, -13.279, -13.199, -13.118, -12.781, -12.600, -12.124, -11.633, -10.921, -10.009, -8.810, + -13.343, -13.417, -13.279, -13.199, -13.118, -12.781, -12.600, -12.124, -11.633, -10.921, -10.009, -8.810, -7.661, -6.272, -5.163, -4.228, -3.519, -3.055, -2.599, -2.512, -2.197, -2.096, -2.087, -1.959, -1.776, -1.879, -1.894, -1.813, -1.740, -1.810]) # integrated heats of injection (kcal/mol injectant) q_n = q_n * DeltaV * Ls_stated * 1000.0 # convert injection heats to cal/injection beta = 1.0 / (kB * temperature) # inverse temperature 1/(kcal/mol) - + #============================================================================================= # Erensto Freire data for HIV protease inhibitors KNI-10033 and KNI-10075 #============================================================================================= experiments = list() - + # # acetyl pepstatin # - + V0 = 1.4301e-3 # volume of calorimeter sample cell listed in manual (L) V0 = V0 - 0.044e-3; # Tellinghuisen volume correction for VP-ITC (L) sample_cell_concentrations = {'HIV protease' : 20.e-6} @@ -842,7 +845,7 @@ def expected_injection_heats(self, experiment, true_sample_cell_concentrations, experiment.name = "acetyl pepstatin binding to HIV protease" experiment.reference = "Nature Protocols 1:186, 2006; Fig. 1, left panel" experiments.append(experiment) - + # # KNI-10033 # @@ -858,11 +861,11 @@ def expected_injection_heats(self, experiment, true_sample_cell_concentrations, experiment = Experiment(sample_cell_concentrations, syringe_concentrations, injection_volumes, injection_heats, temperature) experiment.name = "KNI-10033 displacement of acetyl pepstatin binding to HIV protease" experiments.append(experiment) - + # # KNI-10075 # - + sample_cell_concentrations = {'HIV protease' : 8.8e-6, 'acetyl pepstatin' : 510.e-6} # initial sample cell concentrations (M) syringe_concentrations = {'KNI-10075' : 55.e-6} Ls_stated = 55.e-6 # KNI-10033 syringe concentration (M) @@ -873,19 +876,17 @@ def expected_injection_heats(self, experiment, true_sample_cell_concentrations, experiment = Experiment(sample_cell_concentrations, syringe_concentrations, injection_volumes, injection_heats, temperature) experiment.name = "KNI-10075 displacement of acetyl pepstatin binding to HIV protease" experiments.append(experiment) - + #============================================================================================= # MCMC inference #============================================================================================= - #model = TwoComponentBindingModel(Ls_stated, P0_stated, q_n, DeltaV, temperature, V0) + #model = TwoComponentBindingModel(Ls_stated, P0_stated, q_n, DeltaV, temperature, V0) model = CompetitiveBindingModel(experiments, 'HIV protease', V0, verbose=True) - + niters = 10000 # number of iterations nburn = 1000 # number of burn-in iterations nthin = 1 # thinning period - + model.mcmc.sample(iter=niters, burn=nburn, thin=nthin, progress_bar=True) pymc.Matplot.plot(model.mcmc) - - diff --git a/examples/two-component-binding/acetyl-pepstatin/acetyl-pepstatin.py b/examples/two-component-binding/acetyl-pepstatin/acetyl-pepstatin.py index 99e6314b..841bb717 100644 --- a/examples/two-component-binding/acetyl-pepstatin/acetyl-pepstatin.py +++ b/examples/two-component-binding/acetyl-pepstatin/acetyl-pepstatin.py @@ -19,20 +19,19 @@ run = Run() run.title = "Direct titration of acetyl pepstatin" -run.syringe.append('acetyl pepstatin' : 200 * Units.uM) -run.cell.append('HIV-1 protease' : 20 * Units.uM) +run.syringe['acetyl pepstatin'] = 200 * Units.uM +run.cell['HIV-1 protease'] = 20 * Units.uM run.injection_volume = 10 * Units.uL run.heat_per_injection = [7.347, 7.219, 7.024, 6.537, 5.470, 3.911, 2.258, 1.231, 0.744, 0.497, 0.302, 0.160, 0.087, 0.062, 0.062, 0.040, 0.017, -0.013, -0.013] experiment.runs.append(run) run = Run() run.title = "Competition of acetyl pepstatin by KNI-764" -run.syringe.append('KNI-764' : 250 * Units.uM) -run.cell.append('HIV-1 protease' : 20 * Units.uM) -run.cell.append('acetyl pepstatin' : 200 * Units.uM) +run.syringe['KNI-764'] = 250 * Units.uM +run.cell['HIV-1 protease'] = 20 * Units.uM +run.cell['acetyl pepstatin'] = 200 * Units.uM run.injection_volume = 10 * Units.uL experiment.runs.append(run) # Run the analysis. analyze(experiment) - diff --git a/setup.py b/setup.py index fc31d087..203d834c 100644 --- a/setup.py +++ b/setup.py @@ -7,6 +7,8 @@ import subprocess from setuptools import setup +from basesetup import (find_packages, write_version_py, build_ext, + StaticLibrary, CompilerDetection) DOCLINES = __doc__.split("\n") @@ -117,8 +119,8 @@ def find_package_data(data_root, package_root): url='https://github.com/choderalab/bayesian-itc', platforms=['Linux', 'Mac OS-X', 'Unix'], classifiers=CLASSIFIERS.splitlines(), - package_dir={'bitc': 'bitc'}, - packages=['bitc'], + package_dir={'bitc' : 'bitc', 'bitc.scripts' : 'scripts'}, + packages=find_packages(), package_data={'bitc': find_package_data('examples', 'bitc')}, # NOTE: examples installs to bitc.egg/examples/, NOT bitc.egg/bitc/examples/. You need to do utils.get_data_filename("../examples/*/setup/"). zip_safe=False, install_requires=[