From 15de286b85ef4af9ac2de7f074c5971f47a52fd8 Mon Sep 17 00:00:00 2001 From: Aaron Stern Date: Mon, 30 Nov 2020 07:12:47 -0800 Subject: [PATCH 01/10] Update inference.py --- inference.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inference.py b/inference.py index 27862d9..7e52801 100644 --- a/inference.py +++ b/inference.py @@ -350,7 +350,7 @@ def traj_wrapper(theta,timeBins,N,freqs,z_bins,z_logcdf,z_logsf,ancGLs,ancHapGLs print('========') print('epoch\tselection') for s,t,u in zip(S,timeBins[:-1],timeBins[1:]): - print('%d-%d\t%.5f'%(t,u,s)) + print('%d-%d\t%.7f'%(t,u,s)) # infer trajectory @ MLE of selection parameter From d9d93e92bc153be025084a5bfca577d83a70ce5f Mon Sep 17 00:00:00 2001 From: Aaron Stern Date: Mon, 25 Jan 2021 11:37:23 -0800 Subject: [PATCH 02/10] likelihood function approximation for PALM --- inference.py | 44 +++++++++++++++++++++++++++++++++----------- 1 file changed, 33 insertions(+), 11 deletions(-) diff --git a/inference.py b/inference.py index 7e52801..283c0c6 100644 --- a/inference.py +++ b/inference.py @@ -90,9 +90,12 @@ def parse_args(): parser.add_argument('--tSkip',type=int,default=1) parser.add_argument('--df',type=int,default=150) parser.add_argument('--betaParam',type=float,default=0.5) + parser.add_argument('--lik',action='store_true',help='saves likelihood function, used by PALM. **Only compatible when a single selection coefficient is estimated.**') + parser.add_argument('--w',type=float,default=0.01,help='width used to estimate likelihood function (--lik)') return parser.parse_args() + def load_normal_tables(): # read in global Phi(z) lookups z_bins = np.genfromtxt('utils/z_bins.txt') @@ -214,7 +217,6 @@ def likelihood_wrapper(theta,timeBins,N,freqs,z_bins,z_logcdf,z_logsf,ancGLs,anc return np.inf sel = Sprime[np.digitize(epochs,timeBins,right=False)-1] - tShape = times.shape if tShape[2] == 0: t = np.zeros((2,0)) @@ -308,6 +310,10 @@ def traj_wrapper(theta,timeBins,N,freqs,z_bins,z_logcdf,z_logsf,ancGLs,ancHapGLs S0 = 0.0 * np.ones(T-1) opts = {'xatol':1e-4} + if args.lik and T > 2: + print('ERROR: Option --lik incompatible with estimating >1 selection coefficient.') + raise ValueError + if T == 2: Simplex = np.reshape(np.array([-0.05,0.05]),(2,1)) elif T > 2: @@ -323,28 +329,44 @@ def traj_wrapper(theta,timeBins,N,freqs,z_bins,z_logcdf,z_logsf,ancGLs,ancHapGLs opts['initial_simplex']=Simplex #for tup in product(*[[-1,1] for i in range(3)]): - logL0 = likelihood_wrapper(S0,timeBins,Ne,freqs,z_bins,z_logcdf,z_logsf,ancientGLs,ancientHapGLs,epochs,noCoals,currFreq,h,sMax,changePts) + logL0 = -likelihood_wrapper(S0,timeBins,Ne,freqs,z_bins,z_logcdf,z_logsf,ancientGLs,ancientHapGLs,epochs,noCoals,currFreq,h,sMax,changePts) print('Optimizing likelihood surface using Nelder-Mead...') if times.shape[2] > 1: print('\t(Importance sampling with M = %d Relate samples)'%(times.shape[2])) print() minargs = (timeBins,Ne,freqs,z_bins,z_logcdf,z_logsf,ancientGLs,ancientHapGLs,epochs,noCoals,currFreq,h,sMax,changePts) + + res = minimize(likelihood_wrapper, - S0, - args=minargs, - options=opts, - #bounds=bounds, - method='Nelder-Mead') + S0, + args=minargs, + options=opts, + #bounds=bounds, + method='Nelder-Mead') S = res.x L = res.fun - #Hinv = np.linalg.inv(res.hess) - #se = np.sqrt(np.diag(Hinv)) + + if args.lik: + print('Fitting likelihood function...') + print() + Svec = np.linspace(S[0]-args.w,S[0]+args.w,20) + Lvec = [] + for s in Svec: + l = likelihood_wrapper(np.array([s]),timeBins,Ne,freqs,z_bins,z_logcdf,z_logsf,ancientGLs,ancientHapGLs,epochs,noCoals,currFreq,h,sMax,changePts) + Lvec.append(l) + Lvec = np.array(Lvec) + S = np.array([Svec[np.argmax(-Lvec)]]) + L = np.max(-Lvec) + p = np.polyfit(Svec,Lvec,deg=2) + if args.out != None: + np.save(args.out+'.quad_fit.npy',p) + print('#'*10) print() - print('logLR: %.4f'%(-res.fun+logL0)) + print('logLR: %.4f'%(L-logL0)) print() print('MLE:') print('========') @@ -356,7 +378,7 @@ def traj_wrapper(theta,timeBins,N,freqs,z_bins,z_logcdf,z_logsf,ancGLs,ancHapGLs # infer trajectory @ MLE of selection parameter print(noCoals) - post = traj_wrapper(res.x,timeBins,Ne,freqs,z_bins,z_logcdf,z_logsf,ancientGLs,ancientHapGLs,epochs,noCoals,currFreq,h,sMax,changePts) + post = traj_wrapper(S,timeBins,Ne,freqs,z_bins,z_logcdf,z_logsf,ancientGLs,ancientHapGLs,epochs,noCoals,currFreq,h,sMax,changePts) if args.out != None: out(args,epochs,freqs,post) From b25597811a78ca0cefdde2d6aadb61838e019922 Mon Sep 17 00:00:00 2001 From: Aaron Stern Date: Mon, 25 Jan 2021 12:01:08 -0800 Subject: [PATCH 03/10] sign error --- inference.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/inference.py b/inference.py index 283c0c6..7fc9df6 100644 --- a/inference.py +++ b/inference.py @@ -359,8 +359,9 @@ def traj_wrapper(theta,timeBins,N,freqs,z_bins,z_logcdf,z_logsf,ancGLs,ancHapGLs Lvec = np.array(Lvec) S = np.array([Svec[np.argmax(-Lvec)]]) L = np.max(-Lvec) - p = np.polyfit(Svec,Lvec,deg=2) + p = np.polyfit(Svec,-Lvec,deg=2) if args.out != None: + print(p) np.save(args.out+'.quad_fit.npy',p) From e6baec08066e7d4fe16403892b1a4669940e1be8 Mon Sep 17 00:00:00 2001 From: Aaron Stern Date: Mon, 25 Jan 2021 12:02:05 -0800 Subject: [PATCH 04/10] sign error --- inference.py | 1 - 1 file changed, 1 deletion(-) diff --git a/inference.py b/inference.py index 7fc9df6..ca72fa3 100644 --- a/inference.py +++ b/inference.py @@ -361,7 +361,6 @@ def traj_wrapper(theta,timeBins,N,freqs,z_bins,z_logcdf,z_logsf,ancGLs,ancHapGLs L = np.max(-Lvec) p = np.polyfit(Svec,-Lvec,deg=2) if args.out != None: - print(p) np.save(args.out+'.quad_fit.npy',p) From 20b012db5a35c1eac8939eaa3538b1577bb6af34 Mon Sep 17 00:00:00 2001 From: Aaron Stern Date: Tue, 1 Mar 2022 20:16:02 -0800 Subject: [PATCH 05/10] default beta setting --- inference.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inference.py b/inference.py index ca72fa3..abea7e5 100644 --- a/inference.py +++ b/inference.py @@ -89,7 +89,7 @@ def parse_args(): parser.add_argument('--sMax',type=float,default=0.1) parser.add_argument('--tSkip',type=int,default=1) parser.add_argument('--df',type=int,default=150) - parser.add_argument('--betaParam',type=float,default=0.5) + parser.add_argument('--betaParam',type=float,default=1.0) parser.add_argument('--lik',action='store_true',help='saves likelihood function, used by PALM. **Only compatible when a single selection coefficient is estimated.**') parser.add_argument('--w',type=float,default=0.01,help='width used to estimate likelihood function (--lik)') return parser.parse_args() From 8a6af17c256c3d7a149a083c5136f225484164e3 Mon Sep 17 00:00:00 2001 From: Aaron Stern Date: Tue, 1 Mar 2022 20:39:11 -0800 Subject: [PATCH 06/10] testing --- inference.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/inference.py b/inference.py index b118b60..abea7e5 100644 --- a/inference.py +++ b/inference.py @@ -90,6 +90,8 @@ def parse_args(): parser.add_argument('--tSkip',type=int,default=1) parser.add_argument('--df',type=int,default=150) parser.add_argument('--betaParam',type=float,default=1.0) + parser.add_argument('--lik',action='store_true',help='saves likelihood function, used by PALM. **Only compatible when a single selection coefficient is estimated.**') + parser.add_argument('--w',type=float,default=0.01,help='width used to estimate likelihood function (--lik)') return parser.parse_args() From cfdb4ea0c3e3c9c2ad0133b82a77934189a7b208 Mon Sep 17 00:00:00 2001 From: Evan Irving-Pease Date: Fri, 8 Jul 2022 14:02:15 +0100 Subject: [PATCH 07/10] added `environment.yaml` to create a minimal conda environment for running CLUES --- environment.yaml | 9 +++++++++ 1 file changed, 9 insertions(+) create mode 100644 environment.yaml diff --git a/environment.yaml b/environment.yaml new file mode 100644 index 0000000..73bd326 --- /dev/null +++ b/environment.yaml @@ -0,0 +1,9 @@ +channels: + - conda-forge +dependencies: + - conda-forge::python=3.7.6 + - conda-forge::biopython=1.76 + - conda-forge::matplotlib=3.1.2 + - conda-forge::numba=0.47.0 + - conda-forge::numpy=1.17.5 + - conda-forge::scipy=1.4.1 \ No newline at end of file From b45080d7cc721a8404dd73eca398713fb5c1a2a2 Mon Sep 17 00:00:00 2001 From: Evan Irving-Pease Date: Fri, 8 Jul 2022 14:03:01 +0100 Subject: [PATCH 08/10] bug fix for directory path issue when CLUES is run from a different folder --- inference.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/inference.py b/inference.py index abea7e5..7096bac 100644 --- a/inference.py +++ b/inference.py @@ -97,10 +97,12 @@ def parse_args(): def load_normal_tables(): + import os + dname = os.path.dirname(os.path.abspath(__file__)) # read in global Phi(z) lookups - z_bins = np.genfromtxt('utils/z_bins.txt') - z_logcdf = np.genfromtxt('utils/z_logcdf.txt') - z_logsf = np.genfromtxt('utils/z_logsf.txt') + z_bins = np.genfromtxt(os.path.join(dname, 'utils/z_bins.txt')) + z_logcdf = np.genfromtxt(os.path.join(dname, 'utils/z_logcdf.txt')) + z_logsf = np.genfromtxt(os.path.join(dname, 'utils/z_logsf.txt')) return z_bins,z_logcdf,z_logsf def load_times(args): From 5d0fe03b2592b40bbe9619923e52a1f801cfa209 Mon Sep 17 00:00:00 2001 From: Evan Irving-Pease Date: Fri, 8 Jul 2022 14:07:38 +0100 Subject: [PATCH 09/10] bug fix for `Inf` log-likelihood when the optimal `s` param is close to `sMax` --- inference.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inference.py b/inference.py index 7096bac..f7bcefc 100644 --- a/inference.py +++ b/inference.py @@ -353,7 +353,7 @@ def traj_wrapper(theta,timeBins,N,freqs,z_bins,z_logcdf,z_logsf,ancGLs,ancHapGLs if args.lik: print('Fitting likelihood function...') print() - Svec = np.linspace(S[0]-args.w,S[0]+args.w,20) + Svec = np.linspace(S[0]-args.w,min(S[0]+args.w,sMax),20) Lvec = [] for s in Svec: l = likelihood_wrapper(np.array([s]),timeBins,Ne,freqs,z_bins,z_logcdf,z_logsf,ancientGLs,ancientHapGLs,epochs,noCoals,currFreq,h,sMax,changePts) From d8ec8b4f687867eaf266bf53e744dc669e7b0bc8 Mon Sep 17 00:00:00 2001 From: Evan Irving-Pease Date: Fri, 8 Jul 2022 15:46:57 +0100 Subject: [PATCH 10/10] bug fix for the opposite condition, when `s` is very close to `-sMax` --- inference.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inference.py b/inference.py index f7bcefc..e1d60b6 100644 --- a/inference.py +++ b/inference.py @@ -353,7 +353,7 @@ def traj_wrapper(theta,timeBins,N,freqs,z_bins,z_logcdf,z_logsf,ancGLs,ancHapGLs if args.lik: print('Fitting likelihood function...') print() - Svec = np.linspace(S[0]-args.w,min(S[0]+args.w,sMax),20) + Svec = np.linspace(max(S[0]-args.w,-sMax),min(S[0]+args.w,sMax),20) Lvec = [] for s in Svec: l = likelihood_wrapper(np.array([s]),timeBins,Ne,freqs,z_bins,z_logcdf,z_logsf,ancientGLs,ancientHapGLs,epochs,noCoals,currFreq,h,sMax,changePts)