From: westcott Date: Fri, 23 Dec 2011 13:45:29 +0000 (+0000) Subject: added qvalues to metastats. fixed bug with chimera.uchime location. fixed windows... X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=commitdiff_plain;h=d0051dc9939d3477bd92b42c86bcd3eda743b955 added qvalues to metastats. fixed bug with chimera.uchime location. fixed windows bug with tellg() in trim.seqs. --- 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