From d0051dc9939d3477bd92b42c86bcd3eda743b955 Mon Sep 17 00:00:00 2001 From: westcott Date: Fri, 23 Dec 2011 13:45:29 +0000 Subject: [PATCH] added qvalues to metastats. fixed bug with chimera.uchime location. fixed windows bug with tellg() in trim.seqs. --- Mothur.xcodeproj/project.pbxproj | 83 ++++++ bsplvb.f | 102 +++++++ bsplvd.f | 118 ++++++++ bvalue.f | 135 ++++++++++ chimerauchimecommand.cpp | 5 + chimerauchimecommand.h | 8 + daxpy.f | 69 +++++ ddot.f | 90 +++++++ dpbfa.f | 109 ++++++++ dpbsl.f | 83 ++++++ interv.f | 102 +++++++ intrv.f | 117 ++++++++ makefile | 11 +- mothur.h | 3 +- mothurmetastats.cpp | 204 +++++--------- mothurmetastats.h | 8 +- phylotree.cpp | 2 + sgram.f | 142 ++++++++++ sinerp.f | 98 +++++++ spline.cpp | 450 +++++++++++++++++++++++++++++++ spline.h | 38 +++ sslvrg.f | 136 ++++++++++ stxwx.f | 65 +++++ trimflowscommand.cpp | 2 +- trimseqscommand.cpp | 28 +- 25 files changed, 2026 insertions(+), 182 deletions(-) create mode 100644 bsplvb.f create mode 100644 bsplvd.f create mode 100644 bvalue.f create mode 100644 daxpy.f create mode 100644 ddot.f create mode 100644 dpbfa.f create mode 100644 dpbsl.f create mode 100644 interv.f create mode 100644 intrv.f create mode 100644 sgram.f create mode 100644 sinerp.f create mode 100644 spline.cpp create mode 100644 spline.h create mode 100644 sslvrg.f create mode 100644 stxwx.f diff --git a/Mothur.xcodeproj/project.pbxproj b/Mothur.xcodeproj/project.pbxproj index d1cb6d8..3a37569 100644 --- a/Mothur.xcodeproj/project.pbxproj +++ b/Mothur.xcodeproj/project.pbxproj @@ -18,6 +18,7 @@ A727864412E9E28C00F86ABA /* removerarecommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A727864312E9E28C00F86ABA /* removerarecommand.cpp */; }; A73DDBBA13C4A0D1006AAE38 /* clearmemorycommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A73DDBB913C4A0D1006AAE38 /* clearmemorycommand.cpp */; }; A73DDC3813C4BF64006AAE38 /* mothurmetastats.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A73DDC3713C4BF64006AAE38 /* mothurmetastats.cpp */; }; + A74A9A9F148E881E00AB5E3E /* spline.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A74A9A9E148E881E00AB5E3E /* spline.cpp */; }; A74D3687137DAB8300332B0C /* addtargets2.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A74D3655137DAB8300332B0C /* addtargets2.cpp */; }; A74D3688137DAB8400332B0C /* alignchime.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A74D3656137DAB8300332B0C /* alignchime.cpp */; }; A74D3689137DAB8400332B0C /* alignchimel.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A74D3657137DAB8300332B0C /* alignchimel.cpp */; }; @@ -325,6 +326,19 @@ A7E9B98F12D37EC400DA6239 /* whittaker.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A7E9B87F12D37EC400DA6239 /* whittaker.cpp */; }; A7F9F5CF141A5E500032F693 /* sequenceparser.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A7F9F5CE141A5E500032F693 /* sequenceparser.cpp */; }; A7FA10021302E097003860FE /* mantelcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A7FA10011302E096003860FE /* mantelcommand.cpp */; }; + A7FA2AC714A0E881007C09A6 /* bsplvb.f in Sources */ = {isa = PBXBuildFile; fileRef = A7FA2ABC14A0E881007C09A6 /* bsplvb.f */; }; + A7FA2AC814A0E881007C09A6 /* bsplvd.f in Sources */ = {isa = PBXBuildFile; fileRef = A7FA2ABD14A0E881007C09A6 /* bsplvd.f */; }; + A7FA2AC914A0E881007C09A6 /* bvalue.f in Sources */ = {isa = PBXBuildFile; fileRef = A7FA2ABE14A0E881007C09A6 /* bvalue.f */; }; + A7FA2ACA14A0E881007C09A6 /* daxpy.f in Sources */ = {isa = PBXBuildFile; fileRef = A7FA2ABF14A0E881007C09A6 /* daxpy.f */; }; + A7FA2ACB14A0E881007C09A6 /* ddot.f in Sources */ = {isa = PBXBuildFile; fileRef = A7FA2AC014A0E881007C09A6 /* ddot.f */; }; + A7FA2ACC14A0E881007C09A6 /* dpbfa.f in Sources */ = {isa = PBXBuildFile; fileRef = A7FA2AC114A0E881007C09A6 /* dpbfa.f */; }; + A7FA2ACD14A0E881007C09A6 /* dpbsl.f in Sources */ = {isa = PBXBuildFile; fileRef = A7FA2AC214A0E881007C09A6 /* dpbsl.f */; }; + A7FA2ACE14A0E881007C09A6 /* interv.f in Sources */ = {isa = PBXBuildFile; fileRef = A7FA2AC314A0E881007C09A6 /* interv.f */; }; + A7FA2ACF14A0E881007C09A6 /* sgram.f in Sources */ = {isa = PBXBuildFile; fileRef = A7FA2AC414A0E881007C09A6 /* sgram.f */; }; + A7FA2AD014A0E881007C09A6 /* sinerp.f in Sources */ = {isa = PBXBuildFile; fileRef = A7FA2AC514A0E881007C09A6 /* sinerp.f */; }; + A7FA2AD114A0E881007C09A6 /* stxwx.f in Sources */ = {isa = PBXBuildFile; fileRef = A7FA2AC614A0E881007C09A6 /* stxwx.f */; }; + A7FA2B1614A0EBEA007C09A6 /* sslvrg.f in Sources */ = {isa = PBXBuildFile; fileRef = A7FA2B1514A0EBEA007C09A6 /* sslvrg.f */; }; + A7FA2B5B14A0F0C2007C09A6 /* intrv.f in Sources */ = {isa = PBXBuildFile; fileRef = A7FA2B5A14A0F0C2007C09A6 /* intrv.f */; }; A7FC480E12D788F20055BC5C /* linearalgebra.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A7FC480D12D788F20055BC5C /* linearalgebra.cpp */; }; A7FC486712D795D60055BC5C /* pcacommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A7FC486612D795D60055BC5C /* pcacommand.cpp */; }; A7FE7C401330EA1000F7B327 /* getcurrentcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A7FE7C3F1330EA1000F7B327 /* getcurrentcommand.cpp */; }; @@ -333,6 +347,19 @@ A7FFB558142CA02C004884F2 /* summarytaxcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A7FFB557142CA02C004884F2 /* summarytaxcommand.cpp */; }; /* End PBXBuildFile section */ +/* Begin PBXBuildRule section */ + A7D162CB149F96CA000523E8 /* PBXBuildRule */ = { + isa = PBXBuildRule; + compilerSpec = com.apple.compilers.proxy.script; + fileType = sourcecode.fortran; + isEditable = 1; + outputFiles = ( + "$(TARGET_BUILD_DIR)/$(INPUT_FILE_BASE).o", + ); + script = "gfortran -g -m64 -c ${PROJECT_DIR}/${INPUT_FILE_NAME} -o ${TARGET_BUILD_DIR}/${INPUT_FILE_BASE}.o"; + }; +/* End PBXBuildRule section */ + /* Begin PBXCopyFilesBuildPhase section */ 8DD76FAF0486AB0100D96B5E /* CopyFiles */ = { isa = PBXCopyFilesBuildPhase; @@ -368,6 +395,8 @@ A73DDBB913C4A0D1006AAE38 /* clearmemorycommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = clearmemorycommand.cpp; sourceTree = ""; }; A73DDC3613C4BF64006AAE38 /* mothurmetastats.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = mothurmetastats.h; sourceTree = ""; }; A73DDC3713C4BF64006AAE38 /* mothurmetastats.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = mothurmetastats.cpp; sourceTree = ""; }; + A74A9A9D148E881E00AB5E3E /* spline.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = spline.h; sourceTree = ""; }; + A74A9A9E148E881E00AB5E3E /* spline.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = spline.cpp; sourceTree = ""; }; A74D3655137DAB8300332B0C /* addtargets2.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = addtargets2.cpp; sourceTree = ""; }; A74D3656137DAB8300332B0C /* alignchime.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = alignchime.cpp; sourceTree = ""; }; A74D3657137DAB8300332B0C /* alignchimel.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = alignchimel.cpp; sourceTree = ""; }; @@ -1010,6 +1039,19 @@ A7F9F5CE141A5E500032F693 /* sequenceparser.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = sequenceparser.cpp; sourceTree = ""; }; A7FA10001302E096003860FE /* mantelcommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = mantelcommand.h; sourceTree = ""; }; A7FA10011302E096003860FE /* mantelcommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = mantelcommand.cpp; sourceTree = ""; }; + A7FA2ABC14A0E881007C09A6 /* bsplvb.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = bsplvb.f; sourceTree = ""; }; + A7FA2ABD14A0E881007C09A6 /* bsplvd.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = bsplvd.f; sourceTree = ""; }; + A7FA2ABE14A0E881007C09A6 /* bvalue.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = bvalue.f; sourceTree = ""; }; + A7FA2ABF14A0E881007C09A6 /* daxpy.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = daxpy.f; sourceTree = ""; }; + A7FA2AC014A0E881007C09A6 /* ddot.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = ddot.f; sourceTree = ""; }; + A7FA2AC114A0E881007C09A6 /* dpbfa.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = dpbfa.f; sourceTree = ""; }; + A7FA2AC214A0E881007C09A6 /* dpbsl.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = dpbsl.f; sourceTree = ""; }; + A7FA2AC314A0E881007C09A6 /* interv.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = interv.f; sourceTree = ""; }; + A7FA2AC414A0E881007C09A6 /* sgram.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = sgram.f; sourceTree = ""; }; + A7FA2AC514A0E881007C09A6 /* sinerp.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = sinerp.f; sourceTree = ""; }; + A7FA2AC614A0E881007C09A6 /* stxwx.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = stxwx.f; sourceTree = ""; }; + A7FA2B1514A0EBEA007C09A6 /* sslvrg.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = sslvrg.f; sourceTree = ""; }; + A7FA2B5A14A0F0C2007C09A6 /* intrv.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = intrv.f; sourceTree = ""; }; A7FC480C12D788F20055BC5C /* linearalgebra.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = linearalgebra.h; sourceTree = ""; }; A7FC480D12D788F20055BC5C /* linearalgebra.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = linearalgebra.cpp; sourceTree = ""; }; A7FC486512D795D60055BC5C /* pcacommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = pcacommand.h; sourceTree = ""; }; @@ -1212,6 +1254,26 @@ name = uchime; sourceTree = ""; }; + A7D161E7149F7F50000523E8 /* fortran */ = { + isa = PBXGroup; + children = ( + A7FA2B5A14A0F0C2007C09A6 /* intrv.f */, + A7FA2B1514A0EBEA007C09A6 /* sslvrg.f */, + A7FA2ABC14A0E881007C09A6 /* bsplvb.f */, + A7FA2ABD14A0E881007C09A6 /* bsplvd.f */, + A7FA2ABE14A0E881007C09A6 /* bvalue.f */, + A7FA2ABF14A0E881007C09A6 /* daxpy.f */, + A7FA2AC014A0E881007C09A6 /* ddot.f */, + A7FA2AC114A0E881007C09A6 /* dpbfa.f */, + A7FA2AC214A0E881007C09A6 /* dpbsl.f */, + A7FA2AC314A0E881007C09A6 /* interv.f */, + A7FA2AC414A0E881007C09A6 /* sgram.f */, + A7FA2AC514A0E881007C09A6 /* sinerp.f */, + A7FA2AC614A0E881007C09A6 /* stxwx.f */, + ); + name = fortran; + sourceTree = ""; + }; A7E9BA3812D3956100DA6239 /* commands */ = { isa = PBXGroup; children = ( @@ -1795,6 +1857,7 @@ A7E9BA5612D39BD800DA6239 /* metastats */ = { isa = PBXGroup; children = ( + A7D161E7149F7F50000523E8 /* fortran */, A7E9B6E512D37EC400DA6239 /* fisher2.c */, A7E9B6E612D37EC400DA6239 /* fisher2.h */, A7E9B75512D37EC400DA6239 /* metastats.h */, @@ -1803,6 +1866,8 @@ A79234D613C74BF6002B08E2 /* mothurfisher.cpp */, A73DDC3613C4BF64006AAE38 /* mothurmetastats.h */, A73DDC3713C4BF64006AAE38 /* mothurmetastats.cpp */, + A74A9A9D148E881E00AB5E3E /* spline.h */, + A74A9A9E148E881E00AB5E3E /* spline.cpp */, ); name = metastats; sourceTree = ""; @@ -1827,6 +1892,7 @@ 8DD76FAF0486AB0100D96B5E /* CopyFiles */, ); buildRules = ( + A7D162CB149F96CA000523E8 /* PBXBuildRule */, ); dependencies = ( ); @@ -2191,6 +2257,20 @@ A774104814696F320098E6AC /* myseqdist.cpp in Sources */, A77410F614697C300098E6AC /* seqnoise.cpp in Sources */, A754149714840CF7005850D1 /* summaryqualcommand.cpp in Sources */, + A74A9A9F148E881E00AB5E3E /* spline.cpp in Sources */, + A7FA2AC714A0E881007C09A6 /* bsplvb.f in Sources */, + A7FA2AC814A0E881007C09A6 /* bsplvd.f in Sources */, + A7FA2AC914A0E881007C09A6 /* bvalue.f in Sources */, + A7FA2ACA14A0E881007C09A6 /* daxpy.f in Sources */, + A7FA2ACB14A0E881007C09A6 /* ddot.f in Sources */, + A7FA2ACC14A0E881007C09A6 /* dpbfa.f in Sources */, + A7FA2ACD14A0E881007C09A6 /* dpbsl.f in Sources */, + A7FA2ACE14A0E881007C09A6 /* interv.f in Sources */, + A7FA2ACF14A0E881007C09A6 /* sgram.f in Sources */, + A7FA2AD014A0E881007C09A6 /* sinerp.f in Sources */, + A7FA2AD114A0E881007C09A6 /* stxwx.f in Sources */, + A7FA2B1614A0EBEA007C09A6 /* sslvrg.f in Sources */, + A7FA2B5B14A0F0C2007C09A6 /* intrv.f in Sources */, ); runOnlyForDeploymentPostprocessing = 0; }; @@ -2240,9 +2320,11 @@ "VERSION=\"\\\"1.22.0\\\"\"", "RELEASE_DATE=\"\\\"10/12/2011\\\"\"", ); + "GCC_VERSION[arch=*]" = ""; GCC_WARN_ABOUT_MISSING_NEWLINE = YES; GCC_WARN_ABOUT_RETURN_TYPE = YES; GCC_WARN_UNUSED_VARIABLE = YES; + HEADER_SEARCH_PATHS = ""; INSTALL_PATH = ""; MACH_O_TYPE = mh_execute; ONLY_ACTIVE_ARCH = YES; @@ -2257,6 +2339,7 @@ ); PREBINDING = NO; SDKROOT = macosx10.6; + USER_HEADER_SEARCH_PATHS = ""; }; name = Debug; }; diff --git a/bsplvb.f b/bsplvb.f new file mode 100644 index 0000000..64ca95b --- /dev/null +++ b/bsplvb.f @@ -0,0 +1,102 @@ + subroutine bsplvb ( t, lent,jhigh, index, x, left, biatx ) +c implicit none +c ------------- + +calculates the value of all possibly nonzero b-splines at x of order +c +c jout = dmax( jhigh , (j+1)*(index-1) ) +c +c with knot sequence t . +c +c****** i n p u t ****** +c t.....knot sequence, of length left + jout , assumed to be nonde- +c creasing. +c a s s u m p t i o n : t(left) < t(left + 1) +c d i v i s i o n b y z e r o will result if t(left) = t(left+1) +c +c jhigh, +c index.....integers which determine the order jout = max(jhigh, +c (j+1)*(index-1)) of the b-splines whose values at x are to +c be returned. index is used to avoid recalculations when seve- +c ral columns of the triangular array of b-spline values are nee- +c ded (e.g., in bvalue or in bsplvd ). precisely, +c if index = 1 , +c the calculation starts from scratch and the entire triangular +c array of b-spline values of orders 1,2,...,jhigh is generated +c order by order , i.e., column by column . +c if index = 2 , +c only the b-spline values of order j+1, j+2, ..., jout are ge- +c nerated, the assumption being that biatx , j , deltal , deltar +c are, on entry, as they were on exit at the previous call. +c in particular, if jhigh = 0, then jout = j+1, i.e., just +c the next column of b-spline values is generated. +c +c w a r n i n g . . . the restriction jout <= jmax (= 20) is +c imposed arbitrarily by the dimension statement for deltal and +c deltar below, but is n o w h e r e c h e c k e d for . +c +c x.....the point at which the b-splines are to be evaluated. +c left.....an integer chosen (usually) so that +c t(left) <= x <= t(left+1) . +c +c****** o u t p u t ****** +c biatx.....array of length jout , with biatx(i) containing the val- +c ue at x of the polynomial of order jout which agrees with +c the b-spline b(left-jout+i,jout,t) on the interval (t(left), +c t(left+1)) . +c +c****** m e t h o d ****** +c the recurrence relation +c +c x - t(i) t(i+j+1) - x +c b(i,j+1)(x) = ----------- b(i,j)(x) + --------------- b(i+1,j)(x) +c t(i+j)-t(i) t(i+j+1)-t(i+1) +c +c is used (repeatedly) to generate the +c (j+1)-vector b(left-j,j+1)(x),...,b(left,j+1)(x) +c from the j-vector b(left-j+1,j)(x),...,b(left,j)(x), +c storing the new values in biatx over the old. the facts that +c b(i,1) = 1 if t(i) <= x < t(i+1) +c and that +c b(i,j)(x) = 0 unless t(i) <= x < t(i+j) +c are used. the particular organization of the calculations follows +c algorithm (8) in chapter x of the text. +c + +C Arguments + integer lent, jhigh, index, left + double precision t(lent),x, biatx(jhigh) +c dimension t(left+jout), biatx(jout) +c ----------------------------------- +c current fortran standard makes it impossible to specify the length of +c t and of biatx precisely without the introduction of otherwise +c superfluous additional arguments. + +C Local Variables + integer jmax + parameter(jmax = 20) + integer i,j,jp1 + double precision deltal(jmax), deltar(jmax),saved,term + + save j,deltal,deltar + data j/1/ +c + go to (10,20), index + 10 j = 1 + biatx(1) = 1d0 + if (j .ge. jhigh) go to 99 +c + 20 jp1 = j + 1 + deltar(j) = t(left+j) - x + deltal(j) = x - t(left+1-j) + saved = 0d0 + do 26 i=1,j + term = biatx(i)/(deltar(i) + deltal(jp1-i)) + biatx(i) = saved + deltar(i)*term + 26 saved = deltal(jp1-i)*term + biatx(jp1) = saved + j = jp1 + if (j .lt. jhigh) go to 20 +c + 99 return + end diff --git a/bsplvd.f b/bsplvd.f new file mode 100644 index 0000000..47a5586 --- /dev/null +++ b/bsplvd.f @@ -0,0 +1,118 @@ + subroutine bsplvd ( t, lent, k, x, left, a, dbiatx, nderiv ) +c -------- ------ +c implicit none + +C calculates value and deriv.s of all b-splines which do not vanish at x +C calls bsplvb +c +c****** i n p u t ****** +c t the knot array, of length left+k (at least) +c k the order of the b-splines to be evaluated +c x the point at which these values are sought +c left an integer indicating the left endpoint of the interval of +c interest. the k b-splines whose support contains the interval +c (t(left), t(left+1)) +c are to be considered. +c a s s u m p t i o n - - - it is assumed that +c t(left) < t(left+1) +c division by zero will result otherwise (in b s p l v b ). +c also, the output is as advertised only if +c t(left) <= x <= t(left+1) . +c nderiv an integer indicating that values of b-splines and their +c derivatives up to but not including the nderiv-th are asked +c for. ( nderiv is replaced internally by the integer in (1,k) +c closest to it.) +c +c****** w o r k a r e a ****** +c a an array of order (k,k), to contain b-coeff.s of the derivat- +c ives of a certain order of the k b-splines of interest. +c +c****** o u t p u t ****** +c dbiatx an array of order (k,nderiv). its entry (i,m) contains +c value of (m-1)st derivative of (left-k+i)-th b-spline of +c order k for knot sequence t , i=m,...,k; m=1,...,nderiv. +c +c****** m e t h o d ****** +c values at x of all the relevant b-splines of order k,k-1,..., +c k+1-nderiv are generated via bsplvb and stored temporarily +c in dbiatx . then, the b-coeffs of the required derivatives of the +c b-splines of interest are generated by differencing, each from the +c preceding one of lower order, and combined with the values of b- +c splines of corresponding order in dbiatx to produce the desired +c values. + +C Args + integer lent,k,left,nderiv + double precision t(lent),x, dbiatx(k,nderiv), a(k,k) +C Locals + double precision factor,fkp1mm,sum + integer i,ideriv,il,j,jlow,jp1mid, kp1,kp1mm,ldummy,m,mhigh + + mhigh = max0(min0(nderiv,k),1) +c mhigh is usually equal to nderiv. + kp1 = k+1 + call bsplvb(t,lent,kp1-mhigh,1,x,left,dbiatx) + if (mhigh .eq. 1) go to 99 +c the first column of dbiatx always contains the b-spline values +c for the current order. these are stored in column k+1-current +c order before bsplvb is called to put values for the next +c higher order on top of it. + ideriv = mhigh + do 15 m=2,mhigh + jp1mid = 1 + do 11 j=ideriv,k + dbiatx(j,ideriv) = dbiatx(jp1mid,1) + 11 jp1mid = jp1mid + 1 + ideriv = ideriv - 1 + call bsplvb(t,lent,kp1-ideriv,2,x,left,dbiatx) + 15 continue +c +c at this point, b(left-k+i, k+1-j)(x) is in dbiatx(i,j) for +c i=j,...,k and j=1,...,mhigh ('=' nderiv). in particular, the +c first column of dbiatx is already in final form. to obtain cor- +c responding derivatives of b-splines in subsequent columns, gene- +c rate their b-repr. by differencing, then evaluate at x. +c + jlow = 1 + do 20 i=1,k + do 19 j=jlow,k + 19 a(j,i) = 0d0 + jlow = i + 20 a(i,i) = 1d0 +c at this point, a(.,j) contains the b-coeffs for the j-th of the +c k b-splines of interest here. +c + do 40 m=2,mhigh + kp1mm = kp1 - m + fkp1mm = dble(kp1mm) + il = left + i = k +c +c for j=1,...,k, construct b-coeffs of (m-1)st derivative of +c b-splines from those for preceding derivative by differencing +c and store again in a(.,j) . the fact that a(i,j) = 0 for +c i < j is used.sed. + do 25 ldummy=1,kp1mm + factor = fkp1mm/(t(il+kp1mm) - t(il)) +c the assumption that t(left) < t(left+1) makes denominator +c in factor nonzero. + do 24 j=1,i + 24 a(i,j) = (a(i,j) - a(i-1,j))*factor + il = il - 1 + 25 i = i - 1 +c +c for i=1,...,k, combine b-coeffs a(.,i) with b-spline values +c stored in dbiatx(.,m) to get value of (m-1)st derivative of +c i-th b-spline (of interest here) at x , and store in +c dbiatx(i,m). storage of this value over the value of a b-spline +c of order m there is safe since the remaining b-spline derivat- +c ive of the same order do not use this value due to the fact +c that a(j,i) = 0 for j < i . + do 40 i=1,k + sum = 0.d0 + jlow = max0(i,m) + do 35 j=jlow,k + 35 sum = a(j,i)*dbiatx(j,m) + sum + 40 dbiatx(i,m) = sum + 99 return + end diff --git a/bvalue.f b/bvalue.f new file mode 100644 index 0000000..3201c18 --- /dev/null +++ b/bvalue.f @@ -0,0 +1,135 @@ + real function bvalue ( t, bcoef, n, k, x, jderiv ) +c from * a practical guide to splines * by c. de boor +calls interv +c +calculates value at x of jderiv-th derivative of spline from b-repr. +c the spline is taken to be continuous from the right, EXCEPT at the +c rightmost knot, where it is taken to be continuous from the left. +c +c****** i n p u t ****** +c t, bcoef, n, k......forms the b-representation of the spline f to +c be evaluated. specifically, +c t.....knot sequence, of length n+k, assumed nondecreasing. +c bcoef.....b-coefficient sequence, of length n . +c n.....length of bcoef and dimension of spline(k,t), +c a s s u m e d positive . +c k.....order of the spline . +c +c w a r n i n g . . . the restriction k .le. kmax (=20) is imposed +c arbitrarily by the dimension statement for aj, dl, dr below, +c but is n o w h e r e c h e c k e d for. +c +c x.....the point at which to evaluate . +c jderiv.....integer giving the order of the derivative to be evaluated +c a s s u m e d to be zero or positive. +c +c****** o u t p u t ****** +c bvalue.....the value of the (jderiv)-th derivative of f at x . +c +c****** m e t h o d ****** +c The nontrivial knot interval (t(i),t(i+1)) containing x is lo- +c cated with the aid of interv . The k b-coeffs of f relevant for +c this interval are then obtained from bcoef (or taken to be zero if +c not explicitly available) and are then differenced jderiv times to +c obtain the b-coeffs of (d**jderiv)f relevant for that interval. +c Precisely, with j = jderiv, we have from x.(12) of the text that +c +c (d**j)f = sum ( bcoef(.,j)*b(.,k-j,t) ) +c +c where +c / bcoef(.), , j .eq. 0 +c / +c bcoef(.,j) = / bcoef(.,j-1) - bcoef(.-1,j-1) +c / ----------------------------- , j .gt. 0 +c / (t(.+k-j) - t(.))/(k-j) +c +c Then, we use repeatedly the fact that +c +c sum ( a(.)*b(.,m,t)(x) ) = sum ( a(.,x)*b(.,m-1,t)(x) ) +c with +c (x - t(.))*a(.) + (t(.+m-1) - x)*a(.-1) +c a(.,x) = --------------------------------------- +c (x - t(.)) + (t(.+m-1) - x) +c +c to write (d**j)f(x) eventually as a linear combination of b-splines +c of order 1 , and the coefficient for b(i,1,t)(x) must then be the +c desired number (d**j)f(x). (see x.(17)-(19) of text). +c + integer jderiv,k,n, i,ilo,imk,j,jc,jcmin,jcmax,jj,kmax,kmj,km1 + * ,mflag,nmi,jdrvp1 + parameter (kmax = 20) + real bcoef(n),t(n+k),x, aj(kmax),dl(kmax),dr(kmax),fkmj + bvalue = 0. + if (jderiv .ge. k) go to 99 +c +c *** Find i s.t. 1 .le. i .lt. n+k and t(i) .lt. t(i+1) and +c t(i) .le. x .lt. t(i+1) . If no such i can be found, x lies +c outside the support of the spline f , hence bvalue = 0. +c (The asymmetry in this choice of i makes f rightcontinuous, except +c at t(n+k) where it is leftcontinuous.) + call interv ( t, n+k, x, i, mflag ) + if (mflag .ne. 0) go to 99 +c *** if k = 1 (and jderiv = 0), bvalue = bcoef(i). + km1 = k - 1 + if (km1 .gt. 0) go to 1 + bvalue = bcoef(i) + go to 99 +c +c *** store the k b-spline coefficients relevant for the knot interval +c (t(i),t(i+1)) in aj(1),...,aj(k) and compute dl(j) = x - t(i+1-j), +c dr(j) = t(i+j) - x, j=1,...,k-1 . set any of the aj not obtainable +c from input to zero. set any t.s not obtainable equal to t(1) or +c to t(n+k) appropriately. + 1 jcmin = 1 + imk = i - k + if (imk .ge. 0) go to 8 + jcmin = 1 - imk + do 5 j=1,i + 5 dl(j) = x - t(i+1-j) + do 6 j=i,km1 + aj(k-j) = 0. + 6 dl(j) = dl(i) + go to 10 + 8 do 9 j=1,km1 + 9 dl(j) = x - t(i+1-j) +c + 10 jcmax = k + nmi = n - i + if (nmi .ge. 0) go to 18 + jcmax = k + nmi + do 15 j=1,jcmax + 15 dr(j) = t(i+j) - x + do 16 j=jcmax,km1 + aj(j+1) = 0. + 16 dr(j) = dr(jcmax) + go to 20 + 18 do 19 j=1,km1 + 19 dr(j) = t(i+j) - x +c + 20 do 21 jc=jcmin,jcmax + 21 aj(jc) = bcoef(imk + jc) +c +c *** difference the coefficients jderiv times. + if (jderiv .eq. 0) go to 30 + do 23 j=1,jderiv + kmj = k-j + fkmj = float(kmj) + ilo = kmj + do 23 jj=1,kmj + aj(jj) = ((aj(jj+1) - aj(jj))/(dl(ilo) + dr(jj)))*fkmj + 23 ilo = ilo - 1 +c +c *** compute value at x in (t(i),t(i+1)) of jderiv-th derivative, +c given its relevant b-spline coeffs in aj(1),...,aj(k-jderiv). + 30 if (jderiv .eq. km1) go to 39 + jdrvp1 = jderiv + 1 + do 33 j=jdrvp1,km1 + kmj = k-j + ilo = kmj + do 33 jj=1,kmj + aj(jj) = (aj(jj+1)*dl(ilo) + aj(jj)*dr(jj))/(dl(ilo)+dr(jj)) + 33 ilo = ilo - 1 + 39 bvalue = aj(1) +c + 99 return + end diff --git a/chimerauchimecommand.cpp b/chimerauchimecommand.cpp index 0846bbd..40b3f26 100644 --- a/chimerauchimecommand.cpp +++ b/chimerauchimecommand.cpp @@ -987,6 +987,11 @@ int ChimeraUchimeCommand::driverGroups(SequenceParser& parser, string outputFNam int ChimeraUchimeCommand::driver(string outputFName, string filename, string accnos, string alns, int& numChimeras){ try { + + outputFName = m->getFullPathName(outputFName); + filename = m->getFullPathName(filename); + alns = m->getFullPathName(alns); + //to allow for spaces in the path outputFName = "\"" + outputFName + "\""; filename = "\"" + filename + "\""; diff --git a/chimerauchimecommand.h b/chimerauchimecommand.h index f6e0ab5..45dd8a8 100644 --- a/chimerauchimecommand.h +++ b/chimerauchimecommand.h @@ -149,6 +149,10 @@ static DWORD WINAPI MyUchimeThreadFunction(LPVOID lpParam){ try { + pDataArray->outputFName = pDataArray->m->getFullPathName(outputFName); + pDataArray->filename = pDataArray->m->getFullPathName(filename); + pDataArray->alns = pDataArray->m->getFullPathName(alns); + //clears files ofstream out, out1, out2; pDataArray->m->openOutputFile(pDataArray->outputFName, out); out.close(); @@ -494,6 +498,10 @@ static DWORD WINAPI MyUchimeSeqsThreadFunction(LPVOID lpParam){ try { + pDataArray->outputFName = pDataArray->m->getFullPathName(outputFName); + pDataArray->filename = pDataArray->m->getFullPathName(filename); + pDataArray->alns = pDataArray->m->getFullPathName(alns); + int totalSeqs = 0; int numChimeras = 0; diff --git a/daxpy.f b/daxpy.f new file mode 100644 index 0000000..ddc7449 --- /dev/null +++ b/daxpy.f @@ -0,0 +1,69 @@ + SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) +* .. Scalar Arguments .. + DOUBLE PRECISION DA + INTEGER INCX,INCY,N +* .. +* .. Array Arguments .. + DOUBLE PRECISION DX(*),DY(*) +* .. +* +* Purpose +* ======= +* +* DAXPY constant times a vector plus a vector. +* uses unrolled loops for increments equal to one. +* +* Further Details +* =============== +* +* jack dongarra, linpack, 3/11/78. +* modified 12/3/93, array(1) declarations changed to array(*) +* +* ===================================================================== +* +* .. Local Scalars .. + INTEGER I,IX,IY,M,MP1 +* .. +* .. Intrinsic Functions .. + INTRINSIC MOD +* .. + IF (N.LE.0) RETURN + IF (DA.EQ.0.0d0) RETURN + IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN +* +* code for both increments equal to 1 +* +* +* clean-up loop +* + M = MOD(N,4) + IF (M.NE.0) THEN + DO I = 1,M + DY(I) = DY(I) + DA*DX(I) + END DO + END IF + IF (N.LT.4) RETURN + MP1 = M + 1 + DO I = MP1,N,4 + DY(I) = DY(I) + DA*DX(I) + DY(I+1) = DY(I+1) + DA*DX(I+1) + DY(I+2) = DY(I+2) + DA*DX(I+2) + DY(I+3) = DY(I+3) + DA*DX(I+3) + END DO + ELSE +* +* code for unequal increments or equal increments +* not equal to 1 +* + IX = 1 + IY = 1 + IF (INCX.LT.0) IX = (-N+1)*INCX + 1 + IF (INCY.LT.0) IY = (-N+1)*INCY + 1 + DO I = 1,N + DY(IY) = DY(IY) + DA*DX(IX) + IX = IX + INCX + IY = IY + INCY + END DO + END IF + RETURN + END diff --git a/ddot.f b/ddot.f new file mode 100644 index 0000000..373134f --- /dev/null +++ b/ddot.f @@ -0,0 +1,90 @@ +*DECK DDOT + DOUBLE PRECISION FUNCTION DDOT (N, DX, INCX, DY, INCY) +C***BEGIN PROLOGUE DDOT +C***PURPOSE Compute the inner product of two vectors. +C***LIBRARY SLATEC (BLAS) +C***CATEGORY D1A4 +C***TYPE DOUBLE PRECISION (SDOT-S, DDOT-D, CDOTU-C) +C***KEYWORDS BLAS, INNER PRODUCT, LINEAR ALGEBRA, VECTOR +C***AUTHOR Lawson, C. L., (JPL) +C Hanson, R. J., (SNLA) +C Kincaid, D. R., (U. of Texas) +C Krogh, F. T., (JPL) +C***DESCRIPTION +C +C B L A S Subprogram +C Description of Parameters +C +C --Input-- +C N number of elements in input vector(s) +C DX double precision vector with N elements +C INCX storage spacing between elements of DX +C DY double precision vector with N elements +C INCY storage spacing between elements of DY +C +C --Output-- +C DDOT double precision dot product (zero if N .LE. 0) +C +C Returns the dot product of double precision DX and DY. +C DDOT = sum for I = 0 to N-1 of DX(LX+I*INCX) * DY(LY+I*INCY), +C where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is +C defined in a similar way using INCY. +C +C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T. +C Krogh, Basic linear algebra subprograms for Fortran +C usage, Algorithm No. 539, Transactions on Mathematical +C Software 5, 3 (September 1979), pp. 308-323. +C***ROUTINES CALLED (NONE) +C***REVISION HISTORY (YYMMDD) +C 791001 DATE WRITTEN +C 890831 Modified array declarations. (WRB) +C 890831 REVISION DATE from Version 3.2 +C 891214 Prologue converted to Version 4.0 format. (BAB) +C 920310 Corrected definition of LX in DESCRIPTION. (WRB) +C 920501 Reformatted the REFERENCES section. (WRB) +C***END PROLOGUE DDOT + DOUBLE PRECISION DX(*), DY(*) +C***FIRST EXECUTABLE STATEMENT DDOT + DDOT = 0.0D0 + DTEMP = 0.0D0 + IF (N .LE. 0) RETURN + IF (INCX .EQ. INCY) IF (INCX-1) 5,20,60 +C +C Code for unequal or nonpositive increments. +C + 5 IX = 1 + IY = 1 + IF (INCX .LT. 0) IX = (-N+1)*INCX + 1 + IF (INCY .LT. 0) IY = (-N+1)*INCY + 1 + DO 10 I = 1,N + DDOT = DDOT + DX(IX)*DY(IY) + IX = IX + INCX + IY = IY + INCY + 10 CONTINUE + RETURN +C +C Code for both increments equal to 1. +C +C Clean-up loop so remaining vector length is a multiple of 5. +C + 20 M = MOD(N,5) + IF (M .EQ. 0) GO TO 40 + DO 30 I = 1,M + DDOT = DDOT + DX(I)*DY(I) + 30 CONTINUE + IF (N .LT. 5) RETURN + 40 MP1 = M + 1 + DO I = MP1,N,5 + DDOT = DDOT + DX(I)*DY(I) + DX(I+1)*DY(I+1) + DX(I+2)*DY(I+2) + + 1 DX(I+3)*DY(I+3) + DX(I+4)*DY(I+4) + END DO + RETURN +C +C Code for equal, positive, non-unit increments. +C + 60 NS = N*INCX + DO 70 I = 1,NS,INCX + DDOT = DDOT + DX(I)*DY(I) + 70 CONTINUE + RETURN + END diff --git a/dpbfa.f b/dpbfa.f new file mode 100644 index 0000000..95952c0 --- /dev/null +++ b/dpbfa.f @@ -0,0 +1,109 @@ + subroutine dpbfa(abd,lda,n,m,info) + + integer lda,n,m,info + double precision abd(lda,n) +c +c dpbfa factors a double precision symmetric positive definite +c matrix stored in band form. +c +c dpbfa is usually called by dpbco, but it can be called +c directly with a saving in time if rcond is not needed. +c +c on entry +c +c abd double precision(lda, n) +c the matrix to be factored. the columns of the upper +c triangle are stored in the columns of abd and the +c diagonals of the upper triangle are stored in the +c rows of abd . see the comments below for details. +c +c lda integer +c the leading dimension of the array abd . +c lda must be .ge. m + 1 . +c +c n integer +c the order of the matrix a . +c +c m integer +c the number of diagonals above the main diagonal. +c 0 .le. m .lt. n . +c +c on return +c +c abd an upper triangular matrix r , stored in band +c form, so that a = trans(r)*r . +c +c info integer +c = 0 for normal return. +c = k if the leading minor of order k is not +c positive definite. +c +c band storage +c +c if a is a symmetric positive definite band matrix, +c the following program segment will set up the input. +c +c m = (band width above diagonal) +c do 20 j = 1, n +c i1 = max0(1, j-m) +c do 10 i = i1, j +c k = i-j+m+1 +c abd(k,j) = a(i,j) +c 10 continue +c 20 continue +c +c linpack. this version dated 08/14/78 . +c cleve moler, university of new mexico, argonne national lab. +c +c subroutines and functions +c +c blas ddot +c fortran max0,sqrt +c +c internal variables +c + double precision ddot,t, a(10), b(10), temp + double precision s + integer ik,j,jk,k,mu +c begin block with ...exits to 40 +c +c + do i = 1, 10 + a(i) = i + b(i) = 2*i +c PRINT *, 'a = ', a(i), 'i = ', i +c PRINT *, 'b = ', b(i) + end do + + temp = ddot(10,a,1,b,1) + + do 30 j = 1, n + info = j + s = 0.0d0 + ik = m + 1 + jk = max0(j-m,1) + mu = max0(m+2-j,1) + if (m .lt. mu) go to 20 + do 10 k = mu, m + + t = abd(k,j) - ddot(k-mu,abd(ik,jk),1,abd(mu,j),1) + t = t/abd(m+1,jk) + abd(k,j) = t + s = s + t*t + ik = ik - 1 + jk = jk + 1 + 10 continue + 20 continue + + s = abd(m+1,j) - s +c ......exit + if (s .le. 0.0d0) go to 40 + + abd(m+1,j) = sqrt(s) + + 30 continue + info = 0 + 40 continue + return + end + diff --git a/dpbsl.f b/dpbsl.f new file mode 100644 index 0000000..d910dee --- /dev/null +++ b/dpbsl.f @@ -0,0 +1,83 @@ + subroutine dpbsl(abd,lda,n,m,b) + + integer lda,n,m + double precision abd(lda,n),b(n) +c +c dpbsl solves the double precision symmetric positive definite +c band system a*x = b +c using the factors computed by dpbco or dpbfa. +c +c on entry +c +c abd double precision(lda, n) +c the output from dpbco or dpbfa. +c +c lda integer +c the leading dimension of the array abd . +c +c n integer +c the order of the matrix a . +c +c m integer +c the number of diagonals above the main diagonal. +c +c b double precision(n) +c the right hand side vector. +c +c on return +c +c b the solution vector x . +c +c error condition +c +c a division by zero will occur if the input factor contains +c a zero on the diagonal. technically this indicates +c singularity but it is usually caused by improper subroutine +c arguments. it will not occur if the subroutines are called +c correctly and info .eq. 0 . +c +c to compute inverse(a) * c where c is a matrix +c with p columns +c call dpbco(abd,lda,n,rcond,z,info) +c if (rcond is too small .or. info .ne. 0) go to ... +c do 10 j = 1, p +c call dpbsl(abd,lda,n,c(1,j)) +c 10 continue +c +c linpack. this version dated 08/14/78 . +c cleve moler, university of new mexico, argonne national lab. +c +c subroutines and functions +c +c blas daxpy,ddot +c fortran min0 +c +c internal variables +c + double precision ddot,t + integer k,kb,la,lb,lm +c +c solve trans(r)*y = b +c + do 10 k = 1, n + lm = min0(k-1,m) + la = m + 1 - lm + lb = k - lm + t = ddot(lm,abd(la,k),1,b(lb),1) + b(k) = (b(k) - t)/abd(m+1,k) + 10 continue +c +c solve r*x = y +c + do 20 kb = 1, n + k = n + 1 - kb + lm = min0(k-1,m) + la = m + 1 - lm + lb = k - lm + b(k) = b(k)/abd(m+1,k) + t = -b(k) + call daxpy(lm,t,abd(la,k),1,b(lb),1) + 20 continue + return + end + diff --git a/interv.f b/interv.f new file mode 100644 index 0000000..d65a914 --- /dev/null +++ b/interv.f @@ -0,0 +1,102 @@ + subroutine interv ( xt, lxt, x, left, mflag ) +c from * a practical guide to splines * by C. de Boor +computes left = max( i : xt(i) .lt. xt(lxt) .and. xt(i) .le. x ) . +c +c****** i n p u t ****** +c xt.....a real sequence, of length lxt , assumed to be nondecreasing +c lxt.....number of terms in the sequence xt . +c x.....the point whose location with respect to the sequence xt is +c to be determined. +c +c****** o u t p u t ****** +c left, mflag.....both integers, whose value is +c +c 1 -1 if x .lt. xt(1) +c i 0 if xt(i) .le. x .lt. xt(i+1) +c i 0 if xt(i) .lt. x .eq. xt(i+1) .eq. xt(lxt) +c i 1 if xt(i) .lt. xt(i+1) .eq. xt(lxt) .lt. x +c +c In particular, mflag = 0 is the 'usual' case. mflag .ne. 0 +c indicates that x lies outside the CLOSED interval +c xt(1) .le. y .le. xt(lxt) . The asymmetric treatment of the +c intervals is due to the decision to make all pp functions cont- +c inuous from the right, but, by returning mflag = 0 even if +C x = xt(lxt), there is the option of having the computed pp function +c continuous from the left at xt(lxt) . +c +c****** m e t h o d ****** +c The program is designed to be efficient in the common situation that +c it is called repeatedly, with x taken from an increasing or decrea- +c sing sequence. This will happen, e.g., when a pp function is to be +c graphed. The first guess for left is therefore taken to be the val- +c ue returned at the previous call and stored in the l o c a l varia- +c ble ilo . A first check ascertains that ilo .lt. lxt (this is nec- +c essary since the present call may have nothing to do with the previ- +c ous call). Then, if xt(ilo) .le. x .lt. xt(ilo+1), we set left = +c ilo and are done after just three comparisons. +c Otherwise, we repeatedly double the difference istep = ihi - ilo +c while also moving ilo and ihi in the direction of x , until +c xt(ilo) .le. x .lt. xt(ihi) , +c after which we use bisection to get, in addition, ilo+1 = ihi . +c left = ilo is then returned. +c + integer left,lxt,mflag, ihi,ilo,istep,middle + double precision x,xt(lxt) + data ilo /1/ + save ilo + ihi = ilo + 1 + if (ihi .lt. lxt) go to 20 + if (x .ge. xt(lxt)) go to 110 + if (lxt .le. 1) go to 90 + ilo = lxt - 1 + ihi = lxt +c + 20 if (x .ge. xt(ihi)) go to 40 + if (x .ge. xt(ilo)) go to 100 +c +c **** now x .lt. xt(ilo) . decrease ilo to capture x . + istep = 1 + 31 ihi = ilo + ilo = ihi - istep + if (ilo .le. 1) go to 35 + if (x .ge. xt(ilo)) go to 50 + istep = istep*2 + go to 31 + 35 ilo = 1 + if (x .lt. xt(1)) go to 90 + go to 50 +c **** now x .ge. xt(ihi) . increase ihi to capture x . + 40 istep = 1 + 41 ilo = ihi + ihi = ilo + istep + if (ihi .ge. lxt) go to 45 + if (x .lt. xt(ihi)) go to 50 + istep = istep*2 + go to 41 + 45 if (x .ge. xt(lxt)) go to 110 + ihi = lxt +c +c **** now xt(ilo) .le. x .lt. xt(ihi) . narrow the interval. + 50 middle = (ilo + ihi)/2 + if (middle .eq. ilo) go to 100 +c note. it is assumed that middle = ilo in case ihi = ilo+1 . + if (x .lt. xt(middle)) go to 53 + ilo = middle + go to 50 + 53 ihi = middle + go to 50 +c**** set output and return. + 90 mflag = -1 + left = 1 + return + 100 mflag = 0 + left = ilo + return + 110 mflag = 1 + if (x .eq. xt(lxt)) mflag = 0 + left = lxt + 111 if (left .eq. 1) return + left = left - 1 + if (xt(left) .lt. xt(lxt)) return + go to 111 + end diff --git a/intrv.f b/intrv.f new file mode 100644 index 0000000..827b5bd --- /dev/null +++ b/intrv.f @@ -0,0 +1,117 @@ +*DECK INTRV + SUBROUTINE INTRV (XT, LXT, X, ILO, ILEFT, MFLAG) +C***BEGIN PROLOGUE INTRV +C***PURPOSE Compute the largest integer ILEFT in 1 .LE. ILEFT .LE. LXT +C such that XT(ILEFT) .LE. X where XT(*) is a subdivision +C of the X interval. +C***LIBRARY SLATEC +C***CATEGORY E3, K6 +C***TYPE SINGLE PRECISION (INTRV-S, DINTRV-D) +C***KEYWORDS B-SPLINE, DATA FITTING, INTERPOLATION, SPLINES +C***AUTHOR Amos, D. E., (SNLA) +C***DESCRIPTION +C +C Written by Carl de Boor and modified by D. E. Amos +C +C Abstract +C INTRV is the INTERV routine of the reference. +C +C INTRV computes the largest integer ILEFT in 1 .LE. ILEFT .LE. +C LXT such that XT(ILEFT) .LE. X where XT(*) is a subdivision of +C the X interval. Precisely, +C +C X .LT. XT(1) 1 -1 +C if XT(I) .LE. X .LT. XT(I+1) then ILEFT=I , MFLAG=0 +C XT(LXT) .LE. X LXT 1, +C +C That is, when multiplicities are present in the break point +C to the left of X, the largest index is taken for ILEFT. +C +C Description of Arguments +C Input +C XT - XT is a knot or break point vector of length LXT +C LXT - length of the XT vector +C X - argument +C ILO - an initialization parameter which must be set +C to 1 the first time the spline array XT is +C processed by INTRV. +C +C Output +C ILO - ILO contains information for efficient process- +C ing after the initial call, and ILO must not be +C changed by the user. Distinct splines require +C distinct ILO parameters. +C ILEFT - largest integer satisfying XT(ILEFT) .LE. X +C MFLAG - signals when X lies out of bounds +C +C Error Conditions +C None +C +C***REFERENCES Carl de Boor, Package for calculating with B-splines, +C SIAM Journal on Numerical Analysis 14, 3 (June 1977), +C pp. 441-472. +C***ROUTINES CALLED (NONE) +C***REVISION HISTORY (YYMMDD) +C 800901 DATE WRITTEN +C 890831 Modified array declarations. (WRB) +C 890831 REVISION DATE from Version 3.2 +C 891214 Prologue converted to Version 4.0 format. (BAB) +C 920501 Reformatted the REFERENCES section. (WRB) +C***END PROLOGUE INTRV +C + INTEGER IHI, ILEFT, ILO, ISTEP, LXT, MFLAG, MIDDLE + Double precision X, XT + DIMENSION XT(*) +C***FIRST EXECUTABLE STATEMENT INTRV + IHI = ILO + 1 + IF (IHI.LT.LXT) GO TO 10 + IF (X.GE.XT(LXT)) GO TO 110 + IF (LXT.LE.1) GO TO 90 + ILO = LXT - 1 + IHI = LXT +C + 10 IF (X.GE.XT(IHI)) GO TO 40 + IF (X.GE.XT(ILO)) GO TO 100 +C +C *** NOW X .LT. XT(IHI) . FIND LOWER BOUND + ISTEP = 1 + 20 IHI = ILO + ILO = IHI - ISTEP + IF (ILO.LE.1) GO TO 30 + IF (X.GE.XT(ILO)) GO TO 70 + ISTEP = ISTEP*2 + GO TO 20 + 30 ILO = 1 + IF (X.LT.XT(1)) GO TO 90 + GO TO 70 +C *** NOW X .GE. XT(ILO) . FIND UPPER BOUND + 40 ISTEP = 1 + 50 ILO = IHI + IHI = ILO + ISTEP + IF (IHI.GE.LXT) GO TO 60 + IF (X.LT.XT(IHI)) GO TO 70 + ISTEP = ISTEP*2 + GO TO 50 + 60 IF (X.GE.XT(LXT)) GO TO 110 + IHI = LXT +C +C *** NOW XT(ILO) .LE. X .LT. XT(IHI) . NARROW THE INTERVAL + 70 MIDDLE = (ILO+IHI)/2 + IF (MIDDLE.EQ.ILO) GO TO 100 +C NOTE. IT IS ASSUMED THAT MIDDLE = ILO IN CASE IHI = ILO+1 + IF (X.LT.XT(MIDDLE)) GO TO 80 + ILO = MIDDLE + GO TO 70 + 80 IHI = MIDDLE + GO TO 70 +C *** SET OUTPUT AND RETURN + 90 MFLAG = -1 + ILEFT = 1 + RETURN + 100 MFLAG = 0 + ILEFT = ILO + RETURN + 110 MFLAG = 1 + ILEFT = LXT + RETURN + END diff --git a/makefile b/makefile index 44d8a81..90e4957 100644 --- a/makefile +++ b/makefile @@ -15,11 +15,12 @@ USEREADLINE ?= yes CYGWIN_BUILD ?= no USECOMPRESSION ?= no MOTHUR_FILES="\"Enter_your_default_path_here\"" -RELEASE_DATE = "\"10/18/2011\"" +RELEASE_DATE = "\"12/7/2011\"" VERSION = "\"1.22.0\"" +FORTAN_COMPILER = gfortran # Optimize to level 3: -CXXFLAGS += -O3 +CXXFLAGS += -O3 ifeq ($(strip $(64BIT_VERSION)),yes) #if you are using centos uncomment the following lines @@ -89,14 +90,18 @@ endif OBJECTS=$(patsubst %.cpp,%.o,$(wildcard *.cpp)) OBJECTS+=$(patsubst %.c,%.o,$(wildcard *.c)) +OBJECTS+=$(patsubst %.f,%.o,$(wildcard *.f)) -mothur : $(OBJECTS) uchime +mothur : fortranSource $(OBJECTS) uchime $(CXX) $(LDFLAGS) $(TARGET_ARCH) -o $@ $(OBJECTS) $(LIBS) strip mothur uchime: cd uchime_src && ./mk && mv uchime .. && cd .. + +fortranSource: + ${FORTAN_COMPILER} -c *.f install : mothur # cp mothur ../Release/mothur diff --git a/mothur.h b/mothur.h index 57ece67..3b7c459 100644 --- a/mothur.h +++ b/mothur.h @@ -189,7 +189,8 @@ void convert(const string& s, T& x, bool failIfLeftoverChars = true){ throw BadConversion(s); } - +//********************************************************************************************************************** +template int sgn(T val){ return (val > T(0)) - (val < T(0)); } //********************************************************************************************************************** template diff --git a/mothurmetastats.cpp b/mothurmetastats.cpp index a4f879e..ae82225 100644 --- a/mothurmetastats.cpp +++ b/mothurmetastats.cpp @@ -9,6 +9,7 @@ #include "mothurmetastats.h" #include "mothurfisher.h" +#include "spline.h" /***********************************************************/ MothurMetastats::MothurMetastats(double t, int n) { @@ -204,13 +205,8 @@ int MothurMetastats::runMetastats(string outputFileName, vector< vector storage[i][7]=temp[i][1]; storage[i][8]=pvalues[i]; } - cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint); - cout << "pvalues" << endl; - for (int i = 0; i < row; i++){ if (pvalues[i] < 0.0000001) { - pvalues[i] = 0.0; - }cout << pvalues[i] << '\t'; } - cout << endl; - calc_qvalues(pvalues); + + vector qvalues = calc_qvalues(pvalues); // BACKUP checks cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint); @@ -218,9 +214,9 @@ int MothurMetastats::runMetastats(string outputFileName, vector< vector if (m->control_pressed) { return 1; } - if(pvalues[i] < threshold){ - m->mothurOut("Feature " + toString((i+1)) + " is significant, p = "); - cout << pvalues[i]; + if(qvalues[i] < threshold){ + m->mothurOut("Feature " + toString((i+1)) + " is significant, q = "); + cout << qvalues[i]; m->mothurOutJustToLog(toString(pvalues[i])); m->mothurOutEndLine(); } } @@ -240,13 +236,14 @@ int MothurMetastats::runMetastats(string outputFileName, vector< vector //output column headings - not really sure... documentation labels 9 columns, there are 10 in the output file //storage 0 = meanGroup1 - line 529, 1 = varGroup1 - line 532, 2 = err rate1 - line 534, 3 = mean of counts group1?? - line 291, 4 = meanGroup2 - line 536, 5 = varGroup2 - line 539, 6 = err rate2 - line 541, 7 = mean of counts group2?? - line 292, 8 = pvalues - line 293 - out << "OTU\tmean(group1)\tvariance(group1)\tstderr(group1)\tmean_of_counts(group1)\tmean(group2)\tvariance(group2)\tstderr(group2)\tmean_of_counts(group1)\tp-value\n"; + out << "OTU\tmean(group1)\tvariance(group1)\tstderr(group1)\tmean_of_counts(group1)\tmean(group2)\tvariance(group2)\tstderr(group2)\tmean_of_counts(group1)\tp-value\tq-value\n"; for(int i = 0; i < row; i++){ if (m->control_pressed) { out.close(); return 0; } out << (i+1); for(int j = 0; j < 9; j++){ out << '\t' << storage[i][j]; } + out << '\t' << qvalues[i]; out << endl; } @@ -494,9 +491,10 @@ int MothurMetastats::calc_twosample_ts(vector& Pmatrix, int secondGroupi /***********************************************************/ vector MothurMetastats::calc_qvalues(vector& pValues) { try { - vector qvalues; - int numRows = pValues.size(); + int numRows = pValues.size(); + vector qvalues(numRows, 0.0); + //fill lambdas with 0.00, 0.01, 0.02... 0.95 vector lambdas(96, 0); for (int i = 1; i < lambdas.size(); i++) { lambdas[i] = lambdas[i-1] + 0.01; } @@ -508,21 +506,11 @@ vector MothurMetastats::calc_qvalues(vector& pValues) { int count = 0; for (int i = 0; i < numRows; i++){ // for each p-value in order if (pValues[i] > lambdas[l]){ count++; } - pi0_hat[l] = count/(double)(numRows*(1-lambdas[l])); } + pi0_hat[l] = count/(double)(numRows*(1-lambdas[l])); } - //from R code - replacing with spline and splint below - //vector f_spline = smoothSpline(lambdas, pi0_hat, 3); - //double pi0 = f_spline[(f_spline.size()-1)]; // this is the essential pi0_hat value - - //cubic Spline - double pi0, notsure; //this is the essential pi0_hat value - int notSure; - vector resultsSpline(lambdas.size(), 0.0); - spline(pi0_hat, lambdas, notSure, notSure, resultsSpline); - //some sort of loop to get best value?? - splint(pi0_hat, lambdas, notsure, pi0, resultsSpline); + double pi0 = smoothSpline(lambdas, pi0_hat, 3); //order p-values vector ordered_qs = qvalues; @@ -530,14 +518,13 @@ vector MothurMetastats::calc_qvalues(vector& pValues) { for (int i = 1; i < ordered_ps.size(); i++) { ordered_ps[i] = ordered_ps[i-1] + 1; } vector tempPvalues = pValues; OrderPValues(0, numRows-1, tempPvalues, ordered_ps); - - ordered_qs[numRows-1] <- min((pValues[ordered_ps[numRows-1]]*pi0), 1.0); + + ordered_qs[numRows-1] = min((pValues[ordered_ps[numRows-1]]*pi0), 1.0); for (int i = (numRows-2); i >= 0; i--){ double p = pValues[ordered_ps[i]]; - double temp = p*numRows*pi0/(double)i; - + double temp = p*numRows*pi0/(double)(i+1); + ordered_qs[i] = min(temp, ordered_qs[i+1]); - ordered_qs[i] = min(ordered_qs[i], 1.0); } //re-distribute calculated qvalues to appropriate rows @@ -612,41 +599,54 @@ int MothurMetastats::swapElements(int i, int j, vector& p, vector& exit(1); } } -/*********************************************************** -vector MothurMetastats::smoothSpline(vector x, vector y, int df) { +/***********************************************************/ +double MothurMetastats::smoothSpline(vector& x, vector& y, int df) { try { - - cout << "lambdas" << endl; - for (int l = 0; l < x.size(); l++){ cout << x[l] << '\t'; } - cout << endl << "pi0_hat" << endl; - for (int l = 0; l < y.size(); l++){ cout << y[l] << '\t'; } - cout << endl; - - //double low = -1.5; double high = 1.5; double tol = 1e-04; double eps = 2e-08; double maxit = 500; + + double result = 0.0; int n = x.size(); vector w(n, 1); - - //x <- signif(x, 6L) - I think this rounds to 6 decimals places - //ux <- unique(sort(x)) //x will be unique and sorted since we created it - //ox <- match(x, ux) //since its unique and sort it will match - - vector wbar(n, 0); - vector ybar(n, 0); - vector xbar(n, 0.0); + double* xb = new double[n]; + double* yb = new double[n]; + double* wb = new double[n]; double yssw = 0.0; for (int i = 0; i < n; i++) { - wbar[i] = w[i]; - ybar[i] = w[i]*y[i]; - yssw += (w[i] * y[i] * y[i]) - wbar[i] * (ybar[i] * ybar[i]); - xbar[i] = (x[i] - x[0]) / (x[n-1] - x[0]); + wb[i] = w[i]; + yb[i] = w[i]*y[i]; + yssw += (w[i] * y[i] * y[i]) - wb[i] * (yb[i] * yb[i]); + xb[i] = (x[i] - x[0]) / (x[n-1] - x[0]); } - vector knot = sknot1(xbar); + vector knot = sknot1(xb, n); int nk = knot.size() - 4; - - //double ispar = 0.0; double spar = 0.0; double icrit = 3.0; double dofoff = 3.0; - - return y; + + double low = -1.5; double high = 1.5; double tol = 1e-04; double eps = 2e-08; int maxit = 500; + int ispar = 0; int icrit = 3; double dofoff = 3.0; + double penalty = 1.0; + int ld4 = 4; int isetup = 0; int ldnk = 1; int ier; double spar = 1.0; double crit; + + double* knotb = new double[knot.size()]; + double* coef1 = new double[nk]; + double* lev1 = new double[n]; + double* sz1 = new double[n]; + for (int i = 0; i < knot.size(); i++) { knotb[i] = knot[i]; } + + Spline spline; + spline.sbart(&penalty, &dofoff, &xb[0], &yb[0], &wb[0], &yssw, &n, &knotb[0], &nk, &coef1[0], &sz1[0], &lev1[0], &crit, + &icrit, &spar, &ispar, &maxit, &low, &high, &tol, &eps, &isetup, &ld4, &ldnk, &ier); + + result = coef1[nk-1]; + + //free memory + delete [] xb; + delete [] yb; + delete [] wb; + delete [] knotb; + delete [] coef1; + delete [] lev1; + delete [] sz1; + + return result; }catch(exception& e) { m->errorOut(e, "MothurMetastats", "smoothSpline"); @@ -654,105 +654,21 @@ vector MothurMetastats::smoothSpline(vector x, vector y, } } /***********************************************************/ -// This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479 -int MothurMetastats::spline(vector& x, vector& y, int yp1, int ypn, vector& y2) { - try { - - cout << "lambdas" << endl; - for (int l = 0; l < x.size(); l++){ cout << x[l] << '\t'; } - cout << endl << "pi0_hat" << endl; - for (int l = 0; l < y.size(); l++){ cout << y[l] << '\t'; } - cout << endl; - - double p, qn, sig, un; - - int n = y2.size(); - vector u(n-1, 0.0); - - if (yp1 > 0.99e30) { y2[0] = u[0] = 0.0; } - else { - y2[0] = -0.5; - u[0] = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1]-x[0]) - yp1); - } - - for (int i = 1; i < n-1; i++) { - sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]); - p = sig * y2[i-1] + 2.0; - y2[i] = (sig - 1.0)/p; - u[i] = (y[i+1]-y[i]) / (x[i+1]-x[i]) - (y[i] - y[i-1]) / (x[i]-x[i-1]); - u[i] = (6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p; - } - - if (ypn > 0.99e30) { qn=un=0.0; } - else { - qn=0.5; - un=(3.0/(x[n-1]-x[n-2]))*(ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2])); - } - - y2[n-1]=(un-qn*u[n-2])/(qn*y2[n-2]+1.0); - - for (int k=n-2; k>=0; k--) { - y2[k]=y2[k]*y2[k+1]+u[k]; - } - - return 0; - - }catch(exception& e) { - m->errorOut(e, "MothurMetastats", "spline"); - exit(1); - } -} -/***********************************************************/ -// This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479 -int MothurMetastats::splint(vector& xa, vector& ya, double& x, double& y, vector& y2a) { - try { - int k; - double h,b,a; - - int n = xa.size(); - - int klo=0; - int khi=n-1; - while (khi-klo > 1) { - - if (m->control_pressed) { break; } - - k = (khi+klo) >> 1; - if (xa[k] > x) { khi=k; } - else { klo=k; } - } - - h=xa[khi]-xa[klo]; - if (h == 0.0) { m->mothurOut("[ERROR]: Bad xa input to splint routine."); m->mothurOutEndLine(); m->control_pressed = true; return 0; } - a=(xa[khi]-x)/h; - b=(x - xa[klo])/h; - y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0; - - return 0; - - }catch(exception& e) { - m->errorOut(e, "MothurMetastats", "splint"); - exit(1); - } -} -/***********************************************************/ -vector MothurMetastats::sknot1(vector x) { +vector MothurMetastats::sknot1(double* x, int n) { try { vector knots; - int n = x.size(); int nk = nkn(n); - cout << nk << endl; //R equivalent - rep(x[1L], 3L) knots.push_back(x[0]); knots.push_back(x[0]); knots.push_back(x[0]); //generate a sequence of nk equally spaced values from 1 to n. R equivalent = seq.int(1, n, length.out = nk) vector indexes = getSequence(0, n-1, nk); - for (int i = 0; i < indexes.size(); i++) { knots.push_back(x[indexes[i]]); } + for (int i = 0; i < indexes.size(); i++) { knots.push_back(x[indexes[i]]); } //R equivalent - rep(x[n], 3L) knots.push_back(x[n-1]); knots.push_back(x[n-1]); knots.push_back(x[n-1]); - + return knots; }catch(exception& e) { diff --git a/mothurmetastats.h b/mothurmetastats.h index 229ffb6..d80b687 100644 --- a/mothurmetastats.h +++ b/mothurmetastats.h @@ -32,18 +32,14 @@ class MothurMetastats { int permute_matrix(vector&, vector&, int, vector&, vector&, vector&); int permute_array(vector&); int calc_twosample_ts(vector&, int, vector&, vector&, vector&); - vector smoothSpline(vector, vector, int); + double smoothSpline(vector&, vector&, int); vector calc_qvalues(vector&); - vector sknot1(vector); + vector sknot1(double*, int); int nkn(int); int OrderPValues(int, int, vector&, vector&); int swapElements(int, int, vector&, vector&); vector getSequence(int, int, int); - int spline(vector&, vector&, int, int, vector&); - int splint(vector&, vector&, double&, double&, vector&); - - }; #endif diff --git a/phylotree.cpp b/phylotree.cpp index 5438147..e430fb9 100644 --- a/phylotree.cpp +++ b/phylotree.cpp @@ -241,6 +241,8 @@ int PhyloTree::addSeqToTree(string seqName, string seqTaxonomy){ int level = 1; tree[0].accessions.push_back(seqName); + m->removeConfidences(seqTaxonomy); + string taxon;// = getNextTaxon(seqTaxonomy); while(seqTaxonomy != ""){ diff --git a/sgram.f b/sgram.f new file mode 100644 index 0000000..a8cb6de --- /dev/null +++ b/sgram.f @@ -0,0 +1,142 @@ +C Output from Public domain Ratfor, version 1.0 + +C PURPOSE +C Calculation of the cubic B-spline smoothness prior +C for "usual" interior knot setup. +C Uses BSPVD and INTRV in the CMLIB +C sgm[0-3](nb) Symmetric matrix +C whose (i,j)'th element contains the integral of +C B''(i,.) B''(j,.) , i=1,2 ... nb and j=i,...nb. +C Only the upper four diagonals are computed. + + subroutine sgram(sg0,sg1,sg2,sg3,tb,nb) + +c implicit none +C indices + integer nb + DOUBLE precision sg0(nb),sg1(nb),sg2(nb),sg3(nb), tb(nb+4) +c ------------- + integer ileft,mflag, i,ii,jj, lentb + DOUBLE precision vnikx(4,3),work(16),yw1(4),yw2(4), wpt +c +c integer interv +c external interv + + lentb=nb+4 +C Initialise the sigma vectors + do 1 i=1,nb + sg0(i)=0.d0 + sg1(i)=0.d0 + sg2(i)=0.d0 + sg3(i)=0.d0 + 1 continue + + ileft = 1 + do 2 i=1,nb +C Calculate a linear approximation to the +C second derivative of the non-zero B-splines +C over the interval [tb(i),tb(i+1)]. +C call intrv(tb(1),(nb+1),tb(i),ilo,ileft,mflag) + call interv(tb(1), nb+1,tb(i), ileft, mflag) +C Left end second derivatives +C call bspvd (tb,4,3,tb(i),ileft,4,vnikx,work) + call bsplvd (tb,lentb,4,tb(i),ileft,work,vnikx,3) +C Put values into yw1 + do 4 ii=1,4 + yw1(ii) = vnikx(ii,3) + 4 continue + +C Right end second derivatives +C call bspvd (tb,4,3,tb(i+1),ileft,4,vnikx,work) + call bsplvd (tb,lentb,4,tb(i+1),ileft,work,vnikx,3) + +C Slope*(length of interval) in Linear Approximation to B'' + do 6 ii=1,4 + yw2(ii) = vnikx(ii,3) - yw1(ii) + 6 continue + + wpt = tb(i+1) - tb(i) +C Calculate Contributions to the sigma vectors + if(ileft.ge.4) then + do 10 ii=1,4 + jj=ii + sg0(ileft-4+ii) = sg0(ileft-4+ii) + + & wpt*(yw1(ii)*yw1(jj)+ + & (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 + & + yw2(ii)*yw2(jj)*0.3330d0) + jj=ii+1 + if(jj.le.4)then + sg1(ileft+ii-4) = sg1(ileft+ii-4) + + & wpt* (yw1(ii)*yw1(jj) + + * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 + & +yw2(ii)*yw2(jj)*0.3330d0 ) + endif + jj=ii+2 + if(jj.le.4)then + sg2(ileft+ii-4) = sg2(ileft+ii-4) + + & wpt* (yw1(ii)*yw1(jj) + + * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 + & +yw2(ii)*yw2(jj)*0.3330d0 ) + endif + jj=ii+3 + if(jj.le.4)then + sg3(ileft+ii-4) = sg3(ileft+ii-4) + + & wpt* (yw1(ii)*yw1(jj) + + * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 + & +yw2(ii)*yw2(jj)*0.3330d0 ) + endif + 10 continue + + else if(ileft.eq.3)then + do 20 ii=1,3 + jj=ii + sg0(ileft-3+ii) = sg0(ileft-3+ii) + + & wpt* (yw1(ii)*yw1(jj) + + * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 + & +yw2(ii)*yw2(jj)*0.3330d0 ) + jj=ii+1 + if(jj.le.3)then + sg1(ileft+ii-3) = sg1(ileft+ii-3) + + & wpt* (yw1(ii)*yw1(jj) + + * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 + & +yw2(ii)*yw2(jj)*0.3330d0 ) + endif + jj=ii+2 + if(jj.le.3)then + sg2(ileft+ii-3) = sg2(ileft+ii-3) + + & wpt* (yw1(ii)*yw1(jj) + + * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 + & +yw2(ii)*yw2(jj)*0.3330d0 ) + endif + 20 continue + + else if(ileft.eq.2)then + do 28 ii=1,2 + jj=ii + sg0(ileft-2+ii) = sg0(ileft-2+ii) + + & wpt* (yw1(ii)*yw1(jj) + + * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 + & +yw2(ii)*yw2(jj)*0.3330d0 ) + jj=ii+1 + if(jj.le.2)then + sg1(ileft+ii-2) = sg1(ileft+ii-2) + + & wpt* (yw1(ii)*yw1(jj) + + * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 + & +yw2(ii)*yw2(jj)*0.3330d0 ) + endif + 28 continue + + else if(ileft.eq.1)then + do 34 ii=1,1 + jj=ii + sg0(ileft-1+ii) = sg0(ileft-1+ii) + + & wpt* (yw1(ii)*yw1(jj) + + * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 + & +yw2(ii)*yw2(jj)*0.3330d0 ) + 34 continue + + endif + 2 continue + + return + end diff --git a/sinerp.f b/sinerp.f new file mode 100644 index 0000000..061aee1 --- /dev/null +++ b/sinerp.f @@ -0,0 +1,98 @@ +C Output from Public domain Ratfor, version 1.0 + subroutine sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,flag) +c +C Purpose : Computes Inner Products between columns of L^{-1} +C where L = abd is a Banded Matrix with 3 subdiagonals + +C The algorithm works in two passes: +C +C Pass 1 computes (cj,ck) k=j,j-1,j-2,j-3 ; j=nk, .. 1 +C Pass 2 computes (cj,ck) k <= j-4 (If flag == 1 ). +C +C A refinement of Elden's trick is used. +c Args + integer ld4,nk,ldnk,flag + DOUBLE precision abd(ld4,nk),p1ip(ld4,nk), p2ip(ldnk,nk) +c Locals + integer i,j,k + DOUBLE precision wjm3(3),wjm2(2),wjm1(1),c0,c1,c2,c3 +c +c unnecessary initialization of c1 c2 c3 to keep g77 -Wall happy +c + c1 = 0.0d0 + c2 = 0.0d0 + c3 = 0.0d0 +C +C Pass 1 + wjm3(1)=0d0 + wjm3(2)=0d0 + wjm3(3)=0d0 + wjm2(1)=0d0 + wjm2(2)=0d0 + wjm1(1)=0d0 + do 100 i=1,nk + j=nk-i+1 + c0 = 1d0/abd(4,j) + if(j.le.nk-3)then + c1 = abd(1,j+3)*c0 + c2 = abd(2,j+2)*c0 + c3 = abd(3,j+1)*c0 + else if(j.eq.nk-2)then + c1 = 0d0 + c2 = abd(2,j+2)*c0 + c3 = abd(3,j+1)*c0 + else if(j.eq.nk-1)then + c1 = 0d0 + c2 = 0d0 + c3 = abd(3,j+1)*c0 + else if(j.eq.nk)then + c1 = 0d0 + c2 = 0d0 + c3 = 0d0 + endif + p1ip(1,j) = 0d0- (c1*wjm3(1)+c2*wjm3(2)+c3*wjm3(3)) + p1ip(2,j) = 0d0- (c1*wjm3(2)+c2*wjm2(1)+c3*wjm2(2)) + p1ip(3,j) = 0d0- (c1*wjm3(3)+c2*wjm2(2)+c3*wjm1(1)) + p1ip(4,j) = c0**2 + c1**2*wjm3(1) + 2d0*c1*c2*wjm3(2)+ + & 2d0*c1*c3*wjm3(3) + c2**2*wjm2(1) + 2d0*c2*c3*wjm2(2) + + & c3**2*wjm1(1) + wjm3(1)=wjm2(1) + wjm3(2)=wjm2(2) + wjm3(3)=p1ip(2,j) + wjm2(1)=wjm1(1) + wjm2(2)=p1ip(3,j) + wjm1(1)=p1ip(4,j) + 100 continue + + if(flag.ne.0)then + +C ____ Pass 2 _____ + +C Compute p2ip + do 120 i=1,nk + j=nk-i+1 +C for(k=1;k<=4 & j+k-1<=nk;k=k+1) { p2ip(.) = .. }: + do 160 k=1,4 + if(j+k-1 .gt. nk)goto 120 + p2ip(j,j+k-1) = p1ip(5-k,j) + 160 continue + 120 continue + + do 170 i=1,nk + j=nk-i+1 +c for(k=j-4;k>=1;k=k-1){ + if(j-4 .ge. 1) then + do 210 k= j-4,1, -1 + c0 = 1d0/abd(4,k) + c1 = abd(1,k+3)*c0 + c2 = abd(2,k+2)*c0 + c3 = abd(3,k+1)*c0 + p2ip(k,j)= 0d0 - ( c1*p2ip(k+3,j) + c2*p2ip(k+2,j) + + & c3*p2ip(k+1,j) ) + 210 continue + endif + 170 continue + endif + return + end + diff --git a/spline.cpp b/spline.cpp new file mode 100644 index 0000000..f81bd4e --- /dev/null +++ b/spline.cpp @@ -0,0 +1,450 @@ + +#include "spline.h" + + +extern"C" { + void sgram_(double *sg0, double *sg1, double *sg2, + double *sg3, double* knot, int* nk); + void stxwx_(double *xs, double *ys, double *ws, int *n, + double *knot, int *nk, double *xwy, double *hs0, double *hs1, double *hs2, double *hs3); + void sslvrg_(double *penalt, double *dofoff, double *xs, double *ys, double *ws, double *ssw, int *n, + double *knot, int *nk, + double *coef, double *sz, double *lev, double *crit, int *icrit, double *lspar, double *xwy, + double *hs0, double *hs1, double *hs2, double *hs3, + double *sg0, double *sg1, double *sg2, double *sg3, double *abd, + double *p1ip, double *p2ip, int *ld4, int *ldnk, int *ier); +} +/***********************************************************/ +// used as reference - http://www.koders.com/fortran/fid8AA63B49CF22F0138E9B3DBDC405696F4A62C1CF.aspx +// http://www.koders.com/c/fidD995301A8A5549CE0361F4E7FFDFD3CDC4B4E4A3.aspx +/* A Cubic B-spline Smoothing routine. + + // sbart.f -- translated by f2c (version 20010821). + // ------- and f2c-clean,v 1.9 2000/01/13 + // + // According to the GAMFIT sources, this was derived from code by + // Finbarr O'Sullivan. + // + + + The algorithm minimises: + + (1/n) * sum ws(i)^2 * (ys(i)-sz(i))^2 + lambda* int ( s"(x) )^2 dx + + lambda is a function of the spar which is assumed to be between 0 and 1 + + INPUT + ----- + penalt A penalty > 1 to be used in the gcv criterion + dofoff either `df.offset' for GCV or `df' (to be matched). + n number of data points + ys(n) vector of length n containing the observations + ws(n) vector containing the weights given to each data point + NB: the code alters the values here. + xs(n) vector containing the ordinates of the observations + ssw `centered weighted sum of y^2' + nk number of b-spline coefficients to be estimated + nk <= n+2 + knot(nk+4) vector of knot points defining the cubic b-spline basis. + To obtain full cubic smoothing splines one might + have (provided the xs-values are strictly increasing) + spar penalised likelihood smoothing parameter + ispar indicating if spar is supplied (ispar=1) or to be estimated + lspar, uspar lower and upper values for spar search; 0.,1. are good values + tol, eps used in Golden Search routine + isetup setup indicator [initially 0 + icrit indicator saying which cross validation score is to be computed + 0: none ; 1: GCV ; 2: CV ; 3: 'df matching' + ld4 the leading dimension of abd (ie ld4=4) + ldnk the leading dimension of p2ip (not referenced) + + OUTPUT + ------ + coef(nk) vector of spline coefficients + sz(n) vector of smoothed z-values + lev(n) vector of leverages + crit either ordinary or generalized CV score + spar if ispar != 1 + lspar == lambda (a function of spar and the design) + iter number of iterations needed for spar search (if ispar != 1) + ier error indicator + ier = 0 ___ everything fine + ier = 1 ___ spar too small or too big + problem in cholesky decomposition + + Working arrays/matrix + xwy X'Wy + hs0,hs1,hs2,hs3 the diagonals of the X'WX matrix + sg0,sg1,sg2,sg3 the diagonals of the Gram matrix SIGMA + abd (ld4,nk) [ X'WX + lambda*SIGMA ] in diagonal form + p1ip(ld4,nk) inner products between columns of L inverse + p2ip(ldnk,nk) all inner products between columns of L inverse + where L'L = [X'WX + lambda*SIGMA] NOT REFERENCED + */ + +int Spline::sbart +(double *penalt, double *dofoff, + double *xs, double *ys, double *ws, double *ssw, + int *n, double *knot, int *nk, double *coef, + double *sz, double *lev, double *crit, int *icrit, + double *spar, int *ispar, int *iter, double *lspar, + double *uspar, double *tol, double *eps, int *isetup, + int *ld4, int *ldnk, int *ier) +{ + try{ + /* A Cubic B-spline Smoothing routine. + + The algorithm minimises: + + (1/n) * sum ws(i)^2 * (ys(i)-sz(i))^2 + lambda* int ( s(x) )^2 dx + + lambda is a function of the spar which is assumed to be between 0 and 1 + + INPUT + ----- + penalt A penalty > 1 to be used in the gcv criterion + dofoff either `df.offset' for GCV or `df' (to be matched). + n number of data points + ys(n) vector of length n containing the observations + ws(n) vector containing the weights given to each data point + NB: the code alters the values here. + xs(n) vector containing the ordinates of the observations + ssw `centered weighted sum of y^2 + nk number of b-spline coefficients to be estimated + nk <= n+2 + knot(nk+4) vector of knot points defining the cubic b-spline basis. + To obtain full cubic smoothing splines one might + have (provided the xs-values are strictly increasing) + spar penalised likelihood smoothing parameter + ispar indicating if spar is supplied (ispar=1) or to be estimated + lspar, uspar lower and upper values for spar search; 0.,1. are good values + tol, eps used in Golden Search routine + isetup setup indicator [initially 0 + icrit indicator saying which cross validation score is to be computed + 0: none ; 1: GCV ; 2: CV ; 3: 'df matching' + ld4 the leading dimension of abd (ie ld4=4) + ldnk the leading dimension of p2ip (not referenced) + + OUTPUT + ------ + coef(nk) vector of spline coefficients + sz(n) vector of smoothed z-values + lev(n) vector of leverages + crit either ordinary or generalized CV score + spar if ispar != 1 + lspar == lambda (a function of spar and the design) + iter number of iterations needed for spar search (if ispar != 1) + ier error indicator + ier = 0 ___ everything fine + ier = 1 ___ spar too small or too big + problem in cholesky decomposition + + Working arrays/matrix + xwy XWy + hs0,hs1,hs2,hs3 the diagonals of the XWX matrix + sg0,sg1,sg2,sg3 the diagonals of the Gram matrix SIGMA + abd (ld4,nk) [ XWX + lambda*SIGMA ] in diagonal form + p1ip(ld4,nk) inner products between columns of L inverse + p2ip(ldnk,nk) all inner products between columns of L inverse + where LL = [XWX + lambda*SIGMA] NOT REFERENCED + */ + +#define CRIT(FX) (*icrit == 3 ? FX - 3. : FX) + /* cancellation in (3 + eps) - 3, but still...informative */ + +#define BIG_f (1e100) + + /* c_Gold is the squared inverse of the golden ratio */ + static const double c_Gold = 0.381966011250105151795413165634; + /* == (3. - sqrt(5.)) / 2. */ + + /* Local variables */ + static double ratio;/* must be static (not needed in R) */ + + double a, b, d, e, p, q, r, u, v, w, x; + double ax, fu, fv, fw, fx, bx, xm; + double t1, t2, tol1, tol2; + double* xwy = new double[*nk]; + double* hs0 = new double[*nk]; + double* hs1 = new double[*nk]; + double* hs2 = new double[*nk]; + double* hs3 = new double[*nk]; + double* sg0 = new double[*nk]; + double* sg1 = new double[*nk]; + double* sg2 = new double[*nk]; + double* sg3 = new double[*nk]; + double* abd = new double[*nk*(*ld4)]; + double* p1ip = new double[*nk*(*ld4)]; + double* p2ip = new double[*nk]; + int i, maxit; + + /* unnecessary initializations to keep -Wall happy */ + d = 0.; fu = 0.; u = 0.; + ratio = 1.; + + /* Compute SIGMA, X' W X, X' W z, trace ratio, s0, s1. + + SIGMA -> sg0,sg1,sg2,sg3 + X' W X -> hs0,hs1,hs2,hs3 + X' W Z -> xwy + */ + + /* trevor fixed this 4/19/88 + * Note: sbart, i.e. stxwx() and sslvrg() {mostly, not always!}, use + * the square of the weights; the following rectifies that */ + for (i = 0; i < *n; ++i) + if (ws[i] > 0.) + ws[i] = sqrt(ws[i]); + + if (*isetup == 0) { + /* SIGMA[i,j] := Int B''(i,t) B''(j,t) dt {B(k,.) = k-th B-spline} */ + sgram_(sg0, sg1, sg2, sg3, knot, nk); + stxwx_(xs, ys, ws, n, + knot, nk, + xwy, + hs0, hs1, hs2, hs3); + /* Compute ratio := tr(X' W X) / tr(SIGMA) */ + t1 = t2 = 0.; + for (i = 3 - 1; i < (*nk - 3); ++i) { + t1 += hs0[i]; + t2 += sg0[i]; + } + ratio = t1 / t2; + *isetup = 1; + } + /* Compute estimate */ + + if (*ispar == 1) { /* Value of spar supplied */ + *lspar = ratio * pow(16.0, (*spar * 6.0 - 2.0)); + sslvrg_(penalt, dofoff, xs, ys, ws, ssw, n, + knot, nk, + coef, sz, lev, crit, icrit, lspar, xwy, + hs0, hs1, hs2, hs3, + sg0, sg1, sg2, sg3, abd, + p1ip, p2ip, ld4, ldnk, ier); + /* got through check 2 */ + return 0; + } + + /* ELSE ---- spar not supplied --> compute it ! --------------------------- + + Use Forsythe Malcom and Moler routine to MINIMIZE criterion + f denotes the value of the criterion + + an approximation x to the point where f attains a minimum on + the interval (ax,bx) is determined. + */ + ax = *lspar; + bx = *uspar; + + /* INPUT + + ax left endpoint of initial interval + bx right endpoint of initial interval + f function subprogram which evaluates f(x) for any x + in the interval (ax,bx) + tol desired length of the interval of uncertainty of the final + result ( >= 0 ) + + OUTPUT + + fmin abcissa approximating the point where f attains a minimum + */ + + /* + The method used is a combination of golden section search and + successive parabolic interpolation. convergence is never much slower + than that for a fibonacci search. if f has a continuous second + derivative which is positive at the minimum (which is not at ax or + bx), then convergence is superlinear, and usually of the order of + about 1.324.... + the function f is never evaluated at two points closer together + than eps*abs(fmin) + (tol/3), where eps is approximately the square + root of the relative machine precision. if f is a unimodal + function and the computed values of f are always unimodal when + separated by at least eps*abs(x) + (tol/3), then fmin approximates + the abcissa of the global minimum of f on the interval ax,bx with + an error less than 3*eps*abs(fmin) + tol. if f is not unimodal, + then fmin may approximate a local, but perhaps non-global, minimum to + the same accuracy. + this function subprogram is a slightly modified version of the + algol 60 procedure localmin given in richard brent, algorithms for + minimization without derivatives, prentice - hall, inc. (1973). + + Double a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w + Double fu,fv,fw,fx,x + */ + + /* eps is approximately the square root of the relative machine + precision. + + - eps = 1e0 + - 10 eps = eps/2e0 + - tol1 = 1e0 + eps + - if (tol1 > 1e0) go to 10 + - eps = sqrt(eps) + R Version <= 1.3.x had + eps = .000244 ( = sqrt(5.954 e-8) ) + -- now eps is passed as argument + */ + + /* initialization */ + + maxit = *iter; + *iter = 0; + a = ax; + b = bx; + v = a + c_Gold * (b - a); + w = v; + x = v; + e = 0.; + *spar = x; + *lspar = ratio * pow(16.0, (*spar * 6.0 - 2.0)); + sslvrg_(penalt, dofoff, xs, ys, ws, ssw, n, + knot, nk, + coef, sz, lev, crit, icrit, lspar, xwy, + hs0, hs1, hs2, hs3, + sg0, sg1, sg2, sg3, abd, + p1ip, p2ip, ld4, ldnk, ier); + fx = *crit; + fv = fx; + fw = fx; + + /* main loop + --------- */ + while(*ier == 0) { /* L20: */ + + if (m->control_pressed) { return 0; } + + xm = (a + b) * .5; + tol1 = *eps * fabs(x) + *tol / 3.; + tol2 = tol1 * 2.; + ++(*iter); + + + /* Check the (somewhat peculiar) stopping criterion: note that + the RHS is negative as long as the interval [a,b] is not small:*/ + if (fabs(x - xm) <= tol2 - (b - a) * .5 || *iter > maxit) + goto L_End; + + + /* is golden-section necessary */ + + if (fabs(e) <= tol1 || + /* if had Inf then go to golden-section */ + fx >= BIG_f || fv >= BIG_f || fw >= BIG_f) goto L_GoldenSect; + + /* Fit Parabola */ + + r = (x - w) * (fx - fv); + q = (x - v) * (fx - fw); + p = (x - v) * q - (x - w) * r; + q = (q - r) * 2.; + if (q > 0.) + p = -p; + q = fabs(q); + r = e; + e = d; + + /* is parabola acceptable? Otherwise do golden-section */ + + if (fabs(p) >= fabs(.5 * q * r) || + q == 0.) + /* above line added by BDR; + * [the abs(.) >= abs() = 0 should have branched..] + * in FTN: COMMON above ensures q is NOT a register variable */ + + goto L_GoldenSect; + + if (p <= q * (a - x) || + p >= q * (b - x)) goto L_GoldenSect; + + + + /* Parabolic Interpolation step */ + d = p / q; + u = x + d; + + /* f must not be evaluated too close to ax or bx */ + if ((u - a < tol2 || b - u < tol2)) { + d = abs(tol1) * sgn(xm - x); + } + + goto L50; + /*------*/ + + L_GoldenSect: /* a golden-section step */ + + if (x >= xm) e = a - x; + else/* x < xm*/ e = b - x; + d = c_Gold * e; + + + L50: + u = x + ((fabs(d) >= tol1) ? d : (abs(tol1)*sgn(d))); + /* tol1 check : f must not be evaluated too close to x */ + + *spar = u; + *lspar = ratio * pow(16.0, (*spar * 6.0 - 2.0)); + sslvrg_(penalt, dofoff, xs, ys, ws, ssw, n, + knot, nk, + coef, sz, lev, crit, icrit, lspar, xwy, + hs0, hs1, hs2, hs3, + sg0, sg1, sg2, sg3, abd, + p1ip, p2ip, ld4, ldnk, ier); + fu = *crit; + + if(isnan(fu)) { + fu = 2. * BIG_f; + } + + /* update a, b, v, w, and x */ + + if (fu <= fx) { + if (u >= x) a = x; else b = x; + + v = w; fv = fw; + w = x; fw = fx; + x = u; fx = fu; + } + else { + if (u < x) a = u; else b = u; + + if (fu <= fw || w == x) { /* L70: */ + v = w; fv = fw; + w = u; fw = fu; + } else if (fu <= fv || v == x || v == w) { /* L80: */ + v = u; fv = fu; + } + } + }/* end main loop -- goto L20; */ + + L_End: + + *spar = x; + *crit = fx; + + //free memory + delete [] xwy; + delete [] hs0; + delete [] hs1; + delete [] hs2; + delete [] hs3; + delete [] sg0; + delete [] sg1; + delete [] sg2; + delete [] sg3; + delete [] abd; + delete [] p1ip; + delete [] p2ip; + + return 0; + /* sbart */ + + }catch(exception& e) { + m->errorOut(e, "Spline", "sbart"); + exit(1); + } +} + +/***********************************************************/ diff --git a/spline.h b/spline.h new file mode 100644 index 0000000..94a0c64 --- /dev/null +++ b/spline.h @@ -0,0 +1,38 @@ +#ifndef SPLINE +#define SPLINE + + +/* + * spline.h + * Mothur + * + * Created by westcott on 12/6/11. + * Copyright 2011 Schloss Lab. All rights reserved. + * + */ + +//This class was translated into c++ using c and fortran code from www.koders.com as a reference. + +#include "mothurout.h" + +class Spline { + + public: + Spline() { m = MothurOut::getInstance(); } + ~Spline() {} + + int sbart(double *penalt, double *dofoff, + double *xs, double *ys, double *ws, double *ssw, + int *n, double *knot, int *nk, double *coef, + double *sz, double *lev, double *crit, int *icrit, + double *spar, int *ispar, int *iter, double *lspar, + double *uspar, double *tol, double *eps, int *isetup, + int *ld4, int *ldnk, int *ier); + + private: + MothurOut* m; +}; + +#endif + + diff --git a/sslvrg.f b/sslvrg.f new file mode 100644 index 0000000..221beaf --- /dev/null +++ b/sslvrg.f @@ -0,0 +1,136 @@ +C Output from Public domain Ratfor, version 1.0 + subroutine sslvrg(penalt,dofoff,x,y,w,ssw, n, knot,nk,coef, + * sz,lev, crit,icrit, lambda, xwy, hs0,hs1,hs2,hs3, + * sg0,sg1,sg2,sg3, abd,p1ip,p2ip,ld4,ldnk,info) + +C Purpose : +C Compute smoothing spline for smoothing parameter lambda +C and compute one of three `criteria' (OCV , GCV , "df match"). +C See comments in ./sbart.f from which this is called + + integer n,nk,icrit,ld4,ldnk,info + DOUBLE precision penalt,dofoff,x(n),y(n),w(n),ssw, + & knot(nk+4), coef(nk),sz(n),lev(n), crit, lambda, + * xwy(nk), hs0(nk),hs1(nk),hs2(nk),hs3(nk), + * sg0(nk),sg1(nk),sg2(nk),sg3(nk), abd(ld4,nk), + & p1ip(ld4,nk),p2ip(ldnk,nk) + + EXTERNAL bvalue + double precision bvalue +C local variables + double precision vnikx(4,1),work(16) + integer i,ileft,j,mflag, lenkno, ilo + double precision b0,b1,b2,b3,eps, xv,rss,df, sumw +c +c integer interv +c external interv + + lenkno = nk+4 + ileft = 1 + eps = 1d-11 + ilo = 1 + +C compute the coefficients coef() of estimated smooth + + do 1 i=1,nk + coef(i) = xwy(i) + abd(4,i) = hs0(i)+lambda*sg0(i) + 1 continue + + do 4 i=1,(nk-1) + abd(3,i+1) = hs1(i)+lambda*sg1(i) + 4 continue + + do 6 i=1,(nk-2) + 6 abd(2,i+2) = hs2(i)+lambda*sg2(i) + + do 8 i=1,(nk-3) + 8 abd(1,i+3) = hs3(i)+lambda*sg3(i) + +c factorize banded matrix abd: + call dpbfa(abd,ld4,nk,3,info) + if(info.ne.0) then +C matrix could not be factorized -> ier := info + return + endif +c solve linear system (from factorize abd): + call dpbsl(abd,ld4,nk,3,coef) + +C Value of smooth at the data points + do 12 i=1,n + xv = x(i) + 12 sz(i) = bvalue(knot,coef,nk,4,xv,0) + +C Compute the criterion function if requested + + if(icrit .eq. 0)then + return + else +C --- Ordinary or Generalized CV or "df match" --- + +C Get Leverages First + call sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,0) + do 16 i=1,n + xv = x(i) + call interv(knot(1), nk+1, xv, ileft, mflag) + if(mflag .eq. -1) then + ileft = 4 + xv = knot(4)+eps + else if(mflag .eq. 1) then + ileft = nk + xv = knot(nk+1) - eps + endif + j=ileft-3 +C call bspvd(knot,4,1,xv,ileft,4,vnikx,work) + call bsplvd(knot,lenkno,4,xv,ileft,work,vnikx,1) + b0=vnikx(1,1) + b1=vnikx(2,1) + b2=vnikx(3,1) + b3=vnikx(4,1) + lev(i) = ( + & p1ip(4,j)*b0**2 + 2.d0*p1ip(3,j)*b0*b1 + + * 2.d0*p1ip(2,j)*b0*b2 + 2.d0*p1ip(1,j)*b0*b3 + + * p1ip(4,j+1)*b1**2 + 2.d0*p1ip(3,j+1)*b1*b2 + + * 2.d0*p1ip(2,j+1)*b1*b3 + p1ip(4,j+2)*b2**2 + + & 2.d0*p1ip(3,j+2)*b2*b3 + p1ip(4,j+3)*b3**2 + & )*w(i)**2 + 16 continue + +C Evaluate Criterion + + if(icrit .eq. 1)then +C Generalized CV + rss = ssw + df = 0d0 + sumw = 0d0 +c w(i) are sqrt( wt[i] ) weights scaled in ../R/smspline.R such +c that sumw = number of observations with w(i) > 0 + do 24 i=1,n + rss = rss + ((y(i)-sz(i))*w(i))**2 + df = df + lev(i) + sumw = sumw + w(i)**2 + 24 continue + + crit = (rss/sumw)/((1d0-(dofoff + penalt*df)/sumw)**2) +c call dblepr("spar", 4, spar, 1) +c call dblepr("crit", 4, crit, 1) + + else if(icrit .eq. 2) then +C Ordinary CV + crit = 0d0 + do 30 i = 1,n + 30 crit = crit + (((y(i)-sz(i))*w(i))/(1-lev(i)))**2 + crit = crit/n +c call dblepr("spar", 4, spar, 1) +c call dblepr("crit", 4, crit, 1) + else +C df matching + crit = 0d0 + do 32 i=1,n + 32 crit = crit+lev(i) + crit = 3 + (dofoff-crit)**2 + endif + return + endif +C Criterion evaluation + end diff --git a/stxwx.f b/stxwx.f new file mode 100644 index 0000000..5322228 --- /dev/null +++ b/stxwx.f @@ -0,0 +1,65 @@ +C Output from Public domain Ratfor, version 1.0 + subroutine stxwx(x,z,w,k,xknot,n,y,hs0,hs1,hs2,hs3) + +c implicit none + integer k,n + DOUBLE precision x(k),z(k),w(k), xknot(n+4),y(n), + & hs0(n),hs1(n),hs2(n),hs3(n) +C local + DOUBLE precision eps,vnikx(4,1),work(16) + integer lenxk, i,j, ileft,mflag +c +c integer interv +c external interv + + lenxk=n+4 +C Initialise the output vectors + do 1 i=1,n + y(i)=0d0 + hs0(i)=0d0 + hs1(i)=0d0 + hs2(i)=0d0 + hs3(i)=0d0 + 1 continue + +C Compute X' W^2 X -> hs0,hs1,hs2,hs3 and X' W^2 Z -> y +C Note that here the weights w(i) == sqrt(wt[i]) where wt[] where original weights + ileft=1 + eps= .1d-9 + + do 100 i=1,k + call interv(xknot(1), n+1, x(i), ileft, mflag) +C if(mflag==-1) {write(6,'("Error in hess ",i2)')mflag;stop} +C if(mflag==-1) {return} + if(mflag.eq. 1)then + if(x(i).le.(xknot(ileft)+eps))then + ileft=ileft-1 + else + return + endif +C else{write(6,'("Error in hess ",i2)')mflag;stop}} + endif + call bsplvd (xknot,lenxk,4,x(i),ileft,work,vnikx,1) + + j= ileft-4+1 + y(j) = y(j)+w(i)**2*z(i)*vnikx(1,1) + hs0(j)=hs0(j)+w(i)**2*vnikx(1,1)**2 + hs1(j)=hs1(j)+w(i)**2*vnikx(1,1)*vnikx(2,1) + hs2(j)=hs2(j)+w(i)**2*vnikx(1,1)*vnikx(3,1) + hs3(j)=hs3(j)+w(i)**2*vnikx(1,1)*vnikx(4,1) + j= ileft-4+2 + y(j) = y(j)+w(i)**2*z(i)*vnikx(2,1) + hs0(j)=hs0(j)+w(i)**2*vnikx(2,1)**2 + hs1(j)=hs1(j)+w(i)**2*vnikx(2,1)*vnikx(3,1) + hs2(j)=hs2(j)+w(i)**2*vnikx(2,1)*vnikx(4,1) + j= ileft-4+3 + y(j) = y(j)+w(i)**2*z(i)*vnikx(3,1) + hs0(j)=hs0(j)+w(i)**2*vnikx(3,1)**2 + hs1(j)=hs1(j)+w(i)**2*vnikx(3,1)*vnikx(4,1) + j= ileft-4+4 + y(j) = y(j)+w(i)**2*z(i)*vnikx(4,1) + hs0(j)=hs0(j)+w(i)**2*vnikx(4,1)**2 + 100 continue + + return + end diff --git a/trimflowscommand.cpp b/trimflowscommand.cpp index 8be4a6c..c2ea534 100644 --- a/trimflowscommand.cpp +++ b/trimflowscommand.cpp @@ -202,7 +202,7 @@ TrimFlowsCommand::TrimFlowsCommand(string option) { } catch(exception& e) { - m->errorOut(e, "TrimFlowsCommand", "TrimSeqsCommand"); + m->errorOut(e, "TrimFlowsCommand", "TrimFlowsCommand"); exit(1); } } diff --git a/trimseqscommand.cpp b/trimseqscommand.cpp index 89ea3f9..b804c8f 100644 --- a/trimseqscommand.cpp +++ b/trimseqscommand.cpp @@ -1019,33 +1019,7 @@ int TrimSeqsCommand::setLines(string filename, string qfilename, vector