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 */; };
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 */; };
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;
A73DDBB913C4A0D1006AAE38 /* clearmemorycommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = clearmemorycommand.cpp; sourceTree = "<group>"; };
A73DDC3613C4BF64006AAE38 /* mothurmetastats.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = mothurmetastats.h; sourceTree = "<group>"; };
A73DDC3713C4BF64006AAE38 /* mothurmetastats.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = mothurmetastats.cpp; sourceTree = "<group>"; };
+ A74A9A9D148E881E00AB5E3E /* spline.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = spline.h; sourceTree = "<group>"; };
+ A74A9A9E148E881E00AB5E3E /* spline.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = spline.cpp; sourceTree = "<group>"; };
A74D3655137DAB8300332B0C /* addtargets2.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = addtargets2.cpp; sourceTree = "<group>"; };
A74D3656137DAB8300332B0C /* alignchime.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = alignchime.cpp; sourceTree = "<group>"; };
A74D3657137DAB8300332B0C /* alignchimel.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = alignchimel.cpp; sourceTree = "<group>"; };
A7F9F5CE141A5E500032F693 /* sequenceparser.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = sequenceparser.cpp; sourceTree = "<group>"; };
A7FA10001302E096003860FE /* mantelcommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = mantelcommand.h; sourceTree = "<group>"; };
A7FA10011302E096003860FE /* mantelcommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = mantelcommand.cpp; sourceTree = "<group>"; };
+ A7FA2ABC14A0E881007C09A6 /* bsplvb.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = bsplvb.f; sourceTree = "<group>"; };
+ A7FA2ABD14A0E881007C09A6 /* bsplvd.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = bsplvd.f; sourceTree = "<group>"; };
+ A7FA2ABE14A0E881007C09A6 /* bvalue.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = bvalue.f; sourceTree = "<group>"; };
+ A7FA2ABF14A0E881007C09A6 /* daxpy.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = daxpy.f; sourceTree = "<group>"; };
+ A7FA2AC014A0E881007C09A6 /* ddot.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = ddot.f; sourceTree = "<group>"; };
+ A7FA2AC114A0E881007C09A6 /* dpbfa.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = dpbfa.f; sourceTree = "<group>"; };
+ A7FA2AC214A0E881007C09A6 /* dpbsl.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = dpbsl.f; sourceTree = "<group>"; };
+ A7FA2AC314A0E881007C09A6 /* interv.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = interv.f; sourceTree = "<group>"; };
+ A7FA2AC414A0E881007C09A6 /* sgram.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = sgram.f; sourceTree = "<group>"; };
+ A7FA2AC514A0E881007C09A6 /* sinerp.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = sinerp.f; sourceTree = "<group>"; };
+ A7FA2AC614A0E881007C09A6 /* stxwx.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = stxwx.f; sourceTree = "<group>"; };
+ A7FA2B1514A0EBEA007C09A6 /* sslvrg.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = sslvrg.f; sourceTree = "<group>"; };
+ A7FA2B5A14A0F0C2007C09A6 /* intrv.f */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.fortran; path = intrv.f; sourceTree = "<group>"; };
A7FC480C12D788F20055BC5C /* linearalgebra.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = linearalgebra.h; sourceTree = "<group>"; };
A7FC480D12D788F20055BC5C /* linearalgebra.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = linearalgebra.cpp; sourceTree = "<group>"; };
A7FC486512D795D60055BC5C /* pcacommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = pcacommand.h; sourceTree = "<group>"; };
name = uchime;
sourceTree = "<group>";
};
+ 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 = "<group>";
+ };
A7E9BA3812D3956100DA6239 /* commands */ = {
isa = PBXGroup;
children = (
A7E9BA5612D39BD800DA6239 /* metastats */ = {
isa = PBXGroup;
children = (
+ A7D161E7149F7F50000523E8 /* fortran */,
A7E9B6E512D37EC400DA6239 /* fisher2.c */,
A7E9B6E612D37EC400DA6239 /* fisher2.h */,
A7E9B75512D37EC400DA6239 /* metastats.h */,
A79234D613C74BF6002B08E2 /* mothurfisher.cpp */,
A73DDC3613C4BF64006AAE38 /* mothurmetastats.h */,
A73DDC3713C4BF64006AAE38 /* mothurmetastats.cpp */,
+ A74A9A9D148E881E00AB5E3E /* spline.h */,
+ A74A9A9E148E881E00AB5E3E /* spline.cpp */,
);
name = metastats;
sourceTree = "<group>";
8DD76FAF0486AB0100D96B5E /* CopyFiles */,
);
buildRules = (
+ A7D162CB149F96CA000523E8 /* PBXBuildRule */,
);
dependencies = (
);
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;
};
"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;
);
PREBINDING = NO;
SDKROOT = macosx10.6;
+ USER_HEADER_SEARCH_PATHS = "";
};
name = Debug;
};
--- /dev/null
+ 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
--- /dev/null
+ 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
--- /dev/null
+ 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
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 + "\"";
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();
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;
--- /dev/null
+ 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
--- /dev/null
+*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
--- /dev/null
+ 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
+
--- /dev/null
+ 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
+
--- /dev/null
+ 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
--- /dev/null
+*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
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
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
throw BadConversion(s);
}
-
+//**********************************************************************************************************************
+template <typename T> int sgn(T val){ return (val > T(0)) - (val < T(0)); }
//**********************************************************************************************************************
template<typename T>
#include "mothurmetastats.h"
#include "mothurfisher.h"
+#include "spline.h"
/***********************************************************/
MothurMetastats::MothurMetastats(double t, int n) {
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<double> qvalues = calc_qvalues(pvalues);
// BACKUP checks
cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
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();
}
}
//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;
}
/***********************************************************/
vector<double> MothurMetastats::calc_qvalues(vector<double>& pValues) {
try {
- vector<double> qvalues;
- int numRows = pValues.size();
+ int numRows = pValues.size();
+ vector<double> qvalues(numRows, 0.0);
+
//fill lambdas with 0.00, 0.01, 0.02... 0.95
vector<double> lambdas(96, 0);
for (int i = 1; i < lambdas.size(); i++) { lambdas[i] = lambdas[i-1] + 0.01; }
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<double> 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<double> 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<double> ordered_qs = qvalues;
for (int i = 1; i < ordered_ps.size(); i++) { ordered_ps[i] = ordered_ps[i-1] + 1; }
vector<double> 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
exit(1);
}
}
-/***********************************************************
-vector<double> MothurMetastats::smoothSpline(vector<double> x, vector<double> y, int df) {
+/***********************************************************/
+double MothurMetastats::smoothSpline(vector<double>& x, vector<double>& 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<double> 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<double> wbar(n, 0);
- vector<double> ybar(n, 0);
- vector<double> 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<double> knot = sknot1(xbar);
+ vector<double> 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");
}
}
/***********************************************************/
-// This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
-int MothurMetastats::spline(vector<double>& x, vector<double>& y, int yp1, int ypn, vector<double>& 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<double> 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<double>& xa, vector<double>& ya, double& x, double& y, vector<double>& 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<double> MothurMetastats::sknot1(vector<double> x) {
+vector<double> MothurMetastats::sknot1(double* x, int n) {
try {
vector<double> 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<int> 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) {
int permute_matrix(vector<double>&, vector<double>&, int, vector<double>&, vector<double>&, vector<double>&);
int permute_array(vector<int>&);
int calc_twosample_ts(vector<double>&, int, vector<double>&, vector<double>&, vector<double>&);
- vector<double> smoothSpline(vector<double>, vector<double>, int);
+ double smoothSpline(vector<double>&, vector<double>&, int);
vector<double> calc_qvalues(vector<double>&);
- vector<double> sknot1(vector<double>);
+ vector<double> sknot1(double*, int);
int nkn(int);
int OrderPValues(int, int, vector<double>&, vector<int>&);
int swapElements(int, int, vector<double>&, vector<int>&);
vector<int> getSequence(int, int, int);
- int spline(vector<double>&, vector<double>&, int, int, vector<double>&);
- int splint(vector<double>&, vector<double>&, double&, double&, vector<double>&);
-
-
};
#endif
int level = 1;
tree[0].accessions.push_back(seqName);
+ m->removeConfidences(seqTaxonomy);
+
string taxon;// = getNextTaxon(seqTaxonomy);
while(seqTaxonomy != ""){
--- /dev/null
+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
--- /dev/null
+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
+
--- /dev/null
+
+#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);
+ }
+}
+
+/***********************************************************/
--- /dev/null
+#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
+
+
--- /dev/null
+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
--- /dev/null
+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
}
catch(exception& e) {
- m->errorOut(e, "TrimFlowsCommand", "TrimSeqsCommand");
+ m->errorOut(e, "TrimFlowsCommand", "TrimFlowsCommand");
exit(1);
}
}
#else
fastaFilePos.push_back(0); qfileFilePos.push_back(0);
- //get last file position of fastafile
- FILE * pFile;
- unsigned long long size;
-
- //get num bytes in file
- pFile = fopen (filename.c_str(),"rb");
- if (pFile==NULL) perror ("Error opening file");
- else{
- fseek (pFile, 0, SEEK_END);
- size=ftell (pFile);
- fclose (pFile);
- }
- fastaFilePos.push_back(size);
-
- //get last file position of fastafile
- FILE * qFile;
-
- //get num bytes in file
- qFile = fopen (qfilename.c_str(),"rb");
- if (qFile==NULL) perror ("Error opening file");
- else{
- fseek (qFile, 0, SEEK_END);
- size=ftell (qFile);
- fclose (qFile);
- }
- qfileFilePos.push_back(size);
-
+ fastaFilePos.push_back(1000); qfileFilePos.push_back(1000);
return 1;
#endif
projects / mothur.git / commitdiff