Permalink
Browse files

C: linear solver can now return the residuals.

  • Loading branch information...
1 parent 7b384b9 commit 798e742ee80aa89bf4ce8ec43c198d8ee890d21c @gaborcsardi gaborcsardi committed Jan 5, 2012
Showing with 108 additions and 34 deletions.
  1. +1 −0 pysplicing/src/.gitignore
  2. +1 −0 splicing/src/.gitignore
  3. +16 −8 splicing/src/Rsplicing.c
  4. +30 −0 src/blas.c
  5. +47 −17 src/solve.c
  6. +6 −2 src/splicing.h
  7. +6 −6 src/splicing_lapack.h
  8. +1 −1 tools/getlapack.sh
@@ -20,3 +20,4 @@ util.c
miso_paired.c
lapack.c
complexity.c
+blas.c
View
@@ -105,3 +105,4 @@ vcf.c
lapack.c
splicing_lapack.h
complexity.c
+blas.c
View
@@ -1658,7 +1658,7 @@ SEXP R_splicing_solve_gene(SEXP pgff, SEXP pgene, SEXP preadLength,
int scale=LOGICAL(pscale)[0];
int i, noReads=GET_LENGTH(pcigar);
splicing_matrix_t match_matrix, assignment_matrix;
- splicing_vector_t expression;
+ splicing_vector_t expression, residuals;
double rnorm;
SEXP result, names;
@@ -1669,27 +1669,31 @@ SEXP R_splicing_solve_gene(SEXP pgff, SEXP pgene, SEXP preadLength,
splicing_matrix_init(&match_matrix, 0, 0);
splicing_matrix_init(&assignment_matrix, 0, 0);
splicing_vector_init(&expression, 0);
+ splicing_vector_init(&residuals, 0);
cigarstr = (const char**) R_alloc(noReads, sizeof(char*));
for (i=0; i<noReads; i++) {
cigarstr[i] = CHAR(STRING_ELT(pcigar, i));
}
splicing_solve_gene(&gff, gene, readLength, overhang, &position, cigarstr,
&match_matrix, &assignment_matrix, &expression,
- scale);
+ &residuals, scale);
- PROTECT(result=NEW_LIST(3));
+ PROTECT(result=NEW_LIST(4));
SET_VECTOR_ELT(result, 0, R_splicing_matrix_to_SEXP(&match_matrix));
splicing_matrix_destroy(&match_matrix);
SET_VECTOR_ELT(result, 1, R_splicing_matrix_to_SEXP(&assignment_matrix));
splicing_matrix_destroy(&assignment_matrix);
SET_VECTOR_ELT(result, 2, R_splicing_vector_to_SEXP(&expression));
splicing_vector_destroy(&expression);
+ SET_VECTOR_ELT(result, 3, R_splicing_vector_to_SEXP(&residuals));
+ splicing_vector_destroy(&residuals);
- PROTECT(names=NEW_CHARACTER(3));
+ PROTECT(names=NEW_CHARACTER(4));
SET_STRING_ELT(names, 0, mkChar("match"));
SET_STRING_ELT(names, 1, mkChar("assignment"));
SET_STRING_ELT(names, 2, mkChar("expression"));
+ SET_STRING_ELT(names, 3, mkChar("residuals"));
SET_NAMES(result, names);
R_splicing_end();
@@ -1719,7 +1723,7 @@ SEXP R_splicing_solve_gene_paired(SEXP pgff, SEXP pgene, SEXP preadLength,
double normalVar=REAL(pnormalVar)[0];
double numDevs=REAL(pnumDevs)[0];
splicing_matrix_t match_matrix, assignment_matrix;
- splicing_vector_t expression;
+ splicing_vector_t expression, residuals;
double rnorm;
SEXP result, names;
@@ -1733,6 +1737,7 @@ SEXP R_splicing_solve_gene_paired(SEXP pgff, SEXP pgene, SEXP preadLength,
splicing_matrix_init(&match_matrix, 0, 0);
splicing_matrix_init(&assignment_matrix, 0, 0);
splicing_vector_init(&expression, 0);
+ splicing_vector_init(&residuals, 0);
cigarstr = (const char**) R_alloc(noReads, sizeof(char*));
for (i=0; i<noReads; i++) {
cigarstr[i] = CHAR(STRING_ELT(pcigar, i));
@@ -1743,20 +1748,23 @@ SEXP R_splicing_solve_gene_paired(SEXP pgff, SEXP pgene, SEXP preadLength,
isNull(pfragmentprob) ? 0 : &fragmentprob,
fragmentstart, normalMean, normalVar, numDevs,
&match_matrix, &assignment_matrix, &expression,
- scale);
+ &residuals, scale);
- PROTECT(result=NEW_LIST(3));
+ PROTECT(result=NEW_LIST(4));
SET_VECTOR_ELT(result, 0, R_splicing_matrix_to_SEXP(&match_matrix));
splicing_matrix_destroy(&match_matrix);
SET_VECTOR_ELT(result, 1, R_splicing_matrix_to_SEXP(&assignment_matrix));
splicing_matrix_destroy(&assignment_matrix);
SET_VECTOR_ELT(result, 2, R_splicing_vector_to_SEXP(&expression));
splicing_vector_destroy(&expression);
+ SET_VECTOR_ELT(result, 3, R_splicing_vector_to_SEXP(&residuals));
+ splicing_vector_destroy(&residuals);
- PROTECT(names=NEW_CHARACTER(3));
+ PROTECT(names=NEW_CHARACTER(4));
SET_STRING_ELT(names, 0, mkChar("match"));
SET_STRING_ELT(names, 1, mkChar("assignment"));
SET_STRING_ELT(names, 2, mkChar("expression"));
+ SET_STRING_ELT(names, 3, mkChar("residuals"));
SET_NAMES(result, names);
R_splicing_end();
View
@@ -0,0 +1,30 @@
+
+#include "splicing.h"
+#include "splicing_error.h"
+
+int splicingdgemv_(char *trans, int *m, int *n, double *alpha,
+ double *a, int *lda, double *x, int *incx,
+ double *beta, double *y, int *incy);
+
+int splicing_dgemv(int transpose, double alpha,
+ const splicing_matrix_t* a, const splicing_vector_t* x,
+ double beta, splicing_vector_t* y) {
+ char trans = transpose ? 'T' : 'N';
+ int m, n;
+ int inc = 1;
+
+ m = splicing_matrix_nrow(a);
+ n = splicing_matrix_ncol(a);
+
+ if (splicing_vector_size(x) != (transpose ? m : n)) {
+ SPLICING_ERROR("Invalid matrix-vector sizes", SPLICING_EINVAL);
+ }
+ if (splicing_vector_size(y) != (transpose ? n : m)) {
+ SPLICING_ERROR("Invalid matrix-vector sizes", SPLICING_EINVAL);
+ }
+
+ splicingdgemv_(&trans, &m, &n, &alpha, VECTOR(a->data), &m,
+ VECTOR(*x), &inc, &beta, VECTOR(*y), &inc);
+
+ return 0;
+}
View
@@ -284,10 +284,10 @@ int splicing_solve_gene(const splicing_gff_t *gff, size_t gene,
splicing_matrix_t *match_matrix,
splicing_matrix_t *assignment_matrix,
splicing_vector_t *expression,
- int scale) {
+ splicing_vector_t *residuals, int scale) {
splicing_matrix_t A;
- splicing_vector_t match;
+ splicing_vector_t match, omatch;
splicing_vector_long_t index;
size_t no_classes;
size_t no_reads=splicing_vector_int_size(position);
@@ -339,6 +339,22 @@ int splicing_solve_gene(const splicing_gff_t *gff, size_t gene,
SPLICING_CHECK(splicing_matrix_transpose(&A));
SPLICING_CHECK(splicing_vector_long_init(&index, 0));
SPLICING_FINALLY(splicing_vector_long_destroy, &index);
+ SPLICING_CHECK(splicing_vector_init(&omatch, no_classes));
+ SPLICING_FINALLY(splicing_vector_destroy, &omatch);
+ splicing_vector_update(&omatch, &match);
+
+ SPLICING_CHECK(splicing_nnls(&A, &match, expression, &rnorm, &index,
+ &nsetp));
+
+ if (residuals) {
+ int i, n=splicing_vector_size(&match);
+ splicing_dgemv(/*transpose=*/ 1, /*alpha=*/ 1.0, myass_matrix,
+ expression, /*beta=*/ 0.0, &match);
+ SPLICING_CHECK(splicing_vector_resize(residuals, no_classes));
+ for (i=0; i<n; i++) {
+ VECTOR(*residuals)[i] = VECTOR(omatch)[i] - VECTOR(match)[i];
+ }
+ }
if (!assignment_matrix) {
splicing_matrix_destroy(myass_matrix);
@@ -348,17 +364,16 @@ int splicing_solve_gene(const splicing_gff_t *gff, size_t gene,
splicing_matrix_destroy(mymatch_matrix);
SPLICING_FINALLY_CLEAN(1);
}
-
- SPLICING_CHECK(splicing_nnls(&A, &match, expression, &rnorm, &index,
- &nsetp));
-
+
+ splicing_vector_destroy(&omatch);
splicing_vector_long_destroy(&index);
splicing_matrix_destroy(&A);
splicing_vector_destroy(&match);
- SPLICING_FINALLY_CLEAN(3);
+ SPLICING_FINALLY_CLEAN(4);
if (scale) {
- splicing_vector_scale(expression, 1.0/splicing_vector_sum(expression));
+ double fac=1.0/splicing_vector_sum(expression);
+ splicing_vector_scale(expression, fac);
}
return 0;
@@ -374,10 +389,10 @@ int splicing_solve_gene_paired(const splicing_gff_t *gff, size_t gene,
splicing_matrix_t *match_matrix,
splicing_matrix_t *assignment_matrix,
splicing_vector_t *expression,
- int scale) {
+ splicing_vector_t *residuals, int scale) {
splicing_matrix_t A;
- splicing_vector_t match;
+ splicing_vector_t match, omatch;
splicing_vector_long_t index;
size_t no_classes;
size_t no_reads=splicing_vector_int_size(position);
@@ -442,6 +457,22 @@ int splicing_solve_gene_paired(const splicing_gff_t *gff, size_t gene,
SPLICING_CHECK(splicing_matrix_transpose(&A));
SPLICING_CHECK(splicing_vector_long_init(&index, 0));
SPLICING_FINALLY(splicing_vector_long_destroy, &index);
+ SPLICING_CHECK(splicing_vector_init(&omatch, no_classes));
+ SPLICING_FINALLY(splicing_vector_destroy, &omatch);
+ splicing_vector_update(&omatch, &match);
+
+ SPLICING_CHECK(splicing_nnls(&A, &match, expression, &rnorm, &index,
+ &nsetp));
+
+ if (residuals) {
+ int i, n=splicing_vector_size(&match);
+ splicing_dgemv(/*transpose=*/ 1, /*alpha=*/ 1.0, myass_matrix,
+ expression, /*beta=*/ 0.0, &match);
+ SPLICING_CHECK(splicing_vector_resize(residuals, no_classes));
+ for (i=0; i<n; i++) {
+ VECTOR(*residuals)[i] = VECTOR(omatch)[i] - VECTOR(match)[i];
+ }
+ }
if (!assignment_matrix) {
splicing_matrix_destroy(myass_matrix);
@@ -450,18 +481,17 @@ int splicing_solve_gene_paired(const splicing_gff_t *gff, size_t gene,
if (!match_matrix) {
splicing_matrix_destroy(mymatch_matrix);
SPLICING_FINALLY_CLEAN(1);
- }
-
- SPLICING_CHECK(splicing_nnls(&A, &match, expression, &rnorm, &index,
- &nsetp));
-
+ }
+
+ splicing_vector_destroy(&omatch);
splicing_vector_long_destroy(&index);
splicing_matrix_destroy(&A);
splicing_vector_destroy(&match);
- SPLICING_FINALLY_CLEAN(3);
+ SPLICING_FINALLY_CLEAN(4);
if (scale) {
- splicing_vector_scale(expression, 1.0/splicing_vector_sum(expression));
+ double fac=1.0/splicing_vector_sum(expression);
+ splicing_vector_scale(expression, fac);
}
return 0;
View
@@ -352,7 +352,7 @@ int splicing_solve_gene(const splicing_gff_t *gff, size_t gene,
splicing_matrix_t *match_matrix,
splicing_matrix_t *assignment_matrix,
splicing_vector_t *expression,
- int scale);
+ splicing_vector_t *residuals, int scale);
int splicing_solve_gene_paired(const splicing_gff_t *gff, size_t gene,
int readLength, int overHang,
@@ -364,7 +364,7 @@ int splicing_solve_gene_paired(const splicing_gff_t *gff, size_t gene,
splicing_matrix_t *match_matrix,
splicing_matrix_t *assignment_matrix,
splicing_vector_t *expression,
- int scale);
+ splicing_vector_t *residuals, int scale);
typedef struct splicing_gff_converter_t {
size_t noiso;
@@ -462,6 +462,10 @@ int splicing_metropolis_hastings_ratio_paired(
splicing_vector_t *pcJS,
splicing_vector_t *ppJS);
+int splicing_dgemv(int transpose, double alpha,
+ const splicing_matrix_t* a, const splicing_vector_t* x,
+ double beta, splicing_vector_t* y);
+
int splicing_dgesdd(const splicing_matrix_t *matrix,
splicing_vector_t *values);
View
@@ -2,15 +2,15 @@
#ifdef EXTERNAL_LAPACK
#define splicingdcopy_ dcopy_
#define splicinglsame_ lsame_
-#define splicingdrot _ drot _
-#define splicingdgemm _ dgemm _
+#define splicingdrot_ drot_
+#define splicingdgemm_ dgemm_
#define splicingdnrm2_ dnrm2_
-#define splicingdscal _ dscal _
+#define splicingdscal_ dscal_
#define splicingdswap_ dswap_
#define splicingddot_ ddot_
-#define splicingdgemv _ dgemv _
-#define splicingdger _ dger _
-#define splicingdtrmm _ dtrmm _
+#define splicingdgemv_ dgemv_
+#define splicingdger_ dger_
+#define splicingdtrmm_ dtrmm_
#define splicingdtrmv_ dtrmv_
#define splicingdgesdd_ dgesdd_
#define splicingdbdsdc_ dbdsdc_
View
@@ -1,6 +1,6 @@
#! /bin/sh
#
-# ./getlapack.sh dgesdd
+# ./getlapack.sh dgesdd dgemv
#
make

0 comments on commit 798e742

Please sign in to comment.