Actualiser src/global.mk

be consistent across systems (an error would be generated on Ubuntu when
no primers found)

.cproject

@@ -1,151 +1,221 @@
-<?xml version="1.0" encoding="UTF-8"?>
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
 <?fileVersion 4.0.0?>
 
 <cproject>
-<storageModule moduleId="org.eclipse.cdt.core.settings">
-<cconfiguration id="cdt.managedbuild.toolchain.gnu.macosx.base.2134184396">
-<storageModule buildSystemId="org.eclipse.cdt.managedbuilder.core.configurationDataProvider" id="cdt.managedbuild.toolchain.gnu.macosx.base.2134184396" moduleId="org.eclipse.cdt.core.settings" name="MacOSX GCC">
-<externalSettings/>
-<extensions>
-<extension id="org.eclipse.cdt.core.MachO" point="org.eclipse.cdt.core.BinaryParser"/>
-<extension id="org.eclipse.cdt.core.GASErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
-<extension id="org.eclipse.cdt.core.GLDErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
-<extension id="org.eclipse.cdt.core.GCCErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
-<extension id="org.eclipse.cdt.core.MakeErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
-</extensions>
-</storageModule>
-<storageModule moduleId="cdtBuildSystem" version="4.0.0">
-<configuration artifactName="ecoPrimers" buildProperties="" id="cdt.managedbuild.toolchain.gnu.macosx.base.2134184396" name="MacOSX GCC" parent="org.eclipse.cdt.build.core.emptycfg">
-<folderInfo id="cdt.managedbuild.toolchain.gnu.macosx.base.2134184396.1840911077" name="/" resourcePath="">
-<toolChain id="cdt.managedbuild.toolchain.gnu.macosx.base.766054112" name="cdt.managedbuild.toolchain.gnu.macosx.base" superClass="cdt.managedbuild.toolchain.gnu.macosx.base">
-<targetPlatform archList="all" binaryParser="org.eclipse.cdt.core.MachO" id="cdt.managedbuild.target.gnu.platform.macosx.base.2057035265" name="Debug Platform" osList="macosx" superClass="cdt.managedbuild.target.gnu.platform.macosx.base"/>
-<builder id="cdt.managedbuild.target.gnu.builder.macosx.base.783726363" managedBuildOn="false" name="Gnu Make Builder.MacOSX GCC" superClass="cdt.managedbuild.target.gnu.builder.macosx.base"/>
-<tool id="cdt.managedbuild.tool.macosx.c.linker.macosx.base.914103467" name="MacOS X C Linker" superClass="cdt.managedbuild.tool.macosx.c.linker.macosx.base">
-<inputType id="cdt.managedbuild.tool.macosx.c.linker.input.62980206" superClass="cdt.managedbuild.tool.macosx.c.linker.input">
-<additionalInput kind="additionalinputdependency" paths="$(USER_OBJS)"/>
-<additionalInput kind="additionalinput" paths="$(LIBS)"/>
-</inputType>
-</tool>
-<tool id="cdt.managedbuild.tool.macosx.cpp.linker.macosx.base.691108439" name="MacOS X C++ Linker" superClass="cdt.managedbuild.tool.macosx.cpp.linker.macosx.base"/>
-<tool id="cdt.managedbuild.tool.gnu.assembler.macosx.base.695639877" name="GCC Assembler" superClass="cdt.managedbuild.tool.gnu.assembler.macosx.base">
-<inputType id="cdt.managedbuild.tool.gnu.assembler.input.1507665054" superClass="cdt.managedbuild.tool.gnu.assembler.input"/>
-</tool>
-<tool id="cdt.managedbuild.tool.gnu.archiver.macosx.base.1786370580" name="GCC Archiver" superClass="cdt.managedbuild.tool.gnu.archiver.macosx.base"/>
-<tool id="cdt.managedbuild.tool.gnu.cpp.compiler.macosx.base.454329831" name="GCC C++ Compiler" superClass="cdt.managedbuild.tool.gnu.cpp.compiler.macosx.base"/>
-<tool id="cdt.managedbuild.tool.gnu.c.compiler.macosx.base.1928774909" name="GCC C Compiler" superClass="cdt.managedbuild.tool.gnu.c.compiler.macosx.base">
-<inputType id="cdt.managedbuild.tool.gnu.c.compiler.input.330854350" superClass="cdt.managedbuild.tool.gnu.c.compiler.input"/>
-</tool>
-</toolChain>
-</folderInfo>
-</configuration>
-</storageModule>
-<storageModule moduleId="scannerConfiguration">
-<autodiscovery enabled="true" problemReportingEnabled="true" selectedProfileId="org.eclipse.cdt.make.core.GCCStandardMakePerProjectProfile"/>
-<profile id="org.eclipse.cdt.make.core.GCCStandardMakePerProjectProfile">
-<buildOutputProvider>
-<openAction enabled="true" filePath=""/>
-<parser enabled="true"/>
-</buildOutputProvider>
-<scannerInfoProvider id="specsFile">
-<runAction arguments="-E -P -v -dD ${plugin_state_location}/${specs_file}" command="gcc" useDefault="true"/>
-<parser enabled="true"/>
-</scannerInfoProvider>
-</profile>
-<profile id="org.eclipse.cdt.make.core.GCCStandardMakePerFileProfile">
-<buildOutputProvider>
-<openAction enabled="true" filePath=""/>
-<parser enabled="true"/>
-</buildOutputProvider>
-<scannerInfoProvider id="makefileGenerator">
-<runAction arguments="-f ${project_name}_scd.mk" command="make" useDefault="true"/>
-<parser enabled="true"/>
-</scannerInfoProvider>
-</profile>
-<profile id="org.eclipse.cdt.managedbuilder.core.GCCManagedMakePerProjectProfile">
-<buildOutputProvider>
-<openAction enabled="true" filePath=""/>
-<parser enabled="true"/>
-</buildOutputProvider>
-<scannerInfoProvider id="specsFile">
-<runAction arguments="-E -P -v -dD ${plugin_state_location}/${specs_file}" command="gcc" useDefault="true"/>
-<parser enabled="true"/>
-</scannerInfoProvider>
-</profile>
-<profile id="org.eclipse.cdt.managedbuilder.core.GCCManagedMakePerProjectProfileCPP">
-<buildOutputProvider>
-<openAction enabled="true" filePath=""/>
-<parser enabled="true"/>
-</buildOutputProvider>
-<scannerInfoProvider id="specsFile">
-<runAction arguments="-E -P -v -dD ${plugin_state_location}/specs.cpp" command="g++" useDefault="true"/>
-<parser enabled="true"/>
-</scannerInfoProvider>
-</profile>
-<profile id="org.eclipse.cdt.managedbuilder.core.GCCManagedMakePerProjectProfileC">
-<buildOutputProvider>
-<openAction enabled="true" filePath=""/>
-<parser enabled="true"/>
-</buildOutputProvider>
-<scannerInfoProvider id="specsFile">
-<runAction arguments="-E -P -v -dD ${plugin_state_location}/specs.c" command="gcc" useDefault="true"/>
-<parser enabled="true"/>
-</scannerInfoProvider>
-</profile>
-<profile id="org.eclipse.cdt.managedbuilder.core.GCCWinManagedMakePerProjectProfile">
-<buildOutputProvider>
-<openAction enabled="true" filePath=""/>
-<parser enabled="true"/>
-</buildOutputProvider>
-<scannerInfoProvider id="specsFile">
-<runAction arguments="-E -P -v -dD ${plugin_state_location}/${specs_file}" command="gcc" useDefault="true"/>
-<parser enabled="true"/>
-</scannerInfoProvider>
-</profile>
-<profile id="org.eclipse.cdt.managedbuilder.core.GCCWinManagedMakePerProjectProfileCPP">
-<buildOutputProvider>
-<openAction enabled="true" filePath=""/>
-<parser enabled="true"/>
-</buildOutputProvider>
-<scannerInfoProvider id="specsFile">
-<runAction arguments="-E -P -v -dD ${plugin_state_location}/specs.cpp" command="g++" useDefault="true"/>
-<parser enabled="true"/>
-</scannerInfoProvider>
-</profile>
-<profile id="org.eclipse.cdt.managedbuilder.core.GCCWinManagedMakePerProjectProfileC">
-<buildOutputProvider>
-<openAction enabled="true" filePath=""/>
-<parser enabled="true"/>
-</buildOutputProvider>
-<scannerInfoProvider id="specsFile">
-<runAction arguments="-E -P -v -dD ${plugin_state_location}/specs.c" command="gcc" useDefault="true"/>
-<parser enabled="true"/>
-</scannerInfoProvider>
-</profile>
-<profile id="org.eclipse.cdt.managedbuilder.xlc.core.XLCManagedMakePerProjectProfile">
-<buildOutputProvider>
-<openAction enabled="false" filePath=""/>
-<parser enabled="false"/>
-</buildOutputProvider>
-<scannerInfoProvider id="specsFile">
-<runAction arguments="-E -v ${plugin_state_location}/${specs_file}" command="${XL_compilerRoot}/xlc" useDefault="true"/>
-<parser enabled="true"/>
-</scannerInfoProvider>
-</profile>
-<profile id="org.eclipse.cdt.managedbuilder.xlc.core.XLCManagedMakePerProjectProfileCPP">
-<buildOutputProvider>
-<openAction enabled="false" filePath=""/>
-<parser enabled="false"/>
-</buildOutputProvider>
-<scannerInfoProvider id="specsFile">
-<runAction arguments="-E -v ${plugin_state_location}/${specs_file}" command="${XL_compilerRoot}/xlC" useDefault="true"/>
-<parser enabled="true"/>
-</scannerInfoProvider>
-</profile>
-</storageModule>
-<storageModule moduleId="org.eclipse.cdt.core.externalSettings"/>
-</cconfiguration>
-</storageModule>
-<storageModule moduleId="cdtBuildSystem" version="4.0.0">
-<project id="ecoPrimers.null.1292969001" name="ecoPrimers"/>
-</storageModule>
+	<storageModule moduleId="org.eclipse.cdt.core.settings">
+		<cconfiguration id="cdt.managedbuild.toolchain.gnu.macosx.base.2134184396">
+			<storageModule buildSystemId="org.eclipse.cdt.managedbuilder.core.configurationDataProvider" id="cdt.managedbuild.toolchain.gnu.macosx.base.2134184396" moduleId="org.eclipse.cdt.core.settings" name="MacOSX GCC">
+				<externalSettings/>
+				<extensions>
+					<extension id="org.eclipse.cdt.core.MachO" point="org.eclipse.cdt.core.BinaryParser"/>
+					<extension id="org.eclipse.cdt.core.GASErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
+					<extension id="org.eclipse.cdt.core.GLDErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
+					<extension id="org.eclipse.cdt.core.GCCErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
+					<extension id="org.eclipse.cdt.core.GmakeErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
+					<extension id="org.eclipse.cdt.core.CWDLocator" point="org.eclipse.cdt.core.ErrorParser"/>
+				</extensions>
+			</storageModule>
+			<storageModule moduleId="cdtBuildSystem" version="4.0.0">
+				<configuration artifactName="ecoPrimers" buildProperties="" description="" id="cdt.managedbuild.toolchain.gnu.macosx.base.2134184396" name="MacOSX GCC" parent="org.eclipse.cdt.build.core.emptycfg">
+					<folderInfo id="cdt.managedbuild.toolchain.gnu.macosx.base.2134184396.1840911077" name="/" resourcePath="">
+						<toolChain id="cdt.managedbuild.toolchain.gnu.macosx.base.766054112" name="cdt.managedbuild.toolchain.gnu.macosx.base" superClass="cdt.managedbuild.toolchain.gnu.macosx.base">
+							<targetPlatform archList="all" binaryParser="org.eclipse.cdt.core.MachO" id="cdt.managedbuild.target.gnu.platform.macosx.base.2057035265" name="Debug Platform" osList="macosx" superClass="cdt.managedbuild.target.gnu.platform.macosx.base"/>
+							<builder id="cdt.managedbuild.target.gnu.builder.macosx.base.783726363" keepEnvironmentInBuildfile="false" managedBuildOn="false" name="Gnu Make Builder" superClass="cdt.managedbuild.target.gnu.builder.macosx.base"/>
+							<tool id="cdt.managedbuild.tool.macosx.c.linker.macosx.base.914103467" name="MacOS X C Linker" superClass="cdt.managedbuild.tool.macosx.c.linker.macosx.base">
+								<inputType id="cdt.managedbuild.tool.macosx.c.linker.input.62980206" superClass="cdt.managedbuild.tool.macosx.c.linker.input">
+									<additionalInput kind="additionalinputdependency" paths="$(USER_OBJS)"/>
+									<additionalInput kind="additionalinput" paths="$(LIBS)"/>
+								</inputType>
+							</tool>
+							<tool id="cdt.managedbuild.tool.macosx.cpp.linker.macosx.base.691108439" name="MacOS X C++ Linker" superClass="cdt.managedbuild.tool.macosx.cpp.linker.macosx.base"/>
+							<tool id="cdt.managedbuild.tool.gnu.assembler.macosx.base.695639877" name="GCC Assembler" superClass="cdt.managedbuild.tool.gnu.assembler.macosx.base">
+								<option id="gnu.both.asm.option.include.paths.1544375094" name="Include paths (-I)" superClass="gnu.both.asm.option.include.paths" valueType="includePath"/>
+								<inputType id="cdt.managedbuild.tool.gnu.assembler.input.1507665054" superClass="cdt.managedbuild.tool.gnu.assembler.input"/>
+							</tool>
+							<tool id="cdt.managedbuild.tool.gnu.archiver.macosx.base.1786370580" name="GCC Archiver" superClass="cdt.managedbuild.tool.gnu.archiver.macosx.base"/>
+							<tool id="cdt.managedbuild.tool.gnu.cpp.compiler.macosx.base.454329831" name="GCC C++ Compiler" superClass="cdt.managedbuild.tool.gnu.cpp.compiler.macosx.base"/>
+							<tool id="cdt.managedbuild.tool.gnu.c.compiler.macosx.base.1928774909" name="GCC C Compiler" superClass="cdt.managedbuild.tool.gnu.c.compiler.macosx.base">
+								<option id="gnu.c.compiler.option.include.paths.823251305" superClass="gnu.c.compiler.option.include.paths" valueType="includePath">
+									<listOptionValue builtIn="false" value="/usr/include"/>
+								</option>
+								<inputType id="cdt.managedbuild.tool.gnu.c.compiler.input.330854350" superClass="cdt.managedbuild.tool.gnu.c.compiler.input"/>
+							</tool>
+						</toolChain>
+					</folderInfo>
+				</configuration>
+			</storageModule>
+			<storageModule moduleId="org.eclipse.cdt.core.externalSettings"/>
+			<storageModule moduleId="org.eclipse.cdt.core.language.mapping"/>
+			<storageModule moduleId="scannerConfiguration">
+				<autodiscovery enabled="true" problemReportingEnabled="true" selectedProfileId="org.eclipse.cdt.make.core.GCCStandardMakePerProjectProfile"/>
+				<profile id="org.eclipse.cdt.make.core.GCCStandardMakePerProjectProfile">
+					<buildOutputProvider>
+						<openAction enabled="true" filePath=""/>
+						<parser enabled="true"/>
+					</buildOutputProvider>
+					<scannerInfoProvider id="specsFile">
+						<runAction arguments="-E -P -v -dD ${plugin_state_location}/${specs_file}" command="gcc" useDefault="true"/>
+						<parser enabled="true"/>
+					</scannerInfoProvider>
+				</profile>
+				<profile id="org.eclipse.cdt.make.core.GCCStandardMakePerFileProfile">
+					<buildOutputProvider>
+						<openAction enabled="true" filePath=""/>
+						<parser enabled="true"/>
+					</buildOutputProvider>
+					<scannerInfoProvider id="makefileGenerator">
+						<runAction arguments="-f ${project_name}_scd.mk" command="make" useDefault="true"/>
+						<parser enabled="true"/>
+					</scannerInfoProvider>
+				</profile>
+				<profile id="org.eclipse.cdt.managedbuilder.core.GCCManagedMakePerProjectProfile">
+					<buildOutputProvider>
+						<openAction enabled="true" filePath=""/>
+						<parser enabled="true"/>
+					</buildOutputProvider>
+					<scannerInfoProvider id="specsFile">
+						<runAction arguments="-E -P -v -dD ${plugin_state_location}/${specs_file}" command="gcc" useDefault="true"/>
+						<parser enabled="true"/>
+					</scannerInfoProvider>
+				</profile>
+				<profile id="org.eclipse.cdt.managedbuilder.core.GCCManagedMakePerProjectProfileCPP">
+					<buildOutputProvider>
+						<openAction enabled="true" filePath=""/>
+						<parser enabled="true"/>
+					</buildOutputProvider>
+					<scannerInfoProvider id="specsFile">
+						<runAction arguments="-E -P -v -dD ${plugin_state_location}/specs.cpp" command="g++" useDefault="true"/>
+						<parser enabled="true"/>
+					</scannerInfoProvider>
+				</profile>
+				<profile id="org.eclipse.cdt.managedbuilder.core.GCCManagedMakePerProjectProfileC">
+					<buildOutputProvider>
+						<openAction enabled="true" filePath=""/>
+						<parser enabled="true"/>
+					</buildOutputProvider>
+					<scannerInfoProvider id="specsFile">
+						<runAction arguments="-E -P -v -dD ${plugin_state_location}/specs.c" command="gcc" useDefault="true"/>
+						<parser enabled="true"/>
+					</scannerInfoProvider>
+				</profile>
+				<profile id="org.eclipse.cdt.managedbuilder.core.GCCWinManagedMakePerProjectProfile">
+					<buildOutputProvider>
+						<openAction enabled="true" filePath=""/>
+						<parser enabled="true"/>
+					</buildOutputProvider>
+					<scannerInfoProvider id="specsFile">
+						<runAction arguments="-E -P -v -dD ${plugin_state_location}/${specs_file}" command="gcc" useDefault="true"/>
+						<parser enabled="true"/>
+					</scannerInfoProvider>
+				</profile>
+				<profile id="org.eclipse.cdt.managedbuilder.core.GCCWinManagedMakePerProjectProfileCPP">
+					<buildOutputProvider>
+						<openAction enabled="true" filePath=""/>
+						<parser enabled="true"/>
+					</buildOutputProvider>
+					<scannerInfoProvider id="specsFile">
+						<runAction arguments="-E -P -v -dD ${plugin_state_location}/specs.cpp" command="g++" useDefault="true"/>
+						<parser enabled="true"/>
+					</scannerInfoProvider>
+				</profile>
+				<profile id="org.eclipse.cdt.managedbuilder.core.GCCWinManagedMakePerProjectProfileC">
+					<buildOutputProvider>
+						<openAction enabled="true" filePath=""/>
+						<parser enabled="true"/>
+					</buildOutputProvider>
+					<scannerInfoProvider id="specsFile">
+						<runAction arguments="-E -P -v -dD ${plugin_state_location}/specs.c" command="gcc" useDefault="true"/>
+						<parser enabled="true"/>
+					</scannerInfoProvider>
+				</profile>
+				<scannerConfigBuildInfo instanceId="cdt.managedbuild.toolchain.gnu.macosx.base.2134184396;cdt.managedbuild.toolchain.gnu.macosx.base.2134184396.1840911077;cdt.managedbuild.tool.gnu.c.compiler.macosx.base.1928774909;cdt.managedbuild.tool.gnu.c.compiler.input.330854350">
+					<autodiscovery enabled="true" problemReportingEnabled="true" selectedProfileId="org.eclipse.cdt.managedbuilder.core.GCCManagedMakePerProjectProfileC"/>
+					<profile id="org.eclipse.cdt.make.core.GCCStandardMakePerProjectProfile">
+						<buildOutputProvider>
+							<openAction enabled="true" filePath=""/>
+							<parser enabled="true"/>
+						</buildOutputProvider>
+						<scannerInfoProvider id="specsFile">
+							<runAction arguments="-E -P -v -dD ${plugin_state_location}/${specs_file}" command="gcc" useDefault="true"/>
+							<parser enabled="true"/>
+						</scannerInfoProvider>
+					</profile>
+					<profile id="org.eclipse.cdt.make.core.GCCStandardMakePerFileProfile">
+						<buildOutputProvider>
+							<openAction enabled="true" filePath=""/>
+							<parser enabled="true"/>
+						</buildOutputProvider>
+						<scannerInfoProvider id="makefileGenerator">
+							<runAction arguments="-f ${project_name}_scd.mk" command="make" useDefault="true"/>
+							<parser enabled="true"/>
+						</scannerInfoProvider>
+					</profile>
+					<profile id="org.eclipse.cdt.managedbuilder.core.GCCManagedMakePerProjectProfile">
+						<buildOutputProvider>
+							<openAction enabled="true" filePath=""/>
+							<parser enabled="true"/>
+						</buildOutputProvider>
+						<scannerInfoProvider id="specsFile">
+							<runAction arguments="-E -P -v -dD ${plugin_state_location}/${specs_file}" command="gcc" useDefault="true"/>
+							<parser enabled="true"/>
+						</scannerInfoProvider>
+					</profile>
+					<profile id="org.eclipse.cdt.managedbuilder.core.GCCManagedMakePerProjectProfileCPP">
+						<buildOutputProvider>
+							<openAction enabled="true" filePath=""/>
+							<parser enabled="true"/>
+						</buildOutputProvider>
+						<scannerInfoProvider id="specsFile">
+							<runAction arguments="-E -P -v -dD ${plugin_state_location}/specs.cpp" command="g++" useDefault="true"/>
+							<parser enabled="true"/>
+						</scannerInfoProvider>
+					</profile>
+					<profile id="org.eclipse.cdt.managedbuilder.core.GCCManagedMakePerProjectProfileC">
+						<buildOutputProvider>
+							<openAction enabled="true" filePath=""/>
+							<parser enabled="true"/>
+						</buildOutputProvider>
+						<scannerInfoProvider id="specsFile">
+							<runAction arguments="-E -P -v -dD ${plugin_state_location}/specs.c" command="gcc" useDefault="true"/>
+							<parser enabled="true"/>
+						</scannerInfoProvider>
+					</profile>
+					<profile id="org.eclipse.cdt.managedbuilder.core.GCCWinManagedMakePerProjectProfile">
+						<buildOutputProvider>
+							<openAction enabled="true" filePath=""/>
+							<parser enabled="true"/>
+						</buildOutputProvider>
+						<scannerInfoProvider id="specsFile">
+							<runAction arguments="-E -P -v -dD ${plugin_state_location}/${specs_file}" command="gcc" useDefault="true"/>
+							<parser enabled="true"/>
+						</scannerInfoProvider>
+					</profile>
+					<profile id="org.eclipse.cdt.managedbuilder.core.GCCWinManagedMakePerProjectProfileCPP">
+						<buildOutputProvider>
+							<openAction enabled="true" filePath=""/>
+							<parser enabled="true"/>
+						</buildOutputProvider>
+						<scannerInfoProvider id="specsFile">
+							<runAction arguments="-E -P -v -dD ${plugin_state_location}/specs.cpp" command="g++" useDefault="true"/>
+							<parser enabled="true"/>
+						</scannerInfoProvider>
+					</profile>
+					<profile id="org.eclipse.cdt.managedbuilder.core.GCCWinManagedMakePerProjectProfileC">
+						<buildOutputProvider>
+							<openAction enabled="true" filePath=""/>
+							<parser enabled="true"/>
+						</buildOutputProvider>
+						<scannerInfoProvider id="specsFile">
+							<runAction arguments="-E -P -v -dD ${plugin_state_location}/specs.c" command="gcc" useDefault="true"/>
+							<parser enabled="true"/>
+						</scannerInfoProvider>
+					</profile>
+				</scannerConfigBuildInfo>
+			</storageModule>
+			<storageModule moduleId="org.eclipse.cdt.internal.ui.text.commentOwnerProjectMappings"/>
+		</cconfiguration>
+	</storageModule>
+	<storageModule moduleId="cdtBuildSystem" version="4.0.0">
+		<project id="ecoPrimers.null.1292969001" name="ecoPrimers"/>
+	</storageModule>
 </cproject>

.gitignore

@@ -0,0 +1,15 @@
+*.o
+/src/*.P
+/src/*.a
+/src/ecoPrimer
+/src/libecoPCR/*.P
+/src/libecoPCR/*.a
+/src/libecoPCR/*.o
+/src/libecoprimer/*.P
+/src/libecoprimer/*.a
+/src/libecoprimer/*.o
+/src/libthermo/*.P
+/src/libthermo/*.a
+/src/libthermo/*.o
+phylonorway.*
+src/ecoPrimers

.project

@@ -5,45 +5,26 @@
 	<projects>
 	</projects>
 	<buildSpec>
+		<buildCommand>
+			<name>org.python.pydev.PyDevBuilder</name>
+			<arguments>
+			</arguments>
+		</buildCommand>
 		<buildCommand>
 			<name>org.eclipse.cdt.managedbuilder.core.genmakebuilder</name>
 			<triggers>clean,full,incremental,</triggers>
 			<arguments>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.fullBuildTarget</key>
-					<value>all</value>
-				</dictionary>
 				<dictionary>
 					<key>?name?</key>
 					<value></value>
 				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.enableAutoBuild</key>
-					<value>false</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.enableFullBuild</key>
-					<value>true</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.enableCleanBuild</key>
-					<value>true</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.cleanBuildTarget</key>
-					<value>clean</value>
-				</dictionary>
 				<dictionary>
 					<key>org.eclipse.cdt.make.core.append_environment</key>
 					<value>true</value>
 				</dictionary>
 				<dictionary>
-					<key>org.eclipse.cdt.make.core.contents</key>
-					<value>org.eclipse.cdt.make.core.activeConfigSettings</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.useDefaultBuildCmd</key>
-					<value>true</value>
+					<key>org.eclipse.cdt.make.core.autoBuildTarget</key>
+					<value>all</value>
 				</dictionary>
 				<dictionary>
 					<key>org.eclipse.cdt.make.core.buildArguments</key>
@@ -54,13 +35,37 @@
 					<value>make</value>
 				</dictionary>
 				<dictionary>
-					<key>org.eclipse.cdt.make.core.autoBuildTarget</key>
+					<key>org.eclipse.cdt.make.core.cleanBuildTarget</key>
+					<value>clean</value>
+				</dictionary>
+				<dictionary>
+					<key>org.eclipse.cdt.make.core.contents</key>
+					<value>org.eclipse.cdt.make.core.activeConfigSettings</value>
+				</dictionary>
+				<dictionary>
+					<key>org.eclipse.cdt.make.core.enableAutoBuild</key>
+					<value>false</value>
+				</dictionary>
+				<dictionary>
+					<key>org.eclipse.cdt.make.core.enableCleanBuild</key>
+					<value>true</value>
+				</dictionary>
+				<dictionary>
+					<key>org.eclipse.cdt.make.core.enableFullBuild</key>
+					<value>true</value>
+				</dictionary>
+				<dictionary>
+					<key>org.eclipse.cdt.make.core.fullBuildTarget</key>
 					<value>all</value>
 				</dictionary>
 				<dictionary>
 					<key>org.eclipse.cdt.make.core.stopOnError</key>
 					<value>true</value>
 				</dictionary>
+				<dictionary>
+					<key>org.eclipse.cdt.make.core.useDefaultBuildCmd</key>
+					<value>true</value>
+				</dictionary>
 			</arguments>
 		</buildCommand>
 		<buildCommand>
@@ -73,5 +78,6 @@
 		<nature>org.eclipse.cdt.core.cnature</nature>
 		<nature>org.eclipse.cdt.managedbuilder.core.ScannerConfigNature</nature>
 		<nature>org.eclipse.cdt.managedbuilder.core.managedBuildNature</nature>
+		<nature>org.python.pydev.pythonNature</nature>
 	</natures>
 </projectDescription>

.pydevproject

@@ -0,0 +1,7 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<?eclipse-pydev version="1.0"?>
+
+<pydev_project>
+<pydev_property name="org.python.pydev.PYTHON_PROJECT_INTERPRETER">Python 2.6</pydev_property>
+<pydev_property name="org.python.pydev.PYTHON_PROJECT_VERSION">python 2.6</pydev_property>
+</pydev_project>

VERSION

@@ -0,0 +1 @@
+0.5

src/Makefile

@@ -1,4 +1,4 @@
-EXEC=ecoPrimer
+EXEC=ecoPrimers
 
 PRIMER_SRC= ecoprimer.c
 PRIMER_OBJ= $(patsubst %.c,%.o,$(PRIMER_SRC))
@@ -6,11 +6,12 @@ PRIMER_OBJ= $(patsubst %.c,%.o,$(PRIMER_SRC))
 
 SRCS= $(PRIMER_SRC) 
 
-LIB= -lecoPCR -lapat -lKMRK -lz -lm
+LIB= -lecoprimer -lecoPCR -lthermo -lz -lm 
 
-LIBFILE= libapat/libapat.a \
-         libecoPCR/libecoPCR.a \
-         libKMRK/libKMRK.a
+LIBFILE= libecoPCR/libecoPCR.a \
+         libecoprimer/libecoprimer.a \
+         libthermo/libthermo.a \
+	 
 
 
 include global.mk
@@ -26,8 +27,8 @@ all: $(EXEC)
 
 # executable compilation and link
 
-ecoPrimer: $(PRIMER_OBJ) $(LIBFILE) 
-	$(CC) $(LDFLAGS) -o $@ $< $(LIBPATH) $(LIB)
+ecoPrimers: $(PRIMER_OBJ) $(LIBFILE) 
+	$(CC) $(LDFLAGS) -O5 -o $@ $< $(LIBPATH) $(LIB)
 	
 
 ########
@@ -36,15 +37,14 @@ ecoPrimer: $(PRIMER_OBJ) $(LIBFILE)
 #
 ########
 
-libapat/libapat.a:
-	$(MAKE) -C libapat
-	
 libecoPCR/libecoPCR.a:
 	$(MAKE) -C libecoPCR
 	
-libKMRK/libKMRK.a:
-	$(MAKE) -C libKMRK
+libecoprimer/libecoprimer.a:
+	$(MAKE) -C libecoprimer
 	
+libthermo/libthermo.a:
+	$(MAKE) -C libthermo
 	
 ########
 #
@@ -54,10 +54,26 @@ libKMRK/libKMRK.a:
 
 clean:
 	rm -f *.o	
+	rm -f *.P
 	rm -f $(EXEC)
-	$(MAKE) -C libapat clean
 	$(MAKE) -C libecoPCR clean
-	$(MAKE) -C libKMRK clean
+	$(MAKE) -C libecoprimer clean
+	$(MAKE) -C libthermo clean
 		
 
-		
+		
+########
+#
+# clean for k2 to remove .o and .P files
+#
+########
+
+k2clean:
+	rm -f *.o	
+	rm -f *.P
+	rm -f libecoPCR/*.o
+	rm -f libecoPCR/*.P
+	rm -f libecoprimer/*.o
+	rm -f libecoprimer/*.P
+	rm -f libthermo/*.o
+	rm -f libthermo/*.P

src/global.mk

@@ -1,9 +1,12 @@
-MACHINE=MAC_OS_X
-LIBPATH= -Llibapat -LlibecoPCR -LlibKMRK
+MACHINE=LINUX
+LIBPATH= -LlibecoPCR -Llibecoprimer -Llibthermo -L/usr/local/lib
 MAKEDEPEND = gcc -D$(MACHINE) -M $(CPPFLAGS) -o $*.d $<
 
 CC=gcc
-CFLAGS= -W -Wall -O2 -g
+CFLAGS= -W -Wall -O5
+#CFLAGS= -W -Wall -O5 -m64 -g
+#CFLAGS= -W -Wall -O0 -m64  -g
+#CFLAGS= -W -Wall -O5 -fast -g
 
 default: all
 
@@ -15,4 +18,4 @@ default: all
 	@sed 's/\($*\)\.o[ :]*/\1.o $@ : /g' < $*.d > $@; \
 	rm -f $*.d; [ -s $@ ] || rm -f $@
 
-include $(SRCS:.c=.P)
+include $(SRCS:.c=.P)

src/libKMRK/KMRK.h

@@ -1,309 +0,0 @@
-#ifndef _KMRK_H_
-#define _KMRK_H_
-
-/********************************************
- ********************************************
- **
- ** Declaration of struct
- **
- ********************************************
- ********************************************/
-
-#include <stdint.h>
-#include "repseek_types.h"
-#include "KMRK_mask.h"
-
-/** 
- * Structure used to manipulate simultanously 
- * the v and n vector 
- * 
- */ 
-
-typedef struct {
-
-	int32_t  size;           /**< size of vectors */
-	int32_t  seqCount;       /**< count of concatenated sequences */
-	int32_t  complement;     /**< if seqCount > 1 and complement !=0
-                            *   then second sequence is the inverted complement
-                            *   strand of first one */
-	int32_t  symbmax;
-	int32_t* v;               /**< sequence vector */
-	int32_t* n;	              /**< linked list vector */
-	int32_t  limits[1];       /**< array of end limits of concatenated
-	                           *   sequences in v (array size is seqCount) */
-} vn_type;
-
-
-
-/********************************************
- ********************************************
- **
- ** Declaration of public macro
- **
- ********************************************
- ********************************************/
-
-// Return a pointer to a vector from a vn_type structure
-#define GETVECTOR(x,vector)  (((x)->vector) - 1)
-
-#define IS_MARKED(x,i)  ((x)[i] < 0)
-#define MARK(x,i)       ((x)[i]) = -ABS((x)[i])
-#define UNMARK(x,i)     ((x)[i]) = ABS((x)[i])
-
-#define SET(x,i,v)      ((x)[i]) = (v)
-
-// set and mark in one operation
-#define SETMARKED(x,i,v) ((x)[i]) = -(v)
-
-//internal macro
-#define GET(x,i)        ABS((x)[i])
-
-// get symbole at position i in vector x
-#define SYMBOLE(x,i)    ((IS_MARKED((x),(i))) ? (i): (GET(x,i)))
-
-
-
-/** 
- * Macro used to declare a pointer to a quorum function.
- * 
- * @param name name of the pointer
- * 
- */
-
-#define KMRK_QUORUM_FUNC_PTR(name) int32_t (*name)(vn_type* x,                  \
-                                                   int32_t pos,                 \
-                                                   int32_t count,               \
-                                                   int32_t countSeq)
-
-
-
-
-/** 
- * Macro used to declare a pointer to an initialisation function.
- * 
- * @param name name of the pointer
- * @param quorum name used for the quorum assiciated function
- *
- * @see KMRK_QUORUM_FUNC_PTR
- * 
- */
-
-#define KMRK_INIT_FUNC_PTR(name,quorum)  vn_type* (*name)(char *sequence,              \
-                                                          int32_t complement,          \
-                                                          int32_t count,               \
-                                                          int32_t countSeq,            \
-                                                          int32_t *k,                  \
-                                                          KMRK_QUORUM_FUNC_PTR(quorum),\
-                                                          masked_area_table_t *mask)
-
-
-/********************************************
- ********************************************
- **
- ** Declaration of public functions
- **
- ********************************************
- ********************************************/
-
-
-
-
-
-/** 
- * Initialise a vn_type record from one sequence to run KMRK algorithm
- * 
- * @param sequence pointer to a C string containing the sequence
- * @param complement != 0 means that seq one and two are the two strands of 
- *                   the same sequence.
- * @param count parameter count passed to the quorun function
- * @param countSeq parametter countSeq passed to the quorun function
- * @param k length of the word represented by each symbole of v. 
- *        k is an output parametter
- * @param quorum pointer to a quorum function
- * 
- * @return a pointer to vn_type structure correctly initialized 
- *          to be used by KMRK_RunKMRK
- *
- * @see KMRK_HashOneSequence
- */
-
-vn_type* KMRK_InitOneSequence(char *sequence,              
-                              int32_t complement, 
-                              int32_t count,
-                              int32_t countSeq,
-                              int32_t *k,
-                              KMRK_QUORUM_FUNC_PTR(quorum),
-                              masked_area_table_t *mask);
-
-
-
-
-
-/** 
- * Initialise a vn_type record from one sequence to run KMRK algorithme.
- * In more than KMRK_InitOneSequence, KMRK_HashOneSequence construct 
- * word of len k with an hash algorithm. k used is a function of
- * sequence size and alphabet size. If calculed k is superior to lmax
- * then k = lmax. 
- * 
- * @param sequence   pointer to a C string containing the sequence
- * @param complement != 0 means that seq one and two are the two strands of 
- *                   the same sequence.
- * @param count      parametter count passed to the quorun function
- * @param countSeq   parametter countSeq passed to the quorun function
- * @param k          maximum length of the created word (input)
- *                   length of the word represented by each symbole 
- *                   of v (output). 
- * @param quorum     pointer to a quorum function
- * 
- * @return a pointer to vn_type structure correctly initialized 
- *          to be used by KMRK_RunKMRK
- *
- * @see   KMRK_InitOneSequence
- */
-
-vn_type* KMRK_HashOneSequence(char *sequence,              
-                              int32_t complement, 
-                              int32_t count,
-                              int32_t countSeq,
-                              int32_t *k,
-                              KMRK_QUORUM_FUNC_PTR(quorum),
-                              masked_area_table_t *mask);
-
-
-/** 
- * An example of quorum function testing than a factor is
- * present a least two times. Because of definition of this
- * quorum function count and countSeq parametter have no meanning
- * in this instance of quorum function
- * 
- * @param x   a pointer to vn_type structure to check
- * @param pos position in n of the beginning of the linked list to test
- * @param count minimal number of occurence of factor
- * @param countSeq minimal number of sequences concerned
- * 
- * @return 1 if quorum is ok 0 otherwise.
- */
-
-int32_t KMRK_CoupleQuorum(vn_type* x, 
-			  int32_t pos,
-			  int32_t count, 
-			  int32_t countSeq);
-
-
-/** 
- * An example of quorum function testing than a factor is
- * present a least two times in the direct strand of a sequence or
- * at least one time in the direct strand and one time in the reverse
- * strand. Because of definition of this
- * quorum function count and countSeq parametter have no meanning
- * in this instance of quorum function
- * 
- * @param x   a pointer to vn_type structure to check
- * @param pos position in n of the beginning of the linked list to test
- * @param count minimal number of occurence of factor
- * @param countSeq minimal number of sequences concerned
- * 
- * @return 1 if quorum is ok 0 otherwise.
- */
-
-int32_t KMRK_DirInvQuorum(vn_type* x, 
-			  int32_t pos,
-			  int32_t count, 
-			  int32_t countSeq);
-
-/** 
- * An example of quorum function testing than a factor is
- * present a least one time in the direct strand and one time in the reverse
- * strand. Because of definition of this
- * quorum function count and countSeq parametter have no meanning
- * in this instance of quorum function
- * 
- * @param x   a pointer to vn_type structure to check
- * @param pos position in n of the beginning of the linked list to test
- * @param count minimal number of occurence of factor
- * @param countSeq minimal number of sequences concerned
- * 
- * @return 1 if quorum is ok 0 otherwise.
- */
-
-int32_t KMRK_InvQuorum(vn_type* x, 
-		       int32_t pos,
-		       int32_t count, 
-		       int32_t countSeq);
-
-int32_t KMRK_Int12Quorum(vn_type* x, 
-			 int32_t pos,
-			 int32_t count, 
-			 int32_t countSeq);
-
-int32_t KMRK_IntInv12Quorum(vn_type* x, 
-			    int32_t pos,
-			    int32_t count, 
-			    int32_t countSeq);
-
-int32_t KMRK_IntDirInv12Quorum(vn_type* x, 
-			       int32_t pos,
-			       int32_t count, 
-			       int32_t countSeq);
-/** 
- * realize one cycle of KMR.
- * 
- * @param x a pointer to vn_type created by an intialisation
- *           function or returned by this function.
- * @param k step used to join two words
- * @param count parametter count passed to the quorun function
- * @param countSeq parametter countSeq passed to the quorun function
- * @param KMRK_QUORUM_FUNC_PTR quorum pointer to a quorum function
- */
-
-void KMRK_RunKMRK(vn_type *x, 
-		  int32_t k, 
-		  int32_t count,
-		  int32_t countSeq,
-		  KMRK_QUORUM_FUNC_PTR(quorum));
-
-
-
-
-/** 
- * realises serveral run of KMR cycle to make from a sequence
- * a vn_type structure describing sequences of factors of a precise size.
- * 
- * @param sequence pointer to a C string containing the sequence
- * @param size word size to construct
- * @param count parametter count passed to the quorun function
- * @param countSeq parametter countSeq passed to the quorun function
- * @param quorum pointer to a quorum function
- * @param init pointer to a initialisation  function
- * 
- * @return a vn_type pointer to a structure containing sequences of factors
- */
-
-vn_type *KMRK_RunTo(char *sequence,
-                    int32_t size,
-                    int32_t complement,
-                    int32_t count,
-                    int32_t countSeq,
-                    KMRK_QUORUM_FUNC_PTR(quorum),
-                    KMRK_INIT_FUNC_PTR(init,quorum),
-                    masked_area_table_t *mask);
-
-
-
-/** 
- * free memory associated to a vn_type pointer
- * 
- * @param x a pointer to vn_type structure
- */
-
-void KMRK_FreeVN(vn_type *x);
-
-
-int32_t KMRK_upperCoupleCount(vn_type *x);
-int32_t KMRK_upperInvertedCount(vn_type* x,int32_t wordsize);
-int32_t KMRK_upperInterCount(vn_type* x,int32_t seq1,int32_t seq2,int32_t wordsize);
-
-void KMRK_markStart(vn_type* x);
-
-#endif //_KMRK_H_

src/libKMRK/KMRK_Seeds.c

@@ -1,908 +0,0 @@
-#include "KMRK_Seeds.h"
-#include "memory.h"
-#include <stdlib.h>
-#include <string.h>
-#include "sequence.h"
-
-
-static void SetMultipleLenInvSeeds(SmallSeed_type* seeds,
-				   int32_t         nseeds,
-				   int32_t         wordSize,
-				   int8_t          same,
-				   AllSeeds_type *PtrAllSeeds);
-
-
-/*
-	Concatenate
-	DirSeq\0InvSeq\0
-*/
-static char* makeDirInvSeq(char* seq, int32_t size)
-{
-  char *SeqInv;
-
-	seq = (char *)MyRealloc( (void *)seq, (size*2+2)*sizeof(char),
-                            (size+1)* sizeof(char),  "Cannot allocate space for reverse sequence");
-	
-	SeqInv= seq + size + 1;
-	seq[size]=0;
-	invseq(seq, SeqInv);
-	seq[size]='@';
-	SeqInv[size]=0;
-
-  return seq;
-}
-
-/*
-	Merge the seq1 with seq2
-*/
-static char *merge2seq(char* seq1, char* seq2,
-		       int32_t size1, int32_t size2)
-{
-  char * dest;
-
-  seq1 = (char *)MyRealloc((void *)seq1,  (size1+size2+2) *sizeof(char),
-           (size1+1)*sizeof(char), "Cannot allocate space for reverse sequence");
-
-  dest = seq1 + size1 + 1;
-  seq1[size1]='@';
-  memcpy(dest,seq2,size2);
-  dest[size2]=0;
-
-  return seq1;
-}
-
-static int32_t dirDelta(SmallSeed_type* seed)
-{
-    return seed->pos2 - seed->pos1;
-}
-
-void KMRK_SetMultipleLenDirSeeds(SmallSeed_type* seeds,
-                              int32_t         nseeds,
-                              int32_t         wordSize,
-                              AllSeeds_type   *PtrAllSeeds)
-{
-    
-    int32_t i,j;                                /* dummy counters j=kept seeds ; i = current seed */
-    int32_t curLen=wordSize;                    /* Length of the current seed */  
-    int32_t delta;
-    SmallSeed_type *mainSeed;
-    SmallSeed_type *curSeed;
-    int32_t add;
-
-/*    fprintf(stderr,"New Version\n");*/
-
-    KMRK_sortSeeds(seeds,nseeds,KMRK_cmpDeltaSeedsPos);
-
-    for(j=0,mainSeed=seeds ; 
-        j < nseeds; 
-        j=i,mainSeed=curSeed)
-    {    
-        /* foreach seed */
-        delta = dirDelta(mainSeed);
-        curLen=wordSize;
-        
-        for (i=j+1,curSeed=mainSeed+1; 
-             i < nseeds && 
-             dirDelta(curSeed)==delta &&
-             (curSeed->pos1 - mainSeed->pos1) <= curLen;
-             i++,curSeed++)
-        {
-            add=wordSize - mainSeed->pos1 - curLen + curSeed->pos1;
-            if (add < 0) add = 0;
-            curLen+=add;
-        }   
-        
-        KMRK_pushSeed(PtrAllSeeds,
-                      seeds[j].pos1,seeds[j].pos2,
-                      curLen,
-                      1);
-    }
-}
-
-static int32_t invDelta(SmallSeed_type* seed)
-{
-    return seed->pos2 + seed->pos1;
-}
-
-static void SetMultipleLenInvSeeds(SmallSeed_type* seeds,
-				   int32_t         nseeds,
-				   int32_t         wordSize,
-				   int8_t          same,
-				   AllSeeds_type *PtrAllSeeds)
-{
-
-    int32_t i,j;                                /* dummy counters j=kept seeds ; i = current seed */
-    int32_t curLen=wordSize;                    /* Length of the current seed */  
-    int32_t delta;
-    int32_t pos2;
-    SmallSeed_type *mainSeed;
-    SmallSeed_type *curSeed;
-    int32_t add;
-
-/*    fprintf(stderr,"New Version\n");*/
-   
-    KMRK_sortSeeds(seeds,nseeds,KMRK_cmpDeltaInvSeedsPos);
-    
-    for(j=0,mainSeed=seeds ; 
-        j < nseeds; 
-        j=i,mainSeed=curSeed)
-    {    
-        /* foreach seed */
-        delta = invDelta(mainSeed);
-        curLen=wordSize;
-        
-        for (i=j+1,curSeed=mainSeed+1; 
-             i < nseeds && 
-             invDelta(curSeed)==delta &&
-             (curSeed->pos1 - mainSeed->pos1) <= curLen;
-             i++,curSeed++)
-        {
-            add=wordSize - mainSeed->pos1 - curLen + curSeed->pos1;
-            if (add < 0) add = 0;
-            curLen+=add;
-        }   
-        
-        if ( same && (seeds[j].pos1+curLen>= seeds[i-1].pos2))
-        {
-            curLen = (seeds[i-1].pos2 - seeds[j].pos1 + 1) * 2;
-            pos2 = seeds[j].pos1;
-        }
-        else
-             pos2 = seeds[i-1].pos2;
-             
-        KMRK_pushSeed(PtrAllSeeds,
-                      seeds[j].pos1,pos2,
-                      curLen,
-                      0);
-    }
-}
-
-AllSeeds_type *KMRK_allocSeeds(AllSeeds_type *AllSeeds, 
-			       int32_t size, 
-			       int8_t opt_dir, 
-			       int8_t opt_inv)
-{
-
-  AllSeeds_type *reponse;
-
-  if(opt_inv != 1 && opt_dir != 1)
-    {
-      fprintf(stderr,"AllocSeeds: requiere at least "
-	      "one of opt_dir and opt_inv to be 1\n");
-      exit(4);
-    }
-
-  reponse = AllSeeds;
-
-  if (!reponse)
-     reponse = MyCalloc( 1, sizeof(AllSeeds_type),"KMRK_allocSeeds: cannot allocate new data structure");
-
-	if(opt_dir)
-	{
-		if (reponse->dirSeeds==NULL)
-			reponse->dirSeeds = (Seed_type *)MyCalloc( size,sizeof(Seed_type),"KMRK_allocSeeds: cannot allocate new data structure");
-      else
-		{
-			if(size)
-				reponse->dirSeeds = (Seed_type *)MyRealloc( (void *)reponse->dirSeeds,
-				                     size*sizeof(Seed_type),  reponse->cDirSeeds*sizeof(Seed_type),
-											"allocSeeds: cannot reallocate data structure" );
-			else
-	  		{
-				MyFree(reponse->dirSeeds, reponse->cDirSeeds*sizeof(Seed_type));
-	    		reponse->dirSeeds=NULL;
-	    		reponse->cDirSeeds=0;
-			}
-		}
-
-		reponse->cDirSeeds=size;
-	}
-	
-	
-	if(opt_inv)
-	{
-		if (reponse->invSeeds==NULL)
-			reponse->invSeeds = (Seed_type *)MyCalloc( size, sizeof(Seed_type),"KMRK_allocSeeds: cannot allocate new data structure");
-      else
-		{
-			if(size)
-				reponse->invSeeds = (Seed_type *)MyRealloc( (void *)reponse->invSeeds,
-				                     size*sizeof(Seed_type),  reponse->cInvSeeds*sizeof(Seed_type),
-											"allocSeeds: cannot reallocate data structure" );
-			else
-	  		{
-				MyFree(reponse->invSeeds, reponse->cInvSeeds*sizeof(Seed_type));
-	    		reponse->invSeeds=NULL;
-	    		reponse->cInvSeeds=0;
-			}
-		}
-
-		reponse->cInvSeeds=size;
-	}
-	
-      
-  return reponse;
-}
-
-void KMRK_freeSeeds(AllSeeds_type *AllSeeds)
-{
-  if (!AllSeeds)
-    return;
-
-  if (AllSeeds->dirSeeds)
-    MyFree(AllSeeds->dirSeeds, AllSeeds->cDirSeeds*sizeof(Seed_type) );
-	AllSeeds->dirSeeds=NULL;
-
-  if (AllSeeds->invSeeds)
-    MyFree(AllSeeds->invSeeds, AllSeeds->cInvSeeds*sizeof(Seed_type) );
-	AllSeeds->dirSeeds=NULL;
-
-  MyFree(AllSeeds, 1*sizeof( AllSeeds_type ) );
-}
-
-void KMRK_compactSeeds(AllSeeds_type *AllSeeds)
-{
-  if (AllSeeds)
-    {
-      if (AllSeeds->dirSeeds)
-	KMRK_allocSeeds(AllSeeds,
-			AllSeeds->nDirSeeds,
-			1,
-			0);
-
-      if (AllSeeds->invSeeds)
-	KMRK_allocSeeds(AllSeeds,
-			AllSeeds->nInvSeeds,
-			0,
-			1);
-
-    }
-}
-
-
-void KMRK_pushSeed(AllSeeds_type *AllSeeds,
-		   int32_t       pos1,
-		   int32_t       pos2,
-		   int32_t       length,
-		   int8_t        dir)
-{
-  Seed_type* stack;
-  int32_t    maxcount;
-  int32_t    index;
-
-  if (dir)
-    {
-      dir      = 1;
-      stack    = AllSeeds->dirSeeds;
-      maxcount = AllSeeds->cDirSeeds;
-      index    = AllSeeds->nDirSeeds;
-    }
-  else
-    {
-      dir      = 0;
-      stack    = AllSeeds->invSeeds;
-      maxcount = AllSeeds->cInvSeeds;
-      index    = AllSeeds->nInvSeeds;
-    }
-
-  if (index == maxcount)
-    {
-      (void) KMRK_allocSeeds(AllSeeds,
-			     maxcount * 2,
-			     dir,
-			     !dir);
-      
-      if (dir)
-	stack    = AllSeeds->dirSeeds;
-      else
-	stack    = AllSeeds->invSeeds;
-    }
-
-  stack+=index;
-
-  stack->pos1   = pos1;
-  stack->pos2   = pos2;
-  stack->length = length;
- 
-  if (dir)
-    AllSeeds->nDirSeeds++;
-  else
-    AllSeeds->nInvSeeds++;
-
-}
-
-
-
-
-AllSeeds_type* KMRK_enumerateDirectCouple(AllSeeds_type* Seeds,
-					  int32_t expected,
-					  int32_t wordSize,
-					  vn_type* stack,
-					  int32_t seq)
-{
-
-  int32_t xmin;
-  int32_t xmax;
-  int32_t i;
-  int32_t j;
-  int32_t k;
-  int32_t next;
-  int32_t* n;
-  int32_t nseeds;
-
-  SmallSeed_type *list;
-
-  list = (SmallSeed_type *)MyCalloc( expected, sizeof(SmallSeed_type) ,
-	      "KMRK_enumerateDirectCouple: cannot allocate DirectCouple memory");
-
-  nseeds = 0;
-
-  n    = GETVECTOR(stack,n);
-  
-  if (seq)
-    xmin = stack->limits[seq-1];
-  else
-    xmin = 0;
-
-  xmax = stack->limits[seq];
-
-  for (i=1; i <= xmax; i++)
-	  if (IS_MARKED(n,i))     /* Check begining of chained list */
-      {
-				/* Look for begining of sequence of interest */   
-		  for( ;(i <= xmin) && (i != GET(n,i));
-			   i=GET(n,i));
-		  
-				/* for each factor in  sequence of interest */   
-		  for (j=i; 
-			   (j != GET(n,j)) && (j <= xmax);
-			   j = GET(n,j))
-		  {
-			  next = GET(n,j);
-			  if (next <= xmax)
-				  do
-				  {
-					  k = next;
-					  next = GET(n,k);
-					  list[nseeds].pos1 = j-1;
-					  list[nseeds].pos2 = k-1;
-					  nseeds++;
-				  } while ((k!=next) && (next <= xmax));
-		  }
-      } ;
-  
-  Seeds = KMRK_allocSeeds(Seeds,
-			  expected/20+1,
-			  1,0);
-
-/*  fprintf(stderr,"Expected direct couple : %d\n",expected);*/
-  
-  KMRK_SetMultipleLenDirSeeds(list,nseeds,wordSize,Seeds);
-  MyFree(list, expected*sizeof(SmallSeed_type) );
-  KMRK_compactSeeds(Seeds);
-	  
-  return Seeds;
-}
-
-
-/*
-	From KMR-K Stacks to SmallSeeds
-*/
-AllSeeds_type* KMRK_enumerateInvertedCouple(AllSeeds_type* Seeds,
-					    int32_t expected,
-					    int32_t wordSize,
-					    vn_type* stack)
-{
-  int32_t xmax;
-  int32_t invmax;
-  int32_t posinv;
-  int32_t i;
-  int32_t j;
-  int32_t k;
-  int32_t memk;
-  int32_t* n;
-  int32_t next;
-
-  int32_t nseeds;        /* number of seeds */
-
-  SmallSeed_type *list;   /* seed list with only pos1 and pos2 --simple output from kmrk */
-
-
-  list = (SmallSeed_type *)MyCalloc( expected, sizeof(SmallSeed_type) ,
-	      "KMRK_enumerateInvertedCouple: cannot allocate InvertedCouple memory");
-
-  nseeds = 0;
-
-  n    = GETVECTOR(stack,n);
-  
-  xmax = stack->limits[0];
-  invmax = stack->limits[1];
-
-  for (i=1; i <= xmax; i++)
-	  if (IS_MARKED(n,i))
-      {
-		  for(memk=i ;
-			  (memk < xmax) && 
-			  memk != GET(n,memk) &&
-			  (memk <= invmax);
-			  memk=GET(n,memk));
-		  
-		  if ((memk > xmax) && (memk <= invmax))
-			  for (j=i;
-				   (j <= xmax) && (j != GET(n,j));
-				   j=GET(n,j))
-			  {
-				  next = memk;
-				  do
-				  {
-					  k = next;
-					  next = GET(n,k);
-					  posinv = 2 * xmax - k -wordSize + 3;
-					  if (j <= posinv)
-					  {
-						  list[nseeds].pos1=j-1;
-						  list[nseeds].pos2=posinv-1;
-						  nseeds++;
-					  }
-				  } while ((next <= invmax) && (k != next));
-			  }
-      }
-		  
-		  Seeds = KMRK_allocSeeds(Seeds,
-			  expected/20+1,
-			  0,1);
-
-/*  fprintf(stderr,"Expected inverted couple : %d\n",expected);*/
-
-  SetMultipleLenInvSeeds(list,nseeds,wordSize,1,Seeds);   /* turn the Small seeds into merged seeds */
-  MyFree(list, expected*sizeof(SmallSeed_type) );
-  KMRK_compactSeeds(Seeds);
-
-
-  return Seeds;
-}
-
-AllSeeds_type* KMRK_enumerateInterCouple(AllSeeds_type* Seeds,
-					 int32_t seq1,
-					 int32_t seq2,
-					 int32_t expected,
-					 int32_t wordSize,
-					 vn_type* stack)
-{
-  int32_t xmin;
-  int32_t xmax;
-  int32_t ymax;
-  int32_t ymin;
-  int32_t pos1;
-  int32_t pos2;
-  int32_t i;
-  int32_t j;
-  int32_t k;
-  int32_t memj;
-  int32_t memk;
-  int32_t* n;
-  int32_t next;
-
-  int32_t nseeds;
-
-
-  SmallSeed_type *list;
-  nseeds=0;
-
-  list = (SmallSeed_type *)MyCalloc( expected, sizeof(SmallSeed_type) ,
-	      "KMRK_enumerateInterCouple: cannot allocate InterCouple memory");
-
-  n    = GETVECTOR(stack,n);
-  
-  if (seq1==0)
-    xmin=0;
-  else
-    xmin = stack->limits[seq1-1];
-
-  xmax = stack->limits[seq1];
-  ymin = stack->limits[seq2-1];
-  ymax = stack->limits[seq2];
-
-  for (i=1; i <= xmax; i++)
-	  if (IS_MARKED(n,i))
-      {
-		  for(memj=i ;
-			  (memj < xmin) && 
-			  memj != GET(n,memj);
-			  memj=GET(n,memj));
-		  
-		  if ((memj > xmin) && (memj <= xmax))
-		  {
-			  for(memk=memj ;
-				  (memk < ymin) && 
-				  memk != GET(n,memk);
-				  memk=GET(n,memk));
-			  
-			  if ((memk > ymin) && (memk <= ymax))
-				  for (j=memj;
-					   (j <= xmax) && (j != GET(n,j));
-					   j=GET(n,j))
-				  {
-					  next = memk;
-					  do
-					  {
-						  k = next;
-						  next = GET(n,k);
-						  if (seq1 > 0)
-							  pos1 = j - xmin - 1;
-						  else
-							  pos1 = j;
-						  pos2 = k - ymin - 1;
-						  list[nseeds].pos1 = pos1 - 1;
-						  list[nseeds].pos2 = pos2 - 1;
-						  nseeds++;
-					  } while ((next <= ymax) && (k != next));
-				  }
-		  }
-      }
-
-	  
-	  Seeds = KMRK_allocSeeds(Seeds,
-							  expected/20+1,
-							  1,0);
-	  
-	/*  fprintf(stderr,"Expected inter-direct couple : %d\n",expected);*/
-	  KMRK_SetMultipleLenDirSeeds(list,nseeds,wordSize,Seeds);
-	  
-	  MyFree(list, expected*sizeof(SmallSeed_type) );
-	  
-	  KMRK_compactSeeds(Seeds);
-	  
-	  return Seeds;
-}
-
-AllSeeds_type* KMRK_enumerateInterInvertedCouple(AllSeeds_type* Seeds,
-						 int32_t seq2,
-						 int32_t expected,
-						 int32_t wordSize,
-						 vn_type* stack)
-
-{
-  int32_t xmin;
-  int32_t xmax;
-  int32_t ymax;
-  int32_t ymin;
-  int32_t posinv;
-  int32_t pos2;
-  int32_t i;
-  int32_t j;
-  int32_t k;
-  int32_t memj;
-  int32_t memk;
-  int32_t* n;
-  int32_t next;
-
-  int32_t nseeds;
-
-  SmallSeed_type *list;
-
-  list = (SmallSeed_type *)MyCalloc( expected, sizeof(SmallSeed_type) ,
-	      "KMRK_enumerateInterCouple: cannot allocate InterCouple memory");
-
-
-  nseeds = 0;
-  n    = GETVECTOR(stack,n);
-
-  if (seq2 < 2)
-    {
-      fprintf(stderr,"enumerateInterInvertedCouple: seq2 must be differente to 0");
-      exit(4);
-    }
-    
-  
-  xmin = stack->limits[0];
-  xmax = stack->limits[1];
-
-  ymin = stack->limits[seq2-1];
-  ymax = stack->limits[seq2];
-
-  Seeds = KMRK_allocSeeds(Seeds,
-			  expected,
-			  0,1);
-
-  for (i=1; i <= xmax; i++)
-	  if (IS_MARKED(n,i))
-      {
-		  for(memj=i ;
-			  (memj < xmin) && 
-			  memj != GET(n,memj);
-			  memj=GET(n,memj));
-		  
-		  if ((memj > xmin) && (memj <= xmax))
-		  {
-			  for(memk=memj ;
-				  (memk < ymin) && 
-				  memk != GET(n,memk);
-				  memk=GET(n,memk));
-			  
-			  if ((memk > ymin) && (memk <= ymax))
-				  for (j=memj;
-					   (j <= xmax) && (j != GET(n,j));
-					   j=GET(n,j))
-				  {
-					  next = memk;
-					  do
-					  {
-						  k = next;
-						  next = GET(n,k);
-						  posinv = 2 * xmin - j -wordSize + 3;
-						  pos2 = k - ymin - 1;
-						  list[nseeds].pos1=posinv-1;
-						  list[nseeds].pos2=pos2-1;
-						  nseeds++;
-					  } while ((next <= ymax) && (k != next));
-				  }
-		  }
-      }
-
-		  Seeds = KMRK_allocSeeds(Seeds,
-			  expected/20+1,
-			  0,1);
-
-/*  fprintf(stderr,"Expected inter-inverted couple : %d\n",expected);*/
-  SetMultipleLenInvSeeds(list,nseeds,wordSize,0,Seeds);
-  MyFree(list, expected* sizeof(SmallSeed_type) );
-  KMRK_compactSeeds(Seeds);
-
-		  
-  return Seeds;
-}
-
-  
-int32_t KMRK_cmpSeedsPos(SmallSeed_type *s1, SmallSeed_type *s2)
-{
-  if (s1->pos1==s2->pos1)
-    return s1->pos2 - s2->pos2;
-  else
-    return s1->pos1 - s2->pos1;
-}
-
-int32_t KMRK_cmpDeltaSeedsPos(SmallSeed_type *s1, SmallSeed_type *s2)
-{
-    int32_t delta1 =  s1->pos2-s1->pos1;
-    int32_t delta2 =  s2->pos2-s2->pos1;
-    
-    if (delta1==delta2)
-        return s1->pos1 - s2->pos1;
-    else
-        return delta1 - delta2;
-}
-
-int32_t KMRK_cmpDeltaInvSeedsPos(SmallSeed_type *s1, SmallSeed_type *s2)
-{
-    int32_t delta1 =  s1->pos2+s1->pos1;
-    int32_t delta2 =  s2->pos2+s2->pos1;
-    
-    if (delta1==delta2)
-        
-        return s1->pos1 - s2->pos1;
-    else
-        return delta1 - delta2;
-}
-
-
-void KMRK_sortSeeds(SmallSeed_type* seeds,
-		    int32_t         nseeds,
-		    KMRK_SORT_SEEDS_FUNC_PTR(compare))
-{
-
-  qsort(seeds, 
-	nseeds, 
-	sizeof(SmallSeed_type),
-	  (int (*)(const void *, const void *))compare);
-}
-
-AllSeeds_type* KMRK_get_seeds(char **seq,
-                              int32_t SimpleSeqLen,
-                              int16_t Lmin,
-                              int8_t opt_dir,
-                              int8_t opt_inv, 
-                              int8_t opt_verbose,
-                              masked_area_table_t *mask) 
-{
-  AllSeeds_type* AllSeeds;
-  char *SeqDir = *seq;
-  vn_type * stacks;
-  int32_t dirExpect=0;
-  int32_t invExpect=0;
-  KMRK_QUORUM_FUNC_PTR(quorum);
-
-  if(opt_inv != 1 && 
-     opt_dir != 1)
-    {
-      fprintf(stderr,
-	      "get_seeds: requiered at least "
-	      "opt_dir or opt_inv to be 1\n");
-      exit(4);
-    }
-	
-  if(opt_inv)
-      SeqDir = makeDirInvSeq(SeqDir,SimpleSeqLen);   /* create a sequence with "DirSeq\0InvSeq\0" */
-
-	if (opt_inv){                                     /* Are we interested in dir, inv or both ? */
-		if (opt_dir)
-			quorum = KMRK_DirInvQuorum;
-		else
-			quorum = KMRK_InvQuorum;
-	}
-	else
-		quorum = KMRK_CoupleQuorum;
-
-  stacks = KMRK_RunTo(SeqDir,
-		      Lmin,
-		      opt_inv,
-		      2,
-		      1,
-		      quorum,
-		      KMRK_HashOneSequence,
-              mask);
-
-  invExpect =0;
-
-  KMRK_markStart(stacks);
-
-  if (opt_inv)
-    {
-      SeqDir = (char *)MyRealloc( (void *)SeqDir,  (SimpleSeqLen+1)*sizeof(char),       
-		                            (2*SimpleSeqLen+2)*sizeof(char) , "KRMK_get_seeds: Cannot shrink memory");   /* reset mem to a sigle sequence */
-      SeqDir[SimpleSeqLen]=0;
-    }
-
-    if(opt_inv)
-        invExpect = KMRK_upperInvertedCount(stacks,Lmin);
-
-    if(opt_dir)
-        dirExpect = KMRK_upperCoupleCount(stacks);
- 
-
-  AllSeeds = NULL;
-
-	
-  MyFree(stacks->v, stacks->size * sizeof(int32_t));
-  stacks->v=NULL;
-
-
-  if (opt_dir)
-    AllSeeds = KMRK_enumerateDirectCouple(AllSeeds,dirExpect,Lmin ,stacks,0);
-
-  if (opt_inv)
-    AllSeeds = KMRK_enumerateInvertedCouple(AllSeeds,invExpect,Lmin,stacks);
-
-  KMRK_FreeVN(stacks);
-
-  *seq = SeqDir;
-				 
-  return AllSeeds;
-}  
-
-AllSeeds_type* KMRK_get_seeds_2seqs(char **seq1,
-                                    char **seq2,
-                                    int32_t size1,
-                                    int32_t size2,
-                                    int16_t Lmin,
-                                    int8_t opt_dir,
-                                    int8_t opt_inv, 
-                                    int8_t opt_verbose,
-                                    masked_area_table_t *mask) 
-{
-
-  AllSeeds_type* AllSeeds;
-  char *sequence1 = *seq1;
-  char *sequence2 = *seq2;
-  vn_type * stacks;
-  int32_t dirExpect=0;
-  int32_t invExpect=0;
-  KMRK_QUORUM_FUNC_PTR(quorum);
-  int32_t sizef; 
-
-  if(opt_inv != 1 && 
-     opt_dir != 1)
-    {
-      fprintf(stderr,
-	      "get_seeds_2seqs: requiered at least "
-	      "opt_dir or opt_inv to be 1\n");
-      exit(4);
-    }
-
-  sizef = size1;
-
-  if(opt_inv)
-    {
-      sequence1 = makeDirInvSeq(sequence1,size1);
-      sizef+=(1+size1);
-    }
-
-  sequence1 = merge2seq(sequence1,sequence2,sizef,size2);
-
-  if (opt_inv)
-    if (opt_dir)
-      quorum = KMRK_IntDirInv12Quorum;
-    else
-      quorum = KMRK_IntInv12Quorum;
-  else
-    quorum = KMRK_Int12Quorum;
-
-  stacks = KMRK_RunTo(sequence1,
-		      Lmin,
-		      opt_inv,
-		      2,
-		      2,
-		      quorum,
-		      KMRK_HashOneSequence,
-              mask);
-
-  KMRK_markStart(stacks);
-
-  sequence1= (char *)MyRealloc(
-  			(void *)sequence1,
-			(size1+1)*sizeof(char),
-	      (sizef+size2+2)*sizeof(char),
-	       "KMRK_get_seeds_2seqs: shrink memory from 3N to 1N... ???");
-
-  sequence1[size1]=0;
-
-  if (opt_dir){
-    if (opt_inv)
-      dirExpect = KMRK_upperInterCount(stacks,0,2,Lmin);
-    else
-      dirExpect = KMRK_upperInterCount(stacks,0,1,Lmin);
-	}
-  if (opt_inv)
-    invExpect = KMRK_upperInterCount(stacks,1,2,Lmin);
-
-  AllSeeds = NULL;
-  MyFree(stacks->v, stacks->size*sizeof(int32_t));
-  stacks->v=NULL;
-
-  if (opt_dir){
-    if (opt_inv)
-      AllSeeds = KMRK_enumerateInterCouple(AllSeeds,
-					   0,2,
-					   dirExpect,
-					   Lmin ,
-					   stacks);
-    else
-      AllSeeds = KMRK_enumerateInterCouple(AllSeeds,
-					   0,1,
-					   dirExpect,
-					   Lmin ,
-					   stacks);
-
-	}
-  if (opt_inv)
-    AllSeeds = KMRK_enumerateInterInvertedCouple(AllSeeds,
-						 2,
-						 invExpect,
-						 Lmin ,
-						 stacks);
-
-
-  KMRK_FreeVN(stacks);
-
-  *seq1 = sequence1;
-				 
-  return AllSeeds;
-}
-
-#define SIGN(x) (((x)<0) ? -1:(((x)>0) ? 1:0))
-
-static int32_t compareSeedsByPos(Seed_type* s1,Seed_type* s2)
-{
-  if (s1->pos1 == s2->pos1)
-    return SIGN(s1->pos2 - s2->pos2);
-  else
-    return SIGN(s1->pos1 - s2->pos1);
-}
-
-void KMRK_sortSeedsByPos(Seed_type* seeds, int32_t count)
-{
-  qsort(seeds,
-	count,
-	sizeof(Seed_type),
-	(int (*)(const void *, const void *))compareSeedsByPos);
-};

src/libKMRK/KMRK_Seeds.h

@@ -1,126 +0,0 @@
-#ifndef KMRK_Seeds_h
-#define KMRK_Seeds_h
-
-/********************************************
- ********************************************
- **
- ** Declaration of struct
- **
- ********************************************
- ********************************************/
-
-#include <stdint.h>
-#include <stdio.h>
-#include "KMRK.h"
-
-
-
-
-#define KMRK_SORT_SEEDS_FUNC_PTR(name) int32_t (*name)(SmallSeed_type*, \
-						       SmallSeed_type*)
-
-#define KMRK_DELTA_SEEDS_FUNC_PTR(name) int32_t (*name)(SmallSeed_type*)
-
-/********************************************
- ********************************************
- **
- ** Declaration of public functions
- **
- ********************************************
- ********************************************/
-
-
-AllSeeds_type *KMRK_allocSeeds(AllSeeds_type *AllSeeds, 
-			       int32_t size, 
-			       int8_t opt_dir, 
-			       int8_t opt_inv);
-
-void KMRK_SetMultipleLenDirSeeds(SmallSeed_type* seeds,
-                                 int32_t         nseeds,
-                                 int32_t         wordSize,
-                                 AllSeeds_type   *PtrAllSeeds);
-
-
-void KMRK_freeSeeds(AllSeeds_type *AllSeeds);
-
-void KMRK_compactSeeds(AllSeeds_type *AllSeeds);
-
-void KMRK_pushSeed(AllSeeds_type *AllSeeds,
-		   int32_t       pos1,
-		   int32_t       pos2,
-		   int32_t       length,
-		   int8_t        dir);
-
-AllSeeds_type* KMRK_enumerateDirectCouple(AllSeeds_type* Seeds,
-					  int32_t expected,
-					  int32_t wordSize,
-					  vn_type* stack,
-					  int32_t seq);
-
-AllSeeds_type* KMRK_enumerateInvertedCouple(AllSeeds_type* Seeds,
-					    int32_t expected,
-					    int32_t wordSize,
-					    vn_type* stack);
-
-AllSeeds_type* KMRK_enumerateInterCouple(AllSeeds_type* Seeds,
-					 int32_t seq1,
-					 int32_t seq2,
-					 int32_t expected,
-					 int32_t wordSize,
-					 vn_type* stack);
-
-AllSeeds_type* KMRK_enumerateInterInvertedCouple(AllSeeds_type* Seeds,
-						 int32_t seq2,
-						 int32_t expected,
-						 int32_t wordSize,
-						 vn_type* stack);
-
-
-/** 
- * Compare two seeds and return an integer less than, equal to or greater 
- * than zero considering the relative order of the two seeds. This
- * version take into account only pos1 and pos2 of seeds without taking
- * account of the sequences or of the relative direction
- * 
- * @param s1 pointer to seed one
- * @param s2 pointer to seed two
- * 
- * @return a integer less than, equal to or greater than zero
- */
-
-int32_t KMRK_cmpSeedsPos(SmallSeed_type *s1, SmallSeed_type *s2);
-int32_t KMRK_cmpDeltaSeedsPos(SmallSeed_type *s1, SmallSeed_type *s2);
-int32_t KMRK_cmpDeltaInvSeedsPos(SmallSeed_type *s1, SmallSeed_type *s2);
-
-void KMRK_sortSeeds(SmallSeed_type* seeds,
-		    int32_t         nseeds,
-		    KMRK_SORT_SEEDS_FUNC_PTR(compare));
-
-
-AllSeeds_type* KMRK_get_seeds(char **seq,
-                              int32_t SimpleSeqLen,
-                              int16_t Lmin,
-                              int8_t opt_dir,
-                              int8_t opt_inv, 
-                              int8_t opt_verbose,
-                              masked_area_table_t *mask);
-
-AllSeeds_type* KMRK_get_seeds_2seqs(char **seq1,
-                                    char **seq2,
-                                    int32_t size1,
-                                    int32_t size2,
-                                    int16_t Lmin,
-                                    int8_t opt_dir,
-                                    int8_t opt_inv, 
-                                    int8_t opt_verbose,
-                                    masked_area_table_t *mask);
-
-/** 
- * Order an array of seeds by pos1,pos2
- * 
- * @param seeds pointer to an array of Seed_type object to sort
- * @param count count of element in the array
- */
-void KMRK_sortSeedsByPos(Seed_type* seeds, int32_t count);
-
-#endif   /* KMRK_Seeds_h */

src/libKMRK/KMRK_mask.c

@@ -1,259 +0,0 @@
-/*
- *  KMRK_mask.c
- *  repseek
- *
- *  Created by Eric Coissac on 04/12/04.
- *  Copyright 2004 __MyCompanyName__. All rights reserved.
- *
- */
-
-#include "KMRK_mask.h"
-#include <stdio.h>
-#include <stdlib.h>
-#include "memory.h"
-
-#define MASKED_AREA_TABLE_SIZE(seqcount) (sizeof(masked_area_table_t) + (sizeof(masked_area_list_t*) * ((seqcount)-1)))
-#define MASKED_AREA_LIST_SIZE(areacount) (sizeof(masked_area_list_t) + (sizeof(masked_area_t) * ((areacount)-1)))
-
-#define AREA_COUNT_INIT (1000)
-
-static masked_area_table_t *new_masked_area_table(int32_t seqcount, int32_t areacount);
-static masked_area_list_t *new_masked_area_list(int32_t areacount);
-static masked_area_list_t *realloc_masked_area_list(masked_area_list_t *list,int32_t areacount);
-static int32_t push_area(masked_area_table_t* table,int32_t sequence,int32_t begin,int32_t end);
-static void sort_area_table(masked_area_table_t* table);
-static int32_t compare_area(const masked_area_t* v1,const masked_area_t* v2);
-static int32_t search_area(const masked_area_t* v1,const masked_area_t* v2);
-
-static masked_area_list_t *strip_area_list(masked_area_list_t* list);
-static void strip_area_table(masked_area_table_t* table);
-
-
-static masked_area_list_t *new_masked_area_list(int32_t areacount)
-{
-    masked_area_list_t *list;
-    
-    list = MyCalloc(1,MASKED_AREA_LIST_SIZE(areacount),"Not enougth memory for mask table");
-    list->reserved=areacount;
-    
-    return list;
-}
-
-static masked_area_list_t *realloc_masked_area_list(masked_area_list_t *list,int32_t areacount)
-{    
-    list = MyRealloc(list,
-                     MASKED_AREA_LIST_SIZE(areacount),
-                     MASKED_AREA_LIST_SIZE(list->reserved),
-                     "Not enougth memory for mask table");
-
-    list->reserved=areacount;
-    
-    return list;
-}
-
-static masked_area_table_t *new_masked_area_table(int32_t seqcount, int32_t areacount)
-{
-    masked_area_table_t *table;
-    int32_t i;
-    
-    table = MyCalloc(1,MASKED_AREA_TABLE_SIZE(seqcount),"Not enougth memory for mask table");
-    table->seqcount=seqcount;
-    
-    for (i=0;i<seqcount;i++)
-        table->sequence[i]=new_masked_area_list(areacount);
-    
-    return table;
-}
-
-static int32_t push_area(masked_area_table_t* table,int32_t sequence,int32_t begin,int32_t end)
-{
-    masked_area_list_t * list;
-    
-    if (sequence >= table->seqcount)
-        return -1;
-    
-    list = table->sequence[sequence];
-
-    if (list->reserved == list->count)
-    {
-        list = realloc_masked_area_list(list,list->reserved*2);
-        table->sequence[sequence]=list;
-    }
-    
-    list->area[list->count].begin=begin;
-    list->area[list->count].end=end;
-    
-    list->count++;
-    table->total++;
-
-    return table->total;
-}
-
-static int32_t compare_area(const masked_area_t* v1,const masked_area_t* v2)
-{
-    return v1->begin - v2->begin;
-}
-
-static void sort_area_table(masked_area_table_t* table)
-{
-    int32_t i;
-    
-    for (i=0; i<table->seqcount;i++)
-    {
-        qsort(table->sequence[i]->area,
-              table->sequence[i]->count,
-              sizeof(masked_area_t),
-              (int (*)(const void *, const void *))compare_area);
-    }
-}
-
-static masked_area_list_t *strip_area_list(masked_area_list_t* list)
-{
-    int32_t i;
-    int32_t j;
-    int32_t count;
-    int32_t newcount;
-    
-    count = list->count;
-    newcount=count;
-    
-    for (i=1;i<count;i++)
-    {
-       /* fprintf(stderr,"\n%d->%d    %d->%d ==>",list->area[i-1].begin,list->area[i-1].end,list->area[i].begin,list->area[i].end); */
-        if ((list->area[i].begin-1) <= list->area[i-1].end)
-        {
-            /* fprintf(stderr," joined"); */
-            list->area[i].begin=list->area[i-1].begin;
-            list->area[i-1].begin=-1;
-            newcount--;
-        }
-    }
-    list->count=newcount;
-
-    for (i=0,j=0;i<count;i++)
-    {
-        if (list->area[i].begin>=0)
-        {
-            if (i!=j)
-                list->area[j]=list->area[i];
-            j++;
-        }
-    }
-    
-    return realloc_masked_area_list(list,newcount);
-
-}
-
-static void strip_area_table(masked_area_table_t* table)
-{
-    int32_t seq;
-    int32_t oldcount;
-    masked_area_list_t* list;
-    
-    sort_area_table(table);
-    
-    for (seq=0; seq < table->seqcount;seq++)
-    {
-        list = table->sequence[seq];
-        oldcount = list->count;
-        table->sequence[seq]=strip_area_list(list);
-        table->total-=oldcount - table->sequence[seq]->count;
-    }
-}
-
-
-static int32_t search_area(const masked_area_t* v1,const masked_area_t* v2)
-{
-    int32_t pos;
-    
-    pos = v1->begin;
-    
-    if (pos < v2->begin)
-            return -1;
-    
-    if (pos > v2->end)
-            return 1;
-    
-    return 0;
-}
-
-
-masked_area_table_t *KMRK_ReadMaskArea(char *areafile,int32_t seqcount,int32_t complement)
-{
-    FILE* area;
-    char buffer[1000];
-    char *ok;
-    int32_t begin;
-    int32_t end;
-    int32_t sequence;
-    int32_t column;
-    int32_t linecount;
-    masked_area_table_t *table;
-    
-    if (complement > 0)
-        seqcount++;
-    else
-        complement=0;
-    
-    area = fopen(areafile,"r");
-    linecount=0;
-    table=new_masked_area_table(seqcount,AREA_COUNT_INIT);
-    
-    do {
-        linecount++;
-        ok = fgets(buffer,999,area);
-        if (ok)
-        {
-            column = sscanf(buffer,"%d %d %d",&begin,&end,&sequence);
-            if (column > 1 && begin <= end)
-            {
-                begin--;
-                end--;
-                if (column==3)
-                    sequence--;
-                else
-                    sequence=0;
-                
-                if (sequence && complement)
-                    sequence++;
-                
-                push_area(table,sequence,begin,end);
-                
-                if (!sequence && complement)
-                    push_area(table,1,complement -1 - end,complement -1 -begin);
-                
-            }
-            
-            if (column==1)
-                fprintf(stderr,"WARNING in mask file reading line %d\n",linecount);
-            
-        }
-    } while (ok);
-    
-    fprintf(stderr,"\nread %d masked areas from file\n",table->total);
-    strip_area_table(table);
-    fprintf(stderr,"strip to %d non overlaping areas\n",table->total);
-    
-    return table;   
-}
-
-char KMRK_isMasked(masked_area_table_t *mask,int32_t seq, int32_t position)
-{
-    masked_area_t input;
-    int32_t result;
-    masked_area_list_t *list;
-    
-    if (! mask || (seq >= mask->seqcount))
-        return 0;
-    
-    list = mask->sequence[seq];
-    input.begin=position;
-    
-    result = bsearch(&input,
-                     list->area,
-                     list->count,
-                     sizeof(masked_area_t),
-                     (int (*)(const void *, const void *))search_area) != NULL;
-
-        return result;
-}

src/libKMRK/KMRK_mask.h

@@ -1,37 +0,0 @@
-/*
- *  KMRK_mask.h
- *  repseek
- *
- *  Created by Eric Coissac on 04/12/04.
- *  Copyright 2004 __MyCompanyName__. All rights reserved.
- *
- */
-
-#include <stdint.h>
-
-#ifndef _KMRK_MASK_H_
-#define _KMRK_MASK_H_
-
-typedef struct {
-    int32_t begin;
-    int32_t end;
-} masked_area_t;
-
-typedef struct {
-    int32_t reserved;
-    int32_t count;
-    
-    masked_area_t area[1];
-} masked_area_list_t;
-
-typedef struct {
-    int32_t seqcount;
-    int32_t total;
-    
-    masked_area_list_t *sequence[1];
-} masked_area_table_t;
-
-masked_area_table_t *KMRK_ReadMaskArea(char *areafile,int32_t seqcount,int32_t complement);
-char KMRK_isMasked(masked_area_table_t *mask,int32_t seq, int32_t position);
-
-#endif

src/libKMRK/KMRK_merge_seeds.c

@@ -1,123 +0,0 @@
-/**
- * @file   KMRK_merge_seeds.c
- * @author Eric Coissac <coissac@inrialpes.fr>
- * @date   Wed Mar  3 11:15:57 2004
- * 
- * @brief  Merge function of overlapping seeds
- * 
- * 
- */
-
-#include "KMRK_merge_seeds.h"
-
-void KMRK_MergeSeeds(AllSeeds_type *AllSeeds, 
-		     int8_t opt_dir, 
-		     int8_t opt_inv)
-{
-
-  int32_t i;   /* the current seed */
-  int32_t j;   /* the checked seed */
-		 
-  int32_t N;   /* the kept ones */
-  Seed_type* seeds;
-	
-	
-  if(opt_dir){
-    seeds = AllSeeds->dirSeeds;
-    for(i=0, N=0 ;i<AllSeeds->nDirSeeds; i++){
-			
-      if(seeds[i].pos1==-1)                      /* any seed at -1 is removed */
-	continue;
-
-			
-      j=i+1;
-			
-      while( (seeds[j].pos1!=-1) && 
-	     (seeds[i].pos1!=-1) && 
-	     (j < AllSeeds->nDirSeeds)        &&
-	     (seeds[j].pos1 < seeds[i].pos1+ seeds[i].length))
-      {
-			
-	if( 
-	   ((seeds[j].pos2 >= seeds[i].pos2)  &&
-	    (seeds[j].pos2 <  seeds[i].pos2 + seeds[i].length)) ||        /* if the seeds are overlapping */
-	   ((seeds[j].pos2  + seeds[j].length >= seeds[i].pos2) &&
-	    (seeds[j].pos2  + seeds[j].length < seeds[i].pos2 + seeds[i].length))) 
-	  {										
-	    if(seeds[j].length <= seeds[i].length)
-	      seeds[j].pos1=seeds[j].pos2=seeds[j].length=-1;   /* removed the smallest */
-	    else
-	      seeds[i].pos1=seeds[i].pos2=seeds[i].length=-1;
-	  }
-				  
-	j++;
-      }
-
-      if(seeds[i].pos1 !=-1)
-      {                             /* if this seed is not out, store it */
-			
-	seeds[N].pos1   = seeds[i].pos1;
-	seeds[N].pos2   = seeds[i].pos2;
-	seeds[N].length = seeds[i].length;
-	N++;
-      }
-			
-    }
-		
-    AllSeeds->nFilteredDirSeeds += AllSeeds->nDirSeeds-N;		
-    AllSeeds->nDirSeeds=N;
-				
-			
-  }
-
-  if(opt_inv){
-    seeds = AllSeeds->invSeeds;
-	
-    for(i=0, N=0 ;i<AllSeeds->nInvSeeds; i++){
-			
-      if(seeds[i].pos1==-1)
-	continue;
-
-			
-      j=i+1;
-      
-      while( (seeds[j].pos1!=-1 ) && 
-	     (seeds[i].pos1!=-1 ) && 
-	     (j < AllSeeds->nInvSeeds) &&
-	     (seeds[j].pos1 < seeds[i].pos1+seeds[i].length))
-	{
-			
-	  if( 
-	     ((seeds[j].pos2 >= seeds[i].pos2)  &&                             /* if the seeds are overlapping */
-	      (seeds[j].pos2 <  seeds[i].pos2+seeds[i].length)) ||
-	     ((seeds[j].pos2 + seeds[j].length >= seeds[i].pos2) &&
-	      (seeds[j].pos2 + seeds[j].length <  seeds[i].pos2+seeds[i].length))) 
-	    {
-	      
-	      if(seeds[j].length <= seeds[i].length)
-		seeds[j].pos1=seeds[j].pos2=seeds[j].length=-1;   /* removed the smallest */
-	      else
-		seeds[i].pos1=seeds[i].pos2=seeds[i].length=-1;
-	    }
-				  
-	  j++;
-	}
-
-      if(seeds[i].pos1!=-1)
-	{                             /* if this seed is not out, store it */
-			
-	  seeds[N].pos1   = seeds[i].pos1;
-	  seeds[N].pos2   = seeds[i].pos2;
-	  seeds[N].length = seeds[i].length;
-	  N++;
-	}
-			
-		}
-		
-    AllSeeds->nFilteredInvSeeds += AllSeeds->nInvSeeds-N;		
-    AllSeeds->nInvSeeds=N;
-  }	
-
-  KMRK_compactSeeds(AllSeeds);
-
-}

src/libKMRK/KMRK_merge_seeds.h

@@ -1,11 +0,0 @@
-#ifndef KMRK_merge_seeds_h
-#define KMRK_merge_seeds_h
-
-#include "KMRK_Seeds.h"
-
-void KMRK_MergeSeeds(AllSeeds_type *AllSeeds, 
-		     int8_t opt_dir, 
-		     int8_t opt_inv);
-
-
-#endif /* KMRK_MergeSeeds */

src/libKMRK/memory.c

@@ -1,224 +0,0 @@
-/******
-        file     : memory.c
-        function : All about memory of the KMR, Seeds and Repeats
-		             All MyMalloc() series is about follwoing how mauch memory has been used
-		          
-        created  : 19 Sep 2003
-        modif    : Oct 2003, Feb 2004
-        modif    : Dec 2004 <EC> ; Corrected Memory declaration
-        author   : amikezor
-*****/
-
-#include <stdio.h>
-#include <stdlib.h>
-#include "repseek_types.h"
-#include "memory.h"
-
-
-MemUsage Memory;
-
-
-/*
-	Functions to count the memory usage all along
-	dybamic allocation and free
-*/
-void PrintMem(char *Comment){
-
-	extern MemUsage Memory;
-
-	fprintf(stderr,"\n%s\nMemory Usage\n\t* Max is: %d bytes, %.2f Kb, %.2f Mb\n\t* Cur is: %d bytes, %.2f Kb, %.2f Mb\n",
-	Comment,
-	Memory.Max, (float)Memory.Max/1024, (float)Memory.Max/(1024*1024),
-	Memory.Current, (float)Memory.Current/1024, (float)Memory.Current/(1024*1024));
-}
-void PrintMaxMem( void ){
-
-	extern MemUsage Memory;
-
-	if(Memory.Max < 1024)
-		fprintf(stderr,"Max Memory Usage.. %d bytes\n", Memory.Max);
-	else if(Memory.Max < 1024*1024)
-		fprintf(stderr,"Max Memory Usage.. %.2f Kilobytes\n", (float)Memory.Max/1024);
-	else if(Memory.Max < 1024*1024*1024)
-		fprintf(stderr,"Max Memory Usage.. %.2f Megabytes\n", (float)Memory.Max/(1024*1024));
-	else
-		fprintf(stderr,"Max Memory Usage.. %.2f Gigabytes\n", (float)Memory.Max/(1024*1024*1024));
-}
-void Update_Mem(int32_t Value){
-
-	extern MemUsage Memory;
-
-	Memory.Current += Value;
-	Memory.Max = (Memory.Current>Memory.Max)?Memory.Current:Memory.Max;
-}
-void Init_Mem(int32_t Value){
-
-	extern MemUsage Memory;
-
-	Memory.Current = Value;
-	Memory.Max = Value;
-}
-
-
-/*
-	Replace functions of dynamic allocation 
-	to allow the tracking of memory usage
-*/
-void *MyMalloc( int32_t size , char *Error ){
-
-	void *pointer;
-
-	pointer = malloc(size);
-	if(!pointer)fprintf(stderr,"%s\n",Error),exit(3);
-
-	Update_Mem(size);
-	
-	return pointer;
-}
-void *MyCalloc( int32_t number, int32_t TypeSize , char *Error ){
-
-	void *pointer;
-
-	pointer = calloc(number, TypeSize);
-	if(!pointer)fprintf(stderr,"%s\n",Error),exit(3);
-
-	Update_Mem(number*TypeSize );
-	
-	return pointer;
-}
-void MyFree( void *Pointer, int32_t size){
-	free(Pointer);
-	Pointer=NULL;
-	Update_Mem(-size);
-}
-void *MyRealloc( void *Pointer, int32_t newsize, int32_t oldsize, char *Error){
-
-	Pointer = realloc(Pointer,newsize);
-	if(!Pointer)fprintf(stderr,"%s\n",Error),exit(3);
-	Update_Mem( newsize-oldsize );
-	return Pointer;
-}
-
-
-/*
-	Deal with Stacks structure for KMR
-
-void MallocStack(Stacks *s, int32_t number, int32_t *histo, int32_t AllValues){
-
-	int32_t i;
-		
-	s->nStacks = number;
-	s->nValues =  AllValues;
-		
-	s->ppStacks = (int32_t **)MyMalloc( number * sizeof(int32_t *), 
-	                                  "MallocStack: ppStacks  malloc error, bye\n");
-	s->lenStacks = (int32_t *)MyMalloc( number * sizeof(int32_t),
-	                                "MallocStack: lenStacks malloc error, bye\n");
-	s->ppStacks[0] = (int32_t *)MyMalloc( AllValues * sizeof(int32_t),
-	                                  "MallocStack: ppStacks[0] malloc error, bye\n");
-	s->lenStacks[0]=0;
-	
-	for(i=1;i < number; i++){
-		s->lenStacks[i]=0;
-		s->ppStacks[i] = s->ppStacks[i-1] + histo[i] ;	
-	}
-}
-
-void FreeStack( Stacks *p){
-	MyFree(p->ppStacks[0]  , p->nValues*sizeof(int32_t) );
-	MyFree(p->ppStacks  ,    p->nStacks*sizeof(int32_t *) );
-	MyFree(p->lenStacks ,    p->nStacks*sizeof(int32_t));
-}
-*/
-
-
-/*
-	Deal with the Seeds part
-
-
-void free_Seeds(Seeds AllSeeds)    
-{
-
-	if( AllSeeds.nDirSeeds ){
-		MyFree(AllSeeds.DirPos1, AllSeeds.nDirSeeds*sizeof(int32_t));
-		MyFree(AllSeeds.DirPos2, AllSeeds.nDirSeeds*sizeof(int32_t));
-		MyFree(AllSeeds.DirLen, AllSeeds.nDirSeeds*sizeof(int32_t));
-		MyFree(AllSeeds.DirMeanR, AllSeeds.nDirSeeds*sizeof(float));
-	}
-	
-	if(AllSeeds.nInvSeeds ){
-		MyFree(AllSeeds.InvPos1, AllSeeds.nInvSeeds*sizeof(int32_t)); 
-		MyFree(AllSeeds.InvPos2, AllSeeds.nInvSeeds*sizeof(int32_t));
-		MyFree(AllSeeds.InvLen, AllSeeds.nInvSeeds*sizeof(int32_t));
-		MyFree(AllSeeds.InvMeanR, AllSeeds.nInvSeeds*sizeof(float));
-	}
-
-}
-*/
-/*
-	Malloc if it is the first time otherwise readjust
-
-void AllocSeeds(Seeds *AllSeeds, int32_t size, int32_t old_size, int8_t opt_dir, int8_t opt_inv){
-
-	if(opt_inv != 1 && opt_dir != 1)
-		fprintf(stderr,"AllocSeeds: requiere at least one of opt_dir and opt_inv to be 1\n"),exit(4);
-		
-	if(opt_dir == 1){
-		if( AllSeeds->DirPos1 == NULL){
-			AllSeeds->DirPos1 = (int32_t *)MyCalloc(size , sizeof(int32_t),"AllocSeeds: Alloc for DirPos1 failed, bye");
-			AllSeeds->DirPos2 = (int32_t *)MyCalloc(size , sizeof(int32_t),"AllocSeeds: Alloc for DirPos2 failed, bye");
-		}
-		else{
-			AllSeeds->DirPos1 = (int32_t *)MyRealloc(AllSeeds->DirPos1, size * sizeof(int32_t), old_size* sizeof(int32_t),
-			                                       "AllocSeeds: realloc for DirPos1 failed, bye");
-			AllSeeds->DirPos2 = (int32_t *)MyRealloc(AllSeeds->DirPos2, size * sizeof(int32_t), old_size* sizeof(int32_t),
-			                                       "AllocSeeds: realloc for DirPos2 failed, bye");
-		}
-	}
-
-	if(opt_inv == 1){
-		if( AllSeeds->InvPos1 == NULL){
-			AllSeeds->InvPos1 = (int32_t *)MyCalloc(size , sizeof(int32_t), "AllocSeeds: Alloc for InvPos1 failed, bye");
-			AllSeeds->InvPos2 = (int32_t *)MyCalloc(size , sizeof(int32_t), "AllocSeeds: Alloc for InvPos2 failed, bye");
-		}
-		else{
-			AllSeeds->InvPos1 = (int32_t *)MyRealloc(AllSeeds->InvPos1, size * sizeof(int32_t), old_size* sizeof(int32_t),
-			                                       "AllocSeeds: realloc for InvPos1 failed, bye");
-			AllSeeds->InvPos2 = (int32_t *)MyRealloc(AllSeeds->InvPos2, size * sizeof(int32_t), old_size* sizeof(int32_t),
-			                                       "AllocSeeds: realloc for InvPos2 failed, bye");
-		}
-	}
-}
-*/
-
-/*
-	Deal with the Repeats structure
-*/
-Repeats mem_Repeats(int32_t Ndir, int32_t Ninv){
-
-	Repeats AllRepeats;             /* All Repeats structure */
-	
-	AllRepeats.nDirRep = Ndir;      /* set the number of repeats to the number of seeds */
-	AllRepeats.nInvRep = Ninv;      
-	AllRepeats.nDirBadRep = 0;      /* set the "bad" repet (included into another rep) as 0 */
-	AllRepeats.nInvBadRep = 0;
-	                                
-	if(AllRepeats.nDirRep)
-		AllRepeats.DirRep = (Rep *)MyMalloc( (AllRepeats.nDirRep)*sizeof(Rep), "init_Repeats: repdir malloc error" );
-	else
-		AllRepeats.DirRep = NULL;
-
-	if(AllRepeats.nInvRep)
-		AllRepeats.InvRep = (Rep *)MyMalloc( (AllRepeats.nInvRep)*sizeof(Rep), "init_Repeats: repinv malloc error" );
-	else
-		AllRepeats.InvRep = NULL;
-
-	return AllRepeats;
-}
-void free_Repeats(Repeats AllRep)
-{
-	if(AllRep.nDirRep)
-		MyFree(AllRep.DirRep, AllRep.nDirRep*sizeof(Rep));
-	if(AllRep.nInvRep)
-		MyFree(AllRep.InvRep, AllRep.nInvRep*sizeof(Rep));
-}

src/libKMRK/memory.h

@@ -1,105 +0,0 @@
-/**
- * @file   memory.h
- * @author Achaz G
- * @date   April 2004
- * 
- * @brief  header for memory alloc/dealloc
- * modif    : Dec 2004 <EC> ; Corrected Memory declaration
- * 
- * 
- */
-
-#ifndef _MEMORY_H_
-#define _MEMORY_H_
-
-#include "repseek_types.h"
-
-typedef struct  {           /********** Memory Usage structure **************/
-
-	int32_t Max;
-	int32_t Current;
-
-} MemUsage;
-
-#include <stdint.h>
-
-/********** **********
-
-      Global Variable(s)
-		
- ********** **********/
-
-extern MemUsage Memory;  /* Instance of the global variable for memory tracking */
-
-/*
-	Follow memory usage
-*/
-void PrintMem(char *Comment);
-void PrintMaxMem( void );
-void Update_Mem(int32_t Value);
-void Init_Mem(int32_t Value);
-/*
-	All Alloc/Free to follow of memory usage
-*/
-void *MyMalloc( int32_t size , char *Error );
-void *MyCalloc( int32_t number, int32_t TypeSize , char *Error );
-void MyFree( void *Pointer, int32_t size);
-void *MyRealloc( void *Pointer, int32_t newsize, int32_t oldsize, char *Error);
-/*
-	For Stacks
-
-void MallocStack(Stacks *s, int32_t number, int32_t *histo, int32_t AllValues);
-void FreeStack( Stacks *p);
-
-	For Seeds
-
-void free_Seeds(Seeds AllSeeds);
-void AllocSeeds(Seeds *AllSeeds, int32_t size, int32_t old_size, int8_t opt_dir, int8_t opt_inv);
-*/
-
-/*
-	For Repeats
-*/
-Repeats mem_Repeats(int32_t Ndir, int32_t Ninv);
-void free_Repeats(Repeats AllRep);
-
-
-
-
-/*
-	Not used anymore, but just in case
-*/
-
-
-#include <stdio.h>
-#include <stdlib.h>
-
-#define KMRK_MALLOC(var,type,size,message) {       \
-            var = (type*) malloc(size);            \
-            if (!var)                              \
-              {                                    \
-		fprintf(stderr,"%s\n",message);    \
-		exit(4);                           \
-	      }                                    \
-            }
-
-#define KMRK_CALLOC(var,type,length,message) {             \
-            var = (type*) calloc(length,sizeof(type));     \
-            if (!var)                                      \
-              {                                            \
-		fprintf(stderr,"%s\n",message);            \
-		exit(4);                                   \
-	      }                                            \
-            }
-                
-#define KMRK_REALLOC(var,type,size,message) {              \
-            var = (type*) realloc(var,size);               \
-            if (!var)                                      \
-              {                                            \
-		fprintf(stderr,"%s\n",message);            \
-		exit(4);                                   \
-	      }                                            \
-            }
-
-
-#endif  /* _MEMORY_H_*/

src/libKMRK/repseek_types.h

@@ -1,197 +0,0 @@
-/**
- * @file   repseek_types.h
- * @author Guillaume Achaz <gachaz@oeb.harvard.edu>
- * @date   April 2004
- * @modif  July 2004 turn scores into doubles
- * @brief  definition of general types and macros for repseek
- * 
- * 
- */
-
-
-
-#ifndef _REPSEEK_TYPES_
-#define  _REPSEEK_TYPES_
-
-/*
-	Version of the program
-*/
-#define REPSEEK_VERSION "4.2"
-#define REPSEEK_DATE "Nov 2004"
-
-
-
-/********** **********
-
-        General Macros
-
- ********** **********/
-
-/*
-	Macros to compare 2 or three values
-*/	
-#define MAX2( A, B )     ( ((A)>(B))?(A):(B) )
-#define MAX3( A, B, C )  ( ((MAX2(A,B))>(C))?(MAX2(A,B)):(C) )
-#define MIN2( A, B )     ( ((A)<(B))?(A):(B) )
-
-#define MAX( A, B )     MAX2( A, B )
-#define MIN( A, B )     MIN2( A, B )
-
-/*
-	Absolute values
-*/
-#define ABS(x) (( (x)>=0 )? (x) : -(x))
-
-
-
-/********** **********
-
-   All types used in repseek
-
- ********** **********/
-
-#include <stdio.h>     /* The type FILE * is defined here */
-#include <stdint.h>    /* all, the int??_t are defined in there */
-
-
-/**
- * Store informations about one STRICT repeat (seeds)
- *
- */
-typedef struct {      /* the complete seed structure */
-
-  int32_t pos1;       /**< position of the first copy */
-  int32_t pos2;       /**< position of the second copy */
-
-  int32_t length;     /**< length of the strict repeats */
-  float   rmean;      /**< mean repeat leavel */
-  
-} Seed_type;
-
-typedef struct {      /* Just after a KMRK length X, only the 2 pos matter */
-
-  int32_t pos1;       /**< postion of the first copy */
-  int32_t pos2;       /**< postion of the second copy */
-
-} SmallSeed_type;   
-
-
-
-
-/**
- * Store informations about all strict repeat (seeds)
- *
- */
-typedef struct {
-
-  int32_t      cDirSeeds;          /**< currently allocated space in dirSeeds array */
-  int32_t      nDirSeeds;          /**< count of direct strict repeats */
-  int32_t      nFilteredDirSeeds;  /**< ??? */
-
-  Seed_type*   dirSeeds;           /**< array of direct repeats */
-
-  int32_t      cInvSeeds;          /**< currently allocated space in invSeeds array */
-  int32_t      nInvSeeds;          /**< count of inverted strict repeats */
-  int32_t      nFilteredInvSeeds;  /**< ??? */
-
-  Seed_type*   invSeeds;           /**< array of inverted repeats */
-
-} AllSeeds_type;
-
-
-
-
-
-/**
- * Store informations about one GENERIC repeat
- *
- */
-typedef struct{              
-
-  char type[20];               /* its name; i.e. Tandem, Palindrome, etc... */
-  
-  int32_t pos1, pos2,          /* both copies postions */
-        len1, len2,            /* both copies length */
-        seed_pos1,seed_pos2,   /* pos1 and pos2 of the originate seed */
-        seed_len,              /* len  of  the seed	*/
-		  
-  	match, align_len;           /* number of match and length of alignment */
-	
-	double score;		          /* the alignment score */
-
-	float seed_meanR;           /* the seed meanR */
-
-	float meanR;                /* The mean R-level of the repeat */
-	int32_t mainR;              /* its Mode R */
-	float fraction_mainR;       /* the fraction of length containing the Mode R */
-
-
-} Rep;
-
-
-
-/**
- * Store informations about All GENERIC repeats
- *
- */
-typedef struct {  
-			  
-	int32_t nDirRep;	        /* Total Number of Direct Repats in Mem */
-	int32_t nDirBadRep;	     /* Direct repeats set to -1 -- filtered out and co. */
-	Rep *DirRep;           /* The array of Dir Rep */
-
-	int32_t nInvRep;	        /* Total Number of Inverted Repats in Mem */
-	int32_t nInvBadRep;		  /* Inverted Repeats set to -1 -- filtered out and co. */
-	Rep *InvRep;			  /* The array of Inverted Rep */
-
-} Repeats;
-
-
-#define MATRIX_SIZE 26
-
-typedef struct {            /******* The scoring Matrix  ************/
-
-	double matrix[MATRIX_SIZE][MATRIX_SIZE];  /* the matrix of match/missmatch */
-	double gap_open;                          /* value of gap-open */
-	double gap_ext;                           /* value of gap_ext */
-	double expect;
-
-} SCORING;
-
-
-
-typedef struct {            /******* The Results of Alignement by Dynamik programming  ************/
-
-
-	double *scores;                 /* the score strings (+/- 'matrix') */
-	double *pscore;                 /* pointer to the current score */
-	char *traces;                   /* the path matrix - could take values 'l'eft, 'd'iagonal, or 'a'bove or 'u'nknown */
-
-	double *F;                      /* *Above -- needed for memorizing deletion in seq2 */
-	
-	double *pBestScore;             /* pointer to it */
-	int32_t BestScore_row;          /* its row and col */
-	int32_t BestScore_col;   
-
-	char *traceback;                /* all you need for bactracking */
-	char *traceback_f;              /* memory for forward traceback and then check other seeds */
-	char *traceback_b;              /* memory needed for backward traceback - to avoid erasing the forward one */
-
-	int32_t alignment_len;          /* guess ?? */
-	int32_t matches;                /* number of match (score>0 in scoring matrix) */
-
-	int32_t nSegment;               /* number of segment */
-	int32_t *Segment_begin;         /* begin and end of each segment */
-	int32_t *Segment_end;
-   
-	int32_t max_scores;             /* size of the matrices only for memory purposes */
-	int32_t max_col;
-	int32_t max_row;
-	int32_t max_alignlen;
-
-} RESULTS;
-
-
-
-#endif  /* _REPSEEK_TYPES_ */
-

src/libKMRK/sequence.h

@@ -1,25 +0,0 @@
-/**
- * @file   KMRK_sequence.h
- * @author Eric Coissac <coissac@inrialpes.fr>
- * @date   Tue Feb 24 22:22:57 2004
- * 
- * @brief  Header file for sequence utilities
- * 
- * 
- */
-
-#ifndef KMRK_sequence_h
-#define KMRK_sequence_h
-
-#include <stdint.h>
-
-int8_t CheckSeq(char *seq, char *alpha);
-
-void nonACGTXtoN(char *seq);
-
-void UpperSequence(char *seq);
-
-void invseq(char *seqsrc, char *seqdest);
-
-
-#endif /* KMRK_sequence_h */

src/libapat/CODES/dft_code.h

@@ -1,14 +0,0 @@
-/* ----------------------------------------------- */
-/* dft_pat_seq_code.h                              */
-/* default alphabet encoding for alpha             */
-/* ----------------------------------------------- */
-
-        0x00000001 /* A */, 0x00000002 /* B */, 0x00000004 /* C */,
-        0x00000008 /* D */, 0x00000010 /* E */, 0x00000020 /* F */,
-        0x00000040 /* G */, 0x00000080 /* H */, 0x00000100 /* I */,
-        0x00000200 /* J */, 0x00000400 /* K */, 0x00000800 /* L */,
-        0x00001000 /* M */, 0x00002000 /* N */, 0x00004000 /* O */,
-        0x00008000 /* P */, 0x00010000 /* Q */, 0x00020000 /* R */,
-        0x00040000 /* S */, 0x00080000 /* T */, 0x00100000 /* U */,
-        0x00200000 /* V */, 0x00400000 /* W */, 0x00800000 /* X */,
-        0x01000000 /* Y */, 0x02000000 /* Z */

src/libapat/CODES/dna_code.h

@@ -1,71 +0,0 @@
-/* ----------------------------------------------- */
-/* dna_code.h                                      */
-/* alphabet encoding for dna/rna                   */
-/* -----------------------------------------       */
-/* IUPAC encoding                                  */
-/* -----------------------------------------       */
-/* G/A/T/C                                         */
-/* U=T                                             */
-/* R=AG                                            */
-/* Y=CT                                            */
-/* M=AC                                            */
-/* K=GT                                            */
-/* S=CG                                            */
-/* W=AT                                            */
-/* H=ACT                                           */
-/* B=CGT                                           */
-/* V=ACG                                           */
-/* D=AGT                                           */
-/* N=ACGT                                          */
-/* X=ACGT                                          */
-/* EFIJLOPQZ  not recognized                       */
-/* -----------------------------------------       */
-/* dual encoding                                   */
-/* -----------------------------------------       */
-/* A=ADHMNRVW                                      */
-/* B=BCDGHKMNRSTUVWY                               */
-/* C=BCHMNSVY                                      */
-/* D=ABDGHKMNRSTUVWY                               */
-/* G=BDGKNRSV                                      */
-/* H=ABCDHKMNRSTUVWY                               */
-/* K=BDGHKNRSTUVWY                                 */
-/* M=ABCDHMNRSVWY                                  */
-/* N=ABCDGHKMNRSTUVWY                              */
-/* R=ABDGHKMNRSVW                                  */
-/* S=BCDGHKMNRSVY                                  */
-/* T=BDHKNTUWY                                     */
-/* U=BDHKNTUWY                                     */
-/* V=ABCDGHKMNRSVWY                                */
-/* W=ABDHKMNRTUVWY                                 */
-/* X=ABCDGHKMNRSTUVWY                              */
-/* Y=BCDHKMNSTUVWY                                 */
-/* EFIJLOPQZ  not recognized                       */
-/* ----------------------------------------------- */
-
-#ifndef USE_DUAL
-
-        /* IUPAC */
-
-        0x00000001 /* A */, 0x00080044 /* B */, 0x00000004 /* C */, 
-        0x00080041 /* D */, 0x00000000 /* E */, 0x00000000 /* F */, 
-        0x00000040 /* G */, 0x00080005 /* H */, 0x00000000 /* I */, 
-        0x00000000 /* J */, 0x00080040 /* K */, 0x00000000 /* L */, 
-        0x00000005 /* M */, 0x00080045 /* N */, 0x00000000 /* O */, 
-        0x00000000 /* P */, 0x00000000 /* Q */, 0x00000041 /* R */, 
-        0x00000044 /* S */, 0x00080000 /* T */, 0x00080000 /* U */, 
-        0x00000045 /* V */, 0x00080001 /* W */, 0x00080045 /* X */, 
-        0x00080004 /* Y */, 0x00000000 /* Z */
-
-#else
-        /* DUAL  */
-
-        0x00623089 /* A */, 0x017e34ce /* B */, 0x01243086 /* C */, 
-        0x017e34cb /* D */, 0x00000000 /* E */, 0x00000000 /* F */, 
-        0x0026244a /* G */, 0x017e348f /* H */, 0x00000000 /* I */, 
-        0x00000000 /* J */, 0x017e24ca /* K */, 0x00000000 /* L */, 
-        0x0166308f /* M */, 0x017e34cf /* N */, 0x00000000 /* O */, 
-        0x00000000 /* P */, 0x00000000 /* Q */, 0x006634cb /* R */, 
-        0x012634ce /* S */, 0x0158248a /* T */, 0x0158248a /* U */, 
-        0x016634cf /* V */, 0x017a348b /* W */, 0x017e34cf /* X */, 
-        0x017c348e /* Y */, 0x00000000 /* Z */
-#endif

src/libapat/CODES/prot_code.h

@@ -1,51 +0,0 @@
-/* ----------------------------------------------- */
-/* prot_code.h                                     */
-/* alphabet encoding for proteins                  */
-/* -----------------------------------------       */
-/* IUPAC encoding                                  */
-/* -----------------------------------------       */
-/* B=DN                                            */
-/* Z=EQ                                            */
-/* X=any - {X}                                     */
-/* JOU not recognized                              */
-/* -----------------------------------------       */
-/* dual encoding                                   */
-/* -----------------------------------------       */
-/* B=BDN                                           */
-/* D=BD                                            */
-/* E=EZ                                            */
-/* N=BN                                            */
-/* Q=QZ                                            */
-/* X=any - {X}                                     */
-/* Z=EQZ                                           */
-/* JOU not recognized                              */
-/* ----------------------------------------------- */
-
-#ifndef USE_DUAL
-
-        /* IUPAC */
-
-        0x00000001 /* A */, 0x00002008 /* B */, 0x00000004 /* C */,
-        0x00000008 /* D */, 0x00000010 /* E */, 0x00000020 /* F */,
-        0x00000040 /* G */, 0x00000080 /* H */, 0x00000100 /* I */,
-        0x00000000 /* J */, 0x00000400 /* K */, 0x00000800 /* L */,
-        0x00001000 /* M */, 0x00002000 /* N */, 0x00000000 /* O */,
-        0x00008000 /* P */, 0x00010000 /* Q */, 0x00020000 /* R */,
-        0x00040000 /* S */, 0x00080000 /* T */, 0x00000000 /* U */,
-        0x00200000 /* V */, 0x00400000 /* W */, 0x037fffff /* X */,
-        0x01000000 /* Y */, 0x00010010 /* Z */
-
-#else
-        /* DUAL  */
-
-        0x00000001 /* A */, 0x0000200a /* B */, 0x00000004 /* C */, 
-        0x0000000a /* D */, 0x02000010 /* E */, 0x00000020 /* F */, 
-        0x00000040 /* G */, 0x00000080 /* H */, 0x00000100 /* I */, 
-        0x00000000 /* J */, 0x00000400 /* K */, 0x00000800 /* L */, 
-        0x00001000 /* M */, 0x00002002 /* N */, 0x00000000 /* O */, 
-        0x00008000 /* P */, 0x02010000 /* Q */, 0x00020000 /* R */, 
-        0x00040000 /* S */, 0x00080000 /* T */, 0x00000000 /* U */, 
-        0x00200000 /* V */, 0x00400000 /* W */, 0x037fffff /* X */, 
-        0x01000000 /* Y */, 0x02010010 /* Z */
-
-#endif

src/libapat/Gmach.h

@@ -1,97 +0,0 @@
-/* ---------------------------------------------------------------- */
-/* Copyright (c) Atelier de BioInformatique                         */
-/* @file: Gmach.h                                                   */
-/* @desc: machine dependant setup                                   */
-/* @+     *should* be included in all ABI softs                     */
-/*                                                                  */
-/* @history:                                                        */
-/* @+       <Gloup> : Jul 95 : MWC first draft                      */
-/* @+       <Gloup> : Jan 96 : adapted to Pwg                       */
-/* @+       <Gloup> : Nov 00 : adapted to Mac_OS_X                  */
-/* ---------------------------------------------------------------- */
-
-#ifndef _H_Gmach
-
-                        /* OS names             */
-
-#define _H_Gmach
-
-                        /* Macintosh Classic    */
-                        /* Think C environment  */
-#ifdef THINK_C 
-#define MACINTOSH
-#define MAC_OS_C
-#endif
-
-
-                        /* Macintosh Classic    */
-                        /* Code-Warrior         */
-#ifdef __MWERKS__ 
-#define MACINTOSH
-#define MAC_OS_C
-#endif
-
-                        /* Macintosh OS-X       */
-#ifdef MAC_OS_X
-#define MACINTOSH
-#define UNIX
-#define UNIX_BSD
-#endif
-
-						/* LINUX				*/
-#ifdef LINUX
-#define UNIX
-#define UNIX_BSD
-#endif
-
-                        /* other Unix Boxes     */
-                        /* SunOS /  Solaris     */      
-#ifdef SUN
-#define UNIX
-#ifdef SOLARIS
-#define UNIX_S7
-#else
-#define UNIX_BSD
-#endif
-#endif
-
-                        /* SGI Irix             */      
-#ifdef SGI
-#define UNIX
-#define UNIX_S7
-#endif
-
-/* ansi setup                           */
-/* for unix machines see makefile       */
-
-#ifndef PROTO
-#define PROTO 1
-#endif
-
-#ifndef ANSI_PROTO
-#define ANSI_PROTO PROTO
-#endif
-
-#ifndef ANSI_STR
-#define ANSI_STR 1
-#endif
-
-/* unistd.h header file */
-
-#ifdef UNIX
-#define HAS_UNISTD_H <unistd.h>
-#endif
-
-/* getopt.h header file */
-
-#ifdef MAC_OS_C
-#define HAS_GETOPT_H "getopt.h"
-#endif
-
-#ifdef SGI
-#define HAS_GETOPT_H <getopt.h>
-#endif
-
-
-
-#endif

src/libapat/Gtypes.h

@@ -1,104 +0,0 @@
-/* ---------------------------------------------------------------- */
-/* Copyright (c) Atelier de BioInformatique                         */
-/* @file: Gtypes.h                                                  */
-/* @desc: general & machine dependant types                         */
-/* @+     *should* be included in all ABI softs                     */
-/*                                                                  */
-/* @history:                                                        */
-/* @+       <Gloup> : Jan 91 : MWC first draft                      */
-/* @+       <Gloup> : Jul 95 : Gmach addition                       */
-/* ---------------------------------------------------------------- */
-
-#define _H_Gtypes
-
-#ifndef _H_Gmach
-#include "Gmach.h"
-#endif
-
-#ifndef NULL
-#include <stdio.h>                      /* is the official NULL here ?  */
-#endif
-
-/* ==================================================== */
-/* constantes                                           */
-/* ==================================================== */
-
-#ifndef PROTO
-#define PROTO   1                       /* prototypes flag              */
-#endif
-
-#ifdef MAC_OS_C
-#define Vrai    true                    /* TC boolean values            */
-#define Faux    false                   /*                              */
-#else
-#define Vrai    0x1                     /* bool values  = TRUE          */
-#define Faux    0x0                     /*              = FALSE         */
-#endif
-
-#define Nil     NULL                    /* nil pointer                  */
-
-#define kBigInt16       0x7fff          /* plus grand 16 bits signe     */
-#define kBigInt32       0x7fffffff      /* plus grand 32 bits signe     */
-#define kBigUInt16      0xffff          /* plus grand 16 bits ~signe    */
-#define kBigUInt32      0xffffffff      /* plus grand 32 bits ~signe    */
-
-#ifdef MAC_OS_C
-/* ==================================================== */
-/*  Types (for Macintosh ThinK C || MWerks)             */
-/* ==================================================== */
-
-                                        /* --- specific sizes --------- */
-typedef long            Int32;          /* Int32  = 32 bits signe       */
-typedef unsigned long   UInt32;         /* UInt32 = 32 bits ~signe      */
-typedef short           Int16;          /* Int16  = 16 bits signe       */
-typedef unsigned short  UInt16;         /* UInt32 = 16 bits ~signe      */
-typedef char            Int8;           /* Int8   = 8 bits signe        */
-typedef unsigned char   UInt8;          /* UInt8  = 8 bits ~signe       */
-
-                                        /* --- default types ---------- */
-
-typedef Boolean         Bool;           /* booleen                      */
-
-typedef long            Int;            /* 'natural' int (>= 32 bits)   */
-
-typedef void            *Ptr;           /* pointeur                     */
-
-#elif ((defined SUN) || (defined SGI) || (defined UNIX))
-/* ==================================================== */
-/*  Types (for Sun & Iris - 32 bits machines)           */
-/* ==================================================== */
-
-                                        /* --- specific sizes --------- */
-typedef int             Int32;          /* Int32  = 32 bits signe       */
-typedef unsigned int    UInt32;         /* UInt32 = 32 bits ~signe      */
-typedef short           Int16;          /* Int16  = 16 bits signe       */
-typedef unsigned short  UInt16;         /* UInt32 = 16 bits ~signe      */
-typedef char            Int8;           /* Int8   = 8 bits signe        */
-typedef unsigned char   UInt8;          /* UInt8  = 8 bits ~signe       */
-
-                                        /* --- default types ---------- */
-
-typedef int             Bool;           /* booleen  (int for ANSI)      */
-
-typedef int             Int;            /* 'natural' int (>= 32 bits)   */
-
-typedef void            *Ptr;           /* pointeur                     */
-
-#else 
-/* ==================================================== */
-/*  Types (for undefined machines)                      */
-/* ==================================================== */
-
-#error undefined MACHINE <please edit Gmach.h>
-
-#endif
-
-/* ==================================================== */
-/*  special macro for prototypes                        */
-/* ==================================================== */
-
-#if PROTO
-#define         P(s)    s 
-#else
-#define         P(s)    ()
-#endif

src/libapat/Makefile

@@ -1,24 +0,0 @@
-
-SOURCES = apat_parse.c \
-         apat_search.c \
-         libstki.c
-         
-SRCS=$(SOURCES)
-
-         
-OBJECTS= $(patsubst %.c,%.o,$(SOURCES))
-
-LIBFILE= libapat.a
-RANLIB=ranlib
-
-
-include ../global.mk
-
-all: $(LIBFILE)
-
-clean:
-	rm -rf $(OBJECTS) $(LIBFILE)
-
-$(LIBFILE): $(OBJECTS)
-	ar -cr $@ $?
-	$(RANLIB) $@

src/libapat/apat.h

@@ -1,173 +0,0 @@
-/* ==================================================== */
-/*      Copyright (c) Atelier de BioInformatique        */
-/*      Dec. 94                                         */
-/*      File: apat.h                                    */
-/*      Purpose: pattern scan                           */
-/*      History:                                        */
-/*      28/12/94 : <Gloup> ascan first version          */
-/*      14/05/99 : <Gloup> last revision                */
-/* ==================================================== */
-
-#ifndef _H_Gtypes
-#include "Gtypes.h"
-#endif
-
-#ifndef _H_libstki
-#include "libstki.h" 
-#endif
-
-#define H_apat
-
-/* ----------------------------------------------- */
-/* constantes                                      */
-/* ----------------------------------------------- */
-
-#ifndef BUFSIZ
-#define BUFSIZ          1024    /* io buffer size               */
-#endif
-
-#define MAX_NAME_LEN    BUFSIZ  /* max length of sequence name  */
-
-#define ALPHA_LEN        26     /* alphabet length              */
-                                /* *DO NOT* modify              */
-
-#define MAX_PATTERN       4     /* max # of patterns            */
-                                /* *DO NOT* modify              */
-
-#define MAX_PAT_LEN      32     /* max pattern length           */
-                                /* *DO NOT* modify              */
-
-#define MAX_PAT_ERR      32     /* max # of errors              */
-                                /* *DO NOT* modify              */
-
-#define PATMASK 0x3ffffff       /* mask for 26 symbols          */
-                                /* *DO NOT* modify              */
-
-#define OBLIBIT 0x4000000       /* bit 27 to 1 -> oblig. pos    */
-                                /* *DO NOT* modify              */
-
-                                /* mask for position            */
-#define ONEMASK 0x80000000      /* mask for highest position    */
-
-                                /* masks for Levenhstein edit   */
-#define OPER_IDT  0x00000000    /* identity                     */
-#define OPER_INS  0x40000000    /* insertion                    */
-#define OPER_DEL  0x80000000    /* deletion                     */
-#define OPER_SUB  0xc0000000    /* substitution                 */
-
-#define OPER_SHFT 30            /* <unused> shift               */
-
-                                /* Levenhstein Opcodes          */
-#define SOPER_IDT 0x0           /* identity                     */
-#define SOPER_INS 0x1           /* insertion                    */
-#define SOPER_DEL 0x2           /* deletion                     */
-#define SOPER_SUB 0x3           /* substitution                 */
-
-                                /* Levenhstein Opcodes masks    */
-#define OPERMASK  0xc0000000    /* mask for Opcodes             */
-#define NOPERMASK 0x3fffffff    /* negate of previous           */
-
-                                /* special chars in pattern     */
-#define PATCHARS  "[]!#"
-
-                                /* 26 letter alphabet           */
-                                /* in alphabetical order        */
-
-#define ORD_ALPHA "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
-
-                                /* protein alphabet             */
-
-#define PROT_ALPHA "ACDEFGHIKLMNPQRSTVWY"
-
-                                /* dna/rna alphabet             */
-
-#define DNA_ALPHA "ABCDGHKMNRSTUVWXY"
-
-
-/* ----------------------------------------------- */
-/* data structures                                 */
-/* ----------------------------------------------- */
-
-                                        /* -------------------- */
-typedef enum {                          /* data encoding        */
-                                        /* -------------------- */
-        alpha = 0,                      /* [A-Z]                */
-        dna,                            /* IUPAC DNA            */
-        protein                         /* IUPAC proteins       */
-} CodType;
-
-                                        /* -------------------- */
-typedef struct {                        /* sequence             */
-                                        /* -------------------- */
-    char      *name;                    /* sequence name        */
-    Int32     seqlen;                   /* sequence length      */
-    Int32     seqsiz;                   /* sequence buffer size */
-    Int32     datsiz;                   /* data buffer size     */
-    Int32     circular;
-    UInt8     *data;                    /* data buffer          */
-    char      *cseq;                    /* sequence buffer      */
-    StackiPtr hitpos[MAX_PATTERN];      /* stack of hit pos.    */
-    StackiPtr hiterr[MAX_PATTERN];      /* stack of errors      */
-} Seq, *SeqPtr;
-
-                                        /* -------------------- */
-typedef struct {                        /* pattern              */
-                                        /* -------------------- */
-    int    patlen;                      /* pattern length       */
-    int    maxerr;                      /* max # of errors      */
-    char   *cpat;                       /* pattern string       */
-    Int32  *patcode;                    /* encoded pattern      */
-    UInt32 *smat;                       /* S matrix             */
-    UInt32 omask;                       /* oblig. bits mask     */
-    Bool   hasIndel;                    /* are indels allowed   */
-    Bool   ok;                          /* is pattern ok        */
-} Pattern, *PatternPtr;
-
-/* ----------------------------------------------- */
-/* macros                                          */
-/* ----------------------------------------------- */
-
-#ifndef NEW
-#define NEW(typ)                (typ*)malloc(sizeof(typ)) 
-#define NEWN(typ, dim)          (typ*)malloc((unsigned long)(dim) * sizeof(typ))
-#define REALLOC(typ, ptr, dim)  (typ*)realloc((void *) (ptr), (unsigned long)(dim) * sizeof(typ))
-#define FREE(ptr)               free((void *) ptr)
-#endif
-
-/* ----------------------------------------------- */
-/* prototypes                                      */
-/* ----------------------------------------------- */
-
-                                     /* apat_seq.c */
-
-SeqPtr  FreeSequence     (SeqPtr pseq);
-SeqPtr  NewSequence      (void);
-int     ReadNextSequence (SeqPtr pseq);
-int     WriteSequence    (FILE *filou , SeqPtr pseq);
-
-                                   /* apat_parse.c      */
-
-Int32   *GetCode          (CodType ctype);
-int     CheckPattern      (Pattern *ppat);
-int     EncodePattern     (Pattern *ppat, CodType ctype);
-int     ReadPattern       (Pattern *ppat);
-void    PrintDebugPattern (Pattern *ppat);
-
-                                /* apat_search.c        */
-
-int     CreateS           (Pattern *ppat, Int32 lalpha);
-Int32   ManberNoErr       (Seq *pseq , Pattern *ppat, int patnum,int begin,int length);
-Int32   ManberSub         (Seq *pseq , Pattern *ppat, int patnum,int begin,int length);
-Int32   ManberIndel       (Seq *pseq , Pattern *ppat, int patnum,int begin,int length);
-Int32   ManberAll         (Seq *pseq , Pattern *ppat, int patnum,int begin,int length);
-Int32   NwsPatAlign       (Seq *pseq , Pattern *ppat, Int32 nerr , 
-                              Int32 *reslen , Int32 *reserr);
-
-                                   /* apat_sys.c   */
-
-float   UserCpuTime     (int reset);
-float   SysCpuTime      (int reset);
-char    *StrCpuTime     (int reset);
-void    Erreur          (char *msg , int stat);
-int     AccessFile      (char *path, char *mode);
-

src/libapat/apat_parse.c

@@ -1,369 +0,0 @@
-/* ==================================================== */
-/*      Copyright (c) Atelier de BioInformatique        */
-/*      Mar. 92                                         */
-/*      File: apat_parse.c                              */
-/*      Purpose: Codage du pattern                      */
-/*      History:                                        */
-/*      00/07/94 : <Gloup> first version (stanford)     */
-/*      00/11/94 : <Gloup> revised for DNA/PROTEIN      */
-/*      30/12/94 : <Gloup> modified EncodePattern       */
-/*                         for manber search            */
-/*      14/05/99 : <Gloup> indels added                 */
-/* ==================================================== */
-
-#include <ctype.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-
-#include "Gtypes.h"
-#include "apat.h"
-                                /* -------------------- */
-                                /* default char         */
-                                /* encodings            */
-                                /* -------------------- */
-
-static Int32 sDftCode[]  =  {
-
-#include "CODES/dft_code.h"
-
-};
-                                /* -------------------- */
-                                /* char encodings       */
-                                /* IUPAC                */
-                                /* -------------------- */
-
-                                /* IUPAC Proteins       */
-static Int32 sProtCode[]  =  {
-
-#include "CODES/prot_code.h"
-
-};
-                                /* IUPAC Dna/Rna        */
-static Int32 sDnaCode[]  =  {
-
-#include "CODES/dna_code.h"
-
-};
-
-
-/* -------------------------------------------- */
-/* internal replacement of gets                 */
-/* -------------------------------------------- */
-static char *sGets(char *buffer, int size) {
-        
-        char *ebuf;
-
-        if (! fgets(buffer, size-1, stdin))
-           return NULL;
-
-        /* remove trailing line feed */
-
-        ebuf = buffer + strlen(buffer); 
-
-        while (--ebuf >= buffer) {
-           if ((*ebuf == '\n') || (*ebuf == '\r'))
-                *ebuf = '\000';
-           else
-                break;
-        }
-
-        return buffer;
-}
-
-/* -------------------------------------------- */
-/* returns actual code associated to type       */
-/* -------------------------------------------- */
-
-Int32 *GetCode(CodType ctype)
-{
-        Int32 *code = sDftCode;
-
-        switch (ctype) {
-           case dna     : code = sDnaCode  ; break;
-           case protein : code = sProtCode ; break;
-           default      : code = sDftCode  ; break;
-        }
-
-        return code;
-}
-
-/* -------------------------------------------- */
-
-#define BAD_IF(tst)   if (tst)  return 0
-
-int CheckPattern(Pattern *ppat)
-{
-        int lev;
-        char *pat;
-
-        pat = ppat->cpat;
-        
-        BAD_IF (*pat == '#');
-
-        for (lev = 0; *pat ; pat++)
-
-            switch (*pat) {
-
-                case '[' :
-                   BAD_IF (lev);
-                   BAD_IF (*(pat+1) == ']');
-                   lev++;
-                   break;
-
-                case ']' :
-                   lev--;
-                   BAD_IF (lev);
-                   break;
-
-                case '!' :
-                   BAD_IF (lev);
-                   BAD_IF (! *(pat+1));
-                   BAD_IF (*(pat+1) == ']');
-                   break;
-
-                case '#' :
-                   BAD_IF (lev);
-                   BAD_IF (*(pat-1) == '[');
-                   break;
-
-                default :
-                   if (! isupper(*pat))
-                        return 0;
-                   break;
-            }
-
-        return (lev ? 0 : 1);
-}
- 
-#undef BAD_IF
-           
-
-/* -------------------------------------------- */
-static char *skipOblig(char *pat)
-{
-        return (*(pat+1) == '#' ? pat+1 : pat);
-}
-
-/* -------------------------------------------- */
-static char *splitPattern(char *pat)
-{
-        switch (*pat) {
-                   
-                case '[' :
-                   for (; *pat; pat++)
-                        if (*pat == ']')
-                          return skipOblig(pat);
-                   return NULL;
-                   break;
-
-                case '!' :
-                   return splitPattern(pat+1);
-                   break;
-
-        }
-        
-        return skipOblig(pat);                  
-}
-
-/* -------------------------------------------- */
-static Int32 valPattern(char *pat, Int32 *code)
-{
-        Int32 val;
-        
-        switch (*pat) {
-                   
-                case '[' :
-                   return valPattern(pat+1, code);
-                   break;
-
-                case '!' :
-                   val = valPattern(pat+1, code);
-                   return (~val & PATMASK);
-                   break;
-
-                default :
-                   val = 0x0;
-                   while (isupper(*pat)) {
-                       val |= code[*pat - 'A'];
-                       pat++;
-                   }
-                   return val;
-        }
-
-        return 0x0;                     
-}
-
-/* -------------------------------------------- */
-static Int32 obliBitPattern(char *pat)
-{
-        return (*(pat + strlen(pat) - 1) == '#' ? OBLIBIT : 0x0);
-}       
-                
-
-/* -------------------------------------------- */
-static int lenPattern(char *pat)
-{
-        int  lpat;
-
-        lpat = 0;
-        
-        while (*pat) {
-        
-            if (! (pat = splitPattern(pat)))
-                return 0;
-
-            pat++;
-
-            lpat++;
-        }
-
-        return lpat;
-}
-
-/* -------------------------------------------- */
-/* Interface                                    */
-/* -------------------------------------------- */
-
-/* -------------------------------------------- */
-/* encode un pattern                            */
-/* -------------------------------------------- */
-int EncodePattern(Pattern *ppat, CodType ctype)
-{
-        int   pos, lpat;
-        Int32 *code;
-        char  *pp, *pa, c;
-
-        ppat->ok = Faux;
-
-        code = GetCode(ctype);
-
-        ppat->patlen = lpat = lenPattern(ppat->cpat);
-        
-        if (lpat <= 0)
-            return 0;
-        
-        if (! (ppat->patcode = NEWN(Int32, lpat)))
-            return 0;
-
-        pa = pp = ppat->cpat;
-
-        pos = 0;
-        
-        while (*pa) {
-        
-            pp = splitPattern(pa);
-
-            c = *++pp;
-            
-            *pp = '\000';
-                    
-            ppat->patcode[pos++] = valPattern(pa, code) | obliBitPattern(pa);
-            
-            *pp = c;
-            
-            pa = pp;
-        }
-
-        ppat->ok = Vrai;
-
-        return lpat;
-}
-
-/* -------------------------------------------- */
-/* remove blanks                                */
-/* -------------------------------------------- */
-static char *RemBlanks(char *s)
-{
-        char *sb, *sc;
-
-        for (sb = sc = s ; *sb ; sb++)
-           if (! isspace(*sb))
-                *sc++ = *sb;
-
-        return s;
-}
-
-/* -------------------------------------------- */
-/* count non blanks                             */
-/* -------------------------------------------- */
-static Int32 CountAlpha(char *s)
-{
-        Int32 n;
-
-        for (n = 0 ; *s ; s++)
-           if (! isspace(*s))
-                n++;
-
-        return n;
-}
-           
-        
-/* -------------------------------------------- */
-/* lit un pattern                               */
-/* <pattern> #mis                               */
-/* ligne starting with '/' are comments         */
-/* -------------------------------------------- */
-int ReadPattern(Pattern *ppat)
-{
-        int  val;
-        char *spac;
-        char buffer[BUFSIZ];
-
-        ppat->ok = Vrai;
-
-        if (! sGets(buffer, sizeof(buffer)))
-            return 0;
-
-        if (*buffer == '/')
-            return ReadPattern(ppat);
-
-        if (! CountAlpha(buffer))
-            return ReadPattern(ppat);
-
-        for (spac = buffer ; *spac ; spac++)
-            if ((*spac == ' ') || (*spac == '\t'))
-                break;
-
-        ppat->ok = Faux;
-
-        if (! *spac)
-            return 0;
-        
-        if (sscanf(spac, "%d", &val) != 1)
-            return 0;
-
-        ppat->hasIndel = (val < 0);
-        
-        ppat->maxerr = ((val >= 0) ? val : -val);
-
-        *spac = '\000';
-
-        (void) RemBlanks(buffer);
-
-        if ((ppat->cpat = NEWN(char, strlen(buffer)+1)))
-            strcpy(ppat->cpat, buffer);
-        
-        ppat->ok = (ppat->cpat != NULL);
-
-        return (ppat->ok ? 1 : 0);
-}
-
-/* -------------------------------------------- */
-/* ecrit un pattern - Debug -                   */
-/* -------------------------------------------- */
-void PrintDebugPattern(Pattern *ppat)
-{
-        int i;
-
-        printf("Pattern  : %s\n", ppat->cpat);
-        printf("Encoding : \n\t");
-
-        for (i = 0 ; i < ppat->patlen ; i++) {
-            printf("0x%8.8x ", ppat->patcode[i]);
-            if (i%4 == 3)
-                printf("\n\t");
-        }
-        printf("\n");
-}
-

src/libapat/apat_search.c

@@ -1,339 +0,0 @@
-/* ==================================================== */
-/*      Copyright (c) Atelier de BioInformatique        */
-/*      Dec. 94                                         */
-/*      File: apat_search.c                             */
-/*      Purpose: recherche du pattern                   */
-/*               algorithme de Baeza-Yates/Gonnet       */
-/*                             Manber (agrep)           */
-/*      History:                                        */
-/*      07/12/94 : <MFS>   first version                */
-/*      28/12/94 : <Gloup> revised version              */
-/*      14/05/99 : <Gloup> last revision                */
-/* ==================================================== */
-
-#if 0
-#ifndef THINK_C
-#include <sys/types.h>
-#endif
-#endif
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-
-#include "Gtypes.h"
-#include "libstki.h"
-#include "apat.h"
-
-#define POP             PopiOut
-#define PUSH            PushiIn
-#define TOPCURS         CursiToTop
-#define DOWNREAD        ReadiDown
-
-#define KRONECK(x, msk) ((~x & msk) ? 0 : 1)
-#define MIN(x, y)       ((x) < (y)  ? (x) : (y))
-
-/* -------------------------------------------- */
-/* Construction de la matrice S                 */
-/* -------------------------------------------- */
-
-int CreateS(Pattern *ppat, Int32 lalpha)
-{
-        Int32     i, j, indx; 
-        UInt32    pindx, amask, omask, *smat;
-
-        ppat->ok = Faux;
-
-        omask = 0x0L;
-
-        if (! (smat = NEWN(UInt32, lalpha)))
-            return 0;
-
-        for (i = 0 ; i < lalpha ; i++)
-           smat[i] = 0x0;
-
-        for (i = ppat->patlen - 1, amask = 0x1L ; i >= 0 ; i--, amask <<= 1) {
-
-            indx = ppat->patcode[i];
-
-            if (ppat->patcode[i] & OBLIBIT)
-                omask |= amask;
-
-            for (j = 0, pindx = 0x1L ; j < lalpha ; j++, pindx <<= 1)
-                if (indx & pindx)
-                    smat[j] |= amask;
-        }
-
-        ppat->smat = smat;
-
-        ppat->omask = omask;
-
-        ppat->ok = Vrai;
-
-        return 1;
-
-}
-
-/* -------------------------------------------- */
-/* Baeza-Yates/Manber algorithm                 */
-/* NoError                                      */
-/* -------------------------------------------- */
-Int32 ManberNoErr(Seq *pseq, Pattern *ppat, int patnum,int begin,int length)
-{
-        UInt32     pos;
-        UInt32    smask, r;
-        UInt8     *data;
-        StackiPtr *stkpos, *stkerr;
-        UInt32    end;
-        
-        end = begin + length;
-        end = (end <= (size_t)(pseq->seqlen+pseq->circular)) ? end:(size_t)(pseq->seqlen+pseq->circular);
-        
-
-                                        /* create local masks   */
-
-        smask = r = 0x1L << ppat->patlen;
-
-                                        /* init. scan           */
-        data   = pseq->data + begin;
-        stkpos = pseq->hitpos + patnum;
-        stkerr = pseq->hiterr + patnum;
-
-                                        /* loop on text data    */
-                                        
-        for (pos = begin ; pos < end ; pos++) {
-
-            r = (r >> 1) & ppat->smat[*data++];
-
-            if (r & 0x1L) {
-                PUSH(stkpos, pos - ppat->patlen + 1);
-                PUSH(stkerr, 0);
-            }
-
-            r |= smask;
-        }
-        
-        return (*stkpos)->top;  /* aka # of hits        */
-}
-
-/* -------------------------------------------- */
-/* Baeza-Yates/Manber algorithm                 */
-/* Substitution only                            */
-/*                                              */
-/* Note : r array is stored as :                */
-/*    0 0 r(0,j) r(0,j+1) r(1,j) r(1,j+1) ...   */
-/*                                              */
-/* -------------------------------------------- */
-Int32 ManberSub(Seq *pseq, Pattern *ppat, int patnum,int begin,int length)
-{
-        int       e, emax, found;
-        UInt32     pos;
-        UInt32    smask, cmask, sindx;
-        UInt32    *pr, r[2 * MAX_PAT_ERR + 2];
-        UInt8     *data;
-        StackiPtr *stkpos, *stkerr;
-        UInt32    end;
-        
-        end = begin + length;
-        end = (end <= (size_t)(pseq->seqlen+pseq->circular)) ? end:(size_t)(pseq->seqlen+pseq->circular);
-
-                                        /* create local masks   */
-        emax = ppat->maxerr;
-
-        r[0] = r[1] = 0x0;
-
-        cmask = smask = 0x1L << ppat->patlen;
-
-        for (e = 0, pr = r + 3 ; e <= emax ; e++, pr += 2)
-                *pr = cmask;
-
-        cmask = ~ ppat->omask;
-
-                                        /* init. scan           */
-        data   = pseq->data + begin;
-        stkpos = pseq->hitpos + patnum;
-        stkerr = pseq->hiterr + patnum;
-
-                                        /* loop on text data    */
-                                        
-        for (pos = begin ; pos < end ; pos++) {
-
-            sindx  = ppat->smat[*data++];
-
-            for (e = found = 0, pr = r ; e <= emax ; e++, pr += 2) {
-                
-                pr[2]  = pr[3] | smask;
-
-                pr[3]  =   ((pr[0] >> 1) & cmask)       /* sub   */
-                         | ((pr[2] >> 1) & sindx);      /* ident */
-
-                if (pr[3] & 0x1L) {                     /* found */
-                    if (! found)  {
-                       PUSH(stkpos, pos - ppat->patlen + 1);
-                       PUSH(stkerr, e);
-                    }
-                    found++;
-                }
-            }
-        }
-
-        return (*stkpos)->top;  /* aka # of hits        */
-}
-
-/* -------------------------------------------- */
-/* Baeza-Yates/Manber algorithm                 */
-/* Substitution + Indels                        */
-/*                                              */
-/* Note : r array is stored as :                */
-/*    0 0 r(0,j) r(0,j+1) r(1,j) r(1,j+1) ...   */
-/*                                              */
-/* Warning: may return shifted pos.             */
-/*                                              */
-/* -------------------------------------------- */
-Int32 ManberIndel(Seq *pseq, Pattern *ppat, int patnum,int begin,int length)
-{
-        int       e, emax, found;
-        UInt32     pos;
-        UInt32    smask, cmask, sindx;
-        UInt32    *pr, r[2 * MAX_PAT_ERR + 2];
-        UInt8     *data;
-        StackiPtr *stkpos, *stkerr;
-        UInt32    end;
-        
-        end = begin + length;
-        end = (end <= (size_t)(pseq->seqlen+pseq->circular)) ? end:(size_t)(pseq->seqlen+pseq->circular);
-
-                                        /* create local masks   */
-        emax = ppat->maxerr;
-
-        r[0] = r[1] = 0x0;
-
-        cmask = smask = 0x1L << ppat->patlen;
-
-        for (e = 0, pr = r + 3 ; e <= emax ; e++, pr += 2) {
-                *pr = cmask;
-                cmask = (cmask >> 1) | smask;
-        }
-
-        cmask = ~ ppat->omask;
-
-                                        /* init. scan           */
-        data   = pseq->data + begin;
-        stkpos = pseq->hitpos + patnum;
-        stkerr = pseq->hiterr + patnum;
-
-                                        /* loop on text data    */
-                                        
-        for (pos = begin ; pos < end ; pos++) {
-
-            sindx  = ppat->smat[*data++];
-
-            for (e = found = 0, pr = r ; e <= emax ; e++, pr += 2) {
-                
-                pr[2]  = pr[3] | smask;
-
-                pr[3]   =  ((     pr[0]                 /* ins   */
-                              |  (pr[0] >> 1)           /* sub   */
-                              |  (pr[1] >> 1))          /* del   */
-                            & cmask)
-                          | ((pr[2] >> 1) & sindx);     /* ident */
-
-                if (pr[3] & 0x1L) {                     /* found */
-                    if (! found) {
-                        PUSH(stkpos, pos - ppat->patlen + 1);
-                        PUSH(stkerr, e);
-                    }
-                    found++;
-                }
-
-            }
-        }
-
-        return (*stkpos)->top;  /* aka # of hits        */
-}
-
-/* -------------------------------------------- */
-/* Baeza-Yates/Manber algorithm                 */
-/* API call to previous functions               */
-/* -------------------------------------------- */
-Int32 ManberAll(Seq *pseq, Pattern *ppat, int patnum,int begin,int length)
-{
-        if (ppat->maxerr == 0)
-           return ManberNoErr(pseq, ppat, patnum, begin, length);
-        else if (ppat->hasIndel)
-           return ManberIndel(pseq, ppat, patnum, begin, length);
-        else
-           return ManberSub(pseq, ppat, patnum, begin, length);
-}
-
-
-/* -------------------------------------------- */
-/* Alignement NWS                               */
-/* pour edition des hits                        */
-/* (avec substitution obligatoire aux bords)    */
-/* -------------------------------------------- */
-
-Int32 NwsPatAlign(pseq, ppat, nerr, reslen, reserr)
-        Seq     *pseq;
-        Pattern *ppat;
-        Int32   nerr, *reslen, *reserr;
-{
-        UInt8  *sseq, *px;
-        Int32  i, j, lseq, lpat, npos, dindel, dsub,
-               *pc, *pi, *pd, *ps;
-        UInt32 amask;
-
-        static Int32 sTab[(MAX_PAT_LEN+MAX_PAT_ERR+1) * (MAX_PAT_LEN+1)];
-
-        lseq = pseq->seqlen;
-
-        pc = sTab;              /*  |----|----| --> i   */
-        pi = pc - 1;            /*  | ps | pd | |       */
-        pd = pi - lseq;         /*  |----|----| |       */
-        ps = pd - 1;            /*  | pi | pc | v j     */
-                                /*  |---------|         */
-
-        lseq = pseq->seqlen;
-        lpat = ppat->patlen;
-
-        sseq = pseq->data - 1;
-
-        amask = ONEMASK >> lpat;
-
-        for (j = 0 ; j <= lpat ; j++) {
-
-           for (i = 0 , px = sseq ; i <= lseq ; i++, px++) {
-
-              if (i && j) {
-                  dindel = MIN(*pi, *pd) + 1;
-                  dsub   = *ps + KRONECK(ppat->smat[*px], amask);
-                  *pc    = MIN(dindel, dsub);
-              }
-              else if (i)               /* j == 0       */
-                  *pc = *pi + 1;
-              else if (j)               /* i == 0       */
-                  *pc = *pd + 1;
-              else                      /* root         */
-                   *pc = 0;
-
-              pc++;
-              pi++;
-              pd++;
-              ps++;
-           }
-        
-           amask <<= 1;
-        }
-
-        pc--;
-
-        for (i = lseq, npos = 0 ; i >= 0 ; i--, pc--) {
-            if (*pc <= nerr) {
-                *reslen++ =  i;
-                *reserr++ = *pc;
-                npos++;
-             }
-        }
-
-        return npos;
-}

src/libecoPCR/Makefile

@@ -1,6 +1,5 @@
 
-SOURCES = ecoapat.c \
-         ecodna.c \
+SOURCES = ecodna.c \
          ecoError.c \
          ecoIOUtils.c \
          ecoMalloc.c \

src/libecoPCR/ecoIOUtils.c

@@ -92,7 +92,7 @@ FILE *open_ecorecorddb(const char *filename,
 {
     FILE        *f;
 	int32_t      read;
-	
+
 	f = fopen(filename,"rb");
 	
 	if (!f)

src/libecoPCR/ecoMalloc.c

@@ -2,6 +2,8 @@
 #include <stdlib.h>
 
 static int eco_log_malloc = 0;
+//static size_t eco_amount_malloc=0;
+static size_t eco_chunk_malloc=0;
 
 void    eco_trace_memory_allocation()
 {
@@ -13,53 +15,72 @@ void    eco_untrace_memory_allocation()
 	eco_log_malloc=0;
 }
 
+void ecoMallocedMemory()
+{
+	//eco_amount_malloc;
+}
 
-void   *eco_malloc(int32_t chunksize,
+void   *eco_malloc(int64_t chunksize,
                    const char *error_message,
                    const char *filename,
                    int32_t    line)
 {
 	void * chunk;
-	
+
 	chunk = calloc(1,chunksize);
-	
+
+
 	if (!chunk)
 		ecoError(ECO_MEM_ERROR,error_message,filename,line);
-		
+
+	eco_chunk_malloc++;
+
 	if (eco_log_malloc)
 		fprintf(stderr,
-			    "Memory segment located at %p of size %d is allocated (file : %s [%d])",
+			    "Memory segment located at %p of size %lld is allocated (file : %s [%d])",
 			    chunk,
 			    chunksize,
 			    filename,
 			    line);
-		
+
 	return chunk;
 }
 
 void   *eco_realloc(void *chunk,
-                    int32_t newsize,
+                    int64_t newsize,
                     const char *error_message,
                     const char *filename,
                     int32_t    line)
 {
 	void *newchunk;
-	
+
+	if (newsize == 0)
+	{
+		if (chunk)
+			free(chunk);
+		return NULL;
+	}
+
 	newchunk = realloc(chunk,newsize);
-	
+
 	if (!newchunk)
+           {
+        fprintf(stderr,"Requested memory : %lld\n",newsize);
 		ecoError(ECO_MEM_ERROR,error_message,filename,line);
+           }
+	if (!chunk)
+		eco_chunk_malloc++;
 
 	if (eco_log_malloc)
 		fprintf(stderr,
-			    "Old memory segment %p is reallocated at %p with a size of %d (file : %s [%d])",
+			    "Old memory segment %p is reallocated at %p with a size of %lld (file : %s [%d])",
 			    chunk,
 			    newchunk,
 			    newsize,
 			    filename,
 			    line);
-		
-	return newchunk;	
+
+	return newchunk;
 }
 
 void    eco_free(void *chunk,
@@ -68,7 +89,7 @@ void    eco_free(void *chunk,
                  int32_t    line)
 {
 	free(chunk);
-	
+
 	if (eco_log_malloc)
 		fprintf(stderr,
 			    "Memory segment %p is released => %s (file : %s [%d])",
@@ -76,4 +97,6 @@ void    eco_free(void *chunk,
 			    error_message,
 			    filename,
 			    line);
+
+	eco_chunk_malloc--;
 }

src/libecoPCR/ecoPCR.h

@@ -4,46 +4,45 @@
 #include <stdio.h>
 #include <inttypes.h>
 
-#ifndef H_apat
-#include "../libapat/apat.h"
-#endif
-
 /*****************************************************
- * 
+ *
  *  Data type declarations
- * 
+ *
  *****************************************************/
 
 /*
- * 
+ *
  *  Sequence types
- * 
+ *
  */
 
 typedef struct {
-	
+
 	int32_t  taxid;
 	char     AC[20];
 	int32_t  DE_length;
-	int32_t  SQ_length;
-	int32_t  CSQ_length;
-	
+	uint32_t  SQ_length;
+	uint32_t  CSQ_length; /*what is this CSQ_length ? */
+
 	char     data[1];
-	
+
 } ecoseqformat_t;
 
 typedef struct {
 	int32_t taxid;
-	int32_t SQ_length;
+	uint32_t SQ_length;
+	int32_t isexample;
 	char    *AC;
 	char    *DE;
 	char    *SQ;
-} ecoseq_t;
+	
+	int32_t	ranktaxonid;/*TR: taxon id to which the sequence belongs*/
+} ecoseq_t, *pecoseq_t;
 
 /*
- * 
+ *
  * Taxonomy taxon types
- * 
+ *
  */
 
 
@@ -53,7 +52,7 @@ typedef struct {
 	int32_t	 parent;
 	int32_t  namelength;
 	char     name[1];
-	
+
 } ecotxformat_t;
 
 typedef struct ecotxnode {
@@ -67,23 +66,23 @@ typedef struct {
 	int32_t count;
 	ecotx_t taxon[1];
 } ecotxidx_t;
- 
-	
+
+
 /*
- * 
+ *
  * Taxonomy rank types
- * 
+ *
  */
-	
+
 typedef struct {
 	int32_t count;
 	char*   label[1];
 } ecorankidx_t;
 
 /*
- * 
+ *
  * Taxonomy name types
- * 
+ *
  */
 
 typedef struct {
@@ -91,10 +90,10 @@ typedef struct {
 	int32_t  namelength;
 	int32_t  classlength;
 	int32_t  taxid;
-	char     names[1];	
+	char     names[1];
 } econameformat_t;
- 
- 
+
+
  typedef struct {
  	char 	*name;
  	char 	*classname;
@@ -102,7 +101,7 @@ typedef struct {
  	struct ecotxnode  *taxon;
 } econame_t;
 
- 
+
 typedef struct {
 	int32_t count;
 	econame_t   names[1];
@@ -115,61 +114,61 @@ typedef struct {
 	ecotxidx_t   *taxons;
 } ecotaxonomy_t;
 
- 
+
 /*****************************************************
- * 
+ *
  *  Function declarations
- * 
+ *
  *****************************************************/
 
 /*
- * 
+ *
  * Low level system functions
- * 
+ *
  */
 
 int32_t is_big_endian();
 int32_t swap_int32_t(int32_t);
 
-void   *eco_malloc(int32_t chunksize,
+void   *eco_malloc(int64_t chunksize,
                    const char *error_message,
                    const char *filename,
                    int32_t    line);
-                   
-                   
+
+
 void   *eco_realloc(void *chunk,
-                    int32_t chunksize,
+                    int64_t chunksize,
                     const char *error_message,
                     const char *filename,
                     int32_t    line);
-                    
+
 void    eco_free(void *chunk,
                  const char *error_message,
                  const char *filename,
                  int32_t    line);
-                 
+
 void    eco_trace_memory_allocation();
 void    eco_untrace_memory_allocation();
 
 #define ECOMALLOC(size,error_message) \
 	    eco_malloc((size),(error_message),__FILE__,__LINE__)
-	   
+
 #define ECOREALLOC(chunk,size,error_message) \
         eco_realloc((chunk),(size),(error_message),__FILE__,__LINE__)
-        
+
 #define ECOFREE(chunk,error_message) \
         eco_free((chunk),(error_message),__FILE__,__LINE__)
-        
-        
+
+
 
 
 /*
- * 
+ *
  * Error managment
- * 
+ *
  */
- 
-  
+
+
 void ecoError(int32_t,const char*,const char *,int);
 
 #define ECOERROR(code,message) ecoError((code),(message),__FILE__,__LINE__)
@@ -181,26 +180,26 @@ void ecoError(int32_t,const char*,const char *,int);
 
 
 /*
- * 
+ *
  * Low level Disk access functions
- * 
+ *
  */
 
 FILE *open_ecorecorddb(const char *filename,
                        int32_t    *sequencecount,
                        int32_t    abort_on_open_error);
-                       
+
 void *read_ecorecord(FILE *,int32_t *recordSize);
 
 
 
-/* 
+/*
  *   Read function in internal binary format
  */
 
 FILE             *open_ecoseqdb(const char *filename,
                                 int32_t    *sequencecount);
-                                                                
+
 ecoseq_t         *readnext_ecoseq(FILE *);
 
 ecorankidx_t     *read_rankidx(const char *filename);
@@ -211,21 +210,23 @@ econameidx_t     *read_nameidx(const char *filename,ecotaxonomy_t *taxonomy);
 
 	/**
 	 * Read taxonomy data as formated by the ecoPCRFormat.py script.
-	 * 
+	 *
 	 * This function is normaly uses internaly by the read_taxonomy
 	 * function and should not be called directly.
-	 * 
+	 *
 	 * @arg filename  path to the *.tdx file of the reformated db
-	 * 
+	 *
 	 * @return pointer to a taxonomy index structure
 	 */
- 
-ecotxidx_t       *read_taxonomyidx(const char *filename);
+
+ecotxidx_t       *read_taxonomyidx(const char *filename,const char *filename2);
 
 ecotaxonomy_t    *read_taxonomy(const char *prefix,int32_t readAlternativeName);
 
 ecotx_t *eco_findtaxonbytaxid(ecotaxonomy_t *taxonomy, int32_t taxid);
 
+ecotx_t *eco_findtaxonatrank(ecotx_t *taxon, int32_t rankidx);
+
 int eco_isundertaxon(ecotx_t *taxon, int other_taxid);
 
 ecoseq_t *ecoseq_iterator(const char *prefix);
@@ -247,11 +248,11 @@ int32_t delete_taxonomy(ecotxidx_t *index);
 
 int32_t rank_index(const char* label,ecorankidx_t* ranks);
 
-int32_t  delete_apatseq(SeqPtr pseq);
-PatternPtr buildPattern(const char *pat, int32_t error_max);
-PatternPtr complementPattern(PatternPtr pat);
-
-SeqPtr ecoseq2apatseq(ecoseq_t *in,SeqPtr out,int32_t circular);
+//int32_t  delete_apatseq(SeqPtr pseq);
+//PatternPtr buildPattern(const char *pat, int32_t error_max);
+//PatternPtr complementPattern(PatternPtr pat);
+//
+//SeqPtr ecoseq2apatseq(ecoseq_t *in,SeqPtr out,int32_t circular);
 
 char *ecoComplementPattern(char *nucAcSeq);
 char *ecoComplementSequence(char *nucAcSeq);

src/libecoPCR/ecoapat.c

@@ -19,11 +19,11 @@ void UpperSequence(char *seq)
 {
         char *cseq;
 
-        for (cseq = seq ; *cseq ; cseq++) 
+        for (cseq = seq ; *cseq ; cseq++)
             if (IS_LOWER(*cseq))
                 *cseq = TO_UPPER(*cseq);
 }
- 
+
 #undef IS_LOWER
 #undef TO_UPPER
 
@@ -54,7 +54,7 @@ void EncodeSequence(SeqPtr seq)
             data++;
             cseq++;
         }
-        
+
         for (i=0,cseq=seq->cseq;i < seq->circular; i++,cseq++,data++)
             *data = (IS_UPPER(*cseq) ? *cseq - 'A' : 0x0);
 
@@ -74,38 +74,38 @@ SeqPtr ecoseq2apatseq(ecoseq_t *in,SeqPtr out,int32_t circular)
 		{
 			out = ECOMALLOC(sizeof(Seq),
 			                "Error in Allocation of a new Seq structure");
-			          
-	        for (i  = 0 ; i < MAX_PATTERN ; i++) 
+
+	        for (i  = 0 ; i < MAX_PATTERN ; i++)
 	        {
-			   
+
 	           if (! (out->hitpos[i] = NewStacki(kMinStackiSize)))
 	             	ECOERROR(ECO_MEM_ERROR,"Error in hit stack Allocation");
-	
+
 	           if (! (out->hiterr[i] = NewStacki(kMinStackiSize)))
 	            	ECOERROR(ECO_MEM_ERROR,"Error in error stack Allocation");
 	        }
 		}
 
-		
+
 		out->name    = in->AC;
 		out->seqsiz = out->seqlen = in->SQ_length;
 		out->circular = circular;
-		
+
 		if (!out->data)
 		{
 			out->data = ECOMALLOC((out->seqlen+circular) *sizeof(UInt8),
-		    	     			  "Error in Allocation of a new Seq data member");    
+		    	     			  "Error in Allocation of a new Seq data member");
 		   	out->datsiz=  out->seqlen+circular;
 		}
 		else if ((out->seqlen +circular) >= out->datsiz)
 		{
 			out->data = ECOREALLOC(out->data,(out->seqlen+circular),
 			                      "Error during Seq data buffer realloc");
-		   	out->datsiz=  out->seqlen+circular;			                      
+		   	out->datsiz=  out->seqlen+circular;
 		}
 
 		out->cseq = in->SQ;
-		
+
 		EncodeSequence(out);
 
         return out;
@@ -117,7 +117,7 @@ int32_t delete_apatseq(SeqPtr pseq)
 
         if (pseq) {
 
-            if (pseq->data) 
+            if (pseq->data)
             	ECOFREE(pseq->data,"Freeing sequence data buffer");
 
             for (i = 0 ; i < MAX_PATTERN ; i++) {
@@ -126,74 +126,77 @@ int32_t delete_apatseq(SeqPtr pseq)
             }
 
             ECOFREE(pseq,"Freeing apat sequence structure");
-            
+
             return 0;
         }
-        
+
         return 1;
 }
 
+/*
+
 PatternPtr buildPattern(const char *pat, int32_t error_max)
 {
 	PatternPtr pattern;
 	int32_t    patlen;
-	
+
 	pattern = ECOMALLOC(sizeof(Pattern),
 						"Error in pattern allocation");
-						
+
 	pattern->ok      = Vrai;
 	pattern->hasIndel= Faux;
 	pattern->maxerr  = error_max;
 	patlen  = strlen(pat);
-	
+
 	pattern->cpat    = ECOMALLOC(sizeof(char)*patlen+1,
 	                             "Error in sequence pattern allocation");
-	                             
+
 	strncpy(pattern->cpat,pat,patlen);
 	pattern->cpat[patlen]=0;
 	UpperSequence(pattern->cpat);
-	
+
 	if (!CheckPattern(pattern))
 		ECOERROR(ECO_ASSERT_ERROR,"Error in pattern checking");
-		
+
 	if (! EncodePattern(pattern, dna))
 		ECOERROR(ECO_ASSERT_ERROR,"Error in pattern encoding");
 
    	if (! CreateS(pattern, ALPHA_LEN))
 		ECOERROR(ECO_ASSERT_ERROR,"Error in pattern compiling");
-	
+
 	return pattern;
-		
+
 }
 
 PatternPtr complementPattern(PatternPtr pat)
 {
 	PatternPtr pattern;
-	
+
 	pattern = ECOMALLOC(sizeof(Pattern),
 						"Error in pattern allocation");
-						
+
 	pattern->ok      = Vrai;
 	pattern->hasIndel= pat->hasIndel;
 	pattern->maxerr  = pat->maxerr;
 	pattern->patlen  = pat->patlen;
-	
+
 	pattern->cpat    = ECOMALLOC(sizeof(char)*(strlen(pat->cpat)+1),
 	                             "Error in sequence pattern allocation");
-	                             
+
 	strcpy(pattern->cpat,pat->cpat);
-	
+
 	ecoComplementPattern(pattern->cpat);
-	
+
 	if (!CheckPattern(pattern))
 		ECOERROR(ECO_ASSERT_ERROR,"Error in pattern checking");
-		
+
 	if (! EncodePattern(pattern, dna))
 		ECOERROR(ECO_ASSERT_ERROR,"Error in pattern encoding");
 
    	if (! CreateS(pattern, ALPHA_LEN))
 		ECOERROR(ECO_ASSERT_ERROR,"Error in pattern compiling");
-	
+
 	return pattern;
-		
+
 }
+*/

src/libecoPCR/ecofilter.c

@@ -16,4 +16,5 @@ int eco_is_taxid_included(	ecotaxonomy_t *taxonomy,
 			return 1;
 	
 	return 0;
-}
+}
+

src/libecoPCR/ecorank.c

@@ -12,14 +12,14 @@ ecorankidx_t     *read_rankidx(const char *filename)
 	int32_t      i;
 	int32_t      rs;
 	char         *buffer;
-	
+
 	f = open_ecorecorddb(filename,&count,1);
-	
+
 	index = (ecorankidx_t*) ECOMALLOC(sizeof(ecorankidx_t) + sizeof(char*) * (count-1),
 	                                  "Allocate rank index");
-	 
-	index->count=count;                                 
-	          
+
+	index->count=count;
+
 	for (i=0; i < count; i++)
 		{
 			buffer = read_ecorecord(f,&rs);
@@ -27,21 +27,18 @@ ecorankidx_t     *read_rankidx(const char *filename)
 			                                  "Allocate rank label");
 			strncpy(index->label[i],buffer,rs);
 		}
-		
+
 	return index;
 }
 
 int32_t rank_index(const char* label,ecorankidx_t* ranks)
 {
 	char **rep;
-	
+        fprintf(stderr,"Looking for rank -%s-... ",label);
 	rep = bsearch(label,ranks->label,ranks->count,sizeof(char*),compareRankLabel);
-	
+
 	if (rep)
 		return rep-ranks->label;
-	else
-		ECOERROR(ECO_NOTFOUND_ERROR,"Rank label not found");
-		
 	return -1;
 }
 

src/libecoPCR/ecoseq.c

@@ -2,6 +2,7 @@
 #include <stdlib.h>
 #include <string.h>
 #include <zlib.h>
+#include <ctype.h>
 #include <string.h>
 #include <stdio.h>
 
@@ -11,32 +12,32 @@ static FILE *open_seqfile(const char *prefix,int32_t index);
 ecoseq_t *new_ecoseq()
 {
 	void *tmp;
-	
+
 	tmp = ECOMALLOC(sizeof(ecoseq_t),"Allocate new ecoseq structure");
-	
+
 	return tmp;
 }
 
 int32_t delete_ecoseq(ecoseq_t * seq)
 {
-	
+
 	if (seq)
 	{
 		if (seq->AC)
 			ECOFREE(seq->AC,"Free sequence AC");
-			
+
 		if (seq->DE)
 			ECOFREE(seq->DE,"Free sequence DE");
-			
+
 		if (seq->SQ)
 			ECOFREE(seq->SQ,"Free sequence SQ");
-			
+
 		ECOFREE(seq,"Free sequence structure");
-		
+
 		return 0;
-			
+
 	}
-	
+
 	return 1;
 }
 
@@ -49,9 +50,9 @@ ecoseq_t *new_ecoseq_with_data( char *AC,
 	ecoseq_t *tmp;
 	int32_t lstr;
 	tmp = new_ecoseq();
-	
+
 	tmp->taxid=taxid_idx;
-	
+
 	if (AC)
 		{
 			lstr =strlen(AC);
@@ -75,6 +76,9 @@ ecoseq_t *new_ecoseq_with_data( char *AC,
 			                  "Allocate sequence data");
 			strcpy(tmp->SQ,SQ);
 		}
+
+	tmp->isexample=1;
+
 	return tmp;
 
 }
@@ -97,12 +101,14 @@ ecoseq_t *readnext_ecoseq(FILE *f)
 	int32_t  comp_status;
 	unsigned long int seqlength;
 	int32_t  rs;
-	
+	char *c;
+	int32_t i;
+
 	raw = read_ecorecord(f,&rs);
-	
+
 	if (!raw)
 		return NULL;
-		
+
 	if (is_big_endian())
 	{
 		raw->CSQ_length = swap_int32_t(raw->CSQ_length);
@@ -110,44 +116,50 @@ ecoseq_t *readnext_ecoseq(FILE *f)
 		raw->SQ_length  = swap_int32_t(raw->SQ_length);
 		raw->taxid      = swap_int32_t(raw->taxid);
 	}
-	
+
 	seq = new_ecoseq();
-	
+
 	seq->taxid = raw->taxid;
-	
+
     seq->AC    = ECOMALLOC(strlen(raw->AC) +1,
                            "Allocate Sequence Accesion number");
     strncpy(seq->AC,raw->AC,strlen(raw->AC));
 
-                           
+
     seq->DE    = ECOMALLOC(raw->DE_length+1,
                            "Allocate Sequence definition");
     strncpy(seq->DE,raw->data,raw->DE_length);
-    
+
 	seqlength = seq->SQ_length = raw->SQ_length;
-    
+
     compressed = raw->data + raw->DE_length;
 
     seq->SQ = ECOMALLOC(seqlength+1,
                         "Allocate sequence buffer");
-                        
+
+    seq->isexample=1;
+
     comp_status = uncompress((unsigned char*)seq->SQ,
                              &seqlength,
                              (unsigned char*)compressed,
                              raw->CSQ_length);
-    
+
     if (comp_status != Z_OK)
     	ECOERROR(ECO_IO_ERROR,"I cannot uncompress sequence data");
-    
+
+    for (c=seq->SQ,i=0;i<seqlength;c++,i++)
+    	*c=toupper(*c);
+
+//	fprintf(stderr,"seq name : %30s seq size : %d\n",seq->DE,seq->SQ_length);
 	return seq;
 }
 
 /**
  * Open the sequences database (.sdx file)
- * @param	prefix	name of the database (radical without extension) 
+ * @param	prefix	name of the database (radical without extension)
  * @param	index 	integer
- * 
- * @return	file object 
+ *
+ * @return	file object
  */
 FILE *open_seqfile(const char *prefix,int32_t index)
 {
@@ -161,22 +173,21 @@ FILE *open_seqfile(const char *prefix,int32_t index)
 	                           "%s_%03d.sdx",
 	                           prefix,
 	                           index);
-	                           
-			fprintf(stderr,"# Coucou %s\n",filename_buffer);
-	                           
+
+
 
 	if (filename_length >= 1024)
 		ECOERROR(ECO_ASSERT_ERROR,"file name is too long");
 
 	filename_buffer[filename_length]=0;
-		
+
 	input=open_ecorecorddb(filename_buffer,&seqcount,0);
-	
+
 	if (input)
 		fprintf(stderr,"# Reading file %s containing %d sequences...\n",
 				filename_buffer,
 				seqcount);
-				
+
 	return input;
 }
 
@@ -184,40 +195,40 @@ ecoseq_t *ecoseq_iterator(const char *prefix)
 {
 	static FILE    *current_seq_file= NULL;
 	static int32_t current_file_idx = 1;
-	static char    current_prefix[1024];
+	static char    current_prefix[1025];
 	ecoseq_t       *seq;
-	
+
 	if (prefix)
 	{
 		current_file_idx = 1;
 
 		if (current_seq_file)
 			fclose(current_seq_file);
-			
+
 		strncpy(current_prefix,prefix,1023);
 		current_prefix[1024]=0;
-		
+
 		current_seq_file = open_seqfile(current_prefix,
 		 							    current_file_idx);
-		
+
 		if (!current_seq_file)
 			return NULL;
-			
+
 	}
-	
+
 	seq = readnext_ecoseq(current_seq_file);
-	
+
 	if (!seq && feof(current_seq_file))
 	{
 		current_file_idx++;
 		fclose(current_seq_file);
 		current_seq_file = open_seqfile(current_prefix,
 		 							    current_file_idx);
-		 							    
-		
+
+
 		if (current_seq_file)
 			seq = readnext_ecoseq(current_seq_file);
 	}
-	
+
 	return seq;
-}
+}

src/libecoPCR/ecotax.c

@@ -5,28 +5,46 @@
 
 static ecotx_t *readnext_ecotaxon(FILE *f,ecotx_t *taxon);
 
- /** 
- * Open the taxonomy database 
- * @param	pointer to the database (.tdx file) 
- * @return	a ecotxidx_t structure 
+ /**
+ * Open the taxonomy database
+ * @param	pointer to the database (.tdx file)
+ * @return	a ecotxidx_t structure
  */
-ecotxidx_t     *read_taxonomyidx(const char *filename)
+ecotxidx_t     *read_taxonomyidx(const char *filename,const char *filename2)
 {
 	int32_t      count;
+	int32_t      count2;
 	FILE         *f;
+	FILE         *f2;
 	ecotxidx_t *index;
 	int32_t      i;
-	
-	f = open_ecorecorddb(filename,&count,1);
-	
-	index = (ecotxidx_t*) ECOMALLOC(sizeof(ecotxidx_t) + sizeof(ecotx_t) * (count-1),
+
+	f  = open_ecorecorddb(filename,&count,1);
+	f2 = open_ecorecorddb(filename2,&count2,0);
+
+
+
+	index = (ecotxidx_t*) ECOMALLOC(sizeof(ecotxidx_t) + sizeof(ecotx_t) * (count+count2-1),
 	                                  "Allocate taxonomy");
-	 
-	index->count=count;                                         
+
+	index->count=count+count2;
+
+	fprintf(stderr,"Reading %d taxa...\n",count);
 	for (i=0; i < count; i++){
 		readnext_ecotaxon(f,&(index->taxon[i]));
 		index->taxon[i].parent=index->taxon + (int32_t)index->taxon[i].parent;
-	}					
+	}
+
+	if (count2>0)
+		fprintf(stderr,"Reading %d local taxa...\n",count2);
+	else
+		fprintf(stderr,"No local taxon\n");
+
+	for (i=0; i < count2; i++){
+		readnext_ecotaxon(f2,&(index->taxon[count+i]));
+		index->taxon[count+i].parent=index->taxon + (int32_t)index->taxon[count+i].parent;
+	}
+
 	return index;
 }
 
@@ -34,18 +52,18 @@ ecotxidx_t     *read_taxonomyidx(const char *filename)
 int32_t delete_taxonomy(ecotxidx_t *index)
 {
 	int32_t i;
-	
+
 	if (index)
 	{
 		for (i=0; i< index->count; i++)
 			if (index->taxon[i].name)
 				ECOFREE(index->taxon[i].name,"Free scientific name");
-				
+
 		ECOFREE(index,"Free Taxonomy");
-		
+
 		return 0;
 	}
-	
+
 	return 1;
 }
 
@@ -57,32 +75,32 @@ int32_t delete_taxon(ecotx_t *taxon)
 	{
 		if (taxon->name)
 			ECOFREE(taxon->name,"Free scientific name");
-			
+
 		ECOFREE(taxon,"Free Taxon");
-		
+
 		return 0;
 	}
-		
+
 	return 1;
 }
 
 
 /**
- * Read the database for a given taxon a save the data 
+ * Read the database for a given taxon a save the data
  * into the taxon structure(if any found)
  * @param	*f	pointer to FILE type returned by fopen
  * @param	*taxon	pointer to the structure
- * 
- * @return	a ecotx_t structure if any taxon found else NULL 
+ *
+ * @return	a ecotx_t structure if any taxon found else NULL
  */
 ecotx_t *readnext_ecotaxon(FILE *f,ecotx_t *taxon)
 {
-	
+
 	ecotxformat_t *raw;
 	int32_t  rs;
-	
+
 	raw = read_ecorecord(f,&rs);
-	
+
 	if (!raw)
 		return NULL;
 
@@ -91,18 +109,18 @@ ecotx_t *readnext_ecotaxon(FILE *f,ecotx_t *taxon)
 		raw->namelength = swap_int32_t(raw->namelength);
 		raw->parent     = swap_int32_t(raw->parent);
 		raw->rank       = swap_int32_t(raw->rank);
-		raw->taxid      = swap_int32_t(raw->taxid); 
+		raw->taxid      = swap_int32_t(raw->taxid);
 	}
-	
+
 	taxon->parent = (ecotx_t*)raw->parent;
 	taxon->taxid  = raw->taxid;
 	taxon->rank   = raw->rank;
-	
+
 	taxon->name   = ECOMALLOC((raw->namelength+1) * sizeof(char),
 	                          "Allocate taxon scientific name");
-	                          
+
 	strncpy(taxon->name,raw->name,raw->namelength);
-		
+
 	return taxon;
 }
 
@@ -111,24 +129,28 @@ ecotaxonomy_t    *read_taxonomy(const char *prefix,int32_t readAlternativeName)
 {
 	ecotaxonomy_t *tax;
 	char          *filename;
+	char          *filename2;
 	int           buffsize;
-	
+
 	tax = ECOMALLOC(sizeof(ecotaxonomy_t),
 	                "Allocate taxonomy structure");
-	
+
 	buffsize = strlen(prefix)+10;
-	
+
 	filename = ECOMALLOC(buffsize,
 	                     "Allocate filename");
-	
+	filename2= ECOMALLOC(buffsize,
+	                     "Allocate filename");
+
 	snprintf(filename,buffsize,"%s.rdx",prefix);
-	
+
 	tax->ranks = read_rankidx(filename);
-	
+
 	snprintf(filename,buffsize,"%s.tdx",prefix);
-	
-	tax->taxons = read_taxonomyidx(filename);
-	
+	snprintf(filename2,buffsize,"%s.ldx",prefix);
+
+	tax->taxons = read_taxonomyidx(filename,filename2);
+
 	if (readAlternativeName)
 	{
   	   snprintf(filename,buffsize,"%s.ndx",prefix);
@@ -137,7 +159,7 @@ ecotaxonomy_t    *read_taxonomy(const char *prefix,int32_t readAlternativeName)
 	else
 	   tax->names=NULL;
 	return tax;
-	
+
 }
 
 
@@ -148,15 +170,15 @@ int32_t delete_ecotaxonomy(ecotaxonomy_t *taxonomy)
 	{
 		if (taxonomy->ranks)
 			ECOFREE(taxonomy->ranks,"Free rank index");
-			
+
 		if (taxonomy->taxons)
 			ECOFREE(taxonomy->taxons,"Free taxon index");
-			
+
 		ECOFREE(taxonomy,"Free taxonomy structure");
-		
+
 		return 0;
 	}
-	
+
 	return 1;
 }
 
@@ -165,17 +187,17 @@ ecotx_t *eco_findtaxonatrank(ecotx_t *taxon,
 {
 	ecotx_t *current_taxon;
 	ecotx_t *next_taxon;
-	
+
 	current_taxon = taxon;
 	next_taxon    = current_taxon->parent;
-	
+
 	while ((current_taxon!=next_taxon) &&  // I' am the root node
 		   (current_taxon->rank!=rankidx))
 		   {
 		   	current_taxon = next_taxon;
 		   	next_taxon    = current_taxon->parent;
 		   }
-		   
+
 	if (current_taxon->rank==rankidx)
 		return current_taxon;
 	else
@@ -185,19 +207,19 @@ ecotx_t *eco_findtaxonatrank(ecotx_t *taxon,
 /**
  * Get back information concerning a taxon from a taxonomic id
  * @param 	*taxonomy 	the taxonomy database
- * @param	taxid		the taxonomic id 
- * 
- * @result	a ecotx_t structure containing the taxonimic information  
+ * @param	taxid		the taxonomic id
+ *
+ * @result	a ecotx_t structure containing the taxonimic information
  **/
-ecotx_t *eco_findtaxonbytaxid(ecotaxonomy_t *taxonomy, 
+ecotx_t *eco_findtaxonbytaxid(ecotaxonomy_t *taxonomy,
 							  int32_t taxid)
 {
 	ecotx_t    *current_taxon;
 	int32_t     taxoncount;
 	int32_t     i;
-	
+
 	taxoncount=taxonomy->taxons->count;
-	
+
 	for (current_taxon=taxonomy->taxons->taxon,
 	     i=0;
 	     i < taxoncount;
@@ -207,30 +229,30 @@ ecotx_t *eco_findtaxonbytaxid(ecotaxonomy_t *taxonomy,
 	     	return current_taxon;
 	     }
 	 }
-	
-	return (ecotx_t*)NULL;	
+
+	return (ecotx_t*)NULL;
 }
 
 /**
  * Find out if taxon is son of other taxon (identified by its taxid)
  * @param	*taxon son 		taxon
  * @param	parent_taxid 	taxonomic id of the other taxon
- * 
+ *
  * @return 	1 is the other taxid math a parent taxid, else 0
  **/
-int eco_isundertaxon(ecotx_t *taxon, 
+int eco_isundertaxon(ecotx_t *taxon,
 						int other_taxid)
 {
 	ecotx_t *next_parent;
-	
-	next_parent = taxon->parent;	
-	
-	while ( (other_taxid != next_parent->taxid) && 
+
+	next_parent = taxon->parent;
+
+	while ( (other_taxid != next_parent->taxid) &&
 			(strcmp(next_parent->name, "root")) )
 	{
 		next_parent = next_parent->parent;
 	}
-		
+
 	if (other_taxid == next_parent->taxid)
 		return 1;
 	else
@@ -242,16 +264,16 @@ ecotx_t *eco_getspecies(ecotx_t *taxon,
 {
 	static ecotaxonomy_t *tax=NULL;
 	static int32_t		 rankindex=-1;
-	
+
 	if (taxonomy && tax!=taxonomy)
 	{
 		rankindex = rank_index("species",taxonomy->ranks);
 		tax=taxonomy;
 	}
-		
+
 	if (!tax || rankindex < 0)
 		ECOERROR(ECO_ASSERT_ERROR,"No taxonomy defined");
-		
+
 	return eco_findtaxonatrank(taxon,rankindex);
 }
 
@@ -260,16 +282,16 @@ ecotx_t *eco_getgenus(ecotx_t *taxon,
 {
 	static ecotaxonomy_t *tax=NULL;
 	static int32_t		 rankindex=-1;
-	
+
 	if (taxonomy && tax!=taxonomy)
 	{
 		rankindex = rank_index("genus",taxonomy->ranks);
 		tax=taxonomy;
 	}
-		
+
 	if (!tax || rankindex < 0)
 		ECOERROR(ECO_ASSERT_ERROR,"No taxonomy defined");
-		
+
 	return eco_findtaxonatrank(taxon,rankindex);
 }
 
@@ -279,16 +301,16 @@ ecotx_t *eco_getfamily(ecotx_t *taxon,
 {
 	static ecotaxonomy_t *tax=NULL;
 	static int32_t		 rankindex=-1;
-	
+
 	if (taxonomy && tax!=taxonomy)
 	{
 		rankindex = rank_index("family",taxonomy->ranks);
 		tax=taxonomy;
 	}
-		
+
 	if (!tax || rankindex < 0)
 		ECOERROR(ECO_ASSERT_ERROR,"No taxonomy defined");
-		
+
 	return eco_findtaxonatrank(taxon,rankindex);
 }
 
@@ -297,16 +319,16 @@ ecotx_t *eco_getkingdom(ecotx_t *taxon,
 {
 	static ecotaxonomy_t *tax=NULL;
 	static int32_t		 rankindex=-1;
-	
+
 	if (taxonomy && tax!=taxonomy)
 	{
 		rankindex = rank_index("kingdom",taxonomy->ranks);
 		tax=taxonomy;
 	}
-		
+
 	if (!tax || rankindex < 0)
 		ECOERROR(ECO_ASSERT_ERROR,"No taxonomy defined");
-		
+
 	return eco_findtaxonatrank(taxon,rankindex);
 }
 
@@ -315,15 +337,18 @@ ecotx_t *eco_getsuperkingdom(ecotx_t *taxon,
 {
 	static ecotaxonomy_t *tax=NULL;
 	static int32_t		 rankindex=-1;
-	
+
 	if (taxonomy && tax!=taxonomy)
 	{
 		rankindex = rank_index("superkingdom",taxonomy->ranks);
+                if (rankindex < 0) {
+		   rankindex = rank_index("domain",taxonomy->ranks);
+                }
 		tax=taxonomy;
 	}
-		
+
 	if (!tax || rankindex < 0)
 		ECOERROR(ECO_ASSERT_ERROR,"No taxonomy defined");
-		
+
 	return eco_findtaxonatrank(taxon,rankindex);
-}
+}

src/libecoprimer/Makefile

@@ -0,0 +1,39 @@
+
+SOURCES = goodtaxon.c \
+          readdnadb.c \
+          smothsort.c \
+          sortword.c \
+          hashsequence.c \
+          strictprimers.c \
+          aproxpattern.c \
+          merge.c \
+          queue.c \
+          libstki.c \
+          sortmatch.c \
+          pairtree.c \
+          pairs.c \
+          taxstats.c \
+          apat_search.c \
+		  filtering.c \
+		  PrimerSets.c \
+		  ahocorasick.c
+
+SRCS=$(SOURCES)
+         
+OBJECTS= $(patsubst %.c,%.o,$(SOURCES))
+
+LIBFILE= libecoprimer.a
+RANLIB= ranlib
+
+
+include ../global.mk
+
+
+all: $(LIBFILE)
+
+clean:
+	rm -rf $(OBJECTS) $(LIBFILE)
+
+$(LIBFILE): $(OBJECTS)
+	ar -cr $@ $?
+	$(RANLIB) $@

src/libecoprimer/PrimerSets.h

@@ -0,0 +1,58 @@
+#ifndef PRIMERSETS_H_
+#define PRIMERSETS_H_
+
+#include "ecoprimer.h"
+
+#define PRIMERS_IN_SET_COUNT 10
+
+typedef struct {
+	int		*set_wellIdentifiedTaxa;
+	int32_t set_pairs[PRIMERS_IN_SET_COUNT];
+	float set_specificity;
+	float set_coverage;
+	float set_lmean;
+	float set_lcov;
+	float set_score;
+	int32_t set_intaxa;
+	int32_t set_wi_cnt;
+}pairset;
+
+typedef struct{
+	ppair_t* sortedpairs;
+	int32_t sorted_count;
+	pecodnadb_t seqdb;
+	poptions_t options;
+}SetParams;
+
+typedef struct{
+	float t_spc; //specificity contribution
+	float t_cov; //coverage contribution
+	float t_lmd; //link spread difference
+	float len; //length
+	float score; //score
+}primerscore;
+
+void add_pair_in_set (pairset *pair_set, int32_t pset_idx, int32_t prb_idx, SetParams *pparams);
+void get_next_pair_options (int *pair_wi_count_sorted_ids, pairset *pair_set, SetParams *pparams);
+float get_links_distribution (int prb_idx, pairset *prob_set, SetParams *pparams);
+pairset build_primers_set_greedy_spc (SetParams *pparams);
+void get_set_mean_cov_stats (pairset *prob_set, SetParams *pparams);
+void some_other_set_possibilities (pairset *pair_set,
+		ppair_t * sortedpairs, int32_t sorted_count, pecodnadb_t seqdb, poptions_t options);
+void sets_by_SimulatedAnealing (pairset *pair_set,
+		ppair_t * sortedpairs, int32_t sorted_count, pecodnadb_t seqdb, poptions_t options);
+void sets_by_TabuSearch (pairset *pair_set,
+		ppair_t * sortedpairs, int32_t sorted_count, pecodnadb_t seqdb, poptions_t options);
+pairset * sets_by_BruteForce (ppair_t * sortedpairs, 
+		int32_t sorted_count, pecodnadb_t seqdb, poptions_t options);
+pairset * extend_set_randomly (pairset *pair_set, SetParams *params, int extend_to_cnt);
+void build_and_print_sets (ppair_t * sortedpairs, int32_t sorted_count, pecodnadb_t seqdb, poptions_t options);
+int32_t get_next_option_increasing_cov (pairset *pair_set, SetParams *pparams);
+void reset_set_props (pairset *pair_set, SetParams *pparams);
+void primers_graph_graphviz (ppair_t * sortedpairs,
+		int32_t sorted_count, poptions_t options);
+size_t primers_changeSortedArray (ppair_t ** pairs,
+		size_t sorted_count, poptions_t options);
+size_t primers_filterWithGivenLinks (ppair_t ** pairs,
+		size_t sorted_count, poptions_t options);
+#endif 

src/libecoprimer/ahocorasick.c

@@ -0,0 +1,478 @@
+/*
+ * ahocorasick.h
+ *
+ *  Created on: 26 march 2011
+ *      Author: tiayyba
+ */
+#include <inttypes.h>
+#include "hashencoder.h"
+#include "ahocorasick.h"
+
+void ahoc_graphKeywordTree (aho_state *root);
+aho_state *groot = NULL; //just for graph testing
+
+#define BASEATINDEX(w, l, i) (uint8_t)((((w)&(0x3LLU<<(((l)-(i))*2)))>>(((l)-(i))*2)) & 0x3LLU)
+
+void ahoc_addOutputElement (aho_state *node, bool_t isdirect, uint32_t idx)
+{
+	if (!node) return;
+	if (node->output.count == 0)
+		node->output.out_set = ECOMALLOC(sizeof(aho_output),
+                                      "Cannot allocate memory for aho-corasick state output element");
+	else
+		node->output.out_set = ECOREALLOC(node->output.out_set, (node->output.count+1)*sizeof(aho_output),
+                                      "Cannot allocate memory for aho-corasick state output element");
+	node->output.out_set[node->output.count].wordidx = idx;
+	node->output.out_set[node->output.count].isdirect = isdirect;
+	node->output.count++;
+}
+
+//is the passed output element in the set
+bool_t ahoc_isOutputIn (aho_state *node, aho_output ot)
+{
+	uint32_t i;
+
+	for (i=0; i<node->output.count; i++)
+		if (node->output.out_set[i].isdirect == ot.isdirect && node->output.out_set[i].wordidx == ot.wordidx) return TRUE;
+	return FALSE;
+}
+
+//take union of output of the two nodes and put in node1
+void ahoc_unionOutputElements (aho_state *node1, aho_state *node2)
+{
+	uint32_t i;
+
+	for (i=0; i<node2->output.count; i++)
+		if (ahoc_isOutputIn (node1, node2->output.out_set[i]) == FALSE)
+			ahoc_addOutputElement (node1, node2->output.out_set[i].isdirect, node2->output.out_set[i].wordidx);
+}
+
+void ahoc_addKeyword (aho_state *root, word_t w, bool_t isdirect, uint32_t idx, poptions_t options)
+{
+	uint32_t i;
+	aho_state *nextnode = root;
+	uint8_t basecode;
+	static uint32_t state_id = 0;
+	
+	//fprintf (stderr, "%s\n", ecoUnhashWord(w, options->primer_length));
+	for (i=1; i<=options->primer_length; i++)
+	{
+		basecode = BASEATINDEX (w, options->primer_length, i);
+		//fprintf (stderr, "%d", basecode);
+		if (nextnode->next[basecode] == NULL)
+		{
+			//add new state
+			nextnode->next[basecode] =  ECOMALLOC(sizeof(aho_state),
+                                      "Cannot allocate memory for aho-corasick state");
+			nextnode = nextnode->next[basecode];
+			//initialize state
+			nextnode->id = ++state_id;
+			nextnode->next[0]=nextnode->next[1]=nextnode->next[2]=nextnode->next[3]=NULL;
+			nextnode->fail = NULL;
+			nextnode->output.count = 0;
+		}
+		else
+			nextnode = nextnode->next[basecode];
+	}
+	//fprintf (stderr, "\n", basecode);
+	//new pattern addess so add node ouptup element
+	ahoc_addOutputElement (nextnode, isdirect, idx);
+}
+
+void ahoc_buildKeywordTree (aho_state *root, pwordcount_t words, poptions_t options)
+{
+	uint32_t i;
+	if (!root) return;
+
+	//init root
+	root->id = 0;
+	root->next[0]=root->next[1]=root->next[2]=root->next[3]=NULL;
+	root->fail = NULL;
+	root->output.count = 0;
+
+	//now add each word as a pattern in the keyword tree
+	for (i=0; i<words->size; i++)
+	{
+		//add direct word
+		word_t w=WORD(words->words[i]);
+		ahoc_addKeyword (root, w, TRUE, i, options);
+
+		//add reverse word
+		w=ecoComplementWord(w,options->primer_length);
+		ahoc_addKeyword (root, w, FALSE, i, options);
+	}
+
+	//loop on root if some base has no out going edge from roots
+	for (i=0; i<4; i++)
+		if (root->next[i] == NULL)
+			root->next[i] = root;
+}
+
+void ahoc_enqueue (aho_queue *ahoqueue, aho_state *node)
+{
+	queue_node *q;
+	if (node == NULL) return;
+
+	q = ECOMALLOC(sizeof(queue_node),
+            "Cannot allocate memory for aho-corasick queue node");
+	q->state_node = node;
+	q->next = NULL;
+
+	if (ahoqueue->first == NULL)
+	{
+		ahoqueue->first = q;
+		ahoqueue->last = q;
+	}
+	else
+	{
+		ahoqueue->last->next = q;
+		ahoqueue->last = q;
+	}
+}
+
+aho_state *ahoc_dequeue (aho_queue *ahoqueue)
+{
+	aho_state *node = NULL;
+	queue_node *q;
+
+	if (ahoqueue->first == NULL) return node;
+	q = ahoqueue->first;
+	ahoqueue->first = q->next;
+
+	node = q->state_node;
+	ECOFREE (q, "Cannot free memory for aho-corasick queue node");
+	return node;
+}
+
+//set fail links and output sets for the keyword tree
+void ahoc_updateForFailAndOutput (aho_state *root)
+{
+	int32_t i;
+	aho_queue Q;
+	aho_state *node_r;
+	aho_state *node_u;
+	aho_state *node_v;
+
+	//empty queue
+	Q.first = NULL;
+	Q.last = NULL;
+
+	//for us alphabet has 4 elements, A=0, C=1, G=2 and T=3
+	for (i=0; i<4; i++)
+	{
+		if (root->next[i] != root && root->next[i] != NULL)
+		{
+			root->next[i]->fail = root;
+			ahoc_enqueue (&Q, root->next[i]);
+		}
+	}
+
+	//while queue not empty
+	while (Q.first != NULL)
+	{
+		node_r = ahoc_dequeue (&Q);
+		for (i=0; i<4; i++)
+		{
+			if (node_r->next[i] != NULL)
+			{
+				node_u = node_r->next[i];
+				ahoc_enqueue (&Q, node_u);
+				node_v = node_r->fail;
+				while (node_v->next[i] == NULL) 
+					node_v = node_v->fail;
+				node_u->fail = node_v->next[i];
+				ahoc_unionOutputElements (node_u, node_u->fail);
+			}
+		}
+	}
+}
+
+void ahoc_freeKeywordTree (aho_state *node)
+{
+	int i;
+	for (i=0; i<4; i++)
+		if (node->next[i])
+			ahoc_freeKeywordTree (node->next[i]);
+	if (node->output.count > 0)
+		ECOFREE (node->output.out_set, "Free failed for node output");
+	ECOFREE (node, "Free failed for node");
+}
+
+pprimercount_t ahoc_lookforStrictPrimers (pecodnadb_t database, uint32_t seqdbsize,uint32_t exampleCount,
+		                     pwordcount_t words,poptions_t options)
+{
+	aho_state automaton_root;
+	aho_state *curr_state;
+	//uint32_t  inSequenceQuorum;
+	uint32_t  outSequenceQuorum;
+	pprimer_t  data;
+	pprimercount_t primers;
+	uint32_t   i, j, k;
+	int32_t pos;
+	uint32_t lmax;
+	char *base;
+	int8_t code;
+	uint32_t   goodPrimers=0;
+
+
+	//inSequenceQuorum = (uint32_t)floor((float)exampleCount * options->sensitivity_quorum);
+	outSequenceQuorum = (uint32_t)floor((float)(seqdbsize-exampleCount) * options->false_positive_quorum);
+
+	//fprintf(stderr,"  Primers should be at least present in %d/%d example sequences\n",inSequenceQuorum,exampleCount);
+	fprintf(stderr,"  Primers should not be present in more than %d/%d counterexample sequences\n",outSequenceQuorum,(seqdbsize-exampleCount));
+
+	data = ECOMALLOC(words->size * sizeof(primer_t),
+			         "Cannot allocate memory for fuzzy matching results");
+	for (i=0; i < words->size; i++)
+	{
+		data[i].word=WORD(words->words[i]);
+		data[i].inexample = 0;
+		data[i].outexample= 0;
+
+		data[i].directCount=ECOMALLOC(seqdbsize * sizeof(uint32_t),
+		              	"Cannot allocate memory for primer position");
+		data[i].directPos = ECOMALLOC(seqdbsize * sizeof(poslist_t),
+				"Cannot allocate memory for primer position");
+		data[i].reverseCount=ECOMALLOC(seqdbsize * sizeof(uint32_t),
+		 		"Cannot allocate memory for primer position");
+		data[i].reversePos = ECOMALLOC(seqdbsize * sizeof(poslist_t),
+				"Cannot allocate memory for primer position");
+	}
+
+	//build keywords automaton
+	ahoc_buildKeywordTree (&automaton_root, words, options);
+	//set fail links and output sets
+	ahoc_updateForFailAndOutput (&automaton_root);
+
+	//debug; print keywordtree in a gv file
+	//ahoc_graphKeywordTree (&automaton_root);
+
+	//loop on each sequence for its each base and find words
+	for (i=0; i < seqdbsize; i++)
+	{
+		if(database[i]->SQ_length <= options->primer_length) continue;
+
+		lmax = database[i]->SQ_length;
+		if (!options->circular)
+			lmax += options->primer_length-1;
+		curr_state = &automaton_root;
+
+		for (j=0,base=database[i]->SQ; j<lmax; j++,base++)
+		{
+			if (i==(uint32_t)database[i]->SQ_length) base=database[i]->SQ;
+
+			//code = encoder[(*base) - 'A'];
+			code = *base;
+			//if (iii++ < 30)
+			//	fprintf (stderr, "%d:%d,", *base, code);
+			if (code < 0 || code > 3)
+			{
+				//if error char, start from root for next character
+				//+forget any incomplete words
+				curr_state = &automaton_root;
+				continue;
+			}
+			while (curr_state->next[code] == NULL) curr_state = curr_state->fail;
+			curr_state = curr_state->next[code];
+
+			//start position of primer is options->primer_length-1 chars back
+			pos = j-options->primer_length+1;
+			if (pos < 0) pos = database[i]->SQ_length+pos;
+
+			//set output, if there is some output on this state then
+			//+all words in the output set complete here, so increment their
+			//+found properties for current sequence
+			for (k=0; k<curr_state->output.count; k++)
+			{
+				if (curr_state->output.out_set[k].isdirect)
+					data[curr_state->output.out_set[k].wordidx].directCount[i]++;
+				else
+					data[curr_state->output.out_set[k].wordidx].reverseCount[i]++;
+
+				if (options->no_multi_match)
+				{
+					if ((data[curr_state->output.out_set[k].wordidx].directCount[i] + 
+					  data[curr_state->output.out_set[k].wordidx].reverseCount[i]) > 1)
+						//since multimach not allowd, set an indication on 1st seq position that
+						//+ a multimatch was found, so that this word will be filtered out
+						//+ and because of first postion we wont have to search the whole array
+						//+ to find if it voilated nomultimatch constraint for some seq
+						data[curr_state->output.out_set[k].wordidx].directCount[0] = 2;
+					else
+					{
+						if (curr_state->output.out_set[k].isdirect)
+							//direct word found on jth position of ith sequence
+							data[curr_state->output.out_set[k].wordidx].directPos[i].value = (uint32_t)pos;
+						else
+							//reverse word found on jth position of ith sequence
+							data[curr_state->output.out_set[k].wordidx].reversePos[i].value = (uint32_t)pos;
+					}
+				}
+				else
+				{
+					//okay multi match allowed
+					if (curr_state->output.out_set[k].isdirect)
+					{
+						if (data[curr_state->output.out_set[k].wordidx].directCount[i] == 1)
+							data[curr_state->output.out_set[k].wordidx].directPos[i].value = (uint32_t)pos;
+						else
+						{
+							//need to create or extend the positions list
+							if (data[curr_state->output.out_set[k].wordidx].directCount[i] == 2)
+							{
+								//for second element, first was put in .value, so dont forget to copy that in the array too
+								data[curr_state->output.out_set[k].wordidx].directPos[i].pointer = ECOMALLOC(2 * sizeof(uint32_t),
+		                      													"Cannot allocate memory for primer position");
+								data[curr_state->output.out_set[k].wordidx].directPos[i].pointer[0] = data[curr_state->output.out_set[k].wordidx].directPos[i].value;
+								data[curr_state->output.out_set[k].wordidx].directPos[i].pointer[1] = (uint32_t)pos;
+							}
+							else
+							{
+								//for third or greater element
+								data[curr_state->output.out_set[k].wordidx].directPos[i].pointer = ECOREALLOC(data[curr_state->output.out_set[k].wordidx].directPos[i].pointer, 
+																	data[curr_state->output.out_set[k].wordidx].directCount[i] * sizeof(uint32_t),
+		                      													"Cannot allocate memory for primer position");
+								data[curr_state->output.out_set[k].wordidx].directPos[i].pointer[data[curr_state->output.out_set[k].wordidx].directCount[i]-1] = (uint32_t)pos;
+							}
+						}
+					}
+					else
+					{
+						if (data[curr_state->output.out_set[k].wordidx].reverseCount[i] == 1)
+							data[curr_state->output.out_set[k].wordidx].reversePos[i].value = (uint32_t)pos;
+						else
+						{
+							//need to create or extend the positions list
+							if (data[curr_state->output.out_set[k].wordidx].reverseCount[i] == 2)
+							{
+								//for second element, first was put in .value, so dont forget to copy that in the array too
+								data[curr_state->output.out_set[k].wordidx].reversePos[i].pointer = ECOMALLOC(2 * sizeof(uint32_t),
+		                      													"Cannot allocate memory for primer position");
+								data[curr_state->output.out_set[k].wordidx].reversePos[i].pointer[0] = data[curr_state->output.out_set[k].wordidx].reversePos[i].value;
+								data[curr_state->output.out_set[k].wordidx].reversePos[i].pointer[1] = (uint32_t)pos;
+							}
+							else
+							{
+								//for third or greater element
+								data[curr_state->output.out_set[k].wordidx].reversePos[i].pointer = ECOREALLOC(data[curr_state->output.out_set[k].wordidx].reversePos[i].pointer, 
+																	data[curr_state->output.out_set[k].wordidx].reverseCount[i] * sizeof(uint32_t),
+		                      													"Cannot allocate memory for primer position");
+								data[curr_state->output.out_set[k].wordidx].reversePos[i].pointer[data[curr_state->output.out_set[k].wordidx].reverseCount[i]-1] = (uint32_t)pos;
+							}
+						}
+					}
+				}
+				//dont forget to increment inexample or outexample count, but only once for a sequence
+				if ((data[curr_state->output.out_set[k].wordidx].directCount[i] + 
+				  data[curr_state->output.out_set[k].wordidx].reverseCount[i]) == 1)
+				{
+					if (database[i]->isexample)
+						data[curr_state->output.out_set[k].wordidx].inexample++;
+					else
+						data[curr_state->output.out_set[k].wordidx].outexample++;
+				}
+			}
+		}
+	}
+
+	//Only thing that remains is to remove the failed words
+	for (i=0,j=0; i<words->size; i++)
+	{
+		fprintf(stderr,"Primers %5d/%lld analyzed  => sequence : %s in %d example and %d counterexample sequences       \r",
+				i+1,words->size,ecoUnhashWord(data[i].word,options->primer_length),
+				data[i].inexample,data[i].outexample);
+
+		//if (data[i].inexample < inSequenceQuorum || (data[i].directCount[0] == 2 && options->no_multi_match))
+		if (data[i].directCount[0] == 2 && options->no_multi_match)
+		{
+			//bad word, delete from the array
+			for (k=0; k<seqdbsize; k++)
+			{
+				if (data[i].directCount[k] > 1)
+					ECOFREE (data[i].directPos[k].pointer, "Cannot free position pointer.");
+				if (data[i].reverseCount[k] > 1)
+					ECOFREE (data[i].reversePos[k].pointer, "Cannot free position pointer.");
+			}
+			ECOFREE (data[i].directCount, "Cannot free position pointer.");
+			ECOFREE (data[i].directPos, "Cannot free position pointer.");
+			ECOFREE (data[i].reverseCount, "Cannot free position pointer.");
+			ECOFREE (data[i].reversePos, "Cannot free position pointer.");
+		}
+		else
+		{
+			//data[i].good = data[i].inexample >= inSequenceQuorum && data[i].outexample <= outSequenceQuorum;
+			data[i].good = data[i].outexample <= outSequenceQuorum;
+			goodPrimers+=data[i].good? 1:0;
+			if (j < i)
+				data[j] = data[i];
+			j++;
+		}
+	}
+	fprintf(stderr,"\n\nOn %lld analyzed primers %d respect quorum conditions\n",words->size,goodPrimers);
+	fprintf(stderr,"Conserved primers for further analysis : %d/%lld\n",j,words->size);
+
+	primers = ECOMALLOC(sizeof(primercount_t),"Cannot allocate memory for primer table");
+	primers->primers=ECOREALLOC(data,
+			                    j * sizeof(primer_t),
+								"Cannot reallocate memory for fuzzy matching results");
+	primers->size=j;
+
+	//free memory of keyword table
+	for (i=0; i<4; i++)
+		if (automaton_root.next[i] != &automaton_root)
+			ahoc_freeKeywordTree (automaton_root.next[i]);
+
+	return primers;
+}
+
+void ahoc_graphPrintNodesInfo (aho_state *node, FILE* gfile)
+{
+	uint32_t i;
+	fprintf (gfile, "\"%d\"[\n", node->id);
+	fprintf (gfile, "label=\"%d\\n", node->id);
+	for (i=0; i<node->output.count; i++)
+		fprintf (gfile, "%d%c,", node->output.out_set[i].wordidx, node->output.out_set[i].isdirect?'d':'r');
+	fprintf (gfile, "\"\n];\n");
+
+	for (i=0; i<4; i++)
+		if (node->next[i] != NULL && node->next[i] != node)
+			ahoc_graphPrintNodesInfo (node->next[i], gfile);
+}
+
+void ahoc_graphPrintNodesLinks (aho_state *node, FILE* gfile)
+{
+	uint32_t i;
+	static int j=0;
+
+	for (i=0; i<4; i++)
+		if (node->next[i] != NULL && node->next[i] != node)
+		{
+			fprintf (gfile, "\"%d\" -> \"%d\" [\n", node->id, node->next[i]->id);
+			fprintf (gfile, "label=\"%c\"\n];\n", "ACGT"[i]);
+		}
+
+	if (j++ < 40)
+	if (node->fail != NULL && node->fail != groot)
+	{
+		fprintf (gfile, "\"%d\" -> \"%d\" [\n", node->id, node->fail->id);
+		fprintf (gfile, "color= \"red\"\n];\n");
+	}
+
+	for (i=0; i<4; i++)
+		if (node->next[i] != NULL && node->next[i] != node)
+			ahoc_graphPrintNodesLinks (node->next[i], gfile);
+}
+
+void ahoc_graphKeywordTree (aho_state *root)
+{
+	FILE *gfile;
+
+	groot=root;
+	gfile = fopen ("keywordtree.gv", "w");
+	fprintf (gfile, "digraph keywordtree {\n");
+	ahoc_graphPrintNodesInfo (root, gfile);
+	ahoc_graphPrintNodesLinks (root, gfile);
+	fprintf (gfile, "}\n");
+	fclose(gfile);
+}
+

src/libecoprimer/ahocorasick.h

@@ -0,0 +1,43 @@
+/*
+ * ahocorasick.h
+ *
+ *  Created on: 26 march 2011
+ *      Author: tiayyba
+ */
+
+#ifndef H_ahocorasick
+#define H_ahocorasick
+
+#include "ecoprimer.h"
+
+typedef struct aho_output_t{
+	uint32_t wordidx; //index of strict word (dont save the word of 64B)
+	bool_t isdirect; //we need to find both direct and reverse words so we must know which one is it
+}aho_output;
+
+typedef struct aho_output_count_t{
+	uint32_t count;
+	aho_output *out_set;
+}aho_output_count;
+
+typedef struct aho_state_t{
+	int32_t id;
+	struct aho_state_t *next[4]; //for labels A=0,C=1,G=2 and T=3
+	struct aho_state_t *fail;
+	aho_output_count output;
+}aho_state;
+
+typedef struct queue_node_t {
+	aho_state *state_node;
+	struct queue_node_t *next;
+}queue_node;
+
+typedef struct{
+	queue_node *first;
+	queue_node *last;
+}aho_queue;
+
+pprimercount_t ahoc_lookforStrictPrimers (pecodnadb_t database, uint32_t seqdbsize,uint32_t exampleCount,
+		                     pwordcount_t words,poptions_t options);
+#endif /* H_ahocorasick */
+

src/libecoprimer/amplifiatree.c

@@ -0,0 +1,131 @@
+/*
+ * amplifiatree.c
+ *
+ *  Created on: 7 mars 2009
+ *      Author: coissac
+ */
+
+#include  "ecoprimer.h"
+#include <search.h>
+
+static void cleanamplifia(pamplifia_t amplifia);
+static void deleteamplifialist(pamplifialist_t list);
+static int cmpamplifia(const void* p1,const void*p2);
+
+
+static void cleanamplifiatlist(pamplifiacount_t list)
+{
+	if (list->amplifias)
+		ECOFREE(list->amplifias,
+				"Free amplifia list");
+}
+
+static void cleanamplifia(pamplifia_t amplifia)
+{
+   cleanamplifiatlist(&(amplifia->pcr));
+}
+
+static pamplifialist_t newamplifialist(pamplifialist_t parent, size_t size)
+{
+	pamplifialist_t tmp;
+
+	tmp=ECOMALLOC(sizeof(amplifialist_t)+sizeof(amplifia_t)*(size-1),
+			      "Cannot allocate new amplifia list");
+
+	tmp->amplifiaslots=size;
+	tmp->amplifiacount=0;
+	tmp->next=NULL;
+
+	if (parent)
+		parent->next=(void*)tmp;
+
+	return tmp;
+}
+
+static void deleteamplifialist(pamplifialist_t list)
+{
+	size_t i;
+
+	if (list)
+	{
+		if (list->next)
+		{
+			deleteamplifialist(list->next);
+			list->next=NULL;
+		}
+		for (i=0; i < list->amplifiacount; i++)
+			cleanamplifia((list->amplifias)+i);
+
+		ECOFREE(list,"Delete amplifia list");
+	}
+
+}
+
+static int cmpamplifia(const void* p1,const void*p2)
+{
+	pamplifia_t pr1,pr2;
+
+	pr1=(pamplifia_t)p1;
+	pr2=(pamplifia_t)p2;
+
+	if (pr1->p1 < pr2->p1) return -1;
+	if (pr1->p1 > pr2->p1) return  1;
+
+	if (pr1->asdirect1 < pr2->asdirect1) return -1;
+	if (pr1->asdirect1 > pr2->asdirect1) return  1;
+
+	if (pr1->p2 < pr2->p2) return -1;
+	if (pr1->p2 > pr2->p2) return  1;
+
+	if (pr1->asdirect2 < pr2->asdirect2) return -1;
+	if (pr1->asdirect2 > pr2->asdirect2) return  1;
+
+	return 0;
+}
+
+pamplifia_t amplifiaintree (amplifia_t key,
+		            pamplifiatree_t amplifialist)
+{
+	if (!amplifialist->tree)
+		return NULL;
+
+	return *((pamplifia_t*)tsearch((const void *)(&key),
+								&(amplifialist->tree),
+						        cmpamplifia
+							  ));
+}
+
+pamplifia_t insertamplifia(amplifia_t key,
+		           pamplifiatree_t list)
+{
+	pamplifia_t current;
+	pamplifia_t found;
+
+	if (list->last->amplifiacount==list->last->amplifiaslots)
+	{
+		list->last->next=newamplifialist(list,100);
+		list->last=list->last->next;
+	}
+
+	current = list->last->amplifias + list->last->amplifiacount;
+	*current=key;
+
+	found = *((pamplifia_t*)tsearch((const void *)current,
+			                    &(list->tree),
+			                    cmpamplifia));
+	if (found==current)
+		list->last->amplifiacount++;
+
+	return found;
+}
+
+pamplifiatree_t initamplifiatree(pamplifiatree_t tree)
+{
+	if (!tree)
+		tree = ECOMALLOC(sizeof(amplifiatree_t),"Cannot allocate amplifia tree");
+
+	tree->first=newamplifialist(NULL,500);
+	tree->last=tree->first;
+
+	tree->tree=NULL;
+}

src/libecoprimer/apat.h

@@ -0,0 +1,120 @@
+/* ==================================================== */
+/*      Copyright (c) Atelier de BioInformatique        */
+/*      Dec. 94                                         */
+/*      File: apat.h                                    */
+/*      Purpose: pattern scan                           */
+/*      History:                                        */
+/*      28/12/94 : <Gloup> ascan first version          */
+/*      14/05/99 : <Gloup> last revision                */
+/* ==================================================== */
+
+
+#ifndef H_apat
+#define H_apat
+
+
+#include "libstki.h"
+#include "inttypes.h"
+#include "../libecoPCR/ecoPCR.h"
+
+
+/* ----------------------------------------------- */
+/* constantes                                      */
+/* ----------------------------------------------- */
+
+#ifndef BUFSIZ
+#define BUFSIZ          1024    /* io buffer size               */
+#endif
+
+#define MAX_NAME_LEN    BUFSIZ  /* max length of sequence name  */
+
+#define ALPHA_LEN        4     /* alphabet length              */
+                                /* *DO NOT* modify              */
+
+#define MAX_PATTERN       4     /* max # of patterns            */
+                                /* *DO NOT* modify              */
+
+#define MAX_PAT_LEN      32     /* max pattern length           */
+                                /* *DO NOT* modify              */
+
+#define MAX_PAT_ERR      32     /* max # of errors              */
+                                /* *DO NOT* modify              */
+
+#define PATMASK 0x3ffffff       /* mask for 26 symbols          */
+                                /* *DO NOT* modify              */
+
+#define OBLIBIT 0x4000000       /* bit 27 to 1 -> oblig. pos    */
+                                /* *DO NOT* modify              */
+
+                                /* mask for position            */
+#define ONEMASK 0x80000000      /* mask for highest position    */
+
+                                /* masks for Levenhstein edit   */
+#define OPER_IDT  0x00000000    /* identity                     */
+#define OPER_INS  0x40000000    /* insertion                    */
+#define OPER_DEL  0x80000000    /* deletion                     */
+#define OPER_SUB  0xc0000000    /* substitution                 */
+
+#define OPER_SHFT 30            /* <unused> shift               */
+
+                                /* Levenhstein Opcodes          */
+#define SOPER_IDT 0x0           /* identity                     */
+#define SOPER_INS 0x1           /* insertion                    */
+#define SOPER_DEL 0x2           /* deletion                     */
+#define SOPER_SUB 0x3           /* substitution                 */
+
+                                /* Levenhstein Opcodes masks    */
+#define OPERMASK  0xc0000000    /* mask for Opcodes             */
+#define NOPERMASK 0x3fffffff    /* negate of previous           */
+
+
+
+/* ----------------------------------------------- */
+/* data structures                                 */
+/* ----------------------------------------------- */
+
+
+typedef uint32_t pattern_t[ALPHA_LEN], *ppattern_t;
+
+                                    /* -------------------- */
+typedef struct {                    /* pattern              */
+                                    /* -------------------- */
+int      patlen;                    /* pattern length       */
+int      maxerr;                    /* max # of errors      */
+uint32_t omask;                     /* oblig. bits mask     */
+bool_t   circular;                  /* is circular sequence */
+} patternParam_t, *ppatternParam_t;
+
+
+/* ----------------------------------------------- */
+/* macros                                          */
+/* ----------------------------------------------- */
+
+#ifndef NEW
+#define NEW(typ)                (typ*)malloc(sizeof(typ))
+#define NEWN(typ, dim)          (typ*)malloc((unsigned long)(dim) * sizeof(typ))
+#define REALLOC(typ, ptr, dim)  (typ*)realloc((void *) (ptr), (unsigned long)(dim) * sizeof(typ))
+#define FREE(ptr)               free((void *) ptr)
+#endif
+
+/* ----------------------------------------------- */
+/* prototypes                                      */
+/* ----------------------------------------------- */
+
+                                /* apat_search.c        */
+
+int32_t ManberNoErr(pecoseq_t pseq,ppattern_t pat,
+                    ppatternParam_t param,
+                    StackiPtr stkpos);
+
+int32_t ManberSub(pecoseq_t pseq,ppattern_t pat,
+		          ppatternParam_t param,
+		          StackiPtr stkpos);
+
+int32_t ManberAll(pecoseq_t pseq,ppattern_t pat,
+                  ppatternParam_t param,
+                  StackiPtr stkpos);
+
+
+#endif /* H_apat */
+

src/libecoprimer/apat_parse.c

@@ -0,0 +1,65 @@
+/* ==================================================== */
+/*      Copyright (c) Atelier de BioInformatique        */
+/*      Mar. 92                                         */
+/*      File: apat_parse.c                              */
+/*      Purpose: Codage du pattern                      */
+/*      History:                                        */
+/*      00/07/94 : <Gloup> first version (stanford)     */
+/*      00/11/94 : <Gloup> revised for DNA/PROTEIN      */
+/*      30/12/94 : <Gloup> modified EncodePattern       */
+/*                         for manber search            */
+/*      14/05/99 : <Gloup> indels added                 */
+/* ==================================================== */
+
+#include <ctype.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+
+#include "apat.h"
+#include "ecoprimer.h"
+
+
+ /* IUPAC Dna      */
+static int32_t sDnaCode[]  =  {
+		/* IUPAC */
+
+		        0x00000001 /* A */, 0x0000000E /* B */, 0x00000002 /* C */,
+		        0x0000000D /* D */, 0x00000000 /* E */, 0x00000000 /* F */,
+		        0x00000004 /* G */, 0x0000000B /* H */, 0x00000000 /* I */,
+		        0x00000000 /* J */, 0x0000000C /* K */, 0x00000000 /* L */,
+		        0x00000003 /* M */, 0x0000000F /* N */, 0x00000000 /* O */,
+		        0x00000000 /* P */, 0x00000000 /* Q */, 0x00000005 /* R */,
+		        0x00000006 /* S */, 0x00000008 /* T */, 0x00000008 /* U */,
+		        0x00000007 /* V */, 0x00000009 /* W */, 0x00000000 /* X */,
+		        0x0000000A /* Y */, 0x00000000 /* Z */
+};
+
+
+/* -------------------------------------------- */
+/* internal replacement of gets                 */
+/* -------------------------------------------- */
+static char *sGets(char *buffer, int size) {
+
+        char *ebuf;
+
+        if (! fgets(buffer, size-1, stdin))
+           return NULL;
+
+        /* remove trailing line feed */
+
+        ebuf = buffer + strlen(buffer);
+
+        while (--ebuf >= buffer) {
+           if ((*ebuf == '\n') || (*ebuf == '\r'))
+                *ebuf = '\000';
+           else
+                break;
+        }
+
+        return buffer;
+}
+
+/* -------------------------------------------- */
+/* Interface                                    */
+/* -------------------------------------------- */

src/libecoprimer/apat_search.c

@@ -0,0 +1,155 @@
+/* ==================================================== */
+/*      Copyright (c) Atelier de BioInformatique        */
+/*      Dec. 94                                         */
+/*      File: apat_search.c                             */
+/*      Purpose: recherche du pattern                   */
+/*               algorithme de Baeza-Yates/Gonnet       */
+/*                             Manber (agrep)           */
+/*      History:                                        */
+/*      07/12/94 : <MFS>   first version                */
+/*      28/12/94 : <Gloup> revised version              */
+/*      14/05/99 : <Gloup> last revision                */
+/* ==================================================== */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "libstki.h"
+#include "apat.h"
+
+#define POP             PopiOut
+#define PUSH(s,v)       PushiIn(&(s),(v))
+#define TOPCURS         CursiToTop
+#define DOWNREAD        ReadiDown
+
+#define KRONECK(x, msk) ((~x & msk) ? 0 : 1)
+#define MIN(x, y)       ((x) < (y)  ? (x) : (y))
+
+
+/* -------------------------------------------- */
+/* Baeza-Yates/Manber algorithm                 */
+/* NoError                                      */
+/* -------------------------------------------- */
+int32_t ManberNoErr(pecoseq_t pseq,ppattern_t pat,
+        ppatternParam_t param,
+        StackiPtr stkpos)
+{
+        uint32_t     pos;
+        uint32_t    smask, r;
+        uint8_t     *data;
+        uint32_t    end;
+
+        end = (size_t)(pseq->SQ_length);
+
+        if (param->circular)
+        	end+=param->patlen - 1;
+
+
+                                        /* create local masks   */
+        smask = r = 0x1L << param->patlen;
+                                        /* init. scan           */
+        data   = (uint8_t*)(pseq->SQ);
+
+                                        /* loop on text data    */
+        for (pos = 0 ; pos < end ; pos++,data++) {
+        	if (pos==pseq->SQ_length)
+        		data=(uint8_t*)(pseq->SQ);
+
+        	if (*data < 4)
+        		r = (r >> 1) & pat[*data];
+        	else
+        		r=0;
+
+            if (r & 0x1L) {
+                PUSH(stkpos, pos - param->patlen + 1);
+            }
+
+            r |= smask;
+        }
+        return stkpos->top;  /* aka # of hits        */
+}
+
+/* -------------------------------------------- */
+/* Baeza-Yates/Manber algorithm                 */
+/* Substitution only                            */
+/*                                              */
+/* Note : r array is stored as :                */
+/*    0 0 r(0,j) r(0,j+1) r(1,j) r(1,j+1) ...   */
+/*                                              */
+/* -------------------------------------------- */
+int32_t ManberSub(pecoseq_t pseq,ppattern_t pat,
+        ppatternParam_t param,
+        StackiPtr stkpos)
+{
+        int       e, found;
+        uint32_t     pos;
+        uint32_t    smask, cmask, sindx;
+        uint32_t    *pr, r[2 * MAX_PAT_ERR + 2];
+        uint8_t     *data;
+        uint32_t    end;
+
+        end = (size_t)(pseq->SQ_length);
+
+        if (param->circular)
+        	end+=param->patlen - 1;
+
+                                        /* create local masks   */
+        r[0] = r[1] = 0x0;
+
+        cmask = smask = 0x1L << param->patlen;
+
+        for (e = 0, pr = r + 3 ; e <= param->maxerr ; e++, pr += 2)
+                *pr = cmask;
+
+        cmask = ~ param->omask;
+                                        /* init. scan           */
+        data   = (uint8_t*)(pseq->SQ);
+
+                                        /* loop on text data    */
+
+        for (pos = 0 ; pos < end ; pos++,data++) {
+        	if (pos==pseq->SQ_length)
+        		data=(uint8_t*)(pseq->SQ);
+
+            sindx  = (*data==4) ? 0:pat[*data];
+
+            for (e = found = 0, pr = r ; e <= param->maxerr ; e++, pr += 2) {
+
+                pr[2]  = pr[3] | smask;
+
+                pr[3]  =   ((pr[0] >> 1) & cmask)       /* sub   */
+                         | ((pr[2] >> 1) & sindx);      /* ident */
+
+                if (pr[3] & 0x1L) {                     /* found */
+                    if (! found)  {
+                       PUSH(stkpos, pos - param->patlen + 1);
+                    }
+                    found++;
+                }
+            }
+        }
+
+        return stkpos->top;  /* aka # of hits        */
+}
+
+
+/* -------------------------------------------- */
+/* Baeza-Yates/Manber algorithm                 */
+/* API call to previous functions               */
+/* -------------------------------------------- */
+int32_t ManberAll(pecoseq_t pseq,ppattern_t pat,
+        ppatternParam_t param,
+        StackiPtr stkpos)
+{
+        if (param->maxerr == 0)
+       	  return ManberNoErr(pseq,
+		                     pat, param,
+		                     stkpos);
+        else
+          return ManberSub(pseq,
+                             pat, param,
+                             stkpos);
+}
+

src/libecoprimer/aproxpattern.c

@@ -0,0 +1,237 @@
+/*
+ * aproxpattern.c
+ *
+ *  Created on: 20 nov. 2008
+ *      Author: coissac
+ */
+
+
+#include  "ecoprimer.h"
+#include  "apat.h"
+#include  <math.h>
+
+static uint8_t encoder[] = {0,                                            // A
+		                   4,                                           // b
+		                   1,                                            // C
+	       	               4,4,4,                                     // d, e, f
+		                   2,                                            // G
+		                   4,4,4,4,4,4,4,4,4,4,4,4,          // h,i,j,k,l,m,n,o,p,q,r,s
+		                   3,3,                                           // T,U
+		                   4,4,4,4,4};                              // v,w,x,y,z
+
+
+ppattern_t buildPatternFromWord(word_t word, uint32_t patlen)
+{
+	static pattern_t       pattern;
+	uint32_t i;
+
+    for (i = 0 ; i < ALPHA_LEN ; i++)
+       pattern[i] = 0x0;
+
+    for (i=0;i < patlen; i++)
+	{
+		pattern[word & 3LLU] |= 1 << i;
+		word>>=2;
+	}
+
+	return pattern;
+
+}
+
+
+#ifdef IS_UPPER
+#undef IS_UPPER
+#endif
+
+/* -------------------------------------------- */
+/* encode sequence                              */
+/* IS_UPPER is slightly faster than isupper     */
+/* -------------------------------------------- */
+
+#define IS_UPPER(c) (((c) >= 'A') && ((c) <= 'Z'))
+
+void encodeSequence(ecoseq_t *seq)
+{
+     int   i;
+     uint8_t *data;
+     char  *cseq;
+
+     data = (uint8_t*)(seq->SQ);
+     cseq = seq->SQ;
+
+     for (i=0;i<(int)(seq->SQ_length);i++,data++,cseq++)
+     {
+         *data = encoder[(IS_UPPER(*cseq) ? *cseq : 'Z') - 'A'];
+     }
+}
+
+pprimercount_t lookforAproxPrimer(pecodnadb_t database, uint32_t seqdbsize,uint32_t exampleCount,
+		                     pwordcount_t words,poptions_t options)
+{
+	pprimer_t  data;
+	pprimercount_t primers;
+	ppattern_t pattern;
+	patternParam_t params;
+	uint32_t   i;
+	uint32_t   w;
+	uint32_t   j;
+	Stacki     positions;
+	uint32_t   count=1;
+	uint32_t   goodPrimers=0;
+
+	uint32_t  inSequenceQuorum;
+	uint32_t  outSequenceQuorum;
+	bool_t    conserved = TRUE;
+
+	//poslist_t ttt;
+
+
+	inSequenceQuorum = (uint32_t)floor((float)exampleCount * options->sensitivity_quorum);
+	outSequenceQuorum = (uint32_t)floor((float)(seqdbsize-exampleCount) * options->false_positive_quorum);
+
+	fprintf(stderr,"  Primers should be at least present in %d/%d example sequences\n",inSequenceQuorum,exampleCount);
+	fprintf(stderr,"  Primers should not be present in more than %d/%d counterexample sequences\n",outSequenceQuorum,(seqdbsize-exampleCount));
+
+	data = ECOMALLOC(words->size * sizeof(primer_t),
+			         "Cannot allocate memory for fuzzy matching results");
+
+	params.circular = options->circular;
+	params.maxerr   = options->error_max;
+//	params.omask    = (1 << options->strict_three_prime) -1;
+	params.omask    = 0;
+	params.patlen   = options->primer_length;
+
+	positions.val=NULL;
+
+	for (i=0,w=0; i < words->size; i++)
+	{
+		data[w].word=WORD(words->words[i]);
+		data[w].inexample = 0;
+		data[w].outexample= 0;
+		count = 1;
+
+		if (conserved)
+		{
+			data[w].directCount=ECOMALLOC(seqdbsize * sizeof(uint32_t),
+                                      "Cannot allocate memory for primer position");
+			data[w].directPos = ECOMALLOC(seqdbsize * sizeof(poslist_t),
+				                      "Cannot allocate memory for primer position");
+			data[w].reverseCount=ECOMALLOC(seqdbsize * sizeof(uint32_t),
+	                                      "Cannot allocate memory for primer position");
+			data[w].reversePos = ECOMALLOC(seqdbsize * sizeof(poslist_t),
+					                       "Cannot allocate memory for primer position");
+		}
+
+		pattern = buildPatternFromWord(data[w].word,options->primer_length);
+		positions.val=NULL;
+
+		for (j=0; j < seqdbsize && (count < 2 || !options->no_multi_match); j++)
+		{
+			positions.cursor=0;
+			positions.top =0;
+			if (!positions.val)
+			{
+				positions.size=1;
+				positions.val = ECOMALLOC(sizeof(uint32_t),
+                                      "Cannot allocate memory for primer position");
+			}
+
+
+			count = ManberAll(database[j],pattern,&params,&positions);
+			data[w].directCount[j]=count;
+
+
+			if (count>1)
+			{
+				data[w].directPos[j].pointer = (uint32_t*)positions.val;
+				positions.val=NULL;
+			}
+			else
+			{
+				data[w].directPos[j].pointer=NULL;
+				if (count==1)
+					data[w].directPos[j].value = (uint32_t)*(positions.val);
+			}
+
+
+		}
+
+		pattern = buildPatternFromWord(ecoComplementWord(data[w].word,options->primer_length),
+				                       options->primer_length);
+
+		for (j=0; j < seqdbsize && (count < 2 || !options->no_multi_match); j++)
+		{
+			positions.cursor=0;
+			positions.top =0;
+			if (!positions.val)
+			{
+				positions.size=1;
+				positions.val = ECOMALLOC(sizeof(uint32_t),
+                                      "Cannot allocate memory for primer position");
+			}
+
+			count = ManberAll(database[j],pattern,&params,&positions);
+			data[w].reverseCount[j]=count;
+
+			if (count>1)
+			{
+				data[w].reversePos[j].pointer = (uint32_t*)positions.val;
+				positions.val=NULL;
+			}
+			else
+			{
+				data[w].reversePos[j].pointer=NULL;
+				if (count==1)
+					data[w].reversePos[j].value = (uint32_t)*(positions.val);
+			}
+
+			if (database[j]->isexample)
+			{
+				data[w].inexample+=(data[w].directCount[j] || data[w].reverseCount[j])? 1:0;
+			}
+			else
+			{
+				data[w].outexample+=(data[w].directCount[j] || data[w].reverseCount[j])? 1:0;
+
+			}
+
+            count+=data[w].directCount[j];
+		}
+
+		data[w].good = data[w].inexample >= inSequenceQuorum && data[w].outexample <= outSequenceQuorum;
+		goodPrimers+=data[w].good? 1:0;
+
+		fprintf(stderr,"Primers %5d/%lld analyzed  => sequence : %s in %d example and %d counterexample sequences       \r",
+				i+1,words->size,ecoUnhashWord(data[w].word,options->primer_length),
+				data[w].inexample,data[w].outexample);
+
+
+		conserved=data[w].inexample >= inSequenceQuorum;
+		conserved=conserved && (count < 2 || !options->no_multi_match);
+
+		if (conserved)
+			w++;
+	}
+
+	if (positions.val)
+		ECOFREE(positions.val,"Free stack position pointer");
+
+	if (!conserved)
+	{
+		ECOFREE(data[w].directCount,"Free direct count table");
+		ECOFREE(data[w].directPos,"Free direct count table");
+		ECOFREE(data[w].reverseCount,"Free direct count table");
+		ECOFREE(data[w].reversePos,"Free direct count table");
+	}
+
+	fprintf(stderr,"\n\nOn %lld analyzed primers %d respect quorum conditions\n",words->size,goodPrimers);
+	fprintf(stderr,"Conserved primers for further analysis : %d/%lld\n",w,words->size);
+
+	primers = ECOMALLOC(sizeof(primercount_t),"Cannot allocate memory for primer table");
+	primers->primers=ECOREALLOC(data,
+			                    w * sizeof(primer_t),
+								"Cannot reallocate memory for fuzzy matching results");
+	primers->size=w;
+
+	return primers;
+}

src/libecoprimer/debug.h

@@ -0,0 +1,29 @@
+/*
+ * debug.h
+ *
+ *  Created on: 12 nov. 2008
+ *      Author: coissac
+ */
+
+#ifndef DEBUG_H_
+#define DEBUG_H_
+
+#include <stdio.h>
+
+
+#ifdef DEBUG
+
+#define DEBUG_LOG(message,...) { \
+	  char *text; \
+	  (void)asprintf(&text,(message),##__VA_ARGS__); \
+      fprintf(stderr,"DEBUG %s (line %d) : %s\n",__FILE__,__LINE__,(text)); \
+      free(text); \
+	}
+
+#else
+
+#define DEBUG_LOG(message, ...)
+
+#endif
+
+#endif /* DEBUG_H_ */

src/libecoprimer/ecoprimer.h

@@ -0,0 +1,366 @@
+/*
+ * epsort.h
+ *
+ *  Created on: 6 nov. 2008
+ *      Author: coissac
+ */
+
+#ifndef EPSORT_H_
+#define EPSORT_H_
+
+#include <inttypes.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include "ecotype.h"
+#include "../libecoPCR/ecoPCR.h"
+#include "../libthermo/nnparams.h"
+#include "apat.h"
+
+#define DEBUG
+#include "debug.h"
+
+/****
+ * Word format used :
+ *
+ *     bit 63    : bad word -> this word should not be used
+ *     bit 62    : multi word -> this word is not uniq in at least one seq
+ *     bits 0-61 : hashed dna word of max size 31 pb
+ *                 code used for a : 00
+ *                 code used for c : 01
+ *                 code used for g : 10
+ *                 code used for t : 11
+ */
+
+typedef uint64_t word_t, *pword_t;
+
+#define WORD(x)           ((x) & 0x3FFFFFFFFFFFFFFFLLU)
+#define WORD(x)           ((x) & 0x3FFFFFFFFFFFFFFFLLU)
+
+#define ISBADWORD(x)      (((x) & 0x8000000000000000LLU) >> 63)
+#define SETBADWORD(x)     ((x)  | 0x8000000000000000LLU)
+#define RESETBADWORD(x)   ((x)  & 0x7FFFFFFFFFFFFFFFLLU)
+
+#define ISMULTIWORD(x)    (((x) & 0x4000000000000000LLU) >> 62)
+#define SETMULTIWORD(x)   ((x)  | 0x4000000000000000LLU)
+#define RESETMULTIWORD(x) ((x)  & 0xBFFFFFFFFFFFFFFFLLU)
+
+
+#define WORDMASK(s)       ((1LLU << ((s) * 2)) -1)
+#define LSHIFTWORD(x,s)    (((x) << 2) & WORDMASK(s))
+#define RSHIFTWORD(x,s)    (((x) & WORDMASK(s))>> 2)
+#define ERRORMASK(s)       ((int32_t)((1LLU << (s)) -1))
+
+#define RAPPENDBASE(x,s,c) (LSHIFTWORD((x),(s)) | (word_t)(c))
+#define LAPPENDBASE(x,s,c) (RSHIFTWORD((x),(s)) | ((word_t)((~(c)) & 3) << (((s)-1) *2)))
+
+
+#define ECO_ASSERT(x,message) if (!(x)) \
+	                           { \
+	                              fprintf(stderr,"Assertion Error in %s (line %d): %s\n", \
+	                                          __FILE__,\
+	                                          __LINE__,\
+	                                          message\
+	                                     ); \
+	                             exit(ECO_ASSERT_ERROR); \
+	                           }
+
+#define MINI(x,y) (((x) < (y)) ? (x):(y))
+#define MAXI(x,y) (((x) < (y)) ? (y):(x))
+
+#define FWORDSIZE (13)
+#define FWORDMASK WORDMASK(FWORDSIZE)
+#define FILTERWORD(x) ((uint32_t)((x) & FWORDMASK))
+#define CFILTERWORD(x,s) ((uint32_t)(((x) >> (((s)-FWORDSIZE)*2)) & FWORDMASK))
+
+
+
+typedef struct {
+	pword_t    words;
+	uint32_t    *strictcount;
+	uint32_t     inseqcount;
+	uint32_t     outseqcount;
+	uint64_t     size;
+} wordcount_t, *pwordcount_t;
+
+
+typedef union {
+	uint32_t *pointer;
+	uint32_t value;
+} poslist_t, *pposlist_t;
+
+
+/**
+ * primer_t structure store fuzzy match positions for a primer
+ * on all sequences
+ */
+
+typedef struct {
+    word_t     word;              //< code for the primer
+	uint32_t   *directCount;      //< Occurrence count on direct strand
+	pposlist_t directPos;        //< list of position list on direct strand
+
+	uint32_t   *reverseCount;     //< Occurrence count on reverse strand
+	pposlist_t reversePos;       //< list of position list on reverse strand
+
+	bool_t     good;              //< primer match more than quorum example and no
+	                              //  more counterexample quorum.
+
+	uint32_t   inexample;         //< count of example sequences matching primer
+	uint32_t   outexample;        //< count of counterexample sequences matching primer
+} primer_t, *pprimer_t;
+
+/**
+ * primercount_t structure store fuzzy match positions for all primers
+ * on all sequences as a list of primer_t
+ */
+typedef struct {
+	pprimer_t   primers;
+	uint32_t    size;
+} primercount_t, *pprimercount_t;
+
+typedef struct {
+	pprimer_t   primer;
+	uint32_t    position;
+	bool_t      strand;
+} primermatch_t, *pprimermatch_t;
+
+/*TR: Added*/
+typedef struct {
+	pprimermatch_t matches;
+	uint32_t	matchcount;
+} primermatchcount_t, *pprimermatchcount_t;
+
+typedef struct {
+	pecoseq_t  sequence;
+	bool_t     strand;
+	const char *amplifia;
+	int32_t	   length;
+	uint32_t   begin;
+	uint32_t   end;
+} amplifia_t, *pamplifia_t;
+
+typedef struct {
+	pamplifia_t amplifias;
+	uint32_t    ampcount;
+	uint32_t	ampslot;
+} amplifiacount_t, *pamplifiacount_t;
+
+typedef struct {
+	char *amplifia;
+	int32_t *taxonids;
+	uint32_t seqidcount;
+	uint32_t seqidindex;
+} ampseqset_t, *pampseqset_t;
+
+typedef struct {
+	int32_t taxonid;
+	char **amplifia;
+	uint32_t amplifiacount;
+	uint32_t amplifiaindex;
+} taxampset_t, *ptaxampset_t;
+
+typedef struct {
+	pprimer_t 		p1;
+	bool_t			asdirect1;
+	pprimer_t 		p2;
+	bool_t			asdirect2;
+
+	amplifiacount_t pcr;
+
+	uint32_t   		inexample;         //< example sequence count
+	uint32_t   		outexample;        //< counterexample sequence count
+	uint32_t   		intaxa;            //< example taxa count
+	uint32_t   		outtaxa;            //< counterexample taxa count
+	uint32_t		notwellidentifiedtaxa;
+
+	int				*wellIdentifiedSeqs; //< an array having elements equla to total seqs
+									 // values are either 0 or 1, if seq is well identified
+									 // its 1 else 0
+	int				*coveredSeqs; //< an array having elements equal to total seqs, 1 if seq is covered else 0
+
+					// these statistics are relative to inexample sequences
+
+	uint32_t 		mind;				//< minimum distance between primers
+	uint32_t		maxd;				//< maximum distance between primers
+	uint32_t        sumd;				//< distance sum
+	uint32_t        amplifiacount;
+	float			yule;
+	float			quorumin;
+	float           quorumout;
+	float           bs;
+	float           bc;
+	int32_t        refsequence;
+//
+//	uint32_t		taxsetcount;
+//	uint32_t		taxsetindex;
+//	ptaxampset_t	taxset;
+//
+//	uint32_t		oktaxoncount;
+	uint32_t		curseqid;
+	float			p1temp;			//strict primer1 melting temperature
+	float			p1mintemp;		//approx primer1 minimum melting temperature
+	float			p2temp;			//strict primer2 melting temperature
+	float			p2mintemp;		//approx primer2 minimum melting temperature
+} pair_t, *ppair_t;
+
+/*TR: Added*/
+
+typedef struct {
+	size_t  	paircount;
+	size_t      pairslots;
+    void*       next;
+	pair_t 	    pairs[1];
+} pairlist_t, *ppairlist_t;
+
+typedef struct {
+	ppairlist_t first;
+	ppairlist_t last;
+	void       *tree;
+	int32_t     count;
+} pairtree_t, *ppairtree_t;
+
+typedef struct {
+	pword_t     words;
+	uint32_t   *count;
+	uint32_t    push;
+	uint32_t    pop;
+	uint32_t    size;
+	bool_t      empty;
+	bool_t      full;
+} queue_t, *pqueue_t;
+
+typedef struct {
+	pword_t     words;
+	uint32_t   *count;
+	uint32_t    write;
+	uint32_t    read1;
+	uint32_t    read2;
+	uint32_t    size;
+} merge_t, *pmerge_t;
+
+typedef struct {
+	const char 		*amplifia;
+	bool_t     	strand;
+	int32_t	   	length;
+	int32_t		taxoncount;
+	void		*taxontree;
+}amptotaxon_t, *pamptotaxon_t;
+
+typedef struct {
+	int32_t 	taxid;
+	void		*amptree;
+}taxontoamp_t, *ptaxontoamp_t;
+
+typedef struct {
+	bool_t         printAC;
+	bool_t         statistics;
+	bool_t         filtering;
+	uint32_t        lmin;                   //**< Amplifia minimal length
+	uint32_t        lmax;                   //**< Amplifia maximal length
+	uint32_t        error_max;              //**< maximum error count in fuzzy search
+	uint32_t        primer_length;          //**< minimal length of the primers
+	int32_t		  *restricted_taxid;      //**< limit amplification below these taxid
+	int32_t       *ignored_taxid;         //**< no amplification below these taxid
+	int32_t		  *exception_taxid;
+	char          *prefix;
+	char          *reference;
+	pecoseq_t      refseq;
+	uint32_t       refseqid;
+	uint32_t		  circular;
+	uint32_t        doublestrand;
+	float         strict_quorum;
+	float         strict_exclude_quorum;
+	float         sensitivity_quorum;
+	float         false_positive_quorum;
+	uint32_t      strict_three_prime;
+	int32_t		  r;                      //**< count of restrited taxa (restricted_taxid array size)
+	int32_t		  g;                      //**< count of ignored taxa (ignored_taxid array size)
+	int32_t		  e;                      //**< count of ignored taxa (ignored_taxid array size)
+	bool_t        no_multi_match;
+	char		  taxonrank[20];  //TR to count ranks against a pair
+	int32_t		  taxonrankidx;   //TR to count ranks against a pair
+
+	// Some statistics useful for options filters
+
+	int32_t       dbsize;
+	int32_t		  insamples;
+	int32_t		  outsamples;
+	int32_t       intaxa;
+	int32_t       outtaxa;
+	int			  saltmethod;
+	float		  salt;
+	PNNParams	  pnparm;
+	bool_t		  print_sets_of_primers;
+	float		  specificity_threshold;
+	int			  links_cnt;
+	float		  max_links_percent;
+	bool_t		  filter_on_links;
+} options_t, *poptions_t;
+
+typedef ecoseq_t  **pecodnadb_t;
+
+void sortword(pword_t table,uint32_t N);
+
+
+pecodnadb_t readdnadb(const char *name, ecotaxonomy_t *taxonomy, uint32_t *size,poptions_t options);
+
+int isGoodTaxon(ecotaxonomy_t *taxonomy,int32_t taxon,poptions_t options);
+int isExampleTaxon(ecotaxonomy_t *taxonomy,int32_t taxon,poptions_t options);
+int isCounterExampleTaxon(ecotaxonomy_t *taxonomy,int32_t taxon,poptions_t options);
+
+uint32_t ecoWordCount(uint32_t wordsize, uint32_t circular, ecoseq_t *seq);
+pword_t ecoHashSequence(pword_t dest, uint32_t wordsize, uint32_t circular, uint32_t doublestrand, ecoseq_t *seq,uint32_t *size,int32_t *neededWords,uint32_t neededWordCount,
+	    int32_t quorum);
+uint32_t ecoCompactHashSequence(pword_t dest,uint32_t size);
+const char* ecoUnhashWord(word_t word,uint32_t size);
+word_t ecoComplementWord(word_t word,uint32_t size);
+uint32_t ecoFindWord(pwordcount_t table,word_t word);
+
+
+void ecomerge(pwordcount_t data,uint32_t s1,uint32_t s2,uint32_t remainingSeq,uint32_t seqQuorum);
+pwordcount_t initCountTable(pwordcount_t table, uint32_t wordsize, uint32_t circular, uint32_t doublestrand,uint32_t seqQuorum,ecoseq_t *seq,int32_t *neededWords,uint32_t neededWordCount);
+void addSeqToWordCountTable(pwordcount_t table, uint32_t wordsize, uint32_t circular, uint32_t doublestrand,uint32_t exampleCount,uint32_t seqQuorum,ecoseq_t *seq,int32_t *neededWords,uint32_t neededWordCount);
+
+pqueue_t newQueue(pqueue_t queue, uint32_t size);
+pqueue_t resizeQueue(pqueue_t queue, uint32_t size);
+
+void pop(pqueue_t queue);
+void push(pqueue_t queue, word_t word, uint32_t count);
+
+pqueue_t cleanQueue(pqueue_t queue);
+
+pwordcount_t lookforStrictPrimer(pecodnadb_t database, uint32_t seqdbsize,
+		                         uint32_t exampleCount, poptions_t options);
+uint32_t filterMultiStrictPrimer(pwordcount_t strictprimers);
+
+void encodeSequence(ecoseq_t *seq);
+ppattern_t buildPatternFromWord(word_t word, uint32_t patlen);
+
+pprimercount_t lookforAproxPrimer(pecodnadb_t database, uint32_t seqdbsize,uint32_t exampleCount,
+		                     pwordcount_t words,poptions_t options);
+
+void sortmatch(pprimermatch_t table,uint32_t N);
+
+ppairtree_t initpairtree(ppairtree_t tree);
+ppair_t pairintree (pair_t key,ppairtree_t pairlist);
+ppair_t insertpair(pair_t key,ppairtree_t list);
+
+
+/*TR: Added*/
+ppairtree_t buildPrimerPairs(pecodnadb_t seqdb,uint32_t seqdbsize,pprimercount_t primers,poptions_t options);
+
+int32_t counttaxon(int32_t taxid);
+int32_t getrankdbstats(pecodnadb_t seqdb,
+					   uint32_t seqdbsize,
+					   ecotaxonomy_t *taxonomy,
+					   poptions_t options);
+float taxonomycoverage(ppair_t pair, poptions_t options, pecodnadb_t seqdb,uint32_t seqdbsize);
+char ecoComplementChar(char base);
+void taxonomyspecificity (ppair_t pair, pecodnadb_t seqdb,uint32_t seqdbsize);
+
+int32_t *filteringSeq(pecodnadb_t database, uint32_t seqdbsize,
+					 uint32_t exampleCount,poptions_t options,uint32_t *size,int32_t  sequenceQuorum);
+
+void printSeqTest(pecodnadb_t seqdb,uint32_t seqdbsize);
+
+#endif /* EPSORT_H_ */

src/libecoprimer/ecotype.h

@@ -0,0 +1,14 @@
+/*
+ * ecotype.h
+ *
+ *  Created on: 24 nov. 2008
+ *      Author: coissac
+ */
+
+#ifndef ECOTYPE_H_
+#define ECOTYPE_H_
+
+typedef enum { FALSE=0,TRUE=1} bool_t, *pbool_t;
+
+
+#endif /* ECOTYPE_H_ */

src/libecoprimer/filtering.c

@@ -0,0 +1,222 @@
+/*
+ * filtering.c
+ *
+ *  Created on: 12 mai 2009
+ *      Author: coissac
+ */
+
+#include "ecoprimer.h"
+#include <string.h>
+#include <math.h>
+
+#include "hashencoder.h"
+
+static int32_t *ecoFilteringHashSequence(int32_t *dest,
+										 uint32_t circular,
+										 uint32_t doublestrand,
+										 ecoseq_t *seq,
+										 uint32_t *size);
+
+
+
+
+
+static	int32_t *ecoFilteringHashSequence(int32_t *dest,
+										  uint32_t circular,
+										  uint32_t doublestrand,
+										  ecoseq_t *seq,
+										  uint32_t *size)
+{
+
+	/*
+	 * This function aims at building a vector of count for every possible hash codes
+	 *
+	 * The function must be applied once on each sequence
+	 *
+	 * The function allocates memory on the first call for the dest table
+	 * The function also allocates memory for the static temporary table in_last_seq and
+	 * the function must be called with *dest == -1 in order to free this temporary table
+	 *
+	 */
+	static char    *in_last_seq=NULL;
+	uint32_t i=0;
+	uint32_t j;
+	char *base;
+	int8_t code;
+	int32_t error=0;
+	word_t word=0;
+	word_t antiword=0;
+	uint32_t goodword;
+	uint32_t lmax=0;
+
+	// run on the first call;
+
+	if (dest==(void*)-1)
+	{
+		if (in_last_seq) ECOFREE(in_last_seq,"Free in last seq table");
+		return NULL;
+	}
+
+
+	/* FWORDSIZE = 13 => *size =    67 108 864
+	 * FWORDSIZE = 14 => *size =   268 435 456
+	 * FWORDSIZE = 15 => *size = 1 073 741 824
+	 */
+
+	*size = pow(4,FWORDSIZE);
+
+	/*
+	 * in_last_seq is a vector of char as it is just to avoid counting twice (or more) a hash code (a DNA word)
+	 * it is set to zero (with memset) and then filled with ones for each word belonging to the sequence
+	 */
+
+	if (!in_last_seq)
+		in_last_seq = ECOMALLOC(*size*sizeof(char),
+								"Cannot allocate filtering hash table");
+