(* Mathematica Package *)

(* COPYRIGHT (C) 2010 Wynand Verwoerd 
under the terms of the GNU General Public License V3
This program is distributed in the hope that it will be useful,but WITHOUT ANY WARRANTY; 
without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
See the GNU General Public License for more details.
*)

BeginPackage["FileReadWrite`",{"NetSplitter`"}]
(* Exported symbols added here with SymbolName::usage *)  
Nonzeroes::usage="Nonzeroes[y] gives the number of non-zero entries in vector y"
Zerows::usage="Zerows[A] gives a list of zero rows in matrix A"
Nonzerows::usage="Nonzerows[A] gives a list of non-zero rows in matrix A"
GetInputS::usage="GetInputS Calls readers for input files and digests the results"
GetExternals::usage="GetExternals[Exfile] reads external compounds"
GetTargets::usage="GetTargtes[Targetfile] reads target compounds"
SubExport::usage="SubEport creates individual subnets and controls their storage in separate files"
TSVReader::usage="Reads text format stoichiometry input files" 
SBMLReader::usage="Reads stoichiometry input files in Sys Bio Markup Lang format"
TSVExport::usage="Writes text file specification for a single subnet"
SBMLExport::usage="Writes SBML file specification for a single subnet"

(* Variables only shared with the Printout package, defined here to avoid cluttering
   the list of generally shared "semi-global" variables defined in the Netsplitter` context *)
Fluxmess;Externalmess;Targetmess;
RemovedRXN;RemovedComp;DropRXN;
Reversecount;ConnectionExternalrows;Extras;

Begin["`Private`"] (* Begin Private Context *) 

Nonzeroes[y_]:=Count[y,x_/;x!= 0]

Zerows[A_]:=Flatten[Position[Map[Nonzeroes,A],x_/;x==0]]

Nonzerows[A_]:=Flatten[Position[Map[Nonzeroes,A],x_/;x!= 0]]

GetInputS[]:=Module[{Extens,duplis,NewIDs,Extrarows,ExternalIDs,ncomp,nreact,nrows,ncols,falsereverse,nabsent},Off[Remove::"remal"];
	SplitDone=False;Status="Processing the input files";StartTime=TimeUsed[];
Extens=If[ToUpperCase[Infile]=="DEMO",Infile=ToFileName[NotebookDirectory[],"Demo.tsv"];"TSV",ToUpperCase@StringDrop[Infile,{1,Last[StringPosition[Infile,"."]][[1]]}]];
Switch[Extens,"TSV",TSVreader[Infile],"SBML",SBMLreader[Infile],_,DialogInput[DialogNotebook[{TextCell["Unable to read input file - Only .TSV or .SBML files can be processed"],Button["OK",DialogReturn[1]]}]];Return[False]];

(* Compound/Reaction numbers converted to strings, embedded quotes removed *)
Complist=ToString/@Complist;Complist=Map[StringReplace[#,"\""->""]&,Complist];
StoichioExternals=ToString/@StoichioExternals;StoichioExternals=Map[StringReplace[#,"\""->""]&,StoichioExternals];
CompNames=Map[StringReplace[#,"\""->""]&,CompNames];
Reactlist=ToString/@Reactlist;Reactlist=Map[StringReplace[#,"\""->""]&,Reactlist];
Reverselist=ToString/@Reverselist;Reverselist=Map[StringReplace[#,"\""->""]&,Reverselist];

(* Do a few data integrity checks *)
ncomp=Length[Complist]; nreact=Length[Reactlist]; nrows=Dimensions[S][[1]];ncols=Dimensions[S][[2]];
duplis=Length[Complist]-Length@Union[Complist];
If[duplis>0,DialogInput[DialogNotebook[{TextCell["Invalid data - Duplicate metabolite IDs: \n"<>ToString@Take[Commonest[Complist],Min[3,duplis]]<>" \netc. were found. Check that metabolites listed as external \nin the stoichiometry file are not listed under internals as well!"],Button["OK",DialogReturn[1]]}]];Return[False]];
duplis=Length[Reactlist]-Length@Union[Reactlist];
If[duplis>0,DialogInput[DialogNotebook[{TextCell["Invalid data - Duplicate reaction IDs: \n"<>ToString@Take[Commonest[Reactlist],Min[3,duplis]]<>" \netc. were found."],Button["OK",DialogReturn[1]]}]];Return[False]];
If[ncomp!= nrows,DialogInput[DialogNotebook[{TextCell["Invalid data : "<>ToString[ncomp]<>" metabolite names listed but "<>ToString[nrows]<>" rows in S matrix "],Button["OK",DialogReturn[1]]}]];Return[False]];
If[nreact!= ncols,DialogInput[DialogNotebook[{TextCell["Invalid data : "<>ToString[nreact]<>" reaction names listed but "<>ToString[ncols]<>" columns in S matrix "],Button["OK",DialogReturn[1]]}]];Return[False]];
falsereverse=Complement[Reverselist,Reactlist];nabsent=Length[falsereverse];If[nabsent>0,DialogInput[DialogNotebook[{TextCell["Invalid data - listed reversibles "<>ToString@Take[falsereverse,Min[3,nabsent]]<>" etc. not found on the full reaction list "],Button["OK",DialogReturn[1]]}]];Return[False]];


Newlist={};Refusals={};Stepcount=0;RepeatMerge=False;OldMerge={};Exfile=If[StringLength[Exfile]== 0," ",Exfile];
Externalmess=If[ToString[FileType[Exfile]]== "File",GetExternals[Exfile];"\n"<>ToString[Length[Newlist]]<>" external metabolites read from input file: "<>Exfile,"\nNo external metabolites input file was read."];

RemovedRXN={};RemovedComp={};DropComps={};
DropRXN=Position[Reactlist,x_/;StringMatchQ[ToString[x],Map[___~~#~~___&,StringSplit[Dropstring,"|"]]],{1}];
If[DropRXN!= {},RemovedRXN=Reactlist[[Flatten[DropRXN]]];Reactlist=Delete[Reactlist,DropRXN];S=Transpose@Delete[Transpose[S],DropRXN];
	DropComps=Split@Zerows[S];S=Delete[S,DropComps];RemovedComp=Complist[[Flatten[DropComps]]];
	Complist=Delete[Complist,DropComps];CompIDs=Delete[CompIDs,DropComps];CompNames=Delete[CompNames,DropComps];
	Compartlist=Delete[Compartlist,DropComps];Chargelist=Delete[Chargelist,DropComps];
	StoichioExternals=StoichioExternals\[Intersection]Complist;Reverselist=Reverselist\[Intersection]Reactlist];
Reversecount=Length[Reverselist];

Reversibles=Flatten@Map[Position[Reactlist,#]&,Reverselist];(*Keep a copy to use when saving subnets*)
OriginalReversibles=Reversibles;OriginalReverselist=Reverselist;
Fluxfile=If[StringLength[Fluxfile]== 0," ",Fluxfile];
DoFlux=If[ToString[FileType[Fluxfile]]== "File",If[StringMatchQ[Infile,___~~"Block"~~NumberString~~"."~~__],
	ChoiceDialog[Column[{TextCell["ARE YOU SURE THE FLUX FILE SHOULD BE APPLIED?"],TextCell["\nThe input file name "<> Infile<>" corresponds to a subnet file created by Netsplitter in a previous run."],
	TextCell["\nIf so, flux information used in that run is already incorporated and repeating it will cause chaos."]}],{"Ignore flux file"->False,"Load flux file"->True},WindowFloating->True,ButtonBoxOptions->{ImageMargins->20}],True],False];
Fluxmess=If[DoFlux,Fluxes=Which[StringCount[ToUpperCase[Fluxfile],".XLS"]== 1,Import[Fluxfile][[1]],StringCount[ToUpperCase[Fluxfile],".DIF"]== 1,Import[Fluxfile],StringCount[ToUpperCase[Fluxfile],".CSV"]== 1,Import[Fluxfile],StringCount[ToUpperCase[Fluxfile],".TXT"]== 1,Import[Fluxfile,"Table"],StringCount[ToUpperCase[Fluxfile],".TSV"]== 1,Import[Fluxfile,"Table"],True,{{"Invalid ","format"}}][[All,1;;2]];flip=0;noflip=0;If[Fluxes== {{"Invalid ","format"}},"\nWARNING - unable to process flux file format !",Map[(If[MemberQ[Reactlist,#[[1]]],rno=Position[Reactlist,#[[1]]][[1,1]];Which[#[[2]]<0,S[[All,rno]]=-S[[All,rno]];Reversibles=DeleteCases[Reversibles,rno];flip++;,#[[2]]>0,Reversibles=DeleteCases[Reversibles,rno];noflip++;]])&,Fluxes];Reverselist=Reactlist[[Reversibles]];"\nFlux values were read from file "<>Fluxfile<>"\nAs a result, "<>ToString[noflip]<>" reaction directions were confirmed and "<>ToString[flip]<>" were reversed, leaving "<>ToString[Length[Reversibles]]<>" as undetermined reversibles."],"\nNo flux values input file was read."];

CNodes=Dimensions[S][[1]];RNodes=Dimensions[S][[2]];
Externals=StoichioExternals;
Externalrows=Flatten[Table[Position[Complist,Externals[[i]]],{i,1,Length[Externals]}]];StExrows=Externalrows;

If[Convert,InputExternalrows={};OutputExternalrows={};
Inpath=DirectoryName[Infile];CurDir=Directory[];SetDirectory[Inpath];
BareName=StringDrop[Infile,{1,StringLength[DirectoryName[Infile]]}];BareName=StringDrop[BareName,{Last[StringPosition[BareName,"."]][[1]],StringLength[BareName]}];
FileTitle=BareName<>" X format conversion by Netsplitter";
Switch[Extens,"TSV",SBMLExport[FileTitle,Complement[Complist,StoichioExternals]],"SBML",TSVExport[FileTitle,Complement[Complist,StoichioExternals]]];SetDirectory[CurDir];Return[]];

CR=Abs[Sign[S]]-Sign[S];
RC=Transpose[Abs[Sign[S]]+Sign[S]];
(*Append duplicate rows/cols for reversible reactions *) If[ReverseThem,
RC=Join[RC,Transpose@CR[[All,Reversibles]]];CR=Transpose@Join[Transpose[CR],RC[[Reversibles]]]];
InternalComplist=Complist;InternalCNodes=CNodes;
InputExternalrows=Flatten[Position[Total[RC],0]]; OutputExternalrows=Flatten[Position[Total[Transpose[CR]],0]];
StructuralExternalrows=Union[InputExternalrows,OutputExternalrows];Externalrows=Union[Externalrows,StructuralExternalrows];

Connex=Total[Transpose[Abs[Sign[S]]]];ConnectionExternalrows=Ordering[Connex,Count[Connex,x_/;x>ConnexMax],Greater];
If[Refusals!= {},Refusedrows=Flatten[Table[Position[Complist,Refusals[[i]]],{i,1,Length[Refusals]}]];ConnectionExternalrows=DeleteCases[ConnectionExternalrows,x_/;MemberQ[Refusedrows,x]]];BasicExternals=Complist[[Union[StructuralExternalrows,ConnectionExternalrows]]];
Externalrows=Union[Externalrows,ConnectionExternalrows];
Removedrows=Externalrows;Externals=Complist[[Externalrows]];

(* External compounds from file are matched by ID only, name or compartment discarded *)
ExternalIDs=StringSplit[Externals,Whitespace,2][[All,1]];
Newlist=Complement[Newlist,Refusals];Extras={};Extrarows={};
If[Newlist!= {},NewIDs=StringSplit[Newlist,Whitespace,2][[All,1]];
Extras=Complement[NewIDs\[Intersection]CompIDs,ExternalIDs];Extrarows=Flatten[Table[Position[CompIDs,Extras[[i]]],{i,1,Length[Extras]}]];
Externalrows=Join[Externalrows,Extrarows];
Removedrows=Externalrows;Externals=Complist[[Externalrows]]];
Extras=Complist[[Extrarows]];

(* The following code can be uncommented to force metabolites read from the externals file to
be taken as external and stop them from being reincorporated (for debugging purposes.)
This is not easily done in an existing Smatrix.TSV input file because it would require reodering the 
matrix so the newly classified externals come last. For SBML input there is no problem, simply add
"BoundaryCondition=True" to any metabolites that are forced to be external  *)

(*StoichioExternals=Extras; Extras={};*)


TargetCompounds={};Targetfile=If[StringLength[Targetfile]== 0," ",Targetfile];
Targetmess=If[ToString[FileType[Targetfile]]== "File",GetTargets[Targetfile],"\nNo target metabolites input file was read."];

Return[True]
]

GetExternals[Exfile_]:=Module[{fulltext,Nlines,ReadConfig,Configline,configno,config,ReadMerge,Exclusion,takeit,i},
fulltext=Import[Exfile,"Lines"];Nlines=Length[fulltext];
(*Read the configuration from Externals file*)
ReadConfig=False;
Do[If[StringMatchQ[fulltext[[i]],"%%%%%%%%\tConfiguration\t%%%%%%%%"],ReadConfig=True];If[ReadConfig&&StringMatchQ[fulltext[[i]],"% "~~NumberString~~__],Configline=StringReplace[fulltext[[i]],"\t%"~~__->""];Configline=StringSplit[Configline,Whitespace,3];
configno=ToExpression[Configline[[2]]];config=If[Length[Configline]== 3,Configline[[3]],""];Switch[configno,1,ReverseThem=ToExpression[config];,2,Dropstring=config;,3,DistChoice=ToExpression[config],4,LinkChoice=ToExpression[config],5,ConnexMax=ToExpression[config],6,Fluxfile=config]];If[StringMatchQ[fulltext[[i]],"%%%%%%%%\tChoices\t%%%%%%%%"],Break[]],{i,1,Nlines}];

(*Read the merge history from Externals file*)
ReadMerge=False;
Do[If[StringMatchQ[fulltext[[i]],"%%%%%%%%\tMerges\t%%%%%%%%"],ReadMerge=True];If[ReadMerge&&StringMatchQ[fulltext[[i]],"% "~~NumberString~~__],AppendTo[OldMerge,ToExpression[StringSplit[fulltext[[i]],"\t"][[2;;3]]]];If[StringMatchQ[fulltext[[i]],StartOfString~~Except["%"]~~___],Break[]]],{i,1,Nlines}];OldMerge=Map[{If[Length[#[[1]]]>0,#[[1]],{}],If[Length[#[[2]]]>0,#[[2]],{}]}&,OldMerge];
If[OldMerge!= {},ReadSteps=Length[OldMerge];
RepeatMerge=ChoiceDialog[Column[{TextCell["A merge history of "<> ToString[ReadSteps]<>" steps was read from the externals file."],TextCell["\nSteps up to the following number are reapplied to this run."],InputField[Dynamic[ReadSteps],Number,FieldSize->4],TextCell["\nChoosing None allows externals selection to be continued.\nChoosing OK but with merge steps = 0 skips selection to perform a new merging cycle."]}],{"\tOK\t"->True,"\tNone\t"->False},WindowFloating->True,ButtonBoxOptions->{ImageMargins->20}];
OldMerge=OldMerge[[1;;ReadSteps]]];

(*Read the externals and refusals from Externals fiel*)
Exclusion=False;
For[i=1,i<= Nlines,i++,takeit=StringDrop[fulltext[[i]],Flatten[StringPosition[fulltext[[i]],"%"~~___,Overlaps->False],1]];
takeit=StringReplace[takeit, (StartOfString ~~Whitespace) | (Whitespace ~~ EndOfString) :>  ""];If[ToUpperCase[takeit]== "<REFUSALS>",Exclusion=True;takeit=""];If[StringLength[takeit]>0,If[Exclusion,AppendTo[Refusals,takeit],AppendTo[Newlist,takeit]]]];

Newlist=Map[StringReplace[#, {(StartOfString ~~Whitespace) | (Whitespace ~~ EndOfString) :>  "",Whitespace->" "}]&,Newlist];Newlist=Replace[Newlist,x_/;StringLength[x]>50:>StringTake[x,49]<>" #",1];
Refusals=Map[StringReplace[#, {(StartOfString ~~Whitespace) | (Whitespace ~~ EndOfString) :>  "",Whitespace->" "}]&,Refusals];Refusals=Replace[Refusals,x_/;StringLength[x]>50:>StringTake[x,49]<>" #",1];
]

GetTargets[Targetfile_]:=Module[{fulltext,takeit,Newlist={},NewIDs,i},TargetCompounds={};fulltext=Import[Targetfile,"Lines"];
For[i=1,i<= Dimensions[fulltext][[1]],i++,takeit=StringDrop[fulltext[[i]],Flatten[StringPosition[fulltext[[i]],"%"~~___,Overlaps->False],1]];
	takeit=StringReplace[takeit, (StartOfString ~~Whitespace) | (Whitespace ~~ EndOfString) :>  ""];If[StringLength[takeit]>0,AppendTo[Newlist,takeit]]];
Targetrows={};TargetIDs={};
If[Newlist!= {},NewIDs=StringSplit[Newlist,Whitespace,2][[All,1]];
NewIDs=NewIDs\[Intersection]CompIDs;Targetrows=Flatten[Table[Position[CompIDs,NewIDs[[i]]],{i,1,Length[NewIDs]}]]];
TargetCompounds=Complist[[Targetrows]];TargetIDs=NewIDs;
Return["\n"<>ToString[Length[TargetCompounds]]<>" target metabolites read from input file: "<>Targetfile]
]


(* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> *)
SubExport[]:=(Status="Exporting subnetwork files";Off[DeleteDirectory::"nodir"];Off[DeleteDirectory::"dirne"];
Inpath=DirectoryName[Infile];CurDir=Directory[];SetDirectory[Inpath];
DeleteDirectory["Subnetworks",DeleteContents->True];
CreateDirectory["Subnetworks"];SetDirectory["Subnetworks"];
Extlisting=Complement[PreMergeExternals,BasicExternals];
Configuration={StringJoin[{"%%%%%%%%\tConfiguration\t%%%%%%%%\n","% 1 "<> ToString[ReverseThem] <>"\t%(Reversible reactions duplicated)\n","% 2 "<>Dropstring<>"\t%(Reaction omission)\n","% 3 "<>ToString[DistChoice]<>"\t%(Distance function)\n","% 4 "<>ToString[LinkChoice]<>"\t%(Intercluster distance)\n","% 5 "<>ToString[ConnexMax]<>"\t%(Maximal internal connectivity)\n","% 6 "<>Fluxfile<>"\t%(Flux values path)"}]};
DecisionHistory={StringJoin[Prepend[Riffle[Map[Riffle[#,"\t"]&,Map[ToString,DecisionTrace,{2}]],"\n% "],"%%%%%%%%\tChoices\t%%%%%%%%\n% "]]};
MergeHistory={StringJoin[Prepend[Riffle[Map[Riffle[#,"\t"]&,Map[ToString,MergeTrace,{2}]],"\n% "],"%%%%%%%%\tMerges\t%%%%%%%%\n% "]]};
Export[Outfile,Join[{"%%% List of external compounds to recreate the subnetworks found from","%%% Input file "<>Infile<>" calculated at "<>DateString[]},Configuration,DecisionHistory,MergeHistory,Extlisting,{"<Refusals>"},Refusals],{"Text","Lines"}];
If[FinalBlockCount>1,Blocktitle=Table["Block "<>ToString[i]<>" with "<>ToString[Map[Length,FinalBlocks][[i]]]<>" internal compounds",{i,1,FinalBlockCount}];
					Blocktitle[[-1]]="Orphan Block  with "<>ToString[Length[Last[FinalBlocks]]]<>" isolated compounds";
					dest=ChoiceDialog["Which output file format?",{"Text file"->"tsv","SBML"->"sbml","Cancel"->"exit"},WindowFloating->True,ButtonBoxOptions->{ImageMargins->6}];
	Which[dest== "exit",Null,dest== "tsv",MapThread[TSVExport,{Blocktitle,FinalBlocks}],dest== "sbml",MapThread[SBMLExport,{Blocktitle,FinalBlocks}]];];SetDirectory[CurDir];
)

LastWord[stringarr_]:=Map[Reverse,StringReverse[StringSplit[StringReverse[stringarr],Whitespace,2]],{Depth[stringarr]-1}]

TSVreader[filepath_]:=Module[{fulltext,Splitlist},fulltext=StringSplit[Import[filepath,"String"],"<"~~LetterCharacter..~~">"~~Whitespace];
InputName=ImportString[fulltext[[1]],"TSV"][[1,1]];
Reactlist=Flatten[ImportString[fulltext[[2]],"TSV"]];
Reverselist=Flatten[ImportString[fulltext[[3]],"TSV"]];
Complist=Flatten[ImportString[fulltext[[4]],"TSV"]];
Complist=Map[StringReplace[#, {(StartOfString ~~Whitespace) | (Whitespace ~~ EndOfString) :>  "",Whitespace->" "}]&,Complist];
StoichioExternals=Flatten[ImportString[fulltext[[5]],"TSV"]];
StoichioExternals=Map[StringReplace[#, {(StartOfString ~~Whitespace) | (Whitespace ~~ EndOfString) :>  "",Whitespace->" "}]&,StoichioExternals];
Complist=Join[Complist,StoichioExternals];

(*Splitlist=Transpose@LastWord[Complist];*)
Splitlist=Transpose@StringSplit[Complist,Whitespace,2];
{CompIDs,Compartlist}=If[Length[Splitlist]== 1,{Splitlist[[1]],ConstantArray["SYSTEM",Length[Splitlist[[1]]]]},Splitlist];
CompartIDs=Union[Compartlist];
CompartOuters=ConstantArray["Missing",Length[CompartIDs]];
If[Length[CompartIDs]<Length[Compartlist]/2,CompNames=ConstantArray["NoName",Length[CompIDs]],
	CompNames=Compartlist;Compartlist=ConstantArray["SYSTEM",Length[CompIDs]];CompartIDs={"SYSTEM"};CompartOuters={"Missing"};];
Chargelist=ConstantArray["0",Length[CompIDs]];

Complist=Replace[Complist,x_/;StringLength[x]>50:> StringTake[x,49]<>" #"(*<>ToString[n++]*),1];
StoichioExternals=Replace[StoichioExternals,x_/;StringLength[x]>50:>StringTake[x,49]<>" #"(*<>ToString[n++]*),1];
S=ImportString[fulltext[[6]],"MTX"];]

TabsepnBlock[Namelist_]:=StringJoin[Riffle[StringJoin[Riffle[#,"\t"]]&/@Partition[Namelist,5,5,1,{}],"\n"]]

Reconcile[AllItems_,WithAttributes_,Default_:None]:=Module[{flag=1,attribute},Map[(attribute=Default;Which[flag>Length[WithAttributes],Null,MatchQ[WithAttributes[[flag]],{#,_}],attribute=WithAttributes[[flag,2]];flag++;Exit;,True,Null];attribute)&,AllItems]]

TSVExport[Title_,Blocklisting_]:=Module[{BlockIntRows,BlockIntCols,BlockRevCols,BlockRows,BlockExtRows,BlockStoi,Inflows,Outflows,Crossflows,Reactlisting,Reverselisting,Intlisting,Extlisting,Smatlisting,Outfile,i},
If[Length[Blocklisting]>0,BlockIntRows=Flatten[Table[Position[Complist,Blocklisting[[i]]],{i,1,Length[Blocklisting]}]];
	BlockIntCols=Nonzerows[Normal[Transpose[S[[BlockIntRows]]]]];
	BlockRevCols=BlockIntCols\[Intersection](*Original*)Reversibles;
	BlockRows=Nonzerows[Normal[S[[All,BlockIntCols]]]];
	BlockExtRows=Complement[BlockRows,BlockIntRows];
	Inflows=BlockExtRows\[Intersection]InputExternalrows;
	Outflows=BlockExtRows\[Intersection]OutputExternalrows;
	Crossflows=Complement[BlockExtRows,Inflows\[Union]Outflows];
	BlockExtRows=Join[Inflows,Crossflows,Outflows];
	BlockRows=Join[BlockIntRows,BlockExtRows];
	BlockStoi=S[[BlockRows,BlockIntCols]];
	Reactlisting=TabsepnBlock[Reactlist[[BlockIntCols]]];
	Reverselisting=TabsepnBlock[Reactlist[[BlockRevCols]]];
	Intlisting=TabsepnBlock[Complist[[BlockIntRows]]];
	Extlisting=TabsepnBlock[Complist[[BlockExtRows]]];
	Smatlisting=ExportString[SparseArray[BlockStoi],"MTX"]; 
	Outfile=StringJoin[StringSplit[Title][[1;;2]],".tsv"];
	Export[Outfile,{"<Title>",Title,"<Reactions>",Reactlisting,"<ReversibleReactions>",Reverselisting,"<InternalCompounds>",Intlisting,"<ExternalCompounds>",Extlisting,"<Stoichiometry>",Smatlisting},"List"];
]]

SBMLreader[filepath_]:=Module[{InputXML,ExIDs,ExNames,rxnsbml,Srules},
InputXML=Import[filepath,"XML"];
InputName=Cases[InputXML,XMLElement["model",{"id"->modid_,___},_]:>modid,Infinity][[1]];
CompartIDs=Cases[InputXML,XMLElement["compartment",{___,"id"->compartid_,___},{}]:> compartid,Infinity];
CompartOuters=Cases[InputXML,XMLElement["compartment",{___,"id"->compartid_,___,"outside"->outid_,___},{}]:> {compartid,outid},Infinity];
CompartOuters=Reconcile[CompartIDs,CompartOuters,"Missing"];
Reactlist=Cases[InputXML,XMLElement["reaction",{___,"id"->reactid_,___},{___}]:>reactid,Infinity];
Reverselist=Complement[Reactlist,Cases[InputXML,XMLElement["reaction",{___,"id"->reactid_,___,"reversible"->"false",___},{___}]:>reactid,Infinity]];
CompIDs=Cases[InputXML,XMLElement["species",{___,"id"->compid_,___},{}]:> compid,Infinity];
CompNames=Cases[InputXML,XMLElement["species",{___,"id"->compid_,___,"name"->compname_,___},{}]:> {compid,compname},Infinity];
CompNames=Reconcile[CompIDs,CompNames,"NoName"];
Compartlist=Cases[InputXML,XMLElement["species",{___,"id"->compid_,___,"compartment"->compname_,___},{}]:> {compid,compname},Infinity];
Compartlist=Reconcile[CompIDs,Compartlist,"SYSTEM"];
Chargelist=Cases[InputXML,XMLElement["species",{___,"id"->compid_,___,"charge"->compname_,___},{}]:> {compid,compname},Infinity];
Chargelist=Reconcile[CompIDs,Chargelist,"0"];

If[CompIDs[[1]]== "DUMMY",CompIDs=Rest[CompIDs];CompNames=Rest[CompNames];Compartlist=Rest[Compartlist];Chargelist=Rest[Chargelist];];

rxnsbml=Cases[InputXML,XMLElement["reaction",_,_],Infinity];
Srules=Flatten[Map[ReactionParse,rxnsbml],1];S=SparseArray[Srules];
ExIDs=Cases[InputXML,XMLElement["species",{"id"->compid_,___,"boundaryCondition"->"true",___},{}]:> compid,Infinity];
ExNames=Cases[InputXML,XMLElement["species",{"id"->compid_,___,"name"->compname_,___,"boundaryCondition"->"true",___},{}]|XMLElement["species",{"id"->compid_,___,"boundaryCondition"->"true",___,"name"->compname_,___},{}]:> {compid,compname},Infinity];
ExNames=Reconcile[ExIDs,ExNames,"NoName"];
If[CompIDs!= CompNames,StoichioExternals=MapThread[StringJoin,{ExIDs,ConstantArray[" ",Length[ExIDs]],ExNames}];
			Complist=MapThread[StringJoin,{CompIDs,ConstantArray[" ",Length[CompIDs]],CompNames}],
			StoichioExternals=ExIDs];
Complist=Map[StringReplace[#, {(StartOfString ~~Whitespace) | (Whitespace ~~ EndOfString) :>  "",Whitespace->" "}]&,Complist];
Complist=Replace[Complist,x_/;StringLength[x]>50:>StringTake[x,49]<>" #"(*<>ToString[n++]*),1];
StoichioExternals=Map[StringReplace[#, {(StartOfString ~~Whitespace) | (Whitespace ~~ EndOfString) :>  "",Whitespace->" "}]&,StoichioExternals];
StoichioExternals=Replace[StoichioExternals,x_/;StringLength[x]>50:>StringTake[x,49]<>" #"(*<>ToString[n++]*),1];
]

ReactionParse[RxnSBML_]:=Module[{Scol,SrowIn,SrowOut},
	Scol=RxnSBML/.XMLElement["reaction",{"id"->rname_,___},_]:>Position[Reactlist,rname,1][[1,1]];
	SrowIn=Cases[Cases[RxnSBML,XMLElement["listOfReactants",_,_],{2}],XMLElement["speciesReference",{"species"->spec_,("stoichiometry"->stoich_)...},{}]:> {Position[CompIDs,spec,1][[1,1]],N@If[stoich== "",1,ToExpression[stoich]]},{3}];
	SrowOut=Cases[Cases[RxnSBML,XMLElement["listOfProducts",_,_],{2}],XMLElement["speciesReference",{"species"->spec_,("stoichiometry"->stoich_)...},{}]:> {Position[CompIDs,spec,1][[1,1]],N@If[stoich== "",1,ToExpression[stoich]]},{3}];
	Join[Map[{#[[1]],Scol}->-#[[2]]&,SrowIn],Map[{#[[1]],Scol}->#[[2]]&,SrowOut]]]

ReactantSBML[i_]:=Table["<speciesReference species=\"" ~~Reactions[[i,2,j,1]]~~ "\""~~ " stoichiometry=\""~~ ToString[Reactions[[i,2,j,2]]]~~ "\" />",{j,1,Length[ Reactions[[i,2]]]}];

ProductSBML[i_]:=Table["<speciesReference species=\"" ~~Reactions[[i,3,j,1]]~~ "\""~~ " stoichiometry=\""~~ ToString[Reactions[[i,3,j,2]]]~~ "\" />",{j,1,Length[ Reactions[[i,3]]]}];

SBMLExport[Title_,Blocklisting_]:=Module[{BlockIntRows,BlockIntCols,BlockRows,BlockExtRows,Inflows,Outflows,Crossflows,BlockStoi,Header,BlockCompartIDs,BlockCompartOuters,CompartmentSBL,Speclist,SpeciesSBL,BlockReactlist,Stoicol,Reacts,Prods,ReactionSBML,Outfile,i},
Outfile=StringJoin[StringSplit[Title][[1;;2]],".sbml"];If[Length[Blocklisting]>0,BlockIntRows=Flatten[Table[Position[Complist,Blocklisting[[i]]],{i,1,Length[Blocklisting]}]];BlockIntCols=Nonzerows[Normal[Transpose[S[[BlockIntRows]]]]];BlockRows=Nonzerows[Normal[S[[All,BlockIntCols]]]];BlockExtRows=Complement[BlockRows,BlockIntRows];
Inflows=BlockExtRows\[Intersection]InputExternalrows;Outflows=BlockExtRows\[Intersection]OutputExternalrows;Crossflows=Complement[BlockExtRows,Inflows\[Union]Outflows];
BlockExtRows=Join[Inflows,Crossflows,Outflows];BlockRows=Join[BlockIntRows,BlockExtRows];
BlockStoi=S[[BlockRows,BlockIntCols]];
Header={"<?xml version=\"1.0\" encoding=\"UTF-8\"?> <sbml ","xmlns=\"http://www.sbml.org/sbml/level2\"","xmlns:sbml=\"http://www.sbml.org/sbml/level2\" version=\"1\" level=\"2\"","xmlns:math=\"http://www.w3.org/1998/Math/MathML\"","xmlns:html=\"http://www.w3.org/1999/xhtml\"   >","<model id=\"" ~~Title~~ "\" name=\"" ~~InputName ~~"\" >","<notes>Metabolic subnetwork extracted by SubnetSplitter on " ~~DateString[] ~~"</notes>","<notes>Stoichiometry file processed is: " ~~Infile ~~ "</notes>"};

BlockCompartIDs=Union[Compartlist[[BlockRows]]];BlockCompartOuters=CompartOuters[[Flatten[Table[Position[CompartIDs,BlockCompartIDs[[i]]],{i,1,Length[BlockCompartIDs]}]]]];
CompartmentSBL=Table["<compartment id=\""~~ BlockCompartIDs[[i]] ~~"\" spatialDimensions=\"3\""~~If[BlockCompartOuters[[i]]== "Missing",""," outside=\""~~ BlockCompartOuters[[i]]~~ "\""]~~" />",{i,1,Length[BlockCompartIDs]}];

Speclist=Join[CompIDs[[BlockIntRows]],CompIDs[[BlockExtRows]]];
(*Dummy is a workaround for CellnetAnalyzer bug, name field in first record not read and upsets internal/external recognition *)
SpeciesSBL=Join[{"<species id=\"DUMMY\" compartment=\"SYSTEM\" initialAmount=\"10.0\" name=\"Dummy\" />"},Table["<species id=\"" ~~CompIDs[[BlockIntRows[[i]]]]~~ "\" compartment=\"" ~~Compartlist[[BlockIntRows[[i]]]]~~"\" charge=\"" ~~Chargelist[[BlockIntRows[[i]]]]~~"\" initialAmount=\"10.0\""~~If[CompNames[[BlockIntRows[[i]]]]== "NoName",""," name=\""~~ CompNames[[BlockIntRows[[i]]]]~~ "\""]~~" />",{i,1,Length[BlockIntRows]}],Table["<species id=\"" ~~CompIDs[[BlockExtRows[[i]]]]~~ "\" compartment=\"" ~~Compartlist[[BlockExtRows[[i]]]]~~"\" charge=\"" ~~Chargelist[[BlockExtRows[[i]]]]~~"\" initialAmount=\"10.0\""~~ " boundaryCondition=\"true\""~~If[CompNames[[BlockExtRows[[i]]]]== "NoName",""," name=\""~~ CompNames[[BlockExtRows[[i]]]]~~ "\""]~~" />",{i,1,Length[BlockExtRows]}]];
BlockReactlist=Reactlist[[BlockIntCols]];BlockReactlist=Table[{BlockReactlist[[i]],If[Count[(*Original*)Reverselist,BlockReactlist[[i]]]>0,True,False]},{i,1,Length[BlockReactlist]}];
Stoicol=(Transpose@Normal@BlockStoi);
Reacts=Map[Flatten[Position[#,x_/;x<0]]&,Stoicol];Prods=Map[Flatten[Position[#,x_/;x>0]]&,Stoicol];
Reactions=Table[{BlockReactlist[[i]],Transpose[{Speclist[[Reacts[[i]]]],-Stoicol[[i,Reacts[[i]]]]}],Transpose[{Speclist[[Prods[[i]]]],Stoicol[[i,Prods[[i]]]]}]},{i,1,Length[BlockReactlist]}];
ReactionSBML=Flatten@Table[Join[{"<reaction id=\""~~ Reactions[[i,1,1]]~~"\"" ~~If[ Reactions[[i,1,2]],""," reversible=\"false\" "]~~ " name=\"" ~~ Reactions[[i,1,1]]~~ "\" >"},{"<listOfReactants>"},ReactantSBML[i],{"</listOfReactants>","<listOfProducts>"},ProductSBML[i],{"</listOfProducts>","</reaction>"}],{i,1,Length[Reactions]}];
Export[Outfile,Join[Header,{"<listOfCompartments>"},CompartmentSBL,{"</listOfCompartments>","<listOfSpecies>"},SpeciesSBL,{"</listOfSpecies>","<listOfReactions>"},ReactionSBML,{"</listOfReactions>","</model>","</sbml>"}],"List"];
]]

End[] (* End Private Context *)

EndPackage[]