(* ::Package:: *)

(* TransitionMatrix Package *)
(* Version 0.3 *)

(* Last edited on Apr 6th, 2009 *)

BeginPackage["TransitionMatrix`"];


(* Messages *)

TransitionMatrix::usage = "TransitionMatrix[matrix, options] creates a TransitionMatrix Object.";
TransitionMatrixNormalized::usage = "";
RateMatrix::usage = "RateMatrix[matrix, options] creates a RateMatrix Object.";
RateMatrixNormalized::usage = "";
CountMatrix::usage = "CountMatrix[matrix, options] creates a CountMatrix Object.";
ClusterMatrix::usage = "ClusterMatrix[matrix, options] creates a ClusterMatrix Object.";
TimeSeries::usage = "TimeSeries[matrix, options] creates a TimeSeries Object.";
Trajectory::usage = "";

ShowSize::usage = "";

RowSumOneQ::usage = "";
RowSumZeroQ::usage = "";

ReplaceRules::usage = "";
UnionRules::usage = "";
RemoveRules::usage = "";

NormalForm::usage = "";
Normal::usage = "";
Normalize::usage = "";
NormalizedQ::usage = "";

MatrixPlot::usage = "";
Length::usage = "";

VarianceMatrix::usage = "";
MeanMatrix::usage = "";

MaximumLikelihoodMatrix::usage = "";
NewStateProbabilityVector::usage = "";
DrawTransitionMatrix::usage = "";

CountMatrixQ::usage = "";
TransitionMatrixQ::usage = "";
RateMatrixQ::usage = "";
ClusterMatrixQ::usage = "";
TimeSeriesQ::usage = "";
TrajectoryQ::usage = "";

ToRateMatrix::usage = "";
ToTransitionMatrix::usage = "";
ToCountMatrix::usage = "";
ToTimeSeries::usage = "";

Combine::usage = "";
Spectrum::usage = "";
ImpliedTimeScales::usage = "";
Equilibrium::usage = "";
EquilibriumSensitivityTensor::usage = "";
MeanFirstPassageTime::usage = "";
MeanFirstPassageTimeSensitivityTensor::usage = "";
Committor::usage = "";
CommittorSensitivityMatrix::usage = "";
PotentialOfMeanForce::usage = "";
PotentialOfMeanForceBins::usage = "";
PotentialOfMeanForceSensitivityMatrix::usage = "";
PotentialOfMeanForceError::usage = "";
CommittorReactionCoordinate = "";
VisitPredictionMatrix::usage = "";

MaximumLikelihoodMatrix::NYI = "Not Yet Implemented";
MeanMatrix::NYI = "Not Yet Implemented";
VarianceMatrix::NYI = "Not Yet Implemented";
DrawTransitionMatrix::NYI = "Not Yet Implemented";
ToRateMatrix::MethodNotImplemented = "The Method has not yet been implemented. Please try Method -> MatrixLog.";
Likelihood::usage = "";

DirichletDistribution::usage = "";
DiracDistribution::usage = "";
DirichletWithDBDistribution::usage = "";

Symmetrize::usage = "";
ChooseTransition::usage = "";

KMeansClustering::usage = "";
CrispClusterMatrix::usage = "";
RemoveClusterMatrix::usage = "";
RemoveEmptyStatesClusterMatrix::usage = "";
Cluster::usage = "";
HasEmptyStatesQ::usage = "";
ToCrispRules::usage = "";
PCCASoftClusterMatrix::usage = "";
PCCACrispClusterMatrix::usage = "";
Crispify::usage = "";
InvolvedStates::usage = "";
ClosedEnd::usage = "";
RemoveVoid::usage = "";

TransitionMatrixFromMatrix::usage = "";
RateMatrixFromMatrix::usage = "";

PositionStreamPlot::usage = "";
PositionStreamDensityPlot::usage = "";
PositionDensityPlot::usage = "";
PositionListPlot::usage = "";
PositionStatePlot::usage = "";

ReactiveFluxMatrix::usage = "";
NetFluxMatrix::usage = "";
ReactivityVector::usage = "";
StateExitFlowRate::usage = "";
StateExitFlowRateSensitivityMatrix::usage = "";

FieldStyle::usage = "";
Rotation::usage = "";
Gradient::usage = "";

ReadSparseMatrix::usage = "";
WriteSparseMatrix::usage = "";

AssumeSparse = "";
Hummer::usage = "";
MatrixLog::usage = "";
VoidStates::usage = "";

Vertices::usage = "";
VertexStyle::usage = "";
VertexLabel::usage = "";

ToHummerRateMatrix6::usage = "";
ToHummerRateMatrix6Geometric::usage = "";


(* All options *)

Options[TransitionMatrix]            = {Normalized->True,ShowSize->7,Lagtime->1};
Options[RateMatrix]                  = {Normalized->True,ShowSize->7};
Options[CountMatrix]                 = {Lagtime->1,Distribution->DirichletDistribution,Prior->1,TransitionMatrix->None,ShowSize->7};
Options[TimeSeries]                  = {States->1,Lagtime->1,ShowSize->11,VoidStates->{}};
Options[ClusterMatrix]               = {};
Options[Trajectory]                  = {};

Options[ToTimeSeries]                = {AssumeSparse->False};
Options[CrispClusterMatrix]          = {Complete->False,States->0};
Options[KMeansClustering]            = {Iterations->20};
Options[ToRateMatrix]                = {Method->MatrixLog,Lagtime->Inherit};
Options[DirichletWithDBDistribution] = {Iterations->1000,MaximalResidue->0.001};
Options[MeanFirstPassageTime]        = {Lagtime->1};
Options[PositionDensityPlot]         = Join[Options[ListDensityPlot],{Vertices->False,VertexLabel->False,VertexStyle->{}}];
Options[PositionStatePlot]           = {Text->False, TextStyle->{RGBColor[1,1,1]}, Grid->True, GridStyle->{RGBColor[0,0,0],Thickness[0.012`]}, Outline->False, OutlineStyle->{}, Fill->True, FillStyle->{RGBColor[0.88`,0.47`,0]}, Coloring->False};


Begin["`Private`"];


(* Functions to handle options *)

ReplaceRules[rule_,f___] := 
		Apply[Sequence,Join[{rule},Select[List[f],#[[1]]=!=rule[[1]]&]]]
UnionRules[rules1___] := 
		Apply[Sequence,Union[List[rules1][[All,1]]]/.Map[#[[1]]->#&,List[rules1]]]
RemoveRules[rule_,f___] := 
		Apply[Sequence,Select[List[f],#[[1]]=!=rule&]]


(* Matrix functions *)

RowSumOneQ[m_?MatrixQ] := 
		Apply[And,Map[#==1&,Map[Apply[Plus,#]&,m]]]
RowSumZeroQ[m_?MatrixQ] := 
		Apply[And,Map[#==0&,Map[Apply[Plus,#]&,m]]]


(* CountMatrix functions *)

CountMatrix/:MakeBoxes[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]],StandardForm] := 
		If[OptionValue[CountMatrix,ShowSize]>Length[m],
			InterpretationBox[#,CountMatrix[m,opts]]&[RowBox[{SubscriptBox["\[DoubleStruckCapitalC]",Length[m]],ToBoxes[MatrixForm[m]]}]],
			InterpretationBox[#,CountMatrix[m,opts]]&[SubscriptBox["\[DoubleStruckCapitalC]",Length[m]]]
		]
CountMatrix/:CountMatrix[CountMatrix[m_?MatrixQ,opts1:OptionsPattern[CountMatrix]],opts2:OptionsPattern[CountMatrix]] := 
		CountMatrix[m,UnionRules[opts1,opts2]]
CountMatrix/:NormalForm[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]] := 
		{m,opts}
CountMatrix/:Normal[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]] := 
		m
CountMatrix/:Normalize[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]] := 
		TransitionMatrix[Map[#/Apply[Plus,#]&,m]]
CountMatrix/:MatrixForm[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]] := 
		MatrixForm[m]
CountMatrix/:Length[CountMatrix[m_?MatrixQ],opts:OptionsPattern[CountMatrix]] := 
		Length[m]
CountMatrix/:MatrixPlot[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]] := 
		MatrixPlot[m,Frame->True,BaseStyle->{FontFamily->"GillSans",11},FrameTicks->{Automatic,Automatic,None,None},ImageSize->320,ColorFunction->"TemperatureMap"]
CountMatrixQ[x_] := 
		If[Head[x]===CountMatrix&&MatrixQ[x[[1]]],True,False];
CountMatrix/:VarianceMatrix[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]] := 
		VarianceMatrix[OptionValue[CountMatrix,Distribution][CountMatrix[m,opts]]]
CountMatrix/:MeanMatrix[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]] := 
		MeanMatrix[OptionValue[CountMatrix,Distribution][CountMatrix[m,opts]]]
CountMatrix/:MaximumLikelihoodMatrix[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]],optsA___] := 
		MaximumLikelihoodMatrix[OptionValue[CountMatrix,Distribution][CountMatrix[m,opts]],optsA]
CountMatrix/:DrawTransitionMatrix[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]] := 
		DrawTransitionMatrix[OptionValue[CountMatrix,Distribution][CountMatrix[m,opts]]]
CountMatrix/:NewStateProbabilityVector[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]] := 
		Normal[MeanMatrix[OptionValue[CountMatrix,Distribution][CountMatrix[ReplacePart[ArrayPad[m,{{0,1},{0,1}}],{Length[m]+1,Length[m]+1}->1],opts]]]][[Range[1,Length[m]],Length[m]+1]]
Combine[v_?VectorQ] :=
		Module[{cMatrixQ, matchSizeQ},	
			cMatrixQ=Apply[And,Select[v,CountMatrixQ]];
			If[cMatrixQ,
				matchSizeQ=(Min[Map[Length,v]]==Max[Map[Length,v]]);
				If[matchSizeQ,
						Print["okay!"];
(*					CountMatrix[Apply[Plus,v[[All,1]]]]*)
				]
			]
		]
TimeSeries/:ToCountMatrix[TimeSeries[v_?VectorQ,opts:OptionsPattern[TimeSeries]],optsCM:OptionsPattern[Join[Options[CountMatrix],{Lagtime->1, ClosedEnd->True, SeriesLength->All}]]] := 
		Module[{cM,options,nStates,nLagtime,cV},			
			options=UnionRules[opts,States->Max[v]];
			nStates=Max[OptionValue[TimeSeries,States],Max[v]];
			nLagtime=Lagtime/.List[optsCM]/.{Lagtime->1};
			cM=SparseArray[Map[{Mod[#[[1]],(nStates+1)],Floor[#[[1]]/(nStates+1)]}->#[[2]]/nLagtime&,Tally[Drop[Drop[Join[(nStates+1)*v,Table[0,{nLagtime}]]+Join[Table[0,{nLagtime}],v],-nLagtime],nLagtime]]],{nStates,nStates}];
(*
			If[SeriesLength/.optsCM/.{SeriesLength->All}===All
				,
				cvPart=Partition[v,SeriesLength/.optsCM/.{SeriesLength->Length[v]},SeriesLength/.optsCM/.{SeriesLength->Length[v]},1,{}];
				cM=SparseArray[Table[Map[{Mod[#[[1]],(nStates+1)],Floor[#[[1]]/(nStates+1)]}->#[[2]]&,Tally[Drop[Drop[Join[(nStates+1)*cvPart[[nn]],Table[0,{nLagtime}]]+Join[Table[0,{nLagtime}],cvPart[[nn]]],-nLagtime],nLagtime]]],{nn,1,Length[cvPart]}],{nStates,nStates}];

			]
*)
			cV=v;
			If[ClosedEnd/.List[optsCM]/.{ClosedEnd->True},
				While[(Sort[Union[Drop[cV,nLagtime]]]!=Sort[Union[Drop[cV,-nLagtime]]])&&(Length[cV]>nLagtime),
					cV=Drop[cV,-1];
					cM=SparseArray[Map[{Mod[#[[1]],(nStates+1)],Floor[#[[1]]/(nStates+1)]}->#[[2]]&,Tally[Drop[Drop[Join[(nStates+1)*cV,Table[0,{nLagtime}]]+Join[Table[0,{nLagtime}],cV],-nLagtime],nLagtime]]],{nStates,nStates}];
				];
				If[Length[cV]>=nLagtime,
					cM=SparseArray[Map[{Mod[#[[1]],(nStates+1)],Floor[#[[1]]/(nStates+1)]}->#[[2]]/nLagtime&,Tally[Drop[Drop[Join[(nStates+1)*cV,Table[0,{nLagtime}]]+Join[Table[0,{nLagtime}],cV],-nLagtime],nLagtime]]],{nStates,nStates}],
					cM=SparseArray[{},{nStates,nStates}]
				];
			];
			CountMatrix[cM,RemoveRules[ClosedEnd,optsCM]]
		]

CountMatrix/:CountMatrix[TimeSeries[v_?VectorQ,opts:OptionsPattern[TimeSeries]],optsCM:OptionsPattern[Join[Options[CountMatrix],{Lagtime->1}]]] := 
		Module[{cM,options,nStates,nLagtime},			
			options=UnionRules[opts,States->Max[v]];
			nStates=Max[OptionValue[TimeSeries,States],Max[v]];
			nLagtime=OptionValue[TimeSeries,optsCM,Lagtime];
			CountMatrix[Normal[SparseArray[Map[{Mod[#[[1]],(nStates+1)],Floor[#[[1]]/(nStates+1)]}->#[[2]]&,Tally[Drop[Drop[Join[(nStates+1)*v,Table[0,{nLagtime}]]+Join[Table[0,{nLagtime}],v],-nLagtime],nLagtime]]],{nStates,nStates}]]/nLagtime,optsCM]
		]

CountMatrix/:NormalizedQ[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]] := 
			Apply[And,Map[#==0&,Map[Apply[Plus,#]&,m]]
		]



(*
singleSymMatrix6={
		{(-r12-r13-r14-r15-r16),r12,r13,r14,r15,r16},
		{r12*p1/p2,(-r12*p1/p2-r23-r24-r25-r26),r23,r24,r25,r26},
		{r13*p1/p3,r23*p2/p3,-r13*p1/p3-r23*p2/p3-r34-r35-r36,r34,r35,r36},
		{r14*p1/p4,r24*p2/p4,r34*p3/p4,-r14*p1/p4-r24*p2/p4-r34*p3/p4-r45-r46,r45,r46},
		{r15*p1/p5,r25*p2/p5,r35*p3/p5,r45*p4/p5,-r15*p1/p5-r25*p2/p5-r35*p3/p5-r45*p4/p5-r56,r56},
		{r16*p1/(1-p1-p2-p3-p4-p5),r26*p2/(1-p1-p2-p3-p4-p5),r36*p3/(1-p1-p2-p3-p4-p5),r46*p4/(1-p1-p2-p3-p4-p5),r56*p5/(1-p1-p2-p3-p4-p5),(-r16*p1-r26*p2-r36*p3-r46*p4-r56*p5)/(1-p1-p2-p3-p4-p5)}
};
singleSymMatrixFunction6[r12_,r13_,r14_,r15_,r16_,r23_,r24_,r25_,r26_,r34_,r35_,r36_,r45_,r46_,r56_,p1_,p2_,p3_,p4_,p5_]:=
{
		{(-r12-r13-r14-r15-r16),r12,r13,r14,r15,r16},
		{r12*p1/p2,(-r12*p1/p2-r23-r24-r25-r26),r23,r24,r25,r26},
		{r13*p1/p3,r23*p2/p3,-r13*p1/p3-r23*p2/p3-r34-r35-r36,r34,r35,r36},
		{r14*p1/p4,r24*p2/p4,r34*p3/p4,-r14*p1/p4-r24*p2/p4-r34*p3/p4-r45-r46,r45,r46},
		{r15*p1/p5,r25*p2/p5,r35*p3/p5,r45*p4/p5,-r15*p1/p5-r25*p2/p5-r35*p3/p5-r45*p4/p5-r56,r56},
		{r16*p1/(1-p1-p2-p3-p4-p5),r26*p2/(1-p1-p2-p3-p4-p5),r36*p3/(1-p1-p2-p3-p4-p5),r46*p4/(1-p1-p2-p3-p4-p5),r56*p5/(1-p1-p2-p3-p4-p5),(-r16*p1-r26*p2-r36*p3-r46*p4-r56*p5)/(1-p1-p2-p3-p4-p5)}
};
*)

singleSymMatrix6={
		{-r12-r13-r14-r15-r16,r12,r13,r14,r15,r16}/p1,
		{r12,-r12-r23-r24-r25-r26,r23,r24,r25,r26}/p2,
		{r13,r23,-r13-r23-r34-r35-r36,r34,r35,r36}/p3,
		{r14,r24,r34,-r14-r24-r34-r45-r46,r45,r46}/p4,
		{r15,r25,r35,r45,-r15-r25-r35-r45-r56,r56}/p5,
		{r16,r26,r36,r46,r56,-r16-r26-r36-r46-r56}/(1-p1-p2-p3-p4-p5)
};
(*
singleSymMatrix6Out={
		{(-r12-r13-r14-r15-r16)*p1,r12*p2,r13*p3,r14*p4,r15*p5,r16*(1-p1-p2-p3-p4-p5)}/p1,
		{r12*p1,(-r12-r23-r24-r25-r26)*p2,r23*p3,r24*p4,r25*p5,r26*(1-p1-p2-p3-p4-p5)}/p2,
		{r13*p1,r23*p2,(-r13-r23-r34-r35-r36)*p3,r34*p4,r35*p5,r36*(1-p1-p2-p3-p4-p5)}/p3,
		{r14*p1,r24*p2,r34*p3,(-r14-r24-r34-r45-r46)*p4,r45*p5,r46*(1-p1-p2-p3-p4-p5)}/p4,
		{r15*p1,r25*p2,r35*p3,r45*p4,(-r15-r25-r35-r45-r56)*p5,r56*(1-p1-p2-p3-p4-p5)}/p5,
		{r16*p1,r26*p2,r36*p3,r46*p4,r56*p5,(-r16-r26-r36-r46-r56)*(1-p1-p2-p3-p4-p5)}/(1-p1-p2-p3-p4-p5)
};
*)

singleSymMatrixFunction6[r12_,r13_,r14_,r15_,r16_,r23_,r24_,r25_,r26_,r34_,r35_,r36_,r45_,r46_,r56_,p1_,p2_,p3_,p4_,p5_]:=
{
		{-r12-r13-r14-r15-r16,r12,r13,r14,r15,r16}/p1,
		{r12,-r12-r23-r24-r25-r26,r23,r24,r25,r26}/p2,
		{r13,r23,-r13-r23-r34-r35-r36,r34,r35,r36}/p3,
		{r14,r24,r34,-r14-r24-r34-r45-r46,r45,r46}/p4,
		{r15,r25,r35,r45,-r15-r25-r35-r45-r56,r56}/p5,
		{r16,r26,r36,r46,r56,-r16-r26-r36-r46-r56}/(1-p1-p2-p3-p4-p5)
};

likelihoodSymMatrix6[cm_,r12_?NumberQ,r13_,r14_,r15_,r16_,r23_,r24_,r25_,r26_,r34_,r35_,r36_,r45_,r46_,r56_,p1_,p2_,p3_,p4_,p5_] :=
		Module[{matrixN,matrixExp,matrixAbs,ret,noValue=100000.0},
			matrixN=Evaluate[singleSymMatrixFunction6[r12,r13,r14,r15,r16,r23,r24,r25,r26,r34,r35,r36,r45,r46,r56,p1,p2,p3,p4,p5]];
			matrixExp=Chop[MatrixExp[matrixN]];
			If[(Min[matrixExp]<=0.0)&&(Max[matrixExp]>1.0),
				noValue,
(*				matrixAbs=Abs[Sign[Chop[matrixN]]];
				matrixExp=matrixExp+(1-matrixAbs)*0.00000000001;*)
				ret=-Re[Apply[Plus,Flatten[Log[matrixExp]*cm]]];
    			ret=If[NumberQ[ret],ret,noValue]
			]
		]

CountMatrix/:ToHummerRateMatrix6[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]], optsHR___] :=
		Module[{erg,miR=0.0,maR=1.0,miP=0.0001},
			erg=NMinimize[
				{likelihoodSymMatrix6[m,r12,r13,r14,r15,r16,r23,r24,r25,r26,r34,r35,r36,r45,r46,r56,p1,p2,p3,p4,p5],
				r12>miR,r13>miR,r14>miR,r15>miR,r16>miR,r23>miR,r24>miR,r25>miR,r26>miR,r34>miR,r35>miR,r36>miR,r45>miR,r46>miR,r56>miR,
				r12<maR,r13<maR,r14<maR,r15<maR,r16<maR,r23<maR,r24<maR,r25<maR,r26<maR,r34<maR,r35<maR,r36<maR,r45<maR,r46<maR,r56<maR,
				p1>miP,p2>miP,p3>miP,p4>miP,p5>miP, p1+p2+p3+p4+p5<1.0-miP,p5<0.15,p1+p2+p3+p4+p5>0.85},
				{r12,r13,r14,r15,r16,r23,r24,r25,r26,r34,r35,r36,r45,r46,r56,p1,p2,p3,p4,p5},
			optsHR];
			Print[erg];
			RateMatrix[singleSymMatrix6/.erg[[2]]]
		];
TransitionMatrix/:Likelihood[TransitionMatrix[t_?MatrixQ,optsT:OptionsPattern[CountMatrix]],CountMatrix[m_?MatrixQ,optsM:OptionsPattern[CountMatrix]]] :=
		Module[{matrixAbs,t2},
			matrixAbs=Abs[Sign[m]];
			t2=t+(1-matrixAbs);
			-Re[Apply[Plus,Flatten[Log[t2]*m]]]
		];
RateMatrix/:Likelihood[RateMatrix[r_?MatrixQ,optsT:OptionsPattern[CountMatrix]],CountMatrix[m_?MatrixQ,optsM:OptionsPattern[CountMatrix]]] :=
		Module[{matrixAbs,t,t2},
			t=MatrixExp[r];
			matrixAbs=Abs[Sign[m]];
			t2=t+(1-matrixAbs);
			-Re[Apply[Plus,Flatten[Log[t2]*m]]]
		];


singleSymMatrix6Geometric={
		{(-r12-r13-r15-r16),r12,r13,0,r15,r16},
		{r12*p1/p2,(-r12*p1/p2-r24-r25-r26),0,r24,r25,r26},
		{r13*p1/p3,0,-r13*p1/p3-r34-r35-r36,r34,r35,r36},
		{0,r24*p2/p4,r34*p3/p4,-r24*p2/p4-r34*p3/p4-r45-r46,r45,r46},
		{r15*p1/p5,r25*p2/p5,r35*p3/p5,r45*p4/p5,-r15*p1/p5-r25*p2/p5-r35*p3/p5-r45*p4/p5-r56,r56},
		{r16*p1/(1-p1-p2-p3-p4-p5),r26*p2/(1-p1-p2-p3-p4-p5),r36*p3/(1-p1-p2-p3-p4-p5),r46*p4/(1-p1-p2-p3-p4-p5),r56*p5/(1-p1-p2-p3-p4-p5),(-r16*p1-r26*p2-r36*p3-r46*p4-r56*p5)/(1-p1-p2-p3-p4-p5)}
};
singleSymMatrixFunction6Geometric[r12_,r13_,r14_,r15_,r16_,r23_,r24_,r25_,r26_,r34_,r35_,r36_,r45_,r46_,r56_,p1_,p2_,p3_,p4_,p5_]:=
{
		{(-r12-r13-r15-r16),r12,r13,0,r15,r16},
		{r12*p1/p2,(-r12*p1/p2-r24-r25-r26),0,r24,r25,r26},
		{r13*p1/p3,0,-r13*p1/p3-r34-r35-r36,r34,r35,r36},
		{0,r24*p2/p4,r34*p3/p4,-r24*p2/p4-r34*p3/p4-r45-r46,r45,r46},
		{r15*p1/p5,r25*p2/p5,r35*p3/p5,r45*p4/p5,-r15*p1/p5-r25*p2/p5-r35*p3/p5-r45*p4/p5-r56,r56},
		{r16*p1/(1-p1-p2-p3-p4-p5),r26*p2/(1-p1-p2-p3-p4-p5),r36*p3/(1-p1-p2-p3-p4-p5),r46*p4/(1-p1-p2-p3-p4-p5),r56*p5/(1-p1-p2-p3-p4-p5),(-r16*p1-r26*p2-r36*p3-r46*p4-r56*p5)/(1-p1-p2-p3-p4-p5)}
};
singleSymMatrix6Geometric={
		{(-r12-r13-r15-r16)/p1,r12/p2,r13/p3,0,r15/p5,r16/(1-p1-p2-p3-p4-p5)},
		{r12/p1,(-r12-r24-r25-r26)/p2,0,r24/p4,r25/p5,r26/(1-p1-p2-p3-p4-p5)},
		{r13/p1,0,(-r13-r34-r35-r36)/p3,r34/p4,r35/p5,r36/(1-p1-p2-p3-p4-p5)},
		{0,r24/p2,r34/p3,(-r24-r34-r45-r46)/p4,r45/p5,r46/(1-p1-p2-p3-p4-p5)},
		{r15/p1,r25/p2,r35/p3,r45/p4,(-r15-r25-r35-r45-r56)/p5,r56/(1-p1-p2-p3-p4-p5)},
		{r16/p1,r26/p2,r36/p3,r46/p4,r56/p5,(-r16-r26-r36-r46-r56)/(1-p1-p2-p3-p4-p5)}
};
singleSymMatrixFunction6Geometric[r12_,r13_,r15_,r16_,r24_,r25_,r26_,r34_,r35_,r36_,r45_,r46_,r56_,p1_,p2_,p3_,p4_,p5_]:={{(-r12-r13-r15-r16)/p1,r12/p2,r13/p3,0/p4,r15/p5,r16/(1-p1-p2-p3-p4-p5)},{r12/p1,(-r12-r24-r25-r26)/p2,0,r24/p4,r25/p5,r26/(1-p1-p2-p3-p4-p5)},{r13/p1,0,(-r13-r34-r35-r36)/p3,r34/p4,r35/p5,r36/(1-p1-p2-p3-p4-p5)},{0,r24/p2,r34/p3,(-r24-r34-r45-r46)/p4,r45/p5,r46/(1-p1-p2-p3-p4-p5)},{r15/p1,r25/p2,r35/p3,r45/p4,(-r15-r25-r35-r45-r56)/p5,r56/(1-p1-p2-p3-p4-p5)},{r16/p1,r26/p2,r36/p3,r46/p4,r56/p5,(-r16-r26-r36-r46-r56)/(1-p1-p2-p3-p4-p5)}};
likelihoodSymMatrix6Geometric[cm_,r12_?NumberQ,r13_,r15_,r16_,r24_,r25_,r26_,r34_,r35_,r36_,r45_,r46_,r56_,p1_,p2_,p3_,p4_,p5_]:=Module[{matrixN,matrixExp,matrixAbs},
		matrixN=Evaluate[singleSymMatrixFunction6Geometric[r12,r13,r15,r16,r24,r25,r26,r34,r35,r36,r45,r46,r56,p1,p2,p3,p4,p5]];
		matrixExp=MatrixExp[matrixN];
		matrixAbs=Abs[Sign[matrixN]];
		matrixExp=matrixExp+(1-matrixAbs)*0.00001;
		-Re[Apply[Plus,Flatten[Log[matrixExp]*cm]]]
]

CountMatrix/:ToHummerRateMatrix6Geometric[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]], optsHR___] :=
		Module[{erg},
			erg=NMinimize[
				{likelihoodSymMatrix6Geometric[m,r12,r13,r15,r16,r24,r25,r26,r34,r35,r36,r45,r46,r56,p1,p2,p3,p4,p5],
				r12>0.001,r13>0.001,r15>0.001,r16>0.001,r24>0.001,r25>0.001,r26>0.001,r34>0.001,r35>0.001,r36>0.001,r45>0.001,r46>0.001,r56>0.001,
				r12<2,r13<2,r15<2,r16<2,r24<2,r25<2,r26<2,r34<2,r35<2,r36<2,r45<2,r46<2,r56<2,
				p1>0.001,p2>0.001,p3>0.001,p4>0.001,p5>0.001,
				p1<0.999,p2<0.999,p3<0.999,p4<0.999,p5<0.999,p1+p2+p3+p4+p5<0.999,
				p1},
				{r12,r13,r15,r16,r24,r25,r26,r34,r35,r36,r45,r46,r56,p1,p2,p3,p4,p5},
			optsHR];
			RateMatrix[singleSymMatrix6Geometric/.erg[[2]]]
		];


(* RateMatrix functions *)

RateMatrix/:MakeBoxes[RateMatrix[m_?MatrixQ,opts:OptionsPattern[RateMatrix]],StandardForm] := 
		If[OptionValue[RateMatrix,ShowSize]>Length[m],
			InterpretationBox[#,RateMatrix[m,opts]]&[RowBox[{SubscriptBox["\[DoubleStruckCapitalR]",Length[m]],ToBoxes[MatrixForm[m]]}]],
			InterpretationBox[#,RateMatrix[m,opts]]&[SubscriptBox["\[DoubleStruckCapitalR]",Length[m]]]
		]
RateMatrix/:RateMatrix[RateMatrix[m_?MatrixQ,opts1:OptionsPattern[RateMatrix]],opts2:OptionsPattern[RateMatrix]] := 
		RateMatrix[m,UnionRules[opts1,opts2]]
RateMatrixQ[x_] := 
		If[Head[x]==RateMatrix&&MatrixQ[x[[1]]]&&RowSumZeroQ[m],True,False];
RateMatrix/:NormalForm[RateMatrix[m_?MatrixQ,opts:OptionsPattern[RateMatrix]]] := 
		{m,opts}
RateMatrix/:Normal[RateMatrix[m_?MatrixQ,opts:OptionsPattern[RateMatrix]]] := 
		m
RateMatrix/:Spectrum[RateMatrix[m_?MatrixQ,opts:OptionsPattern[RateMatrix]]] := 
		Eigenvalues[N[m]]
RateMatrix/:MatrixForm[RateMatrix[m_?MatrixQ,opts:OptionsPattern[RateMatrix]]] := 
		MatrixForm[m]
RateMatrix/:Length[RateMatrix[m_?MatrixQ],opts:OptionsPattern[RateMatrix]] := 
		Length[m]
RateMatrix/:RowSumZeroQ[RateMatrix[m_?MatrixQ],opts:OptionsPattern[RateMatrix]] := 
		RowSumZeroQ[m]
RateMatrix/:MatrixPlot[RateMatrix[m_?MatrixQ,opts:OptionsPattern[RateMatrix]],optsMP:OptionsPattern[MatrixPlot]] := 
		MatrixPlot[m,optsMP,Frame->True,BaseStyle->{FontFamily->"GillSans",11},FrameTicks->{Automatic,Automatic,None,None},ImageSize->320,ColorFunction->"TemperatureMap"]
RateMatrix/:ArrayPlot[RateMatrix[m_?MatrixQ,opts:OptionsPattern[RateMatrix],optsAP:OptionsPattern[ArrayPlot]]] := 
		ArrayPlot[m,optsAP,Frame->True,FrameTicks->All,DataReversed->True,BaseStyle->{FontFamily->"GillSans",FontSize->12}]
RateMatrix/:N[RateMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]] := 
		RateMatrix[N[m],opts]
RateMatrixQ[x_] := 
		If[Head[x]==RateMatrix&&MatrixQ[x[[1]]],True,False];
RateMatrix/:NormalizedQ[RateMatrix[m_?MatrixQ,opts:OptionsPattern[RateMatrix]]] := 
		RowSumZeroQ[m]
RateMatrix/:Normalize[RateMatrix[m_?MatrixQ,opts:OptionsPattern[RateMatrix]]] := 
		RateMatrix[Map[#-Apply[Plus,#]&,m],ReplaceRules[Normalized->True,opts]]
RateMatrixNormalized[m_?MatrixQ,opts:OptionsPattern[RateMatrix]] := 
		RateMatrix[Map[#-Apply[Plus,#]&,m],ReplaceRules[Normalized->True,opts]]
RateMatrix/:Equilibrium[RateMatrix[m_?MatrixQ,opts:OptionsPattern[RateMatrix]]] := 
		#/Apply[Plus,#]&[Eigenvectors[Transpose[m]][[-1]]]
RateMatrix/:ToTransitionMatrix[RateMatrix[m_?MatrixQ,opts:OptionsPattern[RateMatrix]],nLagtime_] := 
		TransitionMatrix[MatrixExp[m*nLagtime],Lagtime->nLagtime]


(* TransitionMatrix functions *)

TransitionMatrix/:MakeBoxes[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],StandardForm] := 
		If[OptionValue[TransitionMatrix,ShowSize]>Length[m],
			InterpretationBox[#,TransitionMatrix[m,opts]]&[RowBox[{SubscriptBox["\[DoubleStruckCapitalT]",Length[m]],ToBoxes[MatrixForm[m]]}]],
			InterpretationBox[#,TransitionMatrix[m,opts]]&[SubscriptBox["\[DoubleStruckCapitalT]",Length[m]]]
		]
TransitionMatrix/:TransitionMatrix[TransitionMatrix[m_?MatrixQ,opts1:OptionsPattern[TransitionMatrix]],opts2:OptionsPattern[TransitionMatrix]] := 
		TransitionMatrix[m,UnionRules[opts1,opts2]]
TransitionMatrixFromMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]] := 
		TransitionMatrix[Map[#/Apply[Plus,#]&,m]]
TransitionMatrix/:NormalForm[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]] := 
		{m,opts}
TransitionMatrix/:Normal[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]] := 
		m
TransitionMatrixNormalized[m_?MatrixQ, opts:OptionsPattern[TransitionMatrix]] := 
		TransitionMatrix[Map[#/Apply[Plus,#]&,m], opts]
TransitionMatrix/:Spectrum[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]] := 
		Eigenvalues[N[m]]
TransitionMatrix/:MatrixForm[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]] := 
		MatrixForm[m]
TransitionMatrix/:Length[TransitionMatrix[m_?MatrixQ],opts:OptionsPattern[TransitionMatrix]] := 
		Length[m]
TransitionMatrix/:RowSumOneQ[TransitionMatrix[m_?MatrixQ],opts:OptionsPattern[TransitionMatrix]] := 
		RowSumOneQ[m]
TransitionMatrix/:MatrixPlot[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]] := 
		MatrixPlot[m,Frame->True,BaseStyle->{FontFamily->"GillSans",11},FrameTicks->{Automatic,Automatic,None,None},ImageSize->320,ColorFunction->"TemperatureMap"]
TransitionMatrix/:ArrayPlot[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix],optsAP:OptionsPattern[ArrayPlot]]] := 
		ArrayPlot[m,optsAP,Frame->True,FrameTicks->All,DataReversed->True,BaseStyle->{FontFamily->"GillSans",FontSize->12}]
TransitionMatrixQ[x_] := 
		If[Head[x]===TransitionMatrix&&MatrixQ[x[[1]]]&&RowSumOneQ[m],True,False];
TransitionMatrix/:Spectrum[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]] := 
		Eigenvalues[N[m]]
TransitionMatrix/:Chop[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]] := 
		TransitionMatrix[Chop[m],opts]
TransitionMatrix/:NormalizedQ[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]] := 
		OptionValue[TransitionMatrix,Normalized]
TransitionMatrix/:Normalize[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]] := 
		TransitionMatrix[Map[#/Apply[Plus,#]&,m],ReplaceRules[Normalized->True,opts]]
TransitionMatrix/:Committor[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],states_] := 
		Module[{cM,cMEV},
			cM=m;
			Map[(cM[[#]]=Table[If[ii==#,1,0],{ii,1,Length[cM]}])&,states];
			cMEV=Eigenvectors[cM,Min[Length[states]+50,Length[m]]][[Range[Length[states]]]];
			Inverse[cMEV[[All,states]]].cMEV
		]
TransitionMatrix/:PotentialOfMeanForce[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],states_] := 
		Module[{committor, equilibrium, s2,s3,dd},
			committor=Committor[TransitionMatrix[m,opts],states][[2]];
			equilibrium=Equilibrium[TransitionMatrix[m,opts]];
			s2=Sort[Transpose[{committor,equilibrium}],#1[[1]]<#2[[1]]&];
			dd=0;
			s3=Transpose[{s2[[All,1]],Table[dd=dd+s2[[ii,2]],{ii,1,Length[s2]}]}];
			Transpose[{s3[[All,1]],Identity[s3[[All,2]]]}]
		]
TransitionMatrix/:PotentialOfMeanForceBins[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],states_,bins_] := 
		Module[{committor, equilibrium, s2,s3,dd,binning},
			committor=Committor[TransitionMatrix[m,opts],states][[2]];
			equilibrium=Equilibrium[TransitionMatrix[m,opts]];
			s2=Sort[Transpose[{committor,equilibrium}],#1[[1]]<#2[[1]]&];
			Select[Table[{(bins[[nn]]+bins[[nn+1]])/2,-Log[Apply[Plus,Select[s2,(#[[1]]>=bins[[nn]]&&#[[1]]<=bins[[nn+1]])&][[All,2]]]/(bins[[nn+1]]-bins[[nn]])]},{nn,1,Length[bins]-1}],#[[2]]<\[Infinity]&]
		]
TransitionMatrix/:CommittorReactionCoordinate[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],states_,expectations_] := 
		Module[{committor, equilibrium, s2,s3,dd},
			committor=Committor[TransitionMatrix[m,opts],states][[2]];
			s2=Sort[Transpose[{committor,expectations}],#1[[1]]<#2[[1]]&];
			dd=0;
			s3=Transpose[{s2[[All,1]],Table[dd=dd+s2[[ii,2]],{ii,1,Length[s2]}]}];
			Transpose[{s3[[All,1]],s3[[All,2]]/Apply[Plus,expectations]}]
		]
TransitionMatrix/:Equilibrium[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]] := 
		#/Apply[Plus,#]&[Eigenvectors[Transpose[m],Min[20,Length[m]]][[1]]]
TransitionMatrix/:ToRateMatrix[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],optsTRM:OptionsPattern[ToRateMatrix]] := 
		Module[{rM,eVecs,eVals,nLagtime},
			If[OptionValue[ToRateMatrix,optsTRM,Method]===MatrixLog,
				{eVals,eVecs} = Eigensystem[m];
				nLagtime = If[OptionValue[ToRateMatrix,optsTRM,Lagtime]===Inherit,
					OptionValue[TransitionMatrix,Lagtime],
					OptionValue[ToRateMatrix,optsTRM,Lagtime]
				];
				rM=RateMatrix[Transpose[eVecs].DiagonalMatrix[Log[eVals]/nLagtime].Inverse[Transpose[eVecs]],Lagtime->nLagtime];
			];
			rM
		]
TransitionMatrix/:ImpliedTimeScales[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],optsITS:OptionsPattern[{Lagtime->1}]] := 
		Module[{nLagtime},
			nLagtime=Lagtime/.List[UnionRules[optsITS,opts,Lagtime->1]];
			Join[{Infinity},-nLagtime/Log[Abs[Eigenvalues[N[m]][[Range[2,Length[m]]]]]]]
		]
RateMatrix/:ImpliedTimeScales[RateMatrix[m_?MatrixQ,opts:OptionsPattern[RateMatrix]]] := 
		Module[{},
			Join[{Infinity},Reverse[-1.0/Eigenvalues[N[m]][[Range[1,Length[m]-1]]]]]
		]
TransitionMatrix/:MeanFirstPassageTime[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],state_,optsMFPT:OptionsPattern[MeanFirstPassageTime]] := 
		Module[{cM,cMEV,nLagtime},
			nLagtime=Lagtime/.List[UnionRules[optsMFPT,opts,Lagtime->1]];
			cM=m-IdentityMatrix[Length[m]];
			cM[[state]]=Table[If[ii==state,1,0],{ii,1,Length[cM]}];
			Inverse[cM].Table[If[ii==state,0,-1],{ii,1,Length[cM]}]*nLagtime
		]
TransitionMatrix/:MeanFirstPassageTimeSensitivityTensor[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],state_,optsMFPT:OptionsPattern[MeanFirstPassageTime]] := 
		Module[{cM,cMEV,nLagtime,mfpt,cMInv,size},
			size=Length[m];
			nLagtime=Lagtime/.List[UnionRules[optsMFPT,opts,Lagtime->1]];
			cM=m-IdentityMatrix[Length[m]];
			cM[[state]]=Table[If[ii==state,1,0],{ii,1,Length[cM]}];
			cMInv=Inverse[cM];
			mfpt=cMInv.Table[If[ii==state,0,-1],{ii,1,Length[cM]}]*nLagtime;
			Table[If[ii==state,0,1]*cMInv[[kk,ii]]*mfpt[[jj]],{ii,1,size},{jj,1,size},{kk,1,size}]
		]
TransitionMatrix/:PotentialOfMeanForceError[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],state1_,state2_] := 
		Module[{cM,cP,size,cC,cL,cB},
			size=Length[m];
			cM=PseudoInverse[Join[Transpose[m]-IdentityMatrix[size],{Table[1,{size}]}]];
			cP=#/Apply[Plus,#]&[Eigenvectors[Transpose[m],1][[1]]];
			cC=Committor[TransitionMatrix[m,opts],{state1,state2}][[2]];
			cL=Sort[Transpose[{Range[Length[m]],cC}],#1[[2]]<#2[[2]]&];
			Table[
				cB=Map[Apply[Plus,#]&,cM[[All,cL[[Range[nn],1]]]]];
				{cL[[nn,2]],Apply[Plus,Abs[Flatten[Outer[Plus,cB,-cB]]]]}
			,{nn,1,10}]
		]
TransitionMatrix/:PotentialOfMeanForceSensitivityMatrix[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],state1_,state2_,c_] := 
		Module[{cM,cP,size,cC,cL,cB},
			size=Length[m];
			cM=PseudoInverse[Join[Transpose[m]-IdentityMatrix[size],{Table[1,{size}]}]];
			cP=#/Apply[Plus,#]&[Eigenvectors[Transpose[m],1][[1]]];
			cC=Committor[TransitionMatrix[m,opts],{state1,state2}][[2]];
			cL=Select[Transpose[{Range[Length[m]],cC}],(#[[2]]<=c)&][[All,1]];
			cB=Map[Apply[Plus,#]&,cM[[All,cL]]];
			Outer[Plus,cB,-cB]
			(*Table[cP[[aa]]*(cB[[aa]]-cB[[bb]]),{aa,1,size},{bb,1,size}]*)
		]
TransitionMatrix/:EquilibriumSensitivityTensor[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]] := 
		Module[{cM,cP,size},
			size=Length[m];
			cM=PseudoInverse[Join[Transpose[m]-IdentityMatrix[size],{Table[1,{size}]}]];
			cP=#/Apply[Plus,#]&[Eigenvectors[Transpose[m],1][[1]]];
			Table[cP[[aa]]*cM.Table[If[ii==aa,1,0]+If[ii==bb,-1,0],{ii,1,size+1}],{aa,1,size},{bb,1,size}]
		]
TransitionMatrix/:CommittorSensitivityMatrix[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],states_] := 
		Module[{cM,cS,cA,size},
			size=Length[m];
			cS=IdentityMatrix[size][[Complement[Range[size],states]]];
			cM=Inverse[(m-IdentityMatrix[size])[[Complement[Range[size],states],Complement[Range[size],states]]]];
			cA=Transpose[cS].cM.cS;
			Transpose[cA]
		]
TransitionMatrix/:ToTimeSeries[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],state_,num___:1,optsTTS:OptionsPattern[ToTimeSeries]] := 
		Module[{probabilities,probSum,endState,p,size,list,startState,listOfRules},
			If[OptionValue[ToTimeSeries,optsTTS,AssumeSparse]==True,
				listOfRules=Table[
					(#[[All,1]]->#[[All,2]])&[Map[#[[2]]->#[[1,1]]&,Drop[ArrayRules[m[[nn]]],-1]]]
				,{nn,1,Length[m]}];
				list = NestList[RandomChoice[listOfRules[[#]]]&,state,num];
			,
				size=Length[m];
				list=Table[0,{num+1}];
				endState=state;
				list[[1]]=state;
				Table[
					startState=endState;
					probSum=m[[startState,1]];
					p=RandomReal[];
					endState=1;
					While[p>probSum,
						endState++;
						probSum+=m[[startState,endState]];
					];
					endState=Min[endState,size];
					list[[nn]]=endState;
				,{nn,2,num+1}];
			];
			TimeSeries[list]
		]
TransitionMatrix/:ChooseTransition[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],state_,num___:1] := 
		Module[{},
			ToTimeSeries[TransitionMatrix[m,opts],state,num]
		]
TransitionMatrix/:Symmetrize[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]] := 
		Module[{eqDist,size,cM},
			size=Length[m];
			eqDist=(#/Apply[Plus,#]&[Eigenvectors[Transpose[m],1][[1]]]);
			cM=m;
			Table[cM[[ii]]=cM[[ii]]*eqDist[[ii]];,{ii,1,size}];
			(cM+Transpose[cM])/2
		]
TransitionMatrix/:NormalizedQ[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]] := 
		Apply[And,Map[#==1&,Map[Apply[Plus,#]&,m]]]
TransitionMatrix/:VisitPredictionMatrix[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],length_,optsMFPT:OptionsPattern[MeanFirstPassageTime]] := 
		Module[{evals,evecs,nLagtime,nLength},
			nLagtime=Lagtime/.List[UnionRules[optsMFPT,opts,Lagtime->1]];
			nLength=length / nLagtime;
			{evals,evecs} = Eigensystem[m];
			Transpose[evecs].DiagonalMatrix[Map[If[#<1,(#^(nLength+1)-1)/(#-1),(nLength+1)]&,evals]].Transpose[Inverse[evecs]]
		]


(* Distribution Implementations*)

DirichletDistribution/:MeanMatrix[DirichletDistribution[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]] := 
		TransitionMatrixFromMatrix[Map[#/Apply[Plus,#]&,(m+OptionValue[CountMatrix,Prior])]]
DirichletDistribution/:MaximumLikelihoodMatrix[DirichletDistribution[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]] := 
		TransitionMatrixFromMatrix[Map[#/Apply[Plus,#]&,(m)]]
DirichletDistribution/:VarianceMatrix[DirichletDistribution[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]] := 
		Map[#*(Apply[Plus,#]-#)/(Apply[Plus,#]^2*(Apply[Plus,#]+1))&,(m+OptionValue[CountMatrix,Prior])]
DirichletDistribution/:DrawTransitionMatrix[DirichletDistribution[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]] := 
		TransitionMatrix[Map[#/Apply[Plus,#]&,Map[RandomReal[GammaDistribution[#,1],1][[1]]&,m+1,{2}]]]
DiracDistribution/:MeanMatrix[DiracDistribution[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]] := 
		TransitionMatrixFromMatrix[Map[#/Apply[Plus,#]&,(m)]]
DiracDistribution/:MaximumLikelihoodMatrix[DiracDistribution[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]] := 
		TransitionMatrixFromMatrix[Map[#/Apply[Plus,#]&,(m)]]
DiracDistribution/:VarianceMatrix[DiracDistribution[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]] := 
		0*m
DiracDistribution/:DrawTransitionMatrix[DiracDistribution[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]] := 
		TransitionMatrix[Map[#/Apply[Plus,#]&,m]]
DirichletWithDBDistribution/:MeanMatrix[DirichletWithDBDistribution[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]] := 
		Message[MeanMatrix::NYI]
DirichletWithDBDistribution/:VarianceMatrix[DirichletWithDBDistribution[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]] := 
		Message[VarianceMatrix::NYI]
DirichletWithDBDistribution/:DrawTransitionMatrix[DirichletWithDBDistribution[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]] := 
		Message[DrawTransitionMatrix::NYI]
DirichletWithDBDistribution/:MaximumLikelihoodMatrix[DirichletWithDBDistribution[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]],optsML:OptionsPattern[DirichletWithDBDistribution]] := 
		Module[{size,sumVec,symmMat,result,nIterations,sumResVec,nResidue,oRes,ii},
			size=Length[m];
			nIterations=OptionValue[DirichletWithDBDistribution,List[optsML],Iterations];
			sumVec=Map[Apply[Plus,#]&,m];
			symmMat=m+Transpose[m];
			result=Table[1./size,{size},{size}];
			nResidue=OptionValue[DirichletWithDBDistribution,List[optsML],MaximalResidue];
			For[ii=1,ii<=nIterations,ii++,
				sumResVec=Map[Apply[Plus,#]&,result];
				oRes=result;
				result=Table[symmMat[[ii,jj]]/(sumVec[[ii]]/sumResVec[[ii]]+sumVec[[jj]]/sumResVec[[jj]]),{ii,1,size},{jj,1,size}];
				If[Max[Abs[result-oRes]]<nResidue,ii=nIterations];
			];
			TransitionMatrix[Map[#/Apply[Plus,#]&,result]]
		]


(* TimeSeries & Trajectories *)

TimeSeries/:MakeBoxes[TimeSeries[v_?VectorQ,opts:OptionsPattern[TimeSeries]],StandardForm] := 
		If[Length[v]<OptionValue[TimeSeries,ShowSize],
			InterpretationBox[#,TimeSeries[v,opts]]&[RowBox[{SubscriptBox["\[DoubleStruckCapitalS]",RowBox[{Max[OptionValue[TimeSeries,States],Max[v]]}]],"(",ToBoxes[TableForm[v,TableDirections->Row]],")"}]],
			InterpretationBox[#,TimeSeries[v,opts]]&[RowBox[{SubscriptBox["\[DoubleStruckCapitalS]",RowBox[{Max[OptionValue[TimeSeries,States],Max[v]]}]],"\[LeftAngleBracket]",Length[v],"\[RightAngleBracket]"}]]
		]
TimeSeriesQ[x_] := If[Head[x]==TimeSeries&&VectorQ[x[[1]]],True,False];
TimeSeries/:NormalForm[TimeSeries[v_?VectorQ,opts:OptionsPattern[TimeSeries]]] := 
		{v,opts}
TimeSeries/:Normal[TimeSeries[v_?VectorQ,opts:OptionsPattern[TimeSeries]]] := 
		v
TimeSeries/:MatrixForm[TimeSeries[v_?VectorQ,opts:OptionsPattern[TimeSeries]]] := 
		MatrixForm[v]
TimeSeries/:Length[TimeSeries[v_?VectorQ,opts:OptionsPattern[TimeSeries]]] := 
		Length[v]
TimeSeries/:States[TimeSeries[v_?VectorQ,opts:OptionsPattern[TimeSeries]]] := 
		Max[OptionValue[TimeSeries,States],Max[v]]
TimeSeries/:MatrixPlot[TimeSeries[v_?VectorQ,opts:OptionsPattern[TimeSeries]],optsMP:OptionsPattern[MatrixPlot]] := 
		MatrixPlot[{v},optsMP,Frame->True,BaseStyle->{FontFamily->"GillSans",11},FrameTicks->{None,All,None,None},ImageSize->320,ColorFunction->"TemperatureMap"]


(* TimeSeries functions*)

Trajectory/:MakeBoxes[Trajectory[m_?MatrixQ,opts:OptionsPattern[Trajectory]],StandardForm] := 
		InterpretationBox[#,Trajectory[m,opts]]&[RowBox[{SubscriptBox["\[DoubleStruckCapitalM]",Dimensions[m][[2]]],"\[LeftAngleBracket]",Dimensions[m][[1]],"\[RightAngleBracket]"}]]
Trajectory/:ToTimeSeries[Trajectory[m_?MatrixQ,opts:OptionsPattern[Trajectory]],centers_] := 
		Module[{nearest},
			nearest=Nearest[centers->Automatic];
			TimeSeries[Map[nearest[#,1][[1]]&,m],States->Length[centers]]
		]
TrajectoryQ[x_] := If[Head[x]==Trajectory&&VectorQ[x[[1]]],True,False];
Trajectory/:NormalForm[Trajectory[v_?VectorQ,opts:OptionsPattern[Trajectory]]] := 
		{v,opts}
Trajectory/:Normal[Trajectory[v_?VectorQ,opts:OptionsPattern[Trajectory]]] := 
		v
Trajectory/:MatrixForm[Trajectory[v_?VectorQ,opts:OptionsPattern[Trajectory]]] := 
		MatrixForm[v]
Trajectory/:Length[Trajectory[v_?VectorQ],opts:OptionsPattern[Trajectory]] := 
		Length[v]
Trajectory/:KMeansClustering[Trajectory[m_?MatrixQ,opts:OptionsPattern[Trajectory]],nCenters_,optsKM:OptionsPattern[KMeansClustering]] := 
		Module[{centerList,nearest,clusters,nClusters,nIterations,oCL,residuum,ii},
			centerList = If[
				NumberQ[nCenters],
				m[[Round[Length[m]/nCenters*Range[0.5,nCenters-0.5,1.0]]]],
				nCenters
			];
			nIterations = OptionValue[KMeansClustering,optsKM,Iterations];
			For[ii=1,ii<=nIterations,ii++,
				nClusters=Length[centerList];
				nearest=Nearest[centerList->Automatic];
				clusters=Map[nearest[#,1][[1]]&,m];
				oCL=centerList;
				centerList=Table[
					Mean[m[[Select[Transpose[{Range[Length[clusters]],clusters}],#[[2]]==nn&][[All,1]]]]]
				,{nn,1,nClusters}];
				residuum=Apply[Plus,Map[Norm,centerList-oCL]];
				If[residuum==0, ii=nIterations];
			];
			centerList
		]


(* ClusterMatrix functions *)

ClusterMatrix/:MakeBoxes[ClusterMatrix[m_?MatrixQ,opts:OptionsPattern[ClusterMatrix]],StandardForm] := 
		InterpretationBox[#,ClusterMatrix[m,opts]]&[SubscriptBox["\[DoubleStruckCapitalL]",RowBox[{Dimensions[m][[2]],"\[Rule]",Dimensions[m][[1]]}]]]
ClusterMatrix/:Normal[ClusterMatrix[m_?MatrixQ,opts:OptionsPattern[ClusterMatrix]]] := 
		m
ClusterMatrixQ[x_] := If[Head[x]==ClusterMatrix&&MatrixQ[x[[1]]],True,False];
ClusterMatrix/:NormalForm[ClusterMatrix[m_?MatrixQ,opts:OptionsPattern[ClusterMatrix]]] := 
		{m,opts}
ClusterMatrix/:Normal[ClusterMatrix[m_?MatrixQ,opts:OptionsPattern[ClusterMatrix]]] := 
		m
ClusterMatrix/:MatrixForm[ClusterMatrix[m_?MatrixQ,opts:OptionsPattern[ClusterMatrix]]] := 
		MatrixForm[m]
ClusterMatrix/:Length[ClusterMatrix[m_?MatrixQ],opts:OptionsPattern[ClusterMatrix]] := 
		Reverse[Dimensions[m]]
ClusterMatrix/:MatrixPlot[ClusterMatrix[m_?MatrixQ,opts:OptionsPattern[ClusterMatrix]],optsMP:OptionsPattern[MatrixPlot]] := 
		MatrixPlot[m,optsMP,Frame->True,BaseStyle->{FontFamily->"GillSans",11},FrameTicks->{Automatic,Automatic,None,None},ImageSize->320,ColorFunction->"TemperatureMap"]
ClusterMatrix/:InvolvedStates[ClusterMatrix[m_?MatrixQ,opts:OptionsPattern[ClusterMatrix]]] :=
		Drop[ArrayRules[Transpose[Normal[m]]][[All,1,1]],-1]

CrispClusterMatrix[rules_,opts:OptionsPattern[CrispClusterMatrix]] := 
		Module[{nStates,thisRules},
			nStates=Max[OptionValue[CrispClusterMatrix,States],Max[Flatten[rules]]];
			thisRules=If[OptionValue[Complete],Join[rules,Map[{#}&,Complement[Range[nStates],
				Flatten[rules]]]],
				rules
			];
			ClusterMatrix[Map[Table[0,{nStates}]+Apply[Plus,IdentityMatrix[nStates][[Flatten[{#}]]]]&,thisRules]]
		];
RemoveClusterMatrix[states_,size_] := 
		ClusterMatrix[IdentityMatrix[size][[Complement[Range[size],states]]]]
CountMatrix/:RemoveEmptyStatesClusterMatrix[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]] := 
		Module[{cM,mOld,last,emptyLine},
			mOld=Normal[m];
			last=Length[m]+1;
			cM=Range[Length[m]];
			emptyLine=Table[0,{Length[m]}];
			While[Length[cM]<last,
				last=Length[cM];
				Map[(mOld[[#,All]]=emptyLine;mOld[[All,#]]=emptyLine)&,Complement[Range[Length[m]],cM]];
				cM=Select[Table[{nn,Apply[Plus,mOld[[nn]]],Apply[Plus,mOld[[All,nn]]]},{nn,1,Length[m]}],#[[2]]>0&&#[[3]]>0&][[All,1]];				
			];
			ClusterMatrix[SparseArray[MapIndexed[{#2[[1]],#}->1&,cM],{Length[cM],Length[m]}]]
		]

TransitionMatrix/:RemoveEmptyStatesClusterMatrix[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]] := 
		Module[{cM},
			cM=Select[Table[{nn,Apply[Plus,m[[nn]]]},{nn,1,Length[m]}],#[[2]]>0&][[All,1]];
			ClusterMatrix[SparseArray[MapIndexed[{#2[[1]],#}->1&,cM],{Length[cM],Length[m]}]]
		]

ClusterMatrix/:ToCrispRules[ClusterMatrix[m_?MatrixQ,opts:OptionsPattern[ClusterMatrix]]] := 
		Module[{cM},
			Table[nn->#&[Select[Transpose[{Range[Length[m]],m[[All,nn]]}],#[[2]]==Max[m[[All,nn]]]&][[1,1]]],{nn,1,Length[m[[1]]]}]
		]
ClusterMatrix/:Crispify[ClusterMatrix[m_?MatrixQ,opts:OptionsPattern[ClusterMatrix]]] := 
		Module[{cM},
			ClusterMatrix[Normal[SparseArray[Table[{#,nn}->1&[Select[Transpose[{Range[Length[m]],m[[All,nn]]}],#[[2]]==Max[m[[All,nn]]]&][[1,1]]],{nn,1,Length[m[[1]]]}]]]]
		]
TransitionMatrix/:Cluster[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],ClusterMatrix[c_?MatrixQ,optsCM:OptionsPattern[ClusterMatrix]]] := 
		Module[{eqDist,size,cM},
			size=Length[m];
			eqDist=(#/Apply[Plus,#]&[Eigenvectors[Transpose[m],1][[1]]]);
			cM=m;
			Table[cM[[ii]]=cM[[ii]]*eqDist[[ii]];,{ii,1,size}];
			TransitionMatrix[Map[If[Apply[Plus,#]>0,#/Apply[Plus,#],Table[0,{Length[c]}]]&,c.cM.Transpose[c]],RemoveRules[States,opts]]
		];
RateMatrix/:Cluster[RateMatrix[m_?MatrixQ,opts:OptionsPattern[RateMatrix]],ClusterMatrix[c_?MatrixQ,optsCM:OptionsPattern[ClusterMatrix]]] := 
		Module[{eqDist,size,cM},
			size=Length[m];
			eqDist=(#/Apply[Plus,#]&[Eigenvectors[Transpose[m],1][[1]]]);
			cM=m;
			Table[cM[[ii]]=cM[[ii]]*eqDist[[ii]];,{ii,1,size}];
			RateMatrix[Map[If[Apply[Plus,#]>0,#/Apply[Plus,#],Table[0,{Length[c]}]]&,c.cM.Transpose[c]],RemoveRules[States,opts]]
		];
ClusterMatrix/:Cluster[v_?VectorQ,ClusterMatrix[c_?MatrixQ,optsCM:OptionsPattern[ClusterMatrix]]] := 
		Module[{rules},
			rules=ArrayRules[c];
			Table[Map[{#[[1,1]],#[[2]]}&,Select[rules,#[[1,2]]==v[[nn]]&]],{nn,1,Length[v]}]
		];
CountMatrix/:Cluster[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]],ClusterMatrix[c_?MatrixQ,optsCM:OptionsPattern[ClusterMatrix]]] := 
		CountMatrix[c.m.Transpose[c],RemoveRules[States,opts]]
TimeSeries/:Cluster[TimeSeries[v_?VectorQ,opts:OptionsPattern[TimeSeries]],ClusterMatrix[c_?MatrixQ,optsCM:OptionsPattern[ClusterMatrix]]] := 
		Module[{size,cM,rules},
			rules=ToCrispRules[ClusterMatrix[c,optsCM]];
			TimeSeries[v/.rules]
		];
CountMatrix/:HasEmptyStatesQ[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]] := 
		If[Apply[Times,Map[Apply[Plus,#]&,m]]==0,True,False]
TransitionMatrix/:HasEmptyStatesQ[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]] := 
		If[Apply[Times,Map[Apply[Plus,#]&,m]]==0,True,False]
TimeSeries/:RemoveVoid[TimeSeries[v_?VectorQ,opts:OptionsPattern[TimeSeries]],states_] := 
		Module[{rules,vR,state,vC},
(*			vVoidStates=(VoidStates/.List[UnionRules[optsRV,opts]])/.Options[TimeSeries]; *)
			(* Standard Method *)
			vR = v/.Thread[states->0];
			state = 0;
			vC=Table[
				If[vR[[nn]]>0,state=vR[[nn]],state]
			,{nn,1,Length[v]}];
			TimeSeries[Select[vC,#!=0&]]
		];


(* PCCA+ *)

TransitionMatrix/:PCCASoftClusterMatrix[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],level_] := 
		Module[{eqDist,size,cM,eVecs,distTable,maxDist,cornerList,eVecsNorm,list2,ll},
			size=Length[m];
			eVecs=Transpose[Eigenvectors[N[m],level+50]];
			distTable=Flatten[Table[{qq,pp,Norm[eVecs[[pp]]-eVecs[[qq]]]},{pp,2,size},{qq,1,pp-1}],1];
			maxDist=Max[distTable[[All,3]]];
			cornerList=Select[distTable,#[[3]]==maxDist&][[1,{1,2}]];
			For[ll=3,ll<=level,ll++,
				distTable=Table[{pp,Apply[Plus,Map[Norm[eVecs[[pp]]-eVecs[[#]]]&,cornerList]]},{pp,Complement[Range[size],cornerList]}];
				maxDist=Max[distTable[[All,2]]];
				AppendTo[cornerList,Select[distTable,#[[2]]==maxDist&][[1,1]]];
			];
			Print[cornerList];
			eVecsNorm=Map[#-eVecs[[cornerList[[1]]]]&,eVecs][[All,Range[2,level]]];
   		 list2=eVecsNorm.Inverse[Drop[eVecsNorm[[cornerList]],1]];
			ClusterMatrix[Chop[Join[Transpose[list2],{Map[1-Apply[Plus,#]&,list2]}]]]
		];
TransitionMatrix/:PCCACrispClusterMatrix[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],level_] := 
		Module[{eqDist,size,cM,eVecs,distTable,maxDist,cornerList,eVecsNorm,list2,ll},
			size=Length[m];
			eVecs=Transpose[Eigenvectors[N[m],level]];
			distTable=Flatten[Table[{qq,pp,Norm[eVecs[[pp]]-eVecs[[qq]]]},{pp,2,size},{qq,1,pp-1}],1];
			maxDist=Max[distTable[[All,3]]];
			cornerList=Select[distTable,#[[3]]==maxDist&][[1,{1,2}]];
			For[ll=3,ll<=level,ll++,
				distTable=Table[{pp,Apply[Plus,Map[Norm[eVecs[[pp]]-eVecs[[#]]]&,cornerList]]},{pp,1,size}];
				maxDist=Max[distTable[[All,2]]];
				AppendTo[cornerList,Select[distTable,#[[2]]==maxDist&][[1,1]]];
			];
			eVecsNorm=Map[#-eVecs[[cornerList[[1]]]]&,eVecs][[All,Range[2,level]]];
			list2=eVecsNorm.Inverse[Drop[eVecsNorm[[cornerList]],1]];
			Crispify[ClusterMatrix[Chop[Join[Transpose[list2],{Map[1-Apply[Plus,#]&,list2]}]]]]
		];


(* TransitionPathTheory *)

TransitionMatrix/:ReactiveFluxMatrix[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],stateA_,stateB_] := 
		Module[{eqDist,size,cM,tm},
			size=Length[m];
			tm=(m-IdentityMatrix[size]);
			cM=Committor[TransitionMatrix[m,opts],{stateA,stateB}][[2]];
			eqDist=Equilibrium[TransitionMatrix[m,opts]];
			Table[cM[[jj]]*(1-cM[[ii]])*eqDist[[ii]]*tm[[ii,jj]],{ii,1,size},{jj,1,size}]
		];
TransitionMatrix/:NetFluxMatrix[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],stateA_,stateB_] := 
		Module[{eqDist,size,cM,tm},
			size=Length[m];
			tm=(m-IdentityMatrix[size]);
			cM=Committor[TransitionMatrix[m,opts],{stateA,stateB}][[2]];
			eqDist=Equilibrium[TransitionMatrix[m,opts]];
			Map[Max[0,#]&,(#-Transpose[#])&[Table[cM[[jj]]*(1-cM[[ii]])*eqDist[[ii]]*tm[[ii,jj]],{ii,1,size},{jj,1,size}]],{2}]
		];
TransitionMatrix/:ReactivityVector[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],stateA_,stateB_] := 
		Module[{eqDist,size,cM},
			size=Length[m];
			cM=Committor[TransitionMatrix[m,opts],{stateA,stateB}][[2]];
			eqDist=Equilibrium[TransitionMatrix[m,opts]];
			Table[cM[[ii]]*(1-cM[[ii]])*eqDist[[ii]],{ii,1,size}]
		];
TransitionMatrix/:StateExitFlowRate[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],stateA_,stateB_] := 
		Module[{eqDist,size,cM},
			size=Length[m];
			cM=Committor[TransitionMatrix[m,opts],{stateA,stateB}][[2]];
			eqDist=Equilibrium[TransitionMatrix[m,opts]];
			Apply[Plus,Table[If[jj==stateB,0,1]*cM[[jj]]*(1-cM[[stateA]])*eqDist[[stateA]]*m[[stateA,jj]],{jj,1,size}]]
		];
TransitionMatrix/:StateExitFlowRateSensitivityMatrix[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],stateA_,stateB_] := 
		Module[{committor, equilibrium, committorSensitivity, equilibriumSensitivity, eqInverse, size, dummySum1, dummySum2, matrix1, matrix2, matrix3},
			size=Length[m];
			committor=Committor[TransitionMatrix[m,opts],{stateA,stateB}][[2]];
			equilibrium=Equilibrium[TransitionMatrix[m,opts]];
			committorSensitivity=CommittorSensitivityMatrix[TransitionMatrix[m],{stateA,stateB}];
			eqInverse=PseudoInverse[Join[Transpose[m]-IdentityMatrix[size],{Table[1,{size}]}]];
			equilibriumSensitivity=Table[equilibrium[[kk]]*(eqInverse[[stateA,kk]]-eqInverse[[stateA,ll]]),{kk,1,size},{ll,1,size}];
			dummySum1=Apply[Plus,committor*m[[stateA]]]-committor[[stateB]]*m[[stateA,stateB]];
			dummySum2=Table[Apply[Plus,committorSensitivity[[nn]]*m[[stateA]]]-committorSensitivity[[nn,stateB]]*m[[stateA,stateB]],{nn,1,size}];
			matrix1=SparseArray[{},{size,size}];
			matrix1[[stateA]]=committor-committor[[stateA]];
			matrix1[[stateA,stateB]]=-committor[[stateA]];
			matrix2 = Table[dummySum2[[kk]] * (committor[[kk]] - committor),{kk,1,size}];
			matrix3 = equilibrium[[stateA]]*(matrix1 + matrix2);
			equilibriumSensitivity * dummySum1 + Normal[matrix3]
		];


(* Plot functions *)

PositionListPlot[positionList_,valueList_,opts:OptionsPattern[DensityPlot]] := 
		Module[{pointList},
			pointList=If[NumberQ[positionList],
				Table[{Mod[ii,positionList]+1,Floor[ii/positionList]+1},{ii,0,positionList^2-1}],
				positionList
			];
			ListPlot[Map[pointList[[#]]&,valueList],opts]
		]
PositionDensityPlot[positionList_,valueList_,opts:OptionsPattern[PositionDensityPlot]] := 
		Module[{arrayToPlot,verticesLabel,vertexList,epilog,optslist},
			arrayToPlot=If[NumberQ[positionList],
				Partition[valueList,positionList],
				Transpose[Join[Transpose[positionList],{valueList}]]
			];
			vertexList={};
			If[(Vertices /. {opts} /. Options[PositionDensityPlot])===True,
				verticesLabel=VertexLabel/. {opts} /. {VertexLabel->False};
				vertexList = If[(verticesLabel) === False,
	     	           Map[Point,positionList],
						If[!VectorQ[verticesLabel],
							MapIndexed[Text[#2[[1]],#1]&,positionList],
		                    MapIndexed[Text[verticesLabel[[#2[[1]]]],#1]&,positionList]
						]
		            ];
			];
			epilog=Join[{Epilog}/.{opts}/.{Epilog->{}},{VertexStyle}/.{opts}/.{VertexStyle->{}},vertexList];
			optslist = Select[List[opts],!(SameQ[#[[1]],VertexStyle] || SameQ[#[[1]],VertexLabel] || SameQ[#[[1]],Vertices] || SameQ[#[[1]],Epilog])&];
			ListDensityPlot[arrayToPlot,Epilog->epilog,optslist]
		]
PositionStreamDensityPlot[positionList_,valueList_,opts:OptionsPattern[Join[{InterpolationOrder->3, FieldStyle->{Gradient}},Options[StreamDensityPlot]]]] := 
		Module[{arrayToPlot,interpolation,func,x,y,nOrder},
			arrayToPlot=If[NumberQ[positionList],
				Transpose[Join[Transpose[Table[#&[{Mod[ii,positionList]+1,Floor[ii/positionList]+1}],{ii,0,positionList^2-1}]],{valueList}]],
				Transpose[Join[Transpose[positionList],{valueList}]]
			];
			nOrder = InterpolationOrder/.Flatten[Join[List[opts], {InterpolationOrder->3}]];
			interpolation = Interpolation[arrayToPlot, InterpolationOrder->nOrder];
			func=OptionValue[FieldStyle]/.{Gradient->{D[interpolation[x,y],{x}],D[interpolation[x,y],{y}]},Rotation->{D[interpolation[x,y],{y}],-D[interpolation[x,y],{x}]}};
			StreamDensityPlot[Evaluate@func,{x,interpolation[[1,1,1]],interpolation[[1,1,2]]},{y,interpolation[[1,2,1]],interpolation[[1,2,2]]},Evaluate[RemoveRules[InterpolationOrder,RemoveRules[FieldStyle,opts]]],BaseStyle->{FontFamily->"GillSans",FontSize->12}]
		]
PositionStreamPlot[positionList_,valueList_,opts:OptionsPattern[Join[{InterpolationOrder->3, FieldStyle->{Gradient}},Options[StreamPlot]]]] := 
		Module[{arrayToPlot,interpolation,func,x,y,nOrder},
			arrayToPlot=If[NumberQ[positionList],
				Transpose[Join[Transpose[Table[#&[{Mod[ii,positionList]+1,Floor[ii/positionList]+1}],{ii,0,positionList^2-1}]],{valueList}]],
				Transpose[Join[Transpose[positionList],{valueList}]]
			];
			nOrder = InterpolationOrder/.Flatten[Join[List[opts], {InterpolationOrder->3}]];
			interpolation = Interpolation[arrayToPlot, InterpolationOrder->nOrder];
			func=OptionValue[FieldStyle]/.{Gradient->{D[interpolation[x,y],{x}],D[interpolation[x,y],{y}]},Rotation->{D[interpolation[x,y],{y}],-D[interpolation[x,y],{x}]}};
			StreamPlot[Evaluate@func,{x,interpolation[[1,1,1]],interpolation[[1,1,2]]},{y,interpolation[[1,2,1]],interpolation[[1,2,2]]},Evaluate[RemoveRules[InterpolationOrder,RemoveRules[FieldStyle,opts]]],BaseStyle->{FontFamily->"GillSans",FontSize->12}]
		]
PositionStatePlot[pStateList_,positionList_,pOptions:OptionsPattern[PositionStatePlot]] := 
		Module[{pointList, pFill={}, pGrid={}, pText={}, pOutline={}, pColoring={}},
			pointList=If[NumberQ[positionList],
				Table[{Mod[ii,positionList]+1,Floor[ii/positionList]+1},{ii,0,positionList^2-1}],
				positionList
			];
			If[OptionValue[PositionStatePlot,optsKM,Outline],
				{vorvert,vorval}=VoronoiDiagram[data2D];
				pOutline=DiagramPlot[positionList];
			];
			If[Grid/. replaceRules,
				pGrid={RGBColor[0,0,0],Thickness[0.012`],GridStyle/. replaceRules,(Line[{#1[[2]]+zustandsabstand {0,0},#1[[2]]+zustandsabstand {1,0},#1[[2]]+zustandsabstand {1,1},#1[[2]]+zustandsabstand {0,1},#1[[2]]+zustandsabstand {0,0}}]&)/@({#1,GetPositionFromState[#1]}&)/@pStateList}
			];
			If[Text/. replaceRules,
				pText={RGBColor[1,1,1],Style[Function[(Text[ToString[#1[[1]]],#1[[2]]+{0.5`,0.5`} zustandsabstand])]/@({#1,GetPositionFromState[#1]}&)/@pStateList,BaseStyle/. replaceRules]};
			];
			If[Fill/. replaceRules,
				pFill={RGBColor[0.88`,0.47`,0],FillStyle/. replaceRules,({Rectangle[#1[[2]]+zustandsabstand {0,0},#1[[2]]+zustandsabstand {1,1}]}&)/@({#1,GetPositionFromState[#1]}&)/@pStateList};
			];
			If[Coloring/. replaceRules,
				pColoring={({(ColoringFunction/.replaceRules)[(ColoringData/.replaceRules)[[#1[[1]]]]],Rectangle[#1[[2]]+zustandsabstand {0,0},#1[[2]]+zustandsabstand {1,1}]}&)/@({#1,GetPositionFromState[#1]}&)/@pStateList};
			];
			ret=Graphics[Flatten[{pFill,pColoring,pText,pGrid,pOutline}]];
			ret
		];


(* I/O functions *)

ReadSparseMatrix[fileName_,opts:OptionsPattern[Import]]:=Module[
		{readImport},
		readImport=Import[fileName,"Table",opts];
		SparseArray[Map[1+#[[{1,2}]]->#[[3]]&,readImport]]
]
WriteSparseMatrix[fileName_,m_?MatrixQ,optsE:OptionsPattern[Export]]:=Module[
		{arrayToExport},
		arrayToExport=Map[{#[[1,1]]-1,#[[1,2]]-1,N[#[[2]]]}&,Drop[ArrayRules[SparseArray[m]],-1]];
		Export[fileName,arrayToExport,"Table",optsE]
]
WriteSparseMatrix[fileName_,TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],optsE:OptionsPattern[Export]]:=Module[
		{arrayToExport},
		arrayToExport=Map[{#[[1,1]]-1,#[[1,2]]-1,N[#[[2]]]}&,Drop[ArrayRules[SparseArray[m]],-1]];
		Export[fileName,arrayToExport,"Table",optsE]
]
WriteSparseMatrix[fileName_,RateMatrix[m_?MatrixQ,opts:OptionsPattern[RateMatrix]],optsE:OptionsPattern[Export]]:=Module[
		{arrayToExport},
		arrayToExport=Map[{#[[1,1]]-1,#[[1,2]]-1,N[#[[2]]]}&,Drop[ArrayRules[SparseArray[m]],-1]];
		Export[fileName,arrayToExport,"Table",optsE]
]
WriteSparseMatrix[fileName_,CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]],optsE:OptionsPattern[Export]]:=Module[
		{arrayToExport},
		arrayToExport=Map[{#[[1,1]]-1,#[[1,2]]-1,N[#[[2]]]}&,Drop[ArrayRules[SparseArray[m]],-1]];
		Export[fileName,arrayToExport,"Table",optsE]
]


End[ ];


EndPackage[ ];
