[med-svn] [Git][med-team/mptp][master] 5 commits: New upstream version 0.2.5
Étienne Mollier (@emollier)
gitlab at salsa.debian.org
Sun Aug 25 18:24:03 BST 2024
Étienne Mollier pushed to branch master at Debian Med / mptp
Commits:
605a8cc7 by Étienne Mollier at 2024-08-25T19:21:16+02:00
New upstream version 0.2.5
- - - - -
4fed58c1 by Étienne Mollier at 2024-08-25T19:21:16+02:00
Update upstream source from tag 'upstream/0.2.5'
Update to upstream version '0.2.5'
with Debian dir 1ed2b46f8134ce7f97c9a22ada15de4ca56d23d7
- - - - -
9c22b2df by Étienne Mollier at 2024-08-25T19:21:59+02:00
d/control: add myself to uploaders.
- - - - -
b00c0766 by Étienne Mollier at 2024-08-25T19:22:21+02:00
d/control: declare compliance to standards version 4.7.0.
- - - - -
f9dccf2c by Étienne Mollier at 2024-08-25T19:22:59+02:00
Ready for upload to unstable.
- - - - -
11 changed files:
- ChangeLog.md
- README.md
- configure.ac
- debian/changelog
- debian/control
- man/mptp.1
- src/auto.c
- src/dp.c
- src/likelihood.c
- src/mptp.h
- src/output.c
Changes:
=====================================
ChangeLog.md
=====================================
@@ -2,6 +2,13 @@
All notable changes to `mptp` will be documented in this file.
This project adheres to [Semantic Versioning](http://semver.org/).
+## [0.2.5] - 2023-09-11
+### Added
+ - Likelihood ratio test for the multi method
+ - Added implementation for the incomplete gamma function
+### Removed
+ - Dependency for GNU scientific library
+
## [0.2.4] - 2018-05-14
### Fixed
- If we do not manage to generate a random starting delimitation with the
=====================================
README.md
=====================================
@@ -78,9 +78,9 @@ where `DIR` is the directory where bash autocompletion is stored. You can use
and the documentation, use the following commands:
```bash
-wget https://github.com/Pas-Kapli/mptp/releases/download/v0.2.4/mptp-src-0.2.4.tar.gz
-tar zxvf mptp-src-0.2.4.tar.gz
-cd mptp-src-0.2.4
+wget https://github.com/Pas-Kapli/mptp/releases/download/v0.2.5/mptp-src-0.2.5.tar.gz
+tar zxvf mptp-src-0.2.5.tar.gz
+cd mptp-src-0.2.5
./configure
make
make install # as root, or run sudo make install
@@ -110,12 +110,12 @@ To use the pre-compiled binary, download the appropriate executable for your
system using the following commands if you are using a Linux system:
```bash
-wget https://github.com/Pas-Kapli/mptp/releases/download/v0.2.4/mptp-0.2.4-linux-x86_64.tar.gz
-tar zxvf mptp-0.2.4-linux-x86_64.tar.gz
+wget https://github.com/Pas-Kapli/mptp/releases/download/v0.2.5/mptp-0.2.5-linux-x86_64.tar.gz
+tar zxvf mptp-0.2.5-linux-x86_64.tar.gz
```
You will now have the binary distribution in a folder called
-`mptp-0.2.4-linux-x86_64` in which you will find three subfolders `bin`, `man`
+`mptp-0.2.5-linux-x86_64` in which you will find three subfolders `bin`, `man`
and `doc`. We recommend making a copy or a symbolic link to the mptp binary
`bin/mptp` in a folder included in your `$PATH`, and a copy or a symbolic link
to the mptp man page `man/mptp.1` in a folder included in your `$MANPATH`. The
=====================================
configure.ac
=====================================
@@ -2,7 +2,7 @@
# Process this file with autoconf to produce a configure script.
AC_PREREQ([2.63])
-AC_INIT([mptp], [0.2.4], [t.flouris at ucl.ac.uk])
+AC_INIT([mptp], [0.2.5], [t.flouris at ucl.ac.uk])
AM_INIT_AUTOMAKE([subdir-objects])
AC_LANG([C])
AC_CONFIG_SRCDIR([src/mptp.c])
@@ -42,8 +42,6 @@ AC_FUNC_REALLOC
AC_CHECK_FUNCS([memmove memcpy gettimeofday memchr memset pow regcomp strcasecmp strchr strcspn sysinfo])
AC_CHECK_LIB([m],[cos])
-AC_CHECK_LIB([gslcblas], [cblas_dgemm])
-AC_CHECK_LIB([gsl], [gsl_cdf_chisq_P])
# Bash completions
AC_ARG_WITH([bash-completions],
=====================================
debian/changelog
=====================================
@@ -1,3 +1,18 @@
+mptp (0.2.5-1) unstable; urgency=medium
+
+ [ Andreas Tille ]
+ * Fix syntax
+
+ [ Nilesh Patra ]
+ * [skip ci] Update email
+
+ [ Étienne Mollier ]
+ * New upstream version 0.2.5
+ * d/control: add myself to uploaders.
+ * d/control: declare compliance to standards version 4.7.0.
+
+ -- Étienne Mollier <emollier at debian.org> Sun, 25 Aug 2024 19:22:52 +0200
+
mptp (0.2.4-3) unstable; urgency=medium
* Add correct information about tests README
=====================================
debian/control
=====================================
@@ -1,7 +1,8 @@
Source: mptp
Maintainer: Debian Med Packaging Team <debian-med-packaging at lists.alioth.debian.org>
Uploaders: Andreas Tille <tille at debian.org>,
- Nilesh Patra <nilesh at debian.org>
+ Nilesh Patra <nilesh at debian.org>,
+ Étienne Mollier <emollier at debian.org>
Section: science
Priority: optional
Build-Depends: debhelper-compat (= 13),
@@ -10,7 +11,7 @@ Build-Depends: debhelper-compat (= 13),
libgsl-dev,
pkg-config,
bash-completion
-Standards-Version: 4.5.0
+Standards-Version: 4.7.0
Vcs-Browser: https://salsa.debian.org/med-team/mptp
Vcs-Git: https://salsa.debian.org/med-team/mptp.git
Homepage: https://github.com/Pas-Kapli/mptp
=====================================
man/mptp.1
=====================================
@@ -1,6 +1,6 @@
.\" -*- coding: utf-8 -*-
.\" ============================================================================
-.TH mptp 1 "May 14, 2018" "mptp 0.2.4" "USER COMMANDS"
+.TH mptp 1 "Sep 11, 2023" "mptp 0.2.5" "USER COMMANDS"
.\" ============================================================================
.SH NAME
mptp \(em single-locus species delimitation
@@ -359,5 +359,9 @@ Replaced hsearch() with custom hashtable. Fixed minor output error messages.
If we do not manage to generate a random starting delimitation with the wanted
number of species (randomly chosen), we use the currently generated
delimitation instead.
+.TP
+.BR v0.2.5\~ "released Sep 9th, 2023"
+Added likelihood ratio test for the multi method. Added implementation for the
+incomplete gamma function, and removed dependency for GNU scientific library.
.RE
.LP
=====================================
src/auto.c
=====================================
@@ -317,9 +317,15 @@ void detect_min_bl(rtree_t * rtree)
}
if (minfound && n != 1)
+ {
printf("Minimum branch length (--minbr) should be set to %.10f\n", branch_lengths[n-1]);
+ output_minbr(branch_lengths[n-1]);
+ }
else
+ {
printf("Minimum branch length (--minbr) should be set to 0\n");
+ output_minbr(0);
+ }
free(branch_lengths);
=====================================
src/dp.c
=====================================
@@ -22,6 +22,7 @@
#include "mptp.h"
static unsigned int species_iter = 0;
+static unsigned int coal_param_count = 0;
static void dp_recurse(rtree_t * node, long method)
{
@@ -166,6 +167,47 @@ static void backtrack(rtree_t * node,
}
}
+long multi_coalpopedgecount(rtree_t * node)
+{
+ long edges = 0;
+
+ if (node->left)
+ {
+ if (node->left->length > opt_minbr)
+ ++edges;
+ edges += multi_coalpopedgecount(node->left);
+ }
+ if (node->right)
+ {
+ if (node->right->length > opt_minbr)
+ ++edges;
+ edges += multi_coalpopedgecount(node->right);
+ }
+
+ return edges;
+
+}
+void multi_getcoalparamscount(rtree_t * node, int index)
+{
+ dp_vector_t * vec = node->vector;
+
+ if ((vec[index].vec_left != -1) && (vec[index].vec_right != -1))
+ {
+ node->event = EVENT_SPECIATION;
+
+ multi_getcoalparamscount(node->left, vec[index].vec_left);
+ multi_getcoalparamscount(node->right,vec[index].vec_right);
+ }
+ else
+ {
+ node->event = EVENT_COALESCENT;
+
+ long edges = multi_coalpopedgecount(node);
+ if (edges)
+ coal_param_count++;
+ }
+}
+
void dp_ptp(rtree_t * tree, long method)
{
int i;
@@ -234,19 +276,16 @@ void dp_ptp(rtree_t * tree, long method)
/* do a Likelihood Ratio Test (lrt) and return the computed p-value */
species_count = vec[best_index].species_count;
- // only do LRT for PTP, not for mPTP
- lrt_pass = (method == PTP_METHOD_MULTI) ? 1 : lrt(tree->coal_logl,
- vec[best_index].score_single, 1, &pvalue);
+ /* fills the coal_param_count variable with # coalescent pop parameters */
+ coal_param_count = 0;
+ multi_getcoalparamscount(tree,best_index);
-#ifndef HAVE_LIBGSL
- fprintf(stderr, "WARNING: delimit was not compiled with libgsl. "
- "Likelihood ratio test disabled.\n");
-#endif
+ /* likelihood ratio test */
+ unsigned int df = (method == PTP_METHOD_SINGLE) ? 1 : coal_param_count;
+ lrt_pass = lrt(tree->coal_logl,max,df,&pvalue);
-#ifdef HAVE_LIBGSL
- if (!opt_quiet && method == PTP_METHOD_SINGLE)
+ if (!opt_quiet)
fprintf(stdout,"LRT computed p-value: %.6f\n", pvalue);
-#endif
/* initialize file name */
FILE * out = open_file_ext("txt", opt_seed);
=====================================
src/likelihood.c
=====================================
@@ -21,6 +21,70 @@
#include "mptp.h"
+double IncompleteGamma(double x, double alpha, double ln_gamma_alpha)
+{
+ /* returns the incomplete gamma ratio I(x,alpha) where x is the upper
+ limit of the integration and alpha is the shape parameter.
+ returns (-1) if in error
+ ln_gamma_alpha = ln(Gamma(alpha)), is almost redundant.
+ (1) series expansion, if (alpha>x || x<=1)
+ (2) continued fraction, otherwise
+ RATNEST FORTRAN by
+ Bhattacharjee GP (1970) The incomplete gamma integral. Applied Statistics,
+ 19: 285-287 (AS32)
+ */
+ int i;
+ double p = alpha, g = ln_gamma_alpha;
+ double accurate = 1e-10, overflow = 1e60;
+ double factor, gin = 0, rn = 0, a = 0, b = 0, an = 0, dif = 0, term = 0, pn[6];
+
+ if (x == 0) return (0);
+ if (x < 0 || p <= 0) return (-1);
+
+ factor = exp(p*log(x) - x - g);
+ if (x > 1 && x >= p) goto l30;
+ /* (1) series expansion */
+ gin = 1; term = 1; rn = p;
+l20:
+ rn++;
+ term *= x / rn; gin += term;
+ if (term > accurate) goto l20;
+ gin *= factor / p;
+ goto l50;
+l30:
+ /* (2) continued fraction */
+ a = 1 - p; b = a + x + 1; term = 0;
+ pn[0] = 1; pn[1] = x; pn[2] = x + 1; pn[3] = x*b;
+ gin = pn[2] / pn[3];
+l32:
+ a++;
+ b += 2;
+ term++;
+ an = a*term;
+ for (i = 0; i < 2; i++)
+ pn[i + 4] = b*pn[i + 2] - an*pn[i];
+ if (pn[5] == 0) goto l35;
+ rn = pn[4] / pn[5];
+ dif = fabs(gin - rn);
+ if (dif > accurate) goto l34;
+ if (dif <= accurate*rn) goto l42;
+l34:
+ gin = rn;
+l35:
+ for (i = 0; i < 4; i++) pn[i] = pn[i + 2];
+ if (fabs(pn[4]) < overflow) goto l32;
+ for (i = 0; i < 4; i++) pn[i] /= overflow;
+ goto l32;
+l42:
+ gin = 1 - factor*gin;
+
+l50:
+ return (gin);
+}
+
+#define CDFGamma(x,alpha,beta) IncompleteGamma((beta)*(x),alpha,lgamma(alpha))
+#define CDFChi2(x,v) CDFGamma(x,(v)/2.0,0.5)
+
double loglikelihood(long edge_count, double edgelen_sum)
@@ -34,15 +98,11 @@ double loglikelihood(long edge_count, double edgelen_sum)
int lrt(double nullmodel_logl, double ptp_logl, unsigned int df, double * pvalue)
{
-#ifdef HAVE_LIBGSL
double diff = 2*(ptp_logl - nullmodel_logl);
-
- /* http://docs.scipy.org/doc/scipy/reference/generated/scipy.special.chdtr.html */
- *pvalue = 1 - gsl_cdf_chisq_P(diff,df);
+ *pvalue = 1 - CDFChi2(diff,df);
if ((*pvalue) > opt_pvalue)
return 0;
-#endif
return 1;
}
=====================================
src/mptp.h
=====================================
@@ -44,10 +44,6 @@
#include "config.h"
#endif
-#if (defined(HAVE_CONFIG_H) && defined(HAVE_LIBGSL))
-#include <gsl/gsl_cdf.h>
-#endif
-
/* constants */
#define PROG_NAME PACKAGE
@@ -426,6 +422,7 @@ void output_info(FILE * out,
int lrt_result,
rtree_t * root,
unsigned int species_count);
+void output_minbr(double minbr);
FILE * open_file_ext(const char * extension, long seed);
=====================================
src/output.c
=====================================
@@ -62,12 +62,25 @@ void output_info(FILE * out,
(method == PTP_METHOD_SINGLE) ?
"single" : "multi",
logl);
-#ifdef HAVE_LIBGSL
- if (method == PTP_METHOD_SINGLE)
- {
- fprintf(out, "LRT computed p-value: %.6f\n", pvalue);
- fprintf(out, "LRT: %s\n", lrt_result ? "passed" : "failed");
- }
-#endif
+ fprintf(out, "LRT computed p-value: %.6f\n", pvalue);
+ fprintf(out, "LRT: %s\n", lrt_result ? "passed" : "failed");
fprintf(out, "Number of delimited species: %d\n", species_count);
}
+
+void output_minbr(double minbr)
+{
+ FILE * fp_out;
+ char * filename = NULL;
+
+ if (asprintf(&filename, "%s.%s", opt_outfile, "txt") == -1)
+ fatal("Unable to allocate enough memory.");
+
+ fp_out = xopen(filename,"w");
+ free(filename);
+
+ fprintf(fp_out, "Command: %s\n", cmdline);
+ if (minbr == 0)
+ fprintf(fp_out, "Minimum branch length (--minbr) should be set to 0\n");
+ else
+ fprintf(fp_out, "Minimum branch length (--minbr) should be set to %f\n",minbr);
+}
View it on GitLab: https://salsa.debian.org/med-team/mptp/-/compare/af90077caaba893dd5ac6093d33e48b64fb0754b...f9dccf2cb7814b7e287ff004957649d3994c9c36
--
View it on GitLab: https://salsa.debian.org/med-team/mptp/-/compare/af90077caaba893dd5ac6093d33e48b64fb0754b...f9dccf2cb7814b7e287ff004957649d3994c9c36
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20240825/85f55f7b/attachment-0001.htm>
More information about the debian-med-commit
mailing list